CN102902590B - 面向并行数字地形分析海量dem部署与调度方法 - Google Patents

Abstract

本发明公开了一种面向并行数字地形分析海量DEM部署与调度方法，其步骤如下：第一步是并行应用程序的初始化，各计算节点同时获取并行计算环境及DEM的各种参数；第二步、各计算节点根据本地数据块与缓冲区执行串行算法，在必要时对缓冲区的相应数据进行更新操作；第三步、调用并行计算粒度控制模块，对并行计算的数据粒度进行统计与更新，直到并行计算的效率达到最优；第四步、生成结果文件：对并行地形分析应用申请的各种内存进行回收，清空临时内存与无关的外存数据；更新日志文件，统计并记录各节点每次执行不同操作的时间，以便该地形分析算法的重复执行；进行结果融合与后操作，最终输出结果文件。本发明具有广阔的工业化应用前景。

Description

面向并行数字地形分析海量DEM部署与调度方法

技术领域

本发明涉及一种海量数字高程模型数据管理与调度的方法，具体说，是一种基于分布式并行数字地形分析技术与数据并行计算模型的动态数据分配部署方法。 

背景技术

数字地形分析是在数字高程模型（DEM）上进行地形属性计算和特征提取的数字信息处理技术。自上世纪60年代以来，数字地形分析已在制图、土木、水利、交通、水文、农业、规划等地学领域得到广泛应用。随着研究的深入和应用领域的扩展，各应用领域给数字地形分析提出了新的要求。同时，各种新型传感器，如LiDAR、InSAR、SAR等的出现，使DEM在数据获取上更加灵活和迅速，网络通讯技术使得各类DEM，如SRTM DEM、USGS DEM得以共享，特别是空间数据基础设施的兴起，区域级DEM的建库工作已相继完成，DEM呈现多比例尺、多分辨率的系列化的特征。面对如此浩瀚的、多元化的DEM数据，如何将其转化为能在社会经济过程中直接应用的地学信息，并将这些地学信息转化为地学知识，从而推动地学研究和应用的发展，已经成为国内外相关领域研究的主题。这也迫切要求解决在面对如此海量DEM数据时，现有分析方法及硬件资源在计算能力上严重不足甚至是无能为力的问题。并行计算技术为解决这一问题提供了新的思路。并行计算在空间信息科学领域中的应用由来已久，并取得了一些重要成果。但是，人们对DEM数字地形分析的并行处理模式还没有形成统一而有效的认识，没有完整的数字地形分析并行软件；如何运用统一的模式和方法移植和并行化既有的串行DEM数字地形分析算法，传统的分析算法是否都能够进行并行化，如何在并行的环境中实现并行与串行算法与处理模式共用与互操作，这都是摆在众多研究者面前迫切需要解决的问题。 

数字地形分析是地理信息系统（GIS）软件必备的功能之一，地形因子算法的并行化也是所有与地形相关的复杂地理计算必须考虑的基础要素。数字地形分析并行计算的有效实现是与地形相关的地理建模及复杂地理计算的基础与关键技术问题。目前，众多应用领域对DEM数据分辨率及数字地形分析效率的要求越来越高，相应的研究规模也不断增大，计算资源低利用率和海量数据高分析效率要求之间的矛盾日益突出，归纳起来有以下几点局限性： 

（1）数字地形分析并行计算要解决的关键技术问题包括算法的任务分解、数据的有效拆分及高效融合等。现有数字地形分析的并行处理模式多从单一视角出发，未从整体上对数据、任务、计算环境进行综合考量，也未见有全面探讨三者之间相互关系及量化模型的研究成果，这就使得研究者很难把握处理器负载均衡策略并达到并行效率最大化的目的，也在一定程度上制约了数字地形分析并行技术的发展和应用。 

（2）数字地形分析算法的特点是计算密集和数据密集兼而有之，已有并行计算环境支持下的研究，多集中在某一个方面处理方法，孤立了I/O密集和计算密集问题，并未对任务依赖关系、DEM数据块和CPU/GPU处理能力匹配做深入研究。现有的串行算法的并行化需要重点考虑算法依赖关系，由于地形数据的非结构化和计算密集特征，使得任务调度、数据并行平台的搭建等处理相对于结构化数据来说困难得多。 

（3）在目前基于“数据分块→并行计算→结果合并”的并行计算实施策略中，因为其数据调度模式和计算方法，在计算单元内和基于单一数据集的计算方法并没有本质差异，尚未形成高效、成熟、稳健的CPU/GPU任务分配机制，尤其是在计算过程中，还较少有针对动态数据分块的自适应任务、数据调度策略的研究，包括动态地任务加载、任务计算、任务数量预测等。因此，研究适合不同数字地形分析内容的自适应任务分配策略和动态管理机制将是一个新的挑战。 

发明内容

本发明针对以上缺陷，设计了一套新的并行数字高程模型部署与调度的工业技术方案。该技术流程涵盖了具有自适应计算节点性能数据分配技术和负载均衡可达的格网DEM数据调度技术两大关键技术环节，该技术方案能够解决现有并行地形分析算法中DEM低效率管理与调度的问题。 

本发明的基本思想是在并行程序初始化与并行程序执行的过程中，完成不同性质算法数据的统一部署与调度，实现格网DEM数据在并行程序的整个过程中有效地减少通信次数并能够满足各类地形分析算法的数据约数关系。在保证并行程序高效率执行的同时，通过最优并行粒度量化体系的相关指标，实现并行粒度的动态控制，从而挖掘并行应用的潜在计算效率。 

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

面向并行数字地形分析海量DEM部署与调度方法，其步骤如下： 

步骤1，初始化并行应用系统的信息； 

步骤2，各计算节点调用执行并行算法模块，根据本地数据块与缓冲区执行串行算法； 

步骤3，调用并行数字地形分析控制粒度模块，对并行计算的数据粒度进行统计与更新；

步骤4，生成结果文件：对并行地形分析应用申请的各种内存进行回收，清空临时内存与无关的外存数据；更新日志文件，统计并记录各节点每次执行不同操作的时间，以便该地形分析算法的重复执行；进行结果融合与后操作，最终输出结果文件。 

本发明的技术特点及有益效果： 

（1）本发明提出了一套通用的并行地形分析数据部署与调度策略，运用了数据并行策略对不同性质的数字地形分析算法进行并行计算，通过结合缓冲区机制实现了不同依赖程度的 算法的正确实行，采用了“数据并行为主，辅助以缓冲区”技术，实现了“不同并行算法，通用数据部署；不同依赖关系，动态数据调度；数据控制粒度，网络环境无关”的并行地形分析海量DEM的部署与动态调度，具有广阔的工业化应用前景。 

（2）本发明提出了最优并行控制粒度，通过记录每次程序执行的粒度大小与计算时间，对计算效率进行综合量化评价，提升了高性能地学分析的计算效率。该方法具有很好的可行性和稳定性，动态调度策略可望得较为理想的计算效率，有效地缓解了分布式并行计算效率低与通信耗时的瓶颈问题。 

（3）本发明提出的DEM并行部署模块与并行粒度控制模块可以独立使用，也可以集成使用。使用DEM并行部署模块，快速完成串行算法的并行化与数据调度。而通过最优并行粒度量化体系，可以计算出不同算法的最优并行计算粒度，充分保证海量DEM并行处理的高效率进行。二者具有较高的可移植性，为不同的并行环境提供了参考。在未来复杂多变的异构计算环境下，可用于计算的划分和负载均衡控制，从而提供比传统算法更高的任务执行效率和更好的可移植性。 

附图说明

图1是本发明的初始化工作流程图； 

图2是DEM数据划分方式，（a）一维行划分，（b）一维列划分，（c）二维块划分； 

图3是3乘3邻域型算法数据一维行划分缓冲区示意图； 

图4是二维块划分时不同子块的缓冲区示意图，其中（a）表示子块的上下左右方位都有缓冲区，（b）表示子块的右方位没有缓冲区； 

图5是本发明方法执行的主要流程图； 

图6是优化粒度算法执行主要流程； 

图7是结果融合示意图（进程间数据块存在数据依赖关系）； 

图8是生成结果文件的基本流程图； 

图9是本发明实施例中不同计算粒度下算法的执行时间； 

图10是本发明实施例中不同计算粒度下洼地填平算法的加速比。 

具体实施方式

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

本发明中设计和实现的过程中主要的类和接口如表1。 

表1主要类以及功能 

第一步、并行应用程序的初始化：本发明设计了针对各种形式的格网数字高程模型的并行应用程序初始化操作（见图1）。通过该初始化工作流程，各计算节点可以同时获取并行计算环境及DEM的各种参数。 

并行应用程序初始化的具体实现过程为： 

（1.a）进入MPI环境，各进程读取全局DEM和并行环境的各种信息，如：全局DEM左上角X，Y坐标、无值区的标识值（NOData值）、栅格点间距、栅格点数据类型、全局栅格的总行、列数、当前DEM的波段数、当前进程的进程号与并行环境总的进程数目等。 

（1.b）调用CMemControl模块，为了最大程度的保证物理内存的充分使用，初始内存限制设定为所有进程可用物理内存的平均值的十分之一，各节点可用物理内存的均值(FreeMem)是： 

$\mathrm{FreeMem}=\frac{\underset{\mathrm{rank}=0}{\overset{\mathrm{MPISize}}{\mathrm{\Σ}}}\mathrm{FreeMem}}{\mathrm{MPISize}}---\left(1\right)$

其中，MPISize是进程的总数目。然后，对初始内存限制——MAX_MEM_LIMIT进行初始化： 

MAX_MEM_LIMIT=0.1×FreeMem    （2） 

这里以MB为内存调度的基本单位。 

（1.c）选择数据块划分的类型： 

选择数据块划分类型后，各计算节点将执行同样的数据划分方式。这里的数据划分方式分别为：Row_DimenI是一维行划分，以一个栅格行为最小单元进行划分；Column_DimenI是一维列划分，以一个栅格列为最小单元进行划分；Block_DimenII是二维块划分，划分示意 图如图2所示。 

（1.d）根据内存限制、数据划分类型与栅格数据类型，各计算节点根据内存限制申请内存使用，然后动态地计算各进程第一次需要读取的栅格数据，各节点需要计算的数据块栅格点数LocalPixels为： 

$\mathrm{LocalPixels}=\frac{\mathrm{MAX}\_\mathrm{MEM}\_\mathrm{LIMIT}\×2^{20}}{\mathrm{PixelsBYTES}}---\left(3\right)$

其中，PixelsBYTES为栅格点数据类型的字节数。如果按照1维数据行划分，则本进程的行数、列数分别是： 

$\left\{\begin{array}{c}\mathrm{LocalRowNum}=\frac{\mathrm{TotalPixels}}{\mathrm{TotalCol}\×\mathrm{MPIRank}}\\ =\frac{\mathrm{TotalRowNum}}{\mathrm{MPIRank}}\\ \mathrm{LocalColNum}=\mathrm{TotalCol}\end{array}\right.---\left(4\right)$

其中，TotalPixels是总的栅格点数目，TotalCol是总列数，MPIRank是当前进程号，TotalRowNum是总行数。左上角X坐标与左上角Y坐标分别是： 

$\left\{\begin{array}{c}\mathrm{LocalXTopLeft}=\mathrm{XTopLeft}\\ \mathrm{LocalYTopLeft}=\mathrm{YTopLeft}-(\mathrm{MPIRank}-1)\×\mathrm{LocalRowNum}\×\mathrm{YCellSize}\end{array}\right.---\left(5\right)$

其中，XTopLeft是全局DEM的左上角X坐标（东西方向），YTopLeft是左上角Y坐标（南北方向），YCellSize是南北方向上相邻栅格点间距。 

更新本地的控制信息。调用LogFileManagement模块，将相关信息写入到日志文件。 

（1.e）根据地形分析算法对格网DEM的数据依赖程度，设定一定的数据缓冲区。地形分析算法可以分为邻域型算法和全局型算法，邻域型算法的并行化可以通过设定一定的缓冲区来计算本地数据块边界数据的相应地形因子；全局型算法可以通过若干的缓冲区动态地进行交互，实现全局信息的迭代更新。邻域型算法一维行划分的缓冲区示意图如图3所示，二维块划分缓冲区示意图如图4所示。 

（1.f）调用LogFileManagement模块，将初始化的相关参数、内存申请、计算的数据读取时间与数据划分类型等信息写入日志文件，并以“LOG_”+MPIRank+“_”+MPISize+“_”+DTAALGO命名，其中DTAALGO为本次执行的算法名称。该文件由各进程分别创建并存储在网络文件系统（NFS，Network File System）上，属于共享式全局日志文件。 

第二步、各计算节点调用DTAParaCom模块，根据本地数据块与缓冲区执行串行算法。如果算法类型属于全局类算法，则在程序的执行过程中根据算法对数据的不同依赖程度，调用 BufferShare模块动态的更新缓冲区数据，更新DTABuffer数据结构内的缓冲区数据，更新本节点对缓冲区存在数据依赖的计算结果，并行算法执行的主要流程图如图5所示。 

第三步、调用ControlGranuModule模块，对并行计算的数据粒度进行统计与更新。该模块只适用于全局DEM不能一次被所有的进程读取的情况下，各进程需要执行若干次地形分析算法，也即是对不同的数据块执行同样的地形分析算法。 

控制粒度更新的具体实现过程是： 

（3.a）读取各进程创建的日志文件，统计各进程的执行时间CompTime。 

（3.b）计算一次地形分析结束后的控制粒度计算效率： 

$\mathrm{CompEfficiency}=\frac{\mathrm{LocalPixels}}{\mathrm{Comp}{\mathrm{Time}}_{\mathrm{order}}}---\left(6\right)$

主节点记录该粒度计算效率，并写入日志文件。 

（3.c）从各进程执行第三次数据块的地形分析算法开始，对前两次的粒度计算效率进行比较。当一次子任务计算结束后，主节点计算控制粒度计算效率CompEfficiencyorder(简写为E_order)，并判断比较上一粒度效率，如果E_order>E_order-1，则继续增大计算粒度，即：D_orde+1=2×D_orde，这里的D_order代表第order次数据块计算时的栅格总点数；否则D_orde+1=D_orde；如果粒度效率停止增长，D_orde+1=D_orde×1.5，继续计算粒度计算效率。更新控制粒度，并写入日志文件。这里的控制粒度指的本次数据块的总的栅格数目。 

（3.d）如果总的DEM计算尚未结束，主进程向各进程分发新的计算粒度，各进程按照该粒度即总的栅格数目，执行第二步，如果数据全部计算结束，执行下一步骤，如图6所示。 

第四步、生成结果文件：对并行地形分析应用申请的各种内存进行回收，清空临时内存与无关的外存数据。更新日志文件，统计并记录各节点每次执行不同操作的时间，以便该地形分析算法的重复执行。进行结果融合与后操作，最终输出结果文件，数据融合的示意图如图7所示，生成结果文件的基本流程图如图8所示。 

下面以并行填洼算法为例。 

从数字高程模型中提取河网水系需要经过一系列的预处理，其中最复杂的部分是对DEM中洼地的处理。洼地聚合和嵌套式洼地合并过程被认为是洼地填平算法中计算较为复杂的一类操作。洼地涉及到全局，大小不固定。为了保证并行填洼中洼地的正确处理，各节点在判定边界线上是否存在洼地时可能需要访问更多的数据。当利用一定缓冲区无法判断边界线上是否有洼地时，计算节点需要把数据读取范围继续扩大或者交换相应的全局信息，直到处理所有洼地。当计算粒度较大并超过系统的最大内存限制时，程序将无法处理的数据转存在外存。 

该算法的DEM部署与动态调度过程为： 

第一步、并行应用程序的初始化 

（一）进入MPI并行环境，获取DEM的各种信息，存入相应的数据结构。 

（二）调用CMemControl模块，计算初始内存限制。 

（三）选择数据块划分的类型为一维行划分，转换为行数100行。 

（四）计算各节点需要读取的栅格点总数，并行读取栅格数据。 

（五）设定缓冲区——1行，将相应的数据存入到DTABuffer数据结构。 

（六）将相应的计算、内存申请等情况写入日志文件。 

第二步、调用DTAParaCom模块，对各数据块执行串行算法。各节点在处理完一次串行算法后，更新缓冲区内部数据，并对数据块的结果数据进行相应的更新操作。 

第三步、调用ControlGranuModule模块，对不同进程的数据块进行控制粒度评价，对并行粒度进行相应的控制。如果有可更新的控制粒度，下一次数据块大小将更新，并将相应的信息记录到日志文件。 

第四步、生成结果文件，清空临时内存与无关的外存数据，更新日志文件，统计并记录各节点每次执行不同操作的时间，然后进行结果融合与后操作，最终输出结果文件。 

基于以上并行计算流程与动态调度机制，示例算法并行计算优化的效率如图9、10所示。该发明的方法增加了量化的机制来深入讨论数据并行问题，以及同样数据量时不同并行粒度对计算效率的影响，避免处理海量数据时计算时间、效率预测不准确问题。与传统的并行算法相比，经过并行控制粒度的动态调度后，本方法具有很好的可行性和稳定性，动态调度策略可望取得较为理想的计算效率，有效地缓解了分布式并行计算效率低的问题。此方法有待于在更多的领域内应用，以检验其性能。 

Claims (1)

1.面向并行数字地形分析海量DEM部署与调度方法，其步骤如下：

步骤1，初始化并行应用系统的信息；

步骤2，各计算节点调用执行并行算法模块，根据本地数据块与缓冲区执行串行算法；

步骤3，调用并行数字地形分析控制粒度模块，对并行计算的数据粒度进行统计与更新；

步骤4，生成结果文件：对并行地形分析应用申请的各种内存进行回收，清空临时内存与无关的外存数据；更新日志文件，统计并记录各节点每次执行不同操作的时间；进行结果融合与后操作，最终输入结果文件；

所述步骤1的具体过程如下：

步骤11，在MPI环境，各进程读取全局DEM和并行环境的各种信息；

步骤12，调用内存管理模块，为了最大程度的保证物理内存的充分使用，初始内存限制设定为所有进程可用物理内存的平均值的十分之一；

步骤13，选择数据块划分的类型；

步骤14，根据内存限制、数据块划分类型与栅格数据类型，各计算节点根据内存限制申请内存使用，然后动态地计算各进程第一次需要读取的栅格数据，计算各节点需要的数据块栅格点数；

步骤15，根据地形分析算法对格网DEM的数据依赖程度，设定一定的数据缓冲区；地形分析算法采用邻域型算法或者全局型算法，其中，邻域型算法的并行化是通过设定一定的缓冲区来计算本地数据块边界数据的相应地形因子；全局型算法是通过若干的缓冲区动态地进行交互，实现全局信息的迭代更新；

步骤16，调用日志文件管理模块，将初始化的相关参数、内存申请、计算的数据读取时间与数据划分类型信息写入日志文件；

所述步骤3的具体过程如下：

步骤31，读取各进程创建的日志文件，统计各进程的执行时间；

步骤32，计算一次地形分析结束后的控制粒度计算效率；

步骤33，从各进程执行第三次数据块的地形分析算法开始，对前两次的粒度计算效率进行比较；

步骤34，如果总的DEM计算尚未结束，主进程向各进程分发新的计算粒度，各进程按照该粒度即总的栅格数目，执行所述步骤2，如果数据全部计算结束，执行所述步骤4。