基于HDP-HMM的機械設備故障預測方法研究

王 恒，周易文，瞿家明，季 云

(南通大學 機械工程學院，江蘇 南通 226019)

故障預測與健康管理(Prognostics and Health Management,PHM)是近年來提出的一種集故障診斷、故障預測和健康管理能力與一體的新型系統[1-3]，它借助于信息化、智能化手段對關鍵部位的故障進行實時檢測與隔離，并預測裝備關鍵部位的剩余壽命。PHM系統的構建，必須以具有高準確率的故障預測技術為基礎。由于機械設備本身結構和機理的復雜性，影響系統運行的因素復雜多變，機器所表現出的狀態行為具有大范圍不確定性、高度非線性、動態時變性、強關聯性等特征，從而造成通常難以建立精確、完備的機理模型、難以描述系統各部分之間的依賴關系，給精確預測帶來更大的困難和挑戰。

隱馬爾科夫模型(Hidden Markov Model,HMM)具有狀態隱含、觀測序列可見的雙重隨機屬性，很好地描述了設備運行過程中的衰退狀態與觀測到的征兆信號(如振動、轉速和位移等)之間的隨機關系，在機械設備性能退化評估與預測中得到廣泛應用[4-7]。但是，HMM的定義和學習過程中還亟待解決以下問題：① 需要預設設備所經歷的狀態數，設備所經歷的狀態數必須預先設定，這就需要對模型先驗知識比較了解，但在很多實際問題中，缺乏相應的先驗知識，并不能準確地給出模型的退化狀態數，限制了HMM應用的場合。曾慶虎[8]提出了基于最小描述長度(Minimum Description Length,MDL)算法，通過MDL準則優化調整狀態數，但需要預先確定狀態數、計算過程較為復雜；騰紅智等[9]提出了基于K均值算法和交叉驗證相結合的狀態數優化方法；張星輝等[10]建立了一組聚類方法評價指標，利用K均值算法對狀態特征進行聚類，通過指標評定結果從中選取模型的最優狀態數，但K均值聚類仍需要預先確定狀態數。② 在HMM模型訓練中采用極大似然估計算法，沒有考慮模型的復雜度，容易造成過擬合和欠擬合等問題。

針對HMM存在的不足，本文引入了分層狄利克雷過程—隱馬爾科夫模型(Hierarchical Dirichlet Process-HMM,HDP-HMM)，利用HDP模型良好的聚類特性和分層共享原理，優化了HMM模型結構；在此基礎上，提出了基于HDP-HMM的機械設備故障預測算法，利用HDP-HMM所建立的狀態動態轉移關系，實現設備的健康等級評估和故障預測。滾動軸承軸承性能退化評估與故障預測結果表明了本文所提方法的有效性和適用性。

1 基于HMM的機械設備退化過程描述

設備的全壽命周期一般可以看作多狀態之間循序漸進的轉移過程，由正常狀態經過一系列不同程度的退化狀態，最終到達功能性故障。設備漸進性故障演化過程如圖1所示。

圖1 機械設備漸進性故障演化過程示意圖Fig.1 Schematic diagram of mechanic equipment gradual fault evolution process

當設備運行到“早期故障點P”這一臨界點，也稱變點或異常點[11]，其性能開始緩慢地發生退化，健康指數逐漸下降，直到運行至“設備功能故障點F”，這段時間是存在異常行為的設備性能退化過程。設備剩余壽命是指由監測到設備當前狀態點時刻到設備功能故障點的時間間隔。根據設備當前所處的退化狀態，評估健康狀況、預測剩余使用壽命，對于視情維修與健康管理具有重要的意義。

設備的性能退化是一個輸出特性參數可見、狀態隱藏的隨機過程，可以用HMM描述。假設HMM描述的設備全壽命周期共有S={1,2,...,K}個隱藏狀態，{1}為設備的正常狀態，{2,3,...,K-1}分別為設備K-2個依次嚴重程度的退化狀態，{K}為設備故障狀態。到T時刻為止，產生的隱狀態序列為{s1,s2,...,sT}，{x1,x2,...,xT}是不同狀態下出現的觀測值，觀測值xt關于狀態st條件獨立。設備各時刻所屬的狀態可以從K個不同的狀態中取值，不同狀態之間的轉換由轉移概率矩陣決定，其元素為

πij=p(st=j|st-1=i)

(1)

當HMM的參數λ確定的情況下，離散狀態s和觀測序列x的聯合概率密度函數可表示為

(2)

HMM的模型定義和標準估計過程都受到嚴格的限制。首先，在HMM的定義中，狀態數K是需要提前確定的，而且是有限的，即狀態只能在K個狀態之間進行轉移。狀態數K的確定是HMM模型訓練和測試的關鍵，主要根據經驗人為設定，很難將各種類別都考慮到，而且新的數據中也可能有未知類型出現，依靠訓練樣本得到的固定模型結構對新的觀測數據的適用性、涵蓋性不強。其次，為增加HMM的魯棒性和泛化性，一般需要采用多個觀測樣本進行訓練[12]，但在訓練多組觀測值序列時，只能先訓練單一觀測序列，再綜合所有單個觀測值序列得到的訓練結果，計算時間長、計算量較大。

2 基于HDP-HMM的機械設備故障預測算法

為解決HMM在退化評估中存在的問題，提出了一種基于改進HMM的機械設備故障預測算法，引入HDP-HMM[13-14]，該模型的主要優點是：①在HMM-HDP中，狀態類數可以是不確定的，即狀態空間無限，實現了擴展有限空間狀態的隱馬爾科夫模型到無限維度，并利用狄利克雷過程(Dirichlet Process,DP)的性質，實現HMM狀態數自動生成；②HDP是在DP基礎上擴展基礎分布，使多個數據源之間不必滿足獨立同分布條件，利用參數共享，實現多維特征的信息融合和聚類。算法的流程如圖2所示。

該算法的主要步驟為：

步驟1構造HDP作為HMM參數的先驗分布，通過多組觀測數據計算模型的后驗概率，并對HMM的參數進行動態調整，確定設備狀態數K(默認1為正常狀態、K為故障狀態)；

圖2 基于HDP-HMM的機械設備故障預測算法框圖Fig.2 Flow chart of mechanic equipment prognostic algorithm based on HDP-HMM

步驟2采用不同狀態對應的觀測數據，對HMM進行參數估計，訓練出對應于每個狀態的HMMn(λn)模型，建立設備退化狀態識別庫；

步驟3利用設備全壽命數據，訓練涵蓋設備所有狀態的HMM模型，獲得每個退化狀態所持續時間，并確定設備早期故障點P和功能故障點F；

步驟4針對當前觀測序列，利用退化狀態識別庫計算P(O|λn)，識別設備當前所屬狀態，如果處于退化狀態，預測其剩余壽命；如果處于故障狀態就進行故障報警。

2.1 基于HDP-HMM的設備退化狀態數確定

2.1.1 HDP-HMM的構造

在HMM-HDP中，記t時刻的狀態為st，對于HMM中的每一個狀態，它下一個狀態由當前狀態的轉移概率密度函數決定。設πk為第k個狀態的轉移概率密度函數，則有st∶πst-1。在給定狀態st下,觀測值xt的分布為Fφst(xt)，此分布通常稱為發散分布，各個狀態下觀測值的分布組成發散矩陣。面向設備退化評估的HDP-HMM有向圖模型，如圖3所示。

圖3 面向設備退化評估的HDP-HMM有向圖模型Fig.3 Directed graphical model of HDP-HMM for equipment degradation assessment

HDP-HMM的截棍構造(Stick-breaking)為

β|γ∶GEM(γ)
πk|α,β∶DP(α,β)
φk|H∶H(λ′)

(3)

在時間t=1,2,L,T時，狀態之間的轉移分布和觀測序列分布分別為

(4)

2.1.2 HDP-HMM的采樣

馬爾科夫鏈-蒙特卡羅(Markov Chain Monte Carlo,MCMC)統計模擬方法能夠簡化條件概率的計算，Gibbs算法就屬于MCMC算法之一。Gibbs采樣就是用條件分布的抽樣來代替全概率分布的抽樣。本文中已知觀測數據條件分布，通過Gibbs對條件分布循環抽樣新的觀測數據，實現對隱狀態數的更新。為敘述方便，約定當某一變量的上角標或下角標有符號“”時，例如Zi(Zi)為對應的變量集中移出下標(上標)對應的變量，Zi為zi將從Z中移出后由剩余的數據組成的數據集，即Zi={z1,...,zi-1,zi+1,...,zN}。

(1)采樣β

(5)

(2)采樣st

p(st|S ,x1,…,xT,β,α,λ′)∝
p(st|S ,β,α)p(xt|X ,st,S ,λ′)

(6)

式(6)的第1部分，根據馬爾科夫鏈和Dirichlet過程的性質可得

1)當st=k，k為已出現過的狀態

(7)

(8)

式(6)的第2部分，p(xt|X ,st,S ,λ′)即為觀測數據在隱狀態中的分布概率。設觀測數據分布函數F的密度分布函數為f(·|θ)，基分布函數H的密度分布函數為h(·|λ′)，選取F和H為共軛分布，利用貝葉斯公式將分布參數消去后，可得觀測數據的條件分布

(9)

式中：xji為第j組第i個觀測數據；fkxji為觀測數據x中除了xji外，屬于第k類的條件概率分布。

(3)采樣mjk

(10)

式中：mjk為第j組屬于第k類的聚類數；mjk除了第j組的第k類之外的聚類數；nj,k為第j組觀測數據中屬于k類的數據總數；s(nj,k,m)為Stirling數。

通過對β,m的交叉采樣實現參數更新，并在Stick-breaking構造中運用直接分配后驗采樣算法對觀測數據的指示因子st進行采樣，每次只更新一個數據的聚類屬性，當某個類簇中元素個數為0時，K-1，繼續迭代，待聚類數穩定時，停止迭代，獲得最終的隱狀態數K。

(4)更新超參數

模型中的超參數包括Θ={α,γ}，超參數決定了DP過程的離散程度，初始值的選取對聚類結果影響很大，這大大影響了HDP模型的通用性。因此，在進行后驗參數的更新時，需要同時對超參數進行更新，限于篇幅，超參數更新算法不再贅述。

2.2 基于HDP-HMM的設備退化評估與故障預測

(11)

(12)

(3)利用設備全壽命數據，訓練涵蓋設備所有狀態的HMM模型，得到退化狀態轉移路徑，獲得每個狀態的持續時間，并確定設備早期故障點P和功能故障點F。若判斷出設備處于退化狀態，按式(13)估計其剩余壽命。

(13)

式中：TRULn為設備處于第n個退化狀態時的剩余壽命；Tn為設備第n個退化狀態持續的時間。

3 實例分析

本文采用美國USFI/UCR智能維護系統中心的軸承加速壽命實驗進行實例分析，實驗裝置如圖4所示。四個ZA-2115雙列軸承并列安裝在同一軸上，轉速為2 000 r/min，在軸承上加載約26 671 N的徑向載荷，采樣頻率為20 kHz，每隔10 min采集一次數據，軸承1運行了約163 h后出現嚴重故障而失效。

圖4 軸承全壽命實驗裝置示意圖Fig.4 Schematic diagram of the bearing life device

3.1 HDP-HMM算法有效性分析

(1)多觀測序列聚類分析

選取軸承1的均方根值(Xrms)和峭度指標(Kurtosis)作為HDP-HMM的輸入觀測序列，假設基分布函數H的密度分布服從多項式分布、觀測數據分布的參數服從DP分布。隨機設定初始聚類數目K=50，集聚參數α=20，γ=1，迭代次數為M=300，通過Stick-breaking構造HDP過程，實現狀態數的自動聚類，結果如圖5(b)所示。圖5(a)為采用單一觀測序列(峭度指標)獲得的隱狀態數。對比圖5(a)和圖5(b)可知，單一和多觀測序列聚類結果均收斂于5。因此，HDP算法能夠有效實現多觀測序列組合聚類。

圖5 單一和多觀測序列下的隱狀態數Fig.5 Comparison of hidden states using single and multi-group observation data

(2)HDP-HMM參數設置分析

研究超參數α取值對HDP-HMM聚類結果的影響，分別取α=2，α=20，α=2000，其余參數不變，HDP-HMM聚類結果如圖6所示。結果表明，無論α的初始值如何選取，該模型聚類結果均能收斂到相同的值(超參數γ的情況類似)。由于采用了更新算法，HDP-HMM中超參數初始值可以任意選擇，使得算法具有很強的魯棒性。

3.2 基于HDP-HMM的機械設備性能退化評估

利用軸承的全壽命觀測序列進行訓練和建模，基于HDP-HMM的軸承性能退化狀態識別結果，如圖7所示。

由圖7中的峭度觀測序列曲線可知，軸承在觀測序列500～600波動較為平穩，而在觀測序列700～960波動忽高忽低，毫無規律可循，在觀測序列800附近峭度值回落到正常值。因此，僅從特征值指標判斷軸承的運行狀態容易引起很大的誤差，甚至導致故障發現不及時，造成嚴重事故。

圖6 不同α值的HDP-HMM聚類結果 Fig.6 Clustering results of HDP-HMM under different α

圖7 基于HDP-HMM的軸承退化狀態識別Fig.7 Identification of bearing degradation states based on HDP-HMM

由圖7的軸承運行狀態轉移曲線可知，從正常到失效的全壽命歷程中，一共出現了5次(即K=5)狀態轉移，分別代表正常狀態(1)、早期退化狀態(2)、中度退化狀態(3)、嚴重退化狀態(4)和故障狀態(5)。結合圖1可知，通過HDP-HMM模型有效地找出了軸承運行過程中的早期故障點P(576號文件)，便于預測剩余壽命和健康管理；同時能夠確定軸承的設備功能故障點F(930號文件)，方便及時預警，預防重大事故的發生。此外，在HMM-HDP算法中，馬爾科夫鏈只與前一個狀態有關，即t時刻系統狀態的概率分布只與t-1時刻的狀態有關，轉移概率隨時間變化不斷更新，因此狀態是否發生轉移是由轉移概率和發散概率綜合決定，與觀測序列的實際值關系不大。因此，HDP-HMM能識別軸承在運行過程中的一系列不同程度的退化狀態，較好地模擬了漸進性故障演化過程，為設備的早期維護和故障預測提供了理論指導。

采用均方根值觀測序列得到的軸承退化狀態識別結果，如圖8所示。對比圖7和圖8可知，采用峭度和均方根觀測序列所獲得軸承退化路徑基本上是一致的。但就早期故障點來說，峭度指標是概率密度分布尖峭程度的度量，對早期故障有較高的敏感性，其早期故障點相比于均方根值要早一些；有效值是對時間平均的，用來反映信號的能量大小，對早期故障不敏感，但穩定性比峭度指標要好，波動比峭度指標要小。

圖8 基于均方根值的軸承退化狀態識別Fig.8 Identification of bearing degradation states based on RMS

為了驗證HDP-HMM算法的有效性，與文獻[15]中基于K-S檢驗的軸承退化狀態識別方法進行對比分析。

圖9為基于K-S檢驗的軸承退化狀態識別結果，對比圖7和圖9可知，基于K-S檢驗的退化評估得到的軸承早期故障點P在第533序列處，設備功能故障點F在962處，與本文所提算法獲得的結果大體一致。但在軸承后期嚴重磨損階段，基于距離的K-S檢驗算法將振幅跳變劇烈處的狀態劃分過多，在實際工作過程中不利于對軸承狀態的判定。

圖9 基于K-S檢驗的軸承退化狀態識別Fig.9 Identification of bearing degradation states based on K-S test

3.3 基于HDP-HMM的設備故障預測

根據軸承在每個狀態所對應的持續時間，代入式(13)可得不同退化狀態下的剩余壽命，如表1所示。通過對軸承的健康分級和壽命預測，為基于狀態的設備維護和健康管理提供了依據。

表1 不同退化狀態下的軸承剩余壽命Tab.1 Residual life of different degradation states

4 結 論

(1)提出了一種基于HDP-HMM的機械設備故障預測算法，該算法既有HMM 處理序列數據的特點，又具有HDP自動生成聚類數目的功能，獲得了設備退化狀態數，解決了傳統的HMM模型狀態數必須預先設定的不足。

(2)建立了設備不同狀態下和全壽命歷程對應的HMM模型，利用每個退化狀態所對應的持續時間，實現了對設備不同退化狀態下的剩余壽命預測。

(3)滾動軸承全壽命數據建模結果表明，HDP-HMM模型可以有效地找出軸承運行的早期故障點和功能故障點，并能夠識別軸承在運行過程中的一系列不同程度的退化狀態，與基于距離聚類的K-S檢驗算法相比，HDP-HMM更能有效描述設備實際退化過程。