考虑湖泊范围的子流域划分和编码方法研究

关铜垒，刘佳嘉，周祖昊，井 涌，王鹏翔，杜 崇

(1.黑龙江大学 水利电力学院，黑龙江 哈尔滨 150080；2.中国水利水电科学研究院流域水循环模拟与调控重点实验室，北京 100038；3.陕西省水文水资源勘测局，陕西 西安 710068）

1 引言

分布式水文模型计算单元划分方法主要有网格单元、地貌单元和子流域单元［1］3 种,其中:基于坡面漫流模型的子流域单元划分方法应用较为广泛且被成功集成到各类GIS 软件及水文分析软件(如ArcGIS、ARCSWAT、QSWAT 等)中［2-3］。一些学者在前人研究基础上进行改进,以满足不同情形的子流域划分要求,如刘欢等［2］基于传统子流域划分方法,提出一种面向大尺度区域的改进方法,可实现流域出水口的自动识别,并通过多阈值模拟河网融合技术,解决河网提取失真问题所导致的子流域拓扑关系错乱问题；雷晓辉等［4-5］提出一种子流域划分过程中界河、海岸线处理方法,并基于传统子流域划分方法提出一种改进的子流域加密算法,可解决狭长子流域和水库、水文站控制的子流域划分问题；刘佳嘉等［6］提出一种能够描述城市管网汇流作用的子流域划分方法,可有效解决城市管网对子流域划分的影响。然而无论是传统子流域划分方法还是改进后的方法,在进行子流域划分时都未考虑流域内湖泊的相对独立性,一般通过虚拟河网将湖泊划分为若干部分,每个部分分属不同的子流域。湖泊作为流域的重要组成部分,其水量平衡过程对流域水循环过程的影响不可忽视,现有划分方法使得湖泊水量平衡计算变得比较困难,因此有必要将湖泊划分在同一个子流域内进行研究分析,但在将湖泊进行独立划分的时候,会遇到湖泊上游河流汇入和坡面区域编码等问题。鉴于此,笔者提出一种考虑湖泊范围的子流域划分和编码方法,对湖泊整体及其上游汇入区域进行准确划分及编码,并利用流域分布式水文模型(WEP-L 模型)［7-8］进行模拟验证。

2 研究区概况

选取黄河源区(唐乃亥以上区域)作为研究区。研究区总面积12.2 万km2,包含黄河源头两大淡水湖泊扎陵湖和鄂陵湖。这两个湖泊均位于玛多县,其中:扎陵湖东西长35 km、南北宽21.6 km,湖面面积526 km2,湖泊东北部较深、西部较浅,平均水深约9 m,蓄水量46.7 亿m3；鄂陵湖位于扎陵湖之东,南北长32.3 km、东西宽31.6 km,湖面面积610 km2,平均水深约18 m,蓄水量107 亿m3［9］。

本研究所用DEM(见图1)来源于中国科学院信息中心提供的ASTER GDEMV2 高分辨率数字高程数据,空间分辨率为30 m×30 m。为提高计算效率,采用ArcGIS 软件将该数据重采样为1 km×1 km 分辨率数据。流域水系数据来源于国家基础地理信息中心的全国1 ∶250 000 地形数据库。

图1 研究区DEM

3 子流域划分原理及步骤

3.1 设置模拟河网提取及流域栅格属性

模拟河网提取过程采用传统方法进行,即经填洼、流向计算、汇流累计数计算、河道阈值设置、模拟河网提取等操作得到［10-11］。首先,根据河网和湖泊对流域栅格属性进行设置,主要包括A和B两部分,其中:A＝1 代表湖泊、A＝0 代表非湖泊；B＝1 代表河道、B＝0 代表非河道(坡面)。因此,流域栅格属性有湖泊河道(AB＝11)、湖泊非河道(AB＝10)、非湖泊河道(AB＝01)以及非湖泊非河道(AB＝00)4 种。

3.2 子流域划分和编码

根据湖泊范围及汇流关系,将流域分为湖泊区域、湖滨带区域、其他区域3 类区域,其中:湖泊区域为湖泊范围内部区域,湖滨带区域为直接汇入湖泊的无河网坡面集水区域,其他区域则是除湖泊和湖滨带以外的区域。湖泊区域和湖滨带区域采用湖区子流域划分方法进行划分,其他区域采用非湖区子流域划分方法进行划分。两种方法的选择主要由流域栅格属性确定,即在对河段从下游往上游追溯划分的时候,如果该河段起始河道栅格为非湖泊河道栅格(AB＝01),则采用非湖区子流域划分方法进行划分；如果该河段起始河道栅格为湖泊河道栅格(AB＝11),则采用湖区子流域划分方法进行划分。

3.2.1 非湖区子流域划分方法

非湖区子流域划分方法即传统子流域划分方法,具体步骤［5,11-13］:①从流域出口栅格开始,基于栅格流向沿河道从下游向上游追溯,分河段进行编码；②当完成一个子流域划分后,沿河道向上对该河段的上游河道继续追溯,根据流入河道栅格的汇流量确定干支流,将汇流量最大的河道作为干流,将其余的河道作为支流；③先对汇流量小的支流河道进行子流域划分,再对汇流量大的干流河道进行子流域划分,直至河流源头或者遇到湖泊河道栅格结束。

3.2.2 湖区子流域划分方法

湖区子流域划分方法是本研究的创新点,主要思路是将每个湖泊划分为一个独立的湖泊单元,并且根据河网与湖泊进出口节点将汇入该湖泊的湖滨带区域划分为若干个单元。具体步骤如下。

(1)湖区区域预处理。以湖泊范围内的模拟河网与湖泊出入口节点为分割点,对湖泊进行分割处理,将两条模拟河网之间的湖泊范围内栅格作为一个整体,继而将湖泊范围分割为若干个子区域,每一个区域内的栅格相连通且由模拟河网分割开。对分割开的子区域赋予不同的分区编号X,该编号是从1 开始的连续自然数(如图2 的1～4 号),流域内的每个湖泊分割后的子区域编号均从1 开始。

图2 湖泊预处理结果示意

(2)湖泊单元划分。当湖泊下游子流域划分编码完成后,溯源到当前湖泊出口栅格时开始湖区子流域划分。从该栅格出口开始溯源编码,直至遇到非湖泊栅格结束(AB＝00 或AB＝01),将其中所有湖泊栅格(AB＝10 或AB＝11)都赋予相同的编号,如图3 的6 号和12 号湖泊型子流域。

图3 子流域划分结果示意

(3)湖滨带单元划分。根据湖泊预处理结果,以湖泊范围内拥有不同分区编号的非河网栅格为起点,对湖滨带栅格进行溯源编号,将湖滨带流入湖泊内不同分区的湖滨带栅格赋予相应的子流域编号N＋x(N为湖滨带栅格流入的湖泊单元编号,x 为湖滨带栅格最终流入的湖泊栅格所拥有的湖泊子区域编号,如图3 中7～9 号子流域分别对应6＋1、6＋2、6＋3)。重复上述步骤,直至该湖泊上游所有湖滨带栅格都完成子流域划分,至此完成当前的湖区子流域划分,然后根据湖泊上游汇入河道栅格累计数按从小到大的顺序依次执行非湖区子流域划分。

3.2.3 确定子流域类别属性

所有子流域划分完成后,根据子流域内的湖泊范围和河道情况对子流域进行属性划分,主要包括一般类型、湖泊类型和坡面类型3 种。其中:一般类型的子流域包括一条河道栅格以及连接该河道的坡面栅格,无湖泊栅格(如图3 的1～2 号子流域)；湖泊类型的子流域包含该湖泊范围内的所有河道和坡面栅格(如图3 的6 号、12 号子流域)；坡面类型的子流域为连接湖泊区域的坡面栅格(如图3 的7～8 号、13～14 号子流域),不含河道和湖泊栅格。

3.3 子流域上下游拓扑关系确定方法

3.3.1 确定上游流入子流域

本研究中由于存在不包括河网的坡面类型的子流域,无法直接根据河网汇流次序确定子流域间上下游拓扑关系,因此需要根据栅格流向进行溯源遍历(沿模拟河网对流域内所有栅格逐个访问)确定,具体规则及步骤如下。

(1)如果遍历的子流域为一般类型子流域,则根据子流域内的河道栅格进行溯源遍历,查找上游流入的河道所在子流域作为其上游(如图3 一般类型1 号子流域的上游2 号和5 号子流域、一般类型5 号子流域的上游6 号子流域)。

(2)如果遍历的子流域为湖泊类型子流域,则需要先对上游汇入的湖滨带子流域进行遍历,再对汇入的河道所在子流域进行遍历。首先,根据栅格流向,查找湖滨带直接流入湖泊的坡面类型子流域,根据湖泊预处理结果,将湖滨带直接流入湖泊不同子区域的多个湖滨带子流域都作为湖泊类型子流域的上游坡面类型子流域(如图3 湖泊类型6 号子流域的上游7～9 号子流域、12 号子流域的上游13～16 号子流域)；其次,根据河道栅格流向,沿湖泊内虚拟河道向上继续追溯,将汇入的非湖泊河道栅格所在的子流域作为当前湖泊类型子流域的上游(如图3 湖泊类型6 号子流域的上游11号子流域、12 号子流域的上游17～19 号子流域)。

(3)如果遍历的子流域为坡面类型子流域,则直接结束遍历,无上游子流域,如图3 中坡面类型7 号子流域。

3.3.2 建立上下游拓扑关系属性表

在确定子流域上游流入子流域后,需要构建子流域上下游拓扑关系属性表,上下游拓扑关系属性表内容主要包括当前子流域编号、下游子流域编号、上游子流域数目、多个上游子流域编号。

4 实例应用

按照上述方法对黄河源区进行子流域划分。首先采用ArcGIS 软件经过栅格填洼、计算栅格流向、计算栅格汇流量、确定河网阈值等步骤提取模拟河网,然后基于提取的模拟河网设置流域栅格属性,接着根据鄂陵湖和扎陵湖边界及模拟河网对两个湖泊进行预处理,最后对研究区域进行子流域划分以及构建上下游拓扑关系属性表。黄河源区共划分1 871 个子流域,将两个湖泊划分为单独的子流域,其中:鄂陵湖上游有河道汇入的子流域4 个、无河道汇入的湖滨带坡面型子流域5 个,扎陵湖有河道汇入的子流域5 个、无河道汇入的湖滨带坡面型子流域7 个。计算单元划分结果如图4 所示,两个湖区子流域的上下游拓扑关系属性表见表1。

表1 湖区子流域上下游拓扑关系属性

图4 黄河源区考虑湖泊范围的子流域划分结果

4.1 对比分析

为对比本文提出的考虑湖泊范围的子流域划分方法与传统子流域划分方法之间的优缺点,采用ArcGIS软件对黄河源区(唐乃亥以上区域)进行子流域划分,结果见图5(左下角为湖泊范围子流域划分放大图)。对比图4 和图5 可以发现,两种方法非湖泊范围划分的子流域大体相当,但传统方法将湖泊范围划分为若干个子流域,湖泊范围中心处出现密集的异常子流域；而本文提出的方法将湖泊范围单独划分为一个子流域,保证了湖泊的独立性。

图5 黄河源区常规方法子流域划分结果

4.2 水文模拟验证

采用流域分布式水文模型WEP-L 进行模拟验证。为符合本文提出的子流域划分单元的特点,对该模型的汇流模块加以改进,即在汇流模块中识别湖滨带计算单元,因其内部不包含河网,故只对其进行坡面汇流演算,而不进行河道汇流计算。选取玛曲和唐乃亥水文站1956—2000 年月均流量进行水文模拟验证,其中1956—1980 年为率定期、1981—2000 年为验证期,并将Nash 效率系数和相对误差作为水文模拟结果的检验指标,结果如图6 所示。玛曲水文站1956—1980 年率定期的Nash 效率系数为0.764、相对误差为1.6%,1981—2000 年验证期的Nash 效率系数为0.750、相对误差为-2.3%；唐乃亥水文站1956—1980 年率定期的Nash 效率系数为0.806、相对误差为2.1%,1981—2000 年验证期的Nash 效率系数为0.768、相对误差为-1.4%。月径流量模拟结果较好,表明本文提出的子流域划分方法符合分布式水文模型建模的需求,能有效模拟水文站月径流量过程。本文所提出的方法优势在于能通过水文模拟方便地平衡湖泊的水量,为进一步研究湖泊水循环及伴生生态演变提供支持。在改进计算单元划分的前提下,借助WEP-L 模型输出鄂陵湖和扎陵湖的蒸散发量(见图7)和出、入湖水量(见图8)。扎陵湖的年均蒸散发量为597.7 mm,鄂陵湖的年均蒸散发量为746.9 mm；鄂陵湖多年平均出、入湖水量分别为16.1 亿、17.2 亿m3,扎陵湖多年平均出、入湖水量分别为8.3 亿、9.8 亿m3。

图6 玛曲和唐乃亥水文站月径流量模拟结果

图7 鄂陵湖和扎陵湖多年平均蒸散发量

图8 鄂陵湖和扎陵湖多年平均出入湖水量

5 结论

(1)本文提出一种考虑湖泊范围的子流域划分和编码方法,该方法通过对湖泊范围进行预处理,可自动区分湖泊类型子流域上游汇入的河道和坡面区域,将上游汇入的坡面湖滨带划分为若干个子流域,实现对研究区内湖泊区域的自动识别,并将湖泊划分为单独的子流域然后编码。

(2)与现有子流域划分方法相比,本文提出的改进方法可将湖泊划分为单独的子流域,保证湖泊的完整性,便于运用分布式水文模型计算湖泊的进、出水量。

(3)采用流域分布式水文模型WEP-L 进行月均流量模拟验证,结果表明所提出的子流域划分方法能够满足分布式水文模型建模的需求,保证子流域划分过程中湖泊的独立性、完整性和有效性。