气候变化对黄土高原浅层滑坡影响的模拟研究
——以延安宝塔区为例

徐增辉, 金继明,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的年平均累积滑坡面积的相关性分析可知，未来滑坡出现的频次增多是累积滑坡面积增大的主导因素。