CN102915227A - 面向大区域流域提取的并行方法 - Google Patents

Abstract

本发明公开了一种面向大区域海量DEM数据流域提取的并行方法，其步骤如下：第一步、最优划分粒度评价；第二步、按照数据划分、融合策略进行洼地填平并行计算；第三步、按照数据划分、融合策略，基于洼地填平结果进行水流方向并行计算；第四步、按照数据划分、融合策略,基于水流方向数据，进行汇流累积并行计算；第五步、按照数据划分、融合策略,基于水流方向和汇流累积数据，设置汇流阈值，进行河网水系并行计算；第六步、按照数据划分、融合策略,基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。本发明方法能够在并行时充分考虑数据划分的粒度及并行I/O机制，在串行算法并行解析上能够顾及数据自身特征。

Description

面向大区域流域提取的并行方法

技术领域

本发明涉及一种数字高程模型（DEM）的流域提取并行技术，具体说是一种基于面向大区域海量DEM数据的数据划分、融合策略及串行算法并行解析机制的流域提取并行方法。

背景技术

数字高程模型产品（Digital Elevation Models，缩写为DEM）是地理信息系统（Geographical Information System,缩写为GIS）地理数据库中最为重要的空间信息资料和赖以进行地形分析的核心数据系统，其应用遍布测绘、交通、军事、水利、农业、环境、资源管理、规划与旅游等众多领域，如在测绘部门主要用于三维地形建模、可视域分析、工程土方估算、遥感影像的几何矫正等；在地学界主要用于水文分析、蓄洪计算、淹没分析等三维地学分析；在国土资源部门主要用于土地利用类型详察等。

DEM并不是地形的简单表达，真正价值在于DEM所蕴含的各种地学信息。如何进行“数据挖掘”，获取各种反映实际地形地貌的结构信息，是GIS领域一个重要的研究方向。在此需求环境下，数字地学分析技术应运而生。数字地形分析（Digital Terrain Analysis,DTA）是在DEM上进行地形属性计算和特征提取的数字信息处理技术。在复杂的现实世界地理过程中各影响因子和简单、高效、精确和易于理解的抽象和计算机实现中找到平衡，是数字地形分析的核心任务。用来描述地形特征和空间分布的地形参数有很多，根据地形属性的地学应用范畴，地形属性可分为一般地形属性（如高程、坡度坡向等）和水文特征（如地形结构识别、流域提取等）。其中，流域提取是DTA的重要研究内容。流域是指地表水及地下水分水线所包围的集水区域的总称。区域水土保持的重要工作之一是小流域治理，其产流汇流过程的模拟均离不开小流域这一精确划分治理模型及其评价的基本单元。因此，流域的自动提取在实际应用中具有重要意义。

目前，空间数据获取技术的发展，使大范围高分辨率数字地形数据的快速获取成为现实，为基于海量DEM数据的数字地形分析提供了丰富的数据源。众多应用领域对DEM数据分辨率及数字地形分析效率的要求也越来越高，相应的研究规模也不断增大。如在黄土高原地区，覆盖全区的采用5m分辨率DEM数据量约100GB，很难在单机条件下完成流域提取的计算任务。在数据量从少量到海量、分辨率从低到高、分析结果从低精度到高精度的背景下，单机计算能力不足的问题日益突出。这就造成了计算资源低利用率和海量数据高分析效率要求之间的矛盾，也使得DEM数据得不到充分利用。

在解决基于规模DEM数据集的流域提取方面，Lars Arge（Arge L.,et al.Efficient FlowComputation on Massive Grid Terrain Datasets[J].Geoinformatica,2003,7(4):283-313.）等对I/O效率算法进行了研究，其处理的数据规模可以达到2GB。现有的大多数GIS软件中集成的流域特征提取功能都对经典算法进行改进，如ArcGIS和Rivertools等软件的流域特征提取功能都在某种程度采用了I/O高性能算法，以减少计算时间。对流域算法的研究和改进，在一定程度上提高了算法效率和数据处理规模，然而这些分析算法都是串行的，可扩展性差，在处理更大规模的地形数据时仍然存在较大困难。

并行计算技术的发展为突破流域提取中数据规模的限制这一瓶颈提供了良好的契机。近年来，相关流域并行研究一是并行流域提取算法方面，如周海芳等（周海芳,刘光明,郑明玲等.遥感图像自动配准的串行与并行策略研究[J].国防科技大学学报,2004,26(2):56-61.）提出了基于距离的流域变换的一种优化并行流域分割算法，Ortega（Ortega L,Rueda A.Parallel drainage network computation on CUDA[J].Computers &Geosciences,2010,36:171-178.）提出了基于CUDA的流域网络并行特征提取算法；二是集成并行流域计算系统方面，王光谦（王光谦，刘家宏.黄河数字流域模[J].水利水电技术,2006,37(2):15-21.）等在集群环境下构建了黄河数字流域模型的并行计算平台，并对动态调度进行了相关研究；在并行处理自动机上，

A New Algorithm ForExtraction of Continuous Channel Networks without Problematic Parallels fromHydrologically Corrected DEMS.2010,16(1).）采用DRAW（Drainage axis way）算法解决了流域分析中平行线的问题。以上研究注重了并行算法与流域模型的研究，而数字地形分析面临大数据集和处理过程中大量的中间计算数据问题，产生了对数据的高吞吐量处理要求，该方面也是提高计算性能重要环节。Jianya Gong（Jianya Gong,Jibo Xie,Extraction ofdrainage networks from large terrain datasets using high throughoutcomputing[J].Computers & Geosciences,35(2009):337-346.）采用高吞吐率计算，以流域边界为数据划分策略，从大幅地形数据集并行提取流域网络，但从多尺度角度进行数据重采样获取流域边界划分，无疑增加了计算的开销且无法进行负载均衡控制。

美国犹他州立大学的TAUDEM程序采用MPICH2并行库，实现了水文信息的并行提取和分析，但其仅支持TIFF数据格式，且缺少考虑并行环境中数据划分粒度所产生的系统负载均衡和计算安全、稳定等问题，因此该方法常很难直接应用于实际DEM数字地形分析。

以上所述的现有流域提取方法的不足，在面向大区域进行地形分析和工程应用中带来很大困难，在诸如黄土高原水土保持与防治等工程应用上甚至会造成相关部门的决策失误，从而给国家和人民造成巨大的经济损失。

发明内容

本发明所要解决的技术问题，在于克服现有技术存在的缺陷，提出一种能够在数据并行方式下的数据划分、融合策略及流域提取并行计算的技术流程。该发明能够在并行时充分考虑数据划分的粒度及并行I/O机制，在串行算法并行解析上能够顾及数据自身特征。本技术流程能够实现面向大区域流域提取的并行计算，对大规模的DEM的广泛应用有着十分重要的意义。

为实现上述发明目的，本发明采用的技术方案如下：

面向大区域流域提取的并行方法，包括如下步骤：

第一步、数据划分策略：采用行带通信域划分方法对DEM数据集进行划分，并按照以下数据划分粒度模型对划分粒度进行优化：

PGM=(D,C,A,R)

其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系；

第二步、按照数据划分策略和数据融合策略进行洼地填平并行计算，包括对洼地处理和平地处理两个过程；其中，数据融合策略是对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件；

第三步、按照数据划分策略和数据融合策略，基于洼地填平结果进行水流方向并行计算；

第四步、按照数据划分策略和数据融合策略，基于水流方向数据进行汇流累积并行计算；

第五步、按照数据划分策略和数据融合策略，基于水流方向和汇流累积数据，设置汇流阈值，进行河网水系并行计算，包括河网水系提取和河网水系编码两个过程；

第六步、按照数据划分策略和数据融合策略，基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。

本发明提出的面向大区域流域提取的数据划分、融合策略及串行算法并行解析的并行方法具有以下特点：

（1）本发明提出的数据划分粒度模型，在兼顾负载平衡及计算环境安全的同时，可通过对最大行内存限制的调整进行粒度最优评价。基于最优粒度的数据划分，使得计算系统中的资源得到最优效率的应用。该粒度模型为以后的数据地形分析并行计算提供了技术基础。

（2）本发明构建并实现的面向大区域海量DEM数据的流域提取并行方法，在顾及流域提取算法特征的基础上对串行算法进行了并行解析。同时，该方法可在数据划分后任务数大于计算节点的情况下，完成并行计算任务并保证负载均衡。本方法基于并行I/O方法，大幅度降低了因读取和写入数据对计算效率的影响。对上十种栅格数据的支持，扩展了该方法的应用广度。

附图说明

图1是本发明数据划分策略中二维数组在物理层按行存储；

图2是本发明数据划分策略中行带通信域划分示意图；

图3是数据划分粒度优化示例；

图4是本发明数据融合策略中锚点融合示意图；

图5是本发明并行I/O方法；

图6（a）是实施例中洼地示意图，（b）是实施例中平地示意图；

图7是实施例中洼地处理过程中3×3邻域窗口；

图8（a）是实施例中数据并行方式下洼地填平处理详细程序实现流程图，（b）是（a）中初始化A流程，（c）是（a）中洼地处理B流程，（d）是（a）是平地处理C流程（图中洼①指洼地处理算法步骤①，平①指平地处理算法步骤①）；

图9是实施例中水流方向值示意图；

图10是实施例中数据并行方式下水流方向处理流程图；

图11是实施例中流向单元类型；

图12（a）是实施例中数据并行方式下汇流累积处理详细程序实现流程图,(b)是（a）中汇流累积计算A流程，（c）是（a）中汇流累积计算B流程（图中汇①指汇流累积处理算法步骤①）；

图13是实施例中河网节点类型图；

图14是实施例中河网段Morton编码示意图；

图15（a）是实施例中数据并行方式下河网水系确定处理详细程序实现流程图,(b)是（a）中河网水系编码A流程，（c）是（a）中河网水系编码B流程（图中提取①指河网水系提取算法步骤①，编码①指河网水系编码算法步骤①）；

图16（a）是实施例中水流方向矩阵，（b）实施例中流域划分示意图；

图17（a）是实施例中数据并行方式下河网水系确定处理详细程序实现流程图，（b）是（a）中初始化A流程，（c）是（a）中流域划分B流程（图中流①指流域划分算法步骤①）；

图18是实施例中流域并行提取流程图；

图19是洼地填平结果；

图20是实施例中水流方向结果；

图21是实施例中汇流累积结果；

图22是实施例中河网水系结果；

图23是实施例中流域提取结果。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细说明。

本发明中，数据并行方式下的数据划分、融合策略及并行I/O机制：

一、数据划分策略

（1）划分方式

行带通信域划分（Row Band Communication Partition,RBCP）：DEM数据在进行计算时，是一个m行×n列的矩阵形式，可以用二维数组进行组织。基于二维数组的DEM数据在物理层的存储形式为按行存储（如图1所示）且负载均衡是并行计算过程中各计算节点达到最大计算效率的重要因素两个方面的考虑，本发明进行DEM数据划分时，使用行带冗余划分。行带通信域划分是指依据一定的划分粒度，按行对DEM进行包含节点间通信域的带状划分方式。通信域是流域提取并行计算过程中各节点间进行信息交互的通道。行带通信域划分方式中，通信域的实施方案为在划分粒度的基础上，行带块上下各冗余一行。DEM数据划分后的行带块存在上、下两个通信域，其中首个行带块无上通信域，最后行带块无下通信域（如图2所示）。

通信域的本质在于计算数据的通信交互。流域提取算法具有计算密集的特点，故在其并行计算过程中，需要进行更新通信域中的数据。更新方法为：P_i计算节点的最后一行传递给P_i+1的上通信域，P_i+1处理节点的最上一行传递给P_i的下通信域（i为0、1…n）。

（2）划分粒度模型

数据划分粒度模型（Partition Granularity Model,PGM）：数据划分粒度是并行计算中数据划分的基本单位，是综合数据、计算环境、任务属性的结果。本发明针对流域提取的特征，数据划分粒度模型采用如下形式化定义：

PGM=(D,C,A,R)                （1）

其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系。

1）数据粒度。格网DEM是以规则栅格单元来描述地表的，目前规则栅格单元DEM的存储，一般采用格网原点坐标值(数据类型为浮点型，float)，无值区的标识值（数据类型为浮点型，float）以及格网间距(数据类型为整型，float)、行列数(数据类型为整型，int)，同时高程值集按照浮点型予以存储(数据类型为整型，float)。因此，DEM数据粒度可表达为：

G_d=4*sizeof(float)+2*sizeof(int)+row*col*sizeof(float)    （2）

式中，sizeof为内存容量度量函数，row为格网DEM行数，col为格网DEM列数。

2）计算粒度。并行计算系统中，处理能力随着计算节点数的增大而提高。然在处理海量DEM数据时，各计算节点内存资源有限，即所有的计算不能一次性的全部执行，需要按照最大内存限制对各计算节点的数据分发进行控制。结合DEM数据结构和数据划分方式，最大内存限制可转为最大内存行限制以方便计算管理。因此，并行计算粒度可表达为：

$\left\{\begin{array}{c}{G}_c=n*M_{\mathrm{limit}}\\ M_{\mathrm{limit}}=\mathrm{Maxrow}*\mathrm{col}*\mathrm{sizeof}\left(\mathrm{float}\right)\end{array}\right.---\left(3\right)$

式中，n为并行计算系统中计算节点个数，M_limit为最大内存限制，Maxrow为最大内存行限制，col为格网DEM列数，sizeof为内存容量度量函数。

3）任务属性是指DEM数据划分后的各种附加信息，包含分辨率、高程点数据类型、行列数、无值区标识值，通信域及最大内存限制等。

4）划分粒度评价。由于通信及流域提取中算法的全局性问题各节点处理相同大小数据划分粒度需执行的算法迭代次数也各不相同，这就导致同样大小粒度在不同的处理单元、不同地形分析算法的计算时间各不相同。对各节点的数据分配量进行优化可以优化提高系统的总执行效率，各节点按照最优的执行效率并行执行计算任务。系统对整个数据分配的要求主要是各节点计算的划分粒度不能超过各节点的执行能力、最大完成时间，破坏各节点的负载平衡等。划分粒度后的任务数为：

划分粒度后各计算时间估算表达式为：

$\left\{\begin{array}{c}{T}_{\mathrm{al}{l}_i}=T_{\mathrm{rea}{d}_i}+T_{\mathrm{writ}{e}_i}+T_{\mathrm{comput}{e}_i}\\ T_{\mathrm{all}}=\underset{i=1}{\overset{N_m}{\mathrm{\Σ}}}{T}_{\mathrm{alli}}\end{array}\right.---\left(5\right)$

式中T_read为读取文件时间，T_write为写文件时间；T_compute为任务计算时间。

通过对比T_all可以对数据划分粒度进行评价，在最大内存行限制Maxrow范围改变内存行限制Row_limit的值进行优化（如图3所示）。步骤为：

（a）初始化Row_limit∈(0,Max_row)；

（b）当一次任务计算完成后，记录T_alli，如果T_alli<T_alli-1，则Row_limit=Row_limit*2；否则其值不变；

（c）重复上一步，直到T_all停止减少或者达到最大内存行限制。

二、数据融合策略

锚点融合（Anchor Merge,AM）：对大规模DEM数据集进行划分，并经流域提取后，需要对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件。由于采用了行带通信域划分的方式，不同计算节点的子结果数据集在结果文件中具有唯一的空间位置。如图4所示，各计算节点可以根据锚点和子结果数据集的大小进行无缝融合，锚点定义如下：

锚点（Anchor）：当数据按一维行规则划分后，子数据集通过其左上角的起始行列号在大数据集中进行空间定位。此时起始点作用如同船锚的稳固作用，故称之为锚点。

根据行带通信域划分策略，可以对计算节点的锚点(A_col,A_row)进行计算，公式如下：

A_col=0

三、并行I/O方法

数据划分解决了并行化中各计算节点所需数据，但如何将划分后的子数据集分发给相应的计算节点是并行计算过程中重要的方面。并行I/O方法能够大幅度降低主从式因分发通信而带来的各计算节点数据输入时间，从而提高并行性能。

如图5所示，管理节点、计算节点和共享存储系统等部分组成并行I/O的硬件部分，计算节点以共享存储系统中数据的空间分布元数据作为对象与管理节点进行通信，即各计算节点可以通过高速专网访问共享存储系统，在控制流的引导下，获取相应的数据集。在共享存储系统环境下，各计算节点可以按照数据划分策略对数据集进行并行读取。并行读写方法为：依据各计算节点的锚点和划分粒度确定空间位置读写。基于开源库GDAL对12种常用的数据格式进行控制，扩展本发明的应用范围。

四、主要的类设计

本发明中设计和实现的过程中主要的类有7类：RasterClass用来标准不同格式的DEM数据，同时提供了获取DEM数据的头文件信息、局部高程计算等功能；RasterIOClass实现了基于并行I/O方法的DEM数据的并行读取与输出及不同数据格式DEM之间的转换；PartitionClass根据数据划分及融合策略实现了并行数据划分及融合及任务属性管理；FillClass实现了流域提取过程中并行洼地填平功能；FlowDirectionClass执行着流域提取过程中的水流方向并行计算功能；FlowAccumulationClass实现了流域提取过程中的汇流累积并行计算；StreamClass实现了流域提取过程中的河网水系并行计算；WatershedClass实现了并行流域提取过程中的河网水系段上游汇水区并行计算。

本发明面向大区域流域提取并行方法，主要技术流程如下：

第一步、洼地填平：在进行流域提取时，DEM数据中的洼地和平地两种特殊地形必须先处理以确保每个DEM栅格单元拥有确定的水流方向。洼地是被较高高程所包围的DEM局部地形单元，而平地是DEM数据中的平地（如图6所示）。自然条件下，水流向低处流动，遇到洼地，首先将其填满，然后再从该洼地的某一最低出口流出。由于洼地是局部的最低点，所以无法确定该点的水流方向。因此，洼地对确定水流的方向有重要的影响。同理，平地区域的存在对水流方向的确定也有着重要的影响。

洼地处理算法如下：

①初始化洼地处理结果栅格数据FillRaster（RasterClass类）：如果是栅格边界上的点，则赋予原始DEM的高程值，否则赋予一个极高值FLOAT_MAX；

②初步查找潜在洼地栈pitlist1：遍历结果栅格数据FillRaster中的单元（irow，jcol），执行洼地处理操作。

洼地处理操作为：查找栅格单元8个邻域（如图7所示）中的最小值MinNeighbor，如果原始DEM对应单元值不小于MinNeighbor，则FillRaster中该单元值更新为原始DEM对应值；否则该点加入pitlist1，同时如果FillRaster对应单元值大于MinNeighbor，则该单元值更新为MinNeighbor；

③处理潜在洼地栈pitlist1：循环依次取出pitlist1中的潜在洼地单元，对其进行洼地处理操作，直到pitlist1为空。

洼地出理操作为：查找该栅格单元8个邻域中的最小值MinNeighbor，如果原始DEM对应单元值不小于MinNeighbor，则FillRaster中该单元值更新为原始DEM对应值；否则该点加入另一个潜在洼地栈pitlist2，同时如果FillRaster对应单元值大于MinNeighbor，则该单元值更新为MinNeighbor；

④处理潜在洼地栈pitlist2：循环依次取出pitlist2中的潜在洼地单元，对其进行洼地处理操作，直到pitlist2为空。

洼地出理操作为：查找该栅格单元8个邻域中的最小值MinNeighbor，如果原始DEM对应单元值不小于MinNeighbor，则FillRaster中该单元值更新为原始DEM对应值；否则该点加入另一个潜在洼地栈pitlist1，同时如果FillRaster对应单元值大于MinNeighbor，则该单元值更新为MinNeighbor；

⑤重复③、④，直到pitlist1、pitlist2均为空。

平地处理算法如下：

①遍历经过洼地处理后的DEM数据，查找8个邻域栅格高程都不低于该栅格高程的栅格点，标记为平地单元，加入平地栈flatlist；

②循环依次取出flatlist中的单元，给每一个栅格单元都增加一个微小的增量（如栅格高程采样精度的十分之一、千分之一或万分之一）；

③重复上述①、②过程，直到再也搜索不到平地单元。

在数据并行方式下洼地填平处理详细程序实现流程图如图8所示。

第二步、水流方向确定：水流方向是指水流离开格网时的流向。流向确定目前有单流向和多流向两种，但在流域分析中，常是在3×3局部窗口中通过D8算法确定水流方向。算法原理如下：

在3×3局部窗口（图7）中，设中心格网为c，其流向（即水流的流出方向）在其相邻八个格网点i中选择，i满足条件：

Max{d×(z_c-z_i)}    (i＝1,2,...,8)            (7)

其中，当i位于水平方向（东西）或垂直方向（南北）时，d=1；当i为对角线方向时，

此时，水流方向结果栅格该单元FlowDir值由下式确定（如图9所示）：

FlowDir=2^i-1            (8)

在数据并行方式下水流方向处理流程图如图10所示，详细程序实现步骤为：

①数据划分：根据数据划分对洼地填平后的DEM数据进行数据划分，同时创建水流方向结果栅格文件FlowDirRaster；

②循环任务：根据任务数N_m和计算节点数n确定计算次数block，循环计算次数依次计算各节点数据：

(a)遍历FillRaster数据，每个栅格单元根据水流方向计算公式（7）、（8）进行流向值FlowDir；

（b）根据数据融合策略，将结果数据写入FlowDirRaster；

③水流方向完成。

第三步、汇流累积计算：汇流累积是指流向该格网的所有的上游格网单元的水流累计量（将格网单元看作是等权的，以格网单元的数量或面积计），它是基于水流方向确定的，是流域划分的基础。根据水流方向计算方法，假设每个栅格单元为1个单位流量，本发明的汇流累积计算公式为：

$\mathrm{FlowAcc}\left(i\right)=k+\underset{k}{\mathrm{\&Sigma;}}\mathrm{FlowAcc}\left(k\right)---\left(9\right)$

式中，i指当前处理单元，k为8邻域中流入i的栅格。

水流方向确定后，可将流向单元分为三类（如图11所示）：（1）源单元（Source Cell,SC）。该类单元是指无邻域单元流入。（2）径流单元（Runoff Cell,RC）。此类单元是指有一个邻域单元流入但流向不改变。（3）低洼单元（Sink Cell，SK）。除了上述两类，都属于该类单元。

根据汇流累积计算公式可知，只有当所有流入本单元的邻域单元汇流累积计算完毕时，本单元方可进行计算。由上述流向单元分类，本发明汇流累积算法如下：

①初始化汇流累积结果栅格FlowAccRaster：遍历水流方向栅格FlowDirRaster，如果是SC单元，则赋值于-1且加入栈SClist；如果是RC单元，则赋值于-2;如果是SK单元，则赋值-3-m且加入栈SKlist（m为流入SK单元的邻域单元数，值域范围为[1，7]）；

②源单元流径累积：依次遍历栈SClist元素，从SC单元开始计算其汇流累积值，按照水流方向向下游追溯直到遇到SK单元停止，同时将该SK单元FlowAccRaster的FlowAcc值更新为FlowAcc+1，最后将该SC单元从SClist中移除；

③低洼单元流径累积：遍历栈SKlist元素，如果SK单元FlowAcc==-3，则从SK单元开始计算其汇流累积值，按照水流方向向下追溯直到遇到SK单元停止，同时将该SK单元FlowAccRaster的FlowAcc值更新为FlowAcc+1且如果FlowAcc==-3,则将SK加入栈SKlist，最后将该SK单元从SKlist中移除后重新遍历栈SKlist；如果SK单元FlowAcc<>-3，则不处理；

④当SKlist为空时，任务完成。

在数据并行方式下汇流累积处理详细程序实现流程图如图12所示。

第四步、河网水系确定：河网水系是根据汇流累积大于一定阈值提取而来，同时对每一河网段进行分配一个唯一的识别符。如图13所示，河网节点包括源点和汇合点两类。源点是外部河网段的起始点，没有上游流入河网；汇合点是几个河网的交汇点，汇合点接受多个河网的流量。因此，河网水系确定分为河网水系提取和河网水系编码两个部分。本发明采用Morton编码对河网段进行编码，具体方法为：以河网段的起始单元行、列号的Morton作为整个河网段的编码（如图14所示）。

河网水系提取算法为：

①遍历FlowAccRaster数据，如果汇流累积值不小于给定参数阈值threshold，则河网栅格StreamRaster数据对应单元赋值-1，否则赋值Nodata；

②写入文件，完成河网水系提取。

河网水系编码算法为：

①查找源点和汇合点：遍历StreamRaster，3×3邻域窗口中判断河网节点类型。如果该单元值为-1且无上游单元流入，则加入源点栈SNlist；如果该单元值为-1且有2个以上上游单元流入，则加入汇合点栈JNlist；

②源点河网编码：依次取SNlist栈中的单元，根据水流方向FlowDirRaster追溯下游单元，直到遇到汇合点（或下游单元越出范围）为止，然后计算Morton码值M，将该路径上所有单元StreamRaster值更新为M；最后从SNlist移除该单元。

③汇合点河网编码：遍历取JNlist栈中的单元，依次沿水流方向FlowDirRaster追溯汇入该汇合点的路径上游单元，如果遇到该上游已编码M的单元，则反向向下游追溯进行该路径编码直到该汇合点；如果该汇合点的所以上游路径都已编码，则将该单元作为源点进行源点河网编码操作后从JNlist中移除该单元重新遍历JNlist栈，否则遍历执行下一个汇合点。

④当JNlist为空时，任务完成。

在数据并行方式下河网水系确定处理详细程序实现流程图如图15所示。

第五步、流域划分：格网单元分配，即流域划分，是以河网段为依据，把流向该河网段的所有上游格网标识出来，从而实现子流域的划分。每个河网段的汇水区单元都标识为河网段编码值（如图16所示）。定义汇水种子为邻域的下游单元已标记且自身未标记的格网单元。流域划分算法为：

①初始化结果栅格WatershedRaster数据：遍历StreamRaster，将WatershedRaster对应单元赋予StreamRaster值；

②格网单元分配：遍历WatershedRaster，根据水流方向FlowDirRaster数据，判断当前单元（irow，jcol）是否为汇水种子，如果是，则进行汇水区计算：

（a）标识当前单元（irow，jcol）；

（b）依次判断单元（irow，jcol）的8个邻域单元是否为汇水种子，如果是，则进行汇水区计算；

③重复②，直到无法找到汇水种子为止。

在数据并行方式下河网水系确定处理详细程序实现流程图如图17所示。

以上第一～第五步操作实现了流域并行提取技术，完成了面向大区域DEM流域提取应用（如图18所示）。

下面以黄土高原为例。

黄土高原东西千余千米，南北700千米，面积约40万平方千米，沟壑纵横，是我国水土流失最为严重的区域。以小流域为单位对黄土高原水土流失进行治理，具有良好的效果。国家现有的1:10000比例尺DEM规定的采样间距（格网大小）是5米，覆盖整个黄土高原地区数据量约100GB，很难在单机条件下完成流域提取的计算任务。单机计算的受阻直接造成现有DEM在该地区难以展开有效的地学分析与工程应用。实验数据选取覆盖整个黄土高原的1:50000比例尺DEM数据，约4GB，实验计算环境:4台高性能机架服务器，型号为Blade HS22（2*Intel Xeon CPUs 2.8GHz，6*2GB内存）；采用GDAL1.8开源库和MPICH2并行计算平台。

流域提取并行计算过程为：

第一步、最优划分粒度：采用较小数据量（实验数据中的部分图幅数据）的数据，利用本发明设计数据划分粒度模型对计算系统的最优粒度进行评价。

第二步、洼地填平：调用FillClass过程，完成洼地和平地处理过程。

第三步、流向计算：调用FlowDirectionClass过程，完成水流方向计算。

第四步、汇流累积计算：调用FlowAccumulationClass过程，完成汇流累积计算。

第五步、河网水系：调用StreamClass过程，完成河网水系计算。

第六步、流域划分：调用WatershedClass过程，完成实验区域流域提取计算。

提取过程中，洼地填平、水流方向、汇流累积、河网水系及流域提取结果如图19、20、21、22、23所示。

Claims (5)

1.面向大区域流域提取的并行方法，其特征在于，包括如下步骤：

第一步、数据划分策略：采用行带通信域划分方法对DEM数据集进行划分，并按照以下数据划分粒度模型对划分粒度进行优化：

PGM=(D,C,A,R)

其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系；

第二步、按照数据划分策略和数据融合策略进行洼地填平并行计算，包括对洼地处理和平地处理两个过程；其中，数据融合策略是对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件；

第三步、按照数据划分策略和数据融合策略，基于洼地填平结果进行水流方向并行计算；

第四步、按照数据划分策略和数据融合策略，基于水流方向数据，进行汇流累积并行计算；

第五步、按照数据划分策略和数据融合策略，基于水流方向和汇流累积数据，设置汇流阈值，进行河网水系并行计算；

第六步、按照数据划分策略和数据融合策略，基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。

2.根据权利要求1所述的面向大区域流域提取的并行方法，其特征在于，所述第三步的具体计算过程如下：

步骤31.数据划分：根据数据划分策略对经洼地填平后的DEM数据进行数据划分，同时创建水流方向结果栅格文件；

步骤32.循环任务：根据任务数N_m和计算节点数n确定计算次数，循环计算次数依次计算各计算节点数据：（a）遍历初始化过的洼地填平结果栅格数据，每个栅格单元根据水流方向计算公式进行流向值；（b）根据数据融合策略，将各计算节点的结果数据写入水流方向结果栅格文件。

3.根据权利要求1所述的面向大区域流域提取的并行方法，其特征在于，所述第四步的具体计算过程如下：

步骤41.初始化汇流累积结果栅格数据：遍历水流方向结果栅格数据，如果是SC单元，则赋值于-1且加入栈SClist；如果是RC单元，则赋值于-2;如果是SK单元，则赋值-3-m且加入栈SKlist，m为流入SK单元的邻域单元数，值域范围为[1，7]；

步骤42.源单元流径累积：依次遍历栈SClist元素，从SC单元开始计算其汇流累积值，按照水流方向向下游追溯直到遇到SK单元停止，同时将该SK单元的汇流累积值FlowAcc更新为FlowAcc+1，最后将该SC单元从栈SClist中移除；

步骤43.低洼单元流径累积：遍历栈SKlist元素，如果SK单元的汇流累积值FlowAcc==-3，则从SK单元开始计算其汇流累积值，按照水流方向向下追溯直到遇到SK单元停止，同时将该SK单元汇流累积值FlowAcc更新为FlowAcc+1且如果汇流累积值FlowAcc==-3,则将SK加入栈SKlist，最后将该SK单元从栈SKlist中移除后重新遍历栈SKlist；如果SK单元的汇流累积值FlowAcc<>-3，则不处理；

步骤44.当栈SKlist为空时，任务完成。

4.根据权利要求1所述的面向大区域流域提取的并行方法，其特征在于，所述第五步的具体计算过程如下：

步骤51.提取河网水系；

步骤52.查找源点和汇合点：遍历河网水系栅格数据，3×3邻域窗口中判断河网节点类型：如果该单元值为-1且无上游单元流入，则加入源点栈SNlist；如果该单元值为-1且有2个以上上游单元流入，则加入汇合点栈JNlist；

步骤53.源点河网编码：依次取源点栈SNlist中的元素，根据水流方向栅格数据追溯下游单元，直到遇到汇合点或下游单元越出范围为止，然后计算Morton码值M，将该路径上所有单元河网水系栅格数据值更新为M；最后从源点栈SNlist移除该单元；

步骤54.汇合点河网编码：遍历取汇合点栈JNlist中的元素，依次沿水流方向栅格数据追溯汇入该汇合点的路径上游单元，如果遇到该上游已编码M的单元，则反向向下游追溯进行该路径编码直到该汇合点；如果该汇合点的所以上游路径都已编码，则将该单元作为源点进行源点河网编码操作后从汇合点栈JNlist中移除该单元并重新遍历汇合点栈JNlist，否则遍历执行下一个汇合点；

步骤55.当汇合点栈JNlist为空时，任务完成。

5.根据权利要求1所述的面向大区域流域提取的并行方法，其特征在于，所述第六步的具体计算过程如下：

步骤61.初始化流域划分结果栅格数据：遍历河网水系栅格数据，将流域划分结果栅格对应单元赋予河网栅格数据值；

步骤62.格网单元分配：遍历流域划分结果栅格数据，根据水流方向栅格数据，判断当前单元（irow，jcol）是否为汇水种子，如果是，则进行汇水区计算：（a）标识当前单元（irow，jcol）；（b）依次判断单元（irow，jcol）的8个邻域单元是否为汇水种子，如果是，则进行汇水计算；

步骤63.重复步骤62，直到无法找到汇水种子为止。

Images