基于动力过程的单沟泥石流危险性评价方法

2023-11-04谢奎林陈兴长陈慧唐勤冯鑫

科学技术与工程 2023年29期

谢奎林, 陈兴长, 陈慧, 唐勤, 冯鑫

(西南科技大学环境与资源学院天府新区创新研究院, 绵阳 621010)

泥石流是一种含有大量泥沙、石块的特殊洪流,具有突发性强、破坏力大、流速快、历时短等特点,常常造成大量的人员伤亡和财产损失,是影响中国山区建设的主要地质灾害之一[1-2]。2010年8月13日四川省绵竹市清平乡走马岭沟暴发特大泥石流,造成379户农房被毁,致使绵远河改道并淹没清平老场镇,直接经济损失达3 000万元[3];同年8月14日,强降雨诱发汶川县映秀镇红椿沟暴发泥石流,堵断岷江主河,导致河水改道冲入映秀镇新城,造成13人死亡、59人失踪,数千人被迫避险转移[4];2016年5月8日,福建泰宁县金溪左岸的芦庵坑沟暴发特大规模的泥石流灾害,导致36人死亡或失踪[5]。2020年6月17日,四川省甘孜州丹巴县梅龙沟暴发大规模泥石流,阻断小金川河形成堰塞湖,造成严重的次生灾害链[6]。泥石流灾害引起了众多学者的广泛关注,泥石流灾害的研究和危险性评价一直是灾害研究领域的重点。

泥石流危险性评价分为区域评价和单沟危险性评价。区域评价主要是查明区域内泥石流的危险性,筛选重点研究的泥石流沟。单沟评价应以泥石流危害为主要评价内容。目前,这两种尺度的评价方法很多是类似的,如模糊综合评价法、多元回归分析法、信息量法以及灰色关联分析法等[7-10]。这些评价方法虽然各不相同,但其原理基本一致,从泥石流形成条件出发(如泥石流成因,孕灾环境以及导致其暴发的诱发因子等)进行评价[11],但这更像是泥石流易发性评价(泥石流发生的可能性)。泥石流的危险性来自其危害,应是泥石流发生的结果,而不是原因。因此,泥石流危险性评价应是特定降雨频率下泥石流危害范围和大小[12]。

此外,这些方法得到的评价结果过于笼统和概化,往往一条泥石流沟得到一个危险等级。如果泥石流的危险性来自其危害,那么在一个流域内的不同部位,其危险性等级应该是不一样的。从泥石流形成条件出发的危险性评价方法无法对不同部位进行危险性分区。泥石流堆积区往往是山区人民生活和生产的重要场所之一[13],精细化的危险性分区,对泥石流滩地的合理开发和利用,以及泥石流防灾减灾都具有重要意义。

随着泥石流动力学研究和计算机技术的发展,泥石流动力过程的数值模拟得到了快速发展,为精细化的泥石流危险性评价提供了一种全新的手段。韩玫等[14]利用CFX软件模拟再现了桃关沟泥石流泛滥过程,得到了泥石流动力学参数和相应的危险区范围;刘福臻等[15]使用FLO-2D模拟了江达县典型公路泥石流,合理地预测了公路泥石流的影响区域和危害程度。

尽管目前已经有很多学者使用数值模拟的方法对泥石流进行危险性评价,但评价指标和分区标准不统一,有的以泥深或流速作为评价指标来进行危险性分区,也有的以动量、动能为依据进行危险性划分[16-17]。

基于此,现提出基于泥石流动力学参数(流速、泥深、冲击力等)的危险性评价方法。通过对泥石流危险区动力学参数的计算和评价,结合动力学分区指标,实现单沟泥石流危险性的评价分区。利用该方法对2010年甘肃舟曲“8·7”特大泥石流进行危险性评价和分区,以验证方法的可行性和实用性,为泥石流防灾减灾和灾害精细化管理提供支撑。

1 泥石流危害方式

1.1 冲击危害

泥石流的危害方式主要表现为冲击和淤埋,以冲击破坏最为剧烈,往往给流通路径上的建筑物造成毁灭性破坏。如,1981年成昆铁路利子依达沟暴发特大泥石流,强大的冲击力将直径4 m 多的圆柱形桥墩拦腰截断,当时正在行驶的客车跌入泥石流沟中,数百人死亡[18]。冲击危害与泥石流的冲击力有关,冲击力越大,冲击破坏就越严重。泥石流的冲击力由浆体产生的整体动压力和大石块撞击力构成,其中大石块的瞬时冲击会对建筑物结构造成严重破坏,可能导致建筑物整体垮塌。这两种冲击力在运动过程中往往相互叠加,其破坏能力巨大。

为了计算泥石流冲击力,学者们建立了多种计算模型。其中,浆体冲击力计算模型主要有动力学模型[19-20]、静力学模型[19,21]、动量-冲量模型等[22];大石块冲击力计算方法则主要包括梁式计算法[23]和Hertz弹性碰撞理论计算法[24]等。常见的泥石流冲击力计算方法列于表1。可以看出,除静力流体模型外,泥石流的冲击力都与流速密切相关。浆体的冲击力与流速的平方成正比,而大石块的冲击力与流速成正比。研究表明,随着泥石流浆体和大石块运动速度的增加,会导致砌体结构和框架结构建筑物产生严重位移和破坏,其中大石块对简支梁的最大冲击力更是近似线性增长[25-27]。因此流速是泥石流冲击力大小的主要决定性因素,也是衡量冲击破坏能力的重要参考。

表1 泥石流冲击力计算公式Table 1 Calculation formula of debris flow impact force

1.2 淤埋危害

淤埋危害是指泥石流流经建筑物周边或从建筑物的开口进入内部时,在其周围或建筑物内部淤积大量泥沙、石块的现象。淤埋过程中可能会伴随冲击作用,不过一般不会造成建筑物结构损坏,但会降低建筑物的使用寿命或使其丧失功能性。不仅如此,淤埋危害相对于冲击危害的影响范围更广,危害对象更多,包括农田、房屋以及基础设施等,严重时甚至可能淤塞主河,形成堰塞湖,造成严重的灾害链[6]。如2010年映秀镇红椿沟暴发特大泥石流,冲出的泥石流堵断岷江河道形成堰塞湖,导致洪水冲向映秀新城,致使数千群众受灾[4]。

淤埋危害通常发生在泥石流堆积区,其危害程度与泥石流的淤积深度有关[28]。淤积深度过大不仅造成承灾体功能性丧失,还可能影响其整体结构。泥石流冲出总量越大,淤积范围和深度往往越大。一次泥石流冲出总量可通过实测或计算的方法确定,实测方法精度高,但往往不具备测量条件。实际应用中,往往通过计算暴雨洪峰流量,再按配方法计算泥石流峰值流量,最后根据泥石流暴涨暴落的特点,按泥石流历时和最大流量,将其过程线概化成五边形进行计算。具体计算方法如下。

暴雨洪峰流量Qp计算公式为

(1)

式(1)中:Qp为频率为p时的暴雨峰值流量,m3/s;ψ为流体的径流系数;Sp为p频率下的雨力,mm/h;τ为汇流时间,h;n为暴雨参数;F为汇水面积,km2。

泥石流峰值流量Qc计算公式为

Qc=QpDBF

(2)

式(2)中:Qc为与Qp相同重现期的泥石流峰值流量,m3/s;D为堵塞系数;BF为放大因子。

泥石流一次过程总量计算公式为

Q=0.264TQc

(3)

式(3)中:Q为泥石流一次过程总量,m3;T为泥石流历时,s。

对于一条特定的泥石流沟,不同降雨频率诱发的泥石流冲出总量不同,淤埋危害也不相同,而泥深是衡量泥石流淤积能力最直接的指标。

2 单沟泥石流危险性评价方法

2.1 评价指标

泥石流的危险性来自其危害,而泥石流的危害主要表现为对沿途承灾体的淤埋和冲击破坏,而泥石流淤埋危害与泥深有关,冲击危害与流速有关。此外,承灾体被泥石流淤积的过程中往往还伴随有冲击破坏,因此,基于动力过程的泥石流危险性评价指标应包含泥石流泥深、流速,以及二者的组合。

2.1.1 泥石流泥深

泥深可以反映泥石流对建筑物的淤埋能力,泥深越大对建筑物的淤埋危害越严重。对于一条特定的泥石流沟,泥深主要与一次泥石流过程冲出总量有关。对蒋家沟50组泥石流泥深观测数据分析发现[29],泥石流泥深一般在3.5 m以下,并且在2 m以下居多。不同的泥深对建筑物和人员的危害程度是不同的。唐川等[30]将泥深1～2 m划分为高危险区,泥深0.5～1 m划分为中危险区,泥深小于0.5 m划分为低危险区;常鸣[31]则将泥深大于2.5 m定为高危险,泥深0.5～2.5 m定为中危险区,泥深小于0.5 m定为低危险区;Lin等[32]则认为泥石流前端堆积厚度大于1 m为高危险,小于0.2 m为低危险。

故研究将泥深大于普通住宅楼一半高度即1.5 m时,定为高危险区,该泥深范围内,建筑物内外部淤埋严重,难以恢复使用,同时人员一旦陷入泥石流中,将难以逃脱,会造成严重的人员伤亡和财产损失。将泥深0.5～1.5 m划为中危险,泥石流会部分淤埋建筑物,建筑物在一定程度上可以恢复使用。当泥深小于0.5 m时,对建筑物和人员都难以造成危害,但农田淤埋严重。

2.1.2 泥石流流速

流速是对泥石流冲击力的直接体现,流速越大冲击破坏越剧烈。实际观测数据表明,泥石流的流速通常在4～10 m/s,最小的可能不足1 m/s,较大的则可达13～14 m/s,甚至更高[29,33],如舟曲特大泥石流出山口的瞬时速度可能达到27 m/s。流速是影响冲击力的关键,不同的流速对建筑物的破坏程度不同。Wei等[34]研究发现,框架结构和砌体结构建筑物遭受泥石流灾害的极限承压力分别为110.56 kN/m2和18.22 kN/m2,通过冲击力计算公式反算出两种建筑结构破坏时的临界速度分别为4.3 m/s和1.7 m/s。因此,将流速2 m/s和4 m/s作为泥石流造成破坏时的流速分界点,当流速小于2 m/s时对建筑物的冲击危害较小,不易对建筑物造成破坏;当流速在2～4 m/s时,砖砌体结构的建筑物可能会遭受破坏;当流速超过4 m/s时框架结构建筑可能会被破坏。

2.1.3 泥深和流速的组合

泥石流的冲击和淤埋危害往往同时发生,仅用泥深或流速单个指标不能有效地刻画其破坏能力。因此,很多学者以泥深和流速的组合来反映泥石流的危险性。例如,Fiebiger[35]将泥深和流速进行组合,定义出了泥石流强度公式,用强度大小来决定其危险等级。此外,还有以动量最大值划分泥石流危险性的,Hu等[36]通过研究发现,动能或动量的综合指标比泥深或流速等单个指标更适用于泥石流危险性分区,并推导出用流速和泥深乘积的最大值作为评估动量的参数。

结合上述对评价指标的分析,研究用泥深反映泥石流的淤埋能力,用综合评价指标动量反映泥石流的冲击破坏能力,并以泥深和动量两个指标,划分泥石流危险性等级。

2.2 泥石流危险性分级标准

基于泥深和动量的泥石流危险区划分方法,前人已经做了大量研究,但不同学者对泥石流危险性的划分标准并不一致,表2列出了不同学者的泥石流危险性划分标准。

表2 不同学者的泥石流危险性划分表Table 2 Debris flow hazard classification table of different scholars

考虑到中国山区的主要建筑结构有框架结构、砖混砌体结构和木结构,并多以砌体结构和木结构为主。框架结构具有较好的抵抗地质灾害的能力,而砌体结构和木结构抵抗地质灾害的能力较弱,受到泥石流危害时可能导致房屋损坏或倒塌。分析发现,当泥深达到1.5 m,流速达到4 m/s时,会对建筑物造成严重的淤埋危害,降低建筑物使用寿命或使其丧失功能性;同时在该流速下即使是框架结构的建筑物也可能受到破坏,严重时甚至使建筑物倾斜、坍塌。当泥深在0.5～1.5 m,流速在2～4 m/s时,建筑物局部被淤埋,但农田和基础设施等淤埋严重;砖砌体结构的建筑物可能会被泥石流破坏,生活生产场所恢复有一定困难。当泥深小于0.5 m,流速小于2 m/s时,泥石流对建筑物和人员基本无损伤,但可能造成农田道路被淤埋,危害程度相对较轻。

通过分析,并基于Hu等[36]推导出的动量评估参数,将造成危害的最大泥深和最大流速乘积作为动量分界值即6 m2/s和1 m2/s。于是,提出表3的泥石流危险性划分标准,将泥石流危险性划分为高、中、低3个等级。高危险区主要指泥石流泥深(H≥1.5 m)或动量(HV≥6 m2/s)较大的区域,区内建筑物可能会遭受严重破坏,甚至会造成人员伤亡,危害程度大。低危险区主要是泥石流流速(H<0.5 m)和动量(HV<1 m2/s)都比较低的区域,泥石流基本不会对建筑物造成损坏,但可能会淤埋农田、公路等基础设施,危害程度低。除高危害区和低危险区以外的泥石流危害区域为中等危险区,区内建筑物会有个别损坏,几乎不会出现人员伤亡的情况。

表3 基于动力过程的泥石流危险性划分表Table 3 Hazard classification table of debris flow based on dynamic process

2.3 泥石流动力学指标获取方法

传统的动力学参数估算方法可以得到特定位置的流速和泥深,但无法获得整个泥石流危害区的动力学参数。随着对泥石流动力学机理研究的不断深入和计算机技术的快速发展,数值模拟成为泥石流动力学研究的重要手段。数值模拟不仅可以直观地展示泥石流的运动过程和危害范围,而且可以获得泥深、流速等动力学参数及其分布情况,非常适合基于动力过程的泥石流危险性评价。

随着对泥石流动力学机理认识的不断深入,建立了大量的泥石流运动的本构方程[40-41],为泥石流数值模拟奠定了基础。在较大尺度问题分析中,还可以把数值计算模型或结果与地理信息系统软件相结合[42]。目前,常用于泥石流数值模拟的软件主要有FLO-2D软件[43-44],Debris-2D软件[45]和Massflow软件[46]等。Massflow 软件是从广义的深度积分方法出发,将泥石流运动三维问题简化为二维问题,具有计算效率高、支持二次开发的优点。Debris-2D的控制方程是基于浅水假设的质量和动量守恒方程,在模拟时无需输入曼宁系数和降雨水位图,但需要提供底床物质的屈服应力值、物源的分布等。FLO-2D是基于非牛顿流体模型及有限差分的数值方法,用于计算二维洪水和泥石流运动的模拟软件,它的基底摩擦模型默认为曼宁模型;通过输入曼宁系数、体积浓度、剪切应力、黏度以及流量曲线等参数,可得到整个区域的泥深、平均流速、最大泥深和冲击力等动力学参数,具有模拟结果精度高、计算时间短、使用方便等特点。拟采用FLO-2D对2010年甘肃舟曲“8·7”特大泥石流进行模拟,并利用上述方法对其危险性进行分区。

3 舟曲“8·7”特大泥石流危险性评价

3.1 泥石流概况

甘肃省舟曲县地处青藏高原东北缘,西秦岭山脉与岷山山脉的交汇地区,新构造运动和地震活动强烈,是中国四大泥石流活跃区之一[47]。2010年8月7日23:00左右,舟曲县城东北部发生了短时强降雨,导致罗家峪和三眼峪两沟同时暴发大规模泥石流。泥石流冲毁拦挡坝,穿过舟曲县城直冲入白龙江,形成堰塞湖,导致白龙江水位上升10 m左右,淹没了大半个县城[47]。据统计,泥石流冲毁房屋5 508间,淤埋耕地约93 hm2,导致2万余人受灾,1 481人遇难,284人失踪,1 824人受伤,是中华人民共和国成立以来最严重的泥石流灾害事件[48]。此次泥石流激发雨量为百年一遇的水平。

其中,三眼峪泥石流直接冲击县城,造成的危害最为严重,以三眼峪“8·7”泥石流为例,利用前述方法对其危险性进行评价。三眼峪流域面积24.1 km2,主沟长6.2 km,平均比降241‰。流域内最高海拔为3 772 m,最低点1 278 m(图1),山坡坡度大于35°以上的区域占总面积的51.7%以上,整个流域地形切割强烈,山势陡峻。

图1 研究区位置及流域概况图Fig.1 Location and overview map of the study area

3.2 泥石流动力学参数

FLO-2D数值模拟首先需要对地形数据进行处理,然后确定模拟参数和清水流量。

3.2.1 地形数据处理

模拟选取精度为12.5 m的数字高程模型(digital elevation model,DEM)数据,并将高程数据转换成ASCⅡ文件,导入FLO-2D软件;通过FLO-2D的grid工具对网格进行高程点的插值计算,并将流域划分为10 m×10 m的网格,同时圈定模拟的流域边界,完成地形数据处理。

3.2.2 模拟参数选取

模拟参数主要依据文献[49]并参考FLO-2D使用手册[50]对该流域的体积浓度(Cv)、泥沙相对密度(Gs)、黏滞系数参数(α、β)、层流阻滞系数(K)和曼宁粗糙系数(n)进行综合选取。具体计算参数取值见表4。

表4 泥石流模拟参数Table 4 Debris flow simulation parameter

3.2.3 暴雨峰值流量计算

流量过程曲线计算主要参照《四川省中小流域暴雨洪水计算手册》,研究区位于甘肃省南部,与川北部山区距离较近,地质、地貌、气候均有相似性,计算手册中的参数在研究区的降雨模拟中具有较好的适用性,按式(1)计算出三眼峪泥石流沟百年一遇降雨频率下的暴雨峰值流量,结果见表5。

表5 暴雨峰值流量Table 5 Rainstorm peak flow

3.3 模拟结果

流速模拟结果显示(图2),在整个运动过程中上游沟道的流速较高,主要是因为上游沟道狭窄,地形起伏大,最大流速为11.8 m/s,位于上游沟道地形急剧变化处;泥石流运动到出山口后,由于沟道变缓,地势宽阔,流速降低,但堆积范围明显扩大。当泥石流到达县城时,由于沟道被建筑物过分侵占,导致沟道重新变窄,形成一种束口效应,导致流速再次增加,破坏力增加,给县城造成严重破坏。泥深模拟结果显示,沟道中心泥深较大,基本在4 m以上,沿着沟道中心向两侧逐渐递减,最大泥深位于白龙江内,淤积厚度达10.1 m。

图2 P=1%降雨频率下泥石流模拟结果Fig.2 Debris flow simulation results under P=1% rainfall frequency

3.4 结果验证

为了验证模拟结果的可靠性,将模拟所得的泥石流堆积范围和冲出方量与刘传正等[33]、朱立峰等[51]对舟曲泥石流的调查结果进行比较,并根据泥石流模拟堆积面积Sm,泥石流实际堆积面积Sr与重叠面积S0作为数据基础,计算模拟精度I。I越高,模拟结果越合理,计算公式为

(4)

计算结果见表6。结果表明,泥石流淤积范围模拟的结果比实地调查范围偏小,误差约13%; 冲出方量模拟结果大于实际情况,误差约12%;模拟精度为88%,模拟结果可信度较高。

表6 三眼峪泥石流模拟结果与实际情况对照表Table 6 Comparison of simulated results and realities of Sanyanyu debris flow

3.5 泥石流危险性评价

根据表3的泥石流危险性分区标准,利用ArcGIS软件栅格处理工具,完成泥石流危险性分区图(图3)。危险性分区结果表明,高危险区范围较大,主要分布在沟道中心,面积约284 700 m2,占总面积的49%;中危险区主要分布在泥石流形成区和堆积区边缘,面积约170 400 m2,占总面积的29%;低危险区主要分布在沟道两侧和堆积扇翼端,面积约127 600 m2,占总面积的22%。模拟结果与灾后影像对比发现(图4),靠近沟道中心的月圆村基本处于高危险区域内,只有边缘少部分属于中、低危险区,这也就导致月圆村在强大的泥石流冲击下几乎被夷为平地。其他靠近沟道中心的村子如北关村、三眼村、南街村主要位于中、低危险区,部分位于高危险区域,村寨部分被摧毁,破坏程度相对较小。这与胡凯衡等[47]、赵成等[52]的现场实际调查结果基本相符,因此危险性分区结果可靠性较高。