这是用户在 2024-9-5 19:01 为 https://app.immersivetranslate.com/pdf-pro/909e2084-9d82-4fc7-b12f-dcf4cc3a1698 保存的双语快照页面，由 沉浸式翻译 提供双语支持。了解如何保存？
2024_09_05_4cf61db1cfdf7e4a9692g


将低分辨率图像与农业统计数据结合，以生成亚像素作物类型地图和协调的面积估算


琼 , 何银 , 马克·A·弗里德尔 , 梁志友 , 赵亮 , 华军唐  文彬
湖北省地理过程分析与模拟重点实验室/华中师范大学城市与环境科学学院 武汉 430079，中国
波士顿大学，地球与环境系，美国马萨诸塞州波士顿，邮政编码 02215
肯特州立大学地理系，美国俄亥俄州肯特市林肯街 325 号，邮政编码 44242
国际食品政策研究所，华盛顿特区 I 街 1201 号，邮政编码 20005，美国
华中农业大学经济管理学院宏观农业研究所，中国武汉 430070
农业遥感重点实验室（AGRIRS），农业农村部/中国农业资源与区划研究所
中国北京市 100081 农业科学院


文章信息

 关键词：

 作物类型映射
 随机森林回归

迭代区域间隙空间分配 MODIS
 面积估算
 特征选择

 摘要


摘要 可靠的作物类型地图对于农业监测、确保粮食安全和环境可持续性评估至关重要。由于其短重访周期，MODIS 等低分辨率影像被广泛用于作物类型制图，这对于检测不同作物类型的季节动态具有优势。然而，固有的低空间分辨率可能限制其在农业景观异质性地区的作物类型制图的实用性。农业统计数据提供了不同空间和时间尺度的作物面积信息，有潜力改善遥感作物类型制图。然而，以往的研究通常将农业统计数据作为参考数据来评估卫星衍生作物类型地图的准确性，但很少利用它们来改善作物类型分布制图。将农业统计数据与卫星影像结合以生成高精度作物类型地图的实用性鲜有探讨。本研究提出了一种通过整合 MODIS 时间序列和农业统计数据来绘制亚像素作物类型分布的方法。 我们在中国农业生产最高的黑龙江省测试了我们的方法。首先，我们使用了优化的随机森林回归（RF-r）模型，训练样本来源于高分辨率图像（即 SPOT 和 Landsat），以预测 MODIS 时间序列中的亚像素作物类型分布。为了优化 RF-r 模型，使用了 2011 年 8 天的 MODIS 时间序列五个植被指数作为候选自变量，并实施了向后特征消除策略，以选择最佳变量进行模型预测。其次，我们开发了一种迭代区域差距空间分配（IAGSA）方法，以空间上调和从 MODIS 基础地图估算的作物面积与农业统计数据之间的差异。我们发现，MODIS 衍生的作物分数与高分辨率图像衍生的作物分数一致，所有作物类型的差异为 ，但 MODIS 估算的作物面积与农业统计数据之间存在明显差异。 由 IAGSA 调整的亚像素作物类型图不仅与作物面积的农业统计数据一致，而且保留了原始 MODIS 衍生作物比例的空间分布模式。我们的结果表明，整合低分辨率图像和农业统计数据以绘制亚像素作物类型分布的优势，并提供一致的作物面积估算。所提出的方法具有以成本效益高的方式绘制区域大规模作物类型范围的潜力。

https://doi.org/10.1016/j.rse.2021.112365

收到日期：2020 年 4 月 18 日；修订后收到日期：2021 年 2 月 16 日；接受日期：2021 年 2 月 19 日

在线可用日期：2021 年 3 月 12 日

0034-4257/© 2021 爱思唯尔公司。保留所有权利。

 1. 引言


世界人口预计将从 2019 年的 77 亿增加到 2050 年的 97 亿（联合国，2015），在有限的农业用地上养活这一不断增长的人口是一个巨大的挑战（Cole 等，2018；Zabel 等，2019）。极端事件的增加，如洪水、干旱和霜冻，使这一情况更加恶化（Calzadilla 等，2013；Foley 等，2011）。这些挑战表明，迫切需要可持续管理农业用地。制定可持续农业战略的一个关键组成部分是关于农业用地的空间分布和面积估算的精确信息（Dong 等，2015；Fritz 等，2015；Waldner 等，2015）。在这方面，遥感已成为作物监测的重要工具，因为它具有综合的空间视角、频繁的数据获取和相对较低的观测成本（Song 等，2017a；Wang 等，2019；Xiao 等，2005；Yin 等，2020a）。 例如，地球观测组织（GEO）的农业监测实践社区与综合全球观测战略（IGOS）呼吁建立一个全面、操作性强且准确的全球农业监测系统，利用地球观测数据（Whitcraft et al., 2015）。

粗分辨率卫星图像已广泛用于在区域或更大尺度上分类不同的作物类型（Massey et al., 2017; Skakun et al., 2017; Xiao et al., 2006）。粗分辨率图像的高时间可用性使其能够很好地描绘作物类型在生长期间的光谱动态，这对于捕捉不同作物类型之间微妙但关键的物候差异具有良好的前景（Foerster et al., 2012; Zhang et al., 2015）。例如，中等分辨率成像光谱仪（MODIS）数据具有短的重访周期（每天在 30 度以北），丰富的光谱波段（36 个波段）以及在中等空间分辨率下的全球覆盖（ ），在大面积区分作物类型方面具有很大的潜力（Chen et al., 2018; Wardlow and Egbert, 2008; Zhang et al., 2014; Zhong et al., 2016）。

然而，使用低分辨率图像的时间序列进行作物类型制图面临几个挑战。首先，简单地将所有图像叠加在一起可能无法产生最佳的分类精度，因为信息冗余（胡等，2016；邱等，2017）。此外，从时间序列生成的广泛使用的增强植被指数（EVI）和归一化差异植被指数（NDVI）可能不是作物类型制图的最佳指标。其他植被指数，如陆地表面水分指数（LSWI），在分类不同作物类型时可能更为关键，特别是在作物类型多样且相关种植模式复杂的地区（肖等，2006；钟等，2014）。因此，选择最佳的光谱和时间特征进行作物类型分类，以充分利用时间序列图像是非常重要的。其次，低分辨率图像由于亚像素混合可能会在作物面积估算中产生相当大的误差（杜维耶和德福尼，2010；瓦尔德纳和德福尼，2017）。

为了解决“混合像素”问题，以往的研究基于作物比例与遥感变量之间的关系，绘制了亚像素作物覆盖率。例如，Lobell 和 Asner（2004）使用了一种概率线性解混模型，利用 MODIS 时间序列估计墨西哥雅基谷的作物比例及其不确定性。Pan 等（2012）开发了一种基于物候的指数，利用 EVI 时间序列和回归模型估计中国两个农业区内 250 米 MODIS 像素中的冬小麦比例。除了线性光谱解混模型，Atzberger 和 Rembold（2013）使用神经网络研究 NDVI 时间序列与 AVHRR 像素内作物面积之间的非线性关系。然而，这些解混模型完全依赖于地球观测数据，因此它们的准确性不可避免地受到训练样本的可靠性、遥感图像的质量以及模型的学习和泛化能力的限制（Foody, 2002; Verburg 等, 2011）。

将遥感影像与农业统计数据结合提供了一种替代方法，以提高作物类型制图和面积估算的准确性。由于记录的独立性和时间的连续性，农业统计数据已成为基于卫星的土地利用和土地覆盖研究的重要补充和参考数据（Pittman et al., 2010; Ramankutty et al., 2008）。农业统计数据通常作为参考数据，用于评估不同行政级别（如县、地级和国家级）遥感衍生的作物类型地图的准确性（Chang et al., 2007; Xiao et al., 2005）。农业统计数据还有潜力改善遥感作物空间分布模型的预测，但这种潜力很少被探索。此外，基于遥感的面积估算往往与农业统计数据不匹配。这是一个不幸的情况，因为国际组织，如联合国机构（例如，粮农组织）、世界银行和 CGIAR，无法使用基于遥感的地图，除非它们与农业统计数据一致（Erb et al., 2007; Fritz et al., 2015）。 因此，将遥感与农业统计数据相结合不仅有价值，而且是必要的，以提高基于遥感的地图在农业应用中的实用性，特别是考虑到这些地图可能存在较大的分类错误。

在这里，我们提出了一种将低分辨率图像与农业统计数据相结合的方法，以绘制亚像素作物类型分布并生成协调的面积估算。我们在中国黑龙江省测试了这种方法，该省的农业在国家粮食安全中发挥着至关重要的作用。更具体地说，我们的目标是：（1）评估结合 MODIS 时间序列、RF-r 模型和向后特征消除策略在绘制亚像素作物类型分布方面的潜力；（2）引入一种迭代面积差距空间分配（IAGSA）方法，以协调 MODIS 衍生的作物面积与农业统计数据之间的面积差异，从而提高基于 MODIS 的作物类型图的准确性。

 2. 材料


2.1. 研究地点的描述


我们选择了中国东北的黑龙江省作为研究地点（图 1）。黑龙江的耕地面积约为 1593 万公顷，主要集中在西部的松嫩平原和东部的三江平原。由于其肥沃的黑土（即草甸土）以及高水平的农业机械化，黑龙江自 2011 年以来一直是中国最大的粮食生产省份，2019 年生产了约 6018 万吨粮食（中国国家统计局，2019）。该地区气候以温带大陆性季风气候为主，夏季短暂而温暖，冬季漫长而严酷。由于光照时间和积温的限制，该地区以单作模式为主。主要作物类型为水稻、玉米和大豆，这三者合计约占总耕地面积的 （中国国家统计局，2019）。玉米通常在四月底播种，九月中旬收获。水稻的种植时间为四月初至九月底。大豆的生长周期与玉米相似，但具有独特的年度内发育阶段（例如）。作物的季节动态在不同生长阶段（如开花和结荚阶段）存在明显差异，相关的物理和生化特性也可以通过影像时间序列的光谱特征轻松识别。

 2.2. MODIS 时间序列


我们使用了 2011 年的 8 天复合 MODIS/Terra 地表反射率产品（MOD09A1，集合 6）作为亚像素作物类型制图的主要数据。MOD09A1 产品包括七个光谱波段，名义空间分辨率为 （实际大小为 463.3 米），已针对大气气体和气溶胶进行了校正。

图 1. 研究区域的位置。(a) 从 SPOT4、SPOT5 和 Landsat TM/ETM+解读的作物类型参考图（刘等，2014）。(b) 从 GlobeLand30 数据集中提取的黑龙江省农田图（陈等，2015），以及用于生成图 1(a)中作物类型参考图的 2011 年实地调查的作物类型样本。

（Vermote，2015）。我们下载了从朱利安日期 65 到 305 的所有 MOD09A1 图块，这些图块覆盖了黑龙江省所有作物类型的整个生长季节。计算并用于本研究的五个植被指数，即增强植被指数（EVI）、土壤水分植被指数（LSWI）、归一化差异衰老植被指数（NDSVI）、归一化差异耕作指数（NDTI）和绿色植被指数（VIgreen），这些指数由每个时间点的不同光谱波段组成（见表 1）。我们选择这五个植被指数是因为它们已被证明是检测特定作物类型不同物理和生化特性的良好指标。最终生成了 155 个特征（31 个时间点×5 个植被指数），作为 RF-r 模型的候选变量。


2.3. 作物类型参考图和辅助数据


用于训练 RF-r 模型（在 3.1 中描述）和评估其性能的参考数据来源于 2011 年黑龙江省的现有作物类型地图（图 1a）。这些地图由黑龙江省农业科学院生成（刘等，2014）。具体而言，首先对高分辨率影像（即 SPOT4、SPOT5 和 Landsat TM/ETM +，以下简称 HSRI）进行了最大似然分类。然后，基于对作物分布和物候特征的先验知识，进行了专家视觉解读和检查，以纠正明显的分类错误（刘等，2014）。为了评估这些 HSRI 地图的可靠性，我们使用了 287 个作物类型样本（即 95 个玉米样本；120 个水稻样本；72 个大豆样本）在现场收集，以评估 HSRI 地图的制图精度。我们发现玉米、水稻和大豆的生产者精度分别为 和 ，而它们的用户精度分别为 、 和 ，这表明 HSRI 地图的可靠性。 为了获取基于 MODIS 的作物类型制图样本，我们将基于 HSRI 的作物类型地图重投影到 MODIS 正弦投影，并将其叠加在 MODIS 网格上。随后，估算每个 MODIS 网格的作物比例为 种作物类型像素的数量除以 网格内所有 像素的总数。对于每种作物类型，我们从 作物比例图中随机选择 4000 个像素来训练 RF-r 模型，另选 4000 个像素进行准确性评估。

除了基于 HSRI 的作物类型地图，我们还使用了 2010 年 GlobeLand 30 数据集中的耕地层（Chen et al., 2015）来掩盖 MODIS 像素网格中的非农业区域（见图 1b），并设定一个限制条件，即 IAGSA 调整后的作物比例应小于每个 MODIS 像素内的耕地比例（见第 3.2 节）。在这里，我们假设 2011 年的耕地与 2010 年相同，因此可以在本研究区域使用 GlobeLand30 产品作为我们的耕地掩膜。GlobeLand30 是一个全球土地覆盖产品，具有 空间分辨率（http://www.globallandcover.com），其用户和生产者的准确度在 2010 年黑龙江省的耕地上均大于 （Wang et al., 2018）。我们基于 GlobeLand30 产品，使用与基于 HSRI 的作物类型地图相同的方法计算了黑龙江省每个 网格内的耕地比例。
 表 1

本研究使用的植被指数。波段 1-7 指的是 MODIS 表面反射率中的红光（B1， ）、近红外（B2， ）、蓝光（B3， ）、绿光（B4， ）、短波红外（B6，1.628-1.625 ）和短波红外（B7，2.105-2.155 ）波段。
 光谱通道
 植被指数
  公式
物理和生化性质		  参考
 可见光-近红外 EVI
植被生长状况，冠层结构		  沃德洛和埃格伯特，2008
NIR-SWIR LSWI  水分含量，  肖等，2006
 可见短波红外 NDSVI  残留覆盖
植被生长状态、水分含量、残留覆盖		  钟等, 2014
SWIR-SWIR NDTI
非光合成成分，		  钟等, 2014
 可见-可见 VIgreen  残留覆盖
叶绿素浓度，植被状态		  Löw 等, 2013


2.4. 农业统计


我们使用了中国国家统计局发布的《中国农村统计年鉴》中省级的作物统计数据（可在https://data.cnki.net/Yearbook/Single/N2011110092获取），数据涵盖 2010 年、2011 年和 2012 年。对于县级作物统计数据，我们使用了黑龙江省统计局的黑龙江农村统计年鉴中的数据集（2010-2012 年）。这两个数据集分别提供了 80 个县每种作物类型的播种面积和黑龙江省每种作物类型的总面积。为了减少单一年份农业统计数据的不确定性，我们对相邻三年的作物播种面积（即 2010-2012 年）进行了平均，以代表 2011 年每个县的作物类型面积。我们将农业统计数据与 MODIS 衍生的作物类型面积进行了比较，以评估 MODIS 衍生的作物分数图的表现。此外，我们使用农业统计数据计算了统计衍生的作物类型面积与 MODIS 衍生的作物类型面积之间的面积差异，然后进行了空间分配，以提高 MODIS 衍生的作物类型图和面积估算的准确性。

 3. 方法论


所提出的将 MODIS 时间序列与农业统计数据相结合的方法主要包括两个部分（图 2）。第一部分，我们使用 RF-r 模型预测整个研究区域内每种作物类型的亚像素作物比例，通过建立作物比例与 MODIS 光谱特征之间的关系。通过使用最佳光谱-时间特征作为自变量，采用向后消除特征选择策略确定了最佳 RF-r 模型。第二部分，我们通过使用参考统计数据来修正 MODIS 衍生的作物类型地图的偏差。我们开发了一种方法（即 IAGSA），使用迭代线性分配规则对区域偏差进行空间分配，并生成了新的亚像素作物类型分布图。最后，我们评估并比较了这两个部分衍生的作物比例图。


3.1. 作物比例估算


3.1.1. 随机森林回归


随机森林是一种集成机器学习算法，它允许基于使用训练集观察数据估计的模型来预测响应变量（Breiman, 2001）。RF-r 模型算法提供了回归准确性的内部度量（即袋外误差，称为 OOB 误差）和每个变量的相对重要性，这支持了变量贡献的定性分析，并为理解建模的复杂关系提供了有用的信息（Rodriguez-Galiano 等，2012）。

图 2. 整合 MODIS 时间序列和农业统计数据以生成亚像素作物类型图和调整后的面积估算的工作流程。VIs 分别代表植被指数和时间点。

RF-r 模型使用基本回归公式 进行了说明，在本研究中，因变量 是从基于 HSRI 的作物类型图估算的每个 500 米像素内的作物比例，自变量 是光谱-时间特征， 表示候选特征的数量，最初为 155，然后通过我们的向后特征消除策略逐渐减少（见第 3.1.2 节）。RF-r 模型的预测可能会受到“回归到均值”的影响（Huang 等，2016）。为了纠正这种偏差，我们使用了一种线性旋转方法（Huang 等，2016；Zhang 和 Lu，2012），该方法在补充材料的第一部分中进行了说明。使用这种方法，本研究中 RF-r 模型的回归偏差得到了减少（补充材料的第二部分）。


3.1.2. 向后特征消除


在这里，我们使用了一种向后特征消除方法，通过逐步消除最不重要的变量来选择 RF-r 模型的最佳变量集（Guan 等，2013；Yu 和 Liu，2004）。首先，我们使用原始的 155 个光谱时序特征作为自变量来执行初始的 RF-r 模型。其次，我们对所有变量的特征重要性分数进行排名，并移除得分最低的 个特征。第三，使用袋外（OOB）样本（不包含在训练子集中的数据）及其相关的 OOB 误差，在内部交叉验证过程中评估在移除特征后 RF-r 模型的表现（Breiman，2001；Belgiu 和 Drăguţ，2016）。为了确保稳健性，我们重复上述步骤 30 次。获得最佳准确性（即最低 OOB 误差）的特征被确定为 RF-r 模型的最佳变量。随机森林生成的“%InMSE”指标表示当特定变量被删除时均方误差的增加，用于表征特征的重要性分数。 我们将回归树的数量（ntree）设置为 500，终端节点的最小大小（nodesize）设置为 10。在每个节点尝试的不同预测变量数量（mtry）设置为默认值，即总变量数量的三分之一 （mtry ）。由于总变量数量 在迭代中有所不同，因此随着特征通过消除过程被淘汰，mtry 也会有所变化。

图 3. IAGSA 方法的工作流程，用于空间分配 MODIS 衍生的作物类型面积与农业统计之间的面积差距。


3.2. 作物比例调整


我们提出了一种迭代区域间隙空间分配（IAGSA）方法，以空间调整 MODIS 衍生的亚像素作物分数的面积（图 3）。该算法包括以下五个步骤：

首先，我们通过将每个县的像素面积乘以总像素数，按照以下公式计算了每种作物类型的 MODIS 衍生作物面积 MArea ：
MArea

其中 和 分别代表每个县（例如，北安市）和作物类型（例如，玉米）。 代表特定的 MODIS 像素， 代表像素 中作物类型 的作物比例， 代表县 内的 MODIS 像素总数， 代表 MODIS 像素的面积（ ）。

然后，根据以下公式计算县 中 MODIS 衍生的作物类型面积与作物类型 的统计数据之间的面积差 ：
 阿卡普

其中 MArea 表示从 MODIS 衍生的作物分数计算的县级作物类型面积，CArea 是来自县级统计数据的作物播种面积。

第三，我们计算了目标县 中作物类型 的农田比例 以下的像素数量 ，使用的公式如下：


其中 表示作物类型 在像素 内的作物比例。 是像素 内的耕地比例，该比例是从 GlobeLand30 数据集中估算得出的（详细信息见第 2.3 节）。由于耕地比例应大于每个像素内的作物比例， 作为 ptki 的上限。

如果 Agap 小于 或者 等于 0，则不进行调整；否则，我们估计了目标县 中作物类型 的区域调整参数 ，其定义为：
Ap

(5) 如果 小于 1，则子像素裁剪比例乘以 。如果 大于 1，我们首先将子像素裁剪比例乘以 ，然后将值大于上限 的像素的裁剪比例更改为 。步骤（1）到（4）重复进行，直到 Agap 小于 或 等于 0。该过程的输出是目标县 中每种作物类型 的调整后裁剪比例。请注意，如果在调整过程中 等于 0 的条件在 Agap 变小于 之前满足，则遥感与农业统计之间的面积差异将不会完全按预期分配。为了进一步阐述我们的 IAGSA 方法，我们提供了三个示例，说明了 IAGSA 方法如何调整作物比例（补充材料的第三部分）。由于我们旨在绘制整个黑龙江省的主要作物类型，因此变量 和 需要分别包括三种作物类型（即玉米、水稻和大豆）和 80 个县。


3.3. 性能评估


在进行 IAGSA 调整之前，评估 MODIS 衍生的作物分数图的准确性有两种方法。首先，我们使用了 4000 个测试样本，这些样本与用于校准 RF-r 模型的训练样本是独立的，以基于两个指标评估映射准确性：决定系数 和均方根误差（RMSE）。其次，这些未经过 IAGSA 调整的作物分数图通过将每个像素值乘以其面积转换为作物面积，然后将图生成的面积估算汇总并与县级农业统计数据进行比较。这些比较通过估算 和各县的 RMSE 进行了量化。

为了评估我们 IAGSA 方法的性能，我们首先评估了 IAGSA 方法是否能够显著改变未应用 IAGSA 方法的作物分数的原始空间分布模式。为此，我们分析了应用 IAGSA 方法后作物分数的空间分布如何变化，并比较了使用和不使用 IAGSA 调整的 MODIS 数据所得到的作物分数直方图。此外，我们测试了三种低代表性训练样本的场景，以评估在遥感方法较差时 IAGSA 方法的性能。这三种场景中的训练样本在数量上既不足，也在空间上不离散（补充材料的第四部分）。我们基于 和 RMSE 比较了这些场景中使用和不使用 IAGSA 调整的作物分数图的准确性，这些数据来自于独立收集的 4000 个 HSRI 地图验证样本。最后，由于中国的农业统计通常在三个统计单位层级上记录，即。在县、地级和省级层面，我们进行了测试，将县级统计数据替换为省级统计数据，以进行 IAGSA 分析，以了解统计单位空间尺度对 IAGSA 方法性能的影响。在这些测试中，省级 IAGSA 只有一个区域差距，而县级 IAGSA 需要分配 80 个区域差距。然后比较和评估了县级和省级 IAGSA 方法之间作物比例的空间变化程度。

 4. 结果


4.1. 作物类型制图的最佳光谱-时间特征


图 4 显示了通过向后特征消除策略选择的作物类型分类的最佳光谱-时间特征。我们发现玉米、水稻和大豆的最佳特征数量分别为 103、43 和 76，这表明需要更多特征来将玉米和大豆与其他类别区分开，而水稻则需要较少的特征进行分类（图 4）。这一发现与 Yin 等人（2020b）的先前研究一致，该研究也表明，在黑龙江，使用比其他作物类型更少的光谱-时间特征可以实现令人满意的水稻分类准确性。我们还观察到，随着额外特征的去除，回归准确性提高（即 OOB 误差降低），直到达到最佳准确性（图 4）。此外，随着用于分类的特征数量减少，计算时间显著减少。我们发现，与使用所有 155 个特征相比，使用最佳光谱-时间特征使玉米、水稻和大豆的计算时间分别减少了 和 。

特征重要性分数（即，InMSE%）来自 RF-r 模型和向后特征消除方法，见图 5。结果显示，更多高分数的红色网格与 LSWI 和 NDTI 相关，而与其他植被指数的关联较少，这表明 LSWI 和 NDTI 在亚像素玉米制图中比其他植被指数更为重要（图 5）。在所有光谱时间特征中，{LSWI，年中的第 141-181 天}、{NDTI，年中的第 133-181 天}和{VIgreen，年中的第 149-172 天}提供了最重要的玉米制图信息。这可能是因为在黑龙江，年中的第 140 到 170 天是玉米的三叶期，此时观察到的玉米田光谱特征是低绿植被和土壤背景的混合。

图 4. 在向后特征消除方法的每次迭代中，随着不重要特征的移除，回归准确性和计算时间的变化，针对(a)玉米，(b)稻米和(c)大豆。灰色条表示具有最低 OOB 误差的最佳特征数量。

与其他土地覆盖类型不同。对于水稻，最佳的光谱时间特征是{NDSVI, DOY 133-157}、{NDSVI, DOY 和{LSWI, DOY 220-250}。DOY 在 130 到 150 之间是水稻移栽阶段，此时稻田的水分含量远高于其他作物类型。DOY 在 220 到 280 之间是水稻的衰老阶段，黄叶和干枯的叶子与其他作物类型的光谱差异显著。对于大豆，LSWI 比其他植被指数更为重要。在光谱时间特征方面，{LSWI, DOY 133-189}、{LSWI, DOY 253-301}、{EVI, DOY 221-245}和{NDSVI, DOY 是大豆识别中最关键的特征。DOY 在 210 到 240 之间对应黑龙江大豆的灌浆阶段，此时黄叶和独特的种子纹理使得光谱特征与其他作物类型有显著差异。


4.2. MODIS 衍生的作物分布图


黑龙江省的作物分布图是基于优化的 RF-r 模型和未进行 IAGSA 调整的 MODIS 数据生成的，如图 6 所示。我们发现玉米主要种植在松嫩平原和三江平原。稻米主要种植在三江平原，子像素分数为 。与玉米和稻米相比，大豆的种植面积较小，主要分布在西北地区。图 7 显示了在进行 IAGSA 调整之前，基于 MODIS 的作物分数估计与基于 HSRI 的作物分数之间的比较。总体而言，这两种估计结果相互一致。稻米的 最高，RMSE 最低（RMSE ），而玉米的 和斜率（斜率 ）最低，大豆的 RMSE 最高（RMSE =12.89）。我们还发现，MODIS 推导的稻米分数在不同分数区间的变异性高于玉米和大豆分数。这可能是因为黑龙江的稻田比玉米和大豆田大，导致估计的稻米分数较大。
(a)

后向特征消除
(b)

图 5. 从 RF-r 模型和向后特征消除方法得出的特征重要性分数（InMSE%）。(a) 使用所有 155 个特征的第一次迭代结果。(b) 从 30 次迭代中识别出的最佳特征。图表的横轴表示五个植被指数，纵轴表示从 DOY 65 到 DOY 305 的影像日期，间隔为 8 天。网格单元指的是相应光谱时序特征的 InMSE 值。更高的 InMSE%值表示该特征在作物类型分类中的重要性更高。NF 代表特征数量。

图 6. 2011 年黑龙江省基于优化的 RF-r 模型和 MODIS 数据生成的作物分布图。

图 7. 黑龙江省 MODIS 衍生与 HSRI 衍生的作物比例。结果在每 10%的区间（即 等）中进行了总结。点表示每个区间的中位作物比例（%）。条形表示 HSRI（横向）和 MODIS（纵向）的第 25 和第 75 百分位数。

在每个 500 米像素内，MODIS 衍生的玉米和大豆比例在作物覆盖率为 的地区略有高估，而在作物覆盖率为 的地区则被低估。训练数据的不均匀分布模式以及黑龙江的景观异质性可能导致了这些偏差。MODIS 衍生的作物比例直方图显示， 的像素具有作物覆盖率 ，而且少于 的像素具有覆盖率 （补充材料的第 V 部分）。在 4000 个训练样本中也发现了相同的模式，其中具有比例值 的像素数量几乎占所有样本的 （补充材料的第 V 部分）。

图 8 显示了从 MODIS 时间序列估算的农作物面积与县级农业统计数据估算的农作物面积之间的比较。在 80 个县中，所有三种农作物类型的 值 均已获得。与 HSRI 地图的比较一样，与农业统计数据的比较显示，稻米在三种农作物类型中具有最高的准确性， 为 0.98，RMSE 为 。相比之下，与农业统计数据相比，MODIS 推导的大豆比例的准确性低于其他两种农作物类型， 为 0.83，RMSE 为 。我们还发现，与农业统计数据相比，MODIS 推导的玉米面积在大多数县中被低估。


4.3. 通过 IAGSA 方法调整的作物分数图


在三种作物类型中，玉米在应用 IAGSA 方法后空间范围变化最大，黑龙江省有 个像素的玉米比例增加（图 9）。使用 IAGSA 调整后，稻米比例相对稳定。稻米像素的比例（即 ）保持不变的比例远大于其他作物类型（玉米为 ，大豆为 ），稻米比例的变化范围为 到 。稻米的这种适度变化归因于 MODIS 衍生估计与农业统计之间的小面积差异（图 8b）。图 9 还显示西部地区大豆比例下降，而黑龙江东南部地区则普遍增加。总体而言，减少大豆比例的像素数量多于增加的像素数量，如图 9 中的大豆直方图所示。

使用 IAGSA 方法调整的亚像素作物类型图不仅与作物面积的农业统计数据一致，而且与未使用 IAGSA 调整的原始空间分布模式也非常吻合（图 10）。我们发现这两种估计的直方图在列长度和横轴变化模式上都相似。对于像稻米这样面积差距较小的作物类型，相应直方图之间的差异可以忽略不计。


4.4. IAGSA 方法的规模效应


使用 IAGSA 方法在省级调整的三种作物类型的作物分数如图 11 所示。我们发现 IAGSA 的性能依赖于所选的尺度（图 9 和图 11）。由于基于遥感的作物面积与农业统计之间的面积差异因统计单位而异，调整参数（即 ）和迭代次数在 80 个县之间有所不同。如图 9 所示，基于县的 IAGSA 方法导致了从 MODIS 单独估算的作物分数的明显调整。尽管一些县的作物分数变化趋势相似，但它们的具体变化幅度却有所不同（图 9）。

图 8. MODIS 衍生的作物面积与县级农业统计数据的比较。图(a)、(b)和(c)分别代表玉米、水稻和大豆。

图 9. 使用 IAGSA 方法在县级调整的作物比例。第一行的图(a)表示在县级使用 IAGSA 方法后玉米、稻米和大豆的比例。第二行的图(b)表示应用 IAGSA 方法后作物比例变化的空间分布及三种作物类型的相关统计直方图。对于直方图，横轴表示变化的作物比例，间隔宽度为 ，纵轴表示像素比例（具有特定作物比例的像素数量占总作物像素数量的比例）。对于空间地图，灰色标记的列表示值保持不变的像素。

不同区域间的差距。我们还发现，较小的统计单位在调整后的作物比例中产生了更高的空间方差。与县级 IGASA 方法相比，作物比例的变化比率在各县之间空间变化（图 9），而省级 IAGSA 调整后的作物比例变化（无论是增加还是减少）的空间模式在我们研究区域的所有像素中相对均匀（图 11）。水稻比例保持不变，因为遥感推导的水稻面积与农业统计之间的省级面积差异相当小（即，Agap ）。

 5. 讨论

 5.1. 特征选择


我们使用了五个植被指数的时间序列，这些指数是从不同的光谱通道计算得出的（总共 155 个特征），以充分捕捉作物类型光谱特征在生长季节中的细微差异。与之前仅使用单一光谱指数（例如 EVI 或 LSWI）时间序列的研究相比，我们的方法利用了多个植被指数来区分不同作物类型的物候特征。这对于区分具有相似光谱特征的多种作物类型尤为重要（Peña-Barragán 等，2011）。我们的结果证实，单独使用 EVI 或 LSWI 等植被指数是次优的，而综合使用多个光谱指数是有利的。此外，关键时间段的信息对于作物类型分类也很重要（图 5）。

我们的后向特征消除策略也表明了 RF-r 模型在特征选择中的实用性。之前的研究主要使用 RF 衍生的特征重要性分数来解释不同特征对分类的贡献（Grinand 等，2013；Song 等，2017b），而我们的研究表明，RF 衍生的特征重要性也可以用于选择最佳特征，以提高映射准确性和计算效率。我们整合了后向特征消除策略和 RF 衍生的特征重要性分数，以管理原始的 155 个特征，并自动识别最佳特征。

图 10. 从 MODIS 图像中得出的作物分数直方图，带有（用颜色标记）和不带（用灰色标记）IAGSA 调整。对于直方图，横轴对应作物分数，间隔宽度为 ，纵轴表示像素数量。图（a）、（b）和（c）分别代表玉米、水稻和大豆。

图 11. 使用 IAGSA 方法在省级调整的作物比例。第一行的图(a)表示在省级使用 IAGSA 方法后玉米、稻米和大豆的比例。第二行的图(b)表示与未应用 IAGSA 方法的原始作物比例相比，玉米、稻米和大豆比例的变化。

每种作物类型的光谱-时间特征。我们发现，使用 OOB 误差选择的最佳特征与 4000 个独立测试样本的评估结果高度一致（补充材料第六部分）。这不仅确认了我们特征选择方法的可靠性，还表明本研究中使用的训练样本和测试样本具有很高的代表性。

我们的结果表明，结合向后消除特征策略的 RF-r 模型在特征选择方面是有前景的，特别是在处理大型候选特征集时，它不仅对不同特征的相对重要性进行排名，还有效地产生

最高的分类准确率和更低的计算时间（图 4）。这表明我们的方法可以应用于其他地区的土地覆盖类别映射，特别是对计算要求较高的大面积映射。


5.2. 遥感和农业统计不确定性对 IAGSA 的影响


遥感被认为是监测作物面积及其动态的客观有效方法。然而，训练样本、参考数据、图像预处理中的不确定性以及分类器的不完美学习能力可能会给遥感衍生的作物类型图带来误差（Foody, 2002; Verburg et al., 2011）。我们的研究证明了利用农业统计数据来改善遥感衍生的作物类型图和调和面积估算的价值（图 9 及补充材料的第四部分）。

高质量的农业统计数据对于提议的 IAGSA 至关重要。本研究中的农业统计数据来自《中国农村统计年鉴》，这是中国最广泛使用的农业统计资料（中国国家统计局，2019 年）。为了评估我们研究区域农业统计数据的可靠性，我们通过将无偏的作物面积估计与农业统计数据进行比较，进行了额外的测试。我们使用 Olofsson 等人（2014 年）提出的修改方法，从 MODIS 衍生的作物分布图和 HSRI 参考图生成无偏的面积估计及相关的置信区间（补充材料中的第七部分和第八部分）。我们的结果显示，超过 的目标县的农业统计数据落入所有作物类型的无偏面积估计的置信区间内（补充材料中的第八部分），这表明本研究中使用的农业统计数据的可靠性。 此外，我们通过采访当地专家，包括黑龙江省农业科学院的科学家和黑龙江省几个县（如北安市）的官员，对农业统计数据进行了定性评估。这些受访者确认了我们研究区域农业统计数据的可靠性。

然而，我们承认农业统计数据中可能存在不确定性，因为通过家庭调查和与农民的访谈获得的作物面积信息并不总是准确的（刘等，2020；佩尼亚-巴拉甘等，2011）。由于缺乏绝对的真实情况，通常很难评估农业统计数据的质量，同时也很难判断这些统计数据是否比大区域内的遥感作物面积更准确。因此，我们在此讨论以下三个案例，以说明我们 IAGSA 方法的适用性：

  • 案例 1：遥感衍生的作物估计与农业统计之间高度一致。在这种情况下，通过 IAGSA 方法调整的作物比例变化很小或没有变化（例如，水稻制图，图 8b 和图 9）。

  • 案例 2：遥感区域估计的不确定性大于农业统计数据的不确定性。在这种情况下，IAGSA 方法可以改善基于 MODIS 的估计，这通过我们使用低代表性训练样本进行的额外测试得到了证明（补充材料的第四部分）。在训练样本的质量和数量方面的代表性是一个最具挑战性的问题，常常限制基于遥感的作物类型制图的准确性（Wang et al., 2019; Zhong et al., 2014）。我们发现应用 IAGSA 方法后的制图准确性（用 和 RMSE 表示）显著高于未使用 IAGSA 的方法（补充材料的第四部分）。

  • 案例 3：遥感区域估计的不确定性小于农业统计的不确定性。在本研究中，我们提出的 IAGSA 方法通过迭代线性调整规则分配了遥感与农业统计之间的区域差距，这意味着一个县内所有像素的作物比例几乎成比例地发生了变化。因此，使用 IAGSA 方法后，遥感衍生的作物比例的整体空间分布模式并没有显著改变。如图 10 所示，应用和未应用 IAGSA 方法的遥感图像衍生的作物比例直方图彼此非常一致，这表明 IAGSA 能够在不影响全球视角下固有价值分布特征的情况下调整作物比例。在这方面，IAGSA 可以被视为一种后分类程序，其中农业统计提供了重要的补充信息，以细化基于遥感的作物类型图，而分类准确性主要由遥感决定。 此外，根据 IAGSA 的性质，作物面积的差异被空间分配到一个县内的所有像素，因此农业统计数据所带来的不确定性被所有这些像素共享。这意味着地理上较大的县（即覆盖更多像素的县）对农业统计数据的不确定性影响较小。


5.3. 县级与省级 IAGSA


县级与省级结果的差异（第 4.4 节）引发了如何选择适合 IAGSA 方法的农业统计单位的问题，特别是在进行大规模作物类型制图时。在中国，各级农业统计（即县级、地级和省级）以不同的方式和由不同的行政部门收集（中国国家统计局，2019）。小规模（即县级）收集的统计数据被认为比大规模统计数据更准确，因此推荐用于 IAGSA 方法。此外，由遥感引入的作物类型分类错误可能因景观异质性和种植模式的多样性而因地区而异。这种遥感衍生错误的空间变异性可以通过使用小规模农业统计数据的 IAGSA 方法部分解决。相反，大规模农业统计数据（例如。当省级（或更高）数据更可靠或遥感衍生的误差在空间上均匀时，IAGSA 更倾向于使用这些数据。


5.4. 与其他研究的比较


其他研究结合了遥感和农业统计数据，制作了作物类型地图（Frolking 等，2002；Monfreda 等，2008；You 等，2009）。在这些研究中，使用了卫星衍生的土地覆盖数据集（例如，MODIS 土地覆盖产品）来确定耕地的空间分布，随后开发了一个具有多个约束的统计下行模型，以根据农业统计数据在每个网格（例如，10 公里）内空间分解作物类型面积的总量（Frolking 等，2002；You 等，2014）。然而，这些统计模型很少考虑不同地区个别作物类型的光谱特征。此外，全球土地覆盖图中耕地的定义通常与农业统计记录的作物类型不一致，这也可能引入误差（See 等，2015）。相比之下，我们提出的方法可以生成作物分布图，并在遥感图像和农业统计之间提供一致的面积估算。 与上述统计模型相比，我们的方法利用了作物的光谱特征和物候信息，因此在农业统计的不确定性方面更具鲁棒性，并且使得映射结果更易于解释。

 5.5. 未来改进


尽管在本研究中我们主要关注黑龙江省的三种主要作物类型（即玉米、水稻和大豆）使用 MODIS 数据的种植，但所提出的方法可以扩展到使用其他中等或粗分辨率影像在不同地区绘制任何其他作物类型。为了提高我们方法的适用性，可以进行三项改进。首先，我们的方法可以被视为一种数据融合形式，其中作物类型图的最终准确性主要依赖于遥感的表现。因此，收集和编制高质量的训练样本，以及提供重要辅助数据的可靠耕地参考图是未来的高优先事项。其次，在我们的研究中，我们通过使用三个独立的 RF-r 模型得出了三个作物分数图，然后使用 IAGSA 方法对它们进行了单独调整。这可能会导致某些情况下，作物分数的总和超过一个像素内的耕地分数。尽管这种情况在我们的研究中仅占 的像素，但我们提醒，在同时绘制多种作物类型时需要解决这种冲突。 与此同时，随着混合像素中作物类型数量的增加，作物制图的准确性可能会下降。为了解决这个问题，一方面，多输出回归模型（Dapogny 等，2017）有潜力提高不同作物类型的遥感覆盖一致性。另一方面，平衡多个目标和约束是重要的，随着作物类型数量的增加，这些目标和约束也会增加，以改进 IAGSA 方法，同时为多种作物类型分配区域差距。最后，我们的 IAGSA 方法使用迭代线性函数来修正遥感引入的作物面积偏差。其他可以提供关于基于遥感结果的空间不确定性信息或关于作物面积概率的先验知识的数据，如农业适宜性数据（Nicole 等，2019），也可以被纳入，以建立一个非线性函数，使 IAGSA 能够分配在空间上非平稳的制图误差。

 6. 结论


本文提出了一种将粗分辨率图像与农业统计数据相结合的方法，以绘制亚像素作物类型分布。该方法包括两个部分：（1）基于遥感的作物类型制图；（2）基于统计的调整。使用该方法绘制了 2011 年中国黑龙江省的主要作物类型，即玉米、水稻和大豆。在第一部分中，我们展示了反向特征消除策略在识别不同作物类型最佳光谱时空特征方面的实用性。基于 MODIS 的作物比例与基于 HSRI 的结果良好一致，所有作物类型的 值均大于 0.75。为了进一步减少从 MODIS 图像派生的作物类型图中的潜在偏差或错误，在第二部分中，我们使用 IAGSA 方法通过迭代调整的线性规则，空间分配 MODIS 派生面积与农业统计数据之间的面积差距。调整后的亚像素作物类型图不仅与作物面积的农业统计数据一致，而且与未调整的图也显示出良好的空间一致性。 我们还发现 IAGSA 方法的性能依赖于所选的空间尺度。总体而言，我们的方法可以被视为一种数据融合策略，通过农业统计数据来优化基于遥感的结果，这使得它在遥感方法单独表现不佳时，特别适用于大面积的作物类型制图。


利益冲突声明


作者声明，他们没有已知的竞争性财务利益或个人关系，这些关系可能会影响本文所报告的工作。

 致谢


本研究得到了中国国家自然科学基金（41921001 和 41901380）、中国国家重点研发计划（2019YFA0607401）和中央高校基本科研业务费（CCNU20QN032）的资助。我们感谢黑龙江省农业科学院的研究人员提供黑龙江省的参考作物图。我们还感谢四位匿名评审的建设性意见。


附录 A. 补充数据


本文的补充数据可以在线找到，网址为 https://doi. org/10.1016/j.rse.2021.112365。

 参考文献


Atzberger, C., Rembold, F., 2013. 使用 AVHRR NDVI 时间序列和神经网络对冬季作物的空间分布进行亚像素级映射。遥感. 5, .

Belgiu, M., Drăguţ, L., 2016. 随机森林在遥感中的应用：回顾与未来方向。ISPRS 摄影测量与遥感杂志 114, 24-31。

布雷曼，L.，2001。随机森林。机器学习。5-32。

Calzadilla, A., Rehdanz, K., Betts, R., Falloon, P., Wiltshire, A., 2013. 气候变化对全球农业的影响。气候变化 120, 357-374。

Chang, J., Hansen, M.C., Pittman, K., Carroll, M., DiMiceli, C., 2007. 使用 MODIS 时间序列数据集对美国的玉米和大豆进行制图。农业杂志 99, 1654.

陈, J., 陈, J., 廖, A., 曹, X., 陈, L., 陈, X., 何, C., 韩, G., 彭, S., 陆, M., 张, W., 童, X., 米尔斯, J., 2015. 全球土地覆盖图绘制，分辨率为 30 米：一种基于 POK 的操作方法. ISPRS 摄影测量与遥感杂志, 103, 7-27.

陈, Y., 陆, D., 莫兰, E., 巴蒂斯特拉, M., 杜特拉, L.V., 桑切斯, I.D.A., 达席尔瓦, R.F.B., 黄, J., 路易斯, A.J.B., 德奥利维拉, M.A.F., 2018. 使用 MODIS 时间序列数据绘制农田、种植模式和作物类型. 国际应用地球观测与地理信息杂志 69, 133-147.

科尔, M.B., 奥古斯丁, M.A., 罗伯逊, M.J., 曼纳斯, J.M., 2018. 食品安全的科学. npj Sci. Food 2.

Dapogny, A., Bailly, K., Dubuisson, S., Ieee, 2017. 用于面部动作单元检测的多输出随机森林. (2017). 纽约. 在：第十二届 IEEE 国际自动面部和姿态识别会议，页 135-140.

董杰, 肖晓, 寇伟, 秦勇, 张国, 李丽, 金超, 周勇, 王杰, 比拉达尔, 刘杰, 摩尔, B., 2015. 通过时间序列的 Landsat 影像和基于物候的算法跟踪 1986-2010 年水稻种植面积的动态变化. 遥感与环境. 160, 99-113.

杜维耶尔，G.，德福尔尼，P.，2010。一个概念框架，用于定义利用遥感进行农业监测的空间分辨率要求。《遥感与环境》，114，2637-2650。

Erb, K.-H., Gaube, V., Krausmann, F., Plutzar, C., Bondeau, A., Haberl, H., 2007. 一套与国家统计数据一致的 2000 年全球 5 分钟分辨率土地利用数据集。土地利用科学杂志 2, 191-224。

Foerster, S., Kaden, K., Foerster, M., Itzerott, S., 2012. 使用光谱-时间剖面和物候信息进行作物类型映射。计算机电子农业。 .

福利，J.A.，拉曼库提，N.，布劳曼，K.A.，卡西迪，E.S.，杰伯，J.S.，约翰斯顿，M.，穆勒，N.D.，奥康奈尔，C.，雷，D.K.，韦斯特，P.C.，巴尔泽，C.，贝内特，E.M.，卡彭特，S.R.，希尔，J.，蒙弗雷达，C.，波拉斯基，S.，罗克斯特罗姆，J.，希汉，J.，西伯特，S.，蒂尔曼，D.，扎克斯，D.P.，2011 年。为耕种的星球提供解决方案。《自然》 。

Foody, G.M., 2002. 土地覆盖分类精度评估的现状。遥感环境. 80, 185-201.

弗里茨, S., 西, L., 麦卡勒姆, I., 你, L., 本, A., 莫尔查诺娃, E., 杜尔劳尔, M., 阿尔布雷希特, F., 希尔, C., 佩尔格, C., 哈夫利克, P., 莫斯尼尔, A., 索恩顿, P., 伍德西赫拉, U., 埃雷罗, M., 贝克-雷谢夫, I., 贾斯蒂斯, C., 汉森, M., 龚, P., 阿卜杜勒·阿齐兹, S., 奇普里亚尼, A., 库马尼, R., 切基, G., 康切达, G., 费雷拉, S., 戈麦斯, A., 哈法尼, M., 卡伊塔基雷, F., 马兰丁, J., 穆勒, R., 纽比, T., 诺古尔马, A., 奥卢塞贡, A., 奥特纳, S., 拉贾克, D.R., 罗查, J., 谢帕申科, D., 谢帕申科, M., 特列霍夫, A., 天旺, A., 瓦恩库特森, C., 文特鲁, E., 文斌, W., 范德维尔德, M., 邓伍迪, A., 克拉克斯纳, F., 奥伯斯泰因, M., 2015. 全球农田和田块大小的映射. 全球变化生物学, 21, 1980-1992.

Frolking, S., Qiu, J., Boles, S., Xiao, X., Liu, J., Zhuang, Y., Li, C., Qin, X., 2002. 结合遥感和地面统计数据开发中国稻米农业分布的新地图。全球生物地球化学循环 16。

Grinand, C., Rakotomalala, F., Gond, V., Vaudry, R., Bernoux, M., Vieilledent, G., 2013. 利用多日期 Landsat 卫星图像和随机森林分类器估算 2000 年至 2010 年马达加斯加热带湿润和干燥森林的森林砍伐情况。遥感环境. 139, 68-80.

关, H., 李, J., 查普曼, M., 邓, F., 姬, Z., 杨, X., 2013. 使用随机森林将正射影像和激光雷达数据集成用于基于对象的城市主题制图. 国际遥感杂志, 34, 5166-5186.

黑龙江省统计局，2010-2012 年。黑龙江省农村统计年鉴。
    •  通讯作者。

    电子邮件地址：wuwenbin@caas.cn (W. Wu)。