基于Green函數(shù)的兩種STSF算法的比較

馬佳男，楊德森，時勝國，胡 博

(哈爾濱工程大學(xué) 水聲技術(shù)重點實驗室，哈爾濱 150001)

現(xiàn)有的基于Neumann邊界條件空間聲場變換的有限離散化算法主要采用了k-空間抽樣格林函數(shù)法對聲場進(jìn)行重構(gòu)，該算法直接對波數(shù)域下的格林函數(shù)進(jìn)行抽樣，從而得到波數(shù)空間格林函數(shù)的離散陣列。該算法雖然簡單，工程上易于實現(xiàn)，但是由于Neumann邊界條件下k-空間振速-聲壓格林函數(shù)在輻射圓周上存在奇異性，這種奇異性會使得格林函數(shù)的幅值在輻射圓周上具有很大的躍變，從而影響重構(gòu)精度。

1987年，Veronesi等人［1］首先提出實空間積分格林函數(shù)法，2004年，合肥大學(xué)的陳曉東［2］將該算法運用到了基于Dirichlet邊界條件下的聲壓場重構(gòu)，獲得了很好的重構(gòu)效果，并通過數(shù)值仿真證明該方法不僅適用于近場也同樣適用于遠(yuǎn)場。隨著矢量水聽器的發(fā)展，可以通過測量質(zhì)點振速來對聲場進(jìn)行重構(gòu)和預(yù)測［3-4］。由于介質(zhì)質(zhì)點振速方向與波達(dá)方向一致，所以振速提供了聲源方位的信息。在實際測量中更偏重于通過測量質(zhì)點振速來重構(gòu)聲場以獲得更好的重構(gòu)精度和定位信息。

本文將實空間積分格林函數(shù)法運用到了Neumann邊界條件下基于法向質(zhì)點振速測量的聲壓場重構(gòu)中，推導(dǎo)獲得Neumann邊界條件下該有限離散化算法的計算參數(shù)，通過對其傅里葉變換后的角譜求逆得到聲壓-法向質(zhì)點振速聲場重構(gòu)時格林函數(shù)的角譜，并由數(shù)值仿真給出了該算法與k-空間抽樣格林函數(shù)法在聲場重構(gòu)中基于不同重構(gòu)參數(shù)的誤差分析，從而比較兩種格林函數(shù)法在聲場重構(gòu)中的優(yōu)劣。

1 Neumann邊界條件下的格林函數(shù)

理想流體介質(zhì)中小振幅聲波傳播的波動方程［5］:

對于單頻聲波速度勢Φ=φe-jωt，ω為聲波的角頻率，φ是空間分布函數(shù)，它滿足Helmholtz方程:

其中，c為聲速，k為聲波波數(shù)，故波數(shù)空間又稱為k-空間。通過Helmholtz方程，可以推導(dǎo)出Helmholtz公式。此公式用聲場邊界函數(shù)值表示聲場(穩(wěn)態(tài)單頻波動聲場)的積分形式解。當(dāng)點源都集中在某一封閉曲面s內(nèi)時，Helmholtz公式表示為:

n為封閉曲面的外法線方向，φ0為Helmholtz方程的解析解。ejkr/r為所選的輔助函數(shù)，其中，r表示場點o到封閉曲面s的距離，此函數(shù)除在r=0點有奇點外，其它地方皆滿足假設(shè)條件和波動方程。

由式(3)可見，Helmholtz公式用φ和?φ/?n邊界值的面積分來確定聲場中任意一點的速度勢函數(shù)值，因此當(dāng)已知邊界質(zhì)點振速的分布和聲壓的分布值時，就可以用Helmholtz積分求出場中任意點的速度勢函數(shù)值。

格林函數(shù)表示一定邊界條件下點源的場，與邊界條件一一對應(yīng)。單頻聲場中的格林函數(shù)滿足下面的方程［6］:

其中，r'代表聲源的位置，r代表場點的位置，δ函數(shù)表示點源，時間因子取 e-jωt，解得:

式(5)表示的是自由場中的格林函數(shù)，稱為格林函數(shù)的基本解。該式為具有1/r奇點形式的函數(shù)且滿足Helmholtz方程，可以成為式(3)的輔助函數(shù)。利用格林函數(shù)這種可選性可以選擇適當(dāng)?shù)母窳趾瘮?shù)形式來簡化Helmholtz公式。

Helmholtz方程與第二類邊界條件構(gòu)成的定解問題叫做第二邊值問題或Neumann問題。對于式(3)，第二類邊界條件是指?φ/?n在區(qū)域邊界上為給定函數(shù)。相應(yīng)地，該邊界條件下滿足式(5)和Neumann邊界條件的解稱為Neumann格林函數(shù)。

根據(jù)“虛源法”，平面邊界下格林函數(shù)的數(shù)學(xué)表達(dá)式為:

其中，r″表示虛源的位置。式(6)第二項對應(yīng)著“虛源”，表示由平面邊界引起的反射波部分。在平面邊界絕對硬下，Helmholtz公式可以簡化成一項。由“虛源”和實際聲源的相位關(guān)系以及式(5)得到Neumann格林函數(shù)［7］:

gN即為Neumann格林函數(shù)。由歐拉公式:

其中uz(x，y，z)，p(x，y，z)分別為空間點(x，y，z)處為法向質(zhì)點振速、聲壓。

通過式(8)和式(9)可以得到Neumann邊界條件下，實空間域下法向質(zhì)點振速-聲壓(簡稱為振速-聲壓)的格林函數(shù)，:

當(dāng)?φ/?n=0時，對應(yīng)為 Neumann邊界條件，此時對式(3)進(jìn)行二維空間Fourier變換并利用歐拉公式，整理得到k-空間振速-聲壓的格林函數(shù):

以及k-空間聲壓-法向質(zhì)點振速的格林函數(shù):

其中:

2 實空間積分格林函數(shù)法

本文應(yīng)用的實空間積分格林函數(shù)法是對已經(jīng)離散化的實空間格林函數(shù)各子塊的中點(x0，y0)處進(jìn)行二元泰勒級數(shù)的展開，通過取不同階數(shù)的泰勒展開來近似連續(xù)光滑實空間格林函數(shù)在子塊上的積分值，即式(10)的積分值［8］。本文參考文獻(xiàn)［3］，修正了文獻(xiàn)［1］給出的計算參數(shù)，利用取泰勒展開式的前4階得到實空間積分格林函數(shù)的近似結(jié)果。

首先，令全息孔徑為Lx×Ly，將全息面離散為N×N矩陣，其中每個子塊的大小為 Δx×Δy，其中 Δx=Lx/N，Δy=Ly/N，于是有實空間格林函數(shù)的離散表達(dá)式:

其中，m1，m2分別代表x，y方向上離散格林函數(shù)的子塊序號，于是各子塊中心坐標(biāo)x0=m1Δx，y0=m2Δy，dz為重建面與全息面間的距離，這里簡稱為重構(gòu)距離。

對式(13)的各子塊中心做4階泰勒展開，通過整理得到的數(shù)學(xué)表達(dá)式為:

上式給出的系數(shù)ai，是通過泰勒四階展開得到的格林函數(shù)展開系數(shù):

上式b=kx0，d=ky0。令a=kR0，

式(14)是非中心子塊4階泰勒展開的結(jié)果，當(dāng)x0=0，y0=0且dz→0時R0→0，這時格林函數(shù)按照上式的計算方法具有奇異性，不能進(jìn)行直接的近似，所以在計算中心子塊積分值時將該子塊分成兩部分，即圓形區(qū)域和四個角區(qū)域。

于是，帶入式(10)所得圓形區(qū)域的積分可用極坐標(biāo)表示為:

然后，討論正方形中剩余的4個角域上的積分，因為格林函數(shù)的對稱性，所以這4個積分相同，只需求出其中一個即可，于是有:

將中心子塊的格林函數(shù)與非中心子塊格林函數(shù)相加，就得到了基于法向質(zhì)點振速測量Neumann邊界條件下的實空間積分格林函數(shù)。

令采樣點數(shù)N=64，測量孔徑為4 m×4 m，聲速c=1 500 m/s，對比在不同重構(gòu)距離下，k-空間抽樣格林函數(shù)與實空間積分振速-聲壓格林函數(shù)的幅值分別在k空間和實空間域上的分布特點:

圖1 兩種算法下振速-聲壓格林函數(shù)幅值分布圖Fig.1 The amplitude distribution of the two algorithms in the vibration velocity reconstructing the sound pressure

由圖1可以看出:

(1)振速-聲壓實空間積分格林函數(shù)在頻率和全息孔徑不變的條件下，隨著重構(gòu)距離的減小，其幅值變得尖銳，且孔徑中心與邊緣的幅值衰減急劇增加;而當(dāng)重構(gòu)距離保持不變時，該格林函數(shù)幅值隨著頻率的增加而變得陡峭;以上分布特點說明當(dāng)重構(gòu)距離小于0.05m或頻率增大時，由于實空間積分格林函數(shù)曲線的面積僅限于x=y=0附近很窄的范圍內(nèi)，從而導(dǎo)致該區(qū)域下采樣點數(shù)較少，從而影響重構(gòu)精度。

(2)k-空間抽樣格林函數(shù)在輻射圓周上有躍變，隨著重構(gòu)距離的增大，輻射圓周上的躍變變得愈發(fā)尖銳;同時，聲源頻率的增加同樣可以加劇輻射圓周上的躍變。這種躍變會使重構(gòu)面孔徑邊緣處的函數(shù)具有奇異性，增大重構(gòu)誤差，所以由此可以推斷該算法下的格林函數(shù)對重構(gòu)距離和聲源頻率比較敏感，在對聲場進(jìn)行重構(gòu)時會產(chǎn)生較大誤差。

3 數(shù)值仿真分析

3.1 振速－聲壓的聲場重構(gòu)

下面利用實空間積分格林函數(shù)法和k-空間抽樣格林函數(shù)法分別進(jìn)行振速-聲壓的聲場重構(gòu)。取與圖1相同的仿真條件，f=1 000 Hz，zH為全息面與聲源間的距離。

圖2 兩種算法基于法向質(zhì)點振速測量的聲壓場重構(gòu)Fig.2 The two algorithms reconstructing the sound pressure field based on measuring normal particle vibration velocity

圖2給出的是兩種算法下，法向質(zhì)點振速重構(gòu)聲壓場的等高線圖。從圖中可以看出:

(1)基于實空間積分算法的聲壓場重構(gòu)，雖邊緣有起伏，但該算法下得到的聲壓幅值具有良好的衰減特性。而當(dāng)dz和zH增大后，起伏也隨之增大，由于NAH應(yīng)用的條件之一是聲壓幅值在全息面邊緣上具有足夠大的衰減，因此誤差大的這些點對整個聲場的空間變換影響很小，而影響大的區(qū)域在孔徑中心，只要將該區(qū)域的幅值誤差控制在一定的范圍內(nèi)，就能最終保證重構(gòu)精度。

(2)k-空間抽樣格林函數(shù)法在孔徑邊緣不僅波動較大，且在dz、zH很小的情況下，“卷繞誤差”使孔徑邊緣產(chǎn)生了虛源;隨著dz和zH的增大，虛源強(qiáng)度增大從而影響對聲源特性的識別。即使利用漢寧窗和k域濾波，都不能抑制卷繞誤差對重構(gòu)結(jié)果的影響。

從式(9)可以看到，當(dāng)聲源頻率為聲速的整數(shù)倍時，k-空間格林函數(shù)的分子為零，存在奇異性，不能通過k-空間抽樣格林函數(shù)法對聲場進(jìn)行重構(gòu)。所以我們重點考察聲源頻率對實積分格林函數(shù)法和聲壓場重構(gòu)精度的影響。令dz=0.01 m，zH=0.1 m。

由圖(3)可以看出，實積分格林函數(shù)在低頻條件下重構(gòu)聲壓幅值具有較好的收斂性和衰減性;在高頻條件下，重構(gòu)聲壓幅值衰減率減小，且孔徑邊緣有起伏，重構(gòu)誤差增大。

圖3 不同聲源頻率下實積分格林函數(shù)法的重構(gòu)聲壓場Fig.3 Using the real integral Green algorithm to reconstruct the sound field with different frequencies

3.2 聲壓－法向質(zhì)點振速的聲場重構(gòu)

將經(jīng)過FFT變換得到實空間積分法向質(zhì)點振速-聲壓格林函數(shù)的角譜求逆，就可以得到Neumann邊界條件下，聲壓-法向質(zhì)點振速格林函數(shù)的角譜，通過該譜函數(shù)對法向質(zhì)點振速進(jìn)行聲場重構(gòu)。

下面利用實空間積分格林函數(shù)法和k-空間抽樣格林函數(shù)法分別對基于聲壓測量的聲場進(jìn)行聲壓-振速的聲場重構(gòu)。取與圖2相同的仿真條件，令dz=0.02 m，zH=0.15 m。

圖4 兩種算法基于聲壓測量的法向質(zhì)點振速場重構(gòu)Fig.4 The two algorithms reconstructing normal particle vibration velocity field based on measuring the sound pressure

從圖4可以看出相對于聲壓的幅值分布，法向質(zhì)點振速的幅值函數(shù)的收斂性更好，在孔徑中心較為尖銳。所以為了更好地比較兩種算法，圖4給出的是兩種算法的剖面圖。在k域濾波參數(shù)相同的條件下，實空間積分格林函數(shù)法孔徑邊緣波動較小，而k-空間抽樣格林函數(shù)法孔徑邊緣的波動較大。為了更直觀地對比兩種算法的優(yōu)劣，利用下式計算重構(gòu)聲場的累積誤差:

圖5分別給出了兩種算法下，不同的聲源頻率、重構(gòu)距離、全息面距離重構(gòu)聲場的幅值誤差變化曲線。其中:圖5(a)中zH=0.15 m，dz=0.02 m;圖5(b)中 dz保持0.02 m 不變;圖 5(c)中zH=0.15 m。

可以看出:① 相對于k-空間抽樣算法，實空間積分算法不僅可以得到較高的重構(gòu)精度，而且各重構(gòu)參數(shù)的累積誤差隨著參數(shù)的變化相對穩(wěn)定，說明該算法具有較強(qiáng)的適用性;② 圖5(c)虛線表示當(dāng)dz＞0.04 m時，k-空間抽樣算法對法向質(zhì)點振速場重構(gòu)失敗，因為孔徑邊緣的幅值出現(xiàn)的指數(shù)型增長已經(jīng)無法對聲源做出精確的識別。雖然由實空間積分格林函數(shù)在dz＞0.05 m處得到的累積誤差很大，但是其誤差主要集中在孔徑中心，孔徑邊緣處的誤差很小，幾乎不影響對聲源特性的判斷。

圖5 不同重構(gòu)參數(shù)對重構(gòu)精度的影響Fig.5 Different reconstruction parameters on the reconstruction accuracy

4 結(jié)論

通過數(shù)值仿真結(jié)果可以得出如下結(jié)論:

(1)當(dāng)通過測量法向質(zhì)點振速對聲壓場進(jìn)行重構(gòu)時，相對于k-空間抽樣法，實積分格林函數(shù)法有效地抑制了重構(gòu)中常存在的“卷繞誤差”，且該算法對聲源頻率不如k-空間抽樣格林函數(shù)敏感，雖然在高頻條件下，孔徑邊緣有波動，但沒有出現(xiàn)較大的誤差，對識別聲源的影響較小。說明在基于空間聲場變換的全息重構(gòu)中，由實空間積分得到的格林函數(shù)的角譜更理想并更好地改善了法向質(zhì)點振速-聲壓k空間格林函數(shù)在輻射圓周上的奇異性。

(2)振速-聲壓法向質(zhì)點振速進(jìn)行聲場重構(gòu)時，實空間積分格林函數(shù)法有效抑制了由于逆向重構(gòu)產(chǎn)生的“吉布斯”效應(yīng)，且各重構(gòu)參數(shù)的累積誤差隨著參數(shù)的變化相對穩(wěn)定，累積誤差相對較小。

(3)通過對比和分析，Neumann邊界條件下的實空間積分格林函數(shù)法相對于k-空間抽樣格林函數(shù)法具有更好的重建精度，更有效的識別聲源，而其應(yīng)用特點也可為進(jìn)一步的工程實踐提供參考。

［1］ Veronesi W A，Maynard J D.Nearfield acoustic holography(NAH)II.Holographic reconstruction algorithms and computer implementation［J］.J.Acoust.Soc.Am.，1987，81(5):1307-1322.

［2］陳曉東.近場平面聲全息的測量和重構(gòu)誤差研究［D］.合肥:合肥工業(yè)大學(xué)，2004.

［3］張永斌，畢傳興，陳 劍，等.基于質(zhì)點振速測量的近場聲全息技術(shù)［J］.農(nóng)業(yè)機(jī)械學(xué)報，2007，38(9):112-116.

［4］ Jacobsen F，Liu Y.Near field acoustic holography with particle velocity transducers［J］.J.Acoust.Soc.Am.，2005，118(5):3139-3144.

［5］何祚鏞，趙玉芳.聲學(xué)理論基礎(chǔ)［M］.北京:國防工業(yè)出版社，1981.

［6］梁昆淼.數(shù)學(xué)物理方法［M］.北京:高等教育出版社，1995.

［7］ Veronesi W A，Maynard J D.Digital holographic reconstruction of source with arbitrarily shaped surfaces［J］.J.Acoust..Soc.Am.，1989，85(2):588-598.

［8］ Wu S F.Methods for reconstructing acoustic quantities based on acoustic pressure measurements［J］.J.Acoust.Soc.Am.，2008，124(5):2680-2697.