CN106228593A - 一种图像密集匹配方法 - Google Patents

Abstract

本发明公开了一种图像密集匹配方法，包括提取出特征点图像，构建旋转平移不变特征向量；进行同名初匹配，得到具有稳定相对关系的最大同名点对子集；构建两幅相邻图像平面的第一Delaunay三角网，核实各对同名三角形在两幅相邻图像中的拓扑关系是否准确；移除所有拓扑关系不准确的同名点，构建第二Delaunay三角网，并计算仿射变换参数；计算密集匹配在第一图像的各非特征像点位置，判断各非特征像点与第二Delaunay三角网中三角形的拓扑包容关系，并计算多个仿射变换参数在第二图像中的多个初始同名像点位置；建立两幅相邻图像的局部相关矩阵；计算得到密集匹配结果；对立体观测区域进行三维形貌重建。

Description

一种图像密集匹配方法

技术领域

本发明涉及图像三维重建技术，特别地，涉及一种图像密集匹配方法。

背景技术

图像密集匹配是摄影测量或计算机视觉的图像三维形貌重建工作中最重要的工作、也是当前数字城市建设中的核心工作之一。仅靠图像内稀疏不规则的特征点，无法让普通人群理解图像所观测目标场景的三维形貌，必须依靠图像内占绝大部分比例的非特征点位所构成的密集匹配，才能准确全面描述它。而当前各种方法中，对于纹理贫乏、相邻局部区域相似度较高的三维区域，始终存在对图像内非特征点位难以匹配的问题，且匹配时运算量大，匹配精度低。

针对现有技术中图像密集匹配方法在纹理贫乏、相邻局部区域相似度较高的三维区域中匹配运算量大、匹配精度低的问题，目前尚未有有效的解决方案。

发明内容

针对现有技术中图像密集匹配方法在纹理贫乏、相邻局部区域相似度较高的三维区域中匹配运算量大、匹配精度低的问题，本发明的目的在于提出一种图像密集匹配方法，能够在纹理贫乏、相邻局部区域相似度较高的三维区域中进行工作，使用更少的运算量得到更高精度水平的密集匹配结果。

基于上述目的，本发明提供的技术方案如下：

根据本发明的一个方面，提供了一种图像密集匹配方法，包括：

获取立体观测区域的框幅式序列图像；

从每幅序列图像中选择两幅相邻的图像，提取出两幅相邻图像各自的特征点图像，并根据特征点图像分别构建旋转平移不变特征向量，两幅相邻图像包括第一图像与第二图像；

在核线约束下，使用随机抽样一致算法对旋转平移不变特征向量进行同名初匹配，得到具有稳定相对关系的最大同名点对子集；

根据最大同名点对子集构建两幅相邻图像平面的第一Delaunay三角网，计算第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实各对同名三角形在两幅相邻图像中的拓扑关系是否准确；

根据各对同名三角形在两幅相邻图像中的拓扑关系是否准确的核实结果，移除所有拓扑关系不准确的同名点，构建第二Delaunay三角网，并计算第二Delaunay三角网中各对同名三角形之间的仿射变换参数；

获取密集匹配所需的平面像素间距，根据平面像素间距计算密集匹配在第一图像的各非特征像点位置，并根据各非特征像点与第二Delaunay三角网判断各非特征像点与第二Delaunay三角网中三角形的拓扑包容关系，根据拓扑包容关系获得各像点位置的初始点应调用的同名三角形的多个仿射变换参数，并计算多个仿射变换参数在第二图像中的多个初始同名像点位置；

以第一图像当前的各非特征像点位置和第二图像多个初始同名像点位置为中心，选取多个指定半径的同心圆与辐射状点阵，建立两幅相邻图像的局部相关矩阵；

对局部相关矩阵使用辐射几何改正的最小二乘平差算法，迭代计算得到具有最大相关系数的亚像素像点位置的集合，像点位置的集合为两幅相邻图像的密集匹配结果；

使用上述方法计算出框幅式序列图像中任意两幅相邻图像的密集匹配结果，确定密集匹配结果中空间参考平面与立体图像的内外参数，并根据密集匹配结果、空间参考平面与立体图像的内外参数对立体观测区域进行三维形貌重建。

其中，提取出两幅相邻图像各自的特征点图像包括：

将两幅相邻图像的每一幅划分为多个大小相等的像元；

在每一个像元中的相同位置任取一点，计算该点在该点所在的像元中上下左右四个方向的灰度差分绝对值；

获取事先确定的阈值，并根据点在点所在的像元中上下左右四个方向的灰度差分绝对值确定协方差矩阵与误差椭圆，并根据协方差矩阵与误差椭圆确定每个像元中一初选点与一备选点；

根据每个像元中的一初选点与一备选点，确定每个像元中的特征点；

将所有像元中的特征点组合生成两幅相邻图像的两幅特征点图像。

其中，计算第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实各对同名三角形在两图像中的拓扑关系是否准确包括：

根据第一Delaunay三角网中某一对同名三角形Δp₁p₂p₃与Δp′₁p′₂p′₃的三顶点坐标 $\left\{\begin{array}{c}{p}_1(x_1,y_1),p_2(x_2,y_2),p_3(x_3,y_3)\\ p_1^{\′}(x_1^{\′},y_1^{\′}),p_2^{\′}(x_2^{\′},y_2^{\′}),p_3^{\′}(x_3^{\′},y_3^{\′})\end{array}\right.,$ 计算获得顶点坐标构成的平面矢量

根据顶点坐标构成的平面矢量,计算获得该对同名三角形Δp₁p₂p₃与Δp′₁p′₂p′₃的矢量叉积面积与

若判定拓扑同名关系准确、

若判定拓扑同名关系不准确；

依次对第一Delaunay三角网中的每一对同名三角形进行上述计算，并获得判定结果。

并且，据各非特征像点与第二Delaunay三角网判断各非特征像点与第二Delaunay三角网中三角形的拓扑包容关系包括：

过非特征像点p₀(x₀,y₀)作射线

对于第二Delaunay三角网中的某一个三角形Δp₁p₂p₃，依次判断有向线段与射线 $\stackrel{\→}{l}:y=y_0,x\≥x_0$ 是否存在交点；

统计射线与有向线段的交点个数总数并进行判断：

若交点个数总数为奇数，则非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃内、

若交点个数总数为偶数，则非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃外。

并且，依次判断有向线段与射线是否存在交点包括；

对于

首先，当y₀≤y₁且y₀≤y₂时，或者当y₀≥y₁且y₀≥y₂时，判定有向线段与射线 $\stackrel{\→}{l}:y=y_0,x\≥x_0$ 不存在交点；

其次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≤y₂时，计算Δp₁p₂p₀的面积其中

若判定射线与有向线段存在交点、

若表示射线与有向线段不存在交点；

再次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≥y₂时，计算Δp₁p₂p₀的面积其中

若判定射线与有向线段存在交点、

若表示射线与有向线段不存在交点；

对于与与射线是否存在交点的判定方式与相同。

并且，对局部相关矩阵使用辐射几何改正的最小二乘平差算法包括：

计算两幅相邻图像各自的灰度值、均值和标准差，并根据两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正；

计算两幅相邻图像各自的形变函数与随机噪声，并根据两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正，进行最小二乘平差计算。

一方面，计算两幅相邻图像各自的灰度值、均值和标准差，并根据两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正包括：

获取第一图像的灰度值、灰度均值与标准差以及第二图像的灰度值、灰度均值与标准差

根据第一图像的灰度值、灰度均值与标准差以及第二图像的灰度值、灰度均值与标准差使用中心化方法对第一图像与第二图像的灰度值g_l,g_r进行标准化处理，得到 $\left\{\begin{array}{c}{g}_l^{\′}=\frac{g_l-\stackrel{\‾}{g_l}}{{\mathrm{\σ}}_{g_l}}\\ g_r^{\′}=\frac{g_r-\stackrel{\‾}{g_r}}{{\mathrm{\σ}}_{g_r}}\end{array}\right.,$ 其中，g′_l,g′_r分别为标准化处理后的两幅图像的灰度值；

将第一图像的灰度值g_l与第二图像的灰度值g_r使用线性变换进行匹配，有g′_l＝h₁+h₂g′_r，可获得变换系数h₁,h₂，变换系数h₁,h₂可对辐射畸变进行匹配纠正。

另一方面，计算两幅相邻图像各自的形变函数与随机噪声，并根据两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正包括：

获取第一图像的形变函数g_l(x,y)与随机噪声函数n_l，以及第二图像的形变函数g_r(x,y)与随机噪声函数n_r；

根据第一图像与第二图像在匹配后的形变与噪声一致性，可得g_l(x,y)+n_l＝g_r(x,y)+n_r，并计算出误差方程v＝n_r-n_l＝g_l(x,y)-g_r(x,y)；

根据误差方程与几何畸变误差最小原则进行几何畸变的匹配纠正。

从上面所述可以看出，本发明提供的技术方案通过提取出图像特征点并使用旋转平移不变特征向量进行初匹配并构建Delaunay三角网获取仿射变换参数以建立相关矩阵进而得到密集匹配结果的技术方案，能够在纹理贫乏、相邻局部区域相似度较高的三维区域中进行工作，使用更少的运算量得到更高精度水平的密集匹配结果。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为根据本发明实施例的一种图像密集匹配方法的流程图；

图2为根据本发明实施例的一种图像密集匹配方法中第一图像上同名三角形Δp₁p₂p₃的位置图；

图3为根据本发明实施例的一种图像密集匹配方法中第二图像上同名三角形Δp′₁p′₂p′₃的位置图；

图4为根据本发明实施例的一种图像密集匹配方法中判定非特征像点p₀(x₀,y₀)是否位于三角形Δp₁p₂p₃中的关系示意图；

图5为根据本发明实施例的一种图像密集匹配方法中，当y₀≤y₁且y₀≤y₂时，有向线段与射线不存在交点的关系示意图；

图6为根据本发明实施例的一种图像密集匹配方法中，当y₀≥y₁且y₀≥y₂时，有向线段与射线不存在交点的关系示意图；

图7为根据本发明实施例的一种图像密集匹配方法中，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≤y₂、且时，射线与有向线段存在交点关系示意图；

图8为根据本发明实施例的一种图像密集匹配方法中，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≤y₂、且时，射线与有向线段不存在交点的关系示意图；

图9为根据本发明实施例的一种图像密集匹配方法中，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≥y₂、且时，射线与有向线段存在交点的关系示意图；

图10为根据本发明实施例的一种图像密集匹配方法中，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≥y₂、且时，射线与有向线段不存在交点的关系示意图；

图11为根据本发明实施例的一种图像密集匹配方法中选取多个指定半径的同心圆与辐射状点阵的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进一步进行清楚、完整、详细地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员所获得的所有其他实施例，都属于本发明保护的范围。

根据本发明的实施例，提供了一种图像密集匹配方法。

如图1所示，根据本发明实施例提供的立体图像快速三维重建方法包括：

步骤S101，获取立体观测区域的框幅式序列图像；

步骤S103，从每幅序列图像中选择两幅相邻的图像，提取出两幅相邻图像各自的特征点图像，并根据特征点图像分别构建旋转平移不变特征向量，两幅相邻图像包括第一图像与第二图像；

步骤S105，在核线约束下，使用随机抽样一致算法对旋转平移不变特征向量进行同名初匹配，得到具有稳定相对关系的最大同名点对子集；

步骤S107，根据最大同名点对子集构建两幅相邻图像平面的第一Delaunay三角网，计算第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实各对同名三角形在两幅相邻图像中的拓扑关系是否准确；

步骤S109，根据各对同名三角形在两幅相邻图像中的拓扑关系是否准确的核实结果，移除所有拓扑关系不准确的同名点，构建第二Delaunay三角网，并计算第二Delaunay三角网中各对同名三角形之间的仿射变换参数；

步骤S111，获取密集匹配所需的平面像素间距，根据平面像素间距计算密集匹配在第一图像的各非特征像点位置，并根据各非特征像点与第二Delaunay三角网判断各非特征像点与第二Delaunay三角网中三角形的拓扑包容关系，根据拓扑包容关系获得各像点位置的初始点应调用的同名三角形的多个仿射变换参数，并计算多个仿射变换参数在第二图像中的多个初始同名像点位置；

步骤S113，以第一图像当前的各非特征像点位置和第二图像多个初始同名像点位置为中心，选取多个指定半径的同心圆与辐射状点阵，建立两幅相邻图像的局部相关矩阵；

步骤S115，对局部相关矩阵使用辐射几何改正的最小二乘平差算法，迭代计算得到具有最大相关系数的亚像素像点位置的集合，像点位置的集合为两幅相邻图像的密集匹配结果；

步骤S117，使用上述方法计算出框幅式序列图像中任意两幅相邻图像的密集匹配结果，确定密集匹配结果中空间参考平面与立体图像的内外参数，并根据密集匹配结果、空间参考平面与立体图像的内外参数对立体观测区域进行三维形貌重建。

对绝大部分情况下的场景目标，框幅式数字相机基于空中或地面移动各种运动平台连续成像观测时，不论是从地面或空中飞行等方式获取其移动序列图像数据，也不需要事先对待观测区域布设测量标志，只要图像间存在较大幅面的重叠度，并且成像质量达到正常状态，则即使所观测区域的纹理相当贫乏，以常规方式难以匹配重建(如沙漠、戈壁等)，都可通过本发明的技术方案进行处理，得到高自动化水平、亚像素级的密集匹配结果，为三维形貌重建提供可靠基础数据。

其中，提取出两幅相邻图像各自的特征点图像包括：

将两幅相邻图像的每一幅划分为多个大小相等的像元；

在每一个像元中的相同位置任取一点，计算该点在该点所在的像元中上下左右四个方向的灰度差分绝对值；

获取事先确定的阈值，并根据点在点所在的像元中上下左右四个方向的灰度差分绝对值确定协方差矩阵与误差椭圆，并根据协方差矩阵与误差椭圆确定每个像元中一初选点与一备选点；

根据每个像元中的一初选点与一备选点，确定每个像元中的特征点；

将所有像元中的特征点组合生成两幅相邻图像的两幅特征点图像。

这种提取特征点图像的方法相比于现有技术中的灰度方差算法或自相关函数曲率算法，计算速度更快，计算精度可达到亚像素级，且具有良好的抗噪性。

根据特征点图像分别构建的旋转平移不变特征向量能较好的表现出特征点图像所包含的图像局部特征。图像局部特征是视觉领域中特异性的体现，通常用来描述一块区域，使之与其他区域相比具有高区分度。局部特征建立非常重要，局部特征建立的好坏直接会决定着后面分类、识别的精度。旋转平移不变特征向量是一种提取局部特征的算法，旋转平移不变特征向量在尺度空间寻找极值点，提取亚像素级的位置、尺度、旋转、亮度等数据的不变量成分，对视角变化、仿射变换、随机噪声也保持一定程度的稳定性。旋转平移不变特征向量的独特性好，信息量丰富，适用于在海量特征数据库中进行快速、准确的匹配，能够提高匹配速度并降低匹配计算量。

其中，计算第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实各对同名三角形在两图像中的拓扑关系是否准确包括：

图2与图3分别示出的是同名三角形Δp₁p₂p₃与Δp′₁p′₂p′₃在第一图像与第二图像的特征点图像上的位置。如图2与图3所示，根据第一Delaunay三角网中某一对同名三角形Δp₁p₂p₃与Δp′₁p′₂p′₃的三顶点坐标 $\left\{\begin{array}{c}{p}_1(x_1,y_1),p_2(x_2,y_2),p_3(x_3,y_3)\\ p_1^{\′}(x_1^{\′},y_1^{\′}),p_2^{\′}(x_2^{\′},y_2^{\′}),p_3^{\′}(x_3^{\′},y_3^{\′})\end{array}\right.,$ 计算获得顶点坐标构成的平面矢量

根据顶点坐标构成的平面矢量,计算获得该对同名三角形Δp₁p₂p₃与Δp′₁p′₂p′₃的矢量叉积面积与

若判定拓扑同名关系准确、

若判定拓扑同名关系不准确；

依次对第一Delaunay三角网中的每一对同名三角形进行上述计算，并获得判定结果。

核实各对同名三角形在两图像中的拓扑关系是否准确可以排除在步骤S105中进行同名初匹配时发生误匹配的同名点，并基于所有准确匹配的同名点构建第二Delaunay三角网进行下一步计算，避免误匹配的同名点在后续的计算中造成偏差。

并且，据各非特征像点与第二Delaunay三角网判断各非特征像点与第二Delaunay三角网中三角形的拓扑包容关系包括：

过非特征像点p₀(x₀,y₀)作射线

对于第二Delaunay三角网中的某一个三角形Δp₁p₂p₃，依次判断有向线段与射线 $\stackrel{\→}{l}:y=y_0,x\≥x_0$ 是否存在交点；

统计射线与有向线段的交点个数总数并进行判断：

若交点个数总数为奇数，则非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃内、

若交点个数总数为偶数，则非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃外。

图4示出的是判定非特征像点p₀(x₀,y₀)是否位于三角形Δp₁p₂p₃中的关系示意图。在平面p₁p₂p₃内，判断点p₀(x₀,y₀)是否被三角形Δp₁p₂p₃拓扑包容，即判断计算几何中的数学关系。从点p₀(x₀,y₀)引一条射线计算射线与三角形Δp₁p₂p₃各条边所在线段的交点，统计有向线段存在的交点数量，当交点数量为奇数个，表示该点落在三角形Δp₁p₂p₃内；当交点数量为偶数个，表示该点落在三角形Δp₁p₂p₃外。

并且，依次判断有向线段与射线是否存在交点包括；

对于

首先，当y₀≤y₁且y₀≤y₂时(如图5所示)，或者当y₀≥y₁且y₀≥y₂时(如图6所示)，判定有向线段与射线不存在交点；

其次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≤y₂时，计算Δp₁p₂p₀的面积其中

若(如图7所示)，判定射线与有向线段存在交点、

若(如图8所示)，表示射线与有向线段不存在交点；

再次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≥y₂时，计算Δp₁p₂p₀的面积其中

若(如图9所示)，判定射线与有向线段存在交点、

若(如图10所示)，表示射线与有向线段不存在交点；

对于与与射线是否存在交点的判定方式与相同。

步骤S113中选取多个指定半径的同心圆与辐射状点阵，指定半径的同心圆是指以第一图像当前的各非特征像点位置或第二图像多个初始同名像点位置为中心的，多个不同半径的同心圆；辐射状点阵是指过第一图像当前的各非特征像点位置或第二图像多个初始同名像点位置，具有多个角度不同的直线。图11示出的是一个在第一图像中选取3个不同半径的同心圆与8条辐射状点阵的实施例，由图11可知，箭头左图为第一图像，第一图像上选取了3个为等差数列的不同半径，8条辐射状点阵间隔角相同；在第二图像上同名点集生成了对应的图像。这种与现有技术方案中直接选取矩形窗口不同的选取方式能扩大邻域、增加有效特征并降低计算量，建立的局部相关矩阵能较好地具有代表性地覆盖整个图像的像点情况。

并且，对局部相关矩阵使用辐射几何改正的最小二乘平差算法包括：

计算两幅相邻图像各自的灰度值、均值和标准差，并根据两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正；

计算两幅相邻图像各自的形变函数与随机噪声，并根据两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正，进行最小二乘平差计算。

一方面，计算两幅相邻图像各自的灰度值、均值和标准差，并根据两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正包括：

获取第一图像的灰度值、灰度均值与标准差以及第二图像的灰度值、灰度均值与标准差

根据第一图像的灰度值、灰度均值与标准差以及第二图像的灰度值、灰度均值与标准差使用中心化方法对第一图像与第二图像的灰度值g_l,g_r进行标准化处理，得到 $\left\{\begin{array}{c}{g}_l^{\′}=\frac{g_l-\stackrel{\‾}{g_l}}{{\mathrm{\σ}}_{g_l}}\\ g_r^{\′}=\frac{g_r-\stackrel{\‾}{g_r}}{{\mathrm{\σ}}_{g_r}}\end{array}\right.,$ 其中，g′_l,g′_r分别为标准化处理后的两幅图像的灰度值；

将第一图像的灰度值g_l与第二图像的灰度值g_r使用线性变换进行匹配，有g′_l＝h₁+h₂g′_r，可获得变换系数h₁,h₂，变换系数h₁,h₂可对辐射畸变进行匹配纠正。

另一方面，计算两幅相邻图像各自的形变函数与随机噪声，并根据两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正包括：

获取第一图像的形变函数g_l(x,y)与随机噪声函数n_l，以及第二图像的形变函数g_r(x,y)与随机噪声函数n_r；

根据第一图像与第二图像在匹配后的形变与噪声一致性，可得g_l(x,y)+n_l＝g_r(x,y)+n_r，并计算出误差方程v＝n_r-n_l＝g_l(x,y)-g_r(x,y)；

根据误差方程与几何畸变误差最小原则进行几何畸变的匹配纠正。

在一个较佳实施例中，对第一图像与第二图像同时进行辐射畸变与几何畸变的匹配纠正的方法如下所示。以下所有变量中，下标为1的变量意为指代第一图像的对应变量，下标为2的变量意为指代第二图像的对应变量。

首先，立体图像仿射变换包括像对移位与旋转变换如下：

$\left\{\begin{array}{c}{x}_2=a_0+a_{1}x+a_{2}y\\ y_2=b_0+b_{1}x+b_{2}y\end{array}\right.$

其中，x,y为窗口图像中心原点坐标，且x＝x₁,y＝y₁；a₀,a₁,a₂,b₀,b₁,b₂均为变换参数。

对第一图像与第二图像进行兼顾几何畸变与辐射畸变的纠正时应该满足：

g₁(x,y)+n₁(x,y)＝h₀+h₁g₂(x₂,y₂)+n₂(x₂,y₂)

＝h₀+h₁g₂(a₀+a₁x+a₂y,b₀+b₁x+b₂y)+n₂(x,y)

由上式可得误差方程如下：

v＝c₀δh₀+c₁δh₁+c₂δa₀+c₃δa₁+c₄δa₂+c₅δb₀+c₆δb₁+c₇δb₂-Δg

其中， $\left\{\begin{array}{c}{c}_0=1\\ c_1=g_2\\ \frac{{\∂g}_2}{{\∂x}_2}\·\frac{{\∂x}_2}{{\∂a}_0}={\left({\stackrel{.}{g}}_2\right)}_x=c_2\\ \frac{{\∂g}_2}{{\∂x}_2}\·\frac{{\∂x}_2}{{\∂a}_1}={\left({\stackrel{.}{g}}_2\right)}_{x}x=c_3\\ \frac{{\∂g}_2}{{\∂x}_2}\·\frac{{\∂x}_2}{{\∂a}_2}={\left({\stackrel{.}{g}}_2\right)}_{x}y=c_4\end{array}\right.$ $\left\{\begin{array}{c}\frac{{\∂g}_2}{{\∂y}_2}\·\frac{{\∂y}_2}{{\∂b}_0}={\left({\stackrel{\·}{g}}_2\right)}_y=c_5\\ \frac{{\∂g}_2}{{\∂y}_2}\·\frac{{\∂y}_2}{{\∂b}_1}={\left({\stackrel{\·}{g}}_2\right)}_{y}x=c_6\\ \frac{{\∂g}_2}{{\∂y}_2}\·\frac{{\∂y}_2}{{\∂b}_2}={\left({\stackrel{\·}{g}}_2\right)}_{y}y=c_7\end{array}\right.$

其中，δh0与δh1是辐射改正系数；δa0、δa1与δa2分别为x方向的位移、比例与旋转改正系数；δb0、δb1、δb2分别为y方向的位移、比例与旋转改正系数；与分别为第二图像在x、y方向上的偏导数，由于第二图像是规则排列的等间隔数字阵列，求偏导数可以用差分替代以简化运算。

按以上系数建立两图像窗口的误差方程的矩阵表达形式为：

$\underset{n*n,1}{V}=\underset{n*n,8}{C}\·\underset{\mathrm{8,1}}{X}-\underset{n*n,1}{L}$

其匹配模型参数初始值为：

$\left\{\begin{array}{c}{h}_0=0,h_1=1,a_0=x_2-x,a_1=1\\ a_2=0,b_0=y_2-y,b_1=0,b_2=1\end{array}\right.$

其辐射畸变改正参数计算式为：

$\left\{\begin{array}{c}{h}_0^i=h_0^{i-1}+{\mathrm{\δh}}_0^i+h_0^{i-1}{\mathrm{\δh}}_1^i\\ h_1^i=h_1^{i-1}+h_0^{i-1}{\mathrm{\δh}}_1^i\end{array}\right.$

其几何畸变改正参数计算式为：

$\left\{\begin{array}{c}{a}_0^i=a_0^{i-1}+{\mathrm{\δa}}_0^i+a_0^{i-1}{\mathrm{\δa}}_1^i+b_0^{i-1}{\mathrm{\δa}}_2^i\\ a_1^i=a_1^{i-1}+a_1^{i-1}{\mathrm{\δa}}_1^i+b_1^{i-1}{\mathrm{\δa}}_2^i\\ a_2^i=a_2^{i-1}+a_2^{i-1}{\mathrm{\δa}}_1^i+b_2^{i-1}{\mathrm{\δa}}_2^i\\ b_0^i=b_0^{i-1}+{\mathrm{\δb}}_0^i+b_0^{i-1}{\mathrm{\δb}}_1^i+b_0^{i-1}{\mathrm{\δb}}_2^i\\ b_1^i=b_1^{i-1}+b_1^{i-1}{\mathrm{\δb}}_1^i+b_1^{i-1}{\mathrm{\δb}}_2^i\\ b_2^i=b_2^{i-1}+a_2^{i-1}{\mathrm{\δb}}_1^i+b_2^{i-1}{\mathrm{\δb}}_2^i\end{array}\right.$

迭代计算图像匹配的辐射畸变改正参数与几何畸变改正参数(共计8个)，比较前后两次之间两图像窗口改正值中误差是否继续减小，如果减小则继续迭代，否则终止迭代并输出匹配模型参数。

完成了图像目标窗口与搜索匹配窗口的基于辐射几何改正的匹配模型计算，得到所有搜索窗口与目标窗口的最佳匹配窗口后，即可以认为两窗口中心为对应的共轭点，但在高精度影像匹配中，应该认定图像特征最明显的位置为基于图像特征的最佳匹配点所在，而图像梯度越大，特征越丰富，则匹配越可靠。因此，先对目标图像窗口取梯度加权值求得目标窗口的最佳匹配点坐标，得到高精度匹配点位。计算公式如下：

$\left\{\begin{array}{c}{x}_t=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{x}_{i,j}\frac{{\stackrel{\·}{g}}_x^2(i,j)}{\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\stackrel{\·}{g}}_x^2(i,j)}\\ y_t=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{y}_{i,j}\frac{{\stackrel{\·}{g}}_y^2(i,j)}{\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\stackrel{\·}{g}}_y^2(i,j)}\end{array}\right.$

再根据上式得到的对应仿射变化参数，可求得亚像素的同名像点位置。

综上所述，借助于本发明的上述技术方案，通过提取出图像特征点并使用旋转平移不变特征向量进行初匹配并构建Delaunay三角网获取仿射变换参数以建立相关矩阵进而得到密集匹配结果的技术方案，能够在纹理贫乏、相邻局部区域相似度较高的三维区域中进行工作，使用更少的运算量得到更高精度水平的密集匹配结果。具体地，本发明提供的技术方案适用于具有重叠成像区域的两任意光学图像之间非特征点的密集匹配，充分利用了临近局部区域越大同名特征越多、距离中心位置越近同名特征作用越显著在图像相关中的基本概念，提高了对同名点进行密集匹配的效率与质量；在扩大有效相关区域的同时，尽可能降低了实际参与相关计算的二维矩阵维度，降低了计算量，提高了对单位图像进行处理的速度。

从初步实验结果可知，本发明提供的技术方案提出的密集匹配方法无需人工特征干预，即可适应沙漠等纹理贫乏区域；对于一般地形区域的图像，可达到0.3像素的密集匹配精度，匹配准确度可达99.5％，由于计算量降低了，时间消耗反而还比一般的核线搜索和矩形窗口相关的计算要快；对于建筑物密集的区域，有效匹配的数量比常规方法高出近10倍。

所属领域的普通技术人员应当理解：以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种图像密集匹配方法，其特征在于，包括：

获取立体观测区域的框幅式序列图像；

从每幅所述序列图像中选择两幅相邻的图像，提取出所述两幅相邻图像各自的特征点图像，并根据所述特征点图像分别构建旋转平移不变特征向量，所述两幅相邻图像包括第一图像与第二图像；

在核线约束下，使用随机抽样一致算法对所述旋转平移不变特征向量进行同名初匹配，得到具有稳定相对关系的最大同名点对子集；

根据所述最大同名点对子集构建所述两幅相邻图像平面的第一Delaunay三角网，计算所述第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实所述各对同名三角形在所述两幅相邻图像中的拓扑关系是否准确；

根据所述各对同名三角形在所述两幅相邻图像中的拓扑关系是否准确的核实结果，移除所有拓扑关系不准确的同名点，构建第二Delaunay三角网，并计算所述第二Delaunay三角网中各对同名三角形之间的仿射变换参数；

获取密集匹配所需的平面像素间距，根据所述平面像素间距计算密集匹配在所述第一图像的各非特征像点位置，并根据所述各非特征像点与所述第二Delaunay三角网判断所述各非特征像点与所述第二Delaunay三角网中三角形的拓扑包容关系，根据所述拓扑包容关系获得各像点位置的初始点应调用的所述同名三角形的多个仿射变换参数，并计算多个所述仿射变换参数在所述第二图像中的多个初始同名像点位置；

以所述第一图像当前的各非特征像点位置和所述第二图像多个初始同名像点位置为中心，选取多个指定半径的同心圆与辐射状点阵，建立所述两幅相邻图像的局部相关矩阵；

对所述局部相关矩阵使用辐射几何改正的最小二乘平差算法，迭代计算得到具有最大相关系数的亚像素像点位置的集合，所述像点位置的集合为所述两幅相邻图像的密集匹配结果；

使用上述方法计算出所述框幅式序列图像中任意两幅相邻图像的密集匹配结果，确定所述密集匹配结果中空间参考平面与立体图像的内外参数，并根据所述密集匹配结果、空间参考平面与立体图像的内外参数对所述立体观测区域进行三维形貌重建。

2.根据权利要求1所述的一种图像密集匹配方法，其特征在于，提取出所述两幅相邻图像各自的特征点图像包括：

将所述两幅相邻图像的每一幅划分为多个大小相等的像元；

在所述每一个像元中的相同位置任取一点，计算该点在该点所在的像元中上下左右四个方向的灰度差分绝对值；

获取事先确定的阈值，并根据所述点在所述点所在的像元中上下左右四个方向的灰度差分绝对值确定协方差矩阵与误差椭圆，并根据所述协方差矩阵与误差椭圆确定每个所述像元中一初选点与一备选点；

根据所述每个所述像元中的一初选点与一备选点，确定每个所述像元中的特征点；

将所有像元中的特征点组合生成两幅相邻图像的两幅特征点图像。

3.根据权利要求1所述的一种图像密集匹配方法，其特征在于，计算所述第一Delaunay三角网中各对同名三角形的三顶点坐标构成的二维向量的叉积，并根据该叉积的正负核实所述各对同名三角形在两图像中的拓扑关系是否准确包括：

根据所述第一Delaunay三角网中某一对同名三角形Δp₁p₂p₃与Δp₁′p₂′p₃′的三顶点坐标 $\left\{\begin{array}{c}{p}_1(x_1,y_1),p_2(x_2,y_2),p_3(x_3,y_3)\\ p_1^{\′}(x_1^{\′},y_1^{\′}),p_2^{\′}(x_2^{\′},y_2^{\′}),p_3^{\′}(x_3^{\′},y_3^{\′})\end{array}\right.,$ 计算获得所述顶点坐标构成的平面矢量

根据所述顶点坐标构成的平面矢量,计算获得所述该对同名三角形Δp₁p₂p₃与Δp₁′p₂′p₃′的矢量叉积面积与

若判定所述拓扑同名关系准确、

若判定所述拓扑同名关系不准确；

依次对所述第一Delaunay三角网中的每一对同名三角形进行上述计算，并获得判定结果。

4.根据权利要求3所述的一种图像密集匹配方法，其特征在于，据所述各非特征像点与所述第二Delaunay三角网判断所述各非特征像点与所述第二Delaunay三角网中三角形的拓扑包容关系包括：

过所述非特征像点p₀(x₀,y₀)作射线

对于所述第二Delaunay三角网中的某一个三角形Δp₁p₂p₃，依次判断有向线段与所述射线是否存在交点；

统计所述射线与所述有向线段的交点个数总数并进行判断：

若所述交点个数总数为奇数，则所述非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃内、

若所述交点个数总数为偶数，则所述非特征像点p₀(x₀,y₀)位于三角形Δp₁p₂p₃外。

5.根据权利要求4所述的一种图像密集匹配方法，其特征在于，所述依次判断有向线段与所述射线是否存在交点包括；

对于

首先，当y₀≤y₁且y₀≤y₂时，或者当y₀≥y₁且y₀≥y₂时，判定所述有向线段与所述射线不存在交点；

其次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≤y₂时，计算Δp₁p₂p₀的面积其中

若判定所述射线与所述有向线段存在交点、

若表示所述射线与所述有向线段不存在交点；

再次，当y₀≤y₁且y₀≤y₂不成立、且y₀≥y₁且y₀≥y₂不成立、且y₁≥y₂时，计算Δp₁p₂p₀的面积其中

若判定所述射线与所述有向线段存在交点、

若表示所述射线与所述有向线段不存在交点；

对于与与所述射线是否存在交点的判定方式与相同。

6.根据权利要求5所述的一种图像密集匹配方法，其特征在于，对所述局部相关矩阵使用辐射几何改正的最小二乘平差算法包括：

计算所述两幅相邻图像各自的灰度值、均值和标准差，并根据所述两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正；

计算所述两幅相邻图像各自的形变函数与随机噪声，并根据所述两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正，进行最小二乘平差计算。

7.根据权利要求6所述的一种图像密集匹配方法，其特征在于，计算所述两幅相邻图像各自的灰度值、均值和标准差，并根据所述两幅相邻图像各自的灰度值、均值和标准差进行辐射畸变的匹配纠正包括：

获取所述第一图像的灰度值、灰度均值与标准差以及所述第二图像的灰度值、灰度均值与标准差

根据所述第一图像的灰度值、灰度均值与标准差以及所述第二图像的灰度值、灰度均值与标准差使用中心化方法对所述第一图像与所述第二图像的灰度值g_l,g_r进行标准化处理，得到 $\left\{\begin{array}{c}{g}_l^{\′}=\frac{g_l-\stackrel{\‾}{g_l}}{{\mathrm{\σ}}_{\mathrm{gl}}}\\ g_r^{\′}=\frac{g_r-\stackrel{\‾}{g_r}}{{\mathrm{\σ}}_{\mathrm{gr}}}\end{array}\right.,$ 其中，g_l′_,g_r′分别为标准化处理后的所述两幅图像的灰度值；

将所述第一图像的灰度值g_l与所述第二图像的灰度值g_r使用线性变换进行匹配，有g_l′＝h₁+h₂g_r′_，可获得变换系数h₁,h₂，所述变换系数h₁,h₂可对辐射畸变进行匹配纠正。

8.根据权利要求6所述的一种图像密集匹配方法，其特征在于，计算所述两幅相邻图像各自的形变函数与随机噪声，并根据所述两幅相邻图像各自的形变函数与随机噪声计算出几何畸变的误差方程，再按照几何畸变误差最小原则进行几何畸变的匹配纠正包括：

获取所述第一图像的形变函数g_l(x,y)与随机噪声函数n_l，以及所述第二图像的形变函数g_r(x,y)与随机噪声函数n_r；

根据所述第一图像与所述第二图像在匹配后的形变与噪声一致性，可得g_l(x,y)+n_l＝g_r(x,y)+n_r，并计算出误差方程v＝n_r ^-n_l＝g_l(x,y)-g_r(x,y)；

根据所述误差方程与几何畸变误差最小原则进行几何畸变的匹配纠正。