激光雷达的立体体元模型在森林冠层间隙率提取中的应用

曹庆安 冷 鹏，2 何原荣

（1.江西省核工业地质调查院，江西 南昌 330038；2.福州大学 空间数据挖掘与信息共享教育部重点实验室，福建 福州 350116；3.厦门理工学院 数字福建自然灾害监测大数据研究所，福建 厦门 361024）

0 引言

森林是陆地生态系统重要组成部分，具有调节气候、保持水土、防风固沙等多种功能。叶面积指数（leaf area index，LAI）是描述森林冠层结构的重要参数之一，控制着森林冠层的光合作用、呼吸作用等多种生物物理过程［1-2］。而LAI 的计算一般通过冠层间隙率进行反演，为此，准确有效地估算森林冠层间隙率对准确反演LAI具有重要的意义。

森林冠层间隙率的定义是指一束光沿一定的方向照射在森林冠层，未被森林冠层阻挡的概率。目前常用的测量方法有：半球成像（digital hemispherical photography，DHP）、植物冠层叶面积指数测量仪、跟踪辐射和冠层结构测量仪（tracing radiation and architecture of canopies，TRAC）、DCP（digital cover photography）、多光谱冠层成像仪（multi-spectral canopy imager，MCI）等，这些方法受天气的影响均较大［3-7］。而激光雷达（light detection and ranging，LiDAR）作为一种新型遥感技术，有全天候作业、受天气影响小、精度高的优势，越来越多的学者将激光雷达点云数据运用于森林冠层间隙率提取［8-10］。如CHEN Y 等人基于激光雷达回波强度信息来反演间隙率，但该方法受回波强度校正精度的影响较大；ZHENG G 等人提出采用点云投影方法将3D 点云投影成2D 图像提取间隙率，但该方法受采集数据缺失程度影响较大［11］。因此，本文基于激光雷达点云数据，提出一种立体体元模型法进行间隙率提取，该方法主要分为两部分，一是通过设置体元大小，构建以体元为介质的森林场景，然后通过半球成像原理，生成光线集合，采用光线跟踪算法，对光线与体元进行求交运算，最后根据生成的光线输出DHP影像；二是根据生成的DHP影像，统计各个天顶角区间的黑白像素与空白像素，从而求得各个天顶角区间的间隙率。

1 理论与方法

1.1 半球成像方法

半球成像方法是采用极化投影方式生成，极化投影采用式（1）计算，其基本原理如图1 所示，O为光源点，点C为成像平面的任意点有一条光线与其对应，点B为点C在半球S上对应的投影点，则点C对应的出射光线为OB，ϕ为光线OB的方位角，θ为光线OB的天顶角，D为C到O的距离，R为成像平面半径［12］。

图1 半球成像原理

1.2 点云投影方法

点云投影方法是通过模拟半球成像方法的原理生成DHP 影像，然后基于生成的DHP 影像提取森林冠层间隙率。基本原理如下：先对采集的单站点云数据设置投影中心点与生产的DHP影像大小，从而构建像平面；然后将所有点云利用极化投影方法计算各个点云投影后对应像素，并对该像素赋值；最后根据像素值生成对应的影像［13］。

1.3 立体体元模型方法

立体体元模型方法的基本原理是基于点云数据，构建体元为介质的虚拟森林场景，然后基于虚拟森林场景采用光线跟踪算法进行DHP 影像模拟，最后基于DHP 影像提取间隙率。具体原理：首先将计算样地所有点云数据的最小包围盒；然后根据设定的体元大小将最小包围盒体元化，并计算体元内是否包含点云数据，将所有不包含点云数据的体元剔除，保留包含1 个点云数据以上的体元构建森林场景；然后基于构建的森林场景，并根据半球成像原理，将成像平面每个像素生成光线集合，采用光线跟踪算法进行光线与体元相交计算，最后根据光线求交结果生成对应的DHP影像。

1.4 间隙率计算方法

本文间隙率计算方法是基于DHP 影像，采用图像处理技术，统计各个天顶角区间总像素与黑白像素数量，然后根据式（2）计算各个天顶角方向间隙率［14］。

式中，θ为天顶角；Pgap(θ)为θ天顶角区域对应的间隙率；Nblank(θ)为θ天顶角区域空白像素个数；Nsum(θ)为θ天顶角区域像素总数。

2 数据采集与处理

2.1 数据收集与预处理

本文实验样地在河北承德塞罕坝森林公园，样地树种为落叶松，样地内单树随机分布。实验样地大小为25 m×25 m，点云数据采样点均匀布设在样地中，共布设16 个采样点，每个采样点间隔5 m。数据采集使用Riegl VZ-400地面扫描仪，扫描角分辨率为0.04°，测站高为1.2 m，每个采样点均采集0°和90°姿态，点云数据采集时需要无风。DHP 数据使用Canon 6D 相机搭配Sigma 8 mm鱼眼镜头进行采集，采样点与点云数据采样点一致。

点云数据的处理主要包括配准与去噪，各个测站点云数据使用RiSCAN Pro 软件进行配准，以其中一个测站为基准站，先使用反射片模式配准，对不能使用反射片配准的测站，通过手动选择同名点进行配准，配准中最大误差为3.8 mm。配准后，需要进行去噪处理，先采用RiSCAN Pro软件中的Deviation 去噪功能，设置去噪阈值进行噪点删除，对于少量剩余零散的噪点，采用手动删除［15］。

外业实地拍摄的DHP 影像数据需要进行裁剪、阈值化处理［16］，先通过Photoshop 将图像进行有效区域进行裁剪，然后英语伽马校正与阈值化将影像进行分类，处理前后图像如图2所示。

图2 DHP影像处理前后对比图

2.2 基于点云投影的DHP模拟

根据1.2 节中点云投影方法的基本原理，该方法图像模拟具体流程如图3所示。首先加载去噪后的单站点云数据，计算每个点云的方位角与天顶角；其次设置模拟影像的分辨率与投影中心坐标，生成像平面；然后计算每个点云到成像中心点距离以及投影后对应的位置与像素；最后遍历完成所有点云数据，生成DHP 影像，生成的DHP影像如图4所示。

图3 点云投影方法模拟DHP影像流程图

图4 点云投影DHP影像

2.3 基于立体体元模型的DHP模拟

立体体元模型方法与点云投影方法存在较大差异，该方法采用16 个测站配准后的点云数据，由于数据量大，计算量巨大，本文调用点云库（point cloud library，PCL），构建八叉树数据结构进行加速，同时模拟DHP 影像的采样方案与实地采集点云的采样方案一致，共16 个采样点，模拟高度为1.2 m，DHP 影像分辨率3 000 像素×3 000像素，根据相关文献的研究成果与实验测试，本文设置4 种体元大小，分别为0.004、0.006、0.008、0.01 m。

根据1.3 节中立体体元模型方法的基本原理，进行DHP 影像模拟，具体流程如图5 所示。先加载所有已配准、去噪的点云数据，计算森林样地点云数据最小包围盒；然后根据最小包围盒与设置的体元大小，构建立体体元，并判断体元内是否包含1 个及1 个以上的点云，剔除不包含点云的体元，将包含点云的体元构建森林场景；其次，根据模拟相机参数，将每个像素生成一条光线，并将每条光线初始颜色赋值为白色，即（0，0，0）；在光线赋值后，采用光线跟踪算法，追踪每条光线是否跟虚拟森林场景是否有相交，对有相交的光线，将其颜色值修改为黑色，即（255，255，255），最后根据每条光线颜色值（即像素值）生成DHP 影像，如图6 所示。

图5 立体体元模型方法模拟DHP影像流程

图6 基于立体体元模型方法模拟的DHP影像示意图

3 结果与讨论

3.1 结果

基于模拟的DHP 影像，采用1.4 节所述方法提取间隙率，采用立体体元模型方法提取的不同体元大小的间隙率结果如图7所示。通过图中可以发现，不同体元大小提取的间隙率整体走势基本一致，立体方法提取的间隙率结果在天顶角0°～75°，呈现逐渐减小的趋势，在天顶角75°～90°时呈上升再下降的趋势。随着体元大小的增大，提取的间隙率也逐渐减小，且差异在小天顶角区域较小，在大天顶角区域较大。

图7 不同体元大小下提取的间隙率

为验证立体体元模型方法的可行性，本文引入点云投影与DHP 两种较为常用的间隙率提取方法进行分析对比，对比结果如图8 所示。通过图中可以发现，三种方法提取的森林冠层间隙率整体走势大致相同，点云投影方法提取的间隙率结果均是随着天顶角的增大，逐渐减小，直到天顶角80°以后，趋于稳定。立体体元模型提取结果与DHP、点云投影方法提取的结果在大天顶角（80°以上）区域存在较大的差异。

图8 不同方法提取的间隙率对比

在Miller 提出的总面积指数计算模型中，主要采用了0°～90°、0°～80°、10°～65°三种天顶角区间提取的间隙［4］，为此本文分别对该三种天顶角区间提取的间隙结果与DHP 方法提取结果进行相关性分析，相关系数如表1、表2 所示。在点云投影方法中，三种天顶角区间下提取的间隙率均与DHP方法有较高的相关性，相关系数均在0.91以上；立体体元模型方法与DHP 方法相比，在三种天顶角不同的区间下，随着体元大小的增加，相关性也不断增加，且当天顶角区间在10°～65°时，4种体元大小下提取的结果相关性均较高，在0.92 以上。而立体体元模型方法与点云投影方法相比时，当考虑整个天顶区域时，随着体元大小的增加，相关性也不断提高；而当在天顶角在0°～80°与10°～65°两个区间段时，具有较高的相关性，相关系数在0.91以上。

表1 与DHP方法结果的相关系数

表2 与点云投影方法结果的相关系数

3.2 分析与讨论

立体体元模型法通过16个测站点云配准后，因各个测站之前遮挡情况不一样，可以有效地弥补各个测站之间因冠层遮挡引起的点云数据缺失问题，使得获取的森林点云数据更为完整。但各个测站之间配准都有一定的误差，特别是多测站存在误差累积，这些误差都会影响间隙率的提取。DHP 方法是目前提取间隙率较为成熟、使用较多的光学方法，该方法主要受DHP 采集时相机曝光度影像，在天顶方向容易曝光过渡，而在大天顶角区域容易曝光不足，从而使得提取的间隙率精度不高。点云投影的方法是通过单站扫描的点云投影生成的DHP 影像，该方法容易受冠层、单树之间的遮挡的影响，同时由于激光在不同的距离采样精度不一样，可能导致获取的森林点云数据缺失较多。当采用立体体元方法时，由图6 可知，当体元大小越小时，在大天顶角区域，由于数据的不完整，导致存在较多的小间隙，使得体元越小时，与DHP 方法的相关性越低，这是由于在激光束角分辨率一定时，随着距离越远，采集的点云数据越稀疏，为此，当体元大小越小时，大天顶角区域存在较多间隙；而当体元大小一定时，随着使用的天顶角区间不一样，提取的间隙率结果与DHP 的相关性也不一致，天顶角在10°～65°时，相关性最高，这是由于在该天顶角区间，两种方法提取的间隙相差较小，且走势基本一致，均为随着天顶角的增加，间隙率逐渐减小。因此，本文使用的立体体元模型法可用于森林冠层间隙提取。

4 结束语

间隙率是森林冠层重要的结构参数。本文以25 m×25 m 的落叶松林为实验样地，在样地内均匀布设了16个采样点进行激光雷达数据采集，从而开展基于激光雷达森林冠层间隙率提取方法研究。本文采用一种立体体元模型方法，基于立体体元模型构建的虚拟森林场景，利用光线跟踪算法与八叉树数据结果来模拟DHP 影像开展间隙率提取，共设计了4 种体元大小（0.004 m、0.006 m、0.008 m、0.01 m）下的虚拟森林场景间隙率，分0°～90°、0°～80°、10°～65°三种天顶角区间与DHP、点云投影两种常用的间隙率提取方法进行分析对比，结果发现立体体元模型法与DHP、点云投影两种方法结果整体趋势基本相同，有较高的相关性，该方法可用于森林冠层的间隙提取。