基于快速三因子分解和组稀疏正则化的高光谱图像去噪.docx

    基于快速三因子分解和组稀疏正则化的高光谱图像去噪    高小雨，白静远，黄扬智，宁纪锋（1 西北农林科技大学 信息工程学院，杨凌 712100）（2 西北农林科技大学 理学院，杨凌 712100）0 引言使用光谱仪对同一场景获取不同光谱的图像被称为高光谱图像（Hyperspectral Image，HSI），HSI 能提供丰富的光谱信息，被各领域［1-3］广泛应用然而，由于光照条件、传输条件和成像仪器等多种客观因素的限制，数据在采集和传输过程中不可避免地受到各种噪声的污染，并且由于其丰富的光谱维，其受到的噪声干扰可能性相对更大，噪声会降低图像的辨识度，影响图像质量，并且限制了分类［4-5］和解混［6］等后续处理任务的精度因此，作为预处理步骤，HSI 去噪是一个重要的研究课题目前许多针对HSI 的去噪模型［7-8］都是基于HSI 中的两大重要先验信息，即光谱域的低秩性和空间域的稀疏性HSI 的光谱维度存在低秩先验，即HSI 提供了被观测对象同一场景下数十个至数百个连续波段的光谱信息，所以不同波段之间存在着高相关性，然而直接求解秩函数最小化是一个凸优化问题，CANDES E J等［9］提出的鲁棒主成分分析模型（Robust Principal Component Analysis，RPCA）利用非凸的核范数对秩函数进行凸逼近，取得了良好的性能。

ZHANG Hongyan 和HE Wei 等［10］在RPCA 的框架下对图像进行块处理，然后沿光谱维将三维立方体展开成矩阵，最后使用低秩矩阵复原（Low-Rank Matrix Recover，LRMR）模型恢复矩阵，该做法更好地保存了局部细节和纹理信息，取得了较好的去噪效果有研究［11-12］指出，RPCA 模型中的核范数虽然便于计算但是在逼近秩函数时不准确，因此，研究者们致力于开发更精确的低秩近似表示HUANG Xinjian 等［13］使用加权核范数来约束低秩先验，提升对秩函数逼近的性能ZENG Haijin 等［14］直接定义了张量的秩，提出γ范数来描述HSI 局部块的低秩性质这些基于低秩矩阵逼近（Low-Rank Matrix Approximation，LRMA）的方法，是通过引入非凸函数来逼近矩阵的秩，但是这会涉及到奇异值分解（Singular Value Decomposition，SVD），存在计算复杂度高的缺点为解决以上难题，LIU Yuanyuan 等［15］对低秩矩阵进行快速三因子分解（Fast Tri-Factorization，FTF），将矩阵分解为两个正交因子矩阵和一个秩为r的核心矩阵。

FTF 计算复杂度低，速度快，并且还将低秩约束转移到规模更小的核心矩阵上，所以探索矩阵的低秩性有更高的效率LIU Qiang 等［16］在FTF 框架下提出了一种基于QR 分解的L2，1范数最小化方法该方法在对核心矩阵的约束上，用L2，1范数最小化来替换核范数最小化，这两种约束具有相同的最优解，而且使用L2，1范数相较于核范数能够带来计算速度上的显著提升同时，HSI 空间维的稀疏性也是常用于去噪的重要先验信息对于HSI，不同波段的成像场景相同，不同波段的分段平滑结构也应该相同，即所有波段的差分图像在光谱维度上应服从组稀疏结构这启发研究者们考虑用组稀疏正则化［17-18］去表示不同波段的空间差分图像之间的内部结构特征CHEN Yong 等［19］用组稀疏正则化去除条纹，取得了比使用稀疏正则化更好的效果文献［20］利用加权L2，1范数最小化来约束空间差分图像，但是只对两个空间维度上的差分图像进行约束，没有考虑到光谱维度现有的去噪模型对于矩阵秩最小化问题，多采用核范数最小化进行不断迭代求解，而每次迭代都涉及到SVD，所以这些算法有很高的计算复杂度；除此之外，全变差（Total Variation，TV）项无法探索空间差分图像的共享组稀疏模式。

为了更快速地表示低秩，更精确地表示稀疏，同时保证去噪效果的提升，在FTF 与全局组稀疏的启发下，本文提出了基于快速三因子分解和组稀疏正则化的高光谱图像降噪（Fast Trifactorization And Group Sparsity Regularized Hyperspectral Image Denoising，FTFGS）模型，以更有效地表达HSI 的低秩性和全局组稀疏1 相关工作1.1 HSI 的退化模型假设HSI 受到不同种类噪声的污染，该退化模型可表示为式中，O是观测到的含噪图像；L是清晰的图像；S表示包括脉冲噪声、死线、条带等稀疏噪声；N是高斯噪声；M×N表示HSI 的空间大小，p表示波段数1.2 局部RPCA文献［21］指出局部空间内像素属于同一地物的可能性更大，具有更强的相似性，并利用RPCA 模型来探索这种局部低秩先验，提出局部RPCA 模型基于此，定义取块算子Pi，j：=O→Oi，j=Oi，j，表示从高光谱图像O的空间位置(i，j)截取大小为m×n×p的子块并逐波段列化成矩阵Oi，j相应的，Li，j和Si，j分别是清晰HSI 和稀疏噪声张量对应的子矩阵，λ是正则化参数，ε表示高斯噪声的噪声水平。

那么，局部RPCA 去噪模型可表示为文献［8，21］等指出局部RPCA 以分块方式进行研究具有诸多优势首先，在局部区域中探索低秩属性更符合HSI 的物理特性；其次，局部模型可以减少对噪声是独立同分布这一假设的依赖；最后，这种做法避免了直接将HSI 列化后形成病态矩阵，有利于保护局部图块中的细节信息1.3 快速三因子分解文献［15］为了降低核范数最小化问题的计算成本，提出了一种基于QR 分解的FTF 方法对于秩为r的矩阵M∈Rm×n，它的求解模型为式中，A∈Rm×r是列正交矩阵，C∈Rr×n是行正交矩阵，核心矩阵B∈Rr×r，r≪min(m，n)FTF 模型将矩阵M的低秩性转移到核心矩阵B的低秩性，即文献［16］对B的秩进一步约束，即不仅对r进行限制，而且使用L2，1范数最小化来替换核范数最小化，提出模型式中，A和C可通过QR 分解快速求解且计算复杂度仅是SVD 的10%［16］相比于传统的SVD，在FTF 框架下采用基于QR 的L2，1范数不但降低了计算的复杂度，而且利用对核心矩阵的秩和范数的约束可以更准确、快速地体现图像的局部低秩特性1.4 SSTV 的组稀疏表示由于HSI 波段间的高度相关性，不同波段的分段平滑结构也应相同，即，一波段的某一区域是边界，那么相邻波段也是边界；如果该区域平滑，相邻波段也是平滑的。

对HSI 进行梯度操作之后，边界都非零，平滑区域近似为0，这一物理特性就是HSI 的组稀疏特性文献［20］采用加权的L2，1范数探索了HSI 空间上两个维度梯度张量的组稀疏结构，其数学表示为式中，F=[Fx，Fy]，Fx和Fy是空间维度的两个差分算子，Fx L=L(x+1，y，z)-L(x，y，z)，Fy L=L(x，y+1，z)-L(x，y，z)；w=[wx，wy]中wx和wy表示对空间上两个维度添加的权重2 基于FTFGS 的HSI 去噪模型2.1 FTFGS 模型在图1 中，对HSI 的三个梯度张量每个像素处的梯度光谱曲线求L2范数，并画出它们的直方图从直方图分布来看绝大多数L2范数值集中在固定的值，表明HSI 的梯度张量具有组稀疏的特性，而且这种组稀疏特性不仅存在于空间梯度张量中，也存在于光谱方向的梯度张量中为了探索HSI 每个梯度方向上的组稀疏结构，本文提出一个全新的加权空谱组稀疏正则化‖w⊙FL‖2，1，即图1 探索HSI 的组稀疏性Fig.1 Exploring group sparsity in HSI式中，F=[F1，F2，F3]=[bxFx，byFy，bzFz]，约束因子b可以控制不同方向的梯度对模型的影响，Fx、Fy和Fz分别是空间维度和光谱维度的三个差分算子；w=[wx，wy，wz]则是对三个维度添加的权重。

相较于只对空间维进行约束，对HSI 空间和光谱维度的梯度张量都进行组稀疏正则化，会取得更好的效果本文将全局组稀疏和局部FTF 的低秩表示相结合提出FTFGS 模型，即式中，Ci，j、Di，j和Ri，j分别是将Li，j进行FTF 之后得到的列正交矩阵、核心矩阵和行正交矩阵，λ、γ和τ都是正则化参数，以平衡各个正则项达到更好的去噪效果从式（8）可以看出FTFGS 模型的提出侧重于为低秩和稀疏分量分别开发更精确的近似表示，并在保证最优解的同时带来速度的提升模型的整体架构如图2 所示可以清楚地看到，局部模块采用局部RPCA 模型探索HSI 低秩性，不仅符合HSI 的物理特性，避免形成病态矩阵，而且可以更好地保护局部图块中的细节信息在处理小规模矩阵时，不同于使用传统SVD 方法［10］，FTF 可以达到速度更快、复杂度更低的效果在对核心矩阵进行处理时，不仅对r进行限制，而且用L2，1范数最小化在得到最优解的同时也能够带来计算速度上的显著提升在全局上，提出一个全新的加权空谱组稀疏正则化，相较于只在空间维度进行约束的模型［20］，能充分探索HSI 三个维度的组稀疏特性从而取得更好的效果，并且可以有效抑制高斯噪声，消除分块处理带来的人为结构。

FTFGS 模型实现了局部低秩模块和全局组稀疏模块的有效联动，迭代求解模型至收敛条件，最终不仅可以去除各种混合噪声，而且减少了对噪声独立同分布假设的依赖，有效抑制与结构相关的噪声图2 模型架构Fig.2 Model architecture2.2 模型求解采用交替乘子法［22］（Alternating Direction Method of Multipliers，ADMM）求解目标函数，并引入辅助变量J，X∈RM×N×p，U∈RM×N×p×3，式（8）可改写成式中，F(⋅)是空谱TV 算子，U=[u1，u2，u3]=[F1X，F2X，F3X]利用增广拉格朗日乘子法，式（9）优化问题可以转化成式中，ΓOi，j，ΓDi，j，ΓLi，j，ΓX和Γ是拉格朗日乘子，μ是惩罚参数在第k次迭代中，问题的解可以转化为局部低秩分解子问题与基于SSTV 的组稀疏子问题，即2.2.1 固定其他参数，求解子优化问题(C,D,R,L,S,N)1）求解Ci，j，Ri，j和Di，j求解Ci，j，目标式（11）变为基于QR 分解求L2，1范数最小化方法［16］得到按列求解Ci，j得求解Ri，j，目标式（11）变为类似Ci，j求解步骤，可得按行求解Ri，j得求解Di，j，将式（11）改写成2）求解Li，j计算得到3）求解Si，j可利用软阈值收缩算子计算，即4）求解Ni，j计算可得2.2.2 固定其他参数，求解子优化问题(J,X,U)1）求解J凸优化模型得到解析解式中，I∈RM×N×P是元素全为1 的张量，Pi，j是取样算子，PTi，j是其逆算子。

2）求解X利用快速傅里叶求解X得式中，fftn(·)表示快速傅里叶变换，ifftn(·)是其逆变换3）求解U根据文献［23］给出的引理计算U更新拉格朗日乘子和惩罚参数为在每次迭代中组稀疏权重w(i，j)更新为综上所述，得到FTFGS 模型的求解步骤为：1）输入观测到的HSI 图像O，图块大小blocksize，步长stepsize，秩大小r，正则化参数λ，γ，τ，（bx，by，bz），迭代终止条件ε以及参数更新率ρ2）初始化变量L=J=X=0，U=0，ΓOi，j=ΓDi，j=ΓLi，j=ΓX=Γ=0，μ=03）进入重复迭代：基于式（14）、（19）、（17）、（21）、（23）和（25）更新(Ci，j，Di，j，Ri，j，Li，j，Si，j，Ni，j)；基于式。