麻雀算法參數(shù)優(yōu)化VMD聯(lián)合K-SVD滾動軸承故障診斷

褚 惟，王貴勇，劉 韜，王振亞

（1.昆明理工大學 機電工程學院，昆明 650500；2.內(nèi)蒙古第一機械集團有限公司，內(nèi)蒙古 包頭 014000）

滾動軸承作為機械裝備中的重要零部件被廣泛應用于工程領域[1]。軸承的健康狀態(tài)與機械裝備的正常運轉(zhuǎn)緊密聯(lián)系，一旦發(fā)生故障便會影響到生產(chǎn)安全和經(jīng)濟效益。由于滾動軸承的振動信號蘊含了豐富的設備狀態(tài)信息，因此對軸承進行振動分析已成為狀態(tài)監(jiān)測分析的基本手段[2]。

目前基于正交基的信號分解方法，如傅里葉變換[3]，小波變換[4]等，往往需要大量的基函數(shù)才能對故障信號進行完整表達，一定程度上難以簡潔和自適應地進行信號分解。而稀疏分解以過完備字典為基礎，能根據(jù)分析信號自身特點盡可能地優(yōu)選字典基函數(shù)來表示信號，使所得信號更為稀疏簡潔，避免了基函數(shù)的冗余，更好地刻畫信號的內(nèi)在特征。常用的稀疏字典有基于原始原子伸縮、平移等變換所創(chuàng)建的字典，如小波基字典，傅里葉基字典，離散余弦基(Discrete Cosine Transform，DCT)字典等[5]。這類方法的主要思想是設計具有特征相似度的字典，因此在與信號良好匹配的情況下是有效的，但處理某些復雜多變信號時，其性能會降低。另一類方法則是根據(jù)信號本身進行自適應學習，即基于信號本身的學習字典算法，如最佳方向法(Method of Option Directions，MOD)[6]，K-SVD[7]等。K-SVD算法通過奇異值分解進行字典更新，避免了MOD算法的大規(guī)模矩陣求逆過程，大大提高了運算效率[8]。

在工程環(huán)境中，由于工況條件的復雜多變，滾動軸承故障發(fā)生時信號往往具有非線性、非平穩(wěn)、弱周期性、低信噪比等特性，使得經(jīng)典的K-SVD 算法在自身學習過程中易受到噪聲干擾產(chǎn)生虛假原子從而影響表達效果[9]。VMD[10]通過限制信號各分量的帶寬，對不同中心頻率的信號成分進行分離，可有效地去除噪聲干擾。因此，以經(jīng)VMD分解的本征模函數(shù)(Intrinsic Mode Function，IMF)分量進行K-SVD字典學習來稀疏表達能避免K-SVD 學習過程產(chǎn)生的虛假原子，更有效進行故障識別。首先，通過麻雀搜索算法(Sparrow Search Algorithm，SSA)優(yōu)化VMD 的分解層數(shù)k和平衡因子α，并將所得優(yōu)化值代入VMD中進行分解；其次利用平方包絡譜峭度(Kurtosis of Squared Envelope Spectral，KSES)[11]遴選分解出最優(yōu)IMF分量；最后對最優(yōu)IMF分量進行K-SVD字典學習，提取有效故障信息。

1 理論基礎

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

VMD主要分為變分問題的構(gòu)建和求解，即：

對式(1)引入平衡參數(shù)α和拉格朗日乘數(shù)λ使該變分問題無約束化。因此增廣拉格朗日量可表示為：

用交替乘子法可得式(1)解對應式(2)的鞍點。再預先確定k，初始化參數(shù)并通過式(3)、式(4)對和ωk進行迭代更新。

更新和ωk后，再以式(5)對λ進行更新。

當滿足迭代精度ε時停止，即：

式(6)中，參數(shù)k保證了分解模式數(shù)的適當性和準確性；參數(shù)α與信號的重構(gòu)精度相關，二者的選擇對VMD 算法的分解效果尤為重要。當k和α過大時，容易造成模態(tài)混疊，反之會造成有用信息的缺失。因此本文引入麻雀搜索算法對兩者進行優(yōu)化。

1.2 麻雀搜索算法

SSA算法是Xue等[12]于2020年提出的一種新的優(yōu)化算法，可以概括為尋找-跟隨-預警的抽象模型，它模擬了麻雀的覓食過程以獲取待優(yōu)化問題的解。

設M是麻雀覓食的優(yōu)化搜索空間集，存在N只麻雀個體，第x只麻雀處在M集的位置表示為Sx=[sx1,…,sxd,sxM]，x=1,2,3,…,N，sxd表示麻雀x在M集中居d維的位置，尋找者位置更新表達式為：

其中：t為當前迭代數(shù)；K為最大迭代數(shù)；β是屬于區(qū)間(0,1]的均勻隨機數(shù)；P是服從N(0,1)分布的隨機數(shù)；J為1×d的單位矩陣；R2∈[0,1]為預警值，KZ∈[0.5,1]為安全值。而跟隨者位置更新可表達為：

式中：swtd為麻雀群進行t次迭代時處在d維最差的位置；反之sct+1d為t+1次迭代時最好的位置；當x＞n/2，表示適應度較低，需擴大搜索范圍；當x≤n/2時，表示適應度較高，可在sc位置周圍隨機覓食。預警麻雀位置迭代為：

式中：δ為服從N(0,1)的隨機數(shù)；V為[-1,1]的隨機數(shù)，表示麻雀移動的方向，同時控制步長；e是避免分母為零而設的極小值；hx為處于位置x時麻雀的適應度值，hw為當前麻雀的最差適應度值，hg則為最優(yōu)值。通常，預警麻雀個數(shù)占總種群的15%，為兼顧優(yōu)化準確性和計算效率，本文設置種群數(shù)和最大迭代數(shù)[N,M]=[30,20][13]。

在對VMD 的k、α進行優(yōu)化時，需考慮SSA 算法中一關鍵點，即適應度函數(shù)值的構(gòu)建。本文選取包絡熵為麻雀優(yōu)化算法的適應度函數(shù)值，包絡熵[14]可以很好評價信號的稀疏性，反映所研究信號分解情況的概率分布特性。

1.3 平方包絡譜峭度

在故障發(fā)生時，信號中的瞬時能量變化主要受故障沖擊和噪聲脈沖的影響。對于任意的模態(tài)分解分量，假設SE()為平方包絡信號，其方差的變化可很好表現(xiàn)故障信號的瞬時能量波動，表達式為：

式(10)中E(·)為數(shù)學期望，那么平方包絡信號的峭度可表示為：

由式(11)可知，當d()增大時，峭度也會隨之增大，而故障沖擊和噪聲脈沖的峭度值都偏大，當故障信號瞬態(tài)沖擊循環(huán)頻率過高時，其有效值會明顯增大，但信號的瞬時能量變化范圍反而會減小，d()降低，信號的峭度會減小，導致采用傳統(tǒng)峭度指標容易誤判最優(yōu)模態(tài)。由于低信噪比條件下故障特征所具備的包絡譜結(jié)構(gòu)易被干擾或淹沒，使其值與其他分量值相差不大，不利于篩選。因此，為了凸顯故障循環(huán)沖擊成分的有效性，選擇了平方包絡譜峭度作為篩選指標。

綜上分析，當分量信號中的瞬態(tài)沖擊循環(huán)頻率比較高時，信號的有效值會增加，d([n])較小，分量信號具有較小的峭度值，但平方包絡譜峭度卻能呈現(xiàn)較大值，此時信號的瞬時能量變化最大，克服了傳統(tǒng)峭度對單周期瞬態(tài)沖擊敏感但對多周期沖擊響應存在不足的缺陷，有利于檢驗故障信號中的循環(huán)沖擊特性[15]。

1.4 K-SVD算法

K-SVD算法分為稀疏編碼和字典更新兩部分。給定D∈Rm×K，當K大于m時，稱過完備字典。首先初始化并固定字典D，其中字典D常取DCT 字典。對信號矩陣Y=[y1,y2,y3,…,yn]∈Rm×n，稀疏編碼過程可表示為：

式中：δ為逼近誤差閾值(‖·‖q為lq范數(shù))，X=[x1,x2,…,xn]∈RK×n為待計算的稀疏系數(shù)矩陣。信號矩陣Y通常采用Hankel矩陣，其結(jié)構(gòu)表示為：

式中：H為轉(zhuǎn)換算子，m=N-n+1。

當?shù)玫较禂?shù)矩陣后，特征矩陣可重構(gòu)為：

式中：dk、xk分別表示D的第k列和X的第k行。可以看作是分解后的第k個分量。與原Hankel矩陣具有相同的維數(shù)。可以被視為退化的Hankel矩陣，因此H記為H-1的逆運算為：

式(18)中：P等于在給定條件下i和j的總組合數(shù)。當字典更新時，K-SVD利用殘差矩陣中的主成分按順序更新原子。殘差矩陣表示為：

將SVD應用于殘差矩陣，分離不同奇異值對應的分量。式中：Ekj=表示為第j個奇異分量；Δ=diag(σ1,σ2,…，σn)是奇異值的降序?qū)蔷仃?，σj表示第j個奇異值。uj、vj分別是U、V的第j列。在KSVD 中，dk被U的第一列更新，xk被σ1vT1取代，通過往復迭代，字典D中的所有原子都可以依次更新。

2 滾動軸承故障診斷流程

針對經(jīng)典K-SVD 算法在學習過程中易引入虛假原子導致信號稀疏不徹底以及VMD 中參數(shù)難以確定的問題，提出了基于麻雀算法優(yōu)化VMD參數(shù)與K-SVD的聯(lián)合診斷方法，其流程如圖1所示。

圖1 基于SSA-VMD聯(lián)合K-SVD的故障診斷流程圖

具體步驟如下：

步驟1：初始化麻雀算法種群數(shù)、最大迭代參數(shù)[N,M]、VMD 參數(shù)[k,α]優(yōu)化范圍。如滿足迭代條件，轉(zhuǎn)入下一步；否則更新各麻雀位置，繼續(xù)尋優(yōu)；

步驟2：將優(yōu)化參數(shù)組合[k,α]代入VMD 進行分解，得到k個本征模態(tài)分量IMFs；

步驟3：計算各本征模態(tài)的平方包絡譜峭度，選擇平方包絡峭度值最大的IMF分量；

步驟4：選擇步驟3優(yōu)選的IMF分量相空間構(gòu)建Hankel矩陣。初始化DCT字典，設置最大迭代次數(shù)L和編碼閾值δ，對Hankel 矩陣信號進行字典學習，迭代過完備字典和稀疏編碼系數(shù)，重構(gòu)信號并包絡解調(diào)。

3 仿真信號分析

根據(jù)軸承內(nèi)圈故障特點構(gòu)建仿真信號[16]：

設置內(nèi)圈故障仿真信號xi(t)幅值A0=1，采樣頻率fs=12 000 Hz，共振頻率fn=3000 Hz，仿真模型的衰減系數(shù)B=1000，故障特征頻率fi=120 Hz，轉(zhuǎn)頻fr=20 Hz，τi為服從μ=0、δ2=0.5%×fr的正態(tài)分布隨機滑動系數(shù)，同時加入信噪比為-15 dB 的高斯白噪聲n(t)，分析的信號長度為1×4 096。

內(nèi)圈故障仿真信號時域波形如圖2所示，可以看到?jīng)_擊成分被高斯白噪聲完全淹沒，難以獲取到?jīng)_擊規(guī)律和有效信息。其對應的內(nèi)圈故障仿真信號包絡譜如圖3所示，從中雖然能看到特征頻率，但其周圍存在嚴重的噪聲干擾譜線，且轉(zhuǎn)頻和其他特征倍頻也難以發(fā)現(xiàn)，無法有效進行故障特征識別。

圖2 仿真信號原始時域波形

圖3 仿真信號原始包絡譜

基于本文方法，首先初始化麻雀算法參數(shù)，其中k取[2,10]，平衡因子α取[200,3000]進行麻雀尋優(yōu)[17]。SSA迭代的適應度函數(shù)包絡熵變化趨勢如圖4所示，18次迭代后出現(xiàn)了最小包絡熵值3.622 7，此時通過優(yōu)選得到[k,α]=[6,2214]。將該參數(shù)組代入VMD 對仿真信號進行分解，如圖5所示，各IMF分量時域波形差異性不明顯。

圖4 仿真信號適應度函數(shù)迭代曲線

圖5 仿真信號VMD分解結(jié)果

為篩選出包含沖擊信息最多的模態(tài)分量，選擇了平方包絡譜峭度并與常見指標對比，歸一化結(jié)果如表1所示，其中，樣本熵指向模態(tài)1，歐式距離和峭度指向模態(tài)6。但模態(tài)1和模態(tài)6為虛假模態(tài)，無法給出故障特征頻率?；谄椒桨j譜峭度最大指標確定了最優(yōu)模態(tài)4，其頻譜中心頻率與所仿真信號的共振頻率3 000 Hz 相同，可以看出平方包絡譜峭度相較于傳統(tǒng)峭度指標具有更好的性能。

表1 仿真信號IMF分量遴選指標值

選取模態(tài)4分量構(gòu)建Hankel矩陣，初始化字典，設置字典維數(shù)m和原子個數(shù)n分別為140 和200，最大迭代次數(shù)L=10，編碼閾值δ=1.9[18]。利用KSVD 算法進行字典學習，最優(yōu)稀疏波形及其包絡譜如圖6和圖7所示，從包絡譜中可以明顯觀察到仿真信號的轉(zhuǎn)頻20 Hz、內(nèi)圈故障頻率120 Hz及其2倍特征頻率240 Hz。

圖6 仿真最優(yōu)模態(tài)信號K-SVD稀疏時域波形

圖7 本文所提方法仿真信號包絡譜

4 實驗信號分析

實驗數(shù)據(jù)來源于美國辛辛那提大學IMS中心[19]。實驗軸承為ZA2115 雙列滾子軸承，每列滾子數(shù)為16，滾子直徑為8.4 mm，節(jié)徑為71.5 mm，接觸角為15.17°，實驗時轉(zhuǎn)速為2 000 r/min，采樣頻率為20 kHz，采樣間隔為10 min，共984組，每組信號長度為1×20 480。實驗結(jié)束后發(fā)現(xiàn)轉(zhuǎn)軸上的軸承1外圈出現(xiàn)損壞，依據(jù)故障軸承參數(shù)，計算得該軸承外圈故障頻率fo=236 Hz，本文選取實驗數(shù)據(jù)中的第400組振動信號（分析的信號長度為1×5 120）進行驗證分析。原始信號的時域波形如圖8所示，可以看到時域波形故障沖擊幅值小、周期性不明顯，且伴隨大量背景噪聲。其包絡譜如圖9所示，特征頻率230.6 Hz幾乎被淹沒，難以進行故障識別。

圖8 實驗信號原始時域波形

圖9 實驗信號原始包絡譜

采用SSA 算法對原始信號進行迭代搜索，迭代優(yōu)化過程如圖10所示，可看到迭代至第9 次及其之后，最小包絡熵值穩(wěn)定在3.628 5，即得到最終優(yōu)化參數(shù)[k,α]=[5,316]，以此參數(shù)進行VMD信號分解，最終得到5個分解模態(tài)，其時域波形如圖11所示。

圖10 實驗信號適應度函數(shù)迭代曲線

圖11 實驗信號VMD分解結(jié)果

為定量選擇最優(yōu)模態(tài)，計算模態(tài)的篩選指標平方包絡譜峭度，進行歸一化處理并與常見指標對比，得到的模態(tài)篩選指標如表2所示。可以看到樣本熵指向模態(tài)1，歐式距離指向模態(tài)3，經(jīng)過處理發(fā)現(xiàn)模態(tài)1、3為虛假模態(tài)。而峭度和本文所提指標均指向模態(tài)5，以此易篩選得出模態(tài)5為最優(yōu)模態(tài)分量。

表2 實驗信號IMF分量遴選指標值

為了更好說明平方包絡譜峭度較峭度的優(yōu)越性，給出峭度與本文所提指標對比圖如圖12所示，雖然兩指標均指向模態(tài)5，但本文指標所判定的最優(yōu)分量在整體均值線的上方，而峭度指標波動較大，除最佳模態(tài)分量5外同樣存在其他分量在整體均線上方的情況，以此說明了本文所提指標具有更好的穩(wěn)定性且相較于其他常見指標具有更好的篩選能力。

圖12 峭度與本文所提指標對比圖

基于上述VMD預處理結(jié)果，以模態(tài)5分量相空間構(gòu)建Hankel 矩陣。初始化字典，取字典維數(shù)、原子個數(shù)分別為m=100,n=130，并設置最大迭代次數(shù)、編碼閾值分別為L=10，δ=0.1。利用K-SVD學習字典對所構(gòu)造的Hankel 矩陣信號進行稀疏編碼和字典學習，獲得已恢復至時間序列的稀疏重構(gòu)信號時域波形如圖13所示，可以看出經(jīng)字典學習的稀疏信號沖擊特征和周期性明顯，對其進行包絡分析如圖14所示，易觀察到明顯的故障頻率峰值230.6 Hz，證明所提方法提取軸承故障特征的有效性和泛化性。

圖13 實驗最優(yōu)模態(tài)信號基于K-SVD稀疏結(jié)果

圖14 經(jīng)所提方法處理的實驗信號包絡譜

為了說明結(jié)合VMD算法的必要性，現(xiàn)與設置相同參數(shù)的經(jīng)典K-SVD 算法進行對比。原始信號經(jīng)經(jīng)典K-SVD 處理后的包絡譜如圖15所示，由圖15可明顯看出，故障特征頻率的幅值小、包絡譜線峰值不突出，且在低頻段分布有大量的混淆頻率，難以識別故障。由此看出采用經(jīng)典K-SVD 方法提取微弱故障特征存在困難，本文所提方法可有效避免此局限性。

圖15 經(jīng)經(jīng)典K-SVD處理后實驗信號包絡譜

5 結(jié)語

本文針對VMD中模態(tài)分解層數(shù)k和平衡因子α難以選擇的問題，提出了以SSA 算法進行迭代尋優(yōu)方法。同時，針對經(jīng)VMD分解后的最優(yōu)模態(tài)難以選擇問題，引入平方包絡譜峭度遴選最優(yōu)模態(tài)分量，通過對最優(yōu)模態(tài)的K-SVD 字典學習和包絡檢波捕捉低信噪比條件下的軸承故障特征信息。仿真和試驗結(jié)果表明平方包絡譜峭度指標具有良好的適用性，所提采用SSA優(yōu)化VMD聯(lián)合K-SVD的診斷方法能夠很好地在低信噪比環(huán)境中提取滾動軸承故障，具有良好的泛化性和一定的工程實際意義。