利用经验模态分解和多尺度熵的滚动轴承故障诊断研究

秦喜文 ，高中华，董小刚，张邦成

QIN Xi-wen1～3 , GAO Zhong-hua1 , DONG Xiao-gang2, ZHANG Bang-cheng3

（1.长春工业大学 研究生院，长春 130012；2.长春工业大学 基础科学学院，长春 130012；3.长春工业大学 汽车工程研究院，长春 130012）

0 引言

在常见的机械设备中，滚动轴承是整个设备的重要部件，其在运作中的好坏直接影响到设备的生产效率，因此，准确地识别出运作过程中的状态至关重要。滚动轴承出现故障是因为设备内部各部件周而复始的运作外加载荷转速等其他原因，致使滚动轴承出现一些内外圈以及滚动体的故障，而且在故障出现的同时，伴随的振动信号也表现出剧烈的非线性和非平稳性的特征。传统对非平稳和非线性信号的处理方法，比如小波分析、分形维数、支持向量机等方法已大量的被用在滚动轴承故障的诊断领域，这些方法的渗透极大丰富了非平稳信号特征识别的研究体系[1～3]。

HHT是由N.E.Huang提出的一种更具有适应性的处理复杂非线性和非平稳信号的分析方法[4]，此方法通过经验模态分解(EMD)，使信号本身分解出一组互异的基底，且分解结果具有高度的自适应性。在分解过程中，传统的均值包络建立是采用三次样条插值伴随边界效应、过冲和欠冲等问题。S.R.Qin提出采用分段平均幂函数的方法求取信号的上下包络线抑制这种过冲现象[5]。E.Delechelle提出一种利用四阶抛物型偏微分方程解决均值包络问题的方法[6]。

Richman提出了改进的复杂尺度测量新方法即样本熵，样本熵具有得到稳定估计值所需的干扰能力强、抗噪声和参数取值范围内一致性好等显著特点[7]。由于样本熵的尺度单一导致只是衡量振动信号本身尺度下的复杂性。Costa在样本熵的研究基础上，提出了一种新的时间序列复杂度的衡量方法—多尺度熵[8～10]，这一次研究使熵的内容得到了全面的扩充。

上述两大类方法，国内研究人员都分别在故障诊断识别方面应用过，并取得了一定的成效，但对识别的效果及精确度仍有待改善。针对HHT，本文提出基于改进的经验模态分解，即利用小波函数作为插值函数求原始振动信号的包络，经过多次实验模拟，发现可以良好解决欠冲、过冲现象，基于这种研究下，再与多尺度熵结合起来应用到滚动轴承的振动信号当中，给出了一种全新的滚动轴承故障诊断识别方法。

1 基于小波函数的经验模态分解

EMD分解是HHT方法的核心，EMD方法将任意复合信号分解成有限固有模态函数(IMF)之和。EMD分解对一个x(t)信号进行分解采用步骤如下：

1）首先获得数据的最大值和最小值所有点，再利用小波函数作为基函数进行插值，并把最大值最小值同时连接，形成上下包络线，原始数据处于两线之间。

2）计算出上下包络均值m1，用原始数据减去均值：

判断h1是否符文IMF两个条件（信号零点的个数与信号极值点的个数相等或相差1和整个时间序列关于时间轴局部对称），如果符合，则得到第一个分量。

3）若不符合，把h1看做成原始数据，并把以上2步骤再次运行，得到新的均值m11，判断h11=h1-m1是否符合条件，不符合需再做k次处理，h1(k-1)-m1k=h1k，直至h1k符合为止，为C1=h1k，则C1是信号x(t)的第一个分量。

4）将C1从x(t)中分离出来，得到：

将r1看成原始数据，重复上述过程获得x(t)的第二个分量C2，依次做n次处理，即可获得n个分量，即：

5）最终得到一个单调的函数即循环结束。

式中，rn为残余项。

基于小波函数EMD的优点：

1）一层层将数据基于其本质特征筛分的过程，作用类似一组滤波器。

2）模态波形对称，时间尺度Ci从小到大依次分离。

3）有效的抑制边界效应、过冲和欠冲等问题。

2 多尺度熵

多尺度熵的计算是在样本熵的基础上，将原始数据粗粒化得到的在各个尺度上的样本熵值组成的一组数列，即时间序列在不同尺度下的样本熵。如果一个序列和另一个序列在同种尺度下，前者的熵值比后者高，这说明前者的时间序列的复杂性要高于后者，这就体现出了前者和后者两者之间具有强烈的差异性。多尺度熵的计算过程如下：

1）对给定原始信号数据X={x1,x2,x3,L,xN}，长度是N，即N为序列长度，给出嵌入维数m，并定义

2）定义X(i)与X(j)，它们之间距离D(i,j)为两个样本对应元素差值绝对值的最大值。

其中，i,j=1,2,3L,n-m,k=1,2,3Lm-1。

3）计算所有的X(i)与其他在子样本X(j)的距离，给出相似容限r，统计并计算出D(i,j)＜r的个数，标记成countD(i,j)最后算出平均值：

5）重构数据，令嵌入维数等于m+1，重复过程1-4步骤计算出Bm+1(r)。

6）单一尺度上时间序列的样本熵计算：

7）对于长度为N原始数据X根据尺度因子τ的取值不同把原始数据分割成长度为τ的数据单元，再将这N/τ组单元分别组内求均值，利用所求的均值点形成一个新的时间序列，这个过程也称将原始数据按τ的粗粒化过程。而单元组内求均值的相应点为：

8）对于τ的不同取值，也就形成了对原始X的多尺度化，对每个尺度化后新的时间序列依次求样本熵，即得到了多尺度熵。

多尺度熵的计算跟嵌入维数m、相似容限r以及尺度因子τ 都有关联。其中，m的取值一般为2，相似容限r=0.1～0.25*δ（其中δ为原始数据的标准差），对于尺度因子一般取值不超过20。

多尺度熵的优点：

1）特征明确，其熵值的大小直接反映出时间序列产生新模式概率的大小。

2）具有较强的抗干扰性。

3）对原始信号进行多尺度分析，使得分析更具完备性和系统性。

3 基于经验模态分解和多尺度熵的滚动轴承故障诊断

数据是从美国华盛顿凯斯西储大学实验室的滚动轴承数据中截断获取，被测轴承是SKF6205，轴承的损坏是由人为的加工制作完成，之后通过振动加速度传感器获得各种工作状态的振动数据。

本文数据是轴承转速在1797r/min的情况下采集，采样频率为12kHZ，主要分析滚动轴承的四种工作状态，分别是：滚动轴承正常、滚动轴承内圈故障、滚动轴承外圈故障、滚动体故障，每种状态下的数据共6个小样本，每个小样本数据长度为6000，且每个小样本采集用时30s。所有状态共计24个小样本，耗时12min。在进行多尺度熵计算过程中，m取值为2，相似容限r=0.15*δ，尺度因子τ=15。第一次进行样本采集得到的四种工作状态原始数据样本对比以及利用改进的EMD对四种工作状态分解得到的IMF分别如图1～图4所示。

图1 对正常状态进行改进的EMD得到的分解图

图2 对内圈故障状态进行改进的EMD得到的分解图

图3 对滚动体故障状态进行改进的EMD得到的分解图

从图1～图4中可以直观地判断出四种工作状态所描绘的振动信号有差异性，正常的工作状态与滚动体出现故障对比内圈出现故障和外圈出现故障的折射的波动性差异明显，但正常工作状态和滚动体出现故障仍然无法准确进行区别，当把振幅作为识别特征，内圈故障和外圈故障也出现了相互无法明确分离的情况。

将改进的经验模态分解应用到原始振动信号数据得到的IMF1和IMF2，利用多尺度熵进行计算，结果如图5和图6所示。

图4 对外圈故障状态进行改进的EMD得到的分解图

图5 IMF1在四种工作状态下的多尺度熵

图6 IMF2在四种工作状态下的多尺度熵

从图5中可以看出，第一，从四种状态的多尺度熵值聚合程度来看，滚动体故障和外圈故障相比原始数据多尺度熵值相对聚集。第二，图5中的四种工作状态的IMF1多尺度熵，外圈故障可以准确的与其他三种情况分离开来。其他工作状态仍然出现了熵值近似情况。出现这些情况，总的来说原始信号经过改进的EMD分解得到的IMF1虽然振动信号得到了平稳化，而且为原始信号的高频成份，但对于轴承振动而言，这并不能代表原始信号作为特征进行故障识别。而从图6看出，IMF2的多尺度熵值更加不理想，原始信号经过经验模态分解得到的各自IMF是从高频到低频，IMF1的频率高于IMF2，得到的IMF2是原始信号减去IMF1再进行改进的EMD分解而得，对原始信号而言，特征体现次于IMF1。

对得到的固有模态函数IMF1和IMF2取和进行多尺度熵计算，结果如图7所示。

图7 IMF1+IMF2在四种工作状态下的多尺度熵

分析图7可以明显发现四种状态基于多尺度熵值可以很好地分离开来，这说明了基于改进的EMD对信号进行分解后，得到固有模态函数IMF1与IMF2两者的组合具有代表原始数据特征的特性，也就意味着要从数据的整体性把握对振动信号的故障分析，加上后续利用多尺度熵的方法，对重构后的新数据进行计算后，得出结论，尤其是在尺度因子选为12即τ=12时，可以很好地进行滚动轴承振动信号的故障分离与识别。

基于上述研究，利用MATLAB的LIBSVM包对滚动轴承故障状态进行分类，进而实现诊断识别。首先以上述24组四种工作状态的IMF1与IMF2加和的多尺度熵作为训练集，并做出标记，令正常、内圈故障、滚动体故障以及外圈故障分别为1、2、3、4。然后在实验室的滚动轴承数据中，在不含以上数据的数据段中，把随机采取四种状态的6000点作为测试集，测试结果显示，利用改进的经验模态分解得到的IMF1与IMF2加和，其再进行多尺度熵形成的数组，利用LIBSVM之后可以精确的将新数据的四种情况进行分类，分别是正常为1、内圈为2、滚动体为3和外圈为4，从而达到了诊断目的。

4 结论

对于滚动轴承故障的诊断识别，经验模态分解（EMD）和多尺度熵（MSE）分别的对其进行了计算，识别效果有待改善。本文提出的通过改进的经验模态分解与多尺度熵结合再利用某一尺度下得到的固有模态函数熵值进行诊断识别，这种方法，通过实验，不仅有效的提取不同状态下故障特征的多尺度固有模态函数熵值，同时利用这组数值，结合LIBSVM对滚动轴承的故障进行诊断识别，效果显著。

[1]王书涛,李梅梅,张淑清,等.基于多重分形去趋势波动的机械故障诊断新方法[J].计量学报,2012,33(3):232-235.

[2]张超,陈建军,郭迅.基于能量熵和支持向量机的齿轮故障诊断方法[J].振动与冲击,2010,29(10):216-220.

[3]姚成,吴小培.小波变换与生物医学信号处理[J].生物学杂志,2000,17(1):24-29.

[4]Huang N E,Shen Z,Long S R et al.The empirical mode decomposition method and the hilbert spectrum for non-linear and non-stationary time series analysis[J].Royal Society of London Proceedings,1998,A(454):903-995.

[5]S.R.Qin,Zhong Y M.A new envelope algorithm of Hilbert-Huang Transform[J].Mechanical Systems and Signal Processing,2006,20(8):1941-1952.

[6]Delechelle E,Lemoine J,Niang O.Empirical Mode Decomposition:An analytical approach for sifting process[J].IEEE Signal Processing Letters,2005,11(12):764-767.

[7]Richman,Joshua S,Moorman JR.Physiological time-series analysis using approximate entropy and sample entropyAm J Physio[J].Heart Circ Physiol,2000,278(6):H2039-H2049.

[8]Costa M,Goldberger A L,Peng C K.Multiscale Entropy Analysis of Complex Physiologic Time Series[J].Physical Review Letters,2002,89(6):068102-1-068102-4.

[9]Lin J L,Liu J Y C,Li C W.Motor Shaft Misalignment Detection Using Multiscale Entropy with Wavelet Denoising[J].Expert Systems with Applications,2010,37(10):7200-7204.

[10]Litak G,Syta A,Rusinek R.Dynamical Changes during Composite Milling:Recurrence and Multiscale Entropy Analysis[J].The International Journal of Advanced Manufacturing Technology,2011,56(5-8):445-453.