[发明专利]一种编码孔径光谱成像系统的重构方法

权利要求书

1.一种编码孔径光谱成像系统的重构方法，其特征在于：包含以下步骤，

步骤101：输入编码孔径快照光谱成像仪(Coded Aperture Snapshot SpectralImager，CASSI)获取的压缩光谱采样图像Y、标定后的前向响应矩阵H、主函数最大迭代次数Maxiters和正则化系数τ；

步骤102：初始化重构高光谱图像X⁰，优化目标函数的初始值f⁰，以及交替方向乘子法(Alternating Direction Method of Multipliers，ADMM)辅助矩阵S和U，辅助因子ρ；初始化当前迭代次数t＝0；

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

X⁰＝H^TY (1)

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

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

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

||DX||₁＝∑_i，j，λ|x(i，j+1，λ)-x(i，j，λ)|+|x(i+1，j，λ)-x(i，j，λ)| (3)

x(i，j，λ)为高光谱图像X任意一点的像素值；矩阵D为差分矩阵，矩阵D和图像X的作用结果如下：

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

步骤102所述辅助矩阵S和U是指使用ADMM算法求解优化方程(2)时所需要的变量，大小均为2×M×N×Ω，且都初始化为全零矩阵；M表示高光谱图像空间维的长度；N表示高光谱图像空间维的宽度；M×N表示高光谱图像的空间分辨率；Ω表示高光谱图像的频谱分辨率；根据ADMM算法，对公式(2)引入辅助矩阵S＝DX，得到如下带约束的优化方程：

增广拉格朗日方程为：

其中ρ为ADMM辅助因子；

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

步骤103具体实现方法为，

首先更新高光谱图像X；固定矩阵U和S不变，求解使得公式(6)最小时高光谱图像X的值；优化目标为：

其最小二乘解为：

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

利用共轭梯度下降法求高光谱图像X^t+1的近似解，从而完成高光谱图像X^t+1的更新；

步骤104：更新ADMM算法的辅助矩阵S和U；

步骤104具体实现方法为，

辅助矩阵S的更新是在保持X和U不变的情况下，求解使得公式(6)最小时矩阵S的值；其优化目标为：

其最小二乘解，也就是矩阵S的更新公式为：

上式为软阈值收缩函数，S^t+1能够直接求出；

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

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

U^t+1能够直接求得；

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

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

步骤106：根据步骤105计算的结果执行迭代选择策略，完成高光谱图像的重构；

步骤106所述迭代选择策略如下，

如果收敛，即步骤105的目标函数值f^t+1小于上一个目标函数值f^t，则计算目标函数值的相对变化量：

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

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

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

其中ω(λ)表示CCD相机的频谱响应函数；Cu(i，j)表示编码模版函数；φ(λ)表示色散棱镜的波段偏移函数；y(i，j)为压缩光谱采样图像，v(i，j)表示高斯白噪声；式(14)写成矩阵的形式为：

Y＝HX+V (15)

其中Y表示二维压缩光谱采样图像，大小为M×(N+Ω-1)；X表示三维数据立方体，V表示高斯白噪声；H表示对CASSI标定后的前向响应矩阵，为编码模版函数、CCD相机频谱响应和色散棱镜偏移函数的综合作用，大小为M×(N+Ω-1)×L。