快速检索        
  武汉大学学报·信息科学版  2018, Vol. 43 Issue (7): 1000-1007

文章信息

方留杨, 王密, 潘俊
FANG Liuyang, WANG Mi, PAN Jun
CPU和GPU协同的多光谱影像快速波段配准方法
CPU/GPU Cooperative Fast Band Registration Method for Multispectral Imagery
武汉大学学报·信息科学版, 2018, 43(7): 1000-1007
Geomatics and Information Science of Wuhan University, 2018, 43(7): 1000-1007
http://dx.doi.org/10.13203/j.whugis20160218

文章历史

收稿日期: 2017-03-27
CPU和GPU协同的多光谱影像快速波段配准方法
方留杨1 , 王密2 , 潘俊2     
1. 云南省交通规划设计研究院陆地交通气象灾害防治技术国家工程实验室, 云南 昆明, 650041;
2. 武汉大学测绘遥感信息工程国家重点实验室, 湖北 武汉, 430079
摘要:随着遥感影像数据量的飞速增长,传统的串行波段配准方法已无法满足大数据多光谱影像的实时配准需求。针对该问题,提出了一种CPU和GPU协同的多光谱影像快速波段配准方法。首先进行计算量和并行度分析,将同名点匹配和微分纠正映射至GPU执行,仿射变换系数拟合仍驻留在CPU执行。其次通过核函数任务映射和基本设置,使算法步骤在GPU上可执行,并设计了3种性能优化方法(访存优化、指令优化、传输计算堆叠),进一步提高了波段配准的执行效率。在NVIDIA Tesla M2050 GPU和Intel Xeon E5650 CPU组成的实验平台上,对遥感26号卫星多光谱影像的实验表明,使用该方法加速后的波段配准执行时间仅为3.25 s,与传统串行方法相比,加速比达到了32.32倍,可以满足大数据多光谱影像的近实时配准需求。
关键词CPU和GPU协同     波段配准     计算量和并行度分析     核函数任务映射     性能优化    
CPU/GPU Cooperative Fast Band Registration Method for Multispectral Imagery
FANG Liuyang1, WANG Mi2, PAN Jun2     
1. National Engineering Laboratory for Surface Transportation Weather Impacts Prevention, Broadvision Engineering Consultants, Kunming 650041, China;
2. State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China
First author: FANG Liuyang, PhD, engineer, specializes in remote sensing data processing and traffic disaster prevention technology.E-mail:fangliuyang@stwip.com
Foundation support: The National Natural Science Foundation of China, No.41601476;Research Project of Broadvision Engineering Consultants, No.ZL-2015-03
Abstract: With the rapid increase of data size of remote sensing images, the traditional serial band re-gistration method cannot meet the demand for real-time processing of big-data multispectral images. Therefore, a CPU/GPU cooperative fast band registration method for multispectral imagery is proposed in this paper. Firstly, the computational amount and degree of parallelism are analyzed; point matching and differential rectification are ported to GPU to execute while the affine transformation parameter is still calculated on CPU. Secondly, kernel task assignment and basic settings are made to ensure the two above GPU steps executable. Moreover, three performance optimization methods, including memory access optimization, instruction optimization and transmission/computation overlap, are designed to further improve the efficiency of band registration. The experimental results based on NVIDIA Tesla M2050 GPU and Intel Xeon E5650 CPU show that the running time of YG-26 multispectral image band registration is only 3.25 s with our method, which got a speedup ratio of 32.32 compared with the traditional CPU serial method. The proposed method can provide quasi-real-time processing capability for multispectral imagery with big data size.
Key words: CPU/GPU cooperation     band registration     calculation and parallelism analysis     kernel task assignment     performance optimization    

波段配准是多光谱影像处理的核心环节, 对几何定位精度、地物分类准确性等方面都会产生显著的影响[1]。目前常用的波段配准方法为基于影像灰度信息的小面元微分纠正法(以下简称微分纠正法), 其基本思想是先在参考波段上提取一定数量的配准控制点, 并根据控制点的分布将影像划分为若干小面元, 然后在目标波段上按照一定搜索策略匹配同名点, 通过同名点坐标拟合各面元变换公式, 再对各面元进行坐标变换和重采样以获取配准后的影像。

一方面, 由于微分纠正法需要对参考波段和各目标波段进行同名点匹配、面元变换公式拟合、逐像素坐标变换和重采样, 因此计算量很大; 另一方面, 随着对地观测技术的飞速发展, 各类大幅宽、高分辨率光学卫星相继发射, 获取的数据量呈几何级数增长。因此, 传统的串行波段配准方法已无法适用于大数据光学卫星影像的近实时处理需求, 亟需研究新的快速波段配准方法, 以提高光学卫星影像的波段配准效率。近年来, 图形处理单元(graphics processing unit, GPU)已成为当前计算机系统中具备高性能处理能力的部件[2]。GPU计算性能高、体积小、功耗低、软硬件架构通用, 已广泛应用于测绘及相关学科高性能数据处理任务和算法中。典型的算法和应用如正射校正[2]、影像复原[3]、LiDAR点云滤波[4]、SAR影像分割[5]、GIS矢量地图处理[6]等, 均通过由NVIDIA提出的统一计算设备架构(compute unified device architecture, CUDA)映射至GPU进行处理, 处理效率得到了显著提高。在影像配准方面, 陈茜等[7]提出了一种基于GPU的星图配准算法并行程序, 将星点检测和星点坐标计算移植至GPU处理; 朱智超[8]针对多分辨率红外影像与可见光影像配准, 使用CUDA对特征点提取和尺度空间图像构建进行了加速及优化, 但文献[7-8]并没有涉及到GPU性能优化方面的讨论; 周海芳等[9]选取了基于互信息小波分解配准算法进行GPU并行设计, 将存储优化策略应用于遥感图像GPU配准程序, 加速比最大达到了19.9倍, 但未充分进行性能优化, 算法性能仍有进一步提升的空间; 徐如林[10]提出了一种基于GPU的遥感图像全局配准处理算法, 该算法进行了存储和同步优化, 但实验过程中最大仅选取2 048×2 048像素的影像进行测试, 未对大数据遥感影像的配准效率进行验证。虽然上述文献均属于广义的影像配准范畴, 但与多光谱影像波段配准并不完全一致。实际上, 波段配准计算量密集, 具有先天并行性, 非常适合在GPU上处理。

本文基于目前主流的CUDA并行化架构, 提出了一种CPU和GPU协同的多光谱影像快速波段配准方法, 通过对波段配准的各步骤计算量和并行度进行分析, 将其合理映射至CPU和GPU上, 采用多种方法系统全面地优化性能, 显著提高了大数据多光谱影像波段配准的效率。

1 快速波段配准原理与方法 1.1 CPU/GPU协同原理

NVIDIA公司在2007年推出了CUDA[11-12], 其包括GPU通用计算的完整解决方案。CUDA提供了一种串行-并行混合的CPU/GPU协同模型, 即CPU负责处理串行逻辑任务, 包括数据准备、参数设置、显存分配和初始化等; 此外一些计算量不大、并行程度不高的计算任务通常也由CPU处理; GPU则负责处理计算量大、并行程度高的任务。

1.2 CPU/GPU协同的波段配准方法 1.2.1 计算量和并行度分析

波段配准的具体步骤为:

1) 同名点匹配。选定某一波段作为参考波段, 首先确定格网间隔, 将各格网中心点作为配准控制点, 然后确定搜索范围, 对除参考波段以外的其余各目标波段进行相关系数匹配。根据多光谱相机偏移参数消除波段间初始偏移后, 波段间已大致对齐, 故搜索范围通常较小。此外, 在假设影像内部变形为线性的前提下, 对相关系数低于一定阈值的同名点, 可以利用周围同名点的坐标偏移平均值作为该同名点的坐标偏移值, 对匹配结果进行修正优化。

2) 仿射变换系数拟合。匹配完成后, 对每个由同名点作为角点的四边形面元建立仿射变换多项式:

$ \left\{ \begin{gathered} s = {a_0} + {a_1}x + {a_2}y + {a_3}xy \hfill \\ l = {b_0} + {b_1}x + {b_2}y + {b_3}xy \hfill \\ \end{gathered} \right. $ (1)

式中, sl为目标波段像点坐标; xy为参考波段像点坐标; a0a1a2a3b0b1b2b3为仿射变换模型系数, 该系数可由目标波段和参考波段的同名点坐标拟合进行求取。

3) 微分纠正。根据仿射变换多项式系数, 采用间接纠正法方案, 从参考波段出发, 逐面元、逐像素计算参考波段和目标波段的像点坐标对应关系, 并对目标波段进行重采样, 完成波段配准。理论上讲, 四边形面积越小, 配准精度越高。

在上述3个步骤中, 同名点匹配需要在指定搜索范围内逐像素计算相关系数, 微分纠正需要对影像的所有像素进行坐标变换和重采样, 因此这两个步骤的计算量都较大; 仿射变换系数拟合通过4对同名点坐标实现, 可直接求解, 不需迭代, 计算量较小。另外, 同名点匹配和微分纠正逐格网逐像素独立进行, 并行程度高; 但仿射变换系数拟合受格网在影像中的相对位置影响而有所不同, 并行程度较低。因此, 拟将同名点匹配和微分纠正映射至GPU执行。

1.2.2 GPU并行映射

由于同名点匹配和微分纠正执行逻辑连续, 为减少反复调用核函数带来的额外开销, 将这两个步骤分别映射为函数BR_CorPtsMat和BR_Resample进行实现。基本设置包括:

1) 对于BR_CorPtsMat, 每个线程块负责一个同名点的相关系数匹配, 线程块中的每个CUDA线程负责搜索区域内一个像素的相关系数计算; 对于BR_Resample, 每个CUDA线程负责对一个像素进行微分纠正。

2) 将BR_CorPtsMat的CUDA线程块大小设置为与搜索区域大小一致, 此时每个线程块对一个格网控制点进行匹配; 将BR_Resample的CUDA线程块大小设置为256(16×16)。

3) 为减少对全局存储器的多次访问, 使用程序默认的L1 cache缓存模式。

4) 为减少使用循环指令带来的开销, 采用程序默认的循环展开模式。

1.2.3 GPU性能优化

性能优化遵循“渐进式优化”原则, 主要从访存优化、指令优化、传输计算堆叠[13]这3个方面展开, 即首先通过访存优化和指令优化将波段配准的两个核函数(BR_CorPtsMat和BR_Resample)的性能调至最优, 在此基础上再通过传输计算堆叠将整个波段配准算法的全流程性能调至最优。

1) 访存优化

使用3种方法(线程多元素重访、线程块大小和形状最优选取、存储层次性访问)进行访存优化以使存储吞吐量最大化。线程多元素重访策略利用GPU上多次显存访问可被流水化以及线程寻址可重用等特点, 令一个GPU线程同时处理多个元素, 以提高GPU线程的执行效率, 但由于同名点匹配和微分纠正的处理流程和线程组织方式等存在较大差异, 线程多元素重访的性能提升程度并不相同, 且不会随着线程处理元素数量的增多而无限增长, 因此需要通过实验确定各线程处理元素的具体数量, 使算法执行效率达到最优。

在选择线程块大小时, 需满足3个准则[14]:

准则1  线程块大小最好是线程束数量32的整数倍, 以避免在未完全展开的线程束上进行无意义的操作。

准则2  线程块大小应至少在64以上。

准则3  线程块最大不能超过1 024。

此外, 由于线程块的宽度会影响全局存储器的合并访问效率, 因此将各线程块的宽度从1逐渐增加到最大, 并记录核函数的执行时间, 以寻找核函数执行时间最短时所对应的线程块宽度。

在基本实现中, 核函数采用了默认的L1 cache缓存模式缓存数据, 以减少对全局存储器的访问。实际上, 通用计算GPU提供了两级缓存(L1 cache和L2 cache), 通过调整相应的编译选项, 可以屏蔽L1 cache并启用L2 cache缓存数据。此外, BR_CorPtsMat进行相关系数匹配时, 搜索区域内各匹配窗口像素值被反复使用, 因此可将其载入共享存储器; BR_Resample进行影像重采样时, 数据内插具有明显的二维局部性访存特点, 该特点特别适合于使用纹理存储器进行优化。表 1总结了两个核函数的缓存模式。

表 1 核函数缓存模式 Table 1 Illustration of Kernel Caching Mode
核函数 缓存模式
L1 cache L2 cache 共享存储 纹理存储
BR_CorPtsMat
BR_Resample

2) 指令优化

GPU主要负责对计算量大且并行程度高的任务进行处理, 不涉及过多的控制和逻辑操作, 因此其片上控制单元较少。一些高级语言中常用的逻辑操作, 如分支语句、循环语句等, 在GPU上的执行效率较低。在两个核函数中, 由于线程格网大小不一定与待处理影像大小完全相等, 因此在边界像素处理时会引入逻辑分支, 但该分支属于线程块级分支, 不属于线程束内部分支, 不会造成分支语句的重复遍历执行, 故不会带来额外的性能损失。但是采用线程多元素重访策略对算法进行性能优化后, 每个线程处理的元素数可能多于1, 此时会引入循环语句, 为提高算法的执行效率, 需对循环进行展开, 将循环语句变为顺序语句执行。实际上, GPU编译器对程序进行编译时会自动识别循环语句并将其展开, 但是展开的行为(如是否进行展开、具体展开次数)依据GPU型号和程序的不同而有所差别。此时可以使用预编译指令(#pragma unroll)控制循环展开的行为, 并寻找核函数性能最优时所对应的循环展开数量。

3) 传输计算堆叠

目前主流的GPU内部设计了独立的数据拷贝引擎和核函数执行引擎, 可使用流机制将数据拷贝和核函数执行操作分离至各自引擎执行。CUDA流标识了一组GPU操作队列, 该队列中的操作按指定方式执行, 主要包含内存到显存数据拷贝、若干次核函数执行、显存到内存数据拷贝等。在CUDA程序启动后, GPU底层硬件会将数据拷贝操作映射至数据拷贝引擎, 并将核函数执行操作映射至核函数执行引擎。若程序同时指定了多组流操作, 当一组流进行核函数执行操作时, 另一组流可以进行数据拷贝操作, 从而实现传输计算堆叠, 提高GPU的硬件使用效率。

综上所述, 性能优化后CPU和GPU协同的多光谱影像快速波段配准流程如图 1所示。

图 1 性能优化后CPU和GPU协同执行流程 Figure 1 CPU/GPU Execution Flowchart After Performance Optimization
2 实验与分析 2.1 实验数据及平台

选取3景遥感26号卫星标准景多光谱影像进行实验。影像大小为12 288×12 288像素, 10 bit量化, 由蓝、绿、红、红外4个波段组成, 数据量大小约为1.12 GB。

选择NVIDIA Tesla M2050 GPU和Intel Xeon E5650 CPU构建实验平台, 其中NVIDIA Tesla M2050 GPU基于Fermi架构构建, Intel Xeon E5650 CPU包括6个处理核心, 主频2 660 MHz, 内存大小为32 GB。

2.2 算法性能

为定量评价算法性能, 定义两个性能评价指标, 即加速比和性能提升比。加速比表示GPU映射后算法速度与CPU串行算法相比的提升程度, 即:

$ s = {T_{{\text{CPU}}}}/{T_{{\text{GPU}}}} $ (2)

式中, s为加速比; TCPU为CPU串行算法的执行时间; TGPU为GPU映射后算法的执行时间, 该时间包括3个部分:①CPU系统内存和GPU设备显存之间的数据传输时间; ②GPU核函数的执行时间; ③显存分配、核函数启动等步骤的执行时间。

此外, 采用性能提升比表示某一性能优化前后算法的效率提升程度, 即:

$ p = \frac{{({T_{{\rm{bef}}}}-{T_{{\rm{aft}}}})}}{{{T_{{\rm{aft}}}}}} \times 100\% $ (3)

式中, p为性能提升比; Tbef为性能优化前算法的执行时间; Taft为性能优化后算法的执行时间。

首先, 在Intel Xeon E5650 CPU上, 对遥感26号卫星多光谱影像进行波段配准。同名点匹配使用基于相关系数的灰度匹配策略, 搜索范围设定为5×5, 格网大小为500×500, 微分纠正重采样为双线性内插方式。表 2记录了波段配准各步骤的CPU执行时间和时间占比。其中, 执行时间皆取3景多光谱影像执行时间的平均值。

表 2 波段配准各步骤CPU执行时间和占比 Table 2 CPU Running Time and Its Corresponding Percentage for Each Step of Band Registration
算法步骤 CPU执行时间/s 时间占比/%
同名点匹配 33.27 31.67
仿射变换系数拟合 0.048 0.05
微分纠正 71.72 68.28
合计 105.04 100

表 2可以看出, 波段配准的3个步骤中, 同名点匹配和微分纠正占据了整个算法执行时间的绝大部分(31.67%+68.28%=99.95%); 仿射变换系数拟合的时间很短, 仅占全部执行时间的0.05%, 这是因为:①设定的格网尺寸较大(500×500)且多光谱影像较小(12 288 × 12 288), 因此需要求解仿射变换系数的格网数量较少; ②仿射变换为线性变换, 变换系数可直接求解, 不需迭代, 故该步骤计算量不大。因此, 将仿射变换系数拟合驻留在CPU执行, 同名点匹配和微分纠正则映射至GPU执行。

表 3列出了GPU并行映射后算法各步骤的执行时间(其余耗时项包括显存分配、核函数启动、计时开销等)。从表 3可以发现, CPU和GPU之间存在多次数据往返(分别用CPU→GPU和GPU→CPU表示)。其中, CPU→GPU数据传输时间主要由待配准影像的传输时间决定, GPU→CPU数据传输时间主要由配准后影像的传输时间决定。由于GPU计算性能较高, 波段配准的处理效率有了显著提升, 整体加速比达到了29.76倍; 但核函数BR_CorPtsMat的加速比仅为23.27倍, 远低于BR_Resample的311.83倍, 造成加速比悬殊的主要原因是考虑到经过初始偏移调整后, 多光谱影像各波段间的坐标差已经很小, 因此本文将相关系数的搜索范围设置较小(仅为5×5)。在初始设置中, BR_CorPtsMat一个线程块负责一个格网控制点的搜索匹配, 即线程块的大小也仅为5×5, 小于一个完整的线程束中包含的线程数量, 导致线程束上出现一些无效操作, 使得流式多处理器上的线程占有率较低。此外, 每个线程计算完各自像素点的相关系数后, 还需要设置一个线程对相关系数进行比较, 选择并输出最终结果, 其余线程在这个时间段内处于等待状态, 因此这两个因素在一定程度上影响了BR_CorPtsMat的计算效率。

表 3 GPU映射后波段配准各步骤执行时间和加速比 Table 3 GPU Running Time and its Speed-up Ratio for Each Step of Band Registration
算法步骤 GPU执行时间/s 加速比
CPU→GPU传输待配准影像 0.45 -
BR_CorPtsMat 1.43 23.27
GPU→CPU传输匹配点坐标 3.01×10-6 -
仿射变换系数拟合 0.048 -
CPU→GPU传输仿射变换系数 1.35×10-5 -
BR_Resample 0.23 311.83
GPU→CPU传输配准后影像 0.21 -
其余耗时项 1.162 -
合计 3.53 29.76

利用§1.2.3中的优化策略对核函数进行性能优化。首先使用线程多元素重访策略, 令一个CUDA线程处理多个元素, 相应的核函数执行时间如图 2所示。从图 2可以看出, 当BR_CorPtsMat和BR_Resample的每个线程分别处理1、128个元素时, 核函数执行时间最短, 分别为1 430.1 ms和191.5 ms, 即线程多元素重访策略提高了BR_Resample的执行效率, 但却降低了BR_CorPtsMat的执行效率。这是因为BR_CorPtsMat的每个线程块负责一个格网控制点的搜索匹配, 由于搜索范围设定为5×5, 当每个CUDA线程处理一个像素时, 线程块大小也仅为5×5=25, 尚不足一个完整的线程束大小(1个线程束包含32个CUDA线程)。若使用线程多元素重访策略, 令一个CUDA线程处理多个像素, 线程块将会更小; 当每线程处理像素数为32以上时, 线程块甚至缩小至1×1, 无法充分利用GPU的硬件计算性能。

图 2 线程多元素重访优化后核函数执行时间 Figure 2 Kernel Running Time with Multiple Elements After per Thread Access Optimization

使用不同的线程块大小、形状和缓存模式时, BR_Resample的执行时间如图 3所示。其中柱状子图记录了核函数的执行时间(该执行时间为线程块宽度从1增加到最大时的最短执行时间), 折线子图记录了对应的线程块宽度。从图 3中可以发现:

图 3 线程块大小、形状、缓存模式变化时BR_Resample执行时间 Figure 3 BR_Resample Running Time with Different Block Sizes, Shapes and Caching Modes

1) BR_Resample的线程块最大为512, 该值由式(4)决定:

$ {t_s} = \left\lfloor {\frac{{{r_s}/{r_t}}}{{32}}} \right\rfloor $ (4)

式中, $\left\lfloor \cdot \right\rfloor $ 表示向下取整; rs代表流式多处理器中的寄存器数量, 对于NVIDIA Tesla M2050 GPU, 该值为32 768;rt代表每个线程使用的寄存器数量, 其值为63;ts代表流式多处理器中的线程最大值, 将已知量代入可计算其值为512。

2) 对于同一种缓存模式, 当线程块逐渐增大时, 执行时间无明显差异。这是因为在线程块逐渐增大的过程中, 线程占有率最低时, 流式多处理器中驻留的线程束数量也足够隐藏访问延迟。容易发现, 当线程块大小为320时, 线程占有率最低, 此时可根据式(5)求出隐藏访问延迟所需要的计算访存比R:

$ R = \left\lceil {\frac{{{f_p}/({b_w}/{b_l})}}{n}} \right\rceil $ (5)

式中, $\left\lceil \cdot \right\rceil $ 表示向上取整; n是线程束数, 当线程块大小为320时, 其值为10;fp表示GPU的双精度浮点计算性能, 其值为515 Gflops; bw为GPU的显存带宽, 其值为148 GB/s; bl为核函数访存时读取数据类型的字节长度。BR_Resample的访存数据类型有两种, 即双精度仿射变换系数和影像像素灰度, 其bl值分别为8和2, 此处为方便讨论, 将该值设为8。由于式(5)中blR的值成正比, 将该值设为8就意味着需要更高的计算访存比才能隐藏访问延迟, 条件限制更为严格。将上述值代入式(5)中, 计算得到R的值为3, 也就是说, 理论上BR_Resample每执行1次访存就必须要进行3次计算才能隐藏访存所带来的延迟。由于坐标变换和重采样的计算量较大, 纠正一个像素通常需要上百次的浮点运算, 远多于访存次数, 因此该要求完全可以达到。

3) 无论使用哪种缓存模式(纹理存储器、L1 cache、L2 cache), 核函数执行效率几乎相同。虽然BR_Resample中存在着明显的二维局部访存模式, 且该模式非常适合于使用纹理存储器进行加速, 但纹理存储器在性能上却并未体现出对于L1/L2 cache的优势。

对于BR_CorPtsMat, 由于其线程块大小固定为5×5, 因此仅记录使用不同缓存模式时核函数的执行时间, 如表 4所示。由表 4可知, 当使用L1 cache作为数据缓存时, 核函数效率最优, 相应的执行时间为1 430.1 ms。

表 4 缓存模式变化时BR_CorPtsMat执行时间 Table 4 BR_CorPtsMat Running Time with Different Caching Modes
核函数名 缓存模式 执行时间/ms
共享存储 1 447.2
BR_CorPtsMat L1 cache 1 430.1
L2 cache 1 683.4

由于BR_Resample使用了线程多元素重访策略。令一个线程对多个像素进行重采样, 即核函数中存在着循环逻辑。可使用预编译指令#pragma unroll进行循环展开, 并分析循环展开的次数对核函数执行效率的影响, 如图 4所示。随着循环展开次数的增加, 核函数的执行时间逐渐下降; 当循环被完全展开时, 执行时间最短, 与不使用#pragma unroll时基本相同, 说明编译器已默认将循环完全展开。

图 4 循环展开时BR_Resample执行时间 Figure 4 BR_Resample Running Time with Loop Unrolling

表 5总结了经性能优化后核函数的最优配置、执行时间和性能提升比。从表 5可知, BR_Resample的性能提升较为明显, 达到了22.9%;而BR_CorPtsMat的性能未发生变化, 这是因为在进行参数基本设置时, 人工指定和编译器选择的一些参数已在一定程度上考虑了性能优化的因素。

表 5 性能优化后核函数最优配置、执行时间和性能提升比 Table 5 The Optimal Configuration, Running Time and Performance Improvement Percentage of Kernel After Optimization
核函数 每线程处理像素数 线程块大小(宽度×高度) 缓存模式 循环展开次数 执行时间/ms 性能提升比/%
BR_CorPtsMat 1 25 (5 × 5) L1 cache 1 1 430.1 0
BR_Resample 128 128 (8 × 16) L1 cache 128 187.2 22.9

在完成核函数性能优化的基础上, 使用流机制实现传输计算堆叠, 进一步提高波段配准的执行效率。表 6为分别创建2、3、4个流时, 波段配准的执行时间和性能提升比。由表 6可知, 当创建2个流时, 执行时间最短, 对应的性能提升比最高。由于流的创建和销毁会带来额外开销, 因此创建过多的流反而会带来性能损失。此外, 由于仿射变换系数拟合未映射至GPU上执行, 因此BR_CorPtsMat和BR_Resample两个核函数对应的流被分开执行。第1个流(即BR_CorPtsMat驻留的流)主要负责堆叠CPU→GPU待配准影像数据传输和同名点匹配, 第2个流(即BR_Resample驻留的流)主要负责堆叠微分纠正和GPU→CPU配准后影像数据传输。由于CPU→GPU待配准影像数据传输和同名点匹配的执行时间差异较大(分别为0.45 s和1.43 s), 因此堆叠后性能提升的效果有限。

表 6 使用流后波段配准执行时间和性能提升比 Table 6 Running Time and Performance Improvement Percentage of Band Registration with Streams
流的数量 执行时间/s 性能提升比/%
2 3.25 7.30
3 3.29 5.99
4 3.35 4.10

最终, 性能优化后的波段配准执行时间为3.25 s, 与传统串行方法相比, 加速比达到了32.32倍, 可以满足大数据多光谱影像的近实时处理需求。

2.3 正确性

将本文方法和传统CPU串行方法处理的3景影像分别输出, 逐像素计算均方根误差(root mean square error, RMSE):

$ {\rm{RMSE}} = \sqrt {\frac{{\sum\limits_{k = 1}^p {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{(V_{(k, i, j)}^{{\rm{GPU}}}-V_{(k, i, j)}^{{\rm{CPU}}})}^2}} } } }}{{p \times m \times n}}} $ (6)

式中, V(k, i, j)GPUV(k, i, j)CPU分别为本文方法和传统CPU串行方法处理的第k景影像索引为(i, j)的像素值; mn分别为影像的行宽度和列宽度; p为影像数量。

经统计, 使用两种方法处理的3景影像的浮点均方根误差为6.23×10-8, 但由于波段配准最后输出的影像像素值为舍入后的整数值, 因此最终输出的影像整数均方根误差为0, 说明了本文方法的正确性。实际上, NVIDIA Tesla M2050 GPU和Intel Xeon E5650 CPU都支持IEEE 754-2008浮点数精度标准(即支持微小浮点数和4种IEEE舍入模式), 舍入前出现微小浮点均方根误差的原因是NVIDIA Tesla M2050 GPU支持混合乘加操作, 该操作执行乘加运算时仅最后舍入一次, 但Intel Xeon E5650 CPU进行乘加运算时会舍入两次(即乘法运算后舍入一次, 加法运算后再舍入一次), 因此其计算精度优于Intel Xeon E5650 CPU的计算精度[15]

3 结语

本文提出了一种基于CPU和GPU协同处理的多光谱影像快速波段配准方法。首先进行核函数任务映射和基本配置, 使该方法在GPU上可执行, 然后通过访存优化(包括线程多元素重访、线程块大小和形状最优选择、存储层次性访问)、指令优化、传输计算堆叠3个性能优化手段, 进一步提高了方法的执行效率。实验结果表明, 在NVIDIA Tesla M2050 GPU和Intel Xeon E5650 CPU组成的实验平台上, 使用本文方法对遥感26号卫星多光谱影像进行波段配准的时间仅为3.25 s, 加速比达到了32.32倍, 可以满足大数据多光谱影像的近实时处理需求。此外, 将经过本文方法与传统CPU串行算法处理后的影像采用逐像素方式计算RMSE, 结果发现舍入后的整数均方根误差为0, 说明了本文方法的正确性。

参考文献
[1] Yu Haiyang, Gan Fuping, Dang Fuxing. An Experimental Analysis of Band to Band Registration Error in High Resolution Satellite Remote Sensing Imagery[J]. Remote Sensing for Land & Resources, 2007, 19(3): 39–42 ( 于海洋, 甘甫平, 党福星. 高分辨率遥感影像波段配准误差试验分析[J]. 国土资源遥感, 2007, 19(3): 39–42. DOI:10.6046/gtzyyg.2007.03.09 )
[2] Yang Jingyu, Zhang Yongsheng, Li Zhengguo, et al. GPU-CPU Cooperate Processing of RS Image Ortho-Rectification[J]. Geomatics and Information Science of Wuhan University, 2011, 36(9): 1043–1046 ( 杨靖宇, 张永生, 李正国, 等. 遥感影像正射纠正的GPU-CPU协同处理研究[J]. 武汉大学学报·信息科学版, 2011, 36(9): 1043–1046. )
[3] Fang L Y, Wang M, Li D R, et al. MOC-Based Parallel Preprocessing of ZY-3 Satellite Images[J]. IEEE Geoscience and Remote Sensing Letter, 2015, 12(2): 419–423 DOI:10.1109/LGRS.2014.2345419
[4] Hu X Y, Li X K, Zhang Y J. Fast Filtering of LiDAR Point Cloud in Urban Areas Based on Scan Line Segmentation and GPU Acceleration[J]. IEEE Geoscience and Remote Sensing Letter, 2013, 19(2): 308–312
[5] Sui H G, Peng F F, Xu C, et al. GPU-Accelerated MRF Segmentation Algorithm for SAR Images[J]. Computers & Geosciences, 2012, 43(2): 159–166
[6] Cheng Boyan, Liu Qiang, Li Xiaowen, et al. Parallel Rasterization of Vector Polygon Based on CUDA[J]. Bulletin of Surveying and Mapping, 2014(11): 97–101 ( 程博艳, 刘强, 李小文, 等. 利用CUDA实现矢量地图栅格化的并行处理[J]. 测绘通报, 2014(11): 97–101. )
[7] Chen Xi, Qiu Yuehong, Yi Hongwei. Parallel Programming Design of Star Image Registration Based on GPU[J]. Infrared and Laser Engineering, 2014, 43(11): 3756–3761 ( 陈茜, 邱跃洪, 易红伟. 基于GPU的星图配准算法并行程序设计[J]. 红外与激光工程, 2014, 43(11): 3756–3761. DOI:10.3969/j.issn.1007-2276.2014.11.044 )
[8] Zhu Zhichao. Research on GPU-Based Multi-resolution Infrared and Visible Image Registration[D]. Nanjing: Nanjing University of Aeronautics and Astronautics, 2011 ( 朱智超. 基于GPU的多分辨率红外与可见光图像配准研究[D]. 南京: 南京航空航天大学, 2011 )
[9] Zhou Haifang, Zhao Jin. Parallel Programming Design and Storage Optimization of Remote Sensing Image Registration Based on GPU[J]. Journal of Computer Research and Development, 2012, 49(1): 281–286 ( 周海芳, 赵进. 基于GPU的遥感图像配准并行程序设计与存储优化[J]. 计算机研究与发展, 2012, 49(1): 281–286. )
[10] Xu Rulin. Study of Parallel Algorithms for Remote Sensing Image Registration Based on GPU and Implement of Application System[D]. Changsha: National University of Defense Technology, 2014 ( 徐如林. 基于GPU的遥感图像配准并行算法研究及应用系统实现[D]. 长沙: 国防科技大学, 2014 )
[11] Qiu Deyuan. GPGPU Programming Technique:From GLSL, CUDA to OpenGL[M]. Beijing: China Machine Press, 2011 ( 仇德元. GPGPU编程技术:从GLSL、CUDA到OpenCL[M]. 北京: 机械工业出版社, 2011 )
[12] NVIDIA. CUDA C Programming Guide, V5. 0[S]. Santa Clara: NVIDIA Corporation, 2012
[13] NVIDIA. CUDA C Best Practices Guide, V5. 0[S]. Santa Clara: NVIDIA Corporation, 2012
[14] Kirk D, Hwu W M. Programming Massively Parallel Processors[M]. 2nd ed. Massachusetts: Morgan Kaufmann Publishers, 2012
[15] NVIDIA. NVIDIA's White Paper of Precision & Performance: Floating Point and IEEE 754 Comp-liance for NVIDIA GPUs[S]. Santa Clara: NVIDIA Corporation, 2012