李永琪，彭珍瑞

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

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

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

針對(duì)上述問(wè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,α]，再對(duì)信號(hào)進(jìn)行分解，選擇加權(quán)合成峭度值最大IMF 作為最優(yōu)IMF，包絡(luò)譜分析后判別故障類型。經(jīng)軸承內(nèi)圈故障仿真信號(hào)及實(shí)際軸承內(nèi)、外圈數(shù)據(jù)驗(yàn)證，表明本文方法可有效實(shí)現(xiàn)軸承內(nèi)、外圈故障類型的判斷。

1 基本理論

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

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

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

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

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

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

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

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

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

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

1.3 平均樣本熵

Richman 等[12]提出樣本熵概念。假如涵蓋越多故障信息，出現(xiàn)與故障相關(guān)的周期性沖擊脈沖，信號(hào)的干擾和噪聲較弱，呈現(xiàn)出較強(qiáng)的規(guī)律性，樣本熵值越小表明信號(hào)周期性沖擊越明顯[13]。因此運(yùn)用平均樣本熵值最小來(lái)構(gòu)建目標(biāo)函數(shù)來(lái)總體評(píng)估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)合成峭度

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

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

然而，這些基于峭度的方法往往僅偏向于單個(gè)脈沖頻帶，卻不是所期望的缺陷脈沖，因它們僅僅重視了故障脈沖的沖擊性，而忽略了其循環(huán)平穩(wěn)性。合成峭度[10]可有效解決峭度的這個(gè)缺陷，通過(guò)它們之間期望故障脈沖相比較，發(fā)現(xiàn)僅合成峭度對(duì)故障脈沖識(shí)別率較高，可將其用于脈沖故障識(shí)別。考慮到峭度及相關(guān)系數(shù)各自優(yōu)缺點(diǎn)，且為了綜合考慮脈沖信號(hào)周期的沖擊性和循環(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為信號(hào)x與y相關(guān)系數(shù)；和分別為x及y的均值；E[·]為數(shù)學(xué)期望。

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

基于以上分析，采用平均包絡(luò)譜峭度及平均樣本熵可以總體評(píng)價(jià)其故障成分來(lái)優(yōu)化VMD參數(shù)，利用加權(quán)合成峭度可以綜合考慮故障脈沖周期性及循環(huán)平穩(wěn)性來(lái)優(yōu)選IMF并提取故障特征。同時(shí)為避免在優(yōu)化中某一指標(biāo)信息被淹沒(méi)的情況，故將平均包絡(luò)譜峭度、平均樣本熵分別進(jìn)行歸一化處理并求和，即可得到目標(biāo)函數(shù)式，如式(12)所示，從而構(gòu)建綜合評(píng)價(jià)指標(biāo)f。本文以綜合評(píng)價(jià)指標(biāo)最小為原則，表明各項(xiàng)指標(biāo)都較好，由節(jié)1.2 及1.3 可知平均樣本熵以最小來(lái)構(gòu)建目標(biāo)函數(shù)來(lái)總體評(píng)估IMF 分量，而平均包絡(luò)譜峭度以最大來(lái)構(gòu)建目標(biāo)函數(shù)來(lái)總體評(píng)估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]提出的按信息熵最小來(lái)逐個(gè)優(yōu)化VMD 參數(shù)方法，先給定α=2 000，再利用f最小值來(lái)確定K；然后利用選定好的K，接著再對(duì)α進(jìn)行優(yōu)化。整體流程圖如圖1所示。具體過(guò)程如下：

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

(1)采集原始信號(hào)；

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

(3) 接著對(duì)懲罰因子α再進(jìn)行優(yōu)化。由之前優(yōu)化的K再對(duì)α進(jìn)行優(yōu)化，初始化α=1 000，搜索范圍[1 000，3 000]，步長(zhǎng)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ò)譜對(duì)其進(jìn)行分析，提取故障特征并輸出故障類型。

3 仿真分析

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

式中：周期沖擊信號(hào)s(t)；高斯白噪聲n(t)，對(duì)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。加入噪聲后信號(hào)信噪比為-10 dB，其均方誤差為0.300 7。圖2 是s(t)沖擊信號(hào)，圖3 是仿真信號(hào)x(t)，圖4 是圖3 的包絡(luò)譜圖。盡管可以從仿真信號(hào)包絡(luò)譜圖中看到故障特征頻率，但其余周邊均有噪聲圍繞。

圖2 沖擊信號(hào)

圖3 仿真信號(hào)

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

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

對(duì)上述信號(hào)進(jìn)行VMD分解，其模態(tài)數(shù)K與綜合評(píng)價(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ù)組合對(duì)軸承信號(hào)進(jìn)行VMD信號(hào)分解，得到2個(gè)IMF分量，如圖7所示。圖8為IMFs頻譜圖，圖9為IMFs包絡(luò)譜圖。

圖7 IMF波形圖

圖8 IMFs頻譜圖

用所構(gòu)建的加權(quán)合成峭度來(lái)優(yōu)選IMF，可以綜合考慮軸承故障脈沖信號(hào)的沖擊性及循環(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)合成峭度可用來(lái)選取最優(yōu)IMF，且其可以綜合考慮軸承故障脈沖信號(hào)的沖擊性及循環(huán)平穩(wěn)性。在圖9(b)中，可以清晰地讀出故障頻率的一倍頻fi、二倍頻2fi、三倍頻3fi、四倍頻4fi至7fi。所提方法具有一定的有效性。

圖9 IMFs包絡(luò)譜圖

4 試驗(yàn)分析

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

為了驗(yàn)證本文方法合理性，用美國(guó)西儲(chǔ)大學(xué)軸承數(shù)據(jù)[16]進(jìn)行驗(yàn)證。6205-2RS JEM SKF，故障類型為內(nèi)圈故障，電火花線切割所加工深度為0.28 mm、直徑為0.18 mm 輕微故障，電機(jī)負(fù)載為0，對(duì)應(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)波形及加噪信號(hào)，圖12對(duì)應(yīng)包絡(luò)譜圖。

圖10 內(nèi)圈信號(hào)

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

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

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

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

經(jīng)VMD分解后模態(tài)數(shù)與綜合評(píng)價(jià)指標(biāo)f的關(guān)系如圖13 所示，當(dāng)K=2 時(shí)，綜合指標(biāo)f最小，其值為1.146，故取K=2；懲罰因子α與綜合評(píng)價(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ù)組合對(duì)軸承數(shù)據(jù)進(jìn)行VMD 分解，得到2 個(gè)IMF 分量，圖15 是IMFs信號(hào)，圖16 為IMF2包絡(luò)譜圖。所計(jì)算的峭度、合成峭度、相關(guān)系數(shù)加權(quán)及合成峭度指標(biāo)的相關(guān)值如表2所示。

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

圖15 IMFs信號(hào)

圖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。表明所提方法對(duì)軸承內(nèi)圈故障判別具有一定的有效性。

4.2 軸承外圈信號(hào)試驗(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對(duì)上述信號(hào)分解，當(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 分解信號(hào)得到2 個(gè)IMF，計(jì)算其EKCI值，分別為1.151 6 及1.804 6，故IMF2的值較大，為最佳IMF，圖17是對(duì)應(yīng)的包絡(luò)譜圖。

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

圖17 IMF2包絡(luò)譜圖

為了能充分表明VMD 的優(yōu)越性，采用EMD 對(duì)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信號(hào)分解時(shí)域圖

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

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

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

同樣采用EMD對(duì)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方法故障特征信息更為豐富。所提方法針對(duì)外圈故障診斷也具有一定的有效性。

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

5 結(jié)語(yǔ)

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

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

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

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