Automatic Stitching for Hyperspectral Images Using Robust Feature Matching and Elastic Warp

作者：Yujie Zhang; Zhiying Wan; Xingyu Jiang; Xiaoguang Mei
年份：2020
机构：武汉大学
期刊/会议：IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing
原文地址：Automatic Stitching for Hyperspectral Images Using Robust Feature Matching and Elastic Warp

高光谱图像不仅包含空间信息，还包含丰富的光谱信息，已广泛应用于农业、城市规划等领域，但单幅图像难以覆盖大面积。因此，需要对各个部位进行拍照，并应用图像拼接技术来获得全景高光谱图像。当场景的视点变化很大时，传统方法会出现重影问题。为了获得高精度合成全景图，本文提出了一种使用鲁棒特征匹配和弹性扭曲的高光谱图像自动图像拼接算法。我们的方法包含两个阶段。第一阶段是选择一个波段作为参考波段，获取单波段全景图。特别是，我们通过尺度不变特征变换提取特征点。然后，我们提出了一种称为多尺度顶级 K 秩保留算法的有效算法，用于在两组点之间建立稳健的点对应关系。接下来，我们采用稳健的弹性扭曲来获得每个波段的全景图。第二阶段是基于第一阶段得到的变换拼接所有剩余的波段，并将所有波段的信息融合在一起，得到最终的全景高光谱图像。大量实验证明了我们提出的方法的有效性。第二阶段是基于第一阶段得到的变换拼接所有剩余的波段，并将所有波段的信息融合在一起，得到最终的全景高光谱图像。

Ⅰ 介绍

随着无人机（UAV）的发展，高光谱图像在遥感中的应用得到了改进。然而，为了实现高空间分辨率，瞬时视场会很小。单幅图像很难获得大面积的信息。因此，研究高光谱图像的图像拼接是有意义的。

图像拼接包括图像匹配和图像对齐两个部分：

图像匹配包含两个步骤：

基于点、边、角或者其他实体关系建立两幅图像之间的对应关系，常见的算法包括SIFT、SURF，对于高光谱图像拼接，光谱空间SIFT (SS_SIFT) 可以用来提取3D特征，但处理时间不是实时的，对于实时应用不可行。
去除错误匹配，常用的方法包括随机采样共识 (RANSAC) ，除此之外还有最大似然估计样本共识和渐进样本共识，但并不适用于所有情况。除此之外，还有一些非参数的方法，如向量场一致性 (VFC) 、局部保持匹配 (LPM) 但存在较大异常值时，仍不准确。

图像对齐中：

全局变换：为每个输入图像估计全局变换矩阵，使用该矩阵对齐图像，这种方法是不容许视差的。
尽可能投影APAP：局部自适应投影变换，可以在重叠区域取得好的效果，但是仍然存在失真的问题；
SPHP：利用投影变换，逐渐转变为跨图像的全局相似变换。

RGB图像的拼接方法有了很大的改进，但还不能直接应用于高光谱图像。与 RGB 图像不同，高光谱图像的每个像素都有数十或数百个与其波段对应的幅度值。因此，拼接算法应该是高效的，以避免大量的计算。对于基于无人机的高光谱图像，图像可能会发生尺度变化、视差、非自然旋转或非刚性变换。

作者提出的方法：

首先选择一个单波段作为参考波段，获得单波段全景图。利用SIFT提取特征，通过多尺度Top K透视 (mTopKRP) 去除误差匹配。接着利用稳健弹性扭曲重新投影获取无缝全景。
使用第一个参考波段获得矩阵拼接所有剩余波段并融合所有波段信息，获得最终高光谱全景图。

Ⅱ 方法

方法流程图如下：

高光谱图像拼接流程图

单波段图像拼接

使用SIFT算法提取图像 $I_1$ 和 $I_2$ 的特征点，获得 $N$ 对匹配对 $S=\{(x_i,y_i)\}_{i=1}^N$ ，其中 $x_i$ 和 $y_i$ 为像素坐标。为了获得稳健的匹配对，需要去除错误的匹配对。

最佳内点集 $R^*$ 最优解为
$\begin{equation*} R^*=\arg \min C\left(R;S,\lambda \right).\tag{1} \end{equation*}$
根据mTopKRP，特征点的KNN相似度量衡量特征点是否正确匹配，代价函数 $C$ 表示为：
$\begin{equation*} C\left(\mathbf {p}; S,\lambda \right)=\sum _{i=1}^Np_{i}\left(c_{i}-\lambda \right) +\lambda N \tag{2} \end{equation*}$
其中 $c_i=\frac{1}{M}\sum_{m=1}^M{D_{km}(\sigma(x_i),\sigma(y_i))}$ ，向量 ${\bf p}$ 与 $S$ 相关， ${\bf p}\in(0,1)$ ， $p_i=1$ 代表匹配点对是正确的，为内点。反之 $p_i=0$ 表示匹配点对是错误的，为外点。 ${D_{km}(\sigma(x_i),\sigma(y_i))}$ 表示特征点匹配对前 $k$ 的距离。

代价函数用到了多尺度策略，当找到前 $K$ 个最相似的近邻之后，在不同尺度 $K$ 下建立损失函数。定义不同 $K$ 为 $K =\{K_m\}_{m=1}^M$ 。由于代价函数优化时可能会出现零值，需要添加约束项，确保最大化内点数量的同时保持损失值最小。

给定匹配点对，多有特征点的K近邻被计算出来，接着计算损失值 $\{c_i\}_{i =1}^N$ ，正确匹配会导致出现负值使代价函数变小。
$\begin{equation*} {p_i} = \left\lbrace {\begin{array}{c}1,{c_i} \leq \lambda \\ {0,{c_i} > \lambda } \end{array}}. \right. \tag{3} \end{equation*}$
此时，最优内点可以被下式确定：
$\begin{equation*} {R ^*} = \left\lbrace {i|{p_i} = 1,i = 1,2,\ldots,N\rbrace .} \right. \tag{4} \end{equation*}$
接着进行图像图像扭曲变换。直接使用全局变换会导致错误的匹配，出现视差或者鬼影。作者在这里采取了鲁棒弹性变换避免投影变形。假设 ${\bf p}=[x,y,1]^T$ 和 ${\bf q}=[u,v,1]^T$ 为匹配的特征点
$\begin{equation*} \left[ \begin{array}{c}u\\ v\\ 1 \end{array} \right] \sim \left[ \begin{array}{ccc}h_{11}&h_{12}&h_{13}\\ h_{21}&h_{22}&h_{23}\\ h_{31}&h_{32}&h_{33} \end{array} \right] \left[ \begin{array}{c}x\\ y\\ 1 \end{array} \right] \tag{5} \end{equation*}$
(5)式可以表示为
$\begin{equation*} {\bf 0}_{3\times 1} = \left[ \begin{array}{ccc}\bf 0_{1\times 3} & -{\bf p}^{T}& v {\bf p}^{T} \\ {\bf p}^{T}&{\bf 0}_{1\times 3} & -u{\bf p}^{T} \\ -v {\bf p}^{T} &u{\bf p}^{T} &{\bf 0}_{1\times 3} \end{array} \right] \cdot {\bf h=Ah},\Vert {\bf h}\Vert ^{2} =1 \tag{6} \end{equation*}$
其中 ${\bf h}=[h_{11},h_{12},\cdots,h_{33}]^T$ ， ${\bf h}\in{\mathbb R}^{9\times 1}$ ， ${\bf A} =[a_{1},a_{2},\cdots,a_{N}]^T$ ， ${\bf a} \in{\mathbb R}^{2\times 9}$ ， $i=1,2,\cdots,N$ 。

利用直接线性变换求解(6)式
$\begin{equation*} \mathbf {h}=\arg \min _{\mathbf {h}}\sum ^{N}_{i=1}\Vert \mathbf {a}_{i}\mathbf {h}\Vert ^{2}=\arg \min _{\mathbf {h}}\Vert \mathbf {Ah}\Vert ^{2}. \tag{7} \end{equation*}$
将图像 $I_1$ 分为 $C_1\times C_1$ 块。每个区域的中心点为 ${\bf p}^*$ ，任意位置 ${\bf p}_i$ 的局部单应性可以利用MDLT计算
$\begin{equation*} \mathbf {h}_{*}=\arg \min _{\mathbf {h}}\sum ^{N}_{i=1}\Vert \omega _{*}^{i}a_{i}\mathbf {h}\Vert ^{2}=\arg \min _{\mathbf {h}}\Vert W_{*}\mathbf {Ah}\Vert ^{2} \tag{8} \end{equation*}$
其中 $\omega_*^i=\max(\exp(-\|{\bf p}_*-{\bf p}_i\|^2/\sigma^2),\lambda)$ ， $\sigma$ 为尺度因子且 $\lambda\in[0,1]$ ， $W_*={\rm diag}([w_*^1,w_*^1,w_*^2,w_*^2,\cdots,w_*^N,w_*^N])$ .

然而非重叠区域也有可能出现扭曲变形，为了克服这个问题，作者的解决方法是加权优化算法。

图像 $I_2$ 映射到图像 $I_1$ 上需要单应矩阵变换， $x$ 方向的变换方程可以表示为：
$\begin{equation*} {g_\mu }(x,y) = \mu g(H(x,y)),(0 < \mu (x) < 1). \tag{9} \end{equation*}$
$\mu$ 的值从图像左侧到右侧由1到0逐渐减小。 $y$ 方向的变换方程可以表示为：
$\begin{equation*} {h_\mu }(x,y) = \mu h(H(x,y)),(0 < \mu (x) < 1). \tag{10} \end{equation*}$
将扭曲和全局相似性 ${ \bf S}$ 结合
$\begin{equation*} {J_A} = {\sum \limits _{i = 1}^n {\left\Vert {{\mathbf {S}}{{\mathbf {p}}_i} - {{\mathbf {q}}_i}} \right\Vert } ^2} \tag{11} \end{equation*}$
其中 ${ \bf S}=\begin{bmatrix}c & -s & t_x \\ s & c & t_y\end{bmatrix}$ ， $\{{ \bf p}_i,{ \bf q}_i\}$ 为匹配对。

将单应性和相似性矩阵混合：
$\begin{equation*} {{ \mathbf {H}}_q} = {\mu _h}{\mathbf {H}} + {\mu _s}{{\mathbf {H}}_s} \tag{12} \end{equation*}$
$\mu_h$ 在原图 $I_1$ 上由0到1变化， $I_2$ 的变换如下所示：
$\begin{equation*} {{\mathbf {H}}_p} = {{\mathbf {H}}_q}{{\mathbf {H}}^{ - 1}}. \tag{13} \end{equation*}$

算法：使用鲁棒特征匹配和弹性扭曲自动拼接高光谱图像。

输入图像： $I_i$ ， $i=1,2,\cdots,K$ ，每张图由 $K$ 个波段组成。

输出：无视差的高光谱全景图。

选择第 $\lfloor M/2\rfloor$ 个波段作为参考波段
利用SIFT在参考波段提取特征点匹配对 $S=\{(x_i,y_i)\}_{i=1}^{N}$
基于多尺度K近邻获得每个特征点的排序表
利用式(6)计算代价值 $\{c_i\}_{i=1}^{N}$ 获得内点对 $R^*$
计算单应矩阵 ${\bf H}$ ，利用(7)式计算相似矩阵 ${ \bf H }_s$
定义非重叠区域并利用式(9)和(10)计算全局投影变换
利用(12)计算 ${\bf H}$ 和 ${ \bf H }_s$ 的全局扭曲
对于所有点 $i=1$ 到 $K-1$ 重复5-8
利用(14)和(15)将所有图映射到一张图上
对剩下 $M-1$ 个波段使用相同的全局扭曲变换
所有波段图像线性融合
拼接后的图像可以表示为

$\begin{equation*} I=\frac{\sum _{i=1}^Kw_{i}I_{i}}{\sum _{i=1}^Kw_{i}} \tag{14} \end{equation*}$

其中 $w_i$ 为 $I_i$ 的权重
$\begin{equation*} w_{i}(x,y)=1-\frac{d_{i}(x,y)}{d_{i,\max }} \tag{15} \end{equation*}$
其中 $d_i(x,y)$ 为点 $(x,y)$ 到图像 $I_i$ 中心的距离， $d_{i,\max}$ 表示 $I_i$ 中像素值非0点到中心的最大距离。

Ⅲ 实验分析

拼接场景由四川双力光谱成像技术有限公司提供的GaiaSky-mini传感器采集，是一款无人机载高光谱传感器，波长范围为400-1000nm。我们通过选择第 92 (700.2 nm)、第 47 (547.6 nm) 和第 13 (436.5 nm) 波段来合成伪彩色图像。

SIFT与其他特征点检测算法对比

SIFT、SURF和SS-SIFT的特征点匹配，以及由RANSAC去除的错误匹配：使用SIFT和557(44.42%)使用RANSAC进行内部匹配；80个使用SURF和49(61.25%)使用 RANSAC 进行内部匹配；和82个使用 SS-SIFT 的暂定匹配和64(78.04%)使用 RANSAC（从上到下和从左到右）的内层匹配。黑线连接假定的匹配点。蓝色的连接正确的匹配点。

利用F-Score评估性能
$\begin{equation*} \text{F-score}= \frac{2 \times \text{Precision} \times \text{Recall}}{(\text{Precision}+\text{Recall})} \tag {16} \end{equation*}$

RANSAC、VFC、LPM、和我们的 mTopKRP 在合成孔径雷达 (SAR)、彩色红外航空照片 (CIAP) 和 UAV 上的定量比较。Recall、Precision 和 F-score（从左到右）关于累积分布。图例中报告了平均精度 (AP)、平均召回率 (AR) 和平均 F 分数 (AF)。

通过使用整个特征集在 46 个遥感图像对上构建局部相邻结构，相对于累积分布的 F-score。最佳平均 F-score 及其阈值λ 是 AF = 0.9905 和 λ = 0.7.

缝合性能的定性比较。第一行是ANAP、NISwGSP、ELA和我们提出的单波段方法的对应拼接结果。其余行是每个结果图像的三个代表性区域。(a) ANAP 的结果。(b) NISwGSP 的结果。(c) ELA 的结果。(d) 我们扭曲的结果。

最终的高光谱拼接结果：RGB彩色图像与176波段高光谱图像在像素级融合后，通过我们的方法得到彩色显示的高光谱遥感图像。我们通过选择第 92 (700.2 nm)、第 47 (547.6 nm) 和第 13 (436.5 nm) 波段来合成伪彩色图像。

EOG为能量梯度，表示灰度等级的改变；DFT反映图像空间的全局活动；Variance表示图像灰度值相对于均值的离散程度。它们的定义如下：
$\begin{align*} &\ \text{EOG}=\sum _{x}\sum _{y}\lbrace [f(x\!+\!1,y)\!-\!f(x,y)]^{2}\!+\![f(x,y+1)-f(x,y)]^{2}\rbrace \tag{17} \\ &\ \text{DFT}=\frac{1}{MN}\sum _{\mu =0}\sum _{\nu =0}\sqrt{\mu ^{2}+\nu ^{2}}P(\mu,\nu) \tag{18} \\ &\ \text{Variance}=\frac{1}{M\times N}\sum _{x}\sum _{y}(f(x,y)-\lambda)^{2} \tag{19} \end{align*}$

6种算法不同指标的比较

Ⅳ 频谱分析

评价指标：光谱角度距离 (SAD) :
$\begin{equation*} \ \text{SAD} =\arccos \left[\frac{\sum {\mathbf {a}_{i}\mathbf {b}{_{i}}}}{\sqrt{ \sum {\mathbf {a}_{i}^{2}}}\sqrt{\sum {\mathbf {b}_{i}^{2}}}}\right] \tag{20} \end {equation*}$
其中 $i$ 为波段序号。

image

光谱分析

Ⅴ 结论

本文的主要贡献是结合了SIFT算法、mTopKRP算法、鲁棒弹性扭曲和光谱匹配算法，设计了一种用于高光谱图像的自动图像拼接算法。首先，我们提出了一种称为 mTopKRP 的稳健特征匹配的新的错配去除方法，该方法基于特征对应的稳定相邻拓扑关系。接下来，采用鲁棒弹性扭曲来减少图像的失真。在混合一个波段的 18 张图像后，我们对其他 175 个波段使用相同的全局扭曲和融合参数重复此过程。mTopKRP算法保证了特征点提取和匹配的准确性，鲁棒的弹性扭曲避免了一些形状畸变，光谱匹配算法使光谱与原始光谱一致。最后，在真实高光谱图像上的实验结果表明，我们的方法可以获得高精度的拼接结果。

使用鲁棒特征匹配与弹性扭曲的高光谱图像自动拼接