基于性能非線性退化統計模型的發動機剩余壽命預測

郭 慶，鄭天翔，李印龍

（中國民航大學航空工程學院，天津 300300）

0 引言

航空發動機是航空器的“心臟”，其健康狀態對保證航空器正常運轉具有重要意義，但是航空發動機由于長時間在復雜多變的惡劣環境下運行，其性能會隨著服役時間延長出現退化現象，并且在退化到某一閾值時，最終需將發動機下廠翻修。為了準確地評估發動機的健康狀態，需對發動機可能發生的故障及時處理，并及時地配合發動機的運行規劃和合適的維護決策，減少災難性事故的發生，減小因事故造成的損失，因此需要通過監測航空發動機性能退化數據，構建相應的退化模型，進而預測發動機的實時剩余壽命。

隨著傳感器技術和預測技術的迅速發展，從視情維修中發展出了1個重要方向——預測維修，而作為實現預測維修的核心技術，壽命預測近年來已經成為中國研究的熱點之一。主流的研究方法可分為基于模型的方法和基于數據驅動的方法。基于模型的方法需要借助其故障模式和失效機理建立模型來評估其健康狀態，典型的預測方法是利用統計方法對發動機壽命數據進行統計分析，進而得到其壽命分布，此外還有隨機損傷傳播模型等；而基于數據驅動的方法不需要建立發動機的性能模型，通過狀態監測數據和歷史故障數據進行預測，診斷方法包括神經網絡方法、支持向量機法、相關向量機法、自回歸（Auto Regressive，AR）模型等。然而在剩余壽命預測領域，設備的性能退化過程多以關于時間的線性函數為主，鮮有研究非線性退化的情況，對此王小林等利用非線性Wiener過程對產品性能退化過程建模，表明非線性的退化過程更能反映某些復雜設備的性能退化軌跡。在航空發動機研制中受材料、加工條件等因素的影響，使得同類型產品的性能退化過程存在差異性，需要在剩余壽命預測中加以考慮。

針對航空發動機性能退化的非線性、不確定性，本文基于傳感器得到的發動機監測數據，構建非線性性能退化統計模型，利用貝葉斯準則對參數進行實時更新，開展發動機實時剩余壽命預測。

1 非線性性能退化建模

1.1 非線性退化數據統計模型

在航空發動機退化過程中，按照時間順序，根據傳感器所測得的退化數據分析系統主要性能的退化過程，并建立退化軌道，退化軌道是系統性能退化的真實寫照。具體地，性能退化過程表示為

式中：X(t)為系統在時間t處的真實退化值；y(t)=(y(t),…,y(t))為q×1的矩陣，其中y(t)(l=1,2,…,q)是t的已知函數；α為q×1的隨機變量，且E(α)=μ，Cov(α)=∑。

式（1）給出的軌道中關于t的函數形式可以是線性的，也可以是非線性的。

由于受測量儀器精度和環境噪聲影響，發動傳感器監測的衰退數據不可避免地存在測量誤差，因此引入隨機誤差ε描述性能退化過程中的測量不確定性

據此對n臺發動機在m個給定時間點的數據統計模型綜合如下

式中：α為第i個發動機的隨機參數向量，用以描述個體間的差異不確定性，可看作是來自某概率分布α～N(μ,∑)的1個樣本，并且與ε相互獨立。

實際上，由于受初始條件、制造工藝、操作和環境等因素的相互影響，即使同類型的航空發動機，不同個體間的性能退化過程也往往存在著差異，且在某一時刻的性能衰退量也是隨機的，這些隨機波動時大時小，時正時負，很難控制，最后綜合地表現在退化曲線和失效時間上。發動機性能退化符合非線性退化過程，故在發動機性能衰退建模中，考慮發動機退化軌道模型中關于t的函數形式采用指數模型

綜上所述，基于退化軌跡的發動機性能退化模型采用滿足于退化帶有隨機效應的非線性退化模型，且為了不失一般性，初始值Y(0)設為零

1.2 非線性退化模型參數估計

對于式（2），若Rank(X)=q，則μ,σ,∑分別有以下無偏估計。

（1）隨機參數α均值μ的最小二乘估計。

2 發動機剩余壽命預測

2.1 貝葉斯更新參數

由于受到內部因素和外部因素的影響，航空發動機性能退化過程存在著不確定性，為了減小不確定性對航空發動機剩余壽命預測結果準確性的影響，對于具體運行的發動機，根據貝葉斯準則，利用其運行過程中實時的監測數據對模型中的隨機參數的概率分布進行更新，使得退化建模的結果更加符合該發動機實際的運行情況。因此，在 時刻，該發動機所有的傳感數據為={(t ),≤t ≤t }，并且令X =(t )，基于的先驗分布和歷史數據，根據Bayesian鏈式規則得到模型參數α的后驗分布

2.2 個體發動機剩余壽命預測

航空發動機的剩余壽命是指從當前時刻開始，隨著服役時間的延長，表征發動機健康水平的性能衰退狀態{ }

(),≥0首次達到紅線值的時間。基于首達時間定義發動機的壽命過于嚴格，當性能退化量首次達到紅線值時，發動機可能僅僅是性能降低，并未失效，但從安全角度上考慮，此時應該停止使用。

因此由退化軌道()首次達到紅線值時可以確定1個失效分布，這是退化模型參數的分布，發動機的使用壽命可以定義為

式中：T 為隨機變量，為第個發動機的壽命；為發動機的紅線值。

通常紅線值是確定值，但在發動機的退化過程中存在隨機性，不同類型、批次的發動機在真實失效時對應的退化量不盡相同。為了減少不確定性，在相同失效模式和操作條件下，發動機的紅線值可以通過使待定紅線值與不同發動機失效時刻對應的退化量的方差最小值來確定

對于固定的X ，與L 的分布依賴于基本軌道參數(α;)。在退化軌道的函數形式是線性的情況下，有可能寫出()的顯式表達式；在函數形式是非線性的情況下，這種顯式表達式很難得到，采用數值仿真方法估計()的表達式是很好的解決方式。

為了避免上述多重積分計算可能帶來的困難，本文采用蒙特卡洛方法來模擬發動機退化過程，蒙特卡洛仿真流程如圖1所示。利用隨機數來描述發動機性能狀態隨運行時間的變化，將發動機性能退化的確定性問題轉化為隨機問題。在發動機運行過程中，通過模擬總體樣本發動機退化行為的規律和個體隨機發生相互作用，對大量的隨機事件追蹤和統計來獲得發動機的性能狀態。

圖1 蒙特卡洛仿真流程

蒙特卡洛方法解題主要步驟為：

（1）描述或構造概率過程。對本身具有隨機性的問題建立正確描述問題的概率統計模型，模型的概率分布可以很好地反映所求問題。

（2）概率分布抽樣。概率模型確定后，根據模型進行大量隨機試驗獲得隨機變量的試驗值。確定隨機變量的產生方法，包括隨機數或偽隨機數的產生及根據概率統計模型對隨機數進行抽樣，根據試驗規則進行大量相互獨立的隨機抽樣試驗。

（3）得到各統計量的估計值。根據概率模型進行模擬試驗后，分析試驗結果的數學特性。

3 實例驗證

本試驗采用C-MAPSS數據集中FD001數據集進行試驗驗證，C-MAPSS數據集是基于NASA出品的渦扇發動機模擬模型C-MAPSS，通過模擬各種操作條件并加入不同程度的衰退故障生成的數據集，廣泛地應用于各種工況下的發動機性能退化試驗。該數據集模擬了發動機在單工況單故障情形下模擬運行產生的退化以及失效數據，共包含了100臺失效的發動機全壽命退化數據，每組監測數據包括21個傳感器的監測數據。

根據工程實踐可知，航空發動機的性能衰退和故障通常表現在排氣溫度裕度值（Engine Gas Temperature Margin，EGTM）的衰減，因此EGTM是表征發動機健康狀態的1個重要指標。本文使用與EGTM成正相關關系的排氣溫度裕度（Exhaust Gas Temperature，EGT）作為表征發動機退化狀態的性能參數。為驗證非線性退化軌道模型在航空發動機剩余壽命預測中的準確性，選擇100臺全壽命周期的前80臺發動機EGT歷史監測數據作為訓練組，后20臺發動機的EGT監測數據用作實時更新測試。

由于初始數據有著極大的隨機特性，并且在每個監測參數數據中都加入了大量隨機噪點，該數據集存在一定未知程度的初始退化，不滿足于建模研究，需要對數據進行預處理（如圖2所示），即在不改變數據的數字特征的情況下，降低數據的數量級。隨機取出1臺發動機對其性能參數采用局部加權平均法進行平滑處理

圖2 數據預處理

式中：為與當前中心點x 關聯的跨度內的相鄰點；d 為沿橫坐標x 到跨度內最遠的相鄰點的距離。

采用參數估計方法得到的參數估計值的結果見表1。

表1 參數估計結果

以表1中的參數估計值為先驗信息，利用式（23）、（24）對退化模型中的隨機參數的均值和方差進行實時更新（其結果如圖3、4所示），實現對發動機健康情況的實時跟蹤。

圖3 α均值更新結果

圖4 α方差更新結果

從圖3、4中可見，退化模型的參數可以隨著衰退數據的積累動態更新，并且每當獲得1條新的衰退數據，模型的參數可以自適應地估計和更新，從而得到更新后的剩余壽命的概率密度函數PDF（PDF，Probability Density Function），彌補了歷史數據和先驗信息的不足，表明了本文所提方法能夠實時跟蹤發動機的及時健康情況。

隨機參數的均值和方差隨著發動機退化數據的增加不斷更新。從圖3中可見，在前100個循環中發動機性能退化趨勢平緩，隨機參數的更新結果變化較小，發動機健康狀態基本不變。隨著運行時間的增加，發動機退化趨勢變大，隨機參數均值出現一定范圍的波動，并且其變化趨勢隨監測數據變化而出現不同程度地波動；而隨機參數方差的更新演化趨勢表明，隨著退化數據的不斷增加，對隨機參數估計值的不確定性逐漸減弱。

蒙特卡洛仿真結果數據就是發動機剩余壽命數據，對其以柯爾莫哥洛夫-斯米爾諾夫（Kolmogorov-Smirnov，K-S）檢驗標準進行分布擬合，擬合結果見表2。

表2 觀測顯著性水平P值

值是統計檢驗的觀測顯著性水平，＞0.05時，不排除當前擬合的分布形式。表2中數據展示了在發動機運行60、120和180循環時，蒙特卡洛對發動機剩余壽命預測仿真結果。顯然指數分布的＜0.05，則拒絕0.05水平的分布，而正態分布、對數正態分布和Gamma分布的＞0.05且相差不大，以發動機在60循環時的剩余壽命預測為例，3種分布擬合如圖5所示。從圖中可見，3種分布趨勢相近，都是數據良好的模型。

圖5 3種分布擬合

3種分布擬合優度值見表3。從表中可見，在分布擬合相近的3種分布形式，正態分布擬合優度值比其他2種分布的更接近仿真數據分布模型得到的結果。K-S檢驗結果如圖6所示。從圖中可見，樣本數據點基本在1條直線上且在95%置信帶內，因此選擇正態分布來描述發動機的剩余壽命預測PDF，使用剩余壽命仿真結果的期望值作為發動機剩余壽命的預測值。

圖6 K-S檢驗結果

表3 擬合優度值

對發動機剩余壽命預測結果如圖7所示。從圖7（a）中可見，對發動機剩余壽命預測的誤差集中在發動機運行早期，在中后期的相對誤差明顯減小；在發動機運行后期的剩余壽命預測精度較發動機早期和中期的明顯提高。

圖7 剩余壽命預測結果

從圖7（b）中可見150循環后發動機剩余壽命預測的PDF曲線，每條PDF曲線均包含了發動機真實剩余壽命，且隨著運行時間的增加預測剩余壽命和真實剩余壽命逐漸靠近（如圖8所示），意味著對發動機剩余壽命預測的不確定性逐漸減弱。此外，在每個監測點模型參數會根據退化數據的增加自適應地估計和在線更新，因此，所得到的PDF曲線逐漸變窄而尖銳，說明了本文方法可以更好地預測發動機的剩余壽命。

從圖8中可見，隨著監測數據的積累，黑線越來越貼近紅線，剩余壽命預測的精度逐漸提高，并且根據圖7的概率密度函數，可以計算發動機剩余壽命預測的均方誤差，如圖9所示。從圖中可見，隨著時間的推移，剩余壽命預測的均方誤差越來越小，剩余壽命的預測精度逐漸提高。

圖8 剩余壽命預測值與真實值對比

圖9 均方誤差

為了進一步證明本文方法的優越性，與一種由標準布朗運動驅動的線性Wiener過程預測發動機剩余壽命的方法進行對比。2種方法在所有監測點的預測剩余壽命與真實剩余壽命的比較結果如圖10所示。從圖中可見，線性方法在發動機運行早期的預測結果并不理想，在運行后期的剩余壽命預測具有一定的參考性；而本文方法在發動機運行早期就可以為判斷發動機健康狀態提供一定的依據，在發動機運行全過程中，預測結果逐漸向真實剩余壽命靠近，具有較強的預測能力，剩余壽命預測精度更高。

圖10 剩余壽命預測結果

4 結論

（1）基于首達時間的概念，利用統計和蒙特卡羅仿真可以基于非線性退化數據對航空發動機進行退化建模，實現了對個體發動機剩余壽命實時預測。

（2）利用C-MAPSS數據集對本文所述方法進行試驗，表明本文所述方法隨著運行時間的推移，對發動機剩余壽命的預測精度逐漸提高。