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

張一澍　馬宏偉　王星　王浩添

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

以及k?空間聲壓?法向質(zhì)點振速格林函數(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ù)在輻射圓周上具有奇異性，這種奇異性會使得格林函數(shù)的幅值在輻射圓周上具有很大的跳變，從而影響到重建的精度。k?空間積分格林函數(shù)法通過格林函數(shù)在k?空間的積分值來改善函數(shù)在輻射圓周上的奇異性。

設(shè)全息面大小為[Lx×Ly，]測量點數(shù)為[Nx×Ny，]重構(gòu)距離為[dz，]根據(jù)采樣定理，在空間采樣間隔為[Δ]時，全息面聲壓角譜范圍為[-πΔ≤kx≤πΔ，][-πΔ≤ky≤πΔ。]記[kxmax=πΔ，][kymax=πΔ]。在波數(shù)域的點[kx，ky]附近的環(huán)形區(qū)域帶[k2r2≤k2r≤k2r1]上進行積分求格林函數(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所示。

積分分為三個部分，即輻射圓內(nèi)小于[kr2]的傳播波區(qū)域、輻射圓外大于[kr1]的倏逝波區(qū)域，以及傳播波和倏逝波混合的環(huán)帶區(qū)域。于是通過計算可以得到基于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ù)域濾波

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

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

3.2 信噪比估計

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

設(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ū)域有用信號[Vt（kx，ky，zH）]衰減殆盡，幾乎完全由噪聲占據(jù)，因此可以通過[Ω3]區(qū)域上的數(shù)據(jù)估計噪聲信號。結(jié)合歐拉公式，則信號與噪聲之比SNR可近似用下式進行估計：

[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ù)成分，即有用信號；[PH（Ω3）]表示聲壓角譜位于[Ω3]區(qū)域內(nèi)的波數(shù)成分，即噪聲信號。直接求解[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ù)值仿真

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

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

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

5 結(jié) 論

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

參考文獻

[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)計最優(yōu)近場聲全息重建精度和計算速度優(yōu)化方法[J].聲學(xué)學(xué)報，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é)報，2012（9）：8?11.

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

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

[10] 于飛，陳劍，周廣林，等.噪聲源識別的近場聲全息方法和數(shù)值仿真分析[J].振動工程學(xué)報，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].計量學(xué)報，2010，31（6）：537?542.endprint