CN104077782A - 一种星载遥感图像匹配方法 - Google Patents

Abstract

本发明公开了一种星载遥感图像匹配方法，该方法包括以下步骤：获取当前星载遥感图像与历史星载遥感图像的重叠区域，作为输入图像和参考图像；提取输入图像和参考图像的显著目标区域；在显著目标区域中提取显著特征；为输入图像显著目标区域中每个尺度不变角点特征在参考图像上寻找一匹配点，得到候选匹配点对集合；剔除其中的误匹配点对；去除其中的外点并根据匹配点对计算得到星载遥感图像变换参数。本发明基于目标显著区域内的特征进行匹配，过滤了大量冗余信息，提高了匹配性能，大大减少了计算量；利用渐进最小中位方差法去除外点，提高了对复杂图像匹配的鲁棒性。本发明可广泛应用于星载遥感图像配准、目标监测等诸多应用系统中。

Description

一种星载遥感图像匹配方法

技术领域

本发明涉及图像处理和图像配准等技术领域，特别是一种星载遥感图像的匹配方法。

背景技术

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

目前的图像匹配算法大都是针对普通计算机和地面处理系统设计的。然而，由于卫星计算资源及网络带宽的限制，很多紧急任务无法在卫星上实时完成；等数据下传到地面后，由于多源遥感数据的海量性，现有地面处理系统很难及时的将有用、紧急数据推送给相关部门或系统，从而导致大量数据长期积存在数据库里。这种现状不仅远远不能满足实际应用尤其是紧急任务的需求，还大大浪费了存储资源和网络资源。针对星载系统在线应用的特殊需求，设计星载图像在线匹配算法是解决上述问题的关键，但目前还没有有效的相关算法。

星载图像在线匹配与传统的图像匹配原理基本相同，区别主要在于计算资源的有效性、计算效率的高效性和算法的自动性，即星载图像在线匹配需要在比普通计算机更少的计算资源上以更高的效率、以无需人工干预的自动处理方式完成各种复杂图像的可靠匹配。

与传统的图像匹配相比，星载图像在线匹配有一些可以利用的先验信息，如图像的分辨率差异等。另外，有一些约束也可以用以减少星载图像在线匹配的复杂度和计算量。

发明内容

本发明的目的是针对星载图像在线匹配的特殊性，提供一种有效的星载遥感图像匹配方法。

为了实现上述目的，本发明提出一种星载遥感图像匹配方法，该方法包括如下步骤：

步骤S1，根据位置信息获取当前星载遥感图像与历史星载遥感图像的重叠区域，从当前星载遥感图像和历史星载遥感图像中取出的重叠区域图像分别作为输入图像和参考图像；

步骤S2，分别提取所述输入图像和参考图像的显著目标区域；

步骤S3，在所述输入图像和参考图像的显著目标区域中分别提取显著特征；

步骤S4，根据最近邻匹配原则为所述输入图像显著目标区域中的每一个尺度不变角点特征在所述参考图像上寻找一个匹配点，得到候选匹配点对集合；

步骤S5，剔除所述候选匹配点对集合中的误匹配点对；

步骤S6：去除所述候选匹配点对集合中的外点并根据匹配点对计算得到星载遥感图像变换参数。

本发明所述方法对于星载图像处理系统具有重要的意义，其主要优点如下：

本发明只在目标显著区域内进行特征提取、匹配，过滤了大量冗余信息，提高了匹配性能，大大减少了计算量，适合星载遥感图像在线匹配；

利用散列表结构进行特征匹配，加快了匹配速度，减少了内存要求；

在利用匹配点对于显著区域的一致对应关系剔除误匹配的基础上，利用渐进最小中位方差法去除外点，提高了对复杂图像匹配的鲁棒性、稳定性和自动性。

得益于上述优点，本发明极大地提高了现有星载遥感图像匹配的精度，可广泛应用于星载遥感图像目标识别、目标监测等系统中。

附图说明

图1是本发明星载遥感图像匹配方法的流程图。

图2是根据本发明一实施例的尺度不变角点特征选择的示意图。

图3是根据本发明一实施例的角点特征亚像素定位的示意图。

图4是根据本发明一实施例的确定角点特征主方向的示例图。

图5是根据本发明一实施例构造角点特征描述向量的示意图。

具体实施方式

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

图1是本发明星载遥感图像匹配方法的流程图，如图1所示，所述方法包括以下步骤：

步骤S1，根据位置信息获取当前星载遥感图像与历史星载遥感图像的重叠区域，从当前星载遥感图像和历史星载遥感图像中取出的重叠区域图像分别作为输入图像和参考图像；

为减少计算量并提高图像匹配的正确率，首先根据当前星载遥感图像和历史星载遥感图像的位置信息，比如经纬度信息获取两幅图像的大致重叠区域，具体过程如下：取出当前星载遥感图像和历史星载遥感图像的左上角和右下角的经纬度的交集即公共经纬度，然后根据公共经纬度推导出当前星载遥感图像和历史星载遥感图像的重叠区域的像素坐标。从当前星载遥感图像和历史星载遥感图像中取出的重叠区域图像分别作为输入图像和参考图像，后续的显著目标区域及特征提取、特征匹配都基于输入图像和参考图像进行。

步骤S2，分别提取所述输入图像和参考图像的显著目标区域；

该步骤中，以图像灰度的离散余弦变换的符号函数为目标显著性度量标准来提取显著目标区域。具体地，所述步骤S2进一步包括以下步骤：

步骤S21，对于所述输入图像或参考图像I，计算其显著特征图像其中，DCT(·)和IDCT(·)分别表示离散余弦变换和离散余弦逆变换操作；

步骤S22，基于所述显著特征图像计算得到显著特征映射图像其中，ο为Hadamard点乘运算符，g表示高斯核函数，*表示卷积运算；

步骤S23，基于所述显著特征映射图像，分别提取得到所述输入图像和参考图像的显著目标区域。

该步骤中，首先对于所述显著特征映射图像进行直方图统计，在本发明一实施例中，直方图的灰度级数取为3。将直方图中频数较多对应的像素部分认为是背景，频数较少对应的像素组成的区域认为是显著目标区域。

步骤S3，在所述输入图像和参考图像的显著目标区域中分别提取显著特征；

在本发明一实施例中，所述显著特征包括多尺度角点特征及对应的描述向量，所述步骤S3进一步包括以下步骤：

步骤S31，构建所述输入图像和参考图像的高斯金字塔；

所述输入图像或参考图像I的高斯金字塔中第k层高斯图像表示为L(x,y,kσ)＝G(x,y,kσ)*I(x,y),其中*表示卷积运算，σ为平滑参数，在本发明一实施例中，σ＝0.5。

步骤S32，在所述高斯金字塔的每层上提取所述显著目标区域内的角点特征，得到角点特征集合；

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

步骤S321，计算所述输入图像或参考图像I中每个像素点X的互相关矩阵M(X)，所述互相关矩阵用来描述该像素点与其邻域像素之间的关系，计算公式为： $M\left(X\right)={\mathrm{\Σ}}_{Y\∈N\left(X\right)}w(X-Y)(\▿I\left(Y\right)){(\▿I\left(Y\right))}^T,$ 其中，表示图像I在水平、竖直方向的梯度图像，w()表示高斯权重函数，其标准差为σ，N(X)表示以像素点X为中心的邻域，邻域的半径为3σ_k，σ_k＝kσ，k为像素点X所在的高斯金字塔的层数。

步骤S322，根据所述互相关矩阵的特征值选取角点特征，并计算该角点特征的强度，得到角点特征集合。

首先，对于所述互相关矩阵M(X)进行特征值分解，即：

M(X)＝∑_i＝1,2λ_i(X)Γ_i(X)Γ_i(X)^T，

其中，λ₁(X)≤λ₂(X),λ_i(X)为特征值，Γ_i(X)表示与λ_i(X)对应的特征向量，Γ_i(X)^T表示Γ_i(X)的转置。

若特征值比值λ₁(X)/λ₂(X)>t_m，则X为候选角点特征。

对每个像素点进行上述操作，则得到候选角点特征集合C。

其中，t_m是一阈值，t_m越大得到的候选角点特征就越少，但不同的t_m取值对后续的特征匹配影响不大，在本发明一实施例中，t_m＝0.1。计算每个候选角点特征的特征强度m(X)＝λ₁(X)+λ₂(X)。

步骤S33，从所述角点特征集合中选择尺度不变角点特征；

由于角点特征在不同尺度上的稳定性不同，在本发明一实施例中，利用多尺度的特征强度极值点来近似表征角点特征的尺度不变性，其中，所述特征强度利用自相关矩阵的特征值的和来表示。具体过程为：对于高斯金字塔每一层的任一候选角点特征X，若该候选角点特征在该层及上、下相邻两层的N个邻域内强度最大，则将该候选角点特征作为尺度不变角点特征，其中，N为自然数，在本发明一实施例中取为26。图2是根据本发明一实施例的尺度不变角点特征选择的示意图。

步骤S34，计算所述尺度不变角点特征的亚像素级别的空间位置；

该步骤中，假设尺度不变角点特征X＝(x,y,k)，其中x,y,k分别为角点特征X的x、y方向坐标及其所在高斯金字塔的层数，此时三个方向的坐标都为正整数，则可利用尺度不变角点特征X及以X为中心像素的k×k邻域中的各像素点的特征强度作为权重，并根据最小二乘法计算得到所述尺度不变角点特征的精确的亚像素坐标。

令X_i＝(x_i,y_i)为k×k邻域中的第i个像素点，则其权重表示为：

$p_i=\frac{m\left(X_i\right)}{\underset{j=1}{\overset{k*k}{\mathrm{\Σ}}}m\left(X_j\right)},$

其中，m(X_i)＝λ₁(X_i)+λ₂(X_i)表示第i个像素点的特征强度，λ₁(X_i)和λ₂(X_i)为X_i处的互相关矩阵对应的特征值。在本发明一实施例中，k＝2。

则所述尺度不变角点特征X的亚像素坐标表示为：

$\stackrel{\‾}{x}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i{x}_i,\stackrel{\‾}{y}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_i{y}_i.$

其中，n表示k×k邻域中像素点的个数。

图3是根据本发明一实施例的角点特征亚像素定位的示意图。

步骤S35，确定所述尺度不变角点特征的主方向，其中，所述主方向是指与以所述尺度不变角点特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向；

在本发明一实施例中，首先在以尺度不变角点特征(x,y)为中心、以1.5*k*σ为半径的邻域窗口内进行采样，然后计算第k层高斯图像L(x,y,kσ)在上述邻域窗口内的每一像素的梯度方向然后利用直方图统计邻域像素的梯度方向，得到一个梯度方向直方图。所述梯度方向直方图是一种关于梯度方向θ的统计图，其范围是0～360度，其中每10度一个柱，总共36个柱。该梯度方向直方图的峰值代表了该极值点处邻域梯度的主方向，即作为该角点特征的方向。图4是根据本发明一实施例的确定角点特征主方向的示例图，该图显示了采用7个柱为角点特征确定主方向的示例。在该梯度方向直方图中，当存在另一个相当于主峰值80％能量的峰值时，则将这个方向认为是该角点特征的辅方向。

步骤S36，为所述尺度不变角点特征构造描述向量，其中，所述描述向量是用于刻画该尺度不变角点特征周围图像块统计特征的、由梯度方向直方图构成的向量。

在本发明一实施例中，首先将尺度不变角点特征周围图像块的坐标轴旋转为所述尺度不变角点特征的方向，以确保旋转不变性；然后将尺度不变角点特征周围16×16像素的窗口内分成4×4像素的小块，计算每个小块的8个方向的梯度方向直方图，将每个小块的梯度方向直方图连接起来形成128维的特征向量；最后，将所述特征向量归一化到单位长度即得到所述尺度不变角点特征的描述向量。当然，所述窗口以及分得的小块窗口以也可以为其他像素大小，只是最后得到的特征向量的维数有所不同。图5是根据本发明一实施例的为角点特征构造描述向量的示意图。

步骤S4，根据最近邻匹配原则为所述输入图像显著目标区域中的每一个尺度不变角点特征在所述参考图像上寻找一个匹配点，得到候选匹配点对集合；

令I_a和I_b分别表示所述输入图像和参考图像显著目标区域的尺度不变角点特征集合，集合中元素的个数分别为N_a和N_b，对应的描述向量集合为A和B。以集合B作为参考集合，对集合A中的每一个特征向量，通过查询其在参考集合B中的最近邻点和次近邻点来得到参考集合B中相应的匹配点。具体过程如下：

步骤S41，计算集合B在每一维上的数据分布；

该步骤中，首先对于集合B的每一维数据进行快速排序；然后根据快速排序结果，选择m-1个分界值把集合B中的每一维数据划分为m等份。在本发明一实施例中，m＝5。

步骤S42，建立索引结构和散列桶，以通过索引结构将集合A和集合B中相似的特征向量快速投影到同一个散列桶中；

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

步骤S421，确定索引结构的关键维和关键值；

分别为每个索引结构I_m(m＝1,…,k)随机选择n维作为I_m的关键维。在本发明一实施例中，k＝10。对于每一个关键维，从对应的m-1个分界值中随机选定一个值作为该关键维的关键值。在本发明一实施例中，n＝8。

步骤S422，计算集合B中各特征向量点的散列值；

对于集合B中的每一个特征向量点p和索引结构I_m(m＝1,…,k)，如果p在某个随机维上的值大于I_m所对应的关键值，则分配一个值1给特征向量点p，否则分配0。这样，对应于索引结构I_m，p将得到一个长度为L的二进制串，转换该二进制串得到的整数key作为所述特征向量点p的散列值。

步骤S423，将集合B中的数据点投影到散列桶中去。

该步骤中，建立一个线性投影函数f(key)＝key，使得具有相同散列值的特征点被划分到同一个散列桶中，key为散列桶的编号。

步骤S43，对集合A中的每个尺度不变角点特征，根据索引结构进行投影得到散列桶的编号，从集合B中挑选具有相同散列桶编号的尺度不变角点特征子集并在其中查找集合A中尺度不变角点特征的最近邻特征点和次近邻特征点；

对于集合A中的每一个尺度不变角点特征q，在k个索引结构{I₁,I₂,…,I_k}上进行与步骤S422和步骤S423相同的散列、投影操作并得到q所在的散列桶编号，然后只计算q与集合B中具有相同散列桶编号的尺度不变角点特征对应的描述向量的欧式距离。与最小欧式距离对应的集合B中的尺度不变角点特征为最近邻特征点，与第二最小欧式距离对应的集合B中的尺度不变角点特征为次近邻特征点。当集合B中与q具有相同散列桶编号的尺度不变角点特征个数大于T时，从集合B中随机选择T个尺度不变角点特征，只计算q与T个尺度不变角点特征描述向量之间的欧式距离。在本发明一实施例中，T＝200。

步骤S44，在集合B中为集合A中的每一个尺度不变角点特征q搜寻得到候选匹配点对。

根据所述步骤S41-S43可以在集合B中为集合A中的每一个尺度不变角点特征q搜寻得到最近邻特征点p₁和次近邻特征点p₂。设特征点q与最近邻特征点p₁和次近邻特征点p₂之间的描述向量对应的欧式距离分别为d₁和d₂,若d₁/d₂<τ，则(q,p₁)为一候选匹配点对。其中，τ为一小于1的阈值，在本发明一实施例中，τ＝0.8。

步骤S5，剔除所述候选匹配点对集合中的误匹配点对；

设步骤4得到的候选匹配点对集合表示为S＝{(a_i,b_i)|a_i＝(x_i,y_i),b_i＝(X_i,Y_i)}，i＝1,…,m，其中，m为候选匹配点对集合中元素的个数。上述候选匹配点对主要基于特征点的外观相似性，由于遥感图像的噪声，直接基于上述初始匹配利用RANSAC剔除外点有可能导致失败，并且会增加计算时间。因此，在本发明一实施例中，利用匹配点对与显著区域的一致对应关系来进行误匹配点对的快速剔除。具体为：

对于集合A中的两个特征点a₁和a₂及其在集合B中的匹配点b₁和b₂，若a₁和a₂同属于当前图像的某目标显著区域R_A,i，则b₁和b₂也应同属于参考图像的某目标显著区域R_B,j；若a₁和a₂不属于当前图像的同一目标显著区域，则b₁和b₂也应不属于参考图像的同一目标显著区域。根据上述规则，可将初始匹配中违反上述规则的匹配点对去除。

步骤S6：去除所述候选匹配点对集合中的外点并根据匹配点对计算得到星载遥感图像变换参数。

该步骤中，对于所述步骤S5保留下来的匹配点对按照特征描述子欧式距离进行升序排序。欧式距离越小，说明正确匹配的置信度越高。在本发明一实施例中，利用渐进最小中位方差法去除所述候选匹配点对集合中的外点，即从置信度高的前k个匹配点对中随机选择三个匹配点对作为种子匹配点对(在本发明一实施例中，k＝15)。

然后根据随机选择的三个匹配点对利用最小二乘法计算得到初始仿射变换T，根据所述初始仿射变换T计算匹配点对全集中的每一匹配点对(a_i,b_i)的匹配残差r_i＝||b_i-T(a_i)||₂，T(a_i)表示特征点a_i在初始仿射变换T下的坐标为特征点a_i对应的匹配点。计算所述候选匹配点对集合中内点的数目即||r_i||₂≤ε的数目，其中，ε为内点判断阈值，在本发明一实施例中，ε取为||r_i||₂的中值。

上述随机选择种子匹配点对及内点计算过程重复s次，将具有最小残差的样本集合作为最终的匹配点对集合，根据最终的匹配点对集合利用最小二乘法可以求得最终的仿射变换。在本发明一实施例中，s＝10。

本发明基于目标显著区域内的特征进行匹配，过滤了大量冗余信息，提高了匹配性能，大大减少了计算量，适合星载遥感图像在线匹配；利用散列表结构进行特征匹配，加快了匹配速度，减少了内存要求；在利用匹配点对与显著区域的一致对应关系剔除误匹配的基础上，利用渐进最小中位方差法去除外点，提高了对复杂图像匹配的鲁棒性。本发明可以广泛应用于星载遥感图像配准、星载遥感图像目标监测等诸多应用系统中。

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

Claims (10)

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

步骤S1，根据位置信息获取当前星载遥感图像与历史星载遥感图像的重叠区域，从当前星载遥感图像和历史星载遥感图像中取出的重叠区域图像分别作为输入图像和参考图像；

步骤S2，分别提取所述输入图像和参考图像的显著目标区域；

步骤S3，在所述输入图像和参考图像的显著目标区域中分别提取显著特征；

步骤S4，根据最近邻匹配原则为所述输入图像显著目标区域中的每一个尺度不变角点特征在所述参考图像上寻找一个匹配点，得到候选匹配点对集合；

步骤S5，剔除所述候选匹配点对集合中的误匹配点对；

步骤S6：去除所述候选匹配点对集合中的外点并根据匹配点对计算得到星载遥感图像变换参数。

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

步骤S21，对于所述输入图像或参考图像I，计算其显著特征图像其中，DCT(·)和IDCT(·)分别表示离散余弦变换和离散余弦逆变换操作；

步骤S22，基于所述显著特征图像计算得到显著特征映射图像其中，ο为Hadamard点乘运算符，g表示高斯核函数，*表示卷积运算；

步骤S23，基于所述显著特征映射图像，分别提取得到所述输入图像和参考图像的显著目标区域。

3.根据权利要求2所述的方法，其特征在于，所述步骤S23中，通过对于所述显著特征映射图像进行直方图统计得到所述输入图像和参考图像的显著目标区域。

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

步骤S31，构建所述输入图像和参考图像的高斯金字塔；

步骤S32，在所述高斯金字塔的每层上提取所述显著目标区域内的角点特征，得到角点特征集合；

步骤S33，从所述角点特征集合中选择尺度不变角点特征；

步骤S34，计算所述尺度不变角点特征的亚像素级别的空间位置；

步骤S35，确定所述尺度不变角点特征的主方向，其中，所述主方向是指与以所述尺度不变角点特征为中心的邻域内的梯度方向直方图的峰值对应的梯度方向；

步骤S36，为所述尺度不变角点特征构造描述向量，其中，所述描述向量是用于刻画该尺度不变角点特征周围图像块统计特征的、由梯度方向直方图构成的向量。

5.根据权利要求4所述的方法，其特征在于，所述步骤S33中，利用多尺度的特征强度极值点来从所述角点特征集合中选择尺度不变角点特征。

6.根据权利要求1所述的方法，其特征在于，所述步骤S4中，令I_a和I_b分别表示所述输入图像和参考图像显著目标区域的尺度不变角点特征集合，集合中元素的个数分别为N_a和N_b，对应的描述向量集合为A和B，以集合B作为参考集合，对集合A中的每一个特征向量，通过查询其在参考集合B中的最近邻点和次近邻点来得到参考集合B中相应的匹配点。

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

步骤S41，计算参考集合B在每一维上的数据分布；

步骤S42，建立索引结构和散列桶，以通过索引结构将集合A和集合B中相似的特征向量快速投影到同一个散列桶中；

步骤S43，对集合A中的每个尺度不变角点特征，根据索引结构进行投影得到散列桶的编号，从集合B中挑选具有相同散列桶编号的尺度不变角点特征子集并在其中查找集合A中尺度不变角点特征的最近邻特征点和次近邻特征点；

步骤S44，在集合B中为集合A中的每一个尺度不变角点特征q搜寻得到候选匹配点对。

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

步骤S421，确定索引结构的关键维和关键值；

步骤S422，计算集合B中各特征向量点的散列值；

步骤S423，将集合B中的数据点投影到散列桶中去。

9.根据权利要求1所述的方法，其特征在于，利用匹配点对与目标显著区域的一致对应关系来进行误匹配点对的剔除。

10.根据权利要求1所述的方法，其特征在于，利用渐进最小中位方差法去除所述候选匹配点对集合中的外点。