CN103942788B - 高光谱遥感图像特征提取方法及装置 - Google Patents

Abstract

本发明实施例公开了一种高光谱遥感图像特征提取方法及装置，通过CPU和GPU协同运算实现高光谱遥感图像的快速特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

Description

高光谱遥感图像特征提取方法及装置

技术领域

本发明涉及图像处理技术领域，更具体地说，涉及一种高光谱遥感图像特征提取方法及装置。

背景技术

高光谱遥感图像所具有的大量波段为地物信息提取提供了极其丰富的光谱、辐射和几何信息，有助于更加精细的地物分类和目标识别，然而，光谱波段的增多也会导致信息的冗余和数据处理复杂性的提高，因此，特征提取在高光谱遥感图像处理中具有十分重要的作用。

发明人在实现本发明的过程中发现，现有技术通过CPU提取高光谱遥感图像特征，其特征提取速率较低。

发明内容

本发明的目的是提供一种高光谱遥感图像特征提取方法及装置，以提高高光谱遥感图像特征提取速率。

为实现上述目的，本发明提供了如下技术方案：

一种高光谱遥感图像特征提取方法，应用于中央处理器CPU，包括：

将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

将所述第一矩阵的特征向量矩阵加载至所述GPU中；

向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

将所述第二矩阵读入到所述CPU中。

上述方法，优选的，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三维的长度为所述高光谱遥感图像的波段数；

所述每一个图像块通过一个线程进行噪声估计包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

上述方法，优选的，所述每一个图像块通过一个线程进行噪声估计包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

利用各个图像块的回归参数对各个图像块进行参数估计；

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

一种高光谱遥感图像特征提取方法，应用于图形处理器GPU，包括：

接收中央处理器CPU加载的待处理的高光谱遥感图像；

在接收到所述CPU发送的第一指令后，执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量矩阵由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的获取方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

上述方法，优选的，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三位的长度为所述高光谱遥感图像的波段数；

所述每一个图像块通过一个线程进行噪声估计包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

上述方法，优选的，所述每一个图像块通过一个线程进行噪声估计包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

利用各个图像块的回归参数对各个图像块进行参数估计；

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

一种高光谱遥感图像特征提取装置，包括：

第一加载模块，用于将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

第一指令发送模块，用于向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二指令发送模块，用于向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第一读取模块，用于将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

第一计算模块，用于计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

第二加载模块，用于将所述第一矩阵的特征向量矩阵加载至所述GPU中；

第三指令发送模块，用于向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

第二读取模块，用于将所述第二矩阵读入到所述CPU中。

一种高光谱遥感图像特征提取装置，包括：

接收模块，用于接收中央处理器CPU加载的待处理的高光谱遥感图像；

噪声估计模块，用于在接收到所述CPU发送的第一指令后，执行空噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二计算模块，用于在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第三计算模块，用于在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的计算方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

上述装置，优选的，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三位的长度为所述高光谱遥感图像的波段数；所述噪声估计模块包括：

划分单元，用于将所述高光谱遥感图像的各个波段均划分为若干个图像块；

对应关系确定单元，用于依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

通过以上方案可知，本申请提供的一种高光谱遥感图像特征提取方法，通过CPU和GPU协同运算实现高光谱遥感图像的快速特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本申请实施例提供的高光谱遥感图像特征提取方法的一种实现流程图；

图2为本申请实施例提供的线程、二维线程块、和三维线程网格、之间的关系示意图；

图3为本申请实施例提供的高光谱遥感图像特征提取方法的另一种实现流程图；

图4为本申请实施例提供的高光谱遥感图像特征提取方法的又一种实现流程图；

图5为本申请实施例提供的高光谱遥感图像特征提取装置的一种结构示意图；

图6为本申请实施例提供的高光谱遥感图像特征提取装置的另一种结构示意图；

图7为本申请实施例提供的噪声估计模块的一种结构示意图；

图8为本申请实施例提供的方案与现有技术在不同遥感图像数据量情况下的遥感图像特征提取时间的对比图。

说明书和权利要求书及上述附图中的术语“第一”、“第二”、“第三”“第四”等（如果存在）是用于区别类似的部分，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本申请的实施例能够以除了在这里图示的以外的顺序实施。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本申请实施例提供的高光谱遥感图像特征提取方法的一种实现流程图如图1所示，本申请实施例提供的高光谱遥感图像特征提取方法应用于中央处理器CPU（CentralProcessing Unit），具体可以包括：

步骤S11：将待处理的高光谱遥感图像加载至图形处理器GPU（GraphicsProcessing Unit）中，所述GPU能够同时运行若干个线程；

本申请实施例中，CPU对GPU进行初始化，然后将待提取特征的高光谱遥感图像从CPU中拷贝到图形处理器GPU中。

步骤S12：向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

所述第一指令可以是调用噪声估计方法的核函数的指令，也就是说，CPU可以调用GPU中的噪声估计方法的核函数。GPU接收到CPU发送的调用噪声估计方法的核函数的指令后，执行所述噪声估计方法的核函数，具体的，可以将高光谱遥感图像划分为若干个预设大小的图像块，例如，可以将高光谱遥感图像划分为10*10大小的图像块，当然，也可以划分为其它大小的图像块，这里不做具体限定。

GPU中每个线程只处理一个图像块，而由于GPU能够同时运行若干个线程，因此，在GPU端，噪声估计方法可以通过多个线程并行执行。也就是说，本申请实施例中，将需要在CPU端串行执行的大运算量的噪声估计运算转移到GPU中，通过并行运算方式实现，从而减少了噪声估计所需要的时间。

CPU向GPU发送第一指令后，进入等待状态直到GPU中噪声估计方法的核函数执行完毕。

步骤S13：向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

当GPU执行完噪声估计的核函数后，CPU向GPU发送第二指令。GPU在接收到第二指令后，计算高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；其中，

噪声数据的协方差矩阵的计算方法可以包括：

计算原始噪声数据（为便于叙述，本申请实施例中，将步骤S12得到的噪声数据称为原始噪声数据）的各个波段的均值，包括：计算每个波段的原始噪声数据的总和，将总和除以每个波段的象元的个数得到原始噪声数据的各个波段的第一均值。

将每个波段的原始噪声数据减去对应波段的第一均值，得到去中心化后的噪声数据；

将去中心化的各个波段的噪声数据进行向量化，得到去中心化后的噪声数据的各个波段的向量数据；

将去中心化后的噪声数据的各个波段的向量数据构成第四矩阵；

将第四矩阵转置后乘以第四矩阵，得到的结果再除以波段数减一的值就得到噪声数据的协方差矩阵。

原始数据（即CPU拷贝到GPU中的高光谱遥感图像数据）的协方差矩阵的计算方法可以包括：

计算原始数据的各个波段的均值，包括：计算每个波段的原始数据的总和，将总和除以每个波段的象元的个数得到的原始数据的各个波段的第二均值；

将每个波段的原始数据减去对应波段的第二均值，得到各个波段的去中心化后的数据；

将去中心化后的各个波段的数据进行向量化，得到去中心化后的数据的各个波段的向量数据；

将去中心化的数据的各个波段的向量数据构成第五矩阵；

将第五矩阵转置后乘以第五矩阵，得到的结果再除以波段数减一的值就得到原始数据的协方差矩阵。

其中，所述第二指令可以为一条指令，即指示GPU在计算完噪声数据的协方差矩阵后，计算原始数据的协方差矩阵；

所述第二指令也可以包括两条子指令，向GPU发送第二指令时，CPU可以先向GPU发送第一子指令，指示GPU计算噪声数据的协方差矩阵，当GPU计算完噪声数据的协方差矩阵后，CPU再向GPU发送第二子指令，指示GPU计算原始数据的协方差矩阵。

步骤S14：将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

本申请实施例中，CPU可以在GPU计算完噪声数据的协方差矩阵，和原始数据的协方差矩阵后，将噪声数据的协方差矩阵和原始数据的协方差矩阵一起拷贝到CPU中；

CPU可以在GPU计算完噪声数据的协方差矩阵后，先将噪声数据的协方差矩阵拷贝到CPU中，当GPU计算完原始数据的协方差矩阵后，再将原始数据的协方差矩阵拷贝到CPU中。

步骤S15：计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

本步骤在CPU端执行，其中，特征向量的矩阵中的特征向量可以按从小到大的顺序排列，也可以按照从大到小的顺序排列。本申请实施例中，优先选择按照从小到大的顺序排列。

步骤S16：将所述第一矩阵的特征向量矩阵加载至所述GPU中；

CPU将计算得到的特征向量矩阵拷贝到GPU中。

步骤S17：向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

本申请实施例中，在GPU端计算特征向量矩阵与去中心化后的高光谱遥感图像数据的乘积，得到的第二矩阵即是所要提取的高光谱遥感图像的特征。

步骤S18：将所述第二矩阵读入到所述CPU中。

GPU计算得到高光谱遥感图像的特征后，CPU将特征数据读入到CPU中。

CPU在将所提取的特征数据读入到CPU中后，可以将特征数据保存为磁盘图像文件。

本申请提供的一种高光谱遥感图像特征提取方法，通过CPU和GPU协同运算实现高光谱遥感图像的快速特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

另外，本申请实施例中，先对高光谱遥感图像进行噪声估计，再计算原始高光谱遥感图像数据的协方差矩阵，可以只向GPU中加载一次高光谱遥感图像即可，从而避免先计算原始高光谱遥感图像数据的协方差矩阵，再对高光谱遥感图像进行噪声估计导致需要向GPU中加载两次高光谱遥感图像的问题。

上述实施例中，优选的，为了进一步提高高光谱遥感图像的特征提取速率，CPU可以调用GPU中的CUDA（Compute Unified Device Architecture，统一计算设备架构）这一编程模型中的CUBLAS库（基本线性代数计算库）计算噪声数据的协方差矩阵和原始数据的协方差矩阵，以及第二矩阵。

上述实施例中，优选的，为了方便线程快速确定其所处理的图像块，可以定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三维的长度为所述高光谱遥感图像的波段数；

如图2所示，图2为本申请实施例提供的线程（图中所示Thread）、二维线程块（图中所示Block）和三维线程网格（图中所示Grid）之间的关系示意图；

本申请实施例中，对于GPU中可以执行的若干个线程，预先定义二维线程块，其中，每一个二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引，也就是说，本申请实施例中，将GPU中可以执行的若干个线程划分为多个二维线程块，一个二维线程块为一个线程组。例如，二维线程块的第一维和第二维的长度可以均为16，即，一个二维线程块中包括16*16=256个线程。

还预先定义三维线程网格，该三维线程网格中包括若干个二维线程块，在该三维线程网格中，每一个二维线程块对应一个三维数组索引；

所述每一个图像块通过一个线程进行噪声估计可以包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

当一个线程需要确定其所要处理的图像块是高光谱遥感图像中的哪一个图像块时，可以先通过该线程所属的二维线程块的第三维索引（即该线程所属的二维线程块所对应的波段）确定该线程所对应的波段。

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。表示取整运算。

也就是说，本申请实施例中，每一个线程实体都可以通过（blockIdx.x，blockIdx.y，blockIdx.z，threadIdx.x，threadIdx.y）五个全局参数得到标识自身的一个唯一编号，通过这个编号，以及第三公式，线程可以确定自己要处理的位置是图像中哪一个图像块。

其中，blockIdx.z为所述第一线程所属的二维线程块在所述三维线程网格中的第三维索引，即所述第一线程所属的二维线程块所对应的波段的索引，也就是所述第一线程所对应的波段的索引；

需要说明的是，为方便计算，本申请实施例中，各个线程、各个二维线程块的每一维索引的最小值均为零，最大值为相应维度的长度减一后的值。

上述实施例中，优选的，所述每一个图像块通过一个线程进行噪声估计可以包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，可以包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

本方案中，利用多元线性回归方法计算回归参数，具体如何利用多元线性回归方法计算回归参数属于本领域的公知常识，这里不再赘述。

利用各个图像块的回归参数对各个图像块进行参数估计；

本申请实施例中，利用回归参数对各个数据进行参数估计，得到重新估计的图像块，具体如何利用回归参数对各个数据进行参数估计属于本领域的公知常识，这里不再赘述。

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

本申请实施例提供的高光谱遥感图像特征提取方法的另一种实现流程图如图3所示，本申请实施例提供的高光谱遥感图像特征提取方法应用于GPU，具体可以包括：

步骤S31：接收中央处理器CPU加载的待处理的高光谱遥感图像；

本申请实施例中，CPU对GPU进行初始化，然后将待处理的高光谱遥感图像拷贝到GPU中。

步骤S32：在接收到所述CPU发送的第一指令后，执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

所述第一指令可以是调用噪声估计方法的核函数的指令，也就是说，CPU可以调用GPU中的噪声估计方法的核函数。GPU接收到CPU发送的调用噪声估计方法的核函数的指令后，执行所述噪声估计方法的核函数，具体的，可以将高光谱遥感图像划分为若干个预设大小的图像块，例如，可以将高光谱遥感图像划分为10*10大小的图像块，当然，也可以划分为其它大小的图像块，这里不做具体限定。GPU中每个线程只处理一个图像块，而由于GPU能够同时运行若干个线程，因此，在GPU端，噪声估计方法可以通过多个线程并行执行。也就是说，本申请实施例中，将需要在CPU端串行执行的大运算量的噪声估计运算转移到GPU中，通过并行运算方式实现，从而减少了噪声估计所需要的时间。

CPU向GPU发送第一指令后，进入等待状态直到GPU中噪声估计方法的核函数执行完毕。

步骤S33：在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

当GPU执行完噪声估计的核函数后，CPU向GPU发送第二指令。GPU在接收到第二指令后，计算高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；其中，

噪声数据的协方差矩阵的计算方法可以包括：

计算原始噪声数据（为便于叙述，本申请实施例中，将步骤S22得到的噪声数据称为原始噪声数据）的各个波段的均值，包括：计算每个波段的原始噪声数据的总和，将总和除以每个波段的象元的个数得到原始噪声数据的各个波段的第一均值。

将每个波段的原始噪声数据减去对应波段的第一均值，得到去中心化后的噪声数据；

将去中心化的各个波段的噪声数据进行向量化，得到去中心化后的噪声数据的各个波段的向量数据；

将去中心化后的噪声数据的各个波段的向量数据构成第四矩阵；

将第四矩阵转置后乘以第四矩阵，得到的结果再除以波段数减一的值就得到噪声数据的协方差矩阵。

原始数据（即CPU拷贝到GPU中的高光谱遥感图像的数据）的协方差矩阵的计算方法可以包括：

计算原始数据的各个波段的均值，包括：计算每个波段的原始数据的总和，将总和除以每个波段的象元的个数得到的原始数据的各个波段的第二均值；

将每个波段的原始数据减去对应波段的第二均值，得到各个波段的去中心化后的数据；

将去中心化后的各个波段的数据进行向量化，得到去中心化后的数据的各个波段的向量数据；

将去中心化的数据的各个波段的向量数据构成第五矩阵；

将第五矩阵转置后乘以第五矩阵，得到的结果再除以波段数减一的值就得到原始数据的协方差矩阵。

其中，所述第二指令可以为一条指令，即GPU依据该指令在计算完噪声数据的协方差矩阵后，计算原始数据的协方差矩阵；

所述第二指令也可以包括两条子指令，向GPU发送第二指令时，CPU可以先向GPU发送第一子指令，GPU在接收到第一指令后，计算噪声数据的协方差矩阵，当GPU计算完噪声数据的协方差矩阵后，CPU再向GPU发送第二子指令，GPU在接收到第二子指令后才计算原始数据的协方差矩阵。

步骤S34：在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量矩阵由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的获取方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

本申请实施例中，CPU可以在GPU计算完噪声数据的协方差矩阵，和原始数据的协方差矩阵后，将噪声数据的协方差矩阵和原始数据的协方差矩阵一起拷贝到CPU中；

CPU可以在GPU计算完噪声数据的协方差矩阵后，先将噪声数据的协方差矩阵拷贝到CPU中，当GPU计算完原始数据的协方差矩阵后，再将原始数据的协方差矩阵拷贝到CPU中。

其中，特征向量的矩阵中的特征向量可以按从小到大的顺序排列，也可以按照从大到小的顺序排列。本申请实施例中，优先选择按照从小到大的顺序排列。

本申请实施例中，在GPU端计算特征向量矩阵与去中心化后的高光谱遥感图像数据的乘积，得到的第二矩阵即是所要提取的高光谱遥感图像的特征数据。

GPU计算得到特征数据后，CPU可以将特征数据读入到CPU中。

本申请提供的一种高光谱遥感图像特征提取方法，通过CPU和GPU协同运算实现高光谱遥感图像的快速特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

上述实施例中，优选的，为了方便线程快速确定其所处理的图像块，可以定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三维的长度为所述高光谱遥感图像的波段数；

所述每一个图像块通过一个线程进行噪声估计可以包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

上述实施例中，优选的，所述每一个图像块通过一个线程进行噪声估计可以包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，可以包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

本方案中，利用多元线性回归方法计算回归参数，具体如何利用多元线性回归方法计算回归参数属于本领域的公知常识，这里不再赘述。

利用各个图像块的回归参数对各个图像块进行参数估计；

本申请实施例中，利用回归参数对各个数据进行参数估计，得到重新估计的图像块，具体如何利用回归参数对各个数据进行参数估计属于本领域的公知常识，这里不再赘述。

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

下面结合CPU和GPU对本申请实施例提供的高光谱遥感图像的特征提取方法的实现流程进行说明，本申请实施例提供的又一种高光谱遥感图像的特征提取方法的实现流程如图4所示，可以包括：

步骤S41：CPU将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

步骤S42:CPU向GPU发送第一指令；

步骤S43：GPU在接收到所述CPU发送的第一指令后，执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

步骤S44：CPU在GPU完成噪声估计后，向GPU发送第二指令；

CPU可以主动监测GPU是否完成噪声估计，也可以由GPU在完成噪声估计后通知GPU。

步骤S45：GPU在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

步骤S46：CPU将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

步骤S47：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

步骤S48：CPU将所述第一矩阵的特征向量矩阵加载至所述GPU中；

步骤S49：CPU向GPU发送第三指令；

步骤S410：GPU在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

步骤S411：CPU将所述第二矩阵读入到CPU中。

与方法实施例相对应，本申请实施例还提供一种高光谱遥感图像特征提取装置，应用于CPU，本申请实施例提供的高光谱遥感图像特征提取装置的一种结构示意图如图5所示，可以包括：

第一加载模块51，第一指令发送模块52，第二指令发送模块53，第一读取模块54，第一计算模块55，第二加载模块56，第三指令发送模块57和第二读取模块58；其中，

第一加载模块51用于将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

第一指令发送模块52用于向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二指令发送模块53用于在所述GPU噪声估计完成时，向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第一读取模块54用于在所述GPU计算得到噪声数据的协方差矩阵和原始数据的协方差矩阵时，将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

第一计算模块55用于计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

第二加载模块56用于将所述第一矩阵的特征向量矩阵加载至所述GPU中；

第三指令发送模块57用于向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

第二读取模块58用于在所述GPU计算得到第二矩阵时，将所述第二矩阵读入到所述CPU中。

本申请实施例提供的一种高光谱遥感图像特征提取装置，通过CPU与GPU协同运算实现高光谱遥感图像的快速特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

本申请实施例提供的高光谱遥感图像特征提取装置应用于GPU，本申请实施例提供的高光谱遥感图像特征提取装置的另一种结构示意图如图6所示，可以包括：

接收模块61，噪声估计模块62，第二计算模块63和第三计算模块64；其中，

接收模块61用于接收中央处理器CPU加载的待处理的高光谱遥感图像；

噪声估计模块62用于在接收到所述CPU发送的第一指令后，执行空噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二计算模块63用于在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第三计算模块64用于在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的计算方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

本申请实施例提供的一种高光谱遥感图像特征提取装置，通过CPU与GPU协同运算实现高光谱遥感图像的特征提取，其中，在CPU上实现算法的逻辑控制，以及读取、存储数据等功能，而把大规模的运算（如噪声估计以及协方差矩阵的计算）功能在GPU中实现，其中，在进行噪声估计时，将高光谱遥感图像划分为若干个图像块，通过多线程并行执行，从而提高了高光谱遥感图像特征提取速率。

上述实施例中，优选的，为了方便线程快速确定其所处理的图像块，可以定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三维的长度为所述高光谱遥感图像的波段数；

所述噪声估计模块的一种结构示意图如图7所示，可以包括：

划分单元71，用于将所述高光谱遥感图像的各个波段均划分为若干个图像块；

对应关系确定单元72，用于依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；

其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

为了更好的说明本申请的有益效果，下面给出本申请实施例提供的方案与现有技术在不同高光谱遥感图像数据量情况下的特征提取时间的对比图，如图8所示；

图8中，表示现有技术中通过CPU串行运算实现高光谱遥感图像特征提取所需要的时间；“——”表示通过本申请实施例的并行运算实现高光谱遥感图像特征提取所需要的时间。

由图8可以看出，随着高光谱遥感图像数据量的增加，现有技术需要的计算时间急剧加剧，而本申请实施例的并行运算时间增加不明显，本申请实施例提供的高光谱遥感图像特征提取的加速效果显著。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。

Claims (9)

1.一种高光谱遥感图像特征提取方法，应用于中央处理器CPU，其特征在于，包括：

将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

将所述第一矩阵的特征向量矩阵加载至所述GPU中；

向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

将所述第二矩阵读入到所述CPU中。

2.根据权利要求1所述的方法，其特征在于，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三维的长度为所述高光谱遥感图像的波段数；

所述每一个图像块通过一个线程进行噪声估计包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

3.根据权利要求1或2所述的方法，其特征在于，所述每一个图像块通过一个线程进行噪声估计包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

利用各个图像块的回归参数对各个图像块进行参数估计；

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

4.一种高光谱遥感图像特征提取方法，应用于图形处理器GPU，其特征在于，包括：

接收中央处理器CPU加载的待处理的高光谱遥感图像；

在接收到所述CPU发送的第一指令后，执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量矩阵由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的获取方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

5.根据权利要求4所述的方法，其特征在于，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三位的长度为所述高光谱遥感图像的波段数；

所述每一个图像块通过一个线程进行噪声估计包括：

依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；

依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。

6.根据权利要求4或5所述的方法，其特征在于，所述每一个图像块通过一个线程进行噪声估计包括：

每个线程通过空间-光谱维去相关噪声估计方法对图像块进行噪声估计，包括：

获取回归参数，具体包括：如果所要处理的图像块位于所述高光谱遥感图像的第一个波段，则利用多元线性回归方法计算当前波段图像块的数据与第二个波段的相应位置的图像块的数据的回归参数；如果所要处理的图像块位于所述高光谱遥感图像的最后一个波段，则利用多元线性回归方法计算当前波段图像块与前一个波段的相应位置的图像块的数据的回归参数；否则，利用多元线性回归方法计算当前波段图像块的数据与前一个波段和后一个波段的相应位置的图像块的数据的回归参数；

利用各个图像块的回归参数对各个图像块进行参数估计；

利用原始图像块的数据减去估计得到的图像块的数据得到各个图像块的噪声数据。

7.一种高光谱遥感图像特征提取装置，应用于中央处理器CPU，其特征在于，包括：

第一加载模块，用于将待处理的高光谱遥感图像加载至图形处理器GPU中，所述GPU能够同时运行若干个线程；

第一指令发送模块，用于向所述GPU发送第一指令，所述第一指令用于指示所述GPU执行噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二指令发送模块，用于向所述GPU发送第二指令，所述第二指令用于指示所述GPU计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第一读取模块，用于将所述GPU计算得到的噪声数据的协方差矩阵和原始数据的协方差矩阵读入到所述CPU中；

第一计算模块，用于计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列；

第二加载模块，用于将所述第一矩阵的特征向量矩阵加载至所述GPU中；

第三指令发送模块，用于向所述GPU发送第三指令，所述第三指令用于指示所述GPU计算所述特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

第二读取模块，用于将所述第二矩阵读入到所述CPU中。

8.一种高光谱遥感图像特征提取装置，应用于图形处理器GPU，其特征在于，包括：

接收模块，用于接收中央处理器CPU加载的待处理的高光谱遥感图像；

噪声估计模块，用于在接收到所述CPU发送的第一指令后，执行空噪声估计方法的核函数，以对所述高光谱遥感图像进行噪声估计，包括：将所述高光谱遥感图像的各个波段均划分为若干个图像块，其中，每一个图像块通过一个线程进行噪声估计，从而获取所述高光谱遥感图像的各个波段的噪声数据；

第二计算模块，用于在接收到所述CPU发送的第二指令后，计算所述高光谱遥感图像的噪声数据的协方差矩阵，以及所述高光谱遥感图像的原始数据的各个波段的均值和所述原始数据的协方差矩阵；

第三计算模块，用于在接收到所述CPU发送的第三指令后，计算特征向量矩阵与去中心化后的数据的乘积，得到第二矩阵，其中，所述去中心化后的数据的获取方法为：将每一个波段的原始数据减去原始数据的相应波段的均值得到每一个波段的第三矩阵，将每一个第三矩阵向量化得到每一个波段的向量数据，所述去中心化后的数据由各个波段的向量数据构成；

其中，所述特征向量由所述CPU计算得到，并在发送所述第三指令前加载至所述GPU中，所述特征向量矩阵的计算方法包括：CPU计算所述原始数据的协方差矩阵的逆矩阵与所述噪声数据的协方差矩阵的乘积，得到第一矩阵；计算所述第一矩阵的特征向量矩阵，所述特征向量矩阵中的特征向量按预设的大小顺序排列。

9.根据权利要求8所述的装置，其特征在于，定义二维线程块，所述二维线程块中包括至少一个线程，所述二维线程块中，每一个线程对应一个二维数组索引；定义三维线程网格，所述三维线程网格中包括若干个二维线程块，所述三维线程网格中，每一个二维线程块对应一个三维数组索引；所述三维线程网格的第一维的长度依据第一公式确定，所述三维线程网格的第二维的长度依据第二公式确定；所述三维线程网格的第三位的长度为所述高光谱遥感图像的波段数；所述噪声估计模块包括：

划分单元，用于将所述高光谱遥感图像的各个波段均划分为若干个图像块；

对应关系确定单元，用于依据第一线程所属的二维线程块的第三维索引确定所述第一线程所对应的高光谱遥感图像的波段；依据第三公式在所确定波段的高光谱遥感图像中确定所述第一线程所对应的图像块；其中，

所述第一公式为：

所述第二公式为：

所述第三公式为： $\left\{\begin{array}{c}{x}_{\mathrm{Index}}=\mathrm{blockIdx}.x*\mathrm{blockdim}.x+\mathrm{threadIdx}.x\\ y_{\mathrm{Index}}=\mathrm{blockIdx}.y*\mathrm{blockdim}.y+\mathrm{threadIdx}.y\end{array}\right.$

其中，gridDim.x为三维线程网格的第一维的长度；gridDim.y为三维线程网格的第二维的长度；w为所述高光谱遥感图像的宽度；h为所述高光谱遥感图像的高度；a为所述图像块的宽度；b为所述图像块的高度；blockDim.x为二维线程块的第一维的长度；blockDim.y为二维线程块的第二维的长度；x_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第一维索引；y_Index为所需确定的图像块在所确定波段的高光谱遥感图像中的第二维索引；blockIdx.x为所述第一线程所属的二维线程块在所述三维线程网格中的第一维索引；blockIdx.y为所述第一线程所属的二维线程块在所述三维线程网格中的第二维索引；threadIdx.x为所述第一线程在其所属的二维线程块中的第一维索引；threadIdx.y为所述第一线程在其所属的二维线程块中的第二维索引。