基于經驗模態分解的發動機脈動壓力數據分析

任春雷，周小陳，張炳誠，薛小龍，武艷奎，楊懿

（1.國防科技大學空天科學學院，湖南 長沙 410000；2.中國航天科工飛航技術研究院，北京 100074；3.中國航發北京航科發動機控制系統科技有限公司，北京 102200；4.北京控制與電子技術研究所，北京 100038；5.龐巴迪（天津）航空服務有限公司，天津 300300；6.北京航天試驗技術研究所，北京 100074）

0 引言

在液體火箭發動機試驗中，根據所測壓力信號隨時間的變化可分為穩態壓力和脈動壓力兩大類。穩態壓力的測量結果是壓力平均值，用于研究試驗臺和發動機的性能。脈動壓力則主要針對發動機推力室、泵、燃氣發生器、推進劑管路等高壓、高溫、強振動測點進行脈動壓力參數的測量與分析，其主要目的是研究發動機啟動特性、不穩定燃燒、管路特性以及分析發動機的故障原因［1］。由于脈動壓力測量與穩態壓力測量在測量方法、測量系統構成和數據結構等方面有較大的差別。因此采用準確、有效的方法對脈動壓力數據進行分析是非常重要的。

從時間域內壓力變化趨勢的角度進行分析，脈動壓力與穩態壓力最大的區別在于其變化時間序列具有非線性、隨機性的復雜特征［2］。業內一般采用計算一段時間內的平均值或者觀察某個時間點的壓力值的方法對穩態壓力數據進行分析。而對于高頻采樣下的高精度脈動壓力數據，則不僅需要關注其在時域內的幅值信息，而且還需要根據信號時域和頻域的特征，結合實際研究的需要，采用時頻分析技術深度挖掘數據在時域和頻域的特性信息。目前，在信號時頻域分析技術領域，較為常用的方法有傅里葉變換、短時傅里葉變換、快速傅里葉變換和小波分析方法等等。這幾種時頻分析方法各具特點。

傅里葉變換和快速傅里葉變換方法建立了信號在時、頻域的相互轉換關系。通過繪制信號在頻域內的全局圖譜，可以清晰觀察所分析的信號在時域內的頻率-幅值特征。該方法的不足之處在于缺乏局域分析能力，即無法定位所得到的頻率-幅值特性在時域內的發生時刻［3-4］。短時傅里葉變換的核心思想是將一個非線性信號看作若干個短的線性信號的集合，通過設定一個尺度不變的窗函數，將窗函數在時間軸上進行平移，以此來覆蓋整個非線性信號，分析多個短信號局部的時間-頻率信息。該方法的不足之處在于實際應用中需要根據信號的特征選擇對應的窗函數，且分析效果取決于所選定的窗函數，分析的結果具有一定的不確定度性［5］。小波分析方法的實質是通過分析特定時間域內原始信號的特征，人為選定與之相匹配的基函數。將基函數在尺度上進行相應的伸縮，建立一個函數簇。將函數簇在指定的間域內進行多次平移，并對分解得到信息進行濾波等處理后進行重構，能夠得到原始信號在時間域內局部特征信息。小波分析克服傅里葉變換不能兼顧時間和頻率分辨力的缺點，具有局部自適應分析能力，在工程應用中得到了廣泛的應用。其局限性體現在以下兩點：首先，小波分析以傅里葉變換為基礎，雖然提供了一個尺度可變的時間-頻率窗，但是過度調整窗口的尺度會產生信息泄露的問題。其次，由于已經根據原始信號的特征人為選定了與之相匹配的基函數，且在后續的分析過程中無法進行更改，不同的小波基函數的分析結果往往是不同的［6-7］。

為滿足精確描述頻率隨時間變化的需要，美籍華人科學家Norden E.Huang提出了一種新的信號處理方法：經驗模態分解（empirical mode decomposition，EMD）［8］。該方法對于非平穩線性信號的分析具有非常好的效果，其核心思想是對非平穩線性信號進行平穩化處理，逐級分解信號中所包含的不同尺度下的波動或者變化趨勢，將分解得到的一系列不同特征尺度的數據序列稱之為本征模態函數（Intrinsic Mode Function，IMF）。對每一個IMF進行希爾伯特變換（Hilbert Transform，HT），從而得到不同特征尺度下的瞬時振幅和瞬時頻率。該理論的創新之處在于定義了瞬時頻率的概念，引進了EMD和IMF的概念。通過對IMF分量進行希爾伯特變換得到時頻域的Hilbert能量分布頻譜圖，建立了譜分析與時域的對應關系，能夠精確分析原始信號時頻域的分布規律［9-10］。

汪偉等［11］針對跑步機用戶進行行走速度估計過程中出現的加速度漂移問題，根據跑步機用戶行走速度數據的特征，利用EMD方法在不同頻率段內將擺動腳速度進行分解。通過選擇與用戶步頻相近的頻段，摒棄無關頻段，并對本證模態函數進行重構，從而解決了加速度飄逸的問題。采用集成經驗模態分解法對速度積分后得到的擺動腳位移進行分解與重構，有效解決了積分累積誤差引起的模態混疊現象。徐可、李俊等［12-13］充分利用EMD方法在自適應分析方面的優勢，采用將EMD分解后的本征模態函數所包含的能量特征和支持向量機（SVM）遺傳算法相結合的方法對非平穩軸承振動信號處理進行分析處理，分析的結果表明該方法在識別和定位滾動軸承故障方面具有良好的效果。汪朝海等［14］利用EMD方法在主成分分析、信息提取和自適應分析方面的優勢，對數據分析EMD分解。將IMF函數和殘差函數組成的特征向量作為滾動軸承故障識別的主要特征。對振動數據的分析結果顯示，與支持向量機、BP神經網絡等算法相比，該方法在故障識別率、運行時間和操作便捷性方面更具有優勢。楊建華等［15］對降噪后的滾動軸承故障數據進行EMD分解。通過對比分析軸承故障仿真信號和滾動軸承實驗信號的之間的特性信息，能夠對包含強噪聲的滾動軸承數據進行故障分析和定位。晏紅艷等［16］針對傳統相干噪聲壓縮去噪法對含噪信號進行去噪過程中容易消除有效噪聲的問題，結合多通道聯合時頻分析法，將含噪地震信號分解為多個頻率段下的分量。利用EMD方法在特征信號提取和自適應分解方面的優勢，在消除了含噪地震信號中噪聲的同時也有效保存了原始信號中的有用信息。

本文針對發動機試驗脈動壓力數據的非平穩非線性特征，結合經驗模態分解在自適應特征分析、瞬時頻率的精確定位能力、局部瞬時表達能力以及提取信號分量等方面的優勢，對某型號發動機試驗的脈動壓力數據進行分析。

1 脈動壓力測量系統的基本組成

工程應用中，脈動壓力測量系統具有頻率響應快，采集系統采樣率高的特征，能夠精確測量推力室高溫、高壓和劇烈振動環境下的壓力變化信息。測量系統一般由數采和分析系統、供電系統、脈動壓力傳感器、信號轉換裝置、循環冷卻系統和信號轉換裝置等幾個部分組成［8］，測量系統組成原理如圖1所示。

圖1 脈動壓力測量系統組成原理圖Fig.1 Schematic diagram of fluctuating pressure measurement system composition

脈動壓力傳感器的功能主要是完成脈動壓力測點的測量，實現測點壓力值到電信號的轉換。傳感器必須具有耐強振動、高溫和高壓的工作特性。脈動壓力傳感器的安裝以及測壓導管的長度等都有一定的要求［8］。

數采和分析系統主要完成數據的采集和分析，對其采樣率和頻響也有相應的要求，采集系統采樣率一般大于等于10 kHz。

循環冷卻系統主要功能是在傳感器工作過程中通過水冷或者氣冷的方式降低傳感器內部環境的溫度，保障傳感器正常工作溫度環境，延長工作壽命，提高系統可靠性。

2 經驗模態分解的基本理論與優勢

2.1 經驗模態分解的基本理論

對于一個實信號x(t)，可以構建一個復信號z(t)，定義為［10］

式中：a(t)，θ(t)分別為實信號x(t)的瞬時幅值和瞬時相位。

瞬時頻率等于瞬時相位的導數，即

令EMD分解得到的IMF分量為hi(t)，對其進行希爾伯特變換為

根據公式構造解析信號

則幅值函數為

式中：ai(t)為解析信號的瞬時幅值；φi(t)為解析信號的瞬時相位。

瞬時頻率為

希爾伯特譜為

希爾伯特邊際譜為

式中：T為所分析原始信號的總時長。從式（9）可知，希爾伯特譜能夠在時、頻域內精確描述信號的幅值的變化規律。

工程應用中對EMD分解后的征模態函數（IMF）分量進行HT變換時，IMF分量在數學上需要滿足時間歷程中極值點個數和過零點個數相等或者相差1個和局部極大值、極小值相加的均值為0的要求［17］。雖然IMF函數能夠合理地反映信號內部在每一個周期上的波動性，但是并非在所有情況下滿足上述要求的IMF函數都能夠計算出瞬時頻率。其主要的原因在于工程應用中原始信號往往較為復雜。真實信號中包含大量不相關的噪聲分量使得原始信號中包含多個振動模態分量。針對該問題，Huang和相關研究人員對IMF函數及相應原始信號的特征進行了研究和修訂，提出符合對信號進行EMD分解及對IMF函數進行HT變化的原始信號需要滿足以下假定條件［17］：①原始信號至少包含一個極大值點和一個極小值點；②相鄰兩個極值點之間的時間間隔為特征時間尺度；③對于信號中無極值點，僅有拐點的情況，采用對信號進行微分的方法獲得極值點。

EMD分解的步驟主要分為：

1）定位x(t)中所有極大值點和極小值點。采用三次樣條插值法繪制x(t)的上下包絡線（極大值點形成上包絡線，極小值點形成的下包絡線）。計算上、下包絡線的均值m1。x(t)與m1的差記為h1，即

將h1作為新信號x(t)，重復式（11）的步驟。直至所得到的分量hi滿足IMF函數的基本條件［17］時，將滿足條件的第1個一階IMF分量稱之為C1，C1包含信號中的最高頻成分。

2）將一階IMF分量C1從x(t)中分離出來，令剩余的原始信號為r1，令

將r1作為新信號，重復步驟1），直到第n階的殘余信號成為單調函數。

3）數學上x(t)可視為n個IMF分量和一個殘余分量的和，即

式中：rn(t)為信號中的殘量。各IMF分量Cj(t)則代表了信號中所包含的所有不同高、低頻率的成分。

在實際操作過程中需利用計算機編程完成信號的EMD分解過程。分解流程一般如圖2所示。

2.2 經驗模態分解的優勢

2.2.1 自適應時頻分析

從2.1章節的基本理論可知，EMD理論建立了瞬時頻率的物理意義。通過將原始信號分解為多個IMF分量，對各IMF分量進行HT變換得到原始信號在時、頻域的希爾伯特譜。瞬時特征量精確定位了希爾伯特譜在時頻域的特征。相較傳統的時頻分析方法，EMD方法的優點在于定義了一組根據原始信號特征而制定的具有自適應分解特性的基函數，從而使該方法在實際應用中具有良好的自適應分析優勢。

2.2.2 信號局部瞬時特征表征

根據2.1章節的理論可知，通過給出瞬時頻率的定義，即對EMD分解過程中所得到的多個IMF分量的相位進行求導，能夠從物理意義層面對瞬時頻率在局部瞬時特性方面的優勢進行清晰、明確地說明。相較傳統傅里葉變換在頻率混疊、虛假諧波方面的缺陷，EMD分解中多個IMF分量從高階到低階的排序也有效解決了上述問題。多個IMF分量瞬時頻率相位求導的方法則體現了該方法具有對原始信號的局部特征進行精確描述的能力。

2.2.3 主成分分析

在工程應用中，復雜的非線性非平穩隨機信號往往包含了大量真實信號和各種不相關的噪聲分量。對復雜信號進行分析的重要目的是將復雜原始信號中有價值的真實信號提取出來。從2.1章節的介紹可知，EMD方法能夠將原始信號中的頻率分量按照從高到低的順序分離提取出來以滿足信號分析的需要。數據分解的結果能夠較好地體現各頻率分量在原始信號中的比重。因此，該方法也是一種能夠對復雜信號進行主成分分析、提取的重要工具。

3 實例分析

3.1 脈動壓力數據特征

基于脈動壓力測量系統高頻響、高采樣率的特性，脈動壓力數據具有采樣率高、數據容量大的顯著特征［18］。某型號發動機試車中，按照圖1組建脈動壓力參數測量系統。傳感器采用Kulite-CT-375水冷型脈動壓力傳感器。傳感器的性能參數和計量機構檢定結論分別如表1、表2所示。

表1 CT-375型傳感器性能參數Tab.1 Performance parameters of CT-375 sensor

表2 CT-375型傳感器檢定結果Tab.2 Verification results of CT-375sensor

推力室脈動壓力參數名為Pis1，脈動參數采集系統的采樣率為10 kHz。選取試車中穩定段0.2 s（2000個樣本點）的數據進行分析。試驗全程脈動壓力數據如圖3所示。相鄰測量部位穩態壓力測量數據如圖4所示。測點名：Pcio1，穩態參數采集系統的采樣率為1 kHz。由于脈動壓力測點和穩態壓力測點分屬不同的測量系統，所用數據采集設備不同，故二者在坐標軸上顯示的時間有所區別，但發動機總工作時間是相同的。

從圖3和圖4可知，脈動壓力測點數據和穩態壓力測點數據在波形變化趨勢上基本一致，印證了測量數據的合理性。相較高采樣率、高精度脈動壓力測量系統，穩態壓力測量系統的測量數據在波形上較為光滑。而脈動壓力測量系統憑借高采樣率、高頻響和高測量精度的優勢，能夠非常精確地捕捉到推力室內壓力的波動變化情況。

圖3 某型號發動機試驗脈動壓力全程數據圖Fig.3 Full range data diagram of fluctuating pressure of a certain type of engine test

圖4 某型號發動機試驗穩態壓力全程數據圖Fig.4 Full range data diagram of steady-state pressure of a certain type of engine test

脈動壓力數據穩定段0.2 s和0.02 s局部數據分別如圖5和圖6所示。圖5中，在發動機試車穩定段26.10～26.30 s內，推力室壓力在8.614～8.911 MPa之間震蕩波動，變化范圍達到0.297 MPa。圖6清晰地顯示了推力室內壓力突變的情況，具有高采樣率特征的脈動壓力測量數據在精確描述推力室內脈動壓力在時域內波動情況的同時，也為工程技術人員采用多種時頻分析方法對其進行分析奠定了重要的基礎。

圖5 發動機試驗脈動壓力局部數據圖1Fig.5 Local data diagram 1 of engine test fluctuating pressure

圖6 發動機試驗脈動壓力局部數據圖2Fig.6 Local data diagram 2 of engine test fluctuating pressure

3.2 數據分析

根據2.1章節的理論基礎和方法，對穩定段26.10～26.30 s內的脈動壓力數據進行經驗模態分解，分解結果如圖7所示。

圖7 EMD分解結果圖Fig.7 EMD decomposition result

從圖7中EMD分解結果可知，采用EMD分解方法能夠將脈動壓力數據中的頻率分量分離出來。IMF1～IMF4階分量分別顯示了脈動壓力數據所在時域內所包含的頻率成分以及幅值相對大小。其中IMF1階分量在原始信號中所占的比重最大，頻率分量最高。IMF2～IMF4階分量在原始信號中所占比例逐漸降低。IMF2，IMF3和IMF4分量的幅值波動趨勢要大于IMF1分量的幅值波動。R主要為系統中的一些其他不相關低頻分量，變化過程非常緩慢且無規律，對數據的分析基本無影響。

對IMF函數進行HT變換的希爾伯特譜如圖8所示。

圖8 IMF1～IMF4時間頻率曲線圖Fig.8 IMF1～IMF4 time frequency curve

由圖8可知第1層～第4層時間-頻率曲線非常清晰準確地顯示了全時間域內頻率的變化情況。在第1，2層中，3500～5000 Hz左右的中頻帶是數據的主要組成成分，1500 Hz以下的低頻帶所占的比重較少。第2，3和4層時間-頻率曲線頻率值變化情況顯示在該時域內0.03，0.048，0.155和0.188 s存在多個頻率值突變。結合推力室上下游相應的溫度、壓力測點、試驗時序和試驗工況等信息，對于分析判斷推力室內燃燒狀態具有非常重要的作用。

對分解得到的IMF分量進行重構后的時間-幅值特性曲線如圖9所示。脈動壓力數據的FFT頻譜分析圖如圖10所示。

結合圖5、圖8、圖9和圖10進行分析可得到以下幾點結論：

圖9 EMD重構時間-幅值圖Fig.9 EMD reconstruction time-amplitude diagram

圖10 FFT頻譜分析圖Fig.10 FFT spectrum analysis diagram

1）通過對比圖5與圖9可知，圖9中重構后的時間-幅值特性曲線相比原始數據具有更高的分辨力，能夠更加清晰地分辨時間域內的幅值波動情況。不僅可以清晰分辨脈動壓力周期性波動的情況，而且可以精準定位在0.12，0.18 s的幅值振蕩。說明在該試車的穩定段內，隨著推進劑燃燒對壓強振蕩響應，燃燒產生的能量持續注入燃燒室工作系統引起脈動壓力有規律地振蕩。通過對比同時段推力室對應的振動測點的振動數據，脈動壓力有規律的振蕩在振動數據中也得到了體現。此外，根據EMD分解的特點，針對渦輪泵、燃發器等振動劇烈的測點，還可參考文獻［8］，［18］中的方法，結合相應振動測點的數據進行分析。

2）圖8中第1層和第2層的時間頻率特性曲線與圖10中FFT分析結果相吻合。圖8中第3層～第4層時間-頻率特性曲線非常清晰地顯示了低頻分量的變化相對較緩慢的特性，也能非常清晰定位對應的變化時間點。

3）由圖10中FFT分析結果可知，全時域內幅值峰值多集中在低頻區（0～2000 Hz）和中頻區（3500～5500 Hz）。這與圖8中EMD分解時間-頻率曲線的結果相吻合。但FFT的不足之處是分析結果只能顯示全時域數據在頻域內的信息，無法定位某個幅值峰值所發生的時間信息，也無法分析信號中的各頻率成分。

分別采用Haar小波和Db9小波對脈動壓力數據進行7層分解的結果如圖11至圖14所示。

圖11 Haar小波5～7層分解細節系數圖Fig.11 Detail coefficient diagram of Haar wavelet 5～7 level decomposition

圖12 Haar小波1～4層分解細節系數圖Fig.12 Detail coefficient diagram of Haar wavelet 1～4 level decomposition

圖13 Db9小波5～7層分解細節系數圖Fig.13 Detail coefficient diagram of Db9 wavelet 5～7 level decomposition

圖14 Db9小波1～4層分解細節系數圖Fig.14 Detail coefficient diagram of Db9 wavelet 1～4 level decomposition

由小波分析圖11至圖14可知，小波分析的結果與所選小波基函數有關（文獻［18］介紹了多種小波基函數的特征）。由于Haar小波具有緊支撐、單個矩形波、矩形狀階梯變化等特征，分解得到的多層細節分量具有矩形特征，尤其是d4～d7層細節信息尤為明顯。而Db系小波具有對稱性好、頻域的局部化能力強和消失矩光滑性好等特點，分解得到的多層細節分量的光滑性較好。d1～d5層分解信號中均能定位原始信號中幅值產生突變的時間信息。d6～d7層的分解結果均因分解層數高，過濾掉了過多的有用信息而存在一定程度的失真。小波分析的結果可以結合相應的振動測點數據，進一步分析所測部件的動態特性［8］。小波基函數的選擇與分解層數的差異導致小波分析所得到的結果存在一定的差異，也給準確分析數據特征帶來了困難。

4 總結

本文介紹了型號發動機試驗中脈動壓力測量系統的組成和數據特征。分析了采用傳統傅里葉變換、小波分析等方法對非線性非平穩隨機數據進行時頻分析的缺點和不足。重點介紹了EMD的基本理論。從自適應時頻分析、局部瞬時定位分析和信號主成分特征分析和提取三個方面介紹了該理論的優點以及實際工程應用的詳細步驟。分別采用FFT方法、基于不同小波基函數和分解層數的小波分析方法和EMD分解方法對某型號發動機試驗脈動壓力數據進行分析。對比分析了FFT方法和小波分析方法對非平穩非線性數據的不足。脈動壓力數據的EMD分解結果顯示該方法能夠非常清晰地自適應分離原始信號中的高、低頻分量，顯示數據在時間-頻率域的多維特征。該方法可以與其他數據分析方法相結合，為產品型號設計和試驗人員分析非線性非平穩型數據提供重要的參考依據，具有重要的參考和推廣價值。