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

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

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

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 危险性评价基础数据Fig.2 Basic data of risk assessment

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

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

W=Hsdθρw+1-φρrgcosβsinθ, (1)

σN=Hsθρw+1-φρrgcos2β, (2)

τS=csoli+{Hsdθρw+1-φρrgcos2β-χhρwg}tanγ。 (3)

式中: H_sd为土壤厚度, m; θ 为体积含水量, 无量纲; h为毛细管头, m; ρ _w为水的密度, kg/m³; ρ _r为土体密度, kg/m³; β 为坡角, (° ); ϕ 为孔隙率, %; C_soli为土壤黏聚力, kPa; χ 为剪切强度与毛细管压力的比例因子, 无量纲; γ 为土体内摩擦角, (° )。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 六棱柱触发模型基本单元[23]Fig.3 Basic units of triggering model of hexagonal cross section[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

Tab.1

表1(Tab.1)

表1 土体强度参数 Tab.1 Soil strength parameters 黏聚力/Pa 根系强度/Pa 内摩擦角/(° ) 初始含水率/% 11 000 200 11 37

表1 土体强度参数 Tab.1 Soil strength parameters

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

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

表2

Tab.2

表2(Tab.2)

表2 研究区不同频率下暴雨统计参数及累积雨量 Tab.2 Statistical parameters and accumulated rainfall of rainstorm at different frequencies in the study area 频率 1/6 h 1 h 6 h 24 h 平均 值/mm 变异 系数 模比 系数 设计雨 强/mm 平均 值/mm 变异 系数 模比 系数 设计雨 强/mm 平均 值/mm 变异 系数 模比 系数 设计雨 强/mm 平均 值/mm 变异 系数 模比 系数 设计雨 强/mm 1% 15 0.40 2.58 34.63 35 0.43 2.62 85.20 60 0.41 2.95 140.90 68 0.40 2.62 157.00 2% 2.28 31.23 2.32 76.30 2.55 126.80 2.32 141.60 5% 1.89 26.63 1.91 64.40 2.06 107.70 1.91 120.70 10% 1.59 23.30 1.61 55.51 1.68 92.80 1.61 104.30

表2 研究区不同频率下暴雨统计参数及累积雨量 Tab.2 Statistical parameters and accumulated rainfall of rainstorm at different frequencies in the study area

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

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

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

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

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 研究区土层深度分布Fig.4 Soil depth distribution diagram in the study area

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

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 研究区不同频率降雨历程模拟概化曲线Fig.5 Simulation generalization curves of rainfall processes at different frequencies in the study area

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

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 基于无限边坡模型的不同降雨条件下的稳定性指数Fig.6 Stability index under different rainfall conditions based on infinite slope model

根据图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

Fig.7

	Figure Option ViewDownloadNew Window
	图7 不同降雨条件下稳定性指数栅格值频率统计直方图Fig.7 Raster frequency statistics histogram of stability index under different rainfall conditions

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

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 基于K值聚类分析的8· 31山洪泥石流事件危险性分级Fig.8 Risk classification of 8· 31 flood & debris flow event based on K-value clustering analysis

表3

Tab.3

表3(Tab.3)

表3 K值聚类分析 Tab.3 Clustering analysis of K value 栅格值分类 聚类中心 1 2 3 4 栅格值-4类 0.047 4 0.289 8 0.515 8 0.939 7 栅格值-3类 0.047 9 0.311 2 0.932 1

表3 K值聚类分析 Tab.3 Clustering analysis of K value

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

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 研究区2012年8月31日降雨条件计算结果Fig.9 Calculation results of rainfall conditions in the study area on August 31, 2012

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

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

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 滑坡区稳定性指数验证Fig.10 Stability Index verification of landslide

此外, 研究区易发性评价也可采用斜坡单元进行表示, 因此需将基于栅格的危险性评价图运用ArcGIS中的分区统计进行分区。可根据需要选取平均值、众数、最大值等指标, 本文通过对比现有地质灾害分布位置, 采用众数作为危险性分区的指标, 进行斜坡单元的危险性等级划分。