基于DEM的分布式并行通视分析算法研究

张刚，汤国安＊，宋效东，杨坤

（1.南京师范大学虚拟地理环境教育部重点实验室，江苏 南京 210023；2.南京师范大学计算机科学与技术学院，江苏 南京 210023）

0 引言

地形通视性分析（Terrain Inter－Visibility Analysis）是指基于DEM数据判断地形上任意两点之间是否可见的技术方法，通视的条件取决于视点与目标点间是否存在妨碍视线的障碍物。通视分析实质上属于对地形进行最优化处理的范畴［1］。目前，基于规则格网DEM的地形通视性分析已经广泛应用于通信、军事、房地产、考古、景观设计等多个领域。但随着空间数据采集技术的迅速发展，海量、高精度DEM数据的应用越来越广泛，通视分析算法的计算量也呈现出指数级增长趋势，传统的计算机处理技术已经不能有效地提高基于海量DEM数据的通视分析执行效率［2］。因此，如何利用现有计算资源，提高计算效率，降低算法时间复杂度成为通视分析要解决的关键问题。

近年来，随着海量空间数据分析的不断增长以及高性能地学应用的不断推动，并行空间分析方法成为高性能地学计算的发展趋势［3］。并行空间分析方法是将GIS空间分析方法与并行计算技术相融合的过程，旨在通过多种计算资源解决复杂的空间分析问题，提高空间数据处理的速度和质量。通视分析作为地形分析的重要组成部分，也是空间分析不可或缺的内容［4］。研究并行通视分析算法，成为提高海量数据下通视分析效率的有效途径。

针对通视分析的高效计算问题，国内外学者做了一定研究。Floriani等［5］在MIMD架构基础上，通过静态和动态两种并行数据划分策略，提出了一种基于TIN数据的并行可视性分析算法，但是由于TIN数据存储结构的复杂性，算法执行效率有待进一步提高，使用范围有一定的局限性。Kidner等［6］提出一种基于数据并行的反向可视性分析算法，各处理节点同步执行串行通视算法，获取由面到区域的可视域，该算法有效利用了各计算节点的计算能力并取得了较优的加速比，但是计算精度相对不高。由于高分辨率海量DEM数据的点对点通视分析算法计算量巨大，Mills等［7］针对LOS视线设计了数据划分策略，对传统点点通视分析算法进行了并行化处理，由于测试环境较为简单，实验结果有待进一步验证。在此基础上，本文从负载均衡角度，通过分析DEM数据存储的结构特征，提出了一种面向分布并行环境的数据划分策略，构建了并行Bresenham地形通视分析算法，并以全国90mSRTM为数据源，分析了并行通视分析的运行效率，为并行环境下的地形通视分析提供了新的思路。

1 并行通视分析算法设计框架

1.1 数据划分策略

并行计算是指同时使用多种计算资源解决复杂问题的过程，其包括共享存储模型、消息传递模型和数据并行模型3类［8］。由于栅格数据空间分析具有计算数据密集且数据组织规律性强的特点［9］，大部分栅格数据空间分析的并行模式均采用数据并行处理方式［10］。数据并行是将串行任务需要处理的原始大数据集划分为从节点处理的子数据集，各处理节点独立执行串行程序，获取计算结果。数据并行策略的基础是数据划分，针对DEM数据的矩阵式特征，本文从数据并行的角度对基于DEM的分布式并行通视算法进行分析。DEM数据划分主要包括规则划分和不规则划分两种，其中规则划分又分为一维划分和二维划分（图1）。从数据结构角度，DEM数据是一个m行×n列的矩阵，一般通过数组的形式进行组织。因此，规则划分相较于不规则划分更易实现且有利于并行计算系统的负载平衡。实验表明，在规则划分中，一维数据划分程序的执行效率优于二维划分，而DEM数据在计算机物理内存中是按行存储的，故本文选择一维行划分作为其数据划分方法：将数据集按行平均划分为n等份于n个处理节点，如果存在剩余行，则分配给最后一个处理节点。

图1 DEM数据划分方法Fig.1 Regular partition of DEM

针对基于DEM的通视分析算法，其有效数据区为视点与目标点所对应的外接矩形区。如果在数据处理过程中将整个DEM数据读取后再进行并行通视分析操作，会大大增加程序负载量及计算冗余，降低程序的执行效率。因此，在对DEM进行数据划分操作前，首先完成有效数据区的裁切处理（图2）。在数据裁剪过程中，通过视点与目标点位置计算外接矩形及其左上角点的绝对坐标，更新DEM数据头文件信息并获取有效数据区，合并生成新的DEM数据，为数据划分提供数据源。

图2 数据裁切过程Fig.2 Process of extracting DEM

在数据并行过程中，主节点将数据划分后的子数据集发送给从节点，从节点独立完成数据处理，但针对海量DEM数据，每个从节点所能容纳的最大计算量是有限的，为了保证每个计算节点处理的数据不超过其本身所能容纳的最大计算量，本文还设计了内存限制参数，以保证每个计算节点的最大计算效率，以获取最优加速比。根据数据量大小、计算节点数目以及内存限制，可以得到每个处理节点循环计算的次数：

在理想状态下，数据划分后计算节点间的处理过程是相互独立的，不需要进行通信交互。然而，地形通视分析算法需要实时判断每个子数据集中目标点的可视性，以便及时结束运算，提高计算效率，因此，子数据结果集之间需要实时通信，即结果数据通信。

1.2 并行通视处理

1.2.1 并行通视分析原理 地形通视性分析主要是研究视点与目标点是否可视的问题，具有简单复杂性、不可逆性以及可视不变性三大特征［4］。由于目前规则格网DEM的广泛应用，基于DEM的地形通视性算法众多［11－15］，其基本原理为：做从视点V到目标点T的射线，判断视线所经过的地面点高程与相应视线上点的高程之间的关系，如果在视点与目标点之间，有任意一点的高程高于射线高程，则两点不通视，否则通视。基于DEM通视分析算法基本原理大致相同，不同之处主要集中在高程内插方法和可视判断原则方面，其中，最为著名的是Bresenham算法，该算法基于整数运算，将栅格单元看做是同质的，获取栅格高程的内层循环紧密而且高效，避免了其他通视算法的插值运算过程。因此，本文选择Bresenham算法作为并行通视分析的核心算法。

基于DEM的通视分析并行算法采用数据并行的处理方式，按照一维行数据划分策略，根据内存限制，将更新后的DEM的数据划分为不同的子数据集合，各个计算机处理节点分别读取子数据集数据后，同步执行Bresenham通视分析算法，并将计算结果发送到主进程，主进程根据各个进程的返回结果判断视点与目标点之间的可视情况。如果子进程计算结果存在不可视情况，则视点与目标点之间不可视；否则两点可视。并行通视分析原理如图3所示。

图3 并行通视分析原理Fig.3 Principle of parallel visibility analysis

1.2.2 Bresenham通视分析算法 假设直线位于第一象限，从起点（x1，y1）到终点（x2，y2）。直线方程可以表示为y＝kx＋b（k∈［0，1］）。式中：

如图4所示，当直线y＝kx＋b沿x主轴增加1个单位，即xi＋1＝xi＋1时，在数学意义上其y副轴也会相应的增加斜率大小k。由于计算机以栅格化的方式进行图像表达，实际点位在绘制过程中都会存在一定的误差w（w∈［－0.5，0.5］），此时，yi＋1的真值为yi＋1＝yi＋k＋w，而点yi＋1在计算机屏幕进行绘制时，只可能选择yi＋1＝yi＋1或者yi＋1＝yi，选择的原则为k＋w遵从四舍五入，分段函数表达式为：

yi＋1绘制过程中造成的误差记为wnew，则：

在程序执行过程中，为消除浮点数运算，将上述误差公式两边同乘以dx，并令φ＝dx＊w，则：

图4 Bresenham算法原理Fig.4 Principle of Bresenham algorithm

通过Bresenham算法可以得到视线经过的所有栅格单元的行列号，利用DEM及视线所经过的栅格单元行列号获取地表高程TerrainElev；如果当前点高程值为NoDataValue，即空值，则默认此点可视，继续计算下一点，否则，通过比较视线上当前点高程ViewElev与TerrainElev的大小判断是否可视。

综上，Bresenham通视算法思想为：1）计算dx＝x2－x1，dy＝y2－y1，计算误差初始值φ＝0；2）求取直线下一点位置xi＋1＝xi＋1，若φ＋dy≥dx，yi＋1＝yi＋1，否则yi＋1＝yi；3）计算点（xi＋1，yi＋1）的地表高程TerrainElev及视线高程ViewElev，比较两者大小：若ViewElev＞TerrainElev，则视点与目标点不可视，否者执行步骤4；4）更新误差值φ，若φ＋dy≥dx，φ＝φ＋dy－dx，否则φ＝φ＋dy，转步骤2，直至点（x2，y2），程序结束。

2 实验结果分析

笔者在Microsoft Visual Studio 2010环境下，利用C⧺、GDAL1.6开源库及MPICH2并行计算平台实现了并行通视分析算法。实验的硬件环境为分布式并行环境，其主节点服务器CPU为Intel Xeon E5645，主频2.4GHZ，内存32G，核心数24个；从节点服务器CPU为Xeon E5620，主频2.4GHZ，内存16G，核心数8个，集群环境共包含20个从节点，主节点与从节点之间通过千兆网进行连接。采用全国90m分辨率SRTM陆地表面地形数据作为数据源，其数据大小为13.2GB（74 869列×46 784行），从中随机选取3条视线（长度分别为3 836km、4 473km、4 283km）分析并行通视分析代码效率，如图5、图6所示。

在实验过程中，为了对通视分析并行执行效率进行有效的评估，提出了性能加速比的概念模型。设N为实验过程中总进程数目，P为当前进程数目，T（N，P）为当进程数为P时通视分析并行执行时间，T（N，1）为通视分析在进程数为1时的执行时间，也就是串行执行时间。定义Radio为通视分析并行算法执行加速比，则：

图5 并行通视算法性能测试（通视情况）Fig.5 Test of parallelism of inter－visibility analysis algorithm （visible）

图6 并行通视算法性能测试（不通视情况）Fig.6 Test of parallelism of inter－visibility analysis algorithm （invisible）

从通视分析并行化算法性能测试图中可以看到：1）不论视点与目标点之间是否通视，随着处理节点数的增加，各类型并行化算法执行时间的降低速率变小，直至趋向于平缓。处理节点数小于32时，执行时间变化突出，变化速率较高，而当处理节点数大于32时，逐渐趋向于同一水平状态，执行时间变化缓慢甚至存在效率降低的趋势。因此，基于海量地形数据进行分布式并行通视分析时，其计算效率与进程数具有一定的关系。2）地形通视分析并行化后，不通视情况下最大加速比达到7，通视情况下最大加速比达到4，计算效率明显高于单个处理节点，即高于串行算法执行效率。在通视与不通视状况下，并行通视算法加速比存在不一致的情况，这是因为算法的并行性能在一定程度上受到地形数据的影响。从研究区域的地形复杂性角度分析，当研究区域地形起伏较大时，不通视状况下的计算效率体现得尤为突出；从程序设计的角度分析，在程序设计过程中采用了单一进程终止策略，通视分析并行算法在运行过程中，每一个处理节点都会执行相同的通视算法，当任意一个节点获得不通视的结果后程序就会自动结束此节点的运行过程，并将不通视的结果通过MPI通信机制发送给主节点，由于MPI函数库不存在同时终止所有处理节点运行的机制，因此当所有子节点将返回结果发送到主进程后，主进程对其从进程的发送结果进行分析，只要有一个进程的返回结果是不通视，则视点与目标点之间即为不通视，否则两点之间通视。通过这种单一进程终止策略可有效地提高并行算法的执行效率。从图中可以看出，随着处理节点数目的增加，程序的执行效率存在一定的波动，这主要是由于代码执行效率比较高，进程间通信时间不稳定直接影响着并行程序的执行加速比。

3 结论

本文针对基于DEM的分布式并行通视分析算法进行了分析研究，提出了一种通用的有效数据可达的DEM数据划分策略，实现了Bresenham通视分析并行算法设计，经验证，本研究所设计实现的并行算法在计算时间以及加速比方面能够体现出较好的优势，大大提高了海量数据下通视分析的执行效率，为并行环境下的地形通视分析提供了新的思路。但是，如何进一步提高并行I／O的读取效率以及减少进程间通信时间，进一步提高并行通视分析算法加速比，还需要进一步研究。

［1］ 李志林，朱庆.数字高程模型［M］.武汉：武汉大学出版社，2001.

［2］ 吕品，张金芳，鲁敏，等.基于视线的地形可视性分析研究［J］.计算机工程与应用，2006（8）：223－226.

［3］ 宋效东，刘学军，汤国安，等.DEM与地形分析的并行计算［J］.地理与地理信息科学，2012（4）：1－7.

［4］ 周启鸣，刘学军.数字地形分析［M］.北京：科学出版社，2006.249－266.

［5］ FLORIANI L D，MONTANI C，SCOPIGNO R.Parallelizing visibility computations on triangulated terrains［J］.International Journal of Geographical Information Systems，1994，8（6）：515－531.

［6］ KIDNER D B，RALLINGS P J，WARE J A.Parallel processing for terrain analysis in GIS：Visibility as a case study［J］.GeoInformatica，1997，1（2）：183－207.

［7］ MILLS K，FOX1G，HEIMBACH R.Implementing an intervisibility analysis model on a parallel computing system［J］.Computers ＆ Geosciences，1992，18（8）：1047－1054.

［8］ 张林波，迟学斌，莫则尧，等.并行计算导论［M］.北京：清华大学出版社，2006.

［9］ 江岭，汤国安，刘凯，等.局部型地形因子并行计算方法研究［J］.地球信息科学学报，2012，6（14）：761－767.

［10］ CLARKE K C.Geocomputation′s future at the extremes：High performance computing and nanoclients［J］.Parallel Compu－ting，2003，29（10）：1281－1295.

［11］ FLORIANI L D，MARZANO P，PUPPO E.Line of sight communication on terrain models［J］.International Journal of Geographic Information Systems，1994，8（4）：329－342.

［12］ 应申，李霖，梅洋，等.增量法地形可视计算与分析［J］.测绘学报，2007，36（2）：192－197.

［13］ 王智杰，邱晓刚，李革.RSG地形通视性快速算法设计［J］.计算机仿真，2004，21（12）：92－95.

［14］ 鲁敏，张金芳，范植华，等.基于DEM的视域分析与计算［J］.计算机仿真，2006，23（5）：171－175.

［15］ 刘旭红，刘玉树，张国英.利用最大仰角插值技术的通视性分析算法研究［J］.计算机辅助设计与图形学学报，2005，17（5）：971－975.