劉小峰,史長振,黃洪升,柏 林

(重慶大學機械與運載工程學院,重慶 400044)

1 引言

航空發(fā)動機的剩余壽命預(yù)測(remaining useful life,RUL)對于保障發(fā)動機運行安全可靠,降低其維修成本具有重要的意義[1].鑒于航空發(fā)動機性能衰退物理模型的復(fù)雜性,利用運行狀態(tài)監(jiān)測數(shù)據(jù)構(gòu)建健康指數(shù)(health index,HI)并采用智能網(wǎng)外推HI至失效閾值(failure threshold,FT),逐漸成為其RUL預(yù)測主流研究方法.

為了從航空發(fā)動機監(jiān)測數(shù)據(jù)中構(gòu)建能夠表征其性能退化狀態(tài)的HI,許多國內(nèi)外研究學者通過特征的優(yōu)化選擇以單值原始物理監(jiān)測量(如EGTM,DEGT等)作為發(fā)動機的HI[2-3],或者通過狀態(tài)監(jiān)測量的降維約簡、加權(quán)融合、距離度量等方法得到發(fā)動機的綜合HI[4-5].這些方法亦或?qū)ΡO(jiān)測參量進行人工主觀篩選,亦或在HI構(gòu)建中未加入與RUL相關(guān)的約束條件,無法保證HI與RUL的相關(guān)性.隨著機器學習技術(shù)的發(fā)展,基于智能學習網(wǎng)絡(luò)的HI構(gòu)建方法及其RUL預(yù)測模型得到了長足的發(fā)展,如Babu等[6]建立的基于深度卷積神經(jīng)網(wǎng)絡(luò),Nieto等[7]提出的粒子群優(yōu)化SVM,Martha A等[8]采用的貝葉斯層次模型,Yang等[9]開發(fā)的量子行為算法優(yōu)化的極限學習機,彭開香等[10]提出深度置信網(wǎng)絡(luò)與隱馬爾科夫融合模型,Yuan等[11]采用的長短期記憶神經(jīng)網(wǎng)絡(luò)算法等.要指出的是,實際的航空發(fā)動機運行任務(wù)剖面多雜,工況環(huán)境多變,在測試發(fā)動機與訓(xùn)練發(fā)動機退化過程差異性較大的情況下,預(yù)測得到的RUL誤差往往較大.

目前基于HI外推的RUL預(yù)測研究中采用的失效閾值,大多采用的是人工經(jīng)驗法[12-13]、失效閾值分布統(tǒng)計方法[13-14].前者根據(jù)主觀經(jīng)驗設(shè)置,缺乏科學依據(jù),后者需要以大量全壽命周期的衰退數(shù)據(jù)樣本為支撐.Zhang等[15]通過相似性壽命預(yù)測方法對渦扇發(fā)動機仿真數(shù)據(jù)進行研究,預(yù)測結(jié)果遠優(yōu)于神經(jīng)網(wǎng)絡(luò)方法.曹惠玲等[16]通過相關(guān)分析提取與發(fā)動機性能衰退密切相關(guān)參數(shù),提出了一種基于時間序列相似性匹配的航空發(fā)動機剩余壽命預(yù)測方法.張妍等[17]在特征相似性匹配基礎(chǔ)上對航空發(fā)動機進行RUL進行了加權(quán)融合估計.上述這些方法大多采用特征序列間的歐式距離對發(fā)動機退化過程的進行相似估計,尚未涉及失效閾值的設(shè)置.另者,這種以序列距離匹配為基礎(chǔ)的相似度量方法只考慮類序列的全局相似性,容易受到噪聲干擾的影響,且沒有考慮每個退化序列的個體差異性.

針對以上存在的問題,本論文結(jié)合Box-Cox變換、消噪自編碼器(denoising autoencoder,DAE)及多重核典型相關(guān)(multiple kernel canonical correlation,MKCC)技術(shù),強化了HI指數(shù)與飛行循環(huán)間的相關(guān)性;在充分考慮了發(fā)動機性能退化的個體差異性與相似性的情況下,引入了全階時間冪灰色預(yù)測模型(gray forecasting model with full order time power terms,FOTPGM(1,1))及KL散度加權(quán)方法進行FT的設(shè)置.

2 航空發(fā)動機強化健康指數(shù)構(gòu)建

2.1 監(jiān)測物理量的Box-Cox變換

航空發(fā)動機監(jiān)測的物理參數(shù)眾多,如排氣溫度、燃料流量、旁路比率等,這些參數(shù)受環(huán)境工況、負載等影響隨機性較大,規(guī)律性不強,對航空發(fā)動機性能退化的表征能力不強.為了增強監(jiān)測參量與發(fā)動機RUL間的相關(guān)性,采用Box-Cox變換[18]強化傳感器監(jiān)測參數(shù)與發(fā)動機RUL間的相關(guān)性,使得變換后的參量與RUL具有較好相關(guān)性,同時增強參量之間的獨立性.設(shè)航空發(fā)動機監(jiān)測的物理參量集X包含N個樣本,即X={x1,x2,···,xn,···,xN},每個樣本中共有M個監(jiān)測參量,即xn={xn1,xn2,···,xnm,···,xnM},xnm表示第n個飛行循環(huán)下的第m個監(jiān)測參量幅值.X對應(yīng)的飛行循環(huán)量為L={l1,l2,···,ln,···,lN},采用Box-Cox變換參數(shù)λ={λ1,λ2,···,λm,···,λM}得到變換后監(jiān)測參量X(λ)=,其中:

式中a為使監(jiān)測參量X整體平移后大于0的常數(shù).構(gòu)造如式(2)中似然函數(shù):

2.2 性能退化特征提取

航空發(fā)動機物理監(jiān)測參量受環(huán)境工況、負載等影響存在較強的隨機性,本文采用DAE弱化噪聲及環(huán)境工況等不確定因素的影響的同時從航空發(fā)動機多維監(jiān)測參量中提取深層次的退化狀態(tài)信息.將經(jīng)過Box-Cox變換后的航空發(fā)動機監(jiān)測物理參量X(λ),隨機加入噪聲后記為,然后對進行編碼,獲得隱含層表示為[19]

式 中:D(·)為解碼器的非線性轉(zhuǎn)換,g(·)為和之間的映射函數(shù);θD為解碼器的參數(shù).DAE網(wǎng)絡(luò)通過最小化與X(λ)之間的重構(gòu)誤差實現(xiàn)網(wǎng)絡(luò)的訓(xùn)練學習,得到的編碼向量h={h(1),h(2),···,h(p)},p為編碼向量的維度,即為輸入樣本對應(yīng)的發(fā)動機退化特征表達.

2.3 基于多重核典型相關(guān)的HI構(gòu)建

由于DAE提取退化特征時,并沒有加入與發(fā)動機RUL相關(guān)的約束條件,無法保證提取的退化特征與RUL的相關(guān)性.因此,本文采用了多重核典型相關(guān)算法[20]對提取的退化特征進行強化處理,通過核函數(shù)學習一組相關(guān)投影方向,使不同退化特征投影到核空間后與發(fā)動機飛行循環(huán)數(shù)擁有最大的相關(guān)性,進而采用核空間的融合特征值作為發(fā)動機的強化健康指數(shù).將H={h1,h2,···,hn,···,hN}∈Rp×N記為DAE提取的退化特征,對應(yīng)的飛行循環(huán)數(shù)L={l1,l2,···,ln,···,lN}∈R1×N,將兩者組合成(p+1)×N維的樣本集合Y={Y(1),Y(2),···,Y(p+1)},每個子集的維度都為N.采用MKCC將p+1個子集映射到核空間得到:[Φ1,Φ2,···,Φp+1].多重核典型相關(guān)的目標函數(shù)為

3 自適應(yīng)失效閾值設(shè)置

3.1 參考發(fā)動機FT預(yù)測

航空發(fā)動機的性能退化是由多因素引起的,每個因素對發(fā)動機RUL的影響是未知不確定的,其性能退化問題本身就存在于一個灰色系統(tǒng)中,本文采用的全階時間冪灰色預(yù)測模型FOTP-GM(1,1)可以自適應(yīng)地根據(jù)HI序列動態(tài)變化改變模型結(jié)構(gòu)和參數(shù),提高FT的擬合精度與發(fā)動機的RUL預(yù)測精度.

在FOTP-GM(1,1)模型[21]中,x(0)(k)滿足

其中:a稱為發(fā)展系數(shù),bi(i=1,2,···,h)稱為灰色作用量,h稱為時間冪項bikh?i的階數(shù),h最大為4.x(0)(k)為k時刻的原始信號序列值;z(1)(k)為原始信號累加序列的緊鄰均值生成序列值.采用最小二乘法,可估計出參數(shù)序列=(a,b1,b2,···,bh),根據(jù)式(7)即可求得FTOP-GM(1,1)模型的時間響應(yīng)函數(shù)為

1.1.3 試驗動物和細胞。SD大鼠，第四軍醫(yī)大學試驗動物中心；黑色素瘤B16細胞株，中國科學院上海細胞生物研究所。

其中c為常數(shù).為原始信號的累加序列.則x(0)(k)的FTOP-GM擬合值表示為

由于FOTP-GM(1,1)模型具有無差擬合近似非齊次指數(shù)序列的性能,采用FOTP-GM(1,1)模型對訓(xùn)練發(fā)動機的HI曲線進行擬合,能夠有效避免HI曲線在最終壽命時間點上的端點畸變問題,從而降低FT設(shè)置偏差.全壽命參考發(fā)動機l的HI序列Hl(t),t=t1,t2,···,tl,在最終壽命時間點t=tl上的FTOP-GM(1,1)擬合值,即為該參考發(fā)動機的FT,即

3.2 基于KL距離HI相似匹配的FT加權(quán)設(shè)置

在參考發(fā)動機HI的基礎(chǔ)上,如何確定測試發(fā)動機的FT是影響其RUL預(yù)測準確性的關(guān)鍵.由于各個發(fā)動機退化軌跡的差異性與相似性,本文采用基于KL距離壽命調(diào)節(jié)函數(shù)來綜合確定測試發(fā)動機的失效閾值.計算待測發(fā)動機HI與各個全壽命參考發(fā)動機HI之間的KL距離并以此作為發(fā)動機退化軌跡相似度量指標,選取5個最小KL距離構(gòu)造閾值調(diào)節(jié)函數(shù),結(jié)合對應(yīng)的5個全壽命發(fā)動機的FT,對測試發(fā)動機FT進行加權(quán)設(shè)置,即

其 中:l=1,2,···,5,KL(p(HItest)‖p(HIl))表示測試發(fā)動機的HI序列分布與訓(xùn)練發(fā)動機HI序列分布的KL距離,p(·)表示HI概率分布,ωl代表第l個發(fā)動機的失效閾值的加權(quán)系數(shù).

3.3 基于HI外推的RUL預(yù)測

FOTP-GM(1,1)模型具有自適應(yīng)的結(jié)構(gòu)特征,其結(jié)構(gòu)參數(shù)會根據(jù)發(fā)動機性能衰退演化趨勢動態(tài)地自適應(yīng)調(diào)節(jié),不僅可以精確模擬齊次指數(shù)HI序列,而且可以無差模擬近似非齊次指數(shù)HI 序列,因此具有更高的預(yù)測精度,更強的泛化能力,更優(yōu)的自適應(yīng)性能.設(shè)測試發(fā)動機的歷史HI曲線為HItest(l),l=l1,l2,···,lc,根據(jù)式(11)-(13)對其進行建模得到其對應(yīng)的灰色模型FGMtest,則其RUL估計值為

4 基于強化HI閾值設(shè)置的發(fā)動機RUL預(yù)測

本文提出了基于強化HI相似度量的失效閾值設(shè)定方法,如圖1所示.該方法主要步驟包括:1)采用Box-Cox算法對發(fā)動機監(jiān)測物理量進行變換,增強其與飛行循環(huán)間的相關(guān)性;2)將Box-Cox變換后的物理參量輸入DAE,提取發(fā)動機性能退化特征量;3)對DAE提取的退化特征量與對應(yīng)飛行循環(huán)數(shù)進行MKCC運算,構(gòu)建發(fā)動機強化HI指數(shù);4)采用FOTP-G(1,1)模型對訓(xùn)練發(fā)動機的全壽命HI進行擬合,得到訓(xùn)練發(fā)動機的失效閾值;5)計算訓(xùn)練發(fā)動機的HI曲線與測試發(fā)動機的HI曲線的KL距離,在此基礎(chǔ)上對訓(xùn)練發(fā)動機的FT進行加權(quán)平均,得到測試發(fā)動機的FT;6)采用FOTPG(1,1)模型對測試發(fā)動機HI進行外推至其設(shè)定的失效閾值,得到其RUL估計值.

圖1 提出方法流程圖Fig.1 Flowchart of proposed method

5 實驗驗證

5.1 數(shù)據(jù)描述

本文選用美國宇航局C-MAPSS航空發(fā)動機數(shù)據(jù)[20]對提出方法進行應(yīng)用驗證.該數(shù)據(jù)集共包括不同工況和失效模式的4組數(shù)據(jù),每組數(shù)據(jù)包含3個工況參數(shù)和21個傳感器監(jiān)測參數(shù).本文選取FD001數(shù)據(jù)集(即第一組)進行實驗驗證.該數(shù)據(jù)集包含100臺全壽命訓(xùn)練數(shù)據(jù)、100臺部分測試數(shù)據(jù)和測試數(shù)據(jù)對應(yīng)的真實RUL.測試數(shù)中每臺發(fā)動機的起始狀態(tài)及觀測序列長度均不同.關(guān)于發(fā)動機具體模塊的信息以及傳感器的說明詳見文獻[22].

5.2 監(jiān)測參量強化處理

FD001數(shù)據(jù)集中,不同的監(jiān)測物理參量有不同的退化趨勢,在21個監(jiān)測物理參量中,有7個監(jiān)測參量在整個退化周期中幅值未發(fā)生變化,為了達到對發(fā)動機性能退化狀態(tài)更準確的描述,剔除這7個監(jiān)測參量,將剩余14個監(jiān)測參量組成特征集.為改善傳感器監(jiān)測參數(shù)與發(fā)動機飛行循環(huán)序列之間的線性相關(guān)性,采用Box-Cox算法實現(xiàn)傳感器監(jiān)測參數(shù)的轉(zhuǎn)換.以#2號發(fā)動機T50、P30、Nc、Ps30、BPR 5個傳感器參數(shù)為例,參數(shù)λ取值范圍為[?5,5],通過最大似然估計法計算所得Box-Cox變換參數(shù)λ最優(yōu)值以及變換前后各傳感器參數(shù)與飛行循環(huán)序列的Pearson系數(shù)對比如下表1所示.從表中可以看出,經(jīng)過Box-Cox變換后,表中所展示的5個傳感器參數(shù)與飛行循環(huán)之間的線性相關(guān)度均得到了一定程度的提升.

表1 Box-Cox變換前后Pearson系數(shù)對比表Table 1 Comparison of Pearson coefficient before and after Box-Cox transformation

5.3 建立健康指標

DAE網(wǎng)絡(luò)參數(shù)設(shè)置為:輸入層節(jié)點個數(shù)為14,第1層隱含層節(jié)點數(shù)為10,第2層隱含層節(jié)點數(shù)為5,學習率為0.001.按照式(3)-(6)對DAE網(wǎng)絡(luò)進行訓(xùn)練,得到維度為5的退化特征向量,再采用MKCC得到的強化HI曲線.圖2給出了#40,#56,#78,#88號訓(xùn)練發(fā)動機的HI曲線,可以看出采用本文方法提取的發(fā)動機HI曲線退化趨勢明顯,波動幅度小,前期HI變化幅度不大,說明發(fā)動機屬于正常狀態(tài),而后期HI下降明顯,說明發(fā)動機進入衰退狀態(tài).

圖2 訓(xùn)練發(fā)動機的強化HI曲線Fig.2 Enhanced health index of training aero-engines

圖3給出了直接采用DAE算法構(gòu)建的HI曲線,相較于圖2,HI曲線的波動性明顯增大,本文構(gòu)建的強化HI更大程度上擯除了噪聲的干擾,能更有效地表征發(fā)動機性能退化過程.

圖3 基于DAE算法的訓(xùn)練發(fā)動機HI圖Fig.3 Health index of training aero-engines based on DAE algorithm

5.4 失效閾值確定

以FD001中的#43,#72,#87,#99號測試發(fā)動機為分析對象,對100個全壽命參考發(fā)動機的HI指數(shù)曲線進行KL距離匹配,得到5個KL距離最小的參考發(fā)動機如表3所示.對這5個全壽命發(fā)動機的HI進行FOTPGM(1,1)模型擬合后,得到最后失效點對應(yīng)的FT值,繼而按照式(11)-(12)計算得到#43,#72,#87,#99號測試發(fā)動機的失效閾值,結(jié)果如表2所示.

表2 測試發(fā)動機失效閾值Table 2 FT of test aero-engines

5.5 RUL預(yù)測

在測試發(fā)動機HI曲線與設(shè)定閾值的基礎(chǔ)上,采用FOTP-GM(1,1)模型方法進行測試航空發(fā)動機的剩余壽命預(yù)測.部分測試發(fā)動機剩余壽命預(yù)測結(jié)果如圖4所示,具體預(yù)測結(jié)果數(shù)值由下表3給出.圖4(a)與(b)依次為#43,#72號測試發(fā)動機剩余壽命預(yù)測圖,其真實RUL與預(yù)測RUL之間的飛行循環(huán)次數(shù)誤差分別為3和4.由此可見,對于當前飛行循環(huán)次數(shù)較多的測試發(fā)動機由于具有足夠的歷史監(jiān)測數(shù)據(jù),FOTP-GM(1,1)可以較準確地獲得其剩余壽命.圖4(c)與(d)分別是#87,#99號測試發(fā)動機剩余壽命預(yù)測圖,其真實RUL與預(yù)測RUL的誤差為11和9個飛行循環(huán).可見,飛行次數(shù)少的測試發(fā)動機由于歷史數(shù)據(jù)較缺乏,FOTP-GM(1,1)模型獲得的預(yù)測精度盡管有所降低,但仍在可接受范圍之內(nèi).因此FOTP-GM(1,1)同樣適合于短數(shù)據(jù)的長時間預(yù)測.

圖4 測試發(fā)動機RUL預(yù)測圖Fig.4 RUL prediction for test aero-engines

表3 部分測試發(fā)動機RUL預(yù)測結(jié)果Table 3 RUL prediction results of some test aeroengines

5.6 RUL預(yù)測結(jié)果比較分析

圖5(a)給出了100個測試發(fā)動機的RUL預(yù)測結(jié)果與真實RUL的對比分析圖,可見,各個測試發(fā)動機的預(yù)測誤差都各有不同,最小誤差為0,最大誤差為41個飛行循環(huán),這主要和各個測試對象的歷史數(shù)據(jù)長短,波動程度及退化軌跡各異有關(guān).圖5(b)給出了測試RUL誤差(預(yù)測RUL減去真實RUL)的分布情況,76%的測試發(fā)動機RUL預(yù)測誤差集中在?13至10個飛行循環(huán)且提前預(yù)測數(shù)量遠多于滯后預(yù)測,符合工程實際要求,驗證了本文提出方法的準確性.

圖5 100個測試發(fā)動機RUL預(yù)測誤差圖Fig.5 RUL prediction error for 100 test aero-engines

為了進一步驗證本文預(yù)測方法的優(yōu)勢,采用Score指標對預(yù)測結(jié)果進行量化評價,并與卷積神經(jīng)網(wǎng)絡(luò)(CNN)算法[6]和極限學習機模糊聚類(EML-FC)算法[23]預(yù)測RUL所得Score指標進行對比.Score指標[23]計算公式如下:

式中:actRULi,preRULi,Ei和Ai分別是第i個測試發(fā)動機的實際剩余壽命、模型預(yù)測壽命、誤差率和得分;Score代表N個測試發(fā)動機的平均得分,這里N=100.從式(14)-(15)可以看出,Ai與Score越小,說明預(yù)測出的RUL精度越高,方法性能更加穩(wěn)定.經(jīng)計算,本文方法預(yù)測得到的Score值為346,小于EML-FC方法的Score得分1046,CNN方法的得分1287.可見,本文方法對100個發(fā)動機RUL預(yù)測的平均精度大于傳統(tǒng)方法,該方法通過強化HI的失效閾值自適應(yīng)設(shè)置,提高了發(fā)動機RUL預(yù)測精度與可信度.

6 結(jié)論

本文針對航空發(fā)動機監(jiān)測物理參量對發(fā)動機退化狀態(tài)表征能力不強的問題,利用Box-Cox變換、DAE以及MKCC算法,強化了發(fā)動機HI與其飛行循環(huán)數(shù)的相關(guān)性,并采用了基于KL距離的HI相似度量法對其FT進行了自適應(yīng)設(shè)置.試驗結(jié)果表明,本文提出的FT設(shè)置方法,充分考慮了不同發(fā)動機退化模式相似性與差異性,降低了閾值主觀設(shè)置造成的RUL預(yù)測的不確定性,提高了RUL的預(yù)測精度,在航空裝備壽命預(yù)測領(lǐng)域具有較好的應(yīng)用前景.