基于TRIGRS与Scoops3D耦合模型的潜在滑坡稳定性时空动态预测

郑玲静,李秀珍,姚 杰,余文秀,2

(1. 中国科学院、水利部成都山地灾害与环境研究所,四川 成都 610041; 2. 中国科学院大学,北京 100049;3. 成都市城市安全与应急管理研究院,四川 成都 610031)

0 引言

斜坡失稳破坏与内在自身力学特性及外界触发因素密切相关。在降雨、融雪等外界因素诱发下,陡峻山区极易发生斜坡失稳破坏。据初步统计,在全球范围内每年因滑坡灾害造成的经济损失高达百亿美元,更有数以千计的人因滑坡灾害而失去生命[1]。气候变化大背景下,极端强降雨现象加剧了滑坡灾害的发生,使全世界人民生命、财产安全面临更大的滑坡灾害威胁,如2020年8月,强降雨导致巴基斯坦南部山区发生数起山体滑坡,滑坡、洪涝等灾害造成至少90人死亡。据大量资料统计,由降雨诱发的滑坡约占滑坡总数的80%左右[2],降雨是滑坡的主要诱因。

区域降雨滑坡的稳定性评价与预测是目前滑坡研究中的热点问题之一[3]。其研究方法主要包括统计分析法和确定性模型法两大类[4],其中确定性模型法为当前的主要研究方向。此类方法以降雨诱发斜坡失稳的力学原理为基础,结合极限平衡理论建立力学模型,对特定降雨条件下滑坡的稳定性进行评价与预测。物理确定性模型从一维与二维模型逐步向三维模型发展,一维与二维模型以SHALSTAB模型[5]、dSLAM模型[6]、SINMAP模型[7]、Iverson模型[8]、TRIGRS模型[9]为代表,但以上模型忽略了滑坡的滑动方向及三维几何特性,存在过度简化滑坡发生机理与难以反映滑坡力学机制等问题。而已有的广泛使用的FLAC3d、ANSYS和Abacus等商用成熟软件仅适用于三维单体滑坡的稳定性分析及预测。Scoops3D模型是美国地质调查局开发的三维斜坡稳定性分析模型,具有建模简单,可接入多种降雨入渗与水文模型且可以同时考虑复杂孔隙水压力与地形三维分布情况,适用于降雨诱发的区域潜在滑坡的识别及评价研究。该模型不仅在搜索滑面时考虑了栅格间的相互作用,还具有可设置空间纵向分层岩土体参数等优势,目前已成功应用于潜在滑坡稳定性识别与预测中。例如Scoops3D开发者BRIEN与REID[10-11]耦合MODFLOW-2000与Scoops3D模型,建立了考虑三维孔隙水压力的深层潜在滑坡三维危险性分析模型;TRAN等[12]首次耦合了TRIGRS模型和Scoops3D模型,对降雨条件下潜在滑坡的分布范围进行了有效评价;姚杰等[13]也利用TRIGRS模型与Scoops3D模型进行耦合,对拟建川藏铁路沿线典型段潜在滑坡三维稳定性进行了动态识别研究。这些模型的耦合不仅能充分发挥各个单一模型的优势,同时也提高了潜在滑坡识别的研究程度和应用效果。

根据实地考察与遥感解译,查明中巴公路KKH沿线加格洛特(Jaglot)至哈维连(Havelian)段主要为降雨型滑坡区段,选择10 km范围区域为研究对象,模拟计算百年一遇极端降雨情景下雨强,利用 TRIGRS 降雨入渗模型和考虑了岩土体空间纵向分层的 Scoops3D 三维模型耦合,建立能够较为客观反映滑坡孕育条件的潜在滑坡三维稳定性评价模型。研究天然状态下与百年一遇极端降雨情景下潜在滑坡的空间分布及动态变化情况,并对比分析和讨论了一维TRIGRS模型与三维耦合模型的预测结果,以期能够更加客观地对研究区的潜在滑坡进行识别及评价。

1 研究区概况

中巴公路(Karakoram Highway,KKH)地处三大世界最年轻山系的交汇区,是全球气候变化最为敏感和复杂的地区,沿线滑坡灾害十分发育。文中选取加格洛特-哈维连段KKH沿线两侧10 km范围区域作为主要研究区。

1.1 自然环境条件

研究区深大断裂带下切作用强烈,属深切割高山河谷地貌,具有地势险峻、群山耸立、河谷深切、峡谷较多等特点。区内新生代印度板块不断北移与亚欧板块发生强烈撞击,造成该地区整体产生强烈的垂直差异运动,地质构造活跃,复杂的深大断裂带极为发育,主要呈北西走向分布(图1)。区内及其周边发育的活动断裂带主要有主地幔逆冲断裂带(MMT)、主边界逆冲断裂带(MBT)与主喀喇昆仑逆冲断裂带(MKT)三大断裂带。断裂带周边岩石大多会表现出不同程度的破碎与变质,岩性主要以砂岩、砾石、沉积岩为主,区内地质条件复杂。

图1 研究区地层及断裂带分布图Fig. 1 Distribution map of strata and fault zone

研究区属于南亚次大陆亚热带气候区,地形高差大,且受高山降雨阴影效应影响,时空区域分异与气候垂直分带显著,雨季降雨明显易诱发滑坡灾害。

1.2 滑坡灾害概况

滑坡是中巴公路上发生最为频繁的山地灾害之一,一方面研究区高寒冰冻、地表风化严重,岩层变质与破碎现象明显,受剥蚀后形成大量松散物堆积,为滑坡的发生提供内在条件;另一方面区域构造抬升,河床下切侵蚀,山高谷深,雨季降雨集中,为滑坡的启动提供外在动力。受季风影响,雨季降雨量增加明显,导致山体滑坡发生更为频繁。

本研究主要选择KKH沿线降雨型滑坡集中发育的加格洛特到哈维连段为主要研究区段,该段坡体主要为第四系松散堆积物,少部分为三叠系、二叠系和石炭系沉积岩和志留系片岩、板岩等变质岩,大部分属于易滑岩组(图1)。该区段发育的滑坡主要以降雨诱发的堆积层滑坡为主,破坏类型多为圆弧形破坏或平面破坏。经课题组野外调查与遥感解译,滑坡分布如图2(a)所示,图2(b)与(c)分别为研究区加格洛特与吉拉斯(Chilas)附近的典型历史滑坡。

图2 研究区历史滑坡分布图Fig. 2 Historical landslides in the study area

2 方法

2.1 TRIGRS模型

TRIGRS模型是基于无限边坡理论的降雨诱发型边坡危险性评价模型,模型主要有3个模块,分别为降雨入渗模块、水文模块和稳定性分析模块,充分考虑了降雨入渗土体而导致的瞬态孔隙水压力改变情况,其优势在于可分析不同时刻边坡的孔隙水压力、体积含水量和边坡稳定性等动态变化情况。

降雨入渗模块基于IVERSON[14]提出的Richards的线性解,可计算特定时间内土层指定深度的渗流情况,能较为真实地反映了降雨入渗条件下的边坡渗流变化,本研究中假定深层为基岩风化层(与Scoops3D模型对应),渗透系数极小,因此压力水头的计算使用下边界有限深条件,表达式见式(1):

(1)

式中:Ψ(Z,t)为地下水压力水头,与土厚和时间相关;Z为土层垂直方向的厚度,垂直向下为正;t为降雨入渗时间;d为岩土体初始地下水位深度;N为时间序列;InZ为在第n个时间段入渗速率;Ks为垂向饱和渗透系数;m为收敛级数;H(t-tn)为Heavyside阶梯函数,与第n个时刻对应雨强相关,是它的时间函数;tn为第n个雨强出现的时间;dLZ为基底的深度;ierfc(η)为高斯补差函数对变量η的一次积分值,表达式在无穷级数时会迅速收敛。其中,β和D1见式(2)和式(3):

(2)

(3)

式中:δ为坡度;IZLT为初始表面通量;D0为饱和水力扩散系数。

水文模块假定计算栅格在每一个运算时段内水量守恒,即在当前时段内不能下渗的降雨将作为地表径流流向下游栅格,不会无故消失。

稳定性分析模块基于极限平衡理论结合无限边坡模型计算研究区的安全系数,表达式如式(4):

(4)

式中:Fs(Z,s)为安全系数;φ为土层内摩擦角;c为土层凝聚力;γsat为土容重;γw为地下水容重;δ为坡度;ψ(Z,s)为压力水头,是深度和时间的函数。一般而言,内聚力c与φ的值越小、地下水位越高,安全系数越小,坡体越不稳定。

2.2 Scoops3D模型

Scoops3D模型可以系统全面地对研究区的潜在滑坡进行三维搜索及识别,充分考虑复杂地形与地下水等条件对于潜在滑坡稳定性的影响,具有计算范围广、识别效率高等优势。

Scoops3D 模型基于DEM栅格单元将其网格处理为三维柱体,并在模型上方整体空间范围间隔一定距离生成若干个搜索球心,并以球心为原点参考递增的方式生成半径,以此生成不同球型曲面对斜坡体进行切割,切割得到的岩土体利用稳定性计算公式计算安全系数,球型曲面和三维柱体相交的曲面标记为潜在滑面,被球形曲面切割所包含的三维柱体即为潜在滑坡,如图3。滑坡的搜索将设置体积或者面积阈值,当 Scoops3D 模型在搜索过程中达到最大阈值时,将停止增加搜索半径长度,并采用极限平衡理论对所有与球型曲面相交的三维柱体进行计算,判断其稳定性情况,最终记录最小安全系数所对应的就是稳定性最差的潜在滑动面。本研究采用的是 Bishop 理论进行斜坡体危险性计算,具体见式(5):

图3 Scoops3D潜在滑面示意图Fig. 3 Schematic diagram of potential sliding surface in the Scoops3D model

(5)

式中:Ri,j为Scoops3D 进行搜索时的半径;i、j为分别为三维柱体的第i行和第j列;ei,j为搜索球心到三维柱体质心的距离;ci,j为潜在滑坡体的黏聚力;Ai,j为潜在滑坡面的面积;φi,j为潜在滑坡体的内摩擦角;ki,j为水平振动荷载;Wi,j为潜在滑坡体的重力;ui,j为表示三维住体内部的孔隙水压力;αi,j为表示潜在滑坡体的视倾角;βi,j为潜在滑坡体的真倾角。

2.3 TRIGRS和Scoops3D模型的耦合模型

Scoops3D模型无法输入具体降雨数据,TRIGRS模型在降雨入渗水文计算方面理论较为成熟,且应用较广,本研究将2个模型耦合,示意图如图4。

图4 TRIGRS与Scoops3D模型耦合原理图Fig. 4 Schematic diagram of TRIGRS model and Scoops3D model coupling

结合实际野外调查情况,确定研究区上层是较为松散的土壤层,下层是略为坚硬的基岩风化层,文中尝试对研究区的岩土体进行纵向空间分层,以期能够更加真实地对研究区的潜在滑坡进行搜索及识别。耦合过程中,首先利用TRIGRS一维斜坡稳定性分析模型中的入渗模块模拟的降雨渗流,计算并得到在降雨不同历时时刻的二维孔隙水压力图层文件,然后将孔隙水压力与研究区DEM利用Scoops3D模型建模,结合研究区岩土体情况,线性插值获取不同降雨历时下的三维孔隙水压力应力场,对研究区的潜在滑动面进行三维搜索与识别,获取在天然状态和百年一遇降雨情景下中巴经济走廊KKH沿线降雨型滑坡为主的区段潜在滑坡分布及变化情况。

3 应用

3.1 模型参数选取

3.1.1 降雨参数

经调查,收集Kakul雨量站1952—2012年间最大日降雨量数据,该雨量站位于研究区内哈维连与曼塞赫拉(Mansehra)之间,经度为34.2036°E,纬度为73.2820°N。将最大日降雨量数据利用“皮尔逊Ⅲ型频率曲线”计算经验频率[15],结果如图5所示,可得到当重现期为100 a时,最大日降雨量为231.1 mm/d,在本研究中百年一遇极端降雨情景下降雨强度取2.67477E-06 m/s。

图5 P-Ⅲ型频率曲线Fig. 5 P-Ⅲ frequency curve

短时强降雨条件下,浅层土体地下水位会随降雨入渗发生急剧变化,达到饱和。结合实际野外调查情况,确定研究区上层是较为松散的土壤层,下层是略为坚硬的基岩风化层,多为浅层滑坡,TRIGRS模型模拟瞬态降雨入渗的方式采用饱和入渗,并将在各时段的水压分布文件作为输出,与Scoops3D模型进行耦合,其Scoops3D水文模型中采用变化的饱和水压三维空间分布的水文模型。

3.1.2 土层厚度

土层厚度不仅受气候、生物、母质、地形和化学与物理过程等因素影响,此外,还会因为岩性、坡度、曲率和植被覆盖等因素的改变而变化[16]。参考近年来已提出多种估计土厚空间分布的模型,如均一模型、分级模型和线性模型等[17-19]。结合相关研究成果比较了多种常用模型预测的滑坡结果,确定在TRIGRS与Scoops3D模型采用线性土厚模型[20]。

该模型假设土厚与坡度呈线性分布关系,经野外调查与参考其它文献[17-20],最大坡度对应最小土厚(0.1 m),最小坡度对应最大土厚(3.5 m),研究区坡度在0～70.51°之间,因此确定上层土层厚度(y)与坡度(x)之间的函数关系为:y=-0.048 22x+3.5,下层为无限延伸的岩体风化层。借助ArcGIS平台,得到最终上层土层厚度分布如图6所示。

图6 土层厚度空间分布图Fig. 6 Spatial distribution of soil thickness

3.1.3 岩土体物理力学参数

针对已有区域滑坡三维研究中未考虑纵向岩土体岩性存在差异问题,结合本研究区岩土体特征及分布情况,初步尝试将岩土体进行纵向空间分层,将研究区段岩土体纵向分为土壤层和基岩风化层2层。

结合参数反演,采取文献调研以及工程地质类比等手段,分析对比不同参数在研究区的适用效果,最终确定中巴经济走廊KKH降雨影响段的岩土体参数取值[21-23],具体见表1。水文参数参考裴钻对中巴公路的研究[24,综合渗透系数取值为:KS=2.87E-6,根据经验关系可以得到其它两者的取值分别如下[25]:水力扩散系数D0=200Ks=5.74E-4,入渗速率IZ=0.01Ks=2.87E-8。

表1 岩土体参数取值表Table 1 Values of geotechnical parameters

3.1.4 Scoops3D模型搜索参数

研究区面积约为7 259 km2,考虑模型性能与计算效率,DEM分辨率采用90 m×90 m(数据下载自地理空间云)。研究区海拔在465～4 975 m,DEM行列数分别为1 982与2 186,结合研究区现状,Scoops3D模型的搜索球心高程为470～5 000 m,搜索起点坐标为行编号与列编号均为1,相邻球心间隔的栅格数为4。每次搜索通过在DEM上方生成搜索球心矩阵(搜索点阵),使用由粗到细的搜索方式(Coarse-to-fine), 搜索半径以70 m递增,迭代搜索每一个节点分析三维潜在滑动面。模型搜索时会出现在每个节点搜索过程中有无限多潜在滑动面的情况,因此,在文中使用体积限制潜在滑体来确定最终滑动面,其中根据野外调查与遥感解译研究区滑坡灾害方量,确定最终搜索体积阈值范围为105～2×108m3。

3.2 潜在滑坡预测结果

本研究首先基于岩土体纵向空间分层模型,利用Scoops3D模型模拟天然状态下(无降雨)研究区段滑坡稳定性空间分布情况;再基于TRIGRS与Scoops3D耦合模型,预测了百年一遇极端降雨情景下不同降雨历时潜在滑坡动态变化。根据GB/T 32864—2016《滑坡防治工程勘查规范》划分潜在滑坡,见表2,结果见图7与表3。

表2 潜在滑坡分区标准表Table 2 Potential landslide classification standard

图7(a)与表3结果显示,在天然状态下研究区域加格洛特-吉拉斯段与瑟津(Sazin)-达比尔(Dubair)段零星分布着潜在滑坡,其中潜在滑坡区域分布范围占研究区段总面积的4.14%,天然状态下斜坡稳定性较好,历史滑坡灾害数量为4个,占研究区总灾害数量的7.69%,潜在滑坡区主要沿高差悬殊、地势险峻斜坡分布。

表3 潜在滑坡稳定性结果统计表Table 3 Statistical table of stability calculation results of potential landslides

图7(b)～(d)与表3结果显示,在百年一遇极端强降雨情景下,潜在滑坡区相较天然状态下明显增大,潜在滑坡区主要集中分布在地形高差较大的地区,如加格洛特-吉拉斯段中巴公路两侧与瑟津-塔科特(Thakot)段;吉拉斯-瑟津与塔科特-哈维连段相对较为安全,只分布有零星的潜在滑坡区。降雨对研究区潜在滑坡影响较大,且随降雨时间的延续,潜在滑坡区的分布范围逐渐增大,滑坡灾害点数量占比逐步增加,当降雨1 h后,潜在滑坡面积占研究区段总面积的18.44%,约为天然状态下潜在滑坡区总面积(4.13%)的4.5倍,研究区潜在滑坡对降雨响应明显,区内灾害点数为27个,占比为51.92%;降雨4 h后,潜在滑坡分布面积占研究区段总面积的20.57%,区内灾害点数为29个,占比为55.77%;降雨12 h后,潜在滑坡分布面积占研究区段总面积的37.54%,区内灾害点数为37个,占比为71.15%,12 h后,滑坡分区几乎无明显变化。根据降雨入渗理论可知,在恒定降雨条件下,降雨入渗率随降雨时间的延续,会逐渐减小并趋于稳定,渗入坡体的有效降雨量也会随着土体饱和度的增加逐渐减少。当土体达到饱和状态时,降雨入渗量趋近于0,降雨更多地转化为了地表径流,斜坡稳定性也随之愈来愈无明显变化。

图7 研究区潜在滑坡稳定性模拟结果Fig. 7 Stability simulation results of potential landslides in the study area

3.3 结果可靠性评价

在已有研究中往往采用历史滑坡点数量占各分级等级总面积的百分比来评价模型性能[13],但该评价指标总体较为笼统、简单。文中引入点状滑坡性能百分比指数%LRclass对滑坡稳定性预测结果性能进行评价[12],其%LRclass指数计算公式如下:

(6)

(7)

式中:S为在每类Fs中滑坡点数量所占总滑坡数的百分比;A为对应每类Fs的预测面积占总面积的百分比。%LRclass指数评估基于Fs<1与Fs≥1分级,但同时考虑了模型预测结果中在两类分级区中发生滑坡的情况,当安全系数小于1时,指数值越大,模型的预测结果越好。

由表3可以看出,在天然状态下,%LRclass指数在不稳定区(Fs<1)为67.60%,在百年一遇极端降雨情景下,在降雨1、4、12 h后,%LRclass指数在不稳定区均在79%以上,分别为79.49%、79.71%和79.82%,耦合模型计算结果与历史滑坡灾害分布一致性较好,可信度较高。

4 讨论

本研究为了进一步评价TRIGRS和Scoops3D三维耦合模型对降雨滑坡稳定性预测的性能,将三维耦合模型与一维TRIGRS模型结果进行了分析比较。在百年强降雨情景下,2种模型预测的潜在滑坡稳定性空间分布、各稳定性分区的历史滑坡数量占比、面积占比与%LRclass指数结果如图8与表4所示。

图8(a)与(d)对比可知,TRIGRS模型预测的潜在滑坡区空间分布更广泛,其中加格洛特-塔科特段潜在滑坡区分布均较为集中,根据图示历史滑坡分布可以看出吉拉斯-瑟津段中部几乎无滑坡发生,由此可见,耦合模型更符合实际;图8(b)与(e)显示,耦合模型与TRIGRS模型不稳定区的滑坡数量占比分别为65.38%和46.15%,耦合模型预测结果明显优于TRIGRS模型;图8(c)与(f)显示,耦合模型与TRIGRS模型潜在滑坡面积占比分别为37.54%和47.27%,TRIGRS模型对潜在滑坡的预测结果更偏不安全。这与已有的相关研究结果也是一致的[3]。

从图7和图8中可以看出,耦合模型预测的潜在滑坡区呈片状集中分布,TRIGRS模型预测的潜在滑坡区呈细条状零散分布。究其原因应该主要在于,耦合模型是利用球型曲面将坡面切割成若干三维柱体,考虑了柱体间的相互作用,故模拟潜在滑坡区整体性较强,呈片状集中分布;而TRIGRS模型中每个栅格都进行单独计算,因此每个栅格的安全系数均为独立值[10],因此模拟潜在滑坡区整体性较差,总体呈细条状零散分布。

图8 研究区潜在滑坡稳定性一维与三维模拟结果Fig. 8 Results of 1D and 3D simulation of potential landslide stability in the study area

另从表4也可以看出,TRIGRS模型预测的不稳定区的滑坡数量占比(46.15%)明显低于耦合模型(65.38%),TRIGRS模型预测的不稳定区面积(34.26%)略高于耦合模型(32.33%),一维TRIGRS模型与三维耦合模型的%LRclass(Fs<1)指数分别为62.20%和79.81%,耦合模型明显优于TRIGRS模型。2种模型在降雨滑坡预测中均有一定适用性,但一维TRIGRS模型的过度简化致使难以考虑实际边坡结构、载荷和滑动面的三维地形分布,基于岩土体纵向空间分层模型,使用三维耦合模型利用体积搜索滑动面更适合处理复杂地形的滑坡,能获得更合理的预测结果。

表4 研究区一维与三维计算结果统计表Table 4 Statistical table of 1D and 3D calculation results of potential landslide in the study area

5 结论

1)天然状态下加格洛特-吉拉斯段与瑟津-达比尔段零星分布着潜在滑坡,其中潜在滑坡分布面积占研究区段总面积的4.14%,天然状态下斜坡稳定性较好。

2)研究区潜在滑坡对降雨响应明显。百年一遇极端降雨情景下,降雨1 h耦合模型预测的潜在滑坡面积占比(18.44%)约为天然状态下(4.13%)的4.5倍;随时间延续,潜在滑坡区范围逐渐增大,降雨12 h,潜在滑坡区分布面积占研究区段总面积的37.54%,主要沿高差悬殊、地势险峻斜坡分布。

3)TRIGRS模型在潜在滑坡沿加格洛特-塔科特广泛分布,潜在滑坡面积占比达48.27%,模拟中存在过度预测问题;耦合模型预测的潜在滑坡区主要分布在加格洛特-吉拉斯段与瑟津-达比尔段,区内71.15%历史滑坡灾害点被正确预测,更符合历史滑坡灾害分布特征。

4)三维耦合模型在不同时间段%LRclass(Fs<1)指数均在79%以上,预测潜在滑坡效果明显优于一维TRIGRS模型(62.20%)。

TRIGRS模型与考虑纵向岩土体空间分层的Scoops3D模型耦合在该研究区取得了较为理想的结果,相关研究成果可为中巴经济走廊建设提供可靠防灾减灾应急管理科技支撑。在下一步研究中,可进一步划分工程地质区段,结合地质特征确定纵向岩土体分层情况,叠加研究区常见雨型,同时对研究区岩土体及水文等参数的空间不确定性进行进一步深入研究。