Estimation model of Eucalyptus stand volume based on airborne LiDAR Point Cloud

研究区位于广西壮族自治区南宁市兴宁区的国有高峰林场，22°51′~23°02′N，108°06′~108°31′E，该区地处亚热带地区，年平均气温约21 ℃，年平均降水量为1 200~1 500 mm，相对湿度为79%，属丘陵地貌，海拔为100~460 m，坡度为6°~35°，具有较厚的赤红壤，适宜亚热带和热带树种生长，森林覆盖率达87%，主要树种为杉木Cunninghamia lanceolata、巨尾桉Eucalyptus grandis × E. urophylla、马尾松Pinus massoniana等。

调查时间为2018年1—2月，研究区内共设置71块桉树样地(图1)，其中激光雷达覆盖范围内共57块样地。样地大小为20 m×20 m、25 m×25 m和25 m×50 m。采用实时动态差分技术(real-time Kinematic，RTK)进行样地定位，记录样地中心点及样地角点。采用每木检尺的方法，使用胸径尺、激光测高仪和皮尺逐一测量样地内树木的胸径、树高等数据。统计样地调查数据得到样地算数平均树高、算数平均胸径、样地面积(表1)。通过广西地区桉树二元材积表对单木材积量进行计算^[12]，进而计算得到样地尺度的公顷蓄积量值(V_样地)：V_样地=V_公顷S_样地/10 000。其中：V_公顷为通过二元材积公式计算得的每公顷蓄积量，S_样地为桉树样地面积。

Figure 1.  Location of the sample plots distribution

于2018年1月采集机载LiDAR数据，使用有人机并搭载LMS-Q680i激光雷达扫描仪，实际飞行高度为1 000 m，最大扫描角度为30°，波长为1 550 nm，激光脉冲长度为3 ns，采样间隔为1 ns，最大扫描频率为400 KHz，垂直分辨率为0.15 m，点云密度为10 pt·m⁻²。

使用LiDAR 360软件对机载激光雷达点云进行点云拼接、去噪、地面点分类、基于地面点归一化等预处理，结果如图2所示。根据样地角点以中心点坐标对点云进行裁剪，提取出基于样地尺度特征变量共48个，包括37个点云高度参数、10个密度参数、郁闭度等点云特征变量，具体参数见表2。

Figure 2.  Example of point cloud preprocessing results of sample plots

使用MATLAB激光雷达覆盖范围内57个样本进行随机抽样，按照3∶1的比例选取42个作为训练样本，15个作为验证样本。以样地实测公顷蓄积量为因变量，筛选后点云特征为自变量，采用逐步回归、偏最小二乘回归、随机森林回归、支持向量回归模型进行拟合 。

对于逐步回归，采用逐步筛选法对所有特征变量进行变量筛选；再对筛选出特征变量进行多重共线性检验，计算各变量间容忍度或特征变量间方差膨胀因子(F_VI) 2个统计量^[13]，对于容忍度≤0.2或F_VI≥5的变量进行进一步讨论，确定最优特征变量子集。F_VI=1/(1−R²)。其中：R为特征变量间相关系数。

对于偏最小二乘回归、支持向量机回归、随机森林回归等方法，使用随机森林中重要性排序对特征变量进行优选。其主要原理为随机森林算法在构建各决策树时，对某一特征变量进行取舍，若此时均方误差(mean square error, MSE)有较大变化，则该特征变量重要性高，最终得到所有特征变量重要性排序。

逐步回归(stepwise regression, SR)可用于筛选并剔除引起多重共线性的变量，逐步回归建立模型一般形式为：Y=β₀+β₁X₁+β₂X₂$+\cdots + $β_iX_i+ε。逐步筛选法结合了向前选择变量法和向后选择变量法的优点，对i个自变量X分别与因变量Y建立一元回归模型，计算各变量所对应F值，其中，β₀为常数，β_i为回归系数，ε为误差项。在建立逐步回归模型时，选择当前未加入模型的预测变量中F的最大值所对应的$ {X}_{i} $加入模型，再对已选入预测变量逐个进行t检验，若存在已选入预测变量不再显著，则将其剔除。重复以上步骤，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止，将保留下来的最优子集作为建立模型的特征变量，并将候选变量对应系数保留。

偏最小二乘回归(partial least squares regression, PLSR)结合了主成分分析、典型相关分析的优点，实现了数据降维、信息综合^[14]。首先对训练样本进行数据标准化处理，调用plsregress函数提取主成分与自变量和因变量组合。分别计算自变量和因变量提取出成分的贡献率进而计算得到累计贡献率，当度量因子Q²_h≥0.097 5停止主成分的提取^[15]。统计主成分对个数并建立各自变量与因变量的线性表达式，最后根据所建立的各主成分对模型，整理得到PLSR模型。

支持向量机算法(support vector machine, SVM)利用内积核函数代替高维空间的非线性映射, 不涉及大数定律和概率测度等有关问题, 且SVM的决策函数仅有少数支持向量决定，该算法在解决小样本多维度回归和分类问题时泛化能力强，具备较好的“鲁棒性”。VAPNIK^[16]在SVM分类的基础上引入了不敏感损失函数，得到了支持向量机回归算法(support vector regression, SVR)。

使用Libsvm工具箱实现SVR算法，采用网格搜索法(grid search)对常用的4种核函数进行参数寻优^[17]，即线性、多项式、RBF (径向基核函数)、sigmoid (多层感知机核函数)。同时进行十折交叉验证，保证回归模型中惩罚系数(C)与gamma值(g)达到最优。

随机森林回归(RFR)算法采用自助采样法(bootstrap sampling)，在以决策树构建Bagging集成的基础上，对样本和特征变量进行随机选择^[18]。调用TreeBagger函数进行RFR算法建模，通过对决策树数量(n_tree)和最小叶子点数进行参数寻优，直到袋外(out-of-bag，OOB)误差的MSE达到最小，保证模型预测性能达到最优，将寻优结果作为RFR模型的建模参数，用于模型构建。

本研究使用决定系数(R²)、均方根误差(RMSE)和平均绝对误差(MAE)对林分蓄积量估测模型进行定量的精度验证和评价。其计算公式为：

其中：$ {y}_{i} $、$ {\bar{y}}_{i} $、$ {\hat{y}}_{i} $分别为实测蓄积量值、实测蓄积量均值、模型预测蓄积量值；E_RMS为均方根误差；E_MA为平均绝对误差；n为验证样本数量。进行评估模型时，R²越趋近于1，代表拟合程度越高；E_RMS越小，代表真实值与模型预测值之间离散程度越小；E_MA越小，代表真实值与模型预测值之间的误差越小。

在采用逐步筛选法进行偏F检验时，取偏F检验拒绝域的临界值为F_进＞F_出，F_进为选入变量时的临界值、F_出为删除变量时的临界值。本研究设定F_进为0.10，F_出为0.11，筛选特征变量结果为D₉、H_P95、H_max、H_kurtosis。

对筛选出的变量进行双变量相关性分析，根据相关系数(R)计算各变量间方差膨胀因子(F_VI)。计算结果如表3，H_max与H_P95之间F_VI为250.25，F_VI大于5，认为两者之间存在共线性。因此分别以D₉、H_max、H_kurtosis和D₉、H_P95、H_kurtos为建模因子。

对所有的48个特征变量进行重要性排序，调用随机森林中OOBPermutedVarDeltaError参数，得到所有特征变量的重要性。选择重要性大于0.2的变量作为建模因子，筛选结果如图3，包括8个高度变量、2个密度变量，其中H_curt重要性最高，达0.52；D₆重要性相对最低，为0.22。

Figure 3.  Importance ranking of point cloud characteristic variables

根据特征变量优选结果，分别建立以D₉、H_max、H_kurtosis和D₉、H_P95、H_kurtos为建模因子的多元线性模型:

使用SPSS 26对Y₁、Y₂模型进行初步评估，计算2种模型的相关系数(R)、R²以及标准估算误差(SE)，结果如表4。结果表明Y₂中R²、SE均优于Y₁，因此选定Y₂作为逐步回归模型。

利用PLSR提取的主成分F₁、F₂与因变量Y关系分别为F₁=5.840Y+b和F₂=0.271Y+b，b为常数，得到最小二乘回归模型为：

调用meshgrid函数对C、gamma进行参数寻优并采用十折交叉验证，调用svmtrain函数，分别构建4种不同核函数的SVR模型并进行模型训练，训练结果如表5。得到拟合结果最优模型为RBF-SVR，其C为8，gamma为0.125。RBF-SVR模型R²为0.85，RMSE为29.24 m³·hm⁻²，MAE为94.98 m³·hm⁻²，选定该模型作为本研究的SVR模型。

本研究采用穷举法对决策树数量和最小叶子点数进行参数寻优。分别设置决策树数量和最小叶子点数最小值为50和2，每次递增10和1，对寻优过程循环，通过观察OOB误差的MSE变化，直到寻找到本模型最优参数，寻优结果如图4。最终确定决策树数量为90，叶子节点数量为2，将寻优结果作为建模参数，输入训练样本，完成RFR模型构建。

Figure 4.  Evaluation of growing stock volume inversion model

使用验证样本对各模型的预测性能进行评估，结果如表6、图5。可以看出不同的回归模型均表现出较好的拟合结果，其中拟合结果最优模型为RFR模型，模型评价结果R²为0.95，RMSE为12.64 m³·hm⁻²，MAE为8.00 m³·hm⁻²，RBF-SVR模型其次，R²为0.94，RMSE为13.09 m³·hm⁻²，MAE为11.65 m³·hm⁻²。将验证样本带入模型，检验模型的预测能力，结果显示各模型预测能力与模型的拟合效果一致。通过传统方法的划分训练样本与验证样本，存在一定的偶然性与不确定性，为了进一步确保各模型的稳定性与泛化能力，采用留一法交叉验证(leave-one-out cross validation，LOOCV)对本研究中各模型进行再次评估^[19]，如图6所示：RBF-SVR模型表现最优，R²为0.88，RMSE为21.35 m³·hm⁻²，MAE为16.62 m³·hm⁻²，与其他模型相比，R²高出0.03~0.07，RMSE减少2.58~5.17 m³·hm⁻²，MAE减少0.79~3.32 m³·hm⁻²。

Figure 5.  Scatter diagram of training results and verification results of each model

Figure 6.  Ten-fold cross validation results of each model

本研究采用SR、PLSR等参数回归方法与RFR、SVR等非参数回归方法进行广西高峰林场桉树人工林的蓄积量反演模型研建，通过对点云特征变量筛选、建模因子优选、参数寻优等方法保证各模型达到最优性能，在此基础上采用留一法对各模型进行交叉验证，保证了模型的稳定性与泛化能力，并将该结果作为本研究各模型性能评估的最终结果。研究结果如下：①核函数为RBF的SVR在4种模型中预测精度最高(ΔR²为0.03~0.07、ΔRMSE为2.58~5.17 m³·hm⁻²、ΔMAE为0.79~3.32 m³·hm⁻²)，采用SVR模型在解决林业激光雷达领域的回归预测问题已有广泛应用，肖越^[20]结合Landsat 8、高分2号(GF-2)共2种光学遥感数据与极化SAR数据进行旺业甸林场尺度的蓄积量反演，采用SVR模型拟合精度最高(R²为0.48，RMSE为57.27 m³·hm⁻²)，与MLR、RFR等方法相比，ΔR²为0.10~0.11、ΔRMSE为5.30~6.00 m³·hm⁻²。赵勋等^[21]基于机载激光雷达进行广西高峰林场林分平均树高估测中，通过随机森林特征变量筛选，采用SVR相比较RFR方法，R²增加0.01、RMSE减少0.06 m。以上研究结果与本研究一致，进一步表明解决基于LiDAR点云特征变量的林分蓄积量反演问题时，结合随机森林筛选特征变量与支持向量机回归可作为有效建模方法。

②本研究采用2种方法对特征变量进行筛选，逐步回归法保证入选的特征变量显著且各变量之间不存在共线关系，随机森林则通过计算各特征变量对建模时的贡献度排序进而筛选变量。肖越^[20]利用逐步回归法从遥感特征变量中筛选出建模因子用于蓄积量回归模型构建；周蓉等^[22]在基于Landsat 8 遥感影像反演地上生物量的研究中，采用随机森林重要性排序方法从遥感因子中筛选出特征变量用于构建模型。这2种方法均广泛应用于构建模型前的变量筛选，合理地选择模型筛选方法，更有利于回归模型构建。

③本研究旨在选取最优桉树蓄积量反演模型，参数方法SR中逐步筛选法已对变量进行筛选，故SR模型不再采用随机森林进行特征变量筛选。在本研究模型评估中，SR模型评估结果均优于PLSR模型，袁钰娜等^[23]对东北林区4种不同针叶林蓄积量反演所建立的回归模型中，PLSR模型的拟合结果和预测精度均优于SR (ΔR²为0.05~0.15，ΔRMSE为2.6%~4.2%)，与本研究结果不一致，但由于2种回归方法并未使用同一变量筛选方法，因此，不能直接认为在解决预测回归问题时前者更加可靠，对于以上2种回归方法的模型性能与预测能力需进一步讨论。

④本研究在林场尺度通过较少样本进行桉树林分蓄积量反演，非参数方法总体预测精度较高且RBF-SVR模型略优于RFR模型。SVR核心思想为将低维空间中的向量用非线性函数映射到一个高维特征空间，进而寻求线性回归超平面并解决低维空间中的非线性问题^[24]，支持向量机回归解决多维度小样本回归预测问题相较随机森林回归更有优势。在反演尺度较大且训练样本较多的林分蓄积量反演问题时，采用支持向量机回归方法训练效率将明显下降，此时选择随机森林回归将更适用。

⑤本研究采用逐步筛选法筛选出的4个特征变量中，点云密度变量被选择1次；在采用随机森林筛选的10个变量中，点云密度变量被选择2次，且密度变量D₉在重要性排序中位于第3位。点云密度变量作为本研究中重要解释变量之一，该类型变量能够描述林分中树木的水平结构信息，与树高(垂直结构变量)等参数结合可解释林分空间结构信息。孙忠秋等^[25]加入林分郁闭度作为水平结构的解释变量，与仅用点云高度参数相比较，RMSE下降0.27 m³·hm⁻²、MAE下降0.08 m³·hm⁻²，该变量对模型预测性能有一定提升。本研究将郁闭度作为候选变量，但在逐步筛选法中又将其剔除，因为在基于随机森林法的重要性排序中，林分郁闭度参数在48个候选变量中排名第44位。由此可见，在林业激光雷达应用中，将林分密度用于衡量林分水平结构信息更加可靠。

本研究通过2种参数回归方法(逐步回归、偏最小二乘回归)与2种非参数方法(随机森林回归、支持向量机回归)对林分蓄积量进行估测。在采用传统方法划分训练样本与验证样本的结果中，非参数方法模型精度均优于参数回归方法；采用留一法对各模型预测精度进行交叉验证时，表现最优模型RBF-SVR属于非参数回归方法，表明解决林业激光雷达领域中的回归预测问题时，非参数方法相较参数方法更有优势。本研究中蓄积量反演模型已在广西高峰林场内取得较好模型精度，但对于大尺度反演森林蓄积量的适用性有待进一步论证。