周 航, 苏延池, 李占山, 花昀峤
(1.吉林大学 软件学院, 长春 130012；
2.吉林大学 人工智能学院, 长春 130012；
3.吉林大学 计算机科学与技术学院, 长春 130012；
4.吉林大学 资产管理处, 长春 130012)
高光谱图像是一种高维图像, 可反映真实场景的空间信息和光谱信息.但由于图像获取过程中的不可控因素, 例如传感器的灵敏度、光子效应和校准误差等, 图像中会不可避免地存在噪声.这些噪声将严重影响图像信息的传达和后续处理, 限制了其在视觉领域的进一步应用.图像去噪作为视觉领域的预处理步骤, 旨在去除图像中的噪声成分.非局部自相似性和全局光谱相关性是高光谱图像去噪领域内的两个重要先验知识.Othman等[1]首次将二者相结合.Zhang等[2]先将高光谱图像分解为多个互相重叠的全波段三维图像块, 然后沿光谱维度将三维图像块展开成矩阵, 在低秩矩阵恢复框架下进行去噪.为进一步提升去噪性能, 文献[3]提出了低秩矩阵近似模型, 以适应不同波段噪声强度的差异, 显著提高了低秩矩阵恢复的去噪性能.此外, 不同矩阵秩的近似被集成到去噪模型中, 如γ-范数[4]和加权Schatten p范数[5]等.文献[6]基于光谱相关性提出了低秩矩阵分解方法用于图像去噪.上述方法将低秩张量展开成矩阵进行分析, 破坏了图像的内在三维结构, 不能充分利用高光谱图像的先验信息.在图像高维空间上的块匹配操作也导致这些算法的复杂度较高.
针对上述问题, 本文提出一种基于子空间表示和加权低秩张量正则化的方法(subspace representation and weighted low-rank tensor regularization, SWLRTR).首先, 通过噪声图像学习低秩正交基, 将高光谱图像投影到低维子空间中, 得到维度缩减的图像；
然后, 引入加权低秩张量正则化项和l1范数分别约束随机噪声和稀疏噪声；
最后, 针对提出的模型设计相应的求解算法.实验结果表明, 该方法在模拟噪声和真实数据集上效果均较好.
N阶张量X∈I1×I2×…×IN的F范数定义为
l1范数定义为
N阶张量X∈I1×I2×…×IN的模式n向量是一个以in为元素下标, 其他下标{i1,i2,…,iN}in全部固定不变的In维向量, 用符号记作Xi1…in-1: in+1…iN.张量的矩阵化是指将一个N阶张量重新组成一个矩阵形式的变换, 分为横向展开和纵向展开, 将N阶张量X∈I1×I2×…×IN的模式n向量沿水平方向依次排列, 称为张量的模式n横向展开, 用符号X(n)表示,X(n)∈In×(I1…In-1In+1…IN).将张量的模式n向量沿纵向依次排列, 称为张量的模式n纵向展开, 记为X(n),X(n)∈(I1…In-1In+1…IN)×In.一个N阶张量X∈I1×I2×…×IN与一个Jn×In矩阵U(n)的n模式矩阵积用符号X×nU(n)表示, 是一个I1…In-1×Jn×In+1…IN维张量, 其元素定义为
其中:j=1,2,…,Jn；
ik=1,2,…,Ik；
k=1,2,…,N.Tucker分解是一种高阶奇异值分解, 是矩阵奇异值分解向高阶的推广, 每个I1×I2×…×IN维实张量X均可分解为n模式积X=G×1U(1)×2U(2)…×NU(N), 其中U(n)是一个In×Jn半正交矩阵, 与奇异值矩阵不同, 核心张量G不取对角结构, 一般是一个满张量, 表示各模式矩阵U(n)(n=1,2,…,N)之间的相互作用.
1.1 问题描述
噪声高光谱图像Y∈n1×n2×n3(n1,n2分别为空间域上的长和宽,n3为频谱数目)可被建模为干净信号X与噪声之和, 根据噪声强度的分布特点可将噪声分为随机噪声N和稀疏噪声S.因此, 高光谱图像混合噪声去除模型可表示为
Y=X+N+S.
(1)
求解去噪问题的关键是基于高光谱图像的先验信息建立合适的正则化项.基于此, 高光谱图像的去噪模型可表示为
(2)
(3)
1.2 子空间表示
高光谱图像的每个光谱信号均可由少量光谱成分线性表示.假设高光谱图像X存在于一个k维子空间Hk中, 干净的高光谱图像可表示为X=Z×3A, 其中A∈n3×k是包含简化图像Z∈n1×n2×k的基底的正交矩阵.因此, 高光谱图像去噪模型(3)可表示为
He等[7]已经证明, 矩阵A的正交性使得A中每列对应的高光谱图像特征之间具有较大的差异, 有效避免了Y中的噪声信号在计算过程中出现在矩阵A中, 因此维度缩减后的图像Z保留了与Y中相同的噪声分布.此外, 将高光谱图像分解为模-3张量矩阵积的形式保护了结构信息的完整性, 并降低了算法的复杂度.
1.3 加权低秩张量正则化
本文用简化图像的先验信息通过组合相似的全光谱频带图像块得到一个具有明显低秩性的张量.在实际应用中, 为简化计算一般在参考图像块的邻域内寻找.假设图像块的大小为p×p×k, 相似图像块的数目为q, 则构造张量Zi的过程如下：
步骤1) 在参考图像块邻域内间隔特定步长选取图像块, 即搜索范围;
步骤2) 计算参考图像块与邻域内图像块的欧氏距离, 衡量图像块之间的相似程度, 选取欧氏距离最小的前q个图像块组成相似图像块组;
步骤3) 将相似图像块组中的三维图像块沿第3模式纵向展开成p2×k大小的二维矩阵;
步骤4) 将步骤3)得到的q个二维矩阵堆叠为p2×q×k大小的三维张量Zi.
将上述过程描述为Zi=RiZ, 其中i是参考图像块左上角像素点坐标,Ri表示在简化图像Z上提取构造Zi的一系列操作.不同的正则化项会影响去噪效果, 本文引入加权低秩张量正则化项约束Zi的低秩性.根据Zi估计其低秩近似, 即干净张量Li可描述为
(5)
由于Zi的元素之间具有较强的相关性, 其核心张量的奇异值大多数接近于零.较大的奇异值对应Zi主要的投影方向, 较小的奇异值一般大于干净张量, 本质上是噪声导致的扰动.因此应合理地对较大的奇异值施加较少的惩罚, 对较小的奇异值施加较多的惩罚.权重设置为
(6)
其中c>0是一个常数,ε是一个保证除数不为零的较小常数.基于以上分析, 本文提出的基于子空间迭代的加权低秩张量正则化高光谱图像去噪模型形式如下：
(7)
本文引入一种基于迭代最小化策略的算法求解提出的模型.模型(7)的求解过程包括两步：
求解加权低秩张量正则化项和简化图像.
2.1 求解加权低秩张量正则化项
本文从低秩张量Zi中估计张量Li, 忽略模型中与Li无关的变量, 得到如下子问题：
(8)
将Li替换为其Tucker分解的形式, 则式(8)可转化为
(9)
其中∘表示基于元素的乘积,j表示一个三阶张量的模式索引.分别求解U1,U2,U3, 令W1=Zi(1)(U3⊗U2), 其中⊗表示Kronecker积, 则
U1=PQT,
(10)
(11)
得到4个变量后, 用模-3张量矩阵积恢复最优Li.
2.2 求解简化图像Z
固定Li, 则式(7)可转化为下列最小化问题：
(12)
式(12)的增广Lagrange函数为
其中ξ{Ik}(·)是指示函数.采用交替最小化方法迭代优化S,Z,A, 求解式(13), 步骤如下.
1) 更新S：
固定Z和A, 则式(13)可被改写为
(14)
其闭形式的解为Sλ2(Y-Z×3A), 其中Sτ(x)=sgn(x)max{|x|-τ,0}是软阈值操作.
2) 更新Z：
固定S和A, 则式(13)可被改写为
(15)
式(15)是一个二次优化问题, 具有封闭解, 表示为
(16)
3) 更新A：
固定S和Z, 则式(13)可被改写为
(17)
令M=(Y-S)(3)Z(3), 假设M奇异值分解的结果为M=UΛVT, 则
(18)
2.3 迭代正则化
迭代正则化为每次迭代的输入加入一定比例的噪声信号作为对结果的扰动.该策略已被广泛应用以提升去噪效果.在下一次迭代中更新噪声图像的公式为
Yn+1=αXn+(1-α)Y,
(19)
其中α是平衡下次迭代输入的参数,n是迭代次数.子空间大小k的选择与不同的数据集和噪声强度有关, 本文根据HySime[8]初始化k.随着每次迭代噪声强度越来越弱, 得到的去噪图像越来越接近真实图像,k应逐渐增大使维度缩减后的图像Z保留更多的干净图像信息.因此,k随迭代次数变化的公式为
k=k+β×n,
(20)
其中β是一个决定步长的常数.从而每次迭代后An+1可从输入图像中提取更多信息.去噪过程对应的流程如图1所示.
图1 去噪过程Fig.1 Denoising process
2.4 复杂度分析
其中N1是估计S,Z,A时达到收敛所需的迭代次数,N是去噪过程整体的迭代次数.
下面用仿真实验和真实数据集上的实验说明本文方法的有效性.将本文方法与LRTDTV[9],NMoG[10],L1HyMixDe[11],SDeCNN[12],3DLogTNN[13],SNLRSF[14]高光谱图像去噪方法进行比较, 对比方法的参数设置参照源代码和原文献中的要求.所有测试高光谱图像在去噪前都归一化到[0,1]波段, 所有算法在16 GB RAM的Intel(R) Core(TM) i7-8700 CPU 3.20 GHz环境下, 用MATLAB R2018a实现.
3.1 模拟数据实验
在数据集Washington DC(WDC),Pavia University(PaviaU)和Pavia Center(PaviaC)上进行仿真实验.真实场景下的高光谱图像通常被不同类型的噪声污染, 包括高斯噪声、脉冲噪声和条纹等.为模拟真实高光谱图像的噪声分布, 将不同类型的噪声以不同的形式组合加入到高光谱图像中：
第一种是具有相同分布的高斯噪声加入到3个模拟数据集的每个频带中, 其均值为0, 方差为[0.1,0.2]；
第二种是具有不同强度的高斯噪声加入到3个模拟数据集的每个频带中, 其均值为0, 噪声方差σ为[0.1,0.2]；
第三种将与第二种情况相同的高斯噪声加入到3个高光谱数据集的每个频带中, 此外, 随机选取20个频带加入20%的脉冲噪声；
第四种是在原始图像中加入与第三种情况相同的高斯噪声和脉冲噪声.本文从含有脉冲噪声的频带和其他频带中分别选取10个频带加入宽度为1～3的死线.
本文用4个评价指标度量去噪效果：
平均峰值信噪比(mean peak signal-to-noise ratio, MPSNR)、平均结构相似度(mean structural similarity, MSSIM)、相对无量纲全局误差(relative dimensionless global error synthesis, ERGAS)和平均光谱角(mean spectral angle, MSA).表1列出了不同算法在上述3个数据集上的去噪效果, 包括在4种模拟噪声情况下不同评价指标的值.由表1可见, 本文算法的去噪结果在大多数情况下都最好.
表1 不同算法在3个数据集上的定量实验结果
续表1
图2为不同算法在数据集PaviaC频带64上的去噪结果.由图2可见, 本文算法更好地保留了图像的细节信息并去除了噪声.表2列出了不同去噪算法的运行时间.由表2可见, 本文算法在计算时间方面具有优势, 包括不同算法在3个数据集上的平均运行时间.
图2 不同算法对数据集PaviaC在加入第4种噪声后的去噪结果Fig.2 Denoising results of different algorithms on PaviaC dataset after adding the fourth noise
表2 不同去噪算法的运行时间
3.2 真实数据实验
由于真实数据集KSC没有干净的图像作为参考, 因此本文采用带有径向基函数(radial basis function, RBF)核的支持向量机得到的分类结果定量评估本文算法在真实数据集上的去噪性能.表3列出了不同算法每个类别的分类准确率, 并给出了3个评价指标, 分别为总体准确率(overall accuracy, OA)、平均准确率(average accuracy, AA)、Kappa统计量.由表3可见, 本文算法在3个评价指标上都取得了优于其他算法的分类结果, 表明本文算法在真实数据集上同样具有较好的去噪效果.图3为不同算法在数据集KSC频带89上的去噪效果.由图3可见, 本文算法去噪后的图像具有最好的视觉效果.
表3 基于不同去噪算法数据集KSC的分类结果
图3 不同算法在数据集KSC上的去噪结果 Table 3 Denoising results of different algorithms on KSC dataset
综上所述, 本文提出了一种基于子空间表示的加权低秩张量正则化的高光谱图像去噪方法, 针对高光谱图像中的高斯噪声和稀疏噪声进行建模.相比于其他去噪算法, 本文模型较好地保留了张量的固有结构, 提高了高光谱图像先验信息的利用效率, 并且计算复杂度低于多数去噪算法.实验结果表明, 本文方法在模拟噪声图像和真实高光谱图像上都取得了较好的实验结果.
猜你喜欢张量频带集上偶数阶张量core逆的性质和应用数学物理学报(2021年1期)2021-03-29四元数张量方程A*NX=B 的通解五邑大学学报（自然科学版）(2020年4期)2020-12-09Cookie-Cutter集上的Gibbs测度数学年刊A辑(中文版)(2020年2期)2020-07-25Wi-Fi网络中5G和2.4G是什么？有何区别？电脑知识与技术·经验技巧(2020年5期)2020-06-22链完备偏序集上广义向量均衡问题解映射的保序性数学物理学报(2019年6期)2020-01-13分形集上的Ostrowski型不等式和Ostrowski-Grüss型不等式井冈山大学学报(自然科学版)(2019年4期)2019-09-09单音及部分频带干扰下DSSS系统性能分析航天电子对抗(2019年4期)2019-06-02扩散张量成像MRI 在CO中毒后迟发脑病中的应用山西大同大学学报(自然科学版)(2016年2期)2016-12-12调谐放大器通频带的计算及应用电测与仪表(2015年12期)2015-04-09LTE－U通信技术(2015年8期)2015-02-12