基于GWRRK的多年凍土地溫空間分布模擬研究
——以青藏鐵路昆侖山至尺曲谷地段為例

姚明星，趙 銳，楊林川，齊 華，廖 昕，孟祥連，周福軍

（1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756；2.西南交通大學(xué)建筑學(xué)院，四川 成都 611756；3.中鐵第一勘察設(shè)計(jì)院集團(tuán)有限公司，陜西 西安 710043）

0 引言

隨著全球氣候變化加劇，多年凍土面臨退化趨勢(shì)，可能對(duì)寒區(qū)的生態(tài)安全產(chǎn)生影響，誘發(fā)諸如崩塌、滑坡、泥石流等地質(zhì)災(zāi)害問(wèn)題［1-3］。多年凍土對(duì)地溫變化極其敏感［4］，通過(guò)研究其地溫空間分布，可探索凍土活動(dòng)層厚度的變化規(guī)律，分析凍土熱融狀況以及未來(lái)變化趨勢(shì)，為凍土災(zāi)害防治提供科學(xué)依據(jù)。

多年凍土地溫既受海拔、緯度和經(jīng)度的影響，也受氣象、地形地貌、植被、積雪、土壤巖性、含水量等局地因素的影響［5-10］。學(xué)者基于上述影響因素，結(jié)合凍融指數(shù)、熱傳導(dǎo)、地理加權(quán)回歸、支持向量機(jī)、機(jī)器學(xué)習(xí)等方法開(kāi)展了多年凍土地溫的時(shí)空分布模擬研究［11-16］。例如，Jafarov等［11］、Nicolsky等［17］利用GIPL（geophysical institute permafrost laboratory）模型，考慮氣溫、降水、土壤熱特性和土壤含水量的影響，基于含相變的一維非線性熱傳導(dǎo)方程，模擬了阿拉斯加州多年凍土的季節(jié)凍融深度。Liu等［18］利用二維非穩(wěn)態(tài)模型分析了青藏高原多年凍土區(qū)路基溫度，模擬了凍土土壤與大氣的傳熱過(guò)程。Luo等［19］綜合考慮了氣溫、植被、積雪深度和土壤屬性等因素，對(duì)黃河源區(qū)多年凍土地溫的空間分布進(jìn)行了模擬。Qin等［20］、Yin等［21］利用GIPL模型分別對(duì)青藏高原多年凍土活動(dòng)層厚度和青藏高原工程走廊年平均溫度進(jìn)行模擬，前者基于土壤導(dǎo)熱系數(shù)開(kāi)展多年凍土活動(dòng)層厚度模擬，后者側(cè)重于土壤性質(zhì)對(duì)多年凍土地溫的影響。Liu等［22］綜合考慮氣溫、風(fēng)速、風(fēng)向、地表水分蒸發(fā)量、太陽(yáng)輻射等影響因素，基于二維非穩(wěn)態(tài)熱傳導(dǎo)有限元模型，揭示了凍土路基在給定工況下的地溫時(shí)空分布特征。Zhao等［15］利用地理加權(quán)回歸模型對(duì)青藏鐵路沿線多年凍土地溫分布進(jìn)行模擬，同時(shí)指出各因素對(duì)地溫的影響隨時(shí)空尺度變化增大，呈現(xiàn)出典型的空間非平穩(wěn)性特征。此外，各影響因素間可能還存在著交互影響，易引起多重共線性問(wèn)題［15］。

本研究利用地理加權(quán)嶺回歸克里金方法（geographically weighted ridge regression Kriging，GWRRK）構(gòu)建凍土地溫分布模型，以青藏鐵路昆侖山至尺曲谷地段多年凍土覆蓋區(qū)域?yàn)槔?yàn)證模型的有效性，通過(guò)對(duì)該區(qū)域多年凍土地溫分布的揭示，分析凍土融化深度的變化特征。地理加權(quán)嶺回歸克里金方法（GWRRK）既可以充分考慮多年凍土地溫與其影響因素間的空間非平穩(wěn)性，也可消除影響因素之間的多重共線性［23］，同時(shí)還考慮了擬合殘差的影響。昆侖山至尺曲谷地段是青藏鐵路建設(shè)和運(yùn)營(yíng)中凍土變化監(jiān)測(cè)最密集的區(qū)域，通過(guò)揭示該區(qū)域凍土地溫空間分布變化趨勢(shì)，可為鐵路安全運(yùn)營(yíng)及沿線地質(zhì)災(zāi)害防治提供科學(xué)依據(jù)。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于青藏高原中部，隸屬于青海省玉樹(shù)藏族自治州，北至昆侖山越嶺區(qū)，南到尺曲谷地，橫跨楚瑪爾河高平原、可可西里山區(qū)、北麓河盆地、風(fēng)火山山區(qū)，該地區(qū)是典型的高海拔多年凍土區(qū)域。研究區(qū)平均海拔4 900 m，區(qū)內(nèi)山系呈東西走向，山系與河流交替出現(xiàn)，形成高平原、盆地和谷地等地貌。河流主要有楚瑪爾河、北麓河和通天河，由冰川融水和大氣降水補(bǔ)給，流量季節(jié)性變化較大。研究區(qū)具有獨(dú)特的冰緣干旱氣候特征，且隨海拔增高有明顯的氣候垂直分帶性特征。

1.2 數(shù)據(jù)來(lái)源與處理

本文以地溫鉆孔實(shí)測(cè)數(shù)據(jù)為基礎(chǔ)，從地理要素、氣象條件和局地因素三個(gè)維度選取地溫影響因素，構(gòu)建凍土地溫分布模型。其中，地理要素主要考慮緯度和高程［16，24-25］；氣象條件選擇近地表氣溫、風(fēng)速、降水率、地面向下長(zhǎng)波輻射和地面向下短波輻射［26-29］；局地因素則考慮地層巖性、坡度、坡向、植被覆蓋度［30-33］。

青藏高原氣溫往往在7月達(dá)到最大，且由于向下傳遞時(shí)存在一定滯后性［34-35］，本文重點(diǎn)對(duì)研究區(qū)7—9月的地溫空間分布開(kāi)展模擬，分析其對(duì)應(yīng)的凍土活動(dòng)層厚度的空間變化特征。本文實(shí)測(cè)地溫?cái)?shù)據(jù)是來(lái)自2001年修建青藏鐵路工程的鉆孔地溫?cái)?shù)據(jù)（采用熱敏電阻測(cè)溫，測(cè)溫深度為0.2～15 m，數(shù)據(jù)采集間隔為10天至一年不等）。在研究區(qū)域內(nèi)所獲取的鉆孔點(diǎn)位一共237個(gè)，7月至9月的鉆孔地溫?cái)?shù)據(jù)點(diǎn)位分布較為均勻，共307次測(cè)溫?cái)?shù)據(jù)。其中，7月測(cè)溫96次，8月測(cè)溫128次，9月測(cè)溫83次。通過(guò)插值，獲取時(shí)間序列數(shù)據(jù)進(jìn)行對(duì)比分析，可更好了解7月至9月研究區(qū)內(nèi)整個(gè)路段多年凍土地溫變化及凍土融化深度變化情況。鉆孔位置從北向南穿越不同地形區(qū)域，基本反映了青藏高原中低緯度、高海拔的高原凍土特征，樣本點(diǎn)具有較好的代表性。研究區(qū)及鉆孔測(cè)溫點(diǎn)位分布如圖1所示，其中左上角小圖中的多年凍土分布矢量圖來(lái)自Ran等［36］的研究成果。

圖1 研究區(qū)及鉆孔點(diǎn)位分布Fig.1 Study area and distribution of the boreholes

由于研究區(qū)內(nèi)僅有五道梁一個(gè)氣象觀測(cè)站，難以滿足模擬的空間分辨率要求。因此，本文氣象數(shù)據(jù)采用青藏高原科學(xué)數(shù)據(jù)中心公開(kāi)發(fā)布的中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集（1979—2018年），涉及近地面氣溫、地面降水率、近地面全風(fēng)速、地面向下短波輻射和地面向下長(zhǎng)波輻射共5個(gè)指標(biāo)［37-38］。上述數(shù)據(jù)水平空間分辨率為0.1°。本文利用ArcGIS 10.2軟件對(duì)上述數(shù)據(jù)進(jìn)行柵格轉(zhuǎn)點(diǎn)矢量文件處理，并綜合對(duì)比克里金法、反距離權(quán)重法、自然鄰域法和趨勢(shì)面法所獲取的1 km×1 km分辨率氣象數(shù)據(jù)結(jié)果，選用擬合精度最優(yōu)的克里金插值結(jié)果作為氣象基礎(chǔ)數(shù)據(jù)。地層巖性數(shù)據(jù)通過(guò)全國(guó)2002年發(fā)布的1∶100萬(wàn)地質(zhì)圖裁剪獲取。地層巖性因素是繪制不同地層巖性區(qū)域的實(shí)測(cè)地溫?cái)?shù)據(jù)與其深度的擬合曲線，計(jì)算曲線的拐點(diǎn)坐標(biāo)，利用拐點(diǎn)的深度坐標(biāo)值作為地層巖性量化因子。坡向和坡度數(shù)據(jù)是利用ArcGIS 10.2軟件對(duì)研究區(qū)DEM（分辨率30 m）計(jì)算。植被指數(shù)來(lái)源于Terra衛(wèi)星MOD13 Q1全球250 m分辨率植被指數(shù)16天的合成產(chǎn)品，并利用像元二分模型計(jì)算獲取植被覆蓋度［39］。

2 模型構(gòu)建

2.1 地理加權(quán)嶺回歸模型

地理加權(quán)回歸模型（GWR）是在普通線性回歸模型的基礎(chǔ)上，將空間位置信息嵌入到回歸中，以預(yù)測(cè)每個(gè)位置對(duì)應(yīng)因變量的變化，表達(dá)式［40-41］為

式中：(ui，vi)為第i個(gè)采樣點(diǎn)的坐標(biāo)（如經(jīng)緯度）；p為影響因素的個(gè)數(shù)，本文共考慮了11個(gè)影響因子，故p=11；βk(ui，vi)為第i個(gè)采樣點(diǎn)上的第k個(gè)回歸參 數(shù)，是 關(guān) 于 地 理 位 置 的 函 數(shù)，εi～N(0，σ2)，Cov(εi，εj)=0(i≠j)。

為解決地理加權(quán)回歸中可能存在的因子間多重共線性問(wèn)題，Wheeler［42］將地理加權(quán)回歸（GWR）與嶺回歸結(jié)合，提出了地理加權(quán)嶺回歸模型（GWRR），其回歸參數(shù)可估計(jì)為

式中：λ為嶺回歸參數(shù)。本文的嶺回歸參數(shù)λ通過(guò)設(shè)置方差膨脹系數(shù)閾值小于5（VIF＜5）確定得到。

在地理加權(quán)嶺回歸模型中，第i個(gè)點(diǎn)上的因變量擬合值［42］為

2.2 地理加權(quán)嶺回歸克里金模型

地理加權(quán)嶺回歸克里金模型（GWRRK）是將GWRR與普通克里金插值方法（ordinary Kriging，OK）結(jié)合，即在GWRR擬合后，利用OK處理回歸中的殘差問(wèn)題，以此來(lái)提高GWRR的預(yù)測(cè)精度，表達(dá)式［23］為

第i處地理加權(quán)嶺回歸模型（GWRR）的模擬值表述［23］為

樣本點(diǎn)數(shù)值xi(i=1，2，…，n)與插值點(diǎn)的空間區(qū)域變化函數(shù)為Z(xi)，則GWRR第i處局部回歸殘差項(xiàng)的克里金插值模擬值可表達(dá)［23］為

上式必須滿足以下條件，

式中：C(xi，xj)為采樣點(diǎn)xi與xj之間的變異函數(shù)值；μ為拉格朗日常數(shù)。

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

具有外部漂移的克里金模型（Kriging with an external drift，KED）是帶有趨勢(shì)克里金模型的拓展，既考慮了樣本點(diǎn)與待估計(jì)點(diǎn)之間的空間分布幾何特征，又考慮了樣本間空間關(guān)系的結(jié)構(gòu)特征［43］。

設(shè)u處屬性值Z(u)由漂移函數(shù)m(u)和殘差R(u)組成［43］，即

漂移函數(shù)通常可表示為

式中：β0、βl為未知系數(shù)；xl(u)為漂移函數(shù)的輔助變量。

樣本點(diǎn)數(shù)值xi(i=1，2，…，n)與插值點(diǎn)的空間區(qū)域變化函數(shù)為Z(xi)，則殘差R(u)可表達(dá)［23］為

上式必須滿足式（7）。

3 結(jié)果與討論

3.1 模型精度驗(yàn)證

本文采用留一交叉驗(yàn)證方法（leave-one-out cross-validation，LOO-CV）對(duì)模型結(jié)果進(jìn)行驗(yàn)證，其基本思想是每次從個(gè)數(shù)為N的樣本集中抽取出一個(gè)樣本作為驗(yàn)證集，剩下的N-1個(gè)樣本作為訓(xùn)練集，重復(fù)進(jìn)行N次，依次取遍所有N個(gè)數(shù)據(jù)作為驗(yàn)證集，最后將平均N個(gè)數(shù)據(jù)的結(jié)果作為泛化誤差的估計(jì)［44］。該方法較其他驗(yàn)證方法具有訓(xùn)練樣本量大、接近原始樣本分布的特點(diǎn)［45］。

圖2是研究區(qū)2001年7月至9月鉆孔點(diǎn)位在0.2～15 m深度上的留一交叉驗(yàn)證結(jié)果。可以看出，模型具有較高的擬合優(yōu)度，平均擬合優(yōu)度達(dá)0.89。其中，7月的擬合優(yōu)度最大，R2=0.91，且平均絕對(duì)誤差（mean absolute error，MAE）和均方根誤差（root mean square error，RMSE）相對(duì)較小；8月的擬合優(yōu)度最小，R2=0.87；9月擬合結(jié)果的均方根誤差（RMSE）值最小，達(dá)0.694。7月至9月精度驗(yàn)證結(jié)果的平均絕對(duì)誤差整體偏小，最大誤差為0.613℃，說(shuō)明本文所建GWRRK模型具有一定的穩(wěn)健性，可用于案例研究區(qū)域的多年凍土地溫分布模擬分析。此外發(fā)現(xiàn)，模擬地溫與實(shí)測(cè)地溫在-4～2.5℃的溫度區(qū)間范圍內(nèi)重合度較高，說(shuō)明在該溫度區(qū)間內(nèi)模擬的誤差相對(duì)較小，而在大于2.5℃的溫度區(qū)間內(nèi)偏差較大，可能由于不同地形區(qū)域內(nèi)地溫波動(dòng)較大所致。

圖2 GWRRK模型留一交叉驗(yàn)證結(jié)果Fig.2 Scatter distribution of LOO-CV results regarding the GWRRK model：in July（a），August（b）and September（c）

為進(jìn)一步驗(yàn)證模擬結(jié)果的有效性，本文分別從不同地形區(qū)域中隨機(jī)選取1個(gè)鉆孔測(cè)溫點(diǎn)，點(diǎn)位分布如圖3所示，將其對(duì)應(yīng)的模擬結(jié)果與實(shí)測(cè)結(jié)果對(duì)比。由圖4可知，在不同區(qū)域的地溫模擬值與實(shí)測(cè)值的擬合精度R2都大于0.95；均方根誤差（RMSE）最小達(dá)0.149。此外，除圖4（a）昆侖山區(qū)域7月的模擬值與實(shí)測(cè)值有一定偏差外，其他位置不同時(shí)間段的模擬結(jié)果與實(shí)測(cè)結(jié)果的變化趨勢(shì)幾乎一致。這也說(shuō)明GWRRK在該研究區(qū)域尺度內(nèi)的模擬精度較高，后續(xù)可用于地溫場(chǎng)的分布預(yù)測(cè)。

圖3 不同地形區(qū)域中隨機(jī)選取的鉆孔點(diǎn)位分布Fig.3 Distribution of boreholes selected randomly in different terrain areas

3.2 研究區(qū)凍土地溫分布特征

從圖4中可以發(fā)現(xiàn)，研究區(qū)不同地形區(qū)域多年凍土地溫在5～15 m的深度范圍內(nèi)，地溫分布變化不明顯，地溫曲線呈垂直形態(tài)。因此，本文繪制了研究區(qū)域在2001年7月至9月間多年凍土地下0.2～5.0 m的地溫分布情況，如圖5所示。研究區(qū)內(nèi)山區(qū)地溫顯著低于平原和盆地地溫，風(fēng)火山、可可西里山和昆侖山越嶺區(qū)的地溫相對(duì)較低；而尺曲谷地、北麓河盆地以及楚瑪爾河高平原地區(qū)的地溫相對(duì)較高。通過(guò)對(duì)比圖5中7月至9月的地溫分布情況，發(fā)現(xiàn)最北端的昆侖山越嶺區(qū)在0.2～5.0 m的深度區(qū)間內(nèi)，地溫隨時(shí)間推移，降低得最快；而最南端尺曲谷地的地溫則降低得最慢。楚瑪爾河高平原區(qū)域海拔變化不大，地溫變化并不明顯，可能是由于該區(qū)域內(nèi)分布有大量的季節(jié)性河流，對(duì)多年凍土具有良好的保溫作用。

圖5 研究區(qū)域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

在上述結(jié)果基礎(chǔ)上，通過(guò)提取并繪制0℃等溫線（圖6），獲取了研究區(qū)域多年凍土的融化深度，即多年凍土活動(dòng)層厚度。在風(fēng)火山、可可西里山和昆侖山的南坡，0℃等溫線距離較大；而在北坡，0℃等溫線距離較小。研究區(qū)域內(nèi)多年凍土融化深度從7月到9月逐漸增加。其中，風(fēng)火山南北坡向地區(qū)的凍土融化深度變化最為明顯，可能的原因在于該地區(qū)是整個(gè)研究區(qū)內(nèi)植被發(fā)育最好的地區(qū)。植被在一定程度上可以隔絕冷空氣，減少表土層水分蒸發(fā)，且因該地區(qū)坡度較大，接受太陽(yáng)輻射強(qiáng)，易引起凍融變化［46］。相反，在楚瑪爾河地區(qū)，3個(gè)月的0℃等溫線變化較小，基本重合。值得注意的是7月0℃等溫線在五道梁楚瑪爾河等地區(qū)超過(guò)了8月和9月對(duì)應(yīng)的0℃等溫線，可能是由于該地區(qū)海拔下降，坡向向北，7月能接受的太陽(yáng)光輻射強(qiáng)度最大，而8月和9月太陽(yáng)逐漸轉(zhuǎn)回赤道，該地區(qū)能接受到的太陽(yáng)輻射變少所致；而楚瑪爾河至昆侖山區(qū)域，海拔上升，坡向向南，可以持續(xù)地接受太陽(yáng)輻射，凍土融化深度逐漸增加。

圖6 研究區(qū)域2001年7月至9月的0℃等溫線Fig.6 The 0℃isotherm in the study area in July（blue），August（red）and September（green），2001

本文還對(duì)不同地形區(qū)域多年凍土融化深度的平均值進(jìn)行了統(tǒng)計(jì)，結(jié)果如表1所示。昆侖山地區(qū)多年凍土融化深度從7月到9月呈現(xiàn)先增加后減小的趨勢(shì)，而其他地區(qū)多年凍土融化深度均呈現(xiàn)逐漸增大趨勢(shì)。楚瑪爾河高平原和可可西里山兩個(gè)區(qū)域多年凍土融化深度變化不明顯。該結(jié)果與Yin等［21］2018年對(duì)青藏高原工程走廊多年凍土活動(dòng)層厚度模擬結(jié)果相似（即在高山地區(qū)平均ALT小于2.0 m，而河谷的ALT大于4.0 m，高地平原地區(qū)，ALT的范圍通常為3.0～4.0 m）。

表1 研究區(qū)域內(nèi)不同地形地區(qū)多年凍土融化深度變化（單位：m）Table 1 Variation of thawing depth of permafrost in different terrain areas in the study area（unit：m）

3.3 討論

在充分考慮影響因素間存在的多重共線性和地溫與影響因素間的空間非平穩(wěn)性特征的前提下，本文將GWRRK模型與具有外部漂移克里金方法（KED）以及地理加權(quán)嶺回歸模型（GWRR）的模擬結(jié)果進(jìn)行了對(duì)比，如表2所示。GWRRK在3個(gè)指標(biāo)的評(píng)價(jià)結(jié)果均優(yōu)于GWRR和KED，說(shuō)明該方法具有更好的模擬效果。

表2 不同模型的精度對(duì)比Table 2 Accuracy comparison of different models

同時(shí)，本文也與相似區(qū)域的前人研究結(jié)果進(jìn)行了對(duì)比分析。GWRRK模型的整體模擬精度（R2=0.87）略高于Yin等［21］利用GIPL2.0模型對(duì)青藏高原工程走廊年平均地溫的模擬精度（R2=0.85），但略低于Sun等［47］對(duì)青藏高原多年凍土地表溫度的模擬精度（R2=0.91），可能的原因一方面在于建模參數(shù)選擇的差異化，如本文不僅考慮了土壤和氣溫的影響，還充分考慮了植被覆蓋、降雨量和風(fēng)速的影響；另一方面本文研究的區(qū)域尺度相對(duì)較小，且鉆孔點(diǎn)位在區(qū)域內(nèi)分布較為密集，是否可在更大的區(qū)域尺度上應(yīng)用，值得進(jìn)一步探索。此外，Hinzman等［48］利用一維非穩(wěn)態(tài)熱傳導(dǎo)有限元方法觀察凍土土壤溫度剖面與融化深度的空間分布時(shí)，發(fā)現(xiàn)不同深度地溫的模擬結(jié)果與實(shí)測(cè)地溫結(jié)果存在一定差異。本文結(jié)果與之相近，即在深度0～5 m的區(qū)間內(nèi)，模型擬合值與實(shí)測(cè)值相差較大，但隨著深度的增加，擬合度逐步提高。可能原因在于土壤溫度變化主要受熱傳遞影響，淺層土壤同時(shí)還受降雨、植被、風(fēng)速等影響，在短時(shí)間內(nèi)對(duì)地溫的熱傳導(dǎo)作用較大；而在較深層土壤或富含冰的土壤中，熱傳導(dǎo)受上述作用較小［33，49］。

模型回歸系數(shù)可以反映各因素對(duì)多年凍土地溫的影響程度，其正負(fù)表示影響因素與地溫呈正相關(guān)或負(fù)相關(guān)，絕對(duì)值越大則說(shuō)明該因素影響越大。由圖4不難發(fā)現(xiàn)，在5～15 m深度區(qū)間地溫曲線呈垂直形態(tài)，地溫變化很小。為此，在分析地溫影響因素時(shí)，本文重點(diǎn)關(guān)注0～5 m深度區(qū)間。圖7反映了7—9月0.5～5 m深度區(qū)間各變量回歸系數(shù)的統(tǒng)計(jì)平均值，可發(fā)現(xiàn)各影響因素系數(shù)在3 m以下的深度開(kāi)始逐漸收斂。根據(jù)回歸系數(shù)的絕對(duì)值大小發(fā)現(xiàn)，研究區(qū)多年凍土地溫受到的主要影響因素依次是地面向下長(zhǎng)波輻射、地面向下短波輻射、緯度和氣溫。

圖4 研究區(qū)域內(nèi)不同地形位置鉆孔點(diǎn)地溫模擬值與實(shí)測(cè)值對(duì)比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），F(xiàn)enghuo Mountain（e）and Chiqu Valley（f）

圖7 GWRRK模型回歸系數(shù)平均值變化Fig.7 Change of mean value of regression coefficient of GWRRK model：in July（a），August（b）and September（c）

就同一因素，在不同月份對(duì)多年凍土地溫的影響程度也存在顯著差異。如地面向下長(zhǎng)波輻射在7月對(duì)地溫的影響呈正相關(guān)，而在8月的影響呈負(fù)相關(guān)；9月在0.5～3.0 m的深度區(qū)間呈正相關(guān)，而在3.5～5.0 m的深度區(qū)間呈負(fù)相關(guān)。氣溫在7月和9月呈現(xiàn)負(fù)相關(guān)，而在8月呈現(xiàn)正相關(guān)，其可能原因在于7月氣溫在研究區(qū)域達(dá)到最大，而地溫向下傳遞的過(guò)程仍存在一定的滯后性，故呈現(xiàn)負(fù)相關(guān)；8月地溫在持續(xù)升溫中，大氣溫度還未開(kāi)始明顯下降，故地溫與氣溫呈現(xiàn)明顯正相關(guān)；到9月，大氣溫度開(kāi)始下降，地溫仍處于持續(xù)升溫過(guò)程，故9月地溫與氣溫再次呈現(xiàn)負(fù)相關(guān)。

此外，本文更多的是從數(shù)據(jù)驅(qū)動(dòng)建模的角度對(duì)研究區(qū)地溫分布模擬開(kāi)展研究，在模型的應(yīng)用、影響因素的選取和實(shí)驗(yàn)中存在著一定的不確定性。

首先，文中通過(guò)嶺回歸模型消除影響因素間的多重共線性，方差膨脹系數(shù)VIF值小于10即可認(rèn)為變量間不存在多重共線性。本文選取所有變量VIFmax值小于5來(lái)確定模型嶺回歸參數(shù)λ，可以明確模型消除了多重共線性問(wèn)題，但卻無(wú)法量化嶺回歸參數(shù)λ對(duì)模擬結(jié)果的影響。其次，由于地理位置的特殊性，本文的氣象數(shù)據(jù)采用青藏高原科學(xué)數(shù)據(jù)中心公開(kāi)發(fā)布的中國(guó)區(qū)域地面氣象要素驅(qū)動(dòng)數(shù)據(jù)集插值而得到，但未對(duì)插值后的數(shù)據(jù)進(jìn)行評(píng)價(jià)，可能導(dǎo)致其對(duì)模擬結(jié)果產(chǎn)生偏差。此外，在結(jié)果闡釋方面，研究區(qū)內(nèi)高溫多年凍土和低溫多年凍土可能同時(shí)存在，與自然狀況下的多年凍土地溫曲線有所不同，可能與鐵路修建引起的熱擾動(dòng)相關(guān)，但其影響的深層原因、作用機(jī)制等還有待進(jìn)一步探究分析。后續(xù)研究考慮進(jìn)一步優(yōu)化和改進(jìn)模型，將數(shù)據(jù)驅(qū)動(dòng)模型和多年凍土形成的機(jī)理知識(shí)結(jié)合，將不同區(qū)域的模擬結(jié)果整體對(duì)比，以驗(yàn)證模擬結(jié)果的區(qū)域差異；同時(shí)從不同地形環(huán)境、不同時(shí)段和不同層次的角度，對(duì)研究區(qū)多年凍土地溫的影響機(jī)理進(jìn)行深入分析，聚焦地溫的水平與垂直分布特征分析，以揭示不同因子對(duì)地溫空間及深度分布的貢獻(xiàn)率。

4 結(jié)論

本文以青藏鐵路昆侖山至尺曲谷地段為例，利用GWRRK方法對(duì)鐵路沿線多年凍土地溫的空間分布進(jìn)行了模擬，分析了2001年7月至9月研究區(qū)域內(nèi)多年凍土地溫的分布變化特征，獲取了該地區(qū)多年凍土的融化深度。

結(jié)果表明，研究區(qū)域內(nèi)多年凍土地溫總體呈現(xiàn)山區(qū)地溫低于平原和盆地地溫的特征；地溫隨凍土深度增加而降低，地溫在深度0～5 m的區(qū)間內(nèi)溫度變化較大，平均溫差為10.3℃，而在5～15 m的深度區(qū)間基本保持不變，平均溫差為0.2℃。研究區(qū)域最北端的昆侖山越嶺區(qū)地溫下降最快，最南端尺曲谷地的地溫下降最慢；楚瑪爾河高平原區(qū)域受河流湖泊的影響，地溫和凍土活動(dòng)層厚度的變化不明顯，而風(fēng)火山南北坡向地區(qū)的凍土融化深度變化最顯著。本文將GWRRK方法的模擬效果與具有外部漂移克里金方法（KED）和地理加權(quán)嶺回歸（GWRR）進(jìn)行了對(duì)比，反映出GWRRK方法具有較高的精度優(yōu)勢(shì)，可通過(guò)該方法彌補(bǔ)地面監(jiān)測(cè)的信息缺口。