本发明属于水利工程的水循环技术领域，具体涉及一种面向对象的小水库群时空分布式出流计算方法。
背景技术：
在湿润多雨的山丘地区，为了抵御洪水灾害的侵袭，流域内修建大量水利工程特别是蓄水工程，各型水库不仅拦蓄了上游洪水，减轻了所在流域下游的防洪压力，而且为城镇供水提供优质水源，改善了当地水环境，促进经济和社会发展。然而大量修建的小水库数量众多，区域分散，规模主要与所处地形有关，水库之间串并联结构兼而有之，大部分水库建设初期未经水文计算设计，造成各水库拦蓄能力参差不齐，同一地区有的水库频繁泄洪，而部分水库却多年难以蓄满；水库坝体多为土坝，且部分没有闸门控制，很大程度上改变了流域原有的自然面貌，也使得流域的原有水文格局发生了显著变化。目前难以对所有水库调度过程进行实时监控，因而水库群的实际泄水调度过程通常是未知的，因此其对流域洪水过程的影响也是难以确定的，给水利计算及洪水预报带来了较大困难。
虽然目前有关于水库对流域径流影响的研究，但研究对象多针对大江大河上的大中型水库群，由于其管理体系完善，硬件测量技术配套完备，能较精准地掌握其调度运行过程，可以为水文规律的分析及模型建设提供准确的基础资料。但是对于存在众多小水库群的流域来说，由于小水库群实际泄水调度过程无法预知，使得其流域洪水过程的影响也是难以确定的。目前还未见在欠缺实际水库运行资料的前提下，结合地理要素与水文模型，建立时空分布的、水库群差异运行情景下的水库群定量化出流量报道。
技术实现要素：
本发明的目的是为克服现有技术的不足而提供一种面向对象的小水库群时空分布式出流计算方法，该方法可通过建立地理要素和水文模型之间的关系，实现在时空尺度上流域水库群动态设定并参与产汇流计算，从而得到水库群分布式出流量计算模型，能够有效地推动流域水文模型发展。
本发明解决其技术问题是通过以下技术方案实现的：
一种面向对象的小水库群时空分布式出流计算方法，包括以下步骤：
步骤1，流域下垫面基本数据处理：包括采集水库群的地理数据、计算水库群的子流域拓扑信息和收集水库群的基本信息，最终将得到的信息进行数据汇总和整合；
步骤2，核算水库群下泄流量：根据所述步骤1进行数据汇总和整合后得到的信息，计算水库群控制区内的产汇流，再根据水库的库容控制能力，动态调节不同水库的下泄流量；
步骤3，核算时空尺度内水库群状态：包括根据水库建造时间，设立不同时间维度内的水库对象以及根据水库地理位置，建立水库与流域之间的嵌入式管理体系；
步骤4，获得水库群参与的流域水文计算量：计算流域内非水库群控制范围内的产汇流，再将步骤2的结果，参照不同时间维度，带入到步骤3的嵌入式管理体系中，继而进行时空分布式的、有水库群参与的流域径流计算。
作为进一步改进的技术方案，步骤1中所述的采集水库群的地理数据，包括收集流域地形数据、栅格高程数据、栅格土壤特性数据，所述的计算水库群的子流域拓扑信息，包括计算栅格内的坡度、坡向以及河长、划分流域边界、生成子流域拓扑关系与流域水系；所述的水库群的基本信息，包括水库的地理位置、运行时间、控制面积、库容和最大下泄能力。
作为进一步改进的技术方案，步骤2中所述的动态调节不同水库的下泄流量，包括以下步骤：
步骤2.1：计算时段内的水库库容：
V＝V0+V1－V2
式中，V0为时段初库容；V1为计算时段内由于上游来水导致水库库容增加量；V2为计算时段内由于蒸发导致水库库容减少量，
其中，所述V1和V2通过下式求解：
V1＝QΔt
V2＝EAΔt
其中，Q为水库上游流域坡面汇流后得到入库流量，A为水库水面面积，E为单位时段水面蒸发量，Δt为计算时段长；
步骤2.2：根据所述步骤2.1中计算得到的V，计算水库理论出流库容：
其中，Vc为水库溢洪道高程对应的临界库容；
再根据下式计算得到该时段的出流库容：
Vout＝KRES×ΔV
式中，KRES为水库出流系数；
步骤2.3：根据所述步骤2.2计算得到的ΔV，计算水库实际出库流量：
Qout＝min{Qmax，Vout/Δt}
Qmax为每个水库的最大下泄流量，为计算时段长，Qout为Vout/Δt与Qmax之间的最小值。
与现有技术相比，本发明的有益效果如下：
(1)能够依据不同的来水情况和水库自身出流能力，计算小水库群内所有水库的独立出库流量，便于在分布式模型中应用；
(2)二是本发明能够考虑小水库群的建库时间不同，在不同时间维度内设定水库群水量贡献量；
(3)三是本发明考虑小水库群在流域内的空间分布，将水库及其控制区与自然流域通过水量耦合的方式，组建嵌入式管理系统，统一进行径流演算；
(4)本发明所得模型应用于实时洪水预报将大幅度提高有小水库群参与的流域洪水预报的精度，为这类复杂流域的防洪预警等工作提供可靠的科学支撑。
(5)本发明既适用于集总式水文模型，又适用于分布式水文模型，将有效地推动流域水文科学应用研究的深入发展。
附图说明
图1本发明提出的面向对象的水库群时空分布式出流计算方法示意图；
图2格网点水流方向编码示意图；
图3一个简单的DEM及其计算结果示意图；
图4DEM水文分析流程示意图；
图5流域蓄水容量面积分配曲线与降雨径流间关系图；
图6自由水蓄水容量面积分配曲线与各水源关系图；
图7水库控制流域示意图；
图8含水库的流域产汇流模型计算流程；
图9寿昌江流域地形及水库分布图；
图10寿昌江流域水流方向矩阵形式；
图11寿昌江流域累计矩阵形式；
图12寿昌江流域水系；
图13寿昌江流域子流域分布图；
图14流域内水库控制区图；
图15寿昌江流域预报模型建模界面；
图16水文数据管理、导入导出及洪号管理功能界面；
图17站网管理界面；
图18水库管理界面；
图19数据库连接管理界面；
图20模型计算结果显示界面；
图21计算精度统计界面；
图22径流深相对误差；
图23所选场次洪水模拟对比图；
图24 1972年全年洪水模拟；
图25 1972年特大洪水模拟。
具体实施方式
下面通过附图和具体实施例对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。
如图1所示，一种面向对象的小水库群时空分布式出流计算方法，是以GIS(Geographic information system，地理信息系统)生成数字流域水系为基础，以自然流域径流计算和水库下泄能力计算为技术手段，结合空间拓扑关系和时间维度变化，提取各水库出流量的分布情况，最终构建面向对象的小水库群时空分布式出流计算模型。
具体实施时，包括以下步骤：
步骤1：流域下垫面基本数据处理
具体是以DEM(Digital Elevation Model，数值高程模型)数据为基础，计算栅格内坡度、坡向以及河长，划分流域边界，生成子流域拓扑关系与流域水系，其计算遵循以下顺序进行：
(1)洼地处理
被高程较高的区域围绕的洼地是应用数字高程流域水系、生成模型、生成流域水系的一大障碍，因此在决定水流方向以前，必须先将洼地填充。有些洼地是DEM生成过程中带来的数据错误，但另外一些却表示了真实的地形。修改地形数据，对于可能产生的一些不合理的方向阵，处理方法如下：
根据水流特点，通过对洼地区域和平坦区域标示，利用最短流程算法，修改洼地区域和平坦区域高程，使研究区域中的水流能够通过洼地区域和平坦区域。
(2)水流方向矩阵的计算
如图2所示，采用D8(Deterministic Eight-neighbours)算法进行计算，图3(a)中的原始DEM矩阵，每个网格值为格网点高程值，将中心格网点同其最邻近的8个格网点之间的坡降进行比较，其中与落差最大的一个格网点中心之间连线的方向便定义为中心格网点的水流方向，并且规定一个格网点的水流方向用一个特征码表示。有效的水流方向定义为东、东南、南、西南、西、西北、北和东北，并分别用1、2、3、4、5、6、7和8这8个有效特征码表示，表示方法如图2所示。
中心格网点同相邻8个格网点之间单位距离的高差为：
MD＝Z/D
式中：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，且此方向没有指向中心格网，则以此格网的方向值作为中心格网的方向值；
④重复③，直到所有格网都被赋值，得到水流向矩阵如图3(b)所示。
(3)水流累计矩阵的计算
如图3所示，区域水流累计矩阵表示区域地形每一点的流水累积量。其基本思想是，假定以规则格网表示的数字地面高程模型的每一点有一个单位的水量，按照水流从高往低流的规律，根据区域地形的水流方向矩阵，计算每一格网点流过的水量数值，便可得到该区域的水流累计矩阵。下面给出一个从原始DEM矩阵计算出相应水流方向矩阵及水流累计矩阵的范例。
以图3(a)到图3(b)的变换方法为基础，下面说明从图3(b)到图3(c)的变换方法：
图3(b)中，左上单元格标注为“2”的格网点，因为周围没有单元格的水流流入该单元格，该单元格为流域边界，所以在图3(c)中相应填入0，同理可知图3(c)中第一行所填的数均为0；图3(b)中，第二行第二单元格，由于其左上方单元格的水流将流入此单元格，因此，在图3(c)中相应为指填入1；图3(b)中，第二行第四列单元格为例，其左上方和正上方单元格的水流都将流入其中，即有两个单位的水流汇入此单元格，因此，在图3(c)中相应的单元格填入2；再以图3(b)中第三行第四列单元格为例，其左上方和正上方的水流将汇入其中；另外，左上方已经有一个单位的入流，正上方单元格则有两个单位的入流，再加上两个单元格自身的水量，此时，图3(c)中相应单元格所填入的数为5。以此类推，即可由图3(b)得到图3(c)，从而生成水流累积矩阵。
(4)基本水文分析
包括汇水面积计算，流域分水线识别，河网生成，建立流域水系拓扑关系等四部分：
①汇水面积计算
按最大落差规则确定的水流路径可以非常方便地计算出指定格网点以上的流域汇水面积，若以格网数表示面积的话，则该格网点的汇水面积数值是该格网点以上汇入该格网点的格网点数目。在算法上利用递归算法实现，由指定点开始按逆水流方向向上搜索迭代，即可得到该集水流域内任意格网点的汇水面积，其结果如图3(c)表示的水流累积矩阵。
②流域分水线的识别
给定流域的主要入口和出口的断面位置，即它们所在栅格单元的行列坐标。一旦两者的位置确定，根据它们的汇水面积大小，程序可自动搜索从而勾划出流域边界并计算出流域面积。
③河网生成
河网生成分为三个步骤：确定流域边界内的水道；裁减小于某一临界长度的河段；生成河道编码。
首先，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点即定义为水道。
其次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由第一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道、位于河谷两边的凹痕或沟壑的出口，需要将它们裁减除去。
最后，确定河道级数和河段长度。根据流域出口断面确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的编码。同时可以确定各级支流汇入上一级河道的节点，对河网所有节点进行编码，这样定义的节点编码可用于水文汇流计算或河网数据库的构建。
④流域水系拓扑关系
一旦生成联结完好的河网，根据各河网节点，可以确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息。一方面，河网与子流域边界等空间信息是以栅格形式存储，这样易于用GIS软件作可视化显示；另一方面，河段或子流域的拓扑关系以表格形式存储，有利于数字水文模型的调用。详见图4利用DEM进行水文分析的流程示意图。
(5)小水库群基本信息收集
其中包括各水库地理位置、运行时间，控制面积、库容、最大下泄能力等五部分：
①水库地理位置
根据历史资料收集，实地考察等方式，定位流域内每一座小型水库的经纬度坐标，海拔高程，准确定位到流域水系上。将前一步骤生成的流域边界范围进行再分割，得到由水库控制的上游集水区和下游受水区，将原有的流域分割为多个子流域。
②水库运行时间
获取每个水库的建成时间，运行时段。
③水库的控制面积
由水库的地理位置及水库控制的上游集水区范围，生成水库的控制面积。以便核算流域面上降水，入库水量等信息。
④水库的库容
由于流域内的小水库大多没有设计库容，无法获得准确的库容，且属于自由溢洪形式，因此引入“临界库容”概念，即根据资料中显示的水库溢洪道高程，获得对应的库容。
⑤水库最大下泄能力
根据小水库溢流堰的形态，尺寸，计算过堰的最大出流流量。进而分别核算每一个小水库的最大下泄能力。
步骤2：核算小水库群下泄流量
具体是指根据步骤一中得到的下垫面数据，进行水库控制区内的产汇流计算，根据水库的库容、最大下泄能力以及出流能力共同约束，获得水库的实际出库流量，其计算遵循以下顺序进行：
(1)水库控制区产流
根据Horton产流理论以及山坡水文学产流理论，按降雨产生净雨的不同场所，其径流组成主要分为地面径流、壤中流和地下径流等三种。
水库控制区通常处于土壤缺水量不大的湿润地区。在该类地区，一场较大的降雨常易使全流域土壤含水量达到蓄满。倘若一场降雨不能使全流域蓄满，或一场降雨过程中，全流域尚未蓄满之前，流域内也观测到有径流。因此，一个流域的产流过程在空间上是不均匀的，在全流域蓄满前，存在部分地区蓄满而产流。
考虑到流域土壤缺水量空间不均匀，新安江模型采用流域蓄水容量面积分配曲线表征土壤缺水量空间分布的不均匀性。所谓蓄水容量面积分配曲线是：部分产流面积随蓄水容量而变化的累计频率曲线。
流域蓄水容量面积分配曲线线型：
式中f为产流面积(km2)；F为全流域面积(km2)；W’为流域单点的蓄水量(mm)；WMM为流域单点最大蓄水量(mm)；B为蓄水容量面积曲线的指数。
流域蓄水容量面积分配曲线及其与降雨径流相互转换关系见图5流域蓄水容量面积分配曲线与降雨径流间关系图。
(2)水库控制区水源划分
蓄满产流模型计算出的总径流量R中包括了各种径流成分，由于各种水源的汇流规律和汇流速度不相同，相应采用的计算方法也不同。因此，必须进行水源划分。
自由水蓄水库结构考虑了包气带的垂向调蓄作用。按蓄满产流模型计算出的总径流量R，先进入自由水蓄水库调蓄，再划分水源。自由水蓄水库的产流面积上自由水蓄水库设置了两个出口，一个为旁侧出口，形成壤中流RS；另一个为向下出口，形成地下径流RG。根据蓄满产流的概念，只有在产流面积FR上才可能产生径流，而产流面积是变化的，所以，自由水蓄水库的底宽FR也是变化的。还设置了一个壤中流水库，该水库用于壤中流受调蓄作用大的流域，即将划分出来的壤中流再进行一次调蓄计算。
由于饱和坡面流的产流面积是不断变化的，所以在产流面积FR上自由水蓄水容量分布是不均匀的。三水源划分结构是采用类似于流域蓄水容量面积分配曲线的流域自由水蓄水容量面积分配曲线来考虑流域内自由水蓄水容量分布不均匀的问题。所谓流域自由水蓄水容量面积分配曲线是指：部分产流面积随自由水蓄水容量而变化的累计频率曲线。流域自由水蓄水容量面积分配曲线的线型为
式中S’为流域单点自由水蓄水容量(mm)；MS为流域单点最大的自由水蓄水容量(mm)；EX为流域自由水蓄水容量面积分配曲线的方次。
流域自由水蓄水容量面积分配曲线与各水源的关系描述见图6自由水蓄水容量面积分配曲线与各水源关系图。图中，KG为流域自由水蓄水容量对地下径流的出流系数；KI为流域自由水蓄水容量对壤中流的出流系数。
(3)水库控制区汇流
汇流是指流域上地表径流、壤中流和地下径流汇集到流域出口断面的流量过程。
净雨沿坡面从地面和地下汇入河网，然后沿着河网汇集到流域出口断面，这一完整的过程称为流域汇流过程。前者称为坡地汇流，后者称为河网汇流。
为了考虑降雨分布不均和下垫面分布的不均匀性，采用自然流域划分法或泰森多边形法将计算流域划分为N块单元流域，在每块单元流域内至少有一个雨量站；单元流域大小适当，使得每块单元流域上的降雨分布相对比较均匀，并尽可能使单元流域与自然流域的地形、地貌和水系特征相一致，以便于能充分利用小流域的实测水文资料以及对具体问题的分析处理。
对划分好的单元流域运用新安江模型进行产流、汇流计算，得到单元流域出口的流量过程；对单元流域出口的流量过程进行出口以下的河道汇流计算，得到该单元流域在全流域出口的流量过程；将每块单元流域在全流域出口的流量过程线性叠加，即为全流域出口总的流量过程。可见，三水源新安江模型具有三分特性，即流域分单元、蒸散发计算分层、产流量分水源。
(4)出库流量计算
流域中的小水库属于自由溢洪形式，水库上游流域所有的汇流首先进入水库，经过水库调蓄后，超过溢洪道高程的部分水量以不超过下泄能力的流量下泄，下泄的水量再经过流域汇流直到流域出口，完成水库出库流量计算。
以图7水库控制流域示意图为例说明计算方法。水库上游流域按照原有的新安江模型计算，经过流域坡面汇流后得到入库流量Q，在计算时段内则会增加库容V1＝QΔt，而在该时段内由于蒸发导致水库库容减少V2＝EAΔt，A为水库水面面积，则此时段内水库库容净增量为V＝V1－V2，在时段初库容为V0，利用水量平衡原理，获得计算时段内的水库库容V＝V0+V1-V2。水库溢洪道高程对应的临界库容为Vc。当V≤Vc则此时不能泄水；如果V＞Vc时，ΔV＝V-Vc是可以泄流的库容，成为该时段内的自由水。引入水库出流系数KRES，获得该时段的出流库容Vout＝KRES×ΔV，又由于每个水库均有一个最大下泄流量Qmax，因此真实的出流流量应该是Vout/Δt与Qmax之间的最小值。
再将该出流量转化为非水库控制面积上的径流深，加入到三水源分水后流域的地表水径流深中，参与流域坡面的后续汇流，因此包含小水库的流域新安江模型计算框图如图8含水库的流域产汇流模型计算流程。
步骤3：核算时空尺度内小水库群状态
具体是指根据水库建成时间，运行时段，将水库的运行时间与本身建立匹配，并在长序列时间计算中为每一个水库建立一个时间对象，在流域径流计算过程中识别水库的时间对象，动态加载计算时段内相应的运行水库，采用洪号管理的方式，将有时间重叠的水库对象并行计算。再根据水库在流域中的空间位置，以及拓扑关系，建立水库与流域之间的嵌入式管理体系，将水库控制区与非水库控制区流域通过径流交互的方式实现流域内水量一体化管理，其计算遵循以下顺序进行：
(1)水文模型洪号管理
首先选定模型计算的时段，将其根据需要划分为不同的时间洪号，同一洪号内存储相应时段的水文、气象时间序列资料并完成相应的水文计算。由于各水库的建成时间及运行时间都不一致，为每一个水库建立相应的模型对象，为每一个模型对象设定对应的计算时段；再以洪号为依据，动态搜索所属时段内的水库对象，完成识别，参与洪号内计算。
(2)水库对象并行计算
针对流域上同一洪号内并联的水库对象，在水库控制区内进行产汇流计算时，通过并行、互不干扰的方式进行计算，缩减水文模型计算时间，提高运行效率。
(3)构建嵌入式耦合模型
基于水库控制区与非水库控制流域的拓扑关系，将流域边界进行分割，同时将水库控制区上游的流域出口与水库控制区流域水系构建衔接方式；将水库控制区流域出口与其下游流域的水系相衔接，作为集中入流方式，将水库出库流量与下游流域构建耦合衔接方式，完成水库控制区与其他流域的嵌入式搭建。
步骤四：获得水库群参与流域水文计算量
具体是指根据步骤一中得到的下垫面数据，计算非水库控制区的产汇流情况，参考步骤二的结果，将水库群对象按照洪号分组归类，将同一类的水库群对象的径流计算结果与非水库控制区流域的径流计算量通过径流结合在一起，完成流域径流计算，其计算遵循以下顺序进行：
(1)非水库控制区产汇流
参照湿润区蓄满产流模型新安江模型进行流域产汇流计算。
采用流域蓄水容量面积分配曲线表征土壤缺水量空间分布的不均匀性。根据各种水源的汇流规律和汇流速度，采用相应的计算方法。通过水源划分，将产流的总径流量划分为地表、壤中、地下径流。这三种径流汇入流域出口断面的过程即为汇流过程。
(2)洪号控制下的流域嵌入式径流计算
基于步骤三中完成的水文模型洪号管理和嵌入式模型搭建，在设定的计算洪号内，对流域内进行产汇流计算，并将参与洪号计算的水库群对象一并加入到模型计算中。依据流域与水库控制区的拓扑关系，以及上、下游的径流耦合衔接方式，将模型按照空间从上游到下游的计算顺序，即水库上游流域作为上游来水模块参与计算，将水库出库流量作为集中入流，汇入到下游流域内完成全流域的径流演算，实现流域内面相对象的水库群参与的时空分布式出流量计算。
实施例
流域控制站源口水文站位于杭州建德市寿昌江干流，所在流域称“寿昌江流域”，流域控制面积约689km2。选取寿昌江流域为应用实施区域，利用本发明的全流域有小水库群参与的时空分布式计算模型，结合地理数据驱动的流域水系生成、水库控制区划分，嵌入式管理体系以及面向对象的洪号管理方式，再利用寿昌江流域水文气象数据以及小水库群数据等，驱动全流域有小水库群参与的时空分布式计算模型，模拟寿昌江流域水文过程。
步骤1：流域下垫面基本数据处理
(1)采集寿昌江流域DEM资料
寿昌江流域经度118°53′42″～119°19′55″，纬度29°12′33″～29°27′36″。由公开的30米分辨率的ASTGTM数据，采集研究区范围内地形栅格数据。
(2)洼地处理
在处理水流方向之前，根据水流特点，通过对洼地区域和平坦区域标志，利用最短流程算法，修改洼地区域和平坦区域高程，使研究区域中的水流能够通过洼地区域和平坦区域；
(3)水流方向矩阵的计算
采用D8算法进行计算，最终得到寿昌江流域水流方向矩阵形式。
(4)水流累计矩阵的计算
根据寿昌江流域栅格水流方向阵计算结果，计算相应的水流累计矩阵，最终得到寿昌江流域累计矩阵形式如图13所示。
(5)流域水系生成
首先，给定寿昌江流域的主要入口和出口的断面位置(经度119°15’8”，纬度29°26’31”)。根据流域的汇水面积大小，按照流域累计矩阵搜索从而勾划出流域水系；
其次，给定最小河道给养面积阀值，小于该值的汇水面积不可能产生足够的径流而形成水道。流域范围内汇水面积超过该阀值的那些格网点定义为水道。
再次，给定最小河道长度，若一级河道的累积长度小于该长度，则该水道被裁减。由前一步生成的一些河道可能会很短，那些很短的一级水道很可能是伪水道，将它们裁减除去。
最后，确定河道级数：四级，河段长度：14～1307m不等。根据流域出口断面—源口水文站，确定干流河道，流入干流的河道被定为一级支流，流入一级支流的河道定为二级支流，依此类推，确定所有河道的拓扑关系，寿昌江流域水系，如图14所示。
(6)流域水系拓扑关系
生成联结完好的河网水系，根据河网节点，确定相应各支流的流域边界线，从而建立河网节点、河段和子流域的拓扑关系，包括河段坡度、高程值、上游集水面积与侧向集水面积及相互连接的拓扑信息，从而勾划出流域边界并计算出流域面积，寿昌江流域共划分出11个子流域，如图15所示。
(7)水库群基本信息收集
寿昌江发源于杭州市李家镇长林大坑源。在大同镇久山湖村，大同溪与劳村溪汇合成干流，经航头、寿昌、更楼，于罗桐埠汇入新安江。流域内的主要支流有大同溪、劳村溪、小江溪、乌龙溪、南浦溪、童家溪、翠坑溪、山峰溪、石岭溪、石马溪、甘溪；干支河流属典型的山溪性河流，河床比降大、源短流急，其流量受降水量控制十分明显，洪枯变化悬殊，水位暴涨暴落。寿昌江流域共有小(一)型水库14座(包括白岭坑水库、洞山水库、公曹水库、红塘水库、黄家村跃进水库、里诸水库、绿荷塘水库、三楂坞水库、石鼓水库、石郭源水库、石堂水库、武塘水库、牙坑水库、杨家口水库)，小(二)型水库48座，流域及主要水库位置分布见图13，水库相关信息表见表1，其总控制面积约占流域总面积超过23％，足以对流域防洪过程产生重要影响。
表1
步骤2：核算小水库群下泄流量
寿昌江流域现有雨量站多个，其中包括5个雨量资料站和多个雨量遥测站，选择流域内具有长系列资料的大坑源、大同、曲斗、寿昌、源口等5个雨量站的逐日降雨和时段降雨资料；源口蒸发站逐日实测水面蒸发资料；流域出口断面源口站的逐日流量和洪水要素摘录资料。考虑到根据流域选用的5个雨量站分布，采用泰森多边形法将研究流域划分成5个单元。
对于流域蒸发信息，流域上蒸发采用的是1965～2012年，源口站实测多年平均水面蒸发量为847.8mm，其中7～9月占全年的42.4％，日最大蒸发量为11.6mm(1977年7月16日)。
利用流域上的降水，蒸发信息，以及下垫面划分的水库控制区面积、出流系数等地理信息数据，根据流域面积上的水库控制区见图16，结合水库径流计算方法进行计算。
水库上游流域按照原有的新安江模型计算，经过流域坡面汇流后得到入库流量Q，在计算时段内则会增加库容V1＝QΔt，而在该时段内由于蒸发导致水库库容减少V2＝EAΔt，则此时段内水库库容净增量为V＝V1－V2，在时段初库容为V0，利用水量平衡原理，获得计算时段内的水库库容V＝V0+V1－V2。水库溢洪道高程对应的临界库容为Vc。当V≤Vc，则此时不能泄水；如果V＞Vc时，ΔV＝V-Vc是可以泄流的库容，成为该时段内的自由水，获得该时段的出流库容Vout＝KRES×ΔV，又由于每个水库均有一个最大下泄流量Qmax，因此真实的出流流量应该是Vout/Δt与Qmax之间的最小值。再将该出流量转化为非水库控制面积上的径流深，加入到三水源分水后流域的地表水径流深中，参与流域坡面的后续汇流，因此包含小水库的流域新安江模型计算框图如图8。
步骤3：核算时空尺度内水库群状态
以ArcGIS平台为依托，研发的可视化的流域径流模拟系统，界面如图15。为了能将水库出流计算与流域径流模拟在模型计算和结果显示方面嵌入式管理，本次还专门设计了水文数据管理，洪号管理，站网管理、水库管理，数据库连接以及水文数据输入、输出模块等，详见图16－21；由于本流域中还包含众多小水库，在模型中考虑了小水库的影响，需要输入小水库相关的数据，提供界面如图18；针对模型的率定和验证，提供的图形化的界面，以及精度统计界面，分别见图21－22。
步骤4：获得水库群参与流域水文计算量
(1)日模率定与验证
选用寿昌江流域1975至2015年共41年期间的资料对模型进行率定和验证，其中以41年的日资料对模型进行计算，1975～1988年作为日模参数率定期，1990～2015年为日模验证期。率定和验证成果见表2，径流深相对误差分布见图23。
在率定期(1975～1988)，多年平均径流深误差-4.29％，确定性系数均值为0.82；验证期(1989～2015)，多年平均径流深误差3.19％，确定性系数均值为0.86。历年的径流深相对误差分布见图22，只有五年的误差绝对值超过了10％，其他年份均在-10％～10％以内，参数满足精度要求。
表2
(2)次洪率定与验证
资料长度为41年，洪水场次众多，共超过120场，从中选取43场洪峰流量超过500m3/s的洪水过程进行次洪计算。其中1975～1994年中的19场洪水用来率定参数，其后24场用来验证模型，得到的精度表格见表3，率定期的平均径流深相对误差为4.01％，洪峰流量相对误差为-7.45％，确定性系数均值为0.93；验证期平均径流深相对误差为8.58％，洪峰流量相对误差为-2.29％，确定性系数均值为0.93。目前次洪模拟精度确定性系数超过0.9，总体达到甲级精度。所有选取的洪水过程模拟和实测对比见图23，图形排列按照场次顺序左右排列。
表3
(3)“1972”特大洪水
1972年8月初，7号台风在浙江省平阳附近登陆，于黄山、安庆间形成低压槽。8月3日因受7号台风倒槽影响，寿昌江流域发生大暴雨，最大日雨量高达269.4mm，2-3日总降水量：白沙389.7mm，大坑源364.7mm，上邹畈376.9mm，大洲446.7mm，寿昌江出口处的源口水文站洪水调查最大流量为3160m3/s。采用目前率定获得的参数对1972年的特大洪水进行验证，1972年全年洪水及特大洪水的对比图形见图24、图25。对于洪水过程时段(1972-8-2 6:00:00～1972-8-8 5:00:00)内的模拟对比统计值见表4。
表4
1972年模拟的结果与调查洪峰和洪量都有较大差距，可能的原因包括：大洪水的前期降土壤含水量不准确，导致部分雨量补足土壤缺水；前期河道流量也不准确，导致马斯京根汇流也不准确；洪水调查数据也可能偏大，导致结果误差过大。
验证结果表明本发明提供的方法可以准确、有效的完成流域内大量水库群共同参与的径流模拟演算，系统界面设计简洁，操作方便，计算结果直观显示，可以得出以下两个结论：
小水库群算法准确方面：根据本发明之方法得到的有水库群参与的流域水文模拟能够实现对41场历史洪水的准确的模拟及预报；对历史上的大洪水事件也能给出合理的计算结论，可见本发明的计算方法可行并有效。
系统设计方面：采用面向对象的方式，能够动态、灵活的配置参与模型计算的水库；采用洪号管理以及嵌入式系统设计可以直观，准确的实现时间、空间分布式的模型计算，可见本发明的系统设计方法灵活有效。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。