甘 林,張翔宇,金秉寧*,魏少娟,劉佩進
(1.西北工業(yè)大學(xué) 燃燒、熱結(jié)構(gòu)與內(nèi)流場重點實驗室,西安 710072;2.中國航天科技集團有限公司四院四十一所,西安 710025;3.空軍工程大學(xué) 裝備管理與無人機工程學(xué)院,西安 710000)
非線性燃燒不穩(wěn)定問題一直是固體火箭發(fā)動機研制過程中常遇到的棘手問題之一[1-5],其典型特征是在發(fā)動機燃燒中產(chǎn)生非線性壓強振蕩現(xiàn)象[6-10]。通常,這種非線性壓強振蕩中包含一個或多個聲振蕩模態(tài),且各階模態(tài)受增益阻尼因素的影響,能量從低階模態(tài)向高階模態(tài)傳遞,激發(fā)出高階模態(tài),使得不穩(wěn)定問題愈加嚴重[11-15]。然而,目前仍存在對能量傳遞演化機理認識不清的問題,無法對非線性燃燒不穩(wěn)定的預(yù)估和抑制提供理論指導(dǎo)。因此,需要開展能量傳遞演化機理的研究,而進行這項研究的前提則是從非線性振蕩數(shù)據(jù)中準確的得到各階模態(tài)的關(guān)鍵參數(shù),如振蕩頻率和幅值等。
對于固體火箭發(fā)動機非線性不穩(wěn)定振蕩數(shù)據(jù),科研人員多采用商業(yè)軟件平臺進行分析,其特點是功能強大,操作簡單,經(jīng)過長期探索,可以形成一套數(shù)據(jù)分析流程。但在使用過程中發(fā)現(xiàn),商業(yè)軟件中很多模塊方法和參數(shù)不明確,使得數(shù)據(jù)處理結(jié)果出現(xiàn)失真;此外,受主觀經(jīng)驗判斷的影響,分析結(jié)果也存在一定的誤差。在處理類似信號如電力系統(tǒng)低頻振蕩信號、機械系統(tǒng)振動衰減信號等方面,更多的是采用時頻分析方法,如經(jīng)驗?zāi)B(tài)分解(EMD)[16]和變分模態(tài)分解(VMD)[17-18]等,并且取得了較好的效果。綜上,有必要借鑒處理類似信號的分析方法建立一種方便而準確的分析固體火箭發(fā)動機非線性燃燒不穩(wěn)定振蕩數(shù)據(jù)的方法,以獲得非線性不穩(wěn)定振蕩數(shù)據(jù)的關(guān)鍵參數(shù)。
固體火箭發(fā)動機非線性燃燒不穩(wěn)定振蕩數(shù)據(jù)是典型的非線性信號。在非線性信號處理方面,典型的時頻方法有經(jīng)驗?zāi)B(tài)分解(EMD)[19]、局部均值分解(LMD)[20]、本征時間尺度分解(ITD)[21]和變分模態(tài)分解(VMD)[22]等。其中,前三種基于遞歸迭代的思想,能夠?qū)π盘柲B(tài)進行分解。然而,EMD存在著一些問題,如缺乏嚴格的數(shù)學(xué)公式,易受模態(tài)混疊、端點效應(yīng)、過擬合和欠采樣的影響;LMD在抑制端點效應(yīng)和模態(tài)混疊方面優(yōu)于EMD,但其主要問題是計算復(fù)雜度高,仍容易產(chǎn)生模態(tài)混疊;ITD存在末端效應(yīng)和時間分辨率低的問題;而VMD方法基于求解一個完全非遞歸變分問題的思想,具有理論上嚴謹?shù)臄?shù)學(xué)原理和對噪聲及采樣的魯棒性等優(yōu)點[23]。但VMD方法依賴于經(jīng)驗,需要人工預(yù)先設(shè)定模態(tài)分解個數(shù)K、二次懲罰項α和帶寬τ,參數(shù)設(shè)置不當(dāng)可能導(dǎo)致模態(tài)混疊,影響分解結(jié)果的準確性[24]。
綜合考慮,本文引入了VMD[24]算法,并采用標準粒子群優(yōu)化算法(SPSO 2011,以下簡稱SPSO)[25]對VMD算法所需的預(yù)置參數(shù)進行優(yōu)化,建立了SPSO-VMD方法。將SPSO-VMD方法應(yīng)用于典型特征的燃燒不穩(wěn)定聲振蕩信號處理,準確地獲得了各階模態(tài)的頻率和幅值,驗證了該方法的可行性。并將其應(yīng)用在了真實發(fā)動機振蕩數(shù)據(jù)的處理上,準確地分離出了各階模態(tài)。
VMD方法將原始信號X(t)分解為K(模態(tài)分解個數(shù))個緊繞其中心頻率ω的本征模態(tài)函數(shù)。在該算法中,定義本征模態(tài)函數(shù)(IMF)為一個調(diào)幅-調(diào)頻信號,VMD分解過程實質(zhì)上是變分問題的構(gòu)造和求解過程,詳見文獻[24]。
SPSO算法是粒子群算法(Particle Swarm Optimization,PSO)[26]的一種改進。PSO算法作為一種群體智能優(yōu)化算法,具有原理簡單、易于實現(xiàn)和搜索速度快等優(yōu)勢,在PSO算法中,群體中的每個個體稱為一個粒子,代表優(yōu)化問題的一個潛在解。對于一個n維優(yōu)化問題,每個粒子被賦予一個n維速度向量和n維位置向量。在每一次迭代也即搜索的過程中,每個粒子根據(jù)自身的搜索經(jīng)驗和群體中其他粒子的搜索經(jīng)驗動態(tài)調(diào)整各自的運動狀態(tài),從而保證每個粒子逐步向問題的最優(yōu)解方向靠近[27-28]。但是,基本PSO算法存在兩個缺陷:(1)不能很好地動態(tài)調(diào)整粒子的全局和局部搜索能力;(2)粒子容易陷入搜索迭代停滯[29]。RATNAWEERA等[30]為克服PSO算法的第一個缺陷,提出了線性自適應(yīng)參數(shù)調(diào)整法則,即根據(jù)當(dāng)前迭代次數(shù)和最大迭代次數(shù),線性地調(diào)整參數(shù)的值;CLERC等[25]為克服基本PSO算法過早陷入迭代停滯的缺陷,提出了SPSO算法。
在SPSO算法中,首先根據(jù)每個粒子的當(dāng)前位置Xi(t)、粒子的個體最優(yōu)位置Pbesti(t)以及群體的全局最優(yōu)位置Gbest,為每個粒子定義一個超球面;然后,在定義的超球面中,為每個粒子隨機產(chǎn)生一個新的位置點;最后,將這個隨機產(chǎn)生的位置視為擾動添加到粒子的速度更新公式中,所以在迭代初期,SPSO算法可從最大程度上保證每個粒子的飛行速度不為0??梢员M量保證在迭代的初期每個粒子都能以非零速度在解空間中一直搜索,從而避免算法過早陷入迭代停滯。
本文按照線性自適應(yīng)參數(shù)調(diào)整法則的SPSO算法優(yōu)化VMD分解,建立了SPSO-VMD方法,其思路如下:影響VMD算法的三個需要提前輸入的參數(shù)為模態(tài)分解個數(shù)K、二次懲罰因子α和帶寬τ。在處理的過程中,模態(tài)分解個數(shù)K是通過對振蕩信號進行傅里葉變換并濾波從而確定的;二次懲罰因子α和帶寬τ是通過SPSO算法來優(yōu)化的。
SPSO算法優(yōu)化需要一個適應(yīng)度函數(shù)f(x),其定義如下:首先,對振蕩信號進行傅里葉變換并濾波,得到K個模態(tài)頻率及對應(yīng)的K個幅值A(chǔ)k;再使用VMD算法對振蕩信號進行分解得到K個模態(tài)分量,對這些模態(tài)分量進行傅里葉變換,得到K個模態(tài)分量對應(yīng)的K個幅值A(chǔ)k,VMD(k=1,2,…,K)。各模態(tài)幅值差值ΔA計算公式為
ΔAk=|Ak-Ak,VMD|
ΔA={ΔA1,ΔA2,…,ΔAk}
(1)
將Max(ΔA)定義為SPSO算法的適應(yīng)度函數(shù)f(x),也即SPSO優(yōu)化VMD算法是朝著值Max(ΔA)變小的方向去優(yōu)化的。
應(yīng)用SPSO算法優(yōu)化VMD算法流程如下:首先設(shè)定兩個參數(shù)α和τ優(yōu)化的具體范圍,粒子在這個范圍內(nèi)進行全局尋優(yōu),α和τ的值由粒子的位置來表征,每個粒子每次循環(huán)將自身的位置代入VMD算法對振蕩信號進行一次分解,并計算一次適應(yīng)度值,通過對比新粒子的適應(yīng)度值進行位置的更新,使每個粒子向著最優(yōu)解方向更新位置,通過不斷的迭代更新,最終找到使整個種群全局適應(yīng)度最優(yōu)的位置,也即所需優(yōu)化參數(shù)的最佳值。
SPSO-VMD方法流程圖如圖1所示。
固體火箭發(fā)動機不穩(wěn)定振蕩數(shù)據(jù)中包含振蕩增長段、極限環(huán)段以及振蕩衰減段。并且在整個振蕩數(shù)據(jù)中包含多階模態(tài),各模態(tài)之間存在能量傳遞,各模態(tài)幅值變化率也隨時間變化。因此,構(gòu)造符合上述特征的典型的振蕩信號用于校驗SPSO-VMD模態(tài)分解的準確性。此時,各階模態(tài)壓強振蕩幅值Rm滿足的非線性微分方程[13]為
(2)
式中αm為m階模態(tài)壓強振蕩的總增益系數(shù);A為非線性系數(shù);ωm為m階聲振模態(tài)的角頻率;Emnl為表示各階模態(tài)之間的非線性耦合作用。
按上式構(gòu)造含20階模態(tài)的仿真信號,前5階模態(tài)的幅值以及幅值變化率如表1所示,模態(tài)數(shù)m=20。6~20階初始幅值均為0,其幅值變化率滿足αm=(m-3)×α4。在表1中,正值αm表示能量通過線性作用從燃燒場流入振蕩系統(tǒng)中的第m階模態(tài),負值αm則表示能量從振蕩系統(tǒng)中流出,高階模態(tài)的αm的絕對值越大,表明高階模態(tài)越難被激發(fā)以及越容易被耗散。按上式構(gòu)造含20階模態(tài)的仿真信號如圖2所示。從圖2可看出,上升段振蕩幅值出現(xiàn)了明顯線性向非線性轉(zhuǎn)變的過程,即各階模態(tài)的幅值并不是以線性的方式持續(xù)增長,當(dāng)各階模態(tài)能量傳遞達到動態(tài)平衡過程時,出現(xiàn)了極限環(huán)狀態(tài);下降段各階模態(tài)的幅值變化率也是隨時變化各不相同;圖2(b)為整個振蕩段頻譜分析的前五階結(jié)果,前三階振蕩模態(tài)的幅值占主導(dǎo),所以這里僅分解前三階模態(tài)來驗證分解結(jié)果的準確性。
表1 前五階模態(tài)的幅值以及幅值變化率
(a)Amplitude of the simulated signal
對上述信號在給定α和τ范圍為(0,3000]和[0,0.2]內(nèi)進行SPSO-VMD方法的分解,優(yōu)化所得α和τ分別為2983.78和0.000 4。各階模態(tài)幅值分解結(jié)果如圖3(b)所示。從圖3(b)可以看分解后的各階模態(tài)幅值與圖2(c)幅值輪廓線高度一致;表2展示了SPSO-VMD方法分解該特征信號的頻率和幅值結(jié)果。
表2 SPSO-VMD方法分解特征信號頻率、幅值結(jié)果
(a)Iterative process
這里三個模態(tài)的分解所得極限環(huán)幅值是通過統(tǒng)計的方式計算所得,其與標準極限環(huán)幅值誤差均低于0.5%。從表中結(jié)果可以看出,對于非線性燃燒不穩(wěn)定典型特征信號,SPSO-VMD方法能夠準確的分離各個模態(tài)、辨識各階頻率和幅值。比較了以下兩種方法:
(1)采用VMD算法對特征信號進行分解,分別取值α和τ為16和0.01;
(2)采用PSO算法以同樣思路在同樣參數(shù)范圍內(nèi)對特征信號進行優(yōu)化VMD參數(shù)的分解,對比結(jié)果如表3和圖4所示。
圖4 Mode 1極限環(huán)段振蕩某點誤差示意圖
表3 三種算法對比
為求出全局最小適應(yīng)度函數(shù)值f(x),采用VMD法對特征信號進行分解。其中,α和τ的取值范圍分別為(0,3000]和[0,0.2],α的步長為1,τ的步長為0.01,全局求解共耗時約10 h,結(jié)果如圖5所示。
圖5 f(x)的全局求解
如表3和圖4所示,PSO-VMD方法和SPSO-VMD方法分解的誤差遠小于VMD方法。從圖4中可知,SPSO-VMD方法的誤差略小于PSO-VMD方法,但SPSO-VMD方法可以全局優(yōu)化到與最小適應(yīng)度函數(shù)值相對應(yīng)的參數(shù),且SPSO-VMD方法所需時間更短,效率更高,如圖3(a)所示,在本次計算中,SPSO-VMD方法僅用3步就取得了比PSO方法更好的效果。通過兩者誤差的對比分析,可以認為在適應(yīng)度函數(shù)為10-5數(shù)量級時,分解便已經(jīng)達到了較好的效果。
某發(fā)動機在工作過程中發(fā)生了非線性燃燒不穩(wěn)定,其振蕩幅值數(shù)據(jù)如圖6(a)所示,對其進行FFT處理。頻譜結(jié)果表明,該發(fā)動機在工作過程中出現(xiàn)了典型的非線性燃燒不穩(wěn)定現(xiàn)象,存在多階模態(tài),且前8階模態(tài)幅值明顯,通過計算分析,分解前8階模態(tài)重構(gòu)該信號的精度能達到95%,于是這里使用SPSO-VMD方法對其振蕩數(shù)據(jù)進行前8階模態(tài)分解,所用計算資源CPU型號為intel Xeon(R)Platinum 8259CL、48核,迭代在第9步時,適應(yīng)度函數(shù)便達到了7.763 5×10-5,通過上一節(jié)的誤差分析,可以認為該適應(yīng)度值表明分解已經(jīng)達到了較好的效果,對應(yīng)用時約為1880 s。
(a)Oscillation signal
8階模態(tài)的分解結(jié)果如圖7(a)~(h)所示,各分解模態(tài)的FFT結(jié)果表明,各階模態(tài)被準確分離。
(a)Mode 1 and its FFT result (b)Mode 2 and its FFT result
針對于類似第6階模態(tài),頻譜圖上出現(xiàn)兩個臨近頻率峰值的情況,這里作一定的分析:對第6階模態(tài)進行STFT處理,所得結(jié)果如圖8所示,結(jié)果表明第6階模態(tài)在發(fā)動機發(fā)生燃燒不穩(wěn)定的時間內(nèi)出現(xiàn)了頻率的變化,頻譜圖中兩個峰值恰好對應(yīng)STFT結(jié)果中振蕩幅值最大的0.29 s和0.46 s兩個時刻下的頻率,也即這兩個頻率是第6階模態(tài)不同時刻的振蕩頻率,它們屬于一個模態(tài)。這種頻率變化的現(xiàn)象主要是由燃燒室裝藥燃燒,燃面退移,導(dǎo)致燃燒室內(nèi)聲腔構(gòu)型發(fā)生變化造成的。
圖8 第6階模態(tài)STFT處理結(jié)果 圖9 前8階模態(tài)重構(gòu)信號
這里不考慮9階及其高階模態(tài),將得到的前8階模態(tài)重新疊加,與原始的振蕩幅值進行對比,如圖9所示,計算前8階模態(tài)重構(gòu)振蕩幅值的精度為95.2%。如圖10所示,將前8階模態(tài)繪制在同一張圖中,在0.16~0.4 s時間段范圍內(nèi),壓強信號雖然保持極限環(huán)振蕩現(xiàn)象,但是各階模態(tài)的幅值卻在不斷變化,可看到奇數(shù)階模態(tài)幅值衰減到了一個極小值,而偶數(shù)階模態(tài)出現(xiàn)了一定的增長,仿佛出現(xiàn)了偶數(shù)階增長抑制奇數(shù)階的現(xiàn)象,據(jù)文獻[15]所提:“非線性能量可從一階聲模態(tài)傳遞至所有高階聲模態(tài),而從二階聲模態(tài)傳遞只能向偶數(shù)階傳遞,奇數(shù)階模態(tài)均不發(fā)生能量傳遞過程,且高階聲模態(tài)不能向低階聲模態(tài)傳遞能量。”在該發(fā)動機振蕩信號中,由于某種原因,一階模態(tài)幅值衰減到了接近0的值,傳遞向高階聲模態(tài)的能量極少,而二階聲模態(tài)幅值較大,不斷向偶數(shù)階模態(tài)傳遞著非線性能量,導(dǎo)致偶數(shù)階維持一定幅值的振蕩,而奇數(shù)階高階模態(tài)因為輸入的非線性能量較少,不足以維持自身耗散的作用,導(dǎo)致幅值一直處在一個較低的值。這種現(xiàn)象恰好也提示研究人員可以從能量傳遞的角度對發(fā)動機非線性燃燒不穩(wěn)定進行抑制,比如切斷能量傳輸?shù)穆窂交蛘吒淖兡芰總鬏數(shù)牧康鹊?這些都有待研究。
圖10 前8階模態(tài)奇偶數(shù)階幅值演化對比
(1)SPSO-VMD方法能夠準確地分解發(fā)動機非線性燃燒不穩(wěn)定典型特征信號的各階模態(tài),獲得各模態(tài)的頻率和幅值,可應(yīng)用于真實發(fā)動機燃燒不穩(wěn)定振蕩信號的處理,可為研究非線性不穩(wěn)定能量傳遞演化機制提供數(shù)據(jù)支撐。
(2)SPSO-VMD方法改進了PSO-VMD方法陷入局部最優(yōu)解的缺陷,能夠準確地找到全局最優(yōu)解,分解模態(tài)誤差較小,分解效率較高。
(3)SPSO-VMD方法自適應(yīng)強,不依賴主觀經(jīng)驗,集成性好,使用簡單方便,可適用于工業(yè)部門對發(fā)動機燃燒不穩(wěn)定問題的定位和分析。