基于水土耦合物理机制的区域滑坡概率预报技术研究*

徐 辉, 张少杰, 黎 力, 刘敦龙

(1.国家气象中心，北京 100081;2.中国科学院、水利部成都山地灾害与环境研究所，四川 成都 610041;3.重庆市国土资源和房屋勘测规划院，重庆 401121;4.武汉大学 测绘学院，湖北 武汉 430079;5.成都信息工程大学 软件工程学院，四川 成都 610225)

近些年来,我国每年因滑坡灾害造成的经济损失巨大, 并导致重大的人员伤亡。滑坡预报作为一项非工程型措施，可以通过提前发布灾情的方式，避免群死群伤事件的发生，是地质灾害防灾减灾重要的手段之一。

降雨在很大程度上是滑坡发生的最主要激发因素。降雨型滑坡的预报预警可以通过确定诱发滑坡的降水阈值，并结合降水预报进行滑坡的预报预警。在此方面，国内外许多学者进行了大量的研究，Caine根据滑坡资料研究了降雨的强度、持续时间与滑坡灾害关系并最早提出了I-D降水阈值[1]。Cannon 等根据旧金山海湾滑坡和降雨数据建立了滑坡与降雨强度和持续时间临界关系曲线[2]。Brand[3]在 1984 年就对香港地区降雨与滑坡的关系进行了统计分析。鲍其云[4]确定了浙江62个丘陵山区县( 市、区)引发滑坡和泥石流的降雨强度-历时阈值，并提出了基于降雨强度-历时阈值和滑坡泥石流易发程度区划图的预报系统模型。陈君[5]对湖北省恩施市地质灾害进行详查，通过正态函数分析，确定不同性质滑坡的阈值区间，并建立起降雨诱发滑坡灾害的风险评价模型。

如上所述，目前国内外大部分的滑坡预警预报模型主要以降水阈值为基础的统计预报模型，需要依赖大量的历史数据。而且没有考虑下垫面对地质灾害发生所产生的作用，缺乏地质灾害机理框架的支撑，误报率和漏报率高。而基于物理机制的滑坡预报模型可以在一定程度上解决统计预报模型过于依赖降水因素的弊端[6]，也可以从机理层面上解释滑坡随其内部力学参数变化而演变的规律，因此，受到越来越多的学者关注。基于物理机制的滑坡预报模型是以物理模型分析和实验数据为基础，考虑降雨过程、入渗、径流等作用及其与地质环境因素的相互作用，模拟不同降雨条件下地质环境背景条件的变化，进行滑坡地质灾害背景条件的稳定性分析，从而实现数值模拟和降雨-地质灾害变化过程反演，继而开展基于物理机制的降雨型滑坡地质灾害预警预报[7-17]。

滑坡的水土耦合物理机制预报模型多是采用极限平衡法和数值分析法，通过计算区域内每个网格单元的安全系数来判定网格单元是否失稳破坏，在此基础上开展滑坡的确定性预报。然而，无论是极限平衡法还是数值分析法在计算的过程中所用到的岩土参数和计算得到的滑坡稳定性都是确定的数值，忽略了滑坡岩土参数的空间变异性与不确定性，也将滑坡的安全系数考虑的绝对化，这也导致在对某些滑坡进行稳定性评价时算出来的稳定系数存在偏差，进而影响滑坡的预报准确率，因此，有必要用不确定的分析方法来评价滑坡的安全系数[18-20]。

采用基于水土耦合物理机制的数值模拟进行降雨型滑坡地质灾害预警预报将是未来发展趋势，但实际上目前国内鲜有科研或业务人员在此方面进行研究并在业务中进行预报应用，而且大多只是针对滑坡灾害点或小尺度区域开展基于物理机制的滑坡预报和研究[21]。因此，本文在无限边坡物理模型的基础上采用蒙特卡洛MontE-Carlo随机抽样法提出了基于物理机制的滑坡概率预警模型，并首次尝试选取我国的西南地区(四川省、云南省、贵州省和重庆市)建立了基于水土耦合物理机制的区域滑坡概率模型，开展西南地区滑坡预报试验研究，以期为国家级地质灾害风险预警业务提供预报支撑。

1 基于水土耦合物理机制的滑坡概率预报模型

1.1 滑坡概率预报

以求解安全系数Fs为目标的无限边坡模型具有合理的理论假设(式(1))，在一定程度上可以反映出浅层滑坡的机理过程，具有较高的可靠性。所以，在保证模型准确的力学参数(土体粘结力和内摩擦角)输入变量的前提下，便可利用安全系数Fs对区域内每个网格单元的稳定性进行确定性预报。即当Fs≤1时，网格被判定为不稳定，当Fs>1时，网格被视为稳定。由此，便可实现基于物理机制的确定性滑坡预报模式。

(1)

式中：c为土体的粘结力，φ是土体的内摩擦角，ua是为大气压ua=0，φb与基质吸力相关，当基质吸力较低时，该值与内摩擦角φ接近[16]，Hs为土层厚度，ψ=(ua-uw)为基质吸力。

但无论是采用野外采样的方法[8]还是根据预报区内的土地质地分布[22]而获取的预报区内的每个网格单元的力学参数都存在较高的误差[9]。模型中所采用的表面看似确定的土力学参数实际上却存在不确定性[15]。这一客观存在的、且目前无法解决的事实最终会导致每个网格单元的安全系数FS存在不确定性。

物理参数的不确定性可以采用概率密度函数描述[9]，即假定土力学参数在一定的边界范围内服从均匀分布:c=U(cmin,cmax)，φ=U(φmin,φmax)。采用蒙特卡洛随机取值的思想，在每个滑坡预报时间段内，在U(cmin,cmax)和U(φmin,φmax)区间内随机取上述两个土力学参数n次，以此体现力学参数取值不确定性的问题，具体的随机取值方法是：两个土力学参数随机变量的取值以在U(0,1)上均匀分布的随机数(ri=U(0,1))为基础，利用公式(2)和公式(3)分别在其边界内随机取值。

ci=ri(cmax-cmin)+cmin；

(2)

φi=ri(φmax-φmin)+φmin。

(3)

式中：cmin和φmin是分别是土体粘结力和内摩擦角浮动区间的下边界；cmax和φmax分别是土体粘结力和内摩擦角浮动区间的上边界：

(4)

最后生成的矩阵Mi以一种整体的形式代表了参数取值的不确定性这一物理现象，位于区间内U(cmin,cmax)和U(φmin,φmax)的矩阵中的任意元素均具有物理意义。同时，Mi中的每一组[ci,φi]以数学语言的方式将每个网格单元的土力学参数取值的可能性给定量化表述出来，为建立具有不确定性属性的输入变量与输出结果之间的定量化关联，提供了有效的途径。

另外，由公式(2)和(3)可知，若要实现利用蒙特卡洛法则随机提取土体粘结力c和内摩擦角φ，需要确定上边界和下边界值。但是，确定这两个参数的上下边界值是极为困难的，就目前发表的发现资料而言，还尚未有明确的函数来描述。Raia 等[15]以试验确定的土力学参数 (corigin，φorigin) 为基准值，分别上边界和下边界分别按照1%，10%，和100%浮动，并利用大量的计算证明上下边界值调整幅度大的条件下，预报效果最佳。为此，采用了Raia 等[15]方法中的100%的幅值，来调整本文的土力学参数的上边界和下边界值。

crandom[0.5×corigin,2×corigin]；

(5)

φrandom[0.5×φorigin,2×φorigin]。

(6)

式中:crandom和φrandom为某次随机提取的土体粘结力和内摩擦角值，corigin和φorigin为利用岩石力学手册确定的土力学参数值。

由公式(1)可知，在其他参数已知的情况下，每一组[ci,φi]均会生成唯一对应的安全系数Fsi=[Fs1,Fs2,Fs3,…,Fsn]。安全系数数组Fsi代表了某个网格单元可能存在的n个不同的稳定性状态，而在这n种状态中，Fsi≤1出现的次数(SumFs<1)或比重值(P∈[0,1])则表示了则表示了网格单元倾向于失稳破坏的趋势的大小：

(7)

P值越大，表明在具有不确定性的输入变量(土力学参数)的作用下，同样具有不确定性的预报结果越倾向于滑坡发生。P值超过50%预示着网格单元已完全倾向于失稳破坏的一端。所以，比重P可以用来作为物理预报模型中具有不确定属性的输入量与输出结果之间的纽带。比重P以具体的数值定量化这两者之间的联系：即采用比重P来判定在具有不确定性输入参数的作用下，浅层滑坡发生可能性(或危险度)的大小。若是将P∈[0,1]离散成表1中一系列的参考区间，区间从1到5表示滑坡发生的概率值会逐步增加，危险等级也越来越高，并采用不同的预警颜色与预警等级相匹配。如此，便可利用比重P来建立滑坡发生概率(或危险等级)与输入参数不确定性之间的关联，建立起基于物理机制的滑坡概率预报模型。

表1 滑坡概率预报等级划分

1.2 水文过程模拟

由公式(1)可知，概率预警模型的运行需要关键水文参数的驱动，因此需要利用水文过程模拟为概率模型提供输入参数。降水作用于下垫面的水文过程是一个极为复杂的非线性过程。该过程通过降水截留、入渗和蒸散发等一系列的水文环节改变下垫面土体含水量的分布。降水过程中影响土体含水量分布的最为关键的水文环节是入渗过程。降水满足冠层的截留后，便会到达地面。就中国西南地区的降雨特点而言，区域遭遇雨季时的坡面土体一般处于非饱和状态。非饱和土在降水入渗作用下，其内部垂向的水分运动规律可采用一维非饱和Richards[23]水分运动微分方程描述：

(8)

式中：θ为土体的体积含水量(mm3/mm3)；D(θ)=K(θ)/(dθ/dψm)为非饱和土的扩散率(mm3/h)，ψm为非饱和土基质吸力(mm)；z为竖向坐标(mm)，以地表为原点，向下为正；K(θ)为土层的非饱和导水率(mm/h)。

入渗边界(上边界)控制：地表处于非饱和状态时，地表入渗能力强。若某时刻的降水强度小于地表的入渗能力时，降水全部渗入土体，无径流产生。此时的降水的入渗边界由微分方程(9)控制。

(9)

式中：R(t)为降水强度(mm/h)。而当降水强度超过地表的入渗能力后随即转化为径流流走，不考虑洼地积水的有压入渗，此时地表饱和，控制每个网格入渗边界的方程(7)可转化为θ=θs，θs为区域内的某种土壤类型对应的饱和体积含水率。

下边界控制：对于降水对土壤深度40 cm以下的土体含水量影响较小[24]，且地下水位很深[25]的区域水文特性而言，降水入渗作用对地下水位的变化影响较小。所以下边界不考虑下边界与地下水之间的水分交换，是较为合理的一种选择方式。

2 研究区域和数据

西南地区(包括四川、云南、贵州和重庆市)山地分布广泛，地质条件复杂，构造运动剧烈，地质灾害尤为突出，严重威胁着人民生命财产安全，给国民经济造成重大损失。因此本文选取该区域作为预报区域构建基于水土耦合物理机制的滑坡概率预报模型系统，并进行滑坡预报检验。

2.1 静态数据

本文提出的概率预报模型所需的基础数据可以分为两类：一是主要用于水文过程模拟的下垫面数据，包括水文过程模拟所需的数据包括：数字高程模型(DEM)，土地利用类型、土壤类型和植被指数等。土地利用、土壤类型、土层厚度由数据库获取，同样需要利用重采样技术获取相应的网格数据。

在兼顾计算效率和数据精度的基础上，西南地区的DEM栅格精度为500 m。在DEM的基础上，利用ArcGIS导出西南地区的坡度分布图。

土地利用(栅格精度250 m)、土壤类型(栅格精度1000 m)以及植被叶面指数(栅格精度500 m)等数据通过数据库FAO(http://www.fao.org/geonetwork/srv/en/main.home)数据库获取，需要利用重采样技术获取与DEM相一致的网格数据。

另外一类则主要是影响坡体稳定性的土力学参数(土体粘结力和内摩擦角)的获取。根据公式(2)和(3)，若要实现以蒙特-卡洛的思想在一定的区间内随机取值的方式来体现这两个输入参数的不确定性，必须首先设定一个浮动值的基点。U(cmin,cmax)和U(φmin,φmax)区间的上下浮动边界均与每个网格单元的所设定的这个基点值相关联。为此，我们根据预报区域的地质图，获取区域的岩性分布，并依据岩石力学参数手册对不同岩性的力学参数(土体粘结力c、内摩擦角φ)进行赋值，而后利用重采样技术获取相应的网格数据，以每个网格单元的力学参数值作为中心值，获取每个网格单元的土体粘结力和内摩擦角的浮动区间U(cmin,cmax)和U(φmin,φmax)。

图1 西南地区DEM及坡度值分布图

图2 西南地区土力学参数分布图

图3 2017年7月31日20时-8月1日20时西南地区滑坡预报结果图

2.2 降水数据

本文中采用前期有效含水量经验衰减公式确定区域内的初始含水量分布。在确定了滑坡预报初始时刻的土体含水量后，以国家气象中心提供的降水预报为输入，通过水文模型的数值计算，可以为区域内每个网格单元稳定性计算提供关键水文参数。

2.2.1 基于定量降水估测(QPE)的前期有效含水量处理

国家气象中心提供的定量估测降水(QPE)数据为1 km×1 km的栅格数据，数据范围覆盖整个中国，时间尺度为逐1 h。需要依据西南地区500 m的栅格精度，对QPE数据进行重采样以及裁剪等流程，提取与西南地区DEM精度一致的QPE数据。以预报日期为起点，前推15 d，并以天为单位累加其范围内的24 h QPE数据，生成前15 d的前期雨量数据。而后，利用经验衰减公式，确定区域内每个栅格的前期有效雨量值。

(10)

式中：AntR为每个栅格的前期雨量值，n预报当天前推15 d，K为衰减系数，取值0.84，Ri为某个栅格在前i天的降雨值。

2.2.2 基于定量预报降水(QPF)的预报降水处理

国家气象中心提供的定量预报降水(QPF)数据为5 km×5 km的栅格数据。QPF数据面向整个中国地区可以提供未来3 h、6 h、24 h和48 h等不同时间段的数据形式，系统针对QPF的数据形式，可人工设定对西南地区未来不同时间段的滑坡预报。同样的，依据西南地区500 m的栅格精度，对QPE数据进行重采样以及裁剪等流程，提取与西南地区DEM精度一致的QPE数据。

3 应用实例分析

据中国气象局灾情直报系统不完全统计，2017年汛期期间西南地区发生11起造成人员伤亡的滑坡案例(详见表2)。为说明上述方法的可靠性，对2017年西南地区发生的这11起滑坡事件的预报情况进行检验评估，从预报结果来看，该模型成功预报9次，预报地质灾害等级为3级以上，成功率达82%。具有较好的预报准确率。

表2 2017年西南地区滑坡预报情况一览表

为了进一步说明上述方法的可靠性，下文将选取2017年7月31日至8月1日的降水过程进行具体分析。2017年7月31日21：22至8月1日凌晨03：30，云南省普洱市镇沅县者东镇出现暴雨到大暴雨天气，每小时降雨量达52.6 mm，帮海村累积雨量达100.2 mm。在本轮降水过程作用下，云南省普洱市镇沅县者东镇出现大暴雨天气。强降水引发者东镇境内多处河流漫溢、山体滑坡等衍生地质灾害。灾害造成者东镇 15个村185个村民小组14 695人不同程度受灾，死亡3人。本模型对该时段的滑坡进行等级预报，从2017年7月31日20时-8月1日20时的24 h预报结果(图3)来看，显示7月31日夜间云南省普洱市镇沅县者东镇发生滑坡灾害的概率较高，并发出3级黄色预警和4级橙色预警，成功预报出此次滑坡灾害。

4 结论和讨论

本文在无限边坡物理模型的基础上，通过引入蒙特-卡洛随机取值的方法，构建了滑坡概率预报模型，并结合国家气象中心提供的定量降水估测(QPE)和定量降水预报(QPF)数据，搭建了西南地区滑坡概率预报模型系统，以西南地区2017年的滑坡个例进行预报验证，在一定程度上验证了该模型的应用精度和准确率，具有较好的业务应用前景。

(1)基于水土耦合物理机制的滑坡概率预报模型对西南地区的滑坡灾害有较高的预报命中率，在对2017年发生在西南地区的11起典型滑坡案例预报中，该模型成功预报9次，成功率达82%，仅漏报2例，漏报率为18%左右，具有较好的预报准确率。预期可为滑坡防灾减灾提供技术支撑。模型虽然在本文区域个例中预报应用效果较好，但还需要在更多的滑坡预报个例中进一步验证。

(2)驱动基于水土耦合物理机制的滑坡概率预报模型的前端数据主要是降水实况和预报数据，从预报案例结果来看，降水数据的精度是影响滑坡准确性的重要因素。相较于QPE和QPF的栅格精度，滑坡一般发生在较小的空间尺度，因此需进一步精细化降水产品的空间精度，这不但有助于提升局地突发强降雨的精度，也有助于预测小尺度的滑坡灾害。

滑坡成因机理复杂，且影响滑坡的因素众多，因此所有的滑坡预报模型均存在误报率高的问题，这是当前束缚区域滑坡技术的主要瓶颈。通过获得更加精细化的下垫面数据，进一步改进水文过程模拟能力，深入研究复杂地形下水土耦合的滑坡的形成机理和精度是下一步着重改进和研究的方向，从而进一步消除因数据精度造成的滑坡误报和漏报问题。