基于HHT的航空發動機氣動失穩信號時頻分析方法

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

（中航工業沈陽發動機設計研究所，沈陽110015）

基于HHT的航空發動機氣動失穩信號時頻分析方法

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

（中航工業沈陽發動機設計研究所，沈陽110015）

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

時頻分析；希爾伯特-黃變換；航空發動機；氣動失穩

0 引言

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

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

1 傳統的時頻分析方法

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

1.1短時傅里葉變換

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

1.2小波變換

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

1.3 Wigner-Ville分布

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

1.4正弦曲線擬合

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

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

2.1概述

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

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

2.2 EMD分解方法

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

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

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

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

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

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

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

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

……

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

必須指出，過多的重復該處理過程會導致模態分量變成純粹的頻率調制信號，因此，必須確定1個篩選停止準則。比如限制連續2次處理結果之間的標準差σ的大小來實現。

一般選取σ=0.2～0.4。

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

式中：Ci為本征模態函數；rn為殘余量。2.3 Hilbert變換

對所有IMF進行Hilbert變換，則

式中：P為柯西主值。

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

其中

相應可以得到瞬時頻率

顯然可知

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

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

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

2.4 HHT的特點

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

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

2.5 HHT存在的缺陷

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

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

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

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

3 基于Labview的HHT實現

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

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

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

4 仿真信號試驗驗證

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

圖1 Labview平臺實現的HHT程序流程

圖2 仿真疊加信號

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

圖3 仿真信號經HHT得到的時頻譜

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

5 某型發動機氣動失穩數據時頻分析

目前，HHT已經開始應用于航空發動機氣動不穩定性研究，如發動機或壓氣機試驗的旋轉失速和喘振、動態畸變試驗和壓氣機轉子性能等。

某型發動機臺架試車時出現的氣動失穩信號如圖5所示。對其進行普通的傅里葉變換，如圖6所示。

圖5 某型發動機臺架試車時出現的氣動失穩仿真信號

圖6 氣動失穩仿真信號的傅里葉幅頻譜

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

現采用EMD方法對原始信號進行分解。一共分解出7個經驗模態函數（IMF），這里略去后面幾個趨于直流的IMF分量,只給出前3個IMF分量，如圖7所示，并做出Hilbert時頻譜如圖8所示。

圖7 EMD方法得到的IMF分量

圖8 氣動失穩仿真信號的Hilbert時頻譜

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

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

6 結論

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

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

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

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

[1]管霖，吳國沛，黃雯瑩．小波變換在電力設備故障診斷中的應用研究[J]．中國電機工程學報，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]彭志科，何永勇，盧青．用小波時頻分析方法研究發電機碰摩故障特征[J].中國電機工程學報，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程序設計從入門到精通[M].北京：清華大學出版社,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].華南理工大學學報，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信號分析方法邊緣效應的分析[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]毛煒，金榮洪，耿軍平．一種基于改進Hilbert-Huang變換的非平穩信號時頻分析法及其應用[J]．上海交通大學學報，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]于德介，程軍圣，楊宇.機械故障診斷的Hilbert-Huang變換方法[M].北京：科學出版社，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），男，碩士，工程師，從事航空發動機測試工作。