国内外很多学者对滑（边）坡与降雨特征、地下水渗流关系展开研究。张玉成等^[1]认为滑坡发生明显滞后于降雨，但滞后期一般不超过10 d。房玉洁等^[2]、张勇等^[3]研究巴中地区降雨和滑坡之间的关系，认为滑坡灾害和强降雨类型密切相关。李绍红等^[4]等研究基岩型浅层边坡稳定性，认为降雨具有明显的时效性。张桂荣等^[5]通过模拟预测降雨及库水位联合作用下秭归八字门滑坡的稳定性，认为降雨量超过一定限制对滑坡影响较大且具有滞后性。林锋等^[6]研究贵州大榕滑坡形成机制后认为，多年来在堆积区南侧低洼地带持续堆填弃土及长时间较强降雨导致滑坡发生。郁舒阳等^[7]等讨论了不同非饱和土参数在不同雨型下对边坡渗流特性以及稳定性的影响，认为稳定性系数总体上随降雨持续走低。王雪冬等^[8]等研究降雨入渗特征，根据试验结果对经典Green-Ampt 入渗模型进行改进。Hou等^[9]分析汉源县中海村滑坡发生机制，认为是地下水上升诱发了滑坡。易武等^[10]应用Visual Modflow软件模拟滑坡渗流场，认为滑坡地下水位随库水位波动而变化。张杰等^[11]依据Green-Ampt模型对梯田边坡中农作物灌溉入渗进行稳定分析，认为梯田农业灌溉不利于边坡稳定。

四川盆地以东—大巴山褶皱带以西一带普遍分布侏罗系、白垩系红色碎屑沉积岩层，上覆崩坡积、残坡积含碎块石粉质黏土层，为典型红层地区。区内多为中低山地貌，人口密集，耕地资源紧张。受人类生产活动影响，研究区内斜坡地带多开垦为种植水稻的梯田，对斜坡地层结构、微地貌有不同程度的影响。红层地区滑坡比较发育，很多学者对红层滑坡成因机制进行了相关研究^[12-14]，均认为红层滑坡与降雨密切相关。但针对四川境内梯田分布的滑坡渗流分析及成因机理研究较少。

2021年9月28日四川省中巴市通江县新场镇七家沟村二社水头上滑坡（以下简称“水头上滑坡”）发生剧滑，造成10余处房屋、300 m水泥路损毁，44户149人紧急撤离。水头上滑坡9月1—5日遭遇3 d暴雨天气（日降雨量均超80 mm）至发生剧滑，历时近1个月，远超一般降雨诱发滑坡成灾的滞后时限。该区已进行过多次地质灾害排查和详细调查，采用遥感、InSAR、无人机航拍等新技术的地质灾害风险评价工作也已完成，但是并未发现该区域有何异常。因此，该类型滑坡具有高度复杂性和隐蔽性，对此进行研究十分必要和紧迫。故本文选择水头上滑坡为研究对象，通过野外调查、勘查钻孔、室内试验，结合卫星影像、无人机三维航拍、机载LiDAR等数据资料，对滑坡成因机制进行定性分析，将梯田作为单独地层建立滑坡地质模型，基于Fredlund & Xing土水特征曲线数学模型^[15]，采用Geo-Studio中的SEEP/W模块^[16]进行渗流分析，将不同时间阶段分析结果与SLOPE/W模块^[17]耦合，继而得到稳定系数与降雨、时间之间的动态变化关系。并将水头上滑坡还原至未开垦前，即将梯田土层与坡体含碎块石土视为同一地层作为对比，以研究开垦梯田对滑坡稳定性的影响。

1.   地质环境条件

通江县属亚热带秦巴区湿润季风气候类型，降雨集中在5—9月，多年平均降雨量1 119.3 mm，2021年1月1日至10月1日降雨总量为1 522.3 mm，其中9月降雨量为463.7 mm，是历史平均月降雨量164.5 mm的2.8倍。剧滑前日降雨量83.9 mm，9月日降雨分布如图1所示。

图  1  通江县2021年9月日降雨量分布图

Figure  1.  Distribution of daily rainfall in Tongjiang County in September, 2021

下载: 全尺寸图片 幻灯片

研究区位于通江河右岸缓坡中部，属侵蚀剥蚀-中切割中低山地貌，高程650～850 m。滑坡周边出露的地层除第四系松散堆积层之外，主要为中~上侏罗统（J_2-3）、下白垩统（K₁）砂泥岩互层（图2）。

图  2  研究区地质图

Figure  2.  Geological map of the study area

下载: 全尺寸图片 幻灯片

2.   数据来源

2.1   滑前数据

由于滑坡已发生剧滑，滑坡前高精度地形数据已无法获取。本次所用滑坡前地形数据是在2004年5月《通江土地利用现状图（1∶10000）》的基础上，采用滑坡前Google高程数据结合2016年4月16日Google高清卫星影像图、2017年5月12日Google高清卫星影像图及滑动后人工测量高程值进行内插修正获得。

本文收集了通江县气象站2021年9月1日—30日的降雨数据，主要用于滑坡的诱发因素分析。

2.2   滑坡发生后数据

（1）野外调查勘查

滑坡发生后笔者第一时间参与了灾害成因的现场调查分析工作，防范二次灾害发生。同时，技术人员迅速开展滑坡区基础地形、地质资料和历史卫星影像数据的收集以及现场勘探、测绘、室内试验等工作，通过各类数据的分析研究，为滑坡成因机理分析、变形趋势研判、治理方案选择提供支撑。

（2）无人机航拍影像

采用轻型无人机航拍数据精细建模，再通过数值模拟进行定量分析，最后与工程地质分析方法相互验证，是一条高效、实用、可靠的途径^[18]。

2021年9月29日，利用南方测绘公司SF700无人机对滑坡区进行航拍，航拍面积1.23 km²，共获取航拍照片368张。为了动态跟踪对比滑坡发展情况，2021年10月13日再次对滑坡区相同区域进行航拍。数据获取后采用高性能计算机处理计算，获得高分辨率数字地表模型（digital surface model，DSM）和数字正射影像图（digital ortho map，DOM）。

（3）无人机机载LiDAR

2021年10月13日，利用飞马公司机载D-LiDAR2000获取了滑坡区1.38 km² 高精度激光雷达点云数据，制作了大比例尺数字高程模型（digital elevation model，DEM）。通过LiDAR 数据与无人机航拍获取的三维数据进行精度比对，结合现场调查情况，得出无人机载LiDAR 数据在研究区获取的数据具有更高的精度。

2.3   调查结果

通过野外调查，基本查明了滑坡区下伏基岩岩性、结构面发育及产状特征、坡体地下水出露情况、滑坡影响区变形特征等。通过调查访问，掌握了滑坡前坡体变形历史、坡体人类耕种活动、建筑物变形开裂等相关信息。这些调查结果能为滑坡渗流分析及成因机理的定性判断提供依据。

3.   滑坡特征分析

3.1   坡体结构特征

水头上滑坡平面上总体呈箕状，平面图如图3（d）所示，主滑方向42°，主滑纵坡长375～408 m，宽240～320 m，高差最大118 m。通过现场调查，水头上滑坡区地层岩性为耕植土（${\mathrm{Q}}_4^{ml}$）、第四系崩坡积含碎石粉质黏土（${\mathrm{Q}}_4^{col+dl}$）、侏罗系上统蓬莱镇组（J₃p）粉砂质泥岩，岩层产状320°∠23°，强风化—中风化，上覆含碎块石粉质黏土、耕植土。通过主滑剖面6-6’^[19]揭示（图4），滑体主要由耕植土、含碎石粉质黏土组成，滑带土为软塑粉质黏土，滑床为粉砂质泥岩（图3）。耕植土层主要分布在梯田区域，主要成分为流—软塑状粉质黏土，厚0.3～0.7 m，弱透水性。梯田下部为含碎石粉质黏土，碎石含量15%～30%，厚2.7～12.0 m，稍湿—饱和，孔隙发育，中等透水性。强风化粉砂质泥岩，厚2.5～6.1 m，中等透水性。中风化粉砂质泥岩，揭露厚度4.7～14.8 m，微透水性。勘察发现坡体分布有地下水，为潜水，水位埋深0.7～9.5 m。

图  3  滑坡区影像

Figure  3.  Image of the landslide area

下载: 全尺寸图片 幻灯片

3.2   滑坡变形破坏特征

图  4  滑坡区工程地质剖面图（6-6’剖面）

Figure  4.  Engineering geological profile of the landslide area (6-6’ section)

下载: 全尺寸图片 幻灯片

通过对滑坡区多期卫星影像、无人机倾斜摄影、机载LiDAR等数据处理分析，结合现场调查、勘查、取样测试分析等方法，对滑坡特征有了清晰的认识。滑坡区总体上分为主滑区（A区）、滑坡影响区（B区）、中部弱变形区（C区）3个区域（图3）。

（1）主滑区（A区）特征

根据坡体变形破坏特征，结合现场调查访问、工程钻探、滑坡区各类影像资料，主滑区分为启动区（Ⅰ）、滑移区（Ⅱ）、侧向滑移区（Ⅲ），见图3。启动区（Ⅰ）位于滑坡体中上部，高程745～795 m，下部剪出口高程 745～760 m。启动区纵向平均长140 m，宽 245 m，厚 9 m。启动区坡度17°～39°，坡向30°，主滑方向25°，岩性如前文所述。启动区边界条件和空间特征十分清晰。滑体以岩层层面作为底滑面，以含碎块石粉质黏土层底部软塑粉质黏土为滑带土，以坡体侧面下错陡坎或剪切变形裂缝作为两侧边界，滑坡后缘下错明显，形成1.5～3.0 m滑坡壁，如图5（a）所示。该区域自2021年9月以来先后经历多次强降雨天气，一部分降水渗入坡体，坡体地下水位升高，坡体表面有地下水渗出，如图5（b）所示，土体有效重度、强度降低，为启动区进一步变形创造了有利条件。滑动后，于启动区下部形成总体WN走向的弧形鼓丘。

图  5  启动区滑坡壁与地下水特征

Figure  5.  Characteristics of the landslide wall and groundwater in the start-up area

下载: 全尺寸图片 幻灯片

滑移区（Ⅱ）位于滑坡体中下部，坡向60°，高程650～760 m，以整体式滑移为主要特征。当启动区失稳下滑后，以巨大的动能推挤下部斜坡表面原有松散堆积物。下部斜坡主要为梯田（图6），5—9月封闭蓄水种植。由于地表水、连续强降雨下渗，使得坡体已处于蠕滑状态，但此时滑移面尚未全面贯通。在受到启动区下滑推挤作用之下，滑移面全面贯通，滑移区向下滑动。将无人机正摄影像与滑动前卫星影像进行对比（图4），发现主滑区滑动距离40～48 m，主滑方向42°，滑坡剪出口在水泥路下方。由于滑移区下滑挤压，在本区下部发育总体ESE走向的鼓胀区，造成道路开裂变形，农田、房屋受损严重（图7）。另外，随着滑移区急速下滑，与启动区之间形成一条弧形错断边界，该次级滑坡壁高1.5～7.5 m（图8），壁高随高程逐渐降低。本区地下水十分丰富，坡体各部位均有泉眼出露，并在坡体形成径流，汇集于道路内侧排水沟。滑移区中下段西侧低、东侧高，导致滑动方向向西偏转，这一特点与坡体中地下水渗流方向一致。

图  6  滑坡区内梯田（Google卫星影像图，2016年4月16日）

Figure  6.  Terraced fields in the landslide area (Google satellite image, April 16, 2016)

下载: 全尺寸图片 幻灯片

图  7  受损的房屋道路

Figure  7.  Damaged houses and road

下载: 全尺寸图片 幻灯片

图  8  滑移区次级滑坡壁

Figure  8.  Secondary scarp of the landslide

下载: 全尺寸图片 幻灯片

（2）滑坡影响区（B区）

滑坡影响区（B区）如图3所示。滑坡发生后，通过现场调查及无人机高清三维模型影像，发现受主滑体高速滑动和运动过程中牵引、拖拽作用的影响，在滑坡区上坡一侧形成一条总体走向ESE、长220.0 m，宽15.0～20.0 m的牵引区。牵引区分布大小拉张裂缝20余条。拉张裂缝由坡外向坡内依次下错，方向与滑坡边界走向一致。经笔者现场测量，裂缝长7.5～125.0 m，宽7.0～25.0 cm，深25.0～50.0 cm，无地下水充填。

（3）中部弱变形区

中部弱变形区（C区，见图3）为滑移区（Ⅱ）与侧向滑移区（Ⅲ）中间一带的三角形区域，横向宽50 m，纵向长85 m，高程693～728 m，坡度20°～25°，坡向62°。通过钻探、探槽揭露，覆盖层厚度2.0～3.0 m，基岩埋深普遍较浅。该区域与侧向滑移区（Ⅲ）处于同一斜坡，与滑移区（Ⅱ）斜坡交界处形成NW向山脊。该区域在水头上滑坡发生前后均处于基本稳定状态，调查过程中未发现明显变形迹象。

3.3   滑坡形成过程分析

通过现场调查，结合多期卫星影像分析，水头上滑坡长期处于稳定状态。但是，有利的地形地貌、地层岩性为水头上滑坡发生的主要控制因素，连续性降雨为诱发因素。滑坡位于斜坡中部，后缘堆积层厚，下部滑带土（粉质黏土）强度低，梯田陡坎为滑坡剪出提供了良好的临空条件。所在斜坡地段坡度较缓，地层渗透性一般，加之地表为梯田（水田），排水条件差。2021年9月以来，通江县遭遇罕见连续强降雨天气（图1），导致坡体地下水位上升，稳定性持续下降。根据滑坡变形特征，推断水头上滑坡形成经历了蠕滑、滑动、剧滑、牵引4个阶段：

启动区蠕滑-滑动推移挤压阶段：受连续强降雨影响，雨水下渗及地下水渗流导致启动区地下水位逐渐升高，孔隙水压力持续增大，土体有效重度降低，坡体滑动面逐渐贯通，在月末强降雨诱发作用下，由蠕滑发展到滑动阶段，推动挤压前方土体。通过对滑动前遥感影像进行分析，在启动区下部存在一平台，推断该平台起到一定缓冲作用，因此并未持续下滑。

滑移区蠕滑-滑动-剧滑阶段：处于蠕滑阶段的滑移区，在受到上部土体推移挤压和月末强降雨综合作用下，突然失稳滑动。由于下方缺乏有效堆积支撑平台，滑动段持续运动至坡体下部，形成鼓丘状堆积体，才基本稳定下来。

主滑区牵引拉裂阶段：滑移区剧滑再次激活启动区，启动区继续向下滑动，牵引、拖拽主滑区后缘坡体，在影响范围内形成数量、规模可观的拉张裂缝。

4.   数值分析与耦合

4.1   Fredlund & Xing模型

国内外学者研究表明，对于非饱和多孔介质，其渗透系数为饱和度和体积含水率的函数，由于体积含水率和饱和度与基质吸力之间的关系可以用土水特征曲线来体现，渗透系数也是基质吸力的函数。受土水特征曲线迟滞效应的影响，合理应用土水特征曲线相当困难^[20]。Fredlund等^[15]通过对土体孔径分布曲线的研究，用统计分析理论推导出适用于全吸力范围的任何土类的土水特征曲线的表达式，即Fredlund & Xing模型：

式中：θ_w——体积含水率/%；

θ_s——饱和体积含水率/%；

$C_{\varphi }$——函数的修正函数；

e——常数2.71828；

$\varphi$——负孔隙水压力/kPa；

α、m、n——拟合参数，通常α、m、n的值可用拟 合算法和测量数据点的方式获得^[16]。

本文根据Geo-Studio提供的样本函数曲线^[17]和岩土体饱和状态参数^[19]来拟合确定非饱和岩土体渗透函数和土水特征曲线，其中α、m、n取值为α=100 kPa，m=2，n=1。

利用式（1），Fredlund等提出了渗透系数计算方法，即沿整个体积含水率函数进行积分^[21]，其计算公式为：

式中：k_w——负的孔隙水压力对应的渗透系数/（m·d⁻¹）；

k_s——饱和渗透系数/（m·d⁻¹）；

y——虚拟变量；

i——j到N之间的数值间距；

j、N——最小、最大负孔隙水压力/kPa；

θ——第j步的负孔隙水压力/kPa，θ′为方程的起 始值。

4.2   数值模型的建立

Geo-Studio2012软件中的SEEP/W模块^[16]是一种专门用于岩土渗流计算的有限元模块，可以用来还原滑坡发生前后的动态水文响应过程，该模块中已嵌入式（1）。本文采用该模型计算饱和-非饱和土体动态含水率，将结果用于式（2），得到各时步坡体渗流分析运算结果。

6-6’剖面（图4）在水头上滑坡各阶段起决定作用，因此，以该剖面为基础建立模型。假定梯田为粉质黏土层，均匀分布于坡体表面，取平均厚度0.5 m，其下分别为含碎块石粉质黏土、滑带土（粉质黏土）、强风化粉砂质泥岩、中风化粉砂质泥岩。勘察钻孔揭示启动区段滑带土与滑移区段强度存在差异，在此分别定义为滑带土①、滑带土②，以示区别。另外，根据基岩产状与斜坡坡面产状关系、基岩岩体质量等级，推定滑带土以下基岩为不滑动层，其他地层为满足摩尔-库仑（Mohr-Coulomb）本构模型的弹塑性材料，各地层参数见表1。数值模型采用四边形和三角形网格进行有限元网格划分，单元尺寸2.5 m、5.0 m，共计划分9799个节点，3 247个单元（图9）。

图  9  模拟数值模型与有限元网格划分

Figure  9.  Simulation numerical model and finite element meshing

下载: 全尺寸图片 幻灯片

4.3   坡体渗流数值分析

水头上滑坡前的水文地质资料已无从考证。但是滑动后，通过在滑坡区地质钻孔进行稳定水位观测，结合地区经验复原滑动前坡体地下水稳定水位作为初始状态也是一种可行的方法。基于此，在SEEP/W模块中设定瞬态渗流边界条件：以地表坡面为入渗边界，图1 中的日降雨量曲线作为入渗量曲线；以坡体初始水位线作为模型初始水头（图9虚线所示）；模型左边界、右边界和底边界均设置为透水边界，以真实反映岩土材料性质。模型中分别于启动区（图9）设置4个监测点（A1—A4），用于观察孔隙水压力、材料渗透系数变化。数值模型中各层岩土材料物理力学参数按表1取值。模拟以1 d作为计算步长，计算周期为2021年9月1—27日，共计27个时步。

表  1  计算模型中岩土材料物理力学参数取值

Table  1.  Values of physical and mechanical parameters of geotechnical materials in the calculation model

岩土材料	重度 /（kN·m−3）	黏聚力 /kPa	内摩擦角 /（°）	饱和体积 含水率/%	饱和渗透系数 /（m·d−1）
耕植土（粉质黏土）	20.5/20.5	8	10	60	0.05
含碎块石粉质黏土	20.7/22.0	25	21	45	1.52
滑带土①	20.3/21.8	11	15	45	0.05
滑带土②	20.3/21.8	20	16	45	0.05
强风化粉砂质泥岩	22.5/23.4	30	27	20	1.24
中风化粉砂质泥岩	23.7/24.6	125	30	20	0.001
注：“/”号前后分别为天然重度和饱和重度。

下载: 导出CSV 

| 显示表格

通过瞬态叠加运算，得到每一时步地下水渗流场分布，这里选取第0天、第27天的坡体地下水渗流场进行对比分析（图10）。通过对比坡体初始地下水渗流场和滑坡前地下水渗流场发现滑坡前0孔隙水压力等势线明显升高，在坡体多处位置与地面线相交，说明坡体中有地下水渗出，与现场调查结果一致。通过位于非饱和土体（含碎块石粉质黏土）中的监测点A2，分析滑坡前坡体材料属性、孔隙水压力的动态响应过程。如图11（a）所示，监测点A2含水率变化经历了3个上升阶段：降雨开始前5天，随着雨水入渗和9月4—5日强降雨的发生，体积含水率由38.9%急剧上升至43.1%，这个阶段降雨量最大，对含水率影响也最大；随着降雨停歇，含水率趋于稳定，在中下旬随着降雨入渗，含水率再次上升并达到饱和含水率；在此期间，渗透系数也随之增大至1.5 m/d，如图11（b）所示。A1—A4监测点孔隙水压力监测结果显示，坡体内同一截面不同深度，其孔隙水压力总体变化基本一致，即随降雨入渗孔隙水压力增大。雨停后，随着地下水渗流补给量逐渐减小至小于坡体排泄量，地下水位下降，孔隙水压力降低。强降雨结束半个月以后，坡体仍有泉眼出露，剪出口出水量也未见明显减少，这和计算结果推论一致。降雨对浅层碎石土滑坡的作用主要通过增加滑坡体的孔隙水压力，使滑体的下滑力增大^[22]。另外，研究发现滑移区土体材料属性、孔隙水压力动态响应过程与启动区类似，受篇幅所限不再单独说明。

图  10  滑坡前后地下水渗流分析对比

Figure  10.  Analysis and comparison of groundwater seepage before and after landslide

下载: 全尺寸图片 幻灯片

图  11  滑坡土体水文监测结果（监测点A2）

Figure  11.  Hydrological monitoring results in the landslide soil (monitor point A2)

下载: 全尺寸图片 幻灯片

4.4   坡体稳定性及滑坡过程分析

将SEEP/W模块每一时间步计算得到的渗流分析结果耦合到SLOPE/W模块^[16]，利用Morgenstern-Price 方法计算滑坡的稳定性系数变化。SLOPE/W模块模拟计算模型中各岩土材料参数按表1取值。滑面为指定滑面（图9）。如图12（a）所示，启动区的稳定性系数初值为1.198，稳定性系数在降雨后均有降低趋势；在9月中下旬连续降雨后，稳定性系数降至1.008，启动区坡体已经处于蠕滑阶段；在27日降雨后，稳定性系数降至0.97，启动区失稳滑动。滑移区稳定系数初值为1.205，稳定性系数走势与启动区一致，但在27日降雨后，稳定系数1.026，处于蠕动挤压阶段。通过模拟计算结果显示，启动区率先失稳下滑，推移挤压前方土体，触发了本已临界滑动的滑移区，大规模滑动随后展开，直到在剪出口附近堆积趋稳，这与滑坡形成过程定性分析结果基本一致。

图  12  稳定性系数变化

Figure  12.  Variation of stability coefficient with time

下载: 全尺寸图片 幻灯片

为了进一步研究梯田（水田）对坡体稳定性的影响，将坡体复原后渗流分析结果与坡体复原后SLOPE/W模块进行耦合，稳定性计算结果如图12（b）所示，启动区最小稳定性系数1.001，滑移区最小稳定性系数1.043。复原后稳定性系数计算结果均比现状坡体大，说明在该地区斜坡开垦梯田（水田），改变了坡体原生渗流环境，在一定程度上降低了坡体稳定性。

5.   结论与建议

（1）川东红层梯田水头上滑坡为推移式滑坡，连续性降雨导致坡体地下水位逐步升高，稳定性逐步降低，前期连续性降雨和临近暴雨导致滑坡发生。

（2）把梯田纳入进来进行现状坡体渗流分析，其结果与稳定性分析耦合与复原后坡体稳定性分析结果对比显示，该地区斜坡梯田（水田）对坡体稳定性有不利影响。在梯田斜坡稳定性分析中采用这种方法具有一定的创新性。

（3）滑坡渗流分析和稳定性分析结果与坡体勘察、调查、定性分析结论基本吻合，说明基于Fredlund & Xing模型的有限元渗流分析在川东红层梯田滑坡中有很强的适用性。

降雨入渗是水头上滑坡发生的主因，导致不论是专业人员地面调查的传统方法，还是遥感、InSAR等新技术都难以预报，这是个世界性难题。在此，笔者认为可以从以下几个方面做好防灾减灾工作：使用ArcGIS导入精度大于1∶10 000的DEM，结合现场调查，对斜坡单元按坡度、坡向、地层岩性等孕灾条件进行筛选确认，实行分类地质灾害防控，即对潜在危害性大的斜坡，新建或改扩建坡面截排水设施，将降水阻挡或及时排出；对于潜在危害性斜坡，加强汛期人工巡视排查。

另外，通过本文研究，将水田改为旱田有助于坡面排水，对坡体稳定有利，但还需要加强这方面的研究讨论。

本文受数值模拟软件、研究经费和时间等因素影响，并未考虑上方坡体汇水、坡体自然排泄、滑带土体浸泡软化、大气降水蒸发等因素对坡体渗流分析和稳定性分析的影响，可能影响分析的结果准确性，后续的工作研究中持续跟踪这方面研究进展。