基于三维全变差稀疏先验的高光谱解混压缩感知方法

技术领域

本发明涉及一种高光谱解混压缩感知方法，特别涉及一种基于三维全变差稀疏 先验的高光谱解混压缩感知方法。

背景技术

在高光谱图像中，丰富的光谱信息中包含着极大的数据冗余，严重地增加了高光 谱图像采集、传输以及处理过程中的资源消耗。因此，设计一种高性能的高光谱图像 压缩算法是非常必要的。已有的高光谱图像压缩算法主要分为两类，一类是基于信息 编码的压缩方法，主要利用普通的图像压缩方法去除高光谱图像各个波段内部以及波 段之间的冗余性来实现压缩，其中包括聚类的脉冲差分编码，三维小波变换，三维离 散余弦变换等，然而这类压缩方法作用在图像获取之后，仍然需要消耗大量的资源来 采集和存储数据，而且压缩率较低；另外一类是基于压缩感知(Compressive Sensing) 的压缩方法，该类方法作用于数据采集端，通过采集原始稀疏信号的少量样本点实现 数据的大幅度压缩，解压时利用压缩感知理论实现原始稀疏信号的精确重建，大大减 少了采集和传输过程中的资源消耗。

文献“A compressive sensing and unmixing scheme for hyperspectral data processing, IEEE Transactions on Image Processing,2012,21(3):1200–1210”公开了一种联合光谱解 混的高光谱图像压缩感知算法。该方法首先使用随机观测矩阵对原始数据进行随机采 样实现数据压缩；接着，在对应的光谱库中选择合适的端元；之后利用压缩感知理论 对具有梯度稀疏性的丰度值矩阵进行精确重建；最后，结合选择的端元，利用线性混 合模型实现原始高光谱图像的重建。然而，该方法仅仅考虑了丰度值矩阵空间维上的 梯度稀疏性，忽略了丰度值矩阵光谱维上的潜在稀疏性，从而影响了该压缩算法的重 建精度。

发明内容

为了克服现有联合光谱解混的高光谱图像压缩感知算法精度低的不足，本发明提 供一种基于三维全变差稀疏先验的高光谱解混压缩感知方法。该方法采用随机观测 矩阵从原始数据中抽取少量的样本作为压缩数据。重建过程，根据解混压缩感知模型， 从光谱库中选择适当的光谱作为的模型中的端元矩阵，进而引入丰度值矩阵的三维全 变差稀疏先验，通过求解受限的线性优化问题，精确地求解丰度值矩阵。最后使用线 性混合模型重建原始数据。在HYDICE卫星拍摄的urban数据上当压缩比为1:20时， 归一化的均方误差(normalized mean squared error,NMSE)小于0.09，当压缩比为1:10 时，归一化均方误差同样小于0.08，相对于已有的压缩感知类算法精度提升10%以上。

本发明解决其技术问题所采用的技术方案是：一种基于三维全变差稀疏先验的 高光谱解混压缩感知方法，其特点是包括以下步骤：

步骤一、对于高光谱图像其中每一个像素的光谱x_i表示 成所有端元的线性组合如下：

x_i＝Wh_i   （1）

其中，n_p表示空间上包含的像素数目，n_b表示波段数量，为对应的丰度值向量。

整个数据X表示成丰度值矩阵和端元矩阵W的乘积：

X＝WH   （2）

在H中，行方向是光谱维，每一行代表不同像素的光谱在同一个端元上的投影；列方 向是空间维，每一列代表一个像素的光谱在不同端元上的投影。

步骤二、采用满足高斯随机分布的归一化随机观测矩阵对原始数据进行 随机采样，得到压缩数据如下：

F＝AX＝AWH   （3）

其中，m表示对长度为n_b的信号压缩后的长度，m＜n_b。

步骤三、对于有限的成像场景，根据场景信息从光谱库中抽取n_e个光谱组成端元 矩阵W。

步骤四、(1)在H的光谱维上应用一维的全变差稀疏先验，结合H空间维上的稀疏 性，得到H的三维全变差稀疏先验，如下：

其中，e_j和ε_j分别表示和空间中的第j个单位向量。TV(x)描述的是的全 变差，D_i(x)表示x梯度中的第i个分量。公式（4）中的第一部分表示H空间维上的二 维全变差稀疏先验，其中对应的D_i(·)为二维梯度；第二部分表示H光谱维上的一维全 变差稀疏先验，其中对应的D_i(·)为一维梯度。

(2)构建丰度值的其他先验。将线性解混模型中常用的丰度值先验引入，分别是 混合光谱在不同端元上丰度值投影非负且全和为1的限制，如下：

$1_{n_e}^{T}H=1_{n_p}^T,H>0---\left(5\right)$

其中，和是全部元素为1，长度分别为n_e和n_p的向量。

(3)构建丰度值矩阵H的重建模型。结合公式（3）、（4）和（5）得到如下的重建 模型：

$\left(\begin{array}{c}\underset{H}{\min }\underset{j=1}{\overset{n_p}{\Sigma }}\underset{i=1}{\overset{n_e}{\Sigma }}\left|D_i\left({\mathrm{He}}_j\right)\right|+\underset{j=1}{\overset{n_e}{\Sigma }}\underset{i=1}{\overset{n_p}{\Sigma }}\left|D_i\left({ϵ}_j^{T}H\right)\right|\\ s.t.\mathrm{AWH}=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}\right)---\left(6\right)$

为了方便后续求解，向公式（6）中引入分离变量v_ij＝D_i(He_j),得到：

$\left(\begin{array}{c}\underset{H,{\upsilon }_{\mathrm{ij}},u_{\mathrm{ij}}}{\min }\underset{j=1}{\overset{n_p}{\Sigma }}\underset{i=1}{\overset{n_e}{\Sigma }}\left|v_{\mathrm{ij}}\right|+\underset{j=1}{\overset{n_e}{\Sigma }}\underset{i=1}{\overset{n_p}{\Sigma }}\left|u_{\mathrm{ij}}\right|\\ s.t.v_{\mathrm{ij}}=D_i\left({\mathrm{He}}_j\right),\forall i,j;u_{\mathrm{ij}}=D_i\left({ϵ}_j^{T}H\right),\forall i,j;\mathrm{AWH}=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}\right)---\left(7\right)$

(4)求解公式（7）得到丰度值矩阵H的估计具体求解过程如下：

①使用增广拉格朗日方法，根据公式（7）构建H的增广拉格朗日方程

其中，α＝2⁵,κ＝2⁵,β＝2¹³,γ＝2⁵为二次项惩罚系数，λ_ij,π_ij,Π,υ为对应的拉格朗日乘子， 初始化每个乘子的所有元素为0，||·||_F表示Frobenius范数。

②固定拉格朗日乘子和H，更新分离变量v_ij,u_ij。形式如下：

$\left(\begin{array}{c}{v}_{\mathrm{ij}}=\max \left\{\right|D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{\mathrm{ij}}}{\alpha }|-\frac{1}{\alpha },0\}\mathrm{sgn}(D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{\mathrm{ij}}}{\alpha })\\ u_{\mathrm{ij}}=\max \left\{\right|D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{\mathrm{ij}}}{\kappa }|-\frac{1}{\kappa },0\}\mathrm{sgn}(D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{\mathrm{ij}}}{\kappa })\end{array}\right)---\left(9\right)$

③固定拉格朗日乘子和分离变量v_ij,u_ij，采用梯度下降法更新H。假定第k次更新，由 H^k得到H^k+1，形式如下：

其中，为关于H一阶导数，形式如下：

式中，τ为梯度下降步长。其计算分为初始化和细化两步。在初始化过程中，当第一 次更新H⁰时，τ采用最速下降法进行初始化，之后更新H^k,k＝1,2,...时，对τ采用两点 步长梯度法进行初始化。两点步长梯度法具体形式如下：

其中，tr(·)表示矩阵的迹。τ的细化过程具体如下：

(a)代入初始化的τ，根据公式（10）得到H^k+1，设置参数δ＝3.2×10^-4，η＝0.6和 计数器c＝0；

(b)判断H^k+1是否满足如下的条件：

如果不满足，更新计数器c＝c+1;

如果c＜5，缩小步长τ＝τ·η，继续循环判断是否满足（13）；

否则τ由最速下降法确定，然后由公式（13）得到更新的H^k+1；

否则，得到更新的H^k+1。

④固定更新后的v_ij,u_ij和H，使用如下公式更新拉格朗日乘子：

$\left(\begin{array}{cc}{\lambda }_{\mathrm{ij}}^{k+1}={\lambda }_{\mathrm{ij}}^k-\alpha [D_i\left({\mathrm{He}}_j\right)-v_{\mathrm{ij}}],& {\pi }_{\mathrm{ij}}^{k+1}={\pi }_{\mathrm{ij}}^k-\kappa [D_i\left({ϵ}_j^{T}H\right)-u_{\mathrm{ij}}]\\ {\Pi }^{k+1}={\Pi }^k-\beta (\mathrm{AWH}-F),& {\upsilon }^{k+1}={\upsilon }^k-\gamma {(1_{n_e}^{T}H-1_{n_p}^T)}^T\end{array}\right)---\left(14\right)$

⑤循环步骤②、③和④直至收敛，得到的最终估计的丰度值矩阵

步骤五、结合选择的端元矩阵W和线性混合模型公式（2）得到重建的高光谱数 据

$\hat{X}=W\hat{H}---\left(15\right).$

本发明的有益效果是：该方法采用随机观测矩阵从原始数据中抽取少量的样本作 为压缩数据。重建过程，根据解混压缩感知模型，从光谱库中选择适当的光谱作为的 模型中的端元矩阵，进而引入丰度值矩阵的三维全变差稀疏先验，通过求解受限的线 性优化问题，精确地求解丰度值矩阵。最后使用线性混合模型重建原始数据。在 HYDICE卫星拍摄的urban数据上当压缩比为1:20时，归一化的均方误差(normalized  mean squared error,NMSE)小于0.09，当压缩比为1:10时，归一化均方误差同样小于 0.08，相对于已有的压缩感知类算法精度提升10%以上。

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

具体实施方式

本发明基于三维全变差稀疏先验的高光谱解混压缩感知方法具体步骤如下：

在高光谱图像中，同一像素的不同波段的反射值构成了一个离散向量，称为该像 素的光谱。通常，纯净物质在高光谱图像中具有独一无二的光谱，称作端元。由于存 在地物混合和高光谱图像空间分辨率低等因素，像素的光谱往往是多种不同纯净地物 光谱的混合。这种光谱混合现象通常可以使用线性混合模型来描述。该模型认为任何 一个混合光谱是成像场景中所有端元的线性组合。端元在混合光谱中所占的比例称作 丰度值。因此，对于高光谱图像（n_p表示空间上包含的像素数 目，n_b表示波段数量），其中每一个像素的光谱x_i可以表示成所有端元 的线性组合如下：

x_i＝Wh_i   （1）

其中，为对应的丰度值向量。因此整个数据X可以表示成丰度值矩阵 和端元矩阵W的乘积，如下：

X＝WH   （2）

在H中，行方向是光谱维，每一行代表不同像素的光谱在同一个端元上的投影；列方 向是空间维，每一列代表一个像素的光谱在不同端元上的投影。本发明通过压缩感知 技术重建丰度值矩阵，之后结合选择的端元和线性混合模型重建原始高光谱图像。

1、获取压缩数据。

采用满足高斯随机分布的归一化随机观测矩阵对原始数据进行随机采 样，得到压缩数据如下：

F＝AX＝AWH   （3）

其中，m表示对长度为n_b的信号压缩后的长度，m＜n_b。

2、选择端元。

对于有限的成像场景，通常仅仅包含有限个端元。因此，可以根据场景信息从光 谱库（如ASTER光谱库）中抽取n_e个光谱组成端元矩阵W。

3、重建丰度值矩阵。

根据公式（3）求解丰度值矩阵是一个不适定问题，因此需要在公式（3）的基础 上引入丰度值矩阵H的先验约束。本发明针对H设计了三维全变差稀疏先验，引入到 公式（3）中，最后求解得到重建的具体过程如下：

(1)构建H的三维全变差稀疏先验。H的空间维上集中了不同像素在同一端元上 的丰度值分量。因此，图像场景中的空间相似性在H的空间维上得以保持，而这种相 似性可以使用空间上的二维梯度稀疏性来描述，即空间上的二维全变差稀疏先验。此 外，H的光谱维描述了某个像素的光谱在不同端元上的丰度值投影。通常，在图像匀 质区域内部的光谱，在该类物质对应的端元上具有较大的投影，其他端元上的投影近 似为0，具有一定的稀疏性；而在匀质区域边界处的光谱，往往在两种或两种以上的 端元产生相似的丰度值投影，具有梯度稀疏性。因此，为了满足了这两种情形下的稀 疏性要求，本发明在H的光谱维上应用一维的全变差稀疏先验。结合H空间维上的稀 疏性，得到H的三维全变差稀疏先验，如下：

其中，e_j和ε_j分别表示和空间中的第j个单位向量。TV(x)描述的是的全 变差，D_i(x)表示x梯度中的第i个分量。公式（4）中的第一部分表示H空间维上的二 维全变差稀疏先验，其中对应的D_i(·)为二维梯度；第二部分表示H光谱维上的一维全 变差稀疏先验，其中对应的D_i(·)为一维梯度。

(2)构建丰度值的其他先验。由于本发明中引入了线性混合模型，因此需要将线 性解混模型中常用的丰度值先验引入，分别是混合光谱在不同端元上丰度值投影非负 且全和为1的限制，如下：

$1_{n_e}^{T}H=1_{n_p}^T,H>0---\left(5\right)$

其中，和是全部元素为1，长度分别为n_e和n_p的向量。

(3)构建丰度值矩阵H的重建模型。结合公式（3）、（4）和（5）得到如下的重建 模型：

$\left(\begin{array}{c}\underset{H}{\min }\underset{j=1}{\overset{n_p}{\Sigma }}\underset{i=1}{\overset{n_e}{\Sigma }}\left|D_i\left({\mathrm{He}}_j\right)\right|+\underset{j=1}{\overset{n_e}{\Sigma }}\underset{i=1}{\overset{n_p}{\Sigma }}\left|D_i\left({ϵ}_j^{T}H\right)\right|\\ s.t.\mathrm{AWH}=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}\right)---\left(6\right)$

为了方便后续求解，向（6）中引入分离变量v_ij＝D_i(He_j),得到：

$\left(\begin{array}{c}\underset{H,{\upsilon }_{\mathrm{ij}},u_{\mathrm{ij}}}{\min }\underset{j=1}{\overset{n_p}{\Sigma }}\underset{i=1}{\overset{n_e}{\Sigma }}\left|v_{\mathrm{ij}}\right|+\underset{j=1}{\overset{n_e}{\Sigma }}\underset{i=1}{\overset{n_p}{\Sigma }}\left|u_{\mathrm{ij}}\right|\\ s.t.v_{\mathrm{ij}}=D_i\left({\mathrm{He}}_j\right),\forall i,j;u_{\mathrm{ij}}=D_i\left({ϵ}_j^{T}H\right),\forall i,j;\mathrm{AWH}=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}\right)---\left(7\right)$

(4)求解公式（7）得到丰度值矩阵H的估计具体求解过程如下：

①使用增广拉格朗日方法，根据公式（7）构建H的增广拉格朗日方程

其中，α＝2⁵,κ＝2⁵,β＝2¹³,γ＝2⁵为二次项惩罚系数，λ_ij,π_ij,Π,υ为对应的拉格朗日乘子， 初始化每个乘子的所有元素为0，||·||_F表示Frobenius范数。

②固定拉格朗日乘子和H，更新分离变量v_ij,u_ij。形式如下：

$\left(\begin{array}{c}{v}_{\mathrm{ij}}=\max \left\{\right|D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{\mathrm{ij}}}{\alpha }|-\frac{1}{\alpha },0\}\mathrm{sgn}(D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{\mathrm{ij}}}{\alpha })\\ u_{\mathrm{ij}}=\max \left\{\right|D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{\mathrm{ij}}}{\kappa }|-\frac{1}{\kappa },0\}\mathrm{sgn}(D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{\mathrm{ij}}}{\kappa })\end{array}\right)---\left(9\right)$

③固定拉格朗日乘子和分离变量v_ij,u_ij，采用梯度下降法更新H。假定第k次更新，由 H^k得到H^k+1，形式如下：

其中，为关于H一阶导数，形式如下：

式中，τ为梯度下降步长。其计算分为初始化和细化两步。在初始化过程中，当第一 次更新H⁰时，τ采用最速下降法进行初始化，之后更新H^k,k＝1,2,...时，对τ采用两点 步长梯度法进行初始化。两点步长梯度法具体形式如下：

其中，tr(·)表示矩阵的迹。τ的细化过程具体如下：

(c)代入初始化的τ，根据公式（10）得到H^k+1，设置参数δ＝3.2×10^-4，η＝0.6和 计数器c＝0；

(d)判断H^k+1是否满足如下的条件：

如果不满足，更新计数器c＝c+1;

如果c＜5，缩小步长τ＝τ·η，继续循环判断是否满足（13）；

否则τ由最速下降法确定，然后由公式（13）得到更新的H^k+1；

否则，得到更新的H^k+1。

④固定更新后的v_ij,u_ij和H，使用如下公式更新拉格朗日乘子：

$\left(\begin{array}{cc}{\lambda }_{\mathrm{ij}}^{k+1}={\lambda }_{\mathrm{ij}}^k-\alpha [D_i\left({\mathrm{He}}_j\right)-v_{\mathrm{ij}}],& {\pi }_{\mathrm{ij}}^{k+1}={\pi }_{\mathrm{ij}}^k-\kappa [D_i\left({ϵ}_j^{T}H\right)-u_{\mathrm{ij}}]\\ {\Pi }^{k+1}={\Pi }^k-\beta (\mathrm{AWH}-F),& {\upsilon }^{k+1}={\upsilon }^k-\gamma {(1_{n_e}^{T}H-1_{n_p}^T)}^T\end{array}\right)---\left(14\right)$

⑤循环步骤②、③和④直至收敛，得到的最终估计的丰度值矩阵

4、重建高光谱数据。

结合选择的端元矩阵W和线性混合模型公式（2）得到重建的高光谱数据

$\hat{X}=W\hat{H}---\left(15\right).$

本发明采用随机观测矩阵从原始数据中抽取少量的样本作为压缩数据。重建过程， 根据解混压缩感知模型，从光谱库中选择适当的光谱作为的模型中的端元矩阵，进而 引入丰度值矩阵的三维全变差稀疏先验，通过求解受限的线性优化问题，精确地求解 丰度值矩阵。最后使用线性混合模型重建原始数据。在HYDICE卫星拍摄的urban数 据上当压缩比为1:20时，归一化的均方误差(normalized mean squared error,NMSE)小 于0.09，当压缩比为1:10时，归一化均方误差同样小于0.08，相对于已有的压缩感知 类算法精度提升10%以上。