望蓟门
[唐代][祖咏]
燕台一望客心惊,箫鼓喧喧汉将营。
万里寒光生积雪,三边曙色动危旌。
沙场烽火连胡月,海畔云山拥蓟城。
少小虽非投笔吏,论功还欲请长缨。

本文是对Intel关于户外大气散射实现的文章Rendering Outdoor Light Scattering in Real Time（作者为Naty Hoffman等）的翻译，原文请自行搜索。

Introduction

大气散射会对太阳的颜色以及天空的颜色造成影响，且这个影响过程并不是静态的，而是随着24小时动态变化的。另外，大气散射的表现还跟天气，空气状况等因素相关，对于具有不同大气成分的星球的天空表现可能跟地球的表现是不一样的。

本文将会介绍大气对于天空光照等表现的理论知识，并给出传统雾效模型在模拟大气散射表现上的不足，并给出一种PBR的实现方式。

Interaction of Light with Particles

光照跟粒子的交互是图形学中的基本事件，实际上通常我们所处理的光照与物体表面的交互可以看成是光照与粒子交互的特例。光线作为一种电磁波在于粒子发生交互的时候，会出现以下几种情况：光线被吸收；光线被反射；粒子自发光（本文不予考虑）。

Absorption

粒子对于光线的吸收可以用cross section（交叉区域面积） $\sigma^\lambda_{ab}$ 上的吸收量来表示，其中 $\lambda$ 表示的光波波长，表示这个数值是跟光波波长相关的参数，其单位为 $m^2$ 。具体而言指的是单位入射光照度（incident irradiance，单位为 $W/m^2$ ）吸收的光辐射通量（单位为 $W$ ），光辐射照度通过光辐射通量与表面面积相除得到，想要得到单位光辐射照度，就需要根据光辐射强度调整承光表面倾斜角度，此时包围粒子的包围球在这个角度的横截面投射到与光照方向垂直的位置的面积，即可理解为cross section面积。从物理含义来说，这个面积可以看成是粒子吸收光线的吸收率。

出于叙述方便，这里假设粒子将打在包裹粒子的球形区域的光线全部吸收了（球形区域可能比粒子本身体积大或者小），参考上图中这个球的cross section。在cross section上的每一点，都会进行光线的吸收作用，因此整个粒子的吸收的光辐射通量就是这个section上每个点的入射光辐射照度的积分。

一种吸收性的介质由一定的密度的粒子组成，粒子的密度用 $\rho_{ab}$ 表示，其单位为 $m^-~^3$ ，每个粒子的吸收cross section为 $\sigma^\lambda_{ab}$ 。这里定义此介质的吸收系数 $\beta^\lambda_{ab}$ = $\rho_{ab}$ $*$ $\sigma^\lambda_{ab}$ ，其单位为 $m^-~^1$ 。这个系数表示的物理意义为，每单位体积的粒子所对应的cross-section面积。为了给出这个参数的直观意义，以下图作为示意：

假设以固定的方向朝着上述介质射出光线：在光线垂直方向上定义一个横截面积为A的盒子，这个盒子的厚度为dx，那么其体积为A*dx。按照上面的定义，这个盒子的累计absorption cross section面积为 $\beta^\lambda_{ab} * Adx$ 。而如果沿着光照方向随机发射一个光子，其被盒子所吸收的概率可以用累计吸收交叉面积除以盒子横截面积A来计算，即 $\beta^\lambda_{ab} * dx$ （这个结果是没有单位的）。按照这个结论，可以推导出，强度为 $L^\lambda$ 的光线穿过厚度为 $dx$ 的介质后，被吸收的光强大概为其本身强度的 $\beta^\lambda_{ab} * dx$ ，用公式表示可以得到如下结果：

公式1

对公式1进行解析，得到

公式2

在上述公式中0表示的是光线入射路径上某一起始评估点， $L_0^\lambda$ 则是此点的光强，x表示的是从起始评估点继续往前行走的距离。

公式2中将 $\beta^\lambda_{ab}$ 看成是常量来计算，但实际情况中有可能这个数值是一个随着空间位置而变化的变量，此时得到如下的公式：

公式3

另外，假设传播介质中存在n种不同的吸收粒子，每种吸收粒子有着自己不同的吸收系数 $\beta_{ab(i)}^\lambda$ ，那么总的吸收系数可以表示为多种吸收粒子吸收系数之和：

公式4

吸收系数之所以会跟波长有关，是由于粒子本身的吸收特性所导致的，在实现的时候可以使用经验数据来实现对这些数据的模拟，从[Preetham1999a]中可以找到所有大气散射常见粒子的一些相关数据。

Scattering （Out-Scattering）

与吸收作用一样，粒子对于光线的散射作用也可以用一个类似的散射cross section $\sigma^{(\lambda)}_{sc}$ 来表述，这个cross section的定义是：单位入射irradiance情况下的散射radiant通量。

Fig 3. 散射cross section

如图3所示，将所有的输入irradiance积分起来得到的结果就是所有的散射radiant通量。仿照前面一节的处理方式，引入散射系数： $\beta^{(\lambda)}_{sc} = \rho_{sc} * \sigma^{(\lambda)}_{sc}$ 。这个系数可以看成是单位体积下的散射cross section，在这个定义下，radiance值为 $L^{(\lambda)}$ 的入射光在经过 $dx$ 的传播距离后，由于散射（out-scattering）导致的损耗radiance $dL^{(\lambda)}$ 与散射系数之间存在如下关系：

公式5

对上述公式进行解析，可以得到公式6

公式6

如果考虑散射系数随着空间位置而变化的情况，公式6可以转化为更完备的公式7：

公式7

跟absorption一样，总体的散射系数等于各种不同的粒子类型的散射系数之和：

公式8

如果将散射跟吸收的作用放到一起的话，可以将二者的损耗系数加到一起，得到总的损耗系数：

公式9

Scattering (In-Scattering)

前面考虑的out-scattering对于光照的作用于absorption类似，不过这里还少考虑了一件事情，经过吸收后的光损可以直接忽略不予考虑，但是out-scattering散射出去的光线却不能直接忽略，这部分能量可能会窜入其他光线的作用范围，成为这里需要考虑的in-scattering散射。

前面介绍了散射系数，其作用是给出光线在传播过程中的散射损耗，但是却没有给出散射损耗在各个方向上的占比，这里定义一个相函数用来描述这个数值： $\Phi(\theta, \phi)$ 。相函数会给出散射后的光强在各个方向 $(\theta, \phi)$ 上的分布比例（或者说单束散射后的方向的分布概率）。

相函数的概率分布是以单位固体角为刻度来进行衡量的，其单位为 $sr^{-1}$ 。如果散射发生位置的粒子时球状的，或者粒子半径相对于光波波长而言可以忽略的话（在当前考虑的散射中，是满足这个特征的），那么就可以将光线与粒子的作用简化为二维平面交互，从而忽略以光线传播方向与粒子的交点为圆心的360度角度 $\phi$ 的变化，在这种情况下，相函数可以简化为 $\Phi(\theta)$ 。

由于相函数是一个概率分布函数，那么就可以得到如下的公式：

公式10

而为了保证最终的相函数能够具备上述特性，可能需要对相函数进行归一化处理。对于各向同性的情况而言，归一化因子等于 $1/4\pi$ 。另外需要注意的是，这里给出的相函数没有考虑波长的影响，而这个特征并不是对于所有相函数都适用，部分相函数需要考虑波长对于散射分布的影响，而本文后续选择的相函数则是跟波长无关的。

综上，我们可以得到跟散射角度相关的散射系数分布规律： $\beta^{(\lambda)}_{sc}(\theta) = \beta^{(\lambda)}_{sc} * \Phi(\theta)$ ，其中：

公式11

$\beta^{(\lambda)}_{sc}(\theta)$ 的单位为 $m^{-1}sr^{-1}$ .

Fig 4

如图4所示，入射光固体角为 $d\omega$ ，入射光方向与观察方向夹角 $\theta$ ，那么根据前面的描述，当前从入射光方向散射进入观察者眼睛中的in-scattering光概率等于 $\Phi(\theta)*d\omega$ ，将概率乘以入射光Radiance $L^{(\lambda)}_i(\theta,\phi)$ 得到inscattering radiance $L^{(\lambda)}_i(\theta,\phi)*\Phi(\theta)*d\omega$ （这里的 $d\omega$ 用于限定入射光方向，只有这个方向的入射光，其相函数对应的角度为 $\theta$ ），这个公式也可以看成是 $L^{(\lambda)}_i(\theta,\phi)*d\omega$ 跟 $\Phi(\theta)$ 相乘，而前者就等于 $E_{(d\omega)}$ 。对这个公式进行球面积分（为什么要进行球面积分，一个给定的输入光，其方向不是固定的吗？从下面的公式看出，这里的 $L_i$ 指的是inscattering部分的光照数据而非第i盏灯的光照数据，而这里的球面积分指的是从四面八方各个方向传过来的inscattering光照在相函数作用下沿着视线方向散射的结果的累加，因此这个积分是有必要的），就得到了 $L_i$ 在当前粒子上发生的总的in-scattering radiance：

公式12

上述公式只是单点散射导致的in-scattering radiance，而如果考虑入射光线路径上极小距离 $dx$ 的in-scattering radiance，就需要乘上散射系数以及传播距离（这个数值跟out-scattering中的散射系数一样），从而得到由于in-scattering而导致的radiance增量为：

公式13

将损耗跟in-scattering增益结合起来，得到差分公式：

公式14

对上述公式进行解析，得到公式15：

公式15

其中 $0$ 指的是光线与传播媒介相接触的第一个点（通常都是发光点所在位置）， $x$ 指的是从第一个点到人眼这段距离。公式15看起来有点复杂，但是在实际使用中会有很多简化技巧，总的来说，光线传播过程中有两项作用，一项是损耗，一项是增益：

公式16

Scattering Coefficients

下面介绍一下相函数的一些知识。相函数的表现也就是形状取决于散射媒介（也就是说不同的粒子对应于不同的相函数，那么瑞利散射跟米氏散射的相函数应该是不同的），对于散射粒子半径较小（ $r < 0.05\lambda$ ）的情况，相函数 $\Phi(\theta)$ 可以用瑞利散射中的函数来描述：

公式17

注意，瑞利散射中的相函数在光线传播前后方向上的散射分布是对称的。

瑞利散射的散射系数具体公式在[Preetham1999a]中可以找到，这里仅仅关心的是其散射系数跟波长的四次方成反比（散射效果之所以受到波长的影响，是否也是因为瑞利散射本身是由于小尺寸粒子所导致，因此更容易受波长的影响呢？），因此短波比长波容易散射（天空呈现蓝色与傍晚天空呈现红色的原因）

对于尺寸更大一点的散射粒子，瑞利散射模型就不再适用了，此时需要修正为米氏散射。米氏散射的公式比较复杂，不过这里可以做一些简化处理，从而可以使用 Henyey / Greenstein 相函数 $\Phi_{HG}(\theta)$ 来近似模拟米氏散射的相函数：

公式18

公式18实际上是椭圆方程的极坐标形式，其中g是奇异点。实际使用中，如果g取负值会使得散射光线更多的分布在传播方向的前方范围，而取正值则会导致散射光线更多的分布在传播方向的反方向。大多数粒子对应的g值都是小于0的，且随着粒子尺寸的增大，g值的绝对值也会增大。

实现米氏散射模拟的另一个方法就是LUT，根据参数对预先存储的米氏散射测量结果进行取用即可。

Fig 5

图5给出了散射效率 $Q_{sc}$ 在无任何吸收作用的粒子媒介中随着半径与波长之比 $r/\lambda$ 的变化规律，其中散射效率可以看成是散射cross section与粒子的几何cross section之比： $Q_{sc} = \sigma^{(\lambda)}_{sc} / (2*\pi*r^2)$ （物理意义可以看成是单位横截面积下的散射cross section，也可以理解为粒子的反射能力，从图中可以看到，固定输入光波长，在半径与波长相当的时候，散射能力是最强的，之后就进入震荡并最终收敛下来）。

对于细小粒子而言， $Q_{sc}$ 的变化曲线接近于 $r/\lambda$ 的二次方，加上 $Q_{sc}$ 定义中的几何cross section中的 $r^2$ 项，那么 $\sigma^{(\lambda)}_{sc}$ 的变化规律就跟 $r^4 / \lambda^2$ 保持一致，如果固定住r的话，那么这个变化规律就跟 $1/\lambda^2$ 保持一致，怎么跟瑞利散射的 $1/\lambda^4$ 不一致呢？

随着粒子半径的增加， $Q_{sc}$ 对波长的依赖关系逐渐减弱直到 $r/\lambda$ 接近于1为止。从这之后， $Q_{sc}$ 开始随着波长的增加而增加，曲线表现为震荡，直到最终逼近一个极限值2。对于大尺寸粒子而言，其散射表现跟波长的关系基本上可以忽略。米氏散射中的介质粒子，其尺寸通常不会固定为单一的半径数值，而是在一个范围连续分布的，在这种情况下，各个粒子的震荡输出可以直接使用极限值2来进行拟合，这种情况下就可以直接忽略波长的作用。

Earth’s Atmosphere

地球大气中主要有两类粒子，其中一类是小尺寸气体分子（gas molecules），主要在干净的天空中出现，这类粒子是瑞利散射的主要贡献者；另一类是较大尺寸的气溶胶粒子（aerosols），这类粒子是米氏散射的主要贡献者。部分气溶胶粒子还会带有一定程度的吸收作用，且由于气溶胶粒子的半径范围很广，因此前面用2来拟合 $Q_{sc}$ 的假设是成立的。

Applying the Theory - Aerial Perspective

Fig 6

在现实世界中，如果我们观察远方的物体，其表面颜色会随着距离以及天气的变化而变化，这是由于大气散射的作用导致。如图6所示，相机观察方向面向的是一座山，这座山的表面颜色沿着观察方向前进的过程中会由于大气散射的影响而出现out-scattering以及absorption等原因导致的损耗，同时也会由于太阳光在空气中传播而导致的in-scattering输入。如果当前空气比较干净的话，那么此时散射的主导者是瑞利散射，这会导致大气散射中占据主导地位的是短波也就是蓝色光，也就使得大气散射对于入射光线的损耗系数 $\beta^{(\lambda)}_{ex}$ 是偏红色的，而在in-scattering方向偏离太阳光方向较大时，in-scattering是偏蓝色的（瑞利散射），当两者方向接近时，此时会是无波长影响的米氏散射占据主导（为什么？瑞利散射的强度与波长四次方呈反比，意味着蓝光散射概率大，其他光概率相对低，也就意味着其他光维持原方向前进的概率大，而根据瑞利散射的相函数，我们知道，在夹角较小的时候，散射光的占比更高，那么这个时候的光其实更接近于全色光，也就表现为波长影响不大，对应于米氏散射），从而呈现亮白色。

在游戏中，大气散射通常是使用fog的方式来模拟的，其计算公式为： $L^{(\lambda)} = L^{(\lambda)}_0(1-f)+C_{fog}f$ ，总结起来就是，经过大气散射后的物体的颜色等于未经过散射前的颜色与雾气的颜色通过一个权重f进行混合后的结果。不过呢，这种近似模拟的效果显然跟实际表现相去甚远：首先混合的系数f是一个跟波长无关的数值，这跟实际情况中散射与波长的关联关系不一致；其次输出的结果跟观察方向无关，这也与此前的分析相冲突。

公式15

正确的使用方式应该是以公式15为蓝本，根据实际情况加以简化。而这些简化中最常用的一个就是对多种不同散射介质的密度函数的简化处理。通常情况下，相机跟待观测的物体本身的海拔都是很低的，因此这里可以做出一些简化处理（具体可以参考后面的Implementation Note）。

Applying the Theory - Sunlight

太阳光从太空中向着地球表面传播的过程中会发生一些损耗，根据公式15中的损耗计算公式，可以由太阳的位置以及空气的大气散射系数计算出任意点处的太阳辐射radiance，随着太阳位置的下降，其在大气中传播的距离也相应增加，导致瑞利散射中被散射出去的蓝光比例也相应增加，最终使得早晨以及黄昏的时候天空呈现橘红色。

太阳光的传播虽然看起来像是一种大气散射的效应，但实际上需要考虑的东西比大气散射要多一些，比如至少需要考虑大气层的弯曲形状，因此在实际计算中的时候会需要更为精准的大气密度函数。

Applying the Theory - Sky Color

天空颜色是in-scattering的结果，在天空澄澈的时候，对短波友好的瑞利散射使得天空中分布着大量的蓝光从而呈现蓝色；而在空气质量较差的雾霾天气或者阴雨天气时，此时占据主导的是米氏散射，因此呈现银灰色。此外夕阳西下时，1.太阳周边的一圈红色的染色也是米氏散射的作用（如何解释？短波因为长距传输过程中的高额散射消耗而损耗，只剩下长波还沿着原始的方向传播着，但是这不应该是瑞利散射的作用吗？毕竟米氏散射是跟波长无关的）；2.与太阳相对的天空会由于大气的荫蔽作用以及in-scattering的减少而显得较暗。

对于天空而言， $L^{(\lambda)}_0$ 通常为0（纯黑），因此前面用于计算普通物件经过大气散射后的颜色的传播损耗就可以直接移除了，这也是天空颜色渲染中与普通物件渲染不一样的地方，因此如前所述，需要更为精准的大气密度函数。

Conclusion

本文对大气散射的一些基本知识进行了简要的介绍，并给出了户外大气散射计算的基本公式。

Implementation Note

如果我们假设大气密度是恒定的，且相机与观察对象都处于一个低海拔的状态时，可以对原计算公式15进行简化处理，得到的公式可以用满足实时渲染的需要，公式给出如下：

公式19

这个公式可以看成是物体经过大气散射out-scattering之后的radiance加上太阳光经过in-scattering之后的radiance。其中 $L_0$ 指的是物体原有的颜色， $s$ 指的是物体距离相机的距离， $\theta$ 则是对应于当前像素的观察方向与太阳光方向的夹角， $F_{ex}(s)$ 与 $L_{in}(s,\theta)$ 的计算公式给出如下：

散射后剩余比例等于瑞利散射与米氏散射的乘积，可以理解。

$1 - e^{-(\beta_R+\beta_M)s}$ 可以看成是散射比例，而被观察点处的太阳光irradiance $E_{sun}$ 前面的系数可以看成是散射部分中沿着角度 $\theta$ 方向的比例。很显然，这里只考虑了太阳光在被观察位置的单点in-scattering，没有考虑整个观察方向路径上其他点的in-scattering影响。

散射系数

References – for Further Reading

[Blinn1982] J. F. Blinn. Light Reflection Functions for Simulation of Clouds and Dusty Surfaces.
Computer Graphics, 16(3): 21-29, July 1982.
[Dutré2001] P. Dutré. Global Illumination Compendium. Available online at
http://www.cs.kuleuven.ac.be/~phil/GI/
[Henyey1941] L. G. Henyey and J. L. Greenstein. Diffuse Reflection in the Galaxy. Astrophys.
J. 93, 70, 1941.
[Hoffman2001] N. Hoffman and K. J. Mitchell. Photorealistic Terrain Lighting in Real Time.
Game Developer, 8(7): 32-41, July 2001.
[Klassen1987] R. V. Klassen. Modeling the Effect of the Atmosphere on Light. ACM
Transactions on Graphics, 6(3): 215-237, July 1987.
[Mie1908] G. Mie. Bietage zur Optik truber Medien Speziell Kolloidaler Metallosungen.
Annallen der Physik 25(3): 377, 1908.
[Preetham1999a] A. J. Preetham, P. Shirley, B. E. Smits. A Practical Analytic Model for
Daylight. Computer Graphics (Proceedings of SIGGRAPH 1999): 91-100, August 1999.
[Preetham1999b] A. J. Preetham, A Practical Analytic Model for Daylight. M.Sc. thesis,
Department of Computer Science, University of Utah. March 1999.
[Rayleigh1871] J. W. Strutt (Lord Rayleigh). On the light from the sky, its polarization and
colour. Philos. Mag. 41: 107-120, 274-279, April 1871.
[Yee2002] H. Yee, P. Dutré, S. Pattanaik. Fundamentals of Lighting and Perception: The
Rendering of Physically Accurate Images. Proceedings of Game Developer Conference,
March 2002.

【翻译】Rendering Outdoor Light Scattering in Real Time