CN104933685A - 基于三维张量压缩感知的高光谱压缩成像方法 - Google Patents

Abstract

本发明公开了一种基于三维张量压缩感知的高光谱压缩成像方法，主要解决现有技术在多维压缩感知过程中多维数据的结构信息被破坏的问题。其实现步骤为：通过引入张量的方法，在对高光谱图像不进行向量化处理的基础上，先在高光谱图像的三个维度上同时进行压缩采样，获得测量值；然后求出三个维度上的感知矩阵；再利用张量正交匹配追踪算法求出稀疏系数；最后根据高光谱图像的多维稀疏表示，完成高光谱图像的重构。实验结果表明：在相同采样率下，本发明同其他传统的压缩感知方法相比，重构速度快，效果好，可用于遥感图像获取。

Description

基于三维张量压缩感知的高光谱压缩成像方法

技术领域

本发明属于信号处理技术领域，涉及一种高光谱压缩成像方法，可用于遥感图像获取。

背景技术

压缩感知是近年来图像处理技术领域中发展起来的一种新的采样理论，通过利用信号的稀疏特性，可在远小于传统奈奎斯特采样率的条件下，实现信息的精确恢复。压缩感知算法通常都是对一维信号进行操作，当信号维度超过一维时，通常都是对其进行向量化，把信号转换成一维向量，然后再进行压缩感知操作。而多维信号的各个维度间存在结构相关性，若是简单的向量化，将会破坏多维信号各个维度间的结构，从而失去结构信息，给后续的信号处理增加了困难。对于高光谱图像来说，空间和谱间存在极强的相关性，且维度较大，若对其向量化，不但失去了结构信息，而且增加了算法在时间和空间上的复杂度。在对一维信号进行压缩感知恢复时，采用OMP、BP的算法，很容易求解。若是对于多维信号，普通的OMP、BP方法虽然能够恢复出原信号，但是由于多维信号数据量大，因此运算效率低，耗时较长，特别是BP算法，重构的效果并不特别理想。Hanchao Qi和Shannon Hughes等人提出了一种基于核技巧的压缩感知方法，引入核函数，将观测的过程和恢复重构过程映射到核空间上。这种方法在一定程度上提高了重构效果，但是仍旧沿用对图像进行向量化处理的方法，破坏了多维信号间的结构信息，造成了信息丢失，对后续的压缩感知工作带来困难。

发明内容

本发明的目的在于针对上述已有技术的不足，提出了一种基于三维张量压缩感知的高光谱压缩成像方法，以提高运算效率和图像恢复效果。

本发明的技术关键是，通过引入张量，在高光谱图像的各个维度上同时进行压缩采样成像，利用张量正交匹配追踪算法求出稀疏系数，进而完成图像重构。具体步骤包括如下：

(1)输入高光谱图像，将其表示为三维张量Z，其中I₁、I₂、I₃分别为高光谱图像三个维度的大小；

(2)设三个维度的采样率分别为S₁、S₂、S₃，构造一个具有克罗内克结构的观测矩阵Φ：

$\mathrm{\Φ}={\mathrm{\Φ}}_3\⊗{\mathrm{\Φ}}_2\⊗{\mathrm{\Φ}}_1,$

其中，Φ₁、Φ₂、Φ₃是大小分别为J₁×I₁、J₂×I₂和J₃×I₃的高斯随机矩阵，作为三个维度的观测矩阵，第i个维度的观测矩阵Φ_i的行数J_i由第i个维度的采样率S_i和高光谱图像第i个维度的大小I_i来决定，即J_i＝S_iI_i，i＝1,2,3， 表示矩阵的克罗内克积；

(3)根据高光谱图像Z和三个维度的观测矩阵Φ₁、Φ₂、Φ₃，得到低维的测量值M：

M＝Z×₁Φ₁×₂Φ₂×₃Φ₃，

其中，×_i表示高光谱数据Z与第i个维度上的观测矩阵Φ_i的张量i-模乘；

(4)构造一个具有克罗内克结构的字典D：

$D=D_3\⊗D_2\⊗D_1$

其中，D₁、D₂、D₃是大小分别为I₁×I₁、I₂×I₂和I₃×I₃的多维度上的字典，三个维度的字典均取为离散余弦字典；

(5)根据三个维度的观测矩阵Φ₁、Φ₂、Φ₃和三个维度的字典D₁、D₂、D₃，计算三个维度上的压缩感知矩阵：

Q₁＝Φ₁D₁，

Q₂＝Φ₂D₂，

Q₃＝Φ₃D₃，

得到压缩感知矩阵 $Q=Q_3\⊗Q_2\⊗Q_1,$ 其中

(6)根据测量值M和各个维度上的压缩感知矩阵Q₁、Q₂、Q₃，利用张量正交匹配追踪算法求解下式，得出稀疏系数张量β：

M＝β×₁Q₁×₂Q₂×₃Q₃，

其中，×_i表示稀疏系数张量β和第i个维度上的压缩感知矩阵Q_i的张量i-模乘；

(7)根据稀疏系数β和三个维度上的字典D₁、D₂、D₃，得到重构出的原高光谱图像

$\underset{\‾}{\hat{Z}}=\underset{\‾}{\mathrm{\β}}{\×}_1{D}_1{\×}_2{D}_2{\×}_3{D}_3,$

其中，×_i表示稀疏系数张量β和第i个维度上的字典D_i的张量i-模乘，

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

本发明由于引入张量，使得在对高光谱图像进行压缩感知的过程中，多维数据各维度间的结构信息保留下来，这样不仅能够获得更好的恢复效果，而且也大大减小了压缩成像的复杂度，缩短了实验的运行时间。

附图说明

图1是本发明的实现流程图；

图2是用本发明方法在采样率为6.25％时对高光谱图像Moffet的重构效果图；

图3是用本发明方法在采样率为12.5％时对高光谱图像Moffet的重构效果图。

具体实施方法

参照图1，本发明的实现步骤如下：

步骤1.输入高光谱图像，将其表示为三维张量Z，其中，I₁、I₂、I₃分别为高光谱图像三个维度的大小。

步骤2.构造观测矩阵。

设三个维度的采样率分别为S₁、S₂、S₃，构造一个具有克罗内克结构的观测矩阵Φ：

$\mathrm{\Φ}={\mathrm{\Φ}}_3\⊗{\mathrm{\Φ}}_2\⊗{\mathrm{\Φ}}_1,$

其中，Φ₁、Φ₂、Φ₃是大小分别为J₁×I₁、J₂×I₂和J₃×I₃的高斯随机矩阵，分别作为三个维度的观测矩阵，且第i个维度的观测矩阵Φ_i的行数J_i由第i个维度的采样率S_i和高光谱图像第i个维度的大小I_i来决定，即J_i＝S_iI_i，i＝1,2,3， 表示矩阵的克罗内克积。

步骤3.求出测量值。

根据高光谱图像Z和三个维度的观测矩阵Φ₁、Φ₂、Φ₃，得到低维的测量值M：

M＝Z×₁Φ₁×₂Φ₂×₃Φ₃，

其中，×_i表示高光谱图像χ与第i个维度上的观测矩阵Φ_i的张量i-模乘，按如下方式计算：

在第1个维度上按公式计算高光谱图像Z与第1维度的观测矩阵Φ₁的1-模乘积，得到Z×₁Φ₁，是一个张量，计算的过程相当于，高光谱图像Z其他维度不变，第1维度的元素构成一个向量，则第1维度的观测矩阵Φ₁与该向量相乘，其乘积作为测量值M第1维度的元素，其中，表示张量(Z×₁Φ₁)的j₁i₂i₃位置上的元素，表示高光谱图像Z的i₁i₂i₃位置上的元素，表示第1个维度上的观测矩阵Φ₁的j₁i₁位置上的元素，j₁＝1,2,…,J₁，i₁＝1,2,…,I₁，i₂＝1,2,…,I₂，i₃＝1,2,…,I₃；

在第2个维度上按公式计算高光谱图像Z与第2维度的观测矩阵Φ₂的2-模乘积，得到Z×₂Φ₂，是一个张量，计算的过程相当于，高光谱图像Z其他维度不变，第2维度的元素构成一个向量，则第2维度的观测矩阵Φ₂与该向量相乘，其乘积作为测量值M第2维度的元素，其中，表示张量(Z×₂Φ₂)的i₁j₂i₃位置上的元素，表示高光谱图像Z的i₁i₂i₃位置上的元素，表示第2个维度上的观测矩阵Φ₂的j₂i₂位置上的元素，j₂＝1,2,…,J₂，i₁＝1,2,…,I₁，i₂＝1,2,…,I₂，i₃＝1,2,…,I₃；

在第3个维度上按公式计算高光谱图像Z与第3维度的观测矩阵Φ₃的3-模乘积，得到Z×₃Φ₃，是一个张量，计算的过程相当于，高光谱图像Z其他维度不变，第3维度的元素构成一个向量，则第3维度的观测矩阵Φ₃与该向量相乘，其乘积作为测量值M第3维度的元素，其中，表示张量(Z×₃Φ₃)的i₁i₂j₃位置上的元素，表示高光谱图像Z的i₁i₂i₃位置上的元素，表示第3个维度上的观测矩阵Φ₃的j₃i₃位置上的元素，j₃＝1,2,…,J₃，i₁＝1,2,…,I₁，i₂＝1,2,…,I₂，i₃＝1,2,…,I₃；

在三个维度上计算完成后，得到测量值M，

步骤4.构造一个具有克罗内克结构的字典D：

$D=D_3\⊗D_2\⊗D_1$

其中，D₁、D₂、D₃是大小分别为I₁×I₁、I₂×I₂和I₃×I₃的三个维度上的字典，作三个维度的字典均取为离散余弦字典。

步骤5.计算感知矩阵。

根据三个维度的观测矩阵Φ₁、Φ₂、Φ₃和三个维度的字典D₁、D₂、D₃，计算三个维度上的压缩感知矩阵：

Q₁＝Φ₁D₁，

Q₂＝Φ₂D₂，

Q₃＝Φ₃D₃，

得到压缩感知矩阵 $Q=Q_3\⊗Q_2\⊗Q_1,$ 其中，

步骤6.求出稀疏系数。

根据测量值M和各个维度上的压缩感知矩阵Q₁、Q₂、Q₃，利用张量正交匹配追踪算法求解下式，得出稀疏系数张量β：

M＝β×₁Q₁×₂Q₂×₃Q₃，

其中，×_i表示稀疏系数张量β和第i个维度上的压缩感知矩阵Q_i的张量i-模乘；

6a)设初始迭代次数k＝1，残差R ₁＝M，索引集稀疏度为T＝1000，门限值ε＝0.05；

6b)计算第k次循环的残差R _k与三个维度上的压缩感知矩阵Q₁、Q₂、Q₃中未选中的列向量的乘积H：

$\underset{\‾}{H}={\underset{\‾}{R}}_k{\×}_1{Q}_1^T(:,i_1){\×}_2{Q}_2^T(:,i_2){\×}_3{Q}_3^T(:,i_3),$

其中，i₁、i₂、i₃分别表示每个感知矩阵中未选中的列向量的索引，()^T表示矩阵的转置，i₁＝1,2,…,I₁，i₂＝1,2,…,I₂，i₃＝1,2,…,I₃，

6c)对乘积H各个维度上的值取绝对值，找出各个维度上绝对值最大的H所对应的感知矩阵Q₁、Q₂、Q₃中的列向量的索引i₁ ^k、i₂ ^k、i₃ ^k，然后把索引i₁ ^k、i₂ ^k、i₃ ^k添加到索引集Λ_i中，即

6d)用当前索引集Λ_i所对应的压缩感知矩阵Q_i的列向量构成子矩阵：P_i＝Q_i(:,Λ_i)；

6e)利用递归的算法计算测量值M在当前的子矩阵P_i下的稀疏系数β：

$\underset{\‾}{\mathrm{\β}}=\underset{\underset{\‾}{\mathrm{\β}}}{\arg \min }{\left|\right|\underset{\‾}{M}-\underset{\‾}{\mathrm{\β}}{\×}_1{P}_1{\×}_2{P}_2{\×}_3{P}_3\left|\right|}_2^2,$

其中，稀疏系数β满足||β||₀≤T，||||₀表示求张量的0范数；

6f)迭代次数加1，即k＝k+1；

6g)利用下式更新残差R _k：

R _k＝M-β×₁P₁×₂P₂×₃P₃；

6h)判断是否满足|Λ₁||Λ₂||Λ₃|>T，或者||R _k||_F<ε，其中，|Λ_i|表示索引集Λ_i包含的元素数，||||_F表示Frobenius范数，如果不满足，则重复步骤6b)到步骤6g)，如果满足，停止迭代，得出稀疏系数β。

步骤7.重构原高光谱图像。

根据稀疏系数β和三个维度上的字典D₁、D₂、D₃，得到重构出的原高光谱数据

$\underset{\&OverBar;}{\hat{Z}}=\underset{\&OverBar;}{\mathrm{\&beta;}}{\&times;}_1{D}_1{\&times;}_2{D}_2{\&times;}_3{D}_3,$

其中，×_i表示稀疏系数张量β和第i个维度上的字典D_i的张量i-模乘，

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

1)实验条件

本实验所用的高光谱图像为典型的VIRIS高光谱图像Moffett，Moffett图像是由1992年8月AVIRIS传感器对加利福尼亚州的Moffett地区成像所得，具有0.4um～2.5um波长范围内的224个光谱带,去掉所有像素为0和不透明的波段后有200个谱段，空间分辨率为20m，该高光谱图像的大小为512×512×200。

实验仿真环境：采用MATLAB 2012R软件作为仿真工具，CPU是AMDA8-5550M，主频为2.10GHz，内存4G，操作系统为Windows 7旗舰版。

在两个不同的采样率下进行实验：第一种情况，空域的两个维度采样率S₁＝0.5，S₂＝0.5，谱域的采样率设为S₃＝0.25，则总的采样率S'＝S₁·S₂·S₃＝0.5×0.5×0.25＝6.25％；第二种情况，空域的两个维度采样率S₁'＝0.5，S'₂＝0.5，谱域的采样率设为S₃'＝0.5，则总的采样率S”＝S₁'·S'₂·S₃'＝0.5×0.5×0.5＝12.5％。实验中所用到的观测矩阵均为高斯随机矩阵，所用到的字典均为DCT字典。

2)仿真内容

仿真1：在采样率6.25％下，采用本发明方法对高光谱图像Moffett进行高光谱压缩成像仿真实验，实验结果如图2所示，其中：

图2(a)是高光谱图像Moffett第35谱段的原始图像；

图2(b)是采用本发明方法的重构图像，其PSNR＝28.53dB，该谱段的重构效果最好，PSNR最高；

图2(c)是高光谱图像Moffett第144谱段的原始图像；

图2(d)是采用本发明方法的重构图像，其PSNR＝19.58dB，该谱段的重构效果最差，PSNR最低。

仿真2：在采样率12.5％下，采用本发明方法对高光谱图像Moffett进行高光谱压缩成像仿真实验，实验结果如图3所示，其中：

图3(a)是高光谱图像Moffett第14谱段的原始图像；

图3(b)是采用本发明方法的重构图像，其PSNR 29.04dB，该谱段的重构效果最好，PSNR最高；

图3(c)是高光谱图像Moffett第165谱段的原始图像；

图3(d)是采用本发明方法的重构图像，其PSNR＝21.64dB，该谱段的重构效果最差，PSNR最低。

从图2和图3可以看出，本发明方法可以精确的恢复原高光谱图像。

仿真3：在两种采样率6.25％和12.5％下，分别采用本发明方法与现有的仅在空域采样的OMP以及仅在谱域采样的OMP方法对高光谱图像Moffett进行高光谱压缩成像仿真实验，实验结果如表1所示。

表1两种采样率下本发明方法和OMP方法的实验对比

仿真4：在两种采样率6.25％和12.5％下，分别采用本发明方法与现有的仅在空域采样的BP以及仅在谱域采样的BP方法对高光谱图像Moffett进行高光谱压缩成像仿真实验，实验结果如表2所示。

表2两种采样率下本发明方法和BP方法的实验对比

从表1和表2可以看出，随着采样率的提高，不同方法的PSNR都有所升高，在同一个采样率下，同仅在空域采样的OMP、BP方法对比，本发明方法的PSNR较高，恢复效果较好。同时还可以看出，本发明方法大大缩短了实验运行时间，进一步说明了本发明方法在减小算法复杂度上的有效性。

Claims (2)

1.一种基于三维张量压缩感知的高光谱压缩成像方法，包括如下步骤：

(1)输入高光谱图像，将其表示为三维张量Z，其中，I₁、I₂、I₃分别为高光谱图像三个维度的大小；

(2)设三个维度的采样率分别为S₁、S₂、S₃，构造一个具有克罗内克结构的观测矩阵Φ：

$\mathrm{\&Phi;}={\mathrm{\&Phi;}}_3\&CircleTimes;{\mathrm{\&Phi;}}_2\&CircleTimes;{\mathrm{\&Phi;}}_1,$

其中，Φ₁、Φ₂、Φ₃是大小分别为J₁×I₁、J₂×I₂和J₃×I₃的高斯随机矩阵，作为三个维度的观测矩阵，第i个维度的观测矩阵Φ_i的行数J_i由第i个维度的采样率S_i和高光谱图像第i个维度的大小I_i来决定，即J_i＝S_iI_i，i＝1,2,3， 表示矩阵的克罗内克积；

(3)利用高光谱图像Z和三个维度的观测矩阵Φ₁、Φ₂、Φ₃，得到低维的测量值M：

M＝Z×₁Φ₁×₂Φ₂×₃Φ₃，

其中，×_i表示高光谱图像Z与第i个维度上的观测矩阵Φ_i的张量i-模乘积；

(4)构造一个具有克罗内克结构的字典D：

$D=D_3\&CircleTimes;D_2\&CircleTimes;D_1$

其中，D₁、D₂、D₃是大小分别为I₁×I₁、I₂×I₂和I₃×I₃的三个维度上的字典，三个维度的字典均取为离散余弦字典；

(5)根据三个维度的观测矩阵Φ₁、Φ₂、Φ₃和三个维度的字典D₁、D₂、D₃，计算三个维度上的压缩感知矩阵：

Q₁＝Φ₁D₁，

Q₂＝Φ₂D₂，

Q₃＝Φ₃D₃，

得到压缩感知矩阵 $Q=Q_3\&CircleTimes;Q_2\&CircleTimes;Q_1,$ 其中，

(6)根据测量值M和各个维度上的压缩感知矩阵Q₁、Q₂、Q₃，利用张量正交匹配追踪算法求解下式，得出稀疏系数张量β：

M＝β×₁Q₁×₂Q₂×₃Q₃，

其中，×_i表示稀疏系数张量β和第i个维度上的压缩感知矩阵Q_i的张量i-模乘；

(7)根据稀疏系数β和三个维度上的字典D₁、D₂、D₃，得到重构出的原高光谱图像

$\underset{\&OverBar;}{\overset{^}{Z}}=\underset{\&OverBar;}{\mathrm{\&beta;}}{\&times;}_1{D}_1{\&times;}_2{D}_2{\&times;}_3{D}_3,$

其中，×_i表示稀疏系数张量β和第i个维度上的字典D_i的张量i-模乘，

2.根据权利要求1所述的方法，其中所述步骤6中利用张量正交匹配追踪算法求解稀疏系数β，按如下步骤进行：

2a)设初始迭代次数k＝1，残差R ₁＝M，索引集稀疏度为T＝1000，门限值ε＝0.05；

2b)计算第k次循环的残差R _k与三个维度上的压缩感知矩阵Q₁、Q₂、Q₃中未选中的列向量的乘积H：

$\underset{\&OverBar;}{H}={\underset{\&OverBar;}{R}}_k{\&times;}_1{Q}_1^T(:,i_1){\&times;}_2{Q}_2^T(:,i_2){\&times;}_3{Q}_3^T(:,i_3),$

其中，i₁、i₂、i₃分别表示每个感知矩阵中未选中的列向量的索引，( )^T表示矩阵的转置，i₁＝1,2,...,I₁，i₂＝1,2,...,I₂，i₃＝1,2,...,I₃，

2c)对乘积H各个维度上的值取绝对值，找出各个维度上绝对值最大的H所对应的感知矩阵Q₁、Q₂、Q₃中的列向量的索引i₁ ^k、i₂ ^k、i₃ ^k，然后把索引i₁ ^k、i₂ ^k、i₃ ^k添加到索引集Λ_i中，即 ${\mathrm{\&Lambda;}}_i={\mathrm{\&Lambda;}}_i\&cup;\left[i_i^k\right];$

2d)用当前索引集Λ_i所对应的压缩感知矩阵Q_i的列向量构成子矩阵：P_i＝Q_i(:,Λ_i)；

2e)利用递归的算法计算测量值M在当前的子矩阵P_i下的稀疏系数β：

$\underset{\&OverBar;}{\mathrm{\&beta;}}=\underset{\underset{\&OverBar;}{\mathrm{\&beta;}}}{\arg \min }{\left|\right|\underset{\&OverBar;}{M}-\underset{\&OverBar;}{\mathrm{\&beta;}}{\&times;}_1{P}_1{\&times;}_2{P}_2{\&times;}_3{P}_3\left|\right|}_2^2;$

其中，稀疏系数β满足||β||₀≤T，|| ||₀表示求张量的0范数；

2f)迭代次数加1，即k＝k+1；

2g)利用下式更新残差R _k：

R _k＝M-β×₁P₁×₂P₂×₃P₃；

2h)判断是否满足|Λ₁||Λ₂||Λ₃|＞T，或者||R _k||_F＜ε，其中，|Λ_i|表示索引集Λ_i包含的元素数，|| ||_F表示Frobenius范数，如果不满足，重复步骤2b)到步骤2g)，如果满足，停止迭代，得出稀疏系数β。