CN105205808A - 基于多特征多约束的多视影像密集匹配融合方法及系统 - Google Patents

Abstract

本发明提供基于多特征多约束的多视影像密集匹配融合方法及系统，包括根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图；根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合；联合面特征和线特征进行点云优化，生成最终的点云。本发明技术方案能够自动选择合理的立体像对，充分利用多视信息使匹配结果更为精确可靠，能够获取全局意义下的最优多视融合结果，优化生成的点云更加精细。

Description

基于多特征多约束的多视影像密集匹配融合方法及系统

技术领域

本发明涉及多视影像密集匹配技术领域，尤其是涉及一种基于多特征多约束的多视影像密集匹配融合方法及系统。

背景技术

21世纪是一个信息化的时代，信息化和数字化对社会发展产生了深远的影响。在全世界正从数字地球向智慧地球发展的过程中，我国也在推动“数字中国”、“智慧城市”等的快速发展，在为政府部门、企事业单位和社会公众等提供公共管理、突发事件应急、科学决策等服务方面具有重要意义。影像密集匹配技术是主流的地物目标三维信息获取手段之一，将相机搭载于卫星、飞机、无人机、移动车辆等不同高度的遥感平台，通过同名光线对对相交的原理，快速获取目标的三维测绘信息，在城市三维建模、三维变化检测、DEM/DSM和真正射影像制作等领域有着非常广泛的应用前景。

密集匹配一直是摄影测量和计算机视觉领域经久不衰的研究热点。目前，已经出现了众多商业化的密集匹配软件，如Inpho、Intergraph、UltraMap、SimActive、VisualSFM、smart3D、AgisoftPhotoscan、PhotoModeler、StreetFactory、Autodesk123D等，能够实现卫星、航空、低空无人机、车载等多种情况下的立体影像三维重建。但是，这些商业软件仍然存在问题：匹配的点云不够密集(不是逐像素匹配)、三维模型表达不够精细、纹理贫乏或重复纹理区域误匹配现象严重、线状地物(如建筑物)边缘“毛刺”现象严重、特征描述过于单一(只描述点特征)等，无法完全满足高精度三维重建、地理国情监测、智慧城市等应用的需求。可以预见，未来的密集匹配技术将向准实时、精确、稳健、密集、多特征等方向发展，以满足现代社会发展的迫切需求。

发明内容

本发明主要是解决传统密集匹配技术生成点云不够密集、纹理贫乏以及纹理重复区域误匹配现象严重、线状地物边缘“毛刺”现象严重等问题；提供了一种基于多特征多约束的多视影像密集匹配与点云融合优化技术，能够自动选择最优的匹配模型，在密集匹配过程中充分利用多视约束条件，获取全局意义下的最优融合结果，并联合面特征和线特征对点云进行优化，从而得到更为精确可靠的密集匹配产品。

本发明技术方案提供一种基于多特征多约束的多视影像密集匹配融合方法，包括以下步骤：

步骤1，根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束；

步骤2，对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图，实现方式如下，

1)预测高程范围，包括根据特征匹配产生的物方点的高程，预测基准影像所覆盖测区的地表起伏，找到最大高程Z_max和最小高程Z_min；

2)计算高程方向的步距，包括从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面和最小高程面，分别有两个交点M、M'，在步骤1所得待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'，将最大最小高程差与投影点连线长度之间的比值，作为高程步距；

3)构造代价矩阵，包括以基准影像的像平面作为水平面，根据1)所得高程范围，根据2|)所得高程步距在高程方向划分间隔，建立一个三维网格矩阵，作为代价矩阵；所述代价矩阵中，每一个网格代表某像素在对应高程Z_i下的匹配代价；

4)在代价矩阵中，以互信息作为匹配代价，实现方式如下，

首先，计算基准影像的初始高程图，包括针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i，计算对应的物方点坐标(X_i,Y_i,Z_i)；将物方点同时投影到匹配模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像上像点(x_i,y_i)与各待匹配影像上像点(x_i',y_i')之间的相关系数并取平均，作为该像素在高程Z_i下的最终相关系数，在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程；

然后，将基准影像分别与各待匹配影像构成一个立体像对，根据初始高程图，对基准影像上的每个像素计算与各待匹配影像之间的互信息作为匹配代价，对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中；

最后，进行半全局密集匹配，包括根据代价矩阵得到基准影像上每个像素在高程范围内对应的匹配代价，在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，将匹配问题转化为动态规划问题求解，匹配结果为动态规划的最优路径，得到基准影像的高程图；

步骤3，根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合，实现方式如下，

以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网，三维格网的水平范围为测区的范围，三维格网的高程范围为地表的高程起伏，分别将不同视角的模型匹配结果投影到物方，每个独立的网格都包含零至多个物方点，物方点为步骤3针对各匹配模型所得高程图的相应点云，对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上二值化算子Census值的方差作为像方可见性约束，计算每个网格点的高程置信度如下，

conf＝-N/σ_Census

${\mathrm{\σ}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\‾}{Census})}^2/n}$

其中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能见该网格点的影像数目；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值；

所述全局能量函数如下，

$E=min\underset{i\∈G}{\Σ}(Cost\left(H_i\right)+\underset{q\∈N_i}{\mathrm{\Σ}}\mathrm{\σ}|H_i-H_q|/|I_i-I_q\left|\right)$

其中，E表示能量函数，作为融合结果的测度；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值；

步骤4，联合面特征和线特征进行点云优化，生成最终的点云。

而且，步骤4中，基于面特征进行优化的实现方式如下，

对融合后的点云在物方存在无高程的区域，记为“空洞”，首先根据“空洞”的平面位置找出覆盖该“空洞”的基准影像集合以及对应的匹配结果；对基准影像进行图像分割，提取面特征；将融合后的点云反投影到原始的基准影像上；统计每个分割块内的有效匹配点得到有效高程信息；根据有效高程信息，选择一个最优的高程平面参数，模拟分割块内的高程变化；最后，根据最优高程平面参数，计算分割块内“空洞”的高程，并投影到物方。

而且，步骤4中，基于线特征进行优化的实现方式如下，

首先在基准影像上提取直线边缘，根据基准影像与高程图一一对应的关系，在高程图上找到该边缘的位置，以该边缘为中心，建立一个缓冲区；统计缓冲区边缘两侧的高程分布情况，根据少数服从多数的原则，修正边缘两侧的错误高度，剔除毛刺。

本发明提供一种基于多特征多约束的多视影像密集匹配融合系统，包括以下模块：

匹配模型构建模块，用于根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束；

高程图提取模块，用于对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图，包括如下子模块，

高程范围预测子模块，用于根据特征匹配产生的物方点的高程，预测基准影像所覆盖测区的地表起伏，找到最大高程Z_max和最小高程Z_min；

步距计算子模块，用于从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面和最小高程面，分别有两个交点M、M'，在步骤1所得待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'，将最大最小高程差与投影点连线长度之间的比值，作为高程步距；

代价矩阵构造子模块，用于以基准影像的像平面作为水平面，根据高程范围预测子模块所得高程范围，根据步距计算子模块所得高程步距在高程方向划分间隔，建立一个三维网格矩阵，作为代价矩阵；所述代价矩阵中，每一个网格代表某像素在对应高程Z_i下的匹配代价；

匹配代价提取子模块，用于在代价矩阵中，以互信息作为匹配代价，实现方式如下，

首先，计算基准影像的初始高程图，包括针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i，计算对应的物方点坐标(X_i,Y_i,Z_i)；将物方点同时投影到匹配模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像上像点(x_i,y_i)与各待匹配影像上dioxide像点(x_i',y_i')之间的相关系数并取平均，作为该像素在高程Z_i下的最终相关系数，在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程；

然后，将基准影像分别与各待匹配影像构成一个立体像对，根据初始高程图，对基准影像上的每个像素计算与各待匹配影像之间的互信息作为匹配代价，对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中；

最后，进行半全局密集匹配，包括根据代价矩阵得到基准影像上每个像素在高程范围内对应的匹配代价，在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，将匹配问题转化为动态规划问题求解，匹配结果为动态规划的最优路径，得到基准影像的高程图；

融合模块，用于根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合，实现方式如下，

以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网，三维格网的水平范围为测区的范围，三维格网的高程范围为地表的高程起伏，分别将不同视角的模型匹配结果投影到物方，每个独立的网格都包含零至多个物方点，物方点为步骤3针对各匹配模型所得高程图的相应点云，对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上二值化算子Census值的方差作为像方可见性约束，计算每个网格点的高程置信度如下，

conf＝-N/σ_Census

${\mathrm{\σ}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\Σ}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\‾}{Census})}^2/n}$

其中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能见该网格点的影像数目；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值；

所述全局能量函数如下，

$E=min\underset{i\∈G}{\Σ}(Cost\left(H_i\right)+\underset{q\∈N_i}{\mathrm{\Σ}}\mathrm{\σ}|H_i-H_q|/|I_i-I_q\left|\right)$

其中，E表示能量函数，作为融合结果的测度；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值；

点云优化模块，用于联合面特征和线特征进行点云优化，生成最终的点云。

而且，点云优化模块中，基于面特征进行优化的实现方式如下，

对融合后的点云在物方存在无高程的区域，记为“空洞”，首先根据“空洞”的平面位置找出覆盖该“空洞”的基准影像集合以及对应的匹配结果；对基准影像进行图像分割，提取面特征；将融合后的点云反投影到原始的基准影像上；统计每个分割块内的有效匹配点得到有效高程信息；根据有效高程信息，选择一个最优的高程平面参数，模拟分割块内的高程变化；最后，根据最优高程平面参数，计算分割块内“空洞”的高程，并投影到物方。

而且，点云优化模块中，基于线特征进行优化的实现方式如下，

首先在基准影像上提取直线边缘，根据基准影像与高程图一一对应的关系，在高程图上找到该边缘的位置，以该边缘为中心，建立一个缓冲区；统计缓冲区边缘两侧的高程分布情况，根据少数服从多数的原则，修正边缘两侧的错误高度，剔除毛刺。

因此，本发明具有如下优点：能够自动选择合理的立体像对，无需人工干预，对航线的设计没有任何要求；在密集匹配过程中能够充分利用多视信息，减少误匹配，匹配结果更为精确可靠；能够获取全局意义下的最优多视融合结果，避免网格点高程的选择陷入局部极值解；能够联合面特征和线特征对点云进行优化，生成的点云更加精细。本发明在航天摄影测量、航空摄影测量、低空摄影测量和近景摄影测量领域有较好的应用前景。

附图说明

图1为本发明实施例的匹配模型自动选择技术路线示意图；

图2为本发明实施例的特征匹配约束和交会角约束技术路线示意图；

图3为本发明实施例的互信息计算技术路线示意图；

图4为本发明实施例的基于动态规划匹配示意图；

图5为本发明实施例的多方向动态规划示意图；

图6为本发明实施例的物方多视点云融合技术路线示意图；

图7为本发明实施例的物方“空洞”对应影像的示意图；

图8为本发明实施例的分割块内深度图分布示意图；

图9为本发明实施例的整体流程图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述。

具体实施时，本发明所提供方法可由本领域技术人员采用计算机软件方式实现自动运行流程。参见图9，实施例的流程包括以下步骤：

本发明可以分别处理不同遥感平台的影像，包括卫星影像、航空影像、低空影像和近景影像。其中，卫星影像为卫星切片影像，附带RPC参数；航空影像、低空影像以及近景影像一般为框幅式影像，附带内外方位元素。在构成测区影像集合时，要求测区影像尽量为同一遥感平台拍摄的数据，或者在分辨率、拍摄日期等方面不存在较大的差异。本发明提供的技术方案是，一种自动选择合理的匹配模型，根据多视条件，进行多视半全局密集匹配，获取全局意义下的多视融合结果，并结合多种特征对点云进行优化的方法，包括以下步骤：

步骤1.自动选择合理的立体匹配模型：根据多重约束，为每一张基准影像选择若干张待匹配影像，构成匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束。

匹配模型的选择则是密集匹配流程中必不可少的重要环节。选择合理的匹配模型不仅可以避免“盲匹配”等冗余计算，大大加快整个测区的密集匹配速度，还可以保证密集匹配点云的精度。本发明首先根据输入的测区影像集合(设测区影像集合中的影像数目为n)以及对应的内外方位元素，对影像进行灰度化、建立影像金字塔等预处理操作；在金字塔影像顶层，任意定义一张影像为基准影像，进行多约束的立体像对自动选择，包括采用基线约束、像平面法向量约束、特征匹配约束和交会角约束，在剩余影像集合中，自动找出合理的待匹配影像集合；将待匹配影像集合，存入对应的链码结构，方便后续匹配查询；然后在之前未定义为基准影像的影像集合中，定义一张新的基准影像，依次类推，直至遍历完所有影像为止，如图1所示。其中，N表示测区影像集合中的影像数目，测区影像集合中的影像序号记为i，i＝1,2,…,N；基准影像与根据多重约束选择的待匹配影像集合共同组成一个匹配模型。为了自动选择合理的匹配模型，同时兼顾匹配精度和交会精度，且尽可能地满足影像之间有比较大的重叠度，将综合采用基线约束、像平面法向量约束、特征匹配约束和交会角约束。

1)基线约束

通常情况下，相互重叠的影像往往是相邻影像，即影像之间的基线长度较短。因此可以将基线长度作为一个约束，从而大大缩小立体匹配像对集合的搜索范围。具体为，根据影像方位元素统计基准影像与其余影像之间的基线长度，选择基线长度最短的T_z张相邻影像，每张相邻影像分别与基准影像建立一个立体像对，构成一个候选立体匹配像对集合。具体实施时，本领域技术人员可自行预设T_z取值。

例如T_z＝2的情况，设航测过程中依次拍摄的影像为A、B、C、D、E、F、J，如果以E为待匹配的基准影像，则以E为基准影像，影像D、F则T_z张基线长度最短的相邻影像，D、F分别与基准影像E构成一个立体像对，构成根据基线约束的候选立体匹配像对集合。在大重叠度的情况下，同一个测区往往由数十甚至数百张影像覆盖，可以取较大的T_z值纳入更多的候选立体匹配像对，然后再根据后续约束进一步筛选最优立体像对集合。由于采用多重约束一步一步筛选最优立体像对集合，即使取较大的T_z值，也不会影响后续筛选结果。

2)像平面法向量约束

两张影像之间的重叠情况，除了受基线长度的影响，同时也受两张影像的像平面法向量之间夹角的影响。在基线长度合适的情况下，像平面法向量之间的夹角越大，则影像几何变形越大，且遮挡问题严重。上述问题会对后续的密集匹配造成严重的干扰，属于不合格的立体像对。本发明设计了一种像平面法向量约束条件，计算像平面的法向量，然后再根据式(1)，计算法向量之间的夹角；最后根据预设的夹角阈值T_a，在1)所得候选立体匹配像对集合中，进一步筛选更加合理的立体匹配像对。具体实施时，本领域技术人员可自行预设T_a取值。

θ＝arccos(N₁·N₂/|N₁||N₂|)(1)

式(1)中，θ表示像平面法向量之间的夹角；N₁,N₂代表候选立体匹配像对中第1、2幅影像的法向量；|N₁|、|N₂|代表候选立体匹配像对中第1、2幅影像法向量的模。

3)特征匹配约束和交会角约束

在满足一定基线条件和平面法向量条件的情况下，重叠度越高，则特征匹配的像点数目越多。而且特征匹配的同名点数目越多，说明影像纹理信息丰富，有利于后续的密集匹配；反之，则说明该影像质量较差，密集匹配结果不可靠。因此可以将特征匹配的点数作为判断影像之间重叠度和影像信噪比的一个依据，在多视的情况下可以将重叠度较小或信噪比较低的立体像对舍去，避免重复计算。此外，还需要检查影像之间的交会角，如果影像之间的交会角过小，虽然匹配质量很高，但是交会精度很低，无法生成高精度点云，因此交会角过小的立体像对也是不合格像对。针对上述情况，本发明设计基于特征匹配和交会角的约束条件，首先对基准影像和2)所得候选立体像对集合生成影像金字塔；针对候选立体像对集合中的每一个像对在金字塔顶层采用SIFT算子进行特征匹配，并剔除匹配点数较少的立体像对，具体实施时，本领域技术人员可自行预设匹配点数相应阈值取值；根据匹配点计算每个立体像对的平均交会角，去除交会角过小的立体像对，具体实施时，本领域技术人员可自行预设交会角相应阈值取值；在剩下的立体像对中找出匹配点数最多的若干像对，包括设定一个百分比阈值T_p，将其它像对的匹配点数与最大匹配点数进行比较，当所得百分比小于T_p则舍弃，筛选出最终的立体匹配像对集合，即可得到最终的待匹配影像集合，具体过程如图2所示。具体实施时，本领域技术人员可自行预设阈值T_p取值。

步骤2.多视条件下的半全局密集匹配：在单个匹配模型中，充分利用多视约束条件，进行半全局密集匹配，直接生成单个模型的密集匹配结果(即高程图)。

1)预测高程范围。具体实施时，可利用步骤1中产生的特征匹配点，采用前方交会的方法，生成稀疏的物方点。根据物方点的高程，预测基准影像所覆盖测区的地表起伏，找到一个最大高程Z_max和最小高程Z_min。

2)计算高程方向的步距。从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面(Z_max)和最小高程面(Z_min)，分别有两个交点M、M'。在步骤1所得最终的待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'。将最大最小高程差与投影点连线长度之间的比值，作为高程步距。

采用共线方程将物方点M投影到影像上：

$x-x_0=-f\frac{a_1(X-Xs)+b_1(Y-Ys)+c_1(Z-Zs)}{a_3(X-Xs)+b_3(Y-Ys)+c_3(Z-Zs)}$ (2)

$y-y_0=-f\frac{a_2(X-Xs)+b_2(Y-Ys)+c_2(Z-Zs)}{a_3(X-Xs)+b_3(Y-Ys)+c_3(Z-Zs)}$

式中，x、y表示投影点的像点坐标；X、Y、Z表示LiDAR点云的物方坐标；x₀、y₀、f表示相机的内方位元素；Xs、Ys、Zs表示外方位线元素；a₁～a₃、b₁～b₃、c₁～c₃表示旋转矩阵的9个元素。

高程步距的计算公式为：

Z_Step＝(Z_max-Z_min)/len_m,m'(3)

式中，Z_step表示高程步距；len_m,m'表示像点m、m'之间的长度(以像素为单位)。

3)构造代价矩阵。以基准影像的像平面作为水平面，以最大最小高程Z_max、Z_min作为高程范围，根据高程步距在高程方向划分间隔，建立一个三维网格矩阵，即代价矩阵。代价矩阵中，每一个网格代表像素i在对应高程Z_i下的匹配代价。

4)本发明以互信息作为匹配代价。首先计算基准影像的初始高程图：针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i(Z_max≤Z_i≤Z_min)，计算对应的物方点坐标(X_i,Y_i,Z_i)，如式(4)所示：

$\left(\begin{array}{c}X\\ Y\end{array}\right)={\left(\begin{array}{cc}{P}_0-{\mathrm{xP}}_8& P_1-{\mathrm{xP}}_9\\ P_4-{\mathrm{yP}}_8& P_5-{\mathrm{yP}}_9\end{array}\right)}^{-1}\left(\begin{array}{c}x\left(P_{10}Z+P_{11}\right)-\left(P_{2}Z+P_3\right)\\ y\left(P_{10}Z+P_{11}\right)-\left(P_{6}Z+P_7\right)\end{array}\right)---\left(4\right)$

式中，X、Y表示表示物方坐标；Z表示物方高程值；x、y表示像方坐标；P_0～11表示摄影矩阵P的11个元素。

根据公式(2)，将物方点同时投影到模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像像点(x_i,y_i)与待匹配影像像点(x_i',y_i')之间的相关系数。由于在单个匹配模型内，一张基准影像会对应多张待匹配影像，因此一个像素会对应多个相关系数值。将所有相关系数值相加取平均，其结果作为该像素在高程Z_i下的最终相关系数。

在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程。当基准影像所有像素遍历完后，即可生成一张初始的高程图。

然后将基准影像与各待匹配影像分别构成一个立体像对，根据初始高程图，计算基准影像与待匹配影像之间的互信息，如图3所示。

在图3中，首先，根据初始高程图，按照式(4)和式(2)，计算基准影像(设记为I₁)上，每个像素在某待匹配影像(设记为I₂)上的同名点。根据同名点对的灰度，统计一张二维的联合直方图。由二维联合直方图，统计相应的二维灰度概率密度和一维灰度概率密度。

根据灰度概率密度，可以按照式(5)，计算点与点之间的互信息

$h_{I_1,I_2}(i,k)=-1/n\·log\left(P_{I_1,I_2}\right(i,k)\⊗g(i,k\left)\right)\⊗g(i,k)$

$h_I\left(i\right)=-1/n\·log\left(P_I\right(i)\⊗g(i\left)\right)\⊗g\left(i\right)---\left(5\right)$

${\mathrm{mi}}_{I_1,I_2}(i,k)=h_{I_1}\left(i\right)+h_{I_2}\left(k\right)-h_{I_1,I_2}(i,k)$

式中，表示灰度i和灰度k之间的互信息；g(i)表示一维高斯核函数；g(i,k)表示二维高斯核函数；h表示信息熵；表示对应灰度(i,k)的二维灰度概率密度；P_I(i)表示对应灰度i的一维灰度概率密度，h_I(i)代表核线影像I上，灰度i的信息熵。影像I指影像I₁或影像I₂，具体地，代表影像I₁上灰度i的信息熵，代表影像I₂上灰度k的信息熵，P_I1(i)表示影像I₁的灰度概率密度；P_I2(k)表示影像I₂的灰度概率密度；n表示立体像对重叠区域内的像素总数。具体计算过程如图3所示。先将一维/二维直方图采用高斯卷积进行平滑，再采用log运算，最后再对结果进行高斯平滑，即可得到互信息。

在高程范围内，可以根据高程依次计算基准影像上每个像素的“候选同名点”。根据“候选同名点对”，可以按照式(5)计算出对应的互信息，作为匹配代价。每个像素都可以按照高程范围，计算相应的匹配代价。

由于每张基准影像对应多张待匹配影像，因此基准影像上的每个像素，会对应多张待匹配影像上的匹配代价。需要采用一定的策略对多个匹配代价进行融合。本发明根据“像方一致性”和基线约束，实现多个匹配代价的融合。采用二值化算子Census、影像一阶梯度算子Sobel、影像二阶差分算子LOG和基线长度计算权值，当三种测度分别较小或者基线长度较短时，给予较大的权值；反之，则赋较小的权值。将融合后的匹配代价，存入代价矩阵对应的网格中。

Census测度的计算公式为：

C_census(p,d,i)＝Dis_H(R(p),R'(p,d,i))

R(p)＝{ε(p,p+(x,y))}(x,y)∈Np(6)

$\mathrm{\&epsiv;}(p,p^{\&prime;})=\left\{\begin{array}{ccc}0& if& I\left(p\right)<I\left(p^{\&prime;}\right)\\ 1& & otherwise\end{array}\right.$

式中，p表示基准影像上的像素；p'表示p邻域内的某个像素；(x,y)表示p'在p邻域窗口内位置；Np表示像素p的邻域；d表示像素p对应的高程；R(p)表示像素p在基准影像上的census值；R'(p,d,i)表示像素p在高程为d的情况下，在第i张待匹配影像的同名点的census值；Dis_H表示两个census值之间的Hamming距离(汉明距离)；ε(p,p')为二值化运算符；C_census(p,d,i)表示像素p在高程为d的情况下，与第i张待匹配影像上同名点之间的Census测度。

Sobel测度可以表示为：

C_Sobel(p,d,i)＝||G(p)-G'(p,d,i)||

G(p)＝(g_x(p),g_y(p))G'(p,d,i)＝(g_x'(p,d,i),g_y'(p,d,i))(7)

$\begin{array}{cc}{g}_x\left(p\right)=Gx\&CircleTimes;N\left(p\right)& g_y\left(p\right)=Gy\&CircleTimes;N\left(p\right)\end{array}$

$\begin{array}{cc}{g_x}^{\&prime;}(p,d,i)=Gx\&CircleTimes;N^{\&prime;}(p,d,i)& {g_y}^{\&prime;}(p,d,i)=Gy\&CircleTimes;N^{\&prime;}(p,d,i)\end{array}$

式中，C_Sobel(p,d,i)表示像素p在高程为d的情况下，与第i张待匹配影像上同名点之间的Sobel测度；N(p)表示像素p在基准影像上的邻域；N'(p,d,i)表示像素p的同名点在第i张待匹配影像上的邻域；||||表示向量的模；Gx、Gy分别表示Sobel算子模板；G(p)表示像素p在基准影像上的梯度；G(p,d,i)表示像素p的同名点在第i张待匹配影像上的梯度；g_x(p)、g_y(p)分别表示像素p的梯度x方向和y方向的梯度分量；g_x'(p,d,i)、g_y'(p,d,i)分别表示像素p的同名点在第i张待匹配影像上的梯度分量。

同理，LOG测度，也可以通过模板卷积运算得到。LOG模板可以表示为：

LOG测度可以表示为：

C_LOG(p,d,i)＝|L(p)-L'(p,d,i)(8)

$\begin{array}{cc}L\left(p\right)=LOG\&CircleTimes;N\left(p\right)& L^{\&prime;}(p,d,i)=LOG\&CircleTimes;N^{\&prime;}(p,d,i)\end{array}$

式中，C_LOG(p,d,i)表示像素p在高程为d的情况下，与第i张待匹配影像上同名点之间的LOG测度；N(p)表示像素p在基准影像上的邻域；N'(p,d,i)表示像素p的同名点在第i张待匹配影像上的邻域；LOG表示LOG算子模板；L(p)表示像素p在基准影像上的二阶差分；L'(p,d,i)表示像素p的同名点在第i张待匹配影像上的二阶差分。

根据式(6)、式(7)、式(8)以及基线长度，可以计算待匹配影像上匹配代价的权值，并根据权值对多个匹配代价进行加权融合。设第i张待匹配影像记为I_i'，w_i(p,d)表示第i张待匹配影像的权值，权值计算公式为：

w_i(p,d)＝1/(C_census(p,d,i)·C_Sobel(p,d,i)·C_LOG(p,d,i)·baseline(I₁,I_i'))(9)

式中，baseline(I₁,I_i')表示基准影像I₁与待匹配影像I_i'之间的基线长度。实际计算中，需要对多张待匹配影像的权值w_i(p,d)进行归一化操作。

根据式(9)计算的权值，即可对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中：

${\mathrm{mi}}_{I_1}(p,d)=\underset{i=1}{\overset{n}{\&Sigma;}}{w}_i\&CenterDot;{\mathrm{mi}}_{I_1,{I_i}^{\&prime;}}(p,d)---\left(10\right)$

式中，表示p在给定高程d下的最终匹配代价；w_i表示第i张待匹配影像的权值；表示像素p在第i张待匹配影像上的互信息匹配代价。

所述第i张待匹配影像中，i为测区影像集合中的影像序号，i＝1,2,…,N。

最后，半全局密集匹配。根据代价矩阵，可以得到基准影像上每个像素在高程范围内对应的匹配代价。在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，实际上，可以将匹配问题转化为动态规划问题来求解，如图4所示。

动态规划问题的求解方法，如式(11)所示：

$\begin{array}{cc}Lr\left(p,d\right)=C\left(p,d\right)+\min \left\{\begin{array}{c}Lr\left(p-r,d\right),\\ Lr\left(p-r,d-1\right)+P_1,\\ Lr\left(p-r,d+1\right)+P_1,\\ \begin{array}{cc}\underset{i}{\min }& Lr\left(p-r,i\right)+P_2\end{array}\end{array}\right\}-\underset{k}{\min }& Lr\left(p-r,k\right)\end{array}---\left(11\right)$

式中，p表示当前像素；p-r表示前一个阶段的像素；Lr(,)表示动态规划的路径；P₁、P₂表示惩罚系数。r表示动态规划的步距，一般取1；k表示前一个像素的最优视差值；i表示前一个像素的某个视差值。

匹配结果即为动态规划的最优路径。由于一维匹配很容易产生“条纹状”的误匹配，因此采用多个方向的动态规划策略，将匹配路径按照图5所示的方式，划分为8～16个方向，将所有方向的匹配结果进行累加，能够得到更加稳健的匹配结果。图5中，p表示当前的像素；x、y表示像平面坐标的x轴、y轴。

步骤3.全局意义下的多视点云融合：根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合。

单模型的匹配结果不够稳健，易受遮挡、辐射、纹理等因素的影响，需要进一步容纳多视匹配结果，以提高匹配的准确度。本发明通过融合的方法来纳入多视匹配结果，如图6所示。以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网。格网的水平范围即为测区的范围，格网的高程范围即为地表的高程起伏。因为各基准影像的视角不同，为了得到更加精细的融合结果，将所有模型匹配结果分别投影到物方，每个独立的网格都包含零至多个物方点，当测区影像集合中的影像数目为n，物方点即为图6中的匹配高程图1～n相应点云1～n。对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上Census值的方差作为像方可见性约束。一般来说，网格内点数越大、方差越小，则代表对应物方点为最优高程的可能性越大。

根据物方一致性约束和像方可见性约束，计算每个网格点的高程置信度，如下式所示：

conf＝-N/σ_Census

${\mathrm{\&sigma;}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\&OverBar;}{Census})}^2/n}---\left(12\right)$

式中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能够可见该网格点的影像数目，具体数目信息可在物方投影的同时得到；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值。置信度越小，表示该网格点对应的高程为最优高程的可能性越大。

仅仅根据局部孤立网格点的高程置信度来确定最优高程，很容易陷入局部极值解。考虑场景高程是分片连续的特点，网格点的高程之间，存在一定的平滑约束。在确定格网点高程的同时，考虑周围格网点高程的情况，理论上能够得到更加稳健的融合结果。根据高程置信度和高程平滑约束，可以构建如下所示的全局能量函数：

$E=min\underset{i\&Element;G}{\&Sigma;}(Cost\left(H_i\right)+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}\mathrm{\&sigma;}|H_i-H_q|/|I_i-I_q\left|\right)---\left(13\right)$

其中，E表示能量函数，作为融合结果的测度，E的极小值所对应的融合结果，即为全局最优的融合结果；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，相应地H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值。该能量函数充分考虑场景的高程是分片连续的情况，当第i个格网点与其邻域点q的影像灰度相似时，可以认为是属于同一表面上的点，其高程应该是连续平滑的，此时需要H_i和H_q很接近，才能得到较小的能量函数E；反之，当I_i和I_q相差很大时，即使H_i和H_q相差较大，也不会给能量函数带来太大的影响，因此能够保证边缘处高程的阶跃性。根据实验所得基于全局能量函数最小化的多视点云融合结果，可知全局能量函数考虑了格网点之间的高程平滑。

公式(13)是定义域在二维空间的多元函数，其最优解是一个NP难题，可以采用图割算法或者置信传播算法来计算，但无论是图割算法还是置信传播算法，其算法复杂度均很高，运行时间长，效率低下，不适合地理国情监测、智慧城市等实际应用需求。本发明仍然采用多方向动态规划的方法求解上式，即可得到全局意义下的最优格网点。最后采用加权融合的方法，将网格内的点融合为一个点。

网格内，点云融合的公式为：

$\begin{array}{ccc}{X}^{\&prime;}=\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{\&rho;}}_i{X}_i& Y^{\&prime;}=\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{\&rho;}}_i{Y}_i& Z^{\&prime;}=\underset{i=1}{\overset{N}{\&Sigma;}}{\mathrm{\&rho;}}_i{Z}_i\end{array}$

${\mathrm{\&rho;}}_i=P_i/\underset{j=1}{\overset{N}{\&Sigma;}}{P}_j---\left(14\right)$

$P_i=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{2.5}\mathrm{\&sigma;}}\exp \left(-\frac{1}{2{\mathrm{\&sigma;}}^2}\left({\left(X_i-\stackrel{\&OverBar;}{X}\right)}^2+{\left(Y_i-\stackrel{\&OverBar;}{Y}\right)}^2+{\left(Z_i-\stackrel{\&OverBar;}{Z}\right)}^2\right)\right)$

上式采用正态分布函数计算融合的权值。式中，X'、Y'、Z'表示融合后的点坐标；N表示网格内点的总数；X_i、Y_i、Z_i表示网格内第i个点；σ表示标准差；表示网格内物方点坐标的平均值；exp表示以自然常数e为底的指数函数；P_i或P_j表示网格内第i或j个物方点的高斯权值；ρ_i表示网格内第i个物方点经过归一化后的权值。

步骤4.联合面特征和线特征的点云优化：联合面特征和线特征，对点云进行优化，修正纹理贫乏和重复区域的误匹配，精化线状地物的边缘，生成高质量的最终点云。

(1)基于“面特征”的点云优化方法

受遮挡、纹理贫乏等因素的影响，融合后的点云在物方仍然可能存在“无”高程的区域，即“空洞”。在实际作业中“空洞”往往对应湖泊、雪地、房屋楼底等重要地物，是重建DSM不可缺少的元素。首先根据“空洞”的平面位置找出覆盖该“空洞”的基准影像集合以及对应的匹配结果，如图7所示；对基准影像进行图像分割，提取面特征；将融合后的点云反投影到原始基准影像上；统计每个分割块内的有效匹配点(即高程信息)，如图8所示；根据有效高程信息，选择一个最优的高程平面参数，如式(15)所示，模拟分割块内的高程变化；最后，根据最优高程平面参数，计算分割块内“空洞”的高程，并投影到物方。

图7中，测区包括山地、湖泊，虚线代表物方有效高程，叉号代表“空洞”。根据空洞周围的有效视差，即可迅速确定覆盖该“空洞”的影像，并自动找到其对应的立体匹配像对集合，如图7中圆形虚线所示。

参见图8中影像分割结果，在每个分割块内，用圆点表示每个像素所对应的深度，不同灰度代表不同的深度。叉号表示无效深度，即“空洞”。在图8中，需要计算包含“空洞”的分割块内的深度平面方程，如式(15)所示。在同一个分割块内，经常存在不同的有效深度信息(如图8中的包含“空洞”的分割块内，存在两种深度信息)，为了找出最优深度平面方程，要求利用的有效深度信息满足如式(16)所示的条件。

d_i＝ax_i+by_i+c(15)

其中，i表示分割块内第i个像素；d_i表示像素(x_i,y_i)对应的高程；a,b,c表示高程平面方程的参数。

$E=min\underset{i\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}\left(Cost\left(x_i,y_i\right)+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}{P}_{1}T\&lsqb;|D_i-D_q|=1\&rsqb;+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}{P}_{2}T\&lsqb;|D_i-D_q|>1\&rsqb;\right)---\left(16\right)$

其中，Ω表示分割块内所有像素的集合；Cost(x_i,y_i)表示集合Ω中第i个像素(x_i,y_i)对应的匹配代价；P₁、P₂表示惩罚系数，一般情况下P₂不能小于P₁；函数T[]表示当[]内条件为真时，T＝1，当[]内条件为假时，T＝0；D_i表示以高程步距为单位，第i个格网点的高程；D_q表示以高程步距为单位，第q个格网点的高程，q在i的邻域范围内。

(2)基于“线特征”的点云优化方法

三维重建要求房屋点云的边缘要尽量精确。但是现有的密集匹配算法，无论是局部方法还是全局方法，房屋边缘的匹配结果，会出现严重的“毛刺”现象。为了解决这个问题，本发明设计了一种线特征约束的点云优化方法，充分利用影像中的边缘可以被精确提取的特点，利用基准影像中提取的边缘来改正物方点云的边缘。所述线特征约束的点云优化方法实现方式为，首先在基准影像上提取直线边缘，特别是房屋屋顶的边缘。根据基准影像-高程图一一对应的关系，在高程图上找到该边缘的位置，以该边缘为中心，建立一个缓冲区；统计缓冲区边缘两侧的高程分布情况，根据少数服从多数的原则，修正边缘两侧的错误高度，剔除毛刺。

从实验可以看出，原始影像上能高精度地提取直线，但是对应的深度图毛刺现象严重，需要进一步优化。因此可以以精确提取的直线为中心，建立缓冲区，统计缓冲区域内，边缘两侧的深度分布情况，根据少数服从多数原则，修正边缘两侧的“毛刺”，达到精化边缘的目的。

具体实施时，还可以采用模块化方式提供系统。本发明提供一种基于多特征多约束的多视影像密集匹配融合系统，包括以下模块：

匹配模型构建模块，用于根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束；

高程图提取模块，用于对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图，包括如下子模块，

高程范围预测子模块，用于根据物方点的高程，预测基准影像所覆盖测区的地表起伏，找到最大高程Z_max和最小高程Z_min；

步距计算子模块，用于从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面和最小高程面，分别有两个交点M、M'，在步骤1所得待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'，将最大最小高程差与投影点连线长度之间的比值，作为高程步距；

代价矩阵构造子模块，用于以基准影像的像平面作为水平面，根据高程范围预测子模块所得高程范围，根据步距计算子模块所得高程步距在高程方向划分间隔，建立一个三维网格矩阵，作为代价矩阵；所述代价矩阵中，每一个网格代表某像素在对应高程Z_i下的匹配代价；

匹配代价提取子模块，用于在代价矩阵中，以互信息作为匹配代价，实现方式如下，

首先，计算基准影像的初始高程图，包括针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i，计算对应的物方点坐标(X_i,Y_i,Z_i)；将物方点同时投影到匹配模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像上像点(x_i,y_i)与各待匹配影像上dioxide像点(x_i',y_i')之间的相关系数并取平均，作为该像素在高程Z_i下的最终相关系数，在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程；

然后，将基准影像与各待匹配影像构成一个立体像对，根据初始高程图，对基准影像上的每个像素计算与各待匹配影像之间的互信息作为匹配代价，对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中；

最后，进行半全局密集匹配，包括根据代价矩阵得到基准影像上每个像素在高程范围内对应的匹配代价，在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，将匹配问题转化为动态规划问题求解，匹配结果为动态规划的最优路径，得到基准影像的高程图；

融合模块，用于根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合，实现方式如下，

以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网，三维格网的水平范围为测区的范围，三维格网的高程范围为地表的高程起伏，分别将不同视角的模型匹配结果投影到物方，每个独立的网格都包含零至多个物方点，物方点为步骤3针对各匹配模型所得高程图的相应点云，对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上二值化算子Census值的方差作为像方可见性约束，计算每个网格点的高程置信度如下，

conf＝-N/σ_Census

${\mathrm{\&sigma;}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\&OverBar;}{Census})}^2/n}$

其中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能见该网格点的影像数目；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值；

所述全局能量函数如下，

$E=min\underset{i\&Element;G}{\&Sigma;}(Cost\left(H_i\right)+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}\mathrm{\&sigma;}|H_i-H_q|/|I_i-I_q\left|\right)$

其中，E表示能量函数，作为融合结果的测度；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值；

点云优化模块，用于联合面特征和线特征进行点云优化，生成最终的点云。

各模块具体实现可参见相应步骤，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于多特征多约束的多视影像密集匹配融合方法，其特征在于，包括以下步骤：

步骤1，根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束；

步骤2，对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图，实现方式如下，

1)预测高程范围，包括根据特征匹配生成的物方点的高程，预测基准影像所覆盖测区的地表起伏，找到最大高程Z_max和最小高程Z_min；

2)计算高程方向的步距，包括从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面和最小高程面，分别有两个交点M、M'，在步骤1所得待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'，将最大最小高程差与投影点连线长度之间的比值，作为高程步距；

3)构造代价矩阵，包括以基准影像的像平面作为水平面，根据1)所得高程范围，根据2|)所得高程步距在高程方向划分间隔，建立一个三维网格矩阵，作为代价矩阵；所述代价矩阵中，每一个网格代表某像素在对应高程Z_i下的匹配代价；

4)在代价矩阵中，以互信息作为匹配代价，实现方式如下，

首先，计算基准影像的初始高程图，包括针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i，计算对应的物方点坐标(X_i,Y_i,Z_i)；将物方点同时投影到匹配模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像上像点(x_i,y_i)与各待匹配影像上像点(x_i',y_i')之间的相关系数并取平均，作为该像素在高程Z_i下的最终相关系数，在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程；

然后，将基准影像分别与各待匹配影像构成一个立体像对，根据初始高程图，对基准影像上的每个像素计算与各待匹配影像之间的互信息作为匹配代价，对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中；

最后，进行半全局密集匹配，包括根据代价矩阵得到基准影像上每个像素在高程范围内对应的匹配代价，在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，将匹配问题转化为动态规划问题求解，匹配结果为动态规划的最优路径，得到基准影像的高程图；

步骤3，根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合，实现方式如下，

以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网，三维格网的水平范围为测区的范围，三维格网的高程范围为地表的高程起伏，分别将不同视角的模型匹配结果投影到物方，每个独立的网格都包含零至多个物方点，物方点为步骤3针对各匹配模型所得高程图的相应点云，对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上二值化算子Census值的方差作为像方可见性约束，计算每个网格点的高程置信度如下，

conf＝-N/σ_Census

${\mathrm{\&sigma;}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\&OverBar;}{Census})}^2/n}$

其中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能见该网格点的影像数目；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值；

所述全局能量函数如下，

$E=min\underset{i\&Element;G}{\&Sigma;}(Cost\left(H_i\right)+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}\mathrm{\&sigma;}|H_i-H_q|/|I_i-I_q\left|\right)$

其中，E表示能量函数，作为融合结果的测度；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值；

步骤4，联合面特征和线特征进行点云优化，生成最终的点云。

2.根据权利要求1所述的基于多特征多约束的多视影像密集匹配融合方法，其特征在于：步骤4中，基于面特征进行优化的实现方式如下，

对融合后的点云在物方存在无高程的区域，记为“空洞”，首先根据“空洞”的平面位置找出覆盖该“空洞”的基准影像集合以及对应的匹配结果；对基准影像进行图像分割，提取面特征；将融合后的点云反投影到原始的基准影像上；统计每个分割块内的有效匹配点得到有效高程信息；根据有效高程信息，选择一个最优的高程平面参数，模拟分割块内的高程变化；最后，根据最优高程平面参数，计算分割块内“空洞”的高程，并投影到物方。

3.根据权利要求1所述的基于多特征多约束的多视影像密集匹配融合方法，其特征在于：步骤4中，基于线特征进行优化的实现方式如下，

首先在基准影像上提取直线边缘，根据基准影像与高程图一一对应的关系，在高程图上找到该边缘的位置，以该边缘为中心，建立一个缓冲区；统计缓冲区边缘两侧的高程分布情况，根据少数服从多数的原则，修正边缘两侧的错误高度，剔除毛刺。

4.一种基于多特征多约束的多视影像密集匹配融合系统，其特征在于，包括以下模块：

匹配模型构建模块，用于根据多重约束，为每一张基准影像分别选择若干张待匹配影像，得到待匹配影像集合，基准影像和相应匹配影像集合构成一个匹配模型；所述多重约束，包括基线约束、像平面法向量约束、特征匹配约束和交会角约束；

高程图提取模块，用于对各个匹配模型，利用多视约束条件，进行半全局密集匹配，直接生成单个匹配模型的密集匹配结果，得到相应高程图，包括如下子模块，

高程范围预测子模块，用于根据特征匹配产生的物方点的高程，预测基准影像所覆盖测区的地表起伏，找到最大高程Z_max和最小高程Z_min；

步距计算子模块，用于从基准影像的摄影中心出发，经过像主点，引出一条射线。射线与最大高程面和最小高程面，分别有两个交点M、M'，在步骤1所得待匹配影像集合中，选择一张与基准影像距离最远的一张待匹配影像I'，将交点M、M'投影到影像I'上，分别得到对应的像点m、m'，将最大最小高程差与投影点连线长度之间的比值，作为高程步距；

代价矩阵构造子模块，用于以基准影像的像平面作为水平面，根据高程范围预测子模块所得高程范围，根据步距计算子模块所得高程步距在高程方向划分间隔，建立一个三维网格矩阵，作为代价矩阵；所述代价矩阵中，每一个网格代表某像素在对应高程Z_i下的匹配代价；

匹配代价提取子模块，用于在代价矩阵中，以互信息作为匹配代价，实现方式如下，

首先，计算基准影像的初始高程图，包括针对基准影像上每一个像素(x_i,y_i)，根据高程范围内的高程值Z_i，计算对应的物方点坐标(X_i,Y_i,Z_i)；将物方点同时投影到匹配模型内所有的待匹配影像上，得到对应的像点坐标(x_i',y_i')；计算基准影像上像点(x_i,y_i)与各待匹配影像上dioxide像点(x_i',y_i')之间的相关系数并取平均，作为该像素在高程Z_i下的最终相关系数，在高程范围Z_min～Z_max内，找到相关系数最大的同名点，对应的高程为像素(x_i,y_i)的初始高程；

然后，将基准影像分别与各待匹配影像构成一个立体像对，根据初始高程图，对基准影像上的每个像素计算与各待匹配影像之间的互信息作为匹配代价，对多张待匹配影像的匹配代价进行融合，得到最终的匹配代价，并存入相应的代价矩阵中；

最后，进行半全局密集匹配，包括根据代价矩阵得到基准影像上每个像素在高程范围内对应的匹配代价，在基准影像上，将任意方向直线上的每个像素作为一个阶段，像素的匹配代价作为节点，将匹配问题转化为动态规划问题求解，匹配结果为动态规划的最优路径，得到基准影像的高程图；

融合模块，用于根据格网点之间的高程平滑约束，在全局能量函数最小的条件下，对多个匹配模型的密集匹配结果进行融合，实现方式如下，

以物方坐标X、Y、Z方向为轴，建立一个包含整个测区的三维格网，三维格网的水平范围为测区的范围，三维格网的高程范围为地表的高程起伏，分别将不同视角的模型匹配结果投影到物方，每个独立的网格都包含零至多个物方点，物方点为步骤3针对各匹配模型所得高程图的相应点云，对落在同一个网格中的不同视角的物方点进行统计，以网格内物方点的数量作为物方一致性约束，以网格中心在多张影像上二值化算子Census值的方差作为像方可见性约束，计算每个网格点的高程置信度如下，

conf＝-N/σ_Census

${\mathrm{\&sigma;}}_{Census}=\sqrt{\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{Dis}}_H{({\mathrm{Census}}_i,\stackrel{\&OverBar;}{Census})}^2/n}$

其中，conf表示高程置信度；N表示网格内的物方点数目；σ_Census表示网格中心在多张影像上Census测度的方差；n表示能见该网格点的影像数目；Dis_H表示Hamming距离；Census_i表示物方网格点投影到第i张影像上所对应像点的Census值；表示所有影像上投影像点的Census均值；

所述全局能量函数如下，

$E=min\underset{i\&Element;G}{\&Sigma;}(Cost\left(H_i\right)+\underset{q\&Element;N_i}{\mathrm{\&Sigma;}}\mathrm{\&sigma;}|H_i-H_q|/|I_i-I_q\left|\right)$

其中，E表示能量函数，作为融合结果的测度；G表示物方所有格网点的集合；H_i代表第i个格网点的高程，Cost(H_i)表示第i个格网点的对应的高程置信度；N_i表示第i个格网点的邻域格网；σ表示平滑项系数，用于控制格网点之间的平滑强度；I_i表示第i个格网点对应多张影像像点灰度的均值，H_q代表邻域点q的高程，I_q表示邻域点q对应多张影像像点灰度的均值；

点云优化模块，用于联合面特征和线特征进行点云优化，生成最终的点云。

5.根据权利要求4所述的基于多特征多约束的多视影像密集匹配融合系统，其特征在于：点云优化模块中，基于面特征进行优化的实现方式如下，

对融合后的点云在物方存在无高程的区域，记为“空洞”，首先根据“空洞”的平面位置找出覆盖该“空洞”的基准影像集合以及对应的匹配结果；对基准影像进行图像分割，提取面特征；将融合后的点云反投影到原始的基准影像上；统计每个分割块内的有效匹配点得到有效高程信息；根据有效高程信息，选择一个最优的高程平面参数，模拟分割块内的高程变化；最后，根据最优高程平面参数，计算分割块内“空洞”的高程，并投影到物方。

6.根据权利要求4所述的基于多特征多约束的多视影像密集匹配融合系统，其特征在于：点云优化模块中，基于线特征进行优化的实现方式如下，

首先在基准影像上提取直线边缘，根据基准影像与高程图一一对应的关系，在高程图上找到该边缘的位置，以该边缘为中心，建立一个缓冲区；统计缓冲区边缘两侧的高程分布情况，根据少数服从多数的原则，修正边缘两侧的错误高度，剔除毛刺。