基于多尺度窗口的DEM局部填洼方法

徐静波，许捍卫，于艳超

基于多尺度窗口的DEM局部填洼方法

徐静波1，许捍卫1，于艳超2

（1.河海大学 地球科学与工程学院，江苏 南京 210098；2.贵州省水利水电勘测设计研究院，贵州 贵阳 550002）

为了去除DEM中的伪洼地，使水系提取结果更加精确，在M&V算法的基础上提出了一种改进的填洼方法。用局部处理的思想将填洼步骤简化，减小对DEM的计算范围和迭代次数，从而达到优化算法的目的；同时借助Arcpy与GDAL库，将其封装为一个模块集成到ArcGIS软件中，以更高效地解决DEM填洼预处理问题。

填洼；算法；多尺度；局部；评估

DEM是以高程值为基础的空间分布模型，用来反映区域的地形特征及趋势[1]。随着地理信息技术的发展，以及其与水文领域的密切结合，可借助一定算法自动提取研究范围内的基础水系。由于插值误差或采样误差，原始DEM会形成许多伪洼地，与实际地形存在不符[2]。洼地是指地表的低洼处，在实际地形中真实洼地是极少的；伪洼地是局部地区的最低点或区域，在水流方向的计算中会导致水不能顺利流出，从而对水系的提取造成障碍，影响结果的准确性。因此，在利用DEM提取水系前，必须进行填洼处理。

O'Callaghan J F[3]等提出用小窗口滤波平滑的方法处理DEM，该方法可填平较小的洼地但不适于大范围的洼地填平。Jenson K[4]等提出J&D算法，主体思想是先填洼地为平地，再逐步垫高并找出潜在出口，确定流向。但这种方法的实现较复杂，对于低分辨率和数据量较小的DEM处理十分有效；对于高分辨率和大数据量的DEM处理就显得效率低下。另外，在ArcGIS软件中，J&D算法经过集成被应用于水文分析的填洼中[5]。Moran C J[6]等于1993年提出了一种快速填洼思想（M&V算法），2001年由Planchon O等编程实现[7]。该方法思路简单，对高低分辨率的DEM普遍适用，且效率提高明显，但由于算法的可移植性和发展的局限性，M&V算法作为一种先进算法尚未得到广泛应用。本文借助ArcGIS平台，利用Python语言和Arcpy站点包，在M&V算法的基础上提出了一种改进算法，并将其封装为模块集成到ArcGIS软件中，以此来改进填洼算法的应用，旨在更高效、更全面合理地解决DEM填洼预处理问题。

1 M&V算法

1.1 M&V算法结构

M&V算法的主要思想是先用一极大值水面淹没原始DEM，再迭代移除DEM上多余的水，得到的结果就是填洼后的数据[8]。M&V算法的实现比较简单，且适用于各种分辨率的DEM。其具体步骤为：

1）除边界网格保留原值外，对剩余网格赋极大值水面，淹没原始DEM。极大值默认选取网格最大水面高程值加1。

2）以3×3窗口逐行对所有网格（边界除外）进行遍历。以网格C为中间网格，遍历其8邻域：将网格C的当前值与邻域网格值作比较，若邻域网格值加上极小变量小于C的当前值，则将网格C重新赋值为邻域网格值加上极小变量。由于填洼是将洼地水面逐渐抬高的过程，因此还需判断网格原始高程与邻域网格值加上极小变量的大小关系，若前者大于后者，则将网格新值重新赋值为原始高程。其中，极小变量值可在算法实现中酌情设定，默认值为0.1。

3）多次迭代，即重复步骤2），直到所有网格进行遍历时，没有可执行的单元为止，迭代结束。此时网格无限逼近于初始DEM，且洼地被填平。

1.2 算法实现

M&V算法充分考虑了水流流向问题，处理后的DEM满足在不改变原始DEM结构的基础上，对其中任意栅格点都能找到一条高程递减的路径将水顺利排出。算法思想简单完善，效率明显高于当前ArcGIS软件中的填洼工具，该工具使用J&D算法。实现M&V算法伪代码为：

For each grid cell c of the DEM '遍历DEM中每个网格

if c is on the border： '边界网格保持值不变

H(c)=Z(c)

else： '中间网格赋予极大值，默认为DEM所有网格最大值+1

H(c)=max(Z(n))+1

Do

loopFlag=False

For each grid cell c of the DEM

if H(c)!=Z(c)： '中间网格

For each existing neighbor n of c： '遍历8邻域

if H(c)＞(H(n)+min)：

H(c)=H(n)+min

loopFlag=True

if Z(c)＞(H(n)+min)：

H(c)=Z(c)

loopFlag=True

Loop While loopFlag=True

2 改进的M&V算法与算法实现

2.1 Python与ArcGIS

Python作为一种免费开源的通用脚本语言，其代码简洁明了，能与众多第三方开源库进行无缝衔接，编译过程非常简单，能直接对源代码进行调试执行，且Python编写的脚本程序能跨平台运行。ArcGIS软件是一款专门用于地理空间数据处理和分析的软件，Esri公司在ArcGIS 10.0中引入了Arcpy。Arcpy是基于Python语言的站点包，提供了一种用于集成ArcGIS处理流程和开发工具脚本的丰富的动态环境。在基于Arcpy开发的ArcGIS应用程序和脚本中，可借助大量Python模块和第三方开源库，使功能实现更加便捷[9]。

2.2 基于多尺度窗口的局部填洼算法及实施步骤

M&V算法充分利用了边界条件和DEM的初始值，先将洼地抬高，再通过多次迭代和与邻域值进行比较，使栅格中的水都能顺利流出，从而达到填洼的目的。迭代是M&V算法中最耗时的过程，其时间消耗主要包括两个方面：迭代次数和每次迭代遍历的栅格数。

本文从精简迭代步骤的角度出发，先提取出洼地网格，找到该洼地的相对出口；再利用M&V算法，以半径R为窗口对每个洼地网格进行计算，使洼地里的水可以顺利流出，即以洼地局部迭代替代M&V算法中的全局迭代。为了最大程度地保留原始地形，需制定一个运算规则来确定洼地的迭代窗口R。基于多尺度窗口的局部填洼算法流程如图1所示。

图1 基于多尺度窗口的局部填洼算法流程图

1）遍历整个网格，标定洼地。洼地分为入流洼地和无向洼地两种。入流洼地以单个网格形式存在，即其高程值小于其8邻域网格的高程值；无向洼地通常以片区形式存在，在DEM中表现为网格高程值不低于该洼地所有邻域网格高程值。

2）确定每个入流洼地的相对出口，此时的搜索半径R即为该洼地M&V算法的迭代窗口大小。以半径为R的窗口（初始值默认为5）遍历洼地的周围网格，若最小值小于洼地高程值，则标定该网格；若最小值仍高于洼地高程，则以步长为2逐步扩大搜索范围，直到搜索到符合要求的相对出口为止。记录找到相应出口时的R值，即为M&V算法迭代窗口的大小。最终洼地的流向就是标定的相对出口，相对出口的标定可作为算法检验的辅助参数。考虑到算法的可行性和合理性，R不可能无限增大，设定R的极限值为Rlim，当R增至Rlim仍未找到相对出口时，则增加洼地高程。

3）无向洼地迭代窗口大小的确定。根据无向洼地的分布，用一个最小外接矩形包住所有的洼地网格，迭代窗口即为最小外接矩形的长宽均增加一个长度后的窗口。

4）使用M&V算法对每个洼地进行迭代。

2.3 基于Python的ArcGIS模块集成

本文开发的填洼算法工具界面如图2所示。

图2 基于多尺度窗口的局部填洼算法工具界面

3 性能评估

3.1 算法复杂度分析

算法复杂度的计算涉及时间和存储两方面的需求，是一种事前估算，分为时间复杂度和空间复杂度。度量方式是指当问题的尺度以某种方式从1递增到n时，解决这个问题的算法损耗的时间或占用的空间也以某种方式从1递增到n，用函数表示为f (n)。本文采用时间复杂度进行度量。根据问题的不同，应考虑最坏的情况和最乐观的情况，使用O表示：只有存在正整数c和n0使得T(n)≤cg(n)对所有的n≥n0成立，那么称算法的渐进时间复杂度为T(n)= O(g(n))[7]。

对于一个大小为N的DEM，J&D算法需遍历所有网格进行填洼，总的时间复杂度为O(N2)；M&V算法每次迭代都至少有一个网格被赋值，迭代最大长度为对角线长度（20.5N），时间复杂度为O(N1.414)；本文提出的先查找洼地再分别填充洼地的做法，在最坏的情况下，才会达到M&V算法的渐进时间复杂度，其时间复杂度接近O(NlogN)。

3.2 运行效率对比

采用3组数据对M&V算法和改进的M&V算法进行对比，90 m分辨率DEM是从地理空间数据云上下载的长江上游数据；30 m分辨率DEM为长江南京段数据，为水下地形测量数据；5 m分辨率DEM为珠江流域CGCS2000成果（表1）。

表1 算法运行效率

3.3 填洼量累计分析

由于M&V算法对高程值的改变有加有减，因此将增加和减少的绝对量之和作为衡量标准。本文算法的3 组数据累计改变高度分别为3 541、46 832和184 035，而M&V算法累计改变高度分别为3 967、50 236、254 793，在保留洼地周围的地形特征上本文算法略优（表2）。

表2 填洼量累计分析结果

4 结 语

本文从DEM填洼算法需求的本源出发，在理解时下通用算法的基础上，对M&V算法进行了改进，提出了一种新的基于多窗口尺度的局部填洼算法。算法考虑到每个洼地的地形特征，并针对区域性地形采用不同尺度的窗口对洼地作处理，减少了原始M&V算法的迭代次数和每一次迭代的数据量，并最大程度地保留了洼地周围的地形特征。实验证明，该方法比传统填洼算法具有更高的填洼效率。

[1] 任立良,刘新仁. 基于DEM的水文物理过程模拟[J].地理研究, 2000,19(4)：369-376

[2] Martz L W, Garbrecht J Numerical Definition of Drainage Network and Subcatchment Areas from Digital Elevation Models[J]. Computers & Geosciences,1992,18(6)：741-761

[3] O'Callaghan J F, Mark D M. The Extraction of Drainage Networks from Digital Elevation Data[J].Computer Vision Graphics and Image Processing,1984,27(3)：323-344

[4] Jenson K, Domingue O. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis[J]. Sensing,1988,54(11)：1 593-1 600

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

[6] Moran C J, Vezina G. Visualizing Soil Surface and Crop Residues[J]. IEEE Computer Graphics and Applications, 1993,13(2)：40-47

[7] Planchon O, Darboux F. A Simple and Versatile Algorithm to Fill the Depressions of Digital Elevation Models[J]. Catena,2002, 46(2)：159-176

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

[9] 彭海波,向洪普.基于Python 的空间数据批量处理方法[J].测绘与空间地理信息,2011,34(4)：81-82

[10] 殷人昆.数据结构(用面向对象方法与C++语言描述)[M].北京：清华大学出版社,1999

P208

B

1672-4623（2017）05-0098-03

10.3969/j.issn.1672-4623.2017.0053.0

徐静波，硕士研究生，主要从事GIS开发与应用。

2015-08-31。

项目来源：国家自然科学基金资助项目（41101374、41101308）；水利部公益性行业科研专项经费资助项目（201201025）。