发明名称 一种基于GPU加速的多分辨率大规模森林演替过程仿真方法
摘要 一种基于GPU加速的多分辨率大规模森林演替过程仿真方法,包括以下步骤:1),建立树木数据组织与样地数据组织;根据GPU架构对树木数据进行组织,形成利于GPU计算的数数据组织结构,对于整个样地单元建立一与样地网格大小相同的二维数组结构M;2)将每棵树木的遍历转换成cuda内核中单位线程,计算plot[k]的种子分布情况,Y纬度值代表样地,对线程网格的纬度Y=k处,所有X纬度上的线程进行对样地K的种子分布影响计算;3)种子分布内核方法,成年树NCI方法GPU的内核方法以及成年树光照GLI内核方法。本发明提供一种正确性良好、提升运行效率的基于GPU加速的多分辨率大规模森林演替过程仿真方法。
申请公布号 CN102646261A 申请公布日期 2012.08.22
申请号 CN201210046135.2 申请日期 2012.02.27
申请人 浙江工业大学 发明人 范菁;嵇海锋;汤颖;董天阳
分类号 G06T1/00(2006.01)I;G06T1/20(2006.01)I 主分类号 G06T1/00(2006.01)I
代理机构 杭州天正专利事务所有限公司 33201 代理人 王兵;王利强
主权项 一种基于GPU加速的多分辨率大规模森林演替过程仿真方法,其特征在于:所述仿真方法包括以下步骤:1),建立树木数据组织与样地数据组织;根据GPU架构对树木数据进行组织,形成利于GPU计算的数数据组织结构:所有计算的树木数据结构将被设计成数组形式;在多分辨率数据组成中聚类结点第k层与第k+1层的数据表示,对于第k层单位结点数据由m个k‑1层结点数据累加而成,对于每个单位结点数据,由s种树种对应单位成年树结点数据构成;每棵成年树或者幼年树x的树冠被建模成圆柱体,它们在地表的投影都是一个圆C(x),设样地单元边长为r,则样地单元网格的面积即r2,对于投影完全包含k个样地单元网格的树木,则其暴露面积即为kr2,为正确计算不同树木发生遮挡时暴露面积,设树木树冠上顶面离地面距离为H(x),所有样地网格单元标记为Mij,在理想不想交投影区域,则如果某样地网格Mij完全包含在树木投影C(x)中,则将Mij分配给树木x,否者不对其分配;在投影相交区域,对于样地网格Mij被多个投影所覆盖,该样地网格将会被分配给能够完全包含它的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): <mrow> <msub> <mi>R</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>&eta;</mi> </mfrac> <mi>STR</mi> <munderover> <mi>&Sigma;</mi> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mrow> <mo>(</mo> <msup> <mrow> <mo>(</mo> <mfrac> <msub> <mi>DBH</mi> <mi>k</mi> </msub> <mn>30</mn> </mfrac> <mo>)</mo> </mrow> <mi>&beta;</mi> </msup> <msup> <mi>e</mi> <msubsup> <mrow> <mo>-</mo> <mi>Dd</mi> </mrow> <mi>ik</mi> <mi>&theta;</mi> </msubsup> </msup> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>其中,Ri表示样地i的种子数量,μ,STR,D,θ,β都是与树种相关的常数参数,DBH为成年树树干胸径大小,由成年树胸径的大小及树种参数所决定,dik表示第k棵成年树距离样地i的实际物理距离,N为所有具有繁殖能力的成年树总数,通过对树种真实参数的代入可知,随着母树与样地的距离的dik不断增大,母树所对样地产生的种子数量Ri成单调递减趋势;将每棵树木的遍历转换成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): <mrow> <mi>NCI</mi> <mo>=</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>S</mi> </msubsup> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>j</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </msubsup> <msub> <mi>&lambda;</mi> <mi>i</mi> </msub> <mfrac> <msup> <mrow> <mo>(</mo> <msub> <mi>DBH</mi> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mi>&alpha;</mi> </msup> <msup> <mrow> <mo>(</mo> <msub> <mrow> <mi>dis</mi> <mi>tan</mi> <mi>ce</mi> </mrow> <mi>ij</mi> </msub> <mo>)</mo> </mrow> <mi>&beta;</mi> </msup> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>计算公式中,i=1...S为树种数,j=1...N领域内大于某一固定DBH的树木;α是由目标树树种决定的NCIAlpha参数值;β是由目标树树种决定的NCIBeta参数值;DBHjk是领域树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)累加该样地单元面积r2,若该样地网格已被分配,则获取该样地被分配成年树ID的树高H(assign_tree(k)),比较H(assign_tree(k))与H(i)大小,若H(assign_tree(k))大于等于H(i),则不做任何处理,否者,将assign_tree(k)赋值为i,并对A(i)累加该样地单元面积r2,并对A(assign_tree(k))减去该样地单元面积r2,e继续进入步骤c遍历阴影范围内的样地网格,若遍历完,则进入步骤b进行下一棵成年树,直至结束该方法。
地址 310014 浙江省杭州市下城区朝晖六区