在计算机视觉和机器学习领域中，主成分分析(principal component analysis, PCA)一直是备受关注的方法，旨在分析高维数据的低维表示^[1]。PCA是一种非参数分析，能够恢复被稀疏噪声轻微损坏的数据。在采集和传输视觉数据时，不可避免地会引入异常值或大量噪声，而PCA对这些异常值或噪声非常敏感而不够稳定。为解决PCA的低稳定性问题，提出了鲁棒主成分分析(robust principal component analysis, RPCA)^[2]，旨在从被稀疏噪声破坏的观测矩阵中恢复低秩矩阵。其中矩阵的秩是唯一的，由非零奇异值的个数决定。由于矩阵的秩函数是非凸函数难以凸优化求解，所以RPCA采用核范数作为矩阵中秩函数的凸代替。假设观测矩阵$ {\boldsymbol{X}} \in {\mathbb{R}^{{n_1} \times {n_2}}} $可以分解得为低秩矩阵$ {\boldsymbol{L}} $和稀疏(噪声)矩阵$ {\boldsymbol{E}} $，则可通过奇异值阈值(singular value thresholding, SVT)^[3]方法进行有效求解。如今，RPCA和其拓展模型已经应用到很多领域，如图像恢复^[4]及对齐^[5-6]、前景检测^[7]、子空间聚类^[8-9]和人脸识别^[10]等。然而RPCA只能处理二维数据，但在实际应用中，高维数据的应用场景更为广泛且迅速增加。张量是矩阵的拓展，更适合表示高维数据。例如，彩色图像是大小为高×宽×通道数的三维张量，而灰度图视频是具有列、行和时间模式的三阶张量。为了处理这些张量数据，可以对张量的每个正面切片分别应用RPCA方法，但这种方式会忽略张量背后的多维结构信息，因此研究人员考虑将RPCA扩展到张量域中称为张量鲁棒主成分分析(tensor robust principal component analysis, TRPCA)。

假设给定一个观测张量数据$ \mathcal{{\boldsymbol{X}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，同理，TRPCA旨在从张量$ \mathcal{{\boldsymbol{X}}} $中精确地恢复$ \mathcal{{\boldsymbol{L}}} $和$ \mathcal{{\boldsymbol{E}}} $。这里，与矩阵秩是唯一的性质不同，张量秩的定义有很多种，导致张量分解的方法不同。Tucker秩^[11]是由Tucker分解所定义的，表示沿给定张量的每个模态展开的秩向量。Tucker秩只探讨张量的一个模态和其余模态之间的相关性。事实上，每个矩阵$ {{\boldsymbol{A}}^{(i)}} $都是通过将张量沿一个单模态矩阵化而得到的。因此，它是不平衡的，相应的秩不够大，无法捕捉张量中元素之间的全局相关性^[12]。由于直接最小化Tucker秩是多项式复杂程度的非确定性问题（non-deterministic polynomial，NP）困难的，所以核范数之和(sum of nuclear norm, SNN)^[13]被表示为Tucker秩的凸松弛。Huang等^[14]基于SNN提出了SNN-TRPCA模型，利用了张量子空间沿各模态的低秩性。然而由于实际数据中每个模态的低秩程度通常不同，因此很难设置正则化参数$ {\lambda _i} $。例如，一个灰度图视频沿其时间模式的秩比沿其空间模式的秩低得多。此外，单模态展开会破坏张量固有的结构信息。

因此近年来，张量平均秩^[15]受到了更多的关注，其定义为张量奇异值分解(tensor singular value decomposition, t-SVD)^[16]，其中块循环矩阵化以循环的方式排列张量的切片。研究表明，与Tucker秩相比，张量平均秩能保留更多的切片结构信息。由于张量平均秩的最小化同样是NP困难的，因此张量核范数(tensor nuclear norm, TNN)被推导为张量平均秩的凸代替。目前基于t-SVD，研究人员提出了不同的TNN定义^{[15, 17-18]}。特别是Lu等^[15]推导出了一个新的张量核范数，并给出了严格的理论保证：张量核范数是张量谱范数单位球内张量平均秩的凸包络，该模型可以通过乘法器交替方向进行求解^[19]，其中张量奇异值阈值法(tensor singular value thresholding，t-SVT)是求解的关键步骤。为了实现凸性，t-SVT根据阈值相等地缩小每个奇异值，忽略了关于奇异值的先验条件^[5]。此外，TNN-TRPCA有一个潜在的局限性，只是简单地假设整个底层张量是低秩的，对于视觉数据（如自然图像和视频）这样的假设通常难以满足。因此，TNN-TRPCA不能很好地恢复视觉数据中的细节信息，特别是在结构复杂的数据中。

受到重加权$ {l_1} $最小化方法^[20]的启发，研究人员借助非凸正则化器更精确地评估矩阵的奇异值^[21-24]。并且，越来越多的理论和应用证明了非凸优化方法的优越性^[25-27]。Gao等^[28]将加权张量Schatten-p范数最小化(weighted tensor Schatten-p norm minimization, WTSNM)和t-SVD结合起来求解他们提出的模型(enhanced tensor robust principal component analysis, ETRPCA)，并开发了一种高效的算法。然而，张量tubal秩在表征高阶张量之间的相关性方面并不是很有效^[29]。

针对上述问题，本文基于Gao等^[28]的方法进行改进，提出一种增强型TRPCA模型。首先引入加权非凸型张量Schatten-p范数对张量奇异值进行不同等级的收缩，可自适应阈值来收缩不同的张量奇异值，利用张量的全局低秩性更好地保留张量数据的重要信息。针对张量数据受损严重的问题，为了减少噪声干扰并简化数据结构，本文提出了一种可低秩预分离的方法，可在大量噪声干扰下提前将数据分离成近似的低秩结构和近似的稀疏结构，以便后续算法可分别逐次提高噪声分离效果。这样，可解决张量数据受损严重的问题，提高被噪声破坏的观测数据的恢复精确度。同时本文还提出了随机抖动机制，解决了高阶张量之间相关性弱的问题。该机制利用噪声信息的随机性削弱周期性噪声的影响，实现算法的正则化，并约束了模型复杂度。为求解所提模型，本文基于交替方向乘子法(alternating direction method of multipliers, ADMM)进行优化，将多优化变量问题转化为双闭式解问题简化处理。最后，进行了一系列的实验验证，将所提模型应用到高维数据恢复任务中，包括彩色图像恢复、核磁共振图像恢复、高光谱图像(hyperspectral images, HIS)及多光谱图像(multispectral image, MSI)恢复、灰度视频恢复等多个应用领域。在公开使用的标准数据集上对本文所提出的算法进行了评估和比较。实验结果表明，本文所提出的算法明显优于其他先进的TRPCA算法，并且随着受损数据包含噪声程度的不断增加，本文算法的恢复效果更好。

1.   基础知识和符号

1.1   符号说明

本文引入关于张量的符号及相关定义，采用大写粗体花式字母表示张量，如$ \mathcal {\boldsymbol{A}} $；采用大写粗体字母表示矩阵，如$ {\boldsymbol{A}} $；采用小写粗体字母表示向量，如$ {\boldsymbol{a}} $；采用小写字母表示标量，如$ a $。对于一个三维张量$ \mathcal {\boldsymbol{A}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$ {a_{ijk}} $表示$ \mathcal {\boldsymbol{A}} $中的每一个元素，即$ \mathcal {\boldsymbol{A}} $自身。$ \mathcal {\boldsymbol{A}}(i,:,:) $$ \mathcal {\boldsymbol{A}}(:,j,:) $、 $ \mathcal {\boldsymbol{A}}(:,:,k) $分别表示第i个水平切片、第j个侧面切片和第k个正面切片。通常，用$ {{\boldsymbol{A}}^{(i)}} $表示$ \mathcal {\boldsymbol{A}} $的第i个正面切片，$ \mathcal{{\boldsymbol{A}}}(i,j,:) $表示$ \mathcal {\boldsymbol{A}} $的tube。$ \overline {\mathcal{{\boldsymbol{A}}} }$表示$ \mathcal {\boldsymbol{A}} $沿着第三维度的快速傅里叶变换(fast Fourier transform, FFT)，如$ \overline{ \mathcal{{\boldsymbol{A}}}} = {\mathrm{fft}}(\mathcal{{\boldsymbol{A}}},[{\text{ }}],3) $。同样地经过快速反傅里叶变换(fast inverse Fourier transform, IFFT)有，$ \mathcal{{\boldsymbol{A}}} = {\mathrm{ifft}}(\overline{ \mathcal{{\boldsymbol{A}}}},[{\text{ }}],3) $。

定义1　（F范数^[16]）对于三维张量$ \mathcal{{\boldsymbol{A}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$\mathcal{{\boldsymbol{A}}} $的Frobenius范数为$ \parallel \mathcal{\boldsymbol{A}}{\parallel _F} = \sqrt {{\Sigma _{ijk}}|{a_{ijk}}{|^2}} $。

定义2　（对角矩阵^[16]）对于三维张量$ \mathcal{{\boldsymbol{A}}}\in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，存在对角矩阵$ \overline {\boldsymbol{A}} $，其对角线上的每个块是$ \overline {\mathcal{\boldsymbol{A}}} $的正面切片$ {\overline {\boldsymbol{A}}^{(i)}} $，有

定义3　（块循环矩阵^[16]）对于三维张量$\mathcal{{\boldsymbol{A}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，其块循环矩阵是一种张量矩阵化方式，大小为$ {n_1}{n_3} \times {n_2}{n_3} $，形式为

定理1^[17]　块循环矩阵可通过下式块对角化：

式中：$ \otimes $为克罗内克积，$ {{\boldsymbol{F}}_{{n_3}}} $为大小是$ {n_3} \times {n_3} $的离散傅里叶变换(discrete fourier transform，DFT)矩阵，$ {{\boldsymbol{I}}_{{n_1}}} $和$ {{\boldsymbol{I}}_{{n_2}}} $分别为大小是$ {n_1} \times {n_1} $和$ {n_2} \times {n_2} $的单位矩阵。

定义4　（折叠和展开操作^[16]）对于张量$\mathcal{{\boldsymbol{A}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，有

1.2   张量积和张量奇异值分解

定义5　（张量积(t-product)^[16]）令$\mathcal{{\boldsymbol{A}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$ \mathcal{{\boldsymbol{B}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，则张量积$\mathcal{{\boldsymbol{A}}} \mathit{*} \mathcal{{\boldsymbol{B}}}$的运算结果为$\mathcal{{\boldsymbol{C}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $。

$ {n_1} \times {n_2} \times {n_3} $的三维张量可以看作是$ {n_1} \times {n_2} $的二维矩阵，矩阵中的每一个元素沿着第三维度方向，称为tube。由此可见，张量积是矩阵积运算的拓展，使用循环卷积替代元素间的乘法运算。当$ {n_3} = 1 $时，张量积退化为标准矩阵乘法。因此，使得矩阵RPCA模型成为TRPCA模型的一个特例。

$\mathcal{{\boldsymbol{A}}} $与$\mathcal{{\boldsymbol{B}}} $的张量积可以通过以下步骤有效地计算：

1）计算$ \overline {\mathcal{\boldsymbol{A}}} = {\rm{fft}}(\mathcal{\boldsymbol{A}},[{\text{ }}],3) $，$ \overline {\mathcal{\boldsymbol{B}}} = {\rm{fft}}(\mathcal{\boldsymbol{B}},[{\text{ }}],3) $。

2）将$ \overline {\mathcal{\boldsymbol{A}}} $与$ \overline {\mathcal{\boldsymbol{B}}} $的每一个正面切片相乘得到$ \overline {\mathcal{{\boldsymbol{C}}}} $。

3）计算$ \overline{ \mathcal{{\boldsymbol{C}}}} ={\mathrm{ fft}}(\mathcal{{\boldsymbol{C}}},[{\text{ }}],3) $。

定义6　（共轭转置^[16]）对于三维张量$\mathcal{{\boldsymbol{A}}} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，首先将$\mathcal{A} $的每个正面切片进行共轭转置，再将转置运算后序号从2到$ {n_3} $的正面切片按逆序排列（第一个正面切片的位置固定不变）。即得到$\mathcal{A} $的共轭转置为$ {\mathcal{\boldsymbol{A}}^{\rm{T}}} \in {\mathbb{R}^{{n_2} \times {n_1} \times {n_3}}} $。

定义7　（单位张量^[16]）$\mathcal{I} \in {\mathbb{R}^{n \times n \times {n_3}}} $为单位张量，它的第一个正面切片为大小为$ n \times n $的单位矩阵，其余正面切片均为零矩阵。

定义8　（正交张量^[16]）若$\mathcal{Q} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $满足$ {\mathcal{\boldsymbol{Q}}^{\rm{T}}} * \mathcal{\boldsymbol{Q}} = \mathcal{\boldsymbol{Q}} * {\mathcal{\boldsymbol{Q}}^{\rm{T}}} = \mathcal{I} $，则$\mathcal{Q} $是正交的。

定义9　（f-对角张量^[16]）如果三维张量的每一个正面切片都是对角矩阵，则称为f-对角张量。

定理2　（张量奇异值分解(t-SVD)^[16]）令$\mathcal{A} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$\mathcal{A} $可以被分解为

式中：$ \mathcal{U}\in {\mathbb{R}^{{n_1} \times {n_1} \times {n_3}}} $，$ \mathcal{V} \in {\mathbb{R}^{{n_2} \times {n_2} \times {n_3}}} $，二者为正交张量；$ \mathcal{S} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $为f-对角张量。

图1直观地给出了t-SVD分解。在数学计算软件中可调用内置函数，使用傅里叶变换快速计算t-SVD。实现t-SVD的步骤如下：

图  1  ${n_1} \times {n_2} \times {n_3}$张量的t-SVD示意

Fig.  1  Illustration of the t-SVD of an ${n_1} \times {n_2} \times {n_3}$ tensor

下载: 全尺寸图片

1） 对$\mathcal{A} $沿第三维度进行傅里叶变换，得到$ \overline {\mathcal{\boldsymbol{A}}} $。

2） 对$ \overline {\mathcal{\boldsymbol{A}}} $的每个正面切片进行矩阵奇异值分解(SVD)，得到分解后对应的3个张量$ {\overline {\mathcal{U}} ^{(i)}} $、$ {\overline {\mathcal{S}} ^{(i)}} $、$ {\overline {\mathcal{V}} ^{(i)}} $。

3） 对3个张量进行反傅里叶变换，得到t-SVD的最终结果，即$ \mathcal{U} $、$ \mathcal{S} $、$ \mathcal{V} $。

需要注意，矩阵SVD分解是在复数域上进行的，最后的反傅里叶变换后得到的必须是实数。

定义10　（张量tubal秩(Tensor tubal rank)^[15]）对于三维张量$ \mathcal{A} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，其张量tubal秩记$ {\mathrm{rank}}_t(\mathcal{\boldsymbol{A}}) $，定义为$ \mathcal{S}$的非零奇异管道数量，可写为

式中$ \mathcal{S} $来自t-SVD中$ \mathcal{\boldsymbol{A}} = \mathcal{U} * \mathcal{S} * {\mathcal{V}^{\rm{T}}} $。张量tubal秩是非凸的。

定义11　（张量核范数^[15]）给定$\mathcal{X} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$ l = \min ({n_1},{n_2}) $，其核范数为

张量核范数为$ \mathcal{S} $的第一个正面切片的对角线上奇异值之和，该类凸模型的最大优点是可以计算最优解，但是张量核范数只是张量秩的一种凸优化逼近，在多数情况下会导致这种近似可能太松弛。因此，本文考虑引入一种非凸函数来更好地逼近张量秩，以获得更好的局部解。

2.   增强型张量鲁棒主成分分析模型

2.1   模型的建立与目标

TRPCA是一种典型的张量低秩表示方法。它可从被损坏的张量数据 $\mathcal{X} $中精确地恢复低秩张量数据$\mathcal{L} $，通过下述张量核范数模型^[15]：

式中：$\mathcal{X} $为张量噪声数据；$ \lambda $为平衡项的正则化参数；$\mathcal{L} $为低秩恢复张量，其低秩性由$ {l_ * } $范数描述；$\mathcal{E} $为稀疏恢复张量，其稀疏性由$ {l_1} $范数描述。

根据定义11，在处理张量核范数最小化问题中，所得模型(1)对输入的张量数据所有奇异值平等处理^[30-31]，当数据中存在一些异常值或离群值时，会导致这些异常值对结果产生过大影响，降低算法的鲁棒性，使得模型对数据中的异常值或噪声更为敏感。对所有奇异值进行同等收缩可能导致重要信息丢失，不同的奇异值代表了数据中不同的重要特征，同等收缩会使得这些特征之间的差异变得不那么明显，从而难以区分和识别数据中的重要特征。

在实际应用中，较大的奇异值通常与图像的颜色信息等一些突出信息相关联，较小的奇异值通常含有过高的噪声信号或者异常值。而TRPCA在最小化张量核范数时并没有充分利用奇异值的先验信息。因此，它不能很好地保留一些突出的信息，如颜色信息。为解决这一问题，张量奇异值部分和(tensor partial sum of singular values, TPSSV)被提出，其目标函数为

式中$ \parallel \mathcal{L}{\parallel _{p = r}} = \displaystyle\sum\limits_{i = 1}^{{n_3}} \parallel {\overline {\mathcal{L}}^{(i)}}{\parallel _{p = r}} $。由式(2)可知，在解决张量核范数最小化问题时，它不会收缩前r个最大的奇异值，这表明与前r个最大奇异值相关的信息与图像内容无关。然而由于这种假设条件极其严格，在实际当中并不可行。

受上述启发，为了保留主要信息同时减小冗余噪声干扰，应该使较大的奇异值收缩减小，同时使较小的奇异值收缩增大。因此，引入加权张量Schatten-p范数。

定义12　（加权张量Schatten-p范数）对于给定的$ \mathcal{X} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$ h = \min ({n_1},{n_2}) $，$\mathcal{X} $的加权张量Schatten-p范数为

由此，本文提出一种新的TRPCA模型，其目标函数为

式中：$\mathcal{X} $为张量噪声数据；λ为平衡项正则化参数；$\mathcal{L} $为低秩恢复张量，其低秩性由$ {l_{{\mathrm{Schatten}} - p}} $范数描述；$\mathcal{E} $为稀疏恢复张量，其稀疏性由$ {l_1} $范数描述。

2.2   算法优化

本节将介绍基于噪声预分离与随机抖动正则器的加权张量非凸型Schatten-p范数的TRPCA模型算法，该方法可应用于高维数据恢复。这些正则项之和最小化问题的目标是找到一个最理想的损失函数满足惩罚函数与目标函数之间求解出最优闭式解的问题。通常，给定一个受损张量数据$\mathcal{X} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，惩罚函数$ {{\Gamma }}$，目标函数$ {\Gamma _0} $，找到最小目标函数：

式中：损失函数$ \ell $比较了惩罚函数$ \Gamma $和目标函数$ {\Gamma _0} $，常数$ C $权衡损失函数和正则化项。

最小目标函数的求解取决于损失函数和正则化项的选择，具体过程描述如下。

2.2.1   加权张量非凸型Schatten-p范数最小化问题

前面提到本文使用加权张量非凸型Schatten-p范数的TRPCA模型(3),目标函数是需要获取$\mathcal{L} $的tubal秩最小化和$\mathcal{E} $的$ {l_1} $范数最小化。假设它们满足稀疏先验性，但不一定具有低秩先验性。

在TRPCA最小化问题当中，主要研究的问题是如何从一个被稀疏噪声损坏的张量数据中恢复分别得到低秩部分$\mathcal{L} $和稀疏部分$\mathcal{E} $。而需要重点关注低秩$\mathcal{L} $的恢复，其包含原始数据的主要信息，可近似看成对损坏数据的精确恢复数据。下面将介绍本文所提出方法的基本框架，并在下述部分展开描述详细的实现方法。

通过使用定理2中的t-SVD方法，本文给出ADMM的细节来求解式(3)。受增广拉格朗日函数的启发，可通过最小化下式来求解式(3)：

式中：$\mathcal{Y} $为增广拉格朗日乘子，$ \mu $为避免线性约束的惩罚参数，为正标量。

同时优化所有这些变量是很困难的。因此可以通过最小化一个变量而固定其它变量来近似地解决这个优化问题。这个过程就是所谓的交替方向乘子法ADMM。在多块ADMM框架下，通过交替最小化增广拉格朗日目标函数$ {\bf{\Gamma}} $来更新$\mathcal{L} $和$\mathcal{E} $，在特定的权重排序下，式(5)对各变量的优化问题可以等价地转化为独立的非凸$ {l_{{\mathrm{Schatten}} - p}} $范数和凸优化$ {l_1} $范数2个子问题，且2个子模型均有闭式解。利用ADMM得到的2个子模型及迭代增广拉格朗日乘子项为

对此本文给出了解决以上2个子模型的求解细节：

$\mathcal{E} $子模型　固定$ {\mathcal{L}_{k + 1}} $，$ {\mathcal{Y}_k} $，$ {\mu _k} $。根据软阈值算法，可将式(7)改写为

式中$ {D_\tau }(x) $是软阈值算子^[32]，作用是改变稀疏部分$\mathcal{E} $的数值。其定义为

$\mathcal{L} $子模型　同样地，固定$ {\mathcal{E}_k} $，$ {\mathcal{Y}_k} $，$ {\mu _k} $。为求解式(6)，可使用广义软阈值算法(generalized soft-thresholding algorithm，GST)^[33]求解其全局最优解。通过以下引理和定理展开讲解。

定理3^[28]　存在Y$ \in {\mathbb{R}^{{n_1} \times {n_2}}} $,$ r = \min ({n_1},{n_2}) $,$ \tau > 0 $,$ 0 \leqslant {\omega _1} \leqslant {\omega _2} \leqslant \cdots \leqslant {\omega _r} $,令Y的SVD为$ {\boldsymbol{Y}} = {\boldsymbol{{U_Y}{\Sigma _Y}{V_Y}}}^{\rm{T}} $，且$ \Sigma ={\mathrm{ diag}}({\sigma _1},{\sigma _2}, \cdots ,{\sigma _r}) $，假设所有奇异值均为非升序，则下述模型

根据Von-Neumann的迹不等式^[34]，可以找到式(10)的一个闭型全局最优解为

且$ {\Delta _{\tau * \omega [{\boldsymbol{Y}}]}} = {\mathrm{diag}}({\delta _1},{\delta _2}, \cdots ,{\delta _i}) $，其中$ {\delta _i} $由下式给出：

即使将模型(10)变换为简化后的式(11)，由于目标函数的非凸性和非光滑性，加上附加的阶约束(即$ {\delta _i} \geqslant {\delta _j},i \leqslant j $)，求解该问题仍然是困难的。在不失一般性的前提下，由于奇异值$ {\sigma _i} \geqslant 0 $，因此可以放弃非负约束$ {\delta _i} \geqslant 0 $。很显然，如果式(11)中的约束可以被丢弃，那么问题就可以解耦成r个独立的子问题，从而降低求解难度，将通过下述引理进行求解。

引理1　对于下述优化模型：

对于给定的p和$ {\omega _i} $，存在一个特定的阈值：

由此可以得出如下结论：

1)当$ \sigma\leqslant\tau_p^{\mathrm{GST}}(\omega_i) $时，式(12)的最优解$ T_p^{\mathrm{GST}}(\sigma_i,\omega_i) $=0。

2)当$ \sigma > \tau_p^{\mathrm{GST}}(\omega_i) $时，式(12)的最优解$ T_p^{\rm{GST}}({\sigma _i},{\omega _i}) = {\mathrm{sign}}({\sigma _i})S_p^{\rm{GST}}({\sigma _i},{\omega _i}) $，其中$ S_p^{\rm{GST}}({\sigma _i},{\omega _i}) $可由$ S_p^{\rm{GST}}({\sigma _i},{\omega _i}) - \sigma + \omega p{(S_p^{\rm{GST}}({\sigma _i},{\omega _i}))^{p - 1}} = 0 $求解得到。

由此可求得式(10)中$ {\delta _i} $，其中$ {\delta _i} = T_p^{\rm{GST}}({\sigma _i},{\omega _i}) $。

定理4^[28]　假设存在$\mathcal{A} \in {\mathbb{R}^{{n_1} \times {n_2} \times {n_3}}} $，$ r = \min ({n_1},{n_2}) $，$ \tau > 0 $，$ 0 \leqslant {\omega _1} \leqslant {\omega _2} \leqslant \cdots \leqslant {\omega _r} $，令$ \mathcal{\boldsymbol{A}} = \mathcal{U} * \mathcal{S} * {\mathcal{V}^{\rm{T}}} $，给出模型：

然后，由定理3可得模型(13)的全局最优解为

式中：$ {\mathcal{S}_{\tau * \omega }}_{[\overline {\mathcal{\boldsymbol{A}}}]} $为一个三维张量，$ {\mathcal{S}_{\tau * \omega [{{\overline {\mathcal{\boldsymbol{A}}}}^{(i)}}]}} $为$ {\mathcal{S}_{\tau * \omega }}_{[\overline {\mathcal{\boldsymbol{A}}}]} $的第i个正面切片。

证明　在傅里叶域中，式(13)等效为

根据定义1，可得到

式中$ {\overline {\boldsymbol{X}}^{^{(i)}}} $为$ \overline {\mathcal{X}} $的第i个正面切片。

在式(14)中，每一个变量$ {\overline {\boldsymbol{X}}^{^{(i)}}} $都是独立的，因此可将其拆解为$ {n_3} $个独立的子问题考虑。对于第i个(i = 1, 2, ···, $ {n_3} $)子问题，有

由定理3可知，式(14)的全局最优解为

式中全局最优解$ {\overline {\boldsymbol{X}}^{^{(i)}}}^ * $是$ {\overline {\mathcal{X}}^{ * }} $的第i个正面切片。由于已经得到了所有子问题的全局最优解，根据1.1小节中介绍过的反傅里叶变换，可以得到优化模型(13)的全局最优解：

式中：$ {\mathcal{U}_{\mathcal{\boldsymbol{A}}}} = {\rm{ifft}}({\overline {\mathcal{U}}_{\mathcal{\boldsymbol{A}}}},[{\text{ }}],3) $，$ {\mathcal{V}_{\mathcal{\boldsymbol{A}}}} = {\rm{ifft}}({\overline {\mathcal{V}}_{\mathcal{\boldsymbol{A}}}},[{\text{ }}],3) $。

阐述完引理及定理，回到求解模型(6)。根据定理4可知，Λ子模型(6)的全局最优解为

下面使用算法1总结本文提出的加权张量Schatten-p范数TRPCA模型求解方法的整个优化过程。图2直观地给出了所提出模型针对高维数据恢复的流程。

图  2  加权张量Schatten-p范数TRPCA模型概述

Fig.  2  Overview of the weighted tensor Schatten-p norm TRPCA model

下载: 全尺寸图片

算法1　基于张量加权非凸型Schatten-p范数的TRPCA算法

输入　观测张量数据$\mathcal{X} $，权重向量$ \omega $，正则化参数$ \lambda $

输出　低秩恢复$\mathcal{L} $，稀疏恢复$\mathcal{E} $

初始化　$ {\mu _0} = 1 \times {10^{ - 4}} $,$ \rho = 1.1 $,$ k = 0 $,$ {\mu _{\max }} = 1 \times {10^{10}} $,$ \varepsilon =1 \times {10^{ - 8}} $,$ {\mathcal{L}_0} = {\mathcal{E}_0} = 0 $,$ {\mathcal{Y}_0} = 0 $

通过ADMM求解模型(3)：

while不收敛 do

1）根据模型(9)更新$ {\mathcal{E}_{k + 1}} $

2）根据模型(15)更新$ {\mathcal{L}_{k + 1}} $

3）根据模型(8)更新$ {\mathcal{Y}_{k + 1}} $

4）通过$ {\mu _{k + 1}} = \min (\rho {\mu _k},{\mu _{\max }}) $更新$ {\mu _{k + 1}} $

5）判断收敛性

$ \parallel {\mathcal{E}_{k + 1}} - {\mathcal{E}_k}{\parallel _F} \leqslant \varepsilon$,$\parallel {\mathcal{L}_{k + 1}} - {\mathcal{L}_k}{\parallel _F} \leqslant \varepsilon $,$\parallel {\mathcal{L}_{k + 1}} + {\mathcal{E}_{k + 1}} - \mathcal{X}{\parallel _F} \leqslant \varepsilon $

end while

2.2.2   低秩预分离

在上文中已知TRPCA是对受损数据去除噪声或异常值，通过降维方法张量奇异值分解t-SVD从中提取主要信息使得受损数据得以恢复。而恢复后的数据包括低秩分量$\mathcal{L} $和稀疏分量$\mathcal{E} $共2部分。所以本文考虑高维数据可能受到较大程度损坏的情况，对其进行预处理，使得其在进行TRPCA方法之前已经提前获得近似含有低秩信息和稀疏噪声信息的两分支，将噪声预分离后的2部分分别作为受损数据输入，再次使用TRPCA降维处理。预处理阶段通过识别和分离出数据中的噪声成分，从而使数据更加干净，减少噪声对数据分析结果的影响，有利于后续的分析和建模过程，在此基础上再进行降维处理相较之前更具有明确性和鲁棒性。

基于低秩预分离的方法可以预先处理一个含有大量噪声的受损数据对低秩恢复数据及异常数据的提前分离能力，其基本流程如下所述。

1）平滑化处理。为了减少噪声干扰并简化数据结构，本文考虑首先利用滤波算法^[35-36]（如中值滤波、高斯滤波和小波变换滤波等，针对不同类型的噪声数据可分别选择其相对适宜的滤波算法）对受损张量数据$\mathcal{X} $进行分解。可将首次滤波后结果视为低秩恢复数据，即得到噪声预分离后的近似低秩部分$ \mathcal{L}' $，实现初步筛选低秩结构。

2）作差处理。通过将受损张量数据$\mathcal{X} $与近似低秩部分$ \mathcal{L}' $做差法计算，可自动筛选出近似异常数据，即得到噪声预分离后的近似稀疏部分$ \mathcal{E}' $，其中$ \mathcal{E}' = \mathcal{X} - \mathcal{L}' $。

通过上述预分离的思想可初步筛选出受损严重的张量数据$\mathcal{X} $中的近似低秩部分$ \mathcal{L}' $和近似稀疏部分$ \mathcal{E}' $，实现低秩结构的预分离。该方法的具体过程如图3所示。

图  3  低秩预分离流程

Fig.  3  Flowchart of low-rank pre-separation

下载: 全尺寸图片

2.2.3   随机抖动机制

考虑到高维数据中第3维度相邻切片一致性较弱，提出随机抖动机制来增强一致性。如在连续视频帧中，使每一帧均加入随机的轻微抖动(jitter)，则信号或图像特征（如边缘和纹理）在相邻帧之间的位置差异会相对减小。由此，当对多帧进行平均或采用时域滤波时，随机出现且不随时间连续的噪点（其在不同帧之间不具有一致性）会被削弱，使得图像信息得以保留和加强。

张量数据受到噪声破坏时，通过TRPCA模型恢复低秩张量的过程中可能出现过拟合情况，加入随机抖动机制可破坏周期性噪声从而增加数据的随机性减少数据对特定噪声模式的依赖来削弱过拟合情况。jitter可改变噪声在图像中的分布，使得噪声分布不规律从而破坏噪声的可预测性。正是通过利用受损数据中噪声信息的随机性来正则化算法，以减少过拟合并且提高鲁棒性，有助于约束TRPCA模型的复杂度，可简化上文所提出的非凸TRPCA后续恢复过程。

受文献[37-38]的启发，本文提出的随机抖动正则器具有隐式形式。在将输入图像数据传送到TRPCA损失函数$ \ell $之前，可随机选择区域优化。由此可将式(4)改写为

式中：$ \ell $为惩罚函数$ \Gamma ( \cdot ) $和目标函数$ {\Gamma _0} $之间的损失计算；$ E[ \cdot ] $为数学期望；$ \varsigma = ({\varsigma _1},{\varsigma _2}) $为均匀分布在集合$ \left\{ 0 \right., \cdots ,{\left. {T - 1} \right\}^2} $中的离散随机变量，表示最多$ T - 1 $个像素的随机水平或垂直翻转。运算符$ {\mathrm{jitter}}( \cdot ) $对张量数据$\mathcal{X} $进行如下转换和分割：

式中：$ v $为索引图像行，$ u $为索引图像列，$ 1 \leqslant v \leqslant {n_1} - T + 1 $且$ 1 \leqslant u \leqslant {n_2} - T + 1 $。随机抖动正则化机制每次随机索引一个区域进行优化，期间可出现相互重叠情况（每次索引区域可不一致），经多次索引迭代过程后，将输入数据的完整区域全部进行索引，作为末次迭代终止项。随机抖动机制的示意如图4所示。

图  4  随机抖动机制示意

Fig.  4  Illustration of the random jitter mechanism

下载: 全尺寸图片

3.   实验结果与分析

在实验中将选择一些先进的经典TRPCA模型方法作为对比方法，包括基于克罗内克基表示的KBR方法^[39]、基于t-SVD的TNN方法^[15]、基于凹凸分离型t-SVD的CCSVS方法^[40]与本文所提方法进行比较。分别在彩色图像恢复、核磁共振图像恢复、高光谱及多光谱图像恢复、灰度视频恢复4种应用中展开实验。实验将从定量和定性两方面分析实验结果，采用原始数据、受损数据和恢复数据肉眼定性评估。采用峰值信噪比(peak signal-to-noise ratio, PSNR) $ {E_{{\mathrm{PSNR}}}} $（以dB为单位）、结构相似度(structural similarity, SSIM) $ {E_{{\mathrm{SSIM}}}} $、特征相似度(feature similarity, FSIM) ${E_{{\mathrm{FSIM}}}}$、CPU运行时间（以秒为单位）定量质量指标来评价方法的性能。前3个评估指标定义如下：

式(16)中：$ {\mathcal{L}_0} $为观测张量，$ {\mathcal{L}_0}(:,:,i) $为第$ i $帧图像，$ \mathcal{L} $为恢复张量，$ {\max _\mathcal{L}} $为$ \mathcal{L} $的最大像素值。式(17)中：$ {\mu _L} $和$ {\mu _{{L_0}}} $分别为图像$L$和${L_0}$的均值，$ {\sigma _{L{L_0}}} $为$L$和${L_0}$的协方差，$ {\sigma _L} $和$ {\sigma _{{L_0}}} $分别为$L$和${L_0}$的标准方差，$ {c_1} $和$ {c_2} $均为常数。式(18)中${S_L}(x)$是相位一致性特征(phase congruency, PC)和梯度特征(gradient magnitude, GM)的融合的相似度，$P{C_m}(x)$为衡量${S_L}(x)$在2张图像中全局相似性的重要程度，$\Omega $为整个图像空间域。对于PSNR值越高越好，SSIM值和FSIM值越接近1越佳。

3.1   彩色图像恢复

在本实验中，将验证并评估本文方法及其它对比方法对损坏的彩色图像的恢复效果。选择Berkeley Segmentation数据集作为图库，随机选择50张彩色图像进行测试。每张图像的大小为481×321×3像素。对于每张图像，人为添加随机噪声，分别设置随机噪声，使10%、20%、30%的像素被设置为[0, 255]中的随机值，并且损坏像素的位置未知。每个彩色图像包括3个颜色通道，即红、绿、蓝，可以看作是一个三阶张量。所有3个通道的图像都在相同的像素位置进行损坏。对于其他的对比方法，设置的参数与其参考原文相同。由于本文方法是在Gao等^[28]的方法上进行改进，初始化设置正则化参数$ \lambda =1/\sqrt {3\max ({n_1},{n_2})} $，$ {\mu _0} = 1 \times {10^{ - 4}} $，$ \rho = 1.1 $，$ k = 0 $，$ {\mu _{\max }} = 1 \times {10^{10}} $，$ \varepsilon = 1 \times {10^{ - 8}} $，均来自原文献。对于在不同噪声等级下可调参数p的选取，同样参考原文献，并在后续3.2小节中进行了参数敏感度的具体分析。通过仿真分析证明p值选取的合理性，并在具体应用任务当中加以微调，争取达到最优效果。低秩预分离平滑化过程选择中值滤波器，随机抖动正则器每次迭代选择的区域大小均为80像素×80像素，迭代次数为20次。

表1给出了所有对比方法在不同噪声等级下对于50张彩色图像的平均恢复结果。从这些结果可以看出，本文方法的3个评估指标均较对比方法有更好的性能，几乎在所有情况下都达到了最高的PSNR、SSIM及FSIM值。该方法比较为先进的TNN方法在PSNR上平均提高3.37~5.19 dB，优于仅使用低秩信息的一般方法。KBR与CCSVS方法随噪声等级的增加恢复效果明显变差，具有较差的鲁棒性，不如本文方法。但是由于本文方法每次迭代需要更多的计算量，因此在运行时间指标上表现出一般的性能。此外在数据中显示，当张量Schatten-p范数中的p=0.9时，本文方法在处理较低噪声等级的数据时获得了最好的平均结果；当p=0.95时，在较高噪声等级的数据恢复中具有更好的性能。由此可得出结论：随着噪声等级增加，其他方法降噪效果明显下降，本文方法仍可通过调节p的取值来保持较好恢复效果。

图5以50张彩色图像为例，分别列出了各种对比方法的PSNR、SSIM和FSIM值，并定量地证明本文的恢复效果优于所有其他方法。可看出在10%噪声等级下，本文方法p=0.9时3个评估指标均取得最优效果，p=0.95时在SSIM和FSIM评估指标中均取得次优效果。

图  5  彩色图像恢复结果的PSNR值、SSIM值及FSIM值比较(10%噪声等级)

Fig.  5  Comparison of the PSNR, SSIM and FSIM for color image restorations (10% noise level)

下载: 全尺寸图片

图6从50张彩色图像中任意选取5张图像为例，分别为10%噪声等级下本文方法与多种对比方法的原始图像、噪声图像及其恢复图像。从这些恢复图像中不难看出，本文方法的性能比其他方法好得多。原因是其他方法仅利用了低秩先验条件下数据的多维结构，而本文方法除了利用低秩先验之外，还利用了局部平滑先验。这比仅仅在矩阵上使用低秩先验要好。使用本文方法得到的恢复图像视觉更清晰，边缘更锐利，从而从视觉上证明了本文方法在测试数据集上取得良好的效果。总之，本文方法显示出强大的恢复能力，在可视化中产生更令人满意的去噪恢复结果，并且在定量上具有更高的PSNR、SSIM及FSIM指标。

图  6  10%噪声等级下的彩色图像恢复结果比较

Fig.  6  Comparison of color image restorations results at 10%noise level

下载: 全尺寸图片

3.2   MRI图像恢复

为了证明所提出模型针对高维数据恢复的有效性，在本小节中测试了MRI数据集，大小为181×217×40像素，并对其分别加入10%、20%、30%的随机噪声。本实验中本文方法及其他对比方法的参数配置与上一实验相同，对于可调参数p值的选取最终以实验结果和敏感性分析加以确定。

表2列出了不同噪声等级下MRI数据集各切片的平均PSNR、SSIM、FSIM值和运行时间。可以看出本文方法在噪声等级不断增加的情况下仍然保持较好的PSNR、SSIM、FSIM值。但由于本文方法需要较大的计算量，所以在运行时间上没有明显优势，KBR方法具有最短运行时间。从表2中明显可以看到当噪声等级为10%时，本文方法p=0.7取得最优效果；当噪声等级为20%时，本文方法p=0.8取得最优效果；当噪声等级为30%时，本文方法p=0.9取得最优效果。由此可见随着噪声等级的增加，可通过增大p的取值来保持本文方法的优势。

图7为10%噪声等级下，针对MRI数据集中的每个图像，提取10个相邻切片的PSNR和SSIM值对比。可直观地注意到本文方法明显优于其他对比方法。

图  7  MRI图像恢复中相邻10个切片的PSNR值及SSIM值比较(10%噪声等级)

Fig.  7  Comparison of the PSNR and SSIM of 10 adjacent slices of the MRI image restorations (10% noise level)

下载: 全尺寸图片

图8给出了调节参数p在MRI数据集上的灵敏度分析。通过对p值的敏感性分析，可以发现当$ \eta $介于1.1~1.4时，算法能够取得较好的性能表现。随着p值的减小，曲线的峰值有所右移，进一步表明$ \eta = 2.1 - p $的选择是合理的。为了验证这一选择的普适性，对不同p值进行了进一步的实验分析。结果表明，p值在0.7~0.9的设置与灵敏度$ \eta $的变化相得益彰，能够达到最佳的性能。因此，通过敏感性分析，最终确定了p值范围的最优设置。该参数组合在不同数据集和应用场景中均能稳定表现，具有较强的通用性，避免了依赖先验知识的局限性，能够更好地推广到实际应用中。

图  8  MRI参数$ \eta $调节的灵敏度

Fig.  8  Sensitivity for the parameter $ \eta $ tuning on MRI

下载: 全尺寸图片

图9给出了不同方法在噪声等级分别为10%、20%和30%的MRI图像中第23个切片及第40个切片的恢复结果，其中本文方法的恢复图像选取的为p取最优值时的结果。可明显看出KBR方法效果明显不佳，对比之下CCSVS方法逐渐占有优势，但随噪声等级增加，图像恢复后的轮廓及细节清晰度仍不及本文方法。本文方法恢复后的脑图纹路清晰度几乎高度还原于原始MRI图像。随着噪声等级的增加，其他对比方法的恢复效果明显降低，本文方法优势最佳。

图  9  不同噪声等级下MRI图像的恢复结果比较

Fig.  9  Comparison of MRI image restorations at different noise levels

下载: 全尺寸图片

3.3   高光谱及多光谱图像恢复

本实验分别采用一张高光谱图像和一张多光谱图像，选自高光谱遥感数据库中的Pavia Centre和哥伦比亚MSI数据库中的Cloth进行测试，大小分别为$200 \times 200 \times 80$像素和$512 \times 512 \times 31$像素，并对其分别加入10%、20%、30%的随机噪声。本实验中本文方法及其他对比方法的参数配置与上一实验相同。

表3列出了各种方法对高光谱或多光谱图像在不同噪声等级下各波段的平均PSNR、SSIM、FSIM值和运行时间。从表3中可以看出随着噪声等级的增加，本文方法均在PSNR、SSIM、FSIM值中取得最优效果，且增大p的取值可获得更好效果，CCSVS方法取得次优效果，KBR方法效果明显不如其他方法，但在运行时间方面占据显著优势。

图10选取不同噪声等级下的所有对比方法恢复的高光谱或多光谱图像的第23波段，其中本文方法的恢复图像选取的为p取最优值时的结果。可以清楚地看到，本文方法的恢复效果与CCSVS方法相似，观测数据中噪声被去除，并且图像中的纹理和边缘也得到了很好的保留，多光谱图像Cloth中KBR方法在各个噪声等级中肉眼可见的恢复能力不及本文方法。

图  10  不同噪声等级下对高光谱和多光谱图像的恢复结果比较

Fig.  10  Comparison of HIS and MSI restorations at different noise levels

下载: 全尺寸图片

3.4   灰度视频恢复

灰度视频序列本质上是一个三阶张量，每一帧都是一个灰度图像，可以用一个低tubal秩张量来近似。在本实验中，测试视频来源于YUR视频数据库,从中选取了2个视频进行展示,即Road和Suzie，测试视频的大小分别为158像素×238像素×24像素和144像素×176像素×150像素，并对其分别加入10%、20%、30%的随机噪声。本实验中本文方法及其他对比方法的参数配置与上一实验相同。

表4给出了各种方法对2个视频在不同噪声等级下每一帧的平均PSNR、SSIM、FSIM值和运行时间。从表4中可以看出不论哪一种噪声等级下，本文方法均在PSNR、SSIM、FSIM值中取得最优效果，KBR方法在运行时间方面占据显著优势。在不同噪声等级下，KBR方法恢复效果均明显低于其它方法，并且本文可通过调节增加p的取值来获得最好效果。

图11给出了不同噪声等级下的所有对比方法恢复的测试视频的第23帧，其中本文方法的恢复图像选取的为p取最优值时的结果。可以清楚地看到，Road视频中KBR方法恢复后有明显虚影，Suzie视频中KBR方法恢复后的人物轮廓模糊，在噪声等级较大时CCSVS方法仍存在部分噪声，本文方法的人物边缘清晰度优于TNN方法与CCSVS方法，且明显具有比KBR更好的性能特征。

图  11  不同噪声等级下灰度视频的恢复结果比较

Fig.  11  Comparison of grayscale video restorations at different noise levels

下载: 全尺寸图片

4.   结束语

本文提出了一种增强型TRPCA模型，该模型同时利用了低秩先验结构和局部平滑性。本文模型利用加权张量非凸型Schatten-p范数对张量秩进行松弛，通过区别对待奇异值来更大程度地保留重要信息同时去除稀疏噪声干扰。在面对受损张量数据噪声等级较大的情况，针对性地提出了低秩预分离的方法实现受损数据的近似低秩部分和近似稀疏部分的提前分离，同时利用随机抖动机制对受损数据中第三维度增加一致性且对噪声信息的随机性来正则化算法，以减少过拟合并且提高鲁棒性。针对所提出的模型利用ADMM的算法进行求解，并将该模型应用到高维数据恢复实验当中，如彩色图像恢复、核磁共振图像恢复、高光谱或多光谱图像恢复及灰度视频恢复4个应用。

大量实验表明，在多数情况下，本文模型在定性和定量评估指标上均取得不错效果，优于KBR、TNN、CCSVS这3种最先进的方法。但是所需CPU运行时间较长，这是由于基于Schatten-p范数的非凸优化问题需要对每个张量奇异值进行不同等级的收缩，并且常常需要更多的迭代次数来收敛从而产生更多运算量，消耗了较多时间。不过在处理第3通道数目相对较低的张量数据时是可以接受的等待时间范围。考虑到本文方法能够实现优越的图像处理效果，可以平衡所需的时间代价，未来考虑可在运行速度上进一步提高。