基于稀疏贝叶斯学习的高分辨率Patch近场声全息

扈 宇， 胡定玉， 方 宇， 肖 悦

(1. 上海工程技术大学 城市轨道交通学院，上海 201620；2. 南昌工程学院 江西省精密驱动与控制重点实验室，南昌 330099)

基于空间二维Fourier变换的近场声全息[1-2](Near-field Acoustic Holography, NAH)技术能够充分利用包含高空间频率的倏逝波信息实现高分辨率的声场重建，是一种重要的声场可视化与噪声源识别工具。但NAH的极限分辨率取决于全息面上的测点间隔，因此要获得较高的分辨率，就需要极高的测量成本，这极大影响了该技术的实际应用价值。针对这种问题，Harris等[3]提出了利用声压梯度为指导的全息声压的三次Hermite插值法提高声场重建的分辨率，但该方法需要测量质点振速信息，增大了测量工作量。Xu等[4]提出了利用带限信号恢复算法的NAH分辨率增强方法，但该方法需要经过多次迭代计算，增加了计算时间。随后，徐亮等[5]提出了一种基于正交球面波插值的方法，但是对于形状较为复杂的声源插值效果并不理想。张小正等[6]通过波叠加法原理构建虚源面对全息数据插值，一定程度缓解了声源形状对重建结果的影响。最近，毛荣富等[7]利用支持向量回归算法实现了全息数据插值，有效地减少了测量工作量。孙超等[8]提出了一种基于极限学习机的全息数据插值方法，通过神经网络训练数据提高全息分辨率并对全息数据外推。

基于空间二维Fourier变换的NAH技术在应用中遇到的另一个问题是测量孔径对NAH技术重建精度的影响[9-10]。由于NAH的理论是建立在测量面为无穷大的基础上，因此当测量孔径较小时，窗函数的截断效应会导致波数谱泄露等一系列影响重建精度的问题。Saijyou等[11]在研究大尺寸声源重建问题时通过全息面补零和波数域反复迭代滤波有效抑制了测量孔径的影响。Williams等[12]进一步提出了基于空间二维Fourier变换的Patch NAH，以解决小测量孔径下的全息重建问题。徐亮等[13-14]在波数域进行全息数据外推并通过极小化加权范数优化了外推方法，Scholte等[15-16]对重建过程的低通滤波方法进行优化并利用线性边界预测法实现了空间域上的全息数据外推，有效提高了小孔径全息数据的重建精度。在全息数据外推方法不断发展的同时，基于波叠加法[17]、分布源边界点法[18]、统计最优法[19]与边界元法[20]的Patch NAH技术相继被提出，扩大了Patch NAH的适用范围。

本文针对基于空间二维Fourier变换的NAH重建精度受到测量孔径影响较大的问题，提出一种基于稀疏贝叶斯学习(Sparse Bayesian Learning, SBL)算法的高分辨率Patch NAH方法。该方法利用测点的空间坐标构建高斯核函数，通过SBL算法构建全息面声压拟合模型实现全息面声压数据的插值和外推，并进一步将插值和外推的数据用于NAH重建，有效抑制了测量孔径与测点间隔对NAH的限制。本文首先给出基于SBL的全息面声压插值与外推理论模型，然后通过数值仿真对该方法用于Patch NAH重建的有效性、抗噪性与分辨率提升等性能进行分析，最后通过实验进一步验证该方法的有效性。

1 基于空间二维Fourier变换的NAH

基于空间二维Fourier变换的NAH的基本原理是利用均匀分布在靠近声源表面的全息面上有限个测点的声场信息，通过全息面和重建面声压在波数域的传递关系重建声源辐射区域任意平面的声压及其他声学物理量。如图1所示，波数域声压从全息面到重建面的传递关系为

P(kx,ky,zr)=P(kx,ky,zh)e-ikz(zh-zr)

(1)

(2)

由于在实际计算中采用了快速Fourier变换算法，且测量孔径与测量点数有限，因此重建分辨率受到测量点数的影响。另外将式(2)分解可得到

(3)

在重建过程中通常会利用窗函数抑制噪声的放大，其中最常用的经典窗函数是Maynard等[21]提出的指数窗，其形式为

(4)

图1 声源面、重建面与全息面分布图Fig.1 A diagram of NAH

2 基于SBL算法的全息面声压插值与外推

2.1 高斯核函数

高斯核函数是径向基函数的一种，二维高斯核函数的数学表达式为

(5)

2.2 基于SBL算法的全息面声压插值

高斯核函数可以将低维数据映射到高维并对其进行线性求解，在升维过程中保证了以高斯核函数为基的解的稀疏性，理论上可以通过较少测点拟合声压在全息面上的分布，从而实现超高分辨率[24]。由于在全息重建中使用的是复声压，且复声压实部与虚部的误差虽都服从零均值高斯分布，但方差可能会有不同，因此可以分别对复声压实部与虚部插值并将两者合并获取复声压。下面仅以实部插值为例介绍全息面声压插值过程。

通过高斯核函数映射后，全息面上任意一点ri处的声压实部理论值可由高斯核函数表示为

(6)

pre=Φw+ε

(7)

(8)

(9)

式中：α=[α1,α2,…,αN]T，αj为权重系数wj的超参数。

在α与β2皆为未知的情况下， 通过测量的声压求解权重系数w的问题即可转化为最大化联合后验概率Prob(w,α,β2pre)的问题， Prob(w,α,β2pre)最大值对应的w即为最优的权重系数。 但由于Prob(w,α,β2pre)很难通过积分直接获取， 因此将Prob(w,α,β2pre)分解为

Prob(w,α,β2pre)=
Prob(wpre,α,β2)Prob(α,β2pre)

(10)

式中：w的后验概率分布Prob(wpre,α,β2)可以根据贝叶斯定理表示为

(11)

式中： 分母Prob(preα,β2)可由似然函数式(8)对权重系数w进行边缘积分得到，即

(12)

式中：A=diag(α1,α2,α3,…,αN)，C=β2I+ΦA-1ΦT。

将式(8)、式(9)和式(12)代入式(11)可得

Prob(wpre,α,β2)=N (μ,Σ)

(13)

式中：Σ=(β-2ΦTΦ+A)-1，μ=β-2ΣΦTpre， Prob(wpre,α,β2)满足均值为μ、方差为Σ的N维高斯分布。

得到w的后验分布Prob(wpre,α,β2)后， 稀疏贝叶斯学习过程即被转化为超参数后验分布Prob(α,β2pre)∝Prob(preα,β2)Prob(α)Prob(β2)关于α与β2的最大值问题。 在贝叶斯模式中， Prob(preα,β2)被称为边缘似然函数。 对边缘似然函数进行第Ⅱ类型最大似然参数估计可求出α与β2的估计值， 即利用式(12)分别对α与β2求偏导并使导数等于零， 可求出α与β2的更新公式

(14)

(15)

2.3 基于SBL算法的全息面声压外推

全息面声压外推过程和插值过程类似，都是利用高斯核函数实现，区别在于外推点的位置和插值点位置的不同。由于外推点均位于测量孔径外部，而全息面声压是采用以测点位置为中心峰值位置的高斯核函数进行拟合，因此根据高斯核函数的衰减特性，测量孔径外的外推点声压平滑衰减，最终可得到边缘为零或无限近似于零的全息面声压数据，因此利用外推后的数据进行重建，可以极大避免窗函数截断效应的影响。需要说明的是，在应用过程中可以同时进行全息面声压的插值和外推，可以在几乎不增加计算时间的同时提高重建结果的分辨率。

3 数值仿真

(16)

式中：Pcal为插值或重建的声压，Pth为对应的理论值。

3.1 全息面声压插值

为便于描述，定义插值倍数为插值前后分辨率的比值。对全息面声压进行5倍插值，即插值后全息面声压的分辨率为1 cm。图2显示的是频率为900 Hz时全息面上的测量声压(图2(a))、插值后的声压(图2(a))及其理论值(图2(c))。从图中可以看出，插值结果与理论值几乎完全一致，利用式(16)计算出插值后全息数据的插值误差仅为1.88%，证明了本文插值模型的有效性。

图2 900 Hz时全息面声压插值效果对比图Fig.2 The pressure on the hologram at 900 Hz

为研究插值倍数与插值精度之间的关系，对全息面声压进行不同倍数的插值。频率为900 Hz时插值误差与插值倍数的关系如图3所示。从图3中可以看出，虽然插值点间距缩小使得全息数据点数成倍增长，但插值误差的变化很小，说明插值点数对于插值结果的影响很小。图3同时也显示了全息面声压的测量数据与理论值之间的相对误差。比较测量误差和插值误差可以看出，基于SBL算法的全息数据插值方法对测量数据中的噪声有一定的抑制能力。图4给出了频率为500 Hz和900 Hz时，在不同的信噪比(SNR)的情况下插值误差与测量误差的对比。从图中可以看出，在测量条件较差的情况下，使用插值模型可以较好地抑制噪声，证明了插值模型具有较高的工程实用性。

图3 插值误差与插值倍数的关系Fig.3 The relative errors of the measurement and theinterpolation versus the ratios of the spatial resolutionsbefore and after interpolation

图4 插值误差与测量误差随信噪比变化曲线Fig.4 The relative errors of the measurement and theinterpolation versus the SNR

3.2 高分辨率声场重建

对全息面声压进行5倍插值并同时进行外推，并利用插值和外推后的全息面声压重建距板表面0.02 m处的声压。图5显示的是频率为900 Hz时重建结果与理论值的比较。从图中可以看出，重建的声压与理论值吻合地很好，通过式(16)计算两者之间的误差仅为3.32%，由此可见，对全息面声压进行插值和外推处理可以有效地避免有效孔径效应对重建精度的影响，同时可以获取高分辨率的重建结果。

图6显示的是在100～2 000 Hz的频率范围内声压的重建误差随频率的变化。从图中可以看出，声压的重建误差均在10 %以下，证明了本文方法的有效性。

图5 900 Hz时重建结果与重建面处实际声压幅值分布对比图Fig.5 The pressure on the reconstruction plane at 900 Hz

为进一步验证高分辨率Patch HAH对小测量孔径的适应性，在全息面上取三个不同大小的测量孔径，其位置如图7所示。分别对三个孔径的全息数据进行插值和外推，并进一步利用NAH重建距板表面2 cm处的声压。图8显示的是频率为900 Hz时利用三个测量孔径内声压重建的结果与其理论值的比较。从图中可以看出，三个测量孔径的边缘附近均存在较大峰值，测量孔径产生的截断效应会对基于二维空间Fourier变换的NAH的重建精度产生很大影响。而利用本文方法对全息面声压进行插值和外推后，利用三个测量孔径的全息数据都获得了较好的重建结果。其中，三个孔径的重建误差分别为孔径1∶10.95%，孔径2∶15.53%和孔径3∶16.35%。证明了本文插值外推方法对较小的测量孔径仍具有较好的稳定性。

图7 全息面测量孔径分布位置图Fig.7 Positions of the measurement apertures

图8 900 Hz时各孔径重建结果对比图Fig.8 The pressure on the three reconstruction plane at 900 Hz

4 实验验证

实验在半消声室中进行，半消声室本底噪声小于16 dB。实验装置实物照片如图9所示。声源为80×66 cm2的固支钢板，钢板厚度为2 mm。对钢板施加简谐信号产生声场。全息面位于钢板上方5 cm处，测点间隔为5 cm，共15×11个测点。重建面位于钢板上方2 cm处。为检验重建结果，同时以2.5 cm的测量间隔测量重建面上的声压作为其真实值，共31×21个测点。在15×11个测点得到全息面中，选取2个大小不同、测点数不同的测量孔径并利用孔径内的数据进行重建，孔径在全息面上的位置分布如图10所示。算法迭代过程中判断收敛条件与仿真相同。

图9 实验装置实物图Fig.9 A photo of the experimental setup

图11显示的是频率为700 Hz、900 Hz与1 100 Hz时，对两个测量孔径的全息数据进行外推并插值后重建结果与真实值的比较。从图中可以看出，利用本文方法对不同尺寸测量孔径进行重建皆获得了较好的重建效果。不同激励下各测量孔径重建误差见表1。重建结果表明，本文方法可以准确地重建出声场分布，并且有效地提升了重建结果分辨率，证明了利用SBL算法对全息数据插值与外推后可以实现高分辨Patch NAH。

图10 全息面测量孔径位置Fig.10 Positions of the measuring apertures

激励频率/Hz7009001 100孔径1重建误差/%13.8810.3218.69孔径2重建误差/%7.7414.1511.94

图11 不同激励频率下各孔径重建结果对比图Fig.11 Comparison of reconstruction results with the corresponding true values

5 结 论

本文建立了一种基于SBL的高分辨率Patch NAH方法，解决了测量孔径与测点数对基于空间二维Fourier变换的NAH重建效果影响较大的问题.研究结果表明：利用高斯核函数和SBL算法对全息面数据进行插值和外推，可以有效地避免测量孔径对NAH重建结果的影响，同时获取高分辨率的重建结果。另外，利用插值模型还可以在一定程度上减少测量噪声，具有较高的工程应用价值。