基于VMD和排列熵的滾動(dòng)軸承故障診斷研究*

楊 云，張昊宇，薛元賀,2，丁 磊

(1.華東交通大學(xué)電氣與自動(dòng)化工程學(xué)院，南昌 330033；2.中國(guó)鐵路南昌局集團(tuán)有限公司，南昌 330033)

0 引言

滾動(dòng)軸承是旋轉(zhuǎn)機(jī)械的重要組成部分，如果發(fā)生故障，在其工作時(shí)會(huì)造成安全隱患，因此判斷出軸承的當(dāng)前狀態(tài)，并采取相應(yīng)的措施處理十分必要[1]。

軸承的振動(dòng)信號(hào)多為非線性、非平穩(wěn)信號(hào)，經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical mode decomposition，EMD)作為早期的自適應(yīng)分解算法與其他方法結(jié)合激起了人們的研究興趣。王金東等[2]通過(guò)EMD分解軸承的振動(dòng)信號(hào)最后結(jié)合SVM，可以較準(zhǔn)確實(shí)現(xiàn)故障診斷；由于EMD算法本身的缺陷，馬麗華等[3]提出基于集合經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)和GG聚類的方法，最終的結(jié)果也證明了該方法的可行性。2005年由EMD理論更新出了一種新的自適應(yīng)分解算法局部均值分解算法(Local mean decomposition，LMD)，侯高雁等[4]針對(duì)滾動(dòng)軸承的振動(dòng)信號(hào)，通過(guò)LMD形態(tài)學(xué)與EEMD形態(tài)學(xué)在故障提取中的對(duì)比研究，表明LMD對(duì)故障提取有著速度快，清晰度高的特點(diǎn)。人們?cè)谘芯窟^(guò)程中，發(fā)現(xiàn)EMD、LMD兩者算法本身存在局限性，無(wú)法解決端點(diǎn)效應(yīng)和模態(tài)混疊現(xiàn)象。Dragomiretskiy Konstantin在2014年提出了變分模態(tài)分解(Variational modal decomposition，VMD)方法[5]，此方法可以避免端點(diǎn)效應(yīng)、抑制模態(tài)變分模態(tài)分解和樣本熵的特征提取方法,采用支持向量機(jī)進(jìn)行故障識(shí)別。該方法能精確的實(shí)現(xiàn)故障診斷。熵作為一種構(gòu)建特征向量的方法也廣泛用于軸承的故障診斷領(lǐng)域[7-8]。

本文基于滾動(dòng)軸承的故障運(yùn)行機(jī)制提出一種變分模態(tài)分解和基于峭度準(zhǔn)則排列熵構(gòu)建特征向量的方法，通過(guò)SVM進(jìn)行分類診斷。并通過(guò)實(shí)例信號(hào)進(jìn)行分析，結(jié)果表明提出的方法可以實(shí)現(xiàn)滾動(dòng)軸承的故障診斷。

1 變分模態(tài)分解原理及其參數(shù)選取

本節(jié)介紹變分模態(tài)分解原理和步驟，并且分析參數(shù)設(shè)置對(duì)分解結(jié)果的影響。

1.1 變分模態(tài)分解原理及步驟

1.1.1 變分模態(tài)分解原理

變分模態(tài)分解就是尋求K個(gè)估計(jì)帶寬之和最小的模態(tài)函數(shù)，并且要求所有模態(tài)函數(shù)之和為原函數(shù)[6]，約束變分模型表達(dá)式：

(1)

式中，mk表示分解得到的K個(gè)IMF分量；ωk表示量的中心頻率。

求解式(1)變分問(wèn)題的最優(yōu)解，引入增廣lagrange函數(shù)：

(2)

式中，α為罰因子；λ為lagrange乘子。

利用交替方向乘子算法(ADMM)求取上述變分問(wèn)題，最終結(jié)束迭代得到K個(gè)IMF分量。

1.1.2 變分模態(tài)分解步驟

變分模態(tài)分解算法的步驟如下：

(2)n=n+1，進(jìn)入循環(huán)；

(3)根據(jù)更新公式進(jìn)行更新mk,,ωk，直至分解個(gè)數(shù)達(dá)到K時(shí)停止循環(huán)；

(4)根據(jù)公式更新λ；

(5)給定精度ε，若滿足停止條件，則停止循環(huán)，否則返回步驟(2)繼續(xù)循環(huán)。

1.2 參數(shù)設(shè)置對(duì)分解結(jié)果的影響

變分模態(tài)分解算法包含的參數(shù)有分解尺度K、懲罰因子α、噪聲容限和判別精度，研究發(fā)現(xiàn)，噪聲容限和判別精度對(duì)變分模態(tài)分解的結(jié)果影響較小，本小節(jié)通過(guò)定一求二法分析確定變分模態(tài)分解參數(shù)，介紹K和α對(duì)分解的影響。

本小節(jié)采用西儲(chǔ)凱斯大學(xué)軸承數(shù)據(jù)庫(kù)數(shù)據(jù)進(jìn)行分析，以采樣頻率為12 kHz下的驅(qū)動(dòng)端軸承的滾動(dòng)體故障數(shù)據(jù)做分析，圖1為該故障信號(hào)的時(shí)域圖，橫坐標(biāo)為時(shí)間，縱坐標(biāo)為幅值。

圖1 軸承滾動(dòng)體故障時(shí)的振動(dòng)信號(hào)

1.2.1 分解個(gè)數(shù)K對(duì)分解結(jié)果的影響

首先懲罰因子設(shè)定為1000，分別設(shè)置K值為2、3、4、5，最終分解個(gè)數(shù)和中心頻率結(jié)果如表1所示。根據(jù)表1可知當(dāng)K值小于3時(shí)，分解尺度明顯不足；當(dāng)K大于4時(shí)的中心頻率2698和2840比較接近，很可能存在過(guò)分解現(xiàn)象，初步設(shè)定分解個(gè)數(shù)為4。

表1 不同分解個(gè)數(shù)K所對(duì)應(yīng)的中心頻率

1.2.2 懲罰因子α對(duì)分解尺度的影響

根據(jù)上面K值的確定，設(shè)定K值為4，再設(shè)置當(dāng)α等于500、1000、1500、2000、2500時(shí)，各個(gè)分量的中心頻率數(shù)值如表2所示。

表2 不同懲罰因子α對(duì)應(yīng)的中心頻率

根據(jù)表2所示可知α的值小于1000時(shí)對(duì)于分解結(jié)果出現(xiàn)欠分解現(xiàn)象，α的值大于1000時(shí)趨于穩(wěn)定，因此再次設(shè)置α的值為600、700、800、900分析各個(gè)分量的中心頻率發(fā)現(xiàn)α的值在700時(shí)不會(huì)出現(xiàn)欠分解現(xiàn)象，通過(guò)此方法分析滾動(dòng)體故障信號(hào)大約K=4,α=700左右的時(shí)候通過(guò)變分模態(tài)分解算法分解結(jié)果最優(yōu)。

通過(guò)上述定一求二的分析最終得到的4種狀態(tài)參數(shù)組合如表3所示。

表3 定一求二設(shè)置得到的參數(shù)組合

2 基于峭度準(zhǔn)則的排列熵特征向量構(gòu)建方法

考慮到軸承故障時(shí)故障沖擊隨時(shí)間存在周期性，本文提出了基于峭度準(zhǔn)則的排列熵特征向量構(gòu)建方法，下面分別介紹其概念。

2.1 排列熵概念

排列熵是衡量以為時(shí)間序列復(fù)雜程度的熵，具體原理如下[9]：

(1)對(duì)X(i),i=1,2…,n的一個(gè)時(shí)間序列進(jìn)行空間重構(gòu)，得到矩陣如下所示：

(3)

式中，m為嵌入維數(shù)；τ為延遲時(shí)間；Z為重構(gòu)相空間的向量個(gè)數(shù)。

(2)矩陣中每行可看為一個(gè)重構(gòu)分量，將矩陣中的第j重構(gòu)分量{x(j),x(j+τ),…,x(j+(x(j+(m-1)τ)}以升序的方法排列，得到：

x(j+(i1-1)τ)≤x(j+(i2-1)τ)≤…≤

x(j+(im-1)τ)

(4)

式中，i1,i2,…,im為重構(gòu)分量中的每個(gè)元素所在列的索引。

如果重構(gòu)分量存在相等值，那么：

x(j-(ip-1)τ)=x(j-(iq-1)τ)

(5)

則根據(jù)ip,iq的大小排序，如果ip

x(j-(ip-1)τ)≤x(j-(iq-1)τ)

(6)

因此，對(duì)于重構(gòu)矩陣Y的任一重構(gòu)分量Y(j)都將得到一組位置索引序列：

B(j)=(i1,i2,…,im)，j=1,2,3,…,k

(7)

式中，k≤m！，B(j)是符號(hào)序列的其中一種。

(3)算出每一位置索引序列出現(xiàn)的概率p1，p2…，pk，時(shí)間序列X={x(i)，i=1，2，3…，n}的不同位置索引序列的排列熵可以定義為：

(8)

當(dāng)pi=1/m！時(shí)，Ep(m)就達(dá)到最大值lnm！

(4)為了方便各運(yùn)行狀態(tài)下的排列熵比較，通常用lnm!將Ep(m)進(jìn)行歸一化處理，即：

Ep=Ep(m)/lnm!

(9)

2.2 滾動(dòng)軸承故障運(yùn)行狀態(tài)下模態(tài)分量的峭度分析

峭度為描繪波形尖峰度的參數(shù)，其數(shù)學(xué)描述公式為：

(10)

式中，a為信號(hào)的均值；σ為信號(hào)的標(biāo)準(zhǔn)差。

通過(guò)分析VMD分解后模態(tài)分量的峭度來(lái)分析滾動(dòng)軸承故障狀態(tài)下的運(yùn)行機(jī)理，以西儲(chǔ)凱斯大學(xué)數(shù)據(jù)庫(kù)中的內(nèi)圈故障為例進(jìn)行分析。

選取內(nèi)圈故障2048個(gè)采樣點(diǎn)，通過(guò)上一節(jié)對(duì)不同狀態(tài)下VMD參數(shù)組合[6,2000]，經(jīng)過(guò)VMD分解滾動(dòng)體故障信號(hào)得到的6個(gè)分解模態(tài)，根據(jù)峭度表述公式得到的不同模態(tài)峭度值如表4所示。

表4 不同模態(tài)峭度值

根據(jù)上表可以得知經(jīng)過(guò)VMD分解后的模態(tài)分量當(dāng)n=2、4、5、6時(shí)峭度值較高。

選取不同采樣起始點(diǎn)，同樣選取2048個(gè)采樣點(diǎn)，進(jìn)行VMD分解得到的模態(tài)分量峭度值如圖2所示。

圖2 不同采樣起始點(diǎn)下的模態(tài)分量峭度值

根據(jù)上圖分析在不同采樣起始點(diǎn)下的峭度值同樣遵循n=2、4、5、6時(shí)峭度值較高的規(guī)律，此為滾動(dòng)軸承的故障運(yùn)行機(jī)理。

2.3 基于峭度準(zhǔn)則的排列熵特征向量構(gòu)建方法

軸承無(wú)故障運(yùn)行時(shí)，振動(dòng)信號(hào)近似接近正態(tài)分布，此時(shí)的峭度值近似為3；當(dāng)軸承發(fā)生故障時(shí)，其振動(dòng)信號(hào)概率密度增大，信號(hào)幅值即偏離正態(tài)分布，峭度指標(biāo)的絕對(duì)值越大，軸承的故障就越嚴(yán)重。因此，IMF中計(jì)算得到峭度絕對(duì)值越大，含有故障沖擊的成分就越多[10]。

排列熵計(jì)算過(guò)程中，延遲時(shí)間τ和嵌入維度m的選取對(duì)于排列熵的計(jì)算結(jié)果有一定影響[11]，對(duì)于這兩個(gè)參數(shù)的選取，本文跟經(jīng)驗(yàn)選取τ=6，m=1。

由于第1節(jié)中分析得到的參數(shù)組合中滾動(dòng)體故障時(shí)的K值最小，所以分別取經(jīng)過(guò)VMD分解后內(nèi)圈故障、外圈故障和正常狀態(tài)下模態(tài)分量中峭度值較大的4個(gè)分量用以構(gòu)建特征向量，具體做法如下：

(1)通過(guò)VMD算法分解振動(dòng)信號(hào)，得到若干IMF分量；

(2)計(jì)算各個(gè)分量的峭度值，選取較大的4個(gè)分量；

(3)分別計(jì)算4個(gè)分量的排列熵值，以構(gòu)建特征向量。

通過(guò)此方法得到的內(nèi)圈故障用以構(gòu)建特征向量的模態(tài)分量n為2、4、5、6；外圈用以構(gòu)建特征向量的模態(tài)分量n為2、3、4、5；正常狀態(tài)下用以構(gòu)建特征向量的模態(tài)分量n為1、4、5、6。

3 滾動(dòng)軸承狀態(tài)識(shí)別流程

本文采用定一求二的方法確定VMD算法的參數(shù)組合[K,α]，經(jīng)過(guò)分解得到各個(gè)IMF分量，針對(duì)故障信號(hào)的特點(diǎn)，采用峭度準(zhǔn)則和排列熵方法構(gòu)建特征向量，分析不同故障類型下的排列熵值，確定不同IMF分量下的特征向量，最后將構(gòu)建的特征向量輸入SVM中進(jìn)行訓(xùn)練和模態(tài)分類，診斷流程如圖3所示。

圖3 故障診斷流程圖

具體實(shí)現(xiàn)步驟為：

(1)獲取實(shí)驗(yàn)數(shù)據(jù)，載入原始信號(hào)；

(2)通過(guò)定一求二的方法確定VMD的分解參數(shù)，并通過(guò)分解信號(hào)得到的各個(gè)IMF分量；

(3)計(jì)算不同狀態(tài)類型下不同模態(tài)的峭度值，通過(guò)K值的局限，來(lái)確定用以構(gòu)建特征向量的模態(tài)；

(4)通過(guò)計(jì)算不同模態(tài)的排列熵構(gòu)建特征向量；

(5)建立SVM模型，將訓(xùn)練數(shù)據(jù)和待測(cè)數(shù)據(jù)輸入其中得出結(jié)果。

4 實(shí)例分析

本文選用西儲(chǔ)凱斯大學(xué)實(shí)驗(yàn)室的數(shù)據(jù)作為實(shí)例分析。以采樣頻率為12 kHz下的驅(qū)動(dòng)端軸承故障數(shù)據(jù)做分析，分別采用損傷直徑為0.177 8 mm的內(nèi)圈、滾動(dòng)體、外圈故障和正常狀態(tài)下的振動(dòng)信號(hào)，各選取4種狀態(tài)振動(dòng)信號(hào)的60組樣本，其中訓(xùn)練樣本40組，20組作為測(cè)試樣本。

由于定一求二方法得到的參數(shù)組合(見(jiàn)表3)下4種狀態(tài)中最小K值為4的限制，分別選取不同狀態(tài)下的4個(gè)模態(tài)分量計(jì)算排列熵構(gòu)建特征向量。

以不同采樣點(diǎn)為例介紹不同特征向量的構(gòu)建，分別通過(guò)VMD分解內(nèi)圈故障、外圈故障、滾動(dòng)體故障以及正常狀態(tài)下的軸承振動(dòng)信號(hào)，根據(jù)第3節(jié)的特征向量構(gòu)建方法得到的4類特征如表5所示。

表5 特征向量構(gòu)建值表

表5中標(biāo)簽1～4分別表示為正常狀態(tài)、內(nèi)圈故障、外圈故障以及滾動(dòng)體故障滾動(dòng)軸承狀態(tài)類型；由于最小K值為4的限制，特征1～4分別表示為4個(gè)的模態(tài)分量對(duì)應(yīng)的排列熵值。

由于特征向量數(shù)目過(guò)多在此不加以展示，將上述計(jì)算得到的40組樣本輸入到支持向量機(jī)中訓(xùn)練后，輸入測(cè)試樣本得到的識(shí)別結(jié)果如圖4所示。

圖4 支持向量機(jī)測(cè)試數(shù)據(jù)分類圖

根據(jù)上圖所示的結(jié)果可以看出，4種狀態(tài)下測(cè)試數(shù)據(jù)分類的正確率分別為95%、100%、90%和95%，最終診斷的平均正確率為93.75%，因此本文提出的方法可以良好的實(shí)現(xiàn)滾動(dòng)軸承的故障診斷。

5 結(jié)論

本文提出變分模態(tài)分解、排列熵以及支持向量機(jī)結(jié)合的滾動(dòng)軸承故障診斷方法，通過(guò)實(shí)例信號(hào)實(shí)驗(yàn)分析得出結(jié)論如下：

(1)變分模態(tài)分解的參數(shù)設(shè)置對(duì)于分解結(jié)果尤為重要，良好的參數(shù)設(shè)置對(duì)于后期的故障診斷有著重要作用。

(2)通過(guò)峭度準(zhǔn)則分析了滾動(dòng)軸承的運(yùn)行機(jī)理，并根據(jù)最終的分類結(jié)果可以得出利用峭度準(zhǔn)則結(jié)合VMD和排列熵構(gòu)建特征向量的方法，可以較好的實(shí)現(xiàn)故障診斷。