《测绘学报》
构建与学术的桥梁 拉近与权威的距离
倾斜影像自适应初始物方面元优化匹配算法
张春森1
, 张奇源1, 郭丙轩2, 薛万唱2
1. 西安科技大学测绘科学与技术学院, 陕西 西安 710054;
2. 武汉大学测绘遥感信息工程国家重点实验室, 湖北 武汉 430079
收稿日期：2019-01-28；修回日期：2019-08-02
基金项目：国家自然科学基金重大研究计划（91638203）；陕西省自然科学基金（2018JM5103）；国家自然科学基金（91738301）
第一作者简介：张春森(1963-), 男, 博士, 教授, 研究方向为数字摄影测量计算机视觉与遥感应用。E-mail:zhchunsen@aliyun.com
摘要：采用基于物方面元的最小二乘影像匹配方法匹配倾斜影像时，常出现深度不连续或高差较大区域影像连接点度数低或空三点过少问题。针对此问题，本文提出一种基于自适应初始物方面元的倾斜影像匹配算法。算法利用倾斜影像已有的初始内外方位元素及匹配过程中产生的物点信息，采用多片前方交会和物方差分的方法自适应计算物方面元的高程及法向量方向角初值，进而解决采用物方面元最小二乘影像匹配方法匹配地物高差较大区域的倾斜影像时，因初值不准导致在像方匹配同名点困难的问题。分别采用本文算法和物方面元初始状态为水平面元的最小二乘影像匹配方法对两组倾斜影像进行对比匹配验证。试验结果证明了本文算法的有效性。
关键词：倾斜影像匹配 最小二乘 单应矩阵 自适应物方面元
Oblique image matching algorithm based on adaptive initial object patch
ZHANG Chunsen1
, ZHANG Qiyuan1, GUO Bingxuan2, XUE Wanchang2
1. College of Geomatics, Xi'an University of Science and Technology, Xi'an 710054, China;
2. State Key Laboratory of Information Engineering in Surveying, Mapping&Remote Sensing, Wuhan University, Wuhan 430079, China
Foundation support: The Major Research Plan of the National Natural Science Foundation of China(No.91638203); The Natural Science Foundation of Shaanxi Province (No. 2018JM5103); The National Natural Science Foundation of China (No.91738301)
First author: ZHANG Chunsen(1963—), male, PhD, professor, majors in digital photogrammetry computer vision and remote sensing application. E-mail:zhchunsen@aliyun.com.
Abstract: In order to solve the problem of low degree of image connection points and too few aero triangulation points with large depth discontinuity or difference when using the least square image matching method based on object patch to match oblique images, an algorithm of oblique image matching based on adaptive initial object patch is proposed. The algorithm uses the initial values of interior and exterior orientation elements of the oblique image and the object points generated in the matching process, and adaptively calculates the elevation of the object patch and the initial value of the normal vector direction angle by the method of multi-patch forward intersection and object variance partition. Then it solves the problem that the low accuracy initial value of the oblique image with the least squares image matching method of the object patch is not easily matched to corresponding points when matching the oblique image. Finally, two groups of oblique images are matched by the proposed method and the least square image matching method, in which the initial state of the object patch is horizontal object patch, respectively, to verify the effectiveness of the proposed method.
Key words: oblique image matching least-square homograph matrix adaptive initial object patch
对航空影像进行精确匹配是采用摄影测量技术获取地物的空间位置信息[1]和对地形地物进行三维重建的关键步骤[2]。在过去的几十年中，计算机视觉和摄影测量领域的学者已对影像匹配作了大量卓有成效的研究。其中，基于像方的匹配从基于灰度的影像匹配方法发展到了准确性和可靠性更高的基于特征的影像匹配[3]，如基于特征松弛法的航空影像匹配[4]、结合全局信息的SIFT特征匹配算法[5]和基于特征点引导的多视影像择优匹配方法[6]等。基于物方地面元的影像匹配是在传统铅垂线轨迹法(vertical line locus，VLL)的基础上发展而来，如单立体VLL法[7]和基于多视影像匹配的MVLL法[8-10]等。这些影像匹配方法在解决不同类型影像匹配的同时，极大地促进了数字摄影测量向自动化和智能化方向的发展，降低了技术应用的难度并推动了航空摄影测量行业的发展。
此外，国内学者在经典的基于共线条件约束的多片最小二乘影像匹配算法(multi photo geometrically constrained matching，MPGC)[11]的基础上，引入物方面元的概念，消除了各影像间独立的仿射变形，建立了新的最小二乘影像匹配模型，提出了基于物方面元的最小二乘影像匹配方法[12]。该方法降低了MPGC算法的自由度并在影像匹配的同时可生成数字表面模型(digital surface model，DSM)，故又被称为生成DSM的影像匹配方法[13]。文献[12]在用该方法进行影像匹配时将基准影像上的像点p0投影到初始状态为水平，概略高程为Zc的物方面元上得到物点P(Xc, Yc, Zc)的近似坐标，之后将以物点P为中心的物方面元patch反投影至各搜索影像上，采用相关系数法匹配到同名点后，通过多片前方交会得到物点的初值，最后将物点及物方面元法向量方向角的初值代入误差方程，迭代求解物点和像点的精确坐标。
基于物方面元的最小二乘影像匹配方法在用于传统框幅式航空影像匹配时，在匹配精度、匹配效率及算法的多用途方面都取得了很好的效果。但将该方法用于倾斜影像匹配时，由于倾斜影像倾角大、物方条件复杂(存在较多深度不连续和高差较大区域)，导致将基准影像上的像点p0投影到初始状态为水平，概略高程为Zc的物方面元上得到物点坐标的误差较大，在将物点反投影至各搜索影像上后采用相关系数法匹配到的像点残差较大或匹配失败，导致空三处理后在一些深度不连续和高差较大区域空三点过少。而倾斜摄影三维重建对影像连接点的度数及空三点的密度要求较高，当以空三点为基础经密集匹配、构网、纹理映射后模型上会存在空洞。故该算法在用于倾斜影像匹配时需做进一步改进和优化。
针对上述问题，本文在基于物方面元最小二乘影像匹配方法的基础上，提出一种基于自适应初始物方面元的倾斜影像匹配算法。算法利用影像的初始内外方位元素，采用物方差分的方法，根据匹配点位的变化自适应计算物方面元的高程及法向量方向角的初值，以保证在像方可以匹配到准确可靠的同名点，从而提高影像连接点的度数，扩展物方面元最小二乘影像匹配方法的应用场景。
1 基于自适应初始物方面元倾斜影像匹配算法
在将物方面元最小二乘影像匹配方法应用于倾斜影像匹配时，为了解决在深度不连续或高差较大区域像方匹配同名像点困难的问题，本文通过改进匹配初值的计算方式，通过提高初值的准确性以保证在像方采用相关系数法匹配时，匹配到准确可靠的同名像点。本文方法计算物方面元高程初值，通过一次投影两次前方交会的策略完成；计算物方面元的法向量方向角α、β的初值是基于物方局部表面连续、物方面元是局部表面切平面的思想，利用匹配过程中产生的物方信息采用物方差分计算。本文匹配算法的流程如图 1所示。
图 1 本文倾斜影像匹配方法Fig. 1 The proposed oblique image matching method
图选项
1.1 自适应计算物方面元的高程初值
本文自适应计算物方面元高程初值的步骤为：①由影像的初始内外方位元素计算基准影像与各搜索影像间的相对单应矩阵Hi(i=1, 2, …, m－1)；②由相对单应矩阵Hi估算出基准影像上像点p0在各搜索影像上同名像点pi的近似坐标；③采用多片前方交会计算物方面元的高程近似值Zc；④由物点信息采用物方差分法计算物方面元法向量方向角α、β的初值；⑤通过将像点p0投影到高程为Zc，法向量方向角为α、β的物方面元后得到物方面元中心的物点P(Xc, Yc, Zc)；⑥再将物方面元反投影至各搜索影像得到像方区域wi；⑦在像方wi区域内采用相关系数法匹配到同名点pi′的像点坐标；⑧由像点p0, p′ 1, …, p′ m－1采用多片前方交会计算物点P(Xc, Yc, Zc)初值。
在以上步骤中，由H矩阵估计同名像点pi的近似坐标后采用多片前方交会计算面元的近似高程Zc是自适应计算物方面元高程初值的关键，可以避免因Zc的误差较大导致将像点p0投影到物方面元后物点P不准使得将物方面元反投影至搜索影像后，在像方匹配同名点困难的问题。由影像的初始内外方位元素计算影像间的相对单应矩阵H及由H矩阵估计同名像点近似坐标过程如下：
基准影像和搜索影像的内外方位元素和投影矩阵P (反映中心投影像方到物方的映射关系)之间的关系为
(1)
式中，P1、P2、K1、K2、R1、R2、C1、C2分别为基准影像和搜索影像的投影矩阵、内参数矩阵、旋转矩阵和摄影中心的坐标(影像外参的线元素)。令C1、C2的关系为ΔC = C2－ C1，假设在物方平面上有一物点为X，则对基准影像和搜索配影像其像点x1、x2与物方点X的关为
(2)
将基准影像和搜索影像纳入到以基准影像的像空间坐标系为像空间辅助坐标系进行影像匹配。此时式(2)可写为
(3)
联立式(2)和式(3)得到基准影像上像点x1和搜索影像上像点x2的关系为
(4)
即得基准影像与搜索影像间相对单应矩阵H的表达式为
(5)
由式(5)可利用影像的初始内外方位元素计算基准影像与各搜索影像间的相对单应矩阵。在影像匹配时，由基准影像上像点p0的变化估计同名点坐标，最后根据步骤③至步骤⑧自适应计算物方面元的初始高程。
1.2 物方差分法计算物方面元的法向量方向角初值
物方面元的几何意义是地物局部表面的切平面，当匹配像点对应的物点位于水平面上时，以水平面元作为物方面元的初始状态进行倾斜影像匹配可取得较好的匹配效果。但当匹配像点对应的物点位于高差较大的立面(如建筑物墙面)上时，以水平面元作为物方面元的初始状态进行倾斜影像匹配就会因投影误差较大导致在像方匹配不到准确可靠的同名点。在倾斜影像匹配时，根据匹配点的变化为误差方程给出较为准确的物方面元法向量方向角的初值。本文根据物方局部表面连续的原理，采用物方差分法自适应计算物方面元法向量方向角α、β的初值以改善影像匹配的效果。
根据物方局部表面近似连续的条件，由匹配过程中产生的物方点对局部地形进曲面拟合
(6)
然后计算局部曲面f(x, y)在Pi(Xi, Yi, Zi)点处的切平面的法向量
(7)
在求得物方局部曲面在点Pi处切平面的法向量后，采用差分代替微分计算偏导数的近似值，则以物方点Pi(Xi, Yi, Zi)为中心的局部曲面元(也即物方面元)的法向量n =[n1 n2 1]为
(8)
式中，(Xi－1, Yj, Zi－1, j)、(Xi+1, Yj, Zi+1, j)和(Xi, Yj-1, Zi, j-1)、(Xi, Yj+1, Zi, j+1)分别表示在物方X和Y方向上相邻点的坐标。
将法向量n归一化后，求得物方面元法向量方向角
(9)
最后，将上文中求得的物方面元中心点P(Xc, Yc, Zc)的初值及法向量方向角α、β的初值代入误差方程[12](式(10))进行最小二乘平差解算
(10)
式中，g(xi, yi)为待匹配点处影像的灰度；(dxi, dyi)为像点坐标的改正值。方程的未知数为物方面元的高程Zc，物方面元法向量方向角α、β和基准影像与第i幅搜索影像间的辐射变形参数h0i, h1i。其中，像点坐标的改正值与方程未知数间的关系为
(11)
限于篇幅，式(11)的详细推导过程可参见文献[12]，本文不再赘述。
2 试验及精度分析2.1 试验数据
为验证本文算法的有效性，选择两组倾斜影像数据进行试验。第1组是地物高差较大区域的160幅倾斜影像，第2组是地形较为平坦地区的150幅倾斜影像。两组影像数据的详细信息见表 1。基于Win7 64位操作系统，本文在VS 2010中配置GodWork-AT软件开发环境，采用VC++语言编写程序进行算法验证，运行平台为Intel Core 3.3 GHz处理器，16 GB内存的计算机。
表 1 试验数据相关信息Tab. 1 Experimental data information
试验数据影像数量/幅影像大小/像素焦距/mm航向重叠率/(%)旁向重叠率/(%)航高/m地表近似高程/m
下视侧视
第1组32×57592×53043550806530025
第2组30×55472×364810.410.4806518032
表选项
2.2 地物高差较大区域的倾斜影像匹配试验
分别采用文献[12]物方面元初始状态为水平面元的最小二影像匹配方法(以下简称为“初始水平面元法”)和本文提出的方法，对第1组倾斜影像进行匹配试验。试验中通过改变相关窗口的尺寸大小并统计试验结果，进而研究相关窗口尺寸对影像匹配结果的影响。本文随机选取20个像对，每组像对400对同名点在不同窗口尺寸下匹配的结果如表 2所示。为验证影像初始外方位元素的精度对本文方法匹配结果的影响，在用该方法对第1组倾斜影像进行匹配、光束法平差后，统计了有代表性的两幅影像的外方位元素改正数及每幅影像像点匹配时的平均迭代次数(表 3)。为从整体角度验证采用本文方法对倾斜影像匹配后生成空三点的质量，在空三处理后统计了影像连接点的平均度数及有效物点个数等指标(表 4)。为了比较初值准确性对匹配效率的影响，在表 5中给出了部分像点匹配时的初值及平差后的值。两种匹配方法得到连接点的物点(空三点)如图 2所示。以空三点为基础，通过密集匹配、构网、纹理映射后重建的三维模型如图 3所示。
表 2 窗口尺寸对匹配结果的影响Tab. 2 The effect of window size on matching results
方法窗口尺寸/像素平均相关系数(ρ)平均迭代次数
初始水平面元法9×90.60530.210
15×150.63229.812
19×190.70926.302
25×250.73225.656
31×310.74625.036
37×370.71226.122
43×430.65728.672
本文方法9×90.69327.866
15×150.74223.681
19×190.86321.323
25×250.89218.765
31×310.91815.618
37×370.83222.442
43×430.76821.783
表选项
表 3 影像外方位元素初值对匹配结果的影响Tab. 3 The effect of initial value of elements of exterior orientation on matching result
影像名ΔX/mΔY/mΔZ/mΔΩ/(°)Δρ/(°)Δκ/(°)平均迭代次数
DSC1521-0.2490.3500.856-0.3570.5405-0.43615.325
DSC2153-0.1060.1600.437-0.2030.213-0.35714.863
表选项
表 4 匹配结果及匹配效率对比Tab. 4 Matching results and matching efficiency comparison
影像名平差后像点残差均值/像素影像连接点的平均度数有效像点数有效物点数匹配时间/min
初始水平面元法0.793 23.22 768 365990 13817.3
本文方法0.564 85.23 698 6311 148 96313.2
表选项
表 5 初值的准确性对匹配结果的影响Tab. 5 The effect of the accuracy of initial value on matching result
物点名Zc/m
(α、β)/(°)迭代次数
初值平差后的值初值平差后的值
初始面元法本文方法初始面元法本文方法
初始面元法本文方法初始面元法本文方法初始面元法本文方法
T1033712527.827.627.5
(0, 90)(5, 73)(7, 69)(6, 70)1213
T1033722526.626.126.1(0, 90)(6, 89)(9, 85)(9, 87)1110
T1033732531.730.930.6(0, 90)(16, 46)(15, 51)(19, 51)1413
T1033742536.6/38.6(0, 90)(7, 73)/(11, 79)/12
T1033752532.130.831.2(0, 90)(21, 56)(28, 81)(30, 78)1614
T1033762528.327.126.9(0, 90)(11, 67)(12, 63)(12, 63)1111
T1033772536.637.337.7(0, 90)(45, 60)(40, 71)(43, 68)2314
T1033782535.73636.4(0, 90)(35, 55)(34, 55)(39, 51)2118
T1033792541.640.942.8(0, 90)(8, 65)(14, 65)(11, 63)2817
T1033802548.2/47.6(0, 90)(10, 70)/(8, 68)/14
T1033812551.2/53.4(0, 90)(30, 45)/(32, 41)/16
T1033822533.734.332.9(0, 90)(55, 60)(52, 68)(54, 66)1815
T1033832531.632.432.9(0, 90)(23, 75)(25, 76)(24, 79)1517
T1033842537.839.938.7(0, 90)(40, 55)(40, 58)(36, 51)2917
注：匹配时所开窗口大小为25×25像素，/表示匹配失败。
表选项
图 2 空三点云对比Fig. 2 Comparison of point clouds in aerotriangulation
图选项
图 3 以空三点云为基础重建三维模型对比Fig. 3 Comparison of three-dimensional models reconstructed on the basis of aerotriangulation pointclouds
图选项
2.3 地形平坦地区的倾斜影像匹配试验
为了进一步验证本文方法对地形平坦地区倾斜影像匹配的效果，分别采用初始水平面元法和本文方法，对第2组倾斜影像进行匹配试验。在影像匹配光束法平差后，统计了像点残差均值及影像连接点的平均度数、有效物点个数及匹配用时等指标(表 6)。在影像匹配、光束法平差获得空三点后做了三维重建，用于从视觉上对比两种匹配方法在空三处理(图 4)及三维重建(图 5)方面的效果。
表 6 平坦地区倾斜影像的匹配结果及匹配效率对比Tab. 6 Matching results and matching efficiency comparison of smooth area oblique images
影像名平差后像点残差均值/像素影像连接点的平均度数有效像点数有效物点数匹配时间/min
初始水平面元法0.5725.22 106 015796 43614.6
本文方法0.5835.12 187 015802 55315.7
表选项
图 4 空三点云对比Fig. 4 Comparison of pointclouds in aerotriangulation
图选项
图 5 以空三点云为基础重建三维模型对比Fig. 5 Comparison of three-dimensional models reconstructed on the basis of aerotriangulation pointclouds
图选项
2.4 试验结果分析
试验结果分析如下。
(1) 像方所开窗口的尺寸对匹配结果的影响。由于倾斜影像存在较大的倾角，当基准影像发生变化时，将基准影像上以像点p0为中心的窗口区域投影后所得到的物方区域(即物方面元的大小)会发生较大变化。本文基于物方面元最小二乘匹配方法匹配时是在像方区域采用相关系数法匹配得到较为准确的同名像点初值。故分析以基准影像上像点p0为中心所开窗口大小间接研究面元大小对匹配结果的影响。从表 2中可以看出，窗口尺寸在一定范围内时(小于31×31像素)，其与匹配到同名点的概率呈正相关。随着像方所开窗口的变大，最大相关系数随之变大，即匹配到同名点的准确度在提高。匹配到准确的同名点后可计算出物点坐标和α、β更为准确的初值，迭代次数减少。但窗口过大时，窗口内的地物及其纹理信息会变复杂，导致窗口在给定范围内滑动时求得的最大相关系数变小，匹配像点的准确度降低，故窗口尺寸的选择对匹配结果的影响较大。试验中发现，在某一临界值下，如本次试验中窗口尺寸取31×31时，可以更好地兼顾匹配的准确度和匹配效率。
(2) 影像初始内外方位元素的精度对匹配结果的影响。影像初始内外方位元素越精确，采用相对单应矩阵H估计得像点准确性越高。由同名像点的近似坐标采用多片前方交会得到面元中心的高程值Zc越准确，当把物方面元投影到搜索影像后，采用相关系数法可匹配到更为准确的同名像点，进而可提高匹配精度。从表 3两幅有代表性的影像平差后外方位元素改正数及每幅影像匹配像点时的平均迭代次数可以看到影像DSC2153的初始外方位元素较影像DSC1521的外方位元素较为准确，故影像匹配时平差解算的迭代次数较少。
(3) 初值的准确性对匹配结果的影响。从表 4初值的准确性对匹配结果的影响中可以看出，本文方法相比初始水平面元法可以给出更为准确的初值，故匹配时的迭代次数较少。另外，当Zc的初值偏差过大时，造成在像方的反投影区域内匹配不到同名像点，匹配失败，进而导致空三点较少。如图 2(a)中红色椭圆区域，由于初值不准确使得匹配到的同名点少、连接点度数低，导致在空三处理后空三点稀少的现象。在以空三点为基础进行三维重建时由于物点过少(图 3(a)中红色椭圆区域)，导致所构三角网缺乏细节层次的信息故模型上会存在空洞。
(4) 对地形平坦地区倾斜影像匹配的效果分析。第2组影像所在区域的地表高程约32 m，在匹配试验时初始水平面元法Zc的初值取32 m，本文方法的初值算法自动计算得到。从表 6影像连接点的平均度数及有效物点个数等指标可看出，在匹配高差变化不大区域的倾斜影像时，本文方法与初始水平面元法的匹配效果相近。由于本文方法计算初值的算法复杂度较高，故匹配用时多于初始水平面元法。另外从图 4可以看出，两种方法匹配、平差后得到空三点的密度相近。
3 结论
本文针对基于物方面元最小二乘影像匹配应用于倾斜影像匹配，在深度不连续或高差较大区域影像连接点度数低和空三点过少的问题，提出了一种基于自适应初始物方面元的倾斜影像匹配算法。算法采用多片前方交会和物方差分的方法自适应计算匹配初值，实现了将物方面元匹配方法应用到对倾斜影像的匹配中，拓展了物方面元匹配方法的应用场景。
该算法的核心思想是，通过提高计算物方面元高程及法向量方向角初值的准确度，使得在像方可以匹配到准确可靠的同名点，进而解决物方面元最小二乘匹配方法在匹配高差较大区域的倾斜影像时，因初值不准导致在像方匹配同名点困难的问题。本文通过对两组倾斜影像数据进行匹配对比试验，结果显示本文方法在匹配深度不连续或高差较大区域的倾斜影像时可以得到较好的匹配效果。同时，在匹配地形平坦地区的倾斜影像时，本文方法与文献[12]方法的匹配效果相近。
由于本文匹配算法是基于最小二乘算法属二次精匹配，即匹配精度依赖于初值的准确性，虽然文中创新性地给出了自适应计算匹配初值的方法，但未定量确定影像初始内外方位元素的精度与影像匹配精度间的关系，自适应计算匹配初值的算法复杂度有待进一步优化。此外，在倾斜影像匹配中存在的遮挡问题也将是后续研究的内容。
【引文格式】张春森, 张奇源, 郭丙轩, 等. 倾斜影像自适应初始物方面元优化匹配算法. 测绘学报，2020，49(1)：108-116. DOI: 10.11947/j.AGCS.2020.20190050
