基于VMD和排列熵的滾動軸承故障診斷研究*

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

(1.華東交通大學電氣與自動化工程學院，南昌 330033；2.中國鐵路南昌局集團有限公司，南昌 330033)

0 引言

滾動軸承是旋轉機械的重要組成部分，如果發生故障，在其工作時會造成安全隱患，因此判斷出軸承的當前狀態，并采取相應的措施處理十分必要[1]。

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

本文基于滾動軸承的故障運行機制提出一種變分模態分解和基于峭度準則排列熵構建特征向量的方法，通過SVM進行分類診斷。并通過實例信號進行分析，結果表明提出的方法可以實現滾動軸承的故障診斷。

1 變分模態分解原理及其參數選取

本節介紹變分模態分解原理和步驟，并且分析參數設置對分解結果的影響。

1.1 變分模態分解原理及步驟

1.1.1 變分模態分解原理

變分模態分解就是尋求K個估計帶寬之和最小的模態函數，并且要求所有模態函數之和為原函數[6]，約束變分模型表達式：

(1)

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

求解式(1)變分問題的最優解，引入增廣lagrange函數：

(2)

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

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

1.1.2 變分模態分解步驟

變分模態分解算法的步驟如下：

(2)n=n+1，進入循環；

(3)根據更新公式進行更新mk,,ωk，直至分解個數達到K時停止循環；

(4)根據公式更新λ；

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

1.2 參數設置對分解結果的影響

變分模態分解算法包含的參數有分解尺度K、懲罰因子α、噪聲容限和判別精度，研究發現，噪聲容限和判別精度對變分模態分解的結果影響較小，本小節通過定一求二法分析確定變分模態分解參數，介紹K和α對分解的影響。

本小節采用西儲凱斯大學軸承數據庫數據進行分析，以采樣頻率為12 kHz下的驅動端軸承的滾動體故障數據做分析，圖1為該故障信號的時域圖，橫坐標為時間，縱坐標為幅值。

圖1 軸承滾動體故障時的振動信號

1.2.1 分解個數K對分解結果的影響

首先懲罰因子設定為1000，分別設置K值為2、3、4、5，最終分解個數和中心頻率結果如表1所示。根據表1可知當K值小于3時，分解尺度明顯不足；當K大于4時的中心頻率2698和2840比較接近，很可能存在過分解現象，初步設定分解個數為4。

表1 不同分解個數K所對應的中心頻率

1.2.2 懲罰因子α對分解尺度的影響

根據上面K值的確定，設定K值為4，再設置當α等于500、1000、1500、2000、2500時，各個分量的中心頻率數值如表2所示。

表2 不同懲罰因子α對應的中心頻率

根據表2所示可知α的值小于1000時對于分解結果出現欠分解現象，α的值大于1000時趨于穩定，因此再次設置α的值為600、700、800、900分析各個分量的中心頻率發現α的值在700時不會出現欠分解現象，通過此方法分析滾動體故障信號大約K=4,α=700左右的時候通過變分模態分解算法分解結果最優。

通過上述定一求二的分析最終得到的4種狀態參數組合如表3所示。

表3 定一求二設置得到的參數組合

2 基于峭度準則的排列熵特征向量構建方法

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

2.1 排列熵概念

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

(1)對X(i),i=1,2…,n的一個時間序列進行空間重構，得到矩陣如下所示：

(3)

式中，m為嵌入維數；τ為延遲時間；Z為重構相空間的向量個數。

(2)矩陣中每行可看為一個重構分量，將矩陣中的第j重構分量{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為重構分量中的每個元素所在列的索引。

如果重構分量存在相等值，那么：

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

(5)

則根據ip,iq的大小排序，如果ip

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

(6)

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

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

(7)

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

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

(8)

當pi=1/m！時，Ep(m)就達到最大值lnm！

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

Ep=Ep(m)/lnm!

(9)

2.2 滾動軸承故障運行狀態下模態分量的峭度分析

峭度為描繪波形尖峰度的參數，其數學描述公式為：

(10)

式中，a為信號的均值；σ為信號的標準差。

通過分析VMD分解后模態分量的峭度來分析滾動軸承故障狀態下的運行機理，以西儲凱斯大學數據庫中的內圈故障為例進行分析。

選取內圈故障2048個采樣點，通過上一節對不同狀態下VMD參數組合[6,2000]，經過VMD分解滾動體故障信號得到的6個分解模態，根據峭度表述公式得到的不同模態峭度值如表4所示。

表4 不同模態峭度值

根據上表可以得知經過VMD分解后的模態分量當n=2、4、5、6時峭度值較高。

選取不同采樣起始點，同樣選取2048個采樣點，進行VMD分解得到的模態分量峭度值如圖2所示。

圖2 不同采樣起始點下的模態分量峭度值

根據上圖分析在不同采樣起始點下的峭度值同樣遵循n=2、4、5、6時峭度值較高的規律，此為滾動軸承的故障運行機理。

2.3 基于峭度準則的排列熵特征向量構建方法

軸承無故障運行時，振動信號近似接近正態分布，此時的峭度值近似為3；當軸承發生故障時，其振動信號概率密度增大，信號幅值即偏離正態分布，峭度指標的絕對值越大，軸承的故障就越嚴重。因此，IMF中計算得到峭度絕對值越大，含有故障沖擊的成分就越多[10]。

排列熵計算過程中，延遲時間τ和嵌入維度m的選取對于排列熵的計算結果有一定影響[11]，對于這兩個參數的選取，本文跟經驗選取τ=6，m=1。

由于第1節中分析得到的參數組合中滾動體故障時的K值最小，所以分別取經過VMD分解后內圈故障、外圈故障和正常狀態下模態分量中峭度值較大的4個分量用以構建特征向量，具體做法如下：

(1)通過VMD算法分解振動信號，得到若干IMF分量；

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

(3)分別計算4個分量的排列熵值，以構建特征向量。

通過此方法得到的內圈故障用以構建特征向量的模態分量n為2、4、5、6；外圈用以構建特征向量的模態分量n為2、3、4、5；正常狀態下用以構建特征向量的模態分量n為1、4、5、6。

3 滾動軸承狀態識別流程

本文采用定一求二的方法確定VMD算法的參數組合[K,α]，經過分解得到各個IMF分量，針對故障信號的特點，采用峭度準則和排列熵方法構建特征向量，分析不同故障類型下的排列熵值，確定不同IMF分量下的特征向量，最后將構建的特征向量輸入SVM中進行訓練和模態分類，診斷流程如圖3所示。

圖3 故障診斷流程圖

具體實現步驟為：

(1)獲取實驗數據，載入原始信號；

(2)通過定一求二的方法確定VMD的分解參數，并通過分解信號得到的各個IMF分量；

(3)計算不同狀態類型下不同模態的峭度值，通過K值的局限，來確定用以構建特征向量的模態；

(4)通過計算不同模態的排列熵構建特征向量；

(5)建立SVM模型，將訓練數據和待測數據輸入其中得出結果。

4 實例分析

本文選用西儲凱斯大學實驗室的數據作為實例分析。以采樣頻率為12 kHz下的驅動端軸承故障數據做分析，分別采用損傷直徑為0.177 8 mm的內圈、滾動體、外圈故障和正常狀態下的振動信號，各選取4種狀態振動信號的60組樣本，其中訓練樣本40組，20組作為測試樣本。

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

以不同采樣點為例介紹不同特征向量的構建，分別通過VMD分解內圈故障、外圈故障、滾動體故障以及正常狀態下的軸承振動信號，根據第3節的特征向量構建方法得到的4類特征如表5所示。

表5 特征向量構建值表

表5中標簽1～4分別表示為正常狀態、內圈故障、外圈故障以及滾動體故障滾動軸承狀態類型；由于最小K值為4的限制，特征1～4分別表示為4個的模態分量對應的排列熵值。

由于特征向量數目過多在此不加以展示，將上述計算得到的40組樣本輸入到支持向量機中訓練后，輸入測試樣本得到的識別結果如圖4所示。

圖4 支持向量機測試數據分類圖

根據上圖所示的結果可以看出，4種狀態下測試數據分類的正確率分別為95%、100%、90%和95%，最終診斷的平均正確率為93.75%，因此本文提出的方法可以良好的實現滾動軸承的故障診斷。

5 結論

本文提出變分模態分解、排列熵以及支持向量機結合的滾動軸承故障診斷方法，通過實例信號實驗分析得出結論如下：

(1)變分模態分解的參數設置對于分解結果尤為重要，良好的參數設置對于后期的故障診斷有著重要作用。

(2)通過峭度準則分析了滾動軸承的運行機理，并根據最終的分類結果可以得出利用峭度準則結合VMD和排列熵構建特征向量的方法，可以較好的實現故障診斷。