技术领域
[0001] 本
发明涉及地学、农学、林学和植物学等领域，具体涉及一种多层数据结构的虚植物冠层光合有效辐射分布模拟方法。
背景技术
[0002] 光合有效辐射(Photosynthetic Active Radiation,PAR)是
太阳辐射光谱中能够被绿色植物叶绿素吸收、转化为物质积累的那部分辐射，其辐射
波长介于380nm-710nm之间，约占太阳总
能量的40％。自太阳入射于植物冠层的PAR可分为太阳直射PAR、天空散射PAR及
叶片多次散射PAR，其中太阳直射及天空散射PAR是冠层内PAR的主要组成部分。植物叶片在PAR波段具有高吸收率、低反射及透射率的特点，因而植物冠层内的叶片多次散射辐射较小。植物叶器官多次散射辐射占冠层总辐射的比例小于10％，植物冠层PAR分布模拟时常忽略植物冠层的多次散射辐射效应。PAR是植物光合作用的
基础，植物冠层PAR截获及三维分布测量在农业生产管理、林业经营、科学研究等领域均有着广泛的应用。
[0003] 目前植物冠层光合有效辐射三维分布测量方法主要有三种：地面实测法、数学模型法、三维模型模拟法。地面实测法采用光合有效辐射观测设备(如SunScan、AccuPAR、SunFleck、TRAC、光合有效辐射计(光量子计)等测量冠顶及冠层内多个
位置瞬时或持续光合有效辐射值，并以此为基础计算植物冠层PAR三维分布，实地测量法由于受各种实测客观条件的限制，其测量空间
分辨率和
时间分辨率不高，因此该方法通常仅可在样地尺度或冠层各高度区间尺度实现植物冠层PAR三维分布测量。数学模型法则建立在经典的辐射传输模型基础上描述植物冠层内PAR分布规律，如比尔定律、一维模型等，由于数学模型法通常对植物冠层内PAR三维分布规律及树木冠层三维结构进行高度简化，因此其PAR三维分布模拟
精度不高。近年来，随着虚拟植物技术研究的不断深入，高逼真度虚拟植物
几何模型及虚拟植被场景建模成为可能，虚拟植物技术可真实再现植物冠层三维结构，采用
光线跟踪法或辐射度方法可模拟跟踪太阳光线在植物冠层内的反射、吸收及透射全过程，因而其可精确计算冠层内任意三维位置的PAR分布量。与其它两种方法相比，三维模型模拟法具有计算结果精度高、可重复多次使用、可模拟多种植被冠层特征、投入小、产出大等优点。
[0004] 常见的虚拟植物冠层PAR三维分布模拟以
光线跟踪法和辐射度方法为主。辐射度方法基于热辐射工程中的能量传递和守恒理论，其通过求解辐射系统方程，计算虚拟植物冠层内每个三
角面片的辐射分布情况。辐射度方法可完整模拟PAR与植物冠层之间交互的反射、透射、吸收和多重散射全过程，但其需对形状因子及线性方程组进行解算，因此存在计算量巨大、计算效率低、无法支持大场景、高复杂度场景等缺点。
光线跟踪算法是一种以物理光学原理为基础
模拟现实世界中光线传播的
计算机图形学基础算法，该算法具有原理简单、易于实现、真实感图形绘制及精度高等优点，采用
光线跟踪算法可模拟PAR与植物冠层之间的反射、透射和多次散射过程。尽管光线跟踪算法无需对大型线性方程组进行解算，但该方法也存在计算量巨大的不足，其计算量主要用于光线与场景面片求交判断及递归。传统光线跟踪算法开展光线跟踪模拟时其光线与场景要素求交判断需对场景内所有三角面片求交计算，因此计算效率很低。目前常见的四种光线跟踪
加速方法有四种：①.提高光线与场景基本单元的求交计算效率。传统光线跟踪算法有90％以上的计算时间均用于光线与三角面片、
包围盒、球体等场景最小单元的求交计算，因此光线与场景最小单元求交算法的优劣会极大地影响光线跟踪计算性能；②.减少光线计算数量。光线跟踪计算的复杂度不仅与场景复杂度、三角面片数目关系密切，其还与参与光线跟踪计算的光线数量密切相关，而减少光线跟踪数量则可能导致模拟精度降低；③.减少光线与场景基本单元相交计算次数。光线跟踪计算时根据场景基本单元的空间相关性，采用层次包围盒(Bounding Volume Hierarchies,BVH)、Kdtree、八叉树(Octree)及均匀网格等数据结构组织三维场景可快速判断与光线相交的场景对象，从而避免与无交场景基本单元求交判断，减少光线与包围盒、三角面片的相交计算次数，从而大幅提升光线跟踪计算效率。④.并行加速。光线跟踪计算时其每条光线的计算过程相似，因此使用CPU和GPU的并行处理能
力可大幅提高光线跟踪模拟效率。现有的光线跟踪加速方案中，鉴于光线与场景基本单元求交算法已较为完善及成熟，而减少光线计算数量的方法则会一定程度上降低光线跟踪模拟效果，因此目前光线跟踪加速方案以采用数据结构组织三维场景，以及采取并行和GPU技术提高光线跟踪效率较为常见。
[0005] 传统基于数据结构手段开展光线跟踪加速模拟的方法已广泛应用于游戏、动画、电影、地形及普通场景中，但较少有研究提及采用数据结构开展虚拟植被场景下光线跟踪加速模拟，同时虚拟植被场景与传统场景差异显著，因此需根据虚拟植被场景的特点有针对性设计其加速方案，以提高虚拟植被场景下光线跟踪效率。目前尚未有文献涉及采用多级数据结构开展虚拟植被场景下PAR模拟计算加速的报道，尤其针对复杂地形区虚拟植被冠层PAR三维分布模拟的情况。
发明内容
[0006] 本发明的目的在于提供一种基于多层数据结构的植物冠层光合有效辐射三维分布模拟加速方法，以克服
现有技术中存在的
缺陷。
[0007] 为实现上述目的，本发明的技术方案是：一种基于多层数据结构的虚拟植物冠层光合有效辐射三维分布模拟加速方法，按照如下步骤实现：
[0008] 步骤S1：获取典型单株植物结构特征参数和样地调查
数据库，采用虚拟单株植物建模
软件建立高逼真度、自然生长的单株植物精细三维模型；获取目标区域DEM、TIN及地形数据，并基于植被场景构建软件建立虚拟植被场景；
[0009] 步骤S2：根据样地地理位置和
气候条件，采用天文参数计算算法，计算得出不同时刻的冠顶辐射强度，或采用测量仪器实地、实时收集植被冠顶光合有效辐射强度值；
[0010] 步骤S3：对所述虚拟植被场景中除针叶簇单元之外的其它场景要素三角面片化或结构化，采用数据结构分别针对单株树木的几何模型、场景及地形中所包含的三角面片及针叶簇单元构建三级数据结构；
[0011] 步骤S4：创建太阳直射光线矢量平面，采用光线跟踪算法依次跟踪所述太阳直射光线矢量平面上各光线矢量，并对单株树木的几何模型、场景及地形对应的三级数据结构分别进行遍历；当光线与三角面片有交时，将吸收的太阳直射PAR赋予相交的三角面片，采用蒙特卡洛算法分别模拟太阳直射PAR及天空散射PAR与三角面片反射、透射及吸收的全过程，且均基于单株树木的几何模型、场景及地形三级数据结构分别遍历以判断反射、
透射光线是否与三角面片相交；
[0012] 步骤S5：确定天穹剖分分辨率，采用龟壳算法依次对场景内各三角面片循环，分别在各三角面片所在
水平面单位半球形成光线矢量集合，依次循环光线集合每条光线，并对单株树木的几何模型、场景及地形三级数据结构进行遍历，标记与面片有交的光线，统计各三角面片天空可见率，结合冠顶天空散射PAR和三角面片天空可见率计算各三角面片天空散射PAR；
[0013] 步骤S6：场景内各三角面片太阳直射PAR和天空散射PAR总和即为该三角面片光合有效辐射总值。
[0014] 在本发明一
实施例中，在所述步骤S1中，所述单株植物精细三维模型按照如下步骤建立：
[0015] 步骤S111：以所述典型单株植物结构特征参数为基础提炼植物冠型特征参数、枝干系统拓扑结构信息和叶几何信息；
[0016] 步骤S112：根据所述叶几何信息，采用NURBS、Mesh曲面或标准几何体及其组合建立叶或针叶簇精细三维模型；
[0017] 步骤S113：根据所述植物冠型特征参数以及所述枝干系统拓扑结构信息，基于所述虚拟单株植物建模软件，建立植物枝干系统三维模型，并把所述叶或针叶簇精细三维模型按预设角度和分布
密度挂接到枝干系统上，形成所述单株植物精细三维模型；
[0018] 所述典型单株植物结构特征参数包括：单株植物叶器官特征参数以及植物主体形态结构参数；所述单株植物叶器官特征参数包括：叶片倾角分布、方位角分布、展布面积、针叶簇中针叶束数量、针叶束长度、针叶束与针叶簇夹角、针叶束叶序、针叶束形状以及针叶束朝向；所述植物主体形态结构参数包括：树高、冠幅，主干分节数、分节长度、树干半径变化、初始分枝角度等主干参数、枝条层数序号、分枝长度、枝条基径、枝条直径变化、初始分枝角度、角度变化、枝条分段数、枝条分布密度、下一枝条层数序号以及叶片分布范围；
[0019] 所述样地调查数据库包括：胸径、植被高度、枝下高、冠幅以及植被主干地理坐标；
[0020] 所述虚拟单株植物建模
软件包括：xfrog、ParaTree、SpeedTree、AMAP以及L-System。
[0021] 在本发明一实施例中，在所述步骤S1中，按照如下步骤建立所述虚拟植被场景：
[0022] 步骤S121：加载数字地面模型、道路以及河流矢量图以形成地形、道路以及河流地形要素三维模型；
[0023] 步骤S122：从自然分布、均匀分布、随机分布及聚集分布中选择一种作为虚拟植被场景的树木分布模式；
[0024] 步骤S123：确定虚拟植被场景特征参数，包括：林分密度、样地大小、树木种类以及树木数量；
[0025] 步骤S124：结合所述树木分布模式、所述虚拟植被场景特征参数以及所述单株植物精细三维模型将单株植物三维模型按照虚拟植被场景构建方案逐一种植在地形上。
[0026] 在本发明一实施例中，，在所述步骤S2中，采用所述天文参数法计算算法，计算得出不同时刻的冠顶辐射强度包括以下步骤：
[0027] 步骤S211：根据所述样地地理位置，采用太阳天文参数计算公式得出太阳几何参数，包括：太阳高度角和太阳方位角；
[0028] 步骤S212：根据所述太阳几何参数、所述样地气候条件和地形条件，通过计算太阳光能在大气层顶和大气层中的强度，获取各时刻冠顶太阳直射辐射强度和天空散射辐射强度；
[0029] 步骤S213：将获取的各时刻冠顶太阳直射辐射强度和天空散射辐射强度分别乘以对应的比例系数，获取冠顶太阳直射PAR和天空散射PAR。
[0030] 在本发明一实施例中，在所述步骤S2中，采用测量仪器实地、实时收集植被冠顶光合有效辐射强度值，具体实施步骤包括：将植物冠层光合有效辐射观测设备置于冠顶，通过瞬时或持续观测测量所述植被冠顶光合有效辐射强度值；当采用持续观测模式时，所述植被冠顶光合有效辐射强度为观测时间段内各时刻观测值的平均值；所述植物冠层光合有效辐射观测设备包括：SunScan、AccuPAR、SunFleck、TRAC以及光合有效辐射计，且采用上述设备能同时得到太阳直射PAR和天空散射PAR参数值。
[0031] 在本发明一实施例中，其特征在于，在所述步骤S3中，
[0032] 所述虚拟植被场景单株树木几何模型、场景及地形三个尺度数据结构构建步骤如下：
[0033] 步骤S311：选择Kdtree作为所述单株树木几何模型、场景及地形三个尺度数据结构类型，分别将所述单株树木几何模型、场景及地形包围盒作为Kdtree根
节点，采用中分平面分割法确定初始分割平面；
[0034] 步骤S312：基于分割平面把根节点细分为左
子节点、右子节点，判断两子节点是否满足构建终止条件，即节点深度是否大于节点最大深度
阈值、节点的基本对象数是否小于最小对象阈值及节点光线遍历总代价是否大于其作为叶节点的光线遍历总代价；若满足，则转至步骤S313，否则，转至步骤S314；所述节点最大深度阈值根据场景对象总数自适应计算；
[0035] 步骤S313：生成Kdtree叶节点；
[0036] 步骤S314：采用中分平面分割法确
定子节点初始分割平面，并继续对Kdtree内部子节点进行细分，直到所有节点不需再细分为止，最终生成一个完整的Kdtree；
[0037] 所述针叶簇单元数据结构构建流程采用所述单株树木几何模型、场景及地形中三角面片数据结构构建流程，且所述针叶簇单元数据以针叶簇单元中心点为基础构建，其数据结构最小尺度单元为针叶簇。
[0038] 在本发明一实施例中，在所述步骤S4中，所述太阳直射光线矢量平面创建步骤如下：
[0039] 步骤S411：根据样地地理位置以及对应的时刻，采用太阳天文参数计算公式得出太阳几何参数，包括：太阳高度角和太阳方位角；
[0040] 步骤S412：根据太阳
光源三维位置距样地最近的边界的距离，将其三维坐标表示为 由太阳高度角和方位角计算太阳入射方向光线单位矢量基于 以及单位向量(0,1,0)计算太阳直射光线投射平面局
部
坐标系以及相应的
正交基向量 并得到太阳直
射光线投影水平面；
[0041] 步骤S413：将虚拟植物场景最小轴对齐包围盒八个
顶点，沿太阳入射方向反方向平行投影至所述太阳直射光线投影水平面，得到8个投影点，其最小包围长方形即为太阳直射光线投射平面；记太阳直射光线投射平面水平和垂直方向太阳直射光线投射数量分别为Hdres和Vdres，根据太阳直射光线投放密度将太阳直射光线投放平面均匀剖分成一系列均匀格网，以每个格网中心点作为太阳直射光线投射起点， 为投射光线单位矢量，形成光线方程，得到太阳直射光线矢量集合。
[0042] 在本发明一实施例中，在所述步骤S4中，采用光线跟踪算法开展虚拟植物冠层太阳直射PAR模拟步骤包括如下步骤：
[0043] 步骤S421：依次循环太阳直射光线矢量平面上各光线矢量，将各光线矢量逐一投射至虚拟植物场景，当各光线矢量均已投影跟踪结束时，虚拟植物冠层太阳直射PAR模拟结束；
[0044] 步骤S422：采用光线矢量平面上各光线矢量依次分别对场景、单株树木几何模型及地形三级场景数据结构及针叶簇数据结构进行遍历，判断光线与场景级数据结构内部节点最小轴对齐包围盒是否相交；若不相交，则判断光线与地形级数据结构是否相交，若不相交则结束此条光线矢量跟踪，并返回步骤421，若相交，则求解交点，采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程；当光线矢量与场景级数据结构内部节点相交时，确定相交的内部节点；
[0045] 步骤S423：开始对相交的内部节点数据结构遍历，判断光线与内部节点数据结构内单树几何模型是否相交，若不相交，则判断光线与地形级数据结构内部节点是否相交，若不相交，则结束此条光线矢量跟踪，返回步骤421，若相交，则求解交点，采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程；
[0046] 步骤S424：当光线与内部节点内单株树木几何模型三角面片及针叶簇数据结构相交时，确定相交的单株树木的几何模型，并对该单株树木几何模型数据结构遍历，判断光线与单株树木几何模型中三角面片或针叶簇单元是否相交，获取首次相交的三角面片或针叶簇；
[0047] 步骤S425：采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程；确定蒙特卡洛
采样频率参数，采用随机函数确定太阳直射PAR与相交的三角面片交互时的光线传播类型，包括：反射、透射或吸收，并将该太阳直射PAR能量赋予相交的三角面片或采用随机采样法确定光线下一步传播方向，同时计算其传播能量；根据光线传播方向及相交的三角面片建立光线传播方程，依次判断光线与场景及单树两级数据结构内部节点是否相交，若不相交，则返回步骤S421继续循环，直至循环次数达到设定的采用频率参数阈值；若相交，则进一步确定其相交的三角面片，并根据冠层要素吸收率函数计算冠层要素PAR吸收量，当光线能量小于设定的阈值，或者光线与冠层要素交互的次数大于设定的阈值时，则返回步骤S421继续循环，直至循环次数达到设定的采用频率参数阈值；
[0048] 步骤S426：返回步骤S421直至每条光线矢量均已跟踪，若存在首次相交的三角面片或针叶簇，则返回步骤S422继续循环。
[0049] 在本发明一实施例中，在所述步骤S5中，所述天空可见率计算步骤如下：
[0050] 步骤S511：对虚拟植被场景中所有三角面片及针叶簇循环，分别以每个三角面片或针叶簇中心点为半球中心点，对应创建单位上半球，将半球空间在天顶角及方位角方向分别等分为MH和NV个区间，形成MH×NV个均匀、相等大小的立体角网格，沿半球中心点及各立体角网格中心线建立各三角面片及针叶簇光线矢量集合；
[0051] 步骤S512：依次循环各三角面片及针叶簇光线矢量集合，采用光线矢量集合中各光线矢量对虚拟植被场景三角面片及针叶簇数据结构进行遍历，当光线矢量集合循环跟踪结束时，转至步骤S516；
[0052] 步骤S513：判断光线与虚拟植被场景数据结构内部节点最小轴对齐包围盒是否相交，若不相交，则结束此条光线矢量跟踪，并返回步骤S512；当光线矢量与虚拟植被场景三角面片及针叶簇数据结构内部节点相交时，确定相交的内部节点；
[0053] 步骤S514：对相交的内部节点数据结构遍历，判断光线与节点内单株树木的几何模型三角面片及针叶簇数据结构是否相交，若不相交，则返回步骤S513；当光线与节点内单株树木的几何模型三角面片及针叶簇数据结构相交时，确定相交的单株树木的几何模型；
[0054] 步骤S515：对相交的单树几何模型三角面片及针叶簇数据结构遍历，判断光线与单树几何模型数据结构中三角面片及针叶簇是否相交，若相交，则将该光线矢量记为与场景三角面片相交，返回步骤S512；
[0055] 步骤S516：统计每个三角面片及针叶簇光线矢量集合中与场景三角面片及针叶簇不相交的光线数量，该光线数量与光线矢量集合中光线数量的比值即为该三角面片或针叶簇的天空可见率。
[0056] 在本发明一实施例中，在所述步骤S5中，虚拟植被场景内三角面片天空散射PAR为虚拟植物冠顶天空散射PAR与该三角面片天空可见率的乘积
[0057] 相较于现有技术，本发明具有以下有益效果：本发明所提出的一种虚拟植物冠层光合有效辐射三维分布模拟加速方法，克服了传统三维模型模拟方法计算时间较长、资源消耗巨大、无法支持高复杂度虚拟植被场景，尤其是地形区大尺度、大范围针叶林植被场景的缺陷，大幅提高了虚拟植被冠层PAR三维分布模拟效率，并有效降低了资源消耗及需求，同时可支持高复杂度虚拟植被场景，为虚拟植物冠层光合有效辐射三维分布模拟提供稳定、高效、高精度的计算结果。
附图说明
[0058] 图1为本发明基于多层数据结构的虚拟植物冠层光合有效辐射三维分布模拟加速方法的
流程图。
[0059] 图2为本发明采用三层数据结构、两层数据结构和未包含数据结构的虚拟阔叶
林场景PAR三维部分模拟效率对比结果图。
具体实施方式
[0060] 下面结合附图，对本发明的技术方案进行具体说明。
[0061] 本发明提供一种基于多层数据结构的虚拟植物冠层光合有效辐射三维分布模拟加速方法，包括以下步骤实现：
[0062] 步骤S1：收集典型单株植物结构特征参数和样地调查数据库，采用虚拟植物建模软件建立高逼真度、自然生长的单株植物精细三维模型。收集目标区域DEM、TIN及地形数据，基于植被场景构建软件创建虚拟植被场景；
[0063] 步骤S2：根据样地地理位置和气候条件，采用相关天文参数计算算法，计算得出不同时刻的冠顶辐射强度；或者采用测量仪器实地、实时收集植被冠顶光合有效辐射强度值；
[0064] 步骤S3：对虚拟植被场景中除针叶簇单元之外的其它场景要素三角面片化或结构化，采用数据结构分别针对单株树木几何模型、场景及地形中所包含的三角面片及针叶簇单元构建三级数据结构；
[0065] 步骤S4：创建太阳直射光线矢量平面，采用光线跟踪算法依次跟踪光线平面上各光线矢量，并对单树、场景及地形三级数据结构分别进行遍历，当光线与三角面片有交时，将吸收的太阳直射PAR赋予相交的三角面片。采用蒙特卡洛算法模拟太阳直射PAR及天空散射PAR与三角面片反射、透射及吸收的全过程，其均基于单树、场景及地形三级数据结构分别遍历以判断反射、透射光线是否与三角面片相交；
[0066] 步骤S5：确定天穹剖分分辨率，采用龟壳算法依次对场景内各三角面片循环，分别在各三角面片所在水平面单位半球形成光线矢量集合，依次循环光线集合每条光线，并对单树、场景及地形三级数据结构进行遍历，标记与面片有交的光线，统计各三角面片天空可见率，结合冠顶天空散射PAR和三角面片天空可见率计算各三角面片天空散射PAR；
[0067] 步骤S6：场景内各三角面片太阳直射PAR和天空散射PAR总和即为该三角面片光合有效辐射总值；
[0068] 进一步的，在本实施例中，步骤S1中的收集典型单株植物结构特征参数，其包括单株植物叶器官和植物主体形态结构两方面，叶器官特征参数收集包括叶片倾角分布、方位角分布、展布面积、针叶簇中针叶束数量、针叶束长度、针叶束与针叶簇夹角、针叶束叶序、针叶束形状、针叶束朝向等叶片几何建模及分枝特征参数，主体形态结构参数收集主要包括树高、冠幅，主干分节数、分节长度、树干半径变化、初始分枝角度等主干参数，以及枝条层数序号、分枝长度、枝条基径、枝条直径变化、初始分枝角度、角度变化、枝条分段数、枝条分布密度、下一枝条层数序号、叶片分布范围等枝条参数。
[0069] 进一步的，在本实施例中，步骤S1中的样地调查数据库，其主要包括胸径、植被高度、枝下高、幅宽、植被主干地理坐标等虚拟植被环境构建所需的植被样地数据。
[0070] 进一步的，在本实施例中，在步骤S1中，建立高逼真度、自然生长的单株植物精细三维模型包括以下步骤：
[0071] 步骤S111：以收集的单株植物结构特征参数为基础提炼植物冠型特征参数、枝干系统拓扑结构信息和叶几何信息；
[0072] 步骤S112：根据所述叶几何信息，采用NURBS、Mesh曲面或标准几何体及其组合建立叶或针叶簇精细三维模型；
[0073] 步骤S113：根据所述冠层特征参数和枝干系统拓扑结构信息，基于单株植物建模软件如xfrog、ParaTree、SpeedTree、AMAP、L-System等软件平台，建立植物枝干系统三维模型，并把叶或针叶簇精细三维模型按一定的角度和分布密度挂接到枝干系统上，形成高逼真度、自然生长的单株植物精细三维模型。
[0074] 进一步的，在本实施例中，在步骤S1中，建立虚拟植被场景包括以下步骤：
[0075] 步骤S121：加载数字地面模型、道路及河流矢量图以形成地形、道路、河流等地形要素三维模型；
[0076] 步骤S122：从自然分布、均匀分布、随机分布及聚集分布等树木分布模式中选择一种作为虚拟植被场景树木分布模式；
[0077] 步骤S123：确定虚拟植被场景特征参数，如林分密度、样地大小、树木种类、树木数量等；
[0078] 步骤S124：结合树木分布模式、虚拟植被场景特征参数和单株植物精细三维模型库将单株植物三维模型按照虚拟植被场景构建方案逐一种植在地形上。
[0079] 进一步的，在本实施例中，在步骤S2中，采用天文参数法计算冠顶辐射强度包括以下步骤：步骤S211：根据样地地理位置，采用太阳天文参数计算公式得出太阳几何参数，即太阳高度角和太阳方位角；
[0080] 步骤S212：根据太阳几何参数、样地气候条件和地形条件，采用太阳光能在大气层顶和大气层中的计算公式，计算各时刻冠顶太阳直射辐射强度和天空散射辐射强度；
[0081] 步骤S213：各时刻冠顶太阳直射辐射强度和天空散射辐射强度分别乘以比例系数，得出冠顶太阳直射PAR和天空散射PAR。
[0082] 进一步的，在本实施例中，在步骤S2中，采用仪器法测量植被冠顶光合有效辐射强度值，其主要指将植物冠层光合有效辐射观测设备置于冠顶，通过瞬时或持续观测测量植被冠顶光合有效辐射强度值，当采用持续观测模式时，其冠顶光合有效辐射强度为观测时间段内各时刻观测值的平均值；植物冠层光合有效辐射观测设备主要包括SunScan、AccuPAR、SunFleck、TRAC、光合有效辐射计(光量子计)及自主研究设备等，采用上述设备通常可同时得到太阳直射PAR和天空散射PAR参数值。
[0083] 进一步的，在本实施例中，在步骤S3中，虚拟植被场景单株树木几何模型、场景及地形三个尺度数据结构构建步骤如下：
[0084] 以单株树木几何模型、场景及地形三个尺度数据结构类型为Kdtree为例，分别将单树几何模型、地形及虚拟植被场景包围盒作为Kdtree根节点，采用中分平面分割法确定初始分割平面；基于分割平面把根节点细分为左子节点、右子节点，然后判断两子节点是否满足构建终止条件，即节点深度是否大于设定阈值、节点的基本对象数是否小于设定阈值及节点光线遍历总代价是否大于其作为叶节点的光线遍历总代价，若满足，则生成Kdtree叶节点，否则，采用中分平面分割法确定子节点初始分割平面并继续对Kdtree内部子节点进行细分，直到所有节点不需再细分为止，最终生成一个完整的Kdtree；节点最大深度阈值根据场景对象总数自适应计算获得。
[0085] 进一步的，在本实施例中，针叶簇单元数据结构构建流程和单株树木几何模型、场景及地形中三角面片数据其数据结构构建流程一致，但针叶簇单元数据结构构建需分两个步骤分别完成，三角面片数据数据结构构建以三角面片中心点为基础，而针叶簇单元数据则以针叶簇单元中心点为基础构建，其数据结构最小尺度单元为针叶簇。
[0086] 进一步的，在本实施例中，步骤S4中的太阳直射光线矢量平面创建步骤如下：
[0087] 步骤S411：根据样地地理位置和时间，采用太阳天文参数计算公式得出太阳几何参数，即太阳高度角和太阳方位角；
[0088] 步骤S412：设太阳光源三维位置距样地最近的边界一定距离，其三维坐标可表示为 由太阳高度角和方位角计算太阳入射方向光线单位矢量基于 和单位向量(0,1,0)计算太阳直射光线投射平面局部
坐标系及相应的正交基向量 并得到太阳直射光
线投影水平面；步骤S413：将虚拟植物场景最小轴对齐包围盒八个顶点沿太阳入射方向反方向平行投影至太阳直射光线投影水平面，得到8个投影点，其最小包围长方形即为太阳直射光线投射平面；假设太阳直射光线投射平面水平和垂直方向太阳直射光线投射数量分别为Hdres和Vdres，依据太阳直射光线投放密度将太阳直射光线投放平面均匀剖分成一系列均匀格网，以每个格网中心点作为太阳直射光线投射起点， 为投射光线单位矢量，形成光线方程，从而得到太阳直射光线矢量集合。
[0089] 进一步的，在本实施例中，在步骤S4中的采用光线跟踪算法模拟虚拟植物冠层太阳直射PAR步骤如下：
[0090] 步骤S421：依次循环太阳直射光线矢量平面上各光线矢量，将各光线矢量逐一投射至虚拟植物场景，当各光线矢量均已投影跟踪结束时，虚拟植物冠层太阳直射PAR模拟结束；
[0091] 步骤S422：采用光线矢量平面上各光线矢量依次分别对场景、单树及地形三级场景数据结构及针叶簇数据结构进行遍历，判断光线与场景级数据结构内部节点最小轴对齐包围盒是否相交，如不相交，则判断光线与地形级数据结构是否相交，如不相交，则结束此条光线矢量跟踪，返回步骤421，如相交，则求解交点，采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程；当光线矢量与场景级数据结构内部节点相交时，确定相交的内部节点；
[0092] 步骤S423：开始对相交的内部节点数据结构遍历，判断光线与内部节点数据结构内单树几何模型是否相交，如不相交，则判断光线与地形级数据结构内部节点是否相交，如不相交，则结束此条光线矢量跟踪，返回步骤421，如相交，则求解交点，采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程；
[0093] 步骤S424：当光线与节点内单树几何模型三角面片及针叶簇数据结构相交时，确定相交的单树几何模型；开始对相交的单树几何模型数据结构遍历，判断光线与单树几何模型中三角面片或针叶簇单元是否相交，并找出最近的三角面片或针叶簇，采用蒙特卡洛算法模拟太阳直射PAR与相交的三角面片反射、透射及吸收的全过程。
[0094] 步骤S425：确定蒙特卡洛
采样频率参数，采用随机函数确定太阳直射PAR与相交的三角面片交互时的光线传播类型，如反射、透射或吸收，如光线传播类型为吸收则将该太阳直射PAR能量赋予相交的三角面片，否则采用随机采样法确定光线下一步传播方向，同时计算其传播能量；根据光线传播方向及相交的三角面片建立光线传播方程，依次判断光线与场景及单树两级数据结构内部节点是否相交，如不相交，则返回步骤S421继续循环，直至循环次数达到设定的采用频率参数阈值；如相交则进一步确定其相交的三角面片，并根据冠层要素吸收率函数计算冠层要素PAR吸收量，当光线能量小于设定的阈值，或者光线与冠层要素交互的次数大于设定的阈值时，返回步骤S421继续循环，直至循环次数达到设定的采用频率参数阈值。
[0095] 步骤S426：返回步骤S421直至每条光线矢量均已跟踪，若存在最近相交的三角面片或针叶簇，则返回步骤S422继续循环；
[0096] 进一步的，在本实施例中，步骤5中的天空可见率计算步骤如下：
[0097] 步骤S511：对虚拟植被场景中所有三角面片及针叶簇循环，分别以每个三角面片或针叶簇中心点为半球中心点创建一单位上半球(半球直径为1)，将半球空间在天顶角及方位角方向分别等分为MH和NV个区间，形成MH×NV个均匀、相等大小的立体角网格，沿半球中心点及各立体角网格中心线建立各三角面片及针叶簇光线矢量集合；
[0098] 步骤S512：依次循环各三角面片及针叶簇光线矢量集合，采用光线矢量集合中各光线矢量对虚拟植被场景三角面片及针叶簇数据结构进行遍历，当光线矢量集合循环跟踪结束时跳至步骤S516；
[0099] 步骤S513：判断光线与虚拟植被场景数据结构内部节点最小轴对齐包围盒是否相交，如不相交则结束此条光线矢量跟踪，并返回步骤S512；当光线矢量与虚拟植被场景三角面片及针叶簇数据结构内部节点相交时，确定相交的内部节点；
[0100] 步骤S514：开始对相交的内部节点数据结构遍历，判断光线与节点内单树几何模型三角面片及针叶簇数据结构是否相交，如不相交则返回步骤S513；当光线与节点内单树几何模型三角面片及针叶簇数据结构相交时，确定相交的单树几何模型；
[0101] 步骤S515：开始对相交的单树几何模型三角面片及针叶簇数据结构遍历，判断光线与单树几何模型数据结构中三角面片及针叶簇是否相交，如相交则将该光线矢量记为与场景三角面片相交，返回步骤S512；
[0102] 步骤S516：统计每个三角面片及针叶簇光线矢量集合中与场景三角面片及针叶簇不相交的光线数量，该数量与光线矢量集合中光线数量的比值即为该三角面片或针叶簇的天空可见率。
[0103] 进一步的，在本实施例中，虚拟植被场景内三角面片天空散射PAR为虚拟植物冠顶天空散射PAR与该三角面片天空可见率的乘积；
[0104] 为了让本领域技术人员进一步了解本发明所提出的一种基于多层数据结构的虚拟植物冠层光合有效辐射三维分布模拟加速方法，下面结合具体实例进行说明。
[0105] 如图1所示，以一个阔叶林样地为例，开展三层数据结构的虚拟植物冠层光合有效辐射三维分布模拟加速，其具体步骤如下：
[0106] 1、虚拟植被场景构建。
[0107] 1)根据植物几何
参数化建模所需的植物形态结构参数要求，通过野外实测数据和经典文献数据的收集，获取自然生长形冠型的形态结构参数和样地调查数据库，主要包括叶片倾角分布、方位角分布、展布面积、针叶簇中针叶束数量、针叶束长度、针叶束与针叶簇夹角、针叶束叶序、针叶束形状、针叶束朝向等叶片几何建模及分枝特征参数，表1为阔叶林典型叶器官参数表，主体形态结构参数主要包括树高、冠幅，主干分节数、分节长度、树干半径变化、初始分枝角度等主干参数，以及枝条层数序号、分枝长度、枝条基径、枝条直径变化、初始分枝角度、角度变化、枝条分段数、枝条分布密度、下一枝条层数序号、叶片分布范围等枝条参数，并形成典型统计值，表2为阔叶林树种主要特征参数表。
[0108] 表1
[0109]
[0110] 表2
[0111]
[0112] 2)叶器官建模。采用XFrog建模软件和叶器官统计值建立阔叶树叶器官几何模型。采用3dx MAX软件或者自编程序，结合标准几何体，通过组合、复制、旋转等手段建立
云杉林针叶簇几何模型。依据树木主干、分枝分布规律，采用XFrog建模软件建立单株植物枝干精细几何模型，根据叶器官分布规律将叶器官几何模型挂接于枝干系统，调整单树几何模型形态，形成高逼真度、自然生长的单株植物精细三维模型；
[0113] 3)虚拟植被场景构建。在VisForest软件中创建水平地形面，根据样地调查数据库，将各单株树木几何模型按照其实际生长位置种植在地形上，形成虚拟阔叶林植被场景；
[0114] 2、虚拟植被场景冠顶太阳直射和天空散射PAR计算。
[0115] 设定阔叶林样地经纬度为：北纬58°16′39″，东经27°20′30″，假设PAR模拟日期为2014年3月20日，观测时间为正午12:00，根据天文参数计算方法，虚拟阔叶林样地冠顶太阳直射和天空散射PAR值分别为：1.956、663.73μmol*m-2*s-1。
[0116] 3、单树几何模型、地形及虚拟植被场景三级数据结构构建。
[0117] 单株树木几何模型、场景及地形三个尺度数据结构构建方法选择Kdtree，Kdtree中分平面分割方法为启发式排序法。
[0118] 针对虚拟植被场景场景级数据结构构建，其基本数据结构单元为场景内单株树木三角面片最小包围盒及其中心点，数据结构构建则以其三角面片最小包围盒及其中心点为基础。首先将虚拟植被场景三角面片包围盒作为Kdtree根节点，采用启发式排序法确定初始分割平面，基于分割平面把根节点细分为左子节点、右子节点，然后判断两子节点是否满足构建终止条件，即节点深度是否大于设定阈值、节点的基本对象数是否小于设定阈值及节点光线遍历总代价是否大于其作为叶节点的光线遍历总代价(虚拟植被场景节点基本对象数阈值为1，虚拟植被场景数据结构子节点数为2)，若满足，则生成Kdtree叶节点，否则，采用启发式排序法确定子节点初始分割平面并继续对Kdtree内部子节点进行细分，直到所有节点不需再细分为止，最终生成一个完整的以单株树木三角面片最小包围盒及其中心点为基本单元的虚拟植被场景Kdtree数据结构；节点最大深度阈值根据场景对象总数自适应计算获得，其计算公式为：
[0119] Dmax＝k1+k2*log2 N (1)
[0120] 其中Dmax为Kdtree最大构建深度阈值，N为场景对象总数，k1、k2为常量，当数据结构构建对象为虚拟植被场景时，k1、k2值分别设为0.667、0.667；节点和叶节点光线遍历总代价计算方法则可采用局部贪心表面积启发法计算得到。
[0121] 单株树木几何模型中的三角面片其数据结构构建步骤和虚拟植被场景级数据结构构建步骤相同，但采用中分平面分割子节点时终止条件参数值设置不同，即单树几何模型及地形基本节点数阈值为15，单树几何模型及地形数据结构子节点数为2，同时节点最大深度阈值计算时其k1、k2值分别设为1、1。
[0122] 针对虚拟植被场景中的地形三角面片开展数据结构构建，其步骤和虚拟植被场景数据结构构建步骤相同，但采用中分平面分割子节点时终止条件参数值设置不同，即单树几何模型基本节点数阈值为1，地形数据结构子节点数为2，同时节点最大深度阈值计算时其k1、k2值分别设为0.667、0.667。
[0123] 4、虚拟阔叶林样地太阳直射PAR模拟。
[0124] 1)创建太阳直射光线矢量平面。根据样地地理位置和时间，采用太阳天文参数计算公式得出太阳几何参数，即太阳高度角和太阳方位角分别为90°和0°；设太阳光源三维位置距样地最近的边界一定距离，其三维坐标可表示为 由太阳高度角和方位角计算太阳入射方向光线单位矢量 基于
和单位向量(0,1,0)计算太阳直射光线投射平面局部坐标系及太阳直射光线投影水平面；
将虚拟植物场景最小轴对齐包围盒八个顶点沿太阳入射方向反方向平行投影至太阳直射光线投影水平面，得到8个投影点，其最小包围长方形即为太阳直射光线投射平面；假设太阳直射光线投射平面水平和垂直方向太阳直射光线投射数量分别为3000和3000，依据太阳直射光线投放密度将太阳直射光线投放平面均匀剖分成一系列均匀格网，以每个格网中心点作为太阳直射光线投射起点， 为投射光线单位矢量，形成光线方程，从而得到太阳直射光线矢量集合；
[0125] 2)虚拟阔叶林场景太阳直射PAR模拟。
[0126] 步骤A1：依次循环太阳直射光线矢量平面上各光线矢量，将各光线矢量逐一投射至虚拟植物场景，当各光线矢量均已投射跟踪结束时，虚拟阔叶林场景太阳直射PAR模拟结束，虚拟阔叶林场景三角面片太阳直射PAR计算结果之和即为虚拟阔叶林场景太阳直射PAR模拟最终结果，本案例太阳直射PAR模拟未在地形级数据结构部分多次散射模拟；
[0127] 步骤A2：逐一采用光线矢量平面上各光线矢量依次对场景、单树及地形三级场景数据结构Kdtree数据结构进行遍历；
[0128] 步骤A3：判断光线与场景级数据结构是否相交，如相交则开始判断光线与单树级数据结构是否相交，否则开始光线与地形级数据结构相交判断，求解光线与地形级数据结构相交的内部节点，并将吸收的太阳直射PAR能量赋给相交的三角面片，返回步骤A1。
[0129] 步骤A4：开始光线与单树级数据结构是否相交判断，如相交则求解相交的内部节点，否则开始光线与地形级数据结构相交判断，求解光线与地形级数据结构相交的内部节点，并将吸收的太阳直射PAR能量赋给相交的三角面片，返回步骤A1。
[0130] 步骤A5：根据光线与单树级数据结构相交的内部节点求解光线与单树级数据结构相交的交点，判断光线与冠层要素散射次数是否已达到阈值，如已达到阈值则将吸收的太阳直射PAR能量赋给相交的三角面片，返回步骤A1。当光线与冠层要素散射次数未达到阈值时，判断光线能量衰减是否达到阈值，如已达到阈值则将吸收的太阳直射PAR能量赋给相交的三角面片，返回步骤A1。
[0131] 步骤A6：当光线与冠层要素散射时其光线能量衰减未达到阈值，设置冠层要素反射率、透射率及吸收率参数分别为0.04、0.01、0.95，蒙特卡洛采样频率参数设为10000，利用随机函数确定太阳直射PAR与相交的三角面片交互时的光线传播类型，利用随机数产生0-1区间的随机数，当随机数小于0.95时则光线被吸收，并将吸收的太阳直射PAR能量赋给相交的三角面片，返回步骤A1。当随机数大于0.95并小于0.99时光线被反射，否则光线被透射。当光线与三角面片交互被反射或透射时，采用双随机数确定光线下一步传播的方位和天顶方向，同时计算其传播能量；根据光线传播方向及相交的三角面片建立光线传播方程，逐一对场景、单树及地形三级场景数据结构及针叶簇Kdtree数据结构进行遍历，返回步骤A2。
[0132] 5、虚拟阔叶林样地太阳散射PAR模拟。
[0133] 步骤B1：对虚拟植被场景中所有三角面片循环，以每个三角面片中心点为半球中心点创建一单位上半球(半球为1)，将半球空间在天顶角及方位角方向分别等分为90和360个区间，形成90﹡360个均匀、相等大小的立体角网格，沿半球中心点及各立体角网格中心线建立各三角面片光线矢量集合，当虚拟植被场景中所有三角面片单元均已循环结束时，跳至步骤B8，本案例太阳散射PAR模拟未考虑多次散射模拟；
[0134] 步骤B2：依次循环各三角面片单元对应的光线矢量集合，采用光线矢量集合中各光线矢量对虚拟植被场景Kdtree数据结构进行遍历，当光线矢量集合循环跟踪结束时跳至步骤B7；
[0135] 步骤B3：判断光线与虚拟植被场景Kdtree数据结构内部节点最小轴对齐包围盒是否相交，如不相交则结束此条光线矢量跟踪，并返回步骤B2；当光线矢量与虚拟植被场景三角面片Kdtree数据结构内部节点相交时，确定相交的内部节点；
[0136] 步骤B4：开始对相交的内部节点数据结构遍历，判断光线与节点内单树几何模型三角面片Kdtree数据结构是否相交，如相交则确定相交的单树几何模型，如光线与单树几何模型三角面片Kdtree数据结构相交，则开始对相交的单树几何模型三角面片Kdtree数据结构遍历，判断光线与单树几何模型中三角面片是否相交，如相交则将该光线矢量记为与场景三角面片相交，结束此条光线矢量跟踪，若不存在相交的三角面片，则将该光线矢量记为与地形级数据结构相交，返回步骤B3继续循环；
[0137] 步骤B5：统计与三角面片单元、地形级数据结构三角面片相交的光线矢量数量，该数量与光线矢量集合中光线数量的比值即为该三角面片的天空可见率，返回步骤B1；
[0138] 步骤B6：虚拟植被场景中各三角面片单元太阳散射PAR值为冠顶太阳散射PAR数值与该三角面片单元天空可见率的乘积。
[0139] 如图2为采用三层数据结构、两层数据结构和未包含数据结构的虚拟阔叶林场景PAR三维部分模拟效率对比结果图。
[0140] 6、虚拟阔叶林样地冠层PAR三维分布计算。
[0141] 虚拟阔叶林场景内各三角面片单元太阳直射和天空散射PAR总和即为该三角面片单元PAR，从而得到虚拟阔叶林样地PAR三维分布最终结果。
[0142] 上述实例采用Kdtree和启发式排序法分别作为场景数据结构构建的数据结构方法及中分平面分割方法，但本实施例中所采用的数据结构构建方法不受所述实例的限制，如BVH、八叉树、立体体元模型等数据结构，中分平面分割方法则有空间中分法、对象中分法、启发式排序法及启发式BIN分区法等，其它的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
[0143] 上述实例采用平地(地形坡度为0)作为地形级几何模型数据结构建模的对象，但本实施例的具体实施不受地形坡度的影响及限制，其它的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。
[0144] 以上是本发明的较佳实施例，凡依本发明技术方案所作的改变，所产生的功能作用未超出本发明技术方案的范围时，均属于本发明的保护范围。
