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

張一澍　馬宏偉　王星　王浩添

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

關鍵詞： 近場聲全息； 格林函數； 波數域濾波； 截止波數； 信噪比

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

1 Neumann邊界條件下的格林函數

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

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

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

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

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

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

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

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

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

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

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

以及k?空間聲壓?法向質點振速格林函數：

[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?空間積分格林函數

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

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

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

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

3 波數域濾波

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

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

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

3.1 截止波數的選擇

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

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

其中：

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

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

3.2 信噪比估計

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

設全息面上包含誤差干擾的實際測量法向振速為[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變換。由上述分析可知，在高波數區域，[Vt（kx，ky，zH）]隨波數[kr]的增大迅速衰減，在[Ω3]區域有用信號[Vt（kx，ky，zH）]衰減殆盡，幾乎完全由噪聲占據，因此可以通過[Ω3]區域上的數據估計噪聲信號。結合歐拉公式，則信號與噪聲之比SNR可近似用下式進行估計：

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

式中：[·2]表示2范數；[PH（Ω1Ω2）]表示全息面聲壓角譜中位于[Ω1+Ω2]區域內的波數成分，即有用信號；[PH（Ω3）]表示聲壓角譜位于[Ω3]區域內的波數成分，即噪聲信號。直接求解[PH（Ω3）]相對比較困難，可以用全部區域內聲壓信息[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 數值仿真

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

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

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

5 結 論

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

參考文獻

[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] 張永斌，畢傳興，張小正.統計最優近場聲全息重建精度和計算速度優化方法[J].聲學學報，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] 陳曉東.近場平面聲全息的測量和重構誤差研究[D].合肥：合肥工業大學，2004.

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

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

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

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

[10] 于飛，陳劍，周廣林，等.噪聲源識別的近場聲全息方法和數值仿真分析[J].振動工程學報，2003，16（3）：339?343.

[11] 陳心昭，畢傳興.近場聲全息技術及其應用[M].北京：科學出版社，2013.

[12] 梁昆淼.數學物理方法[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].計量學報，2010，31（6）：537?542.endprint