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

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

0 引言

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

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

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

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

1 無(wú)跡粒子濾波算法

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

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

1.1 粒子濾波框架

對(duì)于航空發(fā)動(dòng)機(jī)EGTM 預(yù)測(cè)，可以用如下的數(shù)學(xué)狀態(tài)空間方程來(lái)表示。

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

系統(tǒng)觀測(cè)方程

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

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

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

對(duì)于一般高階非線性、非高斯過(guò)程，式（4）難以求解。所以粒子濾波采用蒙特卡洛采樣來(lái)避免積分難的問(wèn)題。利用一系列從已知分布采樣的粒子來(lái)估計(jì)后驗(yàn)概率

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

然而，因?yàn)闊o(wú)法獲知真實(shí)的狀態(tài)后驗(yàn)概率，所以如果需要通過(guò)后驗(yàn)概率來(lái)獲取粒子分布，則需要從1個(gè)已知分布q(x|y)中采樣得到粒子

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

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

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

1.2 無(wú)跡粒子濾波算法

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

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

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

（1）確定2n+1個(gè)sigma點(diǎn)Si={Xi,Wi}。

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

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

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

（2）通過(guò)非線性函數(shù)傳播sigma點(diǎn)。

（3）計(jì)算y的均值和方差。

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

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

（1）初始化。

（2）按時(shí)間迭代。

a.用上述過(guò)程計(jì)算sigma點(diǎn)。

b.時(shí)間更新。

c.測(cè)量更新。

UKF 為2 階近似，所以在準(zhǔn)確度上要?jiǎng)龠^(guò)1 階近似的擴(kuò)展卡爾曼濾波器（Extended Kalman Filter，EKF），而且在計(jì)算效率上同樣優(yōu)于EKF，因?yàn)閁KF不需要計(jì)算雅克比矩陣和海森矩陣。

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

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

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

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

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

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

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

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

圖1 UPF的流程

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

2.1 EGTM預(yù)測(cè)過(guò)程

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

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

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

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

EGTM的后驗(yàn)概率密度為

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

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

2.2 仿真結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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