基于超統(tǒng)計的多階段航空發(fā)動機(jī)剩余壽命預(yù)測

劉君強，胡東斌，潘春露，雷凡，趙倩茹

（1.南京航空航天大學(xué) 民航學(xué)院，南京211000； 2.南京航空航天大學(xué) 外國語學(xué)院，南京211000）

航空發(fā)動機(jī)的復(fù)雜性與精密性使之成為“工業(yè)之花”，而其安全與高效工作直接決定了航空運輸?shù)陌踩c效益，因此對它進(jìn)行健康管理及剩余壽命預(yù)測至關(guān)重要?，F(xiàn)代航空器的采購費用和使用保障費用日益龐大，據(jù)美軍官方統(tǒng)計，在武器裝備的全壽命周期費用中，使用與保障費用占到了總費用的72%［1］。與使用保障費用相比，維修保障費用在技術(shù)上更具有可壓縮性。而故障預(yù)測與健康管理（Prognostic and Health Management，PHM）正是壓縮維修保障費用的一個重要手段［2］。

在民航發(fā)動機(jī)的運行保障中，傳統(tǒng)的方法是憑借檢查發(fā)動機(jī)以判斷其運行狀況。近年來，越來越多的國內(nèi)外研究員通過先進(jìn)的PHM 技術(shù)對航空發(fā)動機(jī)的健康狀態(tài)進(jìn)行深入研究，主要內(nèi)容包括發(fā)動機(jī)健康監(jiān)測［3］、發(fā)動機(jī)失效機(jī)理［4］、發(fā)動機(jī)壽命預(yù)測等方面［5］。

在剩余壽命（Remaining Useful Life，RUL）預(yù)測模型研究中，主要包括基于數(shù)據(jù)驅(qū)動的方法、基于物理模型的方法和融合的方法。然而，航空發(fā)動機(jī)因設(shè)計結(jié)構(gòu)復(fù)雜導(dǎo)致其物理模型難以確定，基于物理模型的方法和融合的方法，目前的研究較少。

大多研究人員采用基于數(shù)據(jù)驅(qū)動的方法對航空發(fā)動機(jī)的RUL進(jìn)行預(yù)測。任淑紅和左洪福［6］基于排氣溫度裕度（EGTM）對航空發(fā)動機(jī)進(jìn)行RUL預(yù)測，但是僅通過單參數(shù)的RUL預(yù)測不能客觀反映整臺發(fā)動機(jī)的健康狀況。趙廣社等［7］基于多源數(shù)據(jù)對航空發(fā)動機(jī)進(jìn)行RUL預(yù)測。張馬蘭［8］發(fā)現(xiàn)航空發(fā)動機(jī)RUL退化曲線前后階段呈不同的退化模式，采用Kalman濾波（KF）和粒子濾波算法對其進(jìn)行研究。Ram in和M ing［9］采用集成多參數(shù)的方法進(jìn)行發(fā)動機(jī)RUL預(yù)測。綜上可知，發(fā)動機(jī)RUL預(yù)測已經(jīng)從單參數(shù)方法逐漸過渡到多參數(shù)方法。

Baraldi等［10］提出一種基于KF模型的退化預(yù)測系統(tǒng)進(jìn)行RUL的預(yù)測。Cavarzere和Venturini［5］通過4種不同的方法對比對航空發(fā)動機(jī)的性能進(jìn)行預(yù)測。An等［11］使用粒子濾波器來估計退化模型的參數(shù)。

對于航空發(fā)動機(jī)退化過程存在多階段性的問題，張馬蘭［12］將航空發(fā)動機(jī)退化過程分為2個階段；黃亮等［13］基于W iener過程開始對發(fā)動機(jī)進(jìn)行兩階段RUL預(yù)測的研究。此外，國內(nèi)較少有學(xué)者對發(fā)動機(jī)的多階段退化及RUL預(yù)測問題進(jìn)行研究。

因此，本文提出了一個基于超統(tǒng)計的多階段航空發(fā)動機(jī)RUL預(yù)測模型，該模型可通過多源監(jiān)測參數(shù)和突變點檢測，將發(fā)動機(jī)退化過程劃分為若干個退化階段，從而精確預(yù)測發(fā)動機(jī)的RUL，并對該模型的收斂性給予了相應(yīng)證明?；谠撃Ｐ停疚奶岢隽讼鄳?yīng)的算法，該算法先對航空發(fā)動機(jī)各個參數(shù)的時間序列進(jìn)行多階段退化識別；再采用無跡卡爾曼濾波（UKF）對融合的時變參數(shù)進(jìn)行濾波處理；最后通過非線性擬合對發(fā)動機(jī)RUL進(jìn)行預(yù)測。本文通過美國NASA提供的渦扇發(fā)動機(jī)數(shù)據(jù)對該算法的有效性進(jìn)行了實驗驗證，結(jié)果表明，該算法在發(fā)動機(jī)性能退化中的預(yù)測具有較好的適應(yīng)性，能更準(zhǔn)確地預(yù)測發(fā)動機(jī)的RUL。

1 航空發(fā)動機(jī)RUL預(yù)測模型

針對傳統(tǒng)航空發(fā)動機(jī)RUL預(yù)測模型無法客觀描述多階段性能衰退過程及對于RUL預(yù)測精度不高的問題，本文提出了航空發(fā)動機(jī)RUL預(yù)測模型以取得更準(zhǔn)確的發(fā)動機(jī)RUL估計，該模型包含以下4部分：超統(tǒng)計理論、突變點檢測、UKF、非線性預(yù)測。

1.1 模型組成

1）超統(tǒng)計理論

超統(tǒng)計理論作為統(tǒng)計物理學(xué)的一個分支，致力于研究非線性與非平衡系統(tǒng)，通過疊加多個不同的統(tǒng)計模型來描述目標(biāo)模型的特征［14］。本文研究的是分布函數(shù)不變、分布參數(shù)變化的非平穩(wěn)時間序列。將復(fù)雜的航空發(fā)動機(jī)系統(tǒng)抽象為2個部分的疊加，其一是相對微觀的穩(wěn)定平衡系統(tǒng)，即某時刻發(fā)動機(jī)的各個指標(biāo)；另一部分是對應(yīng)的宏觀系統(tǒng)，即發(fā)動機(jī)整體性能，該系統(tǒng)服從一定統(tǒng)計分布F（x），隨著使用時間的增加而緩慢變化。通過微觀系統(tǒng)與宏觀系統(tǒng)模型的疊加，可以描述隨時間變化的復(fù)雜系統(tǒng)的分布模型。本文通過超統(tǒng)計理論對航空發(fā)動機(jī)的健康狀態(tài)進(jìn)行分析建模。

從微觀角度建模。運用統(tǒng)計學(xué)方法對發(fā)動機(jī)指標(biāo)變化情況進(jìn)行分析，一般情況下，在環(huán)境參數(shù)ν一定的情況下，發(fā)動機(jī)健康狀況異常的條件概率密度為

式中：Γ（ν）為Gamma函數(shù)；ν由發(fā)動機(jī)的運行狀況決定，是一個受客流量大小、氣候條件、運營維護(hù)（維修人員的故障檢測、故障排除）［15］、線路設(shè)備故障等影響因素的正值參數(shù)。

從宏觀角度建模。健康指標(biāo)的狀況是隨著時間的推移而不斷變化的，由于在時間維度上環(huán)境參數(shù)ν是變動的，這使得在較大的時間尺度內(nèi)，健康指標(biāo)的受損分布模型是一個與時間變量t相疊加的統(tǒng)計變量分布模型；引入時間維度的參數(shù)δ，可用δt替換x，則有式（2）成立：

那么在較大的時間尺度內(nèi)，宏觀的發(fā)動機(jī)受損分布模型為微觀與宏觀模型的疊加。發(fā)動機(jī)受損的邊際分布為

式中：B（ν，δ）為Beta函數(shù)。因此，發(fā)動機(jī)受損分布服從Beta分布模型。

2）突變點檢測

航空發(fā)動機(jī)在服役期間，RUL逐漸縮短，通過突變點搜索模型根據(jù)時間序列顯著性變化可以獲得發(fā)動機(jī)潛在的性能突變位置，從而完成發(fā)動機(jī)衰退過程的階段性劃分。為了確定發(fā)動機(jī)突變點的準(zhǔn)確位置，本文引入了真假突變點的概念：真突變點是通過突變點搜索與顯著性分析，確認(rèn)該時刻發(fā)動機(jī)性能相較前一時刻確實發(fā)生較大變化且滿足統(tǒng)計顯著性要求的突變點。假突變點是經(jīng)分析，不滿足顯著性要求的突變點，這類假突變點將被剔除。該模型包括合并偏差SD（i）、t檢驗的統(tǒng)計量T（i）、統(tǒng)計顯著性P（Tmax）3部分的計算。

合并偏差SD（i）的計算為

t檢驗的統(tǒng)計量T（i）的計算為

式中：n1和n2分別為分割點i左邊部分和右邊部分的點的總數(shù)；u1（i）和u2（i）分別為分割點i左邊部分和右邊部分的均值；S1（i）和S2（i）分別為i點前半部分和后半部分的標(biāo)準(zhǔn)偏差；T（i）分別為用于量化分割點i兩側(cè)差異程度的t檢驗統(tǒng)計值。

統(tǒng)計顯著性P（Tmax）的計算為

式中：通過Monte-Carlo方法得η=4.19 ln m-11.54，m 為時間序列長度；ν=0.4，δ=m-2；Bx（a，b）為不完全Beta函數(shù)。通常情況下，L0≥25，P可取0.5～0.95，L0為最小分割長度。

3）UKF

UKF使用線性Kalman濾波的框架，對協(xié)方差預(yù)測矩陣使用無跡變換（UT）來處理均值和協(xié)方差的非線性傳遞問題［16］。UT變換即按照某一規(guī)則選取一定數(shù)目的樣本點去近似一個正態(tài)分布進(jìn)行UT變換，獲得與原樣本點具有相同均值與方差的Sigma點［17］，從而有效克服線性Kalman濾波線性化誤差較大的缺點。

針對某一非線性系統(tǒng)：

式中：Xk為k時刻被估計狀態(tài)矩陣；Zk為k時刻被估計觀測序列；Г為系統(tǒng)噪聲驅(qū)動矩陣；W 為過程白噪聲矩陣，具有協(xié)方差矩陣Q；V為觀測噪聲矩陣，具有協(xié)方差矩陣R。

該模型包括建立Sigma點矩陣、時間更新、觀測更新3部分。

Sigma點矩陣及對應(yīng)權(quán)值的計算［17］：

式中：n為隨機(jī)變量維數(shù)；chol（Pk）為Pk的Cholesky分解；｛chol（Pk）｝T1，n為Pk的Cholesky分解所獲得的下三角矩陣的第1～n列；Xi（i=0，1，…，2n）被稱為Sigma點集；λ=α2（n+k）-n為一個縮放的尺度參數(shù)，α的取值范圍一般為10-4≤α≤1，k的取值一般設(shè)置為0；采樣點的權(quán)值ω上標(biāo)m為均值，c為協(xié)方差，下標(biāo)為第i個采樣點；β包含著x的分布信息，當(dāng)x符合正態(tài)分布時，β=2。

時間更新的計算為

式中：Xk｜k-1=f（Xk-1）為Xk-1的每一列向量通過f（x）變換得到Xk｜k-1的每一個列向量；（Xk｜k-1）i為Xk｜k-1的第i列；Zk｜k-1=h（Xk｜k-1）與Xk｜k-1=f（Xk-1）的變換相似。

觀測更新的計算為

4）非線性預(yù)測

考慮到多項式函數(shù)具有逼近任何非線性函數(shù)的能力，濾波后多階段的時變參數(shù)通過非線性擬合的方法對發(fā)動機(jī)RUL進(jìn)行預(yù)測。

1.2 多階段航空發(fā)動機(jī)RUL預(yù)測過程的收斂性

通過多階段航空發(fā)動機(jī)RUL預(yù)測模型準(zhǔn)確估計監(jiān)測參數(shù)，以獲得航空發(fā)動機(jī)真實的運行狀態(tài)與RUL，這要求運用模型的過程具有某種收斂性。

定理1對于發(fā)動機(jī)性能退化數(shù)據(jù)，基于超統(tǒng)計的多階段航空發(fā)動機(jī)RUL預(yù)測模型具有收斂性。

證明思路首先證明UKF具有收斂性；其次證明基于超統(tǒng)計的多階段劃分只能劃分出有限個階段；最后證明兩者相融合后模型仍然具有收斂性。

證明

1）在時間序列中，UKF具有收斂性。

2）基于超統(tǒng)計的多階段劃分只能劃分出有限個階段。

通過對突變點搜索得到了突變位置，由于被劃分的觀測序列的長度有限的，那么進(jìn)行分割后將得到一個長度更小的子序列，且子序列的個數(shù)有限。

因此，基于超統(tǒng)計的多階段突變點搜索，只能劃分出有限個階段。

3）融合后的多階段RUL預(yù)測模型具有收斂性。

由此可證明，對于發(fā)動機(jī)性能退化數(shù)據(jù)，基于超統(tǒng)計的多階段航空發(fā)動機(jī)RUL預(yù)測模型具有收斂性。

證畢

2 算 法

2.1 BS-M SF算法

本文提出了基于超統(tǒng)計理論的多階段分割濾波（Multi-stage Segmentation Filtering based on Super statistics theory，BS-MSF）算法，如圖1所示。

BS-MSF算法包括以下6個過程：①對數(shù)據(jù)進(jìn)行預(yù)處理；②進(jìn)行基于超統(tǒng)計的突變點劃分；③進(jìn)行多參數(shù)的融合；④接著進(jìn)行多階段UKF；⑤通過非線性擬合得到退化模型非時變參數(shù)；⑥進(jìn)行RUL預(yù)測。

圖1 BS-MSF算法Fig.1 BS-MSF algorithm

步驟1數(shù)據(jù)預(yù)處理。采用式（14）對數(shù)據(jù)進(jìn)行歸一化無量綱處理：

式中：min為由Sin經(jīng)過無量綱處理后得到的值；Si0為第i個監(jiān)測參數(shù)的第1個循環(huán)時的監(jiān)測值；Sin為第i個監(jiān)測參數(shù)的第n個循環(huán)時的監(jiān)測值；Si1為第i個監(jiān)測參數(shù)的最后1個循環(huán)時的參數(shù)。

步驟2基于超統(tǒng)計理論的突變點劃分。分別計算合并偏差SD（i）、t檢驗的統(tǒng)計量T（i）、T（i）中的最大Tmax的統(tǒng)計顯著性P（Tmax）對監(jiān)測參數(shù)的時間序列進(jìn)行分割。設(shè)置最小分割長度L0和統(tǒng)計顯著性臨界值P，當(dāng)分割后的序列長度小于L0時，不再繼續(xù)分割；如果P（Tmax）≥P，則當(dāng)前點為突變點，繼續(xù)進(jìn)行分割，否則不分割。根據(jù)突變點顯著性的大小進(jìn)行分析，將滿足顯著性要求的突變點視為真突變點，不滿足顯著性要求的突變點視為假突變點。

步驟3多參數(shù)融合。采用式（15）對各監(jiān)測參數(shù)的各階段數(shù)據(jù)進(jìn)行融合，生成多階段健康指標(biāo)H I［17］。

步驟4.1 無跡濾波處理。令與P0分別為原始數(shù)據(jù)的均值與方差，濾波參數(shù)初始化；采用式（9）對健康參數(shù)H I進(jìn)行無跡變換處理；再采用式（11）、式（12）求解第一階段的時間更新參數(shù)與觀測更新參數(shù)。

步驟4.2 利用上一階段的時間更新參數(shù)與觀測更新參數(shù)計算下一階段的相關(guān)參數(shù)，直至計算出全周期的時變參數(shù)｜k。

步驟5非線性擬合。利用式（13）對發(fā)動機(jī)的時變參數(shù)進(jìn)行非線性擬合，獲得各階段非時變參數(shù)。

步驟6利用非時變參數(shù)對測試集中的數(shù)據(jù)進(jìn)行RUL預(yù)測并判斷模型的預(yù)測效果。

2.2 算法優(yōu)點

1）基于超統(tǒng)計理論，通過對突變點的搜索，對發(fā)動機(jī)的衰退過程進(jìn)行階段性的劃分，以“多階段退化”代替“正常與故障”2種狀態(tài)，更符合發(fā)動機(jī)實際退化特點。

2）采用多參數(shù)的信息融合方法表征發(fā)動機(jī)整體的健康狀況，避免單參數(shù)預(yù)測造成結(jié)果的不穩(wěn)定，可充分利用監(jiān)測參數(shù)所包含的信息。

3）狀態(tài)空間方法通過狀態(tài)轉(zhuǎn)化關(guān)系描述揭示系統(tǒng)的內(nèi)在規(guī)律，并且將狀態(tài)變化與時間變量結(jié)合，建立了時間序列的預(yù)測模型。

4）UKF算法基于帶噪聲的觀測數(shù)據(jù)遞推得到預(yù)測模型的時變參數(shù)，因此預(yù)測結(jié)果精度較高且具有較好的預(yù)測演化過程；另外，該算法僅使用發(fā)動機(jī)的觀測數(shù)據(jù)，不需要大量的失效數(shù)據(jù)，因而符合實際工程應(yīng)用。

3 實驗分析

本文選取了美國NASA渦扇發(fā)動機(jī)數(shù)據(jù)中與溫度、壓強、轉(zhuǎn)速有關(guān)的T24，T30，T50，P30，Nf，Nc，T48 7種氣路性能參數(shù)作為實驗數(shù)據(jù)［19］，參數(shù)的具體含義如表1所示。

步驟1本文以美國NASA發(fā)布的數(shù)據(jù)為訓(xùn)練樣本，研究發(fā)動機(jī)性能衰退的階段性特點，采用式（13）對數(shù)據(jù)進(jìn)行歸一化處理，使其不同的監(jiān)測數(shù)據(jù)經(jīng)變換統(tǒng)一到同一區(qū)間，結(jié)果如圖2所示。

步驟2基于超統(tǒng)計理論的多階段分割算法進(jìn)行階段劃分。在訓(xùn)練樣本集中，本文采用L0的取值為25，超參數(shù)q的取值為0.95進(jìn)行階段分割，例如對T24的階段劃分，結(jié)果如圖3所示，其中與縱軸平行的線段為分割點，陰影區(qū)域為嚴(yán)重突變區(qū)域。

例如，參數(shù)T24初步劃分后的分割點分別為：141、192、219、252、270，相應(yīng)突變點對應(yīng)的顯著性如表2所示。

表1 用于發(fā)動機(jī)RUL預(yù)測的重要監(jiān)測參數(shù)Table 1 M ajor m onitoring param eters for p rediction of engine RUL

圖2 歸一化的監(jiān)測序列Fig.2 Normalized monitoring sequence

將該特征參數(shù)的各分割點根據(jù)其顯著性由高到低進(jìn)行排序，對應(yīng)顯著性水平最高的4個突變點作為該特征參數(shù)的最終突變點。因此，參數(shù)T24劃分后的真分割點為：141、192、219、252循環(huán)時，而第270循環(huán)時為假突變點被剔除。同理得7個特征參數(shù)的多階段分割結(jié)果如表3所示。

圖3 T24的突變點與退化量Fig.3 Mutation point and degradation quantity of T24

表2 T24突變點的顯著性Table 2 Significance of m utation point of T24

表3 各特征參數(shù)的對應(yīng)突變點Tab le 3 Corresponding m u tation poin t of each characteristic param eter

可見，7組特征參數(shù)的4個突變點沒有嚴(yán)格在同一個點上，但也比較接近，對整體的突變點位置求解采用計算平均值的方法，最終確定訓(xùn)練序列中4個突變點的具體位置在：第121循環(huán)時、第188循環(huán)時、第221循環(huán)時、第256循環(huán)時。

步驟3采用式（15）將7個監(jiān)測參數(shù)的數(shù)據(jù)進(jìn)行融合。對5個階段的7個參數(shù)進(jìn)行融合后的綜合健康指標(biāo)如圖4所示。

圖4 多階段非線性融合H I序列Fig.4 Multi-stage nonlinear fusion of H I sequence

由圖5（e）可見，通過多階段的UKF預(yù)測性能比卡爾曼濾波的效果更優(yōu)，具有更小的絕對誤差，能更好地模擬發(fā)動機(jī)真實的退化過程。

步驟5非線性擬合。經(jīng)過3種方式濾波的仿真計算，本文通過非線性擬合的方式，對發(fā)動機(jī)壽命衰退進(jìn)行建模。

將不同的非時變參數(shù)作對比，通過單階段KF濾波擬合狀況最優(yōu)概率的預(yù)測模型為

（a0，a1，a2）=（-1.857，0.02257，-6.164×10-5）

通過多階段KF濾波最終確定擬合狀況最優(yōu)概率的預(yù)測模型：

第4階段：

第5階段：

通過多階段UKF濾波最終確定擬合狀況最優(yōu)概率的預(yù)測模型：

第4階段：

第5階段：

壽命預(yù)測演化過程如圖6所示，其擬合誤差結(jié)果如表4所示。

可見，隨著飛行循環(huán)時的增加，對于發(fā)動機(jī)RUL的擬合越來越準(zhǔn)確。

圖5 濾波結(jié)果比較Fig.5 Comparison of filtering results

對于利用該方法對發(fā)動機(jī)RUL預(yù)測結(jié)果的優(yōu)劣，可以用相對誤差、誤差平方和、均方根誤差進(jìn)行評價［20］，本文將測試樣本集中數(shù)據(jù)通過相同的方法進(jìn)行計算處理，并對以上3種RUL預(yù)測模型做了多次實驗。

圖6 壽命演化過程Fig.6 Life evolution process

表4 擬合誤差Table 4 Fitting errors

步驟6壽命預(yù)測。實驗中絕對偏差與相對誤差的定義為

式中：HIestimate為預(yù)測壽命；HIreal為實際壽命值。

多次實驗的預(yù)測誤差結(jié)果如表5所示。不難發(fā)現(xiàn)，單階段KF的相對誤差比多階段KF的誤差要高得多，而多階段UKF的誤差均值在大多數(shù)情況下是比多階段KF更小的，這證明了UKF的狀態(tài)估計準(zhǔn)確性在多數(shù)情況下優(yōu)于KF的估計，但也并不絕對，如第2次實驗中，多階段KF的偏差均值比多階段UKF更小，也就是說，UKF在處理非線性問題中從概率統(tǒng)計的意義上優(yōu)于KF。

表5 預(yù)測誤差Tab le 5 Prediction er rors

4 結(jié) 論

本文提出了一種多階段航空發(fā)動機(jī)RUL預(yù)測模型與相應(yīng)的BS-MSF預(yù)測算法。通過實驗對本文模型與算法進(jìn)行驗證，結(jié)果證明了該方法的合理性。本文提出的理論模型與算法主要貢獻(xiàn)如下：

1）采用超統(tǒng)計分割算法對航空發(fā)動機(jī)的RUL階段進(jìn)行劃分，在確定了發(fā)動機(jī)退化階段的同時，用退化階段客觀描述發(fā)動機(jī)RUL的長短。

2）UKF濾波算法克服了KF的系統(tǒng)初始值不確定、線性假設(shè)前提帶來濾波效果下降的缺點，UKF濾波算法降低了5.5%誤差，在退化過程中對發(fā)動機(jī)真實性能的把握具有更好的適應(yīng)性。

因此，基于多階段航空發(fā)動機(jī)RUL預(yù)測模型的BS-MSF算法是一種能對發(fā)動機(jī)狀態(tài)進(jìn)行跟蹤與RUL預(yù)測的有效算法。