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

摘要

本发明公开了一种基于三维全变差稀疏先验的高光谱解混压缩感知方法，用于解决现有联合光谱解混的高光谱图像压缩感知算法精度低的技术问题。技术方案是采用随机观测矩阵从原始数据中抽取少量的样本作为压缩数据。重建过程，根据解混压缩感知模型，从光谱库中选择适当的光谱作为的模型中的端元矩阵，进而引入丰度值矩阵的三维全变差稀疏先验，通过求解受限的线性优化问题，精确地求解丰度值矩阵。最后使用线性混合模型重建原始数据。在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)得到如下的重建模型：$\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.AWH=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}---\left(6\right)$为了方便后续求解，向公式(6)中引入分离变量v_ij＝D_i(He_j),得到：$\begin{array}{c}\underset{H,{\upsilon }_{ij},u_{ij}}{\min }\underset{j=1}{\overset{n_p}{\Sigma }}\underset{i=1}{\overset{n_e}{\Sigma }}\left|v_{ij}\right|+\underset{j=1}{\overset{n_e}{\Sigma }}\underset{i=1}{\overset{n_p}{\Sigma }}\left|u_{ij}\right|\\ s.t.v_{ij}=D_i\left({\mathrm{He}}_j\right),\forall i,j;u_{ij}=D_i\left({ϵ}_j^{T}H\right),\forall i,j;AWH=F,1_{n_e}^{T}H=1_{n_p}^T,H>0\end{array}---\left(7\right)$(4)求解公式(7)得到丰度值矩阵H的估计具体求解过程如下：①使用增广拉格朗日方法，根据公式(7)构建H的增广拉格朗日方程其中，α＝2⁵,κ＝2⁵,β＝2¹³,γ＝2⁵为二次项惩罚系数，λ_ij,π_ij,Π,υ为对应的拉格朗日乘子，初始化每个乘子的所有元素为0，||·||_F表示Frobenius范数；②固定拉格朗日乘子和H，更新分离变量v_ij,u_ij；形式如下：$\begin{array}{c}{v}_{ij}=\max \left\{\right|D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{ij}}{\alpha }|-\frac{1}{\alpha },0\}\mathrm{sgn}\left(D_i\left({\mathrm{He}}_j\right)-\frac{{\lambda }_{ij}}{\alpha }\right)\\ u_{ij}=\max \left\{\right|D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{ij}}{\kappa }|-\frac{1}{\kappa },0\}\mathrm{sgn}\left(D_i\left({ϵ}_j^{T}H\right)-\frac{{\pi }_{ij}}{\kappa }\right)\end{array}---\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，使用如下公式更新拉格朗日乘子：$\begin{array}{cc}{\lambda }_{ij}^{k+1}={\lambda }_{ij}^k-\alpha [D_i\left({\mathrm{He}}_j\right)-v_{ij}],& {\pi }_{ij}^{k+1}={\pi }_{ij}^k-\kappa [D_i\left({ϵ}_j^{T}H\right)-u_{ij}]\\ {\Pi }^{k+1}={\Pi }^k-\beta (AWH-F),& {\upsilon }^{k+1}={\upsilon }^k-\gamma {(1_{n_e}^{T}H-1_{n_p}^T)}^T\end{array}---\left(14\right)$⑤循环步骤②、③和④直至收敛，得到的最终估计的丰度值矩阵步骤五、结合选择的端元矩阵W和线性混合模型公式(2)得到重建的高光谱数据$\hat{X}=W\hat{H}---\left(15\right).$