土木建筑与环境工程  2013, Vol. 35 Issue (3): 157-161   PDF    
应用时域有限差分法计算房间脉冲响应和声场参数
黄坤朋1,2, 赵越喆1    
1. 华南理工大学 亚热带建筑科学国家重点实验室, 广州 510640;
2. 国光电器股份有限公司 博士后工作站, 广州 510800
收稿日期：2012-12-20
基金项目：国家自然科学基金（50938003）
作者简介：黄坤朋(1980-), 男, 博士(后), 主要从事室内声学研究, (E-mail)ohmycar@163.com。
赵越喆(通信作者), 女, 教授, 博士生导师, (E-mail)arzhyzh@scut.edu.cn。
摘要：在室内声学时域有限差分（FDTD）计算模型中，将房间的复阻抗边界用数字阻抗滤波器描述，并给出了门、窗和墙壁的数字滤波器复阻抗模型，应用该模型计算房间脉冲响应和声场参数。对一真实房间，将FDTD模拟计算得到的房间脉冲响应和声学参数与实际测量结果相比较，验证了在室内声学FDTD计算中，采用数字滤波器近似的复阻抗边界模型能较好地模拟房间中的声场。
关键词：室内声学    时域有限差分法(FDTD)    声学复阻抗边界    数字滤波器    
Calculation of Room Impulse Reponse and Acoustic Parameters by Finite-difference Time-domain Method
Huang Kunpeng1,2, Zhao Yuezhe1    
1. State Key Laboratory of Subtropical Building Science, South China University of Technology, Guangzhou 510640, P. R. China;
2. Postdoctoral Programme, Guoguang Electric Company Ltd., Guangzhou 510800 P. R. China
Received: 2012-12-20
Abstract: The administrator boundaries in room were modeled by digital filters of complex impedance in finite-difference time-domain (FDTD) method. Doors, windows and walls were described respectively by their complex impedance boundaries with digital filters. Impulse responses in a room were calculated with the FDTD model, and typical room acoustics parameters were obtained from the responses. The high correlation coefficients between the calculated and measured impulse responses and the good agreement in acoustic parameters show that sound field in room can be properly simulated by this FDTD model.
Key Words: room acoustics    finite-difference time-domain (FDTD)    acoustic impedance boundary    digital filter    

使用时域有限差分(FDTD)数值计算技术求解声波方程的关键是建立精确稳定的边界模型。理想刚性边界和全吸收边界[1-2]只能应用于一些特殊的边界情况。在室内空间中更为普遍的情况是，当声波经房间界面反射后，反射声波的幅值和相位都随频率发生变化。早期在房间声场模拟中使用的简化阻抗边界模型[3]只考虑声波的幅值变化，而忽略声波的相位变化。而复阻抗边界同时考虑声波入射到边界时幅值和相位的改变及其与频率的关系。

传统上，FDTD中的复阻抗边界模型由其等效的力学模型近似[4-6]。近年来，Kowalczyk等[7-8]将界面声阻抗用数字阻抗滤波器描述，与前者相比，数字滤波器法可以更加直接地从界面复阻抗的物理表达式求得其FDTD模型。Kowalczyk等[7]给出了房间复阻抗边界无限长脉冲响应数字滤波器(IIR)的FDTD表达式，并通过与理论值对比，证明该方法可准确地模拟声波局部作用表面的声反射。本文尝试将该技术应用于室内声场三维计算机仿真，通过与实测脉冲响应及声场参数的对比，探讨该模型的有效性。

1 复阻抗边界的FDTD计算模型
1.1 FDTD方程

当声波振幅较小，且声波在空气中传播时的能量衰减可被忽略时，均匀的理想空气介质中的声波方程为：

$ \frac{{{\partial ^2}p}}{{\partial {t^2}}} = {c^2}{\nabla ^2}p $ (1)

其中p为声压，Pa；c为声波在介质中的传播速度，m/s；∇为拉普拉斯算符。

将方程(1)离散化，对非交错网格，其三维空间网格节点处声压的FDTD表达式为：

$ \begin{array}{*{20}{l}} {{p^{n + 1}}\left( {i,j,k} \right) = {{\left( {c\frac{{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}t}}{{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}h}}} \right)}^2}}\\ {\left[ {{p^n}\left( {i + 1,j,k} \right) + {p^n}\left( {i - 1,j,k} \right) + {p^n}\left( {i,j + 1,k} \right) + {p^n}\left( {i,j - 1,k} \right) + {p^n}\left( {i,j,k + 1} \right) + {p^n}\left( {i,j,k - 1} \right)} \right]}\\ { + 2\left( {1 - 3{{\left( {c\frac{{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}t}}{{\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}h}}} \right)}^2}} \right){p^n}\left( {i,j,k} \right) - {p^{n - 1}}\left( {i,j,k} \right)} \end{array} $ (2)

其中，pn(i, j, k)表示在网格节点(i, j, k)上n时刻的声压，n是离散化的时间变量。Δh和Δt分别为网格的单元步长和时间步长。某点下一时刻的声压pn+1(i, j, k)决定于这点上一时刻的声压、这点周围6个相邻点(如图 1)的上一时刻声压，以及这点的再上一时刻声压。因此从时间上来看，总共需要3套存储空间，分别用于存储n+1、n以及n-1时刻的声压值。

图 1 三维坐标下的声压网格点分布

1.2 复阻抗边界IIR数字滤波器的设计

对局部作用表面，界面上某一质点的振动速度只与空气中和这个点最相邻的质点声压有关，而与其它质点声压无关[9]。同时，考虑到复阻抗边界为频率的函数，具有幅度响应和相位响应的特性，而滤波器也有幅度和相位的频率响应，因此可以使用滤波器模拟复阻抗边界[7]。

首先，设计边界的法向反射系数R的IIR模拟滤波器：

$ \begin{array}{*{20}{l}} {\;\;\;\;\;R\left( s \right) = }\\ {\frac{{b'\left( 1 \right){s^n} + b'\left( 2 \right){s^{n - 1}} + \ldots + b'\left( n \right){s^1} + b'\left( {n + 1} \right)}}{{a'\left( 1 \right){s^n} + a'\left( 2 \right){s^{n - 1}} + \ldots + a'\left( n \right){s^1} + a'\left( {n + 1} \right)}}} \end{array} $ (3)

其中s为拉普拉斯变量，b’和a’为滤波器系数。边界的特性声阻抗率ξ(s)的IIR滤波器可由下式计算得到：

$ \begin{array}{l} \;\;\;\;\xi \left( s \right) = \frac{{Z\left( s \right)}}{{{\rho _0}c}} = \frac{{1 + R\left( s \right)}}{{1-R\left( s \right)}} = \\ \frac{{b\left( 1 \right){s^n} + b\left( 2 \right){s^{n-1}} + \ldots + b\left( n \right){s^1} + b\left( {n + 1} \right)}}{{a\left( 1 \right){s^n} + a\left( 2 \right){s^{n-1}} + \ldots + a\left( n \right){s^1} + a\left( {n + 1} \right)}} \end{array} $ (4)

利用双线性变换[10]将(4)式离散化，得到其数字形式：

$ \xi \left( z \right) = \frac{{{b_0} + {b_1}{z^{-1}} + \ldots + {b_{n-1}}{z^{-(n - 1)}} + {b_n}{z^{ - n}}}}{{{a_0} + {a_1}{z^{ - 1}} + \ldots + {a_{n - 1}}{z^{ - (n - 1)}} + {a_n}{z^{ - n}}}} $ (5)

其中，a0，ai，b0，bi为数字滤波器的系数，n为滤波器的阶数。z和s的关系为：

$ s = \frac{2}{T}\frac{{z-1}}{{z + 1}} $ (6)

式中T为计算模型中的采样周期，等于公式(2)中的Δt。

将(5)式应用于FDTD计算模型。例如在右边界中，即p(i, j, k)在边界上，p(i+1, j, k)则在边界的右边，边界上声压p(i, j, k)的FDTD表达式为[7-8]：

$ \begin{array}{l} {p^{n + 1}}\left( {i,j,k} \right) = \\ \left\{ {\left[ {2{p^n}\left( {i - 1,j,k} \right) + {p^n}\left( {i,j + 1,k} \right) + {p^n}\left( {i,j - 1,k} \right) + {p^n}\left( {i,j,k + 1} \right) + {p^n}\left( {i,j,k - 1} \right)} \right] + } \right.\\ \left. {{m^2}\frac{{{m^2}}}{{{b_0}}}{g^n} + \left( {\frac{{m{a_0}}}{{{b_0}}} - 1} \right){p^{n - 1}}\left( {i,j,k} \right)} \right\} \cdot {\left( {1 + \frac{{m{a_0}}}{{{b_0}}}} \right)^{ - 1}} \end{array} $ (7)

式中

$ m = c\frac{{\mathit{\Delta }T}}{{\mathit{\Delta }h}} $ (8)
$ {g^n} = \sum\limits_{i = 1}^N {\left( {{b_i}{x^{n-i}}-{a_i}{y^{n-i}}} \right)} $ (9)
$ {x^n} = \frac{{{a_0}}}{{\lambda {b_0}}}\left[{{p^{n + 1}}\left( {i, j, k} \right)-{p^{n-1}}\left( {i, j, k} \right)} \right] -\frac{{{g^n}}}{{{b_0}}} $ (10)
$ {y^n} = \frac{1}{{{a_0}}}\left( {{b_0}{x^n} + {g^n}} \right) $ (11)

式(2)和式(7)分别为声波在房间中和房间界面上的FDTD计算表达式。在模型中加入激励声源，即可计算房间的脉冲响应，并求出声学参数。

对于复阻抗滤波器的设计，若是穿孔板、多孔吸声材料等则可以根据已有的阻抗模型[9]进行设计。近乎刚性的界面则可令其吸声和反射声压相位的变化均近似为零。对复杂界面，可通过驻波管测得其复阻抗或复数形式的反射系数。

2 计算实例
2.1 房间模型

1个实际房间的长、宽、高尺寸分别6.46 m×3.56 m×3.80 m。铝合金玻璃窗和木门的大小分别为1.88 m×2.4 m和1.5 m×3.19 m，位于房间的两端墙上，见图 2。声源点S高出地面1.28 m，6个受声点R1~R6高出地面1.2 m，它们在水平面上的坐标分别为(单位：m)：R1(1.33, 1.14)，R2(1.33, 1.89)，R3(2.45, 1.10)，R4(2.46, 1.88)，R5(3.33, 1.11)，R6(3.33, 1.86)。

图 2 计算的房间及声源点和受声点

窗和木门声阻抗模型可近似表达为[9]：

$ Z\left( s \right) = Ms + {\rho _0}c $ (12)

其中：ρ0为空气的密度，取为1.21 kg/m3，c等于343 m/s，M为窗或木门的面密度，分别为7.5和12 kg/m2，拉普拉斯变量s=jω。

由式(4)，它对应的反射系数为：

$ R\left( s \right) = \frac{{Z/{\rho _0}c-1}}{{Z/{\rho _0}c + 1}} = \frac{{0 + Ms}}{{2{\rho _0}c + Ms}} $ (13)

将R(s)离散化，采用双线性变换的方法将其转化成数字滤波器的形式：

$ R\left( z \right) = \frac{{M-M{z^{-1}}}}{{{\rho _0}cT + M + ({\rho _0}c-M){z^{ - 1}}}} $ (14)

其对应的特性声阻抗率

$ {\xi _w}\left( z \right) = \frac{{{\rho _0}cT + 2M + ({\rho _0}cT-2M){z^{-1}}}}{{{\rho _0}cT + {\rho _0}cT{z^{-1}}}} $ (15)

公式(15)的滤波器系数为

$ {b_0} = {\rho _0}cT + 2M, {b_1} = {\rho _0}cT-2M $ (16)
$ {a_0} = {\rho _0}cT, {a_1} = {\rho _0}cT $ (17)

房间其余墙面为砖墙抹灰，地面为水磨石，近似刚性界面，为此可设计一个低吸声、声阻抗相位响应近似为零的IIR滤波器。对应地面的数字滤波器系数为：b0=1.990 6，b1=-0.702 0，a0=0.009 4，a1=-0.004 3；对墙面，构建的数字滤波器系数为：b0=1.985 1，b1=-0.749 6，a0=0.014 9，a1=-0.006 0。在2 kHz的频率范围，所设计的墙面和地面的数字滤波器所对应的法向吸声系数以及相位响应如图 3所示。吸声系数在0.01~0.03之间，复阻抗相位频率响应接近于零，即反射声波经界面反射后相位基本不发生变化。

图 3 水泥地面和抹灰砖墙的声阻抗率数字滤波器对应的法向吸声系数及相位响应

2.2 计算步长和脉冲声源

使用FDTD模型计算声场时，为使计算过程不发散，Δh应小于最大研究波长的1/10[11]，结合所用的硬件水平，本计算中取Δh=0.05 m。此外，根据三维稳定性条件[4]：

$ \mathit{\Delta }t \le \frac{{\mathit{\Delta }h}}{{\sqrt {3c} }} $ (18)

取时间步长dt=0.084 2 ms。按此参数设置，至少可以确保计算模型适用于700 Hz以下的频率。

高斯脉冲是声场仿真中常用的脉冲声源之一，它的时域表达式为：

$ s\left( t \right) =- A(t- {t_0}){\rm{exp}}\left[{-{{(t-{t_0})}^2}/\tau } \right] $ (19)

其中t0表示高斯脉冲波节出现的时刻，A和τ为常数，分别用于调节脉冲的幅度和宽度。本文的计算中，取A=0.25，τ=20Δt，t0=28Δt。高斯脉冲信号如图 4所示。

图 4 激励高斯脉冲声源的时域波形及其频率响应

由于采用的高斯脉冲总有一定的时间宽度，即其所对应的频响曲线不是覆盖所有频率范围的平直曲线，因此使用高斯脉冲作为激励声源时，需对所得到的脉冲响应进行修正。

房间中某点的响应在时域上可由下式表达：

$ y\left( t \right) = s\left( t \right) * h\left( t \right) $ (20)

其中s(t)为激励声源信号，h(t)为房间脉冲响应，符号*为卷积运算符。由上式有：

$ h\left( t \right) = Is\left( t \right) * y\left( t \right) $ (21)

其中，Is(t)的傅立叶变换为1/S(ω)，S(ω)为声源的频域表达式。因此只要得到声源的修正滤波器函数Is(t)，将其与受声点接收到的脉冲响应y(t)卷积，即可排除声源特性的影响。

计算中所采用的修正滤波器Is(t)的时域波形和频率响应示于图 5。修正补偿滤波器同时包括了截止频率为90 Hz的高通滤波以及截止频率设计为3 kHz的低通滤波。

图 5 修正滤波器的时域波形及其频率响应

2.3 实测

为验证计算结果的正确性，对上述房间的脉冲响应进行了实测。测试所用仪器设备如图 6所示。测试时声源点和受声点的位置与计算模型相同，如图 2。实验采用MLS信号作为激励声源来测试房间的脉冲响应。由计算机(01 dB测试软件)产生数字MLS噪声信号，然后由Symphony分析仪做D/A转换后得到模型信号，经功放(B & K2716)放大后由正十二面体扬声器(B & K4296)发声；受声点处的传声器(B & K4191)录得声音信号后经Symphony分析仪做A/D转换，再由01 dB软件进行解卷积处理，最后得到房间脉冲响应。为与计算结果进行对比，将实测得到的脉冲响应经过了截止频率为90 Hz的高通滤波和截止频率为3 kHz的低通滤波。

图 6 测试系统图

2.4 结果分析

由于脉冲响应的均方根可以直观描述声场的变化情况，由下式计算脉冲响应的每1 ms时长的均方根RMS[6]：

$ RMS\left( \tau \right) = \sqrt {\frac{{\int_\tau ^{\tau + {\tau _0}} {{p^2}\left( t \right){\rm{d}}t} }}{{{\tau _0}}}} $ (22)

式中p(t)为归一化的脉冲响应，τ0=1 ms。由此得到的受声点1的脉冲响应对比如图 7所示，它们的相关系数为0.77。受声点R2~R6的实测脉冲响应和模拟脉冲响应之间的相关系数分别为：0.76，0.79，0.75，0.79和0.76。从各受声点的实测和模拟得到的脉冲响应的形状以及它们之间的相关系数可知，文中的FDTD模型能较好地模拟三维室内声场情况。

图 7 受声点1实测和计算的脉冲幅度变化曲线对比

本文还对重要的室内声学参量，包括早期衰变时间EDT、混响时间T30以及明晰度C80和清晰度D50等进行了实验与计算数据的对比。计算和实测得到的房间6个点的平均声学参数对比如图 8~11所示。图中的2条虚线为±JND(Just noticeable difference)[12]。从图中可以看出，EDT和T30在125 Hz和1 kHz倍频带上，实测值和模拟计算值存在差异，而在其余频带则符合较好。界面声阻抗是影响计算结果精度的决定因素。文中声波从各方向入射到抹灰墙、地板和天花的阻抗均由其法向阻抗近似；木门和玻璃窗由振动吸声模型描述；但阻抗的物理模型[9]不能完全准确描述真实的界面阻抗；此外，计算模型忽略了实际房间的建筑细节，这些都是造成实测值与计算值差别的原因。对于C80和D50，模拟值基本处于实测值的一个JND范围内。

图 8 实测和计算得到的的EDT平均值

图 9 实测和计算得到的T30平均值

图 10 实测和计算得到的C80平均值

图 11 实测和计算得到的D50平均值

计算结果也表明，在上述Δh和Δt的取值情况下，对于中心频率为2 kHz的倍频程及其以下频率，数值计算过程没有发散，这必须在满足稳定性条件即式(18)下进行多次取值计算，找到一组合理的Δh和Δt。

3 结论

声波在房间界面发生反射时，声能衰减的同时其相位也发生改变，且与声波的频率有关。房间界面的声学特性可用复阻抗表达。将复阻抗边界的IIR数字滤波器模型应用于房间三维声场的FDTD模拟，在2 kHz及其以下的频率范围内得到了与实测较为接近的脉冲响应，且声场参数与实测值相符。

参考文献
[1] Mur G. Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations[R]. IEEE Transactions on Electromagnetic Capability, 1981, EMC-23(4):377-382.
[2] Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 1994, 114(2): 185–200. DOI:10.1006/jcph.1994.1159
[3] Yokota T, Sakamoto S, Tachibana H. Visualization of sound propagation and scattering in rooms[J]. Acoustical Science and Technology, 2002, 23: 40–46. DOI:10.1250/ast.23.40
[4] Botteldooren D. Finite-difference time-domain simulation of low-frequency room acoustic problems[J]. Journal of the Acoustical Society of America, 1995, 98(6): 3302–3308. DOI:10.1121/1.413817
[5] Fung K Y, Ju H B. Time-domain impedance boundary conditions for computational acoustics and aeroacoustics[J]. International Journal of Computational Fluid Dynamics, 2008, 18(6): 503–511.
[6] Sakmoto S, Nagatomo H. Calculation of impulse responses and acoustic parameters in a hall by the finite-difference time-domain method[J]. Acoustical Science and Technology, 2008, 29(4): 256–265. DOI:10.1250/ast.29.256
[7] Kowalczyk K, Walstijn M V. Modeling frequency-dependent boundaries as digital impedance filters in FDTD and K-DWM room acoustics simulations[J]. Journal of the Audio Engineering Society, 2008, 56(7/8): 569–583.
[8] Kowalczyk K, Walstijn M V. Formulation of locally reacting surfaces in FDTD/K-DWM modeling of acoustic spaces[J]. ACTA Acustica United with Acustica, 2008, 94: 891–906. DOI:10.3813/AAA.918107
[9] Kuttruff H. Room acoustics[M]. 5th ed. London and New York: Spon Press, 2009: 164-165.
[10] Parks T W, Burrus C S. Digital Filter Design[M]. New York: John Wiley & Sons, 1987: 209-213.
[11] 李太宝. 计算声学:声场的方程和计算方法[M]. 北京: 科学出版社, 2005.
[12] ISO 3382-1:Acoustics-Measurement of room acoustic parameters-Part 1:Performance rooms[S]. Annex A Auditorium measures derived from impulse responses, 2006.
X
    应用时域有限差分法计算房间脉冲响应和声场参数
    黄坤朋, 赵越喆