新型组合近场声全息技术

张子鑫， 肖友洪， 魏富康

(哈尔滨工程大学 动力与能源工程学院，黑龙江 哈尔滨 150001)

近场声全息技术(near-field acoustic holography, NAH)是一种有效的声场可视化技术，可以通过声场中的二维平面上测量得到的声压获得三维空间中的声学量。基于Fourier变换[1](fast Fourier transform, FFT)、边界元法[2](boundary element method, BEM)、Helmholtz最小二乘法[3](Helmholtz equation least squares, HELS)等的NAH已经得到广泛研究。然而当全息面孔径较小时，重建精度受到影响，原有的方法受到限制，局部近场声全息技术逐渐发展起来，如基于Fourier变换[4]和边界元法[5]的局部近场声全息技术。统计最优近场声全息[6](statistical optimal near field acoustic holography，SONAH)同样也属于局部近场声全息技术，通过平面波与倏逝波的线性叠加重建声场，突破了传统近场声全息技术的孔径限制。尽管类似SONAH的局部近场声全息技术不再依赖全息面的孔径，但算法仍受限于Nyquist采样定理，导致在中高频时即使采样面积较小，也需要大量采样点。

压缩感知技术[7]为上述问题提供了解决办法，将压缩感知与近场声全息技术相结合进行声场重建，可实现较少的采样点得到较好的重建效果[8-9]。本文在压缩感知和边界元法结合的基础上，运用SONAH进行全息面上声压的数据外推，并将其与未外推前的重建误差进行了对比分析。通过仿真与实验研究，验证了该组合近场声全息技术的有效性，并对影响因素进行分析，确定了获得较好重建效果时的外推面积大小。

1 组合近场声全息算法

1.1 边界元法与压缩感知

基于边界元法的近场声全息技术的理论基础为Helmholtz积分方程，若已知振动体表面的声压和法向振速，即可得到声场中任意点处的声压。将场点移到声源表面S，建立声源表面声压pS与法向振速vn之间的关系[10]：

HSpS=GSvn

(1)

利用Helmholtz积分方程与式(1)建立场点声压pf与振动体表面法向振速vn之间的关系：

pf=Gvn

(2)

vn=(GHG+λLHL)-1GHpf

(3)

式中λ为正则化参数，通常采用L曲线法[11]，广义交叉验证法[12]等方法进行选取。由全息面声压得到振动体表面振速后，利用式(2)即可得到声场中任意点处的声压。

压缩感知理论应用的前提是信号是稀疏的，或者在某一个基函数下的分解系数是稀疏的，可以通过一个与该基函数不相干的观测矩阵以远低于Nyquist采样定理的采样频率获取信号的全部信息。在边界元法中表面振速是需要获取的信号，选用声辐射模态作为声源表面振速的基函数[13-14]，其最早被用来分析辐射声功率。声源辐射声功率表示为：

(4)

(5)

(6)

式中矩阵U=[u1u2…uN]即可作为表面振速的基函数，此时，表面振速表示为：

(7)

式中：N为声源表面节点数，此时表面振速的求解过程可转化为先对分解系数α进行求解，因为声源表面振速在基函数U下的分解系数α是K稀疏的，且K≪N，所以可通过稀疏正则化进行求逆，将其转化为最小l0范数进行求解：

arg min‖α‖0， s.t.pf=GUα

(8)

同样，可以采用更加简单的最小l1范数代替最小l0范数，直接利用CVX工具箱进行求解。求得α后即可通过式(7)得到声源表面振速，再利用式(2)得到声场中任意点处的声压。

1.2 SONAH外推

SONAH的前提是将声场中任意点声压看作无数个空间波数域的平面传播波和倏逝波的叠加[15]，其数学表示为：

(9)

式中：M表示离散后的波数矢量总数；km=(kxm,kym,kzm)为第m阶的波数矢量；Φ(km,r)=exp[i(kxmx+kymy+kzmz)]为第m阶单元平面波。

SONAH外推的模型示意图如图1所示，假设全息面h上有n个测量点，外推后的全息面H上有N个测量点，根据波叠加原理，H上任意点Φ(km,rE)处的单元平面波Φ(km,rE)可由h上同阶的单元平面波线性叠加得到：

图1 模型布置Fig.1 Diagram of model layout

(10)

式中：ci(rE)为全息面h上第i个测量点处各阶单元平面波的权重系数。离散后的波数矢量总数为M，由此可得到M个线性方程，组成矩阵形式为：

b=Ac(rE)

(11)

通过正则化方法，可得到系数向量的解为：

c(rE)=(AHA+λLHL)-1AHb

(12)

联合式(9)、(10)和式(12)可以将全息面H上一点rE用全息面h上各点声压表示为：

p(rE)=phc(rE)=ph(AHA+λLHL)-1AHb

(13)

本文提出的组合近场声全息方法，即基于外推后的全息面H上各点的声压，并利用边界元法结合压缩感知，实现了较少场点情况下的声场重建。

2 板声场重建仿真

仿真中，选用0.8 m×0.8 m×0.003 m的铝板作为声源，边界施加固支约束，以铝板中心为坐标原点建立模型，分别在(-0.1,-0.18,-0.001 5) m和(0.03,0.165,-0.001 5) m处施加大小为10 N的激励力。全息面h选在距原点高0.04 m处，大小为0.64 m×0.64 m，离散后的点数为81个，阵元间距为0.08 m。重建面S选在距离原点高0.02 m处，大小为0.8 m×0.8 m，重建点个数为121个。

为了定量评价该方法的结果，误差被定义为：

e=‖pr-pt‖2/‖pt‖2

(14)

式中：pr表示算法得到重建面的声压；pt表示仿真得到的重建面声压理论值。

2.1 所提方法的验证

全息面h上场点数为9×9，外推后全息面H大小为0.8 m×0.8 m，场点数为11×11，利用外推后全息面H上的场点声压，通过边界元法进行声场重建，分别比较了边界元法与传统正则化方法和压缩感知方法相结合的重建误差，并将结果与未外推前的重建误差进行比较。结果如图2所示，图中TR与CS分别为传统正则化方法与压缩感知方法的缩写。

图2 采用组合近场声全息前后重建误差的比较Fig.2 Comparison of reconstruction errors before and after using combined near-field acoustic holography

场点选择81个时，使用压缩感知与边界元法相结合的重建误差整体小于传统正则化方法，说明前者对于较少场点的重建效果较好。但是随着频率的增加，4种情况下的重建误差都随之增加，这是因为频率的增加导致波长变短，倏逝波衰减迅速，声场信息减少，从而降低了重建精度。比较利用外推前后的声压重建声场的误差，可以发现经过外推，重建误差明显降低，且在整个分析频段内重建误差均小于30%，验证了外推方法与边界元法结合后的算法的有效性。但是从图中可以看出，利用外推后的场点声压进行声场重建时，传统正则化方法与压缩感知方法的重建误差十分接近，而且在个别频率传统正则化方法的重建误差略小于压缩感知方法的重建误差，说明外推后的声压产生的误差对传统正则化方法的影响更小。

图3分别展示了1 050 Hz时的重建面声压的理论云图与利用组合近场声全息技术前后重建面的声压云图。从云图中可以看出，测量场点较少时，利用边界元法结合传统正则化方法并不能对声场进行有效重建，并且重建的声压幅值也大于理论值，而虽然利用边界元法结合压缩感知改善了重建结果，但在中心位置重建的幅值较小。利用组合声全息方法，先对全息面声压进行外推，可以对场点声压进行有效地补充，可以看出重建结果与理论值吻合较好，较外推前利用传统正则化方法的结果有很大的改善，相较于少场点时的压缩感知方法，外推后的结果在中心位置的幅值也有所增加，验证了该方法的有效性。

图3 1 050 Hz外推前后重建面声压云图比较Fig.3 Sound pressure nephogram of reconstructed plane before and after extrapolation at 1 050 Hz

2.2 外推点数的影响

这一节讨论阵元间距相同，外推点数不同对重建效果的影响。仿真模型与上一节相同，外推后全息面H的尺寸分别选为0.8 m×0.8 m、0.96 m×0.96 m、1.12 m×1.12 m，对应的外推后的场点分数别为11×11、13×13、15×15。首先计算了3种情况下外推得到的声压与仿真得到的声压之间的误差，其结果如图4所示。

图4 不同外推点数的估计误差比较Fig.4 Estimation error of different extrapolation points

图4中可以看出随着外推点数的增加，外推的整体误差也在增加，且在外推点数为15×15时，400 Hz以上的频率外推误差已经超过40%，此时外推方法已经失去应有的作用，这是外推方法的局限性。外推方法本质上是通过已知的场点声压来对未知的场点声压进行估计，外推的距离越大，预测场点离已知场点越远，声压的变化就越大，外推也就失去效果。而外推点数为13×13时，低频时的外推误差较小，而在中高频时外推误差则逐渐增加，这是因为随着频率的增加，板的振动情况逐渐复杂，外推点的声压与已知点的声压变化较大，导致外推误差增加。

将外推后的声压作为输入计算重建面的重建误差，因为外推点数为15×15的外推误差在中高频时已经超过50%，再用其进行重建也没有意义。因此在此只比较外推点数为11×11和13×13的重建误差，重建算法选择压缩感知与边界元法的结合，重建误差的结果如图5所示。在200 Hz以下时外推点数为13×13的重建误差更小，随着频率增加，各点声压的估计误差也增加，此时引入后续重建的估计误差的比重超过了其对场点信息的补充，导致重建误差也随之增加。因此对于外推方法来说，将外推点数限制在沿x，y方向分别外推两个阵元间距时，可得到较好的重建结果，有效降低重建误差。

图5 不同外推点数的重建误差比较Fig.5 Reconstruction error of different extrapolation points

2.3 外推-内插与仅外推方法的比较

前2节主要讨论了外推前后阵元间距不变时，组合近场声全息技术的优势以及影响因素。本节将外推后阵元间距变为0.04 m，则全息面H上的场点数为21×21，大小为0.8 m×0.8 m，此时是在外推的基础上同时进行了内插，二者运用的方法相同，仿真模型布置图如图6所示。

图6 外推-内插仿真模型布置Fig.6 Diagram of extrapolation-interpolation

首先比较外推-内插与仅外推时，全息面H上的声压的估计误差，结果如图7所示，在800 Hz以下时，外推-内插的误差小于只进行外推的误差，而在频率大于1 200 Hz时，外推-内插的误差却急剧增加。这是因为在低频时板的振动情况较为简单，对全息面h上的数据进行外推-内插，虽然需要估计的场点增多，但是各点估计值较为准确，而在高频时板的振动情况变得复杂，各点估计误差累积导致整体误差增大。1 300 Hz时外推-内插的误差达到了63.8%。

图7 外推-内插与仅外推方法的估计误差比较Fig.7 Estimation error comparison between extrapolation-interpolation and extrapolation methods

图8中分别展示了1 300 Hz时只进行外推的全息面H的理论和外推声压云图，以及进行外推-内插时的全息面H的理论与外推声压云图。相较于121个场点，441个场点时声压云图的分辨率更高。

图8 1 300 Hz时全息面理论与外推声压云图Fig.8 Holographic plane theoretical and extrapolated sound pressure nephogram at 1 300 Hz

从81个场点经过外推-内插得到441个场点的声压，虽然云图内幅值较大位置没有发生偏移，但与理论值相比略显粗糙，因此估计误差远大于仅外推时的误差。

利用拓展后的全息面声压分别通过边界元法结合传统正则化方法与压缩感知方法进行声场重建，重建的误差如图9所示。尽管在1 150 Hz以后，相较于仅外推时，外推-内插对全息面声压的估计误差急剧增加，但是其重建误差却与仅外推时的误差基本相同。出现上述现象的原因在于与理论值相比，外推-内插得到的全息面声压幅值较大的位置能够较为准确地显示出来，且重建面的点数为121个，其的分辨率小于全息面的分辨率。而在1 000 Hz以下的中低频，利用外推-内插得到的声压通过边界元法结合压缩感知方法所得的重建误差在4种情况下最低。但外推-内插在丰富了场点信息的同时，也为后续传递矩阵的计算增加了时间成本。

图9 外推-内插与仅外推时重建误差的比较Fig.9 Reconstruction error comparison between extrapolation-interpolation and extrapolation methods

3 板声场重建实验

为验证算法的性能，针对板声源展开实验研究，实验布置如图10所示。将0.8 m×0.8 m×0.003 m的板固定在箱体上，箱体内悬空设置一个激振器对板施加激励，使其振动作为目标声源，激振器的信号频率设置为100～600 Hz，间隔为100 Hz，分析不同频率下算法的重建结果。全息面与重建面的大小均为0.66 m×0.66 m，位置分别距离板声源0.04 m和0.02 m。全息面与重建面声压均由6×6的麦克风阵列分4次测得，点数为144个，另设参考麦克风保证测量声压的相位，相邻麦克风之间的间距为0.06 m。测量结束后，提取全息面中心64个点作为外推前的场点，组合成全息面h，大小为0.42 m×0.42 m。选取外推后全息面H的大小为0.54 m×0.54 m，场点数为100个。保持全息面H的大小不变，减小阵元间距，则外推-内插的场点数为361个。

图10 板实验布置Fig.10 Experimental setup of plate

图11为600 Hz时采用组合近场声全息技术前后的重建面的声压分布云图，并将其与实验测得的云图进行比较。从图中可以看出，选取64个场点时，全息面h面积较小，虽然在板中心的声压分布情况可以重建，但是对于边缘处则由于缺少必要的场点信息而无法得到较好的重建效果。而采用组合近场声全息技术以后，分别对全息面h上的场点声压进行了外推以及外推-内插，从图中可以看出，在重建面边缘位置可以得到较好的重建结果，且相较于仅外推时，采用外推-内插后重建结果的幅值与测量值更为接近，重建效果有所改善。

图11 600 Hz重建面声压云图Fig.11 Sound pressure nephogram of reconstructed surface at 600 Hz

图12给出了各频率下采用组合近场声全息技术前后的重建误差，可以看出除了100 Hz，其他频率处经过外推后重建误差都有所降低，且随着频率增加，误差降低的幅度也随之增加。和前文分析相一致，低频时较小的全息面即可进行声场重建，在中高频时波长变短，场点包含的信息变少，且板的振动情况也随之变得复杂。此时，采用外推方法可对全息面声压进行有效补充，重建误差也随之降低。采用外推-内插则是在外推的基础上更加丰富了场点信息，在中低频时可改善测量结果。

图12 采用组合近场声全息前后重建误差比较Fig.12 Comparison of reconstruction errors before and after using combined near-field acoustic holography

4 结论

1) 组合近场声全息技术针对传统近场声全息技术少场点测量重建精度不高的问题，通过SONAH进行声场外推，将边界元法与压缩感知相结合，有效降低了声场重建误差。

2) 采用组合近场声全息技术，分析外推点数与外推面积的选定，通常沿x，y方向分别外推2个阵元间距，即可得到较好的重建效果。

3) 相较于对全息面仅外推时，外推-内插的重建结果在中低频时有一定的改善，但同时增加了计算的时间。