基于經驗模式分解和排列熵的軸承故障特征提取

王 濤，胡定玉，丁亞琦，廖愛華，師 蔚

（1.上海工程技術大學 城市軌道交通學院，上海201620；2.上海地鐵維護保障有限公司 車輛分公司，上海200235）

滾動軸承作為減少摩擦損失的精密元件，被廣泛用于旋轉機械中，但高速重載的環境導致軸承極易發生故障，影響列車的安全出行，因此對軸承的運行狀況檢測十分必要。當滾動軸承出現故障時，內外滾道與滾動部件之間發生碰撞，會產生一個強迫短脈沖信號，通過對這些信號進行分析可以有效檢測故障從而減少事故的發生。與傳統的振動分析相比，采用聲信號監測軸承的健康狀況更加方便快捷，無需粘貼傳感器，與檢測對象相互獨立，該方法對軸承早期故障的檢測也有較好的效果。

傳統的時域與頻域分析方法難以有效分析軸承故障信號的時變特性，無法準確體現其故障特征。針對該類信號，通常采用小波分解[1]和經驗模式分解[2]等時頻分析方法。時頻分析將時頻兩域結合起來，可以同時考慮到兩個域的性能；在時頻相平面上可以精確地定位出某一時刻出現了哪些頻率分量。小波分解在分析信號之前設置一些參數（小波基函數），將信號分解成不同頻率范圍的子信號[3-4]，不同的小波基會產生不同的分解結果，最佳小波基的選擇成為了較大的難題。經驗模態分解（EMD）[5]可以自適應地將信號分解為固有模態函數和殘差，可以對不同的分量進行獨立分析[6]。

目前EMD 已在滾動軸承故障診斷領域廣泛應用，楊世錫等[7]將經驗模式分解和Hilbert 變換引入旋轉機械的故障診斷中，并與小波時頻分析相比較，研究結果表明基于EMD 的Hilbert 變換具有更好的效果，但出現了一些虛假分量；張志剛等[8]提出一種改進經驗模態分解和譜峭度相結合的算法，剔除了分解過程中出現的虛假分量，但缺點是其受噪聲影響較大。

本文提出一種基于經驗模式分解和排列熵的軸承故障分析的信號處理方法，為EMD分解后IMF分量的選擇提供了新的方法。首先介紹相關基礎理論，然后通過數值仿真與傳統的包絡分析進行比較，最后通過實驗數據的分析驗證該方法的有效性。

1 基礎理論

1.1 經驗模式分解算法

EMD是一種自適應時頻信號處理方法，通過對復雜信號進行分解，獲得一系列用于表征原始信號特性的IMF 分量，篩選一些合適的本征模態分量實現對原始信號所包含的信息特征進行分析。分解出的IMF 必須滿足以下兩個條件[9-11]：(1)在所需分解的數據段中，極點數和過零點數的差值不能超過一個；(2)信號的任意時刻，由極大值點插值所構造的上包絡線和極小值點構成的下包絡線均值為零。

EMD的實施步驟如下：

(1)找出原始信號所有的極大值點和極小值點，利用插值的方法分別對極大值點和極小值點進行連接來擬合原始信號的波形，形成上包絡線和下包絡線。

(2)計算上下包絡線的平均值m1，求得：

若h1滿足上述條件，則為信號的一個IMF分量，記為IMF1，進入下一個循環。

(3)若h1不滿足條件，則用h1替代初始數據x(t)，重復步驟(a)、(b)，得到h2，并判斷h2是否滿足(b)的條件，若還是不滿足條件，則繼續循環迭代直到獲得滿足條件的IMF1分量。

(4)將IMF1分量從信號x(t)中分離出來，得到新的循環信號：

(5)將R1作為原始信號，重復步驟(1)～(5)，當Rn的極值點數小于2，則循環結束，得到n個IMF 分量和殘余分量Rn，信號x(t)可表示為

在EMD的分解過程中，由于需要人為設置迭代次數和終止條件，因此會產生一些無關信息成分，即虛假分量。針對該問題，可以利用相關系數來識別虛假分量[12]。

1.2 排列熵算法

排列熵作為描述復雜系統的算法，對突變的信號較為敏感。當滾動軸承系統發生故障時，聲信號包含的沖擊成分，會對排列熵值造成影響，因此可以利用排列熵對軸承進行早期故障診斷。

給定一個序列{X(K),K=1,2,…,N}，對其進行相空間重構得：

其中：d為嵌入維數，τ為延遲時間，k為重構分量。對矩陣Z中的每一行序列進行升序排列：

式中：t表示索引的列數，R1,R2,…,Rd為X(K)中各元素所處的位置。定義xdi為重構后的任意一組序列，0 ≤Ri≤d!，對于d維相空間映射有d!種排列的可能性，計算每一種序列出現的可能性P1,P2,…,Pi，對于時間序列X(K)的i種排列熵Hp(d)：

當Pi=m!時，Hp(d)就達到最大值In(m!)，為了方便熵值的比較，往往會對排列熵值進行歸一化處理：

Hp的值越小則表示該序列越有序，反之該序列的復雜程度越高。由于發生故障的軸承元件在運行過程中與其他元件發生周期性的沖擊，相較于正常軸承信號，故障信號更加有序，排列熵值較小。因此可以用排列熵值的突變來判斷軸承是否發生故障。

在計算排列熵時，插入維數d過小會導致重構信號狀態量缺失，對沖擊信號的敏感性將大大下降，算法失去有效性；插入維數過大會使得信號在相空間的重構過程中均勻化時間序列，因而一般維數選擇為3～7，本文取5，時間延遲τ=1[13]。

EMD將信號分解為多個IMF分量，高頻IMF分量中包含較多沖擊成分，但包含大量噪聲，低頻IMF分量噪聲成分較少，但可能是虛假分量。為篩選本征模態分量中包含信息量較大的真實分量進行信號的重構，可以利用排列熵值和相關系數來作為篩選指標。排列熵值Hp與原始信號的有效成分呈反相關，相關系數r與原始信號的線性相關程度呈正相關。基于此本文提出聯合系數Co作為IMF 分量的篩選指標，取聯合系數最大的前兩個本征模態分量進行信號的重構：

1.3 快速譜峭度算法

譜峭度為能量歸一化的4 階譜累積量，利用計算每根譜線的峭度值來反映瞬時沖擊的大小，快速譜峭度(Fast-kurtogram)則以STFT作為窗口，根據快速譜峭度圖選擇最佳濾波參數。

在信號x(t)的激勵下系統響應y(t)Wold-cramer分解的頻域表達式為

其中：H(t,f)為y(t)在頻率f處的復包絡信號。y(t)的四階譜累積量為

式中：S2ny(f)為2n階譜瞬時距，用以表示復包絡信號的能量大小，定義為

式中，df表示對頻率f求導。譜峭度定義為

為減少運算時間，可以應用于實際的工程測量，采用塔式快速譜峭度對信號進行2倍頻采樣濾波迭代，將原信號進行L層濾波分解，每層濾波后的信號頻段為2L個，一般情況下L為整數。為了提高信號分解的精度，在L+1 層和L+2 層中插入3×2L個濾波器，即L+1.6層，如當L=0時，第1層和第2層之間分解層數為1.6層時，表明該層將原始信號分解為3個信號頻段。根據快速譜峭度圖上峭度值的大小，即顏色最深的頻段選取最佳濾波中心和濾波帶寬，進行平方包絡分析和故障提取。

2 EMD-PE故障診斷方法流程

應用快速譜峭度算法診斷軸承故障雖然取得了一定的效果，但對信號進行預處理可以進一步提高診斷效果，本文將EMD 和排列熵值相結合，提出一種基于聲信號滾動軸承故障診斷的EMD-PE 方法，可以有效降低噪聲的干擾并確定合適的濾波器參數。該方法流程圖如圖1所示。具體過程如下。

(1)對原始信號進行EMD分解，把非平穩、非線性的聲信號分解為8 個本征模態分量IMF1，IMF2，···，IMF8。

圖1 EMD-PE診斷故障軸承流程圖

(2)計算各個本征模態分量的排列熵值Hp1，Hp2，··· ，Hp8(0 ＜Hpi＜1)和相關系數r1，r2，··· ，r8(0 ＜ri＜1)；

(3)計算聯合系數，根據聯合系數最大原則，篩選系數值最大的前兩個本征模態分量進行信號的重構。

(4)對重構信號進行快速譜峭度分析，根據峭度最大化原則選擇最佳濾波中心和濾波帶寬。

(5)將從中提取到的軸承特征頻率與理論故障特征頻率對比，進行故障診斷。

3 軸承故障信號仿真分析

為了證明該方法在滾動軸承故障診斷中的有效性，現對模擬的軸承故障仿真信號進行分析，設軸承發生單點局部故障，傳聲器接收到信號的信號數學模型為

其中：x(t)為軸承故障仿真信號，A為沖擊信號幅值，s(t)為脈沖沖擊函數，Ti為特征頻率出現的周期，fs為軸承座、傳感器系統的高頻固有頻率，ξ為幅值衰減系數，n(t)為高斯白噪聲信號。

設仿真參數A=1，特征頻率fi=1/Ti=130 Hz，固有頻率fs=5 000 Hz，衰減系數ξ=1000，采樣頻率Fs=51200 Hz，采樣點數N=10 240。仿真脈沖信號如圖2(a)所示，對仿真信號添加高斯白噪聲，SNR=10，如圖2(b)所示。

對含有高斯白噪聲的仿真信號進行EMD分解，對分解后的各個IMF分量進行排列熵值與相關系數的計算，并根據式(8)求得聯合系數，如表1所示。根據聯合系數最大化原則，篩選出兩個最佳的IMF 分量(IMF6，IMF7)進行信號的重組，圖3顯示的是重組信號進行快速譜峭度分析得到快速譜峭度圖。由圖3可以看出，分解層數為1.6，帶寬為8 333.33 Hz，濾波中心為23 333.33 Hz時，譜峭度值達到最大。在濾波后的平方包絡圖中，準確檢測出了130 Hz 的特征頻率及其倍頻程如圖4所示，這與預設的故障特征頻率相符合，于是實現了故障的診斷，說明基于EMD 和排列熵的快速譜峭度的軸承故障診斷是有效的。

表1 仿真數據

圖2 仿真信號時域圖

圖3 仿真信號快速譜峭度圖

圖4 基于EMD-PE算法的軸承故障仿真試驗結果

4 實驗分析

實驗采用SKF 公司生產的滾動軸承故障實驗臺，故障軸承型號為SKF6016深溝球軸承,內圈存在一條裂紋，深約1 mm，其他部位保持完好，表2為該軸承的參數。采樣頻率是51 200 Hz，采樣時間為10s，電機轉速為404 r/min，轉頻為7 Hz，負載為5 N。

表2 滾動軸承參數

對于外圈固定，內圈轉動的滾動軸承，其內圈故障的特征頻率為

將軸承參數代入式(14)中可得軸承內圈故障特征頻率為55.69 Hz。

聲信號的希爾伯特包絡變換的頻譜圖如圖5所示，Hilbert頻譜圖中雖然可以得到故障頻率，但卻不夠明顯。

對信號進行EMD 分解，得到12 個本征模態分量，計算各個分量的排列熵值、相關系數以及聯合系數，如表3所示，由于篇幅有限，這里只顯示前8 個IMF的計算結果，其中原始信號的排列熵值為0.77。

由表3可知IMF4和IMF5的聯合系數Co較大，原始信號中的突變特性成分較多，對IMF4和IMF5分量進行重組，得到的重組信號的排列熵值為0.29，相比于原始故障信號有用信息比例有所提高。

圖5 轉速為404r min時聲信號的希爾伯特包絡變換頻譜圖

表3 各IMF分量排列熵值、相關系數及聯合系數

對合成信號進行快速譜峭度分析得到快速譜峭度圖和平方包絡圖,如圖6(a)、圖6(b)所示。選取分解層數為4，帶通濾波器的濾波中心18 000 Hz，而濾波帶寬為1 200 Hz,在此濾波范圍內峭度值達到最大，對濾波后的信號進行平方包絡分析可以得到較為明顯的滾動軸承內圈故障頻率53.93 Hz 及其倍頻，證明了本文方法的有效性。

圖6 重構信號的快速峭度譜圖和平方包絡譜圖

5 結語

本文提出了一種基于經驗模式分解和排列熵的軸承故障分析的信號處理方法，利用聯合系數代替了傳統的峭度值，對IMF 分量進行篩選重構。此外針對濾波參數較難確定的問題，譜峭度對聲信號中瞬態沖擊較為敏感，可以自動確定帶通濾波器的參數。經驗模式分解和排列熵相結合的滾動軸承故障診斷方法，能夠有效提取出故障特征頻率及其倍頻。仿真及實驗結果驗證了該方法的有效性。