CN107525588A

CN107525588A - 一种基于gpu的双相机光谱成像系统的快速重构方法

Info

Publication number: CN107525588A
Application number: CN201710699197.6A
Authority: CN
Inventors: 黄华; 张仕鹏; 王立志; 付莹
Original assignee: Beijing Institute of Technology BIT
Current assignee: Beijing Institute of Technology BIT
Priority date: 2017-08-16
Filing date: 2017-08-16
Publication date: 2017-12-29
Anticipated expiration: 2037-08-16
Also published as: CN107525588B

Abstract

本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，涉及能够快速获取高分辨率高光谱图像的方法，属于计算摄像学领域。本发明应用于基于编码孔径快照光谱成像+灰度相机的双相机光谱成像系统，将高光谱图像重构问题转化为多个子优化问题，并使用GPU完成各个子问题的求解：使用cuBLAS库和共轭梯度下降法更新高光谱图像；使用软阈值函数更新辅助变量；重复迭代完成高光谱图像的重构。本发明能够高质量地完成双相机光谱成像系统的高光谱图像重构，在保证重建结果具备高空间分辨率和高光谱保真性的同时，大幅度提高高光谱图像重建的效率，扩展高光谱图像的应用范围。本发明可用于载人航天、地质勘测和植被研究等多个领域。

Description

一种基于GPU的双相机光谱成像系统的快速重构方法

技术领域

本发明专利涉及一种用于双相机光谱成像的高光谱图像重构方法，尤其是能够快速获取高分辨率高光谱图像的方法，属于计算摄像学领域。

背景技术

高光谱成像技术是将空间成像技术和频谱成像技术相结合，对目标场景在频谱上进行几十甚至上百个波段的连续测量。该技术得到的图像包含目标场景的二维空间信息和一维光谱信息，被称为数据立方体。相比传统彩色成像，高光谱成像技术能够获得内容更为丰富，细节更为显著的有用信息。该技术已经被应用于植被研究、大气检测，载人航天、医学诊断等多个领域。

近年来，基于压缩感知的计算光谱成像技术被广泛运用到高光谱成像。与传统光谱成像系统相比，计算光谱成像能够获得更高空间分辨率、光谱分辨率和时间分辨率的高光谱图像，具有更为广阔的应用前景。Ashwin Wagadarikar等人提出的编码孔径快照光谱成像仪(Coded Aperture Snapshot Spectral Imager,CASSI)采用编码孔径和色散介质来对目标场景进行调制，通过探测器获取三维图谱数据的二维压缩投影。双相机光谱成像系统则是在CASSI系统的基础上，加入了灰度相机通道。该系统在获取目标场景的二维压缩光谱投影的同时，还获取目标场景的二维灰度投影。与最初的CASSI系统相比，双相机系统能够最大化入射光利用率，大幅度提高高光谱图像的成像质量，因此已成为该领域的重点研究对象。

目前，从双相机系统重构出原始的三维高光谱图像的方法主要有两类：基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/ThresholdingAlgorithms，TwIST)和基于自适应字典学习的重构算法。TwIST算法建立在迭代收缩阈值算法和迭代加权收缩算法的基础上，其核心是使用反向投影函数去噪，并利用前两次的结果进行迭代更新，从而完成高光谱图像的重建。由于其计算复杂度很高，TwIST算法重构高光谱图像时间很长，完全不能满足快速重建的要求。同时，重构的高光谱图像分辨率也有待提高。

基于自适应字典学习的重构算法利用光谱图像和灰度图像在空间维的相似性，采用字典学习和稀疏表示的方法完成高光谱图像的重构。这种方法虽然能够获得较好的重构质量，但是每次重构都需要重新学习字典，导致重构时间比TwIST算法更长，因此也不能满足快速重建的要求，极大地限制了高光谱图像的应用和发展。

CUDA(Compute Unified Device Architecture)，是显卡厂商NVIDIA推出的通用并行计算架构，该架构使用图像处理器GPU作为计算核心完成复杂的计算难题。相比传统CPU处理器，GPU拥有更多的计算单元，能够快速完成运算量大的计算任务，且已经在深度学习和人工智能等多个领域得到应用。

发明内容

针对现有重构算法存在的重构时间长、图像质量低等问题。本发明要解决的技术问题是提供一种基于GPU的双相机光谱成像系统的快速重构方法，能够快速完成双相机光谱成像系统的重构，具有重构速度快，成像质量高等优点。

为达到以上目的，本发明采用以下技术方案：

本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，应用于基于编码孔径快照光谱成像+灰度相加的双相机光谱成像系统，将高光谱图像重构问题转化为多个子优化问题，并使用图形处理器GPU完成各个子问题的求解：使用共轭梯度下降法更新高光谱图像；使用软阈值函数更新中间变量；重复迭代完成高光谱图像的重构。

本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，包括以下步骤：

步骤101：输入双相机光谱成像仪的采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ。

步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(CodedAperture Snapshot Spectral Imager,CASSI)+灰度相机的双相机采样光谱成像仪。双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机等部件构成。目标场景的高光谱图像X大小为M×N×Ω，高光谱图像X上任意一点的像素值为x(i,j,λ)，1≤i≤M,1≤j≤N,1≤λ≤Ω。其中，M×N表示高光谱图像的空间分辨率，Ω表示高光谱图像的频段数。入射光首先到达分光镜将其能量一分为二，一半进入编码孔径快照光谱成像仪CASSI分支，一半进入灰度相机分支。进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码。经编码后的图像到达色散棱镜后，不同频段的图像会沿着竖直方向进行偏移。最后所有频段的图像到达灰度相机后进行叠加，得到压缩的二维混叠光谱图像。编码孔径快照光谱成像仪CASSI成像仪的数学模型为：

其中ω(λ)表示CCD相机的频谱响应函数，Cu(i,j)表示编码模版函数，φ(λ)表示色散棱镜的波段偏移函数，y_c(i,j)为二维压缩光谱采样图像，v₁(i,j)表示高斯白噪声。

进入灰度相机分支的入射光会直接到达灰度相机得到目标场景的二维灰度投影，数学模型为：

将式(1)和式(2)联立，并写成矩阵的形式为：

Y＝HX+V (3)

其中Y＝[Y_c；Y_p]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体，V＝[V_c；V_p]表示高斯白噪声。H表示对编码孔径快照光谱成像仪CASSI+灰度相机标定后的前向响应矩阵，为编码模版函数、色散棱镜偏移函数和CCD灰度相机频谱响应的综合作用。

步骤102：申请主机CPU内存，将前向矩阵H和二维压缩光谱采样图像Y读入到主机CPU内存。

步骤102所述申请主机CPU内存优选C语言的malloc内存分配函数，申请的主机CPU内存为一维线性内存，方便后续处理中将CPU内存读入图像处理器GPU显存的操作。

步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定：如果色散棱镜是沿竖直方向偏移，则以二维采样图像的水平方向为基准进行读入，即逐行读入；如果色散棱镜是沿水平方向偏移，则以二维采样图像的竖直方向为基准进行读入，即逐列读入。目的是方便在后续处理中编码孔径快照光谱成像仪CASSI系统前向模型的计算，提高计算效率。

步骤103：为输入输出、以及后续计算过程所需要的数据申请GPU显存，将二维压缩光谱图像Y和前向响应矩阵H拷贝到GPU显存中。设置CUDA并行计算时需用的线程格/块数。

步骤103所述申请图形处理器GPU显存的方法优选通用并行计算架构CUDA的cudaMalloc函数，拷贝到图形处理器GPU显存的方法优选通用并行计算架构CUDA的cudaMemcpy函数。

步骤104：初始化重构高光谱图像X⁰，优化目标函数的初始值f⁰，以及辅助矩阵S和U，辅助因子ρ。初始化当前迭代次数t＝0。

步骤104所述重构高光谱图像X⁰的初始化方式为：

X⁰＝H^TY (4)

其中H^T表示编码孔径快照光谱成像仪CASSI系统前向响应H的转置，即从二维压缩光谱图像+二维灰度图像反演到三维数据立方体的过程。

步骤104所述优化目标函数为全局优化目标函数。根据自然图像的平滑特性，将高光谱重构问题转化为基于全变差约束的最优化问题，从而得到全局优化目标函数为：

其中运算符表示L2范数的平方。||DX||₁表示图像的全变差值，其定义如下：

||DX||₁＝Σ_i,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)

矩阵D为差分矩阵，它和图像X的作用结果如下：

公式(7)的结果为两个大小和X相同的矩阵，分别表示图像X在水平方向和竖直方向上的差分值。

步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量，大小均为2×M×N×Ω，初始化为全零矩阵。公式(5)无法直接求解，将公式(5)转化为多个子优化问题的求解，原理如下：

引入辅助矩阵S＝DX，得到如下带约束的优化方程：

其增广拉格朗日方程为：

其中ρ为辅助因子，优选2.048。当式(9)最小时，变量X、S的值即为优化问题式(8)的解，因此得到两个子优化问题，分别为：

接下来，对X、S和U进行交替更新、循环迭代即可完成高光谱图像的重构。

步骤105：使用共轭梯度下降算法更新高光谱图像X^t+1。

针对优化问题式(10)，X的最小二乘解为：

X^t+1＝(H^TH+ρD^TD)^-1(H^TY+D^T(U^t+ρS^t)) (12)

由于矩阵H和矩阵D的规模很大，无法直接求其解析解，因此需利用共轭梯度下降法求高光谱图像X^t+1的近似解，从而完成高光谱图像X^t+1的更新。

同时，为了加速程序的执行效率，借助NVIDIA提供的矩阵线性运算库cuBLAS来实现共轭梯度下降法，从而完成高光谱图像X^t+1的更新。

步骤106：更新辅助矩阵S和U。

根据优化问题式(11)，S的最小二乘解为：

上式为软阈值收缩函数，运算符表示点乘，S^t+1能够直接求出。

辅助矩阵U的更新公式如下：

U^t+1＝U^t+ρ(S^t-DX^t) (14)

U^t+1的更新方式优选矩阵线性运算库cuBLAS库，可以提高程序效率。

步骤107：使用更新后的高光谱图像X^t+1计算全局优化目标函数值f^t+1。

步骤107所述全局优化目标函数与步骤103的相同，即：

步骤108：根据步骤107计算的结果执行迭代选择策略，完成双相机光谱成像系统的快速重构。

步骤108所述迭代选择策略如下：如果收敛，即步骤107的目标函数值f^t+1小于上一个目标函数值f^t，则计算目标函数值的相对变化量：

若不满足迭代停止条件，即Tol大于预设阈值或者当前迭代次数t小于最大迭代次数Maxiters，则更新辅助因子ρ＝1.05×ρ，当前迭代次数t＝t+1，目标函数值f^t＝f^t+1，并转至步骤105进行迭代；否则，停止迭代并输出最后一次更新的高光谱图像，从而完成双相机光谱成像系统的快速重构。

如果不收敛，即步骤107的目标函数值f^t+1大于上一个目标函数值f^t，则更新辅助因子ρ＝2.0×ρ。判断更新后的ρ值，若ρ小于预设阈值，则转至步骤105进行迭代；否则，停止迭代并输出最后一次更新的高光谱图像，从而完成双相机光谱成像系统的快速重构。

有益效果：

1、本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，将高光谱重建问题转化为多个子优化问题进行求解，收敛速度和重构效率得到很大提升。

2、本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，使用交替迭代的方式求解优化问题，每一个更新步骤都能够达到最优解，因此能够提高空间分辨率和光谱保真度。

3、本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，使用GPU完成高光谱图像的重建过程，并利用cuBLAS库来提高GPU的利用率，大幅度降低双相机光谱成像仪的重建时间。

4、本发明公开的一种基于GPU的双相机光谱成像系统的快速重构方法，使用预先统一分配的方式管理GPU显存，最小化计算所需显存的大小，对GPU的要求低，使用入门级显卡就能够获得非常显著的加速效果。

附图说明

图1是本发明中用于双相机光谱成像的系统结构图；

图2是本发明公开的基于GPU的双相机光谱成像系统的快速重构方法的总流程图；

图3是本发明和对比算法对测试图片beads进行仿真重构后波长为440nm时的对比图，其中：图3-a为本发明重构的结果，图3-b为对比算法重构的结果。

图4是本发明和对比算法对测试图片beads进行仿真重构后波长为540nm时的对比图，其中：图4-a为本发明重构的结果，图4-b为对比算法重构的结果。

图5是本发明和对比算法对测试图片beads进行仿真重构后波长为640nm时的对比图，其中：图5-a为本发明重构的结果，图5-b为对比算法重构的结果。

具体实施方式

为了更好地说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。

实施例1：

本实施例公开的一种基于GPU的双相机光谱成像系统的快速重构方法，应用于基于编码孔径快照光谱成像+灰度相机的双相机光谱成像系统。编码孔径快照光谱成像系统，即CASSI系统，最早由Ashwin Wagadarikar等人提出(详见Wagadarikar A,John R,WillettR,Brady D.Single disperser design for coded aperture snapshot spectralimaging[J].Applied optics.2008,47(10):B44-B51.)。CASSI系统利用编码模版将光谱信息在空间维上进行随机调制，并使用色散棱镜在光谱维上进行压缩，从而实现了由三维数据立方体到二维数据的压缩。双相机光谱成像系统则是在CASSI系统的基础上，加入了灰度相机通道，其系统结构如图1所示(详见Wang L,Xiong Z,Gao D,et al.Dual-cameradesign for coded aperture snapshot spectral imaging[J].Applied Optics,2015,54(4):848-58.)。双相机光谱成像系统在获取目标场景的二维压缩光谱投影的同时，还获取目标场景的二维灰度投影。与CASSI系统相比，双相机系统能够最大化入射光利用率，同时大幅度提高高光谱图像的成像质量，因此已成为该领域的重点研究对象。

目前，从双相机系统重构出原始的三维高光谱图像的方法主要有两类：基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/ThresholdingAlgorithms，TwIST)和基于自适应字典学习的重构算法。TwIST算法(详见Bioucas-DiasJM,Figueiredo MAT.A New TwIST:Two-Step Iterative Shrinkage/ThresholdingAlgorithms for Image Restoration[J].IEEE Transactions on ImageProcessing.2007,16(12):2992-3004.)建立在迭代收缩阈值算法和迭代加权收缩算法的基础上，其核心是使用反向投影函数去噪(详见Chambolle A.An Algorithm for TotalVariation Minimization and Applications[M].Kluwer Academic Publishers,2004.)，并利用前两次的结果进行迭代更新，从而完成高光谱图像的重建。但是TwIST算法的计算复杂度较高，重构时间很长；同时，其重构质量也有待提高。

基于自适应字典学习的重构算法利用光谱图像和灰度图像在空间维的相似性，采用字典学习和稀疏表示的方法完成高光谱图像的重构(详见Wang L,Xiong Z,Gao D,etal.High-speed hyperspectral video acquisition with a dual-camera architecture[J].2015:4942-4950.)。这种方法虽然能够获得较好的重构质量，但是每次重构都需要重新学习字典，导致重构时间比TwIST算法更长，因此也不能满足快速重建的要求，极大地限制了高光谱图像的应用和发展。

针对现有重构算法存在的重构质量低、速度慢等问题，本实施例将高光谱图像重构问题转化为多个子优化问题，采用交替更新、重复迭代的方式完成高光谱图像的重构。同时，本实施例利用NVIDIA公司提出的通用并行计算架构(Compute Unified DeviceArchitecture，CUDA)来完成重构过程，极大地提高高光谱图像重构的效率。本实施例流程图如图2所示。

本实施例公开的一种基于GPU的双相机光谱成像系统的快速重构方法，包含以下步骤：

步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(CodedAperture Snapshot Spectral Imager,CASSI)+灰度相机的双通道采样光谱成像仪。双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机等部件构成。目标场景的高光谱图像X大小为M×N×Ω，高光谱图像X上任意一点的像素值为x(i,j,λ)，1≤i≤M，1≤j≤N，1≤λ≤Ω。其中，M×N表示高光谱图像的空间分辨率，Ω表示高光谱图像的频段数，入射光首先到达分光镜将其能量一分为二，一半进入编码孔径快照光谱成像仪CASSI分支，一半进入灰度相机分支。进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码。得到编码后的图像到达色散棱镜后，不同频段的图像会沿着竖直方向进行偏移。最后所有频段的图像到达CCD灰度相机后进行叠加，得到压缩的二维混叠光谱图像。CASSI成像仪的数学模型为：

其中ω(λ)表示CCD相机的频谱响应函数。Cu(i,j)表示编码模版函数。φ(λ)表示色散棱镜的波段偏移函数。y_c(i,j)为二维压缩光谱采样图像，v₁(i,j)表示高斯白噪声。

将式(1)和式(2)联立，并写成矩阵的形式为：

Y＝HX+V (2)

其中Y＝[Y_c；Y_p]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体，V＝[V_c；V_p]表示高斯白噪声。H表示对CASSI成像仪+灰度相机标定后的前向响应矩阵，为编码模版函数、CCD相机频谱响应和色散棱镜偏移函数的综合作用。

步骤102：使用malloc函数申请主机CPU内存，将前向矩阵H和二维压缩光谱采样图像Y读入到主机CPU内存。

步骤102所述malloc函数为C语言内存分配函数，申请的主机CPU内存为一维线性内存，方便后续处理中将CPU内存读入图形处理器GPU显存的操作。

步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定：如果色散棱镜是在竖直方向偏移，则以二维采样图像的水平方向为基准进行读入，即逐行读入；如果色散棱镜是在水平方向偏移，则以二维采样图像的竖直方向为基准进行读入，即逐列读入。目的是方便在后续处理中CASSI系统前向模型的计算，提高计算效率。

步骤103：使用cudaMalloc函数为输入输出、以及后续计算过程所需要的数据申请GPU显存，并使用cudaMemcpy函数将二维压缩光谱图像Y和前向响应矩阵H拷贝到GPU显存中。设置通用并行计算架构CUDA计算时需用的线程格/块数。

步骤103所述线程格/块数的设置方法为：线程格/块均为二维布局，单位线程块中的线程个数为(16,16)，单位线程格中线程块的个数则是根据高光谱图像的空间分辨率确定为((M+15)/16,(N+15)/16)。

步骤104所述重构高光谱图像X⁰的初始化方式为：

X⁰＝H^TY (4)

其中H^T表示CASSI系统前向响应H的转置，即从二维压缩光谱图像反演到三维数据立方体的过程。

||DX||₁＝∑_i,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)

矩阵D为差分矩阵，它和图像X的作用结果如下：

步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量，大小均为2×M×N×Ω，且都是用cudaMemset函数将其初始化为全零矩阵。公式(5)无法直接求解，将公式(5)转化为多个子优化问题的求解，原理如下：

引入辅助矩阵S＝DX，得到如下带约束的优化方程：

其增广拉格朗日方程为：

其中ρ为辅助因子，本实施例取2.048。当式(9)最小时，变量X、S的值即为优化问题式(8)的解，因此得到两个子优化问题，分别为：

步骤105：使用cuBLAS和共轭梯度下降算法更新高光谱图像X^t+1。

根据优化问题式(10)，X的最小二乘解为：

X^t+1＝(H^TH+ρD^TD)-¹(H^TY+D^T(U^t+ρS^t)) (12)

由于矩阵H和矩阵D的规模很大，无法直接求其解析解，因此需利用共轭梯度下降法求高光谱图像X^t+1的近似解(详见Hestenes M R,Stiefel E.Methods of ConjugateGradients for Solving Linear Systems[J].Journal of Research of the NationalBureau of Standards,1952,49(6):409-436.)。同时，为了加速程序的执行效率，本实施例借助NVIDIA提供的矩阵线性运算库cuBLAS来实现共轭梯度下降法，从而完成高光谱图像X^t ⁺¹的更新。

步骤106：更新辅助矩阵S和U。

根据优化问题式(11)，S的最小二乘解为：

辅助矩阵U的更新公式如下：

U^t+¹＝U^t+ρ(S^t-DX^t) (14)

使用cuBLAS库，U^t+1能够快速求得。

步骤107所述全局优化目标函数与步骤103的相同，即：

为说明本发明的效果，本实施例将在实验条件相同的情况下对两种方法进行对比。

1.实验条件

本实验的硬件测试条件为：Inter i5 3210，内存8G，Matlab 2015b。GPU为GT630M，显存2G，CUDA 7.0。测试所用高光谱图片来自于哥伦比亚大学的Cave数据集。输入的CASSI分支压缩光谱采样图像大小为512×542；灰度相机分支的灰度采样图像大小为512x512；重构以后得到的高光谱图像大小为512×512×31。编码孔径模版为贝努利随机矩阵。主循环最大迭代次数为10次。辅助因子ρ＝2.048。共轭梯度下降算法的迭代次数为25次，迭代停止门限值为1.0×10^-4。对比算法为基于全变差约束的TwIST算法，其最大迭代次数为200次。两种算法使用的正则化系数均为τ＝0.1，预设迭代收敛阈值均为1.0×10^-6。

2.实验结果

为了验证本发明的可行性，测试在不同噪声方差sigma下，本发明公开的算法和TwIST算法的重构性能。

首先是两种算法的重构PSNR，单位为dB，结果如表1所示：

表1 PSNR/dB

从表1的结果可以看出，本发明公开的算法迭代10次就能够达到非常好的重构效果，在不同强度的噪声影响下其PSNR均高于TwIST算法。

图3-a、图4-a和图5-a分别为测试图片beads使用本发明仿真重构后，波长分别为440nm、540nm和640nm时的结果；图3-b、图4-b和图5-b分别为使用TwIST算法仿真重构后对应波长的结果。可以看出，本发明公开的方法重构的结果更加清楚，在视觉上的效果要优于TwIST算法。

为了对比两种算法的光谱保真度，其均方根误差RMSE结果如表3所示：

表2 RMSE

从表2的结果可以看出，本发明公开的方法得到的均方根误差更小，表明其在频谱维上的保真度更高，其谱线更接近于真实材料的谱线。

接下来比较两种算法的重构时间，结果如表3所示：

表3重构时间/秒

从表3的结果可以看出，本发明公开的方法重构所用时间非常短，重建时间只需要13秒左右，而TwIST算法重建时间长达1833秒，因此重建速度得到非常大的提升。

此外，为了进一步说明本发明公开的方法在重建效率上的提升，本实施例测试了三种不同高光谱图像分辨率下，两种算法的平均重构时间和加速比，结果如下：

表4不同分辨率下平均重构时间对比

从表4可以看出，在不同分辨率下，本发明公开的方法都能够获得非常高的加速比。同时，随着分辨率的提升，重构计算量成倍增大，本发明公开的方法的加速比也在增加。以分辨率为1340×1018×33的重构图像为例，TwIST算法重建的时间长达2小时多，而使用本发明公开的方法只需要1分钟，其加速效果非常显著。更为重要的，本实施例使用的GPU为GT 630M，属于入门级显卡，如果使用更好的显卡将能够获得更高的加速比。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：包括以下步骤，

步骤101：输入双相机光谱成像仪的采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ；

步骤102：申请主机CPU内存，将前向矩阵H和二维压缩光谱采样图像Y读入到主机CPU内存；

步骤103：为输入输出、以及后续计算过程所需要的数据申请GPU显存，将二维压缩光谱图像Y和前向响应矩阵H拷贝到GPU显存中；设置CUDA并行计算时需用的线程格/块数；

步骤104：初始化重构高光谱图像X⁰，优化目标函数的初始值f⁰，以及辅助矩阵S和U，辅助因子ρ；初始化当前迭代次数t＝0；

步骤105：使用共轭梯度下降算法更新高光谱图像X^t+1；

针对优化问题式(10)，X的最小二乘解为：

X^t+1＝(H^TH+ρD^TD)^-1(H^TY+D^T(U^t+ρS^t)) (12)

由于矩阵H和矩阵D的规模很大，无法直接求其解析解，因此需利用共轭梯度下降法求高光谱图像X^t+1的近似解，从而完成高光谱图像X^t+1的更新；

步骤106：更新辅助矩阵S和U；

根据优化问题式(11)，S的最小二乘解为：

<mrow> <msup> <mi>S</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mi>m</mi> <mi>a</mi> <mi>x</mi> <mo>{</mo> <msup> <mi>DX</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mfrac> <msup> <mi>U</mi> <mi>t</mi> </msup> <mi>&rho;</mi> </mfrac> <mo>-</mo> <mfrac> <mi>&tau;</mi> <mi>&rho;</mi> </mfrac> <mo>}</mo> <mo>&CircleTimes;</mo> <mi>s</mi> <mi>i</mi> <mi>g</mi> <mi>n</mi> <mrow> <mo>(</mo> <msup> <mi>DX</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mfrac> <msup> <mi>U</mi> <mi>t</mi> </msup> <mi>&rho;</mi> </mfrac> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>

上式为软阈值收缩函数，运算符表示点乘，S^t+1能够直接求出；

辅助矩阵U的更新公式如下：

U^t+1＝U^t+ρ(S^t-DX^t) (14)

步骤107：使用更新后的高光谱图像X^t+1计算全局优化目标函数值f^t+1；

步骤107所述全局优化目标函数与步骤103的相同，即：

2.如权利要求1所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤101中所述双相机光谱成像仪为基于编码孔径快照光谱成像仪(Coded ApertureSnapshot Spectral Imager,CASSI)+灰度相机的双相机采样光谱成像仪；双相机光谱成像仪主要由分光镜、物镜、编码模版、中继镜、色散棱镜和灰度相机等部件构成；目标场景的高光谱图像X大小为M×N×Ω，高光谱图像X上任意一点的像素值为x(i,j,λ)，1≤i≤M,1≤j≤N,1≤λ≤Ω；其中，M×N表示高光谱图像的空间分辨率，Ω表示高光谱图像的频段数；入射光首先到达分光镜将其能量一分为二，一半进入编码孔径快照光谱成像仪CASSI分支，一半进入灰度相机分支；进入编码孔径快照光谱成像仪CASSI分支的光会到达编码模版进行0-1编码；经编码后的图像到达色散棱镜后，不同频段的图像会沿着竖直方向进行偏移；最后所有频段的图像到达灰度相机后进行叠加，得到压缩的二维混叠光谱图像；编码孔径快照光谱成像仪CASSI成像仪的数学模型为：

<mrow> <msub> <mi>y</mi> <mi>c</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>0.5</mn> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>&lambda;</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>&Omega;</mi> </msubsup> <mi>&omega;</mi> <mrow> <mo>(</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mi>C</mi> <mi>u</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mi>&phi;</mi> <mo>(</mo> <mi>&lambda;</mi> <mo>)</mo> <mo>)</mo> </mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>-</mo> <mi>&phi;</mi> <mo>(</mo> <mi>&lambda;</mi> <mo>)</mo> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>v</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>

其中ω(λ)表示CCD相机的频谱响应函数，Cu(i,j)表示编码模版函数，φ(λ)表示色散棱镜的波段偏移函数，y_c(i,j)为二维压缩光谱采样图像，v₁(i,j)表示高斯白噪声；

<mrow> <msub> <mi>y</mi> <mi>p</mi> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>=</mo> <mn>0.5</mn> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>&lambda;</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>&Omega;</mi> </msubsup> <mi>&omega;</mi> <mrow> <mo>(</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>v</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>

将式(1)和式(2)联立，并写成矩阵的形式为：

Y＝HX+V (3)

其中Y＝[Y_c；Y_p]表示二维压缩光谱采样图像+二维灰度图像；X表示三维数据立方体，V＝[V_c；V_p]表示高斯白噪声；H表示对编码孔径快照光谱成像仪CASSI+灰度相机标定后的前向响应矩阵，为编码模版函数、色散棱镜偏移函数和CCD灰度相机频谱响应的综合作用。

3.如权利要求2所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤102所述申请主机CPU内存优选C语言的malloc内存分配函数，申请的主机CPU内存为一维线性内存，方便后续处理中将CPU内存读入图像处理器GPU显存的操作。

4.如权利要求3所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤102所述将矩阵H和Y读入主机CPU内存的方式是根据编码孔径光谱仪的色散棱镜偏移方向确定：如果色散棱镜是沿竖直方向偏移，则以二维采样图像的水平方向为基准进行读入，即逐行读入；如果色散棱镜是沿水平方向偏移，则以二维采样图像的竖直方向为基准进行读入，即逐列读入；目的是方便在后续处理中编码孔径快照光谱成像仪CASSI系统前向模型的计算，提高计算效率。

5.如权利要求4所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤103所述申请图形处理器GPU显存的方法选通用并行计算架构CUDA的cudaMalloc函数，拷贝到图形处理器GPU显存的方法选通用并行计算架构CUDA的cudaMemcpy函数。

6.如权利要求5所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤104所述重构高光谱图像X⁰的初始化方式为：

X⁰＝H^TY (4)

其中H^T表示编码孔径快照光谱成像仪CASSI系统前向响应H的转置，即从二维压缩光谱图像+二维灰度图像反演到三维数据立方体的过程；

步骤104所述优化目标函数为全局优化目标函数；根据自然图像的平滑特性，将高光谱重构问题转化为基于全变差约束的最优化问题，从而得到全局优化目标函数为：

其中运算符表示L2范数的平方；||DX||₁表示图像的全变差值，其定义如下：

||DX||₁＝∑_i,j,λ|x(i,j+1,λ)-x(i,j,λ)|+|x(i+1,j,λ)-x(i,j,λ)| (6)

矩阵D为差分矩阵，矩阵D和图像X的作用结果如下：

<mrow> <mi>D</mi> <mi>X</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>x</mi> <mrow> <mo>(</mo> <mi>i</mi> <mo>,</mo> <mi>j</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>

公式(7)的结果为两个大小和X相同的矩阵，分别表示图像X在水平方向和竖直方向上的差分值；

步骤104所述辅助矩阵S和U是指求解优化方程(4)时所需要的辅助变量，大小均为2×M×N×Ω，初始化为全零矩阵；公式(5)无法直接求解，将公式(5)转化为多个子优化问题的求解，原理如下：

引入辅助矩阵S＝DX，得到如下带约束的优化方程：

其增广拉格朗日方程为：

其中ρ为辅助因子，优选2.048；当式(9)最小时，变量X、S的值即为优化问题式(8)的解，因此得到两个子优化问题，分别为：

7.如权利要求6所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤106中U^t+1的更新方式选矩阵线性运算库cuBLAS库，能够提高程序效率。

8.如权利要求1或2所述的一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：步骤108所述迭代选择策略如下：如果收敛，即步骤107的目标函数值f^t+1小于上一个目标函数值f^t，则计算目标函数值的相对变化量：

若不满足迭代停止条件，即Tol大于预设阈值或者当前迭代次数t小于最大迭代次数Maxiters，则更新辅助因子ρ＝1.05×ρ，当前迭代次数t＝t+1，目标函数值f^t＝f^t+1，并转至步骤105进行迭代；否则，停止迭代并输出最后一次更新的高光谱图像，从而完成双相机光谱成像系统的快速重构；

如果不收敛，即步骤107的目标函数值f^t+1大于上一个目标函数值f^t，则更新辅助因子ρ＝2.0×ρ；判断更新后的ρ值，若ρ小于预设阈值，则转至步骤105进行迭代；否则，停止迭代并输出最后一次更新的高光谱图像，从而完成双相机光谱成像系统的快速重构。

9.一种基于GPU的双相机光谱成像系统的快速重构方法，其特征在于：应用于基于编码孔径快照光谱成像+灰度相加的双相机光谱成像系统，将高光谱图像重构问题转化为多个子优化问题，并使用图形处理器GPU完成各个子问题的求解：使用共轭梯度下降法更新高光谱图像；使用软阈值函数更新中间变量；重复迭代完成高光谱图像的重构。