基于小流域单元的泥石流易发性评价

2022-11-16陈攀葛永刚孙庆敏梁馨月

科学技术与工程 2022年29期

陈攀，葛永刚, 孙庆敏, 梁馨月

(1.中国科学院水利部成都山地灾害与环境研究所，成都 610041； 2.西藏大学工学院，拉萨 850000； 3.中国科学院山地灾害与地表过程重点实验室，成都 610041； 4.中国科学院大学，北京 100049)

泥石流通常认为是含有大量松散固体物质的固液两相流[1]。因其具有巨大的体积和极快的流速，极易损坏附近桥梁、道路和房屋建筑，威胁当地居民的生命安全。近些年数次重大泥石流灾害事件[2-4]给国家和人民带来巨大的财产损失。进行泥石流易发性评估能识别出研究区域内的泥石流高易发区，为工程建设和防范治理提供决策依据。

评价单元、指标因子和评价模型作为影响泥石流易发性评估结果的重要因素，不同学者已对其进行了讨论研究。邹强等[5]通过对比水文响应单元与栅格单元，发现利用水文响应单元进行泥石流易发性评价具有诸多优点。此后，基本选择小流域单元作为泥石流易发性评价单元[6]。易发性评价指标基本为孕灾环境因子，且主要围绕地形地貌、构造活动、径流水系进行选择[7-8]，由于区域间存在差异性，所以并没有形成统一的指标因子选择标准。泥石流易发性评价模型按评估方式不同可分为四类，包括主观专家分析[9-10]、统计分析[11-15]、机器学习[16-19]和物理模型[20-21]。主观专家分析具有较强主观性，所得评估结果精度较低。统计分析较主观专家分析更加客观，但未反映评价指标与泥石流之间复杂的非线性关系[22]。机器学习对训练样本要求较高，且易存在过拟合。物理模型需要足够的物理力学参数以保证评估的准确性，通常适合小范围的易发性评价。目前，主要通过耦合主观专家分析、统计模型和机器学习的方式提高泥石流易发性评估精度，缺乏对指标因子选择的讨论。对此，现通过对比有无坡体稳定性因子，分别建立考虑坡体失稳破坏的试验组指标体系和仅有孕灾环境因子的对照组指标体系，采用确定性系数模型对安宁河流域进行易发性等级划分，不仅讨论指标因子对泥石流易发性评价结果的影响，同时也耦合统计模型与物理模型，增加泥石流评估结果的物理意义。

1 研究区概况及数据来源

1.1 研究区概况

安宁河流域地处四川省西南部，位于第一阶梯与第二阶梯过渡区内，整个流域约11 050 km2，自北向南依次主要分布于冕宁县、喜德县、西昌市、德昌县、会理县和米易县，海拔高度介于951 ～ 5 137 m。地貌类型包含断陷盆地、河谷平原、构造剥蚀中山和冰蚀高山[23]。安宁河自北向南流域内河谷地带分布大量第四系松散物，河谷东侧主要分布沉积岩，河谷西侧主要分布岩浆岩。受构造影响安宁河流域内发育有控制性的安宁河断裂和则木河断裂，滑坡、崩塌沿断层带分布，为泥石流形成提供充足物源。研究区属亚热带季风季候，年平均降雨量达1 000 mm以上[24]，为泥石流发育提供了水源条件。此外，人类火烧和不合理工程活动降低了坡体失稳和泥石流形成条件，使安宁河流域内泥石流活动频发。

1.2 数据来源

研究所需数据来源如下：①12.5 m分辨率的ASTER DEM数据提取地形地貌和水系发育信息；②1∶250 000地质图提取安宁河流域地层岩性和构造断层信息；③1∶4 000 000 安宁河流域地震峰值加速度区划图；④Landsat 8卫星数据提取植被覆盖度信息；⑤基于坡体稳定性分析获得潜在失稳面积与小流域面积比数据。

2 泥石流易发性评价

2.1 小流域单元提取

小流域包含泥石流形成和运移的全过程，基于小流域获得的评估结果更能满足实际需求。因此，利用ArcGIS软件中的水文分析模块，经过洼地填充、流向提取、流量计算、河网定义和集水区生成5个步骤，划分安宁河流域范围的小流域单元。如图1所示，共得到491个小流域单元，小流域面积介于2.38～179.48 km2，小流域平均面积为22.52 km2。

2.2 评价指标体系构建

2.2.1 评价指标

泥石流易发性评价因子通常围绕物源条件和地形条件进行选择。结合安宁河流域实际情况，选取地表起伏度、地表粗糙度、沟壑密度、距断层距离、地震峰值加速度和潜在物源分布面积与流域面积比作为试验组指标因子，并选择与之对应的孕灾环境因子构建对照组评价指标体系，对照组指标因子包含植被覆盖度、地层岩性、地表起伏度、地表粗糙度、沟壑密度、距断层距离和地震峰值加速度(图2)。所选评价指标可分为孕灾环境因子和坡体稳定性因子。

图2 泥石流易发性评价因子

孕灾环境因子植被覆盖度从水文机制和力学强度两个方面影响坡体稳定性，同时植被覆盖度对水土流失也具有重要影响。不同的地层岩性具有不同的坡体结构强度，是常见的基础地质因子。地表起伏度和地表粗糙度均可表征地表侵蚀程度和地形的起伏变化，前者同时体现小流域物质运移具有的位置势能大小，后者则从微观反映地表单元地势起伏的复杂程度。沟壑密度描述了小流域范围内水道发育程度，反映了小流域内受到一定规模水流冲刷的区域所占比例。距断层距离是内动力地质构造指标因子，通常断层活动会加剧岩体破坏，影响崩塌和滑坡分布。地震峰值加速度反映了地震活动情况，表征了坡体失稳激发因素的影响。

坡体稳定性因子潜在物源分布面积与小流域面积比表示单位小流域面积内的潜在失稳面积。在保证计算精度的前提下合理简化采样及计算过程，通过根-土分区的方式对坡体稳定性进行分析，具体过程见流程图(图3)。考虑植被覆盖度和土体饱和度的影响，利用安全系数计算公式[25]，计算降雨饱和状态下流域内坡体安全系数，以1.5作为稳定坡体与潜在失稳坡体阈值[26]，划分出安宁河流域坡体潜在失稳区。通过对比流域内地质灾害分布情况和坡体潜在失稳区计算结果，发现从区域上看两者具有相同的分布规律。

图3 潜在失稳区计算流程图

2.2.2 多重共线性诊断

泥石流易发性评价指标体系要求评价因子间互不干扰，通常通过共线性检验结果进行判断，当方差膨胀因子(variance inflation factor，VIF)大于10时，表示指标因子之间存在多重共线性，反之，则所选指标因子间相互独立。

通过SPSS软件对所选指标因子进行共线性诊断，容忍度(tolerance, TOL)和方差膨胀因子(VIF)结果如表1所示。由表1可知，试验组和对照组所选各指标因子的VIF值均小于4，说明各指标因子之间基本相互独立，可参与模型进行泥石流易发性预测。

表1 多重共线性诊断结果

2.3 确定性系数模型

确定性系数模型为概率函数，基于泥石流灾害形成条件不变的假设，可计算不同指标因子对泥石流发生的敏感性，现已作为泥石流易发性评价的常见方法之一[27]。其计算公式为

(1)

式(1)中：CF为确定性系数值；PPa为评价因子a中的条件概率，可用评价因子a中的泥石流灾害个数或面积与评价因子a面积比值表示；PPs为泥石流灾害在研究区的先验概率，可用研究区的泥石流灾害个数或面积与研究区面积比值表示。

2.4 易发性评价模型计算结果

从198条泥石流沟中随机选择80%的样本数据作为训练样本，参考常见断点方式和自然断点法利用确定性系数模型计算不同指标因子的分级CF值。计算结果如表2所示。

表2 各指标因子不同分级CF值计算结果

2.5 泥石流易发性评估

基于确定性系数模型计算结果，利用ArcGIS将研究区内各评价因子指标的CF值叠加，利用自然段点法将安宁河流域泥石流易发性划分为高易发、较高易发、中易发、较低易发和低易发5个易发性等级区，结果如图4所示。

图4 安宁河流域泥石流易发性分区

2.6 泥石流易发性评估结果验证

利用频率比值对泥石流易发性结果进行检验，结果如表3所示。由表3可知,试验组指标体系与对照组指标体系的泥石流易发性高易发区与较高易发区的频率比值之和分别占各自总频率比值的88.17%和85.88%，均大于85%，表明两种指标体系均有效评价了安宁河流域泥石流易发性。

表3 泥石流易发性分区频率比值

以泥石流灾害易发性评价中代表研究区内各易发性等级面积累加百分比为横轴(假阳性率)，泥石流灾害易发性评价中代表研究区各易发性等级内真实发生泥石流灾害的累加百分比为纵轴(真阳性率)，绘制的接受者操作特征(receiver operating characteristics，ROC)曲线如图5所示，对照组指标体系的泥石流易发性ROC曲线下面积AUC值为0.752，试验组指标体系的泥石流易发性AUC值达0.767，与对照组指标体系相比预测精度值提升了1.99%。表明考虑坡体稳定性分析的泥石流易发性评估结果较仅考虑孕灾环境因子的结果好。

图5 不同指标体系ROC曲线

图6中，中坝村泥石流沟和老洼沟泥石流沟均在试验组指标体系泥石流易发性评价结果的较高易发区内，表明泥石流易发性评估结果与实际泥石流发育吻合度较高。

图6 实例验证

3 讨论

坡体稳定性指标耦合了地形和力学强度参数，考虑了孕灾环境因子与泥石流之间非线性关系。对照组指标体系中的地表粗糙度与坡度间存在共线性，为保证评估结果的合理性同时与试验组指标体系形成对照，综合考虑选择地表粗糙度作为指标因子。试验组指标体系则将地表粗糙度(地形地貌因素)和坡度(物源稳定性因素)在相互独立的条件下同时引入，使构建的指标体系更加全面。试验组指标体系相较于传统孕灾环境指标体系加入坡体稳定性因子，增加了泥石流形成机理分析，提高了泥石流易发性评估精度。