基于Kalman濾波和譜峭度的滾動軸承故障診斷

劉志川，唐力偉，曹立軍

(軍械工程學院 火炮工程系，石家莊 050003)

滾動軸承廣泛應用于旋轉機械中，是常見的易損件。軸承故障會導致機器出現振動和噪聲，嚴重時會影響整臺設備的正常運行，因此軸承的狀態檢測和故障診斷具有重要意義。

共振解調技術是目前軸承故障診斷中普遍使用的方法，但是存在帶通濾波器參數難以確定的問題，傳統的小波和小波包分析技術均不能自適應選擇最優帶通，導致應用受限。譜峭度[1]的基本原理是通過計算頻域中每根譜線的峭度值，從而在頻域中定位瞬態沖擊。文獻[2]提出將譜峭度值作為短時Fourier窗口函數，通過譜峭度圖選擇帶通濾波器的參數。針對譜峭度圖計算耗時的問題，文獻[3]提出了快速譜峭度圖的概念。文獻[4]則指出在信噪比較高的情況下，譜峭度能夠較好地識別瞬態沖擊信號。軸承工作環境相對復雜，環境噪聲較大，而復雜工況對譜峭度計算值會產生影響，從而給最優帶通濾波的帶寬和中心頻率的計算帶來誤差，影響診斷效果，因此有必要對信號進行降噪預處理。

Kalman濾波是一種動態數據處理方法，屬于遞推式濾波方法[5]，可以克服傳統時頻分析方法不適用于處理非平穩信號的問題，對測量數據逐一處理，由遞推方程隨時給出某一時刻新的狀態估計值，并與該時刻后的數據一起做出下一狀態的最優估計，達到信號降噪目的。在此，針對軸承內、外圈裂紋故障，將Kalman濾波與快速譜峭度圖相結合，以實現軸承的故障診斷。

1 AR模型

1.1 參數估計

對于時間序列{xt}，t=1，2，…,N，AR模型為

xt=φ1xt-1+φ2xt-2+…+φnxt-n+at，

(1)

式中：φ2,…，φn及σa是AR模型[6]的參數;n為模型階次。參數估計即估計出這n+1個數值，由于

at=xt-φ1xt-1-φ2xt-2-…-φnxt-n，

(2)

(3)

1.2 模型階數

AR模型的階次n未知，需要在遞推過程中確定。在遞推過程中，可以給出低階到高階的每一組參數，模型的最小預測誤差功率P是遞減的，當Pn達到最小值或者不再變化時，此時的階次就是AR模型的正確階次。

定階是建立AR模型的關鍵過程，其中AIC準則(Akaike Information Criterion)應用較為廣泛。它由2項構成：第1項隨階數增大而變小，體現模型擬合好壞；第2項隨階數增大而變大，體現模型參數多少。AIC準則表示為

AIC(n)=Nln(Pn)+2n。

(4)

從n=0開始增加模型階數，AIC(n)先減小,后增大。當階數n達到某一值n0時，AIC(n0)最小，則n0即為AR模型的最合適階次。

2 Kalman濾波

Kalman濾波是工程應用較為成熟的一種遞推線性最小方差估計，能夠通過前一時刻的最優值與當前時刻的測量值一起給出當前時刻的最優估計，能夠對測量數據進行實時更新。通過對原始數據建立的AR模型可以構建Kalman濾波的5個基本方程。

設齒輪箱振動真實信號為x(t),t=1,2,3,…，N，觀測噪聲為R，觀測值為y(t)，則觀測方程可表示為

y(t)=Hx(t)+R,

(5)

其中H=[1 0 … 0]1×n。設系統狀態為s(t)=[x(t),x(t-1),…,x(t-n+1)]T，則狀態轉移方程可表示為

s(t)=As(t-1)+a(t)，

(6)

其中In-1為n-1維單位陣。由此可建立Kalman濾波的5個基本方程：

(1)假設系統現在狀態為t，根據上一時刻狀態預測現在的狀態為

s(t|t-1)=As(t-1|t-1)+a(t);

(7)

(2)s(t|t-1)時刻協方差更新為

P(t|t-1)=AP(t-1|t-1)A′+Q;

(8)

(3)更新Kalman增益為

K(t)=P(t|t-1)H′/[HP(t|t-1)H′+R];

(9)

(4)新觀測得到的狀態最優估計值為

s(t|t)=X(t|t-1)+K(t)[y(t)-

Hs(t|t-1)];

(10)

(5)s(t|t)時刻協方差更新為

P(t|t)=[I-K(t)H]P(t|t-1),

(11)

式中：P為協方差；Q為系統過程噪聲；R為測量噪聲；K(t)為系統增益；I為單位矩陣。

3 譜峭度

信號的Wold-Cramer分解在頻域的表達式為

(12)

式中：H(t,f)為t時刻沖擊響應h(t,f)的頻域表示，是一個復包絡函數。

定義Y(t)過程的4階譜累積量為

(13)

2n階的譜瞬時矩為

S2nY(f)E{|H(t,f)dX(f)|2n}/df。

(14)

則譜峭度定義為

譜峭度[7-8]是關于中心頻率和帶寬的函數，當帶寬無限小時，譜峭度為0，當帶寬過大時，譜峭度無法反映頻帶范圍中的瞬態沖擊現象?？焖僮V峭度圖的原理是找到一個中心頻率fs和帶寬Bw，使峭度達到最大值，這時對應的中心頻率和帶寬即最優濾波器的2個參數。快速譜度圖是平面上的二維圖，圖中顏色的深淺代表不同中心頻率和帶寬組合下的峭度值，顏色最深處峭度值最大，對應獲得最優濾波器的中心頻率和帶寬。

4 基于Kalman濾波和譜峭度的故障診斷方法

軸承工作環境常常有較大噪聲，為降低振動信號噪聲，提高故障特征，先采用Kalman濾波技術對振動信號進行降噪預處理，然后利用快速譜峭度圖確定最優濾波器的中心頻率和帶寬。應用能量算子解調法解決傳統Hilbert包絡解調過程中存在的信號兩端調制和計算量大的問題，然后對包絡進行FFT變換得到包絡譜，將包絡譜中峰值較大的頻率成分與軸承故障特征頻率相比較，實現診斷，其流程如圖1所示。

圖1 軸承故障診斷方法流程圖

5 試驗驗證

軸承內、外圈故障實測信號來自二級減速器試驗平臺，試驗軸承為6206型深溝球軸承，鋼球個數Z=9，鋼球直徑Dw=9.5 mm，球組節圓直徑Dpw=46.5 mm。采用線切割分別在軸承的內圈和外圈上加工寬度0.5 mm、深度1 mm的細小裂紋作為故障。采集信號時，振動加速度傳感器B&K4508固定在故障軸承所在的軸承座上，采樣頻率為10 kHz；轉速由JN338型轉速轉矩傳感器測量。

5.1 內圈故障

輸入軸轉速為610.20 r/min，故障頻率f=27.55 Hz，中間軸轉動頻率fr=5.08 Hz，采樣時間為2.5 s。

首先通過原始信號建立系統AR模型，利用AIC準則計算得到的AIC曲線如圖2所示，從曲線中可以看出，階數取50時，AIC有極小值，因此確定模型階數為50。

圖2 AR模型的AIC曲線

Kalman濾波過程中，過程噪聲和觀測噪聲是2個重要參數[9]。觀測噪聲使用小波變換降噪的方法提取，對振動信號進行小波多尺度分解，得到小波系數，設置合適的閾值，大于閾值的小波系數置零，對保留的小波系數進行逆變換即為噪聲信號。過程噪聲在濾波過程中由大量試驗數據確定。軸承內圈故障信號濾波前、后的時域圖如圖3所示，從圖中可以明顯看出，濾波后信號幅值減小，沖擊更加明顯。

圖3 軸承內圈故障信號濾波前、后時域圖

對濾波后的信號進行譜峭度計算，繪制快速譜峭度圖，如圖4所示。

圖4 濾波后內圈故障信號的快速譜峭度圖

由圖4可知，在中心頻率為4 791.67 Hz、帶寬為416.67 Hz的頻帶內，峭度達到最大值35(在3.6層內)，其他頻帶信噪比較低。以此中心頻率和帶寬為帶通濾波器的最優參數對信號進行帶通濾波，并對帶通濾波后信號進行能量算子包絡解調，結果如圖5所示。圖中可以看出與內圈裂紋故障特征頻率接近的27.48 Hz及其倍頻等，而且故障軸承所在軸的轉動頻率(5.09 Hz)非常明顯，從結果可以判斷是內圈故障。

圖5 帶通濾波后內圈故障信號時頻圖

5.2 外圈故障

輸入軸轉速為1 488.5 r/min，故障頻率f=44.41 Hz，中間軸轉動頻率fr=12.4 Hz，采樣時間為1 s。對采集到的信號進行Kalman濾波及譜峭度計算，繪制的快速譜峭度圖如圖6所示。由圖可知，在中心頻率為4 375 Hz、帶寬為1 250 Hz的頻帶內峭度達到最大值13.6(在2層內)。對信號進行帶通濾波和能量算子包絡解調，結果如圖7所示，可以看出與外圈裂紋故障特征頻率接近的44.26 Hz及其倍頻等，而且故障軸承所在軸的轉動頻率(12.52 Hz)非常明顯，從結果可以判斷是外圈故障。

圖6 濾波后外圈故障信號的快速譜峭度圖

圖7 帶通濾波后外圈故障信號時頻圖

6 結束語

為解決傳統包絡分析方法中心頻率和帶寬參數難以確定的問題，提出了基于Kalman濾波和譜峭度圖的軸承故障診斷方法。首先建立系統AR模型，對原始信號進行Kalman濾波降噪預處理，提高信噪比，然后利用譜峭度自適應確定最優濾波器的中心頻率和帶寬，最后利用能量算子解調法對帶通濾波后的信號進行包絡解調分析，結合軸承故障特征頻率實現故障診斷。

對工程信號的分析處理驗證了本方法的可行性，而且在轉速相對較低的情況下也可以成功診斷出軸承內、外圈故障，為低速無載荷軸承故障診斷提供了新的研究方向。