面向汛旱情监测的遥感影像GPU并行处理算法

赵晓晨，吴皓楠，李林宜，孟令奎

(1.武汉大学遥感信息工程学院,武汉 430079；2.珠江水利委员会珠江水利综合技术中心，广州 510611)

0 引言

在汛旱情监测领域，遥感技术以其长期动态监测的优势，成为强有力的技术手段之一。遥感对地观测技术飞速发展的同时，带来了文件体积过大、处理耗时过长的问题，因此，提高数据处理效率、减少耗时成为亟待解决的问题之一。

汛情及洪涝监测应用的关键在于水体范围及水体变化的信息提取，目前常用归一化差异水体指数(normalized difference water index，NDWI)、改进的归一化差异水体指数(modified NDWI，MNDWI)、增强型水体指数(enhanced water index，WI)等方法提取水体信息，取得了一定成效。旱情监测的主要方法之一是利用植被指数研究不同时期作物长势，经典的指数包括归一化植被指数(normalized difference vegetation index，NDVI)、比值植被指数(ratio vegetation index，RVI)、植被状态指数(vegetation condition index，VCI)等。

统一计算架构(compute unified device architecture，CUDA)的基本思想是驱动程序将密集计算数据划分成细小的区块，映射、加载到图形处理器(graphics processing unit，GPU)中并行执行，尽量开发线程级并行并将这些线程在硬件中动态调度和执行[1-2]。CUDA采用主机程序嵌套内核函数的程序设计模式，主机/设备代码分离的编译模式，数据由主机端传输至设备端，并调用设备端的内核函数派生大量独立线程并行执行。GPU与中央处理器(central processing unit，CPU)协同处理技术是将GPU视为高性能并行计算设备，将CPU称为主机端，把GPU看作CPU的协处理器，避免了从通用算法到应用程序接口(application programming interface，API)的复杂映射，从而直接使用GPU的计算能力，实现了异构环境下对GPU资源的灵活调用。

本文在分析了遥感影像处理原理的基础上，将处理过程划分为若干模块，使用CUDA技术对影像处理方法做并行改进，提出了基于CUDA架构的CPU-GPU协同的遥感影像预处理及相关指数计算并行算法，并选取了典型影像数据对该算法进行验证，证明其有效性。

1 遥感影像预处理及指数算法分析

汛旱情监测应用中的遥感影像处理主要包括辐射校正、几何纠正、遥感指数计算等过程，这里对传统算法原理进行分解分析，并确定可优化模块。

1.1辐射校正

辐射校正的目的是将影像像元量化的、无量纲的像元灰度值定标为辐射亮度值或表观反射率，通常定标为表观反射率的计算表达式为：

L=GainDN+Bias，

(1)

式中：L为大气表观反射率；DN为像元灰度值；Gain为辐射增益参数；Bias为辐射偏置参数，Gain和Bias统称为辐射定标参数。

实现辐射定标首先要读取影像头文件，获取辐射增益、辐射偏置、太阳高度角等元数据信息，其次根据式(1)对各波段影像栅格进行波段运算，考虑对波段运算处理过程进行并行优化。

1.2 几何纠正

影像几何纠正的目的是消除图像的几何畸变，该过程是遥感数据预处理中耗时最长、计算量最大的过程之一，也成为遥感影像快速完成处理的瓶颈。这里采用有理多项式系数(rational polynomial coefficient，RPC)模型对影像纠正(图1)，相关研究表明[3]，RPC模型相对于严格成像模型的拟合误差不超过0.04个像元，均方差小于0.01个像元，因此可以看作对严格成像模型的拟合。

图1 RPC模型几何纠正算法Fig.1 Geometric correction algorithm of RPC model

由图1知，RPC纠正算法的实际开发内容主要包括以下4个模块：①影像读取及批处理任务管理；②RPC正解函数构建；③RPC反解函数构建；④纠正后影像栅格重采样。

前3个模块运行效率主要依赖于计算机存储设备性能，重采样是RPC纠正中耗时最长的步骤。使用传统的串行程序设计方式处理时，对所有像素依次循环处理，造成了极大的流程阻塞，成为算法的最大瓶颈。同时由于每个像素的重采样独立进行，算法具有可并行性，因此考虑作为主要并行化模块进行设计。

1.3 遥感指数计算

1.3.1 汛情遥感监测算法

水体范围变化是汛情及洪涝灾害的主要特征，水体信息的有效提取也一直是水利遥感技术研究的重点。常用的水体光谱指数有NDWI[4]和MNDWI[5]等，具体表达式分别为：

NDWI=(RGreen-RNIR)/(RGreen+RNIR)，

(2)

MNDWI=(RGreen-RMIR)/(RGreen+RMIR)，

(3)

式中，RNIR，RGreen和RMIR分别为近红外、绿光和中红外波段反射率。

1.3.2 旱情遥感监测算法

NDVI，VCI和EVI等遥感指数可用于监测植被长势，进而表征干旱程度。具体表达式分别为：

NDVI=(RNIR-RRed)/(RNIR+RRed)，

(4)

VCI=(NDVI-NDVImin)/(NDVImax-NDVImin)，

(5)

(6)

式中：RRed和RBlue分别为红光和蓝光波段反射率；L为土壤调节参数，L=1；C1为红光大气改正参数，C1=6；C2为蓝光大气改正参数，C2=7.5；G为增益系数。

遥感指数计算在算法上无需直方图统计、栅格索引变更等复杂过程，仅包含不同波段影像栅格数据中同名像元的波段计算，这一过程内在并行性强，因此将波段运算模块作为并行化改造的首要任务。

2 CUDA并行算法设计

遥感影像处理算法并行模式如图2所示，遥感影像预处理算法通常以影像数据的管理与解析作为起始功能模块，将影像元数据和光谱数据栅格分别应用于不同的算法模块。其中元数据主要用于相关模型的构建及参数求解，而光谱数据作为影像信息的主要承载，是包括波段运算、栅格重采样、统计学分析、矢量化等多种算法模块的数据来源，其中波段运算和栅格重采样具有较强的区域并行性，适合进行并行化改造。同时不同功能模块间存在复杂的协同关系，从而使算法理论能够准确地应用于光谱数据本身最终完成产品的输出。

图2 遥感影像处理算法并行模式Fig.2 Parallel mode of remote sensing image processing algorithm

2.1 辐射校正并行算法

根据1.1节的原理分析，辐射校正算法的效率瓶颈集中在波段运算环节。相比于遥感指数计算，辐射定标涉及的波段更多，公式更复杂，分支结构更多，计算密集度更高。从影像头文件中获取辐射偏置、辐射增益、太阳高度角等定标参数，也是算法的重要环节。结合GPU存储方式的分析以及定标参数在算法中定位，实际开发过程中使用__constant__将参数映射至常量内存，提高算法效率。

首先采用GDAL(geospatial data abstraction library)库用于数据解析，并通过注册影像驱动的方式，实现Arc/Info ASCII Grid(asc)，GeoTiff (tiff)，Erdas Imagine Images(img)和ASCII DEM(dem)等不同格式的影像文件读取。

影像数据解析为栅格数组后，以数据集和波段2级数据类型管理，针对异构环境下设备端显存空间有限的问题，采用条带划分、分批次处理模式。这种单一的波段运算不涉及空间运算，仅针对同一像元位置不同波段灰度值，因此对栅格数组逻辑分块可以避免多余的空间信息赋值工作，并且根据GPU的需求弹性调整分块大小。

使用核函数编写波段计算是实现并行处理的核心步骤，实际编写中有以下关键点：①使用__global__函数限定符，主机端负责调用，设备端负责执行；②在核函数中完成栅格数组索引与线程逻辑管理索引的正确映射；③算法公式中涉及的变量选取合适的变量类型确保分配足够的显存空间；③在正确表述算法公式外，使用分支结构处理函数异常值。

其中线程管理索引决定了栅格数组元素在线程逻辑结构中的分配方式，其公式为：

I=(byGx+bx)(BxBy)+tyBx+tx，

(7)

式中：I为栅格数组索引；bx和by分别为线程块在网格中x和y方向上的索引，对应内建变量blockIdx.x和blockIdx.y；Gx为一个网格中x方向上线程块总数，对应内建变量gridDim.x；Bx和By分别为一个线程块中x和y方向上线程总数，对应内建变量blockDim.x和blockDim.y；tx和ty分别为线程在线程区块中x和y方向上的索引，对应内建变量threadIdx.x和threadIdx.y。本文结合二维影像数据自身特点，同样使用二维格网、线程块进行管理。

另外，由于核函数中采用了少量分支结构用于波段运算中异常值的处理，可能造成线程块中执行速率不一致。为此可使用CUDA中内建函数cudaThreadSynchronize()保证同一线程块中各线程状态对齐。

综上，针对波段运算的并行算法设计中，主机端负责批处理任务编排、影像文件读取和输出、影像栅格数组划分等耗时较短的串行任务，设备端负责对影像栅格数据并行执行核函数，完成指数算法及异常值处理。这一设计也是RPC几何纠正并行算法的基础。

2.2 基于RPC模型的几何纠正并行算法

通过1.2节中对基于RPC模型的几何纠正算法的原理分析，结合CUDA程序设计模型设计了CPU-GPU协同的RPC纠正算法，如图3所示。

图3 并行化RPC纠正算法Fig.3 Parallel RPC algorithm

该算法实施步骤如下：①影像解析及参数读取，并存储于主机端内存；②根据RPC模型构建正反解函数，正解函数在主机端调用并执行，且仅用于纠正后影像四至点的结算，反解函数将用于设备端影像并行重采样；③矫正后影像栅格初始化，填充重采样灰度数据，同时对该范围进行逻辑分块，以便在有限的设备端存储空间下完成整张影像的重采样，分块后影像条带范围被传入GPU，在全局内存中初始化地面点坐标数组，根据校正后影像条带范围反解得到相对应的原始影像范围，将该范围内原始影像栅格传入GPU；④RPC反解函数中各项系数随第一个影像块传入GPU，并使用__constant__限定符存储为常量内存，然后将反解函数派生为多线程任务，将分块影像条带范围内的地面点坐标并行解算为对应的影像坐标，并根据此坐标对原始影像条带进行重采样；⑤将纠正后影像条带回传至主机端，合并入纠正后影像栅格。最后使用GDAL库中的Geo-TIFF影像驱动封装为影像产品。其中用于并行坐标反解及重采样主要包括并行纠正核函数、RPC模型填充函数以及RPC反解函数等几类。

2.3 遥感指数并行算法设计

遥感指数并行算法与辐射校正并行算法相似，对读取后的影像波段直接进行数学计算。对于多波段影像，相较于辐射定标过程，指数计算涉及的波段数更少，因此也更为简单。

2.4 访存优化

在2.1节波段运算并行算法设计中提到了对栅格数组进行逻辑划分以避免显存溢出，以及栅格数组在多级线程管理结构中的映射方法，这2项算法设计实现了异构环境下高空间分辨率影像分批并行处理，是提升并行算法效率的关键。因此这里结合GPU硬件参数和相关官方资料进行了实验，并对相关算法参数进行优化。

线程束(warp)作为GPU执行任务时的调度单位，由调度器分配依次进入流处理器簇(SM)执行。SM每次只执行线程块(block)中的一个warp，称作active warp，此时其余warp处于待命状态。当active warp在执行中需要等待时(如访问GPU全局内存)，调度器可以直接切换至其他warp执行。通过这种方式GPU隐藏了多线程的延迟等待，实现了大规模的并行。因此在程序设计时，可以在同一block设置尽量多的warp，避免所有warp均处于等待访存状态。但由于SM内寄存器数量有限，同一个block内warp数量过多可能造成单一线程可使用的寄存器数量过少，影响数据访存效率。同时由于warp中线程(thread)数为固定值，因此同一block内thread总数应设计为这一固定值的倍数，避免warp切换时流处理器(SP)的闲置。通过查阅相关文档获得了实验平台所使用GPU相关硬件参数如表1所示。

表1 GPU硬件参数Tab.1 GPU parameters

3 并行优化实验及分析

3.1 实验平台与数据

GPU硬件及CUDA详细参数如表2。实验使用便携计算机平台，CPU版本为Intel(R)Core(TM)i5-2520M CPU @ 2.50 GHz，GPU为NVIDIA GeForce GT 640 M LE，DDR3内存为16 GB。影像数据Landsat8 Level-1产品，OLI和TIRS传感器，共11个波段，覆盖周期为16 d，扫描宽度为170 km×180 km，单波段影像数据量为100 MB。高分一号(GF-1)L1A产品，共9个波段，PMS传感器幅宽为60 km，空间分辨率为2～8 m，WFV传感器幅宽为800 km，空间分辨率为16m，单景影像数据量为1.2 GB。

表2 GPU硬件及CUDA详细参数Tab.2 GPU and CUDA parameters

3.2 实验参数设置

根据3.1节中的硬件参数，在控制实验平台其他硬件参数不变的前提下，使用1景GF-1 WFV影像(数据量1.63 GB)测试了不同线程块尺寸(以32为基准)下NDVI并行处理算法中波段运算环节处理速度，结果如图4所示。

图4 线程块尺寸实验Fig.4 Different thread block sizes

可见在线程块内最大线程数范围内，线程块尺寸对波段运算耗时的影响并不显著。原因可能在于本算法涉及的数据类型较为单一，各线程对全局内存访存延迟较短。

对于可能影响算法运行效率的另一环节——栅格数组逻辑分块，同样控制了其他变量不变，以影像宽度为基准，分别测试了不同行数的条带划分尺度下波段运算耗时，实验结果如图5所示。

图5 栅格划分尺度实验Fig.5 Grid scale of performance

由图5可以看出，栅格数组的划分尺度的确影响了并行波段运算的效率，更大的条带划分可以减少CPU-GPU间数据交换的频数，缩短并行算法的时间。但条带划分仍需考虑GPU显存容量的上限，结合实验和理论分析可以得出最佳条带划分公式为：

(8)

式中：T为最佳条带宽度；H和W分别为二维影像栅格数组的高度和宽度；n为算法使用的栅格数组数量；M为GPU全局内存;m为影像波段数量。

3.3 实验过程与结果分析

3.3.1 汛旱情监测指数

分别使用基于CUDA的并行算法、传统串行算法处理一景GF-1 WFV影像数据，进行NDVI和NDWI指数产品的生产，对各模块计算耗时统计后分别见图6。以鄱阳湖为例计算光谱指数产品如图7所示。

(a)NDVI指数算法耗时对比 (b)NDWI指数算法耗时对比

(a)NDVI (b)NDWI

对于简单波段运算类光谱指数算法，采用基于CUDA的并行算法能够有效提高波段运算这一瓶颈环节的效率，影像读取和产品输出环节的时间消耗与串行算法基本一致。此外由于波段运算核函数中具体内容与串行算法循环体中函数内容完全一致，因此保证了产品结果的正确性。

3.3.2 辐射校正

分别使用基于CUDA的并行算法、传统串行算法对1景Landsat8影像数据进行辐射定标。由于Landsat8各波段影像使用单独的TIFF文件存储，因此使用各波段处理的平均耗时作为对比。结果如表3所示。表3中可见，由于各波段影像差异和计算机存在性能因素，算法耗时存在一定范围的浮动。从整体来看，波段运算环节处理效率有了大幅提升，相比于串行处理方法，并行算法可节约58.9%的时间，但直方图统计环节仍然有一定的优化空间，算法总耗时得到了一定程度的降低。

表3 传统串行和CUDA并行辐射处理耗时对比Tab.3 Processing time of CPU and CUDA for radiometric correction (s)

(续表)

3.3.3 RPC几何纠正

分别使用基于CUDA的并行RPC纠正算法、传统串行RPC纠正算法，对1景GF-1 WFV影像数据进行几何纠正。在影像重采样环节使用了最邻近像元法和双线性内插法，比较不同重采样算法对并行加速比的影响，如图8和表4所示。根据图8和表4可以发现，基于CUDA的并行重采样算法相比传统串行算法可以实现7～10倍的加速比，且不同重采样方式的加速效果较为接近。使用较为复杂重采样方式总处理时间较长、并行处理占比较大，因此总加速比更高。

(a)最邻近像元重采样 (b)双线性内插重采样

表4 RPC并行校正算法重采样环节加速比Tab.4 Speedup radio of parallel RPC algorithm

4 结论

针对汛旱情监测数据处理包括辐射校正、几何纠正及遥感指数计算等过程耗时过长的问题，面向常见遥感数据如Landsat和GF-1等，分析各类处理过程的共性与性能瓶颈，基于CUDA架构提出了基于功能模块划分的并行处理算法。

通过实验发现，调整核心参数进行算法访存优化，可以有效提升并行运算的效率，实验结果表明该算法在保证成果精度的同时，相比传统串行算法效率有了显著提升。辐射校正与指数计算均涉及的波段计算模块，经过并行优化后，较串行模式可减少70%以上的时间；几何纠正不同重采样方式的加速效果较为接近；通过调整核心参数进行算法访存优化，可以有效提升并行波段运算的效率。本算法针对便携化处理平台而设计，可为汛旱情监测提供切实有效的高性能计算解决方案。