基于覆盖土层厚度识别的区域斜坡降雨入渗稳定性定量评价

Regional Hillslope Stability Analysis under Rainfall Based on Characterization of Overburden Soil Layer Thickness

摘要

区域内覆盖土层孕育着斜坡浅表层破坏，是区域斜坡破坏的最主要来源，但由于覆盖土层厚度的空间变化，导致堆积层滑坡灾害发育的底界难以合理确定，进而使稳定性评价结果未能指示实际滑坡状况.以黄冈市九资河镇为例，对该区域的覆盖层厚度进行详细调查，并通过基于地貌过程的估计方法确定了研究区覆盖层厚度的空间分布；在此基础上，运用斜坡稳定性定量模型，分析了研究区斜坡在降雨条件下的稳定性，并据此估计了区域堆积层滑坡的分布情况和可能性.结果显示：所提出的方法有效且具有很强的实用性，基于GIST模型估计的区域覆盖层厚度与实际情况较为接近（集中在0.5~3.0 m）.覆盖层厚度的合理确定使得区域斜坡降雨入渗稳定性评价结果更加精细且具有实际滑坡指示意义；该区基本稳定和稳定斜坡主要分布于阶地平台或河漫滩等覆盖层厚度大但坡度较小区域，欠稳定及不稳定斜坡主要分布在水库及各支流靠岸位置，这些部位覆盖层厚度虽小但坡度较大、地形外凸且受地下水影响明显；短时间的强降雨对斜坡稳定性影响巨大.该研究有助于推动区域斜坡灾害评价向精细化方向发展.

关键词

地貌过程 / 覆盖层厚度 / 区域斜坡稳定性 / 降雨入渗 / 滑坡 / 工程地质.

Key words

geomorphic processes / overburden thickness / regional slope stability / rainfall infiltration / landslides / engineering geology

0 引言

目前，区域尺度的斜坡稳定性评价多是以易发性和危险性着手，方法上也主要是基于经验驱动的定性方法或基于数据驱动的统计分析（半定量）方法.如高克昌等（2006）以重庆万州滑坡为例，运用信息量模型对区域滑坡进行危险性评价；庄建琦等（2010）采用信息量法，结合GIS技术对汶川地震极重灾区域崩塌滑坡进行危险性评价；宋盛渊等（2014）以岩桑树水电站库区潜在滑坡为例，基于突变理论进行滑坡危险性评价；胡涛等（2020）基于逻辑回归模型和3S技术对贵州省思南县进行滑坡易发性评价；胡燕等（2020）以巴东县城为例，基于证据权法开展滑坡易发性评价研究；黄发明等（2022）基于连续概率降雨阈值，结合C5.0决策树和多层感知器模型对江西宁都县进行危险性分析；杨灿等（2022）利用贝叶斯算法对4种常见滑坡易发性评价机器学习模型进行优化并探索该算法的优化效果；黄武彪等（2022）基于层数自适应加权卷积神经网络对川藏交通廊道沿线滑坡进行易发性评价；Millerand Degg（2009）基于多元逻辑回归模型并结合GIS技术对威尔士东北部区域滑坡进行危险性区划；Cascini et al.（2010）基于地下水模型和滑坡位移预测模型分析了降雨与地下水位变化及滑坡位移之间的相关性，从而对意大利中部的降雨型滑坡进行危险性分析；然而，这类方法得到的稳定性结果只能在一般规律或统计意义上显示区域斜坡平均稳定状况.

近些年，越来越多的学者采取基于物理本构模型的力学计算方法开展区域尺度的斜坡稳定性评价定量研究.可用于稳定性分析的物理本构模型包括：Montgomery and Dietrich（1994）改进的SHALSTAB （Shallow Landsliding Stability）模型以及Pack et al.（1998）提出的 SINMAP （Stability INdex Mapping）模型，这两个模型都是以地下水径流和地表水下渗处于稳流状态为前提的浅层水文模型；Baum et al.（2002）基于降雨入渗诱发滑坡理论对区域浅层斜坡在瞬态非饱和入渗条件下的稳定性评价方法展开了研究，提出了优化的TRIGRS 2.0（Transient Rainfall Infiltration and Grid-based Regional Slope-Stability Analysis）模型.

当采取基于物理本构模型的力学计算方法评价区域斜坡稳定性时，模型参数的确定成为影响评价效果的关键环节.在斜坡稳定性的计算中，覆盖土层厚度、地下水埋深以及物理力学参数等均是对计算结果有直接影响的重要参数.特别地，区域覆盖土层孕育着斜坡浅表层破坏，是区域斜坡破坏的最主要来源；覆盖土层厚度的空间变化情况，圈定了堆积层滑坡灾害发育的边界范围，是分析区域斜坡稳定性的关键信息.然而，目前区域尺度的斜坡系统中对覆盖层厚度认识较少且难以全面测量，仅在小尺度范围且已开展详细勘探的场地才有详细的覆盖层厚度数据，因此区域斜坡稳定性评价结果还难以指示实际滑坡可能.

如刘磊等（2016）以三峡库区万州主城区为例，以瞬态降雨入渗的区域斜坡稳定性计算模型为基础，将一定持续降雨期内各网格单元失稳的概率视为滑坡危险性，并根据岩土物理力学参数的不确定性求解此概率，从而评价区域内滑坡的危险性.在其研究中，研究区覆盖土层厚度分布是基于9个滑坡勘察厚度数据，通过厚度与阶地、水文以及地貌环境等地表因素之间的关系，采用工程地质类比的方法应用统计回归模型推估的，即是基于有限数据的一种区域拟合.因此，基于此精度下的覆盖土层厚度的危险性计算结果仅能在一定程度上反映研究区斜坡稳定性的时空分布特征，其效果有待提高.

区域覆盖土层厚度的确定之所以会成为制约区域斜坡稳定性评价效果的关键因素，究其原因，主要是对一整个区域来说，尤其是多山地区，范围大且条件复杂，详细调查覆盖土层厚度既不经济亦较难实现.若有更加精确的区域覆盖土层厚度分布信息，相应的区域斜坡稳定性评价结果将更为真实可靠.因此，在应用确定性模型时，利用现有的数据以及技术手段得到较为真实的模型参数显得尤为重要.

为此，一些学者做了一定的研究，如：曾立峰等（2012）根据钻孔资料、场地条件等因素，针对研究区覆盖层厚度采用了1/4波长法和卓越频率统计法进行反演，从而得到覆盖层厚度的空间分布；Catani et al.（2010）提出将覆盖层厚度与斜坡坡度、斜坡曲率以及斜坡相对位置3个地貌指标联系起来，并基于修正的经验公式估算区域覆盖层厚度，建立了地貌学覆盖层厚度模型（GIST），这为区域覆盖土层分布的确定提供了新的思路.

区域尺度上的覆盖土层厚度分布受到多个因素的影响，例如：基岩风化程度、气候影响、斜坡坡度、斜坡曲率、上坡贡献区域（斜坡剖面内的相对位置）、植被覆盖及外界营力等.Catani et al.（2010）研究发现：覆盖层厚度与曲率之间呈负相关；其次，斜坡坡度也是影响土体累积和流失的重要因素，当斜坡坡度接近某一临界值时，滑坡发生的规模和频率会逐渐增大，进而影响覆盖层厚度分布；此外，影响覆盖层厚度的另一个因素是斜坡剖面内的相对位置.例如，由于外界因素使得坡顶及坡腰的岩土体发生滑移，坡脚处地形凸起实际上是高覆盖层厚度的标志，但位于坡顶位置处凸起则表征岩土体发生侵蚀形态，具有较浅的土层深度.但这两种情况下，斜坡曲率和坡度可能是一致的.为了避免这种情况被忽略，考虑斜坡剖面中的相对位置是非常有必要且符合实际情况，这些因素的复杂性导致了覆盖层厚度的估计极具挑战.

覆盖层厚度是基于降雨‒力学模型定量评价区域斜坡稳定性的重要参数.鉴于此，以湖北省黄冈市罗田县九资河镇为研究区，对该区域的覆盖层厚度进行详细调查，并通过基于地貌过程的方法估计确定了研究区空间分布的覆盖层厚度，在此基础上，运用斜坡稳定性定量模型，分析了研究区斜坡在降雨条件下的稳定性，并据此估计了区域浅层滑坡的分布情况和可能性.

1 研究区地质环境条件

研究区所在的罗田县地处大别山南麓，属中山、低山及丘陵区地貌.受北东东向构造影响，山脉多呈北东展布；地势北东高，南西低.区内地层主要有太古界片麻岩地层、燕山期花岗岩和新生界第四系松散堆积物地层（图1）.

该区处于北半球亚热带季风气候区，气候温和、雨量充沛，年平均降雨量1 334.5 mm，降雨主要集中在4月至8月，日平均降雨量 176.1 mm.降雨具突发性强、强度大、连续集中的特征.区内地下水动态变化受季节影响明显，表现出随季节涨落幅度大、速度快的特征.

调查结果显示，罗田县地质灾害以小型滑坡为主，多发育于构造侵蚀低山丘陵区、构造侵蚀中低山区.由于这些部位堆积条件较好，堆积物往往成为滑坡的物质来源.堆积物地层主要由第四系全新统冲积、冲洪积及残坡积的粉质粘土、砂、砂砾石等组成，结构松散，力学性质差.据统计，区内约90%的滑坡为以残坡积物为主的堆积层滑坡，其受大气降雨周期性控制，表现为雨季灾害多发、频发的特点.

对于以浅表堆积层滑坡为主要地质灾害类型的地区来说，区域覆盖土层厚度分布情况对于该区域斜坡稳定性无疑是关键因素.覆盖层厚度，即覆盖于基岩之上的第四系松散堆积层的厚度，是控制斜坡稳定性系数的关键参数之一（贾琳等，2021）.堆积层滑坡可以在坡表与基覆界面之间的任一深度处发生，最大的滑动面深度就是覆盖层厚度.可见，研究区覆盖土层厚度分布的识别、确定对于该区的斜坡稳定性评价至关重要.

2 区域覆盖层厚度估计方法

2.1　基于地貌学的覆盖层厚度估计方法

覆盖层厚度的估计主要是基于地貌学覆盖层厚度模型（Geomorphologically Indexed Soil Thickness Model，简称GIST），该模型由Catani et al.（2010）提出.该方法以斜坡为基本单元，分别对每一个斜坡的斜坡曲率、斜坡坡度和斜坡相对位置厚度进行转化，得到数值为0~1之间的指标：斜坡曲率贡献指数（C）、坡度临界值贡献指数（ψ）及斜坡相对位置贡献指数（η），详述如下：

2.1.1　斜坡曲率贡献指数C

C的确定是将真实曲率值进行线性归一化处理.真实曲率值有正负，即凸为正、凹为负.在线性归一化时，最大的真实曲率值（正值）对应到1，最小的曲率值（负值）对应到0，曲率贡献指数越大，则覆盖层厚度越小.

本文使用的真实曲率值是基于DEM计算的标准曲率，其综合了剖面曲率与平面曲率.其中剖面曲率影响地表径流的加速和减速，平面曲率影响地表径流的汇聚和分散，因此标准曲率更能表征覆盖层的侵蚀和沉积情况.

2.1.2　坡度临界值贡献指数 \psi

\psi通过下式确定：

\psi = \left\{\begin{array}{cc}{\left(\mathrm{1}+\mathrm{t}\mathrm{a}\mathrm{n}{\theta }_{\mathrm{t}\mathrm{h}}\right)}^{-\mathrm{1}}& \theta >{\theta }_{\mathrm{t}\mathrm{h}}\\ \mathrm{1}& \theta \le {\theta }_{\mathrm{t}\mathrm{h}}\end{array}\right.，（1）

其中，θ _th是坡度临界值.θ _th可以基于不同地层岩性的内摩擦角并考虑粘聚力和植被固结作用综合确定，亦或是通过反分析及观察、统计等方法进行确定.一般情况下坡度临界值稍大于内摩擦角.

2.1.3　斜坡相对位置贡献指数 \eta

斜坡相对位置通俗来讲就是斜坡单元内某点到坡顶位置的水平距离（p），将坡顶位置设为“0”，坡脚位置设为“1”，其值在0~1之间.

基于均一化的斜坡相对位置p，确定与p对应的位置贡献指数η.这需要通过野外实测的覆盖层厚度拟合出p与η之间的函数关系.由于不同坡形斜坡单元形成的地貌过程不同，其覆盖层厚度在斜坡剖面上的分布规律不同，因而不同坡形斜坡p与η之间的函数关系需要分别确定.拟合出的η（0~1之间的纯数字）越大代表该点处的覆盖层厚度值越大.

2.1.4　覆盖层厚度的估计

将以斜坡为基本单元分别确定的C、ψ和η代入以下经验公式，即可确定覆盖层厚度：

h = K_{\mathrm{c}} \cdot \left(\mathrm{1}-C\right) \cdot \eta \cdot \psi，（2）

其中，h为覆盖层厚度估计值；K _c为校准系数.K _c的目的是将（1-C）、ψ和η之积的纯数字无量纲值转换为实际的覆盖层厚度.K _c根据现场实测的覆盖层厚度拟合确定，通常按不同地层岩性分别确定一个值.在不同地层岩性的K _c差别不大时，也可以简单地按不同斜坡坡形类型分别确定，即每类斜坡坡形类型确定一个K _c值.

2.2　露头覆盖层厚度实测调查

2.2.1　野外覆盖层的界定

通过实地调查，研究区的堆积物分布主要分为3种，分别是第四系堆积区、残坡积区和花岗岩风化残积区.为准确测量覆盖层的厚度需要先辨别基覆界面的位置，而不同岩性区的基覆界面判别标准不同，详述如下：

第四系堆积区下部基岩出露较少，大多位于地面以下，需要一定的开挖才能找到，然后利用钢尺测量获得覆盖层厚度（图2c）.

残坡积区可明显见到基岩出露，基覆界面的判别方法有两种：一是根据覆盖层与基岩之间的颜色差异来判别（图2b）；二是通过地质锤敲击岩土进行判别.

花岗岩风化残积区的基岩出露明显，并且基岩颜色与覆盖层颜色有明显突变，基覆界面容易判别（图2d）.

2.2.2　调查方法

首先，划分研究区的斜坡单元及编号，并确定每一个斜坡单元的坡形（图2a）.在地图上布置调查路线，原则上尽可能穿越不同斜坡单元.

其次，将调查路线上的覆盖层断面露头处作为调查点，记录点位坐标、调查点位于斜坡的大致位置（坡顶、坡腰或坡脚）和覆盖层厚度等信息（图2b~2d）.

最后，将点位坐标及其覆盖层厚度值录入ArcGIS系统中，按斜坡单元坡形分组，并按调查点所处斜坡单元相对位置摆放，观察不同坡形分组内的调查点是否能够覆盖该坡形斜坡单元不同位置（图2a），对于覆盖不足的地方进行补充调查.

2.2.3　调查结果

本次调查共获得251个覆盖层厚度实测数据（图1和图2a），不同坡形斜坡单元内厚度调查点分布较分散，同一坡形不同斜坡单元内分布较均匀，可以为后续模型校准以及验证提供基础数据.

3 研究区覆盖层厚度分布的估计

基于覆盖层厚度实测数据，通过地貌学覆盖层厚度模型（GIST），对研究区覆盖层厚度分布进行估计.这里采用10 m×10 m的栅格单元，将研究区划分成1 099 358个栅格.

3.1　覆盖层厚度影响因素贡献指数确定

3.1.1　斜坡曲率贡献指数C

基于ArcGIS软件，获得研究区每个栅格的标准曲率，随后以斜坡为范围将标准曲率进行归一化处理，即每个斜坡单元内的标准曲率值被转换为0~1之间的纯数字（图3a）.

此时，得到的栅格值表示曲率对覆盖层土流失的贡献指数，栅格值越大则表明其对覆盖层厚度贡献值越小，反之则越大.

3.1.2　坡度临界值贡献指数 \psi

目前，对于坡度临界值的确定主要基于岩土体的力学性质，属于经验取值法.

本文开展相关野外筛分实验，获取了3种岩性区碎石土的颗粒级配特征，并基于力学参数的经验公式（李雷雷等，2019）估计了3种岩性的平均内摩擦角和平均粘聚力；同时，统计分析已发生滑坡的斜坡坡度与其力学参数之间的关系.最后综合确定研究区不同岩性区的坡度临界值（表1）.

随后，在ArcGIS软件中将每个栅格的坡度值与其对应的坡度临界值进行对比，并按公式（1）计算坡度临界值贡献指数ψ，得到坡度贡献指数分布图（图3b）.

由图3可知，对于坡度值大于坡度临界值的栅格单元，ψ小于1，实现了对坡度较大处由于覆土运移导致的覆盖层厚度折减.

3.1.3　斜坡相对位置贡献指数 \eta

为确定斜坡相对位置贡献指数，笔者首先要获得斜坡剖面的相对位置，具体步骤如下：

将ArcGIS中包含有DEM高程和流向的斜坡单元属性表导入MATLAB软件中，通过递归算法计算各斜坡单元内每个栅格点高程到坡顶位置（斜坡单元内高程最大点）的最短上坡水文距离和到坡脚位置（斜坡单元内高程最小点）的最短下坡水文距离.为了计算方便，将坡顶位置设为“0”，坡脚位置设为“1”.实质上是计算每个栅格单元最短上坡水文路径与总水文路径的比值，将实际位置进行归一化处理，得到的栅格值即为斜坡单元剖面中的相对位置（p）（图3c）.

其次，要拟合出不同坡形的斜坡相对位置（p）与贡献指数（η）之间的关系，步骤如下：（1）通过ArcGIS软件，将校准模型的157个实测覆盖层厚度数据按斜坡坡形进行分组；（2）将各坡形分组的覆盖层厚度数据分别进行归一化处理；（3）基于不同坡形斜坡单元实际位置（p）和相对位置贡献指数（η）两组数据，利用Origin软件对其进行拟合即可得到p与η之间的关系式（表2和图4）.

最后，基于已计算出的相对位置p，通过不同坡形斜坡单元的p与η的函数关系，计算出全区域的相对位置厚度贡献指数（η）（图3d）.

3.1.4　校准系数 K_c

野外调查数据显示，凹凸形斜坡、凸凹凸形斜坡和凸形斜坡覆盖层厚度的最大值分别为5.0 m、5.5 m和6.0 m，其最小值都为0 m.

根据斜坡坡形进行分组，将各组覆盖层厚度调查值的最大值和最小值与 \left(\mathrm{1}-C\right) \cdot \eta \cdot \psi的最大值（取平均数）和最小值（取平均数）进行对比，即可确定 K_c，见表3.

3.2　覆盖层厚度估计结果

将确定的C、ψ、η以及 K_c代入公式2，即可获得研究区覆盖层厚度分布的估计结果（图5）.

表4给出了不同厚度区间的栅格数占总栅格数的百分比，可以发现0.5~1.5 m厚度区间的栅格数占比最大，约为37.15%；其次为1.5~2.0 m和2.0~3.0 m厚度区间，占比分别为28.58%和21.85%；大于3 m厚度区间占比最小，仅为0.19%.由此可见，研究区覆盖层厚度主要分布在0.5~3.0 m区间范围内.

3.3　覆盖层厚度估计的验证

3.3.1　模型对比

为了验证GIST模型的计算精度，本文利用反距离权重法（IDW模型）对研究区覆盖层厚度进行估算，并与GIST模型的覆盖层厚度计算结果进行了对比.

IDW模型是基于野外实测的覆盖层厚度数据通过空间插值的方法获得整个研究区的覆盖层厚度分布（图6a），并未考虑斜坡形态特征、曲率、陡坡修正和覆盖层形成演化过程.对比发现：基于GIST模型获取的覆盖层厚度（图6b），由于综合了斜坡形态、曲率、坡度等地貌特征，其结果更加符合实际；而IDW模型对于研究区范围小、实测点数量多的情况较适用.

3.3.2　误差分析

为了准确说明两种模型的使用效果，本文将收集的251个厚度调查点数据分为两组，其中一组包含157个厚度数据，用于校准模型参数；另一组包含94个厚度数据，用以验证模型的估计误差.通过将IDW模型和GIST模型获得的厚度估计结果与对应位置的94个覆盖层厚度实测值进行对比，验证覆盖层厚度的估计效果.

均方误差（MSE）是指参数估计值与参数真值之差平方的期望值，可以评价数据的变化程度，MSE值越小，代表预测模型描述实验数据具有更好的精确度，计算公式如下：

M S E = \frac{\mathrm{1}}{m} {{\sum }_{i=\mathrm{1}}^m\left(y_i-{\mathrm{ŷ}}_{\mathrm{i}}\right)}^{\mathrm{2}}，（3）

其中，y_i 为实测值，ŷ_i 为预测值.

图7为两种模型的覆盖层厚度估计值与实测值的散点图，明显看出，基于GIST模型得到的离散点更接近45°虚线，表明估计值与实测值更接近，且离散点均匀地围绕45°虚线分布，系统性误差小；同时将94个实测值和各点位的预测值代入公式（3），计算得到IDW模型的MSE值为2.61，基于GIST模型的MSE值为1.28.因此，基于GIST模型得到的覆盖层厚度估计值更加接近实际情况，预测效果更好.

4 降雨条件下区域斜坡稳定性评价

4.1　区域斜坡降雨入渗稳定性评价方法

TRIGRS （Transient Rainfall Infiltration and Grid-based Regional Slope-Stability Analysis）模型（Baum et al.， 2002）是美国地质调查局通过Fortran语言开发的一种区域栅格入渗‒稳定性计算模型，可以模拟真实降雨过程引起的区域尺度的瞬态孔隙水压力及斜坡稳定性随时间的变化.该模型采用栅格数据进行计算，不仅可以反映斜坡的地形与地质信息，而且计算的结果可以展示斜坡不同部位的稳定性情况.

TRIGRS模型由地表水文模型、入渗模型和斜坡稳定性模型3部分构成：

4.1.1　水文模型

TRIGRS 模型中的地表水文模型解决了降雨形成的地表水的分配问题，即当某栅格单元地表处的雨水来不及入渗时，模型会将多余的地表水根据流向分配到下游栅格.通常假定栅格单元的最大入渗量等于其饱和渗透系数.

4.1.2　入渗模型

TRIGRS中的入渗模型是基于一维无限斜坡模型（图8）建立，分析的是栅格单元（土柱）的渗流，水文模型为此提供了渗流的上边界条件.渗流场只考虑竖直方向的入渗，忽略了水平方向的侧向流动.具体的渗流计算公式详见Iverson（2000）提出的Richards方程线性解，可描述饱和‒非饱和渗流.

4.1.3　斜坡稳定性模型

斜坡稳定性模型同样是基于一维无限斜坡模型（图8）建立的，考虑了入渗模型确定的降雨入渗作用下栅格单元体不同深度的稳定性系数，忽略了栅格单元之间力的传递.稳定性系数是基于莫尔‒库仑破坏准则通过极限平衡分析确定.每个栅格不同时间步均需计算一次稳定性系数，对于任一栅格的任一时间步，稳定性系数 F s的计算公式如下：

F s \left(z,t\right) = \frac{\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{\text{'}}}{\mathrm{t}\mathrm{a}\mathrm{n}a} + \frac{c\mathrm{\text{'}}-\Psi (Z,t){\gamma }_{\mathrm{w}}\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{\text{'}}}{{\gamma }_{\mathrm{s}}\mathrm{s}\mathrm{i}\mathrm{n}a\mathrm{c}\mathrm{o}\mathrm{s}a}，（4）

式中，Fs为稳定性系数；t为入渗时间；a为斜坡坡度；φ'为覆盖土体的有效内摩擦角；c'为覆盖土体的有效粘聚力；ψ为压水头；γ _w为地下水重度；γ _s为土的重度；Z=z/cos a，z为垂直于坡面方向的某一潜在滑动面深度，其最大值为覆盖土层厚度；Z为坡表到潜在破坏面的垂直距离.

斜坡栅格单元的稳定性系数Fs（t）是该栅格单元沿着任一深度潜在滑动面的稳定性系数Fs（z， t）的最小值，表达式如下：

F s ( t ) = \mathrm{m} \mathrm{i} \mathrm{n} ( F s ( z , t ) )，（5）

当Fs（t）小于1表明t时刻斜坡在该栅格单元处不稳定.

4.2　基于覆盖层厚度识别的区域斜坡降雨入渗稳定性定量评价方法

将GIST估计的区域空间覆盖层厚度分布结果输入TRIGRS模型，即可实现基于覆盖层厚度识别的区域斜坡降雨入渗稳定性定量评价（图9）.

4.3　评价模型参数确定

4.3.1　流向

水流流向决定了降雨地表水的径流路径，是TRIGRS模型计算径流量和径流方向的重要参数之一，在ArcGIS软件中，通过DEM文件可以直接生成水流流向.需要注意的是，在计算流向之前需要对DEM进行填洼处理（李俊超等，2012），以最终得到研究区水流流向分布（图10a）.

4.3.2　初始地下水位

本文通过野外实地测量和查阅相关岩土勘察报告的方法获得研究区不同位置的地下水深度，然后利用ArcGIS软件对实测得到的地下水深度数据进行空间插值，得到研究区初始地下水位分布（图10b），并根据自然断点法将地下水深度划分为5个等级：1.2~2.1 m、2.1~2.5 m、2.5~2.8 m、2.8~3.3 m、3.3~4.0 m.

4.3.3　岩土体力学参数

TRIGRS模型中需要输入的力学参数包括：粘聚力、内摩擦角、重度等，前文将研究区划分为3种岩性区，针对不同岩性区收集已有勘察资料及相关文献，同时开展研究区碎石土的筛土实验，并基于经验公式（李雷雷等，2019）获取相关参数，详见表5.

4.3.4　水文参数

水文参数主要包括：水的重度、水力扩散系数、初始地表入渗率、垂直饱和渗透系数、体积含水率等.研究区水文参数较难获取，本文基于野外水文地质调查、水力学经验公式（朱崇辉等，2006；杨何等，2017）、参考相关勘察资料和文献确定了研究区水文参数（表6）.

4.4　计算工况确定

为模拟研究区在不同降雨强度下斜坡的稳定性状况，基于实际降雨数据（罗田县近61年的最大日降雨量为216.2 mm），在日累计降雨量不变的条件下，分别设置了2种降雨工况，各降雨强度及历时如表7所示.

工况一：为了模拟持续性均匀降雨，这里将降雨强度设定为18 mm/h，降雨时长为 12 h，每小时降雨强度相同.

工况二：为了模拟不同降雨强度和不同降雨历时，可分别设定降雨强度为10 mm/h、降雨历时为4 h；降雨强度为25 mm/h、降雨历时为4 h；降雨强度为50 mm/h、降雨历时为1 h；降雨强度为13 mm/h、降雨历时为2 h.四个阶段的降雨历时和降雨强度不尽相同.

经测算，在2.3 GHz 16 GB RAM的笔记本电脑上，通过TRIGRS模型进行一次1 099 359个栅格数据的12 h区域斜坡降雨入渗稳定性分析计算所需时间平均为15 min，可见计算效率较高，这对于模型的实际应用意义重大.

4.5　区域斜坡稳定性评价

根据计算得到的不同栅格稳定性系数大小，参考地质灾害防治工程勘察规范（DB50/143-2003），将区域斜坡稳定性划分为不稳定、欠稳定、基本稳定、稳定4个等级（表8）.

图11为工况一持续降雨12 h后的斜坡稳定性空间分布图，可以看出：持续均匀降雨12 h后，整个研究区基本稳定和稳定斜坡主要分布于阶地平台地区或河漫滩区域，该区域虽然覆盖层厚度普遍较厚，但坡度平缓，不易发育滑坡，其次是分布在地形内凹的地方；欠稳定及不稳定斜坡主要分布在天堂水库及各支流靠岸位置，这些区域虽然覆盖层厚度较小，但坡度较陡且地形外凸，地下水位浅，使覆盖土层容易达到饱和状态，导致斜坡稳定性较差，易形成滑坡.

图12为工况二不同降雨强度和降雨历时下的斜坡稳定性空间分布图，由图可知，随着时间的推移，累计降雨量的增加，不稳定斜坡的数量越来越多；尤其是2-3阶段，在1 h的强降雨作用下（50 mm/h），不稳定斜坡的数量迅速增加（图12b~12c）.

为准确揭示这种变化，计算得到了工况二不同阶段的斜坡稳定情况面积占比（图13），统计显示：2-2阶段不稳定斜坡面积占比为2.55%，而工况2-3阶段增加到15.06%，这主要是由于栅格单元体在短时强降雨条件下，降雨入渗使得岩土非饱和瞬时孔压增大，从而导致栅格单元稳定性降低，使得不稳定斜坡的面积增加.同时也表明短时的强降雨作用对斜坡稳定性影响很大.

同时，本文也计算出降雨结束时，两种工况的不同稳定性分区的面积占比（图14），计算结果显示：工况一（连续降雨12 h）的不稳定斜坡面积占整个研究区的11.86%，而工况二（连续降雨11 h）的不稳定斜坡面积增长到15.68%，说明在降雨总量相同的条件下，降雨历时越短，对斜坡稳定性越不利.

5 滑坡分布情况及效果评述

在覆盖土层厚度识别基础上，研究区斜坡稳定性评价结果较为精细地圈画出了研究区潜在滑坡位置.为了检验评价结果的准确性，本文选取接近实际降雨量的工况2-2阶段的斜坡稳定性评价结果与15处历史滑坡的分布情况和野外实际调查进行对比验证（图15和图16）.

图15显示了15处历史滑坡在研究区的位置，统计发现：2-2阶段的斜坡不稳定区仅占总面积的2.55%，但有13处历史滑坡点位于该区，占历史滑坡总数的86.7%；其余2处滑坡点分布在欠稳定区内.所有历史滑坡点均处于评价结果的不稳定区和欠稳定区内，表明基于GIST模型和TRIGRS模型计算得到的斜坡稳定性结果较为准确.

此外，本文还将野外调查的实际滑坡情况与稳定性评价结果进行对照检验（图16），对照结果都表现出较好的一致性，下面将通过闻家湾滑坡、胡家大湾滑坡和李河滑坡举例说明：

闻家湾滑坡位于花岗岩风化残积区内，覆盖层厚度为0.5~1.2 m，斜坡坡度较大，现场检验时正值降雨期，发现斜坡表层土体受降雨影响稳定性差，表层可见土体溜滑现象，且有大块巨石堆积在斜坡体表层，推测之前发生过大规模滑坡，这与稳定性结果图中不稳定斜坡区域位置一致.

胡家大湾滑坡位于残坡积区内，南侧为天堂水库，覆盖层厚度为1.5~2.0 m，斜坡坡度大，现场调查中坡体已经发生小范围滑移现象，滑坡后缘清晰可见，现场推测强降雨期可能再次发生滑坡，该区域所在位置与稳定性评价结果一致.

李河滑坡位于残积区内，西侧和南侧距离天堂水库较近，覆盖层厚度为1.0~1.5 m，斜坡坡度较大，受强降雨和库水位影响，滑坡中部发生较大变形，现场检验发现坡表局部仍在继续垮塌，坡体植被杂乱无章，该滑坡所在位置符合实际工况条件下稳定性区划结果.

综上表明：基于覆盖土层厚度识别的区域斜坡稳定性评价结果与野外调查结果较为符合，评价结果具有对实际滑坡的指示作用.

6 结论

论文针对区域斜坡稳定性定量评价做出尝试性工作：以湖北省黄冈市罗田县九资河镇为研究区，在调查认识基础上，建立了基于覆盖土层厚度识别的区域斜坡降雨入渗稳定性定量评价方法；利用该方法，通过地貌过程估计确定了研究区覆盖层厚度空间分布，分析了研究区斜坡在降雨条件下的稳定性，并据此圈定了区域浅层滑坡的位置并进行了验证.主要得到以下结论：

（1）基于地貌过程估计的区域空间覆盖层厚度与实际情况较为接近.研究区覆盖层厚度主要分布在0.5~3.0 m区间范围内.

（2）整个研究区基本稳定斜坡和稳定斜坡主要分布于阶地平台地区或河漫滩区域等覆盖层厚度虽大但坡度较小的地方，其次是分布于地形内凹的地方；欠稳定及不稳定斜坡主要分布在天堂水库及各支流靠岸位置等覆盖层厚度虽小但坡度较大、地形外凸且受地下水影响的部位，这些部位容易发育滑坡.

（3）降雨强度越高，累积降雨量越大，区域斜坡稳定性越差；降雨总量相同条件下，降雨历时越短，对斜坡稳定性越不利.短时间强降雨会对斜坡稳定性造成巨大影响.

（4）评价结果与野外调查对照表明，所提出的基于覆盖土层厚度识别的区域斜坡降雨入渗稳定性定量评价方法有效且具有很强的实用性，评价结果具有对实际滑坡的指示意义.

论文研究适用于以堆积层滑坡为主要地质灾害类型的地区，有助于此类地区区域堆积层斜坡灾害评价向精细化方向发展.

参考文献

[1] Baum, R. L., Savage, W. Z., Godt, J. W., 2002. TRIGRS: A Fortran Program for Transient Rainfall Infiltration and Grid-Based Regional Slope-Stability Analysis, Version 2.0. US Geological Survey, Boulder.

[2] Cascini, L., Calvello, M., Grimaldi, G. M., 2010. Groundwater Modeling for the Analysis of Active Slow-Moving Landslides. Journal of Geotechnical and Geoenvironmental Engineering, 136(9): 1220-1230. https://doi.org/10.1061/(asce)gt.1943-5606.0000323

[3] Catani, F., Segoni, S., Falorni, G., 2010. An Empirical Geomorphology-Based Approach to the Spatial Prediction of Soil Thickness at Catchment Scale. Water Resources Research, 46(5): W05508. https://doi.org/10.1029/2008wr007450

[4] Gao, K. C., Cui, P., Zhao, C. Y., et al., 2006. Landslide Hazard Evaluation of Wanzhou Based on GIS Information Value Method in the Three Gorges Reservoir. Chinese Journal of Rock Mechanics and Engineering, 25(5): 991-996 (in Chinese with English abstract).

[5] Hu, T., Fan, X., Wang, S., et al., 2020. Landslide Susceptibility Evaluation of Sinan County Using Logistics Regression Model and 3S Technology. Bulletin of Geological Science and Technology, 39(2): 113-121 (in Chinese with English abstract).

[6] Hu, Y., Li, D. Y., Meng, S. S., et al., 2020. Landslide Susceptibility Evaluation in Badong County Based on Weights of Evidence Method. Bulletin of Geological Science and Technology, 39(3): 187-194 (in Chinese with English abstract).

[7] Huang, F. M., Chen, J. W., Fan, X. M., et al., 2022. Logistic Regression Fitting of Rainfall-Induced Landslide Occurrence Probability and Continuous Landslide Hazard Prediction Modelling. Earth Science, 47(12): 4609-4628 (in Chinese with English abstract).

[8] Huang, W. B., Ding, M. T., Wang, D., et al., 2022. Evaluation of Landslide Susceptibility Based on Layer Adaptive Weighted Convolutional Neural Network Model along Sichuan-Tibet Traffic Corridor. Earth Science, 47(6): 2015-2030 (in Chinese with English abstract).

[9] Iverson, R. M., 2000. Landslide Triggering by Rain Infiltration. Water Resources Research, 36(7): 1897-1910. https://doi.org/10.1029/2000wr900090

[10] Jia, L., Cai, J. S., Yan, E. C., et al., 2021. Assessment of Landslide Susceptibility in Nanzhang County Based on Geological Environment Zoning. Yangtze River, 52(5): 86-94 (in Chinese with English abstract).

[11] Li, J. C., Ma, Q., Tao, J., 2012. Analysis and Application of Hydrological Basin Based on ArcGIS. Geospatial Information, 10(6): 121-123 (in Chinese with English abstract).

[12] Li, L. L., Teng, W. F., Zhang, Y. L., et al., 2019. Optimization of Shear Strength Parameters of Sliding Zone Soil in Huangtupo Landslide. Low Temperature Architecture Technology, 41(3): 75-79, 95 (in Chinese with English abstract).

[13] Liu, L., Yin, K. L., Wang, J. J., et al., 2016. Dynamic Evaluation of Regional Landslide Hazard Due to Rainfall: A Case Study in Wanzhou Central District, Three Gorges Reservoir. Chinese Journal of Rock Mechanics and Engineering, 35(3): 558-569 (in Chinese with English abstract).

[14] Miller, S. A., Degg, M., 2009. GIS Applied to Landslide Hazard Mapping and Evaluation in North-East Wales. EGU General Assembly Conference, Vienna.

[15] Montgomery, D. R., Dietrich, W. E., 1994. A Physically Based Model for the Topographic Control on Shallow Landsliding. Water Resources Research, 30(4): 1153-1171. https://doi.org/10.1029/93WR02979

[16] Pack, R. T., Tarboton, D. G., Goodwin, C. N., 1998. The SINMAP Approach to Terrain Stability Mapping. The 8th Congress of the International Association of Engineering Geology, Vancouver.

[17] Song, S. Y., Wang, Q., Pan, Y. Z., et al., 2014. Landslide Risk Assessment Based on Catastrophe Theory. Rock and Soil Mechanics, 35(S2): 422-428 (in Chinese with English abstract).

[18] Yang, C., Liu, L. L., Zhang, Y. L., et al., 2022. Machine Learning Based on Landslide Susceptibility Assessment with Bayesian Optimized the Hyperparameters. Bulletin of Geological Science and Technology, 41(2): 228-238 (in Chinese with English abstract).

[19] Yang, H., Tang, M. G., Ye, R. Q., et al., 2017. Study on Soil-Water Characteristics of Unsaturated Soil of Landslide Deposits in Three Gorges Reservoir Area. Water Resources and Hydropower Engineering, 48(9): 168-173 (in Chinese with English abstract).

[20] Zeng, L. F., Wu, Z. J., Chen, T., et al., 2012. Inversing Study of Overburden Thickness in Loess Area of Tianshui. Rock and Soil Mechanics. 33(6): 1912-1916 (in Chinese with English abstract).

[21] Zhu, C. H., Wang, Z. H., Liu, J. M., 2006. Study on the Relation between the Permeation Damage Slope and the Grain Composition of Coarse-Grained Soil. China Rural Water and Hydropower, (3): 72-74, 77 (in Chinese with English abstract).

[22] Zhuang, J. Q., Cui, P., Ge, Y. G., et al., 2010. Risk Assessment of Collapse and Landslide in “5.12” Wenchuan Earthquake: Taking Duwen Highway as an Example. Chinese Journal of Rock Mechanics and Engineering, 29(S2): 3735-3742 (in Chinese with English abstract).

基金资助

国家自然科学基金青年项目(41807264)

教育部重点实验室开放基金课题(2020KDZ01)

岩土钻掘与防护教育部工程研究中心开放课题基金(202102)

湖北省教育厅科学技术研究项目(B2019452)

中国地质大学(武汉)杰出人才培育基金项目(CUG170686)