CN103093114B - 一种基于地形和土壤特性的分布式流域缺水量测算方法 - Google Patents

Abstract

本发明属于水利工程的水循环技术领域，特别是涉及一种基于地形和土壤特性的分布式流域缺水量测算方法。该方法包括步骤一流域下垫面基本数据处理、步骤二核算流域范围内的栅格地形指数、步骤三计算流域范围内的土壤特性参数、步骤四定量表达流域缺水量的分布和步骤五绘制流域缺水量累积分布曲线，从而得到分布式流域缺水量测算模型。本发明所得模型应用于实时洪水预报将可大幅度提高流域洪水预报的精度，为流域防洪预警等提供可靠的科学支撑。本发明既适用于集总式水文模型，又适用于分布式水文模型，将有效地推动流域水文科学研究的深入发展。

Description

一种基于地形和土壤特性的分布式流域缺水量测算方法

技术领域

本发明属于水利工程的水循环技术领域，特别是涉及一种基于地形和土壤特性的分布式流域缺水量测算方法。

背景技术

自然界水循环过程对人类经济活动及社会活动都有极其重要的影响，水循环活动造成的洪灾、旱灾以及水资源的分配不均给全世界造成了巨大损失。为了能够了解水循环的过程，进而通过科学手段尽量避免水旱灾害，水文学对水循环过程做出了很多深入研究。现代水文学通过实验、物理过程分析总结水循环规律，并在此基础上结合数学物理方法建立能描述水循环过程的水文学模型，期待通过模型能准确地了解水循环的各相关过程。从二十世纪50年代到现代，水文学模型研究越来越趋近于理论化，模型从最初的统计模型，逐步过渡到半物理的概念性模型以及全动力过程描述的物理模型，通过这些模型的研究有力地推动了人类对水文循环过程的了解，对指定趋利避害的规则提供了大量的科学依据。目前水文模型主要的研究思路是分别研究水文产流和汇流，其中产流过程是核心，它决定了降雨后形成的洪水总量。而在产流过程中流域土壤的缺水量是决定性因素，降雨扣除缺土壤水量后即可得到产流量。基于此，当前的水文模型均引入了土壤缺水量相关的参数，蓄水容量是反映土壤缺水量的有效表达。

流域地形是确定土壤缺水量的重要因子，我国著名的新安江模型提出了反映土壤缺水量空间分布的蓄水容量曲线，隐含地表达了地形变化和土壤下垫面变化对流域蓄水容量的影响；BeVen1979年提出的TOPModel是一个以地形为基础的半分布式水文模型，其主要特征是利用流域地形指数(ln(α/tanβ))作为一个水文相似性的指数，所有具有相同地形指数值的点，对降雨具有相同的水文响应。；VIC模型也提出了类似于流域缺水量计算方法，隐含了地形的影响。流域土壤特征是影响缺水量的另外一个关键指标，土壤性质决定了土壤的含水量范围，进而影响其缺水量，Troch就从理论上推导出了土壤饱和缺水量和土壤势能之间的关系表达式。

虽然目前有关于地形对流域土壤缺水量影响的研究，但研究的成果多采用表现为经验性曲线来描述，比如新安江模型提出了一个抛物型的蓄水容量曲线，其中有经验性的参数需要率定，模型结构理论化程度不高。缺水量本质上与流域地形和土壤质地密切相关，但目前还未见结合地形和土壤特性来共同确定土壤缺水量的报道。

发明内容

本发明的目的是为克服现有技术的不足而提供一种基于地形和土壤特性的分布式流域缺水量测算方法，可通过建立地理地形数据及土壤类型信息与流域栅格土壤缺水量之间的关系，实现在空间尺度上流域缺水量指标动态设定并与实际情况的吻合，从而得到分布式流域缺水量测算模型，将有效地推动流域水文模型的科学发展。

根据本发明提出的一种基于地形和土壤特性的分布式流域缺水量测算方法，其特征在于包括以下步骤：

步骤一，流域下垫面基本数据处理：地理数据采集、计算栅格指标和核算栅格土壤信息；

步骤二，核算流域范围内的栅格地形指数：一是根据流域高程栅格，分别计算栅格坡度角、积水面积，获得任意栅格的地形指数；二是确立流域栅格内地形指数，获得全流域地形指数分布；

步骤三，计算流域范围内的土壤特性参数：一是根据流域栅格土壤信息，分别计算每个栅格的土壤类型及其所占比重；二是确立流域栅格内土壤特性参数信息，建立土壤类型与土壤水分指标的映射关系；

步骤四，定量表达流域缺水量的分布：计算栅格型土壤缺水量指标；

步骤五，绘制流域缺水量累积分布曲线：根据步骤四的结果，计算各子流域内栅格点的缺水量，按照从小到大的顺序排列，再根据栅格面积比绘制各子流域离散化的地形指数累积分布曲线，即流域的缺水量分布曲线，从而得到分布式流域缺水量测算模型。

本发明进一步的优选方案为：

步骤一中所述的地理数据采集，包括收集流域地形数据、栅格高程数据、栅格土壤特性数据；所述的计算栅格指标，包括计算栅格内的坡度、坡向以及河长,划分流域边界,生成子流域拓扑关系与流域水系；所述的核算栅格土壤信息，包括计算栅格内各类土壤的比例，叠加栅格内各类土壤的组成。

步骤四中所述的计算栅格型土壤缺水量指标，过程包括：

（一）利用无降雨补给情况下土壤水分稳定状态分布曲线(以VG模型描述)，积分获得达到田间持水量状态的土壤缺水量：

（二）利用空间尺度转化方法将土壤参数栅格化，结合实验方式获得各类土壤栅格内的土壤参数，通过这种方式获得式(1)中的θ_s、θ_r、Ψ_c、θ_f、α、n等参数；将式(1)中的地下水埋深D与地形指数通过建立如下关联式(2)～(4)：

式(2)中：D为地下水埋深；S_zm为参数；为流域地下水平均埋深；为流域平均地形指数；tp_i为地形指数；

（三）将式(4)代入式(1)计算栅格型土壤缺水量指标，该指标由流域内各栅格点缺水量与地形因子及土壤特征的相关表达式：

流域中各栅格点地形指数不同，土壤特征各异，通过式(5)得到每个栅格点上相应的缺水量指标。

本发明与现有技术相比的显著优点为：一是本发明能够实际计算出流域中各点的缺水量，反映流域内缺水量的空间分布特征，便于在分布式模型中应用；二是本发明能够考虑流域中的地形和土壤的空间分布，下垫面的非均质性将导致不同的缺水量曲线；三是计算参数少，只需要给出流域的平均蓄水容量WM，即可计算出参数S_zm，然后就可确定缺水量曲线，而WM的给定一般有规律可循；四是在基于现代测量技术所得的前期资料支持下，本发明所得模型应用于实时洪水预报将可大幅度提高流域洪水预报的精度，为流域防洪预警等提供可靠的科学支撑。本发明既适用于集总式水文模型，又适用于分布式水文模型，将有效地推动流域水文科学研究的深入发展。

附图说明

图1本发明提出的分布式流域缺水量测算方法的流程示意图。

图2寻找洼地区域或平坦区域出口的方法示意图。

图3最短流程算法示意图。

图4格网点水流方向编码示意图。

图5一个简单的DEM及其计算结果示意图。

图6DEM水文分析流程示意图。

图7无降雨补给情况下土壤水分稳定状态分布曲线示意图。

图8紫罗山流域水流方向矩阵示意图。

图9紫罗山流域累计矩阵形式示意图。

图10紫罗山流域水系示意图。

图11紫罗山子流域分割及流域边界示意图。

图12紫罗山子流域栅格坡度示意图。

图13紫罗山流域内栅格地形指数示意图。

图14流域缺水量与地形指数的关系示意图。

图15不同土壤类型缺水量曲线示意图。

具体实施方式

下面结合附图和具体应用实施例对本发明作进一步的详细说明。

结合图1，本发明提出的一种基于地形和土壤特性的分布式流域缺水量测算方法，总体思路是以GIS(Geographic information system，地理信息系统)生成数字流域水系为基础，以分割流域地形指数栅格和土壤特征参数栅格为技术手段，结合程序分析各子流域土壤和地形指数分布，提取各子流域的缺水量分布，最终构建流域的缺水量测算模型。实施本发明的具体步骤包括：

步骤一，流域下垫面基本数据处理。具体是以DEM(Digital Elevation Model,数值高程模型)数据为基础，计算栅格内坡度、坡向以及河长，划分流域边界，生成子流域拓扑关系与流域水系，其计算遵循以下顺序进行：

(一)洼地处理

被高程较高的区域围绕的洼地是应用数字高程流域水系生成模型生成流域水系的一大障碍，因为在决定水流方向以前，必须先将洼地填充。有些洼地是DEM生成过程中带来的数据错误，但另外一些却表示了真实的地形。一些研究试图通过平滑处理消除洼地，但平滑处理只能处理浅的洼地，较深的洼地仍然无法处理。处理洼地的另一种方法是通过将洼地的每一格网点赋以洼地边缘最小值，从而达到消除洼地的目的。这些方法都需要对地形数据进行修改，可能会产生一些不合理的方向阵，下面将具体的介绍此种情况的处理。

根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域高程，使研究区域中的水流能够通过洼地区域和平坦区域。算法的计算顺序如下：

①结合图2，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找，直到找到这个最低格网；

②结合图3，如果边缘格网中高程最低的格网点不为洼地区域内的格网点，则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

③结合图3，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的高程。

(二)水流方向矩阵的计算

结合图4，采用D8(Deterministic Eight-neighbours)算法进行计算。以图5(a)的原始DEM矩阵为例，每个网格值为格网点高程值，将中心格网点同其最邻近的8个格网点之间的坡降进行比较，其中与落差最大的一个格网点中心之间连线的方向便定义为中心格网点的水流方向，并且规定一个格网点的水流方向用一个特征码表示。有效的水流方向定义为东、东南、南、西南、西、西北、北和东北，并分别用1、2、3、4、5、6、7和8这8个有效特征码表示，表示方法见图4。

中心格网点同相邻8个格网点之间单位距离的高差为：

MD=Z/D (6)

式(6)中：MD为两个格网点之间的单位距离的高差，表示地形坡度；Z为两个格网点之间的高程差；D为两个格网点中心之间的距离。

确定水流方向的具体计算顺序如下：

①对所有DEM边缘的格网，赋以指向边缘的方向值0；

②对所有在第一步中未赋方向值的格网，计算其对8个邻域格网的单位距离的高差值，确定具有最大落差值的格网，执行以下顺序：

A.如果该格网与相邻8个邻域格网的最大落差小于0，则赋以负值以表明此格网方向未定（这种情况在经洼地处理的DEM中不会出现）。

B.如果该格网的高程与相邻8个邻域格网最大落差大于或等于0，且最大落差只有一个，则该格网的水流方向赋以指向最大落差的方向。

C.如果该格网的高程与相邻8个邻域格点最大落差大于0，且最大落差有多个，则该格网的水流方向在逻辑上以查表方式确定，也就是说，如果中心格网在一条边上的三个邻域点有相同的落差，则中间的格网方向被作为中心格网的水流方向，如果中心格网的相对边上有两个邻域格网落差相同，则任选一格网方向作为水流方向。

D.如果该格网的高程与相邻8个邻域格点最大落差等于0，且最大落差有多个，则以这些0值所对应的方向值相加。

③对①、②步没有赋以负值，0，1，2，3，…，8的每一格网，检查对中心格网有最大落差值的邻域格网。如果邻域格网的水流方向值为1，2，3，…，8，且此方向没有指向中心格网，则以此格网的方向值作为中心格网的方向值；

④重复③，直到所有格网都被赋值，得到水流向矩阵如图5(b)所示。

(三)水流累计矩阵的计算

结合图5，区域水流累计矩阵表示区域地形每一点的流水累积量。其基本思想是，假定以规则格网表示的数字地面高程模型的每一点有一个单位的水量，按照水流从高往低流的规律，根据区域地形的水流方向矩阵，计算每一格网点流过的水量数值，便可得到该区域的水流累计矩阵。下面给出一个从原始DEM矩阵计算出相应水流方向矩阵及水流累计矩阵的范例。

以图5(a)到图5(b)的变换方法为基础，下面简要说明从图5(b)到图5(c)的变换方法：

以图5(b)左上单元格标注为“2”的格网点为例，因为周围没有那个单元格的水流流入该单元格（即该单元格为流域边界），所以在图5(c)中相应填入0，同理可知图5(c)中第一行所填的数均为0；以图5(b)中第二行第二单元格为例，由于其左上方单元格的水流将流入此单元格，因此，在图5(c)中相应为指填入1；以图5(b)中第二行第四列单元格为例，其左上方和正上方单元格的水流都将流入其中，即有两个单位的水流汇入此单元格，因此，在图5(c)中相应的单元格填入2；再以图5(b)中第三行第四列单元格为例，其左上方和正上方的水流将汇入其中；另外，左上方已经有一个单位的入流，正上方单元格则有两个单位的入流，再加上两个单元格自身的水量，此时，图5(c)中相应单元格所填入的数应该是5。以此类推，即可由图5(b)得到图5(c)，从而生成水流累积矩阵。

(四)基本水文分析，其中包括汇水面积计算，流域分水线识别，河网生成，建立流域水系拓扑关系等四部分。

①汇水面积计算

按最大落差规则确定的水流路径可以非常方便地计算出指定格网点以上的流域汇水面积，若以格网数表示面积的话，则该格网点的汇水面积数值是该格网点以上汇入该格网点的格网点数目。在算法上可利用递归算法实现，由指定点开始按逆水流方向向上搜索迭代，即可得到该集水流域内任意格网点的汇水面积，其结果如图5(c)表示的水流累积矩阵。

②流域分水线的识别

给定流域的主要入口和出口的断面位置，即它们所在栅格单元的行列坐标。一旦两者的位置确定，根据它们的汇水面积大小，程序可自动搜索从而勾划出流域边界并计算出流域面积。

③河网生成

河网生成分为三个步：确定流域边界内的水道；裁减小于某一临界长度的河段；生成河道编码。

首先，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点即定义为水道。

其次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由第一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道、位于河谷两边的凹痕或沟壑的出口，需要将它们裁减除去。

最后，确定河道级数和河段长度。根据流域出口断面确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的编码。同时可以确定各级支流汇入上一级河道的节点，对河网所有节点进行编码，这样定义的节点编码可用于水文汇流计算或河网数据库的构建。

④流域水系拓扑关系

一旦生成联结完好的河网，根据各河网节点，可以确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息。一方面，河网与子流域边界等空间信息是以栅格形式存储，这样易于用GIS软件作可视化显示；另一方面，河段或子流域的拓扑关系还以表格形式存储，有利于数字水文模型的调用。详见图6利用DEM进行水文分析的流程示意图。

步骤二，核算流域范围内的栅格地形指数。具体是根据流域高程栅格，分别计算栅格坡度角、积水面积，获得任意栅格的地形指数；通过确立流域栅格内地形指数，获得全流域地形指数分布，其计算遵循以下顺序进行：

(一）计算栅格上游集水面积

根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域水流方向矩阵，使研究区域中的水流能够通过洼地区域和平坦区域，算法的计算顺序如下：

①结合图2，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找，直到找到这个最低格网；

②结合图3，如果边缘格网中高程最低的格网点不为洼地区域内的格网点,则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

③结合图3，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的水流方向阵；

④根据区域地形的水流方向矩阵，计算每一格网点流过的水量数值，得到该区域的水流累积矩阵；

⑤由流域内任一栅格出发，向其上游统计所有汇入该栅格的栅格面积，记做积水面积，计算积水面积与所属网格周长之比，得到流经坡面任一点处单位等高线长度的汇流面积。

(二）计算栅格地形坡度

根据累计矩阵，确定由任一栅格向四周八个相邻栅格中，坡度最陡的方向，以此两栅格高度差与水平距离之比，得到栅格地形坡度数据。

(三）计算栅格地形指数

根据地形指数公式：

式(7)中：T_i为地形指数；a_i为集水区面积；β_i为坡度角；i为栅格编号。

由上述步骤(一）中⑤得到栅格集水区面积，再由步骤(二）得到栅格地形坡度，将上述二者代入到式(7)中，计算得到流域内各栅格上的地形指数T_i。

步骤三，计算流域范围内的土壤特性参数。具体是指根据流域栅格土壤信息，分别计算每个栅格的土壤类型及其所占比重；二是确立流域栅格内土壤特性参数信息，建立土壤类型与土壤水分指标的映射关系，其计算遵循以下顺序进行：

(一）土壤性质参数收集

搜集流域范围内土壤的相关资料，包括土壤质地、土壤类型、土壤饱和含水量，土壤孔隙度、土壤导水性质等方面数据，其中如土壤类型及其所占比重的数据是以栅格形式表达。

(二）栅格尺寸同化

由于土壤质地类型及其所占比例等数据一般是按照栅格形式表达，需要将该数据的栅格尺寸与流域地形指数栅格的尺寸进行匹配，以数据同化思路出发，以地形指数栅格为目标，采用线性处理方式重新划分土壤类型比例，完成土壤类型栅格与地形指数栅格尺寸匹配的工作。

(三）土壤性质多层栅格数据叠加

针对研究区同一栅格范围内存在多种土壤类型，且各成分所占比例不同，将多种土壤类型性状参数数据与其对应的成分比例在同一栅格内叠加，得到栅格内综合的土壤性质参数。

步骤四，定量表达流域缺水量的分布。具体是指以栅格为单位，结合地形指数以及土壤特征参数信息，建立包含地形指数和土壤特征参数的土壤缺水量数学模型，计算栅格型土壤缺水量指标，其计算遵循以下顺序进行：

(一）结合图7，利用无降雨补给情况下土壤水分稳定状态分布曲线，积分获得达到田间持水量状态的土壤缺水量式：

式(8)中：Ψ为土水势对应的吸力；Ψ_c为土壤对应的田间持水量时毛管上升高度；θ(Ψ)为土壤体积含水率；θ_s为土壤饱和含水率；θ_f为土壤田间持水量；α、n为土壤水参数；θ_r为土壤凋萎含水率；D为地下水埋深/mm。

将式(8)中各土壤水分参数（如土壤饱和含水率θ_s、土壤凋萎含水率θ_r、土壤水参数α和n等）及地形参数（如地下水埋深D）与栅格地形指数及土壤特征建立联系。

(二）建立地下水埋深与地形指数的关系

式(9)中：D为地下水埋深；S_zm为参数；为流域地下水平均埋深；为流域平均地形指数；tp_i为地形指数。

将式(9)～式(11)带入到式(8)中，将原有公式中的地下水埋深用地形指数替代，将栅格地形指数信息引入到土壤缺水量计算中。

(三）非饱和带土壤水分状态与土壤势能的关系：

式(12)中：Ψ为土壤吸力，表示土壤水势能；θ_r为土壤凋萎含水量；θ_s为土壤饱和含水量。

结合实验测得多种土壤对应的模型参数如表1所述。将表1中相关参数与同化后的栅格土壤类型建立映射关系，结合其所占比例，将原有式(8)中的土壤饱和含水率、土壤凋萎含水率、土壤水参数等用不同质地土壤类型相关参数替代，将栅格土壤类型信息引入到土壤缺水量计算中。

表1不同土壤类型对应的模型参数列表

(四）计算栅格型土壤缺水量指标，由流域内各栅格点缺水量与地形因子及土壤特征的相关表达：

流域中各栅格点地形指数不同，土壤特征各异，通过式(13)得到每个栅格点上相应的缺水量指标。

步骤五，绘制流域缺水量累积分布曲线。具体是指结合流域分割，得到各子流域内栅格缺水量序列，根据栅格面积比绘制各子流域离散化的地形指数累积分布曲线，即流域的缺水量分布曲线，其计算遵循以下顺序进行：

(一）结合式(13)，栅格地形指数以及土壤类型参数，计算流域范围内各栅格的缺水量；

(二）结合流域水系构成、子流域划分范围，以及各子流域的拓扑关系，划分各子流域内所组成的栅格；

(三）统计各子流域内栅格缺水量，按照从小到大顺序排列，根据栅格面积比建立各子流域离散化的缺水量累积分布曲线，得到流域内的缺水量空间分布,最终得到分布式流域缺水量测算模型。

下面以淮河流域紫罗山水文站控制流域为应用实施区域(下称紫罗山流域)实例，进一步对本发明的具体实施方式作如下详细说明。

紫罗山水文站位于淮河流域沙颍河支流北汝河上游，流域控制面积约1800km²，流域处于山丘区，最高点海拔2122m，最低点海拔284m，流域平均高程约820m，地势呈西南高东北低的对角陡坡形态。流域属大陆性气候，多年平均降雨900mm，60～70%的雨水集中在6、7、8三个月，紫罗山水文站设立于1951年4月，具有长序列的水位、流量观测资料，水文站上游没有大中型水利工程，受人类活动影响较小，保持了较好的自然流域状态。

选取紫罗山流域为应用实施区域，利用本发明的全流域分布式缺水量模型，结合DEM数据驱动流域水系生成、子流域范围分割以及地形指数算法，利用紫罗山流域的水文气象数据，土壤类型栅格数据和数字高程数据等，驱动全流域分布式缺水量模型，模拟紫罗山流域水文过程。

步骤一：流域下垫面基本数据处理

(一)采集紫罗山流域DEM资料

紫罗山流域，经度111°56′53″～112°31′30″，纬度34°13′52″～33°41′22″。由公开的30米分辨率的ASTGTM数据，采集研究区范围内地形栅格数据；

(二)洼地处理

在处理水流方向之前，根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域高程，使研究区域中的水流能够通过洼地区域和平坦区域，按如下顺序执行洼地处理：

首先，结合图2，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找,直到找到这个最低格网；

其次，结合图3，如果边缘格网中高程最低的格网点不为洼地区域内的格网点,则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

最后，结合图3，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的高程。

(三)水流方向矩阵计算

采用D8算法按照如下顺序进行计算：

首先，对所有DEM边缘的格网，赋以指向边缘的方向值0；

其次，对所有在前一步中未赋方向值的格网，计算其对8个邻域格网的单位距离的高差值，确定具有最大落差值的格网；

再次，对前两步没有赋以负值，0，1，2，3，…，8的每一格网，检查对中心格网有最大落差值的邻域格网。如果邻域格网的水流方向值为1，2，3，…，8，且此方向没有指向中心格网，则以此格网的方向值作为中心格网的方向值；

最后，重复上述再次的过程，直到所有格网都被赋值，最终得到紫罗山流域水流方向矩阵形式，详见图8所示。

(四)水流累计矩阵计算

根据紫罗山流域栅格水流方向阵计算结果，计算相应的水流累计矩阵，其计算顺序如下：

以图5(b)左上单元格标注为“2”的格网点为例，因为周围没有那个单元格的水流流入该单元格（即该单元格为流域边界），所以在图5(c)中相应填入0，同理可知图5(c)中第一行所填的数均为0；以图5(b)中第二行第二单元格为例，由于其左上方单元格的水流将流入此单元格，因此，在图5(c)中相应为指填入1；以图5(b)中第二行第四列单元格为例，其左上方和正上方单元格的水流都将流入其中，即有两个单位的水流汇入此单元格，因此，在图5(c)中相应的单元格填入2；再以图5(b)中第三行第四列单元格为例，其左上方和正上方的水流将汇入其中；另外，左上方已经有一个单位的入流，正上方单元格则有两个单位的入流，再加上两个单元格自身的水量，此时，图5(c)中相应单元格所填入的数应该是5。以此类推，即可由图5(b)得到图5(c)，从而生成水流累积矩阵，最终得到紫罗山流域累计矩阵形式如图9所示。

(五)流域水系生成

首先，给定紫罗山流域的主要入口和出口的断面位置（经度：112°30′23″，纬度：34°11′29″）。根据流域的汇水面积大小，按照流域累计矩阵搜索从而勾划出流域水系；

其次，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点即定义为水道。

再次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由前一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道、位于河谷两边的凹痕或沟壑的出口，需要将它们裁减除去。

最后，确定河道级数：四级，河段长度：1497～42275m不等。根据流域出口断面—紫罗山水文站，确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的拓扑关系，紫罗山流域水系，如图10所示。

(六)流域水系拓扑关系

生成联结完好的河网水系，根据河网节点，确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息，从而勾划出流域边界并计算出流域面积，紫罗山流域共划分出50个子流域，如图11所示。

步骤二：核算流域范围内的栅格地形指数

以DEM数据分析过程中的填洼处理、累计矩阵计算、栅格坡度以及上游汇流面积计算等为基础，计算栅格地形指数，其计算遵循以下顺序：

(一)计算栅格上游集水面积

首先，结合图2，先找出洼地区域或平坦区域的边缘格网，然后在边缘格网中找出高程最低的格网点。如果不存在最低格网，则将洼地区域或平坦区域向外围扩大一个格网继续寻找，直到找到这个最低格网。结合图3，如果边缘格网中高程最低的格网点不为洼地区域内的格网点,则将边缘格网中高程最低格网点作为洼地区域或平坦区域的水流出口点；

其次，结合图3，以洼地区域或平坦区域的水流出口点为起点，修改洼地区域或平坦区域内高程低于水流出口格网点的水流方向阵。

再次，根据区域地形的水流方向矩阵，计算每一格网点流过的水量数值，得到该区域的水流累积矩阵。

最后，由流域内任一栅格出发，向其上游统计所有汇入该栅格的栅格面积，记做积水面积，计算积水面积与所属网格周长之比，得到流经坡面任一点处单位等高线长度的汇流面积。

(二)计算栅格地形坡度

结合图12，根据累计矩阵，确定由任一栅格向四周八个相邻栅格中，坡度最陡的方向，以此两栅格高度差与水平距离之比，得到紫罗山流域栅格地形坡度数据。

(三)计算栅格地形指数

根据地形指数公式：

式(14)中：T_i为地形指数；a_i为集水区面积；β_i为坡度角；i为栅格编号。

由步骤(一)中最后一步得到栅格集水区面积，再由步骤(二)得到栅格地形坡度，将上述二者代入到式(14)中，计算得到紫罗山流域内各栅格上的地形指数T_i，详见图13。

步骤三：计算流域范围内的土壤特性参数

结合各种不同尺度的土壤栅格信息，同化栅格内土壤数据，核算土壤特征参数，其计算遵循以下顺序：

首先,收集流域土壤数据，结合中科院土壤所公布的我国1:100万的土壤数据库，将将土壤主要分为黏土、壤土和砂土三类；

其次，结合流域土壤数据30m分辨率的栅格数据信息，将土壤类型以及其所占比重叠加到同一范围栅格内；

最后，结合Troch等（TROCH P.Conceptual basin-scale runoff process models forhumid catchments:analysis,synthesis and applications,Ph.D.thesis.GhentUniversity,1992）通过实验多种土壤得到的土壤参数信息如表2、表3所示，将土壤栅格数据转化为相应土壤水分参数。

表2不同土壤类型对应的模型参数列表

表3毛管水上升高度表

步骤四：定量表达流域缺水量的分布

(一)基本原理

在存在饱和地下水面时，Van-Genuchten（VAN GENUCHTEN M T,1980.A closed formequation for predicting hydraulic conductivity in unsaturated soils.Soil Sci.Soc.Am.J.44,892-898）模型（以下简称VG模型）是描述非饱和带中土壤水分特征时的常用模型，Troch(TROCH P.Conceptual basin-scale runoff process models for humid catchments:analysis,synthesis and applications,Ph.D.thesis.Ghent University,1992）得出如下改进式：

式(15)中：Ψ为土水势对应的吸力；θ(ψ)为土壤体积含水率；θ_s为土壤饱和含水率；α、n为土壤水参数。

用图形表示VG模型在静态无降雨补给情况下的土壤水分分布特征，见图7，图中有三个关键含水率：凋萎含水率θ_r、田间持水率θ_f以及饱和含水率θ_s，其中含水率为θ_s时位于地下水位面附近；含水率为θ_r时一般位于地表附近，但在地下水埋深较浅时，土壤中含水率可能均大于θ_r；含水率为θ_f时与VG曲线的交点(θ_f，Ψ_c)即是由于毛管力维持达到田间持水量时的最大高度，Ψ_c以上直到地表含水率均小于田间持水率，其计算公式可以从式(15)中得出如下式：

图7中非饱和带土壤的厚度为D，假设土壤中水分分布符合VG模型，其土壤含水量可以通过式(15)空间积分计算：

土壤吸力Ψ与土壤基质势对应，可以定义为地下水面上的高度，取Ψ=z，基于该定义式(17)积分的最终表达式为：

当地下水埋深为D时，M(D)可以认为是非饱和带的稳定最低含水量，整个土层中需要达到饱和时缺水量为：

由图7可见，地下水位以上Ψ_c的高度内土壤已经达到田间持水量，Ψ_c~D范围内土壤含水量可以计算为：

深度为D的土层中要满足田间持水量时的缺水量可以通过式(20)推算为：

(二)土壤因素影响下得栅格缺水量定量表达

针对缺水量计算公式(21)中土壤参数指标，如饱和含水率θ_s，凋萎含水率θ_r毛管水上升高度Ψ_c及土壤水参数等，均涉及土壤类型。

借助中科院土壤所公布的我国1:100万的土壤数据库栅格土壤类型数据，及其所占比例，结合Troch等（TROCH P.Conceptual basin-scale runoff process models for humidcatchments:analysis,synthesis and applications,Ph.D.thesis.Ghent University,1992）通过实验获得的各种土壤对应的VG模型参数以及Cosby等（COSBY B J,HORNBERGER GM,et al.Estimating catchment water quality response to acid deposition usingmathematical modelsof soil ion exchange processes.1986,38:77-95）对不同土壤饱和含水率与田间含水率值分析结果详见表2和表3，将栅格内缺水量公式中的上述参数按土壤性状确定。

(三)地形指数因素影响下得栅格缺水量定量表达

地形是影响土壤缺水量的另外一个关键因素，根据TOPMODEL的假设，地下水埋深与地形指数存在下列线性关系（BEVEN K J,KIRKBY M J.A physically based variable contributingarea model of basin hydrology.Hydrological Science Bulletin,1979,24(1):43-69）：

式(22)中：D为地下水埋深；S_zm为参数；为流域地下水平均埋深；为流域平均地形指数；tp_i为地形指数。

将式(22)～式(24)代入到式(21)中，得到流域各点与地形因子相关的缺水量表达式：

式(25)为缺水量的定量模型，其中的地形指数反映了地形、地貌对缺水量的影响，而含水率等水分特征参数则反映了土壤的影响，由于流域中各栅格点地形指数不同，栅格内土壤特征也各异。

(四)参数确定

式(25)中的地形指数和土壤水分特征参数可以通过计算或查表取得，而其中C和S_zm仍然是未知量，需要提出相应的确定方法。在单一土壤分布情况下，饱和含水量式(19)和缺水量式(21)均与地下水埋深D正相关，而式(22)～式(24)表明D_i与地形指数tp_i负相关，因此式(25)中的W_f(D)_i与地形指数负相关，这与TOPMODEL中的基本原理一致。在靠近流域下游河槽处地形指数最大的地方，地下水常出露地表形成回归流，地下水埋深为零，依据这一设想，在最大地形指数时，根据式(22)～式(24)得出C：

(26)

将式(26)代入到式(25)可以得到单点的缺水量计算公式：

式(27)中只有S_zm一个自变量，在文献（林凯荣,郭生练,熊立华等.DEM栅格分辨率对TOPMODEL模拟不确定性的影响研究.自然资源学报,2010,25(6):1022-1032）中指出S_zm的取值范围是1～100cm。当给定流域的平均蓄水容量WM后，首先假设S_zm=0cm，由式(25)求取流域内各栅格点的缺水量，进而可以得出一个初步的流域平均蓄水容量WM0，当WM0与WM存在差距时通过调整S_zm再一次计算，反复迭代直到得出的WM0与给定WM在误差允许范围内即可。此时的S_zm即为流域的有效参数，同步得出各栅格内缺水量W_f。

步骤五：绘制流域缺水量累积分布曲线

结合流域分割，得到各子流域内栅格缺水量序列，根据栅格面积比绘制各子流域离散化的地形指数累积分布曲线，其计算遵循以下顺序进行：

(一）结合式(13)，栅格地形指数以及土壤类型参数，计算流域范围内各栅格的缺水量；

(二）结合流域水系构成、子流域划分范围，以及各子流域的拓扑关系，划分各子流域内所属栅格；

(三)统计各子流域内栅格缺水量，按照从小到大顺序排列，根据栅格面积比建立各子流域离散化的缺水量累积分布曲线，得到流域内的缺水量空间分布,最终得到分布式流域缺水量测算模型。

紫罗山流域共划分50个子流域，将其中得到的三个缺水量曲线差别较大的区域，即图11中的第1、2、3子流域，其地形指数对子流域缺水量曲线的影响规律如图14所示。

采用所提方法计算出三个子流域的缺水量曲线见图14a所示，对应子流域的IRDG曲线如图14b所示，两者虽然在形状上不同，但趋势具有一致性，如IRDG曲线中第1子流域曲线位于最上部，坡度最缓，对应的缺水量曲线中第1子流域曲线也坡度最缓，在三个子流域平均蓄水容量WM相同情况下，第1子流域缺水量曲线表现为更明显的下凹趋势，此时流域缺水量分布更均匀。研究区内所有子流域的缺水量曲线与IRDG曲线都具有相同趋势。

土壤类型对流域缺水量曲线也有较大影响，以第1子流域为例，目前已经获得了实际土壤分布状态下的缺水量曲线，另外假设该子流域的土壤类型单一，分别为黏土、壤土以及砂土，在此根据式(27)得出相应的缺水量曲线，几条曲线对比如图15，针对同一个子流域，当平均蓄水容量WM相同，且砂土质土壤情况下流域缺水量分布最均匀，黏土质时缺水量分布差异最大，壤土质时介于二者中间。

结合缺水量指标，针对紫罗山流域进行径流模拟，选取1970～1979年资料率定，1980～1987年汛期水文资料用于模型验证，以紫罗山水文站流量过程为目标，分别对各场洪水的洪量、洪峰以及确定性系数等进行统计，得到的结果见表4，在12场洪水过程中，洪量相对误差绝对值最大为9.53%，洪峰相对误差绝对值最大为5.77%，洪水过程的流量误差均方根在4.62～67.7m3/s之间，有11场确定性系数在0.9以上，其中一场在0.78，分析原因是由于本次降雨在流域上分布极不均匀，而紫罗山流域内的雨量站个数很有限，尤其在流域上游几乎没有雨量站分布，雨量信息采集的不均匀性导致了径流模拟的较大误差。本次模拟12场洪水的结果符合《水文情报预报规范》的甲级方案标准。

表4模拟成果统计表

验证结果表明根据地形和土壤特征定量确定流域缺水量的方法，该方法不仅能计算流域各栅格点的缺水量，而且可以用于确定流域的缺水量曲线。实现将单点的缺水量确定方法成功应用于分布式水文模型中，提高产流模型的精度；并将新安江模型的缺水量曲线从经验性曲线变为了定量曲线，使其具有物理意义，减少模型参数。通过对地形指数、土壤类型与缺水量的相关分析，可以得出以下两个结论：

①地形指数IRDG曲线的坡度与本文确定的缺水量曲线坡度具有正相关性，IRDG曲线坡度缓时对应的缺水量曲线坡度也缓，流域缺水量分布较均匀。

②当确定流域平均蓄水容量后，全为砂土质的流域缺水量空间分布相对最均匀，其次为壤土质流域，黏土质流域缺水量空间分布差异最大。

Claims (1)

1.一种基于地形和土壤特性的分布式流域缺水量测算方法，包括：

步骤一，流域下垫面基本数据处理：地理数据采集、计算栅格指标和核算栅格土壤信息；

步骤二，核算流域范围内的栅格地形指数：一是根据流域高程栅格，分别计算栅格坡度角、积水面积，获得任意栅格的地形指数；二是确立流域栅格内地形指数，获得全流域地形指数分布；

其特征在于还包括以下步骤：

步骤三，计算流域范围内的土壤特性参数：一是根据流域栅格土壤信息，分别计算每个栅格的土壤类型及其所占比重；二是确立流域栅格内土壤特性参数信息，建立土壤类型与土壤水分指标的映射关系；

步骤四，定量表达流域缺水量的分布：计算栅格型土壤缺水量指标，过程包括：

(一)利用无降雨补给情况下土壤水分稳定状态分布曲线，以VG模型描述，积分获得达到田间持水量状态的土壤缺水量：

$W_f\left(D\right)=(D-{\mathrm{\ψ}}_c)({\mathrm{\θ}}_f-{\mathrm{\θ}}_r)-\[\frac{D({\mathrm{\θ}}_s-{\mathrm{\θ}}_r)}{{\[1+{\left(\mathrm{\α}D\right)}^n\]}^{1/n}}-\frac{{\mathrm{\ψ}}_c({\mathrm{\θ}}_s-{\mathrm{\θ}}_r)}{{\[1+{\left({\mathrm{\α\ψ}}_c\right)}^n\]}^{1/n}}\]---\left(1\right)$

(二)利用空间尺度转化方法将土壤参数栅格化，结合实验方式获得各类土壤栅格内的土壤参数，通过这种方式获得式(1)中的θ_s、θ_r、ψ_c、θ_f、α、n参数；其中：θ_s为土壤饱和含水率；θ_r为土壤凋萎含水率；θ_f为土壤田间持水率；ψ_c为田间持水量时毛细作用水柱上升的高度；α为土壤特征曲线的经验参数1；n为土壤特征曲线的经验参数2；

将式(1)中的地下水埋深D与地形指数通过建立如下关联式(2)～(4)：

$D_i=\stackrel{\‾}{D}-S_{zm}({\mathrm{tp}}_i-\stackrel{\‾}{\mathrm{\λ}})=-S_{zm}{\mathrm{tp}}_i+C---\left(2\right)$

$C=\stackrel{\‾}{D}+S_{zm}\stackrel{\‾}{\mathrm{\λ}}---\left(3\right)$

${\mathrm{tp}}_i=ln\frac{a_i}{{\mathrm{tan\β}}_i}---\left(4\right)$

式(2)中：D为地下水埋深；S_zm为参数；为流域地下水平均埋深；为流域平均地形指数；tp_i为地形指数；a_i为第i点上游流域的面积；β_i为第i点的角度坡度；

(三)将式(4)代入式(1)计算栅格型土壤缺水量指标，该指标由流域内各栅格点缺水量与地形因子及土壤特征的相关表达式：

$\begin{array}{c}{W}_f{\left(D\right)}_i=(-S_{zm}{\mathrm{tp}}_i+C-{\mathrm{\ψ}}_c)({\mathrm{\θ}}_f-{\mathrm{\θ}}_r)\\ -\[\frac{(-S_{zm}{\mathrm{tp}}_i+C)({\mathrm{\θ}}_s-{\mathrm{\θ}}_r)}{{\[1+{\left(\mathrm{\α}\right(-S_{zm}{\mathrm{tp}}_i+C\left)\right)}^n\]}^{1/n}}-\frac{{\mathrm{\ψ}}_c({\mathrm{\θ}}_s-{\mathrm{\θ}}_r)}{{\[1+{\left({\mathrm{\α\ψ}}_c\right)}^n\]}^{1/n}}\]\end{array}---\left(5\right)$

流域中各栅格点地形指数不同，土壤特征各异，通过式(5)得到每个栅格点上相应的缺水量指标；

步骤五，绘制流域缺水量累积分布曲线：根据步骤四的结果，计算各子流域内栅格点的缺水量，按照从小到大的顺序排列，再根据栅格面积比绘制各子流域离散化的地形指数累积分布曲线，即流域的缺水量分布曲线，从而得到分布式流域缺水量测算模型。