一种滑坡灾害预测方法及装置与流程

本发明涉及地质灾害风险评估领域，特别是一种滑坡灾害预测方法及装置。

背景技术：

降雨对滑坡的影响很大，由于雨水的大量下渗，导致斜坡上的土石层饱和，甚至在斜坡下部的隔水层上积水，从而增加了滑体的重量，降低土石层的抗剪强度，导致滑坡产生。滑坡若不能及时得到预防和控制，不仅会造成交通堵塞，威胁人员的生命财产安全，还会引起重大次生危害。

考虑到雨水下渗时，土体或者岩体的异质性可能对入渗过程产生显著影响；土体在有些地方是分层的，例如图1中的边坡就是一个典型的具有分层土体的边坡，而不同土层中的水力学参数可能会显著不同，因此建立一个适用于分层土体的入渗分析方法非常必要；此外，由于实际降雨过程中的降雨强度很少是不变的，因此也需要建立一个适用于任意降雨条件的入渗分析方法，而现有基于入渗分析的区域性滑坡预测模型不能同时满足这两个要求，且很少有能评估滑坡体运动的模型。

对于传统滑坡预测方法而已，主要有两类：一类是基于经验性的降雨阀值，当雨量超过降雨阀值，发出预警信息；该类系统能预测滑坡、落石、泥石流灾害在某一降雨条件下是否会发生，但不能预测灾害发生的时间、地点及影响范围。另一类则是基于监测地质变形情况进行滑坡预警，一般监测边坡的变形量和运动加速度，将其经验性的变形阀值作为预警标准；该类系统在区域范围内对分布密集的灾害源的识别能力较弱，且大规模安装嵌入式传感器花费昂贵，不利于普及。传统滑坡预测方法考虑因素单调，预测结果不准确，易发生低判、错判和漏判，不利于灾害风险管理。

技术实现要素：

本发明的目的在于：针对现有技术中滑坡预测方法考虑因素不全面，导致预测结果不准确的问题，提供一种滑坡灾害预测方法及装置，能够较全面反映监测区域在时间和空间上对降雨过程的响应，准确预测滑坡发生的可能性及影响范围，为风险评价提供依据。

为了实现上述目的，本发明采用的技术方案为：

一种滑坡灾害预测方法，包括以下步骤：

步骤一，建立数字地形模型，将研究区域划分为网格，每一个网格视为一个单元，赋予每个单元所需的地质信息、空间雨量分布信息等水文信息；

步骤二，在任意降雨条件下，对边坡进行入渗分析，获取孔隙水压力；

步骤三，基于孔隙水压力，对于浅层滑坡，用无限边坡模型计算安全系数fs，然后求取失稳概率pf，以判断边坡稳定性和可靠性；其中，浅层滑坡为一种滑坡类型，本领域技术人员将滑坡按滑坡体的厚度划分为浅层滑坡、中层滑坡、深层滑坡和超深层滑坡。

步骤四，收集大量历史滑坡数据，通过回归分析得到用于计算滑坡运动距离平均值的经验关系；该经验关系应该要考虑到滑坡体的体积效应，将基于滑坡运动距离的经验公式作为滑坡运动的休止条件，当滑坡体达到最大的运动距离后，滑坡体将停止运动；

步骤五，将在步骤三中至少一边或者一角相连的不稳定单元归为一组，一组视为一个滑坡，确认滑坡的位置和体积；结合步骤四中基于历史滑坡数据建立的经验关系，将其作为滑坡运动休止条件，评估滑坡体的运动轨迹和堆积位置，进而预测滑坡影响范围。

其中，在入渗分析时，利用理查兹方程(richards方程)描述边坡中的入渗过程：

式中k是渗透系数；ψ是孔隙水压力；z^*为土体边坡法向方向上的位移；βs是边坡坡度；t是时间；w是体积含水量；

利用指数型关系描述渗透方程和土水特征曲线：

其中ks是饱和渗透系数；θr和θs分别是残留体积含水量和饱和体积含水量；αi是材料参数；为得到闭合的解析解，材料参数αi在上下两层土体取相同的值；

进而通过拉普拉斯变换(laplace变换)、留数定理及一系列数学变换，得到恒定降雨条件下kb和kt的解；其中下标b和t分别表示下层和上层土体，k表示渗透系数k与饱和渗透系数ks的比值，kb＝k/ksb，kt＝k/kst；

结合单位阶跃函数和叠加原理，进一步得到在任意降雨条件下的kb和kt、孔隙水压力及地下水位通量函数，其中孔隙水压力为：

其中z也是自定义的无量纲数：

-lb≤z*≤0时，z＝αib(cosβs)²z*；0≤z*≤lt时，z＝αit(cosβs)²z*；

l是土体边坡垂直方向的土体厚度。

当到达地下水的入渗水量超过了最大的重力排水量，地下水位会相应上升。对每个单元逐个进行入渗分析，可以得知研究区域对降雨在时间和空间上的响应。

进一步地，对于浅层滑坡，滑坡体的厚度相对其平面尺寸来说较小，因此可以用无限边坡模型计算安全系数fs，并将安全系数fs作为稳定性评价指标。在步骤三计算安全系数的过程中，由于在理查兹方程(richards方程)的解析解中陡峭的湿润锋通常不存在，因此将湿润锋定义为在体积含水量发生显著变化的地方，从地表到湿润锋对fs进行计算，找到最小fs和相应的位置；如果fs小于1，认为滑坡会发生。

孔隙水压力ψ(z,t)为正值时，fs用如下公式计算：

其中c”(kpa)和φ'(°)分别是土体的有效凝聚力和摩擦角；γw(kn/m³)和γs(kn/m³)分别是水和土的单位容重，z'是沿垂直方向的位移。

孔隙水压力ψ(z,t)为负值时，fs用基于扩展的莫尔库仑破坏准则计算：

其中φb(°)是由基质吸力对剪切强度贡献的等效摩擦角；φb设为φ'的三分之二。

进一步地，为了考虑输入参数的不确定性，为风险评估打下基础，同时进行可靠性分析，求取失稳概率pf作为可靠性分析指标，失稳概率pf大于0.5时，视为滑坡将会发生。室外实验和室内实验研究发现c'和tanφ'通常呈负相关，为了避免负值出现，假设c'和tanφ'服从对数正态分布。对用蒙特卡罗(montecarlo)方法产生的fs进行卡方拟合优度检验，结果表明fs服从对数正态分布时，fs的标准差可表示为：

其中xi是一个随机变量，即c'，tanφ'或tanφb；nv是随机变量的数量；ρ'是随机变量之间的相关系数；v(e)是误差方差，用来描述和实际情况的差异，其服从正态分布。

进一步通过下式求取失稳概率pf：pf＝p(lnfs＜0)＝φ(-λ/ξ)，

其中，

λ＝ln[e(fs)]-0.5ξ²

其中e(fs)是fs的期望值，即滑动面的fs值。

进一步地，在对滑坡运动分析时，滑坡体实际上可能会从一个单元同时运动到几个相邻的较低位置的单元，本发明不考虑滑坡体在运动过程中的发散，也不考虑滑坡体的惯性效应和弹跳效应；假设滑坡体会沿着最陡峭的路径运动到相邻较低的位置，可用以下公式计算一个单元和相邻单元之间的坡度：

对于东西南北方向，坡度

对于对角方向，坡度

其中δh是两个单元之间的高程差；ac是单元的宽度。

本发明基于上述方法步骤提供一种滑坡预测装置，该装置包括至少一个处理器，以及与处理器通信连接的存储器；存储器存储有可被至少一个处理器执行的指令，指令被处理器执行，以使处理器能够执行滑坡预测的方法。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明综合考虑影响滑坡发生的各个因素，能全面得知研究区域在时间和空间上对降雨过程的响应，进而可以预测研究区域范围内的滑坡发生情况。

2、本发明中的入渗分析方法适用于任意降雨条件和双层土体边坡，更贴合复杂的地质条件和实际雨量分布情形，大大增强了本发明的实用性。

3、本发明能够很好地实时预测降雨型滑坡，准确度达到92％，为灾后重建方案选择提供客观科学依据，也可对可能发生滑坡区域范围进行综合检测预警、处治或绕避等安全措施，避免二次伤害，保证生命财产安全。

附图说明

图1是典型分层土体边坡的结构示意图。

图2是滑坡灾害预测模型构架示意图。

图3是简化的双层土体边坡示意图。

图4是任意降雨过程示意图。

图5是不同滑坡的log(h/l)和logv图。

图6是研究单元的运动方向示意图。

图7(a)是模型计算的滑坡分布示意图(t＝5小时)。

图7(b)是模型计算的滑坡分布示意图(t＝15小时)。

图7(c)是模型计算的滑坡分布示意图(t＝36小时)。

图7(d)是模型计算的滑坡分布示意图(t＝46小时)。

具体实施方式

下面结合附图，对本发明作详细的说明。

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例1

一种滑坡灾害预测方法，构架滑坡灾害预测模型，如图2，包括以下步骤：

一、获取数据

建立数字地形模型，将研究区域划分为网格，每一个网格视为一个单元，赋予每个单元相应所需的地质信息、空间雨量分布信息等水文信息，如高程、坡度、地质分类、降雨信息、土体厚度、地下水水位和土体参数。

建立数字地形模型时，可通过地形测绘仪器进行实地野外测绘获得密集高程点数据，在此基础上建立dem(digitalelevationmodel，数字高程模型)；或利用航空摄影测量手段，通过航空影像立体相对提取dem；或利用现有的纸质地形图，通过软件数字化地形图矢量方式dem。

二、灾害分析

1、入渗分析模型

基于获取的空间雨量分布，利用入渗分析获取每个单元的孔隙水压力和基质吸力等水文条件，进而得知研究区域在时间和空间上对降雨过程的响应，好进一步预测滑坡发生的区域范围。

在本实施例中将研究区域内所有的边坡都简化为如图3所示的双层土体边坡。假设每层土是均质的：

利用理查兹(richards)方程描述边坡中的入渗过程：

式中k是渗透系数；ψ是孔隙水压力；z^*为土体边坡法向方向上的位移；βs是边坡坡度；t是时间；w是体积含水量。

利用指数型关系描述渗透方程和土水特征曲线：

其中ks是饱和渗透系数；θr和θs分别是残留体积含水量和饱和体积含水量；αi是材料参数。土体的单位容重随含水量的变化而变化，为得到闭合的解析解，材料参数αi在上下两层土体取相同的值，则上述方程(1)-(3)联合则可转变为以下形式：

设初始地下水位在下层土体的下边界。在降雨过程开始之前的长时段内，入渗率qa设为恒定值，地下水位以上的土体处于非饱和状态。在降雨过程中，入渗率是qb(恒定降雨条件)或者qbn(任意降雨条件)。

为了求解式(4)，首先定义几个无量纲数用于后面的推导：

-lb≤z*≤0时：z＝αib(cosβs)²z*，lb＝αib(cosβs)²lb；(5)

0≤z*≤lt时，z＝αit(cosβs)²z*，lt＝αit(cosβs)²lt；(6)

kb＝k/ksb,qab＝qa/ksb,qbb＝qb/ksb,qbbn＝qbn/ksb(7)

kt＝k/kst,qat＝qa/kst,qbt＝qb/kst,qbtn＝qbn/kst(8)

其中下标b和t分别表示下层和上层土体，l是沿垂直方向的土体厚度。

进而公式(4)可表示为：

其中，

下层土体的初始条件为

上层土体的初始条件为

下层土体下边界的孔隙水压力为

上层土体上边界的边界通量为

上下两层分界面的通量为

上下两层分界面的孔隙水压力为

利用拉普拉斯(laplace)变换和留数定理可得到恒定降雨条件下kb和kt的解：

其中ψ0是地下水位的孔隙水压力；rab，rbb，rcb，rat，rbt和rct在推导过程中产生的关于是z和t的函数。

在实际中，一场暴雨过程中的降雨强度会如图4所示随时间变化。小时降雨强度广泛应用于实际情况，如果降雨历时有n小时，那么就有n个间隔。

利用单位阶跃函数来考虑任意降雨条件。利用单位阶跃函数和叠加原理，在任意降雨条件下的kb和kt可表示为

其中h(t-tn)是单位阶跃函数；nt是降雨历时的小时数；tn是间隔n开始的时刻；qbbn和qbtn分别是qbb和qbt在时间间隔n期间的值。

孔隙水压力可用如下公式得到：

地下水位的通量为

当到达地下水的入渗水量超过了最大的重力排水量，地下水位会相应上升。对每个单元逐个进行入渗分析，可以得知研究区域对降雨在时间和空间上的响应。此滑坡灾害预测模型的入渗分析方法适用于任意降雨条件和双层土体边坡，实用性强。

2、边坡稳定性和可靠性评估模型

基于入渗分析，孔隙水压力的上升或者基质吸力的下降会导致安全系数降低，容易发生滑坡；而且研究区域面积很大，区域内单元数量众多(超过180,000个)，必须进一步采用切实可行的方法来评价大量边坡的稳定性和可靠性。

对于浅层滑坡，滑坡体的厚度相对其平面尺寸来说较小，因此可以用无限边坡模型计算安全系数fs，将安全系数作为稳定性评价指标。由于在理查兹(richards)方程的解析解中陡峭的湿润锋通常不存在，因此将湿润锋定义为在体积含水量发生显著变化的地方，从地表到湿润锋对fs进行计算，找到最小fs和相应的位置；如果fs小于1，认为滑坡会发生。

如果孔隙水压力ψ(z,t)为正值，fs用如下公式计算：

其中c”(kpa)和φ'(°)分别是土体的有效凝聚力和摩擦角；γw(kn/m³)和γs(kn/m³)分别是水和土的单位容重，z'是沿垂直方向的位移。

如果孔隙水压力ψ(z,t)为负值，fs用基于扩展的莫尔库仑破坏准则计算：

其中φb(°)是由基质吸力对剪切强度贡献的等效摩擦角；φb设为φ'的三分之二。

进一步地，为了考虑输入参数的不确定性，为风险评估打下基础，同时进行可靠性分析，求取失稳概率pf作为可靠性分析指标，失稳概率pf大于0.5时，视为滑坡将会发生。室外实验和室内实验研究发现c'和tanφ'通常呈负相关。为了避免负值出现，假设c'和tanφ'服从对数正态分布。对用蒙特卡罗(montecarlo)方法产生的fs进行卡方拟合优度检验，结果表明fs服从对数正态分布时，fs的标准差可表示为：

其中xi是一个随机变量，即c'，tanφ'或tanφb；nv是随机变量的数量；ρ'是随机变量之间的相关系数，本实施例中取值-0.5；v(e)是误差方差，用来描述和实际情况的差异，其服从正态分布。

通过下式求取失稳概率pf：pf＝p(lnfs＜0)＝φ(-λ/ξ)，

其中，

λ＝ln[e(fs)]-0.5ξ²，

其中e(fs)是fs的期望值，即滑动面的fs值。

3、滑坡运动分析模型

基于三维数字高程模型可以进行滑坡运动分析，三维数字高程模型包含平面坐标和相应的高程数据。首先结合入渗分析方法确认滑坡的位置和体积，然后基于历史滑坡数据建立关于滑坡运动距离的经验关系；滑坡体沿着山体最陡峭的方向运动，将最大滑坡运动距离作为运动休止条件，可进而获知滑坡体的运动轨迹和堆积位置，预测滑坡影响范围，这对滑坡风险评估非常重要。

a.建立滑坡运动距离经验公式

收集2008年汶川地震后降雨诱发的20个浅层滑坡数据及地震诱发的31个滑坡数据，包括滑坡运动的高差h、运动距离l和滑坡体的体积v。用这51个滑坡数据进行回归分析，通过回归分析得到的用于计算滑坡运动距离平均值的经验关系如下：logl＝0.023+0.023logv+logh，在回归计算时，pearson相关系数r²取值为0.88。方程可写成到达角的形式：log(h/l)＝-0.023-0.023logv。该经验关系考虑了滑坡体的体积效应，可以作为滑坡运动的休止条件，当滑坡体达到最大的运动距离后，滑坡体将停止运动。

以上51个滑坡数据和得到的经验关系如图5中所示。

b.滑坡运动分析

将有至少一边或者一角相连的不稳定单元归为一组，一组视为一个滑坡。在确定每个不稳定单元的滑坡体体积之后，通过相加，可以确定每组总的滑坡体体积。在预测滑坡运行轨迹时，滑坡体的运动受地形控制，在运动过程中滑坡的方向可能会因地形而改变，滑坡体可能会从一个单元同时运动到几个相邻的较低位置的单元，但本实施例不考虑滑坡体在运动过程中的发散，也不考虑滑坡体的惯性效应和弹跳效应。假设滑坡体会沿着最陡峭的路径运动到相邻较低的位置，如图6所示，每个单元有8个可能的运动方向，滑坡体会沿着最陡峭的路径运动到相邻的较低单元。一个单元和相邻单元之间的坡度可用以下公式计算：

对于东西南北方向，坡度

对于对角方向，坡度

其中δh是两个单元之间的高程差；ac是单元的宽度。

滑坡体沿着最陡峭的路径运动，当达到最大的运动距离后，滑坡体将停止运动。

4、模型验证

2010年“8.13”暴雨席卷了映秀及其附近地区，降雨历时长达46小时，诱发了大量滑坡。本实施例采用此次暴雨诱发的滑坡验证模型的正确性。如图7(a)-(d)是模型计算的滑坡分布。下表1是预测滑坡和实际情况对比的混淆矩阵，可从中获知滑坡预测模型的准确度达到92％，能够很好地实时预测降雨型滑坡，为灾后重建方案选择提供客观科学依据，对可能发生滑坡区域范围进行综合检测预警、处治或绕避等安全措施，避免二次伤害，保证生命财产安全。

表1混淆矩阵

注：(a)定义(tp：真阳性；fn：真阴性；fp：假阳性；fn：假阴性；n：单元总数(182417)；acc：预测精度)。

(b)降雨结束时基于安全系数的困惑矩阵。如果安全系数小于1，滑坡灾害预测模型中滑坡将发生。

(c)降雨结束时基于失稳概率的困惑矩阵。如果失稳概率大于0.5，滑坡灾害预测模型中滑坡将发生。

实施例2

本发明基于实施例1中的预测方法提供一种滑坡预测装置，该装置包括至少一个处理器，以及与处理器通信连接的存储器；存储器存储有至少一个可被对应处理器执行的指令，指令被处理器执行，以使处理器能够执行实施例1所述全部或部分方法步骤。

实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的指令程序可以存储于计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤。当本发明上述集成的单元以软件功能单元的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实施例的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机、服务器、或者网络设备等)执行本发明各个实施例所述方法的全部或部分。而前述的存储介质包括：移动存储设备、rom、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。