郭 莹， 闫美辰

(沈阳工业大学 信息科学与工程学院， 沈阳 110870)

基于镜像源方法的室内声场脉冲响应仿真*

郭 莹， 闫美辰

(沈阳工业大学 信息科学与工程学院， 沈阳 110870)

针对基于室内声场脉冲响应描述声场音质特性时其精确度不高的问题，以经典镜像源方法为核心，讨论了基于变换频域和反射系数的改进算法.以矩形空间为具体实例，给定接收点和声源位置，并利用计算机仿真室内声场脉冲响应.在不同参数条件下，对各个算法进行了分析，估计了能量衰减曲线和房间混响时间.结果表明，基于变换频域和反射系数的改进算法具有更高的准确性和可靠性.

室内声场； 脉冲响应； 音质特性； 镜像源方法； 变换频域； 反射系数； 能量衰减曲线； 混响时间

由于封闭声场内声能密度、声压等参数都满足齐次性和叠加性原理，因此，可将封闭的待测室内声学空间看作一个线性时不变系统[1].时域的通道特性可用声源至接收点的声场脉冲响应来描述[2].借助现代信号处理技术，通过对脉冲响应进行测量和分析，可以得到声学空间的各种声学特性参数[3-5]，如混响时间、明晰度和声场强度等，从而准确地描述室内声场的变化规律.

当仿真计算声场脉冲响应时，需要选取合适的声场模型.镜像源方法是以几何声学模型为理论基础的经典算法，可以较准确地模拟声场脉冲响应.自1979年以来，凭借理论简单、易于理解和实现的优点，镜像源方法已被视为声学研究的基础理论，因而在声学、工程学等很多领域中得到了广泛应用[6].可以利用镜像源方法对封闭环境中声音的传播、封闭空间中噪声的控制，以及声场中的语言清晰度和语音传输指数进行分析[7]，且镜像源方法也可应用于双耳可视化和交互系统等设备的研发领域[8].

本文以经典镜像源方法为基础，结合统计声学，讨论了几种声场脉冲响应的具体算法，编制了仿真程序，给出了矩形空间声场的计算实例，并对直接影响模拟结果的几个重要方面进行了讨论.基于仿真得到了室内声场脉冲响应，并估计了能量衰减曲线和室内混响时间.

1 镜像源方法

1.1 基本原理

图1为镜像源方法俯视图.由图1可见，在封闭空间中，由声源S发射出来的声波在空气中传播时，一部分声波被接收点Q直接接收，另一部分声波将被墙面或障碍物反射，经过一定延时后再被Q接收.由于反射角等于入射角，故反射波的反向延长线与声源到墙面的垂线之间必然存在某一交点Si(i=1，2，3，4)，反射声波可以看做是从墙后的另一声源Si发射出来的，因此，Si被称为声源S的镜像源，且镜像源Si与声源S到墙壁的距离相等.当获知全部镜像源的位置与能量(或声压)后，可以把声源对于接收点的能量等效为相应镜像源的能量之和.

图1 镜像源方法俯视图Fig.1 Vertical view of image source method

镜像源数目与反射级数(声源反射次数)和壁面数目均有关系.随着镜像源数目的增加，总镜像源数目将按照指数规律递增.但是并非所有镜像源都会对封闭声场内的接收点产生能量.事实上，产生能量的镜像源只占很小的比例.因此，在镜像源方法中，一个很重要的步骤就是判断镜像源的可见性，即要从全部的镜像源中，确定哪一些镜像源可以对接收位置产生能量.这一过程需要遵循的原则为产生镜像源的壁面必须与镜像源、接收点的连线存在实际交点.

图2为镜像源的可见性判断示意图，其中实线是声音真正的传播路径，而虚线是所引辅助线.如图2所示，对于声源S和接收点Q而言，一级镜像源S′是可见的，而二级镜像源S″则是不可见的.通过镜像源的可见性判断，可以大大降低镜像源方法的计算量.

图2 镜像源的可见性判断Fig.2 Visibility judgment for image source

1.2 经典算法

可以通过一系列可见镜像源产生的声压和来模拟计算声场脉冲响应.在矩形自由场中，接收点接收到的声压(频域声压脉冲响应)可以表示为

(1)

根据镜像源原理，镜像源与接收点是关于某一墙面对称的，因此，声压还可以表示为

exp(-jωt)

(2)

式中：R-为声源到接收点的距离；R+为镜像源到接收点的距离.此外，R-与R+可以分别表示为

(3)

(4)

考虑到房间中有6面墙，且每一个一阶镜像源还会产生高阶镜像源，计算比较复杂，可将式(2)改写为

exp(-jωt)

(5)

式中：Rp为相对于坐标原点而言，声源到镜像源的位置向量；Rr为考虑高阶反射后的虚拟房间尺寸向量.此外，Rp和Rr可以分别表示为

Rp=(x±x′，y±y′，z±z′)

(6)

Rr=2(nLx，lLy，mLz)

(7)

Rp+Rr=(x±x′+2nLx，y±y′+2lLy，

z±z′+2mLz)

(8)

式中：Lx、Ly和Lz分别为房间的长、宽、高；n、l和m为相应参数.

对式(5)进行傅里叶变换后，得到的时域脉冲响应表达式为

(9)

式中，δ(·)为单位脉冲函数.

在以上推导中，均需要假设墙面是坚硬的.在实际情况中，则需要考虑到每一个壁面均存在吸声系数α和反射系数β的情况.根据Sabine公式[8]，可以得到二者的关系为

α=1-β2

(10)

当考虑到房间反射系数后，可以将式(9)改写为

(11)

式中：

Rp=(x-x′+2qx′，y-y′+2vy′，

z-z′+2kz′)

(12)

q、v和k为相应参数，且分别等于0或1.

此外，室内所有表面的声学特性可用反射系数βx，ξ、 βy，ξ和βz，ξ来表征.其中：当ξ=1时，表示紧邻坐标系原点的墙壁；ξ=2时，表示远离坐标系原点的墙壁.

2 改进算法

一般而言，由基本镜像源方法仿真得到的脉冲响应精确度较低，因而无法真实地反应室内脉冲响应.若把基本方法的室内脉冲响应表达式变换到频域进行计算，并将表达式中的反射系数定义为负反射系数，得到新的室内脉冲响应表达式，则会提高室内脉冲响应建模的精确度.

2.1 频域计算

当在时域处理离散信号时，很难把握时间点的选取.若选取离每个镜像源传播时间点最近的整数点，则时间点的选取较为粗略，因而会使脉冲响应产生相当大的误差.当在频域处理信号时，信号本身是连续的，无需再把连续信号离散化，从而避免了时间点的选取问题.由式(11)可以得到声压脉冲响应的频域表达式，即

exp(-jωt)

(13)

此时，将式(13)进行傅里叶逆变换后，可以得到时域状态下室内脉冲响应表达式，即

hf(t)=F-1(Hf(ω))

(14)

式中，F-1为傅里叶逆变换.

2.2 负反射系数

设任何房间的表面吸收系数均由吸声系数α表示，则由式(10)可以得到反射系数β的表达式[9-10]，即

(15)

由式(15)可知，β可以取正、负两个值，在传统镜像源方法中，β值取的是正数.然而，利用传统镜像源方法，得到的脉冲响应与真实测量值相比，失真情况较为严重，无法清晰地反映室内脉冲的真实情况.事实上，利用变换频域与采用负反射系数的方法，仿真得到的脉冲响应与实际情况更为接近.若将式(15)中的β定义为负数，可以实现对室内脉冲响应的精确建模.此时，可以得到新的声压脉冲响应表达式，即

exp(-jωt)

(16)

此时，将式(16)进行傅里叶逆变换后，得到的时域状态下室内脉冲响应表达式为

hne(t)=F-1(Hne(ω))

(17)

3 矩形空间的算法验证

矩形空间是最典型也是最简单的一种三维模型，对其进行研究和讨论，有利于研究更为复杂的情形，因此，本文针对矩形房间进行了仿真.假设矩形房间天花板材质为石膏；地板材质为树脂.房间内材料的吸声系数如表1所示.

表1 房间内材料的吸声系数
Tab.1 Sound absorption coefficients of room materials

室内空间吸声系数四面侧墙0 90地面0 70顶棚0 99

以房间一角为原点，且以长为x轴、宽为y轴、高为z轴建立三维直角坐标系，可以得到仿真室内空间的相应坐标.表2为仿真室内空间的环境设置情况.

为了比较实验结果的准确性，在与实验环境参数相同的房间中，获取真实的室内脉冲响应，使之与测量结果进行比较.首先采用幅频特性以平坦谱的白噪声信号作为激励信号，然后记录由扬声器输出的激励信号在麦克风接收端的脉冲响应，最后通过解卷积的方法求得室内脉冲响应.为了获得较为准确的脉冲响应，可以通过多次记录数值并取平均值的方式得到真实的室内脉冲响应.

表2 仿真室内空间的环境设置
Tab.2 Environment arrangement of simulated room space

空间尺寸声源坐标/m麦克风坐标/m采样频率/kHz声速/(m·s-1)3 2m×4 0m×2 7m(2 0，3 0，2 0)(1 1，1 0，1 2)16343

图3为由不同方法获得的室内脉冲响应对比结果.由图3可知，与真实的室内脉冲响应相比(见图3c)，利用基本镜像源方法获得的室内脉冲响应的失真情况较为严重(见图3a)；利用变换频域与采用负反射系数方法得到的室内脉冲响应的精确度则相对较高(见图3b).因此，利用变换频域和采用负反射系数的镜像源方法可以获得更接近于真实室内脉冲响应的仿真结果.

为了进一步比较经过变换频域和采用负反射系数的镜像源方法与基本镜像源方法的仿真准确度，选取了5个不同位置的接收点进行了分析，且各接收点的位置如表3所示.

表3 各接收点位置Tab.3 Locations of various receiving points m

为了进行性能比较，采用归一化误差作为评价建模准确性的准则，且其表达式为

(18)

表4为由不同方法得到的室内脉冲响应误差对比结果.

表4 室内脉冲响应误差对比Tab.4 Comparison in errors of room impulse response dB

表4中改进方法1是针对频域方法的改进；改进方法2是针对频域与反射系数的改进.由表4可见，无论接收点位置如何改变，与基本镜像源方法相比，改进镜像源方法所得到的室内脉冲响应误差更小，且改进方法2获得的室内脉冲响应更接近于真实情况.

4 应 用

混响时间是厅、堂、馆音响效果设计中一个重要的客观评价参数.利用脉冲响应获得声场能量衰减曲线(EDC)后估计混响时间，已经成为近年来国家标准化组织(ISO)对各声学参数进行标定的标准测量方法.所谓能量衰减曲线是指声音在充满声场后突然停止，某一点的声能随时间的推移逐渐衰减而形成的曲线.依据能量衰减曲线可以估计混响时间[9-10].

根据式(17)并结合Schroeder方法可以得到能量衰减曲线(EDC)的表达式

(19)

式中，E为能量.

根据式(19)可以仿真得到能量衰减曲线，结果如图4所示.图4中3条曲线是在不同的反射系数下得到的仿真结果，且反射系数的具体设置如表5所示.由图4的能量衰减曲线可直接读取出能量从-5 dB变化到-25 dB的时间.

图4 能量衰减曲线Fig.4 Energy decay curve (EDC)

表5 室内各表面的反射系数Tab.5 Reflection coefficients of various room surfaces

表6为分别采用了基于变换频域和负反射系数的脉冲响应方法、Eyring方法与Sabine方法[11]后获得的混响时间数据.表6中的误差率是指估值与理想值作差取绝对值后，再与理想值作比后的数值.由表6可见，由脉冲响应方法获得的混响时间的准确率远高于另外两种方法.

表6 混响时间比较
Tab.6 Comparison in reverberation time

方法理想值/s估值/s误差率/%脉冲响应方法Eyring方法Sabine方法0 050 04840 100 082180 150 14440 050 045100 100 136360 150 212410 050 063260 100 138380 150 21543

5 结 论

本文以经典镜像源方法为基础，讨论了几种声场脉冲响应的具体算法，编制了仿真程序，给出了矩形空间声场的计算实例，并对直接影响模拟结果的几个重要问题进行了讨论，最后基于仿真得到室内声场脉冲响应，估计了能量衰减曲线和房间混响时间.仿真结果表明，利用基于变换频域和负反射系数的镜像源方法获得的脉冲响应的混响时间结果与理想值具有较好的一致性，且具有更高的准确性和可靠性.

[1]成烨.室内声场的脉冲响应分析 [D].太原：太原理工大学，2014.

(CHENG Ye.Impulse response analysis for room sound field [D].Taiyuan：Taiyuan University of Technology，2014.)

[2]马宏伟，董明.基于矩形换能器空间脉冲响应的相控阵次声场研究 [J].机械工程学报，2014，50(18)：36-42.

(MA Hong-wei，DONG Ming.Research on acoustic field of ultrasonic linear phased array based on spatial impulse response of rectangular planar transducer [J].Journal of Mechanical Engineering，2014，50(18)：36-42.)

[3]黄坤朋，赵越喆.应用时域有限差分法计算房间脉冲响应和声场参数 [J].土木建筑与环境工程，2013，35(3)：157-161.

(HUANG Kun-peng，ZHAO Yue-zhe.Calculation of room impulse response and acoustic parameters by finite-difference time-domain method [J].Journal of Civil，Architectural and Environmental Engineering，2013，35(3)：157-161.)

[4]陈妍.长空间的声学特性研究 [D].南京：南京大学，2011.

(CHEN Yan.A study on the acoustic of long enclosures [D].Nanjing：Nanjing University，2011.)

[5]曾向阳.室内声场计算机模拟发展40年(1968～2008) [J].电声技术，2008(10)：15-17.

(ZENG Xiang-yang.Forty years of development of room acoustic computer simulation [J].Audio Engineering，2008(10)：15-17.)

[6]姜锦云.基于麦克风阵列的近场和远场混合声源定位 [D].成都：西南交通大学，2010.

(JIANG Jin-yun.Mixed near-filed and far-filed speech source local localization based on microphone arrays [D].Chengdu：Southwest Jiaotong University，2010.)

[7]焦伟奇.基于虚源法的声矢量场建模与预报及软件实现 [D].哈尔滨：哈尔滨工程大学，2011.

(JIAO Wei-qi.Research on acoustic vector field calculation and analysis based on image source and software implementation [D].Harbin：Harbin Engineering University，2011.)

[8]矫姗姗.厅堂音质控制计算机仿真的误差分析 [D].青岛：中国海洋大学，2013.

(JIAO Shan-shan.The error analysis of computer simu-lation for the hall acoustics quality control [D].Qingdao：Ocean University of China，2013.)

[9]Li K，Lam P.Prediction of reverberation time and speech transmission index in long enclosures [J].The Journal of the Acoustical Society of America，2005，117(6)：3716-3726.

[10]Omoto A.Measurement of sound absorption coefficients in a reverberation room [J].Journal of the Acoustical Society of Japan，2012：68(12)：616-621.

[11]刘延昌.厅堂声场计算机模拟与实际测量结果对比分析研究 [D].邯郸：河北工程大学，2011.

(LIU Yan-chang.Comparative analysis of the acoustic simulation and the filed measurement’s results [D].Handan：Hebei University of Engineering，2011.)

(责任编辑：尹淑英 英文审校：尹淑英)

Simulation for impulse response of room acoustic field based on image source method

GUO Ying, YAN Mei-chen

(School of Information Science and Engineering, Shenyang University of Technology, Shenyang 110870, China)

In order to solve the problem that the accuracy in describing the sound quality characteristics of acoustic field based on the impulse response of room acoustic field is not high, the classical image source method (ISM) was taken as the core, and the improved algorithms based on the transformation frequency domain and reflection coefficient were discussed. In addition, the rectangular space was taken as the specific example, the receiving points and sound source locations were given, and the impulse response of room acoustic field was simulated with the computer simulation. According to different parameters, the various algorithms were analyzed. In addition, the energy decay curve (EDC) and room reverberation time (RT) were estimated. The results show that the improved algorithm based on transform frequency domain and reflection coefficient has higher accuracy and reliability.

room acoustic field; impulse response; sound quality characteristic; image source method; transform frequency domain; reflection coefficient; energy decay curve; reverberation time

2015-12-17.

国家自然科学基金资助项目(61102123)； 辽宁省教育厅优秀人才资助项目(LJQ2014011)； 辽宁省教育厅一般项目(L2014041)； 沈阳市科技局工业攻关项目(F14-073-2-00).

郭 莹(1975-)，女，辽宁铁岭人，副教授，博士，主要从事非高斯信号处理与参数估计等方面的研究.

16∶06在中国知网优先数字出版.

http：∥www.cnki.net/kcms/detail/21.1189.T.20160907.1606.012.html

10.7688/j.issn.1000-1646.2017.01.11

TN 912

A

1000-1646(2017)01-0055-06