基于FVMD的滾動軸承故障特征提取方法

張 爽， 王曉東， 李 祥， 楊創艷

(1.昆明理工大學 信息工程與自動化學院，昆明 650500； 2.昆明理工大學 云南省人工智能重點實驗室，昆明 650500)

滾動軸承作為機械設備的關鍵零部件，常常工作在高溫、高壓和復雜的力學環境中，極易發生故障[1]。軸承一旦發生故障，就會導致設備停機，生產線難以正常運行，造成重大的經濟損失，甚至導致災難性的人員傷亡[2]。因此，提取滾動軸承的故障特征對及時有效判定其工作狀況和保證生產安全具有重大意義。

振動信號分析作為滾動軸承故障診斷中最受歡迎的技術手段。其中，經驗模態分解(empirical mode decomposition,EMD)是最具代表性的方法之一，因其在處理非線性、非平穩信號時有著良好的效果，并在滾動軸承故障診斷領域中得到了廣泛的應用[3-4]。但同時存在模態混疊，端點效應等問題，嚴重影響其分解效果[5]。為此，有學者提出了局部均值分解(local mean decomposition,LMD)方法，可很好地改善EMD所遺留的模態混疊，端點效應等問題，保證故障特征提取正確性。但是，LMD本質上仍然是遞歸模態分解，容易受采樣頻率的影響，難以保證分解精度[6]。為解決EMD和LMD遞歸模態分解的問題，Dragomirestkiy等[7]提出變分模態分解(variational mode decomposition,VMD)算法，該算法對信號進行分解時采用的是一種非遞歸的處理技術，彌補了循環遞歸篩選的缺陷，從而有效地擺脫LMD與EMD存在的端點效應、模態混疊等束縛[8]。

VMD能克服端點效應和模態混淆的弊端，但涉及到3個參數：模態數K、懲罰因子α、模態初始中心頻率ω的預設置。模態初始中心頻率ω在傳統VMD中給出了3種初始化方式，其中ω等于零和ω在被分解信號的頻譜上呈線性分布的方式具有較高的計算效率，另一種ω隨機取值的方式效率較低。Bi等[9]將信號頻譜中幅值最大點對應的頻率作為模態初始中心頻率ω，并把模態數K設為1進行分解，有效的減少VMD的模態中心頻率迭代次數。另一參數K太小會導致分解結果不完全，故障信息被遺漏，而模態數K太大會導致過分解[10]。懲罰因子α越小，分解后的本征模態函數(intrinsic mode function,IMF)分量衰減越慢、帶寬越大，相反，分解后IMF分量衰減越快、帶寬越小[11]。Ma等[12]通過尺度空間頻譜分割的方法來確定模態數K的值，同時把懲罰因子α默認為2 000，實現了參數K的自適應選擇。劉長良等[13]將模態數K設置為2、懲罰因子α設為2 000，并計算VMD分解后IMF的中心頻率，若出現了中心頻率相近的模態則認為出現過分解，模態數K=K-1,否則K=K+1再次計算IMF的中心頻率。Jiang等[14]把模態數K設為1、懲罰因子α設為2 500進行信號的分解，并把殘余信號依次分解，計算分解出模態的峭度值，當峭度值開始下降時停止分解，對峭度值最大的模態進行分析。

上述大部分研究把懲罰因子α設為2 000，對模態數K進行優化，但沒有考慮到懲罰因子α變化帶來的影響。目前，大部分研究采用優化算法對VMD的參數K和α同時進行選擇：張俊等[15]以包絡譜峰值因子作為優化指標，采用粒子群優化算法尋找最優的模態數K和懲罰因子α；Zhang等[16]把權重峭度作為優化指標，通過蚱蜢優化算法搜尋模態數K和懲罰因子α。利用優化算法和VMD相結合，能夠對參數進行有效的設置，但使得整個流程的計算效率降低、忽略了VMD算法的性質。綜合上述存在的問題，本文從對模態初始中心頻率ω的定位入手，提出了一種快速變分模態分解(fast VMD,FVMD)算法。算法對不同ω的模態計算出對應的α值，能同時對K和α值進行優化。相比于原VMD算法的模態初始中心頻率ω定位方式，還能夠增加其計算效率。

1 基本理論

1.1 VMD

(1)

式中：δ(t)為沖擊函數；f(t)為輸入信號； j為虛數單位；?t為對參數t求偏導。VMD算法可以將輸入信號分解為K個以ωk為中心頻率的IMF分量uk(t)。VMD算法可以分為兩個主要部分：第一個部分是變分問題的構造；第二部分是變分問題的求解。約束變分問題的構造可以用式(1)來進行描述。

然后，引用拉格朗日乘子λ和二次懲罰因子α，將約束變分問題轉化為無約束變分問題進行求解，其增廣拉格朗日如式(2)所示。

(2)

采用乘子交替方向法求解增廣拉格朗日表達式，模態uk(t)及模態中心頻率ωk的更新步驟如下所示。

步驟2n=n+1。

(3)

步驟4對拉格朗日算子λ進行了更新如式(4)。

(4)

τ為噪聲容限參數，為了獲得更好的去噪效果，可設τ=0。

重復步驟2～步驟4，直到滿足式(5)則停止迭代。

(5)

式中，ε值通常設置為1×10-6，詳細內容請參考Liu等的研究。

1.2 閾值去噪

閾值去噪是小波閾值去噪的主要部分，進行小波閾值去噪時，首先需要設置一個臨界閾值λ。當小波系數小于λ，認為該系數主要由噪聲引起，去除這部分系數；若小波系數大于λ，則認為此系數主要是由信號引起，予以保留。因此，將此想法應用到頻譜趨勢的計算中，消除噪聲引起的頻譜趨勢波動。

此處采用硬閾值去噪的方法對頻譜趨勢進行去噪處理

(6)

閾值的計算公式如式(7)[17]所示

(7)

1.3 頻譜趨勢分割

快速經驗小波變換(fast empirical mode decomposition,FEWT)[18]通過快速傅里葉變換(fast Fourier transform,FFT)和快速傅里葉反變換(inverse FFT,iFFT)計算信號的頻譜趨勢，對頻譜趨勢使用閾值去噪進行優化，取優化后頻譜趨勢的極小值點作為分割邊界。本文引用頻譜趨勢分割的方法，來確定VMD算法中模態數K。同時，VMD優先對頻譜中能量大的調頻調幅部分進行分解，將故障信號中調頻調幅部分的中心頻率作為VMD模態初始中心頻率能夠降低模態迭代次數。而頻譜趨勢能夠凸顯此部分，因此，以頻譜趨勢極大值對應的頻率作為模態分量的初始中心頻率。其方法流程如圖1所示，結合圖1和式(8)，方法實現過程闡述如下。

步驟1通過FFT計算得到信號的頻譜，如圖1(a)所示。對信號的頻譜再次使用FFT得到信號頻譜的核函數，如圖1(b)所示。

步驟2取核函數的前30個點(圖1(b)虛線左邊的部分)，利用iFFT計算得到頻譜趨勢(圖1(c)頻譜上方曲線)，頻譜趨勢能夠刻畫故障信號在頻譜上的整體走向，并根據其組成成分能量的大小，頻譜趨勢的高低也隨之變化。并把趨勢的極小值作為分割的邊界，如圖1(c)虛線所示。

步驟3利用硬閾值去噪的方法來優化頻譜趨勢，被噪聲干擾的頻譜趨勢部分呈直線，然后以頻譜趨勢極小值點和直線趨勢的端點作為新分割邊界，如圖1(d)虛線所示。將圖1(d)中的兩段直線趨勢所在頻段視為無效頻段(a段和b段)，其余帶有趨勢極大值點的4段作為有效頻段，即模態數K=4。

步驟4通過優化的頻譜趨勢分割圖，可以容易的找到頻譜趨勢極大值點，將極大值點對應的橫坐標頻率作為模態初始中心頻率即ω=[ω1ω2ω3ω4]，如圖1(d)所示。

圖1 頻譜趨勢分割

(8)

式(8)由諧波和周期脈沖組成，同時加入高斯白噪聲來模擬實際信號中的噪聲干擾，諧波信號y1、y2、y3的頻率分別為40、400、900，旋轉頻率fr為25 Hz，故障頻率fm為70 Hz，衰減系因子ζ為1 000，共振頻率fn為2 000 Hz，n(t)為零均值1方差的高斯白噪聲，采樣頻率為10 kHz，采樣點數N為4 096。

1.4 懲罰因子的計算

根據滾動軸承故障振動信號的頻譜分布特征，中低頻區域主要由旋轉頻率及其相關特征頻率(如軸承故障特征頻率和齒輪嚙合頻率等)的諧波組成，而故障沖擊和噪聲干擾大多位于高頻區域[19]。同時諧波信號具有時域持續時間長、頻域上較為緊湊的特點，而沖擊信號具有時域短、頻域寬的特點。

因此，以各模態分量的中心頻率作為確定相應懲罰因子的依據。通過頻譜趨勢分割得出的模態初始中心頻率較小，說明模態提取的部分主要是諧波，懲罰因子應該大；相反，說明模態提取的部分主要是故障沖擊部分，此時的懲罰因子應該選擇較小的值。因此，可以建立式(9)的映射關系

(9)

1.5 有效權重峭度指標

從分解的K個IMF分量中篩選含豐富故障信息的分量是另一個需要解決的核心任務。峭度指標通常作為沖擊測量指標用于監測機械損傷[20]。但是，峭度指標僅僅考慮沖擊信號的分布密度，而忽略振幅較大、分布不均勻的沖擊分量。相關系數可以表征兩個信號的相似性，能夠考慮到沖擊分量的不均勻性，但在沖擊信號檢測中容易受到噪聲的影響。

結合兩個指標的特點，引入有效權重峭度指標來篩選分量，該指標既考慮到脈沖分量的不均勻性，又有效地濾除了噪聲[21]。該指標可以根據模態分量自身的情況，從中選出含沖擊信息多、噪聲少的有效分量，計算公式如式(10)所示。

(10)

式中，Ki和Ci、KEW,i分別為第i個模態分量的峭度值和相關系數、有效權重峭度指標，i,j=[1,2,3，…，K]。將KEW(k)>0的模態判斷為有效模態分量，其余為含噪模態分量。

2 FVMD的故障特征提取方法

根據第1章的理論分析，本文提出了一種FVMD的滾動軸承故障特征提取方法。方法的流程如圖2所示，結合圖2，方法實現過程闡述如下。

圖2 FVMD故障特征提取流程

步驟1收集滾動軸承故障信號，用頻譜趨勢分割法進行預處理。

步驟2通過頻譜趨勢分割，確定VMD算法參數ω、K、α，并完成信號的自適應分解。把有效頻段中的極大值所對應頻率ωk作為模態初始中心頻率，ω=[ω1ω2ω3ω4…ωk],1≤k≤K；通過式(9)計算出每個初始中心頻率ωk對應的懲罰因子αk，α=[α1α2α3α4…αk]，1≤k≤K；改變模態初始中心頻率和懲罰因子初始化設置方式，利用優化后的K、ω和α更新傳統VMD參數，并執行信號的VMD分解。

步驟3采用有效權重峭度準則，選取有效模態分量進行重構，并對重構信號進行包絡解調完成故障特征提取。計算分解后各模態的峭度值和相關系數；根據式(10)計算各模態的有效權重峭度指標，將有效權重峭度指標大于零的模態分量進行重構；包絡解調獲取故障特征。

3 仿真試驗對比分析

為了驗證所提出方法的有效性，對式(8)組成的仿真信號進行分析，其波形如圖3所示。

圖3 仿真信號

根據圖3(c)的頻譜趨勢分割圖可以快速獲得FVMD的ωF=[0,460,937,2 031]和K=4,并根據式(9)計算出對應的懲罰因子αF=[2 500,1 191,871,482]，從而分解得到模態，如圖4(a)所示。并與傳統VMD算法中的兩種模態初始中心頻率設置方式進行對比，論證FVMD方法動態設置ω和K的有效性。①是ω=0，懲罰因子α設置為2 000、模態數K=4,如圖4(b)所示；②是ω在頻譜上呈線性分布，懲罰因子α設置為2 000、模態數K=4，如圖4(c)所示。此外，圖4(d)還給出ω=ωF和懲罰因子α=2 000、K=4的分解結果，論證了懲罰因子α動態設置的必要性。

同時，圖4(b)與圖4(d)分解出4個故障相關模態，而圖4(c)未分解出第二個諧波，驗證了模態初始中心頻率的合理設置的必要性。

圖4 4種VMD分解方式的IMF分量

為了選取有效分量進行后續的頻譜分析，分別計算圖4(a)～圖4(d)中IMF分量的EWK(effective weight kurtosis)，如表1所示。以及對應IMF分量的中心頻率迭代次數，如圖5所示。根據有效權重峭度準則，選取表1中EWK>0的分量進行重構，并計算包絡解調頻譜，如圖6所示。

表1 4種VMD分解方式下IMF分量的EWK

圖6 4種VMD分解方式的重構信號包絡解調

分析圖5～6可知：①圖5(a)FVMD根據頻譜趨勢分割設置的ωF，相比于圖5(b)～圖5(c)傳統VMD中ω的設置方式，模態中心頻率的迭代次數最少，進而提高了VMD的計算效率。在優化ω的基礎上，同時對懲罰因子進行優化，相比α=2 000進行分解的方式(見圖5(d))增加了迭代次數,但是FVMD方法迭代次數遠低于傳統VMD方法。②圖6(a)所提方法成功提取出故障的基頻和2倍～6倍頻。而圖6(b)～圖6(d)中令α=2 000進行VMD分解的方法只提取到故障基頻和2倍～4倍頻。說明用此方法對α進行優化，有利于取到更豐富的故障特征。③同時計算4種方法重構信號的峭度值，可以看出圖6(a)FVMD的峭度值最大，說明了FVMD的有效性。

4 試驗數據對比分析

為了驗證基于FVMD的滾動軸承故障特征提取方法在實際工程中的有效性和優越性，以下試驗使用CWRU(Case Western Reserve University)數據集和NASA(National Aeronautics and Space Administration)數據集并基于表2中的方法進行分析。

表2 試驗方法

4.1 CWRU滾動軸承振動數據分析

CWRU滾動軸承故障模擬試驗平臺如圖7所示，具體參數如下：軸承型號為SKF6205，轉速為1 750 r/min，采樣頻率為12 kHz，采樣點數為2 048。選取滾動軸承內圈、外圈故障振動信號進行試驗分析，滾動軸承外圈理論故障頻率BPFO=107.305 Hz、內圈故障頻率BPFI=162.185 Hz。

圖7 CWRU軸承數據試驗平臺

4.1.1 軸承外圈故障信號試驗對比分析

軸承外圈故障信號與頻譜趨勢分割如圖8所示，可計算得到模態數K=3，初始中心頻率ωF=[375，2 813，3 375]、懲罰因子αF=[1 616，489，372]。將其參數代入VMD中分解得到系列IMF分量，如圖9(a)所示。同時給出第3章所述的其他3種方法的IMF分量，如圖9所示。

圖8 外圈故障信號

圖9 4種VMD分解方式的IMF分量

從頻域上看，4種方法都分解出3個故障相關模態，與頻譜趨勢分割中的有效頻段相對應。但僅通過圖9對比分析無法直觀看出區別，需要進一步計算圖9(a)～圖9(d)中各IMF分量的EWK來對比分析，計算值如表3所示。同時，給出各IMF分量的中心頻率迭代次數，如圖10(a)～圖10(d)所示。由表3可知，將EWK>0的IMF分量視為有效分量予以保留，并進一步對有效分量重構后進行包絡解調，如圖11(a)～圖11(d)所示。

表3 4種VMD分解方式下IMF分量的EWK

對圖10～11進行對比分析可知：①圖10(a)FVMD方法各模態的中心頻率迭代次數最少為14，均小于圖10(b)～圖10(c)傳統VMD方法的兩種初始化方式，故本文所提方法能有效減少VMD算法的迭代次數。②由圖11(a)～圖11(d)可知，4種方法都存在與外圈故障頻率理論值(107.305 Hz)非常接近的頻率105.5 Hz，還存2倍～7倍頻。但是圖11(a)FVMD方法提取到的故障基頻和倍頻幅值都大于其他3種方法。③同時，從圖11中可知，所提方法的峭度值最大，進一步驗證了動態設置α不僅能保留更多的故障信息，而且能突出振動信號的故障沖擊分量。

圖10 4種VMD分解方式的IMF分量中心頻率迭代次數

圖11 4種VMD分解方式的重構信號包絡解調

4.1.2 軸承內圈故障信號試驗對比分析

按照第3章的分析思路，應用頻譜趨勢分割計算得到內圈故障振動信號的模態數K=4、模態初始中心頻率ωF=[562，1 125，2 625，3 563]。根據式(9)計算懲罰因子αF=[1 420，1 045，533，338]，將其參數K、ωF和αF代入進行內圈振動信號的VMD分解，并給出其他3種分解方式的IMF分量。計算分解后各IMF分量的EWK值，將大于零的分量進行重構。并計算其包絡解調，如圖12所示。

圖12 4種VMD分解方式的重構信號包絡解調

由圖12可知：①通過4種方式進行信號分解，都提取到與內圈故障頻率理論值(162.185 Hz)非常接近的頻率164.1 Hz，還存2倍～7倍頻。但是，圖12(a)FVMD方法提取到的故障特征頻率幅值最大、故障特征最明顯。②進一步計算4種方法重構信號的峭度指標峭度值，圖12(a)所提方法的峭度值最大，具有最豐富的故障特征。圖11與圖12對比可知，外圈與內圈故障振動信號試驗取得一致的結論，進一步論證所提方法(FVMD)的有效性和可行性。

4.2 NASA滾動軸承振動數據對比分析

NASA軸承全壽命周期試驗平臺如圖13所示。采樣頻率為20 kHz，驅動電機轉速為2 000 r/min，徑向負荷2 721.554 22 kg，共采集984組數據,每組長度為20 480。試驗采用外圈第874組數據進行分析。

圖13 NASA軸承數據試驗平臺

同理，采用FVMD方法對故障振動信號進行處理，獲得信號模態數K=5、對應的初始中心頻率ωF=[0，937，1 563，3 438，4 375]、懲罰因子αF=[5 000，2 365，1 913，1 139，888]，4種方法的重構信號包絡頻譜如圖14所示。

圖14 4種VMD分解方式的重構信號包絡解調

結合圖14進行分析可知：①圖14(a)～圖14(d)4種分解方法都能提取到故障的基頻和倍頻。但是，圖14(a)FVMD方法提取到的基頻和倍頻幅值均大于圖14(b)～圖14(d)(α=2 000、ω=0或ω=呈線性或ω=ωF)的方法。②進一步計算重構信號的峭度指標峭度值，對比分析可知所提方法(FVMD)取得最大的峭度值，進一步驗證FVMD方法的優越性。

5 結 論

本文通過仿真信號和CWRU、NASA滾動軸承數據的試驗及對比分析，驗證了基于FVMD的滾動軸承故障特征提取方法的有效性，可以得到以下結論：

(1)針對傳統VMD方法需要預先給定分解模態數K和懲罰因子α的局限性，本文結合頻譜趨勢自適應分割方法，實現了VMD分解模態數K、懲罰因子α和模態初始中心頻率ω的自適應選擇。

(2)本文通過對模態初始中心頻率ω和懲罰因子α進行動態優化降低了模態中心頻率的迭代次數，提高了計算效率；同時，所提方法比傳統VMD方法提取到了更豐富故障特征。

(3)針對如何從IMF分量中篩選有效分量的問題。結合峭度指標和相關系數指標建立了有效權重峭度指標篩選規則，實現了有效分量的選取。

(4)通過仿真信號和CWRU、NASA滾動軸承數據的對比分析，有效論證了FVMD滾動軸承故障特征提取方法的可行性和有效性。