0 引　言

近场声全息技术（Nearfield Acoustic Holography, NAH）由Maynard等在20世纪80年代提出^[1-2]，是声学研究中一种具有高分辨率重建精度的噪声源定位和声场可视化方法，在声场重构以及振动噪声控制等领域有着广泛应用。近场声全息通过多个传感器组成的测量面在被测声源物体表面的近场区域获取全息测量数据，建立重建面与测量面之间的空间变换关系，进而重建三维空间的声压场、声强场、质点振速场或者预测远场指向性，从而为声源识别定位、有效控制结构振动噪声提供技术参考^[3]。

在声场逆向重构过程中，位于高波数区的倏逝波成分将与噪声误差一起被逆传递因子按指数规律放大，导致重建结果产生很大误差，严重影响声场重建的准确性^[4]。针对这一问题，在重建过程中，必须对高波数区噪声误差进行滤波处理，以确保重建精度，为此国内外学者进行了广泛研究。Veronesi等^[5]提出一种在波数域加低通滤波窗的方法来抑制高波数成分波对噪声的放大效应，该方法采取指数滤波器对波数域噪声信号进行处理，其算法简单，可操作性强，但该方法没有考虑声波频率、信噪比以及测量面与源面间距离等因素的影响，使用中常会出现因截止波数取值不当以致滤波效果不理想的情况。赵晨等^[6]对指数滤波器中窗函数参数的合理设置进行分析，针对截止波数和窗函数陡度系数2个关键参数，总结了在不同声源频率下取得最优重建效果的截止波数和陡度系数的最优取值，对窗函数的设计优化具有重要意义。张德俊^[7]在系统存在噪声情况下，对指数滤波器的截止波数进行修正，同时还考虑到指数滤波器仅在测量有限距离时有效，随着测量距离增大，滤波窗收效甚微，提出一种带约束条件的最小二乘法滤波函数。罗禹贡等^[8]总结了倏逝波随距离变化的规律并设计了相应的空间频域滤波器。由于指数滤波器仅对被逆传递函数按指数规律放大的噪声误差有效，对其他区域噪声误差过滤效果不佳，何其伟等^[9]提出一种波数域循环加窗逼近法，通过循环加窗筛选最接近理论声压幅值的数据作为重建数据，从而确保声场重构的准确性。莫登沅等^[10-12]将空间域滤波和波数域滤波相结合，提出一种空间域即波数域联合迭代滤波窗，该滤波窗可对重建面边缘误差进行有效消除，并且经过多次迭代，明显提高了声场重建精度。在指数滤波器实际使用中发现，窗函数陡度系数的取值变化决定着窗函数的衰减速率，陡度系数取值越小，指数窗函数衰减越快，相较于截止波数的取值变化，陡度系数取值变化对反演效果影响较小，故其对指数滤波器衰减程度的调整也有限。此外，指数滤波器的陡度系数会同时影响通带内和高频带内滤波曲线的衰减速度，且不能时时兼顾两者之间的平衡，甚至会导致有价值的倏逝波被剔除，影响滤波效果，进而影响声场重建精度。

本文针对上述问题，选取2个脉动球源作为研究对象，采用巴特沃斯数字滤波器（Butterworth filter）作为滤波函数，对脉动球源声场反演过程进行数值仿真研究。基于重构误差最小原则，得到了使声场重建误差最小的巴特沃斯滤波器最优参数，并比较了该参数条件下，巴特沃斯滤波器与常用指数滤波器的效果。

1 基于Fourier空间声场变换的近场声全息

在均匀理想流体介质中，小振幅声波的三维波动方程为：

$ {\nabla ^2}p(x,y,z;t) = \frac{1}{{{c^2}}}\frac{{{\partial ^2}p(x,y,z;t)}}{{\partial {t^2}}} 。$

(1)

式中：${\nabla ^2}$为拉普拉斯算子；$p(x,y,z;t)$为质点声压函数；c为介质内的声速。对波动方程两边进行傅里叶变换，得到不依赖时间变量的单频声波场亥姆霍兹方程为：

$ {\kern 1pt} {\nabla ^2}p(x,y,z) + {k^2}p(x,y,z) = 0 。$

(2)

式中：$p(x,y,z)$为声场空间点$(x,y,z)$复声压；$k = \frac{\omega }{c}$为声波波数；$\omega = 2{\text{π}} f$为声源简谐振动的圆频率。当所有声源均位于$z = 0$平面以下，$ \mathrm{z} > 0 $$ z>0 $的空间为自由声场时，由Green公式可得任意平面上的声压为：

$ \begin{aligned} & \;{\kern 1pt} p\left( {{x_H},{y_H},{z_H}} \right) = \iint_S p\left( {{x_S},{y_S},{z_S}} \right){g_D}。\\ & \left( {{x_H} - {x_S},{y_H} - {y_S},{z_H} - {z_S}} \right){\rm{d}}{x_S}{\rm{d}}{y_S}。\end{aligned}$

(3)

式中：S为用于积分的无穷大平面；${\kern 1pt} p\left( {{x_{{H}}},{y_{{H}}},{z_{{H}}}} \right)$为测量面声压；$p\left( {{x_S},{y_S},{z_S}} \right)$为重建面声压；${g_D}( {x_H} - {x_S},{y_H} - {y_S},$${z_H} - {z_S} )$为满足Dirichlet边界条件的格林函数。其表达式为：

$ {g_D}(x,y,z) = \frac{{z(1 - {{i}}kR){{{e}}^{{{i}}kR}}}}{{2{\text{π}} {R^3}}} 。$

(4)

式中，$R = \sqrt {{{\left( {{x_{{H}}} - {x_{{S}}}} \right)}^2} + {{\left( {{y_{{H}}} - {y_{{S}}}} \right)}^2} + {{\left( {{z_{{H}}} - {z_{{S}}}} \right)}^2}}$。

对式（3）两边取空间二维Fourier变换，将空域卷积转化为波数域乘积可得：

$ P( {{k_x},{k_y},{z_H}} ) = P( {{k_x},{k_y},{z_S}} ){G_D}( {{k_x},{k_y},{z_H} - {z_S}} ) 。$

(5)

式中：$P( {{k_x},{k_y},{z_H}} )$为测量面声压$p( {{x_{{H}}},{y_{{H}}},{z_{{H}}}} )$的空间Fourier变换；$P( {{k_x},{k_y},{z_S}} )$为重建面声压$p( {{x_S},{y_S},{z_S}} )$的空间Fourier变换；${G_D}( {{k_x},{k_y},{z_H} - {z_S}} )$为Green函数${g_D}( {x,y,z} )$的空间Fourier变换。其解析表达式为：

$ {G_D}( {{k_x},{k_y},{z_H} - {z_S}} ) = {{{e}}^{{{i}}\left( {{z_H} - {z_S}} \right){k_z}}}。$

(6)

式中：${z_H}$为测量面z轴坐标，${z_S}$为声源面z轴坐标，${k_x}、{k_y}$分别为x轴和y轴方向波数，${k_z}$取值存在以下2种情况（$k_x^2 + k_y^2{\text{ = }}{k^2}$被称之为辐射圆，辐射圆区分了声波2种不同的传播模式）^[13]。

当$k_x^2 + k_y^2 \leqslant {k^2}$时：

$ {k_z} = \sqrt {\left( {k_x^2 + k_y^2} \right) - {k^2}}。$

(7)

${G_D}( {{k_x},{k_y},{z_H} - {z_S}} )$中包含了许多频率成分，在此情况下，${k_z}$为非负实数，${G_D}( {{k_x},{k_y},{z_H} - {z_S}} )$为波数空间的低频成分。在声波传播过程中，仅相位发生变化，幅值不随空间位置的变化而变化，称其为传播波。

当$ {k}_{x}^{2}+{k}_{y}^{2}>{k}^{2} $时：

$ {k_z} = \pm i\sqrt {\left( {k_x^2 + k_y^2} \right) - {k^2}}。$

(8)

当${k_z} = - i\sqrt {\left( {k_x^2 + k_y^2} \right) - {k^2}} $时，其对应的解没有物理意义，在无穷远处趋向于无穷大，故此时对${k_z}$取值为${k_z} = i\sqrt {\left( {k_x^2 + k_y^2} \right) - {k^2}} $。在此情况下，${k_z}$为纯虚数，${G_D}( {k_x}, {k_y},{z_H} - {z_S} )$为波数空间的高频成分，在声波传播过程中，其幅值在传播方向随距离变化按指数规律衰减，称其为倏逝波，为了获取到倏逝波中包含的重要信息，需在声源近场对其进行测量。

当已知测量面上声压分布$p\left( {{x_{{H}}},{y_{{H}}},{z_{{H}}}} \right)$时，对式（5）进行空间二维Fourier逆变换，可得重建面上声压为：

$ p\left( {{x_s},{y_s},{z_s}} \right) = {{\rm{F}}^{ - 1}}\left[ {P\left( {{k_x},{k_y},{z_H}} \right)G_D^{ - 1}\left( {{k_x},{k_y},{z_H} - {z_S}} \right)} \right]。$

(9)

式中，$G_D^{ - 1}\left( {{k_x},{k_y},{z_H} - {z_S}} \right)$为声压重建的逆向传递因子。

基于空间声场变换的平面近场声全息逆向重构流程，如图1所示。

2 滤波器

在实际声场反演过程中，位于高波数区的倏逝波成分将与噪声误差一起被逆传递函数按指数规律放大，若不对噪声信号进行低通滤波，将导致重建结果产生很大误差，严重影响声场重建的准确性，因此必须在声场逆向重构过程中进行滤波处理。

2.1 指数滤波器

指数滤波器是声场重构过程中常用的滤波方法，其算法简单、可操作性强、应用广泛。截止波数和陡度系数是该方法的2个关键系数，截止波数对参与声场重构过程的测量面声压角谱范围做出了限定，其取值受声波频率、信噪比以及测量面与源面间距离等因素的影响。截止波数取值过大，则达不到滤波效果，重建误差大；截止波数取值偏小，将导致过滤波，影响重建精度。陡度系数决定着窗函数的衰减速率，其取值越小，窗函数在截止波数处截止得越陡峭，且陡度系数会同时影响通带和高频带内滤波曲线的衰减速度，有时无法兼顾两者之间的平衡，甚至会导致有价值的倏逝波被剔除，影响滤波效果，进而影响声场重建准确性。

指数滤波器计算公式为：

$ W\left( {{k_x},{k_y}} \right) = \left\{ {\begin{split} &{1 - \frac{1}{2}\exp \left[ {\left( {{k_r}/{k_c} - 1} \right)/\alpha } \right],{k_r} \leqslant {k_c}} ，\\ & {\dfrac{1}{2}\exp \left[ {\left( {1 - {k_r}/{k_c}} \right)/\alpha } \right],{k_r} > {k_c}} 。\end{split}} \right. $

(10)

式中：${k_r} = \sqrt {\left( {k_x^2 + k_y^2} \right)} $；$ {k_c} $为截止波数；$ \alpha $为窗函数陡度系数。

2.2 巴特沃斯滤波器

巴特沃斯滤波器被称为最大平坦滤波器，其在通带的频率响应曲线符合最大限度平坦特性，是连续衰减的滤波器，该滤波器能够保持增益的平坦特性，尤为适合设计低通滤波器进行滤波处理，消除高频噪声干扰。巴特沃斯滤波器过渡带的陡峭程度与阶次成正比，随着阶次的增加，幅频特性曲线越陡峭，巴特沃斯滤波器更接近于理想低通状态，抑制噪声的能力变强。随着截止频率增加，系统响应速度变快，但系统稳定性也随之下降^[14]。

巴特沃斯滤波器计算公式为：

$ f({\text{C}},N) = \frac{1}{{1 + {{{C}}^2}{{\left( {\dfrac{{{k_r}}}{k}} \right)}^{2N}}}}。$

(11)

式中：$C = \left( {{{10}^{0.1{\alpha _p}}} - 1} \right)$为巴特沃斯滤波器截止频率；${\alpha _p}$为通带最大允许衰减系数；N为巴特沃斯滤波器阶数；${k_r} = \sqrt {\left( {k_x^2 + k_y^2} \right)} $。

3 仿真研究与参数分析 3.1 两脉动球源声场模型

如图2所示，两脉动球源的中心坐标分别为${C_1}{\text{ = }}\left( {0.175,0,0} \right)$和${C_2}{\text{ = }}\left( { - 0.175,0,0} \right)$，半径${r_0}$为0.01 m，两脉动球同相振动，表面振动速度${U_a}$=1 m/s，频率f=1000 Hz。重建面${z_S}$=0.01 m，测量面$ {z_H} $=0.05 m，重建距离为0.04 m，采样间距为0.05 m，测点数量为25×25，在测量面上等间距分布。

脉动球源辐射声压表达式：

$ {{p}} = \frac{A}{r}{e^{ikr}} 。$

(12)

式中：r为单个声源到测点距离；$\dfrac{A}{r}$绝对值为声压振幅。

$ {{A}} = - {{i}}\frac{{k\cdot \rho \cdot{c_0}}}{{4 {\text{π}} }}{\kern 1pt} {Q_0}。$

(13)

式中：$\rho $为空气密度；${c_0}$为空气中声速；${{Q}} = 4{\text{π}} r_0^2 \cdot {U_a}$。

在频率为1 000 Hz时，测量面测得声压分布图如图3所示。在理论计算的基础上，添加高斯白噪声，模拟信噪比为20 dB时，测量面测得声压如图4所示。

3.2 指数滤波器参数设计

为了能够定量描述滤波器滤波效果，定义重建误差为：

$ {\kern 1pt} E = \frac{{\sqrt {{{\left| {{p_s} - p_s^{the{\text{ }}}} \right|}^2}} }}{{\sqrt {{{\left| {p_s^{the{\text{ }}}} \right|}^2}} }} \times 100 \text% 。$

(14)

式中：${p_s}$为重建声压；$p_s^{the{\text{ }}}$为理论声压。

声源频率为1000 Hz时，指数滤波器误差曲线如图5所示。可知，随着截止波数的增加，重建面误差急剧下降，而在该阶段相较于截止波数的取值变化，陡度系数取值变化对重建误差影响较小。当指数滤波器截止波数到达${k_c} = 2.5k$附近时，滤波效果最好，重建误差最小，此时重建效果达到最优。当截止波数继续增大，参与重建过程的测量面角谱范围进一步扩大，被噪声误差淹没的高波数区倏逝波未被滤掉，滤波效果变差。当$ {k}_{c}=2.5k，\alpha \text{=}0.2 $时，重建误差达到最小，选取此时的指数滤波器参数为最优参数，该参数下指数滤波器如图6所示。

3.3 巴特沃斯滤波器参数设计

声源频率为1000 Hz时，巴特沃斯滤波器误差曲线如图7所示。可以看出，在阶次$N \leqslant 3$时，随着通带允许最大衰减系数${\alpha _p}$的取值增大，滤波器截止频率增大$\left( {C = \left( {{{10}^{0.1{\alpha _p}}} - 1} \right)} \right)$，重建误差不断降低，误差曲线存在最低点，此时重建误差最小，重建效果达到最优。随着阶次继续增加，误差曲线不存在拐点，随通带允许最大衰减系数${\alpha _p}$的取值单调递增，${\alpha _p}$越大，截止频率越大，重建过程中，过高波数的倏逝波未被滤掉，将连同噪声误差被逆传递函数按指数规律放大，造成极大的重建误差。当$N{\text{ = 2}}$，${\alpha _p}$=0.6时，重建误差达到最小，选取此时的巴特沃斯滤波器参数为最优参数，该参数下巴特沃斯滤波器如图8所示。

3.4 巴特沃斯滤波器与指数滤波器滤波效果分析

当声源频率为1000 Hz时，分析图5指数滤波器误差曲线可知，指数滤波器选择最优参数：截止波数${k_c} = 2.5k$，窗函数陡度系数$\alpha {\text{ = }}0.2$时，重建误差最小值为14.07%；分析图7巴特沃斯滤波器误差曲线可知，巴特沃斯滤波器在选择最优参数阶次$N{\text{ = 2}}$，通带最大允许衰减系数${\alpha _p}$=0.6时，重建误差最小值为13.31%。经指数窗函数滤波后重建结果如图9所示，经巴特沃斯窗函数滤波后重建结果如图10所示。通过与理论值对比发现，这2个滤波器均较好地实现声场重构，重构误差相对较小，其中巴特沃斯滤波器重建误差较指数滤波器更小。这一结果说明在声场模型相同时，选择最优参数的巴特沃斯滤波器滤波效果比选择最优参数的指数滤波器效果更好，巴特沃斯滤波器在声场重构过程中可进行有效滤波。

不同频率条件下，巴特沃斯滤波器与指数滤波器误差曲线如图11所示。重建距离为0.04 m时，巴特沃斯滤波器全频段平均误差为14.3%，指数滤波器为14.5%，全频段滤波效果相当。从图中可看出，在高频段，巴特沃斯滤波器重建误差明显小于指数滤波器重建误差，滤波效果更好。这是由于指数滤波器的陡度系数会同时影响通带内和高频带内滤波曲线的衰减率，且不能兼顾两者之间的平衡，导致有价值的倏逝波被滤出，从而影响滤波精度。后移测量面至${z_H}$=0.1 m，重建距离为0.09 m时，不同频率条件下，巴特沃斯滤波器与指数滤波器误差曲线如图12所示。此时，巴特沃斯滤波器全频段平均误差为30.6%，指数滤波器为31.6%，巴特沃斯滤波器滤波效果更好。从图中可看出，改变重建距离后，在高频段巴特沃斯滤波器重建误差仍明显小于指数滤波器重建误差，滤波效果更好。

4 结　语

针对指数滤波器无法兼顾通带衰减率和高频带的衰减率，可能导致有价值的倏逝波被剔除，影响声场重建精度问题，本文提出采用巴特沃斯数字滤波器作为滤波函数的近场声全息滤波方法，选取2个同向振动的脉动球源作为研究对象，对两脉动球源声场反演过程进行数值仿真。通过仿真分析得到了使声场重建误差最小的巴特沃斯数字滤波器最优参数，并使用最优参数下该滤波器与最优参数下指数滤波器的滤波效果进行了对比。结果表明，巴特沃斯滤波器可替代指数滤波器在声场重构过程中，进行滤波处理，且在高频段巴特沃斯滤波器重建误差明显小于指数滤波器重建误差，滤波效果更好，重建精度更高。