GPU在遥感影像数据融合算法中的应用

　　1 引言

　　随着现代遥感技术的发展，可获取的遥感影像数据越来越丰富。遥感技术应用的主要障碍不再是影像源的不足，而是如何从中提取更丰富、有用和可靠信息。多源遥感图像数据融合技术是对地观测信息获取与处理技术的重要研究内容，是解决海量遥感图像数据问题的一个有效途径。在处理大数据量的遥感影像数据时，处理速度成为一种融合算法能否普遍适用的一个关键因素。IHS变换可以提高结果图像的地物纹理特性，增强其空间细节的表现能力，因此，在遥感图像融合中，IHS变换法也是常用、基本的方法。但IHS变换融合算法比较复杂，在CPU的串行编程模式中实现较为费时。

　　图形处理单元（GraphicsProcessingUnit,GPU）由于具备强大的并行处理能力和极高的存储器带宽而被抽象成一个“流处理器”，用于科学运算、数据分析、线性代数、流体模拟等需要大量重复的数据集运算和密集的内存存取的应用程序。CPU本质上是一个标量计算模型，计算单元偏少，主要针对复杂控制和低延迟而非高带宽优化。GPU比CPU更适用于流处理计算。

　　2 IHS变换融合算法模型

　　图像处理中经常应用的彩色坐标系统（或称彩色空间）有2种：（1）由红（R）、绿（G）、蓝（B）三原色组成的彩色空间，即RGB空间；（2）IHS，即亮度（Intensity）、色调（Hue）、饱和度（Saturation）。RGB空间和IHS空间相互转换的处理过程称为IHS变换。IHS空间中的3分量I,H,S具有相对独立性，可分别对它们进行控制，并能够准确定量地描述颜色特征。

　　IHS变换的步骤如下：

　　（1）将原始灰度值R0,G0,B0进行归一化R=R0L,0G=GL,0B=BL其中，L为图像的灰度级数。

　　（2）正变换，计算I,H,S的值当R=G=B时，S=H=0,I=R；否则：

　　当G≥B时，H=h；当G<B时，H=2π？h。

　　（3）直方图匹配，将高分辨率全色影像经过灰度拉伸，使其灰度的均值与方差和IHS空间中I分量图像一致。

　　（4）反变换，将拉伸过的高分辨率图像作为新的I分量代入计算，求原始空间的RGB值，反变换如下：

　　IHS变换融合的流程如图1所示。

　　3 IHS变换的GPU实现

　　基于GPU编程的基础是流和核。在流式编程模型中，所有的数据都表现为流。流的定义是具有相同数据类型数据的有序集。核操作输入流上的所有元素而不是独立的元素，核对输入流的每个元素用一个函数进行求值。流式模型用多种方式实现了高效计算，其中重要的一点是流利用了应用程序中的并行性。由于核操作整个流，流元素可以用数据并行硬件并行地处理，有很多元素的长流使数据级的并行性有很高的效率。由此可以看出GPU是利用各个流元素之间的独立性并行地处理多个数据来实现加速。IHS变换融合算法对遥感影像数据中各个像素点的处理是相对独立的，因此，可以将其映射为GPU上的流式计算。

　　3.1 并行的输入流

　　进入GPU的数据分为2类：

　　（1）描述几何形状的顶点数据；

　　（2）片元数据，片元数据通常从纹理中获得。顶点处理器接收顶点数据并转化为屏幕上的窗口坐标，这些顶点组合成几何形状，在光栅化过程中，形状所覆盖的像素被构造成许多片元，由顶点所构成的多变性内部片元可在这一步骤中计算出来。片元传入片元处理器进行并行计算并被终写入帧缓冲区。

　　相对于顶点处理器而言，GPU往往具有更多个片元处理器，因此，片元级的向量计算可获得更高的并行处理能力。

　　本文利用4个顶点生成一个矩形，经过光栅化过程生成片元值。为保证所生成的片元与所处理的遥感图像数据一一对应，将矩形大小定义为纹理大小而不是图像本身的大小。

　　本文实验处理的是经过配准好的高空间分辨率全色影像和低空间分辨率的RGB3波段多光谱影像，假设每幅影像大小为M×N。为充分利用GPU中的RGBA4个通道的计算资源，须对输入数据进行打包或分组。本文将高空间分辨率全色影像和低空间分辨率的RGB3波段多光谱影像分别载入4块纹理之中，尽量减少纹理中空的区域。

　　这样即可将各影像每4个中相邻的像素分为一组，将它们打包进一个RGBA纹理之内，过程如图2所示。

　　为保证数据处理的，本文将纹理内部格式定义为GL_FLOAT_RGBA32_NV，即每个像素中的4个分量各占1个32位浮点空间，将纹理类型设置为GL_TEXTURE_RECTANGLE_ARB，纹理尺寸的大小是任意的，不必要求是2的n次方，纹理坐标不要求单位化，范围不会被限定在0~1之间。

　　3.2 渲染到纹理技术

　　在传统的GPU渲染流水线中，每次渲染运算的终结束点就是帧缓冲区。帧缓冲区其实是显卡内存中的一块，保存在该内存区块中的图像数据会实时地在显示器上显示出来。

　　使用传统的GPU渲染管线作图像处理须将帧缓冲区对应的窗口大小设置为与纹理大小严格一致，这就会把所能处理的图像大小限制在一定大小的范围之内（图像打包后的纹理大小不能超过屏幕大小），否则就会因为系统对图像数据的重采样而造成失真。根据显示器设置的不同，帧缓冲区可以取得32位的颜色深度，即红、绿、蓝、alpha4个颜色通道共享这32位的数据，每个通道占8位。这对于显示器来说已足够，但如果要在浮点数字下工作，用8位来记录一个浮点数，其数学远远不够。为解决这2个问题，本文采用FBO扩展将数据直接渲染进纹理。现在消费级的显卡支持的纹理一般是4096×4096，的显卡会更大，范围远远超过了屏幕窗口的范围，同时在纹理中可将RGBA4个通道的每个通道都设置为32位的浮点型数据，也比帧缓冲区中的高得多，而且这种离屏渲染的模式是通过硬件加速的。

　　在GPU计算中，一个纹理对象不能同时被读写。根据IHS融合算法，本文将一个FBO与6个纹理对象进行绑定（不同的显卡对每个FBO支持的绑定点数目不一样，可以用GL_MAX_COLOR_ATTACHMENTS_EXT进行查询，如果所用显卡所支持的绑定数为4，则须建立2个FBO，每个FBO绑定3个纹理对象），其中3个纹理用来保存正变换获得的I,H,S值，另外3个纹理用来保存反变换获得的R,G,B值，纹理格式和大小都和输入纹理一致。

　　3.3 多渲染目标

　　从IHS融合算法可以看出，正变换同时计算出了I,H,S3个变量的值，反变换同时计算出了R,G,B3个变量的值，传统的渲染方法只能绘制一张纹理，要绘制多张纹理需要不停地切换渲染对象，还须设计不同的片元着色程序多次渲染，这种只绘制一张纹理的方法效率不高。本文采用绘图缓冲扩展实现了多渲染目标，在IHS的正、反变换中分别将I,H,S3个变量及R,G,B3个变量的值分别并行的渲染输出到3个纹理对象之中。

　　3.4 计算过程

　　按3.1节中的方法分别将低空间分辨率的RGB3波段多光谱影像载入3幅纹理中，创建FBO，将其与6个帧缓冲区进行绑定，这6个帧缓冲区的绑定点分别设定为GL_COLOR_ATTACHMENT0_EXT,GL_COLOR_ATTACHMENT1_EXT,?,GL_COLOR_ATTACHMENT5_EXT，每个帧缓冲区对应1个输出纹理。

　　正变换由一个片元着色程序在单遍内实现。本文的着色程序通过OpenGL着色语言（GLSL）实现，将低空间分辨率的RGB3波段多光谱影像纹理数据作为输入纹理，用绘图缓冲扩展将渲染输出的颜色缓冲区个数设置为3，对应GL_COLOR_ATTACHMENT0_EXT,GL_COLOR_ATTACHMENT1_EXT,GL_COLOR_ATTACHMENT2_EXT这3个缓冲区，分别用它们来保存I,H,S3个变量的值。本文根据IHS变换融合算法模型设计了正、反变换2个片元着色程序。

　　第1次渲染用正变换的片元着色程序，目标缓冲区将获得I,H,S3个分量的值。对于在输出统计或量度时须将数据散射保存的算法，并不合适在GPU中处理[3]，因此，本文将直方图匹配的工作交给CPU处理。将保存I分量的缓冲区数据到系统内存，在CPU中将高空间分辨率全色影像与其做直方图匹配，因为I中的数据在正变换过程中经过单位化，所以将I中的数据还原为原始的灰度级，匹配完成后将经过直方图拉伸的高空间分辨率全色影像单位作为新的I分量载入纹理，同时将保存H,S分量的缓冲区作为源缓冲区（将其中的纹理数据与新的I分量数据一起作为第2次渲染的输入流），将GL_COLOR_ATTACHMENT3_EXT,GL_COLOR_ATTACHMENT4_EXT,GL_COLOR_ATTACHMENT5_EXT作为目标缓冲区。第2次渲染用反变换的片元着色程序，将保存变换结果R,G,B值的目标缓冲区中的数据到系统内存中完成IHS融合。整个计算过程如图3所示。

　　4 实验分析

　　实验在一台装有3GHzPentium（R）4和NVIDIAGeForce8800GTS的DELLPC上进行，显卡具有320MB的显存，实验中的真实数据是同一地区的LandsatTM第8波段影像（分辨率为15m）与第2,4,7这3个波段影像（分辨率为30m），实验选取了配准后大小一致的影像，分4组，各组影像大小分别为：256×256,512×512,1024×1024,2048×2048，单位为像素。基于GPU的IHS融合算法使用OpenGLAPI，利用OpenGLShadingLanguage（GLSL）的片元着色程序加以实现，基于CPU的IHS融合算法在MicrosoftVisalC++6.0中实现。

　　效率评价采用算法执行时间进行度量，文中基于GPU的IHS融合算法的实现主要分为以下几步：

　　（1）内存中的低空间分辨率多光谱影像数据载入纹理；

　　（2）在GPU中用正变换的片元着色程序，获得I,H,S3个分量的值；

　　（3）I分量到系统内存；

　　（4）在CPU中将高空间分辨率全色影像与I分量做直方图匹配；

　　（5）将经过直方图匹配的高空间分辨率全色影像载入纹理，替代I分量；

　　（6）在GPU中用反变换的片元着色程序，获得R,G,B的值；

　　（7）将R,G,B值到主内存。将这7步所耗的时间分别记为t1,t2,?,t7。表1为t1~t7在不同影像数据量下的执行时间（括号内为其所占的百分比）。

　　从表1可以看出，t4（在CPU中执行）基本上随图像数据量的增加倍数而相应地成倍增加，而t2和t6（在GPU中执行）的增加速度却较慢，而且与图像数据量增加的倍数没有对应关系。从t1~t7所占的百分比可以看出t4随着图像数据量的增加其所占的比例越来越大，而t2和t6所占的比例越来越小。

　　本文的IHS变换融合算法与基于CPU的IHS变换融合算法在不同数据量下的时间比较如图4所示。

　　从图4可以看出，当图像数据量较小时，基于GPU的算法不一定占优，而当图像数据增大时，基于GPU的IHS融合算法的优势越来越明显，在图像大小达到2048×2048时，基于GPU的算法比基于CPU的算法速度提高超过4倍。这是因为GPU的流处理机制是对输入流的每个元素用同一个函数并行求值的过程，有很多元素的长流使数据级的并行性有很高的效率，所以随着图像数据量的增大，GPU对数据并行处理的优势会更明显。

　　5 结语

　　遥感影像融合往往涉及到大数据量、过程复杂的运算，存在一定的并行性。随着GPU的飞速发展，的图形处理器能够支持32位的颜色通道和更多的显存空间及更复杂的渲染程序，这是实现更多的复杂计算的基础。但基于GPU的程序设计比较复杂，而且有很多限制条件，缺乏像CPU那样的灵活性，这使其应用受到一定的限制。NVIDIA公司的CUDA架构是针对GPU的C语言开发环境，简化了程序设计的复杂程度，可以充分应用流处理单元强大的浮点运算能力，解决复杂的科学运算问题，这也为解决大数据量的处理问题带来了新变革。