CN114663786B - 基于点云数据和计算机图形学的林分辐射通量计算方法 - Google Patents
基于点云数据和计算机图形学的林分辐射通量计算方法 Download PDFInfo
- Publication number
- CN114663786B CN114663786B CN202210276654.1A CN202210276654A CN114663786B CN 114663786 B CN114663786 B CN 114663786B CN 202210276654 A CN202210276654 A CN 202210276654A CN 114663786 B CN114663786 B CN 114663786B
- Authority
- CN
- China
- Prior art keywords
- solar
- stand
- forest
- triangle
- incident
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000004907 flux Effects 0.000 title claims abstract description 134
- 230000005855 radiation Effects 0.000 title claims abstract description 127
- 238000000034 method Methods 0.000 title claims description 104
- 241000218631 Coniferophyta Species 0.000 claims abstract description 59
- 238000004364 calculation method Methods 0.000 claims abstract description 32
- 238000004088 simulation Methods 0.000 claims abstract description 12
- 239000013598 vector Substances 0.000 claims description 35
- 230000005540 biological transmission Effects 0.000 claims description 23
- 238000005070 sampling Methods 0.000 claims description 22
- 238000004422 calculation algorithm Methods 0.000 claims description 19
- 238000002834 transmittance Methods 0.000 claims description 16
- 239000000443 aerosol Substances 0.000 claims description 15
- 241000894007 species Species 0.000 claims description 11
- 238000005286 illumination Methods 0.000 claims description 10
- 230000003993 interaction Effects 0.000 claims description 9
- 238000002310 reflectometry Methods 0.000 claims description 7
- 238000010521 absorption reaction Methods 0.000 claims description 5
- 230000011218 segmentation Effects 0.000 claims description 5
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 claims description 3
- 241000196324 Embryophyta Species 0.000 abstract description 22
- 210000000056 organ Anatomy 0.000 abstract description 9
- 238000012512 characterization method Methods 0.000 abstract description 4
- 238000009826 distribution Methods 0.000 description 30
- 238000010586 diagram Methods 0.000 description 20
- 238000005259 measurement Methods 0.000 description 18
- 230000000694 effects Effects 0.000 description 7
- 238000011160 research Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 5
- 230000006870 function Effects 0.000 description 5
- 238000011835 investigation Methods 0.000 description 5
- 239000000203 mixture Substances 0.000 description 5
- 230000007423 decrease Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000002329 infrared spectrum Methods 0.000 description 3
- 230000003287 optical effect Effects 0.000 description 3
- 230000035790 physiological processes and functions Effects 0.000 description 3
- 230000003595 spectral effect Effects 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 241000218645 Cedrus Species 0.000 description 2
- 241000723346 Cinnamomum camphora Species 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 239000000428 dust Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000011049 filling Methods 0.000 description 2
- 239000007789 gas Substances 0.000 description 2
- 239000011121 hardwood Substances 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000029553 photosynthesis Effects 0.000 description 2
- 238000010672 photosynthesis Methods 0.000 description 2
- 238000009877 rendering Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- DSSYKIVIOFKYAU-XCBNKYQSSA-N (R)-camphor Chemical compound C1C[C@@]2(C)C(=O)C[C@@H]1C2(C)C DSSYKIVIOFKYAU-XCBNKYQSSA-N 0.000 description 1
- 241000251468 Actinopterygii Species 0.000 description 1
- 235000011201 Ginkgo Nutrition 0.000 description 1
- 244000194101 Ginkgo biloba Species 0.000 description 1
- 235000008100 Ginkgo biloba Nutrition 0.000 description 1
- 241000218300 Liriodendron Species 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 244000242564 Osmanthus fragrans Species 0.000 description 1
- 235000019083 Osmanthus fragrans Nutrition 0.000 description 1
- 241000218657 Picea Species 0.000 description 1
- 241000219000 Populus Species 0.000 description 1
- 244000141698 Prunus lannesiana Species 0.000 description 1
- 235000014001 Prunus serrulata Nutrition 0.000 description 1
- 240000003377 Shepherdia canadensis Species 0.000 description 1
- 235000018324 Shepherdia canadensis Nutrition 0.000 description 1
- XUIMIQQOPSSXEZ-UHFFFAOYSA-N Silicon Chemical compound [Si] XUIMIQQOPSSXEZ-UHFFFAOYSA-N 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000008033 biological extinction Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000009395 breeding Methods 0.000 description 1
- 230000001488 breeding effect Effects 0.000 description 1
- 229960000846 camphor Drugs 0.000 description 1
- 229930008380 camphor Natural products 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000001678 irradiating effect Effects 0.000 description 1
- 230000031700 light absorption Effects 0.000 description 1
- 238000007726 management method Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000002156 mixing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000877 morphologic effect Effects 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000243 photosynthetic effect Effects 0.000 description 1
- 230000029264 phototaxis Effects 0.000 description 1
- 230000010399 physical interaction Effects 0.000 description 1
- 230000004962 physiological condition Effects 0.000 description 1
- 230000001863 plant nutrition Effects 0.000 description 1
- 238000001556 precipitation Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000001902 propagating effect Effects 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 229910052710 silicon Inorganic materials 0.000 description 1
- 239000010703 silicon Substances 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000003892 spreading Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000002195 synergetic effect Effects 0.000 description 1
- 230000005068 transpiration Effects 0.000 description 1
- 238000009827 uniform distribution Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 238000001429 visible spectrum Methods 0.000 description 1
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于点云数据和计算机图形学的林分辐射通量计算方法，包括：采集林分点云数据；计算林分水平表面的太阳辐照度；单个树冠分割；利用空间半椭球几何基元模拟阔叶树冠或/和利用空间圆锥几何基元模拟针叶树冠；对几何基元模拟后的树冠进行三角剖分；确定入射太阳光束与三角形的第一个交点，计算林分的入射辐射通量；计算林分中树木间的反射辐射通量；计算林分的透射辐射通量。本发明通过确定覆盖森林冠层表面的多个空间三角形中太阳光束与可量化能量的发射和碰撞，克服了与树木器官精细表征和植物系统中因器官尺度异质性而导致的阳光遮挡相关的并发症。
Description
技术领域
本发明涉及林分辐射通量研究技术领域，具体是一种基于点云数据和计算机图形学的林分辐射通量计算方法。
背景技术
太阳辐射是植物光合作用和蒸腾等生理过程的重要能量来源。植物冠层接收到的太阳辐射是动态的，不仅受入射太阳强度和方向的影响，而且还受冠层结构的影响。冠层辐射通量的获取是森林表型研究的前沿领域，对于高光效率树型的选择和育种以及指导生产性林分的森林管理具有重要意义。太阳辐射状态也影响净初级生产力、森林物种丰度、生态系统响应机制和地点适应特征。
事先了解水平表面接收到的本地和区域太阳辐射对于正确测量目标森林冠层的辐射是至关重要的。一些太阳辐射模型，如ASHRAE算法和Iqbal模型，可以根据云的类型和覆盖、大气浑浊度、可感知的含水量和日照时数(SD)等许多易于测量的量来计算相应的结果。这些模型被工程界和建筑界广泛应用于预测日照时数下的日均全球辐射。这些模型中的一些经验函数，如third-degree和logarithmic，与现有的太阳辐射计传感器相结合，测量特定地点的太阳辐射通量密度。此外，准确表征冠层光廓线或冠层辐射传输模型对于预测生态动力学和森林生产力具有重要意义，已引起各研究领域的广泛关注。因此，在现代光学成像技术和先进的计算机技术的帮助下，许多计算太阳辐射通过冠层传输的方法被创造出来，可分为两大类:基于半球照片(hemispherical photograph，HP)的方法和基于计算机图形学的方法。
基于HP的方法使用水平向上的鱼眼镜头拍摄森林照片等距投影，它记录了植物营养分布和天空可见通过获得太阳辐射传输树冠洞或者裂缝差距分数或阻塞的叶面积指数(LAI)的光。各种商业和非商业软件程序已经可以用于HP处理和分析，并结合现有的太阳辐射方程来计算照片中的辐射量。将基于HP的方法嵌入到仪器测量中是一种相对直接的方法，通过测量冠层顶部(或没有冠层的开放区域)与冠层底部的辐射差，可用于评估太阳活动和评估拦截辐射量。相比之下，基于计算机图形学的方法模拟植物冠层内部的辐射传递，并采用某种假设或简化来近似冠层结构；这些技术加上射线追踪方法来模拟单个光子的路径集群和比尔-朗伯定律计算辐射衰减的定量评估赖和树冠的消光系数。表1是对现有方法的简要描述。
表1、现有冠层太阳辐射估计方法简介：
此外，利用基于HP的方法估算冠层太阳辐射对仪器和周围环境条件的性能、精度和测量过程提出了很高的要求。例如，相机需要在森林树冠下绝对水平，这不可避免地会导致难以控制的相机失准问题。此外，在冠层内获得的结果会受到部分阴影影响的意外偏差、由于相机设备参数的相机曝光值变化以及空气中存在灰尘或水蒸气而导致的不确定分数计算，从而将噪声引入图像。对于大型林分，林地样地的高时空采样分辨率加上对森林接收的太阳辐射进行量化的统计分析降低了结果的普遍性，实验工作耗时且容易产生不稳定的结果。相比之下，如上所述，基于计算机图形的方法采用光线追踪方法来模拟太阳光束在树冠元件上的透射、反射和入射；然而，这些方法必然需要更多的光线来充分采样树冠中的每个小植物元素。此外，基于计算机图形的方法的计算复杂性已被证明与异质树冠模型的几何细节水平和目标场景中发射的太阳光线数量呈指数相关。因此，在详细器官水平模拟植物系统的生理过程时，必须在树冠的表型细节水平和计算成本之间进行权衡。
发明内容
本发明所要解决的技术问题是针对上述现有技术的不足提供一种基于点云数据和计算机图形学的林分辐射通量计算方法，本基于点云数据和计算机图形学的林分辐射通量计算方法，通过确定覆盖森林冠层表面的多个空间三角形中太阳光束与可量化能量的发射和碰撞，克服了与树木器官精细表征和植物系统中因器官尺度异质性而导致的阳光遮挡相关的并发症。
为实现上述技术目的，本发明采取的技术方案为：
基于点云数据和计算机图形学的林分辐射通量计算方法，包括：
步骤1：通过激光雷达采集林分点云数据；
步骤2：计算林分水平表面的太阳辐照度；
步骤3：对采集的林分点云数据进行噪声和离群值的去除，然后，使用双高斯滤波器从林分点云数据中检测每个树冠，并基于能量函数最小化方法进行单个树冠分割；
步骤4：将林分内的树木种类分为阔叶树或/和针叶树，利用空间半椭球几何基元模拟阔叶树冠或/和利用空间圆锥几何基元模拟针叶树冠；
步骤5：对几何基元模拟后的树冠进行三角剖分；
步骤6：每个源点发射一束太阳光束，计算已知方向和源点的模拟太阳光束，并计算每束发射的太阳光束的辐射通量；
步骤7：利用三角剖分后生成的三个顶点对应的三角形和模拟太阳光束，采用光束相交算法，确定几何基元模拟后的树冠在三角形尺度下的局部光照和阴影区域；
步骤8：确定入射太阳光束与三角形的第一个交点，进而计算林分的入射辐射通量；
步骤9：计算入射太阳光被林分冠层不同三角形反射后的反射矢量，通过反射矢量计算林分中树木间的反射辐射通量；
步骤10：计算入射太阳光被林分冠层透射后的透射矢量，通过透射矢量计算林分的透射辐射通量。
作为本发明进一步改进的技术方案，所述的步骤2具体为：
林分水平表面的太阳辐照度Itotal为：
Itotal＝Idirect+Idiffuse (1)；
其中Idirect＝Iscr0τrτoτgτwτasinθa；Isc为太阳常数，取1367W/m2；r0为地球轨道偏心率相关系数；τr、τo、τg、τw和τa分别是瑞利、臭氧、气体、水和气溶胶散射透过率；
太阳高度角的表达式如下：
sinθa＝sinθlsinθd+cosθlcosθdcosθh (2)；
其中θl为当地纬度，θd为当前太阳的赤纬角，θh为当地太阳时的时角；
到达水平面的漫射辐照度Idiffuse为：
Idiffuse＝Idr+Ida+Idm (3)；
其中：
τaa为气溶胶吸收引起的直接辐射的透射率，ωo为气溶胶散射到总衰减的入射能量的单次散射反照率，取0.9，ma是实际压力下的气团；
Fc为前向散射相对于总散射的分数，取0.84，τas为气溶胶散射后透射入射能量的分数，计算公式如下：
ρg为地面反照率，ρa为无云天空的反照率，计算公式如下：
ρa＝0.0685+(1-Fc)(1-τas) (8)。
作为本发明进一步改进的技术方案，所述的步骤4中利用空间半椭球几何基元模拟阔叶树冠的空间半椭球方程为：
利用空间圆锥几何基元模拟针叶树冠的空间圆锥方程为：
其中a为每棵单树在东西方向上的半冠宽度；b为每棵单树在南北方向上的半冠宽度；c是树高减去枝下高；被分割树的每个树冠中心坐标记录为：(xcenter,ycenter,zcenter)，zcenter为相应树的枝下高；heighttree为相应树高。
作为本发明进一步改进的技术方案，所述的步骤5具体为：
将代表每个树冠的圆锥体或/和半椭球几何面从X-Y平面上的投影圆进行垂直投影，并在每个投影圆内定义多个等距的同心圆；在每个同心圆上均匀取相同点数，得到X-Y平面上的采样点；利用二维Delaunay三角剖分算法，基于采样点生成一个三角形网格，得到每个采样点的索引，定义构成树冠表面的每个三角形的顶点；根据采样点结果，对每个树冠表面进行反投影，实现三维Delaunay三角剖分。
作为本发明进一步改进的技术方案，所述的步骤6具体为：
定义源点的坐标为其中N表示源点总数，每个源点发射一束太阳光束；相邻两源点之间的固定间隔设为0.2m；即发射的太阳光束密度为
每个太阳光束从每个源点出发的方向表示为单位矢量
已知方向和源点的模拟太阳光束Lk为：
由于发射的太阳光束密度为则每束发射的太阳光束的辐射通量为
作为本发明进一步改进的技术方案，所述的步骤7具体为：
每个三角形Tj,j＝1,2,...M的法向量可以由其三个顶点来计算：
其中×表示两个向量的外积；结合公式(11)，采用公式(13)确定射线Lk是否与三角形Tj所在平面相交：
如果F＝0，发射的太阳光线平行于三角形的平面，且不与该平面相交；
如果F≠0，交点存在，交点坐标为：
计算交点后，采用以下准则判断交点是否在三角形Tj内：
其中，N1,N2和N3由式(16)导出：
通过以上计算，进而确定每个太阳光束与剖分森林冠层表面的每个三角形的交点。
作为本发明进一步改进的技术方案，所述的步骤8具体为：
假设三角形Tj与所有入射的太阳光束的第一个交点总数为即入射的太阳光束的总数首先被三角形拦截，那么，林分的入射辐射通量或入射辐射功率表示为：
其中为入射太阳光束Ld与被照射三角形Tj法向量的夹角。
作为本发明进一步改进的技术方案，所述的步骤9具体为：
入射太阳光被森林冠层不同三角形反射后的反射矢量为：
假设三角形Tj与所有反射的太阳光束的第一个交点总数为反射太阳光束以给定的反射角继续传播，并与森林冠层的其他三角形Tj进行第二次相互作用，那么，林分中树木间的反射辐射通量表示为：
τreflect是三角面所属不同树种的树冠反射率；为反射的太阳光束与三角形法向量的夹角。
作为本发明进一步改进的技术方案，所述的步骤10具体为：
入射太阳光被林分冠层透射后的透射矢量为：
其中e为折射率，取0.5；
林分的透射辐射通量计算公式如下：
其中为所有透射太阳光束与三角形Tj的第一次交点总数，τtrans为各树种的冠层透光率，是透射太阳光束与被光束照射的三角形法向量之间的夹角。
本发明的有益效果为：
本发明的方法不需要大量的射线来充分取样小的植物元素，并结合激光雷达数据导出的树种冠几何形状变化的树木生长特性，实现了小区级和种级辐射通量之间的平衡。此外，考虑到冠层反射率、冠层透光率、冠层形状和冠层大小都是可调节的，并利用灵活数量的发射太阳光束和覆盖森林冠层的细分三角形，根据太阳方位角和高度角，采用定向射线跟踪算法估计碰撞和再碰撞概率，模拟太阳光束在林冠间的入射、反射和透射。因此，本发明的方法具有更强的实用性和可扩展性，适用于任何林分类型和复杂、结构多样的冠层，因此可以渲染植被样地场景所接收到的辐射。
本发明的主要参数的具体设置如下，在树冠建模过程中，我们利用机载激光雷达数据导出的生长特性，根据标准形状的针叶树为圆锥，阔叶树为半椭球，对森林冠层组成进行建模。因此，利用三角剖分策略将森林冠层在三角尺度上细分为不同的详细区域。根据太阳辐照度计算出集合能量的每束离散太阳光束沿当前太阳高度角和天顶角方向发射，然后与代表冠面的三角形相互作用，克服了与树木器官精细表征和植物系统中因器官尺度异质性而导致的阳光遮挡相关的并发症。
附图说明
图1中(a)为研究场地的位置示意图。
图1中(b)为研究场地的放大图。
图1中(c)为实验场地1、2、3生长的林分的对应照片。
图2中(a)为用于模拟针叶树的锥形模型示意图。
图2中(b)为用于模拟阔叶树的半椭球模型示意图。
图2中(c)为针叶林的机载LiDAR数据。
图2中(d)为针叶林经过计算机建模的3D森林样地场景。
图2中(e)为阔叶林的机载LiDAR数据。
图2中(f)为阔叶林经过计算机建模的3D森林样地场景。
图2中(g)为混交林的机载LiDAR数据。
图2中(h)为混交林经过计算机建模的3D森林样地场景。
图3中(a)为采样点分布在每个树冠上并投影到X-Y平面上以优化三角剖分示意图。
图3中(b)为针叶林模拟几何基元的三角剖分结果图。
图3中(c)为阔叶林模拟几何基元的三角剖分结果图。
图3中(d)为混交林模拟几何基元的三角剖分结果图。
图4中(a)为入射太阳光束照射示意图。
图4中(b)为入射太阳光束被反射示意图。
图4中(c)为入射太阳光束透射示意图。
图4中(d)为入射光束被抽取以描绘太阳光在森林地块之间的传播示意图。
图4中(e)为反射光束被抽取以描绘太阳光在森林地块之间的传播示意图。
图4中(f)为透射光束被抽取以描绘太阳光在森林地块之间的传播示意图。
图5为基于HP的HPEval软件对三个森林地块进行短波辐射测量的HP图像。
图6中(a)为研究地点在不同日期和一天中不同时间的总太阳辐照度的计算结果示意图。
图6中(b)为研究地点在不同日期和一天中不同时间的直接太阳辐照度的计算结果示意图。
图6中(c)为研究地点在不同日期和一天中不同时间的漫射太阳辐照度的计算结果示意图。
图7中(a)为2019年7月15日11:00针叶林样地冠层入射辐射通量分布示意图。
图7中(b)为2019年7月15日11:00阔叶林样地冠层入射辐射通量分布示意图。
图7中(c)为2019年7月15日11:00混交林样地冠层入射辐射通量分布示意图。
图7中(d)为2019年7月15日11:00针叶林样地冠层反射辐射通量分布示意图。
图7中(e)为2019年7月15日11:00阔叶林样地冠层反射辐射通量分布示意图。
图7中(f)为2019年7月15日11:00混交林样地冠层反射辐射通量分布示意图。
图7中(g)为2019年7月15日11:00针叶林样地冠层透射辐射通量分布示意图。
图7中(h)为2019年7月15日11:00阔叶林样地冠层透射辐射通量分布示意图。
图7中(i)为2019年7月15日11:00混交林样地冠层透射辐射通量分布示意图。
图8中(a)为采用本发明方法计算的针叶林样地的入射辐射通量示意图。
图8中(b)为采用本发明方法计算的阔叶林样地的入射辐射通量示意图。
图8中(c)为采用本发明方法计算的混交林样地的入射辐射通量示意图。
图8中(d)为采用本发明方法计算的针叶林样地的反射辐射通量示意图。
图8中(e)为采用本发明方法计算的阔叶林样地的反射辐射通量示意图。
图8中(f)为采用本发明方法计算的混交林样地的反射辐射通量示意图。
图8中(g)为采用本发明方法计算的针叶林样地的透射辐射通量示意图。
图8中(h)为采用本发明方法计算的阔叶林样地的透射辐射通量示意图。
图8中(i)为采用本发明方法计算的混交林样地的透射辐射通量示意图。
图9中(a)为正午时刻林分阴影示意图。
图9中(b)为正午太阳入射光束视角下的针叶林分图。
图9中(c)为正午太阳入射光束视角下的阔叶林分图。
图9中(d)为正午太阳入射光束视角下的混交林分图。
图10为显示程序运行时间和相关性的双y轴图。
具体实施方式
下面根据附图对本发明的具体实施方式作出进一步说明：
森林环境具有高度复杂的辐射状态，3D冠层结构阻碍了太阳辐射进入森林；因此，冠层辐射随时空、季节和气象而变化，使得辐射通量的测量和建模都具有挑战性。在这里，本实施例提供了一种协同方法，即一种基于点云数据和计算机图形学的林分辐射通量计算方法，使用机载LiDAR数据和计算机图形学来模拟森林冠层并计算三个林分(针叶林、阔叶林和混交林)的辐射通量。根据太阳高度角和方位角发射定向入射太阳光束，将森林冠层表面分解为三角形单元。使用光线追踪算法来模拟反射和透射光束在森林树冠内的传播。本实施例的方法准确地模拟了太阳辐射通量，并与基于半球照片的HPEval软件和日射强度计测量的结果表现出良好的一致性(R2≥0.82)。由于最大的太阳高度角(81.21°)和最密集的树冠，6月15日中午针叶林达到最大入射辐射通量；由于针叶树更高的反射率和对反射太阳光束更好的吸收，针叶林也达到了最大反射辐射通量(10.91-324.65kW)。然而，由于阔叶树种的透射率较高，阔叶林接收到更多透射辐射通量(46.42-280.93kW)。本实施例的方法可以直接模拟冠层辐射的详细林分尺度分布，对于研究光依赖的生物生理过程很有价值。
1.1、研究内容和技术路线：
机载光探测和测距(LiDAR)通过获取被称为点云的高分辨率表示，提供了森林的自上而下视角。该技术可以在单个树冠分割算法的帮助下揭示树冠的空间结构。因此，可以导出树冠的许多关键属性，如树冠宽度、树高、枝下高和位置等，从而利用规则的几何基元将森林图简化为一组伪浑浊的树冠介质。因此，在计算机图形技术的帮助下，例如树冠表面的三角剖分、光线追踪算法和反射和折射的太阳光束的物理相互作用，本实施例的主要贡献是提出了一种计算森林冠层辐射通量的方法，通过确定覆盖森林冠层表面的多个空间三角形中太阳光束与可量化能量的发射和碰撞。这种方法克服了与树木器官精细表征和植物系统中因器官尺度异质性而导致的阳光遮挡相关的并发症。最后，利用所得到的定量结果，将林分接收到的辐射分布直观可视化，并通过现有的基于HP的方法和野外仪器测量进行了验证。
2.1、研究区和激光雷达数据：
研究场地位于江苏省东南部南京市(118°49′24.13″N,32°5′7.11″E)的南京林业大学校园内(如图1中(a)所示)。南京是一个温暖湿润的亚热带季风气候，四季分明，雨量充沛，属国家森林城市。年平均降水量1000-1600毫米，年平均气温15.4℃。南京林业大学占地面积0.838平方千米，最大海拔4米。该地区的乔木种群由杉木、水杉、雪松、银杏、无患子、杨树和樟树等7个主要树种组成。如图1中(b)所示，我们在校园内选择了三种实验立地类型进行实验：实验场地1为草木林纯针叶林；实验场地2是由鹅毛楸、榆、白蜡树等多种树种组成的阔叶林；实验场地3为阔叶-针叶林混交林，混合树种有香樟、桂花、日本樱桃、云杉、雪松等。在接下来的试验中，选择3个样地大小为50m×50m的森林类型(样地)各1个子集作为试验点；试验地点的位置和照片如图1中(b)和(c)所示。
图1中(a)为研究地点的位置，即中国江苏省南京市南京林业大学(谷歌地球)。图1中(b)为研究地点的放大图，三个不同的正方形标出了不同实验场地的边界。图1中(c)中从上至下依次为实验场地1、2、3生长的林分的对应照片。
研究区的激光雷达数据于2019年5月20日获得，使用Velodyne HDL-32E激光扫描仪，飞行高度为100米，飞行速度为10公里/小时，航线重叠率为30％。该扫描仪的波长为930nm，脉冲重复频率为21.7kHz，垂直视场(FOV)为-30.67°-10.67°，水平视场为360°。束发散度为2.79mrad。平均接地点距离10.5cm，脉冲密度约为100点/m2。图2中(c)、(e)、(g)分别显示了实验场地1、2、3(针叶、阔叶和混交林)采集的点云。
图2中(a)是具有可变参数的锥形模型，模拟每个针叶树的树冠。图2中(b)是具有可变参数的半椭球模型，模拟每个阔叶树的树冠。图2中(d)为实验场地1经过计算机建模的3D森林样地场景(针叶)，图2中(f)为实验场地2经过计算机建模的3D森林样地场景(阔叶)，图2中(h)为实验场地3经过计算机建模的3D森林样地场景(混交林)。
2.2、研究区太阳辐照度计算：
本实施例采用Iqbal模型来评估研究场地水平表面的太阳辐照度Itotal。包括来自太阳的直接辐照度Idirect，由瑞利散射、气溶胶散射组成的漫射地面辐照度Idiffuse，以及地面与天空之间的多次反射：
Itotal＝Idirect+Idiffuse (1)；
其中Idirect＝Iscr0τrτoτgτwτasinθa；Isc为太阳常数，取1367W/m2；r0为地球轨道偏心率相关系数。τr,τo,τg,τw,和τa分别是瑞利、臭氧、气体、水和气溶胶散射透过率。计算太阳高度角的表达式如下:
sinθa＝sinθlsinθd+cosθlcosθdcosθh (2)；
其中θl为当地纬度，θd为当前太阳的赤纬角，θh为当地太阳时的时角。
到达水平面的漫射辐照度可以表示为:
Idiffuse＝Idr+Ida+Idm (3)；
其中Idr和Ida分别表示大气中主要由空气分子和气溶胶散射的光，并描述了地面和天空之间的多次反射。可表示为:
τaa为气溶胶吸收引起的直接辐射的透射率，由给出。ωo为气溶胶散射到总衰减的入射能量的单次散射反照率，取0.9，ma是实际压力下的气团。
气溶胶引起的漫射辐照度计算如下：
其中Fc为前向散射相对于总散射的分数，取0.84，τas为气溶胶散射后透射入射能量的分数，计算公式如下：
地面与天空之间多次反射产生的漫反射辐照度可确定为：
ρg为地面反照率，ρa为无云天空的反照率，可由以下公式计算:
ρa＝0.0685+(1-Fc)(1-τas) (8)；
通过上述过程，可以对研究区的太阳辐照度进行评估。
2.3、使用规则几何基元进行树冠模拟：
采集真实林分的机载LiDAR数据后，对三种森林样地类型对应的扫描数据进行细化，去除噪声和离群值。然后，使用双高斯滤波器检测每个树冠，并基于能量函数最小化方法进行单个树冠分割。这样，将每个森林地块分解为一组单株树，从而便于检索每棵树的生长特性，包括南北和东西两个垂直方向的树冠宽度、树高、树中心坐标和枝下高。
三片林分由阔叶树和针叶树构成。利用半椭球和圆锥两种几何基元分别模拟阔叶和针叶树冠，并根据从点云中提取的每个独立树冠的生长特性分配自适应参数。
式(9)为阔叶树冠模拟的空间半椭球方程：
式(10)为针叶树冠模拟的空间圆锥方程：
其中a和b分别为每棵单树在东西、南北方向上的半冠宽度；c是树高减去枝下高；被分割树的每个树冠中心坐标记录为：(xcenter,ycenter,zcenter)，zcenter为相应树的枝下高；heighttree为相应树高。使用几何基元建模树冠的示意图如图2中(a)和(b)所示。
利用我们的方法，在选定的3个林分中，分别在样地1(即实验场地1)、样地2(即实验场地2)和样地3(即实验场地3)中获得了60、60和61株单株分段树。然后对每棵树进行空间显式的坐标和生长参数指定(如表2所示)，并以此来调整半椭球模型和圆锥模型来表示每棵树的树冠。图2中(d)、(f)和(h)显示了使用替代几何模型组合显示的三种林分仿真图像。
表2、树木生长特性：
2.4、树模型的三角剖分：
在地块中的每个树冠用表示后，进行三角剖分，将所有树冠表面分解为一组三角形Tj,j＝1,2,...M，并以三角形比例进行明确描述。这一过程确保了树冠的充分采样，三角形覆盖了森林冠层的外表面，并说明了太阳光束在树冠之间的传输和反射的详细相互作用。由于3DDelaunay三角剖分过于复杂且计算开销大，为了减少程序执行时间，本文采用2DDelaunay三角剖分。首先，将代表每个树冠的圆锥体和半椭球几何面从X-Y平面上的投影圆进行垂直投影，并在每个投影圆内定义多个等距的同心圆(图3)。在每个同心圆上均匀取相同点数，得到X-Y平面上的采样点。其次，利用二维Delaunay三角剖分算法，基于这些采样点生成一个三角形网格(图3中(a))，得到每个采样点的索引，定义构成树冠曲面(即树冠表面)的每个三角形的顶点。最后，根据采样点指标(即采样点结果)，对每个树冠表面进行反投影，实现三维Delaunay三角剖分。这种策略明显提高了程序的效率，使我们的树模型不再需要三维三角剖分。三角剖分树模型如图3所示。其中图3中(a)解释了采样点分布在每个树冠上并投影到X-Y平面上以优化三角剖分。图3中(b)为针叶林模拟几何基元的三角剖分结果。图3中(c)为阔叶林模拟几何基元的三角剖分结果。图3中(d)为混交林模拟几何基元的三角剖分结果。
2.5、离散化太阳辐射：
森林冠层上的太阳辐射是基于太阳光束与森林表面的相互作用，可以定义太阳光线在森林样地内的传播，定量评估太阳光束的入射、透射和反射。模拟太阳辐射的能量来自太阳在无穷远处，太阳光从林分上方同一水平面沿着同一方向和一定的间隔发射出入射光线，源点的坐标表示为其中N表示源点总数，每个源点发射一束太阳光束。这里，相邻两点之间的固定间隔设为0.2m；即发射的太阳光束密度为
3个样地的乔木高度在9～31m之间；因此，我们将所有源点的大小，即源点所在的森林地块上的水平面高度，设置为100m，以模拟阳光向下传播到森林冠层上。
每个太阳光束从每个源点出发的方向表示为单位矢量(这里的指的是入射光线的方向，三个参数分别入射光线在x轴、y轴和z轴上的分量。该向量可以通过方位角和高度角来得出)，与当前太阳方位角和高度角相关。因此，已知方向和源点的模拟太阳光束Lk可以表示为：
由于入射太阳光束的密度为加上2.2节中计算的研究场地的太阳辐照度Itotal(W/m2)，则每束发射的太阳光束的辐射通量可表示为
2.5.1、局部受辐射树冠面积的确定：
当太阳光照射到树冠时，必须确定树冠暴露在进行光合作用的入射光下的局部区域。在这里，利用生成的三个顶点和对应的三角形Tj和模拟的太阳光束Lk，可以采用光束相交算法，确定模拟树冠在三角形尺度下的局部光照和阴影区域。
每个三角形的法向量可以由其三个顶点来计算：
其中×表示两个向量的外积。结合入射光公式(11)，采用公式(13)确定射线Lk是否与三角形Tj所在平面相交。
如果F＝0，发射的太阳光线平行于三角形的平面，且不与该平面相交。如果F≠0，交点存在。交点坐标为：
计算交点后，采用以下准则判断交点是否在三角形Tj内：
其中，N1,N2和N3由式(16)导出：
通过以上计算，可以确定每个太阳光束与剖分森林冠层表面的每个三角形的交点。由于整个森林中植物元素的相互遮挡效应，每个光束都被光束遇到的第一个三角形所拦截。因此，可以确定不同太阳角度下直接太阳辐照度对森林冠层的三角尺度面积。这解决了复杂森林冠层中异质树冠之间的重叠问题。图4中(a)和(d)为抽取的入射太阳光束，其中曲线表示当前太阳光束入射角下的相互遮挡区域。具体为：图4中(a)为发射的入射太阳光束照射具有遮挡区域的邻近树木(曲线)。图4中(b)解释了大部分低仰角的入射太阳光束被反射到相邻的树冠上，而高仰角的太阳光束则被反射到天空或下层的灌木丛中。图4中(c)解释了不同的太阳方位角导致不同程度的树冠重叠效应，阻碍了太阳光束的传输路径。图4中(d)为入射光束被抽取以描绘太阳光在森林地块之间的传播。图4中(e)为反射光束被抽取以描绘太阳光在森林地块之间的传播。图4中(f)为透射光束被抽取以描绘太阳光在森林地块之间的传播。
2.5.2、入射辐射通量计算：
在确定了被照射的三角形之后，我们需要进一步量化覆盖森林冠层的每个被照射三角形的入射辐射通量，然后外推整个森林林分的入射辐射通量。
基于上述冠层光模拟，可以确定入射太阳光束与三角形的第一个交点。假设每个三角形Tj与发射的太阳光束的第一个交点总数为即入射的光束的总数首先被三角形拦截。那么，森林图的入射辐射通量或入射辐射功率可表示为：
其中为入射太阳光束Ld与被照射三角形Tj法向量的夹角，表示上述每个太阳光束的辐射通量。
2.5.3、反射辐射通量计算：
反射比入射更复杂。根据几何光学中的反射定律，光线入射到表面的角度的镜面反射等于它被反射的角度。入射太阳光被森林冠层不同三角形反射后的反射矢量可由下式求得：
在计算了反射太阳光束照射森林冠层三角形的反射矢量后，假设三角形Tj与所有反射的太阳光束的第一个交点总数为即假设一些具有反射能的反射太阳光束以给定的反射角继续传播，并与森林冠层的其他三角形Tj进行第二次相互作用。那么，森林图中树木间的反射辐射通量可以用以下公式表示：
公式(19)定义了相交的三角形内点与反射的太阳光束的能量,τreflect是三角面所属不同树种的树冠反射率(阔叶树和针叶树的相应值分别设置为0.17和0.25)；为反射的太阳光束与三角形法向量的夹角；q是一个时间变量，表示在与森林冠层三角形产生第二次碰撞的反射后的每一束入射太阳光。图4中(b)和(e)分别显示了入射高度相对较低和较高的太阳光束所产生的反射光在森林冠层之间产生的各种局部反射。
2.5.4、透射辐射通量计算：
当入射太阳光照射到林冠表面时，部分光能被植物元素反射和吸收，其余光能穿透林冠表面，进入林分中下层。将冠层的内部和外部作为两种不同的介质，进入冠层的透射光束可以看作是折射光束。这里，按照Snell定律模拟折射，即对于给定的一对介质，入射角θ1和折射角θ2的正弦之比等于两种介质之间的折射率。根据几何光学中的折射定律，矢量(这里的指的是透射光线的方向，三个参数分别透射光线向量在x轴、y轴和z轴的分量。)可由以下公式导出：
其中e为折射率，本研究取0.5。
目标林分的透射辐射通量计算公式如下：
公式(22)计算仅与透射太阳光束相交的三角形点，其中为所有透射太阳光束与三角形Tj的第一次交点总数，τtrans为各树种的冠层透光率(阔叶树和针叶树的对应值分别设为0.2和0.15)，是透射(折射)太阳光束与被光束照射的三角形法向量之间的夹角。如图4中(a)、(b)和(c)所示，一些入射、反射和透射的太阳光束沿给定方向在混合树种中传播，这些混合树种对太阳高度角和方位角始终很敏感。
2.6、方法验证：
为了进一步验证我们方法的有效性，使用现有的两种方法，即基于高分辨率HP计算冠层短波辐射传输的HPEval软件和智能辐射计传感器(S-LIB-M003)来测量单位面积的太阳能功率，利用该方法计算了森林冠层截获的太阳辐射，并与我们的方法进行了比较。
HPEval是一个基于matlab的软件工具，旨在评估数字HP并提供局部短波辐射估计。该工具适用于森林冠层以下的应用，可以输出高分辨率的时间序列，并有潜力解决单个太阳斑点。缺省情况下，HPEval输出潜在短波辐射与相应的冠层上方辐射，但也可以计算冠层拦截的真实太阳辐射。整个过程包括图像采集、元数据采集、图像二值化和天空像素分割。HPEval的HPs输入由高分辨率数码相机和鱼眼镜头拍摄。相机放置在森林内的地面上，光轴向上，拍摄树冠的HPs。为了让图像与太阳的天文路径参数相交，相机必须水平并朝向北方。由此产生的照片记录了植物冠层的结构，以及通过冠层上的孔或缝隙可见的天空部分。这些特征可以被HPEval软件精确测量并用来计算通过植物冠层的太阳辐射。
使用日射强度计对三个森林地块进行短波辐射测量。图5为基于HP的HPEval软件对三个森林地块进行短波辐射测量的HP图像。每个林分以12.5m的固定间距捕获10个HP。HP图像中曲线代表图像采集当天的太阳路径。
图5显示，每个林分采集10个HPs。相邻的HP的圆被12.5米的固定间隔隔开。我们用下面的公式来计算森林图截获的太阳(短波)辐射Sintercept：
Sintercept＝Sabove-Sbelow＝Sabove-(Sdif·Vf+Sdir·τdir) (24)；
其中，冠层上的入射太阳(短波)辐射Sabove等于第2.2节中计算的瞬时总太阳辐照度，冠层下的入射太阳短波辐射Sbelow由漫射分量Sdif和直接分量Sdir组成，可由当时的太阳高度角和太阳常数推导。Vf为从HPs计算得到的天景部分，τdir为直接短波辐射的透过率(本研究假设为0.2)。
每个地块(50米×50米)的面积为2500平方米。因此，利用HPEval软件从HPs计算出各林分的平均入射辐射通量可以表示为以下公式：
这里的i是由于我们每片林分都拍摄了10个HP，因此基于这10个HP，我们计算出了10个HP的平均Sintercept，再乘以面积2500m2即为HPEval方法计算出的辐射通量。
图5所示的硅智能日辐射计传感器(S-LIB-M003)的仪器测量方法对本研究也具有重要的指导意义。根据每个地块的采样HPs中心点位置，利用无人机将传感器垂直提升至冠层表面上方，测量冠层上方的太阳辐射，而冠层下的太阳辐射则是在HPs中心点近地面的另一组测量位置进行测量的。竭尽全力确保两组辐射传感器的测量都是水平的。该智能辐射计传感器有一个插件模块连接器，使其可以很容易地添加到一个站。所有校准参数都存储在智能传感器内，无需任何额外的编程、校正或广泛的用户设置，即可自动将其配置信息传输到记录器。光谱范围为300～1100nm，测量范围为0～1280W/m2。采用工厂重新校准程序后，280～2800nm光谱响应的单位面积太阳能功率可以外推为森林小区截获的太阳(短波)辐射。具体计算公式如下：
这里的Iabove指的是仪器测量出的林分上方的辐射度，Ibelow指的是仪器测量出的林分下方的辐射度。作差即为截获辐射度。其单位为W/m2。
研究结果：
3.1、研究区的太阳辐照度计算：
通过太阳辐照度的计算方法，在当地，结合2.2节中阐述的经度，纬度，累积的日子和高度角，以及相应的空气质量参数，校园的太阳能辐照度可以实现在每时每刻被计算。
在一天的不同时间，林分在太阳照射下的暴露情况不同。太阳高度角和方位角随时间变化，从而影响研究区的太阳辐照度。例如在北半球，夏至中午的太阳高度角明显高于冬至中午的太阳高度角。根据2.2节的Iqbal模型，选择2019年3月至10月的每个月15日，在晴空条件下，以1小时为间隔计算8:00～17:00时水平地表太阳辐照度；结果如图6所示。在早上和下午，太阳在天空中位置较低，因此，它的光线穿过大气层的距离更远，从而降低了太阳辐照度。同样，在夏季(6-8月)，太阳高度角大于同期其他月份时，太阳辐照度较大。
图6中(a)计算了研究地点在不同日期和一天中不同时间的总太阳辐照度。图6中(b)计算了研究地点在不同日期和一天中不同时间的直接太阳辐照度。图6中(c)计算了研究地点在不同日期和一天中不同时间的漫射太阳辐照度。
3.2、利用反演的辐射通量对三个林分进行光照渲染：
太阳光束从离地面100m的水平面(z＝100m)发射，光束密度为入射角度与当前太阳方位角和高度角有关。每个光束被分配一个初始能量指向由许多三角形面组成的离散森林表面区域。这些光束与这些三角形的相互作用，例如它们的吸收、反射和透射，如图7所示。
各向异性林型和异质茎杆结构影响了林分入射太阳光束的空间分布；例如，树冠的几何结构变化导致树冠之间相互遮挡，产生局部阴影，进而影响光线跟踪算法反射和透射光束在林分之间的传播轨迹。图7显示了在2019年7月15日11点，入射，反射和透射光束截获的三角形组成的森林的树冠三个示例图，当时的太阳光线入射方位角为121.88°，高度角为73.53°。
换句话说，图7展示了在三角形尺度上不同照明层次的分布，其中太阳光束局部照亮树木，并根据光线跟踪算法通过森林冠层的几何表示在森林林分中传播。
图7中(a)为2019年7月15日11:00针叶林样地冠层入射辐射通量分布。图7中(b)为2019年7月15日11:00阔叶林样地冠层入射辐射通量分布。图7中(c)为2019年7月15日11:00混交林样地冠层入射辐射通量分布。图7中(d)为2019年7月15日11:00针叶林样地冠层反射辐射通量分布。图7中(e)为2019年7月15日11:00阔叶林样地冠层反射辐射通量分布。图7中(f)为2019年7月15日11:00混交林样地冠层反射辐射通量分布。图7中(g)为2019年7月15日11:00针叶林样地冠层透射辐射通量分布。图7中(h)为2019年7月15日11:00阔叶林样地冠层透射辐射通量分布。图7中(i)为2019年7月15日11:00混交林样地冠层透射辐射通量分布。
对于斜光照的入射太阳光，入射辐射通量集中在林冠面对太阳一侧，而对侧辐射通量较低且有较多的阴影区。另外，大树的相互遮挡也会影响周围小树的日照照度。在稀疏辐亮度分布下，反射辐射通量主要出现在林分中下层。太阳光通过林冠隙或植物元素的透射和折射进入林分。采用模拟光束与林冠背面相互作用的分布，计算了林冠内透射光的辐射通量。
无论同质林分的组成如何，图8说明了入射辐射通量在大约中午达到峰值。这主要是由于南京林业大学的太阳高度角在12:00左右达到最大值，这导致了太阳辐照的区域幅照度最高。此外，由于植物趋光性和树木竞争，所有树木垂直生长，树冠在空间上伸展；因此，太阳发射太阳光的高度角与树木的生长方向是一致的，从而消除了森林林分中树冠之间的重叠区域，而不考虑亚冠(下层)植被对光线的吸收。然而，随着太阳高度角的减少，例如，在清晨和傍晚，太阳高度角在(近15°)时候，由于树之间的阻碍，导致入射辐射通量在17:00达到最低值(129.26-756.51kW)。结合表2中列出的树木生长属性，如图8中(a)，(d)和(g)所示，入射辐射通量在针叶林中出现最大值，由于最高的针叶树有相对较大的树冠，从而形成一个紧密的森林冠层和更大的截留面积，用于斜向传播阳光。相比之下，纯阔叶林由于成熟树木最矮，树冠相对较小，因此接收入射辐射通量最小，从而在森林景观中产生更多的开放空间。
图8中(a)为使用本实施例的方法计算出来的针叶林样地的入射辐射通量；图8中(b)为使用本实施例的方法计算出来的阔叶林样地的入射辐射通量；图8中(c)为使用本实施例的方法计算出来的混交林样地的入射辐射通量；图8中(d)为使用本实施例的方法计算出来的针叶林样地的反射辐射通量；图8中(e)为使用本实施例的方法计算出来的阔叶林样地的反射辐射通量；图8中(f)为使用本实施例的方法计算出来的混交林样地的反射辐射通量；图8中(g)为使用本实施例的方法计算出来的针叶林样地的透射辐射通量；图8中(h)为使用本实施例的方法计算出来的阔叶林样地的透射辐射通量；图8中(i)为使用本实施例的方法计算出来的混交林样地的透射辐射通量。
反射和透射辐射通量与入射光束与树冠几何基元碰撞后反射和透射光束的再碰撞次数相关。许多具有新方向的反射光可能传播到林下层的灌木或森林地块外，而不会向附近的林分传递辐射。相反，所有入射光都穿透树冠表面，剩余的能量被树冠另一侧的三角形吸收。如2.5节所述，指定的针叶树冠层反射率(0.25)大于阔叶树的(0.17)，而针叶树的透过率(0.15)小于阔叶树的透过率(0.20)。结合反射或透射光束与目标表面的夹角，在针叶树图中反射辐射通量大于透射辐射通量，在阔叶树图中则相反。相比之下，混交林样地的针叶树(33株)和阔叶树(28株)数量大致相同；因此，估算的入射通量与反射辐射通量近似。
当前森林样地的太阳光照面积与树木密度、生长树木高度、树冠形状、太阳高度和方位角有关(图9中(a))。相对应的角度入射太阳能光方向的三个观察森林呈现在图9中(b)、(c)和(d)，这表明，针叶林(b)有一个较大的区域接收入射太阳光束，其次是混交林(d)和阔叶林(c)在降序排列。
图9中(b)为正午太阳入射光束视角下的针叶林分图；图9中(c)为正午太阳入射光束视角下的阔叶林分图；图9中(d)为正午太阳入射光束视角下的混交林分图；三片森林的视角与2019年8月15日中午的太阳入射光束方向平行，太阳高度角为73.01°，太阳方位角为180.00°。树冠由空间三角形斑块组成。
3.3、现有方法的比较：
表3为我们的方法、基于HP的HPEval软件和日强度辐射仪器三种方法计算得到的林分的辐射通量。三种方法的相关系数均大于0.83。对于我们的方法与基于HP的软件估计和基于仪器的现场测量之间存在偏差的不完美对齐的合理解释是，一个地块中每个采样点的距离和方向需要在现场严格标准化。提高采样点密度会增加两种方法之间的相关系数，但由于人工操作和判断误差造成的数据扰动会产生偏差。对于HPEval软件，错误地设置HPs二值化阈值可能会导致孔隙度被高估或低估，不同强度光照条件下叶片反射率和透射率的均匀分布，并不能很好地反映出图像灰度值变化的实际生理状况。这些因素导致基于HPEval计算的太阳辐射通量与实际值存在一定偏差。此外，为了有效地防止环境因素干扰内部元件，仪器测量的值应该更接近真实值，但传感器在测量过程中需要大量的人力和准备工作，例如定位采样点和进行认可的校准；在强风和不同于阳光的光谱含量条件下，仪器测量也很困难。因此，在仪器数量有限的情况下，在同一时间在一个图的所有采样点收集现场数据是不现实的。因此，为了获得实地测量，需要迅速地将仪器携带到每个采样地点，以追踪太阳辐照度随时间的变化。此外，传感器上的灰尘降低了其精度，并限制了电池寿命中断和限制了设备的连续运行。
表3、利用该方法反演的辐射通量与基于HP的HPEval软件和日射计传感器获得的辐射通量的比较：
备注：O-H：我们的方法与HPEval之间的相关系数；O-I：我们的方法与仪器测量之间的相关系数；O-I：HPEval与仪器测量之间的相关系数。
讨论：
4.1、森林太阳辐射计算的展望：
从机载激光雷达数据中获取树木生长属性，实现单株树冠分割算法后，利用基于计算机图形学的射线追踪模型和模拟树冠结构的基本几何模型，对森林树冠接收的辐射通量进行估算。树冠内不规则分布的微小叶片会产生多用途的镜面和漫反射光，这极大地复杂化了森林冠层中太阳辐射通量的计算。因此，对每个树冠应用一种简化的近似方法，该方法基于入射太阳光照射下的规则几何原元，并伴随着反射和透射率，以获得每个森林样地的辐射通量。日期和时间的变化导致了太阳高度和方位角度的变化，产生不同的照明角度，因此植被的细节元素的变化与太阳的光束从天空中太阳的当前位置。由于北半球中午的方位角始终为180°(南)，因此在不同月份(3月15日至10月15日)，中午的太阳高度角是一个重要因素。如图6所示，中午时分，南京市太阳高度角在6月15日达到峰值(81.21°)，研究区太阳总辐照度为920.26w/m2。因此，6月15日中午针叶林样地的入射辐射通量为2047.2kW；阔叶林1653.63kW；混合林1785.76千瓦)高于中午的辐射通量在另一个月(从10月的1834.66千瓦到7月份的2018.65千瓦的针叶林，从3月的1440.54kW至7月的1589.81kW的阔叶林，从10月的1530.54kW至7月的1746.1kW的混合林)。随着太阳高度角的减小，森林样地的遮挡效应会增强，因为这些效应是由树冠的相互遮挡的侧面表现出来的，这导致了入射辐射通量的下降。此外，方位角的变化伴随着太阳高度角的变化，为目标林分的入射辐射通量和不同光照视野的变化提供了比通常预期更多的可能性。
我们发现，冠层填充和样地空间填充增加了每个林分相应的辐射通量值。针叶林、阔叶林和混交林的冠层总容积分别为25151、13546和15391m3。虽然有些研究复杂的物种组成是潜在影响森林动态速度、生态系统生产力恢复力、冠层结构可塑性的重要因素。这将提高森林冠层密集的光合组织对太阳辐射的截留能力，而这些组织对太阳光线更难以穿透，但实际上，树木的形态特征，例如树冠大小和树高，都受到森林再造林活动的限制。耕作干预强度与生态微环境条件。第一个地块的针叶树年龄较大，树体更大、更深，树冠呈锥形，这增加了森林冠层的垂直分层，有利于太阳光覆盖整个冠层。相反，样地2、3的阔叶混交树树冠较小，占空间能力和树冠间相互作用减弱，导致大量太阳光直接照射地面，而没有充分利用光。
虽然我们的方法在提取树冠间大量小叶片辐射通量和分布之间的对应关系方面表现相对较弱，但我们的方法不需要大量的射线来充分取样小的植物元素，并结合激光雷达数据导出的树种冠几何形状变化的树木生长特性，实现了小区级和种级辐射通量之间的平衡。此外，考虑到冠层反射率、冠层透光率、冠层形状和冠层大小都是可调节的，并利用灵活数量的发射太阳光束和覆盖森林冠层的细分三角形，根据太阳方位角和高度角，采用定向射线跟踪算法估计碰撞和再碰撞概率，模拟太阳光束在林冠间的入射、反射和透射。因此，我们的系统具有更强的实用性和可扩展性，适用于任何林分类型和复杂、结构多样的冠层，因此可以渲染植被样地场景所接收到的辐射。
4.1、方法参数和算法执行：
我们的方法中主要参数的具体设置如下。在树冠建模过程中，我们利用机载激光雷达数据导出的生长特性，根据标准形状的针叶树为圆锥，阔叶树为半椭球，对森林冠层组成进行建模。因此，利用三角剖分策略将森林冠层在三角尺度上细分为不同的详细区域。根据研究点的太阳辐照度计算出集合能量的每束离散太阳光束沿当前太阳高度角和天顶角方向发射，然后与代表冠面的三角形相互作用。原则上，三角形越小，入射光束密度越大，得到的森林辐射通量结果越准确，但程序执行时间越长。
我们改变了邻近的入射太阳光束之间的间距(从0.1米到0.4米)和平均间距相邻顶点的三角形(从1.01米到1.75米通过调优的同心圆投影圆的数量和每个同心的圆上的采样点)，这使得树木最高的针叶林样地的入射太阳光束数为25236～402741，高于阔叶样地(入射太阳光束数为16310～265161)和混交林样地(入射太阳光束数为19900～318003)。针叶林样地(三角形数为14156～56755)的三角形数也大于阔叶样地(三角形数为6705～20495)和混交林样地(三角形数为9178～29159)，因为针叶林样地树木较高，树冠较宽。双y轴图10显示了使用我们的方法对三个森林图估算辐射通量和时间消耗的RMSE随入射光束和三角形顶点间距的变化。三次多项式回归拟合得到的实线和虚线(即代表均方根误差的实线和虚线)表示我们的方法与参考场测量(基于HP的HPEval和日强度辐射仪器的平均值)在三个森林样地的辐射通量估计RMSE呈稳定上升趋势，这说明光束密集的光束分布和精细的三角形表示可以使均方根误差量级降低。相比之下，图10中代表运行时间的实线和虚线表明，我们的辐射通量计算方法的时间消耗先减少后趋于平稳。因此，权衡后设置入射光线间距为0.2米，三角形顶点平均间距1.36米。图10的灰色区域即为综合考虑程序执行时间和较小的RMSE的标记。图10为双y轴图显示了程序运行时间和使用我们的针叶林、阔叶林和混交林地块的方法计算的辐射通量RMSE的分布。线由三次多项式函数拟合。灰色区域代表我们方法中分配的参数设置，综合考虑了计算耗时和RMSE。
在树冠层和叶片水平的光学特性，例如阔叶和针叶树种的反射率和透射率，在之前的许多工作中已经进行了研究。一项开创性的研究表明，在可见(VIS)光谱(波长从380nm到750nm)中，针叶树的反射率(范围从0.04到0.38)高于阔叶树(范围从0.02到0.29)。在近红外光谱(波长为750～1400nm)中，针叶树的反射率在0.14～0.55之间，高于阔叶树(0.15～0.47)。只有在短波红外光谱(波长为1400～3000nm)中针叶树的反射率(0.02～0.28)略低于阔叶树(0.02～0.33)。在光谱透过率方面，针叶树(0.02-0.41)在VIS光谱中高于阔叶树(0.02-0.30)，而针叶树在近红外光谱(针叶树:0.05-0.32，阔叶树:0.18-0.55)和SWIR光谱(针叶树:0.01--0.21,阔叶:0.03--0.51)。因此，根据相关文献，针叶树和阔叶树的反射率分别为0.25和0.17，对应的透射率分别为0.15和0.2。结果表明，针叶林、阔叶林和混交林的最大反射辐射通量分别为324.65kW、164.37kW和217.02kW。13:00针叶树样地最大辐射通量为211.28kW，中午阔叶样地最大辐射通量为280.93kW,10:00混合样地最大辐射通量为209.44kW。在透射和入射辐射通量之间观察到高度相关，因为每一入射太阳光束击中树冠都会产生相应的透射光束，在树冠内传播，产生透射辐射通量。相反，反射和入射辐射通量之间的相关性相对较弱，因为反射光束的方向是不可预测的，这是因为随着时间的推移，不同方向的入射太阳光束照射到局部树冠区域的法向量各向异性；因此，击中亚冠层和在图外传播的反射光束的比例是可变的。
5、总结：
本研究基于三维冠层结构的基本几何表示，目的是开发一种定量分析森林冠层接收到的太阳光的入射、反射和透射的方法。森林冠层由从冠层收集的机载激光雷达数据建模，冠层表面被细分成三角形，允许提出的方法量化入射、反射、并通过太阳辐射模型和基于计算机图形学的射线追踪算法，传输不同森林样地冠层接收到的太阳辐射通量。提出了一种利用几何基元模型模拟森林冠层的光束相交优化算法，以模拟太阳光在林分中的传播。我们的结果表明，一个高度聚集分布的树木具有相对较大的树冠的针叶树小区将获得最大的入射和反射辐照通量。相比之下，阔叶地块由于其较大的透光率，具有较大的透射辐射通量。相关系数(>0.82)表明，我们的结果与现有方法，即HPEval软件和日射计仪器的测量结果具有较高的一致性。通过对参数进行适当的调整，所提出的方法可以应用于任何阔叶或针叶树树种以及任何森林样地类型。随着计算机图形算法和GPU性能的提升，我们的方法既可以根据需要调整光密度，无限接近实际太阳辐照度，又可以根据点云精细模拟树冠中的植物元素，并定量评估它们与每束太阳光线的生物物理相互作用。
本发明的保护范围包括但不限于以上实施方式，本发明的保护范围以权利要求书为准，任何对本技术做出的本领域的技术人员容易想到的替换、变形、改进均落入本发明的保护范围。
Claims (9)
1.基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：包括：
步骤1：通过激光雷达采集林分点云数据；
步骤2：计算林分水平表面的太阳辐照度；
步骤3：对采集的林分点云数据进行噪声和离群值的去除，然后，使用双高斯滤波器从林分点云数据中检测每个树冠，并基于能量函数最小化方法进行单个树冠分割；
步骤4：将林分内的树木种类分为阔叶树或/和针叶树，利用空间半椭球几何基元模拟阔叶树冠或/和利用空间圆锥几何基元模拟针叶树冠；
步骤5：对几何基元模拟后的树冠进行三角剖分；
步骤6：每个源点发射一束太阳光束，计算已知方向和源点的模拟太阳光束，并计算每束发射的太阳光束的辐射通量；
步骤7：利用三角剖分后生成的三个顶点对应的三角形和模拟太阳光束，采用光束相交算法，确定几何基元模拟后的树冠在三角形尺度下的局部光照和阴影区域；
步骤8：确定入射太阳光束与三角形的第一个交点，进而计算林分的入射辐射通量；
步骤9：计算入射太阳光被林分冠层不同三角形反射后的反射矢量，通过反射矢量计算林分中树木间的反射辐射通量；
步骤10：计算入射太阳光被林分冠层透射后的透射矢量，通过透射矢量计算林分的透射辐射通量。
2.根据权利要求1所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤2具体为：
林分水平表面的太阳辐照度Itotal为：
Itotal＝Idirect+Idiffuse (1)；
其中Idirect＝Iscr0τrτoτgτwτasinθa；Isc为太阳常数，取1367W/m2；r0为地球轨道偏心率相关系数；τr、τo、τg、τw和τa分别是瑞利、臭氧、气体、水和气溶胶散射透过率；
太阳高度角的表达式如下：
sinθa＝sinθlsinθd+cosθlcosθdcosθh (2)；
其中θl为当地纬度，θd为当前太阳的赤纬角，θh为当地太阳时的时角；
到达水平面的漫射辐照度Idiffuse为：
Idiffuse＝Idr+Ida+Idm (3)；
其中：
τaa为气溶胶吸收引起的直接辐射的透射率，ωo为气溶胶散射到总衰减的入射能量的单次散射反照率，取0.9，ma是实际压力下的气团；
Fc为前向散射相对于总散射的分数，取0.84，τas为气溶胶散射后透射入射能量的分数，计算公式如下：
ρg为地面反照率，ρa为无云天空的反照率，计算公式如下：
ρa＝0.0685+(1-Fc)(1-τas) (8)。
3.根据权利要求2所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤4中利用空间半椭球几何基元模拟阔叶树冠的空间半椭球方程为：
利用空间圆锥几何基元模拟针叶树冠的空间圆锥方程为：
其中a为每棵单树在东西方向上的半冠宽度；b为每棵单树在南北方向上的半冠宽度；c是树高减去枝下高；被分割树的每个树冠中心坐标记录为：(xcenter,ycenter,zcenter)，zcenter为相应树的枝下高；heighttree为相应树高。
4.根据权利要求3所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤5具体为：
将代表每个树冠的圆锥体或/和半椭球几何面从X-Y平面上的投影圆进行垂直投影，并在每个投影圆内定义多个等距的同心圆；在每个同心圆上均匀取相同点数，得到X-Y平面上的采样点；利用二维Delaunay三角剖分算法，基于采样点生成一个三角形网格，得到每个采样点的索引，定义构成树冠表面的每个三角形的顶点；根据采样点结果，对每个树冠表面进行反投影，实现三维Delaunay三角剖分。
5.根据权利要求4所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤6具体为：
定义源点的坐标为其中N表示源点总数，每个源点发射一束太阳光束；相邻两源点之间的固定间隔设为0.2m；即发射的太阳光束密度为l＝25/m2；
每个太阳光束从每个源点出发的方向表示为单位矢量
已知方向和源点的模拟太阳光束Lk为：
由于发射的太阳光束密度为l＝25/m2，则每束发射的太阳光束的辐射通量为Itotal/l。
6.根据权利要求5所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤7具体为：
每个三角形Tj,j＝1,2,...M的法向量可以由其三个顶点来计算：
其中×表示两个向量的外积；结合公式(11)，采用公式(13)确定射线Lk是否与三角形Tj所在平面相交：
如果F＝0，发射的太阳光线平行于三角形的平面，且不与该平面相交；
如果F≠0，交点存在，交点坐标为：
计算交点后，采用以下准则判断交点是否在三角形Tj内：
其中，N1,N2和N3由式(16)导出：
通过以上计算，进而确定每个太阳光束与剖分森林冠层表面的每个三角形的交点。
7.根据权利要求6所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤8具体为：
假设三角形Tj与所有入射的太阳光束的第一个交点总数为即入射的太阳光束的总数首先被三角形拦截，那么，林分的入射辐射通量或入射辐射功率表示为：
其中为入射太阳光束Ld与被照射三角形Tj法向量的夹角。
8.根据权利要求7所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤9具体为：
入射太阳光被森林冠层不同三角形反射后的反射矢量为：
假设三角形Tj与所有反射的太阳光束的第一个交点总数为反射太阳光束以给定的反射角继续传播，并与森林冠层的其他三角形Tj进行第二次相互作用，那么，林分中树木间的反射辐射通量表示为：
τreflect是三角面所属不同树种的树冠反射率；为反射的太阳光束与三角形法向量的夹角。
9.根据权利要求8所述的基于点云数据和计算机图形学的林分辐射通量计算方法，其特征在于：所述的步骤10具体为：
入射太阳光被林分冠层透射后的透射矢量为：
其中e为折射率，取0.5；
林分的透射辐射通量计算公式如下：
其中为所有透射太阳光束与三角形Tj的第一次交点总数，τtrans为各树种的冠层透光率，是透射太阳光束与被光束照射的三角形法向量之间的夹角。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210276654.1A CN114663786B (zh) | 2022-03-21 | 2022-03-21 | 基于点云数据和计算机图形学的林分辐射通量计算方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202210276654.1A CN114663786B (zh) | 2022-03-21 | 2022-03-21 | 基于点云数据和计算机图形学的林分辐射通量计算方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN114663786A CN114663786A (zh) | 2022-06-24 |
| CN114663786B true CN114663786B (zh) | 2024-03-22 |
Family
ID=82031910
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202210276654.1A Active CN114663786B (zh) | 2022-03-21 | 2022-03-21 | 基于点云数据和计算机图形学的林分辐射通量计算方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN114663786B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115292616B (zh) * | 2022-06-30 | 2023-07-18 | 北京大学 | 基于光谱不变理论的植被蓝天空反照率估算方法及装置 |
| CN115963092A (zh) * | 2022-12-07 | 2023-04-14 | 浙江大学 | 基于浊度补偿和散射宽度估计的自适应瑞利散射处理方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102914501A (zh) * | 2012-07-26 | 2013-02-06 | 南京大学 | 一种利用激光点云计算三维森林冠层消光系数的方法 |
| CN105371789A (zh) * | 2015-10-09 | 2016-03-02 | 南京大学 | 一种利用航空激光点云计算有效叶面积指数的方法 |
| CN105389538A (zh) * | 2015-10-09 | 2016-03-09 | 南京大学 | 一种基于点云半球切片估算森林叶面积指数的方法 |
| CN105701313A (zh) * | 2016-02-24 | 2016-06-22 | 福州大学 | 多层数据结构的虚植物冠层光合有效辐射分布模拟方法 |
| KR102008017B1 (ko) * | 2018-02-12 | 2019-08-06 | 전주비전대학교산학협력단 | 드론을 이용한 태양광발전소 건설부지의 일사량 계산 방법 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AUPR301601A0 (en) * | 2001-02-09 | 2001-03-08 | Commonwealth Scientific And Industrial Research Organisation | Lidar system and method |
-
2022
- 2022-03-21 CN CN202210276654.1A patent/CN114663786B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102914501A (zh) * | 2012-07-26 | 2013-02-06 | 南京大学 | 一种利用激光点云计算三维森林冠层消光系数的方法 |
| CN105371789A (zh) * | 2015-10-09 | 2016-03-02 | 南京大学 | 一种利用航空激光点云计算有效叶面积指数的方法 |
| CN105389538A (zh) * | 2015-10-09 | 2016-03-09 | 南京大学 | 一种基于点云半球切片估算森林叶面积指数的方法 |
| CN105701313A (zh) * | 2016-02-24 | 2016-06-22 | 福州大学 | 多层数据结构的虚植物冠层光合有效辐射分布模拟方法 |
| KR102008017B1 (ko) * | 2018-02-12 | 2019-08-06 | 전주비전대학교산학협력단 | 드론을 이용한 태양광발전소 건설부지의 일사량 계산 방법 |
Non-Patent Citations (4)
| Title |
|---|
| 利用辐射度模拟虚拟枇杷果园冠层光分布;杨怡斐;唐丽玉;崔磊;陈舒炜;;计算机工程与应用;20170627(第11期);全文 * |
| 基于虚拟现实的杉木形态建模及冠层辐射模拟;陈刚;李鑫;;福建电脑;20130625(第06期);全文 * |
| 激光雷达和点云切片算法结合的森林有效叶面积指数估算;路璐;郑光;马利霞;;遥感学报;20180525(第03期);全文 * |
| 鹤山人工林的辐射能环境研究;任海, 彭少麟;生态科学;19970630(第01期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN114663786A (zh) | 2022-06-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zheng et al. | Computational-geometry-based retrieval of effective leaf area index using terrestrial laser scanning | |
| Wei et al. | An assessment study of three indirect methods for estimating leaf area density and leaf area index of individual trees | |
| Widlowski et al. | The fourth phase of the radiative transfer model intercomparison (RAMI) exercise: Actual canopy scenarios and conformity testing | |
| CN114663786B (zh) | 基于点云数据和计算机图形学的林分辐射通量计算方法 | |
| Morisette et al. | Validation of global moderate-resolution LAI products: A framework proposed within the CEOS land product validation subgroup | |
| CN102314546B (zh) | 基于虚拟植物的植物生长生物量变化估算方法 | |
| US8019117B2 (en) | Method and apparatus for evaluating solar radiation amount | |
| Xue et al. | Shortwave radiation calculation for forest plots using airborne LiDAR data and computer graphics | |
| Van Leeuwen et al. | Automated reconstruction of tree and canopy structure for modeling the internal canopy radiation regime | |
| Zhao et al. | A spectral directional reflectance model of row crops | |
| Ma et al. | Retrieving forest canopy clumping index using terrestrial laser scanning data | |
| Ma et al. | Characterizing the three-dimensional spatiotemporal variation of forest photosynthetically active radiation using terrestrial laser scanning data | |
| Wang et al. | Evaluating a three dimensional model of diffuse photosynthetically active radiation in maize canopies | |
| Schraik et al. | Crown level clumping in Norway spruce from terrestrial laser scanning measurements | |
| CN114066966A (zh) | 基于多孔介质理论和计算机图形学的树冠孔隙率估算方法 | |
| Lin et al. | Reconstruction of a large-scale realistic three-dimensional (3-D) mountain forest scene for radiative transfer simulations | |
| Chen et al. | Effect of layer thickness and voxel size inversion on leaf area density based on the voxel-based canopy profiling method | |
| Hovi et al. | Empirical validation of photon recollision probability in single crowns of tree seedlings | |
| CN116881609A (zh) | 一种普适性的森林冠层间隙率的计算方法 | |
| Guo et al. | Apple tree canopy leaf spatial location automated extraction based on point cloud data | |
| Mõttus | Measurement and modelling of the vertical distribution of sunflecks, penumbra and umbra in willow coppice | |
| Kesselring et al. | Diversity of 3D APAR and LAI dynamics in broadleaf and coniferous forests: Implications for the interpretation of remote sensing-based products | |
| Wang et al. | Estimating photosynthetically active radiation distribution in maize canopies by a three-dimensional incident radiation model | |
| Kwak et al. | Estimation of effective plant area index for South Korean forests using LiDAR system | |
| Sakai et al. | FLiES-SIF ver. 1.0: Three-dimensional radiative transfer model for estimating solar induced fluorescence |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |