CN110717231B - 一种基于坡面沟道河道三级结构的子流域汇流模拟方法 - Google Patents

Abstract

本发明公开了一种基于坡面沟道河道三级结构的子流域汇流模拟方法，涉及分布式水文模型汇流模拟技术领域；包括如下步骤：基于DEM数据提取河道和沟道，并将流域栅格按类别分成坡面栅格、沟道栅格以及河道栅格；计算每个栅格到流域出口的汇流路径长度，并根据沟道栅格类别，统计沟道数量以及沟道平均长度；根据每个栅格最终汇入的河道或者沟道栅格，设置流域内各栅格的汇流属性；根据栅格汇流属性，分别确定等高带坡面、沟道、河道径流分配系数，采用运动波方程作为坡面、沟道、河道汇流的基本模拟公式，对子流域和等高带进行汇流模拟。该方法能够考虑沟道汇流过程对水循环的影响，提高汇流过程模拟效果。

Description

一种基于坡面沟道河道三级结构的子流域汇流模拟方法

技术领域

本发明涉及分布式水文模型汇流模拟技术领域，尤其涉及一种基于坡面沟道河道三级结构的子流域汇流模拟方法。

背景技术

分布式水文模型是研究水循环过程和机理的有效手段，也是解决许多水文实际问题的有效工具。采用子流域对研究流域进行离散，反映参数的空间异质性，是比较常用的方法。该方法具有保持对流域上下游关系刻画精度的同时，大幅度减少计算单元个数的优点。然而受限于子流域划分方法的制约，即要求子流域内有且只有1条河道，使得子流域数量和子流域内河网密度相冲突。当采用较大密度的模拟河网作为子流域划分依据时，不可避免地会生成数量众多的子流域，尤其在超大面积研究流域；当划分较少子流域的时候，又不可避免地减小了河网密度，使得原先快速的河道汇流变成了慢速的坡面汇流。另一方面，对于黄土高原而言，沟壑纵横，这些沟壑在坡面产汇流过程中起着重要的快速汇流作用。而通常采用的坡面河道汇流结构并不能很好地对这种情况进行模拟。为了进一步提高流域坡面河道汇流过程模拟精度，有必要在模型汇流模拟过程中增加对沟壑汇流过程的模拟。

发明内容

本发明的目的在于提供一种基于坡面沟道河道三级结构的子流域汇流模拟方法，从而解决现有技术中存在的前述问题。

为了实现上述目的，本发明采用的技术方案如下：

一种基于坡面沟道河道三级结构的子流域汇流模拟方法，包括如下步骤：

S1，采用子流域套等高带的方式作为基本模拟计算单元，基于DEM对子流域和等高带内的栅格进行分析统计，获取子流域、等高带相关参数；

S2，对每个子流域，采用坡面-沟道-河道三级结构作为子流域汇流的汇流结构；

基于DEM数据，提取河道和沟道，并将流域栅格按类别分成3大类：坡面栅格、沟道栅格以及河道栅格；

S3，计算每个栅格到流域出口的汇流路径长度，并根据沟道栅格类别，统计沟道数量以及沟道平均长度；

S4，根据每个栅格最终汇入的河道或者沟道栅格，设置流域内各栅格的汇流属性；包括最低等高带沟道汇流栅格、最低等高带河道汇流栅格、其他等高带沟道汇流栅格、以及其他等高带河道汇流栅格；

S5，根据步骤S4确定的栅格汇流属性，分别确定等高带坡面、沟道、河道径流分配系数，所述径流分配系数用于匹分坡面汇流量，得到进入下一等高带坡面、沟道以及河道部分的径流量；

S6，采用运动波方程作为坡面、沟道、河道汇流的基本模拟公式，对子流域和等高带进行汇流模拟。

优选地，步骤S1前还包括：提取模拟河道过程，其具体步骤如下：

基于DEM数据，采用坡面汇流累积方法提取模拟河道，包括填洼、流向计算、汇流累计数计算、河道阈值确定及模拟河道提取步骤；

步骤S1中子流域的获取方式根据提取的模拟河道将流域划分成数目不等的子流域，其中，要求每个子流域中只有1条无分岔的河段，并将所有汇入该河段的栅格划分为同一个子流域。

优选地，子流域范围确定后，可通过两种方法确定等高带的个数：

方法一：给定最大等高带个数和最小等高带面积阈值参数计算获得；首先采用子流域面积除以最小等高带面积阈值确定等高带个数，再同最大等高带个数求最小值得到最终的等高带个数；

方法二：给定等高带平均面积参数计算获得，采用子流域面积除以等高带平均面积，向上取整得出最终的等高带个数；

根据等高带个数确定等高带范围，主要包括对子流域内所有栅格按高程降序排列，获得最高高程值和最低高程值，采用等高带个数将这个高程区间进行等分，对等分后的高程区间从高到底依次从1开始自然数编号；

遍历子流域内栅格，根据栅格高程值位于的高程区间分别赋值对应的区间编号，从而确定各等高带范围，即相同编号的栅格位于同一个等高带。需要注意的是，一般情况下DEM高程值为整数，如果高程区间内的整数个数并不足以划分出对应个数的等高带，则需要根据高程区间内的整数个数对等高带个数进行修正。

优选地，根据划分的等高带范围，以DEM为基础数据提取等高带的坡度、长和宽参数，其中，等高带坡度是等高带范围内所有栅格坡度的平均值，等高带宽是沿着坡面流向的长度，等高带长是垂直坡面流向的长度，长和宽计算公式如下：

wid_c＝(H_max-H_min)/sin(slope)

len_c＝A_c/wid

其中，wid_c为等高带宽，len_c为等高带长，H_max是等高带范围内的最大高程，H_min是等高带范围内的最小高程，slope是等高带平均坡度，A_c是等高带面积。

优选地，步骤S2中：

沟道提取步骤为：基于DEM数据，采用坡面汇流累积方法提取模拟沟道，包括填洼、流向计算、汇流累计数计算、沟道阈值确定及模拟沟道提取等步骤；需要注意的是，沟道的提取和河道的提取在方法上是完全一样的，区别在于沟道阈值要小于河道阈值。

栅格按类别划分的标准为：通过河道阈值提取出的模拟河网栅格为河道栅格，通过沟道阈值提取出的非河道栅格的栅格为沟道栅格，剩余的栅格为坡面栅格。

优选地，步骤S3中：

采用GIS软件计算每个栅格到流域出口的汇流路径长度为；所述汇流路径长度为该栅格与流域出口的沿程汇流长度；

由于提取的沟道存在着长度过小的情况，认为这些短的沟道属于提取出的伪沟道，需要加以剔除。因此直接由阈值提取出的沟道为备选沟道。

遍历子流域内的河道栅格，统计流入河道栅格的备选沟道栅格，作为备选沟道的出口。对每个备选沟道出口栅格，统计所有流入该栅格的上游备选沟道栅格，并确定汇流路径长度最大值。采用该最大汇流路径长度减去备选沟道出口栅格的汇流路径长度，得到该备选沟道的概化长度。

设定一个沟道最小长度阈值，并认为长度小于该阈值的备选沟道为提取的伪沟道，长度大于等于该阈值的备选沟道为最终使用的沟道，进而确定该子流域内沟道条数N_g，并对所有沟道长度求平均值得到沟道的平均概化长度len_g。

对流域内所有伪沟道出口进行溯源遍历，将流入该出口栅格的上游所有备选沟道栅格类别属性设置为坡面栅格，消除提取的伪沟道的影响。

优选地，步骤S4中：

实施最低等高带指的是子流域内编号最大的那个等高带，如果子流域只有1个等高带则该等高带就是最低等高带，没有其他等高带；

栅格汇流属性的确定方法为对流域内每个河道和沟道栅格进行溯源遍历，如果沟道栅格位于最低等高带，则给所有汇入该沟道栅格的栅格赋值最低等高带沟道汇流栅格属性(F1)，同理，对所有汇入最低等高带河道栅格的栅格赋值最低等高带河道汇流栅格属性(F2)，对所有汇入其他等高带沟道栅格的栅格赋值其他等高带沟道汇流栅格属性(F3)，对所有汇入其他等高带河道栅格的栅格赋值其他等高带河道汇流栅格属性(F4)。

优选地，步骤S5中：

等高带坡面径流分配系数、沟道径流分配系数、河道径流分配系数分为最低等高带和其他等高带，且最低等高带和其他等高带径流分配系数分别计算；

其中最低等高带径流分配系数计算公式如下：

r_p＝0

r_g＝N_F1/(N_F1+N_F2)

r_r＝N_F2/(N_F1+N_F2)

其他等高带径流分配系数计算公式如下：

r_p＝(N_F1+N_F2)/(N_F1+N_F2+N_F3+N_F4)

r_g＝N_F3/(N_F1+N_F2+N_F3+N_F4)

r_r＝N_F4/(N_F1+N_F2+N_F3+N_F4)

其中，r_p，r_g，r_r分别表示坡面、沟道、河道的径流分配系数；N_F1为该等高带内最低等高带沟道汇流栅格个数；N_F2为该等高带内最低等高带河道汇流栅格个数；N_F3为该等高带内其他等高带沟道汇流栅格个数；N_F4为该等高带内其他等高带河道汇流栅格个数。

优选地，步骤S6中：

采用运动波方法模拟等高带的汇流过程：

ql_i＝q_i/len_i；

S_fi＝S_0i；

其中,ql_i为汇流单宽流入量，q_i为参与汇流过程侧向径流量，len_i为汇流长度，A_i为汇流过水断面面积，Q_i为汇流断面出口流量，n_i为曼宁糙率系数，R_i为水力半径，S_0i为坡降，S_fi为摩擦坡降，t为时间，x为沿汇流方向的坐标，i表示坡面(＝pj,j为等高带编号)、沟道(＝g)、河道(＝r)不同汇流情况。

优选地，对子流域汇流方式进行模拟的过程如下：

对于坡面汇流而言，每个等高带都有坡面汇流过程。汇流断面概化为矩形断面，汇流长度等于等高带宽，汇流断面宽等于等高带长,各等高带坡面汇流上游入流量等于Q_p(j-1)*r_p(j-1)，参与汇流的侧向径流量q_pj为当前等高带产流量，坡降S_0p为等高带坡降。

对于沟道汇流而言，只在子流域层面进行汇流演算，且认为子流域内有N_g条长度为len_g的沟壑，对每条沟道分别进行汇流演算得到各沟道汇流量Q_g，再求和得到总的沟道汇流量Q_g*N_g；沟道的断面概化为倒等腰梯形，相关参数根据调查或文献获得；沟壑汇流上游入流量＝0，参与侧向径流的q_g等于各等高带计算的坡面汇流量Q_pj*r_g之和，坡降S_0g为外部给定的沟道坡降参数。

对于河道而言，也只在子流域层面进行汇流演算，只有1个河段，长度为根据DEM提取的len_r,断面也概化为倒等腰梯形，相关参数根据河流上水文实测断面进行概化得出；河道上游入流量等于上游子流域出口径流量，参与侧向径流的q_r为各等高带计算的坡面汇流量Q_pj*r_r以及沟道汇流量Q_g*N_g总和。

本发明的有益效果是：

本发明公开了一种基于坡面沟道河道三级结构的子流域汇流模拟方法，该方法能够考虑沟道汇流过程对水循环的影响，提高汇流过程模拟效果。

附图说明

图1是本发明应用实例计算流程图；

图2是本发明应用实例等高带划分及栅格分类数值图；

图3是本发明应用实例等高带参数信息；

图4是本发明应用实例沟道选择图；

图5是本发明应用实例汇流属性图；

图6是本发明应用实例提取的子流域坡面、沟道、河道径流分配系数图；

图7是本发明应用实例汇流模拟过程示意图A；

图8是本发明应用实例汇流模拟过程示意图B。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不用于限定本发明。

如图1所示，本发明实施例提供了一种基于坡面沟道河道三级结构的子流域汇流模拟方法，包括如下步骤：

S1，采用子流域套等高带作为模型基本计算单元。本实例以1个子流域为例进行说明，如图2所示。其中，假定该子流域划分为3个等高带，并按栅格高程排序分割得到，等高带1高程最高，等高带3高程最低。基于栅格大小90m的DEM数据计算各栅格坡度，并根据等高带范围统计各等高带坡度、长度和宽参数，相关结果如图3所示。

S2，根据水系提取算法，提取沟道和河道栅格，如图2所示。其中，河道是单一的无分岔河段，沟道是汇入河道的多个分岔河段。在应用实例提取中，河道提取阈值取的是5km²，沟道提取阈值取的是0.5km²。子流域内大部分栅格都是坡面栅格，少部分是沟道和河道栅格。

S3，采用GIS软件计算各栅格到流域出口的汇流路径长度(这里没列出)。查找汇入河道的备选沟道出口栅格，本实例中共有10条备选沟道流入河道，根据各栅格汇流路径长度，求10条备选沟道的长度，如图4所示。其中点栅格为对应备选沟道最长的汇流路径。设置200m作为沟道最小长度阈值，可以发现图4中红圈中的备选沟道长度小于该阈值(180m)，作为伪沟道剔除(见图5中沟道栅格点)。因此本应用实例共9条沟道。为了计算机处理方便，这里直接采用沟道栅格个数(155)乘以栅格边长(90m)作为所有沟道长度，再除以沟道数量，得到每条概化沟道的平均长度1550m。

S4，本应用实例中最低等高带指的是等高带3。按如下步骤对各栅格汇流属性进行赋值。首先将所有栅格赋值为F2；其次对等高带3中的沟道栅格赋值F1，并逐个栅格进行溯源遍历，将汇入其中的所有坡面栅格赋值为F1；第三将等高带1和2中的河道栅格赋值F4，并逐栅格进行溯源遍历，将汇入其中的所有坡面栅格赋值为F4；最后将等高带1和2中的沟道栅格赋值F3，并逐栅格进行溯源遍历，将汇入其中的所有坡面栅格赋值为F3，结果如图5所示。

S5，统计各等高带中F1～F4栅格的个数，并根据公式计算坡面、沟道、河道径流分配系数，结果如图6所示。

S6，采用运动波方程以及图6计算所得分配系数对各等高带坡面汇流以及子流域层面的沟道和河道汇流进行演算，按先等高带坡面汇流、再沟道汇流，最后河道汇流的次序依次进行，汇流计算过程示意如图7和图8所示。

Claims (9)

1.一种基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，包括如下步骤：

S1，采用子流域套等高带的方式作为基本模拟计算单元，基于DEM对子流域和等高带内的栅格进行分析统计，获取子流域、等高带相关参数；

S2，对每个子流域，采用坡面-沟道-河道三级结构作为子流域汇流的汇流结构；

基于DEM数据，提取河道和沟道，并将流域栅格按类别分成3大类：坡面栅格、沟道栅格以及河道栅格；

S3，计算每个栅格到流域出口的汇流路径长度，并根据沟道栅格类别，统计沟道数量以及沟道平均长度；

S4，根据每个栅格最终汇入的河道或者沟道栅格，设置流域内各栅格的汇流属性；包括最低等高带沟道汇流栅格、最低等高带河道汇流栅格、其他等高带沟道汇流栅格、以及其他等高带河道汇流栅格；

S5，根据步骤S4确定的栅格汇流属性，分别确定等高带坡面、沟道、河道径流分配系数，所述径流分配系数用于匹分坡面汇流量，得到进入下一等高带坡面、沟道以及河道部分的径流量；

S6，采用运动波方程作为坡面、沟道、河道汇流的基本模拟公式，对子流域和等高带进行汇流模拟；

步骤S5中：

等高带坡面径流分配系数、沟道径流分配系数、河道径流分配系数分为最低等高带和其他等高带，且最低等高带和其他等高带径流分配系数分别计算；

其中最低等高带径流分配系数计算公式如下：

r_p＝0

r_g＝N_F1/(N_F1+N_F2)

r_r＝N_F2/(N_F1+N_F2)

其他等高带径流分配系数计算公式如下：

r_p＝(N_F1+N_F2)/(N_F1+N_F2+N_F3+N_F4)

r_g＝N_F3/(N_F1+N_F2+N_F3+N_F4)

r_r＝N_F4/(N_F1+N_F2+N_F3+N_F4)

其中，r_p，r_g，r_r分别表示坡面、沟道、河道的径流分配系数；N_F1为该等高带内最低等高带沟道汇流栅格个数；N_F2为该等高带内最低等高带河道汇流栅格个数；N_F3为该等高带内其他等高带沟道汇流栅格个数；N_F4为该等高带内其他等高带河道汇流栅格个数。

2.根据权利要求1所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，步骤S1前还包括：提取模拟河道过程，其具体步骤如下：

基于DEM数据，采用坡面汇流累积方法提取模拟河道，包括填洼、流向计算、汇流累计数计算、河道阈值确定及模拟河道提取步骤；

步骤S1中子流域的获取方式根据提取的模拟河道将流域划分成数目不等的子流域，其中，要求每个子流域中只有1条无分岔的河段，并将所有汇入该河段的栅格划分为同一个子流域。

3.根据权利要求2所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，子流域范围确定后，可通过两种方法确定等高带的个数：

方法一：给定最大等高带个数和最小等高带面积阈值参数计算获得；首先采用子流域面积除以最小等高带面积阈值确定等高带个数，再同最大等高带个数求最小值得到最终的等高带个数；

方法二：给定等高带平均面积参数计算获得，采用子流域面积除以等高带平均面积，向上取整得出最终的等高带个数；

根据等高带个数确定等高带范围，主要包括对子流域内所有栅格按高程降序排列，获得最高高程值和最低高程值，采用等高带个数将这个高程区间进行等分，对等分后的高程区间从高到底依次从1开始自然数编号；

遍历子流域内栅格，根据栅格高程值位于的高程区间分别赋值对应的区间编号，从而确定各等高带范围，即相同编号的栅格位于同一个等高带。

4.根据权利要求3所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，根据划分的等高带范围，以DEM为基础数据提取等高带的坡度、长和宽参数，其中，等高带坡度是等高带范围内所有栅格坡度的平均值，等高带宽是沿着坡面流向的长度，等高带长是垂直坡面流向的长度，长和宽计算公式如下：

wid_c＝(H_max-H_min)/sin(slope)

len_c＝A_c/wid_c

其中，wid_c为等高带宽，len_c为等高带长，H_max是等高带范围内的最大高程，H_min是等高带范围内的最小高程，slope是等高带平均坡度，A_c是等高带面积。

5.根据权利要求1所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，步骤S2中：

沟道提取步骤为：基于DEM数据，采用坡面汇流累积方法提取模拟沟道，包括填洼、流向计算、汇流累计数计算、沟道阈值确定及模拟沟道提取步骤；需要注意的是，沟道的提取和河道的提取在方法上是完全一样的，区别在于沟道阈值要小于河道阈值；

栅格按类别划分的标准为：通过河道阈值提取出的模拟河网栅格为河道栅格，通过沟道阈值提取出的非河道栅格的栅格为沟道栅格，剩余的栅格为坡面栅格。

6.根据权利要求1所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，步骤S3中：

采用GIS软件计算每个栅格到流域出口的汇流路径长度为：所述汇流路径长度为该栅格与流域出口的沿程汇流长度；

由于提取的沟道存在着长度过小的情况，认为这些短的沟道属于提取出的伪沟道，需要加以剔除； 因此直接由阈值提取出的沟道为备选沟道；

遍历子流域内的河道栅格，统计流入河道栅格的备选沟道栅格，作为备选沟道的出口；对每个备选沟道出口栅格，统计所有流入该栅格的上游备选沟道栅格，并确定汇流路径长度最大值；采用该最大汇流路径长度减去备选沟道出口栅格的汇流路径长度，得到该备选沟道的概化长度；

设定一个沟道最小长度阈值，并认为长度小于该阈值的备选沟道为提取的伪沟道，长度大于等于该阈值的备选沟道为最终使用的沟道，进而确定该子流域内沟道条数N_g，并对所有沟道长度求平均值得到沟道的平均概化长度len_g；

对流域内所有伪沟道出口进行溯源遍历，将流入该出口栅格的上游所有备选沟道栅格类别属性设置为坡面栅格，消除提取的伪沟道的影响。

7.根据权利要求1所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，步骤S4中：

实施最低等高带指的是子流域内编号最大的那个等高带，如果子流域只有1个等高带则该等高带就是最低等高带，没有其他等高带；

栅格汇流属性的确定方法为对流域内每个河道和沟道栅格进行溯源遍历，如果沟道栅格位于最低等高带，则给所有汇入该沟道栅格的栅格赋值最低等高带沟道汇流栅格属性F1，同理，对所有汇入最低等高带河道栅格的栅格赋值最低等高带河道汇流栅格属性F2，对所有汇入其他等高带沟道栅格的栅格赋值其他等高带沟道汇流栅格属性F3，对所有汇入其他等高带河道栅格的栅格赋值其他等高带河道汇流栅格属性F4。

8.根据权利要求1所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，S6中：

采用运动波方法模拟等高带的汇流过程：

坡面：

ql_pj＝q_pj/len_pj；

S_fpj＝S_0pj；

其中,ql_pj为第j个等高带汇流单宽流入量，q_pj为第j个等高带参与汇流过程侧向径流量，len_pj为第j个等高带汇流长度，A_pj为第j个等高带汇流过水断面面积，Q_pj为第j个等高带汇流断面出口流量，n_pj为第j个等高带曼宁糙率系数，R_pj为第j个等高带水力半径，S_0pj为第j个等高带坡降，S_fpj为第j个等高带摩擦坡降，t为时间，x为沿汇流方向的坐标；

沟道：

ql_g＝q_g/len_g；

S_fg＝S_0g；

其中,ql_g为沟道汇流单宽流入量，q_g为沟道参与汇流过程侧向径流量，len_g为沟道汇流长度，A_g为沟道汇流过水断面面积，Q_g为沟道汇流断面出口流量，n_g为沟道曼宁糙率系数，R_g为沟道水力半径，S_0g为沟道坡降，S_fg为沟道摩擦坡降，t为时间，x为沿汇流方向的坐标；

河道：

ql_r＝q_r/len_r；

S_fr＝S_0r；

其中，ql_r为河道汇流单宽流入量，q_r为河道参与汇流过程侧向径流量，len_r为河道汇流长度，A_r为河道汇流过水断面面积，Q_r为河道汇流断面出口流量，n_r为河道曼宁糙率系数，R_r为河道水力半径，S_0r为河道坡降，S_fr为河道摩擦坡降，t为时间，x为沿汇流方向的坐标。

9.根据权利要求8所述的基于坡面沟道河道三级结构的子流域汇流模拟方法，其特征在于，对子流域汇流方式进行模拟的过程如下：

对于坡面汇流而言，每个等高带都有坡面汇流过程；汇流断面概化为矩形断面，汇流长度等于等高带宽，汇流断面宽等于等高带长，各等高带坡面汇流上游入流量等于Q_p(j-1)*r_p(j-1)，参与汇流的侧向径流量q_pj为当前等高带产流量，坡降S_0p为等高带坡降；

对于沟道汇流而言，只在子流域层面进行汇流演算，且认为子流域内有N_g条长度为len_g的沟壑，对每条沟道分别进行汇流演算得到各沟道汇流量Q_g，再求和得到总的沟道汇流量Q_g*N_g；沟道的断面概化为倒等腰梯形，相关参数根据调查或文献获得；沟壑汇流上游入流量＝0，参与侧向径流的q_g等于各等高带计算的坡面汇流量Q_pj*r_g之和，坡降S_0g为外部给定的沟道坡降参数；

对于河道而言，也只在子流域层面进行汇流演算，只有1个河段，长度为根据DEM提取的len_r，断面也概化为倒等腰梯形，相关参数根据河流上水文实测断面进行概化得出；河道上游入流量等于上游子流域出口径流量，参与侧向径流的q_r为各等高带计算的坡面汇流量Q_pj*r_r以及沟道汇流量Q_g*N_g总和。

Images