0引言高光谱图像(hyperspectral image, HSI)由成像光谱仪收集的地球表面同一区域数百条从紫外线到红外波长的光谱波段图像组成, 具有非常丰富的光谱信息, 广泛应用于环境监测(Frick和Tervooren，2019)、地物分类(冉琼等，2018；Li等，2020；谭琨等，2019)等各个领域。然而, 受实际的成像条件、天气条件或数据传输程序的限制, 实际的HSI往往受到不同类型的污染, 即噪声(杨润宇等，2019；Zeng等，2020)、欠采样(Zhang等，2020a)或缺失的数据(Xie等，2019)。这些类型的污染严重降低了HSI的质量, 限制了后续处理任务的精度, 如HSI解混(Hong等，2019a)、目标检测(Lin等，2018；Zhang等，2020b)、HSI分类(Hong等，2019b, 2021a, b)等。因此, 作为HSI应用的预处理步骤, HSI复原是一个有待深入研究且具有较强热度的研究方向。复原HSI的一种常见方法是将每一个波段视为灰度图像或将HSI列化成2维矩阵, 利用矩阵复原的方法进行处理。秩最小化是解决这类问题的常用策略。但由于秩函数是离散且非凸的, 因此秩函数最小化问题的求解是一个非确定性多项式(nondeterministic polynominal, NP)难问题。为了解决这个问题, 通常使用秩函数的凸包络—核范数去近似秩函数, 将秩最小化问题转化为核范数最小化问题。虽然核范数能取得良好的复原结果, 但由于凸性的限制, 对矩阵所有奇异值施加相同的收缩, 导致对于秩函数的逼近程度受限。为了提高逼近精度, 不少学者提出利用非凸函数去近似秩函数, 例如: Gu等人(2017)提出加权核范数来对不同的奇异值进行不同程度的收缩。Xie等人(2016)利用logsum函数, Ji等人(2017)利用logDet函数, Zeng等人(2021b)利用 $\gamma $ 函数, 以及Zhang(2010)利用MCP(minimax concave penalty)函数来实现秩函数的非凸逼近。这种方式考虑了矩阵大奇异值主要包含图像信息, 而小奇异值主要包含异常值信息的物理特点, 并针对不同的奇异值施加不同程度的惩罚, 从而改进对秩函数的逼近程度。然而, 尽管基于矩阵方法的模型取得了不错的结果(Gu等，2017；Xie等，2016；Zhang，2010；Liu等，2019), 但是HSI是对同一场景在不同光谱下的成像, 每一个波段的数据都具有高度相关性。通过矩阵列化的方式往往会忽略这种相关性, 对于HSI的低秩特性探索不够完全。Lu等人(2020)，Zhang等人(2020a)等研究在处理类似HSI的高阶张量数据方面, 直接对张量建模会优于矩阵的方法。主流的张量处理方法是基于CP (Candecomp/Parafac) 分解和Tucker分解的。张量的CP分解(Hitchcock，1927)是将张量分解成多个秩1张量的和, 那么张量的CP秩可以定义为秩1张量个数的最小值。尽管张量CP秩是矩阵秩的直接推广, 但是最小化CP秩或者CP秩松弛形式(Friedland和Lim，2018)仍然是NP难问题。Tucker分解(Tucker，1966)是将张量分解为核心张量和因子矩阵模态相乘的形式。基于此, Tucker秩可以定义为张量各模态展开矩阵的秩所组成的向量。Liu等人(2013)将张量各模态展开矩阵核范数之和(sum of nuclear norm，SNN)定义为Tucker秩的凸松弛形式。Liu等人(2021)将其成功应用于高光谱融合复原问题中。Kilmer等人(2013)提出了张量奇异值分解(tensor singular value decomposition, T-SVD)方法, 并据此给出了张量的multi-秩和tubal-秩的定义。Semerci等人(2014)将multi-秩凸松弛为张量核范数(tensor nuclear norm, TNN), 即张量模-3傅里叶变换后各正向切片矩阵核范数的平均。由于TNN独特的计算方法和对时间序列特征的高效利用, 因此受到广泛的关注, 并成功应用于HSI复原问题中。然而, TNN的定义继承了矩阵核范数将所有奇异值同等对待的缺陷，并忽略了张量的频率信息: 在张量模-3傅里叶变换域中, 对清晰的高光谱图像添加噪声，图像信息在低频部分变化较小，而在高频部分变化巨大。因此, 其对张量秩的逼近精度有限。针对上述问题, 本文提出一种基于频率加权张量核范数。所提方法的基本构建思想是通过适当减少对低频正向切片矩阵的惩罚, 加大对高频正向切片矩阵的惩罚, 从而提高张量核范数对张量秩的逼近程度, 以达到提高复原精度的目的。同时, 给出各个频率正向切片矩阵的权重自适应计算方法, 增加基于频率加权张量核范数对数据的自适应性, 减少参数调整过程的人为干预。最后，将基于频率加权张量核范数应用于高光谱复原和去噪问题, 并在多个高光谱数据实验上成功验证了所提出模型的有效性。1相关工作1.1频率带和频率成分对于任意实向量 $\mathit{\boldsymbol{v}} = {\left[ {{v_1}, {v_2}, \cdots, {v_N}} \right]^{\rm{T}}}$ , 其快速傅里叶变换记为 $\mathit{\boldsymbol{\bar v = }}\mathit{fft}\left(\mathit{\boldsymbol{v}} \right) \in {{\bf{C}}^N}$ , 相应的快速傅里叶逆变换记为 $\mathit{\boldsymbol{v = }}\mathit{ifft}\left({\mathit{\boldsymbol{\bar v}}} \right)$ 。那么, 频率带(Wang等，2020)的定义为 1 $\{\bar{v}\}_{n-1}= \begin{cases}\left(\bar{v}_{n}\right) & n=1 \text { 或 } n=N-n+2 \\ \left(\bar{v}_{n}, \bar{v}_{N-n+2}\right) & \text { 其他 }\end{cases}$ 当 $N$ 为奇数时, 频率带为 $\left({{{\bar v}_1}} \right), \left({{{\bar v}_2}, {{\bar v}_N}} \right), \cdots, \left({{{\bar v}_{\frac{{N + 1}}{2}}}, {{\bar v}_{\frac{{N + 3}}{2}}}} \right)$ ；当 $N$ 为偶数时, 频率带为 $\left({{{\bar v}_1}} \right), \left({{{\bar v}_2}, {{\bar v}_N}} \right), \cdots, \left({{{\bar v}_{\frac{{N + 2}}{2}}}} \right)$ 。如果单独将频率带 ${\left\{ {\mathit{\bar v}} \right\}_{n - 1}}$ 进行傅里叶逆变换, 那么就得到频率成分 ${\left[ v \right]_{n - 1}} \in {{\bf{R}}^N}$ , 即 2 $\begin{gathered}{[v]_{n-1}=i f f t\left(\{\bar{v}\}_{n-1}\right)=} \\{ifft}\left(\left[0, \cdots, 0, \bar{v}_{n}, 0, \cdots, 0, \bar{v}_{N-n+2}, 0, \cdots, 0\right]^{\mathrm{T}}\right)\end{gathered}$ 在频率带式(1)和频率成分式(2)的定义下，向量 $\mathit{\boldsymbol{v}}$ 可以表示为所有频率成分之和，即 3 $\boldsymbol{v}=\sum\limits_{i=0}^{N_{v}-1} {ifft}\left(\{\bar{v}\}_{i}\right)=\sum\limits_{i=0}^{N_{v}-1}[v]_{i}$ 式中， ${N_v} = \left\lceil {\left({N + 1} \right)/2} \right\rceil $。 针对任意张量 $\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}$ , 其模-3傅里叶变换记为 $\overline {\mathit{\boldsymbol{\tilde X}}} = \mathit{fft}\left({\mathit{\boldsymbol{\tilde X}}, \left[ {} \right], 3} \right)$ , 相应的模-3傅里叶逆变换记为 $\mathit{\boldsymbol{\tilde X}} = i\mathit{fft}\left({\mathit{\boldsymbol{\tilde X}}, \left[ {} \right], 3} \right)$ 。基于频率带式(1)和频率成分式(2)的定义, 张量 ${\mathit{\boldsymbol{\tilde X}}}$ 也可以表示为所有频率成分之和的形式, 即 4 $\tilde{\boldsymbol{X}}=\sum\limits_{i=0}^{N_{v}-1} {ifft}\left(\{\tilde{\boldsymbol{X}}\}_{i}\right)=\sum\limits_{i=0}^{N_{v}-1}[\tilde{\boldsymbol{X}}]_{i}$ 式中, ${\left\{ {\mathit{\boldsymbol{\tilde X}}} \right\}_i}$ 表示张量 ${\mathit{\boldsymbol{\tilde X}}}$ 的第 $i$ 个频率带; ${\left[ {\mathit{\boldsymbol{\tilde X}}} \right]_i}$ 表示张量 $\mathit{\boldsymbol{X}}$ 的第 $i$ 个频率成分; ${N_v} = \left\lceil {\left({N + 1} \right)/2} \right\rceil $ 。1.2秩函数逼近矩阵或者张量的秩函数离散且非凸, 相应的秩最小化问题的求解是NP难问题。为了克服计算上的困难, 通常的做法是将非凸的秩函数近似成凸函数, 然后代入最小化问题中求解。Fazel(2002)证明了矩阵核范数是矩阵秩函数在矩阵谱范数单位球上的凸包络。此后, 核范数作为秩函数的凸替代, 而秩最小化问题转换成核范数最小化问题(nuclear norm minimization, NNM), 广泛应用于高光谱图像处理领域。针对任意矩阵 $\mathit{\boldsymbol{X}} \in {{\bf{R}}^{m \times n}}$ 相应的矩阵范数定义为 5 $\|\boldsymbol{X}\|_{*}=\sum\limits_{i} \sigma_{i}(\boldsymbol{X})$ 式中, $\sigma_{i}(\boldsymbol{X})$ 为矩阵的第 $i$ 大奇异值。针对矩阵核范数最小化问题, Cai等人(2010)提出了核范数逼近方法, 并使用软阈值算法来求解。尽管NNM在不少方面取得了成功(Lin等，2009；Ji和Ye，2009), 但其将所有的奇异值同等对待, 且用相同的阈值进行收缩处理, 从而未充分考虑矩阵奇异值大小的物理意义, 即较大的奇异值蕴含的是图像主体信息, 较小的奇异值蕴含的是图像异常值信息。为了克服这个缺陷, Gu等人(2017)提出了矩阵加权核范数最小化, 即 6 $\|\boldsymbol{X}\|_{w, *}=\sum\limits_{i} w_{i} \sigma_{i}(\boldsymbol{X})$ 式中, $w_i$ 为第 $i$ 大奇异值的权重。对于HSI复原和去噪问题, 所研究的对象可以建模成张量。针对张量而言, 一种简单的方式是将张量展开成矩阵, 然后利用矩阵低秩问题的方法进行处理(He等，2016, 2018；Zeng等，2021a)。例如, 针对任意三阶张量 $\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}$ , Liu等人(2013)在张量Tucker分解的基础上定义了张量的SNN，即 7 $\|\tilde{\boldsymbol{X}}\|_{*}=\sum\limits_{k=1}^{3} \alpha_{k}\left\|\boldsymbol{X}_{(k)}\right\|_{*}$ 式中, $\boldsymbol{X}_{(k)}$ 是张量 ${\mathit{\boldsymbol{\tilde X}}}$ 的模- $k$ 展开矩阵; ${\alpha _k}$ 为权重参数, 满足 $\sum\limits_{k = 1}^3 {{\alpha _k} = 1} $ 。然而, 张量具有内在的高维结构, 将张量展开成矩阵处理的方式会破坏张量内部的相关性, 并导致信息丢失(Zhou等，2018)。为了解决上述问题, 在T-SVD方法(Kilmer等，2013)的基础上, Lu等人(2020)给出了三阶张量 $\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}$ 的TNN定义为 8 $\|\tilde{\boldsymbol{X}}\|_{\mathrm{TNN}}=\frac{1}{n_{3}} \sum\limits_{i=1}^{n_{3}}\left\|\overline{\tilde{\boldsymbol{X}}}{}^{(i)}\right\|_{*}$ 式中, $\overline{\tilde{\boldsymbol{X}}}{}^{(i)}$ 表示 $\overline{\tilde{\boldsymbol{X}}}$ 的第 $i$ 个正向切片矩阵。与张量的SNN相比, TNN直接在张量上定义, 没有破坏张量内在的高维结构, 对张量的秩有更好的逼近能力。但是, TNN对所有奇异值同等对待, 忽略了奇异值对张量信息的表达特征。因此, Jiang等人(2020)提出部分和张量核范数(partial sum of the tensor nuclear norm, PSTNN), 只对较小的奇异值进行惩罚, 放弃对较大奇异值的收缩, 在一定程度上弥补了TNN的上述缺陷，还有一些类似研究，如Zhang(2010)，Xie等人(2016)，Gu等人(2017)，Zeng等人(2021b)等。然而, 这些文献聚焦于研究矩阵内部奇异值大小所蕴含的物理信息, 并没有考虑张量模-3傅里叶域各频率带所蕴含的物理信息。2频率加权张量核范数傅里叶变换在TNN(式(8))中起到关键的作用，本小节将通过真实数据实验的形式分析张量 ${\mathit{\boldsymbol{\tilde X}}}$ 在模-3傅里叶变换域中的性质, 并以此给出频率域加权张量核范数定义及其权重参数的自适应确定方法。2.1频率加权张量核范数在图 1中, 对真实的高光谱遥感Pavia数据集进行模-3傅里叶变换, 计算出变换后各个频率带切片矩阵奇异值, 并绘制出所有奇异值关于频率带的分布曲线。从图 1中可以明显观察到较大的奇异值主要集中在低频率带, 而较小的奇异值主要分布在高频率带中。利用不同的频率切片进行傅里叶逆变换, 得到不同频率成分张量, 计算各个频率成分张量的峰值信噪比(peak signal-to-noise ratio, PSNR), 并将其绘制在图 2中。从图 2中可以明显地观察到低频率成分的PSNR高于高频率成分的PSNR。此外, 针对清晰的高光谱遥感Indian数据添加稀疏噪声形成噪声数据, 然后分离出噪声数据的零频率成分和非零频率成分。图 3显示了所有数据的第15波段的灰度图像，从图 3中可以明显观测到零频率成分包含了清晰数据的绝大部分信息, 而非零频率成分主要包含了噪声信息。上述真实数据的实验分析表明, 张量 ${\mathit{\boldsymbol{\tilde X}}}$ 主要信息集中在 $\overline {\mathit{\boldsymbol{\tilde X}}} $ 低频切片矩阵中, 而异常值信息则集中在 $\overline {\mathit{\boldsymbol{\tilde X}}} $ 的高频切片矩阵中。 图1 各个频率带的奇异值分布情况 The singular value distribution of each frequency bandFig 1 图2 利用不同频率带复原图像的PSNR PSNR of the restored image with different frequency bandsFig 2 图3 Indian数据的第15波段频率带复原情况 Restoration of the 15th frequency band of Indian dataFig 3((a) clear data; (b) noise data; (c) zero frequency band; (d) non-zero frequency band) 在上述实验分析的基础上, 针对高光谱遥感Pavia数据添加强度为0.1, 0.3和0.5的稀疏噪声, 并在图 4中绘制出各频率带切片矩阵核范数的分布曲线。由图 4可知, 噪声对低频率带切片矩阵核范数的影响较小, 而对高频率带切片矩阵核范数影响较大; 同时, 在噪声强度增加的情况下, 低频率带切片矩阵核范数变化较小, 而高频率带切片矩阵核范数变化较大。因此, 在利用TNN对张量进行复原时, 若能减少对低频率带切片矩阵核范数的收缩, 同时加大高频率带切片矩阵核范数的惩罚, 将更符合张量在频率域的物理性质, 帮助提升复原结果的精度。基于上述真实数据实验分析的结果, 在张量 $\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}$ 的TNN定义过程中, 根据高低频率带的物理含义, 对不同频率带切片矩阵核范数加权求和, 提出频率加权张量核范数(frequency weighted tensor nuclear norm, FWTNN)，即 9 $\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}=\sum\limits_{i=1}^{n_{3}} w_{i}\left\|\overline{\tilde{\boldsymbol{X}}}{}^{(i)}\right\|_{*}$ 图4 受噪声污染前后各个频率带核范数的变化情况 The variation of the nuclear norm of each frequency band before and after noise pollutionFig 4式中, ${w_i}, i = 1, 2, \cdots, {n_3}$ 是频率域权重参数, 表示各个频率切片的惩罚程度。合理的频率域权重参数, 不仅需要考虑不同频率切片所携带的原始张量信息, 而且要兼顾异常值对各个频率带切片矩阵核范数的影响。这是加权核范数和非凸核范数(Zhang，2010；Xie等，2016；Gu等，2017；Zeng等，2020)所没有考虑的物理意义。2.2频率域权重参数在频率加权张量核范数式(9)中, 频率域权重参数的取值非常重要。如果权重参数取值均为1, FWTNN就退化为一般的TNN。如果频率域权重参数 $w_i$ 能在携带张量主要信息的低频带取较小的值, 就可以减少对低频带切片矩阵核范数的惩罚, 在去除异常值的条件下保护原始张量的主要信息; 而在包含异常值的高频率带取较大的值, 那么就可以增加对高频带切片矩阵核范数的惩罚, 从而充分地去除张量中的异常值。由图 4可知, 清晰高光谱遥感图像各个频率带核范数与频率域权重参数的分布特点相反, 即在低频带处核范数取值较大, 而在高频率带处核范数取值较小, 并且核范数取值随频率增加而减少。因此, 可以使用各个频率切片矩阵核范数对数的倒数作为频率域权重参数, 即 10 $w_{i}=\frac{1}{\lg \left(\left|\overline{\tilde{X}}{}^{(i)}\right|_{*}\right)}, i=1, 2, \cdots, n_{3}$ 式中, 对数函数的作用是对奇异值进行尺度收缩, 防止最大奇异值的极端影响。上述方法在原始清晰HSI已知的情况下是可行的。然而, 针对HSI复原和去噪问题, 往往仅有含异常值的观测数据, 原始清晰的HSI未知。如果使用观测数据来计算频率域权重, 由图 5可知, 这样的频率域权重与清晰图像计算的频率域权重相差较大, 而复原数据所计算出的频率域权重与清晰图像计算的频率域权重相差较小。进一步地, 由图 6可知, 随着迭代次数 $k$ 的不断增加, 复原图像所计算出的频率域权重会越来越接近清晰图像所计算的频率域权重。据此, 本文提出迭代的权重参数计算方式。假定第 $k$ 次迭代得到 ${{{\mathit{\boldsymbol{\tilde X}}}^k}}$ , 那么第 $k+1$ 次迭代对应的权重参数记为 11 $w_{i}^{k+1}=C_{1} h_{i}^{k}+\frac{C_{2}}{k} \quad i=1, 2, \cdots, n_{3}$ 图5 Pavia数据各个波段的权重 The weight of each band of Pavia dataFig 5 图6 不同迭代次数所对应的频率域权重 Frequency weights corresponding to different iteration timesFig 6式中, $h_i^k = 1/\lg \left( {{{\left\| {{{\overline {{\mathit{\boldsymbol{\tilde X}}}} }^{k\left( i \right)}}} \right\|}_*}} \right);{C_1}$ 为频率归一化后的缩放系数; $C_2$ 为一个较小的常数, 保证在初始迭代时对所有的频率切片矩阵核范数都进行惩罚。由式(11)可知, 随着迭代次数 $k$ 的不断增加, ${C_2}/k$ 的值越来越小, 对所有频率切片统一收缩力度不断减小, 频率权重 ${w^{k + 1}}$ 的主导也渐渐由 ${C_2}/k$ 转化为 ${C_1}{h^k}$ 。这样设置权重的原因是在迭代初期, 所有的频率切片都含有异常值信息, 因此需要一个较小的收缩常数来确保对每一个频率切片进行收缩。随着迭代次数 $ k$ 的增加, 复原得到的图像越来越接近清晰图像, 低频所需的收缩也就越来越小; 而高频部分仍需要依靠 ${C_1}{h^k}$ 主导的权重来收缩。图 7展示了频率加权张量核范数的示意图。 图7 频率加权张量核范数示意图 Schematic diagram of frequency weighted tensor nuclear normFig 7((a) model the hyperspectral data as a third-order tensor ${\mathit{\boldsymbol{\tilde X}}}$ and perform Fourier transform on it; (b) frequency components and singular value distribution of ${\overline {\mathit{\boldsymbol{\tilde X}}} }$ ; (c) weight the frequency components to get ${{{\mathit{\boldsymbol{\tilde X}}}_{{\rm{FWTNN}}}}}$ ) 3本文模型与求解方法FWTNN反映了张量在频率域的物理特征, 能提供对张量秩函数更精确的逼近, 本文将FWTNN应用于HSI复原和去噪问题。3.1基于FWTNN的HSI复原模型针对清晰HSI， $\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}$ , 成像过程中由于天气因素或者传感器破坏等原因, 观测到的高光谱图像 ${\mathit{\boldsymbol{\tilde O}}}$ 往往包含一些缺失信息。记观测索引集为 $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$ ，在 $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$ 中, 有 ${{\mathit{\boldsymbol{\tilde O}}}_\mathit{\Omega }}\mathit{\boldsymbol{ = }}{{\mathit{\boldsymbol{\tilde X}}}_\mathit{\Omega }}$ ; 而在 $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} $ 的补集 ${\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}^C}$ 中, ${{\mathit{\boldsymbol{\tilde O}}}_{{\mathit{\Omega }^C}}} = 0$ 。那么, 基于FWTNN的HSI复原模型可以表示为 12 $\underset{\tilde{\boldsymbol{X}}}{\arg \min }\|\tilde{\boldsymbol{X}}\|_{\text {FWTNN }} \text { s.t. } \tilde{\boldsymbol{X}}_{\varOmega}=\tilde{\boldsymbol{O}}_{\varOmega}$ 当 ${\left\| {\mathit{\boldsymbol{\tilde X}}} \right\|_{{\rm{FWTNN}}}}$ 中的权重 $w_i$ 都为1时, 该模型退化为基于TNN的HSI复原模型。令 13 $I_{\varPhi}(\tilde{\boldsymbol{X}})= \begin{cases}0 & \tilde{X} \in \varPhi \\ \infty & \text { 其他 }\end{cases}$ 式中, $\mathit{\boldsymbol{ \boldsymbol{\varPhi} = }}\left\{ {\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}:{{\mathit{\boldsymbol{\tilde X}}}_\mathit{\Omega }} = {{\mathit{\boldsymbol{\tilde O}}}_\mathit{\Omega }}} \right\}$ 。基于式(13), 将式(12)重写为 14 $\underset{\tilde{\boldsymbol{X}}}{\arg \min } I_{\varPhi}(\tilde{\boldsymbol{X}})+\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}$ 引入辅助变量 ${\mathit{\boldsymbol{\tilde Y}}}$ , 令正则项 ${I_\mathit{\Phi }}\left({\mathit{\boldsymbol{\tilde X}}} \right)$ 中的 ${\mathit{\boldsymbol{\tilde X}} = \mathit{\boldsymbol{\tilde Y}}}$ , 最小化问题(14)可以改写为 15 $\underset{\tilde{\boldsymbol{X}}}{\arg \min } I_{\varPhi}(\tilde{\boldsymbol{Y}})+\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}} \text { s.t. } \tilde{\boldsymbol{Y}}=\tilde{\boldsymbol{X}}$ 利用增广拉格朗日乘子法, 最小化模型(15)中的目标函数可以转化为 16 $\begin{gathered}L_{\mu}(\tilde{\boldsymbol{X}}, \tilde{\boldsymbol{Y}}, \tilde{\boldsymbol{\varGamma}})= \\I_{\varPhi}(\tilde{\boldsymbol{Y}})+\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\langle\tilde{\boldsymbol{\varGamma}}, \tilde{\boldsymbol{X}}-\tilde{\boldsymbol{Y}}\rangle+\frac{\mu}{2}\|\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{Y}}\|_{\mathrm{F}}^{2}= \\I_{\varPhi}(\tilde{\boldsymbol{Y}})+\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\frac{\mu}{2}\left\|\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{Y}}+\frac{\tilde{\boldsymbol{\varGamma}}}{\mu}\right\|_{\mathrm{F}}^{2}\end{gathered}$ 式中, ${\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}$ 是拉格朗日乘子, $\mu $ 是罚参数。对于 $L_{\mu}(\tilde{\boldsymbol{X}}, \tilde{\boldsymbol{Y}}, \tilde{\boldsymbol{\varGamma}})$ 的最小化问题, 可以使用交替方向乘子法(Boyd等，2011)(alternating direction method of multipliers, ADMM)将其转换成以下的子优化问题来求解。已知第 $K$ 次迭代结果 ${{{\mathit{\boldsymbol{\tilde X}}}^k}}$ , ${{{\mathit{\boldsymbol{\tilde Y}}}^k}}$ 和 ${{\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}^k}$ , 那么第 $k+1$ 次迭代结果的求解过程如下：1) 针对 ${\mathit{\boldsymbol{\tilde X}}}$ 的子优化问题 17 $\tilde{\boldsymbol{X}}^{k+1}=\underset{\tilde{\boldsymbol{X}}}{\arg \min }\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\frac{\mu}{2}\left\|\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{Y}}^{k}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right\|_{\mathrm{F}}^{2}$ 可以通过FTSVT(frequency-weighted tensor singulor value threshad)算子(Wang等，2020)求解, 结果表示为 18 $\tilde{\boldsymbol{X}}^{k+1}=F W_{\frac{w^k}{\mu}}\left(\tilde{\boldsymbol{Y}}-\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right)$ 式中， $FW\left(\cdot \right)$ 为FTSVT算子。2) 针对 ${\mathit{\boldsymbol{\tilde Y}}}$ 的子优化问题 19 $\tilde{\boldsymbol{Y}}^{k+1}=\arg \underset{\tilde{\boldsymbol{Y}}}{\min} I_{\varPhi}(\tilde{\boldsymbol{Y}})+\frac{\mu}{2}\left\|\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{Y}}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right\|_{\mathrm{F}}^{2}$ 解为 20 $\tilde{\boldsymbol{Y}}^{k+1}=\frac{1}{\mu}\left(\mu \tilde{\boldsymbol{X}}^{k+1}+\tilde{\boldsymbol{\varGamma}}^{k}\right)_{\varOmega^{C}}+\tilde{\boldsymbol{O}}_{\varOmega}$ 3) 针对拉格朗日乘子 ${\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}$ 的更新 21 $\tilde{\boldsymbol{\varGamma}}^{k+1}=\tilde{\boldsymbol{\varGamma}}^{k}+\mu\left(\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{Y}}^{k+1}\right)$ 算法1给出了基于FWTNN的HSI复原模型式(12)的详细计算过程。算法复杂度为 ${\rm{O}}\left({{n_1}{n_2}{n_3}\log \left({{n_3}} \right) + {n_1}{n_2}{n_3} + {n_{\left(2 \right)}}{n_3}} \right)$ 。算法1    基于FWTNN的HSI复原求解算法输入: 观测数据 ${\mathit{\boldsymbol{\tilde O}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ ; 索引集合 $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$ ; 停止条件 $\varepsilon $ ; 罚参数更新速率 $\rho $ 。输出: 复原HSI, ${\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ 。过程1初始化: ${{\mathit{\boldsymbol{\tilde X}}}^0} = RAND\left({{n_1} \times {n_2} \times {n_3}} \right), \mathit{\boldsymbol{\tilde X}}_\mathit{\Omega }^0 = {{\mathit{\boldsymbol{\tilde O}}}_\mathit{\Omega }}, {{\mathit{\boldsymbol{\tilde Y}}}^0} = {{\mathit{\boldsymbol{\tilde X}}}^0}$ , $\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }} = ZEROS\left({{n_1} \times {n_2} \times {n_3}} \right)$ , 频率权重 $w$ , 罚参数 $\mu $ 。过程2循环迭代:(1) 基于式(18)更新 ${\mathit{\boldsymbol{\tilde X}}}$ ;(2) 基于式(20)更新 ${\mathit{\boldsymbol{\tilde Y}}}$ ;(3) 基于式(21)更新 ${\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}$ ;(4) 基于式(11)更新 $w$ ;(5) 更新罚参数 ${\mu ^{k + 1}} = \rho {\mu ^k}$ ;(6) 检查停止条件: $\left\|\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{X}}^{k}\right\|_{\infty}\varepsilon, \left\|\tilde{\boldsymbol{Y}}^{k+1}-\tilde{\boldsymbol{Y}}^{k}\right\|_{\infty}\varepsilon, $ $\left\|\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{Y}}^{k+1}\right\|_{\infty}\varepsilon。$ 注：RAND()和ZEROS()为MATLAB内置函数，其含义分别为产生随机矩阵和零矩阵。3.2基于FWTNN的HSI去噪模型实际的HSI在成像过程中也会受到不同类型的噪声干扰。而本文的HSI去噪模型主要针对的是含有稀疏噪声的HSI进行去噪处理。针对三阶张量 ${\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ , 由稀疏噪声 ${\mathit{\boldsymbol{\tilde S}}}$ 污染后得到观测张量为 ${\mathit{\boldsymbol{\tilde O}}}$ , 即 $\tilde{\boldsymbol{O}}=\tilde{\boldsymbol{X}}+\tilde{\boldsymbol{S}}$ 那么, 基于FWTNN的HSI去噪模型可以表示为 22 $\underset{\boldsymbol{\widetilde{X}, \widetilde{S}}}{\arg \min }\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\lambda\|\tilde{\boldsymbol{S}}\|_{1} \text { s.t. } \tilde{\boldsymbol{O}}=\tilde{\boldsymbol{X}}+\tilde{\boldsymbol{S}}$ 利用增广拉格朗日乘子法, 优化问题式(22)中的目标函数可以转化为 23 $\begin{gathered}L_{\mu}(\tilde{\boldsymbol{X}}, \tilde{\boldsymbol{S}}, \tilde{\boldsymbol{\varGamma}})=\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\lambda\|\tilde{\boldsymbol{S}}\|_{1}+ \\\langle\tilde{\boldsymbol{\varGamma}}, \tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{S}}\rangle+\frac{\mu}{2}\|\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{S}}\|_{\mathrm{F}}^{2}= \\\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+\lambda\|\tilde{\boldsymbol{S}}\|_{1}+ \\\frac{\mu}{2}\left\|\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{S}}+\frac{\varGamma}{\mu}\right\|_{\mathrm{F}}^{2}\end{gathered}$ 对于式(23)的最小化问题, 可以使用ADMM将其转换成以下的子优化问题来求解。已知第 $k$ 次迭代结果 ${{{\mathit{\boldsymbol{\tilde X}}}^k}}$ , ${{{\mathit{\boldsymbol{\tilde S}}}^k}}$ 和 ${{\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}^k}$ , 那么第 $k+1$ 次迭代结果的求解过程如下：1) 针对 ${\mathit{\boldsymbol{\tilde X}}}$ 的子优化问题 24 $\begin{gathered}\tilde{\boldsymbol{X}}^{k+1}=\arg\underset{\tilde{\boldsymbol{X}}}{ \min }\|\tilde{\boldsymbol{X}}\|_{\mathrm{FWTNN}}+ \\\frac{\mu}{2}\left\|\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}-\tilde{\boldsymbol{S}}^{k}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right\|_{\mathrm{F}}^{2}\end{gathered}$ 可以通过FTSVT算子(Wang等，2020)求解，结果表示为 25 $\tilde{\boldsymbol{X}}^{k+1}=F W_{\frac{w^k}{\mu}}\left(\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{S}}^{k}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right)$ 2) 针对 ${\mathit{\boldsymbol{\tilde S}}}$ 的子优化问题 26 $\tilde{\boldsymbol{S}}^{k+1}=\arg \underset{\widetilde{\boldsymbol{S}}}{\min } \lambda\|\tilde{\boldsymbol{S}}\|_{1}+\frac{\mu}{2}\left\|\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{S}}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right\|_{\mathrm{F}}^{2}$ 可以通过软阈值收缩算子 $R$ (Zeng等，2020)来求解 27 $\tilde{\boldsymbol{S}}=R_{\frac{\lambda}{\mu}}\left(\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}^{k+1}+\frac{\tilde{\boldsymbol{\varGamma}}^{k}}{\mu}\right)$ 式中, 针对 $x \in {\bf{R}}$ , 且 $\Delta 0$ , 软阈值收缩算子为 28 $R_{\Delta}(x)= \begin{cases}x-\Delta & x\Delta \\ x+\Delta & x-\Delta \\ 0 & \text { 其他 }\end{cases}$ 3) 针对拉格朗日乘子 ${\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}$ 的更新 29 $\tilde{\boldsymbol{\varGamma}}^{k+1}=\tilde{\boldsymbol{\varGamma}}^{k}+\mu\left(\tilde{\boldsymbol{O}}-\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{S}}^{k+1}\right)$ 算法2给出了基于FWTNN的HSI去噪模型式(22)的详细计算过程。算法复杂度为 ${\rm{O}}\left({{n_1}{n_2}{n_3}\log \left({{n_3}} \right) + {n_{\left(1 \right)}}n_{\left(2 \right)}^2{n_3} + {n_{\left(2 \right)}}{n_3}} \right)$ , 其中 ${n_{\left(1 \right)}} = \max \left({{n_1}, {n_2}} \right), {n_{\left(2 \right)}} = \min \left({{n_1}, {n_2}} \right)$ 。算法2    基于FWTNN的HSI去噪求解算法输入: 观测数据 ${\mathit{\boldsymbol{\tilde O}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ ; 参数 $\lambda $ ; 停止条件 $\varepsilon $ ; 罚参数更新速率 $\rho $ 输出: 复原数据 ${\mathit{\boldsymbol{\tilde X}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ 过程1 初始化:${{\mathit{\boldsymbol{\tilde X}}}^0} = \mathit{\boldsymbol{\tilde O}}, \mathit{\boldsymbol{\tilde S}} = \mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }} = ZEROS\left({{n_1} \times {n_2} \times {n_3}} \right)$ ; 频率权重 $w$ ; 罚参数 $\mu $过程2  循环迭代:(1) 基于式(25)更新 ${\mathit{\boldsymbol{\tilde X}}}$ ;(2) 基于式(27)更新 ${\mathit{\boldsymbol{\tilde S}}}$ ;(3) 基于式(29)更新 ${\mathit{\boldsymbol{ \boldsymbol{\tilde \varGamma} }}}$ ;(4) 基于式(11)更新 $w$ ;(5) 更新罚参数 ${\mu ^{k + 1}} = \rho {\mu ^k}$ ;(6) 检查停止条件 $\left\|\tilde{\boldsymbol{X}}^{k+1}-\tilde{\boldsymbol{X}}^{k}\right\|_{\infty}\varepsilon, \left\|\tilde{\boldsymbol{S}}^{k+1}-\tilde{\boldsymbol{S}}^{k}\right\|_{\infty}\varepsilon, $ $\left\|\tilde{\boldsymbol{X}}^{k+1}+\tilde{\boldsymbol{S}}^{k+1}-\tilde{\boldsymbol{O}}\right\|_{\infty}\varepsilon。$ 4实验结果与分析为了验证本文模型的有效性, 选取3种相关的HSI复原方法进行对比, 并用峰值信噪比(peak signal-to-noise ratio, PSNR)和结构相似度(structural similarity index, SSIM)(Wang等，2004)两个定量的图像评价指标(quantitative picture quality indices，PQI)来衡量复原结果的精度。其中, PSNR计算为 30 $P S N R=10 \lg \frac{L^{2}}{\frac{1}{n_{1} n_{2}}\left\|\tilde{\boldsymbol{X}}_{R}-\tilde{\boldsymbol{X}}_{T}\right\|}$ 式中, ${{n_1}, {n_2}}$ 为图像空间大小, ${{{\mathit{\boldsymbol{\tilde X}}}_R}}$ , ${{{\mathit{\boldsymbol{\tilde X}}}_T}}$ 分别为复原张量数据和原始张量数据; $L$ 为原始图像的最大像素值。PSNR用来衡量复原张量与原始张量之间的灰度相似性, 值越大, 说明复原张量数据与原始张量数据越接近, 复原效果越好。SSIM计算为 31 $S S I M=\frac{\left(2 \mu_{\tilde{\boldsymbol{X}}_{T}} \mu_{\tilde{\boldsymbol{X}}_{R}}+c_{1}\right)\left(2 \sigma_{\tilde{\boldsymbol{X}}_{T} \tilde{\boldsymbol{X}}_{R}}+c_{2}\right)}{\left(\mu_{\tilde{\boldsymbol{X}}_{T}}^{2}+\mu_{\tilde{\boldsymbol{X}}_{R}}^{2}+c_{1}\right)\left(\sigma_{\tilde{\boldsymbol{X}}_{T}}^{2}+\sigma_{\tilde{\boldsymbol{X}}_{R}}^{2}+c_{2}\right)}$ 式中, ${{\mu _{{{\mathit{\boldsymbol{\tilde X}}}_T}}}}$ 和 ${{\mu _{{{\mathit{\boldsymbol{\tilde X}}}_R}}}}$ 分别表示原始张量数据 ${{{\mathit{\boldsymbol{\tilde X}}}_T}}$ 和复原张量数据 ${{{\mathit{\boldsymbol{\tilde X}}}_R}}$ 的平均值。 ${{\sigma _{{{\mathit{\boldsymbol{\tilde X}}}_T}}}}$ 和 ${{\sigma _{{{\mathit{\boldsymbol{\tilde X}}}_R}}}}$ 分别表示原始张量数据 ${{{\mathit{\boldsymbol{\tilde X}}}_T}}$ 和复原张量数据 ${{{\mathit{\boldsymbol{\tilde X}}}_R}}$ 的方差, $c_1$ 和 $c_2$ 为非常小的正整数。SSIM的评价指标更接近于人眼的视觉感, 用于衡量复原张量数据和原始数据之间的结构相似性, 值越大, 说明复原张量数据质量越好, 结构性越完整。本文将复原HSI所有通道的PSNR和SSIM分别取平均值, 记为MPSNR和MSSIM, 作为复原结果的评价指标。实验环境为: CPU英特尔酷睿i7-8500(主频1.80 GHz)、内存8 GB、测试平台MATLAB R2016a。4.1HSI复原为了验证基于FWTNN的HSI复原模型的有效性, 本文选用了两种类型的HSI数据集进行实验, 第1个数据集是由高光谱数字图像收集实验(hyperspectral digital image collection experiment, HYDICE)传感器采集的Washington DC Mall (以下简称DC)数据集(http://www.ehu.es/ccwintco/index.php/ Hyperspectral_Remote_Sensing_Scenes)。整个空间大小为1 208×307像素, 波段数为191个光谱波段。实验中, 只选取了256×256像素大小的数据集。第2个数据集为通用分类像素相机(generalized assorted pixel camera, GAPC)传感器采集的Stuff数据集(http://www1.cs.columbia.edu/CAVE/databases/multispectral)。其空间大小为512×512像素, 波段范围为400~700 nm, 总计31个光谱波段。在实验中只选取了256×256像素的子块。对每个数据集 ${\mathit{\boldsymbol{\tilde O}} \in {{\bf{R}}^{{n_1} \times {n_2} \times {n_3}}}}$ , 定义采样比率(sampling rate, SR)为 32 $S R=\frac{|\tilde{\boldsymbol{\varOmega}}|}{\prod\limits_{i=1}^{3} n_{i}}$ 式中, $\left| {\mathit{\boldsymbol{ \boldsymbol{\tilde \varOmega} }}} \right|$ 表示索引集 $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$ 中的采样数, $\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}$ 是通过均匀分布从张量 ${\mathit{\boldsymbol{\tilde O}}}$ 中随机选择的。本文将在SR为10%, 20%和30%的情况下, 分别用基于HaLRTC(high accuracy low rank tensor completion)(Liu等，2013)、TNN(Lu等，2020)、PSTNN (Jiang等，2020)和FWTNN的HSI复原模型进行复原。表 1列出了3种采样率下所有复原方法的MPSNR、MSSIM的数值和算法平均运行时间。由表 1可知, 在不同采样率下, 基于FWTNN的HSI复原模型获得了最佳的PSNR及SSIM。图 8展示的是DC数据集在3种采样比率下, 4种复原方法结果的逐切片PSNR和SSIM曲线。从图 8中可以明显看出, 在不同的采样比率下, 本文方法在绝大部分正向切片矩阵中取得了最优的复原效果。图 9展示的是DC数据在采样比率SR=10%的条件下, 4种复原方法复原结果的第114个波段灰度图像及其PSNR。图 10展示的是Stuff数据集在采样比率SR=10%的条件下, 4种复原方法复原结果的第25个波段灰度图像及其PSNR。从直观的灰度结果中可以明显观察到, 相比于其他HSI复原方法, 本文方法能复原出更多的纹理细节, 取得最好的复原效果, 而相应的PSNR也给出数值例证。表1 在不同采样率(SR)下的定量评价结果和运行时间 SR/% 数据 指标 算法 HaLRTC TNN PSTNN FWTNN 10 DC MPSNR 19.785 4 32.219 3 32.551 0 33.536 8 MSSIM 0.423 6 0.887 9 0.888 3 0.907 6 Stuff MPSNR 21.309 3 25.917 6 26.410 2 28.259 9 MSSIM 0.413 2 0.689 6 0.713 9 0.799 1 20 DC MPSNR 27.792 7 37.127 1 37.453 3 39.094 8 MSSIM 0.763 1 0.956 6 0.958 4 0.968 6 Stuff MPSNR 24.136 3 30.774 0 31.631 1 34.287 1 MSSIM 0.598 4 0.873 7 0.881 1 0.923 6 30 DC MPSNR 33.738 1 40.221 3 40.495 1 42.263 2 MSSIM 0.914 0 0.976 5 0.977 3 0.983 1 Stuff MPSNR 26.133 9 33.898 8 34.995 3 37.905 4 MSSIM 0.730 0 0.928 6 0.933 1 0.956 2 DC 时间/s 109.10 590.55 601.51 647.36 Stuff 时间/s 20.19 89.78 84.49 98.53 Quantitative evaluation results and running time at different sampling ratesTable 1 加粗字体为每行最优值。 图8 DC数据复原结果的逐正向切片矩阵的PSNR与SSIM结果 The PSNR and SSIM of the forward slice matrix of the DC data restoration resultsFig 8((a)SR=10%; (b) SR=20%; (c) SR=30%) 图9 SR=10%时4种方法在DC数据的第114波段复原效果 The restoration effect of the four methods on the 114th band of DC data when SR=10%Fig 9((a) original image; (b) missing image; (c) HaLRTC; (d) TNN; (e) PSTNN; (f) FWTNN) 图10 SR=10%时4种方法在Stuff数据的第25波段复原效果 The restoration effect of the four methods on the 25th band of Stuff data when SR=10%Fig 10((a) original image; (b) missing image; (c) HaLRTC; (d) TNN; (e)PSTNN; (f)FWTNN) 4.2HSI去噪为了验证基于FWTNN的HSI去噪模型的有效性, 本文选用了两个高光谱遥感数据集进行实验。第1个数据集是反射光学系统成像光谱仪(reflective optics spectrographic imaging system, ROSIS-03)所拍摄的Pavia Centre数据集, 以下简称Pavia(http://www.ehu.es/ccwintco/index.php/Hyperspectral_Remote), 空间尺寸为1 096×1 096像素。由于该数据集中的一部分波段被噪声污染严重, 不能作为HSI去噪结果的参照, 故除去了这一部分污染严重的波段。在测试中, 只选取了空间尺寸为200×200像素、共80波段的数据进行实验。第2个数据集是机载可见/红外成像光谱仪(airborne visible infrared imaging spectrometer, AVIRIS)所拍摄的Indian数据集(https://engineering.purdue.edu/~biehl/MultiSpec/), 空间大小为145×145像素。在测试中, 选取前100个波段进行实验。针对两个选定的HSI数据集, 分别添加强度为0.5、0.4和0.3的稀疏噪声, 分别用基于TNN、PSTNN和FWTNN的HSI去噪模型进行复原。表 2列出了3种稀疏噪声强度下3种HSI去噪模型复原优劣的数值评判结果和算法平均运行时间。由表 2可知, 基于FWTNN的HSI去噪模型在绝大部分噪声情况下, 对HSI数据都达到最优的复原结果。而且在稀疏噪声强度最高的情况下, 本文模型能将MPSNR提高10 dB左右, 取得了非常好的去噪效果。表2 在不同的噪声情况下的定量评价结果和运行时间 稀疏度 数据 指标 算法 TNN PSTNN FWTNN 0.5 Pavia MPSNR 19.048 8 19.438 7 28.056 6 MSSIM 0.385 3 0.422 2 0.863 2 Indian MPSNR 21.957 6 23.758 8 34.529 9 MSSIM 0.471 4 0.543 6 0.946 5 0.4 Pavia MPSNR 28.486 6 29.829 9 36.502 8 MSSIM 0.844 0 0.874 1 0.959 3 Indian MPSNR 32.471 1 34.721 4 41.063 5 MSSIM 0.905 5 0.947 3 0.980 4 0.3 Pavia MPSNR 39.008 6 41.851 5 43.588 4 MSSIM 0.975 6 0.981 6 0.980 9 Indian MPSNR 37.323 6 42.762 9 45.661 2 MSSIM 0.983 6 0.987 5 0.988 2 Pavia 时间/s 59.11 97.42 100.27 Indian 时间/s 37.50 62.64 74.62 Quantitative evaluation results and running time under different noise conditionsTable 2 加粗字体为每行最优值。在3种强度稀疏噪声污染的情况下, 图 11显示的是Indian数据由3种HSI去噪模型复原结果所有波段的PSNR和SSIM曲线; 图 12显示的是Pavia数据由3种模型复原结果所有波段的PSNR和SSIM曲线。由图 11和图 12可知, 基于FWTNN的HSI去噪模型在复原结果的绝大多数波段都获得较好的数值结果。在噪声强度为0.4的条件下, 图 13第1行显示了所有模型对Pavia数据复原结果的第32个波段的灰度图像及其PSNR值; 图 14第1行显示了所有模型对Indian数据复原结果的第25个波段的灰度图像及其PSNR值。在噪声强度为0.5的条件下, 图 13第2行显示了所有模型对Pavia数据复原结果的第55个波段的灰度图像及其PSNR值; 图 14第2行显示了所有模型对Indian数据复原结果的第51个波段的灰度图像及其PSNR值。从上述直观复原结果的演示中可以明显观察到, 基于FWTNN的HSI去噪模型复原出的结果具有更多的纹理细节, 而且相应的数值评判结果PSNR值也取得最高。 图11 不同方法在Indian数据去噪结果的不同波段上PSNR与SSIM值 PSNR and SSIM values of different methods on different bands of Indian data denoising resultsFig 11((a) sparsity = 0.3; (b) sparsity = 0.4; (c) sparsity = 0.5) 图12 不同方法在Pavia数据去噪结果的不同波段上PSNR与SSIM值 PSNR and SSIM values of different methods on different bands of Pavia data denoising resultsFig 12((a) sparsity = 0.3; (b) sparsity = 0.4; (c) sparsity = 0.5) 图13 不同方法在Pavia数据的去噪结果 The denoising results of different methods on Pavia dataFig 13((a) original image; (b) noisy image; (c) TNN; (d) PSTNN; (e) FWTNN) 图14 不同方法在Indian数据的去噪结果 The denoising results of different methods on Indian dataFig 14((a) original image; (b) noisy image; (c) TNN; (d) PSTNN; (e) FWTNN) 5参数分析针对Indian数据, 分析了强度为0.5的稀疏噪声污染下基于FWTNN去噪模型参数的灵敏性。参数 $C_1$ 反映了频率域权重的缩放比例。固定其他参数, 图 15(a)(b)显示的是 $C_1$ 在[0.1, 4] 内取不同值时, 复原结果的PSNR和SSIM曲线。从图 15中可以看出, 当 $C_1$ 过小时, 频率域权重变小, 噪声更难去除干净; 当 $C_1$ 过大时, 频率域权重变大, 在抑制噪声的同时去除了张量数据的细节信息。此时, 最佳的参数值为 ${C_1}=1.6$ , 在拟制噪声的同时有效保护张量的细节结构。 图15 稀疏噪声强度为0.5时的Indian数据参数情况 Indian data parameters when the sparse noise intensity is 0.5Fig 15((a) PSNR with different $C_1$ values; (b) SSIM with different $C_1$ values; (c) PSNR with different $C_2$ values; (d)SSIM with different $C_2$ values) 参数 $C_2$ 反映的是迭代初始的频率域权重。将其他参数固定, 图 15(c)(d)显示的是 $ $ 在[0.1, 3] 取不同值时, 复原结果的PSNR和SSIM曲线。随着迭代次数的增加, $C_2$ 作用越来越小，不同的 $C_2$ 对模型的复原结果影响不敏感，此时最佳的参数值为 ${C_2}=1.0 $ 。为了让模型更具有普遍性和鲁棒性, 在最终实验中本文选取了两个让模型较为稳定的数值来作为所有去噪模型的参数, 分别为 ${C_1} = 1.5, {C_2} = 1$ 。正则化参数在所有去噪模型中均为 $1/\sqrt {\min \left({{n_1}, {n_2}} \right){n_3}} \;\;, \;{n_1}, {n_2}$ 为数据空间大小, $n_3$为波段数。图 16(a)—(c)是MPSNR、MSSIM和误差随迭代次数的变化情况。从图 16中可以看出, 本文算法是数值收敛的。 图16 稀疏噪声强度为0.5时Indian数据的迭代收敛情况 Iterative convergence of Indian data when the sparse noise intensity is 0.5Fig 16((a) MPSNR change with the number of iterations; (b) MSSIM change with the number of iterations; (c) error change with the number of iterations) 6结论针对高光谱图像复原问题, 本文在张量核范数定义的基础上, 考虑张量在傅里叶域下的物理意义, 即对清晰的高光谱图像添加噪声，图像信息在低频部分变化较小，而在高频部分变化巨大。基于这样的事实, 本文提出了频率加权的张量核范数来更好地近似张量秩。通过理论分析、实验验证, 给出了不同频率切片矩阵的频率权重的精确计算方法, 使得权重参数不仅与频率带携带的信息有关, 而且还可以根据缺失噪声/稀疏噪声的变化情况来自适应改变权重。在HSI复原和去噪两类实验中, 从定量的评价指标及定性的视觉效果两个方面验证了所提模型的有效性。虽然本文在光谱低秩上取得了更接近张量秩函数的结果, 但并未充分探索高光谱图像的空间平滑信息, 因此并不能有效地去除混合噪声。在后续研究中, 针对高光谱图像混合噪声的分布特点，将深入考虑高光谱的空间平滑性、组稀疏性和几何结构等先验，联合本文提出的低秩模型，提升模型的综合性能。此外, 本文模型也可以应用于诸如高光谱图像聚类、视频修复和MRI(magnetic resonance imaging)重构等高维数据处理任务。