杨 勃 徐禄文 王海涛 陈建明 徐大然

(1国网重庆市电力公司电力科学研究院 重庆 401123)

(2西北工业大学航海学院 西安 710072)

(3国网重庆市电力公司信息通信分公司 重庆 401121)

(4重庆市送变电工程有限公司 重庆 400039)

0 引言

声场仿真预测是声学领域中的一项重要研究内容，在噪声控制、声品质设计以及声环境监测等应用中均可发挥重要的作用。室内环境中的声场仿真受声波反射因素的影响，相较于自由场环境下的仿真具有更高的复杂度，特别是对于户内变电站主变室等低频噪声突出且噪声分布频段较宽的室内环境，其仿真还需要对干涉、衍射等波动现象进行模拟，对仿真方法提出了更高的要求；另外，在户内变电站主变室、大型厅堂等场所中，常存在着较大范围的周期排列结构，例如设备散热器、通风百叶窗、声学扩散体等，这些界面上的声散射具有较强的规律性[1−2]，其模拟准确性对于最终的声场仿真精度也起到重要的限制作用。因此，发展一种能够对周期界面声散射进行准确模拟，且具有宽频高精度的声学仿真方法，对于室内声学设计、噪声控制、声学监测等应用的支撑具有重要的意义。

经典的声场仿真预测方法主要可以分为波动声学方法、统计声学方法以及几何声学方法3类。3类方法目前都有较为成熟的研究，但都存在一定的应用限制。其中波动声学方法主要适用于低频段声场仿真，常见方法包括有限元法(Finite element method,FEM)、边界元法(Boundary element method,BEM)等。波动声学方法通过求解波动方程而获得声场分布结果，是一种具有严格数学理论支撑的仿真方法，可以获得声场分布的精确解。但是，由于波动声学方法需要对空间进行网格离散化处理，受网格尺寸与波长关系影响，当处理较高频段的声场仿真问题时，会造成计算规模的急剧上升，特别是当空间中包含周期界面时，需对所有周期子结构进行非常细化的空间离散处理，会进一步影响计算精度并严重降低计算效率[3]。统计声学方法通常特指统计能量分析法(Statistic energy analysis,SEA)，此方法在高声模态密度条件下分析功率流实现对声场的预测，仅适用于高频仿真问题，而在低频段无法适用，而且在实际应用中，统计能量分析法所需参数较为复杂且不容易获取，例如各类耦合因子等，因此应用受限较多[4]。几何声学方法是另一类重要的声学仿真方法，是指通过模拟声波物理传播过程或物理规律而实现对声场仿真的方法，主要包含声线跟踪法(Ray tracing method,RT)、虚声源法(Image source method,ISM)等。其中声线跟踪法由于具有物理描述直观、建模简易等特点，自其在20世纪60年代被提出后就获得广泛关注[5]。

早期的声线跟踪法中，由于声线无法模拟声波的衍射、干涉等波动现象，导致其只能适用于高频声场仿真。为了在简便的几何声学框架下提高对低频段声场的模拟精度，Gensane等[6]提出了一种能考虑不同反射波之间声压干涉叠加的干涉模型，并将其应用于几何声学方法之中，以提高低频仿真精度；之后，Lemire等[7]在对点声源于封闭空间内的声传播研究的基础上，提出了相干几何法，大大提高了几何声学方法在低频段的仿真精度。而且，由于无需对空间几何模型进行频率相关的模型细化，在计算效率方面相较于波动声学方法具有更加明显的优势，因此，相干几何法已成为涵盖低频仿真精度要求的全频段声学仿真常用方法。

虽然相干几何法目前已成为宽频段声学仿真的一种常用方法，但是在复杂声散射存在的情况下，例如周期散射界面存在时仿真精度仍有不足，主要体现在低频段精度差，其原因是：在低频段，声波在周期结构表面实际并不会发生明显散射，但是传统相干几何法中的声线会受周期结构几何形状影响而向无规方向反射，导致声场趋向于扩散场，与实际情况产生偏差。

针对周期散射界面存在条件下的宽频空间声场仿真问题，本文发展了一种基于声线迭代散射模型的相干声线跟踪法。此方法以经典的相干声线跟踪法为基础，在周期散射结构界面反射过程中，声线将会在散射作用下发生迭代分裂，通过格栅方程预判结构尺寸与波长的关系，给出子声线的反射方向，并通过散射系数估计不同子声线反射后的能量；同时，借助于散射模型的引入，可将复杂的周期结构等效为更为简便的模型，从而降低了几何建模难度及运算时间。与边界元法及传统相干几何法的比较验证了此模型在宽频段，特别是在传统几何法精度较差的低频段，具有较高的仿真精度。

1 基于声线迭代散射模型的相干声线跟踪法

1.1 声线迭代散射模型

在实际空间环境中，常存在各类周期结构，例如户内变电站主变室中的散射设备、大型音乐厅中的周期扩散体、各类厂房空间中的通风百叶窗等，均可视为周期结构。假设一包含周期结构的室内环境如图1所示，其中S为声源，R为接收区域，Γ为包含周期结构的壁面。

图1 包含周期散射界面室内环境示意图Fig.1 The schematic of room environment containing the periodic scattering surface

在声线跟踪法中，所有频段的声波均采用声线的形式模拟。对于传统的相干几何法，当空间中存在周期结构时，低频声线可能会受周期结构的轮廓影响而发生如图1所示的S′1方向的反射。但实际上，当低频声波波长较大时，周期结构不会对低频声波产生散射，即低频声波将以镜面反射的方式与界面发生作用，如图1所示的S1反射，此时传统的相干几何法由于对声线方向模拟不准确，导致提高了声场的扩散性，造成精度下降。针对此问题，本文发展了一种声线迭代模型，此模型中，将周期结构在几何上进行简化，同时当声线在周期结构散射界面上发生散射时，将依据周期散射定理发生分裂，分裂的形式与频率、周期结构尺寸等因素密切相关，以下为模型的具体推导。

假设一根声线入射到具有周期轮廓的散射体上，首先仅考虑声线在周期结构垂直切面xOy内的情况，如图2所示。散射体的单元周期长度为L，凹槽长度为l，高度为H，入射波的声压可表示为

式(1)中，θi为入射角，k=2π/λ为波数，λ为波长。与光栅类似，周期型的声学散射体同样具有可用光栅方程所表示的周期散射效应，在xOy平面内，散射波的阶次及方向可根据周期散射定理求得[8−9]：

式(2)中，θs为散射角，λ为入射波的波长，n为散射波的阶次，n=0即代表镜面反射波。

式(2)表明，当一束声波入射到周期散射界面上时，会被界面散射为多束反射波。由于需要满足−1≤sinθs≤1的条件，因此散射波阶次n的取值会受声波频率、周期结构尺寸以及入射角度多重因素影响。当声波频率较低时，波长较大，散射波阶次较低，从而反射波的数量较少，但至少会有一束镜面反射波；而当频率较高时，散射波阶次会较高，反射波的数量也会相应较多。依据此理论，本文提出了声线迭代来模拟声波在周期结构表面上的散射作用。此算法的核心思想是将图2中边界处的周期结构视为平面结构，这样可以降低建模的难度，也降低跟踪过程中的界面处理复杂度，当声线入射到周期结构上时，声线会依据式(2)迭代分裂为多根子声线，此后仍按照经典方法对各个声线继续跟踪。

图2 入射声波在周期散射界面散射示意图Fig.2 Scattering of sound waves incident on the periodic surface

声线经过周期结构界面散射后，子声线的数量与n的取值相关，其中第n阶子声线的散射方向为

根据此散射方向即可得到此子声线跟踪方向的方向向量。

然后，将上述xOy平面内的入射情况拓展到三维情况，如图3所示，当声线在xOy平面以外入射到周期结构上时，将其首先投影在xOy平面以内，如图3中红色虚线所示，求解入射声线与其投影夹角γ；然后，按式(3)确定散射子声线在xOy平面内的反射角θs；以此为基础，确定xOy平面以外的实际散射子声线，其自然坐标系下的反射角度通过(θs,r)求得，γ为反射子声线与其在xOy平面投影的夹角，等于入射声线与其在xOy平面投影的夹角γ。

图3 三维空间内的入射声波在周期界面散射示意图Fig.3 Scattering of sound waves incident on the periodic surface in three dimension space

1.2 相干声线跟踪

通过上述方式确定散射子声线的反射角度后，即确立了周期结构的声线散射模型，此后可采用经典的相干声线跟踪法完成所有声线的跟踪以及统计。对于每一根声线及其子声线，跟踪其传播路径，当其达到所设置的测点时，即将其记录，若一直未到达测点，则在能量低于某一阈值后对其停止跟踪。最终，在测点处所形成的声场可表示为

式(4)中，等式右侧第一项表示直达声的贡献，第二项表示混响声的贡献，其中E0表示声源的源强，ddir表示声源到测点的直达距离，di表示到达测点的第i根声线所经过的距离，j为虚数单位，N表示所有到达测点处声线或子声线的数量，Qi为声线从声源到达测点之间经过的所有界面的声反射系数，可以由声线传输路径中每次声反射对应的单次反射系数依次相乘得到：

式(5)中，M为声线由声源到测点之间所经过的反射次数，Qm为每次反射的反射系数，在Lemire相干模型中，利用无限大界面上球面波反射场的近似解来表示，由于声线在周期结构界面上有可能会发生散射，因此此处将其定义为传统反射系数与某阶散射系数的乘积：

式(6)中，smn为声线发生第m次反射时第n阶散射波的散射系数，即可确定每一个子声线在新的传播路径中所携带的初始能量。散射系数的计算值可根据周期结构中一个子单元凹槽的性质，通过有限元方法求解[10−11]，当周期结构的轮廓为矩形结构时，可采用更为简便的直接估算方法[12]得到其散射系数，此处散射系数的计算也表明，本文所推导的迭代散射模型实际上可适用于具有任意轮廓的周期结构的情况；Rm为第m个反射面上的平面波反射系数，表达式为

式(7)中，θm表示声线入射到第m个反射面上时与其法线的夹角，βm表示第m个反射面的法向比声导纳。式(6)中，F(w)为界面损失因子，其表达式为

式(8)中，erfc(·)为比例辅助误差函数，而w为数值距离参数，与声线所经过的距离、入射角以及相应的边界有关：

通过上述过程，即可实现对空间内的相干声线跟踪。

1.3 适用频段设定

本文所发展的散射模型一个关键处理是在界面上按照散射性质将原始声线进行迭代处理，在此过程中，原始声线有可能经过反射后会迭代为多根子声线，此处分析了声线数量增长对计算效率的影响。

在传统的相干几何法中，选定一根声线，假设此声线不被吸收，当经过m次反射后，其能量衰减到限定能量值以下时，即停止对这根声线的跟踪。此过程可以表示为

式(10)中，E0为声线初始携带能量，ET为跟踪过程中的限定能量值，α为壁面的吸声系数，m为反射次数。对于这根声线，从其发出到能量过小截止跟踪后，其跟踪总次数为

式(11)中，ℜ{·}表示取大于内部函数值的最大正整数。

对于本文方法，假设一个选定声线每次反射都会发生迭代，迭代后的数量假设为固定值β，即一个声线会分裂为β根子声线；且每根子声线所携带的能量相同，即声线反射后的能量可以平均分配在每根子声线上。此时对于最终的一个子声线，其停止跟踪的过程可以表示为

式(12)中，各符号意义与式相同。根据此式，m值可表示为

对于选定初始声线，由于每次反射后均会变为更多的子声线，是一个等比增长的过程，因此，其跟踪总次数为

通过一特例对跟踪次数进行对比。假设限定能量值与初始能量值之比为0.2，那么在吸声系数从0.1～0.5变化条件下的跟踪次数如图4所示，其中本文方法每次反射后的迭代子声线数量β分别为2、3。

图4 不同方法的声线跟踪次数统计Fig.4 Statistics of ray tracing times in different methods

从图4所示跟踪次数来看，本文方法总体上相较于传统相干几何法会产生较多的跟踪次数，特别是当较大且吸声系数较小时，跟踪次数会严重影响计算效率。

为了能够兼顾计算精度与计算效率，依据室内声学理论[13]，本文对散射模型的使用频段进行了限定，适用于迭代模型的频段为

式(15)中，T为空间的混响时间，V为空间体积。根据室内声学理论，当频率高于此频率时，声场将以扩散性质为主，此时，传统的相干几何法也可以得到较为准确的结果。

2 数值验证

为了验证本文所发展模型的正确性，在一个空间中进行了仿真验证试验。此空间的尺寸为2.6 m×2.2 m×2.0 m，由于尺寸较小，会发生明显的波动效应，以此来验证本文方法对于低频波动效应模拟的正确性。在其一个壁面上布置有矩形轮廓周期结构，周期结构轮廓的周期长度L为0.2 m，高度H为0.2 m。本算例中，所有界面的声阻抗设置为80ρ0c0，即法向比声阻抗为80。算例中设置了一个点声源，其位置为图5(a)中橘色圆点所示的(2.4,0.2,0.2)m处，另外设置了9个测点，坐标分别为(0.6,1.8,0.3)m、(0.6,1.6,0.3)m、(0.6,1.4,0.3)m、(0.6,1.8,0.7)m、(0.6,1.6,0.7)m、(0.6,1.4,0.7)m、(0.6,1.8,1.1)m、(0.6,1.6,1.1)m、(0.6,1.4,1.1)m。

图5 算例1封闭空间及其网格划分示意图Fig.5 The schematic of room environment 1 and its mesh for FEM

此算例首先运用Virtual.Lab软件的有限元法完成了计算，网格尺寸为40 mm，如图5(b)所示。根据有限元法网格应满足上限计算频率对应波长1/6的经验定理，此网格条件下，有限元法可达到的上限计算频率为1416 Hz。为了保证参考方法的准确性，本文选择1000 Hz作为上限对比频率。此外，还利用传统的相干声线法计算了声源到所有测点的频率响应。经过计算，有限元法、本文方法以及传统相干几何法在30～1000 Hz所得到的频率响应如图6所示，此处给出了测点1处的对比结果。

图6所示的频响曲线对比表明，相较于有限元法，传统相干几何法在低频段存在较为明显的误差，而本文方法则具有更高的精度。其中，在400 Hz以内的频段，本文方法所得到的结果与有限元法的结果存在细微的频率偏移，即所得模态频率与准确的模态频率存在一定的误差，但是在频率曲线的结构方面吻合度较高，而传统相干几何法则除了较为明显的频率偏移外，曲线的幅值与有限元法结果也存在一定误差，这说明在周期结构存在的空间内，本文方法相较于传统的相干几何法在低频段的精度上有较大的提高。在400 Hz以上的频段，本文方法以及传统相干几何趋向于与有限元法结果一致，误差均较低。此对比说明当空间内存在周期类型结构时，传统相干几何方法由于未能准确地模拟周期结构的散射特性，导致其误差较高，而本文方法对于散射的处理则更符合实际情况，具有较高的精度。

图6 测点1处不同方法的频率响应对比Fig.6 Comparison of frequency responses of different methods at measuring point 1

为了分析所有测点处的误差，统计了9个测点处的模态频率偏移平均误差。需要说明的是，此处之所以主要对比模态频率的误差，是因为在低频段模态频率的位置是一项基础评价指标，只有当模态频率处于同一位置时，不同方法才可以比较幅值的关系，否则会因为频响曲线在频率轴上的偏移造成幅值无法比较。本文使用式(16)方法统计了本文方法以及传统方法相对于有限元法的模态频率平均误差，结果如图7所示。

图7 本文方法与传统相干几何法的平均误差对比Fig.7 Comparison of average errors of the proposed method and the traditional coherent method

从图7所示的结果来看，在各个测点处，本文方法相较于传统相干几何法均体现出了更高的精度，在30～1000 Hz范围内，本文方法在所有测点处的平均误差均小于2 Hz，而传统相干几何法则比本文方法平均误差要高，进一步证明了本文方法的正确性。

为了进一步分析本文方法对空间声场幅值分布预测的正确性，统计了不同方法在各个测点上的声压级分布对比，限于篇幅，此处给出了100 Hz、300 Hz、500 Hz、700 Hz等4个示例性频率的对比图，如图8所示；另外以有限元法为标准参考方法，统计了本文方法以及传统相干几何法在各个测点上的声压级分布对比误差，如图9所示。

从图8所示的测点声压级对比以及图9所示的误差对比来看，在100 Hz时，本文方法所得结果与有限元法比较接近，一方面是在各个测点的变化趋势上具有良好的一致性，另一方面是在幅值上具有良好的一致性，而传统几何方法虽然在变化趋势上较为接近，但是在幅值方面具有明显的误差，主要表现为在各个测点上的起伏性较低。这说明，在低频处，真实声场分布并不是均匀的，即扩散性较低。传统的相干几何法在模拟声线跟踪的过程中，受周期结构影响，在原本不会发生散射的位置模拟发生了声线的散射，在一定程度上增加了声场的扩散性，从而导致其精度较低，而本文方法则能够准确模拟出散射的情况，从而具有较高的精度。当频率升高时，声场扩散性随之升高，例如图8(b)、图8(c)、图8(d)所示的300 Hz、500 Hz、700 Hz的情况，声压级的变化范围趋向于减小，本文方法均达到了良好的仿真效果。当频率较高时，传统方法因为可以模拟声场的这种扩散性，从而具有良好的模拟精度。除上述4个示例性频率之外，从整个频段的各个测点声压级分布来看，本文方法由于在低频段能够更加准确地模拟周期结构的散射规律，因此相较于传统相干几何法具有更高的精度，随着频率的升高，传统方法对于周期结构散射现象的模拟逐渐适用，精度提高，同时本文方法仍具有较好的精度。

图9 100 Hz、300 Hz、500 Hz、700 Hz时处各测点声压级分布误差Fig.9 Errors of sound field distributions of the different methods at 100 Hz,300 Hz,500 Hz and 700 Hz

为了验证本文方法在更加复杂的周期结构场景中的仿真精度，本文在上一个算例模型基础上，在模型空间内部增加周期结构，如图10所示。增加的周期结构处于与原周期结构相对的壁面上，新的周期结构L为0.1 m，高度H为0.1 m。本算例中，所有界面的声阻抗仍设置为80ρ0c0。此算例中设置了一个点声源，其位置为图5(a)中橘色圆点所示的(2.4,0.2,0.2)m处，另外设置了1个测点，坐标为(0.6,1.8,0.3)m。经计算，有限元法、本文方法以及传统相干几何法在30～1000 Hz所得到的频率响应如图11所示。

图10 算例2封闭空间及其网格划分示意图Fig.10 The schematic of room environment 2 and its mesh for FEM

图11 测点处不同方法的频率响应对比Fig.11 Comparison of frequency responses of different methods at the measuring point

图11所示的频响曲线对比表明，与算例1类似，本文方法在低频段相较于传统相干几何法具有更高的精度，虽然本文方法与有限元法对比存在细微的频率偏移，但是在频率曲线的结构方面吻合度较高。在较高频段，本文方法以及传统相干几何趋向于与有限元法结果一致，误差均较低。此对比说明当空间内存在较为复杂的周期类型结构时，本文方法依然能够给出良好的仿真结果。

3 结论

针对周期结构存在条件下的室内宽频声场仿真问题，本文发展了一种基于迭代散射模型的相干声线跟踪法。此方法在一定频段范围内将周期散射结构视为平面，然后基于周期散射定理，依据周期结构的轮廓尺寸以及声波波长、入射角条件确定不同的散射波，此时原始声线将会根据散射波情况迭代分裂为多个子声线，此后继续按照经典相干法对子声线展开跟踪，直到完成所有声线的跟踪。此迭代散射模型可准确模拟声波在周期结构界面上的散射情况，相较于传统的相干几何法在低频段有了较大的精度提高。此方法是对几何声学方法的有效补充，可为周期散射结构存在条件下的室内声场仿真提供简便且准确的新方法。