CN104732566B - 基于非分离稀疏先验的高光谱图像压缩感知方法 - Google Patents

Abstract

本发明公开了一种基于非分离稀疏先验的高光谱图像压缩感知方法，用于解决现有高光谱图像压缩感知方法重建精度低的技术问题。技术方案是采集每个像素光谱的少量线性观测值作为压缩数据，实现大幅数据压缩的同时降低了图像采集过程中的资源需求。重建过程中，利用经验贝叶斯推理，构建稀疏信号的非分离稀疏先验，充分考量了稀疏信号内部非零元素之间的潜在相关性，实现了高光谱图像的高精度重建。由于该方法使用小波正交基作为字典，消除了对端元的依赖性。此外，基于贝叶斯框架的推理，实现了所有未知参数的全自动估计，无需人为调节，适应性广。试验表明，当采样率为0.1时，本发明获得的峰值信噪比相对于背景技术压缩感知方法提升4db以上。

Description

基于非分离稀疏先验的高光谱图像压缩感知方法

技术领域

本发明涉及一种高光谱图像压缩感知方法，特别是涉及一种基于非分离稀疏先验的高光谱图像压缩感知方法。

背景技术

高光谱图像成百上千的波段中存储着场景的光谱信息，有助于遥感地物的检测、分类与识别。然而，丰富的光谱信息导致高光谱图像数据量巨大，图像的获取、传输和处理需要消耗大量资源，制约了高光谱图像的应用。因此，研究高效的高光谱图像压缩算法是高光谱领域的热点问题之一。目前，普通图像的经典压缩算法已经成功推广到高光谱图像，来同时消除高光谱图像中波段内和波段间的冗余。然而，这些压缩算法均作用在图像获取之后，无法减少成像过程中的巨大资源需求。近年来，压缩感知成像理论表明仅需要从场景中采集少量的线性观测值便可以在需要时对原始场景图像进行高精度的重建，大幅减少图像采集过程中的资源消耗。ChengBo Li等人在文献“A compressive sensing and unmixingscheme for hyperspectral data processing,IEEE Transactions on ImageProcessing,2012,21(3):1200–1210”中，利用压缩感知成像技术，仅采集每个波段的少量线性观测，来实现对高光谱数据的大幅压缩。重建过程中，结合光谱线性混合模型，在少量端元的辅助下，引入全变差梯度稀疏约束，重建空间连续的丰度值矩阵，最终线性混合端元和丰度值矩阵得到原始的高光谱图像。然而，该方法存在以下几个问题：首先，提出的稀疏约束并没有考虑稀疏信号内部非零元素之间的关系，重建精度受限；其次，算法性能严重依赖端元的选择，但文中并未提供自适应的端元选择机制；此外，针对不同的数据需要参数调节，缺乏适应性。

发明内容

为了克服现有高光谱图像压缩感知方法重建精度低的不足，本发明提供一种基于非分离稀疏先验的高光谱图像压缩感知方法。该方法采集每个像素光谱的少量线性观测值作为压缩数据，实现大幅数据压缩的同时也进一步降低了图像采集过程中的资源需求。重建过程中，利用经验贝叶斯推理，构建稀疏信号的非分离稀疏先验，充分考量了稀疏信号内部非零元素之间的潜在相关性，实现了高光谱图像的高精度重建。由于该方法使用小波正交基作为字典，消除了对端元的依赖性。此外，基于贝叶斯框架的推理，实现了所有未知参数的全自动估计，无需人为调节，适应性广。在真实的卫星图像URBAN，PAVIA UNIVERSITY以及INDIANA数据集上的试验结果表明，当采样率为0.1时，本发明获得的峰值信噪比相对于背景技术压缩感知方法提升4db以上。

本发明解决其技术问题所采用的技术方案是：一种基于非分离稀疏先验的高光谱图像压缩感知方法，其特点是包括以下步骤：

步骤一、对于包含n_b个波段，每个波段包含n_p个像素的高光谱图像，将每一个波段拉伸成为一个行向量，所有的行向量组成一个二维矩阵X的每一列表示每一个像素对应的光谱，称为光谱维；X的每一行对应一个波段的所有像素值，称为空间维。

步骤二、使用满足高斯分布，列归一化的随机观测矩阵对高光谱图像的光谱维进行随机采样，得到压缩数据m_b表示n_b个波段压缩后的长度，m_b<n_b，定义采样率为ρ＝m_b/n_b。

G＝AX+N (1)

其中，表示压缩感知采样过程中存在的噪声。

步骤三、由于高光谱图像数据本身并不稀疏，通过引入Haar小波基作为字典对高光谱图像的每一个光谱进行稀疏化，X＝DY，其中D为小波正交基，Y为列稀疏的小波系数矩阵。因此，(1)式中的压缩感知模型进一步表示为G＝AX+N＝ADY+N。假设压缩感知采样过程中存在的噪声服从均值为0的矩阵正太分布则得到压缩感知模型对应的似然函数

其中，Σ_n＝diag(λ)表示以λ的元素为对角线元素的对角矩阵，指示噪声的强度。表示对矩阵Q的加权迹范数。对于稀疏信号Y，假设服从如下的均值为0的矩阵正太分布

其中，控制着Y中每一行的稀疏度。当γ_i＝0表示Y的第i行为0。λ和γ均为待估计的参数。

步骤四、由于噪声强度λ和稀疏先验中超参数γ未知，无法采用最大后验估计直接对稀疏信号Y进行估计。因此，根据经验贝叶斯框架，先基于压缩数据G采用MAP过程估计出未知的参数λ和γ，如下

假设p(γ)为常数，表示扁平先验，在(4)式中省略。通过积分，并引入-2log运算，容易得知(4)式等价于最小化的(5)式

其中，Σ_by＝Σ_n+ADΣ_yD^TA^T，为关于λ和γ的代价函数。

步骤五、通过对(5)式进行适当的变形，得到稀疏信号Y的非分离稀疏约束模型。首先，对(5)式的第一份部分进行变形

然后，将(6)式带入到(5)式中，得

接着，引入新的代价方程如下

显然，而且能够证明最小化(5)式，再进行关于稀疏信号Y的MAP估计，与直接最小化(8)式得到关于λ和γ解相同，关于Y的解仅相差一个常量因此，(8)式看作是一个全新的关于稀疏信号的正则化回归模型，其中为稀疏信号的非分离稀疏约束。该约束不能拆分成关于Y中每一行的独立约束。因此，该约束能够同时约束稀疏信号中非零元素，潜在地考量这些元素之间的相关性。

步骤六、已知压缩后数据G和随机观测矩阵A，采用坐标下降法最小化(8)式，每次仅优化一个未知数而固定其他的所有未知数。具体步骤如下：

①初始化，λ⁰和γ⁰均初始化为对应长度的全1向量，计数变量t＝0；

②更新中间变量Σ_n＝diag(λ^t)，Σ_y＝diag(γ^t)，Σ_by＝Σ_n+ADΣ_yD^TA^T；

③固定λ^t和γ^t，根据(8)式得到关于Y的优化形式，如下

求解得到Y的更新规则如下，

④固定Y^t+1和λ^t，得到关于γ的优化形式，如下

求解得到如下的更新形式：

γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T). (12)

其中，γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T)代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，

⑤固定Y^t+1和γ^t+1，得到关于λ的优化形式，如下

求解得到如下的更新形式：

其中，根号运算表示向量每一个元素开方后组成的向量，./运算代表两个向量对应元素相除后组成的向量，代表对角线元素组成的向量。

⑥计算稀疏信号Y更新前后的差异，如下

其中，表示对Y^t+1内的每一个元素乘以||·||_F表示弗罗贝尼乌斯范数，如果t>400或者η<10^-4，则退出循环；否则，循环执行步骤②至⑥。

⑦假设上述循环结束得到的最优稀疏信号为Y_rec，则待重建的高光谱图像X_rec通过如下的方式得到：

本发明的有益效果是：该方法采集每个像素光谱的少量线性观测值作为压缩数据，实现大幅数据压缩的同时也进一步降低了图像采集过程中的资源需求。重建过程中，利用经验贝叶斯推理，构建稀疏信号的非分离稀疏先验，充分考量了稀疏信号内部非零元素之间的潜在相关性，实现了高光谱图像的高精度重建。由于该方法使用小波正交基作为字典，消除了对端元的依赖性。此外，基于贝叶斯框架的推理，实现了所有未知参数的全自动估计，无需人为调节，适应性广。在真实的卫星图像URBAN，PAVIA UNIVERSITY以及INDIANA数据集上的试验结果表明，当采样率为0.1时，本发明获得的峰值信噪比相对于背景技术压缩感知方法提升4db以上。

下面结合具体实施方式对本发明作详细说明。

具体实施方式

本发明基于非分离稀疏先验的高光谱图像压缩感知方法具体步骤如下：

在本发明中，为了便于处理，对于包含n_b个波段，每个波段包含n_p个像素的高光谱图像，将每一个波段拉伸成为一个行向量，所有的行向量组成一个二维矩阵X的每一列表示每一个像素对应的光谱，该方向为光谱维；X的每一行对应一个波段的所有像素值，该方向为空间维。在压缩过程中，本发明随机采样高光谱图像的光谱维，获得少量的线性观测值作为压缩数据；重建过程中，构建贝叶斯压缩感知模型；然后，采用经验贝叶斯推理构建稀疏信号的非分离稀疏先验；最后，在非分离稀疏先验的约束下，估计稀疏信号，重建原始图像。具体步骤如下：

1、获取压缩数据。

使用满足高斯分布，列归一化的随机观测矩阵对高光谱图像的光谱维进行随机采样，得到压缩数据m_b表示n_b个波段压缩后的长度，m_b<n_b，定义采样率为ρ＝m_b/n_b。

G＝AX+N (1)

其中，表示压缩感知采样过程中存在的噪声。

2、建立贝叶斯压缩感知模型。

由于高光谱图像数据本身并不稀疏，本发明通过引入Haar小波基作为字典对高光谱图像的每一个光谱进行稀疏化，X＝DY，其中D为小波正交基，Y为列稀疏的小波系数矩阵。因此，(1)式中的压缩感知模型可以进一步表示为G＝AX+N＝ADY+N。假设压缩感知采样过程中存在的噪声服从均值为0的矩阵正太分布(Matrix Normal distribution)则可以得到压缩感知模型对应的似然函数

其中，Σ_n＝diag(λ)表示以λ的元素为对角线元素的对角矩阵，指示噪声的强度。表示对矩阵Q的加权迹范数。对于稀疏信号Y，假设服从如下的均值为0的矩阵正太分布

其中，控制着Y中每一行的稀疏度。当γ_i＝0表示Y的第i行为0。在本发明中，λ和γ均为待估计的参数。

3、经验贝叶斯推理。

由于噪声强度λ和稀疏先验中超参数γ未知，无法采用最大后验估计(Maximum aposterior estimation，MAP)直接对稀疏信号Y进行估计。因此，本发明根据经验贝叶斯框架，先基于压缩数据G采用MAP过程估计出未知的参数λ和γ，如下

假设p(γ)为常数，表示扁平先验，因此可以在上式子中省略。通过积分，并引入-2log运算，容易得知(4)式等价于最小化的(5)式

其中，Σ_by＝Σ_n+ADΣ_yD^TA^T，为关于λ和γ的代价函数。

4、建立非分离稀疏先验约束模型。

通过对(5)式进行适当的变形，便可以得到稀疏信号Y的非分离稀疏约束模型。首先，对(5)式的第一份部分进行变形

然后，将(6)式带入到(5)式中，得

接着，引入新的代价方程如下

显然，而且可以证明最小化(5)式，再进行关于稀疏信号Y的MAP估计，与直接最小化(8)式得到关于λ和γ解相同，关于Y的解仅相差一个常量因此，(8)式看作是一个全新的关于稀疏信号的正则化回归模型，其中为稀疏信号的非分离稀疏约束。该约束不能拆分成关于Y中每一行的独立约束。因此，该约束能够同时约束稀疏信号中非零元素，潜在地考量这些元素之间的相关性。

5、模型求解。

已知压缩后数据G和随机观测矩阵A，本发明采用坐标下降法最小化(8)式，每次仅优化一个未知数而固定其他的所有未知数。具体步骤如下：

⑧初始化，λ⁰和γ⁰均初始化为对应长度的全1向量，计数变量t＝0；

⑨更新中间变量Σ_n＝diag(λ^t)，Σ_y＝diag(γ^t)，Σ_by＝Σ_n+ADΣ_yD^TA^T；

⑩固定λ^t和γ^t，根据(8)式得到关于Y的优化形式，如下

求解得到Y的更新规则如下，

固定Y^t+1和λ^t，得到关于γ的优化形式，如下

求解得到如下的更新形式：

γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T). (12)

其中，γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T)代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，

固定Y^t+1和γ^t+1，得到关于λ的优化形式，如下

求解得到如下的更新形式：

其中，根号运算表示向量每一个元素开方后组成的向量，./运算代表两个向量对应元素相除后组成的向量，代表对角线元素组成的向量。

计算稀疏信号Y更新前后的差异，如下

其中，表示对Y^t+1内的每一个元素乘以||·||_F表示弗罗贝尼乌斯范数(Frobenius norm),如果t>400或者η<10^-4，则退出循环；否则，循环执行步骤⑨至

假设上述循环结束得到的最优稀疏信号为Y_rec，则待重建的高光谱图像X_rec可以通过如下的方式得到：

Claims (1)

1.一种基于非分离稀疏先验的高光谱图像压缩感知方法，其特征在于包括以下步骤：

步骤一、对于包含n_b个波段，每个波段包含n_p个像素的高光谱图像，将每一个波段拉伸成为一个行向量，所有的行向量组成一个二维矩阵X的每一列表示每一个像素对应的光谱，称为光谱维；X的每一行对应一个波段的所有像素值，称为空间维；

步骤二、使用满足高斯分布，列归一化的随机观测矩阵对高光谱图像的光谱维进行随机采样，得到压缩数据m_b表示n_b个波段压缩后的长度，m_b<n_b，定义采样率为ρ＝m_b/n_b；

G＝AX+N (1)

其中，表示压缩感知采样过程中存在的噪声；

步骤三、由于高光谱图像数据本身并不稀疏，通过引入Haar小波基作为字典对高光谱图像的每一个光谱进行稀疏化，X＝DY，其中D为小波正交基，Y为列稀疏的小波系数矩阵；因此，(1)式中的压缩感知模型进一步表示为G＝AX+N＝ADY+N；假设压缩感知采样过程中存在的噪声服从均值为0的矩阵正态分布则得到压缩感知模型对应的似然函数

其中，Σ_n＝diag(λ)表示以λ的元素为对角线元素的对角矩阵，指示噪声的强度；表示对矩阵Q的加权迹范数；对于稀疏信号Y，假设服从如下的均值为0的矩阵正态分布

<mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>Y</mi> <mo>|</mo> <mi>&amp;gamma;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mi>exp</mi> <mo>{</mo> <mo>-</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>Y</mi> <mo>|</mo> <msubsup> <mo>|</mo> <msub> <mi>&amp;Sigma;</mi> <mi>y</mi> </msub> <mn>2</mn> </msubsup> <mo>}</mo> </mrow> <mrow> <msup> <mrow> <mo>(</mo> <mn>2</mn> <mi>&amp;pi;</mi> <mo>)</mo> </mrow> <mrow> <msub> <mi>n</mi> <mi>b</mi> </msub> <msub> <mi>n</mi> <mi>p</mi> </msub> <mo>/</mo> <mn>2</mn> </mrow> </msup> <mo>|</mo> <msub> <mi>&amp;Sigma;</mi> <mi>y</mi> </msub> <msup> <mo>|</mo> <mrow> <msub> <mi>n</mi> <mi>p</mi> </msub> <mo>/</mo> <mn>2</mn> </mrow> </msup> </mrow> </mfrac> <mo>,</mo> <msub> <mi>&amp;Sigma;</mi> <mi>y</mi> </msub> <mo>=</mo> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>g</mi> <mrow> <mo>(</mo> <mi>&amp;gamma;</mi> <mo>)</mo> </mrow> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>

其中，控制着Y中每一行的稀疏度；当γ_i＝0表示Y的第i行为0；λ和γ均为待估计的参数；

步骤四、由于噪声强度λ和稀疏先验中超参数γ未知，无法采用最大后验估计直接对稀疏信号Y进行估计；因此，根据经验贝叶斯框架，先基于压缩数据G采用MAP过程估计出未知的参数λ和γ，如下

<mrow> <mo>{</mo> <msub> <mi>&amp;lambda;</mi> <mrow> <mi>o</mi> <mi>p</mi> <mi>t</mi> </mrow> </msub> <mo>,</mo> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>o</mi> <mi>p</mi> <mi>t</mi> </mrow> </msub> <mo>}</mo> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mi>&amp;lambda;</mi> <mo>,</mo> <mi>&amp;gamma;</mi> </mrow> </munder> <mi>p</mi> <mrow> <mo>(</mo> <mi>&amp;lambda;</mi> <mo>,</mo> <mi>&amp;gamma;</mi> <mo>|</mo> <mi>G</mi> <mo>)</mo> </mrow> <mo>&amp;Proportional;</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <mi>&amp;lambda;</mi> <mo>,</mo> <mi>&amp;gamma;</mi> </mrow> </munder> <mo>&amp;Integral;</mo> <mi>p</mi> <mrow> <mo>(</mo> <mi>G</mi> <mo>|</mo> <mi>Y</mi> <mo>,</mo> <mi>&amp;lambda;</mi> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>Y</mi> <mo>|</mo> <mi>&amp;gamma;</mi> <mo>)</mo> </mrow> <mi>p</mi> <mrow> <mo>(</mo> <mi>&amp;gamma;</mi> <mo>)</mo> </mrow> <mi>d</mi> <mi>Y</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>

假设p(γ)为常数，表示扁平先验，在(4)式中省略；通过积分，并引入-2log运算，容易得知(4)式等价于最小化的(5)式

其中，Σ_by＝Σ_n+ADΣ_yD^TA^T，为关于λ和γ的代价函数；

步骤五、通过对(5)式进行适当的变形，得到稀疏信号Y的非分离稀疏约束模型；首先，对(5)式的第一份部分进行变形

<mrow> <mi>t</mi> <mi>r</mi> <mrow> <mo>(</mo> <mrow> <msubsup> <mi>n</mi> <msub> <mi>n</mi> <mi>p</mi> </msub> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <msup> <mi>G</mi> <mi>T</mi> </msup> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>b</mi> <mi>y</mi> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mi>G</mi> </mrow> <mo>)</mo> </mrow> <mo>=</mo> <mi>t</mi> <mi>r</mi> <mrow> <mo>&amp;lsqb;</mo> <mrow> <mrow> <mo>(</mo> <mfrac> <msup> <mi>G</mi> <mi>T</mi> </msup> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> </mfrac> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>b</mi> <mi>y</mi> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mrow> <mo>(</mo> <mfrac> <mi>G</mi> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> </mfrac> <mo>)</mo> </mrow> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mo>=</mo> <munder> <mi>min</mi> <mi>Y</mi> </munder> <mo>|</mo> <mo>|</mo> <mi>A</mi> <mi>D</mi> <mi>Y</mi> <mo>-</mo> <mfrac> <mi>G</mi> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> </mfrac> <mo>|</mo> <msubsup> <mo>|</mo> <msub> <mi>&amp;Sigma;</mi> <mi>n</mi> </msub> <mn>2</mn> </msubsup> <mo>+</mo> <mo>|</mo> <mo>|</mo> <mi>Y</mi> <mo>|</mo> <msubsup> <mo>|</mo> <msub> <mi>&amp;Sigma;</mi> <mi>y</mi> </msub> <mn>2</mn> </msubsup> <mo>.</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>

然后，将(6)式带入到(5)式中，得

接着，引入新的代价方程如下

显然，证明了最小化(5)式，再进行关于稀疏信号Y的MAP估计，与直接最小化(8)式得到关于λ和γ解相同，关于Y的解仅相差一个常量因此，(8)式看作是一个全新的关于稀疏信号的正则化回归模型，其中为稀疏信号的非分离稀疏约束；该约束不能拆分成关于Y中每一行的独立约束；因此，该约束能够同时约束稀疏信号中非零元素，潜在地考量这些元素之间的相关性；

步骤六、已知压缩后数据G和随机观测矩阵A，采用坐标下降法最小化(8)式，每次仅优化一个未知数而固定其他的所有未知数；具体步骤如下：

①初始化，λ⁰和γ⁰均初始化为对应长度的全1向量，计数变量t＝0；

②更新中间变量Σ_n＝diag(λ^t)，Σ_y＝diag(γ^t)，Σ_by＝Σ_n+ADΣ_yD^TA^T；

③固定λ^t和γ^t，根据(8)式得到关于Y的优化形式，如下

求解得到Y的更新规则如下，

<mrow> <msup> <mi>Y</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <msub> <mi>&amp;Sigma;</mi> <mi>y</mi> </msub> <msup> <mi>D</mi> <mi>T</mi> </msup> <msup> <mi>A</mi> <mi>T</mi> </msup> <msubsup> <mi>&amp;Sigma;</mi> <mrow> <mi>b</mi> <mi>y</mi> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msubsup> <mfrac> <mi>G</mi> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>

④固定Y^t+1和λ^t，得到关于γ的优化形式，如下

求解得到如下的更新形式：

γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T). (12)

其中，γ^t+1＝diag(V^T+Y^t+1(Y^t+1)^T)代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，

⑤固定Y^t+1和γ^t+1，得到关于λ的优化形式，如下

求解得到如下的更新形式：

<mrow> <msup> <mi>&amp;lambda;</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <msqrt> <mrow> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>g</mi> <mrow> <mo>(</mo> <msup> <mi>QQ</mi> <mi>T</mi> </msup> <mo>)</mo> </mrow> <mo>.</mo> <mo>/</mo> <mi>&amp;alpha;</mi> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>

其中，根号运算表示向量每一个元素开方后组成的向量，./运算代表两个向量对应元素相除后组成的向量，代表对角线元素组成的向量；

⑥计算稀疏信号Y更新前后的差异，如下

<mrow> <mi>&amp;eta;</mi> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msup> <mi>Y</mi> <mrow> <mi>t</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>.</mo> <mo>*</mo> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> <mo>|</mo> <msub> <mo>|</mo> <mi>F</mi> </msub> </mrow> <mrow> <mo>|</mo> <mo>|</mo> <msup> <mi>Y</mi> <mi>t</mi> </msup> <mo>.</mo> <mo>*</mo> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> <mo>|</mo> <msub> <mo>|</mo> <mi>F</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>

其中，表示对Y^t+1内的每一个元素乘以||·||_F表示弗罗贝尼乌斯范数，如果t>400或者η<10^-4，则退出循环；否则，循环执行步骤②至⑥；

⑦假设上述循环结束得到的最优稀疏信号为Y_rec，则待重建的高光谱图像X_rec通过如下的方式得到：

<mrow> <msub> <mi>X</mi> <mrow> <mi>r</mi> <mi>e</mi> <mi>c</mi> </mrow> </msub> <mo>=</mo> <mi>D</mi> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mrow> <mi>r</mi> <mi>e</mi> <mi>c</mi> </mrow> </msub> <mo>.</mo> <mo>*</mo> <msqrt> <msub> <mi>n</mi> <mi>p</mi> </msub> </msqrt> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> <mo>.</mo> </mrow> 3