渦扇發(fā)動機(jī)故障診斷中粒子濾波改進(jìn)方法

盧俊杰，黃金泉，魯 峰

（南京航空航天大學(xué)能源與動力學(xué)院，南京210016）

0 引言

渦扇發(fā)動機(jī)結(jié)構(gòu)復(fù)雜，工作環(huán)境惡劣，屬于1種故障多發(fā)系統(tǒng)。據(jù)統(tǒng)計，渦扇發(fā)動機(jī)氣路部件故障占渦扇發(fā)動機(jī)總體故障的90%以上，因此，實時檢測發(fā)動機(jī)健康狀況并進(jìn)行氣路性能分析是提高發(fā)動機(jī)安全可靠性的重要途徑[1-3]。發(fā)動機(jī)氣路部件的效率變化系數(shù)、流量變化系數(shù)等健康參數(shù)是發(fā)動機(jī)氣路故障的狀態(tài)特征，會直接導(dǎo)致轉(zhuǎn)速、溫度、壓力等發(fā)動機(jī)測量參數(shù)的變化，因此發(fā)動機(jī)氣路部件故障診斷主要采用特定的算法根據(jù)測量參數(shù)的變化來估計出健康參數(shù)[4]，從而分析發(fā)動機(jī)氣路部件的健康狀況。發(fā)動機(jī)故障診斷方法主要有基于模型（如最小二乘方法[5]，卡爾曼濾波方法[6]），基于數(shù)據(jù)（如神經(jīng)網(wǎng)絡(luò)[7]，支持向量機(jī)[8]）和基于知識（如專家系統(tǒng)）的方法。隨著發(fā)動機(jī)部件級模型的精度以及計算機(jī)速度的提高，基于模型的故障診斷方法在工程實際中得到了廣泛的研究和應(yīng)用。

卡爾曼濾波算法受限于線性高斯系統(tǒng)，一些改進(jìn)的卡爾曼算法[9-10]雖然適用于非線性問題，但依然依賴于高斯噪聲的假設(shè)[11]。而粒子濾波算法適用于含非高斯噪聲的非線性系統(tǒng)，于是有學(xué)者提出使用粒子濾波算法進(jìn)行氣路部件健康參數(shù)的估計，并取得較好效果[12]。標(biāo)準(zhǔn)粒子濾波用于發(fā)動機(jī)故障診斷主要存在以下2方面問題：（1）標(biāo)準(zhǔn)粒子濾波直接采用先驗分布作為重要性密度函數(shù)[13]，沒有結(jié)合最新觀測值，導(dǎo)致在突變故障下，需要診斷步數(shù)多；（2）發(fā)動機(jī)健康參數(shù)的維度高，濾波器的抽樣率低，同時出于實時性考慮，粒子數(shù)少，所以有效粒子數(shù)少，導(dǎo)致診斷結(jié)果噪聲水平高，估計精度低。

本文針對粒子濾波應(yīng)用于突變故障診斷時診斷步數(shù)多和診斷結(jié)果噪聲水平高的問題，通過引入偽協(xié)方差與自適應(yīng)調(diào)整似然分布尾部的平坦程度，提出使用偽協(xié)方差和自適應(yīng)似然分布結(jié)合的1種改進(jìn)粒子濾波方法。

1 標(biāo)準(zhǔn)粒子濾波算法

粒子濾波是1種基于遞推貝葉斯估計和蒙特卡羅積分的統(tǒng)計濾波方法，其主要思想是：首先，由系統(tǒng)狀態(tài)向量的經(jīng)驗條件分布在狀態(tài)空間中產(chǎn)生1組隨機(jī)的粒子；其次，根據(jù)新的觀測值調(diào)整粒子的位置和權(quán)值；最后，通過調(diào)整后的粒子信息修正最初的條件分布[14]。實質(zhì)即為根據(jù)粒子及其權(quán)值組成的離散隨機(jī)測度來近似狀態(tài)的后驗概率分布，且這種近似會隨著粒子的數(shù)目趨于無窮多時收斂于真實后驗概率分布。

考慮如下的非線性系統(tǒng)

式中：xk、uk、zk分別為k時刻的狀態(tài)量、控制量和觀測量；ωk、vk分別為k時刻的過程噪聲和觀測噪聲；f（·）、h（·）分別為狀態(tài)轉(zhuǎn)移函數(shù)和觀測函數(shù)。

k時刻目標(biāo)狀態(tài)xk的后驗概率分布p(xk|z1,…zk)可以用{來描述，其中{是權(quán)值分別為{的粒子集，N為采樣粒子數(shù)。設(shè)歸一化權(quán)值為則k時刻系統(tǒng)狀態(tài)的后驗概率分布近似離散加權(quán)為

標(biāo)準(zhǔn)粒子濾波算法的重要性密度函數(shù)選擇為最容易實現(xiàn)的先驗概率密度

由歸一化權(quán)值加權(quán)求和得狀態(tài)的最優(yōu)估計值

2 改進(jìn)粒子濾波算法

高斯粒子濾波是標(biāo)準(zhǔn)粒子濾波的1種實現(xiàn)方法，通過高斯密度函數(shù)來逼近狀態(tài)后驗概率分布[15]，基本思想是利用反映粒子平均數(shù)值水平的數(shù)學(xué)期望和反映粒子離散程度的協(xié)方差進(jìn)行遞推來簡化粒子集的更新過程。本節(jié)在高斯粒子濾波實現(xiàn)方法的基礎(chǔ)上，提出1種改進(jìn)的粒子濾波算法。該方法結(jié)合發(fā)動機(jī)健康參數(shù)估計的特點，引入偽協(xié)方差的代替協(xié)方差，以期望減少突變故障診斷步數(shù)；另外通過對似然分布自適應(yīng)調(diào)整，以期望降低診斷結(jié)果噪聲水平。

2.1 偽協(xié)方差

在發(fā)動機(jī)故障診斷中，在無突變故障發(fā)生的情況下，粒子加權(quán)均值基本穩(wěn)定，可以近似認(rèn)為。粒子集的協(xié)方差∑k為

2.2 自適應(yīng)似然分布

在似然分布位于轉(zhuǎn)移先驗分布尾部或者觀測模型具有很高精度時[16]，很多樣本由于歸一化權(quán)重很小而成為無效樣本，過低的抽樣率會導(dǎo)致粒子濾波失敗，增大粒子數(shù)是其解決方法之一，但隨之而來的巨大計算量使該方法在工程應(yīng)用中無法實現(xiàn)。

發(fā)動機(jī)故障診斷需要估計多個健康參數(shù)，高維度必然導(dǎo)致似然分布處于轉(zhuǎn)移先驗尾部，如圖1所示。從圖中可見，很少的粒子位于B處，絕大部分粒子位于A處成為無效粒子。同時受實時性限制，粒子總數(shù)少，所以濾波器有效粒子數(shù)低，診斷結(jié)果的噪聲水平高。針對該問題，本節(jié)通過對似然分布進(jìn)行自適應(yīng)調(diào)整以提高抽樣率，使之尾部更為平坦，增加先驗和似然的重疊區(qū)，降低診斷結(jié)果的噪聲水平。改變似然函數(shù)的分布就是改變粒子權(quán)值的分布，所以只需對標(biāo)準(zhǔn)粒子濾波算法的權(quán)值更新作自適應(yīng)改進(jìn)。

圖1 狀態(tài)轉(zhuǎn)移函數(shù)和似然函數(shù)的關(guān)系

式中：β為比例因子，根據(jù)權(quán)值離散水平進(jìn)行調(diào)整；κ為擠壓因子。根據(jù)權(quán)值數(shù)值水平進(jìn)行調(diào)整，取β∈(0.6,1.5),κ∈(0.3,0.9)。

2.3 改進(jìn)算法實現(xiàn)

（3）計算各粒子的權(quán)值

（5）將各粒子的權(quán)值歸一化

（6）求狀態(tài)估計值x?k和粒子集偽協(xié)方差

3 發(fā)動機(jī)氣路故障的粒子濾波診斷

為了表示故障發(fā)生后部件性能的變化程度，引入旋轉(zhuǎn)部件的效率變化系數(shù)ΔSEi和流量變化系數(shù)ΔSWi[17]

式中：ηi、Wi為部件的實際效率和流量；為部件效率和流量的理想值。

健康參數(shù)選為風(fēng)扇、壓氣機(jī)、高壓渦輪和低壓渦輪的效率和流量的變化系數(shù)[18]，共8個，定義為

將其增廣到發(fā)動機(jī)狀態(tài)中，則

式中：NL、NH分別為低、高壓轉(zhuǎn)子轉(zhuǎn)速；ΔSE1、ΔSE2、ΔSE3、ΔSE4分別為風(fēng)扇、壓氣機(jī)、高壓渦輪、低壓渦輪效率變化系數(shù)；ΔSW1、ΔSW2、ΔSW3、ΔSW4分別為風(fēng)扇、壓氣機(jī)、高壓渦輪、低壓渦輪流量變化系數(shù)；上標(biāo)T表示將行向量轉(zhuǎn)換為列向量的變換，涉及的控制量有：Wfb為燃油量，A8為尾噴管面積。

基于非線性模型的氣路部件故障診斷原理如圖2所示。通過發(fā)動機(jī)的輸出值與模型預(yù)測值之間的殘差結(jié)合粒子濾波算法對氣路部件健康參數(shù)進(jìn)行估計。發(fā)動機(jī)傳感器參數(shù)選取見表1，相應(yīng)傳感器的噪聲水平根據(jù)文獻(xiàn)[12，17]確定。

圖2 基于粒子濾波的故障診斷原理

表1 傳感器選取及其測量噪聲

4 仿真與分析

高斯分布為統(tǒng)計學(xué)中常用的理想噪聲模型，然而現(xiàn)實應(yīng)用中噪聲通常都是比較復(fù)雜的，為了模擬發(fā)動機(jī)運行時的復(fù)雜噪聲，本文選取由伽馬噪聲以及高斯噪聲簡單混合組成的雙模噪聲，對發(fā)動機(jī)健康參數(shù)突變故障進(jìn)行估計，以驗證粒子濾波在非高斯噪聲下狀態(tài)參數(shù)估計的準(zhǔn)確性。為了比較的公平性，標(biāo)準(zhǔn)粒子濾波算法以及改進(jìn)粒子濾波算法的粒子數(shù)均設(shè)為50。另外，設(shè)～N(μm，)～Γ(μm，)，觀測噪聲，其中 m∈{1,2,…，8}表示對應(yīng)傳感器，和相互獨立。根據(jù)表1中傳感器測量噪聲標(biāo)準(zhǔn)差取R值，令r=1，μ=-，σ=，=mmmm，則 E()=0，D()=Rm。過程噪聲參考文獻(xiàn)[12]選取為對角陣 Q=0.00152×I10×10。

在地面標(biāo)準(zhǔn)狀況的最大工作狀態(tài)下，分別用標(biāo)準(zhǔn)粒子濾波器和改進(jìn)粒子濾波器對發(fā)動機(jī)氣路部件的單部件故障模式和多部件故障模式進(jìn)行仿真估計。參考NASA在MAPSS仿真平臺中渦扇發(fā)動機(jī)完成一定工作循環(huán)數(shù)量后氣路部件性能參數(shù)變化情況的統(tǒng)計數(shù)據(jù)[19]，選擇單部件故障模式：風(fēng)扇效率變化系數(shù)突變-0.04；多部件故障模式：風(fēng)扇效率變化系數(shù)突變-0.04，壓氣機(jī)流量變化系數(shù)突變-0.02。仿真時間為10 s，突變故障發(fā)生時刻均設(shè)定為第2 s，采樣步長為20 ms。綜合估計精度用全部仿真時間10 s內(nèi)的均方根誤差RMSE來衡量，其表達(dá)式為

式中：定義診斷步數(shù)為從故障突變到健康參數(shù)穩(wěn)定在故障值±0.01范圍內(nèi)的仿真步數(shù)差，用來反映診斷速度；Δh為診斷結(jié)果穩(wěn)定后的平均值，σh為診斷結(jié)果的標(biāo)準(zhǔn)差。

在單故障模式下，標(biāo)準(zhǔn)粒子濾波和改進(jìn)粒子濾波的20次運行平均性能見表2，故障診斷值及其標(biāo)準(zhǔn)差見表3并如圖3所示。在多故障模式下，標(biāo)準(zhǔn)粒子濾波和改進(jìn)粒子濾波的20次運行平均性能見表4，故障診斷值及其標(biāo)準(zhǔn)差見表5并如圖4所示。

表2 單故障模式下濾波性能對比

表3 單故障模式下估計值及標(biāo)準(zhǔn)差

圖3 單故障模式診斷結(jié)果

表4 多故障模式下濾波性能對比

表5 多故障模式下估計值及標(biāo)準(zhǔn)差 %

圖4 多故障模式診斷結(jié)果

在表3、5中Δh和σh均為20次運行的均值，而σh的均值為8個健康參數(shù)的標(biāo)準(zhǔn)差的均值，綜合體現(xiàn)診斷結(jié)果噪聲水平。圖3、4中每個健康參數(shù)對應(yīng)的3 個點分別為 Δh-σh，Δh，Δh+σh，圖中線段長度反映該健康參數(shù)估計值的噪聲水平。

由于改進(jìn)粒子濾波采用更能反映突變情況的偽協(xié)方差代替協(xié)方差陣，抽樣得到的粒子能夠更快地逼近真實值，從表2、4中可見，突變故障診斷步數(shù)明顯減少，2種故障模式下診斷步數(shù)均減少約27%；另外由于采用自適應(yīng)的似然分布函數(shù)，提高了先驗分布和似然分布的重疊區(qū)域，增加了有效粒子數(shù)，從表3、5中可見，診斷結(jié)果噪聲水平顯著下降，在2種故障模式下，估計結(jié)果的平均標(biāo)準(zhǔn)差均降低約39%。從圖3、4中可見，改進(jìn)粒子濾波器的8個維度的健康參數(shù)診斷結(jié)果噪聲水平均明顯下降，估計值更接近真實值。由于改進(jìn)粒子濾波在故障診斷速度和診斷結(jié)果噪聲水平這2個方面的改善，綜合估計誤差明顯減小，從表2、4中可見，2種故障模式下RMSE均減小約38%。

5 結(jié)論

標(biāo)準(zhǔn)粒子濾波算法在發(fā)動機(jī)故障診斷領(lǐng)域有著廣泛應(yīng)用，但該算法存在著突變故障診斷步數(shù)多以及診斷結(jié)果噪聲水平高的問題，本文提出使用偽協(xié)方差和自適應(yīng)似然分布結(jié)合的1種改進(jìn)粒子濾波方法。在發(fā)動機(jī)故障診斷中的仿真驗證表明，相比于標(biāo)準(zhǔn)粒子濾波算法，改進(jìn)方法能在相同粒子數(shù)的前提下，顯著減小診斷的均方根誤差，并且減少突變故障的診斷步數(shù)。綜合考慮單故障模式和多故障模式，改進(jìn)算法能使診斷速度提高約27%，診斷精度提高約38%。優(yōu)越的性能表明，該改進(jìn)粒子濾波方法可以進(jìn)一步推廣應(yīng)用于發(fā)動機(jī)故障診斷系統(tǒng)中。