球形正十二面体声源测量近场HRTF的误差*

余光正 陈泽为

(华南理工大学声学研究所，广东广州510640)

自由场条件下，点声源产生的声波经头部、躯干和耳廓等人体生理结构散射后到达双耳，构成一个线性时不变的声学传输系统，并可描述为头相关传输函数(HRTF).当声源与头中心的距离不小于1.0m时，HRTF近似与声源距离无关，称为远场HRTF;当距离小于1.0 m时，HRTF随距离变化明显，称为近场HRTF.HRTF包含声源的空间定位信息，因而对双耳听觉的研究具有重要意义，并在虚拟听觉重放中有重要的应用［1-2］.

实验测量是准确获取HRTF的主要方法.但在近场HRTF测量中，由于声源靠近受试者，声源尺度引起的多重散射和指向性误差将比远场测量时更明显.为减小声源引起的误差，一些研究中已开始采用几种小尺度声源［3-5］.这些小尺度声源虽能减少多重散射和改善无指向特性，但在1 kHz以下的低频特性变差，不能准确测量低频的近场 HRTF，而HRTF的低频特性对声源距离定位却十分重要［6］.

对于一般的扬声器系统，改善低频特性与减少声源尺度经常是矛盾的.笔者采用脉动球形声源和球形头部的组合模型评估了声源多重散射引起的HRTF测量误差，发现当以1.0dB的幅度误差为限，声源距离不小于0.20 m，频率上限为20 kHz时，球形声源等效半径应不大于0.05 m［7］.但理想的脉动声源(球面波声源)并不存在，实际的声源在高频具有一定指向性.为近似实现无指向的球面波声源，笔者设计制作了半径为0.035m的球形正十二面体声源，其工作频率范围为0.35～20.00 kHz，在低于8kHz的频段近似无指向(±3dB)，较已有声源有较明显的改善［8］.

在近场HRTF测量中，声源多重散射和指向性引起的误差是共存的.因此，文中提出球形正十二面体声源和球形头部组合模型，用球谐函数多极展开和球基函数旋转变换方法，计算和分析声源多重散射和指向性对近场HRTF测量的共同影响，并提出了改进办法，最后用自制球形头部模型和球形正十二面体声源的测量结果对理论计算结果进行了验证.

1 HRTF的计算模型及求解

1.1 刚球模型HRTF的定义

文中所采用的坐标系统如图1所示，声源中心位置用球坐标(r0，α0，β0)或空间位置矢量r0表示.其中声源距离r0通常应大于0.2m;仰角α0的取值范围为0°～180°;方位角β0的取值范围为－90°～270°.α0=90°表示水平面;水平面上，β0= － 90°和90°分别表示正左方和正右方.在该坐标系，将图1的球形正十二面体声源替换为点声源，即可得到刚球模型［6］精确的 HRTF，简记为H*.考虑到模型的对称性，左、右耳受声点(即EL和ER)的HRTF统一定义为

式中:f表示频率，P*表示点声源在左、右耳受声点处的声压，P*0表示点声源在刚球头部中心的自由场辐射声压.刚球模型HRTF的计算方法可参考文献［6］.

1.2 计算模型及误差公式

球形头部(文中将半径固定为0.0875 m)和球形正十二面体声源组成的计算模型如图1所示.

图1 刚球头部模型和球形正十二面体声源的组合模型Fig.1 Combined model consisting of a rigid spherical head model and a spherical dodecahedral sound source

图1中，E表示空间任意受声点;在以O和O'为原点的球坐标系中，分别用r和r'表示E点的位置矢量，用(r，α，β)和(r'，α'，β')表示E点的球坐标; 而(x，y，z)和(x'，y'，z') 分别为对应的笛卡尔坐标.对于图1的计算模型，其HRTF(记为H)可参考式(1)定义为

式中，P表示球形正十二面体声源在左、右耳受声点的声压，P0为球形正十二面体声源在球形头部中心的自由场辐射声压.

文中主要分析和讨论由球形正十二面体声源的散射和指向性共同引起的近场HRTF测量误差，因此这里先给出HRTF测量的幅度谱误差定义:

式中:SD表示幅度谱误差，单位为dB.要求解式(3)的SD，须先求解式(2)的组合模型HRTF，即求解式(2)的双耳声压P和自由场声压P0.P0的求解可参考文献［8］，下面简要推导双耳声压P的求解公式.

1.3 双耳声压P的求解

在图1中，任意受声点E的空间位置矢量用r表示(左、右耳的空间位置矢量rL和rR是其特例).受声点E的复数声压P(r，r0，f)应满足 Helmholtz方程:

式中:k=2πf/c，表示波数;c表示空气中的声速.球形头部表面SH(设为刚性声学边界［9］)和球形正十二面体表面SD(局部脉动单元以外的区域设为刚性声学边界)的声压分别满足以下边界条件:

式中:∂/∂n表示法向导数;j为单位复数;ω=2πf，为角频率;ρ0为空气密度;v为球形正十二面体声源表面的法向速度，且有

式中:i表示局部脉动单元的编号，i=1，2，…，12;v0表示局部脉动单元法向振幅;θ0表示局部脉动单元的半张角;αi为图2 所示球坐标系(ri，αi，βi)中的仰角.在无穷远，声压P(r，r0，f)满足 Sommerfeld辐射条件:

式中，∂/∂r表示对距离求偏导数.求解声压P(r，r0，f)也就是求解式(4)在边界条件(5)－(8)下的解.一般情况下，该声压可分解为3部分的和:

式中:P0(r，r0，f)表示声源在受声点E的直达声压;PH(r，r0，f)表示球形头部的散射声压;PD(r，r0，f)表示球形正十二面体声源的散射声压.三者可按球谐函数展开为表达式(11) －(12)［10］.其中，式(11)在球坐标系(r，α，β)中展开，式(10)和(12)在球坐标系(r'，α'，β')中展开(见图 1).

式中，hl是l阶第二类球汉克尔函数，Ylm是复球谐函数.求解声压P(r，r0，f)转而成为求解式(10)－(12)的展开系数Dlm、Alm、Blm.实际计算时，先解出系数Dlm，Alm和Blm则可通过边界条件式(5)、(6)联立方程组求解［7］.

求解式(10)中的展开系数Dlm时，分别将12个局部脉动单元的辐射声压按式(10)展开并线性叠加.当第i个局部圆形脉动单元的中心不经过z'轴时，须引入图 2 所示的球坐标(ri，αi，βi)，它对应的笛卡尔坐标为(xi，yi，zi)，此时，第i个局部圆形脉动单元的辐射声压PDi(ri，f)可表示为［11］

式中，展开系数Dil的求解公式为

aD为球形正十二面体声源半径，h'l是l阶第二类球汉克尔函数的导数，系数Cl满足

Γl是l阶的勒让德多项式，且Γ－1=1.

为了把球坐标(ri，αi，βi)下的表达式(13)用球坐标(r'，α'，β')重新表达，需用到以下的球基函数旋转变换［12］(示意图见图2):

式中，α'i和β'i表示第i个局部脉动单元中心在球坐标系(r'，α'，β')中的方位角和仰角，变换系数T可用下面的公式求解［12］:

图2 球基函数旋转变换Fig.2 Rotation transform of the spherical basis functions

利用公式(13)－(17)，可将任意第i个局部脉动单元的辐射声压表示在球坐标系(r'，α'，β')中，然后对12个局部脉动单元的辐射声压求和，便可得到式(10)的展开系数Dlm:

将Dlm代入式(10)－(12)，可结合边界条件求出各项展开系数，从而求解出双耳声压和HRTF.

2 HRTF的误差分析

2.1 误差的算例

为简单验证计算结果的准确性，先对幅度谱误差SD分步进行求解，表达式如下:

式中，SD1和SD2按下式确定:

H1是忽略式(9)的PD(r，r0，f)项后算得的 HRTF，因而第1步求解得到的SD1只反映声源多重散射的影响，第2步求解的SD2只反映声源指向性的影响.作为算例，球形正十二面体声源半径aD取为0.035m，脉动单元半张角θ0为23.6°，声源位于r0=0.2 m、α0=90°、β0=90°处，算得右耳受声点ER的 HRTF幅度谱误差 SD1、SD2、SD，如图3 所示.

由图3可知:SD1在±0.5 dB的范围内，SD2在±2.5dB的范围;SD恰好是SD1和SD2的叠加.且SD1和SD2的结果与文献［7］和［8］的计算结果基本一致，说明计算结果是准确的.

图3 右耳HRTF的 SD1、SD2和SDFig.3 SD1，SD2and SD of HRTF of right ear

2.2 总体误差与声源改进

随着声源靠近球体，声源散射误差将变得明显［7］，因而这里只取0.2 m(一般作为分析近场HRTF的最小声源距离［2］)进行分析.由图1所示计算模型的对称性，只给出α0=90°(水平面)、声源方位角β0在－90°～90°(前半周)时的右耳 HRTF幅度谱误差，见图4.

aD=0.035m和θ0=23.6°时的计算结果如图4(a)所示.可见，明显的谱误差出现在约8 kHz以上;随着声源偏离受声耳(β0偏离90°)，HRTF的幅度谱误差略有增大，但出现误差的频段基本保持一致.

为了改进现有的声源，另计算分析了两种情况下的声源误差.图4(b)是减小声源尺寸至aD=0.025m且保持脉动单元半张角不变时的计算结果;图4(c)给出了脉动单元的半张角增大到θ0=30°(接近理想的最大值31.7°)时的计算结果.由图中可见，随着声源尺寸减小和局部脉动单元半张角增大，球形正十二面体声源多重散射和指向性引起的误差幅度减小，且出现误差的频率升高，这对近场HRTF的测量是有利的.如图4(b)所示，当声源半径为0.025 m、局部脉动单元半张角为23.6°时，较明显的误差(一般指大于1 dB)出现在约12 kHz以上，除个别方向( －90°≤β0≤ －50°)和频率(f≥10kHz)，误差幅度可控制在 ±3.5 dB.事实上，在10kHz以上的频率，HRTF所包含的谱特性对定位的贡献已逐渐减少，而各种因素(与声源无关)导致测量HRTF的幅度谱误差较大(通常在10 dB左右).因此，可认为图4(b)的声源引起的误差已基本满足近场HRTF的测量要求，为进一步改进声源提供了理论依据.

图4 r0=0.2m，α0=90°时右耳HRTF的SDFig.4 SD of HRTF of right ear when r0=0.2m and α0=90°

3 测量验证

为进一步验证上述计算结果的合理性，使用实际的球形正十二面体声源(半径aD约0.035m，脉动单元半张角取θ0=23.6°)和木质刚球模型(半径约0.0875 m)进行测量验证，如图5所示(声源位于r0=0.2m、α0=90°、β0=90°处).HRTF 的测量环境和方法与文献［8］相同.另一方面，按照图5实际模型的尺寸，计算图1所示组合模型条件下的HRTF.测量和计算的HRTF结果见图6.由图6可见，在0.35～10kHz的频段，测量值与计算值基本一致，声源指向性和散射的总误差约在±2.5 dB以内.测量值在0.35Hz以下明显衰减，这是由声源的频率响应特性所决定的.其它谱细节差异是由实验误差(主要是实际制作的声源和球体与理想模型的差异)引起的.上述实验结果进一步验证了计算结果的合理性.

图5 刚球模型以及球形正十二面体声源Fig.5 Rigid spherical model and spherical dodecahedral sound source

图6 测量和计算的右耳HRTFFig.6 Measured and calculated HRTF of right ear

4 结语

文中提出了球形正十二面体声源和球形头部的组合计算模型，利用球谐函数多极展开和球基函数旋转变换的方法，计算和分析了近场HRTF测量中由球形正十二面体声源指向性和多重散射共同引起的HRTF幅度谱误差.计算结果的准确性通过两步求解的办法得到了验证，计算结果的合理性用实际制作的球形正十二面体声源和木质球体进行了实验验证.文中还在不同方向全面分析了声源对近场HRTF测量的影响，结果表明，随着声源尺寸的减小和局部脉动单元半张角的增大，球形正十二面体声源多重散射和指向性引起的误差幅度减小，且有向高频移动的趋势.当声源距离不小于0.2 m、声源半径不大于0.025 m、局部脉动单元半张角不小于23.6°时，较明显的误差出现在约12 kHz以上，且误差幅度基本在±3.5 dB，可满足近场HRTF的测量要求.文中结果为进一步改进适合于近场HRTF测量的球形正十二面体声源提供了理论依据.

［1］谢菠荪.头相关传输函数与虚拟听觉［M］.北京:国防工业出版社，2008:1-33.

［2］余光正，谢菠荪.近场头相关传输函数及其应用［J］.电声技术，2007，31(7):45-50.Yu Guang-zheng，Xie Bo-sun.Head-related transfer function for nearby sources and its application［J］.Audio Engineering，2007，31(7):45-50.

［3］Nishino T，Hosoe S，Takeda K，et al.Measurement of the head related transfer function using the spark noise［C］∥The 18th International Congress on Acoustics.Kyoto:Science Council of Japan，2004:1437-1438.

［4］Hosoe S，Nishino T，Itou K，et al.Development of micrododecahedral for measuring head-related transfer functions in the proximal region［C］∥Proceedings of the 2006 IEEE International Conference on Acoustics，Speech，and Signal Proceedings.Toulouse:Centre de Congrès Pierre Baudis，2006:329-332.

［5］Brungart D S，Rabinowitz W M.Auditory localization of nearby sources.head-related transfer functions［J］.Journal of the Acoustic Society of America，1999，106(3):1465-1479.

［6］Duda R O，Martens W L.Range dependence of the response of a spherical head model［J］.Journal of the Acoustic Society of America，1998，104(5):3048-3058.

［7］Yu G Z，Xie B S，Rao D.Effect of sound source scattering on measurement of near-field head-related transfer function ［J］.Chinese Physics Letters，2008，25(8):2926-2929.

［8］Yu G Z，Xie B S，Rao D.Directivity analysis on spherical regular polyhedron sound source used for near-field HRTF measurements［J］.Chinese Physics Letters，2010，27(12):124302.

［9］余光正，谢菠荪，饶丹.躯干散射对近场HRTF影响的微扰法分析［J］.华南理工大学学报:自然科学版，2010，38(3):143-147.Yu Guang-zheng，Xie Bo-sun，Rao Dan.Analysis of effect of torso scattering on near-field head-related transfer functions by using perturbation method ［J］.Journal of South China University of Technology:Natural Science Edition，2010，38(3):143-147.

［10］Gumerov N A，Duraiswami R.Computation of scattering from N spheres using multipole reexpansion［J］.Journal of the Acoustic Society of America，2002，112(6):2688-2701.

［11］张海澜.理论声学［M］.北京:高等教育出版社，2007:244-248.

［12］Gumerov N A，Duraiswami R.Fast multipole methods for the Helmholtz equation in three dimensions［M］.Oxford:Elsevier，2004:122-124.