劉志云, 黃 川, 于 暉, 鐘振濤, 崔福慶
(1.長安大學地質工程與測繪學院,陜西西安 710054; 2.中交第一公路勘察設計研究院有限公司西安中交公路巖土工程有限責任公司,陜西西安 710075)
青藏高原被稱為地球的“第三極”,其內(nèi)發(fā)育著世界上中、低緯度帶海拔最高、面積最廣的凍土區(qū)[1]。隨著全球氣候變暖,青藏高原多年凍土產(chǎn)生嚴重的退化現(xiàn)象,主要表現(xiàn)為地溫升高、活動層厚度增加以及多年凍結層消失等[2-4]?;顒訉邮侵概救诨浼緝鼋Y覆蓋于多年凍土之上的土層,是凍土地層內(nèi)水熱交換最主要的區(qū)域[5],受經(jīng)緯度、高程、植被覆蓋度、地表溫度、土壤性質及氣候環(huán)境等諸多因素的影響。活動層厚度的變化將會對寒區(qū)水文、地質、環(huán)境和工程建筑產(chǎn)生一系列影響。因此,建立活動層厚度的預測模型以及研究活動層厚度的分布特征對青藏高原地區(qū)工程構筑物的設計、建造及后期養(yǎng)護具有重要意義。
活動層厚度的現(xiàn)場測量主要有鉆孔、觸探、挖探等機械方法以及監(jiān)測站監(jiān)測等方法。Hon 等[6]利用動態(tài)圓錐貫入儀在活動層與多年凍土層的動力錐穿透指數(shù)的不同,從而測定活動層厚度。更多學者則是通過不同地區(qū)的現(xiàn)場測量數(shù)據(jù)對活動層厚度的年際變化進行研究分析[7-9],Zhao 等[10]根據(jù)青藏高原鉆孔溫度曲線證明,在1967—1997 年期間,其活動層的厚度以平均速度為0.71 cm·a-1增長。Wu 等[11]通過對青藏公路沿線監(jiān)測站監(jiān)測數(shù)據(jù)分析研究,發(fā)現(xiàn)從1995—2007 年期間,青藏高原多年凍土區(qū)活動層厚度正以平均7.5 cm·a-1的速率持續(xù)增加?;顒訉雍穸痊F(xiàn)場測量雖然能夠獲得較高的精確度,但存在著成本高、數(shù)據(jù)樣本少且難以連續(xù)觀測致使其無法刻畫大區(qū)域特征等缺點。
針對以上不足,許多學者發(fā)現(xiàn)經(jīng)驗公式模型、統(tǒng)計模型和數(shù)值模型對活動層厚度大區(qū)域分布研究具有較好的效果。Pang 等[12-13]用Stefan 和Kudryavtsev 公式計算了青藏高原活動層厚度的變化,給出了青藏高原活動層厚度空間分布圖并預測了2049年、2099年的活動層厚度變化情況。Zhao等[14]將Kudryavtsev 公式與Lund-Postam-Jena 模型[15-16]耦合進一步提高了經(jīng)驗公式模型對活動層厚度的預測能力。Ni 等[17]采用統(tǒng)計模型與機器學習結合的方式對青藏高原活動層厚度進行預測,研究發(fā)現(xiàn)模型比以往的研究具有較高的精度。Zhang 等[18]利用改進的GIPL2 模型對青藏高原無人區(qū)活動層厚度進行數(shù)值模擬,并基于實測活動層厚度數(shù)據(jù)驗證了數(shù)值模型的準確性和優(yōu)越性??傮w而言,經(jīng)驗公式、統(tǒng)計模型和數(shù)值模型能夠有效的描繪大區(qū)域活動層厚度分布情況,但對局地因素考慮不全且未能結合實地勘查數(shù)據(jù)進行研究,致使空間分辨率低,難以反映實際情況。
近年來,隨著機器學習方法和遙感技術的快速發(fā)展,目前,諸多學者采用對現(xiàn)場勘查數(shù)據(jù)進行預測模型建立,再結合大范圍遙感數(shù)據(jù)進行區(qū)域分布模擬的方式已經(jīng)廣泛應用于青藏高原凍土的相關研究,包括凍土分布[19-20]、凍土地溫預測[21-22]、凍土滑坡敏感性[23]等方面。鑒于此,本文通過青藏工程走廊沿線300 組活動層厚度鉆孔監(jiān)測數(shù)據(jù),基于年平均地表溫度、平均植被指數(shù)、等效緯度、緯度、高程和含冰量等參數(shù)建立了活動層厚度的經(jīng)驗公式、隨機森林和RBF 神經(jīng)網(wǎng)絡預測模型,通過對比三種預測模型的預測效果,并結合高精度遙感數(shù)據(jù),運用預測精度最高的活動層厚度預測模型繪制了青藏工程走廊多年凍土區(qū)段沿線活動層厚度分布圖。
青藏工程走廊始于格爾木,止于拉薩市,橫穿青藏高原1 120 多公里,穿越多年凍土區(qū)約550 km,是內(nèi)陸進入西藏的重要通道。本文以走廊內(nèi)多年凍土區(qū)段(西大灘—安多)為研究區(qū),研究活動層厚度分布狀況。如圖1所示,研究區(qū)以青藏公路、青藏鐵路為基準線向兩側外延10 km,全長約540 km,地理坐標位于32°~36°N、91°~95°E,海拔介于3 716~6 191 m,該區(qū)地貌類型豐富,包括中高山區(qū)、高平原、低山丘陵、河谷等。
圖1 研究區(qū)監(jiān)測點分布Fig.1 Distribution of monitoring points in the study
活動層是地-氣間水熱交換的主要場所,大的氣候背景決定了大區(qū)域活動層厚度的宏觀分布狀況。但在一定條件下,局地因素的影響將超過氣候的影響,導致區(qū)域內(nèi)相同氣候背景下活動層厚度的分布異常。局地因素主要通過影響太陽輻射、熱對流和熱傳導等過程從而影響活動層厚度的大小。因此,本文擬選用年平均地表溫度、平均植被指數(shù)、緯度、高程、等效緯度和含冰量六類數(shù)據(jù)作為預測模型建立的預測因子。
年平均地表溫度、平均植被指數(shù)和高程遙感數(shù)據(jù)由美國國家航天航空局下載的地表溫度(MOD11A2 H25V05,2000—2016 年)、植被指數(shù)(MOD13Q1 H25V05,2000—2016 年)和SRTMDEM(Shuttle Radar Topography Mission-Digital Elevation Model)數(shù)據(jù)產(chǎn)品中提取。等效緯度是表征太陽輻射對地表的影響,也是判斷坡面走向的重要數(shù)據(jù),可通過由SRTM-DEM 數(shù)據(jù)產(chǎn)品中提取的坡度、坡向和緯度數(shù)據(jù)計算得到:
式中:φ′為等效緯度;l為坡度;h為坡向;φ為緯度。
凍土含冰量是多年凍土的基本特征指標之一,且活動層內(nèi)不同含冰量也對活動層熱量吸收能力有一定影響。含冰量數(shù)據(jù)由現(xiàn)場鉆取多年凍土上限以下8 m 深度內(nèi)的土層,再經(jīng)過現(xiàn)場觀察記錄和室內(nèi)試驗綜合確定。凍土含冰量鉆孔點沿青藏鐵路和青藏公路間隔2 km 布置,確定所有鉆孔點含冰量賦存狀態(tài)后,再利用概率插值得到走廊帶內(nèi)凍土含冰量分布,如圖2所示。目前,含冰量分類標準主要以《凍土工程地質勘察規(guī)范》(50324—2014)為主[24],分為少冰、多冰、富冰、飽冰和含土冰層五類。但在青藏高原地區(qū),多年凍土上限以下土層內(nèi)含冰量變化較為劇烈。因此,本文將根據(jù)多年凍土上限以下8 m 深度內(nèi)土層含冰量主要的賦存類型將含冰量劃分為少冰-多冰、多冰-富冰、富冰-飽冰、飽冰-含土冰層四類。
圖2 研究區(qū)含冰量分布Fig.2 Distribution of ice content in the study area
用于建立預測模型的活動層厚度實測數(shù)據(jù)來源于青藏公路、青藏鐵路沿線監(jiān)測斷面地溫監(jiān)測數(shù)據(jù)(監(jiān)測點如圖1),監(jiān)測工作由中交第一公路勘察設計研究院有限公司高寒高海拔地區(qū)道路工程安全與健康國家重點實驗室格爾木觀測基地完成,數(shù)據(jù)時限為2006—2016 年,數(shù)據(jù)采集采用測溫法測得,精度為±0.05 ℃?;顒訉雍穸缺O(jiān)測點共計300組,在研究區(qū)內(nèi)分布較為平均,基本體現(xiàn)了青藏高原工程走廊帶活動層厚度特征,具有較好的代表性。
表1是各預測因子之間的相關性分析及共線性分析結果。由表可知,部分預測因子間雖然表現(xiàn)出較顯著的相關性,但線性相關關系較弱,說明各預測因子間雖然能夠相互影響,但影響程度較小。并且從共線性分析可以看出,各預測因子與活動層厚度之間的容差皆大于0.2 且膨脹方差因子(VIF)皆小于5,這進一步說明各預測因子之間不存在明顯的共線性。
表1 各因子之間的相關性及共線性分析Table 1 Correlation and collinearity analysis among factors
1.4.1 經(jīng)驗公式
最小二乘法是通過找尋數(shù)據(jù)誤差平方和的最小值,從而確定數(shù)據(jù)的最佳匹配函數(shù)。在本文中,含冰量數(shù)據(jù)作為分類變量,因此,在進行經(jīng)驗公式擬合之前需對含冰量數(shù)據(jù)進行編碼處理。目前,對于分類變量的編碼方式主要有獨熱編碼、虛擬編碼和效應編碼等。其中,效應編碼具有不冗余、易于解釋的優(yōu)點,其原理是在一個具有n個類別的變量中選取一個類別作為參照(賦值為-1),從而創(chuàng)建n-1 個指標變量(賦值為0 或1)。含冰量數(shù)據(jù)編碼處理方法如表2所示。
表2 含冰量數(shù)據(jù)效應編碼Table 2 Effect coding of ice content data
以上編碼代表當含冰量為少冰-多冰時,X1、X2、X3取值為-1、-1、-1;含冰量為多冰-富冰時,X1、X2、X3取值為1、0、0,以此類推。因此基于最小二乘法擬合的活動層厚度經(jīng)驗公式如式(2)所示:
式中:h為活動層厚度;Ts為年平均地表溫度;N為平均植被指數(shù);φ為緯度;H為高程;φ′為等效緯度。
1.4.2 機器學習
機器學習是現(xiàn)代智能技術中的一種重要方法,其主要研究從樣本數(shù)據(jù)中尋找規(guī)律,并根據(jù)這些規(guī)律預測未來或無法觀測的數(shù)據(jù)[25]。機器學習方法已經(jīng)在青藏高原凍土相關研究中得以廣泛的應用并取得了較好的成果[26-28]。其中,隨機森林和徑向基函數(shù)神經(jīng)網(wǎng)絡在處理復雜的非線性問題時具有良好的抗噪能力,因此,本文將選用隨機森林和RBF神經(jīng)網(wǎng)絡建立活動層厚度預測模型。
隨機森林[29]是通過對多個決策樹(每棵樹都擬合到訓練數(shù)據(jù)的Bootstrap 樣本)求平均值來擬合組合模型,每棵樹中的每個拆分都考慮隨機的預測變量子集,通過這種方式,將多個弱模型組合起來生成更為強大的模型,其結構如圖3 所示。本文使用該方法建立了活動層厚度預測模型,并將數(shù)據(jù)集按8:2 比例進行隨機分塊,分別作為訓練集和驗證集,樣本個數(shù)分別為240和60。隨機森林預測模型中決策樹的數(shù)量為100 個,每棵樹的最小拆分大小為10個,最大拆分大小為2 000個。
圖3 隨機森林結構圖Fig.3 Random forest structure diagram
RBF 神經(jīng)網(wǎng)絡是三層前饋神經(jīng)網(wǎng)絡,由輸入層、隱含層和輸出層三部分構成。其中,輸入層用于接收外部信息,隱含層實現(xiàn)參數(shù)間的非線性轉換,輸出層用于輸出最終結果,其結構如圖4 所示。
圖4 神經(jīng)網(wǎng)絡結構圖Fig.4 Neural network structure diagram
k-fold 交叉驗證法是一種能夠有效的利用小樣本數(shù)據(jù)得到最優(yōu)模型的方法。其原理是將原始數(shù)據(jù)分成k組,將每一組數(shù)據(jù)分別做一次驗證集,其余的k-1 組數(shù)據(jù)作為訓練集,從而得到k個模型,最終選擇最優(yōu)的模型。本文的RBF 神經(jīng)網(wǎng)絡模型的結構設置為年平均地表溫度、平均植被指數(shù)、等效緯度、緯度、高程和含冰量作為輸入層,活動層厚度作為輸出層,高斯徑向基函數(shù)為激活函數(shù)。隱含層設置為兩層,每層12 個節(jié)點。k-fold 交叉驗證數(shù)為5,其中4 份用于模型訓練,1 份用于模型驗證,模型訓練和模型驗證樣本數(shù)分別為240和60。
1.4.3 評價指標
活動層厚度預測模型評價指標采用決定系數(shù)(R-Squared,R2)、平均絕對百分比誤差(Mean Absolute Percentage Error,MAPE)、均方根誤差(Root Mean Square Error,RMSE)和相對誤差±15%內(nèi)占比。其中R2、MAPE、RMSE 表達式分別見式(3)~(5)。
圖5 為三種活動層厚度預測模型的預測結果。由圖可知,經(jīng)驗公式方法擬合效果最差,預測值大部分都分布于±15%誤差線外,并且可以發(fā)現(xiàn)在低活動層厚度區(qū)間預測值偏大,高活動層厚度區(qū)間預測值偏小,整體分布極為離散。隨機森林方法擬合效果次之,預測值大部分落在±15%誤差線內(nèi)。但整體來看,樣本點仍較為離散,部分測點偏離實測值較大。說明隨機森林方法雖然能夠較好的預測活動層厚度,但存在一定誤差。RBF 神經(jīng)網(wǎng)絡方法擬合效果最好,預測值基本位于±15%誤差線內(nèi),并且整體更靠近實測值線,尤其高活動層厚度區(qū)間。因此,從三種方法預測結果對比可知,基于RBF 神經(jīng)網(wǎng)絡方法的活動層厚度預測模型預測效果最為精確、穩(wěn)定。
圖5 各活動層厚度預測模型結果對比Fig.5 Comparison of results of prediction models for the thickness of each active layer
表3為三種預測模型的各項評價指標。由表可知,經(jīng)驗公式預測模型的各項評價指標的評價優(yōu)度均較差,R2僅為0.40,RMSE 和MAPE 達到0.75 和24.5%,相對誤差在±15%以內(nèi)的占比也不足一半。隨機森林和RBF 神經(jīng)網(wǎng)絡預測模型的R2、RMSE、MAPE、相對誤差在±15%內(nèi)占比分別為0.72 和0.84、0.42 和0.32、12.7% 和10.5%、69.6% 和74.6%。與經(jīng)驗公式預測模型相比,各項評價指標的評價優(yōu)度有較大提升。說明機器學習方法對活動層厚度的預測更加準確、誤差更小。且可以看出RBF 神經(jīng)網(wǎng)絡預測模型的各項評價指標評價優(yōu)度均好于隨機森林預測模型。說明RBF 神經(jīng)網(wǎng)絡方法在處理非線性關系、捕捉各活動層厚度擬合參數(shù)之間的特征聯(lián)系以及全局逼近和逼近精度都好于隨機森林方法。
表3 三種預測模型的各項評價指標Table 3 The evaluation indicators of the three prediction models
由以上三種預測模型的預測結果對比可知RBF 神經(jīng)網(wǎng)絡模型為預測活動層厚度的最佳模型,因此選擇對RBF 神經(jīng)網(wǎng)絡模型的預測結果進行敏感性分析。敏感性分析是研究輸入?yún)?shù)的變化對輸出結果的影響程度。Sabol 方法[30]是一種全局敏感性分析方法,能夠有效的處理非線性響應和度量非加性系統(tǒng)中相互作用的影響。其計算公式如下:
式中:STi為總效應指數(shù),表示敏感性程度大??;Var表示方差;E表示期望;Y表示輸出變量;Xi代表輸入因子;X~i表示除Xi所有變量的集合。計算結果如表4 所示,由表可知,含冰量對于活動層厚度的影響性最大,其次為年平均地表溫度、高程、等效緯度、緯度和平均植被指數(shù)。
表4 各預測因子敏感性分析Table 4 Sensitivity of each predictor
表5為以往學者對活動層厚度預測模型和本文RBF 神經(jīng)網(wǎng)絡預測模型的評價指標對比。由表可知,本文以RBF 神經(jīng)網(wǎng)絡建立的活動層厚度預測模型的預測效果好于以往的研究成果,R2和RMSE 的提升幅度較為明顯。此外,可以發(fā)現(xiàn)文獻中經(jīng)驗公式預測方法的預測效果并不理想,而采用機器學習方法后,預測效果得到極大的改善,這進一步說明活動層厚度與各預測因子之間具有極強的非線性關系。從預測因子來看,本文與以往文獻相比除了溫度、植被、地形等因子,還考慮了含冰量的影響。由表可知,未考慮含冰量參數(shù)的RBF 神經(jīng)網(wǎng)絡預測模型的R2和RMSE分別為0.71和0.43,與以往的機器學習方法研究結果精度相當。而考慮含冰量參數(shù)的RBF 神經(jīng)網(wǎng)絡預測模型R2和RMSE 分別為0.84 和0.32,預測效果明顯提升,說明含冰量是活動層厚度的重要影響因素之一。
表5 不同活動層厚度預測模型的對比Table 5 Comparison of different active layer thickness prediction models
通過以上模型的對比分析,最終確定以RBF 神經(jīng)網(wǎng)絡作為預測方法建立活動層厚度預測模型。本文采用ArcGIS 10.5 軟件繪制活動層厚度分布圖,具體步驟為:(1)對年平均地表溫度、平均植被指數(shù)、等效緯度、高程和含冰量數(shù)據(jù)的柵格圖進行點數(shù)據(jù)提取,緯度數(shù)據(jù)由ArcGIS軟件計算得到;(2)將點數(shù)據(jù)輸入至RBF 神經(jīng)網(wǎng)絡活動層厚度預測模型中,求取各點的活動層厚度;(3)將各點活動層厚度導入至ArcGIS 中,為使圖層精度提高、平滑性增強,采用克里金插值法進行制圖,即可獲得研究區(qū)內(nèi)多年凍土活動層厚度分布區(qū)劃圖。
圖6 為青藏工程走廊活動層厚度分布區(qū)劃圖。由以上步驟(2)中的預測模型計算結果可知,研究區(qū)活動層厚度分布范圍在0.60~6.30 m,平均活動層厚度為3.55 m。圖7 為研究區(qū)各活動層厚度面積與面積占比。結合圖6 可知,研究區(qū)內(nèi)活動層厚度主要為2~4 m,總面積為5 468.3 km2,面積占比為47.27%,主要分布于楚瑪爾平原至北麓河盆地和唐古拉山區(qū)南部至頭二九山區(qū);活動層厚度大于4 m 次之,總面積為3 382.3 km2,面積占比為29.24%,整體分布偏向南部地區(qū),主要分布于布曲河谷地至頭二九山區(qū);活動層厚度為0~2 m 在研究區(qū)分布較少,面積占比僅達到12.2%,在研究區(qū)內(nèi)分布較零散,主要分布于楚瑪爾河平原和可可西里山區(qū)。
圖6 研究區(qū)活動層厚度分布Fig.6 Active layer thickness distribution in the study area
圖7 各活動層厚度面積與面積占比Fig.7 Thickness area and area ratio of each active layer
各區(qū)域不同活動層厚度面積如圖8所示。由圖可知,昆侖山區(qū)至北麓河盆地的活動層厚度主要在1 m以上,整個青藏工程走廊內(nèi)活動層厚度為1~2 m區(qū)域主要分布于此區(qū)段?;顒訉雍穸葹?~4 m 的區(qū)域也在此區(qū)段內(nèi)占有較大面積比例,如楚瑪爾河平原和北麓河盆地的面積占比達到69.16% 和68.52%。風火山區(qū)至開心嶺山區(qū)的活動層厚度主要在2 m 以上,且此區(qū)段各區(qū)域內(nèi)活動層厚度主要為2~4 m。其中,尺曲谷地和開心嶺山區(qū)活動層厚度大于4 m 的也占有較大面積。通天河盆地至頭二九山區(qū)的活動層厚度主要以2~4 m 和4 m 以上的情況居多,主要分布于布曲河谷地、唐古拉山和頭二九山地區(qū)。
圖8 各區(qū)域不同活動層厚度面積Fig.8 Areas of different active layer thicknesses in each region
為探究研究區(qū)活動層厚度與含冰量和地溫的關系,對研究區(qū)四類含冰量和不同區(qū)間的地溫的活動層厚度數(shù)據(jù)進行統(tǒng)計,得到其概率分布如圖9 所示。由圖9(a)可知,四類含冰量的活動層厚度主要分布區(qū)間(分布概率大于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)可知,不同區(qū)間的地溫的活動層厚度主要分布區(qū)間(分布概率大于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,即隨著地溫溫度升高,活動層厚度增加。這是由于活動層的地溫較低時,活動層內(nèi)含冰量會相應的增加,而當外界溫度變化時,具有高含冰量的活動層內(nèi)將會發(fā)生著大量的冰-水相變過程,從而導致活動層升溫速率過慢,致使活動層厚度較淺。
圖9 不同含冰量、地溫的活動層厚度頻率分布Fig.9 Frequency distribution of active layer thickness with different ice content and ground temperature
通過利用青藏工程走廊監(jiān)測斷面實測數(shù)據(jù)建立活動層厚度預測模型,再結合高精度遙感數(shù)據(jù)繪制青藏工程走廊多年凍土區(qū)的活動層厚度分布圖,得出以下結論:
(1)經(jīng)驗公式、隨機森林和RBF 神經(jīng)網(wǎng)絡預測模型的預測效果對比可知,基于RBF 神經(jīng)網(wǎng)絡方法的活動層厚度預測模型預測效果最為精確、穩(wěn)定,其預測值基本位于±15%誤差線內(nèi),且活動層厚度較大區(qū)間預測結果更接近實測值。
(2)考慮含冰量參數(shù)的RBF 神經(jīng)網(wǎng)絡預測模型R2和RMSE 分別為0.84 和0.32,預測效果明顯提升,說明含冰量是活動層厚度的重要影響因素之一。
(3)研究區(qū)內(nèi)活動層厚度主要為2~4 m,總面積為5 468.3 km2,面積占比為47.27%,主要分布于楚瑪爾平原至北麓河盆地和唐古拉山區(qū)南部至頭二九山區(qū);活動層厚度大于4 m 次之,總面積為3 382.3 km2,面積占比為29.24%,整體分布偏向南部地區(qū),主要分布于布曲河谷地至頭二九山區(qū);活動層厚度為0~2 m 在研究區(qū)分布較少,面積占比僅達到12.2%,在研究區(qū)內(nèi)分布較零散。
(4)活動層厚度隨含冰量增加而減小、隨地溫升高而增加,四類含冰量的活動層厚度主要分布區(qū)間為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。