基于图论的城市河网水系连通方案优选——以清潩河许昌段为例

窦 明，石亚欣，于 璐，李桂秋，贾瑞鹏

（1.郑州大学 水利科学与工程学院，河南 郑州 450001；2.郑州大学 生态与环境学院，河南 郑州 450001；3.郑东新区管理委员会水务局，河南 郑州 450000）

1 研究背景

河网水系的连通状况对区域水资源配置、水旱灾害防御、水体自净能力提高、水生生境维护以及物质能量传递等都有着重要的影响［1-3］。然而，随着城市化进程的加快，人类活动特别是水利工程建设影响了天然水系的结构形态，改变了河网水系连通性［4-5］。黄草［6］、Pracheil［7］、Cui［8］等指出，通过重新设计河网水系结构、恢复水系连通度，能有效维持河流健康发展。然而，常规的水系连通措施多是借助各种工程来改变水系结构，需要大量资金和人力投入，措施不当还会引起对生态环境的负面影响，因此选择适宜的水系连通方案尤为关键。

近年来，一些学者在水系连通理论方面取得了显著性进展，如水系连通概念界定［9-12］、水系连通效果评价［13-16］等，但整体来看目前研究成果缺乏针对水系连通程度与河网水量变化过程内在联系的定量描述，故无法有效比较不同水系连通方案的实施效果。此外，图论法也被广泛应用于水系连通表征方面［4，13-14，17］，但前期研究多简单将图论用于水系连通度的计算，并没有将河网水量分配过程与连通效果衔接起来。

为此，本文依据图论法原理建立了城市水系图模型，嵌入支流、闸门、排污口等边界条件，定量描述河网连接点与河段的水量平衡关系，在此基础上构建了多闸联合调度的水系连通方案优选模型，优选得到在不同工况下的最佳水系连通方案，力求为城市水系连通格局优化提供典型案例。

2 研究方法

2.1 研究思路本文研究思路为：首先，借鉴图论法原理将研究区水系网格化，以描述河网连接点与河段的连通关系；其次，通过节点河段关系分析河网水量分配方式，考虑上游来水和下游排泄，反映整个河网的入流-自然分流-出流过程；再次，考虑污水处理厂水源汇入和闸门调节，建立多闸联合调度下的水量分配关系，并计算各河段的流量数值；最后，以水生态景观面积最大为优化目标，结合各种流量约束、水力条件约束和河网水量分配关系等，建立多闸联合调度的水系连通方案优选模型，以不同的多闸联合调度方案为方案集，求解研究区最佳的多闸联合调度的水系连通方案。

2.2 河网图模型构建河网图模型就是利用图论中的图模型概念，将水系连通状况通过点和边的连接关系表示出来［13，15］。与传统图模型不同，本文除了考虑河网自然结构连通，还考虑了污水处理厂水源汇入和闸门调节作用，并且对节点进一步分类以突出节点功能。

首先，河网概化。将河段交点和河网边界用点集V 表示，河段用边集E 表示，边的方向表示水流方向，得到河网图模型G =（V，E），如图1 所示。其中，V={vi}，i=1，…，N，N 为节点数；E={ei}，i=1，…，M，M 为边数，每个ei对应一个＜u，v＞，u，v∈V；（uv）表示方向为u 到v，连接两节点的边。图模型可以用矩阵表示，邻接矩阵A={aij}，aij表示vi到vj的边数，aij=1，表示水体可以从vi流向vj，否则aij=0；加权邻接矩阵W={wij}，wij表示vi到vj的流量所占vi总流量的比例，wij的值由vi下游各河段宽的比例确定，设vi到vj的河段宽为bij，vi下游所有河段河宽总和为b，则wij=bij/b，当vi到vj不存在河段时，wij=0。因此，A 可以反映河网水流方向，W 可以反映河网水量分配关系。

图1 河网图模型概化示意图

其次，节点分类：

（1）河网水源节点。①自然输入节点：位于河网上游边界，水源为河网上游外部河流来水，对应v1和v5。②污水输入节点：位于河网内污水处理厂排水口处，水源为排水口来水，对应v3。

（2）水量分配节点。①分流节点：下游有多条河段的河网交点，向下游不同河段分水，根据节点处是否存在闸门可分为人工调节分流节点和自然分流节点，前者对应v2，水量分配方案由闸门调度方案确定，具体过程在下节说明；后者对应v4，水量分配方案由下游各河段河宽的比例确定。 ②汇流节点：上游有多条河段的河网交点，汇集上游多条河段来水，对应v6和v4。③过渡节点：上下游都仅连接一条河段的河网交点，上下游流量一致，对应v3。

（3）水量排泄节点。输出节点：位于河网下游边界，承接河网内部水量并将其排出河网，对应v7和v8。

需要说明的是由于某些节点可以承担多个功能，因此可以同时为多种节点。如，v3同时为污水输入节点和过渡节点；v4同时为汇流节点和分流节点，另外本文的河网除输出节点通过边界向下游排泄，不涉及其他取水过程。

2.3 河网水量自然分配关系根据图模型可以确定河段节点连接关系，基于此，在不考虑闸门调节和其他节点汇入的情况下分析河网水量自然分配关系。此时，v2为自然分流节点，v3为过渡节点；水体经边界自然输入节点进入河网，被河网内部自然分流节点、汇流节点和过渡节点重新分配，再由边界输出节点排出河网。

设河网中流入任一节点vi处流量值为Qi，其上游任一节河段＜vs，vi＞流入该节点的流量为Qsi，其流入下游任一河段＜vi，vj＞上通过的流量为Qij。则有：

对于内部节点：式（1）表示分流过程，若vi为自然分流节点，wij表示分配比例，vi的流量按比例分至vi下游多个河段；若vi为过渡节点，wij=1，节点处上下游流量相等，该式仍成立。式（2）表示汇流过程，若vi为汇流节点，则对于该节点，有多个s 使asi=1，vi处流量等于上游各河段的流量和；若vi为过渡节点，则仅有一个s 使asi=1，上下游流量相等，该式仍成立。对于边界节点：上游边界自然输入节点存在分流过程，可用式（1）表示；下游边界输出节点存在汇流过程，可用式（2）表示。因此，式（1）和式（2）可以表示整个河网的节点流量分配。

在建立河网图模型时，已根据河网结构和河宽确定图模型的邻接矩阵A 和加权邻接矩阵W，将河网的节点流量和河段流量也用矩阵形式表示，则可通过矩阵运算表示整个河网的分流汇流过程。

式中：Qin、QT和Q 为N 阶对角矩阵，对角线元素分别表示自然输入节点的流量值，上游河段汇集在节点处的流量值和所有节点处的流量值三种性质的有序数据系列，数据位置与节点编号一一对应。其中，Qin中表示内部节点和边界输出节点位置和QT中表示自然输入节点位置的数据值为0。Qedge={Qij}N×N为边流量矩阵。式（3）表示河网分流过程，式（4）表示河网汇流过程。由于Qin对角线元素只能表示自然输入节点的流量（除自然输入节点对应位置外所有数据都为0），QT中对角线元素表示除自然输入节点以外其他节点的流量，因此两者相加可以表示所有节点流量（式（5））。假设初始状态下河网内部无水，已知Qin便可根据式（3）—（5）计算河网其他节点和边的流量。

2.4 人工调节下的河网水量分配关系在河网水量自然分配关系的基础上，考虑污水处理厂排水输入和内部闸门的调节作用，建立人工调节下的河网水量分配关系。

（1）污水输入节点汇入：在2.3 小节的河网水量分配关系中加入污水处理厂排水，这时v3同时作为污水输入节点和过渡节点。此时河网水源除天然河流外还有污水水源，整个河网总水量增加。污水流量矩阵SQ为N 阶对角阵，节点处的流量值由SQ相应的对角线元素表示。SQ加上2.3 小节中的节点流量矩阵为此时的节点流量矩阵Q，如下式所示：

（2）人工调节的加权邻接矩阵：加权邻接矩阵反映了河网水量分配关系，考虑闸门调节作用后，v2由自然分流节点变为人工调节分流节点，水量关系发生变化，产生新的加权邻接矩阵。定义一个新变量D 表示闸门过水能力，多闸联合调度的水系连通方案可用S={Di}（i=1，…，N，N 为节点数）表示，但Di只针对有闸门节点取值。对于河网系统中的任一节点vi，当节点处有闸门时，Di表示：

式中：qedge为闸门所在河段的自然径流量；q′edge为闸门所在河段经闸门调节后闸门下游的径流量。（下文的计算针对闸门过水能力展开，闸门开度和过水能力的关系现阶段不做分析）；D 的取值范围为［0，1］，为简化计算，对D 按固定步长0.2 取一组值， D∈{0、0.2、0.4、0.6、0.8、1}。

根据水量平衡，闸门节点处被闸门阻挡的水量进入节点下游其他河段，比例为各河段河宽比，由于研究区河网各节点下游至多两个河段，因此W 矩阵每行至多两个元素。当节点vi处有闸门时，加权邻接矩阵W={wij}第i 行两个元素发生变化，产生新的加权邻接矩阵Wnew={w′ij}，两个矩阵元素之间存在以下关系：

节点vi处有闸门时，Wnew中第i 行的两个元素w′ij表示节点vi到vj的流量所占vi总流量的比例，表示节点vi到除节点vj外的另一节点的流量所占vi总流量的比例（闸门位于河段＜vi，vj＞上），以图1（b）中节点v2为例，闸门位于河段＜v2，v6＞进水端，节点v2所连接的另一河段为＜v2，v3＞，w′ij指w′26，表示节点v2到v6的流量所占v2总流量的比例，节点v2到除节点v3外的另一节点的流量所占v2总流量的比例。无闸门时，第i 行元素与之前相同，根据式（8）可得到Wnew。

（3）人工调节的水量分配关系：在2.3 小节河网系统自然分配关系的基础上，加上污水处理厂的水源输入和闸门调节作用后，河网的节点流量矩阵和加权邻接矩阵发生变化，用新的矩阵替代之前的矩阵，建立人工调节下的河网水量分配关系如下式：

给定边界水源输入条件，由式（9）可对不同闸门调节方案下的河网流量进行计算。

2.5 基于河道流量约束下的水系连通方案优选模型

（1）确定优化指标。为表征不同方案下的水系连通效果，定义水生态景观面积（WA），含义为整个河网所有河段的水面面积之和。WA与水深正相关，水深又与流量正相关，因此该指标可以反映河网流量分配情况；其次WA影响水生生物生存空间大小，以及其他如景观、航运等功能，可以反映河网的功能特性。指标越大，连通效果越好。WA的计算式如下：

式中：Arij为河段＜vi，vj＞的水生态景观面积；WA为整个河网的水生态景观面积，是所有河段水生态景观面积的加和，如果河段＜vi，vj＞不存在，则Arij为零。本文研究区内河网经过修整，河道断面可以概化为梯形断面，任一河段＜vi，vj＞的河道断面示意图如图2 所示。

Bij、Mij、和Hij分别是河段＜vi，vj＞平均断面的底宽，边坡系数和水深。设该河段长为Lij，则Arij计算公式如式（11）所 示，矩阵Ar={Arij}，B={Bij}，M={Mij}，H={Hij}，则Ar 计算公式如式（12）所示。

图2 河道断面示意图

式中：*表示两个矩阵的哈达玛积。

（2）建立优选模型。以最大WA为优化目标，考虑流量约束、水力学约束和2.4 节中建立河网水量分配关系等约束条件，建立优化模型如下：

目标函数为：

约束条件：

①河网水量分配关系约束，见式（9）。

②断面流量约束。

③水力学参数约束。

④WA计算公式，见式（11）（12）。

式（9）为河网水量分配关系，用以计算不同闸门调度的水系连通方案下的河网流量。式（14）表示流量约束，目的是为了筛选河网流量在适宜范围内的连通方案。Qij是vi到vj的流量，流向为vi→vj。式（14）-①中Q限1是为保证河段不断流且筛去流量极小的方案而设置的最小流量，Q限2为防止洪水淹没两岸的最大流量；式（14）-②中Q限3i表示第i 个河段的最小环境流量控制目标，针对需要设置环境流量目标的关键控制河段；式（14）-③中Q限4i表示第i 个河段适宜水生生物生存的最大控制流量，针对有水生生物投放的生态修复河段。3 个流量约束之间存在交叉。式（15）表示河段＜vi，vj＞的水力学约束，断面参数如图2 所示，Aij，xij和Rij表示过水断面面积，湿周和水力半径，nH和iH为河道糙率系数和水力坡度，Cij表示谢才系数。需要说明的是，Q限1和Q限2是针对整个河网取的极值，相当于是一个极小值和极大值，在多数情况下都达不到该值，设定该范围的主要是针对一些不太重要的河段，而对于一些敏感河段或重要河段，则给了更具体环境流量、生物生长流量作为约束如Q限3i和Q限4i。

图3 模型计算流程图

优化模型的计算通过MATLAB 程序实现，计算逻辑为：输入河段参数和输入节点流量，利用式（9）计算不同方案的河网流量，利用公式（14）判断河网流量是否满足流量限制，若满足，利用公式（15）计算河段水深，之后利用该结果计算WA，比较各方案的WA值并进行方案优选。方案通过不同的闸门调度方式确定，具体方案在模型应用时进行说明。计算流程图如图3 所示。

需要说明的是，该方法是在河道断面可以概化为梯形的基础上进行计算的，当河道不满足该条件如为复式断面时需要对涉及断面参数的公式进行调整。此外，由于目前研究区河段都经过修整，渗漏量不大，且该区段主要是景观河道，没有大的供水任务，因此，本次研究暂没有考虑地下水和地表水互通和供水问题，后期为了增强模型的适用性，可以加上渗漏量和供水量两个参数，在计算河段流量Qedge时可以减去渗漏量，并考虑加上供水节点类型，在计算节点流量Q 时减去供水量。同时，由于本文研究区降水量较小，河道内的径流量主要依靠北汝河水源以及污水处理厂水源补充，因此本文并未考虑降水对流量的影响，在扩大模型应用范围时，可以考虑将降水的影响纳入模型。

3 应用研究

3.1 研究区概况及河网概化清潩河属淮河流域沙颍河水系，许昌市境内河流总长51.46 km，流域面积1585 km2。清潩河许昌段天然河道和人工河道纵横交错，结构复杂，闸门众多，河流碎片化严重，同时天然径流不足，依靠北汝河水源和污水处理厂排水补充，为了合理分配水量，满足清潩河许昌段各河段的需求，需要调整闸门调度方案，优化水系连通效果。因此本文以清潩河许昌段为例，对其多闸联合调度的城市河网水系连通方案进行优选，以期为清潩河流域水系健康发展提供有利建议。

研究区内闸门众多，本文选取对水量分配有重要影响的8 个分水闸分析。闸门及污水处理厂排水口分布如下图4（a），对河网概化，对节点和河段进行编号得到图模型如图4（b）。

其次对节点进行分类，具体分类如下：

①自然输入节点：对应图4 中的v1、v6、v16、v27和v36。

②污水输入节点：对应图4 中的v18、v19、v22、v14、v28和v37。

③输出节点：对应图4 中的v38和v39。

④人工调节分流节点：对应图4 的节点v2、v3、v4、v5、v6、v31和v20。

⑤自然分流节点、汇流节点和过渡节点同样按照2.2 小节中的方法确定，不再详细说明。

图4 清潩河许昌段城区河网水系图模型

另外需要说明的是，根据研究区的实际情况，河网中还存在环境流量控制河段e1，e23，e21，e22，e39和e33；水生生物生长河段e23，e25，e11，e20和e21。这些河段存在流量限制，具体数据会在下文进行说明。

3.2 方案优选以八月份为典型月，利用第2 节中建立的基于河道流量约束下的水系连通方案优选模型优选多闸联合调度的水系连通方案，分析优化效果，其具体计算流程如图3 所示。多闸联合调度方案集由D 值确定，D∈{0、0.2、0.4、0.6、0.8、1}，D 有6 个取值，研究区内共有8 个闸，对不同闸门根据D 的可能取值进行排列组合，则共有68种多闸联合调度的水系连通方案。根据《清潩河（许昌段）流域河湖水系2017—2018年度水资源优化调配方案》计算自然输入节点流量数据；根据排水口2017年实测数据计算污水输入节点流量数据；断面数据（河宽、河长、边坡系数）为实测数据。

式（14）中各流量限制设置如下：结合流量数据及实际情况设置Q限1为0.02 m3/s；用水生生物生存的流速最大限制值乘河段断面面积得到各个河段的最大限制流量，由于流域内流量较小，一般达不到该限制，取所有河段最大限制流量中的最小值为河网最大流量限制，设置Q限2为19 m3/s；环境流量Q限3i采用Tennant 法计算；水生生物生存最大控制流量Q限4i计算方法如下：取保障水生生物健康生长的最大水深和流速限制值为计算参数，计算限制流量。河网流量输入，Q限3i和Q限4i计算结果如表1 所示。

表1 模型部分输入参数 单位：（m3/s）

由于研究区内的河流较小，之前没有做过专门的水力学计算，参数获取较为困难，且研究区内皆为规整的城区衬砌河道，水力学参数相对差异不太显著，而本文又重点在模型优化，不在水力计算，因此概化后水力坡度和糙率系数取统一值，式（15）中nH和iH分别取0.033 和0.0008。结合表1 参数代入优选模型，最优的多闸联合调度的水系连通方案和河网流量如图5 所示。

图5 多闸联合调度的水系连通优选方案结果

其中为了书写方便，对图中闸门进行编号，图4（a）中破张闸、孙家门闸、水口闸、长店闸、黄龙池闸、饮马河分水闸、护城河分水闸和天宝河进水闸依次编号为闸1-闸8，下文闸门也按此编号。

为了判断优选结果合理性，在分析实际计算结果之前先根据河网结构对不同闸门D 的取值大小趋势进行定性分析：闸1 从颍汝干渠分水到清泥河，下游水系复杂，并在还未接受其他支流汇入前便要在v7处进行再一次分水，因此需要较大分水量；闸7 从清泥河-清潩河连通渠分水至护城河，而护城河没有其他水源，因此闸7 要承担护城河所需全部水量；闸8 从清潩河分水至饮马河，经许扶运河汇入小洪河，小洪河下游有较大环境目标流量限制，但上游来水不足，因此需要较大分水量；闸4 下游的水汇入清泥河，但上游清泥河有多个支流汇入，因此不需要下放大量水。所以可以判断闸1、闸7 和闸8 的D 值应较大，闸4 的D 值应较小。其他闸门上下游水系结构复杂，不容易确定D 值趋势，但应在0.2～1.0 之间。

根据计算结果，最优方案WA为3.813 km3，闸1、闸7 和闸8 的D 值最大为1，闸4 的D 值最小为0.2，其他闸门D 值在0.2～1.0 之间，结果与定性分析结果相符，结果合理。可以得出结论：通过河网水系结构可以对特殊位置的闸门D 值进行定性分析，一般位置则无法判断，但定性分析不能确定具体数值。而本文建立的优选模型可以对复杂河网的闸门联合调度方式进行优选，确定D 的具体数值，并且优选结果符合水系结构特性，可以作为实际工作的参考依据。

3.3 优化效果对比分析为了进一步分析优化效果，设置对照方案1—5，将每个对照方案闸门D 值设为统一值，方案1—5 闸门的D 值分别为0.2、0.4、0.6、0.8 和1.0，方案5 闸门全开相当于自然状态。计算5 个方案的河网流量和WA，与最优方案结果进行对比，如图6 和图7 所示。

图6 对照方案与优化方案流量对比

图7 对照方案优化目标与限制满足情况

图6 表示方案1—5 与最优方案的河网流量对比，图7 表示方案1—5 与最优方案的优化目标和限制满足情况的对比。由对比结果可知，最优方案大部分河段的流量大于对照方案相应河段的流量，河段流量比方案1 相应河段流量大的有58%，比方案2 有60%，比方案3 有68%，比方案4 有78%，比方案5 有80%。D 值增大，比例增大，这可能是由于D 越大，闸门阻挡作用越小，人工调节越弱，对于本文的复杂河网，人工调节作用减小后，小部分河段集中了相对较大水量，其他河段水量减小，水量分配不合理。另外，最优方案所有河段满足流量限制，而对照方案1—5 都有河段超过流量限制，方案5 超限河段比例最大为16%，这也是水量分配不合理造成的。随着D 值变化，对照方案的WA距最优方案差值的变化没有明显趋势，但方案5 的WA距最优方案差值最大，水系连通效果最差。可以得出结论：在水量不足的复杂河网中，人工调节对河网水量合理分配有明显影响，但一般调节方案难以满足整个河网的水量需求，保障较好的连通效果，本文建立的优选模型可以优选多闸联合调度的水系连通方案，合理分配水量，满足河网水量基本需求，提升水系连通效果，相比自然状态，提升最明显，WA值提高了28 425 m2。

3.4年内方案优选利用优选模型对其他月份的水系连通方案进行优选。每月的自然输入节点流量和控制河段的环境目标流量如下表2 所示，污水输入节点流量和水生生物生长流量限制与前文一致。

表2 其他份自然输入节点参数和限制河段的环境目标流量 单位：（m3/s）

图8 每月最优的多闸联合调度的水系连通方案及优化目标值

将表2 的参数代入优选模型，计算每月最优的多闸联合调度的水系连通方案，结合8月份的优选结果，绘制图8，其中闸1—闸8 含义与上文相同。

从优化目标来看，1—6月WA值较小，7月和8月较大，9月又开始减小，其变化规律可能与来水的季节性丰枯有关。对比7月份和8月份的结果，7月份的流量输入参数虽然比8月份大，但是优化目标值却比8月份小，这是因为式（14）-③筛除了一部分河网内部少数河段流量过大的方案。从闸门调度方案来看，闸1、闸7、闸8 在所有月份的最优方案中，D=1.0，而其他几个闸门在不同月份的优化方案中，D 值浮动较大。因此，在实际调度中，闸1、闸7、闸8 应保持常开状态，其他闸门需要根据计算结果调整D 值。

相比清潩河许昌段目前的调度方案，本方案相对更加灵活。清潩河许昌段目前的闸门调度活动较少，由于整体水量不足，为了保证两侧的生态景观，一般情况下闭闸蓄水，通过一侧泄水口下泄保证下游生态流量，整个河网水流不畅，本文的闸门调度方案在保证河道用水需求的同时增加河道水体的流动性，希望本文的研究工作可以为其提供一些参考。另外，本文的计算结果是针对特定年份进行的，还未考虑水文年的影响产生的不同边界条件，不同的丰、平、枯下的来水条件所得到的结果可能存在偏差，这需要更加详细的分析。

另外需要说明的是，本次的计算是先确定每月的输入条件和限制条件，通过模型对每月的闸门调节方案进行单独优化得到每个月的优化结果。另外一种思路是在开始确定好所有的输入条件和限制条件，通过一次运算将所有月份一齐优化，得到所有结果，不过这样运算时间会呈指数增长，极大地影响运算效率，因此本文还未对此进行深入分析，不过可以作为后续工作深入的切入点。

4 结语

水系连通方案设计是改善水体连通状况、消减水资源不合理开发影响的有效措施。本文利用图论法构建了河网图模型，在对水系连通关系描述的基础上，从水量分配的角度考虑了污水处理厂排水和闸门调度对水量分配的影响，相比传统图模型能更好地适用于城市复杂水系的水系连通效果展示。引用图论法的邻接矩阵和加权邻接矩阵，将河网水量分配关系与河网水系连通效果联系起来，简化了河网流量计算方法。将水生态景观面积最大作为评价水系连通效果的目标函数，并考虑河网水量分配关系、河网流量限制以及水力条件，建立多闸联合调度的水系连通方案优选模型，寻求能够满足河网水量需求并且达到最优连通效果的水系连通方案。通过在清潩河许昌段的应用效果显示，本次优选的水系连通方案可为相关部门开展水资源管理和水系连通工程建设提供借鉴。

相比水文模型对水文过程的概化，本文建立的模型侧重于对河流网络结构的概化以及连通性质的分析，适用于结构复杂但形状规整的河网水系，便于利用图论的原理通过编程操作对河网图模型进行寻优计算；同时，由于计算过程中涉及大量的节点、河网以及数十万种调度方案，计算量巨大，通过本文所建模型的简化处理后在寻优时的计算效率大大提升。不过本文的模型还不够完善，例如地表水地下水的连通，用水户取水以及不同水文年的影响等问题，这些都是我们下一步工作的重点。