基于隱馬爾科夫模型的滾動軸承性能衰退評估

郝 剛，金 濤

(1.海軍工程大學 動力工程學院，湖北 武漢 430033；2.武漢城市職業學院 機電工程學院，湖北 武漢 430064)

0 引 言

滾動軸承是回轉機械重要組成部分。以軸承為核心的轉子系統健康狀態監測、故障診斷、剩余壽命預測與健康狀態評估等研究問題成為熱點。軸承失效輕則造成旋轉機械的停機，重則造成系統損壞，甚至引發災難性事故。因此滾動軸承性能衰退評估對于其預防性維修，避免轉子系統停機或事故具有指導意義。

軸承性能衰退過程指的是從健康運行狀態變遷至亞健康狀態的隨機過程。以典型的滾動軸承為例，在循環交變應力的作用下，滾動軸承的滾道可能出現壓痕、剝落、裂紋，保持架出現膠合、裂紋，滾子出現剝落等現象[1]。這些單因素或多因素早期失效的出現雖然不影響滾動軸承的正常運行，但是會出現振動、噪聲增大、軸承發熱等表征現象，意味著軸承進入亞健康狀態，即性能衰退。

滾動軸承性能衰退的研究熱點主要集中在特征參數提取問題、健康狀態監測問題和壽命模型構建及剩余壽命預測問題。圍繞相關熱點，P.K.KANKAR等[2]采集滾動軸承的振動信號，利用小波分解得到不同的調制頻率進行特征參數提取，使用支持向量機(SVM)對性能衰退進行分類識別取得較好效果；吳軍等[3]在小波分析的基礎上基于CEEMDAN的本征模能量特征提取，然后對不同模態分量IMF進行Hilbert-Huang變換，最后經過PCA和斯皮爾曼等級相關方法對信號進行融合后評估性能衰退的等級；史曉雪等[4]提出了一種自適應遺傳粒子濾波(AGPF)算法對滾動軸承性能衰退的趨勢進行預測從而得出評估結果。滾動軸承性能衰退原因復雜，特征多樣化、采樣信號噪聲大等原因使其性能衰退評估對信號處理的準確性提出較高要求，且依賴歷史數據。

隱馬爾科夫模型是在一般的馬爾科夫鏈的基礎上發展起來的隨機過程，它通過隱含和觀測兩個序列表征雙重隨機過程。隱馬爾科夫模型的狀態轉移矩陣可以用來描述與馬爾科夫鏈的各個狀態之間的相互關聯量化關系[5,6]。

滾動軸承性能衰退的評估模型中，隱馬爾科夫模型力圖通過狀態轉移矩陣去感知和量化滾動軸承的健康狀態[7-9]。筆者提出隱馬爾科夫模型的滾動軸承性能衰退評估方法，使用滾動軸承振動信號的均方根值(root mean square，RMS)作為隱馬爾科夫模型的輸入。該特征是振動信號時域分析中常用的特征參數，適合基于三軸加速度信號的在線快速計算，并能初步剔除奇異值。利用對滾動軸承壽命的衰退階段分割，建立性能衰退評估和分類模型。

1 隱馬爾科夫模型概述

隱馬爾科夫模型可以用五元組來表示，即λ=(N,M,π,A,B)，其中，N為隱馬爾科夫鏈的狀態數目；M為模型每個狀態下對應可能的觀測值數目；B為狀態轉移概率的數據集合。隱馬爾科夫模型所表示的隨機過程如圖1。

圖1 隱馬爾科夫模型

隱馬爾科夫模型中針對評估問題、解碼問題和學習問題分別提出了對應的前向-后向算法、Viterbi算法和Baum-Welch算法。對于滾動軸承性能衰退評估方法的應用中，首先提取一組測試數據作為樣本觀測序列，即訓練數據，利用Baum-Welch算法建立隱馬爾科夫模型，即完成學習問題；其次使用另外一組測試數據作為觀測序列，即測試數據，利用已經建立的隱馬爾科夫模型和Viterbi算法，計算最優狀態序列，解算出以建立隱馬爾科夫模型中的狀態轉移函數，即完成解碼問題；再次利用前向-后向算法對觀測序列值計算后驗概率分布，計算得出在依托樣本觀測序列建立隱馬爾科夫模型后的觀測序列概率分布，由此得出性能衰退評估結果，即完成評估問題，如圖2。

圖2 基于隱馬爾科夫的滾動軸承性能衰退評估模型

2 基于隱馬爾科夫模型的滾動軸承性能衰退評估

滾動軸承性能衰退過程是一個從量變到質變的漸變過程，也是一個隨機過程。這一系列狀態變遷過程的成因復雜多變，沒有規律可以遵循。對工況條件下滾動軸承性能衰退評估及健康狀態預測應該是業界關注的重點，其失效的概率估計問題也是難點[10,11]。

利用一組觀測序列值就能通過隱馬爾科夫模型計算得到相應的后驗概率分布，從而得知目前滾動軸承處于哪個健康狀態的概率值最大，這便是識別問題。隱馬爾科夫模型利用Forward-Backward算法來解決識別問題，即計算P(O|λ)的概率[12]。

前向算法中，首先定義前向變量為：

αt(i)=P(O1,O2,…,Ot,qt=θi|λ)

(1)

式中：i=1,2,3,4,5;t=1,2,…,T。

經初始化：

(2)

遞歸計算：

(3)

式中：1≤t≤T；1≤i,j≤N。

終結計算：

(4)

算法示意如圖3。

圖3 前向算法遞歸關系

該圖反映了在t時刻的狀態1≥i≥N是通過怎樣的途徑到達t+1時刻的狀態j。前向算法的思想是通過遞歸方法計算由時刻t的狀態i向t+1時刻狀態j轉移的所有途徑的概率，這些概率值之和為P(O|λ)。

與前向算法類似，后向變量的定義為：

βt(i)=P(Ot+1,Ot+2,…,OT|qT=θi,λ)

(5)

式中：i=1,2,3,4；t=1,2,…,T-1。

取βt(i)=1，后向算法的計算過程如下：

1)初始化：βt(i)=1,1≤i≤N。

2)遞歸計算：

(6)

式中：1≤i≤T；t=T-1,T-2,…,1。

3)終結計算：

(7)

后向算法初始化時對于所有的狀態i定義βT(i)=1。t+1時刻的狀態j向t時刻的狀態i過度的流程如圖4。

圖4 后向算法遞歸關系

通過建立的基于隱馬爾科夫的滾動軸承性能衰退模型，利用Forward-Backward算法可以計算相應的后驗概率值P(O|λ)，即P(O,qt=θi|λ)的概率分布。利用后驗概率值得到的各狀態概率值來實現對滾動軸承健康狀態的預測[13]。

3 考慮狀態變遷的滾動軸承健康狀態分類

滾動軸承性能衰退評估首先要對其定性分類，明確軸承從健康狀態到失效過程的各個階段，引入狀態變遷矩陣即在定性分類的基礎上開展定量的評估[14]。滾動軸承全壽命周期將經歷4個階段：磨合期、正常使用期、性能衰退期和快速失效期。

全壽命周期的4個階段對應隱馬爾科夫鏈的狀態數目，對于t時刻隱馬爾科夫鏈所處的狀態為qt={θ1,θ2,θ3,θ4}，對應狀態如表1。

表1 全壽命周期內隱馬爾科夫模型狀態

記t時刻觀測到的測量值為Ot，其中，Ot∈(v1,v2,v3,v4)。

設置初始狀態的分布向量，通常認為滾動軸承的起始狀態是健康狀態，即滾動軸承處于正常使用期，因此，π=(1,0,0,0)。

通過對樣本數據集的統計，利用極大似然估計，估算狀態轉移概率分布矩陣的參數值：

(8)

通過對樣本數據集的統計，利用極大似然估計，估算觀測值概率矩陣的參數值：

(9)

由此完成了基于隱馬爾科夫模型的滾動軸承性能衰退預測的初始化參數設置。隱馬爾科夫模型經參數初始化后，這個參數模型并不是最優化的且符合實際情況的。需要利用多觀測序列的值，通過Baum-Welch算法訓練得出最優化的隱馬爾科夫模型[15]。

定義ξt(i,j)為給定訓練序列O和模型λ時，隱馬爾科夫模型t時刻處于i狀態、t+1時刻處于j狀態的概率，即：

ξt(i,j)=P(qt=i,qt+1=j|O,λ)

(10)

由向前變量與向后變量的定義，式(10)可以變形如式(11)：

(11)

則隱馬爾科夫模型在t時刻處于i狀態的概率為：

(12)

Baum-Welch算法是一種迭代算法，通過不斷的迭代，是參數以極大似然估計趨向于最優值。最終得到合理的狀態轉移概率分布矩陣A和觀測值概率矩陣B，完成隱馬爾科夫模型的建立[16]。

滾動軸承性能衰退的過程與時間序列相關，信號采樣時間是隱馬爾科夫模型建立狀態轉移矩陣考慮的重要參數。不同安裝位置的滾動軸承，在不同時間段所獲數據建立的隱馬爾科夫模型具有相對獨立性。在實際分析過程中可固定采樣時間段，規避時間不同對模型準確性的影響。

4 實例計算

實例計算實驗數據來自美國威斯康辛大學智能維護系統實驗室的滾動軸承試驗，試驗裝置結構如圖5。圖5參照試驗裝置說明繪制。3號滾動軸承振動加速度數據作為訓練樣本，4號滾動軸承振動加速度數據作為測試樣本。預處理采用均方根值，即計算振動加速度數據的RMS值，以證明隱馬爾科夫模型對滾動軸承的性能衰退評估的有效性。

圖5 IMS實驗室滾動軸承全壽命測試裝置示意

ZA-2115雙列滾子軸承依次安裝在軸上，軸由交流電機通過帶傳動驅動，維持在2 000 rmp的回轉速度，同時通過彈簧機構徑向加載6 000磅的壓力。所有的軸承均采用強制潤滑。

根據觀測狀態序列設置，分為4個觀測狀態，健康狀態遷移如圖6。選取滾動軸承振動加速度信號的RMS值的歸一化值作為模型的輸入值，則振動加速度信號RMS歸一化值在1和0之間變化。從1到0的過程表征著軸承由磨合期、正常使用期、性能衰退期到快速失效期4個階段，如表2。

圖6 滾動軸承健康狀態遷移

表2 滾動軸承全壽命周期健康狀態劃分

設置初始化狀態的分布向量為π=(1,0,0,0)，表示初始化狀態下默認軸承是健康的狀態。

根據式(8)對樣本數據統計得到狀態轉移概率分布矩陣A：

初始化狀態分布矩陣A為由滾動軸承歷史壽命數據估算得到的各狀態間轉移的概率。狀態分布矩陣元素值隨著迭代計算不斷更新，最終在數據驅動下實現轉移模型的精確化，即得到觀測值概率矩陣B。

根據式(9)對樣本數據進行統計，得到觀測值概率矩陣B：

綜上，可以得出隱馬爾科夫模型初始化參數，即λ=(π,A,B)。

利用Baum-Welch算法建立隱馬爾科夫模型，根據樣本數據得到最優化的狀態轉移概率分布矩陣A和觀測值概率矩陣B。

優化后的狀態轉移概率矩陣A：

優化后的觀測值概率B：

得到樣本參數各性能衰退狀態的概率值，由此建立基于隱馬爾科夫模型的滾動軸承性能衰退評估方法。選取滾動軸承加速壽命試驗中第1組試驗中的4號軸承振動加速度傳感器數據。數據計算RMS值并進行歸一化處理，按照時間順序抽樣275組數據(分為11組，每組25個數據)組成訓練樣本數據。

從以上11組數據的臨近數據抽取11組作為模型校驗數據，以驗證模型準確性。校驗結果證明，觀測矩陣B對于滾動軸承4的狀態轉移概率準確度高于90%。

依據時間序列隨機選取的11組數據作為測試數據，通過訓練好的隱馬爾科夫模型，依次計算，選取奇數組別數據結果，如圖7。由第1組和第3組測試數據得出，兩組數據基本屬于正常使用期；由第5組和第7組測試數據計算得出，初期性能衰退期和性能衰退期均占有一定的比率，相比之下性能衰退期所占的比重較大；由第9組測試數據計算得出，在性能衰退期的概率較大；由第11組測試數據計算得出，在快速失效期的概率較大。

圖7 滾動軸承性能衰退狀態遷移概率分布

從測試數據的概率分布結果表明：隨著時間變化，正常工作狀態的概率值從1逐漸趨向于0；初期性能衰退狀態和性能衰退狀態的概率先增加后減小，類似于高斯分布，但是初期性能衰退狀態的峰值出現時間比性能衰退狀態峰值出現要早；快速失效狀態的概率值從0逐漸趨向于1。圖7顯示的各個狀態概率值變化情況符合滾動軸承全壽命周期的實際情況。

5 結 論

筆者針對滾動軸承性能衰退過程中狀態變遷具有漸變性的特征，結合隱馬爾科夫模型的隱式和顯式的雙重隨機過程對其過程進行表征，對滾動軸承性能衰退的過程中存在的磨合期、正常使用期、性能衰退期和快速失效期4個階段的概率分布進行計算，得到對應區間的概率分布值和狀態遷移曲線。

實例計算表明：結合隱馬爾科夫模型的滾動軸承性能衰退預測能夠有效判斷其健康狀態的概率分布值并識別性能衰退現象，但是隱馬爾科夫模型的初始狀態矩陣與歷史數據有關，需要先驗數據估算的初始狀態矩陣，存在一定局限性。該時域分析方法與頻域分析方法和時頻域分析方法相比，準確度略低，但是計算速度較快，適用于快速在線性能衰退評估。綜合分析，筆者所提出的滾動軸承的性能衰退識別方法及剩余壽命預測研究具有一定指導意義。