聲全息聲壓場插值重構(gòu)方法研究

趙玉貴，劉志紅，儀垂杰，宿元亮

（1.機械與汽車工程學(xué)院青島理工大學(xué)，山東 青島 266520；2.教育部重點實驗室“工業(yè)流體節(jié)能與污染控制”青島理工大學(xué)，山東 青島 266520）

1 引言

聲全息技術(shù)是一種可同時識別聲源位置和傳播性能的有效手段，它利用聲源發(fā)出聲波的幅值和相位信息獲取聲場的聲壓空間分布，依據(jù)聲壓幅值大小識別聲源位置，聲位置矢量確定聲傳播性能[1]，其測量過程，如圖1所示。為提高聲壓的測量和重建精度，需在無限大平面上布置多個麥克風(fēng)完成，而實際測量過程只能在有限空間進行，即將無限大平面用有限采樣面截斷，從而導(dǎo)致測量及重建誤差，聲壓幅值精度下降[2]。如何在不增大測量工作量及經(jīng)濟成本情況下，有效減少測量重建誤差提高聲源識別精度是解決聲全息技術(shù)高效應(yīng)用的關(guān)鍵問題，同時該問題得到學(xué)者們的高度關(guān)注。

在眾多專家學(xué)者的研究下，Patch近場聲全息[3?8]以及插值法等多種方法應(yīng)運而生，其中，插值法因其計算便捷和擴展性備受青睞，2006年文獻[9]利用Hermite 插值對全息面聲壓重建大幅提高聲全息聲源的空間分辨率，但該方法除需要每個測點聲壓外還需測量每點的振速，在一定程度上增加了測量和計算的工作量。文獻[10]提出了一種基于帶限信號恢復(fù)算法的NAH分辨率增強方法，利用聲壓的波數(shù)域帶限特性，通過帶限信號恢復(fù)算法（PGA），實現(xiàn)全息面聲壓插值提高聲源識別精度，但該方法需多次迭代計算。2007年文獻[11]，2015年文獻[12]提出了基于正交球面波插值的NAH分辨率增強方法，以實際測量點數(shù)據(jù)為插值條件，通過若干不同階次的球面波源疊加擬合實際聲場實現(xiàn)全息面插值，在減少測量點數(shù)同時等效增加了全息面聲壓數(shù)據(jù)。不過這種方法因采用無限階正交球面波源疊加的有限截斷增加了實際擬合誤差，且插值過程需將球面波源布置在聲源內(nèi)部而往往聲源位置是無法預(yù)知，限制了該方法的實際應(yīng)用。為提高插值法的應(yīng)用性，2010年文獻[13]等提出虛源插值法，通過在全息面附近布置虛源面，利用實測全息面聲壓求得虛源面上簡單源源強，再由簡單源源強對全息面聲壓插值，利用插值數(shù)據(jù)進行聲壓重建。綜述所述，插值法能在不增加測量工作量條件下，有效提高聲源識別精度，但該方法仍存在急需解決的問題：（1）如何高效利用實際少量測量數(shù)據(jù)，高精度插值出任意位置上的聲壓幅值，且降低過程計算量；（2）如何減少測量點位置對聲壓插值的影響，提高全息面上聲壓重建的有效性。

提出一種融合牛頓插值和克希荷夫積分的新插值方法。該方法以少量實測點數(shù)據(jù)為對象，由二維牛頓插值法對虛擬位置聲矢量和聲壓進行插值，將離散的測量值轉(zhuǎn)換為連續(xù)函數(shù)值，再通過重采樣得到大量離散點，由克希荷夫積分得出克希荷夫面的聲勢，最后由上述測量及插值數(shù)據(jù)重構(gòu)全息面聲壓值。

2 聲全息測量與插值重構(gòu)

2.1 聲全息測量原理

聲全息聲源識別原理，如圖1所示。全息面H是麥克風(fēng)布置平面，重建面S是聲源面，Z0為全息面與重建面間距離，R為全息面麥克風(fēng)位置聲矢量，通過測量得到H面上聲壓分布函數(shù)U(x,y)，根據(jù)式（1）計算得出S上的聲壓分布函數(shù)U(a,Z)，由其幅值大小識別出聲源位置[14?15]。

圖1 聲全息測量原理圖Fig.1 Acoustic Hologram Schematic Represented by a Single Right?Angle Coordinate System

重建面S上的聲壓分布函數(shù)U(a,Z)，由式（1）得出：

式中：U(a,Z)—重建面S上任意一點Q(x,y,zS)的聲壓值；k0—波數(shù)；

C—全息系數(shù)為常值，簡化計算這里取1；

r1—重建源面上的任意點Q(x,y,zS)與全息面上的麥克風(fēng)位置點P(xH,yH,zH)間距離：

將式（2）代入式（1），且將U(a，Z)寫作US(x，y)整理：

這里：

r0=點到重建面上任意一點Q(x，y，zS)的距離。

由式（3）可知，精確確定全息面上的聲壓分布函數(shù)UH(x，y)是重建聲壓分布函數(shù)US(x，y)的關(guān)鍵。

2.2 插值重構(gòu)

為獲得全息面聲壓函數(shù)UH(x，y)，首先，保證測量麥克風(fēng)的測試數(shù)據(jù)的準確性即測得大量的聲源有效聲信息。陣列的合理構(gòu)型是保證陣列獲取大量有效信息的前提[16]。

這里采用同軸圓環(huán)陣列作為測量陣列，如圖2所示。其構(gòu)型方式為：第一環(huán)半徑R為入射波最小波長的一半即R=λmin2，同一環(huán)上麥克風(fēng)間距Δr，λmin4 ≤Δr≤λmin2。不同環(huán)間半徑差ΔR，0 ≤ΔR≤λmin2。其次，高效重建非測量位置的聲壓幅值。非測量位置點即全息面上聲壓待建點的數(shù)量直接影響全息面聲壓函數(shù)的重建精度。為此，以同軸圓環(huán)陣列面為全息面，聲壓待建點的確定沿用陣列麥克風(fēng)位置的設(shè)置原則，聲壓值由牛頓插值和克希荷夫積分實現(xiàn)。最后，對所有非測量位置點和測量點上的聲壓計算處理后得出UH(x,y)值帶入式（3）重建聲源面上的聲壓分布函數(shù)實現(xiàn)聲源識別。

圖2 極坐標同軸圓環(huán)全息面示意圖Fig.2 Polar Coordinate Ring Hologram

由Newton插值算法計算同環(huán)上的插值點聲壓和不同環(huán)上的聲壓待建點聲壓。首先，以單環(huán)陣列為測量陣，構(gòu)建同一環(huán)上的插值點位置的聲壓。取第1環(huán)R上的麥克風(fēng)測量聲壓數(shù)據(jù)及位置矢量θi為變量，構(gòu)建1 環(huán)上插值點的聲壓與位置矢量的函數(shù)關(guān)系，即任給一個位置矢量θ可得出該位置的聲壓幅值P。其P(θ)的各階差商計算過程：

函數(shù)P(θ)關(guān)于點θ0,θk的一階差商表示為：

將函數(shù)P(θ)關(guān)于點θ0,θ1,θk的二階差商表示為：

以此類推，得到函數(shù)P(θ)的k階差商：

根據(jù)差商的定義，把θ看作[0，π]上的一點，可得：

只要將（8）中后一式依次代入前一式，得到最終的插值公式：

由式（6），代入麥克風(fēng)測量聲壓數(shù)據(jù)，即位置矢量θi和聲壓值P(θi)，便可得出第一環(huán)R1上任意位置的聲壓幅值，其它環(huán)上的離散點聲壓可由克希荷夫積分式（10）計算得到，此時，設(shè)置第一環(huán)為虛擬克希荷夫面。

式中：S1—虛擬克希荷夫面；PS1—虛擬克希荷夫面S1上的聲壓函數(shù)—法線；G—自由空間格林函數(shù)；PSk—第k環(huán)上的聲壓函數(shù)。

PS1可由第一環(huán)上各插值點均值表示：

式中：Pi—插值后第i點處聲壓值，而ωi為第i點處的加權(quán)函數(shù)，表示各點對整個聲壓的貢獻大小，此處ωi=1n。而式（10）中的格林函數(shù)G則可以表示為：

式中：rk—聲源到測量點的距離。

為便于計算，將式（12）在球坐標系下展開[17]：

式中：jn(kr<)—球貝塞爾函數(shù)，其中，r<是取r,r0中的較小值；hn(kr>)—勒讓德函數(shù)，式中的r>是取r,r0中的較大值。將式（12）代入式（10）得

其次，采用同心圓環(huán)陣為測量陣列，即采用多個環(huán)上的麥克風(fēng)測量值為插值數(shù)據(jù)，提高測試數(shù)據(jù)的有效性和聲壓的重建精度。

其計算過程與單個環(huán)的一致，只是將克希荷夫面數(shù)量變?yōu)槿齻€，并將三個克希荷夫面對面外聲壓待建點的貢獻度進行加權(quán)計算。

最后，全息面重構(gòu)方法流程，如圖3所示。

圖3 全息面全息函數(shù)獲取流程圖Fig.3 Holographic Surface Holographic Function Acquisition Flow Chart

3 數(shù)值仿真與實驗分析

3.1 插值算法數(shù)值仿真分析

為探討上述方法的有效性和關(guān)鍵影響因素進行如下數(shù)值仿真。仿真數(shù)據(jù)：重建面S的中心點布置一點聲源，聲壓級*70dB，在平行于重建面且間距為Z=0.33m 和Z=2m 的位置分別布置一個3環(huán)同軸圓環(huán)陣，陣列中心與重建面中心等同高，環(huán)半徑差取為ΔR=0.3λ，麥克風(fēng)數(shù)量第1環(huán)8個，第2環(huán)16個，第3環(huán)20個。全息面為（1×1）m，布置9個環(huán)，聲壓待建點為243個，入射波頻率為2000Hz。仿真過程一圓環(huán)陣為測量陣列，陣列構(gòu)型及插值點，如圖4所示，全息面上的待建點位置，如圖5所示。

圖4 測量麥克風(fēng)陣列及插值點Fig.4 Microphone Measurement Array and Interpolation Points

圖5 全息面聲壓待建點Fig.5 The Point on the Holographic Surface where the Sound Pressure is to be Reconstructed

實驗1：首先，為驗證插值點聲壓重建精度，由陣列測量數(shù)據(jù)，經(jīng)插值和克希荷夫積分重構(gòu)計算出每個環(huán)上插值點聲壓幅值，然后，在相應(yīng)插值點位置放置麥克風(fēng)測量，將計算聲壓pin與實測聲壓比較，其誤差評價函數(shù)為：

大量仿真統(tǒng)計分析發(fā)現(xiàn)：聲波入射頻率、測量陣列構(gòu)型、麥克風(fēng)數(shù)量和聲壓待建點）數(shù)量是影響全息面上聲壓重建精度的主要因素。

因全息面上聲壓待建點布置方式是按陣列構(gòu)型中麥克風(fēng)位置的布置方式確定，所以關(guān)鍵因素包括同環(huán)上的麥克風(fēng)間距，它直接影響各單環(huán)插值點的牛頓聲壓插值精度，和相鄰環(huán)的環(huán)間距，它改變了聲壓待建點位置及數(shù)量對整個全息面的聲壓重建影響顯著。分別針對1000Hz，2000Hz，3000Hz三個頻段，不同麥克風(fēng)間距Δr和環(huán)間隔ΔR仿真，其結(jié)果，如表1，表2和圖6、圖7所示。

表1 不同頻率條件下各測點間隔對應(yīng)測量麥克風(fēng)數(shù)目Tab.1 The Number of Measuring Microphones Corresponding to Each Measuring Point Interval Under Different Frequency Conditions

表2 不同頻率條件下不同環(huán)間隔對應(yīng)的全息面聲壓待建點數(shù)量Tab.2 Number of Holographic Surface Acoustic Pressure Points to be Built Corresponding to Different Ring Intervals Under Different Frequency Conditions

圖6 不同頻率條件下，Δ與插值誤差關(guān)系Fig.6 The Relationship Between Δ and Interpolation Error Under Different Frequency Conditions

由圖6 可知，隨著麥克風(fēng)間距增大，插值誤差也在增大，當(dāng)d>λ3（Δλ>0.3）時的插值誤差變大速度快，d<λ3時誤差增速平穩(wěn)，尤其低頻1000Hz 時顯著。所以，麥克風(fēng)間距在λ4 ≤d≤λ3范圍選取效果較好，當(dāng)頻率高于1000Hz時，范圍適當(dāng)放寬到λ4 ≤d≤λ2。

由圖7可知，同軸圓環(huán)面上兩相鄰環(huán)間隔過小或過大均引起較大誤差，當(dāng)ΔR在λ/3~λ/2范圍內(nèi)，積分誤差相對較小，所以ΔR一般在該范圍內(nèi)取得。

圖7 不同頻率條件下，ΔR與積分誤差關(guān)系Fig.7 The Relationship Between ΔR and Integral Error Under Different Frequency Conditions

實驗二：其次，聲源重建面的重建結(jié)果。

3.2 插值算法有效性驗證

3.2.1 實驗方法

該全息測量實驗在半消聲室內(nèi)進行，室內(nèi)本底噪聲小于20dB，單環(huán)測量時實驗裝置實物布置，如圖9所示。

圖9 全息測量裝置圖Fig.9 Holographic Measuring Device Diagram

圖8 聲源面上重建聲源（聲源幅值歸一化）Fig.8 Reconstruction of Sound Source on the Sound Source Surface（Normalization of Sound Source Amplitude）

在源面上布置一個聲源，聲源選擇一個尺寸為（146×170×223）（單位mm）的音箱，發(fā)出頻率為f=2000Hz的單頻信號，音箱紙盆中心坐標為（0，0，0），在平行于源面且間距分別為0.33m 和2m的位置布置一測量面作為全息面，重建面距離測量面距離為0.1m，全息面與重建面尺寸均為（1×1）m。本實驗所用測量面為同軸圓環(huán)傳聲器陣列，陣列中心坐標分別為（0，0，0.33）和（0，0，2），平面內(nèi)最小環(huán)為第一環(huán)，直徑為D=0.17m，其他環(huán)間距為ΔR=6cm。在測量面第一環(huán)上均勻布置8個測量麥克風(fēng)，第二環(huán)及第三環(huán)上分別布置16和20個測量麥克風(fēng)（θ=0°除外），圓環(huán)中心布置1個參考麥克風(fēng)。

該實驗分別在不同距離上進行測量，為了對插值算法的有效性進行驗證，兩次測量均進行兩組對照實驗。

第一組實驗對三環(huán)所測數(shù)據(jù)采用Newton插值及克希荷夫積分處理得到非測量環(huán)的其他環(huán)上聲壓待建點的聲壓數(shù)據(jù)；第二組則采用測量麥克風(fēng)對全息面上所有環(huán)上的聲壓待建點進行實際測量測得聲壓數(shù)據(jù)。將得到的數(shù)據(jù)進行聲源聲壓重構(gòu)生成聲壓幅值圖，最后就重建面位置處的聲壓實測值與聲壓計算值進行比較分析。

3.2.2 實驗結(jié)果

（1）0.33m條件下重建結(jié)果

（2）2m條件下重建結(jié)果

在0.33m與2m距離條件下得到的重構(gòu)后聲壓幅值圖，如圖10、圖11所示。圖10（a）是0.33m距離情況下依據(jù)Newton插值結(jié)合克希荷夫積分公式得出的全息面聲壓數(shù)據(jù)重建出的聲源聲壓幅值圖；圖10（b）則是根據(jù)8環(huán)233個麥克風(fēng)實際測量值重建的聲壓。在2m距離情況下的聲壓幅值圖，如圖11所示。

圖10 0.33m處聲全息重建聲壓幅值圖Fig.10 Sound Pressure Amplitude Map of Acoustic Hologram Reconstruction at 0.33m

從圖10與圖11中可以看出，兩種方法重構(gòu)所得的聲壓幅值圖外形差別較?。辉龠M行數(shù)據(jù)比對，對重構(gòu)面上若干點聲壓進行實際測量，再與兩種方法的對應(yīng)點重構(gòu)的計算值進行比較，用式（15）算得各點重構(gòu)計算所得聲壓數(shù)值與實測值的相對誤差，誤差均小于5%，在可接受范圍之內(nèi)，可以用來替代原始數(shù)據(jù)進行重建，且使用Newton插值與克希荷夫積分在重建效果相差不大情況下計算則更為簡單。

圖11 2m處聲全息重建聲壓幅值圖Fig.11 Sound Pressure Amplitude Map of Acoustic Hologram Reconstruction at 2m

上述實驗結(jié)果同樣可以驗證上述插值方法能很好地實現(xiàn)全息面聲壓重建工作，能夠節(jié)約測量時間、簡化測量工作。

3.3 誤差分析

利用牛頓插值與克希荷夫積分公式相結(jié)合的復(fù)合算法實現(xiàn)聲全息聲壓場的重建的誤差主要存在于以下方面[18]：

這里主要采用插值算法來實現(xiàn)聲壓場的重建，主要是為了減少測量的時間及成本，但插值所得數(shù)值與實際數(shù)值存在一定的偏差，這種偏差是各種算法所存的問題，只能夠盡可能地減小這一偏差而不能消除，因此，應(yīng)合理選擇麥克風(fēng)間距及相鄰環(huán)間隔等參數(shù)。

而且，在進行插值時，測量所帶來的測量誤差也會影響聲壓待建點的數(shù)據(jù)，因此，為了減小測量誤差，可以多次測量，選取合理數(shù)值。

4 結(jié)論

為解決聲全息技術(shù)受測量空間的限制而無法得到無限大平面上準確聲壓問題，這里提出了牛頓插值與克希荷夫積分公式相結(jié)合的聲壓獲取方法。該方法通過對聲源聲矢量以及聲壓進行插值，將離散信息點轉(zhuǎn)化為連續(xù)點，然后重采樣出大量離散信息點構(gòu)建克希荷夫面的聲壓，再通過克希荷夫積分公式算得整個全息面上的聲壓。該方法可以在不增加實際測量點的情況下通過部分測量點的信息得到全息面上任意點處的聲壓信息，接著利用仿真對影響參數(shù)進行了選取，最后通過實測實驗對該理論的有效性及準確性進行驗證，證實該方法能夠通過理論計算的方式得到無限大面上的聲壓且誤差在允許的范圍之類，且計算較為簡單，節(jié)省了測量及計算時間，可以有效提高重建的效率。