声全息声压场插值重构方法研究

赵玉贵，刘志红，仪垂杰，宿元亮

（1.机械与汽车工程学院青岛理工大学，山东 青岛 266520；2.教育部重点实验室“工业流体节能与污染控制”青岛理工大学，山东 青岛 266520）

1 引言

声全息技术是一种可同时识别声源位置和传播性能的有效手段，它利用声源发出声波的幅值和相位信息获取声场的声压空间分布，依据声压幅值大小识别声源位置，声位置矢量确定声传播性能[1]，其测量过程，如图1所示。为提高声压的测量和重建精度，需在无限大平面上布置多个麦克风完成，而实际测量过程只能在有限空间进行，即将无限大平面用有限采样面截断，从而导致测量及重建误差，声压幅值精度下降[2]。如何在不增大测量工作量及经济成本情况下，有效减少测量重建误差提高声源识别精度是解决声全息技术高效应用的关键问题，同时该问题得到学者们的高度关注。

在众多专家学者的研究下，Patch近场声全息[3−8]以及插值法等多种方法应运而生，其中，插值法因其计算便捷和扩展性备受青睐，2006年文献[9]利用Hermite 插值对全息面声压重建大幅提高声全息声源的空间分辨率，但该方法除需要每个测点声压外还需测量每点的振速，在一定程度上增加了测量和计算的工作量。文献[10]提出了一种基于带限信号恢复算法的NAH分辨率增强方法，利用声压的波数域带限特性，通过带限信号恢复算法（PGA），实现全息面声压插值提高声源识别精度，但该方法需多次迭代计算。2007年文献[11]，2015年文献[12]提出了基于正交球面波插值的NAH分辨率增强方法，以实际测量点数据为插值条件，通过若干不同阶次的球面波源叠加拟合实际声场实现全息面插值，在减少测量点数同时等效增加了全息面声压数据。不过这种方法因采用无限阶正交球面波源叠加的有限截断增加了实际拟合误差，且插值过程需将球面波源布置在声源内部而往往声源位置是无法预知，限制了该方法的实际应用。为提高插值法的应用性，2010年文献[13]等提出虚源插值法，通过在全息面附近布置虚源面，利用实测全息面声压求得虚源面上简单源源强，再由简单源源强对全息面声压插值，利用插值数据进行声压重建。综述所述，插值法能在不增加测量工作量条件下，有效提高声源识别精度，但该方法仍存在急需解决的问题：（1）如何高效利用实际少量测量数据，高精度插值出任意位置上的声压幅值，且降低过程计算量；（2）如何减少测量点位置对声压插值的影响，提高全息面上声压重建的有效性。

提出一种融合牛顿插值和克希荷夫积分的新插值方法。该方法以少量实测点数据为对象，由二维牛顿插值法对虚拟位置声矢量和声压进行插值，将离散的测量值转换为连续函数值，再通过重采样得到大量离散点，由克希荷夫积分得出克希荷夫面的声势，最后由上述测量及插值数据重构全息面声压值。

2 声全息测量与插值重构

2.1 声全息测量原理

声全息声源识别原理，如图1所示。全息面H是麦克风布置平面，重建面S是声源面，Z0为全息面与重建面间距离，R为全息面麦克风位置声矢量，通过测量得到H面上声压分布函数U(x,y)，根据式（1）计算得出S上的声压分布函数U(a,Z)，由其幅值大小识别出声源位置[14−15]。

图1 声全息测量原理图Fig.1 Acoustic Hologram Schematic Represented by a Single Right−Angle Coordinate System

重建面S上的声压分布函数U(a,Z)，由式（1）得出：

式中：U(a,Z)—重建面S上任意一点Q(x,y,zS)的声压值；k0—波数；

C—全息系数为常值，简化计算这里取1；

r1—重建源面上的任意点Q(x,y,zS)与全息面上的麦克风位置点P(xH,yH,zH)间距离：

将式（2）代入式（1），且将U(a，Z)写作US(x，y)整理：

这里：

r0=点到重建面上任意一点Q(x，y，zS)的距离。

由式（3）可知，精确确定全息面上的声压分布函数UH(x，y)是重建声压分布函数US(x，y)的关键。

2.2 插值重构

为获得全息面声压函数UH(x，y)，首先，保证测量麦克风的测试数据的准确性即测得大量的声源有效声信息。阵列的合理构型是保证阵列获取大量有效信息的前提[16]。

这里采用同轴圆环阵列作为测量阵列，如图2所示。其构型方式为：第一环半径R为入射波最小波长的一半即R=λmin2，同一环上麦克风间距Δr，λmin4 ≤Δr≤λmin2。不同环间半径差ΔR，0 ≤ΔR≤λmin2。其次，高效重建非测量位置的声压幅值。非测量位置点即全息面上声压待建点的数量直接影响全息面声压函数的重建精度。为此，以同轴圆环阵列面为全息面，声压待建点的确定沿用阵列麦克风位置的设置原则，声压值由牛顿插值和克希荷夫积分实现。最后，对所有非测量位置点和测量点上的声压计算处理后得出UH(x,y)值带入式（3）重建声源面上的声压分布函数实现声源识别。

图2 极坐标同轴圆环全息面示意图Fig.2 Polar Coordinate Ring Hologram

由Newton插值算法计算同环上的插值点声压和不同环上的声压待建点声压。首先，以单环阵列为测量阵，构建同一环上的插值点位置的声压。取第1环R上的麦克风测量声压数据及位置矢量θi为变量，构建1 环上插值点的声压与位置矢量的函数关系，即任给一个位置矢量θ可得出该位置的声压幅值P。其P(θ)的各阶差商计算过程：

函数P(θ)关于点θ0,θk的一阶差商表示为：

将函数P(θ)关于点θ0,θ1,θk的二阶差商表示为：

以此类推，得到函数P(θ)的k阶差商：

根据差商的定义，把θ看作[0，π]上的一点，可得：

只要将（8）中后一式依次代入前一式，得到最终的插值公式：

由式（6），代入麦克风测量声压数据，即位置矢量θi和声压值P(θi)，便可得出第一环R1上任意位置的声压幅值，其它环上的离散点声压可由克希荷夫积分式（10）计算得到，此时，设置第一环为虚拟克希荷夫面。

式中：S1—虚拟克希荷夫面；PS1—虚拟克希荷夫面S1上的声压函数—法线；G—自由空间格林函数；PSk—第k环上的声压函数。

PS1可由第一环上各插值点均值表示：

式中：Pi—插值后第i点处声压值，而ωi为第i点处的加权函数，表示各点对整个声压的贡献大小，此处ωi=1n。而式（10）中的格林函数G则可以表示为：

式中：rk—声源到测量点的距离。

为便于计算，将式（12）在球坐标系下展开[17]：

式中：jn(kr<)—球贝塞尔函数，其中，r<是取r,r0中的较小值；hn(kr>)—勒让德函数，式中的r>是取r,r0中的较大值。将式（12）代入式（10）得

其次，采用同心圆环阵为测量阵列，即采用多个环上的麦克风测量值为插值数据，提高测试数据的有效性和声压的重建精度。

其计算过程与单个环的一致，只是将克希荷夫面数量变为三个，并将三个克希荷夫面对面外声压待建点的贡献度进行加权计算。

最后，全息面重构方法流程，如图3所示。

图3 全息面全息函数获取流程图Fig.3 Holographic Surface Holographic Function Acquisition Flow Chart

3 数值仿真与实验分析

3.1 插值算法数值仿真分析

为探讨上述方法的有效性和关键影响因素进行如下数值仿真。仿真数据：重建面S的中心点布置一点声源，声压级*70dB，在平行于重建面且间距为Z=0.33m 和Z=2m 的位置分别布置一个3环同轴圆环阵，阵列中心与重建面中心等同高，环半径差取为ΔR=0.3λ，麦克风数量第1环8个，第2环16个，第3环20个。全息面为（1×1）m，布置9个环，声压待建点为243个，入射波频率为2000Hz。仿真过程一圆环阵为测量阵列，阵列构型及插值点，如图4所示，全息面上的待建点位置，如图5所示。

图4 测量麦克风阵列及插值点Fig.4 Microphone Measurement Array and Interpolation Points

图5 全息面声压待建点Fig.5 The Point on the Holographic Surface where the Sound Pressure is to be Reconstructed

实验1：首先，为验证插值点声压重建精度，由阵列测量数据，经插值和克希荷夫积分重构计算出每个环上插值点声压幅值，然后，在相应插值点位置放置麦克风测量，将计算声压pin与实测声压比较，其误差评价函数为：

大量仿真统计分析发现：声波入射频率、测量阵列构型、麦克风数量和声压待建点）数量是影响全息面上声压重建精度的主要因素。

因全息面上声压待建点布置方式是按阵列构型中麦克风位置的布置方式确定，所以关键因素包括同环上的麦克风间距，它直接影响各单环插值点的牛顿声压插值精度，和相邻环的环间距，它改变了声压待建点位置及数量对整个全息面的声压重建影响显著。分别针对1000Hz，2000Hz，3000Hz三个频段，不同麦克风间距Δr和环间隔ΔR仿真，其结果，如表1，表2和图6、图7所示。

表1 不同频率条件下各测点间隔对应测量麦克风数目Tab.1 The Number of Measuring Microphones Corresponding to Each Measuring Point Interval Under Different Frequency Conditions

表2 不同频率条件下不同环间隔对应的全息面声压待建点数量Tab.2 Number of Holographic Surface Acoustic Pressure Points to be Built Corresponding to Different Ring Intervals Under Different Frequency Conditions

图6 不同频率条件下，Δ与插值误差关系Fig.6 The Relationship Between Δ and Interpolation Error Under Different Frequency Conditions

由图6 可知，随着麦克风间距增大，插值误差也在增大，当d>λ3（Δλ>0.3）时的插值误差变大速度快，d<λ3时误差增速平稳，尤其低频1000Hz 时显著。所以，麦克风间距在λ4 ≤d≤λ3范围选取效果较好，当频率高于1000Hz时，范围适当放宽到λ4 ≤d≤λ2。

由图7可知，同轴圆环面上两相邻环间隔过小或过大均引起较大误差，当ΔR在λ/3~λ/2范围内，积分误差相对较小，所以ΔR一般在该范围内取得。

图7 不同频率条件下，ΔR与积分误差关系Fig.7 The Relationship Between ΔR and Integral Error Under Different Frequency Conditions

实验二：其次，声源重建面的重建结果。

3.2 插值算法有效性验证

3.2.1 实验方法

该全息测量实验在半消声室内进行，室内本底噪声小于20dB，单环测量时实验装置实物布置，如图9所示。

图9 全息测量装置图Fig.9 Holographic Measuring Device Diagram

图8 声源面上重建声源（声源幅值归一化）Fig.8 Reconstruction of Sound Source on the Sound Source Surface（Normalization of Sound Source Amplitude）

在源面上布置一个声源，声源选择一个尺寸为（146×170×223）（单位mm）的音箱，发出频率为f=2000Hz的单频信号，音箱纸盆中心坐标为（0，0，0），在平行于源面且间距分别为0.33m 和2m的位置布置一测量面作为全息面，重建面距离测量面距离为0.1m，全息面与重建面尺寸均为（1×1）m。本实验所用测量面为同轴圆环传声器阵列，阵列中心坐标分别为（0，0，0.33）和（0，0，2），平面内最小环为第一环，直径为D=0.17m，其他环间距为ΔR=6cm。在测量面第一环上均匀布置8个测量麦克风，第二环及第三环上分别布置16和20个测量麦克风（θ=0°除外），圆环中心布置1个参考麦克风。

该实验分别在不同距离上进行测量，为了对插值算法的有效性进行验证，两次测量均进行两组对照实验。

第一组实验对三环所测数据采用Newton插值及克希荷夫积分处理得到非测量环的其他环上声压待建点的声压数据；第二组则采用测量麦克风对全息面上所有环上的声压待建点进行实际测量测得声压数据。将得到的数据进行声源声压重构生成声压幅值图，最后就重建面位置处的声压实测值与声压计算值进行比较分析。

3.2.2 实验结果

（1）0.33m条件下重建结果

（2）2m条件下重建结果

在0.33m与2m距离条件下得到的重构后声压幅值图，如图10、图11所示。图10（a）是0.33m距离情况下依据Newton插值结合克希荷夫积分公式得出的全息面声压数据重建出的声源声压幅值图；图10（b）则是根据8环233个麦克风实际测量值重建的声压。在2m距离情况下的声压幅值图，如图11所示。

图10 0.33m处声全息重建声压幅值图Fig.10 Sound Pressure Amplitude Map of Acoustic Hologram Reconstruction at 0.33m

从图10与图11中可以看出，两种方法重构所得的声压幅值图外形差别较小；再进行数据比对，对重构面上若干点声压进行实际测量，再与两种方法的对应点重构的计算值进行比较，用式（15）算得各点重构计算所得声压数值与实测值的相对误差，误差均小于5%，在可接受范围之内，可以用来替代原始数据进行重建，且使用Newton插值与克希荷夫积分在重建效果相差不大情况下计算则更为简单。

图11 2m处声全息重建声压幅值图Fig.11 Sound Pressure Amplitude Map of Acoustic Hologram Reconstruction at 2m

上述实验结果同样可以验证上述插值方法能很好地实现全息面声压重建工作，能够节约测量时间、简化测量工作。

3.3 误差分析

利用牛顿插值与克希荷夫积分公式相结合的复合算法实现声全息声压场的重建的误差主要存在于以下方面[18]：

这里主要采用插值算法来实现声压场的重建，主要是为了减少测量的时间及成本，但插值所得数值与实际数值存在一定的偏差，这种偏差是各种算法所存的问题，只能够尽可能地减小这一偏差而不能消除，因此，应合理选择麦克风间距及相邻环间隔等参数。

而且，在进行插值时，测量所带来的测量误差也会影响声压待建点的数据，因此，为了减小测量误差，可以多次测量，选取合理数值。

4 结论

为解决声全息技术受测量空间的限制而无法得到无限大平面上准确声压问题，这里提出了牛顿插值与克希荷夫积分公式相结合的声压获取方法。该方法通过对声源声矢量以及声压进行插值，将离散信息点转化为连续点，然后重采样出大量离散信息点构建克希荷夫面的声压，再通过克希荷夫积分公式算得整个全息面上的声压。该方法可以在不增加实际测量点的情况下通过部分测量点的信息得到全息面上任意点处的声压信息，接着利用仿真对影响参数进行了选取，最后通过实测实验对该理论的有效性及准确性进行验证，证实该方法能够通过理论计算的方式得到无限大面上的声压且误差在允许的范围之类，且计算较为简单，节省了测量及计算时间，可以有效提高重建的效率。