基于Google Earth Engine的淮河流域越冬作物种植面积制图

潘 力，夏浩铭,2，王瑞萌，牛文辉，田海峰，秦耀辰,2

基于Google Earth Engine的淮河流域越冬作物种植面积制图

潘 力1，夏浩铭1,2※，王瑞萌1，牛文辉1，田海峰1，秦耀辰1,2

（1. 河南大学地理与环境学院/河南省地球系统观测与模拟重点实验室/黄河文明与可持续发展研究中心暨黄河文明省部共建协同创新中心，开封 475001；2. 黄河中下游数字地理技术教育部重点实验室（河南大学），开封 475001）

越冬作物是中国重要的作物类型，其面积的变化不仅对中国粮食产量和经济产生直接的影响，还潜在地影响中国的粮食安全，因此有必要准确绘制越冬作物种植面积图来为决策制定者提供科学参考。该研究以淮河流域为例，基于Google Earth Engine云平台，融合时间序列Landsat-7/8和Sentinel-2A/B卫星影像，采用S-G（Savitzky-Golay）滤波器重构作物时间序列的归一化差异植被指数（Normalized Difference Vegetation Index，NDVI），根据不同植被类型物候期的差异，选取越冬作物生长旺盛期NDVI最大值、越冬作物播种期和收获期中相应的NDVI最小值和中位数，在像元尺度上构建越冬作物提取算法，绘制淮河流域越冬作物的种植面积。研究结果表明，所构建算法能够精确提取淮河流域越冬作物的种植面积，总体精度为95.8%，Kappa系数为0.912，该研究可为作物面积的提取和监测提供方法参考。

遥感；作物；制图；种植面积；Google Earth Engine；淮河流域

0 引 言

越冬作物为秋季播种，幼苗经过冬季，并在第二年春季或夏季收割的农作物，例如冬小麦、冬大蒜和冬油菜等，是中国北部重要的粮食作物，及时准确地获取越冬作物种植面积对于预测粮食产量、价格和国家的粮食安全至关重要[1-3]，并可用于提高农业生产力和韧性[4]。精准的区域作物种植分布图为作物生长监测和产量预测提供了基础数据，有助于帮助决策者及生产者制定合理的政策和风险管理策略[5-6]。

通常作物种植面积的获取方法包括农业调查法和遥感分类法。传统的农业调查法包括实地问卷调查和访谈等，逐层统计上报，此类方法既耗时又费力，缺少精确的作物空间分布信息，且时效性较差。遥感分类法能够基于一幅或多幅影像精确地进行作物分类，是一种非常有效的植被监测方法，已广泛应用于区域、国家和全球作物制图等[7-9]。随着遥感数据的激增，免费共享和遥感大数据处理方法的发展，利用遥感技术进行作物播种面积制图逐渐成为作物种植数据统计的重要途径。

作物种植面积遥感制图常用的方法是利用卫星数据和分类算法在地表反射特征和作物生长特征之间建立关系[10]。由于中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）具有高时间分辨率的特点，且具有很好的时间序列特征，目前已开发了多种基于MODIS数据的算法来监测中等空间分辨率下越冬作物的空间分布、估产、时空变化以及物候特征等，如周亮等[11]使用MODIS数据构建了基于卷积神经网络的冬小麦估产模型，张锦水等[12]基于MODIS数据构建了2017—2019年冬小麦时间序列数据对河北省冬小麦空间分布进行识别，Ren等[13]基于MODIS数据对黄淮海平原2001—2016年冬小麦进行估产，黄健熙等[14]利用MODIS的叶面积指数产品提取了冬小麦返青期、拔节期、抽穗期和开花期。尽管MODIS数据在进行作物动态监测上得到了很好的应用，但因受到空间分辨率的限制，难以满足破碎化、异质性小块农田的监测需求。随着免费的中高空间分辨率卫星的迅速发展，Landsat影像（30 m）逐渐应用于遥感制图，如Phalke等[15]利用160 000幅Landsat影像绘制了30 m空间分辨率的欧洲、中东、俄罗斯和中亚大部分地区的农田地图，Wang等[16]利用Landsat影像绘制了1999—2018年美国中西部30 m空间分辨率的玉米和大豆地图。Landsat影像解决了MODIS影像空间分辨率不足的问题，但由于Landsat卫星的时间分辨率较低，且影像经常受到云层的干扰，导致很难使用Landsat影像构建作物完整平滑的时间序列曲线。为了减少因Landsat低时间分辨率而导致的时间序列曲线值缺失的影响，许多研究开始结合使用多个传感器的融合影像来进行研究[17-18]，如Li等[19]结合MODIS和Landsat数据绘制了安徽宣城作物种类和轮作的空间分布图，Liu等[8]结合Landsat-7/8和Sentinel-2数据绘制了中国30 m空间分辨率的种植强度图，Lenco等[20]结合Sentinel-1和Sentinel-2在10 m空间分辨率下对土地覆被进行分类。多个传感器数据源的整合可以提高影像的时间分辨率，减少云覆盖等因素的干扰，进而构建作物完整生命周期曲线，利用不同时期的光谱差异来区分作物的类型。

综上，本研究基于Google Earth Engine云平台，通过融合Landsat-7/8和Sentinel-2A/B数据集，首先构建了30 m高时空分辨率时间序列影像数据集，并计算归一化植被指数（Normalized Difference Vegetation Index，NDVI）；其次采用S-G（Savitzky-Golay）滤波器对NDVI时间序列进行平滑，基于平滑后的NDVI时间序列确定越冬作物生长旺盛期、越冬作物播种期和收获期，并分别计算生长旺盛期NDVI最大值、播种期和收获期相应的NDVI最小值和中位值；最后通过构建分类模型，在像元尺度上绘制2018年淮河流域越冬作物种植面积图。本研究可为作物面积的提取和监测提供科学依据和方法参考，为作物管理和产量预测提供数据支撑。

1 材料与方法

1.1 研究区概况

淮河流域位于中国东部南北气候过渡带（111°55′E～121°20′E，30°55′N～36°20′N），横跨河南、安徽、江苏、山东、湖北5 省（图1）。淮河流域西部、西南部及东北部为山区和丘陵区，其余为平原。流域面积为27万km2，山区约占总面积的1/3，平原约占总面积的2/3。淮河流域是中国重要的农业种植区域，广阔的平原为作物提供了适宜的生长环境，耕地面积占整个流域面积的68%，占全国耕地面积的12%，粮食产量约占全国总产量的17%[21]。淮河流域种植模式多为单季种植和两季种植，其复种比例达到58.4%[22]。单季种植耕地多为单季水稻、单季玉米等，两季种植耕地多为冬小麦-玉米、冬小麦-花生等。近年来，淮河流域耕地面积呈递减趋势，越冬作物作为淮河流域主要作物种植类型，为保障国家粮食产量和粮食安全，探明淮河流域现有越冬作物种植状况对越冬作物未来生产潜力的评估和提升十分重要。

1.2 数据来源

1.2.1 Landsat-7/8和Sentinel-2卫星数据及处理

Landsat-7/8卫星携带增强型专题制图仪（Enhanced Thematic Mapper，ETM+）和陆地成像仪（Operational Land Image ，OLI），空间分辨率为30 m，重访周期为15 d[23]。Landsat-7 ETM+和Landsat-8 OLI影像分别对应“LANDSAT / LE07 / C01 / T1_TOA”和“LANDSAT / LC08 / C01 / T1_TOA”。Sentinel-2由2颗卫星组成，携带多光谱成像仪（Multi-Spectral Instrument，MSI)，空间分辨率为10、20和60 m，单颗星重访周期为10 d[24]。Sentinel-2 A/B对应“COPERNICUS / S2”。本研究以越冬作物的一个完整物候期为例，选取时间范围为2017年9月1日至2018年6月30日，合计4 982景影像。

基于Google Earth Engine云平台，利用JavaScript编程获取2017年9月1日至2018年6月30日淮河流域逐像元Landsat-7/8和Sentinel-2A/B卫星影像的总观测影像幅数如图2a和图2b所示。同时，采用Google Earth Engine云平台的CFMask算法来识别云、云阴影、卷云和冰/雪覆盖，去除所有质量差的观察结果[25]，获得研究区逐像元高质量观测影像幅数如图2c和图2d所示。结果表明，在研究期间，研究区单个像元上观测影像幅数最少为60幅，最高为534幅，其中高质量观测影像幅数最少为12幅，最高为311幅。

1.2.2 样本数据

为确保足够数量和质量的样本点数据用于分类及精度验证，本研究基于中国气象网（http://data.cma.cn/）提供的观测数据资料、Google Earth高分辨率影像和野外调查数据（2018年3月15日—2018年7月1日）获取了共621个具有代表性的、典型的像元作为样本点（图 1），分别为越冬作物样本点250个、春播作物样本点87个、落叶森林样本点76个、常绿森林样本点67个、建筑用地样本点120个和水体样本点21个。本研究按照3∶1的比例将样本划分为训练样本（463个）和验证样本（158个）。

Google Earth影像（空间分辨率为1 m）能清晰地显示耕地、河流、道路、房屋村庄等地表特征，可作为本研究制图结果的直接验证数据。同时，本研究还使用2019年河南省统计年鉴[26]来间接验证本研究的算法。由于淮河流域横跨河南省、山东省、安徽省、江苏省和湖北省5省，每个省统计年鉴的统计指标并不一致，仅有河南省统计年鉴细化至县市级别，因此本研究选取淮河流域所覆盖的河南省62个县市的统计数据用于验证，由于其中9个县市位于淮河流域的边界，不完全属于淮河流域，故剔除这9个县市，最后保留河南省53个县市作为验证区域。

1.2.3 坡度数据

本研究选取数字高程模型（Digital Elevation Model，DEM）数据，计算研究区的坡度（°），DEM数据可在Google Earth Engine云平台中免费获取，空间分辨率为30 m。

1.3 研究方法

本研究的主要思路如图3所示，主要包括4个方面：1）数据集构建，包括构建30 m影像数据集、时间序列NDVI计算、插值和重构；2）基于平滑后的时间序列NDVI数据集，结合越冬作物物候期构建越冬作物提取算法；3）逐像元绘制2018年淮河流域越冬作物分布图；4）利用Google Earth高分辨率影像和统计数据验证越冬作物提取算法精度。

1.3.1 NDVI时间序列数据集构建及处理

基于NDVI时间序列曲线识别作物是常见且有效的方法，NDVI的计算如（1）所示。

由于不同传感器所获得影像的时间分辨率不同，扫描区域存在重叠，为了获取相等时间间隔的时间序列影像集，本研究以10 d内所有高质量影像的NDVI最大值作为该10 d的NDVI值[18,27]。如果影像受云、雪或其他因素的影响，在某些地区无法获得10 d内高质量影像的NDVI值[8,17,22]，本研究则采用相邻像元线性插值算法来填补缺失值，以构建完整的NDVI时间序列数据集[28]。

即使经过严格的预处理，NDVI时间序列数据集中仍残留由云、大气和二向性反射引起的噪声[29]。因此，本研究利用S-G滤波器对NDVI时间序列曲线进行平滑，以消除时间序列中的噪声[30]。

S-G滤波器于1964年被首次提出[31]，是一种在给定滑动窗口宽度内基于局部多项式最小二乘拟合的滤波方法，该滤波器最大的特点为在滤除噪声的同时可以保持曲线的变化趋势。使用S-G滤波器对NDVI时间序列进行平滑如式（2）所示。

式中Y为平滑后的NDVI值；Y为滑动窗口内平滑前Y前后第个点的NDVI值；C为滤波系数，即第期Y在滑动窗口中的权重；为滑动窗口的数据个数2+1，其中为滑动窗口宽度。

1.3.2 确定时间窗口

为确定提取越冬作物的最佳时间窗口，本研究选取具有代表性的越冬作物、春播作物、落叶森林、常绿森林、水体和建筑用地的样本点创建NDVI时间序列曲线。本研究设定滑动窗口为90 d，使用S-G滤波对NDVI时间序列曲线逐像元进行平滑。由于水体和建筑用地的NDVI曲线在整个研究期（2017年9月1日至2018年6月30日）保持稳定，水体的NDVI<0，建筑用地的NDVI<0.4，与植被曲线差异较大，因此主要对越冬作物和其他植被的NDVI曲线差异进行研究。

由平滑后越冬作物、春播作物、落叶森林和常绿森林的NDVI曲线可知（图4），越冬作物在10月中下旬开始播种，此时NDVI值最小（NDVI<0.3），11月份初出苗后，NDVI值逐步增大，1月份停止生长进入越冬期后，NDVI值保持相对不变，2月份越冬作物开始返青，3月份进入生长旺盛期，NDVI值恢复增长，并于成熟期达到峰值（NDVI>0.6），此后NDVI值逐步下降，到5月底至6月初完成收割，NDVI变为最小值（NDVI<0.3）；春播作物在2月中下旬播种前，9月份至第二年2月中下旬NDVI值较低（NDVI<0.3），播种后NDVI值逐渐增大，6月份下旬达到最大值（NDVI>0.6）；落叶森林在10月份之前NDVI值较高（NDVI>0.6），10月份后NDVI值逐渐降低，并于12月份达到NDVI最低值（NDVI<0.3），3月份上旬返青后NDVI逐渐增大（NDVI>0.6）；常绿森林在整个时期均保持较高NDVI值（NDVI>0.6）。

由以上分析得出，越冬作物的NDVI曲线与其他植被的NDVI曲线之间存在明显差异，即越冬作物播种后生长旺盛，春播作物的耕地撂荒、落叶森林枯萎（NDVI值较低）。越冬作物在播种和收获时，春播作物、落叶森林和常绿森林则处于生长旺盛期（NDVI值较高）。因此，本研究选择越冬作物生长旺盛期、播种期和收获期来区分越冬作物和其他植被。根据越冬作物生命周期，确立越冬作物生长旺盛期为2018年03月20日至2018年04月20日，越冬作物播种期为2017年10月11日至2017年11月10日，越冬作物收获期为2018年05月20日至2018年06月30日。

1.3.3 分类参数

对于越冬作物生长旺盛期，在Google Earth Engine云平台使用Max函数逐像元求NDVI最大值（NDVImax），获得NDVImax图像。对于越冬作物播种期和收获期，在Google Earth Engine云平台使用Median函数和Min函数逐像元分别求解2个时期的NDVI中位值（NDVImedian）和NDVI最小值（NDVImin），获得NDVImedian和NDVImin图像。本研究使用NDVImedian是由于在部分云阴影影响区，常绿森林的NDVI值较低，会与越冬作物混淆，故使用NDVImedian提高越冬作物分类阈值，避免这类型常绿森林错误地分成越冬作物，提高了分类精度。至此，获得了NDVImax、NDVImedian和NDVImin3个分类参数。

1.3.4 精度验证方法

本研究基于158个验证样本，对本分类模型进行精度评价。首先，以每个验证样本点为中心，生成100 m×100 m的矩形地面样方，共获取样方158个（包括越冬作物65个，其他地类93个）；其次，获取地面样方空间分辨率为1 m的Google Earth卫星影像，对每个地面样方的影像进行目视解译以获得地表参考数据；最后，将地表参考数据与本研究分类结果构建混淆矩阵，来评价本研究分类算法的精度。评价指标包括用户精度（User Accuracy，UA，%）、生产者精度（Producer Accuracy，PA，%）、总精度（Overall Accuracy，OA，%）和Kappa系数，其计算如式（3）～式（6）所示

式中n为混淆矩阵中第行列的值，n为混淆矩阵中第行的和，为混淆矩阵中第列的和，为验证样本总数，为混淆矩阵行列数。

2 结果与分析

2.1 分类模型构建结果

为了构建分类模型，本研究基于463个训练样本（186个越冬作物样本点、65个春播作物样本点、57个落叶森林样本点、50个常绿森林样本点、90个建筑用地样本点和15个水体样本点），绘制了NDVImax、NDVImedian和NDVImin3个分类参数的箱型图如图5所示，通过分析越冬作物与其他地类的特征参数差异来构建决策树分类模型。具体操作如下：

1）由于选取的样本数无法完全代表研究区内所有越冬作物的物候特征，为了确保能提取出所有越冬作物，本研究将阈值的范围向外扩大了10%。由图5a可知，越冬作物的NDVImax的范围在(0.36,0.82)之间，扩大后的范围为(0.33,0.91)，而水体和建筑用地的NDVImax上限小于0.33，因此将越冬作物NDVImax的阈值设置为0.33可将水体和建筑用地的像元排除。

2）去除了水体和建筑用地像元的干扰后，下一步将区分越冬作物和其他植被。由图5b可知，越冬作物NDVImedian的范围为(0.19,0.45)，扩大后的范围为(0.17,0.50)，有75%其他植被的像元值在0.50之上。因此将NDVImedian的阈值设置为0.50去除75%的其他植被像元。由图5a可知，落叶森林和春播作物的NDVImax范围分别在[0.24,0.44]和[0.18,0.43]之间，扩大范围后这2类植被的NDVImax最高不超过0.48，超过0.48的落叶森林和春播作物像元将被去除。由于越冬作物NDVImax的范围在(0.36,0.82)之间，因此分为NDVImax≤0.48和NDVImax>0.48两种情况来处理：1）当研究区内像元的NDVImax≤0.48时，可以通过设置NDVImin阈值来去除常绿森林像元。由图5c可知，常绿森林的NDVImin范围在0.44以上，此时越冬作物像元的NDVImin范围在(-0.11,0.16)之间，扩大后为(-0.12,0.17)，因此设置判别条件-0.12< NDVImin<0.17来去除该条件下常绿森林像元；2）当研究区内像元的NDVImax>0.48时，设置判别条件-0.12 < NDVImin<0.5NDVImax去除剩余的常绿森林像元。

3）本研究还选取高程模型中的坡度数据来提升精确度，越冬作物设置的坡度阈值为10°。

基于上述分析，构建越冬作物决策树分类模型如图6所示。使用决策树模型第一层“坡度小于10°”判定条件去除非耕地像元，使用决策树模型第二层“NDVImax>0.33”判定条件去除水体和建筑用地像元，使用决策树模型第三层“NDVImedian<0.5”判定条件去除落叶森林和春播作物像元。在决策树模型的最后一层，若像元满足“NDVImax≤0.48”和“-0.12< NDVImin<0.17”条件，或满足“NDVImax>0.48”和“-0.12< NDVImin< 0.5NDVImax”条件，则该像元被判定为越冬作物，否则为其他地类。

a. NDVI最大值 a. NDVImaxb. NDVI中位值 b. NDVImedianc. NDVI最小值 c. NDVImin

注：NDVImax为越冬作物生长旺盛期（2018年3月20日至2018年4月20日）NDVI最大值；NDVImedian和NDVImin分别为越冬作物播种期（2017年10月11日至2017年11月10日）和收获期（2018年5月20日至2018年6月30日）中NDVI中位值和最小值。

Note: NDVImaxis the maximum value of NDVI during the period of peak growth stage (from 2018-03-20 to 2018-04-20); NDVImedianand NDVIminare the median and minimum values of NDVI during the period of sowing stage (from 2017-10-11 to 2017-11-10) and harvest stage (from 2018-05-20 to 2018-06-30),respectively.

图5 各地类样本的分类参数箱型图

Fig.5 Classification parameter box plots of samples of various land types

2.2 越冬作物分类结果

淮河流域越冬作物的总面积为8.762×106hm2，其中淮河流域内山东省、河南省、安徽省和江苏省的作物种植面积分别为1.299×106、2.895×106、2.431×106和2.126×106hm2（图7）。在空间分布方面，32°～35°N是越冬作物的主要分布区域，占淮河流域越冬作物总面积的89%；113°～120°E是越冬作物的主要分布区域，占淮河流域越冬作物总面积的97%；淮河流域东北部、西南部和西部越冬作物分布较少。调查发现东北部和西部多为山脉（图1），地势高，土质较差，且农田多为梯田，面积较小、收割较困难、不易大规模种植作物，大多数农田只在夏季种植一季作物（例如玉米、生姜、芋头、番薯）后便撂荒，冬季不再种植。研究区西南地区（例如信阳市）气温高降水足，水塘密集，灌溉充分，多种植水稻。水稻在9月末至10月初收割，水稻收割期通常雨水较多，导致稻田较湿，影响了越冬作物的播种，故大部分区域种植一季，仅有地势较高区域种植两季，淮河流域内部空白处主要为城市区、果园区（如砀山梨园区）、湖泊和河流（洪泽湖和微山湖）等。

2.3 精度验证结果

基于1.3.4节获取的地表参考数据，包括528个越冬作物像元（47.521 hm2）和819个其他像元（73.714 hm2），计算混淆矩阵。混淆矩阵的越冬作物分类结果精度验证如表1所示，本研究构建的模型对研究区越冬作物分类的用户精度为92.6%，生产者精度为97.0%，总精度95.8%，分类精度较高；Kappa系数为0.912，说明分类结果与地表参考数据二者间的一致性较强。

将基于目视解译后的地表参考数据与本研究模型的分类结果进行直观对比如图8所示，所选中的Google Earth影像区域是包含典型地物的矩形区域（图 8a），该区域内除越冬作物外，还包括房屋村庄、鱼塘、河流和森林等非越冬作物。由图8b和图8c对比可知，本研究模型的分类结果与目视解译后的地表参考数据相似度较高，表明本研究构建的模型以较高的精度提取出了研究区的越冬作物。分类误差（包括错误分类和遗漏分类）如图8d所示，越冬作物错误分类现象主要集中在耕地路网上，这是由于农村道路宽度一般在30 m以内，本研究构建的30 m空间分辨率数据集对这类线状地物识别能力较差；越冬作物遗漏分类则主要集中在村落建筑物周围，村落边界的不规则导致镶嵌在建筑物周边的小型越冬作物种植区域无法正确识别。因此，本研究构建的越冬作物提取模型能够满足大区域分类精度的需求。

2019年河南省统计年鉴[34]中县级单元（即本研究选出的河南省53个县市）越冬作物种植面积与本研究模型提取的越冬作物种植面积之间的相关性分析如图9所示。由图9可知，县级年鉴统计的种植面积与本研究提取的种植面积具有显著的线性关系（2=0.86，<0.01），表明本研究模型的分类结果与统计数据有较高的一致性。其中县级年鉴统计数据中验证区域的越冬作物总种植面积为3.030×106hm2，本研究提取的越冬作物在验证区域的总种植面积为2.804×106hm2，精度为92.6%。

表1 越冬作物分类结果精度验证

3 结 论

本研究基于Google Earth Engine云平台，融合Landsat-7/8和Sentinel-2A/B影像获取高时空分辨率影像集，结合越冬作物独特的物候期，构建了一种简单、有效和高精度的越冬作物识别模型，在像元尺度上绘制淮河流域越冬作物种植面积图，主要的结论如下：

1）融合多个传感器影像能够获取高时空分辨率数据集，克服当前大面积作物制图所用数据源时空分辨率的缺陷（例如MODIS空间分辨率低、Landsat时间分辨率低等）。本研究结合了所有可使用的高质量Landsat-7/8和Sentinel-2 A/B卫星影像（合计4 982景影像），极大地增加了每个像元上可用的高质量影像观测数，更好地捕获了越冬作物的物候。

2）通过10 d合成、插值和平滑算法，重构了时间序列归一化植被指数（Normalized Difference Vegetation Index，NDVI）数据集，揭示了越冬作物和其他植被的NDVI曲线差异，确定了越冬作物关键物候期：生长旺盛期（2018年3月20日至4月20日）和播种期（2017年10月11日至11月10日）和收获期（2018年5月20日至6月30日），成功构建了融合越冬作物独特物候特征的分类模型。

3）2018年淮河流域越冬作物种植面积为8.762×106hm2，其精度满足大区域作物制图的精度要求。

淮河流域处于中国东部季风区，受云的影响，很难保证每个像元10 d有一幅高质量卫星影像，本研究采用相邻像元线性插值算法来填补缺失值，但是在拐点处缺失影像时，插值并不能完全反映作物生长的真实位置，存在一定的不确定性。Sentinel-1A/B影像不受天气的影像，能够估算作物的物候期，未来可结合Sentinel-1 A/B卫星的C波段合成孔径雷达（Synthetic Aperture Radar，SAR）数据进行制图。本研究构建的30 m空间分辨率数据集难以提取研究区地块面积小于0.09 hm2的农田，未来可基于高分1/2/3/5/6等卫星获取10 m或更高空间分辨率的作物分布图。此外，云平台的高性能计算和并行处理能力可用于处理大尺度的地理空间数据集，是未来遥感大数据处理研究的趋势。

致谢：本研究得到了河南大别山森林生态系统国家野外科学观测研究站和国家地球系统科学数据中心-黄河中下游分中心（http://henu.geodata.cn）的数据支持，在此表示谢意！

[1] Franch B, Vermote E F, Becker-Reshef I, et al. Improving the timeliness of winter wheat production forecast in the United States of America, Ukraine and China using MODIS data and NCAR growing degree day information[J]. Remote Sensing of Environment, 2015, 161: 131-148.

[2] Shao Y, Campbell J B, Taff G N, et al. An analysis of cropland mask choice and ancillary data for annual corn yield forecasting using MODIS data[J]. International Journal of Applied Earth Observation and Geoinformation, 2015, 38: 78-87.

[3] 王利民，刘佳，杨福刚，等. 基于 GF-1 卫星遥感的冬小麦面积早期识别[J]. 农业工程学报，2015，31(11)：194-201.

Wang Limin, Liu Jia, Yang Fugang, et al. Early recognition of winter wheat area based on GF-1 satellite[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(11): 194-201. (in Chinese with English abstract)

[4] Hunt M L, Blackburn G A, Carrasco L, et al. High resolutions wheat yield mapping using Sentinel-2[J/OL]. Remote Sensing of Environment, 2019, 233, [2019-09-07], https: //doi. org/10. 1016/j. rse. 2019. 111410.

[5] Jiao X F, Kovacs J, Shang J L, et al. Object-oriented crop mapping and monitoring using multi-temporal polarimetric RADARSAT-2 data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 96: 38-46.

[6] Jin Z N, Azzari G, Lobell D. Improving the accuracy of satellite-based high-resolution yield estimation: A test of multiple scalable approaches[J]. Agricultural and Forest Meteorology, 2017, 247: 207-220.

[7] 解毅，张永清，荀兰，等. 基于多源遥感数据融合和LSTM算法的作物分类研究[J]. 农业工程学报，2019，35(15)：129-137.

Jie Yi, Zhang Yongqing, Xun Lan, et al. Crop classification based on multi-source remote sensing data fusion and LSTM algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 129-137. (in Chinese with English abstract)

[8] Liu L, Xiao X M, Qin Y W, et al. Mapping cropping intensity in China using time series Landsat and Sentinel-2 images and Google Earth Engine[J/OL]. Remote Sensing of Environment, 2020, 239, [2019-12-20], https: //doi. org/10. 1016/j. rse. 2019. 111624.

[9] Liu C, Zhang Q, Tao S Q, et al. A new framework to map fine resolution cropping intensity across the globe: Algorithm, validation, and implication[J/OL]. Remote Sensing of Environment. 2020, 251, [2020-09-08], https: //doi. org/10. 1016/j. rse. 2020. 112095.

[10] Tian H F, Huang N, Niu Z, et al. Mapping winter crops in China with multi-source satellite imagery and phenology-based algorithm[J]. Remote Sensing, 2019, 11(7): 820-843.

[11] 周亮，任建强，吴尚蓉，等. 基于卷积神经网络的中国北方冬小麦遥感估产[J]. 农业工程学报，2019，35(15)：119-128.

Zhou Liang, Ren Jianqiang, Wu Shangrong, et al. Remote sensing estimation on yield of winter wheat in North China based on convolutional neural network[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(15): 119-128. (in Chinese with English abstract)

[12] 张锦水，赵光政，洪友堂，等. 基于像元物候曲线匹配的生长季内河北省冬小麦空间分布识别[J]. 农业工程学报， 2020，36(23)：193-200.

Zhang Jinshui, Zhao Guangzheng, Hong Youtang, et al. Spatial extraction of winter wheat in Hebei in growing season using pixel-wise phenological curve[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(23): 193-200. (in Chinese with English abstract)

[13] Ren S J, Guo B, Wu X, et al. Winter wheat planted area monitoring and yield modeling using MODIS data in the Huang-Huai-Hai Plain, China[J/OL]. Computers and Electronics in Agriculture, 2021, 182, [2021-02-07], https:// doi. org/10. 1016/j. compag. 2021. 106049.

[14] 黄健熙，赵剑桥，汪雪淼，等. 基于遥感和积温的冬小麦生育期提取方法[J]. 农业机械学报， 2019，50(2)：176-183.

Huang Jianxi, Zhao Jianqiao, Wang Xuemiao, et al. Extraction method of growth stages of winter wheat based on accumulated temperature and remote sensing data[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(2): 176-183. (in Chinese with English abstract)

[15] Phalke A, Özdoğan M, Thenkabail P, et al. Mapping croplands of Europe, Middle East, Russia, and Central Asia using Landsat, Random Forest, and Google Earth Engine[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2020, 167: 104-122.

[16] Wang S, Tommaso S D, Deines J M, et al. Mapping twenty years of corn and soybean across the US Midwest using the Landsat archive[J/OL]. Scientific Data, 2020, 7, [2020-08-13], https: //doi. org/10. 1038/s41597-020-00646-4.

[17] Pan L, Xia H M, Zhao X Y, et al. Mapping winter crops using a phenology algorithm, time-series Sentinel-2 and Landsat-7/8 images, and Google Earth Engine[J/OL]. Remote Sensing, 2021, 13, [2021-06-21], https: //doi. org/10. 3390/rs13132510.

[18] Griffiths P, Nendel C, Hostert P. Intra-annual reflectance composites from Sentinel-2 and Landsat for national-scale crop and land cover mapping[J]. Remote Sensing of Environment, 2019, 220: 135-151.

[19] Li R Y, Xu M Q, Chen Z Y, et al. Phenology-based classification of crop species and rotation types using fused MODIS and Landsat data: The comparison of a random-forest-based model and a decision-rule-based model[J/OL]. Soil and Tillage Research. 2021, 206, [2020-10-07], https: //doi. org/10. 1016/j. still. 2020. 104838.

[20] Lenco D, Interdonato R, Gaetano R, et al. Combining Sentinel-1 and Sentinel-2 Satellite Image Time Series for land cover mapping via a multi-source deep learning architecture[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 158: 11-22.

[21] 王情，刘雪华，岳天祥. 淮河流域粮食生产潜力空间格局研究[J]. 生态经济，2014，30(7)：24-27.

Wang Qing, Liu Xuehua, Yue Tianxiang. Study on the spatial patterns of food production potential in the Huaihe River Basin[J]. Ecological Economy, 2014, 30(7): 24-27. (in Chinese with English abstract)

[22] Pan L, Xia H M, Yang J, et al. Mapping cropping intensity in Huaihe Basin using phenology algorithm, all Sentinel-2 and Landsat images in Google Earth Engine[J/OL]. International Journal of Applied Earth Observation and Geoinformation, 2021, 102, [2021-05-20], https: //doi. org/10. 1016/j. jag. 2021. 102376.

[23] Wulder M, Masek J, Cohen W, et al. Opening the archive: How free data has enabled the science and monitoring promise of Landsat[J]. Remote Sensing of Environment, 2012, 122: 2-10.

[24] Torres R, Snoeij P, Geudtner D, et al. GMES Sentinel-1 mission[J]. Remote Sensing of Environment, 2012, 120: 9-24.

[25] Foga S, Scaramuzza P, Guo S, et al. Cloud detection algorithm comparison and validation for operational Landsat data products[J]. Remote Sensing of Environment, 2017, 194: 379-390.

[26] 王世炎，夏雨春，等. 河南统计年鉴2019[M]. 北京：中国统计出版社，2019.

[27] Holben B. Characteristics of maximum-value composite images from temporal AVHRR data[J]. International Journal of Remote Sensing, 1986, 7(11): 1417-1434.

[28] Chen Y L, Song X D, Wang S S, et al. Impacts of spatial heterogeneity on crop area mapping in Canada using MODIS data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 119: 451-461.

[29] Cihlar J, Ly H, Li Z Q, et al. Multitemporal, multichannel AVHRR data sets for land biosphere studies-Artifacts and corrections[J]. Remote Sensing of Environment, 1997, 60(1): 35-57.

[30] Wang J, Xiao X M, Liu L, et al. Mapping sugarcane plantation dynamics in Guangxi, China, by time series Sentinel-1, Sentinel-2 and Landsat images[J/OL]. Remote Sensing of Environment. 2020, 247, [2020-06-09], https: //doi. org/10. 1016/j. rse. 2020. 111951.

[31] Savitzky A, Golay M. Smoothing and differentiation of data by simplified least squares procedures[J]. Analytical Chemistry, 1964, 36(8): 1627-1639.

Mapping of the winter crop planting areas in Huaihe River Basin based on Google Earth Engine

Pan Li1, Xia Haoming1,2※, Wang Ruimeng1, Niu Wenhui1, Tian Haifeng1, Qin Yaochen1,2

(1.,-,475001,; 2.(),,475001,)

The winter crop has been one of the important crop types in China. Accurate and timely spatio-temporal distribution of planting area directly determines the grain output and economy, as well as the national food security. Taking the Huaihe River Basin as an example, this study aims to extract the planting areas of winter crops according to the phenology period using the Google Earth Engine cloud platform and the fusion of Landsat-7/8 and Sentinel-2A/B images. Firstly, a dataset of time-series images was constructed with a spatial resolution of 30 m. A CFMask algorithm was selected to preprocess the images, thereby calculating the Normalized Difference Vegetation Index (NDVI). More importantly, the maximum NDVI of all high-quality images within 10 days was used to obtain time-series data with equal time intervals. The linear interpolation was utilized to fill the pixels without high-quality images. Savitzky-Golay (S-G) filtering (a second-order filter with a moving window of 9 observations) was adopted to smooth the NDVI time series for the removal of noise. As such, a smoothed NDVI time series was obtained with a 10-day interval. Secondly, the peak growth, sowing, and harvest periods were determined to select sample points of winter crops with different spatial distributions, according to the NDVI time series. Subsequently, the winter crops were sowed in mid-late October, when the NDVI values were the lowest. The NDVI values gradually increased, after the emergence of seedlings in early November. The crops stopped growing in January during the overwintering period, where the NDVI stayed the same over the whole period. Furthermore, the NDVI resumed growing and gradually reached the peak growth period, when the winter crops turned green in February. After that, the NDVI reached the peak at the heading stage, and then gradually decreased. Correspondingly, the NDVI dropped to the bottom, when the harvest was over from the end of May to June. According to these characteristics in the process of winter crops growth, the peak growth period was determined from March 20, 2018, to April 20, 2018, the sowing period was determined from October 11, 2017, to November 10, 2018, and the harvest period was determined from May 20, 2018, to June 30, 2018. Particularly, the maximum NDVI was achieved in the peak growth period and the minimum and median of NDVI in the sowing and harvest period. Finally, the classification model of a decision tree was constructed, according to the NDVI boxplots of winter crops and non-winter crops at different time periods. The planting area map of winter crops was also generated for the Huaihe River Basin. The results showed that the planting area of winter crops was 8.762×106hm2in the Huaihe River Basin in 2018. Specifically, the user accuracy was 0.926, the producer accuracy was 0.970, the total accuracy was 0.958, and the Kappa coefficient was 0.912. Consequently, the large-scale planting area of winter crops was extracted accurately for the decision-making in similar areas.

remote sensing; crops; mapping; planting area; Google Earth Engine; Huaihe River Basin

潘力，夏浩铭，王瑞萌，等. 基于Google Earth Engine的淮河流域越冬作物种植面积制图[J]. 农业工程学报，2021，37(18)：211-218.doi：10.11975/j.issn.1002-6819.2021.18.025 http://www.tcsae.org

Pan Li, Xia Haoming, Wang Ruimeng, et al. Mapping of the winter crop planting areas in Huaihe River Basin based on Google Earth Engine[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(18): 211-218. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2021.18.025 http://www.tcsae.org

2020-11-10

2021-08-30

河南省科技攻关计划项目（212102310019）；国家自然科学基金项目（41701433）；黄河文明省部共建协同创新中心重大项目（2020M19）；河南省自然科学基金资助项目（202300410531）

潘力，研究方向为农业遥感。Email: panli970611@henu.edu.cn

夏浩铭，博士，副教授，研究方向为定量遥感及其综合应用。Email: xiahm@vip.henu.edu.cn

10.11975/j.issn.1002-6819.2021.18.025

TP79

A

1002-6819(2021)-18-0211-08