论文笔记15：Hyperspectral Image Denoising via Noise-Adjusted Iterative Low-Rank Matrix Approximation

引言

HSI去噪方法：PCA、NAPCs、INAPCA、PCA with wavelet、multivariate multiresolution PCA、RPCA；PARAFAC、VBM3D、BM4D；基于PDE的方法(如SSAHTV)，基于小波分析的方法(如HSSNR)，小波+PCA、小波+多向Wiener滤波、小波+FORP等。

方法

HSI表示：$\mathbf{Y}=\mathbf{X}+\mathbf{N}$，其中$\mathbf{Y}=[{\mathbf{Y}}_1,{\mathbf{Y}}_2,\cdots ,{\mathbf{Y}}_p]\in {\mathbb{R}}^{MN\times p}$为Casorati矩阵(一种矩阵其列包含HSI的矢量化波段)，三维情形就是$u=f+\mathcal{N}$，$W=diag({\sigma }_1^2,{\sigma }_2^2,\cdots ,{\sigma }_p^2)$为噪声协方差矩阵，${\sigma }_i$为第$i$个波段的噪声标准差。

Patchwise LRMA Denoising

从线性光谱混合模型的角度来看，每个光谱特征(X的行(列?))可以由少量纯光谱端元的线性组合来表示。这启发我们利用LRMA方法来解决去噪问题$$\underset{\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{X}{\Vert }_F^2,\text{ s.t }\mathrm{rank}(\mathbf{X})\le r$$

由于$MN\gg p$，Y通常是一个非常薄(thin)的矩阵。在这种情况下，LRMA去噪后的图像可能会导致模糊和细节丢失。因此，我们以空间拼接的方式分析HSI，而不是全局分析。

设$R_b$是一个二进制算子，它从一个矩阵中提取m行n列，该矩阵对应于一个m×n的空间块(patch)，由每个HSI波段内的索引b指定，$R_b^{\ast }$为其逆过程。PLRMA去噪过程可定义为$${\mathrm{PLRMA}}_r(u)=\left(\sum _{b\in \Omega }{R}_b^{\ast }{\mathrm{LRMA}}_r\left(R_{b}u\right)\right)\cdot /\left(\sum _{b\in \Omega }{R}_b^{\ast }{R}_b\right)$$

其中$r$是所有Casorati矩阵的上界(下面详细讨论)，$\Omega$表示平铺(tiles)图像域的重叠块集合。上式对从HSI数据$u$中提取的矩阵族执行LRMA运算，并累计(accumulates)结果的加权和。

Noise-Adjusted Iteration Framework

迭代正则化的基本思想是将每次迭代的输出去噪图像加回到输入噪声图像中作为下一次迭代的输入$$u^{k+1}=(1-\delta )f^k+\delta u^k$$

其中$\delta \in [0,1]$为松弛(relaxation)参数，$u^k$表示第$k$次迭代的输入图像，$f^k$表示在$u^k$上进行PLRMA处理后的输出图像。由于不同的波段有不同的噪声水平，扩展了迭代正则化的思想$$u_i^{k+1}=(1-{\delta }_i)f_i^k+{\delta }_i{u}_i^{k}i=1,2,\cdots ,p.(1)$$

其中${\delta }_i=e^{-c(\mathbf{W}(i,i))}$，$\mathbf{W}(i,i)$是原始噪声图像的第$i$个波段的噪声方差，c是衰减参数，需要预先定义。

（这段比较难理解）随着迭代的开始，只有强信号(具有大的奇异值)能够经受住LRMA处理，并有助于初始的无噪声HSI $\hat{f}$。然而，部分恢复的信号将通过(1)反馈到输入信号，以降低噪声的估计。反过来，可以识别较弱的信号并将其添加到信号估计中。此外，当输入图像的第$i$个波段中的噪声水平较低时，第$i$个波段中较少恢复的信号将被加回到输入信号，反之亦然。这种策略有助于保护高信噪比波段中的较弱信号，这将反馈到低信噪比波段的性能改善。随着迭代的进行，通常观察到估计的噪声方差单调减小。同时，HSI结构逐渐恢复，直至收敛。

疑问：如果是第6步更新$f^{k+1}$那么第4步的$f^{k+1}$是怎么来的？是否要将第6步改成$f^{k+2}$?
其中
(5)：$u_i^{k+1}=(1-{\delta }_i)f_i^k+{\delta }_i{u}_i^{k}i=1,2,\cdots ,p$
(6)：${\delta }_i=e^{-c(\mathbf{W}(i,i))}$，$\mathbf{W}(i,i)$
(3)：${\mathrm{PLRMA}}_r(u)=\left({\sum }_{b\in \Omega }{R}_b^{\ast }{\mathrm{LRMA}}_r\left(R_{b}u\right)\right)\cdot /\left({\sum }_{b\in \Omega }{R}_b^{\ast }{R}_b\right)$

RSVD Algorithm for LRMA

主要的优化任务是求解${\min }_{\mathbf{X}}\Vert \mathbf{Y}-\mathbf{X}{\Vert }_F^2,\text{ s.t }\mathrm{rank}(\mathbf{X})\le r$，奇异值阈值法(SVT)将其转化为松弛(slack)问题：$$\underset{\mathbf{X}}{\min }\Vert \mathbf{Y}-\mathbf{X}{\Vert }_F^2+\lambda \Vert \mathbf{X}{\Vert }_{\ast }$$

其中$\Vert \mathbf{X}{\Vert }_{\ast }$为$\mathbf{X}$的所有奇异值之和，SURE提供了一种选择上式中正则化参数的有原则的自动方法。然而，SVT(SVD计算)的主要步骤是耗时的，包含$O(MN{p}^2)$浮点运算(flops)。由于每个波段的噪声水平不同，输入噪声图像Y应进行预白化(prewhitened)：$\tilde{\mathbf{Y}}=\mathbf{Y}{\mathbf{W}}^{-1/2}$. 但是，这种方法也会改变每个波段的幅度值(magnitude value)。

RSVD：
Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions

RSVD算法通过使用随机投影探索近似矩阵分解，将过程分为两个阶段。在第一阶段，随机抽样用于获得一个简化(reduced)矩阵，其范围近似于Y的范围，在第二阶段，简化的矩阵被分解。这个过程基于以下两个事实：

事实1：对于秩为$r$的矩阵$\mathbf{Y}\in {\mathbb{R}}^{MN\times p}$，$\mathbf{Y}$的范围(由$\mathbf{Y}$的列组成的子空间)可以定义为$\mathbf{\Omega }=\mathbf{Y}\mathbf{G}$，其中$\mathbf{G}\in {\mathbb{R}}^{p\times r}$为高斯随机矩阵。

事实2：如果$\mathbf{Q}\in {\mathbb{R}}^{MN\times r}$为一个具有跨越(span)$\mathbf{Y}$的范围(range)的正交列的矩阵，则
$$\Vert \mathbf{Y}-{\mathbf{Q}\mathbf{Q}}^{\ast }\mathbf{Y}\Vert \approx \underset{rank(\mathbf{X})\le r}{\min }\Vert \mathbf{Y}-\mathbf{X}\Vert$$

疑问：${\mathbf{Q}}^{\ast }$是什么？共轭转置？

从算法2可以看到，RSVD的计算包含$O(MNplog(r)+(MN+p)r^2)$flops，比SVD效率高。

在算法2中，如果$Y$的奇异值逐渐衰减，估计的精度可能会损失。引入一个修正的RSVD，计算$\mathbf{H}=({\mathbf{Y}\mathbf{Y}}^{\ast }{)}^q\mathbf{Y}\mathbf{G}$而不是$\mathbf{H}=\mathbf{Y}\mathbf{G}$，因为$({\mathbf{Y}\mathbf{Y}}^{\ast }{)}^q\mathbf{Y}\mathbf{G}$具有与$\mathbf{Y}\mathbf{G}$相同的奇异向量，而其奇异值衰减得更快。

Adaptive Determination of the Parameters

我们采用了基于多元回归理论的方法来估计噪声，(以下省略加粗)$Y=[Y_1,Y_2,\cdots ,Y_p]$，令$${\mathbf{Y}}_{{\partial }_i}=\left[{\mathbf{Y}}_1,{\mathbf{Y}}_2,\dots ,{\mathbf{Y}}_{i-1},{\mathbf{Y}}_{i+1},\dots ,{\mathbf{Y}}_p\right]$$

为$MN\times (p-1)$矩阵，我们假设$Y_i$可由剩余的$p-1$个波段线性组合表示$$Y_i=Y_{{\partial }_i}{\beta }_i+{\xi }_i$$其中${\beta }_i$是回归向量，大小为$(p-1)\times 1$，${\xi }_i$为噪声向量，大小为$MN\times 1$，由最小二乘估计可知$$\hat{\beta }={\left({\mathbf{Y}}_{{\partial }_i}^T{\mathbf{Y}}_{{\partial }_i}\right)}^{-1}{\mathbf{Y}}_{{\partial }_i}^T{\mathbf{Y}}_i$$

因此噪声估计为${\hat{\xi }}_i={\mathbf{Y}}_i-{\mathbf{Y}}_{{\partial }_i}\hat{\beta }$，相关矩阵为$\mathbf{W}=[{\hat{\xi }}_1,\dots ,{\hat{\xi }}_p{]}^T[{\hat{\xi }}_1,\dots ,{\hat{\xi }}_p]/(MN)$.

我们提供一种简单的基于奇异值分解的上界秩估计方法。对$\mathbf{Y}$和$\mathbf{N}=[{\hat{\xi }}_1,\dots ,{\hat{\xi }}_p]$使用奇异值分解，得到奇异值分别为$[{\sigma }_1,\cdots ,{\sigma }_p],[S_1,\cdots ,S_p]$，现在可以求出参数值$r$，使得${\sigma }_r\ge S_1$，${\sigma }_{r+1}<S_1$。秩的上界设置为$r$，意味着第一个$r$维子空间包含更多的信号。估计的上界秩为7，去噪的下一步是从每个块的7维信号子空间中分离噪声。

NAILRMA算法中只需要指定三个参数：块大小、步长和衰减参数$c$。块的大小设置为$20\times 20\times p$(与LRMR同步)，步长设置为8×8，即我们只在水平和垂直方向上每8个像素选择一个块，以加快计算速度。每个不同波段的噪声变化率不变，我们简单地将衰减参数$c$设置为5(通过实验证明)，忽略绝对值，但集中于迭代参数$\delta$的相对值。

推广到混合噪声去除

噪声退化模型为$Y=X+S+N$，$S$为稀疏噪声，为脉冲噪声、死像素或死线以及条纹的混合。低秩近似去噪然后被扩展为众所周知的RPCA模型：$$\underset{\mathbf{X},\mathbf{S}}{\min }\Vert \mathbf{Y}-\mathbf{X}-\mathbf{S}{\Vert }_F^2\text{ s.t. }\mathrm{rank}(\mathbf{X})\le r,\mathrm{card}(\mathbf{S})\le b$$

这里，$b$代表稀疏噪声的上限基数。LRMR算法利用RPCA模型对HSI逐块去噪，是最先进的HSI混合噪声去除方法之一。为了进一步提高算法的密集去噪性能，还采用了算法1中提出的噪声调整迭代框架，并结合了LRMR算法(NAILRMR)对HSI进行去噪。唯一的修改是利用RPCA模型代替算法1的步骤6中的(3)。

实验

对比方法：SSAHTV、VBM3D、BM4D、SURE-SVT、LRMR

模拟实验数据：Washington DC Mall(256×256×191)，Pavia city center(200×200×80).

情况1：零均值，噪声方差分别为0.02、0.04、0.06、0.08和0.1.

情况2：不同的方差零均值高斯噪声被添加到两个HSI数据集的每个波段。方差值从0到0.1随机选取，平均方差值分别为0.053和0.051.

SSAHTV和LRMR的参数选择：

为了进一步比较所有去噪算法的性能，我们计算了无噪声和去噪后的HSI之间所有光谱特征的平均谱角距离(越小越好)：
$$\mathrm{M}\mathrm{S}\mathrm{A}\mathrm{D}=\frac{1}{MN}\sum _{i=1}^{MN}\frac{180}{\pi }\times \arccos \frac{{\left({\mathbf{X}}^i\right)}^T\cdot \left({\hat{\mathbf{X}}}^i\right)}{\Vert {\mathbf{X}}^i\Vert \cdot \Vert {\hat{\mathbf{X}}}^i\Vert }$$

混合噪声实验参数：r=7，b=6000，步长=8.

在所有模拟和真实数据实验中，停止标准设定为$\epsilon =1e-3$，最大迭代次数为50.

具体实验结果见原论文，这里不再给出。

代码