黄河上游典型流域滑坡稳定性预测及模型应用

2022-06-09孙何生邱海军朱亚茹高祥语

西北大学学报（自然科学版） 2022年3期

孙何生,邱海军,2，3,朱亚茹,刘雅,高祥语

(1.西北大学城市与环境学院/地表系统与灾害研究院,陕西西安 710127;2.西北大学陕西省地表系统与环境承载力重点实验室,陕西西安 710127;3.西北大学陕西省黄河研究院,陕西西安 710127)

黄河中上游地区地质构造复杂,降雨集中,山地灾害频发,每年都会造成大量的人员伤亡和财产损失[1-7]。近年来,随着全球气候变化和人类活动加剧,滑坡灾害出现加重趋势,严重影响了该区域的生产生活和生态安全[8-9]。因此,进行滑坡稳定性分析尤为重要,这不仅可以为滑坡灾害风险防控提供理论与基础数据支撑,还有利于黄河流域的生态保护与高质量发展[10]。

滑坡稳定性分析对于滑坡风险防控具有重要意义。在进行滑坡稳定性分析时,目前主要分为定性分析和定量分析两种。其中,定性分析是通过利用区域现有的地质、水文、地形、土地利用方式等滑坡孕灾因子信息,结合专家经验判断,对研究区的滑坡稳定性进行评定分级,最后确定滑坡的稳定性。这种分析方法最大的缺点就是依赖于专家的经验,主观性较强,不能做到统一标准[11]。而定量分析就很好地解决了这一问题,它是通过一种数学方法来估算滑坡的稳定系数[12-13]。物理确定性模型作为定量分析的一种类型,受到了众多学者的青睐,特别是以SINMAP模型[14-15]、TRIGRS模型[16-19]和LAPSUS模型[20-21]代表,将降水因素、景观因素考虑在内,得到了较好的模拟效果。但是，作为一种二维模型,没有充分考虑地层岩性、地下水以及孔隙水等因素对模型的影响。本文所采用的Scoops3D模型作为一种三维模型,不仅考虑了众多的因素,例如地震等,而且克服了以往许多模型的局限,对于小区域的滑坡稳定性预测有着良好的效果[22]。因此,近年来三维模型越来越受到了学者们的应用,以此来解决实际问题[23-25]。

对于物理确定性模型来说,参数的选取会直接影响模拟的效果。数字高程模型作为模型必要输入条件,其分辨率的大小直接影响模型的预测性能,许多学者都对此做过研究,然而却出现了两种不同的声音。一种是分辨率越高,对于滑坡敏感性制图和模拟效果更高[26-30]；另一种是分辨率并非越高越好,使用精度过高的DEM进行滑坡稳定性分析时,会存在预测准确性不足或者过度预测的问题[31-33]。辛星等人在对黄土浅层滑坡稳定性分析时,分别用3 m、5 m和10 m三种不同的分辨率来试验,最后发现3m的结果更优[12],但对于三维模型Scoops3D来说,DEM分辨率对于模型的影响研究者仍较少。Sarma等人在对TRIGRS模型的研究时,分别输入不同分辨率的DEM和不同强度的降水,最后发现降水对于模型的影响更强,而分辨率不是决定模型预测结果的决定性参数[34]。对于三维模型来说,DEM分辨率是影响模型结果的重要参数,但其他参数对于模型的结果影响同样显著。对于岩土参数相差不大的同一小区域/流域而言,由于对滑坡大小的定义不同,在模型中输入不同的滑坡搜索范围,也会极大地影响模拟结果。因此,本文主要探究DEM分辨率和滑坡搜索范围对于模型结果的影响,通过与实际的滑坡位置进行对比验证,获得最佳的参数组合。通过探究不同的参数对于模型在滑坡稳定性预测方面的影响,可以在利用Scoops3D模型进行防灾减灾过程中提供更精确、科学的建议,大大提高模拟精度。

1 研究区概况

研究区甲加沟位于黄河上游的支流,地处中国西北部的青海省(见图1)。从地貌上来看,属于青藏高原向黄土高原的过渡地带,处于一、二级阶梯的中间地带[35]。在褶皱构造上隶属于祁连山褶皱段,是一种中新生代山间断陷谷地,造成了地震频发,因而成为滑坡的典型发育区[36-37]。从研究区的裸露地层岩性来看,土壤类型主要包括第三纪红黏土和第四纪黄土,土壤层的厚度不一[38]。该研究区的纵长为6 km,横宽为3 km,面积约为18 km2。最低海拔在1 900 m左右,最高海拔达到了2 800 m,平均海拔也有2 500 m。由于深居西北内陆,海洋暖湿气流气流难以到达,再加上地势的抬升,形成了高原大陆性半干旱气候。从黄河流域的气候湿润状态来看,属于黄河干流黄土丘陵干旱区,降水量较少,然而由于山高谷深得天独厚的地形优势,再加之蒸发量少的缘故,形成了以高山草甸和草甸草原为主的植被景观。在研究区内不仅存在着起伏和缓的丘陵,还有面积广阔的高原台地[39]。

图1 研究区位置及滑坡分布图Fig.1 Location and landslide distribution map of the study area

2 Scoops3D模型

2.1 模型原理

Scoops3D是由美国国家地质调查局于2015年研发,基于数字高程模型,利用极限平衡法,通过简化的毕肖普法来对每个三维滑坡面进行计算,获取各滑动面上的栅格安全系数。Scoops3D是以DEM为基础进行的斜坡三维稳定性分析。首先模型以每一个椭球形的潜在滑动面的旋转中心为节点,搜索并分析每一个单元格所对应的安全系数,进而得到不稳定块体的位置和体积[24]。

Scoops3D通过线性的Coulomb-Terzaghi失稳准则来计算每个滑动面的抗剪强度,如公式1所示。

s=c+(σn-u)tanφ

(1)

式中：s为抗剪强度；c为土体的黏聚力；φ为土体的内摩擦角；σn为滑体所受的正应力；u为作用于剪切面的孔隙水压力。定义安全系数F为抗剪强度s和剪切应力τ的比值,如式(2)所示。

(2)

从理论上说,当F<1时,表示斜坡不稳定,发生滑坡。在平衡状态下,剪切应力τ等于抗剪强度s乘比例常数1/F。

本文采用简化的Bishop法来计算安全系数,因此可以用式(3)表示。

(3)

式中:Ai,j为潜在球形滑动面在三维柱体内的表面积;αi,j为滑面的视倾角;βi,j为三维柱体潜在滑面的倾角;Wi,j为三维柱体的重量;ui,j为作用在三维柱体滑面上的孔隙水压力;ki,j为作用在三维柱体中心的水平震动荷载;ci,j和φi,j为三维柱体滑面的抗剪参数[24]。

Scoops3D将水文条件分为有水和无水两种状态,有水的情况下又分成了几种状况。由于本文所选研究区位于黄河上游,属于半干旱气候条件,受地下水影响较小,可以忽略不计[13],因此安全系数F可以表示为

(4)

2.2 模型敏感性分析

物理确定性模型在分析滑坡稳定性时,所输入的参数会直接影响模型的预测性能,因此需要在输入模拟参数前掌握详细又准确的数据资料。然而在实际工作中,往往受制于研究区实际的地形、时间、滑坡发生前的真实情况等条件的限制,不能获取每个参数的详细信息[14]。因此,通过提前对参数对模型的结果影响进行分析,可以确定参数对模拟效果的影响程度来获取关键参数,这样就可以降低成本和时间。

在模型敏感性分析中,首先选择一组容易计算的值,将这组值作为标准值,然后逐一调整变化每一个参数,经过计算安全系数的变化率来反应该参数对模型的敏感性[8，35]。表1反映了Scoops3D模型敏感性分析选用的各个参数标准值及变化范围。

表1 参数的标准值及变化范围Tab.1 Standard parameters and their variation ranges

对于参数敏感性分析的结果有两种情况,一种是参数变化与安全系数变化成正相关,一种是负相关关系。安全系数随着黏聚力、内摩擦角、球体半径、滑动面面积和潜在滑动角的增大而增大,随着滑动方向视倾角、单元重力、水平加速度系数和球体半径垂直分量的增大而减小[11，13]。同时从变化幅度来说,模型计算的安全系数对内摩擦角、球体半径和潜在滑坡面的敏感性相对较低,而对黏聚力、滑动方向和栅格单元重量更敏感[14]。

3 数据获取与分析

3.1 DEM数据

本文选取的数字高程模型(DEM)数据共有3种不同的分辨率,分别是5 m、12.5 m和30 m。其中5 m分辨率是购买的日本ALOS卫星的AW3D 5m DEM产品(https:∥www.aw3d.jp/en/),ALOS PALSAR 12.5m分辨率是从AlaskaSatellite Facility网站免费获取的12.5m的DEM产品(https:∥search.asf.alaska.edu/#/);SRTM 30m是从美国USGS官网免费下载的DEM产品(http:∥gdex.cr.usgs.gov/gdex/)。

3.2 现场调查及建立滑坡编目

首先，基于遥感影像和谷歌影像对研究区的滑坡进行了解译,并对解译的滑坡进行了野外验证。野外调查一方面是对解译的滑坡进行现场确认并收集相关数据,另一方面是对一些现场新发现的滑坡进行定位并收集数据信息。滑坡编目作为一个数据库,记载着滑坡的相关属性,对于研究滑坡的成灾机理和防灾减灾等方面有着重要的价值。此外,滑坡编目的建立也是模型运行前的准备工作,在模型预测结束后作为一个现实标准与模拟结果进行比对,以此来判断模型的性能。因此,滑坡编目的建立起着至关重要的作用。

本研究区内共发现滑坡66处,有平移式滑坡、旋转型滑坡和崩塌等3种类型。由于难以测算滑坡的体积,因此选择了面积这一可行性标准来判定滑坡的大小。这些滑坡的大小各异,最小的滑坡面积仅为778m2,而面积最大的超过了1×105m2。从量级上看,滑坡的大小相差比较大,跨越了多个量级。在运行模型时,针对滑坡的面积大小将是我们下文的重点研究内容。为了对模型就行精度的验证,同时在非滑坡区域又选取了60个非滑坡点进行验证。

3.3 岩土参数

Scoops3D模型在进行滑坡稳定性模拟分析时,所需的输入参数主要包括建立搜索矩阵所需的参数(最小滑坡面积和最大滑坡面积,搜索半径,搜索高度)和模拟计算所需的岩土参数(内摩擦角,黏聚力,土壤容重)。搜索矩阵所需的参数可以由滑坡编目获得,而对于模型中需要的岩土参数,主要通过室内实验获得。首先利用ArcGIS10.5软件将研究区均匀划分为20个网格,然后在每个网格随机取一个样点,这些工作方便了在实际野外采取土样中位置的选择与确定。在实际采样过程中,依据地形地势的原因适当进行调整,在选取的位置附近进行取样(见图2)。在采样过程中,结合浅层滑坡的厚度以及通过对现场滑坡的厚壁高度进行估算,最后选择在1.5～2 m的深度进行采样,每个样点采取4个土样,将现场取得的样品塑封带回,利用直剪仪获取黏聚力和内摩擦角等参数。对于模型输入参数,在去掉极端值后再求取每个参数的平均值,最后本文所选取的岩土参数是内摩擦角9.4°,黏聚力14.3 kPa,土壤容重15 kN/m3。

A 取样工具;B 取样坑;C 环刀取样;D 实验仪器图2 野外岩土取样和室内实验图Fig.2 Sampling and experiments

3.4 滑坡搜索面积的确立

对于Scoops3D而言,输入不同的滑坡搜索面积会大大影响模拟的效果,因此本文结合滑坡编目,发现本区域内的滑坡大小差异较大,过大或过小的搜索范围都会影响最后的结果,因此选择合适的搜索范围对模型至关重要。结合前人的研究和现有的资料[38],根据滑坡的面积来对滑坡进行分类,一类为小型滑坡,一类为中大型滑坡。其中小型滑坡面积在1×104m2以下,从研究区的实际状况来看,这部分滑坡数占总数的39%,而在实际选择时,为了节约模型运行时间,对于低于1×103m2不予讨论,并且从滑坡编目中了解到这一范围仅有一个,可以忽略不计。而对于中大型滑坡来说,它的面积在1×104m2以上,结合滑坡编目与模型运行来看,排除了那些面积特大而数量又少的,这可能是在勾画滑坡范围时将不确定的边界考虑在内导致的,因此选择1×104～4×104m2作为中型滑坡搜索范围,这一区间的滑坡数目占全部的44%(见图3)。因此,综合考虑,选择了两种不同的搜索面积来对模型进行检验,分别是1×103～1×104m2和1×104～4×104m2。

图3 滑坡面积分布等级图Fig.3 Grade diagram of landslide area distribution

4 结果分析与模型验证

4.1 预测结果及分析

为了探究Scoops3D的不同参数下的模拟性能,本文一共采取了3种分辨率(5m,12.5m和30m)和2种不同的搜索面积(1×103～1×104m2,1×104～4×104m2)分别组合共3组6种工况来组合进行分析。第1组为工况1：5m分辨率和1×103～1×104m2、工况2：5m分辨率和1×104～4×104m2,第2组为工况3：12.5m分辨率和1×103～1×104m2,工况4：12.5m分辨率和1×104～4×104m2,第3组为工况5：30m分辨率和1×103～1×104m2,工况6：30m分辨率和1×104～4×104m2。每一组工况下采用同一分辨率、不同的搜索面积,这样方便在后面的比较。根据三维模型的计算结果,将安全系数值分为5级:1.5≤F,表示稳定;1.2≤F<1.5表示潜在稳定;1≤F<1.2,表示潜在不稳定;0.75≤F<1,表示不稳定;F<0.75,表示极其不稳定[12，40]。

从图4可以看出,在6种工况下,工况1、3、5极不稳定区域和不稳定区域分别占研究区总面积的43.10%、33.40%和19.15%,工况2、4、6极不稳定区域和不稳定区域分别占研究区总面积的52.62%、48.76%和41.06%。从分辨率的角度看,分辨率越高,模型计算的安全系数值越低,说明滑坡稳定性越低。特别是在搜索面积较小时,这种由分辨率带来的差别比在搜索面积大时更明显,这在一定程度上说明了分辨率越高的DEM在识别小型滑坡上更加有优势。而在同一组工况中,即分辨率相同分辨率下,当搜索的面积变大时,它所识别的不稳定区域也就越多,分别增加了9.52%、15.36%和21.91%,呈现一种递增趋势,这就说明了DEM分辨率和搜索滑坡面积都会影响最后的模拟效果。在这6种工况中,工况5的比例最低,这可能是由于分辨率比较粗糙,由于搜索的小型滑坡的面积较小,一个滑坡所占栅格数量较少,最后在Scoos3D计算过程中整体的安全系数较高。从图4中可以看出,极不稳定区域和不稳定区域主要分布在沟谷的两侧,由于沟谷两侧的地形比较高陡,滑坡数目较多,而稳定区域分布在研究区的西部地带。在实地调查中，这一区域为梯田,相对平缓的地带滑坡数目较少，这从侧面反应了地形条件是影响滑坡稳定性的重要因素。从而可以得出Scoops3D模型在小区域的滑坡稳定性方面具有良好的预测效果。

通过与滑坡编目的空间位置对比,在同一小型搜索面积不同的分辨率条件下,分别有61、47和24个滑坡点位于极不稳定区和不稳定区,占全部滑坡数目的92.42%、71.21%和36.36%,在同一大型搜索面积不同的分辨率条件下,分别有62、59和52个滑坡点位于极不稳定区域的不稳定区域,分别占滑坡点总数的93.94%、89.39%和78.79%。从分辨率角度来看,正确识别率随着分辨率的升高而升高,而对于搜索面积而言,较大的搜索面积识别的准确性也随之上升。通过对比可以发现，在分辨率较高的情况下,搜索面积的大小对正确识别滑坡的影响不大,而随着分辨率越粗糙,它在识别滑坡的大小方面就有很大的差别。对于30m的分别来说,在不同的搜索面积下,正确识别滑坡的比例相差高达42.43%,而5m分辨率仅相差1.52%。

从图4和上述分析可以得出，DEM分辨率和搜索面积的大小都会影响模型的模拟效果,模型的准确性与分辨率和搜索面积呈正相关关系,分辨率高的DEM不受搜索面积的影响,而分辨率越低,正确识别能力受到滑坡面积大小的影响。

4.2 精度验证

预测精度检验是评价模型预测性能的重要步骤[14]，因此，本文将滑坡编目中的滑坡点和计算机随机选取的非滑坡点作为标准采用可接受者曲线(ROC)中的真阳性率(TPR)和假阳性率(FPR)两个指标来对滑坡点和非滑坡点与模型预测的空间位置进行对比[41]，定量地评估三维模型的预测性能。其中TPR为正确识别的滑坡数目占全部滑坡点的比例,而假阳性率为非滑坡点被错误识别为滑坡的数目占所有非滑坡点的比例。通过这两个指数可以分析实际的滑坡位置点与预测的滑坡的位置是否一致,同时可以比较几种工况下的模型性能。

然而，仅仅通过这两个指标来评价模型的性能是不够客观全面的。对于滑坡敏感性评价来说,当TPR指数过高,说明实际的滑坡点落入不稳定区域的比例较高,然而当整个研究区大量被识别为不稳定区域时,也造成了真阳性率过高的情况,当然在这种情况下,假阳性率也会比较高。因此,为了对结果有一个全面的验证,采用综合性指数LRclass进行结果精度辅助验证,即用每个安全系数等级预测的滑坡点来评价模型的表现能力。LRclass是基于每个安全系数等级中滑坡点的比例,再结合所有的滑坡点的情况和每个等级预测的滑坡面积进行计算[35-36],具体计算方法如式(5)所示:

(5)

式中：S表示每个安全系数等级里滑坡点的比例,A表示每个安全等级的面积占整个研究区的比例。

A 工况1;B 工况2;C 工况3; D 工况4; E 工况5; F 工况6图4 6种工况下的Scoops3D结果及不同稳定性占比图Fig.4 Scoops3D results and proportion diagrams of different stability under six working conditions

为了进一步比较不同分类情况的综合指数,对原公式进行改良,利用%LRclass来计算不同的类在不同的安全系数等级中所占的比例[42-43],可用式(6)计算:

(6)

在进行计算时,为了方便对模型结果的分类,将安全系数值分成了两类,以1为分界线,当F<1代表不稳定坡体,即出现滑坡;当F≥1代表稳定坡体,这种分类方式方便了下文的计算,同时也受到了大家的欢迎。

从图5可以看出所有工况下,所有的点都位于坐标的左上方,即真阳性/假阳性都大于1。即说明了三维模型在模拟浅层滑坡上是可接受的。从纵轴可以看出,分辨率最高的5m,其真阳性率也最高,分别为92.42%和93.94%。这比分辨率最低的30m的真阳性率分别高56.06%和15.15%,说明了DEM分辨率会极大地影响模拟结果的准确性。而从假阳性率来看,DEM分辨率高的错误识别率也越高,5m的分辨率情况下错误识别个数最多,其中在工况2的情况下更是达到了48.33%。这也就是说明在工况2的参数组合下几乎一半的几率会将非滑坡识别为滑坡,这对于一个模型来说并不是一种好的结果。对于高精度的分辨率而言,由于栅格尺寸变小,区域被划分为更多的格子,在计算安全系数时,计算结果偏低,从图4就可以看出其极不稳定区域和不稳定区域面积最多。同时,由于实际地形地貌的影响,一些小型微地貌的存在,例如裂缝、落水洞、侵蚀沟等,导致局部的坡度值增大,在高分辨率情况下将其识别为不稳定区,而正是由于这些因素的存在,使得模型在高精度的DEM计算下安全系数值偏低。这一问题的解决将会大大降低模型的过度识别,提高模型的正确性。而对于同一分辨率下两种不同的工况相比,随着分辨率越来越粗糙,他们的真阳性率相差也越大,而假阳性率差值最大是30m。

评价一个模型的效果,真阳性率越高而假阳性率越低,也就是距离图5中的50%这根线越远越高,而100%的真阳性率与0%的假阳性率这种情况只是理论存在,在实际情况由于误差等原因是不会出现这种情况的。

图5 6种不同工况在Scoops3D模型下的真阳性率和假阳性率Fig.5 True positive rate and false positive rate of six different working conditions underScoops3D model

为了全面评估不同工况下Scoops3D的性能,我们还引入了改良后的综合性指数%LRclass来对模型的正确识别滑坡能力进行比较验证。表2显示,除了工况5,其余的5种工况综合性指数都在80%以上,而5m分辨率的两种工况居然都达到了90%以上,这就说明Scoops3D在滑坡稳定性预测方面具有很好地适用性。工况1的数值最高,且与工况2相差较小,说明在高精度的DEM分辨率下,搜索面积对结果影响不大。而随着分辨率的精度越来越低,搜索面积对模型的影响也越来越大。从分辨率角度来看,毫无疑问,分辨率越高越好,这也与许多人的结论相一致[20-25]。这为以后的研究在选择合适的参数上提供了数字支撑。