基于Rao－Blackwellised粒子濾波的巖石聲發(fā)射信號去噪

凌同華，劉偉宏，朱 亮，吳聯(lián)迎

在巖石聲發(fā)射的信號處理中，信號去噪是最基礎(chǔ)、也是最關(guān)鍵的一項工作。目前，信號去噪處理的方法主要有兩種：頻域和時間－空間域分析方法。假設(shè)噪聲和期望信號成分在頻域范圍內(nèi)相互不重疊，可以通過頻域方法濾除噪聲（如：高通濾波、帶通濾波及帶阻濾波等）。例如：白噪聲分布于整個頻域范圍內(nèi)，傳統(tǒng)的頻域方法已經(jīng)不能清除噪聲成分。相反，一種基于噪聲統(tǒng)計特性的狀態(tài)空間方法被用來去除噪聲成分。在時域上建立起來的常用濾波有：維納濾波、卡爾曼濾波［1－2］及擴展卡爾曼濾波等。維納濾波局限于處理線性平穩(wěn)靜態(tài)過程，并具有很好的效果?？柭鼮V波常被用于處理非靜態(tài)過程的信號。當需要確保濾波精度時，動態(tài)系統(tǒng)需滿足線性且噪聲服從高斯白噪聲分布。另外，擴展卡爾曼濾波已被用來處理非線性系統(tǒng)，但對預(yù)測狀態(tài)的估值問題仍需線性化處理，噪聲要求服從高斯分布。

近幾年，迅速發(fā)展起來的粒子濾波［3－5］（Particle Filter，簡稱為PF）是一種基于蒙特卡羅的非線性、非高斯系統(tǒng)狀態(tài)估計的濾波方法，完全突破了卡爾曼濾波理論框架，對系統(tǒng)的過程噪音和量測噪音不加任何限制。在PF基礎(chǔ)上，作者擬提出Rao－Blackwellised 粒 子 濾 波 （Rao－Blackwellised particle filter，簡稱為RBPF）來加強巖石聲發(fā)射信號采集質(zhì)量。

1 粒子濾波

從動態(tài)系統(tǒng)中獲取信息和估計系統(tǒng)的狀態(tài)，建立一個精確的模型是至關(guān)重要的。大多數(shù)物理現(xiàn)象服從隨機分布，其狀態(tài)可以用有限維向量描述，并通過向量差分方程來模擬。建立PF系統(tǒng)模型，至少需要兩個方程：①狀態(tài)方程，用于描述系統(tǒng)狀態(tài)隨時間演變的過程；②觀測方程，用于將系統(tǒng)在某時刻的輸出和系統(tǒng)的狀態(tài)聯(lián)系起來。

式中：Xk，Yk分別為動態(tài)系統(tǒng)在k時刻的狀態(tài)向量和觀測向量；f 為狀態(tài)函數(shù)；h為觀測函數(shù)；Wk－1，Uk分別為狀態(tài)過程和觀測過程的噪聲。

粒子濾波算法流程如圖1所示，其詳細過程見文獻［3，4］。

圖1 PF算法流程Fig.1 Flow chart of PF algorithm

2 Rao－Blackwellised粒子濾波

2.1 RBPF方法

由于聲發(fā)射信號的采集頻率比較高（通常在兆赫茲以上），粒子濾波需要大量的樣本點來組建后驗概率分布，使計算效率下降。尤其在高維狀態(tài)空間和實時監(jiān)控應(yīng)用中，計算效率下降較明顯。為實現(xiàn)粒子濾波在實時聲發(fā)射監(jiān)控中應(yīng)用，采取具有降維作用的 RBPF［6－9］。

RBPF將標準狀態(tài)空間向量xk分解成兩個子空間，其中一部分x1k可以通過卡爾曼濾波進行分析解計算；而另一部分x2k可通過粒子濾波計算，并使其計算的復(fù)雜程度大大降低。RBPF的目的就是用較少的粒子獲取同樣的估計精度，從而降低計算量，實現(xiàn)其成本的節(jié)約與推廣應(yīng)用?；趦蓚€子空間所構(gòu)成的條件概率，并對其進行推演：

式中：p（x21：k｜y1：k）為x21：k的后驗概率密度函數(shù)。

如果p（x11：k｜x21：k，y1：k）服從條件線性高斯分布，則可通過卡爾曼濾波實現(xiàn)其最優(yōu)估計。因此，由兩個子空間構(gòu)成的后驗概率密度函數(shù)通過運用粒子濾波估計時，僅需考慮與x21：k相關(guān)的一部分。RBPF可以看成是卡爾曼濾波與粒子濾波的結(jié)合體，是一種加強的粒子濾波。

2.2 RBPF實時處理特征

RBPF是一種狀態(tài)空間時域的濾波技術(shù)，每次輸入一個新的數(shù)據(jù)算法立即更新。在更新過程中，RBPF不依賴于全部已輸入的信號序列，僅需考慮前一時刻的數(shù)據(jù)。因此，在采集過程處理信號時，對采集數(shù)據(jù)長度和采樣頻率沒有特殊要求，此外，RBPF在計算過程所需存儲數(shù)據(jù)的空間將大大降低，其運算速度獲得較大的提高。斷鉛試驗常被用來模擬聲發(fā)射源，斷鉛試驗中采集的聲發(fā)射信號如圖2所示。從圖2中可以看出，濾波效果不受采集數(shù)據(jù)長度的影響。RBPF進行信號去噪處理時不受將來采集數(shù)據(jù)的干擾，僅與當前時刻及前一時刻所處的狀態(tài)有關(guān)。傳統(tǒng)的濾波器不具有這一特征，在濾波過程中需要完整的數(shù)據(jù)鏈。由于RBPF有實時在線處理功能，并能處理高維的非線性非高斯系統(tǒng)，而備受研究者的關(guān)注。

圖2 RBPF實時處理信號濾波的特征Fig.2 The real－time processing feature of RBPF based signal filtering

3 RBPF在聲發(fā)射信號去噪中應(yīng)用

3.1 狀態(tài)空間模型的建立

由于地震波信號與聲發(fā)射信號具有相似性，地震波模型可以成為合適的候選對象來代表結(jié)構(gòu)中的聲發(fā)射事件。其模型為：

式中：A1（t－t0）為聲發(fā)射源波形的振幅；t0為波形到達時間；ω為主導(dǎo)角頻率（ω＝2πf）。

在處理地震數(shù)據(jù)中，類似的方法被其他學(xué)者采用，地震波的振幅可用隨機過程X1（t）表示，同時地震波用隨機過程Xw（t）表示，則該過程可以改寫為：

式中：δ（t）為相位角。

由于高斯－馬爾卡夫（Gauss－Markov）過程能用來模擬許多物理現(xiàn)象，其中也包括模擬地震波的振幅X1（t），其離散形式為：

假定狀態(tài)空間的噪聲部分也可用高斯－馬爾卡夫過程代替，其表達式為：

式中：aX2，bX2均為高斯－馬爾卡夫過程的參數(shù)；w2與w1相互獨立。

這些方程中涉及的參數(shù)可通過聲發(fā)射信號的噪聲部分確定。

利用模型式（6），（7），聲發(fā)射的狀態(tài)空間系統(tǒng)可以構(gòu)造為：

此外，其量測方程構(gòu)造為：

式中：ωk為第k時間步長的主頻率；vk為量測噪聲。

式（8），（9）可分別用簡化的矩陣形式表示為：

由于沒有固定的主導(dǎo)角頻率，通過運用頻率粒子集［ωik］Mi＝1，會比運用具體指定的頻率有更高的估計精度。在研究中，每一個粒子的動態(tài)函數(shù)和量測函數(shù)分別為：

從式（12），（13）可知，每一個粒子的狀態(tài)和量測函數(shù)均為線性。從另一個角度來說，本研究提出了用RBPF的方法來估計聲發(fā)射信號，即由標準卡爾曼濾波處理狀態(tài)更新，同時利用粒子濾波處理非線性的量測函數(shù)。

3.2 基于聲發(fā)射信號去噪的RBPF算法

完整RBPF算法的步驟為：

1）卡爾曼濾波初始化：X1（0）＝0，X2（0）＝Y(jié)0；誤差協(xié)方差在實驗過程中，發(fā)現(xiàn)：R選擇相對較大值時，能得到一個較滿意的結(jié)果；遞推過程收斂時，R的選取對濾波結(jié)果影響甚小。

2）狀態(tài)估計預(yù)測和誤差協(xié)方差預(yù)測：

3）粒子更新：

①采樣ωik～p（ωk），zik＝sin（ωikkΔ），i＝1，2，…，M；

②計算外推量測估計和方差：

式中：Rk為量測函數(shù)中觀測噪聲的方差，即vk～（0，Rk）。

③計算預(yù)測密度：

④計算每個粒子的重要性權(quán)重：

⑥更新量測矩陣：G～k＝［1 z～k］。

8）設(shè)定k＝k＋1，回到步驟2）。

針對算法中出現(xiàn)的參數(shù)，需要說明：

Φk－1能夠根據(jù)狀態(tài)函數(shù)的先驗知識設(shè)定初值，Φk－1可分兩種情況考慮：①時變性，即隨時間需要更新或傳遞；②時不變性。在巖石沖擊聲發(fā)射實驗中，由于實驗環(huán)境相對較簡單，且RBPF用于聲發(fā)射信號去噪目前處于探索性階段，因此僅考慮時不變性，設(shè)定，其中：a1，a2均為高斯－馬爾卡夫過程中的參數(shù)，詳細的解釋可參考式（6），（7）。a1，a2合適的取值可以從信號的噪聲部分（沒有目標信號參與進行采集）和附有噪聲的信號序列中獲取。

Qk－1為狀態(tài)過程噪聲的協(xié)方差，其表達式為：

假定噪聲的方差為時不變隨機過程，那么，Qk－1為固定值Q。Q矩陣可通過信號的噪聲部分初始化，確定Q的步驟為：①采集沒有目標信號參與的噪聲信號，即采集周圍噪聲；②計算采集噪聲信號的自相關(guān)性；③利用模型σ2e－β｜τ｜對自相關(guān)數(shù)據(jù)進行擬合，其中：σ，β均通過回歸分析確定。

利用該方法可獲得σX2和βX2。在探討過程中，觀測噪聲的方差假定為相對較小的常數(shù)。

3.3 RBPF仿真實驗

為檢驗RBPF的濾波效果，借助于Matlab對其進行仿真模擬。首先需要構(gòu)造一個不受噪聲污染的原始信號作為標準信號，便于與其加噪處理后進行濾波所得結(jié)果比較。為驗證RBPF處理噪聲信號的效果，需對原始信號（如圖3（a）所示）加入高斯密度分布的噪聲，其附加噪聲的方差為0.2，附加噪聲的信號如圖3（b）所示。對加噪信號運用RBPF進行濾波處理，處理后的信號如圖3（c）所示。

圖3 仿真結(jié)果Fig.3 Simulation results

從圖3中可以看出，RBPF在信號去噪處理中能獲得較好的濾波效果。

3.4 RBPF在聲發(fā)射信號去噪中應(yīng)用

試驗數(shù)據(jù)來源：利用SHPB系統(tǒng)對兩組石灰?guī)r進行沖擊荷載作用下聲發(fā)射試驗，具體試驗裝置如圖4所示。采集兩組聲發(fā)射數(shù)據(jù)X11和X21，在采集過程中有噪聲干擾，利用本研究提出的RBPF處理采集的聲發(fā)射數(shù)據(jù)，對應(yīng)的去噪信號分別為X21和X22，經(jīng)RBPF后的結(jié)果如圖5所示。

圖4 SHPB測試系統(tǒng)Fig.4 SHPB testing systems

圖5 X11和X21為原始信號，X12和X22為其對應(yīng)的去噪信號Fig.5 X11and X21as original signals，X12and X22as corresponding denoising signals

在處理兩組聲發(fā)射數(shù)據(jù)時，粒子總數(shù)取200，噪聲按高斯－馬爾卡夫過程模擬，計算時間常數(shù)Tc＝4.5×10－4ms，方差σ＝6.7×10－3v2。聲發(fā)射源的采樣頻率40MHz，并遵循假設(shè)：開始波形的振幅為零，X1（1）＝0；第一個測量值y1設(shè)定為噪聲，即X2（1）＝y(tǒng)1。

為進一步分析聲發(fā)射數(shù)據(jù)并檢驗RBPF的效果，利用FFT求出圖5中各信號的頻譜圖，如圖6所示。從圖6中可以看出，X11最大振幅所對應(yīng)的頻率分別為0.15，1和3MHz，而X21最大振幅所對應(yīng)的頻率分別為0.15和1MHz。X21在0.15和1MHz對應(yīng)的幅值分別高于X11所對應(yīng)的點，而在3MHz時，兩者相差不大，甚至低于X11的，在3MHz出現(xiàn)的現(xiàn)象由噪聲所致。通過RBPF對原始信號進行去噪，發(fā)現(xiàn)兩者在3MHz的幅值都趨于零?？梢?，RBPF在處理聲發(fā)射信號中具有很大的潛力。

圖6 圖5中對應(yīng)信號的頻譜Fig.6 Spectrogram of the corresponding signal in Fig.5

4 結(jié)論

本研究提出了一種時空濾波方法RBPF，并對聲發(fā)射信號進行去噪處理。通過合理選擇模型參數(shù)，由實驗結(jié)果驗證了RBPF能夠加強聲發(fā)射信號的信噪比，并保留其主要頻率成分。實驗中還發(fā)現(xiàn)，沖擊荷載作用下石灰?guī)r最大振幅所對應(yīng)的頻率為0.15和1MHz，這對爆破與機械震動作用下的聲發(fā)射監(jiān)控有著重要的參考價值。

值得注意的是：論文研究模型最初來源于其他學(xué)者用來處理地震數(shù)據(jù)。由于考慮到聲發(fā)射與地震波具有一些共通性，通過修正其模型參數(shù)，延伸到RBPF處理聲發(fā)射信號中，并取得了不錯的效果。然而，對聲發(fā)射信號最優(yōu)模型的建立和過程，噪聲特征仍需作進一步的研究。

（

）：

［1］ Niri E D，F(xiàn)arhidzadeh A，Salamone S.Nonlinear Kalman filtering for acoustic emission source localization in anisotropic panels［J］.Ultrasonics，2014，54：486－501.

［2］ Baziw E，Jones I W.Application of Kalman filtering techniques for microseismic event detection［J］.Pure and Applied Geophysic，2002，159：449－471.

［3］ Zhang W M，Du G，Zhong S，et al.Study of nonlinear filter methods：Particle filter［J］.Journal of Systems Engineering and Electronics，2006，17（1）：1－5.

［4］ Petar M D，Jayesh H K，Zhang J Q，et al.Particle filtering［J］.IEEE Signal Processing Magazine，2003，20（5）：19－38.

［5］ Zhou C J，Zhang Y F.Particle filter based noise removal method for acoustic emission signals［J］.Mechanical Systems and Signal Processing，2012，28：63－77.

［6］ Menéndez R M，F(xiàn)reitas N，Poole D.Dynamic modelling and control of industrial processes with particle filtering algorithms［J］.Computer Aided Chemical Engineering，2004，18：721－726.

［7］ Freitas N.Rao－blackwellised particle filtering for fault diagnosis［J］.IEEE Aerospace Conference Proceedings，2002，4：1767－1772.

［8］ Mustière F，Bolic＇M，Bouchard M.Rao－Blackwellised particle filters：Examples of application［A］.IEEE CCECE／CCGEI［C］.Ottawa，Canada：［s.n.］，2006：1196－1200.

［9］ 潘宏俠，門吉芳.粒子濾波在軸承故障振動信號降噪中的應(yīng)用［J］.振動、測試與診斷，2011，31（3）：354－356.（Pan Hong－xia，Men Ji－fang.Bearing fault vibration signal noise reduction based on particle filtering［J］.Journal of Vibration，Measurement ＆ Diagnosis，2011，31（3）：354－356.（in Chinese））