利用多源土地覆盖产品和时间序列 Sentinel-2 影像的自动训练样本采集
王燕钊,孙永华,曹旭岳,王一涵,张望宽,程兴路,王若增,宗金坤 引用本文:王燕钊,孙永华,曹旭岳,王一涵,张望宽,程兴路,王若增,宗金坤(2024)利用多源土地覆盖产品和时间序列 Sentinel-2 影像的自动训练样本采集,GIScience & Remote Sensing,61 卷第 1 期,2352957,DOI:10.1080/15481603.2024.2352957 链接到本文:https://doi.org/10.1080/15481603.2024.2352957
© 2024 作者。由 Informa UK Limited 出版,作为 Taylor & Francis 集团运营。
在线发表时间:2024年5月14日。
向本期刊提交您的文章
文章浏览量:456
查看 Crossmark 数据
利用多源土地覆盖产品和时间序列 Sentinel-2 影像的自动训练样本采集
程兴路,王若增,宗金坤
a,b,c,d
a,b,c,d
^("a,b,c,d ") { }^{\text {a,b,c,d }} 首都师范大学,北京,中国;首都师范大学城市环境过程与数字模拟国家重点实验室,北京,中国;
dKey
dKey
^("dKey ") { }^{\text {dKey }} 教育部三维信息获取与应用实验室,北京,中国
摘要
收集可靠的训练样本在提高土地覆盖(LC)制图产品的准确性方面起着关键作用,这些产品是全球环境和气候变化研究的重要基础数据。然而,该过程劳动强度大且耗时,因为它高度依赖人工解译。本文提出了一种利用多源土地覆盖产品和时间序列 Sentinel-2 影像的自动训练样本收集方法(ATSC)。首先,通过加权多数投票(WMV)算法融合多种土地覆盖产品,生成初步样本数据集。其次,应用局部选择性组合并行异常检测(LSCP)算法过滤异常样本。结果显示:(1)中国土地覆盖数据集(CLCD)整体准确率最高(73.22%),ESRI 土地覆盖(ESRI)整体准确率最低(59.93%)。树木覆盖、建筑区和水体在所有产品中准确率较高,而灌木丛和湿地的准确率普遍较低。(2)四个研究区初步训练样本的平均准确率为
95.62
%
95.62
%
95.62% 95.62 \% 。 然而,仍然存在异常样本,如分类错误、一年内的土地覆盖变化和光谱异常。(3)使用 LSCP 算法,去除了
70.10
%
70.10
%
70.10% 70.10 \% 的异常样本,最终各区域的训练样本准确率均超过了
97.95
%
97.95
%
97.95% 97.95 \% 。ATSC 方法为土地覆盖分类提供了更高质量的训练样本,促进了大规模土地覆盖制图工作。
文章历史
收稿日期:2023年11月29日 接受日期:2024年5月6日
关键词
地表覆盖;训练样本;自动样本采集;加权多数投票;异常检测
1. 引言
地表覆盖(LC)指地球表面各种生物或物理覆盖类型,主要指土地的自然属性(Di Gregorio 和 Jansen 2000)。准确且广泛的地表覆盖制图为自然与生物活动之间的关系、地表覆盖的空间格局、生态环境变化的模拟、监测与评估以及人类社会经济发展等科学研究提供基础数据支持(Cao 等 2015;Giri 等 2013;Kayet 等 2016)。随着遥感技术的快速发展及航空和卫星影像获取的便利,低成本且高效的地表覆盖数据获取已成为可能(Z. Chen 等 2019;Congalton 等 2014)。
全球范围内已发布了众多地表覆盖产品。许多早期的全球尺度产品,如国际地球生物圈计划(IGBP)DISCover(Loveland 和 Belward 1997),
马里兰大学(UMD)(Hansen 和 Reed 2000)和 2000 年全球土地覆盖数据集(GLC2000)(T. R. Loveland 等,2000)是利用 AVHRR 传感器数据创建的,导致空间分辨率较低(1 公里)。后来,500 米空间分辨率的 MODIS 和 300 米空间分辨率的 MERIS 成为土地覆盖制图的主要数据源。随后开发了几种中等空间分辨率的产品,如国家测绘组织的全球土地覆盖(GLCNMO)(Shirahata 等,2017)、MODIS 土地覆盖(MCD12Q1)(Friedl 等,2010)和土地覆盖-气候变化倡议(LCCCI)(Bontemps 等,2012)。近年来,大量的 Landsat 和 Sentinel-2 影像被用于生成空间分辨率为 30 米和 10 米的土地覆盖产品。这些高分辨率产品包括全球土地覆盖的更精细分辨率观测与监测(FROM-GLC)(Gong 等,2013)、全球土地覆盖
30 米分辨率制图(GlobeLand30)(J. Chen 等,2017)、Dynamic World(Brown 等,2022)和 ESA WorldCover(Zanaga 等,2021)。需要注意的是,土地覆盖产品的准确性存在差异。在最佳情况下,这些值约为
80
%
80
%
80% 80 \% ,且仍受不确定性的影响。以往研究者在特定区域评估和比较了各种土地覆盖产品,揭示了低准确率的持续挑战,尤其是在灌木地、湿地和草地类别中(L. Liu 等,2021;Z. Wang 和 Mountrakis,2023;J. Wang 等,2022;T. Zhao 等,2023)。有必要进一步努力提升土地覆盖图的整体质量。
许多学者主张通过融合现有的多源土地覆盖(LC)产品来减少 LC 地图的不确定性。Schepaschenko 等人(2011)使用适宜性指数融合多种数据源,生成了俄罗斯地区的 LC 地图。Kinoshita 等人(2014)采用逻辑回归模型融合了六个全球尺度的 LC 产品,揭示了产品数量对融合结果准确性有显著影响。A. Pérez-Hoyos 等人(2012)运用模糊集理论融合了欧洲地区的四个 LC 产品,提升了准确性。因此,随着可用 LC 产品数量的不断增加,多源产品的融合已成为提高地图准确性的关键方法。
训练样本对于训练分类器至关重要,直接影响土地覆盖(LC)制图结果的准确性和可靠性(Foody 和 Mathur 2006)。收集训练样本有两种方法。第一种也是最常用的方法是基于解译的,这种方法能获得高质量的样本,但在大规模制图时需要大量的人工工作(Calderón-Loor、Hadjikakou 和 Bryan 2021;M. Li 等 2022)。第二种方法是从现有的土地覆盖产品中收集训练样本,已被证明具有全自动收集和生成大规模、地理分布广泛的训练数据集等优势(Colditz 等 2011;Hermosilla 等 2022;Hu、Dong 和 Batunacun 2018;Radoux 等 2014;H. K. Zhang 和 Roy 2017;X. Zhang 等 2021)。然而,在许多研究中,采用第二种方法时通常选择空间分辨率较低的土地覆盖产品,且训练样本大多来自单一产品。分类错误和土地覆盖变化可能影响训练样本的可靠性, 特别是在景观高度破碎的地区,可能导致训练样本质量较低。近年来出现了许多高空间分辨率的土地覆盖(LC)产品,但很少有学者将它们融合用于训练样本的收集。单一来源的 LC 产品所获得信息的可靠性通常低于多产品融合所获得的信息(Ran 等,2012)。多种 LC 产品的融合可以弥补单一产品的不足,减少误差,增强所收集训练样本的可信度。具有高空间分辨率的 LC 产品能够提供丰富的空间细节。从高空间分辨率产品中收集的训练样本将更加准确,且适用于高分辨率的 LC 制图。此外,许多以往的方法未考虑在收集样本点位置存在卫星影像缺失的问题。对缺失值敏感的算法可能会影响模型性能(C. Zhang, Zhang, and Tian 2023)。 因此,有必要利用多种高空间分辨率的土地覆盖产品和卫星影像来收集高质量且云量较少的训练样本。
研究人员在自动收集训练样本过程中,由于土地覆盖产品固有的不确定性,采用了多种方法去除异常值。Radoux 等人(2014)通过形态学处理去除边缘区域的像素,并基于马氏距离测量排除光谱异常的像素。Zhang 等人(2021)基于训练样本的光谱统计分布去除异常点。Jin 等人(2022)使用主成分分析(PCA)降低 68 个光谱特征的维度,并通过统计方法去除异常像素。Wen 等人(2022)利用 NDVI 时间序列数据计算月均值(
μ
μ
mu \mu )和标准差(
σ
σ
sigma \sigma ),去除 NDVI 中落在
μ
±
σ
μ
±
σ
mu+-sigma \mu \pm \sigma 范围之外的玉米样本。然而,许多异常值去除方法未考虑时间序列光谱特征,导致无法检测一年内土地覆盖变化引起的异常。大多数方法采用基于统计的异常值去除技术,假设数据来自特定分布(Chandola, Banerjee, 和 Kumar 2009),通常是正态分布。 然而,这一假设在高维真实数据集中往往不成立。 即使统计假设可以合理地被证明,仍然可以应用多种假设检验统计量来检测异常;选择最佳统计量往往不是一项简单的任务。近年来,异常检测受到了机器学习研究人员的更多关注。许多新颖的算法被提出并应用于各个领域,包括金融欺诈检测、网络病毒攻击预警和自然灾害防范(J. Li 等,2023;Nassif 等,2021;Proverbio、Bertola 和 Smith,2018;Zuo 等,2023)。遗憾的是,许多先进的异常检测算法尚未被用于消除土地覆盖训练样本中的异常值。与统计技术相比,基于机器学习的异常检测算法能够适应各种类型的数据分布,包括高维复杂情况。机器学习方法具有更广泛的适用性和更强的鲁棒性,同时还提供异常分数,有助于理解异常的严重程度(Han 等,2022;Nassif 等,2021)。 特别是,基于集成学习的异常检测算法通过减少对单个检测器的依赖,可以提高检测的准确性和鲁棒性(Ouyang et al. 2021;J. Zhang et al. 2019)。
为填补这些研究空白,开发了一种自动收集土地覆盖(LC)训练样本的方法。该方法能够减少自动训练样本收集中的不确定性,生成更高质量且更具代表性的土地覆盖训练样本。该方法适用于细空间分辨率的土地覆盖分类。主要内容包括:(1)使用统一的验证数据集评估多个土地覆盖产品的准确性。(2)根据用户准确率计算类别权重值,并基于加权多数投票(WMV)算法融合多个土地覆盖产品。(3)从融合的土地覆盖图中提取稳定区域,即高置信度且无云区域,随后采用局部自适应策略收集训练样本。(4)实施基于集成的异常检测算法,利用从 Sentinel-2 影像提取的时间序列光谱特征识别并剔除异常样本。通过该方法收集的可靠训练样本可用于特定区域的大尺度和细尺度土地覆盖制图。
2. 研究区域与材料
2.1. 研究区域
中国选择了四个地理和气候各异的地区作为研究区域。四个地区分别是京津冀地区(218,000
km
2
km
2
km^(2) \mathrm{km}^{2} )、黑龙江省
(
473
,
000
km
2
)
473
,
000
km
2
(473,000km^(2)) \left(473,000 \mathrm{~km}^{2}\right) 、广东省(
180
,
000
km
2
180
,
000
km
2
180,000km^(2) 180,000 \mathrm{~km}^{2} ),以及新疆北部(184,000平方公里,包括博尔塔拉蒙古自治州、伊犁哈萨克自治州、塔城地区、霍城县、可克达拉市、克拉玛依市、石河子市和双河市)。
图1显示了这些区域位于中国的不同地区。新疆北部地区位于中国西北部,主要属于温带大陆性气候。该气候特点为冬季寒冷,夏季炎热,昼夜和年温差显著,年降水量相对较低。北京-天津-河北地区位于北方,黑龙江省位于东北,两者均属温带季风气候。黑龙江地处较高纬度,冬季更长且更寒冷。广东省位于南部沿海,属于亚热带季风气候,全年气温温暖。比较不同区域的结果可以全面理解所提方法的泛化能力。
2.2. 数据与预处理
2.2.1. 地表覆盖产品
利用 2020 年的多种土地覆盖产品收集了训练样本。使用了四种土地覆盖产品,分别是 Dynamic World(DW)、欧洲航天局 WorldCover(ESA)、环境系统研究所土地覆盖(ESRI)和中国土地覆盖数据集(CLCD)(Brown 等,2022;Zanaga 等,2021;Karra 等,2021;J. Yang 和 Huang,2021)。所有这些产品均可通过 Google Earth Engine(GEE)访问。空间分辨率为 10 米,CLCD 除外,其分辨率为 30 米。选择这四种土地覆盖产品是因为它们具有较高的分辨率,并提供多个参考年份的地图,支持跨年份的样本收集。然而,这些土地覆盖产品均存在偏差和不确定性。差异在于
图 1. 四个研究区的地理位置。(a) 新疆北部地区;(b) 黑龙江省;(c) 广东省;(d) 京津冀地区。 分类系统(不同产品对某些类别的定义略有不同)是影响不同产品一致性和偏差的重要因素(Hao 等,2023;Kang 等,2022;Venter 等,2022)。例如,DW 和 ESRI 中的稻田和灌溉/淹没农业被归类为淹水植被,而在 ESA 中则被归类为农田。此外,分类方法、数据来源和预处理技术也会影响各产品之间的一致性(Hao 等,2023;J. Wang 等,2022)。具体来说,CLCD 是通过使用来自 GEE 的 335,709 幅 Landsat 图像构建多个时间指标,并将其输入随机森林分类器获得的(J. Yang 和 Huang,2021)。而 DW 和 ESRI 产品均源自 Sentinel-2 数据并采用深度学习方法,但其数据预处理方法和输入特征不同(Venter 等,2022)。还需注意的是,上述产品的精度验证结果是由数据生产者使用不同的验证样本获得的(T. Zhao 等,2023)。 验证数据集的数量和质量也可能导致评估结果中的错误。
LC 产品的预处理是在 GEE 上进行的。鉴于 DW 提供近实时的 LC 数据, 年度 DW 数据通过对 2020 年 1 月 1 日至 12 月 31 日期间所有可用数据的每个像素取众数进行合成。LC 产品最初使用中国边界矢量数据进行裁剪。所有数据均通过重投影操作标准化为 WGS1984 坐标系统,空间分辨率为 10 米。通过重分类过程统一了 LC 产品的分类系统。生成的类别代码(像素值)及名称如下:1-树木覆盖,2-灌木丛,3-草地,4-耕地,5-建筑/不透水区域,6-裸地,7-积雪和冰,8-水体,9-湿地/洪泛植被。值得注意的是,ESA 中的“草本湿地”和“红树林”类别被合并为“湿地”类别。ESA 中存在的“苔藓和地衣”类别在其他三个 LC 产品中缺失,且在中国分布有限,因此该类别被排除。
2.2.2. 土地覆盖验证数据集
由刘等人(2023 年)开发的 2020 年新型分层随机抽样全球验证数据集(SRS_Val 数据集)(https://zenodo.org/records/7846090 ,访问时间 2023 年 6 月 12 日)被用于评估土地覆盖产品的准确性。SRS_Val 数据集是 采用分层等面积随机抽样策略和目视解译方法(赵涛等,2023 年)建立的。与以往的验证数据集相比,该数据集提高了异质景观和稀有土地覆盖类型的样本密度。它采用联合国土地覆盖分类系统(UN-LCCS)的标准化分类体系,确保与各种土地覆盖产品的良好兼容性和一致性。尽管整合了多源遥感影像并实施了严格的质量控制措施以确保 SRS_Val 数据集的高置信度,但部分验证样本仍存在不确定性。需要注意的是,在国家尺度评估 SRS_Val 数据集时,某些类别的样本可能相对较少。
对 SRS_Val 也进行了裁剪和重分类处理。“阔叶林”、“针叶林”等类似类别合并为“树木覆盖”。“雨养农田”和“灌溉农田”合并为“农田”。 “稀疏植被”和“裸地”合并为“裸地”。由于 SRS_Val 旨在进行全球尺度的精度评估,某些类别在中国区域内的样本可能有限。为减轻因样本有限导致的精度评估不确定性,通过谷歌地球影像的目视解译将样本较少的类别(如湿地和雪/冰)数量增加到 300。最终,验证样本共有 7,691 个,如图 2 所示。
2.2.3. 卫星影像
哨兵-2 卫星配备了一种高分辨率多光谱成像仪,称为多光谱成像仪(MSI)。它被广泛用于土地监测,提供植被、土壤和水体覆盖、内陆水道及沿海区域的影像(Drusch 等,2012;Phiri 等,2020)。该卫星
图2. 中国7691个土地覆盖验证样本的空间分布。 轨道高度786公里,捕获13个光谱波段的影像,扫描带宽为290公里。该卫星单颗卫星的重访周期为10天,双卫星互补覆盖的重访周期为5天,导致不同光谱范围内的空间分辨率有所不同。
Sentinel-2 Level-2A 数据作为卫星数据源。选择了 2020 年 3 月至 11 月的 Sentinel-2 数据,并在 GEE 平台上进行了裁剪。平台的云去除功能被用来消除云和云影像素,均值函数用于合成相应的月度影像。
3. 方法
已开发出一种自动训练样本收集方法(ATSC),利用多源数据进行操作
土地覆盖产品和时间序列卫星影像。该方法能够生成高质量且可靠的样本数据集。图 3 展示了 ATSC 方法的技术流程图。多源土地覆盖产品融合、初步样本收集和光谱特征提取均在 GEE 上实现。异常检测算法使用 Python 实现。各组件的详细内容将在以下章节中介绍。
3.1. 多源土地覆盖产品融合
所使用的土地覆盖产品提供了土地覆盖类别的地图,只有 DW 产品额外提供了每个类别的置信度层。因此,WMV 算法是一种简单且适用的数据融合方法。对于简单多数投票算法,所有土地覆盖产品的投票权值相同,均为 1。给定像素的最终分类结果为
图 3. ATSC 方法的技术流程图。
图4. 用于加权多数投票算法的四个土地覆盖产品各类别的权重值。 通过选择获得最多票数的类别来确定。需要注意的是,每个土地覆盖产品的准确性不同,各类别的分类准确率存在差异。通常需要考虑产品之间的差异,以发挥各产品的优势,实现更有效且更高准确度的组合结果。这就需要在投票过程中为每个土地覆盖产品的投票赋予不同的权重,这种算法称为加权多数投票(WMV)(H. Kim 等,2011;Zhu 等,2021)。每个像素的最终预测基于最高加权票数。WMV 算法的公式如下:
R
k
(
x
)
=
∑
i
=
1
C
i
,
k
=
max
j
=
1
{
R
j
(
x
)
=
∑
i
=
1
w
i
,
j
C
i
,
j
}
R
k
(
x
)
=
∑
i
=
1
C
i
,
k
=
max
j
=
1
R
j
(
x
)
=
∑
i
=
1
w
i
,
j
C
i
,
j
R_(k)(x)=sum_(i=1)C_(i,k)=max_(j=1){R_(j)(x)=sum_(i=1)w_(i,j)C_(i,j)} R_{k}(x)=\sum_{i=1} C_{i, k}=\max _{j=1}\left\{R_{j}(x)=\sum_{i=1} w_{i, j} C_{i, j}\right\}
在此上下文中,
w
i
,
j
w
i
,
j
w_(i,j) w_{i, j} 表示土地覆盖产品
C
i
C
i
C_(i) C_{i} 在类别
j
j
j j 的权重值。
所有产品中每个类别的权重值最初设定为 0.25。随后,使用统一的验证数据集(SRS_Val)评估了中国地区四个产品的准确性。初始权重基于用户准确率(UA)进行了调整。选择 UA 的原因是验证数据集中样本数量相对有限。生产者准确率(PA)可能会受到类别间样本不平衡的影响。权重值计算的具体过程详见公式(2)。为了衡量 融合的土地覆盖分类结果的置信度,每个像素的置信值定义为 WMV 算法计算的最高加权投票值。四个土地覆盖产品分类结果的一致性越高,置信值越高。
w
i
,
j
=
U
A
(
C
i
,
j
)
∑
i
=
1
U
A
(
C
i
,
j
)
w
i
,
j
=
U
A
C
i
,
j
∑
i
=
1
U
A
C
i
,
j
w_(i,j)=(UA(C_(i,j)))/(sum_(i=1)^(UA(C_(i,j)))) w_{i, j}=\frac{U A\left(C_{i, j}\right)}{\sum_{i=1}^{U A\left(C_{i, j}\right)}}
其中
U
A
(
C
i
,
j
)
U
A
C
i
,
j
UA(C_(i,j)) U A\left(C_{i, j}\right) 表示土地覆盖产品
C
i
C
i
C_(i) C_{i} 中类别
j
j
j j 的用户准确率,如公式(3)所示。其他使用的准确性评估指标包括生产者准确率(PA)、总体准确率(OA)和卡帕系数,详见公式(4)-(6)。
U
A
(
class
i
)
=
x
i
i
∑
k
=
1
n
x
i
k
U
A
class
i
=
x
i
i
∑
k
=
1
n
x
i
k
UA(" class "_(i))=(x_(ii))/(sum_(k=1)^(n)x_(ik)) U A\left(\text { class }_{i}\right)=\frac{x_{i i}}{\sum_{k=1}^{n} x_{i k}}
P
A
(
class
i
)
=
x
i
i
∑
k
=
1
n
x
k
i
P
A
class
i
=
x
i
i
∑
k
=
1
n
x
k
i
PA(" class "_(i))=(x_(ii))/(sum_(k=1)^(n)x_(ki)) P A\left(\text { class }_{i}\right)=\frac{x_{i i}}{\sum_{k=1}^{n} x_{k i}}
O
A
=
∑
k
=
1
n
x
k
k
∑
i
,
k
=
1
n
x
i
k
O
A
=
∑
k
=
1
n
x
k
k
∑
i
,
k
=
1
n
x
i
k
OA=(sum_(k=1)^(n)x_(kk))/(sum_(i,k=1)^(n)x_(ik)) O A=\frac{\sum_{k=1}^{n} x_{k k}}{\sum_{i, k=1}^{n} x_{i k}}
Kappa
=
N
∑
i
=
1
n
x
i
i
−
∑
k
=
1
n
(
∑
i
=
1
n
x
i
k
∙
∑
i
=
1
n
x
k
i
)
N
2
−
∑
k
=
1
n
(
∑
i
=
1
n
x
i
k
∙
∑
i
=
1
n
x
k
i
)
Kappa
=
N
∑
i
=
1
n
x
i
i
−
∑
k
=
1
n
∑
i
=
1
n
x
i
k
∙
∑
i
=
1
n
x
k
i
N
2
−
∑
k
=
1
n
∑
i
=
1
n
x
i
k
∙
∑
i
=
1
n
x
k
i
" Kappa "=(Nsum_(i=1)^(n)x_(ii)-sum_(k=1)^(n)(sum_(i=1)^(n)x_(ik)∙sum_(i=1)^(n)x_(ki)))/(N^(2)-sum_(k=1)^(n)(sum_(i=1)^(n)x_(ik)∙sum_(i=1)^(n)x_(ki))) \text { Kappa }=\frac{N \sum_{i=1}^{n} x_{i i}-\sum_{k=1}^{n}\left(\sum_{i=1}^{n} x_{i k} \bullet \sum_{i=1}^{n} x_{k i}\right)}{N^{2}-\sum_{k=1}^{n}\left(\sum_{i=1}^{n} x_{i k} \bullet \sum_{i=1}^{n} x_{k i}\right)}
其中
N
N
N N 是所有验证样本的数量;
n
n
n n 是矩阵的行/列数;
x
i
i
x
i
i
x_(ii) x_{i i} 是 混淆矩阵中第
i
i
i i 行第
i
i
i i 列的元素;以及混淆矩阵中第
i
i
i i 行第
k
k
k k 列的元素为
x
i
k
x
i
k
x_(ik) x_{i k} (
x
k
k
x
k
k
x_(kk) x_{k k} 和
x
k
i
x
k
i
x_(ki) x_{k i} 类似)。
3.2. 通过局部自适应策略收集初步训练样本
为了确保收集到高度可靠的训练样本,提取在土地覆盖分类结果中具有高置信度且受云覆盖影响最小的稳定区域是必不可少的。采用了两个约束条件来提取高置信度区域:(1)融合土地覆盖产品后像素的置信值(即最高加权投票)大于 0.7;(2)对于每个像素,至少有四个土地覆盖产品中的三个分类结果一致。这两个约束条件有助于减轻某一特定土地覆盖产品在某一类别中准确率较差所带来的负面影响。该方法还增强了从分类准确率和一致性较差的类别中收集训练样本的可信度。随后应用形态学处理(即腐蚀)去除边缘区域的像素,从而减少边缘效应引起的不确定性(Wen et al. 2022)。为了确保收集的训练样本受云覆盖掩膜影响最小,从时间序列卫星影像中提取了无数据缺口的区域。训练样本从最终提取的稳定区域中收集。
提取稳定区域后,采用局部自适应策略收集初步训练样本。每个研究区最初生成一个网格数据集(每个网格为
5
km
×
5
km
5
km
×
5
km
5kmxx5km 5 \mathrm{~km} \times 5 \mathrm{~km} )。在每个网格内为每个 LC 类别收集一个样本点。选择该策略是考虑到同一区域内相同的 LC 类别可能因地理位置影响表现出光谱差异。因此,通过逐网格收集训练样本,所收集样本的空间分布更加均匀且具有代表性。多次实验表明,将网格大小设置为 5 公里能够收集足够的训练样本,实现大区域和城市级别的 LC 制图。此外,在 GEE 上执行时,由于收集的样本数量较多,该过程不会显著消耗过多的计算时间和资源。
此外,某些类别,如湿地,在该区域内的空间分布可能相对稀疏。从这些类别获得的初步样本数量可能有限。为了解决这一问题,设定了初步样本数量的初始阈值。考虑到不同区域的面积可能存在显著差异,有必要为每个区域设置不同的阈值(见公式(7))。如果某一类别收集的初步样本总数低于设定的阈值,则在每个网格内为该类别收集更多样本。为了避免某一类别的初步样本过于密集或过于接近,从而影响模型的泛化能力和分类准确性,每个网格收集的样本最大数量设为5个。如果某一类别在每个网格中收集的样本超过五个,即使样本总数未达到设定阈值,也将停止该类别的采样。 值得注意的是,该策略中的网格大小和样本数量阈值仍可由用户根据研究领域和具体应用需求进行调整。
X
=
500
+
500
⋅
S
−
10
10
X
=
500
+
500
⋅
S
−
10
10
X=500+500*(S-10)/(10) X=500+500 \cdot \frac{S-10}{10}
X
X
X X 表示初始阈值,
S
S
S S 表示区域面积(以万平方公里为单位),
x
x
x x 表示上取整函数。
3.3. 利用时间序列光谱特征过滤异常样本
使用了 3 月至 11 月的月度合成 Sentinel-2 影像来计算五个额外的光谱指数,包括归一化植被指数(NDVI)、增强植被指数(EVI)、归一化水体指数(NDWI)、归一化建筑指数(NDBI)和裸土指数(BSI)。在收集初步训练样本后,提取了表 1 中显示的每个月的 11 个光谱特征,总计 99 个特征。时间序列光谱特征被用于从初步训练样本中过滤异常样本。
表 1. 从 Sentinel-2 影像中提取用于过滤异常样本的光谱特征。
缩写
公式
参考
B2
蓝色
B3
绿色
B4
红色
B8
NIR
B11
SWIR1
B12
SWIR2
NDVI
(近红外 - 红光)/(近红外 + 红光)
Tucker (1979)
EVI
2.5
×
(
2.5
×
(
2.5 xx( 2.5 \times( (近红外 - 红光)/(近红外
+
6
×
+
6
×
+6xx +6 \times 红光
−
7.5
×
−
7.5
×
-7.5 xx -7.5 \times 蓝光 +1
)
)
) ) )
Huete 等, (1997)
NDWI
(绿光 - 近红外)/(绿光 + 近红外)
高(1996)
NDBI
(短波红外 - 近红外)/(短波红外 + 近红外)
Zha 等,(2003)
BSI
(
(
( ( SWIR2 + 红光)
−
(
−
(
-( -( 近红外 + 蓝光)/(SWIR2 + 红光)
+
(
+
(
+( +( (近红外 + 蓝光
)
)
) )
Diek 等人,(2017)
Abbreviation Formulation Reference
B2 Blue
B3 Green
B4 Red
B8 NIR
B11 SWIR1
B12 SWIR2
NDVI (NIR - Red)/(NIR + Red) Tucker (1979)
EVI 2.5 xx( (NIR - Red)/(NIR +6xx Red -7.5 xx Blue +1) ) Huete et al., (1997)
NDWI (Green - NIR)/(Green + NIR) Gao (1996)
NDBI (SWIR - NIR)/(SWIR + NIR) Zha et al., (2003)
BSI ( SWIR2 + Red) -( NIR + Blue)/(SWIR2 + Red) +( (NIR + Blue ) Diek et al., (2017) | Abbreviation | Formulation | Reference |
| :--- | :--- | :--- |
| B2 | Blue | |
| B3 | Green | |
| B4 | Red | |
| B8 | NIR | |
| B11 | SWIR1 | |
| B12 | SWIR2 | |
| NDVI | (NIR - Red)/(NIR + Red) | Tucker (1979) |
| EVI | $2.5 \times($ (NIR - Red)/(NIR $+6 \times$ Red $-7.5 \times$ Blue +1$)$ ) | Huete et al., (1997) |
| NDWI | (Green - NIR)/(Green + NIR) | Gao (1996) |
| NDBI | (SWIR - NIR)/(SWIR + NIR) | Zha et al., (2003) |
| BSI | $($ SWIR2 + Red) $-($ NIR + Blue)/(SWIR2 + Red) $+($ (NIR + Blue $)$ | Diek et al., (2017) |
卫星图像的时间序列选择了从三月到十一月,经过了仔细考虑。首先,这一时间段涵盖了大多数植被类型(如树木、灌木和草地)的生长季节,能够准确反映不同类别的特征。其次,对于某些土地覆盖类别,如建筑区和裸地,季节变化相对较小,这一时期的数据足以描述其特征。此外,中国中高纬度和高海拔的某些地区冬季可能受积雪覆盖和植被落叶影响,可能无法很好地代表地表特征。因此,选择三月至十一月的数据可以减少此类干扰,更准确地反映地表特征。
3.3.2. 异常检测算法
3.3.2.1. 局部离群因子(LOF)算法。采用了一种集成策略进行异常检测,以识别初步样本中的异常样本。LOF 作为基础异常检测器,通过计算每个数据点与其邻近点之间的密度差异来确定异常程度(Breunig 等,2000)。LOF 算法不仅考虑单个样本点的密度,还考虑邻近点的密度,使其适应于密度和分布各异的数据。因此,它适用于地理环境中同一土地覆盖类别的光谱特征因气候、位置等因素差异而变化的情况。此外,LOF 算法对高维数据集也具有鲁棒性,无需对数据分布做出明确假设。包括全局和局部异常在内的不同类型的异常点都可以被 检测到(D. Kim, Lee, 和 Lee 2020)。LOF 算法的详细描述如下,算法中计算的距离均指派生变量特征空间中的距离。
对于给定的样本点
x
i
x
i
x_(i) x_{i} ,令
D
k
(
x
i
)
D
k
x
i
D_(k)(x_(i)) D_{k}\left(x_{i}\right) 表示
x
i
x
i
x_(i) x_{i} 与其 k 近邻之间的距离,令
L
k
(
x
i
)
L
k
x
i
L_(k)(x_(i)) L_{k}\left(x_{i}\right) 表示在
k
k
k k 近邻距离内的点集。则两个样本点
x
i
x
i
x_(i) x_{i} 和
x
j
x
j
x_(j) x_{j} 之间的可达距离,记为
R
k
(
x
i
,
x
j
)
R
k
x
i
,
x
j
R_(k)(x_(i),x_(j)) R_{k}\left(x_{i}, x_{j}\right) ,计算如下:
R
k
(
x
i
,
x
j
)
=
max
{
dist
(
x
i
,
x
j
)
,
D
k
(
x
j
)
}
R
k
x
i
,
x
j
=
max
dist
x
i
,
x
j
,
D
k
x
j
R_(k)(x_(i),x_(j))=max{dist(x_(i),x_(j)),D_(k)(x_(j))} R_{k}\left(x_{i}, x_{j}\right)=\max \left\{\operatorname{dist}\left(x_{i}, x_{j}\right), D_{k}\left(x_{j}\right)\right\}
当
j
j
j j 位于密集区域,且
x
i
x
i
x_(i) x_{i} 远离
x
j
x
j
x_(j) x_{j} 时,可达距离度量等于实际距离。如果
j
j
j j 位于稀疏区域,可达距离度量将被其 k 近邻距离平滑处理。这使我们能够通过对
x
i
x
i
x_(i) x_{i} 的
k
k
k k 近邻点的可达距离取平均,计算
x
i
x
i
x_(i) x_{i} 的平均可达距离
(
A
R
k
(
x
i
)
)
A
R
k
x
i
(AR^(k)(x_(i))) \left(A R^{k}\left(x_{i}\right)\right) :
A
R
k
(
x
i
)
=
M
E
A
N
j
∈
L
k
(
x
i
)
R
k
(
x
i
,
x
j
)
A
R
k
x
i
=
M
E
A
N
j
∈
L
k
x
i
R
k
x
i
,
x
j
AR^(k)(x_(i))=MEAN_(j inL_(k)(x_(i)))R_(k)(x_(i),x_(j)) A R^{k}\left(x_{i}\right)=M E A N_{j \in L_{k}\left(x_{i}\right)} R_{k}\left(x_{i}, x_{j}\right)
局部离群因子是
A
R
k
(
x
i
)
A
R
k
x
i
AR^(k)(x_(i)) A R^{k}\left(x_{i}\right) 相对于其
k
k
k k 近邻的平均比值
x
i
x
i
x_(i) x_{i} 。
L
O
F
k
(
x
i
)
=
M
E
A
N
y
i
∈
L
k
(
x
i
)
A
R
k
(
x
i
)
A
R
k
(
x
j
)
L
O
F
k
x
i
=
M
E
A
N
y
i
∈
L
k
x
i
A
R
k
x
i
A
R
k
x
j
LOF_(k)(x_(i))=MEAN_(y_(i)inL_(k)(x_(i)))(AR^(k)(x_(i)))/(AR^(k)(x_(j))) L O F_{k}\left(x_{i}\right)=M E A N_{y_{i} \in L_{k}\left(x_{i}\right)} \frac{A R^{k}\left(x_{i}\right)}{A R^{k}\left(x_{j}\right)}
用于异常检测的光谱特征有很多,这些特征可能表现出高度相关性(尤其是在相邻月份使用相同波长波段时)。使用传统的欧氏距离度量可能无法获得满意的结果。因此,LOF 算法采用了马氏距离作为距离度量。
计算马氏距离的公式如下(De Maesschalck、Jouan-Rimbaud 和 Massart 2000):
d
M
(
x
i
,
x
j
)
=
(
x
i
−
x
j
)
T
Σ
−
1
(
x
i
−
x
j
)
d
M
x
i
,
x
j
=
x
i
−
x
j
T
Σ
−
1
x
i
−
x
j
d_(M)(x_(i),x_(j))=sqrt((x_(i)-x_(j))^(T)Sigma^(-1)(x_(i)-x_(j))) d_{M}\left(x_{i}, x_{j}\right)=\sqrt{\left(x_{i}-x_{j}\right)^{T} \Sigma^{-1}\left(x_{i}-x_{j}\right)}
其中
Σ
Σ
Sigma \Sigma 是向量
x
i
x
i
x_(i) x_{i} 和
x
j
x
j
x_(j) x_{j} 的协方差矩阵。 3.3.2.2. 基于集成的异常检测算法:LSCP。作为一种异常检测的集成方法,采用了赵等人(2019)提出的并行异常检测器中的局部选择组合(LSCP)。LSCP 首先为测试实例定义一个局部区域,然后通过测量与伪真值的相似性,识别该局部区域内最有能力的基础检测器。通过这一集成过程可以实现更稳健的预测(Y. Zhao 等,2019)。具体步骤如下: (i) 最初,使用训练样本
X
train
∈
R
n
×
d
X
train
∈
R
n
×
d
X_("train ")inR^(n xx d) X_{\text {train }} \in R^{n \times d} 训练一组
r
r
r r 个模型,得到一个聚合的异常分数矩阵
O
(
X
train
)
O
X
train
O(X_("train ")) O\left(X_{\text {train }}\right) 。在公式 (12) 中,
C
r
(
⋅
)
C
r
(
⋅
)
C_(r)(*) C_{r}(\cdot) 表示第
r
r
r r 个基检测器的分数向量。每个检测器的分数
C
r
(
X
train
)
C
r
X
train
C_(r)(X_("train ")) C_{r}\left(X_{\text {train }}\right) 通过 Z 分数标准化进行标准化(Aggarwal 和 Sathe 2015;Zimek、Campello 和 Sander 2014)。
O
(
X
train
)
=
[
C
1
(
X
train
)
,
…
,
C
r
(
X
train
)
]
O
X
train
=
C
1
X
train
,
…
,
C
r
X
train
O(X_("train "))=[C_(1)(X_("train ")),dots,C_(r)(X_("train "))] O\left(X_{\text {train }}\right)=\left[C_{1}\left(X_{\text {train }}\right), \ldots, C_{r}\left(X_{\text {train }}\right)\right]
(ii) 生成用于评估的伪地面真值。伪地面真值
O
(
X
train
)
O
X
train
O(X_("train ")) O\left(X_{\text {train }}\right) 是通过所有检测器的最大分数生成的(如原文所述)。这在公式 (13) 中进行了推广,其中
φ
φ
varphi \varphi 表示对所有基检测器的聚合。
target
=
φ
(
O
(
X
train
)
)
∈
R
n
×
1
target
=
φ
O
X
train
∈
R
n
×
1
" target "=varphi(O(X_("train ")))inR^(n xx1) \text { target }=\varphi\left(O\left(X_{\text {train }}\right)\right) \in R^{n \times 1}
(iii) 局部区域定义。测试实例
X
test
(
j
)
X
test
(
j
)
X_("test ")^((j)) X_{\text {test }}^{(j)} 的局部区域
ψ
j
ψ
j
psi_(j) \psi_{j} 定义为其 k 个最近训练对象的集合。形式上表示为公式 (14)。该过程涉及随机选择
t
t
t t 组
d
/
2
d
/
2
d//2 d / 2 到
d
d
d d 维的集合。 特征子空间。在每个子空间中,识别出距离
k
k
k k 最近的
X
test
(
j
)
X
test
(
j
)
X_("test ")^((j)) X_{\text {test }}^{(j)} 训练样本。然后,将在
k
N
N
ens
(
j
)
k
N
N
ens
(
j
)
kNN_("ens ")^((j)) k N N_{\text {ens }}^{(j)} 中出现次数超过
t
/
2
t
/
2
t//2 t / 2 的训练对象纳入,从而定义一个局部区域。按照原文,
k
k
k k 的值设定为训练样本的
10
%
10
%
10% 10 \% ,范围限定在[30;100]之间。
ψ
j
=
{
x
i
∣
x
i
∈
X
train
,
x
i
∈
k
N
N
ens
(
j
)
}
ψ
j
=
x
i
∣
x
i
∈
X
train
,
x
i
∈
k
N
N
ens
(
j
)
psi_(j)={x_(i)∣x_(i)inX_("train "),x_(i)in kNN_("ens ")^((j))} \psi_{j}=\left\{x_{i} \mid x_{i} \in X_{\text {train }}, x_{i} \in k N N_{\text {ens }}^{(j)}\right\}
(iv) 模型选择与组合。一旦确定了局部空间,利用训练好的基础检测器获得局部异常分数矩阵
O
(
ψ
j
)
O
ψ
j
O(psi_(j)) O\left(\psi_{j}\right) 。通过从目标数据集中提取与局部区域
j
j
j j 相关的值,获得局部伪真值目标
ψ
j
ψ
j
^(psi_(j)) { }^{\psi_{j}} 。LSCP 通过评估局部伪真值目标
ψ
j
ψ
j
^(psi_(j)) ^{\psi_{j}} 与局部检测器分数
C
r
(
X
train
ψ
j
)
C
r
X
train
ψ
j
C_(r)(X_("train ")^(psi_(j))) C_{r}\left(X_{\text {train }}^{\psi_{j}}\right) 之间的 Pearson 相关性,衡量每个基础检测器的局部能力(Schubert 等,2012)。最后,使用选定的
x
x
x x 个检测器,计算测试数据
X
test
(
j
)
X
test
(
j
)
X_("test ")^((j)) X_{\text {test }}^{(j)} 的异常分数,这些
x
x
x x 个异常分数的平均值作为测试数据的最终异常分数。
O
(
ψ
j
)
=
[
C
1
(
ψ
j
)
,
…
,
C
r
(
ψ
j
)
]
∈
R
|
ψ
j
|
×
R
O
ψ
j
=
C
1
ψ
j
,
…
,
C
r
ψ
j
∈
R
ψ
j
×
R
O(psi_(j))=[C_(1)(psi_(j)),dots,C_(r)(psi_(j))]inR^(|psi_(j)|xx R) O\left(\psi_{j}\right)=\left[C_{1}\left(\psi_{j}\right), \ldots, C_{r}\left(\psi_{j}\right)\right] \in R^{\left|\psi_{j}\right| \times R}
target
ψ
j
=
{
target
x
i
∣
x
i
∈
ψ
j
}
∈
R
|
ψ
j
|
×
1
target
ψ
j
=
target
x
i
∣
x
i
∈
ψ
j
∈
R
ψ
j
×
1
target^(psi_(j))={" target "_(x_(i))∣x_(i)inpsi_(j)}inR^(|psi_(j)|xx1) \operatorname{target}^{\psi_{j}}=\left\{\text { target }_{x_{i}} \mid x_{i} \in \psi_{j}\right\} \in R^{\left|\psi_{j}\right| \times 1}
通过将 LSCP 与各种传统的全局组合框架进行比较,评估了其性能。这些框架包括以下内容(Aggarwal 和 Sathe 2015):(1). 平均组合:基于每个基础检测器生成的平均得分,为每个数据点分配一个异常分数。(2). 最大组合:使用最大得分作为异常分数。(3). 最大平均(AOM)组合:将基础检测器随机划分为预定义的子集,最终得分通过计算每个子集中最大得分的平均值得出。(4). 平均最大(MOA)组合:最终 得分定义为每个子集中平均得分的最大值。
上述所有组合框架均利用相同的单个基础检测器池以确保一致性。较高的最终异常值分数表明更有可能是异常点。通常需要设置阈值以区分异常样本。本研究中使用了曲线下面积百分比(AUCP)算法来计算阈值。曲线下面积用于评估一种非参数方法,以对异常值分数生成的分数进行阈值划分。异常值被设定为核密度估计(KDE)的曲线下面积(AUC)小于总 KDE AUC 的(均值 + 绝对值(均值-中位数))百分比的任何值(Ren 等,2019)。曲线下面积(AUC)定义如下:
A
U
C
=
lim
x
→
inf
∑
i
=
1
n
f
(
x
)
δ
x
A
U
C
=
lim
x
→
inf
∑
i
=
1
n
f
(
x
)
δ
x
AUC=lim_(x rarr i n f)sum_(i=1)^(n)f(x)delta x A U C=\lim _{x \rightarrow \inf } \sum_{i=1}^{n} f(x) \delta x
f
(
x
)
f
(
x
)
f(x) f(x) 是曲线,
δ
x
δ
x
delta x \delta x 是矩形的增量步长,这些矩形的面积将被累加。AUCP 方法使用归一化决策分数的概率密度函数(pdf)在
0
−
1
0
−
1
0-1 0-1 范围内生成曲线。这是通过核密度估计完成的。
召回率和精确率被用来评估这些算法的效果(Ma 等,0000)。两个指标的计算结果如公式(18)和(19)所示:
Recall
=
T
P
T
P
+
F
N
Precision
=
T
P
T
P
+
F
P
Recall
=
T
P
T
P
+
F
N
Precision
=
T
P
T
P
+
F
P
{:[" Recall "=(TP)/(TP+FN)],[" Precision "=(TP)/(TP+FP)]:} \begin{gathered}
\text { Recall }=\frac{T P}{T P+F N} \\
\text { Precision }=\frac{T P}{T P+F P}
\end{gathered}
其中
T
P
T
P
TP T P 是被判断为异常样本的异常样本数量,
F
N
F
N
FN F N 是被判断为正常样本的异常样本数量,FP 是被判断为异常样本的正常样本数量。
4. 结果与讨论
4.1. 多源土地覆盖产品的精度评估
使用 SRS_Val 数据集验证了中国四种土地覆盖产品的分类精度,结果见表 2。ESA 和 CLCD 的总体精度较高,均超过 70%。DW 和 ESRI 的总体精度相对较低,其中 ESRI 最低,为 59.93%,kappa 系数为 0.523。
不同土地覆盖类别的准确性也存在显著差异。由于全国范围内验证样本分布相对稀疏,采用用户准确率(UA)作为不同类别的准确性指标。表 2 显示,树木覆盖、建筑区、裸地、水体以及冰雪等类别的分类准确率相对较高。然而,某些类别如灌木丛和湿地在所有产品中的分类准确率均不理想。草地仅在 ESA 和 CLCD 中具有较高的分类准确率。某些类别的低分类准确率以及不同产品间类别定义的差异可能影响所收集样本的可靠性(Venter 等,2022;H. Yang 等,2017)。因此,融合多源土地覆盖产品对于收集可靠的训练样本是必要的。
表2. 四种土地覆盖产品的准确性统计。
类别
DW
ESA
ESRI
CLCD
用户准确率(UA)(%)
PA(百分比)
UA(百分比)
PA(百分比)
UA(百分比)
PA(百分比)
UA(百分比)
PA(百分比)
树木覆盖率
89.39
70.72
91.53
75.25
86.89
80.21
91.46
79.79
灌木地
18.37
11.09
3.03
47.37
34.00
5.82
6.69
46.51
草地
15.86
77.78
69.71
59.13
8.40
75.76
79.50
64.09
农田
56.02
87.73
66.80
91.25
60.69
90.57
83.59
82.06
建成区
92.21
43.43
70.13
78.55
93.83
38.18
70.78
62.29
裸地
94.43
55.92
88.39
50.17
53.91
66.24
65.93
78.48
雪和冰
72.29
15.92
79.09
99.56
83.67
96.91
82.09
97.98
水体
89.30
61.52
86.33
77.78
89.33
65.53
84.56
81.29
湿地
18.58
94.83
33.22
89.09
29.67
82.41
8.33
100.00
总体准确率 (%)
62.62
73.22
59.93
卡帕系数
0.540
0.666
0.523
Class DW ESA ESRI CLCD
UA (%) PA (%) UA (%) PA (%) UA (%) PA (%) UA (%) PA (%)
Tree cover 89.39 70.72 91.53 75.25 86.89 80.21 91.46 79.79
Shrubland 18.37 11.09 3.03 47.37 34.00 5.82 6.69 46.51
Grassland 15.86 77.78 69.71 59.13 8.40 75.76 79.50 64.09
Cropland 56.02 87.73 66.80 91.25 60.69 90.57 83.59 82.06
Built area 92.21 43.43 70.13 78.55 93.83 38.18 70.78 62.29
Bare land 94.43 55.92 88.39 50.17 53.91 66.24 65.93 78.48
Snow and ice 72.29 15.92 79.09 99.56 83.67 96.91 82.09 97.98
Water 89.30 61.52 86.33 77.78 89.33 65.53 84.56 81.29
Wetland 18.58 94.83 33.22 89.09 29.67 82.41 8.33 100.00
OA (%) 62.62 73.22 59.93
Kappa 0.540 0.666 0.523 | Class | DW | | | ESA | | | ESRI | | | CLCD | |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| | UA (%) | | PA (%) | UA (%) | | PA (%) | UA (%) | | PA (%) | UA (%) | PA (%) |
| Tree cover | 89.39 | | 70.72 | 91.53 | | 75.25 | 86.89 | | 80.21 | 91.46 | 79.79 |
| Shrubland | 18.37 | | 11.09 | 3.03 | | 47.37 | 34.00 | | 5.82 | 6.69 | 46.51 |
| Grassland | 15.86 | | 77.78 | 69.71 | | 59.13 | 8.40 | | 75.76 | 79.50 | 64.09 |
| Cropland | 56.02 | | 87.73 | 66.80 | | 91.25 | 60.69 | | 90.57 | 83.59 | 82.06 |
| Built area | 92.21 | | 43.43 | 70.13 | | 78.55 | 93.83 | | 38.18 | 70.78 | 62.29 |
| Bare land | 94.43 | | 55.92 | 88.39 | | 50.17 | 53.91 | | 66.24 | 65.93 | 78.48 |
| Snow and ice | 72.29 | | 15.92 | 79.09 | | 99.56 | 83.67 | | 96.91 | 82.09 | 97.98 |
| Water | 89.30 | | 61.52 | 86.33 | | 77.78 | 89.33 | | 65.53 | 84.56 | 81.29 |
| Wetland | 18.58 | | 94.83 | 33.22 | | 89.09 | 29.67 | | 82.41 | 8.33 | 100.00 |
| OA (%) | | 62.62 | | | 73.22 | | | 59.93 | | | |
| Kappa | | 0.540 | | | 0.666 | | | 0.523 | | | |
OA:总体精度;UA:用户精度;PA:制图者精度。
图 5. 融合四个土地覆盖产品结果的置信度值。注意,每个像素的置信度值定义为加权多数投票算法计算出的最高加权票数。(a) 新疆北部地区;(b) 黑龙江省;(c) 广东省;(d) 京津冀地区。
4.2 多源土地覆盖产品的融合
4.2.1. 融合四种土地覆盖产品结果的置信度值
如图 4 所示,每个产品的权重值是基于用户准确率(UA)计算的。值得注意的是,使用此方法计算的权重值可能会导致某一类别的用户准确率较低,但其权重值却较高的情况。例如,尽管灌木地的用户准确率较低 在所有产品中,DW 和 ESRI 中的权重值却较高。应用加权多数投票(WMV)算法时,更有可能将分类结果归为灌木地。然而,这一问题并未显著影响最终收集的训练样本的可靠性。这些不确定区域将在收集初步训练样本之前,通过对 WMV 算法结果施加约束来移除。
由于事先选择了中国不同地理位置的四个示例区域,后续实验主要集中在这四个区域。如图 5 所示,每个像素的置信度值定义为 WMV 算法计算的最高加权投票。置信度值分为四个等级:低置信度(值小于 0.4,表示分类结果一致性低)、中低置信度(值在 0.4 到 0.7 之间)、中高置信度(值在 0.7 到 1.0 之间)和高置信度(值为 1,表示所有产品的分类结果完全一致)。图 5 显示,北疆的高置信度区域相对稀疏,主要包括一些南部、东部和北部地区。北疆许多地区被分类为中高置信度,而中低置信度区域更为集中,且倾向于形成条带状。在黑龙江省,高置信度区域主要位于西部和南部,中高置信度区域集中在北部。 黑龙江省北部和西部的一些地区表现出中低置信度。在广东省,除南部沿海等特定区域显示中低置信度外,大部分地区表现出高置信度。京津冀地区的高置信度区域主要集中在东南部,而北部和西部地区则主要表现为中高置信度值。
总体而言,四个区域中低置信度区域较少,大部分区域的置信度值高于0.7,表明宏观尺度上置信度较高。置信度值的空间分布模式与地理环境特征密切相关。农田、树木覆盖区、建筑区、水体和雪/冰在大多数区域表现出较高的置信度,置信度值达到1,表明分类结果可靠。这些区域的土地覆盖稳定,且在遥感影像中易于识别。相反,某些区域在灌木丛、草地、裸地和湿地等类别的置信度值相对较低,表明土地覆盖产品之间的分类结果不一致。这些区域低置信度值主要归因于光谱和纹理特征的相似性以及较高的... 这些类别的语义相似性(H. Wang 等,2022)。不同土地覆盖产品中草地-灌木丛和湿地-草地的定义差异导致分类结果存在混淆(Baig 等,2022)。此外,在具有复杂土地覆盖类别和地表高空间异质性的区域以及土地特征边缘区域,置信度值往往较低。这主要是由于像素识别的不确定性、类别定义的差异以及边缘效应(例如混合像素或区域)所致,导致不同土地覆盖产品之间空间一致性较差(Hao 等,2023;Radoux 等,2014)。置信度较低的区域会影响训练样本采集的质量。建议避免在这些区域内采集训练样本。此外,如果某一特定土地覆盖类别中存在较大比例的低置信度区域,可靠训练样本的数量也会受到影响。
4.2.2. 使用加权多数投票算法的土地覆盖图谱
通过 WMV 算法融合的分类结果对应于获得最高加权票数的土地覆盖(LC)类别。如图 6 所示,使用 WMV 算法生成的土地覆盖图表现出显著的区域差异。新疆北部地区以草地、裸地和农田为主,西部有艾丁湖(新疆最大的盐水湖)和赛里木湖。新疆北部的高海拔地区覆盖有大量的永久性积雪和冰川。黑龙江省以广泛的农田和林地为特征,西部地区集中分布着水体和湿地。广东省大部分地区以林地为主,拥有众多河流和湖泊。广东南部的珠江三角洲地区由于经济发达,建成区密度较高。在京津冀地区,农田最为广泛,主要集中在东南部,而林地和草地则占据北部和中西部。北京和天津的建成区尤为集中且广泛。
4.3. 初步训练样本的收集
局部自适应策略的初步训练样本收集结果如表3所示
图 6. 使用加权多数投票算法的土地覆盖图。(a) 新疆北部地区;(b) 黑龙江省;(c) 广东省;(d) 京津冀地区。 提取稳定区域后,在黑龙江省、广东省和京津冀地区收集了大量的树木覆盖、农田和建筑/不透水区域样本,因为这些类别在这些地区广泛存在。大部分草地和裸地样本收集于新疆北部地区。然而,一些类别的样本较少,尤其是雪和冰样本,仅在新疆北部收集。湿地样本在除三个地区外均相对稀少。
黑龙江省,那里样本更为丰富。在特定区域收集的灌木丛和草地样本相对较少。例如,在黑龙江省仅收集了8个灌木丛样本,在广东省收集了119个草地样本。这些样本稀缺的主要原因是这些类别的空间分布有限且分类结果的不确定性较高。此外,该方法在采样前提取了3月至11月的无云区域,这也产生了影响。
表3. 通过局部自适应策略在四个区域收集的各种土地覆盖类别的初步训练样本数量。
类别
北疆地区
黑龙江
广东
京津冀
树木覆盖率
2874
16274
6683
5571
灌木地
1129
8
1157
2443
草地
3797
2724
119
2953
耕地
2571
13764
5551
8607
建成区
2083
10915
5946
8510
裸地
5689
2472
2439
2339
雪和冰
1406
0
0
0
水体
1352
7472
4949
4070
湿地
386
2700
251
701
总计
21287
56329
27095
35194
Class Northern Xinjiang Heilongjiang Guangdong Beijing-Tianjin-Hebei
Tree cover 2874 16274 6683 5571
Shrubland 1129 8 1157 2443
Grassland 3797 2724 119 2953
Cropland 2571 13764 5551 8607
Built area 2083 10915 5946 8510
Bare land 5689 2472 2439 2339
Snow and ice 1406 0 0 0
Water 1352 7472 4949 4070
Wetland 386 2700 251 701
Total 21287 56329 27095 35194 | Class | Northern Xinjiang | Heilongjiang | Guangdong | Beijing-Tianjin-Hebei |
| :--- | :--- | :--- | :--- | :--- |
| Tree cover | 2874 | 16274 | 6683 | 5571 |
| Shrubland | 1129 | 8 | 1157 | 2443 |
| Grassland | 3797 | 2724 | 119 | 2953 |
| Cropland | 2571 | 13764 | 5551 | 8607 |
| Built area | 2083 | 10915 | 5946 | 8510 |
| Bare land | 5689 | 2472 | 2439 | 2339 |
| Snow and ice | 1406 | 0 | 0 | 0 |
| Water | 1352 | 7472 | 4949 | 4070 |
| Wetland | 386 | 2700 | 251 | 701 |
| Total | 21287 | 56329 | 27095 | 35194 |
收集的样本数量。四个研究区全年受云覆盖掩膜的影响较小,减少了对样本收集数量的影响。然而,在许多其他地区,尤其是中国南方,夏季云量较大,可能导致合成的月度卫星图像出现大量数据缺失。为了解决这个问题,最佳方案是排除数据缺失严重月份的卫星图像。仅利用数据缺失较少月份的合成图像,并提取无云区域进行初步样本收集。使用收集的训练样本进行土地覆盖分类时也适用同样的方法。
训练样本的质量和数量对于训练分类模型至关重要(Mellor 等,2015)。本文描述的样本收集策略能够收集可靠且均匀分布的初步训练样本。然而,初步样本的数量可能因研究区域的不同而显著变化,这可能导致样本不平衡。样本不平衡是指某一类别或多个类别的训练样本数量显著多于或少于其他类别。这可能导致稀有地类相对于丰富类别的代表性不足,从而降低整体分类精度(Estabrooks, Jo, 和 Japkowicz 2004;Mellor 等,2015)。已有研究探索了多数类的下采样(Freeman, Moisen, 和 Frescino 2012)和少数类的过采样(Ling 和 Li 1998)等技术,以缓解训练样本不平衡的问题。因此,用户在应用我们的样本收集策略时,可能需要根据具体研究区域和分类模型进行调整。 例如,用户可以增加每个网格中样本较少类别的采样数量,或调整网格大小(本文设定为5公里)以 增加或减少某些类别的采样数量。
通过对四个区域中随机选取的 6671 个样本进行目视解译,利用 Google Earth 和 Bing Maps 等在线卫星影像验证了所收集初步训练样本的准确性。每个区域随机选取的样本数量在 1606 至 1739 之间。还提取并分析了 2020 年的时间序列光谱曲线以辅助解译。验证结果显示,北疆、黑龙江、广东及京津冀地区初步训练样本的准确率分别为
96.26
%
,
97.01
%
,
94.24
%
96.26
%
,
97.01
%
,
94.24
%
96.26%,97.01%,94.24% 96.26 \%, 97.01 \%, 94.24 \% 和
94.97
%
94.97
%
94.97% 94.97 \% 。四个区域初步训练样本的平均准确率为
95.62
%
95.62
%
95.62% 95.62 \% ,表明整体质量令人满意。
4.4. 通过 LSCP 算法过滤异常样本
4.4.1. 异常样本类型
尽管收集的初步训练样本具有较高的准确性,但仍存在错误和异常样本。异常样本可分为三类:(1)分类错误(标签错误),(2)一年内的土地覆盖变化,以及(3)光谱异常(也称为特征异常,包括混合像元和由环境因素或卫星数据质量问题引起的异常)。
表 4 列出了六种典型异常样本的异常类型、描述、图像和时间序列光谱曲线。样本(a)和(b)代表由分类错误(标签错误)引起的异常。在光谱曲线图中,橙色曲线表示异常样本的光谱曲线,蓝色曲线表示正常样本的光谱曲线。
表4. 六种典型异常样本的异常类型、描述、图像及时间序列光谱曲线。
异常类型
描述
图像
时间序列光谱曲线
分类错误/标签错误
样本 (a):树木被误分类为建筑区。
1 0.2 0 3 4 5 6 7 月 8 9 10 11
1
0.2
0
3
4
5
6
7 Mont
8
9
10
11 | 1 |
| :--- |
| 0.2 |
| 0 |
| 3 |
| 4 |
| 5 |
| 6 |
| 7 Mont |
| 8 |
| 9 |
| 10 |
| 11 |
样本 (b):裸地被误分类为草地。
月份
一年内的土地覆盖变化
样本 (c):水体的季节性变化。
1 个月
光谱异常
样本(d):混合像元。
0.8
0
3
4
5
6
7 8
9
10
I | 0.8 |
| :--- |
| 0 |
| 3 |
| 4 |
| 5 |
| 6 |
| 7 8 |
| 9 |
| 10 |
| I |
样本(e):由环境因素或卫星数据质量问题引起的异常情况
Sample
(e): Anomaly caused by environmental factors or issues with satellite data quality | Sample |
| :--- |
| (e): Anomaly caused by environmental factors or issues with satellite data quality |
0.6 0.4 0.2
−
0.2
−
0.2
-0.2 -0.2
−
0.4
−
0.4
-0.4 -0.4
−
0.6
−
0.6
-0.6 -0.6
−
0.8
−
0.8
-0.8 -0.8 月份 1 0.8 0
0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
Month
1
0.8
0 | 0.6 |
| :--- |
| 0.4 |
| 0.2 |
| $-0.2$ |
| $-0.4$ |
| $-0.6$ |
| $-0.8$ |
| Month |
| 1 |
| 0.8 |
| 0 |
Anomaly type Description Image Time-series spectral curve
Classification error/Label error Sample (a): Trees misclassified as built area. https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=352&width=351&top_left_y=261&top_left_x=850 "1
0.2
0
3
4
5
6
7 Mont
8
9
10
11"
Sample (b): Bare land misclassified as grassland. https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=352&width=350&top_left_y=659&top_left_x=852 Month
LC changes within a year Sample (c): Seasonal changes in water. https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=349&width=345&top_left_y=1058&top_left_x=853 https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=373&width=682&top_left_y=1040&top_left_x=1211
1 Month
Spectral anomalies Sample (d): Mixed pixels. https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=352&width=350&top_left_y=1451&top_left_x=852 "0.8
0
3
4
5
6
7 8
9
10
I"
"Sample
(e): Anomaly caused by environmental factors or issues with satellite data quality" https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-17.jpg?height=357&width=347&top_left_y=1848&top_left_x=852 "0.6
0.4
0.2
-0.2
-0.4
-0.6
-0.8
Month
1
0.8
0" | Anomaly type | Description | Image | Time-series spectral curve |
| :--- | :--- | :--- | :--- |
| Classification error/Label error | Sample (a): Trees misclassified as built area. |  | 1 <br> 0.2 <br> 0 <br> 3 <br> 4 <br> 5 <br> 6 <br> 7 Mont <br> 8 <br> 9 <br> 10 <br> 11 |
| | Sample (b): Bare land misclassified as grassland. |  | Month |
| LC changes within a year | Sample (c): Seasonal changes in water. |  |  |
| | | | 1 Month |
| Spectral anomalies | Sample (d): Mixed pixels. |  | 0.8 <br> 0 <br> 3 <br> 4 <br> 5 <br> 6 <br> 7 8 <br> 9 <br> 10 <br> I |
| | Sample <br> (e): Anomaly caused by environmental factors or issues with satellite data quality |  | 0.6 <br> 0.4 <br> 0.2 <br> $-0.2$ <br> $-0.4$ <br> $-0.6$ <br> $-0.8$ <br> Month <br> 1 <br> 0.8 <br> 0 |
表4.(续)
异常类型
描述
图像
时间序列光谱曲线
样本(f):由环境因素或卫星数据质量问题引起的异常。
⟶
⟶
longrightarrow \longrightarrow 蓝 绿 红
=
=
= = 近红外
https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-18.jpg?height=364&width=700&top_left_y=268&top_left_x=1207
longrightarrow Blue Green Red = NIR |  |
| :--- |
| $\longrightarrow$ Blue Green Red $=$ NIR |
Anomaly type Description Image Time-series spectral curve
Sample (f): Anomaly caused by environmental factors or issues with satellite data quality. "https://cdn.mathpix.com/cropped/2025_05_03_cbd0180b65abdc269174g-18.jpg?height=364&width=700&top_left_y=268&top_left_x=1207
longrightarrow Blue Green Red = NIR" | Anomaly type | Description | Image | Time-series spectral curve | |
| :--- | :--- | :--- | :--- | :--- |
| | Sample (f): Anomaly caused by environmental factors or issues with satellite data quality. | |  <br> $\longrightarrow$ Blue Green Red $=$ NIR | |
样本。NDVI 时间序列曲线的差异使得异常样本得以区分。样本©显示了由一年内 LC 类别变化引起的异常。从六月到八月,样本©的 NDVI 值为正,而 NDWI 值为负。主要原因可能是该地区的季节性水文变化和潮汐影响,导致水位变化(Y. Li 等,2019;Yin 等,2023)。这些因素使得原本被水覆盖的区域暴露出来,加上夏季植被旺盛生长,导致 NDVI 显著增加。样本(d)代表由卫星影像中的混合像元引起的异常(X. Liu, Li, 和 Zhang 2010)。该样本位于建筑物边缘,卫星影像中的像元包含来自多个 LC 类别的光谱信息,如建筑物和树木。因此,其光谱曲线与纯像元不同。样本(e)和(f)显示了由环境因素或卫星数据质量问题引起的异常。 光谱曲线显示,样本(e)在三月份出现极端异常,绿光和近红外波段的值低于其他月份。NDVI 为负值,NDWI 甚至超过 0.4。样本(f)在十月份出现光谱异常, 可见光和近红外波段的值接近零(0.0001),导致计算出的光谱指数如 EVI、NDVI 和 NDWI 均为零。
训练样本中存在异常样本(噪声)会对模型训练产生不利影响,进而影响土地覆盖分类的准确性。在训练过程中,模型努力最小化损失函数,即尝试减少预测值与真实值之间的差异。 训练样本中的真实标签与模型预测结果之间的关系。如果训练样本中存在噪声,模型可能也会试图拟合这些噪声,从而导致过拟合。以往研究者已经评估了噪声对分类准确率的影响。Rodriguez-Galiano 等人(2012)报告称,随机森林分类器对训练样本中高达 20%的故意错误标注表现出相对不敏感,超过此比例后错误率呈指数增长。Mellor 等人(2015)发现,随着噪声比例的增加,分类准确率逐渐下降。此外,某些分类算法的学习模型复杂度可能会受到影响(Pelletier 等人,2017)。例如,在存在类别标签噪声的情况下,训练实例的平均路径长度可能增加,导致计算训练时间增加。因此,去除初步训练样本中存在的异常样本至关重要。
4.4.2. 不同异常检测模型的比较
利用时间序列光谱特征在 Python 中实现了用于异常检测的集成框架,以识别异常样本。异常检测模型针对四个区域的每个类别分别训练和应用。由于黑龙江省灌木地样本稀缺,故预先将其移除。所有集成框架均使用了由 30 个基于 LOF 的基检测器组成的池,以保证性能评估的一致性。为了在基检测器之间引入多样性(Britto, Sabourin, 和 Oliveira 2014;Zimek, Campello, 和 Sander 2014),采用了不同的初始化方法。
表5. 四个区域中不同异常检测模型性能比较(%)。
模型
平均值
最大值
AOM
MOA
LSCP
新疆北部
召回率
71.67
63.33
66.67
68.33
75.00
精确率
36.44
36.89
37.39
35.65
35.16
去除率
7.84
6.82
7.15
7.65
8.37
黑龙江
召回率
69.23
67.31
69.23
67.31
69.23
精确率
50.70
49.30
50.70
50.72
49.32
去除率
2.94
2.69
2.74
2.83
3.01
广东
召回率
76.34
63.44
68.82
75.27
76.34
精确率
47.97
49.58
50.00
48.61
47.97
去除率
4.74
3.84
4.05
4.57
4.84
京津冀
召回率
55.81
46.51
54.65
55.81
60.47
精确率
44.04
43.01
48.96
44.44
44.83
去除率
5.05
4.53
4.60
4.96
5.29
Models Average Max AOM MOA LSCP
Northern Xinjiang Recall 71.67 63.33 66.67 68.33 75.00
Precision 36.44 36.89 37.39 35.65 35.16
Removal ratio 7.84 6.82 7.15 7.65 8.37
Heilongjiang Recall 69.23 67.31 69.23 67.31 69.23
Precision 50.70 49.30 50.70 50.72 49.32
Removal ratio 2.94 2.69 2.74 2.83 3.01
Guangdong Recall 76.34 63.44 68.82 75.27 76.34
Precision 47.97 49.58 50.00 48.61 47.97
Removal ratio 4.74 3.84 4.05 4.57 4.84
Beijing-Tianjin-Hebei Recall 55.81 46.51 54.65 55.81 60.47
Precision 44.04 43.01 48.96 44.44 44.83
Removal ratio 5.05 4.53 4.60 4.96 5.29 | | Models | Average | Max | AOM | MOA | LSCP |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| Northern Xinjiang | Recall | 71.67 | 63.33 | 66.67 | 68.33 | 75.00 |
| | Precision | 36.44 | 36.89 | 37.39 | 35.65 | 35.16 |
| | Removal ratio | 7.84 | 6.82 | 7.15 | 7.65 | 8.37 |
| Heilongjiang | Recall | 69.23 | 67.31 | 69.23 | 67.31 | 69.23 |
| | Precision | 50.70 | 49.30 | 50.70 | 50.72 | 49.32 |
| | Removal ratio | 2.94 | 2.69 | 2.74 | 2.83 | 3.01 |
| Guangdong | Recall | 76.34 | 63.44 | 68.82 | 75.27 | 76.34 |
| | Precision | 47.97 | 49.58 | 50.00 | 48.61 | 47.97 |
| | Removal ratio | 4.74 | 3.84 | 4.05 | 4.57 | 4.84 |
| Beijing-Tianjin-Hebei | Recall | 55.81 | 46.51 | 54.65 | 55.81 | 60.47 |
| | Precision | 44.04 | 43.01 | 48.96 | 44.44 | 44.83 |
| | Removal ratio | 5.05 | 4.53 | 4.60 | 4.96 | 5.29 |
去除率指的是模型标记为异常的初步训练样本所占的比例。
表 6. LSCP 与传统统计方法性能比较(%)。
方法
召回率
精确率
去除率
参考文献
PCA-Lajda 准则
74.91
7.54
43.95
Jin 等人,(2022)
标准差滤波(NDVI 时间序列)
97.59
5.27
84.24
Wen 等人,(2022)
LSCP
70.10
44.93
4.75
Approaches Recall Precision Removal ratio Reference
PCA-Lajda criterion 74.91 7.54 43.95 Jin et al., (2022)
Standard deviation filtering (NDVI time-series) 97.59 5.27 84.24 Wen et al., (2022)
LSCP 70.10 44.93 4.75 | Approaches | Recall | Precision | Removal ratio | Reference |
| :--- | :--- | :--- | :--- | :--- |
| PCA-Lajda criterion | 74.91 | 7.54 | 43.95 | Jin et al., (2022) |
| Standard deviation filtering (NDVI time-series) | 97.59 | 5.27 | 84.24 | Wen et al., (2022) |
| LSCP | 70.10 | 44.93 | 4.75 | |
超参数,即每个 LOF 检测器中使用的邻居数量(MinPts),在
[
5
,
150
]
[
5
,
150
]
[5,150] [5,150] 范围内选择。对于 AOM 和 MOA 框架,基础检测器被分为五个子组;每组包含六个不重复选择的基础检测器。
经过目视解译的训练样本被用来评估不同模型的性能。LSCP 和其他模型的结果如表 5 所示。召回率表示被正确识别的异常样本的比例。与传统模型相比,LSCP 模型表现出更优越且稳健的性能,召回率普遍较高。LSCP 的优越性归因于其能够仅整合测试实例局部区域内表现良好的基础检测器,从而减轻表现不佳检测器的影响。其优势在新疆北部和京津冀地区尤为明显。然而,不同区域的召回率存在显著差异。新疆北部和广东的召回率超过 75%,而京津冀地区的召回率最低,仅为
60.47
%
60.47
%
60.47% 60.47 \% 。因此,异常检测模型的性能在不同研究区域间差异显著。表 5 显示,所有异常检测模型的精确率均较低,表明许多正常样本被误判为异常。 这种现象可能归因于这些样本在特征空间中的分布偏离了 大多数正常样本,尽管它们的标签正确且在视觉解译过程中未发现明显异常。LSCP 模型在北疆、黑龙江、广东和京津冀初步训练样本中检测到的异常比例(剔除率)分别为
8.37
%
,
3.01
%
,
4.84
%
8.37
%
,
3.01
%
,
4.84
%
8.37%,3.01%,4.84% 8.37 \%, 3.01 \%, 4.84 \% 和
5.29
%
5.29
%
5.29% 5.29 \% 。虽然使用异常检测模型可能会剔除一些正常样本,但相对于样本总数的剔除比例非常小,因此整体影响有限。在 LSCP 模型剔除异常样本后,北疆、黑龙江、广东和京津冀地区最终训练样本的准确率提高了
2.05
%
2.05
%
2.05% 2.05 \% 至 4.34%,分别达到 99.04%、99.06%、98.58%和 97.95%。即使经过样本过滤,最终训练样本中仍存在错误和不确定性。然而,最终样本中的错误较小,足以用于模型训练。 值得注意的是,用于评估异常检测算法的视觉解译样本数量可能有限,且异常样本在整体样本集中的比例特别小。因此,本研究的验证结果可能与实际准确率存在偏差。尽管如此,这些结果仍对读者具有一定的参考价值。
分别计算了三种异常类型的总数和召回率,以分析异常检测模型对不同类型异常的检测效果。总计291个。 在视觉解译过程中,四个区域均发现异常样本。共有 170 个样本存在分类错误(标签错误),33 个样本在一年内发生了土地覆盖变化,88 个样本存在光谱异常(特征异常)。使用 LSCP 模型进行异常检测时,分类错误、年度土地覆盖变化和光谱异常的召回率分别为
67.65
%
,
69.70
%
67.65
%
,
69.70
%
67.65%,69.70% 67.65 \%, 69.70 \% 和
75.00
%
75.00
%
75.00% 75.00 \% 。LSCP 模型对不同类型异常样本的识别准确率有所不同。识别由分类错误引起的异常的准确率最低。某些植被类型,尤其是灌木丛、草地和湿地,由于其光谱和纹理特征相似,以及土地覆盖产品中存在的混淆,容易被误分类。因此,许多带有分类错误的异常样本也难以被正确识别。识别一年内土地覆盖变化异常样本以及光谱异常样本的准确率相对较高。 对于具有光谱异常的样本,一个或多个时间段内的异常光谱值导致其时间序列光谱曲线与其他样本显著不同,使其更易被检测到。此外,由于土地覆盖变化是一个缓慢的过程,一年内发生土地覆盖变化的样本数量相对较少,主要与水体的季节性变化有关(如表4中样本©所示)。异常检测算法在一年内发生土地覆盖变化的异常样本上的有效性需要进一步验证。还应注意,异常检测算法中的阈值是基于所有点的异常分数计算得出的,只有异常分数超过阈值的样本才被分类为异常。由于不同算法在阈值计算上的差异,训练样本的剔除比例和异常检测的准确性也会有所不同。总体而言,基于集成机器学习的异常检测算法在过滤异常样本方面是有效的。 未来的检测性能有望通过异常检测算法和集成策略的改进以及其他典型特征的选择得到进一步提升。
4.4.3. 与传统方法的比较
通过将 LSCP 模型的结果与传统统计技术的结果进行比较,评估了不同方法的性能。
具体而言,采用了 Jin 等人(2022)和 Wen 等人(2022)用于异常值去除的方法。两种方法均依赖于数据的正态分布,去除落在
μ
±
σ
μ
±
σ
mu+-sigma \mu \pm \sigma 范围之外的训练样本。Jin 等人进行了主成分分析(PCA)以计算 PC1 和 PC2,利用这些主成分过滤异常样本。Wen 等人利用 NDVI 时间序列数据从平均 NDVI 曲线中剔除异常样本。其他方法,如 Zhang 等人(2021)提出的方法,也基于光谱统计分布去除异常样本。由于这些方法原理相似,因此未单独进行比较。
召回率、精确率和去除率被用来比较不同的方法,这些评估指标基于来自四个区域的所有样本计算。表 6 显示,LSCP 的召回率为
70.10
%
70.10
%
70.10% 70.10 \% ,低于两种统计方法的召回率。特别是,标准差过滤方法达到了
97.59
%
97.59
%
97.59% 97.59 \% 的召回率,表明其几乎能够识别所有异常样本。然而,两种统计方法的精确率显著低于 LSCP,而去除率则明显更高。尽管传统统计方法能够正确识别大多数异常样本,但也去除了大量正常样本。事实上,视觉解译识别的异常样本比例很小,表明仅依赖统计方法可能导致训练样本被过度去除。此外,由于数据维度的降低,PCA-Lajda 准则方法可能忽略许多在一年内表现出土地覆盖变化和光谱异常的样本。 原文采用标准差滤波方法去除异常的玉米样本。在多类别土地覆盖训练样本中使用时,由于同一类别的时间序列光谱特征多样(例如不同作物表现出不同的光谱曲线),大量样本被剔除。需要注意的是,除了样本质量,样本数量对于模型训练同样至关重要。尽管 LSCP 的召回率略低于传统方法,但它更好地平衡了样本质量和数量,避免了过度剔除训练样本。基于机器学习的异常检测方法如 LSCP 不受数据分布影响,并能提供异常分数。
表7. 单时相与时间序列光谱特征在六个典型异常样本上的模型性能比较。
模型
(a)
(b)
(c)
(d)
(e)
(f)
总计
LSCP_4
1
1
0
1
0
0
3
LSCP_7
1
1
1
1
0
0
4
LSCP_10
1
0
0
1
0
1
3
LSCP_时间序列
1
1
1
1
1
1
6
Model (a) (b) (c) (d) (e) (f) Total
LSCP_4 1 1 0 1 0 0 3
LSCP_7 1 1 1 1 0 0 4
LSCP_10 1 0 0 1 0 1 3
LSCP _TS 1 1 1 1 1 1 6 | Model | (a) | (b) | (c) | (d) | (e) | (f) | Total |
| :--- | :--- | :--- | :--- | :--- | :--- | :--- | :--- |
| LSCP_4 | 1 | 1 | 0 | 1 | 0 | 0 | 3 |
| LSCP_7 | 1 | 1 | 1 | 1 | 0 | 0 | 4 |
| LSCP_10 | 1 | 0 | 0 | 1 | 0 | 1 | 3 |
| LSCP _TS | 1 | 1 | 1 | 1 | 1 | 1 | 6 |
LSCP_4:利用四月特征的 LSCP 模型;LSCP_7:利用七月特征的 LSCP 模型;LSCP_10:利用十月特征的 LSCP 模型;LSCP_TS:利用时间序列特征的 LSCP 模型。 对每个样本,帮助用户理解异常程度。因此,基于机器学习的异常检测算法在过滤异常土地覆盖训练样本方面具有广阔的应用前景。
4.4.4. 利用单时相光谱特征模型的比较
使用表现最佳的集成框架 LSCP,比较了基于单时相和时间序列光谱特征的模型对各种异常的有效性。时间序列特征涵盖了从三月到十一月的光谱特征。单时相光谱特征则从四月、七月和十月的合成影像中提取,分别代表春季、夏季和秋季的光谱特征。使用表 4 中的六个典型异常样本验证了不同时间特征下的检测结果。在表 7 中,若模型成功识别异常样本,则标记为“1”,未识别则标记为“0”。“总计”栏显示成功识别的异常样本总数。
表 7 显示,所有模型都能有效识别由分类错误(例如样本(a)和(b))和混合像元(例如样本(d))引起的异常。仅利用单时相特征准确识别异常样本仍存在挑战。例如,样本(b)在仅依赖十月份的光谱特征时未被识别为异常样本。利用单时相特征的模型在检测一年内土地覆盖变化引起的异常(例如样本©)和光谱异常(例如样本(e)和(f))方面表现不足。以样本©为例,水体的季节性和水位变化主要在夏季显现,因此检测依赖于该季节的光谱特征。样本(e)和(f)也遵循同样的原则。因此,采用基于时间序列光谱特征的异常检测被认为更为合适,并且能够带来 最高的准确率,表现出在检测各种类型异常样本方面的强大性能。
4.5. ATSC 方法的泛化能力
尽管本研究的实验范围仅限于中国,但作为一个整体,中国地域辽阔,气候和地理条件在各地区之间存在显著差异。所选的四个实验区位于中国不同的地理位置,气候、土地利用和土地覆盖、地形以及社会经济发展水平均有显著差异。利用 ATSC 方法获得的土地覆盖训练样本在这些地理多样的区域中表现出理想的准确性,表明该方法具有较强的泛化能力。随着全球和区域尺度高分辨率土地覆盖产品的日益普及,以及 Sentinel-2 等卫星影像的全球覆盖,我们预见 ATSC 方法在跨区域应用中的可行性和潜力。
在将 ATSC 方法应用于其他区域时,需要考虑几个关键因素。首先,选择具有适当时间和空间覆盖的合适 LC 产品至关重要。可以选择覆盖目标区域的全球尺度 LC 产品或区域尺度产品。其次,应用 WMV 算法时,需要覆盖目标区域的验证数据集以计算不同产品的类别权重值。验证数据集可以利用公开可用的数据集(如 SRS_Val)或通过目视解译创建。为了尽量减少时间差异的影响,必须确保所选 LC 产品和验证数据集的参考年份尽可能接近。第三,在采用局部自适应采样策略时,需要合理设置网格大小和 采样数量基于具体的研究区域和应用需求。最后,可能需要对像 LSCP 这样的异常检测模型进行参数调整,以有效适应具体环境。此外,必须注意与该方法相关的不确定性。所收集训练样本的可靠性受 LC 产品的数量和准确性以及验证数据集质量等因素的影响。此外,异常检测算法的准确性可能会受到某些区域 LC 光谱复杂性和变异性的影响。
然而,尚未在全球其他地区进行进一步验证。必须承认,当将 ATSC 方法应用于地理环境显著不同的地区时,例如气候与中国截然不同的热带雨林,其性能可能会有所不同。未来的研究应包括在多样的地理和气候条件下进行广泛的案例研究,以评估 ATSC 方法的泛化能力和适用性。
5. 结论
本文提出了一种用于自动收集土地覆盖分类训练样本的新方法。该方法首先使用统一的验证数据集评估多源土地覆盖产品的准确性。根据用户准确率(UA)计算类别权重值,并利用加权多数投票(WMV)算法融合多种土地覆盖产品。然后从融合的土地覆盖图中提取高置信度区域,并通过形态学处理去除边缘区域的像素。随后,利用时间序列 Sentinel-2 影像提取无云区域。采用局部自适应策略自动收集高质量训练样本。验证结果表明,所收集的初步训练样本在四个研究区的准确率范围为
94.24
%
94.24
%
94.24% 94.24 \% 至
97.01
%
97.01
%
97.01% 97.01 \% ,平均准确率为
95.62
%
95.62
%
95.62% 95.62 \% 。随后,采用基于集成的异常检测算法进一步筛选异常样本。利用人工解译样本进行验证,显示不同异常检测模型在识别土地覆盖初步样本数据集中的异常样本方面的有效性。 基于机器学习的异常检测算法,如 LSCP,被发现比传统统计方法更好地平衡了训练样本的质量和数量。其中, 利用时间序列光谱特征的 LSCP 模型表现最佳。结果表明,在去除异常样本后,所有区域最终训练样本的准确率进一步提高至
97.95
%
97.95
%
97.95% 97.95 \% 或更高。
与以往方法相比,ATSC 方法具有多项优势。(1) 训练样本的收集完全自动化,消除了为大规模土地覆盖分类生成足够训练样本的繁重劳动过程。该方法适用于不同地理位置和环境的区域,以及不同年份的应用。(2) 融合了多种现有的高空间分辨率土地覆盖产品,使收集的训练样本更适合使用高空间分辨率卫星影像(如 Sentinel-2)进行土地覆盖分类。(3) 最终的训练样本数据集包含两个指标——置信值和异常值评分,用于衡量样本的可靠性。用户可以综合利用这两个指标,在应用于不同区域和不同分类算法时进一步筛选训练样本。因此,通过我们的方法收集的可靠训练样本可用于高质量的土地覆盖制图。
披露声明
作者未报告任何潜在的利益冲突。
资金
本工作部分得到了北京杰出青年科学家计划 [BJJWZYJH01201910028032] 及国家重点研发计划 [2018YFC1508902, 2017YFC0406006, 2017YFC0406004] 的支持。
ORCID
作者贡献声明
Yanzhao Wang:构思与设计,数据分析与解释,论文撰写。Yonghua Sun:构思与设计,批判性修订以提升学术内容,最终批准发表版本。Xuyue Cao:构思与设计,论文撰写。Yihan Wang:数据分析与解释,论文撰写。Wangkuan Zhang:数据分析与解释,论文撰写。Xinglu Cheng:分析与