基于蝙蝠算法優化VMD參數的滾動軸承復合故障分離方法

張 偉， 李軍霞， 陳維望

(1. 太原理工大學 機械與運載工程學院，太原 030024；2. 太原理工大學 礦山流體控制國家地方聯合工程實驗室，太原 030024)

滾動軸承是現代機械設備的重要組成部分，在工業生產系統中發揮著重要作用，廣泛用于礦山、交通運輸和醫藥制造等行業。由于負荷大、運行環境惡劣，滾動軸承等關鍵部件可能發生復合故障，與任何單一故障相比，其危害性更大。復合故障的特征提取一直以來都是故障診斷領域面臨的難題[1]。因此，提供一種能夠從原始振動信號中提取微弱周期特征信息，準確識別復合故障位置和類型的方法，對機械設備安全高效運行具有重要意義。

在工業生產中，機械振動響應是多個頻率特征信息的疊加，需要通過信號分解和濾波等操作來提取故障特征，以便最終進行故障識別[2]。針對此問題，國外學者已經設計了多種方法，如小波分解(wavelet transform，WT)、小波包分解(wavelet packet decomposition, WPT)、經驗模態分解(empirical mode decomposition，EMD)和局部平均分解(local mean decomposition，LMD)等。小波分解和小波包分解是非自適應信號分析方法，需要提前確定小波基函數。EMD可以將復雜信號自適應分解為有限個本征模函數(intrinsic mode function，IMF)，但由于模態混疊和邊界效應等缺陷限制了其應用范圍。為了進一步提高EMD的性能，相繼提出了如集成經驗模態分解(ensemble empirical mode decomposition, EEMD)[3]、互補集合經驗模態分解(complementary ensemble empirical mode decomposition，CEEMD)和自適應噪聲完備集合經驗模態分解(complete ensemble empirical mode decomposition with adaptive noise,CEEMDAN)等方法，并用于機械故障診斷。這些方法在一定程度上緩解了模式混合等問題，但增加的白噪聲不能有效消除，降低了計算效率。

受EMD的啟發，Dragomiretskiy等[4]在2014年提出了一種具有堅實理論基礎的完全非遞歸方法，稱為變分模態分解(variational modal decomposition, VMD)。VMD在信號處理方面具有優異的性能，在音調檢測、音調分離和噪聲魯棒性方面優于EMD，避免了模態混合、端點效應和噪聲抑制等不利影響。此外，VMD已廣泛應用于工程實踐中。胡愛軍等[5]采用參數優化VMD結合1.5維譜的方式，實現了滾動軸承復合故障的有效分離。萬書亭等[6]提出基于VMD及最大相關峭度解卷積(maximum correlated kurtosis deconvolution，MCKD)的復合故障分離方法，仿真及實測信號表明該方法能夠從復合故障中分離出單一故障特征。以上研究表明，采用VMD方法能夠實現滾動軸承復合故障的診斷，但強噪聲干擾下的復合故障診斷仍需要進一步研究。

然而，VMD的優越性能嚴重依賴模態數量和懲罰因子等參數的選取，這些參數對分解效果有著顯著影響。利用經驗或先驗準則選取參數極易導致分解結果不準確，降低了VMD的分解效率。鄭圓等[7]根據峭度最大原則確定分解層數K，再通過鯨魚算法優化選擇懲罰因子α，但忽略了兩個參數之間的關聯性。一些學者利用灰狼優化算法(grey wolf optimizer，GWO)、蚱蜢優化算法(grasshopper optimization algorithm，GOA)、布谷鳥搜索 (cuckoo search，CS)算法等智能優化算法同時確定模態數量和懲罰因子的最優值，與以上仿生算法相比，蝙蝠算法(bat algorithm, BA)采用頻率調諧手段，增加了搜索到最優參數的概率，能夠在全局尋優和局部尋優過程之間自由切換，避免陷入局部收斂，提高了算法的收斂性。Meng等[8]通過多個基準問題和工程設計表明，蝙蝠算法在全局優化方面具有高效性和穩定性。

在采用仿生算法優化VMD時，選取合適的目標函數至關重要，目標函數將直接影響最終分解的準確性和效率。目前常用的目標函數有局部極小熵值[9]、集合峭度、加權峭度[10]、平均包絡熵[11]等，采用這些目標函數取得了不錯的分解效果，但仍存在信息丟失或模式冗余等問題。針對此問題，本文提出一種新的復合影響指數，該指數能夠準確識別故障并抑制噪聲干擾。

為了自適應獲取變分模態分解的參數，本文將最小平均復合影響指數作為目標函數，采用蝙蝠算法對參數空間進行全局尋優。對分解后的信號進行包絡解調分析，提取故障特征。通過仿真信號和工程數據表明，該方法能夠準確識別強噪聲干擾下的滾動軸承復合故障類型。

1 基本原理

1.1 VMD算法

VMD算法通過迭代搜尋變分模型的最優解，將原始信號分解為帶寬之和最小的K個IMF分量uk(t)。

其變分模型為

(1)

式中：f(t)為初始輸入信號； {uk}，{ωk}為各模態信號以及中心頻率；*為卷積運算。

引入二次懲罰因子α和Lagrange算子式λ(t)，將約束性變分問題變為非約束變分問題，其拉格朗日表達式[12]為

(2)

在變分模態分解中，對中心頻率分別進行零值、均勻分布、隨機分布初始化，采用交替方向乘子法進行多次迭代求取式(2)的最優解，不斷更新其中心頻率，結果發現得到的中心頻率基本保持不變。在誤差允許范圍內，為了減少計算量，初始模態設為0矩陣，VMD算法流程如下：

步驟2進行迭代n=n+1，執行整個循環；

步驟3初始k=0，k=k+1，對所有的ω≥0，更新uk，ωk為

(3)

(4)

步驟4更新λ為

(5)

1.2 蝙蝠算法

蝙蝠算法通過模擬蝙蝠的覓食行為對優化問題的解空間進行搜索，通過自適應調節聲波脈沖響度、頻率，以達到尋優的目的[13]。在蝙蝠算法中，每一個蝙蝠代表一個可行解，而蝙蝠的獵物則代表優化問題的最優解。在整個解空間里，蝙蝠頻率、位置和速度的更新計算式分別為

fi=fmin+(fmax-fmin)β

(6)

(7)

(8)

式中：β∈[0,1]的隨機數；x*為t時刻全局搜索過程中蝙蝠的最優位置。

局部搜索時，從當前最優解集中隨機選取一個解并對其施加一個隨機擾動，然后在該解的鄰域內進行搜索。此過程中，蝙蝠位置更新計算式為

xnew=xold+ηAt

(9)

式中：xold為所選取的最優解；η∈[-1,1]的隨機數；At為所有蝙蝠在時刻為t時發出聲波的平均響度。

覓食過程中，為了平衡全局搜索和局部搜索過程，蝙蝠發射的聲波響度和頻率需要跟隨搜索獵物的進程自動調節。蝙蝠發出聲波的響度Ai和發射速率ri的調節計算式分別為

(10)

(11)

式中：α為聲波響度增加系數；γ為脈沖頻度增強因子。

BA算法的偽代碼如下：

目標方程f(x),x=(xi,…,xd)T

初始化蝙蝠種群位置xi和速度vi(i=1,2,…,n)

定義位置xi的發射頻率fi

初始化脈沖發射率ri和響度Ai

while (t<最大迭代次數)

通過式(6)～式(8)更新蝙蝠的頻率、速度和位置形成新的解

if(rand>ri)

選出全局最優解

圍繞選定的最佳解決方案生成局部解決方案

end if

通過隨機飛行產生一個新的解(見式(9))

if(rand

接受新的解

增大ri并且減小Ai(式(10)和式(11))

end if

對所有蝙蝠排序，找到目前最好的解x*

end while

后處理結果和可視化

2 方法步驟

2.1 復合影響指數

選取合適的評估標準對機械故障診斷非常重要，峭度、包絡譜、功率譜、相關系數，熵等單一指標已被廣泛研究。考慮到各個單一指標關注重點不同，目前已經提出了包絡譜峭度(kurtosis of envelope spectrum, KES)、平方包絡譜峭度(kurtosis of squared envelope spectrum, KSES)和包絡熵值(envelope entropy, EE)等混合指數。然而，在工程實踐中，除了背景噪聲外，采集到的振動信號還包含有隨機沖擊和周期性諧波等干擾信號，這對目前的測量指標是一個相當大的挑戰。

本文基于峭度、功率譜、相關系數等指標對故障診斷的影響，提出了一種新的復合影響指數(composite impact index, CII)，以便及時準確檢測到故障沖擊，其計算過程如下：

對于采集到的信號f(n)，其希爾伯特變換后的信號fH(n)計算為

fH(n)=Hilbert{f(n)}

(12)

對希爾伯特變換后的信號fH(n)進行N次采樣，并進行離散傅里葉變換，進而得到信號的功率譜SH(k)

(13)

(k=0,1,…,N-1)

(14)

計算功率譜峭度(kurtosis of power spectrum, KPS)

(15)

考慮到模式與主信號之間的聯系，引入相關系數，其定義為

(16)

式中：f和uk分別為原始信號和分解后的模態分量； E[·]為數學期望。

由此可得，CII的表達式為

CII=|CC|·KPS

(17)

本文提出的指數，與文獻[14]提出的綜合影響指數相比，該參數計算過程相對簡單，省去了構造解析函數和計算包絡譜等過程，加快了計算速度。該指數通過計算信號的功率譜峭度，將脈沖故障特征放大，有利于檢測到信號的故障信息。考慮到相關系數的優點，該指數能夠實現濾波功能，減少噪聲干擾。

為了評估該指標的性能，根據第3章使用的仿真信號，構建了一個包含軸承故障沖擊、隨機沖擊、周期諧波和高斯白噪聲的模擬信號。本節將CII與一些現有指標進行比較，如包絡譜峭度、平方包絡譜峭度、整體峭度(ensemble kurtosis, EK)和包絡熵值，為了公平的比較這些指標的性能，對結果進行歸一化處理。首先，每個指標的振幅除以所有振幅的總和，使它們在[0,1]內；然后，將第一步結果的平方根作為最終的歸一化振幅，結果如圖1所示。

為了擴大復合影響指數的應用范圍，以下證明該指數對其余周期性沖擊依然有效。現構建一個包含周期沖擊、隨機沖擊、周期諧波和高斯白噪聲的模擬信號，分別計算KES、KSES、EK和EE等指標，并進行歸一化處理，結果如圖2所示。

2.2 自適應變分模態分解

對復雜環境下采集的故障信號進行VMD分解時，需要確定模態數K、懲罰因子α、噪聲容限τ和收斂誤差ε等參數。噪聲容忍度和收斂誤差對分解結果的影響較小，一般采用默認值。模態數K和懲罰因子α需要通過蝙蝠算法對全局進行有效的搜索。在選取目標函數時，應最大程度上考慮分解后的模態分量的故障信息。因此，本文提出最小平均復合影響指數(minimum average composite impact index, MACII)作為目標函數，可表示為

(18)

針對變分模態分解參數取值難以確定，嚴重依賴經驗選取而導致解析力低的問題。本文提出以最小平均復合影響指數為目標函數，將蝙蝠算法對VMD參數進行全局尋優，自適應地實現信號的頻域剖分及各分量的有效分離，具體流程如圖3所示。

自適應變分模態分解實施步驟為[15]：

步驟1初始化蝙蝠算法參數，設置VMD參數K和α的范圍，確定適應度函數；

步驟2將信號進行變分模態分解，根據式(18)計算信號的平均復合影響指數；

步驟3將平均復合影響指數作為適應度函數進行全局搜索，尋找最小值；

步驟4根據式子更新蝙蝠個體的位置和速度；

步驟5重復步驟2～步驟4，直至確定最小復合影響指數或達到所設定的迭代循環次數，輸出最佳蝙蝠個體(K，α)；

步驟6使用具有最佳參數的VMD分解原始信號。

3 仿真分析

為了驗證該方法的可行性，建立滾動軸承復合故障特征模型，進行自適應變分模態分解，對分解后的信號進行包絡解調分析，從包絡譜中識別出復合故障特征。

考慮到機械系統工況復雜，采集到的振動信號信噪比低，存在各種干擾成分，特別是由軸旋轉或齒輪嚙合引起的隨機沖擊，嚴重阻礙了故障脈沖的提取。為了模擬復雜工況下采集到的真實振動信號，在振動信號中加入背景噪聲、隨機沖擊和周期性諧波。以滾動軸承為例，其復合故障模擬信號由以下5部分組成[16]

(19)

在式(19)中：第一項和第二項分別代表軸承外圈故障和內圈故障造成的周期性沖擊；第三項表示由背景和電磁干擾等產生的隨機沖擊；第四項表示周期諧波，最后一項表示高斯白噪聲。在該振動模型中，A，B，R和P是不同項的振幅，T是兩個相鄰脈沖之間的時間間隔，τ是由于滾珠隨機滑動效應而引起的時間滯后，它通常占T的1%～2%，用于模式振幅調制。fh表示諧波干擾的頻率，φh為相位。s(t)為脈沖響應的函數，其表達式為

s(t)=e-βntsin(2πfnt+φ)

(20)

式中：βn和fn分別為阻尼系數和共振頻率；φ為相位。

通常，軸承外圈故障沖擊的振幅(式(19)中的Ai)為常數，軸承內圈故障沖擊振幅計算如下

B=1+Bbsin(2πfrt+φb)

(21)

本文構建了含有高斯噪聲的的復合故障模擬信號，如圖4所示。外圈故障、內圈故障、隨機沖擊、周期諧波、高斯噪聲和混合信號的時域波形如圖4(a)～圖4(f)所示。關鍵參數值如下：采樣速率為12 kHz，采樣數量為7 200。軸旋轉頻率fr為30 Hz，內外圈故障特征頻率fi和fo分別為110 Hz和80 Hz(Ta=1/fo，Tb=1/fi)。外圈故障、內圈故障和隨機沖擊共振頻率分別為3 000 Hz，4 000 Hz，4 500 Hz。其阻尼系數分別為1 000，1 200，1 500。Ai和Bb分別為1.2和1.4，τi和τj分別為0.01Ta和0.01Tb。Rm和Tr由MATLAB中randn函數隨機選擇產生，Rm的最大值為3。周期諧波部分寫為Psin(2πf1t)cos(2πfrt)，其中振幅P為0.01，頻率f1為10 Hz，高斯白噪聲的信噪比為5 dB。

如圖5所示，由于背景噪聲的干擾，頻譜圖中譜線較為混亂，在3 000 Hz和4 000 Hz出現共振頻帶，無法識別出軸承內圈和外圈的故障頻率，得不到所需的有效信息。包絡譜中雖然能夠找到內、外圈故障特征頻率，但是這些脈沖相當隨意，沒有呈現規律變化，干擾譜線較多，加大復合故障分離難度。

為了能夠精確識別出軸承復合故障的類型，通過蝙蝠算法對變分模態分解進行參數尋優，設置蝙蝠迭代次數為10，種群數量為20。平均復合影響指數隨蝙蝠更新代數變化趨勢如圖6所示，從圖6中可以看出，隨著迭代次數增加，適應度值逐漸減小，這說明采用優化后的參數后，分解結果逐漸接近最優值，最后，通過不變的適應度表明，找到了最優參數。當迭代次數為3時，適應度值最小，此時模態數量和懲罰因子分別為(5,100)。選取最優參數對仿真信號進行變分模態分解，結果如圖7所示，從頻譜圖中，頻率從低到高分布，沒有出現頻率混疊的現象。從包絡譜中可以看出，故障脈沖信號分別位于3模式和4模式，其他模式相對于故障診斷是冗余的。從模式3中，可以清楚地找到外圈故障的頻率以及多個倍頻。模式4中同樣可以找到內圈故障的頻率及多個倍頻，在其倍頻兩側均勻分布著譜線，與倍頻間隔都是fr，可見轉頻信號對外圈故障信號起著調制作用。為了消除轉頻對診斷結果的影響，需要結合其他降噪方法進行進一步的研究。

4 案例研究

采用凱斯西儲大學的滾動軸承試驗數據驗證本文所提方法的有效性，試驗臺如圖8所示。試驗臺由驅動電機、扭矩傳感器、控制器和故障軸承等組成，電機驅動端安裝有SKF6205型號的深溝球軸承，加速度傳感器安裝在軸承的正上方[17]。測試軸承通過電火花加工的方式在其內圈，外圈加工直徑為0.177 8 mm，深度為0.279 4 mm的單點模擬故障，試驗條件如表1所示。在本試驗中，采樣頻率為12 kHz，采樣時間為1 s，采用Miao等的試驗方法，將軸承內圈故障信號和外圈故障信號混合得到復合故障信號。考慮到現場噪聲干擾的隨機性，在試驗信號中加入了5 dB的高斯白噪聲。

根據表1中的參數，內圈、外圈故障特征頻率由式(22)和式(23)計算得到

表1 試驗條件Tab.1 Test conditions

fbpfi=5.415 2×fr=162.185 Hz

(22)

fbpfo=3.584 8×fr=107.365 Hz

(23)

如圖9所示，周期脈沖信號雜亂無章難以識別，埋沒在噪聲之中。圖10為試驗信號的包絡譜圖，包絡譜中雖然能找到內外圈故障頻率，由于受到背景噪聲的干擾，周圍存在眾多譜線，影響軸承故障類型的診斷。因此，需要從噪聲信號中進一步提取故障特征。

首先需要對VMD的模態數量和二次懲罰因子進行優化。設置種群數量為20，迭代次數為10。圖11顯示了搜索結果，隨著迭代次數的增加，平均復合影響指數逐漸減小。當迭代次數為4時，得到最優參數，此時K和α分別為(5,429)。采用最佳參數對信號進行變分模態分解，分解信號如圖12所示，對分解后的信號進行包絡解調分析，如圖12(b)所示，頻譜圖中沒有出現模態混合和復制的現象。從圖12(c)中可以看出，故障脈沖信號分別位于模式3和模式4，其他模式相對于故障診斷是冗余的。從模式3中，可以清楚地找到外圈故障的頻率、2倍頻以及3倍頻，模式4中同樣可以找到內圈故障的頻率及2倍頻，同時混雜著轉頻及其倍頻，需要進行一步的研究，去除轉頻的影響。

5 結 論

通過本文的研究，可以得到以下結論：

(1) 本文提出了一種基于蝙蝠算法優化VMD參數的滾動軸承復合故障分離方法。該方法克服了VMD參數需要提前確定的問題，可以自適應地得到與待分解信號匹配的模態數量和懲罰因子的最優參數組合。

(2) 在相關系數、功率譜、峭度等指標基礎上，提出了一種新的復合影響指數CII，將其與現有指標KES、KESE、EE、EK、SII進行比較，該指數對故障的敏感度提升了29.6%，同時能夠避免周期諧波、隨機沖擊、背景噪聲的干擾，是一種比較適合識別故障沖擊的指數。

(3) 改進的自適應VMD方法能夠消除噪聲的干擾，識別出復合故障的特征頻率。但是在強噪聲干擾下，包絡譜中干擾譜線較多，分解效果不夠理想。因此，將改進的自適應VMD方法與其他降噪方法結合以實現故障特征信息與噪聲的有效分離將值得進一步研究。