一种基于三级汇流结构的分布式侵蚀产沙模拟方法与流程

本发明涉及土壤侵蚀、河道输沙模型与计算机数值计算技术领域，尤其涉及一种基于三级汇流结构的分布式侵蚀产沙模拟方法。

背景技术：

在水土保持研究过程中，除了产流过程外，泥沙过程的同步模拟必不可少。然则，国内外水文模型在泥沙方面都采用半经验式的模拟方法，物理过程刻画薄弱，不能满足区域水土管理越来越精细化和精准化的需求。

此外，已有的流域土壤侵蚀模型都只对坡面和河道两个侵蚀过程进行模拟，输沙过程计算时往往忽视了沟道过程的重要性。坡面侵蚀显然不能反映沟道的侵蚀特点，而沟道的比降一般都比河道比降大，河道产输沙也不能包含沟道产输沙过程。尤其，黄土高原沟道高度发育，沟道水流侵蚀量占总水力侵蚀量的比例接近50％。这一特点使沟道水力侵蚀过程在黄土高原土壤侵蚀中占据举足轻重的地位，也对水沙物理过程的模拟提出了更高的要求。

目前，国内主流水文模型之一的大尺度流域水与能量转化过程模型(waterandenergytransferprocessesinlargeriverbasins，简称wep-l)已很好地应用于黄土高原降水-径流研究，但该模型也仅考虑了“坡面-河道”两级汇流系统，并且不支持对流域侵蚀产沙等物理过程的模拟分析，模型不能直接应用于黄土高原水土保持分析与生态治理研究。

技术实现要素：

本发明的目的在于提供一种基于三级汇流结构的分布式侵蚀产沙模拟方法，从而解决现有技术中存在的前述问题。

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

一种基于三级汇流结构的分布式侵蚀产沙模拟方法，包括如下步骤：

s1，地形数据提取：利用gis平台和模拟区dem数据，提取河网和坡度等地形数据；

s2，子流域-等高带单元划分：利用gis平台和s1中的河网提取结果，基于汇水区划分子流域；在每个子流域内，按照一定的高程间隔划分等高带；确定坡面-沟道-河道三级汇流结构；

s3，基础数据库构建：收集模拟区内的水文、气象、土壤、土地利用、水土保持措施、人类活动取用水的模型输入数据；并按照模型需求和s2中的划分结果，对收集的模型输入数据进行空间或时间尺度上的处理，构建基础数据库；

s4，计算坡面产输沙：利用s3所建数据库中各等高带的水文气象要素、土地利用条件、水保措施数据，以等高带为单元，按照上下承接的关系，依次计算各等高带的坡面沟间及沟道侵蚀量和输沙率；

s5，计算沟道产输沙：基于s4中坡面产沙结果，以子流域为计算单元，进行沟道侵蚀、输沙模拟，计算沟道输沙率；

s6，计算河道产输沙：基于s4中坡面、s5中沟道产输沙结果，以子流域为单元，进行河道侵蚀、输沙模拟；

s7，侵蚀产沙量统计：坡面沟间及坡面沟道产沙由等高带汇总到子流域；坡面、沟道、河道各环节的产沙由子流域汇总到模拟区；最终，在模拟区内对所有环节的产沙进行汇总。

优选地，s4包括如下步骤：

s401，以等高带为计算单元，暴雨期以每小时降水为驱动，非暴雨期以日降水为驱动，计算坡面沟间雨滴溅蚀；

s402，以产流计算所得坡面沟间径流为驱动，计算坡面沟间径流侵蚀；

s403，每个等高带坡面沟间所产沙量进入当前等高带的坡面沟道，以坡面沟道内水流挟沙力为上限，判断坡面沟道冲淤情况，当上一等高带来沙和当前等高带来沙这两个沙源总的含沙量小于水流挟沙力时，则当前等高带坡面沟道发生冲刷；否则发生淤积；如此从上到下串联进行，实现产沙在坡面等高带间的传递。

优选地，s401中，按照如下公式计算雨滴溅蚀：

式中，splmod为雨滴溅蚀侵蚀模数(kg/m²/s)；e为降雨侵蚀动能(j/m²)；i为雨强(mm/min)；α为坡面的坡度(°)；a1、b1、c1、a2、b2为经验参数。

优选地，s402中，按照如下公式计算坡面径流侵蚀：

slomod＝a3·(q-b3)·(1+c3·α)

式中，slomod为坡面沟间径流侵蚀模数(kg/m²/s)；q为坡面沟间单宽流量(m²/s)；α为坡面的坡度(°)；a3、b3、c3为经验参数。

优选地，s403中，按照如下公式计算水流挟沙力：

式中，tv为体积比挟沙力(m³/m³)；cv为体积比含沙量(m³/m³)；v为水流流速(m/s)；ρs为土粒密度(kg/m³)；ρm为浑水密度(kg/m³)；g为重力加速度(n/kg)；r为水力半径(m)；ω为浑水中的泥沙沉速(m/s)；h为水深(m)；d50为泥沙的中值粒径(m)，粒径小于该值的泥沙占全部泥沙的50％；κ为浑水的卡门常数；κ0为清水的卡门常数；a4、b4、c4、d4、a5、b5为经验参数。

优选地，s5具体为：以子流域为计算单元，沟道承接来自当前子流域坡面最低等高带的输沙，根据汇流结果进行沟道输沙计算。

优选地，s6具体为：在子流域内，多条沟道中的输沙通过并联方式汇入河道，通过河道输沙模拟到达子流域出口的沙量；不同子流域出口的产输沙又基于上下游关系，通过串并联混合方式最终传递到模拟区出口。

优选地，所述输沙计算方法中，对经典悬移质连续方程进行改进，添加侧向来沙项(qls)：

式中，x为沿河方向到沟/河道入口的距离(m)；cx为对应位置的水流含沙量(kg/m³)；tx为对应位置的水流挟沙力(kg/m³)；α为恢复饱和系数,在一般水力因素条件下，平衡时恢复饱和系数在0.02～1.78之间，平均接近0.5；q为单宽流量(m²/s)，qls为单位长度沟道/河道侧向来沙输沙率(kg/s/m)。

优选地，河段出口的含沙量按照如下公式进行计算：

式中，c和c0分别为沟道或河道出口和入口的水流含沙量(kg/m³)；t和t0分别为沟道或河道出口和入口的水流挟沙力(kg/m³)；cl为沟道/河道单位长度侧向来沙的含沙量(kg/m³)；l为河段长度(m)。

本发明的有益效果是：本发明提供的基于三级汇流结构的分布式侵蚀产沙模拟方法，通过改进wep-l模型，将wep-l原有的“坡面-河道”二级系统改进为“坡面-沟道-河道”三级系统，并基于侵蚀产沙和输移机制增加流域侵蚀与产沙模块，具体的，根据从坡面到流域出口的产输沙过程，按照土壤侵蚀发生的位置将水沙过程划分为坡面、沟道、河道三个环节，其中坡面侵蚀环节分为雨滴溅蚀、坡面径流侵蚀和坡面沟道产输沙过程；以各环节的侵蚀特点为依据，根据产输沙时空格局、承接传递过程等方面的相似性，分别以子流域和坡面等高带为计算单元，实现侵蚀产沙过程的多环节细化模拟，有效的反映了沟道侵蚀严重的特点，真实刻画了其侵蚀过程。

附图说明

图1是本发明提供的基于三级汇流结构的分布式侵蚀产沙模拟方法流程示意图；

图2是具体实施例中涉及到的模拟区示意图；

图3是窟野河流域温家川水文站1968-1992年输沙率模拟结果示意图。

具体实施方式

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

如图1所示，本发明提供了一种基于三级汇流结构的分布式侵蚀产沙模拟方法，包括如下步骤：

s1，地形数据提取：利用gis平台和模拟区dem数据，提取河网和坡度等地形数据；

s2，子流域-等高带单元划分：利用gis平台和s1中的河网提取结果，基于汇水区划分子流域；在每个子流域内，按照一定的高程间隔划分等高带；确定坡面-沟道-河道三级汇流结构；

s3，基础数据库构建：收集模拟区内的水文、气象、土壤、土地利用、水土保持措施、人类活动取用水的模型输入数据；并按照模型需求和s2中的划分结果，对收集的模型输入数据进行空间或时间尺度上的处理，构建基础数据库；

其中，对收集的模型输入数据进行空间或时间尺度上的处理，具体可以为，利用泰森多边形法和降水等值线图，将收集到的降水数据空间展布到等高带；考虑水利部水土保持小流域治理工作等，对数据缺失年份的土地利用数据进行插补。

s4，计算坡面产输沙：利用s3所建数据库中各等高带的水文气象要素、土地利用条件、水保措施数据，以等高带为单元，按照上下承接的关系，依次计算各等高带的坡面沟间及沟道侵蚀量和输沙率；具体的，包括如下步骤：

s401，以等高带为计算单元，暴雨期以每小时降水为驱动，非暴雨期以日降水为驱动，计算坡面沟间雨滴溅蚀；具体的，按照如下公式计算雨滴溅蚀：

式中，splmod为雨滴溅蚀侵蚀模数(kg/m²/s)；e为降雨侵蚀动能(j/m²)；i为雨强(mm/min)；α为坡面的坡度(°)；a1、b1、c1、a2、b2为经验参数。

s402，以产流计算所得坡面沟间径流为驱动，计算坡面沟间径流侵蚀；具体的，按照如下公式计算坡面径流侵蚀：

slomod＝a3·(q-b3)·(1+c3·α)

式中，slomod为坡面沟间径流侵蚀模数(kg/m²/s)；q为坡面沟间单宽流量(m²/s)；α为坡面的坡度(°)；a3、b3、c3为经验参数。

s403，每个等高带坡面沟间所产沙量进入当前等高带的坡面沟道，以坡面沟道内水流挟沙力为上限，判断坡面沟道冲淤情况，当上一等高带来沙和当前等高带来沙这两个沙源总的含沙量小于水流挟沙力时，则当前等高带坡面沟道发生冲刷；否则发生淤积；如此从上到下串联进行，实现产沙在坡面等高带间的传递，其中，按照如下公式计算水流挟沙力：

式中，tv为体积比挟沙力(m³/m³)；cv为体积比含沙量(m³/m³)；v为水流流速(m/s)；ρs为土粒密度(kg/m³)；ρm为浑水密度(kg/m³)；g为重力加速度(n/kg)；r为水力半径(m)；ω为浑水中的泥沙沉速(m/s)；h为水深(m)；d50为泥沙的中值粒径(m)，粒径小于该值的泥沙占全部泥沙的50％；κ为浑水的卡门常数；κ0为清水的卡门常数；a4、b4、c4、d4、a5、b5为经验参数。

s5，计算沟道产输沙：基于s4中坡面产沙结果，以子流域为计算单元，进行沟道侵蚀、输沙模拟，计算沟道输沙率；具体为：以子流域为计算单元，沟道承接来自当前子流域坡面最低等高带的输沙，根据汇流结果进行沟道输沙计算。

s6，计算河道产输沙：基于s4中坡面、s5中沟道产输沙结果，以子流域为单元，进行河道侵蚀、输沙模拟；具体为：在子流域内，多条沟道中的输沙通过并联方式汇入河道，通过河道输沙模拟到达子流域出口的沙量；不同子流域出口的产输沙又基于上下游关系，通过串并联混合方式最终传递到模拟区出口。

所述输沙计算方法中，对经典悬移质连续方程进行改进，添加侧向来沙项(qls)：

式中，x为沿河方向到沟/河道入口的距离(m)；cx为对应位置的水流含沙量(kg/m³)；tx为对应位置的水流挟沙力(kg/m³)；α为恢复饱和系数,在一般水力因素条件下，平衡时恢复饱和系数在0.02～1.78之间，平均接近0.5；q为单宽流量(m²/s)，qls为单位长度沟道/河道侧向来沙输沙率(kg/s/m)。

河段出口的含沙量按照如下公式进行计算：

式中，c和c0分别为沟道或河道出口和入口的水流含沙量(kg/m³)；t和t0分别为沟道或河道出口和入口的水流挟沙力(kg/m³)；cl为沟道/河道单位长度侧向来沙的含沙量(kg/m³)；l为河段长度(m)。

s7，侵蚀产沙量统计：坡面沟间及坡面沟道产沙由等高带汇总到子流域；坡面、沟道、河道各环节的产沙由子流域汇总到模拟区；最终，在模拟区内对所有环节的产沙进行汇总。

具体实施例

以黄河一级支流窟野河流域内的温家川水文站控制区1968-1992年侵蚀产沙过程模拟为例(图2)，详细说明本发明所提供模拟方法的具体实施步骤。

步骤一：利用gis平台和流域dem数据，提取窟野河流域河网和坡度等地形数据；

步骤二：利用gis平台和步骤一的河网提取结果，基于汇水区将温家川控制区划分为85个子流域；在每个子流域内，按照大约20km的高程间隔划分等高带，不同子流域划分出1～10个等高带不等；

步骤三：收集温家川控制区内的水文、气象、土壤、土地利用、水土保持措施、人类活动取用水等模型输入数据，并按照模型需求和步骤二中的单元划分结果，进行空间或时间尺度上的处理；如，以子流域形心点为插值点，利用距离平方倒数法空间展布降水、气温等气象要素；

步骤四：

1)等高带内，利用降水数据，按照如下公式计算雨滴溅蚀：

式中，splmod为雨滴溅蚀侵蚀模数(kg/m²/s)；e为降雨侵蚀动能(j/m²)；i为雨强(mm/min)；α为坡面的坡度(°)；a1、b1、c1、a2、b2为经验参数。

2)利用坡面流量和坡度数据，按照如下公式计算坡面径流侵蚀：

slomod＝a3·(q-b3)·(1+c3·α)

式中，slomod为坡面沟间径流侵蚀模数(kg/m²/s)；q为坡面沟间单宽流量(m²/s)；α为坡面的坡度(°)；a3、b3、c3为经验参数。

3)，每个等高带坡面所产水沙进入当前等高带的坡面沟道，按下式计算坡面沟道水流挟沙力。以挟沙力为上限，判断坡面沟道冲淤情况。当上一等高带来沙和当前等高带来沙这两个沙源总的含沙量小于水流挟沙力时，当前等高带坡面沟道发生冲刷；否则发生淤积；如此从最高等高带依次计算至最低等高带。

式中，tv为体积比挟沙力(m³/m³)；cv为体积比含沙量(m³/m³)；v为水流流速(m/s)；ρs为土粒密度(kg/m³)；ρm为浑水密度(kg/m³)；g为重力加速度(n/kg)；r为水力半径(m)；ω为浑水中的泥沙沉速(m/s)；h为水深(m)；d50为泥沙的中值粒径(m)，粒径小于该值的泥沙占全部泥沙的50％；κ为浑水的卡门常数；κ0为清水的卡门常数；a4、b4、c4、d4、a5、b5为经验参数。

步骤五：子流域计算单元中，利用坡面来沙，按照下式计算单条沟道输沙率；由于多条沟道通过并联方式汇入河道，因此，单条沟道输沙率直接累加即为沟道总输沙率；

步骤六：子流域计算单元中，利用沟道来沙及上游相邻子流域河道来沙，同样按照下式计算河道输沙率；

步骤七：坡面及坡面沟道水沙由等高带汇总到子流域；坡沟、沟道、河道各环节的水沙由子流域汇总到模拟区；最终，在模拟区内对所有环节产输沙进行汇总

式中，c和c0分别为沟道或河道出口和入口的水流含沙量(kg/m³)；t和t0分别为沟道或河道出口和入口的水流挟沙力(kg/m³)；cl为沟道/河道单位长度侧向来沙的含沙量(kg/m³)；l为河段长度(m)。

利用1968-1982年数据进行参数率定，1983-1992年数据进行模型验证，得到月均输沙率模拟结果如图3和表1所示。

表1窟野河流域温家川水文站1968-1992年月均输沙率模拟评价结果

从图3和表1可以看出：月均输沙率模拟结果与实测值比较一致，线性相关系数在0.75以上，纳什效率系数大于0.60。参考国内外研究中相关侵蚀产沙模型评价结果，该评价结果表明基于三级汇流结构的分布式侵蚀产沙模拟方法是稳定且有效的，具有较好的模拟效果，其模拟结果可以反映沟道高度发育地区的侵蚀产沙程度。

通过采用本发明公开的上述技术方案，得到了如下有益的效果：本发明提供的基于三级汇流结构的分布式侵蚀产沙模拟方法，通过改进wep-l模型，将wep-l原有的“坡面-河道”二级系统改进为“坡面-沟道-河道”三级系统，并基于侵蚀产沙和输移机制增加流域侵蚀与产沙模块，具体的，根据从坡面到流域出口的产输沙过程，按照土壤侵蚀发生的位置将水沙过程划分为坡面、沟道、河道三个环节，其中坡面侵蚀环节分为雨滴溅蚀、坡面径流侵蚀和坡面沟道产输沙过程；以各环节的侵蚀特点为依据，根据产输沙时空格局、承接传递过程等方面的相似性，分别以子流域和坡面等高带为计算单元，实现侵蚀产沙过程的多环节细化模拟，有效的反映了沟道侵蚀严重的特点，真实刻画了其侵蚀过程。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。