1 引 言

摄影测量技术被广泛地应用于从二维航空/航天影像自动获取三维空间目标的几何信息、辐射信息和语义信息，进而构建相应的三维模型。其中关键技术是影像匹配。传统的影像匹配受限制于航空传感器的立体成像能力，以“单基线”的影像匹配为主^[1]。自2000年以来，随着国内外新型航空数字传感器(ADS40/80、DMC、UltraCamD/UltraCamX、SWDC等)相继出现，各种基于低空飞机和无人机平台的数字成像系统的快速发展，获得覆盖同一地区大重叠度的数字影像变得越来越容易，为影像匹配提供大量冗余信息。与单立体影像匹配相比，多视影像匹配主要有两方面优势：第一，利用多视影像的冗余信息有利于解决影像匹配中由于重复纹理、断裂特征等引起的错误匹配，提高匹配的可靠性；第二，利用多视影像进行匹配，最大限度地减少信息盲区，有利于解决影像匹配中的遮挡问题^{[2, 3]}。

针对航空数字传感器获取的城区多视影像特点，本文深入研究如何通过多视影像选择性匹配有效地解决匹配中的遮挡问题。在此前提下，建立像方和物方相结合的多视影像匹配模型，实现多匹配基元集成的多视影像密集匹配。基于多影像、多匹配基元、高精度影像匹配算法的半自动地物提取和测量是目前航空数字影像处理的主要研究方向之一^[4]。匹配基元的存在形式有两种，一种是存在于像方，影像上的特征(包括点、线、区域等)或者规则分布的格网点。常用的基于特征(影像上的特征点、特征线)的影像匹配可靠性较高^{[5, 6, 7, 8, 9]}，有利于后续的三维表面模型重建，但本质上属于稀疏的影像匹配(sparse matching)，匹配结果不够密集。现有方法在基于特征匹配后，通过区域增长或者利用影像上规则分布的格网点实现密集匹配^[2]，但是这种完全采用像方匹配基元的方法受限制于参考影像上的成像信息。如果物方的某一区域在参考影像上由于遮挡没有成像，会造成该区域的匹配漏洞。尽管可以依次选用每一张影像都作为参考影像进行匹配^[10]，但会增加计算量。另一种是存在于物方,抽象化的物方“地面元”。文献^[11]对“地面元”的概念进行讨论，文中将“地面元”(groundel,ground element)定义为在物方空间的一个单元，相似于影像空间的像元。为了形象化地面元，通常对物方作业范围平面化，并将平面划分为规则格网。地面元在影像匹配中有着重要的应用。传统的影像相关匹配方法，铅垂线轨迹法(vertical line locus,VLL)从物方地面元出发，地面元平面坐标(X,Y)已知，通过影像相关确定地面元高程值Z的同时确定影像上的同名像点^[12]。文献^{[13, 14, 15]}将该方法应用到多视影像匹配；文献^[16]认为多片最小二乘匹配中不同影像之间的仿射变形系数不是独立的，而是地面模型的反应。并从物方面元的思想出发，利用斜面元的法向量方向角，重新定义多片最小二乘匹配，大大减少未知数个数。这种完全基于地面元的匹配，匹配从物方出发，匹配方式更为灵活，不受影像的成像信息限制，但是放弃了影像上特征信息，不利于三维建模。

鉴于单一匹配基元存在的问题，本文提出一种多匹配基元集成的多视影像密集匹配方法。该方法利用规则格网划分的空间平面作为基础，集成像方特征点和物方面元两种匹配基元，充分利用多视影像上的特征信息和影像的成像信息，采用参考影像不固定的匹配策略，通过对多视影像进行选择性匹配，避免了遮挡对匹配产生的影响。

该方法以规则格网划分的空间平面作为匹配基础。首先通过移动平面到不同高程位置，利用平面元位置约束多视影像上特征点投影范围，实现特征点和平面元的同时匹配。然后结合VLL和基于高度遮挡检测方法对平面上规则分布的平面元进行匹配，加密初始匹配结果。图 1为多匹配基元集成的多视影像密集匹配方法总体流程图。具体包括空间格网平面确定、初始特征点匹配、平面元加密匹配3部分内容。

图 1 总体算法流程图 Fig. 1 Flow chart of the overall algorithm

图选项

在航天/航空摄影测量过程中，常取平面的方向垂直于物方铅垂线的方向^[17]，大小应包含所有参与匹配的影像在物方的覆盖范围XY∈(X_min~X_max,Y_min~Y_max)。利用规则格网将平面划分为小的平面元，使其在平面上沿平行于物方的XY轴方向规则分布，根据匹配要达到的精度或者根据生成DSM的分辨率来确定面元的大小，但不应该小于对应影像分辨率的大小，每个面元中心的平面坐标是固定的。

匹配过程中，可以抽象地将规则格网单元划分后的平面看作栅格化的DSM，平面上每一个平面元都看做DSM的一个像元，初始化平面元不具备任何高程值属性，最理想的结果是通过多视影像匹配使得平面上的所有平面元都被赋予正确的高程值，但是在实际匹配过程中，尽可能使平面上高程赋值的平面元具有密集性及均匀性。

首先对所有影像都进行特征点提取。确定影像覆盖区域地表的最大高程值Z_max和最小高程值Z_min，使平面在Z_min~Z_max范围内，按给定的步长大小，沿着物方的铅垂线方向从上到下移动。平面移动到每个高程值位置时，平面方程为Z=Z_i，初始化平面上每个面元内通过投影光束的数目number=0。反投影所有影像特征点与平面相交，分别统计平面上每个面元内通过投影光束的数目number，如果number＞T(T为阈值)，则认为通过该面元的不同影像上的特征点为候选同名像点。需要进一步选择候选同名像点所在的影像，进行灰度相关匹配，如果相关系数大于一定阈值，确认是同名点。同时，将该平面元进行高程赋值，不再参与后续的匹配，高程值为平面所处位置的高程值。该方法通过移动平面到不同高程值位置，利用平面元位置约束候选同名特征点，因此该方法也称作移动高程平面约束的多视影像特征点匹配。具体算法可参见参考文献^[18]。

特征点匹配后，平面上成功匹配的平面元被赋予正确高程值，如图 1中灰色格网单元所示。此时，平面可以看作初始的DSM，在此基础上，利用平面上规则分布的平面元进行密集匹配。如图 2中圆状点所在的平面元，已知其平面坐标(X,Y)，利用初始的DSM，计算其初始近似高程值Z₀，采用基于高度的遮挡检测和铅垂线约束相结合的方法对其进行匹配，确定其正确的高程值Z和影像上的同名像点。

图 2 特征点匹配后高程赋值格网平面以及待匹配平面元示意图 Fig. 2 The grid plane with height evaluation after matching and the plane elements to be matched

图选项

基于高度的遮挡检测方法(height-based method)是在真正射影像制作过程中提出的^[19]，用于消除由于影像上遮挡引起的真正射影像上的“鬼影”现象。本文将该方法引入到铅垂线约束的多视影像匹配中。

该方法的基本原理：判断某一物方点P在影像上是否可见，在其搜索路径上，从物方点向地底点方向依次比较该搜索路径上的物方点的高程和该物方点对应位置投影光束的高程值大小，如果该搜索路径上任何一点高程高于投影光束高度，则该物点在成像过程中被遮挡，没有成像。投影光束是从投影中心到地面上物点的光线，搜索路径是该投影光束在水平地面上的投影，如图 3所示。投影光束高程值可以用公式i×ΔZ+Z₀计算得到，。其中i为搜索路径上第i点，d为搜索路径上搜索步长，(X_S,Y_S,Z_S)、(X₀,Y₀,Z₀)为投影中心和待匹配点的物方坐标。

图 3 基于高度遮挡检测方法示意图 Fig. 3 The illustration of occlusion detection using height-based method

图选项

对平面上规则分布的任一平面元，首先判断是否已被高程赋值，如果没有，该平面元作为密集匹配待匹配基元。如图 4中以P点为中心的平面元，已知P点的平面坐标(X,Y)，利用初始DSM内插P点初始近似高程值Z₀。采用基于高度的遮挡检测算法对P点进行遮挡检测，判断P点在其中哪些影像上发生遮挡，哪些影像上成像。假设已成像的影像数目为m，如果m＜2，P点在所有待匹配影像都没有成像或者仅在其中一张影像上成像，不进行匹配；如果满足m≥2，选择已成像的影像，采用铅垂线约束对平面元P进行匹配，在确定其物方正确高程值的同时确定影像上的同名像点，同时得到不同影像上的区域遮挡图。铅垂线约束原理详见参考文献^{[12, 13, 14, 15]}。

图 4 密集匹配基本原理 Fig. 4 The principle of dense matching

图选项

图 4为密集匹配示意图，具体操作步骤如下：

(1) 已知P点平面坐标(X,Y)与近似高程Z₀，对所有待匹配的影像，都进行遮挡检测，判断该点未发生遮挡的影像，同时记录该点发生遮挡的影像。如果未发生遮挡影像的数目m≥2，选择未发生遮挡的影像进行后续匹配，执行下面步骤，否则对下一个点进行判断。

(2) 当m＞2，选取所有待匹配影像中投影中心离该格网单元最近的影像作为参考影像，即min(D₁,D₂,D₃，…，D_m)最小值对应的影像作为参考影像，其中D_i=，(XS_i,YS_i)为第i张影像投影中心坐标；当m=2时，不需要确定参考影像，直接进行匹配。

(3) 根据初始的DSM确定过P点铅垂线上的高程搜索范围Z_min~Z_max，确定搜索步距ΔZ。由P点的平面坐标(X,Y)与高程Z_i=Z_min+i×ΔZ(i=0,1,2,…)，分别计算其在m个影像上对应的像平面坐标(x_i^k,y^_ik)k=1,2,…，m。

(4) 固定以平面元P为中心的物方窗口，采用物方到像方的相关窗口变化策略，计算它在每个影像上的相关窗口。

(5) 分别计算参考影像上相关窗口和其他影像上相关窗口之间的灰度相关系数ρ。

(6) 将i值增加1，重复步骤(3)～(5)，得到ρ₀、ρ₁、ρ₂、…、ρ_n，算计ρ_h=max{ρ₀,ρ₁,ρ₂,…,ρ_n}，如果ρ_h大于一定的阈值，则ρ_h对应高程值Z_h=Z_min+h×ΔZ是平面元P的高程值，该高程值对应的影像上的点为同名像点。

(7) 将平面上所有待匹配的平面元都重复上面步骤(1)～(6)，逐面元进行匹配。

试验数据选取某城区同一条航带上的4张UCX航空数码影像，如图 5，像元大小为7.2 μm，相应地面分辨率0.049 m，航线重叠达80%以上。影像大小为14 430像素×9420像素，覆盖区域最大高差约为120 m，覆盖区域面积约为1025 m²×1069 m²。利用VirtuoZo进行自动空三加密，得到每张影像精确的内外方位元素，用于物方多视影像匹配。影像中包含高大建筑物，其对周围地物产生不同形式的遮挡。

图 5 试验影像区域 Fig. 5 Image areas of matching experiments

图选项

首先确定影像区域的地表覆盖范围，作为格网平面的范围；平面共划分为3207×2844个小的平面元，大小为0.3 m×0.3 m；确定平面沿着物方铅垂线移动的高程范围为3～130 m，高程移动步长取1 m。在此基础上，实现基于像方特征点和物方面元相结合的多视影像密集匹配。

本文采用Forstner算子提取所有影像的特征点，图 5中影像从左到右，依次提取的特征点数目为235 347、306 595、308 347、236 805。平面移动到每一个高程位置，分别对number≥2的平面元内通过的特征点进行灰度相似性约束匹配，灰度相似性测度阈值设置0.85。特征点匹配后得到同名像点对共为99 447。图 6为特征点匹配后高程赋值的格网平面，每一个赋值的平面元都对应一组同名像点。

图 6 特征点匹配后高程赋值的格网平面 Fig. 6 The grid plane with height evaluation after feature point matching

图选项

特征点匹配后，将高程赋值的格网平面看做初始DSM。在此基础上，对平面上小面元进行规则窗口划分，窗口大小为5×5，窗口中心的平面元作为密集匹配待匹配基元，数目共为158 989。图 7为密集匹配得到的不同影像上的同名像点，图中不同颜色点表示选择不同影像进行匹配结果，其中不同的影像组合成功匹配同名点对的数目小于100的全部用黑色表示。从图中色彩可以看出，根据影像的重叠度及遮挡检测结果，密集匹配过程中影像的选择具有一定的规律性：在影像4度重叠区域，如果4张影像全部成像则都参与匹配，得到同名像点用红色点表示，否则根据遮挡检测结果选择其中的3张或者两张影像进行匹配；同样在影像3度重叠区域，如果3张影像全部成像则都参与匹配，否则根据遮挡检测结果选择其中的两张影像进行匹配，如图中粉色、黄色、白色点所示区域，它与图 8区域遮挡图上发生遮挡的位置正好相反。即在3度重叠的影像区域，如果选择其中两张影像进行匹配，那么在另一张影像上该区域一定发生了遮挡。

图 7 密集匹配得到的不同影像上的同名像点 Fig. 7 The corresponding image points on different images after dense matching

图选项

图 8 不同影像上的区域遮挡图 Fig. 8 Occlusion maps in the different images

图选项

图 8为平面元在对应的不同影像上的区域遮挡图，其中底图为经过特征点匹配和平面元匹配后高程赋值的格网平面，细节见图 9。红色标记为物方平面元所在位置在该影像成像过程中由于遮挡没有成像，图中蓝色十字丝标记为该影像的投影中心在平面上的投影。从图中可以看出影像上遮挡产生的趋势：以投影中心为中心，向投影光线方向倾斜，越远离地底点位置，越高大的建筑物产生的遮挡越严重。经密集匹配后，成功匹配平面元的数目为61 787。

图 9 密集匹配后高程赋值的格网平面 Fig. 9 The grid plane with height evaluation after dense matching

图选项

通过特征点匹配和物方平面元匹配后，得到高程赋值的格网平面，如图 9所示。与图 6相同区域(如图中红色方框内)进行比较。可以看出，密集匹配在一定程度上加密了初始特征点匹配结果。

图 10为部分放大结果图，图中两栋建筑物在不同影像上表现较大投影差，对其附近地面及其自身产生不同形式的遮挡，图中右下边的建筑物部分屋顶在影像上表现为相似纹理结构。图 10(a)为匹配得到的影像上的同名像点，图 10(b)为不同影像上的区域遮挡图，红色点表示发生遮挡的平面元。从中可以看出，本文方法通过对多视影像进行选择性匹配，可以有效地避免影像上存在遮挡对匹配产生的影响，取得可靠的匹配结果，有效解决多视影像匹配中的遮挡问题。

图 10 部分放大结果图 Fig. 10 Parts of zoom results

图选项

本文在Matlab 2011b平台下编程实现多视影像匹配算法。试验所采用的硬件设备为普通台式机，其基本配置为：CPU3.10 Ghz,Inter(R) Core(TM) i5-2400；硬盘500 GB；内存3.47 GB；显卡NVIDIA GeForce 605；操作系统Microsoft Windows XP。匹配共包括两个阶段，共耗时1.79 h。其中特征点匹配过程中，平面移动效率为每分钟1.46次。该过程原理简单，易于实现，几何约束性强，不需要金字塔影像约束或者金字塔分级匹配，耗时主要是由于平面移动过程中，特征点反复投影导致计算量增多。该过程运行时间主要取决于影像大小和平面移动的高程范围和步长；物方面元匹配过程中，面元处理效率为每秒133个。该过程根据初始DSM内插出近似高程值，缩短高程搜索范围，并根据遮挡检测结果选择未发生遮挡的影像进行匹配，在一定程度上减少了计算量，提高算法的运行效率。

最后对匹配得到的不同影像上的同名像点，利用多基线前方交会计算匹配点的三维坐标^[20]，再利用离散点生成的结果DSM，如图 11所示。通过对结果DSM进行视觉检测可以看出，DSM较好地描述了影像覆盖区地表真实信息。此外，本文通过定量的评价结果DSM的精度来间接评价匹配结果的可靠性及精度。

图 11 离散点生成DSM Fig. 11 Generation of DSM by discrete points

图选项

在相同定向参数条件下，通过VirtuoZo对影像进行立体观测，共选169个检查点，加上两个GPS的控制点，共171个检查点。这些点在影像上均匀分布，分别位于建筑物屋顶、角点和地面。通过比较结果DSM上检查点处高程值Z和实际量测高程值Z′的差值。通过统计计算均方根误差(root mean square error,RMSE)，误差均值和最大误差max来评价DSM的精度。

点击浏览公式

表 1为精度评定结果，从中可以看出，本文的多视影像匹配方法取得高的匹配精度，从而验证了本文提出的多视影像匹配方法的可靠性。

本文针对多视影像匹配存在的优势，鉴于单一匹配基元存在的问题，利用规则格网划分的空间平面作为基础，对像方特征点和物方平面元两种不同形式的匹配基元进行整合。并建立了多视影像的物方匹配模型，该模型能有效避免遮挡对匹配的影响，实现多视影像密集匹配。论文选用同一条航带上的4张序列UCX航空数码影像进行匹配试验，并采用VirtuoZo立体量测结果对匹配得到的DSM进行精度评定，验证本文方法的可靠性，得到结论如下：

(1) 通过多视影像上特征点投影范围和平面元位置相互制约，实现特征点和平面元的同时匹配。为后续匹配提供初始可靠的DSM，提高后续匹配的可靠性及匹配效率；

(2) 将遮挡检测引入到铅垂线约束的多视影像匹配方法中，根据遮挡检测结果选择未发生遮挡的影像进行匹配，有效地避开发生遮挡的影像对匹配产生的影响，提高了匹配的可靠性；

(3) 在像方特征点匹配后，采用物方平面上规则分布的平面元代替像方的格网点进行加密匹配，使得参考影像选择方式更灵活，同时可以充分利用多视影像的成像信息，避免匹配受限于单一影像的成像信息；

(4) 进一步工作将研究特征线的匹配，在现有结果基础上，实现同名点和空间几何约束下的多视影像特征线匹配，将传统像方空间线特征相似性约束转移到物方空间，根据物方空间中直线的几何关系(直线间角度、距离、重叠度等)以及直线与平面关系进一步确定候选同名直线，同时考虑采用分区域匹配优化算法运行效率。