响度三维空间分布重建

金江明， 卢奂采， 李敏宗

(浙江工业大学 机械工程学院，杭州 310014)

近场声全息(Near-field Acoustic Holography, NAH)方法通过傅里叶变换[1-2]和球面波函数逼近[3-4]，可以实现阵列有效测量频带范围内声场声压、粒子速度[5]、声强[6]的三维分布重建[7]，由声场逆运算也可得以到最靠近声源位置的声场信息。球形传声器阵列[8-9]具有经纬角的全指向性和三维对称性，针对声场情况优化设计的球形传声器阵列结合球面近场声全息方法(Spherical Near-field Acoustic Holography, SNAH)，适用于汽车轿厢、高铁车厢等封闭空间声场的声压三维空间分布重构和声源识别定位。

但声场降噪的最终目的是提高人耳的听觉舒适性。因此在近场声全息声学量空间分布重建和噪声源定位的基础上，有必要进一步开展声品质评价分析，重建声品质客观参量的三维空间分布，得到恼人的(而不是声压级高的)声源的位置等信息，才能更有效地实施减振降噪。

目前，Song等[10-11]利用球谐波束形成方法与Eberhard等[12-13]响度模型，开展了响度二维空间分布重建方法的研究。Song等[14]的研究结果表明高声压级噪声源所在的位置不一定响度值就高。波束形成算法[15]计算结果为阵列位置处声场的指向性信息，为获得用于响度计算的声压，Hald[16]利用一个缩放因子，与波束形成算法的二维输出结果相乘得到一个声压贡献值来近似表示声压，但该声压贡献并不是符合严格物理定义的真实声压。

本文开展了基于球面近场声全息计算结果的声场响度三维空间分布的数值仿真和实验验证。第1节推导了响度三维空间分布重建的模型；第2节对单声源声场和双声源声场进行了数值仿真计算，讨论响度三维空间分布重建方法适用的重构距离和频率范围，并分析双声源间夹角变化对声压和响度重建精度以及声源识别效果的影响；第3节给出全消声室内，使用刚性表面球形传声器阵列为测量前端重建两个扬声器声场的声压和响度三维分布的实验结果，检验了实际声源声场中该方法的重建效果；最后给出本文的结论。

1 理论模型

本文在球面近场声全息模型和ANSI_3.4—2007标准[17]给出的响度计算模型的基础上，建立声压分布到响度分布的矩阵映射模型，获得响度三维空间分布。

1.1 声压场分布重构的数学模型

刚性表面球形传声器阵列结合球面近场声全息方法可重构出处于近场空间的声压场分布，图1表示了阵列与声源的位置关系和球面近场声全息方法的有效区域。

均匀无黏性的声传播介质中空间任意一点处的波动方程，经傅里叶变换可转化为Helmholtz方程，方程解可由变量分离及系列球谐函数展开获得，具体可表示为

(1)

(2)

式中： dΩ=sinθdθdφ；hn(·)为第二类球汉克尔函数；符号“’”表示求导；符号“*”表示复数共轭。

图1 声场模型示意图Fig.1 Sound field diagram

(3)

式中：M为传声器数目；p(a,θi,φi)为第i个传声器测量到的声压；wi为相应的权重。

(4)

式中，prec(r,θ,φ,ω)为位置(r,θ,φ)处的重构声压，式(2)中的球谐函数扩展项为无穷多项，由于数值积分的原因被截断为N项，最优截止项数N的选取取决于频率和传声器数量，对于36个传声器球阵列，N=4时声场重构计算精度最高，文中SNAH算法的扩展项数均为4项。

1.2 响度数学模型

这里采用ANSI_3.4—2007标准定义的Moore响度模型计算声场指定位置的响度。响度计算步骤主要包括原始信号的频率分段处理、计算外耳和中耳传递函数、计算耳蜗听觉滤波器响应以及计算总响度。具体计算步骤如下：

步骤1原始信号的处理

若原始信号为多频率成分的复合音或者频率单一的纯音，以频率成分fj以及相应的声压级作为输入计算响度。若原始信号为多个频带组成的混合噪声，声音就等效为一系列离散正弦信号(频带中心频率)的组合。

步骤2外耳和中耳传递函数计算

声音在从外耳传播到内耳传输过程中，耳道等身体结构等会对声音的传递产生增强或者衰减的影响，分别用外耳和中耳传递函数表示这一过程，函数计算值可由ANSI标准给出的一系列离散值插值获得。频率fj声压的声压级加上外耳传递函数和中耳传递函数的计算值就得到有效频率fj及内耳耳蜗处的声压级Lpj， 写成以Pa为单位声压Pj为：Pj=10Lpj/20×20×10-6。

步骤3计算内耳耳蜗处的ERB(Equivalent Rectangular Bandwidth)级

ERB级(LERB)定义为在耳蜗处以fj中心频率的ERB带宽内的所有声音的能量和，其计算公式为

(5)

式中：P0=2.5×10-5为参考声压；ERB为一个ERB带宽，ERB=24.673(0.004 368fj+1)；fj±ERB/2为以第j个有效频率成分为中心的ERB带宽范围；fk为fj-ERB/2～fj+ERB/2以1 Hz频率间隔划分的频率；W(gj)为第j个有效频率成分fj对应的权函数，计算公式为W(gj)=(1+pjgj)e(-pjgj)， 其中gj=|fk-fj|/fj，pj=4fj/ERB。

步骤4计算单个滤波器的输出激励

为了模拟耳蜗对声音信号的滤波处理过程，ANSI标准将该听觉范围(50～15 000 Hz以划分间隔分成372个频带(ERB带宽)。每个听觉滤波器的中心频率由ERBnumber(即1.8，1.9,…,38.7，38.9)确定，具体计算公式为

听觉滤波器的输出激励为

(6)

式中：Ei为第i个听觉滤波器的输出激励；Pj为第j个有效频率处的声压值；E0为参考输出激励值，表示声压级为0 dB的1 000 Hz纯音在中心频率1 000 Hz的听觉滤波器的输出激励值；W(gij)为第i个滤波器在第j个有效频率处的响应值，计算公式为W(gij)=(1+pijgij)e(-pijgij)，其中

(7)

(8)

步骤5计算特征响度，对其求和获得总响度

特征响度计算公式为

(9)

在分别求得372个听觉滤波器的特征响度后，将其求和并除以10即可得到总响度，计算公式为

(10)

1.3 响度三维空间分布重建的计算流程

基于球面近场声全息方法重建响度三维空间分布的计算流程如下：

首先将刚性表面球形传声器阵列上测得的时域信号经过FFT(Fast Fourier Transform)变换转换为频域的声压信号pt，然后依据SNAH理论重构出声场中入射声压pi，并将其作为响度计算模型的输入，与响度映射矩阵相乘，得到响度三维空间分布计算结果。声场中某点响度映射矩阵可写为

(11)

或简写为：PiW=N′

式中：Pi为声场中某位置处各频率下的声压重构值所组成的矢量；W为由372个滤波器组成的听觉滤波器矩阵，表示人耳对可听频带内所有频率的响应；N′为特征响度矢量。由式(11)对特征响度矢量N′中的各项求和便可获得空间特定点的响度，而对声场中所有位置重复这一计算过程便可获得声场响度的三维空间分布。具体计算步骤总结如下：

步骤1球形阵列采集的声场时域信号；

步骤2声场时域信号的FFT变换；

步骤3SNAH理论重构声场声压；

步骤4响度理论模型，响度矩阵映射模型，获得响度三维空间分布计算；

步骤5确定人耳听觉感受最响的噪声源位置。

2 仿真验证

对封闭空间声场中由球面近场声全息方法得到的声压场的三维空间分布，应用Moore响度模型，可以进一步计算得到声场中不同位置的响度值，由此实现了响度三维分布的重建。本节将在单声源辐射形成的声场和两个声源辐射形成的声场中开展响度三维分布重建的数值仿真验证。

单声源声场的仿真结果给出声场响度三维空间分布重建方法中重构距离、频率等参数的适用范围。双点声源声场的仿真计算结果给出了不同声压级、不同频率成分组成对响度三维空间分布计算的重构误差影响。

用于声压分布计算结果和响度分布计算结果评价的二范数误差公式由如下公式定义。声压重构的误差Ep为

(12)

式中：J为声场重构的空间点数；pij为第j个计算点处的入射声压重构值，由式(2)获得;pj为第j个计算点处的根据点声源获得的入射声压理论值。

响度重构误差定义为

(13)

式中：Nij为第j个计算点处的响度重构值；Nj为第j个计算点处的响度理论值， 分别由声压重构值pij和理论值pj计算获得。

2.1 单声源声场仿真及计算结果分析

如图2所示，单极子点声源布置在x轴正方向上1 m的位置，而半径a=0.1 m的球形传声器阵列放置于坐标原点。球形传声器阵列共有36个传声器，并均布于传声器阵列的刚性球表面上。

图2 单声源声场模型Fig.2 Sound field diagram with one monopole source

设定点声源强度Qs=3.6×10-5，频率fi在100 Hz～1 300 Hz内，计算半径r从0.1 m～0.9 m球面上的声压、响度的重构误差。

图3给出了球面上声压、响度的重构误差与频率、重构半径r之间的关系曲线。从图中可以发现，声压、响度的重构值误差均随重构距离增加和声源辐射频率升高而变大。图3(a)表明在重构距离小于0.3 m且频率不大于1 200 Hz时，声压的重构误差小于20%。但是当频率等于1 300 Hz时，声压重构误差快速增加后几乎不随重构距离增加而变化，说明此频率下响度三维空间分布重建已不适用于现有球形传声器阵列。图4(b)中给出不同频率声音的响度的计算误差随重构距离的变化曲线，从图中可以发现，在频率小于700 Hz且重构距离小于0.3 m时，响度的重构误差小于20%，但随着频率增加误差将快速增加。对比分析图3(a)和图3(b)，可以发现，响度重构误差的分布趋势与声压的重构误差分布相似，但数值大于声压的重构误差，这是因为响度误差是有效频带内所有频率的计算结果综合后一个误差。这也表明近场声全息声场重构误差经中耳和外耳传递函数计算、耳蜗滤波处理后被成倍放大，因此用于响度三维空间分布计算的声压场要有较高的重建精度，由此才能得到理想的响度三维空间分布的计算结果。

(a) 声压重构误差随重构半径的变化图

(b) 响度重构误差随重构半径的变化图图3 半径为r的球面上声压、响度重构误差与频率、重构半径r之间的关系Fig.3 The reconstruction errors of the sound pressure, loudness varying with different calculation radius, frequency

这里需要特别说明的是，鉴于基于球形传声器阵列的球面近场声全息方法的有效频带为100 Hz～1 300 Hz，而响度又是表示20 Hz～20 kHz频段内所有频率声音的响亮程度，在目前的声场重建中本方法仅仅计算有效频带范围内所有声音的响度，对于有效频带范围以外频率的声音，由于不能有效实现声场重建，因此本文计算的响度仅包含有效频带的内声音，对有效频带范围以外的声音未予以考虑。对于这些频率的声音，可以采用不同直径的刚性表面球形传声器阵列实现声压场重建，Williams等的研究结果表明,如果采用适当尺寸的传声器阵列，近场声全息方法可以实现最高15 000 Hz频率的声场重建。

图4给出了图3声场中fs=600 Hz，d=1 m时，半径r=0.4 m球面上的声压、响度的重构值与理论值在重构球面上展开的空间分布对比图，为更清楚显示误差在球面上分布情况，这里选择计算一个较大半径(r=0.4 m)球面上声压和响度的分布。从图4中可以看到，声压、响度的重构值实现了点声源的位置定位，相比于理论值的，其声场声压重构计值算误差为误差小于6%, 响度重构计值算误差为小于22%。

图5给出了声压、响度的重构值与理论值在重构球面上计算点处的误差。图4、图5给出的仿真计算结果表明声场响度三维空间分布重建方法有效重构了频率为600 Hz的点声源声场在半径r=0.4 m、中心为原点的球面上声压、响度分布。

图4 声压、响度的重构值与理论值的球面展开图Fig.4 Spherical surface expanding for reconstruct and analytical value of space distribution of the sound pressure and loudness (fs=600 Hz，r=0.4 m)

(a) 声压重构误差分布

2.2 双声源声场仿真及计算结果分析

图6给出了双声源声场的仿真分析的声场模型。图中，声源1固定在x轴正方向d1处，而声源2与声源1之间成θ角并位于x-y平面内，声源2到坐标原点的距离为d2，采用与单声源声场相同的球形传声器阵列，同样放置于坐标原点。

图6 双声源声场模型Fig.6 Sound field diagram with two monopole sources

取声场两个点声源强度分别为Qs1=2.5×10-5和Qs2=2.0×10-5，d1=d2=1 m， 频率f1、f2采用表1所示的频率组合，利用声场响度三维空间分布重建方法，计算半径r从0.1 m～0.9 m球面上的声压、响度的重构误差。

表1 双声源的不同频率组合

图7表示了两声源间夹角θ=60°时，声压、响度的重构误差随重构距离、频率组合的变化曲线。图7的误差曲线变化规律与单声源声场中的重构误差曲线类似，声压、响度的重构误差随着重构距离的增大而增大，在相同频率的低频声源情况下，重构误差也随着双声源声场中声源2频率升高而增大。从图中也可以看出当重构距离大于0.3 m时，声压、响度的误差显著增大，为此，文中响度三维空间分布的计算限于离原点距离小于0.3 m的声场空间内。另外，双声源声场中声源间夹角变化对重构误差的影响较小，说明不同的声源空间分布对误差最终计算结果影响不大。

(a) 声压重构误差随重构半径的变化图

(b) 响度重构误差随重构半径的变化图图7 两声源间夹角θ=60°时，声压、响度重构误差与频率、重构距离间的关系曲线Fig.7 Reconstruction error of sound pressure and loudness varying with frequency and reconstruction distance when θ=60°

在图6所示的双声源声场的仿真分析模型中，保持d1=d2=0.3 m，f1=150 Hz，f2=800 Hz，声源强度为Qs1=2.5×10-5、Qs2=2.0×10-5(Qs1>Qs2)不变， 让θ在0°～180°以30°为间隔内变化，以此考察声源间夹角θ变化对响度三维空间分布计算结果的影响。

图8给出了重构距离分别为0.1 m和0.3 m时的声压、响度的空间分布图。图中声源1位于每幅图的中心位置，声源2则随夹角θ增大沿水平中心线往左边移动，在θ=180°时声源2位于图的两边边线。从图中的对比可以发现，声压最大的空间位置与响度最大的空间位置并不一定重合，在声源夹角θ=180°时，响度定位结果与声压级声源定位结果完全相反。此计算结果说明：依据声压分布图识别的噪声源并不一定是人耳听到最响的噪声源，而是应该根据响度分布才能有效识别噪声源。

另一方面，对比分析不同声源间夹角大小时的声场分布计算结果发现，响度分布计算结果比声压分布计算结果能识别出更小的声源夹角时两个声源所在的位置，而联合分析响度和声压分布计算结果，能得到比仅根据单一声压、响度分布更好得声源识别定位结果，响度三维空间重建方法一定程度上提高了声源识别定位精度。同时，对比分析重构半径0.1 m和0.3 m的响度及声压三维空间分布计算结果，并综合前面声场重构误差计算结果，可以发现，虽然重构距离增加使声压、响度重构误差有一定的增加，但重构距离增加不改变声场整体变化趋势，也没有改变声源的最终定位结果。

图8 两声源间夹角θ变化时，声压、响度重构值空间分布的球面展开图 (f1=150 Hz，f2=800 Hz)Fig.8 Spherical surface expanding for reconstruct and analytical value of space distribution of the sound pressure and loudness varying with θ

3 实验验证

仿真计算都是在理想点声源、无干扰条件下进行的。为进一步检验声场响度三维空间分布重建方法的计算精度，在全消声室内开展了使用刚性表面球形传声器阵列测量两个扬声器辐射形成的声场的实验研究。

3.1 实验设备和实验过程

实验设备在全消声室内的布置如图9所示。球形传声器阵列位于两个扬声器的中间，三者成一条直线(即此时扬声器相对于球形传声器阵列中心成180°布置)，两个扬声器距离阵列中心都为0.2 m，球形传声器阵列采用B&K Type4958 36通道刚性表面球形传声器阵列，阵列直径为0.2 m，扬声器纸盆直径约为7 cm，全消声室的本底噪声为18 dB，最低频率为63 Hz，自由声场大小为1.5 m×1.2 m×1.5 m。实验采用LAN-Ⅺ 36通道数采和PULSE系统来记录时域数据。

图9 实验设备在全消声室内布置图Fig.9 Equipment layout in anechoic chamber

实验进行时，PULSE LAN-Ⅺ系统输出两个单频正弦信号，由B&K 2716功放驱动声源1和声源2依次产生如表2、表3所示的不同频率组合的单频声音，声音信号频率都在刚性表面球阵列的有效频带范围内。PULSE LAN-Ⅺ系统单次测量的数据记录时间为10 s。采集的时域数据观察无异常后，由PLUSE系统进行FFT变换得到单个每个通道的频谱数据，并提供给球面近场声全息算法做声场重构算法。

3.2 实验结果分析

图10、图11给出了根据实验测量数据计算获得的声压和响度的三维空间分布计算结果，其声场重构计算半径为0.2 m，球形传声器阵列中心为计算坐标系的中心。图中结果表明在声压分布图和响度分布图中最大峰值所在位置不同，主要声源位置的定位结果不同，从而表明声压级高的声源并不一定是人耳最容易听见的声源。

图10 声压、响度重构值的空间分布图Fig.10 Sound pressure varying with frequency, reconstruction distance

图11 声压、响度重构值的空间分布图Fig.11 Space distribution of the reconstruction valve of the sound pressure and loudness

表2 在不同声源强度组合下声压、响度的空间分布图中声源1、声源2位置对应的声压和响度的峰值Tab.2 Peak value in the 3D space distribution of the sound pressure and loudness in different intensity level of the source 1, source 2

表2给出了两个扬声器采用同一频率组合但在不同声源强度组合下声源1、声源2所在位置对应的声压和响度的峰值大小，表中声源声压级大小由B&K Type4190 1/2英寸自由场传声器在扬声器表面测得。表5给出结果表明随着较高频率声源2的声压级从70 dB增加74 dB，声源2附近的响度值从5.3 sones增加到6.6 sones，声源2在73 dB时已大于较低频率声源1(400 Hz,74 dB)所在位置的响度值。由于响度值是所有频率成分综合后一个结果，因此声源2辐射声压级增加也使得声源1位置处响度值同步增加，但增加的幅度较小。

表3 在不同频率成分组合下声压、响度重构值的空间分布图中声源1、声源2位置对应的声压和响度的峰值Tab.3 Peak value in the space distribution of the sound pressure and loudness in case of the different frequency combination of the source 1, source 2

表3给出了声源强度相同但频率成分组合不同时声源1、声源2对应的声压和响度的峰值大小，随着声源2的频率变高，声场响度值也同步增加。以上两个表的计算结果揭示频率高低和声压级大小对声学成像图中声压和响度峰值位置的影响。

4 结 论

声场响度三维空间分布重建方法在给出声场声压的三维空间分布的同时，还能给出响度的三维空间分布，实现了识别定位与人耳听觉感受最响的噪声源位置的目的。通过研究，得到如下结论：

(1) 声场响度三维空间分布重建计算结果表明在一些声场中，声压峰值所在空间位置与响度峰值所在空间位置不同。

(2) 计算结果给出了响度值的变化与声源空间位置分布、频率间的相互关系。

(3) 重构计算结果误差随着频率升高和重构半径增加而变大，但两声源间的夹角变化对重构误差的影响较小。

(4) 计算结果表明响度计算放大了声场重构误差，因此 保证SNAH有较高的声场重构精度，才能有较理想的响度三维空间分布重构结果。同时，声场分布的趋势性没有发生大的变化。

(5) 将声压、响度两种空间分布计算结果联合应用于声场的声源识别定位中，提高了多声源声场的声源识别定位精度。

[ 1 ] VERONESI W A, MAYNARD J D. Nearfield acoustic holography (NAH) II. Holographic reconstruction algorithms and computer implementation[J]. The Journal of the Acoustical Society of America, 1987, 81(5): 1307-1322.

[ 2 ] MAYNARD J D, WILLIAMS E G, LEE Y. Nearfield acoustic holography: I. Theory of generalized holography and the development of NAH[J]. The Journal of the Acoustical Society of America, 1985, 78(4): 1395-1413.

[ 3 ] LU H, WU S F. Reconstruction of vibroacoustic responses of a highly nonspherical structure using Helmholtz equation least-squares method[J]. The Journal of the Acoustical Society of America, 2009, 125(3): 1538-1548.

[ 4 ] WANG Z, WU S F. Helmholtz equation-least-squares method for reconstructing the acoustic pressure field[J]. The Journal of the Acoustical Society of America, 1997, 102(4): 2020-2032.

[ 5 ] WILLIAMS E G, HOUSTON B H, HERDIC P C, et al. Interior near-field acoustical holography in flight[J]. Journal of the Acoustical Society of America, 2000, 108(4): 1451-1463.

[ 6 ] WILLIAMS E G, VALDIVIA N, HERDIC P C, et al. Volumetric acoustic vector intensity imager[J]. Journal of the Acoustical Society of America, 2006, 120(4): 1887-1897.

[ 7 ] LU H, LI M. Reconstruction of sound field based on near-field acoustic holography with a rigid spherical microphone array[C]//ACOUSTIC 2012 Hong Kong Conference and Exhibition. HongKong: ACOUSTIC, 2012: 1-11.

[ 8 ] JACOBSEN F, MORENO-PESCADOR G, FERNANDEZ-GRANDE E, et al. Near field acoustic holography with microphones on a rigid sphere[J]. The Journal of the Acoustical Society of America, 2011, 129(6): 3461-3464.

[ 9 ] WILLIAMS E G, TAKASHIMA K. Vector intensity reconstructions in a volume surrounding a rigid spherical microphone array[J]. The Journal of the Acoustical Society of America, 2010, 127(2): 773-783.

[10] SONG W, ELLERMEIER W, HALD J. Using beamforming and binaural synthesis for the psychoacoustical evaluation of target sources in noise[J]. The Journal of the Acoustical Society of America, 2008, 123(2): 910-924.

[11] SONG W, SAITO H, HADDAD K. Improved noise source identification using sound quality metrics mapping in vehicle noise measurements[Z]. 2011.

[12] EBERHARD Z, HUGO F, ULRICH W, et al. Program for calculating loudness according to DIN 45631 (ISO 532B)[J]. Journal of the Acoustical Society of Japan (E), 1991, 12(1): 39-42.

[13] FASTL H, ZWICKER E. Psychoacoustics: facts and models[M]. Berlin: Springer, 2006.

[14] SONG W, ELLERMEIER W, MINNAAR P. Loudness estimation of simultaneous sources using beamforming[Z]. 2006.

[15] TREES H L V. Detection, estimation and modulation theory. Part IV, optimum array processing[M]. New York: John Wiley&Sons, 2002.

[16] HALD J. Estimation of partial area sound power data with beamforming: inter-noise[Z]. 2005.

[17] The acoustical society of america. Procedure for the computation of loudness of steady sounds：ANSI_3.4—2007 [S]. [S.l.]: The American National Standards Institute, 2005.