最小二乘影像匹配的算法设计与实现

最小二乘影像匹配的算法设计与实现
【摘要】：本文对最小二乘影像匹配的原理进行了深入的分析，描述了算法实现的具体步骤，并结合实现过程中出现的问题提出了值得注意的关键点，最后通过实验效果显示算法实现的效果及匹配精度。

关键字：最小二乘匹配；灰度重采样；几何纠正
1 引言
在数字摄影测量时代，随着对数据处理的自动化程度要求的提高，影像匹配的精度成为影响测图精度提高的主要因素之一。

德国的Ackermann教授提出的一种高精度的影像匹配方法——最小二乘影像匹配，充分利用了影像窗口内的信息进行平差计算，使影像匹配可以达到1/10甚至1/100像素的高精度，即影像匹配精度可达子像素等级。

它不仅可以被用于生产一般的数字地面模型，生产正射影像图，而且可以用于控制点的加密及工业上的高精度量测。

另外，最小二乘影像匹配中可以非常灵活地引入各种已知参数的特点使它具有整体平差的功能。

正是由于最小二乘匹配方法具有灵活、可靠和高精度的特点，因此，受到了广泛的重视，得到很快的发展。

本文将结合实验效果对最小二乘匹配算法的原理以及实现过程进行具体的阐述和深入分析。

2 最小二乘匹配的原理
2.1 概念引入
影像匹配中判断影像相似的测度很多，其中有一种是“灰度差平方和最小”。

若将灰度差记为余差v，则上述判断可写为：
vv=min (1)
因此，它与最小二乘法是一致的。

但是在一般情况下，它没有考虑影响灰度中存在的系统误差，而仅认为影像灰度只存在偶然误差，即：
v=g1(x,y)-g2(x,y) (2)
这就是一般的按“灰度差平方和最小”原则进行影像匹配的数字模型。

若在此系统中引入系统变形参数，仍然按“灰度差平方和最小”的原则解求变形参数，就构成了最小二乘影像匹配系统。

2.2 原理分析
影像灰度的系统变形有两大类：一类是辐射畸变；另一类是几何畸变。

由此产生了影像灰度分布之间的差异。

因此，在最小二乘影像匹配匹配的过程中，根据左右影像目标区所有像素的灰度分布差异，通过迭代不断修正引入的变形参数值，使匹配精度逐步提高。

由于最小二乘影像匹配是非线性系统，迭代时必须已知初匹配结果。

左右目标区的确定也即以初匹配结果中的左右同名点为中心像素的矩形区域，一般窗口的尺寸较小，只考虑一次畸变：
x’=a0+a1▪x+a2▪y (3)
y’=b0+b1▪x+b2▪y
此处的坐标是以局部窗口内某点为原点的局部坐标，考虑编程实现，一般以目标区窗口中具有最小行列号的像素为原点。

再考虑影像的辐射畸变，以右方影像相对左方影像的线性灰度畸变予以描述，即得关系式：
g1(x,y)+n1(x,y)=h0+h1▪g2(a0+a1▪x+a2▪y,b0+b1▪x+b2▪y)+n2(x,y) (4)
(4)式经线性化后，可得最小二乘影像匹配的误差方程式：
v=c1▪dh0+c2▪dh1+c3▪da0+c4▪da1+c5▪da2+c6▪db0+c7▪db1+c8▪db2-△g （5）
其中c1,c2,c3,c4,c5,c6,c7,c8是（4）线性化后得到的系数，dh0,dh1,da0,da1,da2,db0,db1,db2是待定参数的改正值，它们的初值分别为：h0=0,h1=1,a0=0,a1=1,a2=0,b0=0,b1=0,b2=1；观测值△g是左右目标区中相应像素的灰度差。

按（5）式在目标区（一般只对右方影像作畸变处理）逐个像元建立误差方程式，经迭代平差至收敛，得到系统误差的改正系数，最后再对左方影像目标区内的像素位置按梯度进行加权，得到左影像目标区的特征点，依据解算的几何纠正参数即可求得右影像目标区中匹配的同名点。

3 最小二乘影像匹配的算法实现
3.1 影像初匹配
采用一种简单快速的匹配算法，获取左右影像同名点的概略位置，匹配精度达到整像素级或几个像素级即可。

初匹配的精度过低将导致最小二乘影像匹配的迭代不收敛，匹配精度高的算子对应的匹配过程较复杂将影响匹配效率。

一般可以用相关系数法作为初匹配算法，而要求更高的情况下可以选用SIFT算子。

本实验程序采用的是矢量夹角相关的相关系数法作为初匹配算法。

先打开左影像并用Harris算子获取特征点，然后打开右影像，利用相关系数法匹配出右影像上相关系数达到0.90以上的同名点。

在匹配过程中由于是全图搜索最佳匹配点，建立金字塔影像可以明显提高匹配效率，这种情况对于图幅较大的影像尤为明显。

3.2 几何变形改正
影像获取的过程中，左右影像不可避免地都存在几何变形，而考虑到影像匹配的条件，一般只处理左右影像中的一幅即可，此处仅对右影像进行几何纠正。

利用初匹配得到的同名点，在右影像上以同名点为中心像素确定一个适当大小的矩形目标区，然后根据最小二乘影像匹配解求的几何变形改正参数a0,a1,a2,b0,b1,b2，将目标区影像窗口的局部影像坐标进行改正变换，具体变换公式可参照（3）式。

3.3 灰度重采样与辐射畸变的改正
原目标区中像素的位置经几何变形改正后，其局部影像坐标一般不为整数，即当欲知不位于矩阵( 采样) 点上的原始函数g( x，y) 的数值时就需要进行内插，此时称为重采样，意即在原采样基础上再一次采样。

通常的重采样方法有：最近邻法、双线性内插法以及双三次卷积法。

综合考虑采样的质量和效率，本实验程序采用双线性内插法：考虑投影点周围4个相邻像元灰度值，并根据投影点与周围4个整像素点的距离作为权重计算输出像元的灰度值。

灰度重采样之后，根据最小二乘影像匹配解求的辐射畸变改正参数h0,h1对上述灰度重采样结果做辐射改正：
g2(x,y)’=h0+h1▪g2(x,y) (6)
此处采样的辐射畸变改正模型是线性的，其依据是考虑算法整体的复杂程度即运行效率，同时，实验证明线性的改正模型能够满足要求的精度。

3.4 最小二乘平差解算变形参数
利用左影像目标区窗口与经过几何、辐射校正的右影像目标区窗口内的灰度阵列g1与
h1▪g2(x,y)+h0，依据式（5）对窗口内各个像素建立平差方程，通过每次迭代得到8个参数的改正值，并将改正后的参数继续代入下一次迭代计算的方程，如此循环，直到收敛，即相关系数小于前一次迭代后所求得的相关系数，则可认为迭代结束。

但是考虑到初匹配中可能存在粗差同名点，当误差较大的同名点作为最小二乘影像匹配的初值时，迭代将不收敛，某次相关系数小于前一次迭代的相关系数也是很有可能的，如果仅以此条件为迭代结束条件，得到的结果将有存在较大粗差的危险。

解决的办法是增加迭代结束条件：几何纠正参数的常数项a0和b0的最低位数必须达到0.1及以下的精度，亦可以对相关系数的值进行要求，实验证明0.970比较合适。

如此设定迭代结束条件后，能够跳出循环的结果几乎都是符合要求的正确匹配点，但是对于不收敛的那些点将一直循环迭代下去，显然这是没有意义的，因此可以在迭代循环的内部对循环次数设定一个阈值，以防止死循环的产生，对于超过迭代次数阈值的初匹配点予以删除，并继续下一对同名点的精匹配。

3.5 计算最佳匹配点位
影像匹配的目的是获取同名点，通常是以待定目标点建立一个目标影像窗口，窗口中心点即为目标点。

但是，在高精度影像相关中，必须考虑目标窗口的中心点是否为最佳匹配点。

根据最小二乘匹配的精度理论可知：匹配精度取决于影像的灰度梯度大小|gx|和|gy|，因此可用梯度的绝对值为权，在左方影像目标区窗口内对坐标作加权平均：
|gx(7)
xi=∑x▪|gx|/∑|
|gy
yi=∑y▪|gy|/∑|
以(xi,yi)作为目标点坐标，它的同名点坐标可由最小二乘影像匹配所求得的几何换参数求得：
x2=a0+a1▪x+a2▪y (8)
y2=b0+b1▪x+b2▪y
其中，(x,y)是左方影像目标区窗口内求得的特征点的局部坐标，按式（8）求得的点(x2,y2)即右影像目标区中求得的匹配点的局部坐标。

4 实验效果
在最小二乘影像匹配中，左右影像的目标区窗口大小的设置关系到程序运行的效率以及匹配的精度。

从实验的效果来看，目标区域选择29×29的区域效果会比较好，选择区域较大或较小时效果都不是很好。

下面在此条件下，结合实验效果对程序中其他参数的设置进行分析说明：
4.1 迭代次数限定值
迭代次数的设定，其目的是防止不收敛情况的发生而使程序进入死循环。

限定的迭代次数太少，则会错误地将部分初匹配精度不太高但足以使精匹配过程收敛的正确点处理为不收敛的粗差点，这样虽然提高了匹配点的精度，但是匹配的点数过少，影响实用性；限定的迭
代次数过多，则容易造成粗差点许多不必要的运算，影响匹配效率。

实验显示迭代次数设为15次较好，如下所示：
图1 迭代次数5，匹配点数3
图2 迭代次数10，匹配点数6
图3 迭代次数15，匹配点数6
4.2 迭代结束条件
最小二乘影像匹配的主要原理很少注重说明代结束的条件，但该部分也是影响匹配效果的重要环节。

为了控制匹配精度，一般考虑每次迭代中几何纠正参数的两个常数项a0和b0的改正数的精度，一般设为0.01至0.05较好，这个值也反映了匹配精度的子像素等级。

值得注意的是，每次迭代后计算的相关系数在迭代条件中的运用，若仅用相关系数小于前一次迭代后所求得的相关系数是非常危险的，即初匹配结果在最小二乘影像匹配中迭代不收敛时，相邻两次迭代后计算的相关系数的变化是跳动的。

因此，不能仅使用这个条件为迭代终止条件。

较为合理的设置是考虑每次迭代后计算得到的相关系数，设置为高于初匹配的匹配相关系数阈值0.9即可，一般为0.93—0.95，不能设置得过高，因为最小二乘影像匹配的原理仍然是基于灰度的，而由于各种误差因子的存在，左右影像使对应的区域的相关系数并不一定是1，该值设置得过高反而得不到已经匹配正确的点。

另外，经过实验发现：用每次迭代的单位权中方差代替相关系数具有更好的效果，因为相关系数考虑了最方影像目标区的灰度分布情况，而单位权方差仅仅只与左右目标窗口内的灰度差异有关。

现以如下实验结果予以说
明：
图4 单位权方差值设为200，匹配点数为2
图4：匹配点位即相关系数：
Left Graph Right Graph Match the correlation coefficient (234.876 391.961) (236.21 420.545) 0.988442
(67.3299 442.936) (68.7959 471.619) 0.986844
(163.065 489.552) (164.193 518.708) 0.991607
图5 单位权方差值设为300，匹配点数为5
图5：匹配点位即相关系数：
Left Graph Right Graph Match the correlation coefficient (159.025 400.411) (158.795 425.97) 0.968341
(234.876 391.961) (236.373 420.552) 0.963173
(67.3299 442.936) (68.7959 471.619) 0.986844
(113.934 487.516) (115.034 516.477) 0.987542
(163.065 489.552) (164.193 518.708) 0.991607
图6 单位权方差值设为400，匹配点数为6
图5：匹配点位即相关系数：
Left Graph Right Graph Match the correlation coefficient (107.729 435.019) (102.813 454.516) 0.972375
(159.025 400.411) (158.795 425.97) 0.968341
(234.876 391.961) (236.373 420.552) 0.963173
(67.3299 442.936) (68.7959 471.619) 0.986844
(113.934 487.516) (115.034 516.477) 0.987542
(163.065 489.552) (164.193 518.708) 0.991607
如上所示：目视情况下，图4和图5匹配的点位是正确的，但是图5匹配的点数更多，而且相同点位的坐标和相关系数是一致的，即更低的相关系数阈值对应更多的匹配点位。

从图6可以看出：蓝色圆圈标记的点位显然是错误匹配点，也就是单位权方差阈值400的条件过于宽松，但是其匹配相关系数并不是最低的，由此可知是单位权方差抑制了该匹配粗差点的出现。

因此，使用单位权方差代替相关系数作为迭代终止条件显得更为合理。

5 结束语
本文针对影像匹配问题，对最小二乘法影像匹配的算法进行了深入研究，阐述该算法的匹配原理，并用程序对该算法的特点进行了测试。

最小二乘影像匹配算法匹配精度高，能够很方便地引入变形参数进行整体平差，从而抵消左、右影像之间的几何差异和辐射差异，既考虑了辐射误差对影像匹配的影响，同时也考虑了几何畸变对影像匹配的影响。

避免了只考虑灰度的匹配方法出现的误匹配，如目标点与最终计算得出的匹配点虽然灰度相似性很大，但它们可能并不是同名点，因为它们的几何位置可能存在很大的差异，而最小二乘影像匹配的方法由于考虑了目标区与匹配区的几何畸变，所以能有效地避免此类问题的产生。

另外，对编程实现过程中迭代结束条件以及粗差点的自动剔除的研究，也在一定程度上保证了匹配的准确度、精度和效率。

【参考文献】：
[1］张剑清，潘励，王树根.摄影测量学［M］．武汉: 武汉大学出版社，2003.6．
[2] 李福金，秦志伟.最小二乘影像匹配算法设计与实现.辽宁锦州：测绘与空间地理信息第34卷第2期.。