洪水淹没分析是进行水文预测预报、洪水灾害评估的一项重要内容，也是GIS地形三维仿真系统中的一个重要功能^{[1, 2, 3]}。洪水淹没分析主要基于DEM数据，结合水文水动力学构建的洪水演进模型或根据设定的静态“水位面”进行计算。前者可以较精确地模拟洪水在不同时间段下的淹没区变化过程，但是洪水演进模型的构建需要较多的水文水力参数条件，其处理较为复杂^{[4, 5, 6, 7, 8]}。基于静水位下的洪水淹没分析是对洪水淹没过程达到最终平衡状态的近似模拟，尽管没有洪水演进模型方法精确，但其计算简单、实用，能够快速地确定洪水淹没区和水深分布，得到了广泛的应用。静水位下洪水淹没分析主要分为无源淹没和有源淹没两种，其中无源淹没相对简单，只需遍历全部网格即可得出淹没范围，而有源淹没需要考虑如环形山、堤坝等障碍，其计算相对复杂，是洪水淹没分析算法研究的重点。现有洪水淹没分析算法的研究主要来源于图像处理中的种子填充法及其变形，通过给定淹没起算点，在DEM网格中按照四邻域或八邻域进行扩散，可以较好地模拟实际的地表径流淹没过程，以及自动提取淹没区。如文献^[9，10]提出的基于GIS的复杂地形淹没分析及其灾害评估；文献^[11]提出的基于GIS格网模型的洪水淹没分析；文献^[12]提出的基于DEM的“环形”洪水淹没算法；文献^[13]提出的基于DEM的洪水淹没分析等。

由于无法预先判断洪水淹没的大致范围，常规种子点填充算法必须将DEM数据一次性载入内存进行计算，其所能载入的数据量依赖于计算机内存配置。一般三维仿真系统中用于地形表达的DEM数据常常达到几GB甚至几十GB，这在目前个人计算机内存配置条件下是难以全部读入内存的，当进行洪水淹没分析时，容易受到内存的限制导致计算失败。另一方面，当DEM数据较大，淹没范围较广时，种子填充算法中大量的递归操作常导致算法效率急剧降低，同时容易出现堆栈溢出等问题，这极大地限制了种子点填充淹没分析算法的应用范围。因此，在进行大规模DEM数据的洪水淹没分析时，如何避免大量耗时的文件读写操作，同时又能在内存中高效地处理大数据量DEM，是实现高效、大范围、高精度和自动化洪水淹没分析的关键。

针对上述问题，研究人员亦提出了一些优化算法，如文献^[11]提出将DEM数据转换为不规则多边形格网和三角网，在多边形和三角形格网的基础上按照种子点扩散方法进行淹没范围计算。^[14]提出为DEM数据构建瓦片金字塔，在大区域场景下用上层金字塔数据进行淹没分析，以实现快速三维实时显示。这类方法通过对DEM数据进行简化，以简化数据实现快速种子点填充计算，在一定程度上降低了算法对内存的依赖，提升了计算速度，但是淹没区范围计算精度会直接受到简化数据影响，同时针对海量DEM数据的处理仍然有限。

一种较为有效的策略是对DEM数据进行条带分块读取，并实时计算淹没区。以分块种子点填充算法为例，其基本思路为：首先对原DEM数据进行条带划分，在每个条带上采用种子点填充算法计算所有潜在淹没区并编号，然后对相邻条带间潜在淹没区执行连通域检测，查找属于同一连通域的潜在淹没区并构建关联表，最后通过初始种子点所在淹没区编号和关联表提取整个淹没区。如图 1所示，DEM数据被划分为3个条带，图中灰色填充网格代表各条带上的潜在淹没区。通过对比相邻条带间对应网格单元所在的淹没区编号，执行连通域检测并构建关联表，如条带2上的1、2号淹没区和条带3上3号淹没区属于同一个连通域，它们是相互关联的。最后，根据种子点所对应的1号淹没区，逐个条带查找关联表信息并提取其他关联淹没区，构建洪水淹没范围。

图 1 分块种子点填充法示意图 Fig. 1 Strip-divide seed filling algorithm

图选项

上述连通域标记方法在执行条带种子点填充算法后可以及时清空分块数据，降低了对内存的依赖，使得处理海量栅格数据成为可能，在遥感及图像处理领域得到了一定的应用^{[15,16,17,18]}。但将其用于淹没分析时，主要弊端在于：①为进行连通域判定，必须在条带中执行种子点填充算法提取所有潜在淹没区，而大部分潜在淹没区都不是实际淹没范围，这一处理过程较为低效；②条带上种子填充计算结果无法保存在内存中，因此需要预先生成一个与原DEM大小相同的空值栅格文件，及时写入分块条带的计算结果，而在进行连通域判定时则要多次读取此淹没结果文件；③当淹没区范围远小于原DEM范围时，大量条带数据的读取、计算和连通域判定是没有必要的。

考虑到条带中每个网格单元只有潜在淹没和未淹没两种状态，若在读取条带数据后直接将每行中高程值低于设定水位的相邻网格单元按照游程编码方法进行压缩，则可以较大程度上降低栅格数据量，进而将潜在淹没区存储进内存。为避免盲目读取所有条带数据，可以先读取种子点所在条带数据，实现压缩存储并以种子点所在压缩单元进行边界追踪标记，然后向上和向下扩展读取新条带并压缩，每读取一个新条带都以其相邻条带顶栅格行或底栅格行中被追踪过的压缩单元两侧作为追踪边，在新条带上执行边界追踪，最后提取栅格场上被追踪过的所有压缩单元即可构建淹没区，可以称之为分块压缩追踪法。

如图 2所示，首先读取条带2并压缩存储，以种子点所在压缩单元两侧追踪得到淹没区1，该条带底栅格行和顶栅格行都有被追踪过的压缩单元，因此分别读取条带1和条带3并进行压缩存储，分别以条带2底栅格行和顶栅格行上被追踪过的压缩单元两侧作为追踪边向下和向上追踪得到条带1的淹没区1和条带3的淹没区3并进而追踪得到条带2上的淹没区2。因条带1的底栅格行和条带3的顶栅格行都没有被追踪过的压缩单元，停止读取新条带，最后提取所有被追踪过的压缩单元生成洪水淹没范围。

与分块种子点填充法相比，本算法的主要优势在于：①对分块条带直接进行栅格压缩存储，避免了效率低下的种子点填充计算，而将潜在淹没区存储于内存比存储于文件效率更高； ②条带上潜在淹没区以压缩单元形式存储于内存，针对某压缩单元的边界追踪可以跨越多个条带，从而规避了复杂的连通域标记过程；③按照向上和向下扩展读取新条带并及时追踪压缩单元，读取的条带数目与实际淹没范围直接相关，没有洪水淹没的条带则不会读入内存，避免了大量无效条带数据的读取与计算。

图 2 分块压缩追踪法示意图 Fig. 2 Strip-divide compression tracing algorithm

图选项

按照上述算法设计思路，分块压缩追踪法的总体算法流程图如图 3所示。

图 3 算法总体流程 Fig. 3 Algorithm flow chart

图选项

条带压缩与边界追踪标记包含了算法流程中的步骤(1)—步骤(3)。DEM数据条带划分示意图如图 4所示，假定某DEM数据有24行、32列，设定条带规格为8行，则DEM数据被划分为3个条带。为便于说明，将网格上的高程值以整数进行表达，以1—8的数字代表高程值，数字越大表示高程值越高。

图 4 栅格条带分块示意图 Fig. 4 Raster strip-divide diagram

图选项

种子点所在网格位于DEM数据的第10行、第14列，属于条带2的第2行第14列，其高程值为3。假定淹没水位为4，图中以灰色填充的网格为种子点所在条带的潜在淹没范围，其高程值均小于或等于设定水位值。

条带数据的压缩方法为：

(1) 在内存中生成一个以链表为元素的数组，数组大小为原DEM栅格行数目，链表的结点为压缩单元，该链表结构用于条带数据压缩后存储对应栅格行的所有压缩单元。

(2) 在条带上按行遍历，依次将栅格行上高程值小于或等于设定淹没水位的相邻网格采用游程编码方法进行压缩，生成一个压缩单元结构，该压缩单元标记了潜在淹没网格在此栅格行上跨越的起始列和终止列。

(3) 每生成一个新的压缩单元，则将该压缩单元插入对应行上压缩单元链表尾部，最后完成整个条带数据的压缩存储。图 5是种子点所在条带压缩存储示意，以图中第9行的压缩单元为例，该压缩单元标记了1和18两个值，表示该压缩单元跨越了第1至第18列。

图 5 条带压缩与边界追踪 Fig. 5 Strip compression and boundary tracing

图选项

条带数据的压缩存储尚未标记每个压缩单元与淹没源种子点的连通性，可以通过对种子点所在的压缩单元进行边界追踪以标记与种子点连通的其他压缩单元。主要步骤为：首先为每个链表设置一个标记数组，数组元素的个数是链表中压缩单元的两倍，数组依次存储了每个压缩单元的左侧和右侧两个标记，默认将所有的标记数组元素设置为false，即表示所有压缩单元两侧都未被追踪过。其次，查找初始种子点所在的压缩单元，分别以该压缩单元左侧和右侧边界作为起始追踪边，按照四方向追踪法^{[19, 20]}进行栅格边界追踪，每追踪到新的压缩单元左侧或右侧，就将该行上对应的标记数组元素值由false设置为true，当追踪到某侧标记为true的压缩单元则停止追踪。最后，循环遍历初始种子点所在栅格行上所有压缩单元，若某压缩单元一侧被追踪过而另一侧未被追踪，则以未被追踪一侧作为起始追踪边进行边界追踪，直到该行上所有被追踪过的压缩单元满足两侧都被追踪。图 5显示了压缩单元边界追踪结果，追踪边界如图中粗线所示。

依条带扩展存储包含了算法流程中的步骤(5)—步骤(12)，此步骤是一个循环过程，其循环结束条件为步骤(13)。因为可以预先知道当前已读取过哪些条带数据，为便于描述，按照每个条带所处的空间位置，将读取的位于最上方的条带称为最上条带，将位于最下方的条带称为最下条带，将每个条带中最上一栅格行称为该条带的顶栅格行，最下一栅格行称为该条带的底栅格行。在完成种子点所在条带的压缩存储后，可以将该条带同时设定为最上条带和最下条带。

此时可以向上和向下扩展读取新的条带，并按照前述方法进行条带压缩与边界追踪。其原理为：对最上条带和最下条带而言，需要继续向上和向下读取新条带的唯一条件就是其顶栅格行或底栅格行上有被追踪过的压缩单元，表示淹没区还未到达实际边界。当最上条带顶栅格行上有被追踪过的压缩单元时，可以将该行所有被追踪过压缩单元的两侧标记为起始追踪边，读取最上条带的上一条带作为新的最上条带，并进行压缩存储，此时以标记的起始追踪边向上追踪即可得到当前最上条带的淹没区，同理对最下条带执行此过程。不断循环执行最上条带和最下条带的扩展存储，当最上条带顶栅格行和最下条带底栅格行都没有被追踪的压缩单元，或都已到达DEM边界则停止循环。

如图 5所示，链表数组中只存储了一个条带的压缩数据，其顶栅格行(第16行)有两个压缩单元且都被追踪过，因此需要从原DEM中读取当前条带的上一条带，并执行游程压缩存储。图 6显示了读取上一条带后的链表数组变化，在第17至第24行上新增加了压缩单元。此时分别以第16行被追踪过边界的两个压缩单元的左侧和右侧为起始追踪边向上进行边界追踪，图中以粗线条显示的为边界追踪结果。当前最下条带底栅格行位于第9行，有一个压缩单元，且两侧已经被追踪过，需要读取DEM中对应的下一条带并进行压缩存储，结果如图 7所示。此时以第9行压缩单元两侧作为追踪边，执行向下追踪，由于第10行上压缩单元与第9行压缩单元并不相连，追踪停止。

图 6 上一条带压缩存储及边界追踪 Fig. 6 Upper strip compression and boundary tracing

图选项

图 7 下一条带压缩存储及边界追踪 Fig. 7 Lower strip compression and boundary tracing

图选项

淹没区孤岛处理对应算法流程图中步骤(14)。在完成扩展条带压缩存储及边界追踪标记后，需要追踪部分未被标记的孤岛。如图 7所示，淹没区中仍然有两个孤岛未被追踪出来，其中一个孤岛横跨第12至第15行，另一个孤岛横跨第19至第22行。可以逐行遍历追踪标记数组，并两两取值对比，该值代表对应压缩单元的左侧和右侧追踪情况。若两个标记值分别为true和false，则需对标记值为false的一侧进行追踪。如图 7中第12行第一个压缩单元，其左侧已经被追踪，但是右侧未被追踪过，此时以该侧为起始追踪边进行边界追踪标记，将孤岛边界追踪出来。按照此方法，最后追踪结果如图 8所示，此时整个淹没区外侧边界和内侧边界已经全部追踪出来。

图 8 淹没区孤岛提取 Fig. 8 Island extracting in flood region

图选项

基于上述算法设计原理，笔者采用C#语言实现了洪水淹没区生成算法，并在“金沙江水文泥沙信息三维可视化系统”中得到应用。图 9是金沙江流域研究区DEM数据，该数据包含的网格行列数为19719行×19454列，数据文件大小为1.43GB，网格宽度25m，图 10显示了淹没源水位上升50米时的洪水淹没区范围。

图 9 研究区DEM Fig. 9 Study area of DEM data

图选项

图 10 淹没范围示意图 Fig. 10 Flood region diagram

图选项

为验证算法效率，笔者分别采用常规种子填充算法、分块种子点填充算法和本文提出的分块压缩追踪法进行了测试，其中常规种子填充算法因为内存不足导致计算失败。表 1给出了分块压缩追踪法与分块种子点填充法的算法耗时对比。测试结果表明，在相同的条带规格下，分块压缩追踪法耗时远小于分块种子填充法耗时。试验中设定淹没水深为定值，将单个条带的行数从50逐步扩展到2000，当一次读入的条带范围变大时，由于减少了分块数量，同时降低了分块间连通性判定的次数，因此分块种子点填充法耗时是随着分块条带数据的逐步变大而降低的，其最佳算法运行条件为分块条带能刚好被计算机一次性读入内存。分块压缩追踪法是依条带进行扩展读取的，其读取的条带总体范围与实际淹没范围直接相关，尽管在条带分块更小时读取的条带数目更多，但是其实际读取的数据范围可能更小，即文件I/O耗时更少，在条带分块较大时这种差别更明显。而压缩单元边界追踪在不同分块情形下的差别并不大，因为边界追踪需要搜索的压缩单元数目是一致的。随着条带行数的逐步增大，算法总耗时增加并不显著，相反在条带行数较少时算法耗时更少，条带行数从50扩大到2000时算法耗时仅增加了26s，而读取每个分块条带数据所占据的内存则从8MB(50行×2万列×8字节)增加到320MB(2000行×2万列×8字节)。在条带分块策略方面，只要分块数据能够一次性读入内存都是可以满足算法要求的，但分块条带更小时可减少文件I/O耗时，且在一定程度上降低了内存峰值，是较优的选择。

在内存占用方面，常规种子点填充算法需要读取整个DEM数据，以本算法测试数据为例，其需要占用内存约3.2GB(2万行×2万列×8字节)的数据量，这在目前个人计算机上是难以做到的。分块种子点填充法与分块压缩追踪法较好地解决了这一问题，在读取条带数据后通过计算潜在淹没区，可以及时清空分块条带数据所占用的内存，其内存占用的峰值主要来自于条带数据的读入，如条带行数为1000时占用的内存也仅160MB。分块压缩追踪法将潜在淹没区都存储在内存中，测试中整个栅格场共生成了258557个压缩单元，按照每个压缩单元存储起始列和终止列两个字段计算，需要8字节的数据量，则总数据量仅为2MB。若按此栅格数据压缩比进行计算，当栅格数据量达到3.2TB(20万行×200万列×8字节)时，需存储的压缩单元数据量约为2GB。随着64位操作系统的普及，个人计算机内存配置将逐步过渡到4~8GB，本算法将具备处理TB级数据的应用潜力。

本文针对常规种子点填充算法难以用于海量DEM数据淹没分析问题，提出采用分块压缩追踪法进行淹没区生成。通过将DEM数据进行条带划分，并将条带中栅格行上连续多个淹没单元进行压缩存储，以降低数据量，最后采用压缩单元边界追踪标记方法提取淹没范围，从而实现了复杂地形条件下的淹没区生成。与常见种子点填充方法相比，本文方法使用较小的内存配置处理了较大的栅格数据量，同时避免了种子填充法中的大量递归判断，提升了计算速度，具有一定的实用性。