CN109146787B - 一种基于插值的双相机光谱成像系统的实时重建方法 - Google Patents

Abstract

本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，涉及能够实时获取高分辨率高光谱图像的方法，属于计算摄像学领域。本发明实现方法如下：应用于基于全色相机的双相机光谱成像系统，将高光谱重建分成低分辨率高光谱重建、上采样插值和高分辨率高光谱重建三个阶段，并根据系统原理和GPU工作特点建立快速插值模型、快速响应模型和快速差分模型来完成上述重建的三个阶段，能够在保证重建结果具备高空间分辨率和高光谱保真性的同时，大幅度提高高光谱图像重建的效率，达到实时重建高光谱图像的目的，极大地扩展高光谱图像的应用范围。本发明可用于地质勘探、生物研究、和人工智能等多个领域。

Description

一种基于插值的双相机光谱成像系统的实时重建方法

技术领域

本发明涉及一种用于双相机光谱成像系统的高光谱图像重建方法，尤其涉及一种能够实时获取高分辨率高光谱图像的方法，属于计算摄像学领域。

背景技术

得益于成像技术和光谱技术的成熟，高光谱成像技术能够获取目标场景的连续光谱曲线。该技术得到的图像包含目标场景的二维空间信息和一维光谱信息，被称为数据立方体。相比传统彩色成像，高光谱成像技术能够获得内容更为丰富，细节更为显著的有用信息。由于不同物质反射出不同的光谱曲线，高光谱图像能够很好地完成分类和识别等任务。目前该技术已经在地质勘探、农业生产和生物医学等众多领域得到广泛的应用。

近年来，快照式光谱成像系统将计算和成像相结合，能够完成高分辨率高光谱图像的获取，具有非常广阔的应用前景。Ashwin Wagadarikar等人提出的编码孔径快照光谱成像仪(Coded Aperture Snapshot Spectral Imager,CASSI)使用编码孔径和色散介质来对三维高光谱图像分别进行空间和光谱调制，通过探测器获取二维混叠投影。基于全色相机的双相机光谱成像系统(Panchromatic camera based dual camera compressivehyperspectral imaging，PDCCHI)则是在获取目标场景的二维混叠投影的同时，使用一个全色相机来获取目标场景的二维灰度投影。与CASSI系统相比，PDCCHI能够最大化入射光利用率，大幅度提高高光谱图像的成像质量，已成为该领域的研究重点。

如何从PDCCHI系统得到的二维压缩图像重建出原始的三维高光谱图像，是该领域非常重要的研究内容。目前常用于PDCCHI系统光谱重建的方法主要是基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/Thresholding Algorithms，TwIST)和自适应稀疏重建。TwIST算法建立在迭代收缩阈值算法和迭代加权收缩算法的基础上，其核心是使用反向投影函数去噪，并利用前两次的结果进行迭代更新，从而完成高光谱图像的重建。由于其计算复杂度很高，TwIST算法重建高光谱图像时间非常长，无法满足实时重建的要求。针对PDCCHI系统，Wang等人指出单个谱带的高光谱图像与全色分支的灰度采样具有很高的空间内容相似性，因此使用灰度采样训练自适应过完备字典，完成PDCCHI系统的自适应稀疏重建，从而提升重建质量。但由于自适应稀疏重建算法包含字典训练和稀疏表示两个阶段，因此重建时间比TwIST算法更长，也不能满足实时重建的要求，极大地限制了高光谱图像的应用和发展。

图像处理器GPU使用成百上千个计算单元并行地处理大规模计算任务，能够获得比传统CPU更为高效的工作性能。CUDA(Compute Unified Device Architecture)，是NVIDIA公式推出的通用并行计算架构，能够为用户提供更加高效和快捷的GPU控制技术，目前已经在图像处理和人工智能等多个领域得到应用。

发明内容

针对现有重建方法存在的重建时间长、无法实时获取高光谱图像等问题。本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法要解决的技术问题是：实时完成基于全色相机的双相机光谱成像系统的重建，具有重建速度快，重建质量高等优点。

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

本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，应用于基于全色相机的双相机光谱成像系统，将高光谱重建分成低分辨率高光谱重建、上采样插值和高分辨率高光谱重建三个阶段，并根据系统原理和GPU工作特点建立快速插值模型、快速响应模型和快速差分模型来完成上述重建的三个阶段，能够在保证重建结果具备高空间分辨率和高光谱保真性的同时，大幅度提高高光谱图像重建的效率，达到实时重建高光谱图像的目的。

本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，包括以下步骤：

步骤101：输入双相机光谱成像系统的采样图像Y、标定后的前向响应矩阵H、正则化系数τ，下采样比例β，低分辨率重建迭代次数I_d，高分辨率重建迭代次数I_h。

步骤101中所述双相机光谱成像仪为基于全色相机的双相机光谱成像系统(Panchromatic camera based dual camera compressive hyperspectral imaging，PDCCHI)。双相机光谱成像系统主要由分光镜、物镜、编码模版、中继镜、色散棱镜和全色相机部件构成，包含一个编码孔径快照光谱成像系统(Coded Aperture Snapshot SpectralImager,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)为二维混叠采样图像。

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

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

Y＝HX (3)

其中Y＝[Y^c；Y^p]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体，H表示对编码孔径快照光谱成像系统CASSI分支+全色分支的前向响应矩阵，为编码模版函数、色散棱镜偏移函数和CCD相机频谱响应的综合作用。

步骤102：使用快速插值模型对采样图像Y和前向响应矩阵H进行下采样，得到低分辨率采样Y_d和低分辨率前向响应矩阵H_d，空间和光谱的下采样因子均为β。

步骤102所述快速插值模型算法为统一化插值算法。线性图像插值过程都能够看作是源图像像素的线性叠加过程，故用统一的数学模型表示为：

其中d和s分别表示目标图像和源图像，j表示像素点索引。K表示插值所需要的源像素点个数。I(j,m)和W(j,m)表示索引矩阵和权值矩阵。一旦图像分辨率和采样比例确定以后，I(j,m)和W(j,m)就能够提前生成并作为查找表存储在内存中，从而完成快速图像插值模型。

步骤102所述对采样图像Y和前向响应矩阵H进行下采样的方法如下：对于采样图像Y，使用最近邻插值对其CASSI混叠采样部分下采样，使用均值滤波对其全色采样部分进行下采样。同时，使用均值滤波对前向响应矩阵进行下采样。

步骤103：初始化低分辨率高光谱图像

以及低分辨率辅助矩阵

和

步骤103所述重建低分辨率高光谱图像

的初始化方式为：

其中

表示下采样得到的编码孔径快照光谱成像系统的前向响应H_d的转置，

表示由下采样二维观测反演到低分辨率高光谱图像。由于空间和光谱均进行了比例为β的下采样，因此低分辨率高光谱图像

的大小为βM×βN×βΩ。

步骤103所述辅助矩阵S_d和U_d为重建低分辨率高光谱图像时需要用到的辅助矩阵，大小均为2×βM×βN×βΩ，初始化为全零矩阵。

步骤104：初始化优化目标函数的初始值

辅助因子ρ。初始化当前迭代次数t＝0。

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

其中运算符

表示矩阵Frobenius范数的平方。||DX_d||₁表示图像的全变差值，其定义如下：

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

公式(7)的结果为两个大小和X_d相同的矩阵，分别表示图像X_d在水平方向和竖直方向上的差分值。公式(6)无法直接求解，将公式(6)转化为多个子优化问题的求解，具体求解方法如下：

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

公式(9)的增广拉格朗日方程为：

由公式(10)得到两个子优化问题，分别为：

其中t表示迭代次数。接下来，对X_d、S_d和U_d进行交替更新、循环迭代即完成高光谱图像的重建。

步骤105：使用快速响应模型和快速差分模型更新高光谱图像X_d。

针对优化问题式(10)，低分辨率高光谱图像X_d的最小二乘解为：

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

步骤105所述快速响应模型是为了加速执行PDCCHI系统的前向响应过程，也就是公式(3)所描述的PDCCHI获取采样图像的过程。使用查找表来实现快速响应模型，查找表的大小为2×(βN+βΩ-1)，其定义方式为：

使用公式(13)建立的查找表来重写CASSI前向模型即公式(1)，则为：

其中ω_d(λ)和Cu_d(i,j)分别为ω(λ)和Cu(i,j)经过下采样的结果。由于在使用GPU对前向模型实现并行加速时，公式(1)会出现边界判断从而大幅度降低并行加速的整体性能。而上述查找表优化方法能够降低这一影响，显著提高并行优化性能。

步骤105所述快速差分模型是为了优化前向差分过程即DX_d。使用共享内存对前向差分所用数据进行预加载，计算DX_d时数据直接从共享内存读取，从而完成快速差分模型。由于在计算DX_d时会频繁访问GPU的全局内存区，会降低并行数据的访存效率。利用共享内存的快速访问特点来对前向差分过程进行优化，能够提高前向差分的实施效率。

步骤106：使用快速差分模型更新辅助矩阵S_d和U_d。

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

上式为软阈值收缩函数，运算符

表示点乘。使用步骤105所述快速差分模型，

能够直接求出。

辅助矩阵

的更新公式如下：

步骤107：更新辅助因子ρ和当前迭代次数t，并转至步骤105进行迭代，直至t＝I₁，完成低分辨率高光谱图像X_d的重建。

步骤108：对X_d，S_d和U_d均使用快速插值模型进行双线性上采样插值，插值比例为1/β，插值结果分别为高分辨率高光谱图像X_u,高分辨率辅助矩阵S_u和U_u。

步骤108所述快速插值模型为步骤102所述的快速插值模型。

步骤109：初始化高分辨率高光谱图像X⁰＝X_u，高分辨率辅助矩阵S⁰＝S_u，U⁰＝U_u，辅助因子ρ，迭代次数t＝0。

步骤110：使用快速响应模型和快速差分模型更新高分辨率高光谱图像X，其更新方式为：

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

步骤111：使用快速差分模型更新高分辨率辅助矩阵S和U，更新方式为：

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

步骤112：更新辅助因子ρ和当前迭代次数t，并转至步骤110进行迭代，进行高分辨率高光谱图像X的重建，直至t＝I₂，从而完成双相机光谱成像系统的实时重建，得到高质量的实时高光谱重建图像。

有益效果：

1、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，将高光谱图像重建分为低分辨率高光谱图像重建、上采样插值和高分辨率高光谱图像重建三个阶段进行，并建立快速插值模型、快速响应模型和快速差分模型，提高重建效率、能够达到实时的高光谱图像重建，得到高质量的实时重建高光谱图像。

2、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，先重建低分辨率高光谱图像能够快速恢复原始分辨率高光谱图像的结构信息，保证重建的整体质量。

3、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，对低分辨率重建结果进行插值并作为恢复原始分辨率高光谱图像的初值，能够加速重建高分辨率高光谱图像的收敛速度。

4、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，建立统一形式的快速插值模型，并使用查找表来提升下采样过程和上采样过程的执行效率。

5、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，根据重建原理和GPU的使用特性，建立快速响应模型和快速差分模型等并行优化，提升重建高光谱图像的整体效率。

6、本发明公开的一种基于插值的双相机光谱成像系统的实时重建方法，重建速度高且重建高光谱图像的质量高，适用于生物研究、地质勘探和人工智能等多个领域。

附图说明

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

图2是本发明公开的基于插值的双相机光谱成像系统的实时重建方法的总流程图；

图3是本发明提出的用于快速差分模型的共享内存结构。

图4是对比算法和本发明提出的算法在1秒之内能够重建不同分辨率高光谱图像的帧数。

图5是本发明和对比算法对测试图片eve-1551进行仿真重建后波长为610nm时的对比图，其中：图5-a为参考图片，图5-b为对比算法的重建结果，图5-c为本发明提出的算法的重建结果。

图6是本发明和对比算法对测试图片lst-0950进行仿真重建后波长为610nm时的对比图，其中：图6-a为参考图片，图6-b为对比算法的重建结果，图6-c为本发明提出的算法的重建结果。

具体实施方式

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

实施例1：

本实施例公开的一种基于插值的双相机光谱成像系统的实时重建方法，应用于基于全色相机的双相机光谱成像系统(Panchromatic camera based dual cameracompressive hyperspectral imaging，PDCCHI)(详见Wang L,Xiong Z,Gao D,etal.Dual-camera design for coded aperture snapshot spectral imaging[J].AppliedOptics,2015,54(4):848-58.)。PDCCHI系统由一个编码孔径快照光谱成像系统(CodedAperture Snapshot Spectral Imager,CASSI)(详见Wagadarikar A,John R,Willett R,Brady D.Single disperser design for coded aperture snapshot spectral imaging[J].Applied optics.2008,47(10):B44-B51.)和一个全色分支组成。CASSI系统使用编码孔径和色散介质来对三维高光谱图像分别进行空间和光谱调制，通过探测器获取二维混叠投影。而PDCCHI在获取二维混叠投影的同时，使用一个全色相机来获取目标场景的二维灰度投影。与CASSI系统相比，PDCCHI能够最大化入射光利用率，大幅度提高高光谱图像的成像质量，因此已成为该领域的研究重点。

如何从PDCCHI系统得到的二维压缩图像重建出原始的三维高光谱图像，是该领域非常重要的研究内容。目前常用于PDCCHI系统光谱重建的方法主要是基于全变差最小化的两步迭代收缩/阈值算法(Two-Step Iterative Shrinkage/Thresholding Algorithms，TwIST)和自适应稀疏重建。TwIST算法(详见Bioucas-Dias JM,Figueiredo MAT.A NewTwIST:Two-Step Iterative Shrinkage/Thresholding Algorithms for ImageRestoration[J].IEEE Transactions on Image Processing.2007,16(12):2992-3004.)使用反向投影函数去噪(详见Chambolle A.An Algorithm for Total VariationMinimization and Applications[M].Kluwer Academic Publishers,2004.)，并利用前两次的结果进行迭代更新，从而完成高光谱图像的重建。由于其计算复杂度很高，TwIST算法重建高光谱图像时间非常长，无法满足实时重建的要求。针对PDCCHI系统，Wang等人指出单个谱带的高光谱图像与全色分支的灰度采样具有很高的空间内容相似性，因此使用灰度采样训练自适应过完备字典，完成PDCCHI系统的自适应稀疏重建(详见Wang L,Xiong Z,GaoD,et al.High-speed hyperspectral video acquisition with a dual-cameraarchitecture[J].2015:4942-4950.)，从而提升重建质量。但由于自适应稀疏重建算法包含字典训练和稀疏表示两个阶段，因此重建时间比TwIST算法更长，因此也不能满足实时重建的要求，极大地限制了高光谱图像的应用和发展。

针对现有重建方法存在的重建时间长、无法实时获取高光谱图像等问题，本实施例公开的一种基于插值的双相机光谱成像系统的实时重建方法，应用于基于全色相机的双相机光谱成像系统，将高光谱重建分成低分辨率高光谱重建，上采样插值和高分辨率高光谱重建三个阶段，并根据系统原理和GPU工作特点建立快速插值模型、快速响应模型和快速差分模型来完成上述重建的三个阶段，能够在保证重建结果具备高空间分辨率和高光谱保真性的同时，大幅度提高高光谱图像重建的效率，达到实时重建高光谱图像的目的。本实施例流程图如图2所示。

本实施例公开的一种基于插值的双相机光谱成像系统的实时重建方法，包括以下步骤：

步骤101：输入双相机光谱成像系统的采样图像Y、标定后的前向响应矩阵H、正则化系数τ，下采样比例β，低分辨率重建迭代次数I_d，高分辨率重建迭代次数I_h。

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

其中ω(λ)表示CCD相机的频谱响应函数，Cu(i,j)表示编码模版函数，φ(λ)表示色散棱镜的波段偏移函数，y^c(i,j)为二维混叠采样图像。

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

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

Y＝HX (3)

其中Y＝[Y^c；Y^p]表示二维压缩光谱采样图像+二维灰度图像。X表示三维数据立方体，H表示对编码孔径快照光谱成像系统CASSI分支+全色分支的前向响应矩阵，为编码模版函数、色散棱镜偏移函数和CCD相机频谱响应的综合作用。

步骤102：使用快速插值模型对采样图像Y和前向响应矩阵H进行下采样，得到低分辨率采样Y_d和低分辨率前向响应矩阵H_d，空间和光谱的下采样因子均为β。

步骤102所述快速插值模型为统一化插值算法。常见的线性图像插值过程都能够看作是源图像像素的线性叠加过程，故用统一的数学模型表示为：

其中d和s分别表示目标图像和源图像，j表示像素点索引。K表示插值所需要的源像素点个数。I(j,m)和W(j,m)表示索引矩阵和权值矩阵。一旦图像分辨率和采样比例确定以后，I(j,m)和W(j,m)就能够提前生成并作为查找表存储在内存中，从而完成快速图像插值模型。

步骤102所述对采样图像Y和前向响应矩阵H进行下采样的方法如下。对于采样图像Y，使用最近邻插值对其CASSI混叠采样部分下采样，使用均值滤波对其全色采样部分进行下采样。同时，使用均值滤波对前向响应矩阵进行下采样。

步骤103：初始化低分辨率高光谱图像

以及低分辨率辅助矩阵

和

步骤103所述重建低分辨率高光谱图像

的初始化方式为：

其中

表示下采样得到的编码孔径快照光谱成像系统的前向响应H_d的转置，

表示由下采样二维观测反演到低分辨率高光谱图像。由于空间和光谱均进行了比例为β的下采样，因此低分辨率高光谱图像

的大小为βM×βN×βΩ。

步骤103所述低分辨率辅助矩阵S_d和U_d为重建低分辨率高光谱图像时需要用到的辅助矩阵，大小均为2×βM×βN×βΩ，初始化为全零矩阵。

步骤104：初始化优化目标函数的初始值

辅助因子ρ。初始化当前迭代次数t＝0。

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

其中运算符

表示矩阵Frobenius范数的平方。||DX_d||₁表示图像的全变差值，其定义如下：

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

公式(7)的结果为两个大小和X_d相同的矩阵，分别表示图像X_d在水平方向和竖直方向上的差分值。公式(6)无法直接求解，将公式(6)转化为多个子优化问题的求解，具体求解方法如下：

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

公式(9)的增广拉格朗日方程为：

由公式(10)得到两个子优化问题，分别为：

其中t表示迭代次数。接下来，对X_d、S_d和U_d进行交替更新、循环迭代即完成高光谱图像的重建。

步骤105：使用快速响应模型和快速差分模型更新高光谱图像X_d。

针对优化问题式(10)，低分辨率高光谱图像X_d的最小二乘解为：

由于矩阵H_d和矩阵D的规模很大，无法直接求其解析解，因此需利用共轭梯度下降法(详见Hestenes M R,Stiefel E.Methods of Conjugate Gradients for SolvingLinear Systems[J].Journal of Research of the National Bureau of Standards,1952,49(6):409-436.)求低分辨率高光谱图像X_d的近似解，从而完成低分辨率高光谱图像X_d的更新。

步骤105所述快速响应模型是为了加速执行PDCCHI系统的前向响应过程，也就是公式(3)所描述的PDCCHI获取采样图像的过程。使用查找表来实现快速响应模型，查找表的大小为2×(βN+βΩ-1)，其定义方式为：

使用公式(13)建立的查找表来重写CASSI前向模型即公式(1)，则为：

其中ω_d(λ)和Cu_d(i,j)分别为ω(λ)和Cu(i,j)经过下采样的结果。由于在使用GPU对前向模型实现并行加速时，公式(1)会出现边界判断从而大幅度降低并行加速的整体性能。而上述查找表优化方法能够降低这一影响，显著提高并行优化性能。

步骤105所述快速差分模型是为了优化前向差分过程即DX_d。使用共享内存对前向差分所用数据进行预加载，如图3所示，计算DX_d时数据直接从共享内存读取，从而完成快速差分模型。由于在计算DX_d时会频繁访问GPU的全局内存区，会降低并行数据的访存效率。利用共享内存的快速访问特点来对前向差分过程进行优化，能够提高前向差分的实施效率。

步骤106：使用快速差分模型更新低分辨率辅助矩阵S_d和U_d。

根据优化问题式(11)，低分辨率辅助矩阵S_d的最小二乘解为：

上式为软阈值收缩函数，运算符

表示点乘。使用步骤105所述快速差分模型，

能够直接求出。

辅助矩阵

的更新公式如下：

步骤107：更新辅助因子ρ和当前迭代次数t，并转至步骤105进行迭代，直至t＝I₁，完成低分辨率高光谱图像X_d的重建。

步骤108：对低分辨率高光谱图像X_d，低分辨率辅助矩阵S_d和U_d均使用快速插值模型进行双线性上采样插值，插值比例为1/β，插值结果分别为高分辨率高光谱图像X_u，高分辨率辅助矩阵S_u和U_u。

步骤108所述快速插值模型为步骤102所述的快速插值模型。

步骤109：初始化高分辨率高光谱图像X⁰＝X_u，高分辨率辅助矩阵S⁰＝S_u，U⁰＝U_u，辅助因子ρ，迭代次数t＝0。

步骤110：使用快速响应模型和快速差分模型更新高分辨率高光谱图像X，其更新方式为：

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

步骤111：使用快速差分模型更新高分辨率辅助矩阵S和U，更新方式为：

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

步骤112：更新辅助因子ρ和当前迭代次数t，并转至步骤110进行迭代，进行高分辨率高光谱图像X重建。直至t＝I₂，从而完成双相机光谱成像系统的实时重建，得到高质量的实时重建高光谱图像。

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

1.实验条件

本实验的硬件测试条件为：Inter i7 6800K，内存32G，Matlab 2015b。GPU为NVIDIX TITAN X，显存12G，CUDA 9.1。测试所用高光谱图片来自于ICVL数据集(详见AradB,Ben-Shahar O.Sparse Recovery of Hyperspectral Signal from Natural RGBImages[J].2016:19-34.)。编码孔径模版为p＝0.5的贝努利随机矩阵，采样因子β＝0.5。共轭梯度下降算法的迭代次数为25次，迭代停止门限值为1.0×10^-4。低分辨率重建阶段和正常分辨率重建阶段的迭代次数均为3，即I₁＝I₂＝3。采样过程中分别加入方差σ＝0和σ＝0.1的高斯白噪声。当σ＝0时，本发明公开的方法的τ＝0.1，ρ＝2.00；当σ＝0.1时，本发明公开的方法的τ＝0.25，ρ＝7.85。对比方法为基于全变差约束的TwIST算法，其最大迭代次数为100次

2.实验结果

为了验证本发明在高光谱重建效率上的提升，测试在不同分辨率下本发明公开的方法和对比方法的重建时间。

表1不同分辨率下平均时间对比(单位：秒)

从表1可以看出，本发明公开的方法都能够获得非常高的加速比。同时，随着分辨率的提升，本发明公开的方法的加速比也随之增加。以分辨率为512×512×30的重建图像为例，TwIST算法重建的时间长达7分钟，而使用本发明公开的方法在CPU上只需要1分钟，使用GPU仅仅只需要0.31秒，其加速效果非常显著。

图4为不同分辨率下，本发明公开的方法在GPU上能够获得的重建帧率。可以看出，在使用GPU并行加速的情况下，本发明公开的方法每秒能够完成12帧大小为256×256×14的高光谱图像的实时重建。

为了验证本发明的可行性，测试在不同噪声方差下，本发明公开的方法和对比方法的重建性能。为了定量的衡量重建结果的质量，使用峰值信噪比(Peaksignal to noiseratio，PSNR)衡量重建结果的空间质量和视觉效果，单位为dB；使用光谱角制图(Spectralangle mapping，SAM)(详见Kruse F A,Lefkoff A B,Boardman J W,et al.The spectralimage processing system(SIPS)—interactive visualization and analysis ofimaging spectrometer data[J].Remote sensing of environment,1993,44(2-3):145-163.)衡量重建结果的光谱保真度。重建结果如下所示。

表2噪声方差σ＝0时的重建结果

表3噪声方差σ＝0.1时的重建结果

从表2和表3的结果可以看出，本发明公开的方法能够在迭代很少的次数下达到非常好的重建质量。无论是空间质量还是光谱保真度，在不同强度的噪声影响下，本发明公开的方法都能够与TwIST算法相近甚至更好的重建结果。

图5-a、图5-b和图5-c分别为测试图像eve-1551在610nm时的参考图、TwIST重建结果和本发明公开算法的重建结果。图6-a、图6-b和图6-c分别为测试图像lst-0950在610nm时的参考图、TwIST重建结果和本发明公开算法的重建结果。可以看出，本发明公开的方法重建的结果在视觉上的效果与TwIST算法相似甚至更好。

综上所述，本发明公开的方法能够获得非常好的空间重建质量和光谱保真度。在不损失重建精度的情况下，本发明公开的方法在GPU并行加速的情况下能够达到每秒12帧256×256×14高光谱图像的实时重建。因此本发明公开的方法具有非常高的实际利用价值和指导意义。

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

Claims (7)

1.一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：包括以下步骤，

步骤101：输入双相机光谱成像系统的采样图像Y、标定后的前向响应矩阵H、正则化系数τ，下采样比例β，低分辨率重建迭代次数I₁，高分辨率重建迭代次数I₂；

步骤102：使用快速插值模型对采样图像Y和前向响应矩阵H进行下采样，得到低分辨率采样Y_d和低分辨率前向响应矩阵H_d，空间和光谱的下采样因子均为β；

步骤102所述快速插值模型为统一化插值算法；线性图像插值过程都能够看作是源图像像素的线性叠加过程，故用统一的数学模型表示为：

其中d和s分别表示目标图像和源图像，j表示像素点索引；K表示插值所需要的源像素点个数；I(j,m)和W(j,m)表示索引矩阵和权值矩阵；一旦图像分辨率和采样比例确定以后，I(j,m)和W(j,m)就能够提前生成并作为查找表存储在内存中，从而完成快速图像插值模型；

步骤103：初始化低分辨率高光谱图像

以及低分辨率辅助矩阵

和

步骤104：初始化优化目标函数的初始值

辅助因子ρ；初始化当前迭代次数t＝0；

步骤105：使用快速响应模型和快速差分模型更新低分辨率高光谱图像X_d；

步骤105具体实现方法如下，

针对优化问题式(10)，低分辨率高光谱图像X_d的最小二乘解为：

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

步骤105所述快速响应模型是为了加速执行PDCCHI系统的前向响应过程，也就是公式(3)所描述的PDCCHI获取采样图像的过程；使用查找表来实现快速响应模型，查找表的大小为2×(βN+βΩ-1)，其定义方式为：

使用公式(13)建立的查找表来重写CASSI前向模型即公式(1)，则为：

其中ω_d(λ)和Cu_d(i,j)分别为ω(λ)和Cu(i,j)经过下采样的结果；由于在使用GPU对前向模型实现并行加速时，公式(1)会出现边界判断从而大幅度降低并行加速的整体性能；而上述查找表优化方法能够降低这一影响，显著提高并行优化性能；

步骤105所述快速差分模型是为了优化前向差分过程即DX_d；使用共享内存对前向差分所用数据进行预加载，计算DX_d时数据直接从共享内存读取，从而完成快速差分模型；由于在计算DX_d时会频繁访问GPU的全局内存区，会降低并行数据的访存效率；利用共享内存的快速访问特点来对前向差分过程进行优化，能够提高前向差分的实施效率；

步骤106：使用快速差分模型更新低分辨率辅助矩阵S_d和U_d；

步骤106具体实现方法如下，

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

上式为软阈值收缩函数，运算符

表示点乘；使用步骤105所述快速差分模型，

能够直接求出；

辅助矩阵

的更新公式如下：

步骤107：更新辅助因子ρ和当前迭代次数t，并转至步骤105进行迭代，直至t＝I₁，完成低分辨率高光谱图像X_d的重建；

步骤108：对X_d，S_d和U_d均使用快速插值模型进行双线性上采样插值，插值比例为1β，插值结果分别为高分辨率高光谱图像X_u,高分辨率辅助矩阵S_u和U_u；

步骤108所述快速插值模型为步骤102所述的快速插值模型；

步骤109：初始化高分辨率高光谱图像X⁰＝X_u，高分辨率辅助矩阵S⁰＝S_u，U⁰＝U_u，辅助因子ρ，迭代次数t＝0；

步骤110：使用快速响应模型和快速差分模型更新高分辨率高光谱图像X；

步骤111：使用快速差分模型更新高分辨率辅助矩阵S和U；

步骤112：更新辅助因子ρ和当前迭代次数t，并转至步骤110进行迭代，进行高分辨率高光谱图像X的重建，直至t＝I₂，从而完成双相机光谱成像系统的实时重建，得到高质量的实时高光谱重建图像。

2.如权利要求1所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤101中双相机光谱成像仪为基于全色相机的双相机光谱成像系统；双相机光谱成像系统包括分光镜、物镜、编码模版、中继镜、色散棱镜和全色相机部件，包含一个编码孔径快照光谱成像系统和一个全色相机分支；目标场景的高光谱图像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)为二维混叠采样图像；

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

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

Y＝HX (3)

其中Y＝[Y^c；Y^p]表示二维压缩光谱采样图像+二维灰度图像；X表示三维数据立方体，H表示对编码孔径快照光谱成像系统CASSI分支+全色相机分支的前向响应矩阵，为编码模版函数、色散棱镜偏移函数和CCD相机频谱响应的综合作用。

3.如权利要求1所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤102所述对采样图像Y和前向响应矩阵H进行下采样的方法如下：对于采样图像Y，使用最近邻插值对其CASSI混叠采样部分下采样，使用均值滤波对其全色采样部分进行下采样；同时，使用均值滤波对前向响应矩阵进行下采样。

4.如权利要求3所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤103所述辅助矩阵

和

为重建低分辨率高光谱图像时需要用到的辅助矩阵，大小均为2×βM×βN×βΩ，初始化为全零矩阵。

5.如权利要求4所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤103重建低分辨率高光谱图像

的初始化方式为：

其中

表示下采样得到的编码孔径快照光谱成像系统的前向响应H_d的转置，

表示由下采样二维观测反演到低分辨率高光谱图像；由于空间和光谱均进行了比例为β的下采样，因此低分辨率高光谱图像

的大小为βM×βN×βΩ。

6.如权利要求5所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤104所述优化目标函数为全局优化目标函数；根据单谱带图像的分段平滑特性，将高光谱重建问题转化为基于全变差约束的最优化问题，从而得到全局优化目标函数为：

其中运算符

表示矩阵Frobenius范数的平方；||DX_d||₁表示图像的全变差值，其定义如下：

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

公式(7)的结果为两个大小和X_d相同的矩阵，分别表示图像X_d在水平方向和竖直方向上的差分值；公式(6)无法直接求解，将公式(6)转化为多个子优化问题的求解，具体求解方法如下：

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

公式(9)的增广拉格朗日方程为：

由公式(10)得到两个子优化问题，分别为：

其中t表示迭代次数；接下来，对X_d、S_d和U_d进行交替更新、循环迭代即完成高光谱图像的重建。

7.如权利要求6所述的一种基于插值的双相机光谱成像系统的实时重建方法，其特征在于：步骤110使用快速响应模型和快速差分模型更新高分辨率高光谱图像X，其更新方式为：

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

步骤111使用快速差分模型更新高分辨率辅助矩阵S和U，更新方式为：

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

Images