CN102306233B - 降雨作用下流域滑坡时空预测方法 - Google Patents

Abstract

本发明涉及流域尺度的滑坡时空预测方法，旨在提供一种降雨作用下流域滑坡时空预测方法。该方法包括：使用InHM水文模型、无限边坡稳定模型、滑坡计算模块计算流域内任意三维位置的安全系数SF，并通过可视化模块拟合为指定时间流域安全系数分布图和指定时间流域滑坡深度图，并显示于显示设备；若存在滑坡，则选择显示流域滑坡时间分布图。本发明可实现对流域滑坡进行快速有效的分析和预测；克服了以往模型在复杂地形需要洼地预处理的不足，保留了自然流域地貌水文响应特性，滑坡体更接近于实际滑坡体形状，降低预测误差。同时提高了模型的精度和适用性。适用于简单边坡、复杂边坡及流域尺度的滑坡预测。

Description

降雨作用下流域滑坡时空预测方法

技术领域

本项目涉及一种流域尺度的滑坡时空预测方法，通过将基于物理概念的水文模型和无限边坡稳定性模型相结合，对流域范围内滑坡灾变过程进行三维时空预测，是一种将土力学、水力学与计算机数值模拟相结合的滑坡预警预报方法。

背景技术

我国滑坡是地质灾害多发区，滑坡每年造成巨大的经济损失和人员伤亡。滑坡作为一种地质灾害，严重地危害当地居民的生产生活。滑坡的预警预报一直是工程界研究的重点之一，其分析理论越来越丰富，但由于滑坡灾害的突发特征及其在运动规律上所表现出的不确定性，有关滑坡预测方面的研究仍没有达到较高预测精度和通用性。滑坡的产生受多种自然因素以及人类活动的影响，其中降雨是诱发滑坡最重要也是最广泛的因素。降雨对滑坡的作用主要表现在，雨水的大量下渗，导致斜坡上的土石层含水量增大，地下水位线抬升，甚至在斜坡下部的隔水层上击水，从而增加了滑体的重量，降低土石层的抗剪强度。因此，如何准确地获取土石层含水量及地下水位线的时空变化，是求解滑坡体稳定性的基础。目前，滑坡的预测方法主要集中在简单边坡如道路边坡、水库边坡领域，在流域尺度边坡稳定性时空分布规律的影响相对较少。

根据滑坡体稳定性计算方法的不同，滑坡预测模型主要可分为确定性模型与非确定性模型两类。确定性模型通过力学极限平衡原理来求解不同位置不同深度下稳定性系数，应用最广的如无限边坡稳定性模型；非确定性模型主要通过数理统计方法对相关参数进行分析，主要包括指标分析模型、概率统计分析预测模型以及模糊集预测模型等。根据所采用的水文模型类型的不同，滑坡预测模型可分为基于经验概念水文模型和基于物理概念水文模型两大类。基于经验概念的水文模型相对简单，发展较早，以此基础建立的滑坡预测模型已经广泛应用。相比之下，基于物理概念的水文模型模拟信息量丰富，能捕捉多种水文要素的时空变化过程，但由于求解复杂，对计算要求较高，以此基础建立的滑坡预测模型发展相对较晚。

不同的滑坡体稳定计算方法与水文模型相结合，产生了多种类型的滑坡预测模型。这些模型均有各自的优缺点，如非确定性模型一般适用于特定的研究地区，很难提取与滑坡体地质、力学及其他外界条件有关的信息，不能进行很好的物理力学解释，现有的确定性模型主要缺点在于很难找到一个适合于具有多种不确定因素的复杂介质模型，其中现有的无限边坡稳定模型采用地表坡度方向作为滑坡体滑动方向，在洼地区域的滑坡稳定性计算趋于安全。传统的经验水文模型，在求解时需要对地形进行洼地填充预处理，以满足水文模型的正常产汇流以及滑坡方向等要求，在运用得当的前提下，能够提供具有相对精度的结果，适用于某些对过程细节没有要求的科学研究及工程实践。

发明内容

本发明要解决的技术问题是，克服现有技术中的不足，提供一种降雨作用下流域滑坡时空预测方法。

为解决技术问题，本发明的解决方案是：

提供一种降雨作用下流域滑坡时空预测方法，包括以下步骤：

(1)根据流域的DEM地形数据及土层厚度建立三维网格，向参数输入模块输入不同位置的土壤、植被、气象数据，生成InHM水文模型所需的边界条件；

(2)根据其边界条件计算流域内的水文要素时空变化，使用InHM水文模型模拟地表水二维流动和地下水三维流动(模型包括地表水二维流动模拟和地下水三维流动模拟，两者共同影响土壤含水量的变化)、任意位置土壤含水量变化和任意位置水力坡度变化的水文要素；

(3)改进的无限边坡稳定模型以水力坡度方向作为滑动面方向建立极限平衡方程，以水力坡度方向作为潜在的滑动方向，并将之用于滑坡计算模块；根据InHM水文模型生成的水文要素，使用滑坡计算模块计算流域内任意三维位置的安全系数SF；若安全系数SF小于1，则为潜在的滑坡位置；记录该位置的深度、时间和体积；

所述任意点的安全系数SF表示为：

$\mathrm{SF}=\frac{c^{\′}+c_r+{\cos }^2\mathrm{\θ}\tan \mathrm{\φ}{\∫}_0^z{\mathrm{\ρ}}^{\′}\mathrm{gdz}}{\sin \mathrm{\θ}\cos \mathrm{\θ}{\∫}_0^z\mathrm{\ρgdz}}---\left(1\right)$

式中：c′、c_r分别为土的有效粘聚力和根系粘聚力，kPa；θ为滑动方向，即水力坡度方向，°；φ为有效内摩擦角，°；ρ、ρ′分别为土壤密度和土壤有效密度，kg/m³；z为土体深度；

而土壤密度表示为含水量的函数，具体如下：

ρ＝n_sS_wρ_w+(1-n_s)d_sρ_w                       (2)

其中d_s为土壤颗粒比重；n_s为土壤孔隙度；S_w为土体含水量；ρ_w为水密度；

式(1)与(2)中，c′、c_r、φ、d_s和n_s由土壤特性数据给定，S_w通过InHM水文模型获得；

水力坡度方向θ通过线性平均方式获得，计算公式如下：

$\mathrm{\θ}=\arctan (\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{(\mathrm{\ψ}+h)}_i}{\mathrm{\Δ}{l}_i}/n)---\left(3\right)$

式中，Δ(ψ+h)_i表示该点与其相连的第i节点总水头差，Δl_i表示该点与其相连的第i节点距离，n表示与该点相连的节点数；

将式(2)和式(3)代入式(1)可获得边坡上任意点沿深度方向的安全系数；当SF＞1，该位置是稳定的，反之则不稳定，表示存在滑坡；

(4)通过可视化模块，将步骤(3)所得的数据拟合为指定时间流域安全系数分布图和指定时间流域滑坡深度图，并显示于显示设备；若存在滑坡，则选择显示流域滑坡时间分布图；

(5)参数输出模块输出流域的滑坡时间、滑坡位置、滑坡体大小统计值和滑坡预测函数，以此建立不同降雨条件下流域边坡稳定性分区图。

本发明中，所述步骤(1)中，建立三维网格时，DEM数据无需进行洼地处理，保留地形的原始特征。

本发明中步骤(2)中，使用里查兹方程描述地下水在饱和或非饱和土体及孔洞中的流动：

$\▿\·f^a\stackrel{v}{q}\±q^b\±q^e=f^v\frac{\∂\mathrm{\φ}{S}_w}{\∂t}---\left(4\right)$

$\stackrel{v}{q}=-k_{\mathrm{rw}}\frac{{\mathrm{\ρ}}_{w}g}{{\mathrm{\μ}}_w}\stackrel{v}{k\▿}(\mathrm{\ψ}+h)---\left(5\right)$

式中：

为达西流量；q^b为各种边界上的输入/输出项；q^e为地表与地下土体间的交换水量；φ为孔隙率；S_w为土壤饱和度；t为时间；f^a为与地表/地下土体对应的面积；f^v为与地表/地下土体对应的体积，k_rw为相对渗透率，ρ_w为水密度，g为重力加速度，u_w为流体动力粘度，为固有渗透率矢量，h为地表高程，ψ为压强水头。

本发明步骤(3)中，对于任意地点，安全系数的计算从表层逐层计算到底层，选取安全系数最小的层位作为潜在的滑动位置。

本发明步骤(5)中，还包括选择不同降雨特性下流域首次进入不稳定状态的时间以及完全恢复稳定状态的时间作为判定标准，通过指数函数拟合降雨过程中和降雨结束后流域安全系数临界曲线，将预测区域分为无条件稳定，稳定以及不稳定区域。

本发明的有益效果是：

本发明将滑坡确定性模型与基于物理概念的水文模型进行了有效的集成，分析了在降雨作用下流域边坡稳定性的时空分布规律以及相互间的定性定量关系，实现对流域滑坡进行快速有效的分析和预测。采用水力坡度方向作为滑坡体滑动方向，克服了以往模型在复杂地形需要洼地预处理的不足，保留了自然流域地貌水文响应特性，滑坡体更接近于实际滑坡体形状，降低预测误差。同时，考虑到土体含水量的空间分布差异，采用积分表达式求解的滑动力和抗滑力，提高了模型的精度和适用性。本发明的滑坡预测方法适用于简单边坡、复杂边坡及流域尺度的滑坡预测。

附图说明

图1为系统流程图；

图2为边坡滑动示意图；

图3为降雨历时三天，不同雨强下安全系数变化趋势图；

图4为降雨历时一天，不同雨强下安全系数变化趋势图；

图5为降雨历时三天，流域稳定性区域随雨强变化趋势图；

图6为降雨历时一天，流域稳定性区域随雨强变化趋势图。

具体实施方式

本发明以基于物理概念的水文模型Integrated Hydrology Model(InHM)为基础，采用改进的无限边坡稳定性模型，考虑水力坡度方向作为滑坡体滑动方向，分析降雨作用下流域边坡稳定性的时空分布规律，对不同降雨特性下流域滑坡进行预警预报。

该方法通过水文模型和滑坡模型的计算机数值模拟，流程如图1所示，主要包括以下部分：

1)参数输入模块，用于选定预测滑坡的流域，建立三维地形，输入土壤、植被、气象等数据。

2)水文计算模块，用于运行InHM水文模型，计算流域内的水文要素时空变化，包括地表水二维地下水三维流动，任意位置土壤含水量变化，任意位置水力坡度变化等。

3)滑坡计算模块，用于根据InHM水文模型所计算的水文要素，计算流域内任意三维位置的安全系数，若安全系数小于1，则为潜在的滑坡位置，记录该位置深度和时间，并计算体积。

4)可视化模块，用于显示指定时间流域安全系数分布图，显示指定时间流域滑坡深度图，若存在滑坡，显示流域滑坡时间分布图。

5)参数输出模块，用于输出滑坡时间、滑坡位置、滑坡体大小统计值，输出滑坡预测函数，建立不同降雨条件下流域边坡稳定性分区图。

具体地说，本发明是采用以下步骤对流域边坡稳定性进行预测：

1)参数输入模块根据DEM地形数据及土层厚度建立三维网格，DEM数据无需进行洼地处理，保留地形的原始特征，输入不同位置的土壤、植被、气象数据，生成模型所需的边界条件。

2)水文计算模块采用基于物理概念的水文模型，使用里查兹方程描述地下水在饱和或非饱和土体及孔洞中的流动，在降雨的驱动下计算土体含水量的变化过程，计算水力坡度的变化过程。

在InHM中，三维地下水在饱和或非饱和土体及孔洞中的流动表示为：

$\▿\·f^a\stackrel{v}{q}\±q^b\±q^e=f^v\frac{\∂\mathrm{\φ}{S}_w}{\∂t}---\left(\text{6}\right)$

$\stackrel{v}{q}=-k_{\mathrm{rw}}\frac{{\mathrm{\ρ}}_{w}g}{{\mathrm{\μ}}_w}\stackrel{v}{k\▿}(\mathrm{\ψ}+h)---\left(7\right)$

式中：

为达西流量；q^b为各种边界上的输入/输出项；q^e为地表与地下土体间的交换水量；φ为孔隙率；S_w为土壤饱和度；t为时间；f^a为与地表/地下土体对应的面积；f^v为与地表/地下土体对应的体积，k_rw为相对渗透率，ρ_w为水密度，g为重力加速度，u_w为流体动力粘度，

为固有渗透率矢量，h为地表高程，ψ为压强水头，模型的求解详见文献(VanderKwaak，1999)。

3)滑坡计算模块，改进无限边坡稳定模型，采用水力坡度方向作为滑动面方向，建立极限平衡方程，充分考虑地下水含水量对土体滑动力与抗滑力的影响，如图2所示。采用水力坡度方向作为潜在的滑动方向，可以减少地表起伏对滑坡方向的影响。安全系数(safety factor，SF)表达式如下：

$\mathrm{SF}=\frac{c^{\′}+c_r+{\cos }^2\mathrm{\θ}\tan \mathrm{\φ}{\∫}_0^z{\mathrm{\ρ}}^{\′}\mathrm{gdz}}{\sin \mathrm{\θ}\cos \mathrm{\θ}{\∫}_0^z\mathrm{\ρgdz}}---\left(8\right)$

式中SF表示任意点的安全系数，c′、c_r分别为土的有效粘聚力和根系粘聚力，kPa；θ为滑动方向(水力坡度方向)，°；φ为有效内摩擦角，°；ρ、ρ′分别为土壤密度和土壤有效密度，kg/m³；z为土体深度，。

土壤密度可表示为含水量的函数，公式如下：

ρ＝n_sS_wρ_w+(1-n_s)d_sρ_w                       (9)

其中d_s为土壤颗粒比重；n_s为土壤孔隙度；S_w为土体含水量；ρ_w为水密度。式(8)与(9)中，c′、c_r、φ、d_s和n_s可由土壤特性数据给定，S_w通过InHM求解得到。

水力坡度采用线性平均方式求解，公式如下：

$\mathrm{\θ}=\arctan (\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{(\mathrm{\ψ}+h)}_i}{\mathrm{\Δ}{l}_i}/n)---\left(10\right)$

式中Δ(ψ+h)_i表示该点与其相连的第i节点总水头差，Δl_i表示示该点与其相连的第i节点距离，n表示与该点相连的节点数，可通过网格结构和InHM模型求解得到。

将式(9)和(10)代入式(8)可求得边坡上任意点沿深度方向的安全系数。对于任意地点，安全系数的计算从表层逐层计算到底层，选取安全系数最小的层位作为潜在的滑动位置，记录时间与深度，计算该点所代表的体积，如图2所示。当SF＞1，该位置是稳定的，反之，则不稳定。

4)可视化模块，显示指定时间流域安全系数等值线图，若存在滑坡，显示流域滑坡时间等值线及滑坡体深度等值线。

5)参数输出模块，输出滑坡时间、滑坡位置、滑坡体大小统计值。选择不同降雨特性下流域首次进入不稳定状态的时间以及完全恢复稳定状态的时间作为判定标准，通过指数函数拟合降雨过程中和降雨结束后流域安全系数临界曲线，将预测区域分为无条件稳定，稳定以及不稳定区域。

另外，需要强调的是，本领域技术人员在全面理解本发明所述技术内容之后，根据其所掌握的土力学、水力学、计算机软硬件、气象预报学等方面的知识，完全能够以各种方式实现本发明所述的参数输入模块、水文计算模块、滑坡计算模块、可视化模块和参数输出模块，并据此再现本发明的内容。因此，本发明对此不再赘述。

具体的实施例子：

选择浙江省一个典型小流域作为对象，研究台风暴雨作用下流域边坡稳定性的时空分布规律和流域稳定性的预报。通过将模型应于该小流域的边坡稳定分析，可以确定在不同降雨特性下流域边坡的稳定性变化情况，并实现不同降雨特性下的流域稳定性预测，为灾害预警预报提供快速有效的方法。

1)流域概况，建立三维地形，确定土壤、植被、气象等数据

1.1流域概况

流域位于浙江省台州市宁溪镇长潭水库上游，流域出口位置坐标为28.573°N，120.971°E，流域面积174367m²，高程范围64.4m至365.3m。流域出口正对居民区，在强降雨作用下潜在的滑坡可能性为当地居民的安全埋下了隐患。

1.2流域地形

采用分辩率为5m的栅格DEM数据作为原始数据，建立研究区域不规则三角网(Triangulated irregular network，TIN)模型。TIN模型水平方向平均分辨率为5m，垂向分辩率为0.1m，共21层，其中1-20层为土层，第21层为基岩层。由于基岩相对于土体为不透水层，基岩的层数对模拟结果影响可忽略，为了提高计算效率，故仅设为一层。TIN模型网格总节点数为508313，三角形单元数为993324。

1.3流域土体及植被

流域内森林覆盖率达90％以上，主要为针叶林，阔叶林以及蕨类植物。研究流域平均土层厚度约2m，底下基岩为花岗岩，岩石渗透率取花岗岩的平均渗透率6.36×10^-7m/s。根据现场采集的土样，并分析其级配曲线可知，土体为粘质粉土，孔隙度50.1％，土壤颗粒比重为2.705。根据粘质粉土相关研究，粘质粉土的饱和渗透率范围是10^-3～10^-4m/s，内摩擦角的范围为：35°～45°。由于粘质粉土的粘聚力比较小，且植被覆盖状况与根系粘聚力之间缺少确切的转化方式，目前广泛采用的方法是将土壤粘聚力、根系粘聚力均设为0kpa，同时，将内摩擦角取值为45°，弥补粉土粘聚力与根系粘聚力对抗滑力的作用。所示。土壤含水量随着气象条件变化，由于研究区域尚未有长序列含水量数据，模拟中采用的初始含水量假设为30％。

1.4流域气象

研究流域属于典型的亚热带季风气候，年降水量1519mm，年蒸发量为1360mm，台风平均每年发生二、三次，最多可达四、五次，影响季节一般为6～9月。本研究以近十年来台风破坏最严重的2004年为参考，降雨数据主要根据2004年14号台风“云娜”以及18号台风“艾利”登陆后该地区的降雨情况来确定。台风“云娜”于8月12日在浙江台州登陆，是1956年以来在我国登陆的最强一次台风，此次台风所带来的降雨强度大，持续时间长，流域三天平均降雨量达581.3mm，其中8月12日一天流域内平均降雨量达到495mm，超过长潭水库50年一遇的一日设计暴雨量(480mm)，局部地区12日一天降雨量就达645mm。第18号台风“艾利”紧跟14号台风之后，在福建福清登陆，虽风力比较小，但水库流域内的过程降雨量也达到201.2mm。根据上述资料，确定本次模拟的降雨强度与降雨历时的组合为以下两类：(1)降雨强度20mm/d至200mm/d，降雨三天；(2)降雨强度200mm/d至550mm/d，降雨一天。

2)设计模拟工况，运行模型，计算流域安全系数时间分布

台风暴雨降雨特性主要包括降雨强度，降雨历时，降雨移动方向，降雨不均匀性等，对于小流域，后两者对水文响应的影响基本可以忽略，本研究考虑降雨强度和降雨历时两种特性组合的影响。根据台风暴雨资料所确定的降雨强度和降雨历时范围，同时考虑到流域滑坡预测函数的需要，确定模拟的降雨强度与降雨历时的组合为以下两类，共18个工况：(1)降雨强度20mm/d至200mm/d，降雨3天；(2)降雨强度200mm/d至550mm/d，降雨1天，具体工况如表1所示。

表1 模拟工况及模拟参数

序号	降雨强度，mm/d	降雨历时，d	其它参数
					1	20	3	网格层数，-	21
2	40	3	土层厚度，m	2
					3	60	3	岩石层厚度，m	0.1
4	80	3	孔隙度，-	0.501
					5	100	3	土壤饱和渗透系数，m/s	5.0×10-4
6	120	3	岩石饱和渗透系数，m/s	6.36×10-7
					7	140	3	初始含水量，％	30

8	160	3	土壤粘聚力，kPa	0
					9	180	3	根系粘聚力，kPa	0
10	200	3	内摩擦角，°	40
					11	200	1
12	250	1
					13	300	1
14	350	1
					15	400	1
16	450	1
					17	500	1
18	550	1

3)可视化输出任意时间的流域安全系数、滑坡体分布等，输出任意位置安全系数随降雨的变化过程

3.1流域边坡安全系数发展过程

为方便比较，选取降雨历时3天的临界雨强(工况4)，降雨历时3天的极端雨强(工况10)以及降雨历时1天的极端雨强(工况18)的模拟结果，对比流域在不同降雨特性下边坡安全系数在第1天，第3天以及第6天的分布情况，边坡安全系数及其分布情况呈现较明显的规律。1)流域边坡不稳定区域由流域上游逐渐扩大到流域下游，由山脊逐渐向山谷和渠道区域扩大；2)流域边坡不稳定出现时间与雨强、历时相关，当雨强较小时，不稳定时间较晚，出现在降雨结束之后一段时间，当雨强较大时，不稳定时间较早，出现在降雨结束时或之前；3)在降雨强度小，降雨历时长的情况下，流域边坡出现不稳定和恢复稳定时间较晚，反之，流域边坡出现不稳定和恢复稳定时间较早；4)滑坡体的深度和滑坡体体积随着降雨过程逐渐增大，潜在滑动面从地表附近逐渐过度到基岩附近。对于整个流域，陡坡处滑坡体滑坡时间较早，深度大，产生的滑坡体体积也较大，缓坡处滑坡体滑坡时间较晚，深度小，产生的滑坡体体积也较小。

3.2流域边坡最不安全系数发展过程

流域内最不稳定位置随降雨特性的变化而变化，由最不稳定位置触发的滑坡可同时引起周边不稳定区域发生滑坡，滑坡体向流域出口位置移动，因此，探索流域内最不稳定位置随降雨特性的变化规律，掌握流域可能发生滑坡的时间以及流域完全恢复稳定的时间，对于防灾具有重要的意义。在降雨历时3天的情况下：1)当雨强小于80mm/d，3天的降雨不足以引起边坡不稳定，当雨强大于80mm/d时，流域将在不同时期处于不稳定状态；2)随着雨强的增大，流域边坡最不稳定的时间提前，当雨强为150mm/d时，流域边坡最不稳定时间基本出现在降雨结束时，当雨强大于150mm/d时，流域边坡最不稳定时间在降雨结束前出现，反之，在降雨结束后出现；3)当雨强大于60mm/d时，流域边坡最不安全系数随时间变化存在一个拐点，使得最不安全系数迅速减小，雨强越大，拐点出现的时间越早；4)雨强较小时，流域边坡出现不稳定和恢复稳定时间较晚，反之，流域边坡出现不稳定和恢复稳定时间较早，如图3所示。当降雨历时为1天时，流域最不安全系数随时间的变化趋势与降雨历时3天的变化趋势基本一致，雨强大于200mm/d时，1天的降雨就能引起边坡不稳定；雨强大于400mm/d时，流域边坡最不稳定时间基本集中在降雨结束附近，如图4所示。

4)推算流域滑坡稳定性预测函数，为流域稳定性变化快速预测提供辅助

通过流域最不安全系数可以计算流域进入不稳定状态以及恢复稳定状态的时间，对流域边坡稳定性进行预测。选择不同降雨特性下流域首次进入不稳定状态的时间以及完全恢复稳定状态的时间作为判定标准，通过指数函数分别拟合降雨过程中和降雨结束后流域安全系数临界曲线，预测区域分为无条件稳定，稳定以及不稳定，预测曲线L1～L4函数表达式如表2所示。无条件稳定区域为对应降雨历时下，雨强小于保持边坡稳定的临界雨强，该流域在降雨历时3天(图5)的和降雨历时1天(图6)情况下，无条件稳定对应的临界雨强分别为70mm/d和170mm/d。当雨强较小时，流域进入不稳定状态时间较晚，但是一旦进入不稳定状态，完全恢复稳定的时间也较晚，随着雨强的增大，流域进入不稳定状态的时间以及降雨结束后完全恢复稳定的时间较快，由此可见，在实际预报中，不能简单地以有限时间内边坡稳定性情况来判断流域是否稳定，需要以边坡稳定性变化趋势作为辅助参考依据。

表2 流域稳定性预测函数

注：A，T，B分别为待拟合的参数，e为自然底数，变量x为时间。

Claims (5)

1.降雨作用下流域滑坡时空预测方法，包括以下步骤：

(1)根据流域的DEM地形数据及土层厚度建立三维网格，向参数输入模块输入不同位置的土壤、植被、气象数据，生成InHM水文模型所需的边界条件；

(2)根据其边界条件计算流域内的水文要素时空变化，使用InHM水文模型模拟地表水二维流动和地下水三维流动、任意位置土壤含水量变化和任意位置水力坡度变化的水文要素；

(3)无限边坡稳定模型以水力坡度方向作为滑动面方向建立极限平衡方程，以水力坡度方向作为潜在的滑动方向，并将之用于滑坡计算模块；根据InHM水文模型生成的水文要素，使用滑坡计算模块计算流域内任意三维位置的安全系数SF；若安全系数SF小于1，则为潜在的滑坡位置；记录该位置的深度、时间和体积；

所述任意点的安全系数SF表示为：

$\mathrm{SF}=\frac{c^{\′}+c_r+{\cos }^2\mathrm{\θ}\tan \mathrm{\φ}{\∫}_0^z{\mathrm{\ρ}}^{\′}\mathrm{gdz}}{\sin \mathrm{\θ}\cos \mathrm{\θ}{\∫}_0^z\mathrm{\ρgdz}}---\left(1\right)$

式中：c′、c_r分别为土的有效粘聚力和根系粘聚力，单位为kPa；θ为滑动方向，即水力坡度方向；φ为有效内摩擦角，°；ρ、ρ′分别为土壤密度和土壤有效密度，单位为kg/m³；z为土体深度；

而土壤密度表示为含水量的函数，具体如下：

ρ＝n_sS_wρ_w+(1-n_s)d_sρ_w      (2)

其中d_s为土壤颗粒比重；n_s为土壤孔隙度；S_w为土体含水量；ρ_w为水密度；

式(1)与(2)中，c′、c_r、φ、d_s和n_s由土壤特性数据给定，S_w通过InHM水文模型获得；

水力坡度方向θ通过线性平均方式获得，计算公式如下：

$\mathrm{\θ}=\arctan (\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{\mathrm{\Δ}{(\mathrm{\ψ}+h)}_i}{{\mathrm{\Δl}}_i}/n)---\left(3\right)$

式中，△(ψ+h)_i表示该点与其相连的第i节点总水头差，Δl_i表示该点与其相连的第i节点距离，n表示与该点相连的节点数；

将式(2)和(3)代入式(1)可获得边坡上任意点沿深度方向的安全系数；当SF＞1，该位置是稳定的，反之则不稳定，表示存在滑坡；

(4)通过可视化模块，将步骤(3)所得的数据拟合为指定时间流域安全系数分布图和指定时间流域滑坡深度图，并显示于显示设备；若存在滑坡，则选择显示流域滑坡时间分布图；

(5)参数输出模块输出流域的滑坡时间、滑坡位置、滑坡体大小统计值和滑坡预测函数，以此建立不同降雨条件下流域边坡稳定性分区图。

2.根据权利要求1所述的降雨作用下流域滑坡时空预测方法，其特征在于，所述步骤(1)中，建立三维网格时，DEM地形数据无需进行洼地处理，保留地形的原始特征。

3.根据权利要求1所述的降雨作用下流域滑坡时空预测方法，其特征在于，步骤(2)中，使用里查兹方程描述地下水在饱和或非饱和土体及孔洞中的流动：

$\▿\·f^a\stackrel{\→}{q}\±q^b\±q^e=f^v\frac{\∂\mathrm{\φ}{S}_w}{\∂t}---\left(4\right)$

$\stackrel{\→}{q}=-k_{\mathrm{rw}}\frac{{\mathrm{\ρ}}_{w}g}{{\mathrm{\μ}}_w}\stackrel{\→}{k}\▿(\mathrm{\ψ}+h)---\left(5\right)$

式中：

为达西流量；q^b为各种边界上的输入/输出项；q^e为地表与地下土体间的交换水量；φ为孔隙率；S_w为土壤饱和度；t为时间；f^a为与地表/地下土体对应的面积；f^v为与地表/地下土体对应的体积，k_rw为相对渗透率，ρ_w为水密度，g为重力加速度，μ_w为流体动力粘度，

为固有渗透率矢量，h为地表高程，ψ为压强水头。

4.根据权利要求1所述的降雨作用下流域滑坡时空预测方法，其特征在于，步骤(3)中，对于任意地点，安全系数的计算从表层逐层计算到底层，选取安全系数最小的层位作为潜在的滑坡位置。

5.根据权利要求1所述的降雨作用下流域滑坡时空预测方法，其特征在于，步骤(5)中，还包括选择不同降雨特性下流域首次进入不稳定状态的时间以及完全恢复稳定状态的时间作为判定标准，通过指数函数拟合降雨过程中和降雨结束后流域安全系数临界曲线，将预测区域分为无条件稳定，稳定以及不稳定区域。

Images