王 冉, 周雁翔, 胡 雄, 陳 進(jìn)

(1.上海海事大學(xué) 物流工程學(xué)院，上海 201306；2.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240)

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

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

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

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

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

1 理論基礎(chǔ)

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

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

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

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

對(duì)任一信號(hào)x(t)使用EMD方法進(jìn)行分解主要包括以下步驟。

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

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

具體過(guò)程為：首先，計(jì)算中間信號(hào)

h1=x(t)-m1

(1)

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

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

1.2 威布爾分布

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

(2)

(3)

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

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

(4)

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

1.3 隱馬爾可夫模型

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

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

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

(2)M：各個(gè)狀態(tài)對(duì)應(yīng)的觀測(cè)值數(shù)目。記觀測(cè)值集合為V={v1,…,vM}，若ot為t時(shí)刻的觀測(cè)值，則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：觀測(cè)值概率矩陣B=(bjk)N×N，其中bjk=P(ot=vk|qt=Sj),1≤j≤N,1≤k≤M。

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

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

(5)

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

(6)

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

(7)

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

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

(8)

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

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

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

圖1 軸承性能退化評(píng)估流程圖

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

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

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

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

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

3 試驗(yàn)驗(yàn)證

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

(a) 試驗(yàn)測(cè)試裝置

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

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

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

3.1 EMD分解與威布爾分布校驗(yàn)

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

圖4 IMF3威布爾分布校驗(yàn)

3.2 特征提取

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

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

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

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

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

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

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

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

(a) 1 291 min時(shí)域波形

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

(a) 時(shí)域特征+HMM

4 結(jié) 論

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

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

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

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