多重散射的天空光照效果建模与实时绘制

艾祖亮， 张立民， 于文龙

（海军航空工程学院，山东 烟台 264001）

多重散射的天空光照效果建模与实时绘制

艾祖亮， 张立民， 于文龙

（海军航空工程学院，山东 烟台 264001）

依据大气散射的物理原理，提出了一种考虑多重散射的天空光照效果建模与实时绘制方法。该方法首先以太阳和天空光作为光源建立了多重散射的天空光照效果模型，然后综合多种大气粒子密度, 采用合理的分段采样策略，对天空颜色模型的积分进行简化，以减少积分近似计算所带来的误差；通过对简化后的模型进行分析提出了采用二维纹理与三维纹理对光学深度预计算的方法，避免了运行时计算光学深度积分的问题；最后该算法在GPU的片段处理器上执行，实现了天空光照效果模型的实时绘制，可以满足各种实时应用需求。

多重散射；天空光照；光学深度；预计算；实时绘制

太阳光在大气中的多重散射对于天空光照效果的逼真模拟具有重要的意义。大多数大气散射模型[1-5]在模拟天空光照效果时只考虑单散射，而不考虑多重散射光。由于单散射模型忽略了天空体本身作为天空光光源的影响，并不适合模拟黎明或黄昏时的天空光照效果，这主要是因为此时太阳位于水平线下方，仅仅离太阳方位近的区域才能够收到直接的太阳光，而大气层其他区域只能收到非直接的多重散射的天空光。

因此，许多学者提出了多重散射模拟方法，Preetham等[6]通过一个分析模型去拟合多散射蒙特卡罗（Monte-Carlo）仿真的结果，但是他们的模型要求视点必须在地面上才有效；Nishita等[7]和Haber等[8]用体积辐射度方法计算多重散射，但是他们的方法很难达到实时性；Bruneton和Neyret[9]提出了一个实时大气散射绘制方法，并考虑了多重散射，该方法显著提高了多重散射模型绘制速度，但在实时性上距离工程实际应用仍然存在一定的距离。

国内，刘世光和彭群生[10]及钱徽和朱淼良[11]分别对大气效果的建模和绘制技术进行了综述；王长波等[12]提出一种考虑大气折射的天空光照模型, 综合考虑大气压强、温度及水汽压等因素,实现了日、月升落中形状及其周围光晕色彩的变化等效果，但是他们的方法非常复杂，不能达到实时性；黄雷等[13]提出一种适用于雨雾天气的天空光照的多粒子散射模型，实现了雨雾场景中大气散射效果及雨后彩虹的真实感绘制，但是该模型仅仅只考虑了单散射；杜芳等[14]和冯玉康等[15]分别基于GPU对Nishita等[16]的大气单散射模型进行了实现，但对于多重散射并没涉及。

本文提出了一种考虑多重散射的天空光照效果建模与实时绘制方法。该方法首先以太阳和天空光作为光源建立了多重散射的天空光照效果模型，然后根据多种大气粒子密度采用合理的分段采样策略对天空颜色模型的积分进行简化；通过对简化后的模型进行分析提出了采用二维纹理与三维纹理对光学深度预计算的方法，避免了运行时计算光学深度积分的问题；最后该算法在 GPU的片断处理器上执行，实现了多重散射的天空光照效果的实时绘制，可以满足各种实时应用需求。

1 天空光照效果模型

1.1 太阳光与天空光模型

多重散射模型以太阳光和天空光为光源，通过借鉴Nishita模型的思想[7]对天空光进行建模，以大气中任意一点p为中心，按照方位角 φj和高度角ψj采样把整个天空体分割成若干个发光带，其中，单位为度。发光带的强度由采样线的强度确定，为了简化计算，模型不考虑地面方向的多重散射光，如图1所示。

图1 太阳光与天空光

入射到大气中任意一点p的光强是太阳光和天空体各发光带入射到该点处光强的总和Ip(λ)，公式如下：

其中，λ为光线波长，天空体发光带数目为 Nsky，θj为天空光发光带采样线入射方向与太阳光入射方向的夹角，也称散射角。

首先计算太阳光 Isun,p(λ)，s为太阳光经过点p与天空体圆顶的交点，Is(λ)为太阳光在进入天空体圆顶前的初始入射光强，根据大气光衰减模型[17]可知：

式(2)中， t(s p,λ) 为光学深度，定义为大气粒子的消光系数 βex（λ） 与密度比ρ的乘积沿着光线传输路径的积分[16]，公式如下：

其中， H0为大气粒子均值高度，不同类型大气粒子的消光系数与密度比是不同的，对于空气分子主要发生Rayleigh散射，消光系数 β(λ)等于散射系数λ)；对于其他大气悬浮粒子主要发生Mie散射，消光系数(λ)等于粒子的散射系数β(λ) 与吸收系数 β之和。公式(3)中的βex(λ)ρ实质上是等于大气中所有存在的粒子消光系数与密度比乘积的和。

然后计算天空光发光带 Isky,p(λ,φj,ψj,θj)， pj为发光带采样线 ajp上任意一点， sj为经过 pj的太阳光与天空体圆顶的交点，则点 pj处的光强为：

其中， F(θj,λ)为散射相位函数[18]，描述了入射光线在某一点散射后进入特定散射方向的光能占整体入射光能的比例，它取决于发生散射时的光线波长、粒子尺寸和散射角。

光强 Ipj(λ,θj)经过衰减到达p处的光强为：

天空光发光带Isky,p(λ,φj, ψj,θj)为 光 强Ipjp(λ,θj)沿着采样线 ajp的积分：

1.2 天空颜色模型

天空光照计算是天空颜色（sky color）的计算，主要针对天空体圆顶（sky dome）的光照效果模型，如图2所示。

图2 天空颜色模型

图2中， pv为视点， pa为视线与天空体圆顶的交点，θ为视线方向与太阳光入射方向的夹角，θsky,j为视线方向与天空光发光带采样线入射方向的夹角，θj为天空光发光带采样线入射方向与太阳光入射方向的夹角。点p处的太阳光与天空光总的入射光强 Ip（λ）可以通过上一节公式得出，则光强 Ip（λ）在点p处发生散射进入视线方向的光强为：

视线方向总的入射光强为 Ippv(λ,θ,θsky,j)沿着视线 pvpa的积分：

上式即为天空颜色模型，通过该公式可以计算出大气中任意一点所有视线方向的天空颜色。

2 天空光照效果的实时绘制

2.1 模型的积分简化

天空颜色模型 Iv(λ,θ,θsky,j)的计算涉及太阳光和天空光沿着视线方向的积分，天空光发光带Isky,p(λ,φj,ψj,θj)也需要计算太阳光沿着发光带采样线方向的积分，而上述积分并无解析表达且计算非常耗时，因此需对该式进行简化。

积分的近似计算通常采用对积分路径分段采样，然后通过矩形、梯形或者抛物线法来求得积分结果，积分路径的分段间隔往往是相同的。但是对于大气散射模型来说，如果采用常量间隔的分段采样方式，积分近似的误差会非常大，这主要是因为大气粒子密度是随着海拔高度以指数方式递减的，理想的情况应该是分段间隔与大气粒子密度成反比，小的间隔对应于低的高度，大的间隔对应于高的高度。

为此，按照海拔高度把大气层离散为具有一定厚度的一层层球壳和 Hi,max分别为每层球壳的内表面与外表面高度，则有）。通过这些球壳对视线方向上的积分路径进行分段采样，每一段内的粒子密度采样值取该段内中间点的密度值，每一层球壳的厚度就等于分段间隔，如图3所示。

图3 大气层球壳

大气层的厚度大约在1000km以上，并没有明显的界限，但在飞行模拟视景仿真中，通常只考虑离地球海拔高度为35km以下的大气层，即假设大气层最大的厚度 Hmax= HN,max= 35km 。之所以这样假设是因为在此高度之上的大气粒子已经非常稀薄，对大气散射光照效果模型计算结果的影响完全可以忽略。

积分路径的分段大小要综合地考虑大气粒子密度的分布情况，使其在段中近似恒定，分段数N可根据绘制的精度来调整。为了得到合理的分段间隔，假设每层球壳中的粒子密度衰减都相同[8]：

由于不同粒子的均值高度 H0不同，空气粒子均值高度通常为7.994km，灰尘粒子的均值高度为1.2km，其他雾粒子、雪粒子、雨粒子的均值高度与具体气象条件有关，为了使得分段间隔更加合理，公式(10)实际采用的均值高度为上述所有存在粒子均值高度的平均值 Ha。

求解公式(10)得到如下：

其中， H1,min= 0, Hmax= 35km , i= 1,… ,N -1。

通过该递推公式可以计算出所有球壳内外表面高度以及球壳厚度，然后用这些球壳对天空颜色模型和天空光模型的积分路径分段采样，每一段内的采样点取该段的中心点，显然每段中心点的高度即是相应球壳的中心高度 (Hi,min+ Hi,max)/2，每一段内的粒子密度就是该中心高度的密度。通过上述分段采样方法，公式(7)与公式(9)可以简化为：

其中，Nsky,j表示天空光发光带采样线方向积分路径分段采样数目， Nv表示天空颜色模型视线方向积分路径分段采样数目， pj,k与 pi分别为相应的采样点。

2.2 模型的预计算与绘制

经过上节的简化解决了天空颜色模型与天空光发光带的视线方向积分计算问题，但是光学深度积分的计算仍然不可避免，O'Neil[19]针对光学深度的计算提出了预计算的方法，用一个二维纹理存储光学深度的计算结果，这张表将光学深度的积分式替换成一张二维纹理查找表，表的两个维度分别是视点高度与视线方向的天顶角。二维纹理查找表方法利用太阳光是平行光这一近似假设，通过查找表就可以获得大气中任何一点到大气层顶部的光学深度，避免了运行时计算光学深度的问题。采用二维纹理查找表对于单散射天空光照效果模型的计算取得了实时绘制效果，但是对于多重散射天空光照效果模型，由于光学深度的计算量巨大，只采用二维纹理查找表仍然难以达到实时绘制要求。因此，针对这个问题，本文提出了一种采用二维纹理和三维纹理对光学深度预计算的方法，并基于 GPU实现了多重散射天空光照效果模型的实时绘制。

下面对天空光照颜色模型公式(13)的计算进行分析，公式第一项任意采样点 pi来自太阳光的光强贡献为：

公式第二项任意采样点 pi来自天空光某一方向发光带的光强贡献为：

t(papv,λ)光学深度的值与视点的高度、视线方向的天顶角有关，可以预先计算为以视点的高度和视线方向的天顶角为参数的一个二维纹理查找表。

F(θ,λ)、 F(θj,λ)散射相位函数的值与散射角有关，可以预先计算为以散射角为参数的一个一维纹理查找表。

ρpi、 ρpj,k大气粒子密度比的值与高度相关，可以预先计算以高度为参数的一个一维纹理查找表。

由于大气中不同种粒子的光学深度 t(S,λ)、散射相位函数 F(θ,λ)、密度比ρ都是不同的，因此每个预计算的纹理查找表都需要存储多个数值。本文的模型主要考虑了5种类型的大气粒子，它们是空气粒子、灰尘粒子、雾粒子、雨粒子、雪粒子，其中空气分子发生Rayleigh散射，后面4种粒子发生Mie散射。每个纹理查找表的值都包括4个通道，可以存储4个float值，前3个通道分别存储与空气分子、灰尘粒子、雾粒子相关的预计算结果，第4个通道存储与雨粒子或者雪粒子相关的预计算结果，即雨粒子、雪粒子同时只能有一个有效。

通过上述4个预计算的纹理查找表，天空颜色模型表达式的每一项值都可以在运行中通过查表获得，然后只需要通过简单加减乘积运算就可以计算出天空颜色模型的值，因此该方法非常适合在GPU上执行。算法的绘制在OpenGL片段处理器执行，通过简单地调用 texure1D、texure2D和texure3D等几个着色器指令函数就可以直接从已经预计算好的4个纹理查找表中获取光学深度、散射相位函数、密度比的值，从而实现天空颜色模型的实时计算与绘制。

3 实验及结果分析

本文算法的测试环境如下：硬件平台为Intel酷睿i7 3770 CPU、4 GB内存、NVIDIA GeForce GTX 660显卡（2GB显存），软件平台为Windows 7、OpenGL 4.1、GLEW 1.9。

以下所有实验只考虑空气分子和灰尘两种大气粒子，不考虑雾粒子、雨粒子、雪粒子。天空光发光带数目 Nsky为8，大气球壳分层数目N为20。空气分子发生Rayleigh散射，均值高度 HR为7.994km，对应于光线3个波长(680nm, 550nm, 440nm)的散射系数为(0.0058, 0.0135, 0.0331)；灰尘粒子发生Mie散射。

下面对3种不同方法在实际应用中的绘制效率进行实验分析。表1为相同场景环境下采用3种不同方法所消耗的绘制时间对比，生成图像分辨率为1024×768。为了使仿真环境更加真实，场景地形数据库为台湾地区，地形图像分辨率为5m，数据格式为TerraPage(.txp)，对txp格式的多线程页面调度支持采用 OSG开源图形引擎代码。从绘制各阶段消耗时间以及帧数结果可以看出，本文的多重散射方法相比Bruneton和Neyret[9]的方法在绘制速度上具有明显的优势，基于颜色插值的方法[4]虽然绘制速度最快，但由于只是一个简单经验方法，没有任何基于大气物理模型的仿真计算，不能模拟出真实的天空光照效果。

下面针对太阳以及灰尘粒子不同的高度绘制相应的天空光照效果。

图4为太阳不同高度下的天空光照效果图，灰尘粒子均值高度 H 为1.2km，散射系数为

M0.02，消光系数为 0.0222，从绘制的效果可以看出本文的算法能够实时模拟不同太阳高度的真实天空光照效果。

图5为灰尘粒子不同均值高度下的天空光照效果图，灰尘粒子均值高度从左至右依次为0.6、1.2、3.0、4.0、5.0、6.0，单位为km，散射系数为0.02，消光系数为0.0222。从绘制的效果可以看出随着粒子均值高度的增加，粒子的密度越来越大，相应对太阳光的吸收和散射效果也越明显。

下面针对灰尘粒子不同的散射系数和消光系数绘制相应的天空光照效果。

图6所示为灰尘粒子不同散射系数下的天空光照效果，灰尘粒子均值高度为1.2km，散射系数从左至右依次为0.005、0.010、0.020、0.030、0.040、0.050。从绘制的效果可以看出随着散射系数的增加，太阳光在空气中散射的效果越来越明显。

图7所示为灰尘粒子不同消光系数下的天空光照效果，灰尘粒子均值高度为1.2km，散射系数为 0.02，消光系数从左至右依次为0.4000、0.2500、0.2000、0.0667、0.0333、0.0222。从绘制的效果可以看出随着消光系数的减少，光线在大气中的衰减越来越小，相应天空越来越明亮，颜色也越来越蓝。

表1 不同方法绘制时间的对比

图4 太阳不同高度下的天空光照效果

图5 灰尘粒子不同均值高度下的天空光照效果

图6 灰尘粒子不同散射系数下的天空光照效果

图7 灰尘粒子不同消光系数下的天空光照效果

4 结 论

本文提出了一种多重散射的天空光照效果建模与实时绘制方法，实验结果证明该方法能够实时模拟不同时间、不同大气条件下的真实感天空光照效果，可以满足各种实时应用需求。

[1] Dobashi Y, Nishita T, Kaneda K, Yamashita H. A fast display method of sky color using basis functions [J]. The Journal of Visualization and Computer Animation, 1997, 8(2): 115-127.

[2] Dobashi Y, Yamamoto T, Nishita T. Interactive rendering of atmospheric scattering effects using graphics hardware [C]//HWWS 2002 Proceedings of the ACM SIGGRAPH/EUROGRAPHICS Conference on Graphics Hardware, 2002: 99-107.

[3] O'Neil S. Accurate atmospheric scattering [M]. GPU Gems 2: Programming Techniques for High-Performance Graphics and General-Purpose Computation. Addison-Wesley, 2005: 253-268.

[4] Abad J A. A fast, simple method to render sky color using gradients maps[EB/OL]. (2006-10-09) [2013-04-03]. http://www.geocities.com/ngdash/whitepapers/skydom ecolor.html.

[5] Schafhitzel T, Falk M, Ertl T. Real-time rendering of planets with atmospheres [C]//WSCG International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision, 2007: 91-98.

[6] Preetham A J, Shirley P, Smits B. A practical analytic model for daylight [C]//SIGGRAPH 99, 1999: 91-100.

[7] Nishita T, Dobashi Y, Kaneda K. Display method of the sky color taking into account multiple scattering [C]//Pacific Graphics 1996, 1996: 117-132.

[8] Haber J, Magnor M, Seidel H P. Physically based simulation of twilight phenomena [J]. ACM Transactions on Graphics, 2005, 24(4): 1353-1373.

[9] Bruneton E, Neyret F. Precomputed atmospheric scattering[J]. Computer Graphics Forum, 2008, 27(4): 1079-1086.

[10] 刘世光, 彭群生. 气象景观的真实感模拟技术综述[J]. 计算机辅助设计与图形学学报, 2008, 20(4): 409-416.

[11] 钱 徽, 朱淼良. 大气效果的绘制技术综述[J]. 计算机辅助设计与图形学学报, 2002, 14(5): 475-482.

[12] 王长波, 王章野, 曾 运, 彭群生. 考虑大气折射的天空场景真实感绘制[J]. 计算机学报, 2005, 28(6): 939-949.

[13] 黄 雷, 王章野, 王长波. 雨雾天气下光线散射效果的实时绘制[J]. 软件学报, 2006, 17(11): 126-137.

[14] 杜 芳, 张 炎, 晁建刚, 王金坤. 基于GPU的地球大气散射现象可视化仿真[J]. 系统仿真学报, 2009, 21(2): 147-150.

[15] 冯玉康, 周圣川, 马纯永, 韩 勇, 陈 戈. 基于GPU的地球大气层和三维体积云仿真[J]. 计算机工程, 2012, 38(19): 218-225.

[16] Nishita T, Sirai T, Tadamura K. Display of the earth taking into account atmospheric scattering [C]// SIGGRAPH 93, 1993: 175-182.

[17] Hoffman N, Preetham A J. Rendering outdoor light scattering in real time [C]//Game Developer Conference, 2002: 337-352.

[18] Cornette W M, Shanks J G. Physical reasonable analytic expression for the single scattering phase function [J]. Applied Optics, 1992, 31(16): 3152-3160.

[19] O'Neil S. Real-time atmospheric scattering[EB/OL]. (2004-05-03)[2013-04-03]. http://www.gamedev.net/ page/resources/_/technical/graphics-programming-an d-theory/real-time-atmospatmos-scattering-r2093.

Modeling and Real-time Rendering of Sky Illumination Effects Taking Account of Multiple Scattering

Ai Zuliang, Zhang Limin, Yu Wenlong
(Naval Aeronautical Engineering Institute, Yantai Shandong 264001, China)

This paper presents a method of modeling and real-time rendering of the sky illumination effects taking account of multiple scattering with atmospheric scattering based on physical models. Firstly, with the sun and the sky as light sources, the model of the sky illumination effects taking account of multiple scattering is built. Then a reasonable segmented sampling strategy based on a variety of atmospheric particle densities is adopted to simplify the integral of the sky color model to reduce the integral approximation calculation error. By analyzing the simplified model, the pre-computed 2D texture and 3D texture for optical depth are proposed to avoid a run-time calculation of optical depth. Finally, the algorithm is executed on the fragment processor of GPU so that the sky illumination effects taking account of multiple scattering are rendered with high frame rate for real-time applications.

multiple scattering; sky illumination; optical depth; pre-compute; real-time rendering

TP 391.9

A

2095-302X (2014)02-0181-07

2013-04-03；定稿日期：2013-06-28

艾祖亮（1980-），男，湖北随州人，博士研究生。主要研究方向为自然现象模拟、实时绘制技术。E-mail：aizuliang@live.cn