氣候變化對黃土高原淺層滑坡影響的模擬研究
——以延安寶塔區為例

徐增輝, 金繼明,2, 蔡耀輝, 楊 濤

(1.西北農林科技大學 水利與建筑工程學院, 陜西 楊凌 712100; 2.中國旱區節水研究院, 陜西 楊凌 712100;3.中國科學院 水利部 水土保持研究所, 陜西 楊凌712100; 4.西北農林科技大學 水土保持研究所, 陜西 楊凌 712100)

黃土高原由于特殊的地形、地貌和地質構造條件，淺層滑坡發生頻繁、分布廣泛，給當地人民的生命財產安全帶來了嚴重威脅[1-3]。寧奎斌等根據對陜西省2000—2016年地質災害時空分布規律及變化趨勢的研究發現黃土高原地區滑坡災害呈增加趨勢[4]。近年來采用模型對滑坡進行研究的數量逐年遞增。TRIGRS模型是由美國地質調查局開發的基于柵格的降雨誘發型斜坡穩定性計算模型[5]。模型主要由入滲模型，水文模型和斜坡穩定性模型組成，可計算降雨入滲引起的瞬時孔隙水壓力的變化，分析土壤在飽和非飽和條件下的邊坡穩定性，屬于瞬態水文—邊坡模型。夏蒙等利用TRIGRS模型成功模擬了山西興縣的黃土滑坡[6]，探究了不同坡度對區域斜坡破壞概率的影響。莊建琦等對比了TRIGRS和SINMAP滑坡預報模型對延安市寶塔區2013年群發滑坡的模擬結果，得出TRIGRS的模擬效果優于SINMAP的結論[7]。2016年Baum等改進開發了新的TRIGRS并行計算版本，可適用于大范圍多次模擬[8]。因此采用TRIGRS模型對黃土高原滑坡進行模擬是可行的。

在全球氣候變化的大背景下，區域性和局地性極端天氣增多，強降雨、干旱和洪澇災害頻發[9-11]。探究未來氣候變化對滑坡分布的影響，對政府制定防災減災政策，推廣群測群防體系有重要的現實意義。國外學者Ciabatta等在2016年將GCM模型與PRESSCA滑坡預警系統結合，得出意大利Umbria區域在2070—2099年間將增加40%以上[12];Massimiliano等2018年利用WRF區域氣候模型與TRIGRS結合，研究得出意大利區域未來滑坡增多主要受降雨事件時長變化影響[13]。但是目前國內學者對于黃土高原淺層滑坡未來變化的預報研究鮮有涉獵。滑坡事件和降雨密切相關，雨季降雨頻次、降雨強度和持續時長的改變，是否會對黃土高原這樣的生態環境脆弱區域，淺層滑坡的發生頻率和分布產生影響亟待研究。

本研究以黃土高原中部地質災害頻發的延安寶塔區為例，結合參與耦合模式比較計劃第五階段(CMIP5)中的3種模式兩種RCP情景和TRIGRS滑坡預報模型，利用統計降尺度的方法和Rosenblueth點估法對延安寶塔區1979—2100年滑坡分布情況進行模擬，研究不同氣候模式和情景下降雨事件變化趨勢，進而對淺層滑坡的影響。本研究可為黃土高原地區未來災害預警等提供參考。

1 試驗材料與方法

1.1 研究區概況

延安寶塔區(36°35′N,109°29′E)位于陜西黃土高原中部，屬典型繼承型和繼承—侵蝕混合型的黃土梁峁溝壑區，地表破碎，溝壑縱橫。中生界三疊系和侏羅系，新生界新近系上新統、全系統底層自下而上分布，地表更新統風積黃土(馬蘭黃土為主)廣泛分布，土質疏松，垂直節理裂隙發育，是黃土地質災害高易發區。地形特點西高東低，平均海拔800～1 400 m。該區屬典型的暖溫帶與中溫帶過渡區的西北干旱氣候，多年平均降雨量507.7 mm，降雨主要集中在6—9月，約占全年降雨的70%。據統計6—10月為滑坡等地質災害的主要發生期。

1.2 降雨數據與處理方法

1.2.1 歷史降雨數據和降尺度降雨數據 本研究選取CMIP5中3種升溫強度不同的GCM(GFDL-ESM2G，MIROC5，IPSL-CM5A-MR)2016—2100年的降雨預報數據作為源降雨數據(粗空間分辨率(1.3°～2.5°))。RCP假設選取RCP4.5和RCP8.5兩種情景作為對比。利用統計降尺度的方法將原始模型輸出的數據轉換為同ITPCAS數據集相同的0.1°分辨率。

選取中國區域地面氣象要素數據集(ITPCAS)[14]中1979—2015年時間分辨率為3 h，水平空間分辨率為0.1°的網格降雨數據作為歷史降雨數據進行模型驗證。

1.2.2 降雨事件提取方法 降雨事件是滑坡發生的最直接誘因，從長時間序列的氣象數據中提取一場可能誘發滑坡的降雨事件對于滑坡模擬預報起決定性作用，降雨事件的確定則依賴于經驗降雨閾值(累計降雨值)的定義[15-16]，以及分離兩場降雨事件中間的無雨期(T)的分離間隔時長選取[17]。本研究在前期研究中測試了間隔時長12 h，24 h，36 h，48 h等4種情景，并最終選取12 h作為最終模擬的分離間隔時長，同時選取累計降雨值30 mm為閾值對提取的降雨事件進行篩選[18]。將整個研究時間序列劃分為3個時期:歷史時期H(1979—2018年)，兩個未來時期F1(2019—2058年),F2(2059—2098年)。

1.3 淺層滑坡TRIGRS模型模擬

1.3.1 TRIGRS模型設置 TRIGRS模型針對飽和或近飽和土壤初始條件下的入滲模塊是基于Richard方程的線性解析解形式建立，可得到不同土層深度和時間下孔隙水壓力的大小:

ψ(Z,t)=[Z-d]β+

(1)

式中:ψ為孔隙水壓力;t為時間;Z為土層深度;d為地下水位深度;δ為坡角;Ks為垂向飽和滲透系數;β=cos2δ-(IZLT/Ks);IZLT為初始表面通量;InZ為第n個時間段的表面通量的強度;D1=D0/cos2δ;D0為飽和水力擴散系數;N為時間段總數;H(t-tn)為Heaviside階梯函數;ierfc(η)為高斯補誤差函數。

研究區內降雨誘發型滑坡的發生一般都有前期降雨存在，故采用土壤初始條件飽和的無限邊坡模型選項。TRIGRS模擬地形控制參數選取25 m分辨率的數字高程模型(DEM)，相關土壤參數參考莊建琦等[7]在延安寶塔區進行的實地試驗，參數選取見表1。

表1 TRIGRS模型土壤參數輸入

1.3.2 邊坡穩定性的計算 邊坡穩定性計算采用適用于淺層滑坡的無限邊坡穩定性分析，安全系數Fs(Factor of Safety)定義為土塊下滑阻力與動力的比值，基于Mohr-Coulomb破壞準則和孔隙水壓力的變化，不同土層在不同時刻的安全系數可以表示為:

(2)

式中:c′為土的有效黏聚力;φ′為土的有效內摩擦角;γw為地下水的容重;γs為土的容重;Fs為安全系數。

斜坡穩定性分析中存在較大的不確定性，建立合理的可靠度模型是解決不確定性的重要方法[19-20]。本研究采用Rosenblueth在1975年提出的點估法來建立淺層滑坡分析的可靠度分析模型，其基本原理為:在隨機變量分布未知時，可利用各種變量的均值和方差來計算出功能函數的一階矩和二階矩，從而得到斜坡的可靠度指標和破壞概率。斜坡穩定的功能函數可設為:

F=g(X1,X2,…,Xn)

(3)

式中:X1～Xn為與斜坡穩定性相關的多種因素;F為穩定系數。設F服從正態分布,且各項因素之間相互獨立,則相應的可靠度指標為:

(4)

式中:β為可靠度指標;EF和DF分別為穩定系數的均值和方差。

利用點估法，在隨機變量分布函數未知的情況下，在區間(xmin,xmax)上對稱地取2個點，如取均值的正負標準差，即:

(5)

考慮斜坡穩定性影響最大的粘聚力(c)和內摩擦角(φ)兩個因素，則可以得到4種不同的強度參數組合。通過不同強度參數的組合可以得到4種不同的穩定系數:

(6)

式中:σ為標準差。

根據式(4)和式(6)，斜坡的失穩概率Pf為:

Pf=1-Φ(β)

(7)

式中:Φ為標準正態分布函數。

采用上述的Rosenblueth點估法，在3個不同的時間段(H,F1,F2)將提取的降雨事件作為TRIGRS模型的降雨驅動參數分別模擬每場降雨的滑坡發生分布情況，其中每場降雨事件根據選取的點估法的4種參數組合需重復進行四次模擬可得出每個網格的四組安全系數Fs，再根據式(7)計算出每個網格點滑坡發生的概率。

對1.2.2中選取的全部降雨事件進行模擬，其中在歷史時期(H)共進行了3 888次模擬，兩個未來時期(F1,F2)共進行了9 300次模擬。一般來說，0.5的失效概率相當于安全系數為1的情況[21]，因此本研究后續分析中著重關注滑坡失穩概率大于0.5的情況。

2 結果與分析

2.1 TRIGRS模擬驗證

選取2013年7月研究區內發生的百年不遇的持續性強降雨誘發的大量淺層滑坡[22]來驗證TRIGRS模型的模擬結果。圖1為模擬的研究區(41.59 km2)滑坡分布情況，可以看出模擬失穩概率為0.5～1.0的模擬區域可以較好包含滑坡觀測點。

圖1 研究區內滑坡模擬分布

受試者工作特征(ROC)曲線是評估滑坡模擬結果的常用方法[7]，取失穩概率大于0.5的區域利用ROC曲線法對模擬結果進行評估。AUC的值越高說明模擬的結果越好，當AUC的值大于0.5時則說明模擬結果具有統計意義。圖2為研究區滑坡模擬結果的ROC曲線圖，曲線下的面積(AUC)為0.71，說明采用TRIGRS模型可以較好的模擬研究區淺層滑坡分布情況。

圖2 2013年7月滑坡模擬結果的ROC曲線

2.2 不同氣候模式和RCP情景下降雨事件變化趨勢

圖3為3種GCM模型產生的源數據、研究區降雨觀測數據和統計降尺度數據的降雨年平均對比情況，可以看出3種模式的降尺度結果與模型的源數據相比均有較大提高。

圖3 降雨降尺度數據年平均時間序列

3種不同的升溫強度氣候模式的降雨變化表現出不同的趨勢。MIROC5模型的降尺度數據降雨增大的趨勢最明顯，在RCP4.5和RCP8.5的情景下，在1979—2100年擬合趨勢線，其上升斜率分別為4.26，4.22，在F2時期的的年平均降雨量與H時期相比，分別增大了73.54%和62.07%。GFDL-ESM2G模型的降尺度數據年降雨量絕對值為3個模型中最大，但上升趨勢低于MIROC5，兩種RCP情景下全時期趨勢線上升的斜率分別為4.02，2.839，在F2和H兩個時期降雨分別增大58.37%和59.87%。IPSL-CM5A-MR模型在RCP4.5的情景下降雨呈微弱上升趨勢，全時間序列上升趨勢的斜率為0.94,F2和H兩個時期降雨增多15.32%，而在RCP8.5的情景下降雨呈下降趨勢，全時間序列降雨下降趨勢斜率為0.22,F2和H兩個時期相比降雨減少3.05%。因此不同的氣候模型對研究區降雨的未來預報表現出不同的趨勢，但是大多數情景下研究區未來降雨均呈現增多的趨勢。

提取的降雨事件累計降雨量的統計結果表明除IPSL-CM5A-MR模型外，H和F2兩個時期相比，未來降雨事件累計降雨量分布均有增大的趨勢。其中MIROC5模式下表現得變化趨勢最明顯，如在RCP8.5的情景下，累計降雨量大于90 mm的降雨事件出現的頻次提高了約10%。不同GCM和RCP情景假設下在不同時期可能誘發滑坡的降雨事件頻數見表2，其中3種不同模式的降雨事件頻數差異較大，GFDL-ESM2G升溫強度最高相對提取的降雨事件的數量也最多，而同一種模式下不同的RCP情景假設提取的降雨事件頻數規律不明顯。

表2 3種氣候模式兩種RCP情景不同時期降雨事件頻數

2.3 未來黃土高原淺層滑坡變化趨勢

在3種不同的全球氣候模式驅動下研究區內的滑坡分布情況出現了不同的趨勢，采用累積滑坡面積對研究區內TRIGRS滑坡模擬的分布情況進行評價，累積滑坡面積為滑坡次數與每次滑坡面積的乘積，其中包含多次重復滑坡的區域面積。圖4為研究區每十年年平均累積滑坡面積時間序列圖。3種模式中，MIROC5模型下累積滑坡面積的變化趨勢最大，在RCP8.5的情景下H時期的年平均累積滑坡面積由7.10 km2上升到F2時期的10.45 km2，增加了47.17%,RCP4.5情景下，H時期的年平均累積滑坡面積由6.82 km2上升到F2時期的9.77 km2，增加了43.16%。GFDL-ESM2G模型的淺層滑坡面積增加趨勢略小于MIROC5，在RCP4.5和RCP8.5的情景下F2時期的年平均累積滑坡面積分別增大23.1%和31.14%。IPSL-CM5A-MR模型滑坡面積基本沒有上升的趨勢，反而在RCP8.5的情景下，由于降雨事件的頻數和每場降雨事件的累計降雨量均略有減少，F2時期的年平均累積滑坡面積比H時期減少了23.5%。

圖4 每十年年平均累積滑坡面積的時間序列

圖5為研究區每十年滑坡發生頻數統計的時間序列圖。圖中的時間序列表明研究區內淺層滑坡發生的頻次在MIROC5和GFDL-ESM2G兩種模型下均呈現明顯的上升趨勢，滑坡發生頻次由21世紀初年均4～5場左右升高到21世紀末每年6～7次的頻率，IPSL-CM5A-MR模型滑坡發生頻次的變化趨勢最小，在RCP8.5的情景下有微弱下降的趨勢。

圖6為研究區單場滑坡面積每10 a的箱線圖，在時間尺度上，單場滑坡面積的變化趨勢并不明顯，單場滑坡面積的中位數處在一個相對穩定的狀態，每十年出現的最大滑坡事件的面積呈略微上升的趨勢，這是由于氣候變化導致研究區內降雨事件累計降雨量的極值變大導致的。由圖5和圖6統計分析可以得出，每十年的年平均累積滑坡面積變化與每十年滑坡發生的次數與有很高的相關性(R2>0.9)，而單場滑坡面積變化與每十年的年平均累積滑坡面積變化的相關系數R2僅有0.6左右。因此與單場滑坡面積增大相比未來滑坡出現的頻次增多是累積滑坡面積增大的主導因素。

圖5 每十年滑坡發生頻數時間序列

圖6 單場滑坡面積每十年箱線圖序列

4 結 論

(1) 采用全球氣候模式的統計降尺度降雨數據驅動TRIGRS滑坡模型，同時利用Rosenbluth點估法解決淺層滑坡可靠度分析問題，可以較好的模擬黃土高原淺層滑坡。

(2) 采用3種氣候模式進行滑坡模擬結果表現出不同的趨勢。MIROC5模型在未來時期F1(2019—2058年)滑坡增多的趨勢最明顯，到F2時期滑坡比歷史時期H(1979—2015年)年增加了大約45%。GFDL-ESM2G模型同時期滑坡面積為3種模型中最大，但滑坡面積增大的趨勢在H和F2時間段中小于MIROC5約為27%。IPSL-CM5A-MR滑坡面積變化總體呈減小趨勢，但減小的趨勢較小，F2時期的累計滑坡面積比H時期僅減小了11%左右。而同一種模型不同的RCP情景下，未來研究區滑坡分布變化表現出相同趨勢。

(3) 通過對每10 a滑坡發生的次數和單場滑坡面積變化與每10 a的年平均累積滑坡面積的相關性分析可知，未來滑坡出現的頻次增多是累積滑坡面積增大的主導因素。