一种基于Linux集群平台的DEM并行等高线生成方法与流程

本发明属于高性能地学计算领域，具体涉及一种基于Linux集群平台的DEM并行等高线生成方法。

背景技术：

众所周知，自从1958年Miller在“Photogrammeric Engineering”杂志上发表题为“The Digital terrain model:theory and application”的论文中首次提出了数字地形表达的概念后，DTM、DEM的研究在随后的几十年里得到深入研究，已经广泛的应用在诸如测绘与遥感、农业规划、土木工程、军事、地学分析等领域。而等高线是地面高程相等的各相邻点所连成的闭合曲线，用于描述地球表面地貌情况。作为表面形态的一种主要表现手段，等高线不但能简单正确地显示地貌的形状，还可以直观和有效地表示空间连续分布的现象，精确反映在垂直和水平方向的强弱差异，具有较好的视觉效果，并能方便地进行图上量测，这些特点使其成为专题地图的重要表示方法，成为GIS中地理过程模拟与仿真的重要表现手段，在诸如地形模拟、水文模型建立、公路铁路勘测规划、自然资源管理、工程分析、计算领域、军事领域有着广泛应用。

目前，获取等高线的手段多种多样，其中通过地形数据自动生成等高线，如通过纸质地图扫描矢量化生成等高线，DEM栅格数据转换成生成等高线等多种途径和方法。其中利用DEM栅格数据通过算法自动生成等高线是获取等高线的最主要手段，该途径既是利用DEM进行地形分析的重要部分，也是GIS自动制图的一个重要研究内容。对于大多数GIS应用程序而言，数字高程模型通常通过栅格数据结构提供，因此DEM这一术语用来描述数字高程模型已经得到广泛接受。虽然推导坡度坡向的算法已经十分成熟，但是DEM数据会影响推导的参数的准确性。许多研究表明坡度坡向精度分析算法与DEM数据错误、数据精度、栅格分辨率、栅格大小或者网格方向相关。

数字高程模型(DEM)是高程点的集合，可以用来描述地理特征，展示三维地形表面的特定区域。栅格数据格式是DEM最常用的数据格式。随着DEM网格分辨率的提高，数字地形分析是GIS中的计算密集型任务之一。然而，某些特定的应用程序要求GIS软件为存储、检索和处理地理空间数据等操作提供高性能操作，在传统的GIS中，这些需求已经逐渐成为计算效率的瓶颈。

随着HPC的快速发展，并行计算被应用于改善GIS应用程序性能。在过去的几十年中，涌现了许多基于异构网络如何使GIS应用并行处理的研究。2003年，Hawick等人提出了处理大规模地理数据的并行算法。

技术实现要素：

针对上述存在问题或不足，为解决现有等高线生成方法中，处理DEM数据时十分耗时的问题；本发明提供了一种基于Linux集群平台的DEM并行等高线生成方法。

具体包含以下步骤：

步骤1、首先进行环境初始化，然后Linux各个节点借助PVFS并行文件系统并发读入输入栅格数据；

步骤2、GDAL(Geospatial Data Abstraction Library)读取输入的栅格数据，获取相应的输入参数及进行校验，然后初始化临时Shapefile结果文件；

步骤3、计算出每个计算节点应该处理的等高线区间：

即在进程号大于等高线条数与进程数目的模(nHeavryContour＝nSumContour MOD numprocs)的情况下获取平均处理的条数(nAverageContour)，否则获取(nAverageContour+1)即平均处理的条数加1；

步骤4、各个计算节点同时进行下述操作：

获取等高线区间内栅格图像中各点的高程值并形成矩阵，然后获取等高线的条数及具体高程值，当高程矩阵有高程值与要生成的等高线等级值相等时，偏离微小量，便于进行(level-Z0)(level-Z1)<0；

步骤5、各个计算节点各自生成各自产生的等高线数据，然后将这些数据写入对应的临时Shapefile文件里面；

步骤6、主节点接收来自各个计算节点的临时Shapefile文件，并将所有临时Shapefile文件合并，生成等高线。

进一步的，还包括一个后续的步骤7、在上述栅格DEM生成等高线并行算法的基础上，利用精灵进程进行离线的临时的结果数据合成方法对步骤6进行优化及测试分析；或/和利用小热点并行化方法对步骤4的微调过程(偏离微小量)进行优化及测试分析。

综上所述，本发明采用Linux集群平台对DEM等高线生成算法进行并行可以取得可靠的加速比，能大幅度减少运行时间，具有较多的模式选择，性能提升稳定可靠。

附图说明

图1栅格格网设置；

图2栅格矩形四边均有等高线穿过的情况；

图3 GRASS中r.contour模块整体实现算法流程；

图4 GRASS中r.contour模块cont函数实现结构图；

图5改进的栅格DEM生成等高线算法实现原理图；

图6改进的栅格DEM生成等高线算法运行参数及平台硬件配置信息；

图7改进的栅格DEM生成等高线算法热点分析(一)；

图8改进的栅格DEM生成等高线算法热点分析(二)；

图9并行算法设计方案；

图10 MPI实现算法部分程序；

图11 MPI并行优化算法；

图12 MPI并行优化算法实现过程；

图13优化方法二部分代码。

具体实施方式

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

一、栅格DEM生成等高线算法原理

目前，计算机自动绘制等高线的方法主要有两类，即网格法和三角形法。它们生成等高线的步骤一致，主要包括：

1)在格网边上内插等高点；

2)追踪等值点，形成某一高程值的若干等值线序列；

3)光滑等高线并输出。

目前，应用得最为广泛的栅格生成等高线算法采用的是等值线传播算法，隶属于网格法这一类，该算法的效率为O(N)。等值线传播算法的基本原理是：为了得到某一高程值的等高线，首先找到一个有对应等高线穿过的格网作为初始格网；然后由这个初始格网开始向邻近的栅格(上、下、左、右)传播这条等高线；当某一条等高线穿过格网时，可通过线性内插方法求出等高点；同时利用下述公理判断其下一个传播方向：

公理1：单元格网的边与等高线交点的个数为0，2或4；

公理2：已知等高线与单元格网的边有4个交点，则在单元格网内中心点和两个顶点同号的对角线与等高线没有交点。

直至所有栅格单元遍历一遍，再利用相同方法进行下一条等高线的获取。一个完整的可实现的栅格数据生成等高线算法的描述如下：

3-1、标志栅格图像的栅格。

设置一个标志，用来说明该栅格是否需要再被搜索，初始值与说明不再对该栅格进行搜索的标志值不同，被标志后直至等高线高程发生变化后，标志才会失效。

3-2、构成栅格矩形。

将待搜索的栅格与其右，右下，下方各栅格一起组成栅格矩形。如图3.3所示。

3-3、判断栅格点间等高线是否穿过。

如图1所示，以点0和1为例，设点0高程为Z₀，点1高程为Z₁，判断它们之间是否有高程值为level的等高线穿过，可根据是否满足式(1-1)进行判断。如果满足(1-1)所列条件，则说明点0和点1之间有等高线穿过；反之则否。

(level-Z₀)(level-Z₁)<0 (1-1)

3-4、插入等高线矢量点。

对于每条等高线，依次搜索栅格图像的各个栅格，对每一栅格，判断该点是否出界及是否已被标记不再搜索。如果均否，则判断相邻两栅格之间是否有等高线穿过(如图1的点01之间，12之间，23之间及30之间)。若有，则按线性内插公式1-2计算并存入对应数组。

线性插入法使用的条件是平面上点的插入，在我们要处理的栅格图像中，四个角点可看作在一个平面内，因此可用线性插入法。

3-5、寻找该栅格可继续搜索的方向。

通过上述提及的公理可知，寻找当前搜索栅格可继续搜索的方向有两种情况：

a)栅格矩形只有两边有等高线穿过的情况。

在这种情况下，上一步中我们已经插入了其中一个点(可称之为第一个交点)，而该栅格部分的等高线必然是从第一个交点指向第二个交点，因此只需将第二个交点插入即可，将栅格矩形移到以这点所在边为边的新矩形上。

b)栅格矩形四边均有等高线穿过的情况。

在这种情况下，我们已经插入了第一个交点，先求四个角点的高程平均值，设为mid(如图2)，那么在栅格矩形中必然存在一个点，其高程与mid对应，再分别按照3-3中方法判断mid与栅格矩形四个角点之间是否有等高线穿过，即可得到该栅格部分的等高线的后续搜索方向，将栅格矩形移到以这点所在边为边的新矩形上。

3-6、搜索完毕进入下一栅格的搜索。

过程重复3-2至3-5步，直至该栅格图像到底或已无栅格可再搜索，此时将进入下一条等高线的搜索。

二、优化后的等高线追踪串行算法实现

当等高线结束在栅格的边缘(如上下左右)而呈现出首尾不相连的形状时，我们称之为开曲线；而等高线在栅格单元内部形成首尾相连的形状时，则为闭曲线。上述描述的等高线追踪算法是没有区分生成的等高线是开曲线或者闭曲线的，而在实际算法中，是把开曲线和闭曲线分开考虑，以有效地提高等高线追踪算法的效率。本发明中的栅格数据生成等高线算法是采用GRASS GIS中的对应的栅格数据生成等高线模块r.contour所采用的等高线追踪算法，它实际上是将等高线追踪过程分为开曲线和闭曲线两部分。

在GRASS中，由栅格图像生成等高线算法即是将一幅栅格图像按等高线条数或等高线步长，提取高程相等的点，生成该图像对应的等高线矢量图。算法模块主体与上述等值线传播算法基本原理描述类似，只是在在该算法处理之前，需要获取生成等高线的序列，然后依次循环所有不同高程值依循上述算法原理生成对应的等高线。

其算法实现的主要过程如下：

1)生成每一条等高线，需先设立栅格标志数组，初始值为0，当处理过并已记录的栅格点均标志为1，以后不再处理，直至下一条等高线开始；

2)对于每一条等高线，依次从栅格图像的上部(前3行)、下部(最下面3行)、左部(最左边3列)、右部(最右边3列)、中间部(除去上下左右3行部分)的第一个栅格开始搜索，对每一个栅格，判断该栅格点是否出界以及是否已被标记为不再搜索；当需要处理时，则按上下、左右、中间的顺序依次判断相邻两栅格之间是否有等高线穿过；如果有，则按线性插入法插入，并记录可继续搜索的方向；

3)继续其他可搜索的方向。直至该栅格已被标记为不再搜索或出界或搜索的当前栅格的相邻栅格之间无等高线穿过，则该栅格搜索结束。同时需要将插入的矢量点予以记录；

4)进入下一个栅格的搜索，过程重复步骤2、步骤3，直至所有栅格均被搜索；

5)开始下一条等高线的生成。

在GRASS GIS中，栅格数据生成等高线算法的对应模块为r.contour，算法整体实现流程如图3所示。

其中，利用获取的高程值数据(z_array[])以及经过处理的等高线序列(lev[])生成具体的等高线数据(分为开曲线、闭曲线)是整个模块的主要部分，由cont()函数完成。该函数的具体实现流程如图4所示。

在后续的研究中，我们实际采用的栅格数据生成等高线算法模块是基于GRASS GIS的等高线生成模块r.contour的，其算法基本原理一致，但为了保持模块的可独立运行，在输入/输出数据格式上有着较大的区别。

因为：

1)该算法模块直接基于GRASS GIS基础模块，需要运行grass 64初始化模块后才能运行，即它与GRASS GIS有非常紧密的依赖关系。如果采用这个模块，将使得实验程序变得非常庞大，而直接从里面剥离独立算法也非常复杂。

2)由于该算法的输入/输出为GRASS GIS内部定义的栅格和矢量数据格式文件，故如果不采用GRASS GIS依赖包，那么需要对该算法的输入/输出重新进行设计。

改进的、能独立运行的栅格数据生成等高线算法的流程图见图5。

三、串行等高线生成算法热点分析

改进后的串行等高线生成算法名为contourgen，运行参数以及进行热点分析的平台硬件配置如图6所示，其中为了说明该算法在不同计算规模下情况，通过调整生成等高线的步长(参数s)来进行的，步长越短，生成的等高线愈密计算规模愈大。

图7，8分别为该串行程序在Linux平台上的热点测试信息，从中可以看出：不管计算规模如何，其中热点部分(contour部分，即生成等高线部分)占了整个运行时间的89.2％左右(在四个实验中依次为86.2％，89.1％，90.0％，91.5％)。其中另外一个耗时的子过程是Dismatrix，即进行栅格单元与要生成等高线高程值的微调过程(如果栅格单元高程值与要生成的某等高线相同，需要加上一个非常小的增量，便于进行(level-Z0)(level-Z1)<0，该热点在四次中依次对应为11.3％，10.3％，9.7％，8.3％。

从上可知，不管数据大小如何，计算规模如何，热点依然集中在生成等高线这部(contour)，那么为了减少算法运算时间，需要从这部分进行入手进行并行化。

四、栅格DEM生成等高线并行算法设计与实现

基于已有栅格数据生成等高线串行算法以及常见MPI并行化方案，本发明选择了一种并行化方案具体实现。

1)MPI并行化方案

该算法采用以下并行方案：每个节点均可以独立读取完整的DEM输入数据，主节点分解不同等级(如等高线为100,150,200米…)的等高线任务给各个节点；每个节点各自生成不同等级的等高线，然后同时生成结果文件输出。

该算法设定的输出格式为shapefile，本发明采用的方式是：各个进程独立生成临时的shapefile文件，然后由主节点将这些临时的shapefile文件合并的方式进行。具体实现示意图如图9所示。

2)MPI并行算法实现

其中，需要重点关注以下几个主要问题：

a)PVFS并行文件系统的读写问题；PVFS并行文件系统的读写与传统文件读写调用函数一致，它在每个节点都有输入文件的备份，可以有效地避免传统方式多个节点读取同一个文件发生冲突的现象发生。

b)计算各个计算节点处理的等高线区间问题；对每个节点计算的等高线区间，可采用下述算法实现均分(图10)。

采用上述算法的结果时，在进程号大于等高线条数与进程数目的模(nHeavryContour＝nSumContour MOD numprocs)的情况下获取平均处理的条数(nAverageContour)，否则获取(nAverageContour+1)。

c)不同计算节点生成临时SHAPEFILE文件以及主节点合成最后生成文件问题。每个节点各自生成各自产生的等高线数据，然后将这些数据写入对应的临时SHAPEFILE文件里面；主节点在合成之前，需要先在生成SHAPEFILE文件前去探测各个临时SHAPEFILE文件实体对象的个数，然后按照总的实体个数进行初始化。

3)MPI并行算法性能优化方法

通过实验分析，该并行算法可以很好解决多个计算节点同时读取输入数据的问题，但注意到一个现实存在的问题：在将各个不同节点生成临时Shapefile文件合成最后的输出文件时，都是由主节点完成。如果参与的节点很多，而且生成的临时Shapefile文件有很大的情况，这部分的耗时也将成为该MPI并行算法的一个瓶颈所在。

调研已有文献发现：可以利用独立分布式数据库解决这个问题。但是由于该算法输出为传统的Shapefile文件格式，需要打开文件，生成对应的文件头，然后写入数据。这一系列操作皆由Shapefile库提供，非普通的读写操作。正是因为这个原因，采用分布式数据库并不是最好的解决办法。本发明采用将最后Shapefile合并操作采用离线方式完成。

另一方面，如果MPI集群中的每个节点CPU都是多核的，那么还可以通过MPI+OpenMP混合编程模型来提高该并行算法的性能。即将不同的等高线序列分配给不同的计算节点进行处理，在每个独立的计算节点上，某部分等高线序列处理时又采用多核进行多线程处理，从而细化了计算力度，充分利用了硬件潜能。这种混合编程方式尤其对于计算节点＝NP值，而不是核心数对应NP值的情形尤其适用。

3-1)优化1：利用精灵进程进行离线的临时的结果数据合成方法

采用离线合并Shapefile文件方式，最简单的一种方法是：将所有节点的生成的临时Shapefile文件合并为一个完整的Shapefile文件的操作单独抽取出来，做成一个可执行程序的方式由主节点最后以离线方式调用，而非由它执行完后才结束整个程序的方式，如图11所示。

为了让调用的可执行程序(这里指Shapefile合并程序)以离线方式运行，在MPI主进程调用的过程中采用如下形式：

mpirun-np 4./contourgen-i../../data/srtm_41_14/srtm_41_14.gif-o/home/fhuang/contourgen/contourGen-opt2-deamon/output/result-s 5

然而在实际运行中，该合并程序虽然是离线运行，但仍然是在返回子进程结果后，即程序运行完毕后MPI主进程才结束运行。很显然，采用这种方式不能达到预期目的。

后来通过调研，可以采用精灵进程的方式实现该目的。即精灵进程常驻内存运行，当满足设定条件(通过写入文件变量变化或socket编程方式)触发，将指定目录的文件调用Shapefile库函数进行合并。而MPI主进程则在合并前已经终止完成。

考虑到实现的简便性，本实施例采用的是通过触发某一特定文件来驱动精灵进程进行合并的方法。具体技术实现过程如图12所示。

值得指出的是，精灵进程需要先于MPI程序运行前启动。而且为了实现MPI并行程序与离线Shapefile合并程序的异步运行，即MPI并行程序在调用合并程序后自动结束，而合并程序仍在运行，因此需要等待30秒后再对结果数据进行操作。

3-2)优化2：小热点并行化方法

从MPI基本算法的各步耗时对比结果可知，当数据量、计算规模达到一定程度时，第二步操作——对获取的Z数组依次与要生成的等高线序列进行比对，如果有高程值与任一将要生成的等高线数组中值相等时，需要将Z数组中的高程值做一个微小的改动。因为在判断是否有要生成的等高线经过该像元时，是通过公式(level-Z0)(level-Z1)<0进行计算的；如果不做一个微小的调整，将会出现等于0的情况。如果出现这种情况，则直接有可能导致公式(1)出现分母为0的情况出现，程序将出现崩溃。故在生成等高线之前，需要进行如下处理(附图13)：

这部分代码函数对应DisplaceMatrix，从VTune分析结果可知，是这个算法的小热点(平均值9.9％)。该部分在数据量不大(影响Z数组大小)以及等高线间距较密(影响计算规模，直接增加这部分的操作次数)时，这部分的性能将成为整个并行算法的新的瓶颈。

通过分析可知，对这个小热点问题可以采用如下方法进行并行化。

a)不同进程只扫描Z数组中的一部分区域，记录下需要微调的行列信息；

b)每个进程通过MPI_Allreduce()函数获取整个Z数组需要微调的总个数；

c)每个进程分别开辟对应空间，通过MPI_Allgather()和MPI_Allgatherv()函数获取整个Z数组需要微调的全部位置信息；

d)每个进程统一对各自的Z数组进行微调。

采用这种方法，理论上可以将时间缩短为原来的1/N，当然一个像元值出现等于生成等高线数组的概率以及最后需要收集改动信息(需要微调位置的数组)大小有关，因为如果需要在各个进程间收集的数据量非常大，将影响这部分的并行效果。

3-3)终极优化：优化方法3-1+优化方法3-2

对于该算法来说，就是将前面的两种优化综合起来，即充分将整个MPI算法进行整体优化，得到性能的整体大幅度提升。