CN110485360A - 一种基于虚拟水库溃决过程的山洪流量计算方法 - Google Patents

Abstract

本发明公开了一种基于虚拟水库溃决过程的山洪流量计算方法，基于水动力学、土力学、泥沙动力学的相互耦合关系并结合以往计算中忽略的洼地储水量、土壤含水量等对流量的贡献；该方法通过将坡面每一个水文地貌单元假定为一个虚拟水库，将山洪流量分为坡面径流和虚拟水库溃决泄流，建立了坡面水文过程和泥沙侵蚀过程的综合求解算法，以临界起动流速作为控制条件判别临界状态前后坡面产汇流的状态变化，并通过双曲型冲刷公式计算山洪的坡面侵蚀过程，由此计算泄流量，再叠加坡面水文地貌单元的坡面径流量得到山洪总流量。该计算方法具有明确的物理机制，提高了山洪物理过程研究的科学性和流量计算的准确性，为山洪灾害分析和预警工作提供技术支持。

Description

一种基于虚拟水库溃决过程的山洪流量计算方法

技术领域

本发明涉及山洪灾害分析评价技术领域，尤其涉及一种诱发滑坡、泥石流的山洪流量分析计算方法。

背景技术

近年来，由于全球气候变化，极端水文事件频发，造成全球范围内山丘区山洪灾害的频次和强度日益突出，由于山洪自身的突发性、水量集中且流速大、冲刷能力强等特点，常常裹挟大量泥沙、石块等形成泥石流和滑坡，对国民经济和人民生命财产安全构成了严重的威胁。

目前，山洪灾害分析计算或基于临界雨量法，利用历史统计资料分析山洪的降雨临界条件；或基于动态临界雨量法，利用水文模型和经典水文学理论分析山洪过程和降雨条件；或基于边坡土力学，进行山洪滑坡的稳定性分析。然而，这些方法往往忽略了山洪突发过程中水动力学、土力学、泥沙动力学的相互耦合关系，尤其忽略了洼地储水量、土壤含水量等水量对山洪流量的贡献，因此造成现有的计算方法可靠性低。

通过准确的确定山洪流量，可更加清楚识别坡面水文过程和山洪诱发的机理，可为山洪灾害预警分析和风险管理提供理论依据。

发明内容

本发明的目的是提供一种基于虚拟水库溃决过程的山洪流量计算方法，通过建立水文地貌单元虚拟水库溃决和产汇流之间的关系，从而解决现有技术中存在的问题。

为解决上述技术问题，本发明提供了一种基于虚拟水库溃决过程的山洪流量计算方法，包括：

1、首先通过野外观测和试验获取研究区基础数据，其中所述基础数据具体包括：研究区地理位置、地形地貌、土壤类型数据；研究区内的降雨历时曲线，对应于每个地形点的每一离散时刻均有对应的降雨强度；研究区内土体的物理及力学参数包括初始含水量θ₀、土壤饱和含水量θ_s、土体中值粒径D₅₀、孔隙率n。

2、对于坡面每一个水文地貌单元，可以假定为一个虚拟水库，虚拟水库的上游来水表示降雨；虚拟水库蓄水表示单元内植被截留、洼地储水和土壤蓄水；虚拟水库的排水相当于坡面径流、壤中流和土壤下渗。

3、对于每一个虚拟水库，坡面山洪的是由于前期降雨使土体含水量达到饱和，土体下岩层裂隙中的压力水体的压力剧增，当遇到暴雨时，能量迅速累积，破坏了原有的土体平衡，土体中的压力水体冲破表面覆盖层，瞬间倾泻形成山洪；即相当于虚拟水库由于前期降雨使水库蓄满，当降雨强度超过某阈值时，虚拟水库溃决。此过程类似于水库漫顶溃决机理，即当水库蓄满后，漫顶水流流速超过坝顶土体的临界起动流速，水库开始发生侵蚀溃决。

4、如前所述，对于坡面的每一个水文地貌单元，若虚拟水库达到总库容，且坡面的水流流速超过了土体的临界起动流速，水流开始携沙运动，逐渐产生泥石流或滑坡。即：

v_泄≥v_c (1)

式中v_泄表示坡面水流流速，v_c为泥沙临界起动流速，与土壤类型、土体中值粒径D₅₀、坡面植被有关，可以通过实验测得，具有明确的物理含义。

5、在坡面流速未达到临界起动流速之前，坡面产流流量计算公式：

式中i(t)表示每个水文地貌单元对应的降雨历时曲线，f_t表示t时刻的下渗率，R_s表示坡面产流量。

6、对于每个水文地貌单元代表的虚拟水库，在本次降雨过程中的蓄水量计算公式：

W＝W_d+∫f_tdt (3)

式中W_d表示洼地储水量，W表示虚拟水库总蓄水量。

7、随着降雨的持续，虚拟水库持续蓄水，当坡面流速达到临界起动流速，水流开始冲刷坡面携带泥沙运动，虚拟水库发生破坏，之后地表产流将转为降雨产生的坡面径流和虚拟水库溃决出流的总和，后者即是以往计算时忽略的水量，坡面总产流量计算公式为：

式中w_t为虚拟水库溃决泄流量，q为坡面总产流量。

8、在极短的dt时间内，虚拟水库溃决泄流量与水文地貌单元被冲走的土体体积有关，对于一个较小的水文地貌单元，其受到水流冲刷的剪应力均匀分布，则dt时间内被冲走的土体体积计算公式：

V＝S·dz (5)

式中S表示水文地貌单元的面积，dz表示dt时间内被冲走的土体高度，V表示dt时间内被冲走的土体体积。

9、采用双曲型冲刷公式计算dt时间内被冲走的土体高度，计算公式：

τ₀＝k(τ_t-τ_c) (7)

式中τ₀表示当前剪应力和临界剪应力的差值，a和b为经验系数，与土壤性质有关可以通过实验测得且有明确的物理含义，k为单位转换因子，τ_t表示当前流速对应的剪应力，τ_c表示临界起动流速对应的剪应力。

10、剪应力的计算公式：

τ＝γ_wRJ (8)

式中τ为剪应力，γ_w为水的重度，R为水力半径，J为水力坡度。

11、水力坡度J的计算公式：

式中h_f为沿程的水头损失，l为水流段长度。

12、沿程水头损失h_f的计算公式：

式中d为当量直径，在圆管时取直径，非圆形截面取4倍的水力半径R，此时取4倍的水力半径R，v为断面的平均流速，λ为沿程阻力系数。

13、将式(9)、式(10)代入式(8)可以得到剪应力的计算公式：

式中ρ_w表示水的密度。

14、对于阻力系数λ的确定一般根据莫迪图，即通过莫迪图找出λ、Re、ε/d三个参数的关系。

式中Re为雷诺数；μ为水的动力学粘滞系数，20℃时取1.004×10^-6m²/s；ε为仪器底面的绝对粗糙度，即粗糙部分的平均高度，一般取试样平均粒径的0.9倍，ε/d即为相对粗糙度。

15、如前所述，可以计算水文地貌单元dt时间内被冲走的土体体积，在dt时间内，dt时间内水库溃决泄流量计算公式：

w_t＝V·η (13)

式中η表示该部分土体的蓄水能力，与孔隙率和洼地体积有关。

16、如前所述，根据式(4)(5)(6)(7)(11)(13)可以得到，坡面流速达到起动流速后，dt时间内每个水文地貌单元的山洪流量计算公式：

式中v_t为当前的坡面流速，q_t+1为dt时间水文地貌单元的出流量。在下一个时刻，以q_t+1解出v_t+1作为输入参数迭代计算，即可得到每个水文地貌单元的山洪流量。

17、对于整个研究区，可划分为多个水文地貌单元，按照流向对总流量进行累加即可得到山洪的流量过程。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明所述一种基于虚拟水库溃决过程的山洪流量计算方法流程示意图；

图2是本发明所述水文地貌单元对应的虚拟水库示意图；

图3是本发明所述双曲型冲刷公式示意图；

图4是本发明所述参考莫迪图。

具体实施方式

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

实施例1：该方法包括以下步骤：

步骤1：通过野外观测和试验获取研究区基础数据，其中所述基础数据具体包括：研究区地理位置、地形地貌、土壤类型数据。

步骤2：收集研究区内的降雨历时曲线，对应于每个地形点的每一离散时刻均有对应的降雨强度。

步骤3：收集研究区内土体的物理及力学参数包括初始含水量θ₀、土壤饱和含水量θ_s、土体中值粒径D₅₀、孔隙率n。

步骤4：对于坡面的每一个水文地貌单元，若水文地貌单元蓄满，且坡面的水流流速超过了土体的临界起动流速，水流开始携沙运动，逐渐产生泥石流或滑坡。即：

v_泄≥v_c (1)式中v_泄表示坡面水流流速，v_c为泥沙临界起动流速，与土壤类型、土体中值粒径D₅₀、坡面植被有关，可以通过实验测得，具有明确的物理含义。

步骤5：在坡面流速未达到临界起动流速之前，坡面产流流量计算公式：

式中i(t)表示每个水文地貌单元对应的降雨历时曲线，f_t表示t时刻的下渗率，R_s表示坡面产流量。

步骤6：对于每个水文地貌单元代表的虚拟水库，在本次降雨过程中的蓄水量计算公式：

W＝W_d+∫f_tdt (3)

式中W_d表示洼地储水量，W表示虚拟水库总蓄水量。

步骤7：随着降雨的持续，虚拟水库持续蓄水，当坡面流速达到临界起动流速，水流开始冲刷坡面携带泥沙运动，虚拟水库发生破坏，之后地表产流将转为降雨产生的坡面径流和虚拟水库溃决出流的总和，后者即是以往计算时忽略的水量，坡面总产流量计算公式为：

式中w_t为虚拟水库溃决泄流量，q为坡面总产流量。

步骤8：在极短的dt时间内，虚拟水库溃决泄流量与水文地貌单元被冲走的土体体积有关，对于一个较小的水文地貌单元，其受到水流冲刷的剪应力均匀分布，则dt时间内被冲走的土体体积计算公式：

V＝S·dz (5)

式中S表示水文地貌单元的面积，dz表示dt时间内被冲走的土体高度，V表示dt时间内被冲走的土体体积。

步骤9：采用双曲型冲刷公式计算dt时间内被冲走的土体高度，计算公式：

τ₀＝k(τ_t-τ_c) (7)

式中τ₀表示当前剪应力和临界剪应力的差值，a和b为经验系数，与土壤性质有关可以通过实验测得且有明确的物理含义，k为单位转换因子，τ_t表示当前流速对应的剪应力，τ_c表示临界起动流速对应的剪应力。

步骤10：计算剪应力：

τ＝γ_wRJ (8)

式中τ为剪应力，γ_w为水的重度，R为水力半径，J为水力坡度。

步骤11：计算水力坡度J：

式中h_f为沿程的水头损失，l为水流段长度。

步骤12：计算沿程水头损失h_f：

式中d为当量直径，在圆管时取直径，非圆形截面取4倍的水力半径R，此时取4倍的水力半径R，v为断面的平均流速，λ为沿程阻力系数。

步骤13：将式(9)、式(10)代入式(8)可以得到剪应力的计算公式：

式中ρ_w表示水的密度。

对于阻力系数λ的确定一般根据莫迪图，即通过莫迪图找出λ、Re、ε/d三个参数的关系。

式中Re为雷诺数；μ为水的动力学粘滞系数，20℃时取1.004×10^-6m²/s；ε为仪器底面的绝对粗糙度，即粗糙部分的平均高度，一般取试样平均粒径的0.9倍，ε/d即为相对粗糙度。

步骤14：如前所述，可以计算水文地貌单元dt时间内被冲走的土体体积，在dt时间内，dt时间内水库溃决泄流量计算公式：

w_t＝V·η (13)

式中η表示该部分土体的蓄水能力，与孔隙率和洼地体积有关。

步骤15：根据式(4)(5)(6)(7)(11)(13)可以得到，坡面流速达到起动流速后，dt时间内每个水文地貌单元的山洪流量计算公式：

式中v_t为当前的坡面流速，q_t+1为dt时间水文地貌单元的出流量。在下一个时刻，以q_t+1解出v_t+1作为输入参数迭代计算，即可得到每个水文地貌单元的山洪流量。

步骤16：对于整个研究区，划分为多个水文地貌单元，按照流向对总流量进行累加即可得到山洪的流量过程。

实施例1：

参见图1，对于某研究区，首先通过野外观测和试验获取研究区基础数据，包括：研究区地理位置、地形地貌、土壤类型数据、降雨历时曲线、土体的物理参数、初始含水量θ₀、土壤饱和含水量θ_s、土体中值粒径D₅₀、孔隙率n。

参见图2，将整个研究区坡面划分为多个水文地貌单元，每个单元对应一个虚拟水库。利用本发明中的计算方法计算坡面产汇流，并根据研究区土体的起动流速判别临界状态。

参见图3，根据双曲型冲刷公式计算每个水文地貌单元的侵蚀过程，并根据本发明提供的计算方法计算虚拟水库泄流量。

将每个水文地貌单元的计算结果按照流向进行叠加计算，即可得到该研究区的山洪流量。

本发明实施例1中所提到的计算方法可靠性高且方便快捷，对山洪灾害的预警分析提供技术支撑。

Claims (6)

1.一种基于虚拟水库溃决过程的山洪流量计算方法，其特征在于，包括

首先通过野外观测和试验获取研究区基础数据，通过将坡面每一个水文地貌单元假定为一个虚拟水库，将山洪流量分为坡面径流和虚拟水库溃决泄流两部分，考虑了以往计算中忽略的洼地储水量、土壤含水量等水量对山洪流量的贡献，建立了坡面水文过程和泥沙侵蚀过程的综合求解算法，以临界起动流速作为控制条件判别临界状态前后坡面产汇流的状态变化，并通过双曲型冲刷公式计算山洪的坡面侵蚀过程，以此为基础计算虚拟水库的泄流量，再叠加坡面水文地貌单元的坡面径流量得到山洪总流量。

2.根据权利要求1所述的基于虚拟水库溃决过程的山洪流量计算方法，基础数据具体包括：研究区地理位置、地形地貌、土壤类型数据；研究区内的降雨历时曲线，对应于每个地形点的每一离散时刻均有对应的降雨强度；研究区内土体的物理及力学参数包括初始含水量θ₀、土壤饱和含水量θ_s、土体中值粒径D₅₀、孔隙率n。

3.根据权利要求1所述的基于虚拟水库溃决过程的山洪流量计算方法，其特征在于，所述坡面每一个水文地貌单元设为一个虚拟水库，虚拟水库的上游来水表示降雨；虚拟水库蓄水表示单元内植被截留、洼地储水和土壤蓄水；虚拟水库的排水相当于坡面径流、壤中流和土壤下渗。

4.根据权利要求1所述的基于虚拟水库溃决过程的山洪流量计算方法，其特征在于，对于每一个虚拟水库，坡面山洪的是由于前期降雨使土体含水量达到饱和，土体下岩层裂隙中的压力水体的压力剧增，当遇到暴雨时，能量迅速累积，破坏了原有的土体平衡，土体中的压力水体冲破表面覆盖层，瞬间倾泻形成山洪；即相当于虚拟水库由于前期降雨使水库蓄满，当降雨强度超过某阈值时，虚拟水库溃决，此过程类似于水库漫顶溃决机理，即当水库蓄满后，漫顶水流流速超过坝顶土体的临界起动流速，水库开始发生侵蚀溃决。

5.根据权利要求1所述的基于虚拟水库溃决过程的山洪流量计算方法，其特征在于，所述的坡面每一个水文地貌单元，若虚拟水库达到总库容，且坡面的水流流速超过了土体的临界起动流速，水流开始携沙运动，逐渐产生泥石流或滑坡，即：

v_泄≥v_c (1)

式中v_泄表示坡面水流流速，v_c为泥沙临界起动流速，与土壤类型、土体中值粒径D₅₀、坡面植被有关，可以通过实验测得，具有明确的物理含义。

6.根据权利要求1所述的基于虚拟水库溃决的山洪流量计算方法，其特征在于，所述的计算方法包括：

在坡面流速未达到临界起动流速之前，坡面产流流量计算公式：

式中i(t)表示每个水文地貌单元对应的降雨历时曲线，f_t表示t时刻的下渗率，R_s表示坡面产流量；

对于每个水文地貌单元代表的虚拟水库，在本次降雨过程中的蓄水量计算公式：

W＝W_d+∫f_tdt (3)

式中W_d表示洼地储水量，W表示虚拟水库总蓄水量；

随着降雨的持续，虚拟水库持续蓄水，当坡面流速达到临界起动流速，水流开始冲刷坡面携带泥沙运动，虚拟水库发生破坏，之后地表产流将转为降雨产生的坡面径流和虚拟水库溃决出流的总和，后者即是以往计算时忽略的水量，坡面总产流量计算公式为：

式中w_t为虚拟水库溃决泄流量，q为坡面总产流量；

在极短的dt时间内，虚拟水库溃决泄流量与水文地貌单元被冲走的土体体积有关，对于一个较小的水文地貌单元，其受到水流冲刷的剪应力均匀分布，则dt时间内被冲走的土体体积计算公式：

V＝S·dz (5)

式中S表示水文地貌单元的面积，dz表示dt时间内被冲走的土体高度，V表示dt时间内被冲走的土体体积；

采用双曲型冲刷公式计算dt时间内被冲走的土体高度，计算公式：

τ₀＝k(τ_t-τ_c) (7)

式中τ₀表示当前剪应力和临界剪应力的差值，a和b为经验系数，与土壤性质有关可以通过实验测得且有明确的物理含义，k为单位转换因子，τ_t表示当前流速对应的剪应力，τ_c表示临界起动流速对应的剪应力；

剪应力的计算公式：

τ＝γ_wRJ (8)

式中τ为剪应力，γ_w为水的重度，R为水力半径，J为水力坡度；

水力坡度J的计算公式：

式中h_f为沿程的水头损失，l为水流段长度；

沿程水头损失h_f的计算公式：

式中d为当量直径，在圆管时取直径，非圆形截面取4倍的水力半径R，此时取4倍的水力半径R，v为断面的平均流速，λ为沿程阻力系数；

将式(9)、式(10)代入式(8)可以得到剪应力的计算公式：

式中ρ_w表示水的密度；

对于阻力系数λ的确定一般根据莫迪图，即通过莫迪图找出λ、Re、ε/d三个参数的关系；

式中Re为雷诺数；μ为水的动力学粘滞系数，20℃时取1.004×10^-6m²/s；ε为仪器底面的绝对粗糙度，即粗糙部分的平均高度，一般取试样平均粒径的0.9倍，ε/d即为相对粗糙度；

可以计算水文地貌单元dt时间内被冲走的土体体积，在dt时间内，dt时间内水库溃决泄流量计算公式：

w_t＝V·η (13)

式中η表示该部分土体的蓄水能力，与孔隙率和洼地体积有关；

根据式(4)(5)(6)(7)(11)(13)可以得到，坡面流速达到起动流速后，dt时间内每个水文地貌单元的山洪流量计算公式：

式中v_t为当前的坡面流速，q_t+1为dt时间水文地貌单元的出流量，在下一个时刻，以q_t+1解出v_t+1作为输入参数迭代计算，即可得到每个水文地貌单元的山洪流量；

对于整个研究区，可划分为多个水文地貌单元，按照流向对总流量进行累加即可得到山洪的流量过程。