基于GWRRK的多年冻土地温空间分布模拟研究——以青藏铁路昆仑山至尺曲谷地段为例

姚明星，赵 锐，杨林川，齐 华，廖 昕，孟祥连，周福军

（1.西南交通大学地球科学与环境工程学院，四川 成都 611756；2.西南交通大学建筑学院，四川 成都 611756；3.中铁第一勘察设计院集团有限公司，陕西 西安 710043）

0 引言

随着全球气候变化加剧，多年冻土面临退化趋势，可能对寒区的生态安全产生影响，诱发诸如崩塌、滑坡、泥石流等地质灾害问题［1-3］。多年冻土对地温变化极其敏感［4］，通过研究其地温空间分布，可探索冻土活动层厚度的变化规律，分析冻土热融状况以及未来变化趋势，为冻土灾害防治提供科学依据。

多年冻土地温既受海拔、纬度和经度的影响，也受气象、地形地貌、植被、积雪、土壤岩性、含水量等局地因素的影响［5-10］。学者基于上述影响因素，结合冻融指数、热传导、地理加权回归、支持向量机、机器学习等方法开展了多年冻土地温的时空分布模拟研究［11-16］。例如，Jafarov等［11］、Nicolsky等［17］利用GIPL（geophysical institute permafrost laboratory）模型，考虑气温、降水、土壤热特性和土壤含水量的影响，基于含相变的一维非线性热传导方程，模拟了阿拉斯加州多年冻土的季节冻融深度。Liu等［18］利用二维非稳态模型分析了青藏高原多年冻土区路基温度，模拟了冻土土壤与大气的传热过程。Luo等［19］综合考虑了气温、植被、积雪深度和土壤属性等因素，对黄河源区多年冻土地温的空间分布进行了模拟。Qin等［20］、Yin等［21］利用GIPL模型分别对青藏高原多年冻土活动层厚度和青藏高原工程走廊年平均温度进行模拟，前者基于土壤导热系数开展多年冻土活动层厚度模拟，后者侧重于土壤性质对多年冻土地温的影响。Liu等［22］综合考虑气温、风速、风向、地表水分蒸发量、太阳辐射等影响因素，基于二维非稳态热传导有限元模型，揭示了冻土路基在给定工况下的地温时空分布特征。Zhao等［15］利用地理加权回归模型对青藏铁路沿线多年冻土地温分布进行模拟，同时指出各因素对地温的影响随时空尺度变化增大，呈现出典型的空间非平稳性特征。此外，各影响因素间可能还存在着交互影响，易引起多重共线性问题［15］。

本研究利用地理加权岭回归克里金方法（geographically weighted ridge regression Kriging，GWRRK）构建冻土地温分布模型，以青藏铁路昆仑山至尺曲谷地段多年冻土覆盖区域为例，验证模型的有效性，通过对该区域多年冻土地温分布的揭示，分析冻土融化深度的变化特征。地理加权岭回归克里金方法（GWRRK）既可以充分考虑多年冻土地温与其影响因素间的空间非平稳性，也可消除影响因素之间的多重共线性［23］，同时还考虑了拟合残差的影响。昆仑山至尺曲谷地段是青藏铁路建设和运营中冻土变化监测最密集的区域，通过揭示该区域冻土地温空间分布变化趋势，可为铁路安全运营及沿线地质灾害防治提供科学依据。

1 研究区与数据

1.1 研究区概况

研究区位于青藏高原中部，隶属于青海省玉树藏族自治州，北至昆仑山越岭区，南到尺曲谷地，横跨楚玛尔河高平原、可可西里山区、北麓河盆地、风火山山区，该地区是典型的高海拔多年冻土区域。研究区平均海拔4 900 m，区内山系呈东西走向，山系与河流交替出现，形成高平原、盆地和谷地等地貌。河流主要有楚玛尔河、北麓河和通天河，由冰川融水和大气降水补给，流量季节性变化较大。研究区具有独特的冰缘干旱气候特征，且随海拔增高有明显的气候垂直分带性特征。

1.2 数据来源与处理

本文以地温钻孔实测数据为基础，从地理要素、气象条件和局地因素三个维度选取地温影响因素，构建冻土地温分布模型。其中，地理要素主要考虑纬度和高程［16，24-25］；气象条件选择近地表气温、风速、降水率、地面向下长波辐射和地面向下短波辐射［26-29］；局地因素则考虑地层岩性、坡度、坡向、植被覆盖度［30-33］。

青藏高原气温往往在7月达到最大，且由于向下传递时存在一定滞后性［34-35］，本文重点对研究区7—9月的地温空间分布开展模拟，分析其对应的冻土活动层厚度的空间变化特征。本文实测地温数据是来自2001年修建青藏铁路工程的钻孔地温数据（采用热敏电阻测温，测温深度为0.2～15 m，数据采集间隔为10天至一年不等）。在研究区域内所获取的钻孔点位一共237个，7月至9月的钻孔地温数据点位分布较为均匀，共307次测温数据。其中，7月测温96次，8月测温128次，9月测温83次。通过插值，获取时间序列数据进行对比分析，可更好了解7月至9月研究区内整个路段多年冻土地温变化及冻土融化深度变化情况。钻孔位置从北向南穿越不同地形区域，基本反映了青藏高原中低纬度、高海拔的高原冻土特征，样本点具有较好的代表性。研究区及钻孔测温点位分布如图1所示，其中左上角小图中的多年冻土分布矢量图来自Ran等［36］的研究成果。

图1 研究区及钻孔点位分布Fig.1 Study area and distribution of the boreholes

由于研究区内仅有五道梁一个气象观测站，难以满足模拟的空间分辨率要求。因此，本文气象数据采用青藏高原科学数据中心公开发布的中国区域地面气象要素驱动数据集（1979—2018年），涉及近地面气温、地面降水率、近地面全风速、地面向下短波辐射和地面向下长波辐射共5个指标［37-38］。上述数据水平空间分辨率为0.1°。本文利用ArcGIS 10.2软件对上述数据进行栅格转点矢量文件处理，并综合对比克里金法、反距离权重法、自然邻域法和趋势面法所获取的1 km×1 km分辨率气象数据结果，选用拟合精度最优的克里金插值结果作为气象基础数据。地层岩性数据通过全国2002年发布的1∶100万地质图裁剪获取。地层岩性因素是绘制不同地层岩性区域的实测地温数据与其深度的拟合曲线，计算曲线的拐点坐标，利用拐点的深度坐标值作为地层岩性量化因子。坡向和坡度数据是利用ArcGIS 10.2软件对研究区DEM（分辨率30 m）计算。植被指数来源于Terra卫星MOD13 Q1全球250 m分辨率植被指数16天的合成产品，并利用像元二分模型计算获取植被覆盖度［39］。

2 模型构建

2.1 地理加权岭回归模型

地理加权回归模型（GWR）是在普通线性回归模型的基础上，将空间位置信息嵌入到回归中，以预测每个位置对应因变量的变化，表达式［40-41］为

式中：(ui，vi)为第i个采样点的坐标（如经纬度）；p为影响因素的个数，本文共考虑了11个影响因子，故p=11；βk(ui，vi)为第i个采样点上的第k个回归参 数，是 关 于 地 理 位 置 的 函 数，εi～N(0，σ2)，Cov(εi，εj)=0(i≠j)。

为解决地理加权回归中可能存在的因子间多重共线性问题，Wheeler［42］将地理加权回归（GWR）与岭回归结合，提出了地理加权岭回归模型（GWRR），其回归参数可估计为

式中：λ为岭回归参数。本文的岭回归参数λ通过设置方差膨胀系数阈值小于5（VIF＜5）确定得到。

在地理加权岭回归模型中，第i个点上的因变量拟合值［42］为

2.2 地理加权岭回归克里金模型

地理加权岭回归克里金模型（GWRRK）是将GWRR与普通克里金插值方法（ordinary Kriging，OK）结合，即在GWRR拟合后，利用OK处理回归中的残差问题，以此来提高GWRR的预测精度，表达式［23］为

第i处地理加权岭回归模型（GWRR）的模拟值表述［23］为

样本点数值xi(i=1，2，…，n)与插值点的空间区域变化函数为Z(xi)，则GWRR第i处局部回归残差项的克里金插值模拟值可表达［23］为

上式必须满足以下条件，

式中：C(xi，xj)为采样点xi与xj之间的变异函数值；μ为拉格朗日常数。

2.3 具有外部漂移的克里金模型

具有外部漂移的克里金模型（Kriging with an external drift，KED）是带有趋势克里金模型的拓展，既考虑了样本点与待估计点之间的空间分布几何特征，又考虑了样本间空间关系的结构特征［43］。

设u处属性值Z(u)由漂移函数m(u)和残差R(u)组成［43］，即

漂移函数通常可表示为

式中：β0、βl为未知系数；xl(u)为漂移函数的辅助变量。

样本点数值xi(i=1，2，…，n)与插值点的空间区域变化函数为Z(xi)，则残差R(u)可表达［23］为

上式必须满足式（7）。

3 结果与讨论

3.1 模型精度验证

本文采用留一交叉验证方法（leave-one-out cross-validation，LOO-CV）对模型结果进行验证，其基本思想是每次从个数为N的样本集中抽取出一个样本作为验证集，剩下的N-1个样本作为训练集，重复进行N次，依次取遍所有N个数据作为验证集，最后将平均N个数据的结果作为泛化误差的估计［44］。该方法较其他验证方法具有训练样本量大、接近原始样本分布的特点［45］。

图2是研究区2001年7月至9月钻孔点位在0.2～15 m深度上的留一交叉验证结果。可以看出，模型具有较高的拟合优度，平均拟合优度达0.89。其中，7月的拟合优度最大，R2=0.91，且平均绝对误差（mean absolute error，MAE）和均方根误差（root mean square error，RMSE）相对较小；8月的拟合优度最小，R2=0.87；9月拟合结果的均方根误差（RMSE）值最小，达0.694。7月至9月精度验证结果的平均绝对误差整体偏小，最大误差为0.613℃，说明本文所建GWRRK模型具有一定的稳健性，可用于案例研究区域的多年冻土地温分布模拟分析。此外发现，模拟地温与实测地温在-4～2.5℃的温度区间范围内重合度较高，说明在该温度区间内模拟的误差相对较小，而在大于2.5℃的温度区间内偏差较大，可能由于不同地形区域内地温波动较大所致。

图2 GWRRK模型留一交叉验证结果Fig.2 Scatter distribution of LOO-CV results regarding the GWRRK model：in July（a），August（b）and September（c）

为进一步验证模拟结果的有效性，本文分别从不同地形区域中随机选取1个钻孔测温点，点位分布如图3所示，将其对应的模拟结果与实测结果对比。由图4可知，在不同区域的地温模拟值与实测值的拟合精度R2都大于0.95；均方根误差（RMSE）最小达0.149。此外，除图4（a）昆仑山区域7月的模拟值与实测值有一定偏差外，其他位置不同时间段的模拟结果与实测结果的变化趋势几乎一致。这也说明GWRRK在该研究区域尺度内的模拟精度较高，后续可用于地温场的分布预测。

图3 不同地形区域中随机选取的钻孔点位分布Fig.3 Distribution of boreholes selected randomly in different terrain areas

3.2 研究区冻土地温分布特征

从图4中可以发现，研究区不同地形区域多年冻土地温在5～15 m的深度范围内，地温分布变化不明显，地温曲线呈垂直形态。因此，本文绘制了研究区域在2001年7月至9月间多年冻土地下0.2～5.0 m的地温分布情况，如图5所示。研究区内山区地温显著低于平原和盆地地温，风火山、可可西里山和昆仑山越岭区的地温相对较低；而尺曲谷地、北麓河盆地以及楚玛尔河高平原地区的地温相对较高。通过对比图5中7月至9月的地温分布情况，发现最北端的昆仑山越岭区在0.2～5.0 m的深度区间内，地温随时间推移，降低得最快；而最南端尺曲谷地的地温则降低得最慢。楚玛尔河高平原区域海拔变化不大，地温变化并不明显，可能是由于该区域内分布有大量的季节性河流，对多年冻土具有良好的保温作用。

图5 研究区域2001年7月（a）、8月（b）和9月（c）的地温空间分布Fig.5 Distribution of ground temperature in the study area in July（a），August（b）and September（c），2001

在上述结果基础上，通过提取并绘制0℃等温线（图6），获取了研究区域多年冻土的融化深度，即多年冻土活动层厚度。在风火山、可可西里山和昆仑山的南坡，0℃等温线距离较大；而在北坡，0℃等温线距离较小。研究区域内多年冻土融化深度从7月到9月逐渐增加。其中，风火山南北坡向地区的冻土融化深度变化最为明显，可能的原因在于该地区是整个研究区内植被发育最好的地区。植被在一定程度上可以隔绝冷空气，减少表土层水分蒸发，且因该地区坡度较大，接受太阳辐射强，易引起冻融变化［46］。相反，在楚玛尔河地区，3个月的0℃等温线变化较小，基本重合。值得注意的是7月0℃等温线在五道梁楚玛尔河等地区超过了8月和9月对应的0℃等温线，可能是由于该地区海拔下降，坡向向北，7月能接受的太阳光辐射强度最大，而8月和9月太阳逐渐转回赤道，该地区能接受到的太阳辐射变少所致；而楚玛尔河至昆仑山区域，海拔上升，坡向向南，可以持续地接受太阳辐射，冻土融化深度逐渐增加。

图6 研究区域2001年7月至9月的0℃等温线Fig.6 The 0℃isotherm in the study area in July（blue），August（red）and September（green），2001

本文还对不同地形区域多年冻土融化深度的平均值进行了统计，结果如表1所示。昆仑山地区多年冻土融化深度从7月到9月呈现先增加后减小的趋势，而其他地区多年冻土融化深度均呈现逐渐增大趋势。楚玛尔河高平原和可可西里山两个区域多年冻土融化深度变化不明显。该结果与Yin等［21］2018年对青藏高原工程走廊多年冻土活动层厚度模拟结果相似（即在高山地区平均ALT小于2.0 m，而河谷的ALT大于4.0 m，高地平原地区，ALT的范围通常为3.0～4.0 m）。

表1 研究区域内不同地形地区多年冻土融化深度变化（单位：m）Table 1 Variation of thawing depth of permafrost in different terrain areas in the study area（unit：m）

3.3 讨论

在充分考虑影响因素间存在的多重共线性和地温与影响因素间的空间非平稳性特征的前提下，本文将GWRRK模型与具有外部漂移克里金方法（KED）以及地理加权岭回归模型（GWRR）的模拟结果进行了对比，如表2所示。GWRRK在3个指标的评价结果均优于GWRR和KED，说明该方法具有更好的模拟效果。

表2 不同模型的精度对比Table 2 Accuracy comparison of different models

同时，本文也与相似区域的前人研究结果进行了对比分析。GWRRK模型的整体模拟精度（R2=0.87）略高于Yin等［21］利用GIPL2.0模型对青藏高原工程走廊年平均地温的模拟精度（R2=0.85），但略低于Sun等［47］对青藏高原多年冻土地表温度的模拟精度（R2=0.91），可能的原因一方面在于建模参数选择的差异化，如本文不仅考虑了土壤和气温的影响，还充分考虑了植被覆盖、降雨量和风速的影响；另一方面本文研究的区域尺度相对较小，且钻孔点位在区域内分布较为密集，是否可在更大的区域尺度上应用，值得进一步探索。此外，Hinzman等［48］利用一维非稳态热传导有限元方法观察冻土土壤温度剖面与融化深度的空间分布时，发现不同深度地温的模拟结果与实测地温结果存在一定差异。本文结果与之相近，即在深度0～5 m的区间内，模型拟合值与实测值相差较大，但随着深度的增加，拟合度逐步提高。可能原因在于土壤温度变化主要受热传递影响，浅层土壤同时还受降雨、植被、风速等影响，在短时间内对地温的热传导作用较大；而在较深层土壤或富含冰的土壤中，热传导受上述作用较小［33，49］。

模型回归系数可以反映各因素对多年冻土地温的影响程度，其正负表示影响因素与地温呈正相关或负相关，绝对值越大则说明该因素影响越大。由图4不难发现，在5～15 m深度区间地温曲线呈垂直形态，地温变化很小。为此，在分析地温影响因素时，本文重点关注0～5 m深度区间。图7反映了7—9月0.5～5 m深度区间各变量回归系数的统计平均值，可发现各影响因素系数在3 m以下的深度开始逐渐收敛。根据回归系数的绝对值大小发现，研究区多年冻土地温受到的主要影响因素依次是地面向下长波辐射、地面向下短波辐射、纬度和气温。

图4 研究区域内不同地形位置钻孔点地温模拟值与实测值对比Fig.4 Comparison between simulation results and borehole ground temperatures in different terrain areas：Kunlun Mountain（a），Chumar River High Plain（b），Hoh Xil Mountain（c），Beilu River Basin（d），Fenghuo Mountain（e）and Chiqu Valley（f）

图7 GWRRK模型回归系数平均值变化Fig.7 Change of mean value of regression coefficient of GWRRK model：in July（a），August（b）and September（c）

就同一因素，在不同月份对多年冻土地温的影响程度也存在显著差异。如地面向下长波辐射在7月对地温的影响呈正相关，而在8月的影响呈负相关；9月在0.5～3.0 m的深度区间呈正相关，而在3.5～5.0 m的深度区间呈负相关。气温在7月和9月呈现负相关，而在8月呈现正相关，其可能原因在于7月气温在研究区域达到最大，而地温向下传递的过程仍存在一定的滞后性，故呈现负相关；8月地温在持续升温中，大气温度还未开始明显下降，故地温与气温呈现明显正相关；到9月，大气温度开始下降，地温仍处于持续升温过程，故9月地温与气温再次呈现负相关。

此外，本文更多的是从数据驱动建模的角度对研究区地温分布模拟开展研究，在模型的应用、影响因素的选取和实验中存在着一定的不确定性。

首先，文中通过岭回归模型消除影响因素间的多重共线性，方差膨胀系数VIF值小于10即可认为变量间不存在多重共线性。本文选取所有变量VIFmax值小于5来确定模型岭回归参数λ，可以明确模型消除了多重共线性问题，但却无法量化岭回归参数λ对模拟结果的影响。其次，由于地理位置的特殊性，本文的气象数据采用青藏高原科学数据中心公开发布的中国区域地面气象要素驱动数据集插值而得到，但未对插值后的数据进行评价，可能导致其对模拟结果产生偏差。此外，在结果阐释方面，研究区内高温多年冻土和低温多年冻土可能同时存在，与自然状况下的多年冻土地温曲线有所不同，可能与铁路修建引起的热扰动相关，但其影响的深层原因、作用机制等还有待进一步探究分析。后续研究考虑进一步优化和改进模型，将数据驱动模型和多年冻土形成的机理知识结合，将不同区域的模拟结果整体对比，以验证模拟结果的区域差异；同时从不同地形环境、不同时段和不同层次的角度，对研究区多年冻土地温的影响机理进行深入分析，聚焦地温的水平与垂直分布特征分析，以揭示不同因子对地温空间及深度分布的贡献率。

4 结论

本文以青藏铁路昆仑山至尺曲谷地段为例，利用GWRRK方法对铁路沿线多年冻土地温的空间分布进行了模拟，分析了2001年7月至9月研究区域内多年冻土地温的分布变化特征，获取了该地区多年冻土的融化深度。

结果表明，研究区域内多年冻土地温总体呈现山区地温低于平原和盆地地温的特征；地温随冻土深度增加而降低，地温在深度0～5 m的区间内温度变化较大，平均温差为10.3℃，而在5～15 m的深度区间基本保持不变，平均温差为0.2℃。研究区域最北端的昆仑山越岭区地温下降最快，最南端尺曲谷地的地温下降最慢；楚玛尔河高平原区域受河流湖泊的影响，地温和冻土活动层厚度的变化不明显，而风火山南北坡向地区的冻土融化深度变化最显著。本文将GWRRK方法的模拟效果与具有外部漂移克里金方法（KED）和地理加权岭回归（GWRR）进行了对比，反映出GWRRK方法具有较高的精度优势，可通过该方法弥补地面监测的信息缺口。