基于格栅效应的室内周期矩形声扩散体散射系数预测方法

摘要

本发明公开了一种基于格栅效应的室内周期矩形声扩散体散射系数预测方法，用于解决现有方法效率低而导致耗时长的技术问题。技术方案是将三维周期矩形扩散体进行降维处理，以二维周期矩形扩散体轮廓为计算模型，根据格栅方程预测在某一入射角度时，入射波发生的散射阶次，得到这些阶次散射波的传播方向；根据散射波的传播方向对进入矩形声扩散体凹槽内的声波进行路径跟踪建模，依据建模结果对每个阶次散射波的能量进行预测，最终根据各个反射波的能量得到方向入射散射系数及无规入射散射系数的数值。本发明由于无需对周期扩散体模型进行离散处理，仅对二维模型展开计算，因此可提高周期矩形扩散体散射系数的预测效率，适用于工程应用中的快速估计。

说明书

技术领域

本发明涉及一种室内周期矩形声扩散体散射系数预测方法，特别是涉及一种基于格栅效应的室内周期矩形声扩散体散射系数预测方法。

背景技术

室内周期矩形声扩散体是一种特殊的界面声散射结构，它是指拥有矩形几何形状的材料在宏观空间上按照周期规则排列的结构。相较于其他形状的扩散体，矩形扩散体由于设计简单、制造方便等优势在音乐厅、录音室、剧院等对声场扩散要求较高的场所具有十分广泛的应用。室内周期矩形声扩散体应用的关键是对其散射系数进行预测，声学散射系数的定义为通过扩散体反射后，非镜面反射声能与全部反射声能的比值，其获取通常有两种方式，一种是测量实验方法，一种是数值计算方法。对于测量实验，直到2004年才建立国际标准，即混响室转台法，该方法还有诸多需要完善的方面，而且面对种类繁多的声扩散体，实验测量需要消耗大量的人力物力；另外，对于处于设计阶段的声扩散体是无法进行实验测量的，因此研究散射系数的数值计算方法是非常必要的。但是由于理论上的限制，散射模拟中的镜面反射声能很难通过计算得到，因此数值计算的发展一直较为缓慢。2000年E.Mommertz提出计算镜面反射声能时可以用一个与周期结构尺寸相同的纯平的结构代替周期结构，由此发展了散射系数的数值算法。此后，对散射系数的数值计算获得了快速发展以及越来越多的关注。

文献“授权公告号是CN102938017B的中国发明专利”公开了一种基于无网格模型计算周期结构板声学散射系数的方法。该方法将三维的周期结构扩散体作为计算模型，通过无网格数值算法将周期扩散体用节点进行离散，在计算中将三维周期扩散体划分为节点模型，同时建立一个参考板节点模型，此参考板纯平，与周期结构板整体尺寸、形状相同，然后利用移动最小二乘法获得形函数，并离散Helmholtz微分方程，求解可计算节点声压差的系统方程，从而得到计算周期扩散体及其参考板上节点处的声压差，进一步得到接收点的声压，最终计算得到某一入射方向的方向散射系数及无规入射散射系数。此发明专利方法存在的一个突出的问题就是计算效率十分低下，当扩散体面积较大时，离散后的单元、网格、节点数量会十分巨大，使得计算过程中的各项矩阵规模臃肿，从而需要耗费大量的计算时间，显然无法满足实际工程应用中对时效性的要求，因此有必要发展一种能够快速预测扩散体散射系数的方法。周期型的室内矩形声扩散体具有一个重要的性质就是格栅效应，它可以将入射声波反射向若干特定方向，这些方向可根据扩散体的几何参数及入射波参数直接求得，而无需繁杂的矩阵插值迭代运算。基于此格栅效应，可建立一种周期矩形声扩散体散射系数的快速预测方法。

发明内容

为了克服现有方法效率低而导致耗时长的不足，本发明提供一种基于格栅效应的室内周期矩形声扩散体散射系数预测方法。该方法将三维周期矩形扩散体进行降维处理，以二维周期矩形扩散体轮廓为计算模型，根据格栅方程预测在某一入射角度时，入射波发生的散射阶次，得到这些阶次散射波的传播方向；根据散射波的传播方向对进入矩形声扩散体凹槽内的声波进行路径跟踪建模，依据建模结果对每个阶次散射波的能量进行预测，最终根据各个反射波的能量得到方向入射散射系数及无规入射散射系数的数值。本发明由于无需对周期扩散体模型进行离散处理，仅对二维模型展开计算，因此可提高周期矩形扩散体散射系数的预测效率，适用于工程应用中的快速估计。

本发明解决其技术问题所采用的技术方案是：一种基于格栅效应的室内周期矩形声扩散体散射系数预测方法，其特点是采用以下步骤：

步骤一、将三维周期矩形扩散体进行降维处理，以二维周期矩形扩散体轮廓为计算模型，获得各阶次散射波的反射方向。假设有一束频率为f，入射角为θ_i的声波入射到周期矩形扩散体上，经过扩散体的反射之后，反射波的散射角度根据格栅方程求得：

$>{\mathrm{sin\theta }}_s={\mathrm{sin\theta }}_i+\frac{n\lambda }{L},n=0,\pm 1,\pm 2,...\left(1\right)$>

式中，n代表反射波的阶次，θ_i为入射角，即入射波与垂直方向之间的夹角，θ_s为对应于反射阶次n的反射波的散射角，即散射波与垂直方向之间的夹角，L为周期扩散体一个周期的长度，波长λ＝c₀/f，c₀为空气中的声速。

步骤二、假设入射波一部分被周期矩形扩散体的上表面反射，记为声线1，一部分被凹槽内部表面反射，记为声线2，那么声线1与声线2之间的声程差为：

$>dis=c_{0}t=H(\frac{1}{{\mathrm{cos\theta }}_i}+\frac{1}{{\mathrm{cos\theta }}_s})---\left(2\right)$>

式中，t为声线2在凹槽内传播所用的时间，H为矩形扩散体的凹槽高度。

根据两个声线的声程差，得到两声线经过扩散体反射之后各自的相位差为：

$>\left(\begin{array}{c}{\phi }_1=0\\ {\phi }_2=\frac{2\pi }{\lambda }H(\frac{1}{{\mathrm{cos\theta }}_i}+\frac{1}{{\mathrm{cos\theta }}_s})\end{array}\right)---\left(3\right)$>

根据相位差，求得周期矩形扩散体的反射参数为：

式中，x表示二维周期扩散体的横向位置，m表示不同位置处的周期，l表示凹槽的长度。

由于扩散体的类型为周期形式，那么反射参数同样为周期函数，其傅里叶级数为：

$>R\left(x\right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{\alpha }_n{e}^{jn\frac{2\pi }{L}x}---\left(5\right)$>

式中，α_n为傅里叶系数，由下式求得：

$>{\alpha }_n=\frac{1}{L}{\int }_0^{L}R\left(x\right)e^{-jn\frac{2\pi }{L}x}dx---\left(6\right)$>

对于一个入射声波：

$>p_i(x,y)=e^{-j\frac{2\pi }{\lambda }({\mathrm{sin\theta }}_{i}x-{\mathrm{cos\theta }}_{i}y)}---\left(7\right)$>

经过周期矩形声扩散体的作用后，反射波的声压为：

$>p_r(x,0)=p_i(x,0)·R\left(x\right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{\alpha }_n{e}^{-j\frac{2\pi }{\lambda }({\mathrm{sin\theta }}_i-\frac{n\lambda }{L})x}---\left(8\right)$>

式(8)表明α_n与反射波的声压幅值相关，因此只要确定了各个阶次反射波声压的α_n，即确定了反射波的声压。根据式(4)、式(6)，得到α_n的解为：

$>\left(\begin{array}{c}{\alpha }_0={\mathrm{\tau e}}^{{\mathrm{j\phi }}_1}+(1-\tau )e^{{\mathrm{j\phi }}_2}\\ {\alpha }_n=\frac{j}{2n\pi }(e^{{\mathrm{j\phi }}_1}-e^{{\mathrm{j\phi }}_2})(e^{-j2n\pi \tau }-1),n\ne 0\end{array}\right)---\left(9\right)$>

式中τ＝l/L为占空比，表示凹槽长度与单个周期长度的比值。获得表示声压幅值的α_n之后，将其作平方处理即得到反射波的能量：

$>\left(\begin{array}{c}{E}_0=|{\alpha }_0{|}^2=1-2\tau (1-\tau )(1-cos\Delta \phi )\\ E_n=|{\alpha }_n{|}^2=\frac{1}{n^2{\pi }^2}(1-cos2n\pi \tau )(1-cos\Delta \phi ),n\ne 0\end{array}\right)---\left(10\right)$>

式中$>\Delta \phi ={\phi }_2-{\phi }_1={\phi }_2=\frac{2\pi }{\lambda }H(\frac{1}{{\mathrm{cos\theta }}_i}+\frac{1}{{\mathrm{cos\theta }}_s}).$>

步骤三、在获得所有阶次反射波的能量之后，获得某个声源入射方向的散射系数：

$>s_d=\frac{\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{E}_n}{E_0+\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{E}_n}---\left(11\right)$>

步骤四、将平面外的入射声波转换到平面内，将其投影作为新的入射波，同时将频率进行偏移处理，即：

式中，θ_i'为入射声波在平面内的投影与垂直方向的夹角，为入射波与平面的夹角。

在求得所有入射方向上的方向入射散射系数之后，求得无规入射散射系数：

式中，为设置在周期扩散体周围半球面上声源点的入射方向，θ为方位角为俯仰角。

本发明的有益效果是：该方法将三维周期矩形扩散体进行降维处理，以二维周期矩形扩散体轮廓为计算模型，根据格栅方程预测在某一入射角度时，入射波发生的散射阶次，得到这些阶次散射波的传播方向；根据散射波的传播方向对进入矩形声扩散体凹槽内的声波进行路径跟踪建模，依据建模结果对每个阶次散射波的能量进行预测，最终根据各个反射波的能量得到方向入射散射系数及无规入射散射系数的数值。本发明由于无需对周期扩散体模型进行离散处理，仅对二维模型展开计算，因此可提高周期矩形扩散体散射系数的预测效率，适用于工程应用中的快速估计。经测试，从200Hz到4000Hz，1/3倍频程上的无规入射散射系数，通过与测量实验相比较发现两者之间非常接近，平均差值仅为0.014，变化趋势也一致，均是随着频率的增加，散射系数随之增大。由于本发明方法仅需要利用周期矩形扩散体的几何参数进行计算，完成上述频段范围内散射系数的计算所耗费的时间仅为20s。

下面结合附图和具体实施方式对本发明作详细说明。

附图说明

图1是本发明方法中周期矩形声扩散体尺寸及入射波的示意图。

图2是本发明方法中入射波转换示意图。

具体实施方式

参照图1-2。本发明基于格栅效应的室内周期矩形声扩散体散射系数预测方法具体步骤如下：

以单个周期长度为L＝0.2m，凹槽长度为l＝0.1m，凹槽高度为H＝0.02m的周期矩形声扩散体为例。

步骤1：将三维周期矩形扩散体进行降维处理，以二维周期矩形扩散体轮廓为计算模型，获得各阶次散射波的反射方向。假设在二维周期扩散体所在平面内有一束频率为f＝2000Hz，入射角θ_i为30°的声波入射到扩散体上，根据格栅方程：

$>{\mathrm{sin\theta }}_s={\mathrm{sin\theta }}_i+\frac{n\lambda }{L},n=0,\pm 1,\pm 2,...\left(1\right)$>

求得散射波共有2个阶次，分别为n＝0阶、n＝-1阶，所对应的散射角分别为30°及20.5°。

步骤2：求解各阶次散射波的能量。上述入射波一部分被周期矩形扩散体的上表面反射，记为声线1，一部分被凹槽内部表面反射，记为声线2。两个声线之间的声程差为：

$>dis=c_{0}t=H(\frac{1}{{\mathrm{cos\theta }}_i}+\frac{1}{{\mathrm{cos\theta }}_s})---\left(2\right)$>

求得声线1与声线2之间的声程差为：0阶散射波dis＝0.0462m，-1阶散射波dis＝0.0444m。根据声程差，求得两声线经过扩散体反射之后各自的相位差为：

根据相位差求得周期扩散体的反射系数为：

由于扩散体的类型为周期形式，那么反射参数同样为周期函数，其傅里叶级数为：

$>R\left(x\right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{\alpha }_n{e}^{jn\frac{2\pi }{L}x}---\left(14\right)$>

傅里叶系数α_n为：

$>{\alpha }_n=\frac{1}{L}{\int }_0^{L}R\left(x\right)e^{-jn\frac{2\pi }{L}x}dx---\left(15\right)$>

对于步骤1中的入射声波：

$>p_i(x,y)=e^{-j\frac{2\pi }{\lambda }({\mathrm{sin\theta }}_{i}x-{\mathrm{cos\theta }}_{i}y)}---\left(16\right)$>

经过周期矩形声扩散体的作用后，其反射波的声压为：

$>p_r(x,0)=p_i(x,0)·R\left(x\right)=\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{\alpha }_n{e}^{-j\frac{2\pi }{\lambda }({\mathrm{sin\theta }}_i-\frac{n\lambda }{L})x}---\left(17\right)$>

式(8)表明α_n与反射波的声压幅值相关，因此只要确定了各个阶次反射波声压的α_n，即确定了反射波的声压。根据式(4)、式(6)，得到α_n的解为：

$>\left(\begin{array}{c}{\alpha }_0={\mathrm{\tau e}}^{{\mathrm{j\phi }}_1}+(1-\tau )e^{{\mathrm{j\phi }}_2}\\ {\alpha }_n=\frac{j}{2n\pi }(e^{{\mathrm{j\phi }}_1}-e^{{\mathrm{j\phi }}_2})(e^{-j2n\pi \tau }-1),n\ne 0\end{array}\right)---\left(18\right)$>

对α₀及α_n进行平方处理，求得0阶及-1阶的散射波能量为：

$>\left(\begin{array}{c}{E}_0=|{\alpha }_0{|}^2=1-2\tau (1-\tau )(1-\cos \Delta \phi )=0.4319\\ E_{-1}=|{\alpha }_{-1}{|}^2=\frac{1}{n^2{\pi }^2}(1-\cos 2n\pi \tau )(1-\cos \Delta \phi )=0.2169\end{array}\right)---\left(10\right)$>

步骤3：求解方向入射散射系数。在获得所有阶次反射波的能量之后，根据散射系数的定义获得步骤1中入射波所在入射方向的散射系数：

$>s_d=\frac{\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{E}_n}{E_0+\underset{n=-\infty }{\overset{+\infty }{\Sigma }}{E}_n}=0.3343---\left(11\right)$>

步骤4：求解无规入射散射系数。以周期矩形声扩散体中心为球心的一个半球面上设置声源点，此半球面的半径为10m。声源点的布置可根据经纬线交点的形式排列，每一个声源点代表这一个入射方向。由于步骤3中求得的散射系数所对应的入射声波处于二维周期矩形扩散体所在的xoy平面内，为了求得二维周期扩散体所在xoy平面以外的入射方向上的散射系数(参照图2)，需要对平面外的入射声波转换到平面内，将其投影作为新的入射波，同时将频率进行偏移处理，即：

在求得所有入射方向上的方向入射散射系数之后，根据下式求解无规入射散射系数，即工程中常用的散射系数。

式中，为设置在周期扩散体周围半球面上声源点的入射方向，θ为方位角为俯仰角。通过此式所求得的散射系数代表所有入射方向散射系数的一种平均结果，普遍应用于实际工程应用中。

在本实例中，计算了此从200Hz到4000Hz，1/3倍频程上的无规入射散射系数，通过与测量实验相比较发现两者之间非常接近，平均差值仅为0.014，变化趋势也一致，均是随着频率的增加，散射系数随之增大。有效地证明了本发明的有效性。在效率方面，由于本发明所涉及的方法仅需要利用周期矩形扩散体的几何参数进行计算，完成上述频段范围内散射系数的计算所耗费的时间仅为20s，证明了本发明具有极高的时效性。