基于參數(shù)自適應(yīng)的VMD滾動軸承故障診斷

李永琪，彭珍瑞

（蘭州交通大學(xué) 機(jī)電工程學(xué)院，蘭州730070）

滾動軸承在旋轉(zhuǎn)機(jī)械運(yùn)行中起著舉足輕重的重要作用，其動力學(xué)系統(tǒng)十分復(fù)雜。實(shí)際振動信號往往呈非線性、非平穩(wěn)性[1]。因而如何從這些故障信號中提取出含豐富故障信息的故障特征已成為軸承故障診斷中關(guān)鍵問題之一。

經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical mode decomposition，EMD）由Huang等[2]提出后，主要針對非平穩(wěn)和非線性信號，該方法對信號分解時(shí)可以實(shí)現(xiàn)自適應(yīng)分解，但過包絡(luò)、模態(tài)混疊等缺陷仍存在其中。Smith[3]為解決EMD 擬合問題提出了局部均值分解（Local mean decomposition，LMD）方法，可避免過包絡(luò)缺陷，而且降低了模態(tài)混淆及端點(diǎn)效應(yīng)的影響，但是對信號分解而言仍屬遞歸模式分解，在分解過程中仍會存在誤差。2014年，Dragomiretskiy 等[4]提出了變分模態(tài)分解（Variational mode decomposition，VMD），該方法在處理模式混疊方面具有堅(jiān)實(shí)的數(shù)學(xué)理論基礎(chǔ)，能有效地避免和抑制EMD、LMD的模態(tài)混疊及端點(diǎn)效應(yīng)且分解效率較高[5]，但是會存在模態(tài)分量個(gè)數(shù)K和懲罰因子α兩個(gè)參數(shù)如何合理選取的問題[6]。多數(shù)研究中使用VMD時(shí)，兩個(gè)參數(shù)取決于人為經(jīng)驗(yàn)選取，會嚴(yán)重影響其分解[7]，甚至可能會產(chǎn)生虛假分量且出現(xiàn)錯(cuò)誤診斷。鄭小霞等[8]通過觀察中心頻率選取K值，α取默認(rèn)值2 000而未優(yōu)化，沒有足夠的依據(jù)而缺乏嚴(yán)謹(jǐn)性。張俊等[9]利用粒子群算法優(yōu)化了其參數(shù)，但會增加計(jì)算時(shí)間且存在局部最優(yōu)現(xiàn)象。

針對上述問題，VMD參數(shù)難以合理確定且為了綜合考慮故障脈沖的周期性及循環(huán)平穩(wěn)性，本文提出了一種平均包絡(luò)譜峭度結(jié)合平均樣本熵優(yōu)化的VMD 及加權(quán)合成峭度提取最優(yōu)IMF 的軸承故障診斷方法，逐個(gè)優(yōu)化VMD 參數(shù)，得到最優(yōu)[K,α]，再對信號進(jìn)行分解，選擇加權(quán)合成峭度值最大IMF 作為最優(yōu)IMF，包絡(luò)譜分析后判別故障類型。經(jīng)軸承內(nèi)圈故障仿真信號及實(shí)際軸承內(nèi)、外圈數(shù)據(jù)驗(yàn)證，表明本文方法可有效實(shí)現(xiàn)軸承內(nèi)、外圈故障類型的判斷。

1 基本理論

1.1 變分模態(tài)分解

VMD方法摒棄了EMD中模態(tài)分量循環(huán)篩選過程，將固有模態(tài)函數(shù)重新定義為幅頻調(diào)制信號，其表達(dá)式如下：

式中：uk(t)被認(rèn)為是諧波信號振幅和頻率；Ak(t)是uk(t)瞬時(shí)幅值；瞬時(shí)頻率為相位φk(t)1 階導(dǎo)數(shù)。VMD的約束模型式公式如下：

式中：原始振動信號f(t)是分解得到的各IMF 的累加；{uk}={u1,…,uK}是指原始信號被分解后得到的模態(tài)分量；{wk}={w1,…,wK}表示分解得到的每個(gè)IMF 所對應(yīng)的中心頻率；?(t)是對t求偏導(dǎo)；δ(t)為脈沖函數(shù)。

為了得到變分問題的最優(yōu)解，增廣拉格朗日函數(shù)方法被引入。其公式如下：

式中：λ(t)為拉格朗日乘子；α為二次懲罰因子。用交替方向乘子算法通過迭代更新un+1，ωn+1，λn+1可求取式(3)的“鞍點(diǎn)”，即式(2)的最優(yōu)解。

詳細(xì)的分解過程參考文獻(xiàn)[4]。

1.2 平均包絡(luò)譜峭度

包絡(luò)譜峭度[10]可用于測量和評價(jià)信號中的循環(huán)平穩(wěn)性，在機(jī)械故障診斷中被廣泛使用。在對VMD方法參數(shù)優(yōu)化時(shí)，通過平均包絡(luò)譜峭度來總體評價(jià)IMF 分量的故障成分。Wang 等[11]以小波包節(jié)點(diǎn)信號的包絡(luò)譜峭度最大來選擇最優(yōu)解調(diào)頻帶。因此運(yùn)用平均包絡(luò)譜峭度最大來構(gòu)建目標(biāo)函數(shù)來總體評估IMF分量，具體計(jì)算公式如下：

式中：ESK為包絡(luò)譜峭度；NSE為包絡(luò)譜的采樣數(shù)；SE為信號包絡(luò)譜；K為VMD 分解模態(tài)個(gè)數(shù)；p1為平均包絡(luò)譜峭度。

1.3 平均樣本熵

Richman 等[12]提出樣本熵概念。假如涵蓋越多故障信息，出現(xiàn)與故障相關(guān)的周期性沖擊脈沖，信號的干擾和噪聲較弱，呈現(xiàn)出較強(qiáng)的規(guī)律性，樣本熵值越小表明信號周期性沖擊越明顯[13]。因此運(yùn)用平均樣本熵值最小來構(gòu)建目標(biāo)函數(shù)來總體評估IMF 分量，樣本熵具體公式如下：

式中：SampEn 為樣本熵；m是重構(gòu)相空間后矩陣維數(shù)；r是相似容限閾值；B(r)m為兩個(gè)序列在相似容限下匹配m個(gè)點(diǎn)的概率；B(r)m+1為兩個(gè)序列在相似容限下匹配m+1 個(gè)點(diǎn)的概率；P為給定的時(shí)間序列{x(n)}=x(1),x(2),…x(P)的個(gè)數(shù)；p2為平均樣本熵；K為VMD分解模態(tài)個(gè)數(shù)。

1.4 加權(quán)合成峭度

峭度是稀疏度測量中較為重要度量指標(biāo)之一，它可以檢測信號周期脈沖性，在故障診斷中使用廣泛。其計(jì)算公式如下：

式中：Kr為峭度；μ是信號xi的均值；L是采樣長度；σ是xi的標(biāo)準(zhǔn)差。

然而，這些基于峭度的方法往往僅偏向于單個(gè)脈沖頻帶，卻不是所期望的缺陷脈沖，因它們僅僅重視了故障脈沖的沖擊性，而忽略了其循環(huán)平穩(wěn)性。合成峭度[10]可有效解決峭度的這個(gè)缺陷，通過它們之間期望故障脈沖相比較，發(fā)現(xiàn)僅合成峭度對故障脈沖識別率較高，可將其用于脈沖故障識別?？紤]到峭度及相關(guān)系數(shù)各自優(yōu)缺點(diǎn)，且為了綜合考慮脈沖信號周期的沖擊性和循環(huán)平穩(wěn)性，將合成峭度與相關(guān)系數(shù)結(jié)合，構(gòu)建了加權(quán)合成峭度EKCI進(jìn)行VMD的故障特征提取的方法。相關(guān)公式如下：

式中：Kr、EK、ESK分別表示峭度、合成峭度、包絡(luò)幅值峭度；EKCI為加權(quán)合成峭度指標(biāo)；C為信號x與y相關(guān)系數(shù)；和分別為x及y的均值；E[·]為數(shù)學(xué)期望。

2 參數(shù)自適應(yīng)的VMD 滾動軸承故障診斷

基于以上分析，采用平均包絡(luò)譜峭度及平均樣本熵可以總體評價(jià)其故障成分來優(yōu)化VMD參數(shù)，利用加權(quán)合成峭度可以綜合考慮故障脈沖周期性及循環(huán)平穩(wěn)性來優(yōu)選IMF并提取故障特征。同時(shí)為避免在優(yōu)化中某一指標(biāo)信息被淹沒的情況，故將平均包絡(luò)譜峭度、平均樣本熵分別進(jìn)行歸一化處理并求和，即可得到目標(biāo)函數(shù)式，如式(12)所示，從而構(gòu)建綜合評價(jià)指標(biāo)f。本文以綜合評價(jià)指標(biāo)最小為原則，表明各項(xiàng)指標(biāo)都較好，由節(jié)1.2 及1.3 可知平均樣本熵以最小來構(gòu)建目標(biāo)函數(shù)來總體評估IMF 分量，而平均包絡(luò)譜峭度以最大來構(gòu)建目標(biāo)函數(shù)來總體評估IMF分量，為了滿足綜合指標(biāo)的整體最小的要求，故在構(gòu)建綜合指標(biāo)時(shí)要取平均包絡(luò)譜峭度的倒數(shù)，再歸一化處理并求和。平均包絡(luò)譜峭度的倒數(shù)公式如式(13)所示：

式中：f是兩個(gè)指標(biāo)整體表現(xiàn)；i=1，f1=1/p1為平均包絡(luò)譜峭度的倒數(shù)；i=2，f2=p2為平均樣本熵。

根據(jù)李華等[14]提出的按信息熵最小來逐個(gè)優(yōu)化VMD 參數(shù)方法，先給定α=2 000，再利用f最小值來確定K；然后利用選定好的K，接著再對α進(jìn)行優(yōu)化。整體流程圖如圖1所示。具體過程如下：

圖1 VMD參數(shù)優(yōu)化及故障診斷流程圖

(1)采集原始信號；

(2)先優(yōu)化K。設(shè)置K=2及其數(shù)值范圍[2，12]，1 是其搜索步長；另外，噪聲容忍度t及α取給定值，α=2 000，t=1×10-7。然后對原信號進(jìn)行VMD 分解，判斷綜合評價(jià)指標(biāo)的值是否為最小值，如果是，則確定模態(tài)個(gè)數(shù)K=K；若不是，則K=K+1，否則繼續(xù)判斷，直至達(dá)到結(jié)束條件為止；

(3) 接著對懲罰因子α再進(jìn)行優(yōu)化。由之前優(yōu)化的K再對α進(jìn)行優(yōu)化，初始化α=1 000，搜索范圍[1 000，3 000]，步長50，判斷條件與之前準(zhǔn)則一致，直至達(dá)到結(jié)束要求為止。循環(huán)結(jié)束可得到最優(yōu)[K，α]；

(4)用前一步優(yōu)化參數(shù)組合[K，α]進(jìn)行VMD 分解，得到K個(gè)IMF分量；

(5)計(jì)算各IMF分量的加權(quán)合成峭度指標(biāo)值，比較后選取最大值即最優(yōu)IMF；

(6)用包絡(luò)譜對其進(jìn)行分析，提取故障特征并輸出故障類型。

3 仿真分析

為了驗(yàn)證本文所提方法，用軸承內(nèi)圈故障仿真信號[15]進(jìn)行驗(yàn)證，其仿真信號如下：

式中：周期沖擊信號s(t)；高斯白噪聲n(t)，對s(t)加入-10 dB的高斯白噪聲，模擬實(shí)際工況條件；Ai為第i次沖擊幅值；初始幅值A(chǔ)0=0.3；衰減系數(shù)g=700；fr為軸的轉(zhuǎn)頻為軸承內(nèi)圈故障頻率；fn為系統(tǒng)共振頻率，fn=3000 Hz；采樣頻率fs=12 kHz；采樣點(diǎn)數(shù)N=9 600。加入噪聲后信號信噪比為-10 dB，其均方誤差為0.300 7。圖2 是s(t)沖擊信號，圖3 是仿真信號x(t)，圖4 是圖3 的包絡(luò)譜圖。盡管可以從仿真信號包絡(luò)譜圖中看到故障特征頻率，但其余周邊均有噪聲圍繞。

圖2 沖擊信號

圖3 仿真信號

圖4 仿真信號包絡(luò)譜圖

3.1 VMD參數(shù)優(yōu)化

對上述信號進(jìn)行VMD分解，其模態(tài)數(shù)K與綜合評價(jià)指標(biāo)f關(guān)系如圖5 所示。當(dāng)K=2，綜合指標(biāo)f最小，最小值是1.148，故取模態(tài)個(gè)數(shù)K=2；懲罰因子α與綜合指標(biāo)f的關(guān)系如圖6所示。當(dāng)α=1 000時(shí)，f取得最小值，其值為2，則α=1 000，所以確定最優(yōu)參數(shù)組合[2，1 000]。

圖5 模態(tài)個(gè)數(shù)與綜合指標(biāo)間的關(guān)系

圖6 懲罰因子與綜合指標(biāo)間的關(guān)系

3.2 VMD分解結(jié)果

用上一步優(yōu)化后參數(shù)組合對軸承信號進(jìn)行VMD信號分解，得到2個(gè)IMF分量，如圖7所示。圖8為IMFs頻譜圖，圖9為IMFs包絡(luò)譜圖。

圖7 IMF波形圖

圖8 IMFs頻譜圖

用所構(gòu)建的加權(quán)合成峭度來優(yōu)選IMF，可以綜合考慮軸承故障脈沖信號的沖擊性及循環(huán)平穩(wěn)性。峭度、合成峭度、相關(guān)系數(shù)及加權(quán)合成峭度值如表1所示。

表1 各IMF計(jì)算的相關(guān)值

據(jù)表1 中可以看出，IMF2峭度、合成峭度、相關(guān)系數(shù)及加權(quán)合成峭度值均大于IMF1，故選IMF2為最優(yōu)IMF。所構(gòu)建的加權(quán)合成峭度可用來選取最優(yōu)IMF，且其可以綜合考慮軸承故障脈沖信號的沖擊性及循環(huán)平穩(wěn)性。在圖9(b)中，可以清晰地讀出故障頻率的一倍頻fi、二倍頻2fi、三倍頻3fi、四倍頻4fi至7fi。所提方法具有一定的有效性。

圖9 IMFs包絡(luò)譜圖

4 試驗(yàn)分析

4.1 軸承內(nèi)圈信號試驗(yàn)分析

為了驗(yàn)證本文方法合理性，用美國西儲大學(xué)軸承數(shù)據(jù)[16]進(jìn)行驗(yàn)證。6205-2RS JEM SKF，故障類型為內(nèi)圈故障，電火花線切割所加工深度為0.28 mm、直徑為0.18 mm 輕微故障，電機(jī)負(fù)載為0，對應(yīng)轉(zhuǎn)速取1 797 r/min，采樣頻率取12 kHz，采樣點(diǎn)數(shù)選8 192。經(jīng)計(jì)算得到軸的轉(zhuǎn)頻為29.95 Hz，其故障頻率為162.19 Hz，為了模擬實(shí)際工況條件，加-5 dB的高斯白噪聲。圖10和圖11是相應(yīng)波形及加噪信號，圖12對應(yīng)包絡(luò)譜圖。

圖10 內(nèi)圈信號

圖11 內(nèi)圈加噪信號

加入噪聲后信號信噪比為-5 dB，而均方誤差為0.262 8。故障特征頻率fi于圖12 可以被看到，但是其余故障特征信息均被噪聲所包圍，難以清晰地看到其倍頻。

圖12 內(nèi)圈加入噪聲的包絡(luò)譜圖

(1)VMD參數(shù)優(yōu)化

經(jīng)VMD分解后模態(tài)數(shù)與綜合評價(jià)指標(biāo)f的關(guān)系如圖13 所示，當(dāng)K=2 時(shí)，綜合指標(biāo)f最小，其值為1.146，故取K=2；懲罰因子α與綜合評價(jià)指標(biāo)f關(guān)系如圖14 所示。當(dāng)α=1 000 時(shí)，f取得最小值，則取a=1 000，所以確定最優(yōu)參數(shù)組合[2，1 000]。

圖13 模態(tài)個(gè)數(shù)與綜合指標(biāo)間的關(guān)系

圖14 懲罰因子與綜合指標(biāo)間的關(guān)系

(2)VMD分解結(jié)果

用前一步優(yōu)化后參數(shù)組合對軸承數(shù)據(jù)進(jìn)行VMD 分解，得到2 個(gè)IMF 分量，圖15 是IMFs信號，圖16 為IMF2包絡(luò)譜圖。所計(jì)算的峭度、合成峭度、相關(guān)系數(shù)加權(quán)及合成峭度指標(biāo)的相關(guān)值如表2所示。

表2 各IMF計(jì)算的相關(guān)值

圖15 IMFs信號

圖16 IMF2包絡(luò)譜圖

在表2 中，IMF2加權(quán)合成峭度值大于IMF1，即IMF2為最優(yōu)IMF。在圖16中，可以清晰地讀出故障頻率的一倍頻fi、二倍頻2fi乃至三倍頻3fi及一倍轉(zhuǎn)頻fr、二倍轉(zhuǎn)頻2fr和邊頻帶fi-2fr、fi+2fr及2fi-2fr。表明所提方法對軸承內(nèi)圈故障判別具有一定的有效性。

4.2 軸承外圈信號試驗(yàn)分析

為了進(jìn)一步驗(yàn)證本文方法合理性，采用上述軸承數(shù)據(jù)中的外圈故障數(shù)據(jù)繼續(xù)進(jìn)行驗(yàn)證。其中外圈故障選擇發(fā)生在6 點(diǎn)鐘方向損傷的故障數(shù)據(jù)，故障尺寸：深度為0.28 mm、直徑為0.18 mm輕微故障，電機(jī)負(fù)載為0，轉(zhuǎn)速為1 797 r/min，采樣頻率取12 kHz，采樣點(diǎn)數(shù)取8 192。軸的轉(zhuǎn)頻為29.95 Hz，經(jīng)計(jì)算得到其故障頻率fo為107.36 Hz，同樣加-5 dB的高斯白噪聲。

VMD對上述信號分解，當(dāng)K=2時(shí)，綜合指標(biāo)f最小，其值為1.030，故取K=2；當(dāng)α=1 000 時(shí)，f取得最小值，所以確定最優(yōu)參數(shù)組合[2，1 000]。采用所優(yōu)化參數(shù)的VMD 分解信號得到2 個(gè)IMF，計(jì)算其EKCI值，分別為1.151 6 及1.804 6，故IMF2的值較大，為最佳IMF，圖17是對應(yīng)的包絡(luò)譜圖。

在圖17中，可以清晰地讀出外圈故障頻率的一倍頻fo、二倍頻2fo、三倍頻3fo、四倍頻4fo至七倍頻7fo及一倍轉(zhuǎn)頻fr。所提方法針對外圈故障診斷也具有一定的有效性。

圖17 IMF2包絡(luò)譜圖

為了能充分表明VMD 的優(yōu)越性，采用EMD 對4.1節(jié)的故障數(shù)據(jù)進(jìn)行分析，同樣條件加入相同噪聲且采用加權(quán)合成峭度優(yōu)選IMF。EMD 分解后IMFs中僅前5 個(gè)分量加權(quán)合成峭度值較大，其余均小于0.3，僅給出前5 個(gè)分量，如圖18 所示，計(jì)算所得各IMF分量的EKCI值如表3所示。

圖18 EMD信號分解時(shí)域圖

表3 EMD分解的IMFs加權(quán)合成峭度值

據(jù)表3可知，IMF1的EKCI值最大，為最佳IMF，圖19 為其包絡(luò)譜圖。從圖中僅可以看到其2fr及fi，其余均被噪聲包圍，僅靠這些信息難以判斷軸承故障類型。相較于圖16，VMD方法故障特征信息更為清晰、豐富。故所提方法可對軸承內(nèi)圈故障進(jìn)行有效診斷。

圖19 IMF1的包絡(luò)譜圖

同樣采用EMD對4.2節(jié)軸承外圈故障數(shù)據(jù)進(jìn)行分析，同樣條件加入相同噪聲且采用加權(quán)合成峭度優(yōu)選IMF。經(jīng)計(jì)算所得各IMF 分量的EKCI值僅IMF1最大，值為1.648 4，為最佳IMF，其包絡(luò)譜圖如圖20所示。

在圖20中，可以清晰地讀出外圈故障頻率的一倍頻fo、二倍頻2fo、四倍頻4fo至6fo，但是難以清晰地看到其轉(zhuǎn)頻fr及三倍頻3fo，而且噪聲的幅值較大。相較于圖17，VMD方法故障特征信息更為豐富。所提方法針對外圈故障診斷也具有一定的有效性。

圖20 IMF1的包絡(luò)譜圖

5 結(jié)語

本文構(gòu)建了一種平均包絡(luò)譜峭度結(jié)合平均樣本熵共同優(yōu)化VMD 參數(shù)及加權(quán)合成峭度進(jìn)行故障特征提取的滾動軸承故障診斷方法，運(yùn)用內(nèi)圈仿真信號、西儲大學(xué)軸承內(nèi)、外圈數(shù)據(jù)分別驗(yàn)證了該方法的有效性，實(shí)現(xiàn)了VMD 參數(shù)自適應(yīng)分解，可從中提取出有效IMF，從而進(jìn)行有效的故障判別。

（1）利用平均包絡(luò)譜峭度結(jié)合平均樣本熵來自適應(yīng)優(yōu)化VMD 參數(shù)組合，運(yùn)用優(yōu)化參數(shù)的VMD 分解信號，取得了一定的效果。

（2）VMD分解出K個(gè)IMFs，構(gòu)建了加權(quán)合成峭度綜合考慮故障脈沖的周期性及循環(huán)平穩(wěn)性，進(jìn)而提取最優(yōu)IMF，包絡(luò)譜分析，其包絡(luò)譜圖可明顯讀出故障特征頻率及相應(yīng)倍頻。

（3）通過將VMD 參數(shù)自適應(yīng)優(yōu)化結(jié)合加權(quán)合成峭度提取故障特征的方法，經(jīng)與EMD 方法相比，所提方法具有一定的可行性。