专利名称：面向大区域流域提取的并行方法
技术领域：
本发明涉及一种数字高程模型(DEM)的流域提取并行技术，具体说是一种基于面向大区域海量DEM数据的数据划分、融合策略及串行算法并行解析机制的流域提取并行方法。
背景技术：
数字高程模型产品(Digital Elevation Models，缩写为DEM)是地理信息系统(Geographical Information System,缩写为GIS)地理数据库中最为重要的空间信息资料和赖以进行地形分析的核心数据系统，其应用遍布测绘、交通、军事、水利、农业、环境、资源管理、规划与旅游等众多领域，如在测绘部门主要用于三维地形建模、可视域分析、工程土 方估算、遥感影像的几何矫正等；在地学界主要用于水文分析、蓄洪计算、淹没分析等三维地学分析；在国土资源部门主要用于土地利用类型详察等。DEM并不是地形的简单表达，真正价值在于DEM所蕴含的各种地学信息。如何进行“数据挖掘”，获取各种反映实际地形地貌的结构信息，是GIS领域一个重要的研究方向。在此需求环境下，数字地学分析技术应运而生。数字地形分析(Digital TerrainAnalysis, DTA)是在DEM上进行地形属性计算和特征提取的数字信息处理技术。在复杂的现实世界地理过程中各影响因子和简单、高效、精确和易于理解的抽象和计算机实现中找到平衡，是数字地形分析的核心任务。用来描述地形特征和空间分布的地形参数有很多，根据地形属性的地学应用范畴，地形属性可分为一般地形属性(如高程、坡度坡向等)和水文特征(如地形结构识别、流域提取等)。其中，流域提取是DTA的重要研究内容。流域是指地表水及地下水分水线所包围的集水区域的总称。区域水土保持的重要工作之一是小流域治理，其产流汇流过程的模拟均离不开小流域这一精确划分治理模型及其评价的基本单元。因此，流域的自动提取在实际应用中具有重要意义。目前，空间数据获取技术的发展，使大范围高分辨率数字地形数据的快速获取成为现实，为基于海量DEM数据的数字地形分析提供了丰富的数据源。众多应用领域对DEM数据分辨率及数字地形分析效率的要求也越来越高，相应的研究规模也不断增大。如在黄土高原地区，覆盖全区的采用5m分辨率DEM数据量约100GB，很难在单机条件下完成流域提取的计算任务。在数据量从少量到海量、分辨率从低到高、分析结果从低精度到高精度的背景下，单机计算能力不足的问题日益突出。这就造成了计算资源低利用率和海量数据高分析效率要求之间的矛盾，也使得DEM数据得不到充分利用。在解决基于规模DEM数据集的流域提取方面，Lars Arge (Arge L.，etal. Efficient Flow Computation on Massive Grid Terrain Datasets[J]. Geoinformatica, 2003,7(4) :283-313.)等对1/0效率算法进行了研究，其处理的数据规模可以达到2GB。现有的大多数GIS软件中集成的流域特征提取功能都对经典算法进行改进，如ArcGIS和Rivertools等软件的流域特征提取功能都在某种程度采用了 1/0高性能算法，以减少计算时间。对流域算法的研究和改进，在一定程度上提高了算法效率和数据处理规模，然而这些分析算法都是串行的，可扩展性差，在处理更大规模的地形数据时仍然存在较大困难。并行计算技术的发展为突破流域提取中数据规模的限制这一瓶颈提供了良好的契机。近年来，相关流域并行研究一是并行流域提取算法方面，如周海芳等(周海芳，刘光明，郑明玲等.遥感图像自动配准的串行与并行策略研究[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.)等在集群环境下构建了黄河数字流域模型的并行计算平台，并对动态调度进行了相关研究；在并行处理自动机上，GilLGEN (Gi)LGEN, G0KG0Z. A NewAlgorithm For Extraction of Continuous Channel Networks without ProblematicParallels from Hydrologically Corrected DEMS. 2010，16 (I) ·)米用 DRAW (Drainage axis way)算法解决了流域分析中平行线的问题。以上研究注重了并行算法与流域模型的研究，而数字地形分析面临大数据集和处理过程中大量的中间计算数据问题，产生了对数据的高吞吐量处理要求，该方面也是提高计算性能重要环节。Jianya Gong (JianyaGong, Jibo Xie, Extraction of drainage networks from large terrain datasets usinghigh throughout computing [J]. Computers & Geosciences, 35 (2009) : 337-346.)米用高吞吐率计算，以流域边界为数据划分策略，从大幅地形数据集并行提取流域网络，但从多尺度角度进行数据重采样获取流域边界划分，无疑增加了计算的开销且无法进行负载均衡控制。美国犹他州立大学的TAUDEM程序采用MPICH2并行库，实现了水文信息的并行提取和分析，但其仅支持TIFF数据格式，且缺少考虑并行环境中数据划分粒度所产生的系统负载均衡和计算安全、稳定等问题，因此该方法常很难直接应用于实际DEM数字地形分析。以上所述的现有流域提取方法的不足，在面向大区域进行地形分析和工程应用中带来很大困难，在诸如黄土高原水土保持与防治等工程应用上甚至会造成相关部门的决策失误，从而给国家和人民造成巨大的经济损失。
发明内容
本发明所要解决的技术问题，在于克服现有技术存在的缺陷，提出一种能够在数据并行方式下的数据划分、融合策略及流域提取并行计算的技术流程。该发明能够在并行时充分考虑数据划分的粒度及并行I/o机制，在串行算法并行解析上能够顾及数据自身特征。本技术流程能够实现面向大区域流域提取的并行计算，对大规模的DEM的广泛应用有着十分重要的意义。为实现上述发明目的，本发明采用的技术方案如下面向大区域流域提取的并行方法，包括如下步骤第一步、数据划分策略采用行带通信域划分方法对DEM数据集进行划分，并按照以下数据划分粒度模型对划分粒度进行优化PGM= (D, C，A, R)其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系；
第二步、按照数据划分策略和数据融合策略进行洼地填平并行计算，包括对洼地处理和平地处理两个过程；其中，数据融合策略是对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件；第三步、按照数据划分策略和数据融合策略，基于洼地填平结果进行水流方向并行计算；第四步、按照数据划分策略和数据融合策略，基于水流方向数据进行汇流累积并行计算；第五步、按照数据划分策略和数据融合策略，基于水流方向和汇流累积数据，设置 汇流阈值，进行河网水系并行计算，包括河网水系提取和河网水系编码两个过程；第六步、按照数据划分策略和数据融合策略，基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。本发明提出的面向大区域流域提取的数据划分、融合策略及串行算法并行解析的并行方法具有以下特点(I)本发明提出的数据划分粒度模型，在兼顾负载平衡及计算环境安全的同时，可通过对最大行内存限制的调整进行粒度最优评价。基于最优粒度的数据划分，使得计算系统中的资源得到最优效率的应用。该粒度模型为以后的数据地形分析并行计算提供了技术基础。(2)本发明构建并实现的面向大区域海量DEM数据的流域提取并行方法，在顾及流域提取算法特征的基础上对串行算法进行了并行解析。同时，该方法可在数据划分后任务数大于计算节点的情况下，完成并行计算任务并保证负载均衡。本方法基于并行I/O方法，大幅度降低了因读取和写入数据对计算效率的影响。对上十种栅格数据的支持，扩展了该方法的应用广度。
图I是本发明数据划分策略中二维数组在物理层按行存储；图2是本发明数据划分策略中行带通信域划分示意图；图3是数据划分粒度优化示例；图4是本发明数据融合策略中锚点融合示意图；图5是本发明并行I/O方法；图6 (a)是实施例中挂地示意图，(b)是实施例中平地示意图；图7是实施例中洼地处理过程中3X3邻域窗口 ；图8 Ca)是实施例中数据并行方式下洼地填平处理详细程序实现流程图，(b)是Ca)中初始化A流程，(C)是(a)中洼地处理B流程，(d)是(a)是平地处理C流程(图中洼①指洼地处理算法步骤①，平①指平地处理算法步骤①);图9是实施例中水流方向值示意图；图10是实施例中数据并行方式下水流方向处理流程图；图11是实施例中流向单元类型；图12 (a)是实施例中数据并行方式下汇流累积处理详细程序实现流程图，(b)是Ca)中汇流累积计算A流程，(c)是(a)中汇流累积计算B流程(图中汇①指汇流累积处理算法步骤①);图13是实施例中河网节点类型图；图14是实施例中河网段Morton编码示意图；图15 Ca)是实施例中数据并行方式下河网水系确定处理详细程序实现流程图，(b)是(a)中河网水系编码A流程，(C)是(a)中河网水系编码B流程(图中提取①指河网水系提取算法步骤①，编码①指河网水系编码算法步骤①);图16 Ca)是实施例中水流方向矩阵，(b)实施例中流域划分示意图；图17(a)是实施例中数据并行方式下河网水系确定处理详细程序实现流程图，(b)是(a)中初始化A流程，(C)是(a)中流域划分B流程(图中流①指流域划分算法步骤①)； 图18是实施例中流域并行提取流程图；图19是洼地填平结果；图20是实施例中水流方向结果；图21是实施例中汇流累积结果；图22是实施例中河网水系结果；图23是实施例中流域提取结果。
具体实施例方式下面结合附图和具体实施例，对本发明作进一步详细说明。本发明中，数据并行方式下的数据划分、融合策略及并行I/O机制一、数据划分策略(I)划分方式行带通信域划分(RowBand Communication Partition, RBCP):DEM 数据在进行计算时，是一个m行Xn列的矩阵形式，可以用二维数组进行组织。基于二维数组的DEM数据在物理层的存储形式为按行存储(如图I所示)且负载均衡是并行计算过程中各计算节点达到最大计算效率的重要因素两个方面的考虑，本发明进行DEM数据划分时，使用行带冗余划分。行带通信域划分是指依据一定的划分粒度，按行对DEM进行包含节点间通信域的带状划分方式。通信域是流域提取并行计算过程中各节点间进行信息交互的通道。行带通信域划分方式中，通信域的实施方案为在划分粒度的基础上，行带块上下各冗余一行。DEM数据划分后的行带块存在上、下两个通信域，其中首个行带块无上通信域，最后行带块无下通信域(如图2所示)。通信域的本质在于计算数据的通信交互。流域提取算法具有计算密集的特点，故在其并行计算过程中，需要进行更新通信域中的数据。更新方法为=Pi计算节点的最后一行传递给Pi+1的上通信域，pi+1处理节点的最上一行传递给Pi的下通信域(i为0、1···η)。(2)划分粒度模型数据划分粒度模型(Partition Granularity Model, PGM):数据划分粒度是并行计算中数据划分的基本单位，是综合数据、计算环境、任务属性的结果。本发明针对流域提取的特征，数据划分粒度模型采用如下形式化定义PGM= (D, C，A, R)(I)其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系。
I)数据粒度。格网DEM是以规则栅格单元来描述地表的，目前规则栅格单元DEM的存储，一般采用格网原点坐标值(数据类型为浮点型，float)，无值区的标识值(数据类型为浮点型，float)以及格网间距(数据类型为整型，float)、行列数(数据类型为整型，int)，同时高程值集按照浮点型予以存储(数据类型为整型，float)。因此，DEM数据粒度可表达为Gd=4*sizeof(float)+2*sizeof(int)+row*col*sizeof(float) (2)式中，sizeof为内存容量度量函数，row为格网DEM行数，col为格网DEM列数。2)计算粒度。并行计算系统中，处理能力随着计算节点数的增大而提高。然在处理海量DEM数据时，各计算节点内存资源有限，即所有的计算不能一次性的全部执行，需要按照最大内存限制对各计算节点的数据分发进行控制。结合DEM数据结构和数据划分方式，最大内存限制可转为最大内存行限制以方便计算管理。因此，并行计算粒度可表达为
{G, = n* Mlimit·I c(3)
I = Maxrow * co/5,< sizeof { float)式中，n为并行计算系统中计算节点个数，Mlimit为最大内存限制,Maxrow为最大内存行限制，col为格网DEM列数，sizeof为内存容量度量函数。3)任务属性是指DEM数据划分后的各种附加信息，包含分辨率、高程点数据类型、行列数、无值区标识值，通信域及最大内存限制等。4)划分粒度评价。由于通信及流域提取中算法的全局性问题各节点处理相同大小数据划分粒度需执行的算法迭代次数也各不相同，这就导致同样大小粒度在不同的处理单元、不同地形分析算法的计算时间各不相同。对各节点的数据分配量进行优化可以优化提高系统的总执行效率，各节点按照最优的执行效率并行执行计算任务。系统对整个数据分配的要求主要是各节点计算的划分粒度不能超过各节点的执行能力、最大完成时间，破坏各节点的负载平衡等。划分粒度后的任务数为
N(^%GC=0)
m \GdiGc\*n+n {Gd % Gc ^ 0)划分粒度后各计算时间估算表达式为
"rT 二 y1I J1* J1
all i — read j write i compute ^\(5)
Tau —艺 ^7/"
,f 二 I式中TMad为读取文件时间，Twrite为写文件时间；T_pute为任务计算时间。通过对比Tall可以对数据划分粒度进行评价,在最大内存行限制Maxrow范围改变内存行限制Rowlimit的值进行优化(如图3所示)。步骤为(a)初始化 Rowlimit e (O, Maxrow)；(b)当一次任务计算完成后,记录Talli,如果UTallH,则Rowlimit=Rowlimit*2 ;否则其值不变； ( c )重复上一步，直到Tall停止减少或者达到最大内存行限制。二、数据融合策略
锚点融合(Anchor Merge, AM):对大规模DEM数据集进行划分，并经流域提取后，需要对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件。由于采用了行带通信域划分的方式，不同计算节点的子结果数据集在结果文件中具有唯一的空间位置。如图4所示，各计算节点可以根据锚点和子结果数据集的大小进行无缝融合，锚点定义如下锚点(Anchor):当数据按一维行规则划分后，子数据集通过其左上角的起始行列号在大数据集中进行空间定位。此时起始点作用如同船锚的稳固作用，故称之为锚点。根据行带通信域划分策略，可以对计算节点的锚点(kcol，KJ进行计算，公式如下Acol=O
权利要求
1.面向大区域流域提取的并行方法，其特征在于，包括如下步骤 第一步、数据划分策略采用行带通信域划分方法对DEM数据集进行划分，并按照以下数据划分粒度模型对划分粒度进行优化 PGM= (D, C，A, R) 其中，D是数据粒度，C是计算粒度，A是任务的属性，R是划分粒度的评价系； 第二步、按照数据划分策略和数据融合策略进行洼地填平并行计算，包括对洼地处理和平地处理两个过程；其中，数据融合策略是对各计算节点的子结果数据集进行融合处理，形成一个完整的结果文件； 第三步、按照数据划分策略和数据融合策略，基于洼地填平结果进行水流方向并行计算； 第四步、按照数据划分策略和数据融合策略，基于水流方向数据，进行汇流累积并行计算； 第五步、按照数据划分策略和数据融合策略，基于水流方向和汇流累积数据，设置汇流阈值，进行河网水系并行计算； 第六步、按照数据划分策略和数据融合策略，基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。
2.根据权利要求I所述的面向大区域流域提取的并行方法，其特征在于，所述第三步的具体计算过程如下 步骤31.数据划分根据数据划分策略对经洼地填平后的DEM数据进行数据划分，同时创建水流方向结果栅格文件； 步骤32.循环任务根据任务数Nm和计算节点数η确定计算次数，循环计算次数依次计算各计算节点数据(a)遍历初始化过的洼地填平结果栅格数据，每个栅格单元根据水流方向计算公式进行流向值；(b)根据数据融合策略，将各计算节点的结果数据写入水流方向结果栅格文件。
3.根据权利要求I所述的面向大区域流域提取的并行方法，其特征在于，所述第四步的具体计算过程如下 步骤41.初始化汇流累积结果栅格数据遍历水流方向结果栅格数据，如果是SC单元，则赋值于-I且加入栈SClist ;如果是RC单元，则赋值于-2;如果是SK单元，则赋值-3-m且加入栈SKlist, m为流入SK单元的邻域单元数,值域范围为[1,7]； 步骤42.源单元流径累积依次遍历栈SClist元素，从SC单元开始计算其汇流累积值，按照水流方向向下游追溯直到遇到SK单元停止，同时将该SK单元的汇流累积值FIowAcc更新为FlowAcc+Ι,最后将该SC单元从栈SClist中移除； 步骤43.低洼单元流径累积遍历栈SKlist元素，如果SK单元的汇流累积值FlowAcc==-3,则从SK单元开始计算其汇流累积值，按照水流方向向下追溯直到遇到SK单元停止，同时将该SK单元汇流累积值FIowAcc更新为FlowAcc+Ι且如果汇流累积值FlowAcc==-3,则将SK加入栈SKlist，最后将该SK单元从栈SKlist中移除后重新遍历栈SKlist ;如果SK单元的汇流累积值FlowAcc〈>-3,则不处理； 步骤44.当栈SKlist为空时，任务完成。
4.根据权利要求I所述的面向大区域流域提取的并行方法，其特征在于，所述第五步的具体计算过程如下 步骤51.提取河网水系； 步骤52.查找源点和汇合点遍历河网水系栅格数据，3X3邻域窗口中判断河网节点类型如果该单元值为-I且无上游单元流入，则加入源点栈SNlist ;如果该单元值为-I且有2个以上上游单元流入，则加入汇合点栈JNlist ； 步骤53.源点河网编码依次取源点栈SNlist中的元素，根据水流方向栅格数据追溯下游单元，直到遇到汇合点或下游单元越出范围为止，然后计算Morton码值M，将该路径上所有单元河网水系栅格数据值更新为M ;最后从源点栈SNlist移除该单元； 步骤54.汇合点河网编码遍历取汇合点栈JNlist中的元素，依次沿水流方向栅格数据追溯汇入该汇合点的路径上游单元，如果遇到该上游已编码M的单元，则反向向下游追溯进行该路径编码直到该汇合点；如果该汇合点的所以上游路径都已编码，则将该单元作为源点进行源点河网编码操作后从汇合点栈JNlist中移除该单元并重新遍历汇合点栈JNlist,否则遍历执行下一个汇合点； 步骤55.当汇合点栈JNlist为空时，任务完成。
5.根据权利要求I所述的面向大区域流域提取的并行方法，其特征在于，所述第六步的具体计算过程如下 步骤61.初始化流域划分结果栅格数据遍历河网水系栅格数据，将流域划分结果栅格对应单元赋予河网栅格数据值； 步骤62.格网单元分配遍历流域划分结果栅格数据，根据水流方向栅格数据，判断当前单元(irow，jcol)是否为汇水种子，如果是，则进行汇水区计算(a)标识当前单元(irow, jcol)； (b)依次判断单元(irow, jcol)的8个邻域单元是否为汇水种子,如果是,贝U进行汇水计算； 步骤63.重复步骤62，直到无法找到汇水种子为止。
全文摘要
本发明公开了一种面向大区域海量DEM数据流域提取的并行方法，其步骤如下第一步、最优划分粒度评价；第二步、按照数据划分、融合策略进行洼地填平并行计算；第三步、按照数据划分、融合策略，基于洼地填平结果进行水流方向并行计算；第四步、按照数据划分、融合策略,基于水流方向数据，进行汇流累积并行计算；第五步、按照数据划分、融合策略,基于水流方向和汇流累积数据，设置汇流阈值，进行河网水系并行计算；第六步、按照数据划分、融合策略,基于河网水系和水流方向数据进行子流域划分并行计算，完成流域提取。本发明方法能够在并行时充分考虑数据划分的粒度及并行I/O机制，在串行算法并行解析上能够顾及数据自身特征。
文档编号G06F9/50GK102915227SQ201210320618
公开日2013年2月6日 申请日期2012年9月3日 优先权日2012年9月3日
发明者刘学军, 汤国安, 江岭, 宋效东, 祝士杰 申请人:南京师范大学