基于PFC3D位移判据的低渗堆积层滑坡临界降雨阈值反演

左双英，蒲 泉，史文兵，王 勇，吴道勇，王安礼

(1.贵州大学 资源与环境工程学院，贵州 贵阳 550025； 2.贵州大学 喀斯特地质资源与环境教育部重点实验室，贵州 贵阳 550025；3.贵州省有色金属和核工业地质勘查局物化探总队，贵州 都匀 558099； 4.贵州省质安交通工程监控检测中心有限责任公司，贵州 贵阳 550081)

贵州省地处云贵高原的斜坡地带，受河流的侵蚀切割，地势起伏大、地表相当破碎，地质灾害频繁。随着城镇化建设的步伐加快，对山地斜坡的开发利用加快，加之贵州夏季多雨的气候条件，各类崩滑事故愈加频繁。近年来，关岭岗乌[1]、福泉[2]、纳雍张家湾[3]、水城鸡场[4]等几起重大地质灾害事故对当地民生危害巨大。上述崩滑事故中，降雨诱发的堆积层滑坡是其中最主要的滑坡类型之一，其特点是堆积体表层土粒径粗、孔隙大、吸水性强，下部为细颗粒粉土或黏土，渗透性低，在短时强降雨作用下迅速吸水入渗，在坡体内部大量聚集来不及排出，导致坡体孔隙水压力增大，有效应力、抗剪强度降低[5]，破坏了原有的平衡状态而产生滑坡。因此，降雨型堆积层滑坡坡体物质组成与降雨入渗地下水复杂作用的特殊性给斜坡传统力学评价预测方法的有效应用带来了极大的困难与不确定性[6-7]。为此，深入研究滑坡的先兆因素，得出易于监测的指标，提高滑坡预警预报的可操作性和准确性是业内亟待解决的重点和关键问题。

滑坡的预警研究自日本学者斋藤率先提出的经验公式至今已有近50余年的发展历程[8]。期间，国内外众多学者针对这一问题展开系列研究，对滑坡的预警判据研究做出了重要贡献。例如，适用于长期预报的判据[9]：稳定性系数(K)、声发射参数以及可靠概率(Ps)；适用于临滑预报的判据[10]：位移加速度、位移矢量角、蠕变曲线切线[11]等。以上判据的提出对于我们在滑坡的理论研究上有重要意义，但是对于实际的预警预报工作却有诸多不便，例如指标不易观测、经济成本高等。许强等[12]强调应加强区域性气象预警和单体滑坡预警工作，将预警的重心转移到对历史数据的统计分析和基于变形、地下水位、雨量等关键指标的预警模型和判据研究。Shruti Naidu等[13]利用事件进行回归分析，建立降雨阈值方程，通过降雨阈值分析滑坡稳定性的Fs值建立安博里滑坡预警系统，但是该法仅对浅层滑坡适用；陆研等[14]利用Liang Kleeman信息流方法研究得出重庆地区滑坡泥石流发生频次与年均降雨量显著相关，建议在暴雨发生时加强对滑坡灾害的监测与预报；吴益平等[15]基于信息—物元模型将传统的滑坡灾害危险性预警延伸到滑坡灾害风险预警预报，为降雨诱发滑坡的预警工作提出了新的思考。综上所述，目前国内外关于降雨型滑坡的预警研究大多属于大面积区域上的预报，对于单体滑坡的降雨参数预警却很少涉及。

文中依托贵州省地质灾害自动化监测预警项目对该坡体展开深入研究，并基于坡体表面堆积松散、内部兼具低渗透性的特性对其预警预报工作开展研究。通过现场调查和无人机图像(UAV)结合数字高程模型(DEM)建立该滑坡的PFC3D数值模型，模拟该坡体在自重和不同地下水位静水压力作用下的变形和运动过程，监测模拟过程中颗粒的位移和速度，将位移突变点(计算不收敛)对应的地下水位作为临界参数；在Pradel & Raad提出的降雨入渗土坡深度的算法基础上，反演临界降雨时长以及降雨强度参数，并验证降雨预警阈值的合理性。

1 研究区工程地质概况

研究区表层堆积体结构松散，部分降水透过松散土层入渗至薄层泥岩风化物。泥岩风化物级配曲线如图1所示，土体颗粒主要为细粒组和砂粒组，粘粒含量约为21%，粉粒约占32%，砂粒组含量47%，基本符合含砂粘质粉土的粒径组成。另外，土体曲率系数Cc=2.63，不均匀系数Cu=64.2为级配良好型土体，其颗粒级配较好在宏观上表现为较小颗粒能够很好的填充在空隙当中，表现出相对较低的渗透能力。

滑坡后缘、中部不断有工程活动建设，各种填方堆载加大了斜坡下滑力(图2)。滑坡初现变形时间为2008年，于2013年4月开始，受强降雨作用影响，滑坡体上的房屋、道路、地面陆续出现开裂变形、沉降等变化。在多雨季节由于上部土体松散，降水入渗坡体内部；而坡体底部土体渗透系数低，造成水体易进难出，使得斜坡内部地下水位线持续升高，自重增加且孔隙水压力上升，坡体的稳定性朝着不利的方向发展(图3)。从2017年10月25日第一次监测开始，测点位移量较前期监测有增大趋势，滑体坡面上有多条新裂缝，据此推断山体蠕动变形速度具有加剧趋势，稳定性较差，发生滑坡的可能性增大，严重威胁当地居民的生命财产安全。

图3 纳缝滑坡工程地质剖面图(3-3剖面)Fig.3 Engineering geological profile of Nafeng landslide

2 离散元数值模拟

2.1 颗粒流计算原理

近年来，颗粒流(PFC3D)作为离散元软件，已经成为边坡工程、地下洞室开挖、路基稳定等岩土工程领域中十分成熟的数值计算方法[16]。区别于有限元，颗粒流中的球体单元较为符合散体材料的特性，在涉及土质边坡的计算中表现出一定的优越性，因此被众多学者沿用[17]。球体之间的接触包括ball-ball接触和ball-wall接触。球体受到重力以及外力的作用失去原有的平衡，在颗粒的动态变化中通过牛顿第二定律计算其位移变化，通过不断更新球体位移计算模型的运动情况。

2.2 微观参数的获取

在离散元的计算中，散体材料的宏观力学特征通过微观强度表现。但是，并没有合适的函数关系可以合理地表示散体材料的宏观参数与微观强度之间的非线性关系。通常情况下，通过匹配数值试验与室内试验的应力-应变曲线获得[18]。结合该研究区的勘查资料以及现场取样，选择室内大型剪切试验进行力学测试(图4)，不同法向力作用下的剪应力-剪切位移曲线见图5，所得宏观力学参数如表1，其中渗透系数是通过现场双环渗透试验测得。

图4 试验装置 图5 不同法向力下的剪应力-应变曲线

表1 土样宏观物理力学参数Table 1 Macro mechanical parameters of soil samples

为了更好的模拟土体中力的作用，选择接触粘结模型。试验所得粘质粉土细颗粒粒径小于0.05 mm，在模拟过程中，数值试验的颗粒粒径值选取很难完全与真实的颗粒粒径保持一致，为了满足计算效率需要，等比例扩大半径生成颗粒。模型尺寸25×25×25 m，颗粒最大最小半径比为1.66，通过内嵌fish语言伺服控制，实现大剪试验三维数值模拟(图6)。模拟过程中反复调整相应微观参数，直到应力-应变曲线(图7)逼近室内试验结果。数值试验所得微观参数如表2所示。

图6 剪切盒与数值试验结果Fig.6 Shear box and shear results

图7 大剪数值试验剪应力-应变曲线Fig.7 Shear stress-strain curve of large shear numerical test

表2 滑坡模型微观参数Table 2 Micro-strength parameters of the landslide PFC model

2.3 纳缝滑坡(PFC3D)模型建立

纳缝滑坡数值模型主要包括地形高程数据获取以及几何模型的生成。地形数字高程通过现场无人机勘测摄像(UAV)生成，利用该数字高程生成几何模型(图8)，模型的相关参数为已获得的微观强度参数。

图8 纳缝滑坡(PFC3D) 数值模型图 图9 滑体地下水静水压力计算示意图

随着降雨时长的增加，入渗到坡体内部的水也会随之增多，而滑面以上的土层透水性能差致使这些水体不断聚集抬高地下水位线，该处的静水压力值将不断增大，在静水压力的作用下滑面滑动的趋势增强，当地下水位线增加到一定程度时，土体内部的抗剪强度将无法承受变大的静水压力值而发生破坏(图9)。基于此，该模型在计算过程中，主要通过改变颗粒不同静水压力大小而全面分析滑坡各个位置的变形情况，模拟过程中分别对滑坡前缘、后缘以及中部9个区域进行全过程监测。

2.4 模拟结果

模型在代入匹配的微观参数后到达最初的平衡，此时对每个颗粒施加垂直向下的静水压力(p=ρhg)，每个颗粒在静水压力平行于坡面的分力作用下开始发生变形。通过反复调试静水压力的大小，找到模型位移不收敛的静水压力区间。随着降雨入渗，土体颗粒的静水压力增大是一个动态增加的过程。当静水压力值在某个界限值以下时，土体颗粒最终会趋于平衡。在该模型的计算过程中，增加颗粒的静水压力值，颗粒的位移曲线会出现阶梯状态的增长，但是随着静水压力值的增大，其平衡的程度会降低，具体表现为平衡阶段曲线的水平程度降低。静水压力值增大到某一个界限值时，颗粒之间的力链联结发生断裂时，则认为坡体发生了破坏。

静水位高度h=1.60 m时坡体监测点的变形如图10所示。在静水压力的作用下，坡体各个部位均发生不同程度的变形，变形速率先增加随后减小最终趋于零，说明位移收敛，坡体趋于平稳。由图可知，监测点1、监测点4附近位移较大，其峰值分别为0.32 m和0.21 m；其余监测点位置变形相对较小，位移变化在0.02 m-0.09 m之间。结合监测点所对应位置分析，坡体左侧的稳定程度相对较低，当外部条件发生改变，其变形破坏迹象明显；而监测点3、6、9位移变化较小。

图10 静水压力高度1.60m时监测点位移-时步曲线Fig.10 Displacement-time curve of monitoring point at 1.60 m hydrostatic pressure

当静水位高度增加到1.70 m时，随着颗粒受到静水压力的增大，变形也呈现出增大的趋势(图11)，此时各个监测点位置位移增长较为明显。但是在整个坡面区域，位移变化却存在较大的差异性，从各个监测点的位移变化来看，监测点1、4所在位置为整个坡面最不稳定区域，监测点3稳定程度次之，监测点2稳定程度最高。整个坡面的后缘位移最大、中部以及前缘位移依次递减，符合推移式滑坡特点。在该静水压力条件下，整个坡面的位移最终也是趋于稳定，但是稳定程度明显降低，当外营力增加时极有可能诱发滑坡。

图11 静水压力高度1.70 m时监测点位移-时步曲线Fig.11 Displacement-time curve of monitoring point at 1.70 m hydrostatic pressure

当静水位高度继续增加到1.75 m时，各个监测点的位移明显变化，局部位置的位移表现出不收敛的变化特性，如图12中监测点1、3、4。总体来看，滑坡的失稳发生在其后缘、中部位置，而滑坡前缘的位移变化均能趋于稳定(图13)。分析认为，这是由于上部颗粒在水压力作用下向下滑落，在前缘位置堆积所致。在该静水压力条件下颗粒之间的力链联结主要表现为拉应力和压应力，数值模型破坏时的力链分布如图14所示。图中蓝色代表压应力，绿色代表拉应力，剖面3-3’监测点1、4、7所对应位置颗粒之间的力链主要为压应力，局部位置为拉应力，斜坡表现为整体下滑；而剖面1-1’监测点3、6、9所在剖面在后缘位置有着明显的张拉破坏，土体颗粒向下运动出现颗粒的堆积；在滑坡的中部位置以土体的剪切破坏为主，有明显的剪切裂缝产生，并伴有部分张拉-剪切复合破坏，破坏沿着垂直于坡面方向逐渐发展贯通程度逐渐加大，最后裂缝贯通的区域出现整体的下滑并进一步挤压前缘位置颗粒；监测点2、5、8的位移变化情况可以看出斜坡2-2’剖面在该静水压力条件下的位移是收敛的，最终仍然趋于稳定，在其中部颗粒最大位移是25.4 cm，分析认为可能是由于1-1’剖面的张拉-剪切复合破坏扩展发展所致，但是并未形成贯通性的破坏，总体表现出较稳定的状态。

图12 静水压力高度1.75 m时监测点位移-时步曲线Fig.12 Displacement-time curve of monitoring point at 1.75 m hydrostatic pressure

图13 位移云图Fig.13 Displacement cloud map

图14 力链分布图Fig.14 Force chain distribution

以上过程表明，纳缝堆积体滑坡稳定性与入渗坡体形成的静水压力大小关系紧密。经过计算，静水位高度为1.60 m时，坡体位移最终会处于稳定状态；当静水位高度增加到1.70 m时，坡体位移曲线平缓程度下降；继续增加静水位高度到达1.75 m时，部分位移不能收敛，其强变形区为N01、N03以及N04。通过现场无人机勘测摄像(UAV)生成数字高程模型(DEM),建立的纳缝滑坡PFC3D数值模型能够全方位、多角度地反应纳缝滑坡的变形破坏过程，其位移轨迹能够较好的吻合现场监测影像。相比较PFC2D，PFC3D具有明显的优越性，更加真实地反应滑坡过程中各个部位变形的强弱，更有利于在预警工作中对强变形区的识别和判断，做出合理的预警决策。

3 滑坡预警阈值

3.1 反演模型

纳缝堆积体滑坡具有低渗透特性，在降雨过程中斜坡表面的土体含有不同粒径的碎石，因此表现出一定的松散特性，加之在不断的蠕动变形过程中斜坡各部位尤其是后缘发育有多条裂隙，这为降雨入渗提供了便利的渗流通道，更多的降水通过这些通道进入到坡体内部。蠕动变形当中再次松动的是泥岩风化土体，其细粒组含量达到53%，大小颗粒之间能够互相填充表现出良好的级配特性，渗透系数为3.01×10-5cm/s，透水性能较差。因此，入渗到土体内部的水体最终会在该处富集提高地下水位线，抬高的水位线高度可近似与入渗的深度相同，水体的富集会使孔隙水压力上升有效应力降低，同时水位线以下部分土体的静水压力会逐渐升高从而诱发滑坡。

降雨对滑坡稳定性的影响可以分为以下2种情况：一是当降雨强度小于岩土体入渗强度时，坡面大部分降水通过表面径流作用排泄，少部分降水通过表面入渗作用浸入岩土体内部。此时岩土体内部必须考虑非饱和渗流问题，根据相应水文地质、工程地质条件，综合考虑其他相关因素，如渗透性、植被、地下水位变化等情况；二是当降雨强度大于岩土体入渗强度时，水通过已有的渗流通道以及表面的入渗作用能够不断浸润坡体，增加土体自重从而诱发滑坡。上覆土层的饱和区高度受降雨强度影响，并随着降雨时长增加而升高。本文仅针对后者，讨论当降雨强度大于岩土体入渗强度时，通过离散元数值模拟所确定的临界静水位高度进而推算出降雨强度以及降雨时长，作为纳缝滑坡的预警参数阈值，可为纳缝滑坡自动化预警提供可用依据。

强降雨工况下，土体入渗深度的推算需要满足的2个基本条件是降雨强度Imin大于土体入渗率v0并满足一定降雨时长(T)。目前针对降雨土坡入渗深度的计算方法较为合理的是国外学者Pradel & Raad[19]提出的经验算法。本文依托离散元计算出的静水位高度，据此公式反演出相对应的临界降雨强度以及降雨时长，并计算出降雨总量。计算方法如下所示：

(1)

(2)

(3)

式中，v0为入渗率；k为渗透系数(cm/s)；h为入渗深度(cm)；t为降雨时长(s)；S为基质吸力(cm)；θs为土体饱和含水率(%)；θ0为土体实测含水率(%)；ds为相对密度，e为孔隙比。

基质吸力通过粘质粉土的土水特征曲线[20]得到，计算方法如下:

θ0×100=-0.0002S2-0.0115S+46.04.

(4)

式(3)反映了土体体积含水率与基质吸引力的动态变化关系。本文研究土体在降雨过程中的临界阈值问题，随着降水增多土体含水量增加，因此在使用式(3)的时候需要代入土体破坏时的含水量。同时，在降雨过程中土体水量增加土体颗粒有一个从塑限到液限的转变过程。当土体颗粒趋于液限时，土体就由可塑态变为流动状态，土体也就丧失了其应有的强度。基于以上两点原因，在通过式(3)推导其基质吸力时，θ0的取值为土体液限值，即θ0=wL=39.5%。

将θ0=39.5%代入到式(4)求解方程组得：

基质吸力S=154.352 cm

降雨强度：

降雨时长：

累积总降雨量Q=I×T=5.73×10-5×8.87×24×3600=439.13mm

通过以上计算得出纳缝滑坡临界降雨强度I为5.73×10-5cm/s，降雨时长T为8.87 d，累计总降雨量为439.13 mm。

3.2 滑坡预警阈值指标拟定

根据滑坡预警指标划分方式采用阈值的60%、70%、80%和90%进行预警指标划分[21]，分别对应蓝色预警、黄色预警、橙色预警和红色预警，各等级判别的临界位移见表3，降雨强度和降雨总量见表4。

表3 纳缝滑坡位移预警指标Table 3 Early warning of displacement of Nafeng landslide cm

表4 纳缝滑坡降雨预警指标Table 4 Rainfall early warning of Nafeng landslide

3.3 实测资料验证

罗甸县气象局2008年至2017年10年降雨量统计，年最大降雨量1 326.2 mm(2014年)，年最小降雨量846.1 mm(2009年)，年平均降雨量1 098.68 mm。罗甸县月平均降水量如表5所示，从降水量的时间分配来看，5月到10月是多雨的季节，集中了全年降水的89%。

表5 研究区2008-2017年每月平均降雨量Table 5 Average monthly rainfall in the study area, 2008-2017

研究区于2017年10月25日开始布设观测点等基础工作，并于10月25日进行第一次观测点测量。截止2018年9月29日，对该滑坡体进行了20次监测，根据各监测点情况，最大累计水平位移量为310.3 mm(图15)，对应剖面3-3’所在区域，其次是1-1’剖面区域，变化最小的是2-2’区域；从整个坡体来看，位移变化量最大的地方是其后缘公路所在位置，坡脚处位移相对较小。

通过数值计算的结果(图12)与实测数据(图15)以及近10年的降雨数据(表5)分析发现：斜坡在2017年10月到2018年5月的位移变化较小，而近十年的月平均降水最少的时间区间也是10月到次年5月，在降水充沛的月份里位移的变化则是相当明显的，这说明了导致该研究区变形破坏的最主要因素是降雨。从表3拟定的滑坡位移预警指标分析来看，监测期间最大位移为3-3’剖面附近的CJ10处310.3 mm，1-1’剖面最大位移为CJ6处190.5 mm，监测值与计算值接近，并且达到蓝色预警等级，证明了离散元数值计算得出的临界位移变化量的合理性；由表5可知，罗甸县每年5月、6月为全年降雨量最大季节，平均降雨量分别

图15 位移监测结果Fig.15 Displacement monitoring results

为213.27 mm，204.73 mm。而滑坡初现时间2008年5月的降雨量为288.9 mm达到了雨量的蓝色预警值；研究区降雨充沛，为该区的低速缓慢活动提供了不利条件。这种低速缓慢运动使得内部岩土体的强度逐渐变弱，在持续的强降雨条件下都有可能诱发滑坡。2017年8月降雨量为246.4 mm，基本达到蓝色预警雨量，而监测位移结果表明其位移已经达到蓝色预警值。

综合理论计算结果和实际位移监测以及气象局所获得降雨量结果来看，本文计算的降雨阈值具有一定的合理性。但是，工程地质条件的复杂多样性和土体强度随着含水率变化的非线性等使得滑坡预警工作在监测过程中存在很多困难。因此，可加强宏观迹象特征的群测群防来提高监测预警工作保障。

4 结语

(1)纳缝坡体裂缝发育为降雨提供便利的渗流通道，坡体内部土体粘粒含量约21%，粉粒约32%，砂粒组含量47%，级配良好，渗透特性较低。分析认为：坡体表面开裂，易于入渗水体不断在坡体内部聚集，土体有效应力降低，并伴随滑带处不断增大的静水压力是滑坡变形的内在原因。

(2)对斜坡的PFC3D数值计算，证明了静水压力的大小与坡体的稳定性关系紧密：随着静水压力增大坡体稳定性逐渐降低，滑坡破坏临界静水压力所对应的临界静水位高度在1.70-1.75m之间。

(3)通过降雨入渗深度反演静水位高度所对应的降雨强度以及降雨时长，最终确定缝滑坡预警参数阈值：降雨强度I=5.73×10-5cm/s，降雨时长T=8.87d，累积降雨总量Q=439.13 mm，并进一步通过实测数据验证了预警参数阈值的合理性。