EN中文

地球信息科学学报 >

2015 , Vol. 17 >Issue 1: 15 - 21

DOI: https://doi.org/10.3724/SP.J.1047.2015.00015

一种以等差分级的DEM填洼算法

  • 高翔 ,
  • 蔡国林 ,
  • 徐柱 , * ,
  • 蔡伟娣
展开
  • 西南交通大学遥感信息工程系,成都 611756
*通讯作者：徐 柱(1972-),男,博士,教授,研究方向为空间分析、复杂网络等。E-mail:

作者简介：高 翔(1989-),男,江苏盐城人,硕士生,研究方向为地图制图学与地理信息系统。E-mail:

收稿日期: 2014-02-11

  要求修回日期: 2014-04-02

  网络出版日期: 2015-01-05

基金资助

教育部“新世纪优秀人才支持计划”(NCET-12-0942)

“2011计划”轨道交通安全协同创新中心西南交通大学先行先试项目

国家测绘地理信息局重点实验室开放基金项目(DM2013SC03)

收起

A New Algorithm to Fill Depressions in Digital Elevation Model Based on Arithmetic Arrays

  • GAO Xiang ,
  • CAI Guolin ,
  • XU Zhu , * ,
  • CAI Weidi
Expand
  • Department of RS and GIS, Southwest Jiaotong University, Chengdu 611756, China
*Corresponding author: XU Zhu, E-mail:

Received date: 2014-02-11

  Request revised date: 2014-04-02

  Online published: 2015-01-05

Copyright

《地球信息科学学报》编辑部 所有
Fold

摘要

DEM的填洼是水系提取中最耗时的过程。真实的地表低高程点通常在空间上相邻,在对低高程栅格进行填洼时,少数栅格的高程值更新会导致级内大量栅格的循环迭代,从而消耗大量运算时间。为了提高DEM的填洼效率,在分级填洼的基础上,提出一种等差分级填洼算法,该算法顾及DEM各级填洼时间呈幂函数逐级递减这一地形规律。算法首先创建一系列数组,数组大小随等差数列递增;然后依次将排序后的DEM数据存放至数组中,由于定义的数组大小逐渐增大,因此,低高程区域被“分割”得很细;最后,当填洼运算从级内转至数组内时,低高程区域细致的划分极大缩短了栅格循环迭代时间,从而使得算法获得了较高的效率。新算法既能有效缩短填洼过程中的迭代时间,又能保证所提取水系的完整性与连贯性。为验证算法的有效性,选用四川省不同地区的SRTM 90 m分辨率DEM作为实验数据,并与已有的DEM分级填洼算法进行对比分析。实验结果显示：当研究区域栅格总数达到两千万个时,本文算法填洼效率提升了50%左右,且随着DEM数据量的增大,算法效率的提升更为明显,与此同时,利用新方法进行填洼后,DEM水文线连续性较强,表明了新方法的有效性。

关键词： DEM; 洼地填充; 等差分级

本文引用格式

高翔 , 蔡国林 , 徐柱 , 蔡伟娣 . 一种以等差分级的DEM填洼算法[J]. 地球信息科学学报, 2015 , 17(1) : 15 -21 . DOI: 10.3724/SP.J.1047.2015.00015

Abstract

Depressions cause errors or breaks when extracting rivers from DEM. Filling the depressions in the initial DEM is important. Low elevated surface areas in space are usually adjacent, and high elevated surface areas in space are always discrete. Therefore, when the grids of low elevated areas are filled, a large number of grids in low elevated areas will have to be recalculated, which makes the process of depression filling more complex. This paper investigates a new ranking algorithm for improving the efficiency of filling the depressions of DEM. First, considering the fact that the depression filling time in different pieces decreases in the form of power function, the algorithm creates a list of variable arrays. Then the data of DEM, which has been ranked using quicksort algorithm, is put into the arrays successively. Finally, grids are filled within the variable arrays, whose lengths are increased with respect to arithmetic sequence. As the lengths of the arrays are increased, the low elevated areas are divided finely. With the increase of elevation, the arrays will contain more numbers of grids, which will reduce the depth of seeking path in iterations for the unfilled grids. The new algorithm keeps the integrity and continuity of the drainage networks extracted by the pits-removed DEM and well reduces the depth of seeking path in iterations to accelerate pits-removing process. To illustrate this approach, the algorithm has been employed in the SRTM DEM at a 90-meter resolution of Chengdu City. The results show a considerable efficiency. The efficiency of the algorithm outperforms the traditional algorithm, as it will improve 50% efficiency compared with previous algorithm when the number of grids of DEM reaches 20 million. And the improvement will be more evident when the amount of DEM data increases.

Key words： DEM (Digital Elevation Model); depression filling; arithmetic arrays

1 引言

洼地在DEM栅格数据中大量存在,在对DEM数据实施水系提取时,洼地使得水流不能正常流出,因此提取的水系存在较大的误差甚至发生断裂。为准确提取DEM中的自然水系,有时需要对初始DEM数据进行洼地填平[1-7]。
许多学者提出了各自的DEM填洼算法：Marta等建议将闭合洼地当作池塘来处理,用水填满直到溢出;Jenson等认为洼地需要进行填平同时平地需要增高,并提出了相应算法;Marta等也提出类似算法,并将洼地分为凹陷型和阻挡型。上述算法都能对洼地进行相应处理,然而其运算效率并不高。Moran与Vezna在1993年提出了一种快速填洼的思想(以下简称M&V算法),Planchon等在2001年对该思想进行了编程实现,并证实了该方法具有简单、高效的特点[8-14]。后人对M&V作了改进[15],发现先对DEM实施基于高程大小的均匀分级再依级填洼可加速M&V算法[16-18]。
分级填洼算法虽简单、高效,但实验数据表明算法的各级填洼时间存在极大的不均匀性,具体表现为：低高程区域的多数栅格在空间上相邻,其级内填洼时间远大于较高区域所需计算时间。为进一步提升填洼效率,本文对分级过程进行了优化,提出了一种等差分级填洼算法,实验证明此方法较传统分级方法有更高的填洼效率。

2 基于等差分级的DEM填洼算法

2.1 基于M&V算法的高程分级填洼

M&V算法思想的提出者是Moran和Vezina,之后Planchon等依据这一思想实现了该算法。算法使用高水面完全淹没除DEM边界以外的所有待填洼区域,并通过遍历将边界高程值传递至其余栅格。一次遍历即可完成填洼操作,但填洼后的DEM与初始DEM偏差较大,因此算法需通过循环迭代无限逼近初始地表。
M&V填洼算法将DEM看作纯粹的数字矩阵单元,在处理时将所有栅格单元加入到算法中一次性求解,而算法的循环迭代过程往往需要消耗大量时间,这一点在DEM数据庞大时显得尤为突出,因此,为提升填洼效率,算法引入栅格分级策略,即将栅格划分为不同区域,每个区域进行各自的填洼计算。
为进一步提升分级填洼效率,从地形角度对算法进行研究。图1中的等高线示意图表明：地形中高程值介于500~1000 m之间的栅格数最多（低高程栅格）,且这些栅格在空间上连续（未被隔断）;高程值介于1000~1500 m之间的栅格数开始变少,且这些栅格被红线分割;随着高程值的增大,栅格数逐渐减少,且在空间上也逐渐离散。图1中等高线模拟的地形仅包含2个山体,地形变化过于简单,而现实地表通常包含大量山体,因此,低高程区域的大量栅格在空间上连续,随着高程值的增大,栅格逐渐离散的地形特征也将表现的更加明显,而这一地形特征降低了传统分级算法的填洼效率,因为在对低高程栅格进行填洼时少数栅格的高程值更新会导致级内大量栅格的循环迭代,从而消耗大量运算时间,随着高程值的增大,栅格在空间上逐渐离散最终彼此独立（如DEM中的大量山顶）。因此,局部栅格的高程值更新不会引起其余栅格的循环运算,最终所需的填洼时间也远小于低高程区域,这导致了传统分级填洼算法在效率上呈现出前低后高的不合理趋势。
Fig. 1 Schematic diagram of contours

图1 等高线示意图

2.2 一种等差分级的填洼算法

为解决这一缺陷,本文提出了一种等差分级的填洼算法。首先,创建一系列大小以等差数列变化的数组;接着依次将排序后的DEM数据存放至数组中,由于定义的数组大小逐渐增大,因此,低高程区域被“分割”得很细;最后,低高程区域细致的划分极大缩短了栅格循环迭代时间,从而使得算法获得了较高的效率。假设原始DEM数据包括S个待填洼栅格,用N个数组对其进行存储,则第i级的数组大小Mi可由式（1）计算得到。
M i = i × 2 S N × ( N + 1 ) （1）
随着高程的增加,数组包含的栅格逐渐增加,覆盖的区域范围也随之增大,因此,未能填洼的遗留栅格能够很快进入高程值较大的区域进行计算,这很好地解决了传统的分级算法遗留栅格迭代效率低下的缺陷。如图2所示,传统分级填洼的1号区域被1、2、3号数组代替,原存放于一级中的所有栅格被分割至1、2、3号数组中,这一过程减少了级内栅格循环次数,因此,原1号低高程区域的填洼时间被极大的缩短。图2中红色区域表示DEM填洼时需迭代的栅格,当未处理栅格被DEM中最高区域环绕时,若采用传统的分级填洼算法,红色区域中的栅格只有迭代至最高区域（N号区域）才能得以填洼;若采用大小随等差数列增加的数组进行填洼计算,由于数组大小随高程的增加而扩大,因此,未处理栅格能够更快地进入高程较大的区域,在图2中具体表现为数组的N号区域比传统分级的N号区域包含更多的栅格总数。数组的分级策略有效缩短了未处理栅格（红色区域）的迭代路径,从而有效地提升了填洼效率。
Fig. 2 Ranking process of arithmetic arrays

图2 等差分级过程

等差分级填洼算法的处理流程为：
（1）对栅格以高程值大小进行排序;
（2）确定分级数,用式（1）创建等差分级数组;
（3）将排序后的栅格依次填充至各级数组中,并从最低级别开始实施填洼;
（4）确定级内遗留栅格（未填洼栅格）,将其迭代至下一级别内并参与该级的填洼运算;
（5）反复运算第（3）、（4）步骤,直至所有栅格完成填洼计算。

3 基于DEM的填洼算法实验与分析

3.1 传统分级算法效率分析

本文选取成都市某区域SRTM 90 m分辨率的原始高程数据（图3）对传统分级填洼,以及等差分级填洼算法进行测试,所选区域共包括1000×1000=1 000 000个栅格点。
Fig. 3 Original elevation data of SRTM of 90 m resolution

图3 SRTM 90 m分辨率的原始高程数据

算法运行平台为MATLAB,测试计算机配置为2.0 GHz CPU,2 GB内存。将分级数定为20,表1显示了算法运行时各级所需填洼时间及循环次数。
Tab. 1 Time cost of filling depressions of ranked algorithm(s)

表1 高程分级填洼耗时(s)

数 时间（s） 循环次数 级数 时间（s） 循环次数
1 22.86 158 11 1.29 15
2 5.21 58 12 1.11 16
3 8.22 83 13 0.80 11
4 6.01 62 14 0.88 14
5 3.71 40 15 0.82 12
6 3.21 26 16 0.90 10
7 3.53 22 17 0.77 7
8 3.91 25 18 0.58 5
9 2.36 19 19 0.54 4
10 3.23 33 20 0.49 4
图4描述了表1中数据的变化趋势,趋势表明填洼算法所需时间集中在低高程区域,随着高程的增加,填洼时间及级内循环次数迅速减少,如第1级需158 s的填洼时间而第20级仅需4 s即可完成。最终发现,填洼时间呈幂函数逐级递减,函数方程如式（2）所示。
y = 26.408 x - 1.23 （2）
Fig. 4 Curve fitting of depression filling time in pieces

图4 级内填洼时间曲线拟合

幂函数拟合的相关系数R2值等于0.89,说明用幂函数对级内填洼时间序列进行拟合取得了很好的效果。
各级填洼时间呈幂函数递减是因为分级填洼算法在对DEM进行高程分级时,将全部栅格平均分配到各个级别内。然而,真实的地表低高程点通常在空间上相邻,在对低高程栅格进行填洼时少数栅格的高程值更新会导致级内大量栅格的循环迭代,从而消耗大量运算时间,随着高程值的增大,栅格在空间上逐渐离散最终彼此独立,因此,局部栅格的高程值更新不会引起其余栅格的循环运算,最终所需的填洼时间也远小于低高程区域。

3.2 等差分级填洼算法效率随分级数的变化规律

本文选取6组四川省SRTM 90 m分辨率的原始高程数据对等差分级填洼算法进行测试,6组DEM的栅格总数如表2所示。为研究算法效率随分级数的变化规律,测试6组DEM在不同分级数下的填洼时间,结果如表3所示。
Tab. 2 Total number of DEMs from six groups

表2 6组DEM的栅格总数汇总

区域1 区域2 区域3 区域4 区域5 区域6
栅格总数（个） 836 740 2 143 750 4 821 320 8 575 000 13 394 910 19 275 000
Tab. 3 Variation between depression filling time and number of the pieces

表3 填洼时间随分级数的变化规律

分级数 区域1 区域2 区域3 区域4 区域5 区域6 总耗时(s)
1 17.53 57.38 187.57 317.86 973.42 2876.58 4430.34
3 9.47 34.62 117.32 193.54 383.60 1134.52 1873.07
9 4.33 11.20 51.14 136.92 314.77 753.69 1272.05
27 1.12 6.54 23.60 87.11 253.94 674.53 1046.84
81 0.49 4.73 17.38 59.29 234.50 613.97 930.36
243 0.13 4.11 14.26 52.05 227.65 603.41 901.61
729 0.14 4.27 15.80 51.73 259.80 587.32 919.06
2187 0.96 4.18 15.66 55.83 232.19 697.62 1006.44
6561 1.50 7.92 18.73 62.40 296.34 774.36 1161.25
图5进一步对填洼时间随分级数的变化规律进行了描述。图5（a）表明填洼时间随分级数的变化呈现出3个阶段：（1）当分级数小于100时,随着级数的增加,填洼时间迅速减小,这是因为分级数的增加有效地控制了栅格迭代次数;（2）分级数介于100-1000时,填洼时间随级数的变化趋于稳定;（3）分级数大于1000时,随着级数的增加,填洼时间有一定的增加,这是因为当级数设定得太大时,填洼时会因为栅格迭代次数过多（低级别内的大量遗留栅格需要向高级别传递）而一定程度地减低了填洼效率。由于实验的分级数跨度较大,为更好地研究第一阶段填洼时间随分级数的变化规律,不考虑分级数大于100的数据,并得到相应的变化折线图,如图5（b）所示。通过对比不同拟合函数的相关系数发现,幂函数最能够刻画该阶段填洼时间随级数的变化规律,拟合函数为 y = 3342.4 x - 0.337 ,相关系数为R2=0.8632,因此,该拟合函数可作为预测函数对不同级数的填洼耗时进行估计。同时,考虑到较多的分级数会给计算机造成较大的内存开销,故在兼顾算法效率的同时,分级数越小越好。图5（b）中拟合函数与变化关系折线（蓝色实线）交点的级数为35,当级数超过35时,填洼时间随级数的变化趋于平缓。为此,在进行分级填洼时,可将35作为默认的最优分级数。
Fig. 5 Variation between depression filling time and number of the pieces

图5 填洼时间随分级数的变化图

3.3 等差分级填洼算法效率分析

对不同地区的DEM数据以相同的计算机配置进行填洼计算（表4）,分级数统一确定为243（虽不是最优分级数,但该级数可获得最小的填洼总耗时）。等差分级填洼与传统分级填洼总耗时随DEM数据量的增大呈现相似的变化趋势,用二阶多项式与指数函数对填洼时间序列进行拟合（表5）,拟合的相关系数R2分别为0.9984与0.9774,故而拟合的多项式可作为填洼时间随栅格总数变化的预测方程。
Tab. 4 Comparison of efficiency between the new algorithm and the traditional algorithm

表4 新算法与传统算法效率对比(SRTM 90 m分辨率)

区域 等差分级填洼耗时（s） 传统分级填洼耗时（s） 新方法栅格迭代总数 传统分级栅格迭代总数 新算法与传统算法时间差（s）
1 0.13 1.49 13 860 13 736 1.36
2 4.11 4.66 82 734 86 063 0.55
3 14.26 15.03 340 983 379 001 0.77
4 52.05 56.63 2 161 678 2 442 419 4.58
5 227.65 238.64 11 484 804 12 599 816 10.99
6 603.41 1207.84 61 173 465 62 142 823 604.43
Tab. 5 Prediction equations of depression filling time

表5 填洼耗时预测方程

洼算法 多项式拟合方程 R2
传统分级填洼 y = 2.0127e3E-07x 0.9774
等差分级填洼 y = 2E-12x2 - 2E-05x + 23.388 0.9984
等差分级与传统算法填洼耗时对比如图6所示。图6（a）中的填洼耗时对比表明：当栅格总数小于500×104个时,等差分级填洼算法与传统算法相比效率无明显提升,算法运算时间差在1 s左右;当栅格总数在（500~1000）×104个时,等差分级填洼效率较传统算法效率差异仍很小,填洼时间差在5 s内;当栅格总数在（1000~1500）×104个时,等差分级填洼算法与传统算法相比填洼时间差产生了一定的差异,时间差在11 s内;当栅格总数大于1500×104个时,两种算法的填洼效率产生了巨大的差异。如当栅格总数接近2000×104个时,2种算法的填洼时间差达到了600 s。图6(b)的填洼耗时拟合曲线表现出了相同的规律：2条拟合曲线在1500×104处相交,当栅格总数小于1500×104个时,两条曲线时间差很小,当栅格总数大于1500×104个时,随着DEM数据量的增大,等差分级填洼算法较传统算法的优势也急剧增大。
Fig. 6 Comparison of depression filling time between variable arrays based algorithm and the traditional algorithm

图6 等差分级与传统算法填洼耗时对比

通过分析可知,基于等差分级的填洼算法效率总是优于传统的分级填洼算法,且随着DEM栅格总数的增加,二者所需的填洼时间差也急剧增大。实验数据总区域最大的DEM总栅格数为1.9275×107个,传统算法的填洼时间为1207.84 s,等差分级填洼时间为603.41 s,填洼总时间缩短了50%左右。故以等差分级的填洼算法较传统算法在效率上有较大的提升。

3.4 填洼结果分析

使用等差分级方法对DEM实施填洼,其中区域最大（第6组数据）的DEM填洼结果如图7所示。图7（a）表明,若不对初始DEM实施填洼,提取出的水文线连续性差,无法满足后续的分析应用需求;图7（b）表明,等差分级方法能对DEM实施正确的洼地填平,以此获得的水文线也非常连续,能够满足后续的应用需求。图7（c）为DEM中被填洼的部分(白色区域)。在实际地形中,洼地存在伪洼地与天然洼地2种,而算法在填洼过程中并未考虑这一点,真实存在的洼地也被填平。因此,在填洼完成后应当将填洼区域与实际地形做比较,若填平的区域确实为实际存在的洼地,则应当将该区域恢复为初始高程。
Fig. 7 Hydrological analysis based on DEM

图7 基于DEM提取水文

4 结论

利用DEM作水文特征分析,最重要的一步是进行区域洼地填平,其中以M&V算法的分级填洼能简单、高效地实现这一过程。分级填洼既能有效缩短填洼的搜索路径,又使得后续能够完整、连贯地提取河网。然而由于地表低高程区域的栅格通常在空间上相邻,在对其进行分级填洼时部分栅格的高程更新会导致大量栅格的循环运算,因此,算法的填洼过程几乎集中在低高程区域。
为解决这一填洼算法效率的不均匀性问题,本文提出以等差分级的DEM填洼算法,将DEM数据保存至长度随等差数组递增的数组中,从而细分低高程区域而粗分高程较大区域,这一方法可解决传统分级的2个缺陷：（1）将运算平均到各级数组中从而减小总的填洼时间;（2）进一步缩短未处理栅格的迭代搜索路径。实验证明,等差分级填洼算法比传统分级算法更加高效,且随着DEM数据量的增大,新算法的优势也逐渐扩大,同时,用新方法进行DEM填洼后,可获得有较强连续性的水文线。这证明了新方法的有效性。

The authors have declared that no competing interests exist.

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
周启鸣,刘学军.数字地形分析[M].北京:科学出版社,2006:1-327.

[2]
O'Callaghan J F, Mark D M. The extraction of drainage networks from digital elevation data. Computer Vision[J]. Graphics and Image Processing, 1984,28:323-344.

[3]
张行南,井立阳,叶丽华,等.基于数字高程模型的水文模拟对比分析[J].水利学报,2005,36(6):759-763.

[4]
刘学军,卢华兴,卞璐,等.基于DEM的河网提取算法的比较[J].水利学报. 2006,37(8):1022-1028.

[5]
Tarboton D.A new method for the determination of flow directions and contributing areas in grid digital elevation models[J]. Water Resources Research, 1997,33:309-319.

[6]
Garbrecht J, Martz L W.The assignment of drainage direction over flat surfaces in raster digital elevation models[J]. Journal of Hydrology, 1997,193:204-213.

[7]
王建平. 数字流域与数字水文模型的集成研究[D].南京:河海大学,2005.

[8]
Tribe A.Automated recognition of valley lines and drain-age networks from grid digital elevation models: A review and a new method[J]. Journal of Hydrology, 1992,139:263-293.

[9]
Martz L W, Garbrecht J.The treatment of flat areas and depressions in automated drainage analysis of raster digital elevation models[J]. Hydrologic Processes, 1998,12:843-855.

[10]
Wang L, Liu H.An efficient method for identifying and filling surface depressions in digital elevation models for hydrologic analysis and modeling[J]. International Journal of Geographical Information Science, 2006,20(2):193-213.

[11]
Planchon O, Darboux F.A fast, simple and versatile algorithm to fill the depressions of digital elevation models[J]. Catena, 2001,46(2/3):159-176.

[12]
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.

[13]
Martz L W, De Jong E.Catch: A FORTRAN program for measuring catchment area from digital elevation models[J]. Computers&Geosciences, 1988,14(5):627-640.

[14]
Martz L W, Garbrechet J.Numerical definition of drainage network and subcatchment areas from digital elevation models[J]. Computers&Geosciences, 1992,18(6):747-761.

[15]
王建平,任立良,吴益.一种新的DEM填洼处理算法[J].地球信息科学,2005,7(3):51-54.

[16]
杨邦,任立良,贺颖庆.基于快速排序的数字高程模型分级填洼算法[J].计算机应用,2009,29(11):3161-3165.

[17]
于森,任立良.基于DEM模型的新填洼算法[J].地球信息科学学报,2009,11(1):50-55.

[18]
徐精文,张万昌,符淙斌.适用于大尺度水文气候模式的DEM洼地填充和平坦区处理的新方法[J].水利学报,2007,38(12):1414-1420.

[19]
殷人昆,陶永雷,谢若阳,等.数据结构(用面向对象方法与C++描述)[M].北京:清华大学出版社,2004:1-402.

[20]
Hoare C A R. Quicksort[J]. The Computer Journal, 1962,5(1):10-15.

Options
文章导航

/

〈 〉