基于CEEMD散布熵和Hjorth參數(shù)的混合特征滾動軸承故障診斷研究*

夏理健,劉小平,王 新,田 笑,張立杰

(1.燕山大學(xué) 河北省重型機(jī)械流體動力傳輸與控制重點(diǎn)實驗室,河北 秦皇島 066004;2.北京航空航天大學(xué) 自動化科學(xué)與電氣工程學(xué)院,北京 100191;3.燕山大學(xué) 先進(jìn)鍛壓成形技術(shù)與科學(xué)教育部重點(diǎn)實驗室,河北 秦皇島 066004)

0 引 言

由于滾動軸承經(jīng)常在重載、沖擊和變載荷等復(fù)雜工況條件下工作,其故障的發(fā)生率極高。滾動軸承發(fā)生故障時,會嚴(yán)重影響機(jī)械設(shè)備的正常運(yùn)行,甚至?xí)l(fā)安全事故[1]。因此,對滾動軸承進(jìn)行狀態(tài)檢測和故障診斷,對機(jī)械設(shè)備的安全運(yùn)行具有重要意義。

當(dāng)滾動軸承發(fā)生局部故障時,其振動信號通常呈現(xiàn)出非線性和非平穩(wěn)性的特征。早期的脈沖信號非常微弱,易受到強(qiáng)環(huán)境噪聲影響,導(dǎo)致故障特征難以提取。因此,有必要通過信號分解的方式來提取其故障特征。

采用經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)及其改進(jìn)算法獲取振動信號特征的方法是當(dāng)前的熱點(diǎn)研究方向之一。劉永強(qiáng)等人[2]利用集合經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)和自相關(guān)函數(shù)峰態(tài)系數(shù),對軸承進(jìn)行了故障診斷,但每次EEMD分解需要加入不同幅值的白噪聲,使各個固有模態(tài)分量(IMF)中有殘留噪聲,影響了序列中的有效信息提取,導(dǎo)致其分解精度不高。TORRES M E等人[3]提出了一種完備總體經(jīng)驗?zāi)B(tài)分解(CEEMD)方法,采用該方法可以減小由白噪聲導(dǎo)致的重構(gòu)信號誤差,更好地抑制模態(tài)混疊現(xiàn)象,提高分解精度。

熵是分析信號動態(tài)變化的有力工具,可以反映信號的混亂程度。信息熵能夠有效地檢測故障振動信號的動力學(xué)突變,因此,利用熵能夠反映信號復(fù)雜程度的特性,可以達(dá)到提取特征信息的目的。其中,排列熵(permutation entropy,PE)應(yīng)用較為廣泛。

陳祥龍等人[4]提出了一種采用改進(jìn)排列熵,來表示滾動軸承振動信號特征的方法。石志煒等人[5]先用CEEMD分解滾動軸承的振動信號,然后求解IMF分量的排列熵,從而得到了其特征向量。雖然PE計算簡單,但其沒有考慮幅值之間的大小關(guān)系。2016年，ROSTAGHI M等人[6]610-611提出了一種新的信號不規(guī)則程度指標(biāo)—散布熵(DE)。DE算法的速度很快,且考慮了幅值間的關(guān)系,具有更好的抗噪能力。李從志等人[7]采用了一種將EMD與DE相結(jié)合的方法,實現(xiàn)了對滾動軸承的故障診斷。

由于單一方面的特征只包含部分信息,利用多個物理特征進(jìn)行混合特征提取可以全面體現(xiàn)故障特征信息,有利于提取更廣泛的故障特征信息。

Hjorth參數(shù)最早由HJORTH B[8]提出,目前已被廣泛應(yīng)用于腦電信號的特征提取和分析處理。Hjorth參數(shù)是一種可以同時描述振動信號在時域、頻域中瞬時特征的統(tǒng)計函數(shù)。CAESARENDRA W等人[9]在對低速回轉(zhuǎn)軸承狀態(tài)進(jìn)行監(jiān)測的研究中,應(yīng)用3個Hjorth參數(shù),但在軸承狀態(tài)監(jiān)測中的應(yīng)用效果并不理想。GROVER C等人[10]將EMD與Hjorth參數(shù)結(jié)合起來,以計算相關(guān)性最高的IMF的Hjorth參數(shù)作為特征向量,來對滾動軸承進(jìn)行故障診斷,取得了不錯的效果。

但是,只選用相關(guān)性最高的IMF來進(jìn)行計算,必然會丟掉部分故障信號的信息。因此,該方法還有待于進(jìn)一步改進(jìn)。

基于上述問題,筆者提出一種基于CEEMD的散布熵與Hjorth參數(shù)的混合特征故障診斷方法。首先,對滾動軸承信號進(jìn)行CEEMD分解,得到若干個IMF分量,并計算各個分量與原始信號的相關(guān)性,選取相關(guān)性較高的前幾個IMF分量;然后,求所選IMF的散布熵和Hjorth參數(shù),形成散布熵特征向量和Hjorth參數(shù)矩陣,對Hjorth參數(shù)矩陣進(jìn)行奇異值分解(SVD),利用奇異值向量與散布熵特征向量形成混合特征向量來代表滾動軸承振動信號;最后,利用基于粒子群優(yōu)化算法(PSO)優(yōu)化的最小二乘支持向量機(jī)(LSSVM),對滾動軸承不同故障特征向量進(jìn)行訓(xùn)練和識別,利用滾動軸承不同故障試驗數(shù)據(jù)對所提方法進(jìn)行驗證。

1 基本理論和方法

1.1 CEEMD算法理論

由于EMD算法是從原信號中提取若干個固有模態(tài)分量IMF,每一階IMF都反映原始信號的動態(tài)特性。但一些有異常干擾的非線性信號會產(chǎn)生模態(tài)混疊,導(dǎo)致多個模擬主導(dǎo)頻率分量同時出現(xiàn)在一個模態(tài)分量中。

EEMD將白噪聲加入到原始信號當(dāng)中,利用白噪聲頻譜的均勻分布,使不同時間尺度的信號自動分布到合適的參考尺度上,可有效抑制模態(tài)混疊;但添加的白噪聲會使IMF分量在重構(gòu)時產(chǎn)生誤差。

為了解決白噪聲的干擾,且在重構(gòu)信號時產(chǎn)生誤差的問題,通常采用CEEMD算法,通過給原始信號添加符號相反的白噪聲,分別對兩組信號進(jìn)行EMD分解,使重構(gòu)誤差明顯減小,添加的白噪聲得到中和。

CEEMD分解的具體過程如下:

(1)初始化添加輔助白噪聲的次數(shù)n和白噪聲幅值k,令i=1;

(2)將第i次添加白噪聲的信號進(jìn)行EMD分解。在原信號中以正負(fù)成對的形式加入輔助白噪聲,從而得到信號:

(1)

式中:x(t)—原始信號;ni(t)—第i次添加的輔助噪聲;Pi(t)—加入正噪聲后的信號;Ni(t)—加入負(fù)噪聲后的信號;n—信號個數(shù)。

對加入白噪聲后的信號Pi(t)和Ni(t)進(jìn)行EMD分解,使每個信號得到q個IMF分量,即:

(2)

其中：j=1,2,…,q。

(3)如果i

(3)

式中:cj(t)—信號經(jīng)CEEMD分解后得到的第j個IMF分量;r(t)—最終殘余分量。

該方法保證了信號分解的完備性,可以較好地解決模態(tài)混疊效應(yīng)問題。需要注意的是,添加高斯白噪聲的次數(shù)n一般取100～300[11],此時噪聲殘留所引起的誤差非常小。

白噪聲幅值k的值通常取為原始信號標(biāo)準(zhǔn)差的0.2～0.5倍,取值也可隨噪聲的強(qiáng)度而適當(dāng)調(diào)整增大[12]。此處所使用的n取100,k取0.2σx(σx為原始信號的標(biāo)準(zhǔn)差)。

1.2 散布熵算法

散布熵算法是一種用來度量時間序列復(fù)雜性和不規(guī)則程度的新算法。散步熵對同步頻率、振幅值和時間序列帶寬的變化很敏感,而且受突變信號影響較小,具有更好的抗噪能力,因此,它比排列熵的計算效率更高。

對于時間序列x={xj,j=1,2,…,N},散布熵的計算步驟如下:

(1)利用正態(tài)累計分布函數(shù),將x映射為y={yj,j=1,2,…,N},yj∈(0,1),其正態(tài)分布函數(shù)為:

(4)

式中:μ—時間序列x的期望;σ—時間序列x的標(biāo)準(zhǔn)差。

(2)通過線性變換把y映射到[1,2,…,c]范圍,即:

(5)

(6)

式中:m—嵌入維數(shù);d—時延。

(5)計算每種散布模式πv0v1…vm-1的概率p(πv0v1…vm-1):

(7)

(6)根據(jù)香農(nóng)熵的定義,對于嵌入維數(shù)為m、時間延遲為d以及類別數(shù)為c的原始時間序列x,其歸一化的散布熵可表示為:

(8)

散布熵與信號的不規(guī)則程度相關(guān),DE值越大,信號不規(guī)則程度越高,反之越低。

在參數(shù)選擇方面,文獻(xiàn)[6]611建議m通常取2或3,c取[4,8]間的整數(shù),時間序列長度大于2 000。時延d的值大于1可能會造成信息的丟失,故筆者取m=3,c=6,d=1,每段時間序列長度為2 048。

1.3 Hjorth參數(shù)

Hjorth參數(shù)是為計算時變信號的平均功率、均方根頻率以及均方根頻率展開等信號特征時,提供的一種方法。使用Hjorth參數(shù)進(jìn)行計算時,由于涉及方差,使得其計算成本非常低。

3個Hjorth參數(shù)分別稱為:活動性(Activity)、移動性(Mobility)和復(fù)雜性(Complexity)。

活動性定義為信號的方差,表示信號的幅度特性,其公式如下:

Activity=σ2

(9)

式中:σ—信號的標(biāo)準(zhǔn)差。

移動性定義為信號的一階差分信號的方差和信號自身方差之比的均方根。具體公式如下:

(10)

式中:σ′—原始信號的一階差分信號的標(biāo)準(zhǔn)差。

復(fù)雜性定義為振動信號一階導(dǎo)數(shù)的遷移率與振動信號的遷移率之比。它給出了信號帶寬的估計值,表明振動信號與純正弦波的相似程度。其公式如下:

(11)

式中:σ″—原始信號二階差分的標(biāo)準(zhǔn)差。

2 基于CEEMD的混合特征提取

基于CEEMD散布熵和Hjorth參數(shù)的混合特征提取步驟為:

(1)對采集到的信號用CEEMD進(jìn)行分解,得到若干個IMF分量。

(2)為了盡可能提取有用的特征信息,并去除虛假信息,需要對虛假的IMF分量進(jìn)行剔除;通過計算所有IMF分量與原信號的相關(guān)性系數(shù)Cr,篩選出相關(guān)系數(shù)較大的前m個IMF分量,來代表原信號中的有效信息。

相關(guān)系數(shù)Cr的計算公式如下:

(12)

(3)計算篩選出的m個IMF分量的散布熵,并構(gòu)成散布熵特征向量E:

(13)

(4)計算篩選出的m個IMF分量的Hjorth參數(shù),并構(gòu)成Hjorth參數(shù)矩陣H:

(14)

式中:Ai,Mi,Ci(i=1,2,…,m)—第i個IMF分量的活動性、移動性和復(fù)雜性。

(5)為了把H矩陣中的信息和散布熵特征向量E中的信息拼接成混合特征向量,必須把H矩陣中的信息轉(zhuǎn)化成向量的形式。奇異值往往對應(yīng)著矩陣中隱含的重要信息,且重要性和奇異值大小正相關(guān)。奇異值由于具有穩(wěn)定性及比例不變性等良好性質(zhì),常被用來作為矩陣固有特征[13]。

故此處對H矩陣進(jìn)行奇異值分解,提取分解后的奇異值向量Sv。SVD的具體公式如下:

(15)

式中:S=diag(θ1,θ2,…,θq)—奇異值組成的對角陣,θ1≥θ2≥…≥θq,θi(i=1,2,…,q)—H矩陣的奇異值;U,V—左右正交矩陣;ui—矩陣U的第i列;vi—矩陣V的第i列。

其中：q=min(m,3)。

把S對角陣轉(zhuǎn)化成能夠代表H矩陣的奇異值特征向量Sv:

(16)

(6)把散布熵特征向量E和奇異值特征向量Sv拼接起來,構(gòu)成能夠代表原信號信息特征的混合特征向量F:

F=[ESv]=[e1e2…emθ1θ2…θq]

(17)

該方法從多個方面提取信號的特征,得到的信息更加全面,能夠更好地突出不同信號之間的不同特征信息。

3 基于PSO-LSSVM的軸承故障診斷

3.1 最小二乘支持向量機(jī)模型

LSSVM是在原有支持向量機(jī)基礎(chǔ)上的一種改進(jìn)算法,在解決小樣本故障數(shù)據(jù)、非線性問題時具有很大優(yōu)勢。LSSVM的具體算法如下:

設(shè)給定的訓(xùn)練樣本{(x1,y1),(x2,y2),…,(xn,yn)},其中,xi為輸入樣本,yi為輸出標(biāo)簽。

線性決策函數(shù)構(gòu)造如下:

f(xi)=ωTφ(xi)+β

(18)

式中:ω—權(quán)重向量;β—偏差量;φ(·)—核函數(shù)。

按照結(jié)構(gòu)最小化原理,可把LSSVM的優(yōu)化問題轉(zhuǎn)化為以下表達(dá)式:

(19)

式中:ei—樣本的擬合誤差;γ—正則化參數(shù)。

為了解決上述問題,筆者構(gòu)建拉格朗日函數(shù):

(20)

式中:λi—拉格朗日算子。

其最優(yōu)化條件為:

(21)

消去參數(shù)ei和ω,經(jīng)整理可得線性方程組:

(22)

式中:矩陣Ω=yiyjφT(xi)φ(xj)=yiyjK(xi,xj);K(·)—高斯核函數(shù);yT—標(biāo)簽矩陣;I—單位矩陣。

最后,LSSVM的決策函數(shù)為:

(23)

LSSVM將最優(yōu)化問題轉(zhuǎn)換成線性方程組問題,提高了其運(yùn)算速度。

3.2 整體流程

基于CEEMD的散布熵和Hjorth參數(shù)的混合特征故障診斷方法主要包括3個階段:故障特征提取、故障分類模型構(gòu)建和故障診斷。

具體3個階段如下:

(1)利用CEEMD分解原信號,篩選出相關(guān)性高的IMF,計算其散布熵特征和Hjorth參數(shù)特征,形成混合特征向量,構(gòu)建初始特征集;

(2)將初始特征集劃分為訓(xùn)練集和測試集兩部分,利用訓(xùn)練集樣本訓(xùn)練最小二乘支持向量機(jī)分類器,構(gòu)建PSO-LSSVM模型;

(3)利用PSO-LSSVM模型對測試集樣本進(jìn)行故障診斷,識別故障類型。

滾動軸承故障診斷方法的流程圖如圖1所示。

圖1 滾動軸承故障診斷方法流程

4 實驗與結(jié)果分析

此處實驗采用的是美國凱斯西儲大學(xué)電氣工程實驗室的滾動軸承數(shù)據(jù)[14]。該數(shù)據(jù)是使用安裝在驅(qū)動端型號為SKF6205-2RS的深溝球軸承上采集得到的。

實驗軸承數(shù)據(jù)信息如下:采樣頻率12 kHz,電機(jī)轉(zhuǎn)速1 730 r/min;故障軸承采用電火花加工出單點(diǎn)損傷,滾動體、內(nèi)圈和外圈的故障直徑各取3種,外圈滾道的故障在6點(diǎn)鐘位置。

正常信號和不同故障程度的故障信號共10種,每種數(shù)據(jù)取50個樣本,共500個樣本。其中,訓(xùn)練集的樣本數(shù)為400,測試集的樣本數(shù)為100,樣本長度為2 048。

數(shù)據(jù)集的故障類別描述如表1所示。

表1 滾動軸承10種故障類別的數(shù)據(jù)集狀態(tài)

滾動軸承的振動信號中,每種狀態(tài)典型樣本的振動圖像如圖2所示(其中,橫軸代表采樣點(diǎn)1～2 048,縱軸表示幅值)。

圖2 滾動軸承的振動信號

首先,在10種狀態(tài)下,各樣本信號先經(jīng)CEEMD分解得到各個IMF分量,然后計算出與各自原始信號相關(guān)系數(shù)的平均值,如圖3所示。

圖3 10種狀態(tài)各個IMF分量與原信號的相關(guān)系數(shù)

圖3中,幾乎所有前3個IMF分量與原始信號的相關(guān)系數(shù)都在0.1以上,說明前3個IMF分量與原始信號相關(guān)性較強(qiáng),涵蓋了信號中的主要故障信息。

綜上所述,筆者選取前3個IMF分量,分別計算其散布熵和Hjorth參數(shù),并對Hjorth參數(shù)矩陣奇異值進(jìn)行分解,把得到的三維奇異值向量和三維散布熵特征向量拼接成混合特征向量,共可得到10個特征集,從而完成對滾動軸承的特征提取。

為了說明DE算法的優(yōu)越性,筆者將DE與PE進(jìn)行對比,即采用兩種算法,分別計算信號的前3個IMF分量的熵值;然后以IMF1的熵值為橫坐標(biāo)值,以IMF2的熵值為縱坐標(biāo)值,以IMF3的熵值為豎坐標(biāo)值,繪制成特征散點(diǎn)分布圖。

其中,3D散布熵特征散點(diǎn)分布圖如圖4所示。

圖4 3D散布熵特征散點(diǎn)分布圖

3D排列熵特征散點(diǎn)分布圖如圖5所示。

圖5 3D排列熵特征散點(diǎn)分布圖

從圖4和圖5中的坐標(biāo)跨度可以看出:不同故障類別DE特征之間的類別距離相對較大,在不同故障之間有更好的區(qū)分度;而不同故障類別的PE特征之間聚集現(xiàn)象更明顯,只有正常信號與所有故障之間的區(qū)分較大。

為了驗證Hjorth參數(shù)矩陣奇異值分解形成的奇異值特征向量與故障之間的相關(guān)性,筆者計算了每個奇異值特征與故障標(biāo)簽之間的最大信息系數(shù)。

最大信息系數(shù)[15]是用來捕捉每個特征與標(biāo)簽之間的任意關(guān)系(包括線性和非線性關(guān)系)的過濾方法,它返回每個特征與目標(biāo)之間的互信息量的估計量,這個估計量的取值在[0,1]之間(其中,0表示2個變量獨(dú)立,1表示2個變量完全相關(guān))。3個奇異值特征與故障標(biāo)簽之間的最大信息系數(shù)分別為0.919 6、0.741 8、0.958 0,均大于0.5,說明Hjorth參數(shù)和故障之間有很好的相關(guān)性。

為進(jìn)一步驗證3個奇異值與故障之間的相關(guān)性,筆者分別使用1～3個奇異值特征進(jìn)行故障診斷,其識別準(zhǔn)確率如表2所示。

表2 基于不同數(shù)量奇異值特征的準(zhǔn)確率

從表2中可以看出:使用奇異值在1～3個特征數(shù)量,其識別準(zhǔn)確率單調(diào)遞增,說明每個奇異值特征都與故障信息相關(guān)。

在不同故障狀態(tài)下,滾動軸承的三維奇異值特征散點(diǎn)圖,如圖6所示。

圖6 3D奇異值特征散點(diǎn)分布圖

由圖6可知:不同故障狀態(tài)下,特征的分布有明顯的差異,而且坐標(biāo)跨度比三維PE特征散點(diǎn)圖大,故可以說明,將奇異值序列作為特征向量，可以較為明顯地區(qū)分開不同的故障狀態(tài)。

筆者將每種故障混合特征樣本數(shù)據(jù)隨機(jī)打亂,分成400個訓(xùn)練集和100個測試集,用訓(xùn)練集訓(xùn)練PSO-LSSVM模型,用測試集來進(jìn)行驗證。為了使診斷結(jié)果可視化,筆者采用混淆矩陣來表現(xiàn)每一類故障的分類情況,并利用精確度P和召回率R[16]對其進(jìn)行評估。

在故障測試集上,滾動軸承的分類結(jié)果如圖7所示(分類錯誤的個數(shù)為0,即所有故障狀態(tài)的精確度和召回率都是100%)。

圖7 軸承測試樣本分類歸一化混淆矩陣

從以上分析中可知,基于散布熵和Hjorth參數(shù)的混合特征能夠較好地反映不同的故障狀態(tài)。

為了對基于散布熵和Hjorth參數(shù)(表示的是Hjorth參數(shù)矩陣分解的奇異值向量)的混合特征方法在故障分類方面的優(yōu)勢進(jìn)行評估,筆者將該方法與只用單一排列熵、散布熵和Hjorth參數(shù)矩陣奇異值特征方法進(jìn)行比較,并將其與基于排列熵和Hjorth參數(shù)的混合特征進(jìn)行比較,結(jié)果如表3所示。

表3 不同方法的診斷結(jié)果

由表3可以看出:

基于CEEMD散布熵和Hjorth參數(shù)的混合特征方法的診斷準(zhǔn)確率最高,達(dá)到了100%,利用排列熵和Hjorth參數(shù)矩陣奇異值的混合特征方法的準(zhǔn)確率略低;同時,利用單一特征的診斷準(zhǔn)確率均低于利用混合特征的準(zhǔn)確率,證明了CEEMD相對于EEMD算法的優(yōu)越性;筆者提出的Hjorth參數(shù)矩陣特征提取算法,在不同的信號分解算法中都保持了較高的準(zhǔn)確率,受到的影響也較小。

以上分析結(jié)果充分驗證了筆者采用的混合特征在軸承故障診斷中的有效性和準(zhǔn)確性。

5 結(jié)束語

筆者提出了一種基于CEEMD散布熵和Hjorth參數(shù)的混合特征滾動軸承故障診斷方法;首先提取出了故障軸承振動信號多個物理方面的特征信息,然后構(gòu)建起了PSO-LSSVM故障診斷模型,實現(xiàn)了對滾動軸承的故障診斷,最后通過實驗驗證了該混合特征提取方法在對滾動軸承多狀態(tài)故障進(jìn)行診斷時的可行性和優(yōu)越性。

研究結(jié)果表明:

(1)利用CEEMD算法和相關(guān)系數(shù)法提取敏感IMF,可以提取出與故障信息相關(guān)的IMF,解決EMD分解中的模態(tài)混疊問題;

(2)把振動信號的散布熵和Hjorth參數(shù)混合特征用于故障診斷,能夠從多個維度反映軸承故障的特征,解決單一物理特征診斷精度低的問題;

(3)對比分析結(jié)果表明,該方法的診斷結(jié)果優(yōu)于基于單一特征的故障診斷方法。

在今后的研究中,筆者將在考慮更多影響模型的因素的基礎(chǔ)上,針對支持向量機(jī)進(jìn)行更深入的理論研究。同時,由于粒子群優(yōu)化算法的初始參數(shù)設(shè)置基于個人經(jīng)驗,不同初始參數(shù)是否會影響故障的識別效果,還有待于筆者做進(jìn)一步的研究。