基于EEMD的振動信號自適應(yīng)降噪方法

陳仁祥，湯寶平，馬婧華

(重慶大學(xué) 機械傳動國家重點實驗室，重慶 400030)

機械設(shè)備發(fā)生早期故障時，故障信號一般很弱，常常淹沒在噪聲信號中。因此，為了準確提取故障信息，需要對原始信號降噪，提高信噪比。由于機械設(shè)備振動信號的非高斯、非平穩(wěn)特性，基于傳統(tǒng)傅里葉變換的降噪方法存在保護信號邊緣和抑制噪聲之間的矛盾，難以正確識別信號中的噪聲并加以去除。基于小波變換的降噪方法［1－2］，對非平穩(wěn)信號降噪，要比傳統(tǒng)的濾波降噪方法效果好，但這種方法存在選擇小波基和確定閾值等問題［3－4］。

Huang等［5］提出一種新的非平穩(wěn)信號處理方法——經(jīng)驗?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)，EMD是一種完全基于數(shù)據(jù)驅(qū)動的自適應(yīng)分解方法，能將信號從高頻到低頻分解成有限個具有物理意義的固有模式函數(shù)(Intrinsic Mode Function，IMF)和余項之和。文獻［6］進一步分析了EMD方法對高斯白噪聲的分解，得出EMD可以作為二進濾波器的結(jié)論。然而，EMD本身存在一些不足，如模式混疊、端點效應(yīng)、停止條件等［7－8］。模式混疊是指1個IMF中包含差異極大的特征時間尺度，或者相近的特征時間尺度分布在不同的IMF中［10］，導(dǎo)致相鄰的2個IMF波形混疊，相互影響，難以辨別。Huang［5］認為引起模式混疊的原因在于間歇現(xiàn)象(intermittency)，引起間歇現(xiàn)象的往往是異常事件(如間斷信號，脈沖干擾和噪聲等)。因此，模式混疊現(xiàn)象使EMD對含有異常事件的信號的降噪效果不佳［9－11］。

為了抑制模式混疊，Wu等［12］提出一種集合經(jīng)驗?zāi)Ｊ椒纸夥椒?，該方法將噪聲輔助分析應(yīng)用于經(jīng)驗?zāi)Ｊ椒纸庵校源龠M抗混分解，有效的抑制了混疊現(xiàn)象。文獻［11］將EEMD應(yīng)用于對大型旋轉(zhuǎn)機械啟動過程中振動信號的降噪，得出了EEMD降噪既可以彌補小波降噪方法對調(diào)頻調(diào)幅信號處理過程中存在的特征波匹配缺陷，又可以克服EMD降噪方法對脈沖干擾下振動信號濾波能力不足的結(jié)論。該文的方法是，對信號進行EEMD分解后計算IMF分量與原始信號互信息值來與事先確定的閾值對比以選取IMF分量重構(gòu)信號，從而達到降噪的目的。但在該文中未明確說明閾值的確定方法。文獻［13］將EEMD方法用于對疲勞應(yīng)變信號降噪，其方法是對信號做EEMD分解后，選取IMF分量來重構(gòu)信號以對疲勞信號降噪，但文中也未說明IMF分量的選取方法。

針對應(yīng)用EEMD降噪時IMF分量選取問題，本文從白噪聲經(jīng)過EMD分解后的IMF分量特性出發(fā)，設(shè)計了自動選取IMF分量重構(gòu)信號的算法，提出了基于EEMD的振動信號自適應(yīng)降噪方法。

1 EEMD的基本原理及抗混分析

1.1 模式混疊產(chǎn)生分析

EMD方法自適應(yīng)的將任意信號的分解成有限個IMF分量和余項之和，每個IMF需要滿足兩個條件［5］:① 在整個數(shù)據(jù)集中，極值點和過零點交替出現(xiàn)且極值點個數(shù)和過零點數(shù)相等或最多相差1個;② 在任意時刻，由局部極大值點形成的上包絡(luò)線和由局部極小值點形成下包絡(luò)線的平均值為零，即上、下包絡(luò)線相對于時間軸局部對稱。文獻［5］中給出了EMD算法的詳細步驟。

EMD過程中首先需要確定信號的局部極值點，然后用三次樣條線將所有的局部極大值和極小值點分別連接起來形成上下包絡(luò)線，再由上下包絡(luò)線得到均值曲線。在求取包絡(luò)線的過程中，當(dāng)信號中存在異常事件時，勢必影響極值點的選取，使極值點分布不均勻，從而導(dǎo)致求取的包絡(luò)為異常事件的局部包絡(luò)和真實信號包絡(luò)的組合。經(jīng)該包絡(luò)計算出的均值，再篩選出的IMF分量就包含了信號的固有模式和異常事件或者包含了相鄰特征時間尺度的固有模式，從而產(chǎn)生了模式混疊現(xiàn)象。

1.2 EEMD基本原理與方法

從前面的分析可知，極值點分布的不均勻?qū)е铝四Ｊ交殳B現(xiàn)象。為此，Wu等［12］將白噪聲加入待分解信號來平滑異常事件，利用白噪聲頻譜的均勻分布來使不同時間尺度的信號自動分布到合適的參考尺度上。同時，運用白噪聲的零均值特性，經(jīng)過多次平均使噪聲相互抵消，從而抑制甚至完全消除噪聲的影響?；谶@種思想，文獻［12］提出了集合經(jīng)驗?zāi)Ｊ椒纸狻?/p>

顯然，集合經(jīng)驗?zāi)Ｊ椒纸獾谋举|(zhì)是一種疊加高斯白噪聲的多次經(jīng)驗?zāi)Ｊ椒纸狻?/p>

EEMD的分解步驟如下:

第一步 在原始信號x(t)中M次(M＞1)加入均值為0，幅值標準差為常數(shù)的高斯白噪聲ni(t)(i=1～M)，即:

第二步 對xi(t)分別進行EMD分解，得到的K個IMF記為cij(t)(j=1～K)，余項記為ri(t)。其中，cij(t)表示第i次加入高斯白噪聲后，分解所得到的第j個IMF。

第三步 利用不相關(guān)隨機序列的統(tǒng)計均值為0的原理，將以上步驟對應(yīng)的IMF進行總體平均運算，以消除多次加入高斯白噪聲對真實IMF的影響，最終得到的EEMD分解后的IMF為:

式中:cj(t)表示對原始信號進行EEMD分解后所得到的第 j個 IMF。

1.3 EMD與EEMD的比較

從EMD和EEMD的原理和算法可知，EEMD本質(zhì)上是EMD的改進，它不僅具有EMD的優(yōu)點而且可以有效抑制EMD分解過程中的模式混疊現(xiàn)象。下面通過仿真實例來說明EEMD抑制模式混疊的能力。

文獻［5］提出了對于信號x(t)的正交性指標(Index of Orthogonality，IO)來評價EMD的精度，IO值越小則精度越高。IO可定義為［5］:

仿真信號s(圖1(d))由Gauspuls脈沖分量干擾s1(圖1(a))、頻率為7Hz的正弦分量s2(圖1(b))和趨勢項s3(圖1(c))組成，信號長度為1024點。

圖2為仿真信號的EMD分解結(jié)果，包括三個IMF分量(c1～c3)和一個余項res。顯然，只有c3和res才具有真實的物理意義，分別代表了原信號中的正弦成分和趨勢項。c1完全失真，失去了物理意義，多出來的c2其實是s1和s2的模式混疊。事實上，c1和c2中都分布有s1和s2的特征時間尺度，所以，在這種情況下EMD無法分析信號本質(zhì)，存在嚴重的模式混疊現(xiàn)象。圖3為EEMD分解結(jié)果，加入高斯白噪聲標準差為0.01，加入次數(shù)為100次。從圖3中可以看出，c1代表了s1，c3代表了s2，res代表了 s3。多出來的c2時間尺度介于c1和c3之間，是由于s1和s2調(diào)制所致。

圖1 仿真信號及其組成Fig.1 Simulation signal and its components

進一步比較EMD和EEMD方法計算效率和IO指標，如表1所示。顯然EEMD的耗時比EMD多，但其IO值得到了很大提高。綜上，EEMD有效地抑制了模式混疊現(xiàn)象，能高質(zhì)量的分解出原信號中的各個組分，比EMD更具優(yōu)勢。

表1 兩種分解方法結(jié)果比較Tab.1 A comparison between results of two decomposition methods

2 EEMD自適應(yīng)降噪原理與方法

EEMD把信號分解成特征時間尺度由小到大即頻率由高到底的一系列IMF分量。對混有隨機噪聲的信號，經(jīng)分解后的高頻IMF分量通常情況下為噪聲，將這些噪聲IMF分量去除，由剩下的IMF分量重構(gòu)信號即可進行降噪。這就是EEMD降噪的原理。

在使用EEMD降噪時需要解決兩個問題:① 應(yīng)用EEMD分解信號時，需要設(shè)置兩個參數(shù)，即加入原信號的白噪聲的幅值系數(shù)k和執(zhí)行EMD的總次數(shù)M;②重構(gòu)信號時，選取IMF分量的方法。針對問題①，根據(jù)文獻［12］，k用待分析信號的標準偏差乘以一個分數(shù)來定義比較合適，M取一、兩百次時，殘留噪聲引起的誤差為處在一個較低水平(小于1%)。經(jīng)過筆者的試驗分析，M 取100～300時，k值選擇0.01～0.5倍信號的標準偏差較為適宜，M、k值可隨噪聲的強度而適當(dāng)增大;針對問題②，由文獻［14］的研究可知，白噪聲序列經(jīng)EMD濾波分解后各IMF分量的能量密度與其平均周期的乘積為一常量，從該性質(zhì)出發(fā)，設(shè)計了一種自動選擇IMF分量的方法。具體步驟是:

第一步對含噪信號進行EEMD分解后得到K個IMF分量，計算出每個IMF分量的能量密度與其平均周期之積:

式中:Ej和為第j個IMF分量的能量密度和平均周期，其計算公式如式(5)和式(6):

式中:N為每個IMF的長度，Aj為第j個IMF分量的振幅，Oj為第j個IMF分量的極值點總數(shù)。

第二步根據(jù)式(7)計算系數(shù)RPj(j≥2)

當(dāng)RPj≥1時，可確定第j個IMF的Pj相對于前面j－1個IMF的Pj的平均值成倍增大，即認為前面j－1個IMF的能量密度與其平均周期之積為一常量。由上文可知，前面的j－1個IMF分量作為噪聲而去除，同時，余項作為趨勢項也應(yīng)去除，由剩下的IMF分量重構(gòu)信號即可對原含噪信號進行降噪。

綜上，該降噪方法計算流程如圖4所示。

圖4 降噪算法流程圖Fig.4 Flow chart of de-noising method

3 仿真試驗分析

仿真信號s_noise由前文中的仿真信號s加上標準差為0.4的白噪聲組成，信號長度為1024點，波形圖如圖5所示。

圖5 仿真信號s_noiseFig.5 Simulation signal s_noise

首先對 s_noise進行 EEMD分解，M=100，k=0.02，結(jié)果如圖6所示，一共分解出7個IMF分量和一個余項。根據(jù)本文提出的降噪算法流程，計算其RP值得出當(dāng)j=6時，RP=5.7401＞1，所以自動選取c6和c7重構(gòu)信號，降噪結(jié)果如圖7(b)所示。圖7(a)為EMD降噪結(jié)果，圓圈處出現(xiàn)明顯的擾動成分，這是由于EMD的模式混疊所引起。

根據(jù)式(8)計算出信噪比，結(jié)果如表2所示。

式中:Ps和Pn分別是信號和噪聲的有效功率。

從表2中可以看出，經(jīng)過降噪處理后，信號的信噪比得到了大幅度提高。同時，本文方法降噪后的信噪比更高，證明了本文提出的降噪方法是可行的，且更為有效。

表2 兩種降噪方法信噪比Tab.2 The SNR of two de-noising methods

4 工程應(yīng)用

將本文提出的降噪方法應(yīng)用于軸承振動數(shù)據(jù)的降噪。試驗臺由電動機、軸承、轉(zhuǎn)子、傳動軸、軸承支架、模擬載荷、電器控制裝置等組成，以輸入端軸承為測試對象。采集部分選用16位的信號采集儀和壓電加速度傳感器套裝，傳感器布置在軸承座的垂直方向。現(xiàn)場試驗圖如圖8所示。

軸承型號為UN205，其內(nèi)徑為25mm，外徑為52mm，滾動體直徑為7.5mm，滾動體數(shù)目為12，接觸角為0°。該軸承外圈存在故障，轉(zhuǎn)速為800 r/min，經(jīng)過計算其外圈故障頻率為64.4Hz。原始數(shù)據(jù)長度為1024點，采樣頻率為12 k，其波形如圖9(a)所示。受噪聲干擾影響，從圖9(a)的時域波形中觀察故障特征較困難。

圖9(b)為采用本文提出的降噪方法的結(jié)果，從圖中可以明顯看出沖擊的周期性，并且周期性脈沖衰減明顯。其沖擊周期T約為15.5 ms，其倒數(shù)即是信號特征頻率為64.5Hz，該頻率和前文計算出的軸承外圈故障特征頻率吻合，由此可判斷軸承出現(xiàn)了外圈故障，與實際相符。

5 結(jié)論

EEMD能有效的抑制EMD中的模式混疊現(xiàn)象，其分解精度也更高。在應(yīng)用EEMD降噪時，給出了EEMD中參數(shù)選取原則，并根據(jù)白噪聲經(jīng)EMD分解后其IMF的特性，設(shè)計了一種自動選取IMF分量重構(gòu)信號的算法，提出了基于EEMD的振動信號自適應(yīng)降噪方法。通過應(yīng)用本文方法對仿真信號和軸承振動降噪，驗證了該方法的可行性和有效性，為機械振動信號降噪提供了一種新方法。

［1］Liu F，Ruan X E.Wavelet-based diffusion approaches for signal de-noising［J］.Signal Processing，2007，87:1138－1146.

［2］Azzalini A，F(xiàn)arge M，Schneider K.Nonlinearwavelet thresholding:A recursive method to determine the optimal denoising threshold［J］.App.l Comput.Harmon.Ana，l 2005，18:177－185.

［3］張吉先，鐘秋海，戴亞平.小波門限消噪法應(yīng)用中分解層數(shù)及閾值的確定［J］.中國電機工程學(xué)報，2004，2(24):118－122.

［4］程發(fā)斌，湯寶平，鐘佑明.基于最有Morlet小波和SVD的濾波消噪方法及故障診斷的應(yīng)用［J］.振動與沖擊，2008，27(2):91－94.

［5］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of the Royal Society of London Series A-Mathematical Physical and Engineering Sciences，1998，454:903－995.

［6］Flandrin P，Rilling G，Goncalves P.Empirical mode decomposition as a filter bank［J］.IEEE Signal Processing Lett，11，2004:112 －114.

［7］Rilling G，F(xiàn)landrin P，Goncalyes P.On empirical mode decomposition and its algorithms［C］//IEEE-EURASIP Workshop on Nonlinear Signaland Image Processing(NSIPL2003)，Grado(1)，June，2003:8 －11.

［8］Rato R T，Ortigueira M D，Batista A G.On the HHT，its problems，and some solutions［J］.Mechanical Systems and Signal Processing，2008，22:1374 －1394.

［9］Norden E，Huang N E.A new view of non-linearwaves:The Hi-lbert spectrum［J］.Annual Review of Fluid Mechanics，1999，31:417－457.

［10］趙進平.異常事件對EMD方法的影響及其解決方法研究［J］.青島海洋大學(xué)學(xué)報，2001，31(6):805－814.

［11］曹沖鋒，楊世錫，楊將新.大型旋轉(zhuǎn)機械非平穩(wěn)振動信號的EEMD 降噪方法［J］.振動與沖擊，2009，28(9):33－38.

［12］WuZH， HuangN E. Ensembleempiricalmode decomposition:a noise assisted data analysis method［J］.Advances in Adaptive Data Analysis，2009，1(1):1 －41.

［13］陳 雋，李 想.運用總體經(jīng)驗?zāi)Ｊ椒纸獾钠谛盘柦翟敕椒ǎ跩］.振動、測試與診斷，2011，31(1):15－19.

［14］Wu Z，Huang N E.A study of the characteristics of white noise using the empirical mode decomposition method［J］.Proc.R.Soc.London A，2004，460:1597 －1611.