论文引用格式：Long Z， Liu Y P， Gou Y X， Zeng S X， Liu J N， Wen F and Zhu C. 2023. Linear tensor subapace model based hyperspectral image denoising. Journal of Image and Graphics， 28（08）：2505-2521（引用格式:龙珍， 刘翼鹏， 苟艺馨， 曾思行， 刘佳妮， 文飞， 朱策. 2023. 高光谱图像的线性张量子空间模型及降噪应用. 中国图象图形学报， 28（08）：2505-2521）［0　引 言高光谱成像技术将成像技术与光谱技术相结合，在给定光谱范围内，同时采集目标的二维空间和一维光谱信息，比如光谱范围可以从400 nm到 2 500 nm。凭借高光谱图像（hyperspectral image，HSI）丰富的光谱信息（Manolakis和Shaw，2002；Zhang等，2019），高光谱图像在许多领域获得关注，如遥感（Bioucas-Dias等，2013）、食品安全（Qin等，2013）和医学诊断（Lu和Fei，2014；Fei，2019；吕蒙 等，2021）。然而，由于热电子学、暗电流以及成像过程中光计数的随机误差，高光谱图像在采集过程中不可避免地会受到严重的噪声影响，这必然会降低图像质量并阻碍其后续应用，例如目标检测（Manolakis等，2003）、高光谱图像分类（Sun等，2020；Liu等，2022；Fu等，2021；程嵩阳 等，2021）和材料识别（Zhang等，2008；张启忠 等，2021）。因此，高光谱图像中的噪声去除已成为一个重要的研究领域（Zhang等，2020；Duan等，2020；Maffei等，2020；Nguyen等，2021；Shi等，2021；Xue等，2019；巫勇 等，2021；刘盛 等，2021）。最新的高光谱图像去噪算法可以分为两类，一类是基于数据驱动的方式；另一类是基于模型驱动的方式。基于深度学习的方法作为数据驱动方法的代表工作，在高光谱图像去噪方面表现出不错的性能。例如，Xie和Li（2017）提出了一种具有可训练非线性函数的深度卷积神经网络（convolutional neural network， CNN） 高光谱图像降噪算法，但该算法仅能对单个波段进行去噪。为了同时考虑空间光谱信息，Yuan等人（2019）提出了组合空间光谱CNN，但该模型需要针对不同的噪声水平进行重新训练。Maffei等人（2020）提出了一种基于高光谱图像中空间光谱相关性的去噪CNN方法。但是，该方法要求所有波段的噪声水平相当。总体来说，现有的基于深度学习的方法相对简单，不能很好地模拟复杂的噪声条件。从基于模型的角度来看，低秩近似的方法在没有训练的情况下表现出色。目前，基于低秩的去噪方法可以分为基于矩阵的去噪方法和基于张量的去噪方法。基于矩阵的去噪方法要么将三维张量展开为矩阵，要么使用传统的图像去噪算法独立处理每个波段。例如，Zhang等人（2014）将高光谱图像分成几个片段，将这些片段重新排列成一个二维矩阵，并使用低秩矩阵恢复（low rank matrix recovery， LRMR）方法恢复每个片段。该方法在恢复混合噪声的高光谱图像方面非常成功，产生了一系列基于LRMR（He等，2022；Zhang等，2020；Wang等，2020）的去噪模型。但是，由于高光谱图像的空间光谱信息被部分破坏，这些二维去噪算法难以达到最佳效果。针对这个问题，Fan等人（2017）提出了低秩张量恢复（low rank tensor recovery，LRTR）高光谱图像去噪算法，同时使用高光谱图像的光谱和空间信息，并且比LRMR取得了更好的降噪效果。随后，一些利用不同张量分解的降噪工作也相应开展。例如，Wang等人（2018）使用低秩Tucker近似来模拟无噪声污染的高光谱图像，使用𝓁1-范数去除稀疏和死线噪声，通过全变分（total variation，TV）最小化（Liu等，2019，2020）对高光谱图像的分段平滑结构进行建模。之后，Fan等人（2018）考虑了一种空间光谱总变分和张量奇异值分解（tensor-singular value decomposition，t-SVD）方法来去除高光谱图像中的混合噪声。除了TV约束，Xue等人（2019）和Gong等人（2020）考虑高光谱图像去噪的非局部相似性，其中高光谱图像的相似块被聚集成一个组并一起处理。此外，Xiong等人（2018）考虑了高光谱图像的非负性，并提出非负低秩矩阵向量张量分解的高光谱图像降噪算法。然而，在上述LRTR降噪方法中，仅将高光谱图像视为三阶低秩张量。事实上，高光谱图像除了全局低秩张量的特点，还符合线性子空间模型。例如，给定高光谱图像𝒳∈RI×J×K，其中I和J是空间维度大小，K是波段数。每个谱向量xi,j∈RK,i=1,⋯,I;j=1,⋯,J均位于低维的线性子空间中（Zhuang和Bioucas-Dias，2018；Zheng等，2020；Cao等，2019）。图 1显示了一个低维子空间表示谱向量的示例，其中彩色向量和图像矩阵分别表示子空间的正交基及其对应的系数，低维子空间由HySime （hyperspectral signal identification by minimum error）算法（Bioucas-Dias和Nascimento，2008）获得。此外，可以观察到每个子空间对应的系数具有低秩结构。10.11834/jig.220306.F001图1高光谱图像的线性子空间解释模型Fig.1An illustration of the linear subspace for HSI受上述启发，本文提出了一种具有正交向量的结构低秩矩阵向量张量分解（matrix-vector tensor factorization，MVTF），并应用于高光谱图像去噪。如图 2所示，新获得的张量分解可以有效地利用高光谱图像的全局空间光谱结构。此外，F-范数和𝓁1-范数分别用于处理优化模型中的高斯噪声和脉冲噪声，总变分项用于模拟高光谱图像的空间局部平滑属性。优化模型在交替方向乘子法（alternating direction method of multipliers，ADMM）框架下，通过将优化问题划分为多个子问题进行求解。基于高光谱图像降噪实验结果表明，所提出的方法在平均峰值信噪比（mean peak signal to noise ratio，MPSNR）、平均结构相似性指数（mean structural similarity index，MSSIM）（Yuan等，2015）及无量纲全局相对综合误差（relative dimensionless global error synthesis，ERGAS）（Wald，2000）等方面优于现有算法。10.11834/jig.220306.F002图2矩阵向量形式的张量分解Fig.2Matrix vector tensor factorization本文的主要贡献可归纳如下：1）提出一种具有正交向量结构的低秩MVTF算法来建模高光谱图像，其中正交向量对应谱向量的基，矩阵因子的核范数最小化意味着每个系数矩阵具有低秩结构。与标准MVTF上的方法相比，所提方法的运算速度要快得多，并且无需提前给定MVTF空间系数矩阵的秩。2）各向异性总变分（anisotropic total variation，ATV）项用于低秩张量恢复的优化模型中，它可以建模高光谱图像空间域的局部分段平滑属性，提高降噪性能。3）在模拟和真实数据上的实验结果表明了所提方法在MPSNR、MSSIM和ERGAS等衡量指标上的优越性。1.1　符号说明在本文中，标量、向量、矩阵和张量分别用小写字母x、粗体小写字母x、粗体大写字母X和花体字母X表示。对于三阶张量X∈RI1×I2×I3，其第{i1,i2,i3}项记为X(i1,i2,i3)及其正向切片记为X(:,:,i3),i3=1,⋯,I3。X*表示X的共轭。矩阵X的𝓁1-范数为X1=∑i1=1I1∑i2=1I2xi1,i2，矩阵X的F-范数定义为XF=∑i1=1I1∑i2=1I2xi1,i221/2，其中，xi1,i2表示X的第(i1,i2)项。秩R矩阵X的核范数定义为X*=∑r=1Rσr(X), 其中σr(X)是X的第r个最大奇异值。矩阵X∈RI×I的迹定义为tr(X)=∑i=1IX(i,i)，N阶张量X∈RI1×I2×⋯×IN与矩阵U(n)∈RJn×In的模-n积定义为(X×nU(n))i1,⋯,in-1,in,in+1,⋯,iN=∑in=1Inxi1,i2,⋯,iNujn,in。N阶张量X∈RI1×I2×⋯×IN的模-n展开为X(n)，其中，X(n)(in,j)=X(i1,⋯,in,⋯,iN),j=i1,⋯,in-1,in+1,⋯,iN¯ 。对于一个三阶的高光谱图像 X∈RI×J×K，空间域是局部分段平滑的，因此高光谱图像的各向异性总变分记为 XTV=𝒟(X)1=𝒟h(X)1+𝒟v(X)1,其中，𝒟(⋅)=𝒟h(⋅);𝒟v(⋅)是一个二维差分算子，𝒟h和𝒟v是沿着两种不同方式的高光谱数据的一阶有限差分算子，定义为𝒟h(X)=X(i+1,j,:)-X(i,j,:)； 𝒟v(X)=X(i,j+1,:)-X(i,j,:)，其中i=1,⋯,I-1； j=1,⋯,J-1。1.2　线性子空间模型给定一个高光谱图像X∈RI×J×K，谱向量xi,j∈RK,i=1,⋯,I;j=1,⋯,J位于低维子空间中，数学式可以写成X=CZ（Zhuang和Bioucas-Dias，2018），其中X∈RK×IJ是高光谱图像的矩阵表示，C=[c1,⋯,cR]∈RK×R是R维子空间的基，矩阵Z∈RR×IJ表示X相对于C的系数。在大多数情况下，CTC=IR表示子空间的基是不相关的，其中，IR代表单位矩阵，C可以通过 HySime算法或奇异值分解（singular value decomposition，SVD）获得，Z的行可以进一步转化为矩阵。1.3　矩阵向量张量分解MVTF将一个张量X∈RI×J×K分解为几个张量因子之和，每个张量因子可以写成一个矩阵Gr∈RI×J和一个向量cr∈RK的外积，表示为X≈∑r=1RGr∘cr=G×3C （1）式中，C∈RK×R，G∈RI×J×R，Gr是G的第r个正向切片， ∘ 表示外积。2　优化模型如上所述，高光谱图像是一个三阶张量，具有两个空间维度和一个光谱维度。在数据采集过程中，不可避免地会受到噪声的污染。考虑到观测图像𝒴可能受到高斯噪声和异常值的污染，建立测量模型如下Y=X+S+N （2）式中，Y是观察到的图像，X是干净的图像，S是稀疏噪声或异常值，N是高斯噪音。高光谱图像去噪的目的是从被噪声污染的观测数据中估计出干净的高光谱图像。因此，本文在考虑高光谱图像子空间先验信息的基础上，提出了一种新的平滑且稳健的低秩部分正交MVTF的降噪算法。首先，基于高光谱图像的线性子空间假设，提出部分正交低秩MVTF近似模型，它可以很好地表示干净的高光谱图像，模型为minGr,r=1,⋯R,C∑r=1RGr* s.t.  𝒳=𝒢×3C,CTC=IR （3）式中，X表示干净的高光谱图像。cr表示R维子空间的第r个正交基，CTC=IR表示子空间的基是不相关的。Gr表示相应的第r个特征图像，往往是低秩的，如图 1所示。由于函数rank（Gr）是非凸函数，常常使用它的凸近似函数Gr*来刻画这种低秩结构。接下来，使用𝓁1-范数将稀疏异常值与观察值分开。此外，为了提高恢复性能，使用总变分项模拟高光谱图像数据的局部平滑结构。最后，得到用于高光谱图像去噪的平滑且鲁棒的低秩张量恢复（smoothing and robust low rank tensor recovery， SRLRTR）的优化模型，具体为minX,S,N,G,𝒞λTV𝒟(X)1+λSS1+λNNF2+λ𝒢∑r=1RGr*s.t.  Y=X+S+N,X=G×3C,CTC=IR （4）为了解决这个问题，需要引入两个辅助变量Z和L，优化模型（4）可以改写为minX,S,N,G,𝒞,Z,ℒλTVL1+λSS1+λNNF2+λG∑r=1RGr* s.t.  Y=X+S+N, Z=X, L=𝒟(Z)  X=G×3C, CTC=IR （5）优化模型（5）在CTC=IR约束下的增广拉格朗日函数表示为ℒ(𝒳,S,N,G,C,Z,L,Λ1,Λ2,Λ3,Λ4)=λTVL1+λSS1+λNNF2+λG∑r=1RGr*+Λ1,Y-X-S-N+Λ2,Z-X+Λ3,L-𝒟(Z)+Λ4,X-G×3C+β12Y-X-S-NF2+β22Z-XF2+β32L-𝒟(Z)F2+β42X-G×3CF2 （6）式中， Λ1， Λ2， Λ3和Λ4表示对偶变量， β1， β2， β3和β4表示惩罚因子。该优化模型主要通过ADMM方法求解，即固定其他变量，依次更新一个变量。2.1　更新因子G固定除G以外的其他变量，此时优化模型（6）中只有变量G，式（6）可以简化为minGr,r=1,⋯,RλG∑r=1RGr*+Λ4,X-G×3C+β42X-G×3CF2 （7）式中，G(:,:,r)=Gr,r=1,⋯,R。合并式（7）中的后两项，可得Λ4,X-G×3C+β42X-G×3CF2=tr(Λ4T(X-G×3C))+β42tr((X-G×3C)T(X-G×3C))=trβ42GTG-2GTX+Λ4β×3CT+const=β42G-X+Λ4β×3CTF2+const （8）式中，const表示常数值。优化问题（7）可以等价为minGr,r=1,⋯,RλG∑r=1RGr*+β42G-ℳF2 （9）式中， M=(X+Λ4 / β4)×3CT， 因此针对每个Gr可以简化为minGr,r=1,⋯,RλG∑r=1RGr*+β42∑r=1RGr-M(:, :, r)F2 （10）式中，Gr=G(:,:,r)。上述问题的解可以通过奇异值算子直接获得。G的详细求解方法总结在算法1中，其中SVD表示奇异值算子， x的软阈值算子定义为 Thrτ(x)=sgn(x)max(|x|-τ,0)。算法1： 更新因子G输入：干净的图像X∈RI×J×K，拉格朗日乘子Λ4，惩罚因子β4，参数λ𝒮和λ𝒢，子空间维数R以及光谱子空间的基C。输出：因子G的估计值G^。1）更新τ=λGβ4。2）M=X+Λ4β4×3CT3）从r = 1开始循环：4） [U,Σ,V]=SVD(M(:,:,r))5） Σ^=Thrτ(Σ)6） Gr=UΣ^VT7） G^(:,:,r)=Gr8） 当r = R时，结束循环。2.2　更新因子C在模型（6）中，固定除C外的其他变量，可以得到优化模型为minCΛ4,X-G×3C+β42X-G×3CF2  s.t.  CTC=IR （11）化简上式，该问题（11）等价于maxCtr(G(3)(Λ4(3)T+β4X(3)T)C)s.t.  CTC=IR （12）令M=G(3)(Λ4(3)T+β4X(3)T)。上述为著名的正交普鲁克问题，通过在M上执行奇异值分解M=UΣVT，C的解为C=VUT。2.3　更新因子X固定除X以外的变量，优化模型（6）可以重写为minXΛ1,Y-X-S-N+Λ2,Z-X+Λ4,X-G×3C+β12Y-X-S-NF2+β22Z-XF2+β42X-G×3CF2 （13）通过目标函数（13）对X求导，并令求导结果为0，可以得到X的解为X=β1(Y-S-N)+Λ1+β2Z+Λ2+β4(𝒢×3C)-Λ4β1+β2+β4 （14）2.4　更新因子Z类似地，Z的求解优化模型（6）可以重写为minZΛ2, Z-X+Λ3, L-𝒟(Z)+β22Z-XF2+β32L-𝒟(Z)F2 （15）通过对Z求导，并令结果等于0，可以获得Λ2-𝒟*(Λ3)+β2(Z-X)-𝒟*(β3(L-𝒟(Z)))=0 （16）该问题的解化简为(β2I+β3𝒟*𝒟)Z=β2X-Λ2+𝒟*(β3L+Λ3)   （17）式中，算子𝒟*𝒟具有块循环结构，可以通过傅里叶变换快速计算，𝒵的快速计算可以写成Z=ifftnfftn(M)β21+β3(fftn(𝒟*𝒟)) （18）式中，M=β2X-Λ2+𝒟*(β3L+Λ3)， fftn和ifftn分别是3D快速傅里叶变换及其逆变换。2.5　更新因子L为更新L，固定优化函数（6）中的所有其余变量，子问题可以表述为minLλTVL1+Λ3,L-𝒟(Z)+β32L-𝒟(Z)F2 （19）类似于G的后两项合并，该问题可以转化为minLλTVL1+β32L-𝒟(Z)-Λ3β3F2 （20）L的解可以通过软阈值算子得到，具体为L=ThrλTVβ3𝒟(Z)-Λ3β3 （21）2.6　更新因子S通过固定优化函数（6）中除S以外的变量，可以获得子问题如下minSλSS1+Λ1,Y-X-S-N+β12Y-X-S-NF2 （22）与上面求解L类似，这个优化模型可以通过以下方式求解S=Thrλ𝒮β1Y-X-N+Λ1β1 （23）2.7　更新因子N固定方程（6）中的其他变量，更新N的优化模型可以表述为minNλNNF2+Λ1, Y-X-S-N+β12Y-X-S-NF2 （24）对N求导，并令其导数为0，可以获得N=β1(Y-X-S)+Λ1β1+2λN （25）2.8　拉格朗日乘子更新拉格朗日乘子更新如下Λ1=Λ1+β1(Y-X-S-N)Λ2=Λ2+β2(Z-X)Λ3=Λ3+β3(L-𝒟(Z))Λ4=Λ4+β4(X-G×3C) （26）算法 2总结了所提出的基于部分正交MVTF的平滑稳健低秩张量恢复（SRLRTR）的高光谱图像降噪方法。当两个连续张量X之间的相对误差小于阈值，达到收敛条件，即X(k)-X(k+1)F2 / X(k+1)F2≤ε，其中X(k)是在第k次迭代恢复的结果，ε是阈值。在本文中设置ε=10-5。算法2 ：部分正交MVTF的平滑稳健低秩张量恢复算法输入：需要去噪的图像Y∈RI×J×K，惩罚因子β1,β2,β3,β4，拉格朗日乘子λS,λG,λN,λTV以及阈值ε，最大迭代次数K。输出：干净的图像X。循环迭代：1）通过算法1更新G。2）利用式（10）更新C。3）利用式（12）更新X。4）利用式（15）更新Z。5）利用式（18）更新L。6）利用式（20）更新S。7）利用式（22）更新N。8）利用式（23）更新对偶变量。检查收敛条件X(k)-X(k+1)F2 / X(k+1)F2≤εk←k+1当满足收敛条件或kK时，结束循环。2.9　计算复杂度分析本文算法的主要复杂度来自G的更新。对于一幅高光谱图像X∈RI×J×K，最耗时的部分是计算M(:,:,r)∈RI×J,r=1,⋯,R的奇异值分解。假设I=J，计算复杂度为O(I3)。因此，总复杂度为O(KRI3)，其中K是迭代次数。3　实验结果实验主要使用两个模拟数据集和两个真实数据集来验证所提算法的有效性。实验前将高光谱图像中所有像素的反射率值归一化在0和1之间。3.1　数据集1）Washington DC Mall （WDC）数据集。该数据集由美国弗吉尼亚光谱信息技术应用中心许可提供。每个场景大小为1 208 × 307像素，包括400～2 400 nm光谱范围内的210个波段。去除大气不透明的波段，在实验中使用191个波段，并从每幅图像中选取256 × 256个像素，形成大小为256 × 256 × 191的高光谱图像。WDC的伪彩色图像（R = 29，G = 19，B = 9）如图 3（a）所示。10.11834/jig.220306.F003图3高光谱图像测试数据集Fig.3Testing HSI dataset（（a） WDC；（b） PaviaU；（c） EO-1 ；（d） Urban）2）Pavia University （PaviaU）数据集。该数据集通过ROSIS（reflective optics system imaging spectrometer）传感器获取了意大利北部帕维亚区域。PaviaU的空间尺寸为610 × 340像素，有103个通道，范围为430～860 nm。在本实验中，考虑300 ×300像素，形成大小为300 × 300 × 103的高光谱图像。PaviaU的假彩色图像（R = 29，G = 19，B = 9）如图 3（b）所示。3）EO-1 Hyperion（EO-1）高光谱数据集。该数据集涵盖了美国印第安纳州的一个农业区，包含1 000 × 1 000的空间像素和350～2 600 nm的242个波段。去除水吸收的波段，保留剩下的166个波段，本实验从原始场景中选取200 × 200空间像素，形成大小为200 × 200 × 166的观测高光谱图像，该数据主要被高斯噪声、条纹和死线噪声污染。EO-1的假彩色图像（R = 2，G = 132，B = 136）如图3（c）所示。4）HYDICE Urban Dataset（Urban）数据集。该数据集从高光谱数字影像采集实验（Mitchell，1995）中获取，包含307 × 307的空间像素，210个波段，其中波长范围为400～2 500 nm。该数据集主要被条纹和死线噪声污染。Urban的假彩色图像（R = 2，G = 109，B = 207）如图 3（d）所示。3.2　对比方法和衡量指标为了与所提方法进行比较，选择了7种相关的先进算法。具体如下：1）NLRCPTD（nonlocal low-rank regularized CP tensor decomposition）（Xue等，2019），一种基于非局部低秩正则化CP分解的高光谱图像去噪方法。2）LRTDTV（TV-regularized low-rank tensor decomposition）（Wang等，2018），一种基于总变分和低秩Tucker分解的高光谱图像混合噪声去除方法。3）RSLRNTF（reweighed sparse low-rank nonnegative tensor factorization）（Xiong等，2018），一种基于矩阵向量张量分解的新方法，具有低秩和非负约束，用于高光谱图像去噪。4）SSTVLRTF（spatial-spectral TV regularized low-rank tensor factorization）（Fan等，2018），一种基于空间光谱变分和低秩t-SVD分解的方法，用于去除高光谱图像混合降噪。5）SNLRSF（subspace-based nonlocal low-rank and sparse factorization）（Cao等，2019），一种基于非局部低秩和稀疏约束的新型子空间方法，用于高光谱图像混合降噪。6）LRTFDFR（doubled-factor-regularized low-rank tensor factorization）（Zheng等，2020），一种基于双因子正则化低秩张量分解模型的子空间方法，用于高光谱图像混合噪声去除。7）LTDL（low-rank tensor dictionary learning）（Gong等，2020），一种基于低秩张量字典学习的高光谱图像噪声去除新方法。比较算法的参数根据其相应文献中的值设置。此外，为了定量评估性能，使用MPSNR、MSSIM、ERGAS和CPU time共4个指标来评估图像的去噪质量。具体为fMPSNR=1K∑k=1K fPSNRkfMSSIM=1K∑k=1K fSSIMk (27)fERGAS=1K∑k=1KX^(:, :, k)-X(:, :, k)F2μk2式中，K是光谱带的数量。μk是X(:,:,k)的平均值，X^(:,:,k)是恢复的数据。所有实验均使用MATLAB 2019b在具有2.4 GHz四核Intel Core i5处理器和16 GB RAM的计算机上进行。3.3　模拟数据实验结果本实验使用WDC和PaviaU高光谱图像作为模拟数据集。将高斯噪声、椒盐（脉冲）噪声、死线和条状噪声添加到模拟的真实数据集以验证所提方法。进行了4组实验，实验设置如下：实验1 ：高斯噪声+脉冲噪声，在这组实验中，不同波段的噪声强度是相同的。在每个波段中，添加相同的零均值，方差为0.2的高斯噪声和相同百分比为0.2的脉冲噪声。实验2：高斯噪声+死线。在这组实验中，高斯噪声以与实验1相同的方式添加到每个波段，其中噪声水平等于0.15。此外，在从波段41～100的光谱段中添加死线噪声，条纹的数量从3到10随机选择，宽度从1到3中随机生成。实验3：高斯噪声+脉冲噪声+死线。在实验2的基础上，每个波段的高斯噪声的方差设置为0.05，并添加百分比为0.1的脉冲噪声。实验4：高斯噪声+脉冲噪声+死线+条纹。本组实验在实验3的基础上，对于WDC和PaviaU数据，分别在10～190波段之间和61～100波段之间随机添加条纹噪声，条纹数量范围为20到40。对于本文方法，参数设置为λTV = 0.000 2，λ𝒮=0.02，λN=0.1， λG=0.1和R=5。在本文讨论部分会给出如何选择这些参数。表 1—表4显示了不同算法在不同实验设置下，在WDC和PaviaU数据集上，使用MPSNR、MSSIM、ERGAS和CPU time等衡量指标的降噪实验结果。可以观察到，在大多数情况下，本文算法在MPSNR、MSSIM和ERGA方面都优于其他方法。具体来说，本文算法在实验1中显示了去除混合高斯噪声和脉冲噪声的优势，意味着该算法可以有效处理高斯和脉冲噪声。同样，当实验3和实验4包含死线和条纹噪声时，本文方法在恢复性能方面优于现有的最先进方法，这可能意味着提出的方法在处理各种条件下的混合噪声方面具有优势。在包含高斯噪声和死线的实验2中，SNLRSF、NLRCPTD和LTDL性能更好，但计算复杂度在CPU时间方面更高。10.11834/jig.220306.T001表1实验1条件下，不同算法在 PaviaU 和 WDC 数据集上的降噪性能比较Table 1Quantitative comparison of different algorithms on PaviaU and WDC datasets denoising in experiment 1方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sRSLRNTF19.260.44424.75723.3619.200.59490.13631.34LRTDTV29.010.79142.19179.9627.980.78162.25327.80SSTVLRTF20.930.51325.18729.2721.440.61380.481.34×103NLRCPTD22.130.57387.031.8×10419.300.56281.503.98×104SNLRSF22.050.695.1×103794.3321.430.685.5×103601.76LTDL20.110.62449.199.1×10319.420.61480.453.8×104LRTFDFR28.410.75154.1948.4828.080.82165.0555.77SRLRTR29.170.75141.38245.4529.6420.85141.34271.30注：加粗字体表示各列最优结果。10.11834/jig.220306.T002表2实验2条件下，不同算法在 PaviaU 和 WDC 数据集上的降噪性能比较Table 2Quantitative comparison of different algorithms on PaviaU and WDC datasets denoising in experiment 2方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sRSLRNTF31.290.819111.53446.231.890.91104.13888.73LRTDTV31.120.851108.186132.5730.270.85124.26272.92SSTVLRTF30.700.746121.21657.6630.710.89107.481.05×103NLRCPTD33.040.88290.212.7×10431.940.80126.204.24×104SNLRSF35.760.9241.06×1032.1×10336.250.9661.05×103590.02LTDL34.070.9080.958.9×10333.860.943385.573.6×104LRTFDFR33.580.8787.2646.8134.040.948282.4053.63SRLRTR32.630.82398.4598.0833.560.92987.58283.91注：加粗字体表示各列最优结果。10.11834/jig.220306.T003表3实验3条件下，不同算法在 PaviaU 和 WDC 数据集上的降噪性能比较Table 3Quantitative comparison of different algorithms on PaviaU and WDC datasets denoising in experiment 3方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sRSLRNTF26.060.709207.99811.50726.050.82208.55894.08LRTDTV33.230.90187.02129.4632.960.9291.81260.91SSTVLRTF32.7330.825109.89572.93433.960.9385.781.46×103NLRCPTD26.420.77181.372.3×10426.040.79185.33.63×104SNLRSF3.0350.8751.9×1031.4×10329.720.92.11×103592.662LTDL27.610.81181.348.3×10326.890.84194.83.6×104LRTFDFR36.080.9265.7346.5835.710.9669.7355.28SRLRTR36.750.93663.12182.6537.850.9758.97287.08注：加粗字体表示各列最优结果。10.11834/jig.220306.T004表4实验4条件下，不同算法在 PaviaU 和 WDC 数据集上的降噪性能比较Table 4Quantitative comparison of different algorithms on PaviaU and WDC datasets denoising in experiment 4方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sRSLRNTF25.750.693213.79717.9125.830.80214.84979.73LRTDTV33.450.90487.93182.0432.760.9193.68263.30SSTVLRTF31.770.819112.32762.2432.990.92697.851.17×103NLRCPTD26.070.751187.813.6×10425.540.76202.74.31×104SNLRSF30.040.8792.10×103790.4929.610.9022.14×103556.23LTDL27.710.81184.859.3×10327.020.85192.53.1×104LRTFDFR36.170.93865.1848.4835.730.9669.6951.88SRLRTR36.560.94463.34211.80737.1810.97152.55314.45注：加粗字体表示各列最优结果。为了研究为什么SNLRSF、NLR-CPTD和LTDL等基于非局部相似性的方法在实验2中表现良好，添加了两个只有高斯噪声或脉冲噪声的比较案例，其中高斯噪声的方差为0.15，脉冲噪声的百分比为0.2。值得注意的是，从表 1—表4可以看出，非局部相似性的方法计算复杂度高，因此在PaciaU和WDC数据上进行了2倍下采样，使用SNLRSF、NLRCPTD、LTDL和SRLRTR进行降噪性能比较，结果如表5和表6所示。可以看出，基于非局部相似性的方法在高斯噪声上表现良好，而本文方法在脉冲噪声上的表现优于它们，这可能意味着非局部相似性的方法适用于去除高斯噪声，即便在去噪优化模型中考虑了脉冲噪声，例如SNLRSF算法。10.11834/jig.220306.T005表5不同方法在高斯噪声上的降噪性能比较Table 5The recovery performance on Gaussian noise using different methods方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sNLRCPTD31.300.861110.572.9×10327.280.85171.642.6×103SNLRSF34.390.9221.2×103180.3133.940.971.34×10335.87LTDL33.490.90886.384.5×10329.3750.91141.262.4×104SRLRTR32.040.863109.9243.1731.570.95108.2913.40注：加粗字体表示各列最优结果。10.11834/jig.220306.T006表6不同方法在脉冲噪声上的降噪性能比较Table 6The recovery performance on impulse noise using different methods方法PaviaU数据集WDC数据集MPSNR/dBMSSIMERGASCPU time/sMPSNR/dBMSSIMERGASCPU time/sNLRCPTD21.400.550358.603.1×10320.560.64393.82783.30SNLRSF26.880.8003.4×103196.6126.230.863.75×10338.59LTDL22.100.620338.754.5×10321.220.72378.002.4×104SRLRTR30.590.860123.4642.0830.760.95120.5313.49注：加粗字体表示各列最优结果。图4显示了WDC和PaviaU数据在实验4环境下通过不同方法恢复的实际波段的PSNR的比较。在实验4中，在不同波段中添加了各种类型的噪声。从图中可以看出，PSNR值在具有高斯噪声和脉冲噪声的1～41波段范围内变化平稳，而PSNR值在具有死线和条纹噪声50～190个波段内存在波动。就PSNR而言，本文方法的恢复性能几乎在每个波段上都优于其他方法，表明所提方法在去除混合噪声方面有一定优势。10.11834/jig.220306.F004图 4实验4环境下，不同降噪方法各波段对应恢复图像的 PSNR 比较结果Fig.4The comparison performance of different methods in terms of PSNR in case 4 for denoising（（a）WDC detaset；（b）PaviaU dataset）图5给出了实验1环境下WDC的第121个波段的图像比较结果，呈现了在混合高斯噪声和脉冲噪声条件下不同方法恢复的图像结果。可以看出，使用本文方法恢复的图像的分辨率高于其他方法恢复的图像的分辨率。10.11834/jig.220306.F005图5在实验1环境下，WDC数据的第121个谱段所对应的复原图像Fig.5The recovered 121st spectral segment of the WDC data in case 1 （（a） clean image； （b） noisy image； （c） LRTDTV；（d） RSLRNTF； （e） SNLRSF； （f） LRTFDFR； （g） NLRCPTD； （h） SSTVLRTF； （i） LTDL； （j） SRLRTR）图6给出了实验3环境下PaviaU的第85个波段的图像比较结果，呈现了混合高斯噪声、脉冲噪声和死线条件下的复原图像。可以看出，LTDL和RSLRNTF这两种方法恢复的图像仍然有死线噪声。10.11834/jig.220306.F006图6在实验3环境下，PaviaU数据的第85个谱段所对应的复原图像Fig.6The recovered 85th spectral segment of the PaviaU data in case 3 （（a） clean image； （b） noisy image； （c） LRTDTV；（d） RSLRNTF； （e） SNLRSF； （f） LRTFDFR； （g） NLRCPTD； （h） SSTVLRTF； （i） LTDL； （j） SRLRTR）3.4　真实数据实验结果在本实验中，使用EO-1和Urban数据集，这两个真实数据集的某些波段上有高斯噪声、条纹和死线严重的污染。本实验的目标是恢复这些污染的图像，同时不影响干净的图像。图 7和图 8分别展示了EO-1高光谱图像和Urban 高光谱图像在实际波段68和208中的降噪结果。从图 7可以看出，原始数据受到条纹噪声的严重污染，本文算法可以在保留高光谱图像的局部细节和结构信息的同时恢复高光谱图像。从图 8可以看出，场景被高斯噪声、条纹噪声和死线严重污染，只有SNLRSF和本文方法在恢复性能方面表现良好。10.11834/jig.220306.F007图7EO-1 数据集上的第68个波段所对应的恢复图像Fig.7The recovered 68th spectral segment of the EO-1 data （（a） noisy image； （b） LRTDTV；（c） RSLRNTF； （d） SNLRSF； （e） LRTFDFR； （f） NLRCPTD； （g） SSTVLRTF； （h） LTDL； （i） SRLRTR）10.11834/jig.220306.F008图8Urban 数据集上的第208个波段所对应的恢复图像Fig.8The recovered 208th spectral segment of the Urban data （（a） noisy image； （b） LRTDTV；（c） RSLRNTF； （d） SNLRSF； （e） LRTFDFR； （f） NLRCPTD； （g） SSTVLRTF； （h） LTDL； （i） SRLRTR）图 9显示了在不同实验环境下，所提方法在WDC和PaviaU两个数据上的收敛性能，验证了算法的收敛性，其中连续两个恢复图像之间的相对误差用来衡量收敛性。相对误差的定义为Xk-Xk-1FXk-1F，其中Xk表示在第k次迭代后的复原图像。从图9可以看出，在所有实验环境下，所提方法的相对误差在迭代的过程中会接近0，最终收敛。10.11834/jig.220306.F009图9不同实验情况下，算法收敛性能Fig.9The convergence performance on different cases（（a）WDC dataset；（b）PariaU dataset）4　讨论4.1　部分正交MVTF与传统MVTF的区别本文介绍了一种新的平滑且鲁棒的低秩部分正交MVTF分解方法，并将其用于高光谱图像降噪。部分正交结果的MVTF对建模干净的高光谱图像至关重要。与正交MVTF（de Lathauwer，2008）和非负MVTF（Qian等，2017）在内的传统MVTF相比，提出的MVTF在两个方面有所不同。1）传统正交MVTF（de Lathauwer，2008）的定义为X≈∑r=1RGr∘cr=∑r=1RDr×1Ar×2Br×3cr。其中，Gr=ArDrBrT是Gr的奇异值分解，cr的长度为1。在这种情况下，矩阵Ar和Br是正交的。与传统的正交MVTF相比，本文提出的部分正交MVTF，约束CTC=IR，C(:,r)=cr,r=1,⋯,R，在高光谱图像的线性子空间模型假设下，部分正交MVTF的因子具有物理意义。2）非负MVTF（Qian等，2017）的定义为𝒳≈∑r=1RGr∘cr=∑r=1RArBrT∘cr,Ar,Br,cr≥0。其中，Gr的秩为Lr, r=1,⋯,R，并且可以分解为ArBrT。这种低秩模型在高光谱图像解混（Qian等，2017）上显示了其优越性，其中cr被视为端元，矩阵Gr是相应的丰度图。基于该分解模型，Ding等人（2021）提出了用丰度矩阵的Schatten-p函数来刻画其低秩特性，并在高光谱超分辨效果上取得提升。另外，Xiong等人（2018）使用基于该分解模型的高光谱图像去噪算法RSLRNTF，但该模型需要预先定义秩Lr,r=1,⋯,R。在实际应用中，准确地确定第r个端元的Gr的秩是不可能的。此外，与RSLRNTF（Xiong等，2018）相比，本文方法在实验中表现出更好的恢复性能。4.2　参数的影响目标函数（4）中有4个参数，但只需要调整3个参数，因为这些参数代表目标函数中不同项的相对权重。在这种情况下，通过固定λG，分析参数λ𝒮、λN和λ TV。下面，以在实验4环境下，对WDC数据集进行实验为例，说明如何选择这些参数。1）为了调整参数λTV、λS和λN，设置λG=0.1。其他3个参数分别与局部平滑属性、稀疏噪声（即脉冲噪声、条纹和死线）和高斯噪声的权重有关。图10将MPSNR值显示为与λTV和λS相关的函数， λN、λTV和λS从集合0.01, 0.018, 0.032,  0.1， 0.000 1, 0.000 17, 0.000 46, 0.001 3, 0.003 6, 0.01和  0.001,0.002 8,0.007 7,0.02,0.05,0.1中选择。很明显，λTV越小，λN越大，MPSNR的值就越大。当λS=0.02时，所提方法可以达到MPSNR的峰值。10.11834/jig.220306.F010图 10权重λTV、λS和λN对 MPSNR 的影响Fig.10The impact of parameters λTV，λS and λN （（a）λN = 0.01； （b）λN = 0.018； （c）λN = 0.032； （d）λN = 0.1）2）参数R的影响。参数R在算法中表示子空间维度的大小。图 12将MPSNR值显示为R的函数，其中R的变化为2～15，步长为1。从图11可以看到，当R = 5时，MPSNR值达到峰值。10.11834/jig.220306.F011图 11参数R对MPSNR的影响Fig.11The impact of parameter R5　结论本文提出了一种新的用于高光谱图像去噪的部分正交低秩MVTF模型，该模型具有较强的线性子空间模型的解释能力，并可以有效地探索高光谱图像的全局空间光谱结构。另外，总变分算子用于探索高光谱图像的局部相似结构，优化模型中𝓁1-范数和F-范数约束分别用于去除稀疏和高斯噪声。实验表明，与现有性能优异的算法相比，所提算法在混合噪声情况下各个波段都具有良好的恢复性能。在处理高斯噪声时，基于非局部相似性的方法恢复性能更好，但是计算复杂度非常高。在未来的工作中，可以在高斯噪声的情况下考虑加入非局部相似性先验信息，提高降噪的性能。同时，开发一种快速的非局部低秩恢复算法也是一个有意义的研究方向。