基于循环相关和LPSO算法的自适应MCKD方法的滚动轴承早期故障特征提取...
第36卷第22期
振动与冲击
JOURNAL OF VIBRATION AND SHOCK Vol.36 No.22 2017
基于循环相关和LPSO算法的自适应MCKD方法的
滚动轴承早期故障特征提取
陈昆弘，刘小峰
(重庆大学机械传动国家重点实验室，重庆400044)
主商要！针对强噪声情况下滚动轴承早期故障信号特征难以提取的问题，提出了 MCKD与对称差分能量算子解调 的特征提取方法。MCKD算法进行滤波时，滤波器长度L和故障周期T对滤波效果的影响至关重要，因此提出基于循环
相关和LPSO算法结合的自适应的MCKD算法，自动搜寻MCKD算法所需最优参数；原信号经滤波后，故障特征被明显突
出，为了剔除剩余噪声，对滤波后信号进一步做对称差分能量算子解调，剔除剩余噪声同时获得解调谱，进而提取滚动轴
承的早期故障。实验分析验证了该方法的有效性。
关键词：循环相关；局部粒子优化；最大相关峭度解卷积；对称能量算子解调；早期故障；特征提取
中图分类号：TH133.3; TH17 文献标志码：A DOI ： 10. 13465/j. cnki. jvs. 2017. 22. 013
Incipient fault diagnosis of rolling element bearing based onadaptive
maximum correlated kurtosis deconvolution
CHENKunhong，LIU Xiaofeng
(The State Key Laboratory of Mechanical Transmissions，Chongqing University，Chongqing 400044，China) Abstract: Aiming at the problem that the feature of incif)ient fault of rolling element bearing is dfficult to be extracted under strong noise background，a fault diagnosis method based on MCKD and sym operator demodulation was p roposed. The filter size L and the period of interesting signal T play an important role in MCKD filtering. An adaptive MCKD filter based on cyclic correlation a
nd LPSO was proposed，which coul best parameters automatically. The feature after filtering was outstanding，but some residual noise was still remained. To reduce the residual n oise and get the demodulation spectrum，symmetrical differencing energy operator demodulation was applied after filtering. The result of incipient f ault diagnosis of rolling element bearing shows that Key words：cyclic correlation ；local particle swarm optimization ( LPSO ) ; maximum correlated kurtosis deconvolution ( MCKD) ；symmetrical differencing energy operator demodulation; incipient fault; feature extraction
滚动轴承出现早期局部故障时，故障产生的冲击成分一般十分微弱，往往会被强大的背景噪声淹没而无法提取。针对该问题，不少学者进行了深入研究。莫代一等[1]提出一种基于双重j因子的故障诊断方法，可将轴承早期故障信号分解成高共振和低共振两个分量，对低共振分量进行分析即可提取出轴承故障特征信息，但该方法的参数过多，如果任一参数设定不合理都会对最终结果产生不利影响；曾庆虎等[2]利用 小波相关滤波的降噪特性，将相关滤波降噪与包络谱
基金项目：国家自然科学基金项目（51475052)；中央高校基本科研业务 费（106112016CDJZR115502)
收稿日期：2016 -05 -24修改稿收到日期：2016 -09 -06
第一作者陈昆弘男，硕士,1991年生
通信作者刘小峰女，博士，教授，博士生导师,1980年生相结合，提出了基于小波相关滤波包络分析的轴承早期微弱故障特征提取，诊断效果明显，然而该方法中小 波基函数的选择缺乏自适应性，选取不同的小波基函数时，诊断结果差异较大；王宏超等[3]将最小熵卷积(Minimum Entropy Deconvolution，MED)用于早期轴承信号的降噪，取得了较好的效果。McDonald等[4]对 MED进行了改进，提出了最大相关峭度解卷积（Maxi­
mum Correlated Kurtosis Deconvolution，MCKD)方法，以相关峭度为评价指标，充分考虑了信号所含冲击成分的周期特性，通过迭代过程实现解卷积运算，突出信号 中被强烈噪声所掩盖的连续脉冲，并用于齿轮故障诊断，取得了良好的效果。
钟先友等[5<；研究发现，采用MCKD对滚动轴承早期振动信号进行降噪，也取得了良好的结果，但有限 冲击响应滤波器长度和信号周期这两个参数决定了滤
第22期陈昆弘等：基于循环相关和LPSO算法的自适应MCKD方法的滚动轴承早期故障特征提取81
波效果，实际过程中，往往采用试验的方法确定，浪费
时间，故需要采用一种自适应的方法，使得M C K D方法
取得较优效果。
本文提出将L P S O(Local Particle Swarm Optimiza­
tion)算法和循环相关结合，对 M C K D 的影响参数进行
寻优，实现自动寻最优参数的自适应M C K D算法，并
结合对称能量算子解调[7]，得到滤波后信号解调谱，分
析故障特征。实例分析表明，该方法可有效提取滚动
轴承的早期故障特征。
1 MCKD原理
为了提取周期故障信号，MCKD通过选定一个有限
冲击滤波器/使周期已知信号滤波后的相关峭度最大，
从而使得信号中的冲击成分最突出。相关峭度的定义为
# (n y〇-n T)2
DG m(,) 2max^-------
f(#K)2
0=1
/ 2 L f，/2,…，/!]T(1 )
式中：K〇为周期信号；，为信号的周期；/为滤波器向
量；！为滤波器长度； >为位移数，增加位移数能增加
算法的周期序列脉冲数。>一般取值为1 ~ 7，当>的
取值大于7时，会使计算超出浮点数范围，因此本文选
取的 >值为7。
为了选取一个最优的滤波器/，使得CG>(,)最大，
根据函数求极值，即令
R f
CG M(,) 2 0, 6= 1，2,…，！
算得的结果及其矩阵表示为
f= 2"(1〇I，) _1 #J n r«n
2 I I$I I™=〇
其中，
I
C1_r%_r
0 ._,._,
0 0 _,
⑵
⑶-0 0 0 0 x N_L_,,
,= [0,，2,，…，n,]，
_1$$$
71-3,(71，71-,，…，^-肌,）
-1/222X
y2-n,(y2，h-,，…，h-n,) %n = .
-K N-n,(K N，K N_,，*，K N_n,)-
K1，K1_,，…，K1_n,
八K，J l-T，…，J l-m T
-7N，y N-T，…，y N-mT-
在M C K D算法中，滤波器长度L和故障信号周期 ，两个因素决定着滤波器滤波效果。本文提出了改进
的局部 PSO(Particle Swarm Optimization)算法和循环相 关结合的方法，自动获得L和，两个参数的最优解，使 得M C K D能够对周期未知信号进行滤波，扩大了M C K D滤波算法的适用范围。
2改进的LPSO算法
P S O算法是一种粒子优化算法，具有较好的全局寻优能力，假设在D维空间中有 >个粒子，即其组 成的粒子空间"=(1,I,…，1>)，每个粒子是一个 D维向量，其第（个粒子的向量为" = (.，.，…，.l)，代表第（个粒子在D维搜索空间中的位置，同时，第（个粒子的速度为M = (81，8，…，)，其个体历 史最优解尽=(:1，：2，…，：…)，整个空间的全局最优 解为％= ()1，）2，…，）…)，其粒子随着迭代的进行，速 度和位置变化更新公式为
8+1=〇X8 5 cl X rand X(N - .9)5
c2 X rand X (G- x6d)(4)
〜=、d+'d(?)式中:〇为惯性因子；（=1，2,…，>;d= 1，2,…，L; 6为当前的迭代次数;P，P为粒子的学习因子。r n d是 介于[0，1]的随机数。
以上就是标准P S O算法。然而，标准P S O算法存 在易陷入局部最优的缺点。故本文是用改进后的局部P S O算法（L P S O)，改进了速度更新公式，将粒子的全局历史最优解）替换为粒子的邻域内粒子
的最优解，使得粒子速度更新不再依赖全局最优，不易陷 入局部最优解。改进后的速度更新公式为
8+12〇X8 + P X rand X (: - .)+
p X r n d X(:net - .d)(6)式中，：nex为第（个粒子邻域内粒子的最优解。
从参考文献[8 ]发现，为使粒子算法中的速度变 化取得更好的收敛效果，从而得到更精确的结果，惯性 因子O按照凹函数变化时优于线性变化，线性变化优于常量，因此本文中O采用凹函数随迭代次数而变化，公式如下
J X(-^n x)X
(2 X—a_)+ 〇m (7)
n
式中:〇m ax为〇的最大值；〇m i为〇的最小值；6为粒子 此次迭代次数；m a g e n为P S O算法的总迭代次数。
在标准粒子算法中，学习因子P'P—般是恒定 的常值，但是依据实验表明，随着迭代而变化的学习因 子可以获得更好的收敛效果，即在粒子进化初期，使粒
82振动与冲击2017年第36卷
子能仔细的在自己的邻域内搜索，防止粒子较快地收敛于局部最优解；而在进化后期，应使粒子更快、更 准确地收敛于全局最优解[>]。即在初期使P的值较 大，而在后期使P的值随着迭代次数增大。故改进的 学习因子公式为
P= 4 - I(8)
c2 = 4 - cl(9)使用粒子进行寻优解时，需要选取一个适应度函数来衡量寻优效果。每当粒子更新自己的位置时，根据适应度函数，计算一次适应度的值，再比较选出最 优解。依据钟先友等的研究，结合实际中的实验表明，信噪比在衡量M C K D滤波效果方面，比各种熵值衡量效果较优，故本文选取信噪比作为粒子的适应度函数，用于评判粒子在M C K D滤波效果方面的优劣。即 对M C K D滤波后信号取包络，计算功率谱，再计算信噪 比
.=10x l〇g( S(’0))
、$(/o W
(10)
R(/0) = |1(/0)|2(11)式中:R(/o)为滤波后信号的功率谱；$(/〇)为原信号的 功率谱；/〇为滤波后信号在故障频率附近的平均值。信噪比越大，表明滤波效果越好。
在本文处理信号时，因为未提前预知信号故障频率/〇，故需要对信号进行预处理。本文提出，通过循环 相关算法，对信号频率进行预筛选，得到可能的信号频 率7 [/1/2/3…]，再通过L P S O算法得到精确值。
3循环相关
M C K D滤波器中，需要确定[!，,]两个变量，但是 在测试中发现局部P S O算法在处理多变量多峰函数时，即含噪声信号时，依然得不到足够精确的解。故本 文在局部P S O算法前，加入循环相关算法，在信号进入 局部P S O算法处理前，利用循环相关，对信号的周期, 做一轮筛选，使得局部P S O算法对变量，的搜索范围 极大的减小，提高整个算法的精确程度。循环相关分析中，所提取的信号周期定义为循环周期，，频率称之 为循环频率％全部循环频率的集合定义为循环频率域[10]。
定义信号的自相关函数为
Q(/&) =F〇[./-f).(/+ f)] (12)
Q(t，T)=(13)
式中，傅里叶系数C Q被定义为循环频率％处的循环自相关函数，表示为
C Q.(&) ,*Q(/&)e<2!%d/(14)
循环自相关函数能够将调幅信号的调制频率解调 到循环频率低频带上，实现与载波频率的分离。根据 循环相关结果图和对应切片图，对信号周期进行初步筛选，得到信号频率/=[/1/2/3…]，进而得到对应的周 期,= [,,,…]，在程序中，所需，的计算方法为,=77
此时，再初始化变量！，将得到的多个周期序列，与变量！组成一个二维变量[!，,，然后由L P S O自动 搜寻最优解。
4对称差分能量算子解调
能量算子解调相对于F F T(Fast Fourier Transform)解调而言，对噪声不敏感，能较好的突出信号特征部分，但是基本的能量算子解调算法在端点处以及突变点处仍比较敏感，为了尽量减少这些因素对信号解调的影响，故本文选取改进后的对称差分能量算子解调对滤波后信号进行处理。
离散信号的瞬时幅值一频率表示为
.(0)= a(0)cos(!(0))=
<(0)cos('( 0)+ () (15)对.0)进行能量算子计算得
#%(0)=.2(0)- .(0-1).(0 + 1) (16)又由于#(<(0) )+0，！+0，式（16)简化为
#(.(0))+a2(0)sin2('(0)) (17)对称差分能量算子解调对能量算子中的差分序列 进行了改进。首先，对.0)的差分序列定义为
k(0)=.051) -.0-1)(18)
K(0)的差分序列定义为
〜,、.0+2)-2.(0) + .(0-2)/ir*、
K(0) = -------------4-------------(19)
得到改进的算子
#1(.( 0) )=
#[.(0 + 1)] +2#[.(0)] +#[.(0-1)](2〇)
式中为共轭运算；&为时间延迟变量。
当./为二阶循环平稳信号时，Q(/&)随时间 变化呈现出周期特征，即Q(/，&) =Q(/ + ,，&)，,=1/a表示循环周期。因此，二阶循环平稳信号的自 相关函数Q(/&)可以用傅里叶级数展开为
使用传递函数求解，得到能量算子的幅、频估计值为
I a(0)I+
________2#[.0) ]_________(21)
槡#[.(0 + 1)- .0-1)]
第22期
陈昆弘等：基于循环相关和LPSO 算法的自适应MCKD 方法的滚动轴承早期故障特征提取
83
2
4 6 8
10
频率/k l 〇2/Hz
图2
复合信号的H H ert 包络解调谱
Fig. 2 The Hilbert envelop spectrum of composite signal
对模拟使用基本L P S O 算法和对称能量算子解调 对模拟信号进行处理，设置算法中惯性因子'=〇. 7,
£；1=1.49，£；2=1.49，变量为[!，,，当迭代至第18次 时，得到最优解[4，,]
=
[131,500]，再将最优解代
人M C K D 算法，得到滤波后信号，再对滤波后信号用对 称差分能量算子解调，处理结果如图3所示。由图3 可知，虽然可以辨认故障频率，但杂频太多，效果并不 明显。
图1
复合信号的时域波形及频谱图
Fig. 1 The time domain and frequency domain of composite signal
°0 123456789 10
频率/k l 〇2/Hz
图3信号进行LPSO 和对称能量算子解调结果
Fig. 3 The result of LPSO and demodulation
利用循环相关对模拟信号进行筛选，得到循环相 关结果及&=1时切片图，如图4所示。由于滚动轴承
故障频率在低频段，故只需截取0 ~ 500 H z 频带进行筛 选，由图明显可知，可选100 3
及其倍频程为所需频
率，根据公式,=///…得到筛选后多个周期，序列，=
[
200，100，67，50 ]，初始化变量！，组合成所需寻优组合 [
!，,
。
将变量[!，，和模拟信号输入改进L P S O 算法
中，自动对[!，，进行筛选，在第6次迭代时，即当！= 343，, = 200时，得到最优解，将最优解[!0，,
]
=
[
343，200 ]输入M C K D 方法进行滤波，得到滤波后信号
如图5所示。再将滤波后信号用对称差分能量算子解 调，得到解调谱，如图6所示。
'(0) +
5
故障分析
在实际分析中，选择合适[!，，参数，突出被噪声
淹没的特征成分，是利用M C K D 进行滚动轴承早期故 障特征分析的关键。利用循环相关对，参数进行第一 轮筛选后，再利用改进L P S O 对[!，,]参数进行第二轮 筛选，既避免了人为设定的主观性，又提高了对[!，, 两参数选择的准确性。在信号经M C K D 滤波后，信噪 比会显著提高，但仍存在部分噪声。利用对称差分能 量算子解调对噪声的不敏感特性，可有效地剔除噪声 并得到解调谱。具体故障分析步骤如下：
步骤1利用循环相关筛选的多个周期合成的周
期数组,
步骤2
设定L P S O 算法的粒子种大小 > 和 迭代次数$，设置惯性因子的最大值〇m a 及最小值 〇m i ，初始化！，并与由步骤1得到的，组成二维变量数 组[!，,]。
步骤3
将[!，,]和信号输入改进L P S O 中，利用
改进L P S O 对M C K D 的两个影响参数进行寻优，得到最 优解[!0，,0]。
步骤4
根据[!0，,]，输人M C K D 滤波器对信号
进行滤波，得到滤波后信号。
步骤5
对滤波后信号根据对称差分能量算子解
调， 得到信号的解调谱。
步骤=根据解调谱分析信号，判断信号故障
类型。
在本文中，默认粒子种大小> =30，迭代次数$
=
20，°max  =0 -9,O m i n  = 0- 4。
6
仿真信号分析
滚动轴承的故障大部分都是周期性冲击信号，故
仿真信号采用周期冲击信号，其模拟表达式为
.(/) = 5 x  e  4(/ 4,) x  sin [2!,(/ - A , ] X
滚动轴承的特点
H (/ - k T )
+ n (t ) k  = 1，2,3，*
(23)
式中:/为信号的固有频率，设为3 k H z ;，为故障冲击 周期，设为1/100；阻尼系数c =0. 1，"（/-/^)为单位 阶跃函数；n (〇为噪声。本次模拟中，采样频率为
20 k H z 。噪声为高斯白噪声，信噪比为-5，其时域波 形和频谱如图1所示。
对信号进行Hlbert 包络解调得到波形，如图2所 示。由图2可知，故障频率可以识别，但是有太多的杂 频，不易快速辨认。
3-0X F W S
A /v o l x F l }
®
84
振动与冲击
2017年第36卷
对信号进行Hilbert 包络解调，得到波形图如图8
所示。由图可知，Hilbert 包络解调已经不能准确识别 信号成分。
循环频率/k 1 〇2/H z  循环频率/"x  1 02/H z (K 循环相关结果
（b)
& = 1时切片图
图4
循环相关处理结果
Fig. 4 The result of cyclic correlation
由图6可知，信号的故障频率为100 3及其倍频 程，如图中/及其倍频程标示，与预设值相符，且其他频 带对故障的识别无影响。表明利用局部PS 0算法和循 环均值相结合的自适应M C K D 滤波器能对模拟信号有 较好的滤波效果。
I
实验验证
本文采用西储大学公开的轴承故障数据进行实验
验证。西储大学在实验中采用一台2马力的电机进行 驱动，加速度传感器采集信号，信号通过16通道的
D A T 记录器记录。本文使用采样率12 k H z 条件下的
107 m a t 文件中的驱动端数据进行验证，并在实验前， 加上-1d B 的高斯白噪声。其中主要参数如下：故障 深度0. 178 m m 英寸，载荷为1.491 4 k W ，实际转速 1 772 r/min ，内圈故障。取2 048个点进行分析，其时 域波形及频谱如图7所示。
馨
时间?Xl 〇-2/s 频率户取/他
(a )时域波形图 （b )频谱图
图7信号时域波形图及其频谱
Fig. 7 The time domain and frequency domain of signal
2
4
6
8 10
时间"1(T2/s 图5
滤波后信号
Fig. 5 The signal after
filtePing
2
4
6
8
10
频率/x i 〇2/Hz 图6复合信号对称差
分能量算子解调谱
Fig. 6 The symmetrical differencing energy operator demodulation spectrum of composite signal
对信号以（！，，为未知量，使用基本LPS 0进行筛 选，再对滤波后信号进行对称差分能量算子解调，其中 设置算法中惯性因子'=0.7,P  =1.49,P  =1.49，在 迭代至第13次时，得到最优解[!0，,]=[ 10，500 ]，将 最优解代人M C K D 方法，得到滤波后信号，再对滤波后 信号进行对称能量算子解调，解调结果如图9所示。 由图9可知，几乎无法准确识别故障特征。
Fig. 9 The result of LPSOand symmetrical differencing
energy operator demodulation
对信号做循环相关，并对得到频率通过切片进行 筛选，循环相关结果在& = 1时切片图，如图10所示。 因大部分的滚动轴承故障都集中在低频段，故只需截
取0 ~ 500 H z 范围选取频率，取幅值大于0 ~ 500 H z 幅 值均值点对应频率作为筛选频率，将筛选后的频率求 得对应，=[187，145，86，76，40，34，24 ]，初始化变量 !，组成二维变量矩阵[!，,，与原始信号一同导入改进
LPS 0算法中自动寻优，在迭代至18次时，得到最优解
, = 76, ！ =496，将最优解[!。，,] = [496,76]输入
M C K D 方法，得到滤波后信号后，如图11所示。再对
信号进行对称差分能量算子解调，得到解调谱如图12 所示。
由图12可知，故障信号为158. 2 H z 及其倍频程， 如图中/及其倍频程所标示，与理论计算结果趋于一
致，且其余信息对故障识别几乎没有影响，效果较优。 说明本文所述方法能较好地对早期滚动轴承故障信号 进行处理。
%
2
4    6 8 10
频率/x l 〇2/Hz
图8
Hilbert 包络处理结果
Fig. 8 The Hilbert envelop of signal
5
5 0
2 2 .—
- 1
b x F ^S 0
•0.0.0.0.0•0•
m s o
8 6 4 2
0-21
-4I
乙

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。