基于希爾伯特-黃與小波變換的壓氣機失速信號分析方法

向宏輝，侯敏杰，楊榮菲，葛寧，劉志剛，吳晨

（1.南京航空航天大學能源與動力學院，南京210016；2.中國燃氣渦輪研究院航空發動機高空模擬技術重點實驗室，四川江油621703）

基于希爾伯特-黃與小波變換的壓氣機失速信號分析方法

向宏輝1，2，侯敏杰2，楊榮菲1，葛寧1，劉志剛2，吳晨2

（1.南京航空航天大學能源與動力學院，南京210016；2.中國燃氣渦輪研究院航空發動機高空模擬技術重點實驗室，四川江油621703）

針對多級跨聲速壓氣機失速過程壁面靜壓測量數據，分別采用希爾伯特-黃變換（HHT）和小波分析方法進行時頻分析。相比于小波分析方法，HHT時頻圖中失速過程的頻率特征變化更為明顯，且利用HHT圖中失速發生后的特征頻率可判斷出失速團的旋轉速度。結合方差分析，進一步探索了獲取失速先兆信號的可行性。結果表明，兩種方法均可用于壓氣機失速先兆預測，其中HHT方法能夠提供更高層次的分解信號，得到的失速先兆更為提前。

壓氣機；失速信號；失速先兆；希爾伯特-黃變換；小波變換；方差分析；時頻分析

1　引言

壓氣機的失速與喘振直接影響航空發動機的穩定工作范圍，主動控制技術是實現發動機擴穩的重要技術手段。如果在失速完全發生后進行主動控制，由于流體慣性作用，往往需要很大的輸入能量及較長的控制時間才能成功退出失速。反之，如果能在失速發生時刻前及時啟動主動控制系統，則能在較短時間以較小能量（如葉尖微射流技術）實現失速抑制，確保航空發動機安全可靠工作。因此，如何盡早地獲取與識別失速先兆信號，是壓氣機失速主動控制研究的一個重要內容。

壓氣機失速信號分析通常包括信號選取和先兆捕捉兩方面內容。其中，選取的信號可以是原始采集信號，也可以是濾波處理信號。原始信號由于混入大量背景噪聲，相對較弱的失速信號容易淹沒其中，影響失速先兆捕捉精度。因此，失速信號分析過程中，首先需對原始采集信號進行濾波降噪處理。小波分析是應用較廣的一種降噪方法，其降噪結果依賴于小波基函數與分解層次的選取，壓氣機失速信號通常選用Daubechies正交小波系，采用4至5層信號分解［1-2］。除此外，黃鍔等［3-4］提出的經驗模態分解（EMD）也是一種有效的濾波降噪方法。將此方法與希爾伯特譜分析相結合，得到了著名的時頻分析方法，即希爾伯特-黃變換（HHT）。該方法相較于小波分析，最大優勢在于分解結果不依賴于基函數，能自適應處理帶噪聲的非線性非穩態信號，頻率分辨率較高。基于此，HHT在多個工程領域得到廣泛采用，如機械故障診斷［5］、地震信號分析［6］等。在壓氣機失速信號分析方面，張勇等［7］利用HHT時頻分析功能研究了小波降噪后離心壓縮機試驗數據，得到了壓縮機喘振頻率；郭貴喜等［8］利用HHT分析了軸流壓氣機失速信號，得到了失速前、后頻率的變化特征；李長征等［9］利用HHT分析了小型壓氣機出口總壓信號，通過判斷時域圖中頻率突變點發生時刻，得到了壓氣機失速頻率及失速起始時間；而符嬈等［10］利用HHT分析了三級軸流壓氣機失速過程的小波降噪信號，通過預先設定閾值來監測失速先兆頻帶與失速起始時刻。上述研究主要集中在時頻分析結果方面，而在失速先兆捕捉方面，將HHT分析方法與失速判斷方法相結合的研究還未見公開報道。

在失速先兆捕捉方面，可分為兩類方法：利用降噪信號的階躍點進行失速判別與利用失速特征參數的階躍點進行失速判別。其中，失速特征參數包括統計參數（如方差［1，11］、自相關系數［12-13］、互相關系數［13-14］）和混沌特征參數（如關聯積分［15-16］、最大Lyapunov指數［17］、模糊熵［18］等）。混沌特征參數計算需首先進行相空間重構，會增加失速計算過程的時間與復雜程度；而統計參數僅針對信號進行直接處理，計算時間短。馬翼等［1］對小波降噪信號進行方差分析提前獲得了壓氣機失速先兆；高曼等［11］對比分析了原始失速信號的方差分析結果與自相關系數分析結果，認為亞聲速壓氣機中自相關分析法得到的失速起始時間早于方差分析結果；李繼超等［13］指出互相關系數法在失速先兆捕捉方面優于自相關系數法。

為探索高負荷多級軸流壓氣機動態壓力信號分析與失速先兆捕捉的適用方法，本文以某六級跨聲速高壓壓氣機失速試驗數據為研究對象，首先對比HHT與小波時頻分析結果，然后結合失速先兆捕捉方法，分析EMD降噪與小波降噪用于提前獲取壓氣機失速先兆的可行性。本文研究結果有望為高負荷跨聲速壓氣機失速先兆信號及時捕獲與主動控制策略實施提供技術基礎。

2　試驗樣本數據

動態壓力試驗樣本數據來源于某帶可調進口導葉的六級軸流高壓壓氣機［19］。該壓氣機設計轉速為10 512 r/min，第1～2級轉/靜子葉片數分別為47、68、57、82。試驗時，壓氣機進口保持均勻來流條件，分別在第1級轉/靜子之間、第2級轉/靜子之間的機匣壁面沿周向布置2支動態壓力傳感器。從壓氣機進口往后看，沿逆時針方向，第1級壁面2點脈動靜壓周向安裝角度分別為50°、140°（分別用S1-P1、S1-P2表示），第2級壁面2點脈動靜壓周向安裝角度分別為42°、86°（分別用S2-P1、S2-P2表示），靜壓測點采樣頻率為5 kHz。

試驗用傳感器選用美國PCB公司生產的M101A05型壓電式動態壓力傳感器（固有頻率為250 kHz），并采用半無限長管測試技術解決管腔的駐波效應。動態數據采集系統為Genesis系統，采樣率6.4 kS/s，分析帶寬2.5 kHz。試驗中，通過調節試驗器排氣節氣門開度改變壓氣機出口反壓，實時記錄壓氣機95%設計轉速下，從穩定工作狀態連續變化到失速工作狀態的機匣壁面脈動靜壓信號。

3　信號處理方法

3.1HHT方法

HHT方法主要包括EMD和希爾伯特譜分析兩部分。其中EMD類似于小波分解，通過篩選算法將原始信號分解為一系列不同尺度的本征模態函數（IMF），原始信號等于所有IMF分量與殘差之和。希爾伯特譜分析類似于小波時頻分析中的譜分析，得到的是信號瞬時頻率，而不是小波譜得到的一小段時間窗口的計算頻率。相比小波分析，希爾伯特譜分析的頻率分辨率更高。

HHT分析的核心是EMD，利用時間序列的上下包絡線平均值對信號進行篩選。篩選的終止條件為：信號極值點的數量與零點數相等或最多相差1個，信號由極大值點和由極小值點構成的上下包絡線平均值為零。EMD具體分解過程如下：

（1）對原始信號X（t）的極大值點、極小值點采用三次樣條函數插值，構造原始信號的上下包絡線XU（t）、XL（t），以及包絡線平均值M1=（XU（t）+XL（t））/2。

（2）考察H1=X（t）-M1是否滿足篩選條件，如果滿足則轉到下一步，否則將H1作為原始信號繼續操作步驟（1），求得M11及H11=H1-M11；依此篩選，直至第k步H1k=H1k-1-M1k滿足篩選條件，則第一個IMF分量C1=H1k。

（3）將第一個殘留值R1=X（t）-C1作為原始信號，按照上述方法繼續篩選，得到第二個IMF分量C2（t）；依此類推，得到一系列IMF分量，直至Rn為單調信號或只存在一個極點位置。

由上述分解過程可知，EMD得到的IMF信號為不同分解層次的降噪信號，類似于小波分解重構信號，理論上可用于壓氣機失速先兆分析。另外，對不同分解層次的IMF信號進行希爾伯特譜分析，還可得到頻率-幅值隨時間的變化。

3.2方差分析法

時間序列Cj（t）的方差計算公式見式（1）。其中，N為計算窗口內的數據長度，tN為失速監測時刻。將Cj（t）進行分段滑動，滑動寬度為ΔN，則得到方差隨時間的變化關系。

針對試驗樣本數據中的S1-P1原始信號，選取N為100，ΔN為5，計算得到方差時間序列如圖1所示（圖中縱坐標V代表方差）。可見，方差在穩定工作狀態時基本為常值，接近失速時出現明顯突變。定義突變點發生時刻為失速起始時刻。本文定義當前方差值為前100個方差均值的2倍時壓氣機開始進入失速，圖中虛線為根據此判據得到的失速起始時刻，可見該方法能有效捕捉失速起始時刻。

圖1　原始信號的方差Fig.1 Variance of raw signal

圖2　S1-P1測量點信號的時頻分析結果Fig.2 Time frequency diagram of Hilbert-Huang transform and wavelet at S1-P1

4　基于HHT與小波方法的失速信號分析

4.1時頻分析

采用HHT方法與小波方法對失速信號進行時頻分析，可得到信號頻率與幅值隨時間的變化關系，圖2給出了S1-P1測量點的時頻分析結果。對比發現：HHT方法得到的頻率分辨率高于小波方法，在失速發展過程中，HHT時頻圖存在幅值占優的一個或幾個頻率，而小波時頻圖中幅值占優的頻率為一定范圍，不能明確與失速過程相關的主要特征頻率，不利于失速發展過程中干擾頻率判斷。HHT時頻圖中504.40 s左右出現一個約10 Hz的低頻脈動，且隨著時間推移，脈動頻率逐漸增加并穩定在78 Hz左右，此頻率增加過程對應著從小尺度低頻擾動逐漸發展為大尺度中頻失速團的演化過程。失速發生后，對比轉子轉動一周產生的166.5 Hz轉動頻率，HHT時頻圖中78 Hz頻率占轉動頻率的46.8%，可認為是失速團通過頻率，由此得到失速團旋轉速度為46.8%轉子轉速。由圖3中S1-P1、S1-P2測量點失速信號峰值點確定的失速團旋轉速度為47.6%轉子轉速（圖中縱坐標ΔP代表壓力脈動幅值），與HHT圖中失速點頻率判斷出的失速團旋轉速度相近。因此，采用單點HHT分析即可獲得失速團旋轉速度，相較于傳統測量方案中利用周向兩點測量確定失速團旋轉速度而言，可以減少測量點數，降低試驗成本。

圖3　不同周向測量位置原始信號的時域圖Fig.3 Time domain diagram of raw signal at S1-P1 and S1-P2

4.2失速先兆監測

首先利用HHT與小波方法對采集信號進行多重分解，獲得失速先兆計算的待分析信號。圖4給出了S1-P1測量點兩種方法的分解結果對比，其中C0表示原始信號，C1～C8分別代表分解得到的IMF分量或小波分解重構信號分量。可見，兩種方法得到的分解信號在高頻率時基本重合，隨著分解層次的增加，信號中的突變點前移，如果將信號劇烈變化的時刻看成失速起始時刻，則信號分解后計算得到的失速先兆時間相較于原始信號提前。但當分解層次過高時，如第8層，分解信號逐漸變為定常時均信號，信號中不再包含低頻周期性失速信號，由此計算得到的失速先兆可能與真實結果相矛盾。

進一步將上述分解信號與失速先兆捕捉方法相結合，計算不同分解信號失速先兆起始時間。圖5給出了各分解信號失速先兆的方差判斷結果，從中可以看出，對于兩種信號分解方法，隨著信號分解層次的增加，方差判斷得到的失速先兆時間提前。對于低分解層次，高頻信號預測得到的失速起始時間相同；對于中等分解層次得到的信號（如C4～C7），基于IMF分量預測得到的失速時間略微早于基于小波分解預測得到的失速時間；對于高層次分解信號（如C8～C11），兩種方法預測得到的失速起始時間呈現出相反的變化趨勢。由圖4可知，C8或更高層次分解信號的脈動逐漸接近于時均信號，這會對具有低頻信號特征的失速先兆預測帶來誤判。

圖4　本征模態函數分量與小波分解重構信號對比Fig.4 Comparison of intrinsic mode function and wavelet decomposition signal

圖5　S1-P1測量點處本征模態函數分量與小波分解重構信號的失速先兆Fig.5 Stall precursor obtained from intrinsic mode function and wavelet decomposition signal at S1-P1

按照相同的分析步驟，對其他測量點（S1-P2、S2-P1和S2-P2）進行失速先兆判斷，得到不同位置測量點的失速先兆隨信號分解層次的變化規律，如圖6所示。可見，隨著分解層次的增加，不同測量位置的失速信號采用HHT方法均能獲得一個最低的失速起始時間，而采用小波分解分析方法獲得的最低起始時間所對應的分解層次，會隨著測量點位置的變化而變化。因此，采用小波分解方法進行失速判斷時，建議取第4～6層分解信號進行失速分析，可避免測量點位置對失速先兆判斷的影響；而采用HHT方法時，建議取第6～8層分解信號進行失速分析，可獲得最早的失速起始時間。相比之下，HHT方法獲得的失速起始時間，要早于小波分析獲得的失速起始時間，該結論在其他轉速下也成立。

圖6　不同測量點位置處本征模態函數分量與小波分解重構信號的失速先兆Fig.6 Stall precursor obtained from intrinsic mode function and wavelet decomposition signal at different measurement location

5　結論

（1）HHT方法相對于小波方法頻率分辨率更高，從HHT時頻圖中可清楚地顯示出壓氣機從小尺度低頻擾動信號向大尺度中頻失速團發展過程的頻率變化。

（2）HHT時頻圖中失速發生后產生的頻率為失速團旋轉頻率，利用該方法判定失速團旋轉速度僅需一個測量點即可，減少了試驗測量布點數量。

（3）采用HHT分解與小波分解重構獲得的不同分解層次信號，結合方差分析法，均能提前預測壓氣機的失速先兆。相比于小波方法，HHT分解信號得到的失速先兆時間更為提前。

［1］馬翼，陸鳳霞.基于Daubechies小波的壓氣機失速分析及預測［J］.機械工程與自動化，2011，40（4）：123—125.

［2］姜濤，賈智偉，李應紅.基于小波分析的壓氣機失速故障檢測［J］.推進技術，2003，24（1）：55—57.

［3］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 London A，1998，454：903—995.

［4］Wu H，Huang N E.On the filtering properties of the empirical mode decomposition［J］.Advances in Adaptive Data Analysis，2010，2（4）：397—414.

［5］Wu F J，Qu L S.Diagnosis of subharmonic faults of large rotating machinery based on EMD［J］.Mechanical Systems and Signal Processing，2009，23：467—475.

［6］楊培杰，印興耀，張廣智.希爾伯特-黃變換地震信號時頻分析與屬性提取［J］.地球物理學進展，2007，22（5）：1585—1590.

［7］張勇，張春梅.用HHT變換處理離心壓縮機喘振試驗數據［J］.流體機械，2012，40（1）：10—12.

［8］郭貴喜，王振華，王洪祥，等.基于HHT的航空發動機氣動失穩信號時頻分析方法［J］.航空發動機，2013，39（3）：62—67.

［9］Li C Z，Xiong B.Application of HHT for aerodynamic instability signal analysis in compressor［C］//.3rdInternational Conference on Signal Processing Systems.Singapore：2011.

［10］符嬈，雷勇.基于HHT的壓氣機氣動失穩檢測技術研究［J］.測控技術，2010，29（S0）：71—75.

［11］高曼，楚武利，吳艷輝，等.亞聲速軸流壓氣機失速監測方法試驗［J］.推進技術，2012，33（3）：398—404.

［12］Nakajima T，Tahara N，Kurosaki M.Stall prediction methodforaxialflowcompressor：US，6715984B2［P］. 2003-01-09.

［13］李繼超，童志庭，聶超群，等.基于互相關分析的前失速先兆檢測分析［J］.航空學報，2013，24（1）：28—36.

［14］Salunkhe P B，Joseph J，Pradeep A M.Active feedback control of stall in an axial flow fan under dynamic inflow distortion［J］.Experimental Thermal and Fluid Science，2011，35（6）：1135—1142.

［15］Tahara N，Nakajima T，Kurosaki M，et al.Active stall control with practicable stall prediction system using autocorrelation［R］.AIAA 2001-3623，2001.

［16］吳云，李應紅，張樸.某渦噴發動機壓氣機氣動失穩過程的非線性分析［J］.航空動力學報，2005，20（5）：890—894.

［17］楊榮菲，侯安平，周盛，等.最大Lyapunov指數用于壓氣機局部性能的分析［J］.航空學報，2009，30（4）：622—624.

［18］劉志剛，楊榮菲，向宏輝.基于模糊熵的壓氣機失速信號分析研究［J］.測控技術，2014，33（2）：53—56.

［19］向宏輝，侯敏杰，葛寧，等.多級環境下徑向總壓畸變影響效應的試驗評價［J］.航空動力學報，2012，27（6）：1270—1277.

Analysis method investigation of axial compressor stall signal based on Hilbert-Huang transform and wavelet transform

XIANG Hong-hui1，2，HOU Min-jie2，YANG Rong-fei1，GE Ning1，LIU Zhi-gang2，WU Chen2
（1.College of Energy and Power Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2.Key Laboratory on Aero-Engine Altitude Simulation Technology，China Gas Turbine Establishment，Jiangyou 621703，China）

A comparative study of Hilbert-Huang transform（HHT）and wavelet was presented through analyzing the wall static pressure acquired from a transonic multistage compressor which operated from stable to stall point continuously.Firstly，the time-frequency chart was compared.It is shown that the changes of frequency get more recognizable by HHT than by the method of wavelet in the development process of compressor stall.And the rotation speed of stall cells can be calculated based on the frequency that taking from HHT time-frequency chart at the time of compressor undergoing stall.Secondly，feasibility study on the stall precursor capture method was carried out which includes two steps of signal decomposition and variances calculation.The results suggest that the decomposition signal by HHT and wavelet can capture stall inception earlier than raw signal.As HHT can provide high levels of decomposition signal for stall prediction，it obtained earlier stall precursor than wavelet decomposition.

compressor；stall signal；stall precursor；Hilbert-Huang transform；wavelet transform；variance analysis；time-frequency analysis

V231.3

A

1672-2620（2015）05-0001-05

2015-04-28；

2015-10-25

江蘇省普通高校研究生科研創新計劃項目（KYLX15_0260）

向宏輝（1979-），男，湖南沅陵人，高級工程師，博士研究生，主要從事葉輪機性能評定與試驗技術研究工作。