高光谱影像独有的高光谱分辨率特点为地理信息获取和地物识别带来了便利和优势^[1]，然而高光谱传感器在获取大量数据的同时易受外界干扰，导致数据污染，降低了数据精度。因此，减少数据受到的噪声污染，得到更多有效的图像信息，增加图像数据的利用率，是保证高光谱图像能够在后续应用中发挥作用的重要环节。

高光谱成像技术将传统的2维图像遥感技术和光谱技术相结合，可同时获得地表图像及其光谱信息。因此，高光谱影像同时包含空间域和光谱域的信息。且高光谱图像在空间域和光谱域的连续性使其在相邻的通道间的相关性很强，即存在着低秩的特性。早期基于低秩先验的主成分分析(PCA)^[2]被引入来实现高光谱信号和噪声的分离。该方法将3维高光谱影像平铺为2维图像，并采用主成分分析变换进行降维，权重大的主成分波段包含了高光谱影像的主要信息，而噪声平均分布在每个成分中。但是主成分分析变换主要针对的是高斯噪声，并且对异常值非常敏感，然而高光谱中包含着诸如条带噪声、死线噪声等大量异常值。此外，这种方法虽然充分利用了光谱信息，但空间信息丢失严重。全变差(TV)模型的提出充分保证了图像的边缘信息，能够用实现空间分段光滑，使得高光谱的空间信息得到更加充分的利用。1992年Rudin等人^[3]根据偏微分方程的思想首先提出了TV模型，各向同性的TV模型对高斯随机噪声的监测和去除有较好的效果，之后在1996年Li等人^[4]提出的各向异性的TV模型则对图像中存在的椒盐、脉冲等噪声有着较好的处理效果。TV模型简洁有效，但是带来的阶梯效应极大影响了图像的质量。2000年以后，出现了各类改进的TV模型，例如考虑图像像素间的相关性会被噪声破坏使得局部算法的效率降低而提出的非局部TV模型^[5]、为了弥补TV模型存在的分块效应而提出的高阶TV模型^[6]等也逐渐引入到高光谱噪声分析中。近几年，结合低秩先验模型与空间先验模型的优势，所建立的光谱低秩全变差正则化模型^[7-9]也取得了良好的恢复结果。

在低秩表示中，核范数即矩阵的特征值之和是矩阵秩函数的最紧密的凸下界，核范数和矩阵秩函数之间的关系与向量l₁范数和l₀范数之间的关系相似。因此，在众多研究中通常使用核范数作为秩函数的凸松弛，如核范数正则化最小二乘法(NNLS)^[10]、鲁棒的主成分分析法(RPCA)^[11]、空谱全变差正则化方法(LRTV)^[7]等。但如今基于核范数的恢复方法，虽然具有较强的理论依据，但在实际应用中得到的是次优的结果。在秩函数中，特征值包含着观测数据的信息，不同的特征值包含的信息不同，其中较大的特征值主要包含原始数据信息，而较小的特征值主要包含观测数据的噪声信息。但是基于核范数最小化的复原方法将所有的特征值同等对待，对应的模型算法利用相同阈值对特征值进行收缩处理，在抑制噪声的前提下丢失了大量图像信息，是该方法的一个重要缺陷。除此之外，核范数的一些理论需求(如不连贯性)在实际中很难令人满意。为了解决上述核范数方法存在的缺陷，Hu等人^[12]最近提出使用截断核范数, 利用观测数据较小特征值之和，并且充分利用上述特征值的性质，克服了核范数的缺点, 在此基础上，建立了基于最小化截断核范数改进矩阵填充模型。Cao等人^[13]应用截断核范数在RPCA模型上。基于截断核范数最小化问题属于非凸优化问题，之前相关的研究者^{[12, 14]}设计通过两阶段算法求解相应模型。作为该模型的进一步推广，Tao等人^[15-17]运用DCA(difference of convex functions algorithm)理论解决截断核范数最小化的问题。上述研究都丰富了截断核范数的理论，但是现存的研究还存在以下问题：1)截断核范数在高光谱去噪中的应用只停留在利用高光谱光谱低秩先验信息对其进行去噪的阶段，没有利用高光谱的空间低秩先验信息；2)现有的方法恢复结果并不能使人满意。

为了将截断核范数的内容充分应用于高光谱图像复原，提高复原精度，本文首先分析了现今主流的高光谱去噪方法，再对截断核范数的模型进行分析。最终提出了基于截断核范数下的低秩先验约束和TV正则化的高光谱遥感图像复原模型，不仅保留了当今主流模型的处理优点，同时充分利用截断核范数的性质，提高了图像的复原质量。

由于高光谱在成像过程中会受到噪声污染，使图像模糊，因此需要对高光谱影像进行噪声去除处理，使图像的细节纹理更加丰富。假设高光谱影像受到噪声的干扰，则其退化模型为

$ \mathit{\boldsymbol{Y}} = \mathit{\boldsymbol{X}} + \mathit{\boldsymbol{N}} + \mathit{\boldsymbol{S}} $

(1)

式中，$ \mathit{\boldsymbol{Y}} \in {{\bf{R}}^{m \times n \times b}}$表示观测噪声图像；$ \mathit{\boldsymbol{X}} \in {{\bf{R}}^{m \times n \times b}}$为无噪图像；$ \mathit{\boldsymbol{N}} \in {{\bf{R}}^{m \times n \times b}}$表示高密度噪声，例如高斯噪声、泊松噪声等；$ \mathit{\boldsymbol{S}} \in {{\bf{R}}^{m \times n \times b}}$表示具有稀疏先验的噪声，例如条带噪声、死线噪声、脉冲噪声等。其中$ m, n$表示数据的空间尺寸，$b $表示波段数。为了对图像进行处理，将每一波段的图像按列排成一列向量，即$ \mathit{\boldsymbol{Y}} \in {{\bf{R}}^{mn \times b}}$，则$\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{N}}, \mathit{\boldsymbol{S}} \in {{\bf{R}}^{mn \times b}} $。

由于高光谱图像不同波段间存在着高度相关性，即高光谱图像的光谱低秩先验。若将高光谱图像逐波段列化为二维矩阵，那么光谱低秩先验就可以用$\mathit{\boldsymbol{X}} $的核范数最小化近似表示。在最大后验估计理论的基础上建立秩约束的RPCA模型为

$ \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{S}}} {{\left\| \mathit{\boldsymbol{X}} \right\|}_*} + \lambda {{\left\| \mathit{\boldsymbol{S}} \right\|}_1}}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{. }}\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{S}}} \right\|_F^2 \le \varepsilon , \mathit{rank}\left( \mathit{\boldsymbol{X}} \right) \le r} \end{array} $

(2)

式中，$\parallel \cdot \parallel $表示核范数；$ rank\left( \cdot \right)$表示矩阵秩；$\lambda $是正则化参数；$ \varepsilon $表示噪声方差水平。

低秩矩阵分解模型只包含光谱约束而缺少空间约束，所以不能有效去除高密度噪声。因此Candès等人^[18]提出基于TV的ROF模型，它可以显著地保护边缘信息和分段光滑结构，然而光谱间的相似性却被忽略。基于以上的研究现状，He等人^{[7, 9]}将几种模型结合提出LRTV模型

$ \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{S}} \in {{\bf{R}}^{m \times n}}} {{\left\| \mathit{\boldsymbol{X}} \right\|}_*} + \lambda {{\left\| \mathit{\boldsymbol{S}} \right\|}_1} + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{{\rm{HTV}}}}}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{. }}\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{S}}} \right\|_F^2 \le \varepsilon , \mathit{rank}\left( \mathit{\boldsymbol{X}} \right) \le r} \end{array} $

(3)

式中，$ \tau $用来平衡核范数和TV范数，$\lambda $用来限制稀疏噪声的稀疏性。在LRTV模型中，当$ \tau $设置为0时，LRTV模型退化为低秩矩阵分解模型。

在LRTV模型中，利用核范数最小化理论与3D全变差正则化相结合的方法对秩函数进行了逼近。但是正如上文所述，核范数对于秩函数的逼近存在一系列的问题。因此，截断核范数的提出能更好地利用到观测数据的先验特征。该方法是利用当矩阵具有低秩结构时，噪声包含在相对较小的特征值中，利用较小的特征值之和来控制噪声的误差，从而更好地表示出数据的特征。但是截断核范数在高光谱图像恢复问题上的应用并不十分广泛。为此将在此节给出截断核范数的相关内容并将其应用在高光谱图像复原上。

给定一个秩为$t $的矩阵$\mathit{\boldsymbol{X}} \in {{\bf{R}}^{m \times n}} $，其中的噪声信息绝大程度上取决于最小的$ \min (m, n)-t$个特征值的非零元个数。它的截断核范数定义为

$ \begin{array}{*{20}{c}} {{{\left\| \mathit{\boldsymbol{X}} \right\|}_{t, * }} = {{\left\| {\sigma \left( \mathit{\boldsymbol{X}} \right)} \right\|}_{t, 1}} = \sum\limits_{i = t + 1}^{\min \left( {m, n} \right)} {{\sigma _i}\left( \mathit{\boldsymbol{X}} \right)} = }\\ {\sum\limits_{i = t + 1}^{\min \left( {m, n} \right)} {{\sigma _i}\left( \mathit{\boldsymbol{X}} \right)} - \sum\limits_{i = 1}^t {{\sigma _i}\left( \mathit{\boldsymbol{X}} \right)} } \end{array} $

式中，${\sigma _i}\left( \cdot \right) $表示矩阵第$i $个最大特征值。

记$\parallel \mathit{\boldsymbol{X}}{\parallel _{{t^c}, *}} = \sum\limits_{i = 1}^t {{\sigma _i}\left( \mathit{\boldsymbol{X}} \right)} $，则

$ {\left\| \mathit{\boldsymbol{X}} \right\|_*} = {\left\| \mathit{\boldsymbol{X}} \right\|_{t, * }} + {\left\| \mathit{\boldsymbol{X}} \right\|_{{t^c}, * }} $

(4)

根据2.1节的结果分析表明，截断核范数是比核范数更准确地表征观测数据的先验特征的度量^[12]，因此本文依据高光谱遥感数据光谱和空间低秩先验信息，提出新模型，即TNN-LRTV模型。

$ \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{L}}, \mathit{\boldsymbol{S}}} {{\left\| \mathit{\boldsymbol{L}} \right\|}_{t, * }} + \lambda {{\left\| \mathit{\boldsymbol{S}} \right\|}_1} + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{{\rm{HTV}}}}}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{. }}\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{S}}} \right\|_F^2 \le \varepsilon , \mathit{rank}(\mathit{\boldsymbol{L}}) \le \mathit{\boldsymbol{r}}} \end{array} $

(5)

式中，$ {\sigma _i}\left( \cdot \right)$表示矩阵的第$i $个最大特征值。

该模型在处理噪声时能同时兼顾稀疏噪声和强高斯噪声。如果参数$ \tau $取值为0，那么TNN-LRTV模型就退化为基于TNN的矩阵低秩和稀疏分解的正则化优化模型(LRSD-TNN)^[14], 如果用核范数代替模型中的截断核范数，那么TNN-LRTV模型就退化为LRTV模型^[7]。如果$ \tau $与$t $取值同时为0，那么模型就退化为秩约束下的RPCA模型^[11]。

ADMM(alternating direction method of multipliers)是一种适用于分布式凸优化的简单而有效的方法，该方法能分解原函数和扩增函数，以便于在更一般的假设条件下并行优化。因此本文选用ADMM方法^[19]对模型(5)进行求解。首先利用中间变量$ \mathit{\boldsymbol{L}}$同时添加$\mathit{\boldsymbol{L}} = \mathit{\boldsymbol{X}} $的约束项，那么可将模型(5)改写为

$ \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{L}}, \mathit{\boldsymbol{S}}} {{\left\| \mathit{\boldsymbol{L}} \right\|}_{t, * }} + \lambda {{\left\| \mathit{\boldsymbol{S}} \right\|}_1} + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{{\rm{HTV}}}}}\\ {{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{S}}} \right\|_F^2 \le \varepsilon , \mathit{rank}\left( \mathit{\boldsymbol{L}} \right) \le r, \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{X}}} \end{array} $

(6)

进一步可将模型(6)变为如下的无约束优化问题。

$ \begin{array}{*{20}{c}} {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{L}}, \mathit{\boldsymbol{S}}, \Lambda , rank(\mathit{\boldsymbol{L}}) \le r} E\left( {\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{L}}, \mathit{\boldsymbol{S}}, \Lambda } \right) = }\\ {\mathop {\arg \min }\limits_{\mathit{\boldsymbol{X}}, \mathit{\boldsymbol{L}}, \mathit{\boldsymbol{S}}, \Lambda , rank(\mathit{\boldsymbol{L}}) \le r} {{\left\| \mathit{\boldsymbol{L}} \right\|}_{t, * }} + \lambda {{\left\| \mathit{\boldsymbol{S}} \right\|}_1} + \tau {{\left\| \mathit{\boldsymbol{X}} \right\|}_{{\rm{HTV}}}} + }\\ {\left\langle {\Lambda , \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{X}}} \right\rangle + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{Y}} - \mathit{\boldsymbol{X}} - \mathit{\boldsymbol{S}}} \right\|_F^2 + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{X}}} \right\|_F^2} \end{array} $

(7)

式中，$\mu $是惩罚参数；$\Lambda $是拉格朗日乘子。使用交替优化方法求解模型(7), 在第$k + 1 $次迭代中，各个变量的更新形式为

$ {\mathit{\boldsymbol{L}}^{k + 1}} = \mathop {\arg \min }\limits_{rank\left( \mathit{\boldsymbol{L}} \right) \le r} E\left( {{\mathit{\boldsymbol{X}}^k}, \mathit{\boldsymbol{L}}, {\mathit{\boldsymbol{S}}^k}, {\Lambda ^k}} \right) $

(8)

$ {\mathit{\boldsymbol{X}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} E\left( {\mathit{\boldsymbol{X}}, {\mathit{\boldsymbol{L}}^{k + 1}}, {\mathit{\boldsymbol{S}}^k}, {\Lambda ^k}} \right) $

(9)

$ {\mathit{\boldsymbol{S}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{S}} E\left( {{\mathit{\boldsymbol{X}}^{k + 1}}, {\mathit{\boldsymbol{L}}^{k + 1}}, \mathit{\boldsymbol{S}}, {\Lambda ^k}} \right) $

(10)

$ {\Lambda ^{k + 1}} = {\Lambda ^k} + \mu \left( {{\mathit{\boldsymbol{L}}^{k + 1}} - {\mathit{\boldsymbol{X}}^{k + 1}}} \right) $

(11)

该复原模型(6)最终转化为3个主要的子优化问题(8)、(9)、(10)来求解。从而只需考虑3个子问题的求解。

针对子优化问题(8)，有

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}^{k + 1}} = \mathop {\arg \min }\limits_{rank\left( \mathit{\boldsymbol{L}} \right) \le r} E\left( {{\mathit{\boldsymbol{X}}^k}, \mathit{\boldsymbol{L}}, {\mathit{\boldsymbol{S}}^k}, {\Lambda ^k}} \right) = }\\ {\mathop {\arg \min }\limits_{rank\left( \mathit{\boldsymbol{L}} \right) \le r} {{\left\| \mathit{\boldsymbol{L}} \right\|}_{t, *}} + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{L}} - \left( {{\mathit{\boldsymbol{X}}^k} - {\Lambda ^k}/\mu } \right)} \right\|_F^2} \end{array} $

由式(4)可知, $\parallel \cdot {\parallel _{t, *}} $可以转换成两个凸函数相减, 即

$ {\left\| \mathit{\boldsymbol{X}} \right\|_{t, * }} = {\left\| \mathit{\boldsymbol{X}} \right\|_ * } - {\left\| \mathit{\boldsymbol{X}} \right\|_{{t^C}, * }} $

结合DCA方法即可得到两个凸函数

$ \left\{ \begin{array}{l} H\left( \mathit{\boldsymbol{L}} \right) = \lambda {\left\| \mathit{\boldsymbol{L}} \right\|_{{t^C}, * }}\\ G\left( \mathit{\boldsymbol{L}} \right) = \lambda {\left\| \mathit{\boldsymbol{L}} \right\|_ * } + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{L}} - \left( {{\mathit{\boldsymbol{X}}^k} - {\Lambda ^k}/\mu } \right)} \right\|_F^2 \end{array} \right. $

结合DCA理论, 得到迭代产生两个序列

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Z}}_j} \in \partial H\left( {{\mathit{\boldsymbol{L}}_j}} \right)}\\ {{\mathit{\boldsymbol{L}}_{j + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} G\left( \mathit{\boldsymbol{L}} \right) - \left\langle {{\mathit{\boldsymbol{Z}}_j}, \mathit{\boldsymbol{L}}} \right\rangle } \end{array}} \right. $

(12)

1) 对${\mathit{\boldsymbol{L}}_j} $进行SVD分解

$ {\mathit{\boldsymbol{L}}_j} = \mathit{\boldsymbol{U}}\left[ {\begin{array}{*{20}{c}} {{\rm{diag}}\left( d \right)}&0 \end{array}} \right]{\mathit{\boldsymbol{V}}^{\rm{T}}}, 取\;{d_i} = \left\{ {\begin{array}{*{20}{c}} 1&{i \le t}\\ 0&{i > t} \end{array}} \right. $

则${\mathit{\boldsymbol{Z}}_j} = \lambda \mathit{\boldsymbol{U}}[{\mathop{\rm diag}\nolimits} (d)\quad 0]{\mathit{\boldsymbol{V}}^{\rm{T}}} $

2) 经过化简后求解迭代(8)中的第2个子问题

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{L}}_{j + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{L}} \right) - \left\langle {{\mathit{\boldsymbol{Z}}_j}, L} \right\rangle = }\\ {\mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} \lambda {{\left\| \mathit{\boldsymbol{L}} \right\|}_ * } + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{L}} - \left( {{\mathit{\boldsymbol{X}}^k} - {\Lambda ^k}/\mu } \right)} \right\|_F^2 - \left\langle {{\mathit{\boldsymbol{Z}}_j}, L} \right\rangle = }\\ {\mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} \lambda {{\left\| \mathit{\boldsymbol{L}} \right\|}_ * } + \frac{\mu }{2}\left\| {\mathit{\boldsymbol{L}} - \left( {\mu {\mathit{\boldsymbol{L}}^k} - {\Lambda ^k} + {\mathit{\boldsymbol{Z}}_j}} \right)/\mu } \right\|_F^2 = }\\ {{D_{\frac{\lambda }{\mu }}}\left( {{\mathit{\boldsymbol{L}}^k} + \left( {{\mathit{\boldsymbol{Z}}_j} - {\Lambda ^k}} \right)/\mu } \right)} \end{array} $

则子优化问题(8)的解为

$ \left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{Z}}_j} \in \partial {{\left\| \mathit{\boldsymbol{X}} \right\|}_{{t^C}, * }}}\\ {{\mathit{\boldsymbol{L}}_{j + 1}} = {D_{\frac{\lambda }{\mu }}}\left( {{\mathit{\boldsymbol{X}}^k} + \left( {{\mathit{\boldsymbol{Z}}_j} - {\Lambda ^k}} \right)/\mu } \right)} \end{array}} \right. $

(13)

令$ \boldsymbol{Q}=\frac{\boldsymbol{S}^{k}-\boldsymbol{Y}-\boldsymbol{L}^{k+1}}{2}-\frac{\Lambda^{k}}{2 \mu}$，问题(9)可表示为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} E\left( {\mathit{\boldsymbol{X}}, {\mathit{\boldsymbol{L}}^{k + 1}}, {\mathit{\boldsymbol{S}}^k}, {\Lambda ^k}} \right) = }\\ {\mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} \tau \sum\limits_{i = 1}^b {{{\left\| {{\mathit{\boldsymbol{X}}_i}} \right\|}_{{\rm{TV}}}}} + \mu \left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Q}}} \right\|_F^2} \end{array} $

若将第$i $个光谱通道记做${Q_i} $，那么有

$ \mathit{\boldsymbol{X}}_i^{k + 1} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{X}} \tau {\left\| {{\mathit{\boldsymbol{X}}_i}} \right\|_{{\rm{TV}}}} + \mu \left\| {{\mathit{\boldsymbol{X}}_i} - {\mathit{\boldsymbol{Q}}_i}} \right\|_2^2 $

(14)

在本文中，利用基于梯度的快速算法(FGP算法)^{[11, 20]}对式(14)进行求解。

针对子优化问题(10)，有

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}^{k + 1}} = \mathop {\arg \min }\limits_\mathit{\boldsymbol{S}} E\left( {{\mathit{\boldsymbol{X}}^{k + 1}}, {\mathit{\boldsymbol{L}}^{k + 1}}, \mathit{\boldsymbol{S}}, {\Lambda ^k}} \right) = }\\ {\mathop {\arg \min }\limits_\mathit{\boldsymbol{S}} \frac{\lambda }{\mu }{{\left\| \mathit{\boldsymbol{S}} \right\|}_1} + \frac{1}{2}\left\| {\mathit{\boldsymbol{S}} - \left( {\mathit{\boldsymbol{Y}} - {\mathit{\boldsymbol{X}}^{k + 1}}} \right)} \right\|_F^2} \end{array} $

利用软阈值的方法直接给出它的解

$ {\mathit{\boldsymbol{S}}^{k + 1}} = \mathit{shrink}\left( {\mathit{\boldsymbol{Y}} - {\mathit{\boldsymbol{X}}^{k + 1}}, \frac{\lambda }{\mu }} \right) $

(15)

式中，$shrink(x, y) = {\mathop{\rm sgn}} (x) \cdot \max (|x| - y, 0) $。

综上所述，模型(6)的求解可以通过迭代运行式(13)、(14)、(15)和式(11)完成。

具体迭代过程参照算法1进行。

算法1 TNN-LRTV

输入: $\mathit{\boldsymbol{M}} \times \mathit{\boldsymbol{N}} \times \mathit{\boldsymbol{P}} $矩阵$\mathit{\boldsymbol{Y}} $

秩$r $

迭代停止参数$\varepsilon_{1}, \varepsilon_{2} $

正则化参数$ \tau $和$\lambda $

输出:复原图像$ \mathit{\boldsymbol{X}}$

初始化: $ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{X}} = \mathit{\boldsymbol{S}} = 0, {\Lambda _1} = {\Lambda _2} = 0, \mu = {10^{ - 2}}$, $ \mu_{\max }=10^{6}, \rho=1.5$以及$ k = 0$

重复迭代直到收敛:

更新$\boldsymbol{L}^{k+1}, \boldsymbol{X}^{k+1}, \boldsymbol{S}^{k+1}, {\Lambda}^{k+1} $利用式(8)(9)(10)

更新参数$ \mu:=\min \left(\rho \mu, \mu_{\max }\right)$

检查收敛条件

$\left\|\boldsymbol{Y}-\boldsymbol{L}^{k+1}-\boldsymbol{S}^{k+1}\right\|_{\mathrm{F}} /\|\boldsymbol{Y}\|_{\mathrm{F}} \leqslant \boldsymbol{\varepsilon}_{1} $和$\left\|\boldsymbol{L}^{k+1}-\boldsymbol{X}^{k+1}\right\|_{\infty} \leqslant \varepsilon_{2} $

为了验证模型的去噪效果以及模型的普适性和推广性，选择了两幅真实采集的高光谱集进行实验。第1种高光谱数据是Pavia Centre，Pavia Centre空间尺寸为1 096×1 096像素，空间分辨率为1.3 m，共含102个有用波段(http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes#opennewwindow)。Pavia Centre中的部分样本不包含任何信息，在分析前进行舍弃，且由于篇幅所限，本文主要对空间尺寸为256×256像素的高光谱图像块进行实验。第2种则是Washington DC National Mall(DC)高光谱图像块，空间尺寸为1 280×307像素，空间分辨率可达1 m，其中含有191个可用波段(https://engineering.purdue.edu/~biehl/MultiSpec/hyperspectral.html)。本文主要针对该数据的3个子集进行实验。其中两个空间图像的尺寸均为256×256像素，分别记为DC_sub1、DC_sub2，为了测试模型在不同尺寸下的效果，再截取尺寸为307×307、360×307、412×307以及512×307像素的图像，分别记为DC_dsize1、DC_dsize2、DC_dsize3、DC_dsize4。首先进行归一化处理，随后将不同种类的模拟噪声添加到数据中，得到模拟含噪图像。实验中选取的噪声分别为高斯噪声、椒盐噪声和死线噪声。本文选取的对比方法有：NAILRMA(noise-adjusted iterative low-rank matrix approximation)方法^[21]、LRMR(low-rank matrix recover)方法^[22]、LRTV (total variation-regularized low-rank Matrix factorization)方法^[7]。为了比对恢复效果，本文选取了SSIM和PSNR两个数值指标来比较不同方法的复原效果。

在Pavia Centre数据中添加等强度的高斯噪声和椒盐噪声。高斯白噪声的方差强度和稀疏椒盐噪声百分比强度分别记为$G \in\{0.04, 0.08, 0.12\} $ $P \in\{0.10, 0.15, 0.20\} $。对这3种不同的遥感数据分别使用本文方法和对比方法进行复原。同时将图像块复原结果每个通道的SSIM和PSNR分别取平均，记为MSSIM和MPSNR，作为复原效果的数值评判标准。表 1即为复原模拟下的数值指标。这两种混合噪声。图 1(a)绘制的是第97波段的原始图像，图 1(b)则是添加$ G$ =0.04, $P $ =0.10噪声强度后的图像，图 1(c)，图 1(d)，图 1(f)和图 1(e)依次是采用LRTV、TNN-LRTV、LRMR和NAILRMA去噪后的图像。可以看出LRTV等方法去噪后的图像上仍然残留有椒盐噪声，而本文方法去噪更彻底。

为了进一步比较去噪效果，图 2绘制出LRTV方法和本文方法TNN-LRTV在各波段的PSNR曲线图。可以看出，在本文方法去噪后几乎所有的通道PSNR值都高于LRTV方法，证明了该方法的优越性。

选取DC数据中的DC_sub1和DC_sub2数据集进行实验。在每个光谱通道中添加等强度的高斯或椒盐噪声或两者的混合噪声。高斯白噪声的方差强度和稀疏椒盐噪声百分比强度均分为3个级别。再次使用4种方法进行模拟复原。

首先选取DC_sub1，向其中分别添加强度不同的高斯白噪声或者稀疏椒盐噪声，模拟复原后得到表 2数据。可以看出，只有单独噪声时，本文方法比其他方法稍有提高。但是在实际观测中，通常是多种噪声混合在一起。为了进一步验证本方法的优越性，在DC_sub1中添加不同强度的高斯白噪声和椒盐噪声，表 3汇总了4种方法在DC_sub1中添加3种混合噪声下的数值指标。

可以看出，由于本文方法对于秩函数有更加精确的逼近，所得的MPSNR值和MSSIM值都超过其他方法。图 3表示在DC_sub1数据的每个通道中添加方差为0.09的高斯白噪声与比率为0.18的椒盐噪声。图 3(a)显示了第104个通道的原始图像，而图 3(b)是其加噪声之后的含噪图像，由于此时噪声比较高，整幅图几乎无法辨认出原有的图形。而图 3(c)-(f)依次为LRTV、TNN-LRTV、LRMR和NAILRMA的去噪后图像。可以看出，传统方法处理后仍残留着高斯白噪声，而本文方法处理后的图像有效抑制了混合噪声。除了图像直观显示，每个通道的PSNR值和SSIM值也证实本文方法在较高的噪声下，对图像的复原能力更强，并且没有异常值出现，见图 4和图 5。因此可以认为利用截断核范数代替LRTV模型中的核范数能够显著提升模型的复原效果。

为进一步比较4种复原方法的性能，对复原前后高光谱图像中某个点的光谱特征进行了对比，选取子图 1在像素点(100, 120)的光谱特征。图 6给出了原始图像在该点的光谱特征，以及该点的噪声图像。模型对比了4种方法恢复所得的该点光谱信息与原始图像中该点的光谱信息。可以看出，在各个波段，TNN-LRTV都具有更好的恢复效果，所得的光谱特征曲线基本和原始光谱曲线一致。而其他方法如LRTV则在某些波段有较大的偏离。本文方法能够更加有效地恢复像素点的光谱特征，从而达到更好的去噪效果。

最后，在DC_sub2和DC_dsize1中添加不同强度的高斯白噪声和椒盐噪声，表 4和表 5分别汇总了4种方法添加3种混合噪声的数值指标并绘制了各方法在各波段的PSNR对比图。此外，为了测试在不同尺寸下该模型处理噪声的能力，表 6汇总了添加方差为0.08的高斯白噪声与比率为0.15的椒盐噪声下的不同尺寸图像各方法得到的数值指标。在新的数据集下实验，本文方法仍然高于其他方法。特别在噪声强度较强、图像尺寸较大时，TNN-LRTV的恢复效果远高于其他的噪声。说明本文方法具有很强的推广性与普适性，同时对处理不同尺寸下的图片也均具有更好的效果。

首先对DC数据中的DC_sub1以及DC_dsize1进行随机噪声添加，使每一通道的噪声强度都不同。

表 7汇总了4种方法在添加随机高斯噪声和随机椒盐噪声不同尺寸下复原所得的数值指标。可以看出，在处理随机噪声的能力上，本文方法与其他3种方法相比，均有显著提高。

为了证明本文方法的有效性，在上述噪声的基础上，再向DC数据每个通道中添加死线噪声。死线噪声的强度是0到20间的随机数，代表图像行数的百分比。表 8汇总了4种方法添加随机死线噪声、高斯噪声和椒盐噪声不同尺寸下复原结果。综合表 7和表 8可以看出本文方法对于在不同尺寸下处理多种混合随机噪声均具有显著优势。

1) 主要研究内容：高光谱影像可以同时提供空间域和光谱域的信息，且高光谱影像在空间域和光谱域的连续性使其在相邻的通道间的相关性很强。基于这一特性，高光谱数据的复原均考虑采用光谱低秩先验或是空间低秩先验来建立复原模型。将两种模型结合使用可以取得更好的图像复原效果。考虑到现今在光谱低秩先验模型中核范数不能很好地利用到特征值的性质，首先引入了截断核范数的相关内容，然后利用截断核范数代替核范数改进LRTV模型并运用ADMM算法求解该模型。2)结果与分析：通过实验发现，基于截断核范数的LRTV模型在重构矩阵的秩和误差的稀疏性方面更具有鲁棒性，且在处理高密度噪声时效果更加明显，针对不同尺寸下的图像处理均能取得很好的重建结果，表明该模型能够复原出高质量的高光谱遥感数据。3)不足与展望：新方法在光谱低秩中更加准确地表示了观测数据的先验特征，但是在空间低秩先验方面只是沿用全变差正则化的模型。对全变差正则化进行进一步改进，再与截断核范数相结合，是今后研究的重点方向。