附录 A2 - 介质的衰减与单次散射

光穿过吸收层时，一部分到达地面，其余被吸收。本例将均匀纯吸收介质、球形叶倾角分布的 Turbid 黑叶片，以及随机放置的离散水平黑叶片放在一起，比较三种表示的透射和吸收；再让均匀介质具有散射，检验不同方向的单次散射辐亮度。

光子发射位置和成像像元位置采用随机化 Sobol 低差异采样；介质内碰撞和后续散射采用独立随机采样。

场景设置

上图展示其中一个离散叶片场景，可拖动旋转，也可以单独打开。线框表示吸收层的边界；均匀介质和 Turbid 冠层使用相同的区域。

条件	设置
水平单元	2 m × 2 m，水平周期边界，最多穿越 10000 次
吸收层	Z = 0.1–1.1 m，厚度 H = 1 m
地表	Z = 0，完全吸收的黑色地面
光照	550 nm 平行太阳光，法向辐照度 1 W m⁻² nm⁻¹，无大气和天空光
太阳方向	天顶角 0°、45°、75°，方位角 90°
连续介质	吸收系数 k = 0、0.1、1、5 m⁻¹，散射系数为 0，关闭热发射
Turbid 冠层	球形叶倾角分布，叶片反射率与透射率均为 0，LAI = 0、0.2、2、10
离散叶片	100 片水平黑色正方形叶片，边长 0.2 m，LAI = 1；位置在区域内独立均匀抽样
场景重复	离散叶片使用 12 组位置，随机种子为 101–112；三个后端共用这些位置
光子数	介质和 Turbid 各用 100 万、200 万、400 万；离散场景各用 100 万和 400 万；最终均采用 400 万

跨越单元边缘的叶片同时放置相应的周期副本，使两侧叶片连续。每个周期单元的实际叶面积仍为 4 m²。

太阳光倾斜时，水平面接收的辐照度为 \(E_0\cos\theta\)，因此入射通量为 \(P_0=4E_0\cos\theta\)。下文均以各自的入射通量归一化。

解析解

均匀介质与球形 Turbid 冠层

对于厚度为 H 的均匀纯吸收层，光程为 \(H/\cos\theta\)。根据 Beer–Lambert 定律，透射率为

\[ T=\exp\left(-\frac{kH}{\cos\theta}\right). \]

统计黑叶片的透射率为

\[ T=\exp\left(-\frac{G(\theta)\,\mathrm{LAI}}{\cos\theta}\right). \]

球形叶倾角分布的投影系数 \(G=0.5\)。取 \(kH=0.5\,\mathrm{LAI}=\tau\)，两种表示就具有相同的透射率。表中的四组参数对应垂直光学厚度 \(\tau=0,0.1,1,5\)，从真空逐步过渡到强吸收。

吸收层的吸收率为 \(A=1-T\)，黑色地表吸收率为 T。将吸收层分成四个等厚子层，还可逐层比较吸收：深度从顶部向下计，某层的顶部、底部光学厚度分别为 \(\tau_1,\tau_2\)，则

\[ A_{\mathrm{layer}}=e^{-\tau_1/\cos\theta}-e^{-\tau_2/\cos\theta}. \]

离散水平叶片

对随机放置的水平叶片，每条光线在一片叶片所在的高度上，与该叶片相交的概率为 \(a/S\)，其中 a 为叶片面积，S 为单元水平面积。叶片位置相互独立，穿过全部 N 片叶片的统计平均透射率为

\[ \mathbb E[T]=\left(1-\frac aS\right)^N. \]

这是由本例的随机放置方式得到的覆盖概率公式。本例 \(a=0.04\,\mathrm{m^2}\)、\(S=4\,\mathrm{m^2}\)、\(N=100\)，因此

\[ \mathbb E[T]=0.99^{100}=0.36603234. \]

该平均值与太阳天顶角无关。一个具体场景的叶片重叠程度会随观测光路变化，因此分别计算 12 个随机场景，再比较平均结果。

当单片叶面积减小、叶片数量增加，并保持 \(\mathrm{LAI}=Na/S\) 不变时，上式趋近于 \(e^{-\mathrm{LAI}}\)。它也是严格水平统计叶片对应的理论极限。本例 LAI = 1，其极限值为 0.36787944；有限叶片的对比参考采用上面的 \(0.99^{100}\)。Turbid 模拟使用球形分布，水平叶片通过显式三角形表示。

模拟结果

曲线为解析透射率，圆圈和加号分别表示均匀介质和 Turbid 冠层。太阳天顶角越大，路径越长，透射率越低。各后端的结果接近，图中的点会相互重合。

左图对比四个子层的吸收分布，右图展示总吸收与到达黑地表的通量。

左图逐步增加光子数，检查介质和 Turbid 的误差。右图比较离散叶片的平均透射率，误差棒为场景间均值标准误差的 4 倍；实线是有限叶片的解析均值，虚线是无穷小水平叶片的理论极限。

误差评价

指标	阈值
透射率、总吸收率与分层吸收率的最大绝对误差	≤ 0.005
连续介质与球形 Turbid 的透射率、总吸收率差异	≤ 0.005
最后两个光子预算间的最大绝对变化	≤ 0.003
通量闭合误差	≤ 0.0001
逃逸通量比例、未完成路径通量比例	各 ≤ 0.00001
离散叶片平均透射率误差	≤ 0.005，且位于场景间均值的 4 个标准误差内

透射率由黑地表的吸收通量除以入射通量得到，吸收层的吸收率由网格记录的吸收通量计算。离散场景的标准误差为 \(s/\sqrt{12}\)，其中 s 是 12 个场景透射率的样本标准差。这样可以分别查看增加光子数的效果和改变叶片位置产生的波动。

后端	介质最大误差	Turbid 最大误差	两种表示最大差异	离散叶片均值最大误差	结论
optix	0.00015934	0.00015354	0.00016901	0.00260122	通过
vulkan	0.00015508	0.00015230	0.00017426	0.00260119	通过
embree	0.00015911	0.00015330	0.00016901	0.00260109	通过

离散水平叶片的平均透射率如下，解析均值为 0.36603234。

太阳天顶角	OptiX	Vulkan	Embree	场景间均值标准误差（最大）
0°	0.36557375	0.36557375	0.36557346	0.00336840
45°	0.36812525	0.36812525	0.36812511	0.00591345
75°	0.36863356	0.36863353	0.36863343	0.00653107

以上误差均为绝对误差，最终结果使用每次 400 万光子。三个后端均满足评价阈值。

有散射介质的角度响应

沿用上面的 1 m 厚介质层和黑地面，将消光分成吸收与散射两部分。 用消光光学厚度 \(\tau=(k_a+k_s)H\) 和单次散射反照率 \(\omega=k_s/(k_a+k_s)\) 描述介质。 这里比较恰好散射一次后离开介质的辐亮度；模拟结果取单次散射分量。

条件	设置
光照、边界和波长	与上面的介质层相同；太阳天顶角 0°、60°
主要角度组	\(\tau=1\)、\(\omega=0.7\)，相函数参数 \(g=-0.5,0,0.5\)
观测方向	太阳主平面内 −60°、−30°、0°、30°、60°；另加天顶角 60°、相对方位角 90°
光学厚度对照	\(\tau=0.1,3\)，\(\omega=0.7\)、\(g=0\)
散射比例对照	\(\omega=0.2,1\)，\(\tau=1\)、\(g=0\)
对照组方向	太阳天顶角 60°，观测 −60°、0°、60°
采样	8 × 8 像素，覆盖宽度 0.1 m；每像素 128、512、2048 个样本，最终采用 2048

带符号的观测角用于标记太阳主平面的两侧：正值与太阳方位角相同，负值相差 180°。 所有观测均在介质层上方。每个后端共比较 48 组输入。

单次散射解析解

采用 Henyey–Greenstein 相函数：

\[ p_g(\cos\Theta)=\frac{1-g^2}{4\pi(1+g^2-2g\cos\Theta)^{3/2}}. \]

\(\Theta\) 是入射光子的传播方向与出射方向之间的夹角。 \(g=0\) 对应各向同性散射，正值偏向前向，负值偏向后向。 设太阳和观测天顶角为 \(\theta_0,\theta_v\)，相对方位角为 \(\Delta\phi\)，则

\[ \cos\Theta=-\cos\theta_0\cos\theta_v -\sin\theta_0\sin\theta_v\cos\Delta\phi. \]

记 \(\mu_0=\cos\theta_0\)、\(\mu=\cos\theta_v\)，从层顶向下的消光光学深度为 \(t\)。 太阳光到达散射点、再到达层顶的衰减分别为 \(e^{-t/\mu_0}\) 和 \(e^{-t/\mu}\)。 因此单次散射辐亮度的严格解析解为

\[ L_1=\frac{E_0\omega p_g}{\mu}\int_0^\tau e^{-t(1/\mu_0+1/\mu)}\,dt =\frac{E_0\omega p_g\mu_0}{\mu_0+\mu} \left[1-e^{-\tau(1/\mu_0+1/\mu)}\right]. \]

以水平面入射辐照度归一化，得到单次散射 BRF：

\[ \mathrm{BRF}_1=\frac{\pi L_1}{E_0\mu_0} =\frac{\pi\omega p_g}{\mu_0+\mu} \left[1-e^{-\tau(1/\mu_0+1/\mu)}\right]. \]

这条解析式同时约束相函数方向、入射和出射路径的衰减，以及散射比例的线性变化。

上排是两个太阳天顶角下的角度响应，黑线为解析值，彩色标记为模拟结果。 下排检查增加采样后的最大 BRF 误差和分组误差。 每幅影像分为 8 个互不重叠的像素组，用组均值计算标准误差； 以 t 分布和 Bonferroni 校正构造覆盖全部 144 项比较的同时置信区间，族错误率为 1%。

指标	阈值
单次散射 BRF 最大绝对误差	≤ 0.005
单次散射辐亮度最大相对误差	≤ 2%
最后两档采样间的 BRF 最大变化	≤ 0.003
解析值的统计一致性	位于上述同时置信区间内
采样收敛	最终分组相对均方根误差、平均组间方差均小于最低采样档

后端	BRF 最大绝对误差	辐亮度最大相对误差	最后两档 BRF 最大变化	结论
optix	0.00015899	0.1893%	0.00240502	通过
vulkan	0.00015779	0.1893%	0.00241821	通过
embree	0.00015897	0.1893%	0.00240501	通过

144 项比较全部满足上述条件。最终每组输入使用 131072 个样本；分组相对均方根误差从最低预算的约 2.28% 降至 0.54%。

下载与复现

• 场景参数与采样计划
• 吸收层边界（OBJ）与水平叶片（OBJ）
• 复现脚本
• 结果数据

将脚本保存在 scripts/，参数和两个 OBJ 文件放在同级的 data/attenuation/。安装 LESS、NumPy、SciPy 和 matplotlib 后运行：

python scripts/appendix_attenuation.py --output build/attenuation

运行需要 OptiX、Vulkan 和 Embree 三个后端均可用。输出包括一份结果数据和五张 SVG 图。若已下载结果数据，只重算散射部分可加 --scattering-only。