孫玲莉 楊超 郭輝 胡定玉

摘 要： 與傳統(tǒng)波疊加法通過逆問題求解獲得源強不同，提出一種基于粒子濾波的波疊加算法來實現(xiàn)正向求解源強。該方法基于波疊加理論建立狀態(tài)空間模型，根據(jù)初始化粒子的先驗信息預(yù)測并更新狀態(tài)向量。通過權(quán)重計算和重采樣，估計等效源的強度及其位置，進(jìn)而重建三維輻射聲場。仿真分析與實驗驗證表明，該算法可避免逆問題求解中的不適定性難題，實現(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)部若干個不同源強的等效源所產(chǎn)生聲場的疊加。利用波疊加法近場聲全息進(jìn)行聲場重建，既能克服基于空間Fourier變換法的近場聲全息[2]只能計算規(guī)則形狀聲源的缺點，又能避免基于邊界元法的近場聲全息[3]所帶來的復(fù)雜插值運算和奇異積分處理困擾，國內(nèi)外學(xué)者對該方法進(jìn)行了深入研究[4-8]。

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

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

1 理論背景

1.1 PF-WSM算法理論

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

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

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

其中，qn表示源強的系統(tǒng)狀態(tài)矢量;sn表示源位置的系統(tǒng)狀態(tài)矢量且不隨遞歸變化;yn表示觀測的聲壓矢量;n表示估計狀態(tài)遞歸次數(shù)。源強的系統(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ù)的估計狀態(tài)遞歸后，會存在粒子退化問題，可通過重采樣[19]方法來解決。該方法是對后驗概率密度的離散近似表示再進(jìn)行一次采樣，繁殖權(quán)值較高的樣本，淘汰權(quán)值較低的樣本，重新生成一個新樣本集合，以克服樣本權(quán)值退化問題。

1.2 PF-WSM算法流程

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

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

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

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

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

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

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

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

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

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

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

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

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

2 數(shù)值仿真

2.1 仿真設(shè)置

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

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

2.2 仿真分析

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

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

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

3 實驗驗證

為進(jìn)一步檢驗算法的可行性，在半消聲室進(jì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)原點，音箱平面法向為z軸建立空間直角坐標(biāo)系。傳聲器陣列平面平行于聲源面，陣列中心距坐標(biāo)原點0.6 m，離地面高度為1.2 m。傳聲器陣列如圖5（b）所示，采用36個傳聲器組成3個同心圓的面陣列，3個同心圓半徑分別為：0.07 m、0.18 m、0.25 m。重建面大小與傳聲器陣列面相同，重建面距坐標(biāo)原點0.05 m。PF-WSM算法中，粒子的先驗信息是基于高斯分布，36個虛源均勻分布在0.25 m×0.25 m的正方形平面上，間隔為0.05 m，距坐標(biāo)原點-0.01 m。將PF-WSM算法和TR-WSM算法重建的聲壓與傳聲器測得的聲壓進(jìn)行對比。圖6為實驗驗證結(jié)果。圖6（a）為傳聲器陣列測得的聲壓分布，圖6（b）和圖6（c）分別為PF-WSM算法和TR-WSM算法重建的聲壓分布。由式（11）計算PF-WSM算法和TR-WSM算法重建誤差，分別為5.89%和13.43%，可以得出結(jié)果：PF-WSM算法的重建效果優(yōu)于TR-WSM算法的重建效果。

4 結(jié)束語

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

參考文獻(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] 楊殿閣， 李兵， 王子騰， 等. 運動聲源識別的動態(tài)波疊加方法研究[J]. 物理學(xué)報， 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é)報， 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]. 計算機學(xué)報， 2014， 37（8）： 1679-1694.

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

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