0   引 言

自20世纪80年代以来，近场声全息^[1]在声源识别和声场可视化等领域发挥着重要作用^[2]，该技术通过近距离测量声源，高效捕捉声场倏逝波，从而提高声场重建性能。随着近场声全息的不断发展，基于空间傅里叶变换法^[3]、边界元法^[4]、统计最优法^[5]、等效源法(equivalent source method, ESM)^[6-7]等多种方法被相继提出。其中，ESM因适用于任意形状声源，鲁棒性好、计算简单等优点^[8]得到广泛应用。但因等效源数量远大于传声器数量，导致ESM的求解过程存在不适定性问题，无法求得精确解。由此，给出不同条件求近似解，成为声场重构的主流研究方向和热点问题。传统方法主要应用L₂范数正则化方法，如吉洪诺夫(Tikhonov)正则化^[9]和截断奇异值分解^[10]等，但该类方法受奈奎斯特(Nyquist)采样定理的限制，重建效果较差。

随着压缩感知(compressed sensing, CS)技术的出现，基于L₁范数最小化条件的稀疏重构方法被引入到ESM中^[11]，其思想是仅用很少的非零项来描述源强，实现信号的完美重构^[12]。2015年Fernandezg-Grande等率先将CS理论应用到基于球型传声器阵列的波展开中^[13]，并详细研究了所提出的压缩等效源法(L₁-ESM)^[14]，该方法在稀疏等效源法框架下，结合L₁范数惩罚约束，利用CVX工具箱进行求解，针对简单源，在较宽频率范围内获得了较好的重构精度。Hald提出了宽带声全息，利用最速下降法求取稀疏解^[15-16]，但这两种方法在低频时，对相干声源的重建精度不佳，这是因为稀疏过程未考虑声源本身特征，稀疏性仅约束解本身，当声源不是稀疏分布时，无法获得最佳重构性能。加权稀疏模型^[17]是考虑声源特性的正则化方法，利用声源的结构特征构建加权约束矩阵，将矩阵附加到L₁范数正则化中，实现促进解向量结构稀疏性的同时增强解向量本身稀疏性。现有的矩阵包括全变分^[18]、广义全变分^[19]等模型。1992年全变分模型由Rudin等^[20]提出后，依赖其利用空间梯度促进分段常数估计的特性，广泛应用于图像、声学^[21]等领域。2018年Fernandez-Grande将上述变分模型引入近场声全息声场重建，该模型利用一个近似二维的拉普拉斯(Laplace)矩阵构建结构矩阵，提高了解空间导数和解本身的稀疏性，扩展了CS在宽带重构中的适用性^[22]。然而约束矩阵的构建依赖于等效声源的结构特征，传统等效声源布置结构特征不明确，导致数据模型庞杂，计算效率低，重构精度达不到要求，所以如何合理布置等效源成了一个关键问题。Bi等^[23]证实了传递矩阵条件数对重构精度的影响。刘延善等^[24]在反复假设多层等效源面不同距离的基础上，保留各源面中较大源强，从而得到最优等效源面布置形式。王庆华^[25]、倪伟航^[26]针对相干声场，通过源强先验确定声源大体位置，再利用目标优化算法重新布置等效源点坐标，从而使等效源更好代替声源，提高重构精度。这些方法的本质都是通过减少传递矩阵条件数，确定优化等效源数量及位置。

本文提出一种声场全变分联合等效源法(total variational joint equivalentsource method, TV-ESM)的优化方法进行相干声场的稀疏重建。该方法以传声器到不同等效源点与最近等效源点之间的距离差作为筛选较大贡献等效源点的依据，并以较低传递矩阵条件数作为目标，确定优化等效源坐标。以全变分矩阵为基础，构建优化等效源结构特征所对应的单声源和相干声源结构约束矩阵，提高源强解的稀疏性。本文结合约束矩阵和L₁范数正则化方程，构建全变分稀疏重构声场模型，利用凸优化求解，实现声场稀疏重构，达到重构高精度相干声场的目的，拓展了稀疏等效源法的应用范围。

1   稀疏等效源法原理

等效源法(ESM)基本原理是任意声源产生的声场可以被多个等效声源所辐射的声场叠加替代，等效声源的源强可以通过全息面声压反求得到，以此实现声场重建。ESM的原理示意图如图1所示。

图  1  ESM原理示意图

Fig. 1  Schematic diagram of ESM principle

下载: 全尺寸图片 幻灯片

如图1所示，假设全息面有$M$个传声器测量点，等效源面有$N$个虚拟等效源点，则第$m$个传声器的测量压力可表示为

式中：${\boldsymbol{r}}_{m}$,${\boldsymbol{r}}_{n}$表示传声器和等效声源的位置矢量；$g\left({r}_{m}\left|{r}_{n}\right.\right)$表示第$n$个等效源与第$m$个测量声压之间的传递函数。在自由空间中，有：

式中：$k=\omega /c$为波数，$\omega$为角频率，$c$为空气中的声速；$\left|{\boldsymbol{r}}_{mn}\right|$为第$m$个测点与第$n$个等效源之间的矢量距离$\left|{r}_{mn}\right|=\left|{\boldsymbol{r}}_{m}-{\boldsymbol{r}}_{n}\right|$。测量声压可以用矩阵形式表示为

式中：$\boldsymbol{p}$为全息面测量声压构成的$M$个列向量；$\boldsymbol{q}$为等效源强构成的$N$个列向量；$\boldsymbol{G}$为等效源面与全息面之间的传递矩阵。

当等效源分布是稀疏或近稀疏时，通常对等效源强进行稀疏约束可获得准确声源信息。在理想情况下，可直接利用${{\mathrm{L}}}_{0}$范数求解，由于这是一个难解决的非凸优化问题，且在此条件下可与具有稀疏性的${{\mathrm{L}}}_{1}$范数等价，故通常采用${{\mathrm{L}}}_{1}$范数代替${{\mathrm{L}}}_{0}$范数进行优化求解：

式中：${\left\|·\right\|}_{1}$表示${L}_{1}$范数，表示解向量的绝对值之和，它的特点是通过减少解向量中非零值的数量提高解的稀疏性，但并不关心解向量的结构特征的稀疏性。最后通过式(5)重建得到声场任意点声压：

式中：${\boldsymbol{G}}_{\mathrm{r}}$表示等效源面与重建面之间的传递矩阵。

2   声场全变分联合等效源优化方法

由于式(4)不关心解向量的结构特征，当等效源分布非稀疏时，无法有效重建声源，因此本文以优化等效源的结构特征作为条件构造约束矩阵，该约束矩阵以一阶和二阶全变分(解的二阶导数)为基础，结合声场空间的稀疏性，直接约束信号结构，促进等效源稀疏解，拓展${{\mathrm{L}}}_{1}$范数对解向量结构特征的限制。其中问题在稀疏分析框架中制定，如：

式中：$\boldsymbol{D}$表示结构约束矩阵。

2.1   等效源优化

等效源求解是利用声场变换实现的，其对应关系由各平面之间的传递矩阵控制，而当传递矩阵条件数较大时，全息声压的细微扰动会使求解不稳定。这是由于等效声源离全息测点的距离增大时，倏逝波成分呈指数缩减，致使传递矩阵中最大奇异值减小，最小奇异值急剧下降，使传递矩阵条件数急剧增加，数据模型庞杂，计算效率低，最终导致求解难度的增加。因此提出将等效源点和测点之间距离引入到等效源布置当中，仅保留占主要贡献、距离传声器测点较近的等效声源，其余等效声源忽略不计。

当等效源模型确定时，到所有测点整体距离最短的等效源点，贡献量最大；而其余等效源点因距离相对较长，贡献量相对较小。经式(7)统计得出中心等效源点距测量面总体最短，所以设定其为最大贡献等效声源。

式中：${r}_{Mn}$代表第n个等效源到全部测点之间的距离。

为方便统计其余等效源点对声场贡献大小，以中心等效源点与其余等效源到阵列传声器距离差均值${r}_{k}$作为判断依据，并对其归一化处理，即：

式中：${r}_{MO}$代表中心等效源点到全部测点的距离；${r}_{k}$代表某等效源点与中心等效源点到各测点的整体距离差，距离差过大则代表此等效源点贡献量占比低，反之则表示此等效源点贡献量占比高。与倏逝波衰退过程一样，各等效源点贡献量随测量距离呈指数形式下降。为直观表现各等效源贡献量的差异，从上述距离与贡献量关系的分析中，推断距离差与贡献量之间满足指数分布形式：

式中：${P}_{{r}_{k}}$为${r}_{k}$所代表的等效源点贡献量的指数分布，选择几个相差较大距离差值作为备选；以式(10)传递矩阵条件数作为指标，选择最低条件数等效源作为优化等效源，从而确定等效源点的位置及数量。

式中：k_G则代表矩阵条件数，${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$是传递矩阵$\boldsymbol{G}$的最大奇异值，${\sigma }_{\mathrm{m}\mathrm{i}\mathrm{n}}$为最小奇异值。

2.2   全变分正则化

2.2.1   全变分数学定义

假设将等效声源按照列排序矩阵化为$\boldsymbol{q}\in {\mathbb{C}}^{{N}_{x}\times {N}_{y}}$，${\boldsymbol{q}}_{j}$表示矩阵中第$j$列元素，$j=1,2,\cdots {N}_{y}$，矩阵中任意列上相邻元素的变分定义为

式中：上述$\nabla {q}_{h}$为一阶水平全变分。相应的二阶全变分定义为

式中：$(\nabla {\boldsymbol{q}}_{h}{)}_{j}$、$({\nabla }^{2}{\boldsymbol{q}}_{h}{)}_{j}$分别表示矩阵$\boldsymbol{q}$的一阶、二阶$x$轴方向变分。

2.2.2   全变分矩阵表达

全变分矩阵由式(11)、式(12)得到一阶、二阶全变分矩阵${\boldsymbol{D}}_{1}$和${\boldsymbol{D}}_{2}$：

式中：${\boldsymbol{D}}_{1}$为$\left(N-1\right)\times N$维带状矩阵，它用任意相邻等效声源的差值作为惩罚约束，因此每行中${d}_{ii}$值为−1，${d}_{i(i+1)}$值为1；${\boldsymbol{D}}_{2}$为$\left(N-2\right)\times N$维带状矩阵，每行仍保留了${\boldsymbol{D}}_{1}$的特征。

2.3   正则化等效源优化模型

全变分约束矩阵依赖等效源结构特性及剩余数量和实际声源类型构建，将其设置为一对角阵，主对角线元素为不同阶次的全变分矩阵：

定义${\boldsymbol{D}}_{1d}$,${\boldsymbol{D}}_{2d}$为${\boldsymbol{D}}_{1}$,${\boldsymbol{D}}_{2}$的$\left(n-1\right)\times n$维子矩阵，$n$表示局部等效源的数量。如$n=4$时：

针对不同声源类型，选择构建不同结构约束矩阵。单声源辐射声场较为简单，本身具有较强的空间稀疏性，无需过多约束，故针对核心等效源部分选用单位矩阵$\boldsymbol{I}$进行适配；相较于简单源，相干声源辐射声场更为复杂，空间稀疏性较差，矩阵构建更为复杂，将优化后等效源分为3部分，此时等效源结构较之前发生明显变化，左右两侧等效源数量大幅减少，等效源从空间连续分布转化为块结构分布，并以此为约束构建结构约束矩阵，但相比于左半块和右半块，核心区域需保留更多有效声信息。即：

式中：${\boldsymbol{D}}_{\mathrm{t}\mathrm{v}1}$为单声源约束矩阵；${\boldsymbol{D}}_{\mathrm{t}\mathrm{v}2}$为相干声源约束矩阵。至此全变分约束矩阵${\boldsymbol{D}}_{\mathrm{t}\mathrm{v}}$构建完毕。

将全变分约束矩阵${\boldsymbol{D}}_{\mathrm{t}\mathrm{v}}$代入式(4)得到全变分正则化等效源优化方程：

3   仿真分析

为检验TV-ESM算法的声场重建性能，对该算法进行仿真分析。仿真时定义空气中的声速为340 m·s⁻¹，空气密度为1.29 kg·m⁻³，目标声源选用脉动球形声源，其半径设为0.01 m，振速为2.5×10⁻² m·s⁻¹。声源面、等效源面、重建面和全息面的位置按图1所示顺序布置，设定等效源面、重建面和全息面分别距离声源0.001、0.02和0.2 m，将声源中心设为坐标原点，等效源面和重建面网格尺寸均为0.5 m×0.5 m，等效源面上分布有21×21个等效源，网格间隔为0.025 m，全息面选用64传声器随机分布阵列，测点位置如图2所示。在仿真过程中添加30 dB高斯白噪声作为干扰噪声。

图  2  64传声器随机分布阵列

Fig. 2  Random distribution array of 64 microphones

下载: 全尺寸图片 幻灯片

3.1   优化等效源

图3给出了在距离差为0.1，0.4，0.7 的优化等效源和原始等效源的结构分布，图3中蓝色代表原始等效源布置，红色代表优化等效源布置。由图3可知，随着距离差变化，等效源面在−0.1<x<0.1、−0.1<y<0.1范围时，等效声源位置和数量变化较小，在其余范围则按上下左右的形式对称变化，变化幅度较大。这表示中心部分的等效声源占据的声场贡献量最大，保留数目最多，四周部分的等效源贡献相对较小，并向外以环形递减。图4给出了在频率为100~3 000 Hz时，不同布置等效源的传递矩阵条件数对比。从图4可看出，当距离差设置为0.4时，优化结果最好，因此选其作为最终优化等效源布置。

图  3  距离差不同时的优化等效源布置

Fig. 3  Arrangements of optimized equivalent source with different distance differences

下载: 全尺寸图片 幻灯片

图  4  不同等效源分布的条件数随率变化曲线

Fig. 4  Variation curves of the number of conditions with frequency for different equivalent source distributions

下载: 全尺寸图片 幻灯片

3.2   重建精度

对单声源和相干声源两种情况进行仿真性能分析。设置单声源位于坐标(0, 0, 0)处，相干声源位于坐标(−0.2 m, 0, 0)和(0.2 m, 0, 0)处。以传统L₁范数稀疏正则化方法（L₁-ESM）作为声场全变分联合等效源优化方法(TV-ESM)的对照对象。另外，为了直观评估方法优劣，定义重建声压幅值误差：

式中：${p}_{{\mathrm{t}}}$和${p}_{{\mathrm{r}}}$分别代表重建面上的理论声压和重建声压幅度，声压幅度误差E越小，重建效果越好。

3.2.1   单声源仿真

图5比较了重建面上L₁-ESM和TV-ESM方法对频率为200 Hz和2500 Hz单声源的重建声压云图。从图5可以看到，在频率为200 Hz和2500 Hz时，两种方法重建得到的声像面积非常接近理论值，重建效果都非常接近理论声压。结果表明两种方法都可以很好地估计单声源辐射声场。

图  5  L₁-ESM 和TV-ESM 对200和2500Hz单声源重建的声压云图

Fig. 5  Acoustic pressure nephograms reconstructed on the reconstruction plane by L₁-ESM and TV-ESM for the single sound source at 200 Hz and 2500Hz

下载: 全尺寸图片 幻灯片

图6为重建面上过单声源中心处所在行的声压幅度对比图，声源频率在200 Hz和2 500 Hz时，两种方法的重建声压幅度接近理论值。但是值得注意的是，L₁-ESM和TV-ESM方法在声源位置产生的峰值处理论值都未完全重合，在高频时，重建声压幅度更接近理论声压。

图  6  单声源声场的重建面上中间行的声压幅度

Fig. 6  Sound pressure amplitudes along the middle line on the reconstruction plane for the sound fields of the single sound source

下载: 全尺寸图片 幻灯片

图7为在100~3 000 Hz频率范围内两种方法的单声源声场重建误差曲线，每隔100 Hz取一个频率点，每个频率处的误差结果皆是计算10次以上求平均后得到。从图7可知，在较宽频率范围内(100~3 000 Hz)，TV-ESM和L₁-ESM两种算法具有相对稳定的性能，重建误差均保持在5%以下，重建精度较高。

图  7  采用L₁-ESM 和TV-ESM方法对单声源声场的重建误差比较

Fig. 7  Comparison of the reconstruction errors of L₁-ESM and TV-ESM for the single sound source

下载: 全尺寸图片 幻灯片

3.2.2   相干声源仿真

相干声源为两个或多个声波传播方向相同、频率相同，具有固定的相位差，在相遇区域内产生干涉现象的声源。图8是频率为200 、2 500 Hz时相干声源重建面声压云图。在200 Hz时，L₁-ESM方法虽识别出了两个声源，但真实声源位置出现了偏移，TV-ESM方法重建效果更好，声源位置被准确识别定位，声像面积与理论声压云图的误差不超过5%。在2 500 Hz时，两种方法都可以准确定位声源位置，且与理论结果十分接近。从图8可以看出，TV-ESM方法可以显著提高低频定位精度，扩大了适用频段的范围。

图  8  L₁-ESM 和TV-ESM方法对200和2500Hz的相干声源声场的重建面声压云图

Fig. 8  Acoustic pressure nephograms reconstructed on the reconstruction plane by L₁-ESM and TV-ESM for the coherent sound source at 200 Hz and 2500Hz

下载: 全尺寸图片 幻灯片

图9为相干声源声场重建面上中间行的声压幅度对比图。当声源频率设置为200 Hz时，L₁-ESM仅仅粗略识别到两个峰值，且两个峰值与理论值未能吻合。相比之下，TV-ESM所获得的声压幅值与理论值几乎完全一致，这也与图8相吻合。而当声源频率设置为2 500 Hz时，两种方法获得的声压幅度皆与理论值相近。

图  9  相干声源声场重建面上中间行的声压幅度

Fig. 9  Sound pressure amplitudes along the middle line on the reconstruction plane for the sound fields of the coherent sound source

下载: 全尺寸图片 幻灯片

图10为100~3 000 Hz频率范围内相干声源的重建误差曲线。由图10可知，在较宽范围内(500~3 000 Hz)，L₁-ESM和TV-ESM的重建误差均在10%以下，且中高频重建效果要优于低频。L₁-ESM在频率低于500 Hz时，重建精度很差，误差最高可达46%。TV-ESM方法在低频时的误差不超过15%，极大弥补了L₁-ESM方法在低频时识别失效的缺陷，具有更好的声场重建效果。

图  10  L₁-ESM 和TV-ESM方法对相干声源声场重构误差比较

Fig. 10  Comparison of the reconstruction errors of L₁-ESM and TV-ESM for the coherent sound source

下载: 全尺寸图片 幻灯片

3.3   重建性能影响因素

3.3.1   信噪比对重建精度的影响

本节主要分析信噪比R_SN对TV-ESM方法重建精度的影响。测量距离设置为0.2 m，TV-ESM方法在不同信噪比时，100~3 000 Hz频率范围内的重建精度如图11所示，各频率处的误差是通过5次计算求平均得到。从图11可知，随着信噪比的增大，TV-ESM方法的重建误差总体呈下降趋势。当信噪比大于30 dB时，TV-ESM方法的重建精度较高且稳定(在较宽的频率范围内，重建误差基本都在12%以下)。因此，使用TV-ESM方法进行稀疏声源声场重构时，为确保重建精度，建议信噪比设为30 dB及以上。

图  11  不同信噪比下重建误差曲线

Fig. 11  Reconstruction error curve under different signal to noise ratios

下载: 全尺寸图片 幻灯片

3.3.2   测量距离对重建精度的影响

图12为TV-ESM方法重建误差随测量距离的变化曲线。仿真中信噪比设置为30 dB，即干扰噪声较少。对比不同频率f时重建精度随测量距离的变化曲线，发现低频时各个重建距离下的重建精度低于高频时的重建精度，随着距离测量的增加，TV-ESM方法在全频段上的重建精度有所下降。

图  12  不同测量距离下重建误差曲线

Fig. 12  Reconstruction error curves at different measurement distances

下载: 全尺寸图片 幻灯片

4   实验验证

为进一步检验TV-ESM方法的声场重建效果，在半消声室进行实验验证。实验采用LMS Test.Lab声学测量系统，声源为扬声器，传声器型号为G.R.A.S 1/4 in 40PL。实验布置如图13所示。与仿真条件一致，设置单声源坐标为(0, 0, 0)，相干声源坐标为(−0.2 m, 0, 0)、(0.2 m, 0, 0)，声源距全息面0.2 m，声源距离地面高度与阵列中心距离地面高度相同均为1.2 m。设置全息面有10×10个采样点，随机选择36个采样点的声压作为输入声压，并将信号采集设备的固有噪声作为干扰噪声。

图  13  实验测试现场

Fig. 13  Pictures of experimental test site

下载: 全尺寸图片 幻灯片

图14为单声源试验重构声压云图。当声源频率为200 Hz和2 500 Hz时，两种方法均可以识别出单个扬声器位置。但是L₁-ESM方法的“鬼影”（虚假声源数）要多于TV-ESM方法，且随着频率的增加，“鬼影”数也在增加。主要原因在于高频分量波长短，声能量衰减快，同时阵列构型的需求更大，加之受到周围环境噪声的影响，“鬼影”是无法避免的，但是全变分约束矩阵能在一定程度上抑制噪声分量，因此TV-ESM的方法重构效果更好。

图  14  L₁-ESM和TV-ESM对200 Hz和2 500 Hz的单声源声场试验重建面的声压云图

Fig. 14  Experimental acoustic pressure nephograms reconstructed on the reconstruction plane by L₁-ESM and TV-ESM for the single sound source at 200 Hz and 2 500 Hz

下载: 全尺寸图片 幻灯片

图15为相干声源试验重构声压云图。当声源频率为200 Hz时，TV-ESM方法能准确识别出两个声源，且识别位置接近真实声源。L₁-ESM在数量和位置上重建效果均低于TV-ESM方法，这与仿真结果相对应。当声源频率设置为2 500 Hz时，两种方法的重建效果相近，均能识别出两个声源。

图  15  L₁-ESM and TV-ESM对200 Hz和2 500 Hz的相干声源声场试验重建面声压云图

Fig. 15  Experimental acoustic pressure nephograms reconstructed on the reconstruction plane by L₁-ESM and TV-ESM for the coherent sound source at 200 Hz and 2 500 Hz

下载: 全尺寸图片 幻灯片

5   结 论

本文提出了一种正则化等效源优化方法TV-ESM用于声场稀疏重构。该方法利用全息面测点到等效源点与测点到中心等效源点间的距离差作为等效源点对声场贡献的评定依据，筛选较大贡献等效源，以传递矩阵条件数为优化目标确定等效源布置形式；在稀疏等效源框架下，结合声源自身特征及等效声源结构特性构造全变分约束矩阵，通过凸优化求解源强。该方法通过优化等效源布置，明确等效声源结构特征，简化数据模型，提高了计算精度；并针对相干声源构建约束矩阵，解决传统稀疏等效源稀疏性不足的问题。数值仿真结果表明，与已有的L₁-ESM方法相比，TV-ESM方法适用于重构单声源和相干声源，扩展了稀疏等效源法在宽带重构中的适用性。