一种缩小多视影像匹配同名像点搜索范围的方法与流程

本发明属于数字摄影测量、地理信息系统和计算机视觉技术，具体涉及一种缩小多视影像匹配同名像点搜索范围的方法。

背景技术：

摄影测量作为利用影像研究被摄物体的形状、位置、大小、特性及相互位置关系的一门科学，一直是三维地理空间信息获取的主要技术手段，并被广泛应用于基础测绘、地理国情监测、资源环境调查、城市规划和智慧城市建设等国民经济建设和社会发展中的众多领域。在摄影测量的数据后处理中，多视影像匹配是一个关键步骤，其目的是自动寻找同一地物在多幅影像上的同名像点。而多视影像匹配获取的同名像点则是自动空中三角测量、地理信息采集、三维场景重建等摄影测量数据处理环节的重要数据基础。另外，多视影像匹配在视频数据处理、多源数据融合、机器人视觉导航等应用中也具有重要的作用。因此，影像匹配一直是数字摄影测量与计算机视觉领域的热点研究问题。

多视影像匹配包含两个基本问题：匹配相似性测度的计算与同名像点搜索范围的确定。匹配测度是判断多幅影像上的不同像点是否为同名像点的依据，会对多视影像匹配的准确率与稳健性产生重要的影响。同名像点搜索范围决定了参与匹配计算的候选同名像点的数量，其值的大小会对多视影像匹配的可靠性与计算效率产生重要的影响。目前，现有影像匹配方法的研究重点是如何利用局部匹配窗口内的像方信息(灰度值，特征向量等)和引入各种约束条件(核线，平面单应性，像点之间的邻接关系等)，来提高匹配测度的稳健性与约束同名像点的搜索范围。但是，像方信息不具有唯一性和不变性特点，无法构建能唯一地表达同名像点的独特特征。另外，由于遮挡、几何变形等因素的影响，匹配窗口中的像点也常常不满足平面单应性、空间关系保持等约束条件。因此，现有影像匹配方法无法克服影像中“同物异谱”或“同谱异物”现象带来的匹配测度歧义性问题，及地物遮挡与透视失真因素造成的像点间空间关系变化问题；在处理影像中纹理重复、弱纹理、视差不连续、几何变形等匹配困难区域时，存在着匹配率低、匹配结果可靠性差等缺陷。

那么，多视影像中的同名像点位置是否存在某种独特的规律？是否有某种数值指标能独特地指示同名像点的位置？如何定量描述同名像点的这种位置规律？这种规律是否可以用于进一步精化多视影像匹配过程中的同名像点的搜索范围或唯一地确定同名像点？这些问题的研究对于解决现有影像匹配方法的缺陷具有重要的意义。

技术实现要素：

发明目的：本发明的目的是针对多视影像匹配过程中的同名像点精确搜索问题，基于测量平差系统的可靠性理论，提供一种缩小多视影像匹配同名像点搜索范围的方法，以得到反映多视影像中同名像点位置特征的数值指标，从而为大幅降低候选同名像点的搜索个数提供方法支撑。

技术方案：一种缩小多视影像匹配同名像点搜索范围的方法，包括如下步骤：

(1)对于基准影像上的待匹配像点通过影像的几何成像模型和物方高程信息，确定其在搜索影像上的具有统一编组规则的初始候选同名像点集合；

(2)以待匹配像点及其初始候选同名像点集合中的所有候选点为观测数据，基于多视前方交会原理，进行针对同名像点未知物方三维坐标的最小二乘平差计算，并根据多视前方交会的平差结果，计算各组候选同名像点的标准化残差；

(3)根据候选像点的标准化残差分布图、及待匹配像点的实际同名像点位置，对初始候选同名像点集合进行迭代精化处理；

(4)根据步骤(3)各次迭代计算的结果，确定反映多视影像中同名像点位置特征的定量化数值指标，以及区间长度大幅缩小的同名像点搜索范围。

进一步的，步骤(1)中所述基准影像为影像集合中的任意一幅影像，所述搜索影像为影像集合中除基准影像后的影像，所述的待匹配像点为基准影像中的任意像点。

步骤(1)包括如下步骤：

(11)输入具有外方位元素的n幅多视影像(若为线阵推扫式成像的航天摄影影像，则可视有理多项式系数为影像的外方位元素)，及影像覆盖区域的物方概略高程范围，所述的物方概略高程范围表示为：[最大高程值zmax，最小高程值zmin]；

(12)对于基准影像i0上待匹配像点p，根据影像区域的概略高程范围，利用影像的几何成像模型确定其在其余各幅搜索影像ik的同名核线上的m^k个候选同名像点q^ki；所述的几何成像模型表达式如下：

(r,c)＝f(x,y,z)

其中，r、c为像点在影像中的行号与列号，x、y、z为像点对应地物点的物方三维坐标；f为表达影像几何成像模型的函数，k＝1,2,…,n-1；i＝1,2,…,m^k；对于中心透视投影的航空或近景影像，使用共线条件方程模型，对于多线阵推扫式的航天卫星影像，使用有理函数模型；

(13)基于物方坐标约束，对各幅搜索影像上的大小不同的候选同名像点数量进行一致化处理，按物方坐标将各幅搜索影像上相互独立的各个候选同名像点统一编组为mz组初始像点集合，mz为m^k中的最大值，所述初始像点集合为：

{[b¹1，b²1，...，b^n-11]，[b¹2，b²2，...，b^n-12]，...，[b¹j，b²j，...，b^n-1j]，...，[b¹mz，b²mz，...，b^n-1mz]}(j∈{1，2，...，mz})。

进一步的，步骤(2)包括如下步骤：

(21)对于步骤(1)所确定的mz组候选同名像点集合，根据其中的各个候选同名像点的行号与列号数值，剔除每幅搜索影像上的重复像点；

(22)以待匹配像点及其剔除重复点后的所有候选同名像点为观测像点，基于影像的几何成像模型，构建求解像点未知物方三维坐标近似值改正数x3×1＝[dxdydz]^t的误差方程矩阵v2m×1＝a2m×3x3×1-l2m×1；

(23)基于最小二乘原理，对前述的误差方程矩阵进行平差计算，并在平差结果的基础上，计算mz组中第j组候选像点的标准化残差值其中，ω(b^kj)为第j组候选点中的第k幅搜索影像ik上的候选同名像点b^kj的标准化残差值。

进一步的，步骤(3)包括如下步骤：

(31)绘制mz组候选同名像点的标准化残差的分布曲线图，再根据待匹配像点的实际同名像点位置，从初始候选同名像点集合中选择出包含真正同名像点的nz组像点集合(nz≤mz)：{[b¹1，b²1，...，b^n-11]，[b¹2，b²2，...，b^n-12]，...，[b¹j，b²j，...，b^n-1j]，...，[b¹nz，b²nz，...，b^n-1nz]}(j∈{1，2，...，nz})，作为新的候选同名像点集合；

(32)若nz≠mz，则以nz组新候选同名像点集合作为新一次平差计算的观测像点，重新进行步骤(2)的候选点集的标准化残差计算，并再次挑选包含真正同名像点的t组新候选点集合；此过程迭代进行，直至新候选同名像点集合中的像点数量不再变化为止；

(33)上述迭代计算过程结束后，输出的t值为2。

有益效果：与现有技术相比，本发明所述方法获得了反映多视影像匹配过程中待匹配像点的真正同名像点与其它候选像点之间分布特征的独特误差指标，非常有利于同名像点搜索范围的进一步精确确定，另一方面该方法有望在算法原理层面为缩小多视影像匹配的同名像点搜索范围提供一种新的研究思路，从而大幅减少影像匹配过程中候选同名像点的搜索个数，提高影像匹配的计算效率和准确率。

附图说明

图1是本发明实施例的方法框架图；

图2(a)是趋向水平形态下的基于物方约束的多视影像初始候选同名像点确定示意图；

图2(b)是趋向竖直形态下的基于物方约束的多视影像初始候选同名像点确定示意图；

图3是本发明实施例的基于物方约束的多视影像初始候选同名像点统一编组示意图；

图4(a)是实施例中位于基准影像的地面上的待匹配像点，及其在两幅搜索影像上的实际同名像点位置图；

图4(b)是位于基准影像的房顶上的待匹配像点，及其在两幅搜索影像上的实际同名像点位置图；

图4(c)是位于基准影像的房屋墙面上的待匹配像点，及其在两幅搜索影像上的实际同名像点位置图；

图5是针对图4(a)中的待匹配像点，计算的所有候选同名像点的误差指标分布曲线图；

图6是针对图4(b)中的待匹配像点，计算的所有候选同名像点的误差指标分布曲线图；

图7是针对图4(c)中的待匹配像点，计算的所有候选同名像点的误差指标分布曲线图。

具体实施方式

为了详细的说明本发明所公开的技术方案，下面结合说明书附图及具体实施例做进一步的阐述。

本发明所公开的是一种缩小多视影像匹配同名像点搜索范围的方法，该方法首先利用影像的成像模型和先验的物方高程信息，确定基准影像上的待匹配像点在各幅搜索影像上的具有统一编组规则的初始候选同名像点集合；其次以待匹配像点及其初始候选同名像点集合中的所有候选点为观测数据，进行针对同名像点未知物方三维坐标的多视前方交会平差计算，并根据平差结果计算各组候选同名像点的标准化残差；最后根据候选像点的标准化残差分布图、及待匹配像点的实际同名像点位置，对初始候选同名像点集合进行迭代精化处理，从而根据各次迭代计算的结果，获得可以反映多视影像中同名像点位置特征的定量化数值指标，以及区间长度大大缩小的同名像点搜索范围。本发明方法可用于多视影像匹配的同名像点精确搜索，从而大幅减少影像匹配过程中候选同名像点的搜索个数，提高多视影像密集匹配的计算效率和准确率。

本发明的基本计算过程为：

(1)从已知内外方位元素的n幅多视航空影像中，人工确定一幅作为基准影像i0，余下的n-1幅作为搜索影像i1、i2、…、in-1，并从基准影像和搜索影像上选择一个待匹配像点及其对应的n-1个同名像点；

(2)对于给定的待匹配像点，基于核线和物方概略高程范围约束，确定其在各幅搜索影像上的候选同名像点搜索范围；并对这些不同的搜索范围进行一致化处理，以使不同搜索影像上的离散候选同名像点统一编组，从而形成初始的m组候选同名像点集合；

(3)利用待匹配像点及其m组候选同名像点的像方像平面坐标，基于多视前方交会原理，进行针对像点的未知物方三维坐标的最小二乘平差计算；再根据平差结果，计算各组候选同名像点的标准化残差，并绘制候选像点标准化残差分布图；

(4)根据标准化残差分布图、及待匹配像点的实际同名像点位置，从m组候选同名像点集合中选择出包含真正同名像点的新的n组候选同名像点集合；

(5)若n≠m，则将新的n组候选同名像点作为初始候选同名像点集合，再次进行步骤(3)的标准化残差计算和步骤(4)的新候选同名像点集合确定；

(6)迭代进行步骤(5)的计算过程，直至新候选同名像点集合中的像点数量不再变化为止；迭代结束时输出的t组候选同名像点，即是根据标准化残差精化后的候选同名像点集合，而且t会远远小于m。

如附图1所示，一种缩小多视影像匹配同名像点搜索范围的方法主要包括三个部分：

(a)确定基准影像上的待匹配像点在各幅搜索影像上的具有统一编组规则的多组候选同名像点集合；

(b)基于测量平差原理的各组候选同名像点的标准化残差计算；

(c)基于标准化残差的候选点集合的迭代精化处理。具体实施步骤为：

第一步：确定待匹配点在各幅搜索影像上的具有统一编组规则的多组候选同名像点集合。

基于核线和物方概略高程范围等物方约束，确定待匹配像点在各幅搜索影像上的具有统一编组规则的多组候选同名像点集合的具体流程如下：

(1)从已知内外方位元素的n幅多视航空影像中，人工确定一幅作为基准影像i0，余下的n-1幅作为搜索影像i1、i2、…ik、…、in-1(k＝1,2,…,n-1)，基准影像和各幅搜索影像的摄影中心分别为s0、sk，相应的外方位元素分别为并从基准影像和搜索影像上选择一个待匹配像点p及其对应的n-1个同名像点(待匹配像点p的像方行列号记为(r,c)，像平面坐标记为(x,y))。

(2)根据先验知识，可确定基准影像所覆盖地区的近似高程范围：最小高程zmin与最大高程zmax(这个高程范围并不需要很精确，只要能包含待匹配像点的实际高程就行)，则像点p对应的地物点p一定位于摄影光线s0p上的线段pminpmax之间。利用影像几何成像模型(r,c)＝f(x,y,z)的逆算式(见式(1)，其以中心投影影像的共线方程模型为例)，由像点的像平面坐标和物方高程计算出点pmin、pmax的物方平面坐标(xmin,ymin)、(xmax,ymax)。

式中，f是内方位元素中的相机主距，是由影像i0的外方位角元素确定的旋转矩阵中的九个方向余弦。

(3)根据影像的几何成像模型(见式(2)，以中心投影影像的共线方程模型为例)，将点pmin(xmin,ymin,zmin)、pmax(xmax,ymax,zmax)往搜索影像ik上进行投影，得到对应像点及其像平面坐标与则待匹配点p在影像ik上的同名像点p^k一定位于线段上(这个线段实际上就是利用外方位元素确定的同名核线，如附图2所示)，此同名核线上的所有m^k个像点即是待匹配像点p在影像ik上的初始候选同名像点。

式中，是由搜索影像ik的外方位角元素确定的旋转矩阵中的九个方向余弦。

利用内方位元素中的像主点坐标与像元尺寸，可计算出像点的像方行列号与从而得到同名核线的直线方程，并由此核线方程可计算核线上第i个候选点的行列号同名核线的直线方程按如下两种情况计算：

①当时，核线趋向于水平形态(如附图2(a))，其直线方程为其中，斜率截距自变量的取值范围是候选同名像点的总数int()是取整函数，abs()是取绝对值函数,min()是取最小值函数，max()是取最大值函数。

②当时，核线趋向于竖直形态(如附图2(b))，其直线方程为其中，斜率截距自变量的取值范围是候选同名像点的总数

(4)在步骤(3)中，待匹配点p在各幅搜索影像上的初始候选同名像点数量m^k是不同的，且不同搜索影像上的候选点也是互相独立，难以根据它们之间的序列号建立联系。为了便于分析候选同名像点的分布规律，采用如下的方法对各幅搜索影像上的候选同名像点数量进行一致化处理，从而将相互独立的各个离散候选同名像点按物方坐标规则统一编组为不同的候选点集合：

①以候选同名像点数量mz＝max(m^k|k＝1,2,…,n-1)最大的一幅搜索影像为主搜索影像(如附图3中的影像i1)，剩余的n-2幅影像为副搜索影像。

②对于主搜索影像中同名核线b¹1b¹mz上的第j(j∈{1,2,…,mz})个候选同名点先利用双像前方交会方法，由基准影像上的光线s0p与主搜索影像上的光线计算出对应地物点pj的物方三维坐标(xj,yj,zj)；再利用式(2)所示的共线方程与第t幅副搜索影像it(t∈{2,…,n-1})的同名核线的直线方程，将点pj往每一幅副搜索影像进行投影(如附图4中的虚线stp1、…、stpj、…、stpmz)，在影像it上得到与点对应的候选同名点则，为所有搜索影像上与待匹配点对应的第j组候选同名像点。

③当主搜索影像上的所有候选同名像点都进行了步骤②的一致化处理，则待匹配像点p在n-1幅搜索影像上的候选同名像点个数都统一为mz，且相同序列号的候选点属于同一组。即，原来各幅搜索影像上数量各不相同的候选点被统一成如下的mz组：{[b¹1,b²1,…,b^n-11],[b¹2,b²2,…,b^n-12],…,

需要注意的是：经一致化处理后的mz组候选同名像点，除了主搜索影像外，其它某幅副搜索影像it上的mz个候选点是存在重复点的。

第二步：基于测量平差原理的各组候选同名像点的标准化残差计算。

对基准影像上的待匹配像点p，及其统一编组后的各组候选同名像点集合{[b¹i,b²i,…,b^n-1i]|i∈{1,2,…,mz}}，按如下过程进行基于测量平差原理的各组候选同名像点的标准化残差计算，从而得到第i组候选同名像点的标准化残差值ωi：

(1)去除重复的候选同名像点。由前述的搜索影像上候选同名像点的确定原理可知，第k幅搜索影像上mz个候选点中存在重复点，且重复候选点只可能出现在序列号连续的相邻点中。所以，从第2个候选点开始，逐个判断其与相邻的前一个点是否重复，若重复，则将其标记为重复点。当所有mz个候选点都进行判别标记后，可得到没有重复的mk个候选点另外，经重复点判别后，还可得到原来重复点集中点的序列号j(j∈{1,2,…,mz})与去重复点后点集中点的序列号j′(j′∈{1,2,…,mk})的对应关系：j′＝cor(j)(函数cor表示第j个点与第j′个点是同一个点)。此时，基准影像和n-1幅搜索影像上的同名像点的总数为

(2)构建所有像点的误差方程。根据影像的几何成像模型，以影像上像点的像平面坐标为观测值，对像点的未知物方三维坐标(x,y,z)进行泰勒级数线性化，可列出式(3)所示的求解物方三维坐标的误差方程(以式(2)所示的中心投影影像的共线方程为例)：

式中，(x，y)是像点的像平面坐标观测值，(vx，vy)为观测值的残差，((x)⁰，(y)⁰)是将像点的近似物方坐标(x⁰，y⁰，z⁰)带入式(2)所计算的近似像平面坐标，(dx，dy，dz)为待求解的物方三维坐标近似值的改正数。而误差方程式的6个系数(a11，a12，a13，a21，a22，a23)按下式计算：

式中，(a1，a2，a3；b1，b2，b3；c1，c2，c3)是影像的旋转矩阵中的九个方向余弦，(xs，ys，zs)是影像摄影中心的物方三维坐标。

根据式(3)，对于基准影像和各幅搜索影像上的所有m个同名像点，按先基准影像再搜索影像的顺序，逐个列出其上像点的误差方程，形成如式(5)所示的误差方程矩阵。

v2m×1＝a2m×3x3×1-l2m×1(5)

式中，v为2m个像平面坐标观测值的残差矩阵(维数为2m行×1列)，x为待求解的物方三维坐标的近似值改正数矩阵(维数为3行×1列)，a、l分别为误差方程的系数矩阵和常数项矩阵(维数分别为2m行×3列、2m行×1列)。这些矩阵的具体形式见式(6)。

式中，矩阵中各个元素的上标表示影像编号，下标表示数据项和点的序号。

(3)计算可靠性矩阵。根据误差方程中系数阵a、常数阵l，及观测值的权阵pll(本文的观测值为相互独立观测，此矩阵为单位矩阵e)，可计算出如下的可靠性矩阵r：

(4)解算误差方程。利用最小二乘原理，解算式(5)所示的误差方程矩阵，得到未知数矩阵x，所有观测值的残差矩阵v，及单位权中误差σ0：

(5)计算m个同名像点中第i"个点的标准化残差ω(i")。由于每个像点有x、y两个像平面坐标观测值，取像点的x、y的标准化残差的平均值作为该点的标准化残差(见式(9))。

(6)计算第k幅搜索影像ik上的第j个候选同名像点的标准化残差值由于点的序列号j对应于去重复点后的序列号j′，且从基准影像至第k-1幅搜索影像，去重复后的同名像点的总数为所以，的计算方式如下：

(7)计算第i组候选同名像点的标准化残差值ωi。由于第i组候选同名像点[b¹i,b²i,…,b^n-1i]是由n-1幅搜索影像上的n-1个候选像点组成，所以，按下式取n-1个候选像点的标准化残差的平均值作为ωi：

第三步：基于标准化残差的候选点集合的迭代精化处理。

在第二步计算出待匹配像点p在n-1幅搜索影像上的mz组候选同名像点集合{[b¹1,b²1,…,b^n-11],[b¹2,b²2,…,b^n-12],…,[b¹i,b²i,…,b^n-1i],…,}中每组像点的标准化残差值ωi(i∈{1,2,…,mz})后，按如下过程进行候选点集合的迭代精化处理：

(1)根据每组像点的标准化残差值ωi，以组号i为横轴，以标准化残差值ωi为纵轴，绘制出mz组候选同名像点的标准化残差分布曲线图。

(2)根据mz组候选同名像点的标准化残差分布曲线图的形状特点，及待匹配像点p在n-1幅搜索影像上的实际同名像点位置，若候选像点的标准化残差分布曲线呈现有最低点的近似v形分布形状，则从mz组候选同名像点集合中选择出包含真正同名像点的nz组像点集合(nz≤mz)：{[b¹1,b²1,…,b^n-11],[b¹2,b²2,…,b^n-12],…,[b¹nz,b²nz,…,b^n-1nz]}(j∈{1,2,…,nz})，作为新的候选同名像点集合。

(3)若nz≠mz，则对nz组新候选同名像点集合，重新进行第二步的候选点集的标准化残差计算，并再次挑选包含真正同名像点的t组新候选点集合；若t≠nz，则继续进行候选点集的标准化残差计算，及挑选下一次的t组新候选点集合；以上过程迭代进行，直至新候选同名像点集合中的像点数量不再变化为止。

(4)上述迭代计算过程结束时输出的t组候选同名像点，即是根据标准化残差迭代精化后的候选同名像点集合。

根据以上步骤，选择了附图4(a)-图4(c)的实验数据进行了本发明方法的实验验证。附图4展示了位于地面上、建筑物顶部、建筑物墙立面上的3个待匹配像点，及其在2幅搜索影像上的同名像点位置(搜索影像上的直线为同名核线)；对附图4中的3个待匹配像点进行本发明的候选同名像点集合的迭代精化计算，各次迭代计算的候选同名像点的标准化残差分布曲线如附图5、6、7所示。

根据附图5、6、7的结果，本发明总结出多视影像中同名像点的位置特征为：①对于多视影像中某一幅影像上的待匹配像点，其在剩余各幅影像同名核线上搜索区间内的所有候选同名像点的标准化残差值呈现“v”字形分布现象。②在所有候选同名像点的标准化残差值分布曲线上，最小残差值点基本处于所有候选点的中间位置，而且，真正的同名像点位于最小残差点和起点(或终点)之间。③从“v”字形分布曲线上，将包含真正同名像点的最小残差点左(或右)侧的局部点集作为新的候选点集合，迭代进行新点集的标准化残差的平差计算，直至新候选点集合内的像点数量不再变化为止。则每次迭代过程中，各候选同名像点的标准化残差值仍服从规律①、②的分布现象；而且，迭代结束后，每幅搜索影像上的候选同名像点的数量会缩减为2个。

因此，从结果可以看出，本发明的方法可获得区分多视影像匹配过程中待匹配像点的真正同名像点与其它候选像点的独特误差指标及其分布特征。经过本发明方法的有限次数的迭代计算，基准影像上的待匹配像点在各幅搜索影像上的候选同名像点的数量最终会缩减为2个，这大幅减少了多视影像匹配过程中候选同名像点的搜索个数，非常有利于提高多视影像匹配(尤其是多视密集匹配)的计算效率和准确率。