声扩散体在音乐厅等观演建筑中被广泛应用于改善声场均匀度、提供丰富的侧向反射声及避免声聚焦等。自Schroeder运用数论原理提出3种Schroder扩散体^[1-3]，最大长度序列(max length sequence，MLS)扩散体、二次剩余序列扩散体(quadratic residual diffuser，QRD)和本原根序列扩散体(primitive root diffuser，PRD)以来，对该类扩散体扩散性能的研究从未停止。如Cox^[4]提出的更便于生产的无翼板QRD扩散体，降低了扩散体的吸声且保持了较好的扩散性能。为减小Schroder扩散体尺寸并保持低频扩散性能，“折叠阱”及主动控制技术被运用于扩散体设计中^[5-6]。基于声超表面概念的Schroder扩散体，更进一步减小了扩散体尺寸，但相对于传统Schroder扩散体，该扩散体中高频扩散性能大幅降低^[7]。由于单个Schroder扩散体的尺度有限，在实际工程应用中，常将扩散体在界面上按一定规律排布安装，因此需考虑排布方式对整体扩散性能的影响。Angus^[8]首次将通信技术中的调制技术应用于Schroder扩散体的排布以改善其整体扩散性能。

基于光学Fraunhofer理论的扩散体反射声场计算公式存在局限性^[9-14]，难以对扩散体反射声场进行准确预测。随着计算技术的发展，将声场数值计算与寻优算法相结合可设计出形式多样、性能更加优越的声学扩散体。有效的声场计算方法、合理的扩散体性能评价指标以及高效的寻优算法是该技术的关键。基于波动声学的时域有限差分法，被应用于多种声场计算中^[15-17]。本文采用时域有限差分法计算扩散体的反射声场。在近似远场范围内，讨论了传声器、声源与扩散体间的距离对扩散系数计算结果的影响。免疫遗传算法具有全局寻优、鲁棒性的优点，文中将时域有限差分法与免疫遗传算法相结合，对六阶阶梯状扩散体形状进行优化，进而分析周期及非周期排列方式对排布后的扩散体整体扩散性能的影响，为该类扩散体的布置提供参考。

1 计算方法 1.1 扩散系数计算模型

依据ISO 17497-2标准^[18]，将扩散体置于全消声室中，37个接收点以5°间隔均匀布置在以扩散体为中心的半圆上，声源位于扩散体中轴线上，如图 1所示。

截取来自扩散体的反射声信号并对其进行1/3倍频程滤波可得到对应频带的声压级L_i，再由式(1)和(2)计算归一化扩散系数。

$ {d = \frac{{{{\left( {\sum\limits_{i = 1}^n 1 {0^{{L_i}/10}}} \right)}^2} - \sum\limits_{i = 1}^n {{{\left( {{{10}^{{L_i}/10}}} \right)}^2}} }}{{(n - 1)\sum\limits_{i = 1}^n {{{\left( {{{10}^{{L_i}/10}}} \right)}^2}} }}.} $

(1)

$ {{d_n} = \frac{{d - {d_r}}}{{1 - {d_r}}}.} $

(2)

式中n=37。d_n为归一化扩散系数，d_r为与扩散体具有相同尺寸的平板的扩散系数。

1.2 反射声场理论计算模型

Schroder基于光学Fraunhofer理论，依据Helmholtz-Kirchhoff积分方程，得到扩散体反射声场计算公式^[1]：

$ p \approx \int_x R (x){{\rm{e}}^{(2\pi {\rm{i}}/\lambda )x}}{\rm{d}}x. $

(3)

式中，p表示反射声压，R(x)表示局部反射因子，λ表示波长。该式为设计Schroder扩散体的理论基础，空间反射声场与局部反射因子R(x)之间存在Fourier变换关系。式(3)成立需满足如下条件：1) 扩散体阱中传播的是平面波；2) 扩散体表面为声波局部作用表面；3) 接收点需在远场；4) 采用Kirchhoff边界条件，将反射声压表示为入射声压与局部反射因子的乘积。

1.3 声场时域有限差分法数值计算模型

理想流体媒质中小振幅声波的连续性方程与运动方程分别为

$ {\frac{{\partial p}}{{\partial t}} + \kappa \nabla \cdot \mathit{\boldsymbol{u}} = 0,} $

(4)

$ {\nabla p + {\rho _0}\frac{{\partial \mathit{\boldsymbol{u}}}}{{\partial t}} = 0.} $

(5)

式中，ρ₀为媒质密度，u为质点振动速度，κ=ρ₀c²，c为声速。令$\boldsymbol{u}=u x \overrightarrow{\mathrm{i}}+u y \overrightarrow{\mathrm{j}}$，则上述方程的时域有限差分交错网格格式为^[17]

$ {ux_{i + \frac{1}{2},j}^{n + 1} = ux_{i + \frac{1}{2},j}^n - \frac{{\Delta t}}{{{\rho _0}}}\left( {\frac{{p_{i + 1}^{n + \frac{1}{2}} - p_{i,j}^{n + \frac{1}{2}}}}{{\Delta x}}} \right),} $

(6)

$ {uy_{i,j + \frac{1}{2}}^{n + 1} = uy_{i,j + \frac{1}{2}}^n - \frac{{\Delta t}}{{{\rho _0}}}\left( {\frac{{p_{i,j + 1}^{n + \frac{1}{2}} - p_{i,j}^{n + \frac{1}{2}}}}{{\Delta y}}} \right),} $

(7)

$ {p_{i,j}^{n + \frac{1}{2}} = p_{i,j}^{n - \frac{1}{2}} - \kappa \Delta t\left( {\frac{{ux_{i + \frac{1}{2},j}^n - ux_{i - \frac{1}{2},j}^n}}{{\Delta x}} + \frac{{uy_{i,j + \frac{1}{2}}^n - uy_{i,j - \frac{1}{2}}^n}}{{\Delta y}}} \right).} $

(8)

式中，Δt为时间步长，Δx、Δy为空间步长。依据Courant条件$ \Delta x / \Delta t=\sqrt{2} c$，本文中Δx=Δy=0.006 25 m。

2 模型计算结果及其与实测结果对比 2.1 反射声场理论计算结果与FDTD数值计算结果对比

分别采用理论公式与声场FDTD法对单周期排布的扩散体的散射声场进行计算。扩散体阱高为[0.2 m, 0 m, 0.2 m, 0.2 m, 0.2 m, 0 m, 0.2 m]，阱宽0.1 m，单周期扩散体的宽度为0.7 m，如图 2所示。计算时，计算点位于半径为50 m的半圆周上，声源位于扩散体中轴线上距扩散体100 m。

图 3给出了反射角为0°~90°的计算结果。从图中可以看出，频率为250 Hz时，理论计算与FDTD计算结果一致性较差。理论计算公式是建立在扩散体满足局部作用表面的假设下，频率越低，波长越长，声波与扩散体各部分之间相互作用越明显，越偏离局部作用表面的假设，导致2种计算结果差异较大。当频率升高至850 Hz时，波长变短，这种相互作用变弱，如图 3b所示，2种计算结果一致性较好。随着频率继续升高，计算点将不再满足远场条件，如图 3c和3d所示，2种计算结果差异增大。因此，理论计算公式仅在中频范围可对扩散体的反射声场准确预测。这与文[13-14]中BEM与理论计算对比得出的结论一致。

2.2 反射声场FDTD数值计算结果与实测结果对比

应用FDTD对一扩散体(见图 4a)的反射声场进行计算，并与全消声室内的测试结果进行对比。扩散体阱高依次为[0.075 m, 0 m, 0.05 m, 0.05 m, 0 m, 0.075 m]，阱宽为0.086 m，翼板宽与高分别为0.012 m与0.1 m。根据ISO 17497-2^[18]，计算点以5°为间隔均匀布置在以扩散体为中心、半径为2.5 m的半圆上，为保证有效测量范围，计算点布置在与扩散体中轴线夹角-60°~60°之间，如图 4b所示。声源位于扩散体中轴线上，距扩散体5 m。FDTD计算时，采用Gauss脉冲声源，利用完美匹配层^[17](perfectly matched layer，PML)边界模型作为截断边界条件模拟自由场，扩散体设置为刚性界面。实测时采用电火花声源。图 5a-5d给出了扩散体反射声场的FDTD计算结果与测量结果对比。从图中可以看出，在测试的频率范围内，FDTD计算结果与实测结果相吻合。FDTD法可有效计算扩散体反射声场。故下文将采用FDTD计算扩散体反射声场。

3 结果与分析 3.1 声源及测点与扩散体的距离对扩散系数的影响

理论上，远场测量要求传声器与扩散体之间距离数百米，满足该尺度的自由场在实际应用中较难实现^[14]。因此ISO 17497-2^[18]建议按近似远场条件进行测量，且将80%测点布置在镜像反射区间外。

为讨论声源及测点与扩散体的距离对扩散系数的影响，对一阱宽为0.062 5 m、阱高为[0.15 m, 0.2 m, 0 m, 0.1 m, 0.05 m, 0.25 m]的PRD扩散体(如图 6所示)，在4种测量距离下利用FDTD计算得到归一化扩散系数曲线，如图 7所示。图中r为测点与扩散体的距离，计算中分别取2、5、10和20 m。相应地，声源与扩散体的距离s分别为4、10、20和40 m。测量距离r=2 m，声源距扩散体s=4 m，已满足近似远场测量条件。从图 7可以看出，随着距离增大，扩散频率特性曲线在100~5 000 Hz的整个测量频率范围内一致性很好，说明在近似远场条件下可对扩散系数作出准确测量和计算。该PRD扩散体的设计频率约为700 Hz，图 7显示在700及其倍频率所在的1/3倍频带上出现扩散系数极大值，与理论扩散频率吻合。

图 8给出了4种测量距离条件下，37个测点在2 000 Hz的反射声压级。图中显示4条声压级曲线变化趋势一致。声源与扩散体之间的距离s加倍，计算点半径r增加一倍时，在各个测点上声压级降低固定值。此时可将扩散体视为能够向半空间各方向辐射不同声能量的特殊点声源。这也解释了图 7中4种测量距离下扩散系数曲线一致的原因。

3.2 单周期扩散体的优化

采用文[19]中的方法，将FDTD与免疫遗传算法相结合，对图 6所示的单周期六阶PRD扩散体阱高进行优化。优化目标函数见式(9)。d和σ分别表示100~5 000 Hz范围内1/3倍频带归一化扩散系数的平均值和标准差。目标函数值越低，则扩散体的扩散性能越好。

$ f = \frac{1}{{\bar d - 0.5\sigma + 1}}. $

(9)

计算时，计算点与扩散体的距离r=5 m，声源距离扩散体s=10 m。优化后的扩散体记为A，其阱高依次为[0.237 5 m, 0.212 5 m, 0.187 5 m, 0.212 5 m, 0.237 5 m, 0.075 m]，如图 9所示。图 10显示单周期扩散体A与单周期PRD扩散体的扩散系数频率特性曲线对比。从图中可以看出，在200~1 000 Hz与3 000~5 000 Hz，扩散体A扩散系数高于PRD，在低频范围优化效果尤其明显，在300 Hz提高值达到0.6，在1 000~3 000 Hz略低于PRD的扩散系数。由此也说明采用式(9)作为免疫遗传算法的目标函数，引入标准差作为惩罚因子，可较好地起到平滑扩散系数频率特性曲线的作用^[19]。优化后的扩散体，其扩散系数频率特性曲线较PRD平直。

3.3 扩散体多周期重复排布及其对扩散性能的影响

将A扩散体与PRD扩散体分别重复扩展6个周期后，扩散体总宽为2.25 m。对扩散体的反射声场计算时，在满足近似远场条件下，令计算半径r=10 m，声源距扩散体s=20 m。在下文的讨论中r、s均采用该值。

图 11-12分别给出了单周期与6个周期重复排布的A扩散体与PRD扩散体的归一化扩散系数频率特性曲线。从图中可以看出，2种扩散体重复周期扩展后，各频带扩散系数均有降低，这可由图 13解释。如果相叠加的扩散体表面反射声场声压相位差接近2π，那么叠加后声压有极大值。假设为图中的C位置，由于周期扩展的扩散体阱深变化具有空间周期性，因此扩散体表面反射声场的声压极大值在空间中周期出现。这就形成了扩散体表面反射声场的能量瓣，反射声场的能量分布不均匀性增加，从而使扩散体重复周期排布后的整体扩散性能降低。

图 14显示单周期与6周期重复排布的扩散体A在37个测点处的850 Hz反射声压级。从中可以看出，单个A扩散体能够较好地将声能向空间扩散，各测点声压级差较小。而6周期重复排布后的A扩散体以“瓣”的形式向空间反射声能，其中3个主能量瓣集中在测点“4” “19”和“33”位置处，主能量瓣之间有多个次能量瓣，瓣与瓣之间存在极小值，导致扩散性能变差。

3.4 同一扩散体正反相非周期调制排布及其对扩散性能的影响

以6周期排布为例，采用二进制序列对扩散体进行非周期调制排布并分析其整体扩散性能。调制码选择6位二进制码，从000000至111111。令二进制码中“1”布置扩散体，“0”布置扩散体的反相扩散体。扩散体A的反相扩散体，其阱高为[0.012 5 m, 0.037 5 m, 0.062 5 m, 0.037 5 m, 0.012 5 m, 0.175 m]。PRD扩散体的反相扩散体阱高为[0.1 m, 0.05 m, 0.25 m, 0.15 m, 0.2 m, 0 m]。对A和PRD扩散体，依据式(9)，分别计算2⁶种调制方式下的f值，f最小值对应最优的调制方式。其中，111111和000000这2种调制模式分别相当于扩散体及其反相扩散体的重复周期排布。经计算，PRD及A扩散体的最优非周期调制排布方式分别为101001与011001。非周期调制排布与重复扩展周期排布的PRD扩散体和A扩散体的归一化扩散系数频率特性曲线对比如图 15所示。非周期调制后的扩散体整体扩散性能提高，其中低频扩散性能改善明显。

PRD扩散体与扩散体A经非周期调制排布后的扩散性能对比如图 16所示。从图中可以看出，A扩散体非周期调制排布后的整体扩散性能优于PRD，但两者的扩散系数频率特性曲线在1 000 Hz附近都存在一个极小值。PRD扩散体与A扩散体分别采用理论式(3)与免疫遗传算法获得，两者设计理论不同。将扩散体与其反相扩散体按照调制码排布时，扩散体阱深之间形成的相位差与其反相扩散体阱深之间形成的相位差为正负相位，PRD扩散体与其反相扩散体扩散性能相似，在PRD扩散体的设计频率处，图 13中C位置不再出现声压极大值，将扩散体A与其反相扩散体按照调制码排布没有类似的性质。因此对PRD扩散体的排布进行调制仅改善了其设计频率附近的扩散性能，而A的反相扩散体扩散性能影响了调制效果，故采用扩散体与其反向扩散体进行非周期调制没有获得一条平直的扩散系数频率特性曲线。

3.5 不同扩散体非周期调制排布及其对扩散性能的影响

为进一步改善整体扩散体的扩散频率特性，采用6位二进制序列对A扩散体与PRD扩散体进行混合调制，令“1”布置A扩散体，“0”布置PRD扩散体。依据式(9)，计算2⁶种调制方式下的f值，f最小值对应最优的调制方式，其对应的序列为001101。经此方式调制后的整体扩散体(记为：PRD0-A1-001101)的归一化扩散系数如图 17所示。图中同时给出了扩散体A经正反相非周期调制排布后(A-011001)的归一化扩散系数。可以看出，PRD0-A1-001101扩散体在1 000 Hz频率附近的扩散性能显著提升，在300~5 000 Hz频率范围内获得平直的扩散频率特性曲线。从图 15a也可以看出，虽正反相非周期调制排布比周期重复排布能获得更优良的整体扩散性能，但所获得的扩散系数频率特性曲线变化趋势较一致，扩散频率范围并没有更好地扩展。不同扩散体非周期混合调制方式中，2种扩散体的阱高不存在确定关系，可拓宽有效扩散频率。图 18显示800 Hz附近时2种调制方式的反射声压级对比，可以看出不同扩散体非周期混合调制排布下，反射声场中各点声压级差变小，声场扩散更加均匀。

4 结论

本文采用时域有限差分法计算扩散体的反射声场。在近似远场条件下，声源、测点与扩散体之间的距离对扩散系数频率特性曲线没有影响，且与实测结果相吻合。将时域有限差分法与免疫遗传算法相结合，对单周期六阶阶梯状扩散体的形体进行优化。优化后的单周期扩散体的扩散性能优于原扩散体，但经周期排布后，两者的扩散性能均出现降低。

为探求扩散体的最优排布方式，将扩散体及其反向扩散体分别看成二进制序列中的1和0，对于连续排布的n个扩散体，计算2ⁿ种调制方式下的目标函数值，寻求最优排布方式，其整体扩散性能优于重复周期排布扩散体，但这种正反相扩散体的调制方式不能得到平直的扩散系数频率特性曲线。文中进而将优化前后的扩散体分别看成二进制序列中的0和1，将2种扩散体按上述二进制序列调制寻优后得到的混合排布方式能获得更好的整体扩散性能。与传统的扩散体多周期重复排布相比，混合排布方式既可提高整体扩散性能，又可提高设计的灵活性和空间变化，从而满足不同空间设计的需要。在扩散体优化设计中，因排布方式对扩散性能有显著影响，因此仅考虑单个扩散体的扩散性能是不够的，还需考虑整体扩散性能。随着计算技术的发展，将声场时域有限差分计算模型与寻优算法相结合，对扩散体的界面进行整体扩散性能优化，寻求满足一定声场条件下的扩散边界形状，将是扩散体优化设计的重要方向。