基于3D-MUSIC和波疊加的聲場(chǎng)重建方法*

孫玲莉 楊 超 郭 輝 胡定玉 顧汝彬

（1.上海工程技術(shù)大學(xué)機(jī)械與汽車工程學(xué)院 上海 201600）（2.上海工程技術(shù)大學(xué)城市軌道交通學(xué)院 上海 201600）（3.中國人民解放軍32128部隊(duì) 濟(jì)南 250000）

1 引言

近場(chǎng)聲全息（Near-filed Acoustic Holography，NAH）作為一種聲源定位和聲場(chǎng)可視化技術(shù)，可為振動(dòng)、噪聲源分析提供豐富的聲源和聲場(chǎng)信息，對(duì)噪聲源控制及聲輻射研究具有重要價(jià)值。近場(chǎng)聲全息技術(shù)包含多種實(shí)現(xiàn)方法［1～4］，其中，波疊加法對(duì)各類聲場(chǎng)的適應(yīng)性較好，且不用處理復(fù)雜積分問題，而得到廣泛運(yùn)用。如：楊殿閣等［5］提出動(dòng)態(tài)波疊加法計(jì)算運(yùn)動(dòng)聲源的聲場(chǎng)；Bai等［6］將粒子濾波算法與波疊加法相結(jié)合，提出從正向求解等效源源強(qiáng)；為減少空間采樣點(diǎn)數(shù)，Bi等［7］提出一種壓縮模態(tài)等效點(diǎn)源法，Hu等［8］在壓縮感知框架下對(duì)等效源法進(jìn)行了改進(jìn)。

利用波疊加法重建聲場(chǎng)，其計(jì)算精度取決于對(duì)等效源強(qiáng)度的求解，這屬于聲學(xué)逆問題，需采用正則化進(jìn)行求解，已提出多種較為成熟的解決方法［9～10］。除此之外，等效源點(diǎn)的布置對(duì)等效源強(qiáng)度的求解起重要作用，影響波疊加法重建聲場(chǎng)的精度。李加慶等［11］提出利用波束形成進(jìn)行聲源定位，指導(dǎo)等效源配置，再用波疊加法重建聲場(chǎng)。該方法是已知聲源與測(cè)量面的距離，在二維平面上進(jìn)行聲源定位，對(duì)于聲源空間位置完全未知的情況，重建精度會(huì)受到影響。楊超等［12］提出用統(tǒng)計(jì)最優(yōu)近場(chǎng)聲全息兩次定位聲源二維平面的位置，再通過三角函數(shù)定理求解聲源z軸距離，實(shí)現(xiàn)聲源的三維定位，再配置等效源。為進(jìn)一步提高重建精度，Yang等［13］提出多次轉(zhuǎn)動(dòng)全息面進(jìn)行聲源三維定位，獲得足夠多的聲源位置信息，用卡爾曼濾波提高聲源定位精度。該方法將多種算法聯(lián)合，提高了波疊加法重建聲場(chǎng)的精度，但其計(jì)算分析較復(fù)雜，在實(shí)際應(yīng)用中受到一定限制。

針對(duì)上述問題，提出一種基于3D-MUSIC算法和波疊加法的局部聲場(chǎng)重建方法。3D-MUSIC算法通過空間譜矩陣的求解及譜峰搜索，獲得聲源的俯仰角、方位角和距離等空間位置數(shù)據(jù)，可直接準(zhǔn)確定位聲源的三維位置。根據(jù)3D-MUSIC算法搜索的位置信息，確定等效源的配置策略，利用波疊加法重建聲場(chǎng)。

2 理論基礎(chǔ)

2.1 3D-MUSIC算法

MUSIC算法是一種基于子空間特征分解的高分辨率譜估計(jì)算法，根據(jù)信號(hào)子空間和噪聲子空間的正交特性，構(gòu)造出空間譜矩陣，通過空間譜矩陣的求解及譜峰搜索估計(jì)出信號(hào)的DOA參數(shù)。但經(jīng)典MUSIC算法［14］僅適用于二維遠(yuǎn)場(chǎng)聲源定位。改進(jìn)的3D-MUSIC算法［15～16］可計(jì)算近場(chǎng)噪聲源的俯仰角、方位角和距離，實(shí)現(xiàn)聲源的三維空間定位。

如圖1所示，假設(shè)坐標(biāo)原點(diǎn)o為參考點(diǎn)，任一聲源sd與參考點(diǎn)o的連線在xoy平面的投影oc與x軸的夾角θd為方位角，聲源sd與參考點(diǎn)o的連線與z軸的夾角ψd為俯仰角，rd為聲源sd到參考點(diǎn)o的距離，rmd為第m個(gè)傳聲器到聲源sd的距離。

圖1 傳聲器陣列接收信號(hào)示意

假設(shè)空間有D個(gè)近場(chǎng)窄帶信號(hào)源(s1(t)，…，sD(t))，輻射到由M個(gè)傳聲器組成的陣列上，第m個(gè)傳聲器接收到的信號(hào)表達(dá)式為

其中，d=1，2，…，D，m=1，2，…，M，nm(t)表示第m個(gè)傳聲器在t時(shí)刻接收到的高斯白噪聲。整個(gè)傳聲器陣列的接收信號(hào)表達(dá)式為

式中，S(t)為聲源，A是陣列方向矢量，N(t)是高斯白噪聲。其中：

由圖1的幾何關(guān)系可知：

其中，rmd代表第m個(gè)傳聲器與第d個(gè)聲源與的距離；ωc為聲源的中心角頻率；τmd為第m個(gè)傳聲器接收到第d個(gè)聲源信號(hào)的相對(duì)時(shí)延。

對(duì)傳聲器陣列接收到的數(shù)據(jù)進(jìn)行協(xié)方差運(yùn)算，根據(jù)信號(hào)和噪聲互相獨(dú)立的特性，協(xié)方差矩陣可以將信號(hào)和噪聲分離。定義測(cè)量聲壓P(t)的協(xié)方差矩陣R(t)為

式中，ARs(t)AH為信號(hào)部分，σ2I為噪聲部分，Rs(t)為原始信號(hào)的協(xié)方差矩陣。對(duì)R(t)進(jìn)行特征分解可得：

式中，Us為信號(hào)子空間，Un為噪聲子空間。由噪聲子空間與信號(hào)子空間相互正交可得：

構(gòu)造空間譜函數(shù)為

3D-MUSIC算法對(duì)空間譜函數(shù)進(jìn)行網(wǎng)格搜索，進(jìn)而輸出空間譜函數(shù)為極大值的ψ，θ，r值，即空間聲源信號(hào)的俯仰角、方位角和距離。根據(jù)該聲源定位信息指導(dǎo)等效源的配置，進(jìn)而利用波疊加法重建聲場(chǎng)。

2.2 波疊加算法

波疊加法的基本思想是將聲源輻射聲場(chǎng)用聲源內(nèi)部分布的一系列等效源產(chǎn)生的聲場(chǎng)疊加替代，由全息面測(cè)得聲壓反求出等效源的強(qiáng)度，從而實(shí)現(xiàn)整個(gè)聲場(chǎng)的重建與預(yù)測(cè)。其中聲源面、重建面和全息面的位置關(guān)系如圖2所示。

圖2 聲場(chǎng)空間位置示意

由式（2）可知，全息面上M個(gè)傳聲器測(cè)量的聲壓信息為P(t)，是時(shí)域信號(hào)，經(jīng)過傅里葉變換為頻域信號(hào)，再進(jìn)行下一步計(jì)算。假設(shè)等效源面上布置N個(gè)等效源，則全息面上的聲壓可表示為

式中，P?=[p(r1)p(r2)…p(rM)]T，Q為各等效源源強(qiáng)組成的列向量。G為等效源到全息面的傳遞矩陣，其中元素：

其中，g(rm，rn)為全息面上第m個(gè)傳聲器到第n個(gè)等效源的格林函數(shù)。

由式（9）可反解出源強(qiáng)：

式中，G+為G的廣義逆矩陣。為了保證解的唯一性，等效源的總數(shù)N應(yīng)小于等于全息面上的測(cè)量點(diǎn)數(shù)M。采用Tikhonov正則化方法求解等效源強(qiáng)度，正則化參數(shù)由Hald經(jīng)驗(yàn)公式［17］確定。

由此，可計(jì)算出重建面聲壓：

式中，H為等效源面到重建面的傳遞矩陣。

3 數(shù)值仿真

3.1 仿真設(shè)置

數(shù)值仿真參數(shù)設(shè)置如下：選取兩個(gè)振動(dòng)頻率為1400Hz的非相干脈動(dòng)球聲源，聲源半徑為0.005m，表面脈動(dòng)速度為0.02m/s。根據(jù)3D-MUSIC算法定位特點(diǎn)，聲源可設(shè)在距傳聲器陣列中心0.01m～1m、俯仰角1°～90°、方位角1°～360°的任意位置，設(shè)置兩個(gè)聲源在笛卡爾坐標(biāo)系下的位置分別為s1=（-0.4m，0.2m，0.6m）、s2=（0.2m，-0.1m，0.6m）。由傳聲器組成的全息面位于坐標(biāo)原點(diǎn)所在平面，傳聲器陣列如圖1中的全息面所示，采用36個(gè)傳聲器組成三個(gè)同心圓的面陣列，三個(gè)同心圓半徑分別為0.07m、0.18m、0.25m。重建面位置在0.59m處，面積為1×1m2，網(wǎng)格點(diǎn)間隔為0.1m。

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

式中，P(i)為聲壓理論值；Pr(i)為聲壓重建值。

3.2 聲源定位

由3D-MUSIC算法定位聲源位置，其俯仰角、方位角和距離隨信噪比變化的均方根誤差曲線如圖3所示。隨著信噪比的增加，均方根誤差越來越小，這是因?yàn)?D-MUSIC算法中噪聲子空間和信號(hào)子空間的分離更加徹底。根據(jù)圖3的定位誤差分析，選擇在信噪比為20dB的條件下進(jìn)行仿真。表1為在笛卡爾坐標(biāo)系下，聲源理論坐標(biāo)與算法定位結(jié)果的數(shù)值對(duì)比，兩個(gè)聲源的x軸定位誤差絕對(duì)值分別為0.009m和0.002m，y軸定位誤差絕對(duì)值均為0.001m，z軸定位誤差絕對(duì)值分別為0.017m和0.011m。由表1可以看出，在xoy平面上誤差在1cm以內(nèi)，z方向上誤差在1cm～2cm之間，滿足定位精度要求。

圖3 均方根誤差曲線

表1 理論位置與定位結(jié)果對(duì)比

3.3 聲場(chǎng)重建

由3D-MUSIC算法定位結(jié)果可知，聲源估計(jì)位置分別為（-0.391m，0.199m，0.583m）和（0.198m，-0.101m，0.611m）。根據(jù)該定位信息配置等效源，文獻(xiàn)［18～19］指出：等效源配置最好能與振動(dòng)體表面共形，且等效源球面半徑與脈動(dòng)球半徑之比在0.8以下，有較好的重建效果，另一方面，當(dāng)?shù)刃г磾?shù)目增加時(shí)，重建誤差會(huì)隨之減小。由此制定等效源配置策略：等效源面是一個(gè)以定位聲源位置為中心、0.01m為半徑的球面，在兩個(gè)聲源的等效源球面上分別布置16個(gè)等效源點(diǎn)。

圖4、5給出了兩種方法的重建效果。圖4是在z軸為0.61m的平面，等效源面大小為0.5×0.5m2，均勻布置6×6個(gè)等效源點(diǎn)的情況下，傳統(tǒng)波疊加法的重建值與理論值的對(duì)比結(jié)果，可以看出傳統(tǒng)波疊加法重建聲場(chǎng)的精度不高，這是因?yàn)榈刃г淳鶆蚍植荚谄矫嫔蠒r(shí)，等效源強(qiáng)度被均分，聲源位置處沒有重點(diǎn)布置等效源，重建值的峰值小于理論值。圖5是本文所提方法的重建值與理論值的對(duì)比結(jié)果，其重建值與理論值幾乎重合，重建精度較高。由此說明：利用3D-MUSIC算法對(duì)聲源位置進(jìn)行估計(jì)，在聲源估計(jì)位置處密集布置等效源，聲場(chǎng)的重建值更接近理論值。圖6是兩種方法的重建誤差值隨頻率的變化趨勢(shì)，傳統(tǒng)波疊加法在1300Hz～3000Hz范圍內(nèi)的相對(duì)誤差均大于15%，而本文所提方法的誤差在10%以下，重建誤差對(duì)比明顯。由此可以看出，本文所提方法可有效降低重建誤差，對(duì)聲場(chǎng)進(jìn)行準(zhǔn)確重建。

圖4 傳統(tǒng)波疊加法重建結(jié)果

圖5 本文所提方法重建結(jié)果

圖6 重建誤差對(duì)比

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

為進(jìn)一步檢驗(yàn)算法的可行性與準(zhǔn)確性，在半消聲室進(jìn)行實(shí)驗(yàn)驗(yàn)證。實(shí)驗(yàn)設(shè)置如圖7所示。半消聲室尺寸為9.8m×8.6m×3.5m，背景噪聲為18dB（A），截止頻率為125Hz。采用兩個(gè)音箱模擬噪聲源，音箱位置分別為s1=（-0.4m，0.2m，0.6m）、s2=（0.2m，-0.1m，0.6m）。以傳聲器陣列的中心位置為坐標(biāo)系原點(diǎn)，確定水平坐標(biāo)x軸、豎直坐標(biāo)y軸和測(cè)量坐標(biāo)z軸，傳聲器陣列離地面高度為1.2m。傳聲器陣列為三個(gè)同心圓組成的面陣列，其參數(shù)與仿真設(shè)置一致。兩個(gè)音箱發(fā)出的聲信號(hào)是頻率為1400Hz的非相干聲源，利用傳聲器陣列和數(shù)據(jù)采集系統(tǒng)采集數(shù)據(jù)，再利用本文所提方法進(jìn)行實(shí)驗(yàn)驗(yàn)證。

圖7 實(shí)驗(yàn)布置

笛卡爾坐標(biāo)系下，聲源理論坐標(biāo)與3D-MUSIC算法定位結(jié)果的數(shù)值對(duì)比如表2所示。由表2可得，實(shí)驗(yàn)中兩個(gè)音箱的算法定位誤差絕對(duì)值在x軸分別為0.019m和0.028m，在y軸分別為0.020m和0.012m，在z軸 分 別 為0.017m和0.004m。利 用3D-MUSIC算法進(jìn)行聲源定位的實(shí)驗(yàn)誤差在可接受范圍內(nèi)，能夠較為準(zhǔn)確地定位出聲源的空間三維坐標(biāo)。以聲源定位坐標(biāo)（-0.381m，0.220m，0.583m）和（0.172m，-0.112m，0.596m）為中心，以0.01m為半徑的球面上，分別配置16個(gè)等效源，利用波疊加法進(jìn)行聲場(chǎng)重建，重建面在z=0.59m處，重建面大小為1×1m2，網(wǎng)格點(diǎn)間距為0.1m。

表2 理論位置與定位結(jié)果對(duì)比

聲場(chǎng)重建結(jié)果如圖8、9所示，圖8為傳統(tǒng)波疊加法重建值與理論值對(duì)比，圖9為本文所提方法重建值與理論值對(duì)比。實(shí)驗(yàn)結(jié)果表明，傳統(tǒng)波疊加法的聲場(chǎng)重建值與理論值相差較大，這是因?yàn)槠矫嫔暇鶆蚺渲玫刃г磿r(shí)，沒有對(duì)聲源位置處有針對(duì)性的布置等效源，使得重建峰值與理論峰值相差較大；而本文所提方法可以較為準(zhǔn)確地重建聲場(chǎng)，說明在聲源估計(jì)位置處密集布置等效源，可提高聲場(chǎng)重建精度。圖10為兩種方法的重建誤差對(duì)比，可以看出，在1300Hz～3000Hz頻率范圍內(nèi)，本文所提方法的重建誤差明顯低于傳統(tǒng)波疊加法的重建誤差，能夠更好地重建聲場(chǎng)，保證重建精度，驗(yàn)證了本文所提方法的可行性與準(zhǔn)確性。

圖9 本文所提方法重建結(jié)果

圖10 重建誤差對(duì)比

5 結(jié)語

為解決等效源配置不確定性問題，提出了基于3D-MUSIC和波疊加的聲場(chǎng)重建方法。數(shù)值仿真和實(shí)驗(yàn)研究結(jié)果表明，該方法在1300Hz～3000Hz頻率范圍內(nèi)，可有效定位聲源三維位置，準(zhǔn)確重建聲場(chǎng)且重建誤差不高于10%。結(jié)論如下：

1）與傳統(tǒng)波疊加法的均勻分布方式相比，基于3D-MUSIC和波疊加的聲場(chǎng)重建方法的等效源布置是非均勻的，重點(diǎn)聲源布置更多的等效源點(diǎn)，能夠更好地模擬輻射聲場(chǎng)。

2）該方法精準(zhǔn)配置等效源，降低了不必要的等效源數(shù)目，降低了求解方程的難度，提升了求解效率。

3）與現(xiàn)有的復(fù)雜三維定位方法相比，該方法利用3D-MUSIC算法直接定位聲源的空間位置信息，在保證準(zhǔn)確重建聲場(chǎng)的基礎(chǔ)上，降低了計(jì)算復(fù)雜度，易于實(shí)現(xiàn)。

此外，本文實(shí)驗(yàn)所用為兩個(gè)音箱模擬輻射聲場(chǎng)，后續(xù)可用板聲源或其他型聲源進(jìn)行深入研究；本文采用36個(gè)測(cè)點(diǎn)進(jìn)行聲場(chǎng)數(shù)據(jù)采集，在保證定位結(jié)果及重建精度的情況下，可進(jìn)一步研究如何減少采樣點(diǎn)數(shù)，如引入壓縮感知理論等。