基于熵权变异系数融合算法和FLO-2D的朱家沟流域泥石流危险性评价

吴季寰，张春山，杨为民，孟华君，郭 涵，万飞鹏

（1.中国地质科学院地质力学研究所，北京 100081；2.自然资源部活动构造与地质安全重点实验室，北京 100081；3.自然资源部陕西宝鸡地质灾害野外科学观测研究站，北京100081；4.中国地质大学（北京），北京 100083）

引言

青藏高原东北缘地处西秦岭构造带西延［1］，位于南北向地震带与祁连山地震带的交汇区，活动断裂发育，地震活动性强，山地灾害频发。坐落在此的山地城镇，如甘肃岷县、武都等地区，多属泥石流高发地［2，3］，泥石流灾害已严重制约区域经济社会发展。

甘肃岷县朱家沟是青藏高原东北缘典型的小流域泥石流沟，历史上多次发生泥石流灾害，如1982年“8·8”朱家沟泥石流堵断洮河中上游，形成长达15 km的回水区，河水越过防洪堤侵入岷县城区，淹没建筑物1 700余栋，毁坏房屋600余间，造成100余人伤亡，并导致下游约148.67 hm2农田绝收减产。近年来，受2008年汶川地震和2013年“7·22”岷县-漳县Ms6.6地震影响，朱家沟内物源剧增，导致该小流域泥石流危险性急剧升高，极有可能再次爆发大规模泥石流灾害事件，极大威胁流域内及城区居民生命财产安全，因此亟需开展相应的泥石流灾害危险性评估，为防灾减灾提供科学依据。

泥石流危险评价评价方法大致分为多因子叠加法［4-8］和数值模拟办法［9-12］两类，二者区别在于：1）评价方式不同：多因子叠加法以泥石流影响因素为数据基础，对流域进行综合评价分区，直接体现流域整体危险度，数值模拟法基于对流域内所爆发泥石流的泥深与流速的模拟，着重于对泥石流动态过程的反演，并以此为依据进行危险性分区；2）适用范围不同：多因子叠加法适用于大范围区域，而数值模拟办法更适用于小流域或具体沟谷；3）评价目的不同：多因子叠加法旨在评价各流域于未来阶段暴发泥石流的可能性，数值模拟法注重将泥石流强度与重现周期相结合［11］，模拟不同频率下泥石流的运动过程，预测泥石流的运动特征及堆积区域。

本文结合以上两类方法，通过详细现场调查和资料分析，采用熵权与变异系数融合的多因子叠加法开展朱家沟泥石流易发性评价；在此基础上，进行强降雨条件（P=1.6%，P=1%）下朱家沟及其高易发支沟的泥石流运动数值模拟，并与野外科考资料对比验证，最后以流速、泥深、堆积范围为基础进行危险性区划，该研究可为当地泥石流灾害防治提供依据。

1 流域泥石流概况

朱家沟坐落于甘肃岷县城区北部，是洮河一级支沟，流域面积约14 km2，主沟长度约6.7 km，与县城一河之隔。流域位于青藏高原东北缘，东昆仑断裂带和西秦岭北缘断裂带的应变传递和构造转换的中间过渡区，属中强震多发区，基本地震动峰值加速度为0.15 g，沟内哈树—麻路滩断裂为高角度逆冲兼右旋走滑断裂，倾向北东，倾角约55～70°，断层上盘以二叠系炭板岩夹石英砂岩、粉砂岩为主，表层中风化，风化壳厚度2～3 m，属中等硬度块状碎屑岩组，下盘以三叠系板岩、砂质板岩、粉砂岩为主，表层强风化，风化壳厚度5～8 m，属软弱层状碎裂岩。强表层岩土体受地震扰动，结构疏松，崩塌滑坡的触发阈值降低，在地震触发效应下，崩滑发育，沟道中上游淤积严重，具备丰富的泥石流固体物源条件。

朱家沟地势北高南低，属西秦岭北部构造侵蚀剥蚀中低山地貌，海拔介于2 276～3 017 m，相对高差741 m，主沟纵比降约72～110‰，沟道上窄下宽，自北到南曲形延伸，沟道断面自上而下由“V”型向“U”型过渡，两岸分布冲沟33条，支沟17条，其中郎家沟、松树沟、朱林沟和哈固沟是朱家沟的主要支沟。流域地处温带半湿润与高寒湿润气候过渡带，多年平均降水量560.8 mm，降水集中于5～9月，占全年降水总量的78%以上，且多为夏季短时强降雨，受其激发泥石流活动频繁，严重威胁城区及流域内居民生命财产安全。

综上，朱家沟一直是岷县实施泥石流防灾减灾的重要地区。在岷县特大冰雹山洪泥石流灾后重建和岷县-彰县地震灾后重建时，分别在朱家沟、松树支沟、郎家支沟、王家支沟沟口修建排导槽，在松树沟修建拦挡坝3道以及在王家支沟修建拦挡坝1道。

2 基于熵权与变异系数融合算法的泥石流易发性评价

2.1 融合算法原理

信息熵是用于度量事物的不确定性［13］，在地质灾害领域表达为影响因子携带信息熵的大小与灾害发生的概率成反比。通过信息熵计算熵权的过程如下：

（1）建立标准化评价矩阵：

式中m、n、x ij分别代表评价对象数、评价因子数以及第i个评价对象的第j个因子的标准值。

（2）确定第j个因子的信息熵权：

式中e j为第j个因子的信息熵，w1j为第j个因子的熵权。

变异系数是一个表达概率分布离散程度的归一化量度，在评价体系中，某因子的变异系数愈大，则该因子愈能代表体系内部的差异性，对评价结果影响愈大。第j个因子变异系数权重如下：

式中Cv j为第j个因子的变异系数，w2j为第j个因子的变异系数权。

相对熵（式4），可用于衡量两组离散向量之间差异。

根据相对熵的定义，引入融合向量W=(w1,w2,…,w n)与待评价向量组Wp×n，使该向量与熵权法权重向量以及变异系数法权重向量距离之和最小［13］。建立距离规划模型以求解融合权重：

式中p、n的物理意义分别为灾害评价模型数以及灾害评价因子数，w i为第j个因子的融合权重。

2.2 泥石流易发程度的影响因子

从泥石流形成条件：降水、沟谷、物源三方面出发，选取影响泥石流发育的主要因素：形状系数、流域面积、沟谷形心至断层距离、沟谷密度、沟道长度、主沟纵比降、沟口单位宽度清水流量峰值、物源分布、相对高差、流域平均坡度、植被覆盖率，浅表层工程岩组共计12个评价因子，各因子实际取值见表1。

表1 研究区各因子实际取值Table 1 The actual value of each factor in the study area

由于研究区范围较小，同期次降雨下各支沟雨量接近，故降雨量不作为评价因子。但因一次性短历时、强降雨形成泥石流过程中，流域内雨水短期入渗、蒸发较少，多以地表径流形式激发物源启动，并最终从沟口流出，所以沟口处单位宽度清水流量峰值可以直观体现该流域产汇流能力。流域附近历史最大日降雨量为61.5 mm（2010年5月10日），据岷县气象局2020年6～8月逐日降雨量资料，7月18日（以下简称7·18）59.3 mm的降雨量为2020年度当地雨季最大降水量。本文基于FLO-2D的Rain模块，降雨分布如图2，选择适用于西部山区小流域的SKK模型［14］，进行朱家沟“7·18”强降雨-径流模拟，以提取各支沟沟口清水流量过程峰值，见表1。

图2 朱家沟“7.18”降雨曲线图Fig.2 Zhujia gully'7.18'rainfall curve

依据行业规范［15］及相关研究成果［6.16］，将各因子分为：强（IV）、中（III）、弱（II）及不发育（I）四级，并将分级结果代入融合算法计算权重。分级标准与权重计算结果见表2，分级结果见图3。

图3 因子易发程度量化分级Fig.3 Quantitative classification of factors’susceptibility

表2 评级因子发育等级及权重Table 2 Development grade and weight of rating factors

2.3 泥石流易发程度评价结果

经多因子叠加计算，朱家沟内微流域泥石流易发程度评价结果介于［0.424，0.624］，利用统计学中的自然间断点法进行分类（图4）：高易发泥石流沟（0.526，0.624］、中易发泥石流沟（0.464，0.526］、低易发泥石流沟［0.424，0.464］。结合2020年8月野外实地调查分析（见图1），流域内主要支沟：哈固沟、朱林沟、松树沟和郎家沟的均为泥石流高易发沟。

图1 朱家沟流域地质简图Fig.1 Geological schematic map of Zhujiagully watershed

图4 朱家沟泥石流易发程度评价结果图Fig.4 Assessment results of debris flow susceptibility in Zhujia gully

3 基于FLO-2D的朱家沟泥石流危险性评价

综合考虑朱家沟内各支沟泥石流易发程度和获取的ALSO12.5 m数字高程数据（DEM）的精度，选取哈固沟、朱林沟、松树沟、郎家沟等高易发泥石流支沟以及主沟开展评价。

3.1 工况设置及计算参数

FLO⁃2D以非牛顿流体和中央有限差分法为基础，迭代计算流体运动控制方程，模拟洪水、泥石流运动过程［12］。计算条件设置如下：

（1）降雨条件：根据国家气象局提供的1981-2010年岷县逐日降水数据，利用指数外推法对岷县地区年降雨量频率曲线（图5）进行拟合（式6），可知雨量59.3 mm的2020年度朱家沟流域最大单日降雨：“7·18”强降雨，年降雨频率为1.6%，属50年一遇暴雨事件，设置此次降雨作为模拟条件。此外，根据式6选取单日雨量99.38 mm作为百年一遇极端降水条件进行模拟。

图5 岷县地区年降雨量频率拟合曲线Fig.5 Fitting curve of annual frequency of rainfall in Minxian County

式中：R为降雨量，P为年降雨频率，拟合公式置信概率为95%。

（2）有无防治工程条件：历经“5·10”特大冰雹和“7·22”岷县-漳县地震后，朱家沟泥石流危险性大幅增加，且紧邻县城风险较大，定西国土局为有效防治地质灾害，在模拟区内修建工程措施见表3。

表3 朱家沟泥石流治理工程参数Table 3 Parameters of debris flow control project in Zhujia gully

组合降雨、防治工程条件，模拟3种工况下高易发沟道泥石流的危险性，分别为：①无防治工程时“7·18”强降雨（50年一遇）泥石流，②有防治工程时“7·18”强降雨（50年一遇）泥石流以及③有防治工程时百年一遇强降雨泥石流爆发。以工况①、②作为对照组评判工程措施效用，以工况②、③作为对照组分析极端气象下泥石流运动趋势，以极端工况③作为目标组进行危险性区划。

计算所需参数设置如下：1）重度及体积浓度：据李国营等在朱家沟的工作成果，主沟、郎家沟、松树沟、朱林沟和哈固沟的泥石流容重γ（kN/m3）分别为［17.46，17.59，17.53，16.16，17.16］，流域平均容重为17.2 kN/m3，平均比重Gs=2.6 g/cm3，经计算流域内平均孔隙比为e=1.228，体积浓度C V=0.463，对应泥砂比Rns=0.70；2）层流阻滞系数：本文参考侯圣山等［10］研究，用工程地质类比法，取层流阻滞系数K=2 280；3）屈服应力及粘滞系数：根据1）中参数，联立粘滞系数η（Pas）和屈服应力τy（kpa）计算公式以及王裕宜等［17］统一的泥砂比—体积浓度—流变参数关系式（式7），可得［α1，α2，β1，β2］为［0.012 2，0.015 6，17.41，18.28］，以及［η，τy］为［30.146，56.992］。

式中：η为粘滞系数，τy为屈服应力，CV为体积浓度，Rn s为泥砂比，α1、α2、β1、β2为实验系数。

（3）曼宁粗糙系数：利用王裕宜等［17］统一阻力糙率系数表征公式，结合实际调查主沟、郎家沟、松树沟、朱林沟和哈固沟的泥深h（m）为［1.7，1.3，1.3，1.5，1.2］，可得各沟曼宁系数nc为［0.29，0.14，0.14，0.22，0.10］；5）集水点及泥石流流量过程曲线：集水点取自各沟的泥石流形成区底部流通区顶部，根据“7·18”降雨-径流模拟提取的集水点处清水流量线和放大系数B F=1/（1-CV）=1.862可得点位处泥石流流量过程曲线（图6）。

图6 沟谷集水点分布及泥石流流量过程曲线Fig.6 The distribution of catchment points of gullies and flow process curves of debris flow

3.2 模拟结果

图7在实际降雨条件下朱家沟泥石流运动过程的模拟结果。由图7可知：（1）“7·18”强降雨时（图7（a）、（b）、（d）、（e）），无工程条件下最大泥速6.03 m/s，最大泥深4.88～8.07 m，位于沟口堆积扇中部以及沟道汇集处；有工程条件下最大泥速4.54 m/s，三道拦挡坝前最大泥深为7.42～8.58 m，郎家沟至主沟排导槽段受后方补给与前方运动受阻影响，淤积较深可达4.2～5.75 m，沟口堆积区最大深度为3.23 m，拦挡坝与排导槽起到拦截回淤、减速、导流控界的作用。（2）百年一遇强降雨时（图7（c）、（f）），水动力条件增强，雨水润湿沟道，朱家沟内泥速提升、泥深增加，泥石流物质冲出沟口，堵塞洮河，直接威胁岷县县城。

图7 多工况下泥深、泥速分布Fig.7 Distribution of mud depth and mud velocity under different conditions

经验证（表4），模拟精度高于74.12%，能反映“7·18”降雨条件下区内泥石流运动特征与堆积范围，且说明百年一遇降雨条件下泥石流预测较为可信。

表4 模拟结果验证表Table 4 Validation table of simulation results

3.3 泥石流危险性评价

实地调查表明，朱家沟泥石流直进性明显，破坏形式以冲毁、淤埋农田、民居、公路和挤压河道，泥速、泥深与危险性呈正相关关系，参考侯圣山等［10］在岷县耳阳河流域的研究（表5），在百年一遇强降雨条件和2020年“7·18”强降雨条件下，对朱家沟泥石流进行危险性分级（图8）。

表5 泥石流危险性分区标准Table 5 Criteria for hazard zoning of debris flow

泥石流危险性评价结果表明：“7·18”泥石流总危险区面积为371 875 m2，高危险区占危险区总面积的45.72%，主沟中下游均处于高危险区，沟口排导槽将高危险性区限制于槽中，抑制堆积扇成形，此外松树沟内三道拦挡坝处由于起到拦挡回淤作用表现为高危险性；中危险区占危险区总面积的30.70%，分布在松树沟、朱林沟以及主沟上游区域、沟口排导槽出口周围；低危险区占危险区总面积的23.58%，分布于哈固沟和郎家沟。P=1%概率下，泥石流总危险区面积为500 625 m2，与“7·18”降雨条件相比，危险区面积增大34.62%。其中高危险区占总面积71.62%，增加了188 525 m2，具体表现为主沟沟头处危险性上升，朱家沟主沟沿岸房屋建筑易受泥石流威胁，以及沟口冲积扇成型，主要威胁洮河左岸沟口处与右岸30 m范围内城区建筑；中危险区面积约占18.34%，包括朱林沟、哈固沟和松树沟，但由于松树沟内拦挡工程淤满，超过极限状态，使其周边危险性较高；低危险区面积约占10.04%，郎家沟整体危险性较低，沟口排导设施仍可持续倒流。

总体而言，朱家沟泥石流危险性大。其每年雨季向洮河泻流，泥石流物质在沟河交汇处淤积，抬高河床，在极端的气象条件下（图8（b））有阻断洮河，形成堰塞坝的风险，进而会迫使上游水位增涨，淹没岷县城区。并且朱家沟沟口下游2.5 km处为龙王台水电站，沟内流体泻入洮河后，存在部分流体顺流向下汇入水电站库区，挤占库容的可能性。

图8 朱家沟泥石流危险性分区图Fig.8 Zhujia gully debris flow risk zoning map

4 结论与讨论

本文在青藏高原东北部山地泥石流灾害调查的基础上，选取具有区域代表性的泥石流孕灾小流域——朱家沟为研究区，用于探索适合青藏高原东北缘小流域泥石流沟的危险性评价方法。所得方法将依托熵权-变异系数融合算法的多因子叠加模型和数值模拟技术相结合，前者通过衡量泥石流形成的地形、地质条件评价小流域内各沟道泥石流发育程度，判别在降雨诱发下，区内易于暴发泥石流并汇入主沟，参与主沟泥石流流体物质补给，最终共同致灾的支沟；后者则结合前者结果，引入泥石流触发因素，通过模拟极端降雨条件下流域内泥石流的爆发过程和运动特征，进行泥石流危险性区划，二者相互补充，充分考虑青藏高原东北缘小流域的泥石流孕灾条件，用于泥石流危险性评价具有充分的合理性。经2020年朱家沟“7·18”泥石流实例验证，该法具有较高的准确度，用其进行朱家沟泥石流危险性评价得到以下结论和认识：

（1）朱家沟内各支沟泥石流易发程度主要受沟口单宽清水流量峰值、主沟长度、浅表层岩组等条件控制，区内17条支沟中高易发沟6条，中易发沟5条，低易发沟6条，其中哈固沟、朱林沟、松树沟和郎家沟等主要支沟均为高易发程度泥石流沟。

（2）泥石流流速和沟谷淤积数值模拟结果表明，在50年一遇降雨强度以下，朱家沟沟口的排导槽以及松树沟内三座拦挡坝对沟内泥石流起到了导流控界和拦截减速的作用，但基本已达到其能力极限。在百年一遇的极端降雨条件下，对沟口居民和房屋、岷县县城与下游龙王台水电站威胁较大，建议在主沟修建拦挡坝等防治工程，并及时清淤。

（3）“7·18”朱家沟泥石流危险区集中分布于主沟、哈固沟、朱林沟、松树沟和郎家沟沟道，高危险区分布于主沟中下游区域。当遭遇百年一遇强降雨条件（P=1%），朱家沟泥石流危险区面积约500 625 m2，较“7·18”泥石流危险性区面积增加34%，其中高危险区分布于朱家沟主沟、沟口处与洮河右岸30 m范围，对岷县县城构成重大威胁。评价结果可为岷县县城的地质灾害防治、国土空间规划提供支撑。

小流域泥石流危险性区划的影响因素众多，从泥石流高易发支沟的选取到模拟时所选定的极端降雨条件均直接影响评价结果。本文虽建立起包含了12个地形地质方面致灾因子的泥石流易发程度评价体系，但忽略了因子相关性，未能将高度相关因子剔除或整合，一定程度上影响了高易发支沟的选择，同时，仅重现了百年一遇降雨条件下的泥石流过程，难以全面反映爆发不同规模泥石流时流域的危险性分布。因此，为得到更准确的评价结果，在下一步的工作中应重点提炼泥石流致灾因子，凝练评价体系，以及开展不同雨强条件下小流域泥石流危险性综合区划。

致谢

本文在撰写过程中使用了岷县气象局提供的朱家沟流域降雨资料以及岷县岷阳镇朱家沟泥石流治理工程勘察资料，在此表示诚挚的感谢。