陳兆峰，高偉鵬，朱丹宸

(1.海軍研究院，北京 100161；2.海軍士官學(xué)校，安徽 蚌埠 233012)

滾動軸承安裝在各類旋轉(zhuǎn)機械中，在工業(yè)生產(chǎn)中得到廣泛使用，并發(fā)揮著重要作用。然而，滾動軸承常工作在高溫、高壓和高轉(zhuǎn)速的環(huán)境下，極易產(chǎn)生故障，給機械設(shè)備的安全穩(wěn)定運行帶來較大的不確定性，因此急需對滾動軸承的故障特點、故障診斷方法進(jìn)行更深入的研究。目前，振動監(jiān)測方法是滾動軸承故障診斷中的常用方法，但也存在一些問題：1)設(shè)備振源多，導(dǎo)致滾動軸承故障信號常常淹沒在強背景噪聲中；2)為了不破壞設(shè)備結(jié)構(gòu)，測點往往只能選取在遠(yuǎn)離故障軸承的機殼外側(cè)，故障點與測試點之間存在復(fù)雜的傳遞路徑，進(jìn)一步加劇噪聲干擾。

事實上，滾動軸承故障特征提取所要解決的就是在強背景噪聲干擾下提取弱故障特征。為解決這個難題，專家學(xué)者們提出了很多有效的方法，如共振解調(diào)[1-2]、經(jīng)驗?zāi)Ｊ椒纸鈁3]、變分模態(tài)分解[4-5]、形態(tài)學(xué)濾波[6-7]等。上述方法都能夠?qū)L動軸承的故障狀態(tài)進(jìn)行判斷，但也存在一定問題，如強背景干擾下共振頻帶的準(zhǔn)確選取問題，分解參數(shù)設(shè)置對變分模態(tài)分解的影響問題，形態(tài)學(xué)算子以及結(jié)構(gòu)參數(shù)的確定問題等。考慮傳遞路徑的影響，為更有效地提取故障源信號，盲解卷積算法近些年得到了快速發(fā)展，常見的包括最小熵卷積(Minimum Entropy Deconvolution, MED)[8]、最大相關(guān)峭度解卷積(Maximum Correlated Kurtosis Deconvolution, MCKD)[9-10]、自適應(yīng)多點最優(yōu)最小熵解卷積(Multipoint Optimal Minimum Entropy Deconvolution Adjusted, MOMEDA)[11]等。這些方法都是以特征指標(biāo)最大為目標(biāo)進(jìn)行解卷積，去除無關(guān)成分的干擾，進(jìn)而獲取故障特征更突出的信號；但這些算法也受到眾多參數(shù)的制約，如MCKD中移位數(shù)M的選取，MOMEDA中解卷積周期的確定，以及廣泛存在的濾波器長度選取問題。針對這些問題，專家學(xué)者們對算法進(jìn)行了眾多改進(jìn)，提升了實際使用效果[12-13]。近來，最大二階循環(huán)平穩(wěn)盲解卷積(Maximum Second-Order Cyclostationarity Blind Deconvolution, CYCBD)算法[14]的提出進(jìn)一步增強了盲解卷積算法的有效性，它將最大二階循環(huán)平穩(wěn)指標(biāo)(Second-Order Cyclostationarity, ICS2)作為分析目標(biāo)，提取包含滾動軸承周期性故障沖擊特征的源信號；然而，基于CYCBD的特點，為保證分析效果，準(zhǔn)確預(yù)估濾波器長度和特征頻率值尤為重要，但如何更為準(zhǔn)確地選取參數(shù)仍然值得進(jìn)行深入探討[15-17]。

綜上所述，基于CYCBD的特征提取能力，本文提出一種增強最大二階循環(huán)平穩(wěn)盲解卷積(Enhanced CYCBD, ECYCBD)的方法以實現(xiàn)強背景噪聲干擾下的滾動軸承弱故障特征提取。ECYCBD以重加權(quán)峭度值作為適應(yīng)度函數(shù)，預(yù)先確定濾波器長度和特征頻率的搜索范圍；利用鯨魚優(yōu)化算法獲得最優(yōu)解，確保CYCBD的分析效果；利用Teager能量算子進(jìn)一步突出故障成分，實現(xiàn)滾動軸承故障特征的準(zhǔn)確提取。

1 理論方法

1.1 最大二階循環(huán)平穩(wěn)盲解卷積

假設(shè)觀測信號為x，目標(biāo)源信號為s0，則盲解卷積算法的計算過程可以表示為

s=x*h≈s0，

(1)

式中：s，h分別為通過解卷積得到的源信號和濾波器；*為卷積算子。為了方便計算，可以用矩陣形式進(jìn)行表示

s=Xh，

(2)

(3)

式中：L和N分別為s和h的長度。

基于此，最大二階循環(huán)平穩(wěn)指標(biāo)可表示為[14]

(4)

(5)

(6)

E=[e1,…,ek,…,eK]，

(7)

(8)

式中：H為矩陣的共軛轉(zhuǎn)置；RXWX，RXX和W分別為加權(quán)相關(guān)矩陣、相關(guān)矩陣和加權(quán)矩陣；diag為構(gòu)建對角矩陣；k為樣本數(shù)；Ts為故障周期的長度。

基于循環(huán)頻率與信號能量波動相關(guān)的特性，對于軸承、齒輪故障信號，循環(huán)頻率值與故障現(xiàn)象、頻率相關(guān)，循環(huán)頻率值可設(shè)定為k/Ts[14]。此時，向量ek任意一項的上標(biāo)中，k(N-1)/Ts可改寫為kfstN-1/Ts，其中tN-1為第N-1個數(shù)據(jù)對應(yīng)的時間，fs為采樣頻率，則fs/Ts即為滾動軸承的故障特征頻率。

求解下列方程

RXWXh=RXXhλ，

(9)

最大的特征值λ即為最大的ICS2值[14]。

由于算法運行過程中h的初始值并不確定，需要預(yù)先設(shè)定以獲得初始的加權(quán)矩陣W。受到初始化過程隨機性的影響，需要進(jìn)行如下迭代過程以獲得最大ICS2值：1)利用尤爾—沃克方程計算自回歸模型的系數(shù)，并作為白化濾波器對h進(jìn)行初始化；2)借助信號x和h計算得到矩陣W；3)通過(9)式計算得到最大特征值λ及其特征向量h；4)返回第2步，將特征向量h代入重新進(jìn)行計算直至收斂。

本文的收斂條件定義為：1)迭代次數(shù)達(dá)到100次；2)|λi-λi-1|/|λi-1|<0.001，其中λi，λi-1分別為本次和上次迭代中獲得的特征值。

上述計算過程的輸入信號和輸出信號均只有一個，為提高分析效果，CYCBD可用于多個信號輸入的情況，從而綜合利用多信號特征。假設(shè)有Q個觀測信號xq，通過Q個濾波器hq同時對相應(yīng)的觀測信號進(jìn)行解卷積計算，則最終的源信號s是多個信號的解卷積結(jié)果之和。

(10)

此時，整個算法的分析過程仍可利用上述步驟，通過(9)式獲得最大的ICS2值。此時獲得的濾波器h是多個與輸出相關(guān)的濾波器系數(shù)的組合[14]，即

h=[h1，…，hq，…，hQ]T，

(11)

第q個輸出對原信號的貢獻(xiàn)可表示為

sq=Xqhq，

(12)

矩陣Xq的構(gòu)造方法與(3)式中一致且與信號xq相關(guān)，最終的濾波信號為

(13)

1.2 鯨魚優(yōu)化算法

鯨魚優(yōu)化算法(Whale Optimization Algorithm,WOA)[18]模擬了鯨魚圍捕獵物的行為，主要包括鯨魚在覓食過程中的3種行為，分別是包圍獵物、搜尋獵物和氣泡網(wǎng)環(huán)形游動驅(qū)趕獵物并進(jìn)行位置更新，這幾種行為可通過數(shù)學(xué)表達(dá)式表示。

1)包圍獵物。第t+1代中備選解的位置X(t+1)可以通過(14)式進(jìn)行更新

(14)

(15)

式中：X*(t)為第t代中的最優(yōu)位置；a為動態(tài)變量，隨迭代次數(shù)的增加從2至0線性遞減；r為[0,1]范圍內(nèi)的隨機數(shù)。

2)搜尋獵物。搜索獵物的過程與包圍獵物的過程非常類似，不同的是，搜尋獵物的過程不再使用當(dāng)次迭代的最優(yōu)解X*(t)，而是隨機選取了一個備選位置Xrandom(t)

(16)

3)螺旋式位置更新。鯨魚在圍捕獵物的同時會噴出氣泡網(wǎng)用以驅(qū)趕獵物，同時不斷更新自身位置。使用氣泡網(wǎng)時，鯨魚位置的更新公式為

(17)

式中：D′為第t代的最優(yōu)位置與任意備選位置的距離；b為常數(shù)；l為[-1,1]范圍內(nèi)的隨機數(shù)。

鯨魚優(yōu)化算法的總過程如下：

1)初始化鯨魚種群位置Xi(i=1,2,…,N)，計算各鯨魚所在位置對應(yīng)的適應(yīng)度值，X*(t)即為最優(yōu)適應(yīng)度值對應(yīng)的位置。

2)設(shè)定參數(shù)a，A，C，l，p的值，本文計算時參數(shù)b設(shè)定為1。

3)若p<0.5，當(dāng)|A|<1時，利用(14)式對鯨魚當(dāng)前的位置進(jìn)行更新，否則利用(16)式對鯨魚當(dāng)前的位置進(jìn)行更新；若0.5≤p<1，則用(17)式更新鯨魚的當(dāng)前位置。

4)利用新位置更新a，A和C，并重新計算更新后鯨魚位置的適應(yīng)度值，重新確定最優(yōu)位置X*(t)。

5)重復(fù)步驟3和4直到達(dá)到最大迭代數(shù)。

1.3 重加權(quán)峭度值

本文使用一種重加權(quán)峭度指標(biāo)[19]解決峭度指標(biāo)容易受到少量沖擊干擾成分影響的問題，該指標(biāo)的計算過程如下：

1)將采集到的信號x共劃分n段，得到的子信號分別記為x1,x2,…,xn。

2)計算各子信號的峭度值k1,k2,…,kn；對n個峭度值進(jìn)行升序排序，得到峭度序列K。

3)計算每個峭度值所占比

(18)

4)重加權(quán)峭度值為

(19)

1.4 Teager能量算子

Teager能量算子本質(zhì)上是對信號動態(tài)成分總能量的估計，包括信號的瞬時值和微分形式，由于其計算簡便，得到了廣泛使用。

對于連續(xù)信號x(t)，可將其Teager能量算子表示為

(20)

對于離散時間信號x(n)，其Teager能量算子為

ψ[x(n)]=x2(n)-x(n-1)x(n+1)。

(21)

由于Teager能量算子可以對瞬態(tài)沖擊成分進(jìn)行增強，因此適用于突出滾動軸承的故障特征成分。

1.5 故障特征提取流程

基于上述理論，本文提出基于ECYCBD的滾動軸承弱故障特征提取方法，適用于強背景噪聲干擾下的滾動軸承故障特征提取問題，該方法的流程如圖1所示，具體過程為：

圖1 基于ECYCBD的滾動軸承弱故障特征提取方法流程

1)獲取滾動軸承故障信號并利用信號的最大幅值(取絕對值)進(jìn)行標(biāo)準(zhǔn)化處理(文中仿真信號和試驗信號均進(jìn)行了量綱一化處理)，計算該滾動軸承故障特征頻率的理論值f。

2)基于故障特征頻率的理論值確定濾波器長度和特征頻率值的搜索范圍，為保證CYCBD的分析效果和計算效率，本文將故障周期的長度記為L，進(jìn)而將濾波器長度的選取范圍確定為[0.75L,3L]，特征頻率的選取范圍確定為[f-5,f+5]。

3)以最大重加權(quán)峭度值為目標(biāo)，利用鯨魚優(yōu)化算法優(yōu)化選取CYCBD的參數(shù)。

4)借助最優(yōu)參數(shù)組合，獲得包含明顯故障特征的解卷積信號。

5)利用Teager能量算子和快速傅里葉變換得到該解卷積信號的Teager能量譜。

6)通過識別Teager能量譜中的特征判斷滾動軸承故障類型。

2 仿真分析

首先利用仿真信號對本文方法的有效性進(jìn)行驗證，仿真信號構(gòu)建時模擬了滾動軸承內(nèi)圈故障的特點。為模擬較為復(fù)雜的背景環(huán)境，該信號包含幅值調(diào)制、諧波信號和白噪聲成分，可以表示為

(22)

式中：Ai為幅值調(diào)制成分；B(t)為諧波成分；s(t)為自由衰減信號；n(t)為隨機白噪聲；fr為轉(zhuǎn)頻，取42 Hz；fm為諧波信號頻率，取100 Hz；fn為系統(tǒng)的固有頻率，取3 200 Hz；r為阻尼系數(shù)，取0.05；T=1/185 s，即該信號的故障特征頻率為185 Hz；ti為第i個周期內(nèi)由于滾動體滑移引起的延遲，ti=0.01T～0.02T。

仿真信號構(gòu)造時，為突出背景成分的干擾，設(shè)定信號的信噪比為-10 dB，同時信號的采樣頻率和時長分別為32 768 Hz和2 s。該仿真信號的時域波形和包絡(luò)譜如圖2所示，時域波形中未能識別出周期性的軸承故障沖擊成分，包絡(luò)譜中僅能勉強識別故障特征頻率及其2倍頻成分，譜圖中存在較強的干擾，說明對于該仿真信號，包絡(luò)分析效果不夠理想。本文的包絡(luò)譜均是對信號直接進(jìn)行希爾伯特變換和快速傅里葉變換獲得。

(a) 時域波形

(b) 包絡(luò)譜

下面利用本文提出的方法處理該仿真信號，在優(yōu)化過程中，若WOA的種群數(shù)量和迭代次數(shù)設(shè)定過大，將導(dǎo)致整個算法的計算效率較低，若參數(shù)設(shè)置過小，會導(dǎo)致最后的尋優(yōu)結(jié)果不理想。為兼顧算法的有效性和運行時間，通過多次嘗試，本文將WOA中的參數(shù)設(shè)定為：種群數(shù)量20，最大迭代次數(shù)30，其余參數(shù)都選取默認(rèn)值；濾波器長度和特征頻率的搜索范圍確定為132～531和180～190 Hz。分析得到的迭代收斂結(jié)果如圖3a所示，經(jīng)過16次迭代，適應(yīng)度函數(shù)取得最大值4.815，此時對應(yīng)的鯨魚位置為[525,185]。將此位置帶入CYCBD進(jìn)行計算，得到的最優(yōu)解如圖3b所示，該信號中的噪聲成分得到明顯去除，諧波成分的干擾得到有效抑制，且能夠體現(xiàn)出一些沖擊成分。該信號的Teager能量譜如圖3c所示， 圖中故障特征頻率185 Hz及其倍頻成分得到了準(zhǔn)確提取，無關(guān)頻率成分的干擾得到明顯抑制。通過上述特征可以判斷該軸承存在內(nèi)圈故障，表明本文提出的方法能夠解決強背景噪聲帶來的干擾，準(zhǔn)確有效提取滾動軸承弱故障特征。

(a) 迭代收斂曲線

(b) 最優(yōu)解卷積信號

(c) Teager能量譜

為進(jìn)一步突出本文方法的分析效果，利用CYCBD和MCKD進(jìn)行對比。

使用CYCBD對仿真信號進(jìn)行分析，對參數(shù)進(jìn)行預(yù)先確定，將濾波器長度和特征頻率值分別設(shè)定為312和185 Hz，分析結(jié)果如圖4所示。解卷積信號中噪聲干擾仍較為明顯，Teager能量譜提取出滾動軸承故障特征成分但幅值較低，尤其185和555 Hz處譜線幅值不突出，無關(guān)頻率成分仍產(chǎn)生較明顯干擾，分析效果欠佳。上述結(jié)果說明在分析過程中，隨機確定的CYCBD參數(shù)并不準(zhǔn)確。

(a) 解卷積信號

(b) Teager能量譜

利用經(jīng)典的MCKD方法分析該仿真信號，并計算Teager能量譜。MCKD的參數(shù)依據(jù)本文優(yōu)化過程得到的最優(yōu)解[525,185]進(jìn)行設(shè)置，其中濾波器長度取為525，沖擊信號的周期T=172。MCKD方法對仿真信號的分析結(jié)果如圖5所示，解卷積信號中故障沖擊成分難以得到識別且噪聲干擾較強，Teager能量譜雖然在特征頻率附近出現(xiàn)了幅值增大現(xiàn)象，但未能提取出明顯的故障特征頻率成分。

(a) 解卷積信號

上述對比結(jié)果進(jìn)一步突出本文算法在確定最優(yōu)參數(shù)，準(zhǔn)確提取故障特征上的有效性。

3 試驗驗證

為進(jìn)一步凸顯本文算法在工程實際中的價值，對某型試驗臺獲得的滾動軸承故障信號進(jìn)行分析，試驗臺結(jié)構(gòu)如圖6所示，可用于模擬滾動軸承在實際工作中的常見故障類型和運行工況。為模擬實際設(shè)備結(jié)構(gòu)中的復(fù)雜傳遞路徑，確定測點時不考慮靠近故障軸承的軸承座測點，而是選取加載軸承端(測點1)和支承軸承端(測點2)作為本次試驗的測點，測得的信號分別記為信號1和信號2，振動測量方向為徑向。此時，軸承故障信號在向外傳遞時需要經(jīng)過多個連接界面，還要受到主軸、正常軸承等結(jié)構(gòu)的干擾，信號較為復(fù)雜，給故障特征提取帶來一定難度。

圖6 試驗臺

試驗軸承型號為6010，內(nèi)、外徑分別為50，80 mm，球組節(jié)圓直徑Dpw為65 mm，球徑Dw為9 mm，球數(shù)Z為13，接觸角α為0°。利用激光切割法在內(nèi)圈溝道上割出一道寬0.2 mm、深0.2 mm的窄縫模擬內(nèi)圈故障。試驗過程中，試驗臺轉(zhuǎn)速為3 000 r/min(轉(zhuǎn)頻fr=50 Hz)，徑向和軸向加載分別為1和2 kN，采樣頻率為32 768 Hz，軸承內(nèi)圈故障特征頻率fi為370 Hz。

為更好提取故障源信號，本文在實測信號分析時采用多輸入單輸出的CYCBD方法，同時使用信號1和信號2作為CYCBD的輸入進(jìn)行解卷積計算。實測信號及其包絡(luò)譜如圖7所示： 由于復(fù)雜傳遞路徑及噪聲成分的干擾，時域波形中故障導(dǎo)致的沖擊成分未能得到識別；包絡(luò)譜中無關(guān)干擾頻率成分較為明顯，故障特征未能得到有效提取：上述結(jié)果說明僅利用包絡(luò)解調(diào)方法難以在強噪聲干擾下獲得理想效果。

(a) 信號1

(b) 信號2

利用本文方法分析該實測信號，結(jié)果如圖8所示。分析時，濾波器長度的搜索范圍為66～265，特征頻率的搜索范圍為365～375 Hz。通過WOA得到的迭代曲線如圖8a所示，經(jīng)過7次迭代可獲得最大適應(yīng)度值為4.742，此時對應(yīng)的最優(yōu)位置為[254，366]，此組參數(shù)也就是CYCBD計算時的最優(yōu)參數(shù)組合。最優(yōu)特征頻率值366 Hz與理論值不一致可能是由于軸承尺寸存在制造誤差以及試驗過程中轉(zhuǎn)速控制不準(zhǔn)確。將優(yōu)化結(jié)果代入CYCBD進(jìn)行計算，得到的最優(yōu)解如圖8b所示，能夠識別出周期性的沖擊成分，噪聲干擾得到明顯抑制。Teager能量譜如圖8c所示，圖中存在轉(zhuǎn)頻成分、故障特征頻率及其倍頻成分，圍繞故障特征頻率成分的調(diào)制邊頻帶都得到了準(zhǔn)確提取， 故障特征較為豐富。通過上述分析結(jié)果，能夠確定該軸承存在內(nèi)圈故障，驗證了本文方法在弱特征信號提取中的有效性。

(a) 迭代收斂曲線

(b) 最優(yōu)解卷積信號

(c) Teager能量譜

與仿真分析類似，同樣使用2種方法進(jìn)行對比分析。

將CYCBD分析時的濾波器長度和特征頻率值分別隨機設(shè)定為321和370 Hz，將信號1和2共同作為CYCBD的輸入， 分析結(jié)果如圖9所示。解

(a) 解卷積信號

(b) Teager能量譜

積信號中難以識別故障沖擊成分且噪聲干擾成分的抑制效果不明顯，Teager能量譜中也未能識別出任何有關(guān)故障的特征，說明隨機確定CYCBD的分析參數(shù)難以獲得較好的效果。

利用最優(yōu)參數(shù)組合[254,366]確定MCKD的參數(shù)，MCKD的濾波器長度和沖擊信號的周期分別取254和90，對2個信號進(jìn)行分析。信號1的MCKD分析結(jié)果如圖10所示，解卷積信號中未能識別出明顯的沖擊成分，且Teager能量譜僅能勉強確定故障特征頻率的2倍頻成分，故障特征并不豐富，噪聲干擾較為明顯，分析效果欠佳；信號2的MCKD分析結(jié)果如圖11所示，信號2的采集位置離故障軸承更遠(yuǎn)， 從其時域波形或Teager能量譜中均不能識別出有關(guān)故障的任何信息。

(a) 解卷積信號

(b) Teager能量譜

(a) 解卷積信號

(b) Teager能量譜

4 結(jié)束語

針對滾動軸承故障信號中噪聲干擾成分強，故障沖擊成分弱的問題，提出了基于ECYCBD的故障特征提取方法，其具有以下優(yōu)點： 1)利用增強的CYCBD方法可以抑制無關(guān)成分干擾，有效提取故障特征，借助Teager能量譜可以進(jìn)一步突出故障成分，提高解調(diào)效果；2)以最大重加權(quán)峭度值為目標(biāo)，借助WOA的參數(shù)尋優(yōu)，保證了CYCBD參數(shù)選取的準(zhǔn)確性，提高算法的分析效果。仿真和實測信號的對比分析結(jié)果凸顯了本文方法在噪聲干擾成分抑制和故障特征成分提取上的有效性。