基于高階隱半馬爾科夫模型的設備剩余壽命預測

劉文溢，劉勤明，周林森

(上海理工大學 管理學院，上海 200093)

0 引言

隱馬爾可夫模型(Hidden Markov Model，HMM)作為一種概率統(tǒng)計方法，具有良好的隨機性表征能力以及潛在結(jié)構關系描述能力，被廣泛用于復雜系統(tǒng)建模領域。BUNKS等[1]首先將HMM從語音識別領域遷移應用于設備故障診斷領域，開創(chuàng)了故障診斷的全新研究方向；馬倫等[2]從理論上闡明了HMM必須強制服從指數(shù)分布才能滿足其馬爾科夫性的不合理性，從而不能直接用于設備的壽命預測。

為了擺脫HMM本身特性的束縛，研究者們提出了改進的隱半馬爾科夫模型(Hidden Semi Markov Model, HSMM)，將宏觀狀態(tài)進一步分為微觀狀態(tài)，并對每個宏觀狀態(tài)的時間分布做出先驗假設，從而成功實現(xiàn)了隱馬氏模型在剩余壽命預測領域的應用。董明等[3]假設HSMM的狀態(tài)駐留為顯式Gauss分布，并對設備進行了健康診斷與剩余壽命預測；LIU等[4]提出了信息融合的自適應隱馬爾科模型，完成了對設備的診斷與壽命的預測，取得了良好的效果；李永朋等[5]等將Erlang分布引入HSMM模型，并對模型進行了重新推導，提高了識別率，降低了剩余使用壽命(Residual Useful Life,RUL)誤差；原媛等[6]將HSMM輸出的狀態(tài)序列與比例故障率模型相結(jié)合，并假設基礎故障率函數(shù)為兩參數(shù)的威布爾分布。但上述研究均對模型關于時間的問題進行了先驗分布假設，然而特定設備的先驗知識通常難以獲取，從而不可避免地帶入了主觀因素，影響了客觀實際。

為了更好地參考歷史統(tǒng)計信息以提高模型識別率，研究者們提出了高階隱馬爾可夫模型(Higher-order HMM，HOHMM)。針對高階模型隨之而來的參數(shù)爆炸問題以及更復雜的推導問題，研究者們建立了一種逐次降階算法(ORED算法)來簡化HOHMM，在此基礎上提出了FIT算法(fast incremental algorithm)來訓練HOHMM模型，并對其三問題分別進行了探討[7]；相比于ORED算法的逐次降階，HARDER等[8]提出了一種基于等價變換思想的HARDER算法，將任意高階模型轉(zhuǎn)化為對應的一階模型。上述研究者為高階模型的研究提供了一般研究方法，意義重大，但HOHMM在降階為HMM模型過程中，都需要對模型的各個變量做相應的重新推導，并基于Baum-Welch算法對參數(shù)進行重估，且隨著模型階數(shù)上升，工作量也會呈現(xiàn)爆炸式的增長。ZHU等[9]建立了HO-HMMAR模型，進行了變量推導，更好地解決了最優(yōu)投資組合問題；XIHONG等[10]建立了日均氣溫演化模型，并假設非對稱分量為一個高階隱馬OU(Ornstein-Uhlenbeck)過程，結(jié)果表明該模型能有效地捕捉各種情況下的溫度數(shù)據(jù)的特征

多項式擬合作為一種非線性擬合的簡單工具，具有高效的數(shù)據(jù)處理能力。劉慧婷等[11]運用多項式擬合解決了EMD處理上下包絡過程的端點問題，取得了良好的效果；VYAS等[12]運用多項式擬合的系數(shù)產(chǎn)生虹膜識別的特征向量，并通過基準IITD數(shù)據(jù)集和CASIA-v4區(qū)間虹膜數(shù)據(jù)驗證了多項式方法的有效性。隨著神經(jīng)網(wǎng)絡以及各種機器學習方法的普及，多項式擬合簡單易用的閃光點在故障診斷領域愈發(fā)暗淡。但值得注意的是，多項式擬合是為數(shù)不多的能直接寫出非線性關系解析式的方法，而這是神經(jīng)網(wǎng)絡做不到的。

參數(shù)估計問題一直是概率類模型最棘手的環(huán)節(jié)，隨著工作量較大的最大期望(Expectation Maximization，EM)算法依賴于初值而易陷入局部最優(yōu)的不足漸漸呈現(xiàn)，越來越多的研究者將目光轉(zhuǎn)向各種智能優(yōu)化算法。WU等[13]將遺傳算法與灰色模型相結(jié)合，數(shù)值結(jié)果顯示，結(jié)合遺傳算法的新方法可極大地提高模型預測精度；KRAUSE等[14]采用基于遺傳算法的函數(shù)逼近器估計感應電機參數(shù)，取得了較其他方法更好的效果；ZHANG等[15]運用人工免疫算法對HMM模型進行優(yōu)化，得到了辨識度最高的初始觀測矩陣；張小強等[16]運用自適應基因粒子群算法優(yōu)化隱馬爾科夫模型，得到了比傳統(tǒng)HSMM更精確的結(jié)果。如何巧妙運用智能仿真算法簡化隱馬氏模型的研究過程，是未來研究中應當被重視的一點。

為了實現(xiàn)在分布未知前提下更加精確實際的預測，本文在上述研究的基礎上，首先建立了基于HSMM的高階隱半馬爾科夫模型(Higher-order HSMM, HOHSMM)框架，在ORED算法和HARDER算法等價變換的啟示下，提出一種基于排列的模型降階方法，利用高階隱馬類模型的定義，使其可轉(zhuǎn)化為相應的一階模型，并讓低階模型三問題的解決方案能夠用于高階復雜模型。同時，對轉(zhuǎn)移概率矩陣與觀測概率矩陣進行相應的變形，使得高階復雜模型中節(jié)點的相互依賴關系信息自然融入到模型參數(shù)之中，達到了簡化模型的效果。其次，定義推導了輔助駐留變量，并利用智能優(yōu)化算法群對攜帶了更多“依賴關系信息”的參數(shù)組進行估計，以極大化觀測出現(xiàn)概率為目標，依靠參數(shù)組表征高階模型的分解依賴關系，在一定程度上將復雜度從模型本身轉(zhuǎn)移至參數(shù)組。緊接著，運用多項式擬合方法擬合各駐留變量序列，在分布未知的情況下完成了對設備剩余壽命的預測。最后，以美國卡特彼勒公司液壓泵數(shù)據(jù)集對模型進行驗證評價，結(jié)果顯示本文方法是可行且有效的。

1 隱馬爾科夫模型基礎

1.1 一階隱馬爾科夫模型

一階隱馬爾科夫模型，即常規(guī)隱馬爾科夫模型(HMM)，可描述為λ=(π,A,B)。由可觀測的節(jié)點與隱狀態(tài)節(jié)點組成，若令Statet為時刻t時節(jié)點所處的狀態(tài)，隱狀態(tài)節(jié)點狀態(tài)之間的轉(zhuǎn)移符合馬爾科夫性即：

Prob(Statet|State1…Statet-1)=

Prob(Statet|Statet-1)。

1.2 高階隱馬爾可夫模型

高階隱半馬爾科夫模型(HOHMM)是一階馬爾科夫模型的推廣，保留了更多的歷史統(tǒng)計信息。其假設所研究對象當前狀態(tài)與之前數(shù)個狀態(tài)有關。以一個n階隱馬爾科夫模型為例，其時刻t(t≥n+1)時的狀態(tài)與前n個時刻的狀態(tài)均有關，表示為：

Prob(Statet|State1…Statet-1)=

Prob(Statet|Statet-n…Statet-1)。

HOHMM模型雖然在保留更多歷史統(tǒng)計信息的基礎上提高了模型識別率，但其依然未能克服HMM模型自身的不足。

2 改進高階隱半馬爾科夫模型

本文在綜合考慮HMM不足以及高階建模優(yōu)點的前提下，提出一種基于HSMM的改進高階隱半馬爾科夫模型(HOHSMM)。

2.1 基于排列映射的模型降階

同樣以一般意義上的二階HOHSMM為例，由于二階模型在結(jié)構上的變動，使得模型參數(shù)與相關算法隨之產(chǎn)生變化。常規(guī)意義上，低階模型三問題各自的解決算法并不適用于高階模型，對此，本文提出一種基于組合映射的模型降階法，實質(zhì)是將二階模型中相鄰的兩個時點所對應的隱狀態(tài)結(jié)點合并為一個節(jié)點，并對合并后的節(jié)點進行馬爾科夫過程建模，如圖2所示。

其中圖2a為在二階HOHSMM模型上對模型的一個劃分方式，將相鄰狀態(tài)節(jié)點合并為一個更大的“新節(jié)點”，此時新節(jié)點內(nèi)部兩個狀態(tài)的關系對整個模型的馬爾科夫性不產(chǎn)生影響，降階后的新模型的馬爾科夫性可由式(1)描述：

Prob(Statet,Statet-1|Statet-1…State1)=
Prob(Statet,Statet-1|Statet-1…Statet-2)(t≥3)。

(1)

(2)

(3)

式中j為狀態(tài)排列(i,j)中的第2個狀態(tài)，以N=4的二階HOHSMM為例，其降階后的概率轉(zhuǎn)移矩陣有效數(shù)據(jù)為20個，由表1稀疏表示給出，其中1和1*表示的位置為有效數(shù)據(jù)，1*所在位置列頭代表的復合狀態(tài)為主狀態(tài)，1所在位置列頭代表的復合狀態(tài)定義為過渡狀態(tài)，0表示的位置為無效數(shù)據(jù)。

表1 狀態(tài)數(shù)為4的二階HOHSMM降階轉(zhuǎn)移概率矩陣稀疏表示

2.2 模型推理

借鑒前向-后向算法的思想，引入LT(linger time)機制并建立輔助變量ξt(i)，描述為式(4)：

ξt(i,d)=Prob(O[1:t],LT((i,j))=d|λ),t≥d。

(4)

意為在給定模型參數(shù)組λ下，截止t時刻產(chǎn)生觀測序列O[1:t]以及在當前狀態(tài)已駐留了d個時間的概率。

2.1節(jié)中給出了一般二階HOHSMM模型轉(zhuǎn)移概率矩陣的稀疏表示，并分別闡述了主狀態(tài)與過渡狀態(tài)的意義。不難得出，常規(guī)二階模型不同主狀態(tài)之間的轉(zhuǎn)移是一個漸變的過程，即

主狀態(tài)→(過渡狀態(tài)序列)→下一主狀態(tài)。

且該轉(zhuǎn)化過程至多需要j-i+2個時點來完整描述，以狀態(tài)主狀態(tài)(0,0)轉(zhuǎn)移至(3,3)為例，轉(zhuǎn)移過程需要5個時點來描述，如圖3所示。

(1)當i=0時，d=t，ξt(i)描述為:

ξt(0,d)=Prob(O[1:t],LT(0)=d|λ);

(5)

遞歸初值為：

(6)

內(nèi)轉(zhuǎn)移遞歸式為:

(7)

(2)當i=1時，ξt(i)描述為:

ξt(1,d)=Prob(O[1:t],LT(1)=d|λ)。

(8)

(3)當i=2時，承接狀態(tài)可為(0,0)以及(1,1)，ξt(i)描述為：

ξt(2,d)=Prob(O[1:t],LT(2)=d|λ)；

(10)

內(nèi)轉(zhuǎn)移遞歸式為：

(11)

情形1當承接狀態(tài)為(0,0)時，(0,0)～(2,2)需要4個時點刻畫，承接值則為:

(12)

相應的?t(0,2,Road)可表示為:

(13)

相應的供歸值為：

(14)

(15)

θt(1,2,Road)為:

(16)

供歸值為:

(17)

(4)當i=3時，承接狀態(tài)可為(0,0)、(1,1)以及(2,2)，ξt(i)描述為:

ξt(3,d)=Prob(O[1:t],LT(3)=d|λ);

(18)

內(nèi)轉(zhuǎn)移遞歸式為:

(19)

情形1當承接狀態(tài)為(0,0)時，(0,0)至(3,3)需要5個時點刻畫，相應的承接值為:

(20)

?t(0,3,Road)可表示為:

(21)

供歸值為相應承接值與中間值的乘積。

情形2當承接狀態(tài)為(1,1)時，(1,1)～(3,3)需要4個時點刻畫，相應的承接值則為:

(22)

相應的中間值?t(1,3,Road)為:

(23)

遞歸中值為:

(24)

情形3當承接狀態(tài)為(2,2)時，(2,2)～(3,3)需要3個時點刻畫，相應的承接值則為

(25)

相應的中間值?t(2,3,Road)為

(26)

遞歸中值為

(27)

則ξt(i,d)的一般遞歸式可由式(28)表示：

(28)

式中ct=t-d-lr+1，其中l(wèi)r為Road的長度。

根據(jù)上述對遞歸式的分類推導，將狀態(tài)較高時的輔助變量ξt分解到較其更低的輔助變量的表達形式，并在遞歸過程中，由計算機生成所有滿足發(fā)生主狀態(tài)轉(zhuǎn)移的轉(zhuǎn)移過程狀態(tài)路徑，遍歷所有可能的轉(zhuǎn)移狀態(tài)路徑便能計算得到任意ξt(i,d)，各主狀態(tài)所包含的轉(zhuǎn)移狀態(tài)路徑由表2給出，其中x不計入路徑長度。

表2 各主狀態(tài)轉(zhuǎn)移過程路徑生成器

在此基礎上，給定模型參數(shù)λ，產(chǎn)生觀測O的概率為：

定義輔助變量τt(index)，意為給定模型參數(shù)與觀測的前提下，時刻t處于Mapping-1(index)的概率，即

τt(index)=

Prob(Statet=Mapping-1(index)|λ,O)。

τt+1(index)=

(29)

其對應的低階模型中，時刻t時的處于單狀態(tài)i(i∈(1,2,3,4))的概率用高階模型表示為：

其中I為不同復合狀態(tài)下第二子狀態(tài)為i時對應索引集合，索引對應關系如表3所示。

表3 索引對應關系

二階段似然函數(shù)為：

運用不同智能仿真算法來優(yōu)化兩階段似然函數(shù)并進行比較，最終選取其中最優(yōu)的結(jié)果。

3 不確定分布下的剩余壽命預測

2.2節(jié)中對各個狀態(tài)駐留時間的遞歸推理，能得到設備在每個狀態(tài)能夠駐留時間以及產(chǎn)生觀測的聯(lián)合概率值。實際上，設備產(chǎn)生的觀測以及設備在單個狀態(tài)的駐留間存在著一定相關關系，該關系可由觀測、狀態(tài)轉(zhuǎn)移矩陣通過特定的模式進行描述。為便于得到設備在各個狀態(tài)駐留時間的邊緣概率，不妨假設設備觀測與駐留之間相互獨立，則對于上述求得的ξt(i,d)，計算其所對應觀測的產(chǎn)生概率Prob(O|λ)，則根據(jù)條件概率公式可得

使用各個時點所代表的時間來計算某一狀態(tài)產(chǎn)生該段駐留的概率，得到一組駐留-概率序列數(shù)。通常在已知離散數(shù)據(jù)點的情況下，假設先驗分布有利于研究數(shù)據(jù)的連續(xù)性特征，但如果錯誤假設數(shù)據(jù)分布，則擬合將存在極大誤差，并且會丟失數(shù)據(jù)原有的性質(zhì)。因此，本文采用一種基于多項式回歸的方法去擬合未知分布下的數(shù)據(jù)，多項式回歸方法可描述為：

(30)

其中D為各個狀態(tài)的最大持續(xù)時間，Prob(St=i)為時刻t時，設備處于狀態(tài)i的概率，Ri,Rj為離散數(shù)據(jù)連續(xù)化后的積分縮放系數(shù)。

4 算例分析

本文通過美國卡特彼勒公司液壓泵的設備健康診斷與壽命預測實例來驗證評價本文提出的模型與方法。實驗室中設備的振動信號是由安裝在與液壓泵旋轉(zhuǎn)軸平行位置的液壓加速計收集。在應用實例中，分別對液壓泵充入20、40、60與80 mg的微塵，并每隔10 min運用一個長度固定的時間窗采集一個約為1 min的振動信號(Pump6)。隨后使用10dB的小波將振動信號分為五層，得到數(shù)組高頻與低頻小波系數(shù)，將經(jīng)過降維后的小波系數(shù)作為DGHMM的輸入特征序列向量。整個實驗過程中，液壓泵的狀態(tài)可分為四種，分別為Baseline、Cont1、Cont2以及Cont3，其相應的代表性采樣振動信號時域振幅情況如圖4所示(Pump6)，其中Cont3狀態(tài)為設備的徹底失效狀態(tài)。整個實驗分析平臺為Python 3，平臺運行環(huán)境為Windows 10。

4.1 部分數(shù)據(jù)預覽

利用來自液壓加速計振動信號監(jiān)測數(shù)據(jù)，對液壓泵進行健康狀態(tài)診斷以及剩余壽命預測。部分經(jīng)過小波變換后的振動數(shù)據(jù)(Pump6)如表4所示。

表4 Pump6部分小波變換數(shù)據(jù)

4.2 智能優(yōu)化算法群參數(shù)估計

考慮最大期望算法對初值有較大依賴的不足，本文在模型框架下分別采用遺傳算法(Genetic Algorithm, GA)、粒子群優(yōu)化(Particle Swarm Algorithm, PSO)算法以及人工魚群算法(Artificial Fish Swarm Algorithm, AFSA)對同一組全觀測下的模型進行參數(shù)估計。最大迭代次數(shù)為300次，種群數(shù)(粒子群中為粒子數(shù)量)為30，其中GA采用大型參數(shù)自適應值編解碼策略[17]；PSO算法采用定慣權重策略，慣性權重取值為0.5，學習因子均取2，；AFSA感知野取值1，擁擠因子取值0.6。以最大化觀測出現(xiàn)概率為目標優(yōu)化模型，迭代過程及結(jié)果如圖5所示，在相同種群數(shù)量以及迭代次數(shù)前提下，遺傳算法以持續(xù)進化以及最大似然而具有較好效果，3種優(yōu)化算法各自所尋得的參數(shù)組的似然值如表5所示，得到降階后的概率轉(zhuǎn)移概率矩陣由表6給出。

表5 三種算法最大似然值

表6 最優(yōu)概率轉(zhuǎn)移矩陣稀疏表示

4.3 剩余壽命預測

基于模型推理，對各狀態(tài)駐留進行分析，在計算過程中運用全概率公式進行等價替換。將不同時點每個狀態(tài)不同駐留值按狀態(tài)分開得到各個狀態(tài)的駐留-概率序列，運用多項式回歸對序列進行擬合得到各自的解析式，在擬合過程中優(yōu)先選取對后續(xù)狀態(tài)剩余壽命預測共振作用較小的階數(shù)n(下稱優(yōu)先選取原則)。以最終時點4個不同狀態(tài)的多項式擬合為例(如圖6)，其中4個狀態(tài)進行擬合的最優(yōu)m值分別為：20、15、19、17，詳細的多項式回歸系數(shù)由表7給出。

表7 最終時點多項式擬合系數(shù)一覽

續(xù)表7

原則上概率的積分不允許負值出現(xiàn)，但多項式回歸呈現(xiàn)明顯的波動特性，即正負值存在一定的抵消作用，且不同狀態(tài)相同時點也存在共振，不難預知預測的剩余壽命值可能存在一定的波動特性。在實驗過程中，為了更好地進行多項式擬合，結(jié)合原始數(shù)據(jù)特點，對原始數(shù)據(jù)點進行了線性插值。圖6中紅色點為原始數(shù)據(jù)點，即自適應后的駐留產(chǎn)生概率值，淺藍色為插值點，黃色線條為多項式擬合線。

最終得到Pump24的剩余壽命預測情況如圖7所示，其中離散預測點已進行插值平滑處理。一號標記處顯示預測值在最終的時點處存在較大的偏差，對應圖6各子圖中多項式回歸式初期波動幅度較大的現(xiàn)象，但波動隨著時點的增加而逐漸變小，是優(yōu)先選取原則的結(jié)果。二號標記處出現(xiàn)了“提前損壞”的現(xiàn)象，對應圖6d中駐留提前降低的現(xiàn)象。

與低階HSMM結(jié)果相比較[4]，本文模型預測的剩余壽命結(jié)果如表8所示，從抽樣的時點來看，基于HOHSMM與多項式擬合的剩余壽命預測方法整體效果明顯優(yōu)于常規(guī)HSMM，算例表明，基于HOHSMM與多項式回歸的壽命預測方法是有效可行的。

表8 預測剩余壽命相對誤差分析

5 結(jié)束語

本文提出一種基于高階隱半馬爾科夫模型與多項式擬合的設備剩余壽命預測方法。其中排列組合的模型降階方法簡單直觀，巧妙地利用了高階隱馬類模型的定義，使得高階模型可通過轉(zhuǎn)換觀察角度的方法轉(zhuǎn)化為相應的一階模型，使低階模型三問題的解決方案能夠用于高階復雜模型。通過智能優(yōu)化算法群對本文模型進行了參數(shù)估計，將高階模型之中復雜的依賴關系信息轉(zhuǎn)移至變形后的參數(shù)組中，有效地簡化了模型，為研究該類模型提供了一個新的思路，最后基于多項式擬合的剩余壽命預測從結(jié)果來看是有效的。

未來的工作重點應是考慮如何求得歷史壽命與剩余壽命的聯(lián)合分布，以便于修正模型所計算的剩余壽命。