张晓娇 焦裕飞 刘佳 李旺林 李传哲

摘 要：针对流域水循环模拟中的子流域划分问题，以大清河流域为例，首先基于ArcGIS软件，以30 m分辨率的DEM作为数据源，进行填洼处理、流向判定、汇流累计量计算等;然后采用河网密度法和水系分维法确定最佳集水面积阈值，提取出与实际最符合的河网，以此更准确地提取流域水系;最后根据水流流向、汇流累计以及出水口的位置划分子流域，并基于MATLAB程序对子流域及其内部栅格进行编码。结果表明：当集水面积阈值为270 km2时，提取的大清河流域河网水系与实际较为符合，且对子流域的划分及编码可以实现快速准确定位每一个栅格，方便计算机的识别及演算。因此，以DEM为基础实现大流域尺度下子流域划分是可行的，对子流域的划分及编码可以为开发构建大尺度流域分布式水文模型提供技术支撑。

关键词：DEM;水系提取;子流域划分;大清河

中图分类号：TV213.4文献标志码：A

doi：10.3969/j.issn.1000-1379.2020.06.003

Study on Method of Sub-Basin Partition of Daqing River Based on DEM

ZHANG Xiaojiao1，2，JIAO Yufei2，LIU Jia2，LI Wanglin1，LI Chuanzhe2

（1. School of Water Conservancy and Environment，University of Jinan，Jinan 250022，China;

2. State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin，China

Institute of Water Resources and Hydropower Research，Beijing 100038，China）

Abstract：Aiming at the sub-basin partition in the basin water cycle simulation， taking Daqing River basin as an example， based on ArcGIS software and taking the DEM with the resolution of 30 m as data source for depression filling， determination of flow direction and calculation of cumulative amount of confluence， this paper adopted the river network density method and water system fractal method to determine the optimal catchment area threshold and extract the most practical river network， in order to extract the drainage system of Daqing River basin more accurately. Finally， the molecular watershed was divided according to the flow direction， confluence accumulation and the location of the outlet， and MATLAB program was used for the sub-watershed coding. The results show that when the threshold of catchment area is 270 km2， the extracted river network system of Daqing River basin is more consistent with the actual situation， and the partition and coding of sub-basins can achieve rapid and accurate positioning to each grid， which is convenient for computer identification and calculation. Therefore， it is feasible to realize subcatchment of large watershed scale based on DEM， and the coding of subcatchments can also provide technical support for the development and construction of distributed hydrological model of large watershed scale.

Key words： DEM;ditch extraction;subcatchments partition;Daqing River

子流域劃分是构建分布式水文模型的基础，分布式水文模型以子流域作为计算单元，可以清晰地了解单元之间的水文过程，并可很好地与传统水文模型相结合，缩短模型开发时间[1-2]。Jong-Gun K.等[3]研究认为，子流域划分对水文模拟结果的影响大于模型参数的影响，而且大尺度流域所划分的子流域的数量不同使模拟精度存在不确定性。目前，基于数字高程模型（DEM）的流域地形特征提取方法已经成为获取流域信息的主要手段。OCallaghan等[4]基于水文学原理确定水流的路径，并引入集水面积阈值的概念提取连续河网。传统方法下作为子流域和河网划分依据的集水面积阈值主要靠人为经验设定，划分结果与实际误差较大。为了提高精度，学者们进行了许多研究，目前确定最佳集水面积阈值的方法主要有河道平均坡降法[5]、河网密度法[6]、流域宽度分布法[7]、水系分形法[8]等。 笔者针对大尺度流域的子流域划分方法进行研究，借助ArcGIS与MATLAB软件，对大清河流域进行河网提取、子流域划分及编码，以期为开发构建大清河流域分布式水文模型提供技术支撑。

1 子流域划分方法

借助ArcGIS中的Hydrology工具箱以及MATLAB程序，对流域原始DEM进行处理，步骤如下。

（1）地形预处理。原始的DEM中大多存在洼地或尖峰，会导致提取的水系出现断裂或者水流方向错误的现象[9]，因此需要对原始的DEM进行填洼处理。

（2）提取水流流向。D8算法[4]分别用数字 1、2、4、8、16、32、64、128 来表示东、东南、南、西南、西、西北、北、东北8个方向（见图1），利用该算法生成流域水流流向。

（3）计算汇流累计量。栅格汇流累计量的大小代表形成地表径流的难易程度，数值越大，越易形成地表径流。

（4）提取河网水系。通过设定一个合理的集水面积阈值，根据水流流向及汇流累计量提取河网水系[10]。采用河网密度法及水系分维法两种方法来确定最佳集水面积阈值：①河网密度法，选取不同的集水面积阈值对流域水系进行提取，计算不同阈值下的河网密度及河源密度，分别绘制集水面积阈值与河网密度、河源密度的关系曲线，曲线趋于平稳的点即为最佳集水面积阈值;②水系分维法，采用盒维数法[11]，使用边长为r的正方形网格来覆盖水系，确定覆盖到水系的网格数N，分别取边长及网格数的对数进行线性回归，拟合直线的斜率即为水系的分维值，由此可以通过计算不同集水面积阈值下的水系分维值，与实际水系的分维值进行比较，来确定最佳的集水面积阈值。

（5）河网分级。采用Strahler分级法对流域提取的水系进行分级[12]，分级的原理是将河网中所有没有支流汇入的河流作为第1级河流，两个1级河流汇成第2级河流……依此类推，一直到河流出水口。在Strahler分级中，有且仅有同级别的两条河流汇成一条河流时，河流的级别才会增大。

（6）确定子流域范围。子流域的划分是基于河流流向栅格和河流分级栅格，子流域内汇流量最大的栅格为子流域的出水口，将所有能流向同一个子流域出水口的栅格划在同一个子流域范围内[13]。

（7）子流域编码。通过一定的编码规则确定栅格之间的拓扑关系，以此来反映栅格的计算顺序，便于计算机识别[14-15]。基于栅格之间的水流流向以及汇流关系，对每个栅格赋予一个5层编码：第一层编码A表示子流域的编码;第二层编码B表示该栅格在子流域中的汇流层数，子流域出水口栅格汇流层数为1，流入出水口栅格的栅格层数为2……依此类推;第三层编码C表示汇流层内的栅格顺序，用来区分同一层内的栅格;第四层编码D表示流入栅格的汇流层;第五层编码E用来区分河道与非河道栅格，其中河道栅格自出水口处从1开始编码，非河道栅格编码为0。

2 应用实例

2.1 研究区概况

大清河流域位于海河流域中部，流域面积4.3万km2，地跨京、津、冀、晋4省（市）。大清河发源于太行山，分为白沟河水系和赵王河水系。该流域属于典型的温带季风性气候区，四季分明，降水量年际变化大、年内分配不均，冬季和春季降雨少，夏季暴雨频繁。

2.2 河网提取

借助ArcGIS软件，对大清河流域原始DEM进行填洼处理（见图2），并利用D8算法提取水流流向（见图3），计算汇流累计量。

河网提取的重点是通过设定一个最佳的集水面积阈值，使得提取的河网与实际河网水系相吻合。对于同一个流域，选取不同的集水面积阈值将得到不同的河网，阈值的选取也有很大的主观性和随意性，因此本文采用两种方法确定最佳集水面积阈值。

（1）河网密度法。分别选取9、90、225、270、315、360、405、450 km2作为集水面积阈值，对大清河流域进行水系提取，统计不同阈值下的河流长度、河源数量，计算出相应的河网密度与河源密度，结果见表1。

将集水面积阈值分别与河网密度、河源密度进行趋势线拟合（见图4、图5）。集水面积阈值与河网密度、河源密度都存在幂函数关系，拟合优度R2都大于0.99。对所得到的幂函数求二阶导数并再次进行趋势线拟合，将拐点作为最佳集水面积阈值（见图6、图7）。当集水面积阈值为270 km2时，河网密度与河源密度二阶导数都趋于0且之后不再变化，因此确定大清河流域水系提取的最佳集水面积阈值为270 km2。

（2）水系分形法。分别用边长r为100、500、1 000、2 000 m的网格对不同集水面积阈值下提取的水系進行覆盖，统计覆盖水系的网格数目N，并用边长及网格数的对数之间的线性关系计算分维值，结果见表2。

分别用边长为100、500、1 000、2 000 m的网格在ArcGIS中对实际水系矢量图进行覆盖，得到覆盖水系的网格数，并用边长及网格数的对数之间的线性关系计算分维值，结果见表3。由表3可知，实际水系分维值为1.032，与集水面积阈值为270 km2时提取的水系分维值相同，因此确定最佳集水面积阈值为270 km2。

综合上述两种方法，可确定对大清河流域进行水系提取的最佳集水面积阈值为270 km2。

2.3 河网分级

大清河流域河流分级见图8。当最佳集水面积阈值为270 km2时，大清河流域提取的河流一共有85条，其中：1级河流7条，占总河流数量的7.6%;2级河流10条，占总河流数量的10.7%;3级河流24条，占总河流数量的29.8%;4级河流44条，占总河流数量的51.9%。由此可见，大清河流域水系密布、发育较好。

2.4 子流域划分

基于DEM提取的大清河流域面积为43 209 km2，与实际收集的数据基本吻合，可以认为该DEM精度能够满足提取流域特征的要求。集水面积阈值为270 km2时，可将大清河流域划分为56个子流域（见图9），其中：子流域面积小于500 km2的占51%，在500～750 km2之间的占21%，在750～1 000 km2之间的占17%，大于1 000 km2的占11%。

2.5 子流域编码

对子流域的编码主要通过MATLAB程序来实现。图10为大清河流域某一子流域内所有栅格的编码，以流域中某一栅格（18，3，5，3，150）为例，第一层编码18表示第18个子流域，第二层编码3为该栅格在第18个子流域中的第3层汇流层，第三层编码5表示栅格为第3层汇流层中的第5个栅格，第四层编码3表示流入该栅格的汇流所在层数为3，第五层编码150表示第150条河道编码。

3 结 论

本文采用河网密度法与水系分维法相互验证，确定大清河流域最佳集水面积阈值为270 km2，该条件下提取的水系与实际的水系吻合度很高。在此基础上，应用ArcGIS软件对流域河网水系进行分级并划分子流域，然后借助MATLAB程序对子流域进行编码，可以实现对流域内每一个栅格快速、精准定位，有利于计算机的识别及演算。

通过对大清河流域水文特征进行提取与研究，发现大尺度流域的子流域划分关键在于确定最佳的集水面积阈值。由DEM数据提取的河网与实际误差较小，可以满足水文模拟的需要。通过基于水流流向与汇流计算的子流域划分及编码，可以充分了解分布式水文模型构建过程中栅格之间水量交换等水文过程，为开发适合大尺度流域的分布式水文模型奠定基础。

参考文献：

[1] 贾仰文，王浩，严登华.黑河流域水循环系统的分布式模拟（I）：模型开发与验证[J].水利学报，2006，37（5）：534-542.

[2] 王中根，郑红星，刘昌明，等.黄河典型流域分布式水文模型及应用研究[J].中国科学E辑：技术科学，2004，34（增刊1）：49-59.

[3] JONG-GUN K，YOUNSHIK P，DONGSUN Y，et al. Development of a SWAT Patch for Better Estimation of Sediment Yield in Steep Sloping Watersheds[J]. Journal of the American Water Resources Association，2009，45（4）：963-972.

[4] OCALLAGHAN J F. The Extraction of Drainage Networks from Digital Elevation Data[J].Computer Vision Graphics and Image Processing，1984，28（3）：323-344.

[5] 熊立华，郭生练.基于DEM的数字河网生成方法的探讨[J].长江科学院院报，2003，27（4）：14-17.

[6] 孔凡哲，李莉莉.利用DEM提取河网时集水面积阈值的确定[J].水电能源科学，2005，23（4）：65-67，93.

[7] 杨邦，任立良.集水面积阈值确定方法的比较研究[J].水电能源科学，2009，27（5）：11-14.

[8] 李丽. 分布式水文模型的汇流演算研究[D].南京：河海大学，2007：32-53.

[9] 徐新良，庄大方，贾绍凤，等.GIS环境下基于DEM的中国流域自动提取方法[J].长江流域资源与环境，2004，13（4）：343-348.

[10] 宋向阳，吴发启，赵龙山，等.基于DEM的延河流域水文特征提取与分析[J].干旱地区农业研究，2012，30（4）：200-206.

[11] SKUBALSKA-RAFAJOWICZ E. A New Method of Estimation of the Box-Counting Dimension of Multivariate Objects Using Space-Filling Curves[J]. Nonlinear Analysis： Theory，Methods & Applications，2005，63（5-7）：e1281-e1287.

[12] MELTON M A. A Derivation of Strahlers Channel-Ordering System[J]. Journal of Geology，1959，67（3）：345-346.

[13] 葉爱中，夏军，王纲胜，等.基于数字高程模型的河网提取及子流域生成[J].水利学报，2005，36（5）：531-537.

[14] 刘佳嘉，周祖昊，贾仰文，等.分布式水文模型子流域编码方法对比分析[J].河海大学学报（自然科学版），2017，45（1）：22-29

[15] 王春，江岭，陈泰生，等.基于Pfafstetter规则的流域编码算法并行化方法[J].地球信息科学学报，2015，17（5）：556-561.

【责任编辑 张华兴】