劉文溢，劉勤明，周林森

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

0 引言

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

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

為了更好地參考歷史統(tǒng)計信息以提高模型識別率，研究者們提出了高階隱馬爾可夫模型(Higher-order HMM，HOHMM)。針對高階模型隨之而來的參數(shù)爆炸問題以及更復雜的推導問題，研究者們建立了一種逐次降階算法(ORED算法)來簡化HOHMM，在此基礎(chǔ)上提出了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]建立了日均氣溫演化模型，并假設(shè)非對稱分量為一個高階隱馬OU(Ornstein-Uhlenbeck)過程，結(jié)果表明該模型能有效地捕捉各種情況下的溫度數(shù)據(jù)的特征

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

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

1 隱馬爾科夫模型基礎(chǔ)

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)計信息。其假設(shè)所研究對象當前狀態(tài)與之前數(shù)個狀態(tài)有關(guān)。以一個n階隱馬爾科夫模型為例，其時刻t(t≥n+1)時的狀態(tài)與前n個時刻的狀態(tài)均有關(guān)，表示為：

Prob(Statet|State1…Statet-1)=

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

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

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

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

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

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

其中圖2a為在二階HOHSMM模型上對模型的一個劃分方式，將相鄰狀態(tài)節(jié)點合并為一個更大的“新節(jié)點”，此時新節(jié)點內(nèi)部兩個狀態(tài)的關(guān)系對整個模型的馬爾科夫性不產(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)移過程路徑生成器

在此基礎(chǔ)上，給定模型參數(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時對應索引集合，索引對應關(guān)系如表3所示。

表3 索引對應關(guān)系

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

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

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

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

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

(30)

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

4 算例分析

本文通過美國卡特彼勒公司液壓泵的設(shè)備健康診斷與壽命預測實例來驗證評價本文提出的模型與方法。實驗室中設(shè)備的振動信號是由安裝在與液壓泵旋轉(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)為設(shè)備的徹底失效狀態(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算法采用定慣權(quán)重策略，慣性權(quán)重取值為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é)束語

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

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