基于相關(guān)向量機的發(fā)動機剩余壽命預測

彭鴻博, 蔣雄偉

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

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

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

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

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

1 退化特征信息融合

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

1.1 基于KPCA的退化特征提取

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

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

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

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

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

(1)

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

K′=K-INK-KIN+INKIN

(2)

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

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

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

(3)

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

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

2,…,p

(4)

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

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

1.2 健康指數(shù)的計算方法

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

y=aeX′BT+c

(5)

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

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

(6)

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

2 基于RVM的發(fā)動機RUL預測

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

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

2.1 RVM模型

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

(7)

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

RVM模型的預測過程如下。

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

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

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

(8)

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

(9)

(10)

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

μtest=μTφ(xtest)

(11)

(12)

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

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

2.2 模型評價指標

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

(13)

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

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

3 方法驗證

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

表1 傳感器數(shù)據(jù)集描述Table 1 Description of sensor data sets

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

表2 #7號發(fā)動機訓練數(shù)據(jù)Table 2 Training data from seventh engine

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

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

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

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

圖3 #7號發(fā)動機健康指數(shù)序列Fig.3 H series of seventh engine

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

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

圖4 #7號發(fā)動機RUL預測結(jié)果Fig.4 RUL prediction result of seventh engine

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

圖6 數(shù)據(jù)集1所有發(fā)動機RUL預測結(jié)果Fig.6 RUL prediction results of all engines from flight data1

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

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

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

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

4 結(jié)論

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