0 引　言

近场声全息技术是一种先进的噪声源识别方法，最早由Williams提出，并根据空间声场变换进行了大量实验应用^[1,2],其主要算法是通过在辐射体的近场区域内测量声场数据重构和预测整个三维空间的声学特性，从而识别和定位声源，获得声源的声辐射特性，达到判断噪声来源进行声源控制的目的。

等效源法兴起于20世纪90年代，其主要思想是：物体辐射的声场可以由置于该辐射体内部的一系列等效源产生的声场叠加代替，等效源声强则可以由全息面测量得到的声场数据进行反向推算得到^[3]。等效源法的主流方法有波叠加法^[4－5]、等效源模拟法^[6]等。根据等效源的声源形式区别可以将等效源法近场声全息分为采用简单源作为等效声源进行近场声全息的声场重构技术^[7]和用球面波作为等效声源进行近场声全息的声场重构技术^[7]。采用简单源作为等效声源进行近场声全息的声场重构技术主要通过全息面测量得到的声场声波推算出辐射源内的等效源，由加权矩阵进行波叠加模拟出外辐射声场，从而达到声场预测或重构的目的。近年较为代表的是陶文俊提出的用压缩感知将声强信息转为稀疏信号的方法^[8]，陈汉涛运用经验模态分解法与等效源法结合，使低信噪比条件下的弱声源得以检测的方法^[9]。而当用球面波作为等效源时，最具代表性的是基于Helmholtz方程的最小二乘法（Helmholtz equation least squares method，HELS）声全息技术^[10]。

得益于等效源法声全息在环境适用性和数值计算优势，该算法被大量应用于不规则形状的噪声检测。考虑到声场重构的计算效率与精度，在基于等效源的近场声全息算法的基础上提出一种基于虚拟声源定位的改进算法，分别在声源频率与重构距离变化情况下比较改进前后算法的重构精度与重构效果，利用仿真分析与实验验证，证明了基于等效源的近场声全息算法在水下运用的可靠性，以及基于虚拟声源定位进行声场重构时在声源分辨与声场还原的优势。

1 基于定位算法的等效源近场声全息

等效源法又称为波叠加法，该种算法通过辐射体内部等效源叠加替代该物体所产生的辐射场，因此适用于任何闭合形状的声源声场预测或重构，且因为其算法中不包含插值运算、特征波数的非唯一解处理与奇异值积分计算，因此具有使用环境适应性强、误差较小等优势。

在声向外辐射问题求解中，设S为辐射源的闭合外表面，对Helmholtz方程式进行求解，辐射声场中任意点r处的声压可以表示如下：

$ p(r) = \int_\Omega {i\rho ckq({r_0})} g(r,{r_0}){\rm{d}}\Omega ({r_0})。$

(1)

式中： $ \Omega $ 为在发声内部区域等效的连续分布声源； $ q({r_0}) $ 为半径在r₀处等效的声场强度； $ g(r,{r_0}) $ 为自由场中声场转换的格林函数。

通过式（1）可以看出，声辐射场内空间或辐射表面的任意一点声压可以用声辐射体内若干个等效声源体向外辐射的声场进行替代。而在实际计算过程中，连续分布的声源无法在计算机内实现，因此采用离散算法在辐射体内分布若干个等效声源来模拟辐射场的近似声场。因此，基于等效源法的计算方式，辐射场中点r处的声压可以表示如下：

$ p(r) = \sum\limits_{j = 1}^N {i\rho ckg(r,{r_{0j}})q({r_{0j}})}。$

(2)

式中：N表示等效声源的总个数； $ q({r_{0j}}) $ 表示第j个等效源的声源强度； $ {r_{0j}} $ 表示第j个等效源的柱坐标。

同理，辐射体表面r_s处的声波声压可以表示如下：

$ p({r_s}) = \sum\limits_{j = 1}^N {i\rho ckg({r_s},{r_{0j}})q({r_{0j}})}，$

(3)

假设计算声场的辐射体表面有M个边界点，且M≥N，那么可以得到M个与式（3）相同的等式，将整理成矩阵形式可以表示为：

$ {{\boldsymbol{P}}_s} = i\rho ck{{\boldsymbol{G}}_{sp}}{\boldsymbol{Q}} ，$

(4)

式中，声源表面的声压向量可以表示为 ${{\boldsymbol{P}}_s} = [p({r_{s1}}), p({r_{s2}}), \cdots ,p({r_{sM}})]^{\rm{T}}$ ，与声场对应的等效源强度向量可以表示为 ${\boldsymbol{Q}} = {[q({r_{01}}),q({r_{02}}), \cdots ,q({r_{0N}})]^{\rm{T}}}$ 。G_sp表示基于等效源法计算的近场声全息算法中等效源与辐射场中表面声压的传递矩阵，可以由下式进行表示：

$ {{\boldsymbol{G}}_{sp}}(m,j) = g({r_{sm}},{r_{0j}})。$

(5)

在声辐射研究中，若已知声源表面的声压强度时，可以通过式（5）反向计算出等效声源体的声源强度向量，通过式（4）计算得到辐射场中的声场声压，从而进行声场中声强与声功率的计算。

当声场推算过程中已知全息面的声压，假设该全息测量面有M个测量点，那么可以建立M个与式（5）相同的等式，通过矩阵形式进行表示时如下式：

$ {{\boldsymbol{P}}_h} = i\rho ck{{\boldsymbol{G}}_{hp}}{\boldsymbol{Q}}，$

(6)

式中，全息面上测量所得的声压向量 ${{\boldsymbol{P}}_h} = [p({r_{h1}}), p({r_{h2}}), \cdots ,p({r_{hM}})]^{\rm{T}}$ ，等效声源与全息面声压的传递矩阵 ${{\boldsymbol{G}}_{hp}}$ 表示为：

$ {{\boldsymbol{G}}_{hp}}(m,j) = g({r_{hm}},{r_{0j}})。$

(7)

因为随着辐射声场变化，每个声场的等效源会发生改变，若要求得等效源的声强向量，可以利用式（6）对传递矩阵 ${{\boldsymbol{G}}_{hp}}$ 求逆。

由于直接定义复杂声源发出的声场较为困难，因此利用若干个点声源模拟实际复杂声源的声场情况，复杂声源的外辐射声场中某点处声压可以表示为：

$ \begin{split} & p({r_j},{\theta _j},{\gamma _j}) = \sum\limits_{k = 1}^n {{p_k}({r_j} - {r_k},{\theta _j} - {\theta _k},{\gamma _j} - {\gamma _k})} , \\ &j = 1,2,3, \cdots ,m 。\\ \end{split} $

(8)

式中， $ p({r_j},{\theta _j},{\gamma _j}) $ 为声场中点 $ ({r_j},{\theta _j},{\gamma _j}) $ 处声压。点 $ ({r_k},{\theta _k},{\gamma _k}) $ 为第k个虚拟点声源坐标。

通过遍历法对虚拟声源位置进行求解，利用矢量水听器检测得到的声场信号与模拟声场做差，当声压误差最小时，此时虚拟声源发出的声场近似于实际工况中的声场情况，虚拟声源对应的位置近似于实际工况中声源的位置。此时最小误差定义为：

$ {e_p} = \min \left\{ {P - \sum\limits_{j = 1}^m {p({r_j},{\theta _j},{\gamma _j})} } \right\}。$

(9)

式中，P表示矢量水听器采集到的声场声压。

由于矢量水听器的信号采集具有指向性，且采集信号分为x，y，z，p四通道，将其转换为球坐标下的参数进行遍历，即可得到最接近工况中声场的模拟声源坐标位置。

经矢量信号反推后声场中有m个模拟声源（m<n），且每个虚拟声源有计算得到的位置信息 $ ({r_m},{\theta _m},{\gamma _m}) $ ，此时传递矩阵 ${\boldsymbol{G}}_{hp}'$ 可以表示为：

$ {\boldsymbol{G}}'_{hp} = {{\boldsymbol{U}}_m}diag({\sigma _1},{\sigma _2}, \cdots ,{\sigma _m}){{\boldsymbol{V}}^H} ，$

(10)

通过式（6）可以计算出等效源声强向量表达式：

$ {\boldsymbol{Q}} = \frac{1}{{i\rho ck}}{\boldsymbol{G}}'^{ + }_{hp}{{\boldsymbol{P}}_h} = \frac{1}{{i\rho ck}}{({\boldsymbol{G}}'^{H}_{hp}{\boldsymbol{G}}'_{hp})^{ - 1}}{\boldsymbol{G}}'^{H}_{hp}{{\boldsymbol{P}}_h} ，$

(11)

将式（10）代入式（11），整理可得基于等效源法的近场声全息声场重构公式为：

$ {\boldsymbol{Q}} = \frac{1}{{i\rho ck}}{\boldsymbol{V}}{\rm{diag}}{({\sigma _1},{\sigma _2}, \cdots ,{\sigma _m})^{ - 1}}{\boldsymbol{U}}_m^H{{\boldsymbol{P}}_h} 。$

(12)

由于等效源法近场声全息中等效源声强Q的求解是求逆问题，若直接计算会因为测量误差与环境噪声干扰等因素导致广义逆求解过程中误差被放大，从而给计算结果带来巨大误差，造成声场重构失败的结果。因此在声全息计算中引入Tikhonov正则化求解来抑制噪声干扰，从而使声场重构的结果更加稳定精确，即通过增加一个约束项，使得求解最小二乘解的病态方程变成非病态方程，如下式：

$ {J_\alpha }({{\boldsymbol{P}}_S}) = \min \left\{ {\left\| {i\rho ck{\boldsymbol{G}}'_{hp}{\boldsymbol{Q}} - {\boldsymbol{P}}} \right\|_2^2 + {\theta ^2}\left\| {{\boldsymbol{L}}{{\boldsymbol{P}}_S}} \right\|_2^2} \right\}。$

(13)

式中： $ \|\cdot\|_2 $ 表示矩阵的二次范数；L表示作为约束的惩罚矩阵，通常为正定矩阵； $ \theta $ 为正则化参数。

对式（13）的泛函数问题进行求解即可得等效源声强的正则化解为：

$ \begin{split} {{\boldsymbol{Q}}_{reg}} = &\frac{1}{{i\rho ck}}{({\boldsymbol{G}}'^{H}_{hp}{\boldsymbol{G}}'_{hp} + {\theta ^2})^{ - 1}}{\boldsymbol{G}}'^{H}_{hp}{{\boldsymbol{P}}_h} =\\ &\frac{1}{{i\rho ck}}\sum\limits_{i = 1}^m {\frac{{\sigma _i^2}}{{\sigma _i^2 + {\theta ^2}}}\frac{{u_i^H{{\boldsymbol{P}}_h}}}{{{\sigma _i}}}{v_i}} 。\end{split} $

(14)

通过上式计算得到等效源声强后，代入式（4）进行声辐射域内任意位置的声压重构，从而达到整个声场重构或预测的目的。

利用矢量水听器测量得到的声场信息对声源位置进行定位，能对等效声源的点数进行缩减，降低求解计算复杂度，并且减少等效源算法中虚拟声源的产生，能够降低算法的声场计算误差。下文将对比该方法与传统等效源算法的重构精度，分析比较2种重构方式在不同工况情况下的精度变化趋势。

2 基于定位算法的等效源仿真

为验证算法的可行性，建立圆柱壳体有限元仿真模型，壳材料设置为钢，仿真环境为声速为1500 m/s，密度为1000 kg/m³的水域，水域外设置完美匹配层以避免声波的反射。模型如图1所示，柱体内壁（0.2 m，0°，0.15 m）与（0.2 m，0°，0.35 m）处设置2个同频同相向外振动的激振点，将柱坐标中心设为2个激振点所在柱体，柱体半径0.2 m，有效振动柱面高0.5 m，重构柱面半径r_s=0.25 m，全息柱面半径r_H1=0.30 m，r_H2=0.35 m，r_H3=0.40 m，信号采集范围与柱面共形且与有效振动柱面等高，轴向采样点数与周向采样点数分别为21个、16个，研究方式采用频域研究。

将仿真的重构误差定义为分贝误差：

$ {e_p} = 20\lg \frac{{\sqrt {\displaystyle\sum\limits_{m = 1}^M {\sum\limits_{n = 1}^N {{{\left| {{p_{s(m,n)}}} \right|}^2}} } } }}{{\sqrt {\displaystyle\sum\limits_{m = 1}^M {\displaystyle\sum\limits_{n = 1}^N {{{\left| {{p_{h(m,n)}}} \right|}^2}} } } }} 。$

(15)

式中： $ {p_{s(m,n)}} $ 表示重构面的声场声压； $ {p_{h(m,n)}} $ 表示全息面的声场声压。

2.1 声源频率对2种算法重构精度的影响

在声源频率为100～1 000 Hz时对2种算法的重构精度进行计算，重构误差随声源频率变化趋势如图2所示，“+”表示改进前的等效源算法，“*”表示改进后等效源算法。

根据图2的变化趋势可以看出，随着声源频率增加，2种算法的重构误差也在增大。但定位算法对声场中的各向同性噪声屏蔽效果较好，能明显提高等效源算法中等效源场的估计，因此基于虚拟声源定位的改进等效源算法有着更高的重构精度。

2.2 重构距离对2种算法重构精度的影响

将声源频率设置为500 Hz，改变仿真中声场全息面的采集距离，重构柱面半径r_s=0.25 m，全息柱面半径为r_H1=0.30 m，r_H2=0.35 m，r_H3=0.40 m。对2种算法的重构精度进行计算，重构误差随重构距离变化趋势如图3所示，“+”表示改进前等效源算法，“*”表示改进后等效源算法。

根据图3变化趋势可以看出，随着重构距离增加，改进前后算法的重构误差均逐渐增大，但定位算法的抗干扰性优势更为明显，误差上升趋势与改进前相比更缓慢，对传统等效源算法的改善效果显著，提高了声场的重构精度。

3 实验与结果

为验证本文算法的可行性以及仿真的准确性，设计加工了内壁加肋圆柱体并在消声水池中进行近场声全息实验。柱体尺寸为 $ \varnothing$ 0.2 m×0.5 m×0.003 m，肋间距设计为0.1m用以稳定框架。桶底焊有激振器平台，通过装配槽固定激振器以保证实验过程中激振器与柱体保持相对静止。激振器通过桶盖上的水密接头开口与信号发生器，功率放大器相连，利用2根紧贴筒壁的激励杆带动柱体向外辐射声波，实现实时信号的发生与调节。将水听器布放在距离桶面 $ d=r_H-r_s$ 处，x通道正对声源方向，通过行走机构调整水听器在水中的位置，由上而下进行列向声场采集，采集得到的数据通过PC进行保存与后处理。当采集完一列数据后，控制固定端的旋转杆，使柱体旋转一定角度进行另一列声场的采集。重复上述步骤，即可完成21×16的柱面声场采集。实验系统示意图如图4所示。实验装置包括信号源、功率放大器、激振器、水听器、信号调理器、采集卡、上位机。

利用本文算法对采集到的声场信息进行重构，并分别计算不同工况下的重构精度进行比较。

3.1 不同声源频率时的重构情况

将激振器发出的声源频率分别设置为100 Hz，500 Hz，1 000 Hz，并采集距柱体中心r_s=0.25 m与r_H1=0.30 m，即重构距离为0.05 m的声场信号，采集密度为轴向21个点，周向16个点。对比改进前后算法的声场重构效果，对比结果如图5所示，重构误差如表1所示。

由图5可以看出，总体上改进前后算法在设定的频域范围内均能对声场信息进行一定程度的还原，但随着声源频率增大，声场中的旁瓣与虚像干扰逐渐严重，声源位置信息变弱，声场重构效果逐渐变差。具体对比改进前后声场重构效果，可以看出基于定位算法测量得到的声场声源位置信息更明显，旁瓣影响更小。通过误差计算可得，改进后算法误差在低频时降低1.97～5.44 dB，在声场重构过程中更有优势。

3.2 不同重构距离时的重构情况

将激振器发出的声源频率设置为500 Hz，分别采集距柱体中心r_s=0.25 m与r_H1=0.30 m，r_H2=0.35 m，r_H3=0.40 m时的声场信号，此时重构距离分别为0.05 m，0.10 m，0.15 m，采集密度为轴向21个点，周向16个点。对比改进前后算法的声场重构效果，对比结果如图6所示, 重构误差如表2所示。

由图6可以看出，总体上2种算法在距离较远处采集到的声场旁瓣会削弱，但重构效果都变差。具体对比改进前后算法，当重构距离逐渐增大时，改进前算法计算得到的声源位置已发生偏移，声源定位效果较差；基于定位算法计算得到的声场，由于算法计算对干扰信号的抑制作用，虽然声源分辨能力降低，但位置并未发生偏移，在近距离有更优秀的声源位置信息重现。通过误差计算可得，改进后算法的误差在重构距离变化时降低1.97～3.99 dB，在声场重构过程中更有优势。

4 结　语

本文在基于等效源近场声全息算法的基础上提出了一种基于虚拟声源定位的改进算法，分别在声源频率与重构距离变化情况下比较改进前后算法的重构精度与重构效果。仿真与实验结果表明：采用定位算法进行重构时，能够保留较高的声源位置信息，且计算时间较短，重构精度与效果有较为明显的提高。但由于算法中需要对声源位置进行预判，因此适用于近距离下，声源位置误差较小的非接触式噪声分析。