一种基于物理过程的红外积云的建模方法与流程

[0001]
本发明属于红外波领域，具体涉及一种基于物理过程的红外积云的建模方法。


背景技术：

[0002]
红外场景仿真是分析和研究红外场景的有效手段，依赖于精确的三维建模和严格的物理过程模拟，使得红外虚拟场景有很强的真实感，可以满足很多军事需求和民用需求。仿真模拟红外空中场景需要对目标和场景进行分析，尤其是一些常见的空中景象，包括云、阳光、雨雪等。云作为典型的空中景象，在飞行模拟和户外场景仿真中经常出现，但由于它们的光学特征、形状等外观表现与高度、气候、天气、阳光以及观察角度等有关，导致云的仿真相当复杂且困难。现有技术通过采用三维空间内堆积球形粒子的建模方法获得具有层次感的积云形状，并将其结果与改进噪声建模方法的结果作比较，在保证渲染效率的基础上对两种方法的优缺点进行分析。由于现有的云的仿真体系中三维模型大多沿用比较传统的建模方法，对云中的物理过程缺乏准确的描述，并且在光照模型方面，采用较为简单的前向散射模型，而云中存在的多种光照现象，比如二次散射、多次散射效果。在整体上，现有的云仿真体系在综合三维建模与光照建模后，难以做到真实性、实时性的兼顾，并且大多适用于可见光领域，缺乏红外波段适配的光照模型。


技术实现要素：

[0003]
为了解决现有技术中存在的上述问题，本发明提供了一种基于物理过程的红外积云的建模方法。本发明要解决的技术问题通过以下技术方案实现：
[0004]
本发明实施例提供的一种基于物理过程的红外积云的建模方法，包括：
[0005]
s1，获取积云的三维网格模型；
[0006]
其中，所述三维网格模型中每个网格存储多个粒子的第一变量信息，所述粒子表示在所述网格内各个位置处的云，所述变量信息包括：液滴密度、温度以及速度；
[0007]
s2，基于液滴密度、温度以及速度，确定按照液滴密度、温度以及速度扩散后的积云；
[0008]
s3，基于在光照照射所述积云的入射方向与散射方向之间夹角，确定散射相函数；
[0009]
s4，基于积云中滴液的尺寸分布的宽度、滴液的尺寸分布的特征半径，确定液滴尺寸分布的有效相函数；
[0010]
s5，针对处于不同波段的积云，根据处于不同波段的散射效率因子，计算平均散射效率因子；
[0011]
s6，将处于不同波段的积云的散射相函数求平均，确定平均散射函数；
[0012]
s7，基于平均散射函数以及平均散射效率因子，当光线入射所述积云后，沿着色器的视线方向以及光线方向，计算光线在积云内多次散射后的第一辐射能量；
[0013]
s8，针对处于远红外波段以及中红外波段的积云，分别计算该积云对应的自发辐射能量；
[0014]
s9，将所述自发辐射能量与第一辐射能量求和，确定所述积云在视线方向上的总能量；
[0015]
s10，将所述总能量分别与图像三原色通道中每个通道的最大像素值相乘，获得积云在视线方向上像素值；
[0016]
s11，按照所述像素值，扩散后的积云在该在视线方向上位置处进行渲染，并重复执行s8至s11，直至渲染完扩散后的积云每一处位置。
[0017]
本发明实施例通过积云的三维网格模型，计算网格中存在的水汽输送、相变特征等过程，得到云的液滴密度、水汽密度、温度等特征的分布，然后，多重散射简化模型，对积云的可见光模型中存在的光照现象按照光照方向以及视线方向进行分类并计算，在确定散射相函数和基础散射理论后，计算多次散射的辐射能量，并计算红外场景中积云的自发辐射能量，从而得到总能量，对可见光模型的辐射能量进行修正，设计自发辐射在不同波段下的计算过程得到积云红外光照模型，然后对云的光照进行实时渲染。从而兼顾红外波段以及非红外波段积云，使得对积云渲染的真实性、实时性更高，同时提高对积云渲染的准确性。
[0018]
以下将结合附图及实施例对本发明做进一步详细说明。
附图说明
[0019]
图1是本发明实施例提供的一种基于物理过程的红外积云的建模方法的流程示意图；
[0020]
图2是本发明实施例提供的实现视线方向以及光线方向的积分，循环处理红外波段的自发散射计算以及散射计算，从而得出辐射能量的示意图；
[0021]
图3是本发明实施例提供的实现视线方向以及光线方向的积分，循环处理红外波段的自发散射计算以及散射计算，从而得出辐射能量的细节流程图；
[0022]
图4是本发明实施例提供的三维网格模型的示意图；
[0023]
图5是本发明实施例提供的半拉格朗日法求解速度、温度、密度因流体速度场产生的输运效应的流程图；
[0024]
图6是本发明实施例提供的积云扩散的流程图；
[0025]
图7是本发明实施例提供的初始积云的示意图；
[0026]
图8是本发明实施例提供的积云的光照示意图；
[0027]
图9是本发明实施例提供的从视线方向以及光线方向设置采样点，计算总能量的示意图；
[0028]
图10是本发明实施例提供的在红外场景中积云的红外光照仿真结果。
具体实施方式
[0029]
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
[0030]
实施例一
[0031]
如图1所示，本发明提供的一种基于物理过程的红外积云的建模方法，包括：
[0032]
s1，获取积云的三维网格模型；
[0033]
其中，所述三维网格模型中每个网格存储多个粒子的第一变量信息，所述粒子表示在所述网格内各个位置处的云，所述变量信息包括：液滴密度、温度以及速度；
[0034]
s2，基于液滴密度、温度以及速度，确定按照液滴密度、温度以及速度扩散后的积云；
[0035]
s3，基于在光照照射所述积云的入射方向与散射方向之间夹角，确定散射相函数；
[0036]
s4，基于积云中滴液的尺寸分布的宽度、滴液的尺寸分布的特征半径，确定液滴尺寸分布的有效相函数；
[0037]
s5，针对处于不同波段的积云，根据处于不同波段的散射效率因子，计算平均散射效率因子；
[0038]
s6，将处于不同波段的积云的散射相函数求平均，确定平均散射函数；
[0039]
s7，基于平均散射函数以及平均散射效率因子，当光线入射所述积云后，沿着色器的视线方向以及光线方向，计算光线在积云内多次散射后的第一辐射能量；
[0040]
s8，针对处于远红外波段以及中红外波段的积云，分别计算该积云对应的自发辐射能量；
[0041]
s9，将所述自发辐射能量与第一辐射能量求和，确定所述积云在视线方向上的总能量；
[0042]
s10，将所述总能量分别与图像三原色通道中每个通道的最大像素值相乘，获得积云在视线方向上像素值；
[0043]
其中，对总能量需要进行量化，使得量化后的总能量值位于0-1区间。
[0044]
s11，按照所述像素值，扩散后的积云在该在视线方向上位置处进行渲染，并重复执行s8至s11，直至渲染完扩散后的积云每一处位置。
[0045]
参考图2，在图2中在获得三维模型后，需要再实现视线方向以及光线方向的积分，循环处理红外波段的自发散射计算以及散射计算，从而得出辐射能量，对辐射能量灰度化进行渲染，后续将对每个步骤的公式以及原理进行更详细的介绍。
[0046]
本发明实施例通过积云的三维网格模型，计算网格中存在的水汽输送、相变特征等过程，得到云的液滴密度、水汽密度、温度等特征的分布，然后，多重散射简化模型，对积云的可见光模型中存在的光照现象按照光照方向以及视线方向进行分类并计算，在确定散射相函数和基础散射理论后，计算多次散射的辐射能量，并计算红外场景中积云的自发辐射能量，从而得到总能量，对可见光模型的辐射能量进行修正，设计自发辐射在不同波段下的计算过程得到积云红外光照模型，然后对云的光照进行实时渲染。从而兼顾红外波段以及非红外波段积云，使得对积云渲染的真实性、实时性更高，同时提高对积云渲染的准确性。
[0047]
实施例二
[0048]
作为本发明一种可选的实施例，所述s2的步骤包括：
[0049]
s21，设置初始积云；
[0050]
s22，针对初始积云中每一个网格，初始化该网格内粒子的液滴密度、温度以及速度；
[0051]
s23，将所述网格中的每个粒子按照所述初始化后的速度移动；
[0052]
s24，基于粒子的风力、浮力以及涡流抑制力，计算移动后的每个粒子的外力；
[0053]
s25，基于积云当前的大气压强，计算移动后粒子的液滴密度以及温度；2.26
[0054]
s26，基于所述外力对已获取流体运动方程进行简化，获得简化后的流体运动方程；
[0055]
s27，对所述简化后的流体运动方程进行求解，获得所述粒子当前单位时间下的速度；
[0056]
s28，将所述当前单位时间下的速度替换所述初始化的速度，重复执行s23至s27的步骤，直至达到迭代次数，获得扩散后的积云。
[0057]
实施例三
[0058]
作为本发明可选的一种实施方式，步骤s26包括：
[0059]
获取纳维-斯托克斯方程作为流体运动方程；
[0060]
设置中间速度场以及时间分割的最小单位，对所述流体运动方程进行简化。
[0061]
实施例四
[0062]
作为本发明可选的一种实施方式，所述s4的步骤包括：
[0063]
获取散射相函数；
[0064]
基于水滴中总体粒子数，水尺寸分布的宽度，伽马分布函数以及尺寸分布的特征半径，尺寸分布对应的有效半径，确定不同液滴尺寸分布；
[0065]
基于所述不同液滴尺寸分布，对不同液滴尺寸的散射相函数进行加权和求和，获得液滴尺寸分布的有效相函数。
[0066]
实施例五
[0067]
作为本发明可选的一种实施方式，s8的步骤包括：
[0068]
针对处于远红外波段的积云，使用普朗克公式计算该积云对应的自发辐射能量；
[0069]
针对中红外波段的积云，基于已设置的辐射发射率，计算自发辐射能量。
[0070]
实施例六
[0071]
作为本发明可选的一种实施方式，在所述s8的步骤之前，所述建模方法还包括：
[0072]
根据光吸收率以及能量守恒定律，确定处于远红外波段的光线的发射率以及中红外波段的辐射发射率。
[0073]
实施例七
[0074]
作为本发明可选的一种实施方式，s7的步骤包括：
[0075]
s71，设置光线参数；光线参数包括：光线投射方向、光照方向，光线投射步长以及光线投射的最大次数；
[0076]
s72，沿着视线方向，按照光线投射步长确定视线方向采样点；
[0077]
s73，确定在所述视线方向当前采样点处的第二变量信息；
[0078]
s74，基于所述第二变量信息，计算视线方向的透过率；
[0079]
s75，沿着当前采样点的光线方向，按照所述光线投射步长确定在所述光线方向的当前采样点处的第二变量信息；
[0080]
s76，基于所述第二变量信息，计算光线方向当前采样点的透过率；
[0081]
s77，对所述透过率进行累计；
[0082]
s78，重复执行s75至s76的步骤，直至计算完成光线方向的所有采样点的透光率；
[0083]
s79，计算沿光线方向，光经过散射作用的一次散射率以及二次散射率；
[0084]
s80，基于所述一次散射率以及二次散射率，计算第二辐射能量；
[0085]
s81，将所述第二辐射能量分别进行累计；
[0086]
s82，将沿视线方向的下一个采样点作为当前采样点，并重复执行s73至s82的步骤，直至轮询完所述视线方向的采样点；
[0087]
s83，沿着所述光线方向，计算多次散射后的辐射比例；
[0088]
s84，基于所述辐射比例以及累计后的第二辐射能量，计算第一辐射能量。
[0089]
参考图3，在图3中显示了轮询视线方向采样点以及光线方向采样点的全部流程，是实施例七的展示流程，具体过程与实施例七的实施过程相同。
[0090]
实施例八
[0091]
作为本发明可选的一种实施方式，s78的步骤包括：
[0092]
基于液滴的有效半径和液滴的粒子数密度，确定消光截面；
[0093]
基于所述消光界面，确定透光率。
[0094]
实施例九
[0095]
作为本发明可选的一种实施方式，s79的步骤包括：
[0096]
基于光线经过积云的云层时的透过率、视线方向以及光线方向，确定所述光线经过积云的投射作用函数；
[0097]
确定光线经过积云的一次散射作用函数；
[0098]
确定光线经过积云的二次散射作用函数；
[0099]
使用所述一次散射作用函数，计算一次散射率以及使用所述二次散射作用函数，计算二次散射率。
[0100]
实施例十
[0101]
作为本发明可选的一种实施方式，s83的步骤包括：
[0102]
当积云处于热平衡状态以及处于多次散射时，确定光线的透过率、一次散射比例、二次散射比例以及多次散射比例的能量守恒方程；
[0103]
当多次散射是均匀时，基于所述能量守恒方程确定单个面元的辐射比例；
[0104]
当多次散射是无方向时，基于所述能量守恒方程确定单个面源上不同方向的辐射比例。
[0105]
下面对本发明实施例中的基础原理以及实施例一至实施例十中步骤中的详细描述。
[0106]
首先描述参数含义。
[0107]
u
表示速度，t表示时间，ρ表示密度，p表示压强，f表示力的作用，u
*
表示中间速度场，θ表示位势温度，q表示粒子密度，t表示开尔文温度，п表示exner函数，p0表示标准大气压，r
d
表示干空气的气体常数，c
p
表示干空气的恒压比热容，c
v
表示干空气的容积比热容，f
wind
表示风力作用，f
buo
表示浮力作用，f
vc
表示涡流抑制力，θ
v
表示虚位势温度，θ
v0
表示参考虚温度，q
c
表示液滴密度，q
v
表示水汽密度，ω表示流体速度场的旋度，l表示网格单位长度，q
vs
表示饱和水汽密度，z表示海拔高度，γ表示气温垂直递减率，l表示水汽化潜热，θ
b
表示中间位势温度，p(θ)表示散射相函数，k表示波数，csca表示散射截面，x表示矢量散射振幅，n(r)表示液滴尺寸分布，γ表示尺寸分布的宽度，γ(γ)表示伽马分布函数，n0表示水滴中总体粒子数，r
n
表示尺寸分布的特征半径，r
e
表示有效半径，p(θ)表示有效相函数，τ(d)表示
透过率，表示视线方向，表示光线方向，q表示散射效率因子，σ表示消光截面，κ表示消光系数，s0表示透射作用，s1表示一次散射作用，s2表示二次散射作用，表示总体透过比例，表示总体一次散射比例，表示总体二次散射比例，表示总体多次散射比例，表示平均有效相函数，m
λbb
表示光谱辐射出射度，c1表示第一辐射常数，c2表示第二辐射常数，α表示吸收率，ε表示发射率，g表示量化灰度。
[0108]
在实施例一中，三维模型包含了三个维度的信息，包含了xyz轴三个方向的信息。为了方便快捷地使用和存取这些信息，需要设计合适的模型。本发明实施例采用一个三维网格来保存这些信息，存储云中的密度、温度及速度等信息，结果保存在一个三维贴图中，以方便gpu计算过程中的使用和存取。
[0109]
如图4，本章中的三维模型是在一个左手坐标系下建立的，其中xyz轴正向如图所示。图中的正方体为存储云模型信息的三维网格，这个网格被划分成均匀的若干份，图中为4*4*4网格，实际的计算过程中网格精细度远大于此。被划分后的每一个均匀小格称为体素(voxel)，每个体素存储了当前位置的云信息。云产生在网格中，为了控制云的形状和观察云的变化，通常云的体积小于网格的体积。当网格划分的足够精细时，云体积内的体素便可以表示云中的粒子，存储当前位置的液滴信息与水汽信息。
[0110]
图4中的坐标系也是本发明实施例仿真虚拟场景的坐标系，观察相机一般位于x轴正方向的轴上，观察方向(-1,0,0)。本发明实施例的一些方位上的描述是相对于云网格而言的，云网格处于坐标轴原点。其中“前”指的是相对于云网格的轴正方向，“上”指的是相对于云网格的y轴正方向，“右”指的是相对于云网格的z轴正方向，“后”“下”“左”分别与之对应。
[0111]
下面详述实施例二步骤的实施过程以及过程中需要的参数。
[0112]
(1)不可压缩条件
[0113]
真实流体都有程度不同的可压缩性，在空气动力学中，气体的密度变化是否可忽略，要根据气体流动的速度来确定。一般来说，当气体流动速度低于0.3马赫，适用低速空气动力学方法，此时气体可以认为不可压缩；当气体流动速度高于0.3马赫时，适用高速空气动力学方法，此时气体应当认为是可压缩流动的。
[0114]
(2)粘性力可忽略
[0115]
雷诺数(reynolds number)是一种可用来表征流体流动情况的无量纲数，当雷诺数小于等于1，即流体粘度很小而相对滑动速度又不大时，绕流物体边界层外，粘性应力是很小的，方程中粘性项可以忽略，即可看成理想流体。对于自然条件下的水、空气等，由于运动速度较小、密度较低，因此满足这一条件，可以忽略粘性力。
[0116]
本发明实施例二中的流体运动方程以及简化过程做出详述。
[0117]
(1)流体运动方程
[0118]
纳维-斯托克斯方程(navier-stokes方程，简称n-s方程)，为描述粘性不可压缩流体动量守恒的运动方程。由于流体扩张速度远低于音速，可以忽略流体的粘力项，这样n-s方程可以简化为欧拉方程，形式如下：
[0119]
[0120]
同时有连续性条件：
[0121][0122]
式中
u
为速度，ρ为密度，t为时间，p为压强，f为力的作用效果，式(2-1)是流体动量守恒的表现，右一项为流体速度场自身的运动更新，第二项为流体压力梯度引起的加速度，最后一项为外力作用引起的加速度，通常包括重力、风力等。式(2-2)为连续性方程，说明流体速度场无散度，保证质量守恒。
[0123]
在实时条件下完全解出欧拉方程是不现实的，因此需要采用一些近似手段。本发明实施例将式(2-1)分成两个部分，第一步先计算一个中间速度场u
*
，并假设时间分割的最小单位为δt，有：
[0124][0125]
对于式(2-3)，等号右边第一项为流体速度随流体运动产生的输运结果，使用半拉格朗日法可以近似解出，同样的方法也可以用来更新其它常量，比如温度和水密度(这里兼指体素内包含的水汽和液滴)：
[0126][0127][0128]
其中θ为流体的位势温度，q代表粒子的密度(气态或液态)。位势温度在绝热近似的高度变化下是恒定的，可以通过计算得到：
[0129][0130][0131]
式(2-6)中t为温度，单位为k，п为exner函数，可以视为气体所受无量纲压力。p0为海平面标准大气压，p为所处高度的压强。r
d
为干空气的气体常数(287jkg-1
k-1
)，c
p
和c
v
分别为干空气恒压比热容和容积比热容。
[0132]
(2)水汽输运效应
[0133]
半拉格朗日法求解速度、温度、密度因流体速度场产生的输运效应的过程主要分五步，如图5所示：
[0134]
以温度为例，首先获取待更新位置pos的温度及速度v，向后倒回半个步长计算出中间位置midpos＝pos-v*0.5*dt，获取中间位置midpos的速度v
mid
，向后倒回一个步长计算出最终位置finalpos＝midpos-v
mid
*dt，最后用finalpos位置的温度来更新pos位置的温度。这样，就完成了一次温度随流体速度场的输运。
[0135]
(3)力的作用和修正
[0136]
完成了速度随速度场的输送后，考虑外力作用引起的加速度f，通常外力作用一般包括以下三项，风力f
wind
、浮力f
buo
以及涡流抑制力f
vc
，因此有式(2-8)。
[0137]
f＝f
vc
+f
buo
+f
wind
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-8)
[0138]
对于风力f
wind
，通常相对于云的仿真是独立的，因此可以直接进行设置或采取仿
真条件下的典型风速。浮力f
buo
的计算包含了重力，实际上是对一个体素内气团竖直方向的受力分析结果。气团中水汽占比越高，浮力越大；液滴占比越高，重力越大，浮力也就越小。另外，浮力的计算还和所处环境的温度、压强有关。具体计算如下：
[0139][0140]
式中g为重力加速度，为竖直向上的单位向量，q
c
为液滴密度，θ
v0
是参考虚温度，通常取值在290k到300k之间。θ
v
为虚位势温度，代表同样密度、压强条件下干空气的位势温度，计算方法如下：
[0141]
θ
v
＝θ(1+0.61q
v
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-10)
[0142]
其中q
v
为水汽密度。
[0143]
式(2-10)为连续方程，实际的计算过程中只能采用离散的网格计算，即使用离散差分方程去逼近连续方程的解。离散差分方程与连续方程在结果上的偏差，以及计算力要求的较为粗糙的网格划分，决定了计算结果在小尺度上的细节难以体现，比如流体的旋转效果等。为了解决这一问题，引入涡流抑制力f
vc
，使用基于速度场的一系列计算，将这些流体细节表现出来。涡流抑制力可用如下方式进行计算，首先计算流体速度场的旋度ω：
[0144][0145]
再计算旋度场ω的梯度方向并单位化：
[0146][0147]
最后计算得出涡流抑制力：
[0148]
f
vc
＝ηl(n
×
ω)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-13)
[0149]
式(2-13)中l为网格单位长度，η为涡流抑制系数，决定涡流抑制力的强度，即决定积云运动过程中的小尺度作用补偿效果，一般取值在0.1到10之间，视具体情况而定。
[0150]
(4)压强场的解算
[0151]
计算力的修正后，我们得到了一个中间速度场，并且在忽略密度场梯度的影响后该速度场满足方程：
[0152][0153]
上式右侧已知，因此是一个泊松方程，运用迭代法可以解出压强p，运用压强对速度进行修正，可以得到最终的速度场：
[0154][0155]
至此，计算出流体在该单位时间下的运动特征。
[0156]
(5)粒子的相变特征
[0157]
云的相变是指云中出现的水的状态变化，即在液滴、水汽、冰晶之间的转化过程。在本发明实施例中，假设相变的发生是以云中体素为单位的，即每一个体素内，水分子状态变化的过程。通常，水汽密度超过当前条件下的饱和点，水汽会转化为液滴，进而形成之后看到的云。空气中水汽密度的饱和点与当前位置的温度t、大气压强p有关，可用如下方法进
行计算：
[0158][0159]
式(2-16)中压强需要考虑当前位置的压强修正：
[0160][0161]
式(2-17)中z为当前位置的海拔高度，p0为标准大气压，可以取100kpa，t0为地表温度，可以取280k-310k，γ取10k/km，表示每上升1km，气温下降10k。
[0162]
对于云体，必须要考虑在相变过程中水量的守恒，在积云中，满足连续条件如下：
[0163][0164]
式(2-18)中q
v
为水汽密度，q
c
为液滴密度。c为常数，即发生相变的水分子总量是不变的。
[0165]
当水汽密度大于饱和点时，发生相变，水汽密度和液滴密度发生变化：
[0166]
δq
v
＝-c＝min(q
vs-q
vb
,q
cb
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-19)
[0167]
q
v
＝q
vb
+δq
v
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-20)
[0168]
q
c
＝q
cb-δq
v
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2-21)
[0169]
式(2-19)中q
vb
、q
cb
是指前一状态下体素内的水汽密度和液滴密度。
[0170]
虽然绝热近似是水蒸气未饱和的空气运动的有效近似，但不能假定饱和空气的位势温度是恒定的。如果体素内水汽密度超过饱和点，水蒸气就会凝结并释放潜热，使体素气团升温。如果由于冷凝和蒸发而产生的潜热和冷却是唯一的非绝热变化，则热力学第一定律成立，有下式：
[0171][0172]
式(2-22)中l为水汽化潜热，在0摄氏度取2.501j/kg。上式说明位势温度的变化除了随速度场的输送变化，还与当前位置的相变有关。所以在考虑输送变化后，有下式：
[0173][0174]
式(2-23)中θ
b
表示位势温度随速度场输送之后的中间位势温度。至此，解算出粒子相变的变化特征。
[0175]
计算渲染一个单位时间下(在仿真过程中表现为一帧)的积云的密度场，除了需要计算每个体素的运动、相变、热能转换等，随后还需要对偏微分方程进行迭代解算以求出压强，并用此压强分布对速度场进行修正。如此复杂的计算需要相当大的计算次数，但是仿真的实时性的要求对每一帧的计算时间有限制。为了能够实现实时仿真(帧频至少达到30)，本发明实施例将所有的仿真计算全部转移到gpu中进行，利用gpu的高并行性对重复的计算进行优化，从而在保证高真实性的同时，提高帧频以达到实时仿真要求。积云内部各体素在一帧时间下的计算过程如图3所示，图6展示了实施例二实现s2的过程，在初始时，初始积云如图7，白色区域表示初始积云，方框表示一个网格，需要对每个网格的粒子进行扩散，才能
得到扩散后的积云。
[0176]
下面对实施例一至实施例十涉及到的太阳光线照射积云发射的过程内容进行详述。
[0177]
如图8所示，为积云的光照示意图，太阳光与视线方向不相同，太阳光经过散射进入视点。针对云中的一个点p，考虑点p进入视点的光线。首先，光照透过部分云层进入p点，经由p点散射后沿视线方向出射一次散射光线；其次，光照透过部分云层进入其它点，例如p1、p2，光照经过二次、三次散射后进入视线；最后，云层表面的漫反射以及自身辐射等进入视线。另外，对于不同的点，视线的终点不同，因此还需要考虑天空背景辐射，例如图8中需要考虑沿视线方向透过的天空辐射。当辐射仿真以天空为背景时，可以将天空辐射用背景辐射的计量算法来处理。当视线方向与光照方向相同时，需要考虑沿视线方向透过的太阳辐射。
[0178]
(a)积云中的散射现象
[0179]
(1)散射计算方法
[0180]
计算云中的散射效果，首先需要计算散射相函数。散射相函数描述的是入射能量经粒子散射后的能量角分布。设入射方向与散射方向之间的夹角θ为散射角,引入一个无量纲参数来描述散射能量随散射角的变化，称之为相函数p(θ)。对于相函数，有归一化条件：
[0181][0182]
米氏散射相函数(下称mie函数)可以通过求解麦克斯韦方程得到，但是过程太过复杂，在此不再赘述。最终球形粒子在单波长下的散射相函数p可以表示为：
[0183]
p＝2π(|s1|2+|s2|2)/(k2csca)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3-2)
[0184]
其中，k为波数，csca＝∫
4π
x|/k2dω是散射截面，x为矢量散射振幅，s1和s2为复振幅，可以通过求解贝塞尔方程得到，本发明实施例中不作详细说明。
[0185]
低空的积云是由不同大小的水滴混合而成的，假设水滴为球形，此时水滴的尺寸对云的光学特征有巨大的影响：对于一定数量的水，水滴越小，云越不透明。大约尺寸缩小一半，消光效率就会加一倍，mie散射相函数也随着液滴的尺寸而变化。要分析这一情况，首先需要对液滴尺寸进行建模。积云中液滴尺寸分布(dsd)模型近似于修正的伽马分布。定义为：
[0186][0187]
其中n0＝300cm-3
为水滴中总体粒子数，γ为尺寸分布的宽度，通常取γ＝2。γ(γ)为伽马分布函数，r
n
为尺寸分布的特征半径。r
e
＝(γ+2)r
n
为这个分布对应的有效半径，本发明实施例取r
e
＝7μm。确定参数后，可以得到液滴尺寸分布。
[0188]
对应于液滴尺寸分布的有效相函数p(θ)可以通过对不同液滴尺寸的相位函数进行加权和求和得到，求和公式如下：
[0189]
[0190]
最终可得指定波长下对应的有效相函数。
[0191]
由于有效相函数整体计算流程颇为复杂，很难将其嵌入到实时计算当中，因此大多采用的是预计算方法，将其保存在纹理等工具内供实时计算直接取值使用，或者可以通过高阶函数模拟简化计算，以保证仿真过程的实时性。
[0192]
确定液滴尺寸模型和相函数后，可以考虑具体的散射过程。假设视线方向为竖直向下，太阳光方向竖直向下透射，透射是指光线未经过粒子散射直接透过云层最终进入视点的过程。忽略云层的自发辐射时，在这个过程中只需要考虑云层对光线衰减作用，不用考虑粒子的散射作用，最终可表达为：
[0193][0194]
式中τ(d)为光线经过云层时的透过率，当视线方向与光线方向不同时不需要考虑透射作用。
[0195]
假设积云内部单个体素为均匀分布，此时体素的透过性可以用其液滴的有效半径r
e
和液滴的粒子数密度n来描述。利用这两个参数可以定义消光截面：
[0196][0197]
其中q为消光效率因子，是波长的函数。由消光截面可得消光系数：
[0198]
κ＝nσ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3-7)
[0199]
由此可得云的透过率：
[0200][0201]
式中x为指定的光线路径，可见透过率τ是在该路径上的积分结果。考虑路径通过的所有体素，可得到按体素累加的透射结果，其中l代表体素网格对应的长度。从频率角度考虑，这是一束光线经过一段路径后的衰减结果；从概率角度考虑，这是一条光线经过一段路径不碰到粒子的概率，也就是经过一段路径不发生散射的概率。因此在这段路径x内命中粒子发生散射的概率为1-τ，对于其导函数：
[0202]
s(x)＝κe-∫κdx
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3-9)
[0203]
可以理解为在路径x上的某一位置命中粒子的概率密度，因此s(x)是恰好在长度为x的自由路径的末端命中的概率密度。
[0204]
(2)一次散射过程
[0205]
一次散射是指光线只经过一个粒子散射随后透过云层最终进入视点的过程。忽略云层的自发辐射时，在这个过程中需要考虑云层对光线衰减作用，以及粒子的一次散射作用，最终可表达为：
[0206]
s1＝s(d
l
)p(θ
vl
)τ(d
v
)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3-10)
[0207]
等号右边第一项为光线透过云层直接命中指定粒子(例如粒子p)的概率；第二项为相函数，表示的是在当前角度下粒子沿指定方向(例如视线方向v)散射的概率；第三项为光线离开粒子，在未发生碰撞的情况下直接进入视点的概率。
[0208]
(3)二次散射过程
[0209]
二次散射是指光线经过两个粒子散射后透过云层最终进入视点的过程。忽略云层的自发辐射时，在这个过程中需要考虑云层对光线衰减作用，以及粒子的一次、二次散射作
用，最终可表达为：
[0210][0211]
对于二次散射过程，一条确定的散射路径包含两个散射粒子：第一次散射粒子以及第二次散射粒子。光线透过云层命中第一次散射粒子，之后命中第二次散射粒子，最终透射出云层最终进入视点。和处理一次散射过程相同，考虑散射出射粒子，也就是第二次散射粒子。散射出射粒子会散射所有命中该粒子的一次散射光线，例如对于出射点p而言，p1与p2以及其他粒子都会作为p的一次散射粒子，对其提供一次散射出射光线。
[0212]
事实上，积云中所有粒子都可以作为一次散射粒子，区别只有一次散射光线的角度和亮度。但在实时仿真过程中，由于计算成本过高，无法使每一个出射粒子的散射都考虑所有可能的一次散射粒子，因此需要采用一些近似算法来进行化简。
[0213]
如式(3-11)右一项，对于不同角度不同距离的光线进行积分，由于积分太过复杂，因此需要对其进行简化计算，右二、三式是对其进行的近似计算。其中，一式到二式基于同一散射角度下不同路径的辐射权重相同假设，二式到三式则是将二次散射粒子的消光系数代替一次散射粒子的消光系数以简化计算。第一个简化计算主要建立在路径差越大，透过率越低，权重越小的基础上；第二个假设计算主要建立在云中液滴距离较近时粒子数分布相差不大的基础上。
[0214]
在近似积分后，结果与一次散射点的位置无关，最终的计算结果只与二次散射点的位置有关，与一次散射计算类似。因此，可以在计算一次散射的同时，计算该散射点作为二次散射点的散射效果，从而避免多次计算，减轻计算压力。
[0215]
(4)多次散射计算
[0216]
多次散射是指光线经过三个或三个以上粒子散射后透过云层最终进入视点的过程。忽略云层的自发辐射时，在这个过程中需要考虑云层对光线衰减作用，以及粒子的三次及以上的散射作用。
[0217]
可以采用类似与二次散射的方法来近似计算云朵中的多次散射，但无论是计算的效率还是计算得到结果的精确度都并不令人满意。这主要是因为以下原因：1.多次散射作用理论上包含3次到无限次散射过程的计算，即便对于一个确定散射出射点，也无法逐一计算散射作用，最终依然要采取近似算法忽略或者简化一部分计算。即便这样，得到的结果也受近似算法的影响。2.二次散射的近似计算始终包含着方向性、对称性的特点，方向性是指二次散射近似结果与一次散射近似，在入射方向的小角度上，包含出射的大部分能量；对称性是指二次散射近似结果是对称的，即正负散射角度下的散射结果相同。对于多次散射，由于光线的多次散射使得散射的方向性明显下降，随着散射次数的增多，最终散射的方向性完全消失，在不同角度上的散射能量分布相同。因此，如果采用二次散射的近似算法来处理多次散射，会使计算结果也会包含很强的方向性，与实际结果相差较大。
[0218]
对于多次散射计算结果，采用如下假设进行近似处理：1.忽略散射次数对散射结
果的影响，统一考虑三次散射级以上散射效果。2.忽略入射光方向对多次散射效果的影响，假设多次散射效果是均匀的。3.忽略散射效果的方向性，以漫反射模型近似处理。
[0219]
当积云处于热平衡状态时，基于假设1以及能量守恒定律，可以得到如下公式：
[0220][0221]
式中代表光线的透过率，代表总体的一次散射、二次散射比例，代表多次散射比例，有
[0222][0223][0224][0225]
因此可得总体的多次散射比例，有
[0226][0227]
基于假设2可以得到单个面元a的辐射比例，有
[0228][0229]
基于假设3可以得到单个面源上不同方向的辐射比例，有
[0230][0231]
其中θ为观察方向与出射位置面片法向量的夹角。
[0232]
当忽略积云形状的影响，以最简单的均匀球体模型进行近似计算，可得到一个较为快速的计算方法。这种方法中的变量仅有入射光的方向、云的尺寸以及云的平均粒子数。在一个仿真场景中，太阳光的方向一般是固定的，而要仿真的云在预设形状确定时，最终的稳定形状也是固定的，因此可以对这三个变量进行预设以及预计算，从而加快仿真效率，达到实时要求。
[0233]
(5)天空辐射
[0234]
由于天空辐射效果远小于太阳光光照效果，因此不考虑天空辐射进入积云之后的散射效果，将其当作背景辐射处理。如图9所示，观察方向与天空的交点所在的小面片作为对应视线方向下天空辐射的主体发出辐射能量，经由积云的透过进入视点，最终影响观察效果。因此，天空辐射的能量与小面片所处的位置、温度、时段、能量经过大气路径下的衰减、能量经过积云的衰减等有关，在忽略大气影响，并假设仿真时段内温度与太阳位置保持不变的条件下，天空辐射只与小面片的位置、能量经过积云的衰减有关。小面片的位置，尤其是小面片位置对应视点的天顶角决定了天空辐射的初始能量，采用预处理的方法，将位置对应的能量计算完成后保存在纹理中，在计算辐射能量时直接在纹理上采样获得；经过积云的衰减可以使用积云的透过率来表示，这与积云的形状、密度分布有关，可以在之后的计算中获得。
[0235]
下面需要对在红外场景中的自辐射过程进行进一步详述。
[0236]
(a)红外光照模型
[0237]
红外场景中与可见光场景中的积云在仿真上存在许多差异，其中最重要的差异在于对云的自发辐射的处理，其次是因为波段不同导致的散射计算问题。
[0238]
(1)散射计算修正
[0239]
红外波段的散射计算需要在可见光波段的计算基础上进行两个修正：1.由于波段不同导致的散射效率因子变化；2.红外波段通需要考虑波段宽度，并以此修正散射相函数。
[0240]
对于第一点，只需在进行式计算时，代入波段内效率因子平均值即可，在波段内效率因子差距不大时，可以使用中间波长对应的效率因子作为平均值使用。在4μm波长下，液滴的折射率为1.351，粒子有效半径为7μm时，可以得到消光效率因子大约为2.6。
[0241]
对于第二点，则需要计算对应波段内每个波长对应的散射相函数。在根据式(3-4)计算出单波长下的散射相函数后，计算整个波段对应的平均相函数。
[0242][0243]
式(4-1)可以得到在对应波段下有效相函数的精确结果，但是较为复杂，事实上，在一个较短的波段内(比如中红外波段3-5μm)，有效相函数的变化可以忽略，因此可以使用波段内中间波长对应的有效相函数作为平均值来计算整个波段内的有效相函数。
[0244]
(2)自发辐射修正
[0245]
在可见光场景中，虚拟场景的渲染是以rgb三色的运算为基础的，也就是在可见光的红绿蓝三个波段下的辐射亮度为计算对象，但根据普朗克辐射定律，云在这三个波段内的自发辐射能量极小，相比太阳光的光照辐射以及光照导致的漫反射效果来说，完全可以忽略。相反，在红外波段下，尤其是中红外波段(3-5μm)与远红外波段(8-12μm)下积云的自发辐射能量较为集中，相比于光照引起的散射效果不能忽略，因此需要考虑在红外场景仿真中云的自发辐射效果。
[0246]
从仿真的真实性和实时性两个方面考虑，处理红外场景仿真中云的自发辐射效果的传统方法通常是采用以下两种计算模型：黑体近似模型和灰体近似模型。
[0247]
黑体近似模型用于处理远红外波段(8-12μm)下的云红外仿真问题，在这个波段下，太阳辐射强度较小，光照效果对云的仿真影响较小，此时云的红外辐射以自发辐射为主。因此在只考虑太阳为唯一光源，天空作为仿真背景时，可以将云视为黑体，此时的自发辐射能量可以通过式(4-2)普朗克公式直接计算。
[0248][0249]
其中，m
λbb
为光谱辐射出射度，单位w/(m2·
μm)，λ为波长，单位μm，t为热力学温度，单位k，h为普朗克常数，c为光速，c1为第一辐射常数，c2为第二辐射常数，c2＝hc/k
b
≈1.43879
×
104μm
·
k。
[0250]
灰体近似模型用于处理中红外波段(3-5μm)下的云红外仿真问题，在这个波段下，太阳辐射占主导地位，云的辐射效果以散射太阳辐射为主，自发辐射为次。因此在只考虑太阳为唯一光源，天空作为仿真背景时，可以将云视为灰体，辐射发射率一般在0.7-0.9之间。
[0251]
这两种模型都是基于实测数据进行的近似假设，尽管比较简单，但是仿真结果贴
近真实红外图像。参考这两种简单近似模型，对第三章提出的光照模型进行自发辐射修正，以保证光照模型的准确。首先，无论是中红外波段或是远红外波段，需要考虑积云对红外辐射的吸收，假设吸收率为α，因此根据能量守恒定律，式(3-12)应当变为：
[0252][0253]
对于黑体和灰体的吸收率，由基尔霍夫定律可知，
[0254]
α＝ε
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4-4)
[0255]
即吸收率等于发射率，特别地，黑体的发射率等于1。
[0256]
基于灰体近似模型，在中红外波段吸收率为0.7-0.9，代入式(4-3)可以计算得到当前的多次散射能量比例，根据计算结果多次散射比例不足10％。并且本发明实施例的多次散射使用的是类似于漫反射的处理方法，因此完全可以将多次散射效果忽略，将这部分能量使用自发辐射方法进行处理，因此有
[0257][0258]
在远红外波段，由于消光效率因子更小，吸收更强，多次散射效果完全可以忽略。在物理上表现为云中吸收较强导致云中太阳辐射散射次数降低、散射强度变弱，以至于难以探测。另外，对于远红外波段下的一次散射和二次散射效果，由于相应波长下相函数、散射效率因子以及太阳辐射更小，在误差允许的范围内相比自发辐射同样可以忽略，因此在远红外波段，积云展现出黑体辐射特征。
[0259]
给出自发辐射的计算方法后，需要确定相关的仿真方案，将其包含在之前的仿真方法中。类似于光在粒子上发生散射后的透过效果，自发辐射从粒子出发，按照观察路径进行传播，经过路径的吸收和散射，最终进入视点，产生观察结果。如图9所示，p2处粒子产生自发辐射，辐射效果沿p2a进入视点，同样p1、p处粒子产生自发辐射，辐射效果沿p1a、pa进入视点，区别在于路径长度不同导致的透过率不同。由自发辐射的这个特点可以确定仿真方案，与光照模型类似，自发辐射特性放在三维模型确定后的渲染着色器中进行，在由观察路径进行积分时，累加自发辐射积分结果，最终输出该路径下的总自发辐射量。
[0260]
如图2，在计算透过、散射、自发辐射能量后，得到该路径下的总能量。与可见光仿真直接用rgb三色进行计算不同，红外仿真需要一定的量化尺度。在确定的尺度内，对辐射输出的结果进行量化，输出最终的量化结果，可以定量的分析仿真实验。假定一个物体的辐射亮度输出为l
output
，仿真场景量化上限l
max
，下限l
min
，可得最终的量化结果g，有
[0261][0262]
式(4-6)的得到的量化结果g一般由灰度表示，在常见的0-255灰度场景中，可得辐射亮度输出为l
output
对应的灰度g
·
255。
[0263]
在进行展示之前，首先确定仿真场景的量化尺度与温度设置，场景量化尺度为辐射亮度尺度，仿真波段为中红外波段(3-5μm)，设置下限为270k对应的黑体辐射亮度，即0.6064w/(sr
·
m2)，上限为330k对应的黑体辐射亮度，即5.3518w/(sr
·
m2)。地面场景温度为300k，积云底部高度1000m，因此积云网格底部温度290k，在网格中心设置水汽生成区域。设置太阳光光照方向为(1,-1,0)，观察方向(-1,0,0)，太阳表面温度取5500k，当做黑体近似处理。
[0264]
在红外光照仿真中应当存在三种光照效果，分别是一次散射、二次散射以及自发辐射效果。与可见光光照仿真不同，由于积云在红外波段的自发辐射不可忽略，所以集中在一次散射、二次散射和多次散射中的辐射能量会因为自发辐射的存在而变弱。散射衰减的能量被粒子吸收，并通过自发辐射形式进行输出。结合第一节中提到的仿真方案进行场景渲染，结果分强方向性光照(一次散射、二次散射)以及自发辐射两个方面进行展示。
[0265]
图10展示的是在红外场景中积云的红外光照仿真结果，从上到下依次为仿真光照的综合效果、强方向性光照(一次散射、二次散射)以及自发辐射效果。从图中可以很明显的看出，相比于可见光仿真中太阳光所占的主导地位，在红外仿真场景中，强方向性光照(一次散射、二次散射)所占的比例明显下降，但仍然是光照能量的主体，并且由于太阳光光照方向的原因，在当前仿真视角内产生的光照效果方向性极强；仿真结果中部分能量来自自发辐射，并且自发辐射会因为积云三维模型的变化而变化，但是自发辐射的能量同样会被散射和吸收，这导致自发辐射结果随时间变化，即随着积云的发散产生明显的变化。
[0266]
分析稳定状态的仿真结果，其中最高灰度255，出现在积云上边缘，主要是因为发生了强前向散射，散射光在小角度内包含了大量的太阳能量，因此辐射亮度很大；自发辐射效果中平均灰度35，对应辐射亮度1.3508w/(sr
·
m2)，对应黑体温度287k，略低于云体所处高度的大气温度，符合积云发射率的实测结果。
[0267]
由图10可见，自发辐射效果没有多次散射效果的均匀性，反而是因为积云三维模型的变化产生明显的变化。这主要是因为在仿真方案中，自发辐射的计算对象是云中的云粒子，因此自发辐射从粒子出发，在传播一段路径后才会进入视点，在路径上依然会产生吸收和散射等效果，这导致仿真结果的有明显的动态性。
[0268]
另外，图中的仿真结果还有一个明显的特点，在仿真过程中，积云中心的部分区域自发辐射效果比其他区域要强。这主要是因为中心区域是水汽供给的位置，所以会导致两个结果：1.中心区域产生相变过程比其他区域剧烈，水汽凝结产生的热量导致温度升高；2.中心区域的浓度相比其他区域较高，因此透过率相对较低。由于这种结果只会对自发辐射产生影响，因此在可见光仿真中不明显，而在红外仿真中较为明显。在仿真过程中，可以通过降低水汽初始浓度和减少水汽供给区域来有效减弱这种现象，但是会对三维模型的最终仿真结果产生影响。
[0269]
此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，“多个”的含义是两个或两个以上，除非另有明确具体的限定。
[0270]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。此外，本领域的技术人员可以将本说明书中描述的不同实施例或示例进行接合和组合。
[0271]
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在
不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。