自適應VMD算法在滾動軸承故障診斷方面的應用*

王 杰，郭世偉

（西南交通大學機械工程學院，成都 610031）

0 引言

滾動軸承是機械設備的關鍵部件之一，也是易損傷的零部件之一，其工作狀態很大程度影響整臺機械設備的性能[1]。及時診斷出滾動軸承的問題所在，保障機械設備安全可靠地運行，具有十分重要意義。

變分模態分解（VMD）是由Dragomiretskiy.K[2-3]在2014年提出的一種新的自適應時頻分析法，該方法是通過迭代更新模型中的最優解來確定每個模態分量的中心頻率和帶寬，從而自適應地實現將各個分量從頻域分離；唐貴基[2]、馬洪斌[4]提出了以模態分量的包絡熵為適應度函數，并分別采用粒子群算法、蛙跳算法優化模態分量個數k和懲罰因子α；宋坤駿[5]提出以局部最大峭度值所對應模態分量的峭度為目標函數，采用網格搜索法優化k和α這兩個參數。

對于VMD參數選取問題，本文根據最大峭度指標和最小包絡熵組合成復合指標，結合粒子群優化算法對VMD的最佳參數進行搜尋，提出了自適應VMD算法的滾動軸承故障診斷方法。通過仿真和實驗數據的分析，本文提到的方法能有效提取故障特征信息，準確判斷出滾動軸承的故障類別。

1 變分模態分解（VMD）基本原理

假設VMD每個模態函數都有不同中心頻率的帶寬，而中心頻率的帶寬隨著算法的迭代而更新[3]。在各個模態分量互相迭代后能還原成原始信號的前提條件下，求解出k個估計帶寬之和為最小的模態函數uk(t)。

式中：δ(t)為單位沖激函數，{uk}={u1,u2,…,uk}為分解得到的各個模態函數，f(t)為原始信號，{ωk}={ω1,ω2,…,ωk}為各個模態函數的中心頻率；為所有模態函數之和。

Dragomiretskiy和Zosso通過引入二次懲罰因子α和拉格朗日乘子λ(t)，將約束變分求解問題轉化為非約束變分求解問題。將式（1）轉化為增廣拉格朗日表達式：

使用交替方向乘子法（ADMM）依次迭代uk、ωk、λ，直到尋找到增廣拉格朗日鞍點。將uk、ωk、λ的迭代公式轉換到頻域內再計算，最終得到模態函數和中心頻率的更新方法。

2 自適應VMD算法

2.1 粒子群優化算法

粒子群優化算法[6]（PSO）是由J.Kennedy和R.C.E-berhart在1995年提出的群體智能算法，用來模擬鳥群在覓食過程中的遷徙和聚集行為。其具有快速收斂的能力，同時也具有良好的全局尋優能力。具體算法過程如下。

（1）步驟一：初始化粒子群的規模N，學習因子和慣性因子的調整參數，粒子的速度v，粒子位置x，粒子的邊界條件，迭代次數t=1，最大迭代次數Tmax。

（2）步驟二：計算每個粒子的適應度函數ff(t)，取ff(t)的最小值作為全局最優位置gbest0以及局部最優位置pbest0的初始值。

（3）步驟三：根據式（5）～（9），逐步更新學習因子、慣性因子、粒子速度、粒子位置。

式中：R1、R2、R3、R4分別為學習因子c1、c2的調整參數；ωmax、ωmin為慣性因子ω的調整參數；t 為當前迭代次數；r1、r2為服從[0,1]均勻分布的隨機數。

（4）步驟四：計算每個粒子的適應度函數ff(t)，取ff(t)的最小值作為本次迭代的局部最優位置pbest1，與上一次迭代的全局最優位置gbest0進行比較。若pbest1

（5）步驟五：若當前迭代次數t 大于總的迭代次數Tmax，則輸出全局最優位置以及全局最優位置的粒子P；否則，返回步驟三。

2.2 VMD參數自適應優化

峭度指標Kr對沖擊信號敏感，其值不受外界環境影響而改變，非常適用于滾動軸承表面損傷類故障診斷[7]。

式中，N為在這段時間內采集點數；xˉ為均值，

包絡熵Ep反映原始信號的稀疏特性[2]。

式中：a(j)為原始信號經希爾伯特解調后的包絡信號，j=0,1,…,N。

根據最大峭度指標和最小包絡熵這兩個參數共同優化VMD的模態分量個數k和懲罰因子α，使分解得到的模態分量含有較多的周期沖擊性信號成分，這些模態分量盡可能表現出滾動軸承故障特征信息。VMD分量中的最小包絡熵Ep與VMD分量中的最大峭度指標Kr相除，構造出同時反映周期性沖擊性信號成分的變化情況的復合指標KEP。若KEP 值越大，原始信號中的周期性和沖擊性信號比重越低，越不能判定該信號為故障信號；若KEP 值越小，原始信號中的周期性和沖擊性信號比重越高，該信號為故障信號的可能性就越大。

KEP 作為PSO 的適應度函數，對VMD 算法中懲罰因子α和模態分量個數k 進行自適應優化。將k 和α組合成一個粒子P=[k,α]，更新粒子局部最優值和全局最優值，自動篩選出最優粒子從而獲得k與α最優解的組合。該過程稱為自適應VMD算法，程序設計流程如圖1所示。

圖1 自適應VMD算法程序設計流程

3 基于自適應VMD算法的滾動軸承故障診斷

首先根據自適應VMD算法對待測信號進行分解，然后計算各個模態分量的峭度指標，選取最大的峭度指標作為最佳模態分量，最后應用Hilbert 算法對最佳模態分量進行包絡譜分析[8]，判斷包絡譜中故障特征頻率及其倍頻處的幅值是否突出，再根據故障特征頻率判斷其故障類別。

圖2 滾動軸承故障診斷流程

4 仿真分析

為了驗證本文采用方法的可行性，選用恒定轉速的軸承故障仿真模型來驗證[9]。

式中：x(t)為軸承故障模擬信號；A 為幅值，取1；故障特征頻率fo=1/T=155 Hz；B 為衰減系數，取800；fn為共振頻率，取3500 Hz；n（t）為高斯白噪聲。

假設仿真信號的采樣頻率fs=12 kHz，分析數據點數為4800 個，得到軸承故障模擬信號，如圖4 所示。通過圖3～4 對比可知，軸承故障模擬信號淹沒在噪聲背景環境中，無法觀察到軸承故障沖擊信號。

圖3 軸承故障沖擊信號

圖4 軸承故障模擬信號

圖5 KEP隨粒子種群變化曲線

圖6 軸承故障模擬信號VMD分解

圖7 IMF6包絡譜

圖8 軸承故障模擬信號包絡譜

首先根據自適應VMD 算法，對VMD 算法中的k 和α進行尋優。PSO 的初始參數按照參考文獻[10]進行選取，如表1 所示。基于PSO 優化KEP 的迭代過程，如圖5 所示。當迭代次數t 為11時，KEP 達到最優值0.9875，粒子達到最優位置P=[8,3690]；當k=8 和α=3690時，VMD 獲得最優參數的組合。然后設置VMD算法的優化參數，再對該信號進行分解，得到8個模態分量的時域信號，如圖6所示。最后計算8個模態分量的峭度指標，IMF6的峭度指標最大，并對IMF6進行Hilbert包絡解調處理，得到包絡譜，如圖7所示。從包絡譜中可以清晰觀察到155 Hz 以及二倍頻處的幅值，與故障特征頻率fo=155 Hz一致，能判斷出其故障類型。

應用Hilbert 算法對軸承故障模擬信號進行包絡解調，得到包絡譜，如圖8 所示，故障特征頻率155 Hz 處的幅值凸出顯示出來，但倍頻處的幅值淹沒在噪聲環境中，未清晰反映出來，弱于本文采用方法的故障特征信息提取效果。

表1 PSO的基本參數設置情況

5 實驗分析

本實驗選用滾動軸承故障診斷實驗臺，如圖9所示，來繼續驗證本文采用方法的可行性。在恒定轉速n=1200 r/min 條件下，使用某公司研發的YD-61D 壓電加速度傳感器以及配套的采集分析系統，對N205外圈故障的振動信號進行采集與存儲。在該轉速下，N205外圈故障特征頻率fo1=109.64 Hz。YD-61D 安裝在故障軸承端的軸承座的豎直和水平方向，取豎直方向上的連續采樣時間為0.4 s的外圈故障振動信號進行分析，如圖10所示。由圖可知，外圈故障信號淹沒在噪聲中，無法明顯觀察到周期性沖擊信號成分，很難判別該軸承的故障類型。

圖9 滾動軸承故障診斷實驗臺

圖10 外圈故障信號時域圖

圖11 KEP 隨粒子種群迭代變化曲線

首先根據自適應VMD算法，對VMD參數進行優化。PSO的初始參數按照參考文獻[10]進行選取，如表1 所示。基于PSO 優化KEP 的迭代過程，如圖11 所示。當迭代次數t 為10時，KEP 達到最優值3.491，粒子達到最優位置P=[9,697.6]；當k=9和α=697.6時，VMD獲得最優參數組合。然后應用優化參數后的VMD 算法對該信號進行分解，得到9 個模態分量的時域信號，如圖12 所示。最后計算9 個模態分量的峭度指標，IMF8 的峭度指標最大，并對IMF8 進行Hilbert 包絡解調處理，得到包絡譜，如圖13 所示。從包絡譜中可以觀察到110 Hz 以及2～5 倍頻處的幅值凸出顯示出來，與故障特征頻率fo1=109.64 Hz 十分接近，可以判斷出該故障類型為外圈故障。

應用Hilbert 算法對外圈故障信號進行包絡解調，得到包絡譜，如圖14所示。通過圖13～14對比可知，采用本文方法取得較好的故障特征信息提取效果，有利于軸承故障類型的判別。

圖12 外圈故障信號VMD分解

圖13 IMF8包絡譜

圖14 外圈故障信號包絡譜

6 結束語

根據故障軸承的周期沖擊性振動信號特征構造出KEP 參數，采用KEP與PSO相結合的參數優化方法，對VMD參數進行尋優，避免了盲目性選擇k和α這兩個重要參數；采用自適應VMD算法對故障信號進行分解，根據最大峭度指標選取最佳模態分量，并進行包絡譜分析，能準確地判斷出故障類型。通過仿真與實測實驗結果的分析，驗證了該方法在軸承故障診斷方面的可行性，為軸承故障診斷提供了指導意見。