CN107944111B - 基于改进图论和水文模拟的河网水系连通度计算方法 - Google Patents

Abstract

本发明公开了一种基于改进图论和水文模拟的河网水系连通度计算方法，包括：步骤A、提取河网数字水系图并将其概化为拓扑结构，依据图论基本原理建立描述河网结构的邻接矩阵；步骤B、建立研究区HEC‑HMS水文模型，计算典型降雨过程下河网各节点处出现的流量极大值；步骤C、通过构造连通因子将上述两步骤中求得的邻接矩阵与流量值相耦合，生成改进邻接矩阵；步骤D、基于邻接矩阵和改进邻接矩阵计算河网水系连通度。本发明从河网的静态结构和动态径流过程两个角度，重新设计了图论中的连通程度计算方法，使本发明的应用范围较传统方法有明显扩大。

Description

基于改进图论和水文模拟的河网水系连通度计算方法

技术领域

本发明涉及河网水系连通性评价领域，尤其涉及一种基于改进图论和水文模拟的河网水系连通度计算方法。

背景技术

水系连通的实践可以追溯到公元前2400年，如邗沟、灵渠等我国早期的连通工程，在兴利除害、治国安邦中发挥了十分重要的作用。虽然实践由来已久，但“河湖水系连通”的概念直到2009年才首次被明确提出，其理论研究目前尚处于起步阶段，关于水系连通评价方法的研究也很不充分。

长江水利委员会在《维护健康长江，促进人水和谐研究报告》中将“水系连通”定义为河道干支流、湖泊及其他湿地等水系的连通情况，反映水流的连续性和水系的连通状况。其包含两个基本要素：(1)水流的连接通道；(2)能满足一定需求的保持流动的水流。相应的，评价连通程度的好坏也取决于两个条件：(1)连接通道是否保持畅通；(2)水流在满足一定需求情况下的连续性。

目前我国评价水系连通的方法大多局限于水流通道的连接情况或者是水系的结构特征，而忽略了水量问题。例如《基于图论的河道—滩区系统连通性评价方法》(赵进勇，2011)一文，将图论应用于水系连通性评价，但图论仅能研究水流通道是否连通，并未考虑通道的承载能力；又如与本专利申请技术水平最接近的《基于水流阻力与图论的河网连通性评价》(徐光来，2012)一文，以河道水流阻力倒数表征水流畅通度，对图论法进行了改进，但是水流阻力由河道断面和水深推导而来，忽略了降雨、蒸发、渗漏、人类活动等对水量的影响，且该方法仅适用于平原地区水位变化较小的流域，因此用水流阻力来表征水量交换能力不够全面与准确。

发明内容

本发明所要解决的技术问题是针对背景技术中所涉及到的缺陷，提供一种基于改进图论和水文模拟的河网水系连通度计算方法。

本发明为解决上述技术问题采用以下技术方案：

基于改进图论和水文模拟的河网水系连通度计算方法，包括以下具体步骤：

步骤A)，提取研究区河网数字水系图并将其概化为拓扑结构，依据图论基本原理建立描述河网拓扑结构的邻接矩阵；

步骤B)，建立研究区HEC-HMS水文模型，计算典型降雨过程下河网各节点处出现的流量极大值；

步骤C)，通过构造连通因子ω_ij将所述邻接矩阵与河网各节点处出现的流量极大值相耦合，生成改进邻接矩阵；

步骤D)，基于邻接矩阵和改进邻接矩阵计算河网水系连通度。

作为本发明基于改进图论和水文模拟的河网水系连通度计算方法进一步的优化方案，所述步骤A)具体包括：

步骤A.1)，收集研究区DEM高程数据，使用ArcGIS软件中水文分析工具箱，对数据执行填洼、流向、流量、栅格计算器、栅格河网矢量化的窗口命令，提取数字河流网络，并对照Google Earth影像图进行适当修正；

步骤A.2)，将研究区内河流网络按以下规则概化为拓扑结构：

河道概化为线；

水库、湖泊、河道间交点概化为点；

入海口、入江口与河流源头忽略；

步骤A.3)，将拓扑结构图中的点依次编号为V₁，V₂，…，V_n-1，V_n，并定义邻接矩阵A＝(a_ij)_n×n，其中，n为拓扑结构图中点的数目。

作为本发明基于改进图论和水文模拟的河网水系连通度计算方法进一步的优化方案，所述步骤B)具体包括：

步骤B.1)，收集研究区降雨、径流、气温、太阳辐射的水文气象数据；

步骤B.2)，构建流域模块、气象模块、控制运行模块以及时间序列数据模块的HEC-HMS模型组件，并进行地表径流模拟、基流模拟、直接径流模拟和河道汇流模拟；

步骤B.3)，选用五十年一遇的典型降雨过程作为模型输入，输出流域径流过程，提取河网各节点处的流量极大值，用以表征节点的水量传输能力。

作为本发明基于改进图论和水文模拟的河网水系连通度计算方法进一步的优化方案，所述步骤C)如下：

构造连通因子其中，e为自然底数，q_i为顶点V_i处河道峰值流量；

并将邻接矩阵A中的元素a_ij替换为ω_ij得到改进邻接矩阵B。

作为本发明基于改进图论和水文模拟的河网水系连通度计算方法进一步的优化方案，所述步骤D)的具体步骤如下：

步骤D.1)，根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间的连接路径数矩阵S(s_ij)_n×n：

式中，k为大于等于1且小于等于n-1的自然数，s_ij为从顶点V_i到顶点V_j长度为1,2,…,n-1的路径数之和；

步骤D.2)，根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间所有连接路径的连通因子之和的矩阵T(t_ij)_n×n：

式中，t_ij为从顶点V_i到顶点V_j之间长度为1,2,…,n-1的连通因子之和，b_ij为改进邻接矩阵B中第i行第j列元素；

步骤D.3)，根据以下公式计算顶点V_i到顶点V_j之间所有连接路径的平均连通度矩阵：

步骤D.4)，根据以下公式计算研究区整体连通度Z，即河网水系连通度：

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1.将河网结构特性与河道实际水文特性相结合，使评价结果更趋于实际；

2.更加符合水系连通定义、应用范围更加广泛。

附图说明

图1(a)图1(b)图1(c)图1(d)图1(e)图1(f)分别为秦淮河流域数字水系提取过程中的DEM填洼效果图、计算流向效果图、计算流量效果图、阈值15000效果图、Google Earth修正效果图、数字水系图；

图2为秦淮河流域河网结构的邻接矩阵；

图3为秦淮河流域HEC-HMS模型构建成果图；

图4(a)、图4(b)、图4(c)、图4(d)、图4(e)、图4(f)、图4(g)分别为秦淮河流域率定期1996年洪水、率定期1998年洪水、率定期2002年洪水、验证期1991年洪水、验证期1999年洪水、验证期2003年洪水、验证期2004年洪水HEC-HMS模型径流过程模拟值与实测值对比图；

图5为秦淮河流域部分节点流量过程极大值示意图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本发明公开了一种基于改进图论和水文模拟的河网水系连通度计算方法，包括下述步骤：

步骤A、提取河网数字水系图并将其概化为拓扑结构，依据图论基本原理建立描述河网结构的邻接矩阵；

步骤B、建立HEC-HMS水文模型，计算典型降雨过程下河网各节点处出现的流量极大值；

步骤C、通过构造连通因子ω_ij将所述邻接矩阵与河网各节点处出现的流量极大值相耦合，生成改进邻接矩阵；

步骤D、基于邻接矩阵和改进邻接矩阵计算河网水系连通度。

所述步骤A具体包括：

(1)收集研究区DEM高程数据，使用ArcGIS软件中水文分析工具箱，对数据执行填洼、流向、流量、栅格计算器、栅格河网矢量化等窗口命令，提取数字河流网络，并对照Google Earth影像图进行适当修正。

(2)将河流网络按以下规则概化为拓扑结构：河道概化为线；水库、湖泊、河道间交点概化为点；入海口、入江口与河流源头忽略。

概化后的拓扑图可记为一个有序二元组G(V,E)，其中，V中元素称为顶点，E中元素称为边，E联接V中的两个点，e_ij为顶点V_i与V_j之间的连线。一条路径所经过的边的总数(包括边重复出现的次数)称为该路径的长度。

(3)将拓扑结构图中的点依次编号为V₁，V₂，…，V_n-1，V_n，定义邻接矩阵A(a_ij)_n×n，其中：

将A的k次方记为A^k(a_ij)，由图论基本原理可知，从顶点V_i到顶点V_j长度为k的路径数恰好为为矩阵A^k(a_ij)中的第i行、第j列元素的值。在n个顶点构成的图中，如果顶点V_i和顶点V_j连通，则必能在两点间找到一条长度不超过n-1的路径，即不全为0；否则，认为两点不连通。由此可得到图是否连通的判定准则：对于判断矩阵如果矩阵S中的元素全部为非零元素，则图G为连通图，否则图G为非连通图。对于连通图，可以进一步将“使该连通图变得不连通所需要破坏的最少顶点数”定义为连通度，以此评价连通图的连通程度。

由于图论仅能判断河网节点之间是否存在连接通道，没有考虑水量问题，而对于实际河网来说，河道“连”而不“通”的现象在水系中普遍存在，即顶点间有河道相通，但河道的调蓄能力无法满足行洪的需要。因此有必要对图论加以改进，进一步分析通道的水量传输能力。

为此，本发明选用HEC-HMS水文模型进行降雨径流模拟。HEC-HMS在短期、单场次的降雨径流模拟中能保证较高精度的特性恰好满足本发明要求，且从应用范围来看，该模型在平原、丘陵、高原、山地等各种地形中都有较好的模拟效果。因此选用HEC-HMS模型来计算连通性指标对地形的限制较小，克服了传统评价方法仅适用于平原河网的局限性。

所述步骤B具体包括：

(1)收集研究区降雨、径流、气温、太阳辐射等水文气象数据；

(2)构建HEC-HMS模型组件，包括流域模块、气象模块、控制运行模块以及时间序列数据模块。根据研究区的情况，选用合适方法进行地表径流模拟、基流模拟、直接径流模拟、河道汇流模拟。

其中，地表径流模拟常用的方法有初损稳渗法、SCS曲线法等；基流模拟常用的方法有退水曲线法、月恒定流法；直接径流模拟常用方法有运动波法、经验单位线法、snyder单位线法、Clark单位线法等；河道汇流模拟常用方法有马斯京根法、滞后演算法等。

(3)对上述选用方法所需参数进行率定及验证。

以SCS曲线法为例，需要计算或率定的参数为径流曲线数CN，用以反映累计降雨量、土壤质地、土地利用方式以及前期土壤含水量等因素；

以退水曲线法为例，需要计算或率定的参数为衰减指数k；

以Snyder单位线法为例，需要计算或率定的参数有转化系数T_p和峰值系数C_p；

以马斯京根法为例，需要计算或率定的参数有蓄量常数K和流量比重因子x值。

(4)根据《防洪标准》(GB 50201-2014)中规定，一般城市与重要乡村的防洪要求为抵御五十年一遇洪水，假定设计洪水与设计暴雨同频率，则选用五十年一遇的典型降雨过程可最大程度反映流域的排涝能力。因此本发明选用研究区五十年一遇的典型降雨过程作为模型输入，输出流域径流过程，提取河网各节点处的流量极大值，用以表征节点的水量传输能力。

所述步骤C具体包括：

(1)计算评价指标

图论中邻接矩阵所包含的元素“0”和“1”仅能表示顶点间是否有通道连接，但对通道的连接性能缺乏判别能力，因此本发明借助水文模型对水量问题的分析能力来弥补图论这一局限性。

将步骤B中表示水量传输能力的流量极值与步骤A中表示河网结构的邻接矩阵相耦合，生成改进邻接矩阵，耦合方法如下：构造连通因子ω_ij，将邻接矩阵A(a_ij)_n×n中的元素a_ij替换为ω_ij得到改进邻接矩阵，其中：

式中，e为自然底数；q_i为顶点V_i处河道峰值流量。

当q_i＝0时，ω_ij＝0，表示顶点V_i与顶点V_j完全不连通；当q_i→+∞时，ω_ij＝1，表示顶点V_i与顶点V_j完全连通；q_i越大则ω_ij越接近于1。

(2)确定连通度计算方法

设拓扑图G有n个顶点，A(a_ij)为其邻接矩阵，A的k次方记为A^k(a_ij)，由图论基本原理可知，从顶点V_i到顶点V_j长度为k的路径数为矩阵A^k(a_ij)中的第i行、第j列元素的值。则根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间的连接路径数矩阵S(s_ij)_n×n：

式中，k为大于等于1且小于等于n-1的自然数，s_ij为从顶点V_i到顶点V_j长度为1,2,…,n-1的路径数之和。

类比可知，改进邻接矩阵也有相似的性质，矩阵B^k(b_ij)中第i行、第j列元素的值表示从顶点V_i到顶点V_j路径长度为k的连通因子之和。故根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间所有连接路径的连通因子之和的矩阵T(t_ij)_n×n：

式中，t_ij为从顶点V_i到顶点V_j之间长度为1,2,…,n-1的连通因子之和。

根据以下公式计算顶点V_i到顶点V_j之间所有连接路径的平均连通度矩阵：

根据以下公式计算研究区整体连通度Z，即河网水系连通度：

连通度Z即为本计算方法得到的反映河网水系连通程度的数值结果。

下面以秦淮河流域河网水系连通度计算为例，对具体技术的实施步骤进行详细说明，其中关键步骤的实现可参阅相关附图。应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

所述步骤A具体包括：

登录中国科学院计算机网络信息中心下载水平分辨率为30米的2010年秦淮河流域ASTER GDEM数据，导入ArcGIS软件并在ArcToolbox工具箱中依次打开Spatial AnalystTools和Hydrology条目，使用其中的Fill、Flow Direction、Flow Accumulation、RasterCalculator(阈值设为15000)、Stream to Feature等工具处理DEM数据得到水系模拟图，将其转换为kml格式后再导入Google Earth与实际河网进行比对并适当修正，最终得到研究区数字水系图，各项工具处理效果图参阅图1(a)至图1(f)。

将水库、湖泊、河道间交点依次编号为V₁，V₂，V₃，…，V₄₀，根据顶点间是否有河流沟通，构造邻接矩阵A(a_ij)_40×40，参阅图2。

所述步骤B具体包括：

收集秦淮河流域水文、气象数据。本实例中，收集到的水文数据包括1986～2006年21年间流域内的赵村水库、武定门闸等8个站点的日降雨数据和秦淮新河闸、武定门闸两个站点的日径流量数据；气象数据包括对应时段的秦淮河流域逐日降雨、气温数据。

构建HEC-HMS模型组件，包括流域模块、气象模块、控制运行模块以及时间序列数据模块，参阅图3。

流域模块：产流计算方法选择SCS曲线法；直接径流计算方法选择snyder单位线法；基流计算方法选择退水曲线法；河道洪水流量演算选择马斯京根法。

气象模块：降雨方法选择用户自定义法，将各雨量站按泰森多边形分配给就近区域；秦淮河流域地处亚热带湿润季风气候区，蒸散发量较少，故此处不考虑流域蒸散发。

运行控制模块、时间序列模块：将洪水期间各雨量站的雨量时间序列数据及同期观测的流量数据输入可视化数据存储系统HEC-DSSVue中，利用其数学函数功能将不规则时间序列数据进行时间插值，得到每场降雨日累积雨量和观测流量数据序列，供HEC-HMS水文模型中的时间序列数据库调用。

按上述方法建立的HEC-HMS模型主要需要确定的参数有CN值、流域滞时及马斯京根模型中的蓄量常数K、流量比重x。CN值、流域滞时根据流域土地利用情况、河道情况等计算得出，需要率定的参数主要为马斯京根模型中的参数K、x。采用模型自带的单变量梯度搜索法进行参数优化，优化结果如表1所示。

表1秦淮河流域HEC-HMS模型参数

子流域编号	CN	流域滞时(min)	K(h)	x
					1	77	3000	2.3	0.3
2	70	3000	2.5	0.3
					3	70	3000	2.4	0.3
4	70	3000	2.6	0.3
					5	62	3000	2.1	0.3
6	65	3000	2.9	0.3
					7	65	3000	2.8	0.3
8	73	3000	1.7	0.3
					9	82	3000	1.8	0.3
10	78	3000	1.6	0.3
					11	65	3000	2.9	0.3
12	82	3000	3	0.3
					13	80	3000	1.7	0.3
14	80	3000	1.8	0.3

使用秦淮河流域7场典型洪水过程对模型参数进行率定与验证，其中3场用于率定，4场用于验证。

Nash系数用于反映模型的整体效率，Nash系数大于0.8时说明模拟结果较好；相关系数用于反映模拟结果与实测数据的变化趋势是否一致，相关系数大于0.8时说明模拟结果较好；洪量误差及洪峰误差在20％以内认为模拟结果较好。率定与验证统计结果如表2所示，降雨径流过程模拟值与实测值对比图参阅图4(a)至图4(g)。

表2秦淮河流域HEC-HMS模型模拟评价指标

率定期3场洪水的洪量及洪峰相对误差均在20％以内，Nash系数及相关系数均大于0.8，故认为满足要求。验证期4场洪水的洪量及洪峰相对误差基本控制在20％以内，Nash系数及相关系数均大于0.8，故认为满足要求。以上结果表明，HEC-HMS分布式水文模型适用于秦淮河流域的洪水模拟。

通过水文频率计算获取五十年一遇降雨过程并导入HEC-HMS模型，模拟对应径流过程，部分节点流量过程极大值参阅图5。

所述步骤C具体包括：

根据步骤B中得到的流量极大值，计算河网各节点处连通因子ω_ij，计算公式为：

用ω_ij替换邻接矩阵A中的元素得到改进邻接矩阵B。

将矩阵A、B代入下列方程组求得Z＝0.0029。

本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.基于改进图论和水文模拟的河网水系连通度计算方法，其特征在于，包括以下具体步骤：

步骤A)，提取研究区河网数字水系图并将其概化为拓扑结构，依据图论基本原理建立描述河网拓扑结构的邻接矩阵；

步骤A.1)，收集研究区DEM高程数据，使用ArcGIS软件中水文分析工具箱，对数据执行填洼、流向、流量、栅格计算器、栅格河网矢量化的窗口命令，提取数字河流网络，并对照Google Earth影像图进行适当修正；

步骤A.2)，将研究区内河流网络按以下规则概化为拓扑结构：

河道概化为线；

水库、湖泊、河道间交点概化为点；

入海口、入江口与河流源头忽略；

步骤A.3)，将拓扑结构图中的点依次编号为V₁，V₂，…，V_n-1，V_n，并定义邻接矩阵A＝(a_ij)_n×n，其中，n为拓扑结构图中点的数目；

步骤B)，建立研究区HEC-HMS水文模型，计算典型降雨过程下河网各节点处出现的流量极大值；

步骤B.1)，收集研究区降雨、径流、气温、太阳辐射的水文气象数据；

步骤B.2)，构建流域模块、气象模块、控制运行模块以及时间序列数据模块的HEC-HMS模型组件，并进行地表径流模拟、基流模拟、直接径流模拟和河道汇流模拟；

步骤B.3)，选用五十年一遇的典型降雨过程作为模型输入，输出流域径流过程，提取河网各节点处的流量极大值，用以表征节点的水量传输能力；

步骤C)，通过构造连通因子ω_ij将所述邻接矩阵与河网各节点处出现的流量极大值相耦合，生成改进邻接矩阵，具体如下：

构造连通因子其中，e为自然底数，q_i为顶点V_i处河道峰值流量；并将邻接矩阵A中的元素a_ij替换为ω_ij得到改进邻接矩阵B；

步骤D)，基于邻接矩阵和改进邻接矩阵计算河网水系连通度：

步骤D.1)，根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间的连接路径数矩阵S(s_ij)_n×n：

式中，k为大于等于1且小于等于n-1的自然数，s_ij为从顶点V_i到顶点V_j长度为1,2,…,n-1的路径数之和；

步骤D.2)，根据以下公式计算拓扑结构图中顶点V_i到顶点V_j之间所有连接路径的连通因子之和的矩阵T(t_ij)_n×n：

式中，t_ij为从顶点V_i到顶点V_j之间长度为1,2,…,n-1的连通因子之和，b_ij为改进邻接矩阵B中第i行第j列元素；

步骤D.3)，根据以下公式计算顶点V_i到顶点V_j之间所有连接路径的平均连通度矩阵：

步骤D.4)，根据以下公式计算研究区整体连通度Z，即河网水系连通度：