基于Sentinel-2时序多特征的植被分类

郭文婷， 张晓丽

（1.北京林业大学 省部共建森林培育与保护教育部重点实验室，北京100083；2.北京林业大学 精准林业北京市重点实验室，北京100083）

森林资源调查可以统计不同林地的生物学特性、生态特征和经营价值，为森林的合理经营规划提供有效的数据支持［1-2］。由于中国各地区之间自然条件差异较大，森林类型多种多样，林分结构较为复杂，森林资源调查工作难度很大［3］。植被分类是研究森林资源状况和动态变化规律的基础，利用遥感手段可以更加快速、准确地识别植被类型［4］。因此传统的依靠人力的地面调查已经逐渐转变为基于遥感手段的植被识别和植被分类［5-6］。以往许多研究学者基于植被的光谱特征进行了植被分类［7-8］。由于许多地物光谱特征相似度高及环境因素的影响，许多地物难以分离。针对这一问题，有学者采用神经网络［9-10］、支持向量机［11-12］、随机森林［13-14］、决策树等［15-17］分类方法对植被进行分类并取得较好的分类结果。但是植被分布杂乱无序，且存在 “同物异谱”和 “同谱异物”的现象［18］，只利用光谱特征进行分类很难取得高精度的分类结果。有的学者采用结合纹理特征进行分类的方法并得到了更高的分类精度，证明了纹理特征有助于提高分类精度［19］。近年来，为了进一步提高植被的分类精度，有学者对植被的物候信息进行研究，发现有些植被的物候特征差异很大，利用多时相的遥感数据可以提高物候差异明显的植被区分能力［20］。不同植被的物候特征、纹理特征、光谱特征各不相同，因此其归一化植被指数（normalized difference vegetation index，NDVI）时间序列、纹理和光谱响应存在一定的差异。利用单一特征进行植被分类，精度往往较低。由于过去遥感数据源的限制，数据很难同时满足高时间分辨率、高空间分辨率和丰富的光谱信息这3个特点，例如MODIS数据具有较高的时间分辨率和丰富的光谱信息，但空间分辨率较低，只适合大尺度范围的植被分类［21］；Landsat系列数据的空间分辨率也较低，且其光谱信息中缺少植被十分敏感的红边波段［22］；而Sentinel-2数据则具有较高的时间分辨率、空间分辨率和丰富的光谱信息［23-24］，尤其是其多个红边波段对植被区分具有很大帮助。由于不同植被的物候期存在一定的差异，因此其NDVI时间序列变化特征也存在一定差异，利用多期遥感影像构建NDVI时间序列有助于将物候期差异明显的植被区分。本研究基于多时相的Sentinel-2遥感影像，首先选取植被生长旺盛时期的影像，计算NDVI并设定合适的阈值将研究区内的植被提取出来，剩余部分归为非植被。然后选取NDVI时间序列、最佳时相的Sentinel-2数据中10个波段的光谱反射率特征和主成分分析前3个分量的纹理特征作为分类特征，将研究区内的植被类型分为耕地、草地、常绿针叶林、落叶针叶林和落叶阔叶林五大类。

1 研究区概况

旺业甸实验林场（图1）位于内蒙古赤峰市喀喇沁旗西南部，地属燕山山脉北麓，属茅荆坝次生林区的一部分，为浅山丘陵地貌，地势西南高东北低，一般坡度为10°～30°，平均海拔为1300 m，最高峰是翠云峰， 海拔为 1890.9 m。41°21′～41°39′N， 118°09′～118°30′E，东西长 41 km，南北宽 20 km。林场总面积为253.07 km2，有林地面积为219.6 km2，林木总蓄积量为157万m3，森林覆被率为93.24%。在有林地面积中：天然次生林为115.33 km2，人工林为104.27 km2。其中人工落叶松Larix gmelinii面积为 45.33 km2， 油松Pinustabuliformis为 56.47 km2， 樟子松Pinus sylvestrisvar.mongholica为 2.15 km2，云杉Picea asperata为0.33 km2。旺业甸林区的生物资源复杂多样，具有很强的代表性、典型性、特有性及脆弱性，在全国的多样性保护中占有重要地位。林场有林地有落叶松，云杉，樟子松，红松Pinuskoraiensis，油松，黑桦Betula dahurica，白桦Betula platyphylla，山杨Populus davidiana，核桃楸Juglansmandshurica等12个优势树种。

2 研究数据

2.1 遥感数据与预处理

使用的Sentinel-2数据是从欧洲航天局（European Space Agency，ESA）的数据共享网站（https://scihub.copernicus.eu/dhus/#/home）免费下载。Sentinel-2A是欧洲航天局（ESA）于2015年6月23日发射的环境监测卫星，重访周期为10 d，其双星Sentinel-2B于2016年发射后，将重访周期缩短到5 d。Sentinel-2A数据的具体波段信息见表1。由于植被自身的物候特征，植被的光谱反射率在时间序列上有一定的变化规律，因此，在2017年内具有可用影像的月份均选取一景云量较少、质量较高的Sentinel-2影像，共9景，成像时间分别 为 2017-01-15， 2017-02-14， 2017-04-15， 2017-05-25， 2017-06-14，2017-09-22， 2017-10-27， 2017-11-21，2017-12-31。 利用 ESA提供的SNAP软件对数据进行辐射定标、大气校正、几何校正和图像裁剪。将像元较大的波段重采样为10 m×10 m的像元，从而使各个波段的像元大小一致。

图1 研究区位置示意图Figure 1 Location of the study area

表1 Sentinel-2波段信息表Table 1 Band information of Sentinel-2

2.2 外业数据

本研究中的外业数据来源于2017年9月研究区的实地调查，共选取300个适合本次研究的样点数据。实地调查遵循样本具有代表性、样本均匀分布的原则，在研究区内选取具有代表性的植被类型进行采样，且保证采集的样本点均匀分布在整个研究区内。样本点的数据属性包括树种类别和用手持GPS定位的经纬度坐标。由于采集得到的每个类别样本点数量不同，因此分别选取每个类别的训练样本，使其能够均匀分布于研究区内。共选取训练样本数据122个，其中非植被为18个，耕地20个，常绿针叶林22个，落叶阔叶林24个，落叶针叶林22个，草地16个。剩余178个样点数据作为验证样本，其中非植被为19个，耕地13个，常绿针叶林46个，落叶阔叶林41个，落叶针叶林43个，草地16个。

3 研究方法

采用以下3种方法对研究区进行植被分类：①最大似然法，以外业采集数据作为样本，基于Sentinel-2数据的 10个波段（B02， B03， B04， B05， B06， B07， B08， B08a， B11， B12）进行分类， 且各个类别的样本数据光谱特征均服从正态分布。②根据不同地物的归一化植被指数（NDVI）时序变化特征和光谱特征（NDVI时序＋光谱特征），选取最佳分类阈值，建立基于决策树的分层分类模型，将研究区分为非植被、耕地、草地、常绿针叶林、落叶针叶林、落叶阔叶林。③NDVI时序特征＋光谱特征＋纹理特征（Sentinel-2时序多特征），采用支持向量机分类器，经过多次实验核函数选取高斯函数，核函数系数γ为0.2，惩罚因子C为7，对在植被生长旺盛期，利用NDVI阈值提取出的植被进行植被类型分类，将其分为耕地、草地、常绿针叶林、落叶针叶林、落叶阔叶林。

NDVI对植被反应敏感，能够反映出植被与其他地物在可见光、近红外波段的明显差异。因此，首先选取植被生长茂盛时期的影像，计算其NDVI，设置合适的阈值将植被提取出来。NDVI计算公式如下： INDV=（B08-B04）/（B08＋B04）。 其中： INDV为归一化植被指数， B08为近红外波段的反射率，B04为红光波段的反射率。

3.1 构建NDVI时间序列

不同时间的NDVI能够反映植被的物候信息，且两者有一定的规律性，是区分植被类型的重要特征。例如，常绿针叶林常年保持绿色，NDVI在1 a中的变化较小；农作物在6-9月为生长期，此时耕地的NDVI也呈增长趋势，到9月NDVI达到最大值；落叶类植被一般在5月开始长叶10月开始落叶，NDVI也随之在5月开始上升10月开始下降。因此，分别计算9个时相影像的NDVI，构建6种地物的NDVI时序曲线，描述6种地物的动态变化过程。

3.2 最佳时相光谱反射率特征分析

根据地物的物候特征可知：在5-9月均为5种地物的最佳时相，因此选取这段时间内质量最好的2017-06-14影像绘制光谱反射率曲线进行单波段光谱分析。对于NDVI时间序列曲线相似、难以区分的地物，利用不同地物的光谱反射率差异可以增加这种相似地物的区分度，但是仍可能存在混分现象。例如落叶针叶林和草地，两者的NDVI时序变化规律相似，光谱反射率的变化规律也相似。因此，对纹理特征进行分析后，在最终分类时应用纹理特征以提高分类精度。

3.3 纹理特征分析

纹理特征是一种全局特征，它描述了物体表面的特性。与基于像素点的颜色特征不同，纹理特征需要在包含多个像素点的像素区域中进行统计计算。而且纹理特征常具有旋转不变性，对于噪声也具有较强的抵抗能力。本研究共选取了Sentinel-2数据中的10个波段进行分类研究，若提取每个波段的纹理特征并参与最终分类，将会产生大量数据，计算量太大。因此，在提取纹理特征前首先对数据进行主成分分析处理，然后选取前3个主成分分量用来进行后续的特征提取。1973年HARALICK提出的灰度共生矩阵是遥感图像中常采用的纹理特征提取方法［25-26］。灰度共生矩阵是通过研究图像上像元灰度的空间相关特性来描述纹理。通过研究分析选取5×5大小的移动窗口，利用灰度共生矩阵分别计算前3个分量的8种纹理特征，即均值、方差、同质性、对比度、差异性、熵值、二阶矩、相关性。

综上所述，共选取 43个特征用于分类，包括9个时相（2017-01-15，2017-02-14，2017-04-15，2017-05-25，2017-06-14，2017-09-22， 2017-10-27， 2017-11-21， 2017-12-31）的 NDVI特征，10 个波段（B02，B03，B04，B05，B06，B07，B08，B08a，B11，B12）的光谱反射率特征和24个纹理特征。

4 结果与分析

4.1 特征提取结果

由图2分析可得：①常绿针叶林在1 a中都是绿色植被，NDVI变化相对较小。在1，2，3，11，12月其他植物多处于非生长期，常绿针叶林的NDVI明显高于其他地物，容易被区分；②研究区内的耕地在5-9月属于生长期，NDVI处于上升趋势，但同一时期的其他植物也处于生长旺盛期，5种地物的NDVI均较高，相比之下耕地的NDVI值反而是最低的，容易被区分；③落叶阔叶林在5-9月属于1 a中生长最旺盛的时期，其NDVI值达0.80以上，高于同一时期的其他地物，可以进行区分。

通过分析图3发现：①耕地在B11波段的反射率上升相比于其他地物的反射率下降具有显著差异，利用B11波段很容易将耕地与其他地物区分出来；②草地在B06，B07，B08，B08a波段的反射率与其他地物均存在差异，可以利用这几个波段将草地与其他地物区分；③落叶针叶林在B05波段的反射率低于落叶阔叶林，而在B06波段的反射率却高于落叶阔叶林，可以利用这一特征将两者进行区分。

4.2 分层分类阈值设置

通过分析典型地物的NDVI时间序列曲线与典型地物的光谱反射率曲线，得出结合NDVI时序与光谱特征分类方法的最佳分类阈值。在2017-09-22的影像中，INDV＜0.37为非植被；在2017-12-31的影像中，INDV＞0.39为常绿针叶林；在2017-06-14的影像中，B11波段的反射率B11＞25%为耕地（图3）；在2017-06-14的影像中，B08波段的反射率B08＞21%为草地；在2017-05-25的影像中，INDV＞0.8为落叶阔叶林；在2017-05-25影像中，0.7＜INDV＜0.8为落叶针叶林。

图2 典型地物的NDVI时间序列曲线Figure 2 NDVI temporal profiles of typical land cover types

图3 典型地物的光谱反射率Figure 3 Spectral reflectance curve of typical land cover types

4.3 分类结果

通过分析比较3种方法的分类结果（图4和表2）发现：①最大似然法。非植被与耕地在分类结果中的面积明显偏少，存在大量的漏分；草地在分类结果图中没有明显的表现出来，也存在严重的漏分；常绿针叶林、落叶针叶林和落叶阔叶林在分类结果完全是分散分布且交错分布的，这与实地调查情况不符，存在明显的混分。②NDVI时序+光谱特征方法。非植被分类情况与基于多特征的分类结果一致，且明显可以看出包含了一些道路，与实际调查情况相符；耕地和草地的分类结果与其他2种方法相比面积较多，尤其是草地的面积过多，与常绿针叶林、落叶针叶林和落叶阔叶林产生混分。③NDVI时序+光谱特征+纹理特征方法。非植被与耕地的面积和分布情况均符合实际情况；常绿针叶林主要分布在耕地周边，集中在研究区的中部；落叶阔叶林主要分布在研究区的西部和西南位置，少数分散分布于整个研究区，与落叶针叶林交错分布；草地在研究内呈现为一些散碎的斑块，且面积相对较小。整体的分类结果与实地调查情况有较高的一致性。

4.4 精度评价

图4 3种方法的分类结果示意图Figure 4 Classification results of three methods

利用实地调查采集的178个样点数据作为验证样本，其中非植被19个，耕地13个，常绿针叶林46个，落叶阔叶林41个，落叶针叶林43个，草地16个。对比3种方法的分类结果与实地调查结果，并计算出每种方法的类别分类精度及总体分类精度。由表3可见：基于NDVI时序+光谱特征+纹理特征的多特征分类总体精度达87.64%，分别比最大似然法和NDVI时序＋光谱特征的分类总体精度提高了15.73%和14.61%。而且基于NDVI时序+光谱特征+纹理特征的多特征分类中单个类别的分类精度也均较高，其中常绿针叶林的分类精度最高达到了95.65%，这是由于利用了NDVI时序特征，可以通过分析不同植物在1 a中的NDVI变化趋势，选择出最佳的分类时相将常绿针叶林与其他地物区分。耕地的分类精度也比较高，其NDVI在1 a中的变化具有一定的规律性，且耕地在B11波段的反射率与其他地物有明显差异，利用这些特征可与其他地物很好地区分。NDVI时序＋光谱特征的分类精度与NDVI时序+光谱特征+纹理特征的分类精度相比，非植被这一类别的分类精度相同，但其他类别的分类精度则较低。这是由于NDVI时序变化和光谱特征代表的仅仅是一个平均值，并不能完全代表地物的特征变化，因此，根据平均值来选取阈值分层分类的精度较低。

表2 3种分类方法的地类面积Table 2 Land area of three methods

表3 3种分类方法的精度评价Table 3 Accuracy evaluation of three methods

5 结论与讨论

本研究基于多时相的Sentinel-2遥感影像，选取NDVI时间序列、最佳时相的Sentinel-2数据中10个波段的光谱反射率特征和主成分分析前3个分量的纹理特征作为分类特征，利用支持向量机分类器，将研究区内的植被类型分为耕地、草地、常绿针叶林、落叶针叶林和落叶阔叶林。得出结论如下：①基于多特征的分类总体精度达87.64%，Kappa系数为0.85，分别比最大似然法的分类结果提高了15.73%和0.20。因此，结合多种特征的分类方法有助于提高分类精度。②通过构建NDVI时间序列，分析地物的NDVI变化趋势，选择常绿针叶林提取的最佳时相进行提取，可以大大提高常绿针叶林的分类精度。NDVI的时序特征对于植被的区分具有很大帮助。③耕地在Sentinel-2影像的B11波段上的反射率与其他地物有显著差异，而且耕地在NDVI时间序列上也与其他地物有明显差异，利用这些特征极大提高了耕地的分类精度。④采用分层分类的思想，首先将研究区内的植被提取出来，然后再对植被进行类型识别，可以排除非植被因素的干扰，有效提高植被类型的分类精度。

分类结果中常绿针叶林与耕地均得到了很高的分类精度，这是由于本研究考虑了植被的物候信息，构建了不同植被的NDVI时间序列变化曲线，并充分利用植被的物候特征进行分类。分类结果中也显示草地的分类精度是所有类别中精度最低的。一方面原因是，在研究区内草地分布比较分散而且比较破碎，在实地调查中采集的草地样本点相比其他类别样本数量较少，分类时训练样本的数量相比其他类别也比较少，从而导致了分类精度较低。另一方面，草地的NDVI时间序列变化规律与落叶针叶林、落叶阔叶林相似，其光谱反射率的变化规律与落叶针叶林也相似，从而导致草地极易与落叶针叶林、落叶阔叶林混分，使其分类精度较低。本研究只选取了NDVI时序特征、光谱特征和纹理特征3种特征作为分类特征，而地形特征中也包含了丰富的信息，可以充分利用地形特征进行深入研究。