CN107644393A - 一种基于gpu的丰度估计算法的并行实现方法 - Google Patents

Abstract

本发明公开了一种基于GPU的丰度估计算法的并行实现方法，包括：1、将高光谱混合像元数据H存储至GPU端；2、获取端元向量m_i，并将所获取的端元向量m_i逐一分配至CPU端的p个并行线程上；3、获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,...,m_i‑1]的向量4、获得端元集M所对应的正交基并计算对应的内积值；5、计算像元的光谱向量r在方向上的投影；6、计算各端元m_p对应的丰度结果α_p，并存储至GPU端。本发明通过GPU并行处理的核函数计算丰度值，且该并行机制比在传统CPU串行机制上运行快得多，大大提高了实时丰度估计的运行效率。

Description

一种基于GPU的丰度估计算法的并行实现方法

技术领域

本发明属于图像处理技术领域，具体为涉及一种面向高光谱图像的丰度值估计方法的GPU并行处理方法，以提高实时高光谱丰度的快速计算。

背景技术

高光谱遥感技术旨在以数十个到数百个连续的精细光谱波段获取地物信息，目前已被广泛地应用于环境、国防、农业以及深空等众多领域。丰度估计作为高光谱技术中最重要的主题之一，其目的是在端元提取完成之后，获得相应端元在每一个混合像元中所占的对应百分比(即丰度值)。传统算法(如最小二乘误差和正交子空间投影)的实现通常涉及大量矩阵求逆和乘法运算过程，使得其在软件运行上较慢，并且在硬件上也难以实现，因此不适用于对实时需求高的应用领域。

发明内容

本发明针对以上问题的提出，而研制一种基于GPU的丰度估计算法的并行实现方法。

本发明的技术手段如下：

一种基于GPU的丰度估计算法的并行实现方法，包括如下步骤：

步骤1、将高光谱混合像元数据H存储至GPU端，所述高光谱数据表示为H＝{B₁,B₂,,…B_i,…B_l,}，其中B_i＝(r₁ ^l,r₂ ^l,…r_N ^l)^T，N为像元数目，l为波段数目，r为高光谱数据中的高光谱影像中某个混合像元的光谱向量；

步骤2、自所述高光谱数据中的端元集M中获取每一个端元所对应的端元向量m_i同时将第i个端元m_i置于端元集M的最后位置，即端元集并基于CPU端的openMp机制，将所获取的端元向量m_i逐一分配至CPU端的p个并行线程上；

步骤3、基于p个并行线程，对端元向量m_i进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,…,m_i-1]的向量即

步骤4、获得端元集M所对应的正交基并计算对应的内积值即计算

步骤5、计算像元的光谱向量r在方向上的投影即计算以获取每个端元所对应的投影数据并存储至CPU端保存；

步骤6、计算各端元m_p对应的丰度结果α_p，并存储至GPU端。

进一步的，所述步骤1包括：

步骤11、依次读取高光谱混合像元数据H中第i个波段的所有像元并存储至GPU端，其中1<＝i<＝l；

步骤12、依次读取第i个端元m_i的所有波段值并存储至GPU端，其中1<＝i<＝p。

进一步的，所述步骤3包括：

步骤31、令第一个端元m₁的正交向量为端元向量即

步骤32、基于p个并行线程，依次对端元向量m_i对进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,…,m_i-1]的向量即m～_i，对应的计算公式为

步骤33、判断i是否等于p；是则执行步骤5；否则使得i←i+1，并执行步骤32。

进一步的，所述步骤6是指基于OVP算法计算各端元m_p的丰度结果α_p即基于下述公式计算，所述公式为：丰度结果α_p对应的函数

与现有技术相比，本发明的有益效果：

为了实现对高光谱图像丰度估计的有效性，本发明提供一种能够在不影响准确性的前提下，高光谱图像丰度估计实时性的并行处理方法，具体的其根据其适用要求为并行度高、高光谱数据量大的特点，基于并行编程模式，从算法逻辑、线程映射模型、数据读取存储传输方式和正确性等方面进行了并行化设计，大大提高了实时丰度估计的运行效率。

附图说明

图1是本发明所述并行方法所对应的核心步骤示意图；

图2是本发明所述实验采用模拟图像数据；

图3是表示模拟图像在5个端元情况下第1个端元在CPU串行平台下对应的丰度情况的灰度图；

图4是表示在本发明所述并行方法下第1个端元对应的丰度情况的灰度图；

图5是表示模拟图像在5个端元情况下第3个端元在CPU串行平台下对应的丰度情况的灰度图；

图6是表示在本发明所述并行方法下第3个端元对应的丰度情况的灰度图；

图7是表示模拟图像在5个端元情况下第5个端元在CPU串行平台下对应的丰度情况的灰度图；

图8是表示在本发明所述并行方法下第5个端元对应的丰度情况的灰度图；

图9本发明所述实验采用真实图像数据；

图10是表示真实图像在5个端元情况下第1个端元在CPU串行平台下对应的丰度情况的灰度图；

图11是表示真实图像在5个端元情况下第1个端元本发明所述并行方法下对应的丰度情况的灰度图；

图12是表示真实图像在5个端元情况下第3个端元在CPU串行平台下对应的丰度情况的灰度图；

图13是表示真实图像在5个端元情况下第3个端元本发明所述并行方法下对应的丰度情况的灰度图；

图14是表示真实图像在5个端元情况下第5个端元在CPU串行平台下对应的丰度情况的灰度图；

图15是表示真实图像在5个端元情况下第5个端元在本发明所述并行方法下对应的丰度情况的灰度图；

图16是表示OVP并行设计在模拟图像上加速比随着端元个数的增加而变化的情况；

图17是表示OVP并行设计在真实图像上加速比随着端元个数的增加而变化的情况；

图18是表示并行OVP、LSE以及LSOSP在模拟图像上端元个数从3到20的比较情况；

图19是表示并行OVP、LSE以及LSOSP在真实图像上端元个数从3到20的比较情况。

具体实施方式

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

如图1所示的一种基于GPU的丰度估计算法的并行实现方法，包括如下步骤：

步骤1、将高光谱混合像元数据H存储至GPU端即通过通信函数cudamemcpy将其拷贝到GPU显存中，所述高光谱数据表示为H＝{B₁,B₂,,…B_i,…B_l,}，其中B_i＝(r₁ ^l,₂ ^l,…r_N ^l)^T，N为像元数目，l为波段数目，r为高光谱数据中的高光谱影像中某个混合像元的光谱向量；

步骤2、因为施密特正交化(Gram-Schmidt)在处理每个端元互不干涉，为了求解第i个端元即m_i,则需要进行下述：即自所述高光谱数据中的端元集M中获取每一个端元所对应的端元向量m_i同时将第i个端元m_i置于端元集M的最后位置，即端元集并基于CPU端的openMp机制(处理器上编写并行程序的应用编程接口)，将所获取的端元向量m_i逐一分配至CPU端的p个并行线程上；即使得端元m_i分配到第i(1<＝i<＝p)个并行线程上；

步骤3、基于p个并行线程，对端元向量m_i进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,…,m_i-1]的向量即m～_i；

步骤4、从不同的端元开始，不同的端元顺序进行正交化获得端元集M空间上的一组正交基，它代表该向量正交于其它向量扩展空间的分量，则所述端元集M所对应的正交基并计算对应的内积值即计算具体GPU端建立GPU核函数1(Kernel)计算的内积值的步骤为：41、利用Cuda框架下的通信函数cudamemcpy建立GPU与CPU之间的联系，将CPU的数据即通过内存通信拷贝(利用cudaMalloc函数)至GPU显存中；42、利用CPU端的openMp机制给CPU分配p个线程,并行处理内积的结果,将得到的结果放入数组Q中；43、利用通信函数cudamemcpy建立GPU端与CPU端的联系，将计算结果数组Q通过GPU端返回到CPU端，在CPU段利用数组Q’(Q’＝Q)保存。

步骤5、计算像元的光谱向量r在方向上的投影即计算以获取每个端元所对应的投影数据并存储至CPU端保存；GPU端建立GPU核函数2计算计算的步骤为：51、利用cudaMalloc函数将数据拷贝到GPU内存中；

52、分配线程数目，设定线程数为图像端元数p与图像大小N之积，即p×N个线程；53、在GPU端的使得每个线程分别计算一个像元r在方向上的投影，并将线程数目分为p组，及GPU端的第i组线程(1<＝i<＝p)计算的是所有像元在方向上的投影，将投影值存储为C数组的第i个位置，即C＝[q_i1,q_1i2 … q_iN]；54、将计算结果通过GPU端返回到CPU端，此时结果集中存放着端元个数与图像大小之积个数据，同时在CPU端利用数组C’(C’＝C)保存。

步骤6、计算各端元m_p对应的丰度结果α_p，并存储至GPU端。具体的GPU端建立GPU核函数3计算计算丰度值的步骤为：61、分配线程数目，设定线程数目为图像大小与端元个数之积个线程数目，即p×N；62、利用GPU端每个线程i(0<＝i<p×N)计算Q数组和C数组行号和列号对应下标的值的乘积运算q_m×c_mn，m和n为二维数组的行号和列号；63、每N个线程的计算结果表示第m_i(1<＝i<＝p)个端元的解混结果所对应的丰度f_i，将f_i进行存储在F数组的第i个位置；64、通过GPU端返回到CPU端，此时结果矩阵F’中存放着各个端元对应的丰度值。

其中，方法所涉及的数据包括：知高光谱混合像元数据表示为H＝{B₁,B₂,,…B_i,…B_l,}，其中B_i＝(r₁ ^l,r₂ ^l,…r_N ^l)^T，N为像元数目，l为波段数目。r为高光谱影像中某个混合像元的光谱向量；M为l×p的端元矩阵，p为端元个数，l为高光谱波段数；若端元集M＝[m₁,m₂,…,m_p-1,m_p],则其中m_i＝[m_i1,m_i2,…m_il]^T，其对应的该空间上的一组正交基为其中利用F数组存储各端元的解混结果F＝[f₁,f₂,…f_i…f_p]，其中f_i表示第i个端元的丰度估计值。

进一步的，所述步骤1包括：

步骤11、依次读取高光谱混合像元数据H中第i个波段的所有像元并存储至GPU端，其中1<＝i<＝l；

步骤12、依次读取第i个端元m_i的所有波段值并存储至GPU端，其中1<＝i<＝p。

进一步的，所述步骤3包括：

步骤31、令第一个端元m₁的正交向量为端元向量即

步骤32、基于p个并行线程，依次对端元向量m_i对进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,…,m_i-1]的向量即对应的计算公式为

步骤33、判断i是否等于p；是则执行步骤5；否则使得i增加1即i←i+1，并执行步骤32。

进一步的，所述步骤6是指基于OVP算法计算各端元m_p的丰度结果α_p即基于下述公式计算，所述公式为：丰度结果α_p对应的函数

下面以具体的实例进行进一步说明：所述实例主要是在模拟图像以及真实图像上进行验证：

实验平台及实验数据

相关的实验平台搭建：实验机器为ThinkPad T430，硬件配置处理器为Intel(R)Core(TM)i5-2520、4核CPU、主频2.50GHz、内存6GB；显卡信息为NVIDIA NVS 5400M，显卡内存为1664M；软件平台中操作系统为Windows 10，开发环境Visual Studio 2013以及CUDA7.5。

相关的实验数据：本实验采用两组图像数据，其大小分别为350×350和400×400，各有189个波段；模拟图像如图2所示，真实图像如图9所示。

相关的实验条件设定：选择端元个数为依次为10、40和70时，测试OVP算法的CPU串行时间与CUDA+OpenMP并行时间的实验结果，具体如图3—图19，如下：

选择端元个数为依次为10、40和70时，测试OVP算法并行与LSE以及OSP并行时间的实验结果如下：

综上，为了实现对高光谱图像丰度估计的有效性，本发明提供一种能够在不影响准确性的前提下，高光谱图像丰度估计实时性的方法。本发明根据算法并行度高、高光谱数据量大的特点，依据并行思想,从算法逻辑、线程映射模型、数据读取存储传输方式和正确性等方面进行并行化设计,并用模拟高光谱数据以及真实高光谱数据进行了测试,可见在取得较好的丰度估计效果的同时,达到了较高加速比。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (4)

1.一种基于GPU的丰度估计算法的并行实现方法，包括如下步骤：

步骤1、将高光谱混合像元数据H存储至GPU端，所述高光谱数据表示为H＝{B₁，B_2,,…B_i,…B_l,}，其中B_i＝(r₁ ^l,r₂ ^l,…r_N ^l)^T，N为像元数目，l为波段数目，r为高光谱数据中的高光谱影像中某个混合像元的光谱向量；

步骤2、自所述高光谱数据中的端元集M中获取每一个端元所对应的端元向量m_i同时将第i个端元m_i置于端元集M的最后位置，即使得端元集并基于CPU端的openMp机制，将所获取的端元向量m_i逐一分配至CPU端的p个并行线程上；

步骤3、基于p个并行线程，对端元向量m_i进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,...,m_i-1]的向量即

步骤4、获得端元集M所对应的正交基并计算对应的内积值即计算

步骤5、计算像元的光谱向量r在方向上的投影即计算以获取每个端元所对应的投影数据并存储至CPU端保存；

步骤6、计算各端元m_p对应的丰度结果α_p，并存储至GPU端。

2.根据权利要求1所述的并行实现方法，其特征在于所述步骤1包括：

步骤11、依次读取高光谱混合像元数据H中第i个波段的所有像元并存储至GPU端，其中1<＝i<＝l；

步骤12、依次读取第i个端元m_i的所有波段值并存储至GPU端，其中1<＝i<＝p。

3.根据权利要求1所述的并行实现方法，其特征在于所述步骤3中步骤如下：

步骤31、令第一个端元m₁的正交向量为端元向量即

步骤32、基于p个并行线程，依次对端元向量m_i对进行施密特正交化处理以获得端元向量m_i其所对应的正交于向量空间[m₁,m₂,...,m_i-1]的向量即对应的计算公式为

步骤33、判断i是否等于p；是则执行步骤5；否则使得i←i+1，并执行步骤32。

4.根据权利要求1所述的并行实现方法，其特征在于所述步骤6包括：

所述步骤6是指基于OVP算法计算各端元m_p的丰度结果α_p即基于下述公式计算，所述公式为：丰度结果α_p对应的函数