粒子濾波在機械密封聲發射信號降噪中的應用*

(西南交通大學機械工程學院 四川成都 610031)

機械密封作為旋轉機械設備中轉軸密封的主要形式，在石油化工、泵用機械、船舶、航空航天等領域均得到了廣泛的應用。據統計，國外工業發達國家中95%左右的旋轉機械轉軸密封都采用了機械密封[1]。在機械密封裝置運行過程中，若能實時地獲取反映其運行狀態的有效信息，對系統運行狀態進行實時監測，對整個系統的性能和生產安全具有重要的意義。

機械密封裝置在運行過程中產生的聲發射(Acoustic Emission,AE)信號比振動、溫度、壓力等蘊含更為豐富的端面狀態信息，通過對機械密封裝置運行過程中的聲發射信號進行分析和識別，就可以實現對機械密封端面接觸情況以及整個密封裝置工作狀態的監測。實際應用中使用聲發射傳感器采集機械密封端面產生的聲發射信號時，往往會受到其他機械部件所產生的非線性、非平穩信號的干擾，因此對采集的聲發射信號還需要進行進一步地濾波降噪以便后續的分析處理。

傳統的聲發射信號降噪方法有經驗模態分解、小波閾值降噪、小波包降噪等。李曉暉等[2]提出一種基于Elman神經網絡的粒子濾波降噪算法，然而當神經網絡輸入層的神經元數目過大時，其粒子權值的遞推過程將變得很復雜。葛貞笛等[3]提出基于ARIMA模型的自校正卡爾曼濾波技術，但其不適用于非高斯系統且難以得到準確的模型參數。劉東瀛等[4]提出基于EMD與自相關原理的聲發射降噪算法，但其不能解決EMD的模態混疊問題且未保留高頻IMF分量中的有效信息，使信號產生損失。本文作者提出一種基于CEEMD與PF的聲發射信號降噪方法，該方法克服了傳統EMD算法的模態混疊現象和使用小波閾值降噪方法時閾值選擇困難的問題，且不受系統模型和噪聲模型的限制。仿真和實驗結果表明，該方法能很好地濾除聲發射信號中的背景噪聲，為后續的研究提供了良好的數據支撐。

1 CEEMD原理

經驗模態分解(Empirical Mode Decomposition，EMD)是HUANG等[5]提出的一種新的信號時頻處理方法，能夠很好地分析處理非線性非平穩信號。但EMD存在模態混疊的問題，使得分解的準確性降低。為了克服EMD的模態混疊現象，WU和HUANG[6]提出了總體經驗模態分解(Ensemble Empirical Mode Decomposition,EEMD)方法。該方法通過向原始信號中多次添加不同的白噪聲，然后對信號經EMD分解后的結果求取平均值，以此來抵消引入噪聲的影響。EEMD雖然解決了模態混疊問題，但添加噪聲的次數越多，相應的計算時間也越長。YEH和SHIEH[7]對EEMD進行了改進，提出了互補集合經驗模態分解(Complete Ensemble Empirical Mode Decomposition，CEEMD)方法，該方法通過向原始信號中添加n組正、負成對的輔助白噪聲，很好地消除了重構信號中存在的殘余輔助噪聲；同時通過選擇加入較低的噪聲集合次數，計算效率得到了一定提高。CEEMD的具體實現過程[8-9]如下：

(1)向原始信號x(t)中加入n組正、負成對的輔助白噪聲，由此得到的集合信號的個數為2n。輔助噪聲是均值為0，幅值系數k為常數的高斯白噪聲。當n取100～300時，k值取0.01～0.5倍信號的標準差。

(2)依次對集合中的每一個信號作EMD分解，每一個信號經過EMD分解后都可以得到一組IMF分量，記第i個信號的第j個IMF分量為VIMFij。

(3)將得到的2n組IMF分量求平均，就可以得到信號x(t)經CEEMD分解后的各階IMF分量：

VIMFj=12n∑2ni=1VIMFij

(1)

文中根據信號自相關函數的統計特性對CEEMD分解后的信噪分量進行區分[10]。如圖1所示，由于隨機噪聲在各個時刻的弱關聯性，其歸一化自相關函數值在零點處最大，在其他點處迅速衰減到很小值。而一般信號因為存在著強關聯性，雖然在零點處其歸一化自相關函數值也最大，但是在其他點處并沒有很快衰減到很小值，而是隨著時間差不斷地變化。

圖1 隨機噪聲與一般信號及其歸一化自相關函數

2 粒子濾波算法

粒子濾波(Particle Filter，PF)是通過在狀態空間中尋找一組離散的隨機樣本(即粒子集合)來近似系統隨機變量的概率密度函數p(xk|zk)，以樣本均值代替積分運算，從而獲得狀態的最小方差估計的過程。建立系統的數學模型是濾波的先決條件，但系統模型往往具有復雜的非線性和非高斯分布的特性。傳統的Kalman濾波受到高斯分布條件的限制，不能用于非高斯分布的場合。粒子濾波通過抽取大量的隨機狀態粒子近似表示后驗概率分布，能有效地應用和處理任何線性或非線性的系統模型，在解決非線性、非高斯問題上表現出極大的優越性。

粒子濾波算法的主要步驟有：重要性采樣、權值更新和重采樣[11]。

2.1 重要性采樣

重要性采樣基于蒙特卡洛模擬方法，它是從后驗概率分布中分別采集帶權重的粒子集，用粒子集近似表示后驗概率分布，將積分問題轉化為求和問題。但在實際應用中，難以直接從后驗概率分布中得到采樣粒子，為了解決這一問題，重要性采樣引入一個已知的且容易采樣的重要性概率密度函數q(xk|zk)，后驗概率分布就可以用從重要性概率密度函數中獲得的采樣粒子集的加權和來近似。

在估計后驗概率分布時貝葉斯重要性采樣需要利用所有的觀測數據，因此每當得到新的觀測數據，都需要重新計算每一個樣本的重要性權值，這就增大了計算量并占用很大的存儲空間。為了解決這一問題，提出了序貫重要性采樣(Sequential Importance Sampling,SIS),它在時間t+1時從重要性概率密度函數中采樣，不改變過去的狀態序列樣本集，并采用遞推的形式求得相應粒子的重要性權值，最終以粒子加權和來近似表示后驗概率密度。重要性概率密度函數可以寫成如下形式[12]：

q(X0∶k|Z1∶k)=q(X0∶k-1|Z1∶k-1)q(Xk|Z0∶k-1,Z1∶k)

(2)

2.2 權值更新

假設系統狀態符合一階馬爾科夫過程，且在給定的系統狀態下，各次的觀測量相互獨立，即

p(X0∶k)=p(X0)∏kj=1p(Xj|Xj-1)

(3)

p(Z1∶k|X0∶k)=∏kj=1p(Zj|Xj-1)

(4)

將式(2)—(4)代入權值遞推公式：

ωk=p(Z1∶k|X0∶k)p(X0∶k)q(X0∶k|Z1:k)=

p(Z1∶k|X0∶k)p(X0∶k)q(X0∶k-1|Z1:k-1)q(Xk|Z0∶k-1,Z1∶k)

(5)

又由式(5)有

ωk-1=p(Z1∶k-1|X0∶k-1)p(X0∶k-1)q(X0∶k-1|Z1:k-1)

(6)

由式(5)、(6)可得

ωk=

ωk-1p(Z1∶k|X0∶k)p(X0∶k)p(Z1∶k-1|X0∶k-1)p(X0∶k-1)q(Xk|X0∶k-1,Z1∶k)=

ωk-1p(Zk|Xk)p(Xk|Xk-1)q(Xk|X0∶k-1,Z1∶k)

(7)

在給定重要性參考分布q(Xk|X0∶k-1，Z1∶k)的條件下，式(7)給出了遞推計算粒子重要性權值的方法。如何選取重要性參考分布q(Xk|X0∶k-1，Z1∶k)是計算的關鍵問題，最優的選取方法是參考分布等于真實分布，即

q(Xk|X0∶k-1,Z1∶k)=p(Xk|X0∶k-1,Z1∶k)

(8)

但在實際情況下真實分布很難得到，因此常見的參考分布為先驗分布，即

q(Xk|X0∶k-1,Z1∶k)=p(Xk|Xk-1)

(9)

2.3 重采樣

在進行序貫重要性采樣時，經過多次迭代，只有少數粒子的權值較大，其余粒子的權值較小，可以忽略不計，為無效粒子。隨著狀態空間中無效采樣粒子數目的增加，大量的計算會浪費在對估計后驗概率分布幾乎不起作用的粒子更新上，使得粒子濾波的效率降低，這就是粒子權值退化問題。粒子權值的退化程度可以用有效粒子數Neff來衡量，Neff定義為

Neff=N/(1+var(ωik))

(10)

式中：N表示重采樣次數;var(ωik)為第i次重采樣時第k個粒子權值的方差值。

有效粒子數越小，表示粒子權值退化現象越嚴重。通常采用近似公式計算有效粒子數：

N︿eff=1/∑Ni=1(ωik)2

(11)

將重采樣引入到粒子濾波中，在一定程度上解決了SIS的粒子權值退化問題。粒子集合經過N次重采樣后，將之前的粒子集{X(i)k,ω(i)k}Ni=1更新為{X(i)k,1/N}Ni=1，即重采樣后每個粒子的權值相同，均為1/N。常用的重采樣方法有隨機重采樣、多項式重采樣、系統重采樣和殘差重采樣。

3 ARIMA模型與狀態空間模型

粒子濾波算法基于系統的狀態空間模型，模型包括狀態方程與觀測方程，即

Xk=f(Xk-1,Wk)

(12)

Zk=h(Xk,Vk)

(13)

式中：f(·)和h(·)分別表示狀態方程與觀測方程;Xk、Xk-1分別表示系統在k時刻和k-1時刻的狀態;Zk表示Xk的觀測值;Wk為系統過程噪聲;Vk為傳感器的測量噪聲，Wk和Vk相互獨立。

由于試驗條件的限制以及機械密封系統的復雜性，難以建立準確的聲發射信號狀態空間模型。因此文中通過對聲發射信號建立ARIMA模型，并利用小波分解重構的方法[10]提取背景噪聲，獲得其統計特性，再將ARIMA模型轉化為對應的狀態空間模型。

3.1 ARIMA模型原理

差分自回歸滑動平均模型(Autoregressive Integrated Moving Average Model,ARIMA)，又稱為Box-Jenkins模型，是20世紀70年代初由BOX和JENKINS提出的一種時間序列預測方法，其思想是將非平穩序列通過差分平穩化，然后采用ARMA模型分析處理差分后的序列。

差分自回歸滑動平均模型ARIMA(p，d，q)具有如下結構[13]：

{Φ(B)dxt=Θ(B)εt

E(εt)=0,var(εt)=σ2ε,E(εtεs)=0,s≠t

E(Xsεt)=0,?s

(14)

式中：B表示延遲算子，Bnxt=xt-n；為差分算子，d=(1-B)d表示差分運算；Φ(B)=1-Σpi=1φiBi為ARMA(p，q)模型的p階自回歸系數多項式;φ1,......,φp為自回歸系數，φp≠0；Θ(B)=1-Σqi=1θjBj為ARMA(p，q)模型的q階滑動平均系數多項式；θ1,......,θq為滑動平均系數，θq≠0；εt為服從N(0,σ2ε)的白噪聲序列。

3.2 ARIMA模型建模過程

(1)時間序列平穩性判斷。根據時間序列{xt}的時序圖、自相關函數(Autocorrelation Function，ACF)圖、偏相關函數(Partial Autocorrelation Function，PACF)圖，或者利用增廣迪基-福勒(Augmented Dickey-Fuller，ADF)檢驗方法對時間序列的平穩性進行判斷。如果序列為非平穩序列，則需要先對序列進行平穩化處理再進行模型識別。

(2)非平穩序列平穩化處理。如果時間序列{xt}是非平穩的，則需要利用差分運算對序列進行平穩化。對每次差分后得到的新序列都要進行平穩性判斷，如果仍不平穩，則需要進行多次差分直至序列平穩。設非平穩時間序列{xt}經過d階差分運算后得到的平穩時間序列為{yt}，則

yt=dxt=(1-B)dxt=∑di=0(-1)iCidxt-i

(15)

(3)白噪聲檢驗。在對平穩的時間序列建立模型之前，還要對序列進行白噪聲檢驗，以確定得到的平穩時間序列是否具有繼續研究的價值。白噪聲序列各序列值之間不存在任何的相關信息，沒有必要繼續進行分析研究。

(4)模型識別。如表1所示，根據平穩非白噪聲時間序列的ACF和PACF的截尾與拖尾性特征，確定要建立的模型類型并初步確定模型的p、q值，其中p和q分別為偏相關函數和自相關函數顯著不為0的最高階數。

(5)階數確定與參數估計。根據赤池信息量準則(Akaike Information Criterion，AIC)或貝葉斯信息準則(Bayesian Information Criterion，BIC)確定最優p、q值。采用極大似然準則估計模型的參數。AIC或BIC由下式計算得出：

VAIC=-2lnL+2k

(16)

VBIC=-2lnL+lnn·k

(17)

式中：k是參數個數；L是似然函數；n為樣本數量。

(6)模型檢驗。得到模型后，還需對模型殘差進行自相關檢驗。利用Ljung-Box檢驗Q統計量判斷模型的殘差是否為高斯白噪聲序列，只有當模型的殘差為高斯白噪聲序列時，建立的ARIMA(p，d，q)模型才認為是有效的。

表1 ARMA模型識別

3.3 模型的轉化

對得到的ARIMA(p，d，q)模型，需要將其轉化為如下形式的狀態空間模型：

[xt

xt-1

xt-p+1]=[φ1φ2......φp

0

E

0][xt-1

xt-2

xt-p]+

[1-θ1......-θq

0

O

0][εt

εt-1

εt-q]

(18)

zt=[10......0][xt

xt-1

xt-p+1]+μt

(19)

式中：φ1,......,φp、θ1,......,θp是ARIMA模型的各項系數；εt表示過程噪聲；μt表示通過小波變換提取的測量噪聲；E為單位陣；O為零矩陣。

4 仿真分析

選擇頻率為0.25 Hz與0.5 Hz的余弦信號疊加而成的信號作為原始信號，信號長度為N=5 000。在原始信號中添加均值為0，方差為5的高斯白噪聲后將其作為仿真信號，信噪比為-9.076 0。圖2所示為原始信號與加噪聲信號的時域波形圖。

圖2 仿真信號時域波形圖

圖3 不同方法的降噪結果

對仿真信號分別進行CEEMD小波閾值、標準粒子濾波、CEEPF降噪，降噪結果如圖3所示。為了評價各種方法降噪的效果，采用信噪比SNR、平均絕對誤差MAE、均方根誤差RMSE作為評價指標，可由式(20)、(21)、(22)計算得到。

VSNR=10lg(∑x2∑(z-x)2)

(20)

VMAE=1N∑Ni=1|z-x|

(21)

VRMSE=1N∑Ni=1(z-x)2

(22)

式中：x表示純凈原始信號;z表示降噪后的信號。

表2所示為仿真信號經3種方法降噪后各指標的數據對比。綜合以上的數據分析，可以看出CEEPF方法在聲發射信號降噪上較CEEMD小波閾值和標準粒子濾波算法有更好的效果。

表2 3種降噪方法評價指標對比

5 實驗驗證

為驗證CEEPF降噪方法的可行性與有效性，采用在四川日機密封件股份有限公司的機械密封試驗平臺上采集的一組連續型聲發射信號進行實驗分析，采樣頻率為10 kHz，時間為1 s，采樣點數為10 000，信號時域波形圖如圖4所示。

圖4 機械密封聲發射信號時域波形圖

對聲發射信號進行CEEMD分解，得到的各階IMF分量波形圖如圖5所示，其中n=200，k=0.1。

圖5 IMF分量時域波形圖

圖6所示為各階IMF分量的自相關函數圖，可以看出，前8階IMF的歸一化自相關函數值在零點處取得最大，在其他點處迅速衰減至很小值，但不為0。這是因為信號經CEEMD分解后白噪聲的統計特性遭到了破壞，各高頻IMF分量已不再是真正意義上的白噪聲。因此可以判定前8階IMF分量為噪聲分量，后7階IMF分量為信號分量。

圖6 IMF分量自相關函數

重構前8階IMF分量，并對重構后的信號進行ADF檢驗，判斷出信號是非平穩的。對信號進行一階差分后得到如圖7所示的差分信號，且ADF檢驗的結果為平穩，滿足ARMA模型的建模條件。差分信號的ACF與PACF圖如圖8所示，可以看出其PACF呈負指數形式衰減，在滯后時間為8 s之后落入置信區間內；ACF呈余弦振蕩形式衰減，在滯后時間為6 s之后落入置信區間內。信號的自相關函數與偏相關函數均表現出截尾性，根據ARMA模型的識別準則，對差分信號建立ARMA(p，q)模型，且初步確定p=8，q=6。

圖7 一階差分信號

圖8 一階差分信號的ACF與PACF

根據赤池信息量準則，在得到的48個備選模型中，最小的AIC值為-3.029 2×104，對應的p、q值分別為7和5，因此最終確定的最優模型為ARMA(7,5)。采用極大似然準則估計的ARMA模型參數值如表3所示。對模型殘差進行Ljung-Box檢驗，殘差為白噪聲，模型可用。

表3 一階差分信號ARMA模型參數表

根據得到的ARMA模型參數，可計算模型ARIMA(7,1,5)的表達式為

xt=-0.448 64xt-1+0.431 46xt-2-0.096 3xt-3+

0.287 516 4xt-4-0.211 404xt-5+0.222 696 7xt-6+

0.063 073 8xt-7-0.074 366 5xt-8+εt+0.311 835εt-1-

0.506 312εt-2+0.043 497 4εt-3-0.335 813εt-4-

0.513 207εt-5

(23)

其中εt～N(0,0.002 823 01)。

為了對信號進行粒子濾波分析，需要將信號的ARIMA模型轉化為對應的狀態空間模型。根據前文的分析，直接將ARIMA模型作為信號的狀態方程，εt作為狀態噪聲；利用小波分解重構的方法提取信號的背景噪聲，將其作為信號的測量噪聲，并和ARIMA模型一起作為信號的測量方程。提取的背景噪聲方差為0.003 0，選取粒子數N=1 000。差分信號經粒子濾波后的時域波形圖如圖9所示，將其與9至15階的IMF分量進行重構，就可以得到降噪后的信號。

圖9 粒子濾波結果

對原始聲發射信號分別進行CEEMD小波閾值、標準PF、CEEPF方法降噪，降噪效果如圖10所示。采用信號峰值、均值和標準差作為評價指標對3種方法的降噪效果進行定量分析，3種評價指標可由以下各式計算得到：

Xpeak=12)max(xi)-min(xi))

(24)

X=1N∑Ni=1xi

(25)

σx=∑Ni=1(xi-X)2

(26)

圖10 3種方法的降噪效果對比

如表4所示，機械密封聲發射信號經3種方法降噪之后，信號的最大瞬時幅值都得到了降低。降噪前后信號的均值幾乎保持不變，即機械密封聲發射信號的真實信號并未受到降噪帶來的影響。從信號的標準差來看，降噪后的標準差明顯要小于降噪之前，說明經過降噪處理聲發射信號的隨機噪聲得到了有效的抑制。CEEPF方法的標準差顯然小于CEEMD小波閾值和標準PF 2種方法，表明CEEPF方法能夠更好地濾除背景噪聲，達到不錯的降噪效果。

表4 不同降噪方法評價指標對比

6 結論

(1)利用相關系數原理和標準粒子濾波算法對CEEMD算法進行改進，提出了基于CEEMD與PF的機械密封聲發射信號降噪方法，該方法克服了EMD的模態混疊現象，彌補了小波閾值降噪方法的不足，能最大程度地避免有效信息的損失，使CEEMD方法能夠達到更好的降噪效果。

(2)相對于CEEMD小波閾值、標準PF降噪方法，基于CEEMD與PF的機械密封聲發射信號降噪方法能最大程度地濾除信號中的背景噪聲以及保留原始信號中的有效信息，為后續的特征提取、狀態識別奠定了良好的數據基礎。