陳 璐， 郭世旭， 王月兵， 鄭慧峰， 徐遨璇

(中國計量大學(xué) 計量測試工程學(xué)院， 浙江 杭州 310018)

1 引 言

近場聲全息(near-field acoustical holography,NAH)是聲學(xué)研究中的一種極其重要的聲場推算方法，由Williams提出并根據(jù)空間聲場變換進(jìn)行大量實(shí)驗(yàn)應(yīng)用[1～3]。近場聲全息技術(shù)在噪聲源定位、聲場重構(gòu)及可視化方面有重要應(yīng)用。通過在輻射體的近場區(qū)域內(nèi)測量聲場數(shù)據(jù)重建和預(yù)測整個三維空間的聲學(xué)特性，從而識別和定位聲源，獲得聲源的聲輻射特性，達(dá)到判斷噪聲來源進(jìn)行聲源控制的目的。目前工程中基于近場聲全息的噪聲分析識別算法主要包括空間聲場變換聲全息(spatial transformation of sound field,STSF)，統(tǒng)計最優(yōu)聲全息(statistically optimized NAH,SONAH)，等效源法聲全息(equivalent source model,ESM)等算法。在算法的選擇和應(yīng)用上，需要根據(jù)實(shí)際需求以及應(yīng)用場合進(jìn)行選取。

在技術(shù)研究方面，首先興起的是STSF，Arteaga通過對全息聲場數(shù)據(jù)進(jìn)行線性預(yù)測邊界填充外推法，使STSF算法可在較小孔徑下仍具有理想的重建精度[4]。姬慶對比了不同格林函數(shù)對聲場重建的影響[5]。Jiang等通過結(jié)合正交球面波疊加法和數(shù)據(jù)擴(kuò)展法來改進(jìn)算法，提高聲場重構(gòu)的精度[6]。郭世旭結(jié)合二維平面聲壓構(gòu)建技術(shù)提出了高效率的聲場測量方法[7]。由于重建非穩(wěn)態(tài)聲場的需求，基于等效源的近場聲全息算法受到重視，Pereira在球型陣列的等效源法中引入權(quán)矩陣對函數(shù)進(jìn)行迭代求解，達(dá)到了提高低頻近場重建精度的目的[8]。王成利用等效源算法對高頻換能器聲場進(jìn)行有效分析[9]。陶文俊提出了一種用壓縮感知將聲強(qiáng)信息轉(zhuǎn)為稀疏信號的方法，進(jìn)而優(yōu)化該算法[10]。陳漢濤運(yùn)用經(jīng)驗(yàn)?zāi)B(tài)分解法與等效源法結(jié)合，使低信噪比條件下的弱聲源得以檢測[11]。統(tǒng)計最優(yōu)算法研究起步最晚，張永斌提出了一種平面波優(yōu)化方法，對傳遞函數(shù)進(jìn)行改進(jìn)，提升了重建精度[12]。趙報川通過多種波函數(shù)組合得到聲場傳遞函數(shù)提高了聲場重建的精度[13,14]。熊久鵬提出二維平面聲壓重構(gòu)技術(shù)來提高聲場測量效率[15]。

在實(shí)際測量中，近場聲全息算法的精度受多方面影響，主要有聲源頻率、重建距離、測量面采樣間距、正則化參數(shù)選取等，本文在這些研究的基礎(chǔ)上，將仿真和實(shí)驗(yàn)相結(jié)合，互相驗(yàn)證，分析總結(jié)各個算法受上述因素的影響規(guī)律，并對運(yùn)算效率進(jìn)行分析比較，為不同需求的水下近場聲全息算法選擇及參數(shù)選取提供了一定參考。

2 近場聲全息算法原理

2.1 基于空間聲場變換的聲全息技術(shù)

當(dāng)格林函數(shù)滿足聲源表面均勻邊界條件時，Helmholtz方程的積分形式如下表示：

(1)

式中：x是位置矢量；r是聲源和對應(yīng)點(diǎn)所在面的法向距離；下標(biāo)S表示聲源面的性能評價；p表示聲壓。如果聲源面和全息面共形，它們之間可以通過測量數(shù)據(jù)和格林函數(shù)的卷積運(yùn)算得到。波數(shù)域測量面聲壓通過和格林函數(shù)乘積后再經(jīng)由逆Fourier變換得到空間重建面聲場，計算公式如下：

p(kx,ky,zS)=F-1(p(kx,ky,zH)e-ikz(zH-zS))

(2)

式中：kx、ky和kz分別表示在x、y和z軸方向的波數(shù)；p(kx,ky,zS)與p(kx,ky,zH)分別表示全息面與測量面的聲壓；zH和zS分別表示全息面和測量面與聲源的距離；F-1表示傅里葉逆變換。

2.2 基于統(tǒng)計最優(yōu)的聲全息技術(shù)

基于統(tǒng)計最優(yōu)的近場聲全息技術(shù)是局部近場聲全息算法的一種，核心算法是通過賦予空間聲場的復(fù)聲壓權(quán)值來計算重建面的聲壓，自由聲場里任意點(diǎn)的復(fù)聲壓都可以用多個空間聲場中的傳播波和倏逝波疊加計算得到：

(3)

式中：p(rH,n)為測量面上的聲壓。復(fù)聲壓權(quán)重系數(shù)C(rS)只和位置有關(guān)，即只要重建面與全息面的相對位置關(guān)系不變，不同聲壓都可以用同樣的權(quán)重系數(shù)矩陣。用標(biāo)準(zhǔn)Tikhonov正則化求解出系數(shù)C(rS)代入式(3)即可得基于統(tǒng)計最優(yōu)的重建面聲壓計算公式：

p(rs)=pT(rH)(AHA+θ2I)-1AHα(rS)

(4)

式中：H表示矩陣的共軛轉(zhuǎn)置；I表示單位對角矩陣；θ表示正則化參數(shù)；A表示2個平面間的傳遞矩陣；α(rS)為距聲源rS處平面波組成的列向量。

2.3 基于等效源的聲全息技術(shù)

等效源近場聲全息技術(shù)的基本思想是用輻射場內(nèi)部等效源疊加替代該物體所產(chǎn)生的輻射場。單個等效源強(qiáng)度可定義為[16]：

(5)

式中：G表示自由場格林函數(shù)；q表示虛擬聲源。

為避免聲場重建過程中奇異值帶來的誤差，計算虛擬表面聲強(qiáng)時對聲源表面進(jìn)行反運(yùn)算，由測量面聲場信息倒推出等效源后再重建出真實(shí)的重建面聲場。計算公式如下所示：

pS=GVS(pH(GHV)-1)

(6)

3 正則化方法

3.1 Tikhonov正則化

將測量面的聲壓信息pH中的誤差假設(shè)為空間不相關(guān)、方差為r2的高斯噪聲εp，那么噪聲和重建聲壓有如下關(guān)系：

pH=HpS,pH=pH0+εp, var|pH|=r2

(7)

式中：var為方差。

標(biāo)準(zhǔn)Tikhonov正則化分解通過增加一個約束項，使得求解最小二乘解的病態(tài)方程變成非病態(tài)方程。

(8)

標(biāo)準(zhǔn)Tikhonov正則化解為：

C(rs)=(AHA+θ2I)-1AHα(rS)

(9)

式中：H表示矩陣的共軛轉(zhuǎn)置；I表示單位對角矩陣。

正則化參數(shù)θ有多種計算方式，主要有廣義交叉驗(yàn)證法、L曲線法、固定參數(shù)法等，一般情況下固定參數(shù)的計算公式為[17]：

(10)

式中：d表示全息面和重建面間的距離；SNR表示信噪比。

3.2 廣義交叉驗(yàn)證法(GCV)

GCV法的主要思想是當(dāng)在聲全息重建過程中缺失一個聲場信息時，用此時的聲場模型的正則化解預(yù)測缺失的聲場信息。在標(biāo)準(zhǔn)Tikhonov正則化求解中GCV函數(shù)可以用下式表示：

(11)

式中：Tr為矩陣的跡。GCV函數(shù)達(dá)到極小值時對應(yīng)的θ就是所求的正則化參數(shù)。

3.3 L-曲線法(L-curve)

(12)

式中：ρ′、ρ″、ζ=、ζ″分別為ρ和ζ的一階、二階導(dǎo)數(shù)。該函數(shù)的曲率達(dá)到極大值時對應(yīng)曲線的拐角系數(shù)為正則化參數(shù)θ。

4 仿真模型及仿真結(jié)果

首先建立有限元仿真模型，聲源放置于坐標(biāo)原點(diǎn)，測量面與重構(gòu)面處于z方向法線上，且平行于xy平面，兩個面大小相等，見圖1。聲源頻率設(shè)置為500 Hz，測量面距聲源0.3 m，重建面距聲源0.25 m，測量面和重建面的大小都為1 m×1 m，采樣點(diǎn)數(shù)為11×11個。仿真環(huán)境為水域，聲速c=1 500 m/s，水密度ρ=1 000 kg/m3，模型采用邊界元，因此忽略實(shí)驗(yàn)環(huán)境中的水面和池壁反射，為了精確解析整個聲學(xué)域的壓力梯度變化，使用二次單元進(jìn)行離散化處理，在進(jìn)行頻域仿真時，對其采用大小為λ/6的細(xì)化網(wǎng)格，通過有限元仿真獲得測量面是理想值。

圖1 聲場測量示意圖Fig.1 Schematic diagram of sound field measurement

4.1 聲源頻率對重建精度的影響

在上述有限元仿真模型的基礎(chǔ)上進(jìn)行參數(shù)調(diào)整：修改聲源頻率為300 Hz到1 500 Hz，分別計算3種近場聲全息算法隨頻率變化產(chǎn)生的重建誤差。

由圖2可以看出，在低頻區(qū)域(<700 Hz)，ESM算法的誤差低于6%，其重建精度較高；當(dāng)頻率大于700 Hz時，3種算法的重建誤差隨著頻率增大而線性增大，且SONAH>ESM>STSF；同時隨著頻率增大，3種算法的重建精度變化趨勢相對穩(wěn)定。

圖2 誤差隨頻率的變化曲線圖Fig.2 Variation curve of error with frequency

4.2 重建距離對重建精度的影響

對部分仿真條件進(jìn)行調(diào)整，聲源頻率為500 Hz，獲取距聲源0.26～0.53 m的聲場數(shù)據(jù)，計算測量面與重建面距離變化時3種近場聲全息算法的重建誤差。

由圖3可以看出，ESM與SONAH在0.3 m的范圍內(nèi)重建誤差較小，STSF在超過0.1 m后重建誤差呈指數(shù)型增長，結(jié)果發(fā)生嚴(yán)重失真?？傮w上，3個算法的重建誤差都隨著重建距離增大而增大，重建誤差STSF>SONAH>ESM。

圖3 誤差隨重建距離的變化曲線Fig.3 Variation curve of error with reconstruction of distance

4.3 采樣點(diǎn)數(shù)對重建精度的影響

修改仿真條件，聲源頻率為500 Hz，測量面和重建面分別距聲源0.3 m、0.25 m，分別以 5×5到19×19采樣點(diǎn)數(shù)對測量聲場進(jìn)行采樣，計算比較算法重建誤差變化。

由圖4可以發(fā)現(xiàn)，3種算法的重建誤差隨著采樣點(diǎn)數(shù)的增加而減小，且STSF>SONAH>ESM。當(dāng)采樣點(diǎn)數(shù)大于250后，采樣點(diǎn)數(shù)變化對重建精度的貢獻(xiàn)率基本不變。其中隨著采樣點(diǎn)數(shù)的變化，ESM與STSF的誤差浮動范圍為10%，SONAH的誤差浮動范圍為25%。

圖4 誤差隨采樣點(diǎn)數(shù)的變化曲線Fig.4 Variation curve of error with the number of sampling points

由于算法的運(yùn)算方式不同，聲場數(shù)據(jù)量對計算效率的影響會有所差異。聲場數(shù)據(jù)量的大小取決于采樣點(diǎn)數(shù)，通過選取不同采樣點(diǎn)數(shù)，對3種重建算法的重建運(yùn)算時間進(jìn)行計算分析。由于基于統(tǒng)計最優(yōu)的近場聲全息算法運(yùn)算時間較長，因此對各算法的運(yùn)算時間取對數(shù)后進(jìn)行對比。由圖5可以看出運(yùn)算時間：STSF> ESM>SONAH。其中SONAH算法由于測量面聲場信息的矩陣維度提升而大幅降低重建速度，因此在計算高維度聲場且對重建速度要求不高時可選用ESM算法，當(dāng)對速度有較高要求，對重建精度要求不高時可選用STSF算法。

圖5 運(yùn)算時間隨采樣點(diǎn)數(shù)的變化曲線Fig.5 Variation curve of operation time with the number of sampling points

5 實(shí)驗(yàn)與結(jié)果

5.1 實(shí)驗(yàn)系統(tǒng)搭建

為驗(yàn)證算法的重建效果，在開闊大型水域進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)環(huán)境水域深度足夠，湖面平穩(wěn)，無過往船只及人員干擾。測量系統(tǒng)如圖6所示，陣元為相幅一致的水聽器。

圖6 測量系統(tǒng)示意圖Fig.6 Schematic diagram of measurement system

水聽器陣列尺寸為3 m×2.875 m，由168個水聽器排列成7行24列。選用柱形換能器，高度為0.08 m，直徑為0.06 m，實(shí)驗(yàn)中可視為點(diǎn)聲源。固定水聽器陣列，換能器與水聽器陣列中心位置相對。首先激勵換能器工作，在開闊水域中形成聲場分布，然后利用水聽器陣列對水域聲場進(jìn)行測量，輸出信號同步呈現(xiàn)在示波器上，最終將數(shù)據(jù)讀取并存儲在電腦中。實(shí)驗(yàn)器材：低頻換能器、功率放大器、信號源、示波器、水聽器陣列、采集卡、上位機(jī)。

5.2 算法結(jié)果比較

首先從聲源頻率、重建距離、采樣點(diǎn)數(shù)及正則化方式4個方面采集實(shí)驗(yàn)數(shù)據(jù)，然后利用各個算法進(jìn)行聲場重建，最后分析比較，驗(yàn)證仿真結(jié)果。實(shí)驗(yàn)環(huán)境中信噪比接近60 dB，前3者比較中皆選取Tikhonov正則化固定參數(shù)法進(jìn)行重建。

為了便于算法重建效果比較，采用分貝誤差定義聲場重建誤差：

(13)

式中：ps(m,n)和ph(m,n)分別為理論面和重建面點(diǎn)(m,n)處的聲場信息。

5.2.1 根據(jù)頻率變化的重建精度比較

將距離換能器0.05 m處的聲場面作為重建理論面，距離為0.15 m處作為測量面。用換能器分別發(fā)射頻率為500 Hz、1 000 Hz、2 000 Hz的聲波，得到不同頻率下重建面和理論面的聲場信息。用STSF算法、SONAH算法與ESM算法進(jìn)行聲場重建，重建效果與重建誤差如圖7與表1所示。

圖7 3種算法在不同頻率下的重建效果Fig.7 The reconstruction effects of three algorithms at different frequencies

表1 不同頻率下的重建誤差Tab.1 The reconstruction error of each frequency dB

可以看出，在低頻區(qū)域，SONAH算法和ESM算法保留了更多的聲壓信息，對聲場信號有較好的還原效果，重建精度優(yōu)于STSF算法。ESM算法對聲場信息還原最好，但容易產(chǎn)生虛像，造成噪聲源定位偏差。STSF算法具有較好的聲源定位效果，但忽略了其他微弱的聲場信息，重建誤差較大，且受聲源頻率變化的影響較為明顯。當(dāng)聲源頻率增大時算法的重建精度逐漸降低，且ESM>SONAH>STSF。

5.2.2 根據(jù)距離變化的重建精度比較

將距離換能器0.05 m處聲場面作為重建面，發(fā)出頻率為1 kHz的聲波，分別測量距離為0.15 m、0.25 m、0.35 m處的聲場分布作為測量面。計算不同算法對應(yīng)各個重建距離的重建誤差，重建效果與重建誤差如圖8與表2所示。

圖8 3種算法在不同重建距離下的重建效果Fig.8 Reconstruction effects of three algorithms at different reconstruction distances

表2 不同重建距離下的重建誤差Tab.2 The reconstruction error of each distances dB

通過對比可以看出隨著重建距離增大，重建誤差都逐漸變大。SONAH與ESM具有較強(qiáng)的重建效果， STSF算法除了有較好的定位能力， 對聲壓信息還原較差，受重建距離變化的影響最大。ESM與STSF由于算法中格林函數(shù)受距離變化影響較大，從而在重建計算中產(chǎn)生虛像和旁瓣，因此誤差增加明顯，其中STSF算法尤為明顯。

5.2.3 根據(jù)采樣點(diǎn)數(shù)變化的精度比較

對水聽器陣列進(jìn)行調(diào)整，分析對比水聽器采樣點(diǎn)數(shù)為7×8個、7×12個、7×24個情況下的聲場重建效果。重建效果與重建誤差如圖9與表3所示。

圖9 3種算法在不同采樣點(diǎn)數(shù)下的重建效果Fig.9 The reconstruction effect of three algorithms at different number of sampling points

表3 不同采樣點(diǎn)數(shù)時的重建誤差Tab.3 The reconstruction error of each sampling points number dB

經(jīng)過比較可以看出隨著測量點(diǎn)數(shù)的減少， 聲場信息也在減少，各算法的重建效果區(qū)別變得不明顯。當(dāng)測量點(diǎn)數(shù)增多時，SONAH算法有比較理想的精度，隨著測量點(diǎn)數(shù)減少，ESM算法表現(xiàn)出較大的優(yōu)勢。

測量點(diǎn)數(shù)的變化造成聲場數(shù)據(jù)量增多，使矩陣運(yùn)算維度增大引起計算效率差異，為對比采樣率對算法的影響，對3種算法的計算效率進(jìn)行比較，不同采樣率下算法的運(yùn)算時間如表4所示。由表4可以看出STSF算法始終保持著高速的計算效率，而 SONAH算法由于涉及矩陣運(yùn)算，當(dāng)采樣率變大時聲場重建運(yùn)算效率降低，在3種算法中計算速度最慢。ESM算法速度較為穩(wěn)定，受采樣率的影響較小。

表4 不同采樣點(diǎn)數(shù)時的計算時間Tab.4 Operation time at different number of sampling points s

5.2.4 正則化參數(shù)比較

正則化的主要作用是減少噪聲信號對信號源分析的影響，對測量面的聲場信息人為加入方差為36，均值為0的高斯噪聲，再用混入噪聲的測量面聲信號進(jìn)行聲場重建。

本次實(shí)驗(yàn)所用聲源頻率為1 000 Hz，測量面距聲源0.15 m，重建面距聲源0.05 m，采樣點(diǎn)數(shù)選取 7×24個，選取不同正則化參數(shù)進(jìn)行聲場重建，效果對比如圖10所示。

圖10 3種正則化參數(shù)選取方法對重建效果的影響Fig.10 The reconstruction effect of three regularization parameter selection methods

可以看出，在3種正則化方法中廣義交叉驗(yàn)證法的降噪效果最好，L-曲線法對噪聲源的還原能力最強(qiáng)。SONAH算法與L-曲線法結(jié)合的近場聲全息既能有效保留聲場信息，又能對噪聲進(jìn)行抑制，在該工況下具有最好的重建效果。

6 結(jié) 論

通過聲源頻率、重建距離、傳感器陣列的采樣點(diǎn)數(shù)及正則化參數(shù)的選取這4種對聲場重建影響較大的因素分析3種常用的近場聲全息算法在不同工況下的優(yōu)劣，分別對各自算法進(jìn)行了仿真趨勢分析和實(shí)驗(yàn)重建效果驗(yàn)證，得到以下結(jié)論：

1) 隨著聲源頻率降低，重建距離減小，采樣率增加，算法的重建精度提高；

2) STSF算法速度最快，對聲源有較好的定位效果，然而聲場重建精度受重建距離影響明顯，微弱信號聲場信息丟失嚴(yán)重，因此該算法適用于對噪聲定位具有高速且準(zhǔn)確定位的應(yīng)用場合，如海洋軍事與航天航空領(lǐng)域雷達(dá)定位、水域環(huán)境探測等；

3) SONAH算法重建效果良好，且能保留較多聲場信息，配合L-曲線正則化法后對噪聲能進(jìn)行有效抑制，但受測量點(diǎn)數(shù)影響較大，計算速度緩慢，適用于對計算速度要求較低且需要聲場有效還原的應(yīng)用領(lǐng)域，如房間聲場噪聲分析等；

4) ESM算法有較好的重建精度和速度，結(jié)合正則化法后對聲場信息保留情況較好，但在計算過程中對聲源位置有預(yù)測要求，因此適用于近距離下，聲源位置誤差較小的非接觸式噪聲分析，如機(jī)械結(jié)構(gòu)噪音定位等。