基于哨兵2时间序列组合植被指数的作物分类研究

doi:10.11873/j.issn.1004-0323.2020.3.0702

基于哨兵2时间序列组合植被指数的作物分类研究

谷祥辉^,¹, 张英², 桑会勇², 翟亮^,², 李少军³

1.山东科技大学测绘科学与工程学院，山东青岛 266590

2.中国测绘科学研究院, 北京 100830

3.新疆维吾尔自治区测绘科学研究院, 新疆乌鲁木齐 830002

Research on Crop Classification Method based on Sentinel-2 Time Series Combined Vegetation Index

Gu Xianghui^,¹, Zhang Ying², Sang Huiyong², Zhai Liang^,², Li Shaojun³

1.College of Geomatics, Shandong University of Science and Technology, Qingdao 266590, China

2.Chinese Academy of Surveying and Mapping, Beijing 100830, China

3.Xinjiang Academy of Surveying and Mapping, Xinjiang Uygur Autonomous Region Urumqi, Urumqi 830002，China

通讯作者: 翟亮（1980-），男，河南安阳人，博士，研究员，主要从事地理国情监测与分析方面的研究。E⁃mail:zhailiang@casm.ac.cn

收稿日期: 2019-07-01 修回日期: 2020-04-11 网络出版日期: 2020-07-03

基金资助:

中国测绘科学研究院基本科研业务费项目. 7771728. 7771820. AR1928. AR1931
自然资源部地球观测与时空信息科学重点实验室开放基金. AR191902

Received: 2019-07-01 Revised: 2020-04-11 Online: 2020-07-03

作者简介 About authors

谷祥辉（1993-），男，山东枣庄人，硕士研究生，主要从事地理国情遥感监测、定量遥感等方面的研究E⁃mail：1090441654@qq.com , E-mail：1090441654@qq.com

摘要

时间序列是一种常用的物候研究方法。为充分利用哨兵2数据在红边波段的丰富信息，本文利用多种植被指数组合成时间序列进行作物分类。将NDVI、EVI、红边NDVI三种植被指数进行组合，构建时序植被指数图像，然后使用支持向量机、随机森林、CART决策树和最大似然4种不同的算法对四种作物、三种林草、裸露地表、水体进行分类。原始分类结果中，总体精度最高的随机森林为87.92%，最低的最大似然为80.07%，在分类细节上，随机森林和支持向量机的边界最清晰，4种分类结果中，农作物的分类精度均高于其他地类，仅次于水体的精度，误差主要来自三种林草的混分，表明时间序列组合植被指数用于农作物分类是可行的。

关键词： 植被指数 ; 时间序列 ; 遥感 ; 农作物分类 ; 哨兵2

Abstract

Time series is a widely used phenological research method. A new time series vegetation indices which takes full advantage of the red edge information of Sentinel 2 data were used for crop classification to improve the classification accuracy. The NDVI, EVI, and red edge NDVI were combined to construct a time series vegetation index image. Then, four different algorithms (support vector machine, random forest, CART decision tree and maximum likelihood) were used to classify four crops, three forest grasses, bare land, and water bodies. Among the original classification results, the random forest with the highest overall accuracy is 87.92%, and the maximum likelihood with the lowest overall accuracy is 80.07%. In the classification details, the boundaries of random forest and support vector machine are the clearest. Among the four classification results, the classification accuracy of crops is higher than other land types, just smaller than water body. The error mainly comes from the mixture of three forests. It indicates that the time series combined vegetation index is feasible and accurate for crop classification.

Keywords： Vegetation index ; Time series ; Remote sensing ; Crop classification ; Sentinel 2

PDF (5926KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

谷祥辉, 张英, 桑会勇, 翟亮, 李少军. 基于哨兵2时间序列组合植被指数的作物分类研究. 遥感技术与应用[J], 2020, 35(3): 702-711 doi:10.11873/j.issn.1004-0323.2020.3.0702

Gu Xianghui, Zhang Ying, Sang Huiyong, Zhai Liang, Li Shaojun. Research on Crop Classification Method based on Sentinel-2 Time Series Combined Vegetation Index. Remote Sensing Technology and Application[J], 2020, 35(3): 702-711 doi:10.11873/j.issn.1004-0323.2020.3.0702

1 引言

随着城市化的进行，农作物种植面积、熟制均发生剧烈变化，山区耕地利用边际化现象以及华北平原冬小麦休耕现象尤为明显^[1,2,3]。准确监测农作物种植面积、种植结构及空间分布信息对于粮食安全、农业生态系统变化以及政策制定等都具有重要的现实意义。传统普查与调查手段费时费力且难以准确反映农作物的空间分布信息。遥感数据其快速大范围的优势使得遥感手段成为监测作物种植信息的有效手段，且发展出许多监测方法^［4］。

多时相遥感分析方法应用广泛，不同农作物具有不同的生长周期，因此物候特征在农作物识别中扮演着重要角色。植被指数反映了植被对光谱波段的响应特征，作物的植被指数在不同的生长阶段由于色素和结构的不同而不同；不同的植被指数侧重点不同，对同一种作物有不同的反映；不同作物整个生长期的植被指数变化也不同；因此，多种植被指数结合时间序列信息能够区分作物的生长期以及生长状态，进而区分不同作物^[5]。已有学者利用MODIS时间序列植被指数进行农作物遥感识别，并取得了较好的分类结果，验证了利用时序植被指数进行农作物识别与提取的可行性，应用的植被指数包括归一化植被指数（NDVI）^{[6,7,8,9,10,11,12,13]}和增强型植被指数（EVI）^{[14,15,16,17]}，EVI加入了蓝光波段对土壤的影响进行改正，更好地识别植被，但空间分辨率降低到500 m。由于MODIS数据空间分辨率较低的限制，影像中存在较多混合像元，难以识别小面积种植的作物，并且仅使用NDVI或者EVI并不能充分利用MODIS高光谱的优势。

由于MODIS空间分辨率的限制，有学者使用环境星、高分1号WFV、Landsat 8等中高分辨率影像进行作物识别^[18,19]。张荣群等使用5个时期的Landsat 8影像及其他辅助数据，通过分析待分类作物 NDVI 曲线图的变化特征，确定分类知识规则，提取出各类作物。但Landsat 8重访周期较长，再加上云层的影响，难以获取完整的时序影像，因此，有学者探索使用相似度较好的多源数据来提高时间分辨率^[20,21,22]。

哨兵2卫星每5天过境一次，具有较高的时间分辨率，容易获得完整的时间序列影像，空间分辨率较高，并且在红边范围内具有较多波段，适用于植被监测^[23]，国内外均有学者使用单景哨兵2数据实现作物识别^[24,25]。Vrieling等^[26]对比了相机测量计算出的绿色色度坐标序列和哨兵2数据植被指数序列，证明了哨兵2数据提供植被物候信息的能力^[26]。

本研究主要使用哨兵2时间序列植被指数进行作物分类，以评估哨兵2时间序列组合植被指数在作物精细分类方面的潜力及精度。具体来说，使用哨兵2数据丰富的光谱信息构建不同的植被指数，并对每个时相选取合适的植被指数进行组合，利用多种分类方法对时间序列组合植被指数对作物进行分类。为研究精细分类的能力，本文从两个方面进行试验，一是研究区地块较为零散，作物种植尺度较小，其次作物分类较细，将草木类经济作物纳入分类体系。

2 研究区与数据处理

2.1　研究区

如图1所示，本文研究区行唐县（38°20′34″~38°42′39″ N，114°09′56″~114°41′52″ E）位于河北省石家庄市，属太行山东麓浅山区，主要有低山、丘陵、平原三种地形，属于暖温带半湿润大陆性季风气候区，雨热同季，年均降水量450.4 mm。县域面积966 km²,2017年石家庄市统计年鉴显示行唐县总播种面积619 km²，主要作物有冬小麦、玉米、大豆、花生及牧草。行唐县内85%左右的地表为种植土地和林草覆盖，种植土地面积超过600 km²，林草覆盖超过200 km²，房屋建筑、道路、构筑物、人工堆掘地等人工地表的面积约为100 km²，荒漠与裸露地表面积不足5 km²，水域面积不足15 km²。本文获取研究区内2017年2月至10月的哨兵2数据，选择对应物候期的植被指数进行组合，利用组合植被指数对作物的类别和熟制进行提取，以期减少作物提取对样本的依赖性，提高在地形复杂、作物类型丰富地区的分类精度。

图 1

新窗口打开| 下载原图ZIP| 生成PPT

图 1 研究区概况

Fig.1 Location of the study area

2.2　数据获取及预处理

使用的哨兵2多光谱数据来自欧洲航天局（ESA）（https://scihub.copernicus.eu/dhus/#/home），哨兵2包括哨兵2A和哨兵2B，单颗卫星的重访周期为10天，双星提高到5天，空间分辨率为10 m；哨兵2数据具有13个波段，覆盖可见光到近红外部分，除了常用的可见光近红外波段外，还有4个红边波段，且波段宽度较小，光谱分辨率较高（表1）。

表 1 哨兵2数据空间与光谱分辨率

Table 1 The spatial and spectral resolution of Sentinel 2

波段	S2A		S2B		空间分辨率/m	备注
波段	中心波长/nm	波段宽度/nm	中心波长/nm	波段宽度/nm	空间分辨率/m	备注
2	492.4	98	492.1	98	10
3	559.8	45	559	46
4	664.6	38	664.9	39
8	832.8	145	832.9	133
5	704.1	19	703.8	20	20	窄波段红边
6	740.5	18	739.1	18
7	782.8	28	779.7	28
8a	864.7	33	864	32
11	1 613.7	143	1 610.4	141	20
12	2 202.4	242	2 185.7	238	20
1	442.7	27	442.2	45	60
9	945.1	26	943.2	27
10	1 373.5	75	1 376.9	76

新窗口打开| 下载CSV

本研究共获取了2017年2月至10月共15期无云影像（表2），涵盖全部作物的生长季，影像的时间间隔最大为30天，最小为5天，每个月至少有1期影像，作物生长期关键月份5~8月、10月均有2期影像，影像的具体信息如表2所示。使用ESA提供的sen2cor工具对影像进行辐射定标、大气校正、重采样、格式转换以及波段合成等预处理，生成10 m分辨率的多光谱数据，并在ENVI软件中进行空间配准和矢量裁剪，为方便进行图像处理操作，裁剪时保留矢量轮廓外的部分。

表 2 主要农作物发育期

Table 2 The development periods of major crops

获取日期	儒略日*	农作物所处发育期
获取日期	儒略日*	冬小麦	夏玉米	春玉米	大豆	花生
2017/2/27	58 d	返青
2017/3/9	68 d	起身
2017/3/29	88 d	起身
2017/4/18	108 d	拔节		播种
2017/5/18	138 d	开花		出苗
2017/5/28	148 d	乳熟
2017/6/7	158 d	成熟		抽雄
2017/6/27	178 d	收割	播种	吐丝	播种	播种
2017/7/7	188 d		出苗		出苗	出苗
2017/7/12	193 d		拔节	乳熟		开花
2017/8/6	218 d		抽雄		花芽	下针
2017/8/16	228 d		吐丝	成熟	结荚	结荚
2017/9/20	263 d		乳熟	收割	成熟	收割
2017/10/20	293 d	播种	收割		收割
2017/10/30	303 d	出苗

注：*与标准天文记法略有不同，在研究时间序列时，为更好地体现时间延续性，取当年的1月1日记为起始日期，逐日向后顺延，直至365或366，意为该年中的第几天。

新窗口打开| 下载CSV

3 方法研究

3.1　样本选取

本研究使用的样本数量较少，但样本分类较细且分类精确。选择训练样本时充分利用多时相数据和不同组合的假彩色数据，配合作物物候信息进行判断，对于物候期相近且图像表现相似的作物，参照Google Earth 高分辨率遥感影像进行判断，为提高样本和分类精度，在山地和平原地区表现不同的作物，作为两种样本进行选取，分类后再进行合并。样本设置在均一的农田上，以获取纯净像元。

3.2　构建时序植被指数

根据数据源的中心波长及波段宽度（BW）（表 1）选择合适的植被指数，对比以下三种植被指数（表3），选择最优指数进行组合，这三种植被指数的值域范围均为-1~1，无需再次进行归一化处理，便于组合。NDVI属于宽波段植被指数，对红光波段和近红外波段的波长范围要求较小，适用范围较广，但在植被茂密时容易饱和；EVI加入了蓝光波段对土壤背景进行校正，有助于植被信息的提取^{[14,15,16,17]}；红边植被指数（NDVI705）是一种高光谱植被指数，要求输入窄波段，该指数使用植被反射率急剧上升的红—近红外波段范围（红边），对植被状态高度敏感^[27,28]。

表3 3种植被指数及公式

Table 3 Vegetation index and expression

指数	计算方法	波幅宽度要求
NDVI	（ρNIR-ρR）/（ρNIR+ρR）	宽
EVI	2.5*（ρNIR-ρR）/（ρNIR+6ρR-7.5ρB+1）	宽
NDVI705	（ρ750-ρ705）/（ρ750+ρ705）	窄

新窗口打开| 下载CSV

计算方法中，ρNIR为近红外波段的反射率，ρR为红光波段的反射率，ρB为蓝光波段的反射率，ρ750和ρ705分别为高光谱波段中心波长为750 nm和705 nm波段的反射率。由于所需处理数据较多，植被指数计算在ENVI/IDL环境下批量完成。

将三种指数进行组合以充分应用数据的光谱信息，同时为减少不同指数间的冗余，每个时相选择最优指数进行组合。具体而言，分别构建三种植被指数的时序数据，应用所选纯净样本点对三种植被指数进行筛选，样本在三种时序植被指数中的整体趋势相似，但全部样本的NDVI值均高于NDVI705和EVI，EVI对8月植被和非植被的区分度高于其他指数，NDVI在所有时相的标准差均大于另外两种指数，NDVI705在大多数时相的标准差均较小，EVI在10月30日标准差最小。因此，使用8月和10月30日的EVI数据和其他时相的NDVI705数据构建时序植被指数。

3.3　植被指数曲线及可分性

研究区内的典型地表类型包括耕地（冬小麦-玉米大豆双熟制耕地、冬小麦-花生双熟制耕地、夏玉米单熟制耕地、春玉米单熟制耕地）、林草类（牧草、天然草地、林地）、非植被地表（裸露地表、水体），其中非植被地表的植被指数全年变化较小，裸露地表的植被指数序列接近0，水体远小于0。林草类生长期较长，在88 d进入绿色生长期，至138 d前后已比较茂盛，该时期其植被指数逐渐增加，随后保持较高的植被指数不变，直到263 d前后开始枯萎，植被指数逐渐减少，其中林地整体高于天然草地，牧草由于人工种植，在第二个时期高于天然草地，接近林地，但生长最茂盛的时期由于收割而急剧下降再急剧上升，如图2所示。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 林草类时序曲线

Fig.2 Time series curve of forest and grass

耕地的植被指数波动明显，随种植作物类型与熟制而变化，如图3所示，双熟制耕地具有双峰，在冬小麦时期具有相同的变化曲线，秋粮收获时期则不同，夏玉米收割时间比花生略晚，夏玉米单熟制耕地在178 d之前植被指数较低且变化较小，其后与冬小麦-夏玉米表现一致，春玉米生长季较夏玉米整体前移，138 d前后进入生长季，263 d已完成收割。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 农作物时序曲线

Fig.3 Time series curve of crops

在组合植被指数序列中，耕地、林草、非植被地表特征鲜明，识别度极高；不同的耕地其峰值数量、快速上升期和快速下降期存在差异，容易识别；林草中牧草存在谷值，容易从林草中识别。综上所述，依据组合植被指数对作物进行分类是可行的。

3.4　分类方法

本研究分别采用4种算法对组合后的植被指数进行分类，包括最大似然法（MLC）和3种机器学习算法（支持向量机（SVM）、CART决策树、随机森林（RF）），以验证组合植被指数用于作物提取的可行性。

最大似然算法是基于统计的点估计法，假设每一个波段的每一类统计都呈正态分布，用判别函数计算像元属于每个类别的可能性，将其划分到可能性最大的那一类里边，可能性小于阈值的像元不分入该类。该方法分类速度较快，人为干预少，分类结果稳定^[29]。

支持向量机是一种基于统计学习理论和结构风险最小化原则的机器学习算法，它将非线性分类问题转化为高维线性问题，并在高维特征空间中构造线性判别函数，同时引入核函数减少运算量。支持向量机可以自动寻找对分类有较大区分能力的向量，具有较高的分类精度，本文中的SVM使用径向基函数（RBF）作为核函数，Gama in Kernel Function参数设置为0.067。

分类回归树算法（CART）是一种决策树（DT）构建算法，由Breiman于1984年提出，该算法将样本集划分为两个子集，向下构建两个分支，最终生成一个简洁的二叉树，CART算法中采用基尼系数来选择最佳分割阈值，使得每个子节点中的要素都具有相同的特征^[11]。CART算法生成的决策树没有人工干预，简化了构建DT的步骤，减少了主观误差。

随机森林（RF）由若干个DT组成，随机森林随机地从样本中选择一部分确定分类树节点的阈值，随机地生成成百上千个DT，然后选择重复度最高的DT作为最终结果。RF有两个重要的参数，一是树节点的预测变量数目 Number of Feature，设置为 $\sqrt[]{p}$ ，其中p是输入向量的维度。另一个是决策树（DT）的数量，本文中样本特征明显，将其设为100。

4 结果与分析

4.1　精度验证

将研究区划分为623个格网，使用随机函数选取9个格网作为样方进行验证，样方位置如图4所示。通过高分影像目视判读、实地调查等方式解译出样方的地表类别，以Kappa系数、总体精度、制图精度、用户精度对4种分类结果进行验证。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 样本区位置

Fig.4 Distributions of ground samples in study area

4.2　农作物分类结果及精度评价

4种原始结果的分类精度如表4所示，4种分类结果的总体精度关系为RF>SVM>CART>MLC，RF和SVM的精度相近，RF略高于SVM，但在主要农作物类别上，SVM的制图精度整体高于RF方法。由于波段较多，数据存在冗余，MLC分类的精度小于其他分类方法，主要表现在林地、天然草地、牧草等类别上，主要农作物仍具有较准确的分类结果。4种分类方法的原始分类结果中，农作物、水体的精度均较高，裸露地表略低，林草类地物最低，表明组合时间序列植被指数能够有效地用于农作物分类。

表 4 原始分类结果分类精度

Table 4 Classification accuracy of original classification results

类别	指标	分类方法
类别	指标	SVM	RF	CART	MLC
冬小麦玉米大豆	制图精度	89.68	89.54	84.73	79.84
冬小麦玉米大豆	用户精度	90.24	89.88	89.89	92.84
冬小麦花生	制图精度	88.39	88.21	79.30	83.00
冬小麦花生	用户精度	95.60	93.19	86.74	93.70
春玉米	制图精度	93.12	91.48	89.29	91.06
春玉米	用户精度	82.88	87.61	79.56	67.65
夏玉米	制图精度	87.79	87.44	87.04	79.04
夏玉米	用户精度	91.17	92.04	87.03	95.26
林地	制图精度	85.03	77.70	76.08	77.70
林地	用户精度	67.61	67.49	62.16	43.91
天然草地	制图精度	88.86	87.71	83.85	76.21
天然草地	用户精度	91.24	93.72	92.15	94.18
牧草	制图精度	82.34	81.76	72.76	86.26
牧草	用户精度	78.02	70.71	61.64	40.54
水体	制图精度	98.26	99.21	99.30	90.64
水体	用户精度	99.93	99.99	97.58	100.00
裸露地表	制图精度	79.63	85.54	79.99	75.51
裸露地表	用户精度	87.38	81.90	78.15	84.56
总体精度		87.90%	87.92%	84.25%	80.07%
Kappa系数		0.857 9	0.858 1	0.815 4	0.769 3

新窗口打开| 下载CSV

从分类结果细节来看，RF方法分类结果聚合程度高，边界清晰。如图5所示，裸露地表、林地、牧草、冬小麦玉米大豆的边界十分整齐；SVM方法中裸露地表地表、林地、冬小麦玉米大豆的边界较十分整齐，但林地和牧草相互交错，部分林地误分为牧草；CART方法牧草、冬小麦玉米大豆、林地的边界较清晰，但裸露地表中包含许多零散的地类；MLC方法聚合程度较高，但边界不清晰，林地有错分为牧草的现象。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 分类细节

Fig.5 Classification details

原始分类结果中存在孤岛，即分类中出现少数其他类别的现象，实际地物往往集聚在一起，极少出现孤岛，通常是噪声引起的，为去除孤岛，应用卷积模板对原始结果进行处理，将窗口内像元数量少于2个的类别归入背景类别中，以消除孤岛。消除孤岛后，4种分类结果的总体精度和Kappa系数如表5所示。

表 5 消除孤岛后的精度

Table 5 Classification accuracy after islands removed

指标	分类方法
指标	SVM	RF	CART	MLC
总体精度	88.37%	88.52%	85.89%	81.03%
Kappa系数	0.863 3	0.865 0	0.834 4	0.780 1

新窗口打开| 下载CSV

原始分类结果如图6所示，4种分类方法的地类分布整体一致，行唐县城周边平原地区的作物以冬小麦玉米大豆为主；冬小麦花生主要位于行唐县城的东部及东北部，分布比较零散；春玉米主要分布于东部边界处，面积较小；秋玉米主要分布在平原向山区过渡地区，该地区大量弃耕冬小麦，只种夏玉米；牧草呈块状分布于村镇周边。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 原始分类结果

Fig.6 Original classification results

5 结语

本文基于哨兵2号时间序列数据，对行唐县主要农作物进行分类提取。首先根据哨兵2数据的波段信息选取了三种归一化的植被指数，并根据实测样本信息逐时相选取最佳植被指数，构成时间序列组合植被指数。鉴于该时序数据时间间隔不等，且序列较少，改用不变目标法对时序曲线进行校正。最后通过4种分类方法进行分类，提取农作物类型。结果表明：

在时序数据中序列较少，间隔不等的情况下，采用不变目标法进行时序数据校正是可行的，能够在高度保留特征信息的情况下修正时序曲线。

在组合植被指数序列中，耕地、林草、非植被地表的特征明显，区分度极高；耕地中的4种主要作物之间的时序曲线在峰值、快速上升期和快速下降期差异明显。

运用4种分类方法对组合时序数据进行分类，RF和SVM方法在作物上的分类精度最高，若考虑林草，则RF方法高于SVM方法，CART方法精度次之， MLC方法精度最低。分类细节上，RF方法分类结果聚合程度最高，边界清晰。

4种分类方法的原始分类结果总体精度均较高，最低的MLC为80.07%，最高的RF为87.92%，对分类结果进行消除孤岛后，MLC提高到81.03%，RF提高到88.52%。影响总体精度的主要为裸露地表、林地、牧草之间的错分漏分现象。表明时间序列组合植被指数用于农作物分类是可行的。

行唐县的作物在平原地区以冬小麦玉米大豆为主，在浅山区则为夏玉米，春玉米和花生较少且分布零散，牧草呈块状分布于村镇周边。

本研究主要利用哨兵2时间序列数据的植被指数进行作物分类研究，探究了在多种植被指数组合、时间序列不规则的情况下进行作物分类的能力，为减少数据冗余，对每个时相选择最优指数，但从分类结果来看，仍有部分数据冗余，可能存在于非关键时相中，有待进一步研究。其次，本文基于像元进行分类，在平原地区对作物分类时能较好地区分，但涉及到山区的林草时，精度有所降低，可以尝试应用面向对象分类进一步提高。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Wang

Xue

, Li

Xiubin

, Xin

Liangjie

, et al.

Driving Factors on Winter Wheat Abandonment by Farmers in North China Plain and Its Policy Implication

[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(9):248-257.