顾及高阶结构特征的异源遥感影像同名点匹配方法与流程

本发明属于测绘科学与技术领域，涉及异源遥感影像同名点匹配方法，尤其涉及一种顾及高阶结构特征的异源遥感影像同名点匹配方法。

背景技术：

随着对地观测技术的飞速发展，多种类型的传感器如雨后春笋般出现，它们以更高的时间分辨率、空间分辨率记录着地球表面的变化，也使得覆盖同一地区的影像越来越多样。异源遥感影像的目标联合定位已成为摄影测量与遥感领域的重要研究方向之一。然而，无论是单点三维定位还是区域网平差，均需要快速而准确地找寻影像同名点，并量测其影像坐标。

对于影像匹配方法，相关学者进行了大量研究，主要分为三种：基于特征的匹配、基于区域的匹配以及基于整体影像的匹配。基于特征的匹配方法，理论上可以减少影像灰度差异的不利影响，具有更好的适应性。近年来，相关学者提出，如小波特征、形状上下文、梯度径向夹角金字塔直方图、SIFT(尺度不变特征转换，Scale-invariant feature transform，SIFT)、SURF(加速鲁棒特征，Speed Up Robust Features，SURF)等多种特征描述算子。然而，基于该类特征算子的匹配方法，具有计算量大、匹配成功率不高、匹配点分布不均匀等缺点；一些学者提出的基于边缘、区块特征的匹配算法，受影像辐射特性差异以及噪声的影响，提取同名特征十分困难，普适性不强，均不能满足遥感影像空中三角测量作业中对同名点位置分布均匀的要求。

基于区域的匹配方法通常将整个影像或影像子区的灰度或梯度信息进行匹配。这类算法是以重叠图像上含有相应图像的目标区和搜索区中的灰度作为匹配基础，利用某种相似性测度判定两个影像块的匹配程度。然而，无论是基于灰度区域的影像匹配，还是大部分的基于特征的影像匹配，都是基于单点或局部的影像匹配，无法有效保证匹配结果的可靠性。

基于整体影像的匹配方法，通过考虑匹配结果间固有的映射关系或者邻近匹配间的相容一致性，来达到提高匹配可靠性的目的。陈飒等提出基于Contoulet域Hausdorff距离和粒子群的多源遥感图像匹配方法，利用最小二乘Hausdorff距离建立点集之间的对应关系，增强匹配的可靠性。李雨谦等提出基于水平集的图像配准方法，通过映射函数将SAR图像特征和光学特征结合，构造能量泛函模型，采用水平集方法求解曲线演化方程，来实现异源影像配准。但是，对于遥感影像，受影像覆盖范围宽广、地形起伏等因素影响，同名点之间的关系难以用函数化映射关系进行有效表达，且由于匹配点的空间稀疏性，将造成相容一致假设不再完全成立，为此，直接将这些整体匹配方法引入到异源遥感影像匹配并不合适。

在遥感影像匹配方法方面，专利号是ZL201010242888.1的发明专利(一种稳健的高分辨率卫星影像连接点自动匹配方法)以及专利号是ZL201110091756.8的发明专利(基于RFM模型的多源星载SAR影像自动匹配方法)均涉及遥感影像同名点匹配。但其研究重点是将基于RFM模型的光束法区域网平差方法融入到同名点匹配过程中，利用RFM模型所提供的几何约束，提高特征点初始点位预测精度、提供核线几何约束缩小匹配搜索空间以及基于RFM模型区域网平差剔除错误匹配点，但其本质上仍属于单点匹配，匹配过程中未引入匹配点之间的高阶结构特征。

技术实现要素：

为了解决背景技术中存在的上述技术问题，本发明的目的在于提供一种能够提高异源遥感影像匹配的可靠性和成功率的顾及高阶结构特征的异源遥感影像同名点匹配方法，该方法将RFM模型所提供的几何约束以及匹配点之间的高阶结构特征引入到匹配过程中，实现同名点的自动、可靠、精确匹配，大大提高了同名点获取的自动化程度和工作效率，极大降低了人工工作量。

为实现上述目的，本发明采用如下技术方案：

一种顾及高阶结构特征的异源遥感影像同名点匹配方法，其特征在于：所述顾及高阶结构特征的异源遥感影像同名点匹配方法包括以下步骤：

1)异源影像匹配前的数据准备：所述异源影像匹配前的数据包括生成金字塔影像；所述金字塔影像包括原始影像层以及非原始影像层；所述非原始影像层以及原始影像层自上而下依次叠加并形成金字塔结构；所述原始影像层为第0层金字塔影像；

2)在金字塔影像上的原始影像层进行特征点簇的划分与提取：将原始影像层的整幅原始影像均匀分布一定大小的子区，在各子区内部，进一步划分若干均匀分布的格网；在每个格网范围内，采用摄影测量中常用的Forstner特征提取算子点的方法在每个格网内提取最佳特征点，形成特征点点簇；

3)特征点匹配候选点的确定：在每层金字塔影像上进行特征点初始点位的预测、建立核线几何约束方程、匹配窗口影像的几何粗纠正以及基于归一化互信息测度的匹配，确定每个特征点对应的匹配候选点；

4)顾及高阶结构特征的同名点匹配：待步骤3)中的每个特征点对应的匹配候选点的确定完成后，在每层金字塔影像上以特征点点簇为单元，同时考虑每个特征点所对应的匹配候选点的灰度相似性特征以及高阶的几何特征，通过基于超图的整体影像匹配方法确定每个特征点对应的同名点；所述每个特征点所对应的匹配候选点的灰度相似性特征是一阶特征；所述高阶的阶数不低于两阶；

5)在非原始影像层的每层金字塔影像中进行基于RFM模型的区域网平差，剔除错误的匹配候选点以及计算匹配候选点所对应物方坐标；

6)自上而下逐层完成金字塔影像的匹配直到完成原始影像层，最终实现异源遥感影像同名点的自动可靠匹配。

作为优选，本发明所采用的步骤4)的具体实现方式是：

4.1)构建超图模型；

4.2)对步骤1)所得到的超图模型进行超边采样；

4.3)根据采样结果对超图模型进行解算，由解算结果确定每个特征点对应的同名点。

作为优选，本发明所采用的步骤4.1)的具体实现方式是：

设某个特征点点簇为P¹，特征点个数为N₁，所述特征点所对应的匹配候选点集合为P²，特征点个数为N₂，表示集合P¹和P²的第i个特征点；集合P¹和P²对应的超图模型为G＝(V,E,A)和G′＝(V′,E′,A′)；

其中：

V和V′分别为对应的特征点顶点集合，即

为超边的集合，d为超边包含特征点个数；

A和A′为超边所对应的属性集；

设特征点点簇P¹对应的匹配候选点的顶点响应集合为C的一个k元组如公式(1)所示：

c_s1＝(v₁,v′₁)，…，c_s1＝(v_k,v′_k) (1)

其中：

(v₁,…,v_k)∈V；

(v′₁,…,v′_k)∈V′；

分别令超边e_1,…,k表示包含了顶点v₁,…,v_k的边，超边e′_1,…,k表示包含了顶点v′₁,…,v′_k的边；对于一个k阶图的超图匹配，通过比较两个度为k的超边的相似性来衡量k元组的相似性；指定k维相似性度量函数f_k，所述k维相似性度量函数f_k的参数为超图属性集合A和A′中的属性向量；可以将更低阶的超边同时考虑，此时k阶超边的相似度张量为：

T_si⁽¹⁾＝f₁(a_i,a′_i) (2)

式中：

参数γ^(k)表示k阶超边相似度的加权系数；

T的右上标(k)表示张量的维数；

a_1,…,k,a′_1,…,k分别表示k阶超图属性集合中的向量参数；

设超边最大阶数为δ，相似性张量T^δ为高阶张量，包含了所有不同度超边的相似性信息。

作为优选，本发明所采用的步骤4.2)的具体实现方式是：

对特征点点簇及特征点所对应的匹配候选点，采用k-d树(k-dimensional树)建立起空间索引，对特征点点簇的每个特征点，采用随机采样的方式，随机抽取一定数量的包含该特征点的3元组作为其超边；而在特征点匹配候选点所在影像上，并非采样所有的三元组，而是仅考虑由该特征点采样超边所包含特征点对应匹配候选点构成的3元组，且3元组中任意两个点均不重复。

作为优选，本发明所采用的步骤4.3)的具体实现方式是：

超图匹配，是在顶点响应集合C中，寻找最优的响应关系，用一个二进制赋值矩阵描述，其一般满足双向约束，即对二进制赋值矩阵X向量化，用x表示，则对于x，其对应的影像匹配总相似度为：

score(x)取得最大值所对应的x^*即为最优二进制赋值矩阵所对应的向量；

公式(4)可用张量积表示为：

超图匹配利用张量幂迭代算法进行求解，获取特征点点簇中每个特征点所对应的同名点。

作为优选，本发明所采用的步骤3)中基于归一化互信息测度的匹配的具体实现方式是：

在每层金字塔影像上以特征点的初始点位为中心，开辟一个矩形窗口，生成搜索窗口影像；利用所建立的核线几何约束方程，对搜索窗口影像的每个像素遍历，首先计算搜索窗口影像中当前遍历的像素到核线方程的距离，当该像素到核线方程的距离小于给定的阈值时，以该像素为中心开辟匹配窗口；当匹配窗口的几何变形超过指定阈值时，对匹配窗口中的影像进行几何粗纠正，所述匹配窗口的几何变形包括旋转角度变形参数以及缩放比例变形参数；然后，利用归一化互信息测度，计算待匹配特征点与搜索窗口中各像素的相似度；对搜索窗口各像素的归一化互信息测度值，采用非极大抑制算法提取局部最大值，归一化互信息测度局部最大值大于指定阈值所对应的像素即为特征点所对应的匹配候选点；当特征点所对应的匹配候选点个数大于5个时，对特征点所对应的匹配候选点按归一化互信息测度值进行大小排序，取前5个；最后，将特征点所对应的匹配候选点的结果按照特征点点簇为单元进行数据组织。

作为优选，本发明所提供的异源影像匹配前的数据还包括根据需要计算异源遥感影像的RPC参数。

作为优选，本发明所提供的异源遥感影像的RPC参数的获取方式是：

采用与地形无关的方法利用影像的严格几何模型，基于不同高程面生成的密集且均匀分布的虚拟控制格网，利用虚拟控制点按最小二乘原理进行解算，得到影像的RPC参数。

与现有技术相比，本发明具有以下显著优点和效果：

本发明提供了一种顾及高阶结构特征的异源遥感影像同名点匹配方法，首先进行异源影像匹配前的数据准备，根据需要计算异源遥感影像的RPC(有理多项式系数，Rational Polynomial Coefficient,RPC)参数并生成金字塔影像；然后，进行特征点簇的划分与提取，并在每个金字塔影像层，利用影像RPC参数预测待匹配点初始点位、建立近似核线几何、几何粗纠正匹配窗口影像，进行基于归一化互信息测度的匹配，寻找特征点的候选匹配点，再以特征点点簇为单元，进行顾及高阶结构特征的同名点匹配，通过基于超图的整体影像匹配方法，以确定每个特征点对应的最优同名点。对每层金字塔影像匹配结果，利用融合选权迭代法的RFM模型区域网平差方法，实现影像RPC参数的精化、错误匹配点的自动剔除以及匹配点对应物方坐标的计算。逐层精化匹配结果直到原始影像层，最终实现异源遥感影像同名点的自动可靠匹配。本发明在借鉴现有技术的基础上，对特征点提取策略进行改进，在影像划分出若干均匀分布的子区，子区中提取特征点点簇，综合运用有理函数模型和超图影像匹配模型，将几何约束特征和高阶结构特征引入到逐层金字塔影像匹配中，对特征点点簇进行整体影像匹配，可同时顾及匹配点的灰度相似性以及匹配点之间的高阶结构特征，针对影像匹配的特征，仅对特征点点簇中特征点超边进行采样，对匹配候选点超边全部考虑，既保证了超图模型的稀疏性，也使超边采样带来的信息损失降到最小；本发明有效删除匹配中错误的匹配点，提高了异源遥感影像匹配的可靠性和成功率，有效降低同名点人工量测的工作量。

附图说明

图1是本发明所提供的顾及高阶结构特征的异源遥感影像同名点匹配方法的流程示意图。

具体实施方式

下面结合附图对本发明所提供的技术方案做进一步详细描述：

实施例1：

参见图1，本发明提供了一种顾及高阶结构特征的异源遥感影像同名点匹配方法，各步骤详细阐述如下：

步骤1)异源影像匹配前的数据准备：

异源影像匹配前的数据准备，主要包括：

1)计算异源遥感影像的RPC参数：此步骤为可选步骤，仅当遥感影像只提供严格几何模型时进行。采用与地形无关的方法(可参考张永生、巩丹超等著的《高分辨率遥感卫星应用——成像模型、处理算法及应用技术》)，利用影像的严格几何模型，基于不同高程面生成的密集且均匀分布的虚拟控制格网，利用虚拟控制点按最小二乘原理进行解算，得到影像的RPC参数。

2)金字塔影像生成：采用实用、简单的3×3像元平均法，将异源遥感影像生成3级金字塔影像，原始影像层为第0层金字塔影像。

步骤2)特征点簇的划分与提取：

为了确保影像特征点均匀分布性，以及可用于后继的顾及匹配点高阶结构特征的整体影像匹配，本实施例首先对原始影像层中整幅原始影像均匀分布4×4或5×5个一定大小的子区，在各子区内部，进一步划分为7×7或10×10个均匀分布的格网。在每个格网范围内，采用摄影测量中常用的Forstner特征提取方法(可参考张祖勋、张剑清编著的《数字摄影测量学》)在每个格网内提取一个最佳特征点。如果某格网内的特征信息不明显，则将网格中心点作为特征点，从而形成特征点点簇，每一个子区为一个独立的特征点点簇单元。同时，将特征点点簇存储为文件形式，供后继影像匹配使用。

步骤3)特征点匹配候选点的确定：

特征点匹配候选点的确定，在每层金字塔影像上进行，主要包括特征点初始点位的预测、建立核线几何约束方程、匹配窗口影像的几何粗纠正以及基于归一化互信息测度的匹配。其中，特征点初始点位的预测、建立核线几何约束方程、匹配窗口影像的几何粗纠正具体实现方法，可参考专利号是ZL201110091756.8的发明专利(基于RFM模型的多源星载SAR影像自动匹配方法)。

对于基于归一化互信息测度的匹配，在每层金字塔影像上以特征点的初始点位为中心，开辟一个矩形窗口，生成搜索窗口影像；利用所建立的核线几何约束方程，对搜索窗口影像的每个像素遍历，首先计算该像素到核线方程的距离，当该像素到核线方程的距离小于给定的阈值时，以该像素为中心开辟匹配窗口。当匹配窗口的几何变形超过指定阈值时，进行匹配窗口影像的几何粗纠正，匹配窗口的几何变形包括旋转角度变形参数以及缩放比例变形参数。然后，利用归一化互信息测度，计算待匹配特征点与搜索窗口中各像素的相似度。对搜索窗口各像素的归一化互信息测度值，采用非极大抑制算法提取局部最大值，归一化互信息测度局部最大值大于指定阈值所对应的像素即为匹配候选点。当匹配候选点个数大于5个时，对匹配候选点按归一化互信息测度值进行大小排序，取前5个。最后，将特征点匹配候选点的结果按照特征点点簇为单元进行数据组织。

步骤4)顾及高阶结构特征的同名点匹配：

顾及高阶结构特征的同名点匹配，在每层金字塔影像上进行，特征点匹配候选点的确定后，以特征点点簇为单位逐一进行。具体步骤如下：

1)超图模型的构建：设某个特征点点簇为P¹，特征点个数为N₁，所对应的特征点匹配候选点集合为P²，包含特征点个数为N₂，表示集合P¹和P²的第i个特征点。集合P¹和P²对应的超图模型为G＝(V,E,A)和G′＝(V′,E′,A′)。其中，V和V′为对应的特征点顶点集合，即为超边的集合，d为超边包含特征点个数。A和A′为超边所对应的属性集。设特征点点簇P¹对应的候选同名点顶点响应集合为C的一个k元组如公式(1)所示：

c_s1＝(v₁,v′₁)，…，c_s1＝(v_k,v′_k) (1)

其中，(v₁,…,v_k)∈V，(v′₁,…,v′_k)∈V′。分别令超边e_1,…,k表示包含了顶点v₁,…,v_k的边，超边e′_1,…,k表示包含了顶点v′₁,…,v′_k的边。对于一个k阶图的超图匹配，可以通过比较两个度为k的超边的相似性来衡量k元组的相似性。指定k维相似性度量函数f_k，它的参数为超图属性集合A和A′中的属性向量。可以将更低阶的超边同时考虑，此时k阶超边的相似度张量为：

T_si⁽¹⁾＝f₁(a_i,a′_i) (2)

式中，参数γ^(k)表示k阶超边相似度的加权系数，T的右上标(k)表示张量的维数。a_1,…,k,a′_1,…,k分别表示k阶超图属性集合中的向量参数。设超边最大阶数为δ，相似性张量T^δ为高阶张量，包含了所有不同度超边的相似性信息。

本实施例中，根据异源遥感影像匹配的特点，仅选取1阶和3阶超边构建超图相似度张量。1阶超边所采用的相似性度量函数为归一化互信息测度，即考虑特征点匹配的灰度相似性。所采用的3阶超边相似性度量函数，主要顾及特征的几何相似不变性，计算公式(3)如下：

其中，ε为常数项，分别为3阶超边所对应三角形的内角。

2)超边的采样：本实施例中，为了降低计算复杂度，同时保留节点间主要连接关系，保证超图匹配的正确匹配率，具体步骤如下：对特征点点簇及特征点匹配候选点，采用k-d树建立起空间索引，对特征点点簇的每个特征点，采用随机采样的方式，随机抽取一定数量的包含该特征点的3元组作为其超边，本实施例中取50个，即每个特征点采样的超边数为50个，在特征点匹配候选点所在影像上，并非采样所有的三元组，而是仅考虑由该特征点采样的超边中，所包含特征点对应匹配候选点构成的3元组，且3元组中任意两个点均不重复。该方法既有利于超图模型稀疏性，同时利用影像匹配的特点，将所有相关联的超边均纳入考虑，避免对特征点匹配候选点构成超边进行随机采样所带来的信息损失。

3)超图模型的解算：超图匹配，就是在顶点响应集合C中，寻找最优的响应关系，可用一个二进制赋值矩阵描述，其一般满足双向约束，即对二进制赋值矩阵X向量化，用x表示，则对于x，其对应的影像匹配总相似度为：

为此，score(x)取得最大值所对应的x^*即为最优二进制赋值矩阵所对应的向量。

公式(4)可用张量积表示为：

为此，超图匹配问题可利用张量幂迭代算法进行求解，获取特征点点簇中每个特征点所对应的同名点。

对每个特征点点簇逐一处理，完成当前金字塔影像层特征点的匹配。

步骤5)基于RFM模型区域网平差删除错误匹配点：

本实施例采用像方仿射变换模型进行卫星影像区域网平差，辅之于选权迭代法，在平差过程中，利用合理的调整各观测值的权值，保证错误匹配点不影响平差结果，并实现其自动检测和定位，具体方法请参见专利号是ZL201110091756.8的发明专利(基于RFM模型的多源星载SAR影像自动匹配方法(SAR，合成孔径雷达，Synthetic Aperture Radar))。

利用上层匹配结果约束下层金字塔影像匹配，重复步骤3)、步骤4)、步骤5)，判断是否完成原始影像层影像匹配，直至输出匹配结果。