尼泊尔卡利甘达基河流域历史及未来永久冻土分布的机器学习评估

《Research in Cold and Arid Regions》：Machine Learning Assessment of Historical and Future Permafrost Distribution in the Kali Gandaki River Basin, Nepal

　　Earina Sthapit | Binod Pokharel | Dibas Shrestha | Deepak Aryal | Arnab Singh | Anita Tuitui
尼泊尔特里布万大学中央水文与气象系，Kirtipur 44618

**摘要**
永久冻土是指至少连续两年保持温度在0°C以下的土壤，对山区生态系统的稳定性和水分平衡至关重要。其退化会增加滑坡和其他冰冻圈灾害的风险。本研究利用再分析数据、岩冰川目录和机器学习算法，为喜马拉雅中部卡利甘达基河流域（KGRB）绘制了高分辨率（30米）的永久冻土概率图。我们使用了三种模型——梯度提升（GB）、随机森林（RF）和支持向量机（SVM），并结合WorldClim、ERA5和CMIP6数据集来分析永久冻土随时间和空间的变化。根据数据，WorldClim显示1961-1970年间永久冻土面积增加了6.60%（119平方公里），而ERA5则报告在同一时期减少了26.65%（724平方公里）。未来的CMIP6预测显示，到2091-2010年，在低排放情景（SSP1-2.6）下永久冻土面积可能增加0.47%（9平方公里），在中等排放情景（SSP2-4.5）下减少9.29%（179平方公里），在高排放情景（SSP5-8.5）下减少40.94%（822平方公里）。持续存在的永久冻土非常脆弱，研究结果表明加强气候缓解措施对于维持其稳定性和减少碳循环至关重要。1961-2020年间，年平均气温每年升高0.011–0.019°C，预计到2010年这一增幅将达到0.053°C。所有数据都表明，永久冻土逐渐向更高海拔（5000–6000米）集中，这反映了相同的气候避难区的形成。

**1. 引言**
永久冻土是指至少连续两年保持温度在0°C或以下的土壤（Brown, 1995），它覆盖了地球陆地表面的大约四分之一（Anisimov和Nelson, 1996；Dobinski, 2011）。尽管其覆盖范围广泛，但它对气候调节、生态系统动态和人类活动仍起着关键作用（Nelson等人，2002；Jorgenson等人，2006；Walter等人，2006）。永久冻土的融化会导致温室气体（Anisimov和Nelson, 1996）和有机碳（Hugelius等人，2014）的释放，同时还会引发土壤变化（Wang等人，2006）、地面沉降（Nelson等人，2002）以及水文变化（Frey和McCllland, 2009；Karlsson等人，2015）。此外，永久冻土的融化还对基础设施构成威胁，并增加地质灾害的风险（Hassan等人，2021）。

1998年，国际永久冻土协会（IPA）发布的《北极永久冻土和地面冰状况地图》首次全面展示了北半球的永久冻土分布情况（Brown和Haggerty, 1998）。后续研究进一步细化了我们对永久冻土分布的理解。例如，Gruber（2012）提供了全球永久冻土区的高分辨率估计，Schmid等人（2015）利用岩冰川作为指标，首次预测了兴都库什-喜马拉雅（HKH）地区的永久冻土情况。值得注意的是，喜马拉雅地区的冰川退缩速度比其他地区更快（Anthwal等人，2006），而这些地区的岩冰川存在是永久冻土存在的重要标志（Baral等人，2020）。多种方法，包括统计和经验模型以及卫星数据，已被用于绘制永久冻土分布图。例如，Landsat 8提供的地表温度数据被用于印度帕尔瓦蒂谷的永久冻土测绘（Pradhan和Shukla, 2022）。机器学习（ML）算法在永久冻土测绘中越来越受欢迎，Deluigi等人（2017）在阿尔卑斯山区进行了研究，Haq和Baral（2019）在印度山区也进行了类似研究，他们估计该区域60%的剖面存在永久冻土。类似地，一项针对青藏高原东北部的研究结合了岩冰川目录和实地观测与ML模型，准确估计了永久冻土的分布（Hu等人，2024）。

永久冻土研究在欧洲阿尔卑斯山区（Cremonese等人，2011）和亚洲阿尔卑斯山区（Gruber, 2012；Schmid等人，2015）变得越来越重要。在尼泊尔的昆布地区，永久冻土出现在4,900–5,000米的高度（Fujii和Higuchi, 1976），到2004年其下边界下降到5,400–5,000米（Fukui等人，2007）。坎钦琼嘉地区的调查显示，由于气候变化，2009至2039年间永久冻土的下限可能上升188米（Mayer等人，2012）。研究表明尼泊尔东部的永久冻土比西部更多（Maharjan, 2018），并且整个尼泊尔的永久冻土下边界每年因气候变化而下降1.3–2.3米（Chauhan和Thakuri, 2017）。

喜马拉雅山的永久冻土不仅受温度影响，还受地震和构造活动（Jakob, 1992）以及坡度和降水模式（Fort, 2000；Bell等人，2021）等因素的影响。全球约8%–10%的永久冻土位于尼泊尔喜马拉雅地区（Gruber, 2012）。尽管其重要性显而易见，但运用ML模型评估该地区永久冻土分布的研究相对较少。准确绘制永久冻土地图对于了解其范围、特征及其因气候变化而可能发生的变化至关重要，同时也利于预测与其退化相关的地质灾害和生态影响。

本研究采用ML模型生成了KGRB的十年期永久冻土分布图（图1）。参考Johns（2015）的分类方法，岩冰川被分为完整型和遗迹型：完整型岩冰川表示存在永久冻土，遗迹型岩冰川表示永久冻土的缺失。研究将岩冰川分类作为响应变量，地形和气候变量作为预测因子来建模永久冻土分布。本研究采用了Baral等人（2020）的方法，其中年平均气温（MAAT）与岩冰川数据结合，作为ML模型中的关键预测变量。选择这一方法的原因有：（1）确保与通常依赖温度指数作为首要控制因素的大规模永久冻土模型的方法一致性；（2）由于整个研究区域内缺乏其他关键变量（如积雪深度、土壤性质和地表特征）的一致、高分辨率空间数据；（3）建立一个清晰、可复制的基线模型。我们明确承认这种简化忽略了其他重要因素，这是研究的局限性之一。虽然岩冰川是研究区域内确认的永久冻土特征，但可能还存在其他形式的永久冻土，包括地下含冰沉积物、基岩永久冻土以及仅通过遥感或地表观测难以识别的不太明显的表现形式。

**图1.** (a) 尼泊尔地图，显示了KGRB的流域及其出口Devghat；(b) KGRB地图，展示了海拔、山峰、完整岩冰川（黄色点）和遗迹岩冰川（粉红色+符号）的空间分布；(c) 著名山峰Dhaulagiri、Tukuche、Nilgiri和Annapurna的位置图，以及KGRB上游的雪地和冰川补给区。

**2. 研究区域**
尼泊尔约有6,000条岩冰川（RGs），分布在海拔3,325至5,676米的范围内，构成了该地区冰冻圈的重要组成部分（Harrison等人，2024）。这些岩冰川估计储存了160亿至250亿立方米的水，对下游生态系统和社区的水资源至关重要（Jones等人，2018）。选择KGRB作为岩冰川测绘和永久冻土分布分析的研究区域。岩冰川的测绘遵循Jones等人（2018）建立的方法框架。地理上，该流域从南部平原延伸至中部喜马拉雅地区，海拔梯度从64米到8,141米，总面积达11,851平方公里（图1）。卡利甘达基河发源于 Mustang区的喜马拉雅高地，向南流，最终在Chitwan的Devghat与Trisuli河汇合。气候上，该流域具有极端温度变化，高海拔地区的冬季温度可降至-25°C，而低海拔地区的夏季温度可升至35°C（Mishra等人，2014）。图1详细展示了流域的海拔梯度、冰川范围和岩冰川分布，为分析永久冻土动态和冰冻圈过程提供了基础。

该流域的冰川和永久冻土带受到显著山峰的影响，包括Tukuche（6,920米）、Dhaulagiri（8,167米）、Nilgiri（7,061米）和Annapurna（8,091米）。值得注意的是，卡利甘达基谷拥有世界上最深的峡谷，Dhaulagiri与Annapurna之间的垂直落差达6,000–7,000米。研究表明，随着气候变化，2009至2039年间永久冻土的下限可能上升188米（Mayer等人，2012）。尼泊尔东部地区的永久冻土比西部更多（Maharjan, 2018），并且整个尼泊尔的永久冻土下边界每年因气候变化而下降1.3–2.3米（Chauhan和Thakuri, 2017）。

喜马拉雅山的永久冻土不仅受温度影响，还受地震和构造活动（Jakob, 1992）的影响，以及坡度和降水模式等因素的影响（Fort, 2000；Bell等人，2021）。尼泊尔喜马拉雅地区的永久冻土约占全球永久冻土的8%–10%（Gruber, 2012）。尽管如此，应用ML模型评估该地区永久冻土分布的研究较少。准确绘制永久冻土地图对于理解其范围、特征及潜在变化至关重要，同时也有助于预测与其退化相关的地质灾害和生态影响。

**3. 数据**
3.1. 平均年平均气温（MAAT）
本研究使用WorldClim版本2（Fick和Hijmans, 2017）提供的1960至2020年的月最高和最低温度数据，水平分辨率为2.5分钟（赤道附近约为21平方公里）。WorldClim版本2数据适用于锡金和北阿坎德喜马拉雅地区的永久冻土研究（Haq和Baral, 2019；Baral等人，2020）。同样，也使用了ERA5-Land（Buontempo等人，2022）提供的月平均温度数据，其水平分辨率为9公里。ERA5在统计评估中表现良好，但在估计永久冻土面积方面存在偏差，尤其是在加拿大和阿拉斯加（Cao等人，2020）。

起初，我们使用了四个代表性的CMIP 6模型来评估不同气候条件下的永久冻土分布：(a) 温暖干燥型：NESM 3_r1i1p1f1_gn（Cao等人，2018）；(b) 温暖湿润型：UKESM 1-0-LL_r1i1p1f2_gn（Sellar等人，2019）；(c) 寒冷湿润型：EC-Earth 3_r1i1p1f1_gr（D?scher等人，2022）；(d) 寒冷干燥型：MPI-ESM1-2-LR_r1i1p1f1_gn（Mauritsen等人，2019）。然而，为了进行长期多情景分析（SSP1-2.6, SSP2-4.5, SSP5-8.5），我们采用了MIROC6作为主要模型，因为它在冰冻圈过程模拟方面具有更强的能力。MIROC6采用了复杂的多层土壤方案，深度达14米，能详细模拟土壤热传导性、冻融过程和地面冰含量——这些在其他模型中被简化或忽略（Tatebe等人，2019；Hajima等人，2020）。比较验证研究表明，MIROC6在再现观测到的永久冻土范围、活动层厚度和土壤温度剖面方面的能力更强（NSE >0.75），优于NESM3（NSE=0.58）、UKESM 1-0-LL（NSE=0.62）、EC-Earth 3（NSE=0.64）和MPI-ESM 1-2-LR（NSE=0.61）（Burke等人，2020；Yokohata等人，2020）。尽管这些模型在捕捉温度-降水不确定性方面有效，但它们缺乏模拟永久冻土融化阈值和退化速率所需的垂直土壤分辨率和热力学复杂性（Peng等人，2020）。此外，MIROC6在永久冻土-碳循环研究中也得到了广泛应用，能够真实反映高纬度地区的雪绝缘效应、冬季土壤热制度和季节性冻融周期，从而控制永久冻土的稳定性（Guo和Wang, 2016；McGuire等人，2018）。该模型具有适度的平衡气候变化敏感性（2.6°C），并能准确模拟不同升温情景（从SSP 1-2.6（约1.5–2°C）到SSP 5-8.5（约4–5°C）下的永久冻土响应（Zelinka等人，2020）。处理四个模型的多个集合成员和三种SSP情景下的温度变量（最高和最低温度）需要超出本研究范围的计算资源，包括大量数据存储、空间降尺度处理、偏差校正和永久冻土模型耦合。因此，采用MIROC6不仅代表了在过程表示方面的科学改进，而且在实际应用中更加注重选择适合特定目的的模型，而不是关注模型集的大小——这种方法在冰冻圈影响评估中越来越受到推崇（Biskaborn等人，2019年；Andresen等人，2020年）。基于Burke等人（2020年）用于估算全球永久冻土范围的方法和Hu等人（2024年）用于评估北半球永久冻土温度条件的技术，CMIP6情景气候模型被用来预测未来的永久冻土分布。

3.2. 岩石冰川分类和不确定性评估
岩石冰川的提取遵循Jones等人（2018年）的方法，采用手动数字化和半自动分类光学卫星图像的技术。工作流程包括：（1）观察特征性的地貌特征，包括具有脊状和沟槽图案的叶状和舌状结构；（2）分析表面纹理和碎屑分布；（3）评估地形位置和高程限制。提取的不确定性使用Schmid等人（2015年）提出的不确定性指数进行量化，该指数评估了边界划分、积雪覆盖的影响以及流动结构的识别。根据活动状态，岩石冰川被分为两类：完整岩石冰川（IRGs）；显示出活跃的运动迹象，包括陡峭的前缘坡度、明显的脊线和清晰的流动特征，表明当前存在永久冻土以及正在进行的蠕动过程；遗迹岩石冰川（RRGs）；表现出退化的、圆形的地貌特征，流动特征不明显或不存在，前缘坡度较缓，暗示曾经的永久冻土条件可能已不再存在。在KGRB栅格块中，共识别出408个岩石冰川，其中362个被分类为完整岩石冰川，46个为遗迹岩石冰川。这表明该地区存在大量永久冻土的潜力。

与岩石冰川清单相关的不确定性通过多方面的框架进行了严格量化。每个绘制出的特征都分配了一个不确定性指数，该指数根据诊断性地貌特征的清晰度进行评分（中等≤5分，高6-10分，非常高≥11分），结果有93%的清单被归类为高或非常高的确定性。与独立的全球永久冻土分区指数进行空间验证时显示出强一致性，87%的绘制出的岩石冰川位于模型预测的永久冻土区域内，只有约4%位于其边界之外。最后，通过纳入冰层厚度变化和冰含量的不确定性（在40%-60%之间），对推导出的水文属性的不确定性进行了统计限定。这种放大后的结果是全国水体积为3.81 ± 0.84立方千米，为后续的水文和气候应用建立了稳健、可量化的置信区间。

3.3. 高程网格
我们使用了一个空间分辨率为30米的数字高程模型（DEM），该模型来自Shuttle Radar Topography Mission（SRTM）30米数据集（https://dwtkns.com/srtm30m/）。这个高程网格被用来利用焦点线性回归方法将再分析和气候模型的输出降尺度到30米分辨率的MAAT网格。降尺度后的MAAT是通过将高程差栅格与重采样的露点率层相乘，然后将其结果添加到基础MAAT层中计算得出的（Barel等人，2020年）。最终的30米高程网格也被用作通过机器学习模型生成永久冻土概率分布图的输入变量。

4. 方法论
4.1. 温度降尺度
从粗分辨率气候数据集（WorldClim、ERA5和CMIP6）将MAAT降尺度到细分辨率（30米）是使用基于环境露点率原理的方法进行的。控制这种降尺度方法的基本方程式如方程（1）所示：
(1)
这里表示细分辨率（30米）下的降尺度温度，表示粗分辨率下的温度（°C），表示环境露点率（°C/m），表示细分辨率下的高程（m），表示粗分辨率下的高程（m）。这个公式基于大气科学的基本原理，即在对恒定露点率（?0.65 °C/100 m）的情况下，温度随着高度的增加而降低（Hay等人，2000年；Tabor和Williams，2010年）。
环境露点率γ是使用焦点线性回归方法在3×3像素窗口内计算的，该方法遵循地形气候降尺度中的既定技术（Daly等人，1994年；Hijmans等人，2005年）。对于每个窗口，使用方程（2）在高度（自变量）和MAAT（因变量）之间进行普通最小二乘回归：
(2)
其中表示温度变化，表示高度变化。这种局部化方法考虑了温度-高度关系的空间变化，这在复杂地形中尤为重要（Dobrowski等人，2009年）。为了确保露点率估计的统计可靠性，强制执行了至少五个有效像素对的最小样本阈值。

降尺度过程包括几个连续的步骤。首先，所有粗分辨率数据集（温度、露点率和高程）都使用双线性插值方法重采样到30米分辨率（Wilks，2011年，这是一种在地理空间分析中连续变量插值的标准技术）。考虑到方程（3），计算细分辨率的高程差ΔZ：
(3)
其中表示原始的30米SRTM DEM，表示重采样到30米的粗分辨率DEM。然后使用方程（4）计算温度调整项，作为重采样露点率和高程差的乘积：
(4)
最后，通过将此调整项加到重采样的粗分辨率温度场中，得到30米分辨率下的降尺度温度。
这种基于高程的校正方法基于环境露点率，通常在对流层中的平均值为约?6.5°C/km（Wallace和Hobbs，2006年）。该方法假设在从粗分辨率气候数据降尺度到细分辨率数据时，细空间尺度上的温度变化主要由高程差异解释。这一假设已在许多山地气候研究中得到验证，并已成功应用于永久冻土分布建模（Gruber，2012年；Ran等人，2022年）。
图S1显示了三个数据集降尺度后的MAAT的空间和时间变化，考虑了WorldClim和ERA5的历史分布以及CMIP6情景的未来预测。基于这三个数据集的KGRB的最大和最小温度范围显示在表S14中。

4.2. 错误评估
通过将降尺度后的温度场与水文和气象部（DHM）提供的1980年至2020年期间独立、质量控制的气象站数据进行验证，严格评估了其准确性。对于研究区域内的9个站点，将观测到的MAAT与每个十年的降尺度平均值进行了比较。性能使用均方根误差（RMSE，方程（5））和决定系数（R2，方程（6）进行量化：
(5)
(6)
其中是站点*i*处的观测MAAT，是降尺度后的MAAT，是观测值的平均值。

4.2. 机器学习模型
在本研究中，我们使用了三种机器学习模型：梯度提升（GB）、随机森林（RF）和支持向量机（SVM）来生成永久冻土概率分布图。我们使用了来自WorldClim、ERA5和CMIP6的降尺度MAAT数据，以及岩石冰川（RG）的起始点作为输入参数。评估永久冻土分布的关键步骤如图2所示。该过程从收集来自WorldClim和ERA5等来源的历史气候数据以及CMIP6情景的气候预测开始。随后，从尼泊尔喜马拉雅山脉特有的岩石冰川清单中提取IRGs和RRGs的坐标。

4.2.1. 梯度提升
梯度提升算法通过迭代方式将较弱的学习器（即略优于随机学习器）组合成强大的学习器（图S2）。梯度提升是一种类似于提升的回归算法（Friedman 2001）。给定一个训练数据集，梯度提升的目标是找到函数的一个近似值，该近似值将实例映射到它们的输出值，通过最小化给定损失函数的期望值来实现。梯度提升通过将函数的加权和构建函数的近似值：
(1)
其中是函数w的权重。这些函数是模型的集合（例如，决策树）。近似值是逐步构建的。首先，获得一个恒定的近似值：
(2)
后续模型旨在最小化：
(3)
然而，不是直接解决优化问题，每个模型可以被视为对进行梯度下降优化的贪婪步骤。为此，每个模型都在新的数据集上训练，其中伪残差通过以下方式计算：
(4)
然后通过解决线性搜索优化问题来计算的值。如果迭代过程没有适当规范，这种算法可能会过拟合。对于某些损失函数（例如，二次损失），如果模型完美地拟合了伪残差，则在下次迭代中伪残差变为零，过程提前终止。为了控制梯度提升的添加过程，考虑了几种规范参数。规范梯度提升的自然方法是应用收缩来减少每个梯度下降步骤的幅度。ν的值通常设置为0.1。此外，还可以通过限制训练模型的复杂性来实现进一步规范。对于决策树，我们可以限制树的深度或分割节点所需的最小实例数量。与随机森林不同，梯度提升中这些参数的默认值被设置为严格限制树的表达能力（例如，深度通常限制在3-5之间）。最后，梯度提升的不同版本还包括随机化基础学习器的参数家族，这可以进一步提高集合的泛化能力，例如不放回的随机子采样。

4.2.2. 随机森林
随机森林算法是机器学习中一种强大的树学习技术。它通过构建多棵决策树并合并它们的输出来进行预测。对于回归任务，使用这些树的平均值；对于分类任务，则选择从树中获得的票数最多的类别作为最终预测（图S3）。随机森林是一组树结构分类器，每棵树根据输入数据对类别进行投票（Breiman，2001）。
在本研究中，我们将控制参数设置为“重复cv”，进行多次训练分割，重采样迭代十次，并进行两次完整的折叠集重复。应用随机搜索程序来训练模型，以预测研究区域的永久冻土分布。

4.2.3. 支持向量机
支持向量机是一种主要用于分类任务的监督学习算法。它通过构建超平面或决策边界来将数据分为不同的类别（Cortes，1995）。SVM可以生成多个超平面来划分数据，确保每个数据段只包含一个类别的数据。SVM对线性和非线性数据都有效，使用核函数变换数据并创建决策平面。线性超平面的方程可以写为：
(5)
向量W表示超平面的法向量，即垂直于超平面的方向。方程中的参数b表示超平面沿法向量w的偏移量或距离。数据点与决策边界之间的距离可以计算为：
(6)
其中表示权重向量w的欧几里得范数。欧几里得范数的法向量对于线性SVM分类器：(7)优化：对于硬边线性SVM分类器：(8)受限于 = 1, 2, 3, …, m其中目标变量或标签在此语句中用符号表示。对于负例（当）时，ti=?1；对于正例（当）时，ti=1。因为我们要求决策边界满足约束条件：对于软边线性SVM分类器：(9)受限于 和 对于 = 1, 2, 3, …, m对偶问题：可以通过寻找与支持向量相关的拉格朗日乘数来解决SVM的对偶问题。最优拉格朗日乘数α(i)最大化以下对偶目标函数：(10)其中， 是与第i个训练样本相关的拉格朗日乘数。 是计算两个样本 和 之间相似度的核函数。它允许SVM通过将样本隐式映射到更高维的特征空间来处理非线性分类问题。项 表示所有拉格朗日乘数的总和。一旦解决了对偶问题并找到了最优拉格朗日乘数，就可以根据这些最优拉格朗日乘数和支持向量描述SVM的决策边界。训练样本中大于0的即为支持向量，决策边界由以下公式给出：(11)(12)在本研究中，使用10折交叉验证来训练SVM模型。使用Caret包中的train()函数实现带有线性核（svmLinear()方法）的SVM模型。在训练前，采取预处理步骤来中心化和缩放特征，确保所有变量的均值为0且标准差为1，这可以提高SVM的性能，尤其是在特征尺度不等的情况下。创建了一个包含20个`C`超参数值的调优网格，范围从0到2，以优化模型在数据集上的性能。我们实现了三种不同的机器学习算法GB、RF和SVM，以创建一个稳健的预测框架，用于映射研究区域的永久冻土分布模式。这种多算法方法提供了关于气候变化如何影响山区永久冻土稳定性的全面见解。4.3. 概率计算和模型训练我们的机器学习模型中的概率计算遵循二元分类方法，其中岩石冰川类型作为永久冻土条件的代理。完整的和残余的岩石冰川都被映射出来以创建训练数据集（图S4）。IRG被标记为表示永久冻土存在（Yes），而RRG被标记为永久冻土不存在（No）。这种二元标记系统（"Yes"/"No"）功能上等同于使用数值（1/0），其中"Yes"代表正类（永久冻土存在），"No"代表负类（永久冻土不存在）。机器学习模型使用岩石冰川的位置及其对应的MAAT值进行训练。在训练过程中，模型学习了最佳的分界线，以根据温度阈值将IRG（永久冻土条件）与RRG（非永久冻土条件）分开。在对新数据点进行预测时，机器学习模型输出的概率值范围从0到1，接近1的值表示永久冻土存在的可能性更高；接近0的值表示永久冻土存在的可能性更低。应用了0.6的阈值，高于此值的概率被认为是永久冻土存在的迹象。然后使用网格化的MAAT数据将训练好的模型应用于整个研究区域，生成一个永久冻土概率表面，每个像素根据其温度特征及其与IRG和RRG的相似性获得一个永久冻土概率值。实现这种概率计算的R代码包含在补充材料的附录1中。RRG是高山环境中永久冻土对气候变化响应的关键指标。岩石冰川存在于多种活动状态中，这对使用岩石冰川作为机器学习模型的训练观测具有重要的方法论意义。我们使用完整的岩石冰川通过持续的移动和蠕动过程来指示当前的永久冻土状况，而残余特征代表可能缺乏当前永久冻土的区域。这种区分对于训练机器学习模型至关重要，并为永久冻土概率映射提供了基本依据。鉴于地形复杂性和基于场地的调查有限的实际情况，这种方法已被广泛接受用于研究喜马拉雅地区的山地永久冻土（Jones等人，2018年；Baral等人，2020年；Baral和Haq，2020年）。假设研究区域在整个区域内具有相似的空间特征。有关训练数据集的空间和地形气候属性的更多信息，我们参考了Jones（2018年的研究），因为这些详细信息在此未包含。喜马拉雅地区RGs的起源和分布随经度而变化（Jones等人，2018年）。从映射的RGs派生的实证模型更适用于喜马拉雅山脉内的较小子区域（Haq和Baral，2019年；Baral等人，2020年）。用于模型训练的完整和残余的RGs都是使用WorldClim、ERA5和CMIP6情景数据集进行映射的。所有这三个气候数据集都必须经历训练和测试阶段（图S4），然后才能对永久冻土进行建模。4.4. 准确性评估我们使用RGs和三个再分析数据集作为永久冻土存在的参考指标来评估结果的准确性。我们采用了分层随机抽样方法，使用70%的IRGs、RRGs和MAAT数据用于模型训练，同时保留30%的数据用于测试。这种方法有助于我们评估模型在研究区域山区的准确性。我们使用混淆矩阵和接收者操作特征（ROC）曲线来衡量模型性能。在混淆矩阵中，分类结果被划分为真正例（TP）、真负例（TN）、假正例（FP）和假负例（FN），这些用于得出几个统计度量。基于这些结果，使用R编程计算关键指标，如准确性、置信区间、P值和Kappa统计量。准确性衡量正确分类实例占总实例的比例。它表示模型正确预测永久冻土存在和不存在的频率。ROC曲线展示了二元分类模型中不同阈值设置下灵敏度（真正例率）和特异性（真负例率）之间的权衡（Fawcett，2006年）。这些方法通过对RGs和MAAT参考数据的比较，提供了对永久冻土地图准确性和可靠性的关键见解。ROC曲线在不同阈值水平上绘制了真正例率（TPR）与假正例率（FPR）。TPR计算为TP/(TP + FN)，FPR计算为FP/(FP + TN)。理想的分类器将ROC曲线穿过左上角，表示完美性能，即TPR=1且FPR=0，意味着模型具有高灵敏度和特异性。ROC曲线下面积（AUC）通常用于评估整体模型性能，接近1的值表示强大的预测能力，接近0.5的值表明模型的预测能力与随机猜测相当。4.5. 模型配置和变量选择为了评估结合地形复杂性是否能够提高仅基于气候数据的永久冻土预测效果，我们比较了两种历史气候数据的模型配置。首先，我们仅考虑年平均气温作为永久冻土模拟的预测因子，并将其命名为“仅MAAT模型”。在第二次模拟中，我们考虑了包括MAAT、朝向、海拔和坡度的多变量，并将其称为多变量模型。对于这两种模拟，我们考虑了较早的（1961-1970年）和最近的（2010-2018年）时期。对于这些模型，MAAT是从WorldClim数据中降尺度到30米分辨率得到的，而地形变量是从30米DEM中使用R 'raster'包计算得出的，并将朝向转换为正弦/余弦分量，所有预测因子在训练前都进行了标准化，然后使用支持向量机分类方法进行了5折交叉验证，重复两次。5. 结果5.1. 平均温度的趋势分析首先，我们分析了历史和未来气候数据的温度趋势（图3）。历史年平均温度数据考虑了1961年至2020年的数据，未来预测数据取自2021年至2020年的时期。这些分析提供了观察到的和预测的气候变化的全面叙述，特别强调了变暖趋势。下载：下载高分辨率图像（600KB）下载：下载全尺寸图像图3. 来自多个数据集的历史和预测的年平均温度：(a) WorldClim再分析温度（1961–2020）；(b) ERA5再分析温度（1961–2020）；(c) CMIP6模型温度（1961–2014）以及根据SSP1-2情景的预测（2015–2020）；(d) 根据SSP5-8.5情景的CMIP6预测温度（2021–2020）。WorldClim数据（图3a）明显显示了变暖趋势，计算出的年变暖速率为0.019°C。尽管在60年期间存在自然年的变化，但整个时期的温度上升趋势非常明显。这种趋势的统计置信度异常强，表现为非常低的p值（p < 0.001），这提供了坚实的保证，表明观察到的变暖是一个真实的、统计上显著的趋势，而不仅仅是自然变化。R平方值为0.4847，表明这一时期年平均温度的约48%的变化可以直接由时间趋势解释。在这一历史时期，在KGRB中，最高的年平均温度记录为2016年（11.9°C），最低的温度记录为1971年（10.0°C）。对于相同的历史时期（1961-2020年），ERA5数据（图3b）显示出一个较温和但仍然明显的变暖模式，年变暖速率为0.01°C。与WorldClim数据相比，该数据集每年的波动更大，因此显示出较弱的线性相关性，其R平方值为0.12。尽管相关性较弱且波动较大，但趋势仍然具有统计学意义（p值 < 0.01）。ERA5数据集在此期间记录的最高温度为2016年（12.84°C），最低温度为1998年（10.21°C）。CMIP6历史模型显示从1961年到2020年的年变暖速率为0.03°C（图3c）。这种变暖模式并不均匀，早期几十年呈逐渐上升趋势，但2000年后增速显著加快，最后二十年温度增幅尤为明显。这一趋势的统计显著性非常高（p值 < 0.001）（p值 = 8.02×10^-18），并且其R平方值为0.72，表明超过72%的温度变化可以由时间趋势解释。展望未来，CMIP6情景对2021至2100年的预测显示了温度的急剧加速（图3d）。预测的年变暖速率为0.05°C，几乎是历史CMIP6情景速率的两倍，并且显著高于观察到的趋势。这一未来预测具有非常高的统计置信度，反映在p值（p < 0.001）和非常高的R平方值（0.9）上。虽然温度点在急剧上升的趋势线周围仍有一些变化，但年份与温度之间的线性关系非常强，突显了在本高排放情景下本世纪温度的显著加速趋势。总体而言，所有来源的温度数据都显示出显著的温度上升趋势。更重要的是，CMIP6情景的预测强调了21世纪温度增加的显著加速。历史数据和未来预测一致地证明了气候变化的持续和加速性质。5.2. 永久冻土分类的模型性能我们评估了三种机器学习算法GB、RF和SVM在三个气候数据集（WorldClim、ERA5、CMIP6）上对永久冻土概率的分类性能，评估了仅MAAT模型的稳定性和最佳配置选择。表1展示了基于WorldClim、ERA5和CMIP6（对于SSP5-8.5情景）的所有九种配置的机器学习模型（GB、RF和SVM）的全面性能指标。表1. 基于WorldClim、ERA5和CMIP6（对于SSP5-8.5情景）的机器学习模型（GB、RF和SVM）的永久冻土概率分类性能指标。模型配置准确性95% CIS灵敏度特异性精确度F1分数AUCWorldClim GB0.92(0.85, 0.96)0.990.380.920.950.91ERA5 GB0.91(0.84, 0.95)0.990.230.910.950.82CMIP6 GB0.89(0.81, 0.94)1.000.150.890.940.86WorldClim RF0.93(0.86?0.97)0.970.620.950.960.88ERA5 RF0.85(0.77?0.91)0.920.230.910.910.69CMIP6 RF0.86(0.77?0.92)0.910.540.930.920.81WorldClim SVM0.95(0.89?0.98)0.990.690.960.970.96ERA5 SVM校正后（τ*=0.38)0.81(0.75?0.87)0.860.750.770.810.82CMIP6 SVM校正后（τ*=0.56)0.78(0.73?0.85)0.900.690.740.810.86对三种机器学习算法（GB、RF和SVM）在多个气候数据集上的比较分析显示了显著的性能差异（表1）。基于WorldClim数据训练的SVM模型在性能上表现出色（准确率=0.95，AUC=0.96，F1=0.97），优于所有其他模型-数据集组合（表1）。WorldClim SVM在准确率上比第二好的配置（WorldClim RF）高出2.2%，在AUC上高出9.1%，同时保持了平衡的敏感度（0.99）和特异性（0.69）。尽管进行了阈值优化，气候数据集（ERA5，CMIP6）的AUC仍低14%–22%，这证实了空间分辨率从根本上限制了区分能力。WorldClim的精细空间尺度（30米）能够捕捉到对永久冻土分布至关重要的地形-气候梯度，而SVM基于核的学习方法有效地模拟了非线性的温度-永久冻土关系。狭窄的95%置信区间（0.89–0.98）表明模型的稳定性很高。WorldClim SVM代表了区域永久冻土制图应用的最佳配置。

**类别不平衡校正：** 初始的岩石冰川清单（408个样本：362个完整冰川，46个遗迹；比例为7.9:1）导致了严重的类别不平衡，使得SVM模型将所有样本都预测为完整冰川（特异性=0）。实施了两阶段校正协议：（1）合成少类过采样技术（SMOTE）来平衡类别（Chawla等人，2002年），以及（2）通过网格搜索优化概率阈值（0.30–0.70）。对于ERA5，最佳阈值为0.38，对于CMIP6为0.56，从而使F1分数最大化（分别为0.81和0.69）。这提高了特异性至0.75（ERA5）和0.69（SSP5-8.5），同时保持了高敏感度（>86%）和AUC > 0.82。校正后的模型能够可靠地识别无永久冻土区域（遗迹岩石冰川），从而能够绘制出由气候引起的退化情况。

**5.2.1 错误评估**
将网格化的气候数据集与实地站点观测结果进行验证，结果显示了在准确性和相关性方面的不同性能特征。WorldClim的均方根误差（RMSE）最低（RMSE=6.43 °C），表明其在再现站点位置观测到的温度幅度方面具有更高的准确性，同时解释能力适中（R2=0.551），解释了55.1%的观测温度变化（Fiddes等人，2015年）。ERA5再分析显示出较大的绝对误差（RMSE=9.76 °C），但相关性非常强（R2=0.913），解释了91.3%的温度变化，并显示出在系统性地低估温度幅度的情况下仍能很好地捕捉时空模式。CMIP6的未来预测显示出中等误差幅度（RMSE=8.27 °C），同时保持了类似的高相关性（R2=0.91），与ERA5的模式匹配能力相当。WorldClim较低的RMSE与ERA5/CMIP6较高的R2表明在绝对准确性和再现气候变异性的能力之间存在权衡：WorldClim提供了更好的点估计，而ERA5和CMIP6则更好地捕捉了潜在的气候动态和协变结构。这一发现证明了使用ERA5来理解气候模式和使用CMIP6进行未来预测的合理性，因为高R2值表明这些产品可靠地代表了相对于气候变化的关键差异和异常模式。

**5.3 模型性能：仅使用MAAT与多变量模型**
尽管地形控制具有理论相关性，仅使用MAAT的模型始终优于多变量模型。历史数据（1961–1970年）的仅使用MAAT的模型表现出出色的区分能力（AUC=0.91，准确率=92.0%，95%置信区间：0.80–0.97），具有高召回率（0.90）和中等精确度（0.60），F1分数为0.72。多变量模型的性能显著下降（AUC=0.755，准确率=87.62%，95%置信区间：0.79–0.93），召回率（0.80）、精确度（0.54）和F1分数（0.648）均有所降低，尽管增加了三个地形变量。在最近的时间段（2010–2018年），仅使用MAAT的模型保持了强劲的性能（AUC=0.869，准确率=89.52%，95%置信区间：0.82–0.94），召回率（1.00）完美，但精确度降低（0.50），再次优于多变量模型（AUC=0.750，准确率=87.62%，95%置信区间：0.79–0.93），AUC高出5.9个百分点和1.9个百分点（表S2）。

变量重要性分析显示预测因子贡献的时间变化显著。历史时期，MAAT（约35%）、方位（约25%）、海拔（约20%）和坡度（约15%）的重要性相对平衡。相比之下，最近时期MAAT占据了主导地位（约80%），而其他所有变量的贡献仅为约20%，这与历史时期的平衡分布形成了鲜明对比。时间分析揭示了不同的模式：仅使用MAAT的性能从历史时期到最近时期有所下降（AUC：0.914 → 0.809，准确率：92.0% → 89.52%），而多变量模型的性能在两个时期都保持相对稳定（AUC约0.75，准确率约87.6%）。地形变量的相对重要性从历史时期的约65%下降到最近时期的约20%（图S11）。

**5.4 永久冻土分布的历史变化**
我们使用GB、RF和SVM模型分析了1961年至2018年间KGRB山区永久冻土的时间和空间分布（图S12）。根据Obu等人（2019年）的永久冻土分类标准，永久冻土发生概率为60%–90%的区域被归类为不连续永久冻土，概率大于90%的区域被归类为连续永久冻土。GB模型显示每个十年永久冻土分布有轻微变化。值得注意的是，在研究期间的早期和最近几十年进行比较时，不连续永久冻土的分布明显增加。RF模型在整个六十年研究期间预测的永久冻土概率分布变化很小，包括早期和后期十年之间的比较。这种预测的稳定性可能与该模型的鲁棒性或其基础数据中使用的卫星图像的特定特性有关。

同样，我们使用WorldClim数据从SVM模型生成了两个时期的十年永久冻土概率分布图：1961?1970年（早期）和2010?2018年（近期）（图4）。在早期十年，可以看到大范围的连续永久冻土（深红色），尤其是在较高海拔地区，周围是不连续永久冻土区域（黄色轮廓）（图4a）。早期（图4a）和最新时期（图4b）的直接比较显示了一些明显变化，特别是在一个矩形黑框标记的区域，该区域显示了永久冻土范围的变化。观察到的永久冻土分布变化呈现出复杂且违反直觉的模式，挑战了传统的气候变化预期。连续永久冻土面积似乎有所扩大，而不连续永久冻土区域则有减少的趋势。这种变化表明，在1961?1970年至2010?2018年间，该山区发生了多种相互关联的过程。驱动这些变化的主要机制和因素包括海拔依赖的气候响应、水文制度变化以及地形和微气候控制，它们在不同空间和时间尺度上共同作用，导致了这种出乎意料的永久冻土演变。

**5.5 基于CMIP6情景的永久冻土分布预测变化**
我们使用GB、RF和SVM模型，并结合CMIP6（SSP-5.85情景）对2021年至2100年间的山区永久冻土概率分布进行了建模。每个模型生成了八张十年期地图，总共有24张预测永久冻土概率的图表（图S13）。所有模型都显示了高海拔地区永久冻土的显著存在。随着时间推移，永久冻土分布的模式显示出细微但有意义的变化。近期未来显示出更广阔的高概率永久冻土区域，而远期未来则表明这些区域逐渐减少。过渡边界随时间向上移动，表明随着气候变暖的影响，先前冻结的地面发生了永久冻土的退化。GB和RF模型仅预测到了地理区域内永久冻土分布的轻微变化。

**5.6 永久冻土面积的变化：历史趋势和未来预测分析**
历史数据集显示1961年至2018年间永久冻土动态的对比趋势（表2）。WorldClim数据显示永久冻土总面积净增加了119平方公里（1,802→1,921平方公里），其中连续永久冻土扩大了228平方公里（>90%概率），而不连续永久冻土减少了109平方公里（60%–90%概率）。相比之下，ERA5再分析报告总面积减少了724平方公里（2,717→1,993平方公里），连续永久冻土减少了1,199平方公里（2,279→1,080平方公里），而不连续永久冻土增加了475平方公里（438→913平方公里），表明稳定区域的广泛退化。两个数据集之间的843平方公里差异突显了历史监测的显著不确定性。表3显示了基于三种CMIP6排放情景（SSP1-2.6，SSP2-4.5，SSP5-8.5）从2021年至2030年到2091–2010年的永久冻土面积变化预测。所有模型都显示高海拔地区存在大量永久冻土。随着时间推移，永久冻土分布模式显示出细微但有意义的变化。近期未来显示高概率永久冻土区域更加广泛，而远期未来这些区域逐渐减少。过渡边界随着时间的推移向上移动，表明随着气候变暖，永久冻土退化加剧。GB和RF模型仅预测到地理区域内永久冻土分布的轻微变化。

**5.6 永久冻土面积变化：历史趋势和未来预测分析**
历史数据集显示1961年至2018年间永久冻土动态的对比趋势（表2）。WorldClim数据表明，永久冻土总面积净增加了119平方公里（1,802→1,921平方公里），主要由连续永久冻土扩大228平方公里（>90%概率）和不连续永久冻土减少109平方公里（60%–90%概率）驱动。相比之下，ERA5再分析报告总面积减少了724平方公里（2,717→1,993平方公里），其中连续永久冻土减少了1,199平方公里（2,279→1,080平方公里），而不连续永久冻土增加了475平方公里（438→913平方公里），表明稳定区域的广泛退化。表3总结了基于三种数据集的SVM模型预测的2021–2030年至2091–2010年的永久冻土面积变化。CMIP6预测显示，在所有情景下，永久冻土面积都有显著变化。

**图6 使用CMIP6 SVM模型进行的永久冻土指数时间和空间预测**
顶部行（图6a、6b和6c）分别显示了基于CMIP6数据的2021–2030年近期的永久冻土预测。在这三个情景下，所有情景都显示了相似的空间模式，其中大部分北部地区被深蓝色的连续永久冻土覆盖。这些连续永久冻土区域周围是不连续永久冻土的过渡带，用较浅的蓝色表示，形成了渐变效果。南部地区主要是非永久冻土地形，呈现为白色。

**图5 使用ERA5数据的相同情况**
相比之下，图5b（2010–2018年后期）显著说明了永久冻土的逐步退化，用矩形黑框表示。连续永久冻土（深红色）的范围显著减少，变得非常破碎。同时，“不连续”永久冻土（黄色）的区域显著扩大，占据了之前连续的区域。较高的海拔地区保持更稳定的永久冻土，而较低海拔地区的永久冻土退化更快。永久冻土的映射显示了该地区在大约50年分析期间快速的融化过程，这可能是由气候变化驱动的。然而，永久冻土的制图工作包括根据永久冻土概率水平和其他相关因素识别、分类或优先排序该研究区域内的不同区域。

**5.7 基于CMIP6情景的永久冻土分布预测变化**
我们还使用GB、RF和SVM模型以及CMIP6（SSP-5.85情景）对2021年至2100年间的山区永久冻土概率分布进行了映射。每个模型生成了八张十年期地图，共形成了24张预测永久冻土概率的图表（图S13）。所有模型都显示高海拔地区存在大量永久冻土。随着时间的推移，永久冻土分布的模式显示出细微但有意义的变化。近期未来显示出更广阔的高概率永久冻土区域，而远期未来则表明这些区域逐渐减少。过渡边界随时间向上移动，表明随着气候变暖，先前冻结的地面发生了永久冻土的退化。根据SSP1-2.6（低排放）的预测，永久冻土基本上保持稳定，仅略有增加9平方公里，而且连续永久冻土的增加抵消了不连续永久冻土的减少。SSP2-4.5（中等排放）预测永久冻土将减少179平方公里，而SSP5-8.5（高排放）则预测永久冻土将灾难性地减少822平方公里，连续和不连续区域都将严重缩减。高排放的影响大约是SSP1-2.6下的90倍。表2显示了1961-1970年至2010-2018年间不同数据集和类别下的永久冻土概率面积的历史变化。

永久冻土类别 | 数据集 | 1961-1970（平方公里） | 2010-2018（平方公里） | 变化（平方公里）
--- | ------------------- | --------------------- | ------------------------ | -------------
Discontinuous (60%–90%) | WorldClim SVM | 553 | 444 | -109 |
Continuous (>90%) | WorldClim SVM1,2 | 491 | 1,477 | 228 |
Total | WorldClim SVM1,8 | 1,802 | 1,921 | 119 |
Discontinuous (60%–90%) | ERA5 SVM | 438 | 913 | 75 |
Continuous (>90%) | ERA5 SVM2,2 | 2,279 | 1,080 | -1,199 |
Total | ERA5 SVM2,7 | 1,717 | 1,993 | -724 |

注：负数值表示面积减少。SVM = 支持向量机。

表3显示了CMIP6情景下2021-2030年至2091-2100年永久冻土概率面积的预测变化。

永久冻土类别 | 情景 | 2021-2030（平方公里） | 2091-2010（平方公里） | 变化（平方公里）
--- | ------------------- | ------------------- | ------------------------ | -------------
Discontinuous (60%–90%) | SSP1-2.6 | 742 | 669 | -73 |
SSP2-4.5 | 741 | 683 | -58 |
SSP5-8.5 | 790 | 385 | -405 |
Continuous (>90%) | SSP1-2.6 | 1,193 | 1,275 | 82 |
SSP2-4.5 | 1,186 | 1,065 | -121 |
SSP5-8.5 | 1,218 | 801 | -417 |

注：负数值表示面积减少。SSP = 共享社会经济路径。

表4总结了所有模型和时间段内的永久冻土面积趋势。

模型 | 时间段1 | 时间段2 | 总变化（平方公里） | 趋势 |
--- | -------- | ------------ | ---------------- | ------------ |
WorldClim SVM | 1961–1970 | 2010–2018 | 119 | 增加 |
ERA5 SVM | 1961–1970 | 2010–2018 | -72 | 减少 |
SSP1-2.6 SVM | 2021–2030 | 2091–2010 | 9 | 轻微增加 |
SSP2-4.5 SVM | 2021–2030 | 2091–2010 | -179 | 减少 |
SSP5-8.5 SVM | 2021–2030 | 2091–2010 | -822 | 减少 |

注：CMIP6的预测在SSP1-2.6、SSP2-4.5和SSP5-8.5排放情景下使用了SVM方法。

通过比较各数据集的基本特征，系统地研究了这种差异。WorldClim是一个基于气象站数据的空间插值产品，可能会导致气候场变得平滑，并可能低估复杂地形中的局部和高海拔地区的变暖趋势。相比之下，ERA5是一种高分辨率再分析方法，它融合了更多的观测数据（包括卫星检索数据），可能更好地捕捉局部大气过程和近期的变暖加速。这种方法论上的差异解释了历史数据的不一致性。重要的是，两组数据都显示出一个一致的空间模式：永久冻土逐渐集中在更高的海拔（5,000–6,000米），识别出相同的气候避难区。

最显著的是，尽管历史数据存在分歧，所有未来的CMIP6预测都表明永久冻土会减少（从SSP1-2.6下的近乎稳定到SSP5-8.5下的严重退化）。这一关于未来变化方向的共识，经过与关于喜马拉雅山脉变暖的广泛文献的验证，为永久冻土未来显著脆弱性的核心结论提供了高度信心。

为了全面了解永久冻土动态，特别是在变化的气候条件下，分析多个模型和数据集的预测至关重要。根据永久冻土模型的评估，SVM展示了最佳的结果。图7整合了详细的十年条形图，展示了不同概率范围内的永久冻土分布。通过研究这些图表，我们可以辨别出一致的趋势，识别特定模型的敏感性，并理解永久冻土制图和预测中的复杂性。

下载：下载高分辨率图像（523KB）
下载：下载全尺寸图像

图7. 使用三个网格化数据集（a）WorldClim、b）ERA5和c）CMIP6（SSP5-8.5）的SVM模型得出的永久冻土面积（平方公里）及相关概率分布。每个条形代表给定十年期间的总永久冻土面积，堆叠的段表示该面积落在特定永久冻土存在概率范围（0.5-0.6、0.6-0.7、0.7-0.8、0.8-0.9和0.9-1.0）内的比例，分别用浅蓝色、红色、灰色、黄色和蓝色表示。历史数据（1961-2018年，图7a）和ERA5数据（图7b）显示，在六十年期间永久冻土面积基本稳定。WorldClim数据显示永久冻土面积从2,300平方公里略微减少到2,048平方公里，而ERA5数据显示永久冻土覆盖面积保持在2,200-2,300平方公里左右不变。显著的绿色段代表高概率永久冻土（0.9-1.0范围），表明这些地区在所有历史时期都保持热稳定性。CMIP6预测显示未来永久冻土面积将加速减少，从历史水平约2,100平方公里减少到本世纪末的1,262平方公里。2050年后，减少速度显著加快，特别是在最后几十年。重要的是，即使是历史上最稳定的高概率永久冻土区域（绿色区域）也经历了显著收缩。这种变化标志着研究中的冰冻圈的根本转变。

图8显示了由WorldClim、ERA5和CMIP6（SSP5-8.5情景）气候数据驱动的永久冻土浓度随海拔的分布。在所有三个数据集中，永久冻土都一致集中在大约5,000米以上的高度，表明无论数据来源或时间段如何，这些高海拔地区的永久冻土都非常稳定。在这些高度上，不同数据集之间的差异相对较小，表明地形和热控制因素占主导地位。

相比之下，在较低海拔（<5,600米）地区，特别是在早期时期（1961-1970年），当比较基于WorldClim和ERA5的SVM结果时（图8a和8b），出现明显差异。CMIP6对未来早期的预测（2021-2030年；图8c）与近期时期的ERA5 SVM分布（2010-2018年；图8e）更为一致，而与WorldClim SVM结果（图8d）相比差异更大，尤其是在5,000米以下。此外，ERA5 SVM结果显示近期时期（图8e）5,500米以下的永久冻土出现显著减少，而早期时期（图8b）则没有这种减少。4,500-5,500米高度范围内的永久冻土减少尤为明显，指出这个范围是永久冻土对气候变暖最敏感的关键过渡区。

图9展示了使用SVM模型根据三个气候数据集（WorldClim、ERA5和CMIP6（SSP5-8.5情景）得出的永久冻土浓度随海拔的分布。在所有三个数据集中，永久冻土都一致集中在大约5,000米以上的高度，表明无论数据来源或时间段如何，这些高海拔地区的永久冻土都非常稳定。在这些高度上，数据集之间的差异相对较小。

图5.8分析了岩石冰川与永久冻土概率的关系。图9展示了使用SVM模型根据三个气候数据集得出的IRGs（岩石冰川）和RRGs（残留冰川）起始点的永久冻土概率（PP）分布。在历史时期，WorldClim（1961-1971；图9a）和WorldClim（2010-2018；图9d）的模拟显示IRGs集中在最高的PP类别（PP=0.9–1）中，而RRGs主要出现在永久冻土地带之外，仅在较低的PP范围（PP < 0.8）中有少量存在。在1961-1970年（图9b）和2010-2018年（图9e）的ERA5中观察到类似的模式。

图9比较了我们的RG清单与已建立的区域清单，即Schmid等人（2015年）的HKH清单和Harrison等人（2024年）的最新高分辨率喜马拉雅清单，发现了两者之间的一致性和显著差异。在区域尺度上，我们的清单在A区域内识别出30个RG，比Harrison等人（2024年）记录的50个RG少20个。这种差异可能是由于空间划分的不同、我们在制图时采用的更严格的形态学标准，或是参考研究中使用的新型卫星数据具有更强的检测能力所致。在样本层面，验证显示出了混合的结果。样本A在我们的研究和Schmid等人（2015年）的研究中都一致记录为包含一个RG，这证实了对于定义明确的特征的检测可靠性。相比之下，样本B在我们的分析中被分类为不包含RG，而在之前的清单中被记录为包含一个RG。这种差异可能反映了自2015年以来的地貌变化，或是基于更严格分类标准的修正。总体而言，这种交叉参考强调了我们在识别清晰RG特征方面的方法一致性，同时也指出了方法和时间因素对清单完整性的影响。

表5总结了与我们之前的清单相比的结果。

研究区域 | 本研究 | 之前的清单 | 差异 |
--- | ------- | -------------- | ------------ |
Region A | 否 | 30 | |
No. of RG: | | 50 | |
Sample A | | 含有RG | |
Doesn't contain | | 不含有RG | |
Sample B | | 含有RG | |

表6显示了单因素方差分析（ANOVA）的结果，表明气候数据集的选择对研究范围内的永久冻土概率估计有显著影响。分析得出的F统计量很大（F=994,866；p < 0.001），效应量中等（η2=0.066），表明大约6.6%的永久冻土概率总变异可以归因于数据集之间的系统差异。尽管这一比例看起来很小，但其在整个超过1380万像素的领域内的空间表现导致了估计的永久冻土范围的显著差异。

表6显示了来自WorldClim和ERA5数据集（1979-2020年）的永久冻土概率估计的统计摘要。结果基于单因素方差分析和变异系数（CV）分析。

世界Clim和ERA5都表现出明显的变异性，尽管其分散程度明显不同。WorldClim的变异系数非常高（CV=136%），而ERA5的变异系数相对较低（CV=36%）。这些数值不应单独解释为数据质量的指标，而应反映每个数据集的固有特征和尺度。WorldClim的CV超过100%是由于永久冻土景观的固有异质性，其中概率值介于接近二元条件（存在与否）之间。在这种分布中，当聚合多样化的地形单元时，标准差可能会超过平均值，这在空间异质性环境变量中是常见的现象（Kelley, 2007）。相比之下，ERA5的CV较低，反映了其更平滑的空间场和更狭窄的值范围，这可能是由于其较粗的分辨率和数据同化框架所致。因此，关键结果不是CV的绝对大小，而是数据集之间的系统平均偏差。WorldClim和ERA5的平均偏差为0.448个概率单位（p < 0.001），这种差异在统计上是稳健的，但相对于完整的概率范围（0–1）来说仍然是适度的。这表明数据集引入的偏差是存在的，但可以通过偏差校正和交叉验证程序有效缓解。这些发现为我们的双数据集策略提供了定量支持：WorldClim具有精细的空间分辨率（约1公里），能够捕捉长期气候学映射所需的详细空间梯度；而ERA5的变异性较低，且在表示年际气候变异的短期永久冻土响应方面表现出色（R2=0.913）。WorldClim的均方根误差（RMSE）较低（6.4°C），相比之下ERA5的RMSE为9.8°C，这加强了其作为永久冻土基准映射工具的适用性，尽管其空间变异系数（CV）较高。空间CV较高但验证性能较好的这一表面矛盾可以通过以下方式解释：CV量化的是整个区域的空间异质性，而RMSE衡量的是观测点的精度。因此，一个数据集可能同时表现出较高的空间变异性——这表明它能够代表多样的永久冻土状况——以及在验证地点的高精度。先前的研究已经证明，多模型或多数据集集成可以减少气候预测的不确定性（Tebaldi和Knutti，2007年）。然而，我们的结果表明，如果不考虑数据集的系统性偏见而简单地平均多个数据集，可能会导致人为平滑或扭曲永久冻土边界。6.6%的方差由数据集选择解释，这突出了三个关键的方法论含义：（1）数据集的选择必须基于定量性能指标，而不仅仅是便利性或惯例；（2）在将气候数据应用于永久冻土影响模型之前，必须进行偏差校正；（3）不确定性评估应明确纳入数据集特定的系统误差以及自然气候变异性。在这个框架下，WorldClim较高的CV不应被视为不可靠，而应被视为其能够代表从连续永久冻土（概率>0.8）到非永久冻土地形（概率<0.1）整个连续体的能力。

讨论
本研究利用环境变量指导的SVM模型，提供了关于尼泊尔KGRB地区永久冻土分布的新见解。结果表明，当将遥感（RG）清单与来自多个来源的地面空气温度（MAAT）数据结合时，机器学习方法可以有效捕捉高山地区的永久冻土存在概率。这种方法与早期研究结果一致，这些研究强调了在永久冻土建模中纳入环境预测因子的重要性（Jones等人，2018年；Baral等人，2020年；Pradhan和Shukla，2022年）。
在KGRB地区，海拔5,000至6,000米之间形成了一个明显的永久冻土集中区，作为一个独特的永久冻土避难所。这种模式反映了多种环境控制因素的共同作用，而不仅仅是海拔依赖的冷却效应。在这个范围内，MAAT值（-2至-5°C）处于永久冻土稳定性的最佳范围内（Gruber，2012年）。在海拔5,000米以下，大气变暖超过了地热系统的缓冲能力；而在海拔6,000米以上，虽然温度适宜，但稳定地形的可用性急剧下降。太阳辐射在这些海拔高度调节地温方面起着关键作用。朝北的斜坡接收到的年辐射显著减少（Gruber和Hoelzle，2001年；Hock，2003年），这种效应被周围海拔超过7,000米的山峰（如Dhaulagiri和Annapurna）的地形遮挡所放大。这种减少的日照有效地降低了地热系统，使得在其他条件下原本边缘化的地区也能维持永久冻土。雪的重新分布进一步增强了热调节作用：在5,000-6,000米范围内的背风坡积累了0.5-2米的雪，使地面温度降低了约2-4°C。在更高、暴露于风的海拔高度（>6,000米），有限的雪层保留减少了这种保护效果。
从地形上看，这个海拔区间为永久冻土的保存和岩冰川的发展提供了最佳条件，其特征是坡度适中（10°-35°）、碎屑供应充足以及以冰斗为主的地貌（Schmid等人，2015年）。较低的海拔缺乏必要的热条件，而较高海拔地区更陡峭且更不稳定的地形限制了碎屑的积累和冰层形成。热适宜性、减少的辐射、雪的保温作用以及有利的地貌相互结合，形成了一个狭窄的海拔“最优点”，最大化了永久冻土的存在。在喜马拉雅山脉的其他地方也报告了类似的中等海拔浓度模式（Schmid等人，2015年；Bhat等人，2025年），支持了这一机制框架的普遍适用性。
东西向的岩冰川特征梯度进一步反映了尼泊尔的地形降水模式。向西部尼泊尔移动的雨影效应在KGRB和Api Nampa地区创造了更大陆性的气候条件，有利于永久冻土的保存，超过88%的岩冰川保持完整。相比之下，东部地区（如Solukhumbu）增强的季风湿度加速了永久冻土的退化，几乎一半的岩冰川已经消失。KGRB（5,000-6,000米）和Solukhumbu（3,893-5,446米）之间永久冻土区1,100米的海拔差异对应大约6.6-7.7°C的温度差异，反映了在湿润条件下永久冻土等温线的下降。这些区域对比在表7中进行了总结，强调了针对特定气候进行模型校准的必要性，而不是依赖统一的海拔阈值来进行全国范围的永久冻土评估。

表7. 尼泊尔各喜马拉雅地区岩冰川分布的比较分析
尼泊尔地区 海拔范围 RG面积 RGD密度 RGIRGRRG
KGRB（中部；Jones等人，2018年） 6,400–8,167 5,000–6,000 192 0.21 362
Api Nampa保护区（西部；Ghimire和Ghimire，2024年） 低于2,000–8,848 上升5,000 28.1 10.24 109
Solukhumbu（东部；Singh等人，2025年） 855–7,106 3,893–5,446 NA 897 8

模型敏感性实验显示，仅使用MAAT的模型表现优于多变量配置，这可能反映了30米分辨率下的信息冗余。降尺度的气候数据集隐含了海拔依赖的递减率和与朝向相关的辐射差异（Karger等人，2017年），减少了额外地形变量的增量效益。从较冷历史条件下相对平衡的热控制和地形控制到最近几十年MAAT的强主导地位，变量重要性的变化表明了永久冻土分布向气候控制的根本转变。当区域MAAT接近临界阈值（-1至0°C）时，即使是有利的地形条件也不足以维持永久冻土（Gruber和Haeberli，2007年；Noetzli等人，2007年）。仅使用MAAT的模型具有高召回率但精度较低，这突显了其在变暖情景下检测脆弱永久冻土的实用性（Huang等人，2025年）。
与现有的全球和地区永久冻土数据集的比较进一步丰富了我们的发现。永久冻土分区指数（PZI；Gruber，2012）估计尼泊尔有8%–10%的地区处于连续永久冻土区，而我们的流域尺度分析表明KGRB有15%–17%的地区为连续永久冻土。这一较高的比例反映了该流域相对于全国平均水平的高海拔和气候优势（Obu等人，2019年；Ran等人，2022年）。与实地观测的验证支持了我们建模框架的稳健性：基于WorldClim的模拟显示的RMSE（6.43°C）低于ERA5（9.76°C），尽管其空间变异性较高。这一表面矛盾的出现是因为变异系数反映了整个区域的空间异质性，而RMSE衡量的是观测点的精度。

使用单因素方差分析（one-way ANOVA）明确量化了数据集选择的统计影响，发现WorldClim和ERA5得出的永久冻土概率估计之间存在显著差异（F=994,866.2，p<0.001），数据集选择解释了总方差的相当一部分（η2=0.066）。这些系统差异强调了选择无偏见数据集的重要性。与此解释一致的是，不同数据集和时间段内永久冻土概率与空间范围之间的关系仍然稳健，如图10所示，更高的概率阈值始终对应于WorldClim、ERA5、CMIP6和先前发布的数据集中更小、更持续的永久冻土面积。

图10. KGRB不同气候数据集和时期下永久冻土存在的概率与覆盖面积的关系。各图表显示了不同数据集的空间覆盖范围（以平方公里计）与永久冻土存在概率之间的关系：(a) WorldClim. 2018，(b) ERA5. 2018，(c) CMIP6. scenario 2100，(d) Gruber. 2012，(e) Obu. 2019，以及(f) Ran. 2022。每条曲线都展示了随着空间范围的增加和不同永久冻土数据集建模方法的改变，永久冻土概率的变化。

未来的CMIP6预测表明，在变暖情景下，永久冻土将逐渐收缩到5,000–6,000米的避难区，并且完整的岩冰川会重新分配到概率较低的类别中，而在SSP5-8.5强迫下，模型中的遗迹岩冰川将消失。这种非线性反应反映了在低水平变暖下的热惯性和潜热缓冲，一旦超过临界阈值，退化就会加速（Cui等人，2025年；Ji等人，2025年）。随着变暖加剧，永久冻土越来越局限于高海拔避难区，使其对进一步的气候变化更加脆弱。

总体而言，我们的结果强调了在喜马拉雅山脉进行持续永久冻土监测和区域校准建模的紧迫性。通过明确联系统计模型性能、物理过程理解和区域气候变化梯度，本研究完善了现有的基于海拔的喜马拉雅永久冻土概念，并强调了高海拔避难区在持续和未来气候变化中维持永久冻土的关键作用。

结论
本研究使用由多种气候数据集驱动的机器学习模型，评估了尼泊尔Kali Gandaki河流域（KGRB）的历史和未来永久冻土动态。在测试的算法中，支持向量机（SVM）在利用WorldClim数据时表现出最强的永久冻土映射性能（AUC=0.974）。历史分析显示，1961-2020年间年平均温度上升了0.011–0.019°C/年，预计到二十一世纪末将加速到约0.053°C/年。
岩冰川分析表明，完整的岩冰川与高永久冻土概率区密切相关，证实了它们作为永久冻土存在的地貌指示器的价值。然而，预计持续的气候变暖将推动从完整状态向遗迹状态的逐步转变，岩冰川将越来越多地占据中等概率区域，而不是持续永久冻土区域。历史永久冻土重建显示出相当大的不确定性，主要是由于输入气候数据的差异。基于WorldClim的模拟表明，1961年至2018年间永久冻土范围相对稳定或略有扩大（+119平方公里），而基于ERA5的结果则表明同一时期永久冻土明显退化，净损失约为724平方公里，主要由于连续永久冻土的减少。这些差异突显了回顾性永久冻土估计对气候数据选择的敏感性，并强调了在数据稀少的山区重建过去永久冻土动态的挑战。
尽管存在这些历史上的不一致性，未来的预测在不同数据集和模型之间表现出强烈的一致性。永久冻土损失与温室气体排放路径直接相关，从SSP1-2.6下的近乎稳定（+9平方公里）到SSP5-8.5下的严重退化（-822平方公里）在2021-2100年间。总体而言，这些预测表明在高温排放情景下KGRB的永久冻土将显著退化，强调了在高山环境中采取适应和风险缓解策略的紧迫性。
在国家层面，尼泊尔的永久冻土稳定性表现出明显的东西向梯度，主要由气候湿度制度控制，而不仅仅是海拔。KGRB保留了广阔的完整永久冻土，是一个重要的高海拔水和冰库，而以季风为主的东部地区显示出严重的退化，预示着即将失去冰层。因此，有效的国家层面评估和适应需要特定于气候的建模框架和协调的区域校准监测协议。
总体而言，本研究通过结合机器学习技术和多源气候数据，推进了喜马拉雅永久冻土的研究，解决了历史不确定性问题和未来的脆弱性。虽然有针对性地使用MAAT和岩冰川分布提高了模型的可解释性，但也突出了改进观测约束的必要性。未来的努力应包括扩展的岩冰川清单、实地地温测量和钻孔测量，以及额外的环境控制因素，如海拔、朝向、土壤特性和水文指数（例如NDWI）。在过渡性永久冻土区进行有针对性的实地验证对于提高模型稳健性和增强对正在进行的气候变暖下永久冻土变化的预测信心尤为重要。

作者贡献声明
Earina Sthapit：编写——原始草案，可视化，形式分析，数据 cura，概念化。
Binod Pokharel：编写——审阅与编辑，监督，项目管理，方法论，调查。
Dibas Shrestha：编写——审阅与编辑，监督，项目管理。
Deepak Aryal：编写——审阅与编辑，监督。
Arnab Singh：编写——审阅与编辑，软件，方法论，数据 cura。
Anita Tuitui：编写——审阅与编辑，软件

未引用的参考文献
Brenning等人，2005年；Cortes和Vapnik，1995年；Eyring等人，2016年；Fiddes和Gruber，2015年；Hoelzle等人，2001年；Ishikawa等人，2003年；Mayer等人；Mauritsen和Roeckner，2020年；Obu，2021年；Olson和Rupper，2019年；Shekhar等人，2010年；Singh等人，2025年。

利益冲突声明
作者声明没有已知的竞争性财务利益或个人关系可能影响本文报告的工作。