劉前進(jìn)，高丙朋，宋振軍，王維慶

(新疆大學(xué) 電氣工程學(xué)院，烏魯木齊 830047)

軸承早期故障產(chǎn)生的沖擊十分微弱，易被噪聲掩蓋[1]，機(jī)械運(yùn)行時(shí)的轉(zhuǎn)速變化也隱藏一些軸承的故障特征信息；因此，對(duì)變轉(zhuǎn)速工況下軸承故障信號(hào)進(jìn)行預(yù)處理并提取故障特征，是變轉(zhuǎn)速軸承故障診斷的關(guān)鍵[2]。

近年來(lái)，針對(duì)滾動(dòng)軸承故障診斷問(wèn)題， 相關(guān)學(xué)進(jìn)行了大量的研究。從信號(hào)去噪的角度，文獻(xiàn)[3]提出基于峭度的經(jīng)驗(yàn)?zāi)B(tài)分解，利用峭度對(duì)沖擊成分的敏感特性提取軸承微弱故障特征，但峭度對(duì)瞬時(shí)突變特征提取能力有限，無(wú)法很好的診斷軸承動(dòng)態(tài)故障；文獻(xiàn)[4]提出基于變分模態(tài)分解(Variational Mode Decomposition，VMD)與多尺度熵相結(jié)合的軸承故障診斷方法，成功提取了軸承故障特征；文獻(xiàn)[5]借助包絡(luò)熵衡量軸承故障沖擊的稀疏程度，對(duì)變分模態(tài)分解的參數(shù)進(jìn)行優(yōu)化，提高了軸承故障特征的豐富性：上述方法普遍選取恒轉(zhuǎn)速軸承進(jìn)行研究，而旋轉(zhuǎn)機(jī)械實(shí)際運(yùn)行過(guò)程中的轉(zhuǎn)速存在波動(dòng)，軸承轉(zhuǎn)速波動(dòng)與故障沖擊響應(yīng)一樣都會(huì)造成特征模式的變化，但故障特征能量與轉(zhuǎn)速波動(dòng)相比，微小的故障特征被覆蓋，隱藏在與軸角相關(guān)的周期中[6]。階次跟蹤能將非平穩(wěn)時(shí)域信號(hào)轉(zhuǎn)化成循環(huán)平穩(wěn)的角域信號(hào)，為變轉(zhuǎn)速工況下的軸承故障診斷提供了研究方向[7]。

核互信息(Kernel-Based Mutual Information， KEMI)用來(lái)表示不同隨機(jī)向量之間的獨(dú)立程度[8]。文獻(xiàn)[9]提出基于特征選擇的核互信息，從筆跡樣本中預(yù)測(cè)作家的性別，能有效減少數(shù)據(jù)冗余和沖突。文獻(xiàn)[10]對(duì)哺乳動(dòng)物細(xì)胞環(huán)境的重建研究表明，核互信息算法對(duì)網(wǎng)絡(luò)重構(gòu)具有很強(qiáng)的復(fù)原能力。

卷積神經(jīng)網(wǎng)絡(luò)具有很強(qiáng)的自學(xué)習(xí)能力，在軸承故障診斷中有著廣泛的應(yīng)用[11-12]，但其模型搭建需要龐大的數(shù)據(jù)作為支撐。SVM是一種經(jīng)典分類(lèi)算法，利用核函數(shù)將非線性樣本映射到高維特征空間后將高維空間問(wèn)題轉(zhuǎn)化為二次規(guī)劃問(wèn)題，通過(guò)凸優(yōu)化實(shí)現(xiàn)樣本分類(lèi)，在處理非線性小樣本數(shù)據(jù)時(shí)具有獨(dú)特的優(yōu)勢(shì)[13]。

熵是描述粒子混亂程度的概念，也常用來(lái)對(duì)軸承故障進(jìn)行診斷，如近似熵、樣本熵、模糊熵等[14]。多尺度樣本熵作為衡量時(shí)間序列復(fù)雜性的特征指標(biāo)，時(shí)間序列越復(fù)雜，模式更新的機(jī)會(huì)就越大，對(duì)應(yīng)的熵值就會(huì)增大。對(duì)于變轉(zhuǎn)速工況下軸承產(chǎn)生的周期性沖擊，信號(hào)波動(dòng)會(huì)更加劇烈，因此采用多尺度樣本熵能有效區(qū)分不同故障信號(hào)的復(fù)雜程度。

針對(duì)旋轉(zhuǎn)機(jī)械運(yùn)行過(guò)程中軸承轉(zhuǎn)速變化，微弱故障特征不易提取等問(wèn)題，提出基于階次跟蹤(Order Tracking，OT)和參數(shù)優(yōu)化變分模態(tài)分解的特征提取方法，并使用模態(tài)分量的多尺度樣本熵作為特征向量，通過(guò)SVM分類(lèi)器對(duì)軸承故障進(jìn)行分類(lèi)。

1 階次跟蹤

階次跟蹤利用相位信息將振動(dòng)信號(hào)由時(shí)間采樣轉(zhuǎn)換成角度采樣，通過(guò)階次跟蹤反映與參考軸轉(zhuǎn)速相關(guān)的振動(dòng)特征。對(duì)于轉(zhuǎn)速變化的旋轉(zhuǎn)機(jī)械，其振動(dòng)和噪聲信號(hào)的頻率會(huì)隨著轉(zhuǎn)速變化而變化，為確定重采樣時(shí)間，假設(shè)參考軸的角加速度恒定，則軸角θ可表示為

θ=a0+a1t+a2t2，

(1)

式中：a0，a1，a2分別為初始位置、角速度和角加速度，可由轉(zhuǎn)速脈沖信號(hào)確定。

設(shè)定t1，t2，t3為3個(gè)連續(xù)的脈沖，對(duì)應(yīng)角度為0，δ，2δ，設(shè)置參考軸旋轉(zhuǎn)一圈(即δ=2π)，產(chǎn)生一個(gè)鍵脈沖，可以得到

(2)

根據(jù)已知的參數(shù)a0，a1，a2，參考軸角θ與時(shí)間t的關(guān)系為

(3)

假設(shè)恒定角度增量為Δθ，角度采樣方程可定義為

(4)

式中：Tn為第n個(gè)采樣點(diǎn)的采樣時(shí)間標(biāo)記。

階次跟蹤能將非平穩(wěn)時(shí)域信號(hào)轉(zhuǎn)換成平穩(wěn)的角域信號(hào)。因此，對(duì)變轉(zhuǎn)速工況下軸承故障信號(hào)進(jìn)行階次跟蹤處理，能有效消除因轉(zhuǎn)速變化帶來(lái)的不利影響。

2 變分模態(tài)分解及其參數(shù)優(yōu)化

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

變分模態(tài)分解的基本結(jié)構(gòu)是一個(gè)變分問(wèn)題，將信號(hào)x分解為K個(gè)離散的模態(tài)uk并確定每個(gè)模態(tài)的中心頻率ωk和帶寬，通過(guò)迭代搜索變分模型中的最優(yōu)解使每個(gè)模態(tài)的估計(jì)帶寬之和最小[15]。

變分模態(tài)分解方法將各模態(tài)分量定義為調(diào)幅-調(diào)頻(AM-FM)信號(hào)，即

uk(t)=Ak(t)cosφk(t)；k=1,2,…,K，

(5)

每個(gè)模態(tài)分量存在一個(gè)中心頻率ωk，將所有模態(tài)分量相加等于原信號(hào)作為約束條件，則約束變分模型可描述為

(6)

式中：uk為分解后的K個(gè)模態(tài)分量，uk={u1,…,uK}；ωk為各模態(tài)分量的中心頻率，ωk={ω1,…,ωK}。

在求解約束變分模型時(shí)，引入懲罰因子α和拉格朗日算子λ(t)，將約束變分問(wèn)題轉(zhuǎn)換為非約束問(wèn)題，即

(7)

(8)

借助帕賽瓦爾-傅里葉等距變換求解后，可以得到該優(yōu)化問(wèn)題的最優(yōu)解，即

(9)

(10)

(11)

(12)

式中：γ為拉格朗日算子系數(shù)。

變分模態(tài)分解方法運(yùn)算步驟如下：

1)初始化參數(shù)

λ1，α，n；

2)n=n+1；

4)根據(jù)(12)式更新λ；

5)設(shè)置精度ε，若滿足精度條件

(13)

結(jié)束循環(huán)，否則跳回步驟2繼續(xù)循環(huán)。

2.2 基于變異麻雀算法的變分模態(tài)分解參數(shù)優(yōu)化

參數(shù)α和K的選取會(huì)影響變分模態(tài)分解的效果[5]，故采用變異麻雀算法對(duì)變分模態(tài)分解參數(shù)進(jìn)行優(yōu)化。

2.2.1 麻雀算法

麻雀算法是根據(jù)麻雀的覓食和反捕食行為提出的群智能優(yōu)化算法，具有調(diào)整參數(shù)少，收斂速度快的優(yōu)點(diǎn)，但其全局搜索能力較弱，易陷入局部最優(yōu)[16]；因此，基于遺傳變異算法的思想對(duì)其進(jìn)行改進(jìn)，構(gòu)造變異麻雀算法。

麻雀種群分為發(fā)現(xiàn)者和加入者，警戒者從整個(gè)種群中隨機(jī)選取。在D維搜索空間中，麻雀種群為X，種群中麻雀?jìng)€(gè)數(shù)為n，即X=[x1,x2,…,xi],i=1,2,…,n，則種群中第i麻雀的搜索位置在D維空間中可表示為XD=[xi,1,xi,2,…,xi,d]。

1)發(fā)現(xiàn)者位置更新公式為

(14)

式中：t為當(dāng)前迭代次數(shù)，R2為預(yù)警值，R2∈[0,1]；RST為安全值，RST∈[0.5,1]；Q為服從正態(tài)分布的隨機(jī)數(shù)。

2)加入者位置更新公式為

(15)

式中：wi,j為當(dāng)前種群搜索出的最差位置；bi,j為當(dāng)前種群麻雀的最優(yōu)位置。

3)警戒者是從整個(gè)麻雀種群中隨機(jī)選取20%，其位置更新公式為

(16)

式中：β為服從正態(tài)分布的隨機(jī)數(shù)；κ為均勻隨機(jī)數(shù)，κ∈[-1,1]；fi為當(dāng)前麻雀?jìng)€(gè)體的適應(yīng)度；fg為全局最優(yōu)位置的適應(yīng)度。

2.2.2 核互信息函數(shù)

在變分模態(tài)分解參數(shù)優(yōu)化時(shí)，選用核互信息[17]作為變分模態(tài)分解效果的評(píng)判標(biāo)準(zhǔn)。核互信息定義為2個(gè)變量之間的依賴(lài)關(guān)系，當(dāng)互信息非常大時(shí)，表明2個(gè)變量之間存在彈性依賴(lài)，當(dāng)不存在依賴(lài)時(shí)，其值為零?？杀硎緸?/p>

(17)

式中：p(fi)，p(li)分別為變量f和l的邊緣概率密度函數(shù)；p(fi,li)為f和l的聯(lián)合概率密度函數(shù)。在變分模態(tài)分解中，各模態(tài)分量之間的互信息應(yīng)為最小，否則會(huì)導(dǎo)致模態(tài)混疊，而所有模態(tài)分量與原始信號(hào)的互信息之和應(yīng)為最大，避免在分解過(guò)程中丟失有用信息。因此，提出各模態(tài)之間互信息的和與原始信號(hào)與模態(tài)之間互信息之比，可以表示為

(18)

基于核互信息函數(shù)可以保證變分模態(tài)分解后的模態(tài)具有原始信號(hào)的最少互信息和最大信息，從而提取出原始信號(hào)中微弱的故障特征。因此，將(18)式作為算法的適應(yīng)度函數(shù)，則變異麻雀算法的參數(shù)優(yōu)化流程如圖1所示。

圖1 基于變異麻雀算法的參數(shù)優(yōu)化流程

將麻雀種群分為5個(gè)小組進(jìn)行參數(shù)搜索，分別求出每組最小的核互信息作為局部最優(yōu)解；然后，將每組局部最優(yōu)解作為全局最小核互信息的一部分以尋找全局中最優(yōu)參數(shù)，即全局最小核互信息對(duì)應(yīng)的參數(shù)。

minIKEMI={minIkemi}+βNi。

(19)

在局部最小核互信息minIkemi的基礎(chǔ)上，加入迭代次數(shù)Ni，構(gòu)建全局適應(yīng)度函數(shù)，其中β是適應(yīng)度函數(shù)的量化因子[5]，選取1/1 000。

3 變轉(zhuǎn)速工況下軸承故障診斷流程

與近似熵相比，樣本熵不依賴(lài)數(shù)據(jù)長(zhǎng)度，具有更好的一致性，因此具有更強(qiáng)的抗干擾能力。廣泛應(yīng)用于病理分析和故障診斷研究，樣本熵的主要參數(shù)有重構(gòu)維數(shù)m和閾值r，m一般設(shè)置為1或2，r設(shè)置為0.10～0.25倍的時(shí)間序列標(biāo)準(zhǔn)差[18]，計(jì)算步驟如下：

1)對(duì)于一段數(shù)據(jù)長(zhǎng)度為N的時(shí)間序列{x(n)}={x(1),x(2),…,x(N)}，已知重構(gòu)維數(shù)為m，則能夠得到的m維重構(gòu)向量為

Xm(i)={x(i),x(i+1),…,x(i+m-1)}；

1≤i≤N-m+1。

(20)

2)定義d[Xm(i),Xm(j)]為Xm(i)與Xm(j)兩向量中對(duì)應(yīng)元素最大差值的絕對(duì)值，即

d[Xm(i),Xm(j)]=

max(|x(i+k)-x(j+k)|)；k=0,1,…,m-1。

(21)

(22)

(23)

(24)

(25)

5)Am(r)和Bm(r)是2個(gè)序列在相似容限下匹配到m+1個(gè)點(diǎn)和m個(gè)點(diǎn)的概率，則樣本熵ESE可表示為

(26)

樣本熵只是在單一尺度上對(duì)時(shí)間序列進(jìn)行復(fù)雜度估計(jì)，而變轉(zhuǎn)速工況下的軸承故障特征通常分布在多個(gè)尺度上，使用單尺度樣本熵極易丟失特征信息，因此選用多尺度樣本熵作為特征指標(biāo)。多尺度樣本熵是以樣本熵為基礎(chǔ)，將時(shí)間序列粗粒化再進(jìn)行樣本熵計(jì)算，計(jì)算過(guò)程如下。

在一段時(shí)間序列{x(i),i=1,2,…,N}中，以長(zhǎng)度為τ的窗口對(duì)原序列進(jìn)行粗?；幚?，其表達(dá)式為

(27)

則多尺度樣本熵EMSE可以表示為

(28)

式中：τ為尺度因子；j為粗?；蛄袀€(gè)數(shù)。

變轉(zhuǎn)速工況下軸承故障診斷流程如圖2所示，具體實(shí)現(xiàn)步驟如下：

圖2 變轉(zhuǎn)速工況下軸承故障診斷流程圖

1)獲取變轉(zhuǎn)速工況下健康狀態(tài)、內(nèi)圈故障、外圈故障的軸承振動(dòng)信號(hào)；

2)使用階次跟蹤方法對(duì)原始信號(hào)進(jìn)行角域重采樣，得到循環(huán)平穩(wěn)的角域信號(hào)；

3)以核互信息為適應(yīng)度函數(shù)，通過(guò)變異麻雀算法搜最優(yōu)參數(shù)[α0,K0]，設(shè)置參數(shù)tau=0，DC=0，INIT=0，tol=1×10-7。變異麻雀算法的參數(shù)見(jiàn)表1(lb，ub邊界分別為變分模態(tài)分解變量α，K的約束)；

表1 變異麻雀算法參數(shù)設(shè)計(jì)

4)采用最優(yōu)參數(shù)變分模態(tài)分解方法分解角域信號(hào)，并計(jì)算各模態(tài)分量的多尺度樣本熵；

5)將多尺度樣本熵作為特征向量輸入到SVM分類(lèi)器構(gòu)建SVM診斷模型，選用徑向基函數(shù)RBF作為核函數(shù)，采用網(wǎng)格交叉驗(yàn)證方法對(duì)懲罰參數(shù)c和核函數(shù)參數(shù)g進(jìn)行優(yōu)化；

6)將測(cè)試樣本的特征向量輸入SVM模型得到軸承故障診斷結(jié)果。

4 仿真分析

為驗(yàn)證該方法的可行性，在變轉(zhuǎn)速條件下進(jìn)行仿真試驗(yàn)，綜合考慮轉(zhuǎn)速變化、外界干擾以及噪聲的影響，建立軸承外圈故障仿真模型[19]如下

(29)

t0=mod(k/fs,1/fm)；k=0,1,…，2 047，

f(t)=4sin(4t)+20，

式中：x1(t)為故障沖擊信號(hào)；x2(t)為外界干擾信號(hào)；r(t)為信噪比-5 dB的高斯白噪聲；β為衰減系數(shù)，取100；f1，f2分別取200，30 Hz；f(t)為變轉(zhuǎn)速軸承的轉(zhuǎn)速頻率；fs為采樣頻率，取1 024 Hz；fm為調(diào)制頻率，取16 Hz。

軸承轉(zhuǎn)速頻率變化，故障沖擊信號(hào)和軸承仿真故障信號(hào)如圖3所示，對(duì)比分析可知在軸承轉(zhuǎn)速波動(dòng)和干擾噪聲的影響下，故障沖擊信號(hào)的幅值和頻率信息全部丟失。

(a)軸承轉(zhuǎn)速頻率

根據(jù)軸承轉(zhuǎn)頻對(duì)外圈故障仿真信號(hào)進(jìn)行階次跟蹤處理得到角域信號(hào)，其中變分模態(tài)分解最優(yōu)參數(shù)為[3 270，2]，結(jié)果如圖4所示，角域信號(hào)分解得到2個(gè)主要頻率：30 Hz為干擾信號(hào)頻率，200 Hz為沖擊信號(hào)頻率。說(shuō)明該方法能夠在變轉(zhuǎn)速工況下有效提取軸承外圈故障信號(hào)中的沖擊特征。

(a)角域信號(hào)分解波形

5 試驗(yàn)分析

為充分驗(yàn)證本文所提方法對(duì)變轉(zhuǎn)速工況下軸承故障診斷的可靠性，選用加拿大渥太華大學(xué)機(jī)械工程實(shí)驗(yàn)室的數(shù)據(jù)集進(jìn)行試驗(yàn)[20]。試驗(yàn)軸承型號(hào)為ER16K型球軸承，其中包括通過(guò)加速度傳感器和增量編碼器采集不同軸承(正常及內(nèi)、外圈故障軸承)的振動(dòng)信號(hào)和轉(zhuǎn)速脈沖信號(hào)，采樣頻率為200 kHz，轉(zhuǎn)速范圍12～29 r/s，試驗(yàn)軸承的基本參數(shù)見(jiàn)表2(fr為軸的轉(zhuǎn)頻)。

表2 試驗(yàn)軸承基本參數(shù)

模態(tài)分量的能量占比越大，表示該分量的特征信息越多，包含故障特征的概率就越大[21]。不同K值時(shí)對(duì)3種軸承信號(hào)進(jìn)行變分模態(tài)分解并計(jì)算前3個(gè)模態(tài)分量與原信號(hào)能量占比的平均值，結(jié)果見(jiàn)表3，前3個(gè)模態(tài)分量的能量占比都在96%以上，說(shuō)明剩余模態(tài)分量的特征信息較小，可以忽略；因此，選取變分模態(tài)分解前3個(gè)模態(tài)分量的多尺度樣本熵作為特征向量。

表3 不同K值時(shí)變分模態(tài)分解前3個(gè)模態(tài)分量的能量占比

選定軸承旋轉(zhuǎn)一圈作為一組基本樣本，3種軸承信號(hào)各選取180組樣本，其中訓(xùn)練集120組，測(cè)試集60組，使用變異麻雀算法搜索變分模態(tài)分解最優(yōu)參數(shù)組合[α0,K0]，結(jié)果見(jiàn)表4。使用最優(yōu)參數(shù)變分模態(tài)分解訓(xùn)練樣本，選取前3個(gè)模態(tài)分量的多尺度樣本熵作為特征向量，將訓(xùn)練集3×120個(gè)特征向量輸入SVM分類(lèi)器中訓(xùn)練得到預(yù)測(cè)模型，最后利用SVM分類(lèi)模型對(duì)測(cè)試集進(jìn)行軸承故障分類(lèi)。

表4 滾動(dòng)軸承在不同工況下的變分模態(tài)分解最優(yōu)參數(shù)

加速運(yùn)行過(guò)程中不同故障軸承振動(dòng)信號(hào)前3個(gè)模態(tài)分量的多尺度樣本熵分布如圖5所示，由圖可知多尺度樣本熵能很好地區(qū)分不同故障類(lèi)型。

圖5 各模態(tài)分量的多尺度樣本熵分類(lèi)

采用參數(shù)優(yōu)化變分模態(tài)分解和多尺度樣本熵方法能較為準(zhǔn)確地診斷軸承故障類(lèi)型，為突出其在變轉(zhuǎn)速工況下的優(yōu)越性，采用不同方法對(duì)加速運(yùn)行過(guò)程中的軸承振動(dòng)信號(hào)進(jìn)行處理，結(jié)果見(jiàn)表5。其中，傳統(tǒng)變分模態(tài)分解方法中的參數(shù)α=2 000，K=3；OT表示階次跟蹤處理；VMD-BLS表示基于包絡(luò)熵的變分模態(tài)分解方法；VMD-KEMI是基于核互信息的變分模態(tài)分解方法。由表5可知：階次跟蹤方法能明顯提高變轉(zhuǎn)速工況下軸承故障診斷的準(zhǔn)確率，變分模態(tài)分解參數(shù)優(yōu)化也有效提高了故障識(shí)別的準(zhǔn)確率；核互信息方法雖然與包絡(luò)熵方法的平均診斷準(zhǔn)確率一樣，但穩(wěn)定性與一致性方面較好，說(shuō)明本文所提方法具有很強(qiáng)的應(yīng)用價(jià)值。

表5 不同方法對(duì)加速運(yùn)行過(guò)程中軸承故障的診斷結(jié)果

為進(jìn)一步證明本文所提方法在變轉(zhuǎn)速工況下的可靠性，選用相同型號(hào)軸承采集加速、減速、加速再減速、減速再加速這4種工況下的變轉(zhuǎn)速軸承信號(hào)，其速度譜如圖6所示。同時(shí)采用傳統(tǒng)變分模態(tài)分解方法(α=2 000，K=3)進(jìn)行對(duì)比分析，如圖7所示：在先加速再減速狀態(tài)下，傳統(tǒng)變分模態(tài)分解的診斷準(zhǔn)確率略高于參數(shù)優(yōu)化變分模態(tài)分解，但在其余狀態(tài)下均是參數(shù)優(yōu)化變分模態(tài)分解的診斷效果更好，尤其在軸承減速運(yùn)行過(guò)程中，傳統(tǒng)變分模態(tài)分解的故障診斷準(zhǔn)確率明顯下降，這是由于隨著軸承轉(zhuǎn)速降低，故障沖擊能量減小，故障特征更難提取。

圖6 不同工況下軸承運(yùn)行速度譜

①正常狀態(tài)；②內(nèi)圈故障；③外圈故障

在變轉(zhuǎn)速滾動(dòng)軸承故障診斷中，與傳統(tǒng)變分模態(tài)分解方法相比，參數(shù)優(yōu)化變分模態(tài)分解方法無(wú)論是故障診斷的準(zhǔn)確率還是故障識(shí)別的穩(wěn)定性都略勝一籌。變分模態(tài)分解參數(shù)一般憑借經(jīng)驗(yàn)選取，若選取不當(dāng)，振動(dòng)信號(hào)的特征信息不能被完整提取，造成軸承故障診斷不準(zhǔn)確，而通過(guò)變異麻雀算法對(duì)變分模態(tài)分解參數(shù)優(yōu)化，有效提取出故障的特征信息，提高了軸承故障診斷的準(zhǔn)確率，同時(shí)證明了該方法在變轉(zhuǎn)速滾動(dòng)軸承故障中能有效提取微弱故障特征，實(shí)現(xiàn)軸承故障的動(dòng)態(tài)診斷。

6 結(jié)束語(yǔ)

提出了一種基于階次跟蹤和參數(shù)優(yōu)化變分模態(tài)分解的軸承故障診斷方法并應(yīng)用于變轉(zhuǎn)速工況中的軸承故障診斷。該方法能有效地提取變轉(zhuǎn)速軸承故障特征，提高軸承故障診斷的準(zhǔn)確率。但只是針對(duì)單一故障問(wèn)題進(jìn)行診斷，對(duì)于變轉(zhuǎn)速軸承的復(fù)合故障還需要更深入的研究。