黄土高原淤地坝溃决型泥石流易发性评价

马建全， 吴钶桥， 李识博*， 郑苗苗， 陈建平， 高帅

(1.西安科技大学地质与环境学院， 西安 710054； 2.青海省青藏高原北部地质过程与矿产资源重点实验室， 西宁 810012；3.矿山地质灾害成灾机理与防控重点实验室， 西安 710054； 4.陕西省地质环境监测总站， 西安 710054)

淤地坝是黄土高原小流域水土流失综合防治中最为简单有效的工程措施[1-3]，具有拦截泥沙、滞洪削峰、淤地成田的作用，不仅增加了农民收入，还促进了生态环境建设与修复[4]。但淤地坝溃决形成的泥石流灾害对农业生产、人民生活及生命财产安全也造成极大危害，如1994年7—8月陕北地区7 347座淤地坝遭到水毁破坏，严重影响当地居民正常生活[5]。导致众多库淤坝毁事件发生的最根本的原因在于突发性强降雨[6]。随着全球气候变化加剧，气候敏感和生态环境脆弱的黄土高原局部地区极端气候频发，为防止水土流失问题加剧，解决淤地坝溃决型泥石流问题迫在眉睫[7]。

在黄土高原淤地坝溃决型泥石流研究方面，众多学者的焦点在淤地坝建设问题、溃决洪水分析、淤地坝拦沙贡献率等，如杨媛媛等[8]提出一种基于地理信息系统(geographic information system，GIS)和土壤侵蚀模数的淤地坝逐年拦沙量计算方法，分析流域内淤地坝拦沙能力；Tung等[9]以动态和静态两种情况研究水文和水力条件的不确定性，并提出了一种防漫顶风险模型。美国陆军工程师团采用相对风险指数对大坝的综合风险管控进行研究[10]。Augutis等[11]将故障树、事件树、森林模型的理论引进大坝安全分析研究中。张振伟等[12]、蔡荨等[13]、胡建平等[14]依据不同理论建立了土石坝综合安全评价指标体系。目前这些研究的关注点在大坝的安全及风险评价、防洪拦沙作用，而对于淤地坝溃决形成泥石流方面的关注度相对较少。当前对黄土高原淤地坝溃决型泥石流的研究与溃决型泥石流的研究并没有区分开来，如Liu等[15]从溃决型泥石流的成因和分布特征方面进行了分析，为溃决型泥石流灾害防控提供了技术支持；王跃等[16]利用遥感技术分析典型冰湖溃决泥石流灾害调查成因机理，建立了冰湖溃决泥石流危险性评价模型，评价结果与实际调查结果一致；曹春然等[17]对不同因素作用下的泥石流堰塞坝溃决过程进行研究分析，为泥石流堰塞坝的灾害防治提供依据；廖立业等[18]通过影像解译和现场调查，分析了泥石流形成的降雨激发条件，改进了泥石流启动与降雨强度之间的关系式；黄勋等[19]利用物理模型和数值模拟，定量分析了强震区泥石流的动力特性；刘春玲等[20]通过遥感分析和实地调查，对喜马拉雅山地区泥石流发育特征进行了研究；彭铭等[21]、石振明等[22]、朱兴华等[23]对堰塞坝的研究进行了总结和展望，提出了今后的研究方向与重点。

综上所述，针对黄土高原淤地坝溃决型泥石流易发性方面研究内容相对较少，大多是黄土高原淤地坝和溃决型泥石流两个方面的研究。而淤地坝库满溃决时所形成的泥石流与洪水又有所区别，因此，现以黄土高原典型小流域中的淤地坝溃决型泥石流易发性为研究对象，从工程地质条件、产沙量和降水量等入手，建立淤地坝溃决型泥石流易发性评价模型，开展泥石流易发性评价研究，旨在提供一种基于暴雨净雨水量、侵蚀产沙量的泥石流易发性定量评价方法和研究思路。

1 易发性评价模型及其计算方法

黄土高原属于干旱半干旱地区，其地表地形破碎、土质疏松、降雨稀少且植被覆盖度低，水土流失问题严重。在黄土高原建筑淤地坝就是为了拦沙截泥。依据已有文献[24-26]，淤地坝之所以发生溃决就是因为超标准暴雨及径流产生的入库水量和入库沙量太多而导致的，其中以库满水沙漫顶型的模式为主(图1)。降雨产生的入库水量可以用产流模型得出，径流产生的入库沙量可以用产沙模型得出。因此，以黄土高原典型小流域为例，开展黄土高原淤地坝溃决型泥石流易发性评价与分析。

图1 淤地坝溃决型泥石流孕灾过程Fig.1 The formation process of debris flow induced by check-dam outburst

1.1 易发性评价建模过程

易发性评价建模过程(图2)主要如下。

(1)以研究区工程地质条件和降水量为基础，采用产流模型[27-29]进行产流计算，采用产沙模型[30]进行产沙计算。

(2)由产流和产沙计算结果得出入库洪量，根据原有水量和出库洪量得出水库总水量。

(3)根据研究区水库总水量，通过库容关系曲线找出对应水位高程，再减去坝底高程得出库水位高差。

(4)以库水位高差和极限库水位高差的比值进行淤地坝溃决型泥石流易发性评价。

1.2 产流计算方法

降雨时产生的淤地坝入库水量主要来自两个方面，一是库区本身水面直接降雨量，二是库区辐射范围内的产流量。产流入库水量计算公式为

(1)

图2 易发性评价分析过程示意图Fig.2 Process diagram of susceptibility assessment and analysis

式(1)中：Rn为产流计算得到的净雨量，mm；Ac为淤地坝控制面积，m2。

流域产流量通常用产流模型计算出，结合实际情况，为方便计算流域因降雨所产生的净雨量，采用初损后损法[31]进行计算，计算公式为

(2)

1.3 产沙计算方法

暴雨发生时，受到侵蚀的大量泥沙随同径流一起进入淤地坝。根据产沙模型，入库沙量Ws计算公式为

(4)

由于因子K、L、S、C、P是淤地坝控制流域固有属性，在短时间内不会发生明显变化。所以，模型认为在每年雨洪季节来临前，这5个因子基本固定。于是令In=KLSCP，In为淤地坝控制流域系数，则入库沙量Ws计算公式为

(5)

因此暴雨发生时，实时求出降雨侵蚀力因子R就能求出淤地坝控制流域实时入库沙量Ws。

1.4 库水位高差的计算

库水位高差Δh受淤地坝入库洪量、剩余库容、排水能力的综合影响，计算公式为

Δh=zw-zb

(6)

式(6)中：zw为淤地坝坝前洪水位，m；由水库总库容W查库容曲线得到；zb为坝底高程，m。

水库总库容W计算公式为

W=W0+Win-Wout

(7)

Win=Ww+Ws

(8)

式中：W0为淤地坝初始库容，m3；Wout为淤地坝出库洪量，m3；Win为淤地坝库容的增加物理量。

1.5 易发性等级确定

库水位高差的最大值为空库状态时库区最低点和坝顶最高点高程差，记为Δhmax，将库水位高差Δh与Δhmax的比值定义为库水位坝高比Rrd。以Rrd值的大小来对淤地坝溃决型泥石流易发程度进行评判，如表1所示。

表1 淤地坝溃决型泥石流易发程度评分表

1.5.1 单坝易发性等级

当淤地坝坝系从上游到下游只有一个坝时，即单坝情形。对于单坝，此时淤地坝溃决型泥石流易发性计算只考虑自身控制范围产生的洪水入库洪量Win所带来的危害，直接根据库水位坝高比Rrd值大小来对淤地坝溃决型泥石流易发程度进行评判。

1.5.2 串坝易发性等级

当淤地坝处于坝系结构下游，而上游还存在其他淤地坝时，即串坝情形。在串坝溃决易发性评价时，当上游淤地坝库水位坝高比Rrd≥0.6时即易发性等级达到高易发性及以上时，认为上游淤地坝随时可能会发生溃决，此时，将上游淤地坝水库库容Wu计入下游淤地坝入库洪量。此时下游淤地坝入库洪量Win计算公式为

Win=Win_d+Wu

(9)

式(9)中：Win_d为下游淤地坝原控制流域面积产生的洪水入库总量；Wu为上游淤地坝水库库容。

求出新的入库洪量后，按照式(7)计算水库库容W，然后重新查水位库容曲线，求出库水位坝高比Rrd，再按表1评定相应易发等级。

图3 研究区位置及淤地坝分布图Fig.3 The location of the research area and distribution of check-dam

2 研究区概况

2.1 研究区地理位置

流岷峪沟地处陕西省榆林市米脂县境内，米脂县城以东约12 km处(图3)。流域面积7.39 km2，主沟道长度3.92 km，沟床平均比降2.80%。流岷峪沟流域是典型的黄土丘陵沟壑区地貌，以水力侵蚀和重力侵蚀为主，按侵蚀类型可以划分为沟间地侵蚀和沟谷地侵蚀，一般沟谷地侵蚀比沟间地侵蚀更为剧烈，流域侵蚀产沙主要来源于这一区域。

在流岷峪沟流域从上游到下游共选取3座完整的淤地坝组成淤地坝系(图3)，并以此进行淤地坝溃决型泥石流易发性评价，这3座坝分别是流岷峪沟2号坝、流岷峪沟3号坝、峰塔沟2号坝，此坝系包含了大型坝、中型坝、小型坝，同时具备串联、并联以及混联这三种坝系联合方式，具备一定的代表性。流岷峪沟2号坝上游除了流岷峪沟3号坝和峰塔沟2号坝外，还有1座峰塔沟1号坝，但该坝通过历史遥感影响资料发现坝体有缺口，默认为该坝失去拦蓄能力，故在计算中不予考虑。

2.2 降雨量

根据《榆林市实用水文手册》，设计点暴雨如表2所示。榆林地区初始降雨过程中某时段的降雨量Pa和I0之间相关关系不明显，但初损I0主要分布在20 mm以内，故可认为该地区降雨量小于20 mm的降雨基本不产流，在设计条件下，I0取为20 mm。

表2 设计暴雨表

2.3 淤地坝控制面积

依据图2，按分水岭得到各淤地坝控制面积，如表3所示。

表3 淤地坝控制面积

3 淤地坝溃决型泥石流易发性评价与分析

3.1 淤地坝水位库容曲线

以淤地坝控制范围掩模提取DEM，经分析得到水位库容关系，作图分析得到水位库容关系曲线(图4)。

图4 水位库容曲线Fig.4 Water level and storage capacity curve

3.2 实时入库水量推求

3.2.1 设计暴雨流推求

由表2，根据皮尔逊III(P-III)型频率曲线查表得出设计频率对应的变差系数，则设计历时不同频率对应的设计点暴雨量计算公式为

(10)

计算结果如表4所示。

表4 不同历时各频率设计暴雨成果表

3.2.2 设计暴雨时程分配

依据《榆林市实用水文手册》来进行设计暴雨时程分配，结果如表5所示。

表5 流岷峪沟流域设计暴雨时程分配

3.2.3 实时净雨推求

模拟实时降雨，将设计暴雨按照30 min时段随着时段累加试算法推求设计净雨，结果如表6所示。

表6 逐时段累积净雨

3.2.4 实时入库水量推求

入库水量Ww按照式(1)计算，实时入库水量计算结果如表7所示。

表7 实时入库水量

3.3 入库沙量计算

3.3.1 降雨侵蚀力因子

降雨侵蚀力因子R按照30 min时段根据各频率设计暴雨模拟实时降雨，结合净雨计算结果，当净雨为0时，降雨侵蚀因子R为0，当净雨不为0时，降雨侵蚀力因子R计算结果如表8所示。

表8 实时降雨侵蚀力因子计算结果表

3.3.2 淤地坝控制范围的计算

由于淤地坝控制范围In=KLSCP，其中土壤可蚀因子K采用EPIC[32]侵蚀-生产影响力模型中的计算方法，计算得到K值为0.34 t·hm2·h/(hm2·MJ·mm)，坡长因子L和坡度因子S从数字高程模型(digital elevation model，DEM)中提出(图5、图6)，植被覆盖因子C用归一化植被指数(normalized difference vegetation index，NDVI)数据来表示(图7)，水土保持措施因子P参考相关研究[33]按坡度赋值得出(图8)，最后用栅格计算器计算得到In分布图(图9)，取各淤地坝控制范围In(图10)的平均值，统计结果如表9所示。

3.3.3 入库沙量计算

图5 坡长因子L分布图Fig.5 Distribution of slope length factor L

图6 坡度因子S分布图Fig.6 Distribution of slope factor S

图7 植被覆盖物因子C分布图Fig.7 Distribution of vegetation coverage factor C

图8 水土保持措施因子P分布图Fig.8 Distribution of soil and water conservation factor P

3.4 实时水库总库容计算

不考虑这三座淤地坝出库水量，并且假设淤地坝在汛前初始水量为0，因此水库实时总库容W与入库总量Win相等。入库总量Win按式(8)计算，结果如表11所示。

图9 In分布图Fig.9 In distribution

图10 In示意图Fig.10 Schematic diagram of In

表9 各淤地坝控制范围的平均值计算结果表

表10 实时入库沙量表

3.5 易发性评价结果分析

3.5.1 不考虑坝系影响的淤地坝溃决型泥石流易发性评价

依据表11，查询各淤地坝库容曲线，计算库水位坝高比Rrd，并按表1对各淤地坝进行泥石流易发性评价，评价结果如表12～表14所示。由表12可知，当降水平率为5%、时序条件4.5 h时，流岷峪沟2号坝溃决形成泥石流易发性等级变为高易发，但仍未达到极高易发，在一般降水情况下均为不易发和低易发。由表13可知，当降水频率为20%、时序条件4.5 h时，流岷峪沟3号坝溃决形成泥石流易发性等级变为高易发，当降水频率为1%、时序条件4 h时，达到极高易发，在一般降水情况下均为不易发和低易发。由表14可知，在各频率降水情况下，降雨1.5 h以后，易发性为高易发。

表11 实时入库总量表

表12 流岷峪沟2号坝溃决型泥石流易发性等级表

表13 流岷峪沟3号坝溃决型泥石流易发性等级表

表14 峰塔沟2号坝溃决型泥石流易发性等级表

3.5.2 考虑坝系影响的淤地坝溃决型泥石流易发性评价

当流岷峪沟3号坝和峰塔沟2号坝溃决型泥石流易发性等级为高易发时，其溃决可能性较大，随时可能溃决形成泥石流，此时，当该两个淤地坝的库容并入流岷峪沟2号坝，重新计算库水位高差比，并进行泥石流易发性评价(表15)。由表15可知，考虑坝系影响后，当降水为频率10%、时序条件4 h以上时，流岷峪沟2号坝溃决形成泥石流易发性等级就变为高易发，整体易发性评价等级有所提高。

4 讨论

评价模型及计算方法仅在选取的流域内开展了研究，其适用性和科学性有待探讨。从以下几个方面开展了讨论。

(1)易发性评价模型中，从水文要素(降雨)出发点进行了模型构建和计算，并没有考虑坝体本身工程因素，如坝体渗漏、坝肩渗漏、坝基渗漏、填筑质量、坝体稳定性等，而这些因素也是造成淤地坝溃坝的重要原因，因此，在应用本文模型时，应首先确定淤地坝溃坝的因素。

(2)降雨量对易发性评价结果影响至关重要，但降雨等气象资料一般为区域值，对于小流域范围的降雨量和产流量可能会存在一定程度的误差，一定程度上影响易发性评价结果的精度。

(3)产沙模型需要对多因子进行叠加分析，而各因子之间具有一定相关性，且受区域气候和当地治理工程措施影响较大，具有年际变化规律，如土壤可蚀性因子K、植被覆盖与管理因子C、水土保持措施因子P，三者之间相互影响，因子取值随不同年份、季节可能会发生改变，因此，在易发性评价模型应用过程中应考虑年际变化情况。

表15 流岷峪沟2号坝考虑坝系的溃决型泥石流易发性等级表

5 结论

(1)基于产流模型、产沙模型，建立了淤地坝溃决型泥石流易发性评价模型及其计算方法，定义了库水位坝高比Rrd并开展了考虑坝系结构的淤地坝溃决型泥石流易发性评价。

(2)选取了米脂县流岷峪沟淤地坝溃决型泥石流为例。研究结果表明：小型坝峰塔沟2号坝出现2 h连续降水时就容易发生溃决型泥石流灾害，中型坝流岷峪沟3号坝出现连续3 h以上小于10%概率降水时容易发生溃决型泥石流灾害，大型坝流岷峪沟2号坝在连续降水4.5 h以上时才容易发生溃决型泥石流灾害。最后考虑坝系结构时，位于串坝下游的大型坝流岷峪沟2号坝发生溃决型泥石流灾害的易发性增大。

(3)淤地坝溃决型泥石流易发性受短时集中暴雨的影响比较大，坝系的结构对淤地坝溃决型泥石流易发性也有一定的影响。