2019-2021年河南省越冬作物分布数据集

马晓岩，张兴旺，乔龙鑫，郑泽琳，潘力，彭洁，杨涵璐，夏浩铭,2,3,4*

1.河南大学地理与环境学院，河南开封 475004

2.河南大别山森林生态系统国家野外科学观测研究站，河南开封 475004

3.河南省地球系统观测与模拟重点实验室，河南开封 475004

4.黄河中下游数字地理技术教育部重点实验室（河南大学），河南开封 475004

引 言

越冬作物（如冬小麦、油菜籽、大蒜）是中国北方重要的粮食来源，种植越冬作物可以减少冬季闲置土地，提高粮食总产量[1-2]。精准的作物分布图可以为作物生长监测和产量预测提供基础数据。因此，有必要准确、快速绘制越冬作物种植面积图来帮助决策者和生产者制定合理的政策和风险管理策略[3-4]。传统的农业调查方法包括实地调查、逐级统计报告、问卷调查、访谈等[5]，该方法费时费力，且缺乏准确的作物空间分布信息，时效性差[6-8]。遥感分类法基于遥感影像快速准确地对作物进行分类，是一种非常有效的植被监测方法[9]，已广泛应用于区域、国家或全球作物制图[10-15]。

由于中低空间分辨率卫星（例如MODIS影像）在监测破碎化农田或高度异质性地区中的局限性[16-17]，Landsat（30 m）和Sentinel（10 m、20 m和60 m）等高空间分辨率卫星影像被越来越多地用于陆地植被监测。然而，这些高空间分辨率卫星通常时间分辨率较低。此外，影像质量容易受云或云阴影的影响，导致可用高质量影像数量不足，难以构建完整的植被时间序列曲线。因此，许多研究开始融合多源影像数据，进而显著提高影像的时间分辨率，并减少由异常值（由云、云阴影等其它因素所导致）带来的干扰[18-21]。本研究基于谷歌地球引擎(Google Earth Engine, GEE)，通过融合Landsat 7/8和Sentinel-2A/B 卫星影像，开发了一种新的基于物候学的算法来提取越冬作物。本数据集提供了大区域尺度上越冬作物年的空间分布和年际间的面积变化，有利于政府和决策者制定合理的政策和风险管理策略，也可为相关领域的科研人员提供数据参考。

1 数据采集和处理方法

本文的技术路线如图1所示，包含以下处理步骤：

图1 基于Landsat 7/8和Sentinel 2A/B时间序列数据识别和绘制越冬作物的工作流程Figure 1 The workflow for identifying and mapping winter crops fields by time series based on Landsat 7/8 and Sentinel 2A/B data

1.1 数据来源及其预处理

1.1.1 Landsat-7/8和Sentinel-2A/B数据及其预处理

Landsat 7/8分别于1999年4月和2013年2月发射，搭载了增强型专题制图仪(ETM+)和陆地成像仪(OLI)，空间分辨率为30 m，单颗卫星重访周期为15天[22]。Sentinel-2由两颗卫星组成，分别于2015年6月(Sentinel-2A)和2017年3月(Sentinel-2B)发射，携带空间分辨率为10 m、20 m和60 m的多光谱成像仪(MSI)，两颗卫星的重访周期为5天[23]。本研究所用遥感影像来源于GEE平台，GEE提供了大气校正后的地表反射率(SR)数据，与大气天顶反射率相比，地表反射率数据对陆地表面的变化更加敏感[24]。本文使用了研究期间河南省所有可用的Landsat 7/8和Sentinel-2A/B的SR数据作为初始影像集（表1）。

表1 2019-2021年不同物候阶段可用遥感影像数量Table 1 Number of available images at different phenological stages during 2019-2021

为了去除云和云阴影等对影像表面反射率的影响，我们使用GEE中的CFmask算法检测Landsat 7/8影像中质量较差的观测值。该方法是根据云的物理特性识别潜在的云像素，通过近红外波段的变暗效应检测出云阴影，再通过匹配云和云阴影生成劣质观测掩膜。对于Sentinel-2A/B影像数据，我们使用元数据中的质量波段来识别影像中的劣质观测。元数据将那些包含卷云和不透明云的观测识别为劣质观测，并将它们作为NODATA存储在图像中，以生成劣质观测掩膜。通过云掩膜去除云和云阴影得到高质量观测影像。

由于传感器在轨道、空间和光谱配置上的差异，会影响影像的物理测量值和辐射特性[25]，并且每个传感器的波长可能略有不同。因此，有必要对遥感反射率数据进行归一化以获得可比较的结果。本研究使用普通最小二乘回归系数将ETM+和MSI波段转换为OLI波段标准。需要处理的两个波段是ρRED和ρNIR。本研究计算的两个波段的回归系数分别为ETM+的 0.9175、0.9281和MSI的0.9103、0.9701[26-27]。

1.1.2 验证数据

本研究在2021年3月到6月，在河南省内进行了一系列实地采样，最终得到了180个越冬作物验证样本和60个其他作物验证样本。

1.1.3 土地覆盖数据

本文使用的土地覆盖数据为中国地表覆盖数据(CLCD, https://doi.org/10.5281/zenodo.4417810)。CLCD数据集空间分辨率为30 m，包含耕地等九个土地覆被类型，CLCD的总体准确率为79.31%。基于5131个第三方测试样本的进一步评估表明，CLCD的整体准确度优于MCD12Q1、ESACCI_LC、FROM_GLC和GlobaLand30[28]，CLCD得出的耕地面积与除西藏以外的大多数省区统计数据可以很好地吻合，特别是在华北平原，具有很高的一致性[29]。因此，本研究使用CLCD数据集中的Cropland像元对影像进行掩膜。

1.2 数据处理方法

1.2.1 初始植被指数时间序列构建

归一化差异植被指数（NDVI）时间序列曲线是识别农作物常用且有效的方法[30-31]。本文使用Landsat-7/8和Sentinel-2A/B中的近红外波段和红色波段计算并构建了NDVI时间序列数据集来表征植被的生长动态。其计算公式如下：

式中ρNIR为近红外波段反射率，ρRED为红色波段反射率。

1.2.2 高时空分辨率时间序列数据集构建

考虑到不同传感器之间的时间分辨率不同且拍摄区域重叠。为了获得等时间间隔的时间序列数据集，本研究计算了每个像元每10天内NDVI的最大值。由于部分像元受到云、雪或其他因素的影响被当作劣质观测去除，本研究基于空隙处前后的高质量观测值通过线性插值进行填充。此外，对于数据集中由云残留、数据传输错误或冰雪覆盖等引起的噪声[32]，本研究使用 Savitzky-Golay滤波器（S-G滤波器）对时间序列曲线进行平滑，进而降低噪声并构建高质量的NDVI时间序列曲线。

1.2.3 物候提取模型

本研究通过实地调查，从地表参考数据中选择了越冬作物样本点。基于这些样本点和高时空分辨率时间序列数据集构建了NDVI时间序列曲线（图2）。

图2 2020年9月1日至2021年9月1日越冬作物NDVI时间序列及5个物候指标Figure 2 NDVI time series of winter crops from Sep.1, 2020 to Sep.1, 2021, and five phenological indicators

图2显示了越冬作物NDVI值的时间序列变化情况，越冬作物在10月中下旬播种，此时NDVI值逐渐增加，冬季停止增长。2-3月返青后进入快速绿化期，NDVI值迅速增加，4-5月达到生长高峰期，此时NDVI值达到峰值。之后，6月至7月越冬作物进入成熟和收获期，NDVI值逐渐下降。当越冬作物的生长周期结束时，NDVI值再次跌至谷底。根据以上信息，本研究提取了用于分类的物候指标分别为生长季节始期(SOS)、生长高峰期(SDP)、生长季末期(EOS)、绿化速度(GUS)和生长季长度(GSL)。结合NDVI时间序列的峰谷和阈值法计算SOS、SDP和EOS。GUS为快速绿化期和高峰期之间NDVI的变化量与日期跨度的比值，GSL为EOS和SOS之间的差异。具体步骤如下：

（1）识别NDVI时间序列中的波峰和波谷。基于NDVI时间序列，算法会对像元进行迭代，当时间序列中某一时刻的值大于（小于）其上一个值和下一个值时，定义为波峰（波谷），每个像元上的波峰和波谷值以DOY（一年中的第几天）的形式存储。

（2）筛选有效波峰和波谷。然后将波峰的阈值设置为0.6以消除由杂草或收获不全面导致的额外波峰值[33]。

（3）确定SOS、SDP和EOS的日期[34-35]。在本研究中，通过计算NDVI变化的最大速率来确定越冬作物的SOS和EOS。公式如下：

式中：NDVImin和NDVImax——为前一年9月1日到该年6月30日NDVI的最小值和最大值。

从作物出苗来看，NDVIratio随着作物NDV I值的增加而增加。NDVIratio为0和1的时刻分别代表时间序列中的裸土期和SDP。在本研究中，SOS定义为NDVIratio首次增加到 0.2时的DOY，EOS定义为NDVIratio首次下降到0.2时的DOY[36]，SDP定义为出现波峰的DOY。通过以上步骤，本研究得到了像元的SOS、SDP和EOS。

（4）绿化速度(GUS)和生长季长度(GSL)的计算公式如下：

式中：NDVIpeak——生长周期内的 NDVI峰值; NDVIrevive——返青期开始时的第一个 NDVI值；DOYpeak和DOYrevive——NDVIpeak和NDVIrevive对应的日期。

在本研究中，返青期定义为2019年3月1日。GUS 定义为返青期NDVI的日均增长速度。GSL被定义为从SOS到EOS的时间跨度。

（5）根据地面实地调查越冬作物的样本点，得出了不同作物物候信息的差异，如图3所示，图中NDVI/D表示平均每天NDVI的变化量，以此来衡量GUS的变化速率。本研究发现越冬作物的SOS通常发生在10月底至11月底，SDP发生在3月下旬至5月初，EOS发生在5月下旬至7月初。在快速绿化期，NDVI日增长率超过0.001，完整生长周期的时间跨度超过200天。基于这些物候指标，本研究利用决策条件 260 ＜ SOS ＜ 330、80 ＜ SDP ＜ 140、130 ＜ EOS ＜ 190、GUS ＞ 0.001 和200 ＜ GSL ＜ 250 在像元尺度上进行分类，生成河南省越冬作物分布图。

1.2.4 精度评价方法

在高分辨率影像中对验证样本点定位，然后参考样本点所在耕地边界，进行矢量化，得到验证多边形，多边形的大小取决于验证点所在的耕地的面积。最终得到用于验证越冬作物的多边形为180个（3565 像元），其他作物的多边形为60 个（1547 像元）。然后，用混淆矩阵对本分类模型进行精度评价。评价指标包括用户精度(User Accuracy, UA, %)、生产者精度(Producer Accuracy, PA, %)、总精度(Overall Accuracy, OA, %)和Kappa系数。

图3 研究区越冬作物与其他作物物候指标统计Figure 3 Statistics of phenological indexes of winter crops and other crops in the study area

2 数据样本描述

2.1 数据组成

本数据集包含2019年到2021年河南省越冬作物分布数据。数据集被保存在一个名为“2019年到 2021年河南省越冬作物分布数据.7z”的压缩文件内，压缩文件大小为 34.4 MB。该数据集包括2019年、2020年、2021年的河南省越冬作物分布图共3张，为3个*.tif文件，文件以对应年份进行命名，影像中DN值等于1为越冬作物，DN值等于0为其他。

2.2 数据样本

越冬作物主要分布在河南东部和北部的平原地区（图4），占河南省越冬作物面积的80%以上。由于河南省西部主要为山地，地形起伏大且地面坡度较大，不适宜种植越冬作物。河南省南部（如信阳市）降水充沛，水田分布密集，所以多种植水稻等作物，而水稻常在9月到10月收割，收割后稻田较湿，多为单季，不适宜种植越冬作物。

图4 2019（a）、2020（b）、2021（c）研究区越冬作物分布图Figure 4 Distribution maps of winter crops in the study area during 2019 (a), 2020 (b), and 2021 (c)

3 数据质量控制和评估

数据的质量控制和评估是遥感影像信息提取中尤为重要的一环，图5给出了本文数据质量控制和评估的技术路线。数据质量控制包含影像质量控制和样本质量控制。数据质量评估主要是在执行完算法后，采用验证样本集对越冬作物分布结果的分类精度进行评估。

图5 数据质量控制和评估的工作流程Figure 5 The workflow for data quality control and evaluation

本研究将实地采样调查得到的样本点与本文物候提取模型得到的结果进行比较，计算混淆矩阵如表2所示，其中用户精度、生产者精度和总体精度分别为98.00%，98.36%和97.43%，Kappa系数为0.94，说明分类结果与实地采样调查之间的一致性较强。

表2 基于2021年越冬作物和其他作物地面参考数据的精度评估Table 2 Accuracy assessment based on the ground reference data for winter crops and others in 2021

利用《河南省统计年鉴》[37-39]中越冬作物播种面积与本数据集进行了对比。2021年越冬作物分类结果面积与统计年鉴面积差异为0.88%；2020年为3.80%；2019年为11.54%（表3）。其中2019年的面积差异较大，可能的原因是2018年冬季缺少Sentinel-2A/B影像，导致其可使用的高质量影像数量较少。此外，研究区处于中国东部的季风区，受云的影响较大，很难保证每个像元每10 天有一幅高质量卫星影像，本研究采用该像元相邻时间的像元值进行线性插值算法来填补缺失值，但是在拐点处缺失影像时，插值低估作物生长的真实位置，存在一定的不确定性。同时，由于本研究构建的30 m空间分辨率数据集分辨率有限，对于研究区内地块面积小于900 m2的农田识别较为困难，也可能会导致部分小块农田缺失。由于传统的农业普查难以保证精度，统计数据往往会被高估，因此本数据集满足大范围制图精度要求。

表3 本文越冬作物分类结果与统计数据的面积差异Table 3 Area differences between classification results in the study and statistical data for winter crops

4 数据使用方法和建议

本数据集是公开发布的河南省2019年到2021年越冬作物分布专题数据，且具有较高的空间分辨率和分类精度，可为河南省越冬作物生长监测和产量预测提供基础数据，帮助决策者和生产者制定合理的政策和风险管理策略。同时，该数据集也具有较高的科学意义和实用价值。高分1/2/6卫星数据已成功用于作物制图，未来在数据完全放开的情况下，在后续的研究中，我们会尝试结合高分卫星数据进行更高空间分辨率的遥感制图。