未知載荷條件下機(jī)械系統(tǒng)剩余壽命預(yù)測方法

許丹，肖小琦，馮至昕

（北京航空航天大學(xué) 可靠性與系統(tǒng)工程學(xué)院，北京 100083）

傳感器技術(shù)的發(fā)展為復(fù)雜的工程系統(tǒng)的狀態(tài)檢測與運行環(huán)境檢測提供了條件。其中，許多退化部件的傳感器信號被用于退化建模以預(yù)測退化部件的剩余壽命。對于設(shè)備退化建模與壽命預(yù)測，現(xiàn)有研究大多集中在單個傳感器信號在單個運行條件下的分析。具體來說，通常先根據(jù)歷史監(jiān)測數(shù)據(jù)建立隨機(jī)退化模型，再根據(jù)現(xiàn)場傳感器實時采集到的監(jiān)測數(shù)據(jù)對模型參數(shù)進(jìn)行更新。周紹華等［1］使用非參數(shù)方法建立了稀疏數(shù)據(jù)單傳感器信號退化模型，其中工作狀態(tài)參數(shù)由極大似然估計（maximum likelihood estimate，MLE）算法獲取。在基于多傳感器信息退化建模方面［2-3］，通常先構(gòu)建健康因子再進(jìn)行退化過程建模。任子強(qiáng)等［4］將多個傳感器數(shù)據(jù)融合成一個復(fù)合健康指標(biāo)，并采用線性維納過程對復(fù)合健康指標(biāo)進(jìn)行退化建模，進(jìn)而得到系統(tǒng)的預(yù)測壽命。近年來，學(xué)者們開始研究基于機(jī)器學(xué)習(xí)的融合方法。馬奇友等［5］利用深度學(xué)習(xí)和長短期記憶組合來構(gòu)造DLSTM網(wǎng)絡(luò)，將多個傳感器信號數(shù)據(jù)進(jìn)行融合處理，得到系統(tǒng)剩余壽命預(yù)測。Chen等［6］提出了一種基于編碼器-解碼器框架的遞歸神經(jīng)網(wǎng)絡(luò)，融合多個傳感器信號數(shù)據(jù)并預(yù)測HI值，通過線性回歸獲得最終的RUL值。但是，以上研究大多是考慮系統(tǒng)在單一工況下工作或者假設(shè)系統(tǒng)未來載荷與歷史載荷服從獨立同分布。目前，基于多傳感器信息變載荷退化建模的研究較少。Javed等［7］提出了一種基于小波極限學(xué)習(xí)機(jī)與減法最大熵模糊聚類的方法對多維變載荷條件下的退化數(shù)據(jù)進(jìn)行降維，并利用降維后的特征來描述退化過程。Yan等［8］以最小化退化模型不確定度、最小化失效閾值不確定度及最大化信息范圍為目標(biāo)進(jìn)行線性優(yōu)化，得到考慮變載荷情況下的傳感器數(shù)據(jù)加權(quán)系數(shù)，構(gòu)建健康因子。而在變載荷退化建模與預(yù)測中，未來載荷的不確定性也是研究的一個難點。Flory等［9］假設(shè)變載荷狀態(tài)為一個馬爾可夫過程，通過蒙特卡羅抽樣與經(jīng)典貝葉斯信息準(zhǔn)則的方法從傳感器信息中估計環(huán)境狀態(tài)參數(shù)，作為變換融合模型的輸入，估計退化過程的首穿時分布。

在許多情況下，系統(tǒng)通常會在多個載荷下工作，且載荷變化的下一個狀態(tài)常常與此前的工作狀態(tài)無關(guān)，即可以假設(shè)為具有馬爾可夫的無后效特性。而隱藏的半馬爾可夫鏈具有用于近似復(fù)雜概率分布的隱馬爾可夫鏈的靈活性和用于表示時間的半馬爾可夫鏈的靈活性，常用于設(shè)備的故障診斷及預(yù)測［10-11］。因此，本文提出一種基于隱半馬爾可夫模型（hidden semi-Markov model，HSMM）的隨機(jī)載荷空間構(gòu)建和預(yù)測的方法。首先，采用HSMM對隨機(jī)載荷空間進(jìn)行建模。然后，結(jié)合多傳感器信息的剩余壽命預(yù)測方法，得到針對未來載荷未知條件下的基于多傳感器信息的變載荷下剩余壽命預(yù)測方法。

1 隨機(jī)載荷建模

1.1 隱半馬爾可夫模型

在常規(guī)的馬爾可夫模型（Markov model，MM）中，模型在某狀態(tài)停留一定時間的概率為非顯式的概率模型，也就是說，沒有參數(shù)化的表達(dá)形式。HSMM通過在隱馬爾可夫（hidden Markov model，HMM）中添加顯式時間組件（temporal component）來改進(jìn)HMM模型中不準(zhǔn)確的時間建模。研究人員提出，用精心選擇的概率函數(shù)替換持續(xù)時間概率函數(shù)［12］，而具有這種精心選擇的狀態(tài)時間概率函數(shù)的HMM模型成為HSMM模型。

HMM過 程 由π、A、B等 元 素 構(gòu) 成［13］。與HMM模型不同，HSMM模型會在一個狀態(tài)中產(chǎn)生一系列的觀測值，HSMM模型參數(shù)可以記為

式中：π為初始狀態(tài)分布；A為狀態(tài)轉(zhuǎn)移模型；B為觀測概率模型；St表示t時刻下的隱藏狀態(tài)；O為觀測序列；D為狀態(tài)停留時間分布。

1.2 基于HSMM的隨機(jī)載荷空間建模

系統(tǒng)在工作中，常常需要經(jīng)歷多個工作狀態(tài)，如發(fā)動機(jī)在不同的馬赫數(shù)、油門解算器角度、飛行高度時，會對應(yīng)不同的工作狀態(tài)，即變載荷狀態(tài)空間。在載荷隨機(jī)變化的過程中，可以將其視為具有N個不同序貫狀態(tài)的過程。令di表示載荷i下的停留時間，T表示該系統(tǒng)歷經(jīng)所有載荷的總時長，則有

假設(shè)HSMM模型中的狀態(tài)序列χ包含了N個部分，令qi為一個部分中結(jié)束的時間點，則該狀態(tài)序列可以表示為表1。載荷空間變化過程如圖1所示。

表1 HSMM模型特征量Table 1 Characteristic variables in HSMM

圖1 HSMM載荷空間變化示意圖Fig.1 Sketch map of load space change of HSMM

為實現(xiàn)上述基于HSMM模型的隨機(jī)變載荷空間模型，對變載荷空間進(jìn)行預(yù)測，采用改進(jìn)的前向-后向變量與Viterbi算法，并通過動態(tài)規(guī)劃方法對模型進(jìn)行求解。

首先，將時間t時狀態(tài)i下觀測狀態(tài)o1，o2，…，ot的聯(lián)合概率分布定義為αt（i），即前向變量。同樣，可以推出后向變量βt（i）。前-后向變量計算過程如圖2所示。對于HSMM模型λ：

圖2 前向-后向變量算法示意圖Fig.2 Sketch map of forward-backward variable algorithm

O在給定HSMM模型的情況下的概率模型可以記為

為對模型的其他參數(shù)進(jìn)行求解，繼續(xù)定義了3個變量：第1個變量為載荷處于狀態(tài)i時停留時間di，再跳轉(zhuǎn)到狀態(tài)j的平均概率值?t，t′（i，j）；第2個變量為時刻t處于狀態(tài)i和時刻t′處于狀態(tài)j時的部分觀測值o1，o2，…，ot′的聯(lián)合概率αt，t′（i，j）；第3個 變 量 為 在 給 定 觀 測 值 序 列o1，o2，…，oT載荷處于狀態(tài)i時的停留時間為di=t′-t，再跳轉(zhuǎn)到狀態(tài)j的概率值ξt，t′（i，j）。

Viterbi算法的流程可以參考文獻(xiàn)［12］。令λ表示部分HSMM過程的參數(shù)集，則對于觀測序列O與對應(yīng)的狀態(tài)序列Z可以寫為

式中：θdi與θfi為對應(yīng)分布的參數(shù)集。

通過EM算法，可以對λ進(jìn)行求解，其求解過程可以總結(jié)為如下過程。

M步驟：設(shè)定新的λnew使得對數(shù)極大似然函數(shù)最大：

2 未知載荷條件下的剩余壽命預(yù)測

2.1 基于HSMM的隨機(jī)載荷預(yù)測

1）初始狀態(tài)分布

在已知O時，模型的初始狀態(tài)即在狀態(tài)i為初始狀態(tài)的概率：

2）狀態(tài)轉(zhuǎn)移概率

狀態(tài)轉(zhuǎn)移概率的估計為從狀態(tài)i轉(zhuǎn)移到狀態(tài)j的期望與從狀態(tài)i轉(zhuǎn)移的總期望值之間的比值：

3）觀測值分布

當(dāng)使用顯式狀態(tài)停留時間模型時，狀態(tài)停留時間密度反映在觀測值中的重要性在于可以顯著提高模型的質(zhì)量。然而，使用顯式狀態(tài)停留模型當(dāng)作狀態(tài)停留時間密度有一定的缺點：①會使與狀態(tài)停留時間相關(guān)的變量計算量大增；②包含了除基本的HMM模型參數(shù)外的與每個狀態(tài)都相關(guān)的大量參數(shù)d的模型。因此，HSMM模型的參數(shù)估計過程比標(biāo)準(zhǔn)的HMM過程更加困難。解決問題的方案是使用參數(shù)化狀態(tài)停留時間密度函數(shù)而不是非參數(shù)密度函數(shù)。同時，綜合考慮高斯分布、威布爾分布等分布的應(yīng)用情況和在訓(xùn)練時計算量的大小，本文采用高斯分布對狀態(tài)停留時間概率進(jìn)行建模。狀態(tài)停留時間參數(shù)通過訓(xùn)練集進(jìn)行估計，狀態(tài)停留時間的高斯分布的均值與方差可通過如下步驟進(jìn)行估計。

1）訓(xùn)練HSMM模型識別N個不同載荷，即建立對載荷進(jìn)行分類的模型。已知N組觀測序列，對N組不同的HSMM模型（如HSMM1，HSMM2，…，HSMMN）中的每一組都進(jìn)行特征建模。

2）對HSMM模型進(jìn)行分類。在給定觀測序列時，N組HSMM模型中的每一個都通過同一序列呈現(xiàn)。根據(jù)極大似然估計，可以對這些序列進(jìn)行分類。

式中：D（Zi）表示載荷Zi的狀態(tài)停留時間分布；μ（Zi）和σ2（Zi）分別為均值和方差；ρ為相應(yīng)的系數(shù)。

具體的狀態(tài)停留時間模型預(yù)測步驟如下：

步驟1 通過HSMM訓(xùn)練程序（如參數(shù)估計），從而可以獲取HSMM的狀態(tài)轉(zhuǎn)移概率。

步驟2 通過HSMM的參數(shù)估計，每一個載荷下的狀態(tài)停留時間的概率分布函數(shù)可以獲取。

步驟3 通過分類，可以識別實時的載荷。

圖3 基于HSMM的隨機(jī)載荷預(yù)測Fig.3 Random load prediction based on HSMM

2.2 基于隨機(jī)預(yù)測載荷的剩余壽命預(yù)測

基于本文提出的隨機(jī)載荷預(yù)測方法，可以得到每一個HSMM模型中每一個狀態(tài)的停留時間分布D及狀態(tài)之間的轉(zhuǎn)移概率A。基于所求解的HSMM模型，已知tk-1時的狀態(tài)概率Zi，遞推在時間tk時狀態(tài)為Zj的概率為

基于這一遞推公式，可以得到每一時刻的載荷與處于該載荷的概率。將這些概率代入文獻(xiàn)［14］中所述的預(yù)測壽命的方法，便可計算得出預(yù)測的剩余壽命：

式中：yi為系統(tǒng)/部件的失效時間；Ω（·）為累積分布函數(shù)；σ為尺度參數(shù)；π（si，p，l（t））為位置參數(shù)；si，p，l（t）為載荷l下第i個系統(tǒng)第p個傳感器的退化信號。

系統(tǒng)剩余壽命為式（14）得到的壽命估計值減去實時的時間值。基于隨機(jī)預(yù)測載荷的剩余壽命預(yù)測流程如圖4所示。

圖4 基于隨機(jī)預(yù)測載荷的剩余壽命預(yù)測Fig.4 Remaining life prediction based on random load prediction

得到剩余壽命后，對其進(jìn)行誤差分析來驗證本文方法的準(zhǔn)確性。剩余壽命相對誤差計算方法如下：

式中：Emean為均值的相對誤差；Evar為方差的相對誤差；meanes和meansim分別為預(yù)測得到的剩余壽命分布的均值和原始數(shù)據(jù)庫中的剩余壽命分布均值；Vares和Varsim分別為預(yù)測得到的剩余壽命分布的方差和原始數(shù)據(jù)庫中的剩余壽命分布的方差。

3 發(fā)動機(jī)應(yīng)用案例

商用模塊化航空推進(jìn)系統(tǒng)仿真（C-MAPSS）數(shù)據(jù)庫的多傳感器退化數(shù)據(jù)為基于物理退化模型的飛機(jī)渦扇發(fā)動機(jī)仿真數(shù)據(jù)［15］。本文選取CMAPSS數(shù)據(jù)庫中發(fā)動機(jī)在多種載荷條件下的性能退化數(shù)據(jù)進(jìn)行應(yīng)用案例驗證。圖5展示了CMAPSS系統(tǒng)中的發(fā)動機(jī)主要部件。所使用的數(shù)據(jù)基于如下假設(shè)：發(fā)動機(jī)退化是由于單個部件（如高壓腔與高壓壓氣機(jī)）在使用過程中特定載荷下的磨損引起的。同時，在數(shù)據(jù)庫中，特定飛機(jī)上的損傷累計是無法直接基于飛行時間與飛行載荷計算的。因此，只能使用飛行過程與飛行剛剛結(jié)束時多傳感器記錄的退化信號對飛機(jī)渦扇發(fā)動機(jī)的剩余壽命進(jìn)行計算。

圖5 商用渦扇發(fā)動機(jī)結(jié)構(gòu)簡圖Fig.5 Structure of commercial turbofan engine

該仿真數(shù)據(jù)庫包含260個樣本，挑選其中的80%（208個）作為歷史數(shù)據(jù)，剩余20%（52個）作為在線測試樣本。其中，包含了21個性能傳感器參數(shù)，同時包含了3個表征載荷狀態(tài)的傳感器。載荷狀態(tài)空間如表2所示，其不同載荷主要由3個指標(biāo)表征，即馬赫數(shù)、高度、油門解算器角度（TRA）。從表2中可知，這3個表征載荷的指標(biāo)綜合構(gòu)成6個不同的載荷狀態(tài)，該仿真數(shù)據(jù)在這6種載荷狀態(tài)下仿真得來的。由于渦扇發(fā)動機(jī)整體結(jié)構(gòu)十分復(fù)雜，難以建立準(zhǔn)確的物理模型，采用基于數(shù)據(jù)的傳感器融合方法進(jìn)行處理，并基于訓(xùn)練集數(shù)據(jù)對測試集樣本進(jìn)行未來載荷未知情況下的剩余壽命預(yù)測研究。

表2 載荷狀態(tài)空間Table 2 Space of load states

圖6展示了6個不同載荷條件隨時間變化的過程?？梢钥闯?，在渦扇發(fā)動機(jī)模擬仿真過程中，6個應(yīng)力狀態(tài)隨時間變化，且綜合考慮數(shù)據(jù)庫中所有的數(shù)據(jù)，載荷變化呈現(xiàn)出一定的隨機(jī)性。也就是說，3個表征載荷狀態(tài)的參數(shù)呈現(xiàn)出隨時間變化的特點且具有一定的隨機(jī)性，同時，可以認(rèn)為后續(xù)狀態(tài)只與前一個狀態(tài)相關(guān)而與之前的其他狀態(tài)無關(guān)，即具備馬爾可夫空間的無后效性。同時，表征3個指標(biāo)的傳感器可以視為觀測狀態(tài)，而不同的觀測狀態(tài)的組合構(gòu)成了載荷狀態(tài)，因此可以根據(jù)基于HSMM模型的剩余壽命預(yù)測方法對CMAPSS仿真數(shù)據(jù)中的隨機(jī)載荷進(jìn)行建模，將這6個載荷空間中不同載荷狀態(tài)下的停留時間視為一個隨機(jī)分布過程，并使用高斯過程對狀態(tài)停留時間進(jìn)行擬合。從而得到基于HSMM模型的載荷空間構(gòu)建，通過2.1節(jié)中提出的模型參數(shù)方法對載荷空間的參數(shù)進(jìn)行求解，從而在未來載荷未知的情況下對載荷變化進(jìn)行預(yù)測，為未來載荷未知時的渦扇發(fā)動機(jī)剩余壽命預(yù)測提供基礎(chǔ)。

圖6 六個載荷狀態(tài)隨時間變化示意圖Fig.6 Sketch map of six load states changing with time

首先，基于2.1節(jié)提出的基于HSMM模型的變載荷空間建模與參數(shù)預(yù)測方法，可以求解得到HSMM模型的相應(yīng)參數(shù)，如表3和表4所示。假設(shè)初始狀態(tài)為狀態(tài)Ⅴ，未來載荷預(yù)測結(jié)果如表5所示?；谶@一載荷模型可以對未來載荷的變化情況進(jìn)行預(yù)測，結(jié)合文獻(xiàn)［14］中提出的預(yù)測方法與式（13），可以對載荷未知的條件下的剩余壽命進(jìn)行預(yù)測?；贖SMM模型的變載荷空間建模與參數(shù)預(yù)測方法是基于過去已知載荷預(yù)測未來載荷并代入壽命預(yù)測模型預(yù)測壽命，案例中同時也使用發(fā)動機(jī)從開始到失效真實經(jīng)歷的載荷代入文獻(xiàn)［14］的方法進(jìn)行壽命預(yù)測，即未來載荷已知的情況。

表3 六個載荷狀態(tài)之間的轉(zhuǎn)移概率矩陣Table 3 Transition probability matrix of six load states

表4 六個載荷狀態(tài)各自停留的時間均值與方差Table 4 Average and variance of settling time under six load states

表5 未知載荷預(yù)測結(jié)果Table 5 Pr edicted results for unknown loads

使用不同使用信息占比下的部分剩余壽命預(yù)測結(jié)果如表6與表7所示。為了分析預(yù)測結(jié)果的準(zhǔn)確性，再根據(jù)式（15）和式（16）計算剩余壽命相對誤差的均值與方差值，如表8與表9所示。

表6 未來載荷已知時的剩余壽命Table 6 Remaining life when future loads are known

表7 未來載荷未知時的剩余壽命Table 7 Remaining life when future loads are unknown

表8 剩余壽命相對誤差均值Table 8 Average relative error of remaining life

表9 剩余壽命相對誤差方差Table 9 Variance of relative er ror of remaining life

從表8和表9可以看出，預(yù)測模型達(dá)到了一定的準(zhǔn)確度，即對均值與方差的估計值與載荷已知條件下的誤差相差小于5%。因此，使用HSMM模型對隨機(jī)載荷空間進(jìn)行建模與預(yù)測具有一定的有效性。通過表8和表9可以看出，在使用數(shù)據(jù)占比越大時，對剩余壽命的均值與方差估計值的估計誤差越小。這是由于大量的數(shù)據(jù)提高了模型參數(shù)估計值精度提高，從而提高整個模型的預(yù)測精度。

4 結(jié) 論

針對運行在多工況下的機(jī)械系統(tǒng)，本文提出的基于隱半馬爾可夫模型的壽命預(yù)測方法解決了未來載荷未知情況下的壽命預(yù)測問題。通過案例應(yīng)用，驗證了本文提出的預(yù)測方法的有效性和正確性。

1）相比于當(dāng)前壽命預(yù)測方法僅考慮單一工況或者系統(tǒng)未來載荷與歷史載荷服從獨立同分布的假設(shè)，本文將預(yù)測的未來載荷信息代入壽命預(yù)測模型進(jìn)行壽命預(yù)測，更符合實際情況。

2）與發(fā)動機(jī)真實載荷代入計算的結(jié)果進(jìn)行比較，結(jié)果表明本文提出的基于隨機(jī)預(yù)測載荷的剩余壽命預(yù)測方法具有一定的準(zhǔn)確性。由于載荷預(yù)測模型引入了預(yù)測誤差，這種情況下的剩余壽命預(yù)測值比載荷已知時的剩余壽命預(yù)測值的相對誤差更大，但是依然保持在較低的水平。