CN102419788A

CN102419788A - 基于渗蓄一体化动态产流机制的分布式水文模型的设计方法

Info

Publication number: CN102419788A
Application number: CN2011103215374A
Authority: CN
Inventors: 张万昌; 张东
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2010-12-16
Filing date: 2011-10-21
Publication date: 2012-04-18
Also published as: CN102034001A

Abstract

本发明公开了基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，涉及一种涉及水文预报模型设计的方法。该方法的主要步骤为；A)分布式参数获取；B)流域水文过程概化；C)渗蓄一体化产流过程设计；D)汇流过程设计。本发明提出的方法基于渗蓄一体化思想的动态产流机制，协调水量在水文子过程中的分配，动态模拟栅格的产流过程，它不仅可以为其它相关的水文生态研究如流域产、输沙，营养物输移、污染物扩散等提供先进的计算和模拟平台，而且可以为人类减少旱涝灾害、合理开发和利用水资源提供科学的参考依据，为我国的水资源管理和调度提供科学的参考依据和合理的建议。

Description

基于渗蓄一体化动态产流机制的分布式水文模型的设计方法

技术领域

本发明涉及一种涉及水文预报模型设计的方法，更具体地说是一种有效结合遥感(RS)、地理信息系统(GIS)、数字高程模型(DEM)技术建立的具有一定物理机制的通用产汇流方案，并在此基础上构建以栅格为最小模拟单元的、能够反映流域产汇流时空变化机制的分布式生态-水文模型。

背景技术

分布式水文模型作为一种探索水文循环机制、定量评估水资源动态变化的有力工具，在理论上向人们合理解释水文循环的物理基础，并说明水是如何在土壤-植被-大气传输系统中运动并驱动物质和能量的循环。先进技术的应用可以为分布式水文模型带来方法和手段上的更新，但并不能代表模型本身的高效性，需要通过分布式模拟技术加上先进的产汇流机理，才能构建出一个好的实用的分布式水文模型。

从产流机制来看，最早的产流理论是由Horton在1933年提出的著名的超渗产流理论。该理论认为当雨强小于土壤下渗能力时，所有降雨都被土壤吸收；当雨强大于下渗能力时，土壤吸收率等于下渗能力，其余部分形成地表径流。实践证明超渗产流理论适用于干旱、半干旱气候条件下，植被较少或已被破坏，土层较薄，表层土壤易结皮或渗透能力较差的流域的产流计算。与此对应，在湿润的气候条件下，植被茂密，土层比较厚而且渗透能力很大，这时雨强不起作用，所有的降雨都被土壤吸收，当净雨量满足田间持水量后产生饱和坡面流，这就是我们所知的蓄满产流机制。这两种产流过程构成了流域产流过程的两端，被大量的概念性水文模型采用，并通过模型的模拟和应用，验证了这两种产流机制的客观性和合理性，其余的产流理论可以认为是代表了这两端中间不同水源比例的各种情况。

实际上流域的产流是一个动态的过程，具有时空动态变化的特征，即流域的产流不仅在相同时段不同地域单元上客观存在产流类型的差别，而且在相同地域单元不同时段上也同样会存在产流类型之间的差别，因此，产流过程应该是超渗和蓄满产流机制之间的有机结合和动态转换，只不过何时采用何种产流机制，两者之间的比例如何分配，起决定性作用的应该是流域的降雨量特征和下垫面土壤含水量特征的对比，而不是根据气候条件的干旱和湿润与否，在模型建模前决定采用那种产流机制。

超渗产流通过雨强来判断是否产流，需要对降雨过程采用数分钟至半小时的时间步长进行雨量观测，这些数据可用于暴雨过程的模拟，而对日步长以上的长时段水文过程模拟来说，输入的观测数据为日降雨数据，这些数据只能反映降雨量的数值大小，而不能反映雨强，那么如何利用这些日降雨量数据来恰当地估算雨强，从而在长时段水文过程模拟中考虑超渗产流的影响？或者干脆不能计算？同时降雨量观测的时间步长对超渗产流机制的影响到底有多大？

从空间上来讲，湿润地区可通过蓄满产流来模拟，干旱地区可通过超渗产流来模拟，那么，半干旱、半湿润地区用哪种产流机制来模拟？随着气候条件和下垫面条件的变化，处于蓄满和超渗产流过程中间的产流如何计算？两者的比例如何划分？从时间上来讲，在干旱地区，一个以超渗产流为主的流域，遇到长期连绵的低强度降雨，可能发生蓄满产流；在湿润地区，以蓄满产流为主的流域，久旱后遭遇雨强特大的暴雨，也会发生超渗产流；甚至在同一场洪水过程中，也会发生前期为超渗产流，到后期变为蓄满产流的情况。因此流域的产流是一个动态变化的过程，那么如何来表达产流的这种动态变化？

以上产流机制都是基于流域平面来讲的，实际上流域地形可能变化很大，起伏明显，在坡脚以及凹坡等流线辐合的地方容易形成局部饱和坡面流，而且该面积随着时间的变化而改变，这些和流域平面不同的特性主要是由于流域地形影响而造成的，如何考虑地形对产流过程的影响，使产流更加符合实际情况？

从汇流机制来看，概念性水文模型如何与DEM技术紧密结合，实现流域分布式水文过程模拟，不仅得到流域出口的流量模拟结果，而且得到流域内各模拟单元的水文过程分量，值得深入研究。从目前大多数模型的汇流方案来看，不同的模型有的采用简单的单位线法、线性水库法、马斯京根法，有的采用较复杂的运动波、扩散波和动力波模型，将流域的出流汇集到流域出口。不同的汇流方案有各自的适应性，会产生不同的模拟结果，对出口流量的过程描述差别也较大。对不同的模拟流域，该选择何种合适的汇流方案来提高模型的模拟精度，应该给出具体的方案。

发明内容

1.要解决的问题

本发明的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，针对现有技术的不足，有效结合遥感(RS)、地理信息系统(GIS)、数字高程模型(DEM)技术，建立一种具有一定物理机制的通用产汇流方案，并在此基础上构建以栅格为最小模拟单元的、能够反映流域产汇流时空变化机制的分布式生态-水文模型，该模型能够广泛适应不同气候条件、不同下垫面条件下的流域生态-水文过程模拟，在不同的时空尺度上(时间尺度：分钟/小时/日/月/年，空间尺度：30m～1000m)真实模拟和刻画流域生态水文循环过程。

2.技术方案

本发明的目的是通过以下技术方案来实现的：

基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，主要包括以下步骤：

A)分布式参数获取；

B)流域水文过程概化；

C)渗蓄一体化产流过程设计；

D)汇流过程设计。

所述步骤A)中，获取包括气象、土壤、植被和水系分布式参数：

(1)气象参数缺测条件下的获取

本发明以数理统计、随机模拟和弱稳定理论为基础，使用WGEN气象发生器，实现降雨量、气温、太阳辐射、相对湿度和平均风速五类气象参数的随机解集计算，用于站点气象资料不全时的填补、气象数据缺测时的模拟以及未来气象参数的预测。

a.气象参数获取

气象发生器的模拟原理为数理统计、随机模拟和弱稳定过程理论，通过气象资料月平均数据(如月平均值、标准差等)的统计规律加上随机噪声矩阵来模拟周期性的气象变化过程。

通过求解数据变量的月平均值、标准差以及变量的残差元序列来模拟气温和太阳辐射数据的随机解集。相对湿度定义为空气中实际水汽压与同一温度下饱和水汽压的比值。根据风速频率直方图的分布形态，采用指数方程结合随机模拟来实现平均风速的随机解集求解。降雨量随机解集采用一阶马尔科夫链-偏斜正态分布降水模型来计算。

b.气象参数空间离散

使用邻近站点法、反距离加权法、泰森多边形法、复杂的梯度距离平方反比法或基于DEM的PRISM模型法用于模型的参数空间离散，通过GIS空间离散，得到每一个栅格的气象参数，以便进行栅格水文过程模拟。

(2)土壤属性参数库的建立

对于每一种土壤类型，需要外部输入的土壤参数主要分为两类，一类是土壤分层参数，描述具体每一分层的土壤水力物理特性；另外一类是土壤剖面参数，描述该土壤类型下与本栅格所有土层有关的土壤剖面信息。

土壤分层参数：土壤深度、土壤湿容重、有效持水量、饱和导水率、稳定下渗率系数、初始下渗率系数、有机物含量和土壤粒径组成。

土壤剖面参数：土壤类型名称、土壤分层数、植被根区深度、地表反照率、霍顿下渗曲线系数和初始土壤水比例。

对不同的土壤类型建立包含土壤分层参数和剖面参数的土壤属性参数库。

(3)植被属性参数库的建立

植被参数分为土地利用/植被覆被参数和植被生态参数两类。植被类型参数的获取，是由国家或地方林业部门提供流域的土地利用/覆被图，通过GIS数字化以及矢量～栅格格式转换的方式，得到栅格植被类型分布图；植被生态参数是利用遥感反演结合实地调查验证的方法获得的；

土地利用/植被覆被参数：常绿针叶林、常绿阔叶林、落叶针叶林、落叶阔叶林、针阔混交林、林地、林质草地、封闭灌丛、开放灌丛、草地、农田、裸露地和城镇居民地；

植被生态参数：最大叶面积指数、最小叶面积指数、莱峰偏移、地表噪度系数、月平均叶面积指数、月平均反射率和月平均覆盖度。

对不同的植被类型建立包含土地利用/植被覆被参数和生态参数的植被属性参数库，将每一类植被类型代码和相应的植被生态参数之间建立关联，在栅格水文模拟时，通过栅格上的植被类型读取和推算该模拟时段的植被生态参数，用于栅格能量和水量平衡计算；

(4)基于DEM的流域水系参数提取

以D8算法即最大落差算法(芮孝芳，水文学原理，中国水利水电出版社，ISBN：7508421647)为基础，采用基于图论的广度优先搜索算法实现流域所有栅格的遍历，实现模型所需的水系参数的自动提取。基于DEM的流域水系参数提取的步骤包括：水流流向矩阵计算、水流累积矩阵计算、河网自然水系提取和流域边界计算。

a.水流流向矩阵计算

原始的DEM首先通过填洼处理，得到填洼后的DEM，在此基础上进行水系和汇流参数提取。

b.水流累积矩阵计算

在水流方向矩阵的基础上，通过对任意一个栅格利用基于图论的广度优先搜索算法遍历所有与本栅格有入流关系的上游栅格，可以建立水流累积矩阵。

c.河网自然水系提取

通过给定合适的上游给水区面积阈值，将水流累积量矩阵数据中高于此阈值的格网连接起来，在此基础上搜寻出水系的所有栅格点，形成完整的自然水系网络。

d.流域边界计算

设定流域出口断面所在的DEM数据的行列号，通过广度优先搜索算法搜索对流域出口有贡献的栅格，赋为1，其余赋为0，可以自动勾画出流域边界。

所述步骤B)中，流域的陆面水文过程涉及到水量在大气-植被-土壤之间的循环和分配，是一个复杂的时变非线性系统，在建立分布式参数获取方法基础上，建立流域水文过程概化模型。

在不同的空间地域单元内，由于下垫面的空间变异性，水文循环过程差异很大，但是对于同一个地域单元，可以是子流域、高程带、水文响应单元或者栅格，如果空间尺度足够精细，可以认为单元内部下垫面特征是均一的，在相同的植被和土壤类型下，有相同的水文响应特征。因此，对于概念性分布式水文模型，模型的产流过程是在地域单元上实现的。在垂直方向，大气水以降雨或降雪的方式降落，首先接触到植被，植被冠层对雨水进行截留以满足自身的需水和蒸发、蒸腾需要；多余的水量降落到地表，产生地表水下渗，下渗到土壤中的水在土壤和植被根系的作用下在土壤中进行再分配处理，在合适的情况下产生向上运动的蒸散发和侧向流动的地表径流、壤中径流和地下径流，完成单元的产流过程；单元产流产生的径流分量在水头的作用下，沿坡面方向产生运动，完成单元间的水量交换，通过河道汇流运动到流域的出口。对陆相水文循环过程进行简化，将耦合的水文子过程线性化处理，分解为相对独立的栅格单元水文子过程，以能量平衡和水量平衡为基础，建立概念性分布式流域水文模型。

其中：冠层截留模拟采用改进的概念性Aston指数模型，精度较高。融雪径流计算采用常用的能量平衡法或气温指标法。潜在蒸散发计算采用Penman-Monteith(PM)法、Pristley-Taylor(PT)法或Hargreaves法。地表水下渗计算包括霍顿模型(Horton R.E.)、科斯加可夫模型(Kosjiakov，A.N.)、霍尔坦模型(Holton)、菲利普模型(Philip)或格林-安普特模型(Green-Ampt)。

根据土壤分层水量调节参数土层孔隙率、田间持水量时植被可吸收水量、土层每日实际可供水量、土壤水饱和时土层可供水量、田间持水量时土层含水量、枯萎点时土层含水量和土壤水力传导率系数，动态计算每一分层对入渗水的调节状态。根据土壤剖面水量调节参数土壤剖面平均湿容重、土壤剖面平均孔隙率、土壤剖面平均饱和导水率、土壤剖面稳定下渗率、土壤实际含水量、土壤饱和时土壤含水量、田间持水量时土壤含水量和枯萎点时土壤含水量，进行土壤剖面水量调节计算。

在步骤C)中，除了常规的栅格水量在垂直方向上完整的输入和输出，还考虑了栅格水平方向的出流，本发明同时顾及栅格垂直方向和水平方向出流，设计了一个完整的渗蓄一体化栅格产流方案，具体如下：

(1)动态产流方案设计，如表1

表1：栅格通用产流类型表

根据霍顿下渗理论、蓄满产流理论和山坡水文学原理，按照前期土壤水蓄满与否，建立完整的栅格通用产流类型；

表中：D-土壤缺水量(mm)；P-到达地面的净雨(mm)；E-土壤实际蒸散发(mm)；I-雨强(mm/hr)；FC-土壤表层稳定下渗率(mm/hr)；FP-土壤表层时段下渗率(mm/hr)。

按前期土壤水是否蓄满进行预分类的目的在于土壤下渗率是个时变的值，如果前期土壤水已经蓄满，则在计算时段内采用土壤稳定下渗率来和雨强作对比，否则采用土壤时段下渗率来和雨强作对比，以决定是否发生超渗产流。

(2)栅格产流计算(计算流程如图3)

取栅格为水文模拟单元，在高空间分辨率条件下，近似认为栅格内的土壤、植被特征是均一的，针对表1中不同的产流方案和水量平衡理论，按不同的前期土壤水状态，可以分别计算出各部分的水量分配情况。

计算过程中的变量定义如下：R_T-总产流；R_S1-超渗地表径流；R_S2-蓄满地表径流；R_L-壤中径流；R_G-地下径流；PE＝P-E；IF_C＝I-F_C；IF_P＝I-F_P；Δw-土壤含水量的变化量。

I.地表径流计算

●产流类型①——不产流：

在产流类型①的条件下，雨强小于土壤稳定下渗率，同时降雨量不够土壤和植被蒸发，不足部分水量从土壤前期含水量中补充，因此，不会发生产流事件，土壤含水量变化Δw为负，各部分产流关系如下：

R_T＝0，Δw＝PE＜0 (1)

R_S1＝R_S2＝0 (2)

R_L＝R_G＝0 (3)

●产流类型②——超渗产流：

产流类型②的雨强大于土壤稳定下渗率，满足超渗产流条件，所以虽然总的降雨量不够土壤和植被蒸发，但是净雨首先产生超渗地表径流，然后不够蒸发的部分再从土壤水中补充。各部分产流关系如下：

R_T＝R_S，Δw＝PE-R_T (4)

R_S1＝IF_C·Δt (5)

R_S2＝0 (6)

R_L＝R_G＝0 (7)

●产流类型③——蓄满产流：

在雨强小于土壤稳定下渗率的条件下，前期土壤水蓄满后，外部供水除了满足土壤和植被蒸发，还有多余的水分下渗，满足蓄满产流的条件，此时的产流为蓄满产流。由于土壤水始终处于饱和状态，所以土壤含水量变化Δw为0。各部分产流关系如下：

R_T＝PE，Δw＝0 (8)

R_S1＝0 (9)

R_G+R_L＝F_C·Δt (10)

R_S2＝R_T-(R_G+R_L) (11)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf (12)

式中：d_Surf为地表径流R_S2与地表以下径流R_G+R_L的分流比(0≤d_Surf≤1)。设置该变量的目的在于如果Δt取得不合适而导致R_S2小于0时，可以通过d_Surf对R_S2进行调整：如果模拟流域地表径流影响比较大，则d_Surf取较大值；如果模拟流域的径流补给以地下水为主，则d_Surf取较小值。此时，壤中流和地下径流的关系调整如下：

R_G+R_L＝R_T·(1-d_Surf) (13)

●产流类型④——超渗产流+蓄满产流：

前期土壤水蓄满后，雨强大于土壤稳定下渗率，同时外部供水能够满足土壤和植被的蒸发以及土壤水的下渗，此时，将会发生超渗产流和蓄满产流同时存咱的状况。在这种情况下，在控制总水量平衡时，首先满足地表产流，然后计算地下径流和壤中流。各部分产流关系如下：

R_T＝PE，Δw＝0 (14)

R_S1＝IF_C·Δt (15)

R_G+R_L＝F_C·Δt (16)

如果R_S1+R_G+R_L＞R_T，则有：

R_G+R_L＝R_T-R_S1 (17)

R_S2＝0 (18)

否则，有：

R_S2＝R_T-R_S1-F_C·Δt (19)

在式4-74中，如果R_G+R_L＜0，则取

R_S1＝R_T·d_Surf (20)

R_G+R_L＝R_T·(1-d_Surf) (21)

●产流类型⑤——不产流：

在产流类型⑤的条件下，雨强小于土壤时段下渗率，降雨量不够土壤和植被蒸发，不足部分水量从土壤前期含水量中补充，因此，在该条件下，不产生地表径流。由于前期土壤水未蓄满，所以计算不足部分水量的时候，需要计算该水量是否超过土壤前期含水量与枯萎点水量之差。如果超过，则取此时土壤能够提供的自由水量来代替不足水量。各部分产流关系如下：

R_T＝0，Δw＝PE (22)

R_S1＝R_S2＝0 (23)

R_L＝R_G＝0 (24)

如果PE＜SP_Sw-SP_Wp，则：

Δw＝SP_Sw-SP_Wp (25)

●产流类型⑥——超渗产流：

产流类型⑥和产流类型②类似，不同之处在于产流类型⑥的前期土壤水没有蓄满，所以超渗地表径流用时段下渗率来计算，而产流类型②用稳定下渗率来计算。各部分产流关系如下：

R_T＝R_S，Δw＝PE-R_T(26)

R_S1＝IF_P·Δt (27)

R_S2＝0 (28)

R_L＝R_G＝0 (29)

●产流类型⑦——蓄满产流：

和产流类型③相比，由于前期土壤水未蓄满，所以如果产生蓄满产流，必须首先补充土壤包气带缺水量，因此，土壤含水量变化Δw为土壤缺水量D，计算如下：

D = (1 - \frac{SP_Sw}{SP_Fc}) \cdot SP_Por \cdot Σ_{i = 1}^{n} SL_Z (i) - - - (30)

当外部供水大于土壤缺水量D时，剩余水量用于产流计算，各部分产流关系如下：

R_T＝PE-D，Δw＝D (31)

R_S1＝0 (32)

R_G+R_L＝F_P·Δt (33)

R_S2＝R_T-(R_G+R_L) (34)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf (35)

此时，壤中流和地下径流的关系调整如下：

R_G+R_L＝R_T·(1-d_Surf) (36)

●产流类型⑧——超渗产流+蓄满产流：

在产流类型⑧中，外部供水满足了土壤缺水量D，有多余的水量用于产生蓄满地表径流；同时雨强大于土壤地表时段下渗率，满足超渗产流的条件，因此，超渗产流和蓄满产流同时存在。各部分产流关系如下：

R_T＝PE-D，Δw＝D (37)

R_S1＝IF_P·Δt (38)

R_G+R_L＝F_P·Δt (39)

如果R_S1+R_G+R_L＞R_T，则有：

R_G+R_L＝R_T-R_S1 (40)

R_S2＝0 (41)

否则，有：

R_S2＝R_T-R_S1-F_P·Δt (42)

如果R_G+R_L＜0或者R_S1＞R_T，则对径流分量调节如下：

R_S1＝R_T·d_Surf (43)

R_G+R_L＝R_T·(1-d_Surf) (44)

II.壤中径流和地下径流计算

已知地表径流和地表以下径流的条件下，土壤剖面含水量计算如下：

SP_Sw_i＝SP_Sw_i-1+PE_i-R_S1，i-R_S2，i-(R_G，i+R_L，i)(45)

其中，下标i和i-1分别表示当天和前一天的水量值(mm)。通过对时间步长的循环计算，可以得到每一个计算时段的土壤含水量状况。

在地表径流计算过程中，计算得到的是壤中径流和地下径流的总和，通过对土壤水量下渗的处理，来划分两者的比例关系。

土壤水下渗计算分为两种情况，一种是土壤含水量超过田间持水量时的下渗，另外一种是土壤含水量在枯萎点水量和田间持水量之间时的正常下渗。

当土壤含水量大于田间持水量时，计算超过田间持水量部分的土壤水下渗W_P如下，下渗水量随时间呈指数趋势衰减：

W_{P} = (SP_Sw - SP_Fc) \cdot [1 - \exp (- \frac{Δt}{Δ t_{P}})] - - - (46)

式中：Δt为计算时段长度，根据模拟的时间步长具体确定；Δt_P为土壤水下渗时间，它是土壤饱和含水量、土壤田间持水量状态时的含水量以及土壤饱和导水率的函数，计算如下：

Δ t_{p} = \frac{(SP_Sat - SP_Fc)}{SP_Sat_K} - - - (47)

此时的壤中径流定义为下渗水量的函数，有

R_L＝W_P·d_L (48)

式中：d_L为外部率定的壤中径流出流因子，一般取0≤d_L≤1。

当土壤含水量在枯萎点水量和田间持水量之间时，首先根据时段下渗率计算土壤下渗水量W′_P如下：

W′_P＝F_P·Δt (49)

因为Δt取值的不确定性，有可能造成计算出的W′_P大于土壤可供的最大自由水量的情况，此时，取W′_P为土壤自由水量，即

W′_P＝SP_Sw-SP_Wp(50)

在这种状况下，由于土壤水不可能在Δt时段内完全按照时段下渗率来下渗，因此，考虑根据当前土壤含水量值按指数关系衰减，计算实际的土壤水下渗量W_P如下：

W_{P} = W_{P}^{'} \cdot [1 - \exp (- \frac{W_{P}^{'}}{SP_Sw})] - - - (51)

壤中径流量同样由上式计算得到。这样，当计算出壤中径流后，由壤中径流和地下径流之和减去壤中径流，就可以分离出地下径流量。

至此，在渗蓄一体化思想下，构建出栅格地表产流的动态判别机制，计算地表产流，然后根据水量平衡理论决定壤中流和地下径流量，并根据土壤水的指数退水规律进行水源划分，建立通用的产流方案，计算出栅格的地表径流、壤中径流和地下径流流量，完成栅格的产流计算。

在栅格上计算出各径流分量后，通过坡地和河道汇流处理，可以将水量汇集到流域的出口，从而模拟得到流域的流量过程。在步骤D)中，针对次洪过程和长时段降雨～径流过程的不同特征以及坡地汇流和河道汇流的差别，分别设计了马斯京根-康吉法、滞时演算法和河道分段马斯京根法，用于不同条件下模型的汇流过程演算。

(3)汇流计算方法

a.运动波汇流法

运用马斯京根-康吉法进行全栅格运动波汇流模拟。

b.滞时演算法

当栅格地表径流、壤中径流和地下径流的汇流过程不满足洪水演进条件时，采用滞时演算法来进行汇流处理，通过考虑不同径流类型的沿程传输损失，获得流域出口的流量过程。

线性法和指数法计算得到的传输损失系数非常接近，基本上可以相互替代。但是从两种方法的取值来看，由于指数法的参数数值很小，因此汇流计算结果对参数的敏感性比较强，在参数选择上不如线性法来得稳定，所以在本模型中采用线性法来计算汇流的沿程损失。

c.分段马斯京根法

当流域比较大时，对全流域采用全栅格汇流方式将所有栅格的产流量统一汇集到流域出口是不合适的，可以采用划分子流域的方法进行流域汇流演算。

河道的马斯京根汇流演算分为先演后合法和先合后演法两种。先合后演法更符合实际的流量过程。在利用先合后演法进行河道马斯京根分段连续演算时，没有采用等距分段的原则，而是根据自然河段的长度，分别设置每一河段的时间差分权重系数和洪水波在河段中的传播时间，以体现不同河段之间的差别。

(4)汇流模拟策略

本模型的汇流策略采用模块搭配的方式来解决。首先模型提供基本的汇流演算方法，有基于栅格水系的马斯京根-康吉运动波模拟法、滞时演算法、马斯京根法(先演后合、先合后演)或基于河道分段的马斯京根分段汇流演算法；然后根据径流类型和模拟的时段长度来选择不同的汇流演算算法，得到适合的汇流方案组合；最后根据所选用的汇流方案组合，模拟得到流域的汇流过程。

本发明提出了基于渗蓄一体化思想的动态产流机制，摒弃了传统的产流方案对流域气候条件的依赖，通过实时判断降雨特征与下垫面土壤水变化特征的对比，动态选择水文模拟单元的产流类型，以此实现产流机制的气候普适性，改变了目前水文模型不能同时适应干旱区和湿润区水文过程模拟的状态。

3.有益效果

本发明建立一个先进的概念性分布式流域水文模型。

(1)充分耦合先进的RS、GIS和DEM技术，结合地表观测数据，获取模型所需的水文、气象和流域特征参数；

(2)基于渗蓄产流一体化思想，构建通用的渗蓄一体化动态产流模型，真实反映产流过程的时空变化，使模型具有良好的气候、下垫面条件适应性，能够满足不同流域的水文过程模拟；

(3)提供基于栅格和河道的多种汇流演算策略，不仅可以得到流域出口断面的流量过程，而且可以得到任意模拟时段、任意空间位置的水文子过程的空间分布；

(4)可以模拟短期的暴雨径流过程和长期的降雨径流过程，定量分析和预测不同时间尺度下流域的水文水资源动态变化；

本发明建立的分布式流域水文模型不仅可以为其它相关的水文生态研究如流域产、输沙，营养物输移、污染物扩散等提供先进的计算和模拟平台，而且可以为人类减少旱涝灾害、合理开发和利用水资源提供科学的参考依据，为我国的水资源管理和调度提供科学的参考依据和合理的建议。

说明书附图

图1为本发明的模型体系架构图；

图2为栅格水文过程概化图；

图3为栅格通用产流计算流程图；

图4为江口流域次洪模拟结果与实测结果对比图；

图5为江口流域日降雨径流模拟结果与实测结果对比图(1981～1985)；

图6为江口流域月、年平均径流模拟结果与实测结果对比图(1981～1985)。

具体实施方式

本发明的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，模型体系架构图如图1，栅格水文过程概化如图2，主要步骤为：

A)分布式参数获取；

B)流域水文过程概化；

C)渗蓄一体化产流过程设计；

D)汇流过程设计。

以下通过实施例结合附图对本发明作进一步描述：

实施例一

温润区水文过程模拟与精度评价。

1.研究区概况

本研究区选择在汉江上游区的左侧支流——褒河上游的江口流域，流域的集水面积约2413km²，经度106°48′15″～107°25′34″，纬度33°38′03″～34°11′08″，海拔高程约900～3400m，流量控制站点为江口站。

江口流域的气候属于亚热带季风气候，具有明显的垂直地带性，全年降雨多在750～1000mm之间，降雨量年内分配不均，主要集中在6～9月。流域内植被覆盖度高，自然植被保护良好。流域的土壤类型以黄棕壤、棕壤和淋溶褐土为主，土壤分布的空间差异较大。

2.研究区基础数据

A.水文、气象数据

江口流域有雨量站11个，气象站3个。雨量站提供有1981～1985年的实测日降雨量资料，气象站有同期的气温、太阳辐射、相对湿度和平均风速资料，其中，江口站有同期的日蒸散发资料。流域的流量验证站点为流域出口——江口站。

为了模拟和验证暴雨过程，另外摘录了1980-1985年期间的6场典型洪水过程的降雨、蒸发和流量数据。次洪过程模拟的时间步长取为1小时，因此，摘录的暴雨资料通过时间平均，均匀分配到每一个整点时间段，对应的流量过程通过流量资料线性插值方式获得，计算得到的日潜在蒸散发和实际蒸散发则通过简单正弦关系离散到整点时段。

B.DEM数据

江口流域的DEM数据由1∶25万比例尺的地形图数字化生成，空间分辨率为60m，采用UTM投影，WGS84坐标系，第48分带。

C.土壤数据

江口流域的土壤图由1∶25万比例尺的陕西省土壤质地图数字化生成，为了和DEM相匹配，采用与DEM相同的空间分辨率、投影类型和坐标系。土壤属性数据来自《陕西土壤》(郭兆元，1992)，结合模型的土壤粒径转换程序包整理得到。如表2。

表2：汉江流域(温润区)土壤参数类型表

暗棕壤	黄棕壤	棕壤性土	红土性土	石灰性褐土
					潮土	黄棕壤性土	棕壤	粘化黑卢土	紫色土
黄褐土	淋溶褐土	新积土	褐色土	漂洗棕褐

[0203]

黑卢土

嵝土

D.植被数据

流域的植被类型图由2001年7月28日的Landsat+ETM(Pass36，Raw128)遥感影像经过监督分类解译得到，采用野外调查的土地利用数据以及1996年发布的1km精度的全国土地利用/覆盖图作为影像监督分类时选取训练区的参考和分类结果的评价，最终得到的植被类型为6类：林质草地、落叶阔叶林、封闭灌丛、常绿阔叶林、针阔混交林、落叶针叶林。

3.次洪过程模拟和精度评价

江口流域共摘录了1980～1985年间的6次典型洪水过程，以日期为标识，洪号分别为19800628次、19810810次、19830925次、19831003次、19840709次和19850911次，模拟的时段长度从流域降雨开始，到洪峰过程退水结束为止。

对于江口流域的次洪过程模拟，模型采用的运行方案如下表3所示：

表3：江口流域次洪过程模型运行方案表

按照渗蓄一体化动态产流机制，模拟了江口流域日降雨～径流过程的栅格产流类型，通过观察产流类型的空间分布和随时间的变化规律，来验证渗蓄一体化动态产流思想的正确性。模拟的降雨～径流过程时段为从1981年1月1日到1985年12月31日，站点降雨数据通过GIDW法离散到每一个栅格，在对日降雨量和同期暴雨过程的雨强资料对比分析的基础上，率定出日降雨量的平均降雨时长，推算出日降雨的雨强。在整个模拟时段中截取了1981年8月14日至8月24日共11天的模拟结果进行分析，选取该时段的原因在于在此期间有一个典型的强降雨过程，流量过程变化明显，可以清楚地分析流域产流类型的空间分布和时序变化。

综合降雨～径流过程图和栅格产流类型图，可以看出栅格产流类型的变化具有以下基本规律：

1)在土壤前期含水量比较低的情况下，流量过程相对于降雨过程有一个滞后性，显示出流域下垫面对降雨水量的调蓄作用，随着土壤含水量的增加，该滞后性逐渐减短；

2)产流范围的分布与流域的下垫面特征密切相关，由于土壤类型和植被类型分别通过土壤类型图数字化和遥感反演得到，因此图中的产流区域明显具有与土壤类型图对应的连续斑块特征和与遥感影像对应的像元斑块特征；

3)产流总是在河道附近容易发生，而且对于不同Strahler等级的河道来说，高等级的下游河道比低等级的上游河道更容易产流。这一方面是因为河道附近地势低，地下水位相对较高，土壤含水量较大；另外一方面降雨以后水分向河道方向汇集，造成河道附近区域先产流；

4)超渗产流客观上来讲需要大的雨强，而降雨的空间离散采用考虑高程影响的反距离加权GIDW法，所以超渗产流的区域往往容易出现以降雨站点位置为中心向外扩展的圆形形状；

5)从产流类型的空间分布来看，根据降雨和下垫面特征空间分布形态的不同，任一时刻产流类型可以是8种产流类型的全部或者其中的一部分；

6)从产流类型的时序变化来看，如果有持续的降雨补充，前期土壤水不足条件下的四种产流状态开始向土壤水饱和条件下的四种产流状态转变，同时上个时段不产流时，下个时段容易变成蓄满产流，而上个时段超渗产流时，下个时段容易变成不产流或者超渗产流和蓄满产流同时发生；如果没有足够的降雨补充，则产流状态发生反向转变，原先产流的区域可能变成不产流，土壤含水量重新变成不饱和，直到下一个降雨～径流循环过程，重新开始上述的变化状态。

所以可以看出，基于渗蓄一体化动态产流机制模拟(公式51)得到的栅格时序产流类型基本上真实地反映了产流的动态变化，可以应用于分布式产流过程模拟。

汇流方案详细演算方法：

A.马斯京根-康吉法(运动波)

该方法采用四点带权偏心差分格式，空间差分权重取0.5，并考虑了区间侧向入流的影响，计算公式如下所示：

Q_{j + 1}^{n + 1} = C_{1} \cdot Q_{j}^{n} + C_{2} \cdot Q_{j}^{n + 1} + C_{3} \cdot Q_{j + 1}^{n} + C_{4} \cdot q_{n + 1}

其中：

C_{1} = \frac{K \cdot X + 0.5 \cdot Δt}{K \cdot (1 - X) + 0.5 \cdot Δt}

C_{2} = \frac{0.5 \cdot Δt - K \cdot X}{K \cdot (1 - X) + 0.5 \cdot Δt}

C_{3} = \frac{K \cdot (1 - X) - 0.5 \cdot Δt}{K \cdot (1 - X) + 0.5 \cdot Δt}

C_{4} = \frac{Δt \cdot Δl}{K \cdot (1 - X) + 0.5 \cdot Δt}

j为空间步长下标，n为时间步长下标。因此，

和

分别为栅格前一时段的总入流(m³/s)、栅格当前时段的总入流(m³/s)、栅格前一时段的出流(m³/s)和栅格当前时段的出流(m³/s)；q_n+1为栅格区间当前时段的侧向单宽入流(m³/m)；X为时间差分权重系数；K为洪水波在河段长为Δl中的传播时间，对于栅格来讲为坡面径流流过栅格长度或栅格对角线长度的时间(min)。

由于该演算格式为显式有限差分格式，由Von Neumann稳定性分析方法可知，其稳定性条件为

X≤0.5

将基于一维河道洪水波模拟的马斯京根-康吉法运用于栅格水系，需要将二维的栅格矩阵进行降维处理，变成一维栅格矩阵，在此基础上才能利用上述的四点偏心加权差分离散关系，在时间维t和空间维x上对流量进行连续演算。

降维的方法是建立流域内栅格上下游汇流次序表和栅格汇流最优演算次序矩阵。通过建立栅格上下游汇流次序表，可以知道与当前栅格相邻的8个栅格中那些上游栅格会有水量流入本栅格，同时当前栅格又流向哪一个下游相邻栅格；通过建立栅格汇流最优演算次序矩阵，可以知道流域内不同水系所属的栅格之间的先后演算次序。

栅格汇流最优演算次序矩阵可以从DEM上导出。按照DEM高程由高到低的方式，栅格之间的水流从流域的边界逐渐向下游汇集，一直到达流域的出口节点，因此流域出口节点具有最大的汇流演算次序，依次向上游递减，直到流域的边界。这样，汇流时从最小汇流演算次序的栅格向下游汇流，一级一级向前推进，相当于模拟出洪水向前演进的过程。

当提取出运动波汇流参数矩阵后，进行降维处理。将二维的栅格汇流通过栅格之间的高差关系，汇集或合并成一维的栅格入流来处理，达到栅格水系汇流的降维目的。在此基础上，由马斯京根-康吉法计算出的出流，作为另外栅格的入流之一。这样，通过栅格时序演算和不同栅格之间流量的推演，就可以得到任意栅格任意时段的流量值，从而实现栅格汇流的运动波模拟。

B.滞时演算法

栅格水系滞时演算法汇流的关键在于确定栅格的汇流时间。栅格汇流时间定义为栅格产流流量汇集到(子)流域出口所需要的时间，根据定义，每个栅格的汇流平移时间是栅格的水流路径长度和水流路径平均坡度的函数，两者可以通过流域DEM来计算得到。

栅格水流路径长度指流域内每一个对流域出口有贡献的栅格按照流向汇集到出口时流过的栅格长度的总和。根据流向不同，在一个栅格上，流过的长度有栅格宽度或栅格对角线长度的区别。在具体的实际操作中，计算步骤和方法如下：

a)从流域出口栅格向上游按广度优先算法遍历所有对出口有流量贡献的栅格，压入队列；

b)按队列的先后顺序计算栅格的水流路径长度。首先按照栅格水流流向，计算本栅格的水流流经长度；然后按照队列的先后顺序，后面的栅格(上游栅格)的水流路径长度为栅格本身的水流流经长度与下游栅格水流路径长度之和；

c)遍历完队列中的所有栅格，完成流域的栅格水流路径长度计算。

依照栅格水流路径长度矩阵的计算方法，可以方便地计算出栅格水流路径平均坡度矩阵。方法有以下两种：

a)取栅格水流路径上从流域出口栅格到本栅格的各栅格坡度之和的平均值，作为本栅格的水流路径平均坡度；

b)通过栅格水流路径首尾栅格高程差除以栅格水流路径长度，得到本栅格的水流路径平均坡度。

经计算检验，两种方法得到的栅格水流路径平均坡度值结果一致。

当已知了栅格的水流路径长度和水流路径平均坡度后，栅格汇流时间可以通过下式来计算：

{RT}_{i, j} = \frac{L_{i, j}}{V_{i, j}} = \frac{L_{i, j}}{dKv \cdot \sqrt{S_{i, j}}}

式中：i和j为栅格的行列号，RT为栅格平均汇流时间(hr)；L为栅格水流路径长度(m)；V为栅格水流平均流速；S为栅格水流路径平均坡度。由于V随降雨量、降雨强度、土壤下垫面性质以及时间而变，难以确定，因此利用

来代替计算(王兆礼，陈晓宏，杨涛；东江流域径流系数变化特征及影响因素分析；自然资源学报；《自然资源学报》；2010年第08期)。dKv为速度参数，包含了糙率、水力半径等因素对水流的影响，运算时需要通过率定得到。下表4显示的是江口流域次洪过程模拟的速度参数率定值：

表4：江口流域栅格汇流时间参数表

径流类型	地表径流	壤中径流	地下径流
				dKv	4000	3500	3000

根据率定得到的速度参数，分别计算出江口流域次洪过程模拟中的地表径流、壤中径流和地下径流的栅格汇流时间。

在水源沿汇流路径到达(子)流域出口的途中会产生沿程的传输损失，包括沿程的下渗和蒸散发。传输损失按以下公式计算：

Q_{out} = \frac{Q_{i, j}}{K \cdot {RT}_{i, j}^{n}}

式中：Q_out为(子)流域出口处的汇流流量(m³/s)；Q_i，j为分水源后的栅格径流流量(m³)；K为汇流传输损失系数，通过率定得到；n为汇流时间的指数，指数越大，沿程损失越厉害。通常取n＝1，表示汇流流量与汇流时间成倒数关系。下表5列出了江口流域所率定的线性法的汇流传输损失系数：

表5：江口流域线性法汇流传输损失系数表

径流类型	地表径流	壤中径流	地下径流
				K	60	400	1600

不同水源的沿程传输损失还可以用指数关系来计算：

Q_{out} = Q_{i, j} \cdot [1 - \exp (\frac{- K}{{RT}_{i, j}})]

式中的变量含义同上，参数率定如下表6：

表6：江口流域指数法汇流传输损失系数表

径流类型	地表径流	壤中径流	地下径流
				K	0.017	0.0025	0.0006

C.分段马斯京根法

河道的马斯京根汇流演算分为先演后合法和先合后演法两种。先演后合法假定各子流域出口的流量在主河道中满足线性关系，因此通过分段马斯京根法将各子流域出口的流量分别演算到流域出口，然后叠加得到流域出口的总流量。先合后演法则按照子流域出口流量汇入主河道的顺序，在河道节点处进行流量合并，然后再汇入下一河段进行计算。先合后演法更符合实际的流量过程。

在利用先合后演法进行河道马斯京根分段连续演算时，没有采用等距分段的原则，而是根据自然河段的长度，分别设置每一河段的时间差分权重系数和洪水波在河段中的传播时间，以体现不同河段之间的差别。

模型的汇流策略采用模块搭配的方式来解决。首先模型提供基本的汇流演算方法，有基于栅格水系的马斯京根-康吉运动波模拟法、滞时演算法、马斯京根法(先演后合、先合后演)和基于河道的马斯京根分段汇流演算法；然后根据径流类型和模拟的时段长度来选择不同的汇流演算算法，得到适合的汇流方案组合；最后根据所选用的汇流方案组合，模拟得到流域的汇流过程。

模拟过程中率定得到江口流域次洪过程的模型参数如下表7。

表7：江口流域次洪过程模型率定参数表

表中：SurfQ-表示地表径流；LatQ-表示壤中径流；BaseQ-表示地下径流，下同。

根据以上的模型运行方案和率定得到的模型参数，模拟得到江口流域的暴雨径流过程，如图4。

对6次洪峰过程的模拟结果进行精度分析如表8所示。

表8：江口流域次洪过程误差分析表

从模拟结果的对比分析来看，6场洪水过程模拟的相关系数均较高，除了19840709次洪水的相关系数为0.86以外，其余5场洪水过程的相关系数都达到0.90以上，说明模型模拟的洪峰流量过程曲线与实际的流量过程曲线之间有很好的相似性，能够正确描述洪峰过程的变化趋势。洪峰流量的误差在-23.3％～+10.9％之间，峰现时间误差3场洪水为0，3场洪水为2小时，均达到洪峰预报的要求。

模型的确定性系数最低为0.73，最高为0.96，其中4场洪水过程模拟的确定性系数在0.80以上，模拟精度普遍较高。确定性系数的高低与降雨分布有关，当降雨量大、降雨过程集中、洪峰流量较大时，模型的确定性系数较高；而当降雨量比较分散、洪峰流量较小时，确定性系数稍低。以确定性系数为标准，有2场的有效性等级达到甲等，4场达到乙等，表明模型在湿润区和半湿润区流域进行洪水预报时可以达到作业预报的要求。

4.日、月、年尺度水文过程模拟和精度评价

日、月、年尺度的水文过程与暴雨水文过程有不同的水文响应特性，模型需要设计不同的运行方案并率定与所选方案相关的模型运行参数。

江口流域长时段降雨～径流过程具体的模型运行方案表如表9所示：

表9：江口流域长期降雨径流过程模型运行方案表

针对以上运行方案，根据江口流域1981～1985年的水文、气象资料，率定出长时段尺度下降雨～径流过程模拟的模型参数如表10所示。

表10：江口流域长期降雨径流过程模型率定参数表

江口流域日降雨径流模拟结果具体见图5所示。

对日模拟结果按月、年进行统计平均，得到逐月流量图(图6上)、月平均流量图(图6中)和年平均流量图(图6下)。

长时段降雨～径流过程的精度分析和评价如表11所示：

表11：江口流域日、月、年降雨径流过程误差分析表

模拟时间段	确定性系数	相关系数
			日模拟	0.75	0.87
月模拟	0.89	0.95
			年模拟	0.95	0.99

从误差分析表可知，长时段水文过程模拟的精度比次洪过程模拟有一定程度的降低。1981～1985年的日过程模拟确定性系数为0.75，经分析认为，影响模拟精度提高的主要原因有以下两点：

1)流域土壤对水量的调节作用不是很明显，导致模型对降雨过程的响应比较灵敏，一个稍大的降雨过程就有一个对应的流量峰值过程。而实际上在土壤前期含水量较高的流域，由于土壤对入渗水量的调蓄作用，流量过程应该比较平稳。

2)降雨时间的长短将决定雨强的大小，进而影响超渗产流的计算。模型考虑了超渗产流和蓄满产流的同时作用，因此在长时段降雨～径流过程中根据流域的降雨特性来综合确定一个比较合适的日降雨时长。由于降雨时间长度的均化作用和汇流过程的均化作用弱化了超渗产流量在整个栅格产流中的贡献，所以会影响到最终的径流模拟精度。而对于次洪过程来讲，由于时间步长较短，雨强相对接近实际，所以能够得到更高的模拟精度。

从月模拟和年模拟来看，确定性系数分别达到0.89和0.95，说明模型完全能够在月和年尺度上达到径流过程模拟和预测的要求。

实施例二

寒旱区水文过程模拟与精度评价。

1.研究区概况

流域内景观垂直分布明显，森林主要分布于中山地带，灌木和牧草分布在流域各处。以海拔3600m为高山冰雪冻土带和山区植被带的分界线，4500m以上为永久冰川积雪带，以下依次为高山草甸与灌丛、山间盆地、中山森林和中山草甸草原带、低山荒漠带(黄清华，张万昌；SWAT分布式水文模型在黑河干流山区流域的改进及应用[J]；南京林业大学学报(自然科学版)；2004年02期)。流域年降水量约400mm左右，年蒸发能力约1600mm左右，降水量年内分配极为不匀，夏季7～8月占年降水量的70％以上，冬季降水小于5％，属于典型的干旱、半干旱区流域。

2.研究区基础数据

A.水文、气象数据

系统收集了黑河干流山区流域附近10个雨量站的实测降雨量资料和8个蒸发站的实测蒸散发资料，用于次洪过程模拟。最终的流量过程验证站点为流域出口——莺落峡站。

日降雨～径流过程模拟的时间段为1990-2000年，其中1990-1995年的资料用于率定模型参数，1996-2000年的资料用于模型验证。对照11年的流量过程线，结合实际数据情况，摘录了6场典型洪峰过程，洪号分别为19900724次、19960726次、19960820次、19980703次、19980713次和19980915次。

B.DEM数据

黑河山区流域的DEM数据由1∶25万比例尺的黑河流域地形图数字化生成，空间分辨率为120m，采用UTM投影，WGS84坐标系，第47分带。

C.土壤数据

流域的土壤图来源于1∶100万比例尺黑河流域数字化土壤图，为了和DEM相匹配，采用与DEM相同的空间分辨率、投影类型和坐标系，如表12。

表12：黑河流域土壤类型分类表

暗栗钙土	高山寒漠土	淡栗钙土	高山草甸土	栗钙土
					高山草甸草原土	黑钙土	灰钙土	盐化灰漠土

D.植被数据

流域的植被类型图来源于1∶100万比例尺黑河流域数字化土地利用/覆盖图。

3.次洪过程模拟和精度评价

对黑河山区流域摘录的6次典型洪水过程，模型采用的运行方案与汉江上游区江口流域一致(表13)，但是率定的次洪过程模型参数有所不同(表14)。

表13：黑河山区流域次洪过程模型运行方案表

表14：黑河山区流域次洪过程模型率定参数表

根据以上率定的模型参数，模拟得到黑河山区流域的暴雨径流过程。

对6次洪峰过程的模拟结果进行精度分析如表15所示：

表15：黑河山区流域次洪过程误差分析表

以确定性系数为标准，有2场的有效性等级达到甲等，4场达到乙等，表明模型在干旱区、半干旱区的次洪过程模拟有较高的模拟精度，可以达到洪水作业预报的要求。

从两个流域各自的流量过程曲线来看，江口流域对暴雨～径流过程响应的时间比黑河山区流域短，发生暴雨事件后紧接着就会发生产流事件；而黑河山区流域的径流过程对于暴雨过程而言有相对较长的滞后期，体现出首先土壤对降雨有一个较长时间的吸收和调蓄，然后再产流的过程。两个流域的流量过程表现出湿润区超渗产流和干旱区蓄满产流的特征，这是和传统的产流模型理论相违背的。但是在采用了渗蓄一体化动态产流机制后，可以同时模拟栅格超渗产流和蓄满产流之间的转换，也就是说，在同一个栅格上可以同时模拟出超渗产流和蓄满产流，并确定出两者在地表产流中的比例。如果在干旱区能够模拟出较大的蓄满产流，就可以模拟土壤对降雨的吸收和调蓄作用，从而产生径流事件的滞后效应；而如果在湿润区模拟出较大的超渗产流，就可以缩短径流过程相对于暴雨过程的响应时间，出现在暴雨事件后紧接着发生产流事件的现象，这是单纯采用干旱区超渗产流、湿润区蓄满产流所模拟不出来的，也从实际上证明了渗蓄一体化动态产流机制的先进性和正确性。

4.日、月、年尺度水文过程模拟和精度评价

黑河山区流域的日、月、年尺度水文过程模拟的模型运行方案与江口流域一致。根据流域1990～1995年的水文、气象资料，率定出模型运行参数如表16所示。因为黑河山区流域由于属于高寒山区，冬季径流以深层地下基流补给为主，所以在率定参数上比江口流域多增加了一个深层地下基流调节量参数，其作用类似于慢速地下径流。

表16：黑河山区流域日、月、年降雨径流过程模型率定参数表

A.模型率定期水文过程模拟精度分析

采用表16率定的模型参数，模拟得到黑河山区流域率定期日降雨～径流过程。

对日模拟结果按月、年进行统计平均，得到逐月流量图、月平均流量图和年平均流量图。

流域率定期长时段降雨～径流过程精度分析如下表17所示：

表17：黑河山区流域率定期日、月、年降雨径流过程误差分析表

模拟时间段	确定性系数	相关系数
			日模拟	0.68	0.84

[0331]

月模拟	0.92	0.96
			年模拟	0.69	0.85

B.模型验证期水文过程模拟精度分析

为了验证模型在寒旱区的适应性，采用率定期得到的模型运行参数，对流域1996～2000年的数据进行了验证。

同样对日模拟结果按月、年进行统计平均，可得到逐月流量图、月平均流量图和年平均流量图。

流域验证期长时段降雨～径流过程精度分析如表18所示。月、年流量的相对误差分析结果如表19所示。

表18：黑河山区流域验证期日、月、年降雨径流过程误差分析表

模拟时间段	确定性系数	相关系数
			日模拟	0.70	0.85
月模拟	0.90	0.96
			年模拟	0.81	1.00

验证期的日、月、年水文过程模拟的确定性系数和相关系数总体来看和率定期接近，但是稍有提高。月、年流量相对误差的量值和误差分布规律也和率定期基本一致，体现出所率定的参数在黑河山区流域有一定的代表性，同时模型对寒旱区的水文过程模拟有良好的稳定性。如表19。

表19：黑河山区流域验证期月、年流量相对误差分析表

Claims

1.基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，主要步骤为：

A)分布式参数获取；

B)流域水文过程概化；

C)渗蓄一体化产流过程设计；

D)汇流过程设计。

2.根据权利要求1所述的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，其特征在于步骤A)中，获取包括气象、土壤、植被和水系在内的分布式参数：

(1)气象参数缺测条件下的预测

使用WGEN气象发生器，实现降雨量、气温、太阳辐射、相对湿度和平均风速五类气象参数的随机解集计算，用于站点气象资料不全时的填补、气象数据缺测时的模拟以及未来气象参数的预测；

使用邻近站点法、反距离加权法、泰森多边形法、复杂的梯度距离平方反比法或于DEM的PRISM模型法用于模型的参数空间离散，通过GIS空间离散，得到每一个栅格的气象参数，以便进行栅格水文过程模拟；

(2)土壤属性参数库的建立

对于每一种土壤类型，需要外部输入的土壤参数主要分为两类，一类是土壤分层参数，描述具体每一分层的土壤水力物理特性；另外一类是土壤剖面参数，描述该土壤类型下与本栅格所有土层有关的土壤剖面信息；

土壤分层参数：土壤深度、土壤湿容重、有效持水量、饱和导水率、稳定下渗率系数、初始下渗率系数、有机物含量和土壤粒径组成；

土壤剖面参数：土壤类型名称、土壤分层数、植被根区深度、地表反照率、霍顿下渗曲线系数和初始土壤水比例；

(3)植被属性参数库的建立

植被参数分为土地利用/植被覆被参数和植被生态参数两类，植被类型参数的获取，是由国家或地方林业部门提供流域的土地利用/覆被图，通过GIS数字化以及矢量～栅格格式转换的方式，得到栅格植被类型分布图；植被生态参数是利用遥感反演结合实地调查验证的方法获得的；

植被生态参数：最大叶面积指数、最小叶面积指数、莱峰偏移、地表噪度系数、月平均叶面积指数、月平均反射率和月平均覆盖度；

(4)基于DEM的流域水系参数提取

以D8算法即最大落差算法为基础，采用基于图论的广度优先搜索算法实现流域所有栅格的遍历，实现模型所需的水系参数的自动提取。

3.根据权利要求2所述的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，所述基于DEM的流域水系参数提取的步骤包括：水流流向矩阵计算、水流累积矩阵计算、河网自然水系提取和流域边界计算；

a.水流流向矩阵计算

原始的DEM首先通过填洼处理，得到填洼后的DEM，在此基础上进行水系和汇流参数提取；

b.水流累积矩阵计算

在水流方向矩阵的基础上，通过对任意一个栅格利用基于图论的广度优先搜索算法遍历所有与本栅格有入流关系的上游栅格，可以建立水流累积矩阵；

c.河网自然水系提取

通过给定合适的上游给水区面积阈值，将水流累积量矩阵数据中高于此阈值的格网连接起来，在此基础上搜寻出水系的所有栅格点，形成完整的自然水系网络；

d.流域边界计算

4.根据权利要求2所述的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，所述步骤B)中，流域的陆面水文过程涉及到水量在大气-植被-土壤之间的循环和分配，是一个复杂的时变非线性系统，在不同的空间地域单元内，由于下垫面的空间变异性，水文循环过程差异很大，但是对于同一个地域单元，可以是子流域、高程带、水文响应单元或者栅格，如果空间尺度足够精细，可以认为单元内部下垫面特征是均一的，在相同的植被和土壤类型下，有相同的水文响应特征；因此，对于概念性分布式水文模型，模型的产流过程是在地域单元上实现的，在垂直方向，大气水以降雨或降雪的方式降落，首先接触到植被，植被冠层对雨水进行截留以满足自身的需水和蒸发、蒸腾需要；多余的水量降落到地表，产生地表水下渗，下渗到土壤中的水在土壤和植被根系的作用下在土壤中进行再分配处理，在合适的情况下产生向上运动的蒸散发和侧向流动的地表径流、壤中径流和地下径流，完成单元的产流过程；单元产流产生的径流分量在水头的作用下，沿坡面方向产生运动，完成单元间的水量交换，通过河道汇流运动到流域的出口；对陆相水文循环过程进行简化，将耦合的水文子过程线性化处理，分解为相对独立的栅格单元水文子过程，以能量平衡和水量平衡为基础，建立概念性分布式流域水文模型；

其中：冠层截留模拟采用改进的概念性Aston指数模型；融雪径流计算采用常用的能量平衡法或气温指标法，潜在蒸散发计算采用Penman-Monteith(PM)法、Pristley-Taylor(PT)法或Hargreaves法；地表水下渗计算包括霍顿模型(Horton R.E.)、科斯加可夫模型(Kosjiakov，A.N.)、霍尔坦模型(Holton)、菲利普模型(Philip)或格林-安普特模型(Green-Ampt)；根据土壤分层水量调节参数土层孔隙率、田间持水量时植被可吸收水量、土层每日实际可供水量、土壤水饱和时土层可供水量、田间持水量时土层含水量、枯萎点时土层含水量和土壤水力传导率系数，动态计算每一分层对入渗水的调节状态；根据土壤剖面水量调节参数土壤剖面平均湿容重、土壤剖面平均孔隙率、土壤剖面平均饱和导水率、土壤剖面稳定下渗率、土壤实际含水量、土壤饱和时土壤含水量、田间持水量时土壤含水量和枯萎点时土壤含水量，进行土壤剖面水量调节计算。

5.根据权利要求2所述的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，其特征在于步骤C)中，除了常规的栅格水量在垂直方向上完整的输入和输出，还考虑了栅格水平方向的出流，本发明同时顾及栅格垂直方向和水平方向出流，设计了一个完整的渗蓄一体化栅格产流方案，具体如下：

(1)动态产流方案设计

假定D-土壤缺水量(mm)；P-到达地面的净雨(mm)；E-土壤实际蒸散发(mm)；I-雨强(mm/hr)；FC-土壤表层稳定下渗率(mm/hr)；FP-土壤表层时段下渗率(mm/hr)，

I.当D＜0时(前期土壤水已蓄满)

若P-E＜0，I＜F_C，产流方式：不产流； ①

若P-E＜0，I＞F_C，产流方式：超渗产流； ②

若P-E＞0，I＜F_C，产流方式：蓄满产流； ③

若P-E＞0，I＞F_C，产流方式：超渗产流+蓄满产流； ④

II.当D＞0时(前期土壤水未蓄满)

若P-E＜D，I＜F_P，产流方式：不产流； ⑤

若P-E＜D，I＞F_P，产流方式：超渗产流； ⑥

若P-E＞D，I＜F_P，产流方式：蓄满产流； ⑦

若P-E＞D，I＞F_P，产流方式：超渗产流+蓄满产流； ⑧

按前期土壤水是否蓄满进行预分类的目的在于土壤下渗率是个时变的值，如果前期土壤水已经蓄满，则在计算时段内采用土壤稳定下渗率来和雨强作对比，否则采用土壤时段下渗率来和雨强作对比，以决定是否发生超渗产流；

(2)栅格产流计算

取栅格为水文模拟单元，在高空间分辨率条件下，近似认为栅格内的土壤、植被特征是均一的，这样，针对8种不同的产流方案和水量平衡理论，按不同的前期土壤水状态，可以分别计算出各部分的水量分配情况；计算过程中的变量定义如下：R_T-总产流；R_S1-超渗地表径流；R_S2-蓄满地表径流；R_L-壤中径流；R_G-地下径流；PE＝P-E；IF_C＝I-F_C；IF_P＝I-F_P；Δw-土壤含水量的变化量；

I.地表径流计算

●产流类型①——不产流：

R_T＝0，Δw＝PE＜0 (1)

R_S1＝R_S2＝0 (2)

R_L＝R_G＝0 (3)

●产流类型②——超渗产流：

产流类型②的雨强大于土壤稳定下渗率，满足超渗产流条件，所以虽然总的降雨量不够土壤和植被蒸发，但是净雨首先产生超渗地表径流，然后不够蒸发的部分再从土壤水中补充；各部分产流关系如下：

R_T＝R_S，Δw＝PE-R_T (4)

R_S1＝IF_C·Δt (5)

R_S2＝0 (6)

R_L＝R_G＝0 (7)

●产流类型③——蓄满产流：

在雨强小于土壤稳定下渗率的条件下，前期土壤水蓄满后，外部供水除了满足土壤和植被蒸发，还有多余的水分下渗，满足蓄满产流的条件，此时的产流为蓄满产流；由于土壤水始终处于饱和状态，所以土壤含水量变化Δw为0；各部分产流关系如下：

R_T＝PE，Δw＝0 (8)

R_S1＝0 (9)

R_G+R_L＝F_C·Δt (10)

R_S2＝R_T-(R_G+R_L) (11)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf (12)

式中：d_Surf为地表径流R_S2与地表以下径流R_G+R_L的分流比(0≤d_Surf≤1)；设置该变量的目的在于如果Δt取得不合适而导致R_S2小于0时，可以通过d_Surf对R_S2进行调整：如果模拟流域地表径流影响比较大，则d_Surf取较大值；如果模拟流域的径流补给以地下水为主，则d_Surf取较小值；此时，壤中流和地下径流的关系调整如下：

R_G+R_L＝R_T·(1-d_Surf) (13)

●产流类型④——超渗产流+蓄满产流：

前期土壤水蓄满后，雨强大于土壤稳定下渗率，同时外部供水能够满足土壤和植被的蒸发以及土壤水的下渗，此时，将会发生超渗产流和蓄满产流同时存咱的状况；在这种情况下，在控制总水量平衡时，首先满足地表产流，然后计算地下径流和壤中流；各部分产流关系如下：

R_T＝PE，Δw＝0 (14)

R_S1＝IF_C·Δt (15)

R_G+R_L＝F_C·Δt (16)

如果R_S1+R_G+R_L＞R_T，则有：

R_G+R_L＝R_T-R_S1 (17)

R_S2＝0 (18)

否则，有：

R_S2＝R_T-R_S1-F_C·Δt (19)

在式4-74中，如果R_G+R_L＜0，则取

R_S1＝R_T·d_Surf (20)

R_G+R_L＝R_T·(1-d_Surf) (21)

●产流类型⑤——不产流：

在产流类型⑤的条件下，雨强小于土壤时段下渗率，降雨量不够土壤和植被蒸发，不足部分水量从土壤前期含水量中补充，因此，在该条件下，不产生地表径流；由于前期土壤水未蓄满，所以计算不足部分水量的时候，需要计算该水量是否超过土壤前期含水量与枯萎点水量之差；如果超过，则取此时土壤能够提供的自由水量来代替不足水量；各部分产流关系如下：

R_T＝0，Δw＝PE (22)

R_S1＝R_S2＝0 (23)

R_L＝R_G＝0 (24)

如果PE＜SP_Sw-SP_Wp，则：

Δw＝SP_Sw-SP_Wp (25)

●产流类型⑥——超渗产流：

产流类型⑥和产流类型②类似，不同之处在于产流类型⑥的前期土壤水没有蓄满，所以超渗地表径流用时段下渗率来计算，而产流类型②用稳定下渗率来计算；各部分产流关系如下：

R_T＝R_S，Δw＝PE-R_T (26)

R_S1＝IF_P·Δt (27)

R_S2＝0 (28)

R_L＝R_G＝0 (29)

●产流类型⑦——蓄满产流：

D = (1 - \frac{SP_Sw}{SP_Fc}) \cdot SP_Por \cdot Σ_{i = 1}^{n} SL_Z (i) - - - (30)

R_T＝PE-D，Δw＝D (31)

R_S1＝0 (32)

R_G+R_L＝F_P·Δt (33)

R_S2＝R_T-(R_G+R_L) (34)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf (35)

此时，壤中流和地下径流的关系调整如下：

R_G+R_L＝R_T·(1-d_Surf) (36)

●产流类型⑧——超渗产流+蓄满产流：

在产流类型⑧中，外部供水满足了土壤缺水量D，有多余的水量用于产生蓄满地表径流；同时雨强大于土壤地表时段下渗率，满足超渗产流的条件，因此，超渗产流和蓄满产流同时存在；各部分产流关系如下：

R_T＝PE-D，Δw＝D (37)

R_S1＝IF_P·Δt (38)

R_G+R_L＝F_P·Δt (39)

如果R_S1+R_G+R_L＞R_T，则有：

R_G+R_L＝R_T-R_S1 (40)

R_S2＝0 (41)

否则，有：

R_S2＝R_T-R_S1-F_P·Δt (42)

如果R_G+R_L＜0或者R_S1＞R_T，则对径流分量调节如下：

R_S1＝R_T·d_Surf (43)

R_G+R_L＝R_T·(1-d_Surf) (44)

II.壤中径流和地下径流计算

SP_Sw_i＝SP_Sw_i-1+PE_i-R_S1，i-R_S2，i-(R_G，i+R_L，i) (45)

其中，下标i和i-1分别表示当天和前一天的水量值(mm)；通过对时间步长的循环计算，可以得到每一个计算时段的土壤含水量状况；

在地表径流计算过程中，计算得到的是壤中径流和地下径流的总和，通过对土壤水量下渗的处理，来划分两者的比例关系；

土壤水下渗计算分为两种情况，一种是土壤含水量超过田间持水量时的下渗，另外一种是土壤含水量在枯萎点水量和田间持水量之间时的正常下渗；

W_{P} = (SP_Sw - SP_Fc) \cdot [1 - \exp (- \frac{Δt}{Δ t_{P}})] - - - (46)

Δ t_{p} = \frac{(SP_Sat - SP_Fc)}{SP_Sat_K} - - - (47)

此时的壤中径流定义为下渗水量的函数，有

R_L＝W_P·d_L (48)

式中：d_L为外部率定的壤中径流出流因子，一般取0≤d_L≤1；

W′_P＝F_P·Δt (49)

W′_P＝SP_Sw-SP_Wp (50)

W_{P} = W_{P}^{'} \cdot [1 - \exp (- \frac{W_{P}^{'}}{SP_Sw})] - - - (51)

壤中径流量同样由上式来计算得到；这样，当计算出壤中径流后，由壤中径流和地下径流之和减去壤中径流，就可以分离出地下径流量；

6.根据权利要求4所述的基于渗蓄一体化动态产流机制的分布式水文模型的设计方法，在栅格上计算出各径流分量后，通过坡地和河道汇流处理，可以将水量汇集到流域的出口，从而模拟得到流域的流量过程；在步骤D)中，针对次洪过程和长时段降雨～径流过程的不同特征以及坡地汇流和河道汇流的差别，设计了马斯京根-康吉法、滞时演算法和河道分段马斯京根法，用于不同条件下模型的汇流过程演算；汇流策略采用模块搭配的方式来解决；首先模型提供基本的汇流演算方法，有基于栅格水系的马斯京根-康吉运动波模拟法、滞时演算法、马斯京根法或基于河道的马斯京根分段汇流演算法；然后根据径流类型和模拟的时段长度来选择不同的汇流演算算法，得到适合的汇流方案组合；最后根据所选用的汇流方案组合，模拟得到流域的汇流过程。