CN106570215B - 一种基于时变气象数据的热带气旋动态仿真方法 - Google Patents

Abstract

本发明是一种基于时变气象数据的热带气旋动态仿真方法，属于计算机图像学相关技术。主要包括基于气象数据的热带气旋建模方法，多次前向散射的云绘制算法方法，和基于位置的流体的热带气旋仿真算法方法三部分。首先对热带气旋的时变气象数据进行分析，抽取速度场数据备流体仿真使用，提取云强度相关数据，同时建立云的粒子模型，建模成适合计算机图形学绘制的图元数据；再次是实现对图元数据的逼真展现；最后，依据建立的图形模型，使用基于位置的流体仿真算法，实现对气象数据中热带气旋的仿真展示。

Description

一种基于时变气象数据的热带气旋动态仿真方法

技术领域

本发明涉及一种基于时变气象数据的热带气旋动态仿真方法，属于计算机图形学，特别是云的建模领域。

背景技术

热带气旋(Tropical Cyclone,简称TC)指的是一种快速旋转的风暴天气系统，有一个低气压中心和螺旋结构的云体，同时伴随着强风、雷电、暴雨等现象，也被称作台风、飓风、热带风暴、气旋风暴、热带低压等。一个成熟的热带气旋通常有完整的结构，主要包括风眼，地面低压，暖心，中心密集云层区，风眼墙，螺旋雨带，外散环流等构成。

热带气旋的模拟和仿真技术一直是气象、海洋等领域演技研究的重点和难点，而在计算机图形图像学领域，由于其独特的结构，热带气旋也是作为一个重要的研究内容。其中主要包括研究热带气旋的可视化效果，云的真实感绘制技术，以及热带气旋的模拟技术。同时能够较真实的模拟热带气旋的运动以及真实的展现其结构，对气象研究、海洋研究、气象预测等具有极其重要的意义。

在计算机图形图像领域，对热带气旋方面的研究主要集中在对云的建模和仿真。随着虚拟场景渲染技术的发展，云的模拟表现也成为计算机图形学和虚拟现实技术中的热门课题。云的真实感模拟不仅能有效提高场景逼真度，更能传达丰富的气象信息。云作为一种常见的自然现象，由于其形状千变万化，形成、发展和消散的过程又极其复杂，且具有水汽粒子的半透明特征。能够实现即满足气象学应用，同时实现逼真展现的云场景模拟是一项很艰巨的任务。早期在图形图像学领域云建模大体分为两类，即基于过程的云建模和基于物理的云建模。前者侧重于利用噪声、纹理或者交互式的手段对云进行建模，通常需要经历繁琐的参数调整；后者则通过求解简化的NS方程，模拟云生成的物理过程，这种方法耗时大，且不能有效生成预期形状的云。

基于数据驱动的云建模，在一定程度上能够克服之前建模方法的缺点，且由于其建模过程中参考了真实数据，能够在一定程度上反映数据信息、表达气象属性，因此具有很强的应用意义，逐渐成为研究的热点。通过研究和分析云的相关数据(自然图像、观测数据、数值模拟数据、卫星云图等)，进一步从数据中发掘能够指导云建模的信息，构建逼真的三维云场，同时与真实数据在形状、气象信息、规模上具有一定相关性。相应的成果可以为天气数值预测、军事仿真、卫星气象等领域提供可视化工具和环境。

热带气旋的仿真过程是一个流体仿真的过程。其中传统的流体仿真算法主要分为两类，一类是基于有网格的欧拉方法模拟流体运动，一类是基于无网格的拉格朗日方法模拟流体运动。其中网格方法是计算力学中求解工程问题的主要数值计算方法，在求解流体运动，拍击等设计到特大变形的问题时,有限元网格可能会产生严重地扭曲,这样不仅需要网格重构,而且严重地影响计算精度,对于流体运动、拍击、高速撞击等动态问题,显式时间积分的步长取决于有限元网格的最小尺寸,因而网格的扭曲将使时间积分步长过长,大幅度地增加了计算工作量；对裂纹的动态扩展问题,由于裂纹的扩展方向不能事先确定,因而在计算过程中需要不断地重新划分网格以模拟裂纹的动态扩展过程；有限元近似基于网格,因此必然难于处理与原始网格线不一致的不连续性和大变形；复杂三维结构的有限元网格生成也是极具有挑战性的问题。

对于无网格法的研究可以追溯到20世纪70年代对非规则网格有限差分法的研究,由于当时有限元法的巨大成功,这类方法没有受到重视。Lucy和Gingold&Monaghan等分别提出了光滑质点流体动力学(SPH)方法,并且成功地应用于天体物理领域。Johnson等提出了归一化光滑函数算法,使其通过分片试验,可以正确模拟常应变状态,提高了SPH的精度；Danijela等提出了克服零能量模态的具体方案；Monaghan对SPH法进行了总结；SPH法已经被应用于冲击波模拟、流体动力学、水下爆炸仿真模拟、高速碰撞等材料动态响应的数值模拟等等领域。

基于时变气象数据的热带气旋的仿真技术，相比于传统的图形图像学建模方法，由于使用的是数值模拟数据(Weather Research and Forecasting Model，简称WRF)进行建模，模型本身就具有真实性，同时通过时序数据的矫正,使得仿真的数据气象学的意义，符合气象学用户对数据准确性的要求；同时又在时域内使用简单NS方程模拟云的运动过程，使得模拟速度能够达到实时应用的需求。针对以上欧拉方法的不足，本文选择使用拉格朗日方法。

发明内容

本发明技术解决问题：由于SPH在计算过程中对步长的要求比较苛刻，因此本方法提出了一种融入PBD方法的PBF算法，通过求解一组位置约束，从而强制保证密度不变。此方法可以实现类似SPH的仿真过程，同时保留了几何关系的稳定性，最重要的是保持满足实时应用中对长时间步长的要求。通过引入一个人工压力用来改善粒子的分布，创建表面张力并且降低传统SPH对邻居的需求。最后通过应用涡旋约束作为速度处理来解决能量损失的问题。本方法同事使用热带气旋数据作为输入，通过对热带气旋气象数据的分析，建立热带气旋的绘制和仿真所需要的模型，同时提供一种基于位置的流体方法，模拟热带气旋的动态展示。实验表明，本发明提出的方法能够在实现实时逼真的展现热带气旋的运动过程。

本发明技术解决方案：一种基于时变气象数据的热带气旋动态仿真方法，实现步骤如下：

步骤(1)、分析热带气旋的时变气象数据，验证数据的完整性，并获取热带气旋的时变气象数据中的变量，建立适合步骤(2)的绘制和步骤(3)的仿真热带气旋的粒子模型，实现从时变气象数据到热带气旋粒子模型的建模过程；

步骤(2)、使用GLSL实现的基于多次前向散射光照模型的热带气旋绘制方法绘制步骤(1)建立的热带气旋的粒子模型，实现热带气旋的绘制；

步骤(3)、使用基于位置的流体仿真(Position Based Fluids，简称PBF)方法去模拟热带气旋的运动，计算并更新步骤(1)建立的热带气旋粒子模型的下一时刻的粒子属性，供下一次仿真循环中的绘制使用，实现热带气旋的仿真计算；

步骤(4)、重复步骤(2)和步骤(3)，使用步骤(2)的绘制方法绘制由步骤(3)仿真计算的方法更新的热带气旋粒子模型，直到到达热带气旋的气象数据的下一个时刻，再重复步骤(1)到步骤(4)直到用户终止程序。

进一步的，所述步骤(1)中时变气象数据建模热带气旋粒子模型的方法步骤如下：

步骤(1a)、获取热带气旋的时变气象数据，验证数据的完整性，首先验证维度数据是否包含经度、纬度、高度、时间、字符长度，然后验证时变气象数据中是否存在与维度相关的变量数据(包括经度变量，纬度变量，海拔变量，时间变量)，风场变量数据(U-经向风，V-纬向风，W-垂直风)，温度(T)，压强(P)，湿度(H)，密度(D)，云水比(QCLOUD)，冰水比(QICE)，雨水比(QRAIN)，雪水比(QSNOW)，霰水比(QGRAUP)，重力加速度(G)变量，若存在则可以用于建模；

步骤(1b)、步骤(1a)验证完整性后，提取维度变量数据即经度、纬度、海拔初始化热带气旋粒子模型并计算粒子的坐标信息，根据粒子的坐标信息，计算粒子的绘制半径和光滑核半径，其中绘制半径是粒子在绘制过程中，计算粒子的光照对周边粒子的影响，光滑核半径用于仿真过程中计算邻居粒子间以及计算邻居间粒子间相互作用力的大小；

步骤(1c)、在步骤(1b)建立的热带气旋粒子模型基础上，提取风场数据，插值计算每一个粒子的速度，作为粒子的初速度，用于步骤(3)的仿真计算；

步骤(1d)、在步骤(1b)建立的热带气旋粒子模型基础上，提取获取温度，压强，密度，由理想气态方程，计算粒子的质量及质量的倒数，用于步骤(3)的仿真计算；

步骤(1e)、在步骤(1b)建立的热带气旋粒子模型基础上，提取云水比(QCLOUD)，冰水比(QICE),雨水比(QRAIN)，雪水比(QSNOW)，霰水比(QGRAUP)数据，计算出粒子的消光系数，用于步骤(2)的绘制；

步骤(1f)、在步骤(1b)建立的热带气旋粒子模型基础上，提取粒子的重力加速度数据，用于步骤(3)中的仿真计算。

进一步的，所述步骤(2)中使用GLSL实现的基于多次前向散射光照模型的热带气旋绘制方法绘制步骤(1)建立的热带气旋的粒子模型，实现热带气旋的绘制的步骤如下：

步骤(2a)、使用基于多次前向散射的光照模型模拟阳光照射到云粒子，计算粒子吸收的阳光强度和邻居粒子散射的光强，计算粒子的光照强度；

步骤(2b)、根据步骤(2a)计算的粒子吸收的光强和步骤(1e)计算的消光系数，使用Mie相位函数，计算每一个粒子对显示器的作用，通过启动OpenGL的颜色融合来实现云的绘制。

进一步的，所述步骤(3)中热带气旋的仿真计算的方法的步骤如下：

步骤(3a)、设置空气的密度，光滑核半径，以及时间步长，加载步骤(1)建立的热带气旋的粒子模型，初始化仿真计算状态；

步骤(3b)、由时间迭代步长以及热带气旋的粒子模型中的速度，根据粒子的受力情况，计算每一个粒子的欧拉运动位置即欧拉位置；

步骤(3c)、由步骤(3b)计算的粒子的欧拉位置，查找每一个粒子的邻居粒子，构成一个邻居粒子表；

步骤(3d)、对每一个粒子，使用步骤(3c)中计算的邻居粒子表中的邻居粒子，以及粒子的属性，使用PBD(基于位置的动力学模型)原理，求解每一个粒子由应力而产生的位置偏移量，更新粒子的位置，使仿真过程中每一个粒子满足流体的不可压缩性，即求解密度约束；

步骤(3e)、根据步骤(3d)更新的粒子位置，去计算每一个粒子下一次计算所需要的粒子的速度，计算粒子；

步骤(3f)、应用XSPH粘度机制，将粘滞阻力作用到粒子的速度上，使仿真更加真实。

进一步的，所述步骤(3d)中求解每一个粒子的位置的密度约束的方法的步骤如下：

步骤(3d1)根据每一个粒子的邻居粒子，使用SPH(光滑粒子流体动力学)中密度计算公式计算粒子的密度；

步骤(3d2)、根据步骤(3d1)计算的粒子的密度，使用PBD原理，计算每一个粒子受邻居粒子影响的位移权重；

步骤(3d3)、根据步骤(3d2)计算的位移权重，计算每一个粒子由应力而产生的位移，并更新每一个粒子的位置；

步骤(3d4)、循环迭代步骤(3d1)到步骤(3d2)更新粒子的位置，直到满足密度约束或者达到用户指定的迭代次数。

本发明与现有技术相比的优点在于：

(1)本发明输入是时变气象数据，数据本身是通过数值计算以及站点数据同化的得到的，具有一定的准确性，不需要设置初始化参数，用户操作简单方便，同时满足一定的准确性。

(2)传统的仿真过程是基于网格的欧拉法和基于无网格的拉格朗日方法去模拟流体的运动。其中热带气旋在气象学中使用基于网格的欧拉法去仿真其变化过程，但对于流体运动、拍击、高速撞击等动态问题,时间积分的步长取决于有限元网格的最小尺寸,因而网格的扭曲将使时间积分步长过长,大幅度地增加了计算工作量。但是存在计算量大，不能在小型计算机以及微型计算机上应用，而本发明所使用的输入数据就是通过数值计算模拟的输出数据。基于无网格的拉格朗日方法中最成熟的方法为光滑质点流体动力学(SPH，Smoothed Particle Hydrodynamics)方法，但由于SPH方法每一次仿真都是通过计算每一个质点的加速度，通过加速度计算粒子的速度，在通过速度计算粒子的位置，然后执行碰撞检测等过程，使得当时间步长比较大时候，出现穿越现象，不能正确求解粒子的位置。因此SPH方法受到仿真步长的限制，不能应用到大时间步长的应用中。本发明使用基于位置的流体仿真方法，计算的过程中，首先计算粒子有外力产生的加速度，然后将加速度作用到粒子的运动上计算出粒子的位置，然后迭代矫正粒子的位置，直到粒子的位置满足应力(包括碰撞，摩擦阻力，黏性阻力，NS方程)作用，最后在根据粒子的位置计算流体的速度，为下一次计算使用。因本发明是基于无网格的拉格朗日方法，即不需要考虑网格扭曲后带来的时间积分步长过长的问题，同时具有计算量小的优势。同时本发明与SPH方法相比，也不受时间步长的约束。在本发明中，由于有热带气旋的时变气象数据，同时需要实现实时的仿真过程，因此，本发明很适合此种实时应用的场合。

(3)热带气旋的主要绘制方法是基于参与介质的绘制方法，参与介质的云绘制方法主要使用的是基于参与介质的光照模型。参与介质的光照过程中，以热带气旋的粒子为例，主要接受两部分光，一部分是来自太阳等光源和来自天空，背景，地面反射等产生的光，被称作外部光源，一部分就是来自周边粒子散射过来的光，被称作散射光，同时，热带气旋的粒子自身吸收一部分光，同时向周边散射自己接受的光，而传输过程中也会衰减一部分光强，最终计算得到的就是呈现在用户眼中的光强。在参与介质的光照模型中，常用的是简单光照模型、单散射模型，多散射模型以及其他用于近似估计的多散射模型。其中简单光照模型不计算粒子的散射，而单散射模型中，只考虑粒子的一次散射效果，绘制效果没有考虑多散射模型，绘制的效果不好。而对于多散射模型绘制效果好，但是存在的问题就是计算量大，不能满足实时计算的需求。本发明使用一种近似多散射模型的方法，使用的光照模型，是基于假设：参与介质中粒子的散射大部分都是向前散射。由假设，在计算介质光强的时候，总是只考虑太阳前面粒子对后面粒子的影响，而在绘制过程中，也总是考虑相机前面的粒子被相机后面的粒子的影响，这种方法既能实现较真实的绘制效果，同时实现较快的绘制速率。

附图说明

图1为本发明的基本流程图；

图2为三种参与介质的光照模型图；

图3为本发明基于位置的热带气旋模拟仿真流程图；

图4为本发明绘制的一帧热带气旋绘制效果图。

具体实施方式

下面结合附图与实例，对本发明做进一步详细描述：

本发明实时过程中包括三个主要的步骤：

步骤一：分析时变热带气旋气象数据，建立适合后续绘制和仿真使用的热带气旋模型：

从时变的热带气旋气象数据是基于网格数据，通过提取网格数据，建立基于粒子的粒子模型，这样可以在绘制和仿真过程中使用统一数据模型。通过分析绘制(见步骤三)和仿真(见步骤二)过程，设计的粒子模型包括以下属性：

1、粒子的位置仿真和绘制过程中最基础的变量，用于计算粒子的光滑核中心位置，同时计算粒子的下一次位置；

2、粒子的半径粒子的光滑核半径，用于计算粒子的作用范围，在流体仿真和绘制过程中都需要的重要变量；

3、粒子的速度光滑核的运动速度，在仿真中用于计算下一次粒子的估计位置；

4、粒子的质量光滑核粒子的质量，仿真过程中计算粒子的位移，速度等；

5、粒子的加速度粒子所受外力的加速度，通过受力分析计算求得；

6、粒子的入射光在绘制过程的第一阶段通过模拟计算太阳、环境光和邻居粒子对其照射而产生的光照强度；

7、粒子的出射光经过粒子的消光系数，将入射光衰减以后的光照强度；

8、粒子的初始位置粒子的初始位置，用于重设系统；

9、粒子质量倒数粒子质量的倒数，方便计算；

10、消光系数吸收光经过粒子后光照的衰减率；

11、上一帧的位置上一个时间点粒子的位置；

12、上两帧的位置上两帧时刻粒子的位置；

13、到相机的位置绘制过程中绘制粒子的优先级。

其中建模过程中需要确定的信息包括粒子的位置、速度、半径、质量、质量倒数、消光系数、上一帧和上两帧的位置。其中最难获取的两个参数是粒子的消光系数和粒子的速度，下面主要正对消光系数的计算和粒子的速度计算展开。

1、消光系数的计算

热带气旋主要是由于水汽的循环以及收到科里奥利力作用，行程螺旋体的云体，而对于用户来说，可见的主要是空气中的水汽凝结所形成的云系。因此云是由大量的水汽、水滴或冰晶等小颗粒组成的。而在云的绘制过程中，起决定作用的隐私是这些小颗粒聚集在一起的粒子的消光系数。本发明在对云进行建模的过程中使用粒子系统来表征云的属性，而这些信息中粒子衰减系数值对云的绘制是十分重要，该值描述了入射光射入粒子后的衰减比率，不但能够反映网格点上是否含有云，而且在计算云的光照时也会对光线穿透云的光强产生影响。衰减系数越大，粒子能够吸收的光强越大，出射的光强越小，散射的光强越小；衰减系数越小，粒子能够吸收的光强越小，出射的光强越大，散射的光强越大；而当衰减系数小于一定阈值则该网格点处对穿过的光线近似于没有影响，则可以近似认为该网格点处不含有云。因此通过对WRF数据的参量进行分析，计算出粒子的衰减系数值对下一步的工作至关重要。本发明首先选取WRF数据与云粒子的衰减系数值相关的数据参量作为输入，这些参量含义如表1所示：

依据所选参量，参照以下方程计算每一个网格点的衰减系数值：

其中R_f表示不同参量的平均半径，表示单个参量的消光系数，β_ext(S)表示网格点的衰减系数，ρ_air(S)表示空气的密度，ρ_f(S)表示参量的密度，H_f(S)表示输入数据大小。其中不同参量的平均半径如表1所示。依据式3便可以得到网格点的衰减系数值，从而将包含复杂气象参量的WRF数据转化为仅存储衰减系数值的三维网格，使得之后建模工作的输入数据更加简洁。

表1WRF中相关的气象参量表与云粒子属性关系

2、粒子的速度的计算

由于气象数据中的坐标系不同，在计算粒子的初始风速的时候需要考虑风场数据的插值。其中主要使用到了四个坐标系，依次被称作基准网格(base-grid)、U网格(U-grid)、V网格(V-grid)，W网格(W-grid)。其中U、V、W网格依次只得参数U、V、W三个变量的扩充网格。其维度为X*Y*Z，U、V、W网格依次在经向，纬向，海拔三个方向减少半个格点的距离，同时增加一个网格，而使得其基准网格被包含在其之中，来解决边界风场信息，但是在粒子系统的初始化过程中，每一个粒子参照基准网格生成粒子，因此需要在相应的网格中差值，求得每一个粒子各个方向的风速，在方便的仿真热带气旋。

步骤二：多次前向散射的热带气旋的绘制算法

通过之前的工作，已经将气象数据转化为粒子结构的云数据，绘制时只要对热带气旋的云粒子数据进行相应的调度之后就可以进行绘制，此阶段主要求解光线和粒子的相互作用，同时确定粒子和粒子之间的相互作用，确定粒子入射光照强度以及粒子出射光照强度。然而云作为一种参与介质，其内部的物理光照过程比较复杂，为了能够描述这种光线穿透云的光学现象，

图2展示了不同光线对云粒子的入射光作用以及云粒子出射光对人眼的作用，构成了基本的云光照模型。针对任意一个云中的粒子P_i，其光照强度来源主要包括四种情况：太阳光直接照射、其他粒子散射、天空背景光直接照射、地面反射光直接照射。针对太阳光直接照射的光照强度，当太阳光光线以一定的入射角度进入云层后，每穿透一个粒子，都会产生一定衰减，直至到达当前粒子。针对其他粒子散射的光照强度，是由于光线在穿透不均匀参与介质时，由于光学传播较复杂，发生了光线方向偏离的现象，经过多次的不规则作用，到达了当前粒子。针对天空背景光以及地面反射光的光照强度，其计算过程与太阳光直接照射相似，只是入射角度和入射强度有差别。

通过上述过程，云中的所有粒子所具有的光照强度可以计算得到，而人眼看到的云的形态则可以描述为云中所有的粒子都可以发射光线至人眼。因此沿着人眼的视线方向，将所有该条视线上的粒子光照强度进行叠加，同时考虑叠加过程中的衰减，就可以得到三维云的最终绘制效果。

具体实现云的三维仿真通常包括两个步骤，light和render。light是指确定云体中每点的能量，即计算光源(太阳光)以及周边邻居粒子经过散射照射到云粒子上，使得云粒子具有的光照强度，render过程是指确定云最终屏幕成像的过程，即将每一个云粒子看作是一个光源，投影到屏幕上产生的光强作用。Light和render过程中需要考虑光线的散射情况，如图2所示。

1)不考虑粒子之间光照的影响，被称作简单光照模型，其特点就时计算量少，但问题也就是其绘制效果差，不具有真实请。

2)只考虑周边粒子的一次散射效果的光照模型被称作单次散射模型。相对于简单的光照模型，单次散射模型就存在计算量比较大，同时绘制效果要比单次散射光照模型效果好一些。

3)考了不止一次的粒子散射效果的光照模型被称作多次散射模型。多次散射是在单次散射的基础上使云体内的能量趋向于均匀的过程，事实上，多次散射更加符合实际情况。但是实现的复杂度很高，一般都是采用一些近似的算法。

本方法使用一种近似多次散射模型，即多次前向散射模型。假设光能在传输过程中绝大多数的光是沿着光照的前面传输的，在后面传输的光照比较少，本方法忽略粒子向后传出情况，只考虑向前的传输方法，来近似模拟多次散射模型，因此本方法可以通过深度剥离的方法实现一定的加速计算。综上分析本方法能实现快速模拟多次散射光照模型。

步骤三：基于PBF热带气旋仿真算法

本发明采用基于PBF(Position-basis fluid)的流体仿真算法。在流体仿真领域，强制不可压缩性是一个很重要的特性，同时也是十分耗时的计算过程。虽然最近的很多工作已经提高了效率，但对于实时应用来说还是不切实际的，主要是在仿真的过程作用对时间步长的约束比较苛刻，同时在基于力学系统的动力学模型中，无法解决显示积分的过积分现象，等问题。PBS是Miles和Matthias在2013年提出了一种融入PBD(Position BasedDynamics)方法的迭代求解密度的方法，通过求解一组位置约束，从而强制保证密度不变。此方法可以实现类似SPH(Smoothed Particle Hydrodynamic)的仿真过程，同时保留了几何关系的稳定性，最重要的是保持满足实时应用中对长时间步长的要求。通过引入一个人工压力用来改善粒子的分布，创建表面张力并且降低传统SPH对邻居的需求。最后通过应用涡旋约束作为速度处理来解决能量损失的问题。

不同于传统的速度级别的动力学模拟，直接基于位置的动力学模拟可以直接控制物体的位置，不会由于时间步长过大而导致速度由于线性化而产生的不准确，简单的说基于速度方法是求解出速度层次的约束然后进行积分得到位置，而Position Based Dynamic则是先进行速度积分，得出新的位置然后直接在位置层次求解约束，同时将结果和上一个位置的差作为本步的速度。

PBF的流程图参见图3；

在每一次的循环过程中

首先正对每一个粒子，根据外力作用计算粒子的加速度，跟新粒子的速度，再根据粒子的速度，计算预计的粒子位置。

再由粒子的预计位置计算每一个粒子的邻居粒子，得到一个二维的邻居表，用于后面计算粒子正确的位置使用。

然后在迭代计算求解有应力产生的位移，使得满足流体的不可压缩性，(即ρ_i＝ρ₀，其中ρ₀表示大气压强，ρ_i表示粒子的压强)。转换为PBD仿真过程的约束为：

其中P_i表示粒子i，C_i(P₁,P₂,…P_n)表示粒子i需要满足的约束函数，此处称作密度约束，同时ρ_i使用标准的SPH密度估计方法表示，即

其中i和j表示粒子的编号，m_i表示粒子P_i的质量，p_i表示粒子P_i的位置，W(p_i-p_j,h)表示为光滑核函数，P_i表示核函数的半径。而求解上面的约束，需要使用到基于位置的动力学方法求解，其具体方法如下：

对于粒子P_i，期望求得该粒子的校正位置△P满足约束C(P+△P)＝0，而该等式近似为：

而现实△P与同向，这样就映入一个比例系数λ使得△P满足：

这里对粒子的所有属性，使用SPH方法，即：

根据粒子k是不是邻居粒子而有两种形式，即：

即可求得：

在根据λ求得△P。

在计算△P的过程中执行碰撞检测，并更新每一个粒子的调整位置，递归计算知道满足约束，或者满足迭代测试。

最后在在对所有粒子计算粒子的速度和位置，以及执行涡旋约束和XSPH粘性。

在步骤二的热带气旋的绘制方法和步骤三的热带气旋仿真的循环模拟热带气旋的运动效果图如图4所示，其中第一张原气象数据直接绘制的效果，其他为仿真出来的数据绘制的效果。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (4)

1.一种基于时变气象数据的热带气旋动态仿真方法，其特征在于包括以下步骤：

步骤(1)、分析热带气旋的时变气象数据，验证数据的完整性，并获取热带气旋的时变气象数据中的变量，建立适合步骤(2)的绘制和步骤(3)的仿真热带气旋的粒子模型，实现从时变气象数据到热带气旋粒子模型的建模过程；

步骤(2)、使用GLSL实现的基于多次前向散射光照模型的热带气旋绘制方法绘制步骤(1)建立的热带气旋的粒子模型，实现热带气旋的绘制；

步骤(3)、使用基于位置的流体仿真PBF方法去模拟热带气旋的运动，计算并更新步骤(1)建立的热带气旋粒子模型的下一时刻的粒子属性，供下一次仿真循环中的绘制使用，实现热带气旋的仿真计算；

步骤(4)、重复步骤(2)和步骤(3)，使用步骤(2)的绘制方法绘制由步骤(3)仿真计算的方法更新的热带气旋粒子模型，直到到达热带气旋的气象数据的下一个时刻，再重复步骤(1)到步骤(4)直到用户终止程序；

所述步骤(3)中热带气旋的仿真计算方法的步骤如下：

步骤(3a)、设置空气的密度、光滑核半径及时间步长，加载步骤(1)建立的热带气旋的粒子模型，初始化仿真计算状态；所述光滑核半径用于仿真过程中计算邻居粒子间以及计算邻居间粒子间相互作用力的大小；

步骤(3b)、由时间迭代步长以及热带气旋的粒子模型中的速度，根据粒子的受力情况，计算每一个粒子的欧拉运动位置即欧拉位置；

步骤(3c)、由步骤(3b)计算的粒子的欧拉位置，查找每一个粒子的邻居粒子，构成一个邻居粒子表；

步骤(3d)、对每一个粒子，使用步骤(3c)中计算的邻居粒子表中的邻居粒子，以及粒子的属性，使用基于位置的动力学模型PBD原理，求解每一个粒子由应力而产生的位置偏移量，更新粒子的位置，使仿真过程中每一个粒子满足流体的不可压缩性，即求解密度约束；

步骤(3e)、根据步骤(3d)更新的粒子位置，去计算每一个粒子下一次计算所需要的粒子的速度，计算粒子；

步骤(3f)、应用XSPH粘度机制，将粘滞阻力作用到粒子的速度上，使仿真更加真实。

2.根据权利要求1所述的一种基于时变气象数据的热带气旋动态仿真方法，其特征在于：所述步骤(1)中具体步骤如下：

步骤(1a)、获取热带气旋的时变气象数据，验证数据的完整性，首先验证维度数据是否包含经度、纬度、高度、时间、字符长度，然后验证时变气象数据中是否存在与维度相关的变量数据，风场变量数据，温度T，压强P，湿度H，密度D，云水比QCLOUD，冰水比QICE，雨水比QRAIN，雪水比QSNOW，霰水比QGRAUP，重力加速度G变量，若存在则可以用于建模；所述变量数据包括经度变量，纬度变量，海拔变量，时间变量；所述风场变量数据包括U-经向风，V-纬向风，W-垂直风；

步骤(1b)、步骤(1a)验证完整性后，提取维度变量数据即经度、纬度、海拔初始化热带气旋粒子模型并计算粒子的坐标信息，根据粒子的坐标信息，计算粒子的绘制半径和光滑核半径，其中绘制半径是粒子在绘制过程中，计算粒子的光照对周边粒子的影响，光滑核半径用于仿真过程中计算邻居粒子间以及计算邻居间粒子间相互作用力的大小；

步骤(1c)、在步骤(1b)建立的热带气旋粒子模型基础上，提取风场数据，插值计算每一个粒子的速度，作为粒子的初速度，用于步骤(3)的仿真计算；

步骤(1d)、在步骤(1b)建立的热带气旋粒子模型基础上，提取获取温度，压强，密度，由理想气态方程，计算粒子的质量及质量的倒数，用于步骤(3)的仿真计算；

步骤(1e)、在步骤(1b)建立的热带气旋粒子模型基础上，提取云水比QCLOUD，冰水比QICE,雨水比QRAIN，雪水比QSNOW，霰水比QGRAUP数据，计算出粒子的消光系数，用于步骤(2)的绘制；

步骤(1f)、在步骤(1b)建立的热带气旋粒子模型基础上，提取粒子的重力加速度数据，用于步骤(3)中的仿真计算。

3.根据权利要求2所述的一种基于时变气象数据的热带气旋动态仿真方法，其特征在于：所述步骤(2)中使用GLSL实现的基于多次前向散射光照模型的热带气旋绘制方法绘制步骤(1)建立的热带气旋的粒子模型，实现热带气旋的绘制的步骤如下：

步骤(2a)、使用基于多次前向散射的光照模型模拟阳光照射到云粒子，计算粒子吸收的阳光强度和邻居粒子散射的光强，计算粒子的光照强度；

步骤(2b)、根据步骤(2a)计算的粒子吸收的光强和步骤(1e)计算的消光系数，使用Mie相位函数，计算每一个粒子对显示器的作用，通过启动OpenGL的颜色融合来实现云的绘制。

4.根据权利要求1所述的一种基于时变气象数据的热带气旋动态仿真方法，其特征在于：所述步骤(3d)中求解每一个粒子的位置的密度约束方法的步骤如下：

步骤(3d1)根据每一个粒子的邻居粒子，使用光滑粒子流体动力学SPH中密度计算公式计算粒子的密度；

步骤(3d2)、根据步骤(3d1)计算的粒子的密度，使用PBD原理，计算每一个粒子受邻居粒子影响的位移权重；

步骤(3d3)、根据步骤(3d2)计算的位移权重，计算每一个粒子由应力而产生的位移，并更新每一个粒子的位置；

步骤(3d4)、循环迭代步骤(3d1)到步骤(3d2)更新粒子的位置，直到满足密度约束或者达到用户指定的迭代次数。