本文仅做简单理论探讨，学艺不精，如有错误，恳请各位老师同仁斧正。
一、用几何学方法建模声音
在声学领域，当声波的波长显著小于房间的几何特征尺寸、表面凹凸程度及整体空间尺度时，声能的传播表现出类似于光射线的直线特性 。这种物理近似被称为几何声学（Geometrical Acoustics, GA），它通过将复杂的声学波动方程简化为沿直线传播的能量射线，从而在计算资源可控的前提下，实现了对复杂建筑空间声学特性的高精度预测 。上一期我们介绍了最简单的几何学建模方法——ISM镜像声源法，镜像声源法虽然直接，但需要超大运算量以支撑完整的混响尾音，这显然十分不经济。为了解决这一问题，Pyroomacoustics库中提供了Hybird方法，利用射线追踪法进行完整的混响尾音计算。
1.什么是射线追踪法？
射线追踪法（Ray Tracing Method, RT）的理论也相当简单，此方法将声波简化为沿直线传播、携带特定份额能量的“声射线” 。在传播过程中，射线会根据表面的物理属性进行镜面反射、扩散反射或被吸收，能量随之衰减 。最终，系统利用具有一定体积的接收器模型（如探测球）收集这些射线携带的时间和能量信息，生成能量直方图，进而合成房间脉冲响应（RIR）并导出 RT60、清晰度等声学评价指标 。
2.射线追踪法的开始：选取射线分布
射线追踪法首先需要选取需要追踪的射线，为了真实模拟全向声源或具有特定指向性的声源，不能完全依赖随机选取射线，必须在初始时刻确定射线的空间分布规律 。
在计算实现中，最基本的方法是利用蒙特卡罗（Monte Carlo）方法在单位球面上生成均匀分布的初始射线方向 。为了确保空间采样的公平性，避免在球极地区产生射线堆积，算法通常采用准随机序列（如索博尔序列）或等面积细分算法（如二十面体细分）来确定发射角 。每条射线初始携带的能量 E_0 被定义为声源总功率 P_w 与总射线数 N 的商 ：
$$E_0 = \frac{P_w}{N}$$
射线数量的选择直接关系到模拟的信噪比和收敛精度。在大型建筑空间（如音乐厅或剧院）的模拟中，通常需要发射数百万甚至数亿条射线，以确保空间中的接收器体积能够捕获到足够的采样样本 。
3.射线的传播路径与几何衰减
声射线在均匀介质中沿直线传播，其位置 $$P(t)$$ 随时间 $$t$$ 的变化遵循运动学方程
$$P(t) = P_0 + \vec{d} \cdot c \cdot t$$
其中 $ 为声速 。在传播过程中，射线主要经历两种形式的能量耗散。首先是介质吸收，即空气对声能的粘滞耗散和分子弛豫吸收，这在长距离传播和高频段尤为显著 。其次是几何衰减，但在射线追踪模型中，我们仅对每条射线而非声源进行追踪，这使得我们不需要计算平方反比定律 。
4.射线在边界上的反射/散射/吸收
当声射线与房间边界（如墙体、地板或家具表面）发生碰撞时，其能量状态和传播方向将依据界面的声学物理属性发生转变。界面的吸收系数 α 决定了声能损耗。射线碰撞后的剩余能量为
$$E_{out} = E_{in} \cdot (1 - \alpha)$$
在最理想的情况下，射线遵循斯涅尔定律（Snell's Law）进行镜面反射，即反射角等于入射角。这种处理方式适用于大面积、平整且刚性的表面 。
斯涅尔定律（Snell's Law）是描述波（如声波或光波）在两种介质界面处传播方向变化的物理定律。其实就是我们初中就学过的光的折射和反射。
然而，现实中的建筑表面往往具有粗糙度或微观几何结构，这导致声波在反射时不仅产生镜面分量，还会产生向各个方向散开的扩散分量。声学建模中引入散射系数 $s$ 来量化这一现象 。
目前主流的射线追踪引擎采用基于权重的随机决策机制 。在射线与表面交互时，系统生成一个随机数。若随机数小于 s，则射线执行扩散反射；否则执行镜面反射 。扩散反射的方向通常遵循朗伯余弦定律（Lambert's Cosine Law），其反射强度与表面法线夹角的余弦值成正比，从而模拟出声能向半球空间均匀散布的效果 。
朗伯余弦定律（Lambert's Cosine Law）是描述理想扩散表面（漫反射面）反射特性的基本物理定律 。该定律指出，从理想扩散表面反射出的声能，其在某一方向上的强度与该方向与表面法线（垂直线）之间夹角 $\theta$ 的余弦值成正比 。
为了更精确地模拟具有特定结构（如一维扩散器）的表面，一些先进算法引入了定向散射模型。该模型不再简单地将反射方向均匀化，而是根据表面的拓扑特征（如槽深或排列方向）调整散射概率包络，使得散射能量在特定平面内具有更高的密度 。
5.如何接收计算出的射线
想必读到此处，各位一定会有一个疑问：如果我们将接收器也定义为一个点，除非我们采用相当数量的射线，否则射线恰好抵达接收器所在点的概率微乎其微，根本不可能实现我们需求的完整混响。为此，我们将接收器定义为一个体，凡是经过这个体的射线最终都被我们合并到RIR的计算当中去。
这也带来一个问题：接收器通常被定义为一个具有特定半径 r 的探测球（Detection Sphere）。半径的选择存在精度与效率的平衡：如果半径过大，会引入空间平均误差，导致时间分辨率模糊；如果半径过小，则需要海量的射线数才能获得统计稳定的结果 。在符合人耳听觉特性的模拟中，半径通常设定为 0.1 米至 0.2 米，大致对应人头的物理尺寸 。
当射线进入探测球范围时，相邻时间内射线携带的能量、频率及时间延迟信息被记录并累加到对应的频带直方图中 。能量直方图是射线追踪结果的原始呈现形式，通常以倍频程（Octave Bands）为单位划分，涵盖 63Hz 到 8kHz（甚至 16kHz）的听觉范围 。
直方图的时间槽（Bins）宽度决定了计算的分辨率。为了获得高质量的脉冲响应，时间槽通常设定在毫秒级甚至微秒级。直方图中的每一个点代表了在该时间窗口内到达的所有射线能量之和，反映了房间的能量衰减特性，即大致包络 ，但缺乏声波所需的精度 。
6.如何计算最终的RIR响应
射线追踪法的直接产物是能量包络，而不是可听的声压信号。将能量直方图转化为具有正确相位和精细结构的房间脉冲响应（Room Impulse Response, RIR）是射线追踪建模的核心难点之一 。
由于直方图记录的是能量 I，而 RIR 是声压 P 随时间的变化，合成的第一步是提取声压包络 。通常采用开方计算，即声压包络正比于能量的平方根：
$$P \propto \sqrt{I}$$
为了将包络转换为可听的声音，我们需要引入一个泊松随机过程（Poisson Process），首先，我们根据房间体积 V 和声速 c，计算单位时间内反射脉冲的理论密度：
$$\mu(t) = \min\left(10000, \frac{4\pi c^3 t^2}{V}\right)$$
然后根据该密度产生一系列随机分布的时间点（狄拉克脉冲），模拟反射声随时间推移而呈平方级逐渐变密的物理现象 。
此时，我们又遇到了另一个问题：在几何声学近似下，射线的相位信息已经丢失。为了模拟真实房间中扩散声场的特性，合成程序会为上述生成的每个随机脉冲赋予一个在 $$[-\pi, \pi]$$ 范围内均匀分布的随机相位 。这种随机化处理能确保合成的脉冲响应在听觉上具有自然的、无人工痕迹的“空气感”，避免产生由于周期性导致的颤动回声 。
这些计算通常分倍频程进行，现在，我们得到了以倍频程为基础的脉冲响应，此时根据射线追踪得到的各频带能量直方图，对上述随机脉冲序列进行权重调节 。然后每个频带的信号通过一组特定的滤波器。为了防止信号在时间轴上产生错位（即群延迟问题），必须使用零相位滤波器（Zero-phase filters），以保证不同频率的脉冲能够严丝合缝地在时间轴上对齐 。
最后，将所有频带合成的声压信号进行累加，得到一个宽频带的房间脉冲响应（RIR） 。
7.如何使用Pyroomacoustics库进行ISM和RT混合计算
为了利用ISM的确定性以及RT的计算简便性，现代应用场景中通常以时域结合两种方法，为了保证 RIR 在过渡点不出现明显的能量跳变，Pyroomacoustics内部会计算射线追踪生成的能量响应，并将其与 ISM 产生的离散脉冲进行匹配。此时我们拥有一个经验公式（Perceptual Mixing Time）：
$$t_{mix} \approx \sqrt{V} \text{ (ms)}$$
其中V是房间体积。
下面是一段简单的代码，展示了如何使用Pyroomacoutics进行Hybird模拟
import pyroomacoustics as pra
import numpy as np
import matplotlib.pyplot as plt
# 1. 基础配置
fs = 48000 # 采样率
room_dim = [7, 5, 3] # 房间尺寸 (长, 宽, 高) 单位：米
# 2. 设置吸声系数和混合模型参数
# 我们定义一个材料，包含吸声系数和散射系数
material = pra.Material(energy_absorption=0.2, scattering=0.1)
# 创建房间
# max_order: ISM 的最大阶数（处理早期反射）
# ray_tracing: 设置为 True 以启用射线追踪（处理晚期混响）
room = pra.ShoeBox(
room_dim,
fs=fs,
materials=material,
max_order=3, # ISM 阶数，通常 3-5 阶足够
ray_tracing=True # 开启射线追踪
)
# 3. 配置射线追踪引擎的具体参数
# 这些参数决定了晚期混响的质量
room.set_ray_tracing(
receiver_radius=0.5, # 接收器（麦克风）的虚拟半径
n_rays=10000, # 发射射线的数量，越多越精确但越慢
energy_thres=1e-7 # 能量阈值，低于此值停止追踪
)
# 4. 添加声源和麦克风
source_pos = [1, 1, 1.5]
mic_pos = np.array([[5, 3, 1.5]]).T # 必须是 (维度, 数量) 的形状
room.add_source(source_pos)
room.add_microphone(mic_pos)
# 5. 计算 RIR
# 该步骤会自动结合 ISM 和 Ray Tracing
room.compute_rir()
# 6. 可视化结果
rir = room.rir[0][0] # 获取第一个声源到第一个麦克风的 RIR
plt.figure(figsize=(10, 4))
plt.plot(np.arange(len(rir)) / fs, rir)
plt.title("Hybrid RIR (ISM + Ray Tracing)")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()
8.总结
至此，我们已经介绍完毕Pyroomacoustics中支持的，也是行业内最常用的两种几何声学模拟方法：ISM和RT。通过两种方法的结合，我们得以较为完整的对空间进行建模——仅限高频。这种基于高频近似的建模方法存在其天然的适用边界，即所谓的施罗德频率（Schroeder Frequency） 。施罗德频率：
$$f_s \approx 2000 \sqrt{\frac{T_{60}}{V}}$$
施罗德频率划定了房间声场从低频简正模式（Normal Modes）向高频扩散声场（Diffuse Field）转化的临界区域。在施罗德频率以下，声波的波动特性（如驻波、干涉和显著的衍射）占据主导地位，此时射线追踪法的准确性会大幅下降，因为离散的驻波会导致能量衰减出现非线性跳跃 。通常需要结合有限元法（FEM）或边界元法（BEM）等波动学方法进行杂交建模 。
我们将在之后的文章中介绍波动学方法的空间建模。
发表回复