許丹，肖小琦，馮至昕

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

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

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

1 隨機(jī)載荷建模

1.1 隱半馬爾可夫模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3）觀測(cè)值分布

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

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

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

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

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

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

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

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

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

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

基于本文提出的隨機(jī)載荷預(yù)測(cè)方法，可以得到每一個(gè)HSMM模型中每一個(gè)狀態(tài)的停留時(shí)間分布D及狀態(tài)之間的轉(zhuǎn)移概率A?；谒蠼獾腍SMM模型，已知tk-1時(shí)的狀態(tài)概率Zi，遞推在時(shí)間tk時(shí)狀態(tài)為Zj的概率為

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

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

系統(tǒng)剩余壽命為式（14）得到的壽命估計(jì)值減去實(shí)時(shí)的時(shí)間值?；陔S機(jī)預(yù)測(cè)載荷的剩余壽命預(yù)測(cè)流程如圖4所示。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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