一种水土保持工程减水减沙效益定量评价方法与流程

本发明涉及一种水土保持工程减水减沙效益定量评价方法，是一种水文计算方法，是一种使用计算机技术的水文模拟和计算方法。

背景技术：

现有的水土保持水文水资源效应研究所采取的方法均为“水文法”和“水保法”。由于这 两种方法都不能系统地揭示水土保持对水循环影响的物理机制,并且采用的水资源评价口径。只是单一的狭义水资源,使得研究成果不能真实反映水土保持的水文水资源效应。以狭义水资源和广义水资源为评价口径的基于物理机制的分布式 水文模型则能克服水文法和水保法的严重弊端。

“水文法”也叫统计分析法，利用流域水文泥沙观测资料分析水土保持措施蓄水拦沙作用的一种方法。基本原理是基于对降雨产流产沙规律的分析，建立水土保持治理前流域产流产沙关系式，然后将治理后的降雨条件带入关系式，计算得到该降雨条件时流域未治理条件下的产水产沙量，再与治理后的实测产水产沙量比较，从而求得水土保持措施对流域水沙的影响量。

“水保法”也叫成因分析法，是根据不同水土保持措施减水减沙作用实测成果，考虑产沙在河道中的输移及冲淤变化，以及新增水土流失数量等，计算水土保持减水减沙效益的一种方法。一般直接采用具体水土保持措施的面积和该项措施的蓄水保土定额之积作为该项水土保持措施的减水减沙效益。

上述两种方法的问题都在于不能精确的计算或模拟坡面的水沙过程，因此，大大的影响了计算或模拟的精确度。

技术实现要素：

为了克服现有技术的问题，本发明提出了一种水土保持工程减水减沙效益定量评价方法。所述的方法基于坡面或流域水沙物理过程，或者将不同分辨率的栅格作为计算单元，或者将栅格整合为等高带作为计算单元，建立计算单元条件下的降雨-产流-汇流和土壤侵蚀-输移机制的水沙过程，用以水土保持工程减水减沙效益定量评价的模拟方法。

本发明的目的是这样实现的：一种水土保持工程减水减沙效益定量评价方法，所述方法的步骤如下：

地形数字化处理的步骤：用于DEM数据为基础，以ArcGIS水文模块为主要工具，采用D8法，建立以栅格或等高带为基本的计算单元的坡面水文过程计算所需的地形数据、汇流关系数据文件，以这些文件为基础，建立坡面侵蚀与泥沙输移过程计算所需的地形数据文件；

水土保持措施数据收集的步骤：用于收集模拟区域内不同时间的坡面植被、土壤、土地利用情况，以及水土保持措施情况，形成不同时间段的水土保持数据集；

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

水沙计算的步骤：用于输入模型率定参数或应用参数，进行坡面水文过程计算、坡面侵蚀与泥沙输移过程计算：

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

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

判断的步骤：根据汇流关系文件给定的计算关系， 进行过程计算和决定是否输出当前计算结果，如果“是”进入“计算偏差和率定参数的步骤”，如果“否”则进入“传递变量参数的步骤”；

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

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

设置水土保持情景的步骤：用于使用水土保持措施数据，设置多个年代的水土保持措施情景；

生成不同水土保持情景下的水沙过程数据的步骤：用于将设计的水土保持情景数据运行计算，得到多个水土保持措施的水沙过程数据；

分析计算不同水土保持情景的保水保土效益的步骤：用于对比不同水土保持情景下的水沙过程变化量进行分析；

比较并结束的步骤：用于对各个情景水土保持情景的水土保持数据集所对应的水沙过程时间段数据进行比较，并输出比较结果，结束计算过程。

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

式中，F_W、F_U、F_SV、F_IR、F_NI分别为计算单元内水域、不透水域、裸地-植被域、灌溉农田及非灌溉农田的面积率；E_W、E_SV、E_U、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为沟道总长度。

本发明产生的有益效果是：本发明通过收集坡面的水土保持措施的数据以及其他水文、气象数据，利用不同分辨率的栅格DEM数据或等高带数据为平台，以山坡真实水沙物理过程为图景，实现了坡面降雨-径流和侵蚀-泥沙过程的耦合模拟计算能力，在精确水沙过程计算的前提下，对不同时期的坡面水土保持状况进行比较分析，十分精确的获得水土保持措施的效果。本发明由于精确的模拟和计算水沙过程与已有技术相比，对坡面的水土保持状况的分析更加精确，更加接近实际状况，对如何设计水土保持措施提供了更加切合实际的科学依据。

附图说明

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

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

图2是本发明的实施例一所述应用实例中泾河流域降雨-径流量年代变化；

图3是本发明的实施例一所述应用实例中泾河流域降雨-输沙量年代变化；

图4是本发明的实施例一所述应用实例中不同情景下流域单位径流输沙量；

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

具体实施方式

实施例一：

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

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

一、地形数字化处理的步骤：用于DEM数据为基础，以ArcGIS水文模块为主要工具，采用D8法，建立以栅格为基本计算单元的坡面水文过程计算所需的地形数据、汇流关系数据文件，以这些文件为基础，建立坡面侵蚀与泥沙输移过程计算所需的地形数据文件。

（1）以高分辨率栅格DEM数据为基础，以ArcGIS水文模块为主要工具，采用D8法，建立以栅格为基本计算单元的坡面水文过程计算所需的地形数据、汇流关系数据文件。

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

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

二、水土保持措施数据收集的步骤：用于收集模拟区域内不同时间的坡面植被、土壤、土地利用情况，以及水土保持措施情况，形成不同时间段的水土保持数据集。

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

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

数据处理：

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

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

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

四、水沙计算的步骤：用于输入模型率定参数或应用参数，进行坡面水文过程计算、坡面侵蚀与泥沙输移过程计算：

所述的模型率定参数是指在模型建立过程中，首先要对模型是否符合实际请进行率定，因此先要输入一些结果已知的已有参数对模型进行调整和校验，并不断的进行偏差计算，不断的获取模型计算结果与实际结果的偏离度，使模型的计算结果与实际结果的偏离度最小，以此建立十分精确的模型，这些率定模型的已知结果的参数称之为“模型率定参数”。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由于水流在坡面流动具有从上往下并不断扩大流动范围的特点，水沙过程的计算中具有下游的计算单元对该计算单元的上游计算单元，以及上游计算单元相邻的计算单元有数据承传的关系，因此，当一个计算单元完成水沙过程的模拟后，通常需要将水沙过程的计算结果向下游计算单元传递，除非达到坡底的汇流。

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

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

模型参数率定过程是按照模型框架需求，对模型参数进行调整使模型准确反映研究对象真实世界客观规律的过程。模型验证则是对模型参数率定结果的可靠性和准确性检验。验证标准一般有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为模拟均值。

到此为止，水沙耦合模拟模型已经建立完成，如果参数经过率定，就可以进行水体保持计算了。

八、设置水土保持情景的步骤：用于使用水土保持措施数据，设置多个年代的水土保持措施情景。由于水沙耦合模拟模型是基于水沙物理过程建立起来的，其发生和发展过程通过输入数据、计算公式进行了确定了严格的计算关系。水土保持治理通过植树、种草、封育、建设淤地坝等措施改变流域的水沙过程，对应于模拟模型主要表现为土地利用、水利水保工程的变化。因此，通过改变模型用于计算的输入数据和参数的变化，重新运行模拟模型，即可以获得不同水土保持治理条件下的流域水沙过程，通过比较和计算不同的模拟结果，定量评价不同水土保持治理取得的减水减沙效益。

九、生成不同水土保持情景下的水沙过程数据的步骤：用于将设计的水土保持情景数据运行计算，得到多个水土保持措施的水沙过程数据。

十、分析计算不同水土保持情景的保水保土效益的步骤：用于对比不同水土保持情景下的水沙过程变化量进行分析。

十一、比较并结束的步骤：用于对各个情景水土保持情景的水土保持数据集所对应的水沙过程时间段数据进行比较，并输出比较结果，结束计算过程。

汇集指定计算单元或流域出口的水沙数据作为一个水土保持数据集所对应的水沙过程时间段数据，并对多个水沙过程时间段数据进行比较，并输出比较结果，结束计算过程。由于坡面的水土保持随着年代的变化也在不断的变化，如：某一段时间（几年）修建鱼鳞坑，在过几年在鱼鳞坑的基础上栽种树木，在过几年树木被移去而建立梯田等，总之随着时间的推移，坡面在不断的变化，为此，可以按照不同的时间段所收集的坡面水土保持数据，形成这一时段的水土保持数据集，利用这些数据集分别计算各个时段的坡面水沙过程，并对这些不同时段的数据集所对应的坡面水沙过程所形成的水沙过程时间段数据进行比较，以获取不同时段的水土保持措施所取得的效果。

水土保持情景实例：

利用构建好的泾河流域分布式水沙耦合模型，以1956-2000年降雨情景等数据为基础，分析了泾河流域人类活动（下垫面变化）影响下的流域水沙过程变化趋势，以及梯田、坝地等水利水保工程措施减水减沙效益。

情景设置：

为了尝试说明人类活动对流域产水产沙过程的影响，及水土保持工程措施水保效益等进行水土保持治理过程中人们关心的问题，本实例共设定4种不同的人类活动情景：

情景1：历史下垫面情景，即真实下垫面变化情景下模拟1956-2000年泾河流域水沙变化过程，用于还原真实下垫面变化情况下流域水沙变化情况。

情景2：1985年土地利用下垫面情景，即下垫面保持1985年土地利用和水土保持工程措施状态不变条件下1956-2000年模拟泾河流域水沙变化过程。与情景1对比可以说明1985年水土保持措施对流域水沙关系的影响。

情景3：2000年土地利用下垫面情景，即下垫面保持2000年土地利用和水土保持工程措施状态不变条件下1956-2000年模拟泾河流域水沙变化过程。与情景2对比可以说明，1985-2000年间人类活动减水减沙与增水增沙过程的对比关系。

情景4：2000年土地利用-水土保持措施下垫面情景，即下垫面保持2000年土地利用但不含水土保持工程措施状态，模拟泾河流域1956-2000年水沙变化过程。与情景3对比可以说明不同降雨情景下梯田、坝地等工程类水土保持措施的水土保持效益。

使用本实施例建立的模型对人类活动对泾河流域水沙过程演变的影响分析：

研究采用代表时段方法，分别对1956-1959年、1960-1969年、1970-1979年、1980-1989年、1990-2000年、1956-1979年、1980-2000年等7个情景时段的河道水沙量变化进行对比分析。

从表1不同下垫面情景河道径流流量变化可以看出，泾河流域历史下垫面情景（情景1）下1956-1959、1960-1969、1970-1979、1980-1989、1990-2000的流域平均年径流量分别为17.91、24.07、19.96、17.35、8.20亿m³。说明流域径流量与降雨的变化趋势一致，自60年代起流域径流量开始明显减少。根据径流的变化可以以1980年为界，可以明显的分为丰枯两个时段，其中1956-1979的多年平均径流量达21.33亿m³，是1956-2000平均年径流量的1.24倍；而1980-2000年的多年平均径流量为12.56亿m³，仅为1956-2000平均年径流量的72.85%。

与历史下垫面情景（情景1）相比，1985年土地利用下垫面（情景2）1956-1959、1960-1969、1970-1979、1980-1989、1990-2000年对应情景时段的平均年流域径流量分别变化2%、1%、1%、0%、-1%，变化幅度均在±2%以内，说明情景2相对于情景1人类活动对流域的产汇流条件影响较小，但可以明显看出1985年下垫面情况下，增大了1956-1979年的流域横向水通量，说明该时期人类活动造成的增洪量大于减洪量。相对于情景1，情景2的1956-1979年多年平均径流量增加0.27亿m³，而1980-2000年年均径流量减少0.07亿m³，在1956-2000年时段总体上表现为径流量增加；情景3在不同时段的径流量均表现为增加趋势，说明1956-2000年以来人类活动造成的径流增量大于水土保持减水量。对比情景3和情景2可以发现，人类活动在扣除1985-2000年之间新增的水土保持措施减水效益的基础上，在1956－1959、1960－1969、1970－1979、1980－1989、1990－2000时段等不同情况下仍然分别增加0.40、0.76、0.51、0.52、0.25 亿m³。

2000年下垫面（情景3）及2000年下垫面减去工程水保措施（情景4）与情景1相比，在不同时段均表现出径流增加趋势。说明在不同时段，2000下垫面情景现的人类活动增水量大于水土保持措施的减水量。同时对比情景3和情景4可以看出，工程类水土保持措施减水效益明显：1956-1959、1960-1969、1970-1979、1980-1989、1990-2000对应时段分别减水0.76、1.09、1.3、0.74、0.98亿m³，其中由于1956-1979时间段降雨较多，平均减水效益为1.13亿m³，而在1980-2000年时间段为0.87亿m³，总体上年均减水效益为1亿m³。

从图2中可以看出，泾河流域出口的径流量与降雨量的年代变化增减变化趋势一致，说明自然条件对泾河流域径流的影响仍占主导地位。但不同情境下的响应关系不同，其中情景4在1956-1959年到1960-1969年期间的径流量增长率远大于降雨增长率，而从1970-1979年开始，在降雨减少速率较小的情况下径流量迅速减少，说明情景4条件下对缓解气候的变化能力最弱，其次为情景3和情景2。表明梯田、坝地等水土保持工程措施对调节流域降水-产流关系起到了重要作用。

表2可以看出，泾河流域历史下垫面（情景1）下1956－1959、1960－1969、1970－1979、1980－1989、1990－2000年的平均年流域输出沙量分别为2.22、2.53、1.94、2.07、0.62亿t。相对于1956－1959年时段，1960－1969、1970－1979、1980－1989、1990－2000年时段的平均年输出沙量分别变化14%、-13%、-7%、-72%。说明自60年代起流域输出泥沙开始明显减少，但与流域降雨、径流量不同的是，在1980-1989年时段的侵蚀量较1970-1979年时段多出0.13亿t，说明该时期的人类活动的增沙量大于减沙量。以1980年为界，流域在1956-1979年的多年平均输沙量2.23亿t，1980-2000年多年平均输沙量1.31亿t，分别是1956-2000年多年平均的1.24和0.73倍。

1）相比，1985年土地利用下垫面（情景2）1956－1959、1960－1969、1970－1979、1980－1989、1990－2000年对应时段的平均年输沙量分别变化2%、0.3%、1%、1%、-3%，说明在1956-1985年间人类活动的增沙量大于减沙量，在1990-2000年人类活动的增沙量也大于减沙量。这种情况在2000年下垫面（情景3）中体现更为直接，与历史下垫面相比，分别变化5%、6%、5%、5%、3%，说明不同时期的人类活动增加的产沙量大于其减少的产沙量。对比情景3和情景2可以发现，在扣除1985-2000年之间新增的水土保持措减沙效益的基础上，1956－1959、1960－1969、1970－1979、1980－1989、1990－2000年时段等不同情况下人类活动造成的流域增加产沙量分别为0.07、0.13、0.08、0.10、0.04亿t。

对比情景3和情景4可以看出，工程类水土保持措施减沙效益明显：1956－1959、1960－1969、1970－1979、1980－1989、1990－2000年对应时段分别减沙0.23、0.32、0.28、0.16、0.12亿t，其中1956-1979年时段平均减沙效益为0.29亿t，而在1980-2000年时间段为0.14亿t，总体上年均减沙效益为0.22亿t。

从图3，可以看出1956-1959、1960-1969和1970-1979年三个年代间的流域产沙量变化趋势与降雨变化相同。而自1980-1989开始，水土保持措施等人类活动开始对流域减沙产生显著的正面影响。其中情景4的减沙速率明显较大，与情景2对比发现，自1980年以来的非工程水保措施对流域减沙的贡献较大。

从图4可以看出，单位径流输沙量中以情景4最大，而情景3最小。其中情景3在1956-1959、1960-1969、1970-1979年三个年代的单位径流输沙量显著低于其他情景，说明1980年以来的水土保持措施使人类活动显著改变了泾河流域的水沙关系，使流域加速侵蚀的状况得到一定程度的缓解。

实施例二：

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

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

(1)

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

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

(2)

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

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

(3)

E_IR=E_i1+E_i2+E_tr3+Es (4)

E_NI= E_i1+E_i2+E_tr4+Es (5)

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

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

(6)

(7)

(8)

(9)

(10)

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

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

(11)

(12)

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

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

(13)

(14)

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

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

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

(15)

(16)

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

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

(17)

(18)

(19)

(20)

(21)

(22)

(23)

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

实施例三：

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

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

(24)

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

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

(25)

(26)

(27)

(28)

(29)

(30)

(31)

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

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

(32)

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

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

(33)

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

实施例四：

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

地表径流：

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

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

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

(34)

(35)

式中，P为降水量；H_SV为裸地－植生域的洼地储蓄；H_SVMax为最大洼地储蓄深；E_SV为蒸散发；f_SV为由通用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)

式中，H_s为洼地储蓄；H_smax为最大洼地储蓄；Veg₁、Veg₂为裸地－植生域的高植生和低植生的面积率；Rr₁、Rr₂为从高植生和低植生的叶面流向地表面的水量；Q为重力排水；QD_j,j+1为吸引圧引起的j 层与j+1 层土壤间的水分扩散；E₀为洼地储蓄蒸发；E_s为表层土壌蒸发；E_tr为植被蒸散（第一下标中的1表示高植生、2表示低植生；第一下标表示土壤层号）；R₂为壤中径流；k(q)为体积含水率q对应的土壤导水系数；y(q)为体积含水率q对应的土壤吸引力；d为土壤层厚度；W为土壤的蓄水量，w=qd）；W₁₀为表层土壤的初期蓄水量。另外、下标0、1、2、3分别表示洼地储蓄层、表层土壤层、第2土壤层和第3土壤层。

壤中径流：

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

(46)

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

实施例五：

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

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

运动波方程：

（连续方程） (50)

（运动方程） (51)

（Manning公式） (52)

式中，A为流水断面面积；Q为流量；q_L为网格单元或河道的单宽流入量（包含网格内的有效降雨量、来自周边网格及支流的水量）；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方向的分量。

实施例六：

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

（56）

（57）

式中，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之间）。

实施例七：

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

（58）

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

其中雨动能E的计算：

（59）

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

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

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

（60）

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

实施例八：

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

（61）

（62）

（63）

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

实施例九：

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

（64）

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

（65）

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

实施例十：

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

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

（72）

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

重力侵蚀量计算公式为：

（73）

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

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

，n=78，r=0.61 （74）

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

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