基于DEM的流域河网信息自动提取算法

刘永和张万昌（ 河南理工大学资环学院，焦作 454000； 中国科学院遥感与数字地球研究所，北京 00094）

基于DEM的流域河网信息自动提取算法

刘永和1张万昌2
（1 河南理工大学资环学院，焦作 454000；2 中国科学院遥感与数字地球研究所，北京 100094）

分布式陆面水文过程的模拟，除了经典水文模型所必需的流域气象、水文信息外，还需要研究流域的详细地形、水系信息，方可实现流域内产流和汇流的时空演算。以往通常需要借助商业软件在研究区DEM上提取这些水系信息，不仅耗时、操作不便，还使得模型结构松散。在自主研发的一种分布式水文模型的基础上，开发了一套用于提取水系信息的模块，并与该水文模型以同一种语言紧密耦合，程序完全采用面向对象的方式自主编写，大部分数据交换均在内存中直接进行，而无需占用磁盘空间，运行速度快，易于今后不断完善和扩展。详细介绍了模块中填洼与平坦区域处理、流向与累积流向矩阵生成、Strahler河道等级的确定、子流域生成、汇流次序和流程长度信息的生成等多种流域河网信息提取的具体算法及程序实现。本模型系统完全采用自主方式开发，克服了以往使用商业软件提取河网信息时的限制，使得分布式水文模式的流域模拟及分析功能更强大齐全和今后进一步完善和扩展。

分布式水文模型，数字高程模型，水系提取，耦合模块

0　引言

模拟河流和土壤中水流的分布式水文模型被广泛应用于水利工程及环境影响评价[1-5]。分布式模型可用于估计与水流有关属性的空间变化，如化学物质迁移、土壤侵蚀、水分循环的陆面过程以及洪水状况。由于流域水文状况受天气和气候条件影响，人们期望在实现数值天气预报和气候变化预估的同时，能够实现流域的水文预报或未来水文条件预估。由于水流严格地受制于流域的地形特征，只有预先提取出流域的水系信息，才能实现分布式模拟计算地表径流和壤中流的汇流。基于数字高程模型（DEM）提取流域水系是目前一种较为高效且可靠的手段。目前通过现代雷达测量手段已能获得高精度的地面高程信息，如美国地质调查局发布的SRTM DEM的精度已达到90m的水平分辨率，利用该数据能够较为精准地提取水系河网信息。因此，从DEM中由计算机自动提取汇水网络成为水文模拟研究中广泛应用的手段[6]。

通常由DEM提取的水系信息包括像元间的流向、河段以及对应的子流域[7-8]。有关的研究已有很多文献发表[3,7,9-13]。基于矩形规则格网的DEM是流域提取时应用最多的数据[7]，但在一些文献中也会基于不规则三角网（TIN）DEM[14-15]。本文主要指基于规则格网结构的DEM。采用D8算法提取水系已经接近成熟，具体细节见文献[6, 16]。目前在一些成熟的商业软件中已集成了有关基于D8方法的水系提取工具，如ArcGIS中的空间分析工具、ENVI软件的插件River Tools等。然而基于商业软件的水系提取工具其源代码不公开，人们无法知道其算法是否足够高效且健壮，更重要的是商业软件仅能作为一种数据处理工具，用户无法对其处理过程和结果数据进行自由操控，限制了它在自主开发的分布式水文模型中的应用。 分布式水文模型通常表现为一种执行水文计算任务的计算机程序，由于缺少良好的人机交互界面而难以使用。目前一些重要的分布式水文模型（如SWAT）总是被集成于ArcGIS、GRASS等GIS软件工具中，以充分借助GIS软件所具有的DEM空间分析和可视化功能，增强其易用性。然而这种软件集成方法对于分布式水文模型的研究者而言难度较大，影响了分布式水文模型的应用。因此，开发一套专门的水系信息提取工具，以解决当分布式水文模型应用中存在的问题和不足具有一定的实际意义。

基于上述原因，我们针对分布式水文模型自主开发了一套用于提供水系信息的耦合模块。由于所建的分布式水文模型及水系提取模块采用基于同一种语言（C#）的紧密耦合方式，大部分数据交换均在内存中直接进行，而无需占用磁盘空间，且运行速度快；程序完全采用面向对象的方式自主编写，便于建立批量处理执行任务，易于不断完善和扩展。

1　耦合系统的功能与耦合方案

1.1系统的功能

本文所述的流域河网信息自动提取系统及其与分布式水文模型的耦合系统具有以下三类功能：

1）类似于GIS软件中各种位图和简单矢量图显示与漫游等可视化的功能。考虑到水文模型的主要任务是进行流域水文模拟，本系统既具有可显示单幅DEM或其他地理资料矩阵的功能，也可任意打开并显示用ARCGIS里所指定的ASCII二维数据格式或其他用ASCII码矩阵方式存储的二维数据。当光标在地图上移动时，系统可自动计算并显示当前光标处像元的值、行列号、经纬度、UTM坐标值。本系统可以作为查看任意兼容格式数据内容的基本工具，如查看分布式水文模拟的降水、蒸发、径流的空间分布状况。为使显示的DEM底图或其他属性图像与矢量地图相联系，该系统还提供了设置和显示点类型和线类型矢量图的功能。矢量图显示的数据源为ARCGIS软件导出的文本矢量格式。

2）自动提取大量流域河网信息的功能。本系统中的流域河网提取模块具有计算流向矩阵、累积流向矩阵、水系与河道等级矩阵、流域边界、子流域划分、汇流等级矩阵及汇流次序表、流程长度矩阵和汇流时间矩阵等功能，为分布式水文模拟提供必要的信息。

3）分布式水文模型的模拟功能。本系统还具有按逐日和逐时进行径流和洪水模拟的功能，可按马斯京根法实现河道洪水演算以及按等流时线滞时方法进行汇流。

1.2系统的耦合方案

水文模型在运行过程中所需要的气象、土壤、植被、地形等空间分布数据均采用与所用DEM完全一致的分辨率。为方便使用该系统，在运行过程中系统将以程序配置文件中所记录的河流出口断面所在的像元位置（行号和列号）作为运行起始点，一次性完成所有河道信息提取过程。水文模拟过程中所需的所有河道信息以公用数据结构变量的形式保存在内存中，供水文模型中的产汇流模块使用。

2　流域河网信息提取的算法

文中对流域河道信息的提取采用基于D8方法的思想。在每个DEM栅格上，待提取的信息包括流向、累积流向、水系与河道等级、流域边界、子流域信息、汇流次序、汇流时间、流程长度矩阵。在计算流向前，需要对DEM进行填洼（depression filling）和处理平坦区域，以便能够使DEM每个像元通过流向与主河道连接起来。

2.1填洼和平坦区域的处理

DEM中的洼地是指局部位置的像元高程低于外围所有像元，使计算出的流向呈内流现象，即流向链条不能与外围的主河道相连接。DEM上的洼地可能是地表上真实的洼地，如小的湖泊，也可能是由于测量误差以及分辨率较低所造成的假性洼地。在水文研究中，无论流域中真实存在的局部洼地，还是假洼地，都被视作流域中排水网络的组成部分，即必须通过流向将这些区域与水系主河道相连接。填洼曾经是DEM水系提取中最难解决的问题，有关处理该问题的算法较少，用传统的一些算法（如Janson等[11]在1988年提出的方法）对一幅1000行×1000列的DEM进行填洼需花费数小时的时间，且因易出现错误而使运行不稳定，但Planchon等2002年提出的快速算法已解决了该问题[12]，本文使用了这种新算法，一般仅需数秒的时间即能完成，且运行稳定。

平坦区域中，由于所有像元的高程都相等，不存在任何流向，阻碍了水系的生成。目前常用做法是用不同程度的低于分辨率微小数值抬高平坦区域中像元的高程，使抬高的像元与其邻接像元存在高程差，便于计算流向，并通过流向与外侧水系相连接。

2.2水系生成与河道等级的确定

生成水系时需要完成流向计算、流向累积矩阵计算和确定河道等级。

在本系统中，我们采用一种新的邻域位置编码方法（图1）。在这种编码方法中，编码值排列为三行。从左到右，第一行中的编码设为0、1、2，第二行中的编码设为3、4、5，第三行的编码设为6、7、8。需要注意的是，在被设为4的3×3的编码窗口中心的编码不代表任何邻域位置，而是代表了当前指定像元的位置。由这些方向编码很容易算出当前指定像元的邻域像元的行列号。假定当前像元的位置为（i，j），该像元处的水流流向编码为m，则其下游邻接像元的行列号为（i+m/3-1，j+m%3-1），其中符号%表示取余运算。如果一个流向编码为n的邻域像元的流向指向本像元，则其流向编码应为8-n。可见，这种流向编码方法非常有用，且很直接，这是因为当前感兴趣像元的流出和流入像元的位置能够高效计算出来。

本系统中采用单流向算法来计算流向，算法步骤是：（1）首先对DEM边界上像素的流向设为指向DEM外部；（2）对除边界外的所有内部像素，比较其与周围8个邻接像元的高程差构成的梯度，选择高程低于当前像素且梯度最大的邻接像元所在的方向作为流向。这里的梯度是以当前像元与邻接像元的高程差除以像元距离D来计算。对上左下右四个紧邻像元，像元距离D为1，而对角线上的四个像元，像元距离D为。

流向累积矩阵中每个像元上的累积流向值被定义为流过该像元的所有上游像元的总数，它代表了流过指定像元的上游汇水面积。常用的统计当前像元累积流向的方法有两种，一种是递归找出当前像元的上游像元，统计这些找到的像元总数，即为累积流向值。另一种方法是逐个遍历DEM中的每个像元，当遍历到某个像元时，该像元的所有下游像元的累积流向值加1。当DEM中的像元全部遍历一遍时，便得到了累积流向矩阵。

在得到流向矩阵和累积流向矩阵的基础上，可进一步对河流进行基于Strahler分级方法[17]的水系分级，即把位于顶端的不再有分支的河流称为第1级河流，由两个以上的同级河流汇合组成更高一级河流，如果多条河流同时汇聚，但最高级别的支流只有一条，则汇合后的河流级别与其最高级别的支流的级别相同。文献[18]对Strahler河流分级方法进行了较详细的描述，采用从最上游的像素开始向下游搜索，较为可行，但可能有一定的搜索难度。而本文则提出了另一种不需要从最上游像素开始搜索的方法。具体步骤如下：

1）给定一个累积流向（汇流面积）阈值，认为大于该阈值的像元位于河流上，不大于该阈值的不属于河流；

2）在DEM中找出所有累积流向值（汇流面积）大于指定阈值的所有像元，即找出所有位于河网上的像元。为便于快速检索，将这些河网像元放入一个以行号和列号为关键字的字典集合中；

3）对字典集合中的所有河网像元逐个遍历，对找到的每个像元执行第（4）～（8）步操作；

4）若当前像元（用C1表示）已被分配子流域标志，则返回第（3）步，即跳过当前像元，对遍历到的下一个像元进行处理；

5）从当前像元的8个邻接像元中找出属于C1上游的所有像元，若没有找到这样的像元，则将C1的河流级别设为1（属初级河流）并返回第（3）步执行下一趟循环。若找到了这样的像元，执行后续步骤；

6）对C1的上游像元，统计最大河流级别值m，以及拥有最大河流级别值的像元个数n，并判断是否C1的全部上游像元已被标记过河流级别，如果已被标记则执行后续步骤，否则返回第（3）步执行下一趟循环；

7）如果在前面已找出的C1的上游像元个数为1时，则C1的河流级别设为与其上游像元相同的值；如果C1的上游像元多于1个时，且有两个或更多像元级别具有相同的上游最大级别m时，C1的河流级别设为m+1；如果C1的上游像元多于1个，但只有一个上游像元具有最大级别m，则C1的河流级别设为m；

8）返回第（3）步；

9）如果第（2）～（8）步中没有任何像元被标记过，则退出；否则从字典集合中删除已标记过的像素记录，并返回第（2）步继续执行；

10）结束后退出。

2.3子流域划分

子流域一般是指流域面积超过一定面积的分支流域。文献[19]中提出了一种子流域划分方法，但算法不够具体。文献[20]中提出了具体的算法，但略显复杂。本文采用思路与文献[20]基本一样，但采用下面的算法来完成：

1）生成一个空队列Q，将代表河流出口处的像元设为子流域最低编号（如0或1），并将该像元的记录加入队列；

2）若队列不空，从队列中推出一个像元C1，执行后续步骤，否则退出执行；

3）从像元C1的8个邻接方向中逐个查找C1的上游像元；

4）若找到C1的一个上游像元为C2，将C2加入队列；

5）判断C2的累积流向（流域面积）值大于给定阈值，且C2与C1之间的累积流向值之差也大于给定阈值的条件是否成立，若成立则将C2标记为新子流域（其标记值为C1的标记值加1），若不成立，则C2的标记与C1相同；

6）返回第（2）步继续执行。

上面算法第（5）步中的条件十分关键。C2作为C1的上游像元，其上游流域面积必须大于阈值才能成为一个子流域。另外，还需要保证C1存在除C2外的较大分支才能使C2成为新子流域的像元，这是因为，如果除C2外的其他分支的汇流面积过小，小分支无法独立成为子流域，而与C2、C1属同一子流域。步骤（5）中的条件“C2与C1之间的累积流向值之差也大于给定阈值”正是为了保证C1存在除C2外的较大分支。

2.4汇流等级矩阵、汇流次序表、流程长度与汇流时间矩阵的生成

在分布式水文模型的洪水模拟中，需要根据流域中每个像元间的流向关系从上游向下游采用马斯京根方法进行河道汇流。仅仅根据流向矩阵无法得出正确的汇流结果，仍需要根据汇流路径的远近对所有栅格划分汇流顺序，从而得到整个流域或子流域的汇流等级矩阵和汇流次序表。对于同样属最低级别的河流像元，模拟时距离流域出口像元远的像元应该先计算，距离流域出口近的像元应后计算。汇流等级矩阵中的元素值为理论上根据到流域出口的距离来确定的汇流次序，等级值小的像元要先计算汇流，相同汇流等级的多个像元应同时计算汇流。但从实际计算的角度来看，汇流等级相同的多个像元仍是分别来计算汇流的，但计算顺序可是任意的。为方便，一般要导出汇流次序表，用以代替汇流等级矩阵。汇流次序表以汇流的实际计算次序存放流域中所有的像元的行列号，而不是采用矩阵形式存放。

流程长度矩阵是与流域DEM完全对应的二维矩阵（数组），其每个元素存放流域中对应像元通过流向拓扑关系确定的通向流域出口的空间距离。汇流时间矩阵与流程长度矩阵类似，只是其每个元素存放的是对应像元的“水流”流往流域出口所花费的时间。汇流时间矩阵可以通过将流程长度矩阵直接除以平均流速粗略得到，也可以通过将各个像元的水流流速与相邻上下游像元间的梯度建立关系来精确导出。

汇流次序表可以通过汇流等级矩阵或流向矩阵导出。从流向矩阵导出的方法是从河流出口逆向上溯的方法搜索。由于流程长度矩阵的计算以及通过梯度精确计算汇流时间矩阵时也需要采用同样的逆向上溯方法来计算，可将三个计算任务用同一个算法来完成。算法仍借助队列来实现，具体步骤如下：

1）将流域出口像元入队列Q，并往汇流次序表中添加该像元，对应该像元的流程长度矩阵设为1（对汇流时间矩阵可采用类似的方法赋值）；

2）当队列不空时，从中推出一个像元C，否则跳到第（6）步执行；

3）从8个邻接像元中查找属于C的上游流入像元；

4）对找到的所有上游流入像元的流程长度值进行赋值，对斜对角方向的流入像元，流程长度值为C的流程长度加1.414，对上下左右四方向的流入像元，流程长度值为C的流程长度加1。用类似的方法可对汇流时间矩阵进行赋值。把找到的上游流入像元按任意顺序添加入汇流次序表的末尾；

5）返回第（2）步执行；

6）将汇流次序表反转；退出。

3　其他地面信息的生成

除与水系有关的信息外，分布式水文模型仍需要借助地面坡度和坡向信息来计算汇流速度和水分蒸散发量。计算坡度和坡向的方法可参见文献[21]。

为使流域位置上DEM的地形状况更加直观，可生成漫反射光照明暗图。在计算机图形学中，漫反射光照采用Lambert余弦定理来实现[22]。Lambert余弦定理是根据入射光向量与反射面的法向来计算的，因此先要计算出DEM中的每个像元所处位置上的法向。由几何学知道，空间中的平面法向为通过平面内的任意两个向量的向量积。显然，只要已知平面中的三点，就能确定该平面的法向，因此三角形是最适于光照处理的面片。而规则格网DEM中每相邻四个像元都能拼成两个三角形（图2）。需要指出的是，在拼三角形时需要保证三角形的顶点顺序都为逆时针或都为顺时针。设三角形的三个顶点分别为则三角形某两边上的矢量

当每个DEM像元的三角形法向被算出后，再给定一束来自无穷远处的虚拟入射光，即可算出像元上的光强。

4　流域河网信息提取实例

本文截取了山东省境内的沂河流域范围的DEM数据，分辨率为240m，在此基础上生成水文模拟时所用的各种与流域和河网有关的信息。图3显示的是系统的主窗体，其主要工作区中显示了DEM图像以及该区域内的河流、县界及水文测站的矢量信息。图4显示了DEM原图及由其导出的Strahler河流等级、子流域分区、流程长度、汇流次序等信息的图像。图4b显示的由DEM导出的河网与现有的矢量主河道基本一致，这表明流域河网信息提取模块中的填洼及计算累积流向的算法是较为可靠的。图4c显示了以沂河流域内临沂水文站为出口点的面积阈值超过2000像素的子流域划分，由图可见，划分的子流域基本合理，但仍存在不同子流域面积差别较大的问题。图4d显示了由DEM导出的沂河流域内各像素至河流出口处的流程长度，图4e显示了由DEM导出的各像元的汇流次序，这两幅图表明对应的算法计算结果是正确的。图4f显示了由DEM直接生成的虚拟化光照明暗图，该图能够反映研究区域地形起伏的状况。

5　结语

本文介绍了在分布式水文模型中基于DEM提取流域河网信息的耦合模块的算法设计以及耦合系统的功能和耦合方案，对水系生成和河道等级的确定、子流域生成、汇流次序及流程长度信息的生成以及光照明暗图的生成等算法进行了详细介绍。耦合模块采用基于内部数据结构的紧密耦合方式与分布式水文模型耦合成一体化系统，对水文模型所需的各种流域河网信息采用自动批量式生成，节省了对磁盘空间的占用，提高了运行速度，同时也简化了对水文模拟时流域信息的准备过程，能够较大幅度的提高使用效率。耦合模块中的各种算法结果基本合理或正确，运行速度快，十分稳定。

本水文模型系统的开发完全采用自主方式，完全避免了使用商业GIS软件在基于DEM准备各种河网信息时的限制和不足，更有利于今后的进一步完善和扩展。

[1]Yang X J. Use of LIDAR elevation data to construct a highresolution digital terrain model for an estuarine marsh area. International Journal of Remote Sensing, 2005, 26(23): 5163-5166.

[2]Gumbo B, Munyamba N, Sithole G, et al. Coupling of digital elevation model and rainfall-runoff model in storm drainage network design. Physics and Chemistry of the Earth, 2002, 27(11-22): 755-764.

[3]Wang X H, Yin Z Y. A comparison of drainage networks derived from digital elevation models at two scales. Journal of Hydrology, 1998, 210(1-4): 221-241.

[4]Wolock D M, Mccabe G J. Differences in topographic characteristics computed from 100-and 1000-m resolution digital elevation model data. Hydrological Processes, 2000, 14(6): 987-1002.

[5]Dodov B A, Foufoula-Georgiou E. Floodplain morphometry extraction from a high-resolution digital elevation model: A simple algorithm for regional analysis studies. IEEE Geoscience and Remote Sensing Letters, 2006, 3(3): 410-413.

[6]Tribe A. Automated recognition of valley lines and drainage networks from grid digital elevation models - A review and a new method. Journal of Hydrology, 1992, 139(1-4): 263-293.

[7]Turcotte R, Fortin J P, Rousseau A N, et al. Determination of thedrainage structure of a watershed using a digital elevation model and a digital river and lake network.Journal of Hydrology, 2001, 240(3-4): 225-242.

[8]Paz A R, Collischonn W. River reach length and slope estimates for large-scale hydrological models based on a relatively hill highresolution digital elevation model. Journal of Hydrology, 2007, 343(3-4): 127-139.

[9]Colombo R, Vogt R V, Soille P, et al. Deriving river networks and catchments at the European scale from medium resolution digital elevation data. CATENA, 2007, 70(3): 296-305.

[10]Ahamed T, Rao K G, Murthy J. Automatic extraction of tank outlets in a sub-watershed using digital elevation models. Agricultural Water Management, 2002, 57(1): 1-10.

[11]Jenson S K, Domingue J O. Extracting topographic structure from digital elevation data for geographic information-system analysis. Photogrammetric Engineering and Remote Sensing, 1988, 54(11): 1593-1600.

[12]Planchon O, Darboux F. A fast, simple and versatile algorithm to fill the depressions of digital elevation models. CATENA, 2002, 46(2-3): 159-176.

[13]Jana R, Reshmidevi T V, Arun P S, et al. An enhanced technique in construction of the discrete drainage network from low-resolution spatial database. Computers & Geosciences, 2007, 33(6): 717-727.

[14]Tucker G E, Lancaster S T, Gasparini N M, et al. An objectoriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks. Computers & Geosciences, 2001, 27(8): 959-973.

[15]刘学军, 王永君, 任政, 等. 基于不规则三角网的河网提取算法.水利学报, 2008, 39(1): 27-34.

[16]Mcmaster K J. Effects of digital elevation model resolution on derived stream network positions. Water Resources Research, 2002, 38(4): 13-1-13-8.

[17]Strahler A N. Quantitative analysis of watershed geomorphology. Transactions AGU, 1957, 38(6): 913-920.

[18]张渭军, 王文科, 孔金玲, 等. 基于数字高程模型的水系自动生成. 大地测量与地球动力学, 2006, 26(4): 41-44.

[19]沈中原, 李占斌, 李鹏, 等. 基于DEM的流域数字河网提取算法研究. 水资源与水工程学报, 2009, 20(1): 20-28.

[20]叶爱中, 夏军, 王纲胜, 等. 基于数字高程模型的河网提取及子流域生成. 水利学报, 2005, 36(5): 531-537.

[21]汤国安, 刘学军, 闾国年. 数字高程模型及地学分析的原理与方法. 北京: 科学出版社, 2005.

[22]唐泽圣, 周嘉玉, 李新友. 计算机图形学. 北京: 清华大学出版社, 1994.

A DEM Based Drainage Networks Extraction Module Seamlessly Integrated in a Distributed Hydrological Model

Liu Yonghe1, Zhang Wanchang2
(1 Institute of Resources and Environment, Henan Polytechnic University, Jiaozuo 454000 2 Institute of Remote Sensing and Digital Earth (RADI), Chinese Academy of Sciences, Beijing 100094)

For simulation of the hydrology in a watershed by distributed hydrological models, besides the meteorological and hydrological data needed, the detailed information such as topography and watershed networks is also necessary. Traditionally, the drainage networks were derived from digital elevation models (DEM) using commercial software, which is time consuming and inconvenient for operation in hydrological modeling. Based on a distributed hydrological model (DHM) developed by us, an automatic watershed information extraction module was developed which is able to be integrated into the DHM seamlessly using the same computer language of C#. Therefore, most of the data transferring can be finished in memory with no occupation on hard disks, so it is highly efficient when running. The algorithms for the watershed information extraction were introduced, such as the algorithms of removing of depression and flatten areas, generation of flow direction and accumulating flow direction, obtaining Strahler river order, dividing of sub watershed and calculation the sequence of flow order and water-drainage length. This system overcomes the limitations of using commercial software for extracting watershed information, and also favors the modeling system for convenience of modifying and extending the model easily in the future.

distributed hydrological models, digital elevation models, drainage networks extraction, coupled module

10.3969/j.issn.2095-1973.2015.02.009

2013年12月12日；

2014年3月11日

刘永和（1976—），Email: sucksis@163.com

资助信息：国家自然科学基金项目（41105074和41275108）；中科院数字地球重点实验室开放基金项目（2011LDE010）；河南理工大学博士基金（B2011-038）