楊 云,張昊宇,薛元賀,2,丁 磊
(1.華東交通大學(xué)電氣與自動化工程學(xué)院,南昌 330033;2.中國鐵路南昌局集團(tuán)有限公司,南昌 330033)
滾動軸承是旋轉(zhuǎn)機(jī)械的重要組成部分,如果發(fā)生故障,在其工作時會造成安全隱患,因此判斷出軸承的當(dāng)前狀態(tài),并采取相應(yīng)的措施處理十分必要[1]。
軸承的振動信號多為非線性、非平穩(wěn)信號,經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical mode decomposition,EMD)作為早期的自適應(yīng)分解算法與其他方法結(jié)合激起了人們的研究興趣。王金東等[2]通過EMD分解軸承的振動信號最后結(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]針對滾動軸承的振動信號,通過LMD形態(tài)學(xué)與EEMD形態(tài)學(xué)在故障提取中的對比研究,表明LMD對故障提取有著速度快,清晰度高的特點(diǎn)。人們在研究過程中,發(fā)現(xiàn)EMD、LMD兩者算法本身存在局限性,無法解決端點(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í)現(xiàn)故障診斷。熵作為一種構(gòu)建特征向量的方法也廣泛用于軸承的故障診斷領(lǐng)域[7-8]。
本文基于滾動軸承的故障運(yùn)行機(jī)制提出一種變分模態(tài)分解和基于峭度準(zhǔn)則排列熵構(gòu)建特征向量的方法,通過SVM進(jìn)行分類診斷。并通過實(shí)例信號進(jìn)行分析,結(jié)果表明提出的方法可以實(shí)現(xiàn)滾動軸承的故障診斷。
本節(jié)介紹變分模態(tài)分解原理和步驟,并且分析參數(shù)設(shè)置對分解結(jié)果的影響。
1.1.1 變分模態(tài)分解原理
變分模態(tài)分解就是尋求K個估計(jì)帶寬之和最小的模態(tài)函數(shù),并且要求所有模態(tài)函數(shù)之和為原函數(shù)[6],約束變分模型表達(dá)式:
(1)
式中,mk表示分解得到的K個IMF分量;ωk表示量的中心頻率。
求解式(1)變分問題的最優(yōu)解,引入增廣lagrange函數(shù):
(2)
式中,α為罰因子;λ為lagrange乘子。
利用交替方向乘子算法(ADMM)求取上述變分問題,最終結(jié)束迭代得到K個IMF分量。
1.1.2 變分模態(tài)分解步驟
變分模態(tài)分解算法的步驟如下:
(2)n=n+1,進(jìn)入循環(huán);
(3)根據(jù)更新公式進(jìn)行更新mk,,ωk,直至分解個數(shù)達(dá)到K時停止循環(huán);
(4)根據(jù)公式更新λ;
(5)給定精度ε,若滿足停止條件,則停止循環(huán),否則返回步驟(2)繼續(xù)循環(huán)。
變分模態(tài)分解算法包含的參數(shù)有分解尺度K、懲罰因子α、噪聲容限和判別精度,研究發(fā)現(xiàn),噪聲容限和判別精度對變分模態(tài)分解的結(jié)果影響較小,本小節(jié)通過定一求二法分析確定變分模態(tài)分解參數(shù),介紹K和α對分解的影響。
本小節(jié)采用西儲凱斯大學(xué)軸承數(shù)據(jù)庫數(shù)據(jù)進(jìn)行分析,以采樣頻率為12 kHz下的驅(qū)動端軸承的滾動體故障數(shù)據(jù)做分析,圖1為該故障信號的時域圖,橫坐標(biāo)為時間,縱坐標(biāo)為幅值。
圖1 軸承滾動體故障時的振動信號
1.2.1 分解個數(shù)K對分解結(jié)果的影響
首先懲罰因子設(shè)定為1000,分別設(shè)置K值為2、3、4、5,最終分解個數(shù)和中心頻率結(jié)果如表1所示。根據(jù)表1可知當(dāng)K值小于3時,分解尺度明顯不足;當(dāng)K大于4時的中心頻率2698和2840比較接近,很可能存在過分解現(xiàn)象,初步設(shè)定分解個數(shù)為4。
表1 不同分解個數(shù)K所對應(yīng)的中心頻率
1.2.2 懲罰因子α對分解尺度的影響
根據(jù)上面K值的確定,設(shè)定K值為4,再設(shè)置當(dāng)α等于500、1000、1500、2000、2500時,各個分量的中心頻率數(shù)值如表2所示。
表2 不同懲罰因子α對應(yīng)的中心頻率
根據(jù)表2所示可知α的值小于1000時對于分解結(jié)果出現(xiàn)欠分解現(xiàn)象,α的值大于1000時趨于穩(wěn)定,因此再次設(shè)置α的值為600、700、800、900分析各個分量的中心頻率發(fā)現(xiàn)α的值在700時不會出現(xiàn)欠分解現(xiàn)象,通過此方法分析滾動體故障信號大約K=4,α=700左右的時候通過變分模態(tài)分解算法分解結(jié)果最優(yōu)。
通過上述定一求二的分析最終得到的4種狀態(tài)參數(shù)組合如表3所示。
表3 定一求二設(shè)置得到的參數(shù)組合
考慮到軸承故障時故障沖擊隨時間存在周期性,本文提出了基于峭度準(zhǔn)則的排列熵特征向量構(gòu)建方法,下面分別介紹其概念。
排列熵是衡量以為時間序列復(fù)雜程度的熵,具體原理如下[9]:
(1)對X(i),i=1,2…,n的一個時間序列進(jìn)行空間重構(gòu),得到矩陣如下所示:
(3)
式中,m為嵌入維數(shù);τ為延遲時間;Z為重構(gòu)相空間的向量個數(shù)。
(2)矩陣中每行可看為一個重構(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òu)分量存在相等值,那么:
x(j-(ip-1)τ)=x(j-(iq-1)τ)
(5)
則根據(jù)ip,iq的大小排序,如果ip x(j-(ip-1)τ)≤x(j-(iq-1)τ) (6) 因此,對于重構(gòu)矩陣Y的任一重構(gòu)分量Y(j)都將得到一組位置索引序列: B(j)=(i1,i2,…,im),j=1,2,3,…,k (7) 式中,k≤m!,B(j)是符號序列的其中一種。 (3)算出每一位置索引序列出現(xiàn)的概率p1,p2…,pk,時間序列X={x(i),i=1,2,3…,n}的不同位置索引序列的排列熵可以定義為: (8) 當(dāng)pi=1/m!時,Ep(m)就達(dá)到最大值lnm! (4)為了方便各運(yùn)行狀態(tài)下的排列熵比較,通常用lnm!將Ep(m)進(jìn)行歸一化處理,即: Ep=Ep(m)/lnm! (9) 峭度為描繪波形尖峰度的參數(shù),其數(shù)學(xué)描述公式為: (10) 式中,a為信號的均值;σ為信號的標(biāo)準(zhǔn)差。 通過分析VMD分解后模態(tài)分量的峭度來分析滾動軸承故障狀態(tài)下的運(yùn)行機(jī)理,以西儲凱斯大學(xué)數(shù)據(jù)庫中的內(nèi)圈故障為例進(jìn)行分析。 選取內(nèi)圈故障2048個采樣點(diǎn),通過上一節(jié)對不同狀態(tài)下VMD參數(shù)組合[6,2000],經(jīng)過VMD分解滾動體故障信號得到的6個分解模態(tài),根據(jù)峭度表述公式得到的不同模態(tài)峭度值如表4所示。 表4 不同模態(tài)峭度值 根據(jù)上表可以得知經(jīng)過VMD分解后的模態(tài)分量當(dāng)n=2、4、5、6時峭度值較高。 選取不同采樣起始點(diǎn),同樣選取2048個采樣點(diǎn),進(jìn)行VMD分解得到的模態(tài)分量峭度值如圖2所示。 圖2 不同采樣起始點(diǎn)下的模態(tài)分量峭度值 根據(jù)上圖分析在不同采樣起始點(diǎn)下的峭度值同樣遵循n=2、4、5、6時峭度值較高的規(guī)律,此為滾動軸承的故障運(yùn)行機(jī)理。 軸承無故障運(yùn)行時,振動信號近似接近正態(tài)分布,此時的峭度值近似為3;當(dāng)軸承發(fā)生故障時,其振動信號概率密度增大,信號幅值即偏離正態(tài)分布,峭度指標(biāo)的絕對值越大,軸承的故障就越嚴(yán)重。因此,IMF中計(jì)算得到峭度絕對值越大,含有故障沖擊的成分就越多[10]。 排列熵計(jì)算過程中,延遲時間τ和嵌入維度m的選取對于排列熵的計(jì)算結(jié)果有一定影響[11],對于這兩個參數(shù)的選取,本文跟經(jīng)驗(yàn)選取τ=6,m=1。 由于第1節(jié)中分析得到的參數(shù)組合中滾動體故障時的K值最小,所以分別取經(jīng)過VMD分解后內(nèi)圈故障、外圈故障和正常狀態(tài)下模態(tài)分量中峭度值較大的4個分量用以構(gòu)建特征向量,具體做法如下: (1)通過VMD算法分解振動信號,得到若干IMF分量; (2)計(jì)算各個分量的峭度值,選取較大的4個分量; (3)分別計(jì)算4個分量的排列熵值,以構(gòu)建特征向量。 通過此方法得到的內(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。 本文采用定一求二的方法確定VMD算法的參數(shù)組合[K,α],經(jīng)過分解得到各個IMF分量,針對故障信號的特點(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ù),載入原始信號; (2)通過定一求二的方法確定VMD的分解參數(shù),并通過分解信號得到的各個IMF分量; (3)計(jì)算不同狀態(tài)類型下不同模態(tài)的峭度值,通過K值的局限,來確定用以構(gòu)建特征向量的模態(tài); (4)通過計(jì)算不同模態(tài)的排列熵構(gòu)建特征向量; (5)建立SVM模型,將訓(xùn)練數(shù)據(jù)和待測數(shù)據(jù)輸入其中得出結(jié)果。 本文選用西儲凱斯大學(xué)實(shí)驗(yàn)室的數(shù)據(jù)作為實(shí)例分析。以采樣頻率為12 kHz下的驅(qū)動端軸承故障數(shù)據(jù)做分析,分別采用損傷直徑為0.177 8 mm的內(nèi)圈、滾動體、外圈故障和正常狀態(tài)下的振動信號,各選取4種狀態(tài)振動信號的60組樣本,其中訓(xùn)練樣本40組,20組作為測試樣本。 由于定一求二方法得到的參數(shù)組合(見表3)下4種狀態(tài)中最小K值為4的限制,分別選取不同狀態(tài)下的4個模態(tài)分量計(jì)算排列熵構(gòu)建特征向量。 以不同采樣點(diǎn)為例介紹不同特征向量的構(gòu)建,分別通過VMD分解內(nèi)圈故障、外圈故障、滾動體故障以及正常狀態(tài)下的軸承振動信號,根據(jù)第3節(jié)的特征向量構(gòu)建方法得到的4類特征如表5所示。 表5 特征向量構(gòu)建值表 表5中標(biāo)簽1~4分別表示為正常狀態(tài)、內(nèi)圈故障、外圈故障以及滾動體故障滾動軸承狀態(tài)類型;由于最小K值為4的限制,特征1~4分別表示為4個的模態(tài)分量對應(yīng)的排列熵值。 由于特征向量數(shù)目過多在此不加以展示,將上述計(jì)算得到的40組樣本輸入到支持向量機(jī)中訓(xùn)練后,輸入測試樣本得到的識別結(jié)果如圖4所示。 圖4 支持向量機(jī)測試數(shù)據(jù)分類圖 根據(jù)上圖所示的結(jié)果可以看出,4種狀態(tài)下測試數(shù)據(jù)分類的正確率分別為95%、100%、90%和95%,最終診斷的平均正確率為93.75%,因此本文提出的方法可以良好的實(shí)現(xiàn)滾動軸承的故障診斷。 本文提出變分模態(tài)分解、排列熵以及支持向量機(jī)結(jié)合的滾動軸承故障診斷方法,通過實(shí)例信號實(shí)驗(yàn)分析得出結(jié)論如下: (1)變分模態(tài)分解的參數(shù)設(shè)置對于分解結(jié)果尤為重要,良好的參數(shù)設(shè)置對于后期的故障診斷有著重要作用。 (2)通過峭度準(zhǔn)則分析了滾動軸承的運(yùn)行機(jī)理,并根據(jù)最終的分類結(jié)果可以得出利用峭度準(zhǔn)則結(jié)合VMD和排列熵構(gòu)建特征向量的方法,可以較好的實(shí)現(xiàn)故障診斷。2.2 滾動軸承故障運(yùn)行狀態(tài)下模態(tài)分量的峭度分析
2.3 基于峭度準(zhǔn)則的排列熵特征向量構(gòu)建方法
3 滾動軸承狀態(tài)識別流程
4 實(shí)例分析
5 結(jié)論