基于自适应EEMD及改进滑动峰态解调的滚动轴承故障提取方法

2019-07-22杨伟才廖英英刘永强

轴承 2019年1期

杨伟才，廖英英，刘永强

(石家庄铁道大学 a.机械工程学院；b.土木工程学院，石家庄 050043)

滚动轴承所处工作环境非常复杂，早期故障的冲击信号极其微弱，且容易湮没在强背景噪声中，故障特征很难被直接识别出来[1]。因此，众多学者对振动信号降噪方法和故障提取进行了大量的研究。

小波降噪的优点在于具有多分辨率性，但其消噪过程中需选择基函数和阈值，缺乏自适应性[2]。经验模态分解[3](EMD)适合处理非线性和非平稳信号且具有自适应性，但严重的端点效应和模态混叠现象影响分析结果的正确性和精确性。因此，利用高斯白噪声具有频率均匀分布的特征，通过在EMD的基础上加入高斯白噪声提出了集成经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)[4]，可以有效解决EMD的模态混叠问题。依据互相关系数和峭度值筛选IMF分量并进行信号重构对原信号具有一定的降噪作用[5-6]。

传统的包络解调很难提取出弱冲击成分。零均值、零时滞的四阶累积量即峰态对信号的冲击成分非常敏感[7]，基于峰态特性提出的滑动峰态算法(Sliding Kurtosis Demodulation，SKD)在机械故障的弱冲击特征提取得到了成功应用[8-9]。

综上分析，提出了一种基于自适应EEMD和改进滑动峰态解调的滚动轴承故障提取的方法。首先对故障信号进行EEMD处理得到IMF分量；然后依据相关峭度、信息熵和能量加权准则筛选IMF分量并进行重构，对重构信号进行谱峭度带通滤波后计算其滑动峰态时间序列；最后通过Fourier变换分析识别弱冲击故障信息。

1 理论基础

1.1 EEMD原理

EEMD方法实质上是对EMD算法的一种改进，在原信号中加入足够多组不同的白噪声后，仍然进行EMD处理，再根据白噪声均值为零的随机特性对分解结果进行平均处理，以消除白噪声影响。白噪声的加入可以为EMD处理提供一个相对一致的参照尺度分布，保证每个模态函数时域的连续性从而减小模态混叠现象。

设信号为x(t)，具体的分解步骤如下：

1)在x(t)中加入随机白噪声ni(t)，构成新信号xi(t)，即

xi=x(t)+ni(t)；i=1,2,…,M。

(1)

2)对待测信号xi(t)进行EMD处理，得到

(2)

式中：N为分解的IMF数量；ci,n(t)为分解得到的IMF分量；ri,n(t)为余量。

3)每次加入不同的白噪声序列，重复前2个步骤M次，对分解出的IMF分量求总体平均，消除白噪声造成的影响，最终得到EEMD处理的IMF分量cn(t)，即

(3)

1.2 IMF优选准则

1.2.1 相关峭度准则

由于边界效应、包络误差的存在，EEMD处理会产生与原始信号无关的虚假分量，仅仅用峭度准则或者互相关系数准则优选IMF分量很难得到满意的结果。相关峭度(Correlated Kurtosis, CK)兼备峭度和相关函数的特征，是反映振动信号中周期脉冲信号强度的参数[10]，其表达式为

(4)

式中：xn为信号；T为感兴趣脉冲信号的周期；N为采样长度;M为偏移周期个数。

一般情况下，感兴趣周期就是冲击信号周期，预先设置好冲击信号周期之后，根据相关峭度值就可以判断出冲击信号的强弱，其值越大，冲击性越强。

1.2.2 信息熵准则

信息熵是描述系统混乱程度的物理量，当轴承正常工作时，其振动信号杂乱无章，时频分布聚集性差，熵值越大；当发生故障时，振动信号表现为有规律的周期振动，时频分布聚集性强，熵值越小[11-12]。计算IMF分量的信息熵步骤如下：

1) 设ci(t) 为振动信号c(t)的第i个IMF 分量，寻找ci(t) 中的最大值和最小值，分别定义为cmax和cmin。

2) 设置区间[cmin，cmax]中的N个等分点Ai，[cmin，A1]，…，(Ai，Ai+1]，…， (AN-1，cmax]，将这N个区间作为特征量的N个离散值域B∈{B1，B2，…，BN}。当某个样本的属性值落在(Ai，Ai+1]区间时，就认为该样本在该属性上具有相应的离散属性值Bi。

3)IMF分量的时域点个数为n，c(t)落在第i个区间的样本点个数为mi，则其在第i个区间的概率Pi=mi/n。由此可得该IMF 的信息熵为

(5)

1.2.3 能量加权准则

当滚动轴承发生故障时，滚子在通过故障时会引发强烈的冲击，故障引起的冲击能量必然比正常运行时的能量大。EEMD可以将信号分解为不同的IMF分量，能量大的IMF分量，其故障程度越大。

振动信号的能量是指信号幅值绝对值的平方和，即

(6)

式中：En为振动信号能量值；xi为i时刻的幅值。

根据相关峭度、信息熵筛选出IMF分量，对其能量进行求和可得

(7)

式中：E为筛选出的IMF分量总能量；Ei表示第i个IMF的能量值；N为筛选的IMF阶数。

为了更好地表征各IMF 分量的影响，引入能量权系数λi表征不同模态分量对总能量的影响[13]，即

(8)

1.2.4 自适应优选准则

根据相关峭度、信息熵和能量加权各自的性质，将三者进行结合，能够选出最有效的IMF分量，不易丢失重要故障信息，还能够增强重构信号的冲击性，提高信噪比，从而更加容易提取故障特征频率。自适应优选准则如下：

1)计算各个IMF分量的相关峭度值，并且从大到小进行排序。

2)计算出信息熵最小时所对应的重构个数(重构个数i：各个IMF分量相关峭度从大到小前i阶，即当i=1，选取相关峭度最大的IMF分量进行重构，当i=2，选取相关峭度最大和次大值的IMF分量进行重构，以此类推)。

3)确定IMF分量后，计算各个IMF分量的能量加权系数，将其作为IMF分量的权系数。

1.3 滑动峰态算法

设平稳随机信号为x(n)，n=-2,-1,0,1,2,…,其k阶累积量为

ckx(τ1,τ2,τ3,…,τk-1)=

cum[x(n)x(n+τ1)…x(n+τk-1)]，

(9)

则零均值平稳随机过程x(n)的k阶累积量为

ckx(τ1,τ2,τ3,…,τk-1)=

mean[x(n)x(n+τ1)…x(n+τk-1)]-

mean[g(n)g(n+τ1)…g(n+τk-1)] ，

(10)

式中：cum[*]为累积量；mean[*]为均值；g(n)为与x(n)具有相同二阶统计量的Gauss随机过程。

根据(10)式进一步推导可得零均值、零时滞的4阶累积量为

C4x=mean[x4(n)]-3{mean[x2(n)]}2，

(11)

式中：C4x为峰态，可以描述随机序列概率分布曲线的尖峭或平稳程度，能有效反映信号偏离Gauss分布的程度。

基于上述分析，文献[8]提出一种时间序列的滑动峰态算法。对时间序列x(n)，n=-2,-1,0,1,2,…，通过求其滑动峰态序列CL(ti)，以突显信号的冲击特征。滑动峰态序列定义为

(12)

滑动峰态算法原理如图1所示。

图1 滑动峰态算法原理(L=3)

1.4 轴承故障诊断流程图

为解决强背景噪声下滚动轴承微弱故障特征提取问题，提出了一种基于自适应EEMD和改进滑动峰态解调的滚动轴承故障提取方法，其算法流程如图2所示，具体实施步骤如下：

图2 基于自适应EEMD和改进滑动峰态解调的故障诊断流程图

1)确定轴承测点，采集轴承振动信号。

2)采用EEMD对振动信号进行分解，综合运用相关峭度和信息熵准则选出最优IMF分量，进行信号重构以实现降噪；利用能量加权准则增强信号中的冲击成分。

3)利用快速谱峭图方法确定带通滤波器的最佳带宽和中心频率等参数，对信号进行带通滤波处理。

4)利用滑动峰态算法计算峰态时间序列，实现信号的包络处理。

5)对峰态时间序列进行Fourier变换，提取故障特征频率，进行故障诊断。

2 仿真信号分析

为验证新方法对轴承微弱故障诊断的有效性，构建由故障冲击信号、噪声信号和谐波信号3部分构成的仿真信号进行验证。单点损伤振动模型表达式为

(12)

式中：α为衰减率，α=800；A为冲击幅值，A=1 m·s-2；B为噪声幅值，B=2 m·s-2；C为谐波幅值，C=0.2 m·s-2；t为仿真时长，t=1 s；f1为冲击导致的共振频率，f1=3 kHz；z为随机数；f2为转轴转频，f2=35 Hz；fm为冲击信号的频率(即故障特征频率)，fm=120 Hz。

仿真信号的时域和频域波形如图3所示。由于增加噪声的缘故，导致冲击成分被噪声湮没，图3b中也找不出故障特征频率。将仿真信号进行EEMD处理得到13个IMF分量，各IMF分量的相关峭度如图4所示。重构个数与信息熵关系如图5所示。综合考虑图4和图5，选取IMF2，IMF3，IMF4，IMF5和IMF6并计算其能量权系数，结果见表1。利用能量权系数进行信号重构，结果如图6所示，与图3a对比可知信号的信噪比得到了增强。另外，由图7可知，滑动窗宽L=2时对应的信息熵值最小，峰态时间序列规律性最强；因此，选取L=2对重构信号进行滑动峰态解调处理，结果如图8所示。从图中可以清晰看出故障特征频率及其倍频，该仿真分析结果验证了改进方法的正确性和有效性。

图3 仿真信号的时域、频域图

图5 重构个数与信号熵关系

表1 IMF分量的能量权系数

图6 重构信号的时域图

图7 不同窗宽的信息熵值

图8 重构信号的滑动解调谱

3 试验信号验证

为进一步验证新方法的有效性，应用QPZZ-Ⅱ型旋转机械故障试验平台进行验证。试验轴承型号为6205，其主要参数见表2。采用激光加工技术在内圈沟道上加工一个边长0.2 mm，深0.1 mm的正方形凹槽模拟表面损伤类微弱故障。转速设定为882 r/min，采样频率为10 240 Hz，计算得内圈故障特征频率为79.6 Hz。

表2 6205型轴承主要参数

应用加速度传感器采集1 s的数据，故障信号的时域和频域波形如图9所示。将故障信号进行EEMD处理得到15个IMF分量，前5个IMF分量如图10所示。各IMF分量的相关峭度值如图11所示，重构个数与信息熵的关系如图12所示，综合考虑选取IMF2，IMF3，IMF4，IMF5和IMF11分量计算能量权系数，并将其作为重构IMF分量的权系数以增强信号的信噪比，结果见表3。

图9 试验信号的时频域波形

图10 试验信号经EEMD处理后的IMF分量

图11 相关峭度值

图12 重构个数与信息熵关系

表3 IMF分量的权系数

采用互相关系数-峭度方法与改进方法分别处理后，重构信号的时域波形如图13所示。从图中可以看出：改进方法较互相关系数-峭度方法得到的重构信号具有更加明显的冲击成分。

图13 本文方法和传统方法重构时域图

不同窗宽的信息熵值如图14所示，从图中可以看出，滑动窗宽L=5时，峰态时间序列规律性最强。因此，选取L=5对重构信号进行滑动峰态解调处理，结果如图15所示。从图中可以看出：对于未经降噪处理的信号，从其滑动解调谱中只能找到1倍故障特征频率，没有边频带和转频，受噪声的影响较大，整体诊断效果不佳；对于经过互相关系数-峭度处理的重构信号，其滑动解调谱中可找到1～4倍的故障特征频率，但转频只能找到1倍频，谱线不明显；经改进方法处理后，重构信号的滑动解调谱中存在明显的特征频率(78.69 Hz)及其倍频，这些特征频率周围还存在各自的边频带和转频及其倍频。特征频率与内圈故障频率79.6 Hz十分接近，且其频率特性符合内圈故障的频率分布，由此判断该轴承存在内圈故障，与实际情况相符。