基于無跡粒子濾波算法的航空發(fā)動機排氣溫度預(yù)測

余 臻 ，劉 洋 ，魏 芳 ，劉利軍 ，

（1.廈門大學(xué)航空航天學(xué)院，福建廈門 361005；2.中國航發(fā)商用航空發(fā)動機有限責(zé)任公司，上海 200241；3.中國船舶航海保障技術(shù)重點實驗室，天津 300450）

0 引言

航空發(fā)動機的工作環(huán)境惡劣，多為高溫、高腐蝕、高轉(zhuǎn)速條件，其安全工作壽命難以準(zhǔn)確得知[1]。發(fā)動機性能指標(biāo)包括排氣溫度裕度（Exhaust Gas Temperature Margin，EGTM）、高壓壓氣機轉(zhuǎn)子角速度偏差量（Delta N2，DN2）、平均滑油消耗率（Average Oil Consumption，AOC）等，其中EGTM 是最重要的指標(biāo)之一，是發(fā)動機在海平面壓力、拐點溫度條件下全功率起飛時，排氣溫度與標(biāo)準(zhǔn)規(guī)定的排氣溫度紅線值之間的差值。隨著發(fā)動機飛行循環(huán)數(shù)的增加，發(fā)動機各零部件磨損老化程度增加，排氣溫度持續(xù)升高，使得EGTM值逐漸降低，達到標(biāo)準(zhǔn)規(guī)定的EGTM閾值。

在20世紀90年代初建立的基于發(fā)動機性能參數(shù)的性能衰退失效分析方法，為高可靠性設(shè)備的可靠性分析提供了新的分析方法。Gertsbackh 等[2]率先提出基于性能退化數(shù)據(jù)的設(shè)備可靠性模型的預(yù)測方法；Nelson[3]發(fā)現(xiàn)發(fā)動機部件的失效特性會在發(fā)動機衰退信息中有所體現(xiàn)，這些信息能具體反映出零部件的衰退狀態(tài)；Meeker 等[4]利用退化失效模型解決了傳統(tǒng)可靠性理論無法很好處理的一些實際問題；Wu 等[5]基于數(shù)據(jù)驅(qū)動的退化量統(tǒng)計模型，采用模糊聚類方法研究了零部件退化失效分布。隨著發(fā)動機制造技術(shù)的進步，發(fā)動機可靠度提高，發(fā)動機失效數(shù)據(jù)大大減少，使得樣本數(shù)量減少。Erto[6]利用發(fā)動機零部件常見的先驗信息，分別對經(jīng)典的概率分布模型進行分析，對參數(shù)進行了估計；Cox[7]給出了針對單系統(tǒng)的剩余壽命分布的混合估計方法；王鑫等[8]提出了LSTM 預(yù)測模型參數(shù)優(yōu)選算法，提高了對故障時間序列預(yù)測的準(zhǔn)確度；張營等[9]在過程神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上提出了優(yōu)化算法，結(jié)合混沌粒子群算法對發(fā)動機的排氣溫度進行預(yù)測，獲得較高的預(yù)測精度。

發(fā)動機的性能衰退主要表現(xiàn)為其性能衰退參數(shù)值呈非線性減小趨勢[10-11]，可采用能夠處理非線性時間序列的方法進行預(yù)測。

本文提出了一種基于粒子濾波[12-14]（Particle Filter，PF）的航空發(fā)動機剩余大修壽命預(yù)測方法。提出的無跡粒子濾波（Unscented Particle Filter，UPF）算法首次應(yīng)用在航空發(fā)動機EGTM 預(yù)測領(lǐng)域，UPF 以無跡卡爾曼濾波（Unscented Kalman Filter，UKF）的結(jié)果作為建議分布[15]，引入最新觀測值進行預(yù)測修正。

1 無跡粒子濾波算法

粒子濾波基于蒙特卡羅和遞歸貝葉斯濾波方法[16]。為了解決粒子權(quán)重差異大的問題，有2 種策略，即重采樣和選擇合理的建議分布[17]。重采樣會導(dǎo)致粒子貧化，選擇合理的建議分布是解決退化問題的較優(yōu)選擇。

粒子濾波是由表示未知狀態(tài)空間的采樣值的1組粒子的近似。目前，粒子濾波已經(jīng)成為了解決非線性、非高斯系統(tǒng)參數(shù)估計和狀態(tài)濾波問題的主流方法。

1.1 粒子濾波框架

對于航空發(fā)動機EGTM 預(yù)測，可以用如下的數(shù)學(xué)狀態(tài)空間方程來表示。

系統(tǒng)狀態(tài)方程

系統(tǒng)觀測方程

式中：xt∈Rn，為系統(tǒng)的狀態(tài)變量；yt∈Rnz，為系統(tǒng)在t時刻的測量值；ft為通過上一時刻的狀態(tài)值預(yù)測下一時刻狀態(tài)值的函數(shù)；ht為通過當(dāng)前時刻狀態(tài)值預(yù)測輸出值的函數(shù)；mt、nt分別為系統(tǒng)的過程噪聲和觀測噪聲，二者之間相互獨立。

貝葉斯濾波包含預(yù)測和更新2 部分。預(yù)測部分是利用系統(tǒng)的狀態(tài)方程，在系統(tǒng)測量值未知的前提下，預(yù)測狀態(tài)的先驗概率密度

而更新部分則是在已經(jīng)得到最新的觀測值yt的條件下，對先驗的概率分布進行修正，得到xt的后驗概率

對于一般高階非線性、非高斯過程，式（4）難以求解。所以粒子濾波采用蒙特卡洛采樣來避免積分難的問題。利用一系列從已知分布采樣的粒子來估計后驗概率

式中：δ()為狄拉克函數(shù)；N為粒子個數(shù)；i為粒子序號。在得到后驗概率后，當(dāng)前時刻的后驗期望值為。

然而，因為無法獲知真實的狀態(tài)后驗概率，所以如果需要通過后驗概率來獲取粒子分布，則需要從1個已知分布q(x|y)中采樣得到粒子

上述粒子權(quán)重的計算采用直接計算方式，每進行1 次迭代就需要重新計算，效率低。所以需要重新對權(quán)重公式進行遞推

至此描述了完整的粒子濾波算法。從上述算法得知，建議分布q(xt|y1:t)的選取對于算法的效果非常重要。上述涉及到建議分布的內(nèi)容主要包含2 部分：一部分是粒子是從建議分布中選取的，另一部分是粒子的權(quán)重計算需要用到建議分布。同時可知，當(dāng)建議分布選取為系統(tǒng)先驗分布p(xt|xt-1)時，權(quán)重的計算復(fù)雜度會大幅降低，即為

這極大地促進了粒子濾波的實際應(yīng)用。但粒子濾波在易于實現(xiàn)的同時也會導(dǎo)致其他后果：（1）導(dǎo)致了更高階的蒙特卡洛協(xié)方差，進而使預(yù)測效果變差；（2）使用先驗概率作為建議分布之后，忽略了最新的系統(tǒng)觀測值，可能導(dǎo)致粒子嚴重退化，即大多數(shù)粒子的權(quán)值都趨近于0，尤其是在似然函數(shù)達到峰值而預(yù)測的狀態(tài)卻在似然函數(shù)的邊緣時。建議分布的取值在理論上有無窮多種可能性，但是最優(yōu)的分布-權(quán)值的方差最小，即使存在，求解也十分困難。因此，選擇1 個合理的建議分布，對于避免粒子退化、得到良好的濾波結(jié)果具有深刻意義。

1.2 無跡粒子濾波算法

上文已經(jīng)指出了使用先驗概率作為建議分布會導(dǎo)致系統(tǒng)缺陷，而改進建議分布最明顯的方法是納入當(dāng)前觀測的數(shù)據(jù)。于是學(xué)者設(shè)計出很多的卡爾曼濾波器來改進建議分布[18]，但是這些濾波器的表現(xiàn)多基于各自不同的近似假設(shè)。迄今為止，無跡卡爾曼濾波器是針對非線性系統(tǒng)表現(xiàn)最為優(yōu)秀的卡爾曼濾波器。通過采用無跡卡爾曼濾波器產(chǎn)生合適的建議分布，可以將通用粒子濾波轉(zhuǎn)化成高性能無跡粒子濾波。下文將討論無跡卡爾曼濾波的基礎(chǔ)——無跡變換，并且給出完整的無跡粒子濾波算法。

在許多現(xiàn)實應(yīng)用中，都需要估計經(jīng)過非線性變換y=g(x)的隨機變量的低階統(tǒng)計量，例如均值和協(xié)方差。無跡變換能夠準(zhǔn)確地計算出g(x)的泰勒級數(shù)展開不超過2階的均值和方差。

假定n為系統(tǒng)變量x的維度，x的均值為，方差為Px，則通過無跡變換計算得到y(tǒng)=g(x)的均值和方差。

（1）確定2n+1個sigma點Si={Xi,Wi}。

式中：κ為尺度參數(shù)，用來控制sigma 點與均值的距離；α為正的尺度參數(shù)，用以控制來自非線性函數(shù)g(x)的更高階部分的影響；β為用來控制初始的sigma點權(quán)重的參數(shù)。

縱觀現(xiàn)有的研究仍然存在著研究視角單一，對于組織邊界、組織域的關(guān)注較少，缺乏對于組織機制及運作邏輯的結(jié)構(gòu)性分析和建構(gòu)過程的討論。對于體育社會組織的制度約束、資源困境等關(guān)注過多，忽視了對于包括文化傳統(tǒng)、道德習(xí)俗等非正式制度對于組織運作機制的形塑過程的分析，缺乏對于中國傳統(tǒng)社會“合群立會”等歷史視角的社會學(xué)想象力。從研究方法上看，個案研究過多，定量研究較少，從理論的解釋力來看，理論與經(jīng)驗的不匹配、理論和方法背后的邏輯推演的不清晰，引進的西方理論的適用性等問題仍存在于目前的研究當(dāng)中。

對于標(biāo)量情況來說，最優(yōu)值為α=0,β=0,κ=2。請注意，在計算平均值和協(xié)方差時，初始的sigma點的權(quán)重是不同的。

（2）通過非線性函數(shù)傳播sigma點。

（3）計算y的均值和方差。

保證y的均值和協(xié)方差精確到泰勒級數(shù)2階展開。

無跡卡爾曼濾波可以通過擴展?fàn)顟B(tài)空間以包含噪聲分量來實現(xiàn)，即。而Na=Nx+Nm+Nn則是擴展之后的狀態(tài)空間的維度，其中Nm和Nn是噪聲mt和nt的維度，Q和R則是噪聲mt和nt的協(xié)方差。整個UKF過程如下。

（1）初始化。

（2）按時間迭代。

a.用上述過程計算sigma點。

b.時間更新。

c.測量更新。

UKF 為2 階近似，所以在準(zhǔn)確度上要勝過1 階近似的擴展卡爾曼濾波器（Extended Kalman Filter，EKF），而且在計算效率上同樣優(yōu)于EKF，因為UKF不需要計算雅克比矩陣和海森矩陣。

至此已討論傳統(tǒng)的粒子濾波及其不足之處，隨后引入了無跡卡爾曼濾波。然而，無跡卡爾曼濾波對狀態(tài)分布做了高斯假設(shè)。另外，對于粒子濾波，可以對任意分布進行建模，但是將新的觀測加入到建議分布中并不容易。本文提出了改進辦法，在航空發(fā)動機領(lǐng)域，利用UKF 來為粒子濾波生成建議分布，引進最新的觀察來更新狀態(tài)估計，即為無跡粒子濾波。對于每個粒子的取值概率分布明確為

值得注意的是，即使在估計后驗分布p(xt|xt-1,y0:t)時高斯分布假設(shè)是不現(xiàn)實的，但對于每個粒子都有1個不同的均值和方差都不是問題。而且，使用UKF得到的均值和方差均達到2 階展開，所以系統(tǒng)的非線性特性得以保存。將式（32）代入到傳統(tǒng)粒子濾波算法中，即可得到UPF。完整的預(yù)測步驟如下。

（1）利 用式（21）～（31）的UKF 算法 來更 新,i=1,...,N，得到粒子的均值和協(xié)方差。

（3）利用式（5）計算出各粒子的權(quán)重。

（4）利用式（6）來歸一化權(quán)重。

（5）計算有效的粒子規(guī)模大小S。

（6）如果有效粒子數(shù)少于有效數(shù)閾值，則對應(yīng)生成N個等權(quán)粒子。

（7）利用式（7）來計算系統(tǒng)狀態(tài)后驗分布的期望值。xt的條件均值可以通過xt+1=f(xt)計算得到，而xt的協(xié)方差可以通過計算得到?？梢允褂糜嬎愕玫降臈l件均值和協(xié)方差來得到發(fā)動機當(dāng)前時刻排氣溫度的概率密度分布[19-21]。UPF 的流程如圖1所示。

圖1 UPF的流程

2 EGTM的預(yù)測算法及實現(xiàn)

2.1 EGTM預(yù)測過程

系統(tǒng)的測量方程和狀態(tài)方程為

式中：a，b，c均為系數(shù)。

式（34）中的yt即為t時刻系統(tǒng)的觀測值，是航空發(fā)動機EGT 的測量值與EGT 紅線值的正差值，系統(tǒng)噪聲和測量噪聲均假定遵循高斯分布。根據(jù)粒子濾波方法，EGTM值估計為

所以在t時刻的p步預(yù)測表示為

EGTM的后驗概率密度為

當(dāng)設(shè)定的EGTM 閾值為Y時，達到閾值的循環(huán)為L，則

所以在t時刻預(yù)測發(fā)動機在現(xiàn)有EGTM 條件下得到循環(huán)閾值的概率密度分布為

2.2 仿真結(jié)果分析

在本次仿真試驗中，每次迭代中所取的粒子數(shù)均為100。采用無跡粒子濾波方法在18 和30 循環(huán)處對數(shù)據(jù)進行預(yù)測，狀態(tài)方程的參數(shù)見表1。

表1 根據(jù)歷史數(shù)據(jù)得到的狀態(tài)方程參數(shù)

35 循環(huán)處傳統(tǒng)粒子濾波方法的預(yù)測結(jié)果如圖2所示。

圖2 傳統(tǒng)粒子濾波在35循環(huán)處的預(yù)測

由于沒有考慮到最新的觀測值，采用傳統(tǒng)粒子濾波方法得到的數(shù)據(jù)與真實的發(fā)動機水洗循環(huán)數(shù)相差較多，預(yù)測概率密度是假定閾值為50 ℃時的數(shù)據(jù)，而當(dāng)閾值為35 ℃時，預(yù)測的循環(huán)數(shù)在120 循環(huán)之后，與實際值相差較多。

在18 循環(huán)時刻的無跡粒子濾波預(yù)測結(jié)果如圖3所示。

圖3 無跡粒子濾波在18循環(huán)處的預(yù)測

前18 循環(huán)的航空發(fā)動機EGTM 實際數(shù)據(jù)將用來更新模型。從圖中可見，在達到閾值35 ℃時，預(yù)測的最終循環(huán)數(shù)為60，而實際的飛行循環(huán)數(shù)為59，預(yù)測誤差為1 循環(huán)，同時預(yù)測的EGTM 在閾值時的循環(huán)數(shù)的寬度為16。

前30 循環(huán)的發(fā)動機EGTM 實際數(shù)據(jù)將用來更新模型，如圖4 所示。在達到閾值35 ℃時，預(yù)測的最終循環(huán)數(shù)為59，而實際的飛行循環(huán)數(shù)為59，預(yù)測誤差為0 循環(huán)，同時預(yù)測的EGTM 在閾值時的循環(huán)數(shù)的寬度為7。從上述試驗結(jié)果可見，因為有更多的實際數(shù)據(jù)用以更新模型，所以模型的準(zhǔn)確度得到了提升，模型的擬合優(yōu)度統(tǒng)計見表2。預(yù)測的循環(huán)數(shù)的寬度收窄了，說明預(yù)測的精度更高了。

圖4 無跡粒子濾波在30循環(huán)處的預(yù)測

表2 擬合優(yōu)度統(tǒng)計

3 結(jié)束語

隨著民航產(chǎn)業(yè)的飛速發(fā)展，做好航空發(fā)動機性能衰退預(yù)測，對于發(fā)動機有效及時恢復(fù)性能和節(jié)約成本具有重要意義。針對發(fā)動機領(lǐng)域的特點，本文所提出的無跡粒子濾波算法作為一種改進的粒子濾波方法，充分利用了最新觀測數(shù)據(jù)對系統(tǒng)模型進行修正，也降低了重采樣對樣本貧化的影響。因此，該方法不僅有抑制粒子退化的優(yōu)點，而且在一定程度上保持了粒子的多樣性。通過對無跡粒子濾波預(yù)測結(jié)果和傳統(tǒng)粒子濾波計算結(jié)果的比較研究可知，無跡粒子濾波方法在預(yù)測的準(zhǔn)確度（均值與真實值的接近程度）和預(yù)測的置信度（概率密度分布的狹窄程度）上都顯示出巨大的優(yōu)勢。