EN中文

地球信息科学学报 >

2017 , Vol. 19 >Issue 4: 457 - 466

DOI: https://doi.org/10.3724/SP.J.1047.2017.00457

地球信息科学理论与方法

景观指数的并行计算方法

  • 刘洋 ,
  • 关庆锋 , *
展开
  • 1. 中国地质大学（武汉）国家地理信息系统工程技术研究中心,武汉 4300742. 中国地质大学（武汉）信息工程学院,武汉 430074
*通讯作者：关庆锋（1977-）,男,四川绵阳人,教授,主要从事高性能地理计算和地理计算智能研究。E-mail:

作者简介：刘 洋（1991-）,女,陕西西安人,硕士生,研究方向为高性能地理计算。E-mail:

收稿日期: 2016-07-19

  要求修回日期: 2016-11-01

  网络出版日期: 2017-04-20

基金资助

教育部高等学校博士学科点专项科研基金项目（20130145120013）

收起

A Parallel Algorithm for Landscape Metrics

  • LIU Yang ,
  • GUAN Qingfeng , *
Expand
  • 1. National Engineering Research Center of GIS, China University of Geosciences, Wuhan 430074, China2. School of Information Engineering, China University of Geosciences, Wuhan 430074, China
*Corresponding author: GUAN Qingfeng, E-mail:

Received date: 2016-07-19

  Request revised date: 2016-11-01

  Online published: 2017-04-20

Copyright

《地球信息科学学报》编辑部 所有
Fold

摘要

随着地理信息科学和系统的发展,GIS数据的时空分辨率和数据量呈现爆炸式的增长趋势。传统的基于个人计算机的景观指数计算软件难以有效快速地完成海量数据的空间分析。针对该问题,本文提出了一个高效的景观指数并行计算方法。首先对原有的并查集连通域标记算法进行了2点改进：① 在第2次遍历数据时,增加了计算斑块面积、周长等斑块基本信息的功能,为景观指数的计算提供必要参数;② 在第2次遍历过程中,增加了重新标记连续序号的功能,减少了原有算法在合并操作后造成的序号不连续,需要重新遍历数据的开销。在此基础上,本文利用MPI并行编程库,采用数据分割和主从进程协同的并行计算模式实现了景观指数的并行计算。实验表明,在保证计算正确性的基础上,本文的并行算法大幅度提高了景观指数的计算性能,为快速分析大规模数据的景观形态和格局提供了有效手段。

关键词： 景观指数; 连通区域算法; MPI; 并行计算

本文引用格式

刘洋 , 关庆锋 . 景观指数的并行计算方法[J]. 地球信息科学学报, 2017 , 19(4) : 457 -466 . DOI: 10.3724/SP.J.1047.2017.00457

Abstract

Landscape metrics have been widely used to quantitatively evaluate the spatial patterns of landscapes and to analyze the temporal dynamics of landscapes and their effects. However, when dealing with massive amounts of data, the calculation of landscape metrics requires large amount of computing time and extremely large memory size, which greatly decreases the feasibility in real-world applications. This paper presents a parallel algorithm to improve the performance of landscape metric calculation. First, the classical two-pass connected component labeling (CCL) algorithm was modified: (1) the calculations of some basic geometrical metrics of patches, such as areas and perimeters, are integrated into the second pass of the data for the calculation of landscape metrics; and (2) continuous patch IDs are generated along the second pass, to reduce the overheads for re-labeling. Then, a parallel algorithm consisting of a master process and a set of worker processes is designed and implemented using the C++ programming language and Message Passing Interface (MPI). In our parallel algorithm, the whole spatial domain is decomposed into multiple sub-domains and assigned to a set of concurrent processes. Each process uses the modified CCL algorithm to identify the patches within its assigned sub-domain and calculates the basic geometrical metrics of the patches. After gathering and merging the basic metrics from other processes, the master process calculates the final landscape metrics. The experiments using the land-use dataset of California showed that the computing time of landscape metrics was largely reduced using multiple processes. In conclusion, our parallel algorithm provides a high-performance solution for landscape metric calculation using massive and high-resolution datasets.

Key words： landscape metrics; connected component labeling algorithm; MPI; parallel computing

1 引言

景观指数是一种定量描述地理现象空间形态和格局的指标体系。通过计算和比较景观指数可以分析地理现象的时空特征差异性[1],探索城市扩展的动态变化规律[2-3],研究人类对河流、植被、气候地形等自然环境以及其它生物的影响[4-5]。常用的景观指数计算软件包括Fragstats、ArcGIS的Patch Analyst工具,以及IDRISI中的Pattern模块等。
随着地理信息科学和系统的发展,地理空间数据的时空分辨率和数据量呈现爆炸式的增长趋势。传统的基于个人计算机的景观指数计算软件存在着数据量和计算时间的限制,难以有效快速地完成海量数据的空间分析。例如,Fragstats是由美国俄勒冈州立大学森林科学系开发的一个景观指数计算软件[6],集成了60多种景观指数的计算功能,是科研工作中常用的指数计算软件,但即使在64位操作系统中进行计算,Fragstats的计算内存也仅限3GB（还有1GB用作运行程序）。因此,在使用传统的景观指数软件之前,常需要对图像数据做裁切处理以限制运算数据量。这不仅增加了计算景观指数的任务,还增加了计算误差,降低运行结果的精准度。
为了解决传统景观指数计算软件的局限性,本文针对海量数据的景观指数计算,提出了一种进程级别上的景观指数并行算法,并利用Message Passing Interface（MPI）并行编程库,在分布式内存环境下实现了景观指数的高性能计算。

2 景观指数的计算方法

2.1 连通区域标记算法

景观指数是以定量的方法来描述景观的空间格局。在景观生态学中,景观指数可以结合斑块、廊道、基质3种景观要素构成的模型来描述景观的空间结构[7]。斑块是指在景观中属性不同于周围环境的、相对均质的区域。廊道是指不同于两侧基质的狭长地带,可以看做是线状或带状的斑块。基质是指景观中面积最大、连通性最好的景观要素类型,如广阔的草原、沙漠等。在图像处理中,廊道和基质只是具有特定形状的斑块。也就是说,可以认为景观格局中所有的要素都是由斑块这一种景观要素组成的。在计算景观指数时,按照研究的空间尺度,可以划分为斑块（Patch）、类（Class）和景观（Landscape）3个级别,Fragstats软件在这3种景观级别上分别实现了计算多种景观指数的功能。参考其帮助文档可知[6],常用景观指数（如斑块面积、个数、密度、形状,分维数,景观多样性指数等）的计算参数主要包括斑块的面积、周长、个数等基本信息。由此,计算景观指数的基本步骤是先从分类图像中提取斑块的基本信息作为参数,再将其代入公式计算特定的景观指数。
为了从分类图像中提取斑块的基本信息,首先要在图像中识别出斑块。在分类栅格数据中,一个斑块实质上是由多个相互连通（相邻）并具有相同属性的栅格单元组成,即图像处理领域中的连通域（Connected Component）。因此,提取斑块的方法是在栅格数据中查找连通区域,即连通域标记过程（Connected Component Labeling, CCL）。连通区域标记算法的基本思想是将栅格图像中相互连通的区域标记为同一个标识值,以区别不同区域。如图1所示,连通域标记算法利用邻域模板对图像像元依次扫描,并计算相应的标识值。
Fig. 1 Scanning the pixels of image using a mask in the connected component labeling algorithm

图1 连通区域标记算法利用邻域模板扫描图像

最早的连通域标记算法是由Rosenfeld和Pfaltz于1966年提出[8]。然而,该算法需要遍历图像多次,计算复杂度较大。为了降低计算消耗的内存,减少遍历次数,Haralick改进了邻域模板的操作,使标记图像从上、下2个方向交替进行。其计算时间主要取决于图像的复杂度[9]。Dillencourt等将四叉树和并查集算法（Union-Find）应用于邻域模板操作中[10]。Frorio和Gustedt提出了一种结合并查集数据结构的连通域算法[11],只需遍历图像2次。然而,并查集算法采用的是链表存储结构,扫描需耗费大量时间。Suzuki等在总结前人的研究成果基础上对传统算法做了优化,增加标记局部合并最小值数组来提高计算效率[12]。但该算法的缺点是多次遍历后的标识序号在合并时会产生冲突,且图像信息复杂度会影响计算时间。同时,Wu等也优化了前人的算法,并提出将Suzuki等的2次遍历算法结合以数组为存储结构的并查集算法。这样既解决了合并冲突问题,也减小了原先链表结构的并查集连通域算法的内存开销[13-14]。此外,近年也有其他一些连通区域标记算法的研究进展,如Chang等提出的轮廓追踪算法,在处理二值图像时能达到较高的性能[15-16]。
Wu等提出的并查集连通区域算法的具体思路是：构建一个与分类图像（Image数组）一样大小的标识序号数据集（Label数组）,共遍历图像数据集2次[13-14]。第1次遍历利用邻域模板顺序扫描图像,并给标识序号数据集赋标识值。标记过程中,将当前像元P的值与邻域内1、2、3、4所在位置的像元值比较（图1）。若P与其它位置的像元值都不同,就在标识序号数据集上对应P的位置设置新标识值（按顺序的序号值）。在第1次遍历后,对于复杂的图像,同一个连通区域在标识序号数据集中会出现不同标识值。因此,扫描中需将同一连通域的等值标识记录在一个并查集Quick-Union树结构的合并数组中[17],其特点是可以快速找到等值标识中的最小标识值。第2次遍历是在标识序号数据集中对每个标识序号顺序扫描。同时,查找合并数组,将标识序号数据集中等值标识的位置替换成合并数组中的等值最小标识值（图2）。
这个算法的缺点是第2次遍历产生的最小等值标识值并不是一个连续标识序列值。在应用这样的标识序号值时存在许多不便,需要重新标识。本文改进了原有的算法,解决了不连续序号问题,并在原先的连通域算法中增加了计算斑块基本信息的过程。在此基础上,设计并实现了景观指数的并行算法以提高计算效率。
Fig. 2 CCL with two passes

图2 2次遍历的连通区域标记算法

2.2 改进的并查集连通区域标记算法

本文定义了一个FindPatch函数来封装改进的并查集连通区域标记算法。该算法的主要改进包括：① 增加了斑块（连通区域）的面积和周长等基本信息的计算过程,为后续的景观指数计算提供参数;② 在遍历过程中赋予相邻连通域以连续的标识值,以解决合并标识后产生的序号不连续现象,因此省去了为赋值新序号值而重新遍历并查集数组的时间开销。这一改进不仅应用于连通区域标识序号的合并过程中,还应用于并行算法的数据汇总过程。
2.2.1 斑块面积和周长的计算
为记录斑块的面积和周长等基本信息,根据已标记过斑块标识序号的图像,建立一个记录斑块面积、周长等基本信息的数组（Patch Array）。数组的序号与斑块标识序号相互对应,其中每个斑块的初始面积和周长均为0。在并查集连通区域算法的第2次遍历过程中,计算各斑块的面积和周长并更新斑块信息统计数组中相应的记录。以下是面积和周长的计算方法：
（1） 一特定斑块的面积实际是该斑块（连通区域）的栅格单元（即像元）个数乘以栅格单元的面积。计算方法是在逐行遍历图像像元时,根据每个像元对应的斑块标识序号,给信息统计数组中对应斑块的面积数值增加一个像元面积。在后面的并行算法中合并不同从进程中的同一斑块时,面积总值等于各个从进程中统计的斑块面积之和。
（2） 一特定斑块的周长是该斑块的所有边界像元的边界边长度之和。所谓边界边,是指斑块中某像元的一特定边为该斑块边界的一部分。计算方法是在遍历图像过程中,针对一特定像元,分别判断该像元的每条边邻接的像元属性值是否与该像元的属性值相同。若不相同,则该边被记为斑块的一条边界边,将其长度累加到斑块信息统计数组中相应斑块的周长值。如图3所示,该斑块的外边界为12个单位长度（假设像元为正方形）,内边界为4个单位长度,总周长为外边界与内边界的和,即16个单位长度。在后面的并行算法中合并不同从进程中的同一斑块时,周长总值的计算不是简单的相加,而是需要遍历子数据块的边界像元。如果分别属于2个子数据块的2个像元实际是同一斑块中的2个相邻像元,即这2个像元的共享边不是该斑块的边界边,应将这2部分数据的周长相加再减去共享边的2倍长度。
Fig. 3 Boundary edges of a patch and its perimeter calculation

图3 斑块的边界边和斑块周长计算

2.2.2 连续斑块标识序号的生成
Wu的并查集连通区域标记算法[13-14]中使用的并查集数组操作是Quick-Union算法[17],在2次遍历后合并斑块信息时会造成斑块标识序号不连续的现象,不利于后续斑块统计和计算工作。如图4所示,Quick-Union算法会将第1次扫描对应的第1、2、3号元素下存储的标识合并为一棵以1为根节点的树,即在数组中将原先的2号标识赋值为1,将原先的3号标识赋值为2。在第2次遍历时并查集依次查找每个标识对应的树根结点,即1、2、3号标识赋值为1,4号标识值不变。但这样会造成最终的斑块标识（1和4）并不连续。在计算和统计斑块基本信息（如面积和周长）时,为了处理不连续的标识序号,在斑块统计数组中除了存储斑块的面积和周长数值,还需要额外增加存储斑块标识序号的内存开销,以及查找斑块序号才可以获取对应基本信息的时间开销。为了降低查找时间,直接对不连续标识序号进行重标记生成连续序号,则需要重新遍历图像和合并数组,也会增加计算时间成本,且内存和时间开销与图像的破碎度呈正比。针对该问题,本文认为在并查集连通区域标记算法的第2次遍历中,可根据并查集合并数组按最小标识值生成的特性（最小堆的特性）,在第2次遍历图像的同时进行重标记,从而直接生成连续的斑块标识序号。
Fig. 4 Label merging by Quick-Union after the second pass

图4 图像2次遍历后Quick-Union并查集合并标记

考虑到在按行顺序遍历图像时,并查集合并数组中标识序号值是从小到大顺序生成的,且总是选择标识序号最小值作为根节点这一特质（最小堆的特性）,本文提出的改进算法将Quick-Union转换为Quick-Find并查集数组[17],在重标记的同时提高查找效率。其前提是在顺序遍历并查集数组的情况下重新标记序号,一定会先重标记每个连通域中最小标识序号值,即总是先重标记最小堆中树的根节点值,然后才会遇到已标记过根节点的树的子节点,局部看这个节点又是一个最小堆的根节点。
重标记的过程是按标识序号的顺序从小到大遍历并查集数组,设当前标识序号为i,新标识n初始为1,若该标识i未曾遇到过（第一次遇到）,则查找Quick-Union并查集数组第i个元素下取值j（j为i的父节点）,和第j个元素下取值k（k为j的父节点）。若i和k相等,说明遇到了该标识对应的树根节点。将第i个元素取值j修改为新标识值n,同时n=n+1;否则,k一定小于j,而j小于当前遍历标识i（k是j的父节点,j是i的父节点;因为算法总是按大小顺序合并取最小标识序号值作为存储的标识,所以k<j<i）。此时第j个元素下取值k已在之前遍历标识j时重标记过。为了将Quick-Union转化为Quick-Find,将第i个元素标识为k,即对应树根结点取值。若标识i已经遇到过了,则并查集中第i个元素的取值的就是新的重标记值,该值将会小于等 于i,不需更改;扫描下一个标识i=i+1,重复以上过程。完成扫描后,该并查集数组就从Quick-Union转化为Quick-Find思想,即从树型结构转换成为高效搜索的并查集模式（图5）。
Fig. 5 The relabeling process in the Union-Find Array

图5 重标记并查集数组的流程图

在连通区域标记算法中,2次遍历均是逐行从左至右的顺序遍历图像。被扫描到的标识序号的顺序也是从小到大,在并查集数组中查找标识序号时会按照从小到大的顺序进行。这符合以上并查集改进算法的前提假设特性。因此,可根据第1次遍历结束后得到的连通域（即斑块）总个数,定义斑块数组大小。在第2次图像扫描时,查找待合并标识序号时直接对斑块标识重新编号,并统计对应斑块的信息,就可实现从图像中计算各斑块面积和周长,以及统计斑块的数量。
在第2次逐行扫描图像（Image）,对应的标识序号数组（Label）,并查集数组（UF）,定义新标识序号初始值n=1,初始扫描位置i=0。重标记合并序号操作步骤如图6所示。
Fig. 6 The relabeling process in the second pass of CCL

图6 连通域标记算法在第2次遍历时重标记序号的合并流程图

3 并行算法的设计与实现

并行计算将计算任务分解为多个子任务,并利用多个处理器协同工作,是一种提高计算速度和处理能力的有效手段。近年来,并行计算已经被应用于连通域的计算中。马益杭等[18]研究了OpenMP和MPI并行技术下的连通域标记并行算法。其中,基于OpenMP的算法适用于单机多核的高性能计算机;基于MPI的算法更适合于分布式存储架构、集群以及超级计算机等硬件环境,它在大数据量和复杂图像上的并行效率具有显著优势。本文研究目的是为了突破景观指数计算在大数据量及复杂图像中计算耗时的瓶颈,因此采用MPI并行技术实现景观指数的并行算法。同时,本文在算法上对马益杭等的研究工作做了改进,将原算法中使用的链表并查集存储结构优化为本文改进的数组结构并查集算法（2.2.2节）,能够更高效地完成查询和并行数据汇总任务。Harrison等[19]将并查集连通域的MPI并行算法应用在三维格网图像分析中,并运用二叉空间分割（Binary Space Partitioning, BSP）树分配三维格网数据,实现并行算法的负载均衡。参考以上并行算法的研究成果,本文针对分类栅格图像,设计并实现了一个将景观指数计算方法并行化的策略,该算法可以提高在大数据量和复杂图像中的计算效率。
本文的景观指数并行算法采用了主从（Master-Worker）并行计算模式。如图7所示,主进程读入栅格数据,按行划分成多个子数据块,并分配给多个从进程。从进程在各自子数据块内利用并查集连通域标记算法,标识连通区域（斑块）并计算其基本信息（包括面积和周长等）,再汇总到主进程中合并信息,计算最终的景观指数。
Fig. 7 The basic idea of the parallel algorithm for landscape metrics

图7 景观指数并行算法的总体思路

图8展示了并行算法进程间的4次通讯过程。根据通讯的目的划分为3类,分别是分发数据（第1次通讯,橘色空心箭头）、组内消息传递（第2、3次 通讯,红色虚线箭头）和数据汇总（第4次通讯,紫色实心箭头）。
Fig. 8 Framework of parallel algorithm

图8 并行算法总体计算过程

并行算法中进程间的第1次通讯,是主进程将读取的子数据块,分发给从进程的过程。数据分割和分配如图7中（a）到（b）的过程,图像数据按行划分并分配给各个从进程。但在实际运算中,涉及边缘数据的邻域计算和最终数据的合并,因此会为第一个子数据块加一行下边缘行数据,为最后一个子数据块增加一行上边缘数据,为其他子数据块进行增加上、下各一行边缘数据。
从进程接收到分配的子数据块后,用改进的并查集连通域标记算法对数据进行斑块标识,并计算每个斑块的面积、周长等信息,记录到一个斑块信息统计数组（Patch Array）中。由于各个从进程的斑块ID都是从数字1开始标识,最后统计所有进程斑块信息时会产生标识冲突。为此,需要在进程间实现第2次的通信,步骤如下：
（1）各个从进程统计的斑块数为Ni,设第1个从进程计算的斑块总数为Np=N1;
（2）按进程编号顺序,第i个从进程将斑块总数Np传递给下一个从进程i+1,从进程i+1接收Np值,修改自己的斑块标识序号为原始斑块序号值加上Np,新的斑块总数值Np=Np+Ni+1,将Np传递给下一个从进程i+2;
（3）以此类推直到传递到最后一个从进程n,从进程n将斑块标识数值Nn加上Np,将新的Np值作为斑块总和传递给主进程,至此实现了所有斑块的全局标记。
第2次通讯如图7中（c）到（d）的过程。如果有3个从进程,各个从进程标识出的斑块个数分别为2、2、1,则经过以上操作后各个从进程的斑块标识变为全局唯一标识（Global Label）,主进程得到的斑块总和Np值为5（图7（d））。
以上操作得到的斑块总数Np大于等于整幅图像的所有连通区域的个数,图7的图像实际只有1个连通区域,而Np为5。这是因为原图像中的连通区域被分割成多个子数据块,分发到不同的从进程中单独计算斑块信息。为了合并原图像的同一连通域,需要根据子数据块的边缘数据信息,提取出属于同一个连通域的斑块标识序号,再将这些等值标识合并为一个斑块标识序号。从进程间的第3次通信过程如下：
（1）所有从进程将各自子数据块下边缘的行 标识以数组结构（LastRowLabel）传递给下一个从进程。
（2）第i个从进程接收到从进程i-1传递来的行标识作为从进程i数据块的上边缘标识行,再用图9中灰色模板扫描该上边缘行的图像数据。若存在需要合并的图像数据（边缘行中存在连通域）,则将待合并的标识序号值记录在一个等值合并数组中（Merge Array）,该数组每次只记录1对（2个,分别来自Last Row Label和从进程i中原本第一行标识）待合并标识序号值。
Fig. 9 Finding the patches to merge at the margins of sub-datasets

图9 记录子数据块边缘斑块信息

通过以上步骤,每个从进程完成了各自分配的数据处理任务,得到了2个数组,分别是斑块信息统计数组（Patch Array）和等值合并数组（Merge Array）。主进程获取了全局标识总个数Np,此时主进程需要汇总所有从进程的计算结果,获取斑块的基本信息。进程间的第4次通讯过程如下：
（1）所有从进程发送各自的2个数组：等值合并数组（Merge Array）和斑块信息统计数组（Patch Array）。
（2）主进程根据Np建立一个并查集数组（UF）,再接收所有进程发来的等值合并数组（Merge Array）,根据并查集算法将待合并的邻域标识序号全部合并为该邻域标识中的最小标识值。此时斑块标识符取值并不连续,使用我们改进的并查集算法可以直接给各个斑块重新赋以新的连续的标识值,这样节省了以往算法要重新扫描并查集数组再重标记的时间。
（3）主进程接收各个从进程发来的斑块信息统计数组（Patch Array）,再根据并查集数组UF替换掉原先的斑块标识,并将统一斑块的面积数值、周长等数值合并（见2.2.1节）。从而完成了并行计算斑块信息的过程。
在并行算法的第4次通讯中,主进程按照从进程标号顺序接收（Patch Array）数组,并用并查集数组UF合并斑块的统计信息（图8实心箭头）。在合并过程中,并查集数组UF是按照标识序号从小到大的顺序遍历的,因此使用图5的改进方法对斑块重新进行标识。
本文利用C++编程语言和多进程通信的并行计算模式（Message Passing Interface,MPI）实现了以上的斑块基本信息和景观指数的并行计算方法。数据I/O接口采用了开源的栅格数据I/O库Geospatial Data Abstraction Library（GDAL）。

4 景观指数并行算法实验

4.1 实验与实验环境

实验内容分为2个部分：① 通过比较本文的并行算法的计算结果与Fragstats软件的计算结果,验证了算法的正确性;② 以计算斑块面积指数为例,比较了串并算法的计算性能。实验数据采用美国加利福尼亚州1992年土地利用分类数据,图像大小3.61 G（共23 851像元×40 460像元）,分辨率为30 m,斑块总数13 705 118个。
实验采用的计算机集群包含452个计算节点,每个节点有2颗Intel Xeon E5-2670 CPU,共16核,内存64 G。计算节点间通过QDR Infiniband连接,每秒10千兆比特通讯速率。实验使用Scientific Linux 6.4操作系统下的G++4.7编译器,配合OpenMPI 1.6的MPI编程库,以及GDAL 1.11库。

4.2 景观指数计算正确性的验证

本文将并行算法计算出的景观指数值与Fragstats软件计算的景观指数做对比以验证计算结果的正确性。由于Fragstats软件限制了指数计算的图像大小[6]（内存最大限度3 GB,格网单元分辨率大于0.001 m）,因此裁取原图像中一部分数据作为实验数据,包括582像元×4173像元。
实验计算的景观指数按照斑块、类别和景观3个层次分类,计算了包括面积、统计指标、景观百分数、最大斑块指数、斑块密度、斑块周长、斑块形状指数、斑块分维数以及景观多样性指数等指标。在Fragstats的运算计算时间是37.75 s,而我们的并行程序利用12个进程在3.18 s内完成相同的计算。Fragstats计算结果与并行算法的计算结果完全一致,验证了并行算法的正确性。

4.3 并行算法性能评估

该组实验以加州1992年土地利用分类图像为输入数据（3.61 G）,分别测试了串行算法和并行算法的计算时间。如图10所示,串行计算的时间约为242.22 s。在并行计算实验中,当进程个数为2时,计算时间到达了最大值（247.44 s）,这是因为该算法是主从模式的并行计算。当只有2个进程时,主进程和从进程各有1个。从进程不仅要实现串行时的所有运算,还要和主进程通讯,这导致整体计算时间比串行还大。当使用64个进程时,计算时间缩短到10.18 s。
Fig. 10 Computing time of the parallel algorithm for landscape metrics

图10 并行景观指数算法的计算时间

为了进一步衡量算法的性能,本文使用并行计算常用的性能指标：加速比（Speed-up）和效率（Efficiency）。加速比的计算方法为式（1）,Tserial是串行计算时间,Tparallel是并行计算时间。理想状况下,其值等于并行计算的进程数p,此时认为该并行程序具有“线性加速比”。但实际情况下,并行算法往往会增加一些额外开销,因此效率（式（2））常被用于来反映随进程个数的变化并行算法的加速的程度。本文并行算法属于主从模式,当有p个进程时,实际参与计算的进程仅有p-1个（主进程并不参与从图像中计算指数的任务,p-1是从进程个数）,因此实验中的加速比和效率实际是用p-1值计算的。
Speedup = T serial T parallel (1)
Efficiency = Speedup p (2)
如图11、12所示,本文的并行算法在进程数p为4、8、16时（从进程数p-1分别为3、7、15）出现了加速比大于从进程个数p-1的情况。此时对应进程数下的效率值大于100%,称为“超线性加速比”状况,即p-1个从进程的并行运行速度是串行运行速度的p-1倍甚至p-1倍以上。并行程序在实际情况下往往由于通讯延迟、负载均衡等增加了运行开销,很难达到理想速度。但超线性加速比的情况也是存在的。当单机上串行运算需要大量内存空间时,会增加随机访问的寻址时间。而当物理内存空间不够时,还会使用硬盘作为虚拟内存,这都将降低串行运算的速率。而并行计算利用多个计算单元,其总体内存容量也大大增加,从而可出现超线性加速的情况。
Fig. 11 Speed-ups of the parallel algorithm for landscape metrics

图11 并行景观指数计算的加速比

Fig. 12 Efficiencies of the parallel algorithm for landscape metrics

图12 并行景观指数计算的效率

图13为并行计算中各个部分所花费时间的分步百分比堆积图,展示了并行算法在数据分发、数据运算（从进程提取斑块基本信息的计算和通信过程）和数据汇总上分别占用的时间百分比数。通过比较不同进程数下各步骤分别的耗费时长,可以看出算法的运行时间主要取决于分发时间和数据运算时间,汇总时间所占的比例一直保持一个比较低的水平,说明本文改进并查集算法具有较好的可扩展性。随着进程数的增加,数据分发时间的比例逐渐增加,并逐渐大于计算时间所占的比例。这是因为随着进程数的增加,由主进程向其它进程分发数据的通讯耗费越来越多,时间开销也随之增大。因此,将并行I/O方法结合到景观指数并行计算中将是我们下一步研究工作的重点之一,如采用GeoTIFF的并行I/O编程库pGTIOL[20]。
Fig. 13 Percentages of three parts in the parallel algorithm

图13 并行计算的分步百分比堆积时间图

5 结语

本文提出的景观指数的并行计算方法,首先对原有的并查集连通域标记算法进行了2点改进：① 在第2次遍历中,增加了计算斑块面积、周长等斑块基本信息的功能,为景观指数的计算提供必要参数;② 在第2次遍历过程中,增加了重新标记连续序号的功能,减少了原有算法在合并操作后造成的序号不连续,需要重新遍历数据的开销。在此基础上,本文利用MPI并行编程库,采用数据分割和主从进程协同的并行计算模式实现了景观指数的并行计算。实验表明,在保证计算正确性的基础上,本文的并行算法大幅度提高了景观指数的计算性能,为快速分析大规模数据的景观形态和格局提供了有效手段。

The authors have declared that no competing interests exist.

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
[1]
罗娅, 杨胜天,刘晓燕,等.黄河河口镇一潼关区间1998-2010年土地利用变化特征[J].地理学报,2014,69(1):42-53.运用土地利用变化重要性指数(Ci)、土地利用变化面积比重(D)和林草植被变化指数(I)等3个指标,基于干湿条件、地貌类型与坡度坡向,由宏观到微观,分析黄河河口镇—潼关区间1998-2010年土地利用变化的主要类型、变化的剧烈程度以及林草植被变化程度,以期为评价退耕还林(还草)政策在黄河河口镇—潼关区间的实施效果提供参考。结果表明:①1998-2010年,土地利用变化面积占总面积的19.19%,高覆盖度草地、有林地、其它林地面积明显增加,低覆盖度草地、农用地面积明显减少,土地利用变化十分显著;②空间上,土地利用变化主要发生在35°~38°N的区域,包括黄河以西的泾河支流马莲河流域、北洛河流域、延河流域、清涧河流域、无定河上中游区,以及黄河以东的沿黄区域;③土地利用变化过程中,低覆盖度草地转化为中覆盖度草地、中覆盖度草地转化为高覆盖度草地、农用地转化为其它林地、灌木林地转化为有林地是土地变化的主要类型;④土地利用变化呈现明显的干湿、地貌和坡度差异,但无明显的坡向差异。

DOI

[Luo Y, Yang S T, Liu X Y, et al.Land use change in the reach from Hekouzhen to Tongguan of the Yellow River during 1998-2010[J]. Journal of Geographical Sciences, 2014,69(1):42-53. ]

[2]
阿里木江·卡斯木,唐兵,古丽克孜·吐拉克.基于遥感和GIS的新疆绿洲城市扩展时空动态变化分析[J].冰川冻土,2013,35(4):1056-1064.以新疆16座城市为研究对象,运用多时相卫星遥感资料,利用遥感和GIS技术提取城市的空间特征信息,并通过城市扩展强度、扩展速率以及城市空间形态紧凑度和分形维数等指标定量分析了1990-2010年以来新疆城市扩展时空动态变化的过程和特征.结果表明:新疆城市的空间分布由其独特的自然环境和绿洲经济所决定,许多城市的发展依托于绿洲,城镇的形成具有"亲水性"的特点,由于新疆地域辽阔,交通干线附近往往形成城镇相对集中的城市带或城市群.南北疆城市之间在扩展速度及强度上存在较明显的差异,位于北疆天山北坡经济带上的城市大多属于扩展速度相对较快、扩展强度相对较大的扩张类型,而南疆城市的扩展速度则相对较慢、扩展强度也相对较小.通过对城市紧凑度和分维数变化的分析发现,新疆城市平均紧凑度呈下降趋势,而分形维数呈增长的趋势,说明整体上新疆城市空间结构有趋于松散化的趋势,城市的空间结构不够合理,有待改善.

DOI

[Kasimu A, Tang B, Tulake G.Analysis of the Spatial-Temporal Dynamic Changes of Urban Expansion in Oasis of Xinjiang Based on RS and GIS[J]. Journal of Glaciology and Geocryology, 2013,35(4):1056-1064. ]

[3]
吴琳娜,杨胜天,刘晓燕,等. 1976年以来北洛河流域土地利用变化对人类活动程度的响应[J].地理学报,2014,69(1):54-63.土地利用变化是反映人类活动程度的重要因子,分析土地利用时空变化规律,是揭示人类活动程度的有效方式。在遥感和地理信息系统技术的支持下,采用人机交互图像处理方法获取北洛河流域1976年、1998年和2010年土地利用数据;从土地利用变化速度、转移方向和土地利用程度方面,全面分析1976年以来北洛河流域人类活动作用下的土地利用时空变化规律。结果表明:①1976-1998年和1998-2010年间北洛河流域综合土地利用动态度从0.61增加到6.66;耕地和草地面积减少,减小速度分别从2.00%和2.69%增加到26.20%和23.33%,林地和城乡工矿居民用地面积增加,前者增加速度从5.93%增加到59.68%,后者增加速度从6.59%减小到3.52%。②两个时段土地利用类型转移方向都呈现出耕地和草地主要转化为林地,少部分耕地转化为城乡工矿居民用地的特点;③流域各县土地利用程度综合变化指数从-2~1之间,扩大为-27~4之间。证明流域内人类活动对自然环境的影响逐渐增大;主要表现为耕地、草地、林地和城乡工矿居民用地的相互转化;影响区域主要分布在流域上游吴旗县、富县、甘泉、黄陵县和洛川县;而人类活动作用导致的林地面积增大程度远远大于耕地和草地面积减小和城乡工矿居民用地面积增加程度。

DOI

[Wu L N, Yang S T, Liu X Y, et al.Response analysis of land use change to the degree of human activities in Beiluo River basin since 1976[J]. Journal of Geographical Sciences, 2014,69(1):54-63. ]

[4]
黄浦江,刘艳芳,刘畅,等. 基于RS与GIS的武汉城市湖泊演化研究[J].生态环境学报,2012(9):1588-1593.以武汉市主城区为例,利用1995和2005年Landsat 5两期TM影像与2000和2010年Landset 7两期ETM+影像的解译结果,基于景观分形理论与GIS的空间分析功能相结合的研究方法,分别构建湖泊变化强度指数和湖泊分形维数变化指数.从湖泊面积变化和湖泊形态变化,以及湖泊水域和其他土地利用类型的转移变化特点与影响进行综合分析.来丰富湖泊演化的分析方法,并总结高速城市化背景下湖泊的变化规律,从而更深入地了解和认识人类活动因素与湖泊水域动态变化之间的响应关系,同时提出若干城市湖泊治理与管理的方式.研究结果表明,(1)15年间的湖泊水域面积总量的变化呈现萎缩的趋势,年变化量在逐渐减小,湖泊萎缩的速度得到了一定的控制.(2)1995—2010年湖泊分形维数也呈逐期较小趋势,说明湖泊几何形状趋于简单化,人为活动对湖泊的影响加大.(3)15年来主城区28个主要湖泊的变化强度指数均为负值,呈萎缩趋势.但从2000年开始,少数湖泊变化强度指数为正值,萎缩趋势得到初步的控制.(4)1995—2010年间湖泊的水域面积主要转移成为建设用地和耕地.其转移面积占转移总面积的87.02%.围湖造田和城市化率的迅速提高对城市湖泊的演化影响深刻.(5)城市化背景下湖泊面积与形态变化的规律并非完全与城市化发展速率呈正比的关系,当城市化率上升到一定的阶段后,城市内部对生态环境的保护与要求也不断提高.为改善城市环境和城市内部的生态用地,湖泊应得到立法等强制性保护.

[Huang P J, Liu Y F, Liu C, et al. Study on evolution of urban lakes in Wuhan based on RS/GIS[J]. Ecology and Environment, 2012(9):1588-1593. ]

[5]
淡永利,王宏志,张欢,等. 2000-2010年武汉市中心城区湖泊景观变化[J].生态学报,2014,34(5):1311-1317.武汉雅称"百湖之市",湖泊是武汉市的重要名片。以2000、2005和2010年三期Landsat TM数据为数据源,采用归一化差异水体指数(NDWI)提取湖泊水体信息,建立了各个时期武汉市中心城区湖泊矢量图层,计算了湖泊的面积萎缩率、斑块分维数和破碎度等景观指数,对2000—2010年武汉市中心城区湖泊变化特征进行了分析。研究表明,武汉市各个湖泊均有不同程度的萎缩,湖泊面积萎缩率大小和其所隶属的环线及政策因素有较大的关系;各湖泊斑块分维数在1-1.3之间,并越来越接近于1,表明在人类活动持续影响下,其形状变得越来越规则;同时斑块数目增加,湖泊在面积萎缩的同时,也变得越来越破碎,研究显示主要是道路修建所致导致的湖泊分割,湖泊被分割后,自净能力下降,会导致水体污染而最终被填埋。总之,10年武汉市湖泊景观环境朝着不良方向发展,应该制定更严格的政策进行水域管理。

DOI

[Dan Y L, Wang H Z, Zhang H, et al.Lakes evolution of central Wuhan during 2000 to 2010[J]. Acta Ecologica Sinica, 2014,34(5):1311-1317. ]

[6]
Mcgarigal K S, Cushman S A, Neel M C, et al.FRAGSTATS: Spatial pattern analysis program for categorical maps[J]. 2002.Abstract pdf

[7]
Godron M, Baudry J, Forman R T.Thermodynamic foundation and information theory in understanding landscape heterogeneity[J]. Chinese Journal of Ecology, 1995.theoretical framework based on the dynamics and information theory is outlined. Thermodynamic principles explain the flows of energy that.in turn,produce heterogeneity in thedistribution of energy and material.We describe the( 1)trend toward a logarithmic vertical distribution of biomass. (2) thermedynamic linkage across scales (3 )importance of entropy in isolated,closed and open systems,as a landscape changes,and(4) three resistances opposing the law of zonality or horizontal homogenization. Together they provide insight into heterogeneity at the landscape scale. The broad theory of information is used to measure and understand structure,i.e. spatial configuration. Based on probability and conditional algebras,a simple method for measuring the amount of information in a landscape map is described.

[8]
Rosenfeld A, Pfaltz J L. Pfalz, J L.Sequential operations in digital picture processing[J]. Journal of the Acm, 1966,13(4):471-494.

[9]
Haralick R M.Some Neighborhood Operators[M]// Real-Time Parallel Computing. 1981:11-35.

[10]
Dillencourt M B.A general approach to connected-component labeling for arbitrary image representations[J]. Journal of the Acm, 1992,39(2):253-280.

[11]
Fiorio C, Gustedt J.Two linear time Union-Find strategies for image processing[J]. Theoretical Computer Science, 1996,154(2):165-181.ABSTRACT We consider Union-Find as an appropriate data structure to obtain two linear time algorithms for the segmentation of images. The linearity is obtained by restricting the order in which Union's are performed. For one algorithm the complexity bound is proven by amortizing the Find operations. For the other we use periodic updates to keep the relevant part of our Union-Find-tree of constant height. Both algorithms are generalized and lead to new linear strategies for Union-Find that are neither covered by the algorithm of Gabow and Tarjan (1984) nor by the one of Dillencourt et al. (1992).

DOI

[12]
Suzuki K, Horiba I, Sugie N.Linear-time connected-component labeling based on sequential local operations[J]. Computer Vision & Image Understanding, 2003,89(1):1-23.ABSTRACT This paper presents a fast algorithm for labeling connected components in binary images based on sequential local operations. A one-dimensional table, which memorizes label equivalences, is used for uniting equivalent labels successively during the operations in forward and backward raster directions. The proposed algorithm has a desirable characteristic: the execution time is directly proportional to the number of pixels in connected components in an image. By comparative evaluations, it has been shown that the efficiency of the proposed algorithm is superior to those of the conventional algorithms.

DOI

[13]
Wu K, Otoo E, Shoshani A.Optimizing connected component labeling algorithms[J]. Proceedings of SPIE - The International Society for Optical Engineering, 2005,5747(2):165-170.This paper presents two new strategies that can be used to greatly improve the speed of connected component labeling algorithms. To assign a label to a new object, most connected component labeling algorithms use a scanning step that examines some of its neighbors. The first strategy exploits the dependencies among them to reduce the number of neighbors examined. When considering 8-connected components in a 2D image, this can reduce the number of neighbors examined from four to one in many cases. The second strategy uses an array to store the equivalence information among the labels. This replaces the pointer based rooted trees used to store the same equivalence information. It reduces the memory required and also produces consecutive final labels. Using an array instead of the pointer based rooted trees speeds up the connected component labeling algorithms by a factor of 5 {approx} 100 in our tests on random binary images.

DOI

[14]
Wu K, Otoo E, Suzuki K.Optimizing two-pass connected-component labeling algorithms[J]. Formal Pattern Analysis & Applications, 2009,12(2):117-135.We present two optimization strategies to improve connected-component labeling algorithms. Taking together, they form an efficient two-pass labeling algorithm that is fast and theoretically optimal. The first optimization strategy reduces the number of neighboring pixels accessed through the use of a decision tree, and the second one streamlines the union-find algorithms used to track equivalent labels. We show that the first strategy reduces the average number of neighbors accessed by a factor of about 2. We prove our streamlined union-find algorithms have the same theoretical optimality as the more sophisticated ones in literature. This result generalizes an earlier one on using union-find in labeling algorithms by Fiorio and Gustedt (Theor Comput Sci 154(2):165 181, 1996). In tests, the new union-find algorithms improve a labeling algorithm by a factor of 4 or more. Through analyses and experiments, we demonstrate that our new two-pass labeling algorithm scales linearly with the number of pixels in the image, which is optimal in computational complexity theory. Furthermore, the new labeling algorithm outperforms the published labeling algorithms irrespective of test platforms. In comparing with the fastest known labeling algorithm for two-dimensional (2D) binary images called contour tracing algorithm, our new labeling algorithm is up to ten times faster than the contour tracing program distributed by the original authors.

DOI

[15]
Chang F, Chen C J, Lu C J.A linear-time component-labeling algorithm using contour tracing technique[J]. Computer Vision & Image Understanding, 2004,93(2):206-220.A new linear-time algorithm is presented in this paper that simultaneously labels connected components (to be referred to merely as components in this paper) and their contours in binary images. The main step of this algorithm is to use a contour tracing technique to detect the external contour and possible internal contours of each component, and also to identify and label the interior area of each component. Labeling is done in a single pass over the image, while contour points are revisited more than once, but no more than a constant number of times. Moreover, no re-labeling is required throughout the entire process, as it is required by other algorithms. Experimentation on various types of images (characters, half-tone pictures, photographs, newspaper, etc.) shows that our method outperforms methods that use the equivalence technique. Our algorithm not only labels components but also extracts component contours and sequential orders of contour points, which can be useful for many applications.

DOI

[16]
Chang F, Chen C J.A Component-labeling algorithm using contour tracing technique[C]// International Conference on Document Analysis and Recognition, 2003. Proceedings. IEEE, 2003:206-220.

[17]
Sedgewick R, Wayne K. Algorithms, Fourth Edition (Deluxe): Book and 24-Part Lecture Series[M]. Addison-Wesley Professional, 2015.

[18]
马益杭,占利军,谢传节,等.连通域标记算法的并行化研究[J].地理与地理信息科学,2013,29(4):67-71.连通域标记算法在地理栅格数据分析中有广泛应用，当面对大规模地理栅格数据时，连通域标记串行算法十分耗时，亟须算法并行化。但目前连通域标记算法还缺乏并行化，更缺乏对不同并行技术实现时的性能对比。该文对常用的连通域标记两遍扫描法进行了并行化设计，并分别利用OpenMP和MPI两种并行技术实现了不同版本的并行算法，以适用于单机多核、多机多处理器等不同的并行计算硬件环境。对所实现的并行算法在单节点、多节点的不同测试环境下，以不同数据规模和不同连通域复杂度情况的数据进行效率测试，结果表明：该算法均大幅缩短了运行时间；在数据量较小且连通域数目较少的情况下更适合使用OpenMP版本的并行算法；若图像数据规模较大时，MPI并行算法更快、更高效，但是在多节点的集群环境中，如果连通域情况复杂，进程数的增多并不能保证获得更好的加速效果。

DOI

[Ma Y H, Zhan L J, Xie C J, et al.Parallelization of connected component labeling algorithm[J]. Geography and Geo-Information Science, 2013,29(4):67-71. ]

[19]
Harrison C, Childs H, Gaither K P.Data-Parallel Mesh Connected Components Labeling and Analysis[C]// Eurographics Symposium on Parallel Graphics and Visualization, EGPGV 2011, Llandudno, Wales, UK, 2011. Proceedings. 2011:131-140.

[20]
胡树坚,关庆锋,龚君芳,等. pGTIOL:GeoTIFF数据并行I/O库[J].地球信息科学学报,2015,17(5):575-582.在地理栅格并行计算处理中,数据I/O已成为制约计算性能的主要瓶颈之一。本文针对该问题,首先分析广泛应用于GIS栅格数据存储的Geo TIFF格式,重点研究数据的2种存储模式(即条带存储与块状存储),并根据这2种存储方式,分别构建了栅格数据从逻辑结构向物理存储结构的映射模型。然后,针对地理空间并行计算的需要,提出了栅格数据的并行读写框架,并利用MPI并行I/O技术的文件视图方法,实现了Geo TIFF数据并行I/O库(p GTIOL)。结果表明,对比开源栅格空间数据转换库(GDAL)的主从I/O模式,本文提出的p GTIOL准确读写数据,具有更高的性能。该库隐藏了底层并行I/O的细节,提供简单易用的并行读写Geo TIFF栅格数据的接口,支持多数据类型和多种空间分割,实现了对条带存储与块状存储数据的异步并行读写,从而满足动态负载均衡的需求。

DOI

[Hu S J, Guan Q F, Gong J F, et al.pGTIOL: A Parallel GeoTIFF I/O Library[J]. Journal of Geo-Information Science, 2015,17(5):575-582. ]

Options
文章导航

/

〈 〉