青藏工程走廊活動層厚度預測模型與分布特征研究

劉志云， 黃 川， 于 暉， 鐘振濤， 崔福慶

（1.長安大學地質工程與測繪學院，陜西西安 710054； 2.中交第一公路勘察設計研究院有限公司西安中交公路巖土工程有限責任公司，陜西西安 710075）

0 引言

青藏高原被稱為地球的“第三極”，其內發育著世界上中、低緯度帶海拔最高、面積最廣的凍土區［1］。隨著全球氣候變暖，青藏高原多年凍土產生嚴重的退化現象，主要表現為地溫升高、活動層厚度增加以及多年凍結層消失等［2-4］。活動層是指暖季融化冷季凍結覆蓋于多年凍土之上的土層，是凍土地層內水熱交換最主要的區域［5］，受經緯度、高程、植被覆蓋度、地表溫度、土壤性質及氣候環境等諸多因素的影響。活動層厚度的變化將會對寒區水文、地質、環境和工程建筑產生一系列影響。因此，建立活動層厚度的預測模型以及研究活動層厚度的分布特征對青藏高原地區工程構筑物的設計、建造及后期養護具有重要意義。

活動層厚度的現場測量主要有鉆孔、觸探、挖探等機械方法以及監測站監測等方法。Hon 等［6］利用動態圓錐貫入儀在活動層與多年凍土層的動力錐穿透指數的不同，從而測定活動層厚度。更多學者則是通過不同地區的現場測量數據對活動層厚度的年際變化進行研究分析［7-9］，Zhao 等［10］根據青藏高原鉆孔溫度曲線證明，在1967—1997 年期間，其活動層的厚度以平均速度為0.71 cm·a-1增長。Wu 等［11］通過對青藏公路沿線監測站監測數據分析研究，發現從1995—2007 年期間，青藏高原多年凍土區活動層厚度正以平均7.5 cm·a-1的速率持續增加。活動層厚度現場測量雖然能夠獲得較高的精確度，但存在著成本高、數據樣本少且難以連續觀測致使其無法刻畫大區域特征等缺點。

針對以上不足，許多學者發現經驗公式模型、統計模型和數值模型對活動層厚度大區域分布研究具有較好的效果。Pang 等［12-13］用Stefan 和Kudryavtsev 公式計算了青藏高原活動層厚度的變化，給出了青藏高原活動層厚度空間分布圖并預測了2049年、2099年的活動層厚度變化情況。Zhao等［14］將Kudryavtsev 公式與Lund-Postam-Jena 模型［15-16］耦合進一步提高了經驗公式模型對活動層厚度的預測能力。Ni 等［17］采用統計模型與機器學習結合的方式對青藏高原活動層厚度進行預測，研究發現模型比以往的研究具有較高的精度。Zhang 等［18］利用改進的GIPL2 模型對青藏高原無人區活動層厚度進行數值模擬，并基于實測活動層厚度數據驗證了數值模型的準確性和優越性。總體而言，經驗公式、統計模型和數值模型能夠有效的描繪大區域活動層厚度分布情況，但對局地因素考慮不全且未能結合實地勘查數據進行研究，致使空間分辨率低，難以反映實際情況。

近年來，隨著機器學習方法和遙感技術的快速發展，目前，諸多學者采用對現場勘查數據進行預測模型建立，再結合大范圍遙感數據進行區域分布模擬的方式已經廣泛應用于青藏高原凍土的相關研究，包括凍土分布［19-20］、凍土地溫預測［21-22］、凍土滑坡敏感性［23］等方面。鑒于此，本文通過青藏工程走廊沿線300 組活動層厚度鉆孔監測數據，基于年平均地表溫度、平均植被指數、等效緯度、緯度、高程和含冰量等參數建立了活動層厚度的經驗公式、隨機森林和RBF 神經網絡預測模型，通過對比三種預測模型的預測效果，并結合高精度遙感數據，運用預測精度最高的活動層厚度預測模型繪制了青藏工程走廊多年凍土區段沿線活動層厚度分布圖。

1 數據與研究方法

1.1 研究區概況

青藏工程走廊始于格爾木，止于拉薩市，橫穿青藏高原1 120 多公里，穿越多年凍土區約550 km，是內陸進入西藏的重要通道。本文以走廊內多年凍土區段（西大灘—安多）為研究區，研究活動層厚度分布狀況。如圖1所示，研究區以青藏公路、青藏鐵路為基準線向兩側外延10 km，全長約540 km，地理坐標位于32°～36°N、91°～95°E，海拔介于3 716～6 191 m，該區地貌類型豐富，包括中高山區、高平原、低山丘陵、河谷等。

圖1 研究區監測點分布Fig.1 Distribution of monitoring points in the study

1.2 數據源

活動層是地-氣間水熱交換的主要場所，大的氣候背景決定了大區域活動層厚度的宏觀分布狀況。但在一定條件下，局地因素的影響將超過氣候的影響，導致區域內相同氣候背景下活動層厚度的分布異常。局地因素主要通過影響太陽輻射、熱對流和熱傳導等過程從而影響活動層厚度的大小。因此，本文擬選用年平均地表溫度、平均植被指數、緯度、高程、等效緯度和含冰量六類數據作為預測模型建立的預測因子。

年平均地表溫度、平均植被指數和高程遙感數據由美國國家航天航空局下載的地表溫度（MOD11A2 H25V05，2000—2016 年）、植被指數（MOD13Q1 H25V05，2000—2016 年）和SRTMDEM（Shuttle Radar Topography Mission-Digital Elevation Model）數據產品中提取。等效緯度是表征太陽輻射對地表的影響，也是判斷坡面走向的重要數據，可通過由SRTM-DEM 數據產品中提取的坡度、坡向和緯度數據計算得到：

式中：φ′為等效緯度；l為坡度；h為坡向；φ為緯度。

凍土含冰量是多年凍土的基本特征指標之一，且活動層內不同含冰量也對活動層熱量吸收能力有一定影響。含冰量數據由現場鉆取多年凍土上限以下8 m 深度內的土層，再經過現場觀察記錄和室內試驗綜合確定。凍土含冰量鉆孔點沿青藏鐵路和青藏公路間隔2 km 布置，確定所有鉆孔點含冰量賦存狀態后，再利用概率插值得到走廊帶內凍土含冰量分布，如圖2所示。目前，含冰量分類標準主要以《凍土工程地質勘察規范》（50324—2014）為主［24］，分為少冰、多冰、富冰、飽冰和含土冰層五類。但在青藏高原地區，多年凍土上限以下土層內含冰量變化較為劇烈。因此，本文將根據多年凍土上限以下8 m 深度內土層含冰量主要的賦存類型將含冰量劃分為少冰-多冰、多冰-富冰、富冰-飽冰、飽冰-含土冰層四類。

圖2 研究區含冰量分布Fig.2 Distribution of ice content in the study area

用于建立預測模型的活動層厚度實測數據來源于青藏公路、青藏鐵路沿線監測斷面地溫監測數據（監測點如圖1），監測工作由中交第一公路勘察設計研究院有限公司高寒高海拔地區道路工程安全與健康國家重點實驗室格爾木觀測基地完成，數據時限為2006—2016 年，數據采集采用測溫法測得，精度為±0.05 ℃。活動層厚度監測點共計300組，在研究區內分布較為平均，基本體現了青藏高原工程走廊帶活動層厚度特征，具有較好的代表性。

1.3 預測因子間相關性分析

表1是各預測因子之間的相關性分析及共線性分析結果。由表可知，部分預測因子間雖然表現出較顯著的相關性，但線性相關關系較弱，說明各預測因子間雖然能夠相互影響，但影響程度較小。并且從共線性分析可以看出，各預測因子與活動層厚度之間的容差皆大于0.2 且膨脹方差因子（VIF）皆小于5，這進一步說明各預測因子之間不存在明顯的共線性。

表1 各因子之間的相關性及共線性分析Table 1 Correlation and collinearity analysis among factors

1.4 研究方法

1.4.1 經驗公式

最小二乘法是通過找尋數據誤差平方和的最小值，從而確定數據的最佳匹配函數。在本文中，含冰量數據作為分類變量，因此，在進行經驗公式擬合之前需對含冰量數據進行編碼處理。目前，對于分類變量的編碼方式主要有獨熱編碼、虛擬編碼和效應編碼等。其中，效應編碼具有不冗余、易于解釋的優點，其原理是在一個具有n個類別的變量中選取一個類別作為參照（賦值為-1），從而創建n－1 個指標變量（賦值為0 或1）。含冰量數據編碼處理方法如表2所示。

表2 含冰量數據效應編碼Table 2 Effect coding of ice content data

以上編碼代表當含冰量為少冰-多冰時，X1、X2、X3取值為-1、-1、-1；含冰量為多冰-富冰時，X1、X2、X3取值為1、0、0，以此類推。因此基于最小二乘法擬合的活動層厚度經驗公式如式（2）所示：

式中：h為活動層厚度；Ts為年平均地表溫度；N為平均植被指數；φ為緯度；H為高程；φ′為等效緯度。

1.4.2 機器學習

機器學習是現代智能技術中的一種重要方法，其主要研究從樣本數據中尋找規律，并根據這些規律預測未來或無法觀測的數據［25］。機器學習方法已經在青藏高原凍土相關研究中得以廣泛的應用并取得了較好的成果［26-28］。其中，隨機森林和徑向基函數神經網絡在處理復雜的非線性問題時具有良好的抗噪能力，因此，本文將選用隨機森林和RBF神經網絡建立活動層厚度預測模型。

隨機森林［29］是通過對多個決策樹（每棵樹都擬合到訓練數據的Bootstrap 樣本）求平均值來擬合組合模型，每棵樹中的每個拆分都考慮隨機的預測變量子集，通過這種方式，將多個弱模型組合起來生成更為強大的模型，其結構如圖3 所示。本文使用該方法建立了活動層厚度預測模型，并將數據集按8：2 比例進行隨機分塊，分別作為訓練集和驗證集，樣本個數分別為240和60。隨機森林預測模型中決策樹的數量為100 個，每棵樹的最小拆分大小為10個，最大拆分大小為2 000個。

圖3 隨機森林結構圖Fig.3 Random forest structure diagram

RBF 神經網絡是三層前饋神經網絡，由輸入層、隱含層和輸出層三部分構成。其中，輸入層用于接收外部信息，隱含層實現參數間的非線性轉換，輸出層用于輸出最終結果，其結構如圖4 所示。

圖4 神經網絡結構圖Fig.4 Neural network structure diagram

k-fold 交叉驗證法是一種能夠有效的利用小樣本數據得到最優模型的方法。其原理是將原始數據分成k組，將每一組數據分別做一次驗證集，其余的k-1 組數據作為訓練集，從而得到k個模型，最終選擇最優的模型。本文的RBF 神經網絡模型的結構設置為年平均地表溫度、平均植被指數、等效緯度、緯度、高程和含冰量作為輸入層，活動層厚度作為輸出層，高斯徑向基函數為激活函數。隱含層設置為兩層，每層12 個節點。k-fold 交叉驗證數為5，其中4 份用于模型訓練，1 份用于模型驗證，模型訓練和模型驗證樣本數分別為240和60。

1.4.3 評價指標

活動層厚度預測模型評價指標采用決定系數（R-Squared，R2）、平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）、均方根誤差（Root Mean Square Error，RMSE）和相對誤差±15%內占比。其中R2、MAPE、RMSE 表達式分別見式（3）～（5）。

2 結果與討論

2.1 模型結果對比

圖5 為三種活動層厚度預測模型的預測結果。由圖可知，經驗公式方法擬合效果最差，預測值大部分都分布于±15%誤差線外，并且可以發現在低活動層厚度區間預測值偏大，高活動層厚度區間預測值偏小，整體分布極為離散。隨機森林方法擬合效果次之，預測值大部分落在±15%誤差線內。但整體來看，樣本點仍較為離散，部分測點偏離實測值較大。說明隨機森林方法雖然能夠較好的預測活動層厚度，但存在一定誤差。RBF 神經網絡方法擬合效果最好，預測值基本位于±15%誤差線內，并且整體更靠近實測值線，尤其高活動層厚度區間。因此，從三種方法預測結果對比可知，基于RBF 神經網絡方法的活動層厚度預測模型預測效果最為精確、穩定。

圖5 各活動層厚度預測模型結果對比Fig.5 Comparison of results of prediction models for the thickness of each active layer

表3為三種預測模型的各項評價指標。由表可知，經驗公式預測模型的各項評價指標的評價優度均較差，R2僅為0.40，RMSE 和MAPE 達到0.75 和24.5%，相對誤差在±15%以內的占比也不足一半。隨機森林和RBF 神經網絡預測模型的R2、RMSE、MAPE、相對誤差在±15%內占比分別為0.72 和0.84、0.42 和0.32、12.7% 和10.5%、69.6% 和74.6%。與經驗公式預測模型相比，各項評價指標的評價優度有較大提升。說明機器學習方法對活動層厚度的預測更加準確、誤差更小。且可以看出RBF 神經網絡預測模型的各項評價指標評價優度均好于隨機森林預測模型。說明RBF 神經網絡方法在處理非線性關系、捕捉各活動層厚度擬合參數之間的特征聯系以及全局逼近和逼近精度都好于隨機森林方法。

表3 三種預測模型的各項評價指標Table 3 The evaluation indicators of the three prediction models

2.2 敏感性分析

由以上三種預測模型的預測結果對比可知RBF 神經網絡模型為預測活動層厚度的最佳模型，因此選擇對RBF 神經網絡模型的預測結果進行敏感性分析。敏感性分析是研究輸入參數的變化對輸出結果的影響程度。Sabol 方法［30］是一種全局敏感性分析方法，能夠有效的處理非線性響應和度量非加性系統中相互作用的影響。其計算公式如下：

式中：STi為總效應指數，表示敏感性程度大小；Var表示方差；E表示期望；Y表示輸出變量；Xi代表輸入因子；X～i表示除Xi所有變量的集合。計算結果如表4 所示，由表可知，含冰量對于活動層厚度的影響性最大，其次為年平均地表溫度、高程、等效緯度、緯度和平均植被指數。

表4 各預測因子敏感性分析Table 4 Sensitivity of each predictor

2.3 以往文獻結果對比

表5為以往學者對活動層厚度預測模型和本文RBF 神經網絡預測模型的評價指標對比。由表可知，本文以RBF 神經網絡建立的活動層厚度預測模型的預測效果好于以往的研究成果，R2和RMSE 的提升幅度較為明顯。此外，可以發現文獻中經驗公式預測方法的預測效果并不理想，而采用機器學習方法后，預測效果得到極大的改善，這進一步說明活動層厚度與各預測因子之間具有極強的非線性關系。從預測因子來看，本文與以往文獻相比除了溫度、植被、地形等因子，還考慮了含冰量的影響。由表可知，未考慮含冰量參數的RBF 神經網絡預測模型的R2和RMSE分別為0.71和0.43，與以往的機器學習方法研究結果精度相當。而考慮含冰量參數的RBF 神經網絡預測模型R2和RMSE 分別為0.84 和0.32，預測效果明顯提升，說明含冰量是活動層厚度的重要影響因素之一。

表5 不同活動層厚度預測模型的對比Table 5 Comparison of different active layer thickness prediction models

3 活動層厚度分布

通過以上模型的對比分析，最終確定以RBF 神經網絡作為預測方法建立活動層厚度預測模型。本文采用ArcGIS 10.5 軟件繪制活動層厚度分布圖，具體步驟為：（1）對年平均地表溫度、平均植被指數、等效緯度、高程和含冰量數據的柵格圖進行點數據提取，緯度數據由ArcGIS軟件計算得到；（2）將點數據輸入至RBF 神經網絡活動層厚度預測模型中，求取各點的活動層厚度；（3）將各點活動層厚度導入至ArcGIS 中，為使圖層精度提高、平滑性增強，采用克里金插值法進行制圖，即可獲得研究區內多年凍土活動層厚度分布區劃圖。

圖6 為青藏工程走廊活動層厚度分布區劃圖。由以上步驟（2）中的預測模型計算結果可知，研究區活動層厚度分布范圍在0.60～6.30 m，平均活動層厚度為3.55 m。圖7 為研究區各活動層厚度面積與面積占比。結合圖6 可知，研究區內活動層厚度主要為2～4 m，總面積為5 468.3 km2，面積占比為47.27%，主要分布于楚瑪爾平原至北麓河盆地和唐古拉山區南部至頭二九山區；活動層厚度大于4 m 次之，總面積為3 382.3 km2，面積占比為29.24%，整體分布偏向南部地區，主要分布于布曲河谷地至頭二九山區；活動層厚度為0～2 m 在研究區分布較少，面積占比僅達到12.2%，在研究區內分布較零散，主要分布于楚瑪爾河平原和可可西里山區。

圖6 研究區活動層厚度分布Fig.6 Active layer thickness distribution in the study area

圖7 各活動層厚度面積與面積占比Fig.7 Thickness area and area ratio of each active layer

各區域不同活動層厚度面積如圖8所示。由圖可知，昆侖山區至北麓河盆地的活動層厚度主要在1 m以上，整個青藏工程走廊內活動層厚度為1～2 m區域主要分布于此區段。活動層厚度為2～4 m 的區域也在此區段內占有較大面積比例，如楚瑪爾河平原和北麓河盆地的面積占比達到69.16% 和68.52%。風火山區至開心嶺山區的活動層厚度主要在2 m 以上，且此區段各區域內活動層厚度主要為2～4 m。其中，尺曲谷地和開心嶺山區活動層厚度大于4 m 的也占有較大面積。通天河盆地至頭二九山區的活動層厚度主要以2～4 m 和4 m 以上的情況居多，主要分布于布曲河谷地、唐古拉山和頭二九山地區。

圖8 各區域不同活動層厚度面積Fig.8 Areas of different active layer thicknesses in each region

4 活動層厚度與含冰量、地溫關系

為探究研究區活動層厚度與含冰量和地溫的關系，對研究區四類含冰量和不同區間的地溫的活動層厚度數據進行統計，得到其概率分布如圖9 所示。由圖9（a）可知，四類含冰量的活動層厚度主要分布區間（分布概率大于10%）為2.21～4.00 m 與4.42～5.03 m、1.62～3.60 m、1.33～3.10 m 和0.75～2.52 m，均值為3.48 m、2.55 m、2.28 m 和1.78 m。四類含冰量活動層厚度概率曲線隨著含冰量增加，分布明顯左偏，說明活動層厚度隨著土層含冰量增加而減小。由圖9（b）可知，不同區間的地溫的活動層厚度主要分布區間（分布概率大于10%）為2.08～3.97 m 與4.62～4.94 m、1.42～3.34 m、1.00～3.00 m 和0.75～2.34 m，均值為3.44 m、2.41 m、1.92 m 和1.48 m，即隨著地溫溫度升高，活動層厚度增加。這是由于活動層的地溫較低時，活動層內含冰量會相應的增加，而當外界溫度變化時，具有高含冰量的活動層內將會發生著大量的冰-水相變過程，從而導致活動層升溫速率過慢，致使活動層厚度較淺。

圖9 不同含冰量、地溫的活動層厚度頻率分布Fig.9 Frequency distribution of active layer thickness with different ice content and ground temperature

5 結論

通過利用青藏工程走廊監測斷面實測數據建立活動層厚度預測模型，再結合高精度遙感數據繪制青藏工程走廊多年凍土區的活動層厚度分布圖，得出以下結論：

（1）經驗公式、隨機森林和RBF 神經網絡預測模型的預測效果對比可知，基于RBF 神經網絡方法的活動層厚度預測模型預測效果最為精確、穩定，其預測值基本位于±15%誤差線內，且活動層厚度較大區間預測結果更接近實測值。

（2）考慮含冰量參數的RBF 神經網絡預測模型R2和RMSE 分別為0.84 和0.32，預測效果明顯提升，說明含冰量是活動層厚度的重要影響因素之一。

（3）研究區內活動層厚度主要為2～4 m，總面積為5 468.3 km2，面積占比為47.27%，主要分布于楚瑪爾平原至北麓河盆地和唐古拉山區南部至頭二九山區；活動層厚度大于4 m 次之，總面積為3 382.3 km2，面積占比為29.24%，整體分布偏向南部地區，主要分布于布曲河谷地至頭二九山區；活動層厚度為0～2 m 在研究區分布較少，面積占比僅達到12.2%，在研究區內分布較零散。

（4）活動層厚度隨含冰量增加而減小、隨地溫升高而增加，四類含冰量的活動層厚度主要分布區間為2.21～4.00 m 與4.42～5.03 m、1.62～3.60 m、1.33～3.10 m 和0.75～2.52 m，均值為3.48 m、2.55 m、2.28 m和1.78 m。