CN112556848A - 双相机压缩测量协同张量表示的融合计算成像方法 - Google Patents

Abstract

本发明公开了一种双相机压缩测量协同张量表示的融合计算成像方法，包括：建立双光路压缩测量装置，一路为编码孔径快照光谱测量装置，另一路为全色相机测量装置；通过伯努利随机矩阵生成的压缩编码孔径、色散棱镜引起的偏移、相机的光谱响应构造前向压缩采样过程；两路测量信息经过三阶张量表示，综合各向异性加权空‑谱全变差先验约束，建立双相机压缩测量协同张量表示正则化模型；采取交替方向乘子方法进行模型迭代优化求解。本发明方法中的张量表示能充分地挖掘高光谱图像的立方体图谱特性，并能有效地利用全色图像提供的空间细节信息，通过融合计算，大幅度地提升高光谱成像质量和空‑谱特征保持能力。

Description

双相机压缩测量协同张量表示的融合计算成像方法

技术领域

本发明属于计算摄像学领域，特别是一种双相机压缩测量协同张量表示的融合计算成像方法。

背景技术

高光谱图像是由几十甚至上百连续光谱波段组成的三维数据集，相比于RGB彩色图像，高光谱图像具有丰富的光谱信息和细节结构。在环境和军事监测、农业规划、矿物勘察和医学检测等领域具有广泛的应用。

传统的高光谱成像方法，如推扫型高光谱成像仪和摆扫式高光谱成像仪，一次只能捕获一个一维或二维的光谱立方体子集，需要采用时序的扫描策略才能捕获一个完整的三维光谱数据。因而这种成像方法是耗时的，且不能捕获动态的场景。近几年，随着计算成像技术和压缩感知理论的发展，压缩成像测量方法以远低于奈奎斯特采样数完成对目标场景的压缩采样，并利用相关计算成像算法完成对高光谱图像的复原。最典型的压缩成像测量装置是由美国杜克大学DavidBrady等人研制的编码孔径快照光谱成像测量装置(M.Gehm,R.John,D.Brady,R.Willett,and T.Schulz,“Single-shot compressivespectral imaging with a dual-disperser architecture,”Opt.Exp.,vol.15,no.21,pp.14013–14027,Oct.2007)。为了提升重构的质量，Wang等人提出双相机压缩成像测量装置(L.Wang,Z.Xiong,D.Gao,G.Shi,and F.Wu,“Dual-camera design for coded aperturesnapshot spectral imaging,”Applied Optics,vol.54,no.4,pp.848–58,2015.)，它将编码孔径快照光谱成像测量装置和全色相机结合，实现了光谱信息和空间信息的互补。

从这两路采样中融合计算成像原始的高光谱图像是一个欠定的难题。这个压缩感知反问题通常被转化为求解一个l₁范数最优化问题。通过增加稀疏约束，假定高光谱图像在小波基下是可稀疏表示的，该问题可通过基于稀疏重构的梯度投影算法(gradientprojection for sparse reconstruction,GPSR)(Figueiredo M AT,Nowak R D,Wright SJ.Gradient Projection for Sparse Reconstruction:Application to CompressedSensing and Other Inverse Problems[J].IEEE Journal of Selected Topics inSignal Processing,2008,1(4):586-597.)求解。以及基于分段平滑的全变差约束法(2Dtotal variation,2DTV)(Chambolle A.An Algorithm for Total VariationMinimization and Applications[J].Journal of Mathematical Imaging&Vision,2004,20(1-2):89-97.)。但这类算法很容易导致重构的图像出现过平滑的现象，丢失了很多纹理信息。Peng等人提出一种E-3DTV模型(孟德宇，彭江军，谢琦等.一种基于E-3DTV正则的高光谱图像修复方法.CN109191404A.2019)，是对传统3DTV的增强和改进，对三个方向的梯度图做低秩矩阵分解并对基矩阵加稀疏约束，得到E-3DTV正则，更能降低噪声的干扰。Fu等人提出一种低秩矩阵约束模型(low-rank matrix approximation，LRMA)(Fu Y,Zheng Y,SatoI,et al.Exploiting Spectral-Spatial Correlation for Coded Hyperspectral ImageRestoration[C]//Computer Vision&Pattern Recognition.IEEE,2016.)，能够有效地利用高光谱图像的空间-光谱相关性，计算成像质量有一定提升。但这些方法都基于矩阵-向量形式的观测模型，把三维高光谱图像展开成一个一维向量，破坏了高光谱图像原有的内在属性，导致重构的质量有限，从而限制了实际的应用。王立志等人提出一种维度可区分的张量低秩约束模型(dimension-discriminative low-rank tensor regularization)(王立志，张仕鹏，黄华.一种张量低秩约束的快照式光谱成像系统的重构方法，CN110501072A.2019)，利用了高光谱图像的高维结构特性，并能表达高光谱图像在多个维度的相关性，有利于提升高光谱成像的精度。这些方法都把两路观测模型转化为一项压缩测量保真项，未能提升压缩测量项与全色图像两者之间的自适应，且未能完全有效地利用全色图像的空间信息。

发明内容

本发明的目的在于提供一种双相机压缩测量协同张量表示的融合计算成像方法，该方法中的张量表示能充分地挖掘高光谱图像的立方体图谱特性，将两路观测模型转化为压缩测量光谱保真项和全色约束保真项，通过协同张量Tucker分解与全色图像建立耦合关系，综合各向异性加权空-谱全变差，通过融合计算，大幅度地提升高光谱成像质量和空-谱特征保持能力。

实现本发明目的的技术解决方案为：一种双相机压缩测量协同张量表示的融合计算成像方法，包括如下步骤：

步骤1，建立双光路压缩测量装置，形成空-谱信息压缩测量；

步骤2，构造前向压缩采样算子，描述两路压缩测量的观测过程；

步骤3，建立双相机压缩测量协同张量表示正则化模型；

步骤4，交替方向乘子方法迭代优化求解，输出融合计算成像结果。

本发明与现有技术相比，其显著特点在于：(1)本发明建立了一种新的张量表示的双相机压缩测量观测模型，为后续更好地挖掘高光谱图像的立方体图谱特性；(2)将两路观测模型转化为压缩测量光谱保真项和全色约束保真项之和，提升了这两项的自适应；(3)通过协同张量Tucker分解，更能充分地利用高光谱图像波段之间的全局空间和光谱相关性；(4)在全色约束保真项中，将高光谱图像协同张量Tucker分解中的核心张量和空间因子矩阵与全色图像分解中的核心张量和空间因子矩阵建立耦合关系，从而能更好地利用全色图像的空间结构信息；(5)引入各向异性加权空-谱全变差，能够有效地保持高光谱图像全局空间分段平滑性和光谱连续性特性；(6)本发明方法采用交替方向乘子法交替迭代优化求解，有较好的收敛性和鲁棒性。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明双相机压缩测量协同张量表示的融合计算成像方法的流程图。

图2是本发明中交替方向乘子方法求解的流程图。

图3是本发明中用于双相机压缩成像测量装置的结构流程图。

图4(a)是本发明在双相机压缩成像测量装置下，探测器的光谱响应函数曲线，图4(b)为p＝0.5的伯努利矩阵生成的编码孔径。

图5(a)和图5(b)是本发明在双相机压缩成像测量装置下，分别对测试图像flowers仿真得到的压缩测量值和全色图像。

图6(a)～图6(e)是本发明和对比方法在双相机压缩成像测量装置下，对测试图像flowers仿真计算成像后在波长为600nm时的结果图及局部放大图，其中图6(a)为参考图像，图6(b)、图6(c)、图6(d)和图6(e)分别为2DTV算法、GPSR算法、E-3DTV算法和本发明的计算成像结果。

图7(a)和图7(b)是本发明在双相机压缩成像测量装置下，分别对测试图像toy仿真得到的压缩测量值和全色图像。

图8(a)～图8(e)是本发明和对比方法在双相机压缩成像测量装置下，对测试图像toy仿真计算成像后在波长为600nm时的结果图及局部放大图，其中图8(a)为参考图像，图8(b)、图8(c)、图8(d)和图8(e)分别为2DTV算法、GPSR算法、E-3DTV算法和本发明的计算成像结果。

具体实施方式

结合图1，下面详细说明本发明的实施过程，步骤如下：

步骤1：建立双光路压缩测量装置，形成空-谱信息压缩测量；

如图3所示，所述的双光路压缩测量装置，其中一路透过分光镜进入编码孔径快照光谱成像测量装置，通过建立编码孔径快照光谱成像的数据采集测量光路，该光路包括依次设置的物镜、编码孔径、带通滤波器、中继镜、色散棱镜和探测器，其中使用的编码孔径模板如图4(b)所示；另一路经分光镜反射进入全色相机光路，包括带通滤波器和探测器。在双相机压缩成像测量装置下，进行仿真实验时，假设两路光路的相机具有相同的如图4(a)所示光谱响应曲线，图5(a)和图5(b)是对原始高光谱图像flowers仿真得到的压缩测量值和全色图像；图7(a)和图7(b)是对原始高光谱图像toy仿真得到的压缩测量值和全色图像。

步骤2：构造前向压缩采样算子，描述两路压缩测量的观测过程；

前向压缩采样算子是针对编码孔径快照光谱成像测量装置的压缩采样过程，进入编码孔径快照光谱成像测量装置的光源，先经过随机伯努利矩阵的编码孔径的编码；通过带通滤波器进行滤波，限制光谱的范围；编码和滤波后的图像经由中继镜传播到色散棱镜后，不同波段的图像会沿着竖直方向作不同程度的偏移；最后，所有波段的图像在探测器上进行叠加，得到压缩采样的二维图像。将目标场景的高光谱图像表示为一个三维张量

大小为M×N×L，高光谱图像

上任意一点的像素值为x(i,j,k),1≤i≤M，1≤j≤N，1≤k≤L，其中i和j表示空间坐标大小，k表示光谱坐标大小。其数学模型可以表示为：

其中，ω_c(k)表示该探测器的光谱响应函数，T(i,j)表示编码孔径的传播函数，S(k)表示由色散棱镜引起的偏移函数，通常假定为线性色散。用Φ来表示编码孔径快照光谱成像测量装置的前向压缩采样算子，因而它包含编码孔径、色散棱镜和探测器光谱响应函数的共同作用。建立张量表示的观测模型为：

其中g为该测量装置的测量值，大小为M×(N+L-1)，v₁表示高斯白噪声。

另一路全色相机光路，直接到达相机，其数学模型表示为：

其中，ω_p(k)表示全色相机的光谱响应函数。

同样，建立张量表示的全色图像的观测模型为：

其中，全色图像Y是高光谱图像的光谱下采样形式，大小为M×N，R为全色相机的光谱响应函数，大小为1×L，v₂为高斯白噪声。

步骤3：建立双相机压缩测量协同张量表示正则化模型；

所述双相机压缩测量协同张量表示正则化模型包含四项：压缩测量光谱保真项、全色约束保真项、核心张量能量约束项和各向异性加权空-谱全变差，其计算模型为：

其中，运算符

表示Frobenius范数，λ为正权重参数，用于提升压缩测量光谱保真项和全色约束保真项的自适应，β>0，τ>0为正则项参数，

为张量

的Tucker分解中的核心张量

与空间和光谱方向的因子矩阵

且满足C_m<M，C_n<L和C_l<L，

为光谱维因子矩阵的下采样，通过增加核心张量能量约束项，记该张量Tucker分解为协同张量Tucker分解，各向异性加权空-谱全变差(spatial-spectral total variation，SSTV)表示为

其中x_i,j,k为高光谱图像

中的第(i,j,k)元素，ω₁，ω₂，ω₃为三个权重参数，本实施例中设置ω₁＝ω₂＝1，仅改变ω₃的值。

步骤4：交替方向乘子方法迭代优化求解，输出融合计算成像结果；

引入辅助变量

和

将协同张量正则化模型转化为以下等价形式：

其中，D＝[ω₁D_m；ω₂D_n；ω₃D_l]表示所谓的加权三维差分算子，其中D_m，D_n和D_l分别是高光谱图像在空间和光谱方向的一阶差分算子；

上述带约束的最优化问题能被交替方向乘子方法有效求解，其对应的增广拉格朗日函数为：

其中

为三个拉格朗日乘子，ρ>0，γ>0为惩罚参数。

交替方向乘子方法通过迭代优化求解增广拉格朗日函数极小值，输出融合计算成像结果，流程图如图2所示。

具体的交替迭代优化过程如下：

设置迭代的初始值t＝0，

以及随机初始化M⁰，N⁰，L⁰和通过求解目标函数(20)来初始化化

迭代停止条件为：t>t_max或

1)更新

固定其它变量，提取含

的项，得到

对(8)求导，将含

的式子移到等式一边，可得以下线性方程组

其中，Φ^*是Φ的伴随算子，I为单位矩阵。

该等式可以由预处理共轭梯度法求解。

2)更新

固定其它变量，提取含

的项，得到

可得如下线性等式：

该等式有如下闭式解：

其中D^*是D的伴随算子，

表示三维傅里叶变换，

是其逆变换。

3)更新

固定其它变量，提取含

的项，得到

可得如下闭式解：

其中软阈值收缩算子soft(x,y)＝sign(x)*max(|x|-y,0)。

4)更新M，提出含M的项，得到如下子问题：

根据张量模-n展开的性质：

能得到下列等同的目标函数：

其中，Y₍₁₎,X₍₁₎和U₍₁₎分别是张量

和

的模-1展开，

对目标函数求导，取零，M有如下闭式解：

5)更新N，与求解M类似，可得如下闭式解：

其中，Y₍₂₎,X₍₂₎和U₍₂₎分别是张量

和

的模-2展开，

6)更新L，与求解M和N类似，可到如下广义Sylvester方程：

其中，Y₍₃₎,X₍₃₎和U₍₃₎分别是张量

和

的模-3展开，

可采用共轭梯度法求解该等式。

7)更新

固定其它变量，可得关于

的目标函数：

对于

采用ADMM方法求解，引入变量

可得如下约束形式：

其对应增广拉格朗日函数为：

其中v为拉格朗日乘子，μ是惩罚因子，

y，c，x，j和u₁分别是Y，

和

的向量化形式。

循环更新c，j和v如下:

对于变量c，提取其目标函数为：

其闭式解为：

由于矩阵

太大，直接计算它的逆矩阵太困难，可得

其中P_i，Λ_i(i＝1,2,3)分别由M^TM，N^TN和(RL)^T(RL)的特征值分解得到，且P_i为酉矩阵，Λ_i为对角矩阵。所以

为对角矩阵，很容易计算它的逆矩阵。

对于变量j，可得它的闭式解为

采取和求c相同的策略求解j。

更新乘子v,v＝v-(v-j)

8)更新拉格朗日乘子和参数：

ρ＝ξ₁×ρ,γ＝ξ₂×γ (28)

其中，ξ₁>0，ξ₂>0。

本发明的效果可通过以下仿真实验进一步说明：

(1)仿真条件

仿真实验采用的高光谱数据为公开的CAVE数据集，空间分辨率为512×512，共有31个波段，其中波长范围是400-700nm，间隔为10nm。在我们的实验中，裁剪图像大小为256×256×31，且将图像的强度值归一化到[0,1]。编码孔径快照光谱成像测量装置中所用的编码孔径是概率p＝0.5的随机伯努利矩阵，色散棱镜的色散为线性等距色散。仿真得到大小为256×286的压缩测量值和大小为256×256的全色图像，均添加方差为0.1的高斯白噪声。设置最大迭代次数t_max＝100和迭代收敛阈值∈＝1×10^-6。所有测试图像的参数都统一设置为：C_m＝C_n＝70,C_l＝3,权重参数λ＝2，β＝10^-3,ρ＝0.1，γ＝10^-4，ω＝[1,1,0.1],ξ₁＝1.1,ξ₂＝1.2。对比方法为2DTV算法、GPSR算法和E-3DTV算法。仿真实验均在Windows 10操作系统下采用MATLAB R2016a完成，硬件测试条件：Inter i7 6700,内存8G。

为了衡量高光谱成像的质量，本发明采用的评价指标是峰值性噪比(peak singalto noise ratio，PSNR)、均方根误差(root mean square error，RMSE)和相对无量纲全局误差(relative dimensionless global error in synthesis，ERGAS)，其中PSNR值越大，RMSE值越小，ERGAS值越小，重构质量越好。计算成像结果图像和参考图像均映射到[0,255]进行指标的计算。

(2)仿真实验结果分析

针对CAVE数据集下7幅图像的实验结果如表1。

表1不同方法的计算成像质量指标对比

通过表1可以看出，2DTV算法、GPSR算法和E-3DTV算法具有比较接近的结果，这些方法都需要把三维高光谱图像展开成一维向量进行求解，破坏了高光谱图像原有的三维结构特性，因而降低了计算成像的质量。而本发明公开的双相机压缩测量协同张量表示的融合计算成像方法获得了最好的成像结果，说明了本发明能够大幅度地提升高光谱成像的质量和空-谱特征保持能力。

图6(a)～图6(e)和图8(a)～图8(e)分别是测试图像flowers和toy在双相机压缩成像测量装置下的600nm处的计算成像结果及局部放大图，从结果图可以看出基于二维全变差的2DTV算法，很容易导致图像的过平滑，丢失许多细节信息。E-3DTV算法相比2DTV和GPSR算法，虽然能较好地保留纹理信息，但有时会导致图像过于锐化，与参考图像差异较大。而本发明公开方法相比这几种方法，成像结果是最好的，能够有效地复原和保留图像的空间细节和纹理信息，再一次说明了本发明公开方法的有效性。

Claims (6)

1.一种双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，包括以下步骤：

步骤1，建立双光路压缩测量装置，形成空-谱信息压缩测量；

步骤2，构造前向压缩采样算子，描述两路压缩测量的观测过程；

步骤3，建立双相机压缩测量协同张量表示正则化模型；

步骤4，交替方向乘子方法迭代优化求解，输出融合计算成像结果。

2.根据权利要求1所述的双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，步骤1中所述的双光路压缩测量装置，其中一路透过分光镜进入编码孔径快照光谱成像测量装置，通过建立编码孔径快照光谱成像的数据采集测量光路，该光路包括依次设置的物镜、编码孔径、带通滤波器、中继镜、色散棱镜和探测器；另一路经分光镜反射进入全色相机光路，全色相机光路包括带通滤波器和探测器。

3.根据权利要求2所述的双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，步骤2中所述前向压缩采样算子是针对编码孔径快照光谱成像测量装置的压缩采样过程，进入编码孔径快照光谱成像测量装置的光源，先经过随机伯努利矩阵的编码孔径的编码；通过带通滤波器进行滤波，限制光谱的范围；编码和滤波后的图像经由中继镜传播到色散棱镜后，不同波段的图像会沿着竖直方向作不同程度的偏移；最后，所有波段的图像在探测器上进行叠加，得到压缩采样的二维图像；

将目标场景的高光谱图像表示为一个三维张量

大小为M×N×L，高光谱图像

上任意一点的像素值为x(i,j,k),1≤i≤M，1≤j≤N，1≤k≤L，其中i和j表示空间坐标大小，k表示光谱坐标大小；其数学模型可表示为：

其中，ω_c(k)表示该探测器的光谱响应函数，T(i,j)表示编码孔径的传播函数，S(k)表示由色散棱镜引起的偏移函数；用Φ表示编码孔径快照光谱成像测量装置的前向压缩采样算子，它包含编码孔径、色散棱镜和探测器光谱响应函数的共同作用；建立张量表示的观测模型为：

其中g为该测量装置的测量值，大小为M×(N+L-1)，v₁表示高斯白噪声；

另一路全色相机光路，直接到达相机，其数学模型表示为：

其中，ω_p(k)表示全色相机的光谱响应函数；

同样，建立张量表示的全色图像的观测模型为：

其中，全色图像Y是高光谱图像

的光谱下采样形式，大小为M×N，R为全色相机的光谱响应函数，大小为1×L，v₂为高斯白噪声。

4.根据权利要求1所述的双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，步骤3中所述双相机压缩测量协同张量表示正则化模型包含四项：压缩测量光谱保真项、全色约束保真项、核心张量能量约束项和各向异性加权空-谱全变差，其计算模型为：

其中，运算符

表示Frobenius范数，λ为正权重参数，用于提升压缩测量光谱保真项和全色约束保真项的自适应，β>0，τ>0为正则项参数，

为张量

的Tucker分解中的核心张量

与空间和光谱方向的因子矩阵

且满足C_m<M，C_n<L和C_l<L，

为光谱维因子矩阵的下采样，通过增加核心张量能量约束项，记该张量Tucker分解为协同张量Tucker分解，各向异性加权空-谱全变差表示为：

其中x_i,j,k为高光谱图像

中的第(i,j,k)元素，ω₁，ω₂，ω₃为三个权重参数。

5.根据权利要求4所述的双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，设置ω₁＝ω₂＝1，仅改变ω₃的值。

6.根据权利要求1所述的双相机压缩测量协同张量表示的融合计算成像方法，其特征在于，步骤4中引入辅助变量

和

将协同张量正则化模型转化为以下等价形式：

其中，D＝[ω₁D_m；ω₂D_n；ω₃D_l]表示所谓的加权三维差分算子，其中D_n，D_n和D_l分别是高光谱图像在空间和光谱方向的一阶差分算子；

上述带约束的最优化问题能被交替方向乘子方法有效求解，其对应的增广拉格朗日函数为：

其中

为三个拉格朗日乘子，i＝1,2,3；ρ>0，γ>0为惩罚参数；

交替方向乘子方法通过迭代优化求解增广拉格朗日函数极小值，输出融合计算成像结果。

Images