基于Sentinel-2A影像的森林类型提取研究

郑振灿，陈文惠，2*，林莉平，刘育圳

（1.福建师范大学 地理科学学院，福建 福州 350007；2.湿润亚热带山地生态国家重点实验室培育基地，福建 福州 350007）

森林被誉为“地球之肺”能够净化空气，同时也是天然的“绿色水库”和“天然氧吧”，具有涵养水源、保持水土、调节气候等生态功能，不仅是野生动植物的良好栖息地，还为人类提供大量的木材、药材、松脂等林产品[1-2]。我国森林资源虽然丰富，但由于我国是人口大国，森林资源仍处于资源总量相对不足、质量不高、分布不均的状况[3]。因此，及时、高效、准确地获取森林资源现状信息对促进森林相关产业的可持续性发展以及更加科学地管理森林资源具有重要意义。

我国获取森林资源信息的传统方式主要以地面调查为主。地面调查不仅需要消耗大量的人力、物力和财力，而且调查持续时间比较长，在这个过程中林地面积可能产生变化[4]；同时，有些样地在深山老林、高寒地带或者悬崖峭壁上，受地形地貌的影响，样地可达性大打折扣，无法进行实地调查，会使最终的估测结果产生误差。

遥感技术以其探测范围广、获取信息速度快、周期短、手段多、且受地面条件限制少等技术优势，已被广泛应用于林业、农业、水文、测绘、气象等众多领域[5]。中低分辨率遥感影像（如Landsat TMETM+OLI、MODIS）由于受空间分辨率和光谱分辨率的限制，整体分类精度和估测精度往往不尽如人意。高空间分辨率遥感影像能带来强烈的视觉体验，同时空间分辨率的提高也使不同森林类型异质性增强，表观类型增多，引起“同物异谱”和“同谱异物”现象，给遥感影像分类带来挑战。Sentinel-2A 遥感影像具有较高的空间分辨率（最高为10 m）和较丰富的光谱分辨率（0.443～2.190 μm），使得地物的纹理特征和结构信息更加丰富，具有提高森林类型分类精度的潜力；数据处理简单且免费使用，在土地覆盖一类调查中已被广泛应用，而用于森林类型精细分类方面的研究还比较少。

因此，本研究充分利用Sentinel-2A遥感影像丰富的光谱信息以及较高空间分辨率的产品优势，提取影像的光谱特征、指数特征和纹理特征，结合森林二类调查数据、野外调查数据和DEM数据等辅助数据，在对所选特征进行森林类型可分离性分析的基础上，确定最佳的特征参与森林类型分类，利用支持向量机、决策树、随机森林以及组合分类器的分类方法，依次逐层提取研究区的地表覆盖以及森林树种信息，以期利用遥感影像的多元特征以及组合分类器的优势提高研究区森林类型提取精度，为Sentinel-2A遥感影像在森林类型信息提取方面的应用提供参考。

1 研究区与数据

1.1 研究区

本研究以福建省三明市尤溪县为研究区，地理坐标处于25°50'N—26°26'N、117°48'E—118°40'E 之间，如图1所示。尤溪县东邻闽清县和永泰县，南接德化县，西连大田县和沙县，北毗邻南平市。全县总面积为3 463 km2，总人口为45余万人，地貌以低山、丘陵和盆地为主，地势东西高、中部低，平均海拔1 472 m，属于亚热带季风性湿润气候，温暖湿润，四季分明，年平均气温19.2 ℃，年平均降水量1 650 mm，土壤以沙黄壤和沙红壤为主。

图1 研究区地理位置及其标准假彩色影像（黄色点为野外样本点）Figure 1 The geographic location of the study area and its standard false colour image(the yellow dots are the sample points in the field)

尤溪县森林资源丰富，林业用地面积达2 840.17 km2，占全县国土面积的83.0%，森林覆盖率达到76.94%。森林类型主要以针叶林和针阔混交林为主。物种资源丰富，林木生长旺盛，林相整齐[6]。

1.2 数据来源与处理

1.2.1 遥感影像数据

Sentinel-2A卫星携带有一台光谱成像仪（MSI），卫星高度为786 km，覆盖有13个波段，分别为空间分辨率为10 m的3个可见光波段和1个NIR波段、空间分辨率为20 m的4个Red-edge波段和2个SWIR波段以及空间分辨率为60 m的海岸、水汽和卷积云波段，主要用于陆地植被、土壤和水资源等方面的全球陆地观测并且可用于紧急救援服务等。本研究选取天气晴朗、无云、覆盖研究区范围的4景Sentinel-2A影像数据，影像获取时间均为2018年4月19日。由于该数据为已经过正射校正和几何精校正的大气表观反射率数据，因此采用Sen2cor软件对其进行大气校正，转化为地表反射率数据，再将其空间分辨率为20 m的波段重采样成空间分辨率为10 m，与可见光和NIR波段进行波段叠加，得到研究区10 m空间分辨率的10个波段的Sentinel-2A数据。

1.2.2 辅助数据

辅助数据主要有森林二类调查数据、DEM和Google Earth 0.5 m影像。森林二类调查数据来源于林业部门调查数据，野外样本点大多沿道路分布。为使样本点分布较为均匀，选取部分森林二类调查点作为训练样本和精度验证样本。Google Earth 0.5 m影像主要用于对野外调查和森林二类调查数据点的筛选提供高空间分辨率影像。DEM主要用于地形因子的提取（如坡度、坡向等），其数据来源于中国科学院计算机网络信息中心地理空间数据云平台（http://www.gscloud.cn）。

2 影像分割与分类

2.1 影像分割

目前，面向对象分割算法大致可以分为三类：基于阈值分割算法、基于边缘检测分割算法和基于区域生长分割算法[7]。本研究采用eCognition 软件提供的多尺度分割算法，该方法是一种自上而下的区域合并算法，属于区域生长分割算法范畴。采取实验法探索不同地物类型的最优分割参数。

首先，通过控制变量法的原理，随意选一个分割尺度为30，保持形状因子不变，分析不同紧致度因子对不同地物样本的分割效果，寻找最优的紧致度因子。同理，依次按照该方法寻找最优的形状因子，最后保持形状因子和紧致度因子不变，设置起始分割尺度为10，间隔为5，从下往上依次分割至分割尺度为100。不同分割尺度的分割结果如图2所示，采用视觉检验法对分割效果进行评价。

从图2中可以看出，当分割尺度为100时，大面积的林地虽被分割开来，但林地与其它地物类型不同程度地混合在一起，处于欠分割的状态，对象内部的同质性和不同对象间的异质性均不能达到分类的要求；当分割尺度为50时，道路两边和城市区内的林地还存在欠分割现象，不能达到分类提取的需求；当分割尺度为30时，大面积的林地、道路两边的林地以及城区内部和水域附近的林地分割边界比较清晰，对象内部同质性较好，为最佳的分割尺度；当分割尺度为10时，各个地物类型的边界虽然很清晰，但是分割对象过于破碎化，处于过分割的状态，导致相同地物类型对应的特征值不均一，不利于地物的提取。通过不同地物的对比分析，最终确定最佳分割参数形状因子为0.2（光谱因子为0.8）、光滑度因子为0.5（紧致度因子为0.2），耕地、建设用地、水域等一级类分割尺度为30，杉木、马尾松、竹林和其他林地分割尺度为35。

图2 形状因子为0.2，紧致度因子为0.5下不同尺度分割结果Figure 2 The segmentation results of different scales with shape index as 0.2 and compactness as 0.5

其次，探索分割特征组合对分割效果的影响。在面向对象分割方法中，不同分割波段组合具有不同的分割效果，主成分变换尽可能集中影像所有波段的优势，在减少数据冗余的同时增强影像分割所需的信息[8]。为研究主成分变换对影像多尺度分割效果的影响，本研究在最佳分割尺度参数的基础上，对原始光谱波段加入主成分变换的前三个分量参与分割，分割结果如图3所示。

图3 不同方案分割结果对比Figure 3 Comparison of segmentation results of different schemes

从图3中可以看出，相比于原始光谱波段，加入主成分前三个分量对于城市区域的小面积林地、林地边界、裸土边界以及小区域林地边界有更好的分割结果，说明主成分变换的前三个分量的加入对于林地与其它地物类型有更好的分割结果，有利于本研究中林地的提取。最终得到最佳分割参数组合如表1所示。

表1 最佳分割参数组合Table 1 The best segmentation parameters

2.2 特征分析

2.2.1 光谱特征光谱特征是多光谱遥感影像进行树种识别的最基本的特征之一。基于Sentinel-2A影像上典型地类和主要森林类型的光谱特征，分析各地物类别在影像光谱上的差异和光谱区域的可分性。图4（a）、（b）分别为Sentinel-2A影像上典型地类光谱曲线和主要森林类型光谱曲线。

由图4（a）可见，Sentinel-2A影像上耕地和林地在Green波段（Band 2）和Red波段（Band 3）表现出明显的植被光谱响应“峰谷”特征。然而由于时相的原因耕地和林地在各个波段上的光谱特征表现出一致性，特别是在前4个波段的光谱特征值比较接近，在后6个波段的值域范围相差比较大，仅仅依赖影像光谱特征难以区分，需要加入其它的特征。而水体和建设用地与耕地和林地的光谱曲线的变化趋势的差异比较明显，水体在前3个波段与林地值域范围比较接近，其它波段范围的值域范围相差比较大，建设用地在各个波段内值域范围相差比较大。由图4（b）可见，Sentinel-2A影像上马尾松、杉木、竹林和其它林地表现出明显的植被光谱响应“峰谷”特征。杉木林和马尾松林都为针叶林，在Blue波段（Band 3）、Red波段（Band 3）、Red-edge 1波段（Band 4）以及SWIR 1波段（Band 9）和SWIR 2波段（Band 10）光谱响应几乎完全重叠，在Red-edge 2波段（Band 5）、Red-edge 3波段（Band 6）、Red-edge 4波段（Band 8）以及NIR波段（Band 7）值域范围有所差异，但并不能满足森林类型精细化提取的要求，若能加入其它显著特征（如植被指数、纹理及其它数据等）则有可能使得两类完全分开；竹林在SWIR两个波段（Band 9、Band 10）上具有明显特征，可能是竹林本身的生物学特征使其随着季节的变化而发生显著的变化；其它林地在Red-edge 2波段（Band 5）、Red-edge 3波段（Band 6）、Red-edge 4波段（Band 8）以及NIR波段（Band 7）上值域范围与杉木林、马尾松林和竹林上有所差异。

图4 不同地物类型光谱曲线Figure 4 Spectral curves of different features

图5 反映以Sentinel-2A 影像外业调查获取的马尾松、杉木、竹林以及其它林地样点在各个波段上的像素值分布。由图5可知，不管在Sentinel-2A影像哪个波段上，杉木和马尾松的特征值的混淆程度都比较高，单单依靠影像光谱特征很难将其完全区分开来；从图5（b）、（d）、（e）、（f）、（g）、（h）中可以看出，在Green（Band 2）、Red-edge 1（Band 4）、Red-edge 2（Band 5）、Red-edge 3（Band 6）、NIR（Band 7）、Red-edge 4（Band 8）波段其它林地与竹林的特征值混淆比较严重，但与杉木、马尾松的特征值具有一定的差异，虽存在不同程度的混淆，加入其它特征可使其完全区分开来；从图5（i）和（j）中可以看出，在SWIR 1（Band 9）、SWIR 2（Band 10）两个波段上竹林与马尾松林、杉木林和其它林地特征值有显著的差异；相比SWIR 1波段，SWIR 2波段中竹林更不易与其它林地的特征值混淆，更有利于竹林的提取。

此外，除了影像波段均值（Mean）外，还选取标准差（Standard Deviation）、亮度（Brightness）和贡献率（Ratio）作为Sentinel-2A影像主要森林类型的光谱特征参与森林类型分类。

2.2.2 指数特征

结合植被指数和水体指数进行分类可在一定程度上改善影像分类精度[9]。本研究通过对常用的NDVI、RVI、DVI、MNDWI、LSWI、ARVI、SAVI、NDI45、MCARI、S2REP、IRECI、PSSRa、NDVIre1、NDVIre2、NDVIre2等15个指数在Sentinel-2A影像外业调查获取的马尾松、杉木、竹林以及其它林地样点在各个指数上的像素值分布进行研究，经过反复测试，最终选取NDVI、RVI、DVI、LSWI、MNDWI 等5 个常用指数以及MCARI、NDI45、PSSRa等3个红边波段光谱指数作为有效特征参数参与森林类型分类。

图6反映了MNDWI和LSWI指数在杉木、马尾松、竹林和其它林地等森林类型在各个指数上的像素值分布。从图6（a）、（b）中可以看出，指数特征的加入使得杉木、马尾松、竹林和其它林地的特征值峰值出现差异，特别是LSWI指数使杉木、竹林和马尾松的峰值差异比较明显。从图5（g）、（i）中可以看出，虽然杉木、马尾松在这两个波段各自的像元值分布并没有明显的区分性，但是在NIR波段中杉木像元值相比于马尾松像元值大，而在SWIR 波段中正好相反，经过LSWI 指数计算后反而拉大了二者之间的像元值差距，因而使杉木、马尾松的峰值差异比较明显。但是单单依靠指数光谱特征很难将杉木、马尾松、竹林和其它林地完全区分开来，还需加入其它有用的特征进行区分。相比影像波段光谱特征，波段之间的某种线性或者非线性的组合使得其杉木和马尾松之间更具有区分性。

图5 Sentinel-2A影像森林类型光谱值分布Figure 5 Spectral value distribution of forest types of Sentinel-2A image

图6 森林类型指数特征值分布Figure 6 Index feature value distribution of forest types

2.2.3 纹理特征

纹理特征描述了图像或者图像区域所对应目标的表面性质，对于具有规则变化结果的地物的提取是一个非常重要的特征（如茶园、果园、人工培育的杉木林等）[10]。纹理特征可用灰度共生矩阵进行表达。基于Sentinel-2A影像丰富的波段信息，通过对比不同纹理特征在影像样本点纹理特征上的差异，经过反复测试，选取均值（Mean）、相异度（DIS）、同质性（HOM）、熵（ENT）、相关性（COR）、对比度（CON）等6个显著的纹理特征参与影像分类，选取的窗口大小为9×9，选取的方向为135°。

图7反映了利用Sentinel-2A影像NIR（Band 7）波段和SWIR 1（Band 9）波段分别提取的均值（Mean）和相异度（DIS）纹理特征在各地类样本点的像素分布。由图7可知，均值（Mean）纹理特征参量的引入有利于建设用地、水体、耕地和林地之间具有更好的区分性，而相异度（DIS）纹理特征的引入能够使得建设用地与其它土地覆盖类型具有明显的差异。

图7 地类纹理特征值分布Figure 7 Texture feature value distribution of land

利用Sentinel-2A影像Red-edge 1（Band 4）波段和SWIR 2（Band 10）波段提取的均值（Mean）以及NIR波段提取的同质性（HOM）、相异度（DIS）纹理特征在各森林类型样本点的像素分布如图8所示。图8（a）表明Sentinel-2A 影像Red-edge 1（Band 4）波段所提取的均值（Mean）纹理特征的引入使得竹林、其它林地与杉木、马尾松之间有较为明显的可区分性，而杉木和马尾松表现出相同的纹理特征值，难以区分。图8（b）表明，Sentinel-2A影像SWIR 2（Band 10）波段所提取的均值（Mean）纹理特征与光谱特征一样，也在提取竹林中表现出明显的可区分性，验证了短波红外波段是提取竹林的有利波段。图8（c）、（d）表明，Sentinel-2A影像NIR（Band 7）波段提取的同质性（HOM）、相异度（DIS）纹理特征在其它林地与杉木、马尾松和竹林之间有明显的可区分性。

图8 森林类型纹理特征值分布Figure 8 Texture feature value distribution of forest types

2.2.4 地形因子

地形因子特征主要包括坡度和高程。在山地地区，由于坡度、坡向和高程的原因造成影像传感器太阳入射角不同，使得传感器所获取的太阳辐射能量有很大的差别[11]。研究区高程差异比较大，优势树种的生物光学特性、群落分布与地形因子的相关性比较强。同时，通过查阅相关文献资料和实地野外调查发现，研究区的耕地大多数分布于道路两边以及坡度比较小的区域，坡度是提取耕地中比较有用的特征之一[12]。因此，本研究利用获取的DEM数据计算坡度数据，将计算得到的坡度数据作为研究区地物类型提取的有效辅助特征，以期进一步提高森林类型分类精度。

2.3 分类策略

通过查阅相关研究区森林类型文献资料，参考《国家森林资源连续清查技术规定》[13]中土地类型划分标准、植被类型划分标准和《森林资源规划设计调查主要技术规定》[14]中林地分类系统，结合研究区森林植被特点及Sentinel-2A影像森林类型区分能力，建立研究区森林类型分类系统，如表2所示。

表2 研究区森林类型分类系统Table 2 Classification system of forest types of study area

鉴于研究区植被覆盖状况和Sentinel-2A影像特点，在充分分析Sentinel-2A影像所提取特征的森林类型可分离性以及特征组合对影像分割与分类的效果基础上，本研究总体上采取层次化信息提取的策略，由易到难逐层提取各类型信息，使层与层之间的误差以及不确定性降到最低，确保每层信息类型的分类精度，整体层次化分类技术流程图如图9所示。

图9 分类流程图Figure 9 Flow chart of classification

首先，对建设用地、水体、耕地等3类进行信息提取。对Sentinel-2A影像采用建设用地、水体和耕地3类地物的最佳分割参数进行影像分割，利用归一化植被指数NDVI＜0.56提取非林地的建设用地和水体，结合研究区实际状况和野外调查结果，通过设定坡度Slope＜10°提取耕地，避免耕地与森林类型发生混淆，将其与建设用地和水体合并作为Layer 1。在分析地物类型光谱特征可分离性的基础上，选取的窗口大小为9×9，选取的方向为135°，引入对建设用地具有很好区分性的Mean、DIS、HOM、ENT、COR、CON 纹理特征。针对1ayer1影像上三种地物类别采用在小样本中具有独特优势的支持向量机（SVM）分类方法[15]。

其次，对杉木、马尾松、竹林以及其它林地4类森林类型进行信息提取。对Sentinel-2A影像上非Layer 1区域进行合并，将其作为Layer 2，采用杉木、马尾松、竹林和其它林地的最佳分割参数进行影像分割。引入Mean、DIS、HOM、ENT、COR、CON 等纹理特征，NDVI、RVI、DVI、LSWI、MNDWI、MCARI、NDI45、PSSRa 等指数特征以及坡度地形因子特征，采用决策树和随机森林分类器分别进行森林类型分类。组合分类器能够实现各个分类器优势互补，具有取得比单个分类器算法更好的、精度更高的分类结果的可能[16，17]。采用改进的基于先验知识的投票法组合分类，即在进行分类之前，利用精度评价的混淆矩阵中的用户精度和制图精度的平均值作为投票的加权值，使用多数投票法对不同分类器的输出结果进行组合决策并输出分类结果。采用计算混淆矩阵的方法计算分类结果的生产者精度（Producer's Accuracy）、用户精度（User's Accuracy）、总体精度（Overall Accuracy）以及Kappa系数对分类结果进行精度评价[18]。

2.4 结果分析

利用上述面向对象的分层次森林类型分类方法得到研究区的森林类型分布。利用上述森林类型验证点对分别采用决策树、随机森林和组合分类器得到的森林类型分类结果进行精度验证，得到研究区森林类型精度的混淆矩阵，如表3 所示。从表3 中可以看出，组合分类器具有最高的分类精度，总体分类精度达79.80%，Kappa 系数达0.725；其次是随机森林，总体分类精度达77.67%，Kappa 系数达0.696；决策树分类器分类精度最低，总体分类精度为76.48%，Kappa系数为0.678。对于决策树分类模型来说，造成总体分类精度较低的原因主要是杉木和其它林地的错分，由于决策树模型较为简单，而杉木和马尾松都属于针叶林树种，其像元值和特征都比较接近，所以难以区分；而其它林地与杉木类型的错分也比较多，其它林地类型种类多且复杂，与杉木林影像特征较为接近，容易错分。虽然同样的错分现象也出现在随机森林和组合分类模型中，但是相比于决策树分类模型，随机森林模型比较复杂，在马尾松和竹林的分类精度有所提高，但是在杉木和其它林地的分类精度有所降低；组合分类器总体分类精度最高，其中杉木、马尾松、竹林的生产精度分别达到79.60%、73.86%、96.00%，用户精度分别达到80.67%、76.47%、88.07%，充分发挥了决策树和随机森林二者的优势，取长补短，总体分类精度分别提高了3.32%和2.13%。

表3 森林类型分类结果精度验证Table 3 Accuracy validation of forest types classification results

3 结论

本研究以Sentinel-2A影像为主要数据源，结合森林二类调查数据、野外调查数据和DEM 等辅助数据，探索基于面向对象的Sentinel-2A影像森林类型层次化提取方法。结果表明，运用多尺度分割算法建立分层多尺度的不同地类分割参数，逐层次提取地物类型有利于避免非林地地物类型对于森林类型信息提取的干扰，提高分类精度。主成分变换前三个分量的加入有利于提高地类的分割效果和森林类型的分类精度。通过建立分层多尺度提取规则，采用决策树和随机森林算法对研究区主要森林类型进行提取，总体分类精度分别为76.48%和77.67%，Kappa系数分别为0.678和0.696。基于分类精度混淆矩阵中的生产精度和用户精度平均值进行决策树和随机森林两种分类器的决策级投票融合，组合分类器的总体分类精度达79.80%，Kappa系数为0.725。相比于单一的决策树或随机森林分类器，组合分类器总体分类精度分别提高了3.32%和2.13%。

本研究通过对Sentinel-2A影像多元特征进行综合分析，探索符合Sentinel-2A影像特点和内在规律的森林类型信息提取方法，取得了较好的研究成果，对于Sentinel-2A影像在森林资源调查中的应用具有一定参考价值。但本研究也存在一些不足，以下两方面还有待进一步研究。

（1）分类精度的提高。本研究针对研究区杉木、马尾松和竹林的特点，采用分层提取的策略，其中竹林的生产精度和用户精度比较高，能满足林业应用的需求，但是杉木和马尾松存在一定的混淆，分类精度比较低，不能满足森林类型精细化提取的要求。因此，在今后的研究中可通过引进杉木和马尾松具有较高区分度的数据或者对不同数据源进行融合对其进行森林类型信息提取，同时对该提取方法是否具有迁移性进行研究。

（2）森林类型信息提取的进一步精细化。本研究旨在对研究区的杉木、马尾松和竹林进行信息提取，对于其它林地并未能进行更加精细化的划分，还达不到当今精细化的林业应用需求，在今后的研究中可适当对阔叶林树种或者经济林树种进行划分，做到更加精细化的分类。