基于无人机LiDAR的橡胶树单木地上生物量估测

包栎炀　王祥军　李少达　谭骏祥　黄肖　宋佳倩　贾文惠　吴满意

关键词：无人机LiDAR；橡胶树；地上生物量；单木结构参数

中图分类号：S758.1 文献标识码：A

橡胶树（Hevea brasiliensis）原产于南美洲亚马逊河流域的热带雨林，其产出的天然乳胶是一种重要的战略资源与工业原料，因此我国高度重视橡胶树栽培与胶园培育工作，目前我国橡胶树种植面积已达114.53 万hm2[1-2]。橡胶林生态系统是海南重要的人工生态系统，其生物量和生产力受到科研人员的广泛关注[3-5]。橡胶树生物量表示在一定时间内单位空间中橡胶树所产有机质的累积量，是反映橡胶树生产力、固碳能力和碳储量的一个重要指标[6-7]，目前，橡胶树生物量研究主要集中在地上生物量（above-ground biomass，AGB）。另外，橡胶树单木生物量研究有助于动态理解橡胶树特定的生长过程，并从单木尺度上评估胶园的森林生物量构成[8]。目前比较流行的AGB 估测方法是树种异速生长方程，即通过测定易测的单木因子（树高、胸径等）来推算难以测量的AGB[9]，但传统的获取单木因子的方法费时费力，而且由于品种和生境差异，林木生长存在较大的异质性，异速生长方程的泛化能力较低，因此如何快速估测橡胶树AGB 以及橡胶树AGB的动态监测仍存在较大挑战。

随着遥感技术的发展，激光雷达（light detectionand ranging，LiDAR）技术因其对植被空间结构和地形探测能力的优势，广泛应用于林业研究[10-12]。在样地尺度和单木尺度， 无人机LiDAR 和地基LiDAR 也已经成为一种关键的遙感技术，广泛用于小尺度的单木参数提取和AGB估计[13-16]。无人机LiDAR 操作灵活、使用方便，且在冠层结构提取方面优于地基LiDAR，能够快速获取区域尺度的森林AGB 空间分布，而且它所获取的单木因子精度理想，可以大幅提高AGB的估测精度。因此，通过无人机LiDAR 技术获取的橡胶树点云数据提取单木结构参数，进而建立橡胶树单木AGB 估测模型是可行的。

AGB 估测模型常用的预测因子是树高、胸径以及冠幅等，但机载LiDAR 数据很难提取到树木的胸径值[17]，而树冠作为树木进行光合作用和呼吸作用的重要场所，除了直接的单木因子冠幅外，衍生的单木参数如树冠面积和树冠体积等也能在一定程度上反映树冠形态大小，已有研究表明引入树冠投影面积和树冠体积可以提高AGB 模型的拟合精度和预测精度[18]。除了要选择合适的输入变量外，模型的选择也很关键，AGB 估测可分为参数化方法与非参数化方法两大类，参数化方法通常采用统计回归，假设AGB 与预测因子之间可以通过参数预先确定明确的模型表达式，简单且便于计算[19]。AGB 与预测因子之间通常呈现幂、指数或对数等非线性关系，因此非线性回归模型常被用于生物量建模[20]。但由于林木生长的区域差异性，参数化方法得到的AGB 估测模型往往泛化能力不强，模型难以有效和完整表达，现阶段AGB 的估测模型已向非参数化模型拓展[16， 21]。随机森林（random forest， RF）算法是一种抗过拟合的非参数集成建模方法，具有较高的精度和泛化能力，目前广泛用于生物量估测[22]。因此，本研究以海南省儋州市中国热带农业科学院橡胶研究所的橡胶树育种试验林地作为研究对象，获取试验林地的无人机LiDAR 点云数据，在单木分割的基础上，提取树高、冠幅、树冠投影面积、树冠体积等单木结构参数作为输入变量，对比多元非线性回归和随机森林回归2 种方法来估测橡胶树单木AGB 的精度，从而为橡胶树生物量的育种改良和橡胶林生态系统监测提供一定基础。

1 材料与方法

1.1 材料

橡胶树试验林段位于海南省儋州市中国热带农业科学院橡胶研究所育种基地（ 19°11′～19°52′N，108°56′～109°46′E）。儋州市是我国最重要的天然橡胶生产基地，位于海南西北部，属热带季风气候，年均气温23.5～24.1 ℃，年均降雨量1800 mm[23]。实验树于2013 年定植，生长状况良好，林相整齐。

1.2 方法

1.2.1 数据采集 本研究使用MIRACLE 3 无人机系统搭载CBI-120P 激光雷达（成都奥伦达科技有限公司）对橡胶树试验林段进行LiDAR 数据获取。雷达激光波长为905 nm，扫描视场角为360°，扫描频率为5～20 Hz，最大有效扫描速率为每秒64 万点，有效探测距离0～120 m。数据采集时间为2021 年11 月13 日，天气状况良好。由于试验林段橡胶树株行距较小（约1.5 m×3.0 m），林分郁闭度较高，无人机激光雷达很难得到树干的回波点，因此我们对试验林段进行了间伐，以期获得更加完整的橡胶树整株点云。间伐后株行距约3.0 m×6.0 m，林相通透。间伐前后林相如图1 所示，获取的LiDAR 点云数据如图2 所示。

LiDAR 数据获取完成后，随机选择85 株进行AGB 实测，并测定其离地1.5 m 处的茎围。实验树伐倒后，分离树干、树枝和树叶，分别称取鲜重，并对各组织分别取样，3 次重复，在85 ℃恒温下烘至恒重，称取组织样品的干重，计算不同组织含水率，由此计算整株组织的干重。树干、树枝和树叶的干重总和即为橡胶树AGB。

1.2.2 点云数据处理 为了更好地利用LiDAR数据，需要对获取的点云数据进行预处理，主要包括3 个步骤：首先采用离群点检测算法去除孤立的噪点，然后采用布料滤波算法[24]将点云划分为地面点和非地面点，最后利用反距离权重法对地面点进行局部内插处理生成数字高程模型（digital elevation model， DEM），并对非地面点基于DEM 进行高程归一化。

点云数据预处理完成后，即可进行单木点云分割。首先根据点云密度垂直分布曲线对橡胶树林段进行垂直分层，识别林下植被层，然后采用分层DBSCAN（density-based spatial clustering ofapplication with noise）聚类算法提取橡胶树的基本骨架点，并根据骨架点采用自下而上的区域生长算法将橡胶树单木点云分割出来[25]。对于少部分密集林分，由于树冠粘连，会出现单木树冠错分的问题，由此导致后续提取的单木参数值与实际测量值不符，所以对这部分结果进行手动处理，以纯化单木分割结果（图3）。

得到纯化的单木后，将每株橡胶树Z 值最大的点作为单株树木的最高点，最高点的Z 值即为树高。徐志扬等[26]提出了一种树冠角度分区提取冠幅的方法，首先根据点云分层密度（本研究设阈值为80）划分树干与树叶填充部分树冠的高度边界，以树木顶点为中心，计算每个角度分区（本研究设为20°）的最远点求均值的二倍进而得到冠幅。本研究采用此方法，对树冠部分进行投影，采用二维凸包算法[27]求得树冠边缘的凸包点，再以这一系列凸包点的中心点作为分区中心求取单株树木的冠幅，在此基础上，以凸包边缘点构建凸多边形，求取多边形的面积作为树冠投影面积；树冠体积计算采用几何模拟法[28]，将树冠分成1个圆锥体与多个截断圆锥体，对各层采用二维凸包算法求其投影面积，树冠体积等于圆锥体体积和若干个截断圆锥体的体积总和。

1.2.3 地上生物量的多元非线性回归估测 生物量模型最具代表性的是CAR（constant allometricratio）模型[27]，但传统的CAR 模型自变量较少，可能导致估测精度低，而单木参数之间满足相对生长关系，可在CAR 模型的基础上引入相关的单木因子来建立生物量模型[29]。本研究基于最小二乘回归，将提取的单木树高、冠幅、树冠投影面积以及树冠体积作为输入变量建立单木AGB 模型。模型公式如下：

式中，AGB 表示地上生物量；H 表示树高；C 表示冠幅；S 表示树冠投影面积；V 表示树冠体积。a、b、c、d、e 为模型参数。

1.2.4 地上生物量的隨机森林回归估测 随机森林回归是由多个决策树组成的集成模型，广泛应用于分类和回归问题，该模型利用了最小均方误差的原理，能有效解决决策树过拟合和精度低的问题[30]。随机森林采用自助采样法（bootstrapsampling）构建独立决策树的训练子集，这种有放回的抽样使得在原始数据集中有部分数据未被抽中，这些数据被称为袋外数据（out of bag，OOB）。由于OOB 不包含在训练集中，因此可以将其作为决策树的验证集对模型的泛化性能进行评估。另外，由于不同决策树的训练样本不同，从而保证了回归树的差异性。随机森林回归拟合AGB 过程通过使用R 软件的Random Forest 包实现，其中需要确定2 个重要参数：Ntree 和Mtry，Ntree 为决策树的数量，Mtry 为决策树生长分裂时被测试的特征数量。这2 个参数决定着模型的预测能力，根据交叉验证结果选取最优参数，确定Ntree 取值为600，Mtry 取值为4。每棵回归树的叶子结点输出估测值，将所有回归树输出的估测值取平均即可得到最终的橡胶树AGB 估测值。

1.2.5 估测模型的精度评价 为了评价橡胶树AGB 模型的估测精度，本研究采用决定系数（R2）、均方根误差（RMSE）、相对均方根误差（rRMSE）、平均绝对误差（MAE）、总相对误差（TRE）等统计指标来评估模型的拟合精度以及预测能力，公式如下：

1.3 数据处理

采用Alundar Platform v2.4.2 软件直接量测单木点云的树高和冠幅值作为真实参数， 利用Microsoft Visual Studio 2015 进行点云数据预处理、单木分割以及参数提取的算法编写与实现，利用IBM SPSS Statistics 26 、RStudio 以及Microsoft Excel 2019 等软件进行数据处理、分析、模型建立以及图表制作等。

2 结果与分析

2.1 单木茎围分布

橡胶树一般量测胸高处（离地1.5 m 高）的茎围，由于树木的AGB 绝大部分来源于树干，故橡胶树的树围与AGB 之间存在很强的相关关系。结合实测AGB 与实测树围对样地内的橡胶树进行分析，样地内的橡胶树皆为同一时期定植，橡胶树茎围范围为17～63 cm，采用层次聚类构建茎围的层次结构，以欧式距离（euclidean distance）计算变量间的距离，类平均法（average）进行聚类，将茎围聚类为3 个层次，得到聚类的树状图（图4）以及对应的茎围分布信息（表1）。由表1 可知，样地内的橡胶树AGB 随着树围的增长而增加。整体来看，树围在17～30 cm 的橡胶树AGB较小，共有27 个样本；树围在31～49 cm 范围内的橡胶树共47 个样本；树围在50～63 cm 的橡胶树AGB 较大，共11 个样本。

2.2 单木参数提取结果

基于点云数据，求取单木点云的Z 值最大点估测单木树高， 以角度分区法估测冠幅， 以Alundar Platform 软件量测的树高和冠幅值作为实测值，冠幅以东西方向和南北方向量测的平均值作为实测冠幅。分别对估测的树高和冠幅值与实测的树高和冠幅值进行回归分析，散点图见图5。

由图5A 可见， 估测树高与实测树高的Pearson 相关系数（r）为0.999，RMSE 为0.109 m，说明估测树高与实测树高的相关性与一致性极高；由图5B 可见，提取冠幅与实测冠幅的Pearson相关系数（r）为0.951，RMSE 为0.452 m，说明估测冠幅与实测冠幅的相关性与一致性极高，表明LiDAR 点云能够很好地提取单木树高与冠幅，具有较高的可靠性。

分离树冠点与树干点后，基于二维凸包算法对单木树冠点进行投影计算树冠投影面积，并基于改进几何模拟法求得橡胶树单木的树冠体积。样木的树冠投影面积以及体积数据见表2。

2.3 橡胶树生物量估测

2.3.1 橡胶树单木结构参数与地上生物量的相关性分析 将从LiDAR 点云提取的橡胶树单木结构参数与实测AGB 进行Pearson 相关性分析，结果如表3 所示。从表3 可以看出，提取的4 个参数（树高、冠幅、树冠投影面积、树冠体积）与实测AGB 的Pearson 相关系数分别为0.699、0.846、0.855、0.904，相关性均达到极显著水平（P<0.01），说明4 个单木结构参数对AGB 的解释性很好，可用其建立橡胶树AGB 模型。

2.3.2 橡胶树地上生物量估算结果 以基于点云提取的树高、冠幅、树冠投影面积、树冠体积等单木结构参数为自变量，实测AGB 为因变量，分别采用多元非线性回归和随机森林回归2 种方法建立橡胶树AGB 估测模型。2 种模型的五折交叉验证精度评价结果见表4。基于机载LiDAR 点云数据进行AGB 估测可以获得较高的精度，多元非线性回归与随机森林回归2 种模型的R2 分别为0.82 和0.85，均大于0.80，TRE 趋近于0，拟合效果良好。相较而言，随机森林回归拥有更好的拟合效果，其模型拟合优度更高，R2 较多元非线性回归提高了3.64%，预测性能（rRMSE）提升2.66%。采用相同的五折交叉验证策略，即样本划分相同，每次所用的训练样本和验证样本一致时，随机森林回归模型的泛化能力与预测能力也优于多元非线性回归模型，RMSE 减少1.44 kg，MSE减少1.74 kg，rRMSE 减小2.66%。

将2 种模型的估测值与实测值进行比较，得到图6 所示的散点图。由图6 可以看出，2 种模型对于小生物量样本（茎围<50 cm）的拟合优度高于生物量较大的样本（茎围>50 cm），较大生物量样本AGB 估测时的偏差明显增大。这可能由以下3 个原因造成：（1）根据“2.1”茎围与AGB分析结果，由于研究区小生物量的样本较多，在进行五折交叉检验时，选择的训练样本更大几率为小生物量样本，用小生物量样本拟合的生物量模型估算大生物量样本，可能会导致估测结果出现偏高或偏低的情况；（2）由表3 可知，AGB 与冠幅呈极显著正相关，所以小生物量样本往往冠幅较小，树冠交错现象少，单木分割结果准确，提取的参数更准确，估测精度更高。但是，生物量较大样本由于冠幅较大，单木之间易存在树冠交错现象，冠层分割可能不准确，单木参数提取的准确度降低，导致AGB 估测的偏差增大；（3）由表3 可知，AGB 与树冠体积呈极显著正相关，所以生物量较大的样本往往树冠较大，冠层更密，激光受到的遮挡更大，树冠基部点云稀疏，导致树冠体积计算的误差增大，进而导致AGB 估测的偏差增大。

3 讨论

本研究使用无人机激光雷达获取的橡胶树点云数据，对其进行单木分割并提取结构参数来建立橡胶树的单木AGB 模型。传统的模型一般以树高和胸径作为输入参数来构建AGB 的异速生长方程，能得到较高的拟合优度。冯仲科等[31]利用LiDAR 点云数据获取的树高与全站仪测量的树高之间存在显著的相关性，说明使用机载LiDAR 点云数据测量树高是可行的。但是，由于精度缺陷，机载LiDAR 点云数据很难提取到精确的单木胸径，对于该问题，目前解决方法主要有2 种：一是基于统计模型利用提取的单木结构参数来预测胸径[13]，但胸径的估测偏差可能会给AGB 模型带来更大的误差；二是建立不依赖胸径的AGB模型[15， 27]。本研究采取第二种方法解决机载LiDAR 点云难以提取胸径的问题，将树高、冠幅以及衍生的单木参数－树冠投影面积和树冠体积作为输入变量估测橡胶树AGB。结果表明，在没有胸径作为变量构建AGB 模型的前提下，仍然能得到较优的拟合效果，这与刘浩然等[27]的研究结果一致。

传统生物量建模一般选择参数化建模，得到具体的异速生长方程，这能很好地解释单木结构参数与预测变量之间的相关关系。相比于参数化方法，非参数化方法不需要先验假设，模型构建更加便捷，机器学习算法是典型的非参数化方法，尽管不能得出具体模型，但算法的预测精度较高。为充分利用样本信息，文献[32-35]采用留一交叉验证、十折交叉验证等方法来验证AGB 模型的预测精度，得到较高精度的AGB 预测模型，其中卜帆[35]通过留一交叉验证比较不同回归模型的生物量估算精度，得到误差最小的模型作为最优模型估测研究区AGB，有较好的预测效果。因此，本研究采用多元非线性回归和随机森林回归2 种方法估测橡胶树AGB，以五折交叉验证方法对比2 种模型的评价指标，得到随机森林方法构建的AGB 模型泛化能力更强，精度更高，能够更加准确地估算和预测橡胶树AGB。随机森林作为一种优势方法广泛应用于森林AGB 模型的构建[16， 30]，有效解决了决策树过拟合和精度低的问题，适合处理特征高度相关的数据。用于估测林木AGB 的其他常用非参数机器学习方法还包括支持向量机、神经网络、K 最邻近法等。吴宇峰[36]对比了多元线性回归、BP 神经网络和随机森林算法来建立杭州市淳安县的森林生物量模型，结果表明随机森林的反演精度最高，这与本研究结果相同。但是，L?PEZ-SERRANO 等[19]对比了随机森林和支持向量机方法估测墨西哥西北部山脉森林生物量的性能，结果表明支持向量机回归的估测精度更高；邱布布等[37]利用线性回归方法和随机森林法估测杭州森林生物量，结果表明逐步回归的估测精度更优，这与本研究结论不同，说明非参数方法并非在任何情况下都优于参数化方法。所以后续继续对比研究不同模型以选择最优的橡胶树AGB 估测模型是有必要的。

另外，由“2.3.2”模型评价结果可知，本研究所构建的模型对于小生物量样本的估测精度比生物量较大样本好，主要原因在于样本群体中小生物量样本数量较多，大生物量样本较少。未来可以考虑按照径级进行分阶建模，或者加大样本量，构建径级分布更加均衡的建模群体。此外，本研究在模型构建时仅考虑了树高、冠幅、树冠投影面积和树冠体积等4 个单木结构参数，而机器学习往往需要尽可能多的特征来判别，未来可以考虑选择更多的单木结构参数作为机器学习的输入变量以更精确地估算橡胶树AGB。

4 结论

本研究使用无人机激光雷达获取的橡胶树点云数据，对其进行单木分割并提取单木结构参数，将树高、冠幅、树冠投影面积和树冠体积作为输入变量估测橡胶树AGB，结论如下：

（1）基于算法提取的单木树高和冠幅與基于点云人工量测的参数值具有较高的一致性，相关性良好，说明基于无人机LiDAR 点云提取的单木结构参数可以用于单木AGB 估测模型的建立。

（2）基于无人机LiDAR 点云提取的树高、冠幅、树冠投影面积和树冠体积等参数与橡胶树AGB 相关性显著，说明橡胶树的单木结构参数能很好地解释其AGB，可用其作为预测因子估测橡胶树AGB。

（3）采用多元非线性回归和随机森林回归2种方法建立的AGB 模型均具有较高的拟合优度，但是随机森林回归模型的估测精度高于多元非线性回归模型，表明随机森林回归模型的泛化能力更强，能够更加准确估算橡胶树单木AGB。