CN102842104B - 面向海量dem数据的高精度河道洪水淹没区生成方法 - Google Patents

Abstract

一种面向海量DEM数据的高精度河道洪水淹没区生成方法，其针对海量DEM数据进行条带划分，并将条带中栅格行上连续多个淹没单元进行游程压缩存储，以降低数据量，最后采用压缩单元边界追踪方法提取淹没范围，从而实现了复杂地形条件下的淹没区生成。本发明将传统淹没源种子点填充方法改为栅格压缩存储和栅格边界追踪计算，使得较小的内存配置处理了较大的栅格数据量，同时避免了种子点填充法中的大量递归判断，提升了计算速度。本发明充分运用了数字地形分析和数据压缩存储技术，避免了种子填充法进行海量DEM数据淹没分析产生的问题，较好的解决了河湖复杂地形环境下的连通域分析，实现大范围、高精度河湖洪水淹没区的生成。

Description

面向海量DEM数据的高精度河道洪水淹没区生成方法

技术领域

本发明涉及水利信息技术领域，具体是一种面向海量DEM数据的高精度河道洪水淹没区生成方法。

背景技术

我国是一个洪水灾害发生十分频繁的国家，每年因洪水导致的损失难以估计。产生洪水灾害有多种原因，平坦区域大面积降雨或河道上游来水下泄都可以导致洪水发生，其中以河道洪水淹没最为常见，这类淹没分析通常称为有源淹没分析。河道洪水淹没分析要解决的一个核心问题是要确定出错综复杂的连通区，尤其是河网水系发育繁茂的流域，这通常也是洪水发生频率较高的地区。

随着空间信息技术的发展，基于数字化地形的洪水淹没预警与灾害评估已成趋势。数字化的洪水淹没分析是以数字高程模型(DigitalElevation Model，DEM)为基础，通过判断DEM中某网格高程值高于或低于设定的水位来确定当前网格是否被淹。若要精确地计算出流域内洪水淹没区，必须针对高分辨率DEM进行分析。

实现DEM数据洪水淹没区生成的一种有效方法是直接扫描整个地形区域，设置或预处理获取与淹没有关的连通域，进而形成淹没区，但该方法操作繁琐，难以实现淹没范围自动提取，对于拓扑结构复杂的河流地形，人工处理工作量较大。另一种方法来源于图像处理中的种子填充法，通过给定的淹没源种子点按照四邻域或八领域进行扩散，较好的模拟了实际的地表径流淹没过程，可以自动提取相关的淹没连通区，较为成熟可靠，这也是目前使用的主流方法。

DEM是以栅格形式存储的地形数据，其数据量与所存储的栅格场网格行列数相关，网格行列数越大则所要存储的数据量越大。当所存储的地形范围越大，或者其地形范围内数据精度越高，则DEM数据越大。以三峡库区淹没处理的规划方案为例，涉及面积约7.9万平方公里，假如以1米分辨率地形进行计算，则DEM数据大小为79GB，这在当前主流个人计算机内存配置条件下难以一次性载人内存。

当前大多数的洪水淹没分析方法都是基于种子点填充方法或其变形，该类方法的主要缺点是：

(1)必须将地形数据一次性载入内存，这将消耗大量内存资源

由于无法预先自动判断洪水淹没范围，传统的基于种子点填充算法的淹没分析必须将DEM地形数据全部载入内存，而以栅格形式存储的地形数据极耗内存资源，以目前个人计算机内存配置条件，难以一次全部读入，严重限制了种子点填充淹没分析方法的实用性；

(2)难以处理海量DEM数据

由于受内存限制较明显，淹没源种子点填充的淹没分析方法难以处理海量DEM数据，因此无法处理大范围的淹没分析，或者无法处理高精度的地形数据，降低了淹没分析方法的适用范围；

(3)即使对DEM数据进行分块处理，也需要多次重复读写数据块，效率低下

受问题1和问题2的约束，一些改进的淹没源种子点填充方法是对海量DEM数据进行分块处理，即将海量DEM数据分块读取，在块内完成淹没源种子点扩散算法，将淹没源种子点扩散结果写入淹没分析结果文件，最后对块间各自淹没区进行连通性判断，其主要缺点为：必须逐块遍历或额外建立块索引表，由于无法知道淹没区范围大小，每个数据块都会遍历到；必须进行繁琐的块间淹没区连通性判断；由于淹没分析结果数据没有保存在内存，多次重复文件读写难以避免，导致计算效率较低。因此，在进行海量DEM数据下的河道淹没分析时该方法仍然不是一个较优的选择。

内陆河道洪水淹没分析必须同时满足高效、大范围、高精度和自动化，而当前的各种淹没分析方法是难以达到这一要求的。考虑到处理海量DEM数据的关键是既要避免大量耗时的文件读写操作，又要能在内存中有效的处理大数据量DEM，因此，本发明提供一种面向海量DEM数据的高精度河道洪水淹没区生成方法，以满足水利行业科研和产业的需求。

发明内容

本发明提供一种面向海量DEM数据的高精度河道洪水淹没区生成方法，可以解决上述技术问题，满足内陆河道洪水淹没分析的高效、大范围、高精度和自动化的要求。

一种面向海量DEM数据的高精度河道洪水淹没区生成方法，包括如下步骤：

第一步、设定淹没源种子点、淹没水位和DEM条带划分行数，计算DEM条带总数和每个条带在DEM中的具体位置，构建一个空的栅格场结构；

第二步、读取淹没源种子点所在条带的DEM数据，对条带中可能被淹没的网格进行压缩存储进栅格场；查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记；

第三步、依条带向上和向下扩展，完成新条带压缩存储和边界追踪；

第四步、若栅格场当前最上条带顶栅格行和最下条带底栅格行中都没有压缩单元被追踪过，或都已到达原DEM数据边界栅格行，则停止追踪，否则返回执行第三步；

第五步、依条带逐栅格行扫描，查找只有一侧被追踪过的压缩单元，以该压缩单元另一侧为追踪起点进行追踪，实现洪水淹没区中孤岛的提取；

第六步、将栅格场上被追踪的压缩单元生成一副淹没范围栅格图。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第一步中设定淹没源种子点和淹没水位时，淹没源种子点坐标位于原DEM数据范围内，且淹没源种子点位置所对应的原DEM数据中的网格单元高程值小于设定的淹没水位。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第一步中DEM条带划分行数远小于原DEM数据的总栅格行数，当划分到最后一个条带时，若原DEM上剩余栅格行数小于设定的条带行数，则将条带上不足的栅格行以空值补齐。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第二步中读取淹没源种子点所在条带DEM数据具体方法为：根据淹没源种子点所在的空间坐标x，y值，以及原DEM的最小最大坐标和网格行列数，计算淹没源种子点位于原DEM所在网格的行列号，然后通过条带参数计算得到该网格所属条带以及位于条带中的实际行列位置，最后读取DEM数据文件中位于此条带上的网格高程值数据。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第二步中对条带中可能被淹没的网格进行压缩存储进栅格场的步骤具体为：

(1)在第一步构建的栅格场上，为每个栅格行构建一个空的压缩单元链表，该链表的结点为压缩单元结构；

(2)在淹没源种子点所在的条带上，按行进行遍历，依次将每个网格的高程值与给定淹没水位进行对比，若网格高程值小于给定淹没水位则表示该网格可能属于淹没区，称之为可能淹没网格；同时，对连续的几个可能淹没网格采用游程压缩编码技术进行压缩，生成一个压缩单元，该压缩单元结构标记了可能淹没网格在此栅格行上跨越的起始列和终止列；

(3)在执行步骤(2)时，每生成一个新的压缩单元，则将该压缩单元插入到栅格场中对应行的压缩单元链表尾部，完成条带数据的栅格压缩存储。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，第二步中查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记的步骤具体为：

首先，为当前条带每个栅格行设置一个标记数组，数组元素的个数是栅格行上压缩单元的两倍，每个压缩单元有左侧和右侧两个追踪参数，默认将所有的标记数组元素设置为false，即表示压缩单元左侧和右侧都未被追踪过；

其次，计算出淹没源种子点位于DEM数据中的网格单元，根据该网格单元的行列号判定其所属的压缩单元，分别以该压缩单元左侧和右侧作为起始追踪边，按照四方向追踪法进行边界追踪，每追踪到某压缩单元的左侧或右侧，就将追踪到的压缩单元的左侧或右侧标记由false设置为true，当追踪到某侧标记为true的压缩单元，则停止追踪；

最后，循环遍历淹没源种子点所在栅格行上所有压缩单元，若某单元一侧被追踪过而另一侧未被追踪，则以未被追踪边作为起始追踪边进行边界追踪，直到该栅格行上所有被追踪过的压缩单元满足左侧和右侧都被追踪过。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第三步的具体步骤为：在栅格场中查找当前最上条带的顶栅格行，从该顶栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最上条带的上方，进行游程压缩存储进栅格场，以提取的被追踪压缩单元所在的追踪边向上进行边界追踪，最后，循环遍历上述顶栅格行中的压缩单元，对仅有一侧被追踪过的压缩单元，对其未追踪的另一侧进行向上边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪；同理，在栅格场中查找当前最下条带的底栅格行，从该底栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最下条带的下方，进行游程压缩存储，以提取的被追踪压缩单元所在的追踪边向下边界追踪，最后，循环遍历底栅格行压缩单元，对仅有一侧被追踪过的压缩单元，对其未追踪的某侧进行向下边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪。如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第五步中淹没区孤岛边界的提取方法为：对栅格场中每个栅格行进行遍历，两两提取追踪标记值，若两个追踪值一个为true，一个为false，则表示该标记所代表的压缩单元有一侧未被追踪，此时以该压缩单元未被追踪的侧为起始追踪边进行边界追踪。

如上所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，所述第六步中生成淹没范围栅格图的方法为：按照栅格场读入的条带总范围设置淹没区栅格图的基本参数，包括最小最大坐标、网格行列值、网格宽度；逐行遍历栅格场中的压缩单元，将未被追踪的压缩单元和压缩单元之间的“空隙”所在的网格值设置为nodata，将被追踪过的压缩单元所在的网格值设置为淹没水位值，最后生成淹没区栅格图。

本发明针对海量DEM数据进行条带划分，并将条带中栅格行上连续多个淹没单元进行游程压缩存储，以降低数据量，最后采用压缩单元边界追踪方法提取淹没范围，从而实现了复杂地形条件下的淹没区生成。本发明将传统淹没源种子点填充方法改为游程压缩存储和栅格边界追踪计算，使得较小的内存配置处理了较大的栅格数据量，同时避免的种子填充法中的大量递归判断，提升了计算速度。

本发明的有益效果：

(1)采用游程压缩编码思想实现海量DEM数据的压缩存储，有效解决了常规种子填充方法难以读取大数据量DEM的问题。考虑到海量DEM数据中每个网格单元只有“淹没”和“未淹没”两种状态，这是一种典型的二值栅格数据，可以在将DEM条带数据读入内存的同时进行二值栅格化，并实现游程编码压缩。该压缩单元数据与全栅格数据相比，有着较大的压缩比，对于高程值大于设定淹没水位的网格则可以直接剔除，不用存储，而栅格行上连续多个高程值小于设定淹没水位的网格单元则可以进行游程压缩存储。

(2)在基于分块DEM的种子填充方法中，由于计算结果没有保存到内存，涉及到淹没区提取和块间淹没区连通性判定时，需要多次重复读取文件，而文件读写是一种较为耗时的操作。本发明有效避免了多次重复的文件读写操作，除读取原DEM和生成淹没区栅格图外，所有的淹没分析计算都在内存中进行。

(3)由于无法判断淹没区范围大小，无论是常规种子填充方法还是分块种子填充方法都需要读取整个DEM数据，当淹没区范围远小于DEM数据的范围时，对整个DEM进行读取显得冗余，效率低下。本发明在读取每个条带时都对淹没区进行边界追踪，当发现最上条带顶栅格行和最下条带底栅格行上没有被淹没的网格时即停止读取新的条带。这一操作有效降低了保留在内存中的数据量，也提升了淹没区生成的速度，避免了淹没区范围的大量无效判断。

(4)本发明对于小范围的DEM数据仍然适用。因为条带数目是通过原DEM行数除以单个条带行数得到，当原DEM行数较小，则条带数目也较少，当原DEM行数小于设定的单个条带行数，则总的条带数目为一个，本发明对于这一变化同样适用。同时，在单个条带下，基于压缩单元的边界追踪方法效率仍然优于常规淹没源种子点填充方法。

附图说明

图1是本发明面向海量DEM数据的高精度河道洪水淹没区生成方法其中一个实施例的流程示意图；

图2是某地区DEM数据示例；

图3是给定某淹没源种子点及其淹没水深下的洪水淹没分析结果图；

图4是淹没区的局部放大图；

图5是DEM数据及其条带划分示意图；

图6是淹没源种子点及其与条带上栅格单元示意图；

图7是读取条带数据示例；

图8是条带上栅格压缩示意图；

图9是压缩栅格边界追踪结果示意图；

图10是淹没源种子点所在栅格行循环遍历追踪示意图；

图11是上一条带读取及其压缩结果示意图；

图12是向上一条带边界追踪结果示意图；

图13是下一条带读取及其压缩结果示意图；

图14是淹没区孤岛边界追踪示意图。

具体实施方式

下面将结合本发明中的附图，对本发明中的技术方案进行清楚、完整地描述。

图1所示为本发明面向海量DEM数据的高精度河道洪水淹没区生成方法其中一个实施例的流程示意图，图2给出了某地区的DEM数据示例，该DEM数据包含的网格行列数位19719行、19454列，数据量大小为1.43GB，网格宽度25米。图3为给定某淹没源种子点及其淹没水深下的洪水淹没分析结果图，图4为局部放大图，可以从图中清晰的看到淹没区所形成的“空白”洞存在，这是由于此处的淹没孤岛所形成。

本发明实施例面向海量DEM数据的高精度洪水淹没区生成方法，如图1所示，其包括如下步骤：

第一步、设定淹没源种子点、淹没水位和DEM条带划分规格参数(即DEM条带划分行数)，计算DEM条带总数和每个条带在DEM中的具体位置，构建一个空的栅格场结构；

DEM数据一般都以符合国际标准的文件格式进行存储，如GeoTiff格式等，这些DEM格式数据文件存储了该DEM基本规格信息，如DEM最大、最小坐标、网格行数、网格列数、网格宽度、无效值等。图5是DEM数据及其条带划分示意图，该DEM有24行、32列，其左下角和右上角分别是DEM所在区域的最小和最大空间坐标，共有24×32＝768个网格单元，且每个网格单元的宽度是固定的。对于DEM栅格场上无值的网格，常以无效值进行标识，如设置为一个负数值，如图2所示DEM数据左上角大片黑色区域代表无值区。

其中，设定淹没源种子点和淹没水位的原则为：该淹没源种子点坐标必须位于原DEM数据范围内，且淹没源种子点位置所对应的原DEM数据中的网格单元高程值必须小于设定的淹没水位，否则计算将终止。

DEM条带划分行数设置方法为：原DEM数据所在网格为M行N列组成，在划分条带时设定单个条带的列数与原DEM数据网格列数相同，而条带的行数远小于原DEM条带的行数，如设置K行，则原DEM数据将会划分为M/K个条带。当划分到最后一个条带时，若原DEM上剩余栅格行数小于设定的条带行数，则将条带上不足的栅格行以空值补齐。因为原DEM标记了其所在矩形网格区域的最大最小坐标、网格行列数、网格宽度等参数，当对DEM数据按照从下往上划分条带后，每个条带对应的最小最大坐标、网格行列数和网格宽度等参数及其在原DEM数据中的位置可以通过上述参数计算得到。在内存构建一个不含数据的空的栅格场结构，该栅格场行列数目、最小最大坐标等参数与原DEM数据规格参数相同。

在进行条带划分时，条带的行数是可以自行设定的，如图5中设定的条带为8行，则对原DEM数据按照从下往上划分，可以分为3个条带。因为每个条带的栅格行数是已知的，因此每个条带的空间坐标范围以及在原DEM上所处的位置是可以计算得到的。

第二步、读取淹没源种子点所在条带的DEM数据，对条带中可能被淹没的网格进行压缩存储进栅格场；查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记；

其中读取淹没源种子点所在条带DEM数据具体方法为：根据淹没源种子点所在的空间坐标x，y值，以及原DEM的最小最大坐标和网格行列数，计算淹没源种子点位于原DEM所在网格的行列号，然后通过条带参数计算得到该网格所属条带以及位于条带中的实际行列位置，最后读取DEM数据文件中位于此条带上的网格高程值数据。

为便于图示说明，在图6中，将DEM网格上的高程值以整数进行表达，以1到8的数字代表高程值，数字越大表示高程值越高。淹没源种子点所对应的原DEM网格行列号可以通过淹没源种子点空间x，y坐标反算得到，图中涂成黑色的网格代表淹没源种子点所在网格，其位于DEM数据的第10行、第14列，位于第2个条带的第2行第14列。

假定设定的淹没水位值为4，此处淹没源种子点所在的网格高程值为3，在将淹没源种子点所在条带的高程值数据读入内存后，可以对其进行压缩存储，对条带中可能被淹没的网格进行压缩存储进栅格场的步骤具体为：

(1)在第一步构建的栅格场上，为每个栅格行构建一个空的压缩单元链表，该链表的结点为压缩单元结构。

(2)在淹没源种子点所在的条带上，按行进行遍历，依次将每个网格的高程值与给定淹没水位进行对比，若网格高程值小于给定淹没水位则表示该网格可能属于淹没区，称之为可能淹没网格。同时，对连续的几个可能淹没网格采用游程压缩编码技术进行压缩，同时生成一个压缩单元，该压缩单元结构标记了可能淹没网格在此栅格行上跨越的起始列和终止列。

(3)在执行步骤(2)时，每生成一个新的压缩单元，则将该压缩单元插入到栅格场中对应行的压缩单元链表尾部。因为栅格行上会有多个压缩单元标记的淹没网格集合，将压缩单元采用链表串联，插入到上述栅格场中对应的栅格行中去，完成条带数据的压缩存储。

本发明将栅格数据采用游程编码技术进行压缩存储，每个压缩单元记录了其两端位于栅格行中的列值，并采用链表的形式组织栅格行上的压缩单元集合，在初始条件下，栅格场上每行压缩单元数为零，因此，每行的压缩单元链表为空。附图8是压缩压缩单元示意图，以图中第9行的压缩单元为例，该压缩单元标记了1和18两个值，表示该压缩单元跨越了第1至第18列。与原DEM上按照行和列存储的二维全栅格数据相比，压缩存储极大的降低了数据量。

在第二步中，查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记的步骤具体为：首先为当前条带每个栅格行设置一个标记数组，数组元素的个数是压缩单元的两倍，每个压缩单元有左侧和右侧两个追踪参数，默认将所有的标记数组元素设置为false，即表示未被追踪过。其次，根据淹没源种子点计算出该淹没源种子点位于DEM数据中的网格单元，根据该网格单元的行列号判定其所属的压缩单元，分别以该压缩单元左侧和右侧作为起始追踪边，按照四方向追踪法进行边界追踪，每追踪到某压缩单元的左侧或右侧，就将追踪到的压缩单元的左侧或右侧标记由false设置为true，当追踪到某侧标记为true的压缩单元，则停止追踪。最后，循环遍历淹没源种子点所在栅格行上所有压缩单元，若某单元一侧被追踪过而另一侧未被追踪，则以未被追踪边作为起始追踪边进行边界追踪，直到该栅格行上所有被追踪过的压缩单元满足左侧和右侧都被追踪过。

下面举例说明：在进行边界追踪时，首先找到淹没源种子点所在的压缩单元，该压缩单元位于第10行上(如图8所示)。以该压缩单元左侧开始追踪，按照目前较为成熟的栅格边界追踪法，采用上、下、左、右四个方向追踪边界线，其边界追踪结果如图9所示。在执行边界追踪后，可以确定某些压缩单元的左侧或右侧边界已经被追踪过，如图中粗线所示。在读取条带并完成栅格单元压缩后，为每个压缩单元左侧和右侧设定标记，该标记初始值都为false，表示压缩单元左侧和右侧都未追踪。在追踪过程中，每当追踪到一个压缩单元某侧时，若该侧标记值为false，即将其设置为true，否则停止追踪。该标记是通过为每个栅格行设置数组完成的，每行标记数组元素的个数是该行压缩单元个数的两倍，且该数组按照从小到大两两查找就可以判断每个压缩单元的左侧和右侧追踪情况。

由于淹没源种子点所在压缩单元两侧进行追踪后，淹没源种子点所在栅格行有部分压缩单元会出现一侧被追踪，但另一侧未被追踪的情形，而未被追踪的侧向也是淹没区的边界之一，需要即时追踪出来。因此，可以循环遍历淹没源种子点所在栅格行上所有压缩单元，只要某压缩单元两侧追踪标记值不同，则将压缩单元标记值为false的一侧作为起始追踪边进行边界追踪，直到不存在满足上述情形的压缩单元为止。如图10所示，淹没区外边界为C1，内边界分别为C2、C3和C4，淹没源种子点位于条带X，淹没源种子点所在栅格行上有压缩单元A、B、C、D，追踪标记数组元素为1-8。如图10所示，假定淹没源种子点所在的压缩单元为B，则对压缩单元B两侧执行边界追踪，此时标记2、3、4、5所对应的压缩单元左侧或右侧都被追踪过，执行第一次循环，按照A、B、C、D压缩单元两两提取追踪标记值，发现压缩单元A和压缩单元C的两侧标记值不相等，此时做压缩单元A的左侧和压缩单元C的右侧执行边界追踪，此时标记1、6、7、8都被设置成true，其中标记1、8由压缩单元A左侧追踪设定，标记6、7由压缩单元C右侧追踪设定。执行第二次循环时，发现没有出现两侧标记值不相等的压缩单元，此时跳出循环。

第三步、依条带向上和向下扩展，完成新条带压缩存储和边界追踪，其主要步骤为：

在栅格场中查找当前最上条带的顶栅格行，从该顶栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最上条带的上方，进行游程压缩存储进栅格场，以提取的被追踪压缩单元所在的追踪边向上边界追踪，最后，循环遍历顶栅格行压缩单元，对仅有一侧被追踪过的压缩单元，对其未追踪的某侧进行向上边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪；同理，在栅格场中查找当前最下条带的底栅格行，从该底栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最下条带的下方，进行游程压缩存储，以提取的被追踪压缩单元所在的追踪边向下边界追踪，最后，循环遍历底栅格行压缩单元，对仅有一侧被追踪过的压缩单元，对其未追踪的某侧进行向下边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪。

因为原DEM数据是按照条带读取的，条带数据在进行游程压缩后一次性存储到栅格场中对应的各个栅格行中，可以预先知道当前已经有哪些条带数据存储到栅格场中。为便于描述，按照每个条带所处的空间位置，将位于栅格场中最上方的条带称为最上条带，将位于最下方的条带称为最下条带，将每个条带中最上一栅格行称为该条带的顶栅格行，最下一栅格行称为该条带的底栅格行。特别的，在完成步骤二中种子点所在条带的压缩存储后，栅格场上最上条带和最下条带代表同一个条带。

如图9所示，栅格场上只存储有一个条带的数据(即淹没源种子点所在的条带)，其顶栅格行(第16行)有两个压缩单元，且两个压缩单元的左侧和右侧都被追踪过，因此需要从原DEM中读取当前条带的上一条带，并按照第二步中的压缩步骤执行游程压缩存储。图11显示了读取上一条带后的栅格场变化，在栅格场的第17至第24行上新增加了压缩单元。此时分别以第16行被追踪过边界的两个压缩单元的左侧和右侧为起始追踪边向上进行边界追踪，最后，按照第二步中的方法循环遍历顶栅格行(16行)压缩单元，当没有出现压缩单元一侧被追踪而一侧未被追踪的情况时停止遍历。其边界追踪结果如图12所示，图中压缩单元边界以粗线条显示的为追踪边界，此时淹没区已经追踪出了一个孤岛，该孤岛横跨16-18栅格行。

图9中，当前条带底栅格行位于第9行，有一个压缩单元，且两侧已经被追踪过。此时需要读取原DEM中对应的下一条带，即图6所示的条带1。在读取条带1后，进行游程压缩存储，其压缩后结果如图13所示。此时分别以第9行所在压缩单元左侧和右侧作为追踪边，执行向下边界追踪，由于第10行上压缩单元与第9行压缩单元并不相连，追踪停止。按照第二步所述的方法循环遍历底栅格行(9行)压缩单元，没有出现压缩单元一侧被追踪一侧未被追踪的情况，停止遍历。

第四步、若栅格场当前最上条带顶栅格行和最下条带底栅格行中都没有压缩单元被追踪过，或都已到达原DEM数据边界栅格行(表示淹没区范围止于该条带，不需要再向上或向下读取新条带)，则停止追踪，否则返回执行第三步；

如图13所示，由于当前最上条带顶栅格行位于第24行，最下条带底栅格行位于第1行，都位于原DEM栅格行边界，因此停止读取新的条带。

第五步、依条带逐栅格行扫描，查找只有一侧被追踪过的压缩单元，以该压缩单元另一侧为追踪起点进行追踪，实现洪水淹没区中孤岛的提取。淹没区孤岛边界的提取方法为：对栅格场中每个栅格行进行遍历，两两提取追踪标记值，若两个追踪值一个为true，一个为false，则表示该标记所代表的压缩单元有一侧未被追踪，此时以该压缩单元未被追踪的侧(即标记为false的侧)为起始追踪点进行边界追踪。

如图13所示，原淹没区中仍然有两个孤岛未被追踪出来，其中一个孤岛横跨栅格场中第12至15行，另一个孤岛横跨第19至22行。因为每个栅格行都存储了压缩单元左侧和右侧追踪标记，可以遍历每个行上的追踪标记数组，并逐个两两取值，该值即代表对应压缩单元的左侧和右侧追踪情况。若两个标记值都为false，则表示该压缩单元与淹没区不连通，不属于淹没区；若两个标记值都为true，则表示该压缩单元左侧和右侧都已经追踪过，不需处理。仅需要处理两个标记分别为true和false的情况。如图13中第12行第一个压缩单元，其左侧已经被追踪，但是右侧未被追踪过，此时以该压缩单元右侧为起始追踪边进行边界追踪，将孤岛边界追踪出来。按照此方法，最后追踪结果如图14所示，此时整个淹没区外侧边界和内侧边界已经全部追踪出来。

第六步、将栅格场上被追踪的压缩单元生成一副淹没范围栅格图。生成淹没范围栅格图的方法为：按照栅格场读入的条带总范围设置淹没区栅格图的基本参数，包括最小最大坐标、网格行列值、网格宽度等；逐行遍历栅格场中的压缩单元，将未被追踪的压缩单元和压缩单元之间的“空隙”所在的网格值设置为nodata，将被追踪过的压缩单元所在的网格值设置为淹没水位值，最后生成淹没区栅格图。特别的，nodata是DEM栅格数据标记的一种通用做法，即将网格中的无效值以一个无意义的数值进行标记。

一般情形下，淹没区范围所跨的栅格行数小于原DEM所跨的栅格行数，如图14所示，该淹没区范围所跨的栅格行从第9行至第24行。因此，可以按照淹没区的实际范围进行输出，此时淹没区的行数为16行，最小坐标与条带2的最小坐标相同，最大坐标与条带3的最大坐标相同。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何属于本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (6)

1.一种面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于包括如下步骤：

第一步、设定淹没源种子点、淹没水位和DEM条带划分行数，计算DEM条带总数和每个条带在DEM中的具体位置，构建一个不含数据的空的栅格场结构，所述栅格场结构用来标记栅格场的最小坐标、最大坐标、栅格行数、栅格列数；

第二步、读取淹没源种子点所在条带的DEM数据，对条带中可能被淹没的网格进行压缩存储进栅格场；查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记；其中读取淹没源种子点所在条带DEM数据具体方法为：根据淹没源种子点所在的空间坐标x，y值，以及原DEM的最小最大坐标和网格行列数，计算淹没源种子点位于原DEM所在网格的行列号，然后通过条带参数计算得到该网格所属条带以及位于条带中的实际行列位置，最后读取DEM数据文件中位于此条带上的网格高程值数据；

第三步、依条带向上和向下扩展，完成新条带压缩存储和边界追踪；

第四步、若栅格场当前最上条带顶栅格行和最下条带底栅格行中都没有压缩单元被追踪过，或都已到达原DEM数据边界栅格行，则停止追踪，否则返回执行第三步；

第五步、依条带逐栅格行扫描，查找只有一侧被追踪过的压缩单元，以该压缩单元另一侧为追踪起点进行追踪，实现洪水淹没区中孤岛的提取；

第六步、将栅格场上被追踪的压缩单元生成一幅淹没范围栅格图；

所述第二步中对条带中可能被淹没的网格进行压缩存储进栅格场的步骤具体为：

(1)在第一步构建的栅格场上，为每个栅格行构建一个空的压缩单元链表，该链表的结点为压缩单元结构；

(2)在淹没源种子点所在的条带上，按行进行遍历，依次将每个网格的高程值与给定淹没水位进行对比，若网格高程值小于给定淹没水位则表示该网格可能属于淹没区，称之为可能被淹没的网格；同时，对连续的几个可能被淹没的网格采用游程压缩编码技术进行压缩，生成一个压缩单元，该压缩单元结构标记了可能被淹没的网格在此栅格行上跨越的起始列和终止列；

(3)在执行步骤(2)时，每生成一个新的压缩单元，则将该压缩单元插入到栅格场中对应栅格行压缩单元链表尾部，完成条带数据的栅格压缩存储。

2.如权利要求1所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于：所述第一步中设定淹没源种子点和淹没水位时，淹没源种子点坐标位于原DEM数据坐标范围内，且淹没源种子点位置所对应的原DEM数据中的网格单元高程值小于设定的淹没水位。

3.如权利要求1所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于：第二步中查找淹没源种子点所在的压缩单元，进行淹没区边界追踪，同时对追踪过的压缩单元进行标记的步骤具体为：

首先，为当前条带每个栅格行设置一个标记数组，数组元素的个数是栅格行上压缩单元的两倍，每个压缩单元有左侧和右侧两个追踪参数，默认将所有的标记数组元素设置为false，即表示压缩单元左侧和右侧都未被追踪过；

其次，计算出淹没源种子点位于DEM数据中的网格单元，根据该网格单元的行列号查找其所属的压缩单元，分别以该压缩单元左侧和右侧作为起始追踪边，按照四方向追踪法进行边界追踪，每追踪到某压缩单元的左侧或右侧，就将追踪到的压缩单元的左侧或右侧标记由false设置为true，当追踪到某侧标记为true的压缩单元，则停止追踪；

最后，循环遍历淹没源种子点所在栅格行上所有压缩单元，若某单元一侧被追踪过而另一侧未被追踪，则以未被追踪边作为起始追踪边进行边界追踪，直到该栅格行上所有被追踪过的压缩单元满足左侧和右侧都被追踪过。

4.如权利要求3所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于：所述第三步依条带向上和向下扩展，完成新条带压缩存储和边界追踪的具体步骤为：在栅格场中查找当前最上条带的顶栅格行，从该顶栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最上条带的上方，进行游程压缩存储进栅格场，以提取的被追踪压缩单元所在的追踪边向上进行边界追踪，最后，循环遍历上述顶栅格行中的压缩单元，对仅有一侧被追踪过的压缩单元，将其未追踪的另一侧进行向上边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪；同理，在栅格场中查找当前最下条带的底栅格行，从该底栅格行中提取所有被追踪过的压缩单元，从原DEM中读取新的条带，该条带位于当前最下条带的下方，进行游程压缩存储，以提取的被追踪压缩单元所在的追踪边向下边界追踪，最后，循环遍历底栅格行压缩单元，对仅有一侧被追踪过的压缩单元，对其未追踪的某侧进行向下边界追踪，直到所有单侧追踪压缩单元都已经完成了双侧追踪，该过程将完成新读入条带中压缩单元的边界追踪。

5.如权利要求4所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于：所述第五步中淹没区孤岛边界的提取方法为：对栅格场中每个栅格行进行遍历，两两提取追踪标记值，若两个追踪值一个为true，一个为false，则表示该标记所代表的压缩单元有一侧未被追踪，此时以该压缩单元未被追踪的侧为起始追踪边进行边界追踪。

6.如权利要求1所述的面向海量DEM数据的高精度河道洪水淹没区生成方法，其特征在于：所述第六步中生成淹没范围栅格图的方法为：按照栅格场读入的条带总范围设置淹没区栅格图的基本参数，包括最小最大坐标、网格行列值、网格宽度；逐行遍历栅格场中的压缩单元，将未被追踪的压缩单元和压缩单元之间的″空隙″所在的网格值设置为nodata，将被追踪过的压缩单元所在的网格值设置为淹没水位值，最后生成淹没区栅格图。