CN106599471B - 一种分布式流域水沙过程耦合模拟方法 - Google Patents

Abstract

本发明涉及一种分布式流域水沙过程耦合模拟方法，包括：地形数字化处理；建立等高带；建立等高带侵蚀地形数据文件；数据收集与处理；流域水文过程计算；坡面侵蚀与泥沙输移过程计算；沟道或河道或水库的水沙过程计算；判断；传递变量参数；模型参数的率定和模型验证；结束。本发明针对山坡尺度以等高带作为计算平台，以山坡真实水沙物理过程为图景，实现了等高带坡面的降雨‑径流和侵蚀‑泥沙过程的耦合模拟计算能力，在减少计算量的前提下，保证计算精度，有效还原和预测任意栅格尺度的坡面水沙过程。与已有技术相比，本发明以等高带为基本计算单元，实现地形适应性和尺度适应性；同时基于坡面水沙物理过程模拟，能够更加深入的认识坡面水沙过程规律，支撑更复杂的实践应用。

Description

一种分布式流域水沙过程耦合模拟方法

技术领域

本发明涉及一种分布式流域水沙过程耦合模拟方法，是一种水文泥沙过程计算方法，是一种使用计算机技术的水文泥沙过程模拟计算方法。

背景技术

流域的水沙过程，最终都将归结到坡面和河流的水沙过程。坡面水沙过程是复杂的多过程集合，包括降雨、植被冠层截流、地表覆被物截流、地表产流、土壤入渗、坡面汇流、河道汇流和水库水量过程，同时伴生过程包括雨滴溅蚀、薄层水流侵蚀、沟道水流侵蚀、重力侵蚀，以及泥沙输移与沉积等过程。

现有的流域水沙过程模拟多为经验法，即以水土保持试验区观测数据为基础，按照不同水土保持措施条件下降雨、径流、泥沙过程观测资料，经过统计分析，给出不同降雨条件下的径流、泥沙响应关系。或者以通过水土保持治理与非治理条件下的对比，给出不同水土保持措施的治理效益，然后按照坡面水土保持措施类型，计算不同条件下坡面水沙过程。这种以水沙过程观测资料为基础，通过统计分析等直接给出降雨、径流、泥沙过程的经验公式，然后外推应用。

流域水土保持治理是开展区域生态环境治理的基础。随着经济、社会的发展，区域水土管理越来越精细化和精准化，需要根据山坡水沙过程分异的具体条件，精准配置水土保持治理措施，才能有效改善区域水土条件，取得较好的治理成效。在水土流失区，流域水沙过程模拟也是开展大规模水土保持治理及其相关的资源、环境、生态，以及经济社会发展模式等重大决策的基础支撑。

根据已有研究，即使是相同的降雨等驱动条件下，流域水沙过程随着地形、土地利用、治理措施、治理年限等的不同，差异也较大。因此，上述水沙过程计算经验方法，没有考虑水沙过程的具体物理机理，时间和空间尺度较大，难以适应当前区域水土保持治理的需求。

与此同时，流域的水沙过程由于其巨大计算量，往往需要十分巨大的计算机硬件资源才能实现计算，甚至无法达到要求，因此，基于DEM格栅进行计算，在多数情况下运算量较大，对流域水沙过程的计算非常耗时。

发明内容

为了克服现有技术的问题，本发明提出了一种分布式流域水沙过程耦合模拟方法。所述的方法基于流域水沙物理过程，将不同分辨率的格栅整合为等高带条件下的降雨-产流-汇流和土壤侵蚀-输移机制建立一种基于等高带的分布式流域水沙过程耦合模拟方法。

本发明的目的是这样实现的：一种分布式流域水沙过程耦合模拟方法，所述方法的步骤如下：

地形数字化处理的步骤：用于以DEM数据为基础，建立等高带化地形数据；

建立等高带的步骤：用于根据Pfafstetter规则合并格栅，建立以等高带为基本计算单元的坡面水文过程计算所需的汇流关系数据文件和相应的地形数据文件；

建立等高带侵蚀地形数据文件的步骤：用于对每个等高带，根据汇流条件将计算单元划分为由平面、细沟、浅沟和切沟构成的亚计算单元地貌形态；

数据收集与处理的步骤：用于收集模拟区域内及其附近的水文、气象、植被、土壤、土地利用、水土保持措施数据，其中水文气象数据包括降雨、平均风速、平均气温、日照时数、空气相对湿度，对收集到的数据按照模拟计算所需的时间步长进行处理和空间匹配处理；

水沙计算的步骤：用于输入模型率定参数或应用参数，进行流域水文过程计算、流域侵蚀与泥沙输移过程计算和沟道或河道或水库的水沙过程计算：

所述的流域水文过程计算为：从坡面顶端的等高带开始，进行等高带中坡面水文过程计算，包括：蒸发蒸腾、入渗、地表径流、壤中径流、坡面汇流、积雪融雪计算，如果存在上游等高带则接收上游等高带的参数，加入到计算中；

所述的流域侵蚀与泥沙输移过程计算为：从坡面顶端的等高带开始，进行等高带中坡面侵蚀与泥沙输送和移动过程模拟和计算，所述的模拟和计算包括：雨滴溅蚀模拟、薄层水流侵蚀模拟、股流侵蚀过程模拟、重力侵蚀过程模拟、风力侵蚀模拟，如果存在上游等高带则接收上游等高带的参数，加入到计算中；

所述的沟道或河道或水库的水沙过程计算为：按照汇流关系数据文件对沟道或河道、水库的水沙过程进行从上游到流域出口的过程计算；

判断的步骤：用于根据汇流关系文件给定的计算关系，判断是否继续进行过程计算或结束过程计算，如果“否”则进入“计算偏差和率定参数的步骤”，如果“是”则进入“传递变量参数的步骤”；

传递变量参数的步骤：用于收集“水沙计算的步骤”所获得的变量参数，并携带这些变量参数回到“水沙计算的步骤”，继续进行水沙计算；

模型参数的率定和模型验证的步骤：如果输入的参数为模型率定参数，则通过计算偏差对水沙计算的准确度进行评价，以率定模型是否符合真实情况，并不断的调整参数并回到“水沙计算的步骤”进行计算，同时不断的进行偏差计算和比较，使模型计算与实际情况的偏差达到最小，并进行模型的验证，如果输入的参数为应用参数则跳过本步骤；

结束的步骤：用于汇集流域各个等高带的数据并输出，结束计算过程。

进一步的，所述的蒸发蒸腾计算：

式中，F _W 、F _U 、F _SV 、F _IR 、F _NI 分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的面积率；E _W 、E _U 、E _SV 、E _IR 、E _NI 分别为计算单元内水域、不透水域、裸地－植被域、灌溉农田及非灌溉农田的蒸发量或蒸发蒸腾量。

进一步的，所述的入渗计算采用Green-Ampt铅直一维入渗模型模拟降雨入渗及超渗坡面径流，和通用Green-Ampt模型进行计算。

进一步的，所述的地表径流和壤中径流的计算，

地表径流为：地表径流等于降雨减去降雨时的蒸发损失，

不透水域的地表径流按上述公式所述的地表径流和壤中径流的计算，

地表径流为：地表径流等于降雨减去降雨时的蒸发损失，

不透水域的地表径流： ；

及

计算，

式中，P为降雨，H _u 为洼地储蓄，E _u 为蒸发，R _u 为表面径流，H _umax 为最大洼地储蓄深，E _umax 为潜在蒸发，c为城市建筑物在不透水域的面积率，下标1表示城市建筑物，下标2表示城市地表面；

裸地－植被域的地表径流则根据降雨强度是否超过土壤的入渗能力分以下两种情况计算：霍顿坡面径流和饱和坡面径流，

壤中径流R2由下式计算：

式中，k(θ)为体积含水率q对应的沿山坡方向的土壤导水系数，slope为地表面坡度，L为计算单元内的河道长度，d为不饱和土壤层的厚度。

进一步的，所述的坡面汇流计算为：

运动波方程：

式中，A为流水断面面积，Q为断面流量，q _L 为网格单元或河道的单宽流入量，n为曼宁糙率系数，R为水力半径，S ₀ 为网格单元地表面坡降或河道的纵向坡降，S _f 为摩擦坡降，

动力波方程：

式中，V为断面流速，V _x 为单宽流入量的流速在x方向的分量。

进一步的，所述的积雪融雪过程计算：

式中，SM为融雪量，M _f 为融化系数，T _a 为气温指标，T ₀ 为融化临界温度，S为积雪水当量，SW为降雪水当量，E _snow 为积雪升华量。

进一步的，所述的雨滴溅蚀模拟和计算：

式中，D ₁为雨滴击溅侵蚀量，E _drop 为雨滴动能，I为雨强，J ₁为地表坡度，k ₁，α ₁，β ₁为经验参数。

其中雨动能E的计算：

式中，E _unit 为单位降雨动能，k ₁＇,α ₁＇为经验参数；

当水深大于雨滴直径3倍以上时，取水深大于0.6cm时，雨滴溅蚀作用消失；

由雨滴溅蚀增加的土壤侵蚀输沙能力计算公式为：

式中， q _s1为单宽流量输沙能力，k ₂，α ₂，β ₂为经验参数。

进一步的，所述的薄层水流侵蚀模拟：

式中， Dc为水流剥离土壤速率，k ₃为土壤的可蚀性参数，τ _f 为水流对土壤颗粒的剪切应力，τ _c 为土壤的临界抗剪切应力，Dr为细沟水流剥蚀率，q为单宽流量，c为水流泥沙含量， Tc为水流的挟沙能力，k ₄，α ₄为经验常数。

进一步的，所述的股流侵蚀过程模拟，股流的挟沙能力T _SE 计算公式为：

，

式中，k ₅为浅沟水流携沙能力系数，m为侧向汇流影响常数，ω _u 为单位水流功率；

式中，S _E 为股流侵蚀量， Q _E 为流量， D _r 为上游来沙量。

进一步的，所述的重力侵蚀过程模拟，重力侵蚀量V _g 计算公式为：

式中，k ₆为发生重力侵蚀的沟长系数，L _gully 为沟道总长度。

本发明产生的有益效果是：本发明针对山坡尺度以等高带作为计算平台，以山坡真实水沙物理过程为图景，实现了等高带坡面的降雨-径流和侵蚀-泥沙过程的耦合模拟计算能力，在减少计算量的前提下，保证了模拟计算精度，可以有效还原和预测任意等高带尺度的坡面水沙过程。与已有技术相比，本发明以等高带为基本计算单元，实现地形适应性和尺度适应性；同时基于坡面水沙物理过程模拟，能够更加深入的认识坡面水沙过程规律，支撑更多和更复杂的实践应用。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的实施例一所述方法的流程示意图；

图2是本发明的实施例十所述重力侵蚀原理示意图。

具体实施方式

实施例一：

本实施例是一种分布式流域水沙过程耦合模拟方法。本实施例以基于过程的水文模拟为基础，通过系统识别地形特征与水流动力学特征的耦合关系，建立以地形“面（片）蚀-细沟侵蚀-浅沟侵蚀-切沟侵蚀”土壤侵蚀链为平台，以“薄层水流-股流”为典型水动力条件的土壤侵蚀与泥沙输移过程模拟方法。本实施例首先通过对DEM数据的处理形成坡面计算单元拓扑关系和汇流沟道网，然后此为平台对水文气象、水土保持措施（水土保持林草、梯田、水平沟、鱼鳞坑）、土地利用、地形等信息进行处理，将处理后的数据输入模型，接着进行产汇流计算和侵蚀产沙输沙计算，最后利用已有的水文泥沙数据对模型进行率定和校核。从而实现坡面水沙过程模拟。所述实施例综合构建侵蚀地貌形态和水动力学过程特点的流域分布式水沙耦合模型。针对坡面股流侵蚀和重力侵蚀研究相对不足的现状，通过室内试验资料分析，以及黄土抗剪强度变化规律野外试验对黄土区股流侵蚀过程输沙能力以及重力侵蚀中土壤抗剪强度变化规律两个基本问题进行了研究。从而形成了物理机制相对完善的分布式流域水沙耦合模型。

所述方法的步骤如下（流程示意图见图1）：

一、地形数字化处理的步骤：用于以DEM数据为基础，建立格栅化地形数据。

二、建立等高带的步骤：用于根据Pfafstetter规则合并格栅，建立以等高带为基本计算单元的坡面水文过程计算所需的汇流关系数据文件和相应的地形数据文件。

以地形数据、汇流关系数据文件为基础，按照下述方法建立坡面侵蚀与泥沙输移过程计算所需的地形数据文件。

对于超大型流域，如黄河流域 (面积79.5万km²)，直接应用网格单元型分布式水文模型将面临严重的尺度及计算量问题。如果采用数km以内的计算网格单元，目前的PC计算机的存储量和计算速度都难以满足要求。如果采用10km以上的计算网格单元，将使河网水系和子流域边界与实际相比严重失真，带来计算时的水均衡失真和汇流路径失真问题，并造成各类水利工程难于在计算区域内合理体现。因此，对超大型流域，为克服采用小网格单元带来的计算灾难，以及采用过粗网格单元产生的计算失真问题，本实施例采用“子流域内的等高带”为基本计算单元。其编码规则是改进型的Pfafstetter规则。通过该编码规则建立了包含拓扑信息的河网生成、流域划分与编码系统，对子流域及其相应河道进行编码和计算排序，而子流域内的等高带则按照由高至低顺序进行计算。

不同侵蚀地貌的坡面分布密度运用不同侵蚀输沙规律的基础。细沟是黄土坡面分布最广的沟蚀类型之一。由于细沟密度和深度在坡面的分布随坡面长度的增大呈现多峰变化的规律，分布密度最大可以达到6%，深度变化在0-14cm之间。可以利用不同土地利用条件下细沟密度概念进行等高带内浅沟侵蚀计算。具体模拟过程中假定面（片）蚀水流在整个计算单元内发生，计算单元内依据土地利用的不同，分别给出相对于裸地的衰减系数；细沟侵蚀过程模拟中细沟尺寸依据不同土地利用类型按面积比进行概化。浅沟和切沟根据具体的地形条件决定是否发生侵蚀及沟道数量。概化的典型侵蚀地貌单元的断面参数及其发生的临界地形和水动力学条件。

三、建立等高带侵蚀地形数据文件的步骤：对于每个等高带，根据汇流条件将计算单元划分为由平面、细沟、浅沟和切沟构成的亚计算单元地貌形态。根据已有的野外调查数据，首先概括出不同侵蚀地貌的“典型侵蚀形态”。

然后确定不同侵蚀地貌的坡面分布密度：细沟分布密度0~6%，侵蚀过程模拟中细沟尺寸依据不同土地利用类型按面积比进行概化。浅沟、切沟的在对应计算单元上分布的密度采用下列公式计算：

式中，、分别为计算单元分布的浅沟和切沟数量，分别为浅沟和切沟数量计算的参数和指数。J2为计算单元比降，A为计算单元面积。

四、数据收集与处理的步骤：用于收集模拟区域内及其附近的水文、气象、植被、土壤、土地利用、水土保持措施数据，其中水文气象数据包括降雨、平均风速、平均气温、日照时数、空气相对湿度，对收集到的数据按照模拟计算所需的时间步长进行处理和空间匹配处理，同时进入“坡面水文过程计算的步骤”和“坡面侵蚀与泥沙输移过程计算的步骤”。

数据收集：收集模拟区域内及其附近的水文、气象、植被、土壤、土地利用、水土保持措施数据，其中水文气象数据包括降雨、平均风速、平均气温、日照时数、空气相对湿度；植被数据包括植被覆盖度、叶面积指数；土壤数据包括土壤类型、土壤厚度；土地利用数据包括模拟时段内多期土地利用数据，一般5年需要一期；水土保持数据包括鱼鳞坑、水平沟、梯田、水土保持林草、淤地坝等水土保持措施实施的时间、地点、范围。

数据处理：

（1）时间过程处理：对上述收集到的数据按照模拟计算所需的时间步长进行处理。对于不同数据和具体的时间过程有不同的处理方法，相关文献非常多。

（2）空间匹配处理：首先将上述数据进行空间差值和展布，形成覆盖模拟区域的空间分布数据，然后与数字化坡面的计算单元即栅格进行匹配。相关文献非常多。对于坡面侵蚀过程，根据汇流条件将计算单元划分为由平面、细沟、浅沟和切沟构成的亚计算单元地貌形态。根据已有的野外调查数据，首先对细沟、浅沟和切沟的形态进行研究，概括出不同侵蚀地貌的“典型侵蚀形态”概念。

以上步骤是针对整个模拟区域进行的地形的数字化和数据的处理。以下步骤则是针对一个等高带进行的计算和模拟。当一个等高带计算完成后则对该等高带下游的等高带进行计算和模拟，从坡顶一直计算到坡底。一个子流域可以有多个这样从坡顶到坡底的等高带序列，可以同时进行各个等高带序列的计算，也可以在计算完一个等高带序列后在计算相邻的等高带序列，以便应用两个相邻等高带需要的变量参数传递。

五、水沙计算的步骤：用于输入模型率定参数或应用参数，进行流域水文过程计算、流域侵蚀与泥沙输移过程计算和沟道或河道或水库的水沙过程计算：

所述的模型率定参数是指在模型建立过程中，首先要对模型参数是否符合实际请进行率定，因此先要根据对流域水文泥沙过程的认识输入一套模型参数，将计算出的径流泥沙过程与实际观测的径流泥沙过程进行对比，并在此基础上有针对性的对模型参数进行调整，从而不断的减小模型偏差，使模型的计算结果与实际结果的偏离度最小，以此建立合适的模型参数，这些率定模型的已知结果的参数称之为“模型参数的率定”。

模型的验证：经过率定的模型参数还需要选取其他时段进行流域水沙过程计算，然后再次将计算结果与实际水文泥沙过程进行对比，如果模型运算结果仍然能保持精度，则模型建立完成。如果不能保证，则还需要返回上一步进行模型参数的率定。

所述的应用参数是在模型经过率定后，已经到达了与实际十分接近的程度，则输入一些预测或其他应用的参数，以获取一个未知的结果，这些参数称之为“应用参数”。

所述的流域水文过程计算为：从坡面顶端的等高带开始，进行等高带中坡面水文过程计算，包括：蒸发蒸腾、入渗、地表径流、壤中径流、坡面汇流、积雪融雪计算，如果存在上游等高带则接收上游等高带的参数，加入到计算中。

所述的流域侵蚀与泥沙输移过程计算为：从坡面顶端的等高带开始，进行等高带中坡面侵蚀与泥沙输送和移动过程模拟和计算，所述的模拟和计算包括：雨滴溅蚀模拟、薄层水流侵蚀模拟、股流侵蚀过程模拟、重力侵蚀过程模拟、风力侵蚀模拟，如果存在上游等高带则接收上游等高带的参数，加入到计算中。

所述的沟道或河道或水库的水沙过程计算为：按照汇流关系数据文件对沟道或河道、水库的水沙过程进行从上游到流域出口的过程计算。

流域水文过程计算：从坡面顶端的等高带开始，进行等高带中坡面水文过程计算，包括：蒸发蒸腾、入渗、地表径流、壤中径流、坡面汇流、积雪融雪计算、河道或沟道或水库泥沙输移过程模拟，如果存在上游等高带则接收上游等高带的参数，加入到计算中。

坡面水沙过程的模拟是基于物理机制分布式水沙耦合模型的基础。从侵蚀动力学角度，黄土区坡面典型水沙过程主要包括雨滴溅蚀过程、薄层水流侵蚀过程、股流侵蚀过程和重力侵蚀过程。其中面（片）蚀和细沟侵蚀适用于薄层水流侵蚀模拟，浅沟侵蚀和切沟侵蚀适用于股流侵蚀。目前对雨滴溅蚀和薄层水流侵蚀过程的模拟研究已经较为成熟，而对股流侵蚀和重力侵蚀过程的机理描述尚不完善。本实施例采用雨滴溅蚀和薄层水流侵蚀过程模，同时为了建立完善的坡面水沙过程模拟模型，对股流侵蚀过程和重力侵蚀过程在物理图景概化的基础上，通过实验研究建立了基于物理机制的模拟子模型。

蒸发蒸腾的计算：可以采用“马赛克”结构考虑网格单元内的土地利用变异问题，每个网格单元的蒸发蒸腾可能包括植被截留蒸发、土壤蒸发、水面蒸发和植被蒸腾等多项。参照土壤-植被-大气通量交换方法(SVATS)中的ISBA模型，采用Penman公式或Penman-Monteith公式等进行计算。同时，由于蒸发蒸腾过程和能量交换过程客观上融为一体，为计算蒸发蒸腾，地表附近的辐射、潜热、显热和热传导的计算不可缺少，而这些热通量又均是地表温度的函数。为减轻计算负担，模型对热传导及地表温度的计算采用强制复原法。

入渗的计算：可以采用Green-Ampt铅直一维入渗模型模拟降雨入渗及超渗坡面径流。同时，考虑到由自然力和人类活动(如农业耕作)等引起的土壤分层问题，模型采用Jia和Tamai提出的实际降雨条件下的多层Green-Ampt模型，以下称通用Green-Ampt模型。

地表径流的计算：水域的地表径流等于降雨减去降雨时的蒸发损失，裸地－植被域（透水域）的地表径流则根据降雨强度是否超过土壤的入渗能力分以霍顿坡面径流饱和坡面径流进行计算。

壤中径流计算：在山地丘陵等地形起伏地区，同时考虑坡向壤中径流及土壤渗透系数的各向变异性。

坡面汇流计算：在河网水系生成和流域划分的基础上，根据网格单元DEM及土地利用等基本信息准备各子流域等高带的属性表（包括面积、长度、宽度、平均高程、坡度和曼宁糙率等），采用一维运动波法从上游等高带至下游等高带计算坡面汇流，并将下游等高带的坡面汇流输入给所在子流域内的河道。

积雪融雪计算：可以采用“温度指标法”计算积雪融雪的日变化过程，或其他类似的方法。

风力侵蚀过程模拟：风力侵蚀输沙量对于宁蒙等黄河流域典型风沙区的风力侵蚀过程采用下列公式计算：

式中，Q为输沙量；和分别为气流摩阻风速和临界起沙风速；为空气密度；为重力加速度；为沙粒直径；为标准沙粒直径；为系数，取4.2。

坡面侵蚀与泥沙输移过程计算：从坡面顶端的等高带开始，进行等高带中坡面侵蚀与泥沙输送和移动过程模拟和计算，所述的模拟和计算包括：雨滴溅蚀模拟、薄层水流侵蚀模拟、股流侵蚀过程模拟、重力侵蚀过程模拟，如果存在上游等高带则接收上游等高带的参数，加入到计算中。

坡面水沙过程的模拟是基于物理机制分布式水沙耦合模型的基础。从侵蚀动力学角度，黄土区坡面典型水沙过程主要包括雨滴溅蚀过程、薄层水流侵蚀过程、股流侵蚀过程和重力侵蚀过程。其中面（片）蚀和细沟侵蚀适用于薄层水流侵蚀模拟，浅沟侵蚀和切沟侵蚀适用于股流侵蚀。目前对雨滴溅蚀和薄层水流侵蚀过程的模拟研究已经较为成熟，而对股流侵蚀和重力侵蚀过程的机理描述尚不完善。

雨滴溅蚀模拟：对于土壤性质基本一致的流域，影响降雨溅蚀的主要因素为降雨动能、地表坡度。等高带内可依据土地利用的不同，分别给出相对于裸地的衰减系数。其中对于裸土地雨滴溅蚀的计算选取适合于黄土地区，可运用综合考虑降雨及坡度的模型—吴普特模型进行溅蚀计算。

薄层水流侵蚀模拟：可采用Foster和Meyer年提出的关系式描述薄层水流土壤分离速率与输沙率之间的关系。该关系式被用于WEPP模型，且经实验验证这一假定符合黄土区薄层水流侵蚀过程规律。

股流侵蚀过程模拟：股流挟沙能力是指水流能量全部用于泥沙输移时的水流最大挟沙能力。如果使用经过筛分除去石砾的松散土壤，土壤的抗蚀能力下降较自然土壤明显下降，由于黄土的易蚀性，水流的侵蚀耗能很少减少，在这种情况下水流含沙量可以近似看做股流挟沙能力。

重力侵蚀过程模拟：坡面水沙过程中的重力侵蚀，其影响因素相对单一。对于坡面坡度大于黄土休止角的土体，其是否发生重力侵蚀取决于土体的下滑力G _x 与土体抗剪强度τ _c 之间的大小关系。当G _x ＞τ _c 重力侵蚀发生。

沟道或河道或水库的水沙过程计算：按照汇流关系数据文件对沟道或河道、水库的水沙过程进行从上游到流域出口的过程计算。水流过程计算采用动力波方程计算，所述的泥沙过程采用沿水深积分后的一维恒定水流泥沙扩散方程计算。其中对于水流挟沙能力的计算，在沟道中采用费祥俊公式，在河道和水库中采用张红武公式。

河道或沟道、水库泥沙输移过程模拟，河道泥沙过程采用沿水深积分后的一维恒定水流泥沙扩散方程：

式中：S _x ，S _x* 分别为断面的平均含沙量和水流挟沙力，q为单宽流量，α为恢复饱和系数，ω _s 为河流泥沙沉速。上式的积分形式为河段出口含沙量的显式计算公式为：

式中S、S ₀分别为出口断面和进口断面的平均含沙量；S _* 、S _0* 分别为出口断面和进口断面的水流挟沙力；L为河段长度，采用下式计算：

式中，d ₉₀为泥沙的上限粒径，μ为浑水粘度，ρ _m 为含沙水流密度，S _v 为体积比含沙量，S _vm 为极限体积比含沙量，μ ₀为水的粘度，k ₈为泥沙固体浓度修正系数，d _l ，P _l 分别为某一粒径级的平均直径及其相应的重量百分比。

对于沟道水流挟沙能力的计算，采用费祥俊公式计算：

式中，U为断面平均流速，f为达西系数；R为水力半径；ω ₉₀为上限粒径在一定浓度下的沉速。

式中，k _s 为河床糙度，取为2d ₉₀；R _e 为计算的雷诺数，μ为浑水粘度；γ _m 为含沙水流容重。

对于河道和水库水流挟沙能力计算，采用张红武公式计算：

式中，为断面流速；γ _s 为泥沙容重；κ为浑水卡门常数，取0.4。

六、判断的步骤：用于根据汇流关系文件给定的计算关系，判断是否继续进行过程计算或结束过程计算，如果“否”则进入“计算偏差和率定参数的步骤”，如果“是”则进入“传递变量参数的步骤”。

由于等高带DEM数据是从坡面到坡底顺序排列，上下游等高带之间有变量参数传递的关系，因此，计算开始时，首先应当从坡顶开始计算，然后顺序向下计算，并不断的将变量参数向下一个等高带传递，所以一个等高带计算完毕后就要判断一下，是否达到了流域出口，如果没有达到流域出口，就要继续进行计算。如果已经达到了流域出口则不再继续进行计算。

七、传递变量参数的步骤：用于收集“水沙计算的步骤”所获得的变量参数，并携带这些变量参数回到“水沙计算的步骤”，继续进行水沙计算。

为进行下一个等高带的水沙过程计算，需要将本次等高带的计算所获得的变量参数汇集并传递到下一个等高带的水深过程计算中，即将本次计算中的等高带水沙的影响传递到下一个等高带中。

八、模型参数的率定和模型验证的步骤：如果输入的参数为模型率定的已知参数（有结果记录的参数），则通过计算偏差对水沙计算的准确度进行评价，以率定模型是否符合真实情况，并不断的调整参数并回到“水沙计算的步骤”进行计算，同时不断的进行偏差计算和比较，使模型计算与实际情况的偏差达到最小，并进行模型的验证，如果输入的参数为应用参数则跳过本步骤。

模型参数的率定和模型验证是一个反复进行的过程，需要不断的进行偏差计算和调整、修正参数。参数的调整十分复杂，需要调整的参数往往很多，调整哪些些参数和调整的幅度控制多少，才能尽快准确的达到目的，需要智慧和耐力。

模型参数率定过程是按照模型框架需求，对模型参数进行调整使模型准确反映研究对象真实世界客观规律的过程。模型验证则是对模型参数率定结果的可靠性和准确性检验。验证标准一般有Nash效率、相关系数、相对误差等。评价模型校验的好坏，主要采用指标为： Nash效率、相关系数、相对误差等。为了使模型能够更好的模拟水沙过程，进行调参数时优先保证Nash效率、相关系数相对较高。

Nash-Sutcliffe效率：

Nash与Sutcliffe在1970年提出了模型效率系数（也称确定性系数）来评价模型模拟结果的精度，它更直观地体现了实测过程与模型模拟过程拟合程度的好坏，公式如下：

式中：Nash为Nash-Sutcliffe效率系数，其值越接近于1表示实测与模拟流量过程拟合得越好，模拟精度越高；为模拟值，为实测值，为实测平均值。

相关系数：

相关系数是对两个变量之间关系的量度，考查两个事物之间的关联程度。相关系数的绝对值越大，相关性越强，相关系数越接近于1和－1，相关度越强，相关系数越接近于0，则相关度越弱。其计算公式如下：

式中：rxy为相关系数；n为系列的样本数；X、Y分别代表实测系列和模拟系列的数值。通常情况下：|rxy|在0.8-1.0之间为极强相关，在0.6-0.8之间为强相关，在0.4-0.6之间为中等程度相关，在0.2-0.4之间为弱相关，在0-0.2之间为极弱相关或无相关。

相对误差：

相对误差是整个模拟期模拟值与实测值的差值与实测值的百分比，径流量误差绝对值越接近于零越好。

式中：Dv为模拟相对误差（％）；F0为实测值均值；R为模拟均值。

九、结束的步骤：用于汇集流域各个等高带的数据并输出，结束计算过程。输出数据并结束计算过程。

实施例二：

本实施例是实施例一的改进，是实施例一关于蒸发蒸腾计算的细化。本实施例所述的蒸发蒸腾计算：

计算单元内（等高带内）的蒸发蒸腾包括来自植被湿润叶面（植被截留水）、水域、土壤、城市地表面、城市建筑物等的蒸发，以及来自植被干燥叶面的蒸腾。计算单元的平均蒸发蒸腾量模型采用下式算出：

式中，F _W 、F _U 、F _SV 、F _IR 、F _NI 分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的面积率（%）；E _W 、E _U 、E _SV 、E _IR 、E _NI 分别为计算单元内水域、不透水域、裸地－植被域、灌溉农田及非灌溉农田的蒸发量或蒸发蒸腾量。

水域的蒸发量（Ew）由下述Penman公式算出：

式中，RN为净放射量；G为传入水中的热通量；Δ为饱和水蒸气压对温度的导数；δe为水蒸气压与饱和水蒸气压的差；ra为蒸发表面的空气动力学阻抗；ρa为空气的密度；Ｃp为空气的定压比热；λ为水的气化潜热，γ=Ｃp/λ。

裸地－植被域蒸发蒸腾量（ESV）、灌溉农田域（EIR）、和非灌溉农田域（ENI）分别由以下公式计算:

(3)

EIR=Ei1+Ei2+Etr3+Es (4)

ENI= Ei1+Ei2+Etr4+Es (5)

式中，Ei为植被截留蒸发（来自湿润叶面）；Etr为植被蒸腾（来自干燥叶面）；Es为裸地土壤蒸发。另外，下标1表示高植被（森林、城市树木），下标2表示草，下标3表示灌溉农作物，下标4表示非灌溉农作物。

各类植被的截留蒸发（Ei）使用ISBA模型计算：

(6)

(7)

(8)

(9)

(10)

式中，Veg为植被覆盖度；d为湿润叶面积占总叶面积的比例；Ep为可能蒸发量（由Penman方程式计算）；Wr为植被截留水量；P为降雨量；Rr为植被流出水量；Wrmax为最大植被截留水量；LAI为叶面积指数。

植被蒸腾由Penman－Monteith公式计算。

(11)

(12)

式中，RN为净放射量；G为传入植被体内的热通量；rc为植物群落阻抗（canopyresistance）。蒸腾属于土壤、植物、大气连续体（SPAC：Soil－Plant－AtmosphereContinuum）水循环过程的一部分，受光合作用、大气湿度、土壤水分等的制约。这些影响通过式（12）中的植物群落阻抗（rc）来考虑。

植被蒸腾是通过根系吸水由土壤层供给。假定根系吸水率随深度线性递减、根系层上半部的吸水量占根系总吸水量的70%，则可得下式：

(13)

(14)

式中， Etr为蒸腾量；为根系层的厚度；z为离地表面的深度； Sr(z)为深度处的根系吸水强度；Etr (z)为从地表面到深度z处的根系吸水量。

根据以上公式，只要给出植物根系层厚，即可算出其从土壤层各层的吸水量（蒸腾量）。在本研究中，认为草与农作物等低植物的根系分布于土壤层的一、二层，而树木等高植物的根系分布于土壤层的所有三层。结合土壤各层的水分移动模型，即可算出各层的蒸腾量。

裸地土壤蒸发由下述修正Penman公式计算：

(15)

(16)

式中，b为土壤湿润函数或蒸发效率；q为表层（一层）土壤的体积含水率；qfc为表层土壤的田间持水率；qm为土壤单分子吸力（约1000~10000个大气压）对应的土壤体积含水率。

不透水域的蒸发及地表径流由下述方程式求解：

(17)

(18)

(19)

(20)

(21)

(22)

(23)

式中，P为降雨；Hu为洼地储蓄；Eu为蒸发；Ru为表面径流；Humax为最大洼地储蓄深；Eumax为潜在蒸发（由Penman公式计算）；c为城市建筑物在不透水域的面积率；下标1表示城市建筑物，下标2表示城市地表面。

实施例三：

本实施例是上述实施例的改进，是上述实施例关于入渗计算的细化。本实施例所述的入渗计算采用Green-Ampt铅直一维入渗模型模拟降雨入渗及超渗坡面径流，和通用Green-Ampt模型进行计算。

当入渗湿润锋到达第m土壤层时入渗能力由下式计算：

(24)

式中，f为入渗能力；F为累积入渗量；km、Am-1、Bm-1见后述。累积入渗量F的算出方法，视地表面有无积水而不同。

如果自入渗湿润锋进入第m-1土壤层时起地表面就持续积水，那么累积入渗量由式（25）计算；如果前一时段tn-1地表面无积水，而现时段tn地表面开始积水，那么由式（26）计算。

(25)

(26)

(27)

(28)

(29)

(30)

(31)

式中，SW为入渗湿润锋处的毛管吸力；k为土壤层的导水系数；qs为土壤层的含水率；q0为土壤层的初期含水率；t为时刻；Fp为地表面积水时的累积入渗量；tp为积水开始时刻；Ip为积水开始时的降雨强度；tm-1为入渗湿润锋到达第m层与第m-1层交界面的时刻；L为入渗湿润锋离地表面的深度；Li为第i层的厚度；Dq为qs－q0。

土壤水分吸力关系采用Havercamp公式：

(32)

式中，q为土壤体积含水率；qs为饱和含水率；qr为残留含水率；j为吸引压，cm水柱；a和b为常数。

土壤水分导水系数关系采用Mualem公式：

(33)

式中，Ks为土壤饱和导水系数，cm/s；k(q)为含水率q对应的导水系数，cm/s；n为常数。

实施例四：

本实施例是上述实施例的改进，是上述实施例关于地表径流和壤中径流计算的细化。本实施例所述的地表径流和壤中径流的计算。

地表径流：

水域的地表径流等于降雨减去降雨时的蒸发损失，不透水域的地表径流按上述公式（20）及（22）计算，而裸地－植被域（透水域）的地表径流则根据降雨强度是否超过土壤的入渗能力分以下两种情况计算。

（1）霍顿坡面径流（Hortonian overland flow）

当降雨强度超过土壤的入渗能力时将产生这类地表径流R1ie，即超渗产流，由下式计算：

(34)

(35)

式中，P为降水量；HSV为裸地－植生域的洼地储蓄；HSVMax为最大洼地储蓄深；ESV为蒸散发；fSV为由通用Green-Ampt模型算出的土壤入渗能力。

（2）饱和坡面径流（Saturation overland flow）

对于河道两岸及低洼的地方，由于地形的作用，土壤水及浅层地下水逐渐汇集到这些地方，土壤饱和或接近饱和状态后遇到降雨便形成饱和坡面径流（蓄满产流）。此时，Green-Ampt模型已无能为力，需根据非饱和土壤水运动的Richards方程来求解。为减轻计算负担，地表洼地储留层按连续方程、表层土壤分成3层按照Richards方程（积分形式）进行计算：

地表洼地储留层

(36)

(37)

土壤表层

(38)

土壤中层

(39)

土壤底层

(40)

(j：1、3) (41)

(42)

(43)

(j：1、2) (44)

(j：1、2) (45)

式中，Hs为洼地储蓄；Hsmax为最大洼地储蓄；Veg1、Veg2为裸地－植生域的高植生和低植生的面积率；Rr1、Rr2为从高植生和低植生的叶面流向地表面的水量；Q为重力排水；QDj,j+1为吸引圧引起的j 层与j+1 层土壤间的水分扩散；E0为洼地储蓄蒸发；Es为表层土壌蒸发；Etr为植被蒸散（第一下标中的1表示高植生、2表示低植生；第一下标表示土壤层号）；R2为壤中径流；k(q)为体积含水率q对应的土壤导水系数；y(q)为体积含水率q对应的土壤吸引力；d为土壤层厚度；W为土壤的蓄水量，w=qd）；W10为表层土壤的初期蓄水量。另外、下标0、1、2、3分别表示洼地储蓄层、表层土壤层、第2土壤层和第3土壤层。

壤中径流：

在山地丘陵等地形起伏地区，同时考虑坡向壤中径流及土壤渗透系数的各向变异性。从不饱和土壤层流入河道的壤中径流由下式计算：

(46)

式中，k(θ)为体积含水率q对应的沿山坡方向的土壤导水系数；slope为地表面坡度；L为计算单元内的河道长度；d为不饱和土壤层的厚度。

实施例五：

本实施例是上述实施例的改进，是上述实施例关于坡面汇流计算的细化。本实施例所述的坡面汇流计算为：

坡面汇流。模型则在河网水系生成和流域划分的基础上，根据网格单元DEM及土地利用等基本信息准备各子流域等高带的属性表（包括面积、长度、宽度、平均高程、坡度和曼宁糙率等），采用一维运动波法从上游等高带至下游等高带计算坡面汇流，并将下游等高带的坡面汇流输入给所在子流域内的河道。

运动波方程：

（连续方程） (50)

（运动方程） (51)

（Manning公式） (52)

式中，A为流水断面面积；Q为流量；qL为网格单元或河道的单宽流入量（包含网格内的有效降雨量、来自周边网格及支流的水量）；n为Manning糙率系数；R为水力半径；S ₀ 为网格单元地表面坡降或河道的纵向坡降；S_f为摩擦坡降。

动力波方程（Saint Venant 方程）：

（连续方程） （53)

（运动方程） (54)

（Manning公式） (55)

式中，A为流水断面面积；Q为断面流量；q _L 为网格单元或河道的单宽流入量（包含网格内的有效降雨量、来自周边网格及支流的水量）；n为Manning糙率系数；R为水力半径；S ₀ 为网格单元地表面坡降或河道的纵向坡降；S _f 为摩擦坡降；V为断面流速；V _x 为单宽流入量的流速在x方向的分量。

实施例六：

本实施例是上述实施例的改进，是上述实施例关于积雪融雪过程计算的细化。本实施例所述的积雪融雪过程计算：

式中，SM为融雪量，mm/day；M _f 为融化系数或称“度日因子”，mm/°C/day；T _a 为气温指标，°C；T ₀ 为融化临界温度，°C；S为积雪水当量，mm；SW为降雪水当量，mm；E _snow 为积雪升华量，mm。

“度日因子”既随海拔高度和季节变化，又随下垫面条件变化，常做为模型调试参数对待。一般情况下在1 mm/°C/day和7 mm/°C/day之间，且裸地高于草地，草地高于森林。气温指标通常取为日平均气温。融化临界温度通常在-3°C和0°C之间。另外，为将降雪与降雨分离，还需要雨雪临界温度参数（通常在0°C和3°C之间）。

实施例七：

本实施例是上述实施例的改进，是上述实施例关于雨滴溅蚀模拟和计算的细化。本实施例所述的雨滴溅蚀模拟和计算：

式中，D ₁为雨滴击溅侵蚀量，g/m²；E _drop 为雨滴动能，J/m²；I为雨强，mm/min；J ₁为地表坡度，°；k ₁，α ₁，β ₁为经验参数。

其中雨动能的计算：

式中，E _unit 为单位降雨动能，J/(m²·mm)；k ₁＇,α ₁＇为经验参数；

当水深大于雨滴直径3倍以上时，取水深大于0.6cm时，雨滴溅蚀作用消失；

由雨滴溅蚀增加的土壤侵蚀输沙能力计算公式为：

（60）

式中，为单宽流量输沙能力，kg/m²；为经验参数

实施例八：

本实施例是上述实施例的改进，是上述实施例关于薄层水流侵蚀模拟的细化。本实施例所述的薄层水流侵蚀模拟：

（61）

（62）

（63）

式中，为水流剥离土壤速率；为土壤的可蚀性参数；为水流对土壤颗粒的剪切应力；为土壤的临界抗剪切应力；为细沟水流剥蚀率；为单宽流量；为水流泥沙含量； 为水流的挟沙能力；为经验常数。

实施例九：

本实施例是上述实施例的改进，是上述实施例关于股流侵蚀过程模拟的细化。本实施例所述的股流侵蚀过程模拟，股流的挟沙能力计算公式为：

（70）

式中，为浅沟水流携沙能力系数；为侧向汇流影响常数；

（71）

式中，为股流侵蚀量；为流量；为股流携沙能力、为上游来沙量。

实施例十：

本实施例是上述实施例的改进，是上述实施例关于重力侵蚀过程模拟的细化。本实施例所述的重力侵蚀过程模拟，重力侵蚀量计算公式为：

坡面水沙过程中的重力侵蚀，其影响因素相对单一。对于坡面坡度大于黄土休止角的土体，其是否发生重力侵蚀取决于土体的下滑力与土体抗剪强度之间的大小关系。当重力侵蚀发生。如图2所示，土体下滑力可用下式计算：

（72）

式中， Δh为发生重力侵蚀土体深度，细沟沟坡取0.1m，浅沟沟坡取0.5m，切沟沟坡取10m；θ为黄土天然休止角；γ _s 为土壤湿容重。

重力侵蚀量计算公式为：

（73）

式中，为发生重力侵蚀的沟长系数，为沟道总长度，其它参数意义同上。

土壤的临界抗剪强度公式为：

，n=78，r=0.61 （74）

由于坡面侵蚀过程中重力侵蚀发生的尺度相对较小，影响因素相对单一。采用土体的下滑力与土体的抗剪强度大小作为判定条件。其中土体的下滑力由式（72）计算，土体的抗剪强度由式（74）计算，重力侵蚀量由式（73）计算。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如公式的运用、步骤的先后顺序等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (10)

1.一种分布式流域水沙过程耦合模拟方法，其特征在于，所述方法的步骤如下：

步骤1：地形数字化处理的步骤：用于以DEM数据为基础，建立格栅化地形数据；

步骤2：建立等高带的步骤：用于根据Pfafstetter规则合并格栅，建立以等高带为基本计算单元的坡面水文过程计算所需的汇流关系数据文件和相应的地形数据文件，通过该规则建立了包含拓扑信息的河网生成、流域划分与编码系统，对子流域及其相应河道进行编码和计算排序，而子流域内的等高带则按照由高至低顺序进行计算；

步骤3：建立等高带侵蚀地形数据文件的步骤：用于对每个等高带，根据汇流条件将计算单元划分为由平面、细沟、浅沟和切沟构成的亚计算单元地貌形态；

步骤4：数据收集与处理的步骤：用于收集模拟区域内及其附近的水文、气象、植被、土壤、土地利用、水土保持措施数据，其中水文气象数据包括降雨、平均风速、平均气温、日照时数、空气相对湿度，对收集到的数据按照模拟计算所需的时间步长进行处理和空间匹配处理；

步骤5：水沙计算的步骤：用于输入模型率定参数或应用参数，进行流域水文过程计算、流域侵蚀与泥沙输移过程计算和沟道或河道或水库的水沙过程计算；

所述的流域水文过程计算为：从坡面顶端的等高带开始，进行等高带中坡面水文过程计算，包括：蒸发蒸腾、入渗、地表径流、壤中径流、坡面汇流、积雪融雪计算，如果存在上游等高带则接收上游等高带的参数，加入到计算中；

所述的流域侵蚀与泥沙输移过程计算为：从坡面顶端的等高带开始，进行等高带中坡面侵蚀与泥沙输送和移动过程模拟和计算，所述的模拟和计算包括：雨滴溅蚀模拟、薄层水流侵蚀模拟、股流侵蚀过程模拟、重力侵蚀过程模拟、风力侵蚀模拟，如果存在上游等高带则接收上游等高带的参数，加入到计算中；

所述的沟道或河道或水库的水沙过程计算为：按照汇流关系数据文件对沟道或河道、水库的水沙过程进行从上游到流域出口的过程计算；

步骤6：判断的步骤：用于根据汇流关系数据文件给定的计算关系，判断是否继续进行过程计算或结束过程计算，如果“否”则进入“模型参数的率定和模型验证的步骤”，如果“是”则进入“传递变量参数的步骤”；

步骤7：传递变量参数的步骤：用于收集“水沙计算的步骤”中本次等高带的计算所获得的变量参数，并携带这些变量参数回到“水沙计算的步骤”，从而将本次计算中的等高带水沙的影响传递到下一个等高带中；

步骤8：模型参数的率定和模型验证的步骤：如果步骤5“用于输入模型率定参数或应用参数”中所获得的参数为模型率定参数，则通过计算偏差对水沙计算的准确度进行评价，以率定模型是否符合真实情况，并不断的调整参数并回到“水沙计算的步骤”进行计算，同时不断的进行偏差计算和比较，使模型计算与实际情况的偏差达到最小，并进行模型的验证，如果步骤5“用于输入模型率定参数或应用参数”中所获得的参数为应用参数则跳过本步骤；

步骤9：结束的步骤：用于汇集流域各个等高带的数据并输出，结束计算过程。

2.根据权利要求1所述的模拟方法，其特征在于，所述的蒸发蒸腾计算：

E＝F_WE_W+F_UE_U+F_SVE_SV+F_IRE_IR+F_NIE_NI

式中，F_W、F_U、F_SV、F_IR、F_NI分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的面积率；E_W、E_U、E_SV、E_IR、E_NI分别为计算单元内水域、不透水域、裸地－植被域、灌溉农田及非灌溉农田的蒸发量或蒸发蒸腾量。

3.根据权利要求1所述的模拟方法，其特征在于，所述的入渗计算采用Green-Ampt铅直一维入渗模型和通用Green-Ampt模型进行计算，其中所述的Green-Ampt铅直一维入渗模型用于模拟降雨入渗及超渗坡面径流。

4.根据权利要求1所述的模拟方法，其特征在于，所述的地表径流和壤中径流的计算，

地表径流为：地表径流等于降雨减去降雨时的蒸发损失，

不透水域的地表径流：

及

计算，

式中，P为降雨，H_u为洼地储蓄，E_u为蒸发，R_u为地表径流，H_umax为最大洼地储蓄深，E_umax为潜在蒸发，下标1表示城市建筑物，下标2表示城市地表面；

裸地－植被域的地表径流则根据降雨强度是否超过土壤的入渗能力分以下两种情况计算：霍顿坡面径流和饱和坡面径流，

壤中径流R2由下式计算：

R2＝k(θ)sin(slope)Ld

式中，k(θ)为体积含水率q对应的沿山坡方向的土壤导水系数，slope为地表面坡度，L为计算单元内的河道长度，d为不饱和土壤层的厚度。

5.根据权利要求1所述的模拟方法，其特征在于，所述的坡面汇流计算为：

运动波方程：

S_f＝S₀

式中，A为流水断面面积，Q为断面流量，q_L为网格单元或河道的单宽流入量，n为曼宁糙率系数，R为水力半径，S₀为网格单元地表面坡降或河道的纵向坡降，S_f为摩擦坡降，

动力波方程：

式中，V为断面流速，V_x为单宽流入量的流速在x方向的分量。

6.根据权利要求1所述的模拟方法，其特征在于，所述的积雪融雪过程计算：

SM＝M_f(T_a-T₀)

式中，SM为融雪量，M_f为融化系数，T_a为气温指标，T₀为融化临界温度，S为积雪水当量，SW为降雪水当量，E_snow为积雪升华量。

7.根据权利要求1所述的模拟方法，其特征在于，所述的雨滴溅蚀模拟和计算：

式中，D₁为雨滴击溅侵蚀量，E_drop为雨滴动能，I为雨强，J₁为地表坡度，k₁，α₁，β₁为经验参数；

其中雨动能的计算：

式中，E_unit为单位降雨动能，k₁＇,α₁＇为经验参数；

由雨滴溅蚀增加的土壤侵蚀输沙能力计算公式为：

式中，q_s1为单宽流量输沙能力，k₂，α₂，β₂为经验参数。

8.根据权利要求1所述的模拟方法，其特征在于，所述的薄层水流侵蚀模拟：

Dc＝k₃(τ_f-τ_c)

式中，Dc为水流剥离土壤速率，k₃为土壤的可蚀性参数，τ_f为水流对土壤颗粒的剪切应力，τ_c为土壤的临界抗剪切应力，Dr为细沟水流剥蚀率，q为单宽流量，c为水流泥沙含量，Tc为水流的挟沙能力，k₄，α₄为经验常数。

9.根据权利要求1所述的模拟方法，其特征在于，所述的股流侵蚀过程模拟，股流的挟沙能力T_SE计算公式为：

式中，k₅为浅沟水流携沙能力系数，m为侧向汇流影响常数，ω_u为单位水流功率；

式中，S_E为股流侵蚀量，Q_E为流量，D_r为上游来沙量。

10.根据权利要求1所述的模拟方法，其特征在于，所述的重力侵蚀过程模拟，重力侵蚀量V_g计算公式为：

式中，G_x为土体的下滑力，τ_c为土体抗剪强度，k₆为发生重力侵蚀的沟长系数，L_gully为沟道总长度，Δh为发生重力侵蚀土体深度，θ为黄土天然休止角。