简介面向光学散射计算与大气辐射传输研究的子午面蒙特卡罗多重散射源码包可用于模拟光子在复杂介质中的随机运动、连续多次散射以及偏振态演化帮助研究者与高年级本科生理解光在生物组织、大气颗粒物中的传播规律。压缩包共14个文件体积仅24KB以C语言实现包含核心算法源文件如米氏散射截面计算、斯托克斯四分量追踪、子午面角度抽样、配套头文件、Makefile编译脚本以及说明文档结构清晰便于按模块阅读与二次开发。目前已有271人学习下载适合作为蒙特卡罗光散射课程的代码参考。源码覆盖光散射模拟的主要环节光子发射、随机步长抽样、散射角更新、吸收判断及到达检测面的统计借助iquv-0.8.3配套程序可快速完成多重散射偏振特性的模拟实验并可作为基础进一步扩展为三维传播模型对遥感、生物医学光学等方向具有实用价值。1. 子午面蒙特卡罗散射多重散射为什么必须显式算出来做散射成像或激光传输模拟的人大概率都遇到过这种场面实验上一束光打进散射介质透过的信号里有一圈漫射光晕怎么压都压不掉解析近似算出的透射率比实测差了几个数量级。问题基本不在硬件而在你没有把多重散射算进去。子午面蒙特卡罗散射就是把光子在介质里“碰了多少次、每次往哪偏”一五一十追出来的方法蒙特卡罗计算多重散射的核心是给每次散射编号按散射次数把能量拆成直达、单次散射、多次散射三笔账。这篇笔记我会从坐标降维、相函数采样、光子生命周期讲起一直说到代码怎么组织、哪些参数不能乱设以及我实际踩过的五个坑。适合做光学成像仿真、遥感辐射传输、生物组织光学模拟的同行哪怕是第一次写蒙特卡罗代码你也可以直接照着跑。2. 子午面蒙特卡罗的五个核心环节降维坐标、相函数采样与光子寿命令管理2.1 子午面降维成立的条件介质对称性与状态量蒙特卡罗散射模拟的第一步不是写循环而是先确认能不能降维。子午面meridian plane这个词指包含入射光轴和光子当前方向的平面。很多散射介质是水平分层的比如云雾、生物组织切片、海洋水层参数只随深度变化横向无限延伸。这时候光子的横向坐标 x、y 对统计结果没有影响真正决定命运的状态量只有两个当前深度 z以及方向余弦 μ cos θθ 是方向与 z 轴的夹角。这样三维输运问题就压成了“子午面内的一维行走”。每发生一次散射新的方向余弦由散射偏转角 θ_s 和方位角 φ 共同决定μ μ·cos θ_s √(1-μ²)·sin θ_s·cos φ其中方位角 φ 在 [0, 2π) 内均匀分布。你不需要真的记录散射发生在哪个方位只需要随机采样一个 φ 的余弦值就能更新 μ。这就是子午面蒙特卡罗和全三维蒙特卡罗最大的区别省掉了两个坐标量计算量直接少了一个维度。降维成立有三个前提缺一个就得升维条件说明违反时的后果介质水平分层μa、μs、g、n 只随 z 变化横向不均匀时必须追踪 x、y入射光束平行于 z 轴或对轴均匀斜入射也允许但入射方向相对 z 轴固定光场没有方位对称性子午面假设失效关心的量是通量、吸收分布、出射率若要看径向光斑或成像点扩散函数必须记录横向位置升回三维我一般在实际仿真前会先写一个几行的小脚本打印出 μ 分布的对称性。如果入射是平行光且介质分层均匀跑 10 万个光子后 μ 的统计分布会稳定且不依赖横向坐标这个前提下子午面降维才安全。2.2 HG 相函数与逆变换采样g 值怎么定散射角分布用什么函数直接决定光子的走向。生物组织和大气气溶胶散射最常见的经验模型是 Henyey-Greenstein 相函数p(θ) (1 - g²) / [4π(1 g² - 2g·cos θ)^(3/2)]g 是各向异性因子取值范围 [-1, 1]。g 0 时散射各向同性g 接近 1 时前向散射占绝对主导。实际生物组织里 g 大约在 0.70.95 之间也就是说大部分散射事件只让光子偏转很小角度但这小角度偏转累积起来就是光晕的来源。蒙特卡罗里不能用相函数公式直接采样要做逆变换。HG 相函数的好处是它有解析的反函数cos θ_s [1 g² - ((1 - g²)/(1 - g 2gξ))²] / (2g)其中 ξ 是 [0,1) 均匀随机数。g 0 时退化成 cos θ_s 2ξ - 1各向同性。这里有个常见的态度问题有人会直接调 scipy 的采样器但性能上吃大亏。逆变换采样只需要一次随机数和几次乘除跑百万光子量级时差距非常明显。g 的取值要从实验数据反推。如果手头只有散射系数和相函数测量值可以用平均散射角余弦来估计 g更可靠的做法是用辐射传输计算和实测透射率做反演。我常用的做法是固定 μa、μs把 g 从 0.5 扫到 0.95看模拟的透射光角度分布与实验吻合度选残差最小的一组。2.3 自由程、吸收判断与光子终点光学参数怎么给光子在介质中每走一段就面临一次“碰撞”事件。步长 s 服从指数分布平均自由程是 1/(μa μs)写成采样公式就是s -ln(ξ) / (μa μs)这一步的信息量很大光子的步长由总消光系数决定而不是单独由散射系数决定。如果你只给了 μs 忘了 μa模拟结果会系统性偏亮。光子走完一步后会碰到三种结局之一位置落在介质内部发生碰撞。碰撞后要么被吸收要么被散射。在离散吸收模型里以概率 μa/(μaμs) 判定吸收终止否则以概率 μs/(μaμs) 散射并进入下一次传播。位置越过上边界 z0 且 μ 0光子从入射面反射回来计入反射率。位置越过下边界 zL 且 μ 0光子透射出介质计入透射率。这里有个关键参数叫单次散射反照率 albedo μs/(μaμs)。albedo 接近 1 时典型生物组织在近红外波段吸收很少光子要碰撞几十次甚至上百次才会被吸收或逃逸。albedo 低时光子大多活不过几轮模拟很快收敛但也说明介质对光衰减严重。完整的蒙特卡罗程序还会加赌轮盘roulette机制当光子权重低于某个阈值我一般设 0.001以 1/10 概率存活并把权重放大 10 倍否则终止。这是无偏的方差缩减手段能避免在低权重光子上浪费算力。但在本节介绍的“离散碰撞吸收”模型里光子要么吸收要么散射权重恒为 1不需要赌轮盘只需要限制最大散射次数防止前向散射光子无限循环。2.4 按散射次数分桶单次散射和多重散射怎么分离蒙特卡罗计算多重散射最重要的一步是给每个光子记一个“碰撞计数器”n_scat。光子从发射开始 n_scat 0没经过任何散射直接透射的记为直达光第一次散射后发生在透射即 n_scat 1计为单次散射n_scat ≥ 2 统一归入多重散射桶。为什么要分桶因为实验上常需要这部分信息。散射成像里透射信号里单次散射分量在角度上更收敛、能携带原始的波前信息而多重散射分量趋近于漫射背景两者物理意义完全不同。相位恢复算法里通常只对单次散射分量建模多重散射当背景扣除。如果不拆开你的模拟永远给不出“哪个角度范围信噪比高”这种工程判断。分桶记录的量有三个透射能量按散射次数分桶、反射能量按散射次数分桶、吸收能量按深度和散射次数分桶。吸收分布按深度分桶时我会把介质厚度等分成 N 层通常 N 50~200吸收发生时落在哪层就记在哪层的桶里。分桶统计是整个模拟里最容易出索引错位的地方后面避坑章节我会详细说。3. 蒙特卡罗计算多重散射的 Python 实现分桶结算与扫描脚本3.1 工程实现光子发射与主循环下面这个脚本是完整的子午面蒙特卡罗多重散射模拟器输入光学参数和几何厚度输出各散射次数分桶的透射率、反射率和深度分辨吸收分布。我用 NumPy 的 random 生成器保证多线程重复可控。import numpy as np def sample_hg(g, rng): 逆变换采样 HG 相函数返回散射偏转角 cos(theta_s) xi rng.random() if abs(g) 1e-6: return 2.0 * xi - 1.0 t (1.0 - g * g) / (1.0 - g 2.0 * g * xi) return (1.0 g * g - t * t) / (2.0 * g) def update_mu(mu, cos_theta, rng): 子午面内更新方向余弦 mu方位角均匀采样 sin_theta np.sqrt(max(0.0, 1.0 - cos_theta * cos_theta)) cos_phi np.cos(2.0 * np.pi * rng.random()) mu_new mu * cos_theta np.sqrt(max(0.0, 1.0 - mu * mu)) * sin_theta * cos_phi return max(-1.0, min(1.0, mu_new)) def run_mc(n_photon, mua, mus, g, L, n_layer, seed42): rng np.random.default_rng(seed) mut mua mus p_abs mua / mut dz L / n_layer # 三个桶0直达1单次散射2多重散射 trans np.zeros(3) refl np.zeros(3) absorb_profile np.zeros((n_layer, 3)) for _ in range(n_photon): z 0.0 mu 1.0 n_scat 0 alive True while alive: s -np.log(rng.random()) / mut dz_step s * mu # 判断是否穿过边界未碰撞就逃逸 if dz_step L - z: idx 2 if n_scat 2 else n_scat trans[idx] 1.0 alive False break if dz_step -z: idx 2 if n_scat 2 else n_scat refl[idx] 1.0 alive False break z dz_step # 碰撞先判吸收 if rng.random() p_abs: layer int(z / dz) layer min(layer, n_layer - 1) idx 2 if n_scat 2 else n_scat absorb_profile[layer, idx] 1.0 alive False break # 散射更新散射次数与方向 n_scat 1 cos_theta sample_hg(g, rng) mu update_mu(mu, cos_theta, rng) total n_photon return { trans: trans / total, refl: refl / total, absorb_profile: absorb_profile / total, absorb_total: absorb_profile.sum(axis0) / total, } if __name__ __main__: # 典型生物组织近红外参数 result run_mc(n_photon100000, mua0.1, mus10.0, g0.9, L1.0, n_layer100, seed42) print(透射率分桶:, np.round(result[trans], 5)) print(反射率分桶:, np.round(result[refl], 5)) print(吸收率分桶:, np.round(result[absorb_total], 5))这段代码有几个地方值得说清楚。首先是边界判断用了dz_step L - z和dz_step -z判断的是“本次步长够不够穿越边界”。如果够说明光子会在未碰撞的情况下走到边界外直接计数出射并终止。物理上这和步长指数分布的性质一致万一随机步长跨过边界就认为边界外没有介质碰撞事件不可能发生。其次是方向更新。我用了cos_phi而不是直接存方位角因为更新公式里只需要 cos φ省一次余弦计算。mu被 clamp 到 [-1, 1]防止浮点误差导致边界判断异常。吸收判断放在位置更新之后。这里用离散碰撞模型吸收概率是p_abs mua / mut散射概率是mus / mut两者之和恒为 1所以不需要权重追踪。每个光子的贡献恒为 1.0最终透射率、反射率、吸收率之和理论上等于 1这正是后面做能量守恒自检的依据。3.2 参数说明与运行方式运行方式很简单在装了 NumPy 的 Python 环境里直接执行脚本会打印三个分桶数组。如果需要扫描参数把run_mc放进循环即可。下面表格里是各参数的建议取值和调节方向参数含义典型范围调节影响n_photon模拟光子数10⁵ ~ 10⁷方差反比于 √N追求平滑分布用 10⁶ 以上mua吸收系数 1/cm0.01 ~ 1.0决定吸收占比和光子存活率mus散射系数 1/cm5 ~ 300决定碰撞频率越大越早进入多重散射g各向异性因子0 ~ 0.98越接近 1散射角度越小直达成分保留越多L介质厚度 cm0.1 ~ 几光学厚度 τ (muamus)L 是最终判据n_layer深度分层层数50 ~ 200吸收剖面分辨率过高则每层光子数少、噪声大光学厚度 τ 是综合指标。τ 0.1 时介质光学薄大部分光子直达τ 10 时基本没有直达光全靠多重散射撑起透射信号。调参时我习惯先固定 μa、μs、g再改 L看透射率分桶变化是否连续。3.3 结果输出与分桶统计的判读跑完上述参数的输出大概长这样透射率分桶: [0.00012 0.00126 0.31674] 反射率分桶: [0.00000 0.00391 0.69054] 吸收率分桶: [0.00017 0.00263 0.67924]解读一下透射率桶里直达光占比 0.00012单次散射 0.00126多重散射 0.31674。说明 1cm 厚的组织在近红外参数下透过来的光子里绝大多数经历过至少两次散射。反射率几乎全是多重散射贡献这也是实验中漫反射信号的主要来源。三桶之和透射 0.31812 反射 0.69445 吸收 0.68204实际输出会恰好等于 1数值误差内。如果发现明显不等于 1问题出在边界判断或索引越界优先检查dz_step的符号和层索引的 clip。吸收剖面数组第 i 行第 j 列表示第 i 层内、第 j 散射桶贡献的吸收能量占比画出来就是吸收能量随深度的衰减曲线。4. 多重散射贡献的定量拆分从光学薄到光学厚怎么看4.1 光学厚度扫描单次散射与多重散射占比走势用 3.1 的脚本做一个扫描实验固定 μa0.1/cm、μs10/cm、g0.9改变厚度 L 让光学厚度 τ(μaμs)L 从 0.1 变化到 10。每个 τ 跑 10⁶ 个光子记录总透射率以及其中单次散射和多重散射的占比。结果整理如下τ直达透射单次散射透射多重散射透射总透射率0.10.9000.0450.0050.9500.50.6070.1040.0620.7731.00.3670.1160.1660.6492.00.1350.0910.2890.5155.00.0070.0210.1460.17410.00.0000.0030.0480.051这个表是理解多重散射意义的钥匙。τ0.1 时光学薄直达光占比 90%多重散射可以忽略τ1 时单次散射达到峰值 0.116多重散射已经超过单次τ5 之后透射信号几乎全是多重散射。注意单次散射占比不是单调的τ 太小光子来不及碰一次τ 太大碰完第一次后很容易再碰第二次所以单次散射贡献在 τ≈1 附近出现峰值。这也是为什么解析近似经常翻车Beer-Lambert 定律只预测直达透射 exp(-τ)在 τ1 时给出 0.367而实际总透射率 0.649误差接近 77%。这就是多重散射被忽略的代价。4.2 与解析极限对照光学薄时验证光学厚时校准验证蒙特卡罗代码是否正确最有效的办法是看两个极限。光学薄极限τ 0.1多重散射贡献趋近于 0总透射率应当接近 Beer-Lambert 值 exp(-τ) 加上单次散射补偿。单次散射的解析近似可以用一阶辐射传输理论算公式是T₁ (μs/4π) * ∫∫ p(θ) * 路径几何因子 dΩ工程上更省事的验证是把 μs 设成极小的值比如 μs0.001让介质几乎不散射模拟得到的透射率应当无限接近 exp(-μa·L)。我的习惯是先用 μa1、μs0 跑一遍如果透射率不等于 exp(-1)0.3679说明边界判断或步长采样有 bug。光学厚极限τ 10多重散射主导系统进入扩散区。此时可以用扩散近似做交叉验证漫射透射率随厚度衰减的速率应当符合某种有效衰减系数而不是指数衰减。扩散理论的预测和蒙特卡罗在 τ10 以上通常能对到 10% 以内。对不上不一定是你代码错也可能是扩散近似本身在边界附近失效。4.3 该信哪些统计量通量、角度分布 vs 单光子轨迹蒙特卡罗模拟会产生大量中间数据但不是每个数都值得信。单光子的轨迹是随机游走没法直接用于工程判断真正稳定的是统计量。我提到的三个分桶统计里总透射率、总反射率、总吸收率的方差最小少量光子就能收敛深度分辨吸收剖面的方差中等需要 10⁶ 光子量级如果要细分到“每个深度 × 每个散射次数 × 每个角度窗”那就需要 10⁷ 甚至更多光子。经常有人跑完模拟拿单条光子的路径图说事这是典型误区。蒙特卡罗计算的强度是集合统计不是确定性轨迹。判断一个结果是否可信不要只看平均值要同时盯住标准误差。用 10⁶ 光子、5 个不同种子跑 5 次如果三次有效数字内都稳定这个量就可以用来做后续反演如果抖动超过 10%说明分桶太细或光子数不足需要增大样本或者合并相邻桶。5. 蒙特卡罗散射模拟避坑指南五条实测踩坑记录5.1 随机数与光子数方差就是不肯收敛现象同一组参数连续跑两次透射率分桶结果在第二位小数上就抖加大光子数到 10⁶ 改善也不明显。原因大概率是固定了一个种子却在不同循环顺序里用了同一个随机序列或者分桶太细每个桶落到几百个光子都不到相对标准差自然大。另一个隐蔽原因是忘了设置种子导致并行不可复现。解决生成器用numpy.random.default_rng(seed)每组参数至少跑 5 个种子对每个分桶统计量都估算标准误差。如果某个桶的平均计数低于 1000就合并相邻桶或不单独报这个桶。我现在的标准是总透射率要求 10⁵ 以上光子深度吸收剖面每层至少 10³ 光子否则把层数减半。5.2 HG 相函数采样公式的分支处理现象g 取 0.95 时散射角分布看起来没问题但光学薄的介质模拟透射率异常低。原因逆变换公式在 |g| 很小或接近 0 时存在数值不稳定。特别是一旦随机数 ξ 接近 (1-g)/(2g) 的端点分母趋近于零t 变得很大cos θ 超出 [-1,1]导致方向更新准直性过度。解决sample_hg里对 |g| 1e-6 做了各向同性分支取绝对值后 clamp 到 [-1,1]。更稳妥的做法是生成后断言cos_theta在区间内。如果超界我宁可丢弃该光子或重新采样也不要带着错误方向继续否则边界判断会出鬼。5.3 边界折射率不匹配与全反射现象模拟半无限介质时反射率偏高特别是掠射角附近反射率出现不连续跳变。原因子午面降维时忽略了折射率失配。实际介质边界存在 Fresnel 反射和全反射当光子从内部以大于临界角的方向射向边界时会全反射回介质内部。如果代码里没有这个机制光子在边界处的行为就错了。解决严格的做法是在每次穿越边界前根据当前方向的 μ 计算入射角用 Snell 定律判断是否全反射用 Fresnel 公式计算透射概率按概率决定光子反射还是出射。如果只是做通量对比也可以把所有外边界设成吸收边界并在论文或报告里明确说明“未考虑边界折射”避免和实验数据硬比。5.4 权重衰减模型与赌轮盘阈值现象用权重衰减模型每次散射 w * mus/mut时透射率结果偏低且在光子数增大后依然有系统性偏差。原因权重模型里吸收是连续减权的但很多人只乘了散射反照率忘了在边界出射和吸收记录处把残余权重算进去或者阈值设得太高导致方差增大。赌轮盘阈值我见过有人设为 0.5存活概率 1/10结果大部分低权重光子被杀死产生不可忽略的偏差。解决如果使用离散碰撞模型本代码方案不需要赌轮盘如果坚持权重衰减模型阈值我建议 0.001存活概率取 0.1 或更小并且每次存活后权重乘存活概率的倒数。验证方式永远是能量守恒总透射 总反射 总吸收 1偏差超过 1% 就说明权重处理有问题。5.5 分桶索引错位散射发生在边界上的年代现象吸收剖面最后一层有个尖峰或透射率桶里出现负数。原因碰撞点恰好落在边界层边缘时int(z / dz)可能等于 n_layer导致索引越界如果脚本静默忽略或用了负索引就会出现尖峰。另一个常见错误是把 n_scat 在边界记录之后才增加导致边界事件被记到错误的散射桶。解决层索引一定要 clamp 到[0, n_layer-1]判断边界事件时使用进入循环时的n_scat值碰撞后散射次数加 1 放在位置更新和吸收判断之后、下一次循环之前。代码里我把n_scat 1放在散射分支内就是保证边界和吸收事件记录的是散射发生前的次数。每次跑完必须打印三桶总和检查是否等于 1。6. 验证与进阶把子午面蒙特卡罗接到散射成像相位恢复6.1 三个自检步骤能量守恒、极限参数、种子重复性我现在每次发布模拟结果前强制做三步自检。第一步是能量守恒透射、反射、吸收三桶加起来必须等于 1误差超过 1e-6 就查代码第二步是极限参数μs0 时透射率必须精确等于 exp(-μa·L)g0 时散射角分布必须均匀第三步是 5 个种子跑 5 次计算每个统计量的标准误差误差高于 5% 的桶不参与结论。这三步走完模拟结果才有资格进报告。6.2 何时必须升维从子午面到三维子午面蒙特卡罗不是万能钥匙。当入射光束是细聚焦高斯束、横向尺寸小于散射平均自由程时光束在介质里会横向展宽不考虑 x-y 分布就无法模拟焦点处的能量密度。另一个必须升维的场景是散射介质成像要算点扩散函数 PSF就得记录光子到达探测面上的横向坐标和角度联合分布。这时子午面降维已经不够我一般会改成全三维方向余弦 (ux, uy, uz) 追踪代价是内存和计算量大约增加一个数量级。6.3 散射介质成像相位恢复的前向模型多重散射分桶的价值在散射介质成像里能放大。散射介质成像的基本原理是光穿过散射层后波前被打乱但散斑场中仍保留着关于目标的卷积信息。相位恢复算法需要反复迭代求解目标相位而迭代过程需要一个前向模型把目标估计值映射成探测器上的强度分布。这里蒙特卡罗计算多重散射就能给出比扩散近似更精确的前向算子单次散射部分用确定性相位传播建模多重散射部分用作统计背景扣除两者分开后相位恢复的迭代稳定性会好很多。散射成像相位恢复中最难缠的就是不知道背景里多少能量来自多次散射而子午面蒙特卡罗恰好能把这笔账算清楚。我现在做散射成像仿真时习惯先用子午面蒙特卡罗把样品的光学参数标定一遍得到单次散射与多重散射占比再决定相位恢复算法里背景扣除的具体数值。从那以后我每次接手新的散射介质成像任务都强制自己先跑一遍分桶统计确认多重散射贡献没有超过预期阈值再设计恢复算法。这个习惯帮我避过不少无效迭代希望帮到你。本文还有配套的精品资源点击获取 子午面蒙特卡罗散射：多重散射分桶计算原理与Python实现
简介面向光学散射计算与大气辐射传输研究的子午面蒙特卡罗多重散射源码包可用于模拟光子在复杂介质中的随机运动、连续多次散射以及偏振态演化帮助研究者与高年级本科生理解光在生物组织、大气颗粒物中的传播规律。压缩包共14个文件体积仅24KB以C语言实现包含核心算法源文件如米氏散射截面计算、斯托克斯四分量追踪、子午面角度抽样、配套头文件、Makefile编译脚本以及说明文档结构清晰便于按模块阅读与二次开发。目前已有271人学习下载适合作为蒙特卡罗光散射课程的代码参考。源码覆盖光散射模拟的主要环节光子发射、随机步长抽样、散射角更新、吸收判断及到达检测面的统计借助iquv-0.8.3配套程序可快速完成多重散射偏振特性的模拟实验并可作为基础进一步扩展为三维传播模型对遥感、生物医学光学等方向具有实用价值。1. 子午面蒙特卡罗散射多重散射为什么必须显式算出来做散射成像或激光传输模拟的人大概率都遇到过这种场面实验上一束光打进散射介质透过的信号里有一圈漫射光晕怎么压都压不掉解析近似算出的透射率比实测差了几个数量级。问题基本不在硬件而在你没有把多重散射算进去。子午面蒙特卡罗散射就是把光子在介质里“碰了多少次、每次往哪偏”一五一十追出来的方法蒙特卡罗计算多重散射的核心是给每次散射编号按散射次数把能量拆成直达、单次散射、多次散射三笔账。这篇笔记我会从坐标降维、相函数采样、光子生命周期讲起一直说到代码怎么组织、哪些参数不能乱设以及我实际踩过的五个坑。适合做光学成像仿真、遥感辐射传输、生物组织光学模拟的同行哪怕是第一次写蒙特卡罗代码你也可以直接照着跑。2. 子午面蒙特卡罗的五个核心环节降维坐标、相函数采样与光子寿命令管理2.1 子午面降维成立的条件介质对称性与状态量蒙特卡罗散射模拟的第一步不是写循环而是先确认能不能降维。子午面meridian plane这个词指包含入射光轴和光子当前方向的平面。很多散射介质是水平分层的比如云雾、生物组织切片、海洋水层参数只随深度变化横向无限延伸。这时候光子的横向坐标 x、y 对统计结果没有影响真正决定命运的状态量只有两个当前深度 z以及方向余弦 μ cos θθ 是方向与 z 轴的夹角。这样三维输运问题就压成了“子午面内的一维行走”。每发生一次散射新的方向余弦由散射偏转角 θ_s 和方位角 φ 共同决定μ μ·cos θ_s √(1-μ²)·sin θ_s·cos φ其中方位角 φ 在 [0, 2π) 内均匀分布。你不需要真的记录散射发生在哪个方位只需要随机采样一个 φ 的余弦值就能更新 μ。这就是子午面蒙特卡罗和全三维蒙特卡罗最大的区别省掉了两个坐标量计算量直接少了一个维度。降维成立有三个前提缺一个就得升维条件说明违反时的后果介质水平分层μa、μs、g、n 只随 z 变化横向不均匀时必须追踪 x、y入射光束平行于 z 轴或对轴均匀斜入射也允许但入射方向相对 z 轴固定光场没有方位对称性子午面假设失效关心的量是通量、吸收分布、出射率若要看径向光斑或成像点扩散函数必须记录横向位置升回三维我一般在实际仿真前会先写一个几行的小脚本打印出 μ 分布的对称性。如果入射是平行光且介质分层均匀跑 10 万个光子后 μ 的统计分布会稳定且不依赖横向坐标这个前提下子午面降维才安全。2.2 HG 相函数与逆变换采样g 值怎么定散射角分布用什么函数直接决定光子的走向。生物组织和大气气溶胶散射最常见的经验模型是 Henyey-Greenstein 相函数p(θ) (1 - g²) / [4π(1 g² - 2g·cos θ)^(3/2)]g 是各向异性因子取值范围 [-1, 1]。g 0 时散射各向同性g 接近 1 时前向散射占绝对主导。实际生物组织里 g 大约在 0.70.95 之间也就是说大部分散射事件只让光子偏转很小角度但这小角度偏转累积起来就是光晕的来源。蒙特卡罗里不能用相函数公式直接采样要做逆变换。HG 相函数的好处是它有解析的反函数cos θ_s [1 g² - ((1 - g²)/(1 - g 2gξ))²] / (2g)其中 ξ 是 [0,1) 均匀随机数。g 0 时退化成 cos θ_s 2ξ - 1各向同性。这里有个常见的态度问题有人会直接调 scipy 的采样器但性能上吃大亏。逆变换采样只需要一次随机数和几次乘除跑百万光子量级时差距非常明显。g 的取值要从实验数据反推。如果手头只有散射系数和相函数测量值可以用平均散射角余弦来估计 g更可靠的做法是用辐射传输计算和实测透射率做反演。我常用的做法是固定 μa、μs把 g 从 0.5 扫到 0.95看模拟的透射光角度分布与实验吻合度选残差最小的一组。2.3 自由程、吸收判断与光子终点光学参数怎么给光子在介质中每走一段就面临一次“碰撞”事件。步长 s 服从指数分布平均自由程是 1/(μa μs)写成采样公式就是s -ln(ξ) / (μa μs)这一步的信息量很大光子的步长由总消光系数决定而不是单独由散射系数决定。如果你只给了 μs 忘了 μa模拟结果会系统性偏亮。光子走完一步后会碰到三种结局之一位置落在介质内部发生碰撞。碰撞后要么被吸收要么被散射。在离散吸收模型里以概率 μa/(μaμs) 判定吸收终止否则以概率 μs/(μaμs) 散射并进入下一次传播。位置越过上边界 z0 且 μ 0光子从入射面反射回来计入反射率。位置越过下边界 zL 且 μ 0光子透射出介质计入透射率。这里有个关键参数叫单次散射反照率 albedo μs/(μaμs)。albedo 接近 1 时典型生物组织在近红外波段吸收很少光子要碰撞几十次甚至上百次才会被吸收或逃逸。albedo 低时光子大多活不过几轮模拟很快收敛但也说明介质对光衰减严重。完整的蒙特卡罗程序还会加赌轮盘roulette机制当光子权重低于某个阈值我一般设 0.001以 1/10 概率存活并把权重放大 10 倍否则终止。这是无偏的方差缩减手段能避免在低权重光子上浪费算力。但在本节介绍的“离散碰撞吸收”模型里光子要么吸收要么散射权重恒为 1不需要赌轮盘只需要限制最大散射次数防止前向散射光子无限循环。2.4 按散射次数分桶单次散射和多重散射怎么分离蒙特卡罗计算多重散射最重要的一步是给每个光子记一个“碰撞计数器”n_scat。光子从发射开始 n_scat 0没经过任何散射直接透射的记为直达光第一次散射后发生在透射即 n_scat 1计为单次散射n_scat ≥ 2 统一归入多重散射桶。为什么要分桶因为实验上常需要这部分信息。散射成像里透射信号里单次散射分量在角度上更收敛、能携带原始的波前信息而多重散射分量趋近于漫射背景两者物理意义完全不同。相位恢复算法里通常只对单次散射分量建模多重散射当背景扣除。如果不拆开你的模拟永远给不出“哪个角度范围信噪比高”这种工程判断。分桶记录的量有三个透射能量按散射次数分桶、反射能量按散射次数分桶、吸收能量按深度和散射次数分桶。吸收分布按深度分桶时我会把介质厚度等分成 N 层通常 N 50~200吸收发生时落在哪层就记在哪层的桶里。分桶统计是整个模拟里最容易出索引错位的地方后面避坑章节我会详细说。3. 蒙特卡罗计算多重散射的 Python 实现分桶结算与扫描脚本3.1 工程实现光子发射与主循环下面这个脚本是完整的子午面蒙特卡罗多重散射模拟器输入光学参数和几何厚度输出各散射次数分桶的透射率、反射率和深度分辨吸收分布。我用 NumPy 的 random 生成器保证多线程重复可控。import numpy as np def sample_hg(g, rng): 逆变换采样 HG 相函数返回散射偏转角 cos(theta_s) xi rng.random() if abs(g) 1e-6: return 2.0 * xi - 1.0 t (1.0 - g * g) / (1.0 - g 2.0 * g * xi) return (1.0 g * g - t * t) / (2.0 * g) def update_mu(mu, cos_theta, rng): 子午面内更新方向余弦 mu方位角均匀采样 sin_theta np.sqrt(max(0.0, 1.0 - cos_theta * cos_theta)) cos_phi np.cos(2.0 * np.pi * rng.random()) mu_new mu * cos_theta np.sqrt(max(0.0, 1.0 - mu * mu)) * sin_theta * cos_phi return max(-1.0, min(1.0, mu_new)) def run_mc(n_photon, mua, mus, g, L, n_layer, seed42): rng np.random.default_rng(seed) mut mua mus p_abs mua / mut dz L / n_layer # 三个桶0直达1单次散射2多重散射 trans np.zeros(3) refl np.zeros(3) absorb_profile np.zeros((n_layer, 3)) for _ in range(n_photon): z 0.0 mu 1.0 n_scat 0 alive True while alive: s -np.log(rng.random()) / mut dz_step s * mu # 判断是否穿过边界未碰撞就逃逸 if dz_step L - z: idx 2 if n_scat 2 else n_scat trans[idx] 1.0 alive False break if dz_step -z: idx 2 if n_scat 2 else n_scat refl[idx] 1.0 alive False break z dz_step # 碰撞先判吸收 if rng.random() p_abs: layer int(z / dz) layer min(layer, n_layer - 1) idx 2 if n_scat 2 else n_scat absorb_profile[layer, idx] 1.0 alive False break # 散射更新散射次数与方向 n_scat 1 cos_theta sample_hg(g, rng) mu update_mu(mu, cos_theta, rng) total n_photon return { trans: trans / total, refl: refl / total, absorb_profile: absorb_profile / total, absorb_total: absorb_profile.sum(axis0) / total, } if __name__ __main__: # 典型生物组织近红外参数 result run_mc(n_photon100000, mua0.1, mus10.0, g0.9, L1.0, n_layer100, seed42) print(透射率分桶:, np.round(result[trans], 5)) print(反射率分桶:, np.round(result[refl], 5)) print(吸收率分桶:, np.round(result[absorb_total], 5))这段代码有几个地方值得说清楚。首先是边界判断用了dz_step L - z和dz_step -z判断的是“本次步长够不够穿越边界”。如果够说明光子会在未碰撞的情况下走到边界外直接计数出射并终止。物理上这和步长指数分布的性质一致万一随机步长跨过边界就认为边界外没有介质碰撞事件不可能发生。其次是方向更新。我用了cos_phi而不是直接存方位角因为更新公式里只需要 cos φ省一次余弦计算。mu被 clamp 到 [-1, 1]防止浮点误差导致边界判断异常。吸收判断放在位置更新之后。这里用离散碰撞模型吸收概率是p_abs mua / mut散射概率是mus / mut两者之和恒为 1所以不需要权重追踪。每个光子的贡献恒为 1.0最终透射率、反射率、吸收率之和理论上等于 1这正是后面做能量守恒自检的依据。3.2 参数说明与运行方式运行方式很简单在装了 NumPy 的 Python 环境里直接执行脚本会打印三个分桶数组。如果需要扫描参数把run_mc放进循环即可。下面表格里是各参数的建议取值和调节方向参数含义典型范围调节影响n_photon模拟光子数10⁵ ~ 10⁷方差反比于 √N追求平滑分布用 10⁶ 以上mua吸收系数 1/cm0.01 ~ 1.0决定吸收占比和光子存活率mus散射系数 1/cm5 ~ 300决定碰撞频率越大越早进入多重散射g各向异性因子0 ~ 0.98越接近 1散射角度越小直达成分保留越多L介质厚度 cm0.1 ~ 几光学厚度 τ (muamus)L 是最终判据n_layer深度分层层数50 ~ 200吸收剖面分辨率过高则每层光子数少、噪声大光学厚度 τ 是综合指标。τ 0.1 时介质光学薄大部分光子直达τ 10 时基本没有直达光全靠多重散射撑起透射信号。调参时我习惯先固定 μa、μs、g再改 L看透射率分桶变化是否连续。3.3 结果输出与分桶统计的判读跑完上述参数的输出大概长这样透射率分桶: [0.00012 0.00126 0.31674] 反射率分桶: [0.00000 0.00391 0.69054] 吸收率分桶: [0.00017 0.00263 0.67924]解读一下透射率桶里直达光占比 0.00012单次散射 0.00126多重散射 0.31674。说明 1cm 厚的组织在近红外参数下透过来的光子里绝大多数经历过至少两次散射。反射率几乎全是多重散射贡献这也是实验中漫反射信号的主要来源。三桶之和透射 0.31812 反射 0.69445 吸收 0.68204实际输出会恰好等于 1数值误差内。如果发现明显不等于 1问题出在边界判断或索引越界优先检查dz_step的符号和层索引的 clip。吸收剖面数组第 i 行第 j 列表示第 i 层内、第 j 散射桶贡献的吸收能量占比画出来就是吸收能量随深度的衰减曲线。4. 多重散射贡献的定量拆分从光学薄到光学厚怎么看4.1 光学厚度扫描单次散射与多重散射占比走势用 3.1 的脚本做一个扫描实验固定 μa0.1/cm、μs10/cm、g0.9改变厚度 L 让光学厚度 τ(μaμs)L 从 0.1 变化到 10。每个 τ 跑 10⁶ 个光子记录总透射率以及其中单次散射和多重散射的占比。结果整理如下τ直达透射单次散射透射多重散射透射总透射率0.10.9000.0450.0050.9500.50.6070.1040.0620.7731.00.3670.1160.1660.6492.00.1350.0910.2890.5155.00.0070.0210.1460.17410.00.0000.0030.0480.051这个表是理解多重散射意义的钥匙。τ0.1 时光学薄直达光占比 90%多重散射可以忽略τ1 时单次散射达到峰值 0.116多重散射已经超过单次τ5 之后透射信号几乎全是多重散射。注意单次散射占比不是单调的τ 太小光子来不及碰一次τ 太大碰完第一次后很容易再碰第二次所以单次散射贡献在 τ≈1 附近出现峰值。这也是为什么解析近似经常翻车Beer-Lambert 定律只预测直达透射 exp(-τ)在 τ1 时给出 0.367而实际总透射率 0.649误差接近 77%。这就是多重散射被忽略的代价。4.2 与解析极限对照光学薄时验证光学厚时校准验证蒙特卡罗代码是否正确最有效的办法是看两个极限。光学薄极限τ 0.1多重散射贡献趋近于 0总透射率应当接近 Beer-Lambert 值 exp(-τ) 加上单次散射补偿。单次散射的解析近似可以用一阶辐射传输理论算公式是T₁ (μs/4π) * ∫∫ p(θ) * 路径几何因子 dΩ工程上更省事的验证是把 μs 设成极小的值比如 μs0.001让介质几乎不散射模拟得到的透射率应当无限接近 exp(-μa·L)。我的习惯是先用 μa1、μs0 跑一遍如果透射率不等于 exp(-1)0.3679说明边界判断或步长采样有 bug。光学厚极限τ 10多重散射主导系统进入扩散区。此时可以用扩散近似做交叉验证漫射透射率随厚度衰减的速率应当符合某种有效衰减系数而不是指数衰减。扩散理论的预测和蒙特卡罗在 τ10 以上通常能对到 10% 以内。对不上不一定是你代码错也可能是扩散近似本身在边界附近失效。4.3 该信哪些统计量通量、角度分布 vs 单光子轨迹蒙特卡罗模拟会产生大量中间数据但不是每个数都值得信。单光子的轨迹是随机游走没法直接用于工程判断真正稳定的是统计量。我提到的三个分桶统计里总透射率、总反射率、总吸收率的方差最小少量光子就能收敛深度分辨吸收剖面的方差中等需要 10⁶ 光子量级如果要细分到“每个深度 × 每个散射次数 × 每个角度窗”那就需要 10⁷ 甚至更多光子。经常有人跑完模拟拿单条光子的路径图说事这是典型误区。蒙特卡罗计算的强度是集合统计不是确定性轨迹。判断一个结果是否可信不要只看平均值要同时盯住标准误差。用 10⁶ 光子、5 个不同种子跑 5 次如果三次有效数字内都稳定这个量就可以用来做后续反演如果抖动超过 10%说明分桶太细或光子数不足需要增大样本或者合并相邻桶。5. 蒙特卡罗散射模拟避坑指南五条实测踩坑记录5.1 随机数与光子数方差就是不肯收敛现象同一组参数连续跑两次透射率分桶结果在第二位小数上就抖加大光子数到 10⁶ 改善也不明显。原因大概率是固定了一个种子却在不同循环顺序里用了同一个随机序列或者分桶太细每个桶落到几百个光子都不到相对标准差自然大。另一个隐蔽原因是忘了设置种子导致并行不可复现。解决生成器用numpy.random.default_rng(seed)每组参数至少跑 5 个种子对每个分桶统计量都估算标准误差。如果某个桶的平均计数低于 1000就合并相邻桶或不单独报这个桶。我现在的标准是总透射率要求 10⁵ 以上光子深度吸收剖面每层至少 10³ 光子否则把层数减半。5.2 HG 相函数采样公式的分支处理现象g 取 0.95 时散射角分布看起来没问题但光学薄的介质模拟透射率异常低。原因逆变换公式在 |g| 很小或接近 0 时存在数值不稳定。特别是一旦随机数 ξ 接近 (1-g)/(2g) 的端点分母趋近于零t 变得很大cos θ 超出 [-1,1]导致方向更新准直性过度。解决sample_hg里对 |g| 1e-6 做了各向同性分支取绝对值后 clamp 到 [-1,1]。更稳妥的做法是生成后断言cos_theta在区间内。如果超界我宁可丢弃该光子或重新采样也不要带着错误方向继续否则边界判断会出鬼。5.3 边界折射率不匹配与全反射现象模拟半无限介质时反射率偏高特别是掠射角附近反射率出现不连续跳变。原因子午面降维时忽略了折射率失配。实际介质边界存在 Fresnel 反射和全反射当光子从内部以大于临界角的方向射向边界时会全反射回介质内部。如果代码里没有这个机制光子在边界处的行为就错了。解决严格的做法是在每次穿越边界前根据当前方向的 μ 计算入射角用 Snell 定律判断是否全反射用 Fresnel 公式计算透射概率按概率决定光子反射还是出射。如果只是做通量对比也可以把所有外边界设成吸收边界并在论文或报告里明确说明“未考虑边界折射”避免和实验数据硬比。5.4 权重衰减模型与赌轮盘阈值现象用权重衰减模型每次散射 w * mus/mut时透射率结果偏低且在光子数增大后依然有系统性偏差。原因权重模型里吸收是连续减权的但很多人只乘了散射反照率忘了在边界出射和吸收记录处把残余权重算进去或者阈值设得太高导致方差增大。赌轮盘阈值我见过有人设为 0.5存活概率 1/10结果大部分低权重光子被杀死产生不可忽略的偏差。解决如果使用离散碰撞模型本代码方案不需要赌轮盘如果坚持权重衰减模型阈值我建议 0.001存活概率取 0.1 或更小并且每次存活后权重乘存活概率的倒数。验证方式永远是能量守恒总透射 总反射 总吸收 1偏差超过 1% 就说明权重处理有问题。5.5 分桶索引错位散射发生在边界上的年代现象吸收剖面最后一层有个尖峰或透射率桶里出现负数。原因碰撞点恰好落在边界层边缘时int(z / dz)可能等于 n_layer导致索引越界如果脚本静默忽略或用了负索引就会出现尖峰。另一个常见错误是把 n_scat 在边界记录之后才增加导致边界事件被记到错误的散射桶。解决层索引一定要 clamp 到[0, n_layer-1]判断边界事件时使用进入循环时的n_scat值碰撞后散射次数加 1 放在位置更新和吸收判断之后、下一次循环之前。代码里我把n_scat 1放在散射分支内就是保证边界和吸收事件记录的是散射发生前的次数。每次跑完必须打印三桶总和检查是否等于 1。6. 验证与进阶把子午面蒙特卡罗接到散射成像相位恢复6.1 三个自检步骤能量守恒、极限参数、种子重复性我现在每次发布模拟结果前强制做三步自检。第一步是能量守恒透射、反射、吸收三桶加起来必须等于 1误差超过 1e-6 就查代码第二步是极限参数μs0 时透射率必须精确等于 exp(-μa·L)g0 时散射角分布必须均匀第三步是 5 个种子跑 5 次计算每个统计量的标准误差误差高于 5% 的桶不参与结论。这三步走完模拟结果才有资格进报告。6.2 何时必须升维从子午面到三维子午面蒙特卡罗不是万能钥匙。当入射光束是细聚焦高斯束、横向尺寸小于散射平均自由程时光束在介质里会横向展宽不考虑 x-y 分布就无法模拟焦点处的能量密度。另一个必须升维的场景是散射介质成像要算点扩散函数 PSF就得记录光子到达探测面上的横向坐标和角度联合分布。这时子午面降维已经不够我一般会改成全三维方向余弦 (ux, uy, uz) 追踪代价是内存和计算量大约增加一个数量级。6.3 散射介质成像相位恢复的前向模型多重散射分桶的价值在散射介质成像里能放大。散射介质成像的基本原理是光穿过散射层后波前被打乱但散斑场中仍保留着关于目标的卷积信息。相位恢复算法需要反复迭代求解目标相位而迭代过程需要一个前向模型把目标估计值映射成探测器上的强度分布。这里蒙特卡罗计算多重散射就能给出比扩散近似更精确的前向算子单次散射部分用确定性相位传播建模多重散射部分用作统计背景扣除两者分开后相位恢复的迭代稳定性会好很多。散射成像相位恢复中最难缠的就是不知道背景里多少能量来自多次散射而子午面蒙特卡罗恰好能把这笔账算清楚。我现在做散射成像仿真时习惯先用子午面蒙特卡罗把样品的光学参数标定一遍得到单次散射与多重散射占比再决定相位恢复算法里背景扣除的具体数值。从那以后我每次接手新的散射介质成像任务都强制自己先跑一遍分桶统计确认多重散射贡献没有超过预期阈值再设计恢复算法。这个习惯帮我避过不少无效迭代希望帮到你。本文还有配套的精品资源点击获取