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

摘要

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

主权项

一种基于流形结构化稀疏先验的高光谱图像压缩感知方法，其特征在于包括以下步骤：步骤一、针对包含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,\lambda )=\frac{\exp \{-\frac{1}{2}||A\Psi Y-F|{|}_{{\Sigma }_n}^2\}}{{\left(2\pi \right)}^{m_b{n}_p/2}|{\Sigma }_n{|}^{n_p/2}}---\left(2\right)$其中，Σ_n＝diag(λ)表示以λ为对角线元素的对角矩阵；表示Q矩阵的加权迹范数；除过列稀疏性，高光谱图像空间像素之间的相似性使得Y中不同列的稀疏信号位于一个结构未知的流形结构上；为了充分描述Y的特性，假设Y服从如下的矩阵正太分布$p\left(Y\right|{\Sigma }_{ry},{\Sigma }_{cy})=\frac{\exp \{-\frac{1}{2}tr\left({\Sigma }_{cy}^{-1}{Y}^T{\Sigma }_{ry}^{-1}Y\right)\}}{{\left(2\pi \right)}^{n_b{n}_p/2}\left|{\Sigma }_{ry}{|}^{n_p/2}\right|{\Sigma }_{cy}{|}^{n_b/2}},---\left(3\right)$为描述Y中列信号的稀疏性，令Σ_y＝diag(γ)表示以γ为对角线元素的对角矩阵，$\gamma ={[{\gamma }_1,...,{\gamma }_{n_b}]}^T;$令$\kappa ={[{\kappa }_1,...,{\kappa }_{n_b}]}^T,$假设γ服从如下的伽玛分布$p\left(\gamma \right|\kappa )=\underset{i=1}{\overset{n_b}{\Pi }}Gamma(1,\frac{2}{{\kappa }_i})=\underset{i=1}{\overset{n_b}{\Pi }}\frac{{\kappa }_i}{2}\exp (-\frac{{\kappa }_i{\gamma }_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,\lambda )=\frac{\exp \{-\frac{1}{2}||f-\Phi y|{|}_{I\otimes {\Sigma }_n}^2\}}{{\left(2\pi \right)}^{m_b{m}_p/2}|I\otimes {\Sigma }_n{|}^{1/2}}---\left(6\right)$同样，模型(3)中关于Y的先验等价于$p\left(y\right|\gamma ,{\Sigma }_{cy})=\frac{\exp (-\frac{1}{2}{y}^T{\Sigma }_y^{-1}y)}{{\left(2\pi \right)}^{n_b{n}_p/2}|{\Sigma }_y{|}^{1/2}},{\Sigma }_y={\Sigma }_{cy}\otimes {\Sigma }_{ry}.---\left(7\right)$根据公式(6)、公式(7)，所有的未知参数通过求解如下的优化问题得到$\underset{\lambda ,\gamma \ge 0,\kappa ,{\Sigma }_{cy},\Theta }{max}p(\lambda ,\gamma ,\kappa ,{\Sigma }_{cy},\Theta |f)\propto \int p\left(f\right|y,\lambda )p\left(y\right|\gamma ,{\Sigma }_{cy})p\left(\gamma \right|\kappa )p\left({\Sigma }_{cy}\right|\Theta ,l)dy---\left(8\right)$通过积分，并引入‑2log运算，容易得知式子(8)，等价于最小化如下的式子其中，tr(·)表示迹范数，对式子(9)的第一项做如下变形$f^T{\Sigma }_{by}^{-1}f=\underset{y}{min}\left|\right|\Phi y-f|{|}_{I\otimes {\Sigma }_n}^2+y^T{\Sigma }_y^{-1}y---\left(10\right)$将式子(10)带入到式子(9)中，得到如下等价于式子(8)的正则化回归模型$\begin{array}{c}\underset{\lambda ,\gamma \ge 0,\kappa ,{\Sigma }_{cy},\Theta }{max}\left|\right|\Phi y-f|{|}_{I\otimes {\Sigma }_n}^2+y^T{\Sigma }_y^{-1}y+\log |\left|{\Sigma }_{by}\right|+\underset{i=1}{\overset{n_b}{\Sigma }}({\kappa }_i{\gamma }_i-2{\mathrm{log\gamma }}_i)\\ +tr\left({\mathrm{\Theta \Sigma }}_{cy}^{-1}\right)+(n_p+l+1)\log \left|{\Sigma }_{cy}\right|-l\log \left|\Theta \right|\end{array}---\left(11\right)$该模型将信号重建、稀疏先验学习和噪声估计统一到一个框架下；步骤四、为提升算法效率，引入如下的近似关系，${(I\otimes {\Sigma }_n+{\mathrm{\Phi \Sigma }}_y{\Phi }^T)}^{-1}\approx {\Sigma }_{cy}^{-1}\otimes {({\Sigma }_n+{\mathrm{A\Psi \Sigma }}_{ry}{\Psi }^T{A}^T)}^{-1}.---\left(12\right)$基于关系(12)，采用坐标下降法将式子(11)分解为若干个子问题进行迭代求解，每个子问题中仅优化一个变量而固定剩余的其他变量；具体步骤如下：①初始化λ，γ，κ为对应长度的全1向量，Σ_cy＝I，计数变量t＝0；②学习参考协方差矩阵Θ；定义关于观测值矩阵F的权值矩阵MM_ij为M的i行j列的元素，表示空间中以第i个像素为中心，大小为k＝3的邻域窗口中的所有光谱的观测值；||·||_F表示弗罗贝尼乌斯范数，σ＝0.7；参考协方差矩阵Θ＝(D‑M)^‑1，D为对角阵，D_ii＝∑_jM_ij；③固定λ和γ，根据式子(11)得到关于Y的子问题，如下${\min }_y\left|\right|\Phi y-f|{|}_{I\otimes {\Sigma }_n}^2+y^T{\Sigma }_y^{-1}y---\left(14\right)$基于近似关系(12)，求解得到Y的更新规则如下，Y＝Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^‑1F         (15)④固定Y，λ，κ和Σ_cy，利用近似关系(12)得到关于γ的子问题，如下${\min }_{\gamma }\underset{i=1}{\overset{n_b}{\Sigma }}\frac{Y_{i.}{\Sigma }_{cy}^{-1}{Y}_{i.}^T}{{\gamma }_i}+n_{p}log|{\Sigma }_n+{\mathrm{A\Psi \Sigma }}_{ry}{\Psi }^T{A}^T|+\underset{i=1}{\overset{n_b}{\Sigma }}{\kappa }_i{\gamma }_i---\left(16\right)$其中，Y_i.表示Y的第i行，γ_i为γ的第i个元素，求解得到如下的更新形式：${\gamma }_i=(\sqrt{4{\kappa }_i(Y_{i.}{\Sigma }_{cy}^{-1}{Y}_{i.}^T+n_p{\alpha }_i)+n_p^2}-n_p)/\left(2{\kappa }_i\right)---\left(17\right)$其中，α＝diag[Σ_ry‑Σ_ryΨ^TA^T(Σ_n+AΨΣ_ryΨ^TA^T)^‑1AΨΣ_ry]，与之前不同，此处diag(·)表示取矩阵对角线元素组成向量，α_i为α的第i个元素；⑤固定Y和γ，利用近似关系(12)得到Σ_cy的子问题${\min }_{{\Sigma }_{cy}}\underset{i=1}{\overset{n_b}{\Sigma }}\frac{Y_{i.}{\Sigma }_{cy}^{-1}{Y}_{i.}^T}{{\gamma }_i}+\mu log\left|{\Sigma }_{cy}\right|+tr\left({\mathrm{\Theta \Sigma }}_{cy}^{-1}\right)---\left(18\right)$μ＝m_b+n_p+l+1，求解得到Σ_cy的更新形式，如下：${\Sigma }_{cy}=(Y^T{\Sigma }_{ry}^{-1}Y+\Theta +\eta I)/\mu ---\left(19\right)$为提升噪声鲁棒性，令$\mu =\left|\right|Y^T{\Sigma }_{ry}^{-1}Y+\Theta +\eta I|{|}_F;$⑥固定Y和γ，利用近似关系(12)得到关于λ的优化子问题，如下${\min }_{\lambda }\left|\right|A\Psi Y-F|{|}_{{\Sigma }_n}^2+n_{p}log|{\Sigma }_n+{\mathrm{A\Psi \Sigma }}_{ry}{\Psi }^T{A}^T|---\left(20\right)$求解得到如下的更新形式：${\lambda }_i=\sqrt{\left(Q_{.i}^T{Q}_{.i}\right)/\left(n_p{\upsilon }_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 }_{\kappa }\underset{i=1}{\overset{n_b}{\Sigma }}({\kappa }_i{\gamma }_i-2{\mathrm{log\kappa }}_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。