基于k?空間格林函數(shù)的近場聲全息濾波方法研究

張一澍　馬宏偉　王星　王浩添

摘 要： 基于Neumann邊界條件下的空間聲場變換主要采用k?空間格林函數(shù)法，為保證該算法的穩(wěn)定性與可靠性，聲場重構(gòu)過程必須采取波數(shù)濾波處理。針對固定截止波數(shù)不能適應(yīng)濾波要求的局限性，分析波數(shù)域內(nèi)的噪聲分布，對信噪比進(jìn)行估計(jì)，進(jìn)而選取合適的截止波數(shù)，有效減少了攜帶聲場細(xì)節(jié)信息的高波數(shù)成分的損失。數(shù)值仿真結(jié)果驗(yàn)證了該濾波方法的可行性和正確性。

關(guān)鍵詞： 近場聲全息； 格林函數(shù)； 波數(shù)域?yàn)V波； 截止波數(shù)； 信噪比

中圖分類號(hào)： TN911.73?34； TB532 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2017）21?0158?04

Study on near?field acoustic holography filtering method based on k?space Green function

ZHANG Yishu， MA Hongwei， WANG Xing， WANG Haotian

（College of Mechanical Engineering， Xian University of Science and Technology， Xian 710054， China）

Abstract： The k?space Green function method is used in space transformation of sound field under the Neumann boundary condition. In order to ensure the stability and reliability of the algorithm， the k?space filtering processing must be adopted in the sound field reconstruction process of the algorithm. Since the fixed cutoff wavenumber exists the limitation for filtering requirement， the noise distribution in k?space is analyzed， and the signal?to?noise ratio （SNR） is estimated to select the suitable cutoff wavenumber to effectively reduce the loss of the high wavenumber component with the detail information of the sound field. The numerical simulation results show that the filtering method is feasible and valid.

Keywords： NAH； Green function； k?space filtering； cutoff wavenumber； SNR

近場聲全息（NAH）是一種用于聲源識(shí)別定位和聲場可視化的先進(jìn)技術(shù)，受到國內(nèi)外學(xué)者的普遍關(guān)注[1?2]。該技術(shù)通過測量包含倏逝波的近場數(shù)據(jù)，不僅可以重建出聲源的表面聲壓和法向振速，而且還能對整個(gè)三維聲場中任意點(diǎn)處的聲壓、質(zhì)點(diǎn)速度、聲強(qiáng)、聲源輔射聲功率和遠(yuǎn)場指向性等一系列聲學(xué)量進(jìn)行重建和預(yù)測[3]。近年來形成了眾多的NAH重建算法，其中基于空間Fourier變換的NAH技術(shù)理論上易于理解，算法和實(shí)驗(yàn)容易實(shí)現(xiàn)，計(jì)算速度快，因此在目前的實(shí)際應(yīng)用中最為廣泛。在基于Neumann邊界條件下對聲場數(shù)據(jù)進(jìn)行重構(gòu)的過程中，現(xiàn)有的空間聲場變換的有限離散化算法主要采用k?空間格林函數(shù)法，它由文獻(xiàn)[4]首先提出。文獻(xiàn)[5]將文獻(xiàn)[4]算法運(yùn)用到基于Dirichlet邊界條件下的聲壓場重構(gòu)，獲得了很好的重構(gòu)效果。文獻(xiàn)[6]修正了該格林函數(shù)法在輻射圓周上的計(jì)算公式。文獻(xiàn)[7?9]利用在Neumann邊界條件下的k?空間積分格林函數(shù)對聲場進(jìn)行重構(gòu)，有效地改善了k?空間抽樣格林函數(shù)在輻射圓周的奇異性問題。然而在聲場重建過程中，高波數(shù)域的噪聲誤差會(huì)被逆?zhèn)鬟f因子[G-1N]（Neumann邊界條件下格林函數(shù)的Fourier變換的逆）成千上萬倍地放大，從而影響重建精度。因此必須采取波數(shù)域?yàn)V波措施，以去除高波數(shù)誤差成分對重建結(jié)果的影響[10]。針對該問題，學(xué)者們通常采用指數(shù)濾波器進(jìn)行波數(shù)域?yàn)V波，從而抑制誤差的放大，而對于指數(shù)濾波器中的關(guān)鍵參數(shù)——空間截止波數(shù)（kc）值僅憑借經(jīng)驗(yàn)公式選取。本文結(jié)合文獻(xiàn)[11]提出的[kc]估計(jì)方法，研究了基于k?空間積分格林函數(shù)在聲場重建過程中波數(shù)域低通濾波的問題，對點(diǎn)源聲場逆向重構(gòu)進(jìn)行了數(shù)值仿真，并給出相應(yīng)的結(jié)論。

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

用聲場邊界函數(shù)值表示穩(wěn)態(tài)單頻波動(dòng)聲場的積分形式解，當(dāng)點(diǎn)源都集中在某一封閉曲面[s]內(nèi)時(shí)，Helmholtz公式表示為：

[（s）ejkrr???ns-?s??nejkrrds=-4π?0] （1）

式中：[n]為封閉曲面的外法線方向；[?]是空間分布函數(shù)；[?0]為Helmholtz方程的解析解；[ejkrr]為所選輔助函數(shù)，其中[r]表示場點(diǎn)[o]到封閉曲面[s]的距離，該函數(shù)除在[r=0]點(diǎn)有奇點(diǎn)外，其他地方均滿足假設(shè)條件和波動(dòng)方程。自由場中的格林函數(shù)滿足下面的方程[12]：

[g（r，r）=4πr-r-1expjkr-r] （2）

式中：[r]代表聲源的位置；[r]代表場點(diǎn)的位置。式（2）為具有[1r]奇點(diǎn)形式的函數(shù)并且滿足Helmholtz方程。endprint

Helmholtz方程與第二類邊界條件構(gòu)成的定解問題叫做第二邊值問題或Neumann問題。對于式（1），第二類邊界條件是指[???n]在區(qū)域邊界上為給定函數(shù)。相應(yīng)地，該邊界條件下滿足式（2）和Neumann邊界條件的解稱為Neumann格林函數(shù)。根據(jù)“虛源法”和實(shí)際聲源的相位關(guān)系得到Neumann格林函數(shù)[13]：

[g（r，r）=14πr-rejkr-r+14πr-rejkr-r] （3）

式中[r]表示虛源的位置。當(dāng)封閉曲面[s]為平面時(shí)，有[R=r-r=r-r]。通過歐拉公式和式（3）可以得到Neumann邊界條件下實(shí)空間域下法向質(zhì)點(diǎn)振速?聲壓的格林函數(shù)：

[gup（x，y，z）=-jρckejkR2πR] （4）

當(dāng)[???n=0]時(shí)，對應(yīng)為Neumann邊界條件。對式（4）進(jìn)行二維空間Fourier變換，整理得到k?空間下法向質(zhì)點(diǎn)振速?聲壓格林函數(shù)：

[Gkup=ρck?exp（-jdzkz1）kz1-jρck?exp（dzkz2）kz2] （5）

以及k?空間聲壓?法向質(zhì)點(diǎn)振速格林函數(shù)：

[Gkpu=kz1?exp（-jdzkz1）ρck-jkz2?exp（dzkz2）ρck] （6）

其中：

[kz1=k2-（k2x+k2y）， k2≥k2x+k2y] （7）

[kz2=（k2x+k2y）-k2， k2

2 k?空間積分格林函數(shù)

基于Neumann邊界條件下的格林函數(shù)在輻射圓周上具有奇異性，這種奇異性會(huì)使得格林函數(shù)的幅值在輻射圓周上具有很大的跳變，從而影響到重建的精度。k?空間積分格林函數(shù)法通過格林函數(shù)在k?空間的積分值來改善函數(shù)在輻射圓周上的奇異性。

設(shè)全息面大小為[Lx×Ly，]測量點(diǎn)數(shù)為[Nx×Ny，]重構(gòu)距離為[dz，]根據(jù)采樣定理，在空間采樣間隔為[Δ]時(shí)，全息面聲壓角譜范圍為[-πΔ≤kx≤πΔ，][-πΔ≤ky≤πΔ。]記[kxmax=πΔ，][kymax=πΔ]。在波數(shù)域的點(diǎn)[kx，ky]附近的環(huán)形區(qū)域帶[k2r2≤k2r≤k2r1]上進(jìn)行積分求格林函數(shù)的平均值，以克服輻射圓周上的奇異性。記[kr=（k2x+k2y）12，]積分環(huán)帶內(nèi)徑[kr2=（k2x+k2y）12-Δkr，]外徑[kr1=（k2x+k2y）12+Δkr，]其中[Δkr=][2Δkx，]為環(huán)帶寬度的一半。容易證明，積分帶寬[d=2Δkr=22πLx，]全息面波數(shù)域半徑[D=kxmax=πΔ=][NxπLx，]顯然有[D?d。]為便于分析，將k?空間格林函數(shù)積分原理與下文提到的聲壓角譜分布繪制在一張圖中，如圖1所示。

積分分為三個(gè)部分，即輻射圓內(nèi)小于[kr2]的傳播波區(qū)域、輻射圓外大于[kr1]的倏逝波區(qū)域，以及傳播波和倏逝波混合的環(huán)帶區(qū)域。于是通過計(jì)算可以得到基于Neumann邊界條件下的k?空間積分格林函數(shù)：

[GN（kx，ky，z）=-2jρckejdzk2-k2r2-ejdzk2-k2r1（k2r2-k2r1）dz， k2>k2r2-2jρckedzk2r2-k2-ejdzk2-k2r1（k2r2-k2r1）dz， k2r1

3 波數(shù)域?yàn)V波

波數(shù)域?yàn)V波窗函數(shù)最早由Veronesi和Maynard提出，其算法簡單可操作性較強(qiáng)，應(yīng)用最為廣泛。該函數(shù)通過在截止波數(shù)處采取平滑處理，獲得了很好的濾波效果，其窗函數(shù)的表達(dá)式為：

[∏（kx，ky）=1-12ekrkc-1α，kr≤kc12e1-krkcα，kr>kc] （10）

式中：[kc]為空間截止波數(shù)；[kr=k2x+k2y；][α]為可調(diào)參數(shù)，表示濾波器阻帶上的衰減率。

3.1 截止波數(shù)的選擇

截止波數(shù)[kc]的取值決定了參與重建過程的全息面聲壓角譜范圍。角譜中高波數(shù)的倏逝波成分對應(yīng)聲場中隨距離迅速變化的細(xì)節(jié)信息。為了獲得高分辨率的重建圖像，在重建過程中需要盡可能多地包含有效的倏逝波，這就要求選取較大的[kc；]然而，由于倏逝波衰減迅速，過高波數(shù)的倏逝波若未被濾掉，很容易被各種噪聲誤差所淹沒，在重建過程中這些倏逝波連同各種誤差將被逆?zhèn)鬟f因子按指數(shù)規(guī)律急劇放大，從而產(chǎn)生較大的重建誤差，從這一點(diǎn)分析[kc]又不能取得太大。因此，需要確定一個(gè)合適的[kc，]在保證重建精度的前提下以獲得盡可能高的空間分辨率。[kc]的取值通常根據(jù)經(jīng)驗(yàn)公式確定：[kc=0.6πΔ，]公式中固定地將截止波數(shù)取為當(dāng)前采樣間隔下最大波數(shù)成分的60%，由于沒有考慮信噪比、全息面與聲源面間的距離以及聲波頻率等因素的影響，使用中經(jīng)常會(huì)出現(xiàn)選取不當(dāng)而導(dǎo)致濾波失效的情況。為了彌補(bǔ)經(jīng)驗(yàn)公式存在的不足，綜合考慮信噪比、全息面與聲源面間的距離以及聲波頻率等因素，文獻(xiàn)[11]給出一種[kc]的選取公式：

[kc=kq，min（kxmax，kymax）， kq≤min（kxmax，kymax）kq>min（kxmax，kymax）] （11）

其中：

[kq=k2+SNRln1020dz2≥k2x+k2y] （12）

根據(jù)全息面同聲源面的距離[dz，]按式（11）和式（12）即可求解出截止波數(shù)的值，利用該值進(jìn)行濾波處理可以使得重建出的任意波矢量的有用信號(hào)將不會(huì)被噪聲淹沒。然而在實(shí)際測量時(shí)，式（12）中信噪比SNR一般未知，下面通過全息面聲壓信號(hào)角譜進(jìn)行估計(jì)。

3.2 信噪比估計(jì)

估計(jì)全息面聲壓信號(hào)SNR，需要知道全息面聲場信息中有效信號(hào)成分與噪聲的分布情況。全息面聲壓角譜范圍如圖1所示，將其化為[Ω1，Ω2]和[Ω3]三個(gè)區(qū)域。[Ω1]對應(yīng)輻射源以內(nèi)的區(qū)域，即[kr

設(shè)全息面上包含誤差干擾的實(shí)際測量法向振速為[v（x，y，zH）]，全息面理論法向振速為[vt（x，y，zH）]，噪聲誤差成分為[ve（x，y，zH）]。[V（kx，ky，zH），][Vt（kx，ky，zH）]和[Ve（kx，ky，zH）]分別為[v（x，y，zH），][vt（x，y，zH）]和[ve（x，y，zH）]的離散空間Fourier變換。由上述分析可知，在高波數(shù)區(qū)域，[Vt（kx，ky，zH）]隨波數(shù)[kr]的增大迅速衰減，在[Ω3]區(qū)域有用信號(hào)[Vt（kx，ky，zH）]衰減殆盡，幾乎完全由噪聲占據(jù)，因此可以通過[Ω3]區(qū)域上的數(shù)據(jù)估計(jì)噪聲信號(hào)。結(jié)合歐拉公式，則信號(hào)與噪聲之比SNR可近似用下式進(jìn)行估計(jì)：

[SNR≈10lgρckVt（kx，ky，zH）22ρckVe（kx，ky，zH）22=10lgPH（Ω1Ω2）22PH（Ω3）22] （13）

式中：[·2]表示2范數(shù)；[PH（Ω1Ω2）]表示全息面聲壓角譜中位于[Ω1+Ω2]區(qū)域內(nèi)的波數(shù)成分，即有用信號(hào)；[PH（Ω3）]表示聲壓角譜位于[Ω3]區(qū)域內(nèi)的波數(shù)成分，即噪聲信號(hào)。直接求解[PH（Ω3）]相對比較困難，可以用全部區(qū)域內(nèi)聲壓信息[PH（kx，ky，zH）]去掉[PH（Ω1Ω2）]得到，故式（13）可改寫為：

[SNR≈10lgρckVt（kx，ky，zH）22ρckV（kx，ky，zH）-Vt（kx，ky，zH）22=10lgPH（Ω1Ω2）22PH（kx，ky，zH）-PH（Ω1Ω2）22] （14）

4 數(shù)值仿真

按照上述理論分析估計(jì)全息面聲壓信號(hào)信噪比。給出SNR估計(jì)算例：點(diǎn)源中心位于坐標(biāo)原點(diǎn)，全息面的大小為2 m×2 m；測量點(diǎn)數(shù)為[50×50]；振動(dòng)頻率1 000 Hz；重建面到聲源面的距離為[λ10]。圖2（a）表示全部區(qū)域內(nèi)聲壓信息，包含有用信號(hào)和噪聲信號(hào)；圖2（b）表示[Ω1+Ω2]區(qū)域內(nèi)聲壓信息，即有用信號(hào)。由式（14）計(jì)算獲得全息面聲壓信噪比為25.8 dB，結(jié)合式（11）和式（12）解得[kc]的值為51.2 rad/m。另一方面，根據(jù)經(jīng)驗(yàn)公式得到[kc]的值為47.1 rad/m，從而舍掉了部分有效的倏逝波信息。

分別利用基于k?空間格林函數(shù)的兩種重構(gòu)算法，按照上述測量環(huán)境進(jìn)行法向振速的逆向重構(gòu)，仿真結(jié)果如圖3和圖4所示。其中圖3為不加波數(shù)域?yàn)V波的重構(gòu)結(jié)果，圖4為根據(jù)上文計(jì)算得到的[kc]并代入濾波函數(shù)進(jìn)行濾波后的重構(gòu)結(jié)果。

由圖3和圖4可知，在法向振速的逆向重構(gòu)中，兩種格林函數(shù)算法下的聲場重建都獲得了較好的重建分辨率，而且使用k?空間積分格林函數(shù)明顯改善了由輻射圓上的奇異性引起的重建孔徑上的波動(dòng)，使得重建聲壓曲面變得平滑。但是k?空間格林函數(shù)的幅值在輻射圓之外呈指數(shù)增長，由圖3法向振速幅值在重構(gòu)孔徑邊緣處引起的突變說明了這一點(diǎn)。隨著重構(gòu)距離的進(jìn)一步增大，重構(gòu)孔徑邊緣誤差增加十分迅速，甚至?xí)⒅鞣邃螞]。采用濾波函數(shù)減小了波數(shù)區(qū)間周邊處[G-1N]的幅值，觀察圖4濾波后的結(jié)果，可見邊緣誤差明顯被抑制。

5 結(jié) 論

本文將k?空間格林函數(shù)法運(yùn)用到了Neumann邊界條件下法向質(zhì)點(diǎn)振速測量的聲場重構(gòu)中，由數(shù)值仿真的結(jié)果比較了k?空間抽樣格林函數(shù)法與k?空間積分格林函數(shù)法在聲場重構(gòu)中的異同點(diǎn)。在對重構(gòu)結(jié)果進(jìn)行波數(shù)域?yàn)V波時(shí)，信噪比的估計(jì)對于截止波數(shù)的選取起到了關(guān)鍵作用。考慮到重建過程中應(yīng)盡量保留倏逝波信息，因此未將較低波數(shù)區(qū)域內(nèi)的波數(shù)成分估計(jì)為噪聲。對信噪比進(jìn)行估計(jì)后，得到了較精確的截止波數(shù)，確定了參與重構(gòu)過程的全息面聲壓角譜范圍。選擇適宜的截止波數(shù)進(jìn)行濾波，能夠完善k?空間格林函數(shù)法的重構(gòu)效果，提高全息算法的抗噪聲干擾能力。該濾波方法的處理過程雖然采取了一些近似，但是簡單方便，其精度一般滿足重構(gòu)要求，可為進(jìn)一步的工程應(yīng)用提供參考。

參考文獻(xiàn)

[1] ZHU J， CHRISTENSEN J， JUNG J， et al. A holey?structured metamaterial for acoustic deep?subwavelength imaging [J]. Nature physics， 2011， 7（1）： 52?55.

[2] ZHOU Y， LU M H， FENG L， et al. Acoustic surface evanescent wave and its dominant contribution to extraordinary acoustic transmission and collimation of sound [J]. Physical review letters， 2010， 104（16）： 1?4.

[3] 張永斌，畢傳興，張小正.統(tǒng)計(jì)最優(yōu)近場聲全息重建精度和計(jì)算速度優(yōu)化方法[J].聲學(xué)學(xué)報(bào)，2014（2）：191?198.

[4] VERONESI W A， MAYNARD J D. Nearfield acoustic holography （NAH）Ⅱ. Holographic reconstruction algorithms and computer implementation [J]. Journal of the acoustical society of America， 1987， 81（5）： 1307?1322.

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

[6] 金莉萍.基于格林函數(shù)的典型聲場反演技術(shù)[D].哈爾濱：哈爾濱工程大學(xué)，2008.

[7] 邵光輝，牛悅嬌，馬佳男.基于K?空間積分格林函數(shù)的近場聲全息技術(shù)[J].赤峰學(xué)院學(xué)報(bào)，2012（9）：8?11.

[8] 馬佳男，楊德森，時(shí)勝國，等.基于Green函數(shù)的兩種STSF算法的比較[J].振動(dòng)與沖擊，2012，31（6）：155?159.

[9] 馬佳男.基于格林函數(shù)的近場聲全息技術(shù)[D].哈爾濱：哈爾濱工程大學(xué)，2012.

[10] 于飛，陳劍，周廣林，等.噪聲源識(shí)別的近場聲全息方法和數(shù)值仿真分析[J].振動(dòng)工程學(xué)報(bào)，2003，16（3）：339?343.

[11] 陳心昭，畢傳興.近場聲全息技術(shù)及其應(yīng)用[M].北京：科學(xué)出版社，2013.

[12] 梁昆淼.數(shù)學(xué)物理方法[M].北京：高等教育出版社，1995.

[13] VERONESI W A， MAYNARD J D. Digital holographic reconstruction of sources with arbitrarily shaped surfaces [J]. Journal of the acoustical society of America， 1989， 85（2）： 588?598.

[14] 辛雨，張永斌，畢傳興，等.基于空間傅里葉變換的平面近場聲全息中信噪比估計(jì)方法研究[J].計(jì)量學(xué)報(bào)，2010，31（6）：537?542.endprint