郭明威,倪世宏,朱家海,張志鵬
(空軍工程大學 工程學院,西安 710038)
Hilbert-Huang變換(HHT)是 Huang等[1]提出的一種適應于非線性、非平穩信號的時頻處理方法。HHT的提出被認為是近年來以傅里葉變換為基礎的線性和穩態譜分析的一個重大突破,經過十幾年的發展,HHT理論體系不斷完善,已成功應用于地震信號、圖像處理、故障診斷等學科領域[2-14]。但在HHT應用中,端點效應問題直接影響到信號HHT分析結果,已成為阻礙其工程實用的難點問題[1-12]。
針對端點效應問題,國內外學者做了大量研究,目前提出的抑制端點效應的方法大致分為三類:波形延拓法[2-4]、數據預測延拓法[5-9]、極值延拓法[3,10-11]。在波形延拓方面,波形匹配延拓法[2]利用信號內部變化趨勢來延拓信號,結果更接近原始數據特征,該方法延拓后會出現數據“斷裂”現象,引入“異常事件”而影響分解效果;波形鏡像延拓法[3-4]在具有對稱性的極值所在位置放置鏡子,實現鏡內數據的閉合延拓,該方法對信號的對稱性要求較高,且沒有從數據上給出最佳的鏡面位置。在數據預測延拓方面,利用神經網絡建模對數據進行預測延拓[5],借以抑制端點效應,但該方法對非線性、非平穩信號延拓效果不佳;文獻[6]采用自回歸(AR)模型處理EMD方法中的邊界問題,取得了一定的效果,但AR模型只適應于平穩序列;文獻[7]采用時間序列建模與預測方法對信號數據進行延拓,對端點效應進行了探討,但對于非平穩性較強的振動信號預測效果不佳;基于支持矢量回歸機的端點效應問題處理方法[9],比神經網絡和AR模型預測效果好,但該法用于學習訓練的數據點多,預測步數過多,致使計算量龐大,并且延拓誤差隨多步預測而積累變大,預測精度下降。在極值延拓方面,包絡極值延拓方法[3]實質就是對極值點的等幅、等距延拓,不能反應信號趨勢,并且此方法未考慮到特殊端點的處理;端點優化對稱延拓方法[10]比ARMA延拓方法效果要好,但端點優化對稱延拓方法對端點附近數據處理不夠平滑,可能引起數據的階躍;基于相似極值延拓[11],其時間的延拓為等距延拓,延拓極值點的大小取值為已知極值點序列的均值,導致其自適應性不強,延拓效果不理想。
本文提出一種基于時間尺度的LS-SVM端點延拓方法,通過仿真信號分析,并將其應用到航空發動機的振動數據處理中,驗證了該方法能有效地抑制端點效應。
經驗模態分解(Empirical Mode Decomposition,EMD)和與之相應的Hilbert譜信號分析方法統稱為Hilbert-Huang變換。HHT的端點效應現象在EMD和Hilbert變換中都存在,這兩種失真的疊加,造成了HHT的端點部分無法正確反映信號所包含的信息。
EMD本質是通過數據的特征時間尺度來獲得固有波動模式,其基本思想是:對一給定信號,先獲得信號極值點,通過插值獲得信號包絡,得到均值,與均值的差得到分解的一層信號,如此重復,直到將信號分解成有限個基本模式分量(IMF)和殘余項rn(t)的組合,表示為:

EMD是基于數據極點信息的分解過程,或者說極點是信號的EMD分解的關鍵因素,但信號在端點處往往為非極值點,更不可能同時為極大值點和極小值點,以致三次樣條曲線容易在端點處形成較大的擺動,端點附近插值求包絡過程中便會出現誤差,并且隨著分解過程的進行,重復進行插值求包絡線,導致誤差不斷累積并向內傳播污染數據,從而影響分解結果,這就是EMD方法的端點效應問題。如果數據序列較短,端點效應將嚴重影響分解結果,尤其對于低頻的IMF分量,其時間尺度大,極值間隔距離大,同時極值點相對較少,端點效應現象更為突出[6],誤差的積累與傳播導致分解出的低頻IMF完全扭曲,喪失物理意義。因此,為提高EMD的分解質量,對其進行端點效應抑制是十分必要的。
目前,解決EMD端點效應的方法主要有:① 改進插值方法,盡可能降低端點處的擬合誤差;② 對數據進行延拓,形成新的極值點,將擬合誤差置于原始數據段之外。其中第二種方法對EMD端點效應抑制效果明顯,成為當前研究的主要方向。
通過Hilbert變換計算解析信號的瞬時頻率,得到信號Hilbert譜和邊際譜,理論上是理想的,但在實現上是采用快速傅里葉變換實現的。Hilbert變換的端點效應產生于求取90°共軛信號的過程,共軛信號的求取是通過“傅里葉變換-雙邊譜對折為單邊譜-傅里葉逆變換”獲得的,在數字方法實現時由于數據的非完整周期采樣,進行傅里葉變換便會出現頻譜泄露現象[7],時頻圖的兩端也會出現嚴重的突變和振蕩,即Hilbert變換的端點效應,影響對數據的分析。
克服Hilbert變換的端點效應問題最理想的方法就是對數據進行整周期采樣。但在工程實際中即使是旋轉機械產生的信號由于旋轉體頻率脈動及信號采集過程中不可避免的隨機干擾、噪聲等問題的存在,信號不存在嚴格的周期性,所以數據的整周期采樣是不現實的。目前,最常用的抑制Hilbert變換端點效應的方法就是對待變換的數據進行延拓,將端點效應影響盡可能的置于有效數據之外。
本文綜合數據預測延拓和極值延拓兩種延拓方法,提出基于時間尺度的LS-SVM端點延拓方法,利用最小二乘支持向量機(Least Squares Support Vector Machine,LS-SVM)良好的回歸預測性能對信號極值點幅值進行預測延拓,同時根據信號極值時間尺度設計了一種自適用時間尺度延拓算法,然后結合埃爾米特插值方法對數據進行延拓。該算法既參考了端點附近極值點的變化趨勢對延拓極值的影響,又考慮了內部極值信息對延拓極值的影響,充分利用了已知極值點信息,具有很強的自適應性。利用這種方法對原始信號只需進行一次延拓后即可同時解決兩種端點效應的影響,簡化了數據的計算。算法具體描述如下:
(1)極值點幅值的延拓:波動信號是有許多不同時間尺度的基本波動模式疊加而成,因此基于極值點信息的時間尺度并不是隨機的、無規律的,特別對于旋轉機械而言,其振動數據都有一定的規律性[9],因此可以利用LS-SVM對極值點幅值進行預測延拓。
(2)極值點時間的延拓:對實際振動信號而言,在某一時間點一定范圍內振動形式是相對穩定的,而在整個振動信號序列中卻是非穩定的,因此需要一種自適應的估計算法來預測極值點的發生時刻。
設信號x(t)所有極大值點時刻為{tmax(i),i=1,…,n},幅值為{Amax(i),i=1,…,n},相鄰極大值的時間間隔為{Δtmax(i),i=1,…,n-1},n為極大值點數;同理所有極小值點時刻為{tmin(i),i=1,…,m},幅值為{Amin(i),i=1,…,n},相鄰極小值的時間間隔為{Δtmin(i),i=1,…,m-1},m為極小值點數。
以數據右端點極大值點的一步延拓為例,定義端點處時刻為tend,幅值為Aend,描述端點延拓算法為:

p為第n個極大值點左側參考極大值點的個數,k為加權系數,且k∈(0,2)。
① 當 Δtmax(i+1)>tend-tmax(i)時,

② 當 Δtmax(i+1)<tend-tmax(i)時,

① 當 Δtmax(i+1)>tend-tmax(i)時,

② 當 Δtmax(i+1)<tend-tmax(i)時,

經過大量實驗發現,根據信號變化劇烈程度P一般取3~6個預測效果較好(如極值點數小于P,則P可取1,但在實際工程中振動信號極值點較為密集,一般不存在此問題)。
(3)插值延拓:對延拓后的極值點進行插值,從而達到對原始信號進行延拓的目的。由于延拓預測的信息是極值點信息,所以進行插值時必須保證插值函數在極值點處取到相應的極值,這就成為一個帶導數的差值問題,本文采用埃爾米特插值方法,可保證插值函數的極值與預測的極值信息相對應,同時保證了延拓信號的平滑性。
工程實際中測得的振動信號,是各設備振動和各振動模式的疊加,一般都是非平穩過程,特別是出現故障或負荷、轉速發生變化時,非平穩現象加劇。采用某典型振動信號的時域波形如圖1所示,以該信號右端延拓為例,驗證基于極值時間尺度的LS-SVM端點延拓方法的有效性。

圖1 原始振動信號圖Fig.1 Wave of original vibration signal
對該信號EMD分解時得到的包絡線未包含所有的數據,出現失真現象,產生的端點效應會污染整個數據序列。采用最小二乘支持向量機對數據進行延拓,就是在數據序列的兩端增加兩個極大值和極小值,對延拓的信號段進行埃爾米特插值,即可達到對原始信號進行延拓的目的。
首先求取信號所有極值點幅值信息,得到信號的極大值和極小值。用信號的第37~58的極大值作為最小二乘支持向量機學習訓練樣本,一步預測得第59個極大值點,兩步預測得第60個數據點,并與真實的極大值信息作比較,結果如圖2所示。可見,基于LSSVM的方法對該信號極值點幅值的預測精度較高,準確反映了極值點運動趨勢。

圖2 極值點幅值預測Fig.2 Extremum point amplitude forecast
然后以極值點所處位置的數據點為計算單位,對該信號極值點的時刻進行預測計算,根據式(2)~式(6)的計算方法和判斷條件,可得極值點時刻預測結果和極值點幅值預測結果。最后對延拓后的信號埃爾米特插值,即可得到延拓信號。圖3為原始信號EMD分解結果,圖4為基于極值時間尺度的LS-SVM算法實現的EMD分解結果。在圖3中,data1表示原始信號,data2表示EMD分解結果,可以看出EMD分解后有許多失真,尤其對波段較少的低頻分量IMF2來說,端點效應現象較明顯,對比圖4可以看出,經基于極值時間尺度的LS-SVM算法延拓后的信號進行EMD分解后,無論是IMF1還是IMF2都很好地解決了失真與端點效應問題。

圖3 原始信號EMD分解Fig.3 Original signal by EMD

圖4 基于LS-SVM算法的EMD分解Fig.4 Signal EMD based on LS-SVM
圖5為信號延拓處理后IMF瞬時頻率。可以看出,經本方法延拓后的信號EMD分解,IMF1、IMF2較好的反映了信號兩個基本波動模式,對各IMF進行Hilbert變換后求取瞬時頻率,時頻圖的端點未出現大的波動,IMF1、IMF2的瞬時頻率與信號兩個基頻頻率嚴格對應。

圖5 信號延拓處理后各IMF瞬時頻率Fig.5 IMF instantaneous frequency of signal end expansion
為了對比該方法的延拓效果,取振動信號的一段數據向后延拓100個數據點,效果對比如圖6所示,其中data1為原始信號,data2為基于極值時間尺度的LS-SVM端點延拓后信號,data3為利用神經網絡方法延拓后信號。從圖6可以看出,本文算法明顯優于神經網絡算法,實際運算中,對原始信號進行LS-SVM延拓時,預測到第1 009個數據點處出現第一個延拓極值,為極小值,幅值為-1.270 2,而實際幅值為-1.299 7;再經過預測出現第二個極值,為極大值,幅值為0.820 5,而實際幅值為 0.833 9,可見精度和準確度都比較高。利用神經網絡對信號進行延拓,極值存在較大的誤差,經過45步預測后的數據已經完全喪失了原信號的特征。表1為兩種不同延拓方法的性能比較,可以看出,無論運算時間還是延拓信號與真實信號誤差值,本文算法明顯優于神經網絡算法。

圖6 信號端點延拓效果比較Fig.6 The comparison of signal end expansion effect

表1 不同方法性能比較Tab.1 The comparison of different methods effect
在發動機狀態監控或故障診斷中,轉子運轉狀況最能反映發動機狀況,對轉子振動信號的采集和分析可以有效監控發動機狀況。在試車臺上采集某型渦噴發動機的振動數據,采樣頻率為4 000 Hz,采樣時間為1 s,其時域波形如圖7所示。

圖7 航空發動機振動信號Fig.7 Wave of aeroengine original vibration signal
對發動機振動信號分析,就是從數據序列中找出或分離出能夠反映轉子振動特性的信號。由于航空發動機內部振動模式復雜,信號非平穩性強,為了提取有用的高頻信號,需要對信號EMD分解,在EMD分解時生成的 IMF較多,獲得每個 IMF要經過多次,導致EMD端點效應向內傳播嚴重。所以在用EMD對發動機振動信號進行處理時必須抑制端點效應,否則分解結果會毫無意義。

圖8 延拓后振動信號EMD分解Fig.8 Vibration signal after extended decomposed by EMD
圖8是對原始信號延拓后進行EMD分解結果(這里只列出前五個IMF分量)。可以看出:① 適當增多信號延拓長度,有利于獲得更好的分解效果;② 信號延拓后對兩種端點效應都起到了良好的抑制作用;③經EMD分解得到的IMF2分量具有明顯的調幅特征。
對各IMF進行Hilbert變換,得到Hilbert時頻譜,去除端點延拓部分得到各IMF的瞬時頻率,如圖9所示。

圖9 各IMF的瞬時頻率Fig.9 Instantaneous frequency of IMF
可以看出:① IMF1處在多倍轉頻的高頻段,能量占主導地位,為機體設備產生的高頻背景噪聲;② 各IMF經Hilbert變換后求取瞬時頻率,除低頻分量IMF7外,各瞬時頻率端點處未出現大的波動;③ 計算各IMF的歸一化互信息可判定IMF7、IMF8、IMF9為虛假低頻分量;④ IMF2與2倍轉頻較為接近,IMF3與1倍轉頻較為接近,并且它們的幅值分布也表現出一定的規律性;⑤ IMF2、IMF3、IMF4、IMF5處于發動機轉子工作頻率段,頻率范圍比較集中,頻率成分較為單一,與轉子的實際運轉頻率有著較好的對應關系,反映了轉子基頻振動特性,因此可選定這幾個分量進行轉子信號特征提取,作為狀態監控與故障診斷的參考。
因此,經基于極值時間尺度的LS-SVM端點延拓技術改進后的EMD時頻分析方法,能夠較好的降低EMD分解時產生的模式混疊現象,減少虛假模式產生,有效抑制EMD分解和Hilbert變換的端點效應問題,提高航空發動機振動信號特征提取的精度與準確度,為發動機轉子狀態監控和故障診斷提供了有力的技術支撐。
本文針對HHT/EMD算法中信號兩端發散所引起的端點效應問題,提出一種基于極值時間尺度的LSSVM端點延拓方法,該方法既考慮信號端點處的運動特點,又兼顧信號整體的發展趨勢,仿真結果表明該方法適應性強,有效延拓距離長,延拓速度快。將該方法應用于航空發動機振動信號分析上,在抑制端點效應的同時,較好的分解出反映轉子振動特性的信號,為轉子狀態監控與故障診斷提供數據支持,具有較高的實用價值。
[1] 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,1998,454:903-995.
[2]胡愛軍.Hilbert-Huang變換在旋轉機械振動信號分析中的應用研究[D].保定:華北電力大學,2008.
[3]黃大吉,趙進平,蘇紀蘭.希爾伯特-黃變換的端點延拓[J].海洋學報,2003,25(1):1-11.
[4] Zhao J P,Huang D J.Mirror extending and circular spline function for empirical mode decomposition method[J].Journal of Zhejiang University,2001,2(3):247-252.
[5]鄧擁軍.EMD方法及Hilbert變換中邊界問題的處理[J].科學通報,2001,46(3):257-263.
[6]張郁山.應用自回歸模型處理EMD方法中的邊界問題[J].自然科學通報,2003,13(10):1054-1059.
[7]楊建文,賈民平.希爾伯特-黃譜的端點效應分析及處理方法研究[J].振動工程學報,2006,19(2):283-288.
[8]程軍圣,于德介,楊 宇.Hilbert-Huang變換端點效應問題的探討[J].振動與沖擊,2005,24(6):40-42.
[9]程軍圣,于德介,楊 宇.基于支持矢量回歸機的Hilbert-Huang變換端點效應問題的處理方法[J].機械工程學報,2006,42(4):23-31.
[10]曹沖鋒,楊世錫,楊將新.一種抑制EMD端點效應新方法及其在信號特征提取中的應用[J].振動工程學報,2008,21(6):588-592.
[11]沈 路,周曉軍,張志剛,等.Hilbert-Huang變換中的一種端點延拓方法[J].振動與沖擊,2009,28(8):168-171.
[12]楊永鋒,吳亞鋒,任興民,等.基于最大Lyapunov指數預測的EMD 端點延拓[J].物理學報,2009,58(6):3742-3745.
[13]徐曉剛,徐冠雷,王孝通,等.經驗模式分解(EMD)及其應用[J].電子學報,2009,37(3):581-585.
[14] Qi K,He Z G,Zi Y Y.Cosine window-based boundary processing method for EMD and its application in rubbing faultdiagnosis [J]. MechanicalSystems and Signal Processing,2007,21(7):2750-2760.