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

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

1 引 言

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

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

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

本文采用FFT方法和小波變換奇異點檢測方法對發(fā)動機(jī)脈動壓力數(shù)據(jù)進(jìn)行分析，提取特征信息，研究發(fā)動機(jī)的不穩(wěn)定燃燒狀態(tài)。

2 脈動壓力測量系統(tǒng)組成和數(shù)據(jù)特征

2.1 系統(tǒng)組成

脈動壓力測量系統(tǒng)一般由脈動壓力傳感器、傳感器循環(huán)冷卻系統(tǒng)、供電電源、信號轉(zhuǎn)換裝置、數(shù)據(jù)采集和分析系統(tǒng)組成，系統(tǒng)組成圖如圖1所示。

圖1 脈動壓力測量系統(tǒng)組成圖Fig.1 Fluctuating pressure measurement system diagram

2.2 數(shù)據(jù)特征

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

3 脈動壓力數(shù)據(jù)分析方法

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

4 數(shù)據(jù)分析方法的理論基礎(chǔ)

4.1 FFT理論基礎(chǔ)

在數(shù)字信號分析領(lǐng)域，F(xiàn)FT是分析非周期信號頻域信息的重要工具。該方法利用傅里葉積分能夠?qū)⒎侵芷诘暮瘮?shù)在無限區(qū)間上分解的性質(zhì)，對連續(xù)時域信號進(jìn)行抽樣和截斷，實現(xiàn)數(shù)字信號計算機(jī)分析處理。

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

4.2 小波分析的理論基礎(chǔ)

小波分析法通過小波基函數(shù)的伸縮和平移實現(xiàn)信號時頻兩域自適應(yīng)分析，憑借其強(qiáng)大的挖掘信號細(xì)節(jié)信息的能力，近年來在工程應(yīng)用領(lǐng)域得到了廣泛的應(yīng)用[5]。

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

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

4.2.1小波變換與分解重構(gòu)

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

(1)

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

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

(2)

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

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

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

(3)

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

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

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

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

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

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

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

(4)

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

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

(5)

式中：j——二進(jìn)尺度參數(shù)；A——常數(shù)。

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

小波奇異點定位分析的步驟為：對信號進(jìn)行小波離散化處理；對信號進(jìn)行小波分解；對利用小波函數(shù)對信號進(jìn)行去噪；對信號進(jìn)行奇異點定位、分析；對分解信號進(jìn)行小波重構(gòu)。

4.2.3小波基函數(shù)的選擇

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

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

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

5 發(fā)動機(jī)試驗脈動壓力數(shù)據(jù)分析

5.1 FFT分析

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

圖2 Pfs全程脈動壓力數(shù)據(jù)圖Fig.2 Pfs pressure conduit system diagram

圖3 Pfs測點0.2s脈動壓力數(shù)據(jù)圖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中穩(wěn)定段0.2s脈動數(shù)據(jù)可以看出，燃燒室內(nèi)脈動壓力幅值大致在(5.0～5.4)MPa范圍內(nèi)波動且呈現(xiàn)周期性的壓強(qiáng)振蕩特征。從圖4和圖5的FFT數(shù)據(jù)分析頻譜圖可以看出幅值峰值主要集中在(0～5000)Hz的頻率段內(nèi)，(0～1000)Hz的中、低頻率段的幅值波動不明顯，(1000～5000)Hz高頻段出現(xiàn)多個幅值峰值，且振幅的幅值呈遞增趨勢。說明隨著推進(jìn)劑燃燒對壓強(qiáng)振蕩響應(yīng)，燃燒產(chǎn)生的能量持續(xù)注入燃燒室工作系統(tǒng)。由此可以推斷發(fā)生高頻聲振的可能性較大。

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

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

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

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

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

5.2 小波奇異點分析

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

圖7 Pios1脈動壓力全程數(shù)據(jù)圖Fig.7 Pios1 fluctuating pressure data chart

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

5.2.1啟動段數(shù)據(jù)分析

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

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

圖9 小波分析相似系數(shù)圖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層分解細(xì)節(jié)系數(shù)圖Fig.12 db9 wavelet 7-level decomposition detail coefficient diagram

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

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

通過上述分析可以判斷在發(fā)動機(jī)試車的啟動段，燃燒室存在不穩(wěn)定燃燒，在燃燒室切向產(chǎn)生高頻聲振，導(dǎo)致在啟動段脈動壓力的高頻段存在多個峰值奇異點。

5.2.2關(guān)機(jī)段數(shù)據(jù)分析

關(guān)機(jī)段原始數(shù)據(jù)如圖14所示。采用模極大值法，選取db9小波對原始信號進(jìn)行7層分解重構(gòu)，分解尺度為1～32。小波變換系數(shù)如圖15所示。

圖14 Pios1關(guān)機(jī)段脈動壓力局部數(shù)據(jù)圖Fig.14 Local data diagram of pulsating pressure in Pios1

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

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

6 結(jié)束語

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