用剛性表面球形傳聲器陣列重構(gòu)入射波聲場(chǎng)與識(shí)別定位聲源*

李敏宗，盧奐采，2*，金江明

（1.浙江工業(yè)大學(xué)機(jī)械工程學(xué)院，特種裝備制造與先進(jìn)加工技術(shù)教育部/浙江省重點(diǎn)實(shí)驗(yàn)室，杭州310014；2.浙江工業(yè)大學(xué)，浙江省信號(hào)處理重點(diǎn)實(shí)驗(yàn)室，杭州310014）

用剛性表面球形傳聲器陣列重構(gòu)入射波聲場(chǎng)與識(shí)別定位聲源*

李敏宗1，盧奐采1，2*，金江明1

（1.浙江工業(yè)大學(xué)機(jī)械工程學(xué)院，特種裝備制造與先進(jìn)加工技術(shù)教育部/浙江省重點(diǎn)實(shí)驗(yàn)室，杭州310014；2.浙江工業(yè)大學(xué)，浙江省信號(hào)處理重點(diǎn)實(shí)驗(yàn)室，杭州310014）

剛性表面球形傳聲器陣列是測(cè)量三維聲場(chǎng)常用的前端，由于剛性球面對(duì)入射波聲場(chǎng)有散射影響，直接測(cè)量值不是原入射波聲場(chǎng)的聲壓值，因而不能使用在自由聲場(chǎng)條件下建立的聲場(chǎng)模型和近場(chǎng)聲全息方法來(lái)重構(gòu)入射波聲場(chǎng)。通過(guò)剛性球體表面聲波散射的數(shù)學(xué)模型，建立起入射波聲壓與發(fā)生散射后總聲壓之間的關(guān)系，進(jìn)而通過(guò)剛性表面球形傳聲器陣列測(cè)量聲壓，計(jì)算入射波聲場(chǎng)分布，并進(jìn)行聲源識(shí)別定位。通過(guò)仿真與實(shí)驗(yàn)，檢驗(yàn)了球面入射波遇到剛性表面球體發(fā)生散射后，球體表面處聲壓分布的變化。由球面近場(chǎng)聲全息方法，重構(gòu)球形陣列表面及外部空間的入射波聲場(chǎng)，檢驗(yàn)了波數(shù)、重構(gòu)距離參數(shù)變化對(duì)聲場(chǎng)重構(gòu)精度的影響。結(jié)果表明：采用剛性表面球形傳聲器陣列測(cè)得的聲場(chǎng)總聲壓，運(yùn)用根據(jù)聲波散射模型建立的球面近場(chǎng)聲全息方法，可以以一定的精度重構(gòu)出入射波聲場(chǎng)的三維空間分布。

球形傳聲器陣列；入射波聲場(chǎng)重構(gòu)；球面近場(chǎng)聲全息；聲波散射；剛性表面球體

球形傳聲器陣列，由于其球面所具有的經(jīng)緯角全指向性和三維對(duì)稱(chēng)性，能夠測(cè)量聲壓場(chǎng)的三維空間分布，因而對(duì)于封閉空間（如飛機(jī)、汽車(chē)和房間內(nèi)部）的聲場(chǎng)重構(gòu)、聲源識(shí)別定位等應(yīng)用，具有獨(dú)特優(yōu)勢(shì)［1-5］。球形傳聲器陣列是指在一個(gè)固定半徑的球面上按照一定的幾何規(guī)律和優(yōu)化準(zhǔn)則［6］，布置一組傳聲器而形成的聲場(chǎng)全息聲壓測(cè)量陣列。按照結(jié)構(gòu)可以分為空心球形陣列與剛性表面球形陣列，空心球形陣列指的是將各個(gè)傳聲器固定在指定的位置，傳聲器的表面形成一個(gè)虛擬球面的若干個(gè)節(jié)點(diǎn)，它的聲場(chǎng)理論模型假設(shè)聲場(chǎng)不受陣列影響，然而在應(yīng)用時(shí)由于聲場(chǎng)中引入了傳聲器支架與導(dǎo)線(xiàn)等因素，實(shí)際上（尤其在聲音頻率較高時(shí)）滿(mǎn)足不了此假設(shè)條件［3］。剛性表面球形陣列是指將傳聲器鑲嵌到具有剛性表面的球面殼體上，由于剛性表面對(duì)入射聲波產(chǎn)生的散射，球形陣列上傳聲器測(cè)量的聲壓值，不是原入射波聲場(chǎng)的聲壓值，因此不能直接將其應(yīng)用于自由聲場(chǎng)條件下建立的聲場(chǎng)模型和近場(chǎng)聲全息方法，來(lái)重構(gòu)入射波聲場(chǎng)，而需要考慮球形傳聲器陣列的剛性表面在聲場(chǎng)中發(fā)生的散射現(xiàn)象，建立考慮散射現(xiàn)象的聲場(chǎng)的數(shù)學(xué)模型。所幸的是，剛性表面球體對(duì)入射波聲場(chǎng)的散射在聲學(xué)領(lǐng)域是一個(gè)經(jīng)典問(wèn)題，其理論分析與數(shù)值計(jì)算結(jié)果在理論聲學(xué)相關(guān)專(zhuān)著中［7-9］均具有闡述，本文以入射波在剛性球形表面產(chǎn)生散射的數(shù)學(xué)模型為基礎(chǔ)，推導(dǎo)對(duì)應(yīng)的球面近場(chǎng)聲全息算法，最終由剛性表面球形陣列上傳聲器所測(cè)得的聲場(chǎng)總聲壓數(shù)據(jù)，重構(gòu)出陣列周?chē)臻g的入射波聲壓的三維空間分布，進(jìn)而進(jìn)行聲源識(shí)別定位。

目前，剛性表面球形傳聲器陣列多用于遠(yuǎn)場(chǎng)測(cè)試，結(jié)合波束形成方法［10-12］進(jìn)行三維空間的聲源識(shí)別定位。由于剛性表面球形傳聲器陣列在近場(chǎng)可結(jié)合近場(chǎng)聲全息方法進(jìn)行聲場(chǎng)重構(gòu)，在球狀或類(lèi)似于球狀的內(nèi)部封閉空間中的聲源識(shí)別定位有獨(dú)特優(yōu)勢(shì)，故在該方向的研究也有一定進(jìn)展。2008年，F(xiàn)inn Jacob?sen［2］在球面近場(chǎng)聲全息方法中，考慮了剛性表面球形傳聲器陣列的聲波散射影響，在仿真中給出了陣列中心與聲源連線(xiàn)上的重構(gòu)聲壓分布，將其與基準(zhǔn)聲壓進(jìn)行了對(duì)比，并在實(shí)驗(yàn)中比較了空間中單個(gè)點(diǎn)位置處的重構(gòu)聲壓與測(cè)量聲壓的幅值，但沒(méi)有進(jìn)行三維空間的聲場(chǎng)重構(gòu)誤差的分析與聲源識(shí)別定位的驗(yàn)證；2010年，Earl G.Williams［3］使用剛性表面球形傳聲器陣列，進(jìn)行了聲場(chǎng)重構(gòu)的研究，在不同半徑的重構(gòu)球面上，進(jìn)行了三維空間的聲場(chǎng)重構(gòu)誤差分析，并在點(diǎn)聲源的聲場(chǎng)環(huán)境下，進(jìn)行了聲場(chǎng)重構(gòu)與聲源識(shí)別定位的驗(yàn)證，其論文也通過(guò)考慮剛性球面的散射模型，得出了修正的傳播因子，但根據(jù)此修正傳播因子計(jì)算出的聲場(chǎng)為陣列表面及外部空間聲場(chǎng)的總聲壓分布，而不是僅僅由入射波形成的聲場(chǎng)的聲壓分布。本文的研究?jī)?nèi)容為：在球面近場(chǎng)聲全息方法中應(yīng)用剛性球面的聲波散射模型，計(jì)算出重構(gòu)的入射波聲壓在整個(gè)三維空間的分布，并在不同波數(shù)下，進(jìn)行不同半徑的重構(gòu)球面上聲場(chǎng)重構(gòu)誤差分析；通過(guò)重構(gòu)出的三維空間入射波聲壓分布，進(jìn)行三維空間的聲源識(shí)別定位，并使用剛性表面球形傳聲器陣列進(jìn)行實(shí)驗(yàn)驗(yàn)證。

本文首先在第1節(jié)中根據(jù)剛性球面發(fā)生聲波散射的Neumann邊界條件，建立聲波散射數(shù)學(xué)模型，進(jìn)而得出剛性球面散射前入射波聲壓與散射后總聲壓之間的數(shù)學(xué)關(guān)系。通過(guò)仿真，檢驗(yàn)球面入射波遇到剛性表面球形傳聲器陣列時(shí)發(fā)生散射后的聲場(chǎng)變化規(guī)律。然后在第2節(jié)中，將檢驗(yàn)由球形傳聲器陣列測(cè)量的散射后總聲壓，通過(guò)球面近場(chǎng)聲全息方法，重構(gòu)入射波形成的聲場(chǎng)，以及參數(shù)對(duì)重構(gòu)精度的影響。最后在第3節(jié)中給出結(jié)論。

1 剛性球體表面聲波散射的數(shù)學(xué)模型

1.1 剛性表面球體引起聲波散射的數(shù)學(xué)模型

在各向同性聲傳播介質(zhì)中，入射聲波可由Helmholtz方程的解的形式表示為：

式中：（r，θ,?）為聲場(chǎng)中任意一點(diǎn)的球坐標(biāo)，pi為此點(diǎn)處入射波的聲壓，ω表示角頻率，k為波數(shù)（k=ω/c，c為聲音在空氣介質(zhì)中的傳播速度），t表示時(shí)間。jn(kr)為球貝塞爾函數(shù)，表示聲場(chǎng)徑向即r方向的解；為球諧函數(shù)，表示聲場(chǎng)θ和?方向的解；Amn為波動(dòng)方程特解的系數(shù)。對(duì)于穩(wěn)態(tài)聲場(chǎng)，可將時(shí)間依賴(lài)性因子e-iωt略去。

假設(shè)在聲場(chǎng)空間中放置一個(gè)半徑為a的剛性表面球體，其球心與坐標(biāo)系原點(diǎn)重合，則剛性球體表面會(huì)引起入射波的散射現(xiàn)象，并且球面所發(fā)出散射聲波的頻率與入射波的頻率相同。聲場(chǎng)總聲壓 pt由入射波的聲壓pi與散射波的聲壓ps疊加而成：

剛性球體表面向外輻射的散射聲波可表示為：

其中hn(kr)為第一類(lèi)球漢克爾函數(shù)，Cmn為球諧函數(shù)的系數(shù)。

由于剛性球體表面不會(huì)吸收入射波的能量，故將其完全散射出去。同時(shí)，剛性球面也不會(huì)由于散射而變形，因此，球面處的法向介質(zhì)粒子速度即總聲壓的梯度為零，剛性球面發(fā)生散射的Neumann邊界條件［7-9］可由下式表示：

其中，ut(a)表示剛性球面上的法向介質(zhì)粒子速度，將其結(jié)合于公式（1）與公式（3）可得出：

式中：‘′’為求導(dǎo)符號(hào)。從而散射波的聲壓ps，可以表示為：

于是，根據(jù)公式（1）可得到聲場(chǎng)總聲壓（pt）為：

在這里，根據(jù)剛性球面的Neumann邊界條件得出散射波與入射波的關(guān)系，從而進(jìn)一步得出經(jīng)過(guò)剛性表面球體散射之后的總聲場(chǎng)模型。對(duì)于鑲嵌在剛性表面球形陣列上的傳聲器，在應(yīng)用中測(cè)得的是散射后的聲場(chǎng)總聲壓，根據(jù)式（1）與式（7）即可計(jì)算出發(fā)生散射前的入射波聲場(chǎng)聲壓。

1.2 剛性表面球體引起聲波散射后聲場(chǎng)的變化

對(duì)于剛性表面球形傳聲器陣列，各傳聲器都鑲嵌在球體的剛性表面上，測(cè)量與記錄聲場(chǎng)的聲壓。因此，研究球形傳聲器陣列表面處發(fā)生的散射現(xiàn)象，對(duì)于運(yùn)用球形傳聲器陣列測(cè)量聲壓，進(jìn)行近場(chǎng)聲全息或其它用途的聲學(xué)計(jì)算分析具有重要意義。本小節(jié)對(duì)球面入射波聲場(chǎng)進(jìn)行分析，先給出球面波聲場(chǎng)聲壓的球諧函數(shù)展開(kāi)形式的解析公式，然后檢驗(yàn)球面入射波在遇到剛性表面球形傳聲器陣列時(shí)，散射發(fā)生前后球體表面處聲壓分布的變化。

假設(shè)球面波由（rs，θs，?s）處單極子點(diǎn)聲源發(fā)出，則空間中任意一點(diǎn)（r，θ，?）處的聲壓可以表示為［13］：

其中符號(hào)‘*’表示復(fù)數(shù)共軛，ρ為介質(zhì)的密度，Qs為球面波聲源強(qiáng)度，N為球諧函數(shù)的截止項(xiàng)數(shù)，在本文中取N=30。式中r＜表示r與rs兩者中取值較小者，r＞表示r與rs兩者中取值較大者。

圖1 剛性表面球形傳聲器陣列與球面入射波示意圖

圖2 球面入射波聲場(chǎng)中球形傳聲器陣列剛性表面（r=a）處聲壓分布圖。入射波為來(lái)自位于點(diǎn)（r0，θ，φ）=（0.4 m，π/2，π）發(fā)出的球面波。（ka=5，a=0.2 m）

設(shè)置陣列半徑a為0.2 m，球面波聲場(chǎng)由一個(gè)位于（rs，θs，?s）=（0.4 m，π/2，π）處的單極子聲源生成，無(wú)量綱頻率（dimensionless frequency）ka為5。在上述參數(shù)下，得出球面的聲壓分布如圖2所示。在空間中不存在障礙物時(shí)，虛擬球面（r=a）的聲壓分布如圖2（a）），聲壓值在正對(duì)聲源處最大，在背對(duì)聲源處最小，球面上聲壓值的變化較為平緩；而當(dāng)剛性表面球形傳聲器陣列存在時(shí)，球面聲壓分布如圖2（b）），在正對(duì)聲源的球面區(qū)域，總聲壓增大且有最大峰值，其它處的聲壓值則迅速減小，并且在背對(duì)聲源處的聲壓分布有一個(gè)突起。

2 入射波聲場(chǎng)的重構(gòu)

本節(jié)將剛性球體表面散射模型應(yīng)用于球面近場(chǎng)聲全息方法中，通過(guò)仿真和實(shí)驗(yàn)檢驗(yàn)由剛性表面球形陣列上傳聲器測(cè)量的散射后總聲壓，重構(gòu)由入射波單獨(dú)形成的聲場(chǎng)，以及參數(shù)對(duì)重構(gòu)精度影響。

2.1 球面近場(chǎng)聲全息方法

球面近場(chǎng)聲全息方法［1-5］，由布置于球形陣列上傳聲器所測(cè)得的聲壓，通過(guò)聲場(chǎng)逆向空間變換，重構(gòu)出處于近場(chǎng)空間或聲源結(jié)構(gòu)表面的聲學(xué)量分布。聲場(chǎng)重構(gòu)公式如下：

式中：prec(r,θ,?)為位置(r,θ,?)處的重構(gòu)聲壓，Nt為基函數(shù)截止項(xiàng)數(shù)，取值與傳聲器數(shù)目相關(guān)，A~mn為未知系數(shù)。

對(duì)于剛性表面球形傳聲器陣列，基于考慮剛性球面對(duì)入射波的散射影響的聲場(chǎng)模型，可重構(gòu)出散射前的入射波聲場(chǎng)。根據(jù)第1.1節(jié)中入射波聲場(chǎng)與總聲場(chǎng)的關(guān)系，即公式（7），可得出系數(shù)A~mn為：

式中：p(a,θi,?i)為第i個(gè)傳聲器坐標(biāo)處的聲壓值；wi為該點(diǎn)相應(yīng)的權(quán)重，由球面數(shù)值積分方法確定；M為傳聲器的數(shù)目。

2.2 入射波聲場(chǎng)重構(gòu)的仿真

上節(jié)中介紹了使用剛性表面球形傳聲器陣列測(cè)量的聲場(chǎng)總聲壓，重構(gòu)入射波聲場(chǎng)的球面近場(chǎng)聲全息方法和數(shù)學(xué)模型，本節(jié)將檢驗(yàn)各參數(shù)對(duì)聲場(chǎng)重構(gòu)精度的影響。

設(shè)置聲場(chǎng)環(huán)境由一個(gè)位于（r0，θ，?）=（0.5 m，π/2，π）處的單極子聲源生成；剛性表面球形傳聲器陣列位于坐標(biāo)系原點(diǎn)，陣列半徑a=0.2 m，傳聲器數(shù)目為M=64。球形陣列上各傳聲器在根據(jù)公式（9）測(cè)得總聲壓數(shù)據(jù)，引入信噪比為40 dB的隨機(jī)測(cè)量噪聲，再使用這些帶有隨機(jī)測(cè)量噪聲的總聲壓數(shù)據(jù)，運(yùn)用球面近場(chǎng)聲全息理論，進(jìn)行聲場(chǎng)重構(gòu)計(jì)算。由式（10）與式（11），可重構(gòu)出空間中任意一點(diǎn)處的入射波聲壓prec(r,θ,?)，再將其與相同位置的聲壓理論值pthe(r,θ,?)（由公式（8）得出）進(jìn)行比較。在本文中將重構(gòu)一個(gè)以球形陣列中心為球心的球面（a≤R≤rs，其中R為重構(gòu)球面的半徑）上的聲壓分布，在此球面上均勻的取出Q點(diǎn)，聲壓在這些點(diǎn)的歸一化最小二乘（normalized least-squares）重構(gòu)誤差［14］定義為ε：

式中：‖?‖2表示二范數(shù)。

圖3為在波數(shù)k=15下，球形陣列表面上重構(gòu)的入射波聲壓分布，其與入射波聲壓理論值的誤差為ε=0.28%?？梢?jiàn)，采用考慮了球面散射模型的球面近場(chǎng)聲全息方法，在陣列表面處重構(gòu)的聲場(chǎng)聲壓分布，可以達(dá)到非常高的重構(gòu)精度。

圖3 球面近場(chǎng)聲全息方法重構(gòu)出的球形陣列表面入射波聲壓分布（a=0.2 m，k=15，M=64，R=0.2 m SNR=40 dB）

由圖3的結(jié)果得出，球面近場(chǎng)聲全息方法可重構(gòu)出球形陣列表面處聲壓分布。而在實(shí)際應(yīng)用中，需要重構(gòu)球形陣列外部空間的入射波聲壓分布，以作為聲源識(shí)別定位的依據(jù)。圖4為波數(shù)分別為k=5、k=15、k= 25時(shí)，重構(gòu)誤差ε隨著重構(gòu)半徑R增大的變化。

圖4 在不同波數(shù)下聲壓重構(gòu)誤差隨重構(gòu)球面半徑的變化（a=0.2 m，M=64，SNR=40 dB）

由圖4可見(jiàn)，在考慮了剛性球體表面聲波散射的入射波聲場(chǎng)重構(gòu)計(jì)算中，當(dāng)聲壓測(cè)量信號(hào)的信噪比為SNR=40 dB和球形陣列的傳聲器數(shù)目M＝64時(shí)，重構(gòu)半徑R由0.2 m增大至0.25 m的范圍內(nèi)，在波數(shù)k分別為5、15、25的情形下，重構(gòu)誤差ε均在10%以下，這說(shuō)明在球形陣列外部較近的區(qū)域，聲場(chǎng)重構(gòu)能取得較高的精度。這是因?yàn)榍蛐侮嚵袦y(cè)得了足夠多的由重構(gòu)球面?zhèn)鞑?lái)的聲場(chǎng)近場(chǎng)信息。隨著重構(gòu)半徑R的增大，即重構(gòu)球面與聲源之間距離的減小，重構(gòu)誤差ε也隨之增大，特別是當(dāng)波數(shù)k較大時(shí)，重構(gòu)誤差ε隨重構(gòu)半徑R增大而急劇增大，比如在波數(shù)k為25的情形下，重構(gòu)半徑R由0.2 m增大至0.5 m時(shí)，重構(gòu)誤差由2.3%增大至90%，這是因?yàn)橹貥?gòu)半徑R越大，球形陣列采集到由重構(gòu)球面?zhèn)鞑?lái)的近場(chǎng)信息即高波數(shù)信息越少，因而重構(gòu)誤差ε越大。而當(dāng)波數(shù)k較小時(shí)，重構(gòu)誤差ε則隨重構(gòu)半徑R增大較緩慢地增大。如圖4所示，當(dāng)重構(gòu)半徑R由0.2 m增大至0.5 m范圍內(nèi)，在波數(shù)k為15時(shí)，重構(gòu)誤差由0.4%增大至69.2%，而在波數(shù)k為5時(shí)，重構(gòu)誤差僅由0.1%增大至23.7%，這是因?yàn)楫?dāng)波數(shù)較小時(shí)，即使重構(gòu)半徑R較大，球形陣列依然采集到了足夠的由重構(gòu)球面?zhèn)鞑?lái)的近場(chǎng)信息，因而重構(gòu)誤差ε相對(duì)較小。

圖5為波數(shù)k=15，重構(gòu)半徑R=0.5m，即聲源所在重構(gòu)球面處的入射波聲壓分布圖，圖5（a））為由球面近場(chǎng)聲全息方法重構(gòu)的入射波聲壓分布，圖5（b））為根據(jù)公式（12）得出的入射波聲壓理論值的分布。此情形下，由于波數(shù)k及重構(gòu)半徑R較大，因而重構(gòu)誤差ε較大，為69.2%。但是，依然可以通過(guò)重構(gòu)的入射波聲壓分布圖，明顯識(shí)別出聲源的方位為（π/2，π）。

上述研究結(jié)果表明，使用球面近場(chǎng)聲全息方法進(jìn)行入射波聲場(chǎng)的重構(gòu)，在剛性表面球形陣列較近的區(qū)域，有較高的重構(gòu)精度，在球形陣列表面上的重構(gòu)精度最高。當(dāng)聲源輻射含有較高的波數(shù)k的情形下，在球形陣列測(cè)量總聲壓的信噪比SNR和陣列傳聲器數(shù)目M不變時(shí)，可將球形陣列盡量置于聲源的附近，即縮短重構(gòu)半徑R，聲源近場(chǎng)的入射波聲場(chǎng)重構(gòu)的誤差ε，將會(huì)明顯減小。

2.3 入射波聲場(chǎng)重構(gòu)的實(shí)驗(yàn)驗(yàn)證

在實(shí)驗(yàn)中，使用剛性表面球形傳聲器陣列測(cè)量總聲壓，運(yùn)用球面近場(chǎng)聲全息方法，重構(gòu)球形陣列表面及外部空間區(qū)域的入射波聲場(chǎng)。實(shí)驗(yàn)布置如圖6所示，剛性表面球形傳聲器陣列上，各傳聲器處的聲壓，通過(guò)B&K-LANXI設(shè)備采集，并通過(guò)PULSE軟件進(jìn)行FFT處理得到聲壓頻域數(shù)據(jù)，進(jìn)而可進(jìn)一步進(jìn)行聲場(chǎng)重構(gòu)的計(jì)算。實(shí)驗(yàn)在全消聲室內(nèi)進(jìn)行，球形陣列的球心設(shè)為坐標(biāo)系中心；在球形陣列外部空間，布置一個(gè)揚(yáng)聲器作為聲源，揚(yáng)聲器與功放設(shè)備控制模塊相連接，發(fā)出指定頻率的聲音。

圖6 實(shí)驗(yàn)布置照片

圖5 聲源所在球面處聲場(chǎng)聲壓重構(gòu)值與理論值比較（a=0.2 m，k=15，M=64，R=0.5 m，SNR=40 dB）

2.3.1 實(shí)驗(yàn)I

在本實(shí)驗(yàn)中，使用一個(gè)剛性表面球形傳聲器陣列作為測(cè)量前端，球形陣列的半徑a為0.15 m，陣列上傳聲器的數(shù)目M為64。按照上述實(shí)驗(yàn)設(shè)置，將揚(yáng)聲器放置在距離陣列中心0.2 m處，聲源頻率f為800 Hz，即波數(shù)k為14.8，使用球形陣列測(cè)量聲場(chǎng)總聲壓，然后重構(gòu)球形陣列表面入射波聲壓分布，如圖7所示。

圖7 波數(shù)k為14.8時(shí)，由球面近場(chǎng)聲全息方法重構(gòu)出的球形陣列表面處入射波聲壓分布（陣列半徑a=0.15 m，聲源離陣列中心距離d=0.2 m，傳聲器數(shù)目M=64）

將各傳聲器位置處重構(gòu)聲壓與各傳聲器所測(cè)得聲壓進(jìn)行比較，結(jié)果如圖8所示，由公式（12）可得它們之間的偏差達(dá)20.2%，說(shuō)明聲場(chǎng)分布在經(jīng)過(guò)剛性球面散射后變化很大，因而分析由剛性球面引起的聲波散射現(xiàn)象，對(duì)由剛性表面球形傳聲器陣列重構(gòu)入射波聲場(chǎng)的研究有著重要意義。

圖8 各傳聲器位置處聲壓的測(cè)量值與入射波聲壓重構(gòu)值比較（陣列半徑a=0.15 m，聲源離陣列中心距離d=0.2 m，傳聲器數(shù)目M=64，波數(shù)k=14.8）

將重構(gòu)面半徑設(shè)定為R=0.2 m，根據(jù)傳聲器測(cè)量的總聲壓，運(yùn)用球面近場(chǎng)聲全息方法，可得出聲源所在球面處的入射波聲壓分布，如圖9所示，通過(guò)此聲壓分布圖可以明顯得出，聲源的方位為（1.5 rad，-1.1 rad）。

2.3.2 實(shí)驗(yàn)II

在本實(shí)驗(yàn)中，使用一個(gè)B&K公司的剛性表面球形傳聲器陣列［15］作為測(cè)量前端，陣列的半徑a為0.097 5 m，陣列上傳聲器的數(shù)目M為36。按照實(shí)驗(yàn)I的布置，作為聲源的揚(yáng)聲器放置于（0.13 m，1.6 rad，4.7 rad）處。在剛性表面球形陣列外部空間中，放置一個(gè)1/4英寸檢驗(yàn)傳聲器，位于（0.107 5 m，1.6 rad，5.2 rad）處，其方位對(duì)應(yīng)于球形陣列表面處編號(hào)為7的傳聲器，兩者相距1 cm，用于檢驗(yàn)入射波聲場(chǎng)重構(gòu)的精度。

圖9 球面近場(chǎng)聲全息方法重構(gòu)出的聲源所在球面處入射波聲壓分布（陣列半徑a=0.15 m，聲源離陣列中心距離d=0.2 m，傳聲器數(shù)目M=64，波數(shù)k=14.8）

實(shí)驗(yàn)步驟如下：首先分別設(shè)置聲源頻率f為200 Hz、500 Hz、800 Hz、1 500 Hz、2 500 Hz、4 000 Hz，由球形陣列與聲場(chǎng)中傳聲器進(jìn)行聲壓信號(hào)測(cè)量。然后，將球形陣列移開(kāi)，則聲場(chǎng)中只有聲源與檢驗(yàn)傳聲器，保持它們的位置不變，再次分別在聲源頻率f為200 Hz、500 Hz、800 Hz、1 500 Hz、2 500 Hz、4 000 Hz的情形下，由此檢驗(yàn)傳聲器測(cè)量入射波聲場(chǎng)的聲壓信號(hào)。由實(shí)驗(yàn)步驟可知，前一次球形陣列與檢驗(yàn)傳聲器，測(cè)量的信號(hào)為總聲壓值，而后一次檢驗(yàn)傳聲器的測(cè)量值為入射波聲場(chǎng)的聲壓值。根據(jù)第2節(jié)介紹的入射聲場(chǎng)重構(gòu)理論，可采用前一次球形陣列測(cè)量到的總聲壓值，重構(gòu)出檢驗(yàn)傳聲器位置的入射波聲場(chǎng)聲壓值，并與檢驗(yàn)傳聲器直接測(cè)量的入射波聲場(chǎng)聲壓減小對(duì)比，以檢驗(yàn)其聲場(chǎng)重構(gòu)的計(jì)算精度。在上述各實(shí)驗(yàn)頻率下，在檢驗(yàn)傳聲器位置處，直接測(cè)量總聲壓值、直接測(cè)量入射波聲場(chǎng)聲壓值以及運(yùn)用近場(chǎng)聲全息方法重構(gòu)得出的入射波聲場(chǎng)聲壓值，如圖10所示。

圖10 各頻率下，檢驗(yàn)傳聲器位置處總聲壓、入射波聲場(chǎng)聲壓和重構(gòu)的入射波聲場(chǎng)聲壓的比較（陣列半徑a=0.097 5 m，聲源離陣列中心距離d=0.13 m，傳聲器數(shù)目M=36）

由圖10可得，對(duì)于檢驗(yàn)傳聲器位置處，總聲壓測(cè)量值與入射波聲場(chǎng)聲壓測(cè)量值均隨著頻率f的增大而增大，在各頻率下總聲壓測(cè)量值大于入射波聲場(chǎng)聲壓測(cè)量值，這符合第1.2節(jié)中的分析，因?yàn)闄z驗(yàn)傳聲器位于聲源與陣列之間的聲場(chǎng)區(qū)域。由于檢驗(yàn)傳聲器位置處的總聲壓測(cè)量值與入射波聲場(chǎng)聲壓測(cè)量值分別為兩個(gè)實(shí)驗(yàn)步驟下對(duì)應(yīng)的測(cè)量值，因而兩者只能在幅值大小上進(jìn)行比較，實(shí)驗(yàn)結(jié)果為，在檢驗(yàn)傳聲器位置處，各頻率下散射后的總聲壓級(jí)較散射前的入射波聲場(chǎng)聲壓級(jí)要大約4 dB左右，這說(shuō)明聲源輻射的聲場(chǎng)在經(jīng)過(guò)球形陣列散射后發(fā)生了較大的變化。同時(shí)，由圖10可得，在各頻率下，對(duì)于檢驗(yàn)傳聲器位置處，聲場(chǎng)重構(gòu)得出的入射波聲場(chǎng)聲壓級(jí)及與測(cè)量的入射波聲場(chǎng)聲壓級(jí)很接近，誤差在2 dB以?xún)?nèi)，這說(shuō)明運(yùn)用基于散射聲場(chǎng)模型的近場(chǎng)聲全息方法，可以以一定的精度重構(gòu)出入射波聲場(chǎng)。其中，此誤差不僅僅是聲場(chǎng)重構(gòu)計(jì)算方法引起的，其它影響誤差的因素主要有傳聲器測(cè)量誤差以及檢驗(yàn)傳聲器坐標(biāo)誤差等。

上述兩組實(shí)驗(yàn)表明，采用剛性表面球形傳聲器陣列測(cè)得的總聲壓，基于考慮聲場(chǎng)散射的聲場(chǎng)模型和球面近場(chǎng)聲全息方法，可以重構(gòu)出入射波聲場(chǎng)，消除由于球形陣列剛性表面散射造成的計(jì)算入射波聲場(chǎng)的誤差。

3 結(jié)論

本文通過(guò)應(yīng)用剛性球體表面的聲波散射的數(shù)學(xué)模型，建立起聲場(chǎng)入射波聲壓與發(fā)生散射后總聲壓之間的關(guān)系。針對(duì)不同的波數(shù)，檢驗(yàn)了球面波遇到剛性表面球體發(fā)生散射后，球體表面及外部空間聲場(chǎng)聲壓分布的變化規(guī)律。通過(guò)仿真與實(shí)驗(yàn)檢驗(yàn)了由剛性表面球形陣列上傳聲器測(cè)得的總聲壓，運(yùn)用球面近場(chǎng)聲全息方法，在球形陣列表面及外部空間中入射波聲場(chǎng)重構(gòu)的精度。得出如下結(jié)論：①剛性球體表面處聲壓分布在散射前與散射后發(fā)生了較大的變化，因而在入射波聲場(chǎng)重構(gòu)時(shí)，需要在聲場(chǎng)建模中考慮剛性球面的聲波散射。②通過(guò)建立剛性球面對(duì)聲波散射的數(shù)學(xué)模型，并應(yīng)用于球面近場(chǎng)聲全息方法中，可由剛性表面球形傳聲器陣列測(cè)量的總聲壓，實(shí)現(xiàn)入射波聲場(chǎng)的重構(gòu)。在入射波聲場(chǎng)重構(gòu)中，重構(gòu)誤差會(huì)隨著重構(gòu)半徑的增大而增大，在球形陣列外部的較近區(qū)域，能取得較高的重構(gòu)精度。同時(shí)，在波數(shù)較大時(shí)，由于聲源近場(chǎng)信息衰減較快，重構(gòu)誤差隨重構(gòu)半徑增大而急劇增大；而在波數(shù)較小時(shí)，重構(gòu)誤差隨重構(gòu)半徑增大則較緩慢地增大，即使在重構(gòu)半徑較大的情形下，重構(gòu)誤差也相對(duì)較小。此外，當(dāng)重構(gòu)面為聲源所在球面時(shí)，即使在總的重構(gòu)誤差較大的情形下，通過(guò)重構(gòu)出的入射波聲壓分布圖，依然可以很好地識(shí)別出聲源的方位。

本文研究在多個(gè)方面還需要進(jìn)一步研究：①本文以完全剛性球體的邊界條件，得出了聲波散射的聲場(chǎng)數(shù)學(xué)模型。但在實(shí)際中，剛性表面球形傳聲器陣列的球面上鑲嵌有多個(gè)傳聲器，在傳聲器位置波動(dòng)方程解并不滿(mǎn)足剛性面散射的Neumann邊界條件。為提高入射波聲場(chǎng)的重構(gòu)精度，需要考慮建立更精確的剛性表面球形傳聲器陣列的散射邊界條件模型。②在入射波聲場(chǎng)的重構(gòu)中，由測(cè)量誤差或環(huán)境噪聲引起的微小誤差會(huì)在聲場(chǎng)重構(gòu)逆向變換的過(guò)程中放大。引入如Tikhonov法或Landweber迭代法等正則化方法，可進(jìn)一步提高入射波聲場(chǎng)的重構(gòu)精度，由于論文篇幅的限制，本文沒(méi)有進(jìn)行探討。

［1］Williams E G，Nicolas Valdivia，Peter C Herdic，et al.Volumetric Acoustic Vector Intensity Imager［J］.Journal of the Acoustical So?ciety of America，2006，120：1887-1897.

［2］Jacobsen F，Moreno G，Grande E F，et al.Near Field Acoustic Ho?lography with Microphones Mounted on a Rigid Sphere［C］//Pro?ceedings of Inter-Noise.2008.

［3］Williams E G，Takashima K.Vector Intensity Reconstructions in a Volume Surrounding a Rigid Spherical Microphone Array［J］.Jour?nal of the Acoustical Society of America，2010，127（2）：773-783.

［4］Jacobsen F，Moreno G，Grande E F，et al.Near Field Acoustic Ho?lography with Microphones on a Rigid Sphere［J］.Journal of the Acoustical Society of America，2011，129（6）：3461-3464.

［5］Lu H，Li M.Reconstruction of Sound Field Based on Near-Field Acoustic Holography with a Rigid Spherical Microphone Array［C］//Proceedings of the Hong Kong：ACOUSTIC 2012 Hong Kong Conference and Exhibition.2012：3258-3258.

［6］Hardin R H，Sloane N J A.McLaren’s Improved Snub Cube and other New Spherical Designs in Three Dimensions［J］.Discrete Computational Geometry，1995（15）：429-441.

［7］莫爾斯P M，英格德U.理論聲學(xué)［M］.楊訓(xùn)仁，呂如榆，譯.北京：科學(xué)出版社，1984.

［8］杜功煥，朱哲民，龔秀芬.聲學(xué)基礎(chǔ)［M］.第3版.南京：南京大學(xué)出版社，2012.

［9］張海瀾.理論聲學(xué)［M］.北京：高等教育出版社，2007.

［10］褚志剛，楊洋，蔣忠翰.波束形成傳聲器陣列性能研究［J］.傳感技術(shù)學(xué)報(bào)，2011，24（5）：665-670.

［11］丁浩，李春曉，金江明，等.可識(shí)別聲源深度的三維聲聚焦波束形成方法［J］.傳感技術(shù)學(xué)報(bào)，2013，26（2）：175-181.

［12］劉夢(mèng)然，張國(guó)軍，簡(jiǎn)澤明，等.管道內(nèi)檢測(cè)器聲定位技術(shù)研究［J］.傳感技術(shù)學(xué)報(bào)，2014，27（4）：500-504.

［13］Franz Zotter.Analysis and Synthesis of Sound-Radiation with Spherical Arrays［D］.Austria：Institute of Electronic Music and Acoustics University，2009.

［14］Lu Huancai.Reconstruction of Vibroacoustic Responses Using Helmholtz Equation Least-Squares Method［D］.Doctoral Disserta?tion，Detroit：Wayne State University，2007.

［15］Haddad K，Hald J.3D Localization of Acoustic Sources with a Spherical Array［J］.Journal of the Acoustical Society of America，2008，123（5）：3311-3311.

李敏宗（1988-），男，浙江工業(yè)大學(xué)博士研究生，主要研究方向?yàn)樵肼曉醋R(shí)別定位及傳聲器陣列信號(hào)處理，limin?zong_svlab@163.com；

盧奐采（1962-），女，教授，博導(dǎo)，美國(guó)聲學(xué)學(xué)會(huì)會(huì)員，美國(guó)機(jī)械工程師學(xué)會(huì)會(huì)員，2007年于美國(guó)Wayne State Universi?ty獲得機(jī)械工程博士學(xué)位，1996年于浙江大學(xué)獲得工學(xué)博士學(xué)位，主要從事3D傳聲器陣列理論、聲全息理論及聲場(chǎng)可視化和聲學(xué)圖像處理等方面的研究，huancailu@zjut.edu.cn。

Reconstruction of Incident Sound Field and Identification of Sound Source Using a Rigid Spherical Microphone Array*

LI Minzong1，LU Huancai1，2*，JIN Jiangming1
（1.Key Laboratory of Equipment&Manufacturing，Ministry of Education&Zhejiang Province，College of Mechanical Engineering，Zhejiang University of Technology，Hangzhou 310014，China；2.Zhejiang Key Laboratory for Signal Processing，Zhejiang University of Technology，Hangzhou 310014，China）

Rigid spherical microphone array is an ideal measurement front-end for three-dimensional sound field re?construction.However，as the rigid surface of the array scatters incident sound waves，the directly measured pres?sures and near-field acoustic holography approach based on the math model of free field could not be employed to reconstruct the incident sound field.In this paper，the math model of scattered sound field based on Neumann boundary condition was applied；the relationship between the total sound pressure after scattering and the incident sound pressure was developed and analyzed.Then the incident sound field was reconstructed by using the input da?ta of measured total sound pressure with a rigid microphone array.The variation of sound pressure distribution on the surface of rigid sphere was examined through simulations and experiments.And the reconstruction accuracy was examined when the reconstruction parameters，such as frequency and reconstruction radius，varied.The reconstruc?tion results show that the incident sound field can be reconstructed with certain accuracy with a rigid spherical mi?crophone array，based on spherical near-field acoustic holography with the math model of scattered sound field.

spherical microphone array；reconstruction of incident sound field；spherical near-field acoustic holog?raphy；sound wave scattering；rigid surface sphere

TB877.2

A

1004-1699(2015)10-1459-08

??7320Q；7810

10.3969/j.issn.1004-1699.2015.10.007

項(xiàng)目來(lái)源：國(guó)家自然科學(xué)基金項(xiàng)目（51205354，51275469）；浙江省國(guó)際科技合作專(zhuān)項(xiàng)項(xiàng)目（2013C14014）

2015-03-12 修改日期：2015-08-12