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

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

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

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

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

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

1 試驗材料與方法

1.1 研究區(qū)概況

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

1.2 降雨數(shù)據(jù)與處理方法

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

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

1.2.2 降雨事件提取方法 降雨事件是滑坡發(fā)生的最直接誘因，從長時間序列的氣象數(shù)據(jù)中提取一場可能誘發(fā)滑坡的降雨事件對于滑坡模擬預報起決定性作用，降雨事件的確定則依賴于經(jīng)驗降雨閾值(累計降雨值)的定義[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為垂向飽和滲透系數(shù);β=cos2δ-(IZLT/Ks);IZLT為初始表面通量;InZ為第n個時間段的表面通量的強度;D1=D0/cos2δ;D0為飽和水力擴散系數(shù);N為時間段總數(shù);H(t-tn)為Heaviside階梯函數(shù);ierfc(η)為高斯補誤差函數(shù)。

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

表1 TRIGRS模型土壤參數(shù)輸入

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

(2)

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

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

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

(3)

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

(4)

式中:β為可靠度指標;EF和DF分別為穩(wěn)定系數(shù)的均值和方差。

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

(5)

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

(6)

式中:σ為標準差。

根據(jù)式(4)和式(6)，斜坡的失穩(wěn)概率Pf為:

Pf=1-Φ(β)

(7)

式中:Φ為標準正態(tài)分布函數(shù)。

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

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

2 結(jié)果與分析

2.1 TRIGRS模擬驗證

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

圖1 研究區(qū)內(nèi)滑坡模擬分布

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

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

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

圖3為3種GCM模型產(chǎn)生的源數(shù)據(jù)、研究區(qū)降雨觀測數(shù)據(jù)和統(tǒng)計降尺度數(shù)據(jù)的降雨年平均對比情況，可以看出3種模式的降尺度結(jié)果與模型的源數(shù)據(jù)相比均有較大提高。

圖3 降雨降尺度數(shù)據(jù)年平均時間序列

3種不同的升溫強度氣候模式的降雨變化表現(xiàn)出不同的趨勢。MIROC5模型的降尺度數(shù)據(jù)降雨增大的趨勢最明顯，在RCP4.5和RCP8.5的情景下，在1979—2100年擬合趨勢線，其上升斜率分別為4.26，4.22，在F2時期的的年平均降雨量與H時期相比，分別增大了73.54%和62.07%。GFDL-ESM2G模型的降尺度數(shù)據(jù)年降雨量絕對值為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%。因此不同的氣候模型對研究區(qū)降雨的未來預報表現(xiàn)出不同的趨勢，但是大多數(shù)情景下研究區(qū)未來降雨均呈現(xiàn)增多的趨勢。

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

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

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

在3種不同的全球氣候模式驅(qū)動下研究區(qū)內(nèi)的滑坡分布情況出現(xiàn)了不同的趨勢，采用累積滑坡面積對研究區(qū)內(nèi)TRIGRS滑坡模擬的分布情況進行評價，累積滑坡面積為滑坡次數(shù)與每次滑坡面積的乘積，其中包含多次重復滑坡的區(qū)域面積。圖4為研究區(qū)每十年年平均累積滑坡面積時間序列圖。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的情景下，由于降雨事件的頻數(shù)和每場降雨事件的累計降雨量均略有減少，F(xiàn)2時期的年平均累積滑坡面積比H時期減少了23.5%。

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

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

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

圖5 每十年滑坡發(fā)生頻數(shù)時間序列

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

4 結(jié) 論

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

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

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