基于EMD多尺度威布爾分布與HMM的軸承性能退化評估方法

王 冉, 周雁翔, 胡 雄, 陳 進

(1.上海海事大學(xué) 物流工程學(xué)院，上海 201306；2.上海交通大學(xué) 機械系統(tǒng)與振動國家重點實驗室，上海 200240)

滾動軸承是旋轉(zhuǎn)機械中一種不可或缺的重要部件，對滾動軸承采用基于狀態(tài)的預(yù)測性維護策略是提高旋轉(zhuǎn)設(shè)備安全性與可靠性的重要手段[1-2]。滾動軸承從正常運行到最終失效往往要經(jīng)歷一系列性能退化的過程。與故障診斷技術(shù)不同，性能退化評估更側(cè)重于動態(tài)地揭示設(shè)備在運行過程中健康狀態(tài)退化的規(guī)律，因此更有利于開展預(yù)測性維護，避免機械設(shè)備嚴(yán)重故障的發(fā)生。

滾動軸承性能退化評估的主要步驟包括性能退化特征提取以及評估模型的建立，其核心目標(biāo)在于構(gòu)建一種能夠有效反映軸承性能退化趨勢的性能指標(biāo)，且該指標(biāo)需具有良好的單調(diào)性、敏感性與魯棒性的性能指標(biāo)[3]。常用的性能退化特征包括時域統(tǒng)計特征、頻域特征、時頻域特征等[4]。在評估模型方面，常用的評估模型包括支持向量機、模糊C均值[5-6]、支持向量數(shù)據(jù)描述[7]、隱馬爾可夫模型(hidden Markov model，HMM)[8-10]等。其中，HMM的雙隨機鏈結(jié)構(gòu)能夠較為準(zhǔn)確地描述軸承在退化過程中隱含的衰退狀態(tài)與觀測變量之間的關(guān)系，因此被廣泛應(yīng)用于軸承的性能退化評估中。例如：Jiang等[11]將冗余屬性投影與HMM相結(jié)合建立了性能退化評估模型。Yu[12]建立了一種基于自適應(yīng)HMM的在線學(xué)習(xí)框架用于量化評估軸承的性能退化狀態(tài)。這些方法驗證了HMM在軸承性能退化評估上的有效性與優(yōu)越性。然而，上述方法在魯棒性以及對早期故障的敏感性方面仍然存在著一些不足。在實際工程應(yīng)用中，軸承振動難免受到環(huán)境噪聲、退化過程的隨機性等因素的干擾，使得時、頻域退化特征出現(xiàn)隨機波動，從而影響軸承性能退化評估結(jié)果的穩(wěn)定性和準(zhǔn)確性。尤其在軸承的早期退化階段，軸承故障信息往往淹沒在復(fù)雜的噪聲中[13]，導(dǎo)致難以從中提取出有效的退化特征。

退化特征提取是軸承性能退化評估的關(guān)鍵。與常用的時、頻域等特征提取方法不同，從數(shù)據(jù)的統(tǒng)計分布出發(fā)來提取設(shè)備的退化特征更適用于復(fù)雜的工程數(shù)據(jù)，具有更強的魯棒性。在常見的統(tǒng)計分布模型中，威布爾分布是一種適用于機械零部件磨損累計失效的概率分布模型，更接近實際工程信號的長尾分布特點，在可靠性工程及壽命預(yù)測領(lǐng)域中被廣泛應(yīng)用[14]。Ben Ali等[15]利用威布爾分布對均方根、峭度等軸承時域特征進行擬合，以抑制原始時域特征中的隨機波動和異常干擾，然后把威布爾分布處理后的特征輸入神經(jīng)網(wǎng)絡(luò)實現(xiàn)軸承的剩余壽命預(yù)測。陳昌等[16]將威布爾分布與最小二乘支持向量機結(jié)合對滾動軸承的退化趨勢進行了預(yù)測。侯美慧等[17]將威布爾分布用于岸橋鉸點的退化特征提取中，用實際工況下的岸橋全壽命鉸點振動數(shù)據(jù)驗證了威布爾分布在復(fù)雜工程應(yīng)用中的有效性和可行性。然而，上述方法均采用原始信號的威布爾分布參數(shù)作為衡量設(shè)備性能狀態(tài)變化的指標(biāo)，反映的故障特征信息不夠充分，且原始數(shù)據(jù)包含噪聲干擾，有必要對其進行預(yù)處理以增強故障特征信息。

滾動軸承的振動信號是一種典型的非平穩(wěn)信號，對其進行多尺度分解有助于從不同尺度上挖掘軸承隱藏的性能退化特征。經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)是一種自適應(yīng)的多尺度信號分解方法，能夠把復(fù)雜信號分解為一系列本征模態(tài)分量(intrinsic mode function，IMF)分量之和，非常適用于處理非平穩(wěn)信號，在軸承及齒輪箱的故障診斷中已得到成功應(yīng)用[18-19]。

綜上，為了提高軸承性能退化結(jié)果的準(zhǔn)確性與穩(wěn)定性，針對軸承振動信號非平穩(wěn)性的特點，并且鑒于威布爾分布參數(shù)對性能退化的敏感性以及HMM在性能退化評估中的有效性，本文提出一種基于EMD多尺度數(shù)據(jù)威布爾分布與HMM的軸承性能退化評估方法。該方法通過EMD對軸承振動信號進行多尺度分解，并篩選出包含明顯故障信息的IMF分量，然后從各IMF分量中提取威布爾形狀參數(shù)作為退化特征，最后結(jié)合HMM模型達到準(zhǔn)確評估軸承性能退化狀態(tài)的目的。

1 理論基礎(chǔ)

1.1 經(jīng)驗?zāi)B(tài)分解

EMD是一種適用于非平穩(wěn)、非線性信號的多尺度信號分解方法。它可以將信號自適應(yīng)地分解為若干個本征模態(tài)函數(shù)之和。分解得到的IMF分量須滿足以下兩個條件：

(1) 該函數(shù)局部極值點以及過零點的數(shù)量最多只能相差一個。

(2) 任意時刻，該函數(shù)的上包絡(luò)線與下包絡(luò)線關(guān)于x軸局部對稱。

對任一信號x(t)使用EMD方法進行分解主要包括以下步驟。

步驟1確定該信號的局部極大值以及極小值點，并分別擬合形成上包絡(luò)線及下包絡(luò)線。

步驟2計算上下包絡(luò)線的平均值m1，將均值擬合，得到均值包絡(luò)線。

具體過程為：首先，計算中間信號

h1=x(t)-m1

(1)

然后，判斷h1是否滿足上述IMF分量的兩個條件，若滿足，則h1為原始信號x(t)的第一個IMF分量；若不滿足則以該中間信號作為初始信號，然后重復(fù)步驟1～步驟2得到新的IMF分量。

步驟3當(dāng)函數(shù)hn滿足單調(diào)序列或是常數(shù)時，循環(huán)結(jié)束，記r=hn，r稱為殘余分量。最終，原始信號可表示為多個IMF分量及殘余分量的線性疊加。

1.2 威布爾分布

常用的三參數(shù)威布爾分布的概率密度函數(shù)f(x)及累計密度函數(shù)F(x)可表示為

(2)

(3)

式中:β為形狀參數(shù)；η為尺度參數(shù)；γ為位置參數(shù)。在威布爾分布的3個參數(shù)中，形狀參數(shù)決定了分布曲線的基本形狀，尺度參數(shù)可以將曲線進行放大和縮小，而位置參數(shù)決定了函數(shù)坐標(biāo)系中的起始位置。由此可見，位置參數(shù)與軸承性能退化趨勢無關(guān)，因此主要分析數(shù)據(jù)威布爾分布的形狀參數(shù)與尺度參數(shù)。

本文采用極大似然估計方法對威布爾分布的主要參數(shù)進行估計。給定觀測序列X=x1,…,xN，則在形狀參數(shù)β與尺度參數(shù)η下觀測序列的極大似然函數(shù)為

(4)

式中，N為數(shù)據(jù)長度。根據(jù)極大似然原理及牛頓迭代法可求解得到形狀參數(shù)β及尺度參數(shù)η。

1.3 隱馬爾可夫模型

隱馬爾可夫模型是一種雙重隨機過程，不僅狀態(tài)轉(zhuǎn)移過程是隨機的，且每一狀態(tài)的觀測值也是隨機的。其中，描述狀態(tài)轉(zhuǎn)移的馬爾科夫過程不能被直接觀測，但可利用觀測過程進行估計。

一個典型的離散HMM可用以下參數(shù)描述

(1)N：模型中隱含狀態(tài)的數(shù)目。假定狀態(tài)集合為S={S1,…,SN}，記t時刻模型所處狀態(tài)為qt，則qt∈{S1,…,SN}。

(2)M：各個狀態(tài)對應(yīng)的觀測值數(shù)目。記觀測值集合為V={v1,…,vM}，若ot為t時刻的觀測值，則ot∈V。

(3) π：初始狀態(tài)的概率向量π={π1,…,πN}，其中πi=P(q1=Si)，1≤i≤N。

(4)A：狀態(tài)轉(zhuǎn)移概率矩陣A=(aij)N×N，其中aij=P(qt+1=Sj|qt=Si),1≤i,j≤N。

(5)B：觀測值概率矩陣B=(bjk)N×N，其中bjk=P(ot=vk|qt=Sj),1≤j≤N,1≤k≤M。

典型的HMM可由以上5個參數(shù)表示，常被簡記為λ=(π,A,B)。在HMM的訓(xùn)練過程中，常采用期望最大法對模型參數(shù)進行估計。然后，在模型參數(shù)已知的條件下，采用前向-后向算法對觀測序列O={o1,…,oT}的概率P(O|λ)進行計算。在前向算法中，定義前向變量為1～t時刻觀測序列與t時刻模型處于狀態(tài)Si的聯(lián)合概率，即

αt(i)=P(o1,…,ot,qt=Si|λ)

(5)

記前向變量的初始值α1(i)=πibi(o1)，則前向變量的迭代公式為

(6)

迭代結(jié)束后，觀測序列的概率可表示為

(7)

為了防止數(shù)據(jù)下溢，通常對式(7)的概率取對數(shù)，采用對數(shù)似然概率log(P(O|λ))作為給定模型參數(shù)情況下觀測值出現(xiàn)的概率，簡寫為LLP。

需要說明的是：上述HMM的觀測值序列及相應(yīng)的觀測概率矩陣B是離散的，然而軸承的振動信號是連續(xù)的，需要用連續(xù)的概率密度函數(shù)來描述觀測值的概率分布。由于混合高斯模型可以無限逼近任意分布，因此常采用混合高斯模型來擬合任意狀態(tài)Sj下觀測值ot的概率密度函數(shù)，即

(8)

式中:Mj為狀態(tài)Sj的高斯分量數(shù)目;wjm為第m個高斯模型的權(quán)重;μjm與Vjm分別為第m個高斯分布的均值向量與協(xié)方差矩陣。

2 基于EMD多尺度威布爾分布與HMM的軸承性能退化評估方法

本文提出了基于EMD多尺度威布爾分布與HMM的軸承性能退化評估方法。該方法首先采用經(jīng)驗?zāi)B(tài)分解將軸承振動數(shù)據(jù)分解到不同尺度的IMF分量中，并以峭度為標(biāo)準(zhǔn)選取IMF分量；然后提取每個IMF分量的威布爾分布形狀參數(shù)構(gòu)建退化特征集，將其輸入HMM模型中進行訓(xùn)練，最終得到性能退化評估曲線。該方法分為離線訓(xùn)練與在線評估兩個階段，具體流程如圖1所示，主要步驟包括以下。

圖1 軸承性能退化評估流程圖

步驟1利用EMD對采集的軸承振動數(shù)據(jù)進行多尺度分解。

步驟2分別計算各個IMF分量的峭度值，篩選得到故障信息明顯的多個IMF分量。

步驟3對各個IMF分量分別對一個滑動時間窗口內(nèi)的數(shù)據(jù)進行威布爾分布擬合，提取威布爾分布的形狀參數(shù)構(gòu)建多維性能退化特征集。

步驟4將軸承正常狀態(tài)下的數(shù)據(jù)通過步驟1～步驟3得到的退化特征集作為訓(xùn)練數(shù)據(jù)輸入HMM模型，得到性能退化評估模型λ。

步驟5在線評估階段，將待測數(shù)據(jù)的退化特征輸入評估模型λ，利用前向算法中的式(7)計算得到待測數(shù)據(jù)滿足正常模型λ的對數(shù)似然概率LLP。若LLP的值較大，說明待測數(shù)據(jù)處于正常狀態(tài)的概率較高；反之，若LLP的值變小，則說明軸承當(dāng)前待測狀態(tài)開始偏離正常狀態(tài)。因此，可根據(jù)LLP的變化趨勢及時發(fā)現(xiàn)軸承的早期故障，評估其退化狀態(tài)。

3 試驗驗證

本文使用的滾動軸承全壽命周期振動信號來自杭州軸承試驗中心的ABLT-1A軸承壽命強化試驗臺。試驗裝置如圖2(a)所示，主要包括臺架、加速度傳感器、測試軸承、采集系統(tǒng)等，可同時進行4個軸承的加速疲勞試驗。4個測試軸承及3個加速度傳感器(通道1～3)的安裝位置如圖2(b)所示。

(a) 試驗測試裝置

試驗軸承采用型號為6307的單列深溝球軸承采樣頻率為25.6 kHz，每一分鐘采集一次，采樣時間為0.8 s。本文采用最先失效的B12軸承對所提出的方法進行驗證。試驗從開始采集直至軸承完全失效共經(jīng)歷2 469 min，得到2 469組振動數(shù)據(jù)。

軸承全壽命信號的時域波形及有效值曲線如圖3所示。從圖3可以區(qū)分出軸承的正常階段及失效階段，但對軸承早期故障的識別較為困難。同時，如圖3(b)所示，有效值指標(biāo)存在波動和異常值，不宜直接作為軸承的性能指標(biāo)。

(a) 軸承全壽命數(shù)據(jù)

3.1 EMD分解與威布爾分布校驗

首先，對軸承振動信號進行EMD分解，然后根據(jù)峭度值的大小篩選出主要的IMF分量，得到高于各分量平均峭度值的6個IMF分量，即IMF2、IMF3、IMF10、IMF11、IMF13、IMF17。對各IMF分量進行威布爾分布校驗，檢驗結(jié)果表明各分量數(shù)據(jù)均近似滿足威布爾分布。其中，IMF3分量整體數(shù)據(jù)的威布爾分布校驗結(jié)果如圖4所示。

圖4 IMF3威布爾分布校驗

3.2 特征提取

以IMF3為例，取窗口大小為10 min，步長為1 min，每分鐘含20 480個數(shù)據(jù)點，對該分量分別提取威布爾尺度參數(shù)、形狀參數(shù)等特征指標(biāo)，結(jié)果如圖5所示。可見，尺度參數(shù)不隨故障發(fā)展呈單調(diào)變化趨勢，且無法有效反映軸承早期故障，不宜作為評估指標(biāo)。相比之下，形狀參數(shù)與軸承退化趨勢較為一致，單調(diào)性較好，且在軸承的退化初期有所變化，因此本文選取形狀參數(shù)作為軸承的退化特征。

(a) IMF3尺度參數(shù)

對6個IMF分量分別提取威布爾形狀參數(shù)，組成多維退化特征集W=[W1,W2,W3,W4,W5,W6]。

3.3 HMM模型訓(xùn)練與性能退化評估

在HMM訓(xùn)練階段，選擇軸承正常狀態(tài)下的前300組數(shù)據(jù)作為訓(xùn)練集，將其退化特征集作為HMM的輸入，建立HMM模型進行訓(xùn)練。試驗中，混合高斯模型數(shù)M取2，隱含狀態(tài)數(shù)N取5。

在評估階段，將軸承全壽命數(shù)據(jù)作為測試集進行性能退化評估，得到的性能退化曲線如圖6所示。可見，第1 291 min前HMM輸出概率值較高且較為平穩(wěn)，表明軸承運行狀況良好；在1 291 min后，性能指標(biāo)開始下降，表明軸承早期故障的出現(xiàn)；在2 330 min后，曲線出現(xiàn)劇烈下降，表明軸承性能大幅退化，進入故障失效階段。

圖6 基于EMD多尺度數(shù)據(jù)威布爾分布與HMM的性能退化評估結(jié)果

為了驗證所提出方法的準(zhǔn)確性，對兩個變化點即第1 291 min和2 330 min的數(shù)據(jù)進行包絡(luò)譜分析。本試驗中，軸承轉(zhuǎn)頻fr=50 Hz，內(nèi)圈故障特征頻率fi=246 Hz。1 291 min與2 330 min的時域波形及包絡(luò)譜如圖7所示。如圖7(a)所示，1 291 min的時域波形中存在微弱的沖擊成分。同時，如圖7(c)所示，從其包絡(luò)譜中可以觀測到轉(zhuǎn)頻fr及內(nèi)圈故障頻率fi成分，說明此時軸承出現(xiàn)早期故障，發(fā)生輕微點蝕。在第2 330 min，圖7(b)的時域波形振幅明顯增大，沖擊信號增多。同時，如圖7(d)所示，包絡(luò)譜中故障特征頻率分量fi較為明顯，且受轉(zhuǎn)頻調(diào)制影響，在fi附近出現(xiàn)明顯的邊頻帶，表明軸承出現(xiàn)內(nèi)圈嚴(yán)重故障。上述分析結(jié)果驗證了評估結(jié)果的正確性與有效性。

(a) 1 291 min時域波形

為了進一步驗證方法的優(yōu)越性，將所提方法與以下6種方法進行對比分析。采用常用的11個時域特征指標(biāo)、時域特征與主成分分析(principal components analysis，PCA)降維、時域特征與局部保持投影(locality preserving projections，LPP)、EMD后各分量的能量指標(biāo)、原始單一尺度信號的威布爾形狀參數(shù)、EMD多尺度數(shù)據(jù)的威布爾尺度參數(shù)等退化特征提取方法分別與HMM相結(jié)合，得到的性能退化曲線如圖8所示。可見，圖8(a)～圖8(d)的性能指標(biāo)曲線在軸承正常階段均存在異常波動，容易導(dǎo)致誤判，且無法識別出軸承的早期故障。圖8(e)的性能退化曲線在正常階段雖然比較平穩(wěn)，但是仍然難以識別軸承的早期故障。圖8(f)所采用的方法與本文最為相近，從圖8可以看出該方法雖然對軸承的早期故障有一定反映，但是整體趨勢不明顯，波動較大，說明本文采用的威布爾分布的形狀參數(shù)比圖8(f)采用的尺度參數(shù)更適合作為性能退化特征。根據(jù)圖8劃分的具體退化階段評估結(jié)果如表1所示。可見，僅有本方法能及時識別出軸承的早期故障，并能準(zhǔn)確區(qū)分出軸承的早期故障與嚴(yán)重失效階段。

(a) 時域特征+HMM

4 結(jié) 論

本文針對現(xiàn)有軸承性能退化評估方法在魯棒性和早期故障敏感性上的不足，提出一種基于EMD多尺度威布爾分布與HMM的軸承性能退化評估方法，并結(jié)合滾動軸承全壽命試驗數(shù)據(jù)對該方法進行對比分析與驗證，結(jié)論如下：

(1) 采用EMD將軸承振動數(shù)據(jù)中的退化趨勢分解到不同尺度的分量中，與單一尺度相比能夠更有效地挖掘軸承的故障特征信息。

(2) 采用威布爾分布的形狀參數(shù)作為性能退化指標(biāo)，與傳統(tǒng)時域特征作為退化指標(biāo)相比能夠?qū)S承早期故障做出敏感反映。

(3) 與其他常用軸承性能退化評估方法的對比發(fā)現(xiàn)，本文所提出的評估方法能夠準(zhǔn)確檢測到早期故障的發(fā)生，評估結(jié)果具有較高的準(zhǔn)確性及魯棒性。