罗 琼 韩 志 陈希爱 唐延东

1(中国科学院沈阳自动化研究所机器人学国家重点实验室 辽宁 沈阳 110016) 2(中国科学院机器人与智能制造创新研究院 辽宁 沈阳 110016) 3(中国科学院大学 北京 100049)

0 引 言

真实场景的辐射通常遍及一个很广泛的谱段范围。高光谱图像(HSI)包含特定场景各种谱段强度的积分,相比于其他图像,HSI图像能够为真实的场景提供更多有用的信息,因而有利于计算机视觉中的很多任务,比如分类[1-3]、分割[4-5]、检测[6-7]、检索[8]。

然而,在实际情况下,HSI总会受到传感器精度、光照条件或者大气环境所引起的噪声。噪声的存在会严重影响高光谱图像后续任务的性能。因此,HSI的降噪是一个备受关注的问题。

局部平滑先验是HSI去噪任务中使用最广泛的先验之一。空间局部平滑先验是指相似的目标或者场景通常与相似的谱波相邻分布;而光谱平滑性是指HSI的相邻光谱通常是用相似的传感器参数收集而导致具有相似的值。因此,HSI所持有的局部平滑先验等同于沿着HSI的空间和谱段维度的梯度图的稀疏性,这也就是我们平时所声称的3D全变分正则或者3DTV。3DTV正则已经在各种HSI处理问题中大放异彩[9-11],一些甚至获得了最先进的性能。

即使被成功地应用HSI的各种任务,3DTV正则仍然没有充分考虑到HSI的梯度图背后的深刻的稀疏性结构知识。具体来说,遍及梯度图所有频谱的稀疏性被隐式地假设为相似的和不相关的。然而,这一点被Peng等[12]指出背离了事实。一方面,梯度图的不同谱段并不是同分布的;另一方面,HSI梯度图的不同谱段有明显的相关性,这一点源于原始HSI的波段相关性。这种与实际HSI的先验知识偏差使得这种有用的正则化项的能力仍有很大的空间需要进一步加强。因此,Peng等提出了增强版的3DTV正则E-3DTV以解决常规3DTV的不足。

进一步,一个基于E-3DTV的去噪模型被Peng等提出,他们采用1范数来拟合实际噪声。尽管凭借E-3DTV的优势取得了令人印象深刻的效果,然而,这个模型对于噪声的处理过于粗糙,1范数的噪声建模通常只能有效地处理拉普拉斯噪声。而HSI中包含的实际噪声是未知而复杂的,只针对特定类型的噪声分布建模只能获得次优的结果。

为了弥补单一范数的噪声拟合,已经有一些工作在这方面做了尝试。其中Wen等[13]提出用p范数来替代流行的1范数去拟合噪声。相比于只适合处理拉普拉斯噪声的1范数,p范数有一个更加宽广的适应范围,它把普遍应用的1范数和2范数视为其中的一种特例,并能在其中取得较好的均衡,尤其在具有非常严重的拖尾分布的高度脉冲噪声情形。令p<1可能是一种更加明智的选择。

因此,本文针对HSI中多变而复杂的实际噪声提出一种基于p范数的E-3DTV去噪模型,以提高模型在实际噪声场合中的适应能力。

1 符号说明和3DTV正则

(1)

代表在位置(i,j,k)处沿着空间高度、空间宽度和光谱维度三种不同的操作。三种在HSI上不同模式的差分运算可被计算如下:

(2)

(3)

这个线性操作编码了HSIX∈Rhw×s和它的梯度图Gn∈Rhw×s之间的关系。

在梯度上最普遍使用的稀疏度量,也就是1范数标记为:

通过这个稀疏度量,在HSI中广泛使用的局部平滑正则项,也就是所谓的3DTV[14-15]正则构造如下:

2 E-3DTV正则

为了避免3DTV中稀疏性独立同分布的假设,Peng等提出了E-3DTV。与常规3DTV相似,E-3DTV也是在长度、宽度、波段三个维度上操作。不同的是,E-3DTV的稀疏约束不施加在梯度图本身,而是施加在梯度图的子空间基上。每一组基由梯度图向量的线性组合来获得:GnVn。其中Vn是尺寸为s×r的转换矩阵,其中r<

如前所述,对于一个给定的HSI,其梯度图Gn的稀疏性并非是独立而是有相关性的,这意味Gn也可以表达为如下的矩阵分解形式:

Un=GnVn

(8)

E-3DTV正则每一个Gn(n=1,2,3)如下:

式中:Vn可以看作将不同梯度图关联在一起的一个系数矩阵。

可以看出,E-3DTV实际上度量的是Gn经过线性转换后的稀疏性,而转换矩阵Vn由输入数据自动决策。相比于常规的3DTV,E-3DTV在度量梯度图稀疏性同时考虑了梯度之间的相关性,充分利用了在Gn存在一些稀疏的低秩基,但是能用一个系数矩阵Vn表达Gn所有信息的事实,并且当HSI被破坏时,这些基的计算比梯度图本身更稳定。

由于式(9)不容易被直接解出。我们重写式(9)为如下等价形式:

Un∈Rhw×rVn∈Rs×r

通过沿着梯度图的空间高度,宽度和光谱维度应用式(11)的稀疏度量,E-3DTV正则项有如下表达式:

它有如下等价形式:

Un∈Rhw×r,Vn∈Rs×r,n=1,2,3

3 基于p范数的E-3DTV去噪模型

去噪任务是将干净的数据和噪声从带噪声的数据中分开。在许多降噪任务中,通常假定噪声为高斯噪声,因此经常利用干净数据与噪声数据之间的2范数作为损失函数。但是,基于最小二乘的2范数对观测值中的异常值高度敏感,当测量结果包含较大的误差或脉冲噪声时,1范数损失函数与2范数损失函数相比,可以显著改善性能[16-17],因此,Peng等提出基于1范数的E-3DTV去噪模型如下:

s.t.Y=X+E

式中:Y是观测图像;E是观测图像中包含的噪声;τ是权衡因子。

然而,在实际的高光谱场景中,噪声的类型远不止高斯分布和拉普拉斯分布,例如稀疏噪声、条纹噪声、死线、像素缺失等[18-19]。我们因此提出基于p范数的E-3DTV去噪模型,其数学表达式如下:

s.t.Y=X+E,0≤p≤2

p范数是1范数与2范数一种推广,它可以通过调节p值的大小而改变对异常值的敏感性,以适应复杂多变的实际噪声。例如,在没有异常值或异常值较小的情况下,可能要求p≥1,而对于具有非常严重的拖尾分布的高度脉冲噪声,可能要求p≤1。而固定类型的1范数与2范数都不具备这种泛化性。

为了对p范数有效地求解,Wen等通过将p范数函数的近端算子合并到增强拉格朗日方法的框架中,有效地解决p范数中的凸和非凸情况,并针对凸和非凸情况作了收敛性分析,本文可以直接利用其结论。考虑罚数为η的函数g(x):x∈Rm的近端操作的一般形式如下:

情形1:当p=0,近端操作成为众所周知的硬阈值操作:

式中:ti是向量t的第i个元素,i=1,2,…,m。

情形2:当0

情形3:当p=1。近端操作有闭式解:

proxg,η(t)i=Sa/η(t)i=

sign(ti)max{|ti|-a/η,0}

(18)

对于i=1,2,…,m,Sa:Rm→Rm是众所周知的软阈值或者收缩算子。

情形4:当1

proxg,η(t)i=sign(ti)zi

(19)

式中:zi是式(20)的解。

h2(z)=apzp-1+ηz-η|ti|=0z≥0

(20)

注意对于z≥0,h2(z)是一个递增的非凸函数,并且当ti≠0时,h2(z)<0,h2(|ti|)>0。因此,当ti≠0,式(20)满足0

Wen等论证了求解噪声E的子问题即可看成近端操作形式中的一种。接下来,我们将使用ADMM算法以及上述公式对式(14)求解。

4 ADMM算法求解

基于式(14),我们可以有如下等价表达式:

s.t.Y=X+E,

Vn∈Rs×r,Un∈Rhw×r,n=1,2,3

求解这个模型等价于最小化其增广的拉格朗日函数如:

式中:Mn(n=1,2,3)和Γ是拉格朗日乘子;μ是一个大于0的标量。

在ADMM框架下,我们需要交替固定其他向量并优化涉及式(23)中的每一个变量。

更新X。提取式(23)所有关于X的项,我们能获得如下子问题:

优化式(24)等价于求解如下线性系统:

(26)

更新Un,n=1,2,3。提取式(23)关于Un的所有项,我们可以得到:

这个子问题能够使用众所周知的软阈值算子求解[23]。

式中:x∈R且Δ>0是阈值。然后式(28)有如下闭式解:

更新Vn,n=1,2,3。提取式(23)关于Vn的全部项,有如下子问题:

这个子问题的全局解有闭式解[21]:

式中:B、D、C是执行svd分解后的结果。

更新E。提取式(23)所有关于E的项,关于E的子问题近端操作形式中的一种[13],它可以通过如下方式有效的解决:

(32)

基于ADMM准则,乘子的更新方式如下:

完整的ADMM求解过程如算法1所示。算法功能流程如图1所示。

图1 算法功能流程

算法1基于p范数的E-3DTV去噪模型

输出:fold(X)∈Rh×w×s。

初始化:初始化X,E,Un,Vn,Mn,Γ。

循环:

1.通过式(26)和式(32)更新X,E

2.通过式(27)和式(31)更新Un,Vn

3.通过式(33)更新Mn,Γ

4.i=i+1

5.核对收敛准则

或者i=50

结束循环

5 HSI图像去噪实验

这个部分,我们做了广泛的实验论证本文方法的性能,我们把本文方法与4种经典的HSI去噪方法相比较。它们分别是CWM[24]、LRMR[25]、LRTV[26]、E-3DTV[12]。为了全面地评价实验结果,我们采用三种量化指标,它们分别是PSNR、SSIM[27]、ERGAS[28]。PSNR和SSIM是两种经典的基于空间的图像恢复评价指标,而ERGAS是基于光谱的评价指标。PSNR和SSIM越大越好,而ERGAS越小越好。

5.1 实验设置

案例1:每一幅图像都加入零均值、0.1方差的高斯噪声。

案例2:首先加入类似于案例1中的高斯噪声到每一幅图像上,然后从第91幅至第130幅加入条纹噪声。

案例3:每一幅HSI加入零高斯噪声和脉冲噪声。高斯噪声方差设为0.075,脉冲噪声的百分比设为0.5。

案例4:加入像案例3一样的高斯脉冲混合噪声,加入像案例2一样的条纹噪声。

5.2 实验结果展示

视觉质量方面,我们在图2展示了通过不同方法在四种不同噪声案例下的恢复结果,其中案例1展示的图像是第56幅,案例2展示的图像是第112幅,案例3展示的图像是168幅,案例4展示的图像是第224幅。可以明显地看到E-3DTV相较于其他方法不仅能使恢复的图像残留更少的噪声,而且还能保持更加锐利的细节。而我们采用基于p范数的E-3DTV模型能在此基础上获得性能上的进一步提升。这得益于p范数处理不同噪声的灵活性。

图2 不同方法在不同噪声情形下的比较

在量化比较方面,我们把各种方法恢复结果的三种指标展示在表1至表4。E-3DTV模型在其他三种竞争方法面前仍然保持较大优势。这侧面印证了E-3DTV的优越性,然而,我们的方法在各个指标全面超越了E-3DTV,印证了本文方法的意义和有效性。

表1 不同方法在案例1噪声情形下的比较

表2 不同方法在案例2噪声情形下的比较

表3 不同方法在案例3噪声情形下的比较

表4 不同方法在案例4声情形下的比较

由于本文重点考虑的是p范数相较于单纯的1范数和2范数对E-3DTV去噪模型带来的灵活性。我们重点观察了p参数的设置对实验结果带来的影响。以案例1和案例3为例,在保持其他参数设置不变的情况下,我们列举当p<1、p=1、2>p>1、p=2时对结果的影响在表5和表6,我们加粗了最好的结果。可以看到对于不同的噪声采用不同p值对结果有较大影响,因而单纯地用1范数或者2范数建模噪声是不足够的,此外,最好的效果均不在p刚好为1或刚好为2时取得,进一步说明了p范数建模的意义。

表5 本文的方法在案例1噪声情形下不同p值的比较

表6 本文方法在案例3噪声情形下不同p值的比较

6 结 语

本文提出新的基于p范数的E-3DTV模型用于高光谱图像去噪,得益于p范数在处理不同噪声时的灵活性,本文方法具有比原E-3DTV模型更强大的噪声去除能力和细节恢复能力。下一步的研究会考虑把p范数用于约束梯度图的基矩阵Un,以期望获得适应性更强的HSI去噪模型。

此外,由于HSI采集通常有数百个频谱,这使得HSI需要较大的存储空间。在将它们发送到地面基站时,这将引起低效率和高成本的严重问题。因此,有必要针对HSI压缩感知设计有效的技术,以满足实时传输的需要。由于p范数和E-3DTV正则的结合展现出了良好的重构性能,把这个方法拓展在压缩感知模型上也将作为我们未来研究的方向。