基于时序特征的玉麦轮作制耕地提取方法研究

李百红 彭勃 董超

摘要 [目的]该研究对小麦、玉米轮作制下耕地的面积与分布有重要意义。[方法]基于目视判别收集样本点和GEE平台，分析地物的NDVI、EVI等指数时序特征，开发特征识别算法提取2018年玉麦轮作区面积分布，同时与CART算法监督分类结果进行精度分析。[结果]与基于多光谱的监督分类相比，基于时序特征的识别算法准确率、精准率、F1 Score和Kappa系数分别提高了0.121、0.110、0.136和0.246。[结论]基于时序特征的特征识别算法可有效识别玉麦轮作区，该算法可为轮作区农业生产提供基础数据支持。

关键词 GEE;时序;归一化指数;提取算法

中圖分类号 S 126  文献标识码 A  文章编号 0517-6611（2021）19-0214-04

doi：10.3969/j.issn.0517-6611.2021.19.056

开放科学（资源服务）标识码（OSID）：

Study on Extraction Method of Wheat-Maize Rotation System Based on Time Series Characteristics

LI Bai-hong1， PENG Bo2， DONG Chao2

（1.Taian Natural Resources and Planning Bureau of Shandong Province， Taian， Shandong 270000;2.College of Information Science and Engineering， Shandong Agricultural University， Taian， Shandong 271018）

Abstract [Objective] It is important to study the area and distribution of cultivated land under wheat and maize rotation system. [Method] Based on the collection of sample points and GEE platform， we analyzed the NDVI， EVI and other index time series features of ground objects， developed feature recognition algorithm to extract the 2018 wheat-maize rotation area distribution， and conducted precision analysis with the supervised classification results of CART algorithm. [Result] Compared with the multi-spectral supervised classification， the accuracy， F1 Score and Kappa coefficient of the time-series feature recognition algorithm were improved by 0.121， 0.110， 0.136 and 0.246， respectively. [Conclusion] The feature recognition algorithm based on time series features could effectively identify wheat-maize rotation area， and the algorithm could provide basic data support for agricultural production in rotation area.

Key words GEE;Time series;Normalized index;Extraction algorithm

基金项目 山东省重点研发项目（14032761，140380198）。

作者简介 李百红（1985—），女，山东青州人，工程师，硕士，从事土地资源利用与信息化工作。*通信作者，讲师，博士，从事农业信息技术研究。

收稿日期 2021-02-26

快速准确获取作物种植的分布信息有利于实现对耕地的精确管理、产量估算、制定科学的粮食管理政策[1]。当前，作物识别多基于遥感影像光谱差异，如熊勤学等[2-3]利用NDVI提取小麦的种植面积，路中等[4-6]提取玉米的种植面积，刘佳等[7-8]完成了中国冬小麦制图和大规模春小麦提取研究。受限于传统数据处理工具的低效，这些研究只使用了一种指数提取地物信息，难以发挥高分辨率多光谱遥感数据的优势。多时相的遥感数据进行分类不仅可以提取地物物候特征，而且可有效避免单一时相存在的异物同谱和同谱异物状况[9]。Google Earth Engine（GEE）平台可对大尺度区域进行分析，还支持多源多时相海量数据的协同研究[10]。鉴于此，笔者选取华北平原小麦-玉米轮作区，基于GEE平台，利用多源遥感数据分析地物样本点光谱指数的时序特征，依据各类地物的时序特征设计有效的识别算法，通过精度评价评估识别算法的优劣。

1 研究区域与数据

1.1 研究区概况

华北平原小麦-玉米轮作区（32°16′～39°38′N，112°50′～118°57′E）总面积约221 838 km2，跨河北、河南、山东3省（图1）。研究区大部分平原区为栽培植被，其中以冬小麦、夏玉米轮作最为常见，主要分布在河南北部和山东西北的引黄灌溉区、太行山前平原水源条件较好的区域。研究区内小麦10月中上旬播种，次年3月拔节，6月上旬收割;6月初播种玉米，9月下旬收割[11]。

1.2 数据及预处理

1.2.1 卫星影像选取。

研究所用数据包括：Landsat7 SR（Surface Reflectance，SR）数据集;Landsat8 SR数据集;MOD13Q1.006数据集的植被指数（MOD13Q1）第6版数据;Sentinel-1 SAR GRD数据集。研究影像数据来自GEE平台（https：//earthengine.google.com/）。

1.2.2 影像预处理。

影像预处理包括：①使用研究区对影像统一进行掩膜;②利用Landsat7和Landsat8影像中由CFMASK算法生成的質量评估波段去云[12];③对Landsat7和Landsat8计算归一化光谱指数。NDVI 计算公式如下：

NDVI=（NIR-R）/（NIR+R）（1）

根据NIR、R与B波段计算EVI计算公式如下：

EVI=2.5×（NIR-R）/（NIR+6×R-7.5×B+1）（2）

根据G与MIR波段计算MNDWI，计算公式如式下：

MNDWI=（G-MIR）/（G+MIR）（3）

将Landsat7和Landsat8的影像合成指数的数据集，以弥补Landsat系列卫星回访周期较长的缺点。

1.3 样本点与验证点

基于时间序列将研究区地物类别分为小麦-玉米轮作区、林地、城镇建设用地、水域和其他耕地类型。利用在线高分辨影像及历史时间轴，结合研究区内各区县的相关文献资料和统计报告，通过目视的方法选取地物样本点和验证点（表1）。

2 时序数据分类方法

2.1 时序特征的组合

基于指数数据集生成样本点指数时序曲线，取同类指数的平均值生成各类别指数年积日（Day of year，DOY）时序曲线。指数时序曲线可反映类别特点，研究基于NDVI、EVI和MNDWI共3种指数的时序特征，筛选区分地物类别的时序参数。

2.2 轮作区耕地的识别

采用分层分类加决策树的方法，基于指数的时序特征对研究区中的玉麦轮作区外进行分层掩膜[13]，掩膜顺序为：水体、建筑、森林、轮作区。由于Landsat重访周期较长，受天气影响较大，在识别算法中使用MODIS的NDVI产品来对轮作耕地进行最后的掩膜。

2.3 监督分类

利用目视样本点作为训练数据集，通过GEE提供的云评分算法制作研究区最小云量合成影像，并基于影像的蓝光、绿光、红光、近红外、短波红外等7个波段的合成影像，使用分类回归树（Classification And Regression Tree，CART）算法完成轮作区耕地提取[14]。

2.4 精度评价

利用独立的验证点对结果进行分类评价，将非轮作区耕地的地类视为一类与轮作区耕地进行二分类评价，通过准确率、精准率、F1 Score以及Kappa系数对提取精度进行评价。

3 结果与分析

3.1 时序特征分析

由图2a可知，水域MNDWI数值全年都为正值，且均值大于0.4明显区别于其他地类，基于MNDWI提取水域效果最优。图2b中森林和耕地在EVI上有较为明显的差异，在DOY120～273 d，森林的EVI均值大于2，明显区别与小麦收割后裸地。图2c中建设用地的NDVI曲线较为平缓，全年的数值偏低，且年均值小于0.3，以此来区分建筑用地。图2d中玉麦轮作区在NDVI时序上具有明显双峰特征，明显区别于其他地类。儒略日第79～154天是第1个峰，第198～271天是第2个峰。

3.2 监督分类结果

为比较算法的精度，采用监督分类与识别算法进行比较。二分类结果如图3所示，分类结果中轮作耕地分布过多，明显与调查情况不一致。

二分类混淆矩阵如表2所示，可见分类结果中，有大量轮作耕地被误分为非轮作耕地，同时也有大量非轮作耕地被误分为轮作耕地，2种地类相互混淆，提取精度较低。

监督分类结果中主要城市、河流湖泊分类准确，但北部只有少量其他耕地被正确识别，且分类结果中林地分布偏少，面积为33 651.6 km2，轮作耕地分布过多，面积为121 661.7 km2，与调查情况不符。混淆矩阵如表3所示，林地、其他耕地与轮作区耕地之间错分较多，有较大比例的林地和其他耕地被误分为小麦-玉米轮作区。

3.3 时序分类结果

基于时序数据使用识别算法提取结果精度较高，轮作区耕地面积为68 032 km2，其中东北部（即黄河流域）、中部靠太行山一侧和偏南部分布较多，分布情况基本与调查结果相符。二分类混淆矩阵如表4所示，只有少量轮作耕地和非轮作耕地发生错分，该算法能较为精确地对轮作区进行识别。

3.4 精度对比

基于时序特征的识别算法准确率、精准率、F1 Score和Kappa系数分别为0.987、0.970、0.975和0.956，分别比监督分类高出0.121、0.110、0.136和0.246。准确率、精准率、召回率和F1 Score越高，分类精度越高，这4个指标都是时序特征识别算法占优。Kappa系数表示混淆矩阵的不平衡度，监督分类结果的Kappa系数远低于时序识别算法。

4 结论

该研究基于GEE平台，分别利用识别算法与监督分类，对研究区小麦、玉米轮作区耕地进行提取，并进行精度对比分析，得到以下结论：

（1）该研究实现的提取算法可以快速有效地实现大区域小麦-玉米轮作耕地的提取，其二分类准确率、精准率、F1 Score和Kappa系数分别为0.987、0.970、0.975和0.956，分别比基于多光谱的分类高出0.121、0.110、0.136和0.246。

（2）某些特殊时期，林地、轮作耕地及其他耕地具有相似的光谱特征，传统基于多光谱数据的监督分类难以区分，导致其难以准确地提取轮作耕地。

（3）研究区2018年约有小麦-玉米轮作耕地68 032 km2，约占研究区总面积的30.7%，主要分布在河北东南部、河南东北部和山东西北部区域。

利用遥感数据的各类归一化指数的时序特征，能够比利用多光谱数据的监督分类更准确地区分具有类似光谱特征的地物，如轮作耕地与其他耕地，而在GEE平台的支持下，能快速获取并处理大区域、长时间的遥感数据，从而快速通过时序特征来区分不同地物。

参考文献

[1] 史舟，梁宗正，杨媛媛，等.农业遥感研究现状与展望[J].农业机械学报，2015，46（2）：247-260.

[2] 熊勤学，黄敬峰.利用NDVI指数时序特征监测秋收作物种植面积[J].农业工程学报，2009，25（1）：144-148.

[3] 王菲，李红华，苏洁.基于MODIS-NDVI的邯郸市冬小麦和夏玉米面积提取[J].现代农业科技，2018（24）：3-5.

[4] 路中，雷国平，马泉来，等.基于重构的Landsat 8时间序列数据和温度植被指数的区域旱情监测[J].水土保持研究，2018，25（5）：371-377，384.

[5] 苏伟，姜方方，朱德海，等.基于决策树和混合像元分解的玉米种植面积提取方法[J].农业机械学报，2015，46（9）：289-295，301.

[6] 翟世常，王学强，颉耀文.黑河流域中游盆地玉米作物遥感估产研究[J].安徽农业科学，2015，43（31）：350-353.

[7] 劉佳，王利民，杨福刚，等.基于GF-1卫星影像的中国冬小麦制图研究[J].中国农学通报，2019，35（9）：155-164.

[8] 刘沼辉，柳林，郭慧，等.基于GF1-NDVI时序影像对春小麦进行提取研究[J].北京测绘，2018，32（6）：643-646.

[9] 赵丽花，李卫国，杜培军.基于多时相HJ卫星的冬小麦面积提取[J].遥感信息，2011，26（2）：41-45，50.

[10] 谭深，吴炳方，张鑫.基于Google Earth Engine与多源遥感数据的海南水稻分类研究[J].地球信息科学学报，2019，21（6）：937-947.

[11] 侯学会，隋学艳，姚慧敏，等.中国北方麦区冬小麦物候期对气候变化的响应[J].麦类作物学报，2019，39（2）：202-209.

[12] FOGA S，SCARAMUZZA P L，GUO S，et al.Cloud detection algorithm comparison and validation for operational Landsat data products[J].Remote sensing of environment， 2017，194：379-390.

[13] 向海燕，罗红霞，刘光鹏，等.基于 Sentinel-1A 极化 SAR 数据与面向对象方法的山区地表覆被分类[J].自然资源学报，2017，32（12）：2136-2148.

[14] DIDAN K.MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 [EB/OL].[2020-11-05].https：//doi.org/10.5067/MODIS/MOD13Q1.006.