EMD端點效應的改進型混沌延拓方法及其在機械故障診斷中的應用

蔡艷平，李艾華，石林鎖，何艷萍，趙靜茹

（第二炮兵工程學院五系，西安 710025）

1998 年，Huang［1］提出了基于經驗模態分解的EMD算法，EMD算法和與之相應的Hilbert譜統稱為Hilbert-Huang變換。Hilbert-Huang變換方法在處理非平穩、非線性信號方面具有獨特的優勢，但此方法在實際應用中還存在一些問題，其中端點效應引起的問題最為突出［2－12］。目前已有許多學者開展這方面的研究，并提出了不同的邊界延拓方法。通過分析現有解決EMD端點效應問題的研究成果，本文將其歸納為5類方法：

（1）極值點預測法，如陳忠等人［3］提出了通過添加極值點抑制端點效應的方法，劉慧婷等人［4］提出了基于多項式擬合的端點效應處理方法等；

（2）信號序列的鏡像延拓法，如趙進平等人［5］提出鏡像延拓算法，法國學者Rilling等人［6］提出的邊界極值點鏡像延拓方法等；

（3）信號序列預測法，如鄧擁軍等人［7］提出用神經網絡技術對數據序列進行延拓，楊建文等人［8］提出基于AR模型的時間序列邊界延拓方法，程軍圣等人［9］提出了基于支持矢量回歸機技術的邊界延拓方法等；

（4）波形延拓法，如Huang等人［10］提出了在信號兩端分別添加一組特征波的邊界延拓方法，并將該方法申請了美國的專利。胡愛軍等人［11］提出了波形特征匹配延拓法等；

（5）窗函數法，如Qi等人［12］提出用窗函數法對數據序列進行延拓，主要是在信號中間加矩形窗，在端點附近加漢寧窗、海明窗或者余弦窗，把窗函數和原信號相乘后的信號再做EMD分解。

上述這些方法的提出，為EMD的工程應用奠定了基礎，但這些方法均各有利弊。如極值點預測法，僅以信號端點附近極值點的特征信息決定延拓極值信息，因此，對非平穩信號有時會不能準確反映其真實趨勢，無法達到理想的EMD分解結果；鏡像延拓算法在信號具有較強的不對稱性時，則無論把鏡子放在哪里都不可避免的引入人為影響；支持矢量回歸機的邊界延拓方法具有較好效果，但其預測精度受內積函數、損失函數、核函數、精度參數和懲罰參數經驗選擇的影響；波形特征匹配延拓法，該法更適用于周期信號或準周期信號，對于非線性、非平穩信號具有局限性。窗函數法由于加上了一個窗，所以會改變原信號。因此，如何根據機械振動信號的非線性、非平穩特點，有針對性地提出適用性強的邊界延拓新方法，更好地將EMD應用于機械振動模式分析中，是值得進一步研究的。針對EMD端點延拓問題，本文提出了一種基于改進型最大Lyapunov指數的邊界延拓算法。

1 EMD端點效應問題機理分析

在標準EMD中，EMD包絡曲線擬合是采用三次樣條插值求取的，即先找出信號的所有極大值和極小值點，然后分別用三次樣條插值法來對極大值點、極小值點序列插值，從而獲得上、下包絡，再由上、下包絡線求信號的局部均值，以篩出信號的IMF。為了分析EMD端點效應問題，本文首先分析三次樣條插值給EMD包絡曲線擬合帶來的影響。

設插值對象為平面上n+1個點（xi，yi）（i=0，1，…，n），插值點x0＜x1＜x2＜…＜xn；三次樣條插值函數s（xi）滿足以下3個條件：

（1）s（xi）=yi，i=0，1，…，n；

（2）s（xi）在每個區間［xi－1，xi］（i=0，1，…，n）上是一個三次多項式：

si（x）=ai（x－xi）3+bi（x－xi）2+ci（x－xi）+di，i=0，1，…，n；

（3）s（xi）在整個區間［x0，xn］上有連續的一階和二階導數。

s（xi）還應滿足下列三種邊界條件之一：

（1） 已知s'（x0）=y'0，s'（xn）=y'n；

（2） 已知s″（x0）=y″0，s″（xn）=y″n；

其中第一種和第二中邊界條件還可以互相搭配產生新的邊界條件，則滿足上述條件的s（x）存在且唯一：

其中：hi=xi－xi－1，x∈［xi－1，xi］，i=0，1，…，n

由此可知，對于［x0，xn］以外的外插值點，三次樣條插值方法并未給出明確的規定。一般來說，三次樣條插值方法對［x0，xn］以外的插值點的插值只是邊界點處相鄰二點的三次多項式的延伸。

事實上，在EMD分解中，由于信號的兩個端點不一定是極值點，無法提供三次樣條插值函數所需要條件，所以包絡線在端點會產生發散出現，并且這種發散的結果會逐漸向內“污染”整個數據序列而使得結果嚴重失真，數據越短，其影響越大。

為了便于分析比較現有EMD端點效應解決方法，采用文獻［9］、［13］中使用的仿真信號來分析，該信號的表達式為：

仿真信號的采樣頻率為10 Hz，數據長度為100，EMD的端點效應現象見圖1。從圖1可見，該包絡線在數據兩端出現嚴重失真，尤其是上包絡線，發散現象非常嚴重。

圖1 三次樣條函數產生的端點效應Fig.1 The end effect of EMD generated by cubic spline function

2 現有EMD端點效應解決方法分析

本文主要分析目前三種典型邊界延拓算法：基于鏡像延拓的EMD改進算法、基于SVM延拓的EMD改進算法以及基于波形匹配延拓的EMD改進算法。圖2給出了鏡像延拓及其上下包絡線，圖3給出了SVM延拓及其上下包絡線，圖4給出了波形匹配延拓及其上下包絡線。

圖2 鏡像延拓及其上下包絡線Fig.2 Mirror extension and its upper and lower envelope of EMD

圖3 SVM延拓及其上下包絡線Fig.3 SVM extension and its upper and lower envelope of EMD

圖4 波形匹配延拓及其上下包絡線Fig.4 Waveform matching extension and its upper and lower envelope of EMD

從圖2可見，鏡像延拓方法對x（t）兩端的延拓改善了三次樣條包絡，但與真實樣條包絡線相比仍出現了失真，尤其是其上包絡線，在信號兩端的誤差很大，這勢必將會影響EMD分解的精度；從圖3可見，SVM延拓方法的包絡曲線與真實樣條包絡比較吻合，與鏡像延拓方法相比改善效果更好一些，但其上包絡線與真實樣條包絡線相比也出現了部分失真；從圖4可見，波形匹配延拓方法對x（t）兩端的延拓明顯改善了三次樣條包絡，且與真實樣條包絡線的變化趨勢一致，這說明了波形匹配延拓的有效性。顯然，只要原始信號存在一定的規律性，該方法所延拓數據就能延續原始信號的規律。但在實際信號中，如果信號內部沒有任何一段波形的變化趨勢和邊緣處相似，如無規律的強隨機性非平穩信號，則該方法的優勢將不存在。基于上述分析，本文結合混沌相空間重構理論，提出一種基于改進型Lyapunov指數邊界延拓的EMD端點效應處理方法。

3 EMD端點延拓的Lyapunov指數預測模型

由Packard和Takens提出的相空間重構理論，將混沌理論引入到時間序列分析中，系統任意分量的演化是由與之相互作用的其他分量所決定的，因此可以從某一時間序列中提取和恢復出系統原來的規律。對于特定的時間序列：x1，x2，…，xn，如果能找到合適的嵌入維數m和延遲時間τ，則按照Takens定理可以重構其相空間為：

其中：t=（m－1）τ+1，（m－1）τ+2，…，n，N=n－（m－1）τ為相點個數。

3.1 Lyapunov指數預測模型

Lyapunov指數刻畫了系統在相空間中的線度指數變化率，因此可將Lyapunov指數作為單序列變量的預報參數。設Lk為某時刻兩最近鄰點間的距離，經步長Tk后距離演化為L'k+1，并假設在演化過程中L'k+1/Lk近似為一常值，則：

設預測點Y0（t）及其最近鄰點Ys（t）經時間T后演化為Y0（t+T）和Ys（t+T），則可得預測方程：

當提前預報時間T≤τ時，Y0（t+T）中只有一個分量x0（t+T）是未知的，解預測方程可得：

其中：

3.2 Lyapunov指數預測模型的改進

利用最大Lyapunov指數法進行預測時，最大Lyapunov指數表示系統全體軌道在無窮多步演化條件下的平均發散度，因此，Lyapunov指數預測模型僅僅是軌道［Y0（t），Ys（t）］向［Y0（t+1），Ys（t+1）］演化規律的一個近似模型，而真實模型中應該用軌道［Y0（t），Ys（t）］向［Y0（t+1），Ys（t+1）］演化時的發散度代替最大Lyapunov指數，本文對Lyapunov指數預測模型進行改進，步驟如下：

Step1：用互信息法計算時間序列的最佳延遲時間τ，然后再由Cao氏法計算最佳嵌入維數m，重構相空間；

Step2：以局域發散度代替最大Lyapunov指數（基本預測模型同式5）；

局域發散度的計算方法：設參考點為Y0（t），其k個近鄰點為Yss（t），ss=1，2，…，k，Yss（t+T）表示近鄰點向前演化T（T≤τ）步之后對應的相點。對任意近鄰點對（i，j），i，j=1，2，…，k且i＜j，計算：

在查詢近鄰點時不限制短暫分離，考慮夾角，且近鄰點個數較少，目的是使近鄰點之間的關系能夠近似參考點與其最近鄰點之間的關系。

4 改進Lyapunov預測模型在端點延拓中的仿真

為了檢驗本文方法的有效性，分別采用鏡像延拓的EMD、波形匹配延拓EMD和基于改進型最大Lyapunov指數延拓EMD對式（2）給出的仿真信號進行EMD分解，圖5～圖7分別給出各自的EMD分解結果。

圖5 基于鏡像延拓的EMD分解Fig.5 The EMD decomposition result based on mirror extension

圖6 基于波形匹配延拓的EMD分解Fig.6 The EMD decomposition result based on waveform matching extension

圖7 基于本文方法的EMD分解Fig.7 The EMD decomposition result based on the proposed method

從圖5～圖7可以看出：基于改進型最大Lyapunov指數的邊界延拓效果最好，具有較好的端點效應抑制能力。相對來說，基于波形匹配延拓方法和鏡像延拓法都抑制了端點效應，但二者之間，基于波形匹配延拓方法的結果更勝一籌；基于鏡像延拓方法的效果要稍差一點，端點效應仍有較多誤差傳遞至信號內部。

5 旋轉機械故障診斷實例

為了檢驗本文方法在實測信號中應用的有效性，在WS－ZHT1進行模擬油膜渦動故障試驗，圖8為轉子試驗臺測試系統圖。

一般，當轉子的轉速到達某一轉速時，會產生油膜渦動現象，油膜渦動的振動機理如下：軸頸在軸承中轉動時會受到油膜的擠壓力，當擠壓力的切向分量大于阻尼力時轉軸將產生渦動，轉子繞其軸心線高速旋轉，轉速為V，當V大于失穩轉速時軸頸出現渦動，渦動頻率fr約為fV/2，隨著轉速V上升，渦動頻率fr也近似保持這個關系，這種低頻振動稱為半速渦動，其振幅一般不是很大。油膜渦動的振動特點為油膜渦動的頻率接近于轉子轉速的一半；當轉速上升到二倍于軸系的第一階臨界轉速時，油膜渦動轉變為油膜振蕩，表現為轉子的振幅增大，轉速失穩，且振蕩頻率不再追蹤轉速的變化，穩定在某一特定的值（軸系第一階臨界轉速）［14－15］。

圖8 轉子試驗臺測試系統圖Fig.8 Rotor test system

轉子油膜渦動的振動時域波形圖、頻譜圖如圖9、圖10所示。從頻譜圖10可看出，振動以半頻為主，頻率約為54 Hz，為轉子回轉頻率108.3Hz 的一半左右，是典型的油膜渦動故障。

為了檢驗本文方法在提取旋轉機械故障特征的有效性，限于篇幅，這里用鏡像邊界延拓、波形匹配延拓和最大Lyapunov指數邊界延拓三種改進EMD方法分別對轉子油膜渦動故障振動信號進行分析，見圖11～圖13。

圖11 鏡像邊界延拓方法分解結果Fig.11 The EMD decomposition result based on mirror extension

在圖11中，由上至下依次為原始信號和分解得到的3個IMF分量，第一個內稟模態函數（IMF1）和第二個內稟模態函數（IMF2）的端點產生了一定的發散，使得信號分解結果失真。由于鏡像邊界延拓方法存在的缺陷，致使EMD分解的結果中出現偽分量IMF3。在圖12中，由上至下依次為原始信號和分解得到的2個IMF分量。其中第一個內稟模態函數（IMF1）和第二個內稟模態函數（IMF2）分別對應了信號中1倍頻及0.5倍頻振動模式，但第一個內稟模態函數（IMF1）和第二個內稟模態函數（IMF2）的端點仍有部分失真。在圖13中，IMF1和IMF2分別對應了信號中1倍頻及0.5倍頻振動模式，且信號兩端基本沒有出現任何的擺動，分解結果非常理想，其在端點處的變化趨勢和原始信號的吻合度高，從高頻到低頻的2個IMF分量端點效應得到了較好地抑制，包絡線的形式和幅值都較為準確，并體現了原始信號的初始相位信息，上述對比說明本文提出的改進型Lyapunov指數邊界延拓方法很好地抑制了EMD的端點效應。

為更準確地反映出各EMD改進方法所得的IMF特征信息，對上述各EMD改進方法的IMF分別進行HHT時頻分析。圖14為用本文改進的EMD方法得到的HHT時頻圖，圖15為用本文改進的EMD方法得到的邊際譜圖。從圖14中可以看出：各譜線兩端沒有受到端點效應的影響，且時頻譜圖上主要頻率分布為1倍頻及0.5倍頻成分，同時1倍頻有明顯的調頻現象。從圖15中可以看出：0.5倍頻成分已高于1倍頻，與圖10中的幅值譜相比，半頻成分能量更為集中，幅值比基頻大得多，基頻成分準確表示了頻率調制范圍，在圖10中基頻成分頻率調頻特征基本沒有顯示出來，說明本文改進HHT的時頻分析結果能更好地體現油膜渦動的故障特征。

圖16為端點不作處理的EMD方法得到的HHT時頻圖。從圖中可以看出，其HHT時頻圖基本看不出0.5倍頻振動模式，且在l倍頻信號兩端產生了嚴重的端點效應現象。

圖17為用鏡像延拓方法得到的HHT時頻圖。從圖中可以看出，基于鏡像延拓方法的HHT時頻圖基本上包含了l倍頻和0.5倍頻振動模式，但波形存在失真，且在l倍頻信號兩端產生了一定的端點效應現象。

圖18為用波形匹配延拓方法得到的HHT時頻圖。從圖中可以看出，基于波形匹配延拓方法的HHT時頻圖也基本上包含了l倍頻和0.5倍頻振動模式，但效果不及改進型最大Lyapunov的邊界延拓方法，原因是其0.5倍頻振動模式波形存在部分失真。

圖18 基于波形匹配延拓方法的HHT時頻圖Fig.18 HHT time-frequency figure based on waveform matching extension method

對比圖14、圖16、圖17和圖18可看出四種方法分解結果的差異所在。顯然，本文方法使端點處的延拓更加合理，所延拓數據更加趨于真實。因此，基于本文改進的HHT方法更加準確地提取了轉子油膜渦動的故障特征，也進一步說明本文方法的有效性。

6 結論

在HHT應用中，端點效應是必須解決的問題之一，特別是在短數據序列的分析中，端點效應成為影響該方法精度的主要因素。本文在總結現有邊界處理方法特點的基礎上，結合混沌相空間重構理論，提出一種基于改進型Lyapunov指數邊界延拓的EMD端點效應處理方法。其基本原理就是按混沌動力學過程，通過計算時間序列的最大Lyapunov指數，建立混沌預測模型來預測時間序列端點的發展規律，使EMD篩分過程中不會因為端點原因而污染內部數據。在采用最大Lyapunov指數法進行邊界延拓時，針對Lyapunov指數預測模型所表示的平均發散度僅是軌道演化規律的一個近似模型，提出采用局域發散度代替最大Lyapunov指數進行改進。通過與現有典型邊界延拓方法的仿真對比表明，所提方法的包絡線優于其它方法，由于引入了時間序列混沌預測模型使端點處的延拓更合理，所延拓數據更真實；最后將所提方法應用在轉子油膜渦動故障診斷中，結果表明該方法能更準確地提取轉子油膜渦動故障特征。

需要指出的是，本文所提方法對于無規律具有混沌特征的非平穩信號，在減小HHT的端點效應方面相對其它方法來說具有優勢，但由于算法具有一定復雜性，其在計算時間上與神經網絡預測算法相當。EMD端點效應是EMD算法中一個非常棘手的問題，這個問題沒解決好，EMD算法將毫無用武之地。因此，在實際應用中應視被分析信號的需要來權衡計算精度和計算時間上的利弊進行選擇不同的端點效應解決方法。

［1］ Huang N E，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Procedures of the Royal Society of London，Series A，1998，454：903 －995.

［2］沈 路，周曉軍，張志剛，等.Hilbert-Huang變換中的一種端點延拓方法［J］.振動與沖擊，2009，28（8）：168－171，174.

［3］陳 忠，鄭時雄.EMD信號分析方法邊緣效應的分析［J］.數據采集與處理，2003，18：1114－1118.

［4］劉慧婷，張 旻，程家興.基于多項式擬合算法的EMD端點問題的處理［J］.計算機工程與應用，2004，（16）：84－86.

［5］ Zhao J P，Huang D J.Mirror extending and circular spline function for empirical mode decomposition method［J］.Journal of Zhejiang University（Science），2001，2：247－252.

［6］ Rilling G，Flandrin P，Goncalves P.On Empirical Mode Decomposition and Its Algorithms［A］.In IEEE-Eurasip Workshop on Nonlinear Signal and Image Processing，NSIP－03［C］，Grado（I）.2003.

［7］Deng Y J，Wang W，Qian C C.End-processing-technique in EMD method and Hilbert transform［J］.Chinese Science Bulletin，2001，46（11）：954－961.

［8］楊建文，賈民平.希爾伯特－黃譜的端點效應分析及處理方法研究［J］.振動工程學報，2006，19（2）：283－288.

［9］ Cheng J S，Yu D J，Yang Y.Application of support vector regression machines to the recessing of end effects of Hilbert-Huang transform［J］.MechanicalSystemsand Signal Processing，2007（21）：1197 －1211.

［10］ Huang N E.Empirical mode decomposition apparatus，method and article of manufacture for analyzing biological signals and perfoming curve fitting［P］.United States Patent，2002，6：381.

［11］胡愛軍，安連鎖，廖貴基.Hilbert-Huang變換端點效應處理新方法［J］.機械工程學報，2008，44（4）：154－158.

［12］ Qi K Y，He Z J.Cosine window-based end processing method for EMD and its application in rubbing fault diagnosis［J］.Mechanical Systems and Signal Processing，2007（21）：2750－2760.

［13］程軍圣，于德介，楊 宇.Hilbert-Huang變換端點效應問題的探討［J］.振動與沖擊，2005，24（6）：40－42.

［14］馬 輝，李朝峰，軒廣進，等.轉子系統油膜失穩故障的時頻特征分析［J］.振動與沖擊，2010，29（2）：193－195，198.

［15］曹沖鋒，楊世錫，楊將新.一種基于瞬時能量分布特征的汽輪發電機組轉子故障診斷新方法［J］.振動與沖擊，2009，28（3）：35－39.