基于相關向量機的發動機剩余壽命預測

彭鴻博, 蔣雄偉

(中國民航大學航空工程學院，天津300300)

民航發動機的維修與保障對于飛機的安全運行至關重要，基于發動機剩余壽命(remaining useful life, RUL)預測來確定下發期限能夠為發動機工程師制定科學合理的維修方案提供依據。由于民航發動機使用可靠性高、運行周期長，因此利用傳統的統計方法收集失效數據建立統計模型進行壽命預測比較困難[1-2]。并且航空發動機是囊括了機、電、液的復雜系統，失效模式較多，基于失效機理建立物理模型準確預測剩余壽命的方法只能針對發動機的某個部件，不適用于整機[3-4]。

民航發動機在運行過程中必須進行狀態監控，監測參數的變化情況能夠反映發動機的健康狀態，利用性能退化數據挖掘發動機的性能衰退規律，通過統計推斷和機器學習等方法預測剩余壽命是目前研究的熱點。基于狀態監控的剩余壽命預測通常采用的思路是選取單個或多個監測參數作為發動機性能衰退的退化量，通過預測退化量達到失效閾值的時間來實現發動機剩余壽命的預測。馬小駿等[5]采用排氣溫度裕度(exhaust gas temperature margin, EGTM)作為性能衰退的指標，通過分析發動機EGTM的發展趨勢，采用最小二乘支持向量機(least squares-support vector machine, LS-SVM)建立性能可靠性模型來預測發動機的在翼壽命；黃亮等[6]根據發動機的歷史性能監測數據，在考慮發動機性能退化過程非線性和多階段特點的基礎上，建立Wiener退化模型并采用貝葉斯方法更新模型參數，通過預測排氣溫度(exhaust gas temperature, EGT)超過失效閾值的時間來預測剩余壽命。

與使用單性能參數的預測方法相比，基于多個性能參數預測發動機剩余壽命的方法能夠更有效地利用監測數據，預測的準確性也更高。張頡健[7]根據發動機采集的監測參數為多維離散型數據的特點，通過改進過程神經網絡使其能夠滿足離散輸入與多個參數融合的要求，建立了發動機性能狀態預測模型；Nieto等[8]采用支持向量機(support vector machine, SVM)挖掘性能參數與發動機剩余壽命之間的對應關系，建立了直接預測模型，并針對SVM需要優化的模型參數較多的問題采用粒子群算法進行了參數尋優；Babu等[9]將深度卷積神經網絡引入到復雜設備剩余壽命預測領域中，分析其可行性并使用發動機性能退化仿真數據集進行了驗證；Ordóez等[10]利用自回歸滑動平均模型(autoregressive integrated moving average model, ARIMA)對發動機各監測參數分別進行預測，再以多維參數的預測值作為SVM的輸入，以發動機剩余壽命為輸出，建立了預測模型；宋亞等[11]將自編碼神經網絡(autoencoder)和雙向長短期記憶(bidirectional long short-term memory, BLSTM)神經網絡進行整合，結合兩者的優勢提出了一種混合健康狀態預測模型，并將其用于渦扇發動機的剩余壽命預測。

前人研究中主要依靠機器學習算法較強的非線性映射能力直接建立發動機性能參數與剩余壽命的關系，然而預測過程中監測參數的選取對預測結果的影響較大，監測參數篩除較多容易導致信息遺漏，保留過多的監測參數又會使模型誤差累積過大，預測精度降低。為此，在發動機所有監測參數的基礎上，應用核主元分析(kernel principle component analysis, KPCA)算法提取出包含發動機性能衰退信息的低維特征，并進行非線性融合得到描述發動機退化軌跡的健康指數(health index,H)序列，再利用相關向量機(relevance vector machine，RVM)建立回歸模型，通過外推預測H超過失效閾值的時間，間接獲取發動機剩余壽命的預測值。

1 退化特征信息融合

民航發動機在不考慮維修的情況下，其性能衰退過程可近似地視為不可逆的單調退化過程，退化趨勢通過監測的性能參數變化情況表現，但不同參數對發動機退化程度的敏感程度不同，并且各參數之間存在一定的相關性，因此需要對原始監測數據進行處理，將其融合成能夠反映發動機退化趨勢的健康指數。

1.1 基于KPCA的退化特征提取

核主元分析是一種將核函數與主元分析(principle component analysis, PCA)結合的數據處理方法，主要針對信息融合中的非線性問題，其基本原理是通過非線性映射將原始數據投影到高維特征空間，然后在高維空間中執行主元分析，消除原始數據中的冗余信息，提取出能夠描述原始數據變化趨勢的主元序列[12]。民航發動機運行過程中監測的性能參數較多，且各參數之間存在一定的非線性關系，為此采用KPCA對監測數據進行處理。

假設采集到的發動機性能監測數據為X=[x1,x2,…，xN]∈RM×N，其中N為性能參數的個數，xi=[x1,x2,…,xM]T表示第i個監測參數的時間序列，對其進行核主元分析的步驟如下。

(1)首先對監測數據進行標準化預處理，處理后得到矩陣X*，其均值為0，方差為1。

(2)選擇核函數，并求核函數矩陣K。根據監測數據的特點，采用的核函數為高斯徑向基核函數：

K(xi,xj)=exp(-‖xi-xj‖2/σ)

(1)

式(1)中：σ為核函數寬度。采用交叉驗證的方法選取最佳的核函數寬度值。

K′=K-INK-KIN+INKIN

(2)

式(2)中：IN為N維單位矩陣。

(4)求協方差矩陣的特征向量ν和特征值λi，將主對角線上為特征值的對角陣變換成特征值列向量，特征值按降序排列：λ1≥λ2≥…≥λN。

(5)以主元貢獻率為依據提取前p(p≤N)個主元，對應的特征向量可組成特征向量矩陣V=[ν1,ν2,…,νp]。為保證最大程度保留原有信息量的同時又實現特征空間的降維，設定累計主元貢獻率rCS必須大于85%，計算公式為

(3)

式(3)中：λi為前p個主元的特征值；λj為全部特征值。

(6)計算原始數據在特征向量矩陣V上的投影，輸入樣本x在第k個主元νk上的投影為

2,…,p

(4)

式(4)中：Φ(·)為映射函數；αi為特征向量的系數。

得到發動機性能監測數據的核主元特征為X′=[t1,t2,…,tp]，其中tp=[t1,t2,…,tM]T表示第p個主元的時間序列。

1.2 健康指數的計算方法

經過核主元分析處理后得到的主元序列，既保留了原始數據的特征信息，又在維度上實現了約簡，但其不能夠直接反映發動機性能衰退的趨勢，因此在發動機性能監測數據核主元特征的基礎上，將主元序列融合成健康指數來表征系統的健康程度，其衰減的過程對應發動機性能退化的過程[13]。根據發動機性能衰退的特點，采用非線性模型：

y=aeX′BT+c

(5)

式(5)中：a、B=[b1,b2,…,bp]為非線性融合模型的系數；c為常數項；y為發動機的健康指數序列。

該模型可建立核主元特征X′與健康指數y之間的關系，利用該模型來融合多維主元序列可得到一維的發動機健康指數(health index,H)序列。為了求解模型系數，需要構造訓練樣本集，令發動機剛開始使用時前5個飛行循環H=1，因性能衰退需要下發時后5個飛行循環H=0，得到訓練樣本集：

(6)

式(6)中:X′M為原始數據第M個樣本在特征空間上的投影，使用訓練樣本集Ω來訓練模型，可獲得模型系數a、b1，b2，…，bp和常數項c的值。通過非線性模型融合得到的發動機H序列，不僅保留了發動機性能衰退的信息，而且能夠直接反映發動機的性能衰退過程。

2 基于RVM的發動機RUL預測

基于RVM的發動機剩余壽命預測主要包括訓練過程和測試過程兩個階段。訓練階段從發動機歷史失效數據中獲取完整H序列作為輸入，建立RVM模型；測試階段利用訓練獲得的RVM模型對發動機當前的H序列進行外推預測，獲取H預測值，比較H超過失效閾值的時間和預測起始時刻，計算發動機的剩余壽命，預測方法的框架如圖1所示。

圖1 RUL預測框架Fig.1 RUL prediction frame diagram

2.1 RVM模型

相關向量機是一種將Bayesian理論與SVM結合的機器學習算法。相比SVM，RVM的核函數沒有特殊限制，稀疏性更好，較少的超參數降低了核函數的計算量，且具有概率形式的輸出，可用于時間序列的分析與預測[14]。因此，將RVM用于民航發動機RUL預測。

(7)

式(7)中：n為樣本數；ω=(ω0,ω1,…,ωn)T為模型參數；K(x,xi)為非線性核函數。

RVM模型的預測過程如下。

(1)設定超參數α和噪聲方差σ2的初始值。α為由N+1個超參數αi組成的向量，α=[α0,α1,…,αN]T，其中αi為權值ω對應的超參數。

(2)采用高斯核函數，根據當前的α和σ2計算模型參數ω的后驗分布：

p(ω|y,α,σ2)=N(μ,Σ)

(8)

式(8)中：y=(y1,y2,…,yn)T；α=(α1,α2,…,αn)；N(·)為多變量高斯分布；μ=σ-2ΣΦTy為均值；Σ=(σ-2ΦTΦ+A)-1為協方差，其中A=diag(α1,α2,…,αn)，Φ為n(n+1)維的設計矩陣，Φ=[φ(x1),φ(xn),…,φ(xn)]T，其中φ(xi)=[1,K(xi,x1),…,K(xi,xn)]T；均值μ作為發動機H的預測值；協方差Σ表示模型預測的不確定性。

(9)

(10)

式中：γi=1-αiΣii，Σii為協方差矩陣中第i項對角元素；μi為第i個后驗期望值。

μtest=μTφ(xtest)

(11)

(12)

(5)根據預測均值和方差可以計算預測值的概率密度p(ytest)和分布函數F(ytest)，得到預測值的置信區間。

(6)設定失效閾值作為預測終止條件，截取待測樣本預測起始時到H序列超過失效閾值時的飛行循環數作為發動機的剩余壽命預測值。

2.2 模型評價指標

任何方法應用到實際工程中都需要具備一定的實用價值，使用本文方法對待測發動機進行剩余壽命預測后，可用一個或多個評價指標來顯示該方法的預測效果。通常采用的評價方法是對多個設備進行壽命預測，選擇平均百分比誤差(mean absolute percentage error, MAPE)作為評價預測準確性的指標，其計算公式為

(13)

式(13)中：K為待測樣本數量，Rpre,d、Ract,d分別表示第d個待測樣本剩余壽命的預測值和真實值，Rpre,d-Ract,d表示預測誤差，通常MAPE越低代表模型的預測準確度越高。

采用文獻[15]中定義的準確預測數量，對于多個樣本的預測結果，預測值在真實值周圍±α區域[Ract-αLTTF,Ract+αLTTF]內時可視為準確預測，其中Ract為發動機真實剩余壽命值，LTTF為發動機整個壽命長度，α可根據情況設定，取α=0.05。

3 方法驗證

采用NASA Ames研究中心提供的發動機退化仿真數據集來驗證本文方法，該數據集是基于CMAPSS仿真軟件通過一系列航空渦扇發動機從初始狀態運行至失效的仿真實驗獲得。該數據集分為4組，每組數據由獨立的仿真實驗生成，包括訓練數據集、測試數據集和RUL真實值三個部分，同組數據集中的發動機可視為完全相似的發動機，其中訓練數據和測試數據的格式為26維時間序列，分別記錄了發動機編號、工作循環數、3個發動機工作環境參數(飛行高度、馬赫數和油門桿角度)和21個傳感器測量值，如表1所示。

表1 傳感器數據集描述Table 1 Description of sensor data sets

每臺發動機記錄數據結束時的Cycle值代表該發動機的壽命。采用數據集1來開展實例研究和方法驗證，以其中的#7號發動機為例，訓練數據和測試數據分別如表2、表3所示。

表2 #7號發動機訓練數據Table 2 Training data from seventh engine

觀察21個傳感器測量值序列，發現某些參數為常值，如參數1、5、6、10、16、18、19，在數據集中將這幾個序列刪除，把剩下的序列作為初始數據進行核主元分析，根據交叉驗證的方法確定最佳核函數寬度為3.0，設定最大迭代次數為 1 000，當累計貢獻率達到86.7%時得到前3個主元，其時間序列如圖2 所示。

圖2 #7號發動機的核主元特征Fig.2 KPCA trajectories of seventh engine

選取最初5個飛行循環(令H=1)，以及最后5個飛行循環(令H=0)，將構造的矩陣作為非線性融合模型的訓練樣本，經計算后可得到模型的系數為[a,b1,b2,b3]=[-196.9,0.023 65,-0.011 02,-0.030 3]，常數項c為197.7。

使用上述的非線性模型對訓練數據的3個主元序列進行變換，得到#7號發動機完整的健康指數序列，如圖3所示。

圖3 #7號發動機健康指數序列Fig.3 H series of seventh engine

將由訓練數據集得到的發動機完整H序列作為RVM模型的訓練樣本，經訓練后可得到模型參數。再將由測試數據集得到的H序列作為測試樣本，使用訓練后得到的RVM模型對測試樣本進行外推預測，設定失效閾值為H=0，超過失效閾值時預測結束，截取預測過程的飛行循環數作為RUL預測值，得到的預測結果和RUL的概率分布分別如圖4、圖5所示。

圖5 #7號發動機RUL概率分布Fig.5 Probability distribution of seventh engine’s RUL

圖4 #7號發動機RUL預測結果Fig.4 RUL prediction result of seventh engine

#7號發動機的剩余壽命預測期望ERUL=92，接近其真實壽命(91)。使用本文方法對數據集1中的100個發動機分別進行剩余壽命預測，將預測結果與真實值進行對比，如圖6所示。

圖6 數據集1所有發動機RUL預測結果Fig.6 RUL prediction results of all engines from flight data1

經過對比發現，預測起始時如果發動機H較低，得到的剩余壽命預測值更接近于真實值，預測準確度更高，分析其原因，可能是由于發動機H較高時，現有的H序列反映的發動機退化趨勢不明顯，導致預測的失效時間存在一定程度的滯后，這與發動機投入使用早期難以預測壽命的實際情況符合。

根據預測結果，按照相關指標與SVM模型和過程神經網絡模型的預測結果進行對比。在SVM模型中，選用高斯核函數，通過網格搜索方法優化模型參數，懲罰因子C、核函數參數γ、不敏感損失函數參數ε分別設置為12、0.002、1； 在過程神經網絡模型中，神經網絡層數設置為3，神經元激勵函數選擇Sigmoid函數，各層神經元個數為{10、14、1}。通過相關性分析，選擇2、3、4、7、11、12、14號傳感器數據作為特征參數進行剩余壽命預測。三種方法的預測結果對比如表4所示。

表4 三種預測方法結果對比Table 4 Comparison of the results of three prediction methods

由表4可知，本文方法在各項指標上都明顯優于SVM模型和過程神經網絡模型。傳統的機器學習方法直接使用原始數據來預測發動機剩余壽命，實現多維輸入到一維輸出的直接映射，由于發動機監測數據維數較多，并且各性能參數之間存在一定的相關性，因此預測誤差相對較大。而本文方法利用KPCA對原始數據進行特征提取，有效地消除了冗余信息和參數相關性的干擾，并且特征信息的融合使得RVM模型對發動機剩余壽命的預測更易于實現、預測精度更高。

4 結論

提出了一種基于信息融合與RVM的民航發動機剩余壽命預測方法。該方法是一種間接的RUL預測方法，以發動機失效數據為支撐，能夠融合多維退化特征屬性為一維健康指數，根據健康指數的發展來預測剩余壽命，預測的準確度較高，符合實際應用要求。相比于通過非線性映射能力比較強的機器學習算法(如SVM和神經網絡等)建立監測數據到剩余壽命的關系來實現直接預測的方法，本文方法的預測性能更好，說明其更適用于性能退化情況較為復雜的民航發動機。不過當發動機退化軌跡處于較早時期時，該方法的預測精度會有一定程度的降低，未來需要進一步研究提高預測精度的方法。