基于数字高程模型的混合流向算法

夏誉玲，李小娟，王 涛

1. 首都师范大学北京市成像技术高精尖创新中心，北京 100094； 2. 首都师范大学地球空间信息科学与技术国际化示范学院，北京 100048； 3. 首都师范大学资源环境与旅游学院，北京100048； 4. 首都师范大学三维信息获取与应用教育部重点实验室，北京 100048； 5. 首都师范大学水资源安全北京实验室，北京 100048

从数字高程模型(digital elevation model，DEM)中自动地提取汇水网络和流域范围可以为分布式水文模型及相关应用提供重要的基础信息。在现有提取方法中，基于模拟地表径流累积过程的算法得到了广泛应用，该类算法中，地面各点径流方向和流量分配对提取结果有着直接影响。目前的流向确定算法大致可以分为单流向算法、多流向算法和混合流向算法[1]。

单流向算法[2]是基于模拟地表径流累积过程提取汇水信息的主要算法之一，其基本思想是给定地面点上的积水将流向较矮的最陡坡向相邻地面点，在格网DEM中，给定栅格单元周围有8个相邻栅格单元，该算法又称为D8算法。为了模拟地表流水受地物影响，文献[3]在地形计算中引入随机量，提出了随机八方向算法。D8算法在8个离散的方向选择流向，随着径流的延伸，得到的水流路径与真实径流的偏差会越来越大，为克服该局限性，基于全局搜索策略的D8算法[4]和改进后的GD8算法[5]在D8的基础上扩展了检索范围，并通过定义第2流向以及满足使用第2流向的一系列条件修正了逐渐偏离的水流路径。文献[6]用栅格DEM构建规则三角网，根据三角面的坡向决定流向。单流向算法适于模拟水流的汇流机制，在沟谷地区能够提取出效果良好的河网，然而因其流向单一，无法表现自然环境下的坡面漫流，在平缓区域倾向于产生伪平行流径，在后续水文分析中引入了不确定性因素。

多流向算法适于模拟平缓地形中地表漫流过程。文献[7—8]对D8和Rho8进行了改进，分别提出了水流分配权重为1的FD8和权重值为1.1的FRho8，建立漫散流动模型，使得水可以流向中心栅格周围的所有低高程方向栅格。考虑到不同地形条件下水流分配的不同，文献[9]按FRho8方法计算汇流面积，确定汇流面积阈值，再将水流分配权重值设定为随汇流面积变化的函数，即随着累积汇流面积的增加该算法模拟的水流情况越来越类似于单流向算法。文献[10]提出了基于局域形态确定流向的算法，该方法分析八邻域高程值所反映的地表凸凹几何形态，区分邻域中所有可能被分配水流的独立单元，依据该单元涉及的像素数目决定其分配到的流量。文献[11]认为地形变化对水流分配有着重要的影响，设计了水流分配策略随着下坡坡度变化的多流向算法。

在地表坡面径流过程中，坡度较大的区域径流通常呈汇聚形式，通常表现为单方向流动，因此单流向算法能够较好地模拟这些区域的水流运动；然而，在坡度平缓地表上模拟径流过程时，该算法的结果中会出现不合理的平行流现象[2,12]。在相对平坦的坡面区域水流往往会分流到多个方向，这种情况下多流向算法相比于单流向算法能准确地反映出水流向周围较低高程方向分配的思想。于是出现了一种混合的流向算法，该方法根据地形地貌特征决定适用的流向算法，这样综合了多种算法的优点，避免了无差别的应用单一算法带来的部分问题，能够在不同的地形条件下模拟得到更加合理的水流运动以及水量分配。文献[13]提出了在同一个DEM同时采用不同的单流向算法的方法，该算法用3×3窗口依次对DEM进行曲面拟合，根据该曲面计算得到的正切曲率和流线曲率给中心栅格指定流向算法。文献[14]将研究区分为正地形和负地形，对正地形采用多流向算法，负地形采用D8算法。本文设计实现了一种混合的流向算法，该方法首先对地形进行分类，然后根据地形特征决定适用的流向算法。

最初的地形地貌提取是依靠人工判读勾画，随着数据的丰富和数字地形分析技术和方法的发展，出现了许多地貌形态自动分类方法。现有的地形地貌分类方法有：

(1) 基于地表局部几何形态的方法，该方法的原理是分析局部地表相对起伏度等属性来划分地貌类型。如文献[15]使用了3×3大小的模板扫描DEM，将高程值大于南北或东西向的中心栅格标记为山谷。文献[16]提出了对高程值的下四分位数栅格区域采用自定义直径的圆形模板进行扫描检测山谷线的方法。文献[17]提出的Geomorphons方法，利用地形开放度选取地形特征点来决定与中心栅格高程的相对关系，该算法的统计单元搜索半径不同，具有一定的局域地形条件自适应能力。

(2) 对坡度、曲率等地形因子集合进行聚类的方法[18-20]，即对DEM所派生的地形属性集合进行聚类，之后将聚出的类识别为地形元素类型。

(3) 基于地表整体结构的方法。如文献[21]根据地形表面水流方向提出了基于形态学的沟沿线提取方法。文献[22]通过对等高线进行凹凸段划分，获得等高线特征段，对特征段构建三角网，从而获取特征点，匹配特征段和特征点，完成地形特征线的追踪，生成地形特征线。在山谷或山脊这类陡峭区域，径流通常指向一个方向，由于高程、坡度、剖面曲率等因素，坡面地区容易产生不同程度的漫流。为了简单直接地对单流向和多流向算法进行分配，只需要提取山脊线和山谷线，与坡面地区区分开来。本文采用了基于地表局部几何形态的分类方法，利用3×3窗口内高程的相对差异信息分别提取山谷线和山脊线[15]，对研究区的地形进行分类。

1 混合流向算法

本文提出的混合流向算法实现流程如图1所示。基于地表汇水过程模拟的算法实现过程中，假定地表每一点均有下游点，因此需要对DEM进行洼地填平、平地赋流向等预处理[24-25]。

图1 混合流向算法工作流程Fig.1 Hybrid flow direction algorithm

1.1 地形分类

本文采用了基于局部地表几何形态检测的方法进行地形分类。在3×3模板窗口中，如果南北或者东西相邻栅格的高程值大于中心栅格，则该中心栅格被标记为山谷点[15]；如果南北或者东西相邻栅格的高程值小于中心栅格，则该中心栅格被标记为山脊点；如果同时满足以上两个条件，则标记为鞍部点。当模板完成扫描时，这些被标记的栅格可以连接为山脊线和山谷线。由于消除了平坦区域，经过预处理的DEM中每个栅格单元都有一个坡度值，在模板窗口检测后的DEM中，除了山谷线，山脊线以外的区域，其余区域定义为坡面地区，地形坡度较大区域通常表现为径流汇聚特性，因此利用坡度阈值将坡面区域分为陡坡区和缓坡区。经过计算得到的用于混合流向算法中流向分配的地形图包含了山脊线、山谷线、鞍部、陡坡和缓坡5类。

1.2 流向算法分配

现有的单流向算法中应用最多的是D8，该算法为每个栅格单元指定其较低的邻居栅格单元中最陡坡向为水流方向。D8在模拟山谷、山脊等陡峭区域的径流过程中能够获得较为理想的效果，且相较于其他单流向算法，其算法更加简单，易于实现，计算效率高。这里沿用以往做法，采用2的整数次幂来标记流向，8个流向分别使用1、2、4、8、16、32、64、128，从当前栅格单元左侧栅格开始以顺时针方向的不同流向进行标识，如图2(a)所示。

在缓坡以及鞍部区域，采用MFD-fg算法建立径流模型。MFD-fg以最大下坡坡度的线性函数对水流分配进行加权，水流分配策略随着地形变化而变化，使得陡峭的方向分配的水流量越多，平缓的方向分配的水流量越少。通过这种自适应的水流分配算法计算得到的汇水面积(本文用汇水区栅格数代表汇水面积)与D8算法得到的值衔接相对自然，因此选用MFD-fg算法模拟平缓地区的水流运动。MFD-fg的水流分配策略由式(1)定义

(1)

式中，di为第i个邻域栅格的水流分配比例；tanβi为坡度；f(e)(f(e)>0)为流量分配权重，其中e代表的是最大坡度值；emin、emax分别表示的3×3窗口中e的最大值和最小值[11]。Li为第i个相邻栅格的等高线长度加权因子，如式(2)所示

(2)

混合流向算法采用与单流向算法一致的流向标记方式。如图2(b)所示：具有多个流向的栅格，其流向值由代表各个方向的数值之和表示，如：流向值为3的栅格其水流方向指向东和东南，流向值为255的栅格其水流方向指向周围8个栅格。

图2 混合流向值示意图Fig.2 Schematic plot of hybrid flow direction value

2 试验与结果

本文采用了2000年“航天飞机雷达地形测量计划”获取的DEM数据，分别在山西省永和县的黄土地貌区和重庆市巫溪县的中山切割地貌区各选择了一个试验样区，二者均不位于城区，研究区地貌特征如图3所示，图3(a)样区面积为14.9 km2，平均高程为1 109.1 m，平均坡度为18.6°。图3(b)样区面积为43.6 km2，平均高程为1 399.9 m，平均坡度为32.6°。使用的DEM空间分辨率约3″和1″。

图3 研究区概况Fig.3 The general situation of study area

首先对原始SRTM数据进行投影转换，1″和3″数据分别被转换为高斯投影下30 m和90 m的DEM。考虑到SRTM高程数据为地表高程模型(digital surface model)，同时在试验区内存在低矮建筑、植被等信息，为了减少噪声影响，进行了高斯滤波处理。DEM中的洼地和平坦地区会导致栅格没有流向或者流向不明确，因此，进行了洼地填充处理。

接着，利用几何形态法提取了山谷线和山脊线。根据自然断点分级法，确定坡度阈值。永和样区30 m分辨率的DEM的坡度阈值为17.2°，90 m分辨率DEM的坡度阈值为11.6°；巫溪样区30 m分辨率的DEM的坡度阈值为30.7°，90 m分辨率DEM的坡度阈值为27.2°。然后，通过特征线和坡度阈值将流域内的地形划分为：山脊、山谷、鞍部、陡坡和缓坡5类，对山谷、山脊以及陡坡区域采用单流向的D8算法，对鞍部和缓坡地区采用多流向的MFD-fg算法。如图4所示，在地形分类图中，白色栅格表示山脊，黑色栅格表示山谷，灰色系表示坡地和鞍部；在流向算法分配栅格中，灰色代表单流向算法区，黑色代表多流向算法区。

图4 分类结果Fig.4 The result of classification

分别利用单流向算法D8、Rho8，多流向算法FD8、MFD-fg，以及混合流向算法对两个试验区计算了地表各点汇水面积。如图5和图6所示，D8算法和Rho8算法得到的集水区面积在缓坡地区生成了大量与实际不符的单向径流和平行直线。FD8算法和MFD-fg算法得到的结果在山脊、山谷和陡坡地区产生了不合理的漫散分流的水流运动，MFD-fg算法虽在水量分配上有改进，但依然没有克服此类水流运动在陡坡区域的不合理现象。混合流向算法结合了单流向算法能较好模拟陡峭区域水流运动的特点，以及多流向算法在平缓坡区模拟水流漫散流动的优势，相比于另外4种算法得到的汇水面积栅格图，可以看到陡峭区域水流汇聚，径流宽度相对变窄；平缓坡区的径流分流，根据地形特点水量合理分配到中心栅格周围的低高程栅格。混合流向算法更加合理地反映了地形变化对水流运动以及水量分配的影响。两个样区30 m分辨率的汇水面积结果相比于90 m的结果都更加详细地反映了复杂地形的水量分配。

图5 永和样区5种流向算法的汇水面积图Fig.5 Catchment area results of five algorithms in Yonghe area

图6 巫溪样区5种流向算法的汇水面积结果Fig.6 Catchment area results of five algorithms in Wuxi area

3 分析与讨论

为了进一步验证混合流向算法在传统算法上的改进，分别采用了汇水面积差异图、累积频率图、散点分布图矩阵来定量化分析和描述混合流向算法与传统算法的相似性和差异性。

3.1 汇水面积差异图

用混合流向算法分别与另外4种算法得到的汇水面积进行相减计算，将得到的汇水面积差值结果进一步分层设色表示(见图7、图8)。

图7 永和样区汇水面积差异Fig.7 The difference of catchment area in Yonghe area

图8 巫溪样区汇水面积差异Fig.8 The difference of catchment area in Wuxi area

因为汇流区域的上游栅格数量大，累积了所有上游栅格汇水面积的变化量，所以得到的变化较大的地方主要集中在沟谷区域。4种传统算法与混合流向算法相减的差异从小到大依次排序是：D8、MFD-fg、Rho8、FD8。对比不同分辨率DEM，30 m分辨率DEM计算的汇水面积差值明显大于90 m分辨率DEM，因为高分辨DEM的栅格数多，相应位置汇流区域得到的汇水面积值大于低分辨率DEM；并且高分辨率DEM中的地形分类更加复杂，混合流向算法相应的流向分配策略更加详细，因此在高分辨率栅格中，混合流向算法与其他算法的差异体现得更加明显。

3.2 累积频率

由于汇水面积数值区间较大，为了方便对比，首先对其进行对数变换，在此基础上统计累积频率[23]。图9为两个样区不同分辨率DEM计算所得的汇流面积累积频率分布图。累积频率分布图在ln(catchment area)≤6.4时趋于一致，约占总数的98%。图中单流向算法和多流向算法可以明显区分开来，混合流向算法由于在陡峭地区采用了单流向算法，平缓区域采用了多流向算法，从累积频率图上可以看出，D8和Rho8中汇水面积较小的栅格占比大于多流向栅格，其原因是单流向算法的下游区域完全接受上游汇入，而多流向算法中较高区域径流向多个方向扩散。混合流向累积曲线介于D8和MFD-fg之间，说明汇水面积分布介于单流向算法和多流向算法结果之间，能较好地兼顾不同区域水流汇聚和坡面散流的特性。

3.3 散点图

图10是在30 m和90 m DEM中分别使用4种传统流径算法和混合流向算法所得到的汇水面积的散点图。

在同样分辨率DEM上，混合流向算法的汇水面积值与D8和MFD-fg的分布接近，相关系数最高，其次是FD8和Rho8。4种算法与混合流向算法在较大值处都有部分一致分布。对比不同分辨率DEM，30 m分辨率的DEM上混合流向算法与其他4种算法的相似性低于90 m分辨率的DEM。因为DEM的分辨率越高，所能刻画的地形细节越丰富，混合流向算法的流向分配方案适用性更强，所以得到结果与其他算法的差异比低分辨率DEM得到的结果差异更大，与累积频率分析一致。也就是说，在更高分辨率的DEM上，混合流向算法的改进效果更加明显。对比两个地貌试验区的结果，中山切割地貌水流汇聚特征得到了反映，抑制了多流向算法结果的散流现象，而在黄土地貌坡度平缓区域单流向造成的漫流特征也得到了加强。

图9 ln(catchment area)累积频率分布图Fig.9 Cumulative frequency distribution of ln(catchment area)

图10 ln(catchment area)散点分布图Fig.10 Scatter plot matrix for ln(catchment area)

4 结 论

流向算法直接影响了汇水面积的计算结果，是汇水网络提取以及汇水单元划分的关键步骤，对研究分布式水文模型具有重要的意义。多年来，研究者们提出了各种模拟水流运动以及流量分配的方法。单流向算法虽然擅长模拟沟谷等地形中的水流的汇流机制，但在模拟坡面水流时会受到单流向的局限，无法表现漫流情景，导致结果中存在大量平行流现象。多流向算法中的水能流向所有高程较低的方向，弥补了单流向算法模拟坡面漫流的缺陷，但也失去了单流向算法模拟陡峭地形水流运动优势。

本文研究提出的混合流向算法实现了单流向算法和多流向算法的结合，充分考虑了各种地形变化对水流运动以及水量分配的影响，对比结果显示，该算法得到的汇水面积更加合理，并且在相对较高分辨率的DEM中能够得到更加明显的改进。

地形类别的划分对于流向算法的分配有着至关重要的影响，在后续的工作中，将对地形分类和坡度阈值的判断选择作进一步研究，以便改进算法。

参考文献：

[1] 卢庆辉, 熊礼阳, 蒋如乔, 等. 一种融合Priority-Flood算法与D8算法特点的河网提取方法[J]. 地理与地理信息科学, 2017, 33(4): 40-46.

LU Qinghui, XIONG Liyang, JIANG Ruqiao, et al. A Drainage Networks Extraction Algorithm by Integrating the Characteristics of Priority-flood and D8 Algorithm[J]. Geography and Geo-Information Science, 2017, 33(4): 40-46.

[2] O’CALLAGHAN J F, MARK D M. The Extraction of Drainage Networks from Digital Elevation Data[J]. Computer Vision, Graphics, and Image Processing, 1984, 28(3): 323-344.

[3] FAIRFIELD J, LEYMARIE P. Drainage Networks from Grid Digital Elevation Models[J]. Water Resources Research, 1991, 27(5): 709-717.

[4] PAIK K. Global Search Algorithm for Nondispersive Flow Path Extraction[J]. Journal of Geophysical Research: Earth Surface, 2008, 113(F4): F04001.

[5] SHIN S, PAIK K. An Improved Method for Single Flow Direction Calculation in Grid Digital Elevation Models[J]. Hydrological Processes, 2017, 31(8): 1650-1661.

[6] ZHOU Qiming, PILESJÖ P, CHEN Yumin. Estimating Surface Flow Paths on a Digital Elevation Model Using a Triangular Facet Network[J]. Water Resources Research, 2011, 47(7): W07522.

[7] FREEMAN T G. Calculating Catchment Area with Divergent Flow Based on a Regular Grid[J]. Computers & Geosciences, 1991, 17(3): 413-422.

[8] QUINN P, BEVEN K, CHEVALLIER P, et al. The Prediction of Hillslope Flow Paths for Distributed Hydrological Modelling Using Digital Terrain Models[J]. Hydrological Processes, 1991, 5(1): 59-79.

[9] QUINN P F, BEVEN K J, LAMB R. The in(a/tan/β) Index: How to Calculate It and How to Use It within the Topmodel Framework[J]. Hydrological Processes, 1995, 9(2): 161-182.

[10] PILESJÖ P, ZHOU Qiming, HARRIE L. Estimating Flow Distribution over Digital Elevation Models Using a Form-based Algorithm[J]. Geographic Information Sciences: A Journal of the Association of Chinese Professionals in Geographic Information Systems, 1998, 4(1-2): 44-51.

[11] 秦承志, 李宝林, 朱阿兴, 等. 水流分配策略随下坡坡度变化的多流向算法[J]. 水科学进展, 2006, 17(4): 450-456.

QIN Chengzhi, LI Baolin, ZHU Axing, et al. Multiple Flow Direction Algorithm with Flow Partition Scheme Based on Downslope Gradient[J]. Advances in Water Science, 2006, 17(4): 450-456.

[12] GALLANT J C, WILSON J P. Primary Topographic Attributes[M]∥WILSON J P, GALLANT J C. Terrain Analysis: Principles and Applications. New York: John Wiley & Sons, 2000: 51-85.

[13] 张占阳, 郑红梅, 郭新成, 等. 局部曲面正切曲率与流线曲率的水流路径算法[J]. 测绘科学, 2013, 38(2): 122-124.

ZHANG Zhanyang, ZHENG Hongmei, GUO Xincheng, et al. Flow Paths Algorithm Based on Tangential Curvature and Flow-path Curvature of Partial Surface[J]. Science of Surveying and Mapping, 2013, 38(2): 122-124.

[14] XIONG Liyang, TANG Guoan, YAN Shijiang, et al. Landform-oriented Flow-routing Algorithm for the Dual-structure Loess Terrain Based on Digital Elevation Models[J]. Hydrological Processes, 2014, 28(4): 1756-1766.

[15] JOHNSTON E G, ROSENFELD A. Digital Detection of Pits, Peaks, Ridges, and Ravines[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1975, SMC-5(4): 472-480.

[16] LINDSAY J B. Sensitivity of Channel Mapping Techniques to Uncertainty in Digital Elevation Data[J]. International Journal of Geographical Information Science, 2006, 20(6): 669-692.

[17] JASIEWICZ J, STEPINSKI T F. Geomorphons —A Pattern Recognition Approach to Classification and Mapping of Landforms[J]. Geomorphology, 2013, 182: 147-156.

[18] BURROUGH P A, VAN GAANS P F M, MACMILLAN R A. High-resolution Landform Classification Using Fuzzy K-means[J]. Fuzzy Sets and Systems, 2000, 113(1): 37-52.

[21] 闾国年, 钱亚东, 陈钟明. 基于栅格数字高程模型自动提取黄土地貌沟沿线技术研究[J]. 地理科学, 1998, 18(6): 567-573.

LU Guonian, QIAN Yadong, CHEN Zhongming. Study of Automated Extraction of Shoulder Line of Valley from Grid Digital Elevation Data[J]. Scientia Geographica Sinica, 1998, 18(6): 567-573.

[22] 张尧, 樊红, 李玉娥. 一种基于等高线的地形特征线提取方法[J]. 测绘学报, 2013, 42(4): 574-580.

ZHANG Yao, FAN Hong, LI Yu’e. A Method of Terrain Feature Extraction Based on Contour[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(4): 574-580.

[23] 刘学军, 晋蓓, 王彦芳. DEM流径算法的相似性分析[J]. 地理研究, 2008, 27(6): 1347-1357.

LIU Xuejun, JIN Bei, WANG Yanfang. Similarity Analysis of Flow Route Algorithms for Extracting Drainage Network from Grid-based Terrain Model[J]. Geographical Research, 2008, 27(6): 1347-1357.

[24] JENSON S K, DOMINGUE J O. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis[J]. Photogrammetric Engineering and Remote Sensing, 1988, 54(11): 1593-1600.

[25] 李辉, 陈晓玲, 张利华, 等. 基于三方向搜索的DEM中洼地处理方法[J]. 水科学进展, 2009, 20(4): 473-479.

LI Hui, CHEN Xiaoling, ZHANG Lihua, et al. Depression Removal Method for Grid DEM Based on Three-direction Search[J]. Advances in Water Science, 2009, 20(4): 473-479.