一种基于双层模型驱动的非凸张量方法，融合了低秩性与平滑性特性，用于图像去噪

《Applied Mathematical Modelling》：A two-layer model-driven nonconvex tensor method with low-rankness and smoothness fusion for image denoising

　　高凯欣|孟杰|张翔•我们提出了一种用于去除稀疏噪声的新型TNLS模型。•该模型通过非凸函数将全局低秩性与局部平滑性统一起来。•我们首次为这类非凸张量方法给出了统一的误差界。•我们设计了一种具有闭式解的ADMM算法，并证明了其收敛性。•实验表明，TNLS的性能优于现有的先进方法。1. 引言从含噪图像中恢复清晰图像是计算机视觉领域的核心任务[1]、[2]、[3]、[4]、[5]，这类图像可能在采集或传输过程中受到污染。本质上，这属于一个逆问题，由于问题的固有病态性，找到唯一且稳定的解极具挑战性[6]。该逆问题的基本数学形式如下：(1)y=A(x)+e，其中A:RM1→RM2（M2≤M1）为线性变换，x∈RM1代表我们希望恢复的清晰图像，y∈RM2表示观测到的图像，e∈RM2则为噪声。在大多数去噪任务中，A可视为恒等运算符，即A(x)=x，此时问题转化为从y=x+e中找出x。然而，在方程(1)中，图像的向量化表示通常会破坏其固有的空间结构，比如像素邻接关系以及局部区域的连续性，而这些结构对于准确恢复图像至关重要。为保留这些结构信息，实际研究越来越多地采用张量来表示多维图像[7]、[8]、[9]、[10]。基于张量表示，核心的张量去噪任务变为从Y=X+E中恢复X∈RI1×?×IN，其中张量E∈RI1×?×IN代表噪声，Y∈RI1×?×IN则是观测到的张量。本文聚焦于去除实际场景中常见的稀疏噪声。稀疏噪声的特点在于，其噪声观测值中只有少量元素为非零，这类噪声在视觉效果上往往表现为孤立的异常像素。为解决稀疏噪声去除的病态性问题并获得可靠的恢复结果，关于自然图像的有效先验信息非常重要[11]、[12]、[13]、[14]。低秩性是最常被采用的全球先验之一，它源于自然图像中存在的冗余结构信息。图像中的纹理、边缘及其他特征往往会出现重复，这使得清晰图像的张量表示具有明显的低秩特性。结合这一全球先验与噪声的稀疏性，研究人员开发出了多种有效的去噪方法，其中张量鲁棒主成分分析（TRPCA）就是一种代表性方法[4]、[7]、[15]、[16]。从数学上来看，TRPCA的模型形式为(2)minX,Erank(X)+λ∥E∥0,s.t.Y=X+E，其中rank(X)表示X的秩，∥E∥0为L0范数，即稀疏噪声张量E中非零元素的个数，λ＞0为平衡参数，用于调整两种正则化项的权重。要解决方程(2)，首先需要定义张量秩，与矩阵秩有唯一定义不同，张量秩存在多种不同的定义方式。最常用的定义包括CANDECOMP/PARAFAC（CP）秩[17]、Tucker秩[18]、训练秩[19]、管状秩[20]以及环秩[21]。不过，由于秩函数和L0范数的非凸性和不连续性，方程(2)通常属于NP难问题。应对这一计算难题的经典方法是使用凸替代函数。例如，所有模式展开矩阵的核范数之和[22]以及张量核范数[16]分别被用来近似Tucker秩和管状秩，而L1范数则被用来近似L0范数。Goldfarb和Qin[13]、Lu等人[16]、Liu等人[22]、Yang等人[23]的研究都体现了这种替代函数能让TRPCA模型具备凸性。尽管这些模型在数学上易于处理，但由于凸替代函数倾向于对较大系数施加更重的惩罚，因此并非最优选择[11]。为解决这种不平衡问题，TRPCA中采用了许多非凸函数来近似张量秩或L0范数。例如，加权张量Schatten p范数[15]、Lp（0＜p＜1）范数[4]以及eγ型惩罚项[24]都被用来实现对张量秩或稀疏性的近似。最近，高凯欣等人[7]基于Tucker秩提出了一种广义的两层非凸TRPCA框架，该框架涵盖了众多常用的非凸函数。此外，双层结构还能有效降低奇异值分解（SVD）的计算成本。除了TRPCA之外，Guo等人[8]、Chen等人[25]、Gao等人[26]、Tu等人[27]及其参考文献中还提出了更多张量秩的非凸替代函数。这些研究都证明了在实际应用中，非凸替代方法相比凸替代方法更具优势。虽然低秩性先验能够捕捉图像的全局结构特征，但反映自然图像局部内在特性的局部平滑性先验对于获得高质量图像恢复结果同样重要。总变分（TV）正则化是一种用于描述图像局部平滑性的经典且有效的方法[1]、[28]、[29]、[30]。通过最小化图像的TV值，TV正则化能够在去除噪声的同时保持张量的平滑性。TV正则化可分为两类：各向异性TV和各向同性TV[31]，前者以相邻元素差值的绝对距离作为定义依据，后者则采用平方距离。在实际应用中，还为提升TV正则化的性能设计了其他变体，如空间TV[1]、谱空间TV[29]、高阶TV[28]、[32]以及一阶TV和二阶TV的融合[33]。近年来，基于TV的方法被广泛应用于视觉张量数据恢复任务，而且这些应用大多采用结合了低秩正则化和TV正则化的联合模型[12]、[33]。通常，上述联合模型通过独立处理低秩正则化和TV正则化，能够有效捕捉图像的低秩性和平滑性。然而，这类模型大多缺乏理论支撑，且在调整参数以实现不同正则化方式之间的平衡方面面临挑战。因此，彭等人[9]提出了相关TV（CTV）正则化器，该正则化器能够以统一的方式同时纳入低秩性和平滑性。王等人[34]进一步发展了这种方法，他们在高阶张量SVD框架内将低秩性和平滑性进行了整合。Wang等人[10]、Feng等人[35]、Liu等人[36]也提出了其他类似改进方案。这类方法不仅具有理论依据，而且在图像修复和噪声去除等实际任务中也展现出显著的优势。目前基于CTV及其变体的图像噪声去除模型通常采用L1范数作为凸替代函数，用于表征目标张量分量的低秩性和稀疏性，但这些方法不可避免地存在前文所述的凸替代函数的固有局限性。高凯欣等人[7]建立了一种广义的两层非凸张量框架，该框架能有效降低稀疏噪声去除任务的计算成本。不过，这种非凸框架仅考虑了低秩正则化和稀疏性正则化，忽略了图像固有的局部平滑性先验。表1总结了部分相关方法中的先验和替代函数。可见，一些现有方法要么仅利用低秩性，要么仅利用平滑性；而在同时考虑两种先验的方法中，几乎所有方法都依赖凸替代函数。表1. TNLS与相关方法的先验及替代函数对比。方法先验替代函数低秩性平滑性凸非凸TRPCA [16]????ETRPCA [15]????3DCTV [9]????CTV [34]????NTRPCA [7]????RCCTV [37]????TNLS（本文提出的方法）????鉴于上述局限性，我们进一步研究了基于高凯欣等人[7]所提模型构建的模式展开矩阵。如图1所示，模式展开矩阵不仅能保持固有的低秩性，还具有显著的局部平滑性。受此启发，再加上非凸替代函数的优势以及将低秩性与平滑性相结合的潜力，我们在本文中提出了一种全新的双层非凸张量方法，该方法将低秩性与平滑性相融合，称为TNLS。具体而言，我们首先选取能够有效反映观测张量信息的模式展开矩阵，然后为该矩阵引入一个正则化项，该项同时兼顾低秩性和平滑性。此外，我们还利用一类广义非凸替代函数来表征低秩性和稀疏性，这使得我们的框架能够与大多数常见的非凸替代函数兼容。从理论上，我们为TNLS模型建立了最坏情况误差界，该误差界适用于广泛的非凸函数类别。此外，我们还设计了一种交替方向乘子法（ADMM）算法，该算法中每个变量都可以通过闭式解进行更新，同时我们还证明了该算法能够收敛到Karush-Kuhn-Tucker（KKT）点。在多种成像数据集上进行的大量数值实验都证明了所提方法的有效性。下载：下载高分辨率图像（606KB）下载：下载全尺寸图像图1. 以一张481×321×3色的图像和一张128×128×31的多光谱图像为例，展示了模式展开矩阵的低秩性和平滑性。从上到下，这些图像分别代表了原始图像的模式展开矩阵、对应的奇异值以及所有元素的频率直方图。本文的主要贡献如下：•模型设计。我们首次提出了一种新型TNLS模型，将该模型将低秩性和平滑性整合到了统一的非凸框架中。如表1所示，TNLS在广义非凸松弛条件下统一了低秩性和平滑性先验，而凸替代函数自然可作为所提公式的一种特例。•理论分析。我们建立了通用的理论分析框架，用于推导TNLS模型的最坏情况误差界。与局限于单一形式的分析不同，这一理论框架适用于众多类型的非凸替代函数。据我们所知，这是首个为将低秩性和平滑性统一起来的非凸张量方法提供的理论保障。•算法实现。我们设计了用于求解TNLS模型的ADMM算法，该算法中所有变量都有闭式更新表达式。我们还进行了收敛性分析，证实所提出的ADMM算法生成的序列能够收敛到KKT点。•实验验证。我们在多种类型的成像数据集上针对几种先进的算法进行了大量实验。数值指标和视觉对比结果均表明，TNLS的性能优于其他竞争方法。本文的其余部分结构如下：第2节介绍一些基本的符号和定义；第3节提出TNLS模型，并为其建立最坏情况误差界；第4节详细阐述ADMM算法的设计及其收敛性分析；第5节通过大量的数值实验结果来验证我们方法的有效性；最后，第6节对全文进行总结。2. 基础知识在本节中，我们将介绍本文所使用的所有基本符号、定义以及必要的引理。在本文中，小写字母（如x）表示标量，加粗小写字母（如X）表示向量，加粗大写字母（如X）表示矩阵，而手写体字母（如X）则表示张量。对于任意正整数n，我们将集合[n]定义为{1,2,?,n}。给定一个矩阵X∈RI1×I2，其奇异值分解可表示为X=USV?，其中U∈RI1×r和V∈RI2×r为列正交矩阵，S为对角矩阵，其表达式为diag(σ1(X),σ2(X),…,σr(X))，其中r=min{I1,I2}。奇异值按非递增顺序排列：σ1(X)≥σ2(X)≥…≥σr(X)。通常，diag(x1,…,xI)表示一个对角矩阵，其中对于RI中的向量x，其对角线上的第i个元素为xi。矩阵X的核范数记为‖X‖*，即为该矩阵所有奇异值之和，即‖X‖*:=∑i=1rσi(X)。梯度矩阵?1X,?2X∈RI1×I2的定义分别为?1X:=DI1X和?2X:=(DI2X?)?，其中(?·?)?表示转置运算符，而差分矩阵DI∈RI×I（I∈{I1, I2}）是一个I×I阶的行循环矩阵，其元素为(?1,1,0,…,0)。实际上，?1X和?2X分别用于表示矩阵X的垂直方向和水平方向的差分。对于一个N阶张量X∈RI1×I2×?×IN，我们可以用Xi1i2?iN来表示其第(i1,i2,…,iN)个元素，其中in∈[In]，n∈[N]。当我们沿着第n个模式展开该张量时，得到的矩阵被称为第n阶模式展开矩阵，记为X(n):=unfoldn(X)。然后，我们可以通过折叠操作，即X=foldn(X(n))，从第n阶模式展开矩阵X(n)恢复出原始张量X。对于任意n∈[N]，都有X(n)∈RIn×Kn，其元素(X(n))inj=Xi1i2?iN，其中Kn:=∏k=1,k≠nNIk,j=1+∑k=1,k≠nN(ik?1)Jk，而Jk=∏l=1l≠nk?1Il。对于两个张量X,Y∈RI1×I2×?×IN，它们的内积定义为〈X,Y〉:=∑i1=1I1∑i2=1I2?∑iN=1INXi1i2?iNYi1i2?iN。张量的几种常见范数定义如下：L0范数‖X‖0表示矩阵X中非零元素的个数；L∞范数为‖X‖∞:=maxi1,i2,…,iN|Xi1i2?iN|；X的L1范数和Frobenius范数分别为‖X‖1:=∑i1=1I1∑i2=1I2?∑iN=1IN|Xi1i2?iN|和‖X‖F:=〈X,X〉。设f:R→(?∞,+∞]为一个 proper且下半连续的函数。给定 y∈R 且 λ?>?0，定义 proxλf(y):=argminx∈R{λf(x)+12∥x?y∥22}。最后，我们给出以下引理。引理1 Bullen等人[38]（广义幂平均不等式）：设α1,α2,…,αn为n个正数，满足∑i=1nαi=1。对于任意c1,c2,…,cn≥0，设0??0为惩罚参数。ADMM的迭代格式如下：(23){Xk+1=argminX{L(X,Ek+Ak+Bk,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)},Ek+1=argminE{L(Xk+1,E,Ak,Bk,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)},Ak+1=argminA{L(Xk+1,Ek+1,A,Bk,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)},Bk+1=argminB{L(Xk+1,Ek+1,Ak+1,B,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)},T1k+1=T1k+ξ1k(Xk+1+Ek+1?Y),T2k+1=fold(n?k)((T2k)(n?k)+ξ2k(?1X(n?k)k+1?A(n?k)k+1)),T3k+1=fold(n?k)((T3k)(n?k)+ξ3k(?2X(n?k)k+1?B(n?k)k+1)),ξ1k+1=ρξ1k,ξ2k+1=ρξ2k,ξ3k+1=ρξ3k，其中ρ?>?1，k为迭代次数，n?k=argmax(min)1≤n≤N{∥X(n)k∥*}。接下来，我们给出每个子问题的闭式解。(1) 求解X子问题。将(22)中所有涉及X的项提取出来，可得X子问题的形式为：(24)Xk+1=argminX{ξ1k2∥X?Y+Ek+T1kξ1k∥F2+ξ2k2∥?1X(n?k)?A(n?k)k+(T2k)(n?k)ξ2k∥F2+ξ3k2∥?2X(n?k)?B(n?k)k+(T3k)(n?k)ξ3k∥F2}。该问题可通过以下方程求解：(25)ξ1kX(n?k)+ξ2k?1??1X(n?k)+ξ3k?2??2X(n?k)=Gk，其中Gk=ξ1kY(n?k)?ξ1kE(n?k)k?(T1k)(n?k)+ξ2k?1?A(n?k)k??1?(T2k)(n?k)+ξ3k?2?B(n?k)k??2?(T3k)(n?k)。此处?1?(·)和?2?(·)分别表示?1(?·?)和?2(?·?)的转置运算。方程(25)可通过快速傅里叶变换求解，其闭式解为：(26)X(n?k)k+1=F?1(F(Gk)ξ1k1+F(ξ2k?1??1+ξ3k?2??2))，其中F(·)和F?1(·)分别表示傅里叶变换及其逆变换，1为全1矩阵。最终，Xk+1可通过(27)Xk+1=fold(n?k)(X(n?k)k+1)得到。(2) 求解E子问题。固定其他变量，E子问题可表示为：(28)Ek+1=argminE{λψ(E)+ξ1k2∥E?Y+Xk+1+T1kξ1k∥F2}。该问题的闭式解为：(29)(Ek+1)i1i2?iN=proxλξ1kg((Y?Xk+1?T1kξ1k)i1i2?iN)，其中i?∈?[In]，n?∈?[N]。 (3) 求解A子问题。Ak+1可通过(30)Ak+1=fold(n?k)(A(n?k)k+1)得到，而A(n?k)k+1则由下式确定：(31)A(n?k)k+1=argminA(n?k){?(A(n?k))+ξ2k2∥A(n?k)??1X(n?k)k+1?(T2k)(n?k)ξ2k∥F2}。根据Wen等人[45]中的定理1，问题(31)可通过以下方式求解：(32)A(n?k)k+1=U1k+1S^1k+1(V1k+1)?，其中U1k+1S^1k+1(V1k+1)?=?1X(n?k)k+1+(T2k)(n?k)ξ2k，且(S^1k+1)ii=prox1ξ2kg((S1k+1)ii)，其中i?∈?[rk]，rk:=min{In?k,Kn?k}。 (4) 求解B子问题。与Ak+1类似，Bk+1可通过(33)Bk+1=fold(n?k)(B(n?k)k+1)得到，B(n?k)k+1由下式确定：(34)B(n?k)k+1=argminB(n?k){?(B(n?k))+ξ3k2∥B(n?k)??2X(n?k)k+1?(T3k)(n?k)ξ3k∥F2}。同样依据Wen等人[45]中的定理1，该问题也可通过类似方式求解。本文中未赘述这些证明过程。□引理2 若?或ψ中的函数g选取为L1、Lp、log、SCAD或MCP函数，则序列{T1k}、{T2k}和{T3k}是有界的。证明：首先证明?和ψ中的函数g为L1函数时的情况。根据T1k+1的更新规则，有∥T1k+1∥F=∥T1k+ξ1k(Xk+1+Ek+1?Y)∥F=ξ1k∥Ek+1?(Y?Xk+1?T1kξ1k)∥F=ξ1k∥proxλξ1kg(Y?Xk+1?T1kξ1k)?(Y?Xk+1?T1kξ1k)∥F=ξ1k∑i1=1I1?∑iN=1IN(proxλξ1kg(Qi1?iNk+1)?Qi1?iNk+1)2，其中Qk+1:=Y?Xk+1?T1kξ1k。若g为L1函数，则对于任意i?∈?[In]，n?∈?[N]，有|proxλξ1kg(Qi1?iNk+1)?Qi1?iNk+1|=λξ1k当|Qi1?iNk+1|>λξ1k时，且|proxλξ1kg(Qi1?iNk+1)?Qi1?iNk+1|=|Qi1?iNk+1|≤λξ1k当|Qi1?iNk+1|≤λξ1k时。因此，∥T1k+1∥F≤ξ1k∑i1=1I1?∑iN=1IN(λξ1k)2=λI1?IN，即序列{T1k}是有界的。同理可证∥T2k+1∥F≤rk，因此{T2k}和{T3k}也有界。若?或ψ中的函数g为Lp、SCAD或MCP函数，类似分析也可证明{T1k}、{T2k}和{T3k}的有界性，具体证明可参考Gao等人[7]的论述，此处不再详述。□引理3 若ψ中的函数g选取为L1、Lp或log函数，则序列{Xk,Ek,Ak,Bk}是有界的。证明：对于k=0,1,…，根据(23)给出的迭代规则，有(37)L(Xk+1,Ek+1,Ak+1,Bk+1,T1k+1,T2k+1,T3k+1,ξ1k+1,ξ2k+1,ξ3k+1)≤L(Xk+1,Ek+1,Ak+1,Bk+1,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)。根据增广拉格朗日函数(22)，有(38)L(Xk+1,Ek+1,Ak+1,Bk+1,T1k+1,T2k+1,T3k+1,ξ1k+1,ξ2k+1,ξ3k+1)?L(Xk,Ek,Ak,Bk,T1k,T2k,T3k,ξ1k,ξ2k,ξ3k)=〈T1k+1?T1k,Xk+1+Ek+1?Y〉+〈(T2k+1)(n?k)?(T2k)(n?k),?1X(n?k)k+1?A(n?k)k+1〉+〈(T3k+1)(n?k)?(T3k)(n?k),?2X(n?k)k+1?B(n?k)k+1〉+ξ1k+1?ξ1k2∥Xk+1+Ek+1?Y∥F2+ξ2k+1?ξ2k2∥?1X(n?k)k+1?A(n?k)k+1∥F2+ξ3k+1?ξ3k2∥?2X(n?k)k+1?B(n?k)k+1∥F2=ξ1k+1+ξ1k2(ξ1k)2∥T1k+1?T1k∥F2+ξ2k+1+ξ2k2(ξ2k)2∥T2k+1?T2k∥F2+ξ3k+1+ξ3k2(ξ3k)2。由于序列{T1k,T2k,T3k}是有界的，因此{L(Xk+1,Ek+1,Ak+1,Bk+1,T1k+1,T2k+1,T3k+1,ξ1k+1,ξ2k+1,ξ3k+1)}也是有上界的。由ξ1k+1=ρξ1k且ρ?>?1，可得(41)∑s=0∞ξ1s+1+ξ1s2(ξ1s)2=1+ρ2ξ10∑s=0∞1ρs=(1+ρ)ρ2ξ10(ρ?1)。同理可得∑s=0∞ξ2s+1+ξ2s2(ξ2s)2=(1+ρ)ρ2ξ20(ρ?1)，以及∑s=0∞ξ3s+1+ξ3s2(ξ3s)2=(1+ρ)ρ2ξ30(ρ?1)。(42)式可化为：L(Xk+1,Ek+1,Ak+1,Bk+1,T1k+1,T2k+1,T3k+1,ξ1k+1,ξ2k+1,ξ3k+1)+∥T1k+1∥F22ξ1k+1+∥T2k+1∥F22ξ2k+1+∥T3k+1∥F22ξ3k+1=?(Ak+1)+?(Bk+1)+λψ(Ek+1)+ξ1k+12∥Xk+1?Y+Ek+1+T1k+1ξ1k+1∥F2+ξ2k+12∥?1X(n?k)k+1?A(n?k)k+1+(T2k+1)(n?k)ξ2k+1∥F2+ξ3k+12∥?2X(n?k)k+1?B(n?k)k+1+(T3k+1)(n?k)ξ3k+1∥F2。由于(42)式中各项均为非负值，若选取ψ函数为L1、Lp或对数函数，则{Ek+1}序列的有界性可由{ψ(Ek+1)}的有界性保证；再结合{T1k,T2k,T3k}的有界性，可得出{Xk,Ak,Bk}序列也有界，至此证明完成。□在引理2中，若选?函数为对数函数，则需假设引理3中{T2k,T3k}序列有界；若选ψ函数为对数函数，则需假设{T1k}序列有界。在引理3中，若ψ函数为SCAD或MCP函数，则需假设{Ek}序列有界，才能保证{Xk,Ak,Bk}序列有界。接下来给出定理中的主要结果。定理2：在引理3所述的条件下，若满足limk→∞ξ1k(Ek?Ek+1)=0、limk→∞ξ2k(Ak+1?Ak)=0以及limk→∞ξ3k(Bk+1?Bk)=0，那么序列{Xk,Ek,Ak,Bk,T1k,T2k,T3k}的任意聚点都是问题(21)的KKT点。证明：根据引理2和引理3，可知该序列有界。依据波尔查诺-魏尔斯特拉斯定理，存在至少一个聚点，记为{X*,E*,A*,B*,T1*,T2*,T3*}。在不失一般性的前提下，假设该序列收敛于上述聚点。由于n?k的更新规则，存在kˉ使得当k≥kˉ时，n?k=n?kˉ（记n?*:=n?kˉ）。又因为{T1k,T2k,T3k}有界，且对于任意i∈[3]，都有ξik+1=ρξik，其中ρ>1，因此可得出limk→∞(Xk+1+Ek+1?Y)=limk→∞1ξ1k(T1k+1?T1k)=0，limk→∞(?1X(n?k)k+1?A(n?k)k+1)=limk→∞1ξ2k((T2k+1)(n?*)?(T2k)(n?*))=0，limk→∞(?2X(n?k)k+1?B(n?k)k+1)=limk→∞1ξ3k((T3k+1)(n?*)?(T3k)(n?*))=0，于是有(43)式成立，即Y=X*+E*，?1X(n?*)*=A(n?*)*，?2X(n?*)*=B(n?*)*。此外，根据问题(24)、(28)、(31)和(34)的一阶最优条件，可得0=ξ1k(X(n?k)k+1?Y(n?k)+E(n?k)k)+(T1k)(n?k)+?1?(ξ2k(?1X(n?k)k+1?A(n?k)k)+(T2k)(n?k))+?2?(ξ3k(?2X(n?k)k+1?B(n?k)k)+(T3k)(n?k))=(T1k+1)(n?k)?ξ1k(E(n?k)k+1?E(n?k)k)+?1?((T2k+1)(n?k)+ξ2k(A(n?k)k+1?A(n?k)k))+?2?((T3k+1)(n?k)+ξ3k(B(n?k)k+1?B(n?k)k))。进一步有0∈λ?Ek+1ψ(E)+ξ1k(Xk+1?Y+Ek+1)+T1k=λ?Ek+1ψ(E)+T1k+1，0∈?A(n?k)k+1?(A(n?k))+ξ2k(A(n?k)k+1??1X(n?k)k+1)?(T2k)(n?k)=?A(n?k)k+1?(A(n?k))?(T2k+1)(n?k)，0∈?B(n?k)k+1?(B(n?k))+ξ3k(B(n?k)k+1??2X(n?k)k+1)?(T3k)(n?k)=?B(n?k)k+1?(B(n?k))?(T3k+1)(n?k)。将k令趋于无穷大，结合给定条件，可得(45)式：T1*+?1?(T2*)(n?*)+?2?(T3*)(n?*)=0，0∈λ?E*ψ(E)+T1*，0∈?A(n?*)*?(A(n?*))?(T2*)(n?*)，0∈?B(n?*)*?(B(n?*))?(T3*)(n?*)。根据(43)和(45)式的结果，可最终得出{X*,E*,A*,B*,T1*,T2*,T3*}是问题(21)的KKT点。5. 实验部分本节通过实验验证所提方法的有效性。首先介绍整体实验设置，随后分别阐述详细实验步骤与结果。我们将所提出的TNLS方法与多种先进的张量去噪方法进行比较，这些方法包括：张量鲁棒主成分分析(TRPCA)[16]、改进型TRPCA(ETRPCA)[15]、三维相关总变分正则化方法(3DCTV)[9]、相关总变分方法(CTV)[34]、非凸双层TRPCA方法(NTRPCA)[7]以及典型系数相关总变分方法(RCCTV)[37]。为评估不同方法的性能，采用了三种常用指标：峰值信噪比(PSNR)和结构相似度(SSIM)值越高，相对误差(RE)值越低，表示性能越好。对于原始矩阵A∈RI1×I2和恢复后的矩阵A*∈RI1×I2，PSNR和SSIM[46]的定义分别为PSNR:=10log10(n1n22552∥A?A*∥F2)，SSIM:=(2μAμA*+ζ1)(2σAA*+ζ2)(μA2+μA*2+ζ1)(σA2+σA*2+ζ2)，其中μA和μA*分别是A和A*的均值，σAA*是A和A*的协方差，σA2和σA*2分别是A和A*的方差，ζ1和ζ2是为避免分母过小而引入的常数。PSNR和SSIM的值是通过计算张量所有切片的平均值得到的。对于原始张量A和恢复后的张量A*，相对误差RE的定义为RE:=∥A*?A∥F∥A∥F。算法1在满足??≤?0.001或迭代次数达到200时终止，其中?定义为max{∥Xk+1?Xk∥∞,∥Ek+1?Ek∥∞,∥Xk+1+Ek+1?Y∥∞}。实验中所使用的样本张量是通过在原始张量中加入预设噪声比δ的椒盐噪声生成的，噪声在张量范围内随机分布。所有实验均在配备13代Intel(R) Core(TM) i5-13500处理器（2.50?GHz）和16GB内存的台式机上使用Matlab R2025a进行。实验数据集包括彩色图像、灰度人脸图像以及多光谱图像(MSI)，具体如下：•伯克利分割数据集1[47]：从该数据集中随机选取30幅彩色图像，其尺寸为321?×?481?×?3或481?×?321?×?3。•耶鲁人脸数据库B2[48]：从该数据库中选取3个个体，每个个体有50幅图像，这些图像集分别记为人脸1、人脸2和人脸3，每组图像都被转换为尺寸为96?×?84?×?50的三阶张量。•MSI数据集：使用了两幅尺寸为128?×?128?×?31的MSI图像，分别记为MSI 1和MSI 2。算法1中张量展开的模式选择策略参考了Gao等人的方法[7]：对于彩色图像和灰度人脸图像，nk由nk:=argmax1≤n≤3{∥L(n)k∥*}确定；对于MSI图像，模式选择标准则改为nk:=argmin1≤n≤3{∥L(n)k∥*}。5.1 实验设置本小节将详细介绍实验设置，包括?和ψ函数的形式以及参数设置，并对收敛性进行数值验证。以尺寸为321?×?481?×?3的彩色图像作为测试示例。?和ψ函数的选择：在第3.1节中已经列出了满足假设1的一系列函数。我们会在固定噪声水平参数δ为0.1的情况下，评估不同函数配置下TNLS的性能。对于?和ψ，共考虑两种函数组合方式：第一种是让?和ψ使用相同的非凸函数g，对应的实验案例分别记为L1、log、SCAD、MCP和Lp。具体而言，当?和ψ都采用L1函数时，TNLS模型会退化为3节中提到的凸松弛模型。第二种则是为?和ψ使用不同的函数（?用g1，ψ用g2），共有六种具体组合：(1) gL1和gLp，记为L1?Lp；(2) gLp和gL1，记为Lp?L1；(3) glog和gMCP，记为log-MCP；(4) gLp和gMCP，记为Lp-MCP；(5) gMCP和gSCAD，记为MCP-SCAD；(6) gSCAD和glog，记为SCAD-log。图2展示了在上述不同配置下TNLS的恢复性能。从结果中可以得出两个关键结论：首先，当?和ψ都采用gL1函数时，其性能最差，这证明了非凸函数相比凸函数能更有效地捕捉张量的内在结构特征；其次，log、Lp、L1?Lp、log-MCP和Lp-MCP这几组的恢复性能优于其他组合。在这五组中，Lp和Lp-MCP不仅具有更高的PSNR和SSIM值，相对误差也更低。此外，Lp方案涉及的参数更少，这提升了其实用性，同时也降低了参数调整的复杂度。因此，我们在后续实验中选择Lp方案用于TNLS，并与其他方法进行性能比较。下载：下载高分辨率图像（354KB）下载：下载完整尺寸图像图2. 不同函数g下TNLS的恢复性能对比。参数设置：所提出的TNLS模型包含一些需要调整的参数，包括正则化参数λ、拉格朗日惩罚参数ξ1、ξ2、ξ3，以及Lp函数中的参数p。在实验中，我们将ξ设置为ξ:=ξ1=ξ2=ξ3。为了确定最优参数范围，我们对δ为0.1的彩色图像进行了敏感性分析，即固定另外两个参数，仅调整其中一个参数。图3(a)展示了λ在0.02到0.08之间以0.005为间隔变化时TNLS的PSNR和SSIM值。可以看出，当λ约为0.065时，模型的性能最为稳定且优异。基于这一结果，对于彩色图像，我们将λ的范围设定为[0.01, 0.07]。对于灰度人脸图像和MSI，λ的范围则调整为更小的[0.001, 0.02]。对于参数ξ和p，我们测试了多种组合，ξ的取值范围为ξ∈{0.0001,0.001,0.005,0.01,0.05,0.1,0.15,…,0.6}，p的取值范围为0.36到0.6，间隔为0.02。如图3(b)和(c)所示，当ξ小于0.025且p大于0.4时，模型在大多数情况下都能保持良好的性能。在实验中，ξ的选择范围为{0.001, 0.005, 0.01}，p的范围限制在[0.5, 0.7]。为保证比较的公平性，基准方法TRPCA、ETRPCA、3DCTV、CTV和NTRPCA的参数也按照其对应论文中的建议进行设置。下载：下载高分辨率图像（370KB）下载：下载完整尺寸图像图3. 不同λ、ξ和p值下TNLS的PSNR和SSIM值。算法1的数值收敛性：我们使用δ为0.1、0.2、0.3的彩色图像来评估算法1的数值收敛性。图4展示了PSNR、SSIM和RE随迭代次数的变化情况。可以明显看出，这三项指标的曲线在大约40次迭代后就能迅速达到理想状态（高PSNR/SSIM值和低RE值），并保持稳定。这一现象充分证明了算法4的数值稳定性及其高效的收敛速度。下载：下载高分辨率图像（407KB）下载：下载完整尺寸图像图4. PSNR、SSIM和RE随迭代次数的变化情况。5.2 结果比较本小节将展示TRPCA、ETRPCA、3DCTV、CTV、NTRPCA、RCCTV和TNLS在δ为0.1、0.2、0.3的条件下，对彩色图像、灰度人脸图像以及MSI的恢复效果。同时，我们还会对实验结果进行分析，探讨TNLS相较于其他方法的性能优势。实验结果概述：具体而言，图5展示了30幅不同噪声水平下的彩色图像上TNLS及各类基准方法的性能数值。图6、图7、图8分别展示了在δ为0.1、0.2和0.3时，各方法处理过的原始图像、含噪声图像以及恢复后图像的示例。图中还标注了部分区域的放大图，通过对比这些放大区域与原始图像对应区域的差异，直观显示了不同方法在细节保留和噪声去除方面的优劣。表2总结了图6、图7、图8中所展示的六类图像的定量比较结果，同时给出了30幅彩色图像的PSNR、SSIM和RE的平均值。表3和表4分别展示了各类方法在灰度人脸图像和MSI上的数值结果。从视觉效果来看，图9展示了在δ为0.2的灰度人脸数据集中，不同方法处理10幅图像的恢复效果。图10的布局与图6、图7、图8类似，用于展示MSI图像的处理效果。实验结果分析：从表2、表3、表4的定量结果可以看出，TNLS在大多数数据类型和噪声水平下均表现出优于基准方法的优势。例如，在彩色图像上，当δ为0.2时，TNLS的平均PSNR值为38.45?dB，这一数值分别比TRPCA、ETRPCA、3DCTV、CTV、NTRPCA和RCCTV高9.26?dB、7.56?dB、7.08?dB、3.24?dB、3.30?dB和6.05?dB。TNLS的SSIM值可达0.9800，相比其他方法高出0.0117至0.0712，而相对误差则降至0.0271，比基准方法低0.0112至0.0511。对于灰度人脸图像和医学影像，TNLS在恢复质量方面仍具有显著优势，在大多数情况下，其RE值更低，而PSNR/SSIM值更高，优于所有其他对比方法。图6、图7、图8、图9和图10中的视觉对比结果（涵盖彩色图像、灰度人脸图像和医学影像）均表明，TNLS的恢复图像在细节和整体外观上与原始图像更为一致，视觉质量明显优于其他方法。恢复图像与原始图像之间的差异图以更直观的方式凸显了这一优势：TNLS的差异图中白色区域更多，这直接表明其与原始图像的偏差极小；而其他方法的差异图中则存在更为明显的彩色或黑色斑块，这是由于噪声去除不彻底或图像细节丢失所致。总之，TNLS在性能和细节保留方面均优于基准方法，从而实现了显著提升的视觉恢复质量。

此外，我们还有以下观察结果。首先，非凸方法NTRPCA和ETRPCA的表现始终优于其凸形式TRPCA。其次，在大多数情况下，那些利用梯度矩阵的低秩特性来同时考虑低秩性和平滑性的方法（3DCTV、CTV和RCCTV），其性能优于直接限制低秩性的方法（TRPCA）。第三，非凸方法NTRPCA在大多数场景下表现优于CTV，但在某些情况下可能失效。这些观察结果验证了非凸方法以及结合低秩性与平滑性思路的方法的优势。我们提出的TNLS融合了这两种方法的优势，因而能够实现最佳性能。

为更直观地展示TNLS的性能优势，表5列出了TNLS相较于其他方法的相对改进率。表5中的所有数值都是通过计算相同类型数据集在不同噪声水平δ下的平均结果得出的，相对改进率定义为TNLS与各对比方法之间的绝对性能差异除以相应方法的指标值。从表5可以看出，TNLS在大多数情况下都能持续优于所有对比方法。具体而言，与非凸方法ETRPCA和NTRPCA相比，TNLS因能在非凸框架内同时结合低秩性和平滑性先验，而取得了更好的效果；而ETRPCA和NTRPCA则未能引入平滑性约束。与基于梯度矩阵的方法（3DCTV、CTV和RCCTV）相比，TNLS主要得益于对秩函数和L0范数的非凸替代函数，从而提升了性能。总之，TNLS不仅继承了非凸方法在低秩近似方面的优势，还借助基于梯度的策略来保留图像的精细细节，因此在数值指标和视觉质量方面均优于所有对比方法。

表5. TNLS与其他方法的性能比较及相对改进率
| 方法 | 彩色图像 | 灰度人脸图像 | 医学影像 |
|--------------|----------|--------------|----------|
| PSNR ↑ | 26.67% | 6.91% | 59.17% |
| SSIM ↑ | 27.59% | 11.15% | 27.59% |
| RE ↓ | 60.34% | 60.34% | 11.15% |
| PSNR ↑ | 24.81% | 54.00% | 60.59% |
| SSIM ↑ | 1.00% | 17.86% | 60.59% |
| RE ↓ | 69.01% | 0.35% | 30.41% |

| 方法 | 彩色图像 | 灰度人脸图像 | 医学影像 |
|--------------|----------|--------------|----------|
| CTV | 9.06% | 1.10% | 25.67% |
| NTRPCA | 9.00% | 1.52% | 29.63% |
| RCCTV | 17.62% | 4.21% | 42.56% |

6. 结论与未来工作

在本文中，我们提出了一种名为TNLS的两层非凸张量模型，该模型结合了全局低秩性和局部平滑性先验，用于张量数据去噪。TNLS采用通用的非凸替代函数，同时捕捉张量数据的低秩性和稀疏性特征。在理论分析方面，我们为TNLS的最坏情况误差界建立了通用理论框架，从而保证了所提方法的理论有效性和稳定性。在算法实现方面，我们设计了一种具有闭式更新形式的高效ADMM算法，并证明其能够收敛到KKT点。在多种数据集上的大量实验表明，TNLS的性能优于现有的几种先进方法。

本文提出的方法还可以进一步扩展。例如，我们可以考虑整合张量数据的额外先验信息及潜在结构特征，如局部自相似性和群体稀疏性，以此进一步提升该方法的性能，扩大其应用范围，或处理更复杂的现实世界数据。此外，TNLS将低秩性与平滑性结合在非凸替代函数中的思路，也可应用于其他张量模型，以解决诸如图像恢复和物体检测等更实际的任务。

CRediT作者贡献说明
高凯欣：写作——审稿与编辑，写作——初稿撰写，可视化，验证，软件，方法论，研究，概念设计。孟杰：写作——审稿与编辑，验证，监督。张翔：写作——审稿与编辑，可视化，验证。