CN103400341B - 基于压缩感知的空谱域联合恢复高光谱数据的方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的空谱域联合恢复高光谱数据的方法，主要解决现有技术资源利用率低、高光谱数据的空域和谱域分辨率相互制约、提高空域和谱域分辨率难度大且成本高的问题。其实现步骤为：根据光谱成像系统的原理设计一种易于硬件实现的观测矩阵；将高光谱数据分为高密度采样数据和低密度采样数据，使用观测矩阵分别对高密度采样数据和低密度采样数据进行空谱域联合采样；利用高光谱数据的空域相关性和谱域相关性，对高光谱数据进行空谱域联合恢复。本发明能在很低的采样率下恢复高光谱数据，大幅度降低采样率，提高高光谱数据的空域和谱域分辨率，可用于高光谱遥感成像。

Description

基于压缩感知的空谱域联合恢复高光谱数据的方法

技术领域

本发明属于图像处理技术领域，涉及一种基于压缩感知的空谱域联合恢复高光谱数据方法，可用于高光谱图像的低速率采样和恢复。

背景技术

高光谱遥感技术已经成为遥感观测的一种重要手段，在气象、环境检测、地质勘探、农业生产、航空航天以及军事探测等领域广泛应用。高光谱成像系统可以获取成像地区的地物的几何、辐射和光谱信息，是高光谱遥感的核心技术。

遥感高光谱成像系统一般都用于卫星或机载，使得无论增加谱域采样提高光谱分辨率，还是增加空域采样提高空间分辨率成本都非常高。高光谱数据是一种三维立方体数据，本身数据量很大，提高谱域和空域的采样率会使得数据量大幅度增加，这给数据在卫星通道上的高效传输和存储带来了极大的挑战。

高光谱成像系统主要由望远成像部分和光谱成像部分构成，望远成像部分对地面一定区域成像后通过一个狭缝到达光谱成像系统，光谱成像系统通过分光仪分光，形成不同谱段的像。高光谱成像系统的空域和谱域分辨率是互相制约的，如果想提高空域分辨率就要减小狭缝，这样通光量减小，在总的光能不变的情况下，为了使每一个谱段的光能达到成像的要求，必须减少谱段数，这样就会降低谱域分辨率。如果想提高谱域分辨率，增加成像的谱段数，就必须增加通光量，即增大狭缝宽度，这样就会降低空域分辨率。

传统的采样系统都是按照奈奎斯特采样定理进行采样，将得到的数据进行编码压缩，去除大量的冗余信息，再进行存储和传输。这种高速采样再压缩的方法极大的浪费了采样资源。

压缩感知（CompressiveSensing,CS）是一种全新的理论，它能突破传统的奈奎斯特采样定理的限制，实现信号的高效、高精度采样及重构。压缩感知理论的基本思想是信号在某个变换域满足稀疏性或可压缩性的情况下，利用一个与变换基不相关的矩阵对高维信号进行采样，得到采样数据后通过求解一个优化问题就可以重构原始信号。该方法的数学过程简要描述如下：对于N维信号x∈R^N×1，通过一个M×N的观测矩阵Φ_M×N（M＜＜N）进行采样得到观测数据y＝Φ_M×Nx,y∈R^M×1，最后通过优化方法求解得到x。

核压缩感知（KCS，KernelCompressiveSensing）通过引入核方法，在核空间上进行观测和恢复，可以以更低的采样率实现更好的重构效果，而且不用求解优化问题。

利用压缩感知技术实现高光谱数据的压缩感知采样，可以有效解决高光谱数据量大、传输难度大、成像系统成本高等难题。现有的高光谱压缩感知采样方法，都是单独在空域或者谱域来进行的，没有联合考虑高光谱图像的空域相关性和谱域相关性；另一方面现有的高光谱压缩感知采样方法中，采样矩阵都采用高斯随机矩阵，高斯随机矩阵用硬件难以实现，而且成本高。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于压缩感知的空谱域联合恢复高光谱数据的方法，以大幅度降低采样率，降低高光谱成像系统的实现难度，减少成本，提高图像的空域和谱域分辨率。

实现本发明目的技术思路是：通过联合考虑空域和谱域的相关性，在谱域和空域联合采样并进行恢复，并根据高光谱成像仪的成像原理，设计一种新的采样矩阵，易于硬件实现。具体步骤如下：

1一种基于压缩感知的空谱域联合恢复高光谱数据的方法，包括：

（1）观测矩阵设计步骤：

1a)根据待采样高光谱图像的谱域和空域的维数N和采样数M，得到压缩采样率r＝M/N；

1b)建立一个M行N列的矩阵Φ_M×N，该矩阵的第i行第j列Φ_M×N(i,j)取值为：

（2）空谱域联合压缩感知采样步骤：

2a)给定待采样的高光谱立方体数据X_m×k×l，并将X_m×k×l在空域向量化后表示为X_m×n，m为光谱域维数，k,l为高光谱立方体数据空域的维数，n=k×l为空域向量化后空域的维数；

2b)把向量化后的高光谱立方体数据X_m×n按照谱段分为高密度采样数据X^DG和低密度采样数据X^SG，在X_m×n中取出p行记作X_p×n，剩余部分记作X_q×n，把X_p×n作为高密度采样数据X^DG=X_p×n，把X_q×n作为低密度采样数据X^SG=X_q×n，p为低高度采样数据的维数，q为低密度采样数据的维数，0＜p＜m，0＜q＜m，p+q＝m；

2c)设定高密度采样数据X^DG的空域压缩采样率和谱域压缩采样率设定低密度采样数据X^SG的空域压缩采样率和谱域压缩采样率根据设定的这些数据按照步骤(1),分别构造高密度采样数据的空域观测矩阵谱域观测矩阵低密度采样部分的空域观测矩阵和谱域观测矩阵

2d)分别对高密度采样数据X^DG和低密度采样数据X^SG进行空谱域联合采样，得到高密度采样数据的采样值y^DG和低密度采样数据的采样值y^SG：

$y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{DG}}{\left(X^{\mathrm{DG}}\right)}^T\right)}^T,$ $y^{\mathrm{SG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{SG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{SG}}{\left(X^{\mathrm{SG}}\right)}^T\right)}^T;$

其中，T表示矩阵转置；

（3）对高密度采样数据采样值进行恢复：

3a)定义第一中间变量表示则将高分辨率采样数据X^DG的空谱域联合采样值 $y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{DG}}{\left(X^{\mathrm{DG}}\right)}^T\right)}^T$ 转化为： $y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{y}_{\mathrm{spat}}^{\mathrm{DG}};$

3b)对式进行求解，得到第一中间变量的恢复值

3c)用恢复值替换公式中的项，将该公式重新表示为用BP（BasisPursuit,基追踪）算法求解得到高密度采样数据X^DG的恢复值

（4）对低密度采样数据采样值进行恢复：

4a)定义第二中间变量表示则将低密度采样数据X^SG的空谱域联合采样值转化为用TV（TotalVariationMinimization，最小全变分法）方法求解式得到第二中间变量的恢复值

4b)用第二中间变量恢复值替换公式中的项，将该公式重新表示为用步骤3c）得到的高密度采样数据X^DG的恢复值作为训练数据，用KCS（KernelcompressiveSensing,核压缩感知）算法求解得到低密度采样数据X^SG的恢复值

4c)将高密度采样数据X^DG的恢复值和低密度采样数据X^SG的恢复值合并，得到待恢复数据X_m×n的恢复值

本发明与现有技术相比具有以下优点：

1、本发明根据高光谱成像的原理设计了观测矩阵，所设计的观测矩阵具有固定结构，易于硬件实现，现有的高光谱压缩感知采样方法所采用的观测矩阵都是随机矩阵，随机矩阵难于用硬件实现，并且实现成本很高，本发明很好的克服了这个问题。

2、本发明采用空谱域的联合压缩感知恢复技术，突破了现有方法的空域和谱域分辨率相互制约的限制，可以同时提高高光谱数据的空域和谱域分辨率，极大的降低了采样系统的采样率，提高了采样资源的利用效率。

附图说明

图1是本发明的流程图；

图2是本发明利用的光谱成像原理图；

图3是本发明设计采样矩阵的示意图；

图4是本发明在采样率为4%时的不同波段的恢复图像和原图像；

图5是在采样率为4%时，用本发明对高光谱数据进行恢复和在采样率50%时，用现有的正交匹配追踪算法和基追踪算法分别在空域和谱域对高光谱数据进行恢复的结果对比图。

具体实施方式

参照图1，本发明的具体实施步骤如下：

步骤1：根据光谱成像原理设计采样矩阵：

1.1）给定光谱成像系统如图2所示，其中L1为望远物镜，用来接收地面反射光；S为狭缝，用来控制地面一定区域的反射光进入成像系统；L2为准直镜，用来将通过狭缝的入射光变为平行光；Prism为分光器件，用来把入射光分为不同的谱段；L3为成像物镜，用来对不同谱段的入射光在CCD上成像；

1.2）给定光谱成像系统的工作波长范围为λ₁～λ₂，成像谱段数为K，则光谱成像系统的光谱分辨率为Δλ＝(λ₁-λ₂)/K，光谱成像系统在第i个谱段得到的采样值为其中λ₁为光谱成像系统工作波长的下限，λ₂为光谱成像系统工作波长的上限，f(λ)为光谱成像仪的光谱响应函数，λ为表示进入光谱成像系统的入射光波长的变量，c为第i个谱段的起始波长，i＝{1,2,3,...,K}；

1.3）假设高分辨率光谱成像系统G和低分辨率光谱成像系统D的成像谱段数分别为N和M，工作波长范围都为λ₃～λ₄，光谱响应函数都为f₁(λ)，则高分辨率成像系统G的光谱分辨率为Δλ_h＝(λ₁-λ₂)/N，低分辨率成像系统D的光谱分辨率为Δλ_l＝(λ₁-λ₂)/M，其中λ₃为光谱成像系统工作波长的下限，λ₄为光谱成像系统工作波长的上限；

1.4）高分辨率光谱成像系统G对一固定场景成像的结果为x＝{x₁,x₂,...x_j...,x_N},其中x_j为高分辨率光谱成像系统G的第j个谱段的采样值，j＝{1,2,3,...,N}，根据步骤1.2可得d为高分辨率光谱成像系统G的第j个波段的起始波长；

1.5）低分辨率光谱成像系统D在相同的曝光时间段内对同一场景的成像结果为y＝{y₁,y₂,...y_k...,y_M}，其中y_k为低分辨率光谱成像系统D的第k个谱段的采样值k＝{1,2,3,...,M}，根据步骤1.2可得e为低分辨率光谱成像系统D的第k个波段的起始波长；

1.6）定义压缩采样率r＝N/M，高分辨率光谱成像系统G对一固定场景成像的结果x＝{x₁,x₂,...x_j...,x_N}和低分辨率光谱成像系统D在相同的曝光时间段内对同一场景的成像结果y＝{y₁,y₂,...y_k...,y_M}的关系如附图3所示（这里假设r=3），根据光能守恒，在相同的波长范围内，每一个低分辨率成像仪的采样值等价于3个高分辨率成像仪的采样值之和，这样就可以得到：

记上式为：y＝Φ_M×N*x,其中Φ_M×N是M行N列的观测矩阵，Φ_M×N(i,j)表示矩阵第i行第j列的取值，其取值如下式定义：

步骤2：空谱域联合压缩感知采样：

2.1）给定待采样的高光谱立方体数据X_m×k×l，并将X_m×k×l在空域向量化后表示为X_m×n，m为光谱域维数，k,l为高光谱立方体数据空域的维数，n=k×l为空域向量化后空域的维数；

2.2）把向量化后的高光谱立方体数据X_m×n按照谱段分为高密度采样数据X^DG和低密度采样数据X^SG，在X_m×n中取出p行记作X_p×n，剩余部分记作X_q×n，把X_p×n作为高密度采样数据X^DG=X_p×n，把X_q×n作为低密度采样数据X^SG=X_q×n，p为低高度采样数据的维数，q为低密度采样数据的维数，0＜p＜m，0＜q＜m，p+q＝m；

2.3）设定高密度采样数据X^DG的空域压缩采样率和谱域压缩采样率设定低密度采样数据X^SG的空域压缩采样率和谱域压缩采样率根据设定的这些数据按照步骤1,分别构造高密度采样数据的空域观测矩阵谱域观测矩阵低密度采样部分的空域观测矩阵和谱域观测矩阵

2.4）分别对高密度采样数据X^DG和低密度采样数据X^SG进行空谱域联合采样，得到高密度采样数据的采样值y^DG和低密度采样数据的采样值y^SG：

$y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{DG}}{\left(X^{\mathrm{DG}}\right)}^T\right)}^T,y^{\mathrm{SG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{SG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{SG}}{\left(X^{\mathrm{SG}}\right)}^T\right)}^T;$

其中，T表示矩阵转置。

步骤3：对高密度采样数据采样值进行恢复：

3.1）定义第一中间变量表示则将高分辨率采样数据X_DG的空谱域联合采样值 $y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{\left({\mathrm{\Φ}}_{\mathrm{spat}}^{\mathrm{DG}}{\left(X^{\mathrm{DG}}\right)}^T\right)}^T$ 转化为 $y^{\mathrm{DG}}={\mathrm{\Φ}}_{\mathrm{spec}}^{\mathrm{DG}}{y}_{\mathrm{spat}}^{\mathrm{DG}};$

3.2）对式进行求解，得到第一中间变量的恢复值

式是一个典型的压缩感知问题，通过将其转化为一个零范数优化问题来求解，现有的求解该问题的方法主要有：梯度投影算法（GPSR）、匹配追踪算法MP、正交匹配追踪算法（OMP）、基追踪（BP）、最小全变分法（TV）、贝叶斯压缩感知算法（BCS），本实例中采用基追踪（BP）来求解该问题；

3.3）用恢复值替换公式中的项，将该公式重新表示为用基追踪BP算法求解得到高密度采样数据X^DG的恢复值

步骤4：对低密度采样数据采样值进行恢复：

4.1）定义第二中间变量表示则将低密度采样数据X^SG的空谱域联合采样值转化为用最小全变分法（TV）方法求解式得到第二中间变量的恢复值

4.2）用第二中间变量恢复值替换公式中的项，将该公式重新表示为并用步骤3得到的高密度采样数据X^DG的恢复值作为训练数据，采用核压缩感知KCS算法求解

4.2a）构造核函数K(a,b)＝g(<a,b>)＝(<a,b>+cons)^index，a,b为两个不同的函数变量，cons为常数，cons≥0,index为核函数的指数，index取大于0的奇数；

4.2b）用步骤3得到的高密度采样数据X^DG的恢复值作为训练数据，用4.2a）构造的核函数做核主成分分析（KPCA），得到特征值向量Ψ和特征值向量α；

4.2c）用特征值向量Ψ作为稀疏基，对低密度采样数据的转置(X^SG)^T进行核空间的稀疏表示，得到低密度采样数据的转置(X^SG)^T在核空间的稀疏表示式：f((X^SG)^T)=Ψβ，函数f(.)为核空间的投影函数，β为稀疏表示系数；

4.2d）对式 ${\hat{y}}_{\mathrm{spat}}^{\mathrm{SG}}={\left({\mathrm{\&Phi;}}_{\mathrm{spat}}^{\mathrm{SG}}{\left(X^{\mathrm{SG}}\right)}^T\right)}^T$ 两边求转置得到 ${\left({\hat{y}}_{\mathrm{spat}}^{\mathrm{SG}}\right)}^T={\mathrm{\&Phi;}}_{\mathrm{spat}}^{\mathrm{SG}}{\left(X^{\mathrm{SG}}\right)}^T;$

4.2e）将投影到核空间得到：

$f\left({\left({\hat{y}}_{\mathrm{spat}}^{\mathrm{SG}}\right)}^T\right)=K({\mathrm{\&Phi;}}_{\mathrm{spat}}^{\mathrm{SG}},{\left(X^{\mathrm{SG}}\right)}^T)=K({\mathrm{\&Phi;}}_{\mathrm{spat}}^{\mathrm{SG}},\mathrm{\&Psi;})\mathrm{\&beta;};$

4.2f）通过最小二乘法求解式子得到β的估计值

4.2g）把带入f((X^SG)^T)=Ψβ得到得到(X^SG)^T在核空间的稀疏表示的值

4.2h）对求逆得到低密度采样数据X^SG的恢复值

本实例所采用的核压缩感知算法KCS详见H.QiandS.Hughes,“Usingthekerneltrickincompressivesensing:Accuratesignalrecoveryfromfewermeasurements”,inIEEEInt.ConferenceonAcoustics,SpeechandSignalProc.may2011,pp.3940-3943；

4.3）将高密度采样数据X^DG的恢复值和低密度采样数据X^SG的恢复值合并，得到待恢复数据X_m×n的恢复值

本发明的效果通过以下仿真实验进一步说明：

（1）实验仿真条件：

本实验采用数据是典型的AVIRIS高光谱图像：该数据是美国印第安纳州西北部印第安遥感试验区的高光谱图像，总共有16类地物，成像时间为1992年6月。数据共有200个波段，每个波段图像的大小为145×145，每个像素16bits。本实验在CPU为Intel(R)Pentium(R)4、主频2.2GHz，内存为1G的WINDOWS7系统上采用软件MATLAB2011a进行仿真。

（2）恢复数据评价标准：

（2a）空域评价标准：峰值信噪比PSNR

高光谱图像在计算峰值信噪比时先将数据归一化。

峰值信噪比的定义式如下：

$\mathrm{PSNR}=101g\left\{\frac{1}{\frac{1}{\mathrm{MN}}\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{[f(i,j)-\hat{f}(i,j)]}^2}\right\}$

f(i,j)为原数据，为恢复的数据。PSNR值越大，恢复图像与原始图像的差别就越小，图像质量越高。

（2b）谱域评价标准：相关系数

相关系数的定义如下：

其中u为恢复的数据，v为原始数据。Du，Dv分别为u，v的标准差。相关系数衡量原始数据和恢复数据之间的相关性。相关系数越大，表示相关性越大，相关系数的R的最大值为1，最小值为0。

（3）实验内容：

实验一

在高密度采样数据X^DG的维数p=50,空域压缩采样率谱域压缩采样率低密度采样数据X^SG的维数为q=150，空域压缩采样率谱域压缩采样率总的采样率的条件下，利用本发明设计的观测矩阵对实验数据进行采样和恢复，实验结果如图4所示，其中:

图4（a）是第10个波段的原图，

图4（b）是第10个波段的恢复图，其中Psnr=31.6946，

图4（c）是第50个波段的原始图，

图4（d）是第50个波段的原始图,其中Psnr=31.4021，

从图4可见，本发明可以在很低的采样率下，可以得到很好的空域恢复效果。

分别计算恢复数据中的16类地物和原始数据中的16类地物的谱域相关系数，计算结果如表1所示

表1恢复数据中的16类地物和原始数据中的16类地物的谱域相关系数

类别号	相关系数	类别号	相关系数	平均相关系数
					0	0.9807	9	0.9766	0.9809
1	0.9825	10	0.9828
					2	0.9818	11	0.9805
3	0.9815	12	0.9795
					4	0.9764	13	0.9824
5	0.9828	14	0.9822
					6	0.9820	15	0.9813
7	0.9810	16	0.9798
					8	0.9818

从表1可见，恢复数据和原始数据的谱线相关性很高，本发明在采样率很低的条件下，可以很好的恢复高光谱数据。

实验二

对实验数据在空域进行50%的采样，分别用OMP算法和BP算法进行恢复，对实验数据在谱域进行50%的采样，分别用OMP算法和BP算法进行恢复，与本发明在采样率4%时恢复的结果进行对比，实验结果如图5所示，其中：

图5（a）是第30波段的原始图，

图5（b）是本发明恢复的第30波段的恢复图，其Psnr=30.4146，

图5（c）是OMP方法在谱域恢复的第30波段的恢复图Psnr=27.0826，

图5（d）是BP方法在谱域恢复的第30波段的恢复图Psnr=29.0444，

图5（e）是OMP方法在空域恢复的第30波段的结果图Psnr=23.0682，

图5（f）是BP方法在空域恢复的第30波段的恢复图Psnr=26.6664；

从图5可见，本发明在采样率4%时得到的空域恢复效果优于BP和OMP在采样率50%时得到的空域恢复效果；

将本发明的方法与现有方法恢复结果的数值指标在空域和谱域进行对比，如表2所示。

表2本发明的方法与现有方法恢复结果的数值指标对比

从表2可见，本发明在采样率4%时的谱域和空域恢复效果已经明显的优于BP和OMP在采样率50%时的谱域和空域的恢复效果。

综上，本发明可以大幅度降低高光谱采样设备的采样率，提高高光谱图像的空域和谱域分辨率，提高采样资源利用效率，降低采样系统设计难度和设备制造成本，在高光谱遥感成像领域具有广阔的应用前景。

Claims (1)

1.一种基于压缩感知的空谱域联合恢复高光谱数据的方法，包括：

(1)观测矩阵设计步骤：

1a)根据待采样高光谱图像的谱域和空域的维数N和采样数M，得到压缩采样率r＝M/N；

1b)建立一个M行N列的矩阵Φ_M×N，该矩阵的第i行第j列Φ_M×N(i,j)取值为：

(2)空谱域联合压缩感知采样步骤：

2a)给定待采样的高光谱立方体数据X_m×k×l，并将X_m×k×l在空域向量化后表示为X_m×n，m为光谱域维数，k,l为高光谱立方体数据空域的维数，n＝k×l为空域向量化后空域的维数；

2b)把空域向量化后的高光谱立方体数据X_m×n按照谱段分为高密度采样数据X^DG和低密度采样数据X^SG，在X_m×n中取出p行记作X_p×n，剩余部分记作X_q×n，把X_p×n作为高密度采样数据X^DG＝X_p×n，把X_q×n作为低密度采样数据X^SG＝X_q×n，p为高密度采样数据的维数，q为低密度采样数据的维数，0<p<m，0<q<m，p+q＝m；

2c)设定高密度采样数据X^DG的空域压缩采样率和谱域压缩采样率设定低密度采样数据X^SG的空域压缩采样率和谱域压缩采样率根据设定的这些数据按照步骤(1),分别构造高密度采样数据的空域观测矩阵谱域观测矩阵低密度采样部分的空域观测矩阵和谱域观测矩阵

2d)分别对高密度采样数据X^DG和低密度采样数据X^SG进行空谱域联合采样，得到高密度采样数据的采样值y^DG和低密度采样数据的采样值y^SG：

其中，T表示矩阵转置；

(3)对高密度采样数据采样值进行恢复：

3a)定义第一中间变量表示则将高分辨率采样数据X^DG的空谱域联合采样值转化为：

3b)对式进行求解，得到第一中间变量的恢复值

3c)用恢复值替换公式中的项，将该公式重新表示为用基追踪BP算法求解得到高密度采样数据X^DG的恢复值

(4)对低密度采样数据采样值进行恢复：

4a)定义第二中间变量表示则将低密度采样数据X^SG的空谱域联合采样值转化为用最小全变分法TV求解式得到第二中间变量的恢复值

4b)用第二中间变量恢复值替换公式中的项，将该公式重新表示为用步骤3c)得到的高密度采样数据X^DG的恢复值作为训练数据，用核压缩感知KCS算法求解得到低密度采样数据X^SG的恢复值

4c)将高密度采样数据X^DG的恢复值和低密度采样数据X^SG的恢复值合并，得到待恢复数据X_m×n的恢复值