基于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