在Unity中制作Procedural Sky (Part 1-Atmospheric …

RedZero9 · 发表于 2021-12-24 10:09

我想给自己的钓鱼游戏加上日夜变化，如果游戏中太阳会升起落下，月亮会阴晴圆缺，不可避免，游戏中的天空也要像真实的那样。
现实中，在没有雾霾的天气里，我们能够看到太阳升起时的鱼肚白，能够看到白天时湛蓝的天空，也能够看到日落时布满天边的晚霞。

感谢Unity，我们甚至不需要修改代码，就可以在引擎中看到这样的天空，这个要归功于Unity内置了Shader <Skybox/Procedural>，只要在游戏中使用基于这个材质的天空球，就能看到拟真的天空。详细操作可以看下下面这个视频。
引擎中的效果是这样的
日出时

白天时

这个内置天空球基本能够满足我的需求，但是它有个不足，就是它对夜晚效果的支持基本为零，当太阳下山之后，天空中除了黑色就是黑色，这样可不行，夜晚的天空即使没有月亮，起码也要有星星的嘛。要想支持星星还有月亮的渲染，这时候就需要修改默认的Shader <Skybox/Procedural>。不过在修改原始代码之前，我们还是需要理解一下默认Unity代码是怎么实现的，只有了解原来的逻辑，才有可能在这上面做二次开发，加上其他的效果。
下面文章会分两篇，这一篇会解释一下真实天空从物理上是怎么解释的，以及是如何计算的，下一篇会聊一下Unity默认的程序化天空球是怎么模拟的。
<hr/>1. 原理解释

如果渲染一个房间中的物体，我们常常假定光是真空中传播的，这意味着光在传播的过程中不会被散射吸收，这个假定对于光在空气中短距离的传播，是可以接受的，因为空气密度是比较低的，传播距离相对较短，空气分子对光的散射和吸收基本上可以忽略不计。但是如果将光的传播扩大到整个大气层范围，那么大气对光的散射和吸收，就需要纳入考量了。
对于地球和大气，我们先做一些物理上的简化

假定地球是圆的，半径为6378.1公里。（实际上地球并非完美的球体，地球赤道的半径为6378.1公里）
假定地球的大气层高度为160公里。（实际上大气层没有明确的上界，由于地球引力作用，几乎全部的气体集中在离地面100公里的热层、其中99%在低于25～30公里以内，一般都以100公里这个高度（卡门线）当成外太空的起点，这里设置为160公里是根据unity <Skybox/Procedural> shader中的外圈半径(1.025-1) * 6378.1 = 160算出来的）

这里我们先假定摄像机在外太空，然后去观察地球和地球大气层。

图中地球和大气层的半径并非真实比例

上图中摄像机的视线沿着直线AB通过大气层，因为光是直线传播的（爱因斯坦的广义相对论论证光线在通过强引力场附近时会发生弯曲，这里不做考虑），沿着BA方向的光线最后都有可能到达摄像机，然后在摄像机中形成看到的图像。

由神舟十二号拍摄的太空照片

从上图可以看到地球外面有一层薄薄的蓝色半透明大气层，大气正常是无色的，太阳光一般来说白色偏黄一些，但是为什么太阳光通过大气层到达摄像机就偏蓝色呢？
这个就需要先分析一下太阳光通过大气层都发生了什么事情。
有些化学物质会吸收光能，比如说臭氧层会吸收紫外线，只是紫外线不在光的可见光谱范围之内，因此对最后摄像机看到的效果没有影响，计算时常常不做考虑（大气中的水蒸气和二氧化碳也会吸收光线，不过主要是波长2m以上的远红外线，也不在可见光谱范围）。
大气层对光传播的影响主要在于散射。光线在传播过程中，碰到空气中的分子，会不断被散射到其他方向上，如下图，光传播的距离越远，光的强度也会越小。

如果要计算最后到达摄像机（人眼）的光线，我们主要考虑两种散射过程：

Out-scattering：即摄像机在光线传播方向上，但是光线被空气分子转向(deflect)了，并没有射到摄像机里。

2. In-scattering: 即摄像机不在光线传播方向上，光线被空气分子转向(deflect)后，刚好射到摄像机里。

因此最后摄像机能够接收到的光线如下图

可以看到，摄像机最后不光可以接收到direct light，也能接收到indirect light。上图对光的转向(deflect)次数做了简化，只转向了一次，实际上光线在大气层中可能会被转向(deflect)任意的次数，这意味着真实到达摄像机的光的路径是非常复杂的。如果要模拟非常逼真的光线效果，就需要模拟每条光线的真实路径，这个技术叫raytracing，即光线追踪，支持这项技术的所需计算量对于实时渲染是非常昂贵的，目前台式机N卡20系列以后的显卡已经支持实时光线追踪。
我们这里只考虑光线的单次转向，即single scattering（单次散射）,因为做这样的简化，也能得到相对真实的效果。
总结一下，让我们回到最开始的那张地球示意图，从太阳出来的光线经过大气层到达摄像机的过程应该是下图这样的。

因为相对于地球比较远，上图中的太阳光可以当作directional light处理。

2. 公式推导

下面我们要基于前面的示意图开始计算到达摄像机的光的数量。
透射函数(The Transmittance Function)

我们这里会做一个假定，即太阳光经过真空到达地球大气层之前没有衰减，每条到达地球大气的光的量都是一样的，即  。
光线在大气层中传播一定距离后，到达P点，光的量为  。

这条光线从C点传播到P的过程中，不断被散射，即Out-scattering。如果我们能够算出这个过程中有多少比例的光被Out-scattering，那么就能算出。计算这个比例的函数由  来表示。这个函数后半部分再推导。

散射函数(The Scattering Function)

当光线从C点传播到P点，然后有多少比例的光会偏折(deflect)到PA方向，这个是由散射函数  来表示。其中  表示光的波长， $\theta$ 表示偏折的角度，h表示P点距离地面的高度。这里之所以会有高度，是因为大气不同的高度密度是不同的，大气的密度会影响光的散射。注意光的波长也会影响散射。

如果知道散射函数，我们就能算出朝向A的光数量，这条光线再经过PA这段的距离传播Out-scattering，到达A点。

$I_{p\rightarrow a} = I_{p} S(\lambda,\theta,h)T(\bar{pa})$

前面刚刚有推导，代入到上面公式中

$I_{p\rightarrow a} = I_{c}T(\bar{cp})S(\lambda,\theta,h)T(\bar{pa}) = I_{c}S(\lambda,\theta,h)T(\bar{cp})T(\bar{pa})$
数值积分(The Numerical Integration)

上面的公式推导是表示从C点入射的光线最后有多少到达A点，那么怎么表示A点收到的所有光线呢，这时候就需要  对于P点从A到B的积分。

$I_{a} = \sum_{a}^{b}{I_{p_{i}\rightarrow a}ds}$
将  的推导结果代入到上面公式，将  替换为统一的 $I_{s}$ ，表示统一的太阳光强度，前面有解释过因为地球相对太阳距离较远，可以将太阳光认为是directional light。

$I_{a} = \sum_{a}^{b}{I_{c}S(\lambda,\theta,h) T(\bar{c p_{i}} )T(\bar{ p_{i}a})ds} = I_{s}\sum_{a}^{b}{S(\lambda,\theta,h) T(\bar{c p_{i}} )T(\bar{ p_{i}a})ds}$

公式推导到这里，我们剩下的问题就是怎么计算散射函数和透射函数
2.1 计算散射函数

瑞利散射（Rayleigh scattering）

如果在网上查为什么天空是蓝色的，一定会查到Rayleigh scattering。那这个到底是什么意思呢，Rayleigh scattering表示的是光线碰到半径比光的波长小很多的微小颗粒（例如单个原子或分子）之后的散射行为，比如说氧气分子，氮气分子的大小都比光的波长小很多，光线碰到这些大气中占大部分的分子，就会发生Rayleigh scattering。
另外还有还有一个散射行为叫Mie scattering，这个描述的是光线碰到跟光波长差不多大小的颗粒之后的散射行为，比如说花粉，粉尘，污染物。Mie scattering会导致云显示白色，如果空气污染严重，天空中白茫茫的，这个也是Mie scattering的结果。
这里我们主要想要模拟出天空蓝色的效果，所以只考虑Rayleigh scattering。
Rayleigh scattering的散射函数如下：

$S(\lambda,\theta,h) = \frac{\pi^{2}(n^{2}-1)^2}{2}\cdot\frac{\rho(h)}{N}\cdot\frac{1}{\lambda^{4}}\cdot(1+cos^{2}\theta)$
上面公式中，
n表示空气的折射系数，是个常量1.00029。
h表示高度，  表示密度比率，这个值在海平面处为1，随着高度增加，指数降低。
N也是常量，为标准大气的分子数密度，值为 $2.504\times 10^{25}$ 。

表示入射光线的波长。

表示散射的角度。
当高度h为0时，红光蓝光绿光的散射函数的极坐标图如下：

可以看出，同一个波长的光不同角度的散射是不同的，另外不同波长的光的散射也是不同的，图中可以看出蓝光的散射比红光多的多，导致蓝色的天空。
瑞利散射系数(Rayleigh Scattering Coefficient)

上面的瑞利散射函数只表示有多少光被散射到某一个方向，并没有表示所有方向上到底有多少光被散射。对瑞利散射函数做球面积分，可以计算出光在所有方向上的散射能量。结果如下:

$\beta(\lambda,h) = \int_{0}^{2\pi}\int_{0}^{\pi}S(\lambda,\theta,h)sin\theta d\theta d\phi=\frac{8\pi^{3}(n^{2}-1)^{2}}{3}\cdot\frac{\rho(h)}{N}\cdot\frac{1}{\lambda^{4}}$
这个函数表示光线在碰到一个氧气或者氮气分子发生散射损失的能量，这个常被称为瑞利散射系数(Rayleigh Scattering Coefficient)，这个也可以表示衰减系数，用于计算透射函数
要计算  函数是非常费时的，不过我们可以把这个函数拆成两部分

$\beta(\lambda,h) = \frac{8\pi^{3}(n^{2}-1)^{2}}{3}\cdot\frac{1}{N}\cdot\frac{1}{\lambda^{4}}\cdot\rho(h)$
其中前一部分可以由  表示

$\beta(\lambda) = \frac{8\pi^{3}(n^{2}-1)^{2}}{3}\cdot\frac{1}{N}\cdot\frac{1}{\lambda^{4}}$
这个函数对于某一个波长的光，值是一个常量，可以提前计算，因为高度h为0时， $\rho(h)=1$ ，也被称为海平面高度的瑞利散射系数。下图为不同波长的的值。可以看出瑞利散射系数跟波长的关系，地球大气层对蓝光的散射远远大于绿光和红光。

瑞利相位函数(Rayleigh Phase Function)

瑞利散射函数可以由两部分组成，第一部分为散射系数，控制散射的强度，第二部分为相位函数  ，控制散射的方向。

$S(\lambda,\theta,h) = \beta(\lambda,h)r(\theta)$
因为刚刚已经计算出  ，可以得到

$r(\theta) = \frac{S(\lambda,\theta,h)}{\beta(\lambda,h)} = \frac{\pi^{2}(n^{2}-1)^2}{2}\cdot\frac{\rho(h)}{N}\cdot\frac{1}{\lambda^{4}}\cdot(1+cos^{2}\theta) \cdot \frac{3}{8\pi^{3}(n^{2}-1)^{2}}\cdot\frac{N}{\rho(h)}\cdot\lambda^{4} = \frac{3}{16\pi}(1+cos^{2}\theta)$
从相位函数可以看出，它的值只与  有关，而与光的波长无关。
在之后的计算中我们可以理解为什么将瑞利散射函数拆成两部分。
2.2 计算透射函数

大气密度比率(Atmospheric Density Ratio)

前面公式推导中实际上有涉及到大气密度比率  ，这个是由某个高度大气的密度与海平面大气的密度相除得到的。

$\rho(h) = \frac{density(h)}{density(0)}$
因此可知海平面处的 $\rho(0) = 1$
大气密度越高，意味着光线越有可能跟大气中的分子发生散射。瑞利散射主要发生在海拔60Km以下的对流层，在这个高度范围，大气的密度曲线如下图。

那么大气密度比率的曲线如下（其中的Real density ratio）

因为真实曲线的特点是指数衰减的，我们可以用一个指数函数来拟合。

$\rho(h) = e^{-\frac{h}{H_{0}}}$
上面公式中h还是表示高度， $H_{0}$ 为scale height(均质大气高度)，对于瑞利散射，这个值取8500米，对于米氏散射，这个值取1200米。上图中黑色曲线就是这个拟合的指数函数的曲线。

指数衰减(Exponential Decay)

上文中我们实际上已经推导出瑞利散射系数 $\beta(\lambda,h)$ ，这个系数实际上表示光线发生散射损失的能量，那么初始强度为 $I_{0}$ 的光线在大气中传播一段距离后，剩下的光线强度 $I_{1}$ 可以由下面公式表示

$I_{1} = I_{0} - I_{0}\beta = I_{0}(1 - \beta)$
那么经过一段传播距离后，余下的能量为：

简单情况：如果β为定值，我们可以直接通过微积分得到衰减一段距离x后的剩余能量：

$I_{p} = I_{c} e^{-\beta x}$

复杂情况：β是与波长λ和高度h都相关，即不能简单地当成定值。此时，我们就需要用积分来表示了：

$I _{p} = I_{c} e^{-\sum_{c}^{p}\beta(\lambda,h_{q})ds}$

我们在透射函数部分期望可以求出  ，用来求出

我们现在可以知道

$T(\bar{cp}) = e^{-\sum_{c}^{p}\beta(\lambda,h_{q})ds}$
对于 $\beta(\lambda,h_{q})$ ，前面有说过它可以分成两部分即

$\beta(\lambda,h) = \frac{8\pi^{3}(n^{2}-1)^{2}}{3}\cdot\frac{1}{N}\cdot\frac{1}{\lambda^{4}}\cdot\rho(h) = \beta(\lambda) \cdot\rho(h)$
上面  是常量，那么

$T(\bar{cp}) = e^{-\sum_{c}^{p}\beta(\lambda,h_{q})ds} = e^{-\sum_{c}^{p}\beta(\lambda)\rho(h_{q})ds} = e^{-\beta(\lambda)\sum_{c}^{p}\rho(h_{q})ds}$
上面公式中的 $\sum_{c}^{p}\rho(h_{q})ds$ 被称为光学深度(optical depth)，由 $D(\bar{cp})$ 表示。

$T(\bar{cp}) = e^{-\beta(\lambda)D(\bar{cp})}$
公式推导到这里，散射函数和透射函数都已经知道了，后面就可以对摄像机的View Ray进行采样，计算最后到达摄像机的光线。

2.3 采样View Ray(AB)

回顾一下，前面我们已经推导出

$I_{a} = I_{s}\sum_{a}^{b}{S(\lambda,\theta,h) T(\bar{c p_{i}} )T(\bar{ p_{i}a})ds}$

因为前面已经将  拆分

$S(\lambda,\theta,h) = \beta(\lambda,h)r(\theta) = \beta(\lambda)\rho(h)r(\theta)$
将  代入  得到

$I_{a} = I_{s}\sum_{a}^{b}{\beta(\lambda)\rho(h)r(\theta) T(\bar{c p_{i}} )T(\bar{ p_{i}a})ds}$
对于上图中AB之间要采样的P点，  和  是常量，所以可以得出

$I_{a} = I_{s}\beta(\lambda)r(\theta)\sum_{a}^{b}{\rho(h) T(\bar{c p_{i}} )T(\bar{ p_{i}a})ds}$
前面已经推导出透射函数，那么

$T(\bar{c p_{i}} )T(\bar{ p_{i}a}) = e^{-\beta(\lambda)D(\bar{cp_{i}})} \cdot e^{-\beta(\lambda)D(\bar{p_{i}a})} = e^{-\beta(\lambda)(D(\bar{cp_{i}})+D(\bar{p_{i}a}))}$
在前面的推导中我们知道光学深度也是需要采样才能算出来的

$D(\bar{cp}) = \sum_{c}^{p}\rho(h_{q})ds$
从A点往  采样计算出可以累加得到 $D(\bar{p_{i}a})$

从点往 $C_{i}$ 点开始采样，可以得到 $D(\bar{p_{i}c_{i}})$
这样我们最后能够算出

3. 总结

从上面的公式中可以看出这样计算  是非常费时的，有没有更节省时间的方法来计算  呢？
下一篇文章将会介绍Unity的<Skybox/Procedural>Shader 中所用的算法，是如何达到同样的效果的。

参考文章

这篇文章主要参考了Alan Zucconi的系列文章
也参考了冯乐乐的文章
PZZZB在URP中实现了同样的算法

		自动登录	找回密码
密码			立即注册

[笔记] 在Unity中制作Procedural Sky (Part 1-Atmospheric …

本帖子中包含更多资源