基于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方程。
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<k2x+k2y]
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<k2<k2r2-2jρckedzk2r2-k2-edzk2r1-k2（k2r2-k2r1）dz，
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<k]的波数区域，该区域内通常包含较大能量，其波数成分为幅值不随传播距离衰减的传播波；[ω2]对应两个同心圆之间的环状部分，即[k <kr<km]的波数区域，记[km=maxkxmax，kymax]，该区域内波数成分为倏逝波；[ω3]对应图1中四个边角部分，即[km<kr]的波数区域，其上的波数成分同样为倏逝波，且具有比[ω2]内波数成分更高的波数，在传播过程中衰减得更为迅速。通常情况下，[ω3]内的倏逝波几乎全部衰减到了全息测量系统的噪声水平以下，无法再用于重建。因此噪声误差在[ω3]内占主导地位，可以通过[ω3]内的倏逝波成分来估计全息面声压中噪声能量的大小[14]。
摘 要： 基于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方程。
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<k2x+k2y]
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<k2<k2r2-2jρckedzk2r2-k2-edzk2r1-k2（k2r2-k2r1）dz，
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<k]的波数区域，该区域内通常包含较大能量，其波数成分为幅值不随传播距离衰减的传播波；[ω2]对应两个同心圆之间的环状部分，即[k <kr<km]的波数区域，记[km=maxkxmax，kymax]，该区域内波数成分为倏逝波；[ω3]对应图1中四个边角部分，即[km<kr]的波数区域，其上的波数成分同样为倏逝波，且具有比[ω2]内波数成分更高的波数，在传播过程中衰减得更为迅速。通常情况下，[ω3]内的倏逝波几乎全部衰减到了全息测量系统的噪声水平以下，无法再用于重建。因此噪声误差在[ω3]内占主导地位，可以通过[ω3]内的倏逝波成分来估计全息面声压中噪声能量的大小[14]。
