CN103456022A - 一种高分辨率遥感图像特征匹配方法 - Google Patents

Abstract

本发明公开了一种高分辨率遥感图像特征匹配方法，该方法包括：在参考图像上提取显著目标区域；提取参考图像和输入图像的SIFT特征点；为参考图像的每一SIFT特征点在输入图像的SIFT特征点集合中寻找多个候选匹配点；为参考图像显著目标区域内的SIFT特征点集合在对应的所述候选匹配点组成的候选匹配点集合中搜索最优匹配点；对每个显著目标区域得到其最优匹配点集合，每个显著目标区域对应的SIFT特征最优匹配点集合的并集作为最后的特征匹配集合。本发明基于目标显著区域内的SIFT特征进行匹配，过滤了大量冗余信息，提高了匹配性能；减少了外点，提高了对视角变化、图像形变的鲁棒性。本发明可以广泛应用于图像配准、三维重建等诸多应用系统中。

Description

一种高分辨率遥感图像特征匹配方法

技术领域

本发明涉及遥感图像处理、特征匹配、图像配准等技术领域，特别是一种高分辨率遥感图像特征匹配方法。

背景技术

图像配准是将同一地点、不同时间或不同传感器拍摄的两幅图像（一幅为参考图像，另一幅为输入图像）进行几何匹配的过程。图像配准是图像融合、变化检测等应用的前提和基础，图像配准的精度对后续的应用具有重要的影响。

目前的图像配准算法可分为基于像素的配准方法和基于特征的配准方法。由于噪声、拍摄条件、季节变化、视角变化、平台抖动等因素的影响，基于特征的配准方法更适合遥感图像配准。近年来，SIFT特征（尺度不变特征变换，Scale Invariant Feature Transform）在基于特征的遥感图像配准中得到了广泛应用。特征匹配是影响配准精度的重要因素，传统的SIFT特征匹配一般基于参考图像和输入图像的SIFT特征集合利用k-d tree数据结构和最近邻匹配原则进行匹配。由于遥感图像中重复性结构多，传统的全局匹配方法易受背景影响，错误匹配率较高。

遥感图像中显著目标一般具有独特的几何结构，匹配的SIFT特征集合之间应满足一定的几何拓扑约束。从SIFT特征集合中过滤掉植被、沙漠等背景，综合利用显著目标区域的SIFT特征集合之间的外观、几何拓扑相似性是提高匹配正确性的重要途径，但目前还没有有效的相关算法。

发明内容

本发明的目的是针对遥感图像特征匹配的难点和显著目标区域SIFT特征的几何拓扑特性，提供一种有效的高分辨率遥感图像特征匹配方法。

为了实现上述目的，本发明提出一种高分辨率遥感图像特征匹配方法，该方法包括以下步骤：

步骤S1，在参考图像上提取显著目标区域；

步骤S2，提取所述参考图像和输入图像的SIFT特征点；

步骤S3，根据最近邻匹配原则为所述参考图像的每一个SIFT特征点在输入图像的SIFT特征点集合中寻找多个候选匹配点；

步骤S4，为所述参考图像显著目标区域内的SIFT特征点集合在对应的所述候选匹配点组成的候选匹配点集合中搜索最优匹配点；

步骤S5：对每个显著目标区域按照所述步骤S4得到其最优匹配点集合，每个显著目标区域对应的SIFT特征最优匹配点集合的并集作为最后的特征匹配集合。

本发明所述方法对于提高特征匹配的数目和正确匹配率具有重要的意义，其主要优点如下：

本发明基于目标显著区域内的SIFT特征进行匹配，过滤了大量冗余信息，提高了匹配性能；

利用目标显著区域内的SIFT特征集合之间的外观特征、几何拓扑特征的相似性进行区域匹配，减少了外点；

利用显著目标区域内的SIFT特征的外观特征、几何拓扑特征在输入图像中搜索SIFT特征集合，提高了对视角变化、图像形变的鲁棒性；

得益于上述优点，本发明极大地提高了现有高分辨率遥感图像特征匹配的精度，可广泛应用于图像配准、目标识别、目标跟踪等系统中。

附图说明

图1是本发明提出的高分辨率遥感图像特征匹配方法的流程图；

图2是根据本发明一实施例的遥感图像显著目标区域提取结果示意图；

图3是根据本发明一实施例的DOG金字塔构建示意图；

图4是根据本发明一实施例的SIFT特征提取中极值点提取示意图；

图5是根据本发明一实施例的SIFT特征提取中主方向检测示意图；

图6是根据本发明一实施例的为SIFT特征点构造描述向量的示意图；

图7是根据本发明一实施例的最近邻匹配结果示意图；

图8是根据本发明一实施例的区域匹配结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

图1是本发明提出的高分辨率遥感图像特征匹配方法的流程图，如图1所示，所述方法包括以下步骤：

步骤S1，在参考图像上提取显著目标区域；

该步骤中，以图像灰度的离散余弦变换的符号函数为目标区域显著性度量标准来提取所述显著目标区域；

所述步骤S1进一步包括以下步骤：

步骤S11，对于所述参考图像X，计算其显著特征图像

其中，DCT(·)和IDCT(·)分别表示离散余弦变换和离散余弦逆变换，sign(x)表示符号函数，取值为变量x的符号；

步骤S12，基于所述显著特征图像计算显著特征映射图像

其中，ο表示Hadamard点乘运算符，g表示高斯核函数，*表示卷积运算；

步骤S13，基于所述显著特征映射图像提取所述参考图像的显著目标区域。

该步骤中，首先对所述显著特征映射图像进行直方图统计，直方图的灰度级数可设为3；

然后将直方图的频数较多的部分认为是背景，频数较少的像素组成的区域为显著目标区域。

根据本发明一实施例的遥感图像显著目标区域提取结果的示意图如图2所示。

步骤S2，提取所述参考图像和输入图像的SIFT特征点；

所述步骤S2进一步包括以下步骤：

步骤S21，分别构建所述参考图像和输入图像的DOG（Difference ofGaussian，高斯差分）金字塔；

设图像为I(x,y)，则DOG金字塔中第k层高斯差分图像表示为

其中，L(x,y,σ)=G(x,y,σ)*I(x,y)为高斯平滑图像， $G(x,y,\mathrm{\σ})=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\σ}}^2}},$ *表示卷积运算。

根据本发明一实施例的DOG金字塔构建示意图如图3所示，图3中，图像在5个尺度上进行高斯卷积，得到的DOG金字塔有4个图像。

步骤S22，在构建得到的DOG金字塔的每层上提取极值点，所谓极值点是指在局部邻域内D(x,y,kσ)取值最大的点；

该步骤中，提取极值点的过程为，选择DOG金字塔上的任一点，若该点在该层及上、下相邻两层的多个（比如26个）邻域内不是极值点，则将该点去除，否则将该点作为极值点。根据本发明一实施例的提取极值点的示意图如图4所示，图4中，DOG金字塔第k层图像中标“×”的点与第k-1层、第k+1层的26个标“○“的点进行比较，若标“×”的点对应的D(x,y,kσ)是这26个邻域点中的最大值，则标“×”的点为极值点。

步骤S23，对于所提取的极值点，去除其中局部曲率非常不对称的极值点，保留下来的极值点即为SIFT特征点；

在本发明一实施例中，计算差分图像D的局部Hessian矩阵H，然后去除满足下面条件的极值点：tr(H)²/det(H)>10，其中det(H)表示矩阵H的行列式值，tr(H)表示矩阵H的迹；

步骤S24，计算SIFT特征点的亚像素级别的空间位置、尺度；

假设SIFT特征点X1的坐标参数表示为X=(x,y,σ)，其中x,y,σ分别为特征点X1的x、y方向坐标及尺度参数，此时x、y及尺度这三个方向的坐标都为正数。设SIFT特征点X1的亚像素值为X+ΔX，ΔX为待求的坐标增量，则计算SIFT特征点的亚像素级别的空间位置、尺度的步骤具体为：

首先，对高斯差分图像D(X+ΔX)在特征点X1处按泰勒展开公式展开： $D(X+\mathrm{\ΔX})=D\left(X\right)+\frac{{\∂D}^T}{\∂X}\mathrm{\ΔX}+\frac{1}{2}{\mathrm{\ΔX}}^T\frac{{\∂}^{2}D}{{\∂X}^2}\mathrm{\ΔX},$ 令 $\frac{\∂D(X+\mathrm{\ΔX})}{\∂\mathrm{\ΔX}}=0$ 则求得

其中，

和可由SIFT特征点X1及其邻域的像素点根据差分法计算得到；

然后，按照如下规则对SIFT特征点进行亚像素级插值：如果该SIFT特征点在三个方向上的坐标增量都小于0.5，则该点就是极值点，而X+ΔX为所求的亚像素级极值点的坐标；如果某一方向上的坐标增量大于等于0.5，如假设x方向的偏移量大于0.5像素，则将坐标增量四舍五入到一整数值a，将a与x相加得到新的SIFT特征点X2=(x+a,y,s)，然后将特征点X2按上述步骤进行操作，直到三个方向上的坐标增量都小于0.5。

步骤S25，根据所述SIFT特征点的亚像素级别的空间位置、尺度确定所述SIFT特征点的主方向，所谓主方向是指与以SIFT特征点为中心的邻域内的梯度方向直方图的峰值对应的梯度方向。

在本发明一实施例中，在以SIFT特征点X=(x,y,σ)为中心、以1.5σ为半径的邻域窗口内采样，计算高斯平滑图像L(x,y,σ)在上述邻域窗口内的每一像素的梯度方向并用直方图统计邻域像素的梯度方向，得到一个梯度方向直方图。所述梯度方向直方图是一种关于梯度方向θ的统计图，其范围是0～360度，其中每10度一个柱，总共36个柱。该梯度方向直方图的峰值代表了该极值点处邻域梯度的主方向，即作为该SIFT特征点的方向。根据本发明一实施例的梯度直方图示例图如图5所示，图5中显示了采用7个柱为SIFT特征点确定主方向的示例。在该梯度方向直方图中，当存在另一个相当于主峰值80%能量的峰值时，则将这个方向认为是该SIFT特征点的辅方向。一个SIFT特征点可能会被指定具有多个方向（一个主方向，一个以上辅方向）。

步骤S26，根据所述SIFT特征点的主方向，为SIFT特征点构造描述向量，所谓描述向量是指用于刻画该SIFT特征点周围图像块统计特征的、由梯度方向直方图构成的向量。

在本发明一实施例中，首先将SIFT特征点周围图像块的坐标轴旋转为所述SIFT特征点的主方向，以确保旋转不变性；然后将SIFT特征点周围kn×kn像素的窗口内分成n×n像素的小块，计算每个小块的q个（比如8个）方向的梯度方向直方图，将每个小块的梯度方向直方图连接起来形成k×k×q维的SIFT特征向量；最后，将其归一化到单位长度得到所述SIFT特征点的描述向量。图6是根据本发明一实施例的为SIFT特征点构造描述向量的示意图。

步骤S3，根据最近邻匹配原则为所述参考图像的每一个SIFT特征点在输入图像的SIFT特征点集合中寻找多个候选匹配点；

该步骤具体为：利用所述输入图像的SIFT特征点集合构造k-d树，将所述参考图像的每一个SIFT特征点与所述k-d树进行匹配，即计算所述参考图像的每一个SIFT特征点与所述k-d树中每个结点的描述向量之间的欧式距离，将距离最小的多个比如3个SIFT特征点作为所述参考图像中该SIFT特征点的候选匹配点。k-d树是一种分割高维数据空间的数据结构，k-d树通过对输入图像的SIFT特征点集合进行层次划分并建立索引来加快匹配速度。对所述参考图像某显著目标区域特征点集P_i中的每一个特征点p_k，若令C_k表示特征点p_k在输入图像上的候选匹配点集，则与特征点集P_i对应的候选匹配点集合记为

图7是根据本发明一实施例的最近邻匹配结果示意图。

步骤S4，为所述参考图像显著目标区域内的SIFT特征点集合在对应的所述候选匹配点组成的候选匹配点集合中搜索最优匹配点；

所述步骤S4具体为：

利用SIFT特征点集合之间的相似性进行所述参考图像显著目标区域的SIFT特征点集合与所述候选匹配点集合的匹配，为所述参考图像的每个显著目标区域对应的SIFT特征点集合找到对应的最优匹配点集合；

在本发明一实施例中，所述匹配主要是依据集合之间的外观相似性和几何拓扑相似性来进行，所述几何拓扑相似性主要包括几何形状、位置偏差和位置一致性。

具体地，对于匹配点对(p_k,m_k)，本发明使用

$A(p_k,m_k)=\frac{1}{1+\exp (-{\mathrm{\τ}}_a(\frac{1}{{\mathrm{\μ}}_a}{\left|\right|f_{p_k}-f_{m_k}\left|\right|}_2-1))}$ 来描述其外观相似性，其中，m_k表示参考图像上某特征点p_k的对应的在输入图像上的候选匹配点，

和

表示p_k和m_k对应的SIFT特征向量，τ_a和μ_a是调整外观相似性的参数，在本发明一实施例中，τ_a=1.5，μ_a=0.05。

对于正确的或者说是最优的两个匹配点对(p_k,m_k)和(p_k+1,m_k+1)，其几何拓扑相似性是指直线段与

与

的长度接近且趋于平行，本发明使用

G(p_k,p_k+1,m_k,m_k+1)=||(p_k-p_k+1)-(m_k-m_k+1)||₂=||(p_k-m_k)-(p_k+1-m_k+1)||₂来表示相邻点对(p_k,p_k+1)和(m_k,m_k+1)几何形状的差异。

因为参考图像与输入图像已经配准，最优匹配点对(p_k,m_k)的位置偏差不会太大，本发明使用

来描述两个点的位置偏差，其中，t为一阈值，在本发明一实施例中，t=20。

对于正确的两个匹配点对(p_k,m_k)和(p_k+1,m_k+1)，p_k与p_k+1的位置关系与m_k与m_k+1的位置关系应一致，即若p_k在p_k+1的左上方，则m_k应在m_k+1的左上方。本发明使用 $O(p_k,p_{k+1},m_k,m_{k+1})=\begin{array}{}\left\{\begin{array}{cc}1& S((p_k,p_{k+1}),(m_k,m_{k+1}))=1\\ 0& S((p_k,p_{k+1}),(m_k,m_{k+1}))=0\end{array}\right.\end{array}$ 来描述位置一致性，其中，S((p_k,p_k+1),(m_k,m_k+1))表示位置一致关系，S((p_k,p_k+1),(m_k,m_k+1))=1表示位置一致，S((p_k,p_k+1),(m_k,m_k+1))=0表示位置不一致。

设某区域按行或按列得到的特征点集P_i与候选匹配点集合C_i的最优匹配为

，则M^*可以通过优化如下目标函数求得：

$C(P,M)=\underset{k=1}{\overset{l_i-1}{\mathrm{\Σ}}}{w}_{g}G(p_k,p_{k+1},m_k,m_{k+1})+\underset{k=1}{\overset{l_i}{\mathrm{\Σ}}}{w}_{a}A(p_k,m_k)+\underset{k=1}{\overset{l_i-1}{\mathrm{\Σ}}}{w}_{o}O(p_k,p_{k+1},m_k,m_{k+1})+\underset{k=1}{\overset{l_i}{\mathrm{\Σ}}}{w}_{d}D(p_k,m_k)$ 上述目标函数可以利用动态规划策略求解，其中，G(p_k,p_k+1,m_k,m_k+1)表示相邻点对(p_k,p_k+1)和(m_k,m_k+1)几何形状的差异；A(p_k,m_k)表示匹配点对(p_k,m_k)的外观相似性；O(p_k,p_k+1,m_k,m_k+1)表示两个匹配点对(p_k,m_k)和(p_k+1,m_k+1)的位置关系一致性；D(p_k,m_k)表示最优匹配点对(p_k,m_k)的位置偏差；w_g、w_a、w_o和w_d是上述区域相似性特征的权重，在本发明一实施例中，w_g=1.0，w_a=1.25，w_o=1.5，w_d=4.0。为了得到全局最优解，要求上述目标函数对按行排列的特征集P_i ^R和按列排列的特征集P_i ^C都是最优的。

步骤S5：对每个显著目标区域按照所述步骤S4得到其最优匹配点集合，每个显著目标区域对应的SIFT特征最优匹配点集合的并集作为最后的特征匹配集合。

图8是根据本发明一实施例的区域匹配结果示意图。

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种高分辨率遥感图像特征匹配方法，其特征在于，该方法包括以下步骤：

步骤S1，在参考图像上提取显著目标区域；

步骤S2，提取所述参考图像和输入图像的SIFT特征点；

步骤S3，根据最近邻匹配原则为所述参考图像的每一个SIFT特征点在输入图像的SIFT特征点集合中寻找多个候选匹配点；

步骤S4，为所述参考图像显著目标区域内的SIFT特征点集合在对应的所述候选匹配点组成的候选匹配点集合中搜索最优匹配点；

步骤S5：对每个显著目标区域按照所述步骤S4得到其最优匹配点集合，每个显著目标区域对应的SIFT特征最优匹配点集合的并集作为最后的特征匹配集合。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1进一步包括以下步骤：

步骤S11，对于参考图像X，计算其显著特征图像其中，DCT(·)和IDCT(·)分别表示离散余弦变换和离散余弦逆变换，sign(x)表示符号函数，取值为变量x的符号；

步骤S12，基于所述显著特征图像计算显著特征映射图像

其中，o表示Hadamard点乘运算符，g表示高斯核函数，*表示卷积运算；

步骤S13，基于所述显著特征映射图像提取所述参考图像的显著目标区域。

3.根据权利要求2所述的方法，其特征在于，所述步骤S13中，首先对所述显著特征映射图像

进行直方图统计；然后将直方图的频数较多的部分认为是背景，频数较少的像素组成的区域为显著目标区域。

4.根据权利要求1所述的方法，其特征在于，所述步骤S2进一步包括以下步骤：

步骤S21，分别构建所述参考图像和输入图像的高斯差分DOG金字塔；

步骤S22，在构建得到的DOG金字塔的每层上提取极值点；

步骤S23，对于所提取的极值点，去除其中局部曲率非常不对称的极值点，保留下来的极值点即为SIFT特征点；

步骤S24，计算SIFT特征点的亚像素级别的空间位置、尺度；

步骤S25，根据所述SIFT特征点的亚像素级别的空间位置、尺度确定所述SIFT特征点的主方向，所谓主方向是指与以SIFT特征点为中心的邻域内的梯度方向直方图的峰值对应的梯度方向；

步骤S26，根据所述SIFT特征点的主方向，为SIFT特征点构造描述向量，所谓描述向量是指用于刻画该SIFT特征点周围图像块统计特征的、由梯度方向直方图构成的向量。

5.根据权利要求4所述的方法，其特征在于，所述步骤S22中提取极值点时，选择高斯差分金字塔上的任一点，若该点在该层及上、下相邻两层的多个邻域内不是极值点，则将该点去除，否则将该点作为极值点。

6.根据权利要求4所述的方法，其特征在于，所述步骤S23中，局部曲率非常不对称的极值点为满足下面条件的极值点：tr(H)²/det(H)>10，其中det(H)表示差分图像的局部Hessian矩阵H的行列式值，tr(H)表示矩阵H的迹。

7.根据权利要求4所述的方法，其特征在于，所述步骤S24进一步包括以下步骤：

首先，对高斯差分图像D(X+ΔX)在特征点X1处按泰勒展开公式展开： $D(X+\mathrm{\ΔX})=D\left(X\right)+\frac{{\∂D}^T}{\∂X}\mathrm{\ΔX}+\frac{1}{2}{\mathrm{\ΔX}}^T\frac{{\∂}^2}{{\∂X}^2}\mathrm{\ΔX},$ 令 $\frac{\∂D(X+\mathrm{\ΔX})}{\∂\mathrm{\ΔX}}=0$ 则求得

其中，和

可由SIFT特征点X1及其邻域的像素点根据差分法计算得到，X为特征点X1的坐标参数表示X＝(x，y，σ)，x，y，σ分别为特征点X1的x、y方向坐标及尺度参数，X+ΔX为特征点X1的亚像素值，ΔX为待求的坐标增量；

然后，对SIFT特征点进行亚像素级插值：如果该SIFT特征点在三个方向上的坐标增量都小于0.5，则该点就是极值点，而X+ΔX为所求的亚像素级极值点的坐标；如果某一方向上的坐标增量大于等于0.5，则将坐标增量四舍五入到一整数值a，将a与x相加得到新的SIFT特征点X2=(x+a,y,s)，然后将特征点X2按上述步骤进行操作，直到三个方向上的坐标增量都小于0.5。

8.根据权利要求4所述的方法，其特征在于，所述步骤S26中，首先将SIFT特征点周围图像块的坐标轴旋转为所述SIFT特征点的主方向，以确保旋转不变性；然后将SIFT特征点周围kn×kn像素的窗口内分成n×n像素的小块，计算每个小块的q个方向的梯度方向直方图，将每个小块的梯度方向直方图连接起来形成k×k×q维的SIFT特征向量；最后，将其归一化到单位长度得到所述SIFT特征点的描述向量。

9.根据权利要求1所述的方法，其特征在于，所述步骤S3中，利用所述输入图像的SIFT特征点集合构造k-d树，将所述参考图像的每一个SIFT特征点与所述k-d树进行匹配，即计算所述参考图像的每一个SIFT特征点与所述k-d树中每个结点的描述向量之间的欧式距离，将距离最小的多个SIFT特征点作为所述参考图像中该SIFT特征点的候选匹配点。

10.根据权利要求1所述的方法，其特征在于，所述步骤S4中，利用SIFT特征点集合之间的相似性为所述参考图像的每个显著目标区域对应的SIFT特征点集合找到对应的最优匹配点集合。

Images