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

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

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

1 粒子濾波

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

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

粒子濾波算法流程如圖1所示，其詳細(xì)過程見文獻(xiàn)［3，4］。

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

2 Rao－Blackwellised粒子濾波

2.1 RBPF方法

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

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

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

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

2.2 RBPF實(shí)時(shí)處理特征

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

圖2 RBPF實(shí)時(shí)處理信號(hào)濾波的特征Fig.2 The real－time processing feature of RBPF based signal filtering

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

完整RBPF算法的步驟為：

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

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

3）粒子更新：

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

②計(jì)算外推量測估計(jì)和方差：

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

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

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

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

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

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

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

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

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

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

3.3 RBPF仿真實(shí)驗(yàn)

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

（

）：

［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］ 潘宏俠，門吉芳.粒子濾波在軸承故障振動(dòng)信號(hào)降噪中的應(yīng)用［J］.振動(dòng)、測試與診斷，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））