文章
基于整体一致性差异的农田产品融合方法：以中国为例
收到：2019 年 4 月 8 日；接受：2019 年 4 月 26 日；发表：2019 年 5 月 6 日
摘要
现有的遥感耕地产品之间存在不一致性，其估计耕地面积和空间定位的准确性需要提高。现有的生成协同耕地产品以提高现有产品准确性的方法在融合过程中并未考虑不同产品在每个网格单元中的整体一致性差异。为了减少个别耕地产品异常估计耕地面积对结果的影响，本文提出了一种基于整体一致性差异，通过融合多个现有耕地产品生成协同耕地产品的方法。在该方法中，融合多个现有耕地产品的过程是基于每个网格单元中所有耕地产品的估计耕地面积的整体一致性差异。然后，在确定与耕地统计数据的最佳组合水平后生成协同耕地产品。 在本研究中，我们将 2010 年设定为基准年，并使用所提出的方法对四种遥感耕地产品进行实验：GlobCover 2009、MODIS 耕地、MCD12Q1 和 FOM-GLC，以及国家耕地统计数据。结果表明，所提出的方法生成的协同耕地产品在耕地面积估算和空间定位的准确性上优于通过广义模型和原始产品获得的结果。
关键词：遥感；农田产品；多源数据融合；整体一致性差异
1. 引言
准确的农田分布信息对食品安全和环境可持续性非常重要[1]。卫星遥感影像为我们提供了一个高效的数据源，从中可以获得大面积的农田分布信息及其时空变化[2]。近年来，发布了多种全球和洲际遥感土地覆盖产品，其中包括农田类别，许多产品可免费向公众提供。早期的土地覆盖产品空间分辨率较低，通常为 1 公里，如 IGBP-DISCover [3]、UMD LandCover [4]和 GLC2000 [5]。随着卫星技术和分类方法的改进，土地覆盖产品的空间分辨率已提高到 500 米，如 MODIS Collection 5 [6,7]，然后提高到 300 米，如 GlobCover 2009 [8]。如今，已有几种分辨率为 30 米的农田产品可用，如 GlobeLand30 [9,10]。
来自 GLC [11]。然而，从不同的遥感土地覆盖产品中提取的农田类别之间存在相当大的不一致。此外，这些产品估计的农田面积与官方统计数据相差甚远，空间定位精度较差，因为受到混合像素的限制 [12,13]。
已有多项研究集中于通过多源数据融合提高现有遥感产品的准确性。一些研究采用了位置信息和平滑技术，并通过地理回归分析将地理数据作为训练样本与现有遥感产品进行预测，从而预测耕地的分布[14-19]。其他研究则使用了多源数据融合方法，如邓普斯特-谢弗证据理论[18]、模糊集理论[20]和贝叶斯理论[21,22]。上述方法在一定程度上提高了土地覆盖产品和耕地产品的准确性，但需要大量难以获得的先验知识。一些研究通过分析多源数据之间的一致性建立了融合决策规则，但忽视了它们的质量差异[23-27]。Fritz 等[28]、Lu 等[29]和 Chen 等[19]根据输入多源数据的准确性评估对其进行排名，并相应分配不同的权重，以生成协同耕地产品，并提供了解决上述问题的可行方案[12]。 然而，以往研究中使用的现有通用方法在融合农田产品以生成协同农田产品的过程中，并未考虑每个网格单元中不同产品之间的整体一致性差异，因此受到异常农田面积的低质量产品的影响很大，导致估计结果与实际农田面积之间存在较大差异。
在本文中，为了减少个别产品异常估计耕地面积对结果的影响，我们提出了一种通过基于整体一致性差异融合多个现有产品生成协同耕地产品的新方法。在该方法中，我们首先创建一个融合组合水平表，其中每个网格单元中所有产品的估计耕地面积的整体一致性差异被视为融合的基础。最终，在确定与耕地统计数据的最佳组合水平后生成协同耕地产品。该方法可以快速高效地获得研究区域的耕地分布结果，即协同耕地产品，具有高耕地面积估计精度和空间定位能力，而无需依赖训练样本。 在本研究进行的实验中，我们将 2010 年设定为基准年，并使用了 GlobCover 2009 [8]、MODIS Cropland [7]、MCD12Q1 [6]和 FROM-GLC [11]的遥感耕地产品，以及来自地理信息监测云平台的国家耕地统计数据，生成了中国的新协同耕地产品。该协同耕地产品在耕地面积估算的准确性和空间定位的准确性上，均对每个原始产品和基于广义模型的现有方法生成的结果进行了改进。
本文的其余部分组织如下：第二节中，我们描述了本文使用的四种农田产品和通用模型。第三节中，我们介绍了本文提出的生成协同农田产品的模型。第四节中，我们描述了使用所提方法进行的实验，并对协同农田产品进行质量评估。然后，我们将协同农田产品与通用模型获得的结果以及原始产品进行比较，以分析实验结果。最后，第五节提供了总结和讨论。
2. 背景
2.1. 融合数据源的描述
为了为中国生成最佳的协同耕地产品，在融合过程中，选择了来自地理信息监测云平台的 2010 年中国国家耕地统计数据作为非遥感数据，并使用了四种不同的遥感数据
全球或区域尺度的遥感农田产品被选为遥感数据源，来自四个土地覆盖产品：GlobCover 2009、MODIS Cropland、MCD12Q1 和 FROM-GLC。在提取每个土地覆盖产品的农田类型信息后，使用中国区域矢量地图裁剪所有产品，得到的中国所有土地覆盖产品的农田部分如图 1 所示。表 1 提供了这些遥感农田产品的详细信息。
全球覆盖 2009
c. MCD12Q1
b. MODIS 农田
d. FROM-GLC
图 1. 研究中使用的中国遥感耕地产品 (a) GlobCover 2009; (b) MODIS 耕地; (c) MCD12Q1; (d) FROM-GLC。
表 1. 研究中使用的遥感农田产品。
| 产品 | GlobCover 2009 | MODIS 农田 | MCD12Q1 | FROM-GLC |
| 时间 | 2009 | 2010 | 2010 | 2010 |
| 参考 | 欧洲航天局与伦敦大学学院 | USDA | 波士顿大学 | 清华大学 |
| 卫星 | MERIS | MODIS | MODIS | 陆地卫星 |
| 空间分辨率 | 300 米 | 250 米 | 500 米 | 30 米 |
遥感农田产品在空间分辨率、数据来源、时间覆盖、分类方法和准确性方面存在差异。因此，这些产品在空间分布上存在很大的不一致性。
GlobCover 2009 产品由欧洲航天局（ESA）和鲁汶大学（UCL）制作，主要基于中分辨率成像光谱仪
（MERIS）反射率数据，空间分辨率为 300 米。MODIS 农田产品基于 2000 年至 2008 年的 250 米中等分辨率成像光谱仪（MODIS）数据。对于该产品，使用决策树分类器和 MODIS 指标计算了农田的概率。然后，根据农田面积的统计数据，设定了一个阈值以确定农田面积。波士顿大学开发了一种新的全球土地覆盖产品，称为 MCD12Q1，使用 2000 年至 2012 年的 MODIS 数据，空间分辨率为 500 米。MODIS 1-7 波段的光谱和时间特征以及增强植被指数被用于实施决策树分类。本研究使用了 2010 年的数据。FROM-GLC 产品是世界上第一个综合高分辨率土地覆盖图，由清华大学的研究人员使用 Landsat 主题制图仪（TM）/增强主题制图仪（ETM+）图像创建。它是通过整合多种自动分类算法（随机森林、基于特征空间转换的快速聚类算法等）生成的。 尽管这些输入耕地产品在时间覆盖上略有不同，但它们之间的不一致性对最终协同产品的影响远大于时间差异的影响。
2.2. 生成协同农田产品的广义模型
我们研究的基础是多源遥感数据一致性融合的广义模型[28]，也称为改进模糊一致性评分（MFAS）[19]或分层优化协同方法（HPSA）[29]。基于数据一致性生成协同耕地产品的广义融合模型涉及通过分析现有遥感耕地产品的一致性差异来建立某些决策规则，从而生成协同耕地产品
广义模型的一般逻辑是，在重新采样到一定大小的网格单元后，现有遥感耕地产品中一致性更高的网格单元更有可能真实地包含耕地
2.2.1. 融合组合等级表的创建
由于遥感耕地产品之间的一致性越高，实际耕地面积的接近性就越大，因此可以根据所使用的所有遥感耕地产品之间的整体一致性差异将不同的一致性水平进行划分。因此，对于
当某些一致性水平有
2.2.2. 基于直接算术平均法的农田产品融合
基于直接算术平均法的遥感耕地产品融合涉及简单地计算网格单元中不同产品的平均耕地比例。然后根据融合组合水平表进行融合。
每个网格单元中耕地面积的比例
其中
2.2.3. 通过耕地统计确定最佳组合水平
在根据不同融合级别对每个网格单元进行分层集成后，获得了不同级别的融合结果，这意味着不同的网格单元在其对应的融合级别上进行了融合。当根据不同的融合级别完成融合后，融合结果在不同的组合级别上累积，直到结果最接近耕地面积的统计数据，从而获得最佳组合级别。然后可以生成协同耕地产品。
通过许多研究者的尝试，广义模型可以在一定程度上生成一种协同耕地产品，其耕地面积估算和空间定位的准确性高于原始产品。然而，由于不同组合水平的网格单元的融合耕地比例是基于不同产品的平均耕地比例计算的，广义模型仅适用于整体一致性较高的网格单元。因此，它受到异常产品耕地面积比例的很大影响，这可能导致计算结果与实际耕地面积之间存在较大差异。
3. 基于整体一致性差异融合农田产品生成协同农田产品的模型
为了减少个别异常产品对耕地面积整体比例的影响，我们在广义模型中改进了基于直接算术平均的融合算法，并创新性地提出了一种通过耕地产品融合生成协同产品的模型，基于整体一致性差异。
由于广义模型中使用的原则是，现有产品一致性较高的网格单元更可能包含农田，因此该模型首先需要创建融合组合水平表。接下来，根据整体一致性差异融合原始农田产品。在确定与农田统计数据的最佳组合水平后，生成协同农田产品。
有必要验证模型生成的协同耕地产品的准确性，并将其与广义模型的结果以及原始产品进行比较，以评估模型。所提模型的具体流程图如图 2 所示。
图 2. 提议模型的流程图。
3.1. 农田产品的预处理
在不同遥感耕地产品融合之前，所有产品都统一到相同的地理坐标系统（WGS-1984），并进行了地理注册。此外，原始耕地产品需要根据通用方法的实践进行标准化和重采样，以达到所有原始耕地产品的最小分辨率[28,29]。在将所有产品重采样到
3.2. 融合组合等级表的创建
融合组合水平表可以根据原始产品之间更大一致性的原则创建，这意味着更有可能接近真实的耕地。本文中引入的产品数量为四，因此根据不同网格单元中耕地的一致性，有四个不同的一致性水平。一致性水平的空间分布如图 3 所示。
质量评估显示，GlobCover 2009、MODIS 耕地、MOD12Q1 和 FROM-GLC 在耕地面积估算中的准确性依次降低。最终的融合组合水平表可以根据每个产品的质量差异创建，如表 2 所示（其中“0”表示非耕地，“1”表示包含耕地）。
完成这一部分后，每个网格单元将被划分并标记为不同的一致性水平和融合组合水平。一致性水平高的网格单元是该模型下一步的重点。
图 3. 一致性水平的空间分布。
表 2. 融合组合水平表。
一致性级别 |
组合级别 |
全球覆盖 2009 |
MODIS 农田 |
MOD12Q1 | FROM-GLC |
| I | 1 | 1 | 1 | 1 | 1 |
| II | 2 | 1 | 1 | 1 | 0 |
| 3 | 1 | 1 | 0 | 1 | |
| 4 | 1 | 0 | 1 | 1 | |
| 5 | 0 | 1 | 1 | 1 | |
| 6 | 1 | 1 | 0 | 0 | |
| 7 | 1 | 0 | 1 | 0 | |
| 8 | 1 | 0 | 0 | 1 | |
| 9 | 0 | 1 | 1 | 0 | |
| 10 | 0 | 1 | 0 | 1 | |
| 11 | 0 | 0 | 1 | 1 | |
| 13 | 1 | 0 | 0 | 0 | |
| 14 | 0 | 1 | 0 | 0 | |
| 15 | 0 | 0 | 1 | 0 | |
| 13 | 0 | 0 | 0 | 1 |
3.3. 基于整体一致性差异的遥感耕地产品融合
对于一致性水平较高的网格单元，这意味着网格单元内的大多数农田产品都有值，它们的农田面积比例可能差异很大。具体来说，网格单元
因此，在融合过程中，这是该模型最重要的部分，每个网格单元中产品的耕地面积比例的权重是根据整体一致性的差异来确定的。与整体耕地面积比例一致性较高的产品赋予更高的权重，而与整体一致性较低的产品则赋予较小的权重。根据上述原则，建立了每个网格单元的耕地面积比例。
在基于内容的遥感图像检索（CBIR）中，相似性度量用于指示待检索图像的特征与目标图像特征之间的差异[30,31]。对于各种遥感农田产品的数据特征，我们也可以考虑一种称为整体一致性度量的相似性度量，以表示单一产品与所有产品之间的差异。
在开始时，各种农田产品之间的信心距离被定义为在网格单元
其中
为了标准化所有农作物产品的耕地面积比例之间的相互支持，一致性距离
其中
在考虑不同农田产品的空间分布质量差异的情况下，我们让
其中
表示一致性矩阵。
主要权重
最大特征值
对应于最大特征值的特征向量随后被归一化：
在这个过程中获得的
最后，可以按如下方式计算每个网格单元
这一重要步骤可以减少个别异常产品的影响，这些产品的耕地面积比例与整体产品的耕地面积比例一致性较差。然而，这种方法对耕地产品的数量和质量有一定要求。
通过这个过程，高一致性网格单元的耕地面积比例将很少受到个别产品异常估计耕地面积的影响，更接近真实值，而低一致性网格单元的耕地面积比例将类似于基于直接算术平均的融合结果。
3.4. 最佳组合水平的确定
不同层次的融合结果可以通过基于整体一致性差异融合不同的原始遥感农田产品获得，这意味着在所谓的分层集成后，不同的网格单元在其对应的融合层次上被融合，如图 4 所示。
图 4. 不同组合水平的融合结果。
根据已建立的融合组合层级顺序，融合结果按组合层级叠加，通过使用耕地统计数据作为自适应约束，获得最佳融合结果，如图 5 所示。当这个过程完成后，可以生成协同耕地产品。
图 5. 使用耕地统计数据获得最佳组合水平。
在使用耕地统计数据作为分层叠加的约束条件，通过组合水平获得协同耕地产品的步骤中，当使用耕地统计数据自适应地确定结果时，应根据每个网格单元中耕地面积的比例在融合后更接近真实值的范围进行相应调整，从而实现耕地的真实分布结果。
所提模型将每个网格单元中所有产品估计的耕地面积的整体一致性差异作为融合的基础，并减少个别产品异常估计耕地面积对结果的影响。因此，研究区域的耕地分布结果，即协同耕地产品，能够快速高效地获得高精度的耕地面积估计以及空间定位，而无需依赖训练样本。
4. 实验与分析
4.1. 精度测量
4.1.1. 耕地面积估算的准确性
通过使用耕地统计数据来检验遥感耕地产品的估计耕地面积，从而获得耕地面积估算的准确性。
该省
其中
最后，为了反映每种产品估计耕地面积与耕地统计数据的巧合，计算每种产品耕地面积比例与统计数据耕地面积比例的相关系数，如公式（11）所示：
其中
4.1.2. 空间定位的准确性
我们使用混淆矩阵来评估使用选定测试样本的空间定位准确性。由于研究的类别仅包括耕地和非耕地，混淆矩阵如下：
其中
对于实验结果，分类准确率是基于整体分类准确率和卡帕系数进行评估的，农田的错误率和农田的泄漏率作为辅助评估标准 [34]。
整体准确率（OA）：OA 表示正确样本与测试样本总数的比率，指示产品类别是否与真实的地面真相数据相同[35]。公式如下：
Kappa 系数（kappa）：Kappa 系数是一种离散计算方法，是通过观察混淆矩阵中的主对角线以及行和列的总数来获得概率一致性的统计量[36]：
此外，误分类率（MR）和泄漏率（LR）可以按如下方式计算：
4.2. 实验结果
根据第 2 和第 3 节中描述的方法，使用四种原始耕地产品和耕地统计数据生成协同耕地产品。上述所有地理分析和处理，以及不同地图的展示，均由 ESRI ArcGIS Desktop 软件支持。融合算法的结果和确定与耕地统计数据最佳组合水平的过程是通过 MathWorks MATLAB 软件实现的。通过基于直接算术平均法（MDAA）进行融合的模型获得的协同耕地产品结果，以及通过基于整体一致性差异（MOCD）进行融合的模型获得的协同耕地产品结果如图 6 所示。
C
d
d
图 6. (a) 通过基于直接算术平均法（MDAA）融合生成的协同耕地产品。 (b) 通过基于整体一致性差异（MOCD）融合生成的协同耕地产品。 (c) 测试样本的空间分布。 (d) (a) 和 (b) 的局部放大图。
我们然后使用第 4.1 节中描述的准确性测量来评估获得的实验结果的质量，并将其与原始产品进行比较。至于空间定位精度，现有的中国地区耕地验证数据集仍然难以满足耕地准确性评估的定量要求。因此，我们使用高分辨率遥感图像获取耕地验证的测试样本，并将这些测试样本作为基准，以比较不同产品的空间定位精度。在本实验中，我们选择的测试样本部分是基于四个耕地产品的一致性水平的空间分布图，并随机抽取
4.2.1. 耕地面积估算的准确性
图 7 显示了统计的耕地面积比例的散点图，以及通过直接算术平均法（MDAA）和基于整体一致性差异法（MOCD）融合的协同耕地产品估算的耕地面积。图中的虚线表示 1:1 线。可以看出，MDAA 和 MOCD 估算的耕地面积在所有省份都接近耕地统计数据，因为所有点都分布在 1:1 线附近，这远远优于四个原始产品。
图 7. 来自统计数据的耕地面积散点图与以下估算结果的比较：(a) GlobCover 2009；(b) MODIS 耕地；(c) MCD12Q1；(d) FROM-GLC；(e) MDAA 生成的协同耕地产品；(f) MOCD 生成的协同耕地产品。
至于 MDAA，最大面积比例差异不超过
至于 MOCD，最大面积比例差异仅为
4.2.2. 空间定位的准确性
对于空间定位精度，实验结果使用耕地和非耕地测试样本进行比较。获得的精度指标如表 3 所示，这些产品在不同地理区域的空间定位总体精度（OA）如图 8 所示。
表 3. 原始和协同耕地产品的准确性指标。
| GlobCover 2009 | MODIS 农田 | MCD12Q1 | FROM-GLC | MDAA | MOCD | |
| OA |
|
|
|
|
|
|
| 卡帕 | 0.52 | 0.38 | 0.47 | 0.58 | 0.77 | 0.83 |
| MR |
|
|
|
|
|
|
| LR |
|
|
|
|
|
|
图 8. 不同地理区域农田产品的空间定位总体准确性。
从表 3 和图 8 可以看出，MDAA 和 MOCD 的空间定位精度评估指标与原始耕地产品相比显著提高。
至于 MDAA，其 OA 达到
至于 MOCD，其 OA 高达
4.3. 结果分析
从实验结果可以看出，通过基于整体一致性差异（MOCD）融合生成的协同耕地产品在耕地面积的估计精度和空间定位精度上，相较于通过直接算术平均法（MDAA）融合生成的协同耕地产品以及原始产品都有所提高。
这是因为 MOCD 不仅考虑了不同农田产品质量的差异来建立融合组合水平表，这是协同农田产品更高准确度的重要原因之一，而且还根据每个网格单元农田面积比例的整体一致性差异，确定每个产品的农田面积是否一致
在使用耕地统计数据自适应确定结果时，应相应调整范围，从而实现更接近真实耕地面积的结果。同时，在这个过程中，可以充分利用不同产品之间的互补性，以获得更接近真实耕地空间分布的协同耕地产品。更好的融合。
来自上一步的结果为这个过程提供了更好的数据，使这个过程在使用耕地统计数据以获得更好的协同耕地产品时有更大的调整空间。
通过基于整体一致性差异的融合，以及确定最佳组合水平的过程，MOCD 能够轻松获得一种协同耕地产品，其耕地面积估算的准确性以及空间定位优于 MDAA 和原始耕地产品生成的结果。
5. 结论
通过基于整体一致性差异融合多个现有耕地产品的方法，提出了一种生成协同耕地产品的方法，以获得更高质量的产品，因为现有的通用方法容易受到异常耕地面积的低质量产品的影响。在该方法中，我们首先创建一个融合组合水平表。然后，融合多个现有产品的过程考虑每个网格单元中所有产品的估计耕地面积的整体一致性差异作为融合的基础。最终，在确定与耕地统计数据的最佳组合水平后生成协同耕地产品。该方法可以减少个别异常产品对结果的影响。因此，研究区域的耕地分布结果，即协同耕地产品，可以快速高效地获得，具有高耕地面积估计精度和空间定位，而无需依赖训练样本。
实验通过使用四种遥感产品进行测试：GlobCover 2009、MODIS 农田、MCD12Q1 和 FROM-GLC，以及农田统计数据。在对结果和原始产品的准确性进行测量后，我们得出结论，通过基于整体一致性差异（MOCD）融合多个现有产品获得的协同农田产品在农田估计面积准确性和空间定位准确性上优于通过广义模型获得的结果以及原始产品。
在未来，将收集更多的遥感产品作为遥感数据源，以进一步提高生成的协同耕地产品的准确性，特别是在景观异质性区域。当产品数量增加时，基于整体一致性差异的融合算法将更具优势。同时，融合组合水平表将有更多层次，这将使得得出更接近耕地统计的结果变得更加容易。未来，我们将继续探索减少对先验知识（如耕地统计）依赖的新方法，并考虑小时间差对耕地面积变化的影响。
作者贡献：所有作者对该工作做出了重要贡献。Y.Z.、C.L. 和 X.H. 设计了研究并分析了结果。L.W.、X.W. 和 S.J. 为论文的准备和修订提供了建议。
致谢：本研究得到了中国国家重点研发计划（项目编号：2017YFB0504202）、国家自然科学基金（项目编号：41622107 和 41771385）以及湖北省自然科学基金（项目编号：2016CFA029）的支持。
利益冲突：作者声明没有利益冲突。
参考文献
刘, J.; 刘, M.; 田, H.; 庄, D.; 张, Z.; 张, W.; 唐, X.; 邓, X. 1990-2000 年中国耕地的时空模式：基于 Landsat TM 数据的分析。遥感环境. 2005, 98, 442-456. [CrossRef]
羊，D.W.; 布朗，R.B. 精准农业：作物中杂草的遥感与制图。农业工程研究杂志 2002，78，117-125。[CrossRef]
洛夫兰, T.R.; 里德, B.C.; 布朗, J.F.; 奥伦, D.O.; 朱, Z.; 杨, L. 全球土地覆盖特征数据库的开发及 IGBP 从 1 公里 AVHRR 数据中发现的内容. 国际遥感杂志. 2000, 21, 1303-1330. [CrossRef]