振動信號中HHT/EMD端點延拓方法研究

郭明威，倪世宏，朱家海，張志鵬

(空軍工程大學(xué) 工程學(xué)院，西安 710038)

Hilbert-Huang變換(HHT)是 Huang等［1］提出的一種適應(yīng)于非線性、非平穩(wěn)信號的時頻處理方法。HHT的提出被認(rèn)為是近年來以傅里葉變換為基礎(chǔ)的線性和穩(wěn)態(tài)譜分析的一個重大突破，經(jīng)過十幾年的發(fā)展，HHT理論體系不斷完善，已成功應(yīng)用于地震信號、圖像處理、故障診斷等學(xué)科領(lǐng)域［2-14］。但在HHT應(yīng)用中，端點效應(yīng)問題直接影響到信號HHT分析結(jié)果，已成為阻礙其工程實用的難點問題［1-12］。

針對端點效應(yīng)問題，國內(nèi)外學(xué)者做了大量研究，目前提出的抑制端點效應(yīng)的方法大致分為三類:波形延拓法［2-4］、數(shù)據(jù)預(yù)測延拓法［5-9］、極值延拓法［3，10-11］。在波形延拓方面，波形匹配延拓法［2］利用信號內(nèi)部變化趨勢來延拓信號，結(jié)果更接近原始數(shù)據(jù)特征，該方法延拓后會出現(xiàn)數(shù)據(jù)“斷裂”現(xiàn)象，引入“異常事件”而影響分解效果;波形鏡像延拓法［3-4］在具有對稱性的極值所在位置放置鏡子，實現(xiàn)鏡內(nèi)數(shù)據(jù)的閉合延拓，該方法對信號的對稱性要求較高，且沒有從數(shù)據(jù)上給出最佳的鏡面位置。在數(shù)據(jù)預(yù)測延拓方面，利用神經(jīng)網(wǎng)絡(luò)建模對數(shù)據(jù)進(jìn)行預(yù)測延拓［5］，借以抑制端點效應(yīng)，但該方法對非線性、非平穩(wěn)信號延拓效果不佳;文獻(xiàn)［6］采用自回歸(AR)模型處理EMD方法中的邊界問題，取得了一定的效果，但AR模型只適應(yīng)于平穩(wěn)序列;文獻(xiàn)［7］采用時間序列建模與預(yù)測方法對信號數(shù)據(jù)進(jìn)行延拓，對端點效應(yīng)進(jìn)行了探討，但對于非平穩(wěn)性較強(qiáng)的振動信號預(yù)測效果不佳;基于支持矢量回歸機(jī)的端點效應(yīng)問題處理方法［9］，比神經(jīng)網(wǎng)絡(luò)和AR模型預(yù)測效果好，但該法用于學(xué)習(xí)訓(xùn)練的數(shù)據(jù)點多，預(yù)測步數(shù)過多，致使計算量龐大，并且延拓誤差隨多步預(yù)測而積累變大，預(yù)測精度下降。在極值延拓方面，包絡(luò)極值延拓方法［3］實質(zhì)就是對極值點的等幅、等距延拓，不能反應(yīng)信號趨勢，并且此方法未考慮到特殊端點的處理;端點優(yōu)化對稱延拓方法［10］比ARMA延拓方法效果要好，但端點優(yōu)化對稱延拓方法對端點附近數(shù)據(jù)處理不夠平滑，可能引起數(shù)據(jù)的階躍;基于相似極值延拓［11］，其時間的延拓為等距延拓，延拓極值點的大小取值為已知極值點序列的均值，導(dǎo)致其自適應(yīng)性不強(qiáng)，延拓效果不理想。

本文提出一種基于時間尺度的LS-SVM端點延拓方法，通過仿真信號分析，并將其應(yīng)用到航空發(fā)動機(jī)的振動數(shù)據(jù)處理中，驗證了該方法能有效地抑制端點效應(yīng)。

1 端點效應(yīng)分析

經(jīng)驗?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)和與之相應(yīng)的Hilbert譜信號分析方法統(tǒng)稱為Hilbert-Huang變換。HHT的端點效應(yīng)現(xiàn)象在EMD和Hilbert變換中都存在，這兩種失真的疊加，造成了HHT的端點部分無法正確反映信號所包含的信息。

1.1 EMD分解端點效應(yīng)分析

EMD本質(zhì)是通過數(shù)據(jù)的特征時間尺度來獲得固有波動模式，其基本思想是:對一給定信號，先獲得信號極值點，通過插值獲得信號包絡(luò)，得到均值，與均值的差得到分解的一層信號，如此重復(fù)，直到將信號分解成有限個基本模式分量(IMF)和殘余項rn(t)的組合，表示為:

EMD是基于數(shù)據(jù)極點信息的分解過程，或者說極點是信號的EMD分解的關(guān)鍵因素，但信號在端點處往往為非極值點，更不可能同時為極大值點和極小值點，以致三次樣條曲線容易在端點處形成較大的擺動，端點附近插值求包絡(luò)過程中便會出現(xiàn)誤差，并且隨著分解過程的進(jìn)行，重復(fù)進(jìn)行插值求包絡(luò)線，導(dǎo)致誤差不斷累積并向內(nèi)傳播污染數(shù)據(jù)，從而影響分解結(jié)果，這就是EMD方法的端點效應(yīng)問題。如果數(shù)據(jù)序列較短，端點效應(yīng)將嚴(yán)重影響分解結(jié)果，尤其對于低頻的IMF分量，其時間尺度大，極值間隔距離大，同時極值點相對較少，端點效應(yīng)現(xiàn)象更為突出［6］，誤差的積累與傳播導(dǎo)致分解出的低頻IMF完全扭曲，喪失物理意義。因此，為提高EMD的分解質(zhì)量，對其進(jìn)行端點效應(yīng)抑制是十分必要的。

目前，解決EMD端點效應(yīng)的方法主要有:① 改進(jìn)插值方法，盡可能降低端點處的擬合誤差;② 對數(shù)據(jù)進(jìn)行延拓，形成新的極值點，將擬合誤差置于原始數(shù)據(jù)段之外。其中第二種方法對EMD端點效應(yīng)抑制效果明顯，成為當(dāng)前研究的主要方向。

1.2 Hilbert變換端點效應(yīng)分析

通過Hilbert變換計算解析信號的瞬時頻率，得到信號Hilbert譜和邊際譜，理論上是理想的，但在實現(xiàn)上是采用快速傅里葉變換實現(xiàn)的。Hilbert變換的端點效應(yīng)產(chǎn)生于求取90°共軛信號的過程，共軛信號的求取是通過“傅里葉變換-雙邊譜對折為單邊譜-傅里葉逆變換”獲得的，在數(shù)字方法實現(xiàn)時由于數(shù)據(jù)的非完整周期采樣，進(jìn)行傅里葉變換便會出現(xiàn)頻譜泄露現(xiàn)象［7］，時頻圖的兩端也會出現(xiàn)嚴(yán)重的突變和振蕩，即Hilbert變換的端點效應(yīng)，影響對數(shù)據(jù)的分析。

克服Hilbert變換的端點效應(yīng)問題最理想的方法就是對數(shù)據(jù)進(jìn)行整周期采樣。但在工程實際中即使是旋轉(zhuǎn)機(jī)械產(chǎn)生的信號由于旋轉(zhuǎn)體頻率脈動及信號采集過程中不可避免的隨機(jī)干擾、噪聲等問題的存在，信號不存在嚴(yán)格的周期性，所以數(shù)據(jù)的整周期采樣是不現(xiàn)實的。目前，最常用的抑制Hilbert變換端點效應(yīng)的方法就是對待變換的數(shù)據(jù)進(jìn)行延拓，將端點效應(yīng)影響盡可能的置于有效數(shù)據(jù)之外。

2 基于極值時間尺度的LS-SVM端點延拓

本文綜合數(shù)據(jù)預(yù)測延拓和極值延拓兩種延拓方法，提出基于時間尺度的LS-SVM端點延拓方法，利用最小二乘支持向量機(jī)(Least Squares Support Vector Machine，LS-SVM)良好的回歸預(yù)測性能對信號極值點幅值進(jìn)行預(yù)測延拓，同時根據(jù)信號極值時間尺度設(shè)計了一種自適用時間尺度延拓算法，然后結(jié)合埃爾米特插值方法對數(shù)據(jù)進(jìn)行延拓。該算法既參考了端點附近極值點的變化趨勢對延拓極值的影響，又考慮了內(nèi)部極值信息對延拓極值的影響，充分利用了已知極值點信息，具有很強(qiáng)的自適應(yīng)性。利用這種方法對原始信號只需進(jìn)行一次延拓后即可同時解決兩種端點效應(yīng)的影響，簡化了數(shù)據(jù)的計算。算法具體描述如下:

(1)極值點幅值的延拓:波動信號是有許多不同時間尺度的基本波動模式疊加而成，因此基于極值點信息的時間尺度并不是隨機(jī)的、無規(guī)律的，特別對于旋轉(zhuǎn)機(jī)械而言，其振動數(shù)據(jù)都有一定的規(guī)律性［9］，因此可以利用LS-SVM對極值點幅值進(jìn)行預(yù)測延拓。

(2)極值點時間的延拓:對實際振動信號而言，在某一時間點一定范圍內(nèi)振動形式是相對穩(wěn)定的，而在整個振動信號序列中卻是非穩(wěn)定的，因此需要一種自適應(yīng)的估計算法來預(yù)測極值點的發(fā)生時刻。

設(shè)信號x(t)所有極大值點時刻為{tmax(i)，i=1，…，n}，幅值為{Amax(i)，i=1，…，n}，相鄰極大值的時間間隔為{Δtmax(i)，i=1，…，n-1}，n為極大值點數(shù);同理所有極小值點時刻為{tmin(i)，i=1，…，m}，幅值為{Amin(i)，i=1，…，n}，相鄰極小值的時間間隔為{Δtmin(i)，i=1，…，m-1}，m為極小值點數(shù)。

以數(shù)據(jù)右端點極大值點的一步延拓為例，定義端點處時刻為tend，幅值為Aend，描述端點延拓算法為:

p為第n個極大值點左側(cè)參考極大值點的個數(shù)，k為加權(quán)系數(shù)，且k∈(0，2)。

① 當(dāng) Δtmax(i+1)＞tend-tmax(i)時，

② 當(dāng) Δtmax(i+1)＜tend-tmax(i)時，

① 當(dāng) Δtmax(i+1)＞tend-tmax(i)時，

② 當(dāng) Δtmax(i+1)＜tend-tmax(i)時，

經(jīng)過大量實驗發(fā)現(xiàn)，根據(jù)信號變化劇烈程度P一般取3～6個預(yù)測效果較好(如極值點數(shù)小于P，則P可取1，但在實際工程中振動信號極值點較為密集，一般不存在此問題)。

(3)插值延拓:對延拓后的極值點進(jìn)行插值，從而達(dá)到對原始信號進(jìn)行延拓的目的。由于延拓預(yù)測的信息是極值點信息，所以進(jìn)行插值時必須保證插值函數(shù)在極值點處取到相應(yīng)的極值，這就成為一個帶導(dǎo)數(shù)的差值問題，本文采用埃爾米特插值方法，可保證插值函數(shù)的極值與預(yù)測的極值信息相對應(yīng)，同時保證了延拓信號的平滑性。

3 仿真實例分析

工程實際中測得的振動信號，是各設(shè)備振動和各振動模式的疊加，一般都是非平穩(wěn)過程，特別是出現(xiàn)故障或負(fù)荷、轉(zhuǎn)速發(fā)生變化時，非平穩(wěn)現(xiàn)象加劇。采用某典型振動信號的時域波形如圖1所示，以該信號右端延拓為例，驗證基于極值時間尺度的LS-SVM端點延拓方法的有效性。

圖1 原始振動信號圖Fig.1 Wave of original vibration signal

對該信號EMD分解時得到的包絡(luò)線未包含所有的數(shù)據(jù)，出現(xiàn)失真現(xiàn)象，產(chǎn)生的端點效應(yīng)會污染整個數(shù)據(jù)序列。采用最小二乘支持向量機(jī)對數(shù)據(jù)進(jìn)行延拓，就是在數(shù)據(jù)序列的兩端增加兩個極大值和極小值，對延拓的信號段進(jìn)行埃爾米特插值，即可達(dá)到對原始信號進(jìn)行延拓的目的。

首先求取信號所有極值點幅值信息，得到信號的極大值和極小值。用信號的第37～58的極大值作為最小二乘支持向量機(jī)學(xué)習(xí)訓(xùn)練樣本，一步預(yù)測得第59個極大值點，兩步預(yù)測得第60個數(shù)據(jù)點，并與真實的極大值信息作比較，結(jié)果如圖2所示。可見，基于LSSVM的方法對該信號極值點幅值的預(yù)測精度較高，準(zhǔn)確反映了極值點運(yùn)動趨勢。

圖2 極值點幅值預(yù)測Fig.2 Extremum point amplitude forecast

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

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

圖4 基于LS-SVM算法的EMD分解Fig.4 Signal EMD based on LS-SVM

圖5為信號延拓處理后IMF瞬時頻率。可以看出，經(jīng)本方法延拓后的信號EMD分解，IMF1、IMF2較好的反映了信號兩個基本波動模式，對各IMF進(jìn)行Hilbert變換后求取瞬時頻率，時頻圖的端點未出現(xiàn)大的波動，IMF1、IMF2的瞬時頻率與信號兩個基頻頻率嚴(yán)格對應(yīng)。

圖5 信號延拓處理后各IMF瞬時頻率Fig.5 IMF instantaneous frequency of signal end expansion

為了對比該方法的延拓效果，取振動信號的一段數(shù)據(jù)向后延拓100個數(shù)據(jù)點，效果對比如圖6所示，其中data1為原始信號，data2為基于極值時間尺度的LS-SVM端點延拓后信號，data3為利用神經(jīng)網(wǎng)絡(luò)方法延拓后信號。從圖6可以看出，本文算法明顯優(yōu)于神經(jīng)網(wǎng)絡(luò)算法，實際運(yùn)算中，對原始信號進(jìn)行LS-SVM延拓時，預(yù)測到第1 009個數(shù)據(jù)點處出現(xiàn)第一個延拓極值，為極小值，幅值為-1.270 2，而實際幅值為-1.299 7;再經(jīng)過預(yù)測出現(xiàn)第二個極值，為極大值，幅值為0.820 5，而實際幅值為 0.833 9，可見精度和準(zhǔn)確度都比較高。利用神經(jīng)網(wǎng)絡(luò)對信號進(jìn)行延拓，極值存在較大的誤差，經(jīng)過45步預(yù)測后的數(shù)據(jù)已經(jīng)完全喪失了原信號的特征。表1為兩種不同延拓方法的性能比較，可以看出，無論運(yùn)算時間還是延拓信號與真實信號誤差值，本文算法明顯優(yōu)于神經(jīng)網(wǎng)絡(luò)算法。

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

表1 不同方法性能比較Tab.1 The comparison of different methods effect

4 在航空發(fā)動機(jī)振動信號中的應(yīng)用

在發(fā)動機(jī)狀態(tài)監(jiān)控或故障診斷中，轉(zhuǎn)子運(yùn)轉(zhuǎn)狀況最能反映發(fā)動機(jī)狀況，對轉(zhuǎn)子振動信號的采集和分析可以有效監(jiān)控發(fā)動機(jī)狀況。在試車臺上采集某型渦噴發(fā)動機(jī)的振動數(shù)據(jù)，采樣頻率為4 000 Hz，采樣時間為1 s，其時域波形如圖7所示。

圖7 航空發(fā)動機(jī)振動信號Fig.7 Wave of aeroengine original vibration signal

對發(fā)動機(jī)振動信號分析，就是從數(shù)據(jù)序列中找出或分離出能夠反映轉(zhuǎn)子振動特性的信號。由于航空發(fā)動機(jī)內(nèi)部振動模式復(fù)雜，信號非平穩(wěn)性強(qiáng)，為了提取有用的高頻信號，需要對信號EMD分解，在EMD分解時生成的 IMF較多，獲得每個 IMF要經(jīng)過多次，導(dǎo)致EMD端點效應(yīng)向內(nèi)傳播嚴(yán)重。所以在用EMD對發(fā)動機(jī)振動信號進(jìn)行處理時必須抑制端點效應(yīng)，否則分解結(jié)果會毫無意義。

圖8 延拓后振動信號EMD分解Fig.8 Vibration signal after extended decomposed by EMD

圖8是對原始信號延拓后進(jìn)行EMD分解結(jié)果(這里只列出前五個IMF分量)。可以看出:① 適當(dāng)增多信號延拓長度，有利于獲得更好的分解效果;② 信號延拓后對兩種端點效應(yīng)都起到了良好的抑制作用;③經(jīng)EMD分解得到的IMF2分量具有明顯的調(diào)幅特征。

對各IMF進(jìn)行Hilbert變換，得到Hilbert時頻譜，去除端點延拓部分得到各IMF的瞬時頻率，如圖9所示。

圖9 各IMF的瞬時頻率Fig.9 Instantaneous frequency of IMF

可以看出:① IMF1處在多倍轉(zhuǎn)頻的高頻段，能量占主導(dǎo)地位，為機(jī)體設(shè)備產(chǎn)生的高頻背景噪聲;② 各IMF經(jīng)Hilbert變換后求取瞬時頻率，除低頻分量IMF7外，各瞬時頻率端點處未出現(xiàn)大的波動;③ 計算各IMF的歸一化互信息可判定IMF7、IMF8、IMF9為虛假低頻分量;④ IMF2與2倍轉(zhuǎn)頻較為接近，IMF3與1倍轉(zhuǎn)頻較為接近，并且它們的幅值分布也表現(xiàn)出一定的規(guī)律性;⑤ IMF2、IMF3、IMF4、IMF5處于發(fā)動機(jī)轉(zhuǎn)子工作頻率段，頻率范圍比較集中，頻率成分較為單一，與轉(zhuǎn)子的實際運(yùn)轉(zhuǎn)頻率有著較好的對應(yīng)關(guān)系，反映了轉(zhuǎn)子基頻振動特性，因此可選定這幾個分量進(jìn)行轉(zhuǎn)子信號特征提取，作為狀態(tài)監(jiān)控與故障診斷的參考。

因此，經(jīng)基于極值時間尺度的LS-SVM端點延拓技術(shù)改進(jìn)后的EMD時頻分析方法，能夠較好的降低EMD分解時產(chǎn)生的模式混疊現(xiàn)象，減少虛假模式產(chǎn)生，有效抑制EMD分解和Hilbert變換的端點效應(yīng)問題，提高航空發(fā)動機(jī)振動信號特征提取的精度與準(zhǔn)確度，為發(fā)動機(jī)轉(zhuǎn)子狀態(tài)監(jiān)控和故障診斷提供了有力的技術(shù)支撐。

5 結(jié)論

本文針對HHT/EMD算法中信號兩端發(fā)散所引起的端點效應(yīng)問題，提出一種基于極值時間尺度的LSSVM端點延拓方法，該方法既考慮信號端點處的運(yùn)動特點，又兼顧信號整體的發(fā)展趨勢，仿真結(jié)果表明該方法適應(yīng)性強(qiáng)，有效延拓距離長，延拓速度快。將該方法應(yīng)用于航空發(fā)動機(jī)振動信號分析上，在抑制端點效應(yīng)的同時，較好的分解出反映轉(zhuǎn)子振動特性的信號，為轉(zhuǎn)子狀態(tài)監(jiān)控與故障診斷提供數(shù)據(jù)支持，具有較高的實用價值。

［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變換在旋轉(zhuǎn)機(jī)械振動信號分析中的應(yīng)用研究［D］.保定:華北電力大學(xué)，2008.

［3］黃大吉，趙進(jìn)平，蘇紀(jì)蘭.希爾伯特-黃變換的端點延拓［J］.海洋學(xué)報，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］.科學(xué)通報，2001，46(3):257-263.

［6］張郁山.應(yīng)用自回歸模型處理EMD方法中的邊界問題［J］.自然科學(xué)通報，2003，13(10):1054-1059.

［7］楊建文，賈民平.希爾伯特-黃譜的端點效應(yīng)分析及處理方法研究［J］.振動工程學(xué)報，2006，19(2):283-288.

［8］程軍圣，于德介，楊 宇.Hilbert-Huang變換端點效應(yīng)問題的探討［J］.振動與沖擊，2005，24(6):40-42.

［9］程軍圣，于德介，楊 宇.基于支持矢量回歸機(jī)的Hilbert-Huang變換端點效應(yīng)問題的處理方法［J］.機(jī)械工程學(xué)報，2006，42(4):23-31.

［10］曹沖鋒，楊世錫，楊將新.一種抑制EMD端點效應(yīng)新方法及其在信號特征提取中的應(yīng)用［J］.振動工程學(xué)報，2008，21(6):588-592.

［11］沈 路，周曉軍，張志剛，等.Hilbert-Huang變換中的一種端點延拓方法［J］.振動與沖擊，2009，28(8):168-171.

［12］楊永鋒，吳亞鋒，任興民，等.基于最大Lyapunov指數(shù)預(yù)測的EMD 端點延拓［J］.物理學(xué)報，2009，58(6):3742-3745.

［13］徐曉剛，徐冠雷，王孝通，等.經(jīng)驗?zāi)Ｊ椒纸?EMD)及其應(yīng)用［J］.電子學(xué)報，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.