基于压缩感知和快速波叠加谱方法的近场声全息

张 阳，向 宇，石梓玉，陆 静，王玉江

(1. 广西科技大学广西汽车零部件与整车技术重点实验室，广西柳州545006；2. 广西科技大学机械与交通工程学院，广西柳州545006)

0 引 言

近场声全息技术(Near-Field Acoustic Holography, NAH)是一种噪声源识别、定位及声场可视化的有效技术[1]。经过几十年的发展，各国学者相继提出了空间傅里叶变换算法(Spatial Fourier Transform, SFT)[2]、边界元法(Boundary Element Method,BEM)[3]、等效源法(Equivalent Source Method,ESM)[4]等。其中，基于等效源法的近场声全息因原理简单、适应性强等优点得到了广泛应用[5-7]。为了保证该方法的重建精度，必须要有足够数量的测量点和源强点。但在NAH技术中，通常情况下是通过测量外域声场的信息来重建振动体的近场声学量，受限于测量成本和测量条件，一般很难获得测量点数大于源强点数的超定方程组。此时若采用基于l2范数的最小二乘方法求解欠定方程组，难以得到理想的结果。因此，如何在减少测点的情况下，仍能稳定获得较高的重建精度，一直是NAH技术重点关注的问题。

Patch近场声全息方法[8-9]是解决上述问题的一种有效途径，而压缩感知技术(Compressive Sensing,CS)[10-11]则提供了另外一种手段。该技术利用信号的稀疏性来解决欠定问题，若信号本身稀疏，或在某个变换域稀疏，即可通过一个与变换基无关的观测矩阵，在低于奈奎斯特采样定理的采样频率下高概率地重构出原始信号。其实质上是一种寻找欠定方程组稀疏解的技术，由于利用了基于l0或l1范数的重构算法替代传统基于l2范数的最小二乘法，因此在减少测点的情况下也能保证较高的重建精度。2012年，Chardon等[12]首次将压缩感知应用到NAH中，并通过实验证明了该方法在减少测点数时仍具有较高的计算精度。Fernandez-Grande等[13]将压缩感知理论与等效源法相结合，提出了一种压缩等效源法(Compressive Equivalent Source Method,C-ESM)。Hald等[14]在压缩等效源法的基础上，结合最速下降法提高了计算效率，并对其在高频段下的声场重建效果进行了研究。毕传兴等[15]提出了一种基于功率阻矩阵特征分解的压缩模态等效源法，并通过数值仿真与实验证明了该方法在减少测点数方面优于压缩等效源法。胡定玉等[16]基于等效源法与奇异值分解法，提出了一种用于空间扩展声源的稀疏采样声场重建方法。然而，上述方法都是将虚拟源强进行单元数值离散，这将对最终的求解精度造成一定影响，尤其是在高频声场的重建时，其对精度的影响更加严重。

作者在前期研究中，将等效源法中的虚拟源强和积分核函数在轴对称虚拟面上沿周向和子午线方向利用傅里叶级数进行展开，提出了一种求解轴旋转空穴三维诺伊曼(Neumann)声辐射问题的快速波叠加谱方法[17-18]。由于该方法是一种傅里叶级数展开形式的半解析、半数值方法，因而其计算精度高于传统的离散源强等效源法，同时作者在进一步的研究中发现，待求的傅里叶展开系数向量总是一个仅含少量非零值的稀疏向量，该特性正好满足了压缩感知理论对稀疏性的要求。因此，本文将文献[18-19]的快速波叠加谱方法与压缩感知技术进行有机结合，推广用于NAH技术中，提出了一种基于压缩感知和快速波叠加谱方法的近场声全息方法。文中首先推导等效源强和积分核函数在虚拟面上的傅里叶级数展开形式，并分析待求展开系数向量的稀疏性，然后在此基础上结合压缩感知算法求解，最后通过数值仿真验证所提方法在声场重建中的有效性与优越性。

1 基于压缩感知的等效源法

如图1所示，考虑一任意形状的三维振动体，其中S为振动体边界，E+为振动体外部辐射域，振动体内部布置一虚拟曲面SE，其内部区域记为E。

其中，ε表示一个和噪声相关的误差项。

2 结合压缩感知的快速波叠加谱方法

文献[17-18]中，将源强和积分核函数在轴对称的虚拟面上采用傅里叶级数进行展开，提出了一种求解轴旋转空穴三维 Neumann声辐射问题的快速波叠加谱方法，并取得了较好的声场计算结果。本文借鉴文献[17-18]的方法与思路并将其应用于近场声全息的声场计算中。

假设等效源点rE和场点r在柱坐标系下的位置坐标分别为，如图 2(a)所示。它们在各自子午面内的二维坐标分别为，如图2(b)所示。

图2 场点、源点位置示意图Fig.2 Schematic diagram of field point and source point

声波一般是由结构振动激发声介质产生的，根据结构动力学理论，结构在任意载荷激励下产生的动力响应主要由低阶模态叠加而成，因而结构振动辐射的声波也主要集中在低阶模态。此外，由文献[20]可知，声音在传播的过程中，低阶模态衰减慢、辐射能力强，高阶模态衰减快、辐射能力弱。换言之，随着传播距离的增大，声源对声场的贡献主要集中在低阶模态，高阶模态类似于基于空间傅里叶变换的近场声全息中的倏逝波模态，阶数越高衰减越快，也就是高阶的Cmn越来越小。因此，只要式(21)的求和截断项数mf和nf取得足够大，由该式表示的声场即包含了辐射能力强的低阶模态和足够多的倏逝波模态，因而傅里叶展开系数向量Cmn必然是一个具有一定稀疏性的稀疏向量，且截断项数mf和nf取得越大，其稀疏性越强。

由于Cmn是一个稀疏向量，因此可直接利用压缩感知技术进行求解，即在最小l1范数下求解如下问题：

3 数值算例

3.1 脉动球源声场重建验证

在添加噪声的情况下，利用脉动球源对比不同采样点数时本文方法与压缩等效源法的声场重建效果。其中，压缩等效源法与本文方法皆采用SPGL1工具箱[21]中的基追踪降噪算法(Basis Pursuit De-Noising, BPDN)进行求解。

已知脉动球源半径为rS，周边产生均匀径向振速幅值为v0的脉动球源在距离球心 r处所辐射声压的解析解[22]

式中：ρ为介质密度，c为介质中的声速。

仿真中设置脉动球的半径为 1 m，球心在坐标原点处。等效源点、全息采样点分别布置在半径为0.3、1.1 m且与声源面同心的球面上。在压缩等效源法中设置虚拟等效源离散间隔与全息面一致；在本文方法中取虚拟面周向、子午方向抽样数为N=M= 28，周向、子午线方向求和截断项数分别为nf=30、 mf= 5 0。全息面、等效源面结点划分方式见图 3。在进行计算时，两种重建方法均添加信噪比为20 dB的高斯白噪声。重建脉动球表面的114个离散结点声压，并计算其相对误差。相对误差 e定义为

图3 结点划分图Fig.3 Node partition diagram

图4显示了频率为1 000 Hz、采样点数范围为14～314时本文方法和压缩等效源法的重建误差曲线曲线上的数值为对应的采样点数。可以看出，当采样点数增加到 62时，压缩等效源法的重建误差趋向于 2%。但本文方法的误差曲线始终低于压缩等效源法，即使是在 14个采样点的情况下，其重建误差仍一直低于 2%。同时，本文方法的重建精度也不依赖于采样点数，采样点数的增加并未明显改善本文方法的重建精度，这是因为对于脉动球源，利用本文方法得到的傅里叶展开系数具有良好的稀疏性，因此 14个采样点就能够保证本文方法对于脉动球表面声压的重建精确。

图4 f=1 000 Hz时，不同采样点数下本文方法和压缩等效源法的重建误差曲线Fig.4 Reconstruction error curves of the proposed method and C-ESM under different sampling points at f=1 000 Hz

如图 5对本文方法与压缩等效源法在频率f=100～3 000 Hz范围内的声场重建效果进行了比较。算例中取 14个全息采样点，同样在测量声压时添加信噪比为20 dB的高斯白噪声。

图5 f=100～3 000 Hz时，本文方法和压缩等效源法的重建误差曲线Fig.5 Reconstruction error curve of the proposed method and C-ESM at f=100～3 000 Hz

由图5可以看出，随着频率的增加，两种方法的重建误差都随之增大，但在整个计算频率内，本文方法的重建误差均低于压缩等效源法，并能够保证准确地求解结果，体现了本文方法的优越性。

3.2 长条形声源的声场重建

为进一步验证本文方法的优越性，通过重建长条形声源的表面声压，将本文方法与压缩等效源法进行了比较。

仿真中，长条形声源的长宽高尺寸设为0.4 m×0.4 m×1.2 m，其表面声压由置于内部的 16个不同强度的点声源通过“替代法”生成。这些点声源布置在尺寸为 0.025 m×0.025 m×0.1 m 的长方体框架上，如图6所示。全息面设置为半径rH=0.8 m的球面，共布置366个测点。当使用压缩等效源法计算时，设置等效源面与声源面共形，其尺寸较声源面缩小为 0.16 m×0.16 m×0.48 m(缩小比例为1:2.5)，并在x和y方向均匀布置6个等效源、z方向均匀布置16个等效源，共352个等效源，如图6(a)所示。当使用本文方法计算时，为使等效源面尽可能地与声源面共形，将其设置为与全息面同心的椭球面，椭球面z向的长半轴a取0.3 m，x和y向的短半轴b、c均取0.075 m，如图6(b)所示。其中，抽样数M、N及求和截断项数mf、nf均与第3.1节一致。在仿真计算中，两种重建方法均添加信噪比为30 dB的高斯白噪声。

图6 长条形声源和两种方法的等效源布置示意Fig.6 Layout of long strip source and equivalent source of two methods

图 7给出了两种方法在频率 f=100～3 000 Hz范围内重建长条声源表面声压的误差曲线，其中重建误差可根据式(25)计算得到。从图7中可以看出，当采样点数足够时，两种方法在整个分析频段内均能够保证精确的重建结果，且误差都低于4%。

图7 采样点数为366时，两种方法的重建误差随频率变化的曲线Fig. 7 The curves of reconstruction error versus frequency for two methods with 366 sampling points

分析在采样点数减少的情况下两种方法的声场重建效果。图8给出了在114个采样点的情况下，两种方法在频率 f=100～3 000 Hz范围内重建长条声源表面声压的误差曲线。由图8可以看出，在100～2 000Hz频段内，两种方法均能得到准确的重建结果，且重建误差都低于 20%。随着频率的增加、当f=2 500～3 000 Hz时，压缩等效源法的误差逐渐增大，而本文方法仍保持着相对稳定的重建结果。这是因为本文方法是一种傅里叶级数展开形式的解析计算方法，只要求和截断项数nf、mf取得足够大，由式(21)所表达的声场可无限近似真实声场，同时也使得该式包含了辐射能力强的低阶模态和足够多的倏逝波模态，获得了足够的声场细节信息；待求的傅里叶展开系数向量还具有一定的稀疏性，满足了压缩感知理论中的稀疏性要求。因此在相同采样点数的声场重建中(尤其是在高频时)，本文方法的重建精度优于压缩等效源法的重建精度。

图8 采样点数114时，两种方法的重建误差随频率变化的曲线Fig. 8 The curves of reconstruction error versus frequency for two methods with 114 sampling points

图9给出了采样点数为 114、频率分别为1 500 Hz和2 500 Hz时上述两种方法的重建误差随信噪比变化的曲线。从图9中可以看出，当频率为1 500 Hz时，两种方法均可以保证较好的重建精度。但当频率增大时(f=2 500 Hz)，即使增大信噪比，压缩等效源法也得不到很好的重建结果，但本文方法的误差曲线在整个信噪比分析范围内(10～50 dB)始终低于压缩等效源法，并且重建误差始终保持在20%以下，表明本文方法不仅在相同采样点数的高频声场计算中获得了更高的重建精度，而且在不同的信噪比时均具有良好的抗噪性能。

图9 f=1 500、2 500 Hz，采样点数为114时，两种方法的重建误差随信噪比变化的曲线Fig.9 The curves of reconstruction error versus signal to noise ratio for two methods with 114 sampling points at f =1 500 and 2 500 Hz

4 结 论

本文在等效源法的基础上，将等效源强和积分核函数在轴对称的虚拟面上进行双向傅里叶级数展开，所建立的声压表达式是一种半解析半数值的形式，因而相对于其他数值离散方法具有更高的求解精度。另外，由于待求的傅里叶展开系数向量还具有一定稀疏性，因此在结合压缩感知的求解方法后，即使在少测点的欠定条件下也获得了满意的重建结果。文中利用数值仿真对比了所提方法与传统压缩等效源法的声场重建效果。结果表明：

(1) 对于重建脉动球声源，两种方法在整个计算频段下均有较好的重建效果，但所提方法的重建效果略优于压缩等效源法；

(2) 对于重建长条形声源，两种方法在 100～1 500 Hz的频段内重建效果相当，但在 2 500～3 000 Hz的高频段时，本文方法的重建精度均优于压缩等效源法；

(3) 在1 500 Hz和2 500 Hz频率时，本文方法的抗噪性能优于压缩等效源法。

虽然本文仅在轴对称虚拟面上将源强进行展开，但若进一步采用双向正交曲线坐标下的傅里叶展开，即可将本文方法推广到任意形状封闭声源的近场声全息计算中。