论文阅读——Variational Zero-Shot Multispectral Pansharpening

Abstract

全色锐化的目标是通过融合低空间分辨率的多光谱图像（LRMS）和全色图像（PAN）来生成高空间分辨率的多光谱图像（HRMS）。这项任务最具挑战性的问题在于，只有待融合的LRMS和PAN图像是可用的，而现有的基于深度学习（DL）的方法由于依赖大量训练对而不适合解决此类问题。传统的基于变分优化（VO）的方法非常适合解决此类问题。这些方法侧重于为优化问题精心设计显式的融合规则和正则化项，这些规则和正则化项基于研究人员对图像关系和图像结构的发现。与以往基于VO的方法不同，在本工作中，我们通过参数化项而非手动设计项来探索这种复杂的关系。具体而言，我们提出了一种零样本全色锐化方法，通过在优化目标中引入神经网络来实现。该网络估计HRMS的一个表示分量，主要用于描述HRMS与PAN之间的关系。通过这种方式，网络实现了与所谓的深度图像先验（DIP）类似的目标，因为它通过其固有结构隐式地规范了HRMS和PAN图像之间的关系。我们通过交替最小化直接优化网络参数和期望的HRMS图像来最小化这一优化目标。在各种基准数据集上的大量实验表明，与其他最先进的（SOTA）方法相比，我们提出的方法能够实现更好的性能。

Main Method

我们使用大写字母表示矩阵，例如 $A\in {\mathbb{R}}^{H\times W}$。具有超过两个维度的张量用书法体表示，例如 $\mathcal{X}\in {\mathbb{R}}^{H\times W\times S}$。$\mathcal{X}$ 的第 ((i, j, k)) 个元素表示为 ${\mathcal{X}}_{ijk}$。符号“$\otimes$”表示 $A$ 和 $B$ 之间的卷积运算。符号“$\odot$”表示哈达玛积（Hadamard product），即逐元素相乘。相应地，“$\oslash$”表示逐元素除法。“$A{\downarrow }_r$”表示将 $A$ 下采样，尺度因子为 $r$。矩阵或张量的弗罗贝尼乌斯范数（Frobenius norm）写作“$\Vert \cdot {\Vert }_F$”，即 $\Vert \mathcal{X}{\Vert }_F:={\left({\sum }_{ijk}{\mathcal{X}}_{ijk}^2\right)}^{1/2}$。

Overview to Pansharpening

多光谱全色锐化的目标是将观测到的低分辨率多光谱图像（LRMS）和全色图像（PAN）融合为期望的高分辨率多光谱图像（HRMS）。设 $X\in {\mathbb{R}}^{H\times W\times S}$ 表示HRMS，其中 $H$、$W$ 和 $S$ 分别表示HRMS的高度、宽度和光谱维度。对应的LRMS表示为 $Y\in {\mathbb{R}}^{h\times w\times S}$，全色图像表示为 $P\in {\mathbb{R}}^{H\times W}$。LRMS的下采样因子为 $r$，即 $H/h=W/w=r$。将多光谱全色锐化方法视为逆问题[69]，我们可以为基于变分优化（VO）的通用方法[5]建立以下模型：

$$\underset{X}{\min }{L}_y(X,Y)+{\lambda }_1{L}_p(X,P)+{\lambda }_{2}R(X)(1)$$

其中，$L_y(\cdot ,\cdot )$ 和 $L_p(\cdot ,\cdot )$ 分别表示描述HRMS/LRMS和HRMS/PAN关系的数据保真项，$R(X)$ 表示用于表征 $X$ 结构的正则化项，${\lambda }_1$ 和 ${\lambda }_2$ 是两个权衡参数，用于平衡各项的比例。

在基于VO的全色锐化方法中，最常用的HRMS与LRMS之间的关系如下[71]：

$$Y=(X\otimes K){\downarrow }_r+n_1(2)$$

其中，$K\in {\mathbb{R}}^{k\times k}$ 表示模糊核，$n_1$ 表示小的残差值，通常建模为零均值高斯分布。简而言之，上述表示意味着LRMS被认为是HRMS经过模糊和下采样后的结果。因此，$L_y$ 可以表示为：

$$L_y=\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2(3)$$

根据以往的工作，我们还采用与多光谱传感器的调制传递函数（MTF）匹配的高斯滤波器作为核 $K$ [72]。该核是预设且固定的。可以看出，$L_y$ 项主要衡量空间信息损失，因为模糊和下采样是在空间维度上进行的。

与 $L_y$ 不同，$L_p$ 项的设计更为复杂，因为PAN如何从HRMS退化仍然是一个开放性问题。在早期的工作中[38]，提出了简单的线性模型 $P=X{\times }_{3}p$，随后被许多其他工作采用[39]-[43]。这意味着对 $X$ 的光谱波段进行线性组合，系数向量为 $p\in {\mathbb{R}}^{S\times 1}$。线性HRMS/PAN模型的最大优势之一是在求解优化问题时计算方便，因为它仅包含矩阵乘法操作。此外，由于操作是在光谱维度上进行的，空间信息可以很好地包含在内。然而，线性模型在提高模型性能方面存在一些缺点。首先，即使 $P=X{\times }_{3}p$ 严格成立，由于 $\text{rank}(p)=1$，该方程的解空间仍然很大，这意味着线性模型无法提供更多的图像信息。其次，如果 $P=X{\times }_{3}p$ 的关系不严格成立，这种线性模型不可避免地会带来近似误差。

为了进一步提高模型性能，研究了更复杂的模型来描述HRMS/PAN关系，以增强PAN的引导能力[15]-[25]。例如，一些工作重新采用了基于CS（成分替换）和MRA（多分辨率分析）的方法提出的表示[23]-[25]。其中，最简单的模型将HRMS表示为：

$$X=G\odot \hat{P}(4)$$

公式（4）具有非常简洁和灵活的形式，类似于Brovey变换[28]。该表示仅包含两个部分。扩展的全色图像 $\hat{P}\in {\mathbb{R}}^{H\times W\times S}$ 是从PAN构建的，并且具有与HRMS相同的大小，尤其是空间分辨率。例如，$\hat{P}$ 可以与PAN和LRMS进行直方图匹配。可以看出，在这种表示中，$\hat{P}$ 主要为HRMS提供空间信息。然而，由于PAN只有一个通道，$\hat{P}$ 无法保留与 $X$ 一致的光谱信息。这就是为什么需要一个系数张量 $G$ 来平衡近似。理论上，无论 $X$ 和 $P$ 之间的真实关系如何，公式（4）始终成立。根据公式（4），$L_p$ 项可以表示为：

$$L_p=\Vert X-G\odot \hat{P}{\Vert }_F^2(5)$$

让我们更仔细地分析公式（5）。假设 $G$ 完全未确定，$L_p$ 显然不能用于问题（1），因为无论 $X$ 是什么，我们都可以设置 $G=X\oslash \hat{P}$，使得 $L_p$ 达到最小值，即零。换句话说，如果没有额外的约束，$X=G\odot \hat{P}$ 实际上并未揭示关于 $X$ 的有效图像信息。在以往的工作中，无论是基于CS的方法[28]-[32]、基于MRA的方法[34]-[37]，还是相关的基于VO的方法[23]-[25]，$G$ 都应该是预设的，并且提出了各种巧妙的方法。当 $G$ 通过某种方法确定时，我们可以看到 $G\odot \hat{P}$ 是对 $X$ 的直接近似，而 $L_p$ 直接约束 $X$。

Proposed Model

与现有工作主要关注如何预设 $G$ 不同，我们将公式（4）中的 $G$ 视为一个潜在变量，并在多光谱全色锐化问题（1）中与 $X$ 一起优化。具体来说，我们将公式（1）重新表述为：

$$\underset{X,G}{\min }\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+{\lambda }_1\Vert X-G\odot \hat{P}{\Vert }_F^2+{\lambda }_{2}R(G)(6)$$

在提出的模型（6）中，添加了一个关于 $G$ 的正则化项。如上所述，如果正则化直接针对 $X$，则无论 $X$ 是什么，$L_p=\Vert X-G\odot \hat{P}{\Vert }_F^2$ 始终可以达到零。因此，包含PAN信息的 $\hat{P}$ 对问题的解决没有有效贡献。然而，对 $G$ 的正则化使得空间保真项 $L_2$ 影响问题的解。通过这种方式，问题（6）不会高度不适定，并且期望的HRMS会更合理地受到LRMS和PAN的共同引导。

此外，我们在模型（6）中简单地移除了 $R(X)$。需要注意的是，由于 $G\odot \hat{P}$ 是对 $X$ 的近似，对 $G$ 的正则化可以看作是对 $X$ 结构的间接约束。为了更清楚地理解这一点，我们可以将 $R(G)$ 近似为 $R(X\oslash \hat{P})$。此外，我们有以下定理，表明 $R(X)$ 实际上可以被“吸收”到 $R(G)$ 中。

定理1：对于任何正则化 $R_g(G)$ 和 $R_x(X)$，如果 $(X^{\ast },G_1)$ 是以下问题的最小点：

$$\underset{X,G}{\min }\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+{\lambda }_1\Vert X-G\odot \hat{P}{\Vert }_F^2+{\lambda }_x{R}_x(X)+{\lambda }_g{R}_g(G)(7)$$

那么，至少存在一个 $R(G)$ 和 $G_2$，使得 $(X^{\ast },G_2)$ 是问题（6）的最小点。

$Proof$：

我们需要证明存在一个正则化项 $R(G)$ 和一个 $G_2$，使得问题（6）的最小点 $(X^{\ast },G_2)$ 也是问题（7）的最小点。

设 $R(G):=\Vert G-Z{\Vert }_F^2$，其中 $Z$ 是一个未确定的张量。此时，问题（6）变为：
$$\underset{X,G}{\min }\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+{\lambda }_1\Vert X-G\odot \hat{P}{\Vert }_F^2+{\lambda }_2\Vert G-Z{\Vert }_F^2(19)$$

将张量 $X$、$G$、$Y$、$\hat{P}$ 和 $Z$ 展平为向量：
$$y=\text{vec}(Y),x=\text{vec}(X),g=\text{vec}(G),p=\text{vec}(\hat{P}),z=\text{vec}(Z)$$
则问题（19）可以等价地表示为：
$$\underset{x,g}{\min }\Vert y-Ax{\Vert }_F^2+{\lambda }_1\Vert x-g\odot p{\Vert }_F^2+{\lambda }_2\Vert g-z{\Vert }_F^2(20)$$
其中，$A$ 表示下采样和模糊操作的矩阵。

问题（20）是一个凸优化问题，具有唯一的最小点 $(x^{\ast },g^{\ast })$。最小点满足以下方程组：（即(20)式对于$x$ 和 $g$ 分别求导）
$$\{\begin{array}{l}{A}^T(A{x}^{\ast }-y)+{\lambda }_1(x^{\ast }-g^{\ast }\odot p)=0\\ {\lambda }_{1}p\odot (p\odot g^{\ast }-x^{\ast })+{\lambda }_2(g^{\ast }-z)=0\end{array}$$

通过解上述方程组，我们可以构造 $Z$ 为：
$$Z={\text{vec}}^{-1}\left(\frac{B}{{\lambda }_1{\lambda }_{2}p}\right)$$

其中：
$$B=({\lambda }_1{p}^2+{\lambda }_2)\odot \left[(A^{T}A+{\lambda }_{1}I)\text{vec}(X^{\ast })-A^{T}y\right]-{\lambda }_1^2{p}^2\odot \text{vec}(X^{\ast })$$

此时，问题（19）的最小点为 $(X^{\ast },G_2)$，其中：
$$G_2={\text{vec}}^{-1}\left(\frac{(A^{T}A+{\lambda }_{1}I)\text{vec}(X^{\ast })-A^{T}y}{{\lambda }_{1}p}\right)$$

通过上述构造，我们证明了存在一个正则化项 $R(G)=\Vert G-Z{\Vert }_F^2$ 和一个 $G_2$，使得问题（6）的最小点 $(X^{\ast },G_2)$ 也是问题（7）的最小点。因此，定理得证。

---

根据定理1，我们的模型（6）在 $X$ 方面与模型（7）是等价的。比较这两个模型，似乎 $R_x(X)$ 可以被“吸收”到 $R_g(G)$ 中，形成一个新的 $R(G)$。由于本工作的主要关注点是 $G$，我们移除了正则化 $R(X)$，以保持模型（6）尽可能简洁明了。此外，我们希望使我们的算法易于复现。考虑到这一点，我们更倾向于在本研究中保持模型（6）的形式。

正如Xiao等人[25]所指出的，系数张量 $G$ 包含图像结构，尽管它严格意义上并不是“图像”。图2展示了一个示例。我们可以看到 $G$ 中存在清晰的纹理。这些纹理不仅包含使 $G$ 看起来像图像的空间细节，还包含光谱补偿信息以匹配 $\hat{P}$ 和 $X$。$G$ 中这些结构的存在主要归因于HRMS和PAN之间的高空间相似性。具体来说，由于HRMS和 $\hat{P}$ 的空间分辨率应该相同，HRMS和 $\hat{P}$ 之间的低频区域和高频区域的位置几乎相同。因此，它们的逐元素除法结果，即 $G$，会在相同位置产生类似的“低频”和“高频”结构。基于这一观察，$R(G)$ 的直接选择是用于图像的正则化。在本工作中，我们提出使用神经网络 $f_{\theta }(\cdot )$ 来估计系数张量 $G$。遵循深度图像先验（DIP）[26]的思想，网络结构本身可以通过隐式图像先验对其输出进行隐式正则化。然后，我们推导出所提出的方法，其优化目标如下：

$$\underset{X,\theta }{\min }\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+\lambda \Vert X-f_{\theta }(X,P)\odot \hat{P}{\Vert }_F^2(8)$$

其中，$G=f_{\theta }(X,P)$。正则化 $R(G)$ 被吸收到网络 $f_{\theta }(\cdot )$ 中，并未显式出现。对于DIP[26]的一般设置，$f_{\theta }$ 会以随机噪声作为输入。通过这种方式，网络从“无”中学习构建目标。在问题（8）中，网络 $f_{\theta }$ 被设计为以 $X$ 和 $P$ 作为输入。这有三个原因。首先，表示 $X=G\odot \hat{P}$ 显示了 $X$、$G$ 和 $P$ 之间的明显关系。因此，考虑构建从 $X$ 和 $P$ 到 $G$ 的映射是合理的。其次，除了网络结构外，这两个输入可以为 $f_{\theta }$ 提供额外的信息以构建 $G$，我们通过实验发现这对获得更好的估计非常有用。第三，除了仅优化网络参数外，我们还可以通过调整输入 $X$ 动态和逐步修改 $G$。我们将提出的问题（8）称为PSDip。

Optimization for the PSDip Model

通过深度图像先验（PSDip）进行全色锐化可以通过交替最小化[73], [74]方便地求解。设 $L(X,\theta )$ 表示问题（8）的目标函数：
$$L(X,\theta ):=\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+\lambda \Vert X-f_{\theta }(X,P)\odot \hat{P}{\Vert }_F^2(9)$$
在第 $t$ 步中，$X$ 应通过求解相应的 $X$-子问题来更新，其中网络参数 $\theta$ 在 $L(X,\theta )$ 中是固定的。该子问题没有闭式解，因此我们考虑通过应用一步梯度下降来更新 $X$。此外，我们通过实验发现，如果在 $X$-子问题中将网络输入“$X$”也视为变量（即在此步骤中，网络输入 $X$ 的梯度也用于更新 $X$），算法不会产生最佳结果。具体细节在第 IV-E 节中介绍。因此，我们通过将网络输入 $X$ 固定为上一步更新的值（即 $X_{t-1}$）来进一步简化第 $t$ 步的 $X$-子问题。这样，第 $t$ 步的简化 $X$-子问题的目标函数为：
$$L_X(X,{\theta }_{t-1}):=\Vert Y-(X\otimes K){\downarrow }_r{\Vert }_F^2+\lambda \Vert X-f_{{\theta }_{t-1}}(X_{t-1},P)\odot \hat{P}{\Vert }_F^2(10)$$
然后，$X$ 可以通过梯度下降简单地更新：
$$X_t=X_{t-1}-\alpha {\nabla }_X{L}_X(X,{\theta }_{t-1}){|}_{X_{t-1}}(11)$$
其中 $\alpha$ 是步长。网络参数 $\theta$ 通过求解相应的 $\theta$-子问题来更新，其中 $X$ 在 $L(X,\theta )$ 中是固定的。具体来说，第 $t$ 步的 $\theta$-子问题的目标函数为：
$$L_{\theta }(X_t,\theta ):=\Vert Y-(X_t\otimes K){\downarrow }_r{\Vert }_F^2+\lambda \Vert X_t-f_{\theta }(X_t,P)\odot \hat{P}{\Vert }_F^2(12)$$
我们使用 Adam [75] 来更新 $\theta$，正如大多数基于深度学习的方法所做的那样。与 $X_t$ 类似，${\theta }_t$ 通过一步更新计算：
$${\theta }_t={\theta }_{t-1}-\beta \text{AG}\left({\nabla }_{\theta }L(X_t,\theta ){|}_{{\theta }_{t-1}}\right)(13)$$
其中 $\text{AG}(\cdot )$ 表示 Adam 中的更新方向，$\beta$ 是学习率。

交替最小化的良好初始值 $\theta$ 有助于稳定更新过程，从而帮助我们的模型（6）实现更好的性能。我们看到，在 $X\approx G\odot \hat{P}$ 的两侧应用模糊算子会得到 $X\otimes K\approx G\odot (\hat{P}\otimes K)$ [25]。因此，我们通过以下方式初始化 $\theta$：
$${\theta }_0^{\ast }=\arg \underset{\theta }{\min }\Vert \hat{Y}-f_{\theta }(\hat{Y},P)\odot (\hat{P}\otimes K){\Vert }_F(14)$$
其中 $\hat{Y}\approx X\otimes K$ 表示上采样的 LRMS。由于我们无法访问 $X$，因此 $\hat{Y}$ 也被视为 $f_{\theta }$ 的第一个输入。然后，${\theta }_0^{\ast }$ 被用作问题（8）的交替最小化的初始值。在算法 1 中，我们总结了实现 PSDip 的整个过程。算法 1 的实现细节在第 III-D 节中介绍。此外，图 3 展示了 PSDip 的流程图，以便全面了解 PSDip。

---

Implementation Details

我们采用 PanNet[45] 作为网络 $f_{\theta }$ 的主干网络。该网络主要包含卷积层和跳跃连接。此外，我们在 PanNet 的最后一层添加了 ReLU激活函数，以确保输出始终为正。需要注意的是，网络结构并不是唯一指定的。在第 IV-G 节中，我们还展示了另外两种网络的结果。

为了构建 $\hat{P}$，我们首先进行 直方图匹配 生成 ${\hat{P}}^{\prime }$，即对 ${\hat{P}}^{\prime }$ 的每个波段进行平移和拉伸，使其均值和标准差与 LRMS 的对应波段匹配。为了避免分母为零，我们在 ${\hat{P}}^{\prime }$ 上添加一个小的值（$\epsilon =1\times 10^{-2}$），最终得到 $\hat{P}={\hat{P}}^{\prime }+\epsilon$。上采样的 LRMS $\hat{Y}$ 是通过对 LRMS $Y$ 进行双三次插值得到的。模糊核 $K$ 与多光谱传感器的调制传递函数（MTF）匹配[72]。

对于初始化问题（14），我们使用 Adam 优化 $\theta$，直到目标函数收敛，大约需要 8000 步。学习率固定为 $1\times 10^{-3}$。对于主问题（8）的交替最小化，我们在所有实验中设置 $\alpha =2$ 和 $\beta =1\times 10^{-3}$。在交替最小化的每一步中，$X$ 和 $\theta$ 分别通过一步梯度下降和自适应矩估计（ADAM）进行更新。权衡参数 $\lambda$ 在所有实验中设置为 0.1。当 $L(X,\theta )$ 的值变化缓慢时，迭代完成，大约需要 3000 步。也就是说，在算法 1 中，我们设置 $T=3000$。