应用XGBoost算法对森林地上生物量的机载LiDAR反演1)

李洋 彭道黎 袁钰娜

(森林资源和环境管理国家林草局重点实验室(北京林业大学),北京,100083)(珠江水利委员会珠江流域水土保持监测中心站)

森林生物量是评价森林生态系统固碳和碳平衡能力的重要变量,准确估算森林生物量对于研究大面积陆地生态系统的碳循环尤为重要[1]。森林地上生物量占陆地生态系统总量的90%[2],作为森林结构参数的重要组成部分,间接地反映森林的固碳能力,在森林碳储量评估中充当着重要的评价因子。因此,快速、准确地估测森林地上生物量已经成为全球气候变化研究领域日益关注的问题。

国内外对森林地上生物量反演进行了大量的研究,传统的生物量估计通常是根据固定样地连续观测数据,通过异速生长方程预测样地内单木生物量求和得到样地生物量,从而进一步得到更大尺度的森林生物量[3]。但是这种测量方法存在较大的局限性,不仅耗时多、成本高、工作量大,而且对生态系统具有破坏性。随着遥感与地理信息技术的快速发展,利用主被动遥感影像技术(波段信息、植被指数等)与森林地上生物量建立参数化及非参数化模型研究越来越多[4-7]。然而光学被动遥感影像技术易受森林植被物候、天气等影响,并且存在植被信号饱和问题。随着激光雷达技术的不断发展及日趋成熟,森林参数的定量获取取得了显著突破。

激光雷达依据地面采样点激光回波脉冲相对于发射激光主波之间的时间延迟获得传感器到地面采样点的距离[8],获取点云高度、密度、分布、强度及波形信息,从而得到样地的三维结构参数,并且具有极高的测距分辨能力和抗干扰能力。因此,相比于其他遥感技术,激光雷达在林木及林分的垂直、水平信息获取上更加准确便捷。目前,激光雷达技术已被广泛应用于测量植被高度[9]、叶面积指数[10]、郁闭度[11]、地上生物量[12]等林分特征参数[13]。

20世纪80年代中期,激光雷达技术开始应用于估测森林生物量的研究。激光雷达反演生物量的估测方法主要分为传统的统计回归模型(简单线性回归、多元逐步回归等)和机器学习模型(随机森林、神经网络、支持向量机等)。其中传统的统计回归方法常常需满足一定的假设前提,但森林生长数据的连续观测和层次性难以满足以上假设。随着机器学习技术的快速发展,为森林生长收获预估提供了新的途径。机器学习算法能够深度挖掘数据中的有效信息,通过快速处理自变量与因变量的复杂关系建立预测模型,已经广泛应用于遥感和生态领域[14]。李旺等[15]基于激光雷达数据及实测单木结构信息,分别从样地尺度和单木尺度对森林地上生物量进行估算,结果表明两种尺度下模型的预测结果与地面实测值都具有明显的相关性;庞勇等[16]采用高度变量和密度变量,基于多元线性回归方法对3种不同森林类型的生物量进行估测,结果表明3种森林类型估测的决定系数(R2)均在0.8以上。虽然统计回归方法模型公式较为简单、便于计算,但机器学习算法能更好的发挥遥感大数据量的优势,捕捉生物量与各种特征变量之间的复杂非线性关系[17]。Gleason et al.[18]基于激光雷达数据的高度变量,比较了线性混合效应回归、随机森林、支持向量回归和Cubist模型,结果表明在样地水平上支持向量回归的生物量模型精度最高。Cao et al.[19]采用两种遥感数据源,结合5种机器学习算法建立生物量估测模型,结果发现RF算法的模型最优。由此可见,机器学习方法在林业科学领域的研究有广阔的发展前景。

近年来,作为一种新兴的深度机器学习优化算法,极限梯度提升(XGBoost)算法[20]能适应复杂的非线性关系,模型具有更强的并行处理能力,可以有效解决在机器学习回归模型中出现的过拟合问题。黄宇玲等[21]基于逐步回归的XGBoost方法建立森林蓄积量的估测模型;张亦然等[22]基于XGBoost算法进行了草甸地上生物量的反演研究,该算法都表现出更好的效果。然而,该方法在森林地上生物量的反演建模中的研究尚不充分,并且对于XGBoost模型的解释性并未分析。因此,本文利用机载激光雷达数据和样地实测数据,通过多元逐步回归、随机森林、支持向量机和XGBoost算法分别对根河林业局开拉气林场研究区森林地上生物量进行估测比较,并利用样本中每个特征所分配到的沙普利加和解释(SHAP)值对XGBoost模型的可解释性进行分析,进一步验证XGBoost算法在森林地上生物量模型构建中的可行性和适用性。

1 研究区概况

研究区位于地处内蒙古呼伦贝尔盟东北部开拉气林场,隶属根河市林业局(120°12′～122°55′E,50°20′～52°30′N)(见图1)。地处大兴安岭西坡中段,海拔800～1 000 m。该地区属寒温带湿润型森林气候,并具有大陆季风性气候的特征。年降水量400～550 mm,降水多集中在7—8月份,年平均气温-5.3 ℃,结冻期210 d以上,境内遍布永冻层,个别地段30 cm以下即为永冻层。研究区主要的优势树种为兴安落叶松(Larixgmelinii),伴生树种有白桦(Betulaplatyphylla)和山杨(Populusdavidiana)等。

2 研究方法

2.1 样地数据

研究区共建立了125块圆形样地,每块样地半径为13.82 m,面积约为600 m2,其中包括针阔混交林69块、阔叶林30块、针叶林26块。在每个样地的中心点设立明显标志,并采用差分定位方法,准确定位样地中心位置,保证样地横纵坐标的定位精度达到1 m以内,以确保地面样地数据与激光雷达数据的匹配,准确提取遥感信息。地面样地调查时间为2019年9—11月,调查内容包括3大类:(1)林分因子(优势树种组、起源、龄组、郁闭度等);(2)林木因子(样木编号、立木类型、树种、胸径、树高、方位角、水平距等);(3)其他信息(样地号、位置坐标、海拔、坡度、坡位、坡向、调查员、调查日期等)。

根据样地调查数据应用异速生长方程[23]对样地内每株单木的生物量计算,求和得到样地尺度内森林地上生物量,进而换算成单位面积的地上生物量。样地数据基本信息及研究区主要树种异速生长方程见表1、表2。

表1 调查样地数据基本统计结果

表2 调查样地树种异速生长方程

2.2 机载LiDAR数据

机载激光雷达数据在2019年9月6日—10月16日航摄获取,传感器型号为RIEGL VUX-1UAV机载激光扫描仪,扫描仪最大脉冲发射频率为550 kHz,光束发散角为0.5 mrad,光斑直径为50 mm,平均点密度约4点·m-2,平均地面点距离约1 m,测量精度10 mm;搭载平台为中型旋翼无人机,机长2.2 m,翼长3.3 m,飞行高度约100 m,平飞速度70～110 km·h-1。

图1 根河研究区地面样地点分布

本研究应用国产LiDAR360软件对原始点云数据进行处理。首先,根据航迹线和激光发射器搭载平台的参数信息对LiDAR点云进行航带平差和去噪预处理;其次,对去噪处理后的点云数据进行分类,分为地面点(0 m)、低矮植被点(0～1 m)、中等植被点(1～20 m)、高植被点(>20 m)4类;再次,通过不规则三角网算法(TIN)结合地面点插值生成数值高程模型(DEM)[24],对点云进行高度归一化处理,去除地形影响;最后,根据地面样地中心点和地面样地半径对归一化点云进行裁切,得到与125块地面样地对应的LiDAR点云集。

2.3 特征变量提取和筛选

根据激光雷达提取的数据,特征空间可以分为点云和波形两种,其中点云特征数据可以根据点云高度、密度、强度特征等直接进行分析,计算相关统计或指数特征量(不同的高度分位数、平均值、标准差、峰度、偏度等[4,25])。在林业研究中,高度变量是指点云高程值相关的统计参数,直接反映了树木的垂直结构和高度信息;密度变量反映了树木点云的返回密度;强度变量则体现了树木点云返回激光雷达传感器的能量大小。这些指标已经被广泛应用于林业领域中单木及林分结构参数的估算[26-27]。

在本研究中,使用国产LiDAR360软件对高度归一化的点云数据集进行处理,提取得到样地尺度的LiDAR点云特征变量。主要选取包括主要的16个高度变量、5个密度变量和12个强度变量等共33个特征变量[28]作为建模自变量(见表3)。

表3 激光雷达特征变量

续(表3)

许多变量都可以参与构建生物量模型,但这些变量之间往往存在高度相关性,并且与生物量之间的相关性不高[29]。因此,并非所有的变量都可用于生物量模型的构建。为了避免数据冗余,降低模型复杂度,提高模型精度,需要对原始特征变量数据集进行筛选。

在本研究中,对于多元线性回归统计模型,采用皮尔森相关系数对选取的自变量进行筛选,提取与生物量相关极显著的特征因子,再通过逐步回归方法对每一个引入的自变量进行筛选剔除。对于机器学习模型,采用递归特征消除算法(RFE)进行特征选择[30]。RFE算法是一个递归过程,它根据不同的重要性度量对特征变量进行排序。其核心思想是:在每次迭代中都会评估特征的重要性并删除非重要特征,最终得到用于构建模型的最优特征子集[31]。常用的RFE算法有支持向量机-递归特征消除算法(SVM-RFE)和随机森林-递归特征消除算法(RF-RFE)[32-33]。最后,分别利用多元线性回归(MLR)、随机森林(RF)、支持向量机(SVM)和极限梯度提升(XGBoost)算法构建样地生物量估测模型(算法均通过R4.1.1软件实现)。

2.4 生物量模型建立

多元逐步回归模型:多元逐步回归是一种是以向前引入为主,变量可进可出的采用多个自变量对因变量进行预测的统计方法。综合了向前引入法和向后剔除法的优点,可以修正多重共线性,广泛应用于各领域的回归模型[34]。公式如下:y=β0+β1x1+β2x2+…+βnxn+ε。式中:y为因变量(生物量);xi为自变量(激光雷达特征变量),βi为模型参数,i=0、1、2、…、n;ε为误差项。

激光雷达特征变量能间接反映出生物量的变化,当多个自变量和因变量形成线性关系时,即建立生物量与激光雷达变量的多元线性回归模型。

随机森林模型:随机森林(RF)是一种统计学习理论,它是利用自举(bootsrap)重抽样方法从原始样本中抽取多个样本,对每个自举(bootsrap)样本进行决策树建模,然后组合多棵决策树的预测,通过投票得出最终预测结果[35]。在整个采样过程中,有些样本可能会被多次采样,而有些训练数据不会被采样,这部分训练数据称为“袋外”(OOB)数据;OOB数据不参与模型拟合过程,但用于检查模型的泛化。由于随机性可以有效地降低模型方差,RF算法可以实现良好的泛化性和低方差抗性,而无需对决策树进行额外的“修剪”[36]。利用R软件中的随机森林包构建随机森林模型。

支持向量机模型:支持向量机(SVM)算法的基本思想是将样本空间通过非线性映射到一个高维特征空间,使在特征空间中可以应用线性学习机的方法解决样本空间中的高度非线性分类和回归等问题[37]。支持向量回归利用结构风险最小化的同时优化模型的性能和泛化能力,并且能够找到非线性和唯一的解决方法[38]。本研究应用R软件中的e1071包构建SVR模型,核函数选择径向基核函数。

极限梯度提升模型:极限梯度提升(XGBoost)算法是在传统的推进(Boosting)算法基础上引入正则化项,进一步控制模型的复杂度,也是梯度提升(Gradient Boosting)算法的实现,具有运算速度快、效果好、易于调整参数和海量数据处理等优势[20]。并且相比于其他机器学习算法,具有更强的可解释性。XGBoost的核心算法思想是在逐步添加树到模型的过程中,使得模型整体的预测效果提升。其目标是要使得树群的预测值尽量接近真实值,而且有尽量大的泛化能力。XGBoost目标函数分为损失函数和正则化项,损失函数揭示训练误差(即预测值和实测值的误差),正则化定义复杂度,避免过拟合。应用R软件中的极限梯度提升包建立模型。

2.5 模型验证

本研究采用以下评价指标:均方根误差(RMSE)、决定系数(R2)、平均绝对偏差(MAE)。R2值越接近于1,模型的拟合精度越好,均方根误差和平均绝对偏差的值越小,模型的预测精度越高。3个评价指标计算公式如下:

2.6 机器学习模型超参数优选

算法模型的性能由其超参数决定,超参数设定是构建模型的关键。相同的算法,不同的参数,得到的结果具有差异性,参数调整对模型性能有很大影响,优化算法参数能显著提高模型预测精度[39]。采用网格搜索技术对3种机器学习算法进行超参数调优,对训练集数据进行5次10倍交叉验证法得到稳定的模型结果。利用R软件包对SVR、RF和XGBoost算法进行超参数调优(见表4)。

3 结果与分析

3.1 变量筛选

3.1.1 皮尔森相关系数及逐步回归的特征变量选择

由表5可知,应用SPSS19.0软件,对提取的33个激光雷达特征变量和地面样地的森林地上生物量进行皮尔森相关性分析,得到与地上生物量极显著相关(P<0.01)的自变量特征因子共计27个,其中包括16个高度特征变量、5个点云密度特征变量和6个点云强度特征变量。使用R4.1.1软件包,在皮尔森相关性分析的基础上,再通过逐步回归方法,基于最小信息统计量(AIC)为准则进行自变量再优选,最后得到多元线性回归建模的最佳自变量均为高度变量共4个(包括点云高度平均值、高度偏度、点云高度的二次幂平均、点云高度的三次幂平均)。筛选结果剔除了点云密度变量和强度变量,仅保留点云高度变量,样地的点云高度变量与地上生物量存在极显著的相关性(P<0.01),说明点云高度变量是回归分析最佳建模因子。

表4 3种机器学习算法的超参数调整范围

表5 自变量因子与样地地上生物量的Pearson相关系数

3.1.2 递归特征消除法的特征变量选择

由图2可知,在自变量达到4个时,RMSE的值最小。按其重要值排序依次为:点云高度的二次幂平均、点云高度平均值、点云高度的三次幂平均、高度垂直分布的25%分位数。递归特征消除特征变量的筛选结果与Pearson及逐步回归有三个自变量重合。

图2 递归特征消除方法筛选特征变量

3.2 模型精度评价

由表6可知,按照3∶1的比例随机划分为训练集(93块)和测试集(32块),4种模型在训练集和测试集进行拟合,不同模型的训练集和测试集拟合精度不同。在训练集上,RF模型拟合效果最好(RMSE=9.98 t·hm-2,R2=0.93,MAE=5.69 t·hm-2),其次为XGBoost模型(RMSE=10.80 t·hm-2,R2=0.89,MAE=7.24 t·hm-2),而MLR模型(RMSE=15.92 t·hm-2,R2=0.81,MAE=10.58 t·hm-2)和SVM模型(RMSE=16.49 t·hm-2,R2=0.81,MAE=10.15 t·hm-2)拟合精度相近,4种模型在训练集上的表现都很好,R2都在0.8以上,但这并不能说明模型的泛化能力。在测试集上,XGBoost模型的拟合效果最好(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2),其次是SVR模型(RMSE=12.88 t·hm-2,R2=0.69,MAE=9.31 t·hm-2),而MLR模型的精度相对最低(RMSE=13.99 t·hm-2,R2=0.68,MAE=10.21 t·hm-2)。综合上看,XGBoost模型的两种数据集拟合精度都高且RMSE、R2、MAE差距最小,并且在测试集上明显优于MLR、RF和SVR模型,具有最佳的稳定性和泛化能力。

表6 不同模型的精度评价结果

3.3 XGBoost重要值分析及模型解释

同样应用R软件中的xgboost包功能,进一步分析经递归特征消除算法(RFE)筛选后的4种特征变量与模型生物量的相关性和重要性,并利用SHAP值(样本中每个特征所分配到的数值)解释XGBoost模型。SHAP在2017年提出,用于解释XGBoost等“黑箱”模型[40],其基本定义是对于每个预测样本,模型都产生一个预测值,SHAP值是该样本中每个特征所分配到的数值。相比于传统的特征值重要性排序,SHAP值可以进一步反映出每个样本中特征变量的影响力及正负性[41]。

图3是SHAP特征密度散点图,图上的每个点都是一个特征和一个样本的SHAP值,该值代表了这个特征对单个预测的贡献,点的集合代表了特征整体上对预测结果影响的大小和趋势。y轴上的位置从高到低由特征重要性决定,x轴上的位置由SHAP值决定,颜色从浅到深代表特征值从小到大,该图结合特征重要性对影响地上生物量的因素进行分析。由图3可知,参与建模的4种变量重要值排序依次为:点云高度的平均值、点云高度的三次幂平均、点云高度的二次幂平均、点云高度垂直分布的25%分位数的重要值从大到小依次为11.684、7.058、2.284、1.324。SHAP给出的影响地上生物量的最重要特征为点云高度的平均值。并且该特征与地上生物量呈正相关关系,即随点云高度平均高的增加,估测的林分地上生物量越大。总体上看,4种特征变量与生物量基本都呈正相关关系。

Hm为点云高度平均值;Hsq为点云高度的二次幂平均;Htq为点云高度的三次幂平均;H25为点云高度垂直的25%分位数。

4 结论与讨论

本研究以机载LIDAR数据和林分样地调查数据作为数据源,单位地上生物量为研究对象,根据LIDAR数据源提取的特征变量,分别采用MLR、RF、SVR和XGBoost共两类4种算法建立研究区的森林地上生物量反演模型,探讨了XGBoost算法在森林地上生物量反演模型中的适用能力。

总体而言,XGBoost模型和MLR、RF、SVR模型相比,XGBoost在训练集(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2)和测试集(RMSE=12.20 t·hm-2,R2=0.83,MAE=8.30 t·hm-2)的精度评价指标接近,综合表现最优,而其他三种模型在测试集的评价指标相较于训练集都存在一定程度上的降低,说明XGBoost模型具有更强的泛化能力。张亦然等[22]在利用XGBoost算法建立草地地上生物量模型时,发现整体上利用XGBoost算法要优于MLR和RF算法建立的模型。Li et al.[42]使用Landsat 8和Sentinel-A影像估算湖南省亚热带森林生物量的研究结果表明,在3种遥感数据集下XGBoost模型表现均优于MLR和RF模型。而且与以往其他学者在该地区的地上生物量反演的研究成果相比[43],本研究的XGBoost模型精度要明显优于其基于k-NN算法的研究结果。因此,根据XGBoost算法构建的模型精度最高,具有更好的稳定性和准确性。

在筛选建模因子方面,初始的LIDAR特征变量集包含高度、密度和强度等3种类型变量。而两种不同筛选方法都同时剔除了密度变量和强度变量,只保留高度变量参与模型建立,说明剔除的两种变量与生物量的相关性不强,不适于参与构建生物量模型。这是由于变量易受发射功率、范围、入射角、环境参数和目标结构特征的影响[44],导致在不同情况下同一特征的获取值差异较大,难以反映目标的真实特征。例如,机载LiDAR的回波信号不易穿透高郁闭度林分,从而获取的点云数据多分布于冠层表面[45]。另外,筛选得到的高度变量都包含点云高度的二次幂平均(Hsq)、点云高度平均值(Hm)和点云高度的三次幂平均(Htq),证明三种高度变量都与生物量有显著的相关性,点云高度能够很好地反映林分的平均高度信息[27]。

在模型解释性分析方面,本研究基于SHAP值对XGBoost算法构建生物量模型结果进行简单分析。4种特征变量都对地上生物量有积极影响。但由于研究中样本数较少,数据类型单一,尚未显示出SHAP模型对于多源数据和复杂模型关系的解释分析能力。随着多元遥感数据融合与机器学习算法在林业领域中的应用越来越广泛,可以采用SHAP模型对复杂变量及“黑箱”模型进行解释,为林业数据挖掘与相关分析提供新途径。