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

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

1 引言

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

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

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

2 聲全息測量與插值重構

2.1 聲全息測量原理

聲全息聲源識別原理，如圖1所示。全息面H是麥克風布置平面，重建面S是聲源面，Z0為全息面與重建面間距離，R為全息面麥克風位置聲矢量，通過測量得到H面上聲壓分布函數(shù)U(x,y)，根據式（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)與全息面上的麥克風位置點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)的關鍵。

2.2 插值重構

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

這里采用同軸圓環(huán)陣列作為測量陣列，如圖2所示。其構型方式為：第一環(huán)半徑R為入射波最小波長的一半即R=λmin2，同一環(huán)上麥克風間距Δr，λmin4 ≤Δr≤λmin2。不同環(huán)間半徑差ΔR，0 ≤ΔR≤λmin2。其次，高效重建非測量位置的聲壓幅值。非測量位置點即全息面上聲壓待建點的數(shù)量直接影響全息面聲壓函數(shù)的重建精度。為此，以同軸圓環(huán)陣列面為全息面，聲壓待建點的確定沿用陣列麥克風位置的設置原則，聲壓值由牛頓插值和克希荷夫積分實現(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)陣列為測量陣，構建同一環(huán)上的插值點位置的聲壓。取第1環(huán)R上的麥克風測量聲壓數(shù)據及位置矢量θi為變量，構建1 環(huán)上插值點的聲壓與位置矢量的函數(shù)關系，即任給一個位置矢量θ可得出該位置的聲壓幅值P。其P(θ)的各階差商計算過程：

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

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

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

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

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

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

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

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

式中：Pi—插值后第i點處聲壓值，而ωi為第i點處的加權函數(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)上的麥克風測量值為插值數(shù)據，提高測試數(shù)據的有效性和聲壓的重建精度。

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

最后，全息面重構方法流程，如圖3所示。

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

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

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

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

圖4 測量麥克風陣列及插值點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ù)據，經插值和克希荷夫積分重構計算出每個環(huán)上插值點聲壓幅值，然后，在相應插值點位置放置麥克風測量，將計算聲壓pin與實測聲壓比較，其誤差評價函數(shù)為：

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

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

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

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

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

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

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

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

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

3.2 插值算法有效性驗證

3.2.1 實驗方法

該全息測量實驗在半消聲室內進行，室內本底噪聲小于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），平面內最小環(huán)為第一環(huán)，直徑為D=0.17m，其他環(huán)間距為ΔR=6cm。在測量面第一環(huán)上均勻布置8個測量麥克風，第二環(huán)及第三環(huán)上分別布置16和20個測量麥克風（θ=0°除外），圓環(huán)中心布置1個參考麥克風。

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

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

3.2.2 實驗結果

（1）0.33m條件下重建結果

（2）2m條件下重建結果

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

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

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

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

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

3.3 誤差分析

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

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

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

4 結論

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