运用Sentinel-2遥感影像数据估测森林蓄积量1)

李坤 吴达胜 方陆明 刘建超

(浙江省林业智能监测与信息技术研究重点实验室(浙江农林大学)，杭州，311300) (浙江农林大学暨阳学院)

森林是陆地生态系统的重要组成部分，在全球碳循环中起着至关重要的作用，全面快速地掌握森林信息对国民经济和绿色可持续发展有着积极作用。森林蓄积量是衡量森林质量的一个重要因子，也是反映森林资源水平的一个基本指标，精准估测森林蓄积量有助于监测森林质量，为林业资源管理决策提供科学依据。传统的森林资源调查不仅费时费力，而且难以反映森林的动态发展情况。因此，林业工作者及研究人员一直致力于探索如何利用飞速发展的遥感技术对蓄积量进行估测[1-2]。

目前，遥感技术应用于森林蓄积量估测主要体现在三方面：一是选取光学遥感影像，对其进行主成分分析后选取主成分并提取纹理特征，探讨不同分辨率不同窗口大小提取到的纹理特征与森林蓄积量的相关性；二是结合不同遥感数据源(光学遥感、雷达遥感)构建估测模型；三是使用不同建模方法，从回归模型发展至人工神经网络、k最近邻算法等经典机器学习模型[3]。刘明艳等[4]利用主成分分析法对老秃顶子自然保护区进行蓄积量估测，模型预估精度达到92.18%；杨明星等[5]利用多窗口遥感影像纹理特征对蓄积量进行建模评估，使用随机森林模型对蓄积量进行建模评估，精度达到75.46%，估测效果较好；蒋馥根等[6]基于方差优化K近邻法对蓄积量进行反演，该模型决定系数(R2)达到0.69，相对均方根误差(RRMSE)为32.04%。机器学习模型为非参数模型，相较于传统的经验模型，不仅能够克服特征因子之间的共线性问题，而且具有更好的泛化能力，能够对类别特征进行处理，更适用于处理数据量大、特征繁多的林业数据。

国内外学者通常运用当地的遥感影像和外业调查数据对当前区域的蓄积量进行预测或者反演[7-9]，很少将同一模型应用于不同区域，检验模型的普适性；在选取特征建模时，由于传统经验模型不能直接处理类别特征，需要根据类别特征将数据划分为不同的类别，再根据不同类别分别建立不同的模型[10-11]，这样的模型建立繁琐，且适用性不强。

本研究从非参数模型泛化性能较强的特点出发[6，12]，结合机器学习算法处理类别特征的能力[13]，提出最小绝对收缩和选择算子(Lasso)特征选择方法的堆叠(Stacking)融合模型。根据Sentinel-2遥感影像、二类调查数据与数字高程模型(DEM)数据，首先选择K最近邻算法(K-NN)、梯度提升迭代决策树(GBDT)、极端梯度提升(XGBoost)、梯度增强集成分类器(CatBoost)4种单一模型进行建模估测，然后探讨类别特征的加入对模型指标的影响，最后建立堆叠(Stacking)融合模型。采用十折交叉验证法检验模型性能指标，比较融合模型与单一模型的精度，验证不同模型的普适性，探讨类别特征对模型优化的影响，为森林蓄积量估测中特征选取以及模型构建提供新的思路。

1 研究区概况

淳安县位于浙江省西部，地处北纬29°11′～30°2′，东经118°20′～119°20′，陆域面积4 417.48 km2，属中亚热带季风气候。全县森林覆盖率达65.00%，生物种类繁多、资源丰富，年平均温度为17 ℃，年平均降水量1 430 mm，地形以丘陵山地为主，海拔50～1 500 m，土壤以红壤、黄壤、岩性土、水稻土为主，其中红壤为面积最大的土属。

临海市位于浙江省东南沿海，地处北纬28°40′～29°4′，东经120°49′～121°41′，陆地总面积2 203 km2，属亚热带季风气候，年平均气温17 ℃，年平均降水量1 550 mm，四季分明，适合植被生长，境内森林覆盖率为64.20%。地形以丘陵和山地为主，山地面积占总面积的70.7%，平原面积占22.8%，呈“七山一水两分田”地貌。

2 研究方法

2.1 遥感影像数据和地形数据

本研究使用Sentinel-2遥感卫星影像，成像时间为2017年11月27日，共12景，根据研究区行政区边界裁剪拼接后的遥感影像如图1(a)和图1(c)所示。Sentinel-2号卫星(Sentinel-2A、Sentinel-2B)携带一枚多光谱成像仪(MSI)，可覆盖13个光谱波段(见表1)。具有光学遥感卫星常见的可见光、近红外和短波红外波段，还包含3个应用在高光谱数据中的红边波段，具有不同的空间分辨率。Sentinel-2号数据是在红边范围内的3个波段对监测植被的辐射强度非常有效[14-15]。

表1 Sentinel-2遥感影像波段信息

本研究中的地形数据来源于2017年研究区的二类调查数据，以及30 m分辨率的数字高程模型(DEM)数据。根据研究区行政边界拼接裁剪后的数字高程模型数据如图1(b)和图1(d)所示。

2.2 特征因子提取

由于遥感影像、数字高程模型数据、二类调查数据坐标系不同，在数据整合前需进行坐标系转换，本研究统一使用CGCS_2000_3_Degree_GK_CM_120E坐标系[16]。

光谱特征因子：Sentinel-2共包含13个波段，B1波段与B10波段与本研究无关，故选取剩余11个单波段因子。基于单波段因子进行波段运算，提取12个常见植被指数，即土壤调节植被指数、比值植被指数、非线性指数、归一化水指数、归一化植被指数、改进的归一化植被指数、归一化红外指数、归一化差异绿度指数、增强型植被指数、差值植被指数、类胡萝卜素反射指数、绿光叶绿素指数。根据Sentinel-2特有的红边波段提取7个植被指数：比值红边植被指数、归一化红边植被指数1、归一化红边植被指数2(、改进的归一化红边植被指数、红边叶绿素指数、非线性红边指数、反转红边叶绿素指数[17]。植被指数的具体计算公式见表2。

纹理特征因子：本研究采用灰度共生矩阵法提取遥感影像的纹理特征[19]。以B2、B3、B4和B8等4个10 m分辨率波段组合后进行主成分分析，使用ENVI软件提取第一主成分的8个纹理特征因子：均值、方差、协同性、对比度、相异性、熵、角二阶矩、相关性。

表2 植被指数

地面数据特征因子：从二调数据中获取每个小班的坡向、坡位、土层厚度、腐殖质厚度、植被覆盖度、群落、优势树种、年龄、郁闭度等9个自变量因子，其中坡向、腐殖质厚度、群落、优势树种为类别特征因子。从数字高程模型数据中提取海拔、坡度和坡向等3个自变量因子。

本研究共提取50个特征因子(包含四个类别特征因子)。使用ArcMap软件以森林资源细班为单元整合所有自变量因子，获得淳安县数据133 734条、临海市数据59 636条。删除字段缺失数据，再根据3倍标准差原则去除部分数据，共获得淳安县研究数据40 216条、临海市研究数据19 840条。

2.3 基本模型

K最近邻算法(K-NN算法)：K最近邻算法是比较成熟的非参数方法中最典型的方法之一。该算法通过数据所有的特征计算新数据与数据集中数据点的距离(该距离通常采用欧式距离以及曼哈顿距离)，按照距离递增次序进行排序，返回K个点的加权值作为预测值。该算法没有过多的理论推导与步骤，实现较为方便。

图1 研究区卫星影像图和数字高程模型数据

梯度提升迭代决策树(GBDT)：梯度提升迭代决策树是一种集成学习算法，由决策树和梯度提升组成，GBDT采用基函数的线性组合以及不断减小训练过程产生的残差(预测值与真实值的差值)来完成数据回归。其核心是计算上一轮迭代的负梯度，减少残差的梯度方向上建立新的决策树，通过不断迭代改进来得到更准确的预测结果。该算法很好地降低了多种因素之间的信息冗余程度。

极端梯度提升算法(XGBoost算法)：极端梯度提升算法是一种基于决策树的集成机器学习算法。XGBoost算法的基本思想是首先建立一个基分类/回归器，逐步加入新的分类/回归器，每加入一个分类/回归器后，再次计算其目标函数的值，保证在迭代过程中，目标函数值逐步下降，以此不断提升模型的表达效果。该算法有较强的泛化性能，能降低过拟合，适用性较强，在速度和效率上都优于其他梯度提升算法。

梯度增强集成分类器(CatBoost)：传统的梯度增强对样本梯度计算依赖于样本自身，噪声点会带来预测偏移，最终导致过拟合，为解决这一问题，CatBoost首先对整个数据集进行若干次排序，去除第i条数据，针对前i-1条数据，分别计算损失函数以及梯度，并建立残差树，最后将残差树累加到原模型上，有效避免了预测偏移，减少过拟合。而且CatBoost算法能自动处理类别特征，并根据特征的内在联系将原有类别特征进行组合，丰富特征维度，提升预测结果的精确性。

堆叠法(Stacking)模型融合：Stacking是一种集成学习方法，基于K折交叉验证的思想，通过元分类器或元回归器来融合多个分类或回归模型。首先基于完整的训练集对基础模型进行训练，然后将每个基础模型的预测值作为元模型的训练集，用新的训练集对元模型进行训练。Stacking通常会使用不同的机器学习算法来保证基础模型的多样性，并通过元模型以最优方法对初步回归结果进行整合，可以充分利用基础模型的预测结果，而且可以发现纠正基础模型中的部分偏差，使得准确率得到提高，同时降低过拟合的风险。本研究使用CatBoost和GBDT作为基础模型，XGBoost作为元模型。

图2 Stacking模型融合示意图

Lasso特征选择：最小绝对收缩和选择算子(Lasso)是一种基于L1范式的特征选择方法，属于嵌入式特征选择。Lasso通过加入L1惩罚项，将一些不重要特征的系数压缩为0，即剔除该特征，从而达到稀疏化和特征选择的目的。与其他变量选择方法比较，Lasso方法能够较准确地选择出重要特征，而且还具备特征选择的稳定性。

类别特征编码：本研究共选取四个类别特征。由于K-NN、GBDT和XGBoost不直接支持类别特征，因此需要对类别特征进行编码。根据类别特征属性，本研究采用两种类别编码方式。

标签编码：使用腐殖质厚度、群落为定序变量，采用标签编码，将类别特征依照特质以数字划分。将腐殖质厚度的{“薄”，“中”，“厚”}、群落的{“简单”，“较完整”，“完整”}转换为{0，1，2}。

独热编码：使用坡向取值为{“阴坡”，“阳坡”，“无坡向”}、优势树种取值为{“杉木”，“马尾松”，“阔叶混交林”等8种优势树种}为定类变量，采用独热编码，其方法是使用N位状态寄存器来对N个状态进行编码，每个状态都有独立的寄存器位，在任意时候，其中只有一位有效。坡向编码结果如图3所示。

图3 独热编码图

2.4 模型构建

因本文研究数据涉及到的自变量因子有46个(不包括4个类别特征)，为了找出对森林蓄积量具有更大的影响作用的因子，进一步提升模型的估测精度及其泛化能力。使用Lasso特征选择方法，从繁多的特征因子中选取部分重要特征，建立K-NN、GBDT、XGBoost和CatBoost等4种单一模型，并使用控制变量法探索类别特征对模型指标的影响，然后使用Stacking方法对单一模型进行融合，分析比较融合模型与单一模型的优劣。本文选用的Lasso特征选择以及机器学习模型均通过网格搜索(GridSearchCV)进行参数调优。

建模的具体步骤如下：(1)将所有样本按照7∶3的比例随机划分训练集与测试集；(2)通过网格搜索对4个单一模型进行反复调参，建立最优参数估测模型；(3)使用最优参数模型对测试集进行估测，得到单一模型预测结果；(4)检验类别特征对模型性能指标的影响，分析模型的普适性；(5)使用Stacking方法融合多个单一模型，并与单一模型估测结果进行对比。

2.5 模型精度评价

本研究选用十折交叉验证的方式来评估模型的精度。采用的性能指标有决定系数(R2)、均方误差(MSE)、平均绝对误差(MAE)、平均百分比误差(MAPE)，计算公式如下：

由于使用K-NN算法前需要对数据进行归一化，而树模型并不需要对数据进行归一化，从而导致数据量纲不统一。使用不同量纲下的MAE和MSE不能体现模型之间的差异，故使用均方误差与因变量平均值的比值(相对均方误差)替代性能指标的均方误差，平均绝对误差与因变量平均值的比值(相对平均绝对误差)替代性能指标的平均绝对误差。

3 结果与分析

3.1 数值特征选择

本文共提取数值特征46个，类别特征4个。为了减少信息冗余，降低共线性，使用Lasso方法对46个数值特征进行特征选择。经参数调优后，参数阿尔法取值0.003，系数阈值0.05，并选取系数值大于0.05的自变量因子作为建模特征因子(包括郁闭度、比值红边植被指数、改进的归一化植被指数、年龄、方差、相关性、相异性和熵等特征因子)。

3.2 类别特征对单一模型的影响

为比较类别特征(坡向、腐殖质厚度、群落和优势树种)对于模型性能指标的影响，使用淳安县数据，在未加入类别特征的情况下，经筛选后的8个数值特征进行建模和运用4个类别特征及筛选后的8个数值特征进行建模，模型性能指标见表3。

由表3可知，类别特征的加入对模型各项性能指标的提升有着显著影响。在加入类别指标后，GBDT模型的R2提高了0.12，CatBoost模型的R2提高了0.16；GBDT模型的相对均方误差从55.55%降低到32.19%，CatBoost模型的相对均方误差从54.0%降低到26.7%；GBDT模型的相对平均绝对误差从28.32%降低到20.17%，K-NN模型的相对平均绝对误差从34.02%降低到24.47%；K-NN模型的平均百分比误差从36.74%降低到33.75%，CatBoost模型的平均百分比误差从36.07%降低到22.57%。因此，加入类别特征后，4种模型的各项性能指标均有显著改善，CatBoost模型总体表现最好，各项指标优于其他3种模型。加入类别特征后的4种模型估测结果散点图如图4所示。

表3 淳安县森林蓄积量估测模型的性能指标

由图4可知，CatBoost、XGBoost、GBDT、K-NN等4种模型的R2逐渐下降，预测点与真实值之间差值越来越大，观察点在回归直线附近愈发稀疏，模型的可解释性降低。

3.3 单一模型的普适性

为探究模型跨空间普适能力，并进一步验证类别特征的加入对模型性能指标的影响是否具有普适性。将4种模型应用于临海市森林蓄积量估测，估测结果性能指标如表4所示。该结果进一步验证了类别特征对模型的优化作用具有一定的普适性，且XGBoost和CatBoost模型优于K-NN以及GBDT模型。

图4 4种模型对淳安县森林蓄积量预测的拟合

表4 临海市森林蓄积量估测模型的性能指标

虽然4种模型的R2相较淳安县的森林蓄积量估测结果有所下降，但是相对均方误差，相对平均绝对误差以及平均百分比误差等3项指标均好于淳安的估测结果，说明模型具有一定的跨空间普适性，在不同地区均能取得较好表现。加入类别特征后的蓄积量估测散点图如图5所示。

3.4 融合模型

为验证融合模型与单模型的优劣，基于前述筛选的类别特征及数值特征，以CatBoost和GBDT为基础模型，XGBoost为元模型，进行Stacking模型融合，与相同情况下单模型中表现最好的CatBoost进行对比(见表5)。由表5可知，模型融合后，除相对均方误差外，其余性能指标较CatBoost均有所改善，说明融合模型的综合效果优于单一模型。

图5 4种模型对临海市森林蓄积量预测的拟合

4 结论与讨论

遥感影像包含了更广域的、动态的地面覆盖物信息，为森林蓄积量估测提供了范围更广、实时性更强的数据来源，结合DEM、森林资源二类调查数据，为森林资源动态监测提供更为有效的支持。本研究使用2017年淳安县及临海市的Sentinel-2遥感影像、二类调查以及DEM数据，经过数据预处理以及特征选择，建立了单模型和多模型融合的蓄积量预测模型。通过控制变量法探索类别特征对于模型指标的影响，证明类别特征的加入能有效提高模型性能；实验中，使用临海市森林资源数，进一步据验证模型的间普适性以及类别特征对模型精度的影响；Stacking融合模型与单模型相比，综合指性能标更佳，估测精度更高。

在所选取的特征因子中，获取数值特征46个，Lasso特征选择方法有效减少了特征数量，使建模所用数值特征缩减到8个，加快了模型训练速度、提高了模型泛化能力，并使模型具有更好的解释性；类别特征的加入明显优化了模型，相比于常见的根据类别特征将数据划分为不同类别，建立模型方法更加简便易用；单模型中，CatBoost综合表现最佳，但融合模型Stacking的综合性能指标更优于CatBoost模型；Sentinel-2遥感影像特有的三个红边波段对于森林蓄积量的估测结果有较大的影响，与前人研究结果一致[19-20]。

但受遥感影像精度及与二类调查数据获取时间差异等资源的限制，本研究仍存在一些不足，需要进行优化改进：研究所选取的研究区域空间跨度较大，但是仍未超出省级界限，可选择空间跨度更大的地区作为研究区域，进一步分析模型的普适性；使用了Sentinel-2的单波段与基于部分波段提取的纹理特征因子，但是缺少对纹理特征因子的深入探究，可使用不同窗口大小、不同步长、不同波段组合来提取纹理特征，探究不同纹理特征对模型精度的影响[21]；本研究使用遥感影像分辨率为10 m，若能获取更高分辨率的卫星遥感影像和DEM数据进行研究，则可以进一步提高预测精度；除光学遥感外，雷达遥感影像也常被用来作为蓄积量估测的数据源，雷达遥感影像可穿透森林冠层，获取植被垂直结构参数，可应用于蓄积量预测[19]，进一步探讨雷达影像对蓄积量估测模型精度的影响。