1 引言

近年来,高光谱遥感技术^[1]迅猛发展,极大地提高了高光谱图像(HSI)的空间和光谱分辨率,被广泛应用于城市规划、勘探测绘、农作物监测和生物医学^[2-3]等领域。但HSI在获取或传输过程中,由于热电子、暗流以及在成像过程中光计数的随机误差,会受到各种类型噪声的污染,包括高斯、脉冲、条带和坏行列等^[4-5]。噪声不仅降低图像质量,而且会降低后续HSI处理与信息提取的准确性。

HSI降噪就是在合理消除噪声的同时,尽可能多地保留有效的图像信息。由于高光谱数据具有“图谱合一”^[6]的优势,从降噪处理过程中是否充分利用图谱信息来考虑,现有降噪方法可分为空间域、光谱域和空谱结合三种。空间域是先将每个波段的数据看作灰度图像,再扩展到HSI降噪,但会忽略HSI像元间特有的光谱信息,如块匹配和3D过滤(BM3D)^[7]。光谱域是考虑HSI波段间的光谱信息,但忽略了空间信息,如低秩矩阵恢复(LRMR)^[8]和组低秩表达(GLRR)^[9]。空谱结合是在光谱域降噪的基础上,同时在HSI子立方体上利用局部或全局空间信息对图像进行降噪,但关键是如何将空间与光谱结合,现有方法有利用低秩表达、小波、张量分解和稀疏表示等技术,如双线性低秩矩阵分解(BLRMF)^[10]、主成分分析与小波变换^[11]、联合空间光谱低秩惩罚(SSLR)^[12]和低阶张量恢复(LRTR)^[13]。

光谱的低秩特性是当前的一个研究热点,被应用于降噪和异常目标检测^[14]等领域,因利用HSI高相关性的光谱信息^[15]已被证明其具有更好的降噪结果^[16]。LRMR方法是第一次发现并运用HSI光谱域中潜在的低秩(LR)结构,利用低秩矩阵恢复模型对HSI混合噪声进行去噪。该方法将HSI根据空间信息按先行后列的顺序分成多个重叠的正方形块,依次在块内采用GoDec(Go Decomposition)算法进行降噪,但最优块的大小要通过实验确定。文献[ 17]针对LRMR未考虑整幅谱带噪声方差的变化,提出了噪声校正迭代低阶矩阵近似(NAILRMA)。为了更好地利用空间信息,He等^[18]在LRMR的基础上,考虑了相邻像素间的空间信息,提出全变分正则低秩矩阵分解(LRTV),整合低秩方法和全变分(TV)并应用于混合噪声降噪。以上方法都是针对局部块内的光谱低秩性进行去噪,并未考虑整幅图像的低秩性,同时引入TV导致计算复杂度变高。Fan等^[19]利用超像素将图像分割成不同均匀区域,并对其进行低秩去噪。但基于图像的熵率超像素分割^[20]选择K参数,会对实验结果产生影响,特别是对不同数据集产生的最优K值需反复实验来确定。

在HSI数据中,根据同类像素间的光谱反射率一致性更强的特点^[21],发现相同类别的地物间具有非常相似的光谱特征。如果在光谱域中不加以考虑,这种空间上的相关性极有可能会受到去噪算法的影响。基于低秩方法采用的局部块中全局谱带的差异性较大,在一定程度上破坏了HSI的全局低秩特性,而在同地物分块中可有效增强低秩特性。因此,本文借助HSI的先验知识,即地物类别信息,提出一种结合地物类别和低秩特性(GTLR)的HSI降噪方法。该方法结合地物类别简单而有效地确定分块的大小和个数,先利用块内形成矩阵的空间和光谱低秩特性,再考虑整幅图像的低秩特性,最后采用改进的IALM(Inexact Augmented Lagrange Multiplier)方法^[22]进行模型求解。在公共数据集Washington DC Mall和Indian Pines上分别进行模拟噪声实验和真实图像降噪实验,证明了所提降噪方法的有效性。

2 基于低秩矩阵恢复HSI模型

2.1 LRMR模型

LRMR模型先由Wright等^[23]提出并被理想化为“鲁棒主成分分析(RPCA)”问题。假设一个低秩矩阵L∈ℝ^m×n被一个稀疏误差矩阵S∈ℝ^m×n破坏,则观测数据矩阵D^*∈ℝ^m×n可分解为一个稀疏矩阵和一个低秩矩阵的和,即D^*=L+S。理想的RPCA问题:给定D^*,L和未知的S,目标是恢复L,该优化问题的表达式为

minL,S [rank(L)]+λ‖S‖0,s.t.D*=L+S,(1)

式中:‖ · ‖₀为0范数矩阵;rank为矩阵的秩;λ为平衡两者的正则化参数。(1)式为一个高度非凸优化问题,通过凸松弛将矩阵0范数替换为1范数,将秩替换为核范数得到可计算的优化问题,表达式为

minL,S‖L‖*+λ‖S‖1,s.t.D*=L+S,(2)

式中:‖ · ‖_*为核范数矩阵,即矩阵奇异值(SV)之和;‖ · ‖₁为1范数矩阵,即矩阵绝对值之和;λ为用于平衡核范数和1范数间的相对贡献的正则化参数。

文献[ 24]提出求解(2)式的等价优化问题,即

minL,SL‖*+λ‖S‖1,s.t.‖D*-L-S‖F≤δ,(3)minL,S‖D*-L-S‖F2,s.t.rank(L)≤r,card(S)≤z,(4)

式中:‖ · ‖_F为F范数矩阵;δ为随机噪声的标准差;r为L秩的上界;card为S的基数;z为S中误差水平。

2.2 基于LRMR模型HSI降噪

假设原始HSI的图像立方体为D∈ℝ^w×h×b,其中w、h代表图像空间的行、列维数,b代表图像的光谱维数。按照LRMR的低秩矩阵转换公式,利用全局光谱相关性,将D∈ℝ^w×h×b转换成D∈ℝ^o×b,其中o=w×h,每列由一个波段图像数据组成。其中选取一子数据立方体B,大小为m×n×b,m×n是子数据立方体中的空间信息。

假设HSI被随机和稀疏混合噪声污染,二维矩阵形式下的LRMR退化模型为

D=A+E,(5)

其中:D∈ℝ^o×b为原始的高光谱图像,A∈ℝ^o×b为具有低秩性的无噪声图像;E∈ℝ^o×b 为混合类型噪声,可分为具有稀疏特性的噪声和随机噪声。在等价模型(3)式和(4)式中修正LRMR数据模型,为D=A+F+N,其中F∈ℝ^o×b主要是具有稀疏特性的噪声;N∈ℝ^o×b为随机噪声。因此,可在HSI中采用(2)式及其等价模型。

求解(2)式及其等价模型的较优方法,分别有GoDec^[24]、SSGoDec(Semi-Soft GoDec)^[8]和ALM(Augmented Lagrange Multiplier)^[22],ALM分别包含EALM(Exact ALM)和IALM(Inexact ALM)^[22]。由于在相似的方法中能够快速收敛到最优解,故而采用效率较高的IALM算法进行求解。

3 结合地物类别和低秩特性HSI降噪方法

先通过局部块间的相关性和同物块间的相关性进行实验,分析得到同物分块的可行性,定义同物空谱低秩,再详细描述GTLR方法。

3.1 同物空谱低秩

基于低秩方法,如LRMR等,通过反复实验确定最优降噪结果的正方形局部块的大小和步长,根据超像素分割法确定分块的大小和个数,但目前关于分块大小的选择并没有固定方式。

实验中LRMR使用SSGoDec算法,得到最优解的局部块中各个大小均固定不变,如LRMR方法中采用最优块的大小为20×20×191。在结合地物类别实验中,为了区别局部光谱低秩,此时定义为同物空谱低秩。当采用的像素数小于191(波段数)时,形成矩阵行列间的相关性比LRMR最优块大小形成的矩阵更强,具有更强的低秩特性。当像素数大致为20×20=400时,分别对结合地物类别信息分块和最优局部分块的低秩性探讨如下。

图1中SV(Singular Value)曲线急剧衰减趋势表明,具有相同大小的局部块间的相关性和同物块间的相关性一致。因此,在相同大小的情况下,分析同物选择的立方体与局部立方体间的低秩属性,确实存在与局部块相似的相关性,明确分块大小,证明提出的同物分块法具有可行性。

图 1. HSI低秩特性。(a)局部块;(b)同物块

Fig. 1. HSI low rank feature. (a) Local block; (b) object block

下载图片 查看所有图片

文献[12]采用前k个SV累积能量比(AccEgyR)分析低秩属性,定义XAccEgyR=∑jkσj/∑jnσj,其中σj代表第j个SV,n=S,表示所有SV矩阵的总个数。图1(a) SV曲线中左侧箭头表示前10项SV可达到的比率,为0.9522,右侧箭头表示前56项可达到的比率,为0.9903;图1(b) SV曲线中左侧箭头表示前12项SV可达到的比率,为0.9520,右侧箭头表示前63项可达到的比率,为0.9902。这些定量值表明,局部块的相关性和同物块的相关性具有一致性,即同物空谱低秩与局部光谱低秩都具有较强的相关性。

3.2 GTLR降噪方法

假设图像中共有l种地物,表示为∑={i:i=1,2,…,l},其中i代表某种地物。基于地物类别的低秩矩阵恢复模型表示为

minA,E∑i=0l‖Ai‖*+λ‖Ei‖1),s.t.Di=Ai+Ei。(6)

针对(6)式的约束最优化问题,即

(A,E)=argminA,E∑i=0l‖Ai‖*+λ‖Ei‖1),s.t.Di=Ai+Ei。(7)

(7)式针对每个子立方体B⁽ⁱ⁾构造增广拉格朗日函数,即

L(Ai,Ei,Λ,μ)=‖Ai‖*+λ‖Ei‖1+<Λ,Di-Ai-Ei>+μ2‖Di-Ai-Ei‖F2,(8)

式中:μ为惩罚参数;Λ为拉格朗日乘子。若求解(8)式,先求解A⁽ⁱ⁾,再求解E⁽ⁱ⁾,即

Ak+1i=argminL(Ai,Eki,Λk,μk),(9)Ek+1i=argminL(Ak+1i,Ei,Λk,μk),(10)Λk+1=Λk+μk(Di-Ak+1i-Ek+1i),(11)

式中:k为迭代次数。(9)式和(10)式存在闭式解,可选择采用IALM算法进行求解,分别可得到 Ak+1i=US_μ-1k [S]V^T, Ek+1i=S_λμ-1k [Z],其中(U,S,V)=svd(W),W= Dki- Eki+ μk-1Λ_k,Z= Dki- Ak+1i+ μk-1Λ_k,U、V分别代表左、右SV矩阵,S代表SV对角矩阵,S_ε[x]代表软阈值(收缩)运算符,svd代表奇异值分解,ε代表相关参数。S_ε[x]表达式为

Sε[x]=x-ε,x>εx+ε,x<-ε0,x≤ε。(12)

实验先利用同物空谱低秩特性,根据HSI的地物类别信息,当选择类标号为i的所有像素组成一个B⁽ⁱ⁾时,HSI降噪的模型分为以下两种情况:当以图像块内的像素间的空间相关性为主时,即形成矩阵D⁽ⁱ⁾∈ℝ^p×b的行(row)相关性;当以图像块内的像素间的光谱相关性为主时,即形成矩阵 (Di)T∈ℝ^p×b的列(column)相关性,其中p= Bi代表类标号为i的按先行后列的所有像素数。由于同地物像素都可看作一个不重叠区域,形状确定,因此只需将基于不同地物的区域进行连接,即可获得降噪后的HSI。上述描述如图2(a)过程,其中地物类别信息与HSI数据一样来源于同一数据官方网站,类别标签的精度在实验中并未有特殊要求。

再利用HSI的全局光谱相关性,在图2(a)局部块低秩性的基础上,进一步利用整幅HSI的光谱低秩特性,达到再次降噪的目的,过程如图2(b)所示。

HSI数据的处理过程:实验前经归一化、降噪后各个波段值还原为原始数据范围。关于模型参数λ的取值,在某个B⁽ⁱ⁾∈ℝ^m×n×b中,当以像素间的空间相关性为主时,使用最优参数λ=1/[max(mn, b)]^1/2;当以矩阵中光谱相关性为主时,使用函数中默认的参数λ=1/mn,其中mn代表某个图像块内的像素总数。

综上所述,GTLR算法的主要步骤如下。

输入D∈ℝ^w×h×b,输出A^w×h×b,初始化λ=1/[max(mn, b)]^1/2或λ=1/mn。

1) 借助标签信息对HSI中ℝ^w×h×b进行分块,生成分块B={B⁽¹⁾,B⁽²⁾,…,B^(L)}。

2) 在某个分块B⁽ⁱ⁾∈ℝ^m×n×b中,构造矩阵D⁽ⁱ⁾∈ℝ^b×p。

3) 采用IALM算法,从原始矩阵D⁽ⁱ⁾中依次恢复信号矩阵 A1i。

① 当pw,h)×2时,使用{ Ai1, Ei1}=X_IALM(D⁽ⁱ⁾)。

② 当p> max(w,h)时,使用{ Ai1, Ei1}=X_IALM[ (Di)T]。

③ 以空间低秩为主,采用步骤①;以光谱低秩为主,采用步骤②。

4) 根据标签原始位置,还原低秩部分A⊂ℝ^w×h×b。

5) 将低秩部分A^w×h×b转换为D∈A^o×b,再次恢复信号矩阵 Ao×b2。

6) 得到降噪后的A^w×h×b数据。

图 2. GTLR降噪过程。(a)利用同物块空谱低秩;(b)利用全局图像光谱低秩

Fig. 2. GTLR noise reduction process. (a) Using spatial-spectral low rank of object block; (b) using global image spectral low rank

下载图片 查看所有图片

4 实验结果与分析

为了验证所提方法的降噪性能,在不同实验数据集上分别进行模拟实验和真实图像降噪实验,比较GTLR方法与LRMR方法和NAILRMA方法的实验结果,以证明其优越性。

对模拟实验结果采用图像质量评价指标进行评价,分别计算模拟添加HSI噪声数据与原始数据(第一种情况)、降噪后的HSI数据与原始数据(第二种情况)两种情况下的指标值。指标包含信噪比(SNR)、峰值信噪比(PSNR)、结构相似性(SSIM)和特征相似度(FSIM)。计算第一种情况,只计算SNR和PSNR两个指标,分别代表整幅图像的输入信噪比和输入峰值信噪比。计算第二种情况,当SNR和PSNR值较高时,降噪后的图像在视觉上更接近原始图像。SSIM和FSIM值越接近1,表明降噪后的图像与原始图像越相似,降噪后的结果对图像结构和特征保留更好。评价指标是计算HSI中所有光谱波段的图像,再计算这些波段的平均值,分别表示为MPSNR、MFSIM和MSSIM^[25]。

4.1 实验数据集

实验中共采用两种数据集,其中Washington DC Mall数据集用于模拟实验,实际受噪声污染的Indian Pines数据集用于真实图像的降噪实验。

Washington DC Mall数据集^[26]由图像尺寸为1208 pixel×307 pixel和191个波段组成,数据成像质量较高,故而选择为无噪的参考图像。图1中左侧的三维立方体图是由60、27和17波段的红色、绿色和蓝色制作的假彩色图像,将其作为模拟实验选取的数据集。

Indian Pines数据集^[27]由图像尺寸为145 pixel×145 pixel和224个波段组成,由于存在20条覆盖吸水区域严重的波段,去除波段104~108,150~163和220后,采用最后保留的200个波段。

4.2 模拟降噪实验

4.2.1 模拟噪声

主要针对大部分实际HSI中存在的随机噪声和稀疏噪声的混合情况,按照惯有模拟添加噪声的方式,同时添加以高斯噪声为代表的随机噪声和以椒盐噪声为代表的稀疏噪声,主要考虑上述两种类型噪声的混合噪声。

噪声情况1:选择实验数据集时,对每个波段添加随机高斯方差范围为0~0.1的噪声,随机选取某几个波段添加脉冲噪声。目前在波段20~25和波段140~142共选定9个波段添加脉冲噪声,方差为0.15。最终输入SNR值为6.4483 dB,MPSNR值为14.2516 dB。

噪声情况2:添加高斯方差范围为0~0.1的噪声,HSI光谱波段中随机选取1~2、20~21、50~51、73~74、120~121、142~143和188~189波段添加均值为0(其中波段对应的随机方差与噪声情况1相同),脉冲噪声与噪声情况1一致,最终输入SNR值为15.6367 dB,MPSNR值为1.3644 dB。

噪声情况3:每个波段添加0~0.2随机方差的高斯噪声,随机选择20~30、70~71、73、140~142共17个波段添加密度为0.2的脉冲噪声。输入SNR值为12.7994 dB。单独脉冲噪声的平均SNR值为15.3862 dB。单独高斯噪声的平均SNR值为3.3992 dB。最终输入SNR值为12.7994 dB,MPSNR值为0.7990 dB。

4.2.2 降噪结果

根据对比方法的参数选取,采用最优参数来获得最佳降噪结果。在三种模拟噪声情况下,按照LRMR和NAILRMA对比方法和所提方法GTLR(包含Row GTLR和Column GTLR,分别简写为R-GTLR和C-GTLR)的降噪结果分别列出模拟噪声实验的降噪评价指标值,结果如表1~3所示。对于噪声情况1,所提方法的提升效果不明显,但指标MSSIM和MFSIM值都有所提高,分别最高提高了0.0409和0.0150(所提方法中效果最好的方法与对比方法中效果最好的方法的差)。但在表2和表3中(分别代表模拟噪声情况2和模拟噪声情况3),所提方法的提高比较明显,降噪效果更好。特别在表3中对于混合噪声污染更加严重时,C-GTLR方法的MPSNR值比LRMR和NAILRMA提高4 dB以上,指标MSSIM和MFSIM分别最高提高0.0308和0.0201。

根据噪声情况3对每波段的PSNR、SSIM和FSIM的降噪结果进行对比,如图3~图5所示,所提方法的MPSNR值更高,评价指标MSSIM和MFSIM值越接近于1,则降噪后的图像与原始图像越相似,表明结构和特征更好,降噪效果更优。所提方法不仅可去除高斯噪声,而且对于稀疏噪声去除的效果更好。综上所述,所提方法的噪声去除结果更优。

图 3. 降噪结果的PSNR值对比

Fig. 3. Comparison of PSNR values of noise reduction results

下载图片 查看所有图片

图 4. 降噪结果的SSIM值对比

Fig. 4. Comparison of SSIM values of noise reduction results

下载图片 查看所有图片

图 5. 降噪结果的FSIM值对比

Fig. 5. Comparison of FSIM values of noise reduction results

下载图片 查看所有图片

4.3 真实图像降噪

利用Indian Pines真实存在噪声的图像数据,借助于地物标签信息,共有17种地物类别,其中像素总数小于500有9种,像素总数位于500~5000之间有7种,像素总数大于5000有1种。先对原始HSI数据进行归一化操作,再分别采用LRMR、NAILRMA和模拟实验中最优的C-GTLR方法进行降噪,最后还原结果数据并输出灰度图。其中三个方法的用时分别为24.2249 s、89.2688 s、18.0848 s,所提方法用时最短。

对原始Indian Pines的每个波段视觉进行观察,发现存在混合噪声的波段有1~3、61、75~76、103~107、144~146和198~200,累计共17个波段。在这17个波段中,所提方法的降噪视觉效果更好。由于篇幅限制,故而选择具有代表性波段1~2和103经降噪还原后的输出对比结果,如图6~8所示。从图6~8可以看到,与对比方法相比,所提方法在混合噪声较弱部分,降噪效果一般,但在随机和稀疏噪声较为严重部分,降噪性能更优。

图 6. Indian Pines图像在波段1的降噪结果。(a)原始图;(b) LRMR;(c) NAILRMA;(d)所提方法

Fig. 6. Band 1 noise reduction results of Indian Pines image. (a) Original image; (b) LRMR; (c) NAILRMA; (d) proposed method

下载图片 查看所有图片

图 7. Indian Pines图像在波段2的降噪结果。(a)原始图;(b) LRMR;(c) NAILRMA;(d)所提方法

Fig. 7. Band 2 noise reduction results of Indian Pines image. (a) Original image; (b) LRMR; (c) NAILRMA; (d) proposed method

下载图片 查看所有图片

图 8. Indian Pines图像波段103降噪结果。(a)原始图;(b) LRMR;(c) NAILRMA;(d)所提方法

Fig. 8. Band 103 noise reduction results of Indian Pines image. (a) Original image; (b) LRMR; (c) NAILRMA; (d) proposed method

下载图片 查看所有图片

4.4 参数讨论

求解过程中存在两大问题:一是同物块过大或过小,二是惩罚参数值λ的选取。根据所提方法的思路,选取分块大小最优的结果与结合地物类别标签信息分块的结果进行对比,发现结果相差不大,故而证明结合标签信息分块的正确性。采用结合地物类别信息进行分块时,需考虑在某个相同地物类别中,像素数太大和太小时如何划分图像块。如模拟实验数据集,大小为256×256×191,实验过程中降噪方式以空间为主,当像素数大于256×2=512时,则不采用R-GTLR方法降噪方式;实验过程中降噪方式以光谱为主,当像素数小于行或列值的绝对值时,即小于sqrt(256)=16,则不采用C-GTLR方法。根据实验过程总结,当相同地物下像素数过大时,如Indian Pines数据中标签号为0的总像素值有10776个,则采用以空间为主的GTLR方法,可在该地物类别内重新按照指定默认局部块大小依次选取某个分块B⁽ⁱ⁾∈ℝ^m×n×b中max(m,n)的像素大小;以光谱为主的GTLR方法可无需依次选取分块。

5 结论

基于低秩降噪方法,考虑现有方法存在的不足,针对高光谱遥感数据的特性,分析数据中局部分块形成矩阵间的相关性和同地物分块形成矩阵间的相关性,发现当同物分块像素数与HSI波段数大致相等时,两者具有一致性;当同物分块像素数小于HSI波段数时,具有更强的低秩特性。在利用地物先验信息分块的基础上,基于地物类别建立低秩矩阵恢复模型,利用同物块内低秩性去除混合噪声;再在去噪过程中引入HSI整幅全局光谱低秩性,进一步去除混合噪声。所提方法不仅能够明显去除稀疏噪声,同样可有效去除高密度随机噪声。同时,所提模型在类似均匀区域的同物块内降噪,耗费时间更短。