由于特殊的地理单元和较低的年平均气温，青藏高原约40%的区域赋存多年冻土^[1-3]。受全球暖湿化和人类工程活动的共同影响，青藏高原多年冻土区正处于持续的退化状态，主要表现为地温升高、活动层厚度增加和多年冻土层融化等^[4-8]。青藏高原作为出疆入藏大通道的必经之路，修筑有青藏铁路、青藏公路和共玉高速公路等重大交通基础设施。由于青藏高原暖湿化导致多年冻土持续退化，很多路段出现了较为严重的沉陷病害。

蔡汉成等^[9]分别通过现场监测和理论分析等手段研究了青藏铁路、青藏公路和G214共玉高速公路等设施下伏多年冻土的演变规律，其中1999—2014年内青藏铁路沿线多年冻土平均人为上限与多年冻土平均天然上限退化速度分别为4.5 cm/a和1.9 cm/a；2002—2022年青藏公路多年冻土平均人为上限与多年冻土平均天然上限退化速度分别为15.5 cm/a和4.9 cm/a；2003—2014年G214共玉高速公路，左路肩、路基中心及右路肩下多年冻土人为上限下降速率分别为16.82，25.36，16.73 cm/a^[9-12]。由此可见，青藏高原多年冻土退化已成为不争的事实，但目前鲜有针对青藏高原高速公路路基下伏多年冻土时空演变特征及其影响因素的研究。

综上所述，本研究选取共玉高速公路巴颜喀拉山至清水河段路基下伏多年冻土为研究对象，利用ArcGIS对研究区内公开的气候、地形及钻芯数据进行深入挖掘。分析了高程、坡度、坡向这3个地形地貌影响因子和归一化植被指数(NDVI)、年平均气温(MAT)、年平均地温(MAGT)、平均年降雨量(MAP)和年平均湿度(MAH)这5个时变影响因子作用下高速公路路基下伏多年冻土时空演变特征，并与青藏铁路和青藏公路分析结果进行对比，以期为青藏高原冻土区高速公路的治理改造和选线提供支撑和建议。

研究区位于青藏高原不稳定型多年冻土-非冻土过渡带，沿路线上行方向，地形以山前斜坡和山前平原为主，植被较为发育，冻土类型为少-饱冰冻土，不稳定型多年冻土分布区与非冻土区泾渭分明，不稳定型多年冻土共41.5 km，占路线总长的83%，非冻土区共8.5 km，占路线总长的17%；研究区高程介于4 400~4 700 m之间，沿路线上行方向高程呈现逐渐降低的趋势；沿路线方向多年冻土基础影响因子在不稳定型多年冻土分布区与非冻土区内栅格数据差异性大，规律性强；时变影响因子栅格数据逐年变化显著。

2011年研究区路线K610+000~K611+936.5多年冻土人为上限属于0~1.5 m区间范围，平均上限为1.37 m；路线K611+936.5~K637+600多年冻土人为上限属于2~2.5 m区间范围，平均上限为2.31 m；路线K637+600~K642+100多年冻土人为上限属于2.5~3 m区间范围，平均上限为2.56 m；路线K642+100~K650+000无冻土分布。即沿路线上行方向多年冻土上限逐渐下移，直至无冻土分布。

研究区的基础数据源包括：

(1) 2011年统计汇总的共玉高速公路详细工程地质勘察报告包括钻孔日期、位置和深度信息、工程地质描述和土质分层数据；

(2) 2019年和2022年共玉高速公路病害路段冻土工程地质条件及工程措施建议一览表包括钻孔日期、位置和深度信息及土质分层数据；

(3) 数字高程模型(Digital Elevation Model，DEM)，来自于地理数据云免费公开的基础数据，主要用于高程、坡度和坡向等地形地貌基础环境因子的获取，DEM栅格数据的空间分辨率为30 m；

(4) NDVI，MAT，MAGT，MAP，MAH数据来源于资源环境科学数据注册与出版系统，气候因子和NDVI栅格数据的空间分辨率均为1 km，在研究区范围内栅格数据精度可达到90%以上^[13]；

(5) 多年冻土稳定型数据来源于2005—2015年青藏高原多年冻土稳定性制图^[14]。

利用ArcGIS中“面转栅格”工具将59处多年冻土发育钻孔面和研究区整体分别划分为157 470和195 106个栅格单元，栅格单元中储存有地理信息数据。基础环境因子的有效选取会影响评价结果的可靠性与准确性。其选取原则主要包括环境因子的客观存在性、显著性和继承性等^[15-16]。在这些原则的基础上本研究选取了高程，坡度，坡向，NDVI，MAT，MAGT，MAP，MAH共8类多年冻土基础影响因子。本研究主要运用自然间断点法把选取的高程，坡度，坡向和NDVI因子划分不同属性区间，运用等间隔法把选取的MAT，MAGT，MAP，MAH划分不同属性区间，并得到各基础影响因子的属性区间分级及其频率比值(见表 1)。频率比值的大小决定了该环境因子属性区间是否有利于多年冻土发育，比值>1有利于多年冻土发育，比值＜1则不利于多年冻土发育^[17]。

地形地貌是多年冻土形成和发育的重要影响因素之一，它与气候、植被和水文地质等因素相互作用，对多年冻土的形成、分布和性质产生显著影响。本节综合高程、坡度和坡向共3个地形地貌因子进行频率比和图谱分析。结果表明，当高程属于4 561~5 089 m区间范围，坡度属于0~6.74°区间范围，坡向为平面、东北(阴坡)、北(阴坡)、东(阴坡)和东南(阴坡)时，频率比值>1，有利于多年冻土形成和发育；当高程属于4 440~4 561 m区间范围，坡度属于6.74~14.79°区间范围，坡向为南(阳坡)、西北(阳坡)、西南(阳坡)和西(阳坡)时，频率比值＜1，不利于多年冻土形成与发育。这说明高程与多年冻土形成和发育程度呈正相关，坡度与多年冻土形成和发育程度呈负相关，阴坡较阳坡更有利于多年冻土的形成与发育(见表 1)。

植被覆盖率对多年冻土的分布与发育会产生较为明显的影响。一般情况下，茂密植被覆盖下近地表层冻土的地温较低，温度场平面分布较为连续^[18]，本节对NDVI因子进行频率比和图谱分析。结果表明，当NDVI属于0.6~0.76区间范围，频率比值>1，有利于多年冻土形成和发育；当NDVI属于0.38~0.6区间范围，频率比值＜1，不利于多年冻土形成与发育。这说明NDVI与多年冻土形成和发育程度呈正相关(见表 1)。

多年冻土是气候和环境因素的共同产物，对温、湿度和降雨非常敏感，气候条件决定着多年冻土要素的生存寿命。为揭示青藏高原多年冻土区南界冻土空间分布差异性的内在驱动力，综合MAT，MAGT，MAP，MAH共4个气候因子进行频率比和图谱分析。结果表明，当MAT属于–7~–5 ℃区间范围，MAGT属于–3~–1.5 ℃区间范围，MAP属于450~470 mm区间范围，MAH属于61%~62%区间范围时，频率比值>1，有利于多年冻土形成和发育；当MAT属于–5~–4.5 ℃区间范围，MAGT属于–1.5~–0.5 ℃区间范围，MAP属于470~500 mm区间范围，MAH属于62%~64%区间范围时，频率比值＜1，不利于多年冻土形成与发育。这说明MAT，MAGT，MAP，MAH与多年冻土形成和发育程度呈负相关。降雨会增加土壤含水量、提升土层导热强度、增大地表热传递和土体融化潜热量，从而导致活动层和多年冻土温度升高^[19]；年空气平均湿度越高说明地表水更多地以蒸发的形式耗散，水分蒸发耗热量大，土层中的热量散失多，热通量向下，不利于下层地基土多年冻土的形成与发育^[20](见表 1)。

本研究定义了8个通用参数指标，分别为影响因子年均幅值波动速度v_i：影响因子平均年均幅值波动速度v_i：路基下伏多年冻土退化速度v_d：多年冻土人为上限退化量d_r：多年冻土天然上限退化量d_t：多年冻土人为上限的下移速度v_r：多年冻土天然上限的下移速度v_t：新建路基路面工程后多年冻土上限下移速度增量v_e。在本研究中共有2个统计周期，分别为2011—2019周期和2019—2022周期，这2个统计周期内分别采用不同的参数指标进行计算评估，n为统计周期时长，m为钻孔数量。

2011—2019周期统计结果为：

(1)

(2)

(3)

2019——2022周期统计结果为：

(4)

(5)

(6)

(7)

(8)

多年冻土人为上限h_r指的是人为活动影响导致多年冻土天然上限的变化而形成的新的上限，在本研究中是人为新建路基路面结构导致多年冻土天然上限变化而形成的上限，具体为在共玉高速公路路面结构上进行钻探，钻芯结果中多年冻土上限距原天然大地表层的垂直距离；多年冻土天然上限h_t指的由自然原因形成的上限，具体为在天然大地上进行钻探，钻芯结果中多年冻土上限距天然大地上表层的垂直距离。比较不同年份多年冻土上限数值大小，可分析多年冻土上限随时间的变化规律，将不同年份多年冻土上限做差/比值，可得到v_d，v_r，v_t，v_e等数据。

采用夏皮洛-威尔克检验，验证各要素的正态分布显著性，结果如表 2所示，所有要素的显著性指数α均大于0.05，符合正态分布，可进行Pearson线性相关分析。

采用Pearson函数分别求解要素间相关系数，结果如图 1(a)所示，同时采取配对要素t检验方法，对相关系数显著性进行分析，结果如表 3所示。将地形影响因子中的坡向数据分为阴坡和阳坡2个类别，采用哑变量设置后进行线性回归分析，结果如表 4所示。

图 1 多年冻土影响因子间相关系数 Fig. 1 Correlation coefficients among influencing factors of permafrost

图选项

分析要素间的相关系数，得到各时变影响因子年均幅值波动速度v_i/各地形地貌影响因子数值与v_d间的相关性(见图 1)；分析线性相关的显著性参数，得到坡向与v_d间的相关性(见表 4)。依据v_d数值大小将研究区分为多年冻土低速退化区(v_d＜30 cm/a)与多年冻土高速退化区(v_d≥30 cm/a)，将多年冻土低速退化区与多年冻土低速退化区v_d分别与5个时变影响因子年均幅值波动速度v_i和2个地形地貌影响因子数值(除坡向外)之间进行相关性分析，结果如图 1(b)和图 1(c)所示。

比较不同年份多年冻土上限数值，可分析多年冻土上限随时间的变化规律，将不同年份多年冻土上限分别做差值和比值分析，可得到v_d，v_r，v_t等数据，结果如图 2所示。

图 2 2011，2019，2022年多年冻土上限对比 Fig. 2 Comparison of permafrost upper limits in 2011, 2019 and 2022

图选项

将2011年和2019年多年冻土低速退化区与多年冻土低速退化区影响因子平均值分别做差值和比值分析，研究v_d的空间异质性，结果如图 3所示。

图 3 2011年与2019年多年冻土高、低速退化区影响因子对比 Fig. 3 Comparison of influencing factors in high and low permafrost degradation zones in 2011 and 2019

图选项

分析2019—2022年路基路面结构对多年冻土上限的退化影响，比较v_e与v_t的数值，分析路基路面结构与自然环境因素对多年冻土上限的退化影响程度大小，结果如图 4所示。

图 4 多年冻土上限下移速度对比 Fig. 4 Comparison of descent speeds of permafrost upper limit

图选项

当相关系数R越接近1时，对样本数据的正线性相关程度越强；当相关系数R越接近-1时，对样本数据的负线性相关程度越强；当相关系数R越接近0时，对样本数据的线性相关程度越弱^[21]。由表 3可知，所有配对要素的显著性指数p均小于0.05，则可以拒绝零假设，相关系数是显著的。由图 1(a)可知，研究区内时变影响因子中R[v_i(MAGT), v_d]=0.835，相关性最强；R[v_i(NDVI), v_d]=–0.466，相关性最弱；地形影响因子中R[H, v_d]=–0.695，相关性最强；R[S, v_d]=–0.3，相关性最弱；由图 1(b)可知，多年冻土低速退化区中R[v_i(MAH), v_d]=–0.989，相关性最强，R[v_i(NDVI), v_d]=–0.574，相关性最弱，地形影响因子中R[H, v_d]=–0.811，相关性最强，R[S, v_d]=–0.728，相关性最弱；由图 1(c)可知，多年冻土高速退化区中：R[v_i(MAGT), v_d]=0.802，相关性最强，R[v_i(NDVI), v_d]=–0.039，相关性最弱，地形影响因子中R[H, v_d]=–0.742，相关性最强，R[S, v_d]=–0.035，相关性最弱。阳坡的冻土退化程度显著高于阴坡(p＜0.05)，多年冻土退化速度比阴坡高174.6%(见表 4)。

多年冻土天然上限是评价自然环境下多年冻土稳定性的重要指标。2019—2022年天然大地下伏多年冻土天然上限变化情况见图 2(b)。2019年和2022年对应位置钻孔多年冻土天然平均上限值分别为308 cm和338 cm，这说明2022年天然大地下伏多年冻土上限较2019年出现明显的下移，为2019年天然大地下伏多年冻土上限的1.1倍。由于近年来青藏高原暖湿化进程不断推进，气温和降雨量逐年递增，导致天然大地下伏多年冻土处于持续的退化状态。

多年冻土人为上限是评价多年冻土区高速公路路基稳定性的重要指标。2011—2019和2019—2022年高速公路路基下伏多年冻土人为上限变化情况见图 2。2011年和2019年对应位置钻孔多年冻土人为平均上限值分别为244 cm和532 c m，这说明2019年路基下伏多年冻土上限较2011年出现明显的下移，为2011年多年冻土人为上限的2.18倍。2019年和2022年对应位置钻孔多年冻土人为平均上限值分别为436 cm和502 cm，这说明2022年路基下伏多年冻土上限较2019年出现明显的下移，为2019年路基下伏多年冻土上限的1.15倍。对于高速公路路基下伏多年冻土而言，由于外部环境的逐年变化，以及路基路面结构的吸热、导热作用，导致路基下伏多年冻土的升温或融化，因此，2011—2022年路基下伏多年冻土人为上限呈现出明显的下降趋势。

由图 2(a)可知，2011—2019年研究区多年冻土退化速度呈现出一定的空间异质性。按路基下伏多年冻土退化速度v_d数值大小可将研究区划分为2类，多年冻土低速退化区(v_d＜30 cm/a)和多年冻土高速退化区(v_d≥30 cm/a)，低速退化区起讫点分别为K611+650和K627+100，高速退化区起讫点分别为K627+100和K635+500。低速退化区和高速退化区平均退化速度分别为23.96 cm/a和45.16 cm/a，高速退化区多年冻土平均退化速度为低速退化区的1.88倍。

由图 3(a)和图 3(b)可知，2011年多年冻土高速退化区时变影响因子NDVI，MAT，MAGT，MAH，MAP的平均值分别是多年冻土低速退化区的0.969，0.942，0.823，1.001，0.993倍，其中MAGT平均值相差比例最大，高速退化区比低速退化区高出了0.317 ℃；2019年多年冻土高速退化区时变影响因子NDVI，MAT，MAGT，MAH，MAP的平均值分别是多年冻土低速退化区的1.046，0.937，1.321，0.997，1.003倍，其中MAGT平均值相差比例最大，高速退化区比低速退化区高出了0.353 ℃。2011—2019年研究区多年冻土退化速度呈现出的空间异质性，主要取决于高速退化区和低速退化区内MAGT的空间差异性。

多年冻土人为上限下移速度是反映多年冻土路基沉降变形的重要指标。2019—2022年多年冻土人为上限平均下移速度值、多年冻土天然上限平均下移速度值和新建路基路面结构后多年冻土上限平均下移速度增量值分别为22.0，9.9，12.2 cm/a。故在时间维度上说明2019—2022年路基路面结构加速了多年冻土上限的下移速度，高速公路路基下伏多年冻土上限为天然大地下伏多年冻土上限下移速度的2.22倍，路基路面结构导致的多年冻土上限下移速度为自然环境因素的1.23倍。若研究区高速公路按照现多年冻土退化速率持续退化，那么路基下伏多年冻土的融化将可能导致高速公路产生更大的沉降变形。针对研究区高速公路特殊路基结构中的片块石路基和高填方路基进行研究，2019—2022年片块石路基和高填方路基新建路基路面结构后多年冻土上限平均下移速度增量值分别为11.3 cm/a和12.7 cm/a。这说明高填方路基对多年冻土上限下移速度影响更大，约为片块石路基的1.12倍(见图 4)。

本研究通过对公开的气候、地形及钻芯数据进行深入挖掘，研究青藏高原高速公路路基下伏多年冻土的时空演变特征，得到主要结论如下。

(1) 2011—2019年研究区高速公路路基下伏多年冻土平均人为上限下移量为2.88 m，退化速度为36 cm/a，多年冻土人为上限退化速度分别为同期青藏高原多年冻土核心区青藏铁路和青藏公路的8倍和2.32倍，为同期同区位G214共玉高速公路的1.83倍。

(2) 2011—2019年沿线多年冻土退化呈现出一定的空间异质性。路线K611+650~K627+100为低速退化区，K627+100~K635+500为高速退化区，影响高速退化区多年冻土退化速度的主要因素是MAGT，影响程度为0.802；影响低速退化区多年冻土退化速度的主要因素是MAH，影响程度为0.989；影响整个研究区多年冻土退化速度的主要因素是MAGT，影响程度为0.835。

(3) 2019—2022年多年冻土人为上限v_r、多年冻土天然上限v_t平均下移速度和新建路基路面结构后多年冻土上限平均下移速度增量v_e分别为0.220，0.099，0.122 m/a，故在时间维度上说明路基路面结构加速了多年冻土上限的下移速度，且高速公路路基下伏多年冻土上限下移主要是由路基路面结构修筑引起的，影响程度为自然环境因素的1.23倍。

(4) 在同类型工程的选线设计中，优先选择无多年冻土分布路线，若天然大地赋存多年冻土，则优先选择高海拔阴坡路线，后依次选择低海拔阴坡、高海拔阳坡和低海拔阳坡路线；若路线走向定于后3种情形，需对路基采取延缓多年冻土退化的处置措施；同时选线区域宜具备MAGT，MAT，MAP年均幅值波动速度小，MAH年均幅值波动速度大的条件。在同类型工程的建设与养护中，宜采用高功率强制冷措施降低地温。