1 引言

青藏高原作为全球气候研究的热点区域，冻土分布广泛，其土壤冻融过程对水文、生态、气候和农业等诸多领域产生了重要影响[1-2]。近年来，受气候变暖影响，地表温度升高，部分多年冻土退化，导致地下冰融化，地下水位下降，某些寒区植被无法充足吸收水分，造成植被退化，进而对整个生态系统造成严重威胁[3-4]。同时冻土中含有大量有机碳，土壤冻融交替导致更多的温室气体释放，进而加剧气候变暖[5]。此外，频繁的冻融交替改变了土壤肥力和结构，影响农作物产量[6]。并且土壤松散，土壤稳定性下降，滑坡、泥石流等地质灾害概率增加[7]。因此，研究青藏高原的土壤冻融过程对理解该地区生态系统与气候变化的相互作用、碳循环过程以及对气候变暖的响应机制具有重要意义。

目前，关于土壤冻融过程研究通过站点观测、数值模拟和遥感反演3种方法进行。王康等[8]利用中国境内845个气象站数据，以冻结天数为指标，探讨其时空分布特征及其与气候变化的关系。杨淑华等[9]利用87个气象站的逐日地表温度数据，将土壤冻融状态分为3种状态，分析了1980—2015年青藏高原不同状态的时空变化并探讨了和气温的关系。站点观测虽能提供高精度数据，但受站点数量和分布不均的限制，难以全面反映大尺度区域的冻融状况[10]。为了获取研究区土壤冻融状态的空间分布特征，夏坤等[11

通过使用陆面过程模式CLM3.0，模拟了青藏高原东北部季节性冻土区域，结果显示CLM3.0模式模拟冻结过程好于融化过程。刘火霖等[12]利用高寒气候环境观测研究站资料，对冻融参数进行优化，通过陆面模式CoLM更好地模拟了那曲地区土壤冻融过程特征。数值模拟能够更好地模拟土壤冻融空间分布特征，但模拟参数较多，不同参数会直接影响冻融判别的精度[13]。近年来，遥感技术也被广泛应用于表层土壤冻融过程的研究。Hu等[14]基于MODIS和AMSR-2遥感数据进行融合，分析研究了中国区域2013和2014年的土壤冻融情况，生成5 km分辨率冻融图。沈麒凯等[15]基于MODIS、AMSR-E及AMSR-2遥感数据，使用Fourier非线性模型进行数据融合，生成了1 km高分辨率冻融结果。徐富宝等[16]基于AMSR-2亮温数据和土壤湿度数据，使用LightGBM算法和随机森林算法对2014年7月—2015年6月青藏高原4个典型区域进行土壤冻融监测。遥感技术虽能弥补站点观测的不足，但受空间分辨率和混合像元等因素影响，结果仍存在不确定性[17]。

ERA5-LAND数据集是由多种观测资料进行数据同化得到的再分析产品，其时间序列长、不受地形限制，且经过质量控制，是研究土壤冻融过程的理想数据源[18]。青藏高原作为全球高海拔冻土分布最广泛的地区，其土壤冻融过程研究更具科学意义和实践价值。目前多以站点观测数据研究青藏高原局部地区土壤冻融时空变化，针对整个青藏高原地区的研究较少且仅探讨单因子对土壤冻融过程的影响，在驱动因子对土壤冻融过程影响程度的定量化研究存在不足。因此，本文选择以青藏高原为研究区，采用2003—2022年ERA5-LAND地表温度数据来研究青藏高原土壤冻融过程，选用冻结起始时间、冻结结束时间、冻结持续时间及冻结天数4种参数，探讨了4种参数的时空变化特征，并结合地理探测器分析青藏高原不同区域冻结天数空间分异的主导因素，最后使用相关性分析法分析冻融参数与主要驱动因子的相关性，以期为青藏高原对气候变化的响应提供理论依据。

2 研究区概况、研究方法与数据来源

2.1 研究区概况

青藏高原作为中国最大、世界海拔最高的高原，被称为"世界屋脊"、"第三极"，平均海拔超过4 000 m，地理位置如图1（a）所示，覆盖了西藏、青海、四川、新疆、甘肃和云南等部分地区，其基于地理空间数据云（https://www.gscloud.cn/）提供的2000年SRTMDEM 90 m分辨率原始高程数据制作。青藏高原冻土分布广泛，如图1（b）所示，由国家青藏高原科学数据中心（https://data.tpdc.ac.cn/）提供的青藏高原新绘制冻土分布图（2017）[19]制作，主要分为季节冻土和多年冻土。近年来，伴随着全球变暖的趋势不断增加，青藏高原多年冻土已然发生了极为显著的退化现象，具体呈现为地表温度升高、活动层厚度增加、多年冻土面积不断减少等[20]。

青藏高原研究区概况

图1 青藏高原研究区概况

2.2 研究方法

2.2.1 技术路线

本文技术路线如图2所示。首先进行数据收集与处理，通过ERA5-LAND地表温度数据预处理获取青藏高原逐日最高最低温度数据，站点观测数据对其进行验证。随后通过土壤冻融阈值划分获取表征土壤冻融过程的冻融参数，并与经过预处理的驱动因子数据进行坐标系匹配、分辨率匹配以及像素匹配；通过空间分布状况、像元面积占比和变化速率表征青藏高原土壤冻融时空变化特征；通过地理探测器将冻结天数与驱动因子数据结合进行驱动力分析对土壤冻融参数定量化研究，最终根据得到的主导因素结合冻融参数进行相关性分析。

青藏高原土壤冻融参数时空变化及驱动因子定量分析技术路线

图2 青藏高原土壤冻融参数时空变化及驱动因子定量分析技术路线

2.2.2 土壤冻融阈值划分

为使一个完整冻结期包含在内，同时根据青藏高原特有地理条件及气候情况，本文选取当年7月1日—次年6月30日定为一个循环年。如2003年7月1日—2004年6月30日记为2003年的数据，依此获取2003—2022年共20年数据资料。

本文参考杨淑华等[21]研究选取冻结起始时间、冻结结束时间、冻结持续时间和冻结天数4种参数表征土壤冻融过程。定义方法参考岳书平等[22]及薛华柱等[23]研究，对数据进行5 d滑动平均处理，且以0 °C定为土壤冻融划分的阈值。在一个冻融年内，当首次出现日最低土壤温度在0 °C及以下时，视为冻结起始时间；最后一次出现该种情况时，视为冻结结束时间；冻结持续时间即为前2种时间之间的天数；而冻结天数指的是整个冻融年内实际处于冻结状态的天数，因之间存在部分日期最低温度回升到0 °C以上的情况，所以一般冻结持续时间的值大于冻结天数。此外为了方便统计，本文以7月1日为第一天，7月2日为第二天，以此类推。

2.2.3 Theil-Sen斜率估计法

Theil-Sen斜率估计法是一种常用的非参数统计方法，其计算效率高，对测量误差和利群数据不敏感，常与Mann-Kendall检验法结合用于进行长时间序列的趋势分析[24]。计算公式如下：

β = Median( (xj - xi) / (j - i) ) ∀j > i

式中：xj和xi分别代表第j年和第i年的冻融参数，且j > i；β斜率代表冻融参数的变化趋势；Median代表取中值。若β > 0，则代表上升趋势；β < 0，则代表下降趋势。

2.2.4 Mann-Kendall检验法

Mann-Kendall检验法又称为M-K检验，也作为一种非参数统计方法，可以进行时间序列数据的趋势分析和突变点检验，该方法不需要数据满足特定的分布，且不易受到缺失数据与异常数据的干扰，因此广泛应用于长时间序列的趋势分析[25]。通过比较序列中数据对的大小关系,构建统计量S，进而判断序列是否存在显著趋势。统计量计算公式如下：

S = ∑i=1n-1 ∑j=i+1n sgn(xj - xi)

式中：xj和xi分别代表第j年和第i年的冻融参数，sgn为符号函数，计算公式为：

sgn(xj - xi) = { 1 if xj - xi > 0; 0 if xj - xi = 0; -1 if xj - xi < 0 }

使用检验统计量Z进行趋势检验，计算方法：

Z = { S/√Var(S) if S > 0; 0 if S = 0; (S+1)/√Var(S) if S < 0 }

式中：Var为随机变量的方差函数，其计算公式为：

Var(S) = n(n-1)(2n+5)/18

式中：n表示数据个数。使用双边趋势检验和临界值Z1-α/2来判断趋势显著性。在显著性水平α=0.05下，若|Z|大于1.96，代表趋势显著，并通过95%的显著性检验，具体判断方法见表1。

表1 2003—2022年青藏高原土壤冻融参数趋势类别及显著性划分

β Z 趋势类别 趋势特征
β > 0 1.96 < Z 2 显著增加
β > 0 Z ≤ 1.96 1 不显著增加
β = 0 Z 0 无变化
β < 0 Z ≤ 1.96 -1 不显著减小
β < 0 1.96 < Z -2 显著减小

2.2.5 地理探测器

地理探测器是一种探测数据空间分异的统计方法[26]。本文使用地理探测器的单因子探测及交互因子探测方法分析青藏高原及不同冻土区域冻结天数空间分异的驱动力。由于数据需要进行离散化处理，本文使用自然间断法结合不同数据信息特征，将其重分类为5~9类。地理探测器通过q值度量空间分异，取值范围为[0，1]，其值越大代表自变量因子对研究区冻结天数空间分异的解释力越强。

2.2.6 相关性分析法

相关分析法体现了2个因素之间的关联程度[27]。本文基于像元尺度的4种参数分别与主要驱动因子进行相关性分析，相关系数的计算公式如下：

Rxy = ∑i=1n[(xi - x̄)(yi - ȳ)] / √[∑in(xi - x̄)2 ∑i=1n(yi - ȳ)2]

式中：Rxy表示指标x和y的相关系数；xi表示第i年的冻融参数；yi表示第i年的驱动因子的值；x̄和ȳ分别表示冻融参数和驱动因子n年的平均值。Rxy > 0代表正相关，反之呈负相关，其绝对值越大代表冻融参数与驱动因子的关联程度越强。

2.3 数据来源

本文所使用的地表温度（0~7 cm土层土壤温度）数据，来自欧洲中期天气预报中心（European Centre for Medium-Range Weather Forecasts，ECMWF）提供的ERA5-LAND再分析数据集，空间分辨率达0.1°，时间分辨率为1 h，时间范围为2003年7月1日—2023年6月30日，通过逐小时数据逐像元获取研究区范围土壤的逐日最高最低温度数据。为了检验ERA5-LAND地表土壤温度数据在青藏高原地区适用性，本文采用青藏高原土壤温湿度逐时观测数据集[28-32]与ERA5-LAND再分析数据集进行对比，该数据集包括青藏高原4个典型区域，同样统计其逐日最高最低温度数据，对比结果如图3所示。ERA5-LAND地表温度在4个站点的R2分别为0.75、0.87、0.86、0.88，均值为0.84；均方根误差分别为7.82、3.46、4.05、6.27，均值为5.4 °C。考虑到ERA5-LAND地表温度为0~7 cm深度平均值，且与站点观测数据空间尺度的差异性，认为该数据满足青藏高原地区土壤冻融研究要求。其他辅助数据介绍见表2，所有数据均重采样到9 km，并统一投影至Albers_Conical_Equal_Area坐标系。

ERA5-LAND与站点地表温度散点图

图3 ERA5-LAND与站点地表温度散点图

表2 辅助数据

数据 指标 时间范围 空间分辨率 数据来源
气温 气温 2003年7月—2023年6月 1 km 国家青藏高原科学数据中心
降水 降水 2003年7月—2023年6月 1 km 国家青藏高原科学数据中心
NDVI NDVI 2003年7月—2023年6月 1 km MOD13A2 16d植被指数合成产品
地表反照率 地表反照率 2003年7月—2023年6月 500 m MCD43A4数据集
积雪深度 积雪深度 2003年7月—2023年6月 0.25° 国家青藏高原科学数据中心
太阳辐射 太阳辐射 2003年7月—2023年6月 4 km 全球每月气候和气候水平衡的网格数据集

3 结果与分析

3.1 青藏高原土壤冻融时空变化特征

3.1.1 青藏高原土壤冻融参数的空间分布特征

2003—2022年青藏高原土壤冻融参数的空间分布状况如图4所示，由图可知，青藏高原西北地区冻结起始时间较早，主要发生在7月，结束时间较晚，主要发生在6月，持续时间大部分区域在320 d以上，并且冻结天数多在280 d以上。而东南地区冻结起始时间较晚，主要集中在10月及月初，结束时间较早，主要发生在5月，持续时间及冻结天数大多区域不足200 d。青藏高原整体呈现出冻结起始时间由西北向东南延伸的趋势，西北低，东南高，越靠近东南区域值越高；冻结结束时间、持续时间及冻结天数则由东南向西北延伸的趋势，东南低，西北高，越靠近西北区域值越高。主要是因为西北地区海拔较高，多为连续多年冻土，土壤温度较低；东南地区海拔较低，多为季节冻土，土壤温度较高。在西藏山南市和林芝市部分区域的4种参数值均较小，其所在区域日最低土壤温度普遍高于0 °C，没有出现冻结现象，与冻土分布结果一致。

2003—2022年青藏高原土壤冻融参数的空间分布状况

图4 2003—2022年青藏高原土壤冻融参数的空间分布状况

为了更好地分析2003—2022年青藏高原空间上4种参数的变化情况，结合Theil-Sen斜率估计及M-K检验法对整个青藏高原进行整体趋势分析，结果如图5所示。从空间变化趋势来看，4种参数的变化趋势十分显著。冻结天数变化趋势最为显著，整个研究区绝大区域呈现减小趋势；冻结起始时间变化趋势最小，研究区绝大区域呈现不显著推迟趋势。近20年冻结起始时间显著推迟区域像元面积占比为29.27%，冻结结束时间显著提前区域像元面积占比为44.36%，冻结持续时间显著缩短区域像元面积占比为58.34%，冻结天数显著减小区域像元面积占比为74.05%。

2003—2022年青藏高原土壤冻融参数的变化趋势

图5 2003—2022年青藏高原土壤冻融参数的变化趋势

3.1.2 青藏高原土壤冻融参数的时间变化特征

为探究2003—2022年青藏高原土壤冻融整体的变化趋势，本文研究区剔除了湖泊和冰川区域，将其余部分在下文中称为青藏高原区。在此基础上，结合不同冻土类型将其进行空间平均处理，得到4种参数在不同区域的年际变化趋势（图6）。

2003—2022年青藏高原土壤冻融参数的年际变化情况

图6 2003—2022年青藏高原土壤冻融参数的年际变化情况

结果显示，2003—2022年，4种参数在不同区域随时间均存在明显变化趋势。不同冻土区土壤冻融参数的年际变化趋势与青藏高原区一致，但变化速率略有差异。多年冻土区冻结起始时间变化速率最快，为0.67 d·a-1，推迟12.7 d；季节冻土区变化速率为0.40 d·a-1，青藏高原区变化速率为0.51 d·a-1，分别推迟7.7 d和9.7 d。季节冻土区冻结结束时间下降速率最快，为-0.65 d·a-1，提前12.4 d；多年冻土区变化速率为-0.26 d·a-1，青藏高原区变化速率为-0.48 d·a-1，分别提前4.9 d和9.1 d。季节冻土区冻结持续时间下降速率最快，为-1.06 d·a-1，缩短20.1 d；多年冻土区变化速率为-0.92 d·a-1，青藏高原区变化速率为-0.99 d·a-1，分别缩短了17.5 d和18.8 d。多年冻土区冻结天数下降速率最快，为-1.08 d·a-1，缩短20.4 d；季节冻土区变化速率为-0.73 d·a-1，青藏高原区变化速率为-0.87 d·a-1，分别缩短13.9 d和16.5 d。综上所知，2003—2022年4种参数的明显变化代表地表温度呈显著上升趋势，这种趋势在青藏高原不同区域普遍存在。多年冻土区冻结起始时间变化速率最快，但冻结结束时间和冻结持续时间变化速率均低于季节冻土区。多年冻土区地表温度较低，热惯性较大，对气候变化响应滞后，而季节冻土区土壤热容量较小，对气候变化更为敏感，因此季节冻土区变化更为显著。但冻结天数的缩短代表多年冻土在逐渐退化。

3.2 青藏高原土壤冻融驱动因子定量分析

3.2.1 青藏高原冻结天数驱动力分析

为定量分析驱动因子对土壤冻融过程的影响，本文参考前人经验[34-36]，以冻结天数为例，选取海拔、气温、降水、NDVI、地表反照率、积雪深度和太阳辐射7个与地表温度相关的因子，探究不同因子在不同区域对冻结天数的驱动水平，结果如表3所示。

表3 2003—2022年不同区域冻结天数单因子探测

驱动因子 青藏高原区 多年冻土区 季节冻土区
海拔(X1) 0.474* 0.217* 0.316*
气温(X2) 0.605* 0.291* 0.428*
降水(X3) 0.490* 0.520* 0.449*
NDVI(X4) 0.604* 0.504* 0.589*
地表反照率(X5) 0.408* 0.170* 0.336*
积雪深度(X6) 0.191* 0.020* 0.102*
太阳辐射(X7) 0.182* 0.147* 0.175*

由表3可知，不同区域主导冻结天数的主导因子不同，7个因子在不同区域均通过95%的显著性检验，7个因子均对冻结天数的解释能力显著。气温是青藏高原区冻结天数空间分异的主导因子，其次依序为NDVI、降水、海拔、地表反照率、积雪深度、太阳辐射。降水是多年冻土区冻结天数空间分异的主导因子，其次是NDVI。而在季节冻土区NDVI对冻结天数的解释力最强。

通过交互探测（图7）可知，在青藏高原区，任意驱动因子对冻结天数空间分异的交互作用均为双因子增强，这表明任意2个驱动因子对冻结天数空间分异的交互作用都要大于单一因子作用。在多年冻土区，降水与气温对冻结天数的交互作用最显著（q=0.665），这反映了气候因子对多年冻土区冻结天数的主导调控作用。降水通过调节积雪深度间接影响地表热量传递过程，而气温通过改变降水的形式进一步加剧了降水对冻土热平衡的调节能力。积雪深度受到多种因子调控，未在交互作用中表现为强影响因子。在季节冻土区，NDVI与气温对冻结天数的交互作用最显著（q=0.762），表明植被覆盖状况对冻结天数的调节作用尤为突出。植被不仅通过改变地表热量交换方式影响土壤冻结过程，还与气温变化形成强烈的耦合关系。这种敏感性反映了季节冻土区生态系统对气候变暖的脆弱性，可能加剧冻融过程对农业和水资源的影响。

2003—2022年不同区域驱动因子对冻结天数的交互探测结果

图7 2003—2022年不同区域驱动因子对冻结天数的交互探测结果

3.2.2 青藏高原土壤冻融参数与气候的相关性分析

由表3可知，气候因子对冻结天数的交互作用十分显著。本文首先选取气温通过使用相关性分析法探究青藏高原土壤冻融参数与气温的相关性，如图8所示。图中阴影部分代表显著性≤0.05，即置信度水平达到95%以上，结果显著。由图8可知，冻结起始时间与年均气温的相关系数为-0.8~1.0，冻结结束时间、冻结持续时间、冻结天数与年均气温的相关系数为-1.0~0.8。其中，冻结起始时间与年均气温的相关系数为正值的像元面积为68.71%，31.29%区域呈现负相关关系，但并没有通过95%显著性检验；冻结结束时间与年均气温的相关系数为正值的像元面积为8.04%，通过95.00%显著性检验的像元面积为29.74%，主要分布在多年冻土区和季节冻土东南地区；冻结持续时间与年均气温的相关系数为正值的像元面积为15.30%，通过95%显著性检验的像元面积为21.38%，主要分布在季节冻土区；冻结天数与年均气温的相关系数为正值的像元面积为6.32%，通过95%显著性检验的像元面积为30.24%，主要分布在季节冻土区。

2003—2022年青藏高原土壤冻融参数与年均气温的相关系数分布

图8 2003—2022年青藏高原土壤冻融参数与年均气温的相关系数分布

为了更直观地体现土壤冻融参数与年均气温的相关性，去除显著性>0.05的区域，保留显著性≤0.05的区域（表4）。可以看出，青藏高原土壤冻融4种参数均与年均气温有较好的相关性：冻结起始时间与年均气温为正相关；冻结结束时间、冻结持续时间和冻结天数与年均气温为负相关。冻结天数与年均气温的相关性最大，为-0.55；冻结起始时间、冻结持续时间与年均气温的相关性最小，绝对值为0.48。通过对青藏高原土壤冻融参数整体统计，气温每升高1.0 °C，冻结起始时间推迟13.4 d，结束时间提前12.7 d，持续时间缩短26.0 d，冻结天数减少22.8 d。

表4 2003—2022年青藏高原土壤冻融参数与年均气温的平均相关系数

参数 平均相关系数
冻结起始时间 0.48
冻结结束时间 -0.52
冻结持续时间 -0.48
冻结天数 -0.55

通过土壤冻融参数与降水的相关性（图9）可知，4种参数通过95%显著性检验的像元面积分别为0.83%、0.59%、1.11%、0.75%，土壤冻融参数与降水之间并没有表现出很强的相关性。近些年来，青藏高原的降水量显著增加，但总量仍相对较少，一部分水分被植被吸收以及其它环境因素影响，可能使得降水对土壤冻融的影响被掩盖。

2003—2022年青藏高原土壤冻融参数与年均降水的相关系数分布

图9 2003—2022年青藏高原土壤冻融参数与年均降水的相关系数分布

3.2.3 青藏高原土壤冻融参数与NDVI和海拔的相关性分析

青藏高原土壤冻融参数与年均NDVI的相关性如图10所示。4种参数通过95%显著性检验的像元面积分别为1.25%、5.33%、4.89%、7.36%，这表明4种参数与NDVI之间仍没有显现出很强的相关性。青藏高原土壤冻融参数与NDVI之间的相关系数呈现明显的东西分带特征，通过显著性检验的像元多为东部地区，土壤冻融参数与NDVI的变化和气温保持一致。在高海拔地区冰川丰富，气温较低，植被覆盖度较低，冻结起始时间没有出现明显变化趋势；冻结结束时间与NDVI在青藏高原东部地区呈负相关，而在西部呈正相关，但并没有通过显著性检验，植被变化在气候比较湿润的地区对冻结结束时间的影响更大；冻结持续时间和冻结天数的变化是前两者协同作用的结果。

2003—2022年青藏高原土壤冻融参数与年均NDVI的相关系数分布

图10 2003—2022年青藏高原土壤冻融参数与年均NDVI的相关系数分布

不同海拔之间土壤冻融参数的空间分布不同，本文通过绘制箱线图（图11）分析不同海拔之间土壤冻融参数的差异。由图可知，海拔在2 000 m以下4种参数的时间尺度跨度较大，其主要分布在西藏自治区东南端、新疆维吾尔自治区部分以及青藏高原东南端部分区域。经纬度跨度较大，冻结时间也有所不同。在2 000 m以上，随海拔的升高，冻结起始时间逐渐提前，当海拔＞6 000 m时，起始时间达到最低点，高海拔地区温度更低，在冻融阈值以下的地区较多。其余3种参数分布基本一致，随海拔升高，冻结结束时间推迟，持续时间和冻结天数延长。图中黑点表示异常值，不同参数在不同海拔下均出现异常值，这表明4种参数的变化不单纯取决于海拔高度，其波动可能是由气候、地形以及其它因子等综合影响的结果。

2003—2022年青藏高原土壤冻融参数与海拔的箱线图

图11 2003—2022年青藏高原土壤冻融参数与海拔的箱线图

4 讨论

由上述研究可知,青藏高原2003—2022年冻融过程呈现出冻结起始时间推迟，结束时间提前，持续时间缩短，冻结天数减少的趋势。并且冻融过程存在明显的空间异质性，青藏高原西北地区冻结起始时间较早，结束时间较晚，持续时间和冻结天数较长，东南地区与此相反。这与海拔有关。青藏高原土壤冻融过程在不同区域存在明显的变化趋势，但变化速率略有差异。多年冻土区冻结起始时间变化速率最快，已有研究表明，受气候变化影响，低温冻土区地表温度升温速率高于高温冻土区[37]。多年冻土区地表温度较低，接近冻结点，随气候变化更易冻结。但结束时间和持续时间变化速率低于季节冻土区，这是由于季节冻土区海拔较低，土壤热容量较小，对气候变暖响应更快，因此季节冻土区冻土更易融化。冻结持续时间受起始时间和结束时间共同作用影响，因此多年冻土区变化速率较低。从青藏高原整体来看，青藏高原区冻结起始时间变化速率为0.51 d·a-1，冻结结束时间变化速率为-0.48 d·a-1，冻结持续时间变化速率为-0.99 d·a-1，冻结天数变化速率为-0.87 d·a-1，这与已有研究变化趋势一致[38-39]，但数值不一致。如杨淑华等[21]认为冻结起始时间变化速率为0.72 d·a-1，冻结结束时间变化速率为-0.40 d·a-1，冻结持续时间变化速率为-1.13 d·a-1，冻结天数变化速率为-0.93 d·a-1。数值不一致可能：① 选取青藏高原研究区范围不同，本研究中青藏高原范围剔除了湖泊和冰川的影响；② 选取的数据源有所不同。本研究是基于ERA5-LAND地表温度为基础获取的，而杨淑华等[21]采用87个气象台站资料获取，数据精确但无法全面反映整个青藏高原的冻融状况；③ 研究的时间范围也有所不同，本文以2003—2022年展开研究，而杨淑华等[21]研究范围为1980—2015年，变化速率与所选时间尺度有较大关系，不同时间尺度下对同一事物的变化速率计算有一定差异。

以往研究仅探讨单因子对土壤冻融过程的影响[39]，本文综合采用Theil-Sen斜率估计法、Mann-Kendall检验法、地理探测器模型以及相关性分析等方法，有效避免了单一方法可能导致的局限性，弥补了土壤冻融过程定量化研究不足，实现更全面的冻融参数影响因素分析。但仍有不足之处：

（1）使用地理探测器分析将多种数据统一重采样至9 km，通过自然间断法重分类获取q值，结果显著，可以反映不同驱动因子对冻结天数的驱动力水平[40]。但在模型中，数据离散化和空间尺度效应也带来一定的影响差异，今后在研究中可以将多种数据离散化方法结合空间尺度参数寻找最优解，以实现更为准确的分析。

（2）本文以冻结天数为例，结合地理探测器分析青藏高原不同区域冻结天数空间分异的主导因素，发现气温是青藏高原区冻结天数空间分异的主导因子。近地表土层是对气温最为敏感的界面，地表温度对气温的变化更为敏感。使用相关性分析法探讨了4种冻融参数与气温的相关性，结果显著，气温与冻结起始时间呈正相关，与结束时间、持续时间和冻结天数呈负相关。降水和NDVI对4种冻融参数没有表现出明显的相关性，而不同海拔下4种冻融参数呈现明显梯度差异，可见冻融过程是一个复杂的多因子驱动系统，理解这些因子如何共同作用并影响冻融过程十分重要。但本文使用地理探测器只选取了部分影响因子，想要深入探讨局部影响因素，还需要加入更多影响因子如土壤湿度等参数，结合地理探测器等方法综合考虑。

5 结论

本文使用ERA5-LAND再分析数据分析了2003—2022年青藏高原土壤冻融参数的时空变化特征并结合地理探测器分析冻结天数空间分异的主导因素，最后使用相关性分析法分析主要驱动因子的相关性，结论如下。

（1）青藏高原2003—2022年西北地区冻结起始时间较早，结束时间较晚，持续时间长，冻结天数多；而东南地区与此相反。青藏高原土壤冻融参数空间变化趋势显著，其中冻结起始时间显著推迟区域的像元面积占比为29.27%，冻结结束时间显著提前区域像元面积占比为44.36%，冻结持续时间显著缩短区域像元面积占比为58.34%，冻结天数显著减小区域像元面积占比为74.05%。

（2）不同区域随时间均存在明显变化趋势，多年冻土区冻结起始时间变化速率最快，为0.67 d·a-1，推迟12.7 d；季节冻土区冻结持续时间和冻结天数变化速率最快，分别为-0.65 d·a-1和-1.06 d·a-1，分别推迟和缩短12.4 d和20.1 d；多年冻土区冻结天数变化速率最快，为-1.08 d·a-1，缩短20.4 d。冻结天数的缩短代表多年冻土在逐渐退化。

（3）使用地理探测器分析冻结天数空间分异的主导因素，对青藏高原区的解释力依序为：气温（0.605）＞NDVI（0.604）＞降水（0.490）＞海拔（0.474）＞地表反照率（0.408）＞积雪深度（0.191）＞太阳辐射（0.182）。通过交互探测结果发现，气温与其它因子的交互作用在不同区域的解释力均最强。气温直接影响地表热平衡，控制土壤冻融过程。在多年冻土区降水通过调节积雪深度间接影响地表热量传递过程，对调节冻结天数发挥重要作用。在季节冻土区，植被覆盖状况对冻结天数的影响尤为突出，季节冻土区土壤较薄，热容量较小，对气候变化更为敏感，更容易受到植被覆盖的影响。

（4）结合相关性分析法分析土壤冻融参数与前4个因子的相关性发现，土壤冻融参数与气温变化更为敏感，与降水和NDVI之间没有显现出明显的相关性。受多种驱动因子综合影响，降水和NDVI对土壤冻融的影响被掩盖。土壤冻融参数与海拔表现出明显的梯度差异，随海拔升高，冻结起始时间提前，结束时间推迟，持续时间和冻结天数延长。