多传感器时空融合技术用于30米分辨率的日连续雪盖地图绘制

《Water Resources Research》：Multi-Sensor Spatiotemporal Fusion for 30-m Daily Gapless Snow Cover Mapping

　　摘要

高时空分辨率的遥感数据对于监测异质性山区雪盖至关重要。尽管时空融合为高分辨率雪盖监测提供了一种有前景的方法，但云层污染和观测数据稀疏仍然限制了其在大规模和长期应用中的有效性。为了解决这个问题，我们提出了一种自适应时间序列融合框架，基于增强型空间和时间自适应反射率融合模型（ESTARFM）生成30米分辨率的每日连续雪盖数据。我们整合了多源粗分辨率数据和多源细分辨率数据，以增加有效像素的数量并提高数据密度，从而捕捉雪盖的快速时空变化。此外，我们引入了时间序列图像对来调整ESTARFM方法，该方法克服了细分辨率数据的空间完整性限制，并为每个像素动态选择最接近目标时间的时间信息。综合评估证实了融合结果的高准确性，这与参考数据的一致性（R = 0.776–0.964）证明了这一点。进一步地，基于地面雪盖观测的验证表明，30米分辨率的每日雪盖数据在捕捉雪盖时间动态方面优于广泛使用的500米数据，这一点通过其与地面雪盖物候指标的强一致性得到了证实。此外，该数据还提供了关于山区雪盖空间分布的新见解。在海拔低于3,500米、坡度小于25度或处于阴影区域的地区，30米数据能够捕捉到小规模、稀疏且碎片化的雪盖，这对需要准确雪盖估计的水文研究和实际应用具有重大潜力。

**通俗语言总结**

高时空分辨率的遥感数据对于监测山区雪盖至关重要，但云层污染和观测数据稀疏限制了时空融合技术在大规模和长期应用中的有效性。为了解决这个问题，我们提出了一种自适应时间序列融合框架，该框架整合了多源数据并引入了时间序列粗细图像对。这种方法能够捕捉雪盖的快速时空变化，并克服了输入细分辨率数据在时空融合中的空间完整性限制。值得注意的是，我们融合的30米分辨率每日雪盖数据在捕捉雪盖时间动态方面优于广泛使用的500米分辨率每日雪盖数据，这一点通过与地面雪盖物候指标的强一致性得到了证实。它还提供了关于山区雪盖空间分布的新见解。在海拔低于3,500米、坡度小于25度或处于阴影区域的地区，30米数据能够捕捉到更小、更稀疏且更碎片化的雪盖斑块。这些数据对于依赖准确雪盖估计的水文研究和实际应用具有重大潜力。

**1 引言**

山区雪盖及其融化是全球许多地区的重要水源，特别是在干旱和半干旱地区（Chen等人，2022年；Painter等人，2016年）。然而，雪盖表现出显著的时空异质性，空间分布碎片化且变化迅速（Bormann等人，2018年；Flanner等人，2011年），这使得准确估计雪盖变得困难。这些估计雪盖的困难体现在对细尺度变化的量化上，并可能对水资源管理、农业灌溉等依赖准确雪盖估计的领域产生重大影响（Henderson等人，2018年；Immerzeel等人，2010年）。遥感技术对于雪盖监测至关重要，因为它能够提供长期和大范围的观测（Dong，2018年）。粗分辨率的遥感卫星传感器，如中分辨率成像光谱辐射计（MODIS）、可见光红外成像辐射计套件（VIIRS）和高级非常高分辨率辐射计（AVHRR），提供了百米或千米分辨率的近乎每日覆盖，能够检测到雪盖的快速变化（Huang等人，2022年）。然而，山区雪盖通常在更细的尺度上发生变化，需要更高空间分辨率的数据来探索其时空动态。Landsat（30米，16天周期）和Sentinel-2（10米、20米、60米，10天周期）被广泛用于十米级雪盖制图（Gascoin等人，2019年），但它们的时间分辨率超过一天，不足以有效监测雪盖的快速变化（Slater等人，2013年）。传感器在时间和空间分辨率之间的权衡为获得高时空分辨率雪盖数据带来了挑战。对高时空分辨率雪盖数据的需求推动了数据融合技术的发展（Mityók等人，2018年；Premier等人，2021年；Zhang & Jiang，2022年）。Berman等人（2018年）采用了一种动态时间扭曲技术，将历史Landsat观测与MODIS时间序列对齐，从而生成了30米分辨率的每日雪盖数据。然而，该方法仅考虑无云数据，难以捕捉短期雪盖变化事件。Zakeri和Mariethoz（2024年）提出了一种基于假设的方法，即具有相似气候特征的年份雪模式会重复出现，使用K最近邻算法建立时间序列气象预测因子与可用Landsat/Sentinel-2雪盖图像之间的关系，从而生成30米分辨率的每日雪盖地图。与Zakeri和Mariethoz（2024年）类似，一些研究利用了辅助信息，如气象和地形因素，来降尺度粗分辨率MODIS雪盖数据，以提高估计精度（Li等人，2015年；Mahanthege等人，2024年；Richiardi等人，2023年；Walters等人，2014年）。例如，Rittger等人（2021年）开发了一种两阶段随机森林算法，通过结合辅助数据（包括地形预测因子和地球物理特征）将MODIS雪盖分数（FSC）降尺度到30米分辨率，从而实现了更清晰的边缘和更详细的空间分布。然而，机器学习通常依赖于高质量的训练样本（Zhang等人，2025年），鉴于雪盖的显著时空异质性，确保这些样本的代表性对于模型在大范围和长时间内的有效性至关重要（Luan等人，2022年）。作为一种替代方案，时空融合技术提供了一种有前景的解决方案，旨在生成具有精细空间和时间分辨率的数据（Wang等人，2023年）。它至少需要同一天的一个粗分辨率和一个细分辨率图像对，以及目标时间的粗分辨率图像来预测细分辨率图像（Li等人，2020年）。然而，时空融合本质上涉及估计异质景观和突然变化的不确定性（Ghamisi等人，2019年），这在描述具有快速时空变化的雪盖时会被放大。由于复杂的地形和变化的太阳辐射，山区雪盖分布表现出显著的空间异质性，其特征是薄、稀疏和碎片化的覆盖。此外，雪盖在短时间内会发生快速变化，尤其是在积雪和融化期间（Slater等人，2013年）。在过去的十年中，多项研究应用了各种时空融合方法，结合了来自Landsat、Sentinel-2和MODIS的两种类型的数据，生成了高时空分辨率的雪盖数据（Bousbaa等人，2022年；Liu & Zhang，2022年；Mityók等人，2018年），并缓解了快速时空变化带来的挑战（Gao等人，2022年；H. Wang等人，2022年；Y. Wang等人，2022年；Xiao等人，2024年）。例如，Y. Wang等人（2022年）结合了时空数据融合方法和深度学习，融合了MODIS和Landsat数据，解决了由于雪盖快速变化引起的光谱失真问题，并生成了30米空间分辨率的每日雪盖数据。同样，Gao等人（2022年）将归一化差异雪指数（NDSI）和数字高程模型（DEM）信息整合到增强型空间和时间自适应反射率融合模型（ESTARFM）中，以提高其在雪覆盖山区中的适用性，并减轻由表面突然变化引起的辐射异常。值得注意的是，ESTARFM（Zhu等人，2010年）由于其在该地区的有效性能而被用于雪盖监测。Guo等人（2025年）基于ESTARFM模型，将更高时空分辨率的Sentinel-2数据与MODIS融合，增加了数据密度，更好地捕捉了中国祁连山脉东部的快速雪动态。同时，Dong等人（2024年）比较了两种融合策略：先混合后索引（BI）和先索引后混合（IB），结果表明IB在NDSI方面取得了更好的精度。这一结论也使用ESTARFM方法得到了证实。Xiao等人（2024年）提出了一种复合ESTARFM方法，该方法结合了相邻±2天数据中的时空信息，结合MODIS和被动微波数据，生成了北美地区的1公里分辨率每日FSC估计。扩展时间范围并在时空融合中结合相邻天的雪盖信息，为获得准确且空间完整的雪盖数据提供了机会（Y. Ma等人，2023年）。然而，大多数时空融合研究的共同限制是它们依赖于将单源粗分辨率数据（例如MODIS）与单源细分辨率数据（例如Landsat或Sentinel-2）进行融合（Gevaert & García-Haro，2015年）。云层污染限制了单源粗分辨率数据中有效像素的可用性，而单源细分辨率卫星的不频繁重访频率增加了错过降雪和融雪事件的风险（Stewart，2009年）。此外，大多数现有研究仅依赖一两个无云且空间完整的图像对进行时空融合（Liu等人，2021年）。然而，云层和云阴影经常阻碍高质量图像的获取，特别是在山区，时间序列高分辨率数据中经常观察到长时间的数据间隙（Jia等人，2011年），这进一步限制了它们在大规模和长期雪盖监测中的适用性（Shen等人，2015年；Wang & Atkinson，2018年）。为了解决这些挑战，我们提出了一种自适应时间序列融合框架，将MODIS、VIIRS、Landsat和Sentinel-2数据整合到粗细时间序列中，以提高数据可用性，从而帮助捕捉雪盖的快速时空动态。具体来说，我们引入了时间序列图像对作为基于ESTARFM方法的时空融合框架的输入，该方法不受细分辨率数据空间完整性的限制，并能够为每个像素动态选择最接近目标时间的时间信息。我们使用参考细分辨率数据和地面观测进行了广泛的验证，以评估融合的30米分辨率每日雪盖数据在监测雪盖变化方面的准确性和有效性。此外，我们还定量评估了30米分辨率数据捕捉山区雪盖时间动态和空间分布的能力，探索了雪盖模式的新见解。本研究的主要贡献如下：

我们提出了一种多传感器时空融合方法，整合了多源粗分辨率和多源细分辨率数据，生成30米分辨率的每日雪盖数据。通过有效增加有效像素的数量并提高数据密度，该方法能够准确捕捉动态雪盖变化。我们提出了一种基于ESTARFM的自适应时间序列融合框架，专门设计用于30米分辨率的每日雪盖重建，能够处理稀疏观测和云层覆盖的数据间隙。该方法动态选择并应用与目标时间最接近的观测数据中的时空信息进行每个像素的融合，从而克服了传统上对连续细分辨率输入数据的需求。我们使用参考细分辨率数据和地面观测验证了融合的30米分辨率每日雪盖数据的有效性。结果证明了其在捕捉雪盖快速时间动态方面的优势，并揭示了山区雪盖空间分布的新见解。

**2 研究区域和数据**

**2.1 研究区域**

研究区域是中国新疆的乌鲁木齐（42°54′16″N–44°58′16″N，86°46′10″E–88°59′48″E），面积约为14,193平方公里（图1）。该地区三面被天山山脉环绕，山峰高度可达5,000米。该地区位于中纬度西风带内，具有典型的大陆性干旱至半干旱气候，特征是降水量低、气候干燥以及昼夜温差显著。此外，所有河流均为内流河，主要由融雪补给，在春季灌溉高峰期对下游绿洲至关重要（Wang等人，2011年）。因此，准确监测雪盖分布和动态对于该地区的水资源管理、农业生产及生态可持续性至关重要。图1展示了新疆乌鲁木齐的位置，以及该地区的主要河流、湖泊和气象站。

**2.2 数据**

本研究采用了多种数据集进行多传感器时空融合，包括多源粗分辨率数据和多源细分辨率数据（表1）。所有用于研究的卫星图像数据均来自Google Earth Engine（GEE），并在GEE平台上进行了预处理，覆盖了2022年9月1日至2023年8月31日的水文年度。此外，还引入了辅助数据，包括用于全面分析地形影响的DEM和来自气象站的雪深数据以进行验证。表1。2022年9月至2023年8月期间研究中使用的多传感器数据描述

传感器 数据集 空间分辨率 时间分辨率 图像数量
Terra MODIS MOD09GA.061 500米 1天 365
MOD10A1.061 500米 1天 365
Aqua MODIS MYD09GA.061 500米 1天 365
MYD10A1.061 500米 1天 365
S-NPP VIIRS VNP09GA.001 500米/1000米 1天 365
Landsat-8 OLI LC08/C02/T1_L2 30米 16天 192
Landsat-9 OLI LC09/C02/T1_L2 30米 16天 193
Sentinel-2 A/B MSI S2_SR_HARMONIZED 10米/20米/60米 5天 1,174

2.2.1 多源粗分辨率数据及其预处理
2.2.1.1 MODIS数据
本研究使用的MODIS数据包括每日500米分辨率的MODIS地表反射率产品（MOD/MYD09GA V061）和每日500米分辨率的MODIS积雪覆盖产品（MOD/MYD10A1 V061）。MOD/MYD09GA数据（Vermote & Wolfe, 2021）提供了可见光和短波红外波段的1-7个反射率通道，以及包含云层掩膜的数据。MOD/MYD10A1数据用于积雪覆盖映射算法中获取NDSI数据（Hall et al., 2016）。我们选择了值范围在-1到1之间的原始NDSI层。通过使用MODIS云层掩膜移除了云层和云影像素，以确保数据质量。图像被投影到WGS84椭球坐标系中。为了实现像素级的数据融合，MODIS数据通过立方卷积插值被重采样到30米分辨率，以便与高分辨率数据进行时空融合。此外，自2020年起，全球无云MODIS NDSI数据集（MOD10A1F V061）已在线提供。该数据集是一种新的每日填充云层的积雪覆盖产品，它通过保留前几天的晴空观测数据来填补MODIS Terra NDSI产品中的云层覆盖区域（Hall & Riggs, 2020）。然而，由于中国积雪覆盖具有周期性和短暂性的特点（Jing et al., 2022），该数据集的适用性有限。因此，我们仅使用该数据集与融合结果进行比较。

2.2.1.2 VIIRS数据
VIIRS安装在Suomi国家极轨合作伙伴（S-NPP）卫星上。随着MODIS的老化，VIIRS被认为可以为MODIS提供连续的观测数据（Román et al., 2024）。MODIS和VIIRS仪器的设计相似，都适合用于积雪覆盖的映射（Riggs et al., 2017; Zhang et al., 2020）。本研究选择了每日500米和1公里分辨率的VIIRS地表反射率产品（VNP09GA V001）（Vermote et al., 2023）。同样，对VIIRS数据进行了云层掩膜处理，并从反射率数据中计算了NDSI数据，其他处理步骤与MODIS的数据处理一致。

2.2.2 多源高分辨率数据及其预处理
2.2.2.1 Landsat OLI数据
Landsat-8具有与Terra MODIS相似的轨道参数和过境时间，常被用作时空融合的高分辨率输入数据。Landsat-8操作陆地成像仪（OLI）的地表反射率产品包含4个可见光波段、1个近红外波段和2个短波红外波段，分辨率为30米，重访周期为16天。该产品使用Land Surface Reflectance Code（LaSRC）进行了大气校正。共收集了192张Landsat-8 OLI图像和193张Landsat-9 OLI图像。QA_PIXEL位掩码带用于屏蔽由云层、云影和卷云引起的低质量观测。

2.2.2.2 Sentinel-2 MSI数据
Sentinel-2任务具有宽幅覆盖范围和短的重访周期，单颗卫星的重访周期为10天，两颗卫星组合后的重访周期为5天。Sentinel-2A/B搭载了多光谱仪器（MSI），提供13个地表反射率波段，空间分辨率为10-60米（Jia et al., 2021）。我们收集了1,174张Sentinel-2 MSI图像，并屏蔽了由不透明云层和卷云引起的低质量观测。

2.2.2.3 高分辨率填充间隙的积雪覆盖数据
由Copernicus陆地监测服务（2025）发布的基于Sentinel-1和Sentinel-2的Gap-filled FSC（GFSC）产品整合了最新的观测数据，生成了空间完整的积雪覆盖地图。该产品的空间分辨率为60米，为每个像素提供了积雪覆盖比例（0%-100%），并对受云层覆盖、阴影和无数据区域影响的像素进行了屏蔽（Krebs et al., 2025）。我们将该数据集重采样到30米分辨率，并与我们生成的每日30米积雪覆盖数据进行了比较，以验证所提出方法的有效性。

2.2.3 辅助数据
2.2.3.1 气象站观测数据
收集了2022年9月1日至2023年8月31日期间中国气象局（CMA）的每日气象观测数据，以验证时空融合的性能和融合积雪覆盖数据的准确性。数据遵循特定协议：雪深小于0.5厘米的记录为0厘米，而雪深等于或大于0.5厘米的记录四舍五入到最接近的整厘米。本研究将雪深大于0厘米的观测记录视为积雪覆盖。每日雪深数据来自10个站点，共计3,525条记录（图1）。此外，我们还收集了美国国家海洋和大气管理局（NOAA）下属的国家环境信息中心（NCEI）提供的10个站点的442条雪深观测数据（https://www.ncei.noaa.gov/data/global-summary-of-the-day/archive/），以验证阿尔卑斯山研究区域的积雪覆盖数据（见第4.4节）。这些站点提供雪深记录，但不包括无积雪覆盖的观测数据。

2.2.3.2 数字高程模型（DEM）
研究中使用的是NASADEM（NASA JPL, 2020），用于分析不同地形因素对不同空间分辨率下积雪覆盖的影响。NASADEM是航天飞机雷达地形任务（SRTM）数据的重新处理版本，空间分辨率为30米。它通过整合多个辅助数据集来提高精度并减少空洞。基于NASADEM计算了坡度和朝向，以便后续分析。

3 方法论
3.1 多传感器数据的时空融合
在多传感器数据融合中，我们应用了一种增强的时空融合方法（即ESTARFM），整合了粗分辨率和细分辨率数据的时空信息，以捕捉快速的积雪覆盖变化。该方法需要两对无间隙的粗分辨率和细分辨率图像以及目标时间的1幅粗分辨率图像（Zhu et al., 2010），但由于云层污染（Ju & Roy, 2008），这些数据往往难以获得，进一步限制了其在大规模和长期应用中的适用性。为了解决这个问题，Xiao等人（2024）提出了一种多阶段数据融合框架，通过整合粗分辨率和细分辨率输入来生成无间隙的1公里FSC数据。他们的方法使用了时空邻域立方体来填充被动微波FSC和MODIS FSC数据集中的间隙，然后在5天窗口内进行基于ESTARFM的融合，并进行后处理以填充间隙，确保空间完整性。受启发于在更广泛的时空融合框架中结合ESTARFM和时间窗口的想法，我们的方法适应了更高分辨率数据（例如Landsat和Sentinel-2），生成了30米分辨率的每日积雪覆盖数据。由于细分辨率数据通常具有3-4天的重访周期，并且容易因云层覆盖和卫星过境限制而产生间隙，因此狭窄的时间窗口（例如±2天）不足以识别足够的有效像素。为了解决这个问题，我们扩展了Xiao等人（2024）的概念，实施了自适应时间序列融合策略。在这种方法中，我们引入了多源粗分辨率和细分辨率数据的时间序列对作为输入，从而能够自适应地选择有效的时空信息进行融合。此外，它只要求粗分辨率数据是无间隙的，这确保了数据完整性，因为时空融合模型有严格的输入要求（Goyena et al., 2023），同时允许细分辨率数据中存在空间间隙（支持信息S1中的图S1）。具体来说，该方法分为两个步骤：时间序列数据准备和多传感器数据融合（图2）。在时间序列数据准备阶段，多源粗分辨率数据集（MODIS和VIIRS）和细分辨率数据集（Landsat和Sentinel-2）被协调和整合，以提高有效像素的可用性和数据密度。通过时空插值填充粗分辨率数据中的缺失像素，以确保空间完整性（Abowarda et al., 2021; Chen et al., 2024; Long et al., 2019）。相比之下，细分辨率数据中的空间间隙是可以接受的，因为我们通过时间将细分辨率数据与粗分辨率数据匹配，并将它们作为时间序列纳入融合模型中，从而有效利用不同时间的时空信息。

3.1.1 时间序列数据准备
在数据预处理阶段，为了确保不同数据集之间的一致性，对它们进行了协调处理。鉴于MODIS和VIIRS仪器的设计相似，以及先前的研究已经证明了Landsat和Sentinel-2数据的有效整合（Liu et al., 2020; Mandanici & Bitelli, 2016; Pahlevan et al., 2019; Zhang et al., 2018），我们使用普通最小二乘（OLS）回归对MODIS和VIIRS数据的反射率波段以及OLI和MSI数据进行了标准校准。协调后的结果显示Landsat和Sentinel-2之间，以及MODIS和VIIRS之间有很好的匹配（表2）。OLS回归的R2值均高于0.86，表明回归系数适合于协调粗分辨率和细分辨率传感器之间的地表反射率和NDSI。协调后的Terra、Aqua和VIIRS数据被合并成一个粗分辨率数据集，协调后的Landsat-8 OLI、Landsat-9 OLI和Sentinel-2 MSI数据被合并成一个细分辨率数据集。在合并过程中，如果任何传感器没有提供有效值，则为像素分配一个空值。如果只有一个传感器提供了值，则直接使用该值。如果多个传感器为同一像素报告了值，则计算它们的平均值并分配该值。

为了评估整合MODIS和VIIRS数据集的有效性，我们计算了有效观测像素的增加量，定义为MODIS和VIIRS数据集中的有效像素与MODIS数据集中的有效像素之比。如图3所示，整合显著提高了有效像素的可用性，在雪季相比非雪季有显著改善。同样，多源细分辨率数据集的组合也增加了有效观测的数量。Landsat-8/Landsat-9和Sentinel-2的平均每月有效像素率分别约为10.0%和30.9%，整合后增加到38.3%（图4a）。此外，传感器的互补性质提高了数据密度，这一点通过比较12月的观测结果得到了证明（图4b）。

通过合并MODIS和可见光红外成像辐射计套件的多传感器观测数据，有效观测的比例增加了。图4显示了Landsat和Sentinel-2单独数据集及组合数据集的有效像素比例和数据密度。（a）每月有效像素的平均比例；（b）以12月的数据密度为例。由于云层污染的影响，组合的粗分辨率数据中大约有28.8%的数据间隙。为了获得无间隙的粗分辨率数据，我们使用了时空插值方法来填充这些间隙（Hou et al., 2022; Zhang et al., 2025）。积雪覆盖在短时间内可能会迅速变化，先前的研究表明2-8天的时间窗口可以有效捕捉这些动态变化（Jing et al., 2022）。因此，采用了9天的时间窗口（±4天）来平衡时间代表性和数据可用性。插值基于缺失观测前后时间最接近的有效像素值。线性插值的公式如下：

其中y是缺失时间t的插值像素值，
y1和y2分别是缺失时间之前的和之后的有效像素值，
t分别是缺失时间t1和t2。如果在9天窗口内应用时间插值后仍有缺失像素，则随后进行空间插值来填充间隙。采用了一种逆距离加权方法，使用3×3的空间窗口，其中缺失的像素值是根据其八个相邻像素的值来估计的，权重与距离成反比（Pan等人，2024年；Zhao等人，2022年）。这种时空插值方法生成了无间隙的粗分辨率数据，确保了数据的连续性和完整性。由于在整合多源数据集后仍有大约54.4%的缺失值，因此使用时空插值来填补这些间隙具有挑战性，从而难以获得空间上完整的细分辨率数据。因此，我们通过将细分辨率数据与同一天的粗分辨率数据配对来构建时间序列图像对，用于融合输入，这消除了连续融合对无间隙细分辨率数据的需求，并克服了长期和大规模应用的限制。

3.1.2 多传感器数据融合
在多传感器数据融合阶段，我们采用了ESTARFM时空融合方法（Zhu等人，2010年），该方法通过应用基于粗分辨率影像的时间变化的空间加权方法来估计细分辨率的时间变化。尽管ESTARFM方法最初是为反射率融合开发的，但先前的研究已经证明它也适用于指数融合（Jarihani等人，2014年）。为了生成30米的日雪盖数据，我们应用SNOMAP算法将融合数据转换为二进制雪盖数据，这需要NDSI、绿带和近红外（NIR）波段。因此，我们对NDSI进行了指数融合，对绿带和NIR波段进行了反射率融合。ESTARFM最初需要五张可用图像，包括一张粗分辨率图像（

${C}_{{t}_{p}$）在预测时间（${t}_{p}$），以及在预测时间之前（${t}_{m}$）和之后（${t}_{n}$）获取的两对粗分辨率和细分辨率图像，分别表示为$\left({C}_{{t}_{m}},{F}_{{t}_{m}\right)$和$\left({C}_{{t}_{n}},{F}_{{t}_{n}\right)$。ESTARFM的实现包括四个步骤：（a）在局部搜索窗口内基于相似的光谱特性寻找与中心像素相似的像素；（b）根据光谱和空间特征计算相似像素的权重；（c）计算粗像素和细像素之间的转换系数；（d）估计预测时间中心像素的细分辨率值。主要方程如下：

F$\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{p},B\right)={T}_{m}\times {F}_{m}\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{p},B\right)+{T}_{n}\times {F}_{n}\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{p},B\right)$
F$_{k}\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{p},B\right)=F\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{k},B\right)+\sum\limits _{i=1}^{N}\,{W}_{i}\times {V}_{i}\times \left(C\left({x}_{i},{y}_{i},{t}_{p},B\right)-C\left({x}_{i},{y}_{i},{t}_{k},B\right)$
其中$F$和$C$分别代表细分辨率和粗分辨率值；$F\left({x}_{w/2},{y}_{w/2},{t}_{p},B\right)$是预测时间内搜索窗口大小为$w$的中心像素的细分辨率值的最终预测；$T_{m}$和$T_{n}$分别是时间$t_{m}$和$t_{n}$的时间权重；$F_{k}\left({x}_{w/2},{y}_{w/2},{t}_{p},B\right)$是通过结合所有相似像素的空间权重和基于时间$t_{k}$的粗分辨率值的时间变化计算出的中心像素的细分辨率值。权重$W_{i}$表示每个相似像素对中心像素的影响，转换系数$V_{i}$是使用回归分析为每个相似像素计算的；$N$是相似像素的数量，$(x_{i},{y}_{i})$是相似像素的位置。时间权重$T_{n}$和空间权重$W_{i}$是使用以下方程计算的：

T$_{k}=\frac{\frac{1}{\vert \sum\limits _{j=1}^{w}\,\sum\limits _{l=1}^{w}\,C\left({x}_{j},{y}_{l},{t}_{k},B\right)-\sum\limits _{i=1}^{w}\,\sum\limits _{l=1}^{w}\,C\left({x}_{j},{y}_{l},{t}_{p},B\right)\vert }}{\sum\limits _{k=m,n}(1/\vert \sum\limits _{j=1}^{w}\,\sum\limits _{l=1}^{w}\,C\left({x}_{j},{y}_{l},{t}_{k},B\right)-\sum\limits _{i=1}^{w}\,\sum\limits _{l=1}^{w}\,C\left({x}_{j},{y}_{l},{t}_{p},B\right)\vert },(k=m,n)$
W$_{i}=\frac{\frac{1}{{D}_{i}}}{\sum\limits _{i=1}^{N}\,\left(\frac{1}{{D}_{i}}\right)}$
R$_{i}=\left(1-{R}_{i}\right)\times {d}_{i}$
$R_{i}=\left(1-{R}_{i}\right)\times {d}_{i}$
$F_{i}=\left\{F\left({x}_{i},{y}_{i},{t}_{m},{B}_{1}\right),\mathit{\dots },F\left({x}_{i},{y}_{i},{t}_{m},{B}_{n}\right),F\left({x}_{i},{y}_{i},{t}_{n},{B}_{n}\right)\right\}$
$C_{i}=\left\{C\left({x}_{i},{y}_{i},{t}_{m},{B}_{1}\right),\mathit{\dots },C\left({x}_{i},{y}_{i},{t}_{m},{B}_{n}\right),C\left({x}_{i},{y}_{i},{t}_{n},{B}_{n}\right)\right\}$
$d_{i}=1+\frac{\sqrt{{\left({x}_{\frac{w}{2}}-{x}_{i}\right)}^{2}+{\left({y}_{\frac{w}{2}}-{y}_{i}\right)}^{2}}{\frac{w}{2}}$
$W_{i}$的范围是从0到1，所有相似像素的总权重为1；$D_{i}$结合了相似像素的光谱和地理距离；$R_{i}$是相似像素的细分辨率和粗分辨率之间的光谱相关系数；$F_{i}$和$C_{i}$是包含每个波段在日期$t_{m}$和$t_{n}$的反射率的光谱向量；$d_{i}$是相似像素与中心像素之间的地理距离。相似像素的识别方法如下：

$\begin{array}{c}\vert F\left({x}_{i},{y}_{i},{t}_{k}\right)-F\left({x}_{\frac{w}{2}},{y}_{\frac{w}{2}},{t}_{k}\right)\vert \mathit{\le }\frac{2\sigma }{b}\end{array}$
其中$\sigma$是每个波段中值的标准差，$b$代表土地覆盖类别。我们的多传感器融合方法基于ESTARFM，并通过整合时间序列图像对来扩展时间相邻输入的选择（支持信息S1中的图S1）。它不是依赖于固定的时间对，而是动态地识别每个像素最近的有效细分辨率观测，并使用相应的粗细图像对进行融合。这种自适应策略允许每个像素利用最相关的时空信息，从而获得更可靠和连续的融合预测。具体来说，我们使用了预测时间之前和之后的时间序列图像对（

C$T?i},F_{T?i},i=1,2,\mathit{\dots },m$和之后（C$T+j},F_{T+j},j=1,2,\mathit{\dots },n$）以及预测时间的粗分辨率图像（C$T$），作为增强型ESTARFM的输入（图2）。对于每个像素，该方法动态地识别时间序列中最接近的有效细分辨率观测，并检索相应的图像对（C$T?i},F_{T?i},C$T+j},F_{T+j}$），然后将这些特定于像素的图像对与C$T$整合到ESTARFM框架中，以逐像素为基础估计预测时间的细分辨率值。按照这个程序，我们在研究区域生成了一年的30米日数据，包括NDSI、绿带和NIR反射率波段。

3.2 融合结果的评估
我们使用了参考细分辨率数据，并采用了包括相关系数（R）、平均绝对误差（MAE）和均方根误差（RMSE）等指标来评估融合数据（J. Ma等人，2022年）。此外，我们还使用了R、相对平均绝对误差（rMAE）和相对均方根误差（rRMSE）来评估和比较指数融合和反射率融合之间的差异。进一步地，我们使用气象站的雪深测量数据来评估融合数据在雪盖映射方面的性能。具体来说，我们应用SNOMAP算法对融合数据进行处理，生成了30米的日二进制雪盖数据（即雪/无雪），如果NDSI≥0.4，则像素被视为有雪，绿带的反射率≥0.1，近红外波段的反射率>0.1（Hall等人，1995年）。雪盖映射的性能使用以下指标进行评估：总体准确性（OA）、Kappa系数、精确度、召回率和F1分数，值越接近1表示结果越好。为了进一步评估生成的30米日雪盖数据捕捉雪盖时间动态的能力（Xiao等人，2024年；H. Wang等人，2022年），计算和验证了雪盖指标，包括雪盖天数（SCD）、雪盖开始日期（SCOD）和雪盖结束日期（SCED）。雪盖天数是指在一个水文年内记录的雪覆盖天数（2022年9月1日至2023年8月31日）。SCD的计算方法如下：

SCD=$\sum\limits _{i=1}^{n}\left({S}_{i}\right)$
其中$S_{i}$表示二进制雪盖像素值，1表示有雪，0表示无雪，n表示水文年内的像素数量。SCOD和SCED分别指水文年内连续5天雪覆盖期的第一天和最后一天（Chen等人，2015年）。SCOD和SCED的计算方法如下：

SCOD
=
min
{
t
}
其中
t ∈ [0, m-4]

并且
∑k = 0, k = 0, …, 4
S_t + k = 5

(13)

SCED
=
max
{
t
}
其中
t ∈ [0, m-4]

并且
∑k = 0, k = 0, …, 4
S_t + k = 5

(14)

其中
t 表示一个日期，
m 表示一个水文年中的日期数量，
S_t 表示一个像素的值，1表示有雪，0表示无雪。

4. 融合结果的准确性评估

为了评估多传感器时空融合方法在生成融合数据方面的性能，我们根据雪季的一般特征将雪季分为三个时期：积雪期（9月至11月）、雪稳定期（12月至2月）和融雪期（3月至5月）（Ke等人，2016年；Sun等人，2024年）。我们从这些时期选择了九个场景，分别代表积雪期（2022年11月14日、25日、26日）、稳定期（2022年12月20日、24日、28日）和融雪期（2023年3月6日、16日、21日）。此外，我们还从每个月选择了晴朗的天空场景进行验证，并计算了月平均值。NDSI的融合性能在不同雪季有所不同（图5和表3）。值得注意的是，雪稳定期的准确性最高，其特点是相关性更强（R = 0.926–0.964）和误差更低（MAE = 0.102–0.150；RMSE = 0.151–0.199）。这种更好的性能主要归因于在低温条件下雪盖的持续性，这确保了有利于数据融合的时空一致性。相比之下，在积雪期和融雪期，频繁且快速的降雪和融雪事件导致雪盖发生高度动态和不稳定的变化，从而显著降低了融合准确性。此外，指数融合与反射率融合的比较显示，指数融合通常表现出更高的相关性和更低的rMAE和rRMSE（图6），这与先前的研究结果一致（Dong等人，2024年）。总体而言，评估表明所提出的数据融合框架在生成30米分辨率的每日雪盖数据方面是有效的。

图5展示了30米参考归一化差异雪指数（NDSI）与融合NDSI之间的关系，分别对应积雪期（a–c）、稳定期（d–f）和融雪期（g–i）。表3显示了2022年9月至2023年8月期间融合归一化差异雪指数的月平均准确性验证结果。

表3. 2022年9月至2023年8月期间融合归一化差异雪指数的月平均准确性

图6展示了30米融合归一化差异雪指数、绿色NIR反射率波段的准确性，包括相对平均绝对误差（R）（a）、相对均方根误差（b）和rRMSE（c）。

4.2 多源数据的时空融合有效性

我们使用不同的数据源对融合结果进行了准确性评估（图8）。多传感器数据融合的结果优于仅使用Landsat或MODIS的数据融合结果，这突显了结合多传感器数据源的有效性。特别是，多传感器高空间分辨率数据的结合显著提高了准确性，因为精细分辨率图像的时间密度增加提供了更多接近目标时间的有效像素。

4.3 基于地面观测的二元雪盖映射的准确性比较

我们使用SNOMAP算法生成了30米和500米的融合每日二元雪盖地图，以评估空间分辨率对雪盖映射的影响。我们使用2022年9月至2023年5月期间从10个气象站收集的雪深观测数据对30米和500米的每日二元雪盖进行了验证。如表4所示，30米数据在准确性（OA）、Kappa值、召回率和F1分数方面优于500米数据，分别提高了0.07、0.13、0.14和0.09。这种提高主要是由于30米空间分辨率更高，能够捕捉更细微的细节，从而更准确地识别雪覆盖区域。这在地形复杂的山区尤为重要，因为局部变化（如小雪斑或融化模式）可能显著影响雪盖估计。图11显示，30米数据能够有效捕捉精细的地形特征和小规模的雪覆盖差异。

4.3.3 不同地形因素下山区雪盖的比较

我们从30米和500米的每日雪盖数据中提取了三个雪相指数（SCD、SCOD和SCED），以比较雪盖的时间动态。使用来自8个气象站的2,920个连续观测数据（这些数据具有完整的年度雪深记录和超过30天的SCD），我们评估了这两组数据（图12）。30米数据的SCD、SCOD和SCED与气象站结果非常吻合，平均差异分别为9天、5天和6天。相比之下，500米数据的差异几乎翻倍，分别为17天、11天和13天。这表明30米每日雪盖数据更准确、更一致地捕捉了雪盖的时间变化，与地面观测结果更为一致。

4.3.2 提取的雪相指数差异

我们从30米和500米的每日雪盖数据中提取了三个雪相指数（SCD、SCOD和SCED），以比较不同地形下雪盖的时间动态。我们选择了天山山脉的一个区域作为代表区域，分析了不同坡度、海拔和朝向下的雪覆盖像素计数。如果 ${\mathrm{D}}_{i} > 0$，即第 $i$ 个区域的30米分辨率雪覆盖像素数量超过500米分辨率的数量；如果 ${\mathrm{D}}_{i} < 0$，则500米分辨率的数量更多。当 ${\mathrm{D}}_{i} = 0$ 时，30米和500米分辨率的雪覆盖像素数量相等。我们将海拔高度划分为7个500米间隔的区间，并计算每个海拔区间内30米和500米分辨率雪覆盖数据的雪覆盖像素数量。在3500米以下的区域，30米分辨率数据的雪覆盖像素数量始终高于500米分辨率数据（见图13）。然而，在3500米以上的高度，500米分辨率数据显示出更高的雪覆盖程度，这可能是由于在崎岖的高海拔地区存在混合像素导致的过度估计。30米更精细的分辨率使得雪覆盖的检测更加精确，尤其是在低海拔和中等海拔地区，从而提高了雪覆盖监测的准确性。由于混合像素效应，500米分辨率数据经常显示连续的雪斑，无法准确捕捉复杂地形中的雪分布，特别是在中高海拔地区（见图14c）。在较低海拔地区，许多小的雪覆盖区域也会被遗漏（见图14d）。相比之下，30米分辨率数据能更好地捕捉复杂的雪模式，准确地描绘出复杂地形中的雪覆盖情况，并保留低海拔地区小雪斑的详细记录（见图14a和14b）。

图13展示了整个雪季（a）不同海拔范围内30米和500米日雪覆盖数据平均雪覆盖像素数量的比较，包括雪积累（b）、稳定期（c）和融化期（d）。

图14展示了山区不同海拔范围内30米和500米雪覆盖数据的比较，针对两个代表性场景。我们将坡度划分为7个5°间隔的区间，并计算每个坡度范围内30米和500米雪覆盖数据的雪覆盖像素数量。如图15所示，在缓坡地区（≤25°），30米分辨率数据记录的雪覆盖程度始终高于500米分辨率数据。然而，随着坡度的增加，这种差异减小，当坡度超过25°时，500米分辨率数据显示的雪覆盖程度超过了30米分辨率数据。在陡坡上，雪往往会因风或重力作用而重新分布，从而减少实际雪积累。30米分辨率更高的空间分辨率使得在陡峭地形上更准确地检测到雪覆盖情况，避免了这种过度估计。此外，它还能保留缓坡上零散的雪模式，在融化期间提供更清晰的残余雪分布图。

图15展示了整个雪季（a）不同坡度范围内30米和500米日雪覆盖数据平均雪覆盖像素数量的比较，包括雪积累（b）、稳定期（c）和融化期（d）。我们将朝向划分为8个45°间隔的区间，并计算每个朝向的30米和500米雪覆盖数据的雪覆盖像素数量。在雪积累和稳定期间，30米分辨率数据在所有朝向上都一致地捕捉到更广泛的雪覆盖范围，尤其是在阴凉坡面上差异更为明显。在融化期间，30米分辨率数据显示阴凉坡面上的雪量显著更多，而500米分辨率数据倾向于高估阳光照射坡面上的雪量。由于阴凉坡面因温度较低和融化速度较慢而保持更多雪，30米更精细的分辨率更准确地反映了这种地形驱动的变化。

图16展示了整个雪季（a）不同朝向范围内30米和500米日雪覆盖数据平均雪覆盖像素数量的比较，包括雪积累（b）、稳定期（c）和融化期（d）。

4.4 与高分辨率填补间隙的雪覆盖产品的比较

为了评估所提出的多传感器时空融合框架的有效性，我们将我们的结果与高分辨率雪覆盖产品GFSC进行了比较。由于该数据集覆盖了欧洲，我们选择了阿尔卑斯山区的一个区域作为研究区域（见支持信息S1中的图S2），并生成了2022年9月至2023年5月的30米日雪覆盖数据以进行比较。我们首先使用0.5的阈值将GFSC转换为二值图像，然后使用NCEI提供的气象站雪深观测数据对二值化的GFSC和融合后的30米雪覆盖数据进行了验证。由于雪深观测数据仅包含正样本，因此只能计算遗漏误差（OE）作为评估指标。GFSC产品的遗漏误差（0.11）略低于我们的融合30米雪覆盖数据（0.13）。这可以归因于我们的雪覆盖数据仅依赖于光学观测，而GFSC还结合了合成孔径雷达（SAR）数据，为复杂地形中的可靠检测提供了更多信息。总体而言，尽管存在这种小差异，两种产品在检测雪覆盖方面表现出可比的能力。尽管GFSC产品提供了雪覆盖的空间完整概览，但它仍然包含标记为云或无数据的像素。相比之下，我们提出的方法生成了空间覆盖完整的无间隙雪覆盖数据。为了评估我们的融合30米雪覆盖数据是否能在GFSC缺失的区域提供可靠的估计，我们使用标记为云或无数据的像素对应的雪深观测数据进行了验证。验证结果表明，我们的融合数据具有良好的准确性，遗漏误差为0.15，证明了其在GFSC数据缺失区域可靠估计雪覆盖的能力。图17显示了2022年9月至2023年5月阿尔卑斯山区研究区域内GFSC产品和融合30米雪覆盖数据的日数据完整性百分比。数据完整性百分比是基于研究区域内每天的所有有效像素计算的。GFSC产品的平均数据完整性为50.4%，而融合30米雪覆盖数据实现了完整的空间覆盖，从而得到了无间隙的雪覆盖数据集。此外，如图18所示，我们的融合30米雪覆盖数据与GFSC产品展示了高度一致的空间模式，同时有效填补了其数据间隙。

5 讨论

5.1 多传感器数据融合的优势和局限性

我们整合了多传感器数据，并采用了时间序列时空融合方法，使用粗分辨率和细分辨率图像的对来更好地捕捉快速的雪覆盖动态。MODIS和VIIRS数据的整合增加了有效雪季观测的数量。尽管使用了间隙填补技术来生成无间隙的粗分辨率输入，但大量缺失数据的存在可能会引入插值伪影，从而可能降低融合输出的真实性。在未来的工作中，我们将进一步研究和优化粗分辨率雪覆盖产品的间隙填补策略，并评估它们对融合性能的影响，以提高长期、大尺度高时空分辨率雪覆盖绘制的准确性和鲁棒性。为了解决细空间分辨率下的时间稀疏性问题，我们的研究整合了来自Landsat-8 OLI、Landsat-9 OLI和Sentinel-2 MSI的多源光学观测数据，以及粗分辨率数据，构建连续的时间序列作为融合框架的输入。与Xiao等人（2024年）专注于通过填充MODIS数据中的缺失值来生成无间隙的1公里日FSC不同，我们的研究目标是重建30米日雪覆盖，其中高分辨率传感器的低重访频率和频繁的云污染带来了独特的挑战。我们引入了一种灵活的自适应时间序列融合策略，动态地为每个像素选择最相关的时空信息，而不是像Xiao等人（2024年）那样使用固定的±2天时间窗口。这种方法通过适应稀疏的高分辨率观测和多源数据，扩展了传统的融合框架，而不是依赖于预定的图像对。因此，我们的方法不仅增加了时间密度，解决了空间不完整性问题，还生成了无间隙的30米日雪覆盖地图。与基于单源传感器数据的融合结果相比，基于多源整合的融合输出在捕捉快速雪覆盖动态方面表现出更大的灵活性、有效性和准确性。此外，多传感器数据融合的一个关键挑战是解决不同数据集之间的一致性问题，特别是在雪覆盖具有快速时空变化的情况下。我们使用OLS方法对不同数据集进行了对齐，帮助它们达到一个共同的尺度。然而，考虑到传感器之间的空间分辨率、太阳角度和信噪比差异（Wang & Hu, 2018），需要更全面的方法来解决它们之间的不一致性。因此，未来的工作应专注于开发更稳健的方法来减轻潜在的不确定性，或者考虑使用如Harmonized Landsat Sentinel-2（HLS）数据集（Ju等人，2025年）这样的协调数据集，以更好地解决不同数据源之间的不一致性。除了这些挑战之外，所提出的融合策略的像素级特性还引入了空间一致性的挑战。具体来说，因为不同的像素可能利用来自不同图像对的空间时间信息，所以在快速雪积累和融化期间可能会出现马赛克伪影。这些相邻像素之间的时间不一致性可能导致融合图像中的可见不连续性。未来的研究应优先考虑在充分利用多传感器数据集的时间丰富性的同时保持空间一致性的方法，以进一步推进大尺度、高精度的雪覆盖动态监测。

5.2 30米日数据在山区雪覆盖绘图中的应用潜力

雪覆盖对地形变化非常敏感，因为雪具有平滑效果（Mott等人，2010年）。虽然多项研究表明MODIS雪产品可以有效应对山区的相关挑战（Bormann等人，2018年；Hammond等人，2018年），但雪覆盖和物候趋势表现出显著的区域和海拔变化，以及独特的局部模式（Notarnicola，2020年）。MODIS的粗空间分辨率加上复杂地形限制了雪覆盖估计的准确性。因此，高时空分辨率数据对于更准确地理解山区的雪覆盖时空模式至关重要。尽管有几项研究致力于改进融合模型和提高准确性，但很少有研究使用融合的30米日雪覆盖数据来研究山区雪覆盖的具体模式和动态。在这项研究中，我们对30米和500米日雪覆盖数据进行了比较分析。与地面数据的验证表明，30米日数据表现出更高的准确性，更有效地捕捉了日常雪动态，并与气象站观测到的雪物候更为吻合。在3500米以下或坡度小于25度的地区，30米数据检测到更多的雪覆盖，包括更小、更分散和零散的雪斑。相比之下，500米数据在高海拔和陡坡上倾向于高估雪覆盖，主要是由于混合像素效应。此外，30米数据在所有坡度方向上都提供了更准确的雪覆盖估计，特别是在阴凉坡面上。这一进展对水文建模、气候研究和高山生态系统监测具有重大价值，这些领域都依赖于准确的雪估计。然而，仍然存在挑战。在持续云覆盖或观测输入有限的区域（例如本研究中的条件3），融合数据的可靠性仍需进一步验证。此外，尽管30米数据捕捉到更详细的空间模式，但它们在多样化气候区和雪覆盖制度下的有效性仍有待系统探索。鉴于ESTARFM等时空融合方法的严格输入数据要求和长时间处理时间，我们最初生成了一年的数据（2022年9月至2023年8月）来验证所提出的基于ESTARFM的融合框架是否适合生成大尺度、高时空分辨率的雪覆盖产品。为了考虑气候条件的影响，我们使用标准化降水蒸发指数（SPEI）评估了干旱状况，平均SPEI值为-1.89，表明这是一个干旱年份。因此，基于单一年份的分析可能无法充分捕捉雪的变化，这些变化可能受到年际气候差异的强烈影响（例如，干旱或湿润条件）。尽管存在这一局限性，为期1年的分析仍然证明了所提出的基于ESTARFM框架的有效性，并为未来生成长期、大规模雪盖数据集的工作奠定了基础。此外，未来使用此类高分辨率产品来探索长期趋势的研究还应考虑降水、温度等气候因素，以更准确地理解由气候驱动的雪相变化。

6 结论

在这项研究中，我们提出了一个自适应时间序列融合框架，该框架整合了多传感器数据以捕捉快速的雪盖变化。通过利用粗分辨率和细分辨率数据的时间序列对，该方法有效解决了细分辨率数据在空间完整性方面的局限性所带来的挑战。广泛的验证表明，多传感器融合方法具有很高的相关性（R = 0.776–0.964）和较低的成本误差（MAE = 0.079–0.207，RMSE = 0.151–0.285）。与广泛使用的500米分辨率数据相比，30米分辨率的日雪盖数据在所有指标上都显示出显著改进。它不仅更可靠地捕捉到了雪的变化动态，还为山区雪盖提供了新的见解，特别是在海拔3500米以下、坡度小于25度的地区以及阴坡上，能够检测到更详细的雪模式。未来的研究可以探索融合数据在时间和空间上的一致性，特别是在雪变化快速的区域。此外，使用高时空分辨率数据对不同气候区进行长期趋势分析对于准确监测雪盖至关重要。

致谢

本研究得到了中国国家自然科学基金（项目编号：42471410；42571382）和南京大学人工智能与科学应用项目（项目编号：020914380141）的支持。

利益冲突

作者声明与本研究无关的任何利益冲突。

数据来源声明

我们的研究主要使用了两种类型的数据：遥感数据和气象站观测数据。研究中使用的遥感数据来源于Google Earth Engine (GEE)平台，包括MODIS数据（https://developers.google.com/earth-engine/datasets/catalog/modis）、VIIRS数据（https://developers.google.com/earth-engine/datasets/catalog/NASA_VIIRS_002_VNP09GA）、Landsat OLI数据（https://developers.google.com/earth-engine/datasets/catalog/landsat）、Sentinel-2 MSI数据（https://developers.google.com/earth-engine/datasets/catalog/sentinel）以及DEM数据（https://developers.google.com/earth-engine/datasets/catalog/NASA_NASADEM_HGT_001）（NASA JPL，2020年）。高分辨率的填充间隙FSC数据来源于Copernicus土地监测服务（2025年）。雪深数据的日观测数据来自中国气象局（CMA）（https://data.cma.cn）。此外，阿尔卑斯地区的雪深日观测数据由美国国家海洋和大气管理局（NOAA）全球地表摘要日（GSOD）项目下的国家环境信息中心（NCEI）提供（https://www.ncei.noaa.gov/data/global-summary-of-the-day/archive/）。