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

摘要

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

主权项

一种基于GPU加速的多分辨率大规模森林演替过程仿真方法，其特征在于：所述仿真方法包括以下步骤：1）建立树木数据组织与样地数据组织根据GPU架构对树木数据进行组织，形成利于GPU计算的数据组织结构：所有计算的树木数据结构将被设计成数组形式；在多分辨率数据组成中聚类结点第q层与第q+1层的数据表示，对于第q层单位结点数据由m个q‑1层结点数据累加而成，对于每个单位结点数据，由s种树种对应单位成年树结点数据构成；每棵成年树或者幼年树x的树冠被建模成圆柱体，它们在地表的投影都是一个圆C（x），设样地单元边长为r，则样地单元网格的面积即r²，对于投影完全包含k₁个样地单元网格的树木，则其暴露面积即为k₁r²，为正确计算不同树木发生遮挡时暴露面积，设树木树冠上顶面离地面距离为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（f*n），f为单棵树木最大可能覆盖的样地网格，在1m*1m的样地网格；2）线程分配设计：根据模型公式（1）：其中，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纬度s，,获取线程在X纬度上的值；b判断该线程所计算结点与线程所计算目标样地单元的距离d，根据不同距离尺度标准从而进行判断是否进行细分结点，第二层结点距离尺度标准为ρ2，如果d大于ρ2，则不进行划分，如果小于ρ2，则进行细分结点；c如果结点为叶子结点或不需要细分结点，则直接获取线程所对应母结点聚类数据；通过上一步所获得的母结点与目标样地距离进行种子密度的计算；d如果细分结点，则该线程将获取母聚类结点所对应子节聚类点数据进行计算，重复上面步骤b；成年树NCI方法GPU的内核方法，根据模型公式（6）：$\mathrm{NCI}={\Sigma }_{v=1}^S{\Sigma }_{z=1}^N{\lambda }_v\frac{{\left({\mathrm{DBH}}_{\mathrm{vz}}\right)}^{\alpha }}{{\left({\mathrm{dis}\tan \mathrm{ce}}_{\mathrm{vz}}\right)}^{\beta }}---\left(6\right)$计算公式中，v=1...S为树种数，z=1...N邻域内大于某一固定DBH的树木；α是由目标树树种决定的NCI Alpha参数值；β是由目标树树种决定的NCI Beta参数值；DBH_zp是邻域树z的DBH值；λ_vp是树种z相对于目标树种的NCI Lambda参数；distance_vp是目标树与邻域树之间的距离；成年树NCI方法GPU的内核方法具体流程如下：a获取对线程网格的纬度Y=k₂处，获取线程在X纬度上的值，即线程所计算多分辨率数据根结点层的索引；b判断该线程所计算结点与线程所计算目标树木的距离d，根据距离判断该结点是否在目标树木的邻域影响圈内；c如果在影响圈内，则直接获取最底层叶子结点数据，精确判断叶子结点是否在目标树木的影响圈，如果是，则根据模型计算其NCI值，如果不是，则不计算，d如果线程所对应结点超出影响圈，则该线程任务结束，成年树光照GLI内核方法具体流程如下：a数据初始化：分配线程，将线程划分为g批进行处理，每批M个线程，线程总量为g*M，为成年树与幼年树木数量，动态建立实际分辨率所需样地单元网格二维数组记录结构assign_tree，对其所有元素初始化为‑1表示未分配任何树木ID，初始化所有成年树暴露面积A为0，获取整个森林场景中成年树与幼年树DBH，位置坐标，树种等属性数组，b获取线程号thread_id与块号block_id，通过thread_id与block_id索引树木索引tree_id，对于树木tree_id为l的树木，获取其DBH属性DBH（l），树木所在位置L（l），并通过模型公式由DBH（l）计算其高度H（l），树冠半径R（l），c通过L（l）获取树木l所在样地单元PID（l），以所在样地单元坐标PLocation（PID）为中心，遍历（x，y）坐标在x∈（L（PID）‑R（l），L（PID）+R（l）），y∈（L‑R（l），L（PID）+R（l））内的所有样地网格，若该范围内的样地单元网格k距离L（l）最远顶角坐标小于R（l），则该样地网格完全包含于树木1的投影圆内，进入步骤下面d，否则不作任何处理，d样地单元网格w完全包含于投影圆内，通过记录数组获取网格w所被分配的成年树ID:assign_tree（w），若该样地网格未被分配，则ID为‑1，此时只需将assign_tree（w）赋值为l，并对A（l）累加该样地单元面积r²，若该样地网格已被分配，则获取该样地被分配成年树ID的树高H（assign_tree（w）），比较H（assign_tree（w））与H（l）大小，若H（assign_tree（w））大于等于H（l），则不做任何处理，否则，将assign_tree（w）赋值为l，并对A（l）累加该样地单元面积r²，并对A（assign_tree（w））减去该样地单元面积r²，e继续进入步骤c遍历阴影范围内的样地网格，若遍历完，则进入成年树光照GLI内核方法具体流程中的步骤b进行下一棵成年树，直至结束该方法。