文章
基于整体一致性差异的农田产品融合方法:以中国为例
钟燕飞
, 罗畅
, 胡欣
, 魏丽飞
, 王欣宇
和 金书颖
1 国家测绘、遥感信息工程重点实验室,
武汉大学,中国武汉 430079;zhongyanfei@whu.edu.cn (Y.Z.); whu_huxin@whu.edu.cn (X.H.);wangxinyu@whu.edu.cn (X.W.); jsy@whu.edu.cn (S.J.)
湖北省自然资源遥感监测工程研究中心
武汉大学,中国武汉 430079
湖北大学资源与环境科学学院,中国武汉 430062;weilifei2508@hubu.edu.cn
* 通信: luochang@whu.edu.cn
收到: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]。基于数据一致性生成协同耕地产品的广义融合模型涉及通过分析现有遥感耕地产品的一致性差异来建立某些决策规则,从而生成协同耕地产品
。这本质上等同于对多个现有遥感分类产品的决策级融合。
广义模型的一般逻辑是,在重新采样到一定大小的网格单元后,现有遥感耕地产品中一致性更高的网格单元更有可能真实地包含耕地
。基于这一原则,广义模型可以分为以下三个步骤:(1)根据每个耕地产品的质量差异创建融合组合水平表;(2)根据建立的融合组合水平表,采用直接算术平均法融合多个现有产品;(3)根据耕地统计数据自适应地确定最佳组合水平约束。
2.2.1. 融合组合等级表的创建
由于遥感耕地产品之间的一致性越高,实际耕地面积的接近性就越大,因此可以根据所使用的所有遥感耕地产品之间的整体一致性差异将不同的一致性水平进行划分。因此,对于
个产品,有
个一致性水平。产品之间的一致性越高,网格单元获得的一致性水平就越高,网格单元存在耕地的可能性或概率也就越高。
当某些一致性水平有
种不同的产品组合时,总共有
种不同的产品组合。因此,根据不同原始产品在耕地面积估算中的准确性差异,我们对不同产品赋予不同的权重,以便对这些不同的产品组合进行排序,并制作融合组合水平表,以生成协同耕地产品。
2.2.2. 基于直接算术平均法的农田产品融合
基于直接算术平均法的遥感耕地产品融合涉及简单地计算网格单元中不同产品的平均耕地比例。然后根据融合组合水平表进行融合。
其中
是产品
的耕地比例,
是网格单元中耕地面积比例大于 0 的产品数量;也就是说,协议水平是根据产品一致性来确定的。
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. 融合组合水平表。
|
|
|
|
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]。对于各种遥感农田产品的数据特征,我们也可以考虑一种称为整体一致性度量的相似性度量,以表示单一产品与所有产品之间的差异。
在开始时,各种农田产品之间的信心距离被定义为在网格单元
中多个现有产品之间相互支持程度的测量,具体如下:
其中
是
中的元素,集合表示所有农田产品的农田比例。
的值越大,网格单元中两个产品的面积比例差异越大,这意味着
和
之间的相互支持较弱。置信距离的定义充分利用了多个现有农田产品在每个网格单元中的隐含信息,并减少了对先验信息的要求。
为了标准化所有农作物产品的耕地面积比例之间的相互支持,一致性距离
的定义如下:
其中
代表
中的最大值。如果
中的所有比例完全相等,则网格单元
中的耕地比例可以直接由任何产品
的比例表示。除了这种情况,
中的置信距离并不全为 0,这意味着
。显然,
具有以下两个特征[32]:(1)
与
成反比;(2)
。
在考虑不同农田产品的空间分布质量差异的情况下,我们让
表示
与所有
比例的相似性:
其中
是每个产品的权重系数,根据每个产品的空间精度获得。此外,
与产品
的空间定位精度的 kappa 系数正相关,即
。同时,整体一致性向量
被构建用来表示对
中所有元素的整体支持程度,其中向量
代表一致性,向量
代表权重。
表示一致性矩阵。
主要权重
与 CBIR 研究中相似性度量定义的权重不同。例如,Rui 等人 [30] 直接使用标准差度量的倒数作为权重的估计,这也是数据元素之间信心的一个测量。由于目标对象的不同,我们研究中生成的实际方程、相似性度量的函数和整体一致性
也有所不同。最重要的是,整体一致性向量
被用来分解,以获得对应于最大特征值的特征向量,最终生成不同产品的权重以用于融合过程。
最大特征值
存在于
中,且仅对应于该特征值的特征向量
中的元素均为正,且使得
。
对应于最大特征值的特征向量随后被归一化:
在这个过程中获得的
被视为
的权重,即产品
在融合过程中的耕地比例。在 CBIR 系统中,权重可以通过相关反馈以递归实现动态修改[31]。在本文中,不同网格单元中不同产品的最终权重也可以根据它们的整体一致性通过上述过程自动调整。
最后,可以按如下方式计算每个网格单元
融合后的耕地面积比例: