陆冬华，赵英俊，张东辉，秦 凯

（1. 核工业北京地质研究院 遥感信息与图像分析技术国家级重点实验室，北京 100029）

遥感对地观测作为一项重要的空间探测技术，自20 世纪60 年代以来逐渐发展，已在气象、农业、林业、资源、环境、生态、水利、海洋、大气、水文、灾害、全球变化等民用领域，以及情报侦查、战场监视、导弹制导、隐身等军用领域发挥了积极作用[1]。遥感传感器技术的发展，使得获取的遥感数据空间分辨率和光谱分辨率不断提高，数据量不断增大，因此需要更强的计算资源对其进行处理[2]。美国最先进的WorldView3卫星遥感数据，其全色空间分辨率为0.3 m，多光谱16个波段的空间分辨率为1.24 m[3]，每100 km2融合后的数据约为17 GB；加拿大Itres公司的航空高光谱传感器CASI和SASI，分别拥有32和100个波段[4]，相对飞行高度为1 000 m时，空间分辨率分别为0.49 m和1.22 m，长度为10 km的单条测线获取的数据分别约为8 GB和15 GB。目前大多数的商业软件均未采用多线程数据处理技术，因此对数据量巨大的遥感图像进行分析时需要耗费大量的计算时间。

国内外学者对多线程和并行处理遥感图像方面开展了相关研究工作。QIN C Z[5]等利用多线程的方式对存储在文件系统中的遥感数据进行并行读写，降低了数据读写在遥感影像处理中占用的时间；刘世永[6]等利用MPI技术在超微高性能计算机上实现了多线程原始栅格影像的标准瓦片输出；吴炜[7]等针对集群计算机实现了高分辨率遥感影像的多尺度分割；周海芳[8]等实现了基于GPU的图像配准并行计算；LYU Z H[9]等针对高性能集群实现了遥感图像聚类分析。

遥感图像并行处理面临两个难题[5]：遥感数据量巨大，从计算机硬盘调入内存需耗费大量时间；存储遥感数据的文件格式种类繁多。目前国内外的研究主要侧重于多台计算机或集群的分布式计算以及GPU并行计算，该类方法虽可大幅提高数据的处理效率，但由于大多数科研人员仍以单台计算机运算为主要工作方式，且编程复杂度高，分布式需综合考虑不同计算机之间的数据传输与协同问题，GPU的程序开发需综合考虑GPU内不同类型显存的性能差异以及显存与内存之间传输数据的低效率问题[10]。随着计算机硬件的发展，Intel 和AMD等多核处理器的市场占有率逐步上升，而以往的分类算法一般为串行执行，无法发挥多核处理器的计算优势[11]。因此，本文重点针对单机多核环境下的遥感图像并行处理关键技术进行研究，一方面利用空间数据抽象库（GDAL）开源库解决不同格式数据的输入输出问题，另一方面通过OpenMP解决单台计算机CPU资源的充分利用问题，降低了开发难度，有效提高了图像处理速度。

1 不同类型遥感数据的输入和输出

由于多源遥感数据来自不同的卫星中心和传感器，遥感数据格式缺乏统一标准，不同的传感器数据具有不同的文件格式，甚至同一传感器，因其来源不同，文件格式也略有差异[12]。为了解决这一问题，1998年加拿大的Warmerdam F开始了GDAL项目的编写工作，并得到了许多个人和团体的支持[13]。截至目前GDAL已能支持超过150种空间数据格式，并被ESRI的ArcGIS10、Google Earth，跨平台的GRASS GIS 系统等多家遥感和GIS软件使用[14]。

GDAL 是一个独立的专业开源库，是在X/MIT 许可协议下的开源栅格空间数据转换库[13]。GDAL 的特点为[15]：①可扩展性、可移植性较好，若想支持某种新的数据格式，只需添加该格式驱动程序即可；②API 功能全面，可嵌入到其他程序中；③可跨平台运行；④支持常用的空间数据格式。

GDAL读写栅格数据最主要的核心类如图1所示。GDAL读写文件操作流程如图2所示。

GDALDataset类继承自GDALMajorObject，负责管理存储在硬盘上不同文件格式的遥感数据，利用GDALDataset中的函数可获取遥感数据的元数据信息、投影坐标信息、波段数量等；并可利用该类中的RasterIO函数完成多波段数据读写操作。

GDALDriver类是数据格式的抽象，每种格式均会生成GDALDriver的一个实例，所有属性和方法均来自于其对应的GDALDataset类[16]，通常用于创建不同格式的遥感数据文件。将遥感数据格式的名称传递给GDALDriver类，其将装载相应格式的信息，并创建相应的GDALDataset类，实现对文件的读写。

GDALRasterBand类负责管理GDALDataset类中的单一波段，获取该波段的信息和数据，实现数据的读写操作。

图2 GDAL读写文件操作流程图

2 OpenMP共享存储编程

OpenMP 是由世界上主要的计算机硬件和软件厂商提出的标准，跨平台、可伸缩的模型为并行程序的设计者在桌面电脑和巨型机上开发并行程序提供了简单灵活的接口，规范了一系列编译指导语句、子程序库和环境变量，在不降低性能的情况下克服了机器专用编译指导带来的移植性问题，得到了各界的认可。自发布以来，OpenMP已得到工业界和学术界的大力推动，成为共享主存多处理机上并行编程的事实标准。理论上，OpenMP更好地利用了共享内存体系结构，允许运行时调度，避免了消息传递的开销，并提供了细粒度和粗粒度运行机制[17-18]。

并行计算的主要实现过程为：主进程（处理器）将一个大任务（待处理数据）分派到若干个子进程（处理器）上进行处理，再由主进程负责收集不同子进程的数据处理结果并进行组合，以达到多处理器共同完成某一任务的目的[19]。首先由主线程执行程序的串行部分，再由派生出的其他线程执行程序的并行部分。相对于其他并行算法而言，OpenMP并行算法的开发十分简单，即便是已开发的单线程算法也十分容易改造，无需程序员进行复杂的线程创建、同步、负载平衡和销毁等工作。

OpenMP多线程的编写非常简单，包括两种方式：①在循环体for(int i=0;i＜N;i++)结构的前一行加入#pragma omp parallel for标识；②通过#pragma omp parallel{…} 标记并行块范围。系统将根据当前线程数自动对任务进行分配，利用omp_set_num_threads(int num)可设置并行区中活动的线程个数。使用OpenMP时，需在包含文件中加入函数库＜omp.h＞。常用的函数和编译指导语句如表1[20-21]所示。

表1 OpenMP常用函数和编译指导语句

3 遥感图像单机并行化处理

遥感图像并行化处理的前提是遥感图像各部分计算不存在依赖关系，即在串行计算时，先计算得到的结果不对后计算得到的结果产生影响。遥感所涉及的算法种类繁多，大多可实现并行化处理，如图像变换、图像配准、图像融合、目标提取等；但各方法处理过程完全不同，因此并行化算法的开发也不尽相同。为了降低算法的并行化难度，突出遥感并行化算法的一般过程，本文以Majority滤波和局部均方差计算为例进行阐述。

3.1 遥感图像的分块

由于遥感图像输入输出数据占用了大量的磁盘空间，无法一次性将所有数据读入内存，因此必须将其进行分块处理，每次仅读取数据的一个分块到内存中，处理后再输出到磁盘空间。分块方法通常包括行优先、列优先和平均分块[5]3种。行优先处理方法，根据内存大小，将遥感图像分成行数相近的若干个数据块；列优先处理方法，根据内存大小，将遥感图像分成列数相近的若干个数据块；平均分块处理方法，根据内存大小，将遥感图像分成大小相近的数据块，每个数据块中数据的行数和列数相当。

处理过程中，由于Majority滤波和局部均方差计算均需对M×M窗口内的像素值进行统计和计算，为了保证分块数据边界计算的正确，输入的数据应在分块的基础上，在上下左右各增加(int)(M/2)，即若第i个分块的左上角坐标为(xoff, yoff)，大小为(xBlockSize,yBlockSize)，则读取的数据范围应为(xoff-(int)(M/2),yoff-(int)(M/2))，大小应为（xBlockSize+M－1,yBlockSize+M－1），其中M为基数。

3.2 遥感并行算法的实现

首先利用GDAL打开待处理的图像文件，读取图像的行列数、坐标信息、投影信息等基本信息；再计算分块大小，确定每个分块的左上角坐标以及行数和列数，利用OpenMP开始多线程计算，利用GDAL读取分块数据，并进行Majority计算；然后利用GDAL写入处理后的数据；当所有线程都结束后，关闭遥感数据文件。

以数据分块结果构建for循环体，并在for循环体前加入#pragma omp parallel for编译指导语句，由于在每个线程均存在IO操作，而单机多线程的读写将显著降低读写效率，因此在IO语句前加入#pragma omp critical原子处理编译指导语句，避免多线程同时进行IO操作。for循环体前的编译指导语句为：

//#pragma omp parallel for

for (int iBlock= 0; iBlock ＜ 分块总数 ;++ iBlock)

IO操作前的编译指导语句为：

#pragma omp critical //原子处理

iBand-＞RasterIO(GF_Write, xOff, yOff, xBlockSize,yBlockSize,……);

4 数据实验与分析

4.1 实验算法

1）Majority算法，主要应用于遥感图像分类的后处理。利用M×N的移动窗口，对窗口内的图像像素值进行统计，将出现频率最多的像素值设置为M×N窗口的中心像素值，M和N通常为基数。

2）局部均方差算法，主要应用于图像增强和边缘检测。利用M×N的移动窗口，对窗口内的图像像素值进行标准方差计算。

4.2 数据实验

实验平台为惠普计算机图形工作站Z800，内存为64 GB，2颗Intel Xeon x5690 CPU，主频为3.46 GHz，每颗CPU为6核，开发平台为Viusal Studio，将项目→属性→VC++目录→包含目录和库目录中分别设置为GDAL的头文件目录和LIB文件目录。将GDAL的Bin文件夹设置于环境变量中，或直接拷贝Bin文件夹中的文件到工程数据文件目录，完成对GDAL库编译环境的设置；然后在C/C++的语言选项中设置OpenMP支持为“是(/openmp)”，完成对OpenMP编译环境的设置。

以Majority算法为实验对象，输入图像大小为3 000×2 500的单波段分类数据，输出为Majority滤波结果，分别测试在窗口尺寸为3×3和5×5，并发线程数为1～12情况下的计算性能，如表2所示，可以看出，线程并发数由1增加到4的过程中，计算时间逐渐缩短，当并发线程增加到5之后，计算时间开始变长，达到12时计算时间最长。图3更加清晰地反映了这一过程。

表2 Majority算法不同线程耗时表

图3 Majority算法不同线程耗时曲线

Majority算法实验说明尽管计算机CPU共有12个核心，但计算时间并未随并发数量的增加而逐渐缩短，而是在4个线程并发时达到最高效率。为了进一步研究并行算法的效率，本文继续对局部均方差算法进行实验。

以局部均方差算法为实验对象，输入图像大小为3 000×2 500的单波段16 bit整型灰度数据，输出为方差计算结果，分别测试在窗口尺寸为3×3和11×11，并发线程数为1～12情况下的计算性能，具体如表3所示，可以看出，从线程并发数由1增加到11的过程中，计算时间基本上逐渐缩短，达到12后略有延长。图4更加清晰地反映了这一过程。

表3 局部均方差算法不同线程耗时表

图4 局部均方差算法不同线程耗时曲线

4.3 实验结果分析

实验结果表明，局部均方差算法通过OpenMP的加速效果较好；而当并发线程数为4时，Majority算法的计算效率最高。为了进一步分析导致实验结果的原因，对参与计算的各函数进行逐一计算，分析各函数对并行效率的影响。其代码主要结构为：

omp_set_num_threads(nThread);

#pragma omp parallel for

for(int i= 0;i＜M;i++)

for(int j=0;j＜N;j++)

{

//测试代码

…

}

当测试代码为空或仅包含基本的数学运算（加、减、乘、除、开方运算）时，并发线程由1增加至12时，程序运行时间逐渐缩短。

当测试代码为int a[5000]或int *a=new int[5000/nThread]时，并发线程由1增加至6时，程序运行时间逐渐缩短，增加至7时，运行时间有所延长，但随后又逐渐缩短。

当测试代码为int *a= new int[10]时，程序运行时间与并发线程数的变化与Majority算法的规律相同，即1～4个线程，程序运行时间逐渐缩短，之后随并发数增加程序运行时间逐渐延长。

Majority算法中使用了C++STL map类，该类对并发线程的加速也有较大影响，当测试代码部分为map＜int,int＞ statics变量声明时，1～5个并发线程数时，程序运行时间缩短，随后延长；将map＜int,int＞ statics变量声明放在for循环体之外，测试代码仅为statics.find(2)时，1～4个线程程序运行时间逐渐缩短，之后随并发数增加程序运行时间逐渐延长。

由此可知，影响Majority算法加速效果的主要原因为C++STL map类的变量声明和查找函数；由map类中的哪个函数产生的影响还需进一步研究，但比较明确的是，程序在运行过程中申请内存或新建数组，会在一定程度上影响OpenMP的加速效果。

5 结 语

本文将GDAL开源库与OpenMP相结合，开展了遥感数据在单机多核计算机环境下的处理实验，利用GDAL实现了多种遥感数据的读写操作，利用OpenMP将串行程序改写为并行程序。GDAL读取遥感数据后，应根据内存大小对遥感数据进行分块处理，利用OpenMP将for循环体并行化，可有效提高遥感图像的处理效率。

实验结果表明，并行程序对遥感图像处理的加速效果与处理算法紧密相关，不同的算法需要调整并发线程数才能达到最佳的处理效率，单纯的增加线程有时不但不能提高处理效率，反而会延长处理时间。