文章
基于超像素和多源土地覆盖产品融合方法及统计提取增强分辨率和提高精度
齐金
1
1
^(1) { }^{1}
2
,
∗
(
D
)
2
,
∗
(
D
)
^(2,**(D)) { }^{2, *(D)}
1
1
^(1) { }^{1} 1 地球探测科学与技术学院，吉林大学，长春 130026，中国；jinqi19@mails.jlu.edu.cn (Q.J.); zxq@jlu.edu.cn (X.Z.) 2 土地表面格局与模拟重点实验室，中国科学院地理科学与资源研究所，北京100101，中国 * 通讯作者：xueq@igsnrr.ac.cn
学术编辑：Onisimo Mutanga 和 Lalit Kumar
收到：2022年1月31日 接受：2022年3月19日 发布：2022年3月31日 出版商声明：MDPI 对已发表的地图中的领土主张和机构隶属关系保持中立。
摘要 现有土地覆盖数据中的差异相对较高，表明局部精度低且应用受限。多源数据融合是解决此问题的有效方法；然而，融合过程通常需要重采样以统一空间分辨率，导致空间分辨率降低。为解决这个问题，本研究提出了一种基于精细分辨率训练样本过滤和产品校正的多源产品融合制图方法。基于超像素算法、主成分分析（PCA）和统计提取技术，结合谷歌地球引擎（GEE）平台，获取了可靠的土地覆盖数据。GEE 和机器学习算法将多个产品的不可靠信息校正为新的土地覆盖融合结果。与从现有产品中提取一致像素的常用方法相比，我们提出的方法有效去除了近
38.75
%
38.75
%
38.75% 38.75 \%
85.80
%
85.80
%
85.80% 85.80 \%
11.75
−
24.17
%
11.75
−
24.17
%
11.75-24.17% 11.75-24.17 \%
2.99
%
2.99
%
2.99% 2.99 \%
20.71
%
20.71
%
20.71% 20.71 \%
关键词：土地覆盖；多源信息融合；Google Earth Engine；超像素；大面积遥感产品
1. 引言 精确的大面积土地覆盖地图是探索自然和生物活动与空间模式之间关系的基本数据支持[1,2]；生态和环境变化的模拟、监测和评估[3]；人类社会经济发育[4,5]；以及其他科学研究。粮食安全和耕地面积评估、森林变化监测、城市范围扩张和结构分析、水体面积提取和污染评估也需要及时更新大主题地图，为人类可持续发展战略提供重要指标[6,7]。随着遥感技术的发展和各种卫星数据源的出现，遥感已成为大面积绘制土地覆盖的重要方法[8,9]。
基于 NOAA/AVHRR、MODIS、ENXISAT/MERIS 和其他卫星传感器的图像数据，土地覆盖产品通常具有粗糙的空间分辨率（300-1000 米）。例如，我们拥有来自波士顿的 500 米空间分辨率的 MODIS 全球土地覆盖数据
大学[8]。欧洲航天局（ESA）的 ESA-CCI 数据集具有 300 米的空间分辨率[10]。哥白尼全球土地服务（CGLS）数据集具有 100 米的空间分辨率[11]。先前研究表明，较低的空间分辨率通常会导致较低的精度
[
12
,
13
]
[
12
,
13
]
[12,13] [12,13]
[
16
,
20
]
[
16
,
20
]
[16,20] [16,20] 数据融合可以通过整合多源数据来克服单一数据源相关的有限精度和不确定性[21,22]。一些融合决策方法，如贝叶斯理论、Dempster-Shafer 证据理论和模糊集理论，已在众多研究中得到有效应用[20,23,24]。一些研究引入了多源统计来校准融合产品，这可以提高融合结果的精度。原始产品的精度显著影响融合结果[25,26]，因此，随着输入产品的增加，产品的权重是根据先验知识确定的；否则，很难获得好的结果[27]。由于融合决策方法存在上述局限性，研究人员从先前的地表覆盖产品中选取样本来更新地图。例如，美国地质调查局（USGS）提出了光谱变化监测方法来识别未变化区域，并使用这些区域作为样本来训练决策树分类器，该分类器可以快速更新土地覆盖地图[28]。 然而，仅从单一地表覆盖产品来源获得的信息往往比融合多个产品 [21] 的可靠性要低。
多源数据融合通常需要重采样以统一空间分辨率并获得一致区域 [16]。一致区域被定义为在相同地理位置上，各种土地覆盖产品仍保持在同一类别的区域 [29]。我们通常对保留的每个地面覆盖产品信息有很高的信心 [27]。因此，从多源产品的一致区域中提取有效信息并校正不一致区域被认为是有效提高制图精度的融合方法 [30]。然而，由于重采样，最终结果仍然具有粗糙的空间分辨率。
为了解决上述问题，本研究基于现有方法[31]进一步提出了一种基于超像素和统计提取的多源土地覆盖产品融合方法，该方法实现了精细空间分辨率和高精度融合结果。第一步，我们分析多个产品的数据一致性，将研究区域划分为 300 米空间分辨率（粗）一致区域和不一致区域，并在 GEE 平台上构建 30 米特征图像层。第二步，使用 PCA 和 SNIC（简单非迭代聚类）算法对粗一致区域内的图像进行分割，去除异常像素，获得 30 米空间分辨率（细）一致区域。第三步，通过局部自适应采样方法获得可靠的训练样本，并通过机器学习校正不一致区域，以生成精细空间分辨率和高精度融合结果。然后，我们以东南亚为例应用和验证了所提出的方法。
2. 材料 2.1. 研究区域 研究区域位于
30
∘
N
30
∘
N
30^(@)N 30^{\circ} \mathrm{N}
90
∘
E
90
∘
E
90^(@)E 90^{\circ} \mathrm{E}
120
∘
E
120
∘
E
120^(@)E 120^{\circ} \mathrm{E}
km
2
km
2
km^(2) \mathrm{km}^{2} 图 1. 研究区域和多源土地覆盖产品一致性分析结果。(a) 研究区域地图。数据来自 NASA 地球观测站网站（https:/ / neo.gsfc.nasa.gov （2021 年 10 月 30 日访问）），用于 2015 年净初级生产力产品，以及将复合月度产品转换为年度产品；(b) 300 米空间分辨率（粗）多源产品完全一致性区域，其中一致性区域被定义为在相同地理位置各种土地覆盖产品保持在同一类别中的区域。研究区域被划分为
5
∘
×
5
∘
5
∘
×
5
∘
5^(@)xx5^(@) 5^{\circ} \times 5^{\circ}
b
b
b b 位置差异会导致不同地物的光谱反射率出现巨大变化，对于受气候影响的地物（例如耕地和草地）这种变化更为明显[34]。此外，局部区域内的地物分布与整个研究区域相比存在显著差异，这种差异在选择样本时会对精度产生影响[35]。同时，每次解释的区域面积必须减少以确保 GEE 的运行内存不被超出。因此，本研究将研究区域划分为 31 个
5
∘
×
5
∘
5
∘
×
5
∘
5^(@)xx5^(@) 5^{\circ} \times 5^{\circ} 2.2. 数据来源与预处理 本研究使用了 2015 年覆盖研究区域的八个产品：中分辨率成像光谱仪土地覆盖类型产品（MCD12Q1）、CCI-LC、哥白尼全球土地服务（CGLS）、FROM-GLC、GFSAD30、PALSAR、全球地表水数据（GSWD）和全球人类居住层建成区（GHS-BUILT）。产品详细信息列于表 1。我们将这八个产品统一到同一地理坐标系（WGS84，1984 年世界大地测量系统）以及空间分辨率统一，以促进多源数据融合。在先前的研究中，数据空间分辨率统一通常是通过选择所有产品的粗空间分辨率作为重采样的目标来实现的，以最佳地控制产品精度。考虑到这八个产品中只有一个是 500 米分辨率的（MCD12Q1），其次是 300 米分辨率，重采样到 300 米可以保留最大量的细节。最近邻是用于重采样的技术，这是通过参考先前研究确定的。 由于不同的生产者和传感器，这八种产品的分类系统定义有所不同，我们将每种产品的分类系统进行整合，最终将其分类为九种地面覆盖类别之一（图 1b）。有关分类系统的更多详细信息，请参见文献和 S 1 [31]。
表1。八种土地覆盖产品的详细信息。
产品名称 来源 空间分辨率（米） 传感器 分类方法
MCD12Q1
波士顿大学 500
MODIIS
决策树和神经网络
CCI-LC
ESA
300
MERIS FR/RR, AVHRR, SPOTVGT, PROBA-V
非监督分类和机器学习
CGLS
ECJRC
100
PROBA-V
随机森林
FROM-GLC
清华大学，中国 30
TM ETM +
支持向量机，随机森林
GFSAD30
USGS
30
MODIS
机器学习
PALSAR
JAXA
25
PALSAR
监督分类
GSWD
ECJRC
30
TM ETM + OLI
监督分类
GHS-BUILT
ECJRC
30
TM ETM + OLI
机器学习
Product Name Source Spatial Resolution (m) Sensor Classification Method
MCD12Q1 Boston University 500 MODIIS Decision tree and neural network
CCI-LC ESA 300 MERIS FR/RR, AVHRR, SPOTVGT, PROBA-V Unsupervised Classification and Machine learning
CGLS ECJRC 100 PROBA-V Random forest
FROM-GLC Tsinghua University, China 30 TM ETM + Support vector machine, random forest
GFSAD30 USGS 30 MODIS Machine learning
PALSAR JAXA 25 PALSAR Supervised classification
GSWD ECJRC 30 TM ETM + OLI Supervised classification
GHS-BUILT ECJRC 30 TM ETM + OLI Machine learning | Product Name | Source | Spatial Resolution (m) | Sensor | Classification Method |
| :--- | :--- | :--- | :--- | :--- |
| MCD12Q1 | Boston University | 500 | MODIIS | Decision tree and neural network |
| CCI-LC | ESA | 300 | MERIS FR/RR, AVHRR, SPOTVGT, PROBA-V | Unsupervised Classification and Machine learning |
| CGLS | ECJRC | 100 | PROBA-V | Random forest |
| FROM-GLC | Tsinghua University, China | 30 | TM ETM + | Support vector machine, random forest |
| GFSAD30 | USGS | 30 | MODIS | Machine learning |
| PALSAR | JAXA | 25 | PALSAR | Supervised classification |
| GSWD | ECJRC | 30 | TM ETM + OLI | Supervised classification |
| GHS-BUILT | ECJRC | 30 | TM ETM + OLI | Machine learning |
一致性水平被定义为在图像位置与某一类别一致的陆地覆盖产品数量。较高的水平代表对相应类别的更高置信度[27]。本研究中一致性区域生成方法与刘和徐[31]的方法相同。每个类别的最高一致性水平像素被选为该类别的区域一致区域 (图 1b)。其余部分称为不一致区域（图 1c），约占 55.3%。请注意，这里的一致性区域是 300 米空间分辨率，我们称之为粗一致性区域。因为 30 米空间分辨率的陆地覆盖产品包含更多细节，粗一致性区域忽略了这些细节，这限制了其应用。同时，错误标记的像素将提供错误的样本，限制了采样精度。
2.2.2. Landsat 数据组成 本文选取了 30 米空间分辨率的 Landsat 影像来重新解释不一致区域。在 GEE 平台上对影像进行无云像素筛选和云掩膜处理显示，2015 年的影像不足以覆盖整个研究区域，因此选取了 2015 年相邻年份（2014 年和 2016 年）的 Landsat 影像来填补数据空白。筛选了覆盖研究区域 2014-2016 三年的地球观测卫星 7 ETM+和 8 OLI 传感器的正射校正地表反射率数据，共计 27,223 张影像，包括 11,278 张 Landsat 7 ETM+和 15,945 张 Landsat 8 OLI。利用 GEE 平台上的 FMASK 算法对每张影像中因云、云阴影和雪导致的无效观测进行掩膜处理
[
36
,
37
]
[
36
,
37
]
[36,37] [36,37] Landsat 7 ETM+和 Landsat 8 OLI 使用了六个光谱波段，即蓝光、绿光、红光、近红外（NIR）、短波红外（SWIR1）和 SWIR2，并计算了 11 个光谱指数，每个光谱指数的公式如表 2 所示 [39]。参考先前研究以确保足够的物候信息和图像质量，将一年分为三个时期（时期 1：2015 年 1-120 天，时期 2：2015 年 121-240 天，时期 3：2015 年 241-365 天）[31,35]。我们使用每个时期的时序中值像素作为该时期的光谱特征，因为时序中的中值对物候变化不敏感 [40]。为了增加特征的区分度，我们还计算了每个特征三年的标准差，从而在每个瓦片上通过 GEE 平台合成了 68 个波段特征层（
(
6
(
6
(6 (6
×
3
×
3
xx3 \times 3
+
(
17
)
+
(
17
)
+(17) +(17) 表2. 光谱指数的公式。
光谱指数 公式
归一化植被指数（NDVI）[41] 近红外-红光 近红外+红光 (1)
绿色叶绿素植被指数（GCVI）[42]
NIR
Green
NIR
Green
(" NIR ")/(" Green ") \frac{\text { NIR }}{\text { Green }} (2)
增强型植被指数 (EVI) [43]
2.5
×
NIR-Red
NIR
+
6
×
Red
−
7.5
×
Blue
+
1
2.5
×
NIR-Red
NIR
+
6
×
Red
−
7.5
×
Blue
+
1
2.5 xx(" NIR-Red ")/(" NIR "+6xx" Red "-7.5 xx" Blue "+1) 2.5 \times \frac{\text { NIR-Red }}{\text { NIR }+6 \times \text { Red }-7.5 \times \text { Blue }+1} (3)
标准化燃烧指数 (NBR) [44] 近红外-短波红外 近红外+短波红外 (4)
标准化差异水指数 [45] 绿色-近红外 绿色+近红外 (5)
归一化建置指数 [46] 短波红外1-近红外 短波红外1+近红外 (6)
归一化雪指数 (NDSI) [47] 绿色-SWIR1 绿色+SWIR1 (7)
改进型土壤调整植被指数 (MSAVI) [48]
2
×
N
I
R
+
1
−
(
2
×
N
I
R
)
2
−
8
×
(
NIR
−
Red
)
2
2
×
N
I
R
+
1
−
(
2
×
N
I
R
)
2
−
8
×
(
NIR
−
Red
)
2
(2xx NIR+1-sqrt((2xx NIR)^(2)-8xx(" NIR "-" Red ")))/(2) \frac{2 \times N I R+1-\sqrt{(2 \times N I R)^{2}-8 \times(\text { NIR }- \text { Red })}}{2} (8)
土壤调整总植被指数 (SATVI) [49]
1.5
×
SWIR1-Red
SWIR1
+
Red
+
0.5
−
SWIR1
2
1.5
×
SWIR1-Red
SWIR1
+
Red
+
0.5
−
SWIR1
2
1.5 xx(" SWIR1-Red ")/(" SWIR1 "+" Red "+0.5)-(" SWIR1 ")/(2) 1.5 \times \frac{\text { SWIR1-Red }}{\text { SWIR1 }+ \text { Red }+0.5}-\frac{\text { SWIR1 }}{2} (9)
裸土指数 (BSI) [50]
(
SWIR2
+
Red
)
−
(
NIR
+
Blue
)
(
WIR2
2
Red
)
+
(
NIR
+
Blue
)
(
SWIR2
+
Red
)
−
(
NIR
+
Blue
)
(
WIR2
2
Red
)
+
(
NIR
+
Blue
)
((" SWIR2 "+" Red ")-(" NIR "+" Blue "))/((" WIR2 "2" Red ")+(" NIR "+" Blue ")) \frac{(\text { SWIR2 }+ \text { Red })-(\text { NIR }+ \text { Blue })}{(\text { WIR2 } 2 \text { Red })+(\text { NIR }+ \text { Blue })} (10)
蓝-红 (BR) [51] 蓝 - 红 (11)
Spectral Index Formula
Normalized Difference Vegetation Index (NDVI) [41] NIR-Red NIR+Red (1)
Green Chlorophyll Vegetation Index (GCVI) [42] (" NIR ")/(" Green ") - 1 (2)
Enhanced Vegetation Index (EVI) [43] 2.5 xx(" NIR-Red ")/(" NIR "+6xx" Red "-7.5 xx" Blue "+1) (3)
Normalized Burn Index (NBR) [44] NIR-SWIR NIR+SWIR (4)
Normalized Difference Water Index [45] Green-NIR Green+NIR (5)
Normalized Difference Built-up Index [46] SWIR1-NIR SWIR1+NIR (6)
Normalized Difference Snow Index (NDSI) [47] Green-SWIR1 Green+SWIR1 (7)
Modified Soil-Adjusted Vegetation Index (MSAVI) [48] (2xx NIR+1-sqrt((2xx NIR)^(2)-8xx(" NIR "-" Red ")))/(2) (8)
Soil-Adjusted Total Vegetation Index (SATVI) [49] 1.5 xx(" SWIR1-Red ")/(" SWIR1 "+" Red "+0.5)-(" SWIR1 ")/(2) (9)
Bare Soil Index (BSI) [50] ((" SWIR2 "+" Red ")-(" NIR "+" Blue "))/((" WIR2 "2" Red ")+(" NIR "+" Blue ")) (10)
Blue-Red (BR) [51] Blue - Red (11) | Spectral Index | Formula | |
| :--- | :--- | :--- |
| Normalized Difference Vegetation Index (NDVI) [41] | NIR-Red NIR+Red | (1) |
| Green Chlorophyll Vegetation Index (GCVI) [42] | $\frac{\text { NIR }}{\text { Green }}$ - 1 | (2) |
| Enhanced Vegetation Index (EVI) [43] | $2.5 \times \frac{\text { NIR-Red }}{\text { NIR }+6 \times \text { Red }-7.5 \times \text { Blue }+1}$ | (3) |
| Normalized Burn Index (NBR) [44] | NIR-SWIR NIR+SWIR | (4) |
| Normalized Difference Water Index [45] | Green-NIR Green+NIR | (5) |
| Normalized Difference Built-up Index [46] | SWIR1-NIR SWIR1+NIR | (6) |
| Normalized Difference Snow Index (NDSI) [47] | Green-SWIR1 Green+SWIR1 | (7) |
| Modified Soil-Adjusted Vegetation Index (MSAVI) [48] | $\frac{2 \times N I R+1-\sqrt{(2 \times N I R)^{2}-8 \times(\text { NIR }- \text { Red })}}{2}$ | (8) |
| Soil-Adjusted Total Vegetation Index (SATVI) [49] | $1.5 \times \frac{\text { SWIR1-Red }}{\text { SWIR1 }+ \text { Red }+0.5}-\frac{\text { SWIR1 }}{2}$ | (9) |
| Bare Soil Index (BSI) [50] | $\frac{(\text { SWIR2 }+ \text { Red })-(\text { NIR }+ \text { Blue })}{(\text { WIR2 } 2 \text { Red })+(\text { NIR }+ \text { Blue })}$ | (10) |
| Blue-Red (BR) [51] | Blue - Red | (11) |
3. 方法 图 2 显示了本研究用于实现多源土地覆盖产品融合方法的过程流程图。该方法分为三个主要部分。首先，在 GEE 上检索研究区域的 Landsat ETM 和 OLI 图像，然后进行预处理以创建复合图像层。其次，使用 PCA 技术对 68 个波段进行降维，使用 SNIC 算法对粗一致区域进行图像层分割，并通过统计方法去除异常像素。最后，通过局部自适应采样方法在 30 米空间分辨率一致区域获得大量可靠的训练样本，我们重新解释不一致区域并评估精度。通过高空间分辨率图像的样本点对结果进行验证，并通过 Google Earth 进行目视解译。 图 2。流程图概述。PCA—主成分分析[51]；SNIC—简单非迭代聚类[52]。
3.1. 粗一致区域的主成分分析 大量研究表明，通过光谱特征的异常分布来确定每个像素是否已从现有的土地利用覆盖产品中发生变化[28,53]。本研究将该方法应用于净化粗一致区域（图 1b）中的样本。用于统计分布的光谱特征通常由多个指标组成，例如张等人使用距离 39 个指标均值最小的像素作为样本[35]，因为不存在能够区分多个特征的单一指标，并且对于多个指标的外部检测是一个复杂的问题。在本研究中，我们采用主成分分析（PCA）技术对 68 波段特征层进行特征重建，以再生 PC1-68 波段，其中第一个主成分（PC 1）可以等效地定义为最大化投影数据方差的方向，并依次递减[50]。通过这种方法，每个特征都可以通过较少的波段来区分，实现数据降维。我们对每个瓦片内的粗一致区域分别进行了 PCA，这一步骤在 GEE 上执行。
3.2. 粗一致区域的超像素去除 多源数据预处理和重采样可以实现粗一致区域。事实上，像 Landsat 图像这样的数据已经可以提供更详细的结果。因此， 我们旨在净化粗一致区域并将其提升到精细分辨率。然后，我们从最初获得的精细一致区域中提取有效信息，并应用由 GEE 提供的 30 米空间分辨率 Landsat 图像来纠正多个来源之间的不确定性信息[31]。
然而，逐像素去除会产生大量斑点（“椒盐现象”），这是由于图像像素内的异质性[54]，并且大量像素的计数往往会超过 GEE 内存限制。超像素是由一系列具有相似特征（如颜色、亮度和纹理）的相邻像素组成的小区域。这些小区域中的大部分保留了有效信息，并且通常不会破坏图像中对象的边界信息。通过超像素分割算法将具有某些相似特征的像素分组，然后去除像素组（超像素）将大大减少椒盐现象，同时避免了由于冗余信息引起的计算复杂性[55]。本研究选择了 GEE 提供的简单非迭代聚类（SNIC）超像素算法，它是简单线性迭代聚类算法（SLIC）的改进版本，不需要迭代，效率更高，并且更适合边界保留。
此外，SNIC 需要一个大小参数来确定种子间隔，考虑到一致区域具有 300 米空间分辨率，而 Landsat 图像具有 30 米空间分辨率，我们通过实验确定大小值为 3，这足以对 Landsat 图像在一致区域内进行过度分割，避免分割不足对后期去除异常值的影响[56]。我们在 SNIC 中将紧凑度参数设置为 0，这禁用了空间距离加权，因为我们不希望生成的像素是规则的紧凑正方形；具有相似属性像素的聚合足以满足我们的需求。2015 年年中 Landsat 图像的六个原始波段均值被合成并用作输入分割的波段。SNIC 分割后，计算每个超像素内 PC1 和 PC2 的均值，并将其用作超像素的特征。通过聚合具有相同属性的图像元素，每个超像素的统计特征与真实地理对象更加一致，异常值更容易被检测到[57]。
然后，我们统计了生成的超像素的 PC1 和 PC2 属性。通常，异常值检测的阈值是通过经验或最优阈值搜索算法确定的，并且被移除的超像素覆盖了粗一致区域中的错误标记区域[28]。然而，粗一致区域是多类别的，如果同时进行会导致欠移除或过移除。因此，我们分别对每个类别统计并移除粗一致区域。经过我们的测试，每个特征的 PC1 和 PC2 都服从正态分布，因此我们将 PC1 和 PC2 结合起来，并应用 Lajda 准则构建异常值判别条件：
Class
(
i
)
=
{
consistency area
i
f
(
P
C
1
―
i
−
δ
P
C
1
i
≤
P
C
1
i
,
j
≤
P
C
1
―
i
+
δ
P
C
1
i
)
∩
(
P
C
2
―
i
−
δ
P
C
2
i
≤
P
C
2
i
,
j
≤
P
C
2
―
i
+
δ
P
C
2
i
)
inconsistency area else
Class
(
i
)
=
consistency area
i
f
P
C
1
¯
i
−
δ
P
C
1
i
≤
P
C
1
i
,
j
≤
P
C
1
¯
i
+
δ
P
C
1
i
∩
P
C
2
¯
i
−
δ
P
C
2
i
≤
P
C
2
i
,
j
≤
P
C
2
¯
i
+
δ
P
C
2
i
inconsistency area else
Class(i)={[" consistency area "if( bar(PC1)_(i)-delta_(PC1_(i)) <= PC1_(i,j) <= bar(PC1)_(i)+delta_(PC1_(i)))nn],[( bar(PC2)_(i)-delta_(PC2_(i)) <= PC2_(i,j) <= bar(PC2)_(i)+delta_(PC2_(i)))],[" inconsistency area else "]:} \operatorname{Class}(i)=\left\{\begin{array}{c}
\text { consistency area } i f\left(\overline{P C 1}_{i}-\delta_{P C 1_{i}} \leq P C 1_{i, j} \leq \overline{P C 1}_{i}+\delta_{P C 1_{i}}\right) \cap \\
\left(\overline{P C 2}_{i}-\delta_{P C 2_{i}} \leq P C 2_{i, j} \leq \overline{P C 2}_{i}+\delta_{P C 2_{i}}\right) \\
\text { inconsistency area else }
\end{array}\right. 其中
i
i
i i
j
j
j j
P
C
1
i
,
j
,
P
C
2
i
,
j
P
C
1
i
,
j
,
P
C
2
i
,
j
PC1_(i,j),PC2_(i,j) P C 1_{i, j}, P C 2_{i, j}
P
C
1
,
P
C
2
P
C
1
,
P
C
2
PC1,PC2 P C 1, P C 2
P
C
1
―
i
,
P
C
2
―
i
P
C
1
¯
i
,
P
C
2
¯
i
bar(PC1)_(i), bar(PC2)_(i) \overline{P C 1}_{i}, \overline{P C 2}_{i}
P
C
1
,
P
C
2
P
C
1
,
P
C
2
PC1,PC2 P C 1, P C 2
δ
P
C
1
i
,
δ
P
C
2
i
δ
P
C
1
i
,
δ
P
C
2
i
delta_(PC1_(i)),delta_(PC2_(i)) \delta_{P C 1_{i}}, \delta_{P C 2_{i}}
0.5
−
1.5
0.5
−
1.5
0.5-1.5 0.5-1.5 3.3. 本地适应样本集 根据 2.1 节，研究区域被划分为 31 个瓦片以减少纬度和经度差异对不同地物反射率的影响，并采用了一种局部自适应采样方法[35]。每个瓦片内每个地物的样本必须首先从该瓦片内每个类别的精细一致性区域内随机选择，并在净化后进行。我们仅对每类的样本数量设置了下限（1000），对总样本数量除了满足 GEE 操作的内存限制外没有上限。如果由于某些类别的精细一致性区域缺乏或一致性区域较少，相应瓦片中的样本不足，则从最近的相邻瓦片中补充缺失的样本。选定的样本仅用于训练当前瓦片。样本总数超过 300,000；这组被称为局部自适应样本集。使用这种局部自适应采样方法可以最小化由于空间分布导致的地物类别特征属性的不一致性。随机采样是通过在精确一致性区域内生成随机点来生成的。
3.4. 不一致区域的校正 在本研究中，通过第 3.3 节的本地适应性训练样本对随机森林（RF）分类器进行训练，旨在对不一致区域进行重新解释。RF 是一种包含多棵决策树的机器学习分类器，其输出类别由各棵树的输出类别的多数决定[59]。在遥感中，RF 相对于其他机器学习分类器具有对噪声不敏感、避免过拟合以及实现更高精度的优点。基于先前研究，我们选择了 300 棵树，使用
63
%
63
%
63% 63 \% 3.5. 验证和精度评估 为验证从精细一致性区域提取样本的可靠性，使用谷歌地球中的历史高空间分辨率图像对它们进行验证。我们从每个类别中随机选择了900个样本点（每个类别100个）[62]，统计了每个像素点的类型和九类混淆矩阵，并在移除前后计算了每个点被正确掩膜的准确性。
使用分层随机抽样生成的验证点对不一致区域的校正结果进行了验证和评估，这些验证点也通过谷歌地球进行了判读。验证点的数量为 1507 个（表 3）。还比较和评估了研究区内四个原始空间分辨率的多类别土地覆盖产品（MCD12Q1、CCI-LC、CGLS 和 FROM-GLC）以及四个单类别产品（GFSAD30、PALSAR、GSWD 和 GHS-BUILT）。值得注意的是，我们的精度评估仅在不一致区域进行，因为各种产品的精度差异主要体现在该区域[31]。相比之下，所有产品在精细一致区域的类别分布相同且完全可靠。我们为每个产品计算了混淆矩阵，并计算了常用的指标（生产者精度、用户精度和总体精度）来评估校准结果和四个多类别产品的空间分布精度[63,64]。
表3。研究区内校正结果精度评估的验证点。
类别 测试样本数量
耕地 355
森林 565
草地 155
灌木丛 76
水 86
城市建成区 100
裸地 62
永久性积雪和冰川 57
湿地 51
总计 1507
Class Number of Test Samples
Cropland 355
Forest 565
Grassland 155
Shrubland 76
Water 86
Urban/Built-up 100
Bare land 62
Permanent snow and ice 57
Wetland 51
Total 1507 | Class | Number of Test Samples |
| :--- | :--- |
| Cropland | 355 |
| Forest | 565 |
| Grassland | 155 |
| Shrubland | 76 |
| Water | 86 |
| Urban/Built-up | 100 |
| Bare land | 62 |
| Permanent snow and ice | 57 |
| Wetland | 51 |
| Total | 1507 |
4. 结果 4.1. 30米空间分辨率粗一致区域去除结果 图 3a 显示了超像素移除的结果。在计数像素数量后，大约有
38.75
%
38.75
%
38.75% 38.75 \% 我们验证了表 4 中显示的 900 个随机样本的混淆矩阵（图 5），发现灌木林和永久性雪冰的样本精度仅为
90
%
90
%
90% 90 \%
99
%
99
%
99% 99 \%
93.44
%
93.44
%
93.44% 93.44 \%
10
%
10
%
10% 10 \% 图 3. 移除粗一致区域：结果和统计。 (a) 移除的粗一致区域：结果。黑色代表移除的像素，其他是每个类别的精细一致区域。 (b) 移除结果的统计：红色是移除像素的百分比，蓝色是剩余像素移除的百分比。 图 4. 移除粗一致区域的视觉细节。基准图像是同一地点的高分辨率遥感图像。 (a-i) 分别是耕地、林地、草地、灌丛地、水体、城镇/建筑、裸地、永久性积雪和冰、湿地的一致区域移除结果。黑色对角线区域描绘了移除的图像元素。相应的彩色虚线区域描绘了对应类型净化后的精细一致区域。 图 5. 精细一致区域内随机像素点的验证结果。
表 4. 精细一致区域内随机像素点验证结果的混淆矩阵。
类别 耕地 森林 草地 灌木丛 水 城市/建设用地 裸地 永久性冰雪 湿地 总计
耕地 97
1
1
3
0
1
1
0
5
109
森林 0
98
1
2
0
0
0
0
0
101
草原 1
0
92
4
0
1
4
6
1
109
灌木丛 0
1
3
90
0
0
0
0
0
94
水域 1
0
0
0
99
1
0
0
3
104
城市/建成区 1
0
1
0
1
95
2
0
0
100
裸地 0
0
1
1
0
2
91
4
0
99
永久性积雪和冰川 0
0
1
0
0
0
2
90
0
93
湿地 0
0
0
0
0
0
0
0
91
91
总计 100
100
100
100
100
100
100
100
100
900
Class Cropland Forest Grassland Shrubland Water Urban/BuiltUp Bare Land Permanent Snow and Ice Wetland Total
Cropland 97 1 1 3 0 1 1 0 5 109
Forest 0 98 1 2 0 0 0 0 0 101
Grassland 1 0 92 4 0 1 4 6 1 109
Shrubland 0 1 3 90 0 0 0 0 0 94
Water 1 0 0 0 99 1 0 0 3 104
Urban/Built-up 1 0 1 0 1 95 2 0 0 100
Bare land 0 0 1 1 0 2 91 4 0 99
Permanent snow and ice 0 0 1 0 0 0 2 90 0 93
Wetland 0 0 0 0 0 0 0 0 91 91
Total 100 100 100 100 100 100 100 100 100 900 | Class | Cropland | Forest | Grassland | Shrubland | Water | Urban/BuiltUp | Bare Land | Permanent Snow and Ice | Wetland | Total |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| Cropland | 97 | 1 | 1 | 3 | 0 | 1 | 1 | 0 | 5 | 109 |
| Forest | 0 | 98 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 101 |
| Grassland | 1 | 0 | 92 | 4 | 0 | 1 | 4 | 6 | 1 | 109 |
| Shrubland | 0 | 1 | 3 | 90 | 0 | 0 | 0 | 0 | 0 | 94 |
| Water | 1 | 0 | 0 | 0 | 99 | 1 | 0 | 0 | 3 | 104 |
| Urban/Built-up | 1 | 0 | 1 | 0 | 1 | 95 | 2 | 0 | 0 | 100 |
| Bare land | 0 | 0 | 1 | 1 | 0 | 2 | 91 | 4 | 0 | 99 |
| Permanent snow and ice | 0 | 0 | 1 | 0 | 0 | 0 | 2 | 90 | 0 | 93 |
| Wetland | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 91 | 91 |
| Total | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 900 |
4.3. 不一致区域校正结果和精度评价 不一致区域的校正结果如图 6 所示。基于验证点的 PA、UA、OA 和混淆矩阵如表 5 所示。最终校正结果的总体精度为
85.8
%
85.8
%
85.8% 85.8 \%
91.5
%
91.5
%
91.5% 91.5 \%
73.65
%
73.65
%
73.65% 73.65 \%
94.25
%
94.25
%
94.25% 94.25 \%
68.97
%
68.97
%
68.97% 68.97 \%
5
%
5
%
5% 5 \%
12.25
%
12.25
%
12.25% 12.25 \%
9.47
%
9.47
%
9.47% 9.47 \% 图6。不一致区域的校正结果和验证点的分布。
表5。不一致区域验证点的混淆矩阵。
类别 耕地 森林 草原 灌木丛 水域 城市/建设用地 裸地 永久性积雪和冰川 湿地 总计 PA
OA
卡拉
耕地 307
24
11
0
4
3
0
0
6
355
86.48%
85.80%
0.82
森林 28
517
7
9
2
1
0
0
1
565
91.50%
草原 8
7
121
10
0
4
0
0
5
155
78.06%
灌木丛 7
7
6
56
0
0
0
0
0
76
73.68%
水 3
0
1
0
75
1
1
0
5
86
87.21%
城市/建成区 12
0
3
0
1
82
1
0
1
100
82.00%
裸地 2
1
7
1
2
0
49
0
0
62
79.03%
永久性积雪和冰川 0
0
2
0
0
0
9
46
0
57
80.70%
湿地 5
3
0
2
1
0
0
0
40
51
78.43%
总计 372
559
158
78
85
91
60
46
58
1507
UA
82.53%
92.49%
92.49%
76.58%
71.79%
88.24%
90.11%
81.67%
100.00%
68.97%
Class Cropland Forest Grassland Shrubland Water Urban/BuiltUp BARE LAND Permanent Snow and Ice Wetland Total PA OA Kарра
Cropland 307 24 11 0 4 3 0 0 6 355 86.48% 85.80% 0.82
Forest 28 517 7 9 2 1 0 0 1 565 91.50%
Grassland 8 7 121 10 0 4 0 0 5 155 78.06%
Shrubland 7 7 6 56 0 0 0 0 0 76 73.68%
Water 3 0 1 0 75 1 1 0 5 86 87.21%
Urban/Built-up 12 0 3 0 1 82 1 0 1 100 82.00%
Bare land 2 1 7 1 2 0 49 0 0 62 79.03%
Permanent snow and ice 0 0 2 0 0 0 9 46 0 57 80.70%
Wetland 5 3 0 2 1 0 0 0 40 51 78.43%
Total 372 559 158 78 85 91 60 46 58 1507
UA 82.53% 92.49% 92.49% 76.58% 71.79% 88.24% 90.11% 81.67% 100.00% 68.97% | Class | Cropland | Forest | Grassland | Shrubland | Water | Urban/BuiltUp | BARE LAND | Permanent Snow and Ice | Wetland | Total | PA | OA | Kарра |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| Cropland | 307 | 24 | 11 | 0 | 4 | 3 | 0 | 0 | 6 | 355 | 86.48% | 85.80% | 0.82 |
| Forest | 28 | 517 | 7 | 9 | 2 | 1 | 0 | 0 | 1 | 565 | 91.50% | | |
| Grassland | 8 | 7 | 121 | 10 | 0 | 4 | 0 | 0 | 5 | 155 | 78.06% | | |
| Shrubland | 7 | 7 | 6 | 56 | 0 | 0 | 0 | 0 | 0 | 76 | 73.68% | | |
| Water | 3 | 0 | 1 | 0 | 75 | 1 | 1 | 0 | 5 | 86 | 87.21% | | |
| Urban/Built-up | 12 | 0 | 3 | 0 | 1 | 82 | 1 | 0 | 1 | 100 | 82.00% | | |
| Bare land | 2 | 1 | 7 | 1 | 2 | 0 | 49 | 0 | 0 | 62 | 79.03% | | |
| Permanent snow and ice | 0 | 0 | 2 | 0 | 0 | 0 | 9 | 46 | 0 | 57 | 80.70% | | |
| Wetland | 5 | 3 | 0 | 2 | 1 | 0 | 0 | 0 | 40 | 51 | 78.43% | | |
| Total | 372 | 559 | 158 | 78 | 85 | 91 | 60 | 46 | 58 | 1507 | | | |
| UA | 82.53% | 92.49% | 92.49% | 76.58% | 71.79% | 88.24% | 90.11% | 81.67% | 100.00% | 68.97% | | | |
4.4. 与其他产品比较 4.4.1. 多类别产品比较 将校正后的不一致区域与精细一致区域镶嵌，以获得高空间分辨率土地覆盖融合结果。瓦片边界没有差异，如图 7a 所示。与其它多类别产品相比，所提出方法获得的地表覆盖结果在特征类分布方面有些相似（图
7
b
−
e
7
b
−
e
7b-e 7 \mathrm{~b}-\mathrm{e} 图 7. 校正产品和现有多类别产品的比较：(a) 融合结果；(b) CCI-LC；(c) CGLS；(d) FROM-GLC；(e) MCD12Q1。
然而，空间分辨率和映射方法的差异导致每个产品在细节上存在显著差异，如图 8 所示。由于 CCI-LC、CGLS 和 MCD12Q1 这三个产品的空间分辨率低于 30 米，许多细节被遗漏。例如，一些细小的河流（图 8(a1)）在这些三个产品（图 8(a3,a4,a6)）中被完全忽略，即使在空间分辨率较粗的 MCD12Q1 中，也只有森林类别。城市/建成区类型内部更为复杂（图 8(b1)），在空间分辨率较粗的产品（图 8(b3,b4,b6)）中，城市/建成区类型对内部信息的覆盖更加不完整。永久性积雪和冰、草地和裸地是 在空间分辨率较低的产品中未能准确表现（图 8(c3,c4,c6)）。FROM-GLC 与我们的融合结果具有相同的空间分辨率，尽管视觉上未忽略细节，但我们的土地覆盖产品在表现真实地形方面更准确。例如，图 8(a1,a2)中的许多裸地和草地区域在 FROM-GLC 中被错误地标记为耕地（图 8(a5)），图 8(b1)中的连续水体被错误地标记为建筑物（图 8(b5)），在图 8(c1)中，FROM-GLC 产品错误地将位于阴影中的部分裸地和草地标记为水体（图 8(c5)）。
耕地 森林 草地 灌木丛
水 城市/建成区 裸地 永久性冰雪 湿地
Cropland Forest Grassland Shrubland
Water Urban/ Built up Bare land Permanent snow and ice Wetland | Cropland | Forest | Grassland | Shrubland | |
| :--- | :--- | :--- | :--- | :--- |
| Water | Urban/ Built up | Bare land | Permanent snow and ice | Wetland |
图 8. 融合结果与其他四种多类别产品细节对比：(a-c) 为三个区域；(a1-c1) 为三个区域的遥感影像；(a2-c2) 为三个区域的融合结果；(a3-c3) 为 CCI-LC 在三个区域的类别分布；(a4-c4) 为 CGLS 在三个区域的类别分布；(a5-c5) 为 FROM-GLC 在三个区域的类别分布；(a6-c6) 为 MCD12Q1 在三个区域的类别分布。
四种多类别产品的不一致区域验证点精度评估结果分别显示在表 S2-S5 中，校准结果与四种产品（包括 PA、UA、OA 和 Kappa 系数）的比较评估结果如图 9 所示。图 9 显示，我们不一致区域校正结果中每种类型的 PA 都高于其他四种产品。在 UA 的比较中，除了草地、灌丛地和裸地较为模糊，以及湿地难以识别外，UA 并非最高。其他类型的 UA 均高于其他四种产品。我们的校正结果具有最高的 OA
(
85.80
%
)
(
85.80
%
)
(85.80%) (85.80 \%)
61.63
%
61.63
%
61.63% 61.63 \%
11.75
−
24.17
%
11.75
−
24.17
%
11.75-24.17% 11.75-24.17 \%
0.16
−
0.3
0.16
−
0.3
0.16-0.3 0.16-0.3 图9。与其他四种多类别产品验证结果的校正结果比较。
4.4.2. 单类产品比较 图 10 给出了在标定结果中对耕地、林地、水域和城市/建成区的比较，以及与四种单类产品（GFSAD30、PLASAR、GWSD 和 GHS-BUILT）的比较。校正结果中相应类别的分布仍然与四种产品一致。四种单类产品不一致区域的校正和验证结果混淆矩阵分别显示在表 S6-S13 中，比较评估结果显示在表 6 中。根据不一致区域相应类别准确率的比较，我们的校正结果的 PA 与现有的单类产品更接近，而 UA、OA 和 kappa 系数则高于相应的单类产品。不一致区域中每种类型的 OA 改善程度不同，森林的改善最大，而最小的是 仅用于水体。至于 kappa 系数，城市/建成区的最大改进为 0.56，水体的最小改进为 0.22。结果表明，与单一类别相比，我们的方法略有提高了准确性。 图 10。与四个单一类别产品校正结果的比较：（
a
1
,
a
2
)
a
1
,
a
2
)
a1,a2) \mathbf{a} \mathbf{1}, \mathbf{a} 2)
b
1
,
b
2
b
1
,
b
2
b1,b2 \mathbf{b} \mathbf{1}, \mathbf{b} 2
b
1
b
1
b1 \mathbf{b} 1
a
2
a
2
a2 \mathbf{a} 2
c
1
,
c
2
c
1
,
c
2
c1,c2 \mathbf{c} 1, \mathbf{c} 2
c
1
c
1
c1 \mathbf{c} 1
c
2
c
2
c2 \mathbf{c 2} 表6。与四个单一类别产品的精度评估结果比较。
土地覆盖产品 土地覆盖类型 生产者精度（%） 用户精度（%） 总体精度（%） 卡拉
GFSAD30 正确结果 耕地 86 75
83 52
93 78
0.8 0.47
PLASAR 正确结果 森林 92 66
92 64
94 73
0.87 0.43
正确结果 GWSD 水 87 79
88 59
99 96
0.87 0.65
正确结果 GHS-BUILT 城市/建成区 82 74
90 24
98 83
0.85 0.29
Land Cover Product Land Cover Type Producer's Accuracy (%) User's Accuracy (%) Overall Accuracy (%) Карра
Correct result GFSAD30 Cropland 86 75 83 52 93 78 0.8 0.47
Correct result PLASAR Forest 92 66 92 64 94 73 0.87 0.43
Correct result GWSD Water 87 79 88 59 99 96 0.87 0.65
Correct result GHS-BUILT Urban/Built-up 82 74 90 24 98 83 0.85 0.29 | Land Cover Product | Land Cover Type | Producer's Accuracy (%) | User's Accuracy (%) | Overall Accuracy (%) | Карра |
| :--- | :--- | :--- | :--- | :--- | :--- |
| Correct result GFSAD30 | Cropland | 86 75 | 83 52 | 93 78 | 0.8 0.47 |
| Correct result PLASAR | Forest | 92 66 | 92 64 | 94 73 | 0.87 0.43 |
| Correct result GWSD | Water | 87 79 | 88 59 | 99 96 | 0.87 0.65 |
| Correct result GHS-BUILT | Urban/Built-up | 82 74 | 90 24 | 98 83 | 0.85 0.29 |
5. 讨论 为了解决先前融合方法中由于重采样导致融合结果无法达到精细空间分辨率的问题[27,29,30]，本研究开发了一种基于 SNIC 分割算法和主成分分析技术的多源土地覆盖产品融合方法。结果表明，该方法可以获得具有精细空间分辨率的融合结果，并不同程度地提高了不一致区域的精度。该方法创造性地使用 SNIC 分割算法将粗一致区域内的 Landsat 图像层分割成像素组，并通过主成分分析和统计方法去除异常像素组。因此，粗一致区域可以净化为细一致区域，空间分辨率也从 300 米提高到 30 米，这在先前研究中从未做过。结果表明，在粗一致区域中移除了
38.75
%
38.75
%
38.75% 38.75 \% 同时，通过采样检测精细一致区域，样本准确率可以达到
93.44
%
93.44
%
93.44% 93.44 \% 在 GEE 的大规模数据集和大量样点（超过 30 万）的支持下，不一致区域得到了重新校正。获得了新的精细空间分辨率融合结果。结果表明，在四种可用的多类别产品中，CCILC、CGLS 和 MCD12Q1 的空间分辨率低于 30 米，并且提供的信息明显少于本研究中的融合结果。对于具有相同空间分辨率的 FROM-GLC，我们的融合结果在不一致区域提供了显著更高的精度。在不一致区域的整体精度可以达到
85.80
%
85.80
%
85.80% 85.80 \%
11.75
%
11.75
%
11.75% 11.75 \%
2.99
%
2.99
%
2.99% 2.99 \% 我们发现，在多类别土地覆盖制图中，耕地和草地容易混淆，灌木林是草地和裸地之间的过渡性土地覆盖类型[30]。湿地本质上是一种难以制图的土地覆盖类型，在湿地地区，湿地、耕地和草地也容易发生错分[52]。这些类型的光谱特征非常相似，因此在这些产品中，草地、灌木林和湿地的制图精度相对较低[14]。尽管我们的方法可以提高这些类型在不一致地区的制图精度，但仍然难以实现高精度，这在未来的研究中需要注意，以便更准确地区分这些类型。草地、裸地和永久性积雪和冰川分布在高纬度和高海拔地区，由于季节性变化导致的积雪和冰川融化，它们之间经常被误认为是彼此[65,66]，这也需要在未来的研究中加以注意。
尽管我们的方法能够整合多个土地覆盖产品中的有效信息并显著提高精度，但考虑到当前方法存在的问题，仍有一些方面可以改进，并应在未来继续研究。本文提出的方法需要在存在大量同一年份的土地覆盖数据且遥感影像完全覆盖的条件下进行。本研究选择的地表覆盖产品来自 2015 年，仅能使用 30 米空间分辨率的 Landsat 影像覆盖研究区域。随着土地覆盖数据的更新和高分辨率遥感影像（如 Sentinel 2 和 GF 系列）的开放获取，本文提出的方法仍然适用，并能获得空间更高分辨率的土地覆盖融合产品。然而，仍存在一些问题和改进的可能性。首先，大面积、高空间分辨率影像分割是一个挑战，确保精确分割和操作效率是一个待解决的问题。 其次，研究表明数据立方体可以解决全球遥感数据覆盖不足的问题，因此本文提出的方法也可以使用数据立方体获取高时间频率的土地覆盖产品[66]。最后，深度学习近年来在遥感领域显示出巨大的潜力[67]，与传统模型相比，通过深度学习融合和预测更准确的信息也将是我们的未来方向。
6. 结论 本研究提出了一种基于 SNIC 分割算法和 PCA 技术的多源土地覆盖产品融合方法。该方法解决了先前融合方法因重采样导致融合结果空间分辨率不精细的问题。它创造性地使用 SNIC 分割算法将粗一致区域内的 Landsat 图像层分割成像素组，并通过 PCA 和统计方法去除异常像素组，从而获得精细一致区域。该方法应用并使用 Google Earth 样本点进行了验证，结果表明从精细一致区域提取训练样本的准确率为
93.44
%
93.44
%
93.44% 93.44 \%
85.80
%
85.80
%
85.80% 85.80 \%
11.75
−
24.17
%
11.75
−
24.17
%
11.75-24.17% 11.75-24.17 \%
2.99
−
20.71
%
2.99
−
20.71
%
2.99-20.71% 2.99-20.71 \%
0.22
−
0.56
0.22
−
0.56
0.22-0.56 0.22-0.56 补充材料：以下支持信息可以下载于：https://www.mdpi.com/article/10.3390/rs14071676/s1，表 S1：本文使用的分类系统；表 S2：CCI-LC 的验证点混淆矩阵；表 S3：CGLS 的验证点混淆矩阵；表 S4：FROM-GLC 的验证点混淆矩阵；表 S5：MCD12Q1 的验证点混淆矩阵；表 S6：GFSAD30 评估结果的混淆矩阵；表 S7：PLASAR 评估结果的混淆矩阵；表 S8：GWSD 评估结果的混淆矩阵；表 S9：GHS-BUILT 评估结果的混淆矩阵；表 S10：耕地校正结果的混淆矩阵；表 S11：森林校正结果的混淆矩阵；表 S12：水域校正结果的混淆矩阵；表 S13：城市/建成区校正结果的混淆矩阵。
作者贡献：概念化、方法论，Q.J. 和 E.X.；验证、形式分析，Q.J.；资源，X.Z.；原始草稿撰写、审阅和编辑，Q.J.；监督、项目管理、资金获取，E.X.。所有作者均已阅读并同意已发表的文稿版本。
资助：本研究由第二次青藏高原科学考察与研究计划（STEP）（2019QZKK0603）、中国科学院战略性优先研究项目（XDA20040201）和中国科学院青年创新促进协会（2021052）资助。
机构审查委员会声明：不适用。 知情同意声明：不适用。 数据可用性声明：Google Earth Engine（https:/ / code.earthengine.google.com/ （访问于 2021 年 10 月 30 日）是一个免费且开放的平台。本研究中生成或使用的数据、模型或代码均可向通讯作者申请获取。 利益冲突：作者声明无利益冲突。
参考文献
刘纪远；匡威；张玉书；徐祥德；秦勇；宁家骏；周伟；张胜；李如忠；严春；等。中国自 20 世纪 80 年代末以来土地利用变化的空间时间特征、格局及成因。地理学报 2014，24，195-210。[DOI] 陆地，D.；田，H.；周，G.；葛，H. 使用多传感器遥感数据对中国东南部的人类居住地进行区域制图。遥感环境 2008，112，3668-3679。 [CrossRef] 塞托，K.C.；居纳尔普，B.；胡特拉，L.R. 对 2030 年城市扩张的全球预测以及对生物多样性和碳库的直接影响。美国国家科学院院报 2012，109，16083-16088。 [CrossRef]
