CN105427351A - 基于流形结构化稀疏先验的高光谱图像压缩感知方法 - Google Patents

Abstract

本发明公开了一种基于流形结构化稀疏先验的高光谱图像压缩感知方法，用于解决现有高光谱图像压缩感知方法精度低的技术问题。技术方案是随机采样每个像素光谱的少量线性观测值作为压缩数据，通过流形结构化稀疏先验，同时刻画高光谱图像稀疏化后光谱维中的稀疏性和空间维中的流形结构；通过隐变量贝叶斯模型，将信号重建，稀疏先验学习以及噪声估计统一到一个正则化回归模型进行优化求解。学习得到的稀疏先验既能充分地刻画高光谱图像的三维结构，又具有较强的噪声鲁棒性。利用该稀疏先验，实现了高光谱图像的高精度重建。据测试，当在压缩数据中加入高斯白噪声使得压缩数据信噪比为15db，采样率为0.09时，获得了23db的峰值信噪比。

Description

基于流形结构化稀疏先验的高光谱图像压缩感知方法

技术领域

本发明涉及一种高光谱图像压缩感知方法，特别是涉及一种基于流形结构化稀疏先验的高光谱图像压缩感知方法。

背景技术

高光谱图像包含成百上千的波段，每个像素包含一条连续的光谱。丰富的光谱信息使得高光谱图像在目标检测、识别等方面具有极大的优势，然而，高光谱图像巨大的数据量对图像的采集、传输和处理提出了苛刻的要求，制约了其实际应用。因此，高光谱图像压缩是高光谱领域的热点研究之一。压缩感知成像理论证明仅需要采集少量的线性观测值便可以精确重建原始场景的图像。相对于传统的图像压缩算法，极大地减少了成像过程中的资源消耗。

ChengBoLi等人在文献“Acompressivesensingandunmixingschemeforhyperspectraldataprocessing,IEEETransactionsonImageProcessing,2012,21(3):1200–1210”中公开了一种高效的高光谱图像压缩感知算法。成像过程中，使用单像素相机采集少量线性观测值作为压缩数据。重建过程中，基于线性混合模型，引入适量的端元光谱，重建空间梯度稀疏的丰度值矩阵。最后，通过线性混合重建的丰度值矩阵和引入的端元光谱重建高光谱图像。然而，该算法仅考虑了空间的稀疏性，未能充分利用高光谱图像的三维结构，重建精度受限；其次，该算法的稀疏性约束噪声鲁棒性差；此外，算法性能严重依赖端元光谱的选择，实用性受限。

发明内容

为了克服现有高光谱图像压缩感知方法精度低的不足，本发明提供一种基于流形结构化稀疏先验的高光谱图像压缩感知方法。该方法随机采样每个像素光谱的少量线性观测值作为压缩数据，通过流形结构化稀疏先验，同时刻画高光谱图像稀疏化后光谱维中的稀疏性和空间维中的流形结构；通过隐变量贝叶斯模型，将信号重建，稀疏先验学习以及噪声估计统一到一个正则化回归模型进行优化求解。学习得到的稀疏先验既能充分地刻画高光谱图像的三维结构，又具有较强的噪声鲁棒性。利用该稀疏先验，实现了高光谱图像的高精度重建。在真实的高光谱遥感数据集Urban上的实验结果表明，当在压缩数据中加入高斯白噪声使得压缩数据信噪比为15db，采样率为0.09时，获得了23db的峰值信噪比。

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

步骤一、针对包含n_b个波段，每个波段包含n_r行和n_c列的高光谱图像，将每一个波段拉伸成为一个行向量，重新组成一个二维矩阵，其中，X的每一列对应每个像素的光谱；每一行对应每个波段的所有像素值。将行和列分别称为空间维和光谱维。

获取压缩数据过程中，利用列归一化的高斯随机观测矩阵随机采样高光谱图像X的光谱维，获得压缩数据m_b为压缩后波段长度。

F＝AX+N(1)

其中，表示采样中的噪声。采样率ρ定义为ρ＝m_b/n_b。

步骤二、利用Haar小波基对高光谱图像的每个光谱进行稀疏化，如X＝ΨY，Ψ为小波基，Y为列稀疏的系数矩阵，模型(1)表示为F＝AΨY+N。假设采样过程中噪声N服从的矩阵正太分布，I为对应大小的单位矩阵，模型(1)对应的似然函数定义为

$p\left(F\right|Y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{n}_p/2}|{\mathrm{\Σ}}_n{|}^{n_p/2}}---\left(2\right)$

其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵。表示Q矩阵的加权迹范数。

除过列稀疏性，高光谱图像空间像素之间的相似性使得Y中不同列的稀疏信号位于一个结构未知的流形结构上。为了充分描述Y的特性，假设Y服从如下的矩阵正太分布

$p\left(Y\right|{\mathrm{\Σ}}_{ry},{\mathrm{\Σ}}_{cy})=\frac{\exp \{-\frac{1}{2}tr\left({\mathrm{\Σ}}_{cy}^{-1}{Y}^T{\mathrm{\Σ}}_{ry}^{-1}Y\right)\}}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}\left|{\mathrm{\Σ}}_{ry}{|}^{n_p/2}\right|{\mathrm{\Σ}}_{cy}{|}^{n_b/2}},---\left(3\right)$

为描述Y中列信号的稀疏性，令Σ_y＝diag(γ)表示以γ为对角线元素的对角矩阵， $\mathrm{\γ}={\[{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_{n_b}\]}^T.$ 令 $\mathrm{\κ}={\[{\mathrm{\κ}}_1,\mathrm{...},{\mathrm{\κ}}_{n_b}\]}^T,$ 假设γ服从如下的伽玛分布

$p\left(\mathrm{\γ}\right|\mathrm{\κ})=\underset{i=1}{\overset{n_b}{\Π}}Gamma(1,\frac{2}{{\mathrm{\κ}}_i})=\underset{i=1}{\overset{n_b}{\Π}}\frac{{\mathrm{\κ}}_i}{2}\exp (-\frac{{\mathrm{\κ}}_i{\mathrm{\γ}}_i}{2})---\left(4\right)$

式子(3)中，Σ_cy描述Y中不同列信号之间的相关性，因此，式子(3)隐式地表示Y中不同的稀疏信号之间存在的流形结构。为了更加灵活地学习Σ_cy，进一步假设Σ_cy服从如下的反威沙特分布

其中，l是给定的常量，表示自由度，Γ_np是多变量伽玛函数，为参考协方差矩阵。该先验通过最小化Σ_cy和Θ之间的布雷格曼散度，使得Σ_cy趋近于Θ，从而减轻了Σ_cy学习过程中的过拟合问题。

步骤三、为使得流行结构化稀疏先验能够更好地匹配图像分布并具有较强的噪声鲁棒性，通过隐变量贝叶斯模型对噪声参数λ和先验参数γ，κ，Σ_cy和Θ进行估计。令f＝vec(F)，y＝vec(Y)，n＝vec(N)和vec(Q)表示将矩阵Q拉成列向量，表示克罗内克积，则模型(2)等价于

$p\left(f\right|y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||f-\mathrm{\Φ}y|{|}_{I\⊗{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{m}_p/2}|I\⊗{\mathrm{\Σ}}_n{|}^{1/2}}---\left(6\right)$

同样，模型(3)中关于Y的先验等价于

$p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})=\frac{\exp (-\frac{1}{2}{y}^T{\mathrm{\Σ}}_y^{-1}y)}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}|{\mathrm{\Σ}}_y{|}^{1/2}},{\mathrm{\Σ}}_y={\mathrm{\Σ}}_{cy}\⊗{\mathrm{\Σ}}_{ry}.---\left(7\right)$

根据公式(6)、公式(7)，所有的未知参数通过求解如下的优化问题得到

$\underset{\mathrm{\λ},\mathrm{\γ}\≥0,\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}}{\max }p(\mathrm{\λ},\mathrm{\γ},\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}|f)\∝\∫p\left(f\right|y,\mathrm{\λ})p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})p\left(\mathrm{\γ}\right|\mathrm{\κ})p\left({\mathrm{\Σ}}_{cy}\right|\mathrm{\Θ},l)dy---\left(8\right)$

通过积分，并引入-2log运算，容易得知式子(8)，等价于最小化如下的式子

其中，tr(·)表示迹范数，对式子(9)的第一项做如下变形

$f^T{\mathrm{\Σ}}_{by}^{-1}f=\underset{y}{\min }\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(10\right)$

将式子(10)带入到式子(9)中，得到如下等价于式子(8)的正则化回归模型

$\begin{array}{c}\underset{y,\mathrm{\λ}\≥0,\mathrm{\γ}\≥0,\mathrm{\κ},{\Σ}_{cy},\mathrm{\Θ}}{\min }\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\Σ}_n}^2+y^T{\Σ}_y^{-1}y+\log |{\Σ}_{by}|+\underset{i=1}{\overset{n_b}{\Σ}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\γ}}_i)\\ +tr\left(\mathrm{\Θ}{\Σ}_{cy}^{-1}\right)+(n_p+l+1)\log \left|{\Σ}_{cy}\right|-l\log \left|\mathrm{\Θ}\right|\end{array}---\left(11\right)$

该模型将信号重建、稀疏先验学习和噪声估计统一到一个框架下。

步骤四、为提升算法效率，引入如下的近似关系，

${(I\⊗{\mathrm{\Σ}}_n+{\mathrm{\Φ\Σ}}_y{\mathrm{\Φ}}^T)}^{-1}\≈{\mathrm{\Σ}}_{cy}^{-1}\⊗{({\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T)}^{-1}.---\left(12\right)$

基于关系(12)，采用坐标下降法将式子(11)分解为若干个子问题进行迭代求解，每个子问题中仅优化一个变量而固定剩余的其他变量。具体步骤如下：

①初始化λ，γ，κ为对应长度的全1向量，Σ_cy＝I，计数变量t＝0；

②学习参考协方差矩阵Θ。定义关于观测值矩阵F的权值矩阵M

M_ij为M的i行j列的元素，表示空间中以第i个像素为中心，大小为k＝3的邻域窗口中的所有光谱的观测值。||·||_F表示弗罗贝尼乌斯范数，σ＝0.7。参考协方差矩阵Θ＝(D-M)^-1D为对角阵，D_ii＝∑_jM_ij。

③固定λ和γ，根据式子(11)得到关于Y的子问题，如下

${\min }_y\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(14\right)$

基于近似关系(12)，求解得到Y的更新规则如下，

Y＝Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1F(15)

④固定Y，λ，κ和Σ_cy，利用近似关系(12)得到关于γ的子问题，如下

${\min }_{\mathrm{\γ}}\underset{i=1}{\overset{n_b}{\Σ}}\frac{Y_{i.}{\Σ}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+n_p\log |{\Σ}_n+A\mathrm{\Ψ}{\Σ}_{ry}{\mathrm{\Ψ}}^T{A}^T|+\underset{i=1}{\overset{n_b}{\Σ}}{\mathrm{\κ}}_i{\mathrm{\γ}}_i---\left(16\right)$

其中，Y_i.表示Y的第i行，γ_i为γ的第i个元素，求解得到如下的更新形式：

${\mathrm{\γ}}_i=(\sqrt{4{\mathrm{\κ}}_i(Y_{i.}{\mathrm{\Σ}}_{cy}^{-1}{Y}_{i.}^T+n_p{\mathrm{\α}}_i)+n_p^2}-n_p)/\left(2{\mathrm{\κ}}_i\right)---\left(17\right)$

其中，α＝diag[Σ_ry-Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1AΨΣ_ry]，与之前不同，此处diag(·)表示取矩阵对角线元素组成向量，α_i为α的第i个元素。

⑤固定Y和γ，利用近似关系(12)得到Σ_cy的子问题

${\min }_{{\mathrm{\Σ}}_{cy}}\underset{i=1}{\overset{n_b}{\Σ}}\frac{Y_{i.}{\Σ}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+\mathrm{\μ}log\left|{\mathrm{\Σ}}_{cy}\right|+tr\left({\mathrm{\Θ\Σ}}_{cy}^{-1}\right)---\left(18\right)$

μ＝m_b+n_p+l+1，求解得到Σ_cy的更新形式，如下：

${\mathrm{\Σ}}_{cy}=(Y^T{\mathrm{\Σ}}_{ry}^{-1}Y+\mathrm{\Θ}+\mathrm{\η}I)/\mathrm{\μ}---\left(19\right)$

为提升噪声鲁棒性，令 $\mathrm{\μ}=\left|\right|Y^T{\mathrm{\Σ}}_{ry}^{-1}Y+\mathrm{\Θ}+\mathrm{\η}I|{|}_F.$

⑥固定Y和γ，利用近似关系(12)得到关于λ的优化子问题，如下

${\min }_{\mathrm{\λ}}\left|\right|A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2+n_{p}log|{\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T|---\left(20\right)$

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

${\mathrm{\λ}}_i=\sqrt{\left(Q_{.i}^T{Q}_{.i}\right)/\left(n_p{\mathrm{\υ}}_i\right)}---\left(21\right)$

其中，λ_i为λ的第i个元素，Q＝AΨY-F，Q_.i表示Q的第i列，υ_i为向量υ＝diag[(Σ_n+AΨΣ_ryΨ^TA^T)^-1]的第i个元素，diag(·)运算和④步相同。

⑦固定γ，得到关于κ的优化子问题，如下

${\min }_{\mathrm{\κ}}\underset{i=1}{\overset{n_b}{\Σ}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\κ}}_i)---\left(22\right)$

κ_i为κ的第i个元素。求解得到如下的更新形式

κ_i＝2/γ_i(23)

⑧假设上一次迭代重建得到的稀疏信号为Y′，最新重建的稀疏信号为Y，说计算更新前后的差异，η＝||Y′-Y||_F/||Y′||_F，计数器t加1。如果计数器t≤200并且更新差异η≥10^-4，则循环执行步骤③至⑧；否则，退出循环。

假设最终得到最优估计的Y_rec，则重建高光谱图像，X_rec＝ΨY_rec。

本发明的有益效果是：该方法随机采样每个像素光谱的少量线性观测值作为压缩数据，通过流形结构化稀疏先验，同时刻画高光谱图像稀疏化后光谱维中的稀疏性和空间维中的流形结构；通过隐变量贝叶斯模型，将信号重建，稀疏先验学习以及噪声估计统一到一个正则化回归模型进行优化求解。学习得到的稀疏先验既能充分地刻画高光谱图像的三维结构，又具有较强的噪声鲁棒性。利用该稀疏先验，实现了高光谱图像的高精度重建。在真实的高光谱遥感数据集Urban上的实验结果表明，当在压缩数据中加入高斯白噪声使得压缩数据信噪比为15db，采样率为0.09时，获得了23db的峰值信噪比。

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

具体实施方式

本发明基于流形结构化稀疏先验的高光谱图像压缩感知方法具体包括以下步骤：

针对包含n_b个波段，每个波段包含n_r行和n_c列的高光谱图像，将每一个波段拉伸成为一个行向量，重新组成一个二维矩阵，其中，X的每一列对应每个像素的光谱；每一行对应每个波段的所有像素值。将行和列分别称为空间维和光谱维。本发明主要包含以下四个步骤：

1、获取压缩数据。

压缩过程中，利用列归一化的高斯随机观测矩阵随机采样高光谱图像X的光谱维，获得压缩数据m_b为压缩后波段长度。

F＝AX+N(1)其中，表示采样中的噪声。采样率ρ定义为ρ＝m_b/n_b。

2、建立基于流形结构化稀疏先验的压缩感知模型。

利用Haar小波基对高光谱图像的每个光谱进行稀疏化，如X＝ΨY，Ψ为小波基，Y为列稀疏的系数矩阵。因此，模型(1)可表示为F＝AΨY+N。假设采样过程中噪声N服从的矩阵正太分布，I为对应大小的单位矩阵(下同)。因此，模型(1)对应的似然函数可以定义为

$p\left(F\right|Y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{n}_p/2}|{\mathrm{\Σ}}_n{|}^{n_p/2}}---\left(2\right)$

其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵。表示Q矩阵的加权迹范数。

除过列稀疏性，高光谱图像空间像素之间的相似性使得Y中不同列的稀疏信号位于一个结构未知的流形结构上。为了充分描述Y的特性，假设Y服从如下的矩阵正太分布

$p-\left(Y\right|{\mathrm{\Σ}}_{ry},{\mathrm{\Σ}}_{cy})=\frac{\exp \{-\frac{1}{2}tr\left({\mathrm{\Σ}}_{cy}^{-1}{Y}^T{\mathrm{\Σ}}_{ry}^{-1}Y\right)\}}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}\left|{\mathrm{\Σ}}_{ry}{|}^{n_p/2}\right|{\mathrm{\Σ}}_{cy}{|}^{n_b/2}},---\left(3\right)$

为描述Y中列信号的稀疏性，令Σ_y＝diag(γ)表示以γ为对角线元素的对角矩阵，此外，令假设γ服从如下的伽玛分布

$p\left(\mathrm{\γ}\right|\mathrm{\κ})=\underset{i=1}{\overset{n_b}{\Π}}Gamma(1,\frac{2}{{\mathrm{\κ}}_i})=\underset{i=1}{\overset{n_b}{\Π}}\frac{{\mathrm{\κ}}_i}{2}\exp (-\frac{{\mathrm{\κ}}_i{\mathrm{\γ}}_i}{2})---\left(4\right)$

式子(3)中，Σ_cy描述Y中不同列信号之间的相关性，因此，式子(3)可以隐式地表示Y中不同的稀疏信号之间存在的流形结构。为了更加灵活地学习Σ_cy，进一步假设Σ_cy服从如下的反威沙特分布

其中，l是给定的常量，表示自由度，Γ_np是多变量伽玛函数，为参考协方差矩阵。该先验通过最小化Σ_cy和Θ之间的布雷格曼散度，使得Σ_cy趋近于Θ，从而减轻了Σ_cy学习过程中的过拟合问题。

3、建立正则化的回归模型。

为使得提出的流行结构化稀疏先验能够更好地匹配图像分布并具有较强的噪声鲁棒性，本发明提出了一种隐变量贝叶斯模型对噪声参数λ和先验参数γ，κ，Σ_cy和Θ等进行估计。令f＝vec(F)，y＝vec(Y)，n＝vec(N)和vec(Q)表示将矩阵Q拉成列向量，表示克罗内克积，则模型(2)等价于

$p\left(f\right|y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||f-\mathrm{\Φ}y|{|}_{I\⊗{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{m}_p/2}|I\⊗{\mathrm{\Σ}}_n{|}^{1/2}}---\left(6\right)$

同样，模型(3)中关于Y的先验等价于

$p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})=\frac{\exp (-\frac{1}{2}{y}^T{\mathrm{\Σ}}_y^{-1}y)}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}|{\mathrm{\Σ}}_y{|}^{1/2}},{\mathrm{\Σ}}_y={\mathrm{\Σ}}_{cy}\⊗{\mathrm{\Σ}}_{ry}.---\left(7\right)$

根据公式(6)、(7)，所有的未知参数可以通过求解如下的优化问题得到

$\underset{\mathrm{\λ},\mathrm{\γ}\≥0,\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}}{\max }p(\mathrm{\λ},\mathrm{\γ},\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}|f)\∝\∫p\left(f\right|y,\mathrm{\λ})p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})p\left(\mathrm{\γ}\right|\mathrm{\κ})p\left({\mathrm{\Σ}}_{cy}\right|\mathrm{\Θ},l)dy---\left(8\right)$

通过积分，并引入-2log运算，容易得知式子(8)，等价于最小化如下的式子

其中，tr(·)表示迹范数，对式子(9)的第一项做如下变形

$f^T{\mathrm{\Σ}}_{by}^{-1}f=\underset{y}{\min }\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(10\right)$

将式子(10)带入到式子(9)中，得到如下等价于式子(8)的正则化回归模型

$\begin{array}{c}\underset{y,\mathrm{\λ}\≥0,\mathrm{\γ}\≥0,\mathrm{\κ},{\Σ}_{cy},\mathrm{\Θ}}{\min }\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\Σ}_n}^2+y^T{\Σ}_y^{-1}y+\log |{\Σ}_{by}|+\underset{i=1}{\overset{n_b}{\Σ}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\γ}}_i)\\ +tr\left(\mathrm{\Θ}{\Σ}_{cy}^{-1}\right)+(n_p+l+1)\log \left|{\Σ}_{cy}\right|-l\log \left|\mathrm{\Θ}\right|\end{array}---\left(11\right)$

该模型将信号重建、稀疏先验学习和噪声估计统一到一个框架下。一方面，学习的稀疏先验会根据估计的噪声进行调整；另一方面，基于学习的稀疏先验，噪声估计更加准确。因此，该模型能够从噪声污染的观测值中精确地重建高光谱图像。

4、模型求解。

为提升算法效率，引入如下的近似关系，

${(I\⊗{\mathrm{\Σ}}_n+{\mathrm{\Φ\Σ}}_y{\mathrm{\Φ}}^T)}^{-1}\≈{\mathrm{\Σ}}_{cy}^{-1}\⊗{({\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T)}^{-1}.---\left(12\right)$

基于关系(12)，采用坐标下降法将式子(11)分解为若干个子问题进行迭代求解，每个子问题中仅优化一个变量而固定剩余的其他变量。具体步骤如下：

⑨初始化λ，γ，κ为对应长度的全1向量，Σ_cy＝I，计数变量t＝0；

⑩学习参考协方差矩阵Θ，首先，定义关于观测值矩阵F的权值矩阵M

M_ij为M的i行j列的元素，表示空间中以第i个像素为中心，大小为k＝3的邻域窗口中的所有光谱的观测值。||·||_F表示弗罗贝尼乌斯范数，σ＝0.7。参考协方差矩阵Θ＝(D-M)^-1D为对角阵，D_ii＝∑_jM_ij。

固定λ和γ，根据式子(11)得到关于Y的子问题，如下

${\min }_y\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(14\right)$

基于近似关系(12)，求解得到Y的更新规则如下，

Y＝Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1F(15)

固定Y，λ，κ和Σ_cy，利用近似关系(12)得到关于γ的子问题，如下

${\min }_{\mathrm{\γ}}\underset{i=1}{\overset{n_b}{\Σ}}\frac{Y_{i.}{\Σ}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+n_p\log |{\Σ}_n+A\mathrm{\Ψ}{\Σ}_{ry}{\mathrm{\Ψ}}^T{A}^T|+\underset{i=1}{\overset{n_b}{\Σ}}{\mathrm{\κ}}_i{\mathrm{\γ}}_i---\left(16\right)$

其中，Y_i表示Y的第i行，γ_i为γ的第i个元素，求解得到如下的更新形式：

${\mathrm{\γ}}_i=(\sqrt{4{\mathrm{\κ}}_i(Y_{i.}{\mathrm{\Σ}}_{cy}^{-1}{Y}_{i.}^T+n_p{\mathrm{\α}}_i)+n_p^2}-n_p)/\left(2{\mathrm{\κ}}_i\right)---\left(17\right)$

其中，α＝diag[Σ_ry-Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1AΨΣ_ry]，与之前不同，此处diag(·)表示取矩阵对角线元素组成向量，α_i为α的第i个元素。

固定Y和γ，利用近似关系(12)得到Σ_cy的子问题

${\min }_{{\mathrm{\Σ}}_{cy}}\underset{i=1}{\overset{n_b}{\Σ}}\frac{Y_{i.}{\Σ}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+\mathrm{\μ}log\left|{\mathrm{\Σ}}_{cy}\right|+tr\left({\mathrm{\Θ\Σ}}_{cy}^{-1}\right)---\left(18\right)$

μ＝m_b+n_p+l+1，求解得到Σ_cy的更新形式，如下：

${\mathrm{\Σ}}_{cy}=(Y^T{\mathrm{\Σ}}_{ry}^{-1}Y+\mathrm{\Θ}+\mathrm{\η}I)/\mathrm{\μ}---\left(19\right)$

为提升噪声鲁棒性，本发明令

固定Y和γ，利用近似关系(12)得到关于λ的优化子问题，如下

${\min }_{\mathrm{\λ}}\left|\right|A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2+n_{p}log|{\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T|---\left(20\right)$

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

${\mathrm{\λ}}_i=\sqrt{\left(Q_{.i}^T{Q}_{.i}\right)/\left(n_p{\mathrm{\υ}}_i\right)}---\left(21\right)$

其中，λ_i为λ的第i个元素，Q＝AΨY-F，Q_.i表示Q的第i列，υ_i为向量υ＝diag[(Σ_n+AΨΣ_ryΨ^TA^T)^-1]的第i个元素，diag(·)运算和④步相同。

固定γ，得到关于κ的优化子问题，如下

${\min }_{\mathrm{\κ}}\underset{i=1}{\overset{n_b}{\Σ}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\κ}}_i)---\left(22\right)$

κ_i为κ的第i个元素。求解得到如下的更新形式

κ_i＝2/γ_i(23)

假设上一次迭代重建得到的稀疏信号为Y′，最新重建的稀疏信号为Y，说计算更新前后的差异，η＝||Y′-Y||_F/||Y′||_F，计数器t加1。如果计数器t≤200并且更新差异η≥10^-4，则循环执行③至⑧；否则，退出循环。

假设最终得到最优估计的Y_rec，则可以重建高光谱图像，X_rec＝ΨY_rec。

Claims (1)

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

步骤一、针对包含n_b个波段，每个波段包含n_r行和n_c列的高光谱图像，将每一个波段拉伸成为一个行向量，重新组成一个二维矩阵，(n_p＝n_r×n_c)；其中，X的每一列对应每个像素的光谱；每一行对应每个波段的所有像素值；将行和列分别称为空间维和光谱维；

获取压缩数据过程中，利用列归一化的高斯随机观测矩阵随机采样高光谱图像X的光谱维，获得压缩数据m_b为压缩后波段长度；

F＝AX+N(1)

其中，表示采样中的噪声；采样率ρ定义为ρ＝m_b/n_b；

步骤二、利用Haar小波基对高光谱图像的每个光谱进行稀疏化，如X＝ΨY，Ψ为小波基，Y为列稀疏的系数矩阵，模型(1)表示为F＝AΨY+N；假设采样过程中噪声N服从的矩阵正太分布，I为对应大小的单位矩阵，模型(1)对应的似然函数定义为

$p\left(F\right|Y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{n}_p/2}|{\mathrm{\Σ}}_n{|}^{n_p/2}}---\left(2\right)$

其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵；表示Q矩阵的加权迹范数；

除过列稀疏性，高光谱图像空间像素之间的相似性使得Y中不同列的稀疏信号位于一个结构未知的流形结构上；为了充分描述Y的特性，假设Y服从如下的矩阵正太分布

$p\left(Y\right|{\mathrm{\Σ}}_{ry},{\mathrm{\Σ}}_{cy})=\frac{\exp \{-\frac{1}{2}tr\left({\mathrm{\Σ}}_{cy}^{-1}{Y}^T{\mathrm{\Σ}}_{ry}^{-1}Y\right)\}}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}\left|{\mathrm{\Σ}}_{ry}{|}^{n_p/2}\right|{\mathrm{\Σ}}_{cy}{|}^{n_b/2}},---\left(3\right)$

为描述Y中列信号的稀疏性，令Σ_y＝diag(γ)表示以γ为对角线元素的对角矩阵， $\mathrm{\γ}={\[{\mathrm{\γ}}_1,...,{\mathrm{\γ}}_{n_b}\]}^T;$ 令 $\mathrm{\κ}={\[{\mathrm{\κ}}_1,...,{\mathrm{\κ}}_{n_b}\]}^T,$ 假设γ服从如下的伽玛分布

$p\left(\mathrm{\γ}\right|\mathrm{\κ})=\underset{i=1}{\overset{n_b}{\Π}}Gamma(1,\frac{2}{{\mathrm{\κ}}_i})=\underset{i=1}{\overset{n_b}{\Π}}\frac{{\mathrm{\κ}}_i}{2}\exp (-\frac{{\mathrm{\κ}}_i{\mathrm{\γ}}_i}{2})---\left(4\right)$

式子(3)中，Σ_cy描述Y中不同列信号之间的相关性，因此，式子(3)隐式地表示Y中不同的稀疏信号之间存在的流形结构；为了更加灵活地学习Σ_cy，进一步假设Σ_cy服从如下的反威沙特分布

其中，l是给定的常量，表示自由度，是多变量伽玛函数，为参考协方差矩阵；该先验通过最小化Σ_cy和Θ之间的布雷格曼散度，使得Σ_cy趋近于Θ，从而减轻了Σ_cy学习过程中的过拟合问题；

步骤三、为使得流行结构化稀疏先验能够更好地匹配图像分布并具有较强的噪声鲁棒性，通过隐变量贝叶斯模型对噪声参数λ和先验参数γ，κ，Σ_cy和Θ进行估计；令f＝vec(F)，y＝vec(Y)，n＝vec(N)和vec(Q)表示将矩阵Q拉成列向量，表示克罗内克积，则模型(2)等价于

$p\left(f\right|y,\mathrm{\λ})=\frac{\exp \{-\frac{1}{2}||f-\mathrm{\Φ}y|{|}_{I\⊗{\mathrm{\Σ}}_n}^2\}}{{\left(2\mathrm{\π}\right)}^{m_b{m}_p/2}|I\⊗{\mathrm{\Σ}}_n{|}^{1/2}}---\left(6\right)$

同样，模型(3)中关于Y的先验等价于

$p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})=\frac{\exp (-\frac{1}{2}{y}^T{\mathrm{\Σ}}_y^{-1}y)}{{\left(2\mathrm{\π}\right)}^{n_b{n}_p/2}|{\mathrm{\Σ}}_y{|}^{1/2}},{\mathrm{\Σ}}_y={\mathrm{\Σ}}_{cy}\⊗{\mathrm{\Σ}}_{ry}.---\left(7\right)$

根据公式(6)、公式(7)，所有的未知参数通过求解如下的优化问题得到

$\underset{\mathrm{\λ},\mathrm{\γ}\≥0,\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}}{max}p(\mathrm{\λ},\mathrm{\γ},\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}|f)\∝\∫p\left(f\right|y,\mathrm{\λ})p\left(y\right|\mathrm{\γ},{\mathrm{\Σ}}_{cy})p\left(\mathrm{\γ}\right|\mathrm{\κ})p\left({\mathrm{\Σ}}_{cy}\right|\mathrm{\Θ},l)dy---\left(8\right)$

通过积分，并引入-2log运算，容易得知式子(8)，等价于最小化如下的式子

其中，tr(·)表示迹范数，对式子(9)的第一项做如下变形

$f^T{\mathrm{\Σ}}_{by}^{-1}f=\underset{y}{min}\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(10\right)$

将式子(10)带入到式子(9)中，得到如下等价于式子(8)的正则化回归模型

$\begin{array}{c}\underset{\mathrm{\λ},\mathrm{\γ}\≥0,\mathrm{\κ},{\mathrm{\Σ}}_{cy},\mathrm{\Θ}}{max}\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y+\log |\left|{\mathrm{\Σ}}_{by}\right|+\underset{i=1}{\overset{n_b}{\mathrm{\Σ}}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\γ}}_i)\\ +tr\left({\mathrm{\Θ\Σ}}_{cy}^{-1}\right)+(n_p+l+1)\log \left|{\mathrm{\Σ}}_{cy}\right|-l\log \left|\mathrm{\Θ}\right|\end{array}---\left(11\right)$

该模型将信号重建、稀疏先验学习和噪声估计统一到一个框架下；

步骤四、为提升算法效率，引入如下的近似关系，

${(I\⊗{\mathrm{\Σ}}_n+{\mathrm{\Φ\Σ}}_y{\mathrm{\Φ}}^T)}^{-1}\≈{\mathrm{\Σ}}_{cy}^{-1}\⊗{({\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T)}^{-1}.---\left(12\right)$

基于关系(12)，采用坐标下降法将式子(11)分解为若干个子问题进行迭代求解，每个子问题中仅优化一个变量而固定剩余的其他变量；具体步骤如下：

①初始化λ，γ，κ为对应长度的全1向量，Σ_cy＝I，计数变量t＝0；

②学习参考协方差矩阵Θ；定义关于观测值矩阵F的权值矩阵M

M_ij为M的i行j列的元素，表示空间中以第i个像素为中心，大小为k＝3的邻域窗口中的所有光谱的观测值；||·||_F表示弗罗贝尼乌斯范数，σ＝0.7；参考协方差矩阵Θ＝(D-M)^-1，D为对角阵，D_ii＝∑_jM_ij；

③固定λ和γ，根据式子(11)得到关于Y的子问题，如下

${\min }_y\left|\right|\mathrm{\Φ}y-f|{|}_{I\⊗{\mathrm{\Σ}}_n}^2+y^T{\mathrm{\Σ}}_y^{-1}y---\left(14\right)$

基于近似关系(12)，求解得到Y的更新规则如下，

Y＝Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1F(15)

④固定Y，λ，κ和Σ_cy，利用近似关系(12)得到关于γ的子问题，如下

${\min }_{\mathrm{\γ}}\underset{i=1}{\overset{n_b}{\mathrm{\Σ}}}\frac{Y_{i.}{\mathrm{\Σ}}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+n_{p}log|{\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T|+\underset{i=1}{\overset{n_b}{\mathrm{\Σ}}}{\mathrm{\κ}}_i{\mathrm{\γ}}_i---\left(16\right)$

其中，Y_i.表示Y的第i行，γ_i为γ的第i个元素，求解得到如下的更新形式：

${\mathrm{\γ}}_i=(\sqrt{4{\mathrm{\κ}}_i(Y_{i.}{\mathrm{\Σ}}_{cy}^{-1}{Y}_{i.}^T+n_p{\mathrm{\α}}_i)+n_p^2}-n_p)/\left(2{\mathrm{\κ}}_i\right)---\left(17\right)$

其中，α＝diag[Σ_ry-Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^-1AΨΣ_ry]，与之前不同，此处diag(·)表示取矩阵对角线元素组成向量，α_i为α的第i个元素；

⑤固定Y和γ，利用近似关系(12)得到Σ_cy的子问题

${\min }_{{\mathrm{\Σ}}_{cy}}\underset{i=1}{\overset{n_b}{\mathrm{\Σ}}}\frac{Y_{i.}{\mathrm{\Σ}}_{cy}^{-1}{Y}_{i.}^T}{{\mathrm{\γ}}_i}+\mathrm{\μ}log\left|{\mathrm{\Σ}}_{cy}\right|+tr\left({\mathrm{\Θ\Σ}}_{cy}^{-1}\right)---\left(18\right)$

μ＝m_b+n_p+l+1，求解得到Σ_cy的更新形式，如下：

${\mathrm{\Σ}}_{cy}=(Y^T{\mathrm{\Σ}}_{ry}^{-1}Y+\mathrm{\Θ}+\mathrm{\η}I)/\mathrm{\μ}---\left(19\right)$

为提升噪声鲁棒性，令 $\mathrm{\μ}=\left|\right|Y^T{\mathrm{\Σ}}_{ry}^{-1}Y+\mathrm{\Θ}+\mathrm{\η}I|{|}_F;$

⑥固定Y和γ，利用近似关系(12)得到关于λ的优化子问题，如下

${\min }_{\mathrm{\λ}}\left|\right|A\mathrm{\Ψ}Y-F|{|}_{{\mathrm{\Σ}}_n}^2+n_{p}log|{\mathrm{\Σ}}_n+{\mathrm{A\Ψ\Σ}}_{ry}{\mathrm{\Ψ}}^T{A}^T|---\left(20\right)$

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

${\mathrm{\λ}}_i=\sqrt{\left(Q_{.i}^T{Q}_{.i}\right)/\left(n_p{\mathrm{\υ}}_i\right)}---\left(21\right)$

其中，λ_i为λ的第i个元素，Q＝AΨY-F，Q_.i表示Q的第i列，υ_i为向量υ＝diag[(Σ_n+AΨΣ_ryΨ^TA^T)^-1]的第i个元素，diag(·)运算和④步相同；

⑦固定γ，得到关于κ的优化子问题，如下

${\min }_{\mathrm{\κ}}\underset{i=1}{\overset{n_b}{\Σ}}({\mathrm{\κ}}_i{\mathrm{\γ}}_i-2{\mathrm{log\κ}}_i)---\left(22\right)$

κ_i为κ的第i个元素；求解得到如下的更新形式

κ_i＝2/γ_i(23)

⑧假设上一次迭代重建得到的稀疏信号为Y′，最新重建的稀疏信号为Y，说计算更新前后的差异，η＝||Y′-Y||_F/||Y′||_F，计数器t加1；如果计数器t≤200并且更新差异η≥10^-4，则循环执行步骤③至⑧；否则，退出循环；

假设最终得到最优估计的Y_rec，则重建高光谱图像，X_rec＝ΨY_rec。