CN102646261A - 一种基于gpu加速的多分辨率大规模森林演替过程仿真方法 - Google Patents

Abstract

一种基于GPU加速的多分辨率大规模森林演替过程仿真方法，包括以下步骤：1)，建立树木数据组织与样地数据组织；根据GPU架构对树木数据进行组织，形成利于GPU计算的数数据组织结构，对于整个样地单元建立一与样地网格大小相同的二维数组结构M；2)将每棵树木的遍历转换成cuda内核中单位线程，计算plot[k]的种子分布情况，Y纬度值代表样地，对线程网格的纬度Y＝k处，所有X纬度上的线程进行对样地K的种子分布影响计算；3)种子分布内核方法，成年树NCI方法GPU的内核方法以及成年树光照GLI内核方法。本发明提供一种正确性良好、提升运行效率的基于GPU加速的多分辨率大规模森林演替过程仿真方法。

Description

一种基于GPU加速的多分辨率大规模森林演替过程仿真方法

技术领域

本发明涉及森林动态模拟仿真领域，尤其是一种大规模森林演替过程仿真方法。 

背景技术

对于森林动态模拟仿真，国内外已经提出了很多生长模型，这些模型中，森林的动态仿真大部分都基于以下几个环节：种子分布阶段，萌芽阶段，生长阶段，死亡阶段。在这些阶段，其中包括光照模型的计算，环境资源计算，成年树生长计算，幼年树生长计算等。传统森林仿真不具有科学真实性，森林场景只由主观想象绘制的树木图像简单叠加而成，未充分体现出森林场景的复杂性与客观真实性。另一方面，虽然有很多基于林学数学模型的计算机模拟森林动态仿真，但很多在森林场景的规模上都未达到很大规模，在达到大规模场森林场景的动态仿真时，其方法运行效率不高，通常需要耗费高代价的时间来获得模拟结果，其未能充分发挥现在PC计算机硬件的潜力，尤其是图形硬件GPU处理器上越来越通用的并行计算能力。 

现有方法对上述这些森林动态生长阶段的加速计算主要基于GPU基础上，对于种子分布阶段，本系列方法发明所借助的模型是2000年Philip T.LePage等人提出的WeiBull Function模型，由于对于某块样地的种子数量涉及整片森林场景的树木，因此其计算量是非常庞大的，现在已有的加速方法是基于在CPU上建立四叉树数据结构实现数据优化，进行计算的方法。 

对于成年树阶段，借助的生长模型是2009年K.David Coates等人提出的温带植物资源竞争模型。其中在该成年树生长模型中，对于计算机的模拟主要耗时的是成年树领域内树木竞争的计算与光照资源竞争的计算，在该模型中，其中成年 树领域树的竞争以NCI指数所表示，光照资源以光照指数(GLI)所表示。对于成年树的光照指数(GLI)计算，现有使用基于光线与树冠求交的方法进行实现，即使在GPU上实现求交，由于场景规模的巨大，其计算量与所需线程量也相当庞大。 

现有多分辨率技术主要用于图形图像中的网格简化，通过构造原始网格模型的多个逼近表示，结合硬件资源的绘制能力和绘制误差选择最优的细节层次(level-of-detail，LOD)进行绘制，并在保证绘制速度的前提下尽可能提高绘制质量。 

随着计算机硬件的不断发展，GPU并行计算能力与通用性的提高，另外由于森林场景中树木个体之间具有相同的生长阶段，生长模型公式具有相同性，设计适用于GPU硬件架构的加速方法对于高效大规模的森林动态仿真具有很高的实用价值，对于高效的管理大片森林资源，预测未来发展趋势具有良好的应用价值。 

发明内容

为了克服已有森林动态模拟仿真方法的正确性较差、运行效率较低的不足，本发明提供一种正确性良好、提升运行效率的基于GPU加速的多分辨率大规模森林演替过程仿真方法。 

本发明解决其技术问题所采用的技术方案是： 

一种基于GPU加速的多分辨率大规模森林演替过程仿真方法，所述仿真方法包括以下步骤： 

1)，建立树木数据组织与样地数据组织； 

根据GPU架构对树木数据进行组织，形成利于GPU计算的数数据组织结构：所有计算的树木数据结构将被设计成数组形式；在多分辨率数据组成中聚类结点第k层与第k+1层的数据表示，对于第k层单位结点数据由m个k-1层结点数据累加而成，对于每个单位结点数据，由s种树种对应单位成年树结点数据构成； 

每棵成年树或者幼年树x的树冠被建模成圆柱体，它们在地表的投影都是一个圆C(x)，设样地单元边长为r，则样地单元网格的面积即r²，对于投影完全包含k个样地单元网格的树木，则其暴露面积即为kr²，为正确计算不同树木发生遮挡时暴露面积，设树木树冠上顶面离地面距离为H(x)，所有样地网格单元标记为M_ij，在理想不想交投影区域，则如果某样地网格M_ij完全包含在树木投影C(x)中，则将M_ij分配给树木x，否者不对其分配；在投影相交区域，对于样地网格M_ij被多个投影所覆盖，该样地网格将会被分配给能够完全包含它的H(x)最大的那棵树木x； 

对于整个样地单元建立一与样地网格大小相同的二维数组结构M，所有元素所记录的是所被分配的成年树ID或者未被分配标记；对于具有n棵树木的森林场景，首先遍历树木，遍历中对树木x投影所包含的所有样地网格单元M(i，j)进行分配，如果该网格单元未分配给任何树木，则将其分配给x；如果该网格单元已经被分配给树木t，则通过x树木ID(x)获取树高H(x)，通过记录数组获取t的树木ID(t)，利用ID(t)获取树高H(t)，比较H(t)与H(x)大小，若H(x)大于H(t)，则将该样地单元网格M(i，j)重新分配给树木x，否则不做任何变化；该方法的最大时间复杂度为O(kn)，k为单棵树木最大可能覆盖的样地网格，在1m*1m的样地网格； 

2)线程分配设计：根据模型公式(1)： 

$R_i=\frac{1}{\mathrm{\η}}\mathrm{STR}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left({\left(\frac{{\mathrm{DBH}}_k}{30}\right)}^{\mathrm{\β}}{e}^{{-\mathrm{Dd}}_{\mathrm{ik}}^{\mathrm{\θ}}}\right)---\left(1\right)$

其中R_i表示样地i的种子数量，η，STR，D，θ，β都是与树种相关的参数。DBH为成年树胸径的大小，d_ik表示第k棵成年树距离样地i的实际物理距离，N为所有具有繁殖能力的成年树总数，通过对树种真实参数的代入可知，随着母树与样地的距离的d_ik不断增大，母树所对样地产生的种子数量R_i成单调递减趋 势； 

将每棵树木的遍历转换成cuda内核中单位线程，计算plot[k]的种子分布情况，Y纬度值代表样地，对线程网格的纬度Y＝k处，所有X纬度上的线程进行对样地K的种子分布影响计算； 

对于样地plot(i，j)的种子分布计算，线程网格Y纬度为i+j*PlotWidth的计算其样地的种子分布情况，其中PlotWidth为样地的宽度，对于计算同一样地种子分布的线程划分，是根据邻域树木与样地的距离的不同，线程所计算的结点层次也不同； 

3)种子分布内核方法具体流程如下： 

a获取线程所在块的Y纬度k，，获取线程在X维度上的值； 

b判断该线程所计算结点与线程所计算目标样地单元的距离d，根据不同距离尺度标准从而进行进行判断是否进行细分结点，第二层结点距离尺度标准为ρ2，如果d大于ρ2，则不进行划分，如果小于ρ2，则进行细分结点； 

c如果结点为叶子结点或不需要细分结点，则直接获取线程所对应母结点聚类数据；通过上一步所获得的母结点与目标样地距离进行种子密度的计算； 

d如果细分结点，则该线程将获取母聚类结点所对应子节聚类点数据进行计算，重复上面步骤b； 

成年树NCI方法GPU的内核方法，根据模型公式(6)： 

$\mathrm{NCI}={\mathrm{\Σ}}_{i=1}^S{\mathrm{\Σ}}_{j=1}^N{\mathrm{\λ}}_i\frac{{\left({\mathrm{DBH}}_{\mathrm{ij}}\right)}^{\mathrm{\α}}}{{\left({\mathrm{dis}\tan \mathrm{ce}}_{\mathrm{ij}}\right)}^{\mathrm{\β}}}---\left(6\right)$

计算公式中，i＝1...S为树种数，j＝1...N领域内大于某一固定DBH的树木；α是由目标树树种决定的NCIAlpha参数值；β是由目标树树种决定的NCIBeta参数值；DBH_jk是领域树j的DBH值；λ_ik是树种j相对于目标树种的NCILambda参数；distanceik是目标树与领域树之间的距离； 

成年树NCI方法GPU的内核方法具体流程如下： 

a获取线程所在块的Y纬度k，获取线程在X维度上的值(即线程所计算多分辨率数据根结点层的索引)； 

b判断该线程所计算结点与线程所计算目标树木的距离d，根据距离判断该结点是否在目标树木的邻域影响圈内； 

c如果在影响圈内，则直接获取最底层叶子结点数据，精确判断叶子结点是否在目标树木的影响圈，如果是，则根据模型计算其NCI值，如果不是，则不计算， 

d如果线程所对应结点超出影响圈，则该线程任务结束， 

成年树光照GLI内核方法具体流程如下： 

a数据初始化：分配线程，将线程划分为k批进行处理，每批m个线程，线程总量为k*m，为成年树与幼年树木数量，动态建立实际分辨率所需样地单元网格二维数组记录结构assign_tree，对其所有元素初始化为-1(未分配任何树木ID)，初始化所有成年树暴露面积A为0，获取整个森林场景中成年树与幼年树DBH，位置坐标，树种等属性数组， 

b获取线程号thread_id与块号block_id，通过thread_id与block_id索引树木索引tree_id，对于树木tree_id为i的树木，获取其DBH属性DBH(i)，树木所在位置L(i)，并通过模型公式由DBH(i)计算其高度H(i)，树冠半径R(i)， 

c通过L(i)获取树木i所在样地单元PID(i)，以所在样地单元坐标PLocation(PID)为中心，遍历(x，y)坐标在x∈(L(PID)-R(i)，L(PID)+R(i))，y∈(L-R(i)，L(PID)+R(i))内的所有样地网格，若该范围内的样地单元网格k距离L(i)最远顶角坐标小于R(i)，则该样地网格完全包含于树木i的投影圆内，进入步骤4，否者不作任何处理， 

d样地单元网格k完全包含于投影圆内，通过记录数组获取网格k所被分配 的成年树ID：assign_tree(k)，若该样地网格未被分配，则ID为-1，此时只需将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²，若该样地网格已被分配，则获取该样地被分配成年树ID的树高H(assign_tree(k))，比较H(assign_tree(k))与H(i)大小，若H(assign_tree(k))大于等于H(i)，则不做任何处理，否者，将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²，并对A(assign_tree(k))减去该样地单元面积r²， 

e继续进入步骤c遍历阴影范围内的样地网格，若遍历完，则进入步骤b进行下一棵成年树，直至结束该方法。 

本发明的有益效果主要表现在：在复杂模型下大规模场景的计算高效。对于上面所阐述的方法，已经对其进行了在基于整个生长流程中的系统实现，其与现有普通的计算机森林模拟仿真相比具有以下特点： 

(1)思路新颖。对于多分辨率技术，其主要适用于图形学中应用，虽然有针对其它应用的多分辨率数据组织，但主要针对CPU上数据的组织，对于利用多分辨率组织数据并进行GPU架构上的并行计算该思路具有良好的创新性。 

(2)运行效率高。突破传统CPU加速方法的特点，充分利用当前流行的PC机上GPU硬件并行计算特点，尽可能的发挥了PC机上有限的硬件计算资源，实现模拟阶段中树木的并行计算，对于种子分布与成年树NCI计算，突破单纯GPU简单将程序并行化特点，实现建立不同分辨率的数据进行线程级多分辨率数据选择计算，大大提高了GPU计算效率。 

(3)具有良好扩展性。对于方法在GPU上进行多分辨率数据计算，虽然该系列方法针对大规模森林场景计算机模拟所设计，通过组织多分辨率数据，减轻GPU所开辟线程量，减轻GPU计算负担，但由于该应用所适用对象具有并行性特点，这对于其它具有并行性特点，并且其数据能够根据某一标度(本方法中使用了距离尺度)进行多分辨率组织的应用具有很好的适用性。 

(4)GLI通过暴露面积拟合突破传统光线求交方法，对于林业生长许多光照模型具有参考价值。本方法虽然只对四种树种暴露面积与其光照指数GLI进行了实验数据的拟合，并确立其现线性关系，但该方法对于林学中许多利用光线求交方法的模型具有很好的参考价值，这使得能够在较小误差情况下具有很高的计算效率。 

(5)方法易于实现，通俗易懂。由于该方法是基于NVIDA所推出的CUDA计算架构，CUDA内核是基于C语言基础编程的一个GPU开发平台，因此，对于大部分具有并行编程经验的C语言程序员都能较快速的理解该系列方法，并对其进行实现。 

附图说明

图1是山杨与冷杉两种树种种子数与距离的函数关系图。 

图2是山杨与冷杉两种树种竞争影响因子(CrowdingEffect)与距离的函数关系图。 

图3是组织聚类数据示意图。 

图4是树冠光照指数(GLI)计算示意图。 

图5是成年树部分数据组织结构的示意图。 

图6是聚类结点层次组织结构的示意图。 

图7是样地单元数据组织结构的示意图。 

图8是种子分布CUDA内核设计的示意图。 

图9是成年树NCI CUDA内核设计的示意图。 

图10是基于多分辨率数据CUDA内核线程基本流程的示意图。 

图11是计算暴露面积CUDA内核线程计算流程的示意图。 

具体实施方式

下面结合附图对本发明作进一步描述。 

参照图1～图11，一种基于GPU加速的多分辨率大规模森林演替过程仿真方法所需步骤： 

1)树木与样地数据组织，该步骤主要介绍样地与树木数据在计算机内存中的组织方式。 

树木数据组织结构 

根据GPU架构对数据进行组织，形成利于GPU计算的数数据组织结构：所有计算的数据结构将被设计成数组形式，其中图5为整个场景该成年树部分数据组织示意图，图6为聚类数据组织示意图。 

图5中，左图为树木对象各属性值(树木DBH，树木树种Species，树木所在三维空间位置Location，树木所在样地单元索引号PlotID等)图6为树木各属性值一数组方式进行存储示意图，所示是在多分辨率数据组成中聚类结点第k层与第k+1层的数据表示，对于第k层单位结点数据由m个(图中为4个)k-1层结点数据累加而成，对于每个单位结点数据，由s种树种(图中为4种)对应单位成年树结点数据构成。 

样地数据组织结构 

在整个森林场景中，由于每棵成年树或者幼年树x的树冠都被建模成圆柱体，因此，它们在地表(xy平面)的投影都是一个圆C(x)。若理想情况下树木间不发生遮挡，则树木各自的投影圆面积即可近似为其暴露面积；但由于树木彼此间发生遮挡，不同树木各自投影圆会发生相交(如图4)，就需要合理分配遮挡部分面积。为计算投影圆的面积，本文借助样地网格单元求得圆的面积，即该投影圆的面积近似为其所包含所有样地网格单元面积之和，设样地单元边长为r，则样地单元网格的面积即r²，对于投影完全包含k个样地单元网格的树木，则其暴露面积即为kr²。为正确计算不同树木发生遮挡时暴露面积，设树木树冠上顶面离地面距离为H(x)，所有样地网格单元标记为M_ij。在理想不想交投影区域，则如果 某样地网格M_ij完全包含在树木投影C(x)中，则将M_ij分配给树木x，否者不对其分配；在投影相交区域，对于样地网格M_ij被多个投影所覆盖，该样地网格将会被分配给能够完全包含它的H(x)最大的那棵树木x。例如树木x₁，x₂，x₃同时完全包含样地网格M_uv，但H(x₂)是最大树高，因此，样地网格M_uv将被分配给树木x₂。为正确计算根据树高分配样地网格单元，可以先将场景中的所有树木按树高H(x)按降序排列，其次按照树高从高到低分配样地单元网格，对于已经分配的样地单元，不准进行第二次分配，因此，样地单元总是会被分配给投影能够所完全覆盖它最高的那棵成年树。该方法对于暴露面积精确性取决于样地单元边长r，样地单元网格越精细，其树木暴露面积越准确，但随着样地单元分辨率的提高，计算时间也会伴随着增大，因此选择合适分辨率计算是相当必要的。 

由于对于树高的排序是非常耗时的，为此，在本方法中对于整个样地单元建立一与样地网格大小相同的二维数组结构M，该结构主要记录所对应样地单元各自所分配情况，所有元素所记录的是所被分配的成年树ID或者未被分配标记。为正确分配样地单元网格，在添加该记录结构后，只需遍历一遍场景中所有树木即可。该方法核心思想为对于具有n棵树木的森林场景，首先遍历树木，遍历中对树木x投影所包含的所有样地网格单元M(i，j)进行分配，如果该网格单元未分配给任何树木，则将其分配给x；如果该网格单元已经被分配给树木t，则通过x树木ID(x)获取树高H(x)，通过记录数组获取t的树木ID(t)，利用ID(t)获取树高H(t)，比较H(t)与H(x)大小，若H(x)大于H(t)，则将该样地单元网格M(i，j)重新分配给树木x，否则不做任何变化。对于x投影所包含的所有样地单元网格做同样处理。该方法的最大时间复杂度为O(kn)(k为单棵树木最大可能覆盖的样地网格，在1m*1m的样地网格，对于本文所述的四种树种：冷杉，云杉，黑松，白杨，k＜10，)。图7为方法所采用的记录结构示意图，其左图为样地单元计算暴露面积所需部分属性值，其中Multi_width，Multi_height分 别代表划分样地划分宽度与高度，PLocation代表样地单元中心点所对应三维坐标，assign_tree即在方法中所需要的二维记录表，用于记录对应样地所被分配的树木ID。右图是对PLocation与assign_tree两个属性的说明。其中assign_tree为一整形二维记录表，其元素值代表该样地网格所被分配的树木ID，若该样地网格未被分配，则元素值为-1。 

2)线程分配设计，该步骤主要完成对GPU内核的线程分配设计，从而达到用竟可能少量的线程高效的完成种子分布的计算：

为计算样地plot[k]的种子分布情况，根据模型公式(1)可知，需要对每棵树木进行遍历计算，为提高计算速度，我们将每棵树木的遍历转换成cuda内核中单位线程，这使得种子分布的计算时间大大的提高了。这种情况下，为计算plot[k]的种子分布情况，Y纬度值代表样地，对线程网格(grid)的纬度Y＝k处，所有X纬度上的线程进行对样地K的种子分布影响计算。但是该方法具有一个缺点：会导致X纬度上的线程量过多，从而导致整体计算线程量过多。 

为解决X纬度上线程量过多，可以在一个线程内计算多棵树木，该方法一定程度上解决了线程量过多的问题，但也导致了单位线程运行时间增加，从而导致整个种子分布的计算时间增加。为从根本上解决此问题，下面简单介绍下改进后cuda内核的设计与具体实现方法(图8左面为样地网格，右边为线程网格，黑色圆圈代表聚类结点，圆大小代表聚类结点的大小，图8右图黑色短条虚线单表单位线程)： 

对于样地plot(i，j)的种子分布计算，线程网格Y纬度为i+j*PlotWidth的计算其样地的种子分布情况，其中PlotWidth为样地的宽度。对于计算同一样地种子分布的线程划分，是根据邻域树木与样地的距离的不同，线程所计算的结点层次也不同。如图8所示，对于样地plot(i，j)，与其距离很近的树木(图8左图最内层方块样地区域，样地单元矩形大小最小所在层)，单位线程的计算基本单元是 一棵树木本身的属性值；对于距离相对较远的树木(图8左图中间层样地单元区域)，其单位线程计算的基本单元就是多块样地的树木形成聚类的结点树木属性值；对于非常远的树木(图8左图最外层样地单元区域，样地单元矩形大小最大所在层)，单位线程计算的基本单元就是更大范围的样地所聚类形成的结点的属性值。 

3)内核执行方法步骤设计，该步骤主要对种子分布GPU内核与成年树生长GPU内核的线程内核的设计做了介绍： 

种子分布内核方法具体流程如下： 

a获取线程所在块的Y纬度k(即该线程单元所计算样地k)，获取线程在X维度上的值(即线程所计算多分辨率数据根结点层的索引)。 

b判断该线程所计算结点与线程所计算目标样地单元的距离d。根据不同距离尺度标准(由第一节的距离与种子分布密度的分析所得)从而进行进行判断是否进行细分结点。如，第二层结点距离尺度标准为ρ2，如果d大于ρ2，则不进行划分，如果小于ρ2，则进行细分结点。 

c如果结点为叶子结点或不需要细分结点，则直接获取线程所对应母结点聚类数据(树种参数，聚类DBH属性值等)，通过上一步所获得的母结点与目标样地距离进行种子密度的计算。 

d如果细分结点，则该线程将获取母聚类结点所对应子节聚类点数据进行计算，重复上面步骤b。 

上面阐述了种子分布方法的GPU内核设计，下面阐述成年树NCI方法GPU的内核设计。为计算成年树tree[k]的NCI值，根据模型公式(6)可知，需要对森林中其他所有树木进行遍历计算，同种子分布一样，为提高计算速度，我们将每棵树木的遍历转换成cuda内核中单位线程，这使得种子分布的计算时间大大的提高了。这种情况下，为计算tree[k]的NCI值，对线程网格(grid)的纬度Y＝k 处，所有X纬度上的线程进行对成年树K的竞争影响计算。其方法和种子分布本质是一样的，因此，可以采用种子分布多分辨率GPU计算的策略进行NCI的计算。 

但成年树的NCI与种子分布计算有一个细小的区别：在第上面介绍的模型中，成年树NCI的计算有一最大竞争影响范围圈(即在该影响范围圈外的树木不与目标成年树发生竞争)。为此，为提高成年树NCI多分辨率计算结果的精确性，对于在目标树木影响圈内的树木，在GPU上都是用最高分辨率(即最底层叶子结点)进行计算，在影响圈范围外的结点，都使用最低分辨率计算(因为在影响圈外，所以最后不对最终结果产生影响)。 

对于成年树K的NCI计算，线程网格Y纬度为K的对该成年树进行计算，对于计算同一纬度(Y＝k)线程划分，是根据邻域树木与目标的距离的不同，线程所计算的结点层次也不同。图9中，左面为样地网格，右边为线程网格，黑色圆圈代表聚类结点，圆大小代表聚类结点的大小，右图黑色短条虚线单表单位线程，其中左图不同虚线类型所绘制圆表示不同树种的竞争影响圈，最外层为单棵山杨竞争影响圈，第二层为黑松竞争影响圈，第三层为冷杉竞争影响圈，最内层为云杉竞争影响圈。由图9可知，根据线程所计算树木树种不同，所对应的影响圈也不同，对于在影响圈内的，单位线程计算的是最高分辨率的聚类结点(图9左图最内层样地单元区域，小矩形单元)，在影响圈范围外的，单位线程计算的是最低分辨率的聚类结点(图9左图外层样地单元区域，大矩形单元)数据。 

成年树NCI内核方法具体流程如下： 

a获取线程所在块的Y纬度k(即该线程单元所计算样地k)，获取线程在X维度上的值(即线程所计算多分辨率数据根结点层的索引)。 

b判断该线程所计算结点与线程所计算目标树木的距离d。根据距离判断该结点是否在目标树木的邻域影响圈内。 

c如果在影响圈内，则直接获取最底层叶子结点数据，精确判断叶子结点是否在目标树木的影响圈，如果是，则根据模型计算其NCI值，如果不是，则不计算。 

d如果线程所对应结点超出影响圈，则该线程任务结束。 

对于种子分布与成年树NCI的内核设计虽有不同之处，但其基本流程具有共同性，其两者方法流程可以按以图10进行。 

对于种子分布与成年树方法设计如上所述。下文将对成年树光照GLI的计算方法设计进行阐述。 

光照指数GLI内核方法具体流程如下： 

a数据初始化。分配线程，将线程划分为k批进行处理，每批m个线程，线程总量为k*m，为成年树与幼年树木数量。动态建立实际分辨率所需样地单元网格二维数组记录结构assign_tree，对其所有元素初始化为-1(未分配任何树木ID)，初始化所有成年树暴露面积A为0，获取整个森林场景中成年树与幼年树DBH，位置坐标，树种等属性数组。 

b获取线程号thread_id与块号block_id，通过thread_id与block_id索引树木索引tree_id。对于树木tree_id为i的树木，获取其DBH属性DBH(i)，树木所在位置L(i)，并通过模型公式由DBH(i)计算其高度H(i)，树冠半径R(i)。 

c通过L(i)获取树木i所在样地单元PID(i)。以所在样地单元坐标PLocation(PID)为中心，遍历(x，y)坐标在x∈(L(PID)-R(i)，L(PID)+R(i))，y∈(L-R(i)，L(PID)+R(i))内的所有样地网格，若该范围内的样地单元网格k距离L(i)最远顶角坐标小于R(i)，则该样地网格完全包含于树木i的投影圆内，进入步骤4，否者不作任何处理。 

d样地单元网格k完全包含于投影圆内，通过记录数组获取网格k所被分配的成年树ID：assign_tree(k)。若该样地网格未被分配，则ID为-1，此时只需将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²。若该样地网格已 被分配，则获取该样地被分配成年树ID的树高H(assign_tree(k))。比较H(assign_tree(k))与H(i)大小，若H(assign_tree(k))大于等于H(i)，则不做任何处理。否者，将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²，并对A(assign_tree(k))减去该样地单元面积r²。 

e继续进入步骤3遍历阴影范围内的样地网格，若遍历完，则进入步骤b进行下一棵成年树，直至结束该方法。 

本实施例中，基于大规模森林场景的GPU加速方法主要包括对两个阶段的加速：种子分布多分辨率的GPU加速方法，成年树生长多分辨率的GPU加速方法。其中成年树生长部分包括领域树竞争(NCI指数的计算)与光照计算(GLI指数的计算)两个部分。另外，本系列方法所需模型实验数据以及参数主要使用山杨，冷杉，黑松，云杉四种树种组成。 

首先介绍下本系列加速方法所使用的模型公式基础与多分辨率数据组成原理。对于种子分布，其模型公式如下所示： 

$R_i=\frac{1}{\mathrm{\η}}\mathrm{STR}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left({\left(\frac{{\mathrm{DBH}}_k}{30}\right)}^{\mathrm{\β}}{e}^{{-\mathrm{Dd}}_{\mathrm{ik}}^{\mathrm{\θ}}}\right)---\left(1\right)$

其中R_i表示样地i的种子数量，η，STR，D，θ，β都是与树种相关的参数。DBH为成年树树干胸径大小，由成年树胸径的大小及树种参数所决定，d_ik表示第k棵成年树距离样地i的实际物理距离，N为所有具有繁殖能力的成年树总数。通过对树种真实参数的代入可知，随着母树与样地的距离的d_ik不断增大，母树所对样地产生的种子数量R_i成单调递减趋势。经过真实参数测试可知，当距离大于60m左右时，母树对目标树木产生的种子数趋于0，图1为山杨与冷杉两种树种种子数与距离的函数关系图。横轴为母树与单元样地的距离(单位米)，纵轴为该母树在样地单元中产生的种子密度。因此，对于多分辨率的计算，可以凭借距离为尺度，对于不同距离范围的样地划分成不同分辨率的类簇。由于当多棵树木形成聚类结点时，聚类结点内的所有树木距离近似成聚类结点与目标样地的距离，此时 对于聚类内的所有树木的距离为d_ik是相等的，因此将d_ik近似为d，d为聚类结点到目标样地的距离。则推出对于同一聚类结点内的成年树 

为定值。 

所以，对于一具有m棵成年树的聚类结点对目标样地i种子分布的计算可以表示为： 

$R_i=(\frac{1}{\mathrm{\η}}*\mathrm{STE}*e^{-D{d}^{\mathrm{\θ}}})\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{DBH}}_k}{30}\right)}^{\mathrm{\β}}---\left(2\right)$

由于对于同一聚类结点 

的值是固定的，从而，对于具有m棵成年树的聚类结点数据可以按以下公式表示： 

${\mathrm{NODE}}_i=\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{DBH}}_k}{30}\right)}^{\mathrm{\β}}---\left(3\right)$

成年树的生长模型公式如下所示： 

Growth＝MaxGrowth×SizeEffect×ShadingEffect×CrowdingEffect      (4) 

其中，Growth：为最终DBH的增长量。DBH是成年树树干胸径大小。MaxGrowth：为一棵成年树在最优环境下(无遮挡，无邻域内树木竞争，正处于生长旺盛期)的最大生长量。SizeEffect：主要是当前植株受自身DBH生长的影响因素，主要由当前树木的DBH所决定。ShadingEffect：主要是当前树木受领域内树木遮挡的影响因素。其中，ShadingEffect＝1-GLI/100，GLI为成年树光照影响值。CrowdingEffect是指领域内的树木竞争的影响因素。在该生长模型的计算中，其中CrowdingEffect与ShadingEffect是影响计算效率的两个主要因子。CrowdingEffect其具体计算公式如下(5)所示： 

$\mathrm{CrowdingEffect}=e^{-C*{\mathrm{DBH}}^{\mathrm{\γ}}*{\mathrm{NCI}}^D}---\left(5\right)$

对于上述公式，C，D，γ都是与树种相关的常数参数，DBH也只与自身树木的胸径及其他树种参数相关，是一具体值，在该因素的计算中，竞争影响因子主要体现在NCI的计算上。其中： 

$\mathrm{NCI}={\mathrm{\Σ}}_{i=1}^S{\mathrm{\Σ}}_{j=1}^N{\mathrm{\λ}}_i\frac{{\left({\mathrm{DBH}}_{\mathrm{ij}}\right)}^{\mathrm{\α}}}{{\left({\mathrm{dis}\tan \mathrm{ce}}_{\mathrm{ij}}\right)}^{\mathrm{\β}}}---\left(6\right)$

计算公式中，i＝1...S为树种数，j＝1...N领域内大于某一固定DBH的树木；α是由目标树树种决定的NCIAlpha参数值；β是由目标树树种决定的NCI Beta参数值；DBH_jk是领域树j的DBH值；λ_ik是树种j相对于目标树种的NCI Lambda参数；distance_ik是目标树与领域树之间的距离。由公式(6)可知，NCI值与距离成反比，与CrowdingEffect(其值大于0小于1)成正比。另外，对于成年树生长，每棵成年树有一领域影响范围圈，超出该范围内的树木将不对成年树生长产生影响。经过真实参数测试可知，当距离大于40m左右时，此时NCI的值趋于0，CrowdingEffect趋于1。图2为山杨与冷杉两种树种竞争影响因子(CrowdingEffect)与距离的函数关系图。其中横轴为目标树木与邻域树木的距离(单位m)，纵轴为CrowdingEffect的值，这意味着在该距离外的领域树木对目标树木的竞争影响因素几乎没有。因此，对于NCI的聚类节点可使用上面种子分布所述相同方法以距离为尺度建立多分辨率的数据聚类结点。 

其中α，β，λ为与树种相关的参数，为常数。同种子分布一样，对于同一聚类结点内的成年树，其到目标树木的距离近似相等的，因此将distance_j近似等于distance，其中distance为聚类结点到目标成年树的距离。从而，对于同一聚类结点内的树木，NCI的计算可以近似为： 

$\mathrm{NCI}=\frac{\mathrm{\λ}}{{\left(\mathrm{dis}\tan \mathrm{ce}\right)}^{\mathrm{\β}}}{\mathrm{\Σ}}_{k=1}^m{\mathrm{DBH}}^k---\left(7\right)$

其中m为一聚类结点内的成年树数量。对于同一聚类结点的值 

为定值，从而对于聚类结点i数据属性，可以由该公式表示：

${\mathrm{NODE}}_i={\mathrm{\Σ}}_{k=1}^m{\mathrm{DBH}}^k---\left(8\right)$

在计算目标成年树的NCI值时，可以按上述结点DBH值取代原来的m棵成年树DBH值。 

综上所述，对于种子分布与成年树NCI的计算，其基本聚类策略可以描述为，对于在一定距离之外的树木群体，为计算目标值，可将其群体聚合成一个值由公式(3)或者(8)组成的数据结点，从而计算目标处具体数据值。其示意图如图3，左图为未形成聚类结点前，右图为形成聚类结点后，其中阴影区域为形成聚类结点样地范围，随着距离的不同，所形成聚类结点所包含的样地单元数量也不同。 

对于成年树的光照资源的计算，其上面已经介绍ShadingEffect＝1-GLI/100，其中 

$\mathrm{GLI}\left(x\right)={\mathrm{\Σ}}_{i=1}^m{\mathrm{\Σ}}_{j=1}^{n}L(i,j)---\left(9\right)$

因此，对于光照因子ShadingEffect的计算主要是对GLI的计算。其中GLI(x)为单棵树木x所接收到的光照指数强度，L(i，j)表示在(i，j)方向处光线初始化平均光强，i为光线的经度方向角，j为光线的纬度方向角。在实际自然界中，所有树木所接收的光线不可能完全很理想，很多情况下会被领域树木所遮挡，这就意味着光线会发生衰减。对于越是矮小的树木，在树木密集处，被遮挡的光照量越大。因此，对于树木光照指数并不仅仅是原始光强的累加。为此，对于光照指数的计算，添加一衰减系数λ，当方向为(i，j)的光线与树木发生相交时，其光照强度为L(i，j)*λ。当该方向光线与k棵树木发生相交时，其光照衰减强度即为λ^k。所以对于到达树冠光照强度衰减情况，按以下公式所定义： 

其中P(x，i，j)代表最终光线L(i，j)衰减强度，λ_s代表树种s的透光率，k代表森林场景中与光线L(i，j)相交的树冠集合总数。因此，对于单位个体树冠最终所接收到的光照强度如公式(11)所示： 

$\mathrm{GLI}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}L(i,j)P(x,i,j)---\left(11\right)$

对于GLI的计算原理，可以用示意图4所表示： 

在图4中，对于树冠所获取的光照指数(GLI)值，对于左边相对较矮小树木，由于受到右边树木的遮挡，从右边照射过来的光线发生衰减(光线用虚线表示)。而图中右边树木的光线，未被任何树木所遮挡(光线用实线表示)，因此，它所接收到的光照强度未发生任何衰减。 

虽然利用该求交的方法能够正确计算成年树的GLI值，但由于求交过程的复杂，对于一个10⁵数量级的成年树，为求得每棵成年树的GLI值，在具有100单元方向的光照情况下，将进行10⁵*10⁵*100次的求交，计算量是相当庞大。为此，本发明通过分析GLI计算特点，利用林学上树木暴露面积的概念通过近似的拟合从而计算GLI值。所谓暴露面积即对于一棵树木树冠，光线能够直射到其表面点集合在地表投影的面积，对于树木树冠上的一点p，如果垂直方向光线L(i，j)从p点出发，未与任何树木发生相交，则称该点为暴露点，单棵树木树冠上所有暴露点的集合投影到地表所在三维空间中xy平面上，该投影区域的面积即该树木的暴露面积。 

对于暴露面积与树木光照指数(GLI)，其都有一共同点，就是最终值都受领域树遮挡的影响。对于高大的树木，其被邻域树木所遮挡的面积相对较少，其矮小的树木，被邻域树木所遮挡的树木相对较大。图4中，左边树木相对于右边树木高度低点，导致其从右面方向部分光线被遮挡，也将导致GLI值相对于右边树木偏小。对于两棵树木的暴露面积，将其树冠平面投射至样地平面。从图4中看出，由于左边树木受到其领域右树的遮挡，暴露面积映射投影两圆发生相交。左边灰色阴影区域为较矮小树木暴露面积区域。由图4可知，右边阴影区域将遮挡左边阴影区域，导致左边区域面积减少，这与树木被遮挡导致到达较矮树木获得光照强度衰减原理相同，由此可推断暴露面积的大小与直接利用光线求交计算光照指数可能存在正比关系。本发明已经经过真实森林场景实验数据对暴露面积与树木GLI正比关系进行验证拟合，其关系如表1所示。 

通过上面拟合直线图可以发现，成年树的GLI值与其对应的暴露面积存在着线性直线关系，并且通过拟合直线理想GLI值与实际实验数据值所偏离误差相对较小。通过上述拟合，将成年树GLI值与暴露面积的线性关系设为： 

GLI＝K*Area+B                                        (12) 

下表为四种树种K与B的拟合值。 

树种	K	B
			冷杉Subalpine Fir	3.753549	41.712700
云杉Interior Spruce	3.725126	42.619499
			黑松Lodgepole Pine	1.835209	48.267399
山杨Trembling Aspen	1.239926	48.423000

表1四种树种线性关系斜率与截距拟合值

上面所阐述的为方法所需的计算公式及基本理论准备，下文将详细介绍针对GPU架构的计算的内核方法。对于种子分布，成年树的邻域竞争指数NCI将使用类似的方法，对于成年树生长光照指数GLI将在上文阐述的暴露面积拟合的线性关系基础上进行多分辨率的加速方法计算。 

Claims (1)

1.一种基于GPU加速的多分辨率大规模森林演替过程仿真方法，其特征在于：

所述仿真方法包括以下步骤：

1)，建立树木数据组织与样地数据组织；

根据GPU架构对树木数据进行组织，形成利于GPU计算的数数据组织结构：所有计算的树木数据结构将被设计成数组形式；在多分辨率数据组成中聚类结点第k层与第k+1层的数据表示，对于第k层单位结点数据由m个k-1层结点数据累加而成，对于每个单位结点数据，由s种树种对应单位成年树结点数据构成；

每棵成年树或者幼年树x的树冠被建模成圆柱体，它们在地表的投影都是一个圆C(x)，设样地单元边长为r，则样地单元网格的面积即r²，对于投影完全包含k个样地单元网格的树木，则其暴露面积即为kr²，为正确计算不同树木发生遮挡时暴露面积，设树木树冠上顶面离地面距离为H(x)，所有样地网格单元标记为M_ij，在理想不想交投影区域，则如果某样地网格M_ij完全包含在树木投影C(x)中，则将M_ij分配给树木x，否者不对其分配；在投影相交区域，对于样地网格M_ij被多个投影所覆盖，该样地网格将会被分配给能够完全包含它的H(x)最大的那棵树木x；

对于整个样地单元建立一与样地网格大小相同的二维数组结构M，所有元素所记录的是所被分配的成年树ID或者未被分配标记；对于具有n棵树木的森林场景，首先遍历树木，遍历中对树木x投影所包含的所有样地网格单元M(i，j)进行分配，如果该网格单元未分配给任何树木，则将其分配给x；如果该网格单元已经被分配给树木t，则通过x树木ID(x)获取树高H(x)，通过记录数组获取t的树木ID(t)，利用ID(t)获取树高H(t)，比较H(t)与H(x)大小，若H(x)大于H(t)，则将该样地单元网格M(i，j)重新分配给树木x，否则不做任何变化；该方法的最大时间复杂度为O(kn)，k为单棵树木最大可能覆盖的样地网格，在1m*1m的样地网格；

2)线程分配设计：根据模型公式(1)：

$R_i=\frac{1}{\mathrm{\η}}\mathrm{STR}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}\left({\left(\frac{{\mathrm{DBH}}_k}{30}\right)}^{\mathrm{\β}}{e}^{{-\mathrm{Dd}}_{\mathrm{ik}}^{\mathrm{\θ}}}\right)---\left(1\right)$

其中，R_i表示样地i的种子数量，μ，STR，D，θ，β都是与树种相关的常数参数，DBH为成年树树干胸径大小，由成年树胸径的大小及树种参数所决定，d_ik表示第k棵成年树距离样地i的实际物理距离，N为所有具有繁殖能力的成年树总数，通过对树种真实参数的代入可知，随着母树与样地的距离的d_ik不断增大，母树所对样地产生的种子数量R_i成单调递减趋势；

将每棵树木的遍历转换成cuda内核中单位线程，计算plot[k]的种子分布情况，Y纬度值代表样地，对线程网格的纬度Y＝k处，所有X纬度上的线程进行对样地K的种子分布影响计算；

对于样地plot(i，j)的种子分布计算，线程网格Y纬度为i+j*PlotWidth的计算其样地的种子分布情况，其中PlotWidth为样地的宽度，对于计算同一样地种子分布的线程划分，是根据邻域树木与样地的距离的不同，线程所计算的结点层次也不同；

3)种子分布内核方法具体流程如下：

a获取线程所在块的Y纬度k，，获取线程在X维度上的值；

b判断该线程所计算结点与线程所计算目标样地单元的距离d，根据不同距离尺度标准从而进行进行判断是否进行细分结点，第二层结点距离尺度标准为ρ2，如果d大于ρ2，则不进行划分，如果小于ρ2，则进行细分结点；

c如果结点为叶子结点或不需要细分结点，则直接获取线程所对应母结点聚类数据；通过上一步所获得的母结点与目标样地距离进行种子密度的计算；

d如果细分结点，则该线程将获取母聚类结点所对应子节聚类点数据进行计算，重复上面步骤b；

成年树NCI方法GPU的内核方法，根据模型公式(6)：

$\mathrm{NCI}={\mathrm{\Σ}}_{i=1}^S{\mathrm{\Σ}}_{j=1}^N{\mathrm{\λ}}_i\frac{{\left({\mathrm{DBH}}_{\mathrm{ij}}\right)}^{\mathrm{\α}}}{{\left({\mathrm{dis}\tan \mathrm{ce}}_{\mathrm{ij}}\right)}^{\mathrm{\β}}}---\left(6\right)$

计算公式中，i＝1...S为树种数，j＝1...N领域内大于某一固定DBH的树木；α是由目标树树种决定的NCIAlpha参数值；β是由目标树树种决定的NCIBeta参数值；DBH_jk是领域树j的DBH值；λ_ik是树种j相对于目标树种的NCILambda参数；distance_ik是目标树与领域树之间的距离；

成年树NCI方法GPU的内核方法具体流程如下：

a获取线程所在块的Y纬度k，获取线程在X维度上的值(即线程所计算多分辨率数据根结点层的索引)；

b判断该线程所计算结点与线程所计算目标树木的距离d，根据距离判断该结点是否在目标树木的邻域影响圈内；

c如果在影响圈内，则直接获取最底层叶子结点数据，精确判断叶子结点是否在目标树木的影响圈，如果是，则根据模型计算其NCI值，如果不是，则不计算，

d如果线程所对应结点超出影响圈，则该线程任务结束，

成年树光照GLI内核方法具体流程如下：

a数据初始化：分配线程，将线程划分为k批进行处理，每批m个线程，线程总量为k*m，为成年树与幼年树木数量，动态建立实际分辨率所需样地单元网格二维数组记录结构assign_tree，对其所有元素初始化为-1(未分配任何树木ID)，初始化所有成年树暴露面积A为0，获取整个森林场景中成年树与幼年树DBH，位置坐标，树种等属性数组，

b获取线程号thread_id与块号block_id，通过thread_id与block_id索引树木索引tree_id，对于树木tree_id为i的树木，获取其DBH属性DBH(i)，树木所在位置L(i)，并通过模型公式由DBH(i)计算其高度H(i)，树冠半径R(i)，

c通过L(i)获取树木i所在样地单元PID(i)，以所在样地单元坐标PLocation(PID)为中心，遍历(x，y)坐标在x∈(L(PID)-R(i)，L(PID)+R(i))，y∈(L-R(i)，L(PID)+R(i))内的所有样地网格，若该范围内的样地单元网格k距离L(i)最远顶角坐标小于R(i)，则该样地网格完全包含于树木i的投影圆内，进入步骤4，否者不作任何处理，

d样地单元网格k完全包含于投影圆内，通过记录数组获取网格k所被分配的成年树ID：assign_tree(k)，若该样地网格未被分配，则ID为-1，此时只需将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²，若该样地网格已被分配，则获取该样地被分配成年树ID的树高H(assign_tree(k))，比较H(assign_tree(k))与H(i)大小，若H(assign_tree(k))大于等于H(i)，则不做任何处理，否者，将assign_tree(k)赋值为i，并对A(i)累加该样地单元面积r²，并对A(assign_tree(k))减去该样地单元面积r²，

e继续进入步骤c遍历阴影范围内的样地网格，若遍历完，则进入步骤b进行下一棵成年树，直至结束该方法。

Images