基于非線性Wiener過程的航空發動機性能衰減建模與剩余壽命預測

王璽，胡昌華,*，任子強，熊薇

1. 火箭軍工程大學 導彈工程學院，西安 710025 2. 國防大學 聯合勤務學院，北京 100089

航空發動機作為飛機的“心臟”，不僅為飛機提供動力，同時也為航空事業的發展注入強大動力。然而，航空發動機由于長期運行在復雜多變的惡劣環境下，其性能會隨之產生衰減，并在累積到某種量級時，導致發動機失效[1]。目前，航空發動機在航空領域發揮著重要的作用，一旦發生由失效引起的事故，所造成的后果將難以估計。因此，通過監測的航空發動機性能衰減數據，建立相應的性能衰減模型，進而實現航空發動機的剩余壽命(Remanining Useful Life,RUL)預測是近年來熱門的一項關鍵技術[2]。然而，航空發動機的壽命有多種定義，如首翻壽命、大修壽命、總壽命等，同時又受多種因素約束，如性能衰減、關鍵件壽命等。鑒于此，本文主要研究的是在首翻壽命的定義下，基于性能衰減建模的航空發動機剩余壽命預測。

目前，基于性能衰減建模的剩余壽命預測方法已經廣泛地應用于航空發動機[3-5]、導彈[6-8]、電機[9]等。例如，黃亮等[3]基于多階段Wiener過程建立了航空發動機的性能衰減模型，通過融合發動機的歷史性能衰減數據和實時監測數據，實現對發動機的剩余壽命預測。王浩偉等[6]針對導彈多個部件間的競爭退化，綜合考慮了退化失效和突發失效2種競爭模式，提高了對導彈剩余壽命預測的精度。袁慶洋等[9]針對電機的退化特性，建立了多段Wiener過程模型，實現了電機的剩余壽命預測。

然而，對于航空發動機這類復雜的機械系統，在性能衰減過程中普遍存在著三源不確定性，即時變不確定性[10]、個體差異不確定性[11]、測量不確定性[12]。時變不確定性是指性能衰減過程隨著時間變化的固有不確定性；個體差異不確定性是指同類型設備之間的性能衰減差異，即同類型設備的不同個體在性能衰減過程中雖然具有相似的衰減路徑，但實際的衰減率是不相同的；測量不確定性是指測量的性能衰減數據與實際的性能衰減狀態之間的差異。在實際中，由于受噪聲、干擾、儀器精度等影響，想要對性能衰減狀態的精確測量幾乎是不可能的，因此測量中難免會引入測量誤差，造成測量的不確定性。

目前，有關同時考慮三源不確定性的研究較少。為此，文獻[13]基于線性的Wiener過程建立了同時考慮三源不確定性的性能衰減模型，并使得預測的剩余壽命分布可以隨著性能衰減數據在線更新。然而，對于航空發動機而言，在性能衰減過程中不僅具有三源不確定性，還具有強非線性。鑒于此，文獻[14]在非線性下考慮了三源不確定性，并且證明了同時考慮三源不確定性可以有效地提高剩余壽命預測的精度。文獻[15]在此基礎上，利用粒子濾波實現了對鋰電池的剩余壽命預測。但是，以上有關考慮三源不確定性的工作中存在一個共性問題，即假定有多組同類型設備的歷史性能衰減數據來離線地辨識模型初始參數，并且這些參數一旦確定，就不再隨著新的性能衰減數據進行更新。因此，剩余壽命預測的精度受限于數據量，并且預測的結果在很大程度上依賴于選擇的模型初始參數。然而，針對新研發的航空發動機，往往是缺乏此類相關的歷史性能衰減數據和先驗信息，因此，以上方法具有一定的局限性，且目前相關的研究較少。

對于目前航空發動機剩余壽命預測存在的問題，本文建立了同時考慮三源不確定性和非線性的性能衰減模型，并在首達時間下推導出剩余壽命的分布。然后，針對新研發的航空發動機，提出了一種基于Kalman濾波和條件期望最大化(Expectation Conditional Maximization,ECM)算法的參數估計方法，以克服缺乏歷史性能衰減數據和先驗信息的問題，并當新的監測數據可用時，能夠自適應地估計模型參數，進而更新剩余壽命的分布，實現了對新研發航空發動機剩余壽命的在線預測。最后，通過實例驗證了本文方法可以有效地提高預測的精度并降低預測的不確定性，且更加適用于新研發的航空發動機。

1 模型描述

針對航空發動機在性能衰減過程中存在的非線性和不確定性，本文采用一類一般性的非線性Wiener過程來建模航空發動機的性能衰減過程。具體地，性能衰減過程{X(t),t≥0}由標準Brown運動B(t)驅動，可以表示為

(1)

Y(t)=X(t)+ε

(2)

在以上框架下，本文的主要目標是針對缺乏歷史數據和先驗信息的新研發航空發動機，利用實時監測的性能衰減數據來自適應地預測航空發動機的剩余壽命，并且在獲得新的性能衰減數據后實現對剩余壽命的在線更新。具體地，首先假定獲得的性能衰減數據是在離散時間點0=t0

(3)

2 剩余壽命自適應預測

考慮到表征設備健康水平的潛在性能衰減狀態{X(t),t≥0}一旦首次達到設定的失效閾值ω(一般由工業標準設定)時，則認為設備失效。因此，基于首達時間的概念[1-5]，設備在tn時刻的剩余壽命Ln定義為

Ln=inf{ln>0:X(ln+tn)≥ω}

(4)

相應地，剩余壽命Ln的概率密度函數(PDF)為fLn|Y0:n(ln|Y0:n)。因此，預測航空發動機剩余壽命的關鍵是基于監測的性能衰減數據Y0:n={y0,y1,…,yn}，求解fLn|Y0:n(ln|Y0:n)，并在獲得新的性能衰減數據后，實現剩余壽命分布的在線更新。

為了實現剩余壽命的自適應預測，首先將式(3)進一步改寫為

(5)

(6)

式中：δn～N(0,φn)，且有

定義基于監測數據Y0:n估計的zn的期望和協方差分別為

(7)

同樣地，定義一步前向的期望和協方差為

(8)

狀態估計

(9)

協方差更新

Pn|n=Pn|n-1-KnCPn|n-1

(10)

濾波初值

(11)

根據建立的狀態空間模型式(6)，基于監測的性能衰減數據Y0:n，計算剩余壽命的PDF為

fLn|Y0:n(ln|Y0:n)?

(12)

3 基于ECM算法的模型參數估計

(13)

根據極大似然估計(MLE)理論，性能衰減數據Y0:n關于模型參數向量Θ的對數似然函數為

(14)

(15)

然而，由于將潛在性能衰減狀態xn和隨機參數λn視作“隱含”狀態，因此無法直接最大化式(15)。為解決這一問題，自然會采用期望最大化(EM)算法進行參數估計，通過最大化聯合似然函數p(zn,Y0:n|Θ)去逼近模型參數的MLE。EM算法的估計過程可以分為以下2步。

E步：

(16)

M步：

(17)

首先，為了計算對數聯合似然函數，可以通過條件概率公式得到

2lnp(zn,Y0:n|Θ)=

2[lnp(Y0:n|zn,Θ)+lnp(zn|Θ)]=

(18)

然后，通過E步對式(18)求條件期望，有

(19)

式中：

(20)

顯然，為了推導出E步結果，需要先計算式(20)。這里，本文采用Rauch-Tung-Striebel(RTS)平滑濾波算法[19]計算Ω、Ξ、Φ：

平滑濾波

(21)

協方差初值

Mn|n=(I2×2-KnC)AnPn-1|n-1

(22)

后向迭代

(23)

式中：Mi|n=cov(zi,zi-1|Y0:n)。基于RTS算法，條件期望Ω、Ξ、Φ可以計算為

(24)

(25)

(26)

(27)

(28)

CM2：為了基于性能衰減數據Y0:n估計參數θ，首先將式(25)和式(28)一起代入到式(19)中，便可以通過最大化來估計參數即

(29)

由此可見，與傳統的EM算法相比，本文的ECM算法通過CM1和CM2兩步，可以將原本需要優化的3個參數簡化為一個參數θ，降低了計算復雜度，提高了在線性能。此外，本文提出的參數估計方法可以在監測到一個新的性能衰減數據時，實現對模型參數的自適應估計和在線更新，克服了新研發航空發動機歷史數據和先驗信息不足的限制，這對于工程實際具有重要意義。

4 實例驗證

根據相關文獻可知，航空發動機的性能衰減和故障通常表現在排氣溫度裕度(EGTM)的衰減[22-23]。因此，EGTM是表征航空發動機健康狀態的一個重要指標。EGTM的定義為航空發動機出廠前測定的排氣溫度(EGT)限制值與全推力(或全功率)起飛時實際的EGT之間的差值。隨著航空發動機工作時間的增加，EGT會逐漸升高，反映到EGTM就會逐漸衰退，因此準確預測EGTM對提高航空發動機的可靠性至關重要。

實際中，由于外部環境和內部因素的隨機性，航空發動機在某一時刻的性能衰減量也是隨機的，而且同類型航空發動機EGTM的衰減軌跡雖然具有相似性，但實際上不同個體之間仍然存在著差異性。因此，監測的EGTM數據不僅會受到時變不確定性的影響，還會受到個體差異不確定性的影響。此外，EGT數據來源于航空發動機上安裝的傳感器，經過空地數據鏈、報文解析程序等步驟進入到性能監測系統，在這整個過程中不可避免地會引入測量誤差，造成測量的不確定性。根據EGTM的定義可知，EGTM是基于EGT計算的，因此EGTM數據中必然存在測量誤差。通過上述分析發現，監測的EGTM數據中存在著三源不確定性共存的現象，因此同時考慮三源不確定性是必要的，也是更加符合實際需求的。下面以某型號航空發動機在6年內監測的EGTM數據來驗證本文方法的有效性，收集的性能衰減數據如圖1所示。根據工程經驗，設定EGTM的失效閾值ω=0 ℃，并且考慮以下冪函數形式的性能衰減模型：

圖1 航空發動機的性能衰減路徑

X(t)=x0+λ·tθ+σBB(t)

(30)

這里，根據EGTM數據的特點，選擇式(1)中非線性函數的形式為ρ(t;θ)=θ·tθ-1。到目前為止，此類模型已經在航空發動機、軸承、電池等復雜系統的性能衰減建模中得到了廣泛的應用[3,12,17]。

如圖1所示，監測間隔為100次循環，在監測到第4 000次循環時，EGTM首次達到失效閾值ω，期間共監測了41次。顯然，隨著循環次數的增加，航空發動機的性能衰減路徑整體呈下降趨勢。然而，本文所提出的剩余壽命預測方法是基于性能衰減路徑呈上升趨勢的情況下推導出來的，因此需要對原始性能衰減數據進行轉換。這里，本文利用EGTM的初始值減去其他所有的性能衰減數據，相應地，失效閾值轉換為初始值減去ω。下面，基于轉換后的性能衰減數據來驗證本文方法的有效性和優越性。隨著性能衰減數據的積累，采用所提出的參數估計方法可以對模型參數Θ在每一個監測點進行自適應地估計和更新，結果如圖2所示。

圖2 模型參數的自適應估計過程

圖2的結果表明，模型參數可以隨著性能衰減數據的積累很快地收斂，并且在獲得一個新的性能衰減數據時，模型參數和剩余壽命的PDF都可以進行自適應地估計和更新。因此，克服了歷史數據和先驗信息不足的限制，能夠適用于新研發航空發動機的剩余壽命預測。

為了進一步證明本文方法的優越性，選擇以下2種方法進行對比研究：

圖3給出了3種方法在最后6個監測點預測的剩余壽命的PDF曲線，從圖中可以看出，3種方法得到的PDF曲線均包含了實際的剩余壽命，且隨著循環次數的增加逐漸變窄，這意味著預測的剩余壽命的不確定性越來越低。此外，由于本文方法在性能衰減建模中同時考慮了三源不確定性，并且在每一個監測點可以對模型參數進行自適應地估計和在線更新，因此本文方法得到的PDF曲線更加窄而尖銳，說明本文方法在剩余壽命預測方面具有更高的精度和更好的性能。

圖3 三種方法的比較結果

為了進一步量化比較結果，本文采用可靠性領域中常用的均方誤差(MSE)指標，該指標同時考慮了剩余壽命預測的精度和不確定性，在tn時刻的定義式為

(31)

圖4 所有監測點的MSE

圖5 4個不同監測點的箱形圖(點劃線表示實際剩余壽命)

圖4給出了3種方法在所有監測點的MSE比較結果。由于本文方法和M2方法均不依賴于歷史數據和先驗信息，且模型的初始參數是隨機給定的，因此在監測初期性能衰減數據較少時，剩余壽命預測的MSE高于M1，但隨著性能衰減數據的積累，本文方法的MSE明顯低于M1和M2。為進一步驗證，圖5給出了4個不同監測點的箱形圖。其中，箱子的大小表示預測的剩余壽命的不確定性，若箱子越大，則預測的不確定性越大，精度越低；點劃線表示當前監測點的實際剩余壽命；箱子中的實線表示預測的剩余壽命的期望，實線越接近點劃線，說明預測的精度越高，結果越準確。從圖5(a)中可以看出，此時由于性能衰減數據較少，M1優于本文方法和M2，但是隨著性能衰減數據的積累，本文方法的箱子更小且更加接近實際的剩余壽命，如圖5(b)～圖5(d)所示。此外，引入相對誤差(RE)指標來進一步量化剩余壽命預測的精度，在壽命的20%、45%、70%和95%分位點上分別給出了3種方法RE的比較結果，如表1所示。

通過表1的結果可以看出，本文方法可以有效減少剩余壽命預測的RE，進而提高了剩余壽命預測的精度，尤其在壽命的95%分位點處，剩余壽命預測的RE僅為1.68%，明顯小于其他2種方法。以上的實驗結果表明，本文方法具有更好的預測能力，剩余壽命預測的精度更高，且不依賴于歷史數據和先驗信息，更加適用于新研發的航空發動機。

表1 不同壽命分位點相對誤差的比較結果

5 結 論

1) 本文在三源不確定性下建立了非線性的性能衰減模型，不僅克服了現有方法的潛在假設，而且可以較好地建模航空發動機的性能衰減軌跡，并能夠將三源不確定性的影響納入到預測的剩余壽命分布中。

2) 利用Kalman濾波技術可以對潛在的性能衰減狀態進行實時估計，并基于首達時間的概念，得到了同時考慮三源不確定性的剩余壽命分布。

3) 在獲得一個新的性能衰減數據時，基于RTS平滑和ECM算法對模型參數進行自適應地估計和在線更新，克服了新研發航空發動機缺乏歷史數據和先驗信息的問題，降低了參數估計的不確定性，顯著提高了剩余壽命預測的精度，尤其在壽命的95%分位點處，剩余壽命預測的RE僅為1.68%。