蒙特卡洛光线追踪的(准)蒙特卡洛积分基础知识系列,上一部分讲完了随机变量、概率密度函数、期望、方差这些基本概念,也把蒙特卡洛积分估计子的来龙去脉梳理了一遍。但估计子写出来只是第一步,真正让一个渲染器从"能出图"到"出好图"的,是对方差的理解和控制。这一部分就是接着这个话题往下走。
这篇文章不是零基础科普,而是"基础知识二"——适合已经掌握蒙特卡洛积分基本框架的读者,最好你正在写自己的路径追踪器,或者至少读过几套开源渲染器的核心采样代码。我们会重点拆解三块内容:重要性采样为什么是控制方差的第一武器、(准)蒙特卡洛方法里的低差异序列是怎么工作的、多重重要性采样又是怎么把不同策略组合出更高效率的。每块内容我都会配合实际渲染场景来讲,最后把我调试渲染器时踩过的坑一并倒出来。
1. 方差为什么是渲染画质的隐形天花板
1.1 估计子的收敛速度:为什么砸采样数是最后的笨办法
回到上一部分的基础公式。要估计积分 I = ∫₀¹ f(x)dx,蒙特卡洛估计子的形式是:
F = (1/N) · Σᵢ f(Xᵢ) / p(Xᵢ)
其中 Xᵢ 独立采样自概率密度函数 p。无偏性来自期望的线性:每个 f(Xᵢ)/p(Xᵢ) 的期望恰好等于 I。但无偏不等于好用,关键看方差。在 p 的支撑集覆盖 f 的支撑集(也就是不存在某些区域 p 为零而 f 非零)的前提下,方差可以写成:
Var[F] = (1/N) · Var[f(X)/p(X)]
这个式子看起来平淡,其实是全文最重要的一个公式。它告诉我们两件事:第一,估计子的方差随 N 线性下降;第二,方差的大小直接由比值 f/p 的波动程度决定,与单独的 f 无关。
再往下推一步,中心极限定理告诉我们,当 N 足够大时,估计子 F 的误差大致服从以 0 为均值、标准差正比于 1/√N 的正态分布。换句话说,均方根误差(RMSE)的收敛率是 O(1/√N)。这意味着什么?想让误差降一半,采样数要涨到四倍。想从 1024 个样本压到和 4096 个样本一样的噪点水平,你只能把每像素采样数乘 4,渲染时间也差不多乘 4。渲染器跑一帧要几分钟的时候,这条路基本走不通。
顺带回答一个很多从数值分析转过来的人会问的问题:为什么不用梯形法、辛普森法那种收敛率漂亮得多的确定性求积?答案很简单——维度。梯形法在一维是 O(1/N²),辛普森是 O(1/N⁴),看起来比蒙特卡洛好太多,但要在 d 维空间维持这个收敛率,需要的网格点是 N^d 量级。路径追踪里一条 8 次弹射的路径就要消耗 16 个随机数来描述各段方向,也就是 16 维的积分域。网格法在这个维度下早就变成天文数字了,蒙特卡洛的 O(1/√N) 虽然慢,但它是与维度无关的——不管 2 维还是 32 维,都是同一个收敛率。所以"慢"这个标签只是相对低维求积而言,在高维问题上它是唯一实用的选择。
既然收敛率换不了,你能做文章的就只有方差公式里那个 Var[f(X)/p(X)]。后面所有技术——重要性采样、分层、低差异序列、多重重要性采样——本质上都在和这个系数搏斗。
1.2 偏差与方差:两类错误要用完全不同的方法治
在动手优化之前,必须先分清楚渲染结果里的"不对"到底属于哪一类。偏差(bias)是估计子的期望与真值之间的系统性偏移,方差(variance)是估计子围绕期望的波动程度。两者在画面上的表现截然不同。
我见过不少项目组在 debug 阶段浪费了大量时间,就是没分清这两类问题。举一个很典型的场景:有人实现了某种近似的光泽反射采样,渲染出来的高光区域总是"不够亮"或者"形状不对"。他以为是方差问题,一个劲加采样数,画面确实越来越稳定,但稳定地错——高光还是不对。这就是偏差问题,加多少样本都没用,得从算法本身找原因。
反过来,焦散效果"脏兮兮"的、噪点不随采样数增加而均匀下降的,那大概率是方差爆炸。上面第 1.1 节那个公式可以帮助定量理解:如果你某个 pdf 采样策略在某些路径上产生的 f/p 比值大得离谱,少数几个"幸运样本"就能把整体方差拉爆。
判断技巧很简单:先把采样数从 16 提到 1024,如果噪点明显均匀变小,是方差问题;如果画面几乎不变,只是更稳定地错,是偏差问题。第二招是换随机种子:方差问题的两版渲染会有肉眼可见的差异,偏差问题则趋同于同一个错误结果。这个"先定性再动手"的习惯,能让你的优化工作少走一半弯路。
2. 重要性采样:把采样能量花在贡献最大的地方
2.1 方差公式告诉我们的优化方向
上一节给出了优化目标:压低 Var[f(X)/p(X)]。展开写一下:
Var[f/p] = E[(f/p)²] − (E[f/p])² = E[(f/p)²] − I²
第二项 I² 是常数,所以唯一能优化的是第一项。根据柯西-施瓦茨不等式,E[(f/p)²] ≥ (E[f/p])² = I²,等号成立的条件是 p 与 f 成比例,也就是 p(x) = f(x) / ∫f。代入之后,f(X)/p(X) 恒等于 ∫f,方差直接归零。
这个"理想 p"当然是事后诸葛亮——你已经知道积分值才能构造出它。但它的指导意义极强:p 的形状越贴近 f,方差越低。贴近到什么程度算"够好"?实践中并不需要精确匹配,能把 f 的主要峰值轮廓抓住,方差就能降一两个数量级。
渲染场景里的 f 通常是多个因子的乘积,比如直接光照积分里是 BRDF、余弦项、入射辐射度的乘积。你不可能把所有这些因子都塞进 p,但可以找出主导方差的那一个:如果光源很小、BRDF 很宽,主导方差的是光源项,应该按光源采样;如果光源很大、BRDF 很窄,主导方差的是 BRDF,应该按 BRDF 采样。这个"找到主导因素、再构造对应采样策略"的思路,就是所有实用重要性采样的设计起点。
2.2 三个最常用的采样策略:从漫反射到微表面到光源
余弦加权半球采样是最基础的例子。漫反射材质的 BRDF 是常数,被积函数里唯一的强方向依赖就是 cosθ 项(严格说是辐射度与立体角的投影关系)。这时按 p(θ,φ) = cosθ/π 采样半球方向,f/p 里的 cosθ 被消掉,方差会大幅下降,而且实现极其简单。
CDF 反演的过程值得走一遍:先把概率密度写成球坐标形式 p(θ,φ) = cosθ sinθ/π(sinθ 是球坐标的 Jacobian),然后求累积分布函数。对 θ 积分能得到 C(θ) = sin²θ,对 φ 积分得到 C(φ) = φ/(2π)。分别反演得到 θ = arccos(√(1−ξ₁)),φ = 2πξ₂,对应的 z 方向分量等于 cosθ = √(1−ξ₁)。落到代码里可以写成这样:
// 生成局部坐标系下的余弦加权方向,z 轴为法线方向 struct LocalDir { float x, y, z; // 局部方向 float pdf; // 立体角概率密度 }; LocalDir SampleCosineHemisphere(float xi1, float xi2) { float r = sqrt(xi1); // 圆盘采样半径 float phi = 2.0f * PI * xi2; float x = r * cos(phi); float y = r * sin(phi); float z = sqrt(max(0.0f, 1.0f - xi1)); // z = cosθ return {x, y, z, z / PI}; }这就是所谓的 Malley 方法:先均匀采样单位圆盘,再投影到半球上,结果天然服从余弦加权分布。很多渲染器里本来就带着圆盘采样函数,顺手复用就行,不用重新写逆变换。
GGX 微表面模型的采样稍微复杂一点,但原理完全一样。各向同性 GGX 的法线分布函数是:
D(h) = α² / (π · ((cosθₕ)²(α²−1) + 1)²)
其中 θₕ 是微表面法线与宏观法线的夹角,α 是粗糙度参数。要按 D(h)·cosθₕ 采样微表面法线 h,对 CDF 求逆之后得到:
θₕ = arctan(α·√(ξ₁/(1−ξ₁))),φₕ = 2πξ₂
用这个公式采样出 h,再根据入射方向 wi 做一次反射变换 wi→wo,才是完整的出射方向采样。这里有一个特别容易漏的点:上面求的 pdf 是针对微表面法线 h 的 pdf = D(h)·cosθₕ,但估计子里需要的是出射方向 wo 的 pdf。反射变换本身有 Jacobian,转换关系是 p(wo) = D(h)·cosθₕ / (4|wi·h|)。忘了除以这个 4|wi·h|,直接光的强度就会出现一种"高光越强偏差越明显"的奇怪错误,排查起来非常头疼。
光源面积采样完全换了角度:既然积分里有一个因子是光源辐射度,那就直接在光源表面上采样点,然后对着着色点发射 shadow ray。光源上的均匀采样 pdf 是 1/A(A 是光源面积),但积分是定义在立体角上的,涉及面积→立体角的 Jacobian 变换,即 dω = |cosθ'|/r² · dA。于是估计子里要乘上这个修正项。完整做下来,同样的采样数,小面积光源的阴影噪点会下降得极其明显。
姑且分享一个我实际项目里的数据:某个室内场景里有一盏很小的射灯,漫反射地面用均匀半球采样时,噪声主要来自"能不能命中灯"这个二项事件的方差。换成光源面积采样之后,同样的 64 spp,RMSE 大约降了一个数量级。这种性价比是调渲染参数给不了的。
2.3 两个容易翻车的点:PDF归一化与Jacobian缺失
重要性采样实现上的坑,我按遇到频率排前列的两个说。
第一,采样函数里写的 PDF 必须严格归一化。很多人从公式抄来 p(θ,φ) = cosθ/π 就写进代码,但忘了确认它确实是归一化的概率密度。验证方法很简单:在整个定义域上对 p 做数值积分,结果必须等于 1。如果 p 少了归一化系数,估计子里的 1/p 会整体放大或缩小,画面不是偏亮就是偏暗,而且这个错误不会随采样数增加而消失——因为它是有偏的。
第二,凡是做了坐标变换,都要检查有没有漏 Jacobian。上面 GGX 反射的 4|wi·h|、光源面积采样的 cosθ'/r²,还有各种坐标系旋转涉及的 sinθ 因子,都属于这个范畴。这种 bug 最讨厌的地方在于画面往往只是"略亮一点"或者"略暗一点",不仔细对比根本看不出问题。我的排查习惯是:搭一个解析可积的简单场景(比如一个平行光源照射均匀漫反射平面,正好在某个角度下有手算的精确值),然后用渲染器跑这个场景直接对比。偏差一测一个准。
3. 从伪随机到低差异:准蒙特卡洛凭什么更快
3.1 伪随机点集的"聚簇"问题
重要性采样管的是"往哪采样",但即便采样的方向分布完全正确,样本点在积分域内的几何分布质量也直接影响方差。想象你在二维方格里撒点,伪随机数发生器(比如 PCG 或 xorshift)生成的点从统计检验上看是挺均匀的,但把点画出来看,总会有一些区域点密得扎堆,另一些区域空荡荡。这个"聚簇"现象本质上就是随机波动本身,也就是方差的来源之一。
分层采样是第一个觉悟:把积分域划分成若干格子,每个格子里强制放一个随机点,保证不出现"两个点挤一格、另一格空着"的情况。一维和二维下,这个操作能把方差显著压下去。但分层采样有致命弱点:在 d 维空间里,格子数是维度的指数级,路径追踪一出场就是十几维,你根本不可能做到"每格一点"。
那怎么办?答案是用低差异序列——它不追求统计学上的"随机",而是直接追求几何上的均匀覆盖。
3.2 低差异序列三兄弟:Halton、Sobol、Hammersley
低差异序列的核心指标是差异度(discrepancy),直观理解就是:任意一个子区间里,点的实际数量与理想均匀数量之间,最大能偏出多少。差异度越小,点集的覆盖质量越好。三类最常见的序列:
Halton 序列用一组互质的素数(2、3、5……)做底数,第 i 个点的第 d 维取 i 在对应底数下的基数逆(radical inverse)表示。它的优点是构造简单、不需要预计算,二维三四维下效果很好,很适合学习验证和做测试用例。
Sobol 序列基于二进制方向数(direction numbers)构造,在高维下的分布质量通常优于 Halton,也是现代学术研究和工业渲染器里最主流的选择。常见的开源随机数库里基本都有现成的实现,拿来直接用就好。
Hammersley 点集的形式是 (i/N, φ₂(i), φ₃(i), …),它依赖总点数 N,适合静态、固定预算的场景。路径追踪里每像素要大量样本,动态图像上它对 N 的变化不友好,所以更多用在某一类固定采样数的工具链里。
| 序列 | 构造基础 | 优点 | 适用场景 |
|---|---|---|---|
| Halton | 互质素数基数逆 | 实现简单、无预计算 | 学习验证、低维测试 |
| Sobol | 二进制方向数 | 高维质量好、主流选择 | 路径追踪实际渲染 |
| Hammersley | 依赖总点数 N | 差异度优秀、形式简洁 | 固定采样数的静态问题 |
三者有一个共同点:渐近差异度大约在 O((log N)^d / N) 量级,这比伪随机数的期望差异度 O(1/√N) 要快得多。这就是"免费获得一个维度优势"的来源。
3.3 Koksma-Hlawka 定理:准蒙特卡洛收敛优势的理论根基
如果说蒙特卡洛靠的是大数定律和中心极限定理,那么准蒙特卡洛的精气神就在 Koksma-Hlawka 定理里。这个定理把 QMC 的误差拆成了两个独立项的乘积:
| (1/N)·Σf(xᵢ) − ∫f dx | ≤ V(f) · D*(x₁,…,xN)
左边是 QMC 估计与真值的误差,右边第一项 V(f) 是被积函数在 Hardy-Krause 意义下的全变差,它刻画了函数有多"陡峭";第二项 D* 是点集的星差异。
这个公式的含义很朴素:误差的上界由"函数自身的复杂程度"和"你选的采样点集的均匀程度"共同决定。函数复杂你管不了(那是场景问题),但点集的差异度是可以通过换序列直接优化的。伪随机点的星差异期望约 1/√N,于是误差上界也按这个速率衰减;低差异序列的星差异是 O((log N)^d/N),N 的指数从 1/2 变成了 1,收敛速度就上了一个台阶。
需要提醒一句:Koksma-Hlawka 给的是确定性误差上界,实践中 QMC 的优势未必总是那么理想,尤其是维度升高之后,log N 的幂次会侵蚀一点优势。但工程实测里,随机化 QMC 在绝大多数渲染场景下都优于纯伪随机,这一点几乎没什么争议。
3.4 路径追踪里怎么落地:维度分配与序列随机化
路径追踪的积分域远不止二维:一条长度为 k 的路径,需要 2k 个随机数来描述路径上每个顶点处的出射方向。真把全部随机数扔进一条高维低差异序列里,也不是不行,但效果往往不如"分而治之"。常见做法是:
每个像素维护一个采样索引 i,第一次弹射的方向用二维低差异序列的点(比如 Sobol 的点对)来生成,相当于在像素内做了高质量的分层;后续弹射方向再从序列的高维段去取。这样安排之后,像素内的子样本覆盖质量最好,而路径更深处的随机数质量下降对最终图像的影响相对较小。
还有一个工程上必须做的步骤:序列随机化。纯确定性的低差异序列渲染出来的图像,噪点不是随机颗粒感的,而是带有某种规律性的条纹或网格纹理,动画里还会表现为"爬动"的伪影。解决办法是给序列加一个随机偏移,最经典的是 Cranley-Patterson 旋转:每个样本点加上一个均匀随机向量,再对小数部分取模。这样既保留了低差异的均匀覆盖性质,又让每个像素或每帧拥有独立的随机扰动。
实测下来,我的路径追踪器从纯伪随机换成随机化 Sobol 序列之后,同采样数下噪点大约降了三到五成,这个提升不需要改任何着色代码,纯粹是采样层面的白捡收益,优先级非常高。
4. 多重重要性采样:把多种策略做成组合拳
4.1 单一采样策略的天花板:一个必然顾此失彼的困境
直接光照积分的被积函数可以粗略拆成两部分:BRDF 因子和光源因子。假设场景里有一个极小的点光源,BRDF 却是很宽的漫反射——按 BRDF 采样方向,打到光源上的概率微乎其微,绝大多数样本都落在"没光照"的区域,方差巨大。反过来,如果光源又大又均匀,BRDF 却是窄窄的镜面反射,按光源采样生成的入射方向里能命中镜面波瓣的也少得可怜。
这就是单一策略的宿命:它只能匹配一个因子的形状,必然在某些场景下被另一个因子压制。选择哪一个,本质上是在"大概率命中光源"和"大概率命中 BRDF 波瓣"之间做取舍,但被积函数是两者的乘积,只匹配一边就迟早翻车。
4.2 MIS 估计子的形式与无偏条件
多重重要性采样的思路朴素得让人意外:与其二选一,不如把几种策略都跑一遍,然后用权重把结果组合起来。假设有 m 种采样策略,每种策略 i 产生 Nᵢ 个样本,MIS 估计子写成:
F = Σᵢ (1/Nᵢ) · Σⱼ wᵢ(Xᵢ,ⱼ) · f(Xᵢ,ⱼ) / pᵢ(Xᵢ,ⱼ)
其中 wᵢ 是策略 i 的样本在点 x 处被分配的权重,要求对任意 x,所有权重加起来等于 1(严格说,是所有采样到 x 的概率非零的策略权重之和为 1)。另一个关键条件是"兼容性":如果某个策略在 x 处的 pdf 为零,而另一个策略恰好采到了 x,那么前者的权重必须为零,否则分母为零、无偏性直接崩溃。
MIS 的核心定理是:只要满足上面两个条件,估计子就是无偏的。这个自由度给了我们极大的设计空间——权重不用调参,公式给完就能用,这也是它成为现代渲染器标配的直接原因。
4.3 平衡启发式与幂启发式:两个最常用的权重公式
最常用的是平衡启发式:
wᵢ(x) = pᵢ(x) / Σⱼ pⱼ(x)
直观理解:在某个位置 x,哪个策略的采样概率密度最高,它的样本就越要被保留。这个公式是对称的,自动在策略之间做调和,几乎不需要人工干预。代码极其直白:
float BalanceHeuristic(float pdfA, float pdfB) { return pdfA / (pdfA + pdfB); }第二种是幂启发式,通常取 β=2:
wᵢ(x) = pᵢ²(x) / Σⱼ pⱼ²(x)
幂启发式对"各策略 pdf 形态差异巨大"的情形更激进:它会更狠地压低那些落在另一个策略强分布区域的样本的权重,从而抑制异常大值对结果的冲击。实践里,它的稳定性通常略优于平衡启发式,所以我现在的默认做法是幂启发式,特殊场景下再换回平衡启发式对比。
注意:权重公式里的 pᵢ 应该严格写成 Nᵢpᵢ——当各策略样本数不同时,严格形式是 wᵢ(x) = Nᵢpᵢ(x) / Σₖ Nₖpₖ(x)。如果两组策略采样数差距大而你又忽略了 Nᵢ,组合结果会出现细小的偏差,肉眼未必看得出来,但定量对比时一定会露馅。
4.4 双向路径追踪为什么离不开 MIS
双向路径追踪(Bidirectional Path Tracing)就是把"从相机出发的路径前缀"和"从光源出发的路径后缀"在中间某个顶点处连接起来,每条完整路径对应一种连接方案,而每种连接方案实际上就是一条独立的采样策略。难就难在:一条焦散路径,从光源路径那一端看过去,它的概率密度极小;但从相机路径那一端看过去,概率密度又非常自然。没有 MIS 的话,你只能把每一条路径硬塞给其中一个策略来估计,那结果就像把两个问题简单叠加在一起——方差照样爆炸。
有了 MIS 之后,每一个连接策略都根据自身与对方的 pdf 获得权重,两侧测量自动找到了最优组合。这就是为什么现代双向路径追踪的焦散渲染比单向路径追踪干净那么多的底层原因。
实现 MIS 最需要纪律的地方是记录 pdf:路径构建过程中每个顶点都要保存该顶点的出射方向在各采样策略下的 pdf,而且要跟着坐标变换同步更新。如果路径构建时发生了某种方向压缩变换,变换后的 pdf 必须跟着变,否则 MIS 权重算出来是错的。这个"多存几个 float"的代价微不足道,换来的质量提升绝对值回票价。
5. 工程实录:调试采样与方差的个人经验
5.1 按出现频率排名的三个 Bug
第一个是 PDF 不归一化。任何采样函数写完,第一件事就是验算它的 pdf 在整个定义域上的积分。我习惯在开发环境里写一个数值积分函数,采样网格上对 p 做累加,跑一次单元测试就知道有没有归一化问题。这比盯着公式发呆快得多。
第二个是估计子里的 1/p 与实际采样分布不匹配。改动采样策略时,很容易出现"采样函数换成了新分布、估计子里的 1/p 还写着旧分布"的错位。解决方案是把采样函数和 pdf 计算封装在同一个结构体里,成对出现、成对修改,物理上就杜绝了分离导致的错位。
第三个是漏 Jacobian。GGX 反射的 4|wi·h|、光源面积采样的 |cosθ'|/r²、球坐标系变换的 sinθ,这些都是经典重灾区。这类 bug 最好的排查手段是构造一个解析可积的验证场景,直接用数值积分基准对比渲染结果,偏差立刻暴露,不用瞎猜。
5.2 三步定位:方差问题还是偏差问题
判断阶段我总结了一套三步法,实测非常管用:
- 把每像素采样数从 16 提到 1024。噪点明显均匀变小的,大概率是方差问题;画面几乎不变或者只是更稳定地错,需要考虑偏差。
- 换一个随机种子重新渲染。方差问题下两版图像有明显差异,偏差问题下两版趋近于同一个错误结果。
- 渲染一个解析解已知的简单场景,直接算 RMSE。这能锁定问题到底是采样策略不对、权重公式错了,还是更底层的物理假设有问题。
关于参考图,我想多说一句。很多人习惯用高采样数的路径追踪结果当"金标准"参考图,但要注意:如果参考图本身用的采样策略和你在测的策略不同,两者的偏差来源也不同,拿它当衡量尺子反而会误伤。我会维护两套参考:一套是 4096 或 8192 采样数的同策略结果,用来评估方差降低的幅度;另一套是数学上有闭合解的简化场景,用来评估绝对精度。两套配合,定位问题非常高效。
5.3 封装采样接口的一点纪律性做法
这几个基础工具——重要性采样、低差异序列、MIS——任何一个单独拎出来都不算难,真正的难点在于组合。你在同一个路径追踪器里同时实现余弦加权采样、GGX 采样、光源采样和 MIS,又换了 Sobol 序列,那么每一层的改动都必须保证新的 pdf 与旧的估计子兼容,一个环节脱节,整条链路的方差管理就崩了。
我个人的工程习惯是:把"采样方向、返回 pdf、计算逆 CDF"三件事封装成一组不可分割的接口,并在数据结构里显式记录"这个样本来自哪个策略",一路传递到 MIS 权重计算处,不靠隐式约定。另外,在调试模式下我会给每个策略的样本加一个 tag,渲染结果里可以按策略拆分统计贡献量。这个输出在排查"哪个策略把方差拉爆了"的时候特别有用——有时候根本不是你以为的那个策略在出问题。
最后说一个长期养成的习惯:每次往渲染器里加新的采样策略,我都会保留旧策略的完整实现,在同一个测试场景里并列渲染对比。这看起来笨,其实最省时间。加 QMC 那一次,我以为改动很小不需要对比,结果一个维度分配的错误让我花了整整两天才找到。从那以后,对比测试就是我实现任何采样相关功能的固定流程了。这个系列如果继续写下去,下一步我会展开路径空间公式和更高级的零方差技巧,那些内容全建立在本文这套基础上,先把今天这些吃透再前进。