随着卫星遥感技术的快速发展, 遥感图像的数据量越来越大。比如, SPOT-5卫星图像存储了5.76×10⁸个像素, WorldView-2卫星图像存储了2.62× 10⁸个像素。如此大的数据量给遥感数据处理带来极大的压力和挑战, 主要表现在程序运行时间过长、计算机内存不足。当前计算机CPU主频和内存容量增长比较缓慢, 因此依靠硬件性能的提高解决遥感图像处理时间过长和内存不足的问题比较困难。如何提高遥感图像处理的时间效率和空间效率, 已成为遥感图像处理领域面临的最紧迫的问题之一。

为了提高遥感图像的处理速度, 各种并行计算方法被引入遥感图像处理过程中, 包括基于MPI的遥感图像快速处理方法^[1], 它采用MPI标准并行编程环境和基于集群的多机处理系统对遥感图像进行信息提取; 基于多线程的遥感图像快速处理方法^[2], 其中主线程负责对图像进行分块, 创建子线程, 具体由子线程完成各个分块的处理。在单核处理器中, 这种方法一般能够节省30%的处理时间; 基于GPU与CPU异构模式的遥感图像快速处理方法^[3-5], 在传统计算机系统中加入GPU作为加速部件并配合CPU共同承担计算任务的新型系统, 相比于传统的单纯以CPU作为计算部件的同构计算系统, 异构系统优势明显; 基于GPU的遥感图像快速处理方法^[6-9], 与CPU只有少数几个核(多核)相比, GPU几百个核(众核)具有更大的并行计算潜力, 可用于处理计算密度高、逻辑分支简单的大规模数据并行任务。GPU技术是一种显卡的并行计算技术, 它以大量线程实现面向吞吐量的数据并行计算, 相对于其他技术而言, GPU技术更适合像素数目庞大、逻辑分支简单、迭代次数多的遥感图像处理算法。

为了解决大数据量遥感图像处理内存不足的问题, 最基本的想法就是将图像划分为若干块, 每次只处理一块的图像数据, 这种方法可以显著降低内存的使用量。然而, 由于需要频繁进行I/O操作, 对于需要若干次迭代收敛的算法, 图像处理速度将会下降几倍甚至数十倍。解决这个问题的一种方法是采用分布式计算技术^[10-12], 每个计算节点处理的分块大小应当小于其内存容量, 避免了迭代过程中多次读写磁盘; 另一种方法是采用内存映射文件技术^[13-14]。内存映射文件建立了从磁盘中文件到特定内存区域的映射关系, 通过这种映射关系可以显著提升文件访问速度。

虽然集群、网格和分布式系统都能够提高遥感图像处理算法的时空效率, 但是它们需要特定的计算基础设施, 不仅代价昂贵而且环境配置和使用过程复杂。本文研究了普通计算机环境下的遥感图像处理算法优化方法, 其基本思想是利用GPU并行计算技术提高遥感图像的处理速度, 采用图像分块处理方法降低内存需求, 利用内存映射文件技术提高图像文件访问速度, 最终实现在不增加任何计算资源的情况下提高遥感图像处理算法的时空效率。

1 K-Means算法、CUDA和内存映射文件

在介绍本文的方法之前, 首先回顾一下K-Means聚类算法、CUDA平台和内存映射文件的相关知识。

1.1 K-Means聚类算法

K-Means聚类算法是应用较广的一种非监督分类算法, 其基本思想是将所有样本聚类为$k$个簇, 每个簇内具有较高的相似度。给定训练样本$X = $ $\{ {x^{(1)}}, {x^{(2)}}, \cdots, {x^{(n)}}\} $, 每个${x^{\left( i \right)}} \in {R^p}$, K-Means聚类算法可以描述如下:

(1) 随机选取k个聚类中心: ${\mu _1}, {\mu _2}, \cdots, {\mu _k} \in {R^p}$

(2) 重复下面过程直至收敛: {

1) 对于每一个样本i, 计算其应该属于的类:

$ {c^{\left( i \right)}}: = \arg \;\mathop {\min }\limits_j {\left\| {{x^{\left( i \right)}}-{\mu _j}} \right\|^2}\;\;i = 1, 2, \cdots, n\;\;j = 1, 2, \cdots, k $

(1)

2) 对于每一个类j, 重新计算该类的聚类中心

$ \begin{array}{*{20}{c}} {{\mu ^j}: = \frac{{\sum\limits_{i = 1}^n {1 \cdot \{ {c^{(i)}} = j\} \cdot {x^{(i)}}} }}{{\sum\limits_{i = 1}^n {1 \cdot \{ {c^{(i)}} = j\} } }}}&{j = 1, 2, \cdots, k} \end{array} $

(2)

}

上式中, ${c^{(i)}}$代表样本$i$与$k$个类中距离最近的那个类, ${c^{(i)}}$的值是1到$k$中的一个。如果${c^{(i)}} = j$, 则$\{ {c^{(i)}} = j\} $的值为1;如果${c^{(i)}} \ne j$, 则$\{ {c^{(i)}} = j\} $的值为0。公式(2)中的分母表示第$j$类中样本的数量, 分子表示第$j$类中所有样本向量的总和, 因此除法计算的结果是第$j$类的中心向量。K-Means聚类算法的优点是简单、快速, 对于处理大数据集是相对可伸缩和高效率; 缺点是聚类结果对初始值敏感, 易受噪声和孤立点的影响。K-Means聚类算法的时间复杂度是$O(nkt)$, 其中$n$为样本数目, $k$为聚类中心数目, $t$为迭代次数, 通常$k \ll n$且$t \ll n$; 空间复杂度是$O((n + k) \cdot p)$, 其中$n$为样本数目, $k$为聚类中心数目, $k \ll n$, $p$为样本维数。

1.2 CUDA平台

CUDA(Compute Unified Device Architecture, 统一计算架构)是NVIDIA(英伟达)公司推出的一种并行计算架构^[15-17]。该架构通过利用GPU(Graphic Processing Unit, 图形处理器)的处理能力, 可大幅提升计算性能。用户可使用NVIDIA的GeForce 8及以后的GPU进行计算, 亦可以利用GPU作为C-编译器的开发环境。以GeForce 8800 GTX为例, 其核心拥有128个处理器。利用CUDA技术可以将那些内处理器串通起来, 成为线程处理器解决数据密集的计算。各个内处理器能够交换、同步和共享数据。利用NVIDIA的C-编译器, 就能利用这些功能。

从硬件角度讲, GPU是由多个SM(Streaming Multiprocessor, 流多处理器)组成, 每个SM包含有多个SP(Streaming Processor, 流处理器), SP是最基本的处理单元。从软件角度讲, CUDA是SIMT(单指令多线程)架构, grid(线程网格)和block(线程块)是thread(线程)的组织形式。最小的逻辑单位是一个thread, 若干个thread(128~512)组成一个block, 若干个block组成一个grid对应一个kernel(函数)。一个SM负责多个block的计算任务(一般为8个), 每个SP在任一时刻只负责一个thread。SM中有共享内存、寄存器和L1缓存, 因此block内可以共享公共内存, 每个thread拥有自己的局部内存。warp(线程束)是SM调度和执行的基本单位, SIMT架构使得同一个warp里的线程根据不同的数据执行相同的指令。SM一次只能运算一个block里的一个warp, 如果warp中有线程的数据没有取到, 那么SM调度下一个warp进行运算。

1.3 内存映射文件

内存映射文件(Memory Mapping File, MMF)是由一个文件到一块内存的映射。Win32提供了允许应用程序把文件映射到一个进程中的函数(CreateFileMapping)。内存映射文件技术提供了一种独特的内存管理特征, 它允许应用程序与通过指针访问动态内存相同的方式访问磁盘上的文件。因此, 可以在进程地址空间中将磁盘上的文件部分或全部映射到特定的地址范围, 然后通过指针就可以访问内存映射文件中的内容。一旦某个文件被映射, 就可以访问它, 就像整个文件已经被加载到内存中一样, 从而可以不必对文件执行I/O操作, 这对于大数据文件来说存取的效率极高。

内存映射文件技术比磁盘I/O操作快得多。应用程序可以像读写内存数据那样直接访问或修改映射后的内存文件, 而不需要像普通I/O操作那样从文件开头进行寻址。内存映射文件的另外一个好处是“延迟加载”, 当文件大小远大于内存容量时, 将整个文件加载到内存需要操作系统反复将文件读入内存和写磁盘, 这会引起严重的“抖动”。内存映射文件每次只需要加载较小的页面大小区域, 因此基于内存映射文件的程序能够处理任意大的数据文件。

2 本文方法

为了解决遥感图像处理内存不足的问题, 本文采用了图像分块的思想, 即将遥感图像划分为若干块, 每次只处理一个分块。然而, 分块处理通常需要图像数据在内存和硬盘之间频繁交换。相对于内存访问来说, 磁盘读写速度非常慢, 从而导致遥感图像处理时间比较长。本文将基于CUDA并行计算技术、图像分块处理方法和内存映射文件技术3种手段相结合, 克服内存空间不足和磁盘访问速度慢的瓶颈, 提高遥感图像处理算法的时空效率。

2.1 基于CUDA的遥感图像K-Means聚类算法的并行化

GPU拥有众多的核(SP), 非常适合遥感图像的并行处理。由于SP在任一时刻只能执行一个线程, 处理一个像素, 因此需要对遥感图像进行分块, 每次GPU处理一块图像数据。分块大小通常设置为256×256或512×512, 即一个grid包含256/512个block, 一个block包含256/512个thread, 每个thread处理一个像素。太大或太小的分块都不利于充分发挥GPU的计算性能。CUDA能够将计算任务映射为大量的可以并行执行的线程, 并由硬件动态调度和执行这些线程。因此, 算法只需要将图像划分为若干块, 然后将每个块发送给GPU处理, 不需要关注其内部的线程调度和执行问题。

对于遥感图像的K-Means聚类算法, 假设相对于图像大小内存空间足够大, 基于CUDA平台的并行化处理算法可以描述如下:

(1) 使用RasterIO函数将遥感图像读入到内存中;

(2) 随机选取$k$个聚类中心为${\mu _{^1}}, {\mu _{^2}}, \cdots, {\mu _{^k}} \in {R^p}$;

(3) 将遥感图像划分为若干块, 每一块的大小等于(最后一行或一列可能小于)256×256/512×512;

(4) 重复下面过程直至收敛: {

1) 对于每一个分块, 将其数据从内存拷贝到设备端;

2) 将每一个类的聚类中心从内存拷贝到设备端;

3) 利用公式(1)在设备端对分块图像做并行计算;

4) 将分块计算结果(即样本的类别标记)从设备端拷贝到内存;

5) 处理下一个分块, 直到所有分块处理完成为止;

6) 对于每一个类, 利用公式(2)重新计算该类的聚类中心;

}

基于CUDA的K-Means聚类算法的时间复杂度是$O((2np/C{\rm{ + }}\alpha nk/m) \cdot t)$, 其中, n为样本数目, p为样本维数, C为从内存到设备端的数据拷贝速度, 因为需要双向拷贝, 所以乘以2; $\alpha $为CPU与GPU单核速度比值, k为聚类中心数目, t为迭代次数, m为流处理器数目。基于CUDA的K-Means聚类算法的空间复杂度是$O((n + k) \cdot p)$。

2.2 基于CUDA和MMF遥感图像K- Means聚类算法优化

由于操作系统和其他程序运行占用了较多的内存空间, 遥感图像处理程序可以利用的内存空间通常只有几百兆字节, 处理大小不足百兆字节的遥感图像。为了能够处理大数据量的遥感图像, 必须对图像进行分块处理。分块处理方法需要将各个分块依次读入到内存中进行计算, 然后再将中间计算结果写到磁盘中。对于需要若干次迭代的遥感图像处理算法(如K-Means聚类算法), 频繁地读写磁盘文件将会大大增加程序的执行时间。

为了解决这个问题, 本文首先对遥感图像数据进行重新组织和存储, 然后利用内存映射文件技术提高图像分块的读写速度。传统的遥感图像数据组织格式有三种类型: BSQ、BIP和BIL。BSQ(Band SeQuential Format, 波段顺序格式):每行数据后面紧接着同一波段的下一行数据; BIP(Band Interleaved by Pixel format, 波段按像元交叉格式):按顺序存储第一个像元的所有波段, 接着是第二个像元的所有波段, 然后是第三个像元的所有波段, 交叉存储直到像元总数为止; BIL(Band Interleaved by Line Format, 波段按行交叉格式):先存储第一行的第一个波段, 接着是第一行的第二个波段, 然后是第一行的第三个波段, 交叉存储直到最后一行最后一个波段为止。

BSQ格式适合于对单个波段某个空间范围内的数据存取, BIP格式适合于像元光谱数据的存取, BIL格式适合于存取某些行的像元光谱数据。如果需要读取遥感图像某个分块的所有波段数据, BSQ、BIP和BIL格式都不是最优方案, 因为分块各个波段的数据被分散存储在磁盘的不同区域, 无法利用MMF技术进行存取。为此, 本文提出第四种存储格式BIT(Band Interleaved by Tile Format, 波段按瓦片交叉格式):先存储第一个瓦片的第一个波段, 接着是第一个瓦片的第二个波段, 然后是第一个瓦片的第三个波段, 交叉存储直到最后一个瓦片的最后一个波段为止。每个瓦片对应于遥感图像的一个分块。这种存储格式使同一个瓦片的所有波段连续存储, 因而可以利用MMF技术通过内存指针连续读写瓦片的各个波段数据。

基于CUDA和MMF的K-Means聚类算法具体实现描述如下:

(1) 将遥感图像按行列切分成M个瓦片(Tile), 其中每个瓦片的大小为B(512×512或其他尺寸);

(2) 建立图像文件tempImage保存原始图像中的每个瓦片, 建立类别文件tempClass保存图像中每个像元的类别; 此操作目的是建立内存映射文件;

(3) For i = 1 to M Do Begin

1) 使用RasterIO函数将第i个瓦片读入到内存中;

2) 将第i个瓦片写入到图像文件tempImage中;

End

(4) 随机选取$k$个聚类中心为: ${\mu _{^1}}, {\mu _{^2}} \cdots, {\mu _{^k}} \in {R^p}$;

(5) 打开图像文件tempImage并映射到内存指针ptrImage, 打开类别文件tempClass并映射到内存指针ptrClass;

(6) 重复下面过程直至收敛: {

For i = 1 to M Do Begin

1) 使用内存指针ptrImage读取第i个瓦片的图像数据, 并拷贝到设备端;

2) 将每一个类的聚类中心从内存拷贝到设备端;

3) 利用公式(1)在设备端对瓦片数据做并行计算;

4) 将分块计算结果(即样本的类别标记)从设备端拷贝到内存;

5) 使用内存指针ptrClass将类别信息写到类别文件tempClass中;

End

对于每一个类j, 重新计算该类的聚类中心

$ {\mu ^j}: = \frac{{\sum\limits_{r = 1}^{\text{M}} {\sum\limits_{s = 1}^B {1 \cdot \left\{ {{c^{\left( {r, s} \right)}} = j} \right\} \cdot {x^{\left( {r, s} \right)}}} } }}{{\sum\limits_{r = 1}^{\text{M}} {\sum\limits_{s = 1}^B {1 \cdot \left\{ {{c^{\left( {r, s} \right)}} = j} \right\}} } }}\;\;\;j = 1, 2, \cdots, k $

(3)

其中, ${x^{(r, s)}}$表示通过内存指针ptrImage读取的第r个瓦片第s个像元的光谱向量, ${c^{(r, s)}}$表示通过内存指针ptrClass读取的第r个瓦片第s个像元的类别, 而${\mu ^j}$表示第j类的聚类中心(光谱向量均值)。

}

基于CUDA和MMF的K-Means聚类算法时间复杂度是$\text{O}(2np/D + (2np/C{\rm{ + }}\alpha nk/m) \cdot t)$, 其中: n为样本数目, p为样本维数, C为从内存到设备端的数据拷贝速度, D为磁盘读写速度, $D \ll C$; $\alpha $为CPU与GPU单核速度比值, k为聚类中心数目, t为迭代次数, m为流处理器数目。算法空间复杂度是$O((B + k) \cdot p)$, 其中B为瓦片大小, 通常$B \ll n$, $k$为聚类中心数目, p为样本维数。

在处理大数据量遥感图像时, 如果仅仅采用分块方法而没有使用MMF技术, 那么每一轮迭代中所有遥感图像分块都要读写磁盘各一次, 这种情况下算法的时间复杂度是$\text{O}((2np/D + 2np/C{\rm{ + }}\alpha nk/m) \cdot t)$, $D \ll C$; 算法的空间复杂度是$\text{O}((B + k) \cdot p)$, $B \ll n$。

2.3 复杂度对比与分析

本文已讨论了4种实现策略下的K-Means聚类算法。第1种策略基于CPU和内存实现遥感图像的聚类, 简称“CPU+内存”; 第2种策略基于CUDA和内存实现, 简称“CUDA+内存”; 第3种策略基于CUDA和内存映射文件实现, 简称“CUDA+ MMF”, 第4种策略基于CUDA和图像分块方法实现, 简称“CUDA+分块”。表 1中包含了各种实现策略的算法复杂度。

关于算法的时间复杂度, 由于GPU通常拥有数百个核, 即m值很大, 理论上时间复杂度可以达到原有算法的1/m, 但在实现时由于CUDA复杂的执行模型、CPU与GPU时钟频率的差异以及数据在内存和设备端双向拷贝花费的时间等因素影响, 其实际时间复杂度一般高于理论时间复杂度。“CUDA+ MMF”策略只比“CUDA+内存”策略增加了一次图像读取和写入磁盘的时间, 而“CUDA+分块”策略在每一轮迭代都需要增加一次图像读取和写入磁盘的时间, 因此其时间复杂度大大提高。

关于算法的空间复杂度, “CUDA+MMF”和“CUDA+分块”策略的空间复杂度是相同的。由于B远小于n, 而k可以忽略不计, 因此后两种实现策略的空间复杂度远远低于前两种实现策略。

3 实验与分析

本文中采用的实验数据是广州省雷州湾的SPOT-5多光谱遥感图像, 图像尺寸6 000×6 000像素, 大小137MB, 包含绿、红、近红外、短波红外四个波段, 其假彩色图像如图 1所示。

计算机的硬件环境配置如下: CPU(Intel Core 2 Duo Processor E6550, 4M Cache, 2.33 GHz); 显卡(NVIDIA GeForce GTX 650 Ti, 768个流处理器, GPU核心时钟324MHz, 内存时钟162MHz, 显存1.0GB); 内存2.0GB。软件环境配置如下:操作系统Windows XP; 运行平台CUDA 5.0;开发环境Visual Studio 2010;编程语言C++。

K-Means聚类算法设定的收敛条件是连续两次迭代得到的聚类中心完全相同或很相似, 可以采用绝对距离、欧氏距离、马氏距离、曼哈顿距离、夹角余弦^[18]距离度量公式进行相似性评价。本文采用了夹角余弦法, 因为它的计算结果不依赖于图像的数据类型和取值范围, 鲁棒性较好。当对应聚类中心的夹角余弦值都大于某个给定的阈值(如0.999 9)时, 算法终止。

为了检验K-Means聚类算法在不同实现策略下的性能, 本文使用了一幅裁切后的小图像(600×600像素)和整景图像(6 000×6 000像素)进行了实验。由于K-Means聚类算法的初始聚类中心是随机选择的, 所以每次算法执行需要的迭代次数都太不相同。为此本文做了12组对比实验, 每一组实验随机选择15个初始聚类中心, 然后执行各种K-Means聚类算法的实现策略。实验结果如表 2中所示。

根据K-Means聚类算法的执行时间和迭代次数, 可以求出算法单次迭代需要花费的时间。对于600× 600像素的图像, “CPU+内存”、“CUDA+内存”和“CUDA+MMF”三种策略单次迭代花费的时间如图 2(a)所示。可以看出, “CUDA+内存”和“CUDA+ MMF”策略花费的时间远远小于“CPU+内存”策略花费的时间, 只有后者的3%左右, 其中“CUDA+ MMF”策略花费的时间略高于“CUDA+内存”策略花费的时间。

对于6000×6000像素的SPOT-5整景图像, “CUDA+内存”、“CUDA+MMF”和“CUDA+分块”三种策略单次迭代花费的时间如图 2(b)所示。可以看出“CUDA+内存”和“CUDA+MMF”策略花费的时间基本相同, 并且远远小于“CUDA+分块”策略花费的时间, 也只有后者的3%左右。由此可见MMF技术确实能够极大地提高K-Means聚类算法的时间效率。

图 3中显示了“CUDA+内存”、“CUDA+MMF”和“CUDA+分块”3种策略占用的内存数量。对于600×600像素的小图像, “CUDA+MMF”和“CUDA+分块”策略占用的内存略小于“CUDA+内存”策略占用的内存; 对于3 000×3 000像素的中等大小图像, “CUDA+MMF”策略和“CUDA+分块”策略占用的内存明显小于“CUDA+内存”策略占用的内存, 只有后者的25%左右; 对于6 000×6 000像素的大图像, “CUDA+MMF”和“CUDA+分块”策略占用的内存远小于“CUDA+内存”策略占用的内存, 只有后者的10%左右。由此可见, 遥感图像的尺寸越大, 使用分块处理方法的“CUDA+MMF”和“CUDA+分块”策略对内存数量的节省越明显。

从图 3中还可以看出, 使用“CUDA+内存”策略对6 000×6 000像素的137MB的SPOT-5多光谱图像进行处理时, 占用的内存数量已经超过800MB。如果需要处理更大数据量的遥感图像, 可能会出现内存不足的问题。使用“CUDA+MMF”或“CUDA+分块”策略, 当图像数据量从34MB(3000×3000像素)增到137MB时, 算法占用的内存数量基本保持在70MB左右不变, 因此这两种策略都能够处理特大数据量的遥感图像。

通过上面的讨论可以看出, “CUDA+内存”策略运行速度最快, 但是占用内存数量较大; “CUDA+分块”策略占用内存数量较小, 但是运行速度较慢; “CUDA+MMF”策略运行速度与“CUDA+内存”策略相近, 占用内存数量与“CUDA+分块”策略相同, 显然它是最优的实现策略。“CPU+内存”策略不仅运行速度较慢, 而且占用内存数量较大, 通常只适用于处理较小的遥感图像。

本文还研究了瓦片大小对“CUDA+MMF”实现策略执行时间的影响, 总共做了5组对比试验。每一组实验中瓦片的大小分别设定为128×128、256×256、512×512和1024×1024四种尺寸, 占用的内存数量如表 3所示。可以看出, 随着瓦片尺寸变大, “CUDA+ MMF”策略的执行时间呈下降趋势, 但是下降幅度很小, 通常情况下瓦片大小设定为512×512像素或256×256像素即可。

4 结论

本文提出了一种结合CUDA并行计算技术、分块处理方法和内存映射文件技术的遥感图像处理算法优化方法。内存映射文件技术的引入不仅保留了分块处理方法内存占用少的优点, 而且图像文件访问速度与直接内存访问几乎相同。CUDA并行计算技术可以在不增加计算资源的条件下极大地提高遥感图像的处理速度。本文介绍了K-Means聚类算法的“CUDA+内存”和“CUDA+MMF”两种优化策略的实现方法, 并且分析了各种实现策略的时间复杂度和空间复杂度。实验结果表明, 与传统K-Means聚类算法相比, “CUDA+MMF”实现策略能够提高计算速度30多倍, 而占用的计算机内存数量仅为原有算法的10%左右。此外, 图像分块时瓦片的大小对算法性能影响不大。