考慮三源不確定性的Wiener過程航空發(fā)動機剩余使用壽命預(yù)測

趙廣社，趙春濤

(1.西安工程大學 電子信息學院，陜西 西安 710048；2.西安交通大學 電子與信息學部，陜西 西安 710049)

0 引 言

航空發(fā)動機作為飛機的核心部件，其可靠性和安全性至關(guān)重要。但是，航空發(fā)動機不僅結(jié)構(gòu)復(fù)雜，而且運行環(huán)境復(fù)雜多變，因而航空發(fā)動機的故障具有易發(fā)生、模式多、多故障模式同時存在等特點。預(yù)測與健康管理(prognostics and health management, PHM)是一種降低故障率、保證飛行安全的有效方法[1]，其中RUL預(yù)測是PHM中的關(guān)鍵技術(shù)之一。RUL可以定義為系統(tǒng)或部件可繼續(xù)正常使用的時間長度[2]，即當前時刻與失效時刻之間的時間間隔。RUL預(yù)測通過分析監(jiān)測數(shù)據(jù)或建立合適的性能退化模型對系統(tǒng)或部件的RUL進行預(yù)測。

現(xiàn)有RUL預(yù)測方法分為3類：基于經(jīng)驗的預(yù)測方法[3]、數(shù)據(jù)驅(qū)動的預(yù)測方法[4]和基于物理失效模型的預(yù)測方法[5]?；诮?jīng)驗的預(yù)測方法易理解、易實現(xiàn)[6]，但其預(yù)測精度有限，并且依賴大量歷史數(shù)據(jù)?；谖锢硎Ｐ偷念A(yù)測方法具有預(yù)測精度高的優(yōu)點[7]，但需要構(gòu)建準確的物理失效模型，對航空發(fā)動機這種復(fù)雜系統(tǒng)來說極難實現(xiàn)?；跀?shù)據(jù)驅(qū)動的預(yù)測方法通過挖掘退化數(shù)據(jù)內(nèi)在的信息[8]，以表現(xiàn)系統(tǒng)自身變化及工作環(huán)境影響導致的系統(tǒng)失效，它避免了復(fù)雜的物理失效模型的構(gòu)建過程，且準確率較高。因此，本文主要研究的是基于數(shù)據(jù)驅(qū)動的航空發(fā)動機RUL預(yù)測。

目前，基于數(shù)據(jù)驅(qū)動的RUL預(yù)測是研究的熱門，并且擁有大量成熟的退化模型，如隨機系數(shù)模型[9]、隨機過程模型[10]和濾波模型[11]等。在基于隨機過程的退化模型中，常常采用Wiener過程退化模型。WANG提出了基于EM算法的退化模型未知參數(shù)估計方法[12]。劉君強等通過實驗證明了航空發(fā)動機的退化過程可以劃分為幾個不同的退化階段[13]。黃亮等建立了多階段Wiener過程性能退化模型，并實現(xiàn)航空發(fā)動機的RUL預(yù)測[14]。然而，對于發(fā)動機這種復(fù)雜系統(tǒng)，在退化過程中還普遍存在三源不確定性，即時變不確定性[15]、個體差異不確定性[16]、測量不確定性[17]?，F(xiàn)有關(guān)于三源不確定性的文獻較少。ZHENG等首次將三源不確定性引入到線性Wiener過程退化模型中[18]，然而，該文獻并沒有考慮退化過程中的非線性性，WANG等同時考慮了非線性性和三源不確定性，并且證明了考慮三源不確定性可以有效提高預(yù)測精度[19]。DONG等利用粒子濾波實現(xiàn)了對鋰電池的RUL預(yù)測[20]。王璽等同時考慮三源不確定性和非線性實現(xiàn)對航空發(fā)動機排氣溫度裕度的預(yù)測[21]。然而，上述考慮三源不確定性的研究中，在考慮測量不確定性時都假設(shè)測量誤差服從正態(tài)分布，但正態(tài)分布不適合用來描述極端值(測量誤差較大)出現(xiàn)概率較大的情況。

對于目前航空發(fā)動機RUL預(yù)測存在的問題，本文建立了考慮三源不確定性的非線性Wiener過程性能退化模型，將三源不確定性引入退化模型并且假設(shè)測量誤差服從Logistic分布。提出一種基于最大熵無跡粒子濾波和條件期望最大化算法的參數(shù)估計方法。最后，采用大型渦輪風扇發(fā)動機數(shù)據(jù)集C-MAPSS進行仿真實驗,驗證本文方法的有效性。

1 三源不確定性的剩余壽命預(yù)測方法

1.1 模型描述

性能退化過程{X(t),t>0}被標準Brown運動B(t)驅(qū)動，那么一般的非線性Wiener過程可表示為

(1)

式中：λ為漂移系數(shù)；σB為擴散系數(shù)；ρ(τ，θ)為時間t的非線性函數(shù)；θ為未知參數(shù)；x0為初始退化狀態(tài)，x0=0。如果ρ(τ，θ)=ρ那么該式就變成線性退化模型，可見該退化模型不失一般性。

從實際角度出發(fā)，個體差異不確定性是指同型號同批次的航空發(fā)動機，受到不同的外部工作環(huán)境和內(nèi)部因素的影響，不同個體的性能退化過程也是存在差異的。因此，在退化模型中引入個體差異不確定性是必要的。根據(jù)文獻[22],假設(shè)漂移系數(shù)λ為隨機參數(shù)，來表現(xiàn)個體間差異的不確定性，建立參數(shù)λ的更新機制為

λn=λn-1+γ

(2)

受到傳感器精度、外部運行環(huán)境噪聲等因素的影響，測量誤差是不可避免的。因此，有必要在建立退化模型時引入測量誤差不確定性。為了表示這種測量誤差，令測量過程{Y(t),t>0}為

Y(t)=X(t)+ε

(3)

式中：ε為測量誤差。為了更好地擬合測量誤差極端值出現(xiàn)概率較大情況下的退化過程，假設(shè)ε服從Logistic分布。

Logistic分布與正態(tài)分布形狀相似，其概率密度函數(shù)為

(4)

式中：v為位置參數(shù)；s為尺度參數(shù)。由式(4)可以算出Logistic分布的超額峰度系數(shù)為1.2。通過超額峰度系數(shù)可知Logistic分布具有尖峰、厚尾的特征，因此Logistic分布比正態(tài)分布適合描述測量誤差極端值出現(xiàn)概率較大的情況，Logistic分布也就常被用來在統(tǒng)計推斷中增加推斷的可靠性。同時，假設(shè)ε、λ和B(t)是相互獨立的。

可以得出性能退化數(shù)據(jù)X(x1,x2,…,xn)的考慮三源不確定性的非線性Wiener過程性能退化模型的狀態(tài)空間模型為

(5)

至此，已經(jīng)在非線性退化模型中引入了三源不確定性。將性能退化狀態(tài)xn和隨機參數(shù)λn作為“隱含”狀態(tài)，進一步整理式(5)可得退化模型的狀態(tài)空間模型：

(6)

1.2 關(guān)鍵檢測參數(shù)選取

實際上，每一個監(jiān)測參數(shù)對退化過程造成的影響不同，為了進行準確的RUL預(yù)測，要篩選出對退化過程影響大的監(jiān)測參數(shù)[23]。航空發(fā)動機的退化過程是必然發(fā)生且不可逆轉(zhuǎn)的，對發(fā)動機退化過程貢獻大的參數(shù)應(yīng)具有單調(diào)性的特點。文獻[24-25]采用定性分析監(jiān)測參數(shù)的單調(diào)性來篩選關(guān)鍵參數(shù)，但該方法只是定性分析，并沒有對參數(shù)的變化趨勢進行量化，不能區(qū)分不同參數(shù)的變化程度。因此，采用Spearman系數(shù)計算各個監(jiān)測參數(shù)的單調(diào)趨勢，可表示為

(7)

由式(7)可以看出，Spearman系數(shù)的取值范圍為[-1,1]。當參數(shù)值是常數(shù)或者隨時間變化隨機變化時，那么Spearman系數(shù)值為0；參數(shù)值隨時間變化增加時，Spearman系數(shù)值大于0，參數(shù)值隨時間變化的單調(diào)增加趨勢越，Spearman系數(shù)值越接近于1，當參數(shù)值隨時間完全單調(diào)增加時，Spearman系數(shù)值為1 ；參數(shù)值隨時間變化的單調(diào)減少趨勢越明顯，Spearman系數(shù)值越接近于-1，當參數(shù)值隨時間完全單調(diào)減少時，Spearman系數(shù)值為-1。

基于上述分析，為了使選取的監(jiān)測參數(shù)擁有明顯的單調(diào)性同時能夠盡量全面地表現(xiàn)航空發(fā)動機的退化過程，選取趨勢值超過閾值0.5的參數(shù)進行后續(xù)的航空發(fā)動機RUL預(yù)測。

1.3 最大熵無跡粒子濾波算法

用監(jiān)測數(shù)據(jù)Y0:n估計的zn的期望和協(xié)方差分別為

(8)

一步前向的期望和方差為

(9)

根據(jù)上述假設(shè)，可以用最大熵無跡粒子濾波算法聯(lián)合估計性能退化正態(tài)和隨機參數(shù)組成的zn。

根據(jù)式(6)，基于監(jiān)測數(shù)據(jù)Y0:n，RUL的概率密度函數(shù)可表示為

(10)

式中：

2 實例驗證

2.1 數(shù)據(jù)來源

為驗證所提方法的有效性，在C-MAPSS數(shù)據(jù)集上進行航空發(fā)動機RUL預(yù)測實驗。該數(shù)據(jù)集是美國航空航天局預(yù)測中心提供的一組C-MAPSS工具模擬產(chǎn)生的由正常運行到失效的過程數(shù)據(jù)集。該數(shù)據(jù)集包含4組不同工作環(huán)境與不同故障模式下的數(shù)據(jù)集，包括引擎單元號、時間戳、3種配置以及 21 個傳感器讀取的數(shù)據(jù)，其中21個傳感器數(shù)據(jù)和3個工作環(huán)境配置參數(shù)用來記錄發(fā)動機的整個性能退化過程。每個數(shù)據(jù)集包括訓練集、測試集和RUL預(yù)測子集。第1組訓練集數(shù)據(jù)是在同一種運行環(huán)境下發(fā)生的單一故障的集合，包含100個發(fā)動機的監(jiān)測樣本，本文選取單一故障模式的第1組數(shù)據(jù)進行后續(xù)實驗。

2.2 關(guān)鍵參數(shù)選取

根據(jù)1.2的分析，選擇FD001中監(jiān)測數(shù)據(jù)樣本進行實驗，計算Spearman系數(shù)得到如圖1所示的參數(shù)變化趨勢結(jié)果圖。

圖 1 參數(shù)變化趨勢Fig.1 Parameter variation trend

從圖1可以看出，數(shù)據(jù)集里的21個傳感器參數(shù)有常數(shù)(Spearman系數(shù)為0)、單調(diào)增加(Spearman系數(shù)為正)與單調(diào)減少(Spearman系數(shù)為負)3種情況。為便于后續(xù)計算，選取擁有單調(diào)遞增趨勢的參數(shù)，但本文方法不會受到監(jiān)測數(shù)據(jù)單調(diào)趨勢的影響。符合單調(diào)遞增的參數(shù)序號為2、3、4、8、11、13、15、17，其趨勢值如表1所示。

2.3 航空發(fā)動機RUL預(yù)測

由航空發(fā)動機的歷史性能退化數(shù)據(jù)建立考慮三源不確定性的航空發(fā)動機性能退化模型。首先，利用Spearman系數(shù)對21個監(jiān)測參數(shù)進行篩選。其次，在確定了關(guān)鍵參數(shù)的基礎(chǔ)上，利用最大熵無跡粒子濾波和ECM算法，進行參數(shù)先驗分布的估計。最后，在獲得新的性能退化數(shù)據(jù)后通過模型參數(shù)的在線更新得到參數(shù)后驗分布。RUL預(yù)測的具體步驟如圖2 所示。

圖 2 RUL預(yù)測流程圖Fig.2 Flow chart of remaining useful life time prediction

為驗證所提方法的有效性，選取M1方法進行對比實驗[19]，實驗分為3個方面進行：首先，將本文方法與M1方法分別對航空發(fā)動機性能退化數(shù)據(jù)進行RUL預(yù)測，對比2種預(yù)測方法的預(yù)測精度；然后，使用M1中卡爾曼濾波與ECM算法相結(jié)合的模型參數(shù)估計方法估計模型的未知參數(shù)，分別采用本文模型與M1模型進行航空發(fā)動機RUL預(yù)測；最后，分別用最大熵無跡粒子濾波與ECM算法相結(jié)合的參數(shù)估計方法和M1中卡爾曼濾波與ECM算法相結(jié)合的參數(shù)估計方法對退化模型的未知參數(shù)進行估計，使用本文模型作為航發(fā)動機的退化模型進行航空發(fā)動機RUL預(yù)測對比實驗。

M1方法：同時考慮三源不確定性和非線性性，但其假設(shè)測量誤差服從正態(tài)分布，不能很好地反應(yīng)測量出現(xiàn)極端值的情況。并且，在參數(shù)估計過程中，模型參數(shù)一旦確定，就不會再隨退化數(shù)據(jù)的累積而更新。

為更好地衡量比較結(jié)果，本文選取RUL預(yù)測領(lǐng)域常用的MSE作為指標，該指標可以很好地表現(xiàn)RUL預(yù)測的精度和不確定性，計算公式為

(11)

圖3給出了2種預(yù)測方法在航空發(fā)動機從出廠到失效的所有觀察點的MSE對比結(jié)果。

圖 3 本文方法與M1方法的MSE對比圖Fig.3 MSE comparison between the method in this paper and the M1 method

從圖3可以看出，在25個飛行循環(huán)之前本文預(yù)測方法的誤差明顯高于M1方法。這是因為，由于本文方法退化模型的初始參數(shù)是隨機的，因此在監(jiān)測初期性能退化數(shù)據(jù)較少時，本文方法的MSE高于M1。在第25個飛行循環(huán)與第60個飛行循環(huán)之間，本文預(yù)測方法的MSE值與M1方法的MSE值接近。這是因為，隨著飛行循環(huán)的增加、監(jiān)測數(shù)據(jù)的增多，本文預(yù)方法的模型參數(shù)自適應(yīng)更新使模型參數(shù)更加貼合實際情況。在過了65個飛行循環(huán)之后本文預(yù)測方法的預(yù)測誤差明顯低于M1預(yù)測方法，整體上看，本文預(yù)測方法的MSE值為13.25，M1預(yù)測方法的MSE值為16.12，說明本文預(yù)測方法有更高的預(yù)測精度，驗證了它的有效性。

在此基礎(chǔ)上，為驗證所提模型的有效性，在2種模型上均采用M1方法估計未知參數(shù)，并進行航空發(fā)動機RUL預(yù)測，實驗結(jié)果如圖4所示。

圖 4 方法相同時2種模型的MSE對比圖Fig.4 MSE comparison diagram of the twomodels when the methods are the same

從圖4可以看出，在前10個飛行循環(huán)2種模型的預(yù)測精度接近，但10個飛行循環(huán)之后本文退化模型比M1模型具有更好的預(yù)測效果。應(yīng)注意的是，在75個飛行循環(huán)附近，2種模型的預(yù)測誤差均出現(xiàn)了上升趨勢，這是因為在75個飛行循環(huán)附近，監(jiān)測數(shù)據(jù)出現(xiàn)了異常值，影響了預(yù)測精度。但本文模型的預(yù)測精度高于M1模型的預(yù)測精度，說明本文模型能夠有效地減弱監(jiān)測數(shù)據(jù)異常值對預(yù)測精度的影響。從整體來看，本文模型比M1模型具有更高的預(yù)測精度，說明本文模型可以有效提高預(yù)測準確性。

接下來，為驗證本文參數(shù)估計方法的優(yōu)越性，采用本文模型作為退化模型，參數(shù)估計方法分別采用M1方法和本文方法，實驗結(jié)果如圖5所示。

圖 5 模型相同時2種參數(shù)估計方法的MSE對比圖Fig.5 MSE comparison diagram of the twoparameter estimation methods when the models are the same

從圖5可以看出，在退化模型相同時，使用本文參數(shù)估計方法進行RUL預(yù)測，具有更高的預(yù)測精度，同時，可看出本文參數(shù)估計方法可以提高預(yù)測模型的收斂速度和魯棒性。

3 結(jié) 語

本文提出了一種考慮三源不確定性的Wiener過程航空發(fā)動機RUL預(yù)測方法。該方法將三源不確定性引入退化模型,且假設(shè)測量誤差服從Logistics分布，能對測量誤差較大的情況進行更好的擬合。在進行航空發(fā)動機RUL預(yù)測時，使用基于Spearman系數(shù)的關(guān)鍵參數(shù)篩選方法，篩選出對退化過程貢獻更大的參數(shù)。然后，將基于最大熵無跡粒子濾波與ECM算法結(jié)合估計模型參數(shù)。在C-MAPSS數(shù)據(jù)集上的仿真實驗，證明了該方法可以提高預(yù)測的準確性，可以提高航空發(fā)動機的使用率和可靠性。航空發(fā)動機的退化過程不僅普遍存在三源不確定性和非線性性，還存在多階段的特點，在下一步研究中，應(yīng)考慮航空發(fā)動機退化過程的多階段性。