Simultaneous Image Registration and Fusion in a Unified Framework

Posted on Edited on In
Symbols count in article: 7.4k Reading time ≈ 21 mins.
图像融合

 

Foreword

  • 2015年录用于IEEE Transactions Image Processing(TIP)
  • 第一作者chen chen, specially on low level image to image processing.

 

Conceptual Understanding

 

Abstract

  • 本文提出在同一地理位置上融合高分辨率全色图像(high-resolution panchromatic image)和低分辨率多光谱图像(low- resolution multispectral image)的新方法。
  • 化为凸优化问题，是最小化最小二乘拟合项和动态梯度稀疏性正则化的线性组合。前者是为了保持多光谱图像的准确光谱信息，后者是为了保持高分辨率全色图像的清晰边缘。why?
  • 提出在融合过程中同时配准两幅图像，这是通过动态梯度稀疏特性自然实现的。
  • 大量的实验结果表明，所提出的方法在空间和光谱质量方面明显优于其他方法。

 

Introduction

  • 多光谱(MS)图像在遥感领域应用广泛，高分辨率质谱传感器的设计受到机载存储和带宽传输基础设施限制。具有高空间分辨率的全色(Pan)灰度图像可以更方便地获得，因为它们由数量少得多的像素组成。我们期望通过图像融合获得高空间分辨率和高光谱分辨率的图像，这在文献中也被称为泛锐化。

  • 图像融合是一个典型的逆问题。第一个问题是如何从MS和Pan中保留准确的信息。

  • 许多传统方法使用投影和替换，包括主成分分析(PCA)、强度色调饱和度(IHS)、小波(wavelet)及其组合。这些方法在以下方案中执行融合:上采样、前向变换、强度匹配、组件替换和反向转换。但它们很可能受到光谱失真的影响。

  • 为解决光谱失真，提出一些变分方法，基于稍微弱的假设来制定能量函数，然后优化这样的函数以获得最佳值。这些方法也被称为基于模型的融合。然而，由于缺乏有效的模型来保存空间信息，在融合结果中可能会出现可见的伪影或模糊。此外，所有这些方法通常涉及高计算复杂性，这使得这些方法无法扩展到大规模数据集。

  • 融合中的第二个问题是如何减少失调的影响。以上方法几乎都需要在融合前进行精确的配准。然而，由于输入图像之间的显著分辨率差异，预配准相当具有挑战性。预配准后，当分辨率相差四倍时，多光谱图像上0.5像素的失准对应于全景图像上2像素的失准。

  • 本文提出了一种新的图像同时配准和融合的方法，以首字母SIRF命名。

  • 我们假设下采样后的融合图像应该接近输入的MS图像，该图像被公式化为最小二乘拟合项以保持光谱信息。

  • 受融合图像和输入Pan图像之间的地理关系的驱动，发现、定义并利用动态梯度稀疏特性来提高空间质量。重要的是，我们发现组合模型不违反遥感物理，并且动态梯度稀疏性自然地诱导精确配准，即使在严重的强度失真下。

  • 此外，方法结合了不同波段的内在相关性，这是以前很少考虑的。为了优化整个能量函数，设计了一种基于最近梯度技术的新算法。具体来说，通过分别应用快速迭代收缩阈值算法(FISTA和带回溯的梯度下降法)来有效地解决子问题。整个算法在每次迭代中保持线性计算复杂度，因此可扩展到大规模数据集。该算法可直接应用于真实数据集，无需预过滤、特征提取、训练等。

  • 最后，SIRF只有一个非敏感参数，这是与现有变分方法相比的另一个优点。大量实验结果表明，我们的方法可以显著减少光谱失真，同时保留融合图像中清晰的物体边界。

 

Modeling

notations

  • \(\mathbf{P} \in \mathbb{R}^{m \times n}\)为Pan图像

  • \(\mathbf{M} \in \mathbb{R}^{\frac{m}{c} \times \frac{n}{c} \times s}\)为MS图像。c是常数，当Pan图像的分辨率为0.6m，MS图像的分辨率为2.4m时，c = 4

  • \(\mathbf{X} \in \mathbb{R}^{m \times n \times s}\)是要融合的图像

  • \(\|\cdot\|_{F}\)是Frobenius范数

  • \(\mathbf{X}_{i, j, d}\)表示X的第i行、第j列和第d带中的元素

 

local spectral consistency

  • 许多现有方法对MS图像进行上采样，并从该上采样的MS图像中提取光谱信息。然而，向上采样的图像模糊且不准确。因此，我们只假设下采样后的融合图像接近原始的MS图像。最小二乘拟合用于模拟这种关系: \[ E_{1}=\frac{1}{2}\|\psi \mathbf{X}-\mathbf{M}\|_{F}^{2} \]
  • 其中ψ表示下采样运算符。局部光谱信息被强制与MS每个像素一致。这个函数是物理激励的，因此可以避免结果中的频谱失真。由于极低的下采样率(例如，当c = 4时为1/16)，最小化E1将是一个严重不适定的问题。没有强先验信息，X几乎不可能被准确估计

 

dynamic gradient sparsity

  • Pan图像提供了这样的先验信息。

  • 由于遥感图像通常是分段平滑的，因此它们的梯度往往是稀疏的，非零值对应于边缘。

  • 此外，当图像对齐良好时，这些边缘的位置应与平移图像上的位置相同。根据参考图像，证明了稀疏性不是固定的而是动态的。我们将具有这种性质的数据称为动态梯度稀疏信号/图像。

  • 全变分可以促进梯度的稀疏性。然而，在该正则项中没有整合来自Pan图像的参考信息。与以前的工作不同，我们的方法鼓励动态梯度稀疏性。除了先前方法试图使用的先验信息之外，我们还注意到跨不同波段的内部相关性，因为它们是相同陆地对象的表示。因此，不同波段的梯度应该是组稀疏的(group sparse)。L1范数鼓励稀疏，L2,1范数鼓励群体稀疏。因此，我们提出了一个新的能量函数来同时促进动态梯度稀疏性和组稀疏性: \[ \begin{aligned} E_{2} &=\|\nabla \mathbf{X}-\nabla D(\mathbf{P})\|_{2,1} \\ &=\sum_{i} \sum_{j} \sqrt{\sum_{d} \sum_{q}\left(\nabla_{q} \mathbf{X}_{i, j, d}-\nabla_{q} \mathbf{P}_{i, j}\right)^{2}} \end{aligned} \]

  • 其中q = 1,2，D(P)表示复制P到s波段。有趣的是，当没有参考图像时，即P = 0，结果与矢量全变分(VTV)的结果相同。

  • 情形(a)-(d)分别具有组稀疏数8、4、4、2。

  • 结合这两个能量函数，图像融合问题可以表述为:

\[ \begin{array}{l} \min _{\mathbf{X}}\left\{E(\mathbf{X})=E_{1}+\lambda E_{2}\right. \\ \left.\quad=\frac{1}{2}\|\psi \mathbf{X}-\mathbf{M}\|_{F}^{2}+\lambda\|\nabla \mathbf{X}-\nabla D(\mathbf{P})\|_{2,1}\right\} \end{array} \]

  • 所提出的动态梯度稀疏性仅迫使支持集相同，而梯度的符号以及信号的幅度不需要相同。这些属性使其在对比度反转下不变，并且对光照条件不敏感。它可以应用于不同来源或不同采集时间的图像融合。方法还可以同时联合融合多个波段，这提供了对噪声的鲁棒性。

 

simultaneous image registration

  • 提出在融合过程中同时配准图像。一方面，多光谱图像被锐化到更高的分辨率。这使我们能够更准确地记录图像。另一方面，未对准被逐渐消除，并且图像可以被更精确地融合。我们反复运行这两个过程，直到收敛。
  • 基于配准中使用的特征，现有的图像配准方法可以分为基于特征的配准和基于像素(或基于强度)的配准。基于特征的方法依赖于从图像中提取的标志。更感兴趣的是基于强度的配准，它可以在统一的优化方案中与融合相结合。
  • 图像配准最重要的组成部分之一是能量函数的测量相似度。
  • 使用动态梯度稀疏性来保存空间信息。任何错位都会增加梯度的稀疏性。动态梯度稀疏性可以自然地用作相似性度量。我们可以修改同时进行图像配准和融合的能量函数: \[ E(\mathbf{X}, \mathcal{T})=\frac{1}{2}\|\psi \mathbf{X}-\mathbf{M}\|_{F}^{2}+\lambda\|\nabla \mathbf{X}-\nabla \mathcal{T}(D(\mathbf{P}))\|_{2,1} \tag{a} \]   其中\(\mathcal{T}\)是需要估计的变换
  • 将提出的相似性度量与现有的方法进行了比较。图中的蓝色曲线显示了不同测量的响应，所有这些测量都可以在零平移时找到最佳对准。在增加强度失真和重新缩放之后，图(b)中所示的源图像的外观与原始Lena图像的外观不一致。红色曲线表示的结果表明，只有RC和所提出的方法可以处理这种强度失真。

 

Algorithm

  • 目标是最小化能量函数(a)。先通过固定T来解决关于X的问题，然后通过固定X来解决关于T的问题，对于X子问题: \[ E(\mathbf{X})=\frac{1}{2}\|\psi \mathbf{X}-\mathbf{M}\|_{F}^{2}+\lambda\|\nabla \mathbf{X}-\nabla \mathcal{T}(D(\mathbf{P}))\|_{2,1} \tag{b} \]

  • (b)一个明显的凸函数。第一项是光滑的，而第二项是不光滑的。在FISTA框架中解决这个子问题。已经证明，对于一阶方法，FISTA可以达到最优收敛速度。

  • 关于T的第二子问题可以写成: \[ \min E(\mathcal{T})=\|\nabla \mathbf{X}-\nabla \mathcal{T}(D(\mathbf{P}))\|_{2,1} \tag{c} \]

  • 下图总结了提出的同时进行图像配准和融合(SIRF)的泛锐化(pan-sharpening)算法：.

  • 这里ψ表示ψ的逆算子。L是\(\psi^{T}(\psi \mathbf{X}-\mathbf{M})\)的李普希茨常数。

  • 观察到解决方案是基于Xk 和 Xk-1更新的，而在以前的方法中使用的Bregman方法仅基于Xk更新X，也就是为什么我们的方法收敛得更快的原因。

  • SIRF 子问题1中，L=1， \[ \mathbf{X}^{k}=\arg \min _{\mathbf{X}}\left\{\frac{1}{2}\|\mathbf{X}-\mathbf{Y}\|_{F}^{2}+\lambda\|\nabla \mathbf{X}-\nabla \mathcal{T}(D(\mathbf{P}))\|_{2,1}\right\} \tag{d} \]

  • 取\(\mathbf{Z}=\mathbf{X}-\mathcal{T}(D(\mathbf{P}))\) ，带入式(d)： \[ \mathbf{Z}^{k}=\arg \min _{\mathbf{Z}}\left\{\frac{1}{2}\|\mathbf{Z}-(\mathbf{Y}-\mathcal{T}(D(\mathbf{P})))\|_{F}^{2}+\lambda\|\nabla \mathbf{Z}\|_{2,1}\right\} \]

  • 这个替代问题是一个VTV去噪问题，Xk由\(\mathbf{Z}^{k}+\mathcal{T}(D(\mathbf{P}))\)更新。VTV去噪算法的慢版本是基于FISTA框架来加速求解(d)的，这在Algorithm 2中进行了总结

  • 线性算子定义: \(\mathcal{L}(\mathbf{R}, \mathbf{S})_{i, j, d}=\) \(\mathbf{R}_{i, j, d}-\mathbf{R}_{i-1, j, d}+\mathbf{S}_{i, j, d}-\mathbf{S}_{i, j-1, d}\) ，对应的逆算子： \(\mathcal{L}^{T}(\mathbf{X})=(\mathbf{R}, \mathbf{S})\) with \(\mathbf{R}_{i, j, d}=\mathbf{X}_{i, j, d}-\mathbf{X}_{i+1, j, d}\) and \(\mathbf{S}_{i, j, d}=\mathbf{X}_{i, j, d}-\mathbf{X}_{i, j+1, d} \cdot \mathbb{P}\) 是投射算子确保： \(\sum_{d=1}^{s}\left(\mathbf{R}_{i, j, d}^{2}+\right.\) \(\left.\mathbf{S}_{i, j, d}^{2}\right) \leq 1,\left|\mathbf{R}_{i, n, d}\right| \leq 1\), and $|_{m, j, d}| $

  • 用带回溯的梯度下降法(gradient descent method with backtracking))来解决第二个子问题(c)，好处是在最小化时没有引入新参数变量
  • 设\(\mathbf{r}=\nabla \mathbf{X}-\nabla \mathcal{T}(D(\mathbf{P}))\) 是残差，对(c)的紧密估计： \[ E(\mathcal{T}) \approx \sum_{i} \sum_{j} \sqrt{\sum_{d}\left(\nabla_{1} \mathbf{r}_{i, j, d}\right)^{2}+\left(\nabla_{2} \mathbf{r}_{i, j, d}\right)^{2}+\epsilon} \]
  • 根据链式法则，梯度： \[ \nabla E(\mathcal{T})=-\frac{\partial E(\mathcal{T})}{\partial \mathbf{r}} \nabla \mathcal{T}(D(\mathbf{P})) \frac{\partial \mathcal{T}}{\partial \theta} \]
  • \(\epsilon\)是一个小量,\(\nabla \mathcal{T}(D(\mathbf{P}))\)表示图像亮度梯度，θ表示变换T的参数，在这里，变换主要是仿射或平移
  • \(\partial \mathcal{T} / \partial \theta\)是基于一阶泰勒近似计算的。我们设置初始步长t0= 1，η = 0.8。所有的运算都是线性的。函数值是在两幅图像的重叠区域上计算的。为了避免平凡解，如放大黑暗区域，在这里使用归一化函数值(除以重叠的像素M)。当没有重叠时，函数值将是无穷大。

 

Experiments

simulation

  • 该方法在Quickbird、Geoeye、SPOT和IKONOS卫星的多光谱数据集上进行了验证。Pan图像的分辨率范围为0.41米至1.5米,所对应的MS图像都具有较低的分辨率，c = 4，并且包含蓝、绿、红和近红外波段。为方便起见，仅显示RGB波段。由于缺乏同一场景的多分辨率图像，原始图像被视为 Ground truth，低分辨率图像从 Ground truth下采样得到。

 

visual comparison

  • 结果
  • 为了更好地可视化，图6和图7以相同的比例显示了与地面真实情况相比的误差图像。从这些误差图像中，可以清楚地观察到光谱失真、块状伪影和模糊。

 

quantitative analysis

  • 为了评估不同方法的融合质量，使用四个度量来衡量光谱质量，一个度量来衡量空间质量。光谱度量包括合成中的相对无量纲全局误差(ERGAS) 、光谱角度映射器(SAM) 、通用图像质量指数(Q-average) 和相对平均光谱误差(RASE) 。过滤后的相关系数(FCC) 被用作空间质量度量。此外，峰值信噪比(PSNR)、均方根误差(RMSE)和平均结构相似度(MSSIM) 用于评估与地面真实值相比的融合精度。

 

efficiency comparison

 

translation

  • 在图1所示的平移图像中添加了3个像素的人工水平平移。SIRF的估计翻译如图10 (a)所示。通过大约100次迭代，可以非常准确地恢复真实的translation。

 

real-world datasets

  • 图像是由IKONOS多光谱成像卫星获取，其中包含预配准的以其捕获分辨率拍摄的Pan和MS图像。

 

Conclusion and Discussion

  • 基于动态梯度稀疏性的特性，提出了一种新的变分模型，用于在统一的框架下同时进行图像配准和融合。该模型自然地结合了来自高分辨率全色图像的梯度先验信息和来自低分辨率多光谱图像的光谱信息。
  • 比现有的方法有几个优点:首先，所提出的动态梯度稀疏性可以直接利用来自Pan图像的尖锐边缘;第二，通过合并不同波段的内部相关性来联合锐化图像，而大多数现有的方法是基于波段间的融合；最后，尽管由于空间分辨率的不同，Pan和MS图像之间的配准相当困难，但我们的方法可以在融合过程中同时配准两个输入图像，直接作用于输入图像，而无需任何预滤波和特征提取。
  • 方法被证实在空间和光谱质量方面始终优于现有技术。

 

Code

background knowledge