基于重加权拉普拉斯稀疏先验的高光谱图像压缩感知方法

摘要

本发明公开了一种基于重加权拉普拉斯稀疏先验的高光谱图像压缩感知方法，用于解决现有高光谱图像压缩感知方法重建精度低的技术问题。技术方案是随机采集每个像素光谱的少量线性观测作为压缩数据，建立基于重加权拉普拉斯稀疏先验的压缩感知模型和稀疏正则化的回归模型，对所建立的模型求解。由于随机采集少量线性观测作为压缩数据，减少了图像采集过程中的资源消耗。重加权拉普拉斯稀疏先验准确刻画了高光谱图像中的强稀疏性，克服了传统拉普拉斯稀疏先验对非零元素的非均匀约束，提高了高光谱图像的重建精度。经测试，当采样率为0.15，压缩数据中存在信噪比为10db的强噪声时，本发明的峰值信噪比相对于背景技术方法提升4db以上。

说明书

技术领域

本发明涉及一种高光谱图像压缩感知方法，特别是涉及一种基于重加权拉普拉斯 稀疏先验的高光谱图像压缩感知方法。

背景技术

高光谱图像的光谱信息有助于遥感地物的检测定位与分类识别，然而，高光谱图 像的巨大数据量对图像采集、传输和处理中的软硬资源提出了严格的要求，制约了高 光谱图像的应用。因此，高光谱图像压缩算法是高光谱领域的研究热点之一。目前， 大量普通图像的压缩方法已成功推广应用到高光谱图像上。然而，这类方法只能压缩 已经获取的图像，无法减少成像过程中的巨大资源需求。近年来，压缩感知成像理论 证明，仅需要从场景中采集少量的线性观测值便可以对原始场景的图像进行高精度的 重建。ChengBo Li等人在文献“A compressive sensing and unmixing scheme for  hyperspectral data processing,IEEE Transactions on Image Processing,2012, 21(3):1200–1210”中，针对高光谱图像每个波段使用单像素相机采集少量线性观测值作 为压缩数据。重建过程中，引入少量的端元光谱，结合线性混合模型，在全变差梯度 稀疏的约束下，重建空间连续的丰度值矩阵。最终，利用线性混合模型重建高光谱图 像。然而，该方法使用的稀疏约束忽略了稀疏信号内非零元素之间的关系，重建精度 受限；其次，算法性能严重依赖端元的选择，不当的端元选择可导致重建彻底失败； 此外，算法参数需要调节，适应性差。

发明内容

为了克服现有高光谱图像压缩感知方法重建精度低的不足，本发明提供一种基 于重加权拉普拉斯稀疏先验的高光谱图像压缩感知方法。该方法随机采集每个像素光 谱的少量线性观测作为压缩数据，建立基于重加权拉普拉斯稀疏先验的压缩感知模型 和稀疏正则化的回归模型，对所建立的模型求解。由于随机采集每个像素光谱的少量 线性观测作为压缩数据，减少了图像采集过程中的资源消耗。重建过程中，引入的重 加权拉普拉斯稀疏先验准确刻画了高光谱图像中的强稀疏性，克服了传统拉普拉斯稀 疏先验对非零元素的非均匀约束。通过经验贝叶斯推理，构建噪声鲁棒的非分离稀疏 先验约束，实现了高光谱图像的高精度重建。小波正交基的使用消除了对端元的依赖 性，所有的未知参数均可全自动估计。在真实的卫星图像Urban，Pavia University以 及Indiana数据集上的试验结果表明，当采样率为0.15，压缩数据中存在信噪比为10db 的强噪声时，本发明获得的峰值信噪比相对于背景技术压缩感知方法提升4db以上。

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

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

步骤二、采用列归一化的高斯随机观测矩阵采样高光谱图像X的光谱维， 获得压缩数据m_b表示压缩后的波段长度，ρ＝m_b/n_b为采样率。

G＝AX+N  (1)

其中，表示采样中的噪声。

步骤三、利用Haar小波基稀疏化高光谱图像的每一个光谱，X＝DY，D为小波 正交基，Y为列稀疏的小波系数矩阵。因此，(1)式表示为G＝ADY+N。假设采样过 程中噪声服从的矩阵正太分布，则对应的似然函数为

$p\left(G\right|Y,\lambda )=\frac{\exp \{-\frac{1}{2}{\left|\right|\mathrm{ADY}-G\left|\right|}_{{\Sigma }_n}^2\}}{{\left(2\pi \right)}^{m_b{n}_p/2}{\left|{\Sigma }_n\right|}^{n_p/2}},{\Sigma }_n=\mathrm{diag}\left(\lambda \right),---\left(2\right)$

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

对于稀疏信号Y，由于拉普拉斯分布不是(2)式中似然函数共轭先验，采用级联 分布的方式构造重加权拉普拉斯稀疏先验。首先假设Y服从如下分布

$p\left(Y\right|\gamma )=\frac{\exp \{-\frac{1}{2}{\left|\right|Y\left|\right|}_{{\Sigma }_y}^2\}}{{\left(2\pi \right)}^{n_b{n}_p/2}{\left|{\Sigma }_y\right|}^{n_p/2}},{\Sigma }_y=\mathrm{diag}\left(\gamma \right),---\left(3\right)$

其中，Σ_y＝diag(γ)表示以γ的元素为对角线元素的对角矩阵，控制Y中 每一行的稀疏度，γ_i＝0表示Y的第i行为0。假设则其中任一列y_i服从 高斯分布。假设超参数γ服从以下的伽马分布，

$p\left(\gamma \right|\kappa )=\underset{i=1}{\overset{n_b}{\Pi }}\mathrm{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)$

以上两级先验，等价于重加权拉普拉斯分布，因为对于y_i有

$p\left(y_i\right)=\int p\left(y_i\right|\gamma )p\left(\gamma \right|\kappa )\mathrm{d\gamma }=\frac{1}{2}\underset{j=1}{\overset{n_b}{\Pi }}{\kappa }_j^{1/2}\exp (-{\left|\right|{\mathrm{Ky}}_i\left|\right|}_1)---\left(5\right)$

其中，采用级联先验，λ，γ和κ均为待估计参数。

步骤四、由于λ，γ和κ未知，因此，根据经验贝叶斯框架，先基于压缩数据G利 用MAP估计未知参数λ，γ和κ，如下

$\{{\lambda }_{\mathrm{opt}},{\gamma }_{\mathrm{opt}},{\kappa }_{\mathrm{opt}}\}=\arg \underset{\lambda ,\gamma ,\kappa }{\max }p(\lambda ,\gamma ,\kappa |G)\propto \arg \underset{\lambda ,\gamma ,\kappa }{\max }\int p\left(G\right|Y,\lambda )p\left(Y\right|\gamma )p\left(\gamma \right|\kappa )\mathrm{dY}---\left(6\right)$

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

其中，tr(·)表示迹范数，Σ_by＝Σ_n+ADΣ_yD^TA^T，为代价函数。

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

$\mathrm{tr}\left(n_{n_p}^{-1}{G}^T{\Sigma }_{\mathrm{by}}^{-1}G\right)=\mathrm{tr}\left[\left(\frac{G^T}{\sqrt{n_p}}\right){\Sigma }_{\mathrm{by}}^{-1}\left(\frac{G}{\sqrt{n_p}}\right)\right]=\underset{Y}{\min }{\left|\right|\mathrm{ADY}-\frac{G}{\sqrt{n_p}}\left|\right|}_{{\Sigma }_n}^2+{\left|\right|Y\left|\right|}_{{\Sigma }_y}^2.---\left(8\right)$

然后，将(8)式带入到(7)式中，得

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

显然，而且能够证明，最小化(7)式再对稀疏信号Y进行 MAP估计，与直接最小化(10)式得到的λ，γ和κ相同，关于Y的解仅相差一个常 量因此，(10)式看作是关于稀疏信号Y的正则化回归模型，其中为 稀疏信号的非分离稀疏约束。该约束不能拆分成对于Y中每一行的独立约束，因此该 约束能同时约束稀疏信号中非零元素，潜在地考量非零元素之间的相关性。此外，Σ_by中包含了表征噪声强度的λ，因此，得到的稀疏约束随着估计的噪声强度自适应的变 化，具有噪声鲁棒性。

步骤五、已知压缩后数据G和采样矩阵A，采用坐标下降法最小化(10)式，每 次迭代中仅优化一个变量而固定剩余的其他变量。具体步骤如下：

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

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

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

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

$Y^{t+1}={\Sigma }_y{D}^T{A}^T{\Sigma }_{\mathrm{by}}^{-1}\frac{G}{\sqrt{n_p}}---\left(12\right)$

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

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

${\gamma }_i^{t+1}=\frac{\sqrt{4{\kappa }_i{n}_p{z}_i+n_p^2}-n_p}{2{\kappa }_i}---\left(14\right)$

其中，为γ^t+1的第i个元素，代表V^T+Y^t+1(Y^t+1)^T的对角线元素组成的向量，z_i为z的第i个元素；

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

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

${\lambda }^{t+1}=\sqrt{\mathrm{diag}\left({\mathrm{QQ}}^T\right)./\alpha }---\left(16\right)$

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

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

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

${\kappa }^{t+1}=\frac{2}{\gamma +d}---\left(18\right)$

前式中的加法和除法运算均作用在向量的每一个元素上，得到一个新向量，d＝10_-⁶的 引入是为了确保γ中出现0时，(18)式依然有意义。

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

$\eta =\frac{{\left|\right|Y^{t+1}.*\sqrt{n_p}\left|\right|}_F}{{\left|\right|Y^t.*\sqrt{n_p}\left|\right|}_F}---\left(19\right)$

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

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

$X_{\mathrm{rec}}=D(Y_{\mathrm{rec}}.*\sqrt{n_p})---\left(20\right).$

本发明的有益效果是：该方法随机采集每个像素光谱的少量线性观测作为压缩数 据，建立基于重加权拉普拉斯稀疏先验的压缩感知模型和稀疏正则化的回归模型，对 所建立的模型求解。由于随机采集每个像素光谱的少量线性观测作为压缩数据，减少 了图像采集过程中的资源消耗。重建过程中，引入的重加权拉普拉斯稀疏先验准确刻 画了高光谱图像中的强稀疏性，克服了传统拉普拉斯稀疏先验对非零元素的非均匀约 束。通过经验贝叶斯推理，构建噪声鲁棒的非分离稀疏先验约束，实现了高光谱图像 的高精度重建。小波正交基的使用消除了对端元的依赖性，所有的未知参数均可全自 动估计。在真实的卫星图像Urban，Pavia University以及Indiana数据集上的试验结果 表明，当采样率为0.15，压缩数据中存在信噪比为10db的强噪声时，本发明获得的峰 值信噪比相对于背景技术压缩感知方法提升4db以上。

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

具体实施方式

本发明基于重加权拉普拉斯稀疏先验的高光谱图像压缩感知方法具体步骤如下：

将包含n_b个波段，每个波段包含n_p个像素的高光谱图像的每一个波段拉伸成为一 个行向量，所有行向量组成一个二维矩阵，其中，X的每一列表示每一个 像素对应的光谱，该方向为光谱维；X的每一行对应一个波段的所有像素值，该方向 为空间维。本发明主要包含以下四个步骤，具体如下：

1、获取压缩数据。

采用列归一化的高斯随机观测矩阵采样高光谱图像X的光谱维，获得压 缩数据m_b表示压缩后的波段长度，ρ＝m_b/n_b为采样率。

G＝AX+N  (1) 其中，表示采样中的噪声。

2、建立基于重加权拉普拉斯稀疏先验的压缩感知模型。

本发明利用Haar小波基稀疏化高光谱图像的每一个光谱，X＝DY，D为小波正 交基，Y为列稀疏的小波系数矩阵。因此，(1)式可表示为G＝ADY+N。假设采样过 程中噪声服从的矩阵正太分布，则对应的似然函数为

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

对于稀疏信号Y，由于拉普拉斯分布不是(2)式中似然函数共轭先验，本发明采 用级联分布的方式构造重加权拉普拉斯稀疏先验。首先假设Y服从如下分布

其中，Σ_y＝diag(γ)表示以γ的元素为对角线元素的对角矩阵，控制Y中 每一行的稀疏度，γ_i＝0表示Y的第i行为0。假设则其中任一列y_i服从 高斯分布。假设超参数γ服从以下的伽马分布，

以上两级先验，等价于重加权拉普拉斯分布，因为对于y_i有

其中，在本发明中不直接使用重加权拉普拉斯先验，而是采 用级联先验，保证求解便利，而且λ，γ和κ均为待估计参数。

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

由于λ，γ和κ未知，无法采用最大后验估计(Maximum a posterior esitmation，MAP) 直接对稀疏信号Y进行估计。因此，本发明根据经验贝叶斯框架，先基于压缩数据G利 用MAP估计未知参数λ，γ和κ，如下

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

其中，tr(·)表示迹范数，Σ_by＝Σ_n+ADΣ_yD^TA^T，为代价函数。

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

然后，将(8)式带入到(7)式中，得

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

显然，而且可以证明，最小化(7)式再对稀疏信号Y进行 MAP估计，与直接最小化(10)式得到的λ，γ和κ相同，关于Y的解仅相差一个常 量因此，(10)式可以看作是关于稀疏信号Y的正则化回归模型，其中为稀疏信号的非分离稀疏约束。该约束不能拆分成对于Y中每一行的独立约束，因此 该约束能同时约束稀疏信号中非零元素，潜在地考量非零元素之间的相关性。此外，Σ_by中包含了表征噪声强度的λ，因此，得到的稀疏约束可以随着估计的噪声强度自适应 的变化，具有噪声鲁棒性。

4、模型求解。

已知压缩后数据G和采样矩阵A，本发明采用坐标下降法最小化(10)式，每次 迭代中仅优化一个变量而固定剩余的其他变量。具体步骤如下：

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

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

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

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

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

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

其中，为γ^t+1的第i个元素，代表V^T+Y^t+1(Y^t+1)^T 的对角线元素组成的向量，z_i为z的第i个元素；

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

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

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

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

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

上式子中的加法和除法运算均作用在向量的每一个元素上，得到一个新向量，d＝10_-⁶的引入是为了确保γ中出现0时，(18)式依然有意义。

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

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

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