基于超統計的多階段航空發動機剩余壽命預測

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

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

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

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

在剩余壽命（Remaining Useful Life，RUL）預測模型研究中，主要包括基于數據驅動的方法、基于物理模型的方法和融合的方法。然而，航空發動機因設計結構復雜導致其物理模型難以確定，基于物理模型的方法和融合的方法，目前的研究較少。

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

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

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

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

1 航空發動機RUL預測模型

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

1.1 模型組成

1）超統計理論

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

從微觀角度建模。運用統計學方法對發動機指標變化情況進行分析，一般情況下，在環境參數ν一定的情況下，發動機健康狀況異常的條件概率密度為

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

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

那么在較大的時間尺度內，宏觀的發動機受損分布模型為微觀與宏觀模型的疊加。發動機受損的邊際分布為

式中：B（ν，δ）為Beta函數。因此，發動機受損分布服從Beta分布模型。

2）突變點檢測

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

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

t檢驗的統計量T（i）的計算為

式中：n1和n2分別為分割點i左邊部分和右邊部分的點的總數；u1（i）和u2（i）分別為分割點i左邊部分和右邊部分的均值；S1（i）和S2（i）分別為i點前半部分和后半部分的標準偏差；T（i）分別為用于量化分割點i兩側差異程度的t檢驗統計值。

統計顯著性P（Tmax）的計算為

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

3）UKF

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

針對某一非線性系統：

式中：Xk為k時刻被估計狀態矩陣；Zk為k時刻被估計觀測序列；Г為系統噪聲驅動矩陣；W 為過程白噪聲矩陣，具有協方差矩陣Q；V為觀測噪聲矩陣，具有協方差矩陣R。

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

Sigma點矩陣及對應權值的計算［17］：

式中：n為隨機變量維數；chol（Pk）為Pk的Cholesky分解；｛chol（Pk）｝T1，n為Pk的Cholesky分解所獲得的下三角矩陣的第1～n列；Xi（i=0，1，…，2n）被稱為Sigma點集；λ=α2（n+k）-n為一個縮放的尺度參數，α的取值范圍一般為10-4≤α≤1，k的取值一般設置為0；采樣點的權值ω上標m為均值，c為協方差，下標為第i個采樣點；β包含著x的分布信息，當x符合正態分布時，β=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）非線性預測

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

1.2 多階段航空發動機RUL預測過程的收斂性

通過多階段航空發動機RUL預測模型準確估計監測參數，以獲得航空發動機真實的運行狀態與RUL，這要求運用模型的過程具有某種收斂性。

定理1對于發動機性能退化數據，基于超統計的多階段航空發動機RUL預測模型具有收斂性。

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

證明

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

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

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

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

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

由此可證明，對于發動機性能退化數據，基于超統計的多階段航空發動機RUL預測模型具有收斂性。

證畢

2 算 法

2.1 BS-M SF算法

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

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

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

步驟1數據預處理。采用式（14）對數據進行歸一化無量綱處理：

式中：min為由Sin經過無量綱處理后得到的值；Si0為第i個監測參數的第1個循環時的監測值；Sin為第i個監測參數的第n個循環時的監測值；Si1為第i個監測參數的最后1個循環時的參數。

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

步驟3多參數融合。采用式（15）對各監測參數的各階段數據進行融合，生成多階段健康指標H I［17］。

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

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

步驟5非線性擬合。利用式（13）對發動機的時變參數進行非線性擬合，獲得各階段非時變參數。

步驟6利用非時變參數對測試集中的數據進行RUL預測并判斷模型的預測效果。

2.2 算法優點

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

2）采用多參數的信息融合方法表征發動機整體的健康狀況，避免單參數預測造成結果的不穩定，可充分利用監測參數所包含的信息。

3）狀態空間方法通過狀態轉化關系描述揭示系統的內在規律，并且將狀態變化與時間變量結合，建立了時間序列的預測模型。

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

3 實驗分析

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

步驟1本文以美國NASA發布的數據為訓練樣本，研究發動機性能衰退的階段性特點，采用式（13）對數據進行歸一化處理，使其不同的監測數據經變換統一到同一區間，結果如圖2所示。

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

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

表1 用于發動機RUL預測的重要監測參數Table 1 M ajor m onitoring param eters for p rediction of engine RUL

圖2 歸一化的監測序列Fig.2 Normalized monitoring sequence

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

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

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

表3 各特征參數的對應突變點Tab le 3 Corresponding m u tation poin t of each characteristic param eter

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

步驟3采用式（15）將7個監測參數的數據進行融合。對5個階段的7個參數進行融合后的綜合健康指標如圖4所示。

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

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

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

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

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

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

第4階段：

第5階段：

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

第4階段：

第5階段：

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

可見，隨著飛行循環時的增加，對于發動機RUL的擬合越來越準確。

圖5 濾波結果比較Fig.5 Comparison of filtering results

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

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

表4 擬合誤差Table 4 Fitting errors

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

式中：HIestimate為預測壽命；HIreal為實際壽命值。

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

表5 預測誤差Tab le 5 Prediction er rors

4 結 論

本文提出了一種多階段航空發動機RUL預測模型與相應的BS-MSF預測算法。通過實驗對本文模型與算法進行驗證，結果證明了該方法的合理性。本文提出的理論模型與算法主要貢獻如下：

1）采用超統計分割算法對航空發動機的RUL階段進行劃分，在確定了發動機退化階段的同時，用退化階段客觀描述發動機RUL的長短。

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

因此，基于多階段航空發動機RUL預測模型的BS-MSF算法是一種能對發動機狀態進行跟蹤與RUL預測的有效算法。