烟雾扩散与气体污染的动态仿真研究

唐 勇，甄志华，汪新宇，孙旭东

(1.燕山大学 信息科学与工程学院，河北 秦皇岛 066004；2.河北省计算机虚拟技术与系统集成重点实验室，河北 秦皇岛 066004)

0 引言

气体污染是生活中常见的污染之一，包括燃料燃烧、尾气排放、工业污染等，其中，烟雾扩散导致的气体污染对日常生活造成了极大的影响，因此，模拟气体污染具有重要的现实意义和广泛的应用前景。一方面，由于烟雾随时间变化扩散较快，烟雾细节难以捕获；另一方面，烟雾扩散同气体污染难以建立联系，值得深入研究。

近年来，烟雾模拟一直都是国内外研究的焦点话题，自2003年Stam[1]和Fedkiw等人[2]提出模拟烟雾运动的方法以来，烟雾的模拟应用越来越广泛。采用物理模型的方法通常能真实地展现烟雾细节，Xie等[3]通过将低分辨率烟雾进行神经网络训练，合成具有真实细节的高分辨率烟雾，并将其应用于不同的物理模型，但是其神经网络训练时间较长，无法达到实时；唐勇等[4-5]提出一种改进的空间自适应漩涡限制方法，生成清晰漩涡烟雾细节；在烟雾路径方面，通过采用改进的有限差分法求解N-S方程，加快求解速度，并引入吸引力和驱动力，实现大规模烟雾路径模拟。而基于经验模型的模拟烟雾以及气体污染，虽然能迅速达到实时状态，但牺牲烟雾细节，真实感不足。陆薇等[6]提出了球形雾化渲染模型，结合传统Perlin噪声，通过HDR实时绘制层次雾；Guo等[7]采用Perlin噪声生成异质密度分布纹理，再利用MRF模型结合大气散射透视图评估渲染，绘制异质雾图像。

因此，本文提出一种混合经验模型、物理模型的方法，绘制由烟雾扩散引起的动态气体污染。首先，针对物理模型优化烟雾运动轨迹，并采用基于物理渲染的方式对烟雾粒绘制，提升扩散细节；其次，根据大气污染理论建立烟雾扩散同气体污染之间的联系，并根据优化后的经验模型，模拟具有真实感、动态可变的气体污染。

1 基于物理模型的烟雾模拟

为模拟真实的烟雾扩散，改进半拉格朗日法中的驱动力，并优化漩涡计算过程，结合基于物理的光照模型生成真实烟雾效果。

1.1 N-S方程构建物理模型

采用半拉格朗日法模拟烟雾粒子能更加真实地计算烟雾扩散运动，其动量守恒方程和质量守恒方程为

(1)

∇·u=0，

(2)

其中，u为不可压缩流体速度场，ρ为流体密度，p为压强，v为黏度系数，f表示流体合力项，∇表示梯度算子。

在流体合力项中，通过引入漩涡约束力增加烟雾模拟细节，漩涡场的计算公式为

ω=∇×u，

(3)

其中，×表示卷积操作。将生成的漩涡场代入式(1)得到涡度守恒方程

(4)

根据Biot-Savart公式，通过漩涡场带动速度场更新，实现真实烟雾粒子物理运动。

1.2 引入K-D树提升计算效率

由于引入漩涡后，速度场迭代会受多涡影响，计算量显著增大。通过引入K-D树来降低计算量。相比八叉树，K-D树在空间划分方面具有明显的优越性，模拟效率显著提升。

根据K-D树分割点设置距离阈值，如果粒子到分割点距离小于阈值，计算距离阈值内每个漩涡对粒子的速度场影响。单漩涡对环境影响的计算公式为

(5)

为简化计算，将漩涡视为点，由积分转化为求和方式，速度场可近似表示为

(6)

在计算过程中，若粒子到漩涡的距离大于阈值，则视涡旋为单个涡旋结构，而小于阈值的多个涡旋簇，采用矢量和计算涡旋位置L，

(7)

(8)

1.3 改进外力项提升运动细节

计算外力项过程中，仅靠漩涡力并不能提供更多的细节，本文通过改善风力提升烟雾扩散轨迹。

空间中的风通常由顺风向Vm，横风向Vw以及竖风向Vh组成。采用Kaimal谱作为风速谱表达式，比与高度无关的Davenport功率谱能更好展现风场的自相关特性。因此风速计算中脉动风通过高斯过程求解，得到v(t)表达式

(9)

其中，p为AR模型的阶数，Δt为模拟风速时程的时间步长，ψk为AR模型的自回归矩阵系数，N(t)为独立随机过程向量。根据随机振动理论，利用维纳-辛钦公式求得相关函数后，利用期望操作得到RN，具体计算公式为

(10)

对RN进行Cholesky矩阵分解，可求得具有时间间隔的随机风速向量，将其带入风力项改善烟雾运动轨迹。

1.4 烟雾粒子真实感渲染绘制

烟雾作为非均匀介质媒体，光线穿越其中会产生散射和反射现象，如图1所示。考虑到实时性问题，将真实烟雾纹理结合基于物理的双向反射分布函数渲染更真实的烟雾粒子。

图1 光线的散射和反射Fig.1 The scattering and reflection of light

该模型中，出射辐射率L0(v)等于所有入射方向的辐射率积分和BRDF值，以及余弦值的乘积，出射辐射率的计算公式为

(11)

反射项通常由次表面散射和反射组成。次表面散射的计算公式为

(12)

其中，F=0.5+2r(hI)2。

采用Torrance-Sparrow微面元模型计算高光项

(13)

其中，F(I,h)为菲涅尔反射模型，用于处理反射光和入射光比率，G(I,v,h)为阴影遮掩函数，D(h)为法线分布函数。

2 气体污染动态模拟

考虑到实时性，采用基于经验模型的方法模拟气体污染，并建立烟雾扩散同气体污染之间的联系，绘制真实动态的污染。

2.1 构建烟雾与气体污染关系模型

污染物高度由烟云抬升高度同颗粒物高度组成，其中烟云抬升高度根据烟雾粒子高度赋值，颗粒物高度采用线性高度场，为搭建烟雾扩散与气体污染之间的关系，将烟雾建模同大气污染理论模型结合。本文引入高斯烟羽污染模型

(14)

其中，X(x,y,z,t,H)为时间t下，竖风向x米，横风向y米，距离地面z米的浓度，Q代表污染源强度，H为粉尘的有效高度，u为脉动风速，σy、σz为颗粒物的水平和垂直扩散系数。增加污染因子k并优化为指数函数，实现污染浓度渐变效果。针对不同污染变化，渲染成不同程度的气体污染。

2.2 优化污染衰减公式增加污染细节

考虑到实时性，利用基于经验模型的屏幕技术进行绘制。首先，根据污染浓度衰减公式获得浓度变化梯度，计算公式为

(15)

由于气体污染扩散的不规则性，式(15)渲染产生的细节较少，为此采用基于3D Perlin噪声优化后的分型噪声，并在其中引入风速因子，生成动态变化的气体污染效果，具体计算公式为

(16)

其中，m为倍频数，f为频率，p为振幅，vp(x,y,z,t)代表P处的脉动风场，生成真实的污染衰减公式。

2.3 动态计算环境光

优化生成的污染浓度衰减公式，虽然细节得以提升，但颜色无法随时间动态变化。为解决该问题，引入优化后的时间轴算法，生成动态可变的烟雾，不同时刻环境光的计算公式为

(17)

其中，Cinc表示环境光的插值颜色，T0表示初始时间(即凌晨初始时间)，Tcurrent表示当前时间，Ttotal表示总体时间。根据负余弦函数平滑过渡，计算不同时刻的环境光CTcurrent。

另外，引入透明度系数λ，并将污染颜色同光照颜色作点乘操作，获得最终的真实光照，如公式(18)所示，进而渲染出真实的动态气体污染场景。

C=f(CTcurrent+Clight)+λ(1-f)Cf

(18)

3 烟雾扩散与气体污染仿真实现

3.1 本文整体程序框架

本文整体程序流程图如图2所示。

图2 整体程序流程图Fig.2 Overall program flow chart

3.2 实验结果与分析

本文实验采用基于Windows系统，Unity3d平台开发的动态气体污染仿真系统，硬件环境为：Intel Core i7-4790 CPU 3.60 GHz，16 G RAM，显卡为NVIDIA GeForce GTX 750Ti。

图3(a)、3(b)为文献[3]与本文方法烟雾扩散对比。图3(a)为文献[3]中使用时间相关的生成模型来解决流体流动问题的实验效果图，网格分辨率为256×180×180，虽然仅使用了单个时间步长，但帧率仅有0.000 8 fps，无法实时渲染。本文方法能在实时基础上，生成大量扩散细节。图3(c)、图3(d)为文献[4]与本文方法烟雾对比。图3(c)为文献[4]采用改进的空间自适应漩涡限制方法模拟烟雾，网格分辨率为64×86×64，烟雾扩散细节不够明显。本文利用纹理结合物理渲染模型的方法模拟的烟雾色彩、光照更加逼真。图3(e)、图3(f)为真实图片与本文方法应用至大规模场景的对比，保证实时的状态下，模拟烟雾真实自然。

图3 与文献[3]、[4]以及真实烟雾对比实验Fig.3 Comparison with literature [3],[4] and real smoke

图4展示为导弹烟雾尾迹扩散对比实验。其中图4(a)为真实图片，图4(b)为文献[5]方法，图4(c)为采用本文方法应用到烟雾路径方向。通过对比看出，本文利用优化的物理模型方法比文献[5]生成的更自然灵活，尾迹扩散效果更明显，更接近真实图片效果。

图4 导弹烟雾尾迹对比实验Fig.4 Missile smoke wake comparison experiment

图5为优化后的经验模型模拟的气体污染，通过将不同污染因子代入气体污染模拟效果进行对比。图中分别为污染因子为k=0、0.4、1下的气体污染模拟，可以明显比较出不同污染因子产生不同梯度的气体污染，效果明显。

图5 污染因子对气体污染影响实验Fig.5 Experiment on the influence of pollution factors on smoke pollution

图6展示了清晨T=6、正午T=13、夜晚T=20时刻下，烟雾扩散以及气体污染随着时间轴的动态变化，不同时刻光照插值计算颜色所展现出的气体污染明暗效果。

图6 优化的时间轴算法实验Fig.6 Optimized timeline algorithm experiment

为了检测动态污染模拟效率，表1列出了本文实验以及部分文献数据的对比状况。表中粒子和网格指文献中模拟烟雾效果使用的粒子数或网格大小。其中，可以明显看出本文算法在显著提升模拟的细节同时，保证了实时性。

表1 不同实验场景帧率统计Tab.1 Frame rate statistics for different experimental scenes

4 结论

本文提出一种烟雾扩散动态污染模拟的方法。首先，采用半拉格朗日方法计算粒子运动轨迹，并使用K-D树提升计算效率，引入基于Kaimal谱的脉动风模型并结合真实物理光照，在避免粒子颗粒感的同时，改善烟雾扩散细节；此外，将烟雾扩散同优化的高斯烟羽模型结合，利用改进的经验模型以及Perlin分型噪声生成更加真实的气体污染；采用优化的时间轴算法，解决了污染颜色无法动态变化的问题，大幅度提升污染真实感。实验数据表明，本文方法能实现烟雾扩散下的气体污染实时模拟。在未来工作中，需要进一步对气体污染与环境的交互进行研究。