基于无人机LiDAR 的杉木树冠上部外轮廓模拟与可视化研究

徐志扬，刘浩栋，陈永富，陈巧＊，李华玉,3，王娟,3

(1.中国林业科学研究院资源信息研究所，国家林业和草原局林业遥感与信息技术重点实验室，北京 100091；

2.国家林业和草原局华东调查规划设计院，浙江 杭州 310019；

3.西南林业大学林学院，云南 昆明 650224)

杉木（Cunninghamia lanceolata(Lamb.) Hook.）是我国南方重要用材树种，分布广泛，杉木上冠在杉木全冠占主要地位，决定了杉木的光合作用、生长状况以及冠形特征，因此对杉木树冠上部外轮廓模拟具有重要意义，准确描绘树冠结构可以为树种鉴别提供参考。树冠结构因子包括冠幅、冠长、冠基高、树冠形状等，其中树冠形状可以用树冠轮廓模型表达[1]。国内外进行了很多树冠轮廓模型的研究，有些学者使用规则几何体表达[2-3]，可是很难用简单几何体对自然界各种形态各异的树冠进行特定描述，因此学者们也采用相对复杂的数学模型模拟树冠轮廓。郭艳荣等[4]、吴丹子等[5]采用多项式和指数函数模拟杉木树冠外轮廓。由于树冠上部和下部往往具有不同的几何形状，一些学者认为采用分段形式能够更为精确的描绘树冠外轮廓。Dong等[6]分上冠、下冠、全冠对福建顺昌杉木采用9个模型拟合并选择最佳模型估计树冠形状。高慧淋等[7-9]利用东北林区人工红松（Pinus koraiensisSieb.et Zucc.）、人工长白落叶松（Larix olgensisHenry）、人工樟子松（Pinus sylvestrisvar.mongolicaLitv.）枝解析数据，采用分段抛物线方程、分段幂函数方程、分段单分子式方程、修正Kozak 方程和修正Weibull 方程对树冠轮廓进行模拟，构建3个树种的树冠轮廓非线性混合效应模型。MARSHALL等[2]、卢康宁等[10]、王小明等[11]均将树冠从最大半径位置处分为上部和下部两部分，分别采用线性或非线性形式单独建立方程模拟，取得不错的效果。

树冠轮廓建模变量的测量一般采取树木伐倒再枝解析方式解决，精确度较高，但是效率不够高。主动遥感作为传统遥感的重要补充，为林业调查以及测量提供了新的途径。激光雷达（LiDAR）是一种主动遥感技术，能够穿透树冠获取其三维结构信息，在林分及单木水平上都得到了广泛应用[12-13]。无人机（UAV）作为一种新兴的低空遥感平台，能够灵活、高效地获取遥感数据，无人机和激光雷达的系统集成已成为森林精细化调查的有力支撑[1,14]。利用无人机激光雷达能够快速获取杉木的树冠垂直结构参数，但是，无人机激光雷达在杉木树冠轮廓模拟以及可视化方面的研究却鲜有报道。

本研究采用半自动化的方式，利用无人机激光雷达数据提取杉木单木点云进行杉木树冠上部外轮廓建模，以相对着枝深度为自变量、枝长为因变量，分别采用多项式、幂函数、指数函数建立杉木树冠上部外轮廓模型，并选择最优模型进行可视化，旨在为研究树种识别和单木预测提供参考。

1 研究区概况与数据来源

1.1 研究区概况

研究区位于江西省分宜县境内的中国林业科学研究院亚热带林业实验中心年珠实验林场（27°30′～27°45′ N，114°30′～114°45′ E），属低山丘陵，海拔220～1 092 m，土壤为黄红壤，年均气温16.8℃，日最高气温39.9℃，最低气温−8.3℃，年降水量1 950.9 mm，集中在3—6月，无霜期252 d。主要森林类型为杉木林（Cunninghamia lanceolata(Lamb.) Hook.）、毛竹林（Phyllostachys edulis(Carr.) H.de Lehaie）、阔叶林等。阔叶林中主要树种有樟树（Cinnamomum camphora(L.)Presl.）、大叶锥栗（Castanea henryi(Skan) Rehd.et Wils.）、钩锥（Castanopsis tibetanaHance）、甜槠（Castanopsis eyrei(Champ.ex Benth.)Tutch.）、苦槠（Castanopsis sclerophylla(Lindl.et Paxton) Schottky）、木荷（Schima superbaGardn.et Champ.）、刨花润楠（Machilus pauhoiKanehira）等。

1.2 数据来源

无人机搭载RIEGL VUX-1LR 激光雷达传感器，通过近红外（1 550 nm）激光束和旋转镜330°视场角快速扫描实现数据的高速获取，精度为15 mm，数据采集于2019年植被生长旺季，采用仿地飞行模式，以地形表面为基准设置相对高度，依次在40、70、100、130、160 m 共5个不同高度扫描获取激光雷达数据。此外，无人机搭载Sony ILCE-6000 微单相机采集可见光数据以供后续生成高分辨率正射影像（DOM）、搭载RedEdge M 快照式多光谱相机采集5 波段多光谱数据以供后续生成多光谱影像。使用160 m 航高激光雷达数据，正射影像作为激光雷达数据的辅助参考材料。

地面样地数据的采集也是在2019年植被生长旺季完成，依据研究区范围内森林资源规划设计调查（“二类调查”）数据中的树种信息，在无人机有效飞行范围内设置面积为0.04 hm2的圆形样地，分别记录样地中心点经纬度与3株定位树经纬度，进行每木检尺，记载5 cm 以上树木的每木相对位置、树种、胸径、树高、枝下高、林木分级、枯死情况、东西与南北向冠幅等。选择树种组成为杉木纯林的2个样地（46 号、47 号样地）进行样地的单株杉木提取，样地所在林分为杉木近成熟林。

2 研究方法

2.1 数据预处理

通过数字绿土公司LiDAR360 软件进行激光雷达点云数据预处理。首先去除点云中的噪声点，它一般为孤立点，包括明显高于地表物和低于地表面的激光点，根据绝对高程或阈值去除较为明显的异常点；其次进行点云分类，将点云分为地面点和非地面点，本研究非地面点为森林的激光雷达反射脉冲点；然后，利用Kriging 插值法对地面点进行插值，得到0.5 m 分辨率的数字高程模型（DEM），进而对点云进行归一化，将去噪分类后的点云高程值减去对应DEM 像元值得到点云相对地面点的绝对高度。为减少CHM 孔洞影响，参照相关研究的孔洞填充算法[15-16]，利用LASTools 开源工具对归一化后的点云生成无孔洞的冠层高度模型（Pit-Free CHM）。

对于地面样地数据，需要将属性表示的每木数据转换为矢量数据。首先在ArcGIS 中添加经纬度记录的中心点和定位树位置数据，定义地理坐标系为GCS_WGS_1984，生成经纬度表示的SHP 矢量图层，然后进行投影变换，定义投影坐标系为WGS_1984_UTM_Zone_50N，转换为横纵坐标表示的SHP 矢量数据，以中心点横纵坐标值为原点对样地内每木相对位置进行三角函数解算以得到每木的横纵坐标，最后以横纵坐标为XY值将每木数据添加到ArcGIS 中生成每木矢量位置数据，采用投影变换后的定位树位置对每木位置进行检查控制，从而得到样地内单木矢量数据。

2.2 杉木参数提取

杉木参数提取中，在LiDAR360 软件下完成单木分割，在ArcGIS 下完成单木匹配与补充建模样本，单木点云纯化与参数（树高、冠幅、上部冠长）提取等均采用Python3.6 编程自动化完成。

2.2.1 单木分割与纯化 采用基于冠层高度模型（CHM）种子点的点云分割算法，先利用局部最大值滤波器对生成的无孔洞的冠层高度模型进行单木树顶探测生成种子点，根据种子点完成单木点云分割。

采用选定的2个杉木纯林样地对分割后的单木进行匹配检验。借鉴Reitberger等[17]的匹配原则来匹配分割后的单株树木与实测样木，在ArcGIS 中叠加CHM、检测到的样木、实测样木，经多次对比试验，探测到的树顶与参考树顶距离若小于样地中所有样木平均距离的60%、而且树高差异小于样地中最大的树高20% 的单木为1:1 匹配的样木，若一个参考的树顶对应了多个探测的树顶，则将探测的最小距离的树顶作为1∶1 匹配的样木，然后从匹配的样木中剔除枯死木、下层被压木、非杉木，剩余43株杉木作为单株模拟样本。

为弥补单株杉木样本量不足的缺陷，以DOM和CHM 为底图，结合“二类调查”矢量数据优势树种属性和先验知识，在样地周边以及大片杉木纯林中从单株探测结果矢量中目视选择补充杉木单株样本到150个以上。

对于本研究的密集林分下分割后的单株点云，仍然存在着“欠分割”样本，这将导致提取的单株树冠特征因子异常与无效，因此需要对单株杉木点云再进行纯化。为有效提取树冠的特征因子，引入Kim等[18]单株分割点云纯化思路并改进优化，首先，以单株杉木树顶为中心，在水平方向上，按照90°间隔将单株点云分为4个分区，各分区的点云再按照0.5 m 间隔从树木中心向外逐层分段缓冲，统计在垂直方向上各分区、各段点云的平均高度；其次，以实地调查的单株杉木树高与冠幅数据建立树高-冠幅模型（CW=HT0.3446−1.0645，R2为0.458 3，其中CW为冠幅，HT为树高），利用该模型结合单木CHM 最大值得到一个冠幅预测值，再将该值乘以1.25 倍系数，以此对各分区的分段缓冲进行冠幅约束。正确分割的单株杉木点云各分区中逐层缓冲的高度值从树木中心点向树冠边缘逐渐减小，而当单株点云“欠分割”时，就会发生平均高度减小到一定水平（拐点）然后增大的现象，因此，“欠分割”单株点云的该拐点处应当是可识别的树冠边缘，故舍弃拐点处以外的点云，保留各段内树木中心到拐点分区的点云，形成纯化后的目标单株点云样本。单木点云纯化原理如图1 所示（XY顶视图）。

图1 单木点云纯化示意图Fig.1 Schematic diagram of individual tree point clouds purification

2.2.2 树高提取 单株杉木点云纯化后完整的单木点云为树冠点与树干点的非地面点集合，单株点云的最高点即为树冠顶点，树冠顶点高度就是单木树高。为减少噪声点对树冠顶点的干扰，使用单株杉木CHM 最大值对单株点云的最高点进行检查，二者差异绝对值大于固定阈值则对单株点云输出二维X-Z图像进行目视检查，确定是否存在噪声点，本研究此阈值设置为0.5 m。

2.2.3 冠幅提取 参考Reitberger等[17]的分层处置方法，模拟冠幅的定义提取单木冠幅。在垂直方向上，对纯化后的完整单株点云，自下而上以0.5 m为间隔分层，统计每层点云数量与占整株点云数的百分比；同时，各分层在水平方向上，以树冠顶点垂直投影为中心，按照5°间隔对点云进行分区，共得到72个分区，求取分区中点云距离树木中心的最大值作为分区的树冠半径，取各分区树冠半径均值的2 倍为该层树冠平均直径。经多次实验比较，取自下而上点云累计数量占整株点云数2%以上、高度5 m 以上的最大树冠直径为冠幅。如图2 所示，(a)对单木点云进行分层，(b)为各层点云数量占单木总点云数量比例，(c)为对各分层点云进行分区求取树冠平均直径。

图2 冠幅提取示意图Fig.2 Schematic diagram of crown width extraction

2.2.4 上部冠长提取 从树冠顶点到最大树冠位置的部分即为树冠上部，单木树高减去最大树冠半径所在高度就可得到上部冠长（CLU）。

2.3 树冠上部轮廓模拟

构建树冠上部外轮廓模型需要获取单株树冠的外部轮廓点。本研究参考相关研究采用宽度分位数表示树冠外轮廓的方法[7]，对纯化后的单株树冠点云自最大树冠位置处向上按照0.5 m 间隔分层提取，分别统计各层点云宽度百分位数85%、90%、95%、98%当作相应层外轮廓点，目视对比树冠上部外轮廓模型精度，最终选择95%分位数建立树冠外轮廓模型。

将处理后的单株样木按照4:1 随机分成建模样本与检验样本，并且对建模和检验的树冠外轮廓点样本分别按照相对着枝深度分层，使用3 倍标准差法剔除外部树冠半径异常的树冠上部外轮廓点。以相对着枝深度（RDINC）和外部树冠半径（OR）建立模型，相对着枝深度为树冠轮廓点到树冠顶点在垂直方向上的相对距离，等于着枝深度与上部冠长之比，外部树冠半径等于轮廓点到树冠中心的水平距离。建模变量模型参数见图3 所示。

图3 单木树冠上部外轮廓模型参数示意图Fig.3 Schematic diagram of individual tree upper crown profile model’s parameters

构建模型时，选取树冠轮廓模拟较为常用的二次抛物线、幂函数和对数函数作为基础模型[4-6,19]，方程形式如下：

采用决定系数（R2）和均方根误差（RMSE）评价拟合优度，使用平均偏差（ME）和平均绝对偏差（MAE）检验模型。

树冠上部轮廓模拟时，外轮廓点参数提取、单木样本随机分配和异常轮廓点剔除、模型建立均采用Python3.6 编程完成。

3 结果与分析

3.1 单木探测及空间格局

两个杉木林圆形样地中，46 号样地31株树木中无人机激光雷达探测到23株全部为杉木，47 号样地23株树木中探测出20株全部为杉木，两个样地探测率分别为74.19% 和86.96%，综合探测率79.63%，如图4 所示，可以看出样地中杉木水平分布并不均匀。图4 中的单株树冠边界显示的是提取后的冠幅。

图4 无人机激光雷达单木探测结果Fig.4 Individual tree detection of UAV-LiDAR

3.2 单木树冠特征因子提取结果比较

图5 所示的是43株检测到的杉木样本提取树高与实测树高的散点图与回归关系，回归相关系数R2为0.890 5，单木树高提取精度均值达到95.25%，表明无人机激光雷达能够获得与地面实测相关性较好且精度较高的预测树高值。

图5 提取树高与实测树高的回归关系Fig.5 Regression equation of inversed tree height and measured tree height

图6 所示的是43株检测到的杉木样本提取冠幅与实测冠幅的散点图与回归关系，回归相关系数R2为0.845 6，单木冠幅提取精度均值达到94.32%，表明本研究单木冠幅提取算法较好的反映了树冠真实大小和分布。

图6 提取冠幅与实测冠幅的回归关系Fig.6 Regression equation of the inversed crown and measured crown

3.3 树冠上部外轮廓拟合模型

采用95% 分位数方法提取出141株建模杉木样本的2 538个树冠上部外轮廓点，再按照3 倍标准差法剔除异常外轮廓点得到2 152个树冠上部外轮廓点，据此拟合出了3个杉木树冠上部外轮廓模型，拟合结果如图7 所示。

图7 无人机激光雷达拟合杉木树冠上部轮廓Fig.7 The fitted upper canopy profile of C.lanceolata with UAV-LiDAR

由表1 可看出，参数估计值标准误均较低，可见参数估计值较为稳定。建模样本的R2均达到0.8 以上，二次抛物线和对数函数拟合R2较近，幂函数拟合R2明显高于二者，达到0.817 0，且幂函数的均方根误差（RMSE）均低于二者，从验证样本的平均偏差（ME）和平均绝对偏差（MAE）来看，幂函数亦明显优于二次抛物线和对数函数，因此本研究中幂函数拟合效果最优。表1 中模型参数估计值、建模评价与验证指标编程得到的结果与SPSS 19.0 所得结果一致，估计值标准误从SPSS19.0 得到。

表1 无人机激光雷达杉木树冠上部轮廓拟合与模型检验Table 1 Upper canopy profile model fitting and validating of C.lanceolata with UAV-LiDAR

3.4 树冠上部三维可视化模拟

在Python3.6 下编程读入已提取的单株点云的树高、最大树冠半径和上部冠长，将单木参数代入幂函数外轮廓模型后可以生成一组高度和树冠半径数据，再以杉木中心为原点将树冠半径转入X-Y平面中进行三角函数解算，将外轮廓二维数据转换为三维模拟数据，加载样地内所有样木的三维模拟数据，调用matplotlib 库对两个建模杉木林样地进行样地尺度的三维重建，为减少地理横、纵坐标大量的数字冗余显示，在三维坐标系中保持Z轴为树木高度而在水平面上对样地中样木地理位置的横、纵坐标进行整体平移，在样地范围内将地理坐标原点移动到样地横、纵坐标最小值位置，并用X-Y轴分别表示相对地理横、纵坐标，结果如图8 所示。

图8 2个样地树冠上部三维模型重建Fig.8 Three dimensions model reconstruction of 2 plot’s tree upper canopy

4 讨论

在单木提取方面，本研究2个圆形样地54株树木检测出1∶1 匹配的样本43株（全为杉木），综合探测率为79.63%。与同样采用无人机激光雷达技术进行单木探测结果比较，Wallace等[20]分别采用5种提取算法对4年生桉树林进行单木提取，提取精度均达到90%以上，最高为98.00%，与之相比本研究探测率明显偏低，这除了与点云密度有关外，还与提取对象的林分特征有关，桉树林中树木在水平方向上更为均匀且垂直方向上更为整齐，桉树林的这种林分特征能够非常明显的提高单木提取精度，而本研究的杉木林单株树木分布不均匀，甚至部分树木群团聚集在一起，46 号样地因此现象有5株树木未被检测出，此外，杉木近成熟林中部分杉木树冠衰老、枝叶凋落造成点云质量不佳也影响了单木检测结果；在国内，文献[1]基于CHM和标记控制分水岭法对东北落叶松林进行单木提取，两个样地323株共计检测出1∶1 匹配的树木235株，综合精度为72.76%，本研究的杉木林较落叶松林更为复杂但结果仍然略高于其提取精度，原因主要是本研究采用基于CHM种子点的单木点云分割方法，其次应与样地大小（样本量）有关，样本较少时得出的探测率意义并不强，这也是下一步工作中需要注意的地方。

本研究对象为密集林分的单株对象，单木点云分割后仍然存在“欠分割”现象，为准确提取单株树冠的建模特征因子，研究中改进Kim[18]的单株点云纯化思想，利用理想状态的单木树冠轮廓点高度值呈现出从树木中心点向外逐渐减小的规律，并结合树高-冠幅模型进行树冠点云约束，对单木点云进行筛选纯化。因杉木为针叶树，冠形规律性较强，该算法效果较为明显，但是，对于其他树种则不一定适用。

本研究中3个模型的建模拟合R2均达到了0.80 以上，较文献[1]基于无人机激光雷达的0.752要高，这与对树冠上部外轮廓点样本按照相对着枝深度分层进行3 倍标准差法剔除外部树冠半径异常点所起作用有关；与郭艳荣等[4]基于实测值建模的杉木中龄林拟合R2结果（0.80～0.85）较为接近，略高于吴丹子等[5]基于杉木实测值建模拟合的R2（0.76～0.80），表明在特定树种和林分条件下，采用无人机激光雷达模拟树冠轮廓，理论上可以达到实测值的拟合效果。

本研究仅对树冠上部进行外轮廓模拟，一是为基于无人机多源遥感数据树种分类识别提供树木冠形方面的参考，二是亚热带密集林分条件下获取的树冠上部点云尚能够反映树木的冠形状态，而树冠下部点云因林分条件影响较为稀疏或者对树木冠形描绘无任何作用。

大多数树冠轮廓模型与可视化以实测的树木枝干解析数据为样本材料，往往通过选取平均木伐倒进行枝条测量，精度较高但是效率有限，且人力物力耗费较大，而无人机激光雷达能够提高工作效率，本研究也表明了利用无人机激光雷达表达杉木树冠上部外轮廓冠型结构的可行性。

本研究杉木人工同龄林为近成熟林，而同一树种在不同龄组的树冠形态并不完全一致，今后可以分别龄组开展研究。

5 结论

本研究中单株树木点云外轮廓参数提取是关键，采用单木点云分割与纯化相结合得到纯化后的单木点云数据，进而提取外轮廓参数，程序实现便捷，但适用对象有限，因此研究在高密度林分条件下单木点云有效提取纯化方法对于客观描绘树冠三维形状具有重要价值。基于无人机激光雷达建立的3个杉木树冠上部模型均取得较好的拟合效果，而幂函数模型表现更优，模拟结果反映了杉木的树冠上部形态，它为杉木的分类识别提供了参考依据。