高光谱拼接算法（六）RANSAC 误匹配剔除

本篇算法仓位：RANSAC

1. 特征匹配中的几何约束

上一篇我们梳理了 SIFT 的多种变体与改进方向，现在再回到拼接管线本身：
之前我们利用 SIFT 完成了特征匹配，拿到了初步的匹配点对，但即使经过比率测试，结果中仍然存在误匹配。原因是：

比率测试只关心描述子空间的"数值接近性"，不关心两张图像之间的几何一致性。

比如一对双胞胎站在一起，我们凑近为其中一人拍了一张单人照，然后又拉远镜头拍了一张合照，那么这两幅图像间的特征就很容易匹配到另外一个人上。

换句话说，比率测试只能回答：“这两个局部纹理看起来很相似。”，但无法回答：“这两个点在真实三维空间中是否对应同一个位置？”

这就引出了拼接管线中的下一步：几何验证，其核心思路是：

利用匹配点对之间的几何约束关系，识别并剔除那些"看起来像但实际位置不对应"的误匹配。

在图像拼接场景中，若场景近似平面，这个几何约束通常是单应性矩阵。
而求解单应性矩阵并同时剔除误匹配的最经典方法，是 RANSAC，来自 1981 年的论文 Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography，就和 SIFT 一样，它也拥有丰富的变体和现代化改进。

不过，与 SIFT 不同，RANSAC 的职责十分明确：它并不负责提取或匹配特征，而只是在已有匹配点的基础上完成几何验证。
因此，它的发展并没有形成像 SIFT 那样庞杂的谱系，而是不断结合新技术实现一个目标：如何更高效准确地从包含误匹配的数据中估计正确的几何模型。

具体内容之后会再详细展开。

2. 单应性矩阵

RANSAC 的前置内容是单应性矩阵（Homography），虽然是矩阵形式，但其实使用函数来概括说明可以更好地理解它的作用：

描述同一场景的特征点在两幅图像中对应关系的函数（变换矩阵），输入特征点坐标，输出在另一角度图像中的对应坐标。

大致明确其作用后，再来展开原理：

2.1 什么是单应性？

理论证明：两个视角拍摄同一平面，其对应关系可以由一个单应性矩阵描述。

因此，单应性描述的是同一平面场景在两幅图像之间的投影变换关系。

展开来，假设场景中存在一个平面，比如农田、墙面、地面等，相机从两个不同角度拍摄，那么两幅图像之间的对应点满足：

\[\begin{bmatrix} x' \\ y' \\ 1 \end{bmatrix} \propto H \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} \]

其中 \(H\) 就是 \(3\times3\) 的单应性矩阵：

\[H= \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \]

这里有一个细节：单应性矩阵是齐次坐标下的变换（基础线代），这意味着 \(H\) 乘以任意非零常数后，表示的是同一个变换。因此 \(H\) 的自由度实际上只有 8，最后一个元素通常固定为 1，即 \(h_{33}=1\)。

全是参数看起来十分复杂，但其实这类变换并不突兀，摆几个例子：

变换类型	矩阵形式	自由度	效果
平移	\(\begin{bmatrix}1&0&t_x\\0&1&t_y\\0&0&1\end{bmatrix}\)	2	横向/纵向移动
旋转	\(\begin{bmatrix}\cos\theta&-\sin\theta&0\\\sin\theta&\cos\theta&0\\0&0&1\end{bmatrix}\)	1	绕原点旋转
缩放	\(\begin{bmatrix}s_x&0&0\\0&s_y&0\\0&0&1\end{bmatrix}\)	2	放大缩小
仿射	\(\begin{bmatrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\0&0&1\end{bmatrix}\)	6	平行四边形变换
透视	全 8 参数	8	任意四边形变换

可以看到，平移、旋转、缩放、仿射变换都可以表示为单应性矩阵的特殊形式，这也是为什么单应性矩阵能够描述相机旋转、平移、视角变化等复杂场景下的图像关系。

而结合到这个系列来说：在无人机航拍拼接中，由于地面近似平面，且无人机通常在同一高度飞行，单应性矩阵可以很好地描述相邻航带图像之间的几何关系。

现在，把矩阵运算展开成公式，从 \((x,y)\) 到 \((x',y')\) 的实际映射关系是这样的：

\[x' = \frac{h_{11}x + h_{12}y + h_{13}}{h_{31}x + h_{32}y + 1} \]

\[y' = \frac{h_{21}x + h_{22}y + h_{23}}{h_{31}x + h_{32}y + 1} \]

显然，问题在于求解参数。

2.2 如何求解单应性矩阵？

这部分其实就是矩阵形式的解方程：一对匹配点 \((x,y)\) 和 \((x',y')\) 提供了两个方程。将上面的映射关系改写为线性形式：

\[x' (h_{31}x + h_{32}y + 1) = h_{11}x + h_{12}y + h_{13} \]

\[y' (h_{31}x + h_{32}y + 1) = h_{21}x + h_{22}y + h_{23} \]

再写成矩阵形式，可以发现每对点贡献 \(2\times9\) 的行：

\[\begin{bmatrix} x & y & 1 & 0 & 0 & 0 & -x'x & -x'y & -x'\\ 0 & 0 & 0 & x & y & 1 & -y'x & -y'y & -y' \end{bmatrix} \begin{bmatrix} h_{11} \\ h_{12} \\ h_{13} \\ h_{21} \\ h_{22} \\ h_{23} \\ h_{31} \\ h_{32} \\ h_{33} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} \]

这代表每对匹配点可以提供 2 个约束，\(H\) 有 8 个自由度，理论上只需要 4 对不共线的匹配点即可求解整个单应性矩阵。

这种求解方法叫直接线性变换（Direct Linear Transform，DLT）：

把原本非线性的几何关系，通过齐次坐标改写成关于未知参数的线性方程组，然后利用线性代数求解。

前置内容就到这里，现在问题来了：

我们经过比率测试后有几十甚至上百对匹配点，其中还混有误匹配。用哪些点来解单应性矩阵？

这正是 RANSAC 要解决的问题。

3. RANSAC 算法

RANSAC 的核心思想很朴素，流程也较为简单，但却十分有效：

与其试图一次性用所有数据拟合模型（容易被少量错误数据严重带偏），不如反复随机采样最小子集来拟合模型，然后看每个模型能获得多少其他数据的支持，最终选择支持度最高的那个模型。

3.1 算法流程

RANSAC 的完整步骤如下：

1. 随机采样： 从当前的匹配点集中随机选取 4 对匹配点，进行退化检查（重复点，共线，近似共线等）直至取点合理，使用最小集是为了满足求解条件的同时尽可能防止被误匹配点污染模型。
2. 拟合模型： 用第一步得到的 4 对匹配点通过 DLT 求解单应性矩阵 \(H\)。
3. 内点判定： 用 \(H\) 将所有匹配点中图像 A 的特征点映射到图像 B 的坐标系中，计算映射点与实际匹配点之间的重投影误差，判定内外点并记录。

这里的重投影误差其实就是映射坐标和真实坐标间的距离，之所以称为"重投影"，是因为单应性矩阵根据图像 A 的点重新预测了其在图像 B 中的位置，相当于再次投影到图像 B 上，因此称为重投影误差：

\[e_i = \sqrt{(x'_i - \hat x'_i)^2 + (y'_i - \hat y'_i)^2} \]

其中 \((x'_i, y'_i)\) 是第 \(i\) 对匹配点在图像 B 中的原始坐标，\((\hat x'_i, \hat y'_i)\) 是 \(H\) 映射后的坐标。

而内外点的定义是这样的：我们根据对精度的要求设置距离阈值 \(T\) ，通常取 \(3\sim5\) 像素。
如果 \(e_i < T\) ，则判定该匹配为内点（inlier），否则为外点（outlier）。

最后记录这次采样的内点数量，进入下一步：

4. 重复迭代： 重复步骤 1~3 共 \(N\) 次。\(N\) 次迭代结束后，选取内点数量最多的那次结果对应的 \(H\) 作为最优模型。
5. 最终优化： 对最优模型得到的所有内点，在 \(H\) 基础上使用最小二乘法迭代优化总的重投影误差，重新估计单应性矩阵，得到更精确的 \(H\)，利用更多样本降低随机误差。

RANSAC 的过程其实类似抽签：经过足够多次的随机抽样，抽到一组"全内点"的最小集，从而定位到正确的模型。

3.2 迭代次数与置信度

一个涉及实现的问题是：

需要迭代多少次才能保证找到较准确的模型？

而这和数据本身相关，不能去盲目寻找“官方配置”，一个重要因素是数据中外点的比例，也就是 SIFT 本身得到的匹配质量。
我们可以进行如下推导，假设外点比例为 \(\epsilon\)，则一次采样中 4 个点全部为内点的概率为：

\[P_{good} = (1-\epsilon)^4 \]

如果我们希望以概率 \(p\) 至少找到一组全是内点的采样，则需要进行 \(N\) 次迭代：

\[1 - (1 - (1-\epsilon)^4)^N = p \]

整理得到：

\[N = \frac{\log(1-p)}{\log\left(1-(1-\epsilon)^4\right)} \]

代几个值看看：

外点比例 \(\epsilon\)	单次采样全内点概率	\(N\)（\(p=0.99\)）	\(N\)（\(p=0.999\)）
10%	65.6%	4	6
30%	24.0%	17	25
50%	6.25%	72	107
60%	2.56%	179	268
70%	0.81%	567	850
80%	0.16%	2876	4313
90%	0.01%	45996	68995

可以看到两个结论：

1. 外点比例在 50% 以内时，RANSAC 是高效的，几十到一百多次迭代就能以 99% 的置信度找到较准确内点组合。
2. 一旦外点比例超过 70%，所需的迭代次数会急剧上升。到 90% 时，已经需要近 5 万次迭代。

这说明提高初始匹配质量十分重要。比例测试能够有效降低外点比例，从而显著减少 RANSAC 所需迭代次数，因此二者具有很强的互补性。

在真实实验中，迭代次数少有超过 2000 的设置。
而且，随着 RANSAC 的不断改进，自适应迭代次数的逻辑已经占据主流，这便是下一篇的内容了。

代码的结果图如下，但是因为是我的随手一拍，外点率在 80%，所以结果仍然可以看到不准确的地方，不过在和比例测试的对比中，就可以看到明显的剔除效果：