基于实景植被空间分布粗糙度的精细风场模拟方法与流程

本发明属于地理信息技术和计算机图形学显示的交叉邻域，具体涉及一种基于实景植被空间分布粗糙度的精细风场模拟方法。

背景技术：

目前国内外研究风场的方法有限，除了传统的插值算法外，国内外相对比较完整的检测风的研究方法主要分为集中的现场测风、风洞试验、线性模拟和计算流体力学模拟等。众多文献表明，基于cfd(计算流体力学)的风场模拟是现今能科学模拟风场的有效方法。风在近地表中运行过程中，主要受地形粗糙度和植被分布的影响。现有的模拟大多假设地形的植被是单一或者模拟边界层中的粗糙度是均一的情况。这种模拟的结果与实际的地表覆盖情况下的真实风场存在一定的系统误差。因此本专利公开了一种在复杂地形下，基于实景植被空间分布粗糙度的精细化风场计算方法。一方面，在林相图、遥感影像及实地测量的基础上，本发明采用二维栅格方式存放植被的在复杂地形上的空间分布情况，考虑了植被分布、粗糙度分布、植被高度精细化信息，所以极大的提高了模拟的精度，另一方面，地表粗糙度是研究风场垂直分布风廓线的重要因素，对于准确描绘风廓线的形态起到重要的作用。

地表粗糙度通常有两种理解，一种是从空气动力学角度，因地表起伏不平或者地物本身几何形状的影响，风廓线上风速为零的位置并不在地表，而在离地表一定高度处，这一高度则被定义为地表粗糙度，也称为空气动力学粗糙度。另一种主要是从地形学角度出发，将地面凹凸不平的程度定义为粗糙度，也称为地表微地形。本专利所要讨论的是地表的空气动力学特征，反应地表对风速的减弱作用。空气动力学粗糙度并非像机械加工上仅指物体表面的粗糙程度，而主要是从多相流体力学上，指出物体表面对流经流体的流型、流态及阻滞力影响的一个综合力学参数。空气动力学意义上的地表粗糙度表征地表与大气的相互作用，反应地表对风速的消减作用以及对风沙活动的影响，已经被广泛应用于表征各种地表类型(如沙地、植被、冰雪面、海洋)的空气动力学性质。

当地表被不同植被覆盖时，近地表的气流就会发生变化，则不论是下垫面还是近地面，流场会产生不同的变化。地表粗糙长度是根据地表气流被抬高的高度来预估计的，因此，地表植被的不同高度也会不同，才抬高后产生的气流也会不同，所以本专利提出一种不同植被下最合适的粗糙度，以此来模拟更精细的风场，更细致的描绘不同植被下的风廓线的形成。

本申请在分析了植被的情况下，通过大量的实验来验证粗糙度，并且对比不同的粗糙度参数，采用cfd分别进行计算模拟，并将这种方法运用到未来的森林火灾蔓延过程，为后来的林火预防工作提供一系列的数据支撑，同样可应用于风场发电、污染源扩散以及影视动画特效等方面。

技术实现要素：

本发明的目的在于提供一种基于实景植被空间分布粗糙度的精细风场模拟方法，该方法考虑到地表粗糙度的情况下，在计算机中的虚拟地理环境中，利用计算流体力学来模拟精细的风场，在不同植被下得到更细致、更准确的风廓线。

为实现上述目的，本发明的技术方案是：一种基于实景植被空间分布粗糙度的精细风场模拟方法，包括如下步骤，

s1、选取地形区域：利用三维地形空间分析算法，提取等高线，对坡度、坡向进行分析，提取三维地形数据；

s2、三维地形建模与网格划分：根据三维地形数据进行三维地形的建模，并采用混合网格方法对三维地形进行网格划分；

s3、地形粗糙度参数选定：根据不同植被设置不同粗糙度，进而选定cfd粗糙度参数；

s4、风场cfd计算模拟：根据选定的cfd粗糙度参数，利用cfd模型进行风场模拟计算；

s5、粗糙度影响分析：根据风场模拟计算结果，判断各种植被的粗糙度对于风场模拟的影响。

在本发明一实施例中，所述步骤s2中，所述混合网格方法采用的是契形和四边形混合网格方法。

在本发明一实施例中，所述步骤s3中，根据不同植被设置不同粗糙度长度，进而获得各植被的整体粗糙长度的具体方式为：首先，将各植被按植被种类数分组，对每一组植被进行粗糙度长度为0～0.08m粗糙度长度的模拟实验，即将粗糙度长度分成四个等级，各等级分别为：0～0.02m、0.02～0.04m、0.04～0.06m、0.06～0.08m进行模拟实验；由于选用的参数是粗糙度高度，因此将粗糙度长度转化为相应的粗糙度高度；而后，采用cfd模型进行风场模拟计算，模拟每一种植被下每一种粗糙度长度的风场，并且与实测结果进行对比，以选取每一种植被最合适的粗糙度长度，作为该植被的整体粗糙长度；最后，选取不同植被区域，设置各自的整体粗糙长度，再根据每种植被区域中风速大小来细分区域，风速大的区域粗糙度长度高于整体粗糙长度，风速小的区域粗糙度长度低于整体粗糙长度，其余部分的粗糙度长度等于整体粗糙长度。

在本发明一实施例中，所述步骤s4的风场cfd计算模拟具体实现如下：

s41、湍流模型的选取

三维地形的粗糙度会使得近地面的风场会产生湍流运动，此处湍流模型选用(rng)k-ε模型，其中k方程和ε方程分别为：

式中

s42、边界条件的设置

由于粗糙度所影响的是地形上的风廓线，因此假设风廓线采用指数分布来表示

u(z)表示地面高度为z的平均风速，u*表示摩擦速度，z0表示粗糙度长度，k为vonkarman常数，k为0.4，d为零平面位移；

入口边界设置：根据入口处的植被高度判断，因为模拟出的粗糙度长度z0与粗糙度高度h有如下线性关系，即：

h＝(2m/s)z0＝α×z0(2)

其中，α是地面粗糙度系数，是根据地面覆盖情况而定；α＝2m/s，m代表计算区域内平均单个粗糙度元素所占的水平面积(m²)，s是地面粗糙度元素迎风面平均横截面面积；

通过上述参数通过代入式(1)即可得到入口处的风速；

出口边界设置：出口设置为自由出流；

壁面条件设置：地表植被的分布不同会对地表风速影响很大，主要是由于植被高度的不同影响粗糙度高度的不同；因此，此处将不同植被区域设置的不同的整体粗糙长度，以及每种植被区域中由于风速大小细分区域的粗糙度长度代入式(2)，得到不同植被区域植被的粗糙度高度；

s43、piso算法的设置：由于三维地形复杂且地面起伏度较大，因此采用瞬态piso算法；

s44、风场的计算：根据s41-s43的设置，完成风场模拟计算。

相较于现有技术，本发明具有以下有益效果：本发明所述cfd(计算流体力学)风场模拟的下垫面粗糙度精细化方法，可以克服现有技术模拟风场精度不高、准确性差等缺陷，大大降低了对计算机软硬件的要求，并且对于各种条件下的地形模型更具有适用性，准确性。

附图说明

图1为实验地形植被分布图。

图2为各类植被二维点阵图

图3为粗糙度等级实验图。

图4为植被高度分布图

图5为本发明的流程示意图。

图6为瞬态piso算法流程图。

图7为森林条件下指数风廓线曲线示意图。

具体实施方式

下面结合附图1-7，对本发明的技术方案进行具体说明。

本发明的一种基于实景植被空间分布粗糙度的精细风场模拟方法，目的是为了考虑到地表粗糙度的情况下，在计算机中的虚拟地理环境中，利用计算流体力学来模拟精细的风场，在不同植被下得到更细致、更准确的风廓线。其中本发明中风场的形态要满足三个要求：准确性、真实感和科学性。

在实际场景中，近地面层，一般是边界层下部10％左右，高度在50～100m左右，风速随高度是否遵从对数规律分布主要由当时的大气层结决定。大气层结主要可分为不稳定、稳定和中性三种，在静力中性条件下，风速随高度呈现对数变化规律。因此首先需要判断当时的大气稳定性，再利用有关关系式进行粗糙度的分析，进一步推出风速廓线。本文利用richardson方法来判断大气稳定度，假设ri满足-0.0807<ri<0.0918，根据其标准认为大气处于静力中性状态。

本发明考虑到粗糙度与下垫面的形态结构有关，和粗糙元存在着一一对应关系，主要是有粗糙元的形态大小、空间密度等几何特性决定。以往用于粗糙度计算的方法有lettau、kondo和yamazawa等，但是这些都是运用于粗糙度较均匀的情况下。但是地形上的粗糙度较为复杂，由于植被分布并不均匀，地形复杂，如果采用传统计算粗糙度的方法，计算结果不准确精度不高。所以采用根据植被不同，分布的地理位置不同来预估计较为准确的粗糙度。

在本发明中采用的植被主要是若干种，分别是灌木、树木、草地和岩石等。提取若干种植被的高度数据，根据相应的高度来估计粗糙度。由于粗糙度所影响的是地形上的风廓线，因此假设风廓线采用指数分布来表示：

u(z)表示地面高度为z的平均风速，u*表示摩擦速度，z0表示粗糙度长度，k为vonkarman常数，这里取0.4；根据模拟出的粗糙度长度z0与粗糙度高度h有如下线性关系，即：

h＝(2m/s)z0＝α×z0

其中，α是地面粗糙度系数，是根据地面覆盖情况而定；α＝2m/s，m代表计算区域内平均单个粗糙度元素所占的水平面积(m²)，s是地面粗糙度元素迎风面平均横截面面积；

复杂地形上的风场并不是规则分布，地形、植被分布、温度和湿度的不同，其风场的形成也就不同。每一种植被的粗糙度并不是固定的，而是根据区域划分的不同，其粗糙度也就不同。本专利的目的就是考虑到这些因素，通过实验分析，定义不同种植被最合适的粗糙度，进行分区域分析，根据计算流体力学的知识对每次实验进行模拟分析。

本发明方法的具体实现过程如下。

本发明的一种基于实景植被空间分布粗糙度的精细风场模拟方法，总体上按流程分为选取地形区域并提取数据、三维地形建模与网格划分、地形粗糙度参数选定、风场cfd计算模拟、粗糙度影响分析；具体如下：

第一步：选取地形区域，如图1所示。利用三维地形空间分析算法，提取等高线，对坡度坡向进行分析，并且提取数据。

第二步：导入三维地形实现地形建模，分析地形的几何特点,基于地形的几何属性,提出一种混合网格模型，对地形进行网格划分。

第三步：根据含有实验林区的植被分布以及植被盖度的林相图、高分辨率遥感影像数据，实现对灌木、树木和草地等进行分类，如图2所示，并以一定的网格精度进行标识保存成二维点阵图(栅格数据)。该数据还包括了经纬度坐标，此外，为了提高分类的准确性，可以随机采样几个位置点并到现场进行分类结果验证。

第四步：根据经验将若干种植被分成若干组，每一组植被进行0～0.08m粗糙度长度的模拟实验，将粗糙度分成四个等级，如图3所示。第一个等级粗糙度长度为0～0.02m，后面三个等级分别为0.02～0.04m、0.04～0.06m、0.06～0.08m进行模拟实验。选用的参数是粗糙度高度，所以要将粗糙度长度转化为相应的高度进行设置。此外，通过现场粗糙度实测也可提高粗糙度分类的准确性。

第五步：采用计算流体力学进行模拟实验，模拟流程如图5所示。模拟在每一种植被下的每一种粗糙度的风场，并且与实测结果进行对比，选取每一种植被最合适的粗糙度。

第六步：选取不同植被区域，设置该植被合适的粗糙长度，分析各种植被下风场的效果，与实测值进行对比，分析哪一种植被对于实现风场的准确度影响最大。

第七步：根据地形植被的分布来划分区域，对每一种植被区域先设置整体的粗糙度，此粗糙长度值是模拟实验结果总结出的最合适的粗糙长度，再根据每种植被区域中风速大小来细分区域，风速大的区域粗糙长度高于整体粗糙长度，风速小的区域粗糙长度低于整体粗糙长度，其余部分的粗糙程度等于整体粗糙长度。

第八步：设置湍流模型，地形的粗糙度会使得近地面的风场会产生湍流运动，因此在这里湍流模型选用(rng)k-ε模型，该模型在形式上类似于标准k-ε模型，但是在计算功能上强于标准k-ε模型，它在ε方程中增加了一个fij附加项，使得在计算速度梯度较大的流场时精度更高，并且考虑了旋转效应，因此对强旋转流动计算精度也得到提高。由于地形比较复杂而且倾斜角度较大，所以采用rng能够更好的处理高应变率以及流线弯曲程度较大的流动。(rng)k-ε模型，得到的k方程和ε方程分别为：

式中

第九步：设置边界条件，进行实际地形下的模拟计算，假设入口处的风廓线呈现指数分布，每种植被网格下的粗糙度进行设定，风廓线示意图如图6所示，图中z0代表表面粗糙长度，d代表零平面位移，h代表粗糙度高度，代入分廓线方程中：

u(z)表示地面高度为z的平均风速，u*表示摩擦速度，z0表示粗糙度长度，k为vonkarman常数，k为0.4，d为零平面位移

入口边界设置：假设摩擦速度为u*＝5m/s，根据入口处的植被高度判断，因为根据模拟出的粗糙度长度z0与粗糙度高度h有如下线性关系，即：

h＝(2m/s)z0＝α×z0(2)

其中，α是地面粗糙度系数，是根据地面覆盖情况而定；α＝2m/s，m代表计算区域内平均单个粗糙度元素所占的水平面积(m²)，s是地面粗糙度元素迎风面平均横截面面积，当这里入口处的表面粗糙长度z0＝0.1h时，则d＝0.75h，带入上述公式(1)，可以得到入口处的风速。

出口边界设置：出口设置为自由出流。

壁面条件设置：地表植被的分布不同会对地表风速影响很大，主要是由于植被高度的不同影响粗糙度高度的不同。通常情况下，地表整体的粗糙度参数的选定都只是在一个范围取值，并不准确，精度也不高。这里我们采用通过实验分析，考虑到植被分布不同的情况下，将分组的每种粗糙度长度带入上式(2)，得到粗糙度高度。

第十步：设置算法，由于地形复杂而且地面起伏度较大，所以这里采用piso算法，大多数模拟风场采用的算法都是针对稳态问题，但是多数工程的实际问题是瞬态问题，或者称为非稳态问题。算法流程如图6所示，瞬态问题因场变量与时间有关，因此计算较为复杂，这里采用瞬态问题的piso算法。piso算法原本是基于子分裂技术的求解瞬态问题的非迭代算法，它的精度依赖于所选取的步长，它是一步预测和二步修正。

第一次预测，根据动量方程的离散后的u动量方程和v动量方程：

ai,ju^*i,j＝∑anbu^*nb+(p^*i-1,j-p^*i,j)ai,j+bi,j

ai,jv^*j,i＝∑anbv^*nb+(p^*i,j-1-p^*i,j)ai,j+bi,j

u^*i,j表示在(i,j)处的速度，v^*j,i表示在(j,i)处的速度，ai,j和ai,j表示速度的系数，将系数ai,j和ai,j都增加了a⁰p＝ρ⁰pδv/δt，源项bi,j和bi,j都增加了a⁰pu⁰p和a⁰pv⁰p。

第一修正步和第二修正步。根据离散后的一次和二次压力修正方程：

ai,jp'i,j＝ai+1p'i+1,j+ai-1p'i-1,j+ai,j+1p'i,j+1+ai,j-1p'i,j-1+b'i,j

ai,jp”i,j＝ai+1p”i+1,j+ai-1p”i-1,j+ai,j+1p”i,j+1+ai,j-1p”i,j-1+b”i,j

p′i,j表示在(i,j)处的第一次压力修正，p″i,j表示在(i,j)处的第二次压力修正，其源项都增加了(ρ⁰p-ρp)δv/δt。在这里采用瞬态计算，因此不需要迭代计算，它的精度取决于时间步长，在预测修正过程中压力修正与动量方程计算所达到的精度分别是3(δt³)和4(δt⁴)的量级。当选择足够小的时间步长时，piso算法可以取得较高的精度。

第十一步：每个风场根据其特点在一定范围内随即偏差0.001，初始化设置风场的速度和方向等属性，进行迭代计算。

以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。