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

一种基于栅格数字高程模型(DEM)的流域洼地填充预处理方法，涉及地理信息系统及分布式水文模型中基于栅格DEM的流域数字水系提取方法，尤其是能够应用于大河流域的方法。本发明专利技术利用一个以高程为关键字的优先队列，从DEM边缘最低像元处开始向内部逐个淹没像元，记录各像元的淹没次序，得到淹没次序矩阵。所有与海面相邻的像元加入优先队列，每次由优先队列弹出高程最低的像元，海面高度升至该像元高度，记录该像元的淹没次序，然后将该像元的邻域像元中高程低于海面的所有像元升至海面高度，并加入到优先队列中。由淹没次序矩阵代替原有DEM数据集完成流向提取和河道网络提取等任务，可提取出大河流域的完整水系网络。

【技术实现步骤摘要】
基于栅格DEM的大河流域洼地填充预处理方法
本专利技术涉及地理信息系统以及分布式水文模型中基于栅格数字高程模型的流域数字水系提取方法，尤其是能够应用于大河流域的洼地填充预处理方法。
技术介绍
基于数字高程模型(DEM)提取流域水系信息是实现分布式水文模型数值模拟的重要方法。流域水系信息的提取包括计算流向、累积流向、水系分级以及汇水网络等。目前常用八向方法(简写为D8方法)来获得上述信息，然而这些信息能否成功提取还取决于对洼地及平地像元的预处理。目前已有许多可用于这些任务的软件工具，如著名GIS商业软件ArcGIS中的空间分析工具，ENVI的插件RiverTools，GlobalMapper软件等，这些软件对于较小的规模或者地理空间范围不大的流域可以成功提取出完整的水系。由于大河流域涉及范围较广，数据量较大，DEM中会存在大量高程误差，以及高程完全相同的“平地”像元，即使经过洼地填充处理的DEM仍难以正确提取出完整的水系。事实证明，对于黄河中上游流域，采用上述软件工具仍无法正确提取，其输出结果是一堆“断头河”。采用Plachon&Doubarx(2002)[Planchon,O.andDarboux,F.,2002.Afast,simpleandversatilealgorithmtofillthedepressionsofdigitalelevationmodels.CATENA,46(2-3):159-176.]方法及其各种变形，也得到的是类似的“断头河”。Liu[LiuYonghe,2009.AnotherFastandSimpleDEMDepression-FillingAlgorithmBasedonPriorityQueuetructure.AtmosphericandOceanicScienceLetters,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是该专利技术方法的执行流程。图2marks数组示例(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的值本文档来自技高网...

【技术保护点】
一种基于栅格DEM的大河流域洼地填充预处理方法，其特征在于，包括以下步骤：(1)创建并初始化一个以浮点数为关键字的优先队列对象；(2)将栅格DEM四条外边界上的像元按照高程从低到高的顺序加入优先队列；(3)当优先队列不空时，继续执行第（4）步；否则从第（12）步开始执行；(4)弹出排在最前面的像元，它是优先队列中存放的所有像元中高程最低的，把它作为当前像元，它的高程作为当前海平面高度；(5)借助计数器为当前像元设置淹没次序编号，并将此编号存入淹没次序矩阵中的相应位置上；(6)若当前像元已是淹没状态，则返回第3步执行，否则将其标记为淹没状态后执行第7步；(7)找出当前像元周围8个方向上的所有未被淹没的邻域像元；(8)将低于海平面的邻域像元执行步骤9‑10；高于海平面的邻域像元直接加入优先队列，之后返回至第3步；(9)将尚未淹没的邻域像元加入优先队列，并将其标记为临岸状态；(10)当前邻域高程改设为当前海平面的高度；返回第8步将其它邻域像元完成同样的操作；(11)返回至第（3）步；(12)输出淹没次序矩阵。

【技术特征摘要】
1.一种基于栅格DEM的大河流域洼地填充预处理方法，其特征在于，包括以下步骤：(1)创建并初始化一个以浮点数为关键字的优先队列对象；(2).将栅格DEM四条外边界上的像元按照高程从低到高的顺序加入优先队列；(3).当优先队列不空时，继续执行第（4）步；否则从第（12）步开始执行；(4)弹出排在最前面的像元，它是优先队列中存放的所有像元中高程最低的像元，把它作为当前像元，它的高程作为当前海平面高度；(5).借助计数器为当前像元设置淹没次序编号，并将此编号存入淹没次序矩阵中的相应位置上；(6).若当前像元已是淹没状态，则返回第（3）步执行，否则将其标记为淹没状态后执行第(7)步；(7)找出当前像元周围8个方向上的所有未被淹没的邻域像元；...

【专利技术属性】
技术研发人员：刘永和，李艳利，胡永红，王燕平，李艳粉，
申请(专利权)人：河南理工大学，
类型：发明
国别省市：河南;41