基于Harris-Affine的宽基线立体影像LSM匹配方法

姚国标^{1, 2}，邓喀中¹，张力²，杨化超¹，艾海滨²，杜全叶²

(1. 中国矿业大学 环境与测绘学院，

国土环境与灾害监测国家测绘地理信息局重点实验室，江苏 徐州，221116；

2. 中国测绘科学研究院，北京，100830)

摘要：针对宽基线影像的高畸变特点，提出基于几何一致性约束的Harris-Affine高精度最小二乘匹配(Least Square Matching, LSM)方法。首先，用Harris-Affine算子在影像高斯多尺度空间中提取仿射不变特征，并根据特征区域灰度矩阵与主梯度方位来实现特征的几何归一化，继而提取特征描述符，接着采用一种由粗到精的特征匹配策略来渐进地获取正确率占优的特征匹配。然后，通过LSM迭代方法来补偿同名像点的定位误差，LSM迭代所需的良好的几何畸变初值由同名特征区域灰度矩阵与主梯度方位求取，而良好的辐射畸变初值通过同名区域像元灰度的最小二乘线性拟合法求取。实验结果表明：该算法鲁棒有效，且实现了亚像素匹配精度。

关键词：宽基线立体影像；Harris-Affine算子；特征匹配；最小二乘匹配；亚像素

中图分类号：TP391          文献标志码：A         文章编号：1672-7207(2014)08-2661-08

Least square matching method for wide baseline stereo images based on Harris-Affine features

YAO Guobiao^{1, 2}, DENG Kazhong¹, ZHANG Li², YANG Huachao¹, AI Haibin², DU Quanye²

(1. Key Laboratory for Land Environment & Disaster Monitoring of SBSM,

School of Environment Science and Spatial Informatics, China University of Mining & Technology, Xuzhou 221116, China;

2. Academy of Surveying and Mapping, Beijing 100830, China)

Abstract: Due to the complex geometric and radiometric distortions on widely separated images, the traditional stereo algorithms were not suited for them, thus a novel high-accuracy least square matching (LSM) method based on Harris-Affine features and geometric consensus constraint was proposed. Firstly, Harris-Affine features were detected in multiple scale space of Gaussian images. Local features were normalized separately based on intensity moment and the principle gradient orientation of feature region, and then these normalized results were described respectively. Subsequently, a novel feature matching strategy from coarse to fine was implemented on candidates in order to improve the correctness of correspondences. Secondly, LSM algorithm was performed on these conjugate features so as to obtain the fine matches. As the LSM must require good initial value for iteration convergence, the geometrical distortion initial value was calculated based on intensity moment and the principle gradient orientation of correspondent affine regions, and the radiometric distortions initial value was obtained through least square fitting method which was implemented on normalized correspondent regions. Experiments on wide baseline images of real-world scenes indicate that the proposed algorithm has significant superiority of matching precision.

Key words: wide baseline stereo images; Harris-Affine operator; feature matching; least square matching; sub-pixel

宽基线影像匹配在机器人视觉^[1]、三维重建^[2]、信息融合^[3]等领域发挥着关键作用。常用的匹配方法一般分为基于灰度与基于特征2类。基于灰度的匹配算法，如归一化互相关(normalized cross correlation，NCC)算法^[4]、最小二乘影像匹配(least square matching, LSM)算法^[5]等，具有高精度特点，非常适合于影像的优质匹配，但它们需要已知待匹配区域的几何与辐射(亮度)畸变量，因此加大了其在宽基线影像匹配应用中的难度。基于特征的匹配算法对影像间的各种畸变不敏感，但仅可达到粗匹配精度。总之，上述2类算法优势与挑战并存，如何有效融合二者的精度与鲁棒优势，是目前宽基线影像匹配研究的热点与难点之一^[6-9]。源于机器人视觉领域的特征提取及描述方法的出现，为解决宽基线影像的自动匹配难题提供了思路。特征匹配方法包括3步^[10]：(1) 特征提取；(2) 特征描述符提取；(3) 相似性判断。现有的角点特征提取算子如Frstner^[11]、Harris^[12]等，现有的区域特征提取算子如尺度不变特征变换^[13](Scale Invariant Feature Transform, SIFT)、最大稳定极值区域^[14](Maximally Stable External Regions, MSER)以及Harris-Affine^[15] (Ha-A)等，区域特征较角点特征具有更稳定的抗变形特点，其中Ha-A则又兼具良好仿射不变性与多量性优势。流行的特征描述符，如以特征区域梯度信息为基础的SIFT描述符^[13]、直方图描述符^[16]等，以特征边缘信息为基础的形状描述符^[17]等，文献[18]通过多类影像数据的比较实验表明，SIFT描述符效果最为满意。特征相似性判断主要是根据描述符的距离测度或比率测度进行同名特征判断。为提高特征匹配精度，Xiao和Shah^[19]、杨化超等^[5]分别提出了基于Harris和SIFT的高精度最小二乘匹配算法。对存在复杂畸变的宽基线影像，应用LSM方法的要点是如何获得几何与辐射畸变参数初值。杨化超等^[5]的方法有效利用了SIFT特征尺度与方位信息作为最小二乘匹配几何畸变初值进行匹配迭代优化；而Xiao和Shah^[19]的方法则是在仿射参数的离散化空间内采用逐步搜索的方法来获得。为实现宽基线影像的精确、鲁棒自动匹配，上述方法存在如下缺陷：(1) 最小二乘匹配迭代的几何参数初值获取过程较为繁琐；(2) 在LSM匹配迭代过程中，对应区域的几何畸变与辐射畸变参数具有一定的相关性，即当任一参数不准确时，势必会影响到其余参数甚至算法的收敛可靠性，而上述文献方法均未能提供可靠的辐射畸变初值。与Harris和SIFT相比，仿射不变区域特征Ha-A对存在复杂畸变的宽基线影像具有较好的适应性，且基于Ha-A固有的区域信息易于为LSM提供良好初始值。在此，本文作者提出基于几何一致性约束的Ha-A高精度最小二乘匹配算法。首先，该算法采用一种由粗到精的特征匹配策略获取高正确率的特征匹配。然后，通过LSM算法对特征匹配结果进行精确迭代优化，LSM迭代所需的良好的几何变换初始值由同名特征区域灰度矩阵与主梯度方位求取，而良好的辐射畸变初值根据同名区域像元灰度的最小二乘线性拟合法求取。最后，采用多组宽基线影像实验来验证本文算法的可行性与精度优势。

1  算法

1.1  特征匹配

基于文献[15]算法，在参考影像与待匹配影像中分别提取众多Ha-A特征，特征区域呈椭圆形，而特征描述需在规则区域中(如圆形或矩形)进行，因此，首先需要对特征实施归一化处理，之后再进行特征描述。图1所示为某对Ha-A特征和归一化过程，其中x和x′为特征的质心；U和U′为2阶灰度矩阵，包含了特征区域仿射信息；α和α′表示特征主梯度方位。

特征区域归一化与描述原理如下：(1) 仿射归一化。保持特征区域的形状和倾角不变，然后扩大k(k取3.0)倍，得到特征测量区和，见图1第1列白色椭圆包围的区域。和按(1)式进行仿射变换，得到和，见图1第2列所示结果。(2) 旋转归一化。为克服同名特征间的方向差异，用区间直方图分别统计和区域内各像元梯度方向，并取直方图的峰值作为主梯度方位，则和主梯度方位统计结果为α和α′，再进一步按式(2)构成旋转矩阵R和R′，和按式(3)进行旋转变换，最终得到方向对齐的和，见图2第3列所示结果。(3) 特征描述。SIFT描述符即在归一化区域和中提取。

，     (1)

，  (2)

，       (3)

宽基线立体影像视差变化幅度较大，导致对应区域图像内容差异明显。因而，利用模式识别自动寻找同名特征很容易出现误匹配，为此，采用一种由粗到精的匹配策略来渐进提高匹配正确率：(1) 初始匹配。利用SIFT描述符与最/次近邻距离比率(NNDR)获得初始匹配，并利用K-d树搜索提高匹配效率，NNDR阈值取值0.75。(2) 基于随机采样一致性(Random Sample Consensus, RANSAC)算法逐步排除误匹配。从步骤(1)匹配的结果中随机和递归地抽取8对不同的特征点，利用式(4)与最小二乘法分别计算所有情况的基本矩阵F与单应矩阵H_g。根据式(5)对所得的基本矩阵F与单应矩阵H_g联合进行最小二乘误差(Least Squares Error, LSE)估计。LSE所对应的8对同名特征点就是8对最佳的匹配点，LSE对应的和即是最优的基本矩阵与单应矩阵，然后采用F和H_g的双重几何约束排除误匹配，其中极线约束阈值和单应映射约束阈值分别设置为3.50像素和3.00像素。

        (4)

式中：x和x′为同名特征点坐标；F和H_g分别为特征点应满足的极线约束关系和单应映射关系。

    (5)

式中：E为最小二乘误差估计；m为初始匹配特征数。

图1  特征区域归一化过程

Fig. 1  Normalization of feature regions

1.2  基于特征匹配结果的最小二乘匹配

1.2.1  LSM匹配模型

设和为Ha-A任意特征匹配，以x为中心取大小为(2w+1)×(2w+1)邻域。利用单应变换矩阵H近似对应邻域窗口之间几何变换，则有下式成立：

              (6)

令，，，则。

式中：A，t，v分别为邻域窗口之间的仿射、平移及透视变换参数。令p = [p, q]^T∈[-w, w]，p′=[p′，q′]^T，为p按式(6)进行单应变换后的坐标。则以x和x′为中心的2个相关窗口W和W′可表达为

             (7)

基于式(7)，进一步引入辐射变换参数h₀和h₁，可建立基于单应变换模型的最小二乘匹配条件方程，线性化后可得最小二乘匹配的误差方程式^[5]。其中，需要迭代估计的参数为。

LSM需要较为良好的迭代参数初值。因平移和透视畸变相对较小，其相应的迭代参数初值可分别设为：。

1.2.2  仿射变换初值求取

利用图1说明同名特征仿射畸变初始矩阵A的求取方法。设与为同名特征，分别确定区域和，二者满足如式(8)表达的仿射变换关系，它们分别经过各自的归一化，得到区域和，如图1所示，在此令=，进一步由式(1)与式(3)联立得到式(9)，然后比较式(8)与式(9)，得式(10)。至此，根据式(10)便可计算各对同名特征邻域间的仿射变换矩阵A。

             (8)

         (9)

            (10)

通过多组实验验证，利用A可较好地恢复同名邻域的几何一致性，如图2第1列所示，说明通过同名特征区域灰度矩阵与主梯度方位求取A具有准确性和有效性。

1.2.3  辐射变换初值求取

利用图2来说明辐射畸变初值求取方法。从图2第1列白色四边形内部分别提取像元信息，得灰度分布和，见图2第2列，其中直接提取，则通过双线性插值得到，且此处x=x′∈[-w, w]²，考虑到相对于的线性辐射畸变，则可得

           (11)

图2  同名邻域窗口辐射变换初值求取

Fig. 2  Calculation initialization of radiometric distortions between correspondent regions

又考虑到和的可能存在的位置偏差，故对中各像元灰度的计算采用如下策略：以当前像元为中心取八邻域进行加权平均得到的灰度作为当前像元灰度。由式(11)进一步得到误差方程式

     (12)

其中，方程式系数为：c₁=g(x)，c₂=1；常数项l=，像元数m=9，权值，c为标量且满足，σ为高斯核大小，设定σ=20。按式(12)建立误差方程式，并写成矩阵形式为

               (13)

式中：C表示方程式系数矩阵；B=[h₀, h₁]^T，由最小二乘法解得

            (14)

按式(14)便可计算相对于的辐射畸变初值B，利用B对进行辐射纠正，结果记为g″(x′)，由图2可以看出：g″(x′)与的辐射信息具有较好一致性，表明通过同名区域像元灰度的最小二乘拟合法求解B是有效的。

1.2.4  LSM匹配优化

对任意特征匹配与，便可利用上文方法为LSM匹配提供几何与辐射变换初值，迭代求取最优的单应变换矩阵H以及辐射变换参数h₀和h₁。迭代终止条件可设为W，W′间相关系数r达到阈值条件(本文设置为0.95)。迭代完成后，按式(15)对x′的坐标进行精确补偿。下面随机选取2对同名点，点号分别记为Ⅰ和Ⅱ，来直观说明LSM匹配迭代过程。表1给出了LSM单应变换窗口、匹配误差ε(单位：像素)、相关窗口以及相关系数ρ随迭代次数的变化，其中，匹配误差ε利用已知的单应变换矩阵H₀按式(16)计算得到。为展现几何与辐射变换初值对LSM匹配迭代的影响，图3所示为基于3种初值情况的LSM匹配误差随迭代次数的变化。由图3可知：(1) 仅提供几何变换初值情况：即按仿射变换初值求取方法提供初值，而辐射畸变初值分别设为h₀=0，h₁=1；(2) 仅提供辐射变换初值情况：即按辐射变换初值求取方法提供初值，而几何畸变初始矩阵设为单位矩阵A=E；(3) 提供几何与辐射变换初值情况：即按本文方法同时提供几何与辐射变换初值。

            (15)

            (16)

由表1可以看出：即使2个邻域窗口之间存在显著的几何与辐射畸变，LSM匹配亦能快速收敛，且随着迭代次数递增，相关性迅速提高，同时匹配点误差得到较好的优化补偿，表明：本文方法能够为LSM匹配提供良好的几何与辐射畸变初值；LSM匹配能够显著地提高匹配点的定位精度。图3表明：仅提供几何变换初值时，LSM匹配亦可收敛；而仅提供辐射变换初值，LSM匹配却难以收敛；同时提供几何与辐射变换初值，LSM匹配则能够更为快速、准确地收敛。

表1  基于几何&辐射变换初值的LSM匹配迭代

Table 1  LSM iteration of correspondent neighborhoods and windows based on initialization of geometric & radiometric

图3  基于3种初值情况LSM匹配误差随迭代次数的变化

Fig. 3  LSM deviation changing with iterations based on three kinds of initializations

2  算法测试与分析

2.1  实验数据与精度指标

基于Visual Studio 2010的C++高效编程语言及VLFeat计算机视觉库编写本文的算法程序，实验平台为主频3 GHz，内存2 GB的PC机。选取如图4所示3组数据进行算法测试(数据来源：http://lear.inrialpes.fr/people/mikolajczyk，各组影像真单应变换H₀已知)。(1) Boat影像数据(图4第1行)，分辨率为850×680，影像间存在约90°旋转变换和2倍的尺度变换；(2) Clown影像数据(图4第2行)，分辨率为800×640，影像间存在约60°视角变换；(3) 无人机(Unmanned Aerial Vehicle, UAV)影像数据(图4第3行)，分辨率为850×680，影像间存在显著的视角、旋转与尺度变换。另外，各组影像均存在不同程度的辐射畸变。通过式(17)计算匹配均方根误差 (单位：像素)对算法性能进行客观测试。

        (17)

式中：n为最终匹配点数；和为同名像点。

2.2  实验结果与分析

选取5种方法用于对比试验，分别为ASIFT算法^[20]、MSER算法、Ha-A算法、本文算法Ⅰ和本文算法。ASIFT算法首先按模拟的姿态角对参考影像与待匹配影像分别进行重采样，得到众多影像，然后利用SIFT算法对各像对进行匹配，最后把匹配点坐标归化到原始影像中，ASIFT算法具有全仿射不变性优点，但缺陷为时间开销过大；MSER算法与Ha-A算法均是基于本文提出的由粗到精的特征匹配策略来获取同名像点；本文算法Ⅰ与本文算法的区别在于LSM提供初值方法不同：前者仅提供几何变换初值，而后者同时提供几何和辐射变换初值。

为检验LSM策略对匹配精度的优化效果，从各组数据分别随机抽取19对特征匹配及相应的LSM匹配结果，它们的精度对比如图5所示。基于本文算法的各组数据于各阶段的匹配结果见图4，其中图4(c)中直线为根据最终匹配点估算的同名核线；表2所示为基于5种匹配算法分别针对本文3组影像数据的实验结果，其中t为匹配耗时，n为最终的匹配点数，利用本文算法的匹配点估算单应变换关系并根据式(6)进行影像叠加，它们的对应结果如图6所示。

对存在显著几何与辐射畸变的Boat影像、Clown影像以及UAV影像，本文算法均能获得数量可观、精度为亚像素的匹配点，验证了本文算法对宽基线影像的有效性与精确性。如表2中本文算法针对3组影像的匹配点数分别为294对、152对和304对，相应的匹配精度分别达到0.63像素、0.81像素和0.53像素，各组匹配结果用图4(c)中的十字点表示。从图6可见：本文算法的影像叠加结果相当清晰，图5(c)中青色核线对应性良好，均说明本文算法的匹配准确性。

图4  基于本文算法各组数据的匹配结果

Fig. 4  Matching results with different groups of data based on proposed algorithm

图5  各组数据LSM匹配前后精度对比结果

Fig. 5  Comparison results of accuracy before and after performing LSM based on different data

表2  不同匹配算法实验结果综合对比

Table 2  Comprehensive comparisons using different matching algorithms

本文算法相比于其余算法具有精度优势。表2表明：ASIFT算法、MSER算法和Ha-A算法的匹配精度分别仅能达到像素级，原因在于特征检测与特征匹配的相对独立，使得特征匹配点往往难以精确对应。本文算法即是在特征匹配结果的基础上再进一步通过LSM匹配进行迭代优化，表1和图5均表明特征匹配点精度能够得到较大程度提高。

LSM匹配对几何畸变和辐射畸变的初值要求严格，而本文方法能够较好地满足要求。表2表明：本文方法的匹配效率比本文方法I的高，且匹配精度亦较高。原因是在LSM匹配阶段，本文方法较本文方法I更具备可靠的几何与辐射畸变初值。然而，图5表明，对于部分特征匹配，利用LSM亦不能有效补偿它们的匹配误差，针对这类粗差点，进一步利用RANSAC算法进行剔除，最终得到高精度匹配点如图4(c)所示的十字点。

图6  利用本文算法对各组影像数据叠加结果

Fig. 6  Registration results based on three groups of images by using proposed method

宽基线影像的畸变程度会对匹配结果造成一定影响。从表2可知：针对Clown影像数据，各算法的匹配精度结果均低于其余2组数据。这是因为Clown立体影像存在的仿射(透视)畸变比Boat影像和UAV影像的畸变程度要更大、更复杂。

3  结论

(1) 本文算法具有可行性，能够有效补偿特征点的定位误差，且匹配精度较高。

(2) 由于基础特征所限，本文算法在纹理均匀区域难以得到匹配点，下一步将通过融合与Ha-A特征具有空间分布互补性的特征(如SIFT、MSER等)来提高本文算法的鲁棒性。

参考文献：

[1] Dutta A, Kar A, Chatterji B N. A new approach to corner matching from image sequence using fuzzy similarity index[J]. Pattern Recognition Letters, 2011, 32(5): 712-720.

[2] Gruen A. Development and status of image matching in photogrammetry[J]. Photogrammetric Record, 2012, 26(137): 36-57.

[3] Li H, Ma Y, Liang K. Rapid matching algorithm for hyperspectral image based norm sifting[J]. Chinese Optics Letters, 2012, 10(1): 72-80.

[4] 罗三定, 陈海波. 基于区域增长的自适应窗口立体匹配算法[J]. 中南大学学报(自然科学版), 2005, 36(6): 1042-1047.

LUO Sanding, CHEN Haibo. Stereo matching algorithm of adaptive window based on region growing[J]. Journal of Central South University (Science and Technology), 2005, 36(6): 1042-1047.

[5] 杨化超, 张书毕, 张秋昭. 基于SIFT的宽基线立体影像最小二乘匹配方法[J]. 测绘学报, 2010, 39(2): 187-194.

YANG Huachao, ZHANG Shubi, ZHANG Qiuzhao. Least squares matching methods for wide base-line stereo images based on SIFT features[J]. Acta Geodaetica et Cartographica Sinica, 2010, 39(2): 187-194.

[6] Bastanlar Y, Temizel A, Yardimci Y. Improved SIFT matching for image pairs with scale difference[J]. Electronics Letters, 2010, 46(5): 346-354.

[7] 姚国标, 邓喀中, 杨化超, 等. 基于几何约束的倾斜立体影像匹配[J]. 中国矿业大学学报, 2013, 42(4): 676-682.

YAO Guobiao, DENG Kazhong, YANG Huachao, et al. Matching methods for oblique stereo images on the basis of geometry constraint[J]. Journal of China University of Mining & Technology, 2013, 42(4): 676-682.

[8] 智喜洋, 张伟, 侯晴宇, 等. 影响测绘相机匹配精度的辐射指标量化[J]. 光学精密工程, 2012, 20(2): 387-394.

ZHI Xiyang, ZHANG Wei, HOU Qingyu, et al. Quantification for radiometric specification impacting on matching accuracy of mapping camera[J]. Opt Precision Eng, 2012, 20(2): 387-394.

[9] Chiou Y C, Liang Y T. An effective corner detection method using subpixel edge detector and gaussian filter[J]. Sensor Review, 2010, 30(1): 93-101.

[10] 门朝光, 边继龙, 李香. 基于迭代传播的小基高比立体匹配方法[J]. 中南大学学报(自然科学版), 2012, 43(4): 1362-1369.

MEN Chaoguang, BIAN Jilong, LI Xiang. Iterative diffusion based stereo matching method for small baseline[J]. Journal of Central South University (Science and Technology), 2012, 43(4): 1362-1369.

[11] Frstner W, Glch E. A fast operator for detection and precise location of distinct points, corners and centers of circular features[C]//ISPRS Conference on Fast Processing of Photogrammetric Data. Interlaken, Switzerland, 1987: 281-305.

[12] Harris C, Stephens M J. A combined corner and edge detector[C]//Proceedings of the 4th Alvey Vision Conference. Manchester, Britain, 1998: 147-151.

[13] Sadek R, Constantinopoulos C, Meinhardt E, et al. On affine invariant descriptors related to SIFT[J]. Slam Journal on Imaging Sciences, 2012, 5(2): 652-687.

[14] Matas J, Chum O, Urban M, et al. Robust wide-baseline stereo from maximally stable extremal regions[J]. Image Vision Computing, 2004, 22(10): 761-767.

[15] Mikolajczyk K, Schmid C. Scale & affine invariant interest point detector[J]. International Journal of Computer Vision, 2004, 60(1): 63-86.

[16] 王梅, 屠大维, 周许超. SIFT特征匹配和差分相乘融合的运动目标检测[J]. 光学精密工程, 2011, 19(4): 892-899.

WANG Mei, TU Dawei, ZHOU Xuchao. Moving object detection by combining SIFT and differential multiplication[J]. Opt Precision Eng, 2011, 19(4): 892-899.

[17] 易盟, 郭宝龙, 严春满. 结合优化梯度滤波与投影不变的航拍视频配准[J]. 光学精密工程, 2012, 20(3): 651-660.

YI Meng, GUO Baolong, YAN Chunman. Aerial video registration combining optimal gradient filters and projective invariant[J]. Opt Precision Eng, 2012, 20(3): 651-660.

[18] Mikolajczyk K, Tuytelaars T, Schmid C. A comparison of affine region detectors[J]. International Journal of Computer Vision, 2005, 60(1): 163-186.

[19] Xiao J J, Shah M. Two-frame wide baseline matching[C]//9th IEEE International Conference on Computer Vision. Nice, France, 2003: 603-609.

[20] Morel J M, Yu G S. ASIFT: A new framework for fully affine invariant image comparison[J]. Slam Journal on Imaging Sciences, 2009, 2(2): 438-469.

(编辑  赵俊)

收稿日期：2013-10-15；修回日期：2013-11-19

基金项目：国家自然科学基金资助项目(41272389，41371438)；江苏省基础研究计划(自然科学基金)青年基金资助项目(BK20130174)；江苏省普通高校研究生科研创新计划项目(CXZZ13_0937)；江苏高校优势学科建设工程项目(SZBF2011-6-B35)

通信作者：邓喀中(1957-)，男，四川资中人，教授，博士生导师，从事3S技术集成与数据处理方面研究；电话：15062195909；E-mail：cumtdkz@sina.com