陈舒炜,唐丽玉,尹 丹

(福州大学空间数据挖掘与信息共享教育部重点实验室,地理空间信息技术国家地方联合工程研究中心,福建 福州 350116)

0 引言

果树形态结构决定了冠层内光分布,而光分布影响光合作用效率和冠层内微环境,从而影响果实产量和质量. 目前植物冠层内光环境定量分析方法有地面实测法[1-2]、数学模型法[3-4]和三维结构模型法[5-6]. 地面实测法通过带有传感器的冠层分析仪或三维数字化仪直接测量,但受到天气状况、人力等因素的限制,实际测量的空间和时间分辨率有限. 数学模型法根据辐射传播规律和林冠结构的数学抽象,建立冠层结构参数和生态因子的辐射传输模型,这类模型一般基于实际情况的简化和假设,部分条件在真实植物群落中难以满足. 三维结构模型法以树木空间结构数据为基础,以更精细的面元表达枝干和叶片器官三维信息,建立冠层三维结构模型,模拟任意时刻不同尺度的冠层内光分布情况.

辐射度算法是研究三维植物冠层辐射传输模拟的主要方法之一. 近年来学者们集成了虚拟植物冠层三维模型和辐射度算法,用于模拟冠层内辐射传输过程. Wiechers等[7]从枝叶尺度,利用辐射度方法模拟温室中黄瓜冠层PAR(photosynthetically active radiation)分布. 谢东辉等[8]利用结合真实结构场景的辐射度模型,计算玉米场景的二向性反射率因子. 由于传统辐射度算法中计算形状因子效率较低,Huang等[9]提出RAPID模型模拟复杂场景中的任意植被冠层,减少大型场景计算中巨大的内存需求. 文献[10-11]利用CUDA架构对辐射度算法进行并行计算,提高其计算效率.

植物冠层内辐射仿真模拟是定量分析冠层内光分布的有效新手段[12]. 在分析冠层内光截获效率方面,从叶面积角度考虑,大多数研究引入分析指标STAR(silhouette to total area ratio)量化评价冠层光截获效率[13-15]； 从辐射能量角度考虑,则通常选用冠层内相对光合有效辐射RPAR(relative photosynthetically active radiation)说明冠层光截获能力[16-17].

目前植物冠层光环境模拟和分析大多数集中于株型较小的作物研究,关于虚拟枇杷树的研究报导仍不多见. 枇杷树是亚热带小型常绿乔木,其果实味道甜美、营养丰富,其叶子亦可在中医中入药. 本研究以枇杷树为例,开展冠层内光合有效辐射分布模拟与光截获分析,通过分析与评价枇杷冠层光截获效率,为冠层结构改善和株型优化提供参考依据.

1 虚拟植物精细建模

以国家果树种质福州枇杷圃(26°7′N,119°20′E)为调查样地,于2017年3月1日调查龙早1号枇杷树形态结构参数. 使用皮尺、量角器等工具测量干高、干周、冠幅等主要形态结构参数,统计一级枝干数、二级枝干数、总枝梢数等主要生长特征参数,用于枇杷树枝干建模. 使用Artec Eva手持三维结构光扫描仪扫描枇杷叶簇,获取叶片激光点云数据,用于枇杷叶器官建模.

1.1 基于激光点云数据的叶片器官建模

图1 枇杷叶片器官三维模型Fig.1 Three-dimensional model of loquat leaf organs

利用Artec Studio软件处理叶片点云数据,获取叶片顶点坐标和纹理坐标. 使用Geomagic Studio软件对叶器官模型面片网格简化,减少模拟过程中计算量,同时保持叶片轮廓特征,如图1(a)所示. 为模拟枇杷叶真实的组织排列状态,使用3D Max将多个叶片组合成单簇枇杷叶片三维模型,如图1(b)所示.

1.2 基于L-系统生长规则的树木枝干建模

基于研究团队自主研发的单树建模软件LSTree[18],根据枇杷主要形态结构参数,使用参数化L-系统编写龙早1号枇杷生长规则,构建枝干三维模型,如图2(a)所示. 然后在枝干顶端挂接叶片器官,形成单株枇杷三维精细模型,如图2(b)所示. 三维模型中枝干器官剖分面元数量为12 720个,叶片器官剖分面元数量为390 888个,冠层叶面积约为6.663 m2.

图2 集成枝干和叶片的单株枇杷三维模型Fig.2 Three-dimensional model of single loquat with branches and leaves

1.3 三维模型参数对比

表1给出了枇杷树实际测量的树高、胸径、冠幅、叶幕层厚度等主要参数和三维枇杷模型的相应参数. 从绝对误差和相对误差分析,三维模型的参数与真实值比较接近,能表达冠层形态特征.

表1 三维枇杷模型与真实枇杷果树的主要形态结构参数对比

2 冠层内光分布模拟

图3 冠层内光分布模拟技术流程图 Fig.3 Flow diagram of light distribution simulation in canopy

根据调查样地的地理坐标以及调查日期,利用天文公式计算太阳高度角和太阳方位角[19],从而确定光源位置、入射光方向和入射光强度等必要的环境模拟参数. 在虚拟场景中,构成枝干、叶片、光源的三角形面元是最基本元素. 采用辐射度算法模拟冠层内光分布,模拟主要过程为： 首先加载枇杷三维模型,构建场景中直射光源； 然后构建最小八叉树包围盒并均匀剖分成若干子立方体空间,判定所有三角形面元归属； 最后求解辐射度方程,利用CUDA并行加速技术实现辐射模拟. 具体技术流程见图3.

2.1 计算形状因子

形状因子是求解辐射度模型的关键. 在虚拟场景中,光源与树模型视作由若干个三角形面元组成. 假设不存在遮挡的理想情况下,树模型面元j接受来自光源面元i的辐射能量,其数学表达式为：

Φj=BiFijSi

(1)

其中：Φj为树面元j接受来自光源的辐射能量；Bi为离开光源面元i的总辐射能量；Si为光源面元i的面积；Fij为光源面元i与树面元j之间的形状因子,表示离开面元i而到达面元j的能量占离开面元i的能量的百分比,形状因子的数学表达式为

(2)

其中：θi和θj分别为面元i和面元j的法向量与两面元中心点连线之间的夹角;rij为面元中心点之间的距离; HID为遮挡函数,当面元i与面元j之间存在遮挡时,其值为0,无遮掩时为1.

2.2 光线遮挡判断

实际计算所有面元的形状因子时,必须考虑三维场景中存在大量面元会被遮挡,无法接受辐射能量的情况. 当某个树面元j在以面元i为光源的入射光线方向上被其他树面元遮挡时,只需找到光线路径上第一个与光线相交的面元进行形状因子计算,排除其他面元多余的计算.

采用射线求交方法,对所有树面元和光源面元进行两层遍历. 以光源面元中心点为起点,以树面元中心点与起点连线为入射光方向进行求交计算,得到与该条入射光线有交点的所有树面元,除距离光源面元最近的面元外,其余面元均视为被遮挡. 依次判断每一对树面元和光源面元之间的遮挡情况.

3 冠层内PAR分布

图4 冠层内PAR强度与太阳高度角日变化情况Fig.4 Diurnal change trend of PAR intensity and solar elevation

晴朗无云的天气下,2017年3月1日7:00至17:00各整点时刻冠层内直射PAR、散射PAR和太阳总辐射PAR模拟结果如图4所示. 冠层内直射PAR强度从7:00逐渐增大,于12:00到达最高点后逐渐减小,全天变化趋势呈单峰型曲线,峰值为355.119 μmol·(m2·s)-1； 冠层内散射PAR强度全天变化平稳,日均值为31.429 μmol·(m2·s)-1； 太阳总辐射PAR与直射PAR变化趋势十分接近,日均值为267.179 μmol·(m2·s)-1. 直射PAR是冠层内总PAR的主要组成部分,散射PAR占比仅为11.76%.

太阳高度角与冠层内直射PAR强度变化趋势相似,随时间呈先升后降的单峰型曲线,峰值为56.1°,出现在12:00,此时太阳高度角最大,冠层内直射PAR强度最大. 冠层内直射PAR强度与太阳高度角随时间呈正相关,散射PAR强度与太阳高度角无明显对应关系.

利用OpenGL渲染技术绘制枇杷三维模型模拟结果,图5展示了不同时刻枇杷冠层内太阳直射PAR空间分布情况. 图中假彩色表示树模型面元接收的太阳直射PAR强度,每个面元接收的PAR能量越多,其渲染的颜色越亮,反之颜色越暗表明面元受遮挡情况越严重.

图5 枇杷冠层内不同时刻太阳直射PAR分布情况Fig.5 Distribution ofsolar direct PAR within a loquat canopy at several o’clock

图6 PAR实测值与模拟值相关性分析Fig.6 Correlation between measured values and simulated values of PAR

太阳光源位置随时间推移由东至西,早晨太阳高度角较低,入射光从枇杷树东侧面进入冠层,冠层东部外围及上层区域接收光照较充分,树的背阴面则明显受遮挡严重； 正午时分,入射光从枇杷树侧上方进入冠层,整个冠层上层区域及外围中低层区域能接收光,冠层内膛及中低层区域受到上方枝叶遮挡； 傍晚太阳高度角变小,入射光从枇杷树西侧面进入冠层,冠层西部外围及上层区域接收光照较充分,树的背阴面则遮挡严重,但因为真实的冠层结构不对称,东西两侧枝叶数量和密度均相差较大,东部区域枝叶繁茂而西部区域枝叶较稀疏,所以傍晚入射光穿透冠层的情况好于早晨. 随时间推移,冠层内光主要分布在向阳面的外围及中上层区域.

实地测量过程中,随机选择多个时间点,使用LI-COR太阳辐射传感器采集冠顶与冠层内不同部位的瞬时辐射能量数据,对比计算机程序模拟的瞬时辐射能量数据,验证模拟计算的准确程度. 利用Origin软件对PAR实测值和模拟值进行线性拟合分析,结果如图6所示. 样本数量n为16,线性回归方程y=1.066 2x-6.083 1,决定系数R2为0.968 6,证明实测值与模拟值相关性很高； 均方根误差RMSE为29.527 μmol·(m2·s)-1,误差在可接受范围内,表明模拟值与实测值较吻合.

4 冠层光截获分析

冠层内光截获分析基于冠层内直射PAR模拟结果. 假设将整个冠层空间均匀划分为等大小的若干个0.5 m×0.5 m×0.5 m的小立方格区域,将所有叶器官三角面元根据空间位置归属至各个小立方格中,并将其叶面积、PAR能量等属性视为该立方格的属性. 以小立方格作为基本单位进行不同时刻的辐射模拟,计算得到每个小立方格的属性值,用于分析整个冠层光截获情况.

4.1 冠层内STAR分布

为评价植物冠层光截获效率,从叶面积受光情况考虑,通常引用STAR值来表示[14]. STAR的含义是冠层内叶片被光线照射到的叶面积(PLA)与冠层总叶面积比值(TLA)[15],其数学表达式为：

(3)

结合辐射模拟结果与冠层空间分格,可得到STAR 值空间分布. 表2是以STAR值为评价标准的不同时刻冠层光截获效率情况,其中12:00光截获效率最高,8:00光截获效率最低,下午光截获效率高于上午,冠层内STAR日平均值为0.601. 当某区域STAR值较低时,表明该区域受遮挡情况严重,光截获能力弱. 以区域STAR<0.3为低光截获区的衡量指标,分析冠层内光截获效率.

表2 基于STAR值的不同时刻冠层光截获效率情况

图7(a)为8:00、11:00、14:00、17:00时刻的冠层内STAR<0.3的区域分布情况,图中紫色立方格表示低光截获区,红色箭头表示入射光大致方向. 上午太阳直射光从图中右上角射入,8:00时太阳高度角为19.5°,直射光从东侧面进入冠层,低光截获区主要分布在树背阴面,上、中、下层区域均有； 11:00时太阳高度角为51.6°,直射光从东侧上方进入冠层,低光截获区数量减少,分布在树背阴面中下层区域. 下午太阳直射光从图中左上角射入,14:00时太阳高度角为47.8°,直射光从西侧上方进入冠层,低光截获区较少,零散分布在树背阴面中下层区域； 17:00时太阳高度角为13.4°,直射光从西侧面进入冠层,此时低光截获区增加,集中分布在树背阴面中下层区域.

4.2 冠层内RPAR分布

从辐射能量角度考虑,使用相对光照强度RPAR值,即冠层内某处的光照强度与冠顶光照强度的比值,分析冠层内光截获效率. 当冠层内叶片接受RPAR低于12%时[16],全天净光合总量为0,这类叶片可称为无效叶片,其所在区域为无效光区. 无效光区占冠层空间的比例反映了冠层内光截获情况(见图7(b)).

图7 冠层内不同时刻光截获空间分布情况Fig.7 Spatial distribution of light interception within canopy at several o’clock

图7(b)为8:00、11:00、14:00、17:00时刻的冠层内RPAR<0.12的区域分布情况,其中墨绿色立方格表示RPAR<12%的区域,红色箭头表示入射光大致方向. 需要说明一点： 图7(b)各时刻显示的无效光区数量均多于图7(a)显示的低光截获区,是因为RPAR作为一个相对比值,RPAR<0.12仅供参考, 而不是绝对的衡量标准,其曲线变化趋势更有价值,两者不能从绝对数量上体现差异.

对比图7(b)中各时刻RPAR分布情况,可以发现17:00时冠层内无效光区明显偏少,这与从叶面积角度的无效光区分布随时间变化情况存在区别. 结合真实冠层情况分析,其原因在于冠层西部区域枝叶稀疏,数量明显少于枝叶繁茂的冠层东部区域.

4.3 光截获综合分析

综合STAR与RPAR两项分析指标,分别以STAR<0.3和RPAR<0.12作为低光截获区和无效光区的衡量标准,分析冠层内光截获日变化情况,结果见图8. 图8中STAR平均值曲线最低点出现在8:00,最高点在12:00,表明8:00时冠层内光截获效率最低,而12:00时光截获效率最高. 其变化趋势虽非单峰型曲线,但大致上同太阳高度角变化相似,STAR值随太阳高度角增大而上升,随太阳高度角减小而下降.

无效光区比例曲线最高点出现在8:00,说明此时冠层内无效光区数量最多,光截获效率最低. 随后曲线逐渐下降并保持较稳定状态,至14:00后骤降,于17:00降至最低点. 其中10:00至14:00太阳高度角保持较大值,太阳光源位置长时间处于冠层侧上方位置,冠顶直射PAR强度较大,因此冠层内光截获效率上升,无效光区减少； 而15:00后,虽太阳高度角减小,但冠层西部区域枝叶稀疏,入射光易进入冠层,冠层光截获效率保持较高水平,无效光区较少. 冠层内无效光区变化较大,其原因在于真实冠层结构的复杂性和不对称性,东部区域枝叶繁茂而西部区域枝叶较稀疏.

图8 冠层内光截获能力日变化情况Fig.8 Diurnal change trend of light interception ability within canopy

试验结果表明,该枇杷冠形光截获效率的时空分布具有如下特征： ① 时间角度. 8:00时冠层内无效光区比例最高且STAR平均值最小,此时冠层内遮挡最严重,光截获效率最低； 12:00时因太阳高度角最大,冠顶直射PAR强度最强,冠层STAR平均值最大,光截获效率较高； 因冠层西部区域枝叶稀疏,14:00之后冠层内无效光区比例相对较小,冠层截获效率较高. ② 空间角度. 从冠层顶部至底部,从冠层外围到内膛,随着枝叶遮挡情况逐渐严重,获得光照逐渐减少,无效光区集中在树体的背阴面、冠层的内膛及中下层区域.

5 结论

1) 基于激光点云数据的叶片器官建模方法与基于L-系统生长规则的树木枝干建模方法相结合所建立精细的三维枇杷模型,与真实枇杷果树对比,其主要参数精度较高,能够表达枇杷冠层形态特征.

2) 冠层内直射PAR强度全天变化趋势呈先升后降的单峰型曲线,峰值355.119 μmol·(m2·s)-1出现在12:00,与太阳高度角变化趋势成正相关. 冠层内光主要分布在向阳面的外围及冠层中上层区域.

3) 从叶面积受光情况的角度和辐射能量角度综合分析模拟结果,总结冠层内光截获日变化情况. 冠层内STAR日平均值为0.601,无效光区日平均比例为34.8%,8:00时冠层内光截获效率最低. 冠层内光截获效率的变化,是太阳高度角、入射光强度和枇杷真实冠型结构共同作用的结果.

4) 探索冠层内光分布规律,考虑冠层STAR和RPAR分布变化情况,可用于评价冠层光截获效率,寻找无效光区,为冠型修剪、改善冠层光照微环境提供参考. 通过改变冠层结构、减少总枝量、降低冠层内枝叶相互遮挡程度,可改善光照分布和光能有效利用情况,提高果树产量和品质.

致谢： 感谢福建省农业科学院果树研究所蒋际谋研究员、邓朝军老师在枇杷形态结构调查过程的支持、帮助和指导,感谢研究团队其他成员参与野外调研和数据处理工作.