基于栅格DEM的大河流域洼地填充预处理方法

技术领域

本发明涉及地理信息系统以及分布式水文模型中基于栅格数字高程模型的流域数字水系提取方法，尤其是能够应用于大河流域的洼地填充预处理方法。

背景技术

基于数字高程模型(DEM)提取流域水系信息是实现分布式水文模型数值模拟的重要方法。流域水系信息的提取包括计算流向、累积流向、水系分级以及汇水网络等。目前常用八向方法(简写为D8方法)来获得上述信息，然而这些信息能否成功提取还取决于对洼地及平地像元的预处理。目前已有许多可用于这些任务的软件工具，如著名GIS商业软件ArcGIS中的空间分析工具，ENVI的插件River Tools，Global Mapper软件等，这些软件对于较小的规模或者地理空间范围不大的流域可以成功提取出完整的水系。由于大河流域涉及范围较广，数据量较大，DEM中会存在大量高程误差，以及高程完全相同的“平地”像元，即使经过洼地填充处理的DEM仍难以正确提取出完整的水系。事实证明，对于黄河中上游流域，采用上述软件工具仍无法正确提取，其输出结果是一堆“断头河”。采用Plachon&Doubarx(2002)[ Planchon, O. and Darboux, F., 2002. A fast, simple and versatile algorithm to fill the depressions of digital elevation models. CATENA, 46(2-3): 159-176.]方法及其各种变形，也得到的是类似的“断头河”。Liu [Liu Yonghe, 2009.Another Fast and Simple DEM Depression-Filling Algorithm Based on Priority Queue tructure. Atmospheric and Oceanic Science Letters, 2(4):214-219.]提出了一种洼地填充方法，该方法具有较强的适应性，但对于黄河流域这样的大流域，经过洼地填充同样无法获得完整的水系信息。如何在数据量较大、空间跨度大的DEM中正确提取水系信息仍是一个需要解决的问题。

发明内容

为正确高效的在栅格DEM基础上提取大河流域的水系信息，本发明提供一种以洼地填充步骤为核心的处理方法。其基本思想是借助基于优先队列的洼地填充方法来从DEM边缘最低像元处开始向内逐个“海水淹没”像元，记录每个像元被淹没的次序。由于每个像元可拥有一个不同的次序编号，而且高程愈高的像元，获得的次序编号愈大；高程相同的多个像元，淹没时间愈晚，其次序编号愈大。用淹没次序编号矩阵来代替原DEM高程矩阵来完成流向矩阵、累积流向矩阵、水系分级矩阵，以及其它汇水网络信息。

本发明的技术方案包括以下步骤：

1.创建并初始化一个以浮点数为关键字的优先队列对象；

2.将栅格DEM四条外边界上的像元按照高程从低到高的顺序加入优先队列；

3.当优先队列不空时，继续执行第4步；否则从第12步开始执行；

4.弹出排在最前面的像元，它是优先队列中存放的所有像元中高程最低的，把它作为当前像元，它的高程作为当前海平面高度；

5.借助计数器为当前像元设置淹没次序编号，并将此编号存入淹没次序矩阵中的相应位置上；

6.若当前像元已是淹没状态，则返回第3步执行，否则将其标记为淹没状态后执行第7步；

7.找出当前像元周围8个方向上的所有未被淹没的邻域像元；

8.将低于海平面的邻域像元执行步骤9-10；高于海平面的邻域像元直接加入优先队列，之后返回至第3步；

9.将尚未淹没的邻域像元加入优先队列，并将其标记为临岸状态；

10.当前邻域高程改设为当前海平面的高度；返回第8步将其它邻域像元完成同样的操作；

11.返回至第3步；

12.输出淹没次序矩阵。

本发明的有益效果是基于优先队列的处理方法完成了DEM中填充洼地像元以及无数据像元的同时，且获取了淹没次序号矩阵，用它代替原有DEM高程来完成后续的流向及水系网络信息的提取，可有效避免原DEM中高程相同像元所造成的“平地”像元阻碍水系提取的问题，能获得与真实水系一致的数字水系模型，可用于任何大小的流域提取。

本发明方法已成功的实现了黄河中上游流域的水系正确提取，表明该方法能够正确完成DEM中洼地及无数据区域的填充，并能提取出大河流域的完整水系网络。

附图说明

图1是该发明方法的执行流程。

图2 marks数组示例(1表示未处理像元，0表示被水淹没像元，2表示临岸像元)。

图3是一个小型DEM示例。

图4洼地填充之后的示例DEM(○像元原来为洼地像元)。

图5执行本方法后输出的淹没次序矩阵(即orders数组)。

图6是淹没次序矩阵基础上提取的黄河流域水系。

具体实施方式

    完整的实施流程见图1。

(1)   栅格DEM数据集的准备

栅格DEM数据集中可能存在一些无数据区，需要扫描出这些无数据像元。对由遥感方式获取的DEM数据集，这类像元通常存在于山区或湖泊位置。可以将它们的高程设为一个极小的值，如-9999，这样相当于是将它作为DEM中存在的洼地来对待。在实施本发明方法之后，所有洼地像元被填充，实现了对数据的插补，从而不影响对主要河流的提取。

(2)   优先队列的创建

需要自行编写一个优先队列结构类型，或利用现成优先队列结构类型。优先队列的关键字为浮点型，对应DEM的高程值，每个关键字所关联的存储信息为DEM像元的地址，即行列号位置。优先队列中的关键字以升序排序，处于优先队列中最前面的元素(即像元)对应最小的键值(即高程)。当在优先队列中任意插入一个元素时，应该根据关键字的值来按序插入，以保证优先队列一直保持关键字排序正确。每次从优先队列中弹出的元素总是位于队列最前方、关键字最小的元素。

(3)像元的状态与进出优先队列的关系

像元分初始、淹没、临岸三种状态，分别用值1，0，2来标记。凡是处于优先队列中的像元，对应的状态都应该设置为2，表示处理仍未全部完成；凡是已从优先队列弹出的像元，对应的状态都应该设置为0，这类像元将不再需要处理。当某一像元从优先队列弹出时，意味着它被淹没，这时海面高度应升至该像元的高程上。并且它的那些未被淹没的邻域像元便成为临岸像元。临岸像元若低于海面高度时，其高程需要升至海面高度，这是洼地填充的关键操作。凡是新的临岸像元出现，不仅需要设置临岸状态，还需将其加入优先队列。

(4)初始化数组与计数器

为了标记所有像元的处理状态，需要创建一个与DEM行列数相同的整型数组，下面以名称marks指代该数组； marks是用来存放像元是否被“淹没”状态的数组，其元素值只能为0，1和2。元素值为0时表示像元已被“淹没”，该像元将不需再被处理了；元素值为1时表示像元既没被“淹没”，也不在岸边；元素值为2时表示像元位于岸边，但尚未被“淹没”。marks中所有元素值都需初始化为1。有关marks数组的指示意义可参考图2。

为了记录各像元的淹没次序，还需要创建另一个整型数组，下面以orders指代。为了跟踪淹没次序，还需要再申明一个整型计数器变量DCount，其初始值设为0。每当一个新像元被淹没时，DCount的值就增1，同时将此值写入orders数组中相应当前像元的位置上。这里所谓的“淹没”，是指当一个像元的记录(或地址)从优先队列被弹出时，它的状态将被设置为淹没状态。因此，每次从优先队列中弹出一个像元时，DCount的值增1，用于记录当前像元弹出时的次序，同时orders中对应当前像元的值设为DCount的值。

(5) 对洼地像元的填充

    对洼地像元的填充过程可以看作是通过两步来实现的：第一步是用当前尚未被海水淹没但又临岸的 (marks标记为2)高程最低像元来假定海平面高度；第二步是更新当前临岸像元8个邻域像元的高程，即当某一邻域像元的高程比当前海平面还低，则需要将其高程提升至海面高度，再加入优先队列中等候被淹没，同时该邻域像元的marks标记值应设为2，以确保所有位于优先队列中像元(均为临岸像元)的marks全部为2。

(6) 输出淹没次序矩阵

由本发明方法计算获得的淹没次序矩阵(orders数组)将作为代替原有DEM的数据集，在其基础上可采用常规的D8方法来完成后续的流向、累积流向、河道汇水网络、Strahler河道分级、流域范围等提取任务。

(8)结果验证

    图3是一个小型的DEM数据实例，经过本发明方法处理之后，得到了洼地填充之后的DEM（图4），同时得到了淹没次序矩阵(即orders数组)(图5)，该矩阵中无重复的编号，避免了原DEM中存在高程十分接近的情形，并且矩阵中也不存在洼地像元，因而可以用于提取河道信息。对于覆盖中国北方及黄河流域的大范围DEM，它含有不少无数据像元区块，经本发明方法处理后得到的淹没次序矩阵，它显示无数据区已被填充，淹没次序矩阵也能够反映原DEM的地形对比状况。将淹没次序矩阵代替原DEM，进行常规的流向提取及水系提取过程，获得了黄河流域水系状况（图6），表明所提取的水系是完整的，且与真实的黄河水系十分接近。