基于HHT的航空發(fā)動機(jī)氣動失穩(wěn)信號時頻分析方法

郭貴喜，王振華，王洪祥，高輝國，薛秀生

（中航工業(yè)沈陽發(fā)動機(jī)設(shè)計研究所，沈陽110015）

基于HHT的航空發(fā)動機(jī)氣動失穩(wěn)信號時頻分析方法

郭貴喜，王振華，王洪祥，高輝國，薛秀生

（中航工業(yè)沈陽發(fā)動機(jī)設(shè)計研究所，沈陽110015）

針對航空發(fā)動機(jī)工作范圍內(nèi)存在的氣動失穩(wěn)現(xiàn)象，運用希爾伯特-黃變換（HHT）分析其信號的時變特征；通過對其處理非線性、非平穩(wěn)信號分析新方法與其他時頻分析方法的對比，并對HHT存在的問題進(jìn)行針對性解決，使其適應(yīng)航空發(fā)動機(jī)氣動失穩(wěn)數(shù)據(jù)分析的要求。并在Labview平臺上通過數(shù)字仿真試驗實現(xiàn)和驗證了HHT方法；結(jié)果表明：此算法準(zhǔn)確有效；通過航空發(fā)動機(jī)工程試驗數(shù)據(jù)的處理過程，驗證了HHT在處理此類相關(guān)問題時的可行性、適用性，同時指出其仍然存在的缺陷。

時頻分析；希爾伯特-黃變換；航空發(fā)動機(jī)；氣動失穩(wěn)

0 引言

航空發(fā)動機(jī)在其工作范圍內(nèi)存在“旋轉(zhuǎn)失速”和“喘振”2類氣動失穩(wěn)現(xiàn)象，危害性極大。目前，發(fā)動機(jī)失穩(wěn)現(xiàn)象是世界各國航空科研人員研究的重要課題。發(fā)動機(jī)氣動失穩(wěn)產(chǎn)生的流路內(nèi)壓力脈動本質(zhì)上是一種非線性，非平穩(wěn)信號，特別是在過渡過程中時變性很強(qiáng)。采用普通的處理線性、平穩(wěn)信號的分析方法有非常大的局限性。目前大多數(shù)時頻變換方法存在著不足，如所使用的基函數(shù)固定不變，難以準(zhǔn)確匹配信號中的時頻結(jié)構(gòu)，而且時頻分辨率還受到Heisenberg不確定性原理的限制[1-2]等。

本文應(yīng)用希爾伯特-黃變換（HHT）對某型發(fā)動機(jī)壓力脈動信號進(jìn)行時頻分析，提取信號的幅值和頻率的時變特征，并將分析結(jié)果和傅里葉變換分析方法進(jìn)行對比。

1 傳統(tǒng)的時頻分析方法

自然界中的信號幾乎都是多種信號成分的疊加，既有平穩(wěn)信號也有非平穩(wěn)信號；信號既可能是線性的也可能是非線性的。傳統(tǒng)的基于傅里葉變換的信號處理方法難以避免其局限性。它將信號從時域整個變換到頻域，其譜分析僅能反映一段整體時間內(nèi)所有頻率與其幅值的對應(yīng)關(guān)系，而對頻率隨時間的變化規(guī)律無法分析，也就是說它不能刻畫出信號的頻率等尺度在局部的瞬時特征。針對傅里葉變換不能反映信號瞬時特征的缺點，研究出多種時頻分析方法。

1.1短時傅里葉變換

短時傅里葉變換（STFT）是1種最基本的時頻分析方法。它通過1個滑動時窗，分段對信號進(jìn)行傅里葉變換，實際上就是先對信號進(jìn)行加窗處理，然后再進(jìn)行傅里葉變換。這種方法在一定程度克服了傅里葉變換的缺點，實現(xiàn)了一定的分析信號局部特征的能力。但是窗口一旦選定，就不能任意更改，如果信號在時間或者頻率上的變化小于窗口，窗口內(nèi)信號平穩(wěn)的假設(shè)就不能成立，信號的局部特征就不能被反映。而且它不能在時間和頻率2個方面同時達(dá)到很高的分辨率，因此，短時傅里葉變換只適用于緩變信號的分析處理。

1.2小波變換

小波變換（WT）是20世紀(jì)80年代中后期發(fā)展起來的1種新型線性時頻分析方法。該方法實現(xiàn)了窗口可調(diào)，在高頻時使用窄窗口，在低頻時使用寬窗口。正是基于這種多分辨率的思想，小波變換被譽(yù)為數(shù)學(xué)顯微鏡。但其本質(zhì)上仍是是1種窗口可調(diào)的傅里葉變換，不可避免地仍具有窗函數(shù)的局限性。而且小波變換存在眾多的的小波基函數(shù)，各種基函數(shù)使用的范圍不同，造成小波基選擇的困難，而且小波方法不具自適應(yīng)性，一旦基函數(shù)選定后，就只能用來分析處理所有數(shù)據(jù)，這也是一直困擾小波變換研究應(yīng)用的瓶頸。

1.3 Wigner-Ville分布

Wigner-Ville分布屬于2次時頻分析，是信號在時頻平面上的聯(lián)合功率譜，分辨率很高。該方法雖然克服了短時傅里葉變換的一些缺點，但在分析多分量復(fù)雜信號時，存在嚴(yán)重的交叉干擾，目前雖然有了一些消除干擾的方法，但都是以降低分辨率為代價。

1.4正弦曲線擬合

Choi[3]等首先使用正弦曲線擬合的方法來分析語音信號。該方法是利用1組不同頻率的正弦曲線來擬合一時窗內(nèi)信號，取擬合誤差最小的正弦曲線為該信號的最佳擬合曲線。同時給出頻率、振幅和相位信息，取擬合的剩余曲線重新擬合，直到擬合誤差小到一定程度為止。這種方法可以用于多頻率信號，但是該方法的擬合必須是在1個窗口內(nèi)進(jìn)行，因此時間分辨率不高，而且不同頻率成分的信號會相互干擾，如果提取1個頻率成分信號時通常就會破壞其他信號，只有選擇的窗長是半個周期或整數(shù)倍附近時，才對其他信號影響較小。

2 希爾伯特-黃變換的理論和算法

2.1概述

1998年美國華裔科學(xué)家E.Huang提出了1種新型的非線性、非平穩(wěn)信號的處理方法，即希爾伯特-黃變換（HHT）。該方法從信號自身特征出發(fā)，解決了Hilbert變換不能處理多值頻率的限制，首先使用經(jīng)驗?zāi)B(tài)(EMD)方法把信號分解成一系列的本征模態(tài)函數(shù)(IMF)，然后對這些IMF分量進(jìn)行Hilbert變換，從而在時頻平面上得到能量分布的Hilbert時頻譜圖。

該方法主要由“經(jīng)驗?zāi)B(tài)分解(EMD)”和“Hilbert譜”2部分組成。EMD方法自適應(yīng)地將被分析信號分解成反映信號波動模式的若干個本征模態(tài)函數(shù)(IMF)。IMF沒有固定的表達(dá)式，而是根據(jù)信號自身的波動情況自適應(yīng)決定。對EMD的分解結(jié)果再進(jìn)行Hilbert變換，計算瞬時頻率和幅值，由此得到表征信號時間-頻率-幅值的Hilbert譜。HHT算法執(zhí)行效率高，具有良好的局部時頻聚集性和自適應(yīng)性。目前，該方法已經(jīng)廣泛應(yīng)用于機(jī)械工程、通信電子、地球物理、生物醫(yī)藥等許多涉及信號分析的領(lǐng)域。

2.2 EMD分解方法

HHT主要包含EMD分解和Hilbert變換2個過程。根據(jù)信號的時間尺度，EMD方法自適應(yīng)地通過1個“篩選過程”從被分析信號中提取生成其本身固有的1族本征模態(tài)函數(shù)(IMF)，每個IMF必須滿足以下2個條件：

（1）對于1列數(shù)據(jù)，極值點個數(shù)和過零點個數(shù)相等或至多相差1個；

（2）在任意點處，所有局部極大值點確定的上包絡(luò)和由所有局部極小值點確定的下包絡(luò)的均值為0。

根據(jù)這些條件，EMD的篩選過程主要包括以下3個步驟：

（1）首先找出原始信號X(t)的局部極大值和極小值，將其用3次樣條函數(shù)求取信號的上下包絡(luò)線：U(t)和V(t)；

（2）計算上下包絡(luò)線的平均曲線M(t)=[U(t)+V(t)]/2；

（3）用信號X(t)減去M(t)后得到1個去掉低頻的剩余部分h1(t)=X(t)-M(t)。

一般來講，h1(t)仍然不是1個平穩(wěn)數(shù)據(jù)序列，即分解得到的h1(t)并不一定完全滿足IMF條件[4]，用h1(t)代替X(t)，再求出上下包絡(luò)線U1(t)和V1(t)，重復(fù)上述過程，即

……

直到所有的hk(t)滿足IMF條件。此時分解得到的第1個IMF，C1(t)和r1(t)；即

必須指出，過多的重復(fù)該處理過程會導(dǎo)致模態(tài)分量變成純粹的頻率調(diào)制信號，因此，必須確定1個篩選停止準(zhǔn)則。比如限制連續(xù)2次處理結(jié)果之間的標(biāo)準(zhǔn)差σ的大小來實現(xiàn)。

一般選取σ=0.2～0.4。

對信號的剩余部分r1(t)繼續(xù)進(jìn)行EMD分解，直到剩余部分為1個單一信號或者其值小于預(yù)先設(shè)定值時，分解完成。最終原始信號是所有IMF及殘余量之和

式中：Ci為本征模態(tài)函數(shù)；rn為殘余量。2.3 Hilbert變換

對所有IMF進(jìn)行Hilbert變換，則

式中：P為柯西主值。

通過Hilbert變換，Ci（t）和yi（t）可以組成解析信號Zi（t），即

其中

相應(yīng)可以得到瞬時頻率

顯然可知

這里省略了殘余函數(shù)rn，Re表示取實部。展開式（14）稱為Hilbert幅值譜，記為

這樣可以得到瞬時頻率和幅值，可以在時間頻率平面上根據(jù)幅值的強(qiáng)弱繪制時頻譜圖。如將其沿時域進(jìn)行積分，可以進(jìn)一步定義邊際譜

類似傅里葉變換得到的幅頻譜，邊際譜顯示信號能量在頻域內(nèi)的分布，但它們的物理意義完全不同。

2.4 HHT的特點

作為新發(fā)展的1種時頻分析方法，HHT便于理解、掌握和使用。首先區(qū)別于其他方法，它具有自適應(yīng)性，EMD根據(jù)信號的自身特點，直接從信號本身出發(fā)對信號進(jìn)行分解，自適應(yīng)地將信號分解成數(shù)目有限的IMF分量。與短時傅里葉變換相比，它能從時域頻域2方面分析信號的局部特征。與小波變換相比，它沒有小波基的選擇問題，無需根據(jù)信號特點選擇小波基。

另外，HHT具有良好的時頻聚集性。由EMD方法將復(fù)雜信號分解成單分量信號IMF后，對各IMF再進(jìn)行Hilbert變換。通過對解析信號的相位函數(shù)進(jìn)行微分運算來獲得瞬時頻率，最后得到Hilbert譜。在Hilbert變換時，積分運算不是全局性的，而是通過與1/t進(jìn)行卷積，使得它的結(jié)果相當(dāng)?shù)木植炕?。因此，HHT具備良好的時頻局部化特性，或者說具有了良好的時頻聚集性。

2.5 HHT存在的缺陷

作為1種新興的信號分析方法，HHT在理論上還有不成熟的地方。

（1）模態(tài)混疊。所謂模態(tài)混疊就是1個IMF分量包括了尺度差異較大的信號。產(chǎn)生模態(tài)混疊的原因是由某1或某幾個時間尺度的IMF不連續(xù)造成的。而且EMD分解會產(chǎn)生一些低頻的與原信號相關(guān)性較弱的IMF，這些IMF會在Hilbert譜刻畫中產(chǎn)生某些虛假的頻率分量。對此，一些研究提出了某些解決方法，2003年Rilling[5]提出了局部EMD方法；2005年Ryan Deering[6]提出了掩膜信號法等。這些方法雖然在一定程度上改善了分解結(jié)果，但模態(tài)混疊對分析結(jié)果的影響仍很嚴(yán)重，特別是頻率信息復(fù)雜的信號更為突出，是目前制約HHT發(fā)展的瓶頸。要想從根本上解決，還要從理論上入手。

（2）端點效應(yīng)問題。在對原始數(shù)據(jù)取包絡(luò)時，需要對上極點和下極點進(jìn)行樣條插值擬合，然后再平均，在進(jìn)行樣條插值時，除非數(shù)據(jù)的2個端點處就是數(shù)據(jù)的極值點，否則就會產(chǎn)生擬合誤差。在EMD篩選過程中，由于端點處極值的不確定性導(dǎo)致誤差的不斷向下積累，以此類推，隨著分解的進(jìn)行，誤差就會由端點向內(nèi)逐漸傳播，污染內(nèi)部數(shù)據(jù)。目前，有些方法試圖來改善端效應(yīng)對分析結(jié)果的影響，比如使用比較多的極值延拓法，邊界局部特征延拓法，鏡像延拓法等，但這些方法只能在一定程度上降低或改善端效應(yīng)對分析結(jié)果的影響，不能從根本上消除。

（3）曲線擬合問題。HHT首先通過極值點擬合信號包絡(luò)線，再通過EMD方法篩選IMF，這里就涉及到插值擬合問題。HHT提出使用3次樣條插值法雖然具有良好的光滑性，但也容易造成過沖和欠沖。因此后續(xù)研究者相繼提出了不少的修整和改進(jìn)。比如高次樣條插值法、多項式擬合法、阿克瑪插值法等。

3 基于Labview的HHT實現(xiàn)

Labview是由美國NI公司開發(fā)的在業(yè)界領(lǐng)先的工業(yè)標(biāo)準(zhǔn)圖形化編程工具，主要用于開發(fā)測試、測量和控制系統(tǒng)。它將軟件和各種不同的測量儀器硬件及計算機(jī)集成在一起，建立虛擬儀器系統(tǒng)，形成用戶自定義的解決方案[7]。

基于Labview直觀易學(xué)、接近人類思維、無需編寫繁瑣代碼的特點，本文將上述HHT算法在Labview平臺上實現(xiàn)。系統(tǒng)的程序流程如圖1所示。

直接利用EMD方法分析計算遠(yuǎn)遠(yuǎn)不夠，還要考慮到EMD方法造成的端點效應(yīng)和模態(tài)混疊等問題[8]。本程序采用對稱延拓法來消除端點效應(yīng)[9]，為了抑制模態(tài)混疊和減少虛假IMF分量個數(shù)，本程序還要對每次迭代過程的極值點個數(shù)及距離等進(jìn)行控制和有限插值[10]。

4 仿真信號試驗驗證

現(xiàn)利用某一典型的線性調(diào)頻信號sin 2π·（10+40·t）·t和1個正弦信號sin 2π·100·t，t∈[0s,1s]相疊加，采樣頻率為2 kHz，采樣長度為1 s。對HHT算法程序的EMD分解和時頻分析可靠性進(jìn)行檢驗。顯然，此疊加信號的頻率變化范圍在10～50 Hz，同時還有1個頻率為100 Hz的正弦信號。仿真信號如圖2所示。

圖1 Labview平臺實現(xiàn)的HHT程序流程

圖2 仿真疊加信號

經(jīng)過HHT算法程序得到的仿真信號Hilbert時頻譜如圖3所示。圖中曲線顏色的深淺代表幅值的強(qiáng)弱。與圖3的傅里葉變換結(jié)果相比，HHT將原始仿真信號的頻率和幅值隨時間的變化規(guī)律同時顯示在時頻譜中。而傅里葉變換（如圖4所示）只能單純給出頻域特征，不能同時給出時間和頻率的聯(lián)合信息，對非線性、非平穩(wěn)類信號無能為力。

圖3 仿真信號經(jīng)HHT得到的時頻譜

圖4 仿真信號的傅里葉變換幅頻譜

5 某型發(fā)動機(jī)氣動失穩(wěn)數(shù)據(jù)時頻分析

目前，HHT已經(jīng)開始應(yīng)用于航空發(fā)動機(jī)氣動不穩(wěn)定性研究，如發(fā)動機(jī)或壓氣機(jī)試驗的旋轉(zhuǎn)失速和喘振、動態(tài)畸變試驗和壓氣機(jī)轉(zhuǎn)子性能等。

某型發(fā)動機(jī)臺架試車時出現(xiàn)的氣動失穩(wěn)信號如圖5所示。對其進(jìn)行普通的傅里葉變換，如圖6所示。

圖5 某型發(fā)動機(jī)臺架試車時出現(xiàn)的氣動失穩(wěn)仿真信號

圖6 氣動失穩(wěn)仿真信號的傅里葉幅頻譜

在除去直流分量，只保留交流分量的傅里葉幅頻譜上可觀察到，信號主要含有20 Hz和40 Hz 2個頻率分量，由傅里葉變換的特點，2個頻率分量出現(xiàn)的時間點不得而知。

現(xiàn)采用EMD方法對原始信號進(jìn)行分解。一共分解出7個經(jīng)驗?zāi)B(tài)函數(shù)（IMF），這里略去后面幾個趨于直流的IMF分量,只給出前3個IMF分量，如圖7所示，并做出Hilbert時頻譜如圖8所示。

圖7 EMD方法得到的IMF分量

圖8 氣動失穩(wěn)仿真信號的Hilbert時頻譜

從圖8中可見，此氣動失穩(wěn)仿真信號包含的20Hz和40 Hz 2個主要頻率分量，基本上同時出現(xiàn)，并貫穿整個信號始終，在信號的結(jié)尾0.92 s處開始出現(xiàn)較高頻分量，這從EMD分解的幾個IMF結(jié)果中也可以看出。根據(jù)EMD分解性質(zhì)，各個IMF分量的局部頻率依次降低。本次分解結(jié)果的IMF1基本上是40 Hz的頻率分量，在0.92 s處開始過渡到較高頻分量；同樣在IMF2中基本上是20 Hz的頻率分量，在0.92 s處開始過渡到40 Hz的頻率分量；在IMF3以后的分量中信號的直流成分開始占主導(dǎo)。

在信號的起止端，由于數(shù)據(jù)的截斷造成了時頻譜的端點效應(yīng)問題（又稱端點飛翼），所以時頻譜圖在時間起點和終點處的分析結(jié)果不準(zhǔn)確，如果要求分析此處數(shù)據(jù)，可以適當(dāng)延長被分析數(shù)據(jù)的長度。另外，圖中各頻率分量抖動的部分是由HHT固有的模態(tài)混疊造成，可以通過極值加密、EEMD法、改進(jìn)掩膜信號法等得到改善，但是基于HHT的實際應(yīng)用超前理論的研究和發(fā)展，若不從理論算法上改進(jìn)完善，這些問題是不能從根本上得到解決。

6 結(jié)論

（1）運用圖形化編程語言Labview編制HHT算法，開發(fā)過程更加簡捷，程序界面更加友好。在運行時間上與Matlab程序相比獲得了性能提升，這源于Labview程序所具有的并行特性。

（2）使用HHT方法，定位出不同頻率分量在信號中出現(xiàn)的不同時刻，為氣動失穩(wěn)這種非平穩(wěn)、非線性信號的分析探索出新的途徑，同時還驗證了HHT的一些優(yōu)勢和不足。

（3）針對使用這種新型的信號時頻分析技術(shù)處理實際信號時，HHT存在的固有問題[11]，本文采取了一些有效措施。

（4）HHT擺脫了傳統(tǒng)的傅里葉分析的局限，是一種簡單易行，且有效的信號分析方法。將HHT變換應(yīng)用在發(fā)動機(jī)氣動失穩(wěn)信號的分析具有適用性。但HHT是新興的時頻分析方法，有待進(jìn)一步研究和改進(jìn)。

[1]管霖，吳國沛，黃雯瑩．小波變換在電力設(shè)備故障診斷中的應(yīng)用研究[J]．中國電機(jī)工程學(xué)報，2000，20（10）：46-49.

GUAN Lin,WU Guopei,HUANG Wenying．Studies on the application of wavelet transform in fault diagnosis of electric devices[J]．Journal of Chinese Electrical Engineering Science, 2000,20（10）：46-49．(in Chinese)

[2]彭志科，何永勇，盧青．用小波時頻分析方法研究發(fā)電機(jī)碰摩故障特征[J].中國電機(jī)工程學(xué)報，2003，23（5）：75-79．

PENG Zhike,HE Yongqing,LU Qing．Using wavelet method to analyze fault features of rubrotor in Generator[J]．Journal of Chinese Electrical Engineering Science,2003，23（5）：75-79. (in Chinese)

[3]Choi A.Real-time fundamental frequency estimation by least-square fitting[J].IEEE Trans.Speech Audio Processing，1997（5）：201-205.

[4]HUANGDaji,ZHAOJinping,SUJilan.Pracitical implementation of Hilbert-Huang Transform algorithm[J].Acta Oceanologica Sinica,2003,22（1）：1-14.

[5]Rilling D,Flandrin P,Gonc P,et al.On empirical mode decomposition and its algorithms[C]//Grado,Italy:IEEE，2003：1-5.

[6]Ryan D,James F K.The use of amasking signal to improve empiricalmode decomposition[C]//Philadelphia:International Conference on Acoustics Speech,and Signal Processing(ICA SSP),2005：485-488.

[7]陳錫輝，張銀鴻.Labview程序設(shè)計從入門到精通[M].北京：清華大學(xué)出版社,2007：1-3．

CHEN Xihui,Zhang Yinhong．Labview Program design junior to superior[M].Beijing:Tsinghua University Publishing House, 2007：1-3．(in Chinese)

[8]胡維平，莫家玲，杜明輝.經(jīng)驗?zāi)B(tài)分解中的頻域分辨率及其改進(jìn)法[J].華南理工大學(xué)學(xué)報，2007（5）：2-3．

HU Weiping,MO Jialing,Du Minghui.Method for improving frequency resolution of empirical mode decomposition[J]．Journal of South China University of Technology,2007（5）：2-3.(in Chinese)

[9]陳忠，鄭時雄．EMD信號分析方法邊緣效應(yīng)的分析[J]．?dāng)?shù)據(jù)采集與處理，2003，18（1）：2-3.

CHEN Zhong,ZHENG Shixiong．Analysis on end effects of EMD method[J].Journal of Data Acquisition and Processing, 2003,18（1）：2-3.(in Chinese)

[10]毛煒，金榮洪，耿軍平．一種基于改進(jìn)Hilbert-Huang變換的非平穩(wěn)信號時頻分析法及其應(yīng)用[J]．上海交通大學(xué)學(xué)報，2006（5）：1-2．

MAOWei,JIN Ronghong,GENG Junping．A time-friquency analysis for non-stationary signals based on improved Hilbert-Huang Transform and it's application[J].Journal of ShangHai JiaoTong University,2006（5）：1-2．(in Chinese)

[11]于德介，程軍圣，楊宇.機(jī)械故障診斷的Hilbert-Huang變換方法[M].北京：科學(xué)出版社，2006：46-51．

YU Dejie,CHENG Junsheng,YANG Yu.Hilbert-Huang Transform method formechanical fault diagnosis[M].Beijing: Science Publishing House,2006：46-51.(in Chinese)

Time-frequency Analysis Method of Aeroengine Aerodynam ic Instability Signal Based on Hilbert-Huang Transform

GUO Gui-xi,WANG Zhen-hua,WANG Hong-xiang,GAO Hui-guo, XUEXiu-sheng
（AVIC Shenyang Engine Design and Research Institute,Shenyang 110015,China）

Aiming at aerodynamic instability in aeroengine work scope,the time varying signal characteristics were analyzed by Hilbert-Huang Transform(HHT).The analysis new method of nonlinear and non-stationary signal was compared with other timefrequency analysis methods.The problems of HHT were solved to adapt the data analysis requirements of aeroengine aerodynamic instability.HHT method was implemented and verified on Labview platform by the numerical simulation test.The results show the method is accurate and effective and the project test data processing of aeroengine is presented.The HHT is feasibility and adaptability to solve relative problems.There are some little defection for HHTmethod also.

time-frequency analysis;Hilbert-Huang Transform;aeroengine;aerodynamic instability

2012-12-12

郭貴喜（1980），男，碩士，工程師，從事航空發(fā)動機(jī)測試工作。