快速变化和高速发展的现代社会，土地利用和覆盖变化一直都是人们关注的话题。在众多的土地动态变化检测方法中，遥感技术以其经济、实用、高效的优势，在变化检测领域受到了广泛应用。然而随着遥感技术的迅速发展，遥感图像的数据量急剧增加，数据结构也逐渐复杂，在遥感图像的读取，显示和变化检测等过程中暴露出许多问题，如内存溢出、读取和存储数据耗时长等。另一方面，为了提高遥感图像变化检测的精度，各种算法层出不穷，也越来越复杂，因而计算量和计算耗时也大大增加。由于CPU受体系结构限制，大部分片上资源都用来做控制和缓存，实际的计算单元十分有限，导致基于遥感图像的变化检测效率低下。为了满足日益增长的计算需求，实时地掌握土地的变化信息，就需要一种快捷有效的图像处理方式来检测土地利用和覆盖的变化情况。

目前用于遥感图像处理的并行技术主要是基于计算机集群的并行和基于加速器的并行。基于计算机集群的并行技术有：MPI、PVM、Express等，其中最受欢迎的是基于消息传递的MPI技术，该技术需要多个计算机节点通过网络链接起来，运算时将任务分配到各个节点上运算，最终结果返回某一个节点。国内外学者对基于集群的变化检测做了一些研究：Mubasher等人(2012)将基于背景模型的变化检测嵌入到4台计算机组成的集群上进行研究，使得变化检测的效率提高了 3.4 倍；叶琛(2014)将基于PCA的变化检测嵌入到由3台VPF2组成的集群上进行并行研究，实验表明变化检测的效率提高了3.34倍；胡维(2011)分别利用1—8个节点的计算机集群对基于主成分分析的变化检测进行研究，实验表明随着节点数的增多，加速比基本上成线性增加。用于并行的加速器主要包括：DSP、FPGA、GPU、Intel MIC、APU等，借助于这些硬件设备能够很大程度上提高遥感图像处理的效率。Kockara(2009)利用GPU的并行计算能力，在一定程度上实现了医学影像和遥感影像的实时变化检测；Zhu等人(2012)利用GT630M显卡将基于KL散度的变化检测效率提高了将近20倍。Asano等人(2009)将2维滤波、立体视觉以及K均值聚类算法分别在拥有4核的CPU、FPGA和GPU上进行研究分析，实验表明：对于简单的数值计算GPU表现出良好的并行计算性能，与FPGA相当；而对于复杂的逻辑性计算，CPU和FPGA运算性能优于GPU。可见，对于需要大规模的数值计算的图像处理而言，加速器体现出无法比拟的优势。

数字图像本质上就是一个2维数组，在大多数图像处理过程中，每个像素的运算过程是相互独立而且完全一样的数字计算，这恰恰符合了GPU高性能并行运算的条件，因此，GPU在遥感图像处理领域越来越受到人们的关注。然而直接对GPU操作却非常复杂，且灵活度很差。为了克服这些缺点，NVIDIA推出统一计算设备构架CUDA(Compute Unified Device Architecture)，CUDA是一种专门针对NVIDIA/GPU硬件而设计的并行运算架构，用于图像处理通用并行计算领域。CUDA的推出在学术界和产业界引起了热烈反响，尤其对于具有海量数据的遥感图像处理技术更是起到了巨大的推动作用。Yang等人(2008)将多种典型的遥感图像处理技术嵌入到GPU中研究，结果表明，图像直方图计算可以提高至少40倍，去云处理大约可以提高79倍，离散余弦变换可以提高8倍而边缘检测至少可以提高200倍；方留杨等人(2014)利用Tesla M2050 GPU将资源三号卫星影像的正射校正的效率提高了110倍；Zhang等人(2010)利用GPU分别将Sobel边缘检测以及形态学处理的运算效率提高了25倍和49倍，验证了GPU在图像处理中的可行性；肖汉和张祖勋(2010)设计了一种基于CUDA的快速影像匹配算法，相对于CPU上的串行算法提高了7倍；Hu和Yang(2013)把GPU运用到3D图像处理中，实验表明，GPU并行运算大大提高了3D图像的细化处理的效率。可见，基于GPU并行运算可以有效地提高遥感图像处理的效率。

目前基于GPU/CUDA并行技术已在遥感图像处理中得到了广泛应用，但是基于GPU/CUDA的变化检测并不多见(Zhu等，2012；Asano等，2009；Zhu和Cui等，2012)。为了可以在普通的PC机上能够实时地掌握土地利用和土地覆盖的情况，有必要把基于遥感图像的变化检测嵌入到GPU中进行进一步的研究。在众多的变化检测方法中，本文选用了基于变化矢量分析CVA(Change Vector Analysis)的遥感变化检测方法。该变化检测技术能够同时提取变化强度信息和变化方向信息。它将地物的变化表示为多光谱特征空间中像元的变化矢量，并通过图像处理和图像判读来确定地物变化的范围并进行统计分析，确定地物变化的类型信息。该方法简单、直观、易操作，大大减少了对变化区域做分类时作业人员的工作量，避免了人工操作带来的不便和误差；此外，该方法避免了分类过程引入虚假变化类型，有效的提高了检测精度。本文将CVA遥感变化检测过程嵌入到支持CUDA的NVIDIA GPU上进行研究，主要针对变化强度算子、映射表以及变化方向算子做分析，并对多组不同数据量的图像做实验。

CUDA是一种由NVIDIA推出的通用并行计算架构，通过该架构用户可以方便地对GPU进行操作。通过CUDA来操作GPU，避免了用户直接调用复杂的GPU接口函数，使用户专注于算法的研究。这种优势主要是由CUDA的软件结构决定的，该结构(如图 1所示)主要由以下3个层次组成：CUDA库函数、CUDA运行时API和CUDA驱动API。

如图 1所示，在基于CUDA的编程模型中，最重要的是CUDA C语言，它包含对C语言的最小扩展集和一个运行时库。在CUDA Library中有CUFFT、CUBLAS和CUDPP 3个函数库，它们提供了一些简单高效的常用函数，方便用户调用。而CUDA运行时API和CUDA驱动API提供了设备管理、上下文管理、存储管理、代码块管理、执行控制等应用程序接口，可以访问GPU的接口函数，实现对GPU的操作。

在实际应用过程中，用户只需要针对具体计算任务，确定程序中的并行处理部分；然后通过CUDA C在CPU中编写项目的并行程序。这些并行程序的源文件在运行过程中必须通过NVCC编译器进行编译。然而CUDA C语言编写的并行程序编译得到的只是GPU端代码，而如果需要管理GPU资源，将CUDA C编写的并行处理代码交给GPU来处理，就必须借助CUDA运行时API与CUDA驱动API来实现。需要注意的是，在一个程序中只能使用其中的一种，不能混合使用(张舒和褚艳利，2009)。GPU端的并行处理程序是以内核函数(Kernel)为单位的。然而Kernel函数并不是一个完整的程序，它只是整个CUDA程序中的一个可以被并行执行的步骤。在一个完整的CUDA程序中，除了设备端并行处理的Kernel函数以外，还有主机端的串行程序。在两个核之间运行的就是串行代码。理论上，CUDA核函数以网格(Grid)的形式执行，每个网格由若干个线程块(Block)组成，每一个线程块由若干个线程(Thread)组成。

Malila(1980)最早提出了变化矢量分析法。变化矢量分析法主要是从变化的强度和变化的方向两个角度对遥感影像进行变化检测。其中，变化强度图像用于表示研究区域的变化范围；变化方向图像用于表示研究区域具体的变化类型。传统的基于CVA的变化检测方法对两个不同时相的图像做差值运算，组成一幅新的差值多光谱图像。差值图像中的每一个位置处的像元灰度值都可以组成一个矢量，如式(1)所示。

$\Delta X = \left({\begin{array}{*{20}{c}}{\Delta {x_1}}\\{\Delta {x_2}}\\ \vdots \\{\Delta {x_k}}\\ \vdots \\{\Delta {x_{b and s}}}\end{array}} \right)= \left({\begin{array}{*{20}{c}}{x_{ij}^1\left({{t_2}} \right)- x_{ij}^1\left({{t_1}} \right)}\\{x_{ij}^2\left({{t_2}} \right)- x_{ij}^2\left({{t_1}} \right)}\\ \vdots \\{x_{ij}^k\left({{t_2}} \right)- x_{ij}^k\left({{t_1}} \right)}\\ \vdots \\{x_{ij}^{b and s}\left({{t_2}} \right)- x_{ij}^{b and s}\left({{t_1}} \right)}\end{array}} \right) $ (1)

式中，x_ij^k(t₁)和x_ij^k(t₂)分别表示两个时相的图像在(i，j)位置的像元灰度值，k表示图像的波段索引(其中k=1，2，…，b and s)，b and s是多光谱图像的波段数。下面以2维空间中矢量的大小和方向来说明变化矢量分析法的原理，如图 3所示。图 3中矢量的方向代表了地物的变化类型，矢量的长度代表了地物变化大小程度。

基于CVA的变化检测的强度图像是利用欧氏距离得到的，即变化强度图像的灰度值是由差值多光谱图像上对应位置的像元灰度值的平方和再开方得到的，表达式如式(2)：

$ DX = \sqrt {{{\left({\Delta {x_1}} \right)}^2} + {{\left({\Delta {x_2}} \right)}^2} + \cdots + {{\left({\Delta {x_{b and s}}} \right)}^2}} $ (2)

式中，DX表示某一像元的变化强度。在实际应用中，可根据区域具体情况对变化强度设定一个阈值。若像元的变化强度在此阈值范围内，可以认为该点未发生变化；若超出此阈值范围，则认为发生了变化。通过合理的阈值分割处理，就可以有效地区分变化和未变化区域。变化矢量分析法通常是选用简单差值法(如公式1所示)来计算变化矢量以及变化模式，而变化方向图像是根据变化矢量ΔX的方向进行编码得到的。对于差值图像而言，每个像元都有b and s个不同的值，每个波段上的数值有可能为负值，有可能为零，也可能为正值。这3种可能分别对应3种编码：-1，0，1。若以3个波段为例，如表 1所示。

每一位置的像元对应一个向量，向量的每一个元素是-1，0，1中的一个值。对于3个波段的图像，假设第一波段减弱，则编码为1；第二波段增强，则编码为-1；第三波段不变，则编码为0。这3个编码组成一个变化方向向量，即ΔX_ij=[1，-1，0]，其中(i，j)为当前像元的位置。然后遍历所有情况可以得到3^{b and s}种不同的向量，对应着3^{b and s}种不同的变化。最后建立这些向量与1—3^{b and s}范围上的一一映射关系，这样就实现了不同的方向矢量对应不同的灰度值。对差值图像上的每个像素进行判断，根据这种映射关系就可以得到变化方向图像。3个波段的图像对应的映射表如表 2所示。

对于上述得到的变化方向图像而言，各种类型的变化都是在变化强度确定的变化范围内部的，变化强度确定的未变化区域没有发生任何类型的变化。此外，变化方向图像的量化级为3^{b and s}。理论上，如果每一种变化类型都发生，图像上应该有3^{b and s}种不同的灰度值；但是对于某一特定的研究区域，并不是每一种变化类型都会发生。事实上，大多数的变化类型并没又发生或者变化区域很小，即变化方向图像上不同灰度值的种类远远小于其量化级。通过重新量化，变化方向图像的量化级就会大大较少，并使得每一种变化类型都有对应的变化区域。然后经过合理的伪彩色增强技术就可以把变化区域内不同类型的变化区分开来。

本文的实验主机平台是 Intel(R)Core(TM)i5-3470 CPU，主频率为3.2 GHz，系统内存为4 G，CPU有4个核；显卡采用的NVIDIA GeForce GT 610，GPU型号是GF119，属于Fermi构架。GPU的核心频率为1.62 GHz，显卡内存为1 G，GPU运算核为48个。操作系统为Microsoft Window 8.1 64位，开发环境为Microsoft Visual Studio 2010，整个实验基于GDAL- 1.9.2 、GeForce 340.52 Driver和CUDA6.0 包括：CUDA Toolkit 6.0和CUDA Samples 6.0。

为了实现对土地利用和覆盖变化的实时监测、对海量遥感数据的快速处理，本文结合CUDA编程模型，针对基于CVA的变化检测算法设计了一种并行处理模型。该模型将CVA变化检测中的变化强度检测、映射表构建和变化方向检测3个过程嵌入到GPU中进行研究。本文利用CPU串行、OpenMP并行以及GPU/CUDA并行3种方法对不同数据量遥感进行CVA变化检测实验。结果表明，与另外两种方法相比，基于CUDA C的CVA变化检测效率有明显提高。当图像的数据量相对较小时，结合CUDA的CVA变化检测耗时等于在CPU上对一个像元计算的耗时，一般不会随着数据量的增加而增加，从而大大地提高了变化检测的运算速度；当图像数据量比较大时，就需要对图像进行分块并行处理，运算耗时与分块的个数成正比例关系。根据图 12中GPU耗时随着图像分块大小的变化，可知对于大数据量的遥感影像变化检测而言，在内存和显存不溢出的前提下，增加分块的尺寸，可以减少分块的数量，从而提高变化检测的效率。

总之，在实际处理过程中，针对不同的数据，需要充分地运用全局存储器和各种具有高缓存的存储器来提高数据的调用和传输效率，减少Kernel的个数和图像分块的个数以减少数据在不同设备间传输的次数，并调节Grid和Block的大小以便为并行任务合理地分配处理单元，从而充分利用GPU的并行运算能力，提高遥感图像处理的效率，以满足各种领域对数字图像的处理要求。

附录A

在利用GPU设备进行并行变化检测之前，首先需要获取GPU设备的基本信息，包括GPU个数、每个GPU设备的名称等信息、Block的最大尺寸以及Grid的最大尺寸等信息。GPU设备初始化伪代码如下：

_host_ bool InitCUDA()

Dim count, gpu_id As Integar

    cudaGetDeviceCount(&count);

    For gpu_id =0 To count

        cudaDeviceProp prop;

        cudaGetDeviceProperties(&prop, gpu_id);

//Block的最大尺寸

MaxBlockX=prop.maxThreadsDim[0];

MaxBlockY=prop.maxThreadsDim[1];

MaxBlockZ=prop.maxThreadsDim[2];

//Grid的最大尺寸

MaxGridX=prop.maxGridSize[0];

MaxGridY=prop.maxGridSize[1];

MaxGridZ=prop.maxGridSize[2];

……………………

        EndFor

End InitCUDA

附录B

在进行CVA变化检测之前，我们需要构造映射表，用于对变化矢量进行编码，用于计算变化方向图像。构造映射表的伪代码如下所示：

_global_ void MappingTable_Kernel

(Input bands; Output DirectionVector, DirectionGrayvalue )

Dim index, k, index_k As Integer

index = threadIdx.x+blockIdx.x*blockDim.x;

k = threadIdx.y+blockIdx.y*blockDim.y;

index_k=pow(3.0,(bands-k));

①//反射率减少

If ((index+1)%index_k≤(index_k/3)

AND

(index+1)%index_k !=0)

Then DirectionVector[index][k]=-1

②//反射率不变

Else

If((index_k/3≤(index+1)%index_k≤2*index_k/3))

Then DirectionVector[index][k]=0

③//反射率增加

Else

If (index+1)%index_k>2×index_k/3 or

(index+1)%index_k=0

Then DirectionVector[index][k]=1

Eendif

//与DirectionVector的行一一对应的编码值

DirectionGrayvalue[index]=index+1

End MappingTable_Kernel

其中，DirectionVector表示所有变化矢量组成的矩阵，空间大小为[3^bands，bands]；DirectionGrayvalue表示与变化矢量一一对应的灰度值，空间大小为[3^bands，1]，index对应映射表的行，k对应映射表的列。

附录C

在基于CVA变化检测的程序设计中，最重要的就是结合像元与线程的对应关系，设计变化强度算子和变化方向算子，其对应的Kernel函数的如下所示：

_global_ void ChangeDetection Kerenl (

Input Inputimg1, Inputimg2,

DirectionVector, DirectionGrayvalue, height, width, bands;

Output Magnitudeimg, Directionimg )

Dim Row, Col, offset As Integer

Row = threadIdx.y+blockIdx.y×blockDim.y;

Col = threadIdx.x+blockIdx.x×blockDim.x;

offset = Col+Row×width;

GrayScale_Direction=pow(3.0, bands);

①//变化强度算子

For k= 0 To bands

Magnitudeimg[offset]+=pow((Inputimg2[k][offset]-Inputimg1[k][offset]),2);

EndFor

    Magnitudeimg[offset]=sqrt(Magnitudeimg[offset]);

②//变化方向算子

For k= 0 To bands

     Differimg[k][offset]=Inputimg2[k][offset]-Inputimg1[k][offset];

Dim threshold [bands] As Integer

//反射率减少

If Differimg ≤ -threshold(k)

    Then Differimg[k][offset]=-1

     //反射率不变

     Else

If -threshold(k) ≤ Differimg[k][offset] ≤

threshold(k)

Then Differimg[k][offset]=0

    //反射率增强

Else

If Differimg[k][offset] ≥ threshold(k)

Then Differimg[k][offset]=1;

Endif

//根据当前像素的变化矢量查询对应的编码值

For index=0 To GrayScale_Direction

Dim index_XY As Integer

index_XY=0;

//判断DirectionVector的第index行的每个值与当前像素的矢量的每个值是否相等

For k=0 To bands

If DirectionVector[index][k]=Differimg[k][offset]

Then index_XY++;

Endif

EndFor

//如果当前像素的变化矢量与DirectionVector的第index行相等，则将对应的DirectionGrayvalue的index值赋给变化方向图像的对应当前像素

If index_XY=bands

Then Directionimg[offset]=DirectionGrayvalue[index]

Endif

    EndFor

EndFor

End ChangeDetection Kerenl

其中，Magnitudeimg和Directionimg分别表示变化强度图像指针和变化方向图像指针，它们的空间大小都是[1，height×width]；Differimg表示差值图像指针；Inputimg1和Inputimg2分别表示时相1和时相2图像的指针，空间大小为[bands，height×width]，bands为图像的波段数，height表示图像的行数，width表示图像的列数；GrayScale_Direction表示变化方向图像中最大的理论值，即3^bands；threshold表示各个波段的阈值，threshold[k]为第k个波段的阈值，用于判断每个波段的反射率是否减少、不变或增加。