CN103822616A - 一种图分割与地形起伏约束相结合的遥感影像匹配方法 - Google Patents

Abstract

一种图分割与地形起伏约束相结合的遥感影像匹配方法，包括对基准影像进行分块，得到基准影像块；对每个基准影像块提取特征点，计算所有特征点的局部反差，选取局部反差较大的若干特征点作为种子点；投影得到各基准影像块在待匹配影像上对应的搜索区域；计算基准影像与待匹配影像间的尺度差异因子与旋转矩阵，平衡尺度与旋转差异；引入基于全球SRTM数据，对每个基准影像块基于种子点进行分层匹配，剔除误匹配点；进行影像分割与合并，提取边缘特征与分割区域，获得密集的匹配点云。本发明的优点在于可全自动实现不同分辨率、不同视角的遥感影像匹配，对于纹理贫乏、相似纹理、地表不连续等匹配困难区域，依然能够获得可靠、密集的同名点。

Description

一种图分割与地形起伏约束相结合的遥感影像匹配方法

技术领域

本发明属于测绘科学与技术领域，涉及一种图分割与地形起伏约束相结合的遥感影像匹配方法，主要应用于多源遥感影像的匹配、数字表面模型生成等领域。

背景技术

影像匹配作为数字摄影测量与计算机视觉领域的核心技术，在近十年来一直处于国内外学者研究的前沿。影像匹配即通过一定的匹配算法在两幅或多幅影像之间识别同名点的过程，广泛应用于空中三角测量、地物表面三维恢复、城市三维重建、影像变化检测等各个环节。

目前的影像匹配方法大致可分为三类：基于区域灰度的匹配、基于特征的匹配和基于语义的匹配。基于区域灰度的影像匹配是指在两张影像上比较一定尺寸窗口内像素的灰度分布相似程度，利用相关系数作为相似测度。基于特征的影像匹配是指在两张影像上分别提取特征，然后对特征进行比较来确定同名像点。基于语义的匹配是指根据从两影像中提取的语义信息确定同名像点。

基于区域灰度的匹配是使用最为广泛的方法，它常以相关系数作为同名点的相似性测度。它的优点是算法简单而且匹配精度高，然而，当影像间旋角较大或是不同比例尺影像进行匹配时，基于区域灰度的匹配方法往往结果不甚理想。

基于语义的是一种新型的匹配算法，它理论上可以具有极高的可靠性和精度，是一种理想的匹配算法。但由于它涉及诸如计算机视觉、模式识别和人工智能等许多领域，并且依赖于这些领域中许多理论方面的突破，因此，到目前为止还没有明显进展。

基于特征的匹配方法是在基于区域灰度匹配方法的基础上提出的，主要目的是利用影像的某些特殊纹理特征，从而将单纯基于点的影相匹配提升为基于特征的影像匹配，这些特征可以是点，也可以是线和面。由于对影像进行了进一步解译，基于特征的匹配方法，相对于传统的基于灰度的匹配而言，匹配结果的鲁棒性和精度更强。

1999年D.G.Lowe初步提出SIFT(Scale Invariant Feature Transform)算子，通过高斯滤波提取出在尺度空间内稳定的特征点，并用特征向量予以描述。在2004年D.G.Lowe对算法进行了进一步总结和完善，正式提出了一种基于尺度空间的、对图像缩放、旋转甚至仿射变换保持不变性的图像局部特征描述算子—SIFT算子，即尺度不变特征变换。Silverira等人提出了结合特征匹配与灰度匹配的自动匹配算法，该算法首先利用SIFT算子来提供一组种子点，然后以这些种子点为基础使用区域增长的方式进行最小二乘匹配。该方法适用于平坦地形下同名点的匹配，但对于复杂地形如居民地、山地、丘陵等，该方法匹配效果较差。李芳芳等人提出来利用线特征与点特征相结合的遥感影像匹配方法，该方法利用线特征约束获取高精度点特征，避免了几何差异与辐射差异造成的配准误差。但该算法计算量较大，且线特征无法自动提取。Zhang等人提出了一种新型的基于多基线影像匹配的高精度DEM自动生成方法，但该方法需依赖于较为精确的姿态轨道参数获得大致准确的同名预测点。Zhu等提出一种基于自适应三角形的影像匹配理论，三角网随着匹配过程中产生的新的匹配点的加入而进行实时更新，实现从纹理丰富区域至纹理贫乏区域的传播，但该方法只考虑了兴趣点的匹配，而忽视了重要的边缘特征。Zhang和Gruen采用线特征作为匹配的基本单元，并且将匹配成功的线作为断裂线合并到最终的DEM中。然而该方法中的边缘或者线特征只是用来作为约束，并且和点匹配相互独立，未能充分利用它们对影像匹配的作用。

法国欧空局与德国Inpho公司开发的处理大重叠度影像的高性能航空数字摄影测量处理软件，该软件采用了多像匹配算法，应用于大比例尺城区的空三自动转点,大大提高了自动转点的成功率。法国地理院开发了一套处理大重叠度影像的比较成熟的高性能航空数字摄影测量处理系统——像素工厂(Pixel Factory--PF)，将DPW中的自动化处理与多影像匹配方法相结合，用于高速生产正射影像。

总的来说，现有匹配方法仍存在的一些难以克服的问题，如：①影像必须具有丰富的灰度信息（高反差、高对比度），纹理信息贫乏、相似纹理、地表不连续的区域很难进行匹配；②影像间必须具有相似的比例尺；③影像间的旋转角不宜过大；④影像间的几何与辐射畸变不能太大；⑤影像的阴影或遮蔽范围不能太多等等。

发明内容

针对现有技术存在的不足，本发明提出一种图分割与地形起伏约束相结合的遥感影像匹配方法。

本发明提出了一种图分割与地形起伏约束相结合的遥感影像匹配方法，包括以下步骤：

步骤1，对基准影像进行分块，得到基准影像块；

步骤2，对每个基准影像块提取特征点，计算所有特征点的局部反差，选取局部反差较大的若干特征点作为种子点；

步骤3，投影得到各基准影像块在待匹配影像上对应的搜索区域；

步骤4，计算基准影像与待匹配影像间的尺度差异因子与旋转矩阵；

步骤5，根据尺度差异因子与旋转矩阵，对基准影像与待匹配影像平衡尺度与旋转差异；

步骤6，对每个基准影像块基于种子点进行匹配，匹配方式如下，

步骤6.1，根据基准影像块在待匹配影像上对应的搜索区域，进行分层金字塔构建；

步骤6.2，判断当前是否已处理分层金字塔的顶层，若否则基于全球SRTM数据进行顶层的近似核线生成，若是则以还未处理的分层中最高层为当前层，基于上一层获得的表面高程模型进行当前层的近似核线生成；然后，在近似核线约束下，计算局部畸变参数，进行局部几何和辐射畸变的校正；

步骤6.3，进行种子点匹配，前方交会，然后判断是否处理完原始层，若是则对该基准影像块匹配完成，否则根据当前层的匹配结果，更新表面高程模型，然后返回步骤6.2根据表面高程模型更新下一层近似核线；

步骤7，误匹配点的剔除；

步骤8、影像分割与合并，包括对各基准影像块与对应的搜索区域进行影像分割，将邻近高相关性的区域进行合并，提取边缘特征与分割区域；

步骤9、根据边缘特征与分割区域进行匹配生长，获得密集的匹配点云。

而且，步骤6.2中，进行近似核线生成采用投影轨迹法。

而且，步骤7中，误匹配点的剔除实现方式如下，

当影像摄影光线存在交会角时，则采用前方交会得到的物方三维坐标反投到基准影像与待匹配影像上，计算出像方偏差的标准误差，大于三倍标准误差的点位予以剔除；

当影像间光线平行时，则利用RANSAC算法，以二次多项式作为模型进行迭代解算，最终将误匹配点进行剔除。

而且，利用基于标记的分水岭算法对基准影像块与对应的搜索区域进行影像分割。

而且，步骤9的实现方式为，将每一分割区域看作局部平滑区域，分割区域内有种子点的，根据区域内种子点进行以下处理，分割区域内不存在种子点的，根据邻近分割区域内种子点进行以下处理，

以种子点作为根节点，确定分割区域的视差值，获得特征点和边缘点的候选同名点位置，以候选同名点为圆心，确定圆形搜索窗口；利用距离、角度及归一化相关系数作为相似性测度，在圆形搜索窗口中确定最优同名点，并通过二次曲面拟合获得子像素级点位。

而且，设左影像为基准影像，右影像为待匹配影像，圆形搜索窗口的半径r按以下方式确定，

设视为根节点的种子点在左影像上记为R、在右影像上记为R’，另外两个与根节点相邻的种子点视为起始节点，在左影像上记为P₁、P₂，在右影像上记为P₁’、P₂’，且左影像上起始节点与根节点形成的夹角记为θ，右影像上起始节点与根节点形成的夹角记为θ′，θ和θ′大于0°小于180°，左影像上P₂与R的距离为λ，左影像上P₂’与R’的距离为λ′，则r＝max(|λ-λ′|,|L·tan(θ-θ′)|)

而且，在圆形搜索窗口中确定最优同名点按以下方式实现，

设左影像上某特征点或边缘点记为点T，右影像上点T的候选同名点为点T’，利用距离、角度及归一化相关系数作为相似性测度，得到右影像上圆形搜索窗口内各像素点(x,y)的特征向量与点T的特征向量；

计算圆形搜索窗口内各像素点(x,y)的特征向量与点T的特征向量的欧式距离P(x,y)，取圆形搜索窗口内欧式距离最小的像素点为点T的同名点，

$P(x,y)=\min (\mathrm{xqrt}{\left(\right(R_r-R_l)}^2+{(A_r-A_r)}^2+{\left(\frac{1}{{\mathrm{NCC}}_r-{\mathrm{NCC}}_{\mathrm{Threshold}}}\right)}^2))$

其中，其中R_r与R_l分别为点T与像素点(x,y)的距离，A_r与A_l分别为点T与像素点(x,y)的夹角，NCC_r为点T与像素点(x,y)的归一化相关系数，NCC_Threshold为预设的相关系数阈值。

本发明公开了一种图分割与地形起伏约束相结合的遥感影像匹配方法。该方法引入全球SRTM（Shuttle Radar Topography Mission，航天飞机雷达地形测绘使命）数据作为辅助，首先消除影像间的尺度与旋转差异，并综合利用金字塔影像匹配策略、小面元几何纠正，近似核线约束，进行动态多尺度的种子点匹配，获取稳定可靠的种子点。之后基于影像分割算法，以空间位置、局部纹理反差、灰度相关作为相似性测度，进行特征点与边缘点的匹配生长，最终获得密集的匹配点云。本发明的优点在于可全自动实现不同分辨率、不同视角的遥感影像匹配，对于纹理贫乏、相似纹理、地表不连续等匹配困难区域，依然能够获得可靠、密集的同名点。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的同名点位坐标重心化示意图；

图3为本发明实施例的基于RAG的区域合并示意图；

图4为本发明实施例的加密点预测示意图；

图5为本发明实施例的搜索区域的确定示意图；

图6为本发明实施例的搜索区域内各像素点相似性测度计算示意图；

图7为本发明实施例的搜索区域内各像素点的特征向量示意图；

图8为本发明实施例的同名点的确定示意图。

具体实施方式

本发明技术方案可采用计算机软件实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。

如图1，设左影像为基准影像，右影像为待匹配影像，一般以最接近于垂直摄影的影像作为基准影像，将另外一张影像作为待匹配影像。实施例的流程如下：

步骤1、影像物理分块：针对基准影像，自动检测计算机最大读取影像范围，智能确定影像分块尺寸，得到相互独立的基准影像块，并记录影像块与原始影像的对应关系。

由于影像数据量过大，无法将其加载到内存中进行整体处理。需要通过分块方式将基准影像分割为独立且具有一定重叠关系的基准影像块，以影像块为单位，进行后续处理。在对基准影像进行分块处理时，可通过获取计算机内存大小来自动检测最大读取影像范围，智能确定影像分块尺寸，得到相互独立的基准影像块，可通过建立影像分块索引机制，并记录影像块与原始影像的对应关系，考虑到分块影像边界连续问题，块与块之间在列和行方向均可具有一定重叠。

步骤2、种子点选择：针对基准影像，利用特征提取算子，对基准影像块提取特征点。引入局部纹理反差算子，计算所有特征点的局部反差，选取局部反差较大的若干特征点作为种子点。

实施例利用特征提取算子在每个基准影像块上分别提取特征点作为兴趣点，然后引入局部纹理反差算子LC值，从兴趣点中筛选出局部反差较大的作为种子点，具体实现可参考现有技术，本发明不予赘述。

进行步骤2后，基于基准影像和待匹配影像，进行种子点匹配。首先进入步骤3。

步骤3、搜索区域的确定：

当基准影像被分为一系列块区域后，每一块需要在待匹配影像上找到覆盖同一地理范围的相应区域，在全球SRTM数据的辅助下，通过像地正反算，将基准影像块的各角点投影到待匹配影像上，得到基准影像块在待匹配影像上对应的搜索区域。即通过初始的相关性预测，某块区域的各角点可以被投影到待匹配影像上，一个搜索区域就形成了。

步骤4、初始尺度与旋转差异计算：

为消除基准影像与待匹配影像间的尺度与旋转差异，首先进行本步骤，包括以下两个字步骤：

（1）初始同名点预测。

在对遥感影像进行匹配时，影像间不可避免的会存在由空间分辨率的不同带来的尺度差异，以及由于拍摄条件等的不同带来的旋转差异。这些差异往往会造成匹配出现异常甚至失败，因此在做影像匹配前，需要先消除影像间的尺度与旋转差异。为了消除这些差异，需在两张待匹配的影像上选取初始同名点，以便计算尺度因子与旋转矩阵。在确定初始同名点时，采用的方法是利用影像姿轨文件或者POS数据或RPC文件进行影像空间坐标与地面坐标的正反算。

初始同名点预测主要包括两个过程，像地正向投影和地像反向投影。在像地正向投影中，首先选取接近于垂直摄影的影像作为基准影像，基准影像上的兴趣点通过摄影光线投影到由全球SRTM数据所内插出的真实地面高程面上，进而投影到地面从而获得对应的物方空间坐标。在地像反向投影中，则是利用共线方程多次迭代得到待匹配影像上同名点的像方坐标。

（2）根据得到的初始同名点，计算基准影像与待匹配影像的尺度因子与旋转矩阵，并根据这些参数对两影像进行重采样以消除尺度与旋转差异。首先将得到基准影像与待匹配影像的同名点的像方坐标进行重心化处理，进而得到两张影像的尺度关系等式与旋转关系式，通过解算这些方程，可以得到两张影像的尺度差异因子与旋转矩阵。

为便于实施参考起见，提供实施例的具体实现说明如下：

在初始同名点预测完成后，在左右待匹配影像上可以得到两组同名点集X₁和X₂。将各点像方坐标重心化后可得两组向量X_1G，X_2G，如图2所示，它们之间的关系可见公式1：

x_2Ｇ＝T+λRx_1G  (1)

式中,T是平移向量，λ是尺度参数，R是旋转矩阵。根据公式1，两组重心化后的坐标可以消除平移参数，即可得到公式2，根据均值欧式几何距离，λ可由公式3计算得出。

X_2G＝λRX_1G  (2)

$\mathrm{\λ}\sqrt{\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|X_{2\mathrm{Gi}}{|}^2}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\left|X_{1\mathrm{Gi}}\right|}^2}}$

式中，N为每组点集同名点的个数，i＝1,2,…,N，X_1Gi为第一组同名点集中第i点的重心化后坐标，X_2Gi为第二组同名点集中第i点的重心化后坐标。为了计算旋转矩阵R，矩阵S被构建并进行奇异值分解（SVD，Singular Value Decomposition)。在公式4中Σ是实数对角阵，并且Σ^T=Σ，左右奇异向量矩阵U、V均为标准正交矩阵。因此可以得到公式5、公式6。

$S=X_{1G}{X_{2G}}^T\⇒S=U{\mathrm{\ΣV}}^T\⇒\begin{array}{c}{U}^T\mathrm{SV}=\mathrm{\Σ}\\ V^T{S}^{T}U={\mathrm{\Σ}}^T\end{array}---\left(4\right)$

$U^T\mathrm{SV}=V^T{S}^{T}U\⇒{\mathrm{\λU}}^T{X}_{1G}{X_{1G}}^T{T}^{T}V={\mathrm{\λV}}^{T}R{X}_{1G}{X_{1G}}^{T}U---\left(5\right)$

X_1GX_1G ^T＝(UV^TR)X_1GX_1G ^T(UV^TR)  (6)

公式6成立的必要条件是UV^TR＝(±1)I，这里I是一个单位阵。因为R的行列式为1，并且U和V的行列式也为1，所以可以得到公式7，其中det（）表示求取行列式。

$R=<\begin{array}{}\end{array}\left.\begin{array}{c}{\mathrm{VU}}^T,\det \left({\mathrm{VU}}^T\right)=1\\ -{\mathrm{VU}}^T,\det \left(V{U}^T\right)=-1\end{array}\&DoubleRightArrow;R=V\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& \det \left({\mathrm{VU}}^T\right)\end{array}{U}^T\right.\right]---\left(7\right)$

步骤5、平衡尺度与旋转差异：

根据尺度差异因子，可将高分辨率影像按比例缩小，或将低分辨率影像按比例放大，在进行旋转差异消除时，将待匹配影像根据旋转矩阵进行重采样以消除两张影像间的旋转差异，至此，影像间存在的尺度与旋转差异都被消除。

步骤6，对每个基准影像块分别进行匹配，匹配方式如下，

步骤6.1，进行分层金字塔构建：根据基准影像块在待匹配影像上对应的搜索区域，建立多层金字塔。原始影像数据为最低层，即原始层，建立多层金字塔可参考现有技术。

步骤6.2，判断当前是否已处理分层金字塔的顶层，若否则进行顶层的近似核线生成，若是则对当前层进行近似核线生成，包括设还未处理的所有分层中最高层为第m层，以第m层为当前层，进行第m层的近似核线生成即可。

然后，在近似核线约束下，计算局部畸变参数，进行局部几何和辐射畸变的校正。

初次执行步骤6.2时，基于全球SRTM数据的近似核线生成，得到顶层的近似核线。包括在顶层，利用全球SRTM数据获得影像覆盖范围内的真实高程范围与搜索步距，建立近似核线。在近似核线约束下，计算局部畸变参数，以便进行小面元几何纠正，消除窗口变形对于匹配的影响，获得初始匹配点。

之后每次执行步骤6.2时，对当前层进行近似核线生成则是利用上一层获得的表面高程模型代替全球SRTM数据，自动更新并精化近似核线及局部畸变参数，这样逐层处理直到原始层，逐步获得更为精确的同名点位，更新初始匹配点。

在本方法中引入了全球的SRTM数据以及迭代搜索策略后，可以大大降低近似核线的预测误差，获得了更为精确地地表高程，从而为后续的影像匹配提供了可靠的地形约束。顶层的近似核线生成过程中，对种子点通过像地正算得到其在预测高程面上的平面坐标。根据这个平面坐标，可以得到在SRTM数据中（进行第m层的近似核线生成时是上一层获得的表面高程模型中）该位置的高程值Z_s，其中s表示种子点的编号。然后，再次将种子点通过正算得到其在高程面Z_s+1上的平面坐标，根据新得到的平面坐标，可以得到新的高程值Z_s+1。通过逐步迭代直至偏差小于给定的阈值。因此，通过现有投影轨迹法，基于利用SRTM数据或上一层获得的表面高程模型中邻域内的地形起伏（高程起伏），不断投影到待匹配影像上，可以构建出一条近似核线。

dz＝Z_s+1-Z_s（s＝0,1,2…)  （8）

式中，dz是指高程搜索步距。

再者进行的是局部几何与辐射畸变的校正。在步骤5对影像进行尺度与旋转差异消除后，影像内仍然存在局部畸变，这些畸变的存在可能会导致对应匹配窗口间的低相关性，甚至其对应在搜索影像上的匹配窗口会出现不规则，甚至不连续的现象。因此，本方法计算了匹配窗口间的仿射变换参数以及线性辐射变换参数，以对相关窗口与匹配窗口间的几何和辐射畸变进行补偿。

步骤6.3，进行种子点匹配，前方交会，然后判断是否处理完原始层，若是则对该基准影像块匹配完成，否则根据当前层的匹配结果，更新表面高程模型，然后返回步骤6.2，根据表面高程模型更新下一层近似核线；

步骤7，误匹配点的剔除。|

当影像摄影光线存在交会角时，则采用前方交会得到的物方三维坐标反投到基准影像与待匹配影像上，计算出像方偏差的标准误差，大于三倍标准误差的点位予以剔除。当影像间光线平行时，则利用RANSAC算法，以二次多项式作为模型进行迭代解算，最终将误匹配点进行剔除。

实施例的具体实现如下：

由于遥感影像具有初始的位置参数信息，因此经空间前方交会可以得到各同名点三维物方坐标。对于每一对同名点p_f(x_f,y_f)和p_f′(x_f′,y_f′)，其中f是同名点的编号，f=1,2,3，…F，F为同名点的对数，x_f与x_f′为对应点位列方向的影像坐标，y_f与y_f′为对应点位行方向的影像坐标。通过空间前方交会得到它们的三维物方坐标为((X_f,Y_f,Z_f)，然后物方坐标经后方交会到两张影像上，得到新的影像坐标p_fb(x_fb,y_fb)和p_fb′(x_fb′,y_fb′)，x_fb与x_fb′为后方交会后对应点位在列方向的影像坐标，y_fb与y_fb′为后方交会后对应点位在行方向的影像坐标。点p_f的像方偏差可以由式9计算得到，所有同名点的均方根误差可由式10计算得到。当某点的像方偏差Δp_f大于三倍标准误差RMSE时（式11），认为该点为误匹配，予以剔除。

$\mathrm{\&Delta;}{p}_f=\sqrt{\frac{{(x_f-x_{\mathrm{fb}})}^2+{(y_f-y_{\mathrm{fb}})}^2+{({x_f}^{\&prime;}-{x_{\mathrm{fb}}}^{\&prime;})}^2+{({y_f}^{\&prime;}-{y_{\mathrm{fb}}}^{\&prime;})}^2}{2}}---\left(9\right)$

$\mathrm{RMSE}=\sqrt{\frac{\underset{f=1}{\overset{F}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{p_f}^2}{F}}---\left(10\right)$

$p_f=\{\begin{array}{cc}\mathrm{Mismatch}& {\mathrm{\&Delta;p}}_f\text{\&GreaterEqual;3.RMSE}\\ \mathrm{Correctmatch}& {\mathrm{\&Delta;p}}_f<3.\mathrm{RMSE}\end{array}---\left(11\right)$

Mismatch为误匹配，Correctmatch为正确匹配。

另外，该方法可以用于多视同名点的误匹配检查，并且随着视角的增多，可以得到优于两视的判断效果。但值得注意的是，上述误匹配点剔除的方法适用于影像的投影光线存在交会角的情况，但当影像的投影光线平行时，利用空间前方交会的方法并不能得到物方三维坐标，导致该方法失效。

对于投影光线平行的同名点，因它们分布在左右影像上，把位于同一张影像上的点归为一个点集，这样就形成了两个点集，然后利用RANSAC算法，并且采用二次多项式作为模型去剔除异常值。RANSAC算法和二次多项式模型为现有技术，本发明不予赘述。

步骤8，影像分割与合并：采用分裂合并方法，首先利用基于标记的分水岭算法对基准影像块与对应的搜索区域进行影像分割，并利用LBP/C算子构建纹理表达模型，对初始分割区域进行描述；通过基于纹理特征的区域合并代价准则进行区域的合并，从而提取边缘特征与分割区域

步骤9，匹配生长：将每一分割区域看作局部平滑区域，以种子点作为根节点，确定所属分割区域的视差值，从而获得特征点和边缘的候选同名点位，以候选同名点为圆心，确定圆形搜索窗口。利用距离，角度，以及归一相关系数作为相似性测度，在搜索窗口中确定最优同名点，并通过二次曲面拟合获得子像素级点位。利用以上方法将不存在种子点的分割区域同样进行匹配，最终获得密集的匹配点云，以邻近分割区域内的种子点作为根节点，以邻近分割区域内的其他两个种子点为起始节点即可。

为了得到更加密集的匹配结果，本方法步骤8、9采用一种基于图分割的匹配生长算法。在这种算法下，经图像分割算法得到的每一个区域都作为一个空间约束，并且种子点可以为相应区域提供已知的差异值。通过多参数的相似性测度，实现同名点的增长，最终实现密集的匹配结果。具体过程如下：

首先是图像分割。利用基于分水岭的分割算法实现对基准影像块的初始区域的划分，利用LBP/C算子对初始分割区域构建联合直方图，将联合直方图作为初始区域的表达模型，利用基于纹理特征的区域合并代价准则，对初始区域进行合并，获得最终的分割结果，并且边缘点和分割区域都被提取出来。用LBP/C算子构建联合直方图具体实现为现有技术，可参见Ojala T,Pietikainen M,Maenpaa T.Multiresolution gray-scale and rotationinvariant texture classification with local binary patterns[J].Pattern Analysis andMachine Intelligence,IEEE Transactions on,2002,24(7):971-987.过程如图3所示，不同的区域用不同的标记表示，具有高空间与结构相似的临近区域被合并，错误的分割区域被剔除，例如图3中区域1、2、3、4进行多级合并，区域2合并入区域1，然后区域3合并入区域1。

其次是候选点预测。每一个区域都被当作是一个局部平滑空间，当种子点位于区域内时，其平移参数就被得到。当特征点或者边缘点位于区域内时，候选点就可以通过平移参数预测得到。如图4所示，（x，y）为左影像点位的影像坐标，（△x，△y）为影像坐标平移参数（即视差值），（x，，y，）为预测到右影像上的预测点位影像坐标，x，=x+△x，y，=y+△y。其中有种子点，在x和y方向上的平移参数可以通过计算得到。左影像上剩余还未匹配成功的特征点和边缘点经预测可以得到其在右影像上的对应点。

再者是搜索区域的确定。尽管待匹配影像间的全局性尺度与旋转差异已经被校正，但是因投影差以及地形起伏的原因，在预测点与真实同名点之间仍存在着局部畸变。因此，需要确定一个搜索区域以便得到最佳的同名匹配。如图5所示，区域内的一个种子点，视为根节点（左影像上R，右影像上R’），另外两个与根节点相邻的种子点（左影像上P₁、P₂，右影像上P₁’、P₂’），视为起始节点，且起始节点与根节点形成的夹角（左影像上θ，右影像上θ′）大于0°小于180°，但并不要求起始节点和种子点都位于同一区域内。左影像上P₂与R的距离为λ，左影像上P₂’与R’的距离为λ′。当起始节点确定后，根据角度θ和θ′，以及距离λ和λ′，局部差异就可以计算得出。剩余点之一兴趣点T，其在右影像上的预测点（即候选点、待选点）位为T’，以T’为中心形成一个圆形搜索区域，圆的半径r由式12计算得出，式中L是左影像根节点R到待匹配点T的距离，L，是右影像根节点R，到预测点T，的距离。点T的同名点就位于这个圆形的搜索区域中。

r＝max(|λ-λ′|,|L·tan(|θ-θ′)|)  （12）

最后就是根据约束条件进行匹配实现匹配点生长。在本方法中使用的约束条件是相似性测度，它由三个元素组成，分别是空间距离，角度，以及归一相关系数。首先对圆形搜索区域内各像素点进行相似性测度的计算（如图6），记录下各像素点到根节点R’的距离，记录下各像素点与根节点起始节点形成的夹角，以及各像素点的归一化相关系数（采用现有灰度互相关的相似性测度），这样每个像素点都形成了由距离，角度，归一化相关系数形成的特征向量（如图7）。同样的，在左影像上也形成了关于点T的由距离，角度，归一化相关系数形成的特征向量。计算圆形搜索区域内各像素点(x,y)的特征向量与点T的特征向量的欧式距离P(x,y)，如式13所示，其中R_r与R_l分别为点T与像素点(x,y)的距离向量，A_r与A_l分别为点T与像素点(x,y)的夹角向量，NCC_r为点T与像素点(x,y)的相关系数值，NCC_Threshold为预先设定的相关系数阈值。最小者即为点T的同名点（如图8）。同样的方法对剩余的边缘点和特征点进行匹配生长。

$P(x,y)=\min (\mathrm{xqrt}{\left(\right(R_r-R_l)}^2+{(A_r-A_r)}^2+{\left(\frac{1}{{\mathrm{NCC}}_r-{\mathrm{NCC}}_{\mathrm{Threshold}}}\right)}^2))---\left(13\right)$

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

Claims (7)

1.一种图分割与地形起伏约束相结合的遥感影像匹配方法，包括以下步骤：

步骤1，对基准影像进行分块，得到基准影像块；

步骤2，对每个基准影像块提取特征点，计算所有特征点的局部反差，选取局部反差较大的若干特征点作为种子点；

步骤3，投影得到各基准影像块在待匹配影像上对应的搜索区域；

步骤4，计算基准影像与待匹配影像间的尺度差异因子与旋转矩阵；

步骤5，根据尺度差异因子与旋转矩阵，对基准影像与待匹配影像平衡尺度与旋转差异；

步骤6，对每个基准影像块基于种子点进行匹配，匹配方式如下，

步骤6.1，根据基准影像块在待匹配影像上对应的搜索区域，进行分层金字塔构建；

步骤6.2，判断当前是否已处理分层金字塔的顶层，若否则基于全球SRTM数据进行顶层的近似核线生成，若是则以还未处理的分层中最高层为当前层，基于上一层获得的表面高程模型进行当前层的近似核线生成；然后，在近似核线约束下，计算局部畸变参数，进行局部几何和辐射畸变的校正；

步骤6.3，进行种子点匹配，前方交会，然后判断是否处理完原始层，若是则对该基准影像块匹配完成，否则根据当前层的匹配结果，更新表面高程模型，然后返回步骤6.2根据表面高程模型更新下一层近似核线；

步骤7，误匹配点的剔除；

步骤8、影像分割与合并，包括对各基准影像块与对应的搜索区域进行影像分割，将邻近高相关性的区域进行合并，提取边缘特征与分割区域；

步骤9、根据边缘特征与分割区域进行匹配生长，获得密集的匹配点云。

2.根据权利要求1所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：步骤6.2中，进行近似核线生成采用投影轨迹法。

3.根据权利要求1所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：步骤7中，误匹配点的剔除实现方式如下，

当影像摄影光线存在交会角时，则采用前方交会得到的物方三维坐标反投到基准影像与待匹配影像上，计算出像方偏差的标准误差，大于三倍标准误差的点位予以剔除；

当影像间光线平行时，则利用RANSAC算法，以二次多项式作为模型进行迭代解算，最终将误匹配点进行剔除。

4.根据权利要求1所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：步骤8中，利用基于标记的分水岭算法对基准影像块与对应的搜索区域进行影像分割。

5.根据权利要求1或2或3或4所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：步骤9的实现方式为，将每一分割区域看作局部平滑区域，分割区域内有种子点的，根据区域内种子点进行以下处理，分割区域内不存在种子点的，根据邻近分割区域内种子点进行以下处理，

以种子点作为根节点，确定分割区域的视差值，获得特征点和边缘点的候选同名点位置，以候选同名点为圆心，确定圆形搜索窗口；利用距离、角度及归一化相关系数作为相似性测度，在圆形搜索窗口中确定最优同名点，并通过二次曲面拟合获得子像素级点位。

6.根据权利要求5所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：设左影像为基准影像，右影像为待匹配影像，圆形搜索窗口的半径r按以下方式确定，

设视为根节点的种子点在左影像上记为R、在右影像上记为R’，另外两个与根节点相邻的种子点视为起始节点，在左影像上记为P₁、P₂，在右影像上记为P₁’、P₂’，且左影像上起始节点与根节点形成的夹角记为θ，右影像上起始节点与根节点形成的夹角记为θ′，θ和θ′大于0°小于180°，左影像上P₂与R的距离为λ，左影像上P₂’与R’的距离为λ′，则r＝max(|λ-λ′|,|L·tan(|θ-θ′)|)。

7.根据权利要求5所述图分割与地形起伏约束相结合的遥感影像匹配方法，其特征在于：在圆形搜索窗口中确定最优同名点的按以下方式实现，

设左影像上某特征点或边缘点记为点T，右影像上点T的候选同名点为点T’，利用距离、角度及归一化相关系数作为相似性测度，得到右影像上圆形搜索窗口内各像素点(x,y)的特征向量与点T的特征向量；

计算圆形搜索窗口内各像素点(x,y)的特征向量与点T的特征向量的欧式距离P(x,y)，取圆形搜索窗口内欧式距离最小的像素点为点T的同名点，

$P(x,y)=\min (\mathrm{xqrt}{\left(\right(R_r-R_l)}^2+{(A_r-A_r)}^2+{\left(\frac{1}{{\mathrm{NCC}}_r-{\mathrm{NCC}}_{\mathrm{Threshold}}}\right)}^2))$

其中，其中R_r与R_l分别为点T与像素点(x,y)的距离，A_r与A_l分别为点T与像素点(x,y)的夹角，NCC_r为点T与像素点(x,y)的归一化相关系数，NCC_Threshold为预设的相关系数阈值。

Images