基于粒子濾波的波疊加算法

孫玲莉 楊超 郭輝 胡定玉

摘 要： 與傳統(tǒng)波疊加法通過逆問題求解獲得源強(qiáng)不同，提出一種基于粒子濾波的波疊加算法來實(shí)現(xiàn)正向求解源強(qiáng)。該方法基于波疊加理論建立狀態(tài)空間模型，根據(jù)初始化粒子的先驗(yàn)信息預(yù)測并更新狀態(tài)向量。通過權(quán)重計(jì)算和重采樣，估計(jì)等效源的強(qiáng)度及其位置，進(jìn)而重建三維輻射聲場。仿真分析與實(shí)驗(yàn)驗(yàn)證表明，該算法可避免逆問題求解中的不適定性難題，實(shí)現(xiàn)聲場的高精度重建。

關(guān)鍵詞： 粒子濾波; 波疊加; 聲源定位; 聲場重建

文章編號： 2095-2163（2021）03-0080-05 中圖分類號：TP391.9; TB52 文獻(xiàn)標(biāo)志碼：A

【Abstract】Different from the traditional wave superposition method which is utilized to obtain the source intensity by solving the inverse problem， an innovative wave superposition algorithm based on particle filtering is proposed to realize the forward solving of source intensity. This method presents a state space model based on the wave superposition theory， then updates the state vector according to the prior information prediction of initializing particle. Through weight calculation and resampling， the intensity and the position of the equivalent source can be estimated. Thus， the three-dimensional sound field can be reconstructed. Simulation analysis and experimental verification show that the algorithm can avoid the ill-posed problem of solving the inverse problem and achieve high-precision reconstruction of the sound field.

【Key words】 particle filter; wave superposition; sound source localization; sound field reconstruction

0 引 言

波疊加法由Koopmann等人[1]在1989年提出，是一種求解聲輻射問題的積分方程法，其原理為：任何物體輻射的聲場均可以等效為該輻射體內(nèi)部若干個不同源強(qiáng)的等效源所產(chǎn)生聲場的疊加。利用波疊加法近場聲全息進(jìn)行聲場重建，既能克服基于空間Fourier變換法的近場聲全息[2]只能計(jì)算規(guī)則形狀聲源的缺點(diǎn)，又能避免基于邊界元法的近場聲全息[3]所帶來的復(fù)雜插值運(yùn)算和奇異積分處理困擾，國內(nèi)外學(xué)者對該方法進(jìn)行了深入研究[4-8]。

利用波疊加法進(jìn)行聲場重建時，其關(guān)鍵是獲取等效源強(qiáng)度和配置等效源位置。通常采用Tikhonov正則化方法[9-10]求解聲學(xué)逆問題中的病態(tài)方程，獲取準(zhǔn)確的等效源強(qiáng)度，再合理配置等效源位置，進(jìn)行聲場重建。Antoni[11]提出了貝葉斯正則化準(zhǔn)則的基本理論，Pereira等人[12]進(jìn)一步通過理論和實(shí)驗(yàn)對比了等效源法和貝葉斯法。Chardon等人[13]針對傳統(tǒng)近場聲全息逆問題正則化方法的復(fù)雜性，提出了一種新的基于壓縮采樣的正則化方法。Yang等人[14]提出了一種基于卡爾曼濾波的混合局部近場聲全息技術(shù)，旨在解決等效源配置不確定性問題。

以上所提方法均是從求解聲學(xué)逆問題的角度出發(fā)獲得源強(qiáng)，為避免求解逆問題的不適定性，Bai等人[15-16]提出基于等效源模型的自適應(yīng)濾波方法，對麥克風(fēng)陣列信號進(jìn)行處理，該方法利用遞歸方式估計(jì)未知源振幅，屬于聲學(xué)正問題。本文提出一種基于粒子濾波的波疊加方法（Particle Filter based WSM，PF-WSM），該方法以波疊加法為理論基礎(chǔ)，建立狀態(tài)空間模型，利用粒子濾波算法從正向估計(jì)等效源強(qiáng)度和位置，根據(jù)估計(jì)得到參數(shù)計(jì)算重建面聲壓值。本文安排如下：首先對PF-WSM算法進(jìn)行理論描述，再仿真分析聲源位置已知和未知兩種情況的聲場重建結(jié)果，最后通過實(shí)驗(yàn)驗(yàn)證，證明了該方法的可行性。

1 理論背景

1.1 PF-WSM算法理論

傳統(tǒng)波疊加法通過Tikhonov正則化方法求解聲學(xué)逆問題中的病態(tài)方程，獲得準(zhǔn)確的等效源強(qiáng)度，故被稱為Tikhonov正則化波疊加（TR-WSM）算法。與TR-WSM算法反解源強(qiáng)相比，PF-WSM算法正向估計(jì)源強(qiáng)，避免了求解逆問題的不適定性。

粒子濾波（Particle Filter，PF）是一種基于蒙特卡洛方法進(jìn)行遞推貝葉斯估計(jì)的濾波算法[17-18]。通過采用一組具有相應(yīng)權(quán)值的隨機(jī)粒子估計(jì)狀態(tài)的前驗(yàn)分布，并通過觀測值逼近狀態(tài)的后驗(yàn)分布，再求出這組粒子的加權(quán)求和作為狀態(tài)的估計(jì)值。

在PF-WSM算法中，用等效源強(qiáng)度和位置描述系統(tǒng)狀態(tài)，源強(qiáng)與傳遞矩陣計(jì)算的聲壓描述觀測狀態(tài)，由此，狀態(tài)空間模型表示為：

其中，qn表示源強(qiáng)的系統(tǒng)狀態(tài)矢量;sn表示源位置的系統(tǒng)狀態(tài)矢量且不隨遞歸變化;yn表示觀測的聲壓矢量;n表示估計(jì)狀態(tài)遞歸次數(shù)。源強(qiáng)的系統(tǒng)噪聲和觀測噪聲服從均值為零、方差分別為σ21和σ23的正態(tài)分布，即v1～N0，σ21和v3～N0，σ23;源位置的系統(tǒng)噪聲服從均勻分布，即v2～U-σ2，σ2。G表示等效源到全息面的傳遞矩陣。

第i個粒子的權(quán)值更新方程表示為：

當(dāng)經(jīng)過一定次數(shù)的估計(jì)狀態(tài)遞歸后，會存在粒子退化問題，可通過重采樣[19]方法來解決。該方法是對后驗(yàn)概率密度的離散近似表示再進(jìn)行一次采樣，繁殖權(quán)值較高的樣本，淘汰權(quán)值較低的樣本，重新生成一個新樣本集合，以克服樣本權(quán)值退化問題。

1.2 PF-WSM算法流程

粒子濾波的核心是利用優(yōu)選隨機(jī)粒子來實(shí)現(xiàn)狀態(tài)分布的估計(jì)，為更好地進(jìn)行源強(qiáng)分布估計(jì)，給定粒子初始狀態(tài)分布為高斯分布，這是因?yàn)楦鶕?jù)中心極限定理，隨著粒子數(shù)量的增加，大多數(shù)粒子均值變得相似并越來越接近總體均值。均值為μ和方差為σ2的高斯分布的概率密度函數(shù)為：

對于復(fù)數(shù)域的源強(qiáng)概率密度函數(shù)，式（7）可以看作是一個聯(lián)合概率密度函數(shù)的實(shí)部和虛部，即：

利用PF-WSM算法在2種情況下進(jìn)行聲場重建對比，根據(jù)公式（1）～公式（3）建立狀態(tài)空間模型。對這2種情況可簡述為：

（1）聲源位置已知，其觀測方程中聲壓值的變化僅與粒子的源強(qiáng)有關(guān)。

（2）聲源位置未知，聲壓值的變化依賴于粒子的源強(qiáng)和源位置。

文中給出了PF-WSM算法的流程如圖1所示。具體過程如下：

（1）粒子初始化：第一次遞歸時，從高斯分布Nμ，σ2中隨機(jī)產(chǎn)生K個粒子qini=1，2，…，K，并平均分配權(quán)重。

（2）系統(tǒng)狀態(tài)矢量預(yù)測：根據(jù)式（1）～式（2）所示的狀態(tài)轉(zhuǎn)移方程，計(jì)算粒子qin在n+1次遞歸的預(yù)測粒子qin+1，服從一階馬爾可夫過程。

（3）粒子權(quán)重更新：利用式（3）計(jì)算n+1次粒子預(yù)測值對應(yīng)的觀測聲壓，同時計(jì)算觀測聲壓與n+1次的實(shí)際測量聲壓之間的差值，根據(jù)差值更新每個粒子的權(quán)重（式（4）～式（5）），再將權(quán)重進(jìn)行歸一化。

（4）源強(qiáng)分布估計(jì)：利用當(dāng)前遞歸次數(shù)中每個粒子的權(quán)重win+1以及取值qin+1，估計(jì)系統(tǒng)狀態(tài)矢量（源強(qiáng)分布），計(jì)算公式為：

（5）重采樣：為了快速有效地選擇接近真實(shí)源強(qiáng)分布的粒子，在進(jìn)行下一次系統(tǒng)狀態(tài)估計(jì)之前，需要按照粒子的權(quán)重對粒子進(jìn)行重采樣。重復(fù)上述步驟，即可基于粒子濾波算法構(gòu)建系統(tǒng)狀態(tài)估計(jì)器，實(shí)現(xiàn)源強(qiáng)分布估計(jì)。

（6）聲場重建：根據(jù)上述構(gòu)建的系統(tǒng)狀態(tài)估計(jì)器，得到合理的估計(jì)源強(qiáng)分布，利用源強(qiáng)與重建面之間的傳遞矩陣求得重建聲壓：

其中，H表示等效源面到重建面的傳遞矩陣，Q表示PF-WSM算法估計(jì)源強(qiáng)分布。

2 數(shù)值仿真

2.1 仿真設(shè)置

為驗(yàn)證PF-WSM算法在聲場重建中的可行性，將其與Tikhonov正則化波疊加 （TR-WSM）算法進(jìn)行比較。單極子源的半徑為0.01 m，表面脈動速度為0.2 m/s，由25個單極子聲源組成的輻射聲場布置在xoy面，以聲源面中心為坐標(biāo)原點(diǎn)o，其法面方向?yàn)閦軸，建立空間直角坐標(biāo)系。全息面平行于xoy面，距坐標(biāo)原點(diǎn)0.2 m，在0.8×0.8 m2的平面上均勻布置25個傳聲器，間距為0.2 m。重建面尺寸為1×1 m2，網(wǎng)格間距為0.1 m，距坐標(biāo)原點(diǎn)0.05 m。對全息面的測量數(shù)據(jù)加上隨機(jī)高斯白噪聲，假設(shè)聲場中的信噪比為20 dB。設(shè)置遞歸次數(shù)N為90次，粒子數(shù)K為100個。

仿真過程中，相對誤差的計(jì)算公式如下：

2.2 仿真分析

聲源位置已知時，假設(shè)25個單極子源均勻分布在聲源面上，如圖2（a）所示，其源位置與全息面上傳聲器位置在xoy平面的俯視圖上一一對應(yīng)。根據(jù)1.2節(jié)的算法流程計(jì)算出重建聲壓，圖2（b）比較了100個粒子的PF-WSM算法和TR-WSM算法聲壓重建誤差。源位置均勻分布在網(wǎng)格點(diǎn)上，PF-WSM算法基于源強(qiáng)為高斯分布的先驗(yàn)信息僅估計(jì)源強(qiáng)，與TR-WSM算法相比，PF-WSM算法在聲壓重建中誤差較大。

聲源準(zhǔn)確位置已知的情況下，TR-WSM算法的性能要優(yōu)于PF-WSM算法。但在實(shí)際場景中，通常無法預(yù)先知道聲源的準(zhǔn)確位置，對源位置假設(shè)不正確，導(dǎo)致聲場重建值偏離實(shí)際值，因此假設(shè)聲源位置已知是一種理想化的情況。由此，假設(shè)源位置是未知的，并且存在位置擾動，更符合真實(shí)情況。

聲源位置未知時，假設(shè)25個單極子源分布在聲源面上存在擾動，源位置相對于傳聲器分布位置如圖3（a）所示。圖3（b）是PF-WSM算法和TR-WSM算法的重建誤差對比。在源位置未知且存在擾動的情況下，PF-WSM算法同時估計(jì)源強(qiáng)和源位置，基于源強(qiáng)服從高斯分布和源位置服從均勻分布的先驗(yàn)信息，在不斷遞歸之后，源位置越接近聲源真實(shí)位置，源強(qiáng)越接近真實(shí)源強(qiáng)，PF-WSM算法的聲場重建結(jié)果越接近真實(shí)輻射聲場。而TR-WSM算法在源位置未知的情況下，由于等效源配置不準(zhǔn)確，重建誤差增大。圖3（b）中PF-WSM算法的重建誤差小于TR-WSM算法，說明PF-WSM算法在源強(qiáng)和源位置均未知的情況下，估計(jì)效果更好，而TR-WSM算法在未知源位置時的重建效果較差，其應(yīng)用有一定局限性。圖4給出了2種算法在100 Hz時的聲場重建效果對比，圖4（a）為理論聲壓，圖4（b）為PF-WSM算法的重建聲壓，圖4（c）為TR-WSM算法的重建聲壓。由圖4可以看出，PF-WSM算法的聲場重建效果更好，其熱點(diǎn)最大值與理論聲壓相同，說明粒子濾波可以較準(zhǔn)確地估計(jì)出源強(qiáng)和源位置，而TR-WSM算法的熱點(diǎn)最大值小于理論值，這是因?yàn)椴恢缆曉次恢脮r，等效源配置不準(zhǔn)確，影響聲場重建效果。

3 實(shí)驗(yàn)驗(yàn)證

為進(jìn)一步檢驗(yàn)算法的可行性，在半消聲室進(jìn)行實(shí)驗(yàn)驗(yàn)證。實(shí)驗(yàn)設(shè)置如圖5（a）所示。半消聲室尺寸為9.8 m×8.6 m×3.5 m，背景噪聲為18 dB（A），截止頻率為125 Hz。采用音箱模擬噪聲源，頻率為300 Hz，音箱紙盆中心位置為（0， 0， 0）。以音箱紙盆中心為坐標(biāo)原點(diǎn)，音箱平面法向?yàn)閦軸建立空間直角坐標(biāo)系。傳聲器陣列平面平行于聲源面，陣列中心距坐標(biāo)原點(diǎn)0.6 m，離地面高度為1.2 m。傳聲器陣列如圖5（b）所示，采用36個傳聲器組成3個同心圓的面陣列，3個同心圓半徑分別為：0.07 m、0.18 m、0.25 m。重建面大小與傳聲器陣列面相同，重建面距坐標(biāo)原點(diǎn)0.05 m。PF-WSM算法中，粒子的先驗(yàn)信息是基于高斯分布，36個虛源均勻分布在0.25 m×0.25 m的正方形平面上，間隔為0.05 m，距坐標(biāo)原點(diǎn)-0.01 m。將PF-WSM算法和TR-WSM算法重建的聲壓與傳聲器測得的聲壓進(jìn)行對比。圖6為實(shí)驗(yàn)驗(yàn)證結(jié)果。圖6（a）為傳聲器陣列測得的聲壓分布，圖6（b）和圖6（c）分別為PF-WSM算法和TR-WSM算法重建的聲壓分布。由式（11）計(jì)算PF-WSM算法和TR-WSM算法重建誤差，分別為5.89%和13.43%，可以得出結(jié)果：PF-WSM算法的重建效果優(yōu)于TR-WSM算法的重建效果。

4 結(jié)束語

本文提出一種基于粒子濾波的波疊加算法進(jìn)行聲場重建。該方法在源位置未知且隨機(jī)的情況下，能夠有效利用參數(shù)的先驗(yàn)知識，通過粒子濾波估計(jì)器，達(dá)到求解源強(qiáng)的效果。特別是在聲源位置未知且存在擾動的情況下，源位置是非高斯分布的，傳統(tǒng)波疊加法由于等效源配置不準(zhǔn)確而聲場重建效果不好，PF-WSM算法則利用這一特點(diǎn)，可以更好地重建聲場。PF-WSM算法不僅可以估計(jì)源強(qiáng)，還可以估計(jì)源位置，甚至在源位置未知的情況下，估計(jì)性能更好，即對源位置的估計(jì)越準(zhǔn)確，重建的聲壓就越準(zhǔn)確。另一方面，粒子數(shù)的增加直接影響估計(jì)狀態(tài)向量的精度，粒子數(shù)越多，估計(jì)狀態(tài)越接近真實(shí)狀態(tài)。本文在權(quán)衡運(yùn)算時間和估計(jì)精度之后，選擇100個粒子進(jìn)行估計(jì)，若增加粒子數(shù)（如1 000個粒子），估計(jì)值將越接近真實(shí)值，但運(yùn)算時間相應(yīng)增加，后續(xù)可采用改進(jìn)粒子濾波算法進(jìn)行優(yōu)化，降低采樣復(fù)雜度，提升運(yùn)算速度。

參考文獻(xiàn)

[1] ??KOOPMANN G H， SONG L， FAHNLINE J B. A method for computing acoustic fields based on the principle of wave superposition[J]. The Journal of the Acoustical Society of America， 1989， 86（6）： 2433-2438.

[2] HALD J. Patch near-field holography using a new statistically optimal method[J]. Proceedings of Inter-noise， 2003（975）： 2203-2210.

[3] UCHIDA H，SAIJYOU K. Data extrapolation method for boundary element method-based near-field acoustical holography[J].The Journal of the Acoustical Society of America， 2004， 115（2）： 785-796.

[4] 楊殿閣， 李兵， 王子騰， 等. 運(yùn)動聲源識別的動態(tài)波疊加方法研究[J]. 物理學(xué)報(bào)， 2012， 61（5）：054306.

[5] FERNANDEZ-GRANDE E， XENAKI A， GERATOFT P. A sparse equivalent source method for near-field acoustic holography[J]. The Journal of the Acoustical Society of America， 2017， 141（1）： 532-542.

[6] 胡定玉， 李再幃， 方宇. 非自由聲場中目標(biāo)聲場還原與重建的等效源方法[J]. 聲學(xué)學(xué)報(bào)， 2017，42（4）： 465-475.

[7] XU Zhongming， WANG Qinghua， HE Yangsong， et al. A monotonic two-step iterative shrinkage/thresholding algorithm for sound source identification based on equivalent source method[J]. Applied Acoustics， 2018， 129： 386-396.

[8] BI Chuanxing， LIU Yuan， ZHANG Yongbin， et al. Extension of sound field separation technique based on the equivalent source method in a sparsity framework[J]. Journal of Sound and Vibration， 2018， 442： 125-137.

[9] TIKHONOV A N. Solution of incorrectly formulated problems and the regularization method [J]. Soviet Math Dokl， 1963， 4： 1035-1038.

[10]LECLERE Q. Acoustic imaging using under-determined inverse approaches： Frequency limitations and optimal regularization[J]. Journal of Sound and Vibration， 2009， 321（3-5）： 605-619.

[11]ANTONIO J. Bayesian focusing：A unified approach to inverse acoustic radiation[C]//Proceedings of ISMA 2010. ?Leuven， Belgium： Katholieke Universiteit Leuven， 2010： 35-45.

[12]PEREIRA A， LECLRE Q， ANTONI J. A theoretical and experimental comparison of the equivalent source method and a Bayesian approach to noise source identification[C]//BeBec-2012，2012： 4.

[13]CHARDON G， DAUDET L， PEILLOT A， et al. Nearfield Acoustic Holography using sparsity and compressive sampling principles[J]. Journal of the Acoustical Society of America， 2012， 132（3）： 1521-1534.

[14]YANG C， WANG Y S， Guo H. Hybrid patch near-field acoustic holography based on Kalman filter[J]. Applied Acoustics， 2019， 148： 23-33.

[15]BAI M R， CHEN C C. Kalman filter-based microphone array signal processing using the equivalent source model[J]. Journal of Sound and Vibration， 2012， 331（22）： 4940-4955.

[16]BAI M R， AGARWAL A， CHEN C C， et al. Bayesian approach of nearfield acoustic reconstruction with particle filters[J].The Journal of the Acoustical Society of America， 2013， 133（6）： 4032-4043.

[17]王法勝， 魯明羽， 趙清杰， 等. 粒子濾波算法[J]. 計(jì)算機(jī)學(xué)報(bào)， 2014， 37（8）： 1679-1694.

[18]昝孟恩， 周航， 韓丹， 等. 粒子濾波目標(biāo)跟蹤算法綜述[J]. 計(jì)算機(jī)工程與應(yīng)用， 2019， 55（5）： 8-17.

[19]孔紅山， 李小鵬， 郁濱. SIR粒子濾波的改進(jìn)算法[J]. 計(jì)算機(jī)工程與設(shè)計(jì)， 2020，41（7）： 1899-1904.