基于增強最大二階循環平穩盲解卷積的滾動軸承弱故障特征提取

陳兆峰，高偉鵬，朱丹宸

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

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

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

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

1 理論方法

1.1 最大二階循環平穩盲解卷積

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

s=x*h≈s0，

(1)

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

s=Xh，

(2)

(3)

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

基于此，最大二階循環平穩指標可表示為[14]

(4)

(5)

(6)

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

(7)

(8)

式中：H為矩陣的共軛轉置；RXWX，RXX和W分別為加權相關矩陣、相關矩陣和加權矩陣；diag為構建對角矩陣；k為樣本數；Ts為故障周期的長度。

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

求解下列方程

RXWXh=RXXhλ，

(9)

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

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

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

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

(10)

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

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

(11)

第q個輸出對原信號的貢獻可表示為

sq=Xqhq，

(12)

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

(13)

1.2 鯨魚優化算法

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

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

(14)

(15)

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

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

(16)

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

(17)

式中：D′為第t代的最優位置與任意備選位置的距離；b為常數；l為[-1,1]范圍內的隨機數。

鯨魚優化算法的總過程如下：

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

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

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

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

5)重復步驟3和4直到達到最大迭代數。

1.3 重加權峭度值

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

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

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

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

(18)

4)重加權峭度值為

(19)

1.4 Teager能量算子

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

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

(20)

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

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

(21)

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

1.5 故障特征提取流程

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

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

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

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

3)以最大重加權峭度值為目標，利用鯨魚優化算法優化選取CYCBD的參數。

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

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

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

2 仿真分析

首先利用仿真信號對本文方法的有效性進行驗證，仿真信號構建時模擬了滾動軸承內圈故障的特點。為模擬較為復雜的背景環境，該信號包含幅值調制、諧波信號和白噪聲成分，可以表示為

(22)

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

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

(a) 時域波形

(b) 包絡譜

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

(a) 迭代收斂曲線

(b) 最優解卷積信號

(c) Teager能量譜

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

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

(a) 解卷積信號

(b) Teager能量譜

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

(a) 解卷積信號

上述對比結果進一步突出本文算法在確定最優參數，準確提取故障特征上的有效性。

3 試驗驗證

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

圖6 試驗臺

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

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

(a) 信號1

(b) 信號2

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

(a) 迭代收斂曲線

(b) 最優解卷積信號

(c) Teager能量譜

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

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

(a) 解卷積信號

(b) Teager能量譜

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

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

(a) 解卷積信號

(b) Teager能量譜

(a) 解卷積信號

(b) Teager能量譜

4 結束語

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