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

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

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

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

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

1 理論基礎(chǔ)

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

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

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

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

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

當(dāng)滿足迭代精度ε時(shí)停止，即：

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

1.2 麻雀搜索算法

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

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

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

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

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

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

1.3 平方包絡(luò)譜峭度

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

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

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

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

1.4 K-SVD算法

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

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

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

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

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

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

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

2 滾動(dòng)軸承故障診斷流程

針對(duì)經(jīng)典K-SVD 算法在學(xué)習(xí)過(guò)程中易引入虛假原子導(dǎo)致信號(hào)稀疏不徹底以及VMD 中參數(shù)難以確定的問(wèn)題，提出了基于麻雀算法優(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 進(jìn)行分解，得到k個(gè)本征模態(tài)分量IMFs；

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

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

3 仿真信號(hào)分析

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

設(shè)置內(nèi)圈故障仿真信號(hào)xi(t)幅值A(chǔ)0=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)分布隨機(jī)滑動(dòng)系數(shù)，同時(shí)加入信噪比為-15 dB 的高斯白噪聲n(t)，分析的信號(hào)長(zhǎng)度為1×4 096。

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

圖2 仿真信號(hào)原始時(shí)域波形

圖3 仿真信號(hào)原始包絡(luò)譜

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

圖4 仿真信號(hào)適應(yīng)度函數(shù)迭代曲線

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

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

表1 仿真信號(hào)IMF分量遴選指標(biāo)值

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

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

圖7 本文所提方法仿真信號(hào)包絡(luò)譜

4 實(shí)驗(yàn)信號(hào)分析

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

圖8 實(shí)驗(yàn)信號(hào)原始時(shí)域波形

圖9 實(shí)驗(yàn)信號(hào)原始包絡(luò)譜

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

圖10 實(shí)驗(yàn)信號(hào)適應(yīng)度函數(shù)迭代曲線

圖11 實(shí)驗(yàn)信號(hào)VMD分解結(jié)果

為定量選擇最優(yōu)模態(tài)，計(jì)算模態(tài)的篩選指標(biāo)平方包絡(luò)譜峭度，進(jìn)行歸一化處理并與常見指標(biāo)對(duì)比，得到的模態(tài)篩選指標(biāo)如表2所示?？梢钥吹綐颖眷刂赶蚰B(tài)1，歐式距離指向模態(tài)3，經(jīng)過(guò)處理發(fā)現(xiàn)模態(tài)1、3為虛假模態(tài)。而峭度和本文所提指標(biāo)均指向模態(tài)5，以此易篩選得出模態(tài)5為最優(yōu)模態(tài)分量。

表2 實(shí)驗(yàn)信號(hào)IMF分量遴選指標(biāo)值

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

圖12 峭度與本文所提指標(biāo)對(duì)比圖

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

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

圖14 經(jīng)所提方法處理的實(shí)驗(yàn)信號(hào)包絡(luò)譜

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

圖15 經(jīng)經(jīng)典K-SVD處理后實(shí)驗(yàn)信號(hào)包絡(luò)譜

5 結(jié)語(yǔ)

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