一种基于E-3DTV正则的高光谱图像修复方法与流程

本发明涉及对传统的3dtv的一种增强改进，使得新的正则可以兼顾高光谱图梯度图像的结构的相关性和稀疏性，而不同于传统的3dtv只能刻画梯度图的稀疏性，而忽略了其本质的相关性。具体涉及针对高光谱图像数据的去噪和压缩重构方法。

背景技术：

卫星上携带的光谱成像仪可以在同一时刻对同一地物形成许多连续谱段下的图像，因此在光谱维度每一个像素都会形成一条独特的光谱曲线，这是区分不同物质的一个显著特征。因此，高光谱图像相比传统图像在光谱维度具有更为丰富的信息，因此被广泛应用到地物识别，城市规划，军事勘测等应用问题中。

然而，地端到卫星之间的成像距离通常很长，成像中往往会受到大气气溶胶、光照、水气等因素的污染，使得高光谱图像中往往混杂大量噪音，从而使得部分高光谱图像的结构纹理模糊，物质的光谱曲线被严重破坏，这些现象都严重影响了针对高光谱图像进行应用与处理的后续目标，因此针对高光谱图像的修复具有非常重要的实用价值。除此之外，高光谱图像数据需要被传输到地面来进行进一步的分析与处理任务，由于高光谱图像通常蕴含高谱段的丰富信息，其数据量往往非常巨大。如何有效设计传输重构算法来快速有效的解决高光谱图像传输和储存问题，是针对高光谱图像另一个具有重要意义而挑战性的任务。由于近年来压缩重构理论的发展，压缩重构逐渐倾向于突破奈奎斯特采样的限制，所需传输的采样子个数期待大大减少，高光谱的传输问题也转变为利用压缩重构获得好的重构效果的问题。

以上高光谱的去噪和压缩重构问题为最具代表性的两类高光谱图像修复问题。

高光谱图像修复技术主要基于构建高光谱图像先验信息来进行恢复。根据其先验只是的分类，可以将高光谱图像修复方法分为三类。基于高光谱图像光谱之间的相关性先验知识构建的方法，基于高光谱图像的局部平滑性先验构造的方法，与基于高光谱图非局部自相似性先验构造的方法。实际方法可能会综合采用以上两类及多类先验来构造高光谱图像修复算法。

其中基于光谱相关性先验知识构造的方法，是指利用高光谱图像在光谱维度上的相关性，通过构造一组基及在基下的表达系数来表示高光谱图像。比较常见的方法包括低秩类方法，包括低秩矩阵和低秩张量类方法，以及字典表达类方法。此类方法均为基于高光谱的主体信息可以通过稀疏表达而噪音不能被稀疏表达的基本礼节来实现噪音和干净数据分离的效果。

基于图像的非局部自相似先验类的方法，是在高光谱图像块层面上进行的。将图片分成许多小块后，有许多相似的小块出现。这样就可以将图片的小块分成许多相似的小组，由于每个小组的小块都比较相似，可利用低秩和字典学习等方法对每组块进行共有信息的抽取，从而实现对每组块共同的修复，最终将所有修复组进行合并就可以最终得到修复后的高光谱图。

基于图像的局部平滑性先验类的方法，是基于如下统计规律而构建：相邻的像素一般是同类的物质，因此相邻像素的灰度值差异很小，在不同物质的连接处，像素的灰度值的差异很大，这种灰度值的差异反映了图像的主体结构，即是图像的边缘结构。可以用图片的梯度图的稀疏性来刻画这种灰度值差异的先验，因此一般用全变差正则(tv)来刻画这种先验信息。当对图像使用梯度或者差分算子操作时，梯度图可以显著的提取图像边缘结构，这对高光谱图像的修复会起到很好的作用。

现有的技术尽管是基于对高光谱图的先验知识来构造模型，但其先验并未能充分精细刻画高光谱图像的信息。如基于高光谱图的相似性的方法忽略了高光谱图的局部平滑性，而基于图像的局部平滑性的又不能很好的刻画光谱的相似性，尽管也有一些基于两类甚至三类先验融合的方法，但这些方法的先验知识通常采用简单的叠加，并不能充分的刻画高光谱图丰富的先验结构耦合的信息。因此寻找能够充分表达高光谱图信息先验及其先验编码方式对于推动高光谱图修复方法的性能具有非常重要的作用与意义。

技术实现要素：

本发明的目的在于提供一种基于增强三维全变差正则(enhanced3-dimensionaltotalvariation,简称e-3dtv)的高光谱图像(hyperspectralimage)修复方法，对高光谱图像数据进行去噪和压缩重构，通过对传统3dtv的增强改进，使得新的正则可以兼顾梯度图像的结构相关性和稀疏性，克服传统的3dtv只能刻画梯度图的稀疏性，而忽略了相关性的缺陷。

为了实现上述目的，本发明采用的技术方案是：

一种基于e-3dtv正则的高光谱图像修复方法，包括如下步骤：

步骤s1：将含噪音的原三维高光谱数据沿着光谱维度展开成矩阵得到y∈r^mn×t，其中m,n代表每个谱段下图像的长和宽，t代表高光谱的光谱数，初始化噪音项ε∈r^m×n×t的矩阵形式e∈r^mn×t和待修复的高光谱数据的矩阵表示x∈r^mn×t以及admm框架下的其他模型变量及参数；

步骤s2：对待修复的高光谱数据沿着水平、竖直、光谱三个维度做差分操作，差分算子记为得到的三个不同方向的梯度图，并沿着光谱维度展开成矩阵，记水平、竖直、光谱方向的梯度图分别为i＝1,2,3；

步骤s3：对三个方向的梯度图矩阵分别做低秩uv分解，其中基矩阵u反应反映了梯度图的共有的结构和特征信息，v是系数矩阵，它反映了梯度图矩阵的每一列由基矩阵表达时的系数；对梯度图的基矩阵u加稀疏性约束，得到e-3dtv正则；

步骤s4：根据步骤s3，对待修复的x加e-3dtv正则，写出优化模型，利用admm框架进行迭代求解，到迭代稳定，获得修复图像x和噪音e。

所述步骤s1中，admm框架下的其他模型变量及参数包括admm中需要用到的上升乘子矩阵m和上升乘子常数的初始值μ和每次上升的倍率ρ，m＝0，μ的经验取值在[1e^-3,1e^-2]之间选择，ρ的经验取值在[1.05,1.8]之间选择。

所述步骤s2中三个方向的差分图计算公式定义如下所示：

其中，i,j,k分别为高光谱图像对应的沿空间长度、宽度与谱段方向的坐标位置，表示原高光谱图像在位置(i,j,k)的灰度值，s＝1,2,3分别表示的是原高光谱的在位置的灰度值表示沿着空间长度、宽度与谱段对高光谱图进行差分运算之后得到的梯度图在位置(i,j,k)的取值。分别表示对应高光谱图像的(i,j,k)位置沿着空间长度、宽度与谱段方向分别平移一个位置的像素灰度值。为了保证差分图和待修复图尺寸一致便于计算，对差分图做zero-padding处理，将梯度图沿着光谱维度展开为矩阵，得到如下的公式：

其中unfold表示沿着光谱方向展开。

所述步骤s3中，e-3dtv正则的形式如下：

其中，ui为梯度矩阵的基，由于ui中包含了梯度矩阵的共有结构或者特征，因此给ui加上了一范数‖ui‖1来约束基的稀疏性；vi为梯度矩阵在基矩阵下的系数矩阵。为了保证上述优化问题解的唯一性，给系数矩阵加上vivi^t＝i的限制，这里i表示单位矩阵；r为低秩uv分解中的秩，r<min(mn,t)，蕴含了重建矩阵的低秩性信息。

所述步骤s4中，建立的优化模型如下所示：

s.t.y＝x+e(4)

其中τ为调节参数，作用在于协调e-3dtv正则和误差项的作用折中，y为输入的带噪音数据，e为噪音项，x为待修复的高光谱图数据，等式(4)展开即可得到下面的公式(5)：

所述步骤s4中，利用admm框架进行迭代求解的步骤如下：

s6.1)给出公式(5)的增广拉格朗日函数：

其中，mi,γ是拉格朗日乘子，μ是上升乘子常数，<.，.>表示内积运算，表示矩阵的f范数；

s6.2)建立交替方向乘子法的迭代格式与终止条件：

μ^k+1＝ρμ^k(12)

迭代终止条件为：

s6.3)对问题(6)、(7)、(8)、(9)进行求解，给出迭代的具体算式；

s6.4)设置迭代的初始值为：e⁰＝0，由奇异值分解方法作用在y上产生；

s6.5)进行公式(6)-(12)的迭代运算，直到迭代满足终止条件。

通过所述步骤s1-s4实现高光谱图像去噪；在所述步骤s1-s4的基础上，增加一步由少量采样数据y来计算重构数据的过程，实现高光谱图像压缩重构。

所述高光谱图像压缩重构的具体流程如下：

s8.1)压缩重构问题简述

构造采样算子ψ＝d.h.p，其中d是随机下采样算子，h是随机置换矩阵，p是walshhadamard变换算子，固定采样算子ψ之后，高光谱数据z∈r^mn×t经ψ的作用会得到采样过后的观测子y∈r^l,l＜＜mnt，压缩重构的目的就是通过观测y重构出z。压缩重构蕴含内在等式约束为y＝ψ(z)；

s8.2)e-3dtv正则用以高光谱图像压缩重构的计算方法；

s8.2.1)e-3dtv正则用以高光谱图像压缩重构的优化模型如下

将z分解成待修复高光谱干净数据项x和噪音项e,对x干净数据加e-3dtv正则，得到优化模型如下所示：

s.t.y＝ψ(z),z＝x+e(17)

其中τ为调节参数，作用在于协调e-3dtv正则和误差项的作用折中。，将(17)式展开即得：

s.t.y＝ψ(z)，z＝x+e(18)

s8.2.2)给出公式(18)的增广拉格朗日函数

其中，mi,γ1,γ2是拉格朗日乘子，μ是一个大于0的数；

s8.2.3)建立交替方向乘子法的迭代格式与终止条件：

μ^k+1＝ρμ^k(27)

迭代终止条件为：

s8.2.4)对问题(19)、(20)、(21)、(22)、(23)进行求解，给出迭代的具体算式；

s8.2.5)设置迭代的初始值为：e⁰＝0，由奇异值分解方法作用在ψ^*(y)上产生，其中ψ^*是ψ的共轭算子；

s8.2.6)进行(19)-(27)的迭代运算，直到迭代满足终止条件。

传统的3dtv正则仅考虑了高光谱梯度图的结构稀疏性，并没有更深层次的对梯度图先验知识进行更深入的挖掘。而本发明所提的e-3dtv正则利用子空间学习的方法，对于梯度图的基进行稀疏约束，从而很好的刻画了高光谱梯度图结构的相关稀疏性，这一点是之前tv类方法从未考虑到的。首先，利用差分算子对高光谱图的三个维度进行差分操作，得到原高光谱图像的梯度图，空间上的梯度图可以反映出图片的边缘结构，光谱维度的梯度图可以反映出相邻谱段的相似性。通过算法的不断迭代，梯度图的基及其梯度图在基图下的表示矩阵和噪音矩阵均不断更新，同时基、系数矩阵及其噪音同时作用，不断反解出待修复的高光谱图。本发明通过对梯度图的认知，对3dtv正则进行本质的改进，从而在高光谱图像上大幅度的提升针对此类图像的修复效果。

附图说明

图1为本发明的流程图。

图2为本发明e-3dtv正则的阐释示意图。

图3为本发明在梯度图的基上进行修复能够带来鲁棒性效果的示意图。

图4为本发明基于e-3dtv正则在indianpines高光谱数据上band#220的去噪效果。其中(a)为原图，(b)为噪音，(c)为修复图。

图5为本发明基于e-3dtv正则在dcmall高光谱数据上band#120的去噪效果。其中(a)为原图，(b)为噪音，(c)为修复图。

图6为本发明基于e-3dtv正则在真实的indianpines高光谱数据上band#108的去噪效果。其中(a)为噪音图，(b)为修复图。

图7为本发明基于e-3dtv正则在真实的urban高光谱数据上band#108的去噪效果。其中(a)为噪音图，(b)为修复图。

图8为本发明基于e-3dtv正则在真实的terrian高光谱数据上band#207的去噪效果。其中(a)为噪音图，(b)为修复图。

图9为本发明基于e-3dtv正则在真实的lowal高光谱数据上band#109的去噪效果。其中(a)为噪音图，(b)为修复图。

图10为本发明基于e-3dtv正则在采样率为5％的情况下dcmall高光谱数据上band#5的重构效果。其中(a)为原图，(b)为重构图。

图11为本发明基于e-3dtv正则在采样率为5％的情况下urban高光谱数据上band#160的重构效果。其中(a)为原图，(b)为重构图。

图12为本发明基于e-3dtv正则在采样率为5％的情况下lowal高光谱数据上band#160的重构效果。其中(a)为原图，(b)为重构图。

图13为本发明基于e-3dtv正则在采样率为5％的情况下moffett高光谱数据上band#80的重构效果。其中(a)为原图，(b)为重构图。

具体实施方式

下面结合附图和实施例详细说明本发明的实施方式。

如图1所示，本发明一种基于e-3dtv正则的高光谱图像修复方法，包括如下步骤：

步骤s1：将含噪音的原三维高光谱数据沿着光谱维度展开成矩阵得到y∈r^mn×t，其中m,n代表每个谱段下图像的长和宽，t代表高光谱的光谱数，初始化噪音项ε∈r^m×n×t的矩阵形式e∈r^mn×t和待修复的高光谱数据的矩阵表示x∈r^mn×t以及admm框架下的其他模型变量及参数，其他模型变量及参数包括admm中需要用到的上升乘子矩阵m和上升乘子常数的初始值μ和每次上升的倍率ρ，m＝0，μ的经验取值在[1e^-3,1e^-2]之间选择，ρ的经验取值在[1.05,1.8]之间选择；

步骤s2：对待修复的高光谱数据沿着水平、竖直、光谱三个维度做差分操作，得到的三个不同方向的梯度图，并沿着光谱维度展开成矩阵，记水平、竖直、光谱方向的梯度图分别为i＝1,2,3；

三个方向的差分图计算公式定义如下所示：

其中，i,j,k分别为高光谱图像对应的沿空间长度、宽度与谱段方向的坐标位置，表示原高光谱图像在位置(i,j,k)的灰度值，s＝1，2，3分别表示的是原高光谱的在位置的灰度值表示沿着空间长度、宽度与谱段对高光谱图进行差分运算之后得到的梯度图在位置(i，j，k)的取值。分别表示对应高光谱图像的(i,j,k)位置沿着空间长度、宽度与谱段方向分别平移一个位置的像素灰度值。为了保证差分图和待复原图尺寸一致便于计算，对差分图做zero-padding处理，将梯度图沿着光谱维度展开为矩阵，得到如下的公式：

其中unfold表示沿着光谱方向展开。

步骤s3：对三个方向的梯度图矩阵分别做低秩uv分解，其中基矩阵u反映了梯度图的共有的结构和特征信息，v是系数矩阵，它反映了梯度图矩阵的每一列再由基矩阵表达时的系数；对梯度图的基矩阵u加稀疏性约束，得到e-3dtv正则；

e-3dtv正则的形式如下：

其中，ui为梯度矩阵的基，由于ui中包含了梯度矩阵的共有结构或者特征，因此给ui加上了一范数‖ui‖1来约束基的稀疏性；vi为梯度矩阵在基矩阵下的系数矩阵。为了保证上述优化问题解的唯一性，给系数矩阵加上vivi^t＝i的限制，这里i表示单位矩阵；r为低秩uv分解中的秩，r<min(mn,t)，蕴含了重建矩阵的低秩性信息。

由线性代数的知识可得，每一个谱段的梯度都可以通过基和系数来进行线性表示。相对传统的刻画梯度图稀疏性先验的3dtv而言，本发明的正则不是在上加一范数进行稀疏性约束，而是在基矩阵ui上加一范数进行稀疏性约束。这样做的好处是可以充分的挖掘梯度图之间的相关性，通过稀疏基矩阵ui和正交系数来表达梯度图的稀疏性和关联性，直观有效。

步骤s4：根据步骤s3，对待修复的x加e-3dtv正则，写出优化模型，利用admm框架进行迭代求解，到迭代稳定，获得已修复的图像x和噪音e。

优化模型如下所示：

s.t.y＝x+e(4)

其中τ为调节参数，作用在于协调e-3dtv正则和误差项的作用折中，y为输入的带噪音数据，e为噪音项，x为待复原的高光谱图数据，等式(4)展开即可得到下面的公式(5)：

利用admm框架进行迭代求解的步骤如下：

s6.1)给出公式(5)的增广拉格朗日函数：

其中，mi,γ是拉格朗日乘子，μ是上升乘子常数，<.,.>表示内积运算，表示矩阵的f范数；

s6.2)建立交替方向乘子法的迭代格式与终止条件：

μ^k+1＝ρμ^k(12)

迭代终止条件为：

s6.3)对问题(6)、(7)、(8)、(9)进行求解，给出迭代的具体算式；

s6.4)设置迭代的初始值为：e⁰＝0，由奇异值分解方法作用在y上产生；

s6.5)进行公式(6)-(12)的迭代运算，直到迭代满足终止条件。

步骤s5：步骤s1到s4，给出了e-3dtv正则的一般高光谱图像修复流程，具体到去噪问题中直接套用。针对压缩重构问题，需要多加一步由少量采样数据y来计算重构数据的过程，具体地，高光谱图像压缩重构的具体流程如下：

s8.1)压缩重构问题简述

构造采样算子ψ＝d.h.p，其中d是随机下采样算子，h是随机置换矩阵，p是walshhadamard变换算子，固定采样算子ψ之后，高光谱数据z∈r^mn×t经ψ的作用会得到采样过后的观测子y∈r^l,l＜＜mnt，压缩重构的目的就是通过观测y重构出z。压缩重构蕴含内在等式约束为y＝ψ(z)；

s8.2)e-3dtv正则用以高光谱图像压缩重构的计算方法；

s8.2.1)e-3dtv正则用以高光谱图像压缩重构的优化模型如下

将z分解成待复原高光谱数据项x和噪音项e,对x干净数据加e-3dtv正则，得到优化模型如下所示：

s.t.y＝ψ(z),z＝x+e(17)

其中τ为调节参数，作用在于协调e-3dtv正则和误差项的作用折中。，将(17)式展开即得：

s8.2.2)给出公式(18)的增广拉格朗日函数

其中，mi,γ1,γ2是拉格朗日乘子，μ是一个大于0的数；

s8.2.3)建立交替方向乘子法的迭代格式与终止条件：

μ^k+1＝ρμ^k(27)

迭代终止条件为：

s8.2.4)对问题(19)、(20)、(21)、(22)、(23)进行求解，给出迭代的具体算式；

1:更新噪音项e，从拉格朗日式子中，提取含e的元素，得到

2:更新基矩阵ui,i＝1,2,3，从拉格朗日式子中，提取含ui的项目，得到

3:更新表达系数vi,i＝1,2,3，从拉格朗日式子中，提取含vi的项目，得到

该问题可由矩阵的vonneumann’straceinequality得到显式解:

4:更新x，从拉格朗日式子中，提取含x的项目，得到

对以上算式求导，将含x的式子移到等式一边，可得以下线性方程组。

与去噪模型相似，可算得其显式解形式，

5:更新z，从拉格朗日式子中，提取含z的项目，得到

对(33)求导，将含z的式子移到等式一边，可得以下线性方程组

μψ^*ψz+μz＝μ(x+e)+μψ^*y+ψ^*γ1-γ2(34)

其中ψ^*是ψ的伴随算子。该线性系统可以由现成的共轭梯度算法进行求解。

s8.2.5)设置迭代的初始值为：e⁰＝0，由奇异值分解方法作用在ψ^*(y)上产生，其中ψ^*是ψ的共轭算子；

s8.2.6)进行(19)-(27)的迭代运算，直到迭代满足终止条件。

本发明的原理及效果可参考图2-图13。

图2为本发明e-3dtv正则的阐释示意图，图2中第a行的分别表示的是原高光谱图像数据及其沿着光谱维度展开后的矩阵图。第b行表示的是三个方向的差分操作示意图。第c行表示的是三个方向上得到的梯度图。第d行表示的是梯度图沿着光谱维度展开后得到的矩阵图。第e行表示的是梯度图的数值频数分布和矩阵的奇异值曲线，两条曲线表明梯度图具有很强的稀疏性和相关性，这就说明梯度图本身有很强的冗余性，这种冗余在这里就体现为相关性，而传统的3dtv只是刻画梯度图的稀疏性，并不能刻画这种相关性。第f行表示的对原图做了一个正交变化之后，可以得到一个更紧致的数据表达，这个表达实际上就是梯度矩阵的基。通过频数分布图可得知在梯度图的基上面也是具有稀疏性的。

在上面的描述，第f行得到的数据表达就是梯度图的基的证明如下：

对矩阵gi做低秩矩阵分解，得到gi＝uv^t,其中v^tv＝i。把u记为矩阵的基，v记为系数。等式两边同时乘以v^t，可以得到giv＝uv^tv＝u。

图3为本发明在梯度图的基上面进行修复能够带来鲁棒性效果的示意图，可以阐明在梯度图的基上面进行修复能够带来的鲁棒性的效果。

图3中展示的第二列是干净图像的信息，第一列是加了少许噪音的图像。第a行是原始高光谱图像；第b行是水平方向的差分图信息，第c行是光谱维度上的差分图；通过第b,c两行可知在混杂噪音情形，图像的梯度被严重破坏；第d行是对水平方向高光谱图的梯度图做主成分提取得到的前两个基图；第e行是对光谱方向上的高光谱图的梯度图做主成分提取得到的前两个基图；通过第d,e两行可以知道在噪音情况，梯度图的基仍具有一定的鲁棒性。所以在构造e-3dtv的时候，选择在基图上进行操作是合理的。

图4为展示的e-3dtv方法在仿真实验indianpines数据上band#220的去噪的效果。第一列是原图，第二列是在原图中加上高斯核稀疏噪音得到的噪音图，第三列是e-3dtv方法去噪后的结果，该结果与原图越接近，说明去噪方法越好。为了方便展示本发明方法在高光谱图像上的修复效果，对图片的局部进行了放大。

图5展示的是e-3dtv方法在仿真实验dcmall数据上band#120的去噪的效果。第一列是原图，第二列是在原图中加上高斯核稀疏噪音得到的噪音图，第三列是e-3dtv方法去噪后的结果。同样为了方便展示本发明方法在高光谱图像上的修复效果，对图片的局部进行了放大。

图6展示的是e-3dtv方法在真实实验indianpines数据上band#108的去噪的效果。由于是真实数据，没有真实的无噪音的图片，所以只能通过修复后的视觉效果来进行观察。

图7展示的是e-3dtv方法在真实实验urban数据上band#108的去噪的效果。同样为了方便展示本发明方法在高光谱图像上的修复效果，对图片的局部进行了放大。

图8和图9展示的是e-3dtv方法在真实实验terrian数据上band#207和lowal数据上band#109的去噪的效果。可以看到这两个数据集和前面的数据的结构以及真实的噪音是大不相同，而本发明方法依旧可以很好的去除噪音。同样为了方便展示本发明方法在高光谱图像上的修复效果，对图片的局部进行了放大。

图10到图13展示的是e-3dtv方法在高光谱压缩重构上的修复效果。同样为了展示方便，对修复图片的局部进行了放大。在展示的四幅图中，采样率均为5％，即本发明方法需要从原图的5％的采样数据修修复始图像。修复的图像和原始图像越相近，则修复的效果越好。其中图10展示的是dcmall数据上在band#5下的修复效果，图11展示的是urban数据上在band#160下的修复效果，图12展示的是lowal数据上在band#160下的修复效果，图12展示的是moffett数据上在band#80下的修复效果。