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

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

（1.南京航空航天大學(xué)能源與動力學(xué)院，南京210016；2.中國燃?xì)鉁u輪研究院航空發(fā)動機高空模擬技術(shù)重點實驗室，四川江油621703）

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

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

（1.南京航空航天大學(xué)能源與動力學(xué)院，南京210016；2.中國燃?xì)鉁u輪研究院航空發(fā)動機高空模擬技術(shù)重點實驗室，四川江油621703）

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

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

1　引言

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

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

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

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

2　試驗樣本數(shù)據(jù)

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

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

3　信號處理方法

3.1HHT方法

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

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

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

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

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

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

3.2方差分析法

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

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

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

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

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

4.1時頻分析

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

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

4.2失速先兆監(jiān)測

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

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

圖4　本征模態(tài)函數(shù)分量與小波分解重構(gòu)信號對比Fig.4 Comparison of intrinsic mode function and wavelet decomposition signal

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

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

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

5　結(jié)論

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

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

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

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

［2］姜濤，賈智偉，李應(yīng)紅.基于小波分析的壓氣機失速故障檢測［J］.推進(jìn)技術(shù)，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］楊培杰，印興耀，張廣智.希爾伯特-黃變換地震信號時頻分析與屬性提?。跩］.地球物理學(xué)進(jìn)展，2007，22（5）：1585—1590.

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

［8］郭貴喜，王振華，王洪祥，等.基于HHT的航空發(fā)動機氣動失穩(wěn)信號時頻分析方法［J］.航空發(fā)動機，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的壓氣機氣動失穩(wěn)檢測技術(shù)研究［J］.測控技術(shù)，2010，29（S0）：71—75.

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

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

［13］李繼超，童志庭，聶超群，等.基于互相關(guān)分析的前失速先兆檢測分析［J］.航空學(xué)報，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］吳云，李應(yīng)紅，張樸.某渦噴發(fā)動機壓氣機氣動失穩(wěn)過程的非線性分析［J］.航空動力學(xué)報，2005，20（5）：890—894.

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

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

［19］向宏輝，侯敏杰，葛寧，等.多級環(huán)境下徑向總壓畸變影響效應(yīng)的試驗評價［J］.航空動力學(xué)報，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

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

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