CN104050653A - 基于非负结构稀疏的高光谱图像超分辨率算法 - Google Patents

Abstract

本发明公开了一种基于矩阵结构稀疏非负分解的光谱图像超分辨率重建算法，从低空间分辨率高光谱图像和高分辨率彩色图像联合重建高分辨率高光谱图像，解决现有算法难以精确恢复高分辨率光谱图像的问题。其实现步骤是：(1)输入低分辨率高光谱图像和对应的高分辨率彩色图像；(2)利用光谱图像的局部与非局部自相似性，构建基于矩阵结构稀疏非负分解的光谱重建目标函数；(3)采用交替方向乘子法交替求解得到优化的光谱材料系数和光谱材料基；(4)利用得到的优化光谱材料系数矩阵和光谱材料基矩阵，重建出高分辨率高光谱图像。本发明所恢复的高光谱图像更清晰，图像边缘更锐利，可有效提升光谱图像空间分辨率。

Description

基于非负结构稀疏的高光谱图像超分辨率算法

技术领域

本发明涉及图像处理技术领域，具体涉及一种高分辨率的高光谱图像重建方法，主要用于从低分辨率高光谱图像和高分辨率彩色图像联合重建高分辨率的高光谱图像。

背景技术

高光谱图像在国防科学研究和军事应用领域均有着重要需求。受传统高空间分辨率成像内在固有成像特性和高分辨率传感器限制，现有高光谱成像方法在实现高光谱成像时很难达到高分辨率。所以使用信号处理技术从低分辨率高光谱图像重建出高分辨率高光谱图像成为高分辨率高光谱图像获取的一个重要途径。

为了获得高分辨率的高光谱图像，Kawakami等人在“High-resolutionHyperspectral Imaging via Matrix Factorization，2011”中提出了用低分辨率光谱相机和高分辨率彩色相机获取高分辨率光谱图像的方法，该方法采用信号处理的方法从低分辨率高光谱图像和高分辨率彩色图像重建出高分辨率高光谱图像。由于该光谱重建问题是一个病态逆问题，为了获得高精度光谱图像，Kawakami等人提出了基于光谱稀疏性的重建方法，利用光谱图像在光谱材料基上的稀疏分解特性，联合彩色图像重建高空间分辨率高光谱图像。针对Kawakami等人提出的光谱重建方法中，光谱材料基和光谱材料系数可能为负数而不满足实际物理含义的问题，Eliot Wycoff等人在“A Non-Negative Sparse Promoting Algorithm For High ResolutionHyperspectral Imaging，2013”中提出了一种基于矩阵非负分解(Sparse Non-NegativeMatrix Factorization,SNNMF)的光谱图像重建方法，利用光谱图像和彩色图像的相关性，构建了基于矩阵非负分解的光谱稀疏重建模型，通过交替方向乘子法进行高效优化求解，重建出高空间分辨率高光谱图像。

另一方面，现有的基于SNNMF的光谱重建方法只利用了高光谱图像谱间稀疏性，没有充分利用光谱图像存在的较强的局部和非局部的空间结构相似性，导致难以精确恢复精细的光谱图像的边缘和纹理细节。为了进一步提升光谱图像重建的稳定性和精度，需要进一步挖掘光谱图像的空间结构相关性。

发明内容

本发明的目的在于针对现有基于矩阵非负分解的光谱超分辨率重建算法没有利用光谱图像局部和非局部的相似性的问题，提出一种基于矩阵结构稀疏非负分解的光谱图像超分辨重建算法，从低分辨率光谱图像和高分辨率彩色图像重建出高精度高分辨率光谱图像。通过利用光谱图像的局部和非局部结构相似性，采用矩阵结构范数构建更加稳健、精确的光谱图像矩阵稀疏非负分解模型，有效提升光谱图像的重建精度。

本发明的技术思路是：对光谱图像谱线进行结构聚类；对每一个类中的相似谱线在材料基上采用结构范数同时进行非负分解，使得各相似谱线的分解系数具有结构相似性；构建基于矩阵结构稀疏非负分解的光谱重建优化目标函数，采用交替方向乘子法交替求解光谱材料系数和光谱材料基，重建出高分辨率高光谱图像。具体步骤包括如下：

(1)分别输入低分辨率高光谱图像和高分辨率的彩色图像其中，M_h表示谱段数，L_h表示低分辨率高光谱图像每个谱段的像素数，L_c表示高分辨率高光谱图像每个谱段的像素数，且L_h＜L_c，M_c为彩色图像的谱段数，M_c＝3；

(2)假设高分辨率光谱图像Z可以表示为：其中表示光谱材料基矩阵，表示光谱材料系数矩阵，其每一列表示各谱线在A上的稀疏分解系数。假设低分辨率光谱图像X和彩色图像Y与高分辨率光谱图像Z之间满足如下线性关系X＝ASG+W，Y＝FAS+V，其中表示下采样矩阵，表示彩色图像变换矩阵，和表示加性高斯白噪声。利用高光谱图像的局部与非局部相似性，构建如下基于结构稀疏光谱非负分解目标函数：

$(A,S)=\underset{A\∈R_{+}^{M_h\×N}}{\underset{S\∈R_{+}^{N\×L_C}}{\arg \min }}\frac{1}{2}{\left|\right|X-\mathrm{ASG}\left|\right|}_F^2+\frac{1}{2}{\left|\right|Y-\mathrm{FAS}\left|\right|}_F^2+\mathrm{\λ}\underset{k}{\mathrm{\Σ}}{\left|\right|{\mathrm{SR}}_k\left|\right|}_{\mathrm{1,2}},$

其中R_k表示相似谱线索引矩阵，SR_k表示从稀疏系数矩阵S中取出第k类中相似谱线对应的稀疏分解系数向量，即表示第k类中第p条相似谱线对应的稀疏系数，p＝1,2...n，n为第k类中相似谱线对应稀疏系数的个数，表示矩阵的弗罗贝尼乌斯范数的平方，λ表示正则参数，||·||_1,2表示对矩阵每行求二范数再求和，即w为SR_k矩阵的行数，r为SR_k矩阵的第r行，r＝1,2...w；

(3)交替迭代求解上述两个子目标函数，初始化循环次数t＝0,1,2...T，其中T为最大循环次数，将步骤(2)中的目标函数转化为如下两个交替优化的子目标函数：

$S^{t+1}=\arg \underset{{S\∈R}_{+}^{{N\×L}_c}}{\min }J(A^t,S),$

$A^{t+1}=\arg \underset{{A\∈R}_{+}^{M_h\×N}}{\min }J(A,S^{t+1}),$

其中J(A,S)表示步骤(2)中的目标函数。

交替迭代求解上述两个子目标函数，当光谱材料基矩阵A和光谱材料系数矩阵S的联合目标函数的循环次数达到T次时循环停止，得到优化的光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1；

(4)利用光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1，重建出高空间分辨率光谱图像Z^t+1＝A^t+1S^t+1。

本发明与现有技术相比具有以下优点：

第一，本发明由于在图像重建中，利用到了稀疏系数的结构相关性，能更加精确恢复光谱图像；

第二，本发明由于在参数调节中，对正则参数自适应的选取，使得重构算法更加鲁棒性。

附图说明

图1为本发明的实现流程图；

图2为本发明仿真实验所用图像balloons高分辨率的彩色图；

图3为用现有的MF方法对图像balloons在下采样为16*16时的恢复结果；

图4为用现有的CNMF方法对图像balloons在下采样为16*16时的恢复结果；

图5为用现有的PCA方法对图像balloons在下采样为16*16时的恢复结果；

图6为用现有的SNNMF方法对图像balloons在下采样为16*16时的恢复结果；

图7为用本发明方法对图像balloons在下采样为16*16时的恢复结果。

具体实施方式

参照图1，本发明基于非负结构稀疏表示的高光谱图像超分辨率方法，其实现步骤如下：

步骤1，输入低空间分辨率高光谱图像和高空间分辨率的彩色图像其中M_h表示谱段数，L_h表示低空间分辨率高光谱图像每个谱段的像素数，L_c表示高空间分辨率光谱图像每个谱段的像素数，且L_h＜L_c，M_c为彩色图像的谱段数，M_c＝3。

步骤2，假设高分辨率光谱图像Z可以表示为：Z＝AS，其中表示光谱材料基矩阵，表示光谱材料系数矩阵，其每一列表示各谱线在A上的稀疏分解系数。假设低分辨率光谱图像X和彩色图像Y与高分辨率光谱图像Z之间满足如下线性关系X＝ASG+W，Y＝FAS+V，其中表示下采样矩阵，表示彩色图像变换矩阵，和表示加性高斯白噪声。利用高光谱图像的局部与非局部相似性，构建如下基于结构稀疏光谱非负分解目标函数：

$(A,S)=\underset{A\∈R_{+}^{M_h\×N}}{\underset{S\∈R_{+}^{N\×L_C}}{\arg \min }}\frac{1}{2}{\left|\right|X-\mathrm{ASG}\left|\right|}_F^2+\frac{1}{2}{\left|\right|Y-\mathrm{FAS}\left|\right|}_F^2+\mathrm{\λ}\underset{k}{\mathrm{\Σ}}{\left|\right|{\mathrm{SR}}_k\left|\right|}_{\mathrm{1,2}},$

其中R_k表示相似谱线索引矩阵，SR_k表示从稀疏系数矩阵S中取出第k类中相似谱线对应的稀疏分解系数向量，即表示第k类中第p条相似谱线对应的稀疏系数，p＝1,2...n，n为第k类中相似谱线对应稀疏系数的个数，||？||2F表示矩阵的弗罗贝尼乌斯范数的平方，λ表示正则参数，||·||_1,2表示对矩阵每行求二范数再求和，即 ${\left|\right|{\mathrm{SR}}_k\left|\right|}_{\mathrm{1,2}}=\underset{r=1}{\overset{w}{\mathrm{\Σ}}}\sqrt{{\left(s_{k_1}^r\right)}^2+{\left(s_{k_2}^r\right)}^2+...+{\left(s_{k_n}^r\right)}^2,}$ w为SR_k矩阵的行数，r为SR_k矩阵的第r行，r＝1,2...w；

步骤3，通过交替优化方法求解步骤2中的目标函数：先利用初始光谱材料基矩阵优化求解出光谱材料系数矩阵S，再固定光谱材料系数矩阵S更新光谱材料基矩阵A，如此交替迭代，最终得到优化的光谱材料基矩阵和光谱材料系数矩阵。

(3a)初始化步骤1：利用双三次样条插值法，对低空间分辨率高光谱图像X进行插值，得到初始高空间分辨率高光谱图像Z⁰。结合初始恢复的光谱图像Z⁰和彩色图像Y对谱线进行近邻聚类：

$G_k=\left\{j\right|{\left|\right|{\tilde{z}}_j-{\tilde{z}}_i\left|\right|}_2^2<\mathrm{\&tau;}\},{\tilde{z}}_j={[{\left(z_j^t\right)}^T,y_j^T]}^T,$

其中，表示第t次迭代中高光谱图像的第j条谱线，y_j表示彩色图像的第j条谱线，(·)^T表示矩阵转置，j表示寻找到的相似谱线在三维图像块中的位置，τ为预设的相似度阀值。利用G_k，生成第k类相似谱线选择矩阵R_k；

(3b)初始化步骤2：通过交量最大化方法AVMAX(T.-H.Chan,W.-k.Ma等作者的“A simplex volume maximization framework for hyperspectral endmemberextraction”)得到初始光谱材料基矩阵A⁰；

(3c)初始化循环次数t＝0,1,2...T，通过交替优化将步骤2中的目标函数转化为如下两个子目标函数：

$S^{t+1}=\arg \underset{{S\&Element;R}_{+}^{{N\&times;L}_c}}{\min }J(A^t,S),$

$A^{t+1}=\arg \underset{{A\&Element;R}_{+}^{M_h\&times;N}}{\min }J(A,S^{t+1}),$

其中，S^t+1为第t+1次循环得到的优化光谱材料系数矩阵，A^t+1为第t+1次循环得到的优化光谱材料基矩阵，J(A,S)表示步骤2中的目标函数；

(3d)应用交替方向乘子法(alternating direction method of multipliers，ADMM)，将步骤(3c)中关于光谱材料系数矩阵S的优化问题转换成如下的优化问题：

$(M_k,S,E,Q)=\arg \underset{M_k,S,E,Q}{\min }{L}_{\mathrm{\&mu;},\mathrm{\&gamma;}}(M_k,S,E,Q),$

其中L_μ,γ(M_k,S,E,Q)表示增广拉格朗日函数，定义如下：

$\begin{array}{c}{L}_{\mathrm{\&mu;},\mathrm{\&gamma;}}(M_k,S,E,Q)=\frac{1}{2}{\left|\right|X-\mathrm{AQG}\left|\right|}_F^2+\frac{1}{2}{\left|\right|Y-\mathrm{FAQ}\left|\right|}_F^2+\mathrm{\&lambda;}\underset{k}{\mathrm{\&Sigma;}}{\left|\right|M_k\left|\right|}_{\mathrm{1,2}}+L_{+}\left(E\right)\\ +\frac{\mathrm{\&mu;}}{2}\underset{k}{\mathrm{\&Sigma;}}{\left|\right|M_k-{\mathrm{SR}}_k+\frac{H_1}{\mathrm{\&mu;}}\left|\right|}_F^2+\frac{\mathrm{\&gamma;}}{2}{\left|\right|E-S+\frac{H_2}{\mathrm{\&gamma;}}\left|\right|}_F^2+\frac{\mathrm{\&gamma;}}{2}{\left|\right|E-Q+\frac{H_3}{\mathrm{\&gamma;}}\left|\right|}_F^2\end{array},$

其中H₁、H₂、H₃分别为约束条件M_k＝SR_k、E＝S和Q＝E对应的拉格朗日乘子，μ和γ为惩罚因子，L₊(E)表示E的非负指示函数。

(3d1)初始化循环次数d＝0,1,2,...L，L表示最大循环次数，交替求解下列方程来求解步骤(3d)中的目标函数：

${m_i}^{d+1}=\frac{{g_i}^d}{{\left|\right|{g_i}^d\left|\right|}_2}\max \{{\left|\right|{g_i}^d\left|\right|}_2-\frac{\mathrm{\&lambda;}}{\mathrm{\&mu;}},0\},$ 其中 ${g_i}^d=S^d{R}_k-\frac{{h_1}^d}{\mathrm{\&mu;}},$

$S^{d+1}={(\underset{k}{\mathrm{\&Sigma;}}{\mathrm{\&mu;D}}_k{D_k}^T+\mathrm{\&gamma;I})}^{-1}(\underset{k}{\mathrm{\&Sigma;}}({{\mathrm{\&mu;D}}_k}^T{y^{d+1}}_k+{D_k}^T{h_1}^d)+{\mathrm{\&gamma;p}}^d+{h_2}^d),$ 其中 $D_k={R_k}^T\&CircleTimes;I,$

q^d+1＝(B₁ ^TB₁+γI)^-1(B₁ ^Tz+γp^d+h₃ ^d)，

e^d+1＝[s^d+1+q^d+1-γ^-1(h₂ ^d+h₃ ^d)]₊/2，

H₁ ^d+1＝H₁ ^d+μ(M_k ^d-S^dR_k)，

H₂ ^d+1＝H₂ ^d+γ(E^d-S^d)，

H₃ ^d+1＝H₃ ^d+γ(E^d-Q^d)，

μ＝min(βμ,10⁶)，其中β＝1.01，

γ＝min(βγ,10⁶)，

其中，m_i ^d+1为第d+1次迭代中第k类相似谱线系数矩阵M_k ^d+1的第i行向量，||·||₂表示向量的二范数，max{a,b}表示取a与b中的最大值；s^d+1为第d+1次迭代中光谱材料系数矩阵S^d+1的向量形式，h₁ ^d、h₂ ^d、h₃ ^d分别表示第d次迭代中拉格朗日乘子H₁ ^d、H₂ ^d、H₃ ^d的向量形式；q^d+1和e^d+1分别为第d+1次迭代中矩阵Q^d+1和E^d+1的向量形式， $B_1={\left[{(G^T\&CircleTimes;A)}^T{(I\&CircleTimes;\mathrm{FA})}^T\right]}^T\&Element;R^{(M_h{L}_h+3L_c)\&times;{\mathrm{NL}}_c},$ [·]₊表示取正数；

(3d2)当L_μ,γ的相对变化量小于阀值tol＝10^-5，停止迭代，得到优化的光谱材料稀疏系数矩阵S。

(3e)输入步骤(3d2)中得到的优化的光谱材料稀疏系数矩阵S，使用交替方向乘子法求解光谱材料基矩阵A，将步骤(3c)中的光谱材料基矩阵A的目标函数式改写为：

$(x,a)=\underset{a,x}{\arg \min }{L}_{{\mathrm{\&mu;}}_1}(x,a),$

其中(x,a)表示增广拉格朗日函数，定义如下：

$L_{{\mathrm{\&mu;}}_1}(x,a)=\frac{1}{2}{\left|\right|z-\mathrm{Ba}\left|\right|}_2^2+L_{+}\left(x\right)+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|x-a+\frac{h_4}{{\mathrm{\&mu;}}_1}\left|\right|}_2^2,$

其中x为交替方向乘子法中定义的向量，vec(·)表示向量算子， $B={\left[{({\left(\mathrm{SG}\right)}^T\&CircleTimes;I)}^T{(S^T\&CircleTimes;F)}^T\right]}^T\&Element;R^{(M_h{L}_h+3L_c)\&times;M_{h}N},$ 表示克罗内克积运算，a为光谱材料基矩阵A的向量形式，其中h₄为拉格朗日乘子，μ₁为惩罚因子。

(3e1)初始化循环次数c＝0,1,2,...L，交替求解下列方程来优化求解步骤(3e)中的目标函数：

x^c+1＝(a^c-h₄ ^c/μ₁)₊，

a^c+1＝(B^TB+μ₁I)^-1(B^Tz+μ₁x^c+1+h₄ ^c)，

h₄ ^c+1＝h₄ ^c+μ₁(x^c+1-a^c+1)，

μ₁＝min(βμ₁,10⁶)；

(3e2)当的相对变化量小于设定的阀值tol＝10^-5时，停止迭代，得出优化的光谱材料基矩阵A。

(3f)循环执行步骤(3d)和(3e)，当光谱材料基矩阵A和光谱材料系数矩阵S的联合目标函数的循环次数达到T次时循环停止，得到优化的光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1。

步骤4，利用步骤3f得到的光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1，得到重建的高分辨率高光谱图像Z^t+1＝A^t+1S^t+1。

本发明的效果可以通过如下仿真实验具体说明：

1.仿真条件：

1)仿真实验中的彩色图像变换矩阵F通过Nikon相机D700得到，下采样矩阵G为每16*16个像素求平均为一个像素的下采样方式；

2)仿真实验所用编程平台为MatlabR2012b；

3)仿真实验所用的高空间分辨率彩色图像如图2所示；

4)仿真实验所用的光谱数据为20个场景，谱段从400nm到700nm每10nm一个谱段，共31个谱段，每张高空间分辨率图像大小为512*512；

5)仿真实验中，采用峰值信噪比PSNR指标来评价压缩感知实验结果，峰值信噪比PSNR定义为：

$\mathrm{PSNR}=10{\log }_{10}\left(\frac{1}{M_h}\underset{i=1}{\overset{M_h}{\mathrm{\&Sigma;}}}\frac{\mathrm{MA}{X}_i^2}{{\mathrm{MSE}}_i}\right),$

其中，MAX_i和MSE_i为重构出来的高空间分辨率高光谱图像Z的每个谱段的最大像素值和均方误差。

2.仿真内容：

仿真1，采用现有的CNMF方法，对高光谱图像balloons在下采样为16*16时进行恢复，其恢复结果如图3所示；

仿真2，采用现有的MF方法，对高光谱图像balloons在下采样为16*16时进行恢复，其恢复结果如图4所示；

仿真3是采用现有的PCA方法，对高光谱图像balloons在下采样为16*16时进行恢复，其恢复结果如图5所示；

仿真4，采用现有的SNNMF方法，对高光谱图像balloons在下采样为16*16时进行恢复，其恢复结果如图6所示；

仿真5，采用本发明方法对高光谱图像balloons在下采样为16*16时进行恢复。其恢复结果如图7所示。

从图3-图7所显示的高光谱图像balloons的恢复结果可以看出，本发明的恢复出来的图像比其他方法恢复出来的图像更干净，清晰，图像边缘更锐利，视觉效果更好。

3.峰值信噪比PSNR对比

将现有的CNMF方法、MF方法、PCA方法、SNNMF方法和本发明方法分别对高光谱图像balloons进行重构仿真，得到的峰值信噪比PSNR如表1。

表1恢复图像的峰值信噪比PSNR值(单位dB)

从表1可以看出，本发明的峰值信噪比PSNR比MF和PCA在下采样为16*16的时候要平均高出1.7dB和2.6dB，比SNNMF高出1.5dB。

Claims (2)

1.一种基于非负结构稀疏的高光谱图像超分辨率算法，包括如下步骤：

(1)分别输入低分辨率高光谱图像和高分辨率的彩色图像其中，M_h表示谱段数，L_h表示低分辨率高光谱图像每个谱段的像素数，L_c表示高分辨率高光谱图像每个谱段的像素数，且L_h＜L_c，M_c为彩色图像的谱段数，M_c＝3；

(2)假设高分辨率光谱图像Z可以表示为：Z＝AS，其中表示光谱材料基矩阵，表示光谱材料系数矩阵，其每一列表示各谱线在A上的稀疏分解系数。假设低分辨率光谱图像X和彩色图像Y与高分辨率光谱图像Z之间满足如下线性关系X＝ASG+W，Y＝FAS+V，其中表示下采样矩阵，表示彩色图像变换矩阵，和表示加性高斯白噪声。利用高光谱图像的局部与非局部相似性，构建如下基于结构稀疏光谱非负分解目标函数：

$(A,S)=\underset{A\&Element;R_{+}^{M_h\&times;N}}{\underset{S\&Element;R_{+}^{N\&times;L_C}}{\arg \min }}\frac{1}{2}{\left|\right|X-\mathrm{ASG}\left|\right|}_F^2+\frac{1}{2}{\left|\right|Y-\mathrm{FAS}\left|\right|}_F^2+\mathrm{\&lambda;}\underset{k}{\mathrm{\&Sigma;}}{\left|\right|{\mathrm{SR}}_k\left|\right|}_{\mathrm{1,2}},$

其中R_k表示相似谱线索引矩阵，SR_k表示从稀疏系数矩阵S中取出第k类中相似谱线对应的稀疏分解系数向量，即表示第k类中第p条相似谱线对应的稀疏系数，p＝1,2...n，n为第k类中相似谱线对应稀疏系数的个数，表示矩阵的弗罗贝尼乌斯范数的平方，λ表示正则参数，||·||_1,2表示对矩阵每行求二范数再求和，即w为SR_k矩阵的行数，r为SR_k矩阵的第r行，r＝1,2...w；

(3)交替迭代求解上述两个子目标函数，初始化循环次数t＝0,1,2...T，其中T为最大循环次数，将步骤(2)中的目标函数转化为如下两个交替优化的子目标函数：

$S^{t+1}=\arg \underset{{S\&Element;R}_{+}^{{N\&times;L}_c}}{\min }J(A^t,S),$

$A^{t+1}=\arg \underset{{A\&Element;R}_{+}^{M_h\&times;N}}{\min }J(A,S^{t+1}),$

其中J(A,S)表示步骤(2)中的目标函数。

交替迭代求解上述两个子目标函数，当光谱材料基矩阵A和光谱材料系数矩阵S的联合目标函数的循环次数达到T次时循环停止，得到优化的光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1；

(4)利用光谱材料基矩阵A^t+1和光谱材料系数矩阵S^t+1，重建出高空间分辨率光谱图像Z^t+1＝A^t+1S^t+1。

2.根据权利要求1所述的方法，其特征在于所述步骤(3)中对两个子目标函数的求解，按如下步骤进行：

(3a)应用交替方向乘子法，将步骤(3)中关于光谱材料系数矩阵S的优化问题转换成如下的优化问题：

$(M_k,S,E,Q)=\arg \underset{M_k,S,E,Q}{\min }{L}_{\mathrm{\&mu;},\mathrm{\&gamma;}}(M_k,S,E,Q),$

其中L_μ,γ(M_k,S,E,Q)表示增广拉格朗日函数，定义如下：

$\begin{array}{c}{L}_{\mathrm{\&mu;},\mathrm{\&gamma;}}(M_k,S,E,Q)=\frac{1}{2}{\left|\right|X-\mathrm{AQG}\left|\right|}_F^2+\frac{1}{2}{\left|\right|Y-\mathrm{FAQ}\left|\right|}_F^2+\mathrm{\&lambda;}\underset{k}{\mathrm{\&Sigma;}}{\left|\right|M_k\left|\right|}_{\mathrm{1,2}}+L_{+}\left(E\right)\\ +\frac{\mathrm{\&mu;}}{2}\underset{k}{\mathrm{\&Sigma;}}{\left|\right|M_k-{\mathrm{SR}}_k+\frac{H_1}{\mathrm{\&mu;}}\left|\right|}_F^2+\frac{\mathrm{\&gamma;}}{2}{\left|\right|E-S+\frac{H_2}{\mathrm{\&gamma;}}\left|\right|}_F^2+\frac{\mathrm{\&gamma;}}{2}{\left|\right|E-Q+\frac{H_3}{\mathrm{\&gamma;}}\left|\right|}_F^2\end{array},$

其中H₁、H₂、H₃分别为约束条件M_k＝SR_k、E＝S和Q＝E对应的拉格朗日乘子，μ和γ为惩罚因子，L₊(E)表示E的非负指示函数。

(3a1)初始化循环次数d＝0,1,2,...L，L表示最大循环次数，交替求解下列方程来求解步骤(3a)中的目标函数：

${m_i}^{d+1}=\frac{{g_i}^d}{{\left|\right|{g_i}^d\left|\right|}_2}\max \{{\left|\right|{g_i}^d\left|\right|}_2-\frac{\mathrm{\&lambda;}}{\mathrm{\&mu;}},0\},$ 其中 ${g_i}^d=S^d{R}_k-\frac{{h_1}^d}{\mathrm{\&mu;}},$

$S^{d+1}={(\underset{k}{\mathrm{\&Sigma;}}{\mathrm{\&mu;D}}_k{D_k}^T+\mathrm{\&gamma;I})}^{-1}(\underset{k}{\mathrm{\&Sigma;}}({{\mathrm{\&mu;D}}_k}^T{y^{d+1}}_k+{D_k}^T{h_1}^d)+{\mathrm{\&gamma;p}}^d+{h_2}^d),$ 其中 $D_k={R_k}^T\&CircleTimes;I,$

q^d+1＝(B₁ ^TB₁+γI)^-1(B₁ ^Tz+γp^d+h₃ ^d)，

e^d+1＝[s^d+1+q^d+1-γ^-1(h₂ ^d+h₃ ^d)]₊/2，

H₁ ^d+1＝H₁ ^d+μ(M_k ^d-S^dR_k)，

H₂ ^d+1＝H₂ ^d+γ(E^d-S^d)，

H₃ ^d+1＝H₃ ^d+γ(E^d-Q^d)，

μ＝min(βμ,10⁶)，其中β＝1.01，

γ＝min(βγ,10⁶)，

其中，m_i ^d+1为第d+1次迭代中第k类相似谱线系数矩阵M_k ^d+1的第i行向量，||·||₂表示向量的二范数，max{a,b}表示取a与b中的最大值；s^d+1为第d+1次迭代中光谱材料系数矩阵S^d+1的向量形式，h₁ ^d、h₂ ^d、h₃ ^d分别表示第d次迭代中拉格朗日乘子H₁ ^d、H₂ ^d、H₃ ^d的向量形式；q^d+1和e^d+1分别为第d+1次迭代中矩阵Q^d+1和E^d+1的向量形式， $B_1={\left[{(G^T\&CircleTimes;A)}^T{(I\&CircleTimes;\mathrm{FA})}^T\right]}^T\&Element;R^{(M_h{L}_h+3L_c)\&times;{\mathrm{NL}}_c},$ [·]₊表示取正数；

(3a2)当L_μ,γ的相对变化量小于阀值tol＝10^-5，停止迭代，得到优化的光谱材料稀疏系数矩阵S。

(3b)输入步骤(3a2)中得到的优化的光谱材料稀疏系数矩阵S，使用交替方向乘子法求解光谱材料基矩阵A，将步骤(3c)中的光谱材料基矩阵A的目标函数式改写为：

$(x,a)=\underset{a,x}{\arg \min }{L}_{{\mathrm{\&mu;}}_1}(x,a),$

其中(x,a)表示增广拉格朗日函数，定义如下：

$L_{{\mathrm{\&mu;}}_1}(x,a)=\frac{1}{2}{\left|\right|z-\mathrm{Ba}\left|\right|}_2^2+L_{+}\left(x\right)+\frac{{\mathrm{\&mu;}}_1}{2}{\left|\right|x-a+\frac{h_4}{{\mathrm{\&mu;}}_1}\left|\right|}_2^2,$

其中x为交替方向乘子法中定义的向量，vec(·)表示向量算子， $B={\left[{({\left(\mathrm{SG}\right)}^T\&CircleTimes;I)}^T{(S^T\&CircleTimes;F)}^T\right]}^T\&Element;R^{(M_h{L}_h+3L_c)\&times;M_{h}N},$ 表示克罗内克积运算，a为光谱材料基矩阵A的向量形式，其中h₄为拉格朗日乘子，μ₁为惩罚因子。

(3b1)初始化循环次数c＝0,1,2,...L，交替求解下列方程来优化求解步骤(3e)中的目标函数：

x^c+1＝(a^c-h₄ ^c/μ₁)₊，

a^c+1＝(B^TB+μ₁I)^-1(B^Tz+μ₁x^c+1+h₄ ^c)，

h₄ ^c+1＝h₄ ^c+μ₁(x^c+1-a^c+1)，

μ₁＝min(βμ₁,10⁶)；

(3b2)当的相对变化量小于设定的阀值tol＝10^-5时，停止迭代，得出优化的光谱材料基矩阵A。