基于Green函数的两种STSF算法的比较

马佳男，杨德森，时胜国，胡 博

(哈尔滨工程大学 水声技术重点实验室，哈尔滨 150001)

现有的基于Neumann边界条件空间声场变换的有限离散化算法主要采用了k-空间抽样格林函数法对声场进行重构，该算法直接对波数域下的格林函数进行抽样，从而得到波数空间格林函数的离散阵列。该算法虽然简单，工程上易于实现，但是由于Neumann边界条件下k-空间振速-声压格林函数在辐射圆周上存在奇异性，这种奇异性会使得格林函数的幅值在辐射圆周上具有很大的跃变，从而影响重构精度。

1987年，Veronesi等人［1］首先提出实空间积分格林函数法，2004年，合肥大学的陈晓东［2］将该算法运用到了基于Dirichlet边界条件下的声压场重构，获得了很好的重构效果，并通过数值仿真证明该方法不仅适用于近场也同样适用于远场。随着矢量水听器的发展，可以通过测量质点振速来对声场进行重构和预测［3-4］。由于介质质点振速方向与波达方向一致，所以振速提供了声源方位的信息。在实际测量中更偏重于通过测量质点振速来重构声场以获得更好的重构精度和定位信息。

本文将实空间积分格林函数法运用到了Neumann边界条件下基于法向质点振速测量的声压场重构中，推导获得Neumann边界条件下该有限离散化算法的计算参数，通过对其傅里叶变换后的角谱求逆得到声压-法向质点振速声场重构时格林函数的角谱，并由数值仿真给出了该算法与k-空间抽样格林函数法在声场重构中基于不同重构参数的误差分析，从而比较两种格林函数法在声场重构中的优劣。

1 Neumann边界条件下的格林函数

理想流体介质中小振幅声波传播的波动方程［5］:

对于单频声波速度势Φ=φe-jωt，ω为声波的角频率，φ是空间分布函数，它满足Helmholtz方程:

其中，c为声速，k为声波波数，故波数空间又称为k-空间。通过Helmholtz方程，可以推导出Helmholtz公式。此公式用声场边界函数值表示声场(稳态单频波动声场)的积分形式解。当点源都集中在某一封闭曲面s内时，Helmholtz公式表示为:

n为封闭曲面的外法线方向，φ0为Helmholtz方程的解析解。ejkr/r为所选的辅助函数，其中，r表示场点o到封闭曲面s的距离，此函数除在r=0点有奇点外，其它地方皆满足假设条件和波动方程。

由式(3)可见，Helmholtz公式用φ和∂φ/∂n边界值的面积分来确定声场中任意一点的速度势函数值，因此当已知边界质点振速的分布和声压的分布值时，就可以用Helmholtz积分求出场中任意点的速度势函数值。

格林函数表示一定边界条件下点源的场，与边界条件一一对应。单频声场中的格林函数满足下面的方程［6］:

其中，r'代表声源的位置，r代表场点的位置，δ函数表示点源，时间因子取 e-jωt，解得:

式(5)表示的是自由场中的格林函数，称为格林函数的基本解。该式为具有1/r奇点形式的函数且满足Helmholtz方程，可以成为式(3)的辅助函数。利用格林函数这种可选性可以选择适当的格林函数形式来简化Helmholtz公式。

Helmholtz方程与第二类边界条件构成的定解问题叫做第二边值问题或Neumann问题。对于式(3)，第二类边界条件是指∂φ/∂n在区域边界上为给定函数。相应地，该边界条件下满足式(5)和Neumann边界条件的解称为Neumann格林函数。

根据“虚源法”，平面边界下格林函数的数学表达式为:

其中，r″表示虚源的位置。式(6)第二项对应着“虚源”，表示由平面边界引起的反射波部分。在平面边界绝对硬下，Helmholtz公式可以简化成一项。由“虚源”和实际声源的相位关系以及式(5)得到Neumann格林函数［7］:

gN即为Neumann格林函数。由欧拉公式:

其中uz(x，y，z)，p(x，y，z)分别为空间点(x，y，z)处为法向质点振速、声压。

通过式(8)和式(9)可以得到Neumann边界条件下，实空间域下法向质点振速-声压(简称为振速-声压)的格林函数，:

当∂φ/∂n=0时，对应为 Neumann边界条件，此时对式(3)进行二维空间Fourier变换并利用欧拉公式，整理得到k-空间振速-声压的格林函数:

以及k-空间声压-法向质点振速的格林函数:

其中:

2 实空间积分格林函数法

本文应用的实空间积分格林函数法是对已经离散化的实空间格林函数各子块的中点(x0，y0)处进行二元泰勒级数的展开，通过取不同阶数的泰勒展开来近似连续光滑实空间格林函数在子块上的积分值，即式(10)的积分值［8］。本文参考文献［3］，修正了文献［1］给出的计算参数，利用取泰勒展开式的前4阶得到实空间积分格林函数的近似结果。

首先，令全息孔径为Lx×Ly，将全息面离散为N×N矩阵，其中每个子块的大小为 Δx×Δy，其中 Δx=Lx/N，Δy=Ly/N，于是有实空间格林函数的离散表达式:

其中，m1，m2分别代表x，y方向上离散格林函数的子块序号，于是各子块中心坐标x0=m1Δx，y0=m2Δy，dz为重建面与全息面间的距离，这里简称为重构距离。

对式(13)的各子块中心做4阶泰勒展开，通过整理得到的数学表达式为:

上式给出的系数ai，是通过泰勒四阶展开得到的格林函数展开系数:

上式b=kx0，d=ky0。令a=kR0，

式(14)是非中心子块4阶泰勒展开的结果，当x0=0，y0=0且dz→0时R0→0，这时格林函数按照上式的计算方法具有奇异性，不能进行直接的近似，所以在计算中心子块积分值时将该子块分成两部分，即圆形区域和四个角区域。

于是，带入式(10)所得圆形区域的积分可用极坐标表示为:

然后，讨论正方形中剩余的4个角域上的积分，因为格林函数的对称性，所以这4个积分相同，只需求出其中一个即可，于是有:

将中心子块的格林函数与非中心子块格林函数相加，就得到了基于法向质点振速测量Neumann边界条件下的实空间积分格林函数。

令采样点数N=64，测量孔径为4 m×4 m，声速c=1 500 m/s，对比在不同重构距离下，k-空间抽样格林函数与实空间积分振速-声压格林函数的幅值分别在k空间和实空间域上的分布特点:

图1 两种算法下振速-声压格林函数幅值分布图Fig.1 The amplitude distribution of the two algorithms in the vibration velocity reconstructing the sound pressure

由图1可以看出:

(1)振速-声压实空间积分格林函数在频率和全息孔径不变的条件下，随着重构距离的减小，其幅值变得尖锐，且孔径中心与边缘的幅值衰减急剧增加;而当重构距离保持不变时，该格林函数幅值随着频率的增加而变得陡峭;以上分布特点说明当重构距离小于0.05m或频率增大时，由于实空间积分格林函数曲线的面积仅限于x=y=0附近很窄的范围内，从而导致该区域下采样点数较少，从而影响重构精度。

(2)k-空间抽样格林函数在辐射圆周上有跃变，随着重构距离的增大，辐射圆周上的跃变变得愈发尖锐;同时，声源频率的增加同样可以加剧辐射圆周上的跃变。这种跃变会使重构面孔径边缘处的函数具有奇异性，增大重构误差，所以由此可以推断该算法下的格林函数对重构距离和声源频率比较敏感，在对声场进行重构时会产生较大误差。

3 数值仿真分析

3.1 振速－声压的声场重构

下面利用实空间积分格林函数法和k-空间抽样格林函数法分别进行振速-声压的声场重构。取与图1相同的仿真条件，f=1 000 Hz，zH为全息面与声源间的距离。

图2 两种算法基于法向质点振速测量的声压场重构Fig.2 The two algorithms reconstructing the sound pressure field based on measuring normal particle vibration velocity

图2给出的是两种算法下，法向质点振速重构声压场的等高线图。从图中可以看出:

(1)基于实空间积分算法的声压场重构，虽边缘有起伏，但该算法下得到的声压幅值具有良好的衰减特性。而当dz和zH增大后，起伏也随之增大，由于NAH应用的条件之一是声压幅值在全息面边缘上具有足够大的衰减，因此误差大的这些点对整个声场的空间变换影响很小，而影响大的区域在孔径中心，只要将该区域的幅值误差控制在一定的范围内，就能最终保证重构精度。

(2)k-空间抽样格林函数法在孔径边缘不仅波动较大，且在dz、zH很小的情况下，“卷绕误差”使孔径边缘产生了虚源;随着dz和zH的增大，虚源强度增大从而影响对声源特性的识别。即使利用汉宁窗和k域滤波，都不能抑制卷绕误差对重构结果的影响。

从式(9)可以看到，当声源频率为声速的整数倍时，k-空间格林函数的分子为零，存在奇异性，不能通过k-空间抽样格林函数法对声场进行重构。所以我们重点考察声源频率对实积分格林函数法和声压场重构精度的影响。令dz=0.01 m，zH=0.1 m。

由图(3)可以看出，实积分格林函数在低频条件下重构声压幅值具有较好的收敛性和衰减性;在高频条件下，重构声压幅值衰减率减小，且孔径边缘有起伏，重构误差增大。

图3 不同声源频率下实积分格林函数法的重构声压场Fig.3 Using the real integral Green algorithm to reconstruct the sound field with different frequencies

3.2 声压－法向质点振速的声场重构

将经过FFT变换得到实空间积分法向质点振速-声压格林函数的角谱求逆，就可以得到Neumann边界条件下，声压-法向质点振速格林函数的角谱，通过该谱函数对法向质点振速进行声场重构。

下面利用实空间积分格林函数法和k-空间抽样格林函数法分别对基于声压测量的声场进行声压-振速的声场重构。取与图2相同的仿真条件，令dz=0.02 m，zH=0.15 m。

图4 两种算法基于声压测量的法向质点振速场重构Fig.4 The two algorithms reconstructing normal particle vibration velocity field based on measuring the sound pressure

从图4可以看出相对于声压的幅值分布，法向质点振速的幅值函数的收敛性更好，在孔径中心较为尖锐。所以为了更好地比较两种算法，图4给出的是两种算法的剖面图。在k域滤波参数相同的条件下，实空间积分格林函数法孔径边缘波动较小，而k-空间抽样格林函数法孔径边缘的波动较大。为了更直观地对比两种算法的优劣，利用下式计算重构声场的累积误差:

图5分别给出了两种算法下，不同的声源频率、重构距离、全息面距离重构声场的幅值误差变化曲线。其中:图5(a)中zH=0.15 m，dz=0.02 m;图5(b)中 dz保持0.02 m 不变;图 5(c)中zH=0.15 m。

可以看出:① 相对于k-空间抽样算法，实空间积分算法不仅可以得到较高的重构精度，而且各重构参数的累积误差随着参数的变化相对稳定，说明该算法具有较强的适用性;② 图5(c)虚线表示当dz＞0.04 m时，k-空间抽样算法对法向质点振速场重构失败，因为孔径边缘的幅值出现的指数型增长已经无法对声源做出精确的识别。虽然由实空间积分格林函数在dz＞0.05 m处得到的累积误差很大，但是其误差主要集中在孔径中心，孔径边缘处的误差很小，几乎不影响对声源特性的判断。

图5 不同重构参数对重构精度的影响Fig.5 Different reconstruction parameters on the reconstruction accuracy

4 结论

通过数值仿真结果可以得出如下结论:

(1)当通过测量法向质点振速对声压场进行重构时，相对于k-空间抽样法，实积分格林函数法有效地抑制了重构中常存在的“卷绕误差”，且该算法对声源频率不如k-空间抽样格林函数敏感，虽然在高频条件下，孔径边缘有波动，但没有出现较大的误差，对识别声源的影响较小。说明在基于空间声场变换的全息重构中，由实空间积分得到的格林函数的角谱更理想并更好地改善了法向质点振速-声压k空间格林函数在辐射圆周上的奇异性。

(2)振速-声压法向质点振速进行声场重构时，实空间积分格林函数法有效抑制了由于逆向重构产生的“吉布斯”效应，且各重构参数的累积误差随着参数的变化相对稳定，累积误差相对较小。

(3)通过对比和分析，Neumann边界条件下的实空间积分格林函数法相对于k-空间抽样格林函数法具有更好的重建精度，更有效的识别声源，而其应用特点也可为进一步的工程实践提供参考。

［1］ Veronesi W A，Maynard J D.Nearfield acoustic holography(NAH)II.Holographic reconstruction algorithms and computer implementation［J］.J.Acoust.Soc.Am.，1987，81(5):1307-1322.

［2］陈晓东.近场平面声全息的测量和重构误差研究［D］.合肥:合肥工业大学，2004.

［3］张永斌，毕传兴，陈 剑，等.基于质点振速测量的近场声全息技术［J］.农业机械学报，2007，38(9):112-116.

［4］ Jacobsen F，Liu Y.Near field acoustic holography with particle velocity transducers［J］.J.Acoust.Soc.Am.，2005，118(5):3139-3144.

［5］何祚镛，赵玉芳.声学理论基础［M］.北京:国防工业出版社，1981.

［6］梁昆淼.数学物理方法［M］.北京:高等教育出版社，1995.

［7］ Veronesi W A，Maynard J D.Digital holographic reconstruction of source with arbitrarily shaped surfaces［J］.J.Acoust..Soc.Am.，1989，85(2):588-598.

［8］ Wu S F.Methods for reconstructing acoustic quantities based on acoustic pressure measurements［J］.J.Acoust.Soc.Am.，2008，124(5):2680-2697.