郝 剛，金 濤

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

0 引 言

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

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

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

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

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

1 隱馬爾科夫模型概述

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

圖1 隱馬爾科夫模型

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

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

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

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

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

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

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

(1)

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

經(jīng)初始化：

(2)

遞歸計算：

(3)

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

終結(jié)計算：

(4)

算法示意如圖3。

圖3 前向算法遞歸關(guān)系

該圖反映了在t時刻的狀態(tài)1≥i≥N是通過怎樣的途徑到達t+1時刻的狀態(tài)j。前向算法的思想是通過遞歸方法計算由時刻t的狀態(tài)i向t+1時刻狀態(tài)j轉(zhuǎn)移的所有途徑的概率，這些概率值之和為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)終結(jié)計算：

(7)

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

圖4 后向算法遞歸關(guān)系

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

3 考慮狀態(tài)變遷的滾動軸承健康狀態(tài)分類

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

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

表1 全壽命周期內(nèi)隱馬爾科夫模型狀態(tài)

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

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

通過對樣本數(shù)據(jù)集的統(tǒng)計，利用極大似然估計，估算狀態(tài)轉(zhuǎn)移概率分布矩陣的參數(shù)值：

(8)

通過對樣本數(shù)據(jù)集的統(tǒng)計，利用極大似然估計，估算觀測值概率矩陣的參數(shù)值：

(9)

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

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

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

(10)

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

(11)

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

(12)

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

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

4 實例計算

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

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

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

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

圖6 滾動軸承健康狀態(tài)遷移

表2 滾動軸承全壽命周期健康狀態(tài)劃分

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

根據(jù)式(8)對樣本數(shù)據(jù)統(tǒng)計得到狀態(tài)轉(zhuǎn)移概率分布矩陣A：

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

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

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

利用Baum-Welch算法建立隱馬爾科夫模型，根據(jù)樣本數(shù)據(jù)得到最優(yōu)化的狀態(tài)轉(zhuǎn)移概率分布矩陣A和觀測值概率矩陣B。

優(yōu)化后的狀態(tài)轉(zhuǎn)移概率矩陣A：

優(yōu)化后的觀測值概率B：

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

從以上11組數(shù)據(jù)的臨近數(shù)據(jù)抽取11組作為模型校驗數(shù)據(jù)，以驗證模型準確性。校驗結(jié)果證明，觀測矩陣B對于滾動軸承4的狀態(tài)轉(zhuǎn)移概率準確度高于90%。

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

圖7 滾動軸承性能衰退狀態(tài)遷移概率分布

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

5 結(jié) 論

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

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