基于指數(shù)平滑和XGBoost的航空發(fā)動(dòng)機(jī)剩余壽命預(yù)測(cè)

賴儒杰 范啟富

（上海交通大學(xué)電子信息與電氣工程學(xué)院 系統(tǒng)控制與信息處理教育部重點(diǎn)實(shí)驗(yàn)室）

航空發(fā)動(dòng)機(jī)是飛行器的核心部件， 其穩(wěn)定、可靠地運(yùn)行是飛行安全的基礎(chǔ)。 傳統(tǒng)的發(fā)動(dòng)機(jī)維護(hù)采取簡(jiǎn)單的被動(dòng)和主動(dòng)策略，即當(dāng)發(fā)動(dòng)機(jī)出現(xiàn)故障時(shí)進(jìn)行零件維修或更換，或者發(fā)動(dòng)機(jī)運(yùn)行到默認(rèn)時(shí)限則對(duì)零件統(tǒng)一更新［1］。 由于航空發(fā)動(dòng)機(jī)結(jié)構(gòu)復(fù)雜， 運(yùn)行狀態(tài)受眾多環(huán)境因素的影響，其使用壽命范圍較大，固定的默認(rèn)時(shí)限設(shè)置無法權(quán)衡維護(hù)策略的安全性和經(jīng)濟(jì)性。 近年來針對(duì)航空發(fā)動(dòng)機(jī)的故障預(yù)測(cè)與健康管理 （Prognostics and Health Management，PHM）的研究廣泛開展，通過準(zhǔn)確預(yù)測(cè)剩余壽命 （Remaining Useful Life，RUL），可以對(duì)發(fā)動(dòng)機(jī)的健康狀態(tài)即時(shí)監(jiān)測(cè)，采取合理措施進(jìn)行維護(hù)，從而在保證飛行器可靠飛行的同時(shí)，降低發(fā)動(dòng)機(jī)維護(hù)成本。

PHM領(lǐng)域的RUL預(yù)測(cè)方法主要分為3類：基于物理模型、 基于數(shù)據(jù)驅(qū)動(dòng)和二者結(jié)合的方法［2，3］?；谖锢砟Ｐ偷姆椒ㄐ枰獙?duì)部件進(jìn)行退化過程建模，高度依賴部件原理的相關(guān)知識(shí)；基于數(shù)據(jù)驅(qū)動(dòng)的方法是利用統(tǒng)計(jì)學(xué)、機(jī)器學(xué)習(xí)及深度學(xué)習(xí)等方法對(duì)多源傳感器的數(shù)據(jù)進(jìn)行分析從而得出預(yù)測(cè)模型；二者結(jié)合的方法則在提出物理模型的基礎(chǔ)上，用數(shù)據(jù)驅(qū)動(dòng)方法對(duì)模型參數(shù)進(jìn)行學(xué)習(xí)和調(diào)優(yōu)。

由于航空發(fā)動(dòng)機(jī)復(fù)雜的結(jié)構(gòu)和工作原理使得準(zhǔn)確建模極為困難，對(duì)RUL預(yù)測(cè)采用上述第1類和第3類方法的文獻(xiàn)較少， 基于數(shù)據(jù)驅(qū)動(dòng)的方法則被廣泛研究。 傳統(tǒng)機(jī)器學(xué)習(xí)方面：文獻(xiàn)［3］使用滑動(dòng)時(shí)間窗處理原始信號(hào)的時(shí)間序列，對(duì)時(shí)間窗內(nèi)信號(hào)提取平均值和趨勢(shì)系數(shù)作為特征，最后利用支持向量回歸（SVR）進(jìn)行擬合。 該方法的特征處理簡(jiǎn)單有效，模型參數(shù)少，易于調(diào)整，但在官方評(píng)價(jià)指標(biāo)上預(yù)測(cè)效果欠佳。 深度學(xué)習(xí)方面：文獻(xiàn)［4，5］分別使用卷積神經(jīng)網(wǎng)絡(luò)（CNN）和長(zhǎng)短期記憶網(wǎng)絡(luò)（LSTM）搭建深度模型進(jìn)行特征提取，后接全連接網(wǎng)絡(luò)進(jìn)行回歸；文獻(xiàn)［6］提出一種基于雙邊長(zhǎng)短期記憶網(wǎng)絡(luò)（BiLSTM）的自編碼器（AE）結(jié)構(gòu)計(jì)算健康指數(shù)（HI），再對(duì)HI進(jìn)行線性回歸；文獻(xiàn)［7］則將神經(jīng)網(wǎng)絡(luò)與集成學(xué)習(xí)結(jié)合，使用差分進(jìn)化算法訓(xùn)練出一種以深度信念網(wǎng)絡(luò)（DBF）為基學(xué)習(xí)器的集成模型。 深度學(xué)習(xí)具有強(qiáng)大的非線性擬合能力，但由于其結(jié)構(gòu)復(fù)雜，通常需要大量數(shù)據(jù)才能較好擬合。 文獻(xiàn)［6，7］雖然從模型結(jié)構(gòu)和訓(xùn)練方式上進(jìn)行改進(jìn)， 但受制于數(shù)據(jù)樣本較小，預(yù)測(cè)效果提升有限。

由于航空發(fā)動(dòng)機(jī)退化過程中傳感器變化呈現(xiàn)明顯的趨勢(shì)性，合理的時(shí)域方法可以有效提取特征。 另一方面，一些新興的機(jī)器學(xué)習(xí)方法在小樣本數(shù)據(jù)集上也表現(xiàn)出強(qiáng)大的回歸能力，如Chen T Q和Guestrin C于2016年提出XGBoost，該算法的卓越性能在眾多數(shù)據(jù)競(jìng)賽中得到廣泛認(rèn)可［8］。 因此，筆者將傳統(tǒng)信號(hào)處理方法與XGBoost結(jié)合：首先根據(jù)傳感器信號(hào)的變化趨勢(shì)進(jìn)行特征篩選，使用指數(shù)平滑對(duì)傳感器數(shù)據(jù)降噪處理，利用時(shí)間窗增強(qiáng)特征，最后使用XGBoost進(jìn)行回歸。 筆者選擇在NASA C-MAPSS航空發(fā)動(dòng)機(jī)仿真數(shù)據(jù)集上進(jìn)行試驗(yàn)。

1 指數(shù)平滑

指數(shù)平滑（Exponential Smoothing，ES）常采用指數(shù)加權(quán)移動(dòng)平均（Exponential Weighted Moving Average，EWMA）實(shí)現(xiàn)。 EWMA是一種監(jiān)測(cè)過程平均值的統(tǒng)計(jì)方法，其一般形式為［9］：

其中，xt為t時(shí)刻真實(shí)值，yt為t時(shí)刻觀測(cè)值，衰減因子α∈（0，1）（α為定值）。

筆者采用Python數(shù)據(jù)分析庫(kù)Pandas中提出的一種修正權(quán)重的觀測(cè)值計(jì)算方式，即各項(xiàng)系數(shù)為指數(shù)函數(shù)的一般加權(quán)平均：

式（2）適用于有限長(zhǎng)序列的觀測(cè)值計(jì)算，且滿足t→∞時(shí)與式（1）等價(jià)。

對(duì)于任一給定α，在計(jì)算yt時(shí)，由于距離當(dāng)前時(shí)刻較遠(yuǎn)的xt-i權(quán)重較小， 對(duì)計(jì)算的影響可忽略，因此可理解為近端一定范圍s內(nèi)的xt-i起主要作用，即{xt-i|0≤i≤s-1}。 視s為平滑系數(shù)，α與s關(guān)系如下：

2 XGBoost

XGBoost 是廣義梯度提升決策樹（Gradient Boosting Decision Tree，GBDT） 的一種高效實(shí)現(xiàn)。GBDT以迭代的方式構(gòu)建多棵CART樹，對(duì)每一輪的預(yù)測(cè)殘差逐步進(jìn)行擬合。 XGBoost在GBDT的基礎(chǔ)上對(duì)目標(biāo)函數(shù)引入L2正則項(xiàng)和二階導(dǎo)數(shù)，并從系統(tǒng)設(shè)計(jì)的角度對(duì)運(yùn)算性能進(jìn)行優(yōu)化［8］。

考慮數(shù)據(jù)集D={（xi，yi）}有n個(gè)樣本，每個(gè)樣本m維特征（|D|=n，xi∈Rm，yi∈R）。 訓(xùn)練得到K棵樹的集成模型，其預(yù)測(cè)結(jié)果表示為：

其中，Γ={f（x）=ωq（x）}（q：Rm→T，ω∈RT）為K棵樹的集合，q為樹的結(jié)構(gòu)，T為樹的葉節(jié)點(diǎn)個(gè)數(shù)，ω為葉節(jié)點(diǎn)權(quán)重。

目標(biāo)函數(shù)定義為：

其中，l為損失函數(shù)，Ω為復(fù)雜度函數(shù)， ‖ω‖2為葉節(jié)點(diǎn)權(quán)重的L2范數(shù)，懲罰系數(shù)γ和λ用于調(diào)整模型復(fù)雜度，控制過擬合。

損失函數(shù)可自定義， 但需滿足二階可導(dǎo)，對(duì)于回歸問題通常采用平方誤差：

訓(xùn)練階段，每一輪增加一棵樹對(duì)上一輪的預(yù)測(cè)殘差進(jìn)行擬合，第t輪的目標(biāo)函數(shù)如下：

對(duì)損失函數(shù)進(jìn)行二階泰勒展開，并去除常數(shù)項(xiàng)，目標(biāo)函數(shù)改寫為：

每一輪迭代采用貪心算法構(gòu)建CART樹，詳細(xì)過程見文獻(xiàn)［8］。

3 數(shù)據(jù)預(yù)處理和特征工程

試驗(yàn)數(shù)據(jù)來源于NASA C-MAPSS航空發(fā)動(dòng)機(jī)仿真數(shù)據(jù)集中FD001子集。監(jiān)測(cè)數(shù)據(jù)包括3個(gè)操作變量和21個(gè)傳感器數(shù)據(jù)。 訓(xùn)練集包含100臺(tái)發(fā)動(dòng)機(jī)從啟動(dòng)到失效各個(gè)運(yùn)行周期的監(jiān)測(cè)數(shù)據(jù)，最長(zhǎng)周期數(shù)362cycle，最短周期數(shù)128cycle，共20 631個(gè)樣本； 測(cè)試集包含100臺(tái)發(fā)動(dòng)機(jī)從啟動(dòng)到某一運(yùn)行周期（中斷時(shí)刻）的監(jiān)測(cè)數(shù)據(jù)及對(duì)應(yīng)RUL，最長(zhǎng)周期數(shù)303cycle，最短周期數(shù)31cycle，共13 096個(gè)樣本。

訓(xùn)練集的RUL標(biāo)簽設(shè)置通常有兩種方式：線性退化和分段線性退化。 考慮運(yùn)行初期發(fā)動(dòng)機(jī)性能較穩(wěn)定，退化不明顯，筆者采用分段模型，初始RUL設(shè)置為127cycle［2］，訓(xùn)練集的RUL標(biāo)簽設(shè)置如圖1所示。

傳感器數(shù)據(jù)變化趨勢(shì)見表1。 筆者忽略數(shù)值恒定的傳感器，選取剩余14個(gè)傳感器的信號(hào)作為輸入特征［7］。 由于原始信號(hào)噪聲較大，采用式（2）進(jìn)行指數(shù)平滑處理，使預(yù)測(cè)算法更好地捕捉信號(hào)的趨勢(shì)變化。 圖2展示訓(xùn)練集1號(hào)發(fā)動(dòng)機(jī)3號(hào)傳感器信號(hào)的平滑效果， 可以看出平滑系數(shù)s越大，信號(hào)波動(dòng)越小。

圖1 訓(xùn)練集RUL標(biāo)簽設(shè)置

表1 傳感器數(shù)據(jù)變化趨勢(shì)

圖2 不同平滑系數(shù)時(shí)的平滑效果

針對(duì)多變量時(shí)間序列問題，僅使用單個(gè)樣本點(diǎn)的信息會(huì)使得數(shù)據(jù)利用不充分，因此使用定長(zhǎng)時(shí)間窗增加各傳感器連續(xù)時(shí)間步歷史數(shù)據(jù)［7］，達(dá)到特征增強(qiáng)的目的。 同時(shí)由于運(yùn)行時(shí)長(zhǎng)也是影響RUL的重要因素之一， 運(yùn)行時(shí)長(zhǎng)也被作為一維特征。 設(shè)時(shí)間窗寬度為Ntw（含當(dāng)前時(shí)刻），選取傳感器14個(gè)， 則特征工程處理后每個(gè)樣本點(diǎn)包含14×Ntw+1維特征。 考慮測(cè)試集最短周期數(shù)為31cycle，本試驗(yàn)中Ntw取值為30。

由于各傳感器真實(shí)數(shù)值差異較大，訓(xùn)練前需對(duì)各特征進(jìn)行歸一化處理， 筆者采用Min-Max歸一化，將特征值映射到［0，1］范圍：

4 試驗(yàn)結(jié)果分析

該數(shù)據(jù)集通常采用Score、RMSE和Accuracy指標(biāo)對(duì)預(yù)測(cè)結(jié)果進(jìn)行評(píng)價(jià)［3］（其中Score為該數(shù)據(jù)集的官方指標(biāo)）：

其 中，di為 預(yù) 測(cè) 誤 差，di=RULest，i-RULactual，i，RULest，i和RULactual，i分別為第i個(gè)樣本的RUL預(yù)測(cè)值和 真 實(shí) 值，Cor（di）表 示 當(dāng) 預(yù) 測(cè) 誤 差di屬 于［-13，10］區(qū)間時(shí)認(rèn)為該樣本估計(jì)準(zhǔn)確。

繪制的指標(biāo)評(píng)價(jià)結(jié)果如圖3所示。 Score為非對(duì)稱指標(biāo)，對(duì)過預(yù)測(cè)懲罰較大，對(duì)欠預(yù)測(cè)懲罰較小；RMSE為對(duì)稱指標(biāo)；Accuracy用于統(tǒng)計(jì)在合理誤差范圍內(nèi)的預(yù)測(cè)結(jié)果比例。

圖3 評(píng)價(jià)指標(biāo)結(jié)果

試驗(yàn)使用Python 語言下scikit-learn 接口的XGBoost回歸模型，在經(jīng)指數(shù)平滑和特征工程處理后的數(shù)據(jù)集上進(jìn)行訓(xùn)練和測(cè)試， 訓(xùn)練中損失函數(shù)使用XGBoost 算法默認(rèn)的平方誤差 （式（6））。對(duì)于不同平滑系數(shù)處理后的訓(xùn)練集，分別采用5折交叉驗(yàn)證進(jìn)行超參數(shù)選擇，得到各條件下的模型最優(yōu)參數(shù)， 然后使用訓(xùn)練集所有數(shù)據(jù)對(duì)最優(yōu)參數(shù)下的模型進(jìn)行訓(xùn)練， 最后對(duì)測(cè)試集中斷時(shí)刻的RUL進(jìn)行預(yù)測(cè)， 并結(jié)合其真實(shí)RUL進(jìn)行效果評(píng)估，試驗(yàn)結(jié)果見表2（表中粗體字為該列最佳值）。

表2 不同平滑系數(shù)取值下XGBoost預(yù)測(cè)效果對(duì)比

由表2可知， 隨著平滑系數(shù)的增大平滑效果增 強(qiáng)，Score 整 體 呈 減 小 趨 勢(shì)，RMSE 減 小，Accuracy提高，即算法性能在3個(gè)指標(biāo)衡量下均得以提升，由此推測(cè)XGBoost算法對(duì)噪聲較敏感，使用平滑后的數(shù)據(jù)作為輸入，可使算法更好地捕捉趨勢(shì)特征，取得更好的預(yù)測(cè)效果。 平滑系數(shù)為20時(shí)，XGBoost的參數(shù)列表和預(yù)測(cè)結(jié)果見表3、圖4。

表3 XGBoost參數(shù)列表

圖4 XGBoost預(yù)測(cè)結(jié)果

進(jìn)一步將筆者所提方法與其他主流方法的預(yù)測(cè)效果進(jìn)行比較分析，結(jié)果見表4（表中粗體字為該列最佳值）。 筆者所提方法在RMSE指標(biāo)上取得最好結(jié)果， 考慮XGBoost采用平方誤差作為損失函數(shù)進(jìn)行訓(xùn)練， 可認(rèn)為充分發(fā)揮了XGBoost的擬合能力。 對(duì)比BiLSTM-ED方法，BiLSTM-ED的Score值更小， 而RMSE略高， 且Accuracy明顯較低， 推測(cè)該方法傾向于欠預(yù)測(cè)RUL。 對(duì)比SVR方法，SVR的Accuracy略高于筆者所提方法，而Score則明顯較高，推測(cè)該方法在誤差邊界外的預(yù)測(cè)值偏離真實(shí)值較大。

表4 筆者方法與其他主流方法預(yù)測(cè)效果對(duì)比

綜上所述，筆者所提方法預(yù)測(cè)值相對(duì)于真實(shí)值總體偏差最小， 且預(yù)測(cè)沒有明顯偏向性；在RMSE指標(biāo)上取得最優(yōu)結(jié)果，Score和Accuracy指標(biāo)略差于最好方法，但顯著優(yōu)于其他方法，因此可以認(rèn)為筆者所提方法綜合性能最佳。

5 結(jié)束語

針對(duì)航空發(fā)動(dòng)機(jī)退化過程多源傳感器測(cè)量信號(hào)，使用修正權(quán)值的指數(shù)平滑方法對(duì)有限長(zhǎng)時(shí)間序列數(shù)據(jù)進(jìn)行降噪處理，并添加時(shí)間窗內(nèi)歷史數(shù)據(jù)進(jìn)行特征增強(qiáng)， 使XGBoost的預(yù)測(cè)性能得到大幅提升。 通過與其他前沿算法進(jìn)行對(duì)比，可驗(yàn)證筆者所提方法在對(duì)稱指標(biāo)上預(yù)測(cè)效果最佳，在非對(duì)稱指標(biāo)上僅次于最優(yōu)方法，因此可體現(xiàn)出方法優(yōu)異的綜合性能。 下一步工作可研究其他信號(hào)處理算法進(jìn)行數(shù)據(jù)平滑（如卡爾曼濾波）以及使用非對(duì)稱損失函數(shù)訓(xùn)練XGBoost對(duì)其預(yù)測(cè)性能的影響。