利用經驗模態分解和多尺度熵的滾動軸承故障診斷研究

秦喜文 ，高中華，董小剛，張邦成

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.