CN106599473A - 一种坡面水沙过程耦合模拟方法 - Google Patents

Abstract

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

Description

一种坡面水沙过程耦合模拟方法

技术领域

本发明涉及一种坡面水沙过程耦合模拟方法，是一种水文计算方法，是一种使用计算机技术的水文计算方法。

背景技术

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

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

坡面水土保持治理是开展区域生态环境治理的基本单元。随着经济、社会的发展，区域水土管理越来越精细化和精准化。在水土流失区，坡面水沙过程规律是开展精细化和精准化治理的基础。需要根据山坡水沙过程分异的具体条件，精准配置水土保持治理措施，才能有效改善区域水土条件，取得较好的治理成效。

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

发明内容

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

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

建立数字化坡面的步骤：用于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数据为平台，以山坡真实水沙物理过程为图景，实现了坡面降雨-径流和侵蚀-泥沙过程的耦合模拟计算能力，提升了模拟计算精度，可以有效还原和预测任意格栅尺度的坡面水沙过程。与已有技术相比，本发明以格栅DEM的格栅单元为基本计算单元，实现地形适应性和尺度适应性；同时基于坡面水沙物理过程模拟，能够更加深入的认识坡面水沙过程规律，支撑更多和更复杂的实践应用。

附图说明

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

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

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

具体实施方式

实施例一：

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

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

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

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

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

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

二、建立侵蚀地形数据文件的步骤：用于以地形数据、汇流关系数据文件为基础，建立坡面侵蚀与泥沙输移过程计算所需的侵蚀地形数据文件。对于坡面侵蚀过程，根据汇流条件将计算单元划分为由平面、细沟、浅沟和切沟构成的亚计算单元地貌形态。根据已有的野外调查数据，首先概括出不同侵蚀地貌的“典型侵蚀形态”。

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

式中，N _E 、N _G 分别为计算单元分布的浅沟和切沟数量，k、k＇、α、α＇、分别为浅沟和切沟数量计算的参数和指数，J ₂为计算单元比降，A为计算单元面积。

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

数据收集：收集模拟区域内及其附近的水文、气象、植被、土壤、土地利用、水土保持措施数据，其中水文气象数据包括降雨、平均风速、平均气温、日照时数、空气相对湿度；植被数据包括植被覆盖度、叶面积指数；土壤数据包括土壤类型、土壤厚度；土地利用数据包括模拟时段内多期土地利用数据，一般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为模拟均值。

八、结束的步骤：用汇集坡底各个格栅的数据并输出，结束计算过程。当计算到达坡底的时候，需要将所有坡底格栅的数据汇集，形成坡底的汇流水沙过程。这一数据可以有许多的用途，如流域水沙过程的模拟、进行水土保持分析等。

实施例二：

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

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

(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 ）分别由以下公式计算:

E _SV =E _i1 +E _i2 +E _tr1 +E _tr2 +Es (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 为植物群落阻抗（canopyresistance）。蒸腾属于土壤、植物、大气连续体（SPAC：Soil－Plant－AtmosphereContinuum）水循环过程的一部分，受光合作用、大气湿度、土壤水分等的制约。这些影响通过式（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=q d；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）

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

其中雨动能E的计算：

（59）

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

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

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

（60）

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

实施例八：

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

（61）

（62）

（63）

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

实施例九：

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

（70）

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

（71）

式中，S _E 为股流侵蚀量， Q _E 为流量，T _SE 为股流携沙能力，D _r 为上游来沙量。

实施例十：

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

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

（72）

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

重力侵蚀量计算公式为：

（3-73）

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

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

，n=78，r=0.61 （74）

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

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

Claims (10)

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

建立数字化坡面的步骤：用于DEM数据为基础，以ArcGIS水文模块为主要工具，采用D8法，建立以格栅为基本计算单元的坡面水文过程计算所需的地形数据、汇流关系数据文件；

建立侵蚀地形数据文件的步骤：用于以地形数据、汇流关系数据文件为基础，建立坡面侵蚀与泥沙输移过程计算所需的侵蚀地形数据文件；

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

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

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

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

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

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

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

结束的步骤：用汇集坡底各个格栅的数据并输出，结束计算过程。

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

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

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

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

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

不透水域的地表径流： ；

及

计算，

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

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

壤中径流R2由下式计算：

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

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

运动波方程：

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

动力波方程：

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

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

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

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

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

其中雨动能E的计算：

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

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

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

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

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

式中，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 计算公式为：

式中，k ₆为发生重力侵蚀的沟长系数，L _gully 为沟道总长度，Δh为发生重力侵蚀土体深度，θ为黄土天然休止角。