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

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

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

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

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

1 試驗(yàn)材料與方法

1.1 研究區(qū)概況

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

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

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

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

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

1.3 淺層滑坡TRIGRS模型模擬

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

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

(1)

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

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

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

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

(2)

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

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

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

(3)

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

(4)

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

利用點(diǎn)估法，在隨機(jī)變量分布函數(shù)未知的情況下，在區(qū)間(xmin,xmax)上對(duì)稱(chēng)地取2個(gè)點(diǎn)，如取均值的正負(fù)標(biāo)準(zhǔn)差，即:

(5)

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

(6)

式中:σ為標(biāo)準(zhǔn)差。

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

Pf=1-Φ(β)

(7)

式中:Φ為標(biāo)準(zhǔn)正態(tài)分布函數(shù)。

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

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

2 結(jié)果與分析

2.1 TRIGRS模擬驗(yàn)證

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

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

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

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

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

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

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

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

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

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

2.3 未來(lái)黃土高原淺層滑坡變化趨勢(shì)

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

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

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

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

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

圖6 單場(chǎng)滑坡面積每十年箱線(xiàn)圖序列

4 結(jié) 論

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

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

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