基于Sentinel-2光谱与地形特征的山区森林分类——以武夷山国家公园为例

张春莹，江 洪，，，林敬兰，岳 辉

（1.福州大学 空间数据挖掘和信息共享教育部重点实验室，卫星空间信息技术综合应用国家地方联合工程研究中心，数字中国研究院（福建），福建 福州 350108；2.福建省水土保持试验站，福建 福州 350002；3.长汀县水土保持中心，福建 长汀 366300）

森林是全球陆地生态系统的重要组成部分［1］.传统的森林生态系统调查依靠实地调查来收集相关信息，但是需要大量的人力、物力和财力，并且容易受到区域的限制［1-2］.遥感技术在林业资源管理中具有监测范围广、成像速度快、重访周期短、数据成本低等优势［3-4］.传统基于像元的分类方法仅利用了地物的光谱信息，因此难以区分光谱信息相似的地物，而面向对象分类方法的基本处理单元是对象，可在此基础上充分利用不同对象的光谱、空间和纹理等相关特征［5］.光学图像中丰富的光谱变量是遥感图像分类中最重要的变量，研究表明，近红外和红边波段是植被较为敏感的波段［6］.例如，李丹等［7］通过采用WorldView-2 影像来识别典型乔木树种，研究发现红边波段明显改善了识别结果；Tigges 等［8］研究发现RapidEye 影像的红边波段能显著提升植被的分类结果的准确性.

欧洲空间局（European Space Agency，ESA）分别于2015 年6 月和2017 年3 月成功发射了Sentinel-2 系列卫星，该卫星为森林植被的识别提供了新的数据源，其搭载的多光谱仪（MSI）已被应用于树种分类实验中［9-10］.例如，于婉婉等［11］研究发现不同生长期的优势树种在Sentinel-2数据的红边波段表现有不同差异；黄双燕等［12］利用Sentinel-2影像提取干旱区农作物的分类信息，研究发现红边特征对于识别不同作物间的物候差异表现更为敏感.此外，在地形复杂的山区，山区地物通常形成地带性分布，辅助以地形因子进行地类信息提取可显著提升分类精度.例如，陈元鹏等［13］利用GF-1影像在山地工矿复垦区进行地类信息提取时加入地形因子得到了较高精度提取结果；Dubeau等［14］利用Landsat影像光学、雷达和地形变量的组合方式提取达布斯湿地信息，分类效果较好.考虑到植被生长茂密的区域多位于山地丘陵区，在进行地物提取时复杂的地形条件会对结果造成较大干扰，导致在开展森林植被信息提取时精度难以满足林业资源管理的需求.因此，采用Sentinel-2 影像的红边光谱特征辅助地形因子开展山区森林分类研究很有必要.鉴于此，笔者以武夷山国家公园为研究区，以较高分辨率的Sentinel-2 影像和ALOS PALSAR 数据为数据源，采用多尺度分割方法对多光谱影像进行分割以构建影像对象，提取27个光谱特征变量并计算其最优特征空间，将地形信息与最优特征空间相结合组成3种不同的特征变量子集，在此基础上对武夷山国家公园进行森林分类，以期为地形复杂的山区高精度森林分类提供有效的技术手段.

1 研究区概况与数据源

1.1 研究区概况武夷山国家公园（117 °24 ＇13 ″～117 °59 ＇19 ″ E，27 °31 ＇20 ″～27 °55 ＇49 ″ N）（如图1所示），总规划面积1 001.41 km²，属于中亚热带季风气候区，地貌以山地和丘陵为主，海拔高度差异悬殊，最高处海拔为2 160.8 m，最低处海拔为176.1 m，植被类型多样，森林覆盖率为93.74%.

图1 研究区示意图

1.2 数据源与预处理

1.2.1 数据源本文所用数据包括Sentinel-2 数据、ALOS PALSAR 高程数据、植被分布图，Sentinel-2 数据下载于欧洲航天局（ESA）官方网站（http：//scihub.copernicus.eu/），数据获取时间为2019 年9 月24 日，ALOS PALSAR高程数据下载于美国阿拉斯加卫星设备分布式活动档案中心（https：//search.asf.Alaska.edu/#/）.植被分布图来自国家林业和草原局昆明勘察设计院.所采用的Sentinel-2相关波段参数如表1所示.

表1 Sentinel-2波段介绍

1.2.2 样本数据通过对武夷山国家公园植被分布图的读取，参考国家林业和草原局出台的《国家森林资源连续清查技术规定 2014》，并结合研究需求，得到本文实验区所涵盖的地物类别，包括阔叶林、针叶林、竹林、农业用地、建设用地、水域.通过Google Earth 高分辨率影像（http：//earth.google.com/）选取样本点，遵循样本点的选取尽可能均匀分布的原则，样本点的空间分布如图2 所示，共选取750 个样本点，158个针叶林样本，150个阔叶林样本，183个竹林样本，62个农业用地样本，91个建设用地样本和61个水域样本，其中训练样本533个，验证样本217个.

图2 样本点空间分布

1.2.3 数据预处理笔者使用的是Sentinel-2 影像中除沿海大气气溶胶波段和卷积云波段以外的10 个波段，首先采用欧空局提供的Sen2Cor插件将下载的数据级别由L1C级转化为L2A级，即将10个波段的大气表观反射率处理成大气底层反射率，然后在SNAP软件中全部重采样为10 m，最终得到实验所需的包含10个波段的L2A级Sentinel-2影像数据.

2 研究方法

采用面向对象的方法进行分类，影像分割是面向对象分类方法中的重要步骤，分割是从子像素开始，从上到下合并区域，直到满足异质性标准［15］，特征提取是利用影像的光谱和空间特征等来提取分割对象所需要的类别信息.基于eCognition9.0遥感处理平台对预处理后的影像进行面向对象分类实验.为了检验Sentinel-2 光谱特征和地形因子对森林分类的能力和贡献，构建了3 个特征组合方案：方案1：参考光谱指数；方案2：参考光谱指数+红边指数；方案3：参考光谱指数+红边指数+地形因子.3 个特征组合方案作为随机森林分类器的输入特征，比较不同特征组合方案在分类精度提升上的表现，技术流程如图3所示.

图3 技术流程

2.1 影像分割与最优分割尺度评价影像分割是面向对象分类的一个重要部分，分割尺度是决定分割效果的关键.选用基于可见光和近红外波段的面向对象多尺度分割方法进行影像分割.为避免由于分割尺度不合适所出现的欠分割或过分割的弊端，采用ESP2（Estimation Scale Parameter）工具辅助目视解译来确定影像的最优分割尺度.ESP2 来源于Dragut 等［16］提出的ESP 算法，该工具计算影像对象异质性局部方差的变化率，生成局部方差变化率曲线，当曲线出现峰值时，表示该点变化率最大，则该点相对应的分割尺度为最优分割尺度.

影像多尺度分割在可见光和近红外4 个波段上进行，经过多次试验，最终形状因子确定为0.2，紧致度因子确定为0.5，在此基础上，使用ESP2工具计算最优分割尺度，将分割尺度的步长设为1，起始分割尺度设为30，得到ESP 分割尺度评价结果（如图4 所示），变化率曲线有多个峰值，表示影像的最优分割尺度参数有多个，选择分割尺度40、50、89、111进行分割尺度效果的目视解译对比.

图4 ESP分割尺度评价

分割尺度为40、50、89、111 下的分割效果如图5 所示，图5a 和e 表示分割尺度为40 时，阔叶林、针叶林、农业用地等分布区域面积大的地物类型被分割得较为细碎，产生了过分割现象，不满足实际地物的分割需求；如果分割尺度过大，小面积地物会与大面积地物混杂在一起，导致分割地物的混淆，如图5c～d 和g～h，当分割尺度为89和111时，同一对象内包含了不同的分割地物，6种地物类型都表现出了不同程度的欠分割现象.而在最佳分割尺度下，既合并了同类地物又与其他地物明显区分开来，地物分割都达到了较好的分割效果，如图5b 和f，当分割尺度为50 时，各地物类型分割效果较好，分割对象符合地物实际空间分布边界，可在此基础上进行后续实验.因此，将最优分割尺度选定为50.

图5 不同分割尺度下影像分割效果

2.2 Sentinel-2光谱特征提取根据参考的相关文献［9］，［17-26］，构建27个光谱特征，包含Sentinel-2影像的10个原始波段和17个光谱指数（如表2所示），其中包括11个红边指数.

表2 基于Sentinel-2影像的光谱指数

2.3 地形因子提取地形因子中包含丰富的地形表面形态特征信息，将各种地形因子结合在一起，能够较好地刻画地表起伏的变化.因此，根据本文需求，将高程、坡度和坡向作为地形因子变量，ALOS PALSAR 数据即DEM（Digital Elevation Model）数字高程模型含有该区域海拔高度信息；坡度、坡向特征可使用ArcMap软件从DEM中计算获得.

2.4 特征优化与分类随机森林（Random Forest，RF）算法是Leo 等［27］于2001 年提出的由多个决策树分类器集成的一种机器学习算法.采用RF 模型用来作为分类与计算和选择特征重要性的工具.该算法采用自助法重采样技术，将原始样本中有放回地随机抽取2/3的数据作为训练集，训练集中的每个样本都由多个决策树分类，产生多个分类结果，投票决定每个样本的分类类别.未被抽取的样本称为袋外数据（Out of Bag，OOB），OOB误差越低，RF模型准确性越高.

RF 算法在完成特定分类任务的同时，还可评估模型中每个特征变量的重要性.由于Sentinel-2 影像的原始光谱波段丰富，在后续研究中众多学者研究出了相关光谱指数，面向对象的分类方法可输入众多特征参与分类，但较多特征会导致信息的冗余，降低分类的速度和精度.因此，为充分利用影像对象丰富的光谱与空间特征，采用RF算法中的基尼系数（Gini）评估各特征的重要性，确定最优特征空间.Gini系数数值越大，则表明该特征的重要性也越高.最终根据重要性排序以及OOB 误差进行特征优选，确定最终的特征变量子集.

2.5 精度评价本文中3 种特征组合方案的森林分类精度采用混淆矩阵进行评价，评价指标有总体精度（Overall Accuracy，OA）、Kappa 系数、生产者精度（Producer’s Accuracy，PA）、用户精度（User’s Accuracy，UA）以及调和平均值（F1）.其中，F1的计算公式［28］为

3 结果与分析

3.1 特征重要性排序及评价采用RF Gini 系数计算得到Sentinel-2 影像27 个光谱特征的重要性排序（如图6 所示）.其中重要性排序前五的特征是LSWI、B2（蓝波段）、B4（红波段）、B12（短波红外2 波段）和Ndre1.重要性排名最高的特征为LSWI，为11.79%，该特征由近红外波段B8 与短波红外1 波段B11 计算得到.其次为蓝波段B2和红波段B4，重要性分别为10.35%和7.50%，其中，蓝波段B2的重要性得分远高于其他波段，可能受到其他相关特征的影响［29］.B12、Ndre1和GNDVI的重要性较高，分别为4.48%、4.37%和4.30%，由于研究区植被覆盖度高，植被在此处具有较强的吸收性，所以红光波段达到贡献度高.贡献度最低的为NDVIre3，重要性得分仅有0.81%.

图6 特征变量重要性排序

通过改变特征数量得到27 个光谱特征数量与分类精度的关系如图7a 所示，前期精度值随特征数的增多快速增加，此后持续保持上下浮动的趋势，当特征数为17 时，精度达到最大值0.911 0.因此，在森林分类中选择重要性排序较高的前17个光谱特征构建最优特征空间.对这17个特征中所包含的Sentinel-2原始波段数量进行统计（如图7（b）所示），数量最多的是B4（红波段），8个；其次是B8（近红外波段），5个；B5（红边1 波段），5 个.此外，B6（红边2 波段）和B7（红边3 波段）也较多，进一步表明红边波段在森林分类中具有较高的贡献度.

图7 特征数量与分类精度的关系和重要特征波段统计

3.2 分类结果与精度评价采用面向对象多尺度分割方法结合不同特征组合方案进行森林分类，如图8所示.在方案1中仅利用参考光谱特征分类（图8b）时，“椒盐现象”十分明显；在方案2中加入红边特征参与分类（图8c）后效果较好，各地类的图斑破碎度有所降低；在方案3 中将光谱特征与地形因子结合进行分类（图8d）后效果最好，各地类“椒盐现象”得到改善.此外，将分类结果与Sentinel-2影像对比可以看出，阔叶林、针叶林与竹林交错分布，阔叶林分布面积最大，最连续；竹林与阔叶林镶嵌分布，且与农业用地交错分布；针叶林分布较为集中.

图8 分类结果

不同特征组合方案的分类结果精度对比如图9 所示.方案1、2 和3 的总体精度（OA）与Kappa 系数分别为88.13%、0.856 4，89.50%、0.871 0 和90.87%、0.887 8.对比分析3 种特征组合的森林分类精度，仅使用参考光谱特征分类时，阔叶林、竹林和农业用地的生产者精度（PA）较低，三者的F1分别为0.81、0.88和0.67.加入红边特征进行分类后，阔叶林和竹林的生产者精度（PA）分别提升至0.88和0.90，其F1都提升了0.03，阔叶林和农业用地的用户精度（UA）分别从方案1 的0.76 和0.61 提升至0.79 和0.70.在光谱特征基础上加入地形特征参与森林分类时，竹林和农业用地的生产者精度（PA）分别提升至0.92和0.77，相对于方案1来说分别提升了0.08和0.05，相对于方案2来说提升了0.02和0.15，F1分别提高至0.92和0.79，比只使用参考光谱指数提升了0.03 和0.12，阔叶林和农业用地的用户精度（UA）提升至0.82 和0.81，相对于方案1 提升了0.06 和0.20.由此可见，方案3 的总体分类精度（OA）和Kappa 系数最高，表明将Sentinel-2的光谱特征与研究区地形因子相结合作为分类器的输入变量参与分类时，能在一定程度上提高森林分类精度，得到最佳的分类效果.

图9 不同特征组合的分类精度对比

不同特征组合分类混淆矩阵如图10所示.使用方案1进行分类的效果较差，如图10a所示，特别是农业用地与竹林和阔叶林之间存在一定程度上的地物混淆，导致分类结果整体来看效果欠佳，这是因为竹林、阔叶林和农业用地在Sentinel-2 影像上视觉差较小，因此仅利用参考光谱指数很难将三者区分开.如图10b所示，在方案2中加入红边特征后，竹林精度有了明显提高，说明竹林对于红边数据敏感度更高.方案3中加入地形因子后，如图10c所示，阔叶林与农业用地2类地物相比方案1和方案2来说可更为有效的区分开，说明地物的分布与地形有一定程度上的相关性.

图10 不同特征组合的分类混淆矩阵

从3 种分类特征组合结果可见，方案3 使用Sentinel-2 的光谱特征辅助地形信息的总体分类精度（OA）、Kappa 系数和F1最高，其次为方案2，方案1 的分类精度最低，效果最不好.阔叶林、针叶林、竹林和农业用地在方案3中取得较好的效果.由此可见，由于植被对红边波段的敏感度更高，红边波段对于森林分类也显得更为重要，此外，加入地形因子可提高地形复杂区域地类信息提取的精度，这与陈元鹏等［13］、Dubeau等［14］的研究结论相同.

3.3 分类结果与地形因子分异统计分析为分析各地类在高程、坡向和坡度上的分布特点，通过自然间断法对高程、坡向、坡度3种地形因子重分类［30］，利用ArcMap分区统计工具统计不同地形因子等级内森林类型面积所占比例，不同森林类型在高程、坡向、坡度上的统计结果如图11所示.各森林类型面积比例随着高程的不同而出现明显差异如图11a 所示，其中，针叶林的面积比例随着高程的升高波动起伏较大，在1 427～1 661 m 和1 662～2 162 m 范围内占比分别为63%和57%；农业用地、水域和建设用地在178～446 m范围内占比为46%.各森林类型面积比例在不同的坡向上也有明显的地形响应如图11b所示，其中，农业用地、水域和建设用地在平面上面积占比为79%，阔叶林在东坡、东南坡和南坡的分布面积比例远高于其他森林类型.各个森林类型的面积比例随着坡度的升高而出现明显且较为规律的变化如图11c 所示，阔叶林、竹林和针叶林的面积比例随着坡度的增大逐渐升高，农业用地、建设用地和水域随着坡向的增加面积占比逐渐减小.由此可见，不同森林类型的分布与地形因子之间存在一定的规律性，所以在进行山区森林分类时加入地形因子是提高分类精度的一种有效手段.

图11 分类结果与地形因子统计

4 结束语

利用多尺度分割的方法对Sentinel-2 影像进行分割以构建影像对象，在此基础上提取影像的光谱特征，将其与研究区地形因子相结合，构建3种特征组合方案对武夷山国家公园进行森林分类.结果表明：

1）为避免特征数量过多导致的信息冗余，采用RF Gini系数对特征空间进行重要性评估，在Sentinel-2影像上重要性排名前17 的最优特征空间中，红波段、近红外波段和新增的红边波段在森林分类中的重要性较高；

2）不同特征组合方案参与分类时，光谱特征与地形因子组合变量的分类方案精度最高，总体分类精度为90.87%，Kappa系数为0.887 8.因此，采用Sentinel-2红边光谱特征结合地形因子的多特征组合方法是提高山区森林分类精度的一种有效手段.

随着遥感影像空间、光谱及时间分辨率的提高，充分利用Sentinel-2 影像红边波段光谱信息并将其与地形信息相结合的方式是实现森林高精度分类的一种有效手段.但笔者仅讨论了单一时相下将Sentinel-2影像的光谱特征和地形特征相结合的森林分类方法，下一步研究如何融合多时相遥感影像，将其应用于山区复杂地形条件下的森林制图.