秦绪佳，赵 铮，柯玲玲，郑红波，马 骥

(浙江工业大学 计算机科学与技术学院, 杭州 310032)

E-mail:zhb@zjut.edu.cn

1 引 言

云可分为三大云族，十大云属，是一种自然现象.它们是由聚集在空气中的微尘周围的水分子产生的水滴或冰晶散射产生的.云没有特定的表面，也没有清晰的边界，且晶体结构时时刻刻都在变化，表面呈现出的透明度也不断在改变.每一朵云都是不尽相同的，所以对云的可视化也是十分复杂的.现在存在的一些有不错视觉效果的云仿真方法大多依赖于数学模型，而这些模型都十分复杂且计算成本相对较高.所以对三维云景的绘制要达到高效率和高仿真效果是当前计算机图形学领域面临的一大挑战.

对三维云的可视化在很多日常仿真应用中都是不可或缺的.在很多游戏中，尤其是在大型网游中，云仿真技术的水平对游戏场景的真实感影响极大，真实感较强的云仿真可以使场景更逼真，用户的体验感也会更好；在飞行视景仿真系统中，云层效果的真实以及良好的实时性也是至关重要的；在气象研究中，真实的气象数据可以反映真实的天气情况，对气象的研究、气象的预报等都有极大的意义，云作为气象数据中不可或缺的一部分，对真实气象数据的深入研究和利用都十分重要.

想要得到真实性较高的云的颜色就需要模拟其内部的光照效果.胡香等人[14]提出了一种利用云的多向前方散射特性建立的体云光照模型，效果逼真，与真云一样具有银色的边缘.Xiaoyan Hu等人[15]提出了一种简单的实时云建模与绘制方法，通过模拟云粒子在二维平面上的投影运动而形成云的密度，然后使用一种实时近似散射的算法来渲染.但这类模型忽略光的多向散射，只计算视点方向的外散射，只适用于描述低散射率的光照效果但不能处理云团等高散射率媒介的光散射问题.Zhengwei Sui等人[16]认为云是一个收集来自太阳光的粒子集合，而三维云是用八叉树索引和粒子系统建模的，故使用散射光模型和广告牌技术来渲染它.

开始先使用传统的Ray Casting算法来绘制三维云，为了提高绘制和渲染的速度，还使用了光栅化求交、提前光线终止法和GPU加速等方式.因为所使用的网格数据精度不够高、层级较少使得模拟数据在垂直方向上较为稀少，导致云表面十分粗糙、渲染缺乏真实性.为了改善这些问题，本文在传统体绘制Ray Casting算法的基础上，生成3D Perlin噪声纹理并采样，再和气象数据的原始纹理以不同的比例进行融合.用这样的方式对云表面产生随机扰动，增加所绘制云的真实感.最后，因为涉及到材质的各向异性，本文还使用了一种结合Cook Torrance光照模型和双向反射分布函数(BRDF)的新的光照模型，以此来模拟光的散射，从而进一步增强仿真效果.

2 基于GPU的Ray Casting算法绘制三维云

由于现在常用的光线投射算法三线性插值运算，需要基于GPU,会产生一定时延.所以本文利用GPU高效的浮点算术功能来加快渲染.首先把体数据载入到显存之中，再用一个与其等大的三维纹理表示，其中和这个三位纹理相对应的集合对象被称作代理几何.这样可以很好的解决GPU无法直接访问内存数据的问题.本章中Ray Casting算法的代理几何为立方体，需要将其先转换映射到三维球体上，最后和光线相交计算得到采样点位置.

2.1 气象体数据源及数据格式

1)体数据源

气象数据有结构化、规则的网格数据，也有非结构化的散乱数据.前者是通过天气模型得到的；后者是利用气象观测得到的真实数据.气象观测指的是研究测量和观测地球大气中温度、湿度、大气气压、降水、云量等气象.非结构化的气象观测数据可通过插值、重采样转化成结构化数据.

2)数据格式

本文利用通过气象观测得到的真实数据作为研究数据源.该数据源涉及到相对湿度、降水量、雷达数据等.其中相对湿度使用的数据集是在不同的压强下每3小时进行一次采样获得，以0时为起点，总共获得80个样本.一共有361×181个格点数据，每个格点占据4个字节，存储的是单精度的浮点数.

图1 相对湿度预报格点数据实例Fig.1 An example of grid data for relative humidity prediction

图1为本实验所用到的一个相对湿度预报格点数据实例.前4行为数据头，存储的是数据类型、描述信息、数据采样的相关信息以及数据经纬度的相关信息等.从第5行开始存储数据部分，数据按先纬度向后经度的顺序存放.图1中所示数据是106类数据，该数据的描述信息为“T799模式2017年06月04日20时00分200hPa相对湿度030小时预报(0602)”.该数据是从2017年6月4日20:00开始采样的，并且在200百帕压强的条件一共采样了30小时.数据采样从经度0°，纬度360°开始，到经度-90°，纬度90°，每隔1经度、1纬度采样.将绘制球面分成361×181个格点，从最小等值线值-20到最大等值线值110每隔10个单位长度绘制.之后的数据即每个格点内的相对湿度值.

本实验中需要用到三维体数据，而使用的数据源是多层二维网格数据.故将这些数据垂直逐层叠加形成三维空间的体数据.再将格点数据映射到灰度图上，经过归一化处理后得到的数据将用于下文研究.

2.2 Ray Casting算法实现过程

基于GPU的Ray Casting可以分为5个部分，绘制、归一化、计算、采样最后绘制到网格上.其具体描述如下：

1)绘制六面体包围盒，将所使用的气象体数据源中的网格坐标转化为三维球体上的经纬度.

2)将二维网格数据垂直叠加得到的三维体数据进行归一化处理后做成三维纹理，传到着色器shader中.

3)确认投射光线传播的方向以及光线在该方向上前进的最远传播距离.

首先要找到光线穿过体纹理时和距离视点最近以及最远的两个面的交点，这两个交点分别为光线的出点和入点.如何计算得到出入点是影响整个算法性能的一个关键点，光栅化以及光线/包围盒交叉测试这两种方法都可以求得交点，相比较而言前者的效率会更高，故本文采用这种方式来计算出入点.

4)沿着射线前进方向进行等距离的采样，并累计计算该方向的不透明度以及颜色值.

为了使绘制速度更快，本文参考了Levoy等提出的提前光线终止法[17]，当累计的不透明度达到阈值0.99时就终止计算.

5)最后生成一个二维网格面片，将该面片转换为球坐标后包裹球体，并在该网格面上绘制前面累计的像素值.

3 结合纹理抖动和BRDF的云景增强真实感绘制

3.1 生成3D Perlin噪声纹理

在三维可视化中要实现降噪有很多方法，周天琪等人[18]提出了一种将三维数据分为多个相邻二维切片数据分成降噪的算法.本文使用的是生成3D Perlin噪声纹理，将任意的一点p(x,y,z)输入到已知的三维网格数据空间中，并为附近的网格点都设置伪随机梯度向量.Perlin首先提出了将256个伪随机梯度向量预先存储在一个二维数组G[256][3]之中，然后将0～255这256个数字随机放入一个一维排列数组P[256]中，再用式子G[P[P[x]+y]+z]来得到最终的随机梯度向量值.后来Perlin又提出了两点改进：

1)构造新的缓和曲线函数f(x)=6x5-15x4+10x3代替先前的缓和曲线函数f(x)=3x2-2x3可以得到更加的平滑的噪声值.

2)如图2所示，以p为基准点，围绕p最近的8个网格点形成一个单位立方体，点p指向立方体12条边中点构成12个向量.用这12个向量代替先前方法中的随机梯度向量，最终的随机梯度向量由公式G[p[p[p[x]+y]+z]]计算得到.采用这种方式得到的随机梯度向量虽然不是完全随机的，但能消除梯度中的方向偏差.

图2 三维方格Fig.2 Three dimensional grid

本文在生成伪随机梯度向量时，加入了以上2点后再生成3D Perlin噪声.

在求点积和加权平均值的过程中，首先将x-x0，y-y0，z-z0分别代入缓和函数f(x)中，得到比例系数：

(1)

再以p(x,y,z)为起点，以相邻8个网格端点为终点，分别生成向量：

(2)

以点p为起点，12条边中点为终点，生成12个向量：

(1,1,0),(-1,1,0),(1,-1,0),(-1,-1,0),
(1,0,1),(-1,0,1),(1,0,-1),(-1,0,1),
(0,1,1),(0,-1,1),(0,1,-1),(0,-1,-1)

将公式(2)中的每一个向量作为输入，并与上述伪随机向量中的某个进行点积运算，得到：

(3)

再将这些随机生成的离散点进行插值运算，以得到连续的噪声波.经典Perlin噪声使用的是线性插值，但本文使用的是余弦插值.因为余弦插值函数可以得到更加平滑的噪声，并且计算相对简单，可以很好的替代线性插值函数.

利用余弦插值计算出图2中点A1、A2、A3、A4的值为：

(4)

其中i=1-cos(sy×PI)×0.5.

利用余弦插值计算出图2中点B1、B2的值为：

(5)

其中i=1-cos(sz×PI)×0.5.

利用余弦插值计算出图2中点C的值即为最终结果，其值为：

Noise(x,y,z)=z0×(1-i)+z1i

(6)

其中i=1-cos(sz×PI)×0.5.注意公式(4)-公式(6)中的PI值均约等于3.1415927.

3.2 利用纹理噪声引入随机抖动

将光线投射时在体数据纹理上采样得到的基础值记作Valuesrc，其中体数据纹理是由气象数据生成，记作texVolume_.在这个基础值Valuesrc上加一个随机抖动，首先要在3D Perlin噪声生成的噪声纹理上采样，采样得到的噪声纹理记作texNoise_.再分别获取x、y、z方向上的微量偏移，记作Valuejitter.将Valuejitter和Valuesrc线性叠加，得到采样值Value.如公式(7)所示，其中rf为控制噪声对采样结果影响大小的调节因子：

Value=Valuesrc×(1-rf)+Valuejitter×rf

(7)

最后将Value带入映射函数得到颜色以及透明度，并沿着射线前进方向累加，最后得到合成图像.

3.3 Cook Torrance BRDF 光照模型

常见的简单着色模型都是将需要着色的平面认定为光滑的，这些平面只需要一个法向量就可以表示其方向；但微平面模型是将需要着色的平面认定为由无数个微小光滑平面组成的，这些微小平面远小于入射光线所涉及的区域.显然，微平面模型不可能只用一个法向量来表示.微平面模型需要用一个概率分布函数来描述，这个概率分布函数可以计算微平面中某个方向上着色区域存在的可能性大小.微平面的斜率变化大小可以很好的表现出待着色平面的粗糙程度.Cook-Torrance模型就可以很好的体现微平面的特性.

本文采用Cook-Torrance BRDF模型来近似模拟物体受到光照时由物其自身材质决定的散射的辐射量以及镜面反射的多少.该BRDF函数表示见公式(8)：

fr=kdflambert+ksfcook-torrance

(8)

漫反射部分由flambert表示，镜面反射部分由fcook-torrance表示.用kd来衡量给定物体自身材质表现出扩散行为，值越高其扩散行为越强烈；用ks来衡量物体表现出的行为是否接近镜子，值越高越接近.

由于现实中存在的材质接受能量一定不少于其反射的能量，故所有BRDF的权重和必须小于等于1，见公式(9)：

(9)

并且为了不过多消耗能量，需要始终保证kd+ks≤1.

Ward BRDF经验模型由漫反射和镜面反射组成.其中漫反射用Lambert项表示，镜面反射用椭圆高斯函数表示[19]，见公式(10)所示：

(10)

公式(10)中，单位矢量i表示入射的方向，o表示观察的方向，h表示入射方向和观察方向之间的半角矢量.θi表示相对于法线n的入射角度，θo表示相对于法线n的观察角度，θh表示相对于发现n的半角矢量角度.φh为投射到表面微平面上的半角矢量的方位角，ρd为漫反射率，ρs为镜面反射率.x轴和y轴在表面的切平面中，各向异性由椭圆高斯分布函数表示，并由σx和σy这两个粗糙参数的垂直斜率分布定义该函数.其中表面粗糙大于沿其垂直方向的方向为x轴的方向，即σx≥σy.

本文学习Ward BRDF的经验模型，借鉴其漫反射部分，对fcook-torrance定义见公式(11)：

(11)

式中D为微平面分布函数，该分布函数表示在指定点处微小平面的统计方向；虽然Phong函数和Beckmann函数都可以表示微平面的分布，但相对而言，GGX分布[20]有更强的尾部，从而产生更多的阴影，故本文借鉴Walter等人文章中提出的GGX分布[20]，其定义见公式(12)：

(12)

(13)

其中r表示曲面的粗糙度，范围在0-1之间(包括0和1)；m是要检查的方向矢量.

G表示由于微平面彼此遮挡所造成的光衰减.它近似模拟了绘制过程中在指定的微平面上与邻近几个微平面之间彼此遮掩光线和相互反弹光线的影响，由此计算在到达人眼之前这一系列过程中能量消耗的概率.

在微平面上反射光有以下3种现象，分别为入射光没有受到遮掩，100%的光照都到达观测人员眼睛；入射光被遮掩了一部分；以及反射光被遮掩了一部分.

几何函数G与GGX分布的D配对，见公式(14)所示：

(14)

随着r的值越来越大，即粗糙程度的增大，光线被遮挡的程度也不断增加.

公式(11)中F表示菲涅耳函数，该函数模拟了光在不同角度下和平面相互作用的方式.Shlick C[21]提出了一个可以快速评估得到结果的方式，公式(15)即为其近似模拟公式：

(15)

公式(15)中，θ表示人眼观察方向和半角矢量的夹角；η1和η2分别是两种不同介质的折射率.公式(15)表示光从具有不同折射率的电介质材料之间穿过后的相互作用.

本文设定光线射出介质为空气或真空，借鉴文献[22]提出的近似模拟公式，见公式(16)：

(16)

公式(16)中，θ表示人眼观察方向和半角矢量之间的角度，η表示材质的折射率，k表示材质的吸收系数.图3是运用导体菲涅尔系数得出的实验效果图.

在图3(a)和图3(b)中，η的取值从最左边开始，分别为1.2、1.5、1.8和2.4.图3明显的可以发现，两个函数得到的效果差异较大，由于太阳光照射到云的表面，穿过的介质是大气，所以本文中Cook-Torrance BRDF光照模型采用导体菲涅尔系数.

图3 导体菲涅尔系数近似实验Fig.3 Conducting Fresnel coefficient approximation experiment

因此，Cook-Torrance BRDF光照模型可由式(17)近似表示：

(17)

3.4 本文算法整体实现过程

本文使用了基于GPU的Ray Casting算法绘制三维云，并为了绘制的真实感，结合了纹理抖动和BRDF.首先生成3D柏林噪声纹理，与实测气象体数据纹理进行不同比例的融合，从而实现对云表面的随机抖动效果，得到采样值；然后将该采样值与新的光照模型Cook Torrance BRDF相结合，对颜色和透明度值在射线前进方向进行累计，增加光照的效果.最终绘制出来的三维云相比较先前的效果，云显得更加柔和、更立体和真实.本文改进云绘制算法的详细描述如算法1所示.

算法1.

输入：通过气象观测得到的预报数据，本文使用相对湿度预报格点数据作为样例，记作Vrh

输出：采样点的颜色值记作C，不透明度记作O

步骤1.将多层二维网格数据Vrh，映射到灰度图上并做归一化处理；

步骤2.绘制六面体包围盒，将所使用的气象体数据源Vrh中的网格坐标转化为三维球体上的经纬度；

步骤3.通过光栅化求交计算光线穿过体纹理时和距离视点最近以及最远的两个面的交点，即入点和出点，从而相减得到光线传播方向；

步骤4.将多层相对湿度格点数据Vrh做成的体纹理texVolume_，利用3D Perlin噪声函数生成的噪声纹理texNoise3D_以及在配置文件里读取的颜色数组做成的纹理texTF_一起传入着色器中，其中3个纹理大小一致；

步骤5.沿光线传播方向每隔相同距离进行采样，将从texVolume_采样得到的Valuesrc，和从texNoise3D_采样得到的Valuejitter做线性叠加，并用调节因子rf控制噪声对采样结果的影响，得到最终的采样值Value=Valuesrc×(1-rf)+Valuejitter×rf，还可以通过调整采样间距来控制三维云的细腻程度；

步骤6.根据前面步骤所得到的采样值从texTF_中获取颜色值Csrc，再计算环境光ambient和漫反射光diffuse的效应，并且根据Cook-Torrance光照模型计算出镜面反射效应specular，将ambient、diffuse和specular这3者叠加起来后与Csrc关联，得到相应采样点的颜色值C；

步骤7.对采样点的颜色值C和不透明度O进行累计，并设一个阈值来控制不透明度，本实验中将其值设为0.99，为了进一步提高本实验的效率，本文还设定了采样最大步数的累计值maxSteps来提前终止采样，以防重复采样.

4 实验结果与分析

本实验使用实际观测得到的相对湿度数据为输入，基于GPU绘制三维云景，并且使用了提前光线终止和光栅化求交等加速方法优化得到的效果如图4所示，图4(a)是绘制在3D数字球体上的全局图，图4(b)是相应局部的放大图.

图4 实验所绘三维云效果展示Fig.4 3d cloud effect is shown in the experiment

从图4可见，使用相对湿度数据绘制出的云景，云面没有太多层次感，且缺乏三维感，不够真实.

图5中，图5(a)为没有应用噪声纹理的三维云绘制效果及局部扩大图，图5(b)、图5(c)和图5(d)分别表示在调节因子rf为0.05，0.1，0.15下噪声纹理应用至采样过程的云景绘制效果及局部扩大图.比较图5(a)和图5(b)，可以看到图5(a)中云表面的细节没有图5(b)中云的丰富，但图5(b)中噪声干扰效果并不明显.观察比较图5(b)、图5(c)和图5(d)可以看出，随着调节因子增大，噪声纹理对采样结果值以及三维云绘制效果的影响也不断增强.从图中左下角的局部放大图来看，随着调节因子的逐渐增大，云块的边缘愈加的清晰，云面也愈加的柔和，真实感也随着愈加的强烈了.

图5 不同大小调节因子所绘三维云Fig.5 Rendering effects of different adjustment factors

比较5(d)和图5(e)可以发现，如果调节因子一直增大，直到一定程度时，云面虽然十分柔和，但会丢失不少细节(图中框住的小块区域放大至右上角来看)，导致整体效果反而更缺乏真实性.如当调节因子增加到0.2时(如图5(f))，三维云表面的细节几乎完全丢失.所以应用噪声纹理绘制云景时必须控制好调节因子，才能既保持了较强的真实感，又可以不丢失云面的细节.

本文实验在绘制云过程中加入了随机抖动，对三维云应用了3种光照模型，包括了2种经典光照模型，还有本文中提到的Cook-Torrance BRDF.实验过程中，使用了361×181的规则网格数据，经纬密度间隔为1.0；材质折射率为2.4吸收系数为1.2.截取的三维云局部效果图如图6所示.

图6 3种光照模型效果对比Fig.6 Comparison of rendering effects of three illumination models

图6(a)为没有添加任何光照模型的效果，图6(b)为添加了Blinn Phong光照模型的效果，图6(c)为添加了Phong光照模型的效果、图6(d)为添加了Cook-Torrance BRDF光照模型的效果.对比图6(b)、图6(c)和图6(a)，无光照模型采用单一颜色赋值，色彩层次不足.Phone光照模型根据其反射光线方向和观察发现方向的夹角计算镜面瓣，因此在光照方向和观察方向一致时具有伪影(图6(c)中框选位置).Blinn Phong光照模型利用中间向量和观察方向计算镜面瓣，消除了这种伪影.但是Phone和Blinn Phone光照模型都是基于经验性的，Cook-Torrance BRDF光照模型基于物理的特性使其能够正确的计算光照，而不会引入伪影.可以发现Cook-Torrance BRDF光照模型应用的三维云绘制效果，与应用Phong模型所绘制出来的感觉十分相似，在云的细节上有很好的体现并且所绘制的云景不会给用户太过于艺术的感官体验；除此之外，Cook-Torrance BRDF光照模型还和Blinn Phong模型一样有强立体感.

本文实验在Windows10的Visual Studio2010的编程环境下实现，相关计算机硬件配置如下：8核Intel(R)Core(TM)i5-8300H的处理器，8.0GB内存，GTX 1060的GPU以及4.0GB显存.实验中使用了phong、Blinn Phong以及Cook-Torrance BRDF这3种不同光照模型分别以同一尺度的相对湿度进行绘制，可以发现3种模型下绘制的计算量相差不大，但本文提出的Cook-Torrance BRDF光照模型的绘制效率略有提升，且绘制效果也有一定改善.

5 结 语

本文使用实际探测到的气象数据，首先利用Ray Casting算法绘制三维云.但因为网格数据在垂直方向上比较少，导致所绘制的云效果不佳，立体感不够，在细节上较粗糙.为了改善这些情况，本文将3D Perlin噪声做成纹理对绘制的云面产生随机抖动，并在生成3D Perlin噪声时，加入Perlin的两点改进意见，且为了使插值运算得到的结果更加的平滑，使用了余弦插值来优化.从结果可以看到，纹理噪声抖动之后的云面看起来更加柔和，也更具有真实感.又考虑到云在阳光充足时看上去会更明亮，因而本文在绘制过程中，尝试了不同的光照模型来模拟，并对其产生的效果进行了对比.从结果可以看到Phong光照模型和Blinn Phong光照模型艺术感过强且立体感不够.而本文提出的新光照模型，Cook Torrance BRDF光照模型充分考虑到材质的各向异性，可以使得绘制出的云更具真实感及立体感，且绘制速度并没有降低.