基于GEE平台和自动统计分配算法的大范围冬小麦提取

赵亮，刘莉，司丽丽，赵铁松，黄敬峰

（1.河北省气象灾害防御和环境气象中心/河北省气象与生态环境重点实验室，石家庄 050021；2.浙江大学农业遥感与信息技术应用研究所，杭州 310058）

小麦是全球种植范围最广泛的粮食作物。作为中国第二大粮食作物，小麦种植面积占粮食总面积的25%，产量占粮食总产的22%，在中国粮食构成中占重要地位［1-3］。作为承灾体，许多农业气象灾害监测工作的开展、指标的制定都离不开作物的分布信息［4，5］。因此，在短时间内获取准确的小麦种植信息对保障粮食安全、水资源利用和碳循环研究、指导农业生产和防灾减灾工作都具有重要意义［6-8］。

传统的农业统计方法在很长一段时间一直被用于作物监测（如FAOSTAT），这些统计数据提供了高度准确的土地利用信息，但通常缺乏空间细节［9］，同时存在耗时费力、信息发布滞后等缺点［10，11］。自20世纪80年代后期，遥感技术被广泛用于作物类型识别、面积估算和物候信息获取［12-14］。赵叶等［15］基于2017年11月21日（小麦分蘖期）和2017年12月24日Landsat8 OLI影像，将MIR、NIR和RED波段进行HSV变换，提取河北省中南部越冬前冬小麦的面积。而Sentinel-2卫星具有幅宽大、重访周期短、数据免费、波段信息丰富以及较多的红边波段等特点，这些特点在农作物识别和面积提取方面具有很大的优势［16，17］。尹捷等［18］利用Sentinel-2A、MODIS以及高光谱珠海一号OHS-2A卫星等多源遥感数据，采用支持向量机算法提取河北省雄安新区的小麦信息，比较后发现Sentinel-2A对小麦的识别效果最佳，总体精度为85.57%，Kappa系数0.81。李方杰等［19］利用时序Sentinel-2遥感影像生成的归一化植被指数（Normalized Difference Vegetation Index，ND⁃VI），提出基于复合型混合演化算法和区域作物种植面积总量控制的NDVI时序相似性阈值，提取河北省衡水市武邑县的冬小麦面积总体精度达98.08%。

上述研究小麦信息提取精度普遍较高，但研究区相对较小，提取的也多为单年的信息。目前大面积农作物信息遥感提取困难主要有两个方面。①庞大的数据量以及由此产生的繁杂数据预处理和指数计算工作［20，21］，在利用高时空分辨率数据进行大范围农作物多年信息提取时，这种困难尤其明显。但随着地理大数据与云平台、云计算的发展，Google Earth Engine（GEE）作为一个基于云平台的全球尺度地理空间分析平台，为快速的数据预处理、下载以及遥感地物信息的提取带来了新机遇［22-24］。②同一地区基于不同数据来源的农作物种植面积不一致［9］。由于信息的获取方法不同，不仅来自不同研究者的结果会不一致，基于单纯遥感技术获取的农作物种植面积与农作物统计面积数据之间也会不一致。这些会在一定程度上影响农作物空间分布遥感信息在资源环境、生态、气候变化和粮食安全等地学领域中基于行政单元的空间分析应用与尺度转换［19］。因此，利用多源数据融合，以提高大范围作物遥感识别精度，确保作物面积总量与作物面积统计数据一致性，值得进一步研究。

综上，本研究借助GEE平台，将遥感数据与统计数据协同使用，设计自动统计分配算法，以分析基于不同植被指数合成影像的冬小麦面积及分布信息提取精度；同时利用筛选出的最优植被指数开展河北省多年的冬小麦信息提取工作。研究结果可为大范围小麦快速分布制图提供一定的技术方法参考和思路借鉴，为后续农作物空间分布遥感信息在资源环境、生态、气候变化和粮食安全等地学领域中基于行政单元的空间分析应用打下基础。

1 研究区概况与数据来源

1.1 研究区概况

河北省横跨华北、东北两大地区，位于东经113°27′—119°50′，北 纬36°05′—42°40′，总 面 积18.85万km2，地势西北高、东南低，由西北向东南倾斜。河北省属温带大陆性季风气候，大部分地区四季分明。年日照时数2 126~3 063 h，年无霜期80~205 d，年平均气温2.2~14.6℃，年降水量338.4~688.9 mm，降水量分布特点为东南多西北少，为作物生长提供了良好的环境［25］。

河北省是中国13个小麦主产省份之一，主要种植冬小麦-夏玉米，一年两熟制，另外根据土水条件和气候条件的差异，还存在少部分的单季玉米种植区和两年三熟区。其中冬小麦多集中分布在中南部海拔较低的平原地区（图1），而张家口市、承德市等北部地区鲜少种植。

图1 河北省地理位置和高程分布

1.2 数据获取与预处理

研究使用的数据主要包括哨兵2号卫星数据、耕地分布数据、验证样本数据、面积统计数据、数字高程数据以及行政区划数据等。

1.2.1 哨兵2号卫星数据遥感影像主要来源于哨兵2号（Sentinel-2）卫星。该卫星作为欧洲空间局（European Space Agency，ESA）“哥白尼”计划系列中的重要光学卫星组成部分，分为2A和2B两颗。每颗卫星的重访周期为10 d，两颗互补，形成5 d的重访周期。卫星上携带的宽幅高分辨率多光谱成像仪（MultiSpectral Instrument，MSI），幅宽为290 km，覆盖13个光谱带。不同波段的空间分辨率略有不同，分为10、20和60 m 3个尺度［26-28］；另外有3个QA波段，其中QA60波段包含云掩膜信息，用以去除影像中的卷云和厚云。Sentinel-2数据以其相对较高的时空分辨率和免费获取的特点被广泛应用于地物信息提取中。

本研究使用GEE平台数据库提供的Sentinel-2 L1C产品，该等级产品已经过数字高程模型校正以及亚像素多光谱配准，提供的是经过正射校正的大气顶（Top of Atmosphere，TOA）反射率。通过JavaS⁃cript编程语言，在GEE平台上进行了去云、遥感影像光谱指数计算、最大值选取、拼接、剪裁、转投影、重采样等预处理步骤，得到2016—2019年的30 m分辨率、河北省全覆盖的最大值合成影像，用于后续的小麦信息提取。以下是部分数据预处理方法详细介绍。

1）光学遥感影像云污染像素的去除。光学遥感影像在获取时易受到云、云在地面投影时的阴影以及气溶胶的影响［29，30］，严重影响地面信息的获取。因此本研究利用GEE平台上官方提供的Sentinel-2数据云掩膜算法，设置QA60波段所有像元Bit10和Bit11的值均为0，得到掩膜，对输入的符合时间和空间范围的影像数据集去除有云像元，之后合成目标年份最小云量影像。图2展示了利用GEE平台2018年11月15日至2019年3月31日共622景Sen⁃tinel-2影像、经去云后合成的河北省完整遥感影像。

图2 基于GEE平台合成的河北省Sentinel-2原始波段遥感影像（R/G/B=Band 4/3/2）

2）光谱指数的计算。对于复杂的地物覆盖区使用原始影像单一波段或者几个波段的信息进行单独分析往往具有很大的局限性，而将多光谱数据进行组合运算所得到的指数可将地物的差异放大或突出，特别对某些植被的长势、生物量、健康状况等具有一定的指示意义。这些通过组合运算所产生的指数被称为植被指数。传统的植被指数有比值植被指数（Ration Vegetation Index，RVI）、归一化植被指数（NDVI）、增强型植被指数（Enhanced Vegetation In⁃dex，EVI）［31］和二波段增强型植被指数（2-bands En⁃hanced Vegetation Index，EVI2）［32］等。本研究使用了后三者。其中，NDVI常被用于提取植被的生长状态、植被覆盖度等信息，其计算方法如下。

式中，NIR为近红外波段反射率，R为红光波段反射率。

而EVI相较于NDVI更加复杂，它引入对气溶胶敏感的蓝光波段来校正气溶胶对红光波段的影响，并且引入了冠层调整因子，使其对大多数冠层背景不敏感，其计算方法如下。

式中，G为增益系数，C1、C2为气溶胶阻抗系数，L为冠层背景调整因子，BLUE为蓝光波段反射率。计算过程中G、C1、C2、L分别取值为2.5、6.0、7.5、1.0。

EVI2则是一种于2008年被提出的基于增强型植被指数的变型，其在计算过程中未使用蓝光波段，并且保持了与增强型植被指数的较好相似性。EVI2的计算公式如下。

式中，c为蓝光反射率转化成红光反射率的转换系数，取值为2.08。图3展示了基于GEE平台，利用合成的去云Sentinel-2影像（2018年11月15日至2019年3月31日）计算的NDVI、EVI和EVI2分布图。1.2.2耕地分布数据本研究使用的河北省耕地分布数据来自2019年的地理国情数据，为第三次全国国土调查结果。原始数据由河北省第三测绘院提供，包括耕地、园地、林地、草地、商服用地、工矿仓储用地、住宅用地等共12类。将原始栅格重采样至30 m，与遥感影像保持一致。将栅格值为1的耕地像元提取出来，作为后续小麦信息自动统计分配格点化的基础。

图3 基于GEE平台和去云Sentinel-2影像计算的2019年河北省NDVI、EVI和EVI2最大值合成图

1.2.3 面积统计数据和验证样本数据通过收集和整理河北省农村统计年鉴上的县级小麦面积，得到整个河北省2016—2019年县级小麦面积参考数据。另外为了更好地比较不同植被指数的信息提取效果，评估小麦自动统计分配格点化结果在空间分布上的准确度，本研究利用GEE合成的整个河北省2019年冬小麦生长期最小云量影像（图2），根据不同地物的不同颜色、纹理等差异，基于目视解译和2020年实地考察结果（期间通过走访调查，记录了田块近2年的种植情况），在小麦主要种植区河北省南部和中部地区选取了冬小麦、水体、山体、人工地表和其他（裸地和其他植被）5类地物共907个矢量面（表1），各地类矢量面的分布（图4a）和特征（图4b至图4k）。河北省的山区多分布在西部，其他地类样本随机选取，均匀分布。

图4 基于GEE平台选取的5种样本地类矢量面的分布（a）、各地类原始影像特征（b-f）和各地类EVI2分布特征（g-k）

表1 验证样本数据相关信息

1.2.4 数字高程数据和行政区划数据利用先进的星载热发射和反射辐射仪全球数字高程模型（Ad⁃vanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model，ASTER GDEM）生产的数字高程数据，制作研究区底图，辅助验证样本的判断。该数据是根据NASA的对地观测卫星Terra的详尽观测结果制作完成的，空间分辨率为30 m。另外行政区划数据包括省、市、县三级行政区划数据，来源于国家基础地理信息系统中心，数据格式为ESRI Shapefile。

2 研究方法

2.1 试验设计

本研究的试验设计如下：①利用GEE平台进行Sentinel-2数据预处理和计算得到整个河北省2019年不同植被指数（NDVI/EVI/EVI2）影像（“1.2.1”部分）；②将各影像输入自动统计分配算法中，得到基于不同植被指数的小麦分布图；③利用“1.2.3”部分得到的验证样本（907个矢量面）对小麦分布图进行比较和检验，挑选出最优的植被指数；④基于最优的植被指数影像和自动统计分配算法，得到2016—2019年的河北省小麦分布图，分析近年河北省小麦的分布变化。

2.2 自动统计分配算法原理与方法

由于农作物的生长有着明显而特别的季节性特征，利用不同作物物候及生长发育在时间序列上反映出的植被指数特征进行作物信息识别。该方法已经成为农作物面积提取的一种重要方式。以EVI为例，冬小麦的植被指数在播种期最小，越冬前出现一个小峰值，在越冬后期出现波谷。返青后冬小麦种植区的植被指数会随着小麦生物量的增加而逐步增加，并在抽穗期至乳熟期达到最大，之后随着叶片变黄、麦穗成熟，植被指数又逐步下降［33］。河北省冬小麦约前一年10月下旬出苗结束，当年的4月末开始抽穗，6月初开始成熟收获。而在4月以前，其他作物鲜少种植，因此植被指数越大的像元为小麦的可能性越高［34］。所以自动统计分配算法的原理是自动将统计数据中的小麦面积分配给小麦可能性更高的像素，自适应调整小麦的分布，直到累积小麦面积接近统计数据为止。

研究利用GEE平台筛选前一年11月至当年3月间无云的最大值像元，合成覆盖整个研究区的30 m分辨率的植被指数（NDVI/EVI/EVI2）影像；之后利用耕地分布数据对得到的最大值影像进行掩膜，只保留耕地像元；利用Python编程语言通过分县自动设定并调整阈值，对耕地像元进行筛选，保留冬小麦可能性更高的像元，直到累积小麦面积接近统计数据则停止筛选，以实现冬小麦种植区多年分布信息的提取。

2.3 精度评价方法

利用“1.2.3”部分得到的验证样本（907个矢量面）作为参考数据，对基于不同植被指数影像和自动统计分配算法得到的小麦分布图进行比较和检验，以判断输入不同植被指数影像对小麦信息提取结果的影响。具体方法是统计基于算法得到的小麦像元与各真实地类重合的像元个数并计算百分比。另外，利用统计数据，对算法提取的2016—2019年冬小麦分布结果进行评估。验证指标有平均精度（Aˉx）、平均绝对误差（MAE）和均方根误差（RMSE），其计算方法如下。

3 结果与分析

3.1 最优植被指数挑选结果

统计基于算法得到的小麦像元与各真实地类重合的像元个数，并计算百分比。不同植被指数的小麦信息提取精度评价结果如表2所示。结果表明，基于NDVI、EVI和EVI2的小麦分布结果准确率差别不大。冬小麦真实样本总共有136 312个，其中EVI的用户精度为94.58%，制图精度为91.67%，Kappa系数为0.967 3，均为最大。EVI2的结果与EVI相似，而NDVI的 准确率 相对EVI和EVI2为最低，用户精度、制图精度和Kappa系数分别为94.10%、91.66%和0.966 6。最终确定EVI为自动统计分配算法的最优输入。但不难看出，无论是基于哪个植被指数，自动统计分配算法提取的小麦分布信息准确率均大于94.00%，总体精度均大于98.00%，Kappa系数均大于0.96，提取结果真实可靠。

表2 不同植被指数2019年冬小麦分布结果精度评价混淆矩阵

3.2 河北省冬小麦提取结果精度分析

将GEE平台得到的河北省EVI影像输入自动统计分配算法中，得到2016—2019年河北省30 m分辨率的小麦分布（图5）。由图5可以看出，河北省的小麦主要分布在中南部地区，廊坊市和秦皇岛市的种植密度较低，承德市和张家口市没有冬小麦种植区。

图5 2016—2019年的河北省冬小麦空间分布

分别计算2016—2019年河北省所有小麦种植县冬小麦面积的平均精度，见表3。由表3可知，2016—2019年河北省冬小麦面积平均精度均大于97.00%，算法提取的整个河北省小麦面积与统计面积高度一致。另外计算所有县的MAE和RMSE，并输出散点图，输出结果如图6所示，算法提取结果与统计结果有很好的线性关系。

图6 2016—2019年的河北省冬小麦县级精度验证散点图

表3 冬小麦面积精度验证结果

3.3 河北省冬小麦分布变化

计算每个像元2016—2019年冬小麦种植年份数，得到整个河北省冬小麦分布变化图；另外，统计河北省各个市每年的小麦面积并制作折线图，结果如图7所示。研究发现，2016—2019年，河北省冬小麦面积从224.595万hm2上涨至227.956万hm2，增长幅度不大，2018年达到最高，为231.748万hm2。相较2016年，邢台、衡水和秦皇岛市的小麦种植面积总体呈增长趋势，其中衡水市涨幅最大，2019年相较2016年增长了7.180万hm2，呈由北向南扩张趋势；保定和沧州市4年内的小麦种植面积基本保持不变；而其他市的小麦种植面积总体呈下降趋势，其中石家庄市下降最多，2019年相较2016年减少了2.683万hm2，减少区域主要为石家庄西部地区。

图7 2016—2019年的河北省各市冬小麦面积和分布变化

4 讨论与结论

4.1 研究的优势、不足及下阶段计划

从卫星产品估算的农田面积通常与统计数据不一致，会阻碍农田地图在资源环境、生态、气候变化和粮食安全等地学领域中基于行政单元的空间分析应用。本研究采用遥感数据和统计数据协同使用的自动统计分配算法，制作了河北省多年的30 m分辨率小麦分布图。算法有以下3项优势。①解决了训练样本挑选难的问题。传统的大范围作物分布制图通常需要相对大量的训练样本来训练模型，训练样本的挑选通常需要较大的时间和精力投入，同时训练样本的准确性也决定了算法的准确性，而本算法不需要训练样本，能够快速在其他年份复制使用。②研究基于GEE平台进行数据预处理，极大地提高了制图效率。作为可免费访问、基于云数据的地理空间计算平台，GEE平台为高效处理大量遥感数据提供了新途径［35］，随着数据可用性和强大计算资源的增加，越来越多的遥感信息提取工作在GEE平台开展［36，37］。GEE样本选取方便及精度验证简单等优点让研究者能专注于开发和改进自动统计分配算法本身，提高算法信息提取的精度，而不需要过多关注一些重复的技术性工作。③极大提高了作物分布制图的准确性，使之与官方统计数据一致，为产量预测、水资源利用以及农作物气象灾害监测等后续研究打好基础。

本研究使用的验证样本是矢量面，虽然在采集过程中尽量保证矢量面内地物的一致性，但由于人为主观因素的存在，加上混合像元的问题，验证样本仍然存在一定的不准确性。另外，本研究使用冬小麦生长期内植被指数的变化差异作为基础进行小麦的提取且精度较高，是因为在相同的生长季内（前一年11月至当年3月），河北省鲜少有其他作物大面积被种植。这使得算法在与河北省种植制度相似的地区复刻，同样能获得较好的结果，但在冬季作物种类偏多的地区，例如南方冬小麦和油菜的混合种植区，为获得较高的精度，自动统计分配算法需添加额外的条件，这也是下阶段算法的改进目标。同时，不同遥感数据的联合使用也值得进一步研究。

4.2 结论

为快速准确获取大范围小麦信息，本研究基于Google Earth Engine云平台进行Sentinel-2数据的预处理和植被指数计算，得到整个河北省2019年不同植被指数（NDVI/EVI/EVI2）影像；输入自动统计分配算法，得到基于不同植被指数的小麦分布图并比较，挑选出最优的植被指数；基于最优的植被指数影像和自动统计分配算法，快速得到2016—2019年河北省冬小麦的空间位置分布信息。主要结论如下。

1）利用GEE平台辅助提取作物信息，有数据获取容易、计算速度快、样本选取方便等优点。GEE能够快速完成整个河北省多年影像的去云、镶嵌、裁剪及指数计算等预处理工作，较本地处理有明显优势。

2）NDVI、EVI和EVI2三者的小麦分布制图准确率差别不大，均大于94.00%，其中EVI的精度最高。河北省所有小麦种植县2016—2019年的平均精度均大于97%。相较2016年，衡水市小麦面积4年增长了7.180万hm2，主要呈由北向南扩张趋势；石家庄市减少2.683万hm2，减少区域主要为西部区域。

3）自动统计分配算法不需要依赖训练样本，极大地降低了大面积作物分布制图的成本，且能快速地在其他年份复制使用，提高了大面积作物分布制图效率。另外，算法将遥感数据和统计数据协同使用，获得的小麦面积和分布与真实情况有良好的一致性。本研究可为大范围小麦分布制图提供一定的技术方法参考和思路借鉴，为指导农业生产、农业气象灾害监测和防灾减灾工作提供基础。