ZSL (TPAMI 2023): 基于Zero-Shot方法的光谱锐化模型

Zero-Shot Hyperspectral Sharpening

论文链接：Zero-Shot Hyperspectral Sharpening
代码：

• ZSL_official
• ZSL_format

摘要

​ 近年来CNN在光谱图像融合锐化效果出色，然而这些方法往往缺乏训练数据，泛化能力有限，于是提出了一种用于光谱图象的Zero-Shot学习方法。

​ 具体而言，我们首先提出了一种高精度定量估计成像传感器空间和光谱响应的新方法。在训练过程中，我们根据估计的空间响应对MSI和HSI进行空间抽样，并使用下采样的HSI和MSI来推断原始HSI。这样，我们不仅可以利用HSI和MSI的固有信息，而且训练好的CNN也可以很好地泛化到测试数据中。

​ 此外，我们在HSI上进行了降维，在不牺牲融合精度的情况下减少了模型尺寸和存储使用量。此外，我们为CNN设计了一个基于成像模型的损失函数，进一步提高了CNN的学习性能。

现阶段问题

问题：缺乏训练数据和泛化能力有限

1. 图像数据少，难以生成：在实际应用中，难以获取同一场景的LR-HSI,HR-MSI和HR-HSI。为了解决这个训练数据少的问题，在图像去噪领域，Imamura提出对原始噪声图像进行加噪声生成训练数据。但是在全色锐化领域，无法简单通过这种方式进行数据生成。
2. 成像模型不同导致CNN泛化性差：当训练数据和测试数据的光谱范围、光谱波段数、成像模型不同时，模型表现会受到影响。为了解决这个问题，Xie通过使用各种常见的成像方式生成训练数据，以尽可能的减小成像不同产生影响。

解决思路：

​ 现阶段生成训练数据，都是假设我们预先知道了光谱响应函数(SRF)和传感器的点扩散函数 (PSF)，但是在实际情况中，这两个参数都是难以预先知道的。因此，本文提出一种高精度定量估计光谱和空间响应的新方法。通过从测试数据预先计算SRF和PSF，然后将这两个参数用于生成训练数据，这样就在一定程度保证了训练数据和测试数据在成像模型上保持一致。

主要贡献

1. 提出使用HSI锐化的Zero-Shot学习方法，一种高精度定量估计成像传感器空间和光谱响应的方法，用于Zero-Shot学习方法生成数据。
2. 对HSI进行降维，提高模型的效率
3. 基于子空间表示，设计一个基于成像模型的损失函数，由训练数据损失和成像模型损失组成。

HSI锐化成像模型

​ 常见的将HR-HSI图像变为LR-HSI以及HR-MSI有如下关系，其中B为模糊核(表示空间响应，取决于传感器的点扩散函数 (PSF))，D为下采样因子，R为光谱响应函数(SRF)
$$X_{lr}=(X\ast B)D+N_x,\ast 为卷积$$

$$Y=RX+N_y$$

$$(\mathbf{Y}\ast \mathbf{B})\mathbf{D}={\mathbf{RX}}_{lr}+\mathbf{N},$$

model-based方法

​ 通过(1),(2)，我们可以通过求如下优化问题：
$$\underset{\mathbf{X}}{\min }||{\mathbf{X}}_{lr}-(\mathbf{X}\ast \mathbf{B})\mathbf{D}|{|}_F^2+||\mathbf{Y}-\mathbf{RX}|{|}_F^2+\varphi (\mathbf{X}),$$
​ 其中$\varphi (X)$表示HR-HSI的先验信息，从上述优化问题我们也知道，必须要提前知道$B$,$R$才能进行求解，然而$B,R$一般都是未知的，只能进行定量估计。

cnn-based方法

​ 通常，基于CNN的方法需要同一场景的LR-HSI,HR-MSI和对应的HR-HSI进行训练，然而这种训练数据是很难获取的，一般都是采用data simulation进行空间和光谱退化生成同一场景的HSI和MSI。

​ 然而，这种方式忽略了模拟中的空间和光谱退化算子应该与测试数据一致，以达到良好性能，这种方法也通常是假设B和R在模拟实验中是完全已知，但在实际中是不成立的，因此，对B,R的定量估计也是只管重要的。

求解过程

由式子1，2，可以得到：
$$(Y\ast B)D=R{X}_{lr}+N$$
优化问题可以变为：
$$\underset{B,R}{\min }||(Y\ast B)D-R{X}_{lr}||$$
由于（6）是非凸的，解不唯一，因此利用先验知识规范化这个优化为题，对就模糊核B进行归一化，并让B,R满足非负约束。
$$\begin{array}{c}\underset{\mathbf{B},\mathbf{R}}{\min }||(\mathbf{Y}\ast \mathbf{B})\mathbf{D}-\mathbf{R}{\mathbf{X}}_{lr}|{|}_F^2\\ s.t.\sum _{i,j}{\mathbf{B}}_{i,j}=1,\forall i,j,{\mathbf{B}}_{i,j}\ge 0,{\mathbf{R}}_{i,j}\ge 0.\end{array}$$
优化问题相对于B,R是凸的，采用交替最小化方案(alternating minimization scheme)
$${\mathbf{R}}^{(k)}=\underset{\mathbf{R}}{\mathrm{argmin}}||{\mathbf{H}}^{(k)}-\mathbf{R}{\mathbf{X}}_{lr}|{|}_F^2,s.t.\forall i,j,{\mathbf{R}}_{i,j}\ge 0,$$

$$\begin{array}{rl}{\mathbf{B}}^{(k)}=\underset{\mathbf{B}}{\mathrm{argmin}}\Vert & (\mathbf{Y}\ast \mathbf{B})\mathbf{D}-{\mathbf{E}}^{(k)}{\Vert }_F^2,\\ & s.t.\sum _{i,j}{\mathbf{B}}_{i,j}=1\mathrm{and}\forall i,j,{\mathbf{B}}_{i,j}\ge 0,\end{array}$$

其中$H^k=(Y\ast B^{k-1})D,E^k=R^k{X}_{lr}$

更新R

$$\begin{array}{rl}L(\mathbf{R},\mathbf{K},\mathbf{L})& =||{\mathbf{H}}^{(k)}-\mathbf{R}{\mathbf{X}}_{lr}|{|}_F^2\\ & +\mu {\Vert \mathbf{R}-\mathbf{K}+\frac{\mathbf{L}}{2\mu }\Vert }_F^2+l(\mathbf{K}),\end{array}$$

$L$是拉格朗日算子， $l(\cdot )$ 是指示函数，当$K$的所有元素大于0才为0。

$$l(\mathbf{K})=\{\begin{array}{cc}0& \text{ if }\forall i,j,{\mathbf{K}}_{ij}\ge 0,\\ +\infty & \text{ otherwise}\end{array}$$

方程(10)通过交替最小化更新 $R,K,L$.$R,K$根据如下方程更新
$$\begin{array}{rl}& \mathbf{R}=\underset{\mathbf{R}}{\mathrm{argmin}}L(\mathbf{R},\mathbf{K},\mathbf{L}),\\ & \\ & \mathbf{K}=\underset{\mathbf{K}}{\mathrm{argmin}}L(\mathbf{R},\mathbf{K},\mathbf{L})\end{array}$$

(12)的解如下：
$$\begin{gathered} R=\left({\mathrm{H}}^{(k)}{\mathbf{X}}_{lr}^T+\mu \mathbf{K}-\frac{\mathbf{L}}{2\mu }\right)({\mathbf{X}}_{lr}{\mathbf{X}}_l^T+\mu \mathbf{I}{)}^{-1} \\ K=(R+\frac{\mathbf{L}}{2\mu }{)}_{+} \end{gathered}$$

$(\cdot {)}_{+}$是将K的非0元素置0，$I$为单位矩阵

拉格朗日算子$L$通过如下方式更新
$$\mathbf{L}=\mathbf{L}+2\mu (\mathbf{R}-\mathbf{K})$$

更新B

​ $(Y\ast B)D$表示模糊核B和步长的图像Y各波段的二维卷积,因此可以先将$Y$进行分解，并使用stride为s进行下采样，优化问题变为：
$$\underset{\mathbf{b}}{\min }\sum _{i,j}||{\mathbf{y}}_{ij}^T\mathbf{b}-{\mathbf{E}}_{ij}^{(k)}|{|}_F^2,s.t.\sum _i{\mathbf{b}}_i=1\mathrm{and}\forall i,{\mathbf{b}}_i\ge 0$$
令$v=b$，同时引入了拉格朗日算子g
$$\begin{array}{rl}E(\mathbf{b},\mathbf{v},\mathbf{g})& =\sum _{i,j}||{\mathbf{y}}_{ij}^T\mathbf{b}-{\mathbf{E}}_{ij}^{(k)}|{|}_F^2\\ & +\mu ||\mathbf{b}-\mathbf{v}+\frac{\mathbf{g}}{2\mu }|{|}_F^2+h(\mathbf{v})\end{array}$$

指示函数如下:

$$h(\mathbf{v})=\{\begin{array}{cc}0& \text{ if }\forall i,{\mathbf{v}}_i\ge 0\text{ and }{\sum }_i{\mathbf{v}}_i=1\\ +\infty & \text{ otherwise}\end{array}$$

$$\begin{gathered} b=\arg \underset{b}{\min }E\left(\mathrm{b},\mathrm{v},\mathrm{g}\right) \\ ={\left(\sum _{i,j}{\mathbf{y}}_{ij}{\mathbf{y}}_{ij}^T+\mu \mathbf{I}\right)}^{-1}\left(\sum _{i,j}{\mathbf{E}}_{ij}^{(k)}{\mathbf{y}}_{ij}+\mu \mathbf{v}-\frac{\mathbf{g}}{2}\right) \end{gathered}$$

$$\begin{array}{rl}\mathbf{v}& =\underset{\mathbf{v}}{\mathrm{argmin}}E(\mathbf{b},\mathbf{v},\mathbf{g})\\ & \\ & =\underset{\mathbf{v}}{\mathrm{argmin}}{\left|\left|\mathbf{b}-\mathbf{v}+\frac{\mathbf{g}}{2\mu }\right|\right|}_F^2+h(\mathbf{v}),\end{array}$$

问题（19）下通过下述引理更新求解（对v的更新）

[63]为L. Condat, “Fast projection onto the simplex and the 1 ball,” Math. Program., vol. 158, no. 1/2, pp. 575–585, 2016.

拉格朗日算子更新：
$$g=g+2\mu (b-v)$$

迭代收敛判断

$$\frac{\Vert {\mathbf{B}}^{(k+1)}-{\mathbf{B}}^{(k)}{\Vert }_F}{\Vert {\mathbf{B}}^{(k)}{\Vert }_F}+\frac{\Vert {\mathbf{R}}^{(k+1)}-{\mathbf{R}}^{(k)}{\Vert }_F}{\Vert {\mathbf{R}}^{(k)}{\Vert }_F}\le tol$$

$$\begin{gathered} f(R,B)=Q(R,B)+l(R)+h(B) \\ Q(R,B)=||(Y\ast B)D-R{X}_{lr}|{|}_F^2 \end{gathered}$$

命题 1：

​ 假设方阵 $X_{lr}{X}_{lr}^T$ 和 ${\sum }_{ij}{y}_{ij}{y}_{ij}^T$ 具有满秩，并且 $(B^{(k)},R^{(k)}{)}_{k\in \mathbb{N}}$，即由算法 1 生成的序列是有界的。在这种情况下，序列将收敛到函数 $f$ 的临界点。

证明：

​ 首先，函数 $Q$ 是可微的，并且其梯度具有 Lipschitz 连续性（由于 $(B^{(k)},R^{(k)}{)}_{k\in \mathbb{N}}$的有界性）。进一步假设方阵 $X_{lr}{X}_{lr}^T$和 ${\sum }_{ij}{y}_{ij}{y}_{ij}^T$ 具有满秩，这表明函数 $Q$ 对于 $B$ 和 $R$ 分别是强凸的。除此之外，$l(R)$ 和 $h(B)$ 是凸集的指示函数，而 $f:{\mathbb{R}}^{sS}\times {\mathbb{R}}^m⟶\mathbb{R}$ 满足 Kurdyka-Łojasiewicz（KL）不等式（详见引用 [64] 的2.2Kurdyka-Lojasiewicz inequality）并且具有下界。

基于这些假设和性质，通过直接应用引用 [64] 中的定理 2.8，可以证明序列收敛性。

在这里，假设Algorithm 1是 [64]中Algorithm 1的实例，使用（1.3a）进行更新。

在实际的 HSI 和 MSI 数据中，方阵 $X_{lr}{X}_{lr}^T$和 ${\sum }_{ij}{y}_{ij}{y}_{ij}^T$ 通常具有满秩，这符合上述假设条件。

（1.3a）

子空间表示降维

​ HSI具有较高的光谱分辨率和许多光谱波段，显著提高了CNN的计算负担和卷积滤波器的数量。为了减少计算负担和参数，我们建议将HSI分解为低维子空间和系数。

​ 利用光谱波段相似度引起的HSI的低秩结构，可以很好地逼近许多真实的HSIs
$$X=SA$$

​ 上式具体参见Non-local Meets Global: An Iterative Paradigm for Hyperspectral Image Restoration

​ 因为其低秩结构，采用SVD奇异值分解估计光谱子空间
$$\begin{gathered} Y=U\Sigma V^T \\ S=U(:,1:L) \end{gathered}$$

子空间是半酉的，$S{S}^T=I$

$$\begin{gathered} A_{lr}=S^T{X}_{lr}, \\ {A}_d=S^T{X}_d, \\ {A}_d=(A_{lr}\ast B)\hat{D} \end{gathered}$$

损失函数

$$\begin{array}{rl}loss& =los{s}_g+\lambda los{s}_m\\ & =\frac{1}{M}||{\mathbf{X}}_{lr}-\mathbf{S}{\mathbf{A}}_{lr}|{|}_1+\frac{\lambda }{P}||(\mathbf{S}{\mathbf{A}}_{lr}\ast \mathbf{B})\mathbf{D}-{\mathbf{X}}_d|{|}_1\\ & +\frac{\lambda }{Q}||{\mathrm{RSA}}_{lr}-{\mathrm{Y}}_d|{|}_1\end{array}$$

其中

• $los{s}_g$为训练数据损失，估计的高分辨率 HSI 和地面实况之间的差异损失。

• $los{s}_m$为成像损失，也就是通过计算得到的相应矩阵B,R对图像进行转换时的损失。

训练部分

代码：采样$X_{lr}$

[HRHSI, R] = dataset_input(data2, downsample_factor)

"""计算X_lr"""
# (X*B)D
HSI0 = Gaussian_downsample(HRHSI, PSF, downsample_factor)
# (X*B)D+N_x=X_lr
HSI0 = HSI0 + 0 * np.random.randn(HSI0.shape[0], HSI0.shape[1], HSI0.shape[2])

"""计算y，用来给model提供channel数据"""
MSI0 = np.tensordot(R, HRHSI, axes=([1], [0]))

"""子空间表示"""
HSI3 = HSI0.reshape(HSI0.shape[0], -1)
# 奇异值分解用来降维
U0, S, V = np.linalg.svd(np.dot(HSI3, HSI3.T))
# p为论文中的L，该奇异值分解保留了L个最大奇异值
U0 = U0[:, 0:int(p)]
# 论文中的A_lr
HSI0_Abun = np.tensordot(U0.T, HSI0, axes=([1], [0]))

上方是针对数据模拟的，假设我们先已知了psf和R，然后从HR-HSI图像先去获取LR-HSI，HR-MSI

如果是实际情况，我们先获取到了HR-MSI和LR-HSI，则只需要以下

"""
模拟数据生成，假设已知psf和R，使用dataset_input
对于真实情况，如果是已知HR-MSI和LR-HSI，可以通过调用Kernal_estimation.m近似求出R和psf
然后将R，psf代入到下面的式子中直接进入子空间表示，降维部分，HSI0=LR-HSI

接着对HSI0进行数据增强，然后在生成数据过程中，采样得到X_d,对于Y_d则是通过响应函数直接计算得到。
"""

#HSI0=LR-HSI
#MSI:用来提供model初始化输入channel的数据
#psf，R=读取mat得到


"""子空间表示，降维"""
HSI3 = HSI0.reshape(HSI0.shape[0], -1)
# 奇异值分解用来降维
U0, S, V = np.linalg.svd(np.dot(HSI3, HSI3.T))
# p为论文中的L，该奇异值分解保留了L个最大奇异值
U0 = U0[:, 0:int(p)]
# 论文中的A_lr
HSI0_Abun = np.tensordot(U0.T, HSI0, axes=([1], [0]))

数据生成

先对$X_{lr}$进行augmentation，此处已经使用降维后的数据代替HSI图像，Abun 表示降维后的。

augument = [0]
HSI_aug = []
HSI_aug.append(HSI0)
U = U0
train_hrhs = []
train_hrms = []
train_lrhs = []
for j in augument:
    # 垂直翻转
    HSI = cv2.flip(HSI0, j)
    HSI_aug.append(HSI)

for j in range(len(HSI_aug)):
    HSI = HSI_aug[j]
    #        MSI = MSI_aug[j]
    """计算A_d"""
    HSI_Abun = np.tensordot(U.T, HSI, axes=([1], [0]))
    HSI_LR_Abun = Gaussian_downsample(HSI_Abun, PSF, downsample_factor)

    """计算Y_d"""
    MSI_LR = np.tensordot(R, HSI, axes=([1], [0]))


    """通过stride移位置遍历图像每一个区域，划分得到training_size尺寸的数值"""
    for j in range(0, HSI_Abun.shape[1] - training_size + 1, stride):
        for k in range(0, HSI_Abun.shape[2] - training_size + 1, stride):
            temp_hrhs = HSI[:, j:j + training_size, k:k + training_size]
            temp_hrms = MSI_LR[:, j:j + training_size, k:k + training_size]
            temp_lrhs = HSI_LR_Abun[:, int(j / downsample_factor):int((j + training_size) / downsample_factor), int(k / downsample_factor):int((k + training_size) / downsample_factor)]

            train_hrhs.append(temp_hrhs)
            train_hrms.append(temp_hrms)
            train_lrhs.append(temp_lrhs)

#用于当标签            
train_hrhs = torch.Tensor(np.array(train_hrhs))
#LR-LSI  
train_lrhs = torch.Tensor(np.array(train_lrhs))
#HR-MSI
train_hrms = torch.Tensor(np.array(train_hrms))

matlab代码

​ demo_kernel_estimation.m：假设已知光谱相应R和通过fspecial求解出初始的psf(B),通过使用R,B对HR-HSI进行转换，得到LR-HSI和HR-MSI。再将LR-HSI和HR-MSI通过kernel_estimation.m进行求解得到R,B

​ 实际应用中，我们一般知道LR-HSI和HR-MSI，然后通过这两个输入到kernel_estimation.m就能求解出相应R,B。然后通过R,B去对图像进行采样操作得到训练数据。

结论

本文介绍了一种新的基于 CNN 的 HSI 锐化方法。与现有方法不同的是，我们的方法基于HSI锐化的成像模型生成训练数据。具体来说，我们首先提出了一种新的方法来准确估计观测到的HSI和MSI的光谱和空间响应。然后，我们利用估计的响应来生成训练数据，可以利用HSI和MSI的固有信息并训练光图像自适应CNN。此外，提出了一种新的损失函数，它由训练数据和成像模型组成，以提高融合性能。此外，我们的方法在CNN的输入和输出中采用降维技术，即可以大大减少CNN的参数，提高计算效率。实验证明，我们的方法可以在不牺牲锐化图像的质量的情况下大大减少计算时间。在未来的工作中，我们还可以考虑利用 CNN 中的估计响应，这可能会进一步提高融合精度。