趙錦劍，楊光永，2，周安然，項(xiàng)敏敏

(1．云南民族大學(xué)電氣信息工程學(xué)院，云南昆明　650504；2．華南理工大學(xué)機(jī)械工程學(xué)院，廣東廣州　510640)

0　引言

機(jī)械系統(tǒng)振動(dòng)信號處理及故障診斷方法包括自適應(yīng)濾波[1]、高階譜分析[2]以及獨(dú)立分量分析(ICA)[3]等。自適應(yīng)濾波適用于較復(fù)雜的機(jī)械振動(dòng)信號處理，它具有很強(qiáng)的自學(xué)習(xí)和自跟蹤能力，能夠通過系統(tǒng)辨識來對某一未知的動(dòng)態(tài)系統(tǒng)進(jìn)行逼近。高階譜是分析非平穩(wěn)信號和非高斯信號的有力工具，可以定量地描述信號中與故障密切聯(lián)系的非線性相位耦合。Yang[4]等研究了利用雙譜分析定量地描述信號中與機(jī)械故障密切聯(lián)系的二次相位耦合的方法，進(jìn)而對轉(zhuǎn)子故障進(jìn)行診斷。但常規(guī)的雙譜分析只以單通道信息為研究對象，不能完全地反映旋轉(zhuǎn)系統(tǒng)的非線性特征，從而影響了故障診斷的準(zhǔn)確性。ICA是盲源信號分離方法之一，能夠在無正交限制下抽取信號的統(tǒng)計(jì)獨(dú)立分量，并且對微弱機(jī)械振動(dòng)信號的特征提取有較好的效果。如文獻(xiàn)[5]對齒輪箱的振動(dòng)信號進(jìn)行分離，使微弱故障信息明顯增強(qiáng)。但I(xiàn)CA存在收斂速度較慢以及誤差累積等不足。

旋轉(zhuǎn)機(jī)械系統(tǒng)穩(wěn)定運(yùn)行過程產(chǎn)生的振動(dòng)響應(yīng)信號可以由AR模型來描述，常用的AR模型參數(shù)估計(jì)方法有最小均方(LMS)算法和遞推最小二乘(RLS)算法。LMS算法的特點(diǎn)是運(yùn)算簡單，無需計(jì)算相關(guān)函數(shù)，但收斂速度較慢；RLS算法雖然具有較快的收斂速度，但其得到的數(shù)據(jù)自適應(yīng)濾波器是對每一組輸入數(shù)據(jù)而言的，未能體現(xiàn)信息的完整性[6]。文中提出了一種基于Kalman濾波原理的旋轉(zhuǎn)機(jī)械振動(dòng)信號處理方法，并將其應(yīng)用于旋轉(zhuǎn)機(jī)械系統(tǒng)故障診斷。

1　Kalman濾波原理

Kalman濾波器是一種高效率的遞歸濾波器，它能夠從一系列的不完全及包含噪聲的測量中，估計(jì)動(dòng)態(tài)系統(tǒng)的狀態(tài)。其算法簡述如下[7]：

1．1狀態(tài)方程和觀測方程

狀態(tài)方程：

X(n+1)=F(n+1，n)X(n)+W(n)

(1)

觀測方程：

Y(n)=C(n)X(n)+V(n)

(2)

式中：X(n+1)為系統(tǒng)在離散時(shí)刻n+1的狀態(tài)向量；Y(n)為時(shí)刻n的觀測向量；矩陣F(n+1，n)和C(n)分別為狀態(tài)轉(zhuǎn)移矩陣和觀測矩陣；向量W(n)和V(n)為互不相關(guān)的零均值高斯白噪聲，它們的相關(guān)矩陣分別為Q1和Q2(n)。

1．2Kalman遞推公式

根據(jù)Kalman濾波算法原理，其遞推公式如下：

狀態(tài)向量預(yù)報(bào)方程：

(3)

狀態(tài)預(yù)報(bào)相關(guān)陣：

P(n+1，n)=F(n+1，n)P(n)FT(n+1，n)+Q1(n)

(4)

Kalman增益：

G(n)=P(n，n-1)CT(n)[C(n)P(n，n-1)CT(n)+Q2(n)]-1

(5)

新息過程：

(6)

一步預(yù)報(bào)：

(7)

狀態(tài)濾波誤差相關(guān)陣：

P(n)=P(n，n-1)-G(n)C(n)P(n，n-1)

(8)

Kalman濾波預(yù)估計(jì)就是用前面兩個(gè)時(shí)間更新方程(式(3)～式(4))獲得先驗(yàn)估計(jì)，然后通過后面4個(gè)狀態(tài)更新方程(式(5)～式(8))對先驗(yàn)估計(jì)矯正獲得最優(yōu)估計(jì)。根據(jù)這6個(gè)遞推公式，便可以連續(xù)地求出各個(gè)時(shí)刻的濾波值以及相應(yīng)的估計(jì)誤差。

2　旋轉(zhuǎn)機(jī)械振動(dòng)信號Kalman濾波模型的建立

Kalman濾波最初是為了解決工程控制問題而提出來的，但它的思想完全可用于對旋轉(zhuǎn)機(jī)械振動(dòng)信號進(jìn)行濾波處理。文中根據(jù)旋轉(zhuǎn)機(jī)械振動(dòng)信號的特點(diǎn)，利用Kalman濾波算法對振動(dòng)響應(yīng)信號的AR模型參數(shù)進(jìn)行估計(jì)，然后利用AR模型建立穩(wěn)定機(jī)械系統(tǒng)振動(dòng)響應(yīng)信號對應(yīng)的狀態(tài)空間方程，并結(jié)合Kalman濾波遞推方程組對旋轉(zhuǎn)機(jī)械振動(dòng)信號進(jìn)行降噪濾波處理。

2．1旋轉(zhuǎn)機(jī)械振動(dòng)響應(yīng)AR模型參數(shù)估計(jì)

通過數(shù)據(jù)采集系統(tǒng)得到的旋轉(zhuǎn)機(jī)械系統(tǒng)振動(dòng)信號數(shù)據(jù)屬于時(shí)間序列數(shù)據(jù)，因此可以建立旋轉(zhuǎn)機(jī)械系統(tǒng)振動(dòng)響應(yīng)信號的AR模型并進(jìn)行狀態(tài)預(yù)測。

如果把AR模型看成一步線性最小方差預(yù)報(bào)，已知的輸入向量由x(n-1)，x(n-2)，…，x(n-p)構(gòu)成，將這些輸入向量和p個(gè)參數(shù)ak進(jìn)行加權(quán)運(yùn)算，便可求出下一步預(yù)測向量x(n)的估計(jì)，希望輸出為x(n)。其方程式可以表示如下：

(9)

(10)

式中ε0(n)為期望輸出x(n)的最優(yōu)估計(jì)誤差。

AR模型的階數(shù)p一般事先是未知的，可以根據(jù)信息論準(zhǔn)則和最終預(yù)測誤差準(zhǔn)則來確定。

利用自適應(yīng)Kalman濾波算法對AR(p)參數(shù)進(jìn)行在線估計(jì)時(shí)，考慮到最優(yōu)參數(shù)是隨時(shí)間變化的，因此把參數(shù)估計(jì)過程看作非平穩(wěn)過程，引入過程噪聲ε1(n)，則有：

A(n)=A(n-1)+ε1(n)

(11)

根據(jù)Kalman濾波原理，由式(10)和式(11)可得Kalman濾波器的狀態(tài)方程和觀測方程：

狀態(tài)方程：

A(n)=A(n-1)+w(n)

(12)

觀測方程：

x(n)=X(n)A(n)+v(n)

(13)

式中w(n)和v(n)是相互獨(dú)立的零均值平穩(wěn)隨機(jī)過程。

利用Kalman濾波算法獲得旋轉(zhuǎn)機(jī)械振動(dòng)響應(yīng)信號的p階AR模型的參數(shù)a1，a2，…，ap和最優(yōu)估計(jì)誤差ε0(n)，便可建立式(9)所示的AR模型。

2．2基于AR模型的狀態(tài)空間方程的建立

設(shè)穩(wěn)定旋轉(zhuǎn)機(jī)械系統(tǒng)振動(dòng)響應(yīng)信號的時(shí)間序列數(shù)據(jù)為x(t)，t=1，2，…，n，利用前面所述的Kalman算法獲得穩(wěn)定旋轉(zhuǎn)機(jī)械系統(tǒng)振動(dòng)響應(yīng)信號的p階AR模型的參數(shù)a1，a2，…，ap以及e(t)，建立AR模型方程。則有：

(14)

式(14)可以改寫為：

(15)

若取狀態(tài)向量X(t)=[x(t)，x(t-1)，x(t-2)，……，x(t-p+1)]T，則可得到如式(1)形式的狀態(tài)方程：

X(t)=FX(t-1)+W(t)

(16)

若觀測數(shù)據(jù)為y(t)，則得到如式(2)的觀測方程：

y(t)=CX(t)+v(t)

(17)

3　基于Kalman濾波算法的故障診斷原理

4　仿真實(shí)驗(yàn)

根據(jù)旋轉(zhuǎn)機(jī)械振動(dòng)特點(diǎn)，構(gòu)造仿真信號：

x(t)=20[sin(100πt)+sin(180πt+π/6)]t

=[0，T]

(18)

y(t)=x(t)+n(t)

(19)

式中：x(t)為旋轉(zhuǎn)機(jī)械系統(tǒng)穩(wěn)定運(yùn)行時(shí)的振動(dòng)響應(yīng)原始信號；y(t)為疊加了噪聲的觀測信號；n(t)為干擾噪聲；T為原始信號的取樣長度。

4．1利用Kalman算法對旋轉(zhuǎn)機(jī)械振動(dòng)信號進(jìn)行濾波

對x(t)進(jìn)行仿真取樣，采樣頻率為1 kHz，取樣200個(gè)點(diǎn)，即T=0.2 s．當(dāng)p=6時(shí)，AR(6)模型波形與原始波形x(t)可以達(dá)到很好的擬合度，如圖1所示。此時(shí)AR(6)模型參數(shù)為：a1=-0．140 0，a2=-0．431 9，a3=1．629 4，a4=-0．786 6，a5=-1．986 9，a6=2．652 7。AR模型參數(shù)的收斂速度與RLS算法以及LMS算法的進(jìn)行比較，結(jié)果如圖2所示?？煽闯?，AR參數(shù)具有很好的收斂速度且優(yōu)于RLS算法和LMS算法。

圖1　原始信號波形與AR(6)模型波形

圖2　3種算法AR(6)模型參數(shù)的收斂速度

利用所求的AR(6)模型構(gòu)建狀態(tài)空間方程，結(jié)合Kalman濾波遞推方程組對觀測信號y(t)進(jìn)行濾波。由于機(jī)械結(jié)構(gòu)的復(fù)雜性，往往機(jī)械振動(dòng)信號中疊加了多種類型的噪聲，因此在文中給原始信號x(t)同時(shí)疊加均值為0、方差為15的高斯白噪聲，λ=15的泊松噪聲，以及σ=15的瑞利噪聲進(jìn)行濾波仿真實(shí)驗(yàn)，結(jié)果如圖3所示?？煽闯鯧alman濾波算法能夠有效地濾除觀測信號中的噪聲，在剛開始時(shí)估計(jì)值與實(shí)際值相差較大，但隨著Kalman迭代過程的進(jìn)行，估計(jì)值與實(shí)際值將不斷縮小，最終兩者的差值將趨于系統(tǒng)狀態(tài)噪聲值。

(a)疊加3種噪聲后波形

(b)原始波形和估計(jì)波形

4．2利用Kalman濾波算法對旋轉(zhuǎn)機(jī)械系統(tǒng)進(jìn)行故障診斷

對x(t)進(jìn)行仿真取樣，采樣頻率為1 kHz，取樣1 000個(gè)點(diǎn)，即T=1 s；設(shè)高斯白噪聲n(t)的均值為0、方差為5。在x(t)的1 000個(gè)取樣點(diǎn)中包含微弱突變，圖4為對應(yīng)的Kalman濾波估計(jì)誤差，可以明顯的看到大約在第850個(gè)取樣點(diǎn)處信號發(fā)生了突變。由此可判斷在第850個(gè)取樣點(diǎn)，設(shè)備運(yùn)行異常。

圖4　帶突變點(diǎn)原始信號的Kalman濾波估計(jì)誤差

5　結(jié)束語

文中為旋轉(zhuǎn)機(jī)械系統(tǒng)振動(dòng)信號的處理及故障診斷提供了一種有效可行的方法。該方法能夠有效地濾除旋轉(zhuǎn)機(jī)械振動(dòng)信號中的噪聲，為后續(xù)的機(jī)械振動(dòng)響應(yīng)信號的分析處理和運(yùn)行狀態(tài)的特征提取提供有效的數(shù)據(jù)保障，而且能夠有效地檢測出機(jī)械振動(dòng)信號中的微弱突變信息，對機(jī)械故障進(jìn)行診斷。同時(shí)，其收斂速度方面優(yōu)于RLS算法和LMS算法。

參考文獻(xiàn)：

[1]HAYKIN S．Adaptive filter theory．4th ed．Englewood cliffs：Prentice Hall，2002：436-498．

[2]NIKIAS C L，PETROPULU A P．Higher Order Spectral Analysis：A Nonlinear Signal Processing Framework．Englewood cliffs：Prentice Hall，1993．

[3]HYVARINEN A，KARHUNEN J，OJA E．Independent Component Analysis．New York：Wiley & Sons，2001．

[4]YANG J T，XU J W．Higher order spectral analysis in fault diagnosis of rotors．Chinese Journal of Mechanical Engineering，2001，14(1)：40-44．

[5]黃晉英，畢世華，潘宏俠，等．獨(dú)立分量分析在齒輪箱故障診斷中的應(yīng)用．振動(dòng)、測試與診斷，2008，28(2)：126-130．

[6]VASEGHI S V．Advanced Digital Signal Processing and Noise Reduction．4th ed．Chichester：John Wiley&Sons Ltd，2008：191-224．

[7]EUBANK R L．A Kalman Filter Primer．Boca rat on：CRC Press，2006．

[8]李曼，司頡，張鋒軍．礦井主通風(fēng)機(jī)在線監(jiān)測與故障診斷系統(tǒng)．儀表技術(shù)與傳感器，2013(1)：62-64．