液體火箭發動機試驗脈動壓力數據分析方法研究

楊 懿 蹇 玥 賈志杰 王永鵬 郭亞男

(北京航天試驗技術研究所，北京 100074)

1 引 言

在液體火箭發動機試驗中，壓力是反映發動機和試車臺設備工作過程特性的重要參數[1]。根據所測壓力信號隨時間的變化類型，壓力參數可分為穩態壓力和脈動壓力。穩態壓力一般用于考察和研究發動機、試驗臺工藝系統的能力。脈動壓力一般用于測量燃燒室、推力室內高頻、高壓、高溫壓力，能彌補穩態壓力測量方法在高、低溫和水擊壓力等特殊環境下測量能力的不足。

針對脈動壓力數據的特點，研究行之有效的數據分析方法，挖掘測量數據背后的有價值信息是研究發動機動態性能和燃燒狀態的重要手段。

脈動壓力參數數據采樣率高，系統響應頻率快，僅僅分析數據在時域內的信息是遠遠不夠的，還需要深度挖掘數據的其他信息。商霖[2]等通過研究燃燒室內脈動壓力能量分布的頻率特性與隨機振動的關系，研究發動機的工作原理。劉英元[3]等以數據頻率為研究對象，通過FFT方法和小波變換方法提取發動機狀態特征的方法。費繼友[4]等將小波分析與奇異值檢測理論相結合，應用小波奇異檢測和誤差濾波對瞬態數據進行處理，通過試車試驗驗證該方法有效。

本文采用FFT方法和小波變換奇異點檢測方法對發動機脈動壓力數據進行分析，提取特征信息，研究發動機的不穩定燃燒狀態。

2 脈動壓力測量系統組成和數據特征

2.1 系統組成

脈動壓力測量系統一般由脈動壓力傳感器、傳感器循環冷卻系統、供電電源、信號轉換裝置、數據采集和分析系統組成，系統組成圖如圖1所示。

圖1 脈動壓力測量系統組成圖Fig.1 Fluctuating pressure measurement system diagram

2.2 數據特征

脈動壓力數據需要精確體現燃燒室、推力室等高溫、高壓、劇烈振動環境下的壓力變化信息。測量系統一般具備以下特征：測量系統的頻率響應快，通常在10kHz以上；傳感器耐高溫、高壓和劇烈振動，頻率特性至少應有10k～20kHz的平坦段；采集系統采樣率高，一般具備10kHz/s以上的高采樣速率。因此，脈動壓力數據具有高頻、采樣點多和數據容量大的顯著特征。

3 脈動壓力數據分析方法

本文采用兩種方法分析某型號發動機地面試車的脈動壓力數據：(1)FFT分析方法。結合相應振動測點的數據，分析脈動壓力數據的幅頻特性，為判斷是否發生不穩定燃燒提供依據。(2)采用小波分析法對數據進行去噪、奇異點檢測，研究發動機動態性能和燃燒狀態。

4 數據分析方法的理論基礎

4.1 FFT理論基礎

在數字信號分析領域，FFT是分析非周期信號頻域信息的重要工具。該方法利用傅里葉積分能夠將非周期的函數在無限區間上分解的性質，對連續時域信號進行抽樣和截斷，實現數字信號計算機分析處理。

窗函數的選擇是FFT的重要步驟，直接影響頻譜分析的結果。常見的窗函數有Hamming窗、Hanning窗、矩形窗和Bartlett窗。Hanning窗具有主瓣寬、旁瓣峰值低、衰減速度快和能有效減少頻譜泄漏的優點，本文采用Hanning窗進行數據分析。

4.2 小波分析的理論基礎

小波分析法通過小波基函數的伸縮和平移實現信號時頻兩域自適應分析，憑借其強大的挖掘信號細節信息的能力，近年來在工程應用領域得到了廣泛的應用[5]。

信號的奇異性分析是提取信號特征信息的重要手段。在工程應用中，信號的奇異點包含系統運行過程中很多的重要信息，對奇異點的分析能夠挖掘系統運行的關鍵信息。

采用小波分析法將突變奇異信號進行多尺度分解，根據Lipschitz指數和模極大值在不同尺度上的關系分析信號的突變特征并定位奇異值的位置和大小。通過對信號奇異點突變時間、突變類型和變化的幅值等參數的考查，對發動機的燃燒狀態進行分析[3]。

4.2.1小波變換與分解重構

設φ(t)為一個基本小波，L2(R)表示平方可積的實數空間，φ(t)∈L2(R)，若其傅里葉變換函數Ψ(ω)滿足條件

(1)

式中：R——實域；ω——屬于實域R的變量。

設伸縮系數為a,a>0且連續可變化數值。設平移系數為b，b∈R且連續可變化數值。將基本小波函數φ(t)進行伸縮和平移，令伸縮和平移后的函數φa,b(t)為

(2)

式中：φa,b(t)——依賴于a,b變化的連續小波基函數。

對于將脈動壓力數據視為非周期信號，用函數x(t)來表示，其連續小波變換的表達式為

WTx(a,b)=〈x(t),φa,b(t)〉=

(3)

式中：WTx(a,b)——表示對函數x(t)進行以a,b為系數的連續小波變換。

調整a和b的大小，不僅可以將信號分解到各個指定的頻帶，還可以保留各分量的時域信息。為了便于計算機進行計算分析，需要將連續小波函數、a和b離散化[5]。離散化的方法通常把伸縮系數a和平移系數b取冪級數。在實際的應用中，對a和b進行二進離散時一般取a0=2,b0=1，即a=2j，b=k2j，j,k為整數。

小波變換方法對信號進行分析處理時一般采用Mallat算法將信號進行多重分解與重構。算法的基本原理是：利用小波高通、低通濾波器對函數x(t)在分辨率2j下進行分解，獲取低頻和高頻小波系數。小波重構則是小波分解的逆運算。采用相關算法將高、低頻小波系數進行重構[6]。

4.2.2小波突變點檢測的基本原理

在小波分析方法中，通常用Lipschitz指數α描述信號的奇異性。隨著小波變換尺度的增加，小波變換的模極大值在信號奇異點的衰減速度取決于信號在突變點的指數α。

采用小波變換來描述連續信號f(t)在t0點的奇異指數。設小波基具有n(n為正整數)階消失矩，且n階可微、具備緊支撐特征。令f(t)∈L2(R)。若在t0點的鄰域和所有尺度上，存在常數K滿足

|WTf(s,t)|≤K(sα+|t-t0|α)

(4)

式中：|WTf(s,t)|——對函數f(s,t)進行以s、t為系數的連續小波變換；K——常數；s、t——大于0的常數；t0——數據中的奇異點；α——f(t)在點t0處的Lipschitz指數為α,0<α≤n。

從公式(4)可以看出，信號奇異點分布在模極值線上，其Lipschitz指數α不等于1，突變信號呈現奇異性且指數α>0。當t0為信號f(t)的局部奇異點，t取離散值，則在該點的小波變換取得模極大值。在離散二進小波變換中，公式(4)變換為

(5)

式中：j——二進尺度參數；A——常數。

從公式(5)可以看出，如果信號在t0處的奇異指數大于零，則隨著尺度j的增加，小波變換模極大值也增加。在應用中，通過小波變換后的模極大值在不同分解尺度上的衰減速度衡量信號的局部奇異性。當小波函數可以看作平滑函數的一階、二階導數時，信號小波變換模的局部極值點和過零點分別對應于信號的突變點[7]。因此可建立信號奇異點與模局部極大值之間的關系。但是模極大值與奇異點之間并非一一對應的關系，只有在多個尺度上檢測均為極值點的位置才是真正的突變點所在位置，繼而進一步確定突變點的類型。在進行奇異點定位之前，需要對數據進行去噪處理，排除系統中噪聲的干擾[4]。

小波奇異點定位分析的步驟為：對信號進行小波離散化處理；對信號進行小波分解；對利用小波函數對信號進行去噪；對信號進行奇異點定位、分析；對分解信號進行小波重構。

4.2.3小波基函數的選擇

在工程應用中，常用的小波基函數有：Haar函數、DbN函數、SymN函數、CoifN函數和Dmey函數等幾種。各種基函數都有自身的特點，使用的基函數不同，得到的處理結果也不一樣。在實際的運用過程中需要根據數據和小波基函數自身的特點，選擇合適的小波基函數進行數據處理[3]。

小波基函數的選取主要遵循三點原則：(1)支撐長度是指尺度函數、小波函數及其傅里葉變換的收斂速度；緊支撐基函數的收斂速度快，在描述原始信號時能減少計算量，提高計算的精度。而非緊支撐基函數的收斂速度慢，在描述原始信號時的計算效果和逼近效果要差一些。(2)在選擇基函數時，一般應選擇具備緊支撐特征的基函數。(3)對稱性能夠有效避免信號處理中的移相問題，正則性和消失矩呈正相關關系，主要影響小波系數重構的穩定性。各種小波基函數的特點如表1所示。

表1 各小波基函數特點Tab.1 Characteristicsofwaveletbasisfunctions名稱Haar函數DbN函數SymN函數CoifN函數Dmey函數特點正交、緊支撐、單個矩形波、矩形狀階梯變化、計算簡單正則性好、對稱性好、頻域的局部化能力強、頻帶劃分效果好、緊支撐、消失矩光滑性好DbN函數的一種改進型、具備DbN函數的優點、對稱性比Db函數好對稱性比Db函數好不具備緊支撐性、收斂速度快

5 發動機試驗脈動壓力數據分析

5.1 FFT分析

采用FFT方法分析試驗數據時，一般需要針對發動機不同的工況，如啟動段、穩定段和關機段等工況分段分析。以某次發動機試驗穩定段工況的燃燒室脈動壓力數據為例。選取穩定段0.2s(2000個樣本點)的數據進行FFT頻譜分析，數據采樣率為10kHz/s，測點參數名為Pfs。試驗全程脈動壓力數據如圖2所示，穩定段0.2s數據如圖3所示。為了能夠清晰觀察數據樣本在各頻率段內的幅值變化情況，分別列出(0～5000)Hz和(0～11000)Hz頻率段的頻譜分析圖，如圖4、圖5所示。

圖2 Pfs全程脈動壓力數據圖Fig.2 Pfs pressure conduit system diagram

圖3 Pfs測點0.2s脈動壓力數據圖Fig.3 0.2s fluctuating pressure data of Pfs measuring point

圖4 脈動壓力Pfs(0～11000)Hz頻譜分析圖Fig.4 Spectrum analysis diagram of fluctuating pressure Pfs(0～11000)Hz

圖5 脈動壓力Pfs(0～5000)Hz頻譜分析圖Fig.5 Spectrum analysis diagram of fluctuating pressure Pfs(0～5000)Hz

圖3中穩定段0.2s脈動數據可以看出，燃燒室內脈動壓力幅值大致在(5.0～5.4)MPa范圍內波動且呈現周期性的壓強振蕩特征。從圖4和圖5的FFT數據分析頻譜圖可以看出幅值峰值主要集中在(0～5000)Hz的頻率段內，(0～1000)Hz的中、低頻率段的幅值波動不明顯，(1000～5000)Hz高頻段出現多個幅值峰值，且振幅的幅值呈遞增趨勢。說明隨著推進劑燃燒對壓強振蕩響應，燃燒產生的能量持續注入燃燒室工作系統。由此可以推斷發生高頻聲振的可能性較大。

為研究穩定段的不穩定燃燒聲振的類型，采用FFT分析方法對燃燒室軸、徑、切向三個振動測點數據進行分析。徑向振動測點ab1-2的頻譜分析結果如圖6所示?？梢钥闯鲈?1000～5000)Hz頻率段，徑向振動測點ab1-2存在多個幅值峰值。

圖6 ab1-2振動測點頻譜分析圖Fig.6 Spectrum analysis diagram of ab1-2 vibration measuring point

軸、徑、切向三個振動測點FFT頻譜分析幅值峰值頻率點和脈動壓力數據幅值峰值頻率點的對比分析情況如表2所示。

表2 燃燒室振動測點與脈動壓力幅頻數據表Tab.2 Vibrationmeasurementpointsandamplitudefrequ-encydatatableofpulsatingpressureincombustionchamberPfs幅值峰值對應的頻率(Hz)ab1-1(軸向)ab1-2(徑向)ab1-3(切向)904╱╱╱1611╱√╱2243╱╱╱2884╱√╱3685╱╱╱4822√√╱9643╱√╱ 注:√表示在振動數據FFT分析的幅值峰值的頻率值與Pfs吻合;╱表示在振動數據FFT分析的幅值峰值的頻率值與Pfs不吻合。

通過表2的分析數據可以看出，脈動壓力測點Pfs的幅頻峰值特性與振動徑向測點ab1-2的幅頻峰值特性在1611Hz、2884Hz、4822Hz和9643Hz四個頻率點，即發動機工作頻率的2、6、12倍頻存在耦合關系，而與軸向和切向的耦合點較少，由此可以推斷燃燒室發生高頻橫向聲振。通過上述分析得出結論：脈動壓力數據可以清晰地反映燃燒室內的壓強振蕩；燃燒室內產生不穩定燃燒時，通過對脈動壓力數據進行FFT分析可以獲取其幅頻特性。通過分析脈動壓力和振動測點的頻譜特征，可以為判斷發動機燃燒不穩定性聲振的類型提供強有力的參考依據。

5.2 小波奇異點分析

某次發動機試驗燃燒室脈動壓力測點Pios1全程試驗數據如圖7所示。從圖7可以看出在發動機試車啟動段和關機段的脈動壓力波動異常。

圖7 Pios1脈動壓力全程數據圖Fig.7 Pios1 fluctuating pressure data chart

分別選取啟動段(3.18～3.58)s，樣本點數4000和關機段(12.29～12.39)s，樣本點數1000的數據進行奇異點特征分析。數據采集系統采樣率均為10kHz/s。

5.2.1啟動段數據分析

啟動段(3.18～3.58)s的數據如圖8所示。在奇異點檢測中，選擇小波的正則性非常重要。根據對各類小波的特征分析，分別采用db1和db9小波對啟動段的數據進行相似系數計算、FFT分析和7層分解重構。由于篇幅關系，僅列出數據處理效果較好的db9小波分析結果。db9小波分解7層近似系數結果如圖9所示。(0～10000)Hz和(0～250)Hz的FFT分析局部放大效果圖分別如圖10、11所示。采用db9小波7層分解細節系數結果如圖12所示。

圖8 Pios1啟動段脈動壓力局部數據圖Fig.8 Local data diagram of pulsating pressure in Pios1startup section

圖9 小波分析相似系數圖Fig.9 Similar coefficient graph of wavelet analysis

圖10 FFT幅頻特征圖(0～10000)HzFig.10 FFT amplitude frequency characteristic chart(0～10000)Hz

圖11 FFT幅頻特征圖(0～250)HzFig.11 FFT amplitude frequency characteristic chart(0～250)Hz

圖12 db9小波7層分解細節系數圖Fig.12 db9 wavelet 7-level decomposition detail coefficient diagram

從圖8中可以看出，啟動段脈動壓力存在多個峰值波動，且平均壓強呈遞增的趨勢。從圖9中可以看出，7層分解重構后的相似系數圖在波峰個數以及波峰對應數據點與圖8的原始數據基本一致，表明采用db9小波對數據進行處理很好吻合了數據的變化趨勢。從圖10和圖11的FFT分析結果可以看出，FFT方法能夠分析數據在頻域的特征，由于沒有時間分辨力，無法提取壓力幅值突變點的空間位置和特征信息。圖12列出了采用db9小波對數據進行7層分解后的細節系數。圖中可以看出，d1～d7的7層分解信號中均準確定位了原始信號中幅值突變奇異點的位置。7層分解信號中的噪聲隨著分解尺度的增加迅速減少，第6、7層分解重構后的曲線存在一定程度的失真，主要因為分解層數越高，過濾掉了更多有用的信息。第4、5層分解重構的效果最好，能夠清晰地觀察到信號中的不連續點和幅值突變點。d1～d3的高頻細節系數顯示信號中高頻部分是原始信號中突變的主要成分。通過對該時間段燃燒室軸向、徑向和切向振動測點的進行FFT頻譜分析，發現a3切向振動測點在高頻段(2000～8000)Hz能量較集中，高頻段振幅遠高于低頻段振幅。此外，a3切向振動測點高頻段的振幅均大于同頻率段軸向和徑向的振幅，且a3測點在工作倍頻之外還出現了多個振動峰值。燃燒室a3切向振動測點頻譜圖如圖13所示。

圖13 燃燒室關機段a3切向振動測點頻譜分析圖Fig.13 Spectrum analysis of a3 tangential vibration measuring point in shutdown section of combustion chamber

通過上述分析可以判斷在發動機試車的啟動段，燃燒室存在不穩定燃燒，在燃燒室切向產生高頻聲振，導致在啟動段脈動壓力的高頻段存在多個峰值奇異點。

5.2.2關機段數據分析

關機段原始數據如圖14所示。采用模極大值法，選取db9小波對原始信號進行7層分解重構，分解尺度為1～32。小波變換系數如圖15所示。

圖14 Pios1關機段脈動壓力局部數據圖Fig.14 Local data diagram of pulsating pressure in Pios1

圖15 小波變換系數與分解層數結果圖Fig.15 Results of wavelet transform coefficients and decomposition levels

從圖15可以清楚看出，在第50、876、918、984樣本點，小波系數呈現一個倒錐形的區域，且在2、4、8、16、32層分解的尺度上，小波變換系數均呈現一致間斷的特征(均為白色區域)。由此可以判斷在第50、876、918、984樣本點，信號均為奇異信號。該結論與圖14中原始數據的波動非常一致。從圖14的原始數據中可以觀察到在上述奇異點，脈動壓力的幅值存在(7.88～10.49)MPa幅度的異常波動。

6 結束語

根據脈動壓力數據的特點，本文給出了基于FFT和小波變換奇異點檢測的液體火箭發動機試車脈動壓力數據分析方法。通過對脈動壓力數據和相應振動測點的FFT頻譜分析，能夠在頻域內有效判斷發動機燃燒室的工作狀態和聲振類型。但是由于傅里葉變換在時域內的局限性，無法對信號中的突變點進行定位。采用基于小波變換模極大值和Lipschitz指數的小波奇異點檢測方法對發動機啟動段和關機段的異常數據在多個尺度進行分解分析，不僅能有效定位信號中的突變點位置，還能獲取信號在頻域內的信息。結合相應振動測點的頻譜數據分析結果，能夠為判斷發動機燃燒狀態和聲振類型提供有力的依據。上述兩種方法在其他型號發動機和組合件試驗的脈動壓力數據分析中也具有重要的推廣價值。