彭康健，陳君若，吳智恒

（650500 云南省 昆明市 昆明理工大學(xué) 機(jī)電工程學(xué)院）

0 引言

滾動(dòng)軸承作為機(jī)械設(shè)備中的重要零部件，對(duì)于機(jī)械設(shè)備的運(yùn)行起著至關(guān)重要的作用，據(jù)統(tǒng)計(jì)，在使用滾動(dòng)軸承的機(jī)械設(shè)備中，約有30%的機(jī)械故障是由軸承故障而引起的[1]。作為機(jī)械設(shè)備中的重要零部件，一旦發(fā)生故障，不僅會(huì)影響其正常運(yùn)行，還可能給企業(yè)帶來(lái)巨大的經(jīng)濟(jì)上的損失甚至安全問(wèn)題[2-3]。因此滾動(dòng)軸承的狀態(tài)監(jiān)測(cè)和故障診斷具有重要的現(xiàn)實(shí)研究意義。

軸承故障診斷過(guò)程大致可分為信號(hào)采集、特征提取、故障識(shí)別和診斷分析4 個(gè)步驟，其中特征提取可以說(shuō)是故障診斷過(guò)程中最關(guān)鍵的一步。常用滾動(dòng)軸承故障特征提取的方法是，先通過(guò)傳感器采集軸承的振動(dòng)信號(hào)，這些振動(dòng)信號(hào)與設(shè)備的運(yùn)行狀態(tài)直接相關(guān)，但由于采集的信號(hào)通常是非線性、非平穩(wěn)的，所以并不能直接提取故障特征[4]。于是一些信號(hào)時(shí)頻分析方法應(yīng)運(yùn)而生，并對(duì)設(shè)備的監(jiān)測(cè)監(jiān)控有著重大現(xiàn)實(shí)意義[5]。目前常用的時(shí)頻分析方法有：Huang[6]等人提出的經(jīng)驗(yàn)?zāi)B(tài)分解（Empircial Mode Decomposition）方法；Smith[7]提出的一種針對(duì)非線性非平穩(wěn)信號(hào)的局部均值分解（Local Mean Decomposition）方法；以及Wu[8]等又在現(xiàn)有經(jīng)驗(yàn)?zāi)B(tài)分解方法基礎(chǔ)上提出的集合經(jīng)驗(yàn)?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition）方法。上述方法都在軸承故障診斷分析中取得了良好的效果，但也有著各自的局限性與不足。比如EMD 方法缺少?lài)?yán)謹(jǐn)?shù)臄?shù)學(xué)基礎(chǔ)，容易產(chǎn)生模態(tài)混疊和端點(diǎn)效應(yīng)；LMD 方法計(jì)算量大，可能產(chǎn)生由解調(diào)引起的信號(hào)突變問(wèn)題。而變分模態(tài)分解方法作為一種相對(duì)成熟的算法，目前在軸承故障診斷領(lǐng)域有著不錯(cuò)的效果。

變分模態(tài)分解方法作為一種自適應(yīng)信號(hào)處理方法有著較為成熟的理論基礎(chǔ)，于2014 年由Dragomiretskiy 等人提出[9]。該方法通過(guò)迭代搜索變分模型的最優(yōu)解來(lái)確定各個(gè)分量的頻率中心和帶寬，從而有效地獲取模態(tài)分量，自適應(yīng)地完成軸承信號(hào)的頻域及各分量的有效分離。相對(duì)于傳統(tǒng)EMD 方法，該方法有著良好的數(shù)學(xué)理論支撐，同時(shí)有著較好的抗噪性和自適應(yīng)性，避免了遞歸模式分解造成的包絡(luò)線估計(jì)誤差不斷累積，克服了端點(diǎn)效應(yīng)[10]。

本文針對(duì)變分模態(tài)分解中的模態(tài)數(shù)K 和懲罰因子α兩個(gè)參數(shù)的取值，通過(guò)熵最小值的方法確定最佳參數(shù)組合[K，α]。然后對(duì)信號(hào)進(jìn)行變分模態(tài)分解，得到了K 個(gè)IMF 分量，選取包絡(luò)熵最小的IMF 分量，對(duì)其IMF 分量進(jìn)行包絡(luò)解調(diào)分析，從而提取軸承故障特征頻率信息。

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

變分模態(tài)分解（Variational Mode Decomposition，VMD）是基于維納濾波的自適應(yīng)信號(hào)處理方法，它把信號(hào)分解成一系列固有模態(tài)分量（Intrinsic Modal Function，IMF），每個(gè)IMF 分量表示為

式中：Ak（t）——uk（t）的瞬時(shí)幅值；ωk（t）——uk（t）的瞬時(shí)頻率，ωk（t）=φk'（t）=dφk/dt。

對(duì)于方法變分模型，假設(shè)信號(hào)被分解成K 個(gè)模態(tài)分量（IMF），每個(gè)模態(tài)分量uk(t)都緊密?chē)@在一個(gè)確定的中心頻率ωk（t）周?chē)?，?duì)于一個(gè)一維信號(hào)f，確定其模態(tài)分量uk（t）的頻帶范圍的算法過(guò)程如下：

（1）通過(guò)希爾伯特變換方法得到與各個(gè)模態(tài)分量uk（t）相應(yīng)的解析信號(hào)，從而獲得其單邊頻譜。

（2）加入指數(shù)項(xiàng)調(diào)整各自估計(jì)的中心頻率，將每個(gè)模態(tài)分量的頻譜調(diào)制到基頻帶上。

（3）通過(guò)H1高斯平滑估計(jì)解調(diào)信號(hào)的帶寬，即計(jì)算其梯度二范數(shù)的平方。變分問(wèn)題就變成了尋找K 個(gè)模態(tài)分量uk（t），并使所有模態(tài)分量估計(jì)帶寬最小的問(wèn)題。設(shè)原始信號(hào)經(jīng)過(guò)變分模態(tài)分解后的分量數(shù)為K，則約束變分模型為

式中：{uk}={u1，u1，…，uk}——分解后的k 個(gè)模態(tài)分量；{ωk}={ω1，ω1，…，ωk}——各模態(tài)分量的頻率中心；δ——狄拉克分布；t——時(shí)間；*——卷積；k 越大表示的頻率成分越低。

為確保信號(hào)的絕對(duì)可積性，這里引入一個(gè)二次懲罰因子和增廣 Lagrange 函數(shù)，即

式中：?——懲罰因子；λ——Lagrange 乘子。

根據(jù)上述理論，VMD 計(jì)算流程如下：

Step2：n=n+1,執(zhí)行循環(huán)程序；

Step3：更新{uk}，{ωk}；

Step4：k=k+1，重復(fù)step3，直到k=K；

Step5：更新λn+1=λn+τ（f-∑ukn+1）；

Step6：設(shè)置判別精度ε，直到滿足迭代停止條件，結(jié)束循環(huán)，得到K 個(gè)模態(tài)分量輸出。

2 鯨魚(yú)優(yōu)化算法

鯨魚(yú)優(yōu)化算法[11]是一種啟發(fā)式算法，具有簡(jiǎn)單易實(shí)現(xiàn)、能有效避免局部最優(yōu)、可擴(kuò)展性好等優(yōu)點(diǎn)。目前對(duì)于變分模態(tài)分解參數(shù)優(yōu)化已經(jīng)使用的優(yōu)化算法有：遺傳算法[12]、果蠅優(yōu)化算法、蝙蝠優(yōu)化算法、粒子群優(yōu)化算法[13]等。而鯨魚(yú)優(yōu)化算法相比于以上幾種算法有原理簡(jiǎn)單、易于實(shí)現(xiàn)、調(diào)整參數(shù)少以及魯棒性強(qiáng)等優(yōu)點(diǎn)。

鯨魚(yú)算法的原理主要模擬座頭鯨的捕食行為，這種較為獨(dú)特的捕食行為被叫作泡泡網(wǎng)覓食法，如圖1 所示。

圖1 鯨魚(yú)捕食圖Fig.1 Whale hunting map

鯨魚(yú)優(yōu)化算法的過(guò)程主要包括隨機(jī)搜尋食物、包圍捕食和氣泡捕食3 個(gè)階段。其捕食過(guò)程的數(shù)學(xué)模型如下：

2.1 鯨魚(yú)的包圍捕食

鯨魚(yú)捕食過(guò)程中，先觀察和尋找獵物的位置，然后對(duì)獵物進(jìn)行包圍。假設(shè)待優(yōu)化的問(wèn)題解及問(wèn)題變量為最優(yōu)的鯨魚(yú)所在的位置，其他鯨魚(yú)也會(huì)朝向該鯨魚(yú)的位置移動(dòng)從而更新自身的位置，所以要求的是個(gè)體與最優(yōu)鯨魚(yú)位置（獵物）之間的距離

鯨魚(yú)根據(jù)最優(yōu)鯨魚(yú)的位置更新而更新

2.2 鯨魚(yú)的氣泡捕食階段

鯨魚(yú)的氣泡捕食階段分為兩種方式：

（1）搖擺包圍捕食

（2）螺旋吐氣泡捕食

此階段中，鯨魚(yú)首先計(jì)算出自身到獵物（當(dāng)前為止最好的位置）之間的距離，然后鯨魚(yú)會(huì)在以螺旋形姿勢(shì)向上游動(dòng)的過(guò)程中吐出大小不等的氣泡捕食魚(yú)蝦。該行為的數(shù)學(xué)模型如下：

式中：p——[0，1]之間的隨機(jī)數(shù)。鯨魚(yú)的捕食行為除了氣泡捕食方法外，還會(huì)隨機(jī)尋找視野范圍內(nèi)的獵物，該行為為鯨魚(yú)的全局搜索。

然而，鯨魚(yú)在捕食過(guò)程中會(huì)隨著其他鯨魚(yú)位置改變而更新自身位置，也就是說(shuō)此階段的鯨魚(yú)不再隨著最好的鯨魚(yú)的位置更新自身位置，而是大范圍地隨機(jī)搜索獵物，因此算法是隨著收斂因子A 的變化而進(jìn)行更大范圍的搜索過(guò)程。當(dāng) A>1時(shí)，鯨魚(yú)將進(jìn)行隨機(jī)搜索獵物的行為，這樣可以避免陷入局部最優(yōu)的情況，其數(shù)學(xué)表達(dá)式如下：

（1）初始化參數(shù)：鯨魚(yú)的種群規(guī)模X、最大迭代次數(shù)Tmax，初始化鯨魚(yú)種群位置；

（2）計(jì)算各鯨魚(yú)的適應(yīng)度值，并確定當(dāng)前最優(yōu)鯨魚(yú)個(gè)體X*；

（3）進(jìn)入算法主循環(huán)，如果p<0.5 且|A1|<1，每只鯨魚(yú)個(gè)體按照式（6）更新當(dāng)前位置，否則按照式（12）更新鯨魚(yú)個(gè)體位置。如果p ≥5.0，則每只鯨魚(yú)個(gè)體依據(jù)式（9）更新位置；

（4）再次對(duì)整個(gè)鯨魚(yú)種群進(jìn)行評(píng)價(jià)，找到全局最優(yōu)的鯨魚(yú)個(gè)體及其位置；

（5）若滿足算法的終止條件(最大迭代次數(shù))，則結(jié)束；否則轉(zhuǎn)到步驟（2），繼續(xù)進(jìn)行算法迭代。

（6）輸出全局最優(yōu)解*X 。

3 基于VMD 的軸承故障特征提取方法

VMD 分解過(guò)程中的模態(tài)分量K 和懲罰因子α的選取直接關(guān)系著VMD 的分解效果。如果分量個(gè)數(shù)K 取值過(guò)小，可能造成模態(tài)混疊的現(xiàn)象；K 值過(guò)大，分解的結(jié)果中就可能產(chǎn)生虛假分量，不利于分析。而二次懲罰因子α的取值越小，則分解后的各IMF 分量的帶寬越大；反之，α越大，則各個(gè)IMF 分量的帶寬越小。所以，選擇適當(dāng)?shù)腒、α值對(duì)VMD 分解效果非常重要。

對(duì)于優(yōu)化問(wèn)題，首先需要設(shè)定一個(gè)適應(yīng)度函數(shù)，以此函數(shù)來(lái)確定最佳參數(shù)。故使用鯨魚(yú)優(yōu)化算法對(duì)VMD 進(jìn)行參數(shù)優(yōu)化時(shí)，需設(shè)定一個(gè)適應(yīng)度函數(shù)來(lái)進(jìn)行迭代尋優(yōu)，確定VMD 的最佳參數(shù)組合。文獻(xiàn)[3]提出了包絡(luò)熵的概念，包絡(luò)熵反映了信號(hào)的稀疏特性，若信號(hào)經(jīng)過(guò)VMD 分解后IMF 分量中含有較多的周期性故障沖擊信息時(shí)，則其稀疏特性強(qiáng)，包絡(luò)熵值小；反之，如果分解后IMF 分量中只有較少的周期性故障沖擊信息時(shí)，則其稀疏特性弱，包絡(luò)熵值大。本文中將VMD 分解后的局部極小包絡(luò)熵作為鯨魚(yú)優(yōu)化算法中的適應(yīng)度函數(shù)，搜尋最優(yōu)的VMD 參數(shù)組合[k，α],即模態(tài)分量個(gè)數(shù)和二次懲罰因子。包絡(luò)熵Ep的計(jì)算公式為

式中：pi——a（i）的歸一化形式；a（i）——信號(hào)x（i）經(jīng)Hilbert 解調(diào)后得到的包絡(luò)信號(hào)。

鯨魚(yú)算法優(yōu)化VMD 參數(shù)步驟：

（1）確定適應(yīng)度函數(shù)，設(shè)鯨魚(yú)個(gè)體的解為（α，K），初始化算法參數(shù)；

（2）將信號(hào)進(jìn)行VMD 過(guò)程得到各個(gè) IMF，根據(jù)式子計(jì)算每個(gè) IMF 的包絡(luò)熵；

（3）將最小包絡(luò)熵值作為適應(yīng)度函數(shù)進(jìn)行全局搜索；

（4）根據(jù)式（6）、式（12）更新鯨魚(yú)個(gè)體的位置；

（5）重復(fù)步驟（3）至步驟（5），只要達(dá)到包絡(luò)熵值最小或者達(dá)到設(shè)定的最大迭代次數(shù)，則輸出最佳參數(shù)組合個(gè)體（α，K）；

（6）設(shè)置最優(yōu)參數(shù)（α，K）組合對(duì)信號(hào)進(jìn)行 VMD 分解。

表1 中：X——鯨魚(yú)的種群規(guī)模；Tmax——最大迭代次數(shù)；Lb——分解層數(shù)K 和懲罰系數(shù)alpha的下限；Ub——分解層數(shù)K 和懲罰系數(shù)alpha 的下限和上限。判別精度設(shè)置為1×10-6。

表1 鯨魚(yú)算法各參數(shù)設(shè)置Tab.1 Each parameter setting of whale algorithm

如圖2 所示，對(duì)于軸承信號(hào)處理，首先是獲取原始信號(hào)，對(duì)原始信號(hào)進(jìn)行鯨魚(yú)算法優(yōu)化確定VMD 相關(guān)參數(shù)組合（α，k），然后將信號(hào)進(jìn)行VMD 分解后得到各IMF 分量，并選擇最佳模態(tài)分量進(jìn)行包絡(luò)譜分析，得到故障特征頻率，從而判別故障類(lèi)型，實(shí)現(xiàn)故障診斷。

圖2 軸承故障診斷步驟Fig.2 Step of bearing fault diagnosis

4 實(shí)驗(yàn)驗(yàn)證

采用美國(guó)西儲(chǔ)大學(xué)滾動(dòng)軸承故障模擬實(shí)驗(yàn)臺(tái)的公開(kāi)數(shù)據(jù)進(jìn)行分析驗(yàn)證。實(shí)驗(yàn)軸承安裝在電機(jī)轉(zhuǎn)軸兩側(cè)，其中一端為驅(qū)動(dòng)端，另一端為風(fēng)扇端。這里選用的待測(cè)軸承安裝在轉(zhuǎn)速為 1 797 r/min 驅(qū)動(dòng)端軸承SKF6205，采樣頻率為12 kHz；

表2 6205-2RS JEM SKF 深溝球軸承Tab.2 6205-2RS JEM SKF deep groove bearing

其中：2RS 表示雙面密封，SKF（Svenska Kullager-Fabriken）為斯凱孚公司，軸承滾動(dòng)體個(gè)數(shù)為9 個(gè)，軸承的故障類(lèi)型為電火花加工的單點(diǎn)損傷，損傷直徑0.177 8 mm。

根據(jù)公式計(jì)算得軸承理論的外圈、內(nèi)圈、故障特征頻率分別為：BPO=107.3 Hz、BPI=162.1 Hz，實(shí)測(cè)信號(hào)的時(shí)域波形圖如圖3、圖4 所示。

圖3 滾動(dòng)軸承內(nèi)圈時(shí)域波形圖Fig.3 Time domain waveform of inner ring of rolling bearing

圖4 滾動(dòng)軸承外圈時(shí)域波形圖Fig.4 Time domain waveform of moving bearing outer ring

將上述軸承信號(hào)加載至MATLAB2020a 中，以最小包絡(luò)熵為適應(yīng)度函數(shù)，使用鯨魚(yú)算法優(yōu)化變分模態(tài)的α，K 值。由圖5、圖6 可以得到：對(duì)于內(nèi)圈故障，當(dāng)群體進(jìn)化代數(shù)到第11 代時(shí)得到了局部極小熵值4.887；對(duì)于外圈故障，當(dāng)群體進(jìn)化代數(shù)到5 代時(shí)得到了局部極小熵值4.914。于是可分別得到最佳α，K 組合：內(nèi)圈為 alpha= 1 880，K=4；外圈為alpha=2 596，K=6。

圖5 鯨魚(yú)算法優(yōu)化值隨進(jìn)化代數(shù)的變化曲線（內(nèi)圈故障）Fig.5 Curve of optimal value of whale algorithm with evolutionary algebra (inner)

圖6 鯨魚(yú)算法優(yōu)化值隨進(jìn)化代數(shù)的變化曲線（外圈故障）Fig.6 Curve of optimal value of whale algorithm with evolutionary algebra (outer)

確定參數(shù)后，分別對(duì)兩種故障的振動(dòng)信號(hào)進(jìn)行變分模態(tài)分解和頻譜、包絡(luò)譜分析，分別如圖7—圖9 所示，圖10 為最佳模態(tài)分量包絡(luò)譜。

圖7 模態(tài)分量IMF1-IMF4Fig.7 Modal component IMF1-IMF4

圖8 模態(tài)頻譜Fig.8 Modal frequency spectrum

圖9 IMF1-IMF4 分量包絡(luò)譜Fig.9 IMF1-IMF4 envelope spectrum

圖10 最佳模態(tài)分量IMF4 包絡(luò)譜Fig.10 The best component IMF4 envelope spectrum

外圈故障VMD分解各模態(tài)分量如圖11所示，圖12 為最佳模態(tài)分量。

圖11 IMF1-IMF4 分量包絡(luò)譜Fig.11 IMF1-IMF4 envelope spectrum

圖12 最佳IMF 分量Fig.12 Best IMF

由頻譜圖可見(jiàn)，軸承轉(zhuǎn)頻為60 Hz，各種類(lèi)型的故障頻率161.9，107.7 Hz分別對(duì)應(yīng)軸承內(nèi)圈、外圈的理論故障特征頻率162.1，107.3 Hz。

同樣，對(duì)內(nèi)外圈軸承振動(dòng)信號(hào)進(jìn)行EMD 分解，對(duì)于同樣的故障振動(dòng)信號(hào)，對(duì)于內(nèi)圈故障，當(dāng)EMD 分解到IMF5 時(shí)，就基本看不到故障特征頻率了；對(duì)于外圈故障，當(dāng)EMD 分解到IMF6 時(shí)，就基本看不到故障特征頻率了。在這里我們?nèi)∮行У腎MF 分量進(jìn)行頻譜分析，從圖10 和圖13、圖12 和圖14 中IMF 分量進(jìn)行對(duì)比可以看出，VMD 分解最佳分量相對(duì)EMD 分解分量幅值更大，故障特征頻率相對(duì)更加明顯。

圖13 內(nèi)圈EMD 分解IMF1-IMF4Fig.13 Inner EMD IMF1-IMF4

圖14 外圈EMD 分解IMF1-IMF5Fig.14 Outer EMD IMF1-IMF5

5 結(jié)論

VMD 方法作為一種有著良好數(shù)學(xué)理論支撐的自適應(yīng)分解算法，能夠有效改善信號(hào)分解過(guò)程中的模態(tài)混疊與端點(diǎn)效應(yīng)問(wèn)題，對(duì)非平穩(wěn)非線性的信號(hào)的故障特征頻率的提取有著良好的處理效果。本文針對(duì)VMD 算法中參數(shù) K與α的取值問(wèn)題，提出了基于最小包絡(luò)熵的鯨魚(yú)算法優(yōu)化方法，并證明該方法能夠較為有效地得到參數(shù) K 與α的取值并得到較好分解效果。通過(guò)與EMD 方法進(jìn)行實(shí)驗(yàn)對(duì)比，表明該方法相比EMD 方法更能突出故障特征頻率，從而判別故障類(lèi)型。

但是該方法也存在以下不足：（1）算法參數(shù)選取經(jīng)驗(yàn)性較強(qiáng)，在優(yōu)化過(guò)程中可能存在優(yōu)化結(jié)果不穩(wěn)定，當(dāng)?shù)螖?shù)過(guò)大時(shí)也存在運(yùn)行時(shí)間長(zhǎng)的問(wèn)題；（2）除K 和α值外，對(duì)于VMD 的其他參數(shù)的設(shè)定也有待研究。