利用辐射度模拟虚拟枇杷果园冠层光分布

杨怡斐，唐丽玉，崔 磊，陈舒炜

YANG Yifei1，2,TANG Liyu1，2,CUI Lei1，2,CHEN Shuwei1，2

1.福州大学 空间数据挖掘与信息共享教育部重点实验室，福州 350116

2.福州大学 地理空间信息技术国家地方联合工程研究中心，福州 350116

1.Key Laboratory of Spatial Data Mining&Information Sharing of MOE,Fuzhou University,Fuzhou 350116,China

2.National Engineering Research Center of Geospatial Information Technology,Fuzhou University,Fuzhou 350116,China

1 引言

果树形态结构决定了冠层内的太阳光的分布情况，而光分布影响了光合作用效率和冠层内微环境，进一步影响果实产量和质量。目前冠层内光环境量化分析方法主要有地面实测法[1]、数学模型法[2]和三维植物结构模拟法[3]。地面实测法受实际客观条件限制，空间分辨率和时间分辨率有限；数学模拟法一般对真实冠层进行简化或某种假设，采用一些函数描述叶倾角分布和方位角[4]，在遥感反演得到较多应用，但与真实冠层还有较大误差。基于三维冠层结构的模拟方法已受到学者们的青睐，三维冠层结构可以更精确地描述冠层组分之间的拓扑结构，可以采用更精细面元表达冠层的三维信息。

辐射度算法和光线跟踪算法是三维植物冠层模拟辐射传输的主要方法。例如，文献[5]把冠层三维结构离散化为小立方体集合，同时使用高斯迭代法为每个立方体赋予光属性的方法模拟离散的光线集合。文献[6]利用嵌套辐射度方法，使用蒙特卡洛框架模拟自然光在植物冠层的三维分布。以上模拟已达到植物器官级别或更小单元，因此这类模型可以模拟冠层内任一点的冠层光分布或光强。文献[7]提出了RAPID模型，可以模拟复杂场景中任意的植被冠层，进而快速计算二向反射因子（BRF），适用于不同分辨率的卫星数据。文献[8]从叶片、枝梢、树枝以及整株四个尺度上研究树形结构对冠层光截获的影响，并引入评价光截获效率的STAR（Silhouette to Total Area Ratio）值，即冠层或枝梢叶片被太阳照射的叶面积与总叶面积之比。文献[9-11]分别以黄瓜、玉米、烟草冠层为例分析了植物冠层内的光传输。近年来辐射度模型应用于植物冠层内光分布模拟多见于单株尺度，而果园尺度的研究对象为果树群体，由单株走向群体是模型应用中的必然步骤。在植物群体尺度上，一些学者采用光线跟踪方法[12-13]模拟光分布，文献[12]未考虑冠层内部相互反射（如枝条、果实等）、叶片透射等情况，模拟苹果树群体冠层光分布；文献[13]利用了开源光线跟踪模拟平台GroIMP，且用递归方法计算光线在榉木群体冠层内的多次反射。以上两篇文献的不足之处在于未能考虑植株所处地形的反射光对冠层光分布的影响。另一些学者在此基础上，加入了植物周围环境对群体冠层内光分布的影响，如文献[14]以玉米为例，利用辐射度-图形学结合模型（RGM）计算玉米群体冠层内光分布情况，且考虑了玉米叶片、茎秆面元的反射率与透射率以及地形的影响。不足在于假定地形为平面，与实际情况不符。文献[15]在温室内人工光环境条件下，综合考虑太阳直射、天空散射、人工光源模拟了现代月季群体冠层内光分布与光截获，并计算了温室内地板、墙壁与屋顶的反射光对群体冠层光分布的影响，但其研究对象冠层结构简单，未能应用于复杂的果树冠层。

目前植物群体光环境模拟，大多未考虑实际植物所处地形表面的起伏对光分布的影响。相关研究表明，地形因子对BRF分布有着非常明显的影响[16]。地形是果树的主要载体，地面反射光在冠层内的分布情况是研究虚拟枇杷果园冠层光分布必不可少的组成部分。构建地形格网，剖分为三角形面元参与辐射计算，可进一步提高光辐射模拟的准确度。

本文将采用辐射度算法模拟顾及地形影响的枇杷果园的光分布。

2 枇杷果园场景构建

2.1 数据获取

构建精细枇杷几何模型是进行冠层光分布模拟的基础。本研究团队于2016年1月在福建省云霄县溪口村（117°7′，24°8′）选取若干枇杷树作为调查对象，树龄为3年。调查与测量的主要内容有：（1）使用杆尺、皮尺等进行枇杷形态结构参数测量，具体参数有树高、叶幕层厚度、干高、冠幅、各级枝条分枝长度、周长、角度、分枝数等；（2）使用Artec Eva手持三维激光扫描仪扫描枇杷叶簇获取叶器官点云数据；（3）使用RIGEL VZ400地面三维激光扫描仪获取枇杷果园地形点云数据。

2.2 地形网格生成

使用VZ400地面三维激光扫描仪获取的地形点云数据密度较大，且分布均匀，适合于反距离权重插值法（Inverse Distance Weighted，IDW）的应用条件[17]。假设未知地形点z0受近距离已知地形点的影响比远距离已知地形点的影响更大，z0的高程为：

式中zi为已知点i的高程；di是点i与点0之间的距离；s为用到的已知点的数量；k是确定的幂，在本例中k选值为2。IDW插值的一个重要特性是所有预测点高程值都位于已知点高程值最大最小值之间，插值后地形点高程不会超出原始地形点云覆盖的高程范围，且分布均匀，便于构造三角格网。地形格网生成过程如图1所示。

图1 地形格网生成过程

2.3 单株木冠层模型构建

利用叶器官点云数据，使用Geomagic Studio与3ds Max软件构建枇杷叶片精细三维模型。根据叶片大小不同，每个叶片约由80～200个三角面元构成。为了模拟枇杷树叶真实的组织状态，由叶簇挂接于枝干系统，形成虚拟枇杷树模型，每个叶簇约由3～6片叶片构成，包含的三角面元数约为500～1 000个，如图2所示。

采用研究团队自主研发的ParaTree参数化单树建模软件[18-19]，输入调查获得的各项枇杷形态结构特征参数，建立枝干系统，在枝系统上挂接叶器官模型，构建精细枇杷单树几何模型，如图3所示，模型总面元数为141 217个。

图2 叶簇构建过程

图3 单树枇杷几何模型

2.4 虚拟果园场景生成

虚拟果园场景抽象为由三维地形和三维冠层模型组成。本研究样地的地形呈阶梯状，果树栽种于每一垄。根据实地调查数据，可以计算每株果树种植的地理坐标。精细单树枇杷几何模型数据结构中包含了组成树模型不同器官（包含枝干、叶片等）的大量三角形面元的索引数组，索引数组对应于坐标数组，实现了一一对应关系。因此在虚拟枇杷果园场景中只需使用实例化技术调用树模型对应的类的构造函数，再结合枇杷树种植点位，即可在不同种植点位“种植”枇杷树模型，如图1（4）所示。此时，树模型的坐标数组并未根据“种植”点位不同而存储多次，仅在绘制阶段调用三角形面元索引，获取对应的坐标值，做相应的矩阵平移操作，节约了大量的内存开销，保证了三维果园场景的流畅显示。

3 枇杷果园冠层内光三维分布模拟

3.1 技术流程与实现

果园的光源主要包括太阳直射光和天空散射光。虚拟枇杷果园冠层内光分布模拟总体技术流程如图4所示。

首先根据研究样地的地理经纬度以及模拟的日期、时刻，利用天文计算公式[20]估算样地所在地的太阳高度角和方位角，进而确定光源位置、光线方向、光源面片初始携带能量等参数。对于太阳直射光源，在光源位置构建光投放平行四边形并剖分为多个三角形面元，模拟平行直射光入射植物冠层。而散射光源则构建一个包围整个果园场景的天空球，将此天空球按照经纬度方向剖分为多个三角形面元，从而模拟均质散射光从多个方向离散射入植物冠层。根据枇杷几何模型与光源几何信息，采用三维离散视角因子方法[21]（3D-Discrete View Factor，3D-DVF）计算辐射度模型形状因子。求解辐射度模型，从而计算冠层内太阳直射光分布、天空散射光分布以及地形网格反射初始能量值。将地形网格作为反射光源，与枇杷几何模型再次构建辐射度模型并进行求解，计算枇杷果园冠层内地面反射光分布。将每个面元获得的直射光、散射光、地面反射光的辐射强度进行累加，即得到枇杷果园冠层内每个面元的辐射强度。

图4 虚拟枇杷果园冠层光分布模拟总体技术流程

3.2 虚拟果园场景组织

实地调查样地枇杷树总量为12棵，按每棵树约由14万个三角形面元构成，场景数据量约为170万个面元。以计算地面反射光在冠层内分布为例，为保证精度，将地形面元剖分为1 250个三角形面元，总计算量超过21亿次以上，系统内存消耗巨大，因此高效组织管理虚拟场景非常重要。为了实现树面元与光源面元间的遮挡判断，采用均匀体素剖分场景。根据场景树模型，建立场景最小包围盒，将包围盒划分为64×64×64均匀体素网格。针对场景中每个树面元做归属判定，若某树面元三角形中心点位于某体素内，则认为该面元与此体素具有唯一相关性。将平面直线的Bresenham算法推广至三维条件下，设某空间直线段两端点坐标为(x1,y1,z1)与 (x2,y2,z2)，令 dx=|x1-x2|，dy=|y1-y2|，dz=|z1-z2|，将该空间直线段投影至xoy平面后即可根据平面直线的Bresenham算法求出x坐标与y坐标；同理可求出投影至xoz平面上的x坐标与z坐标，据此可求出该直线与空间体素相交的任意三维点坐标。

光源面元发射出的光线在包围盒内遮挡判断，光线原点为光源三角面元中心点，方向对应于树面元中心点，除每条光线与包围盒相交的第一个体素外，其余相交体素内的树面元均视为被遮挡。CUDA编程模型提供了基于线程网格与线程块的双层结构。采用一维线程格网与一维的线程块，一个线程内计算单个树面元与单个光源面元间的遮挡情况，一个线程块内计算单个树面元与所有光源面元间的遮挡情况。基于三维体素遍历的并行遮挡判断流程如图5所示。

图5 遮挡判断流程图

3.3 形状因子与辐射度求解

形状因子又称视角因子，是仅与面元在场景中所处位置有关的无量纲常数，面元Si与面元Sj间的形状因子可表示为：

式中HID(dAi,dAj)定义为遮挡系数；θi和θj是面元Si和面元Sj法向量与两面元间夹角；rij为面元Si与面元Sj间的距离。采用3D-DVF法求解辐射度模型，具有普适性高、可并行实现的特点。该方法可表示为：

使用CUDA编程模型实现形状因子与辐射度求解，需要迭代3次式（3），第一次将地形面元视为果园场景一部分，与光源面元进行辐射计算。采用一维线程格网与一维线程块双层结构，一个GPU线程计算一对树面元（包含地形面元）与光源面元间的形状因子能量，一个线程块计算一个树面元（包含地形面元）与所有光源面元间形状因子与能量，据此求出太阳直射辐射强度[22]。同理可求出天空散射辐射强度。将地形面元接收到的直射与散射辐射能量乘以地面反射率，取值为0.20[23]，作为地面反射光源，与所有树面元进行辐射计算，CUDA线程结构与前两次类似。

4 分析与讨论

以2016年1月20日于福建省云霄县溪口村样地实地调查所获数据为依据，枇杷树叶无凋落现象，不考虑叶片光透射的情况。估算冠顶辐射强度时，假设研究样地天气晴朗无云。模拟从早上8时开始，间隔为1 h，选取一天中的10个时刻，依次模拟果园群体冠层的光分布。本案例模拟太阳光在400～700 nm波段的辐射，即为光合有效辐射（Photosynthetically Active Radiation，PAR）分布。

果园场景模拟以三角形面元为最小单位，模拟结果可以输出果园冠层每个叶片的每个三角面元截获的PAR能量值，从而可以计算冠层在某时刻获得的PAR能量，亦可进行三维可视化输出。如图6（1）所示为样地正午12时冠层光分布模拟效果图，在接受太阳直射、天空散射与地面反射光源辐射后，模型中颜色越亮的部分代表到达该面元的光辐射能量越多，PAR强度越强。

根据研究样地地形数据和树木的分布情况，构建12株与16株两种不同种植密度的虚拟枇杷果园，其中12株情况为样地真实果树位置，16株情况为增加种植密度后的对照组。分别研究太阳直射、天空散射、地面反射以及总PAR辐射强度。单位面积叶获得PAR瞬时能量值代表该时刻枇杷果园冠层对光的截获能力，以及枇杷植株对光的利用率，进而可判断该种植密度是否合理。

图6（2）中早8时12株密度下获得平均太阳直射PAR强度高于16株密度，是因为早晨太阳高度角较低，光线从植株侧面进入冠层，植株间遮挡明显，更大的种植密度的情况下会有更多的叶片相互遮挡；在正午时刻，光线从冠顶进入冠层，这两种种植密度的冠层内部叶片相互遮挡情况类似，且3年生枇杷还未郁闭，植株间的遮挡不明显，因此，此时太阳直射瞬时PAR能量接近。图6（3）中12株密度下天空散射瞬时能量值在各个时刻均高于16株密度，且差值相近，是因为不同时刻下散射光源天穹发出的光线方向一致，仅有辐射能量大小不同，枇杷果园冠层内叶片相互遮挡情况几乎一致。图6（4）和图6（5）表明果园冠层叶片获得经过地面反射后的光辐射能量约为地形网格初始反射能量的1/3，且正午时刻最高截获能量速率为35 μmol·m-2·s-1。结合图6（6）分析，同一时刻冠层光截获地面反射约占冠层总光截获的1/20～1/10。两种不同种植密度下的虚拟枇杷果园在一天10个时刻内的单位面积的PAR变化趋势一致，早晨和傍晚太阳辐射强度弱，因此冠层获得的PAR值相应变小。对一天10个时刻的模拟值进行平均，12株密度下果园获得的平均PAR值为236.4 μmol·m-2·s-1，16 株密度下则为 224.8 μmol·m-2·s-1，两者相差约为11.6 μmol·m-2·s-1。模拟结果表明，对于3年生的枇杷植株，在样地的生境下，这两种种植密度均未达到郁闭，株间遮挡对光分布影响不是很大。

图6 果园冠层内光分布分析

5 结束语

目前果树冠层光分布模拟多数停留在单株尺度，针对这种研究现状，同时考虑地面反射对冠层内光截获的影响，设计并实现了利用辐射度的枇杷果园冠层光分布模拟方法。利用地面激光雷达扫描获取的地形点云数据构建果园地形网格，手持激光扫描获取的枇杷叶器官点云数据构建叶簇实体三维模型，然后采用单树建模软件ParaTree建立枇杷冠层模型，最后根据植株的分布规则，形成三维果园场景。使用3D-DVF方法求解辐射度模型，开展不同种植密度下果园冠层光分布模拟，可定量化计算冠层内每个组分获得的光的强度，且可直观表达光在冠层内空间分布情况。模拟结果表明，果园冠层内光强度分布日变化呈近似抛物线，早晚较小而正午到达最大值；果园冠层截获的地面反射PAR能量主要作用于冠层内部与下层，这部分光能对于果园内枇杷植株生长发育的作用不容忽视。

参考文献：

[1]高照全，赵晨霞，李志强，等.我国4种主要苹果树光合能力差异研究[J].中国生态农业学报，2013，21（7）：853-859.

[2]Baldocchi D，Hutchison B.On estimating canopy photosynthesis and stomatal conductance in a deciduous forest with clumped foliage[J].Tree Physiology，1986，2：155-168.

[3]Tang L，Hou C，Huang H，et al.Light interception efficiency analysis based on three-dimensional peach canopy models[J].Ecological Informatics，2015，30：60-67.

[4]Ross J.The radiation regime and architecture of plant stands[J].The Quarterly Review of Biology，1982，71（3）：344.

[5]Gastellu-Etchegorry J P，Demarez V，Pinel V，et al.Modeling radiative transfer in heterogeneous 3-D vegetation canopies[J].Remote Sensing of Environment，1996，58（2）：131-156.

[6]Chelle M，Saint-Jean S.Taking into account the 3D canopy structure to study the physical environment of plants：The Monte Carlo solution[C]//4th International Workshop on Functional-Structural Plant Models.Montpellier，France：CIRAD-AMAP，2004：176-180.

[7]Huang H，Qin W，Liu Q.RAPID：A Radiosity Applicable to Porous IndiviDual objects for directional reflectance over complex vegetated scenes[J].Remote Sensing of Environment，2013，132（10）：221-237.

[8]Silva D D，Han L，Costes E.Light interception efficiency of apple trees：A multiscale computational study based on MAppleT[J].Ecological Modelling，2014，290（10）：45-53.

[9]Wiechers D，Kahlen K，Stützel H.Evaluation of a radiosity based light model for greenhouse cucumber canopies[J].Agricultural and Forest Meteorology，2011，151（7）：906-915.

[10]孔娅，劳彩莲，曹素云.利用3D模型模拟天空与叶面散射对玉米冠层截光率的影响[J].农业工程学报，2011，27（5）：248-252.

[11]黄帆，劳彩莲，肖翠霞.基于光线跟踪的冠层光分布模型参数研究[J].中国农业大学学报，2013，18（6）：96-101.

[12]Kang H，Fiser M，Shi B，et al.IMapple-functional structural model of apple trees[C]//2016 IEEE International Conference on Functional-Structural Plant Growth Modeling，Simulation，Visualization and Applications，China，2016：90-97.

[13]Reinhard H，Ole K，Dirk L，et al.The rule-based language XL and the modelling environment GroIMP illustrated with simulated tree competition[J].Functional Plant Biology，2008，35（10）：739-750.

[14]温维亮，孟军，郭新宇，等.基于辐射照度的作物冠层光分布计算系统设计[J].农业机械学报，2009，40（S1）：190-193.

[15]Buck-Sorlin G，de Visser P H，Henke M，et al.Towards a functional-structural plant model of cut-rose：Simulation of light environment，light absorption，photosynthesis and interference with the plant structure[J].Annals of Botany，2011，108（6）：1121-1134.

[16]王新云，郭艺歌，过志峰，等.地形因子对森林冠层BRF模拟的影响分析[J].遥感信息，2012（4）：82-85.

[17]刘光孟，汪云甲，张海荣，等.空间分析中几种插值方法的比较研究[J].地理信息世界，2011（3）：41-45.

[18]Tang L，Chen C，Zou J，et al.OntoPlant：An integrated virtual plant software package for different scale applications[C]//Spatial Data Mining and Geographical Knowledge Services（ICSDM），2011：308-314.

[19]林定，陈崇成，唐丽玉，等.基于颜色编码的虚拟树木交互式修剪技术及其实现[J].计算机辅助设计与图形学学报，2011，23（11）：1799-1807.

[20]王炳忠.太阳辐射计算讲座：第一讲 太阳能中天文参数的计算[J].太阳能，1999（2）：8-10.

[21]赵权，黄运保，孙宇航.CUDA架构下的靶丸辐射能流并行计算[J].计算机辅助设计与图形学学报，2013，25（7）：937-945.

[22]侯璨，唐丽玉，陈崇成，等.基于并行辐射度的虚拟植物冠层内光分布模拟[J].系统仿真学报，2015（10）：2337-2343.

[23]林正云.福建省太阳总辐射和地面辐射平衡的分布[J].太阳能学报，1994（3）：248-256.