基于数值模拟的城镇地质灾害危险性评价方法

2022-09-02龚凌枫铁永波卢佳燕高延超

中国地质调查 2022年4期

龚凌枫，徐伟，铁永波，卢佳燕,2，张玙，高延超

(1.中国地质调查局成都地质调查中心，四川成都 610081； 2.中国地质大学(北京)，北京 100083)

0 引言

随着城镇化水平的不断提升，对重点区城镇地质灾害风险精细化评价需求及评价结果的准确性要求进一步提升[1-2]。当前危险性评价方法主要包括定性评价方法(专家打分法[3]、模糊综合评判[4-7])、半定量评价方法(信息量法、层次分析法[7-8]、证据权法[9]、频率比法、逻辑回归模型[10]、判别分析模型、聚类分析模型[11]、人工智能模型[12-13])及定量评价方法(无限斜坡模型[14]、3D极限平衡方法、斜坡水文模型[15])，上述方法都是地质灾害危险性评价中常用的分析方法。

定性评价方法需丰富的灾害点数据支撑，且对灾害点数据的分布均匀性有较高的要求，否则难以避免因调查范围、精度差异造成的评价结果与灾害点分布的自然属性不符，该评价方案常用于小比例尺的地质灾害危险性评价[16]。定量评价方法是以岩土体力学性质、气象及地表水文、地下水渗流等[6,17-20]与灾害形成直接相关的特性为基础，针对不同区域地质灾害主控因素辅以对应模型进行计算，适用于重点调查区精细化地质灾害危险性评价，可实现大比例尺地质灾害危险性评价。

目前，较多的地质灾害单体危险性评价应用为重要灾害点危险性评价提供了较多的案例，但区域性危险性评价较少采用数值计算方法，主要是受到群发性地质灾害评价方法的制约，而单体地质灾害应用于区域性评价不可避免造成工作量剧增，因而难以推广应用。浅层滑坡、群发性泥石流优势发育地区的数值计算方法较为成熟[15-16]，可以用于精度要求高的重点区和典型城镇地质灾害危险性评价。

本文以典型红层地区浅层滑坡为研究对象，选用滑坡流体力学触发模型(hydromechanical triggering of landslides,LHT)进行重点研究区地质灾害危险性评价。2012年8月31日，四川凉山州喜德县遭遇百年一遇的特大暴雨，孙水河上游米市监测当日最大降雨量达到149.2 mm；强降雨发生后，在孙水河流域两河口以上流域的山坡多处可见“山剥皮”现象，大量植被被破坏，对生态环境产生较大的负面影响。本文以“8·31”群发性滑坡泥石流事件中灾害解译点为验证数据，探索不同降雨频率条件下重点区浅层滑坡危险性评价方法及基于数值模拟的不同降雨频率下1∶1万地质灾害危险性定量计算方法。

1 研究区概况及数据来源

研究重点区为四川省凉山州喜德县洛哈镇的集镇区，面积44.5 km2，其地理位置、地质概况及地质灾害点分布如图1所示。研究区具有冬季干燥无严寒，夏季温凉多雨，四季不分明，气温日差大，年差较小，风多，夜雨多，冰雹多的典型亚热带季风和高原气候特征。

图1 研究区地质概况及地质灾害点分布

1.1 工程地质概况

区域内主要构造为米市向斜，轴线经过米市镇、洛哈镇，呈NE27°，南起额尼村，向北延伸。向斜轴部及两翼地层均为上白垩统雷打树组上段(K2l2)砂质泥岩。东翼倾向北南270°～318°，倾角2°～19°；西翼倾向南东120°～145°，倾角10°～15°。局部地形起伏大，沿路及沿河陡坡区域卸荷裂隙发育，岸坡风化层厚度大，为小规模的崩塌及大量的浅层滑坡提供了良好的地形条件，且由于该区域植被覆盖度较差，尤其是海拔较高的位置山地裸露，为长历时及强降雨作用下浅层滑坡的发生提供了的条件。

研究区岩性为灰黄色质泥质砂岩和紫红色泥岩。河流及筑路切坡处陡坡基岩以泥岩为主，夹薄层砂岩。强风化基岩破碎呈碎块状、中风化基岩呈短柱状，强度较低。人类居住区和耕值区的地形受人类活动影响改变较大，主要表现为农作物耕植过程中在坡体上开挖形成梯田状小土坎，长度一般为8～15 m，坡高一般为0.5～1.5 m。总体来说，区内人类工程活动较强烈，对地形的改变为滑坡的滑动提供了有利的变形空间。

1.2 数据来源

研究区数据主要包括研究区高分二号卫星数据(分辨率为0.8 m，拍摄时间为2020年)、1∶10 000地形图、1∶50 000地质图、1∶50 000地质灾害分布图、综合遥感解译图及野外地质调查资料等，可通过上述资料获取该区域危险性评价基础图件(图2)。通过前期资料收集及野外调查，区内原有灾害点3处，初步光学遥感解译及InSAR识别获取滑坡点17处，通过精细光学遥感解译共获取浅层滑坡点409处，滑坡面积共计1.07 km2，占研究区总面积的2.4%。

图2 危险性评价基础数据

2 研究方法

本文以浅层滑坡发育区地质灾害危险性精细化评价为目标，选取数值模拟方法进行评价。通过对国内外已有成果分析可知，危险性精细化评价主要在于岩土体稳定性计算方法选取及不同频率条件下的降雨历程曲线计算。

2.1 数值计算方法选取

对重点区调查采用基于栅格计算斜坡稳定性的LHT模型[21]。LHT是一种基于力学公式的数值模型，能完整模拟坡体在降雨条件下土体发生渐进性破坏，最终整体失稳的过程。模型把覆盖在基岩上的土体离散成一组机械连接的六角形土柱[22]，通过物理模型将土体的水文状态和力学状态联系起来，以摩擦力、黏聚力、毛细力及植物根强度来保持坡面土体的稳定性。

六角土柱的基本受力模型如图3所示，在降雨条件下，单个六棱柱的力学参数包括土柱体积重度W、正应力σN和不饱和抗剪强度τS,计算公式为

W=Hsd[θρw+(1-φ)ρr]gcosβsinθ，

(1)

σN=Hs[θρw+(1-φ)ρr]gcos2β，

(2)

τS=csoli+{Hsd[θρw+(1-φ)ρr]gcos2β-
χhρwg}tanγ。

(3)

式中：Hsd为土壤厚度，m；θ为体积含水量,无量纲；h为毛细管头，m；ρw为水的密度，kg/m3；ρr为土体密度，kg/m3；β为坡角,(°)；φ为孔隙率,%；Csoli为土壤黏聚力,kPa；χ为剪切强度与毛细管压力的比例因子,无量纲；γ为土体内摩擦角,(°)。

图3 六棱柱触发模型基本单元[23]

该模型将滑坡过程建模为一个瞬时事件，没有去描述滑坡土体运动的具体过程，也没有跟踪滑坡的质量变化，通过由微观到宏观的力学描述，实现了降雨条件下滑坡的触发模拟。本项目采用STEM TRAMM软件进行计算，需要输入的数据包括栅格数据(土层厚度、植被覆盖率、地形数据)和土体强度指标。

(1)土体厚度。以研究区地形数据为基础数据，分别选取高程、坡度(以坡降计)、地形湿度指数(TWI)、距水系距离、方位、曲率，综合计算土体厚度栅格图。在研究区，选取因防火通道开挖而形成的开挖边坡进行线路调查，共选取了162个点统计其土层深度，其中81个点用于栅格深度校准，另81个点作为效果分析点位。

(2)植被覆盖度。植被发育为基于Landsat 8影像提取的归一化植被指数(NDVI)图，利用公式NDVI=(IR-R)/(IR+R)计算NDVI值，其中IR和R分别为近红外(NIR)和红色波段。NDVI数据栅格大小为10 m×10 m。

(3)土体强度。参考文献[24]及区内21份勘察报告实验结果，综合选取用于数值计算的参数，主要有土体黏聚力、根系强度、内摩擦角及初始含水率，对应参数取值见表1。

表1 土体强度参数

2.2 降雨特征曲线

根据前文阐述，降雨条件需考虑10 a、20 a、50 a及100 a一遇降雨工况，但受制于降雨监测数据精度及历史的制约，雨量数据常获取困难，无法根据要求选用不同频率的日降雨量参与分析评价。因此，对降雨数据的分析主要需探索不同频率下的降雨量及降雨历程曲线。

关于降雨量，考虑到四川省内实际情况，各研究区可直接参考《四川省暴雨统计参数图集(2010年)》中所附的暴雨量等值线图，其他地区可参照地方图集或全国暴雨参数图集。表2给出了研究区中米市—洛哈镇1/6 h、1 h、6 h、24 h多年暴雨均值、变异系数等参数，并通过查询皮尔逊Ⅲ型曲线得到不同频率下的模比系数(Kp)，求得不同频率下的雨强值。

表2 研究区不同频率下暴雨统计参数及累积雨量

不同的降雨过程影响地质灾害发育分布特征及灾害规模[25]，且上述数值模型的运用需将降雨过程纳入模型的运算中，因此，评价研究区降雨的过程特征是进行不同区域、不同频率降雨诱发地质灾害研究的基础。根据对地质灾害发生与降雨过程的响应过程分析，在该区域典型的3次灾害事件中，中期集中降雨诱发地质灾害的特征最为明显，且马厂沟1985年7·1泥石流、1973年6·30群发地质灾害及1983年7·1群发地质灾害均为6 h左右。因此，针对该区域降雨特征，将降雨概化曲线设定为降雨历时为6 h的中期集中降雨型事件。

本次评价主要选取累计降雨量(Ra)—持续时间(t)的关系式[26-27](Ra=-1.804 9t3+2.317 3t2+0.471 9t)来代表本区域诱发灾害事件的降雨雨型，累计雨量如表2所示。

此外，也可根据地区特征，选取国家雨量站数据(https://rp5.ru；http://data.cma.cn/)提取3 h或1 h降雨数据产品，经过处理获取月累计平均降雨量、日降雨量等数据，进而通过统计学方法获取目标区域的降雨历程曲线。工作区降雨量在空间上存在显著差异，需进一步考虑该指标随高程变化的空间差异性。由于本研究区高差较小，所以研究中未考虑降雨的空间差异性。

3 计算结果与效果分析

土体深度及降雨过程曲线是数值计算的基础，以下分别对土体深度计算、土体深度数据分析及降雨过程曲线进行分析，并基于数值模拟计算浅层滑坡稳定性指数。

3.1 土体深度计算与数据校准

通过提取野外调查线路上的栅格值与对应调查位置处的土体深度，对比校核获取研究区土层深度的栅格图(图4(a))，并选取剩余81处调查点进行精度分析(图4(c))。根据图4(a)中深度分布栅格图可知，研究区坡体土层深度范围为0.05～10.18 m，采取调查路线上对应点的栅格数据可获取图4(c)中灰色虚线所示的深度值。栅格图中部分点土层深度波动较为剧烈，通过移动平均法绘制黑色趋势线。选取图4(b)中因消防通道开挖揭露良好的土层深度值，结合实测点深度标准误差绘制散点图(图4(c))，结果表明研究区栅格图能够较好地反映研究区坡体处土层的深度分布特征。

(a) 研究区土层深度栅格图 (b) 土层深度调查路线 (c) 土层深度及误差分析

3.2 降雨历程曲线计算

结合安宁河流域降雨诱发灾害及喜德县历次群发性地质灾害降雨历程可知[24, 28]，喜德县群发性浅层滑坡均发生于前期长历时降雨之后的短时强降雨期间。对比前述雨型分布特征的描述，喜德县群发性地质灾害发生条件与式(1)中雨型分布特征更为一致，因此，不同降雨频率下的危险性分析均在此雨型条件下设定。根据研究区不同频率下6 h累计降雨量所计算的降雨分布过程曲线如图5所示。

图5 研究区不同频率降雨历程模拟概化曲线

3.3 浅层滑坡稳定性指数及分析方法

根据前文中所示的降雨历程输入数值计算模型中，结合易发性评价中地形(DEM)、覆盖层厚度、归一化植被指数(NDVI)及岩土体力学参数(表1)进行降雨条件下的稳定性计算，将计算结果中基于栅格的稳定性指数稳定(0)—不稳定(1)作为不同降雨频率下的危险性评价值，评价结果如图6所示。

根据图7所示的评价结果统计可知，研究区在10 a一遇降雨频率下稳定性指数的平均值为0.42，标准差为0.32； 20 a一遇降雨频率下稳定性指数的平均值为0.44，标准差为0.33； 50 a一遇降雨频率下稳定性指数的平均值为0.58，标准差为0.37； 100 a一遇降雨频率下稳定性指数的平均值为0.58，标准差为0.37。可见，在图5所示的降雨频率下，随着降雨强度及降雨量的增加，不稳定指数的平均值也在增加，表明不稳定坡体不断增加。但随着降雨量进一步增加到100 a一遇时，稳定性指数平均值及标准偏差均无变化，同时，从图7中也可以看出评价图的指数均无变化。这说明研究区在50 a一遇降雨条件下即可导致区域不稳定坡体达到较高含水量出现垮塌、溜滑。因此，50 a一遇降雨值可以作为区域浅层滑坡、泥石流及流域山洪的参考阈值。

图7 不同降雨条件下稳定性指数栅格值频率统计直方图

4 讨论

上述基于不同频率降雨的浅层滑坡危险性评价是针对降雨统计的预测值，其准确性可以通过研究区典型的山洪泥石流事件后浅层滑坡分布位置与对应降雨条件下稳定性指数计算进行有效性评价。以下以2012年“8·31”洪水泥石流事件为实例探讨上述方法的准确性。

4.1 “8·31”山洪泥石流事件危险性评价

解明礼等[29]通过收集易发性相关文献中具有代表性的200条成果进行统计，分析不同分级方法的等级划分准确性。由图8可知，基于数值模拟获取的不同降雨条件下稳定性指数具有聚类的特征，运用K值聚类分析方法进行分级也取得了较好的效果(表3、图8)。

图8 基于K值聚类分析的8·31山洪泥石流事件危险性分级

表3 K值聚类分析

由浅层滑坡解译的数据来源可知，本文浅层滑坡主要是在2013年3月的光学影像基础上通过目视解译获取，基本反映了“8·31”孙水河流域山洪泥石流事件中浅层滑坡的发育特征。因此，将2012年8月31日降雨及其形成的山洪泥石流事件作为实例，利用其降雨曲线(图9(a))进行计算，获取如图9(b)所示的此次降雨条件下的地质灾害危险性评价结果。

(a) 降雨历时曲线

(b) 地质灾害危险性评价

4.2 精度分析

合理的危险性评价结果应实现现有地质灾害点分布于高危险区、极高危险区，而非滑坡灾害点主要位于低易发区。量化检验评价结果的方法主要有AUC值、Kappa系数、经验概率值、总体精度、混淆矩阵及历史数据论证等方法。AUC值评价法在危险性评价中较常用，且依赖于2012年“8·31”山洪泥石流事件后的浅层滑坡解译点作为样本，具有较高的可信度。因此，本文选取AUC值探讨“8·31”山洪泥石流事件中浅层滑坡点的评价准确性。

分别提取滑坡解译面文件对应栅格值及非滑坡区随机生成的测试点栅格值，图10给出了滑坡验证成功点的验证曲线，基于“8·31”浅层滑坡灾害点的验证曲线表明该事件中AUC值为0.65，具有较高的准确性。由图9(c)可知，受栅格数据的原始精度、光学影像精度、光学影像拉伸等因素的影响，较大区域数值模拟结果与实际位置对应性受到限制。故通过数值计算表达了区域大部分灾害点发生位置已具有良好的适用性。