CN104318548A - 一种基于空间稀疏度和sift特征提取的快速图像配准实现方法 - Google Patents

Abstract

本发明公开了一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，主要克服经典SIFT特征提取算法中，经常会提取到平滑区域和纹理区域不稳定的特征点而造成特征误匹配的问题。其步骤包括：1)分别提取参考图和待匹配图稀疏区域。2)对参考图和待匹配图的稀疏区域提取SIFT特征点。3)对参考图和待匹配图上提取到的SIFT特征点集进行粗匹配。4)用随机一致性估计算法来滤除粗匹配结果中的误匹配。5)利用参考图和待匹配图最终的匹配点对，通过仿射变换得到仿射变换参数来实现对两幅SAR图像配准。本发明能在保证经典SIFT特征配准算法精度前提下提高配准效率，可用于SAR图像的配准处理。

Description

一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法

技术领域

本发明属于图像处理技术领域，涉及图像配准，具体是一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，可用于SAR影像的变化检测、融合、拼接等方面的前期配准工作。

背景技术

合成孔径雷达(SAR)所成的图像具有全天候、全天时、高分辨率和强大的穿透能力等特点，因此，这种图像被广泛的应用到目标识别，变换检测和水面监视，日益成为了当今很具有代表性的对地观测手段之一。在对SAR图像进行拼接、融合、变化检测等操作前，需将来自同一地区、在相同时间不同视点或同一传感器不同时间获取的图像在空间上进行配准，消除因获取图像的时间、角度、环境和传感器成像机理的不同造成图像间的平移、旋转、伸缩及局部形变等问题。

作为遥感图像处理中的一个基本任务，配准的目标就是将不同时间、不同视角或者不同传感器获得的两幅或多幅图像进行几何对准的处理过程。对于SAR图像，其相干成像过程中形成的斑点噪声对图像灰度分布造成影响，另外成像条件的差异、拍摄时间和季节变化、拍摄场景地物的变化等因素造成图像具有灰度分布复杂、不稳定的特点，使得在光学遥感图像中能够成功应用的算法和技术在应用于SAR图像时却屡屡受挫。

目前，大量配准的方法已经被提出。根据特征空间的不同，配准算法一般可以分为两类：基于灰度的配准方法和基于特征的配准方法。基于灰度的方法一般都是直接利用整幅图像的灰度信息，通过建立某种像素间的相似性测度来衡量两幅图像重叠部分地表反射属性的匹配程度，进而寻找到最优匹配时的平移、旋转和伸缩等变换参数。基于灰度的图像配准方法，其相似性度量函数包括互相关函数、傅里叶函数和互信息函数等，通过使用相似性最大情况下所对应的变换关系进行变换，最终实现图像之间的配准。在利用基于灰度的配准方法对SAR图像的配准时，容易受到SAR图像成像过程中导致的灰度差异及噪声的影响。基于特征的图像配准方法是通过从图像中选择显著和突出的目标进行检测，提取这些特征点的特征，通过这些特征之间的匹配建立图像间的变换关系，最终实现图像配准。基于特征的图像配准方法能够解决存在较大几何畸变及灰度信息差别的图像间的配准问题。常用的图像特征有：特征点(包括角点、高曲率点、拐点等)、直线段、边缘、轮廓、闭合区域、特征结构等。基于特征的配准方法则首先需要对待配准图像提取特征点、边缘等，然后进行对应特征对或特征集之间的搜索匹配，进而得到图像变换参数。基于边缘特征的配准方法需利用一些边缘检测算子进行边缘特征提取，再对提取的边缘进行匹配，但往往由于提取的边缘过于琐碎，不能很好地体现图像中的结构特征，对后续匹配精度和速度造成影响。而点特征相对于线、面等特征有着易于提取、速度快、精度高、稳定性好等许多方面的优势而有着更广泛的应用。

基于特征的图像配准对噪声和形变的抗干扰能力更强，因而具有更高的鲁棒性。最常用的基于特征的配准方法之一是基于SIFT特征提取的配准方法。Lowe将SIFT特征的主要步骤概括为：(1)尺度空间极值点提取；(2)极值点定位；(3)主方向分配；(4)关键点描述。SIFT在较小旋转角度下具有较强的旋转、尺度不变性，归一化后还对光照变化鲁棒。但由于SIFT采用了一个128维的向量来描述特征点，在特征点较多的情况下增加了运算代价和时间复杂度，并且在寻找最佳匹配时需要特征间具有较好的区分度。在利用SIFT特征提取算法处理SAR图像时，不仅在直线、折线、交叉处等稀疏结构的区域提取到较多的SIFT特征点，而且在特征不明显的平滑区域和纹理区域也能提取到很多特征点，但是这些平滑区域和纹理区域的特征点周围存在很多与其相似的点，因而耗时较多且这些点特征不明显，容易形成误匹配。

发明内容

本发明的目的在于克服已有技术的不足，提出了一种联合空间稀疏结构和SIFT特征提取的SAR图像配准方法，以实现对SAR图像中稀疏结构区域集中提取SIFT特征点，在减少算法运算代价和时间效率的基础上降低了匹配错误率，提高图像配准的效率。

实现本发明目的的技术思路是：输入的参考图和待匹配图，分别获取两幅图像的稀疏区域，对参考图和待匹配图的稀疏区域提取SIFT特征点，利用最近次近距离比值准则(NNDR)作为相似度量，得到粗匹配对，并利用随机一致性估计(RANSAC)方法来去除误匹配，最后利用仿射变换来求得仿射变换参数，从而实现参考图和待匹配图的配准。

本发明的技术方案是：一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，包括如下步骤：

(1)输入在不同时间获取的同一地区的两幅多时相SAR图像，分别记作参考图和待匹配图，对参考图和待匹配图的每个像素点计算稀疏度，挑选稀疏度高的像素点所在的区域作为稀疏区域，分别记作Mask_R和Mask_S；

(2)分别对参考图稀疏区域Mask_R和待匹配图稀疏区域Mask_S提取SIFT特征点，得到参考图和待匹配图SIFT特征点集分别记为{R_m}、{S_n}；

(3)对参考图和待匹配图上提取到的SIFT特征点集{R_m}和{S_n}进行特征匹配，得到粗匹配特征对{(pq)_i}，参考图和待匹配图中的粗匹配点集分别为{p_i}和{q_i}，i表示匹配对的索引号；

(4)利用随机一致性估计算法对粗匹配结果去除误匹配对，去除误匹配对后参考图和待匹配图保留下来的特征点对集合记为{(PQ)_i}，此时参考图和待匹配图中各自保留下来特征点分别为{P_i}、{Q_i}，其中i表示匹配对的索引号；

(5)利用参考图和待匹配图的匹配点集{P_i}和{Q_i}，通过仿射变换得到仿射变换参数，最后根据仿射变换参数对两幅图像配准。

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

1)分别对参考图R和待匹配图S进行两层小波分解，小波基取db1，得到参考图和待匹配图的细节分量图，分别记为R_db1、S_db1；

2)对参考图细节分量图和待匹配图细节分量图的每一个像素点计算方差系数，然后以当前像素点为中心、取一个大小为m*m的邻域计算方差，即可得到当前像素点的方差系数CV；设定一个阈值T_cv，它与SAR图像的方差系数呈T_cv＝β·CV的线性关系,其中β为线性系数；由阈值T_cv可得到二值化的方差系数矩阵CV_mask：

$\mathrm{CV}\_\mathrm{mask}=\left\{\begin{array}{cc}1,& \mathrm{CV}>T_{\mathrm{cv}}\\ 0,& \mathrm{CV}\_<=T_{\mathrm{cv}}\end{array}\right.$

可以计算出参考图和待匹配图二值化的方差系数矩阵，分别记为CV_mask_R和CV_mask_S；

3)对参考图的方差系数矩阵细化处理来得到细化矩阵T_R，同样对待匹配图方差系数矩阵细化处理得到细化矩阵T_S；

4)计算参考图和待匹配图的稀疏矩阵；

5)挑选稀疏矩阵中稀疏度值最大的前20％像素点作为稀疏度点，稀疏点膨胀处理即得到了稀疏区域。

上述步骤4)包括如下步骤：

(a)计算细化矩阵中各像素点梯度的大小和方向：

(b)分析各像素点梯度的方向，统计当前像素点的方向特征，得到参考图和待匹配图的方向特征矩阵，分别记为D_R、D_S；

(c)由细化矩阵和方向特征矩阵，可以计算参考图和待匹配图的稀疏度，稀疏度定义为：

$S=\frac{1}{N}\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{D}_{(i,j)}\&CenterDot;{{\text{T}}_{(i,j)}}^2$

其中，(i,j)是当前像素点的索引，Ω是以当前像素点为中心大小为r×r邻域，N为邻域Ω中像素点总数，D为方向特征矩阵，定义T为细化矩阵；由此再得到参考图和待匹配图的稀疏矩阵S_R、S_S。

上述步骤(b)包括如下步骤：

a)构造一个(2r+1)*(2r+1)大小的16方向模板，方向模板各个子方向对应的角度范围分别为i∈[1,16]，r为方向模板半径；

b)用16方向模板分析细化矩阵中各细化边缘点的方向特征；

c)对参考图和待匹配图的方向特征矩阵归一化处理，即可得到参考图和待匹配图的方向特征矩阵D_R、D_S。

本发明的有益效果：本发明是利用SAR图像方差系数判断出图像的平滑区域、纹理区域、稀疏结构区域，排除特征不明显且易造成误匹配的平滑区域和纹理区域，仅对那些折线拐点处、交叉处等稀疏度较高的像素点所在的稀疏区域提取SIFT特征点进行特征匹配，通过求取仿射变换参数来实现图像配准。本发明与现有的技术相比具有以下优点：

1.本发明是采用经典了的SIFT算子，在提取SIFT特征点的过程中的建立高斯尺度空间时，利用高斯函数模糊图像的同时也减少了SAR图像斑点噪声，从而提高了像素点在其邻域的可区分性，在一定程度上减少了误匹配的发生。

2.本发明设计了一种稀疏度定义法，能够比较准确的分析出折线、拐点、交叉点等特征明显的稀疏区域，排除平滑区域、纹理区域等稀疏度低且特征不明显的非稀疏区域。

3.本发明仅对参考图和待匹配图稀疏区域提取SIFT特征点，而排除平滑区域、纹理区域等稀疏度低、特征不明显的非稀疏区域，因而降低了误匹配概率，并且使算法更加快捷。

4.本发明实现简单、结构清晰，对参考图和待匹配图重叠面积要求小，并且SIFT算子对具有旋转、尺度变换的情形或包含变化区域的情形仍具有较高的配准精度。

附图说明

图1是本发明的流程图；

图2是本发明中计算稀疏区域的子流程图；

图3是本发明实验输入的两幅SAR图像，图3(a)为参考图，图3(b)为待匹配图；

图4是参考图和待匹配图的稀疏区域Mask_A、Mask_B；

图5是经典SIFT算法提取到的参考图和待匹配图的特征点分布；

图6是对参考图和待匹配图的稀疏区域Mask_A、Mask_B提取SIFT特征点分布；

图7是根据仿射变换参数得到的参考图和待匹配图的配准结果，用棋盘法显示。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，输入两幅在不同时间获取的同一地区的多时相SAR图像R和S，为了便于描述，我们称图像R为参考图，称图像S为待匹配图，分别计算两幅图像的稀疏区域Mask_R和Mask_S。本步骤参照图2，具体步骤如下：

(1a)分别对参考图R和待匹配图S进行两层小波分解，小波基取db1，得到参考图和待匹配图的细节分量分别记为R_db1、S_db1；

(1b)对参考图细节分量R_db1和待匹配图细节分量S_db1的每一个像素点计算方差系数CV。CV是概率分布离散度的度量，它被定义为标准差和均值的比：

$\mathrm{CV}=\frac{{\mathrm{\&sigma;}}_x}{{\mathrm{\&mu;}}_x}$

其中σ_x表示以当前像素点为中心大小为r×r邻域的标准差，_μx表示以当前像素点为中心大小为r×r邻域的均值。然后以当前像素点为中心、取一个大小为m×m的邻域计算方差，即可得到当前像素点的方差系数CV，本发明实例中r＝3、m＝5。用以上方法可得到参考图细节分量R_db1的方差系数矩阵CV_R；同理对于待匹配图细节分量S_db1有方差系数矩阵CV_S。设定一个阈值T_cv，它与SAR图像的方差系数CV有关(对于强度SAR图像对于幅度SAR图像其中L是SAR图像的视数)，一般经验性地把它看作是与方差系数呈线性关系，即T_cv＝β·CV，本发明中系数常量β值取为0.55。可得到二值化的方差系数矩阵CV_mask：

$\mathrm{CV}\_\mathrm{mask}=\left\{\begin{array}{cc}1,& \mathrm{CV}>T_{\mathrm{cv}}\\ 0,& \mathrm{CV}\_<=T_{\mathrm{cv}}\end{array}\right.$

可以计算出参考图和待匹配图二值化的方差系数矩阵，分别记为CV_mask_R和CV_mask_S。

(1c)对参考图CV_mask_R矩阵细化处理得到细化矩阵T_R，对待匹配图CV_mask_S矩阵细化处理得到细化矩阵T_S。

(1d)分别对参考图和待匹配图的稀疏矩阵。计算稀疏矩阵的具体步骤如下：

1d1)计算细化矩阵中各像素点的梯度大小和方向：

$m_{(x,y)}=\sqrt{{(d_{x+1,y}-d_{x-1,y})}^2+{(d_{x,y+1}-d_{x,y-1})}^2}$

$k_{(x,y)}={\tan }^{-1}\left(\frac{d_{x+1,y}-d_{x-1,y}}{d_{x,y+1}-d_{x,y-1}}\right)$

其中(x,y)表示的细化矩阵中当前像素点的索引，d为细化矩阵中像素点的灰度值，m_(x,y)为梯度大小，k_(x,y)为梯度方向。

1d2)构造15×15大小的16方向模板，方向模板各个子方向对应的角度范围分别为i∈[1,16]，方向模板上各个子方向上的像素点总数目为N，方向模板半径为r。

1d3)利用16方向模板分析细化矩阵每个细化边缘点的方向特征。在细化矩阵中，取以当前细化边缘点为中心大小为15×15的邻域，分析各子方向上每个像素点的梯度方向，统计梯度方向与该模板子方向相一致的像素点数目，记为s_i，统计可得当前细化边缘点归一化后的方向特征向量 $\mathrm{hist}=[\frac{s_1}{\mathrm{\&Sigma;}{s}_i}...\frac{s_i}{\mathrm{\&Sigma;}{s}_i}...\frac{s_{16}}{\mathrm{\&Sigma;}{s}_i}].$

1d4)计算参考图和待匹配图的方向特征矩阵。设定一个阈值N为方向模板上各个子方向上的像素点总数目，r为方向模板半径，本发明中阈值系数β取值为1.8。分析细化边缘点的方向特征，若h_i≥T_hist，则此方向为主方向之一，把每个边缘点的主方向个数作为其方向特征，可得到参考图和待匹配图的方向特征矩阵D_R、D_S。

1d5)计算参考图和待匹配图的稀疏度，定义像稀疏矩阵，

$S=\frac{1}{N}\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{D}_{(i,j)}\&CenterDot;{{\text{T}}_{(i,j)}}^2$

其中，(i,j)是当前像素点的索引，Ω是以当前像素点为中心大小为r×r邻域，N为邻域Ω中像素点总数，D为方向特征矩阵，定义T为细化矩阵。可得到参考图和待匹配图的稀疏矩阵S_R、S_S。

(1e)挑选稀疏矩阵中稀疏度值最大的前20％像素点作为稀疏度点，稀疏点膨胀处理即得到了稀疏区域。

步骤2，分别对参考图R和待匹配图S的稀疏区域提取SIFT特征点，特征点提取的具体步骤如下：

(2a)构建高斯尺度空间。一幅二维图像在不同尺度下的尺度空间可由原图像与高斯核函数卷积得到：

L(x,y,σ)＝G(x,y,σ)*I(x,y)

其中：L(x,y,σ)表示图像的尺度空间，I(x,y)表示原图像，(x,y)表示输入图像当前像素点的索引，*表示卷积运算，G(x,y,σ)为可变尺度高斯核函数，σ是尺度空间因子，核函数定义为：

$G(x,y,\mathrm{\&sigma;})=\frac{1}{2\mathrm{\&pi;}{\mathrm{\&sigma;}}^2}{e}^{-\frac{x^2+y^2}{2{\mathrm{\&sigma;}}^2}}$

对尺度空间进行离散化采样，即取σ的不同倍数下的高斯核与下采样图像卷积，所生成的一系列高斯图像，称为高斯金字塔(LOG)。高斯金字塔(LOG)分为O组，每组包含S+3层尺度空间，本发明中取O＝4、S＝2。

先求第一组高斯尺度空间中的每幅图像。输入图像I₀，经不同大小高斯核的核函数(高斯核分别为σ_i＝kⁱσ)卷积得到输出图像I_i，i＝0,12,3,4，即为高斯尺度空间的第一组图像。

第二组高斯金字塔的每层图像大小分别是第一组图像大小的1/4，其输入图像是由第一组LoG中第S层图像经采样率为2的下采样得到。如当S＝2时，第二组高斯金字塔的输入图像I′₀由第一组高斯金字塔I₁下采样得到，然后经不同大小高斯核的核函数(高斯核为σ_i＝kⁱσ)卷积得到输出I′_i，i＝0,12,3,4，即得到高斯尺度空间的第二组图。同理可得到高斯尺度空间的第三、四组图像。这样就生成了高斯尺度空间。

(2b)构建差分高斯尺度空间(DOG)。将高斯金字塔(LOG)中相邻两层的尺度空间函数进行相减就能够得到高斯差分金字塔(DOG)，用D(x,y,σ)表示：

D(x,y,σ)＝(G(x,y,kσ)-G(x,y,σ)×I(x,y))＝L(x,y,kσ)-L(x,y,σ)

其中，(x,y)表示输入图像中的当前像素点，I(x,y)表示待处理图像，即差分高斯尺度空间D(x,y,σ)由相差一个常数因子k的两层相邻高斯尺度空间函数相减得到。

(2c)极值点检测。将差分高斯尺度空间(DOG)里的每个像素点与同层的8个像素点、上层的9个像素点和下层的9个像素点共26个像素点比较大小，如果当前像素点为极大值或极小值，则把当前像素点定为关键点。

(2d)关键点的精确定位。由于边缘和噪声对差分高斯尺度空间(DOG)的影响较大，因此还需进一步的筛选上一步中检测到的关键点，精确定位具有尺度不变性的特征点，并且需要滤除低对比度点和边缘响应点，提高特征匹配的鲁棒性和稳定性。为了精确确定特征点的尺度和位置，须对检测出得局部极值点进行三维二次函数拟合。在局部极值点(x₀,y₀,σ)处，尺度空间函数D(x,y,σ)的泰勒展开式为：

$D\left(X\right)=D\left(X_0\right)+\frac{{\&PartialD;D}^T}{\&PartialD;X}(X-X_0)+\frac{1}{2}{(X-X_0)}^T\frac{{\&PartialD;}^{2}D}{{\&PartialD;X}^2}{(X-X_0)}^2$

上式中，X＝(x,y,σ)^T表示特征点的位置和尺度信息，X₀表示待选特征点的位置和尺度信息，一阶和二阶偏导计算如下，

$\frac{\&PartialD;D}{\&PartialD;X}=\left[\begin{array}{c}\frac{\&PartialD;D}{\&PartialD;x}\\ \frac{\&PartialD;D}{\&PartialD;y}\\ \frac{\&PartialD;D}{\&PartialD;\mathrm{\&sigma;}}\end{array}\right],\frac{{\&PartialD;}^{2}D}{{\&PartialD;X}^2}=\left[\begin{array}{ccc}\frac{{\&PartialD;}^{2}D}{{\&PartialD;x}^2}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{xy}}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{x\&sigma;}}\\ \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{yx}}& \frac{{\&PartialD;}^{2}D}{{\&PartialD;y}^2}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{y\&sigma;}}\\ \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;x}}& \frac{{\&PartialD;}^{2}D}{\&PartialD;\mathrm{\&sigma;y}}& \frac{{\&PartialD;}^{2}D}{\&PartialD;{\mathrm{\&sigma;}}^2}\end{array}\right]$

然后对此拟合函数进行求导并令一阶导数为零，即可得到x的极值x

$\hat{X}=-\frac{{\&PartialD;}^2{D}^{-1}}{{\&PartialD;}^2{X}^2}\frac{\&PartialD;D}{\&PartialD;X}$

将代入尺度空间函数D(x,y,σ)，可得

$D\left(\hat{X}\right)=D\left(X_0\right)+\frac{1}{2}\frac{{\&PartialD;D}^T}{\&PartialD;X}\hat{X}$

若则表示该点的对比度较低，此时将其剔除。

得到稳定精确的特征点的一个重要步骤，即是过滤掉边缘响应的特征点。过滤边缘响应特征点时，参照如下准则：在高斯差分空间函数D(x,y,σ)的边缘交叉处和峰值处，一个较大的主曲率值存在于图像边缘的特征点上，但在垂直方向主曲率的值较小。其中，求主曲率值可通过2×2的Hessian矩阵H，Hessian矩阵H，

$H=\left[\begin{array}{cc}{D}_{\mathrm{xx}}& D_{\mathrm{xy}}\\ D_{\mathrm{yx}}& D_{\mathrm{yy}}\end{array}\right]$

由于H矩阵的特征值和D的主曲率成比例关系，采用不具体求特征值只求其比例ratio的方法，α、β分别为最大特征值和最小的特征值，γ＝α/β，那么有：

$\mathrm{ratio}=\frac{\mathrm{Tr}{\left(H\right)}^2}{\mathrm{Det}\left(H\right)}=\frac{{(\mathrm{\&alpha;}+\mathrm{\&beta;})}^2}{\mathrm{\&alpha;\&beta;}}=\frac{{(\mathrm{\&gamma;}+1)}^2}{\mathrm{\&gamma;}}$

矩阵H的主对角元素之和Tr(H)的计算公式为：Tr(H)＝D_xx+D_yy＝α+β；矩阵H行列式的值Det(H)的计算公式为：Det(H)＝D_xxD_yy-(D_xy)²＝αβ。

因此我们只需要检测

$\frac{\mathrm{Tr}\left(H\right)}{\mathrm{Det}\left(H\right)}<\frac{{(\mathrm{\&gamma;}+1)}^2}{\mathrm{\&gamma;}}$

其中γ＝10。

(2e)生成特征描述子。为了保证特征点的旋转不变性，需要利用特征点周围领域像素点的梯度分布特性来计算特征点的主方向。每个像素点的梯度大小和方向计算如下：

$m(x,y)=\sqrt{{(L(x+1,y)-L(x-1,y))}^2+{(L(x+1,y)-L(x-1,y))}^2}$

$\mathrm{\&theta;}(x,y)=\arctan \left(\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)}\right)$

其中，特征点在(x,y)处梯度的大小为m(x,y)，方向为θ(x,y)，在金字塔中特征点(x,y)的尺度用L(x,y)来表示。在以当前特征点为中心的邻域内，采用梯度方向直方图统计邻域内像素的梯度方向。梯度方向范围为[0,360°)，其中每10°为一个直方图柱面。共有36柱面。梯度方向直方图的波峰代表了该特征点处邻域梯度的主方向，作为该特征点的主方向。若梯度方向直方图中存在主峰值能力80％以上的峰值，该峰值对应的方向作为辅方向，每个特征点有一个主方向和若干个辅方向。

确定特征点主方向之后，可以生成特征描述子。以特征点为中心取16×16窗口，在4×4的图像小块上计算8个方向的梯度方向直方图，绘制每个梯度方向的累加值。为每个4×4的子块统计梯度信息，16×16的邻域产生4×4个子块，因此特征描述子一共4×4×8＝128维。特征描述子由所有子块的梯度方向直方图构成，最终形成128维的SIFT特征描述向量。

步骤3，对参考图像和待匹配图像的SIFT特征点进行特征匹配，再利用随机一致性估计(RANSAC)方法滤除匹配结果中的误匹配对。具体步骤如下：

(3a)对参考图和待匹配图稀疏区域提取到的SIFT特征点进行特征点匹配，用最近邻距离与次近邻距离比(NNDR)准则来实现。假设A为参考图中的某个特征点，D_A其对应的特征向量。在待配准图中特征点A的最近邻(欧式距离最小)为特征点B和次近邻(欧氏距离次最小)为特征点C，它们对应的特征向量为D_B和D_C，最近邻距离与次近邻距离比(NNDR)准则定义为：

$\frac{|D_A-D_B|}{|D_A-D_C|}>T_{\mathrm{dis}}$

其中T_dis为距离比阈值，D_A表示参考图像中当前特征点A的128维特征描述子，D_B表示待匹配图中与A最相似的特征点B的128维特征描述子，D_C表示待匹配图中与A次相似的特征点C的128维特征描述子，当最近邻欧式距离和次近邻欧式距离之间的比值大于阈值T_dis时，则认为最近邻欧氏距离所对应的一对特征点是匹配特征点对。阈值取值越大，以为着匹配要求越宽松，错误匹配点存在概率增大；阈值取值越小意味着匹配的要求越严格，错误匹配点存在的概率减小，同时会去除一些正确的匹配。因此阈值取值根据所配准图像种类等实际情况设定。匹配完成后参考图和待匹配图中匹配中的点集分别记作{p_i}和{q_i}，即得到了粗匹配对{(pq)_i}，i表示匹配对的索引号。

(3b)因为杂乱背景和噪声等的影响，上一步得到SIFT匹配对中常常存在许多错误匹配点对，它们会引起图像几何校正的误差，我们用随机一致性估计(RANSAC)方法去除粗匹配结果中的错误匹配对，去除误匹配后参考图和待匹配图中匹配对中的点集分别记作{P(x_i,y_j)}和{Q(x′_i,y′_i)}，i表示匹配对的索引号。

步骤4，利用参考图和待匹配图的匹配点集{P(x_i,y_i)}和{Q(x′_i,y′_i)}，通过仿射变换函数得到精确的仿射变换参数，放射变换模型如下：

$\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right]\&CenterDot;\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\end{array}\right]+\left[\begin{array}{c}{c}_1\\ c_2\end{array}\right]$

变换矩阵M为：

$M=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\end{array}\right]$

其中，a₁、a₂、b₁、b₂为旋转变量，c₁为水平偏移量，c₂为垂直偏移量。最后根据求解的仿射变换矩阵对图像进行配准。

本发明效果可以通过以下实验进一步证实：

一.实验条件和内容

硬件平台为：Intel(R)Pentium(R)1CPU2.4GHz；

软件平台为：WindowXPProfessional，MATLABR2010；

实验所使用的输入图像如图3所示，其中图3(a)和图3(b)是SAR图像。

实验内容：

仿真内容1，在上述实验条件下，分别利用经典SIFT算法、本发明中的快速配准算法提取SIFT特征点，再根据最近次近距离比值准则(NNDR)，在不同距离比阈值T_dis的条件下对两种方法提取到的特征点进行特征点匹配。

仿真内容2，在得到的粗匹配结果基础上，用随机一致性估计算法(RANCAC)去除误匹配。利用最后保留下来的匹配对来求仿射变换参数，最后实现图像的配准。

考虑到如果距离比阈值T_dis取得太小，粗匹配结果中正确匹配对的数目就会太少，无法满足仿射变换的要求；如果距离比阈值T_dis取得太大，粗匹配结果中正确匹配对在所有匹配对中所占有的比例太低，随机一致性估计(RANCAC)算法鲁棒性太差。所以实验中距离比阈值T_dis选取如下几个值：0.70、075、0.80、0.85。

二.实验结果

本发明得到的稀疏区域如图4所示，图4(a)为参考图的稀疏区域；图4(b)为待匹配图的稀疏区域。从图中可以看出稀疏区域区域大多集中在折线拐点、交叉处等区域，这些区域稀疏度高、特征明显。

使用经典SIFT特征提取的方法对参考图和待匹配图提取特征点的结果如图5所示，图中白色点表示提取到特征点，图5(a)为参考图的SIFT特征点分布；图5(b)为待匹配图的SIFT特征点分布。从结果图可以看出来，经典SIFT特征提取方法能对图像中线目标区域、交叉处区域等稀疏度高的区域检测出较多的SIFT特征点，对特征不明显的平滑区域和纹理区域也能检测到一些特征点，而这些平滑区域和纹理区域的点周围往往存在很多与其相似的点，因而容易形成误匹配。

使用基于空间稀疏度与SIFT特征提取的快速图像配准方法对参考图和待匹配图提取特征点的结果如图6所示，图中白色点表示提取到特征点，图6(a)为参考图的特征点分布；图6(b)为待匹配图的特征点分布。从结果图可以看出来基于空间稀疏度与SIFT特征提取的快速图像配准方法提取到的特征点仅存在于稀疏度较高且特征明显的折线拐点、交叉处区域，这些区域的点稀疏度高、特征稳定、不易形成误匹配。

使用基于欧式距离比的特征匹配方法对两种方法提取到的特征点进行粗匹配。再用随机一致性估计法RANSAC对粗匹配结果去除误匹配。两种方法粗匹配和去除误匹配的结果列在表1中，用匹配正确率来作为匹配效果好坏的评价指标。匹配正确率的定义：

表1  两种方法特征点匹配的结果

从表1中可以发现，在不同最近次近距离比阈值T_dis＝0.70、T_dis＝0.75、T_dis＝0.80、T_dis＝0.85下，本发明中正确匹配对在所有匹配对中所占的比例都比经典SIFT算法都高。并且，在粗匹配结果中误匹配率较高的情况下本发明的相对于经典SIFT算法的优势更加明显。从匹配正确率统计可以看出，在不同欧氏距离比阈值条件下，本发明匹配正确率比经典SIFT算法的匹配正确效率更高。

通过仿射变换得到仿射变换参数，取阈值为T_dis＝0.75时最终匹配结果得到的仿射变换参数如表2所示：

表2：仿射变换参数结果

从表2中可以发现，本发明方法得到的仿射变换参数与SIFT算法结果相差不大，但是在时间效率上有了很大的改善，本发明中的算法时间由原来的的81.04秒减少到了20.06秒，缩短到了原来的四分之一，大大提高了配准效率，阈值为T_dis＝0.75时，两幅图像的配准结果如图7所示，用棋盘法显示。因为本发明仅对SAR图像的折线、拐点、交叉点等特征明显的稀疏区域提取SIFT特征点，并且排除平滑区域、纹理区域等稀疏度低、特征不明显、易形成错误匹配的区域，因而在降低误匹配率的前提下提高了配准效率。

本发明设计出了一种新的基于空间稀疏度和SIFT特征特征提取的快速图像配准方法，先分析出图像的稀疏度矩阵，对图像中特征明显的稀疏区域提取SIFT特征点，在保证配准精度的同时提高了算法的效率。

综上，本发明是利用SAR图像方差系数判断出图像的平滑区域、纹理区域、稀疏结构区域，排除特征不明显且易造成误匹配的平滑区域和纹理区域，仅对那些折线拐点处、交叉处等稀疏度较高的像素点所在的稀疏区域提取SIFT特征点进行特征匹配，通过求取仿射变换参数来实现图像配准。本发明与现有的技术相比具有以下优点：

1.本发明是采用经典了的SIFT算子，在提取SIFT特征点的过程中的建立高斯尺度空间时，利用高斯函数模糊图像的同时也减少了SAR图像斑点噪声，从而提高了像素点在其邻域的可区分性，在一定程度上减少了误匹配的发生。

2.本发明设计了一种稀疏度定义法，能够比较准确的分析出折线、拐点、交叉点等特征明显的稀疏区域，排除平滑区域、纹理区域等稀疏度低且特征不明显的非稀疏区域。

3.本发明仅对参考图像和待匹配图像稀疏区域提取SIFT特征点，而排除平滑区域、纹理区域等稀疏度低、特征不明显的非稀疏区域，因而降低了误匹配概率，并且使算法更加快捷。

4.本发明实现简单、结构清晰，对参考图和待匹配图重叠面积要求小，并且SIFT算子对具有旋转、尺度变换的情形或包含变化区域的情形仍具有较高的配准精度。

本实施方式没有详细叙述的部分属本行业的公知的常用手段，这里不一一叙述。以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (4)

1.一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，其特征在于：包括如下步骤：

(1)输入在不同时间获取的同一地区的两幅多时相SAR图像，分别记作参考图和待匹配图，对参考图和待匹配图的每个像素点计算稀疏度，挑选稀疏度高的像素点所在的区域作为稀疏区域，分别记作Mask_R和Mask_S；

(2)分别对参考图稀疏区域Mask_R和待匹配图稀疏区域Mask_S提取SIFT特征点，得到参考图和待匹配图SIFT特征点集分别记为{R_m}、{S_n}；

(3)对参考图和待匹配图上提取到的SIFT特征点集{R_m}和{S_n}进行特征匹配，得到粗匹配特征对{(pq)_i}，参考图和待匹配图中的粗匹配点集分别为{p_i}和{q_i}，i表示匹配对的索引号；

(4)利用随机一致性估计算法对粗匹配结果去除误匹配对，去除误匹配对后参考图和待匹配图保留下来的特征点对集合记为{(PQ)_i}，此时参考图和待匹配图中各自保留下来特征点分别为{P_i}、{Q_i}，其中i表示匹配对的索引号；

(5)利用参考图和待匹配图的匹配点集{P_i}和{Q_i}，通过仿射变换得到仿射变换参数，最后根据仿射变换参数对两幅图像配准。

2.根据权利要求1所述的一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，其特征在于：所述的步骤(1)，包括如下步骤：

1)分别对参考图R和待匹配图S进行两层小波分解，小波基取db1，得到参考图和待匹配图的细节分量图，分别记为R_db1、S_db1；

2)对参考图细节分量图和待匹配图细节分量图的每一个像素点计算方差系数，然后以当前像素点为中心、取一个大小为m*m的邻域计算方差，即可得到当前像素点的方差系数CV；设定一个阈值T_cv，它与SAR图像的方差系数呈T_cv＝β·CV的线性关系,其中β为线性系数；由阈值T_cv可得到二值化的方差系数矩阵CV_mask：

$\mathrm{CV}\_\mathrm{mask}=\left\{\begin{array}{c}1,\mathrm{CV}>T_{\mathrm{cv}}\\ 0,\mathrm{CV}<=T_{\mathrm{cv}}\end{array}\right.$

计算出参考图和待匹配图二值化的方差系数矩阵，分别记为CV_mask_R和CV_mask_S；

3)对参考图的方差系数矩阵细化处理来得到细化矩阵T_R，同样对待匹配图方差系数矩阵细化处理得到细化矩阵T_S；

4)计算参考图和待匹配图的稀疏矩阵；

5)挑选稀疏矩阵中稀疏度值最大的前20％像素点作为稀疏度点，稀疏点膨胀处理即得到了稀疏区域。

3.根据权利要求2所述的一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，其特征在于：所述的步骤4)，包括如下步骤：

(a)计算细化矩阵中各像素点梯度的大小和方向：

(b)分析各像素点梯度的方向，统计当前像素点的方向特征，得到参考图和待匹配图的方向特征矩阵，分别记为D_R、D_S；

(c)由细化矩阵和方向特征矩阵，可以计算参考图像和待匹配图像的稀疏度，稀疏度定义为：

$S=\frac{1}{N}\underset{(i,j)\&Element;\mathrm{\&Omega;}}{\mathrm{\&Sigma;}}{D}_{(i,j)}\&CenterDot;{T_{(i,j)}}^2$

其中，(i,j)是当前像素点的索引，Ω是以当前像素点为中心大小为r×r邻域，N为邻域Ω中像素点总数，D为方向特征矩阵，定义T为细化矩阵；由此再得到参考图和待匹配图的稀疏矩阵S_R、S_S。

4.根据权利要求3所述的一种基于空间稀疏度和SIFT特征提取的快速图像配准实现方法，其特征是：所述的步骤(b)，包括如下步骤：

a)构造一个(2r+1)*(2r+1)大小的16方向模板，方向模板各个子方向对应的角度范围分别为i∈[1,16]，r为方向模板半径；

b)用16方向模板分析细化矩阵中各细化边缘点的方向特征；

c)对参考图和待匹配图的方向特征矩阵归一化处理，即可得到参考图和待匹配图的方向特征矩阵D_R、D_S。