CN102034001A

CN102034001A - 一种以栅格为模拟单元的分布式水文模型设计方法

Info

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

Abstract

本发明公开了一种以栅格为模拟单元的分布式水文模型设计方法，简称ESSI。其步骤为：分布式参数获取，矢量数据转换为栅格数据；流域水文过程概化，建立栅格通用产流类型；产流过程设计，按不同的前期土壤水状态，可以分别计算出各部分的水量分配情况；汇流过程设计，分别设计马斯京根－康吉法、滞时演算法和河道分段马斯京根法，用于不同条件下模型的汇流过程演算；模型开发集成。本发明提以模块化、集成化思想为手段，不仅实现可以完成干旱区和湿润区的流域水文过程模拟，而且可以实现流域的短期洪水预报和长期降雨～径流过程模拟和预测，从而为人类深入了解水循环的物理机理，减少旱涝灾害，合理开发和利用水资源提供科学的参考依据。

Description

一种以栅格为模拟单元的分布式水文模型设计方法

技术领域

本发明以涉及一种水文模型的设计方法，更具体的说是一种以栅格为模拟单元的分布式水文模型(简称ESSI)的设计方法。

背景技术

分布式水文模型作为一种探索水文循环机制、定量评估水资源动态变化的有力工具，不仅要在理论上向人们合理解释水文循环的物理基础，说明水是如何在SVAT***中运动并驱动物质和能量的循环；而且更要有实用性和可操作性，因为只有具有良好实用性和可操作性的模型***，才能作为一个有效的水文过程模拟平台，在对实际水问题的分析和模拟过程中为国民经济的发展提供正确、合理的分析结果和建议。所以尽管当前分布式水文模型发展很快，也取得了一定的社会经济效益，但分布式水文模型还存在进一步发展的空间。从以下五个方面来看：

首先，基于水力学计算的分布式水文模型虽然物理机制突出，但是输入数据的准备工作繁琐，专业性强，而且对模型的模拟单元尺度有一定的限制，原因在于当模型运用到一个较大的流域时，计算开销量与模拟单元尺度的大小呈指数关系增长，因此模型不太适用于较大流域的水文过程模拟，而实际上目前中、大尺度流域甚至是全球尺度的水文过程模拟分析更有实用价值。这一类模型目前需要解决的主要问题在于协调网格单元的尺度与计算量以及模拟精度之间的平衡关系：网格精度越高，在一定程度上模型模拟越精确，但同时对分布式的输入数据要求以及计算机开销要求也越高；反之，网格精度降低，网格与网格之间的水力学关联就不那么紧密，相应地削弱了模型的水力学物理基础。因此除了模拟流域具备完整的模拟资料，而且计算机的计算和存储技术有很大改进，那么具有水力学物理基础的分布式水文模型由于尺度效应问题将很难在短期内运用到较大尺度的流域上。从当前的研究成果来看，具有物理基础的分布式水文模型和概念性分布式水文模型在模拟精度上并没有明显的区别，因此，将RS、GIS、DEM技术和概念性水文模型结合，在自然流域或相对较大尺度的栅格上进行分布式水文过程模拟，对水资源进行动态分析和评估，具有良好的应用前景和实际应用价值，是当前分布式水文模型发展的一个重要方向。

其次，在分布式模拟技术和概念性水文模型结合的研究领域，国内外在RS、GIS、DEM等先进技术和概念性水文模型的耦合应用方面所采取的技术路线比较类似，但是，耦合应用在以下方面还可以进一步发展和提高：1)有效解决遥感数据的瞬时性和水文模拟的时序性之间的关系；2)充分利用遥感影像数据源，尽可能多地反演水文模型参数，如初始的土壤湿度场、地表反照率、地表粗糙度、NDVI、LAI等，并有效提高参数反演的精度；3)解决遥感数据尺度和水文模型的模拟单元尺度之间的匹配，提出水文模拟单元分辨率取值的合理建议方案；4)利用GIS技术，结合数学手段，在模型参数的空间离散化方面发挥更大的作用，使模型参数更加符合流域实际的气象、水文空间分布特征；5)GIS技术强大的可视化表达和结果分析功能可以使水文模型更好地分析和显示水文模拟结果。

第三，先进技术的应用可以为分布式水文模型带来方法和手段上的更新，但并不能代表模型本身的高效性，需要通过分布式模拟技术加上先进的产汇流机理，才能构建出一个好的实用的分布式水文模型，因此在模型核心的产汇流机理研究方面还有很大的发展空间。

产流过程应该是超渗和蓄满产流机制之间的有机结合和动态转换，只不过何时采用何种产流机制，两者之间的比例如何分配，起决定性作用的应该是流域的降雨量特征和下垫面土壤含水量特征的对比，而不是根据气候条件的干旱和湿润与否，在模型建模前决定采用那种产流机制。那么，采用何种方式来平衡两种已知的而且经过充分验证的产流机制，解决模型的产流问题，使模型能够适用于各种气候、植被、土壤状况下的水文过程模拟，就对我们提出以下的问题：

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

从空间上来讲，湿润地区可通过蓄满产流来模拟，干旱地区可通过超渗产流来模拟，那么，半干旱、半湿润地区用哪种产流机制来模拟？随着气候条件和下垫面条件的变化，处于蓄满和超渗产流过程中间的产流如何计算？两者的比例如何划分？从时间上来讲，流域的产流是一个动态变化的过程，那么如何来表达产流的这种动态变化？

如何考虑地形对产流过程的影响，使产流更加符合实际情况？这些问题都对我们研制分布式水文模型的通用性产流机制研究提出了很高的要求。

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

第四，在模型的开发和集成，作为一个优秀的模型***，应该有好的用户界面、简单方便的操作性能和良好的***强壮性。现有模型显得规模比较小，功能不够完善，而且***的集成度低，难以满足不同用户的实际需求。

第五，在模型的实用性。现有分布式模型由于起步比较晚，虽然继承了先前模型的开发经验，但是由于开发周期短，普遍存在模型测试不够的问题，所以虽然出现了很多的分布式水文模型，但是到目前为止，基本上还没有比较成熟的或者得到国际上普遍认可的分布式水文模型。

发明内容

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

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

一种以栅格为模拟单元的分布式水文模型设计方法，主要包括以下步骤：

A)分布式参数获取；

B)流域水文过程概化；

C)产流过程设计；

D)汇流过程设计；

E)模型开发集成。

步骤A)中，分布式参数获取包括气象、土壤、植被和水系参数获取。在传统观测技术和先进的空间技术支持下，我们将模型驱动参数分为两种类型，一种是单点数据，另外一种是空间数据，具体如表1所示：

表1：模型驱动参数表

1)气象参数缺测条件下的随机解集计算和精度评估

分布式水文模型中，气象参数是影响水文循环的重要驱动因素，通过它们可以分析陆面水文过程中水量的输入(降雨、降雪)和输出(蒸发、蒸腾)等水分循环过程。ESSI模型中主要的气象驱动参数有：降雨量、气温(最高、最低和平均)、太阳辐射、相对湿度(实际水汽压)、平均风速和潜在蒸散发，它们的获取有三种途径：1)直接采用流域附近水文和气象站点的长期气象观测参数；2)利用气象发生器(Weather GENerator，WGEN)模拟所需的气象参数，主要用于实现站点气象资料不全时的填补、气象数据缺测时的模拟以及对未来气象参数的预测(陈明昌等，1994；S.L Neitsch，et al.，2002；张东等，2004)；3)通过气候学计算的方法，结合气候的季节性变化规律得到时序的气象参数，比如通过日照时数计算太阳辐射，通过气温、太阳辐射、相对湿度、平均风速计算潜在蒸散发等。

本发明以数理统计、随机模拟和弱稳定理论为基础，使用气象发生器(weathergenerator，WGEN)，实现降雨量、气温、太阳辐射、相对湿度和平均风速等五类气象参数的随机解集计算，用于站点气象资料不全时的填补、气象数据缺测时的模拟以及未来气象参数的预测。通过潜在蒸散发计算结果的间接验证，表明随机解集可应用于水文过程模拟。

考虑空间离散方法的可操作性、计算机实现的难易程度、空间离散计算的效率以及离散精度，使用邻近站点法、反距离加权法、泰森多边形法和复杂的梯度距离平方反比法和基于DEM的PRISM模型法，这些不同方法开发成独立模块，供模型组合时调用，用于模型的参数空间离散。

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

土壤参数是分布式水文模型中非常重要的下垫面参数，它控制着土壤对流域水量的调蓄。当降雨发生时，土壤作为一个巨大的蓄水体，充分吸收和存储降雨入渗水量；降雨发生后，土壤调蓄的水量在植被根系吸收、地下水重力势的作用下，缓慢地产生向上、向下和侧向运动，驱动水量在大气-植被-土壤***中的循环。因此，土壤参数的精度好坏将在极大程度上影响模型的水文模拟精度。

在ESSI模型中，模拟的水文单元为栅格，因此对土壤参数的描述分为两个方向：水平方向按土壤类型划分，不同栅格可以具有不同的土壤类型，其主要数据来源为国家和地方土地管理部门调查和收集的土壤类型图，通过GIS数字化以及矢量～栅格格式转换，得到栅格土壤分布图；竖直方向根据土壤特性如土壤孔隙率、湿容重、饱和导水率等物理属性的不同进行分层，每一土层分别具有不同的水力响应特性。***设计了模型水文过程模拟时所需要的土壤属性参数，总结了土壤分层参数和土壤剖面参数的计算方法，开发出基于累积频率曲线统计原理的土壤粒径转换程序包，解决了不同来源土壤资料之间土壤粒径标准不一致的难题。

3)植被生态参数的遥感反演和属性参数库的建立

植被生态参数涉及到植被的生长过程植被以及对水分的分配，主要包括植被指数NDVI、叶面积指数LAI、植被反照度Albedo和植被盖度。在ESSI模型中植被参数分为植被类型参数和植被生态参数两类，植被类型参数作为植被生态参数的索引，其获取有两种手段：一种是由国家或地方林业部门提供流域的土地利用/覆被图，通过GIS数字化以及矢量～栅格格式转换的方式，得到栅格植被类型分布图；另外一种是利用遥感反演结合实地调查验证的方法获得。遥感反演方法的优点在于只要有合适的遥感影像，不仅可以得到高分辨率的植被类型图，而且植被类型可以按模拟时间的要求动态更新。

LAI的遥感反演方法主要分为两类：统计模型法和光学模型法(方秀琴等，2003)。光学模型法建立在辐射传输模型的基础上，是一种基于植被的BRDF(双向反射率分布函数)特征的LAI遥感反演模型，该方法具有很强的物理基础，但是模型需要的参数众多，参数难以获取而且精度难以达到模型要求，因此在实际应用中比较少。统计模型法进行遥感定量反演LAI的理论依据是植被冠层的遥感光谱特征(方秀琴，2004)。因为遥感红光波段反射率包含了植冠顶层叶片的大量信息，近红外波段反射率包含了冠层内叶片的很多信息，所以对这些信息进行组合利用，就能够提取出我们所关心的LAI。

统计模型法以LAI为因变量，以遥感光谱数据或其变换形式(如植被指数NDVI)为自变量，建立估算模型，再通过与地表实测LAI建立统计相关，来计算研究区的LAI。该方法形式简洁，计算简单易行，对输入参数要求也不高，是目前LAI定量遥感估算的主要方法。

ESSI模型中LAI遥感反演的主要数据源为Landsat TM/ETM+数据和MODIS数据，两者在时间分辨率和空间分辨率之间存在互补关系，通过尺度转换，可以得到不同尺度的LAI，应用于中小流域的水文过程模拟。

利用统计模型法计算LAI时，首先将遥感影像进行几何校正和大气校正(Gliabert M.A，et al.，1994；Zhang W C，et al.，2000)，在此基础上，选择SR(比值植被指数)、NDVI(归一化差值植被指数)、RSR(缩小的比值植被指数)、SAVI(土壤修正植被指数)和ARVI(大气阻抗植被指数)等五种植被指数来进行LAI的遥感估算，计算公式如下：

表2：植被指数及计算公式(据方秀琴，2004)

植被指数	计算公式
		SR	R_n/R_r
NDVI	(R_n-R_r)/(R_n+R_r)
		ARVI	(R_n-R_rb)/(R_n+R_rb)，其中R_rb＝R_r-γ(R_b-R_r)
RSR	SR*(1-(R_swir-R_swir，min))/(R_swir，max-R_swir，min)
		SAVI	(R_n-R_r)*(1+l)/(R_n+R_r+l)

其中，R、R_r和R_swir分别为Landsat7ETM+的4、3、5通道的反射率；γ为与大气气溶胶类型有关的参数，在植被覆盖区γ＝1.0，裸土区γ＝0.5(Kaufman & Weigand，1992)；R_swir，min和R_swir，max为5通道反射率的最小值和最大值，分别取5波段反射率直方图两端1％处对应的反射率值(L.Brown，et al.，2000)；l为调整参数，在中等植被覆盖条件下，取l＝0.5(Huete，1988)。

根据计算得到的植被指数，与实测的LAI进行统计相关分析，针对不同的植被类型分别建立LAI遥感估算模型，用于研究区的LAI遥感反演(方秀琴，2004)。需要说明的是，对于不同的研究区，遥感估算模型的参数需要重新率定。

草地：LAI＝0.6557*exp(2.4412*ARVI)

针叶林：LAI＝0.7211*SR+0.03

阔叶林：LAI＝0.9154*RSR+0.0815

农作物：LAI＝1.1674*exp(ARVI*1.906)

灌木：LAI＝0.3224*RSR+1.4904

植被地表反照率Albedo表征地球表面不同植被类型对太阳辐射的反射能力，决定着能量在地-气之间的分配比率，主要用于计算地表能量收支状况和水量平衡(徐兴奎等2002)。由于地表反照率受地表覆盖类型和太阳高度角的影响，具有时空分异性，因此利用遥感资料求取区域地表反照率的方法日益受到重视(陈云浩等，2001)。从目前的研究来看，地表反照率的遥感反演方法主要有两类，分别为BRDF模型法和统计分析模型法。BRDF模型法认为自然的下垫面大多为非朗伯体，遥感观测的波段反射率是在特定的太阳高度角和卫星视角下的值，与方向性有关，可以通过获取双向反射分布函数来求出半球反照率。这种方法的理论性较强，比较符合复杂地表的实际情况，但是需要大量的连续时间序列的遥感观测数据，计算速度较慢，目前还不能满足数值计算模型的需要。统计分析模型法则根据各波段在太阳辐射中所占的权重来反演全波段反照率。因为遥感数据是不连续波段光谱辐射能量的反映，所以通过对不同波段反射率赋予不同的权重进行组合，从不连续的单波段反射率可以推算出连续光谱范围的反照率(刘三超等，2003)。波段组合的方法有两种，一是先对地表覆盖进行分类，然后利用TM影像的2、4、7三波段组合得到植被地表反照率(Brest & Gowand，1987)；另外一种是直接采用TM影像的可见光到中红外所有6个波段的组合，获取全波段的植被地表反照率(Wang J，etal..2000)。

植被盖度定义为观测区域内植被垂直投影面积占地表面积的百分比，它是影响地表植被蒸腾和土壤水分蒸发损失过程的重要控制因子(张云霞等，2003)。植被盖度C_p如果通过NDVI来计算，公式如下(Gutman，1998；张佳华等，2000)：

C_{p} = \frac{{NDVI - NDVI}_{\min}}{{NDVI}_{\max} - {NDVI}_{\min}}

式中：NDVI、NDVI_max和NDVI_min分别为逐日、年最大、年最小植被指数。

如果通过叶面积指数LAI来计算，公式如下(P.Droogers，2000)：

C_{p} = 1 - e^{- α \cdot {LAI}_{day}}

式中：α为地表太阳辐射的消光系数，缺省值为0.45。张佳华等(2000)对冬小麦进行研究，取返青前α为0.28，返青后为0.35。这样，可以根据已经计算出的栅格NDVI或者LAI来推算出计算时段的植被盖度，得到栅格植被盖度分布图。

遥感解译得到的植被参数既可以作为日变化曲线上用于插值的控制点数据，也可以作为验证数据来检验这种逐日插值得到的水文生态参数的可靠性，使水文模型的输入参数更接近实际情况。

结合四维数据同化技术，得到水文模型所需的叶面积指数(LAI)、归一化植被指数(NDVI)、地表反照率、植被盖度的时序变化规律。

在此基础上，建立了模型的植被属性参数库，用于长时段的水文过程模拟。

目前在ESSI模型中一共建立了13类植被类型，如下表所示：

表3：植被参数类型表

代码	英文含义	中文含义
			EGNEEDFR	Evergreen Needleleaf Forest	常绿针叶林
EGBROADF	Evergreen Broadleaf Forest	常绿阔叶林
			DENEEDFR	Deciduous Needleleaf Forest	落叶针叶林
DEBROADF	Deciduous Broadleaf Forest	落叶阔叶林
			MIXEDCOV	Mixed Cover	针阔混交林
WOODLAND	Woodland	林地
			WGRASSLD	Wooded Grassland	林质草地
CLSSHRUB	Closed Shrubland	封闭灌丛
			OPNSHRUB	Open Shrubland	开放灌丛
GRASSLND	Grassland	草地
			CROPLAND	Cropland	农田
BARELAND	Bare Ground	裸露地
			URBANLND	Urban and Build-Up	城镇和居民地

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

●水流流向矩阵计算

原始的DEM首先通过填洼处理，得到填洼后的DEM，在此基础上进行其它的水系和汇流参数提取(Martz & Garbrecht，1992)。对每一栅格格网，水流方向定义为水流离开此格网时的指向。D8算法将格网的8个邻域格网进行编码，水流方向定义为：1-东；2-东南；4-南；8-西南；16-西；32-西北；64-北；128-东北。栅格水流流向为其中的任一值。确定水流方向的方法和步骤如下：

a)假定计算区域是另一更大数据区域的一部分，对所有数据边缘的格网赋以指向边缘的方向值；

b)对第一步中未赋方向值的格网，计算其对8个邻域格网的距离权落差；

c)确定具有最大落差值的格网，执行如下步骤：

如果最大落差值小于0，则赋以负值表示此栅格流向未定；

如果最大落差值大于或等于0，且最大值只有一个，则将对应此最大值的方向值作为中心格网的水流方向；

如果最大落差值大于0，且有一个以上的最大值，则在逻辑上以查表的方式确定水流方向。

●水流累积矩阵计算

在水流方向矩阵的基础上，通过对任意一个栅格利用基于图论的广度优先搜索算法遍历所有与本栅格有入流关系的上游栅格，可以建立水流累积矩阵(张乃孝等，1998)。具体方法说明如下：

a)对水流流向矩阵以先行后列的方法处理每一个栅格；

b)处理过程用到两个数据矩阵，第一个是水流流向矩阵，由上面的计算得到；第二个是水流累积矩阵的初始状态矩阵，用作数据处理标识以及存储累积量计算结果；

c)开辟一个队列存放与当前处理栅格有入流关系的栅格坐标，同时初始化水流累积矩阵，初始值为0；

d)对任意一个栅格，分析和它相邻的8个栅格的入流关系，记录有入流关系的而且未处理过的栅格坐标，压入队列，利用广度优先搜索算法遍历队列中的所有栅格坐标，直到遍历完该栅格的所有上游流域栅格；

e)对队列中的栅格反向计算累积流量，计算包括两部分，第一部分对水流流向矩阵进行累积，得到栅格本身入流的栅格数；第二部分对水流累积矩阵进行累积，得到栅格上游对本栅格的贡献。

f)行、列所有栅格处理完毕，则可生成所需的水流累积矩阵。

●河网自然水系提取

DEM中某一个栅格点如果能够形成水系，必须存在一定规模的上游给水区，因此可以通过给定合适的上游给水区面积阈值，将水流累积量矩阵数据中高于此阈值的格网连接起来，在此基础上搜寻出水系的所有栅格点，形成完整的自然水系网络(李志林等， 2001)。

对整个水系的追索和分级按Strahler的水系分级方案进行，具体算法和实现步骤如下：

a)选定流域出口栅格，按广度优先搜索算法将上游流域所有的栅格按汇流等级压入队列，直到流域水系起始栅格点；

b)从起始栅格点开始，往下游搜索，将所有等级为1的河道遍历完毕，再遍历等级为2的河道，一直到等级最高的河道遍历完毕，完成所有栅格河道的Strahler等级赋值工作；

c)从第2级开始有河道交汇的情况，按以下原则确定交汇栅格的等级：如果同级的河流交汇，所形成的河流级数增加一级；如果不同级的河流交汇，所形成的河流级数取两者中较高者。

●流域边界计算

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

本发明以D8算法为基础，采用基于图论的广度优先搜索算法实现流域所有栅格的遍历，设计流域水系参数提取算法，实现ESSI模型所需的水系参数的自动提取，同时算法和数据结构可以非常方便地改进和进一步扩充。

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

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

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

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

栅格产流类型的变化具有以下基本规律：

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

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

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

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

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

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

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

步骤C)中，在知道栅格的总水量的条件下，根据水量平衡理论来决定各分水源水量的多少。一般来讲，栅格产流由地表径流、壤中流和地下径流组成，地表径流通过降雨和下垫面特征的对比已经可以计算，这样，总水量扣除地表径流和实际蒸散发，剩余的水量就是壤中流和地下径流之和，壤中流通过土壤下渗水量来计算，扣除壤中流水量，剩下的水量就是地下径流量。在这中间，作为土壤对水量调节的重要指标，土壤含水量的变化对栅格产流的分配起到非常重要的作用。

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

1)前期土壤水已蓄满条件下的地表径流计算：

产流类型①——不产流

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

R_T＝0，Δw＝PE＜0

R_S1＝R_S2＝0

R_L＝R_G＝0

产流类型②——超渗产流

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

R_T＝R_S，Δw＝P_E-R_T

R_S1＝IF_C·Δt

R_S2＝0

R_L＝R_G＝0

产流类型③——蓄满产流

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

R_T＝PE，Δw＝0

R_S1＝0

R_G+R_L＝F_C·Δt

R_S2＝R_T-(R_G+R_L)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf

式中：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)

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

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

R_T＝PE，Δw＝0

R_S1＝IF_C·Δt

R_G+R_L＝F_C·Δt

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

R_G+R_L＝R_T-R_S1

R_S2＝0

否则，有：

R_S2＝R_T-R_S1-F_C·Δt

如果R_G+R_L＜0，则取

R_S1＝R_T·d_Surf

R_G+R_L＝R_T·(1-d_Surf)

2)前期土壤水未蓄满条件下的地表径流计算

产流类型⑤——不产流

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

R_T＝0，Δw＝PE

R_S1＝R_S2＝0

R_L＝R_G＝0

如果PE＜SP_Sw-SP_Wp，则：

Δw＝SP_Sw-SP_Wp

产流类型⑥——超渗产流

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

R_T＝R_S，Δw＝PE-R_T

R_S1＝IF_P·Δt

R_S2＝0

R_L＝R_G＝0

产流类型⑦——蓄满产流

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

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

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

R_T＝PE-D，Δw＝D

R_S1＝0

R_G+R_L＝F_P·Δt

R_S2＝R_T-(R_G+R_L)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf

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

R_G+R_L＝R_T·(1-d_Surf)

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

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

R_T＝PE-D，Δw＝D

R_S1＝IF_P·Δt

R_G+R_L＝F_P·Δt

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

R_G+R_L＝R_T-R_S1

R_S2＝0

否则，有：

R_S2＝R_T-R_S1-F_P·Δt

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

R_S1＝R_T·d_Surf

R_G+R_L＝R_T·(1-d_Surf)

壤中径流和地下径流计算如下：

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

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

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

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

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

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

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

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

{Δt}_{p} = \frac{(SP_Sat - SP_Fc)}{SP_Sat_K}

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

R_L＝W_P·d_L

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

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

W′_P＝F_P·Δt

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

W′_P＝SP_Sw-SP_Wp

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

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

当计算出壤中径流后，由壤中径流和地下径流之和减去壤中径流，就可以分离出地下径流量。

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

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

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

步骤E)中，选用Visual C++作为开发工具，实现相关动态库和界面的编程，矢量和栅格数据的支持由Esri提供的二次开发组件MapObjects 2.0实现，图表功能采用Borland公司开发的图表控件TChart实现。对模型的体系架构和开发集成技术，总结如下：

1)模型的开发框架为：

以C++为模型主要算法的开发工具，实现栅格水文模拟总控部分动态库的开发和衔接；

通过MapObjects 2.0组件提供模型对矢量数据的显示、分析和处理；

利用ODBC方式实现属性参数和驱动参数的数据库访问；

由Visual C++进行***总集成，开发出***的人机交互界面，实现研究区数据准备、水文过程模拟、模拟结果精度分析和评价等一系列功能的紧密耦合，成为一个完整的分布式水文模型***；

2)模型分为三个部分：

输入部分提供ArcView Shp格式的矢量数据支持和Esri ASCII文本格式的栅格数据支持，采用多线程同步访问技术对栅格进行分块处理，实现大文本栅格数据的快速访问；

栅格水文模拟总控部分由一系列动态库组成，包括气象参数模拟及空间离散动态库、栅格河道信息提取动态库、栅格水文子过程模拟动态库、栅格水量平衡计算动态库和栅格河道汇流处理动态库。各动态库功能上互相衔接，数据上自动流转，不仅可以实现分布式的栅格水文、生态、能量过程模拟，而且可以得到栅格或(子)流域出口的流量过程；

输出部分不仅可以实现模拟结果的精度自动分析和评价，而且可以通过图表的形式直观显示不同水文过程参量的二维空间分布和流量过程的一维时序变化；

3)围绕模型的栅格水文模拟总控需求，定义了六个主体类：

主体类分别为：气象参数模拟类，栅格参数处理类，下垫面特征参数类，栅格产流计算类，水文模拟总控制类和模型运行参数类；

主体类之间的逻辑处理顺序如下：首先，下垫面特征参数类提供流域土壤物理参数和植被生态参数信息，气象参数模拟类提供实测或模拟得到的站点时序气象参数信息，栅格参数处理类实现气象参数的空间离散以及基于DEM的河道信息提取和汇流参数结构填充；然后，在这三个类的支持下，通过模型运行参数类设定合适的模型运行方案和运行参数，水文模拟总控制类开始控制模型的水文过程运行，利用栅格产流计算类计算栅格的水量分配，调用汇流方案进行栅格间的汇流处理，最终得到流域出口栅格的流量过程；

4)模型采用面向对象方式进行开发，对所有的水文子过程进行封装和模块化处理，通过界面定制的方法，人为选择模拟区域最佳的水文处理模块组合，使模型在使用上具备很强的灵活性和可扩充性；

5)模型的开发集成模式采用了以水文模型为主，将GIS功能嵌入水文模型中紧密集成开发的模式。这种方式的不足之处是底层程序开发工作量大，GIS数据分析管理和可视化功能仅限于二次开发组件提供的基本功能，因此不及商业GIS软件强大。但是它的最大优点是模型的可扩展性很强，可以根据实际需要随时扩充新的数据结构，使用灵活，模型执行效率高，这是采用其它的开发模式所不能及的。

本发明与现有技术相比，其改进主要体现如下：

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

(2)以新提出的渗蓄一体化动态产流机制为核心，以模块化、集成化思想为手段，开发和构建了分布式陆面水文过程模型***ESSI。该***是一个具有自主知识版权的分布式水文模型，具有可扩充性好、可操作性强、界面友好、容错性高的特点，经过更多不同类型流域的测试，可以成功应用于我国的水文、水资源模拟和研究中。

有益效果：

本发明提供了以栅格为模拟单元的分布式水文模型(简称ESSI)的设计方法，可以建立一个先进的概念性分布式流域水文模型、：

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

(2)以栅格为水文模拟单元，在方便与RS、GIS和DEM数据匹配的同时，具备与气候模式嵌套的快速扩展能力；

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

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

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

(6)模型开发基于模块化、集成化思想，使水文子过程具有良好的可扩展性和维护性，通过模块搭配，优选流域合理的产汇流组合；通过集成化，使模型的外部数据输入和不同阶段模拟结果的表达更方便，更高效；

(7)模型产品化，具有自主版权，具备良好的用户界面和可操作性特征，通过不同流域的水文过程测试，验证模型的合理性和实用性。

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

四、说明书附图

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

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

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

图4为江口流域水系和汇流参数提取结果；

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

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

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

具体实施方式

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

实施例

为检验模型效率与模拟结果的精度，模型采用确定性系数(Nash & Sutcliffe，1970)和相关系数r来进行分析和评价。

确定性系数NSC定义如下：

NSC = 1 - \frac{Σ_{i = 1}^{n} {(Q_{obs, i} - Q_{sim, i})}^{2}}{Σ_{i = 1}^{n} {(Q_{obs, i} - {\overset{&OverBar;}{Q}}_{obs})}^{2}}

相关系数r定义如下：

r = \frac{Σ_{i = 1}^{n} (Q_{obs, i} - {\overset{&OverBar;}{Q}}_{obs}) \cdot (Q_{sim, i} - {\overset{&OverBar;}{Q}}_{sim})}{\sqrt{Σ_{i = 1}^{n} {(Q_{obs, i} - {\overset{&OverBar;}{Q}}_{obs})}^{2} \cdot Σ_{i = 1}^{n} {(Q_{sim, i} - {\overset{&OverBar;}{Q}}_{sim})}^{2}}}

式中：Q_obs，i为实测流量(m3·s-1)，Q_sim，i为模拟流量(m3·s-1)，

为实测流量均值(m3·s-1)，

为实测流量均值(m3·s-1)，n为模拟的时段数。确定性系数越接近于1，表明模型效率越高；相关系数越大，表明模拟流量过程和实测流量过程形态越接近。

同时，针对次洪过程，通过检验实测洪峰流量和模拟洪峰流量之间的流量误差和峰现时间误差，来描述ESSI模型对洪峰过程模拟的精确程度；针对长时段降雨～径流过程模拟，特别是月、年以上时间尺度的模拟，在模拟时间段比较短的情况下，通过分析流量过程的平均相对误差，辅助分析模型对水量模拟和预测的精度。

另外，采用《水文情报规范》(SD 138-85)进行水文过程模拟结果的有效性评定和检验，其基本标准为：

根据确定性系数来评定等级有效性等级：

1)甲等：确定性系数≥0.90；

2)乙等：0.90＞确定性系数≥0.70；

3)丙等：0.70＞确定性系数≥0.50。

预报方案达到甲、乙两个等级时，可用于作业预报；达到丙等时可用于参考性预报。丙等以下则不能用于作业预报，只能用作参考性估报。

汉江是长江上游第一大支流，发源于秦岭南麓，流经陕西、湖北两省，全长1577km，流域面积15.9×104km2(杨永德等，1997)。汉江以白河和丹江口水库为界，分为三段，白河以上为汉江上游区，丹江口以下为汉江下游区，在武汉市汇入长江。

汉江上游区河长735km，流域面积59115km2，占全流域面积的37.2％。在地理位置上，北依秦岭，南接大巴山、米仓山，西端与嘉陵江相临，北、西、南三面环山，东面开阔平坦。汉中平原附近河道宽浅，水流平缓，洪水时水面宽可达1～2km，有良好的滞洪作用；当河道出洋县以后，河深变窄，河槽调蓄能力很小。

流域的气候属于北亚热带季风气候，具有四季温暖，雨量充沛，干湿分明的特点，多年平均气温12～16℃，极端最高气温42℃，极端最低气温-13℃，多年平均相对湿度74％，最大风速21m/s，多年平均蒸发量848mm。上游支流呈羽叶状分布于干流两侧，左侧支流源高流长，右侧支流源低流短。这些支流坡降大，汇流速度快，发生暴雨时，暴雨产生的洪水过程能够迅速集中。

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

江口流域的气候属于亚热带季风气候，具有明显的垂直地带性，全年降雨多在750～1000mm之间，降雨量年内分配不均，主要集中在6～9月。流域内植被覆盖度高，自然植被保护良好，以温带落叶阔叶林为主，如栓皮栎林、锐齿栎林、红桦林以及一些杨、柳、槭林等杂木林，还有一部分温带落叶灌丛，如胡枝子灌丛、白刺花灌丛、绣线菊灌丛等，另外还包括少部分的温带和亚热带针叶植物如油松林、华山松林以及一些温带草丛和草甸(朱求安，2005)。流域的土壤类型以黄棕壤、棕壤和淋溶褐土为主，土壤分布的空间差异较大。

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

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

表中：SurfQ——表示地表径流；LatQ——表示壤中径流；BaseQ——表示地下径流对6次洪峰过程的模拟结果进行精度分析如下表所示：

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

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

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

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

长时段水文过程模拟与次洪过程模拟的运行方案的不同点有三个：1)长时段水文模拟涉及到冬春季节的枯季径流过程，冬春季节气温寒冷，存在降雪融雪过程，需要考虑融雪径流的影响；2)在冬季进行蒸散发模拟时，采用互补相关法和Kojima法联合模拟的方案进行计算，两种方法之间通过气温阈值来确定方法的选取；3)次洪过程的坡面汇流方案采用马斯京根-康吉法来模拟洪水波在坡面的演进，而对于长时段水文过程模拟来说，时间步长较长，不容易满足坡面洪水波汇流的条件，因此选择滞时演算法来进行坡面汇流处理。

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

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

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

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

与次洪过程相比，可以看出长时段模拟下的滞时演算汇流系数比次洪过程减小很多，说明长时段下的水文过程模拟需要更长的汇流时间。这些参数可根据流域最远栅格的产流到达流域出口的大致时间来估算得到，因为对于一个流域来说，在重力势和沿程损失的共同作用下，汇流时间可看作是流域固有的特征量。在汇流时间相对固定的情况下，汇流沿程传输损失系数决定了栅格汇集到流域出口的水量贡献，而在湿润区地表径流是主要的径流分量，所以模型特别设计了枯、平、丰水年的地表径流传输损失系数，来表达在长时段水文过程模拟中降雨量的年际变化规律对流量过程的影响，并取得好的模拟效果。

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

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

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

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

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

1)流域土壤对水量的调节作用不是很明显，导致模型对降雨过程的响应比较灵敏，一个稍大的降雨过程就有一个对应的流量峰值过程。而实际上在土壤前期含水量较高的流域，由于土壤对入渗水量的调蓄作用，流量过程应该比较平稳。因为影响土壤对水量的调蓄作用的关键在于土壤属性参数是否正确，所以如果能够进一步提高土壤属性参数的精度，就有可能提高水文过程模拟的精度；

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

从月模拟和年模拟来看，确定性系数分别达到0.89和0.95，说明ESSI模型完全能够在月和年尺度上达到径流过程模拟和预测的要求。但是对于月、年模拟来说，由于模拟时段比较短，数据量少，确定性系数不能全面说明问题，所以进一步对月和年模拟结果进行流量数据的相对误差统计分析，结果如下表所示：

表10：江口流域月、年流量相对误差分析表(1981～1985)

从1981～1985年的逐月平均流量的季节性变化结合月平均流量图来看，在冬春季节，实际流量值比较小，1月～5月间模拟平均流量普遍比实测平均流量高，平均相对误差约为-35％左右；从6月开始降雨量明显增大，外部供水大量补充到土壤包气带中，实测流量接近模拟流量；在8月～11月间，模拟平均流量小于实测平均流量，平均相对误差约为+22％左右。因此ESSI模型对月流量过程模拟的规律为冬春季节模拟流量大于实测流量，夏秋季节模拟流量小于实测流量，在月平均流量图上表现出月模拟水文过程曲线的变化趋势比实测水文过程曲线平缓的状态，降雨的高峰期出现在6月，而径流的高峰期出现在7月，说明存在植被和土壤对外部供水的调节作用。

从年平均流量的年际变化来看，模型对流量大的年份模拟精度较高，1981年的平均相对误差仅为2.78％，而流量相对较小的年份的模拟精度稍低，平均相对误差最大为-26.44％，相对误差的正负分布与流量的大小无关。

ESSI模型在进行水文过程运算的同时，可以同步输出分布式的水文中间变量，不仅可以用于显示水文变量的时空变化规律，而且可用于辅助判断水文计算结果的可靠性。

1981年8月17日是1981～1985年日降雨径流过程模拟中的一个时间片断，选择这一天的模拟结果进行分析的原因在于这一天前后在江口流域内有强降雨过程，因此，各水文子过程模拟结果的变化及空间分布特征比较明显。

江口站的实测蒸散发经折算后为1.2mm，模拟值在1.83～1.886mm之间，在数值上，两者的数值大小接近，模拟结果的数量级正确；在空间分布形态上，河道附近高程低、降雨量大，区域蒸散发比较小，越往山区，随着高程的变大和空气对流的增强，蒸散发增大，这种变化形态与汉江上游区的蒸散发变化规律是一致的。

冠层截留的计算方法由于和流域的植被参数紧密相关，所以空间形态上表现出明显的经遥感影像解译得到的植被分布特征，河谷地区植被较稀疏，山区植被较茂密，冠层截留量分别占降雨量的5％和10％左右。从计算结果来分析，冠层截留量偏小，主要原因在于栅格逐日LAI和植被盖度的取值偏小。

与降雨量空间分布图对照来看，栅格地表产流类型呈以下分布特征：在降雨量大的区域容易发生超渗产流；靠近河道附近，地下水位高，包气带土壤缺水量容易满足，产流以超渗产流和蓄满产流同时发生为主；河道上游的大部分坡面不产流，这些区域降雨的绝大部分补充到土壤中。因此，整个流域的产流区集中在主河道附近，并逐渐向河道两侧的坡地发展，逐步过渡到不产流状态。

在广大的不产流区域，没有地表径流和壤中径流，但是有12～37mm的地下径流。在产流区域，由于降雨量和雨强大，产流量中超渗产流的比例大于蓄满产流，所以最大的地表径流产生在河谷附近超渗和蓄满产流同时发生的区域，径流量大于120mm，沿着河道向上游发展，地表径流量逐渐减小。壤中径流分布比较零散，主要与土壤的下渗能力和雨强有关。地下径流主要集中在流域的左侧支流附近，反映了不同类型土壤的蓄水能力的空间差异。

土壤时段下渗率的空间分布与地下径流的空间分布形态有一定的相似性，但是在数值大小上没有显示出明确的正相关或负相关，这是因为土壤时段下渗率还与土壤的前期含水量、土壤水的饱和程度以及外部供水状况等有关，是上述影响因素综合作用的结果。

土壤剖面含水量表示了当前时段土壤剖面深度内存储的水量的多少，土层平均含水量则反映了各土壤分层剖面含水量与土层厚度之比的加权平均。在分布形态上，河道附近含水量较大，而且越靠近流域出口，土层含水量越大，最大达到535mm左右。另外，从土壤剖面含水量的空间分布图上可以看到反距离加权法离散得到的降雨量空间分布的痕迹，说明外部供水对土壤含水量的变化起到重要的影响。

综上所述，分布形态和流域的客观规律相近，数值不存在数量级上的差异，说明通过水文过程模拟，ESSI模型不仅可以得到不同时间尺度下高精度的流量模拟结果，而且可以用来动态描述水文中间变量的时空变化规律。

Claims

1.一种以栅格为模拟单元的分布式水文模型的设计方法，主要包括以下步骤：

A)分布式参数获取；

B)流域水文过程概化；

C)产流过程设计；

D)汇流过程设计；

E)模型开发集成。

2.根据权利要求1所述的以栅格为模拟单元的分布式水文模型的设计方法，其特征在于步骤A)中，分布式参数有：降雨量、气温、太阳辐射、相对湿度、平均风速和潜在蒸散发，它们的获取是直接采用流域附近水文和气象站点的长期气象观测参数、利用气象发生器模拟所需的气象参数或通过气候学计算的方法，结合气候的季节性变化规律得到时序的气象参数；模拟的水文单元为栅格，对土壤参数的描述分为两个方向：水平方向按土壤类型划分；植被参数分为植被类型参数和植被生态参数两类。

3.根据权利要求2所述的分布式水文模型的设计方法，其特征在于步骤B)中根据霍顿下渗理论、蓄满产流理论和山坡水文学原理，建立栅格通用产流类型，如果前期土壤水已经蓄满，则在计算时段内采用土壤稳定下渗率来和雨强作对比，否则采用土壤时段下渗率来和雨强作对比，以决定是否发生超渗产流。

4.根据权利要求3所述的分布式水文模型的设计方法，其特征在于步骤C)中确定前期土壤水已蓄满条件下的地表径流计算、前期土壤水未蓄满条件下的地表径流计算、壤中径流和地下径流计算：取栅格为水文模拟单元，在高空间分辨率条件下，近似认为栅格内的土壤、植被特征是均一的，根据水量平衡理论，按不同的前期土壤水状态，可以分别计算出各部分的水量分配情况，计算过程中的变量定义如下：R_T-总产流；R_S1-超渗地表径流；R_S2-蓄满地表径流；R_L-壤中径流；R_G-地下径流；PE＝P-E；IF_C＝I-F_C；IF_P＝I-F_P；Δw-土壤含水量的变化量；

1)前期土壤水已蓄满条件下的地表径流计算：

产流类型①——不产流

R_T＝0，Δw＝PE＜0

R_S1＝R_S2＝0

R_L＝R_G＝0

产流类型②——超渗产流

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

R_T＝R_S，Δw＝PE-R_T

R_S1＝IF_C·Δt

R_S2＝0

R_L＝R_G＝0

产流类型③——蓄满产流

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

R_T＝PE，Δw＝0

R_S1＝0

R_G+R_L＝F_C·Δt

R_S2＝R_T-(R_G+R_L)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf

式中：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)

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

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

R_T＝PE，Δw＝0

R_S1＝IF_C·Δt

R_G+R_L＝F_C·Δt

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

R_G+R_L＝R_T-R_S1

R_S2＝0

否则，有：

R_S2＝R_T-R_S1-F_C·Δt

如果R_G+R_L＜0，则取

R_S1＝R_T·d_Surf

R_G+R_L＝R_T·(1-d_Surf)

2)前期土壤水未蓄满条件下的地表径流计算

产流类型⑤——不产流

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

R_T＝0，Δw＝PE

R_S1＝R_S2＝0

R_L＝R_G＝0

如果PE＜SP_Sw-SP_Wp，则：

Δw＝SP_Sw-SP_Wp

产流类型⑥——超渗产流

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

R_T＝R_S，Δw＝PE-R_T

R_S1＝IF_P·Δt

R_S2＝0

R_L＝R_G＝0

产流类型⑦——蓄满产流

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

R_T＝PE-D，Δw＝D

R_S1＝0

R_G+R_L＝F_P·Δt

R_S2＝R_T-(R_G+R_L)

如果R_S2＜0，则取

R_S2＝R_T·d_Surf

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

R_G+R_L＝R_T·(1-d_Surf)

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

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

R_T＝PE-D，Δw＝D

R_S1＝IF_P·Δt

R_G+R_L＝F_P·Δt

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

R_G+R_L＝R_T-R_S1

R_S2＝0

否则，有：

R_S2＝R_T-R_S1-F_P·Δt

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

R_S1＝R_T·d_Surf

R_G+R_L＝R_T·(1-d_Surf)

壤中径流和地下径流计算如下：

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

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

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

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

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

{Δt}_{p} = \frac{(SP_Sat - SP_Fc)}{SP_Sat_K}

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

R_L＝W_P·d_L

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

W′_P＝F_P·Δt

W′_P＝SP_Sw-SP_Wp

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

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

当计算出壤中径流后，由壤中径流和地下径流之和减去壤中径流，就可以分离出地下径流量；

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

5.根据权利要求4所述的分布式水文模型的设计方法，其特征在于步骤D)中在模型中分别设计马斯京根-康吉法、滞时演算法和河道分段马斯京根法，用于不同条件下模型的汇流过程演算。

6.根据权利要求2、3、4或5所述的分布式水文模型的设计方法，其特征在于步骤E)中对模型的体系架构和开发集成技术如下：

1)模型的开发框架为：

利用ODBC方式实现属性参数和驱动参数的数据库访问；

2)模型分为三个部分：

栅格水文模拟总控部分由一系列动态库组成，包括气象参数模拟及空间离散动态库、栅格河道信息提取动态库、栅格水文子过程模拟动态库、栅格水量平衡计算动态库和栅格河道汇流处理动态库；

3)围绕模型的栅格水文模拟总控需求，定义了六个主体类：

4)模型采用面向对象方式进行开发，对所有的水文子过程进行封装和模块化处理，通过界面定制的方法，人为选择模拟区域最佳的水文处理模块组合；

5)模型的开发集成模式采用了以水文模型为主，将GIS功能嵌入水文模型中紧密集成开发的模式。