一种鲁棒SURF无人机彩色遥感图像配准方法与流程

本发明涉及无人机像配准技术，尤其是一种遥感图像配准方法。

背景技术：

surf(speededuprobustfeatures)算法是一种具有鲁棒性的局部特征检测算法，该算法不仅具有重复性高的检测器和可区分性强的描述符，此外具有很强的鲁棒性及较高的运算速度；在精度上能达到亚像素级别，且在尺度缩放、亮度变化、旋转、遮挡和噪声等情况下均能保证从图像中提取出的特征描述子具有良好不变性，目前该算法已被广泛应用于图像配准技术中。但由于无人机遥感图像像幅小、数量多，对于配准精度的要求极其苛刻。如何增强彩色图像特征点显著性，从而进一步提高配准精度，尤其是特征相似度较高的彩色遥感图像精确配准一直是亟待解决的难题。

一些研究人员提出了将rgb分量的强度归一化后再用surf算法来处理彩色遥感图像的方法和基于surf的颜色不变特征的方法，但第一种方法在描述色彩不变特征时有局限性，无法对彩色信息进行充分利用，第二种方法耗时较长，实时性差，均不能满足实际要求。

技术实现要素：

要解决的技术问题

本发明针对无人机彩色遥感图像使用surf算法配准时仅利用图像亮度信息忽视彩色信息使得配准精度降低的问题，提出一种鲁棒surf无人机彩色遥感图像配准方法。

技术方案

一种鲁棒surf无人机彩色遥感图像配准方法，其特征在于步骤如下：

步骤1：对获取的无人机彩色遥感参考图像image1、待配准图像image2进行灰度处理变为灰度图像，对灰度图像进行高斯平滑处理完成尺度空间构建，采用特征检测算子hessian矩阵的行列式进行特征点提取，获取image1、image2的所有特征点分别记为目标集{pi},i＝1,2,……,n和基准集{qj},j＝1,2,……,m；

步骤2：对于得到的目标集{pi}和基准集{qj}，以每个特征点为中心确定边长为20s的正方形区域，其中s表示不同尺度空间下的采样间隔，将其周围邻域信息分成4×4个子区域，每个小区域又分为5×5个采样点，最后用harr小波计算每个小区域垂直和水平方向的响应及模值dx,|dx|,dy,|dy|，并统计5×5个采样点的总响应∑dx,∑|dx|,∑dy,∑|dy|，可得到特征点的64维surf特征描述符：

v＝(∑dx,∑|dx|,∑dy,∑|dy|)

则特征点的surf描述符在水平或垂直方向上定义为：

vgray＝(v1,v2,v3,v4,...,v16)

其中，vi,i＝1,2,3,...,15,16表示16个子区域所对应的描述向量；由此可得目标集{pi}和基准集{qj}中所有的特征点的64维surf特征描述符；

步骤3：对于坐标为(x,y)特征点，其在rgb颜色空间表示为r(x,y)、g(x,y)、b(x,y)，将其进行颜色空间转换，转换到yiq颜色空间为y(x,y)、i(x,y)、q(x,y)；

rgb颜色空间变换到yiq颜色空间表示为：

将其增加到原始算法描述符中，得到新的彩色描述符如下：

对上式进行归一化处理：

其中ik,k＝1,2,3,...,63,64代表经典surf算法中特征点的64维描述向量；由此可得目标集{pi}和基准集{qj}中各个特征点的67维描述向量；

步骤4：特征点双向匹配

步骤4a：使用欧氏距离公式计算参考图像image1中的特征点p1在image1中正向的最近欧式距离dnt和次近欧式距离dsnt；

使用欧氏距离公式计算参考图像image1中的特征点p1在image2中正向的最近欧式距离d1nt和次近欧式距离d1snt的特征点分别为p2nt和p2snt；

步骤4b：计算最近距离d1nt和次最近距离d1snt的比率t1，即

将t1与判断阈值t进行比较：如果t1<t，随之进行步骤4c，否则，自动跳出；

步骤4c：对于p2nt对应特征点，使用欧氏距离公式计算在image1中反向的最近距离d2nt和次最近距离d2snt的特征点分别为p1nt和p1snt；

步骤4d：计算最近距离和次最近距离的比率t2

将t1与判断阈值t进行比较：如果t2<t，同时p1nt所寻找的特征点和原始特征点p1是同一个特征点，则认为是正确匹配点对，否则返回步骤4a对参考图像image1中的其它特征点进行新的特征点判断；

依次遍历参考图像目标集中所有特征点，最终实现双向匹配提纯；

步骤5：采用ransac算法实现变换矩阵参数求取

步骤5a：寻找步骤4完成后的双向匹配提纯后的特征点对mi'←→mi，i＝1,2,......n；

步骤5b：从匹配点对集合中，任意选取4对匹配点对，使用这4对匹配点对的坐标值，求得变换矩阵h；

步骤5c：以欧氏距离为依据，在匹配点对集合中寻找符合d(himi,mi')<t条件的点对，其中，himi表示对image2中特征点进行矩阵变换映射到image1中的位置坐标，mi'表示mi在image1中对应特征点的位置坐标，d(himi,mi')表示两个坐标向量的欧氏距离，将符合条件的点对作为最终的内点，并记录满足hi约束的内点数量；

步骤5d：重复5b和5c两步n次，记录每一次的内点数量；

步骤5e：选取对应内点数量最大的hbest，寻找所有符合d(hbestmi,mi')<t条件的匹配点对，将它们作为最终的内点，即正确的匹配点对，不符合d(hbestmi,mi')<t条件的错误匹配点对，即为外点，予以剔除；

步骤5f：根据随机抽样一致性算法求得n个匹配点对，对这个2n个由x1,y1,x2,y2所构成的矩阵进行标记，记作a，对其进行svd奇异值分解得a＝udv^t，u和v中分别是a的奇异向量，d为对角线上的奇异值按降序排列，v的最后一列重构为3*3矩阵即为所求的透视变换矩阵参数估计；

步骤6：双三次插值实现配准

若参考图像image1上任意一点坐标为(x,y)，其灰度值为g(x,y)，该点在待配准图像image2上的对应点坐标为(u,v)，其灰度值为f(u,v)，[·]表示取整数，待配准图像中对应点4*4邻域内16个像素点组成的矩阵为b：

矩阵b中元素中函数f(u+i,v+j),i＝-1,0,1,2；j＝-1,0,1,2,定义为在待配准图像上的对应点坐标为(u+i,v+j)时的灰度值；

则参考图像中待插值点的灰度值计算公式如下：

g(x,y)＝f(u,v)＝abc

其中：

s(w)为双三次插值的基函数，是对sin(x*π)/x的逼近，通过双三次插值，将image2图像插值到参考图像image1中，实现最终的高精度配准。

所述的步骤4b中的判断阈值t为0.4-0.8。

所述的步骤5c中的t选取4。

所述的步骤5d中的100≤n≤200。

有益效果

本发明提出的一种鲁棒surf无人机彩色遥感图像配准方法，即在特征点描述阶段，对特征点自身进行颜色空间变换，得到特征点自身的彩色信息，随后在描述算子中进行表述。匹配效果良好，图像没有出现交叉线。

附图说明

图1为本发明无人机彩色遥感图像image1和image2，其中图1(a)为参考图像image1，图1(b)为待配准图像image2。

图2为原surf算法的配准结果。

图3为本发明一种鲁棒的surf配准结果。

具体实施方式

现结合实施例、附图对本发明作进一步描述：

由于rgb颜色空间的分量与明亮程度关系紧密，而通常采集到的彩色遥感图像对亮度过于敏感，所以rgb颜色空间不适合进行彩色图像处理。yiq颜色空间(y分量代表图像的亮度信息，i、q两个分量则携带颜色信息)和rgb颜色空间通过一个矩阵相乘可直接相互变换，且为线性变换，实时性高。另一方面，yiq具有能将亮度分量提取出来，同时彩色信息也被分离出来，计算量小，聚类特性好，因此能够有效地用于彩色图像处理。本方法在特征点描述阶段将特征点由rgb空间变换到yiq空间中，对特征点进行yiq模型表示，将特征点自身彩色信息添加到描述符中，在原经典surf的64维描述符基础上增加三维彩色信息向量构成67维描述符以增强彩色图像特征点描述符区分度。随后结合ransac(randomsampleconsensus)、双向匹配、双三次插值算法实现无人机彩色遥感图像精确配准。主要包括如下步骤：

步骤一、利用surf算法进行特征点检测获取目标集和基准集

令无人机彩色遥感图像image1为参考图像，其中的特征点描述子集为目标集{pi},i＝1,2,……,n，另一幅彩色遥感图像image2为待配准图像，其中的特征点描述子集为基准集{qj},j＝1,2,……,m；

对获取的无人机彩色遥感图像image1、image2进行灰度处理，变为灰度图像，进行高斯平滑处理完成尺度空间构建，采用特征检测算子hessian矩阵的行列式进行特征点提取，获取image1、image2的所有特征点分别记为目标集{pi},i＝1,2,……,n和基准集{qj},j＝1,2,……,m，具体描述如下：

即对于灰度处理后所得到的灰度图像i(x,y)中的某一点x＝(x，y)，在尺度σ上的hessian矩阵定义为：

其中：lxx(x,σ)是尺度为σ的高斯二阶微分和图像i(x,y)在点x处卷积的结果：

其中，代表了二维卷积运算，同理，lxy(x,σ)和lyy(x,σ)分别是尺度为σ的高斯二阶微分和i(x,y)在点x处卷积。

hessian矩阵的行列式：

det(h)＝lxx(x,σ)lyy(x,σ)-l²xy(x,σ)(3)

不同尺度σ代表不同空间分辨率，不同尺度σ值进行高斯平滑处理后的图像不同，进行尺度空间构建：

组间尺度：

σoctave＝σ0*2^octave-1(4)

组内尺度：

σ-----尺度空间坐标，即为尺度空间的一个数值表征量；

s-----组内尺度空间坐标，即为组内尺度空间的一个数值表征量；

σ0-----初始尺度，根据相机性能决定，一般选1.6；

s-----每组层数(4层)；

octave-----组数(4组)；

由式(4)、(5)建立尺度空间，然后通过三维非最大值抑制搜索技术在图像的尺度空间中搜索行列式响应的局部极值完成特征点的检测；

步骤二、特征点描述及其添加彩色信息

2.1特征点描述

对于得到的目标集{pi}和基准集{qj}，以每个特征点为中心确定边长为20s(s表示不同尺度空间下的采样间隔)的正方形区域，将其周围邻域信息分成4×4个子区域，每个小区域又分为5×5个采样点，最后用harr小波计算每个小区域垂直和水平方向的响应及模值dx,|dx|,dy,|dy|，并统计5×5个采样点的总响应∑dx,∑|dx|,∑dy,∑|dy|，可得到特征点的64(4×4×4)维surf特征描述符：

v＝(∑dx,∑|dx|,∑dy,∑|dy|)(6)

假设某特征点坐标为(x,y)，则该特征点的surf描述符在一个方向上(水平或垂直)可定义为：

vgray＝(v1,v2,v3,v4,...,v16)(7)

其中，vi(i＝1,2,3,...,15,16)表示16个子区域所对应的描述符。利用上述方法可得目标集{pi}和基准集{qj}中所有的特征点的64维surf特征描述符。

2.2彩色信息添加

将特征点自身的彩色信息加入到原有描述符中，坐标为(x,y)特征点，其在rgb颜色空间表示为r(x,y)、g(x,y)、b(x,y)，将其进行颜色空间转换，转换到yiq颜色空间为y(x,y)、i(x,y)、q(x,y)，rgb颜色空间在日常生活中最为常用。然而，我们通常采集到的彩色遥感图像对亮度过于敏感，而rgb颜色空间的分量与明亮程度关系紧密。所以，rgb颜色空间不适合进行彩色图像处理。而yiq和rgb通过一个矩阵相乘可直接相互变换，成线性变换，实时性高。另一方面，yiq颜色空间具有能将亮度分量提取出来，同时彩色信息也被分离出来，计算量小，聚类特性好，因此能够有效地用于彩色图像处理，这也是我们改进的根源所在。

rgb颜色空间变换到yiq颜色空间表示为：

将其增加到原始算法描述符中，得到新的彩色描述符如下：

vcolor＝(v1,v2,v3,v4,...,v16,y(x,y),i(x,y),q(x,y))(9)

为了使得表述更加清晰，将式(6)进行简单变换

v＝(i1,i2,i3,i4)(10)

随后对公式(9)式进行重新定义如下：

vcolor＝(i1,i2,i3,i4,...,i64,y(x,y),i(x,y),q(x,y))(11)

为了对旋转、尺度、光照有更好的鲁棒性。这里需要对式(11)进行归一化处理：

其中|vcolor|是(11)式的求模式，表达式为：

其中ik(k＝1,2,3,...,63,64)代表经典surf算法中特征点的64维描述符。从公式(11)中我们可以看到，在原始的64维的特征点描述符的基础上，特征点自身的彩色信息也被加载进来，构成了67维的描述符，和原来的算法相比，时间上没有太大变化，但性能上因为彩色信息的加入，使得精度会进一步提升。

以此类推，最终求得目标集{pi}和基准集{qj}中各个特征点的67维描述符。

步骤三、特征点匹配

对于参考图像image1中目标集中每一个特征点，在基准集中查询得到它的最近邻欧氏距离dnt和次近邻欧氏距离dsnt，若满足：

则保留符合公式(14)的目标集中的特征点与其目标集中该特征点在基准集中查询到的最近邻的特征点最近邻构成的匹配，否则剔除这个匹配对，公式中的t是判断阈值，在0.4-0.8之间取值。

通过以上步骤我们得到了最近邻欧氏距离与次近邻欧式距离比值度量准则下的匹配点对。这些匹配点对出现误匹配的可能性依然很高，需要进行提纯处理，这里用双向匹配的方法，操作如下：

为使得表述更加清晰，这里将image1和image2中的描述符式(12)记为v1＝(v11,v12,v13,v14,...,v167)和v2＝(v21,v22,v23,v24,...,v267)，由欧氏距离公式可以得到：

其中，v1、v2两个向量之间的欧氏距离表示为d(v1,v2)，v1i、v2i(i＝1,2,3,...,66,67)表示向量v1、v2的各分量值；

最近距离和次近距离表示为dnt和dsnt，ti是最近距离与次近距离的比率，表示如下：

根据公式(15)和(16)进行双向匹配点对提纯，具体步骤如下：

(1)参考图像image1中的特征点p1，用欧氏距离公式在image2中寻得其最近欧式距离和次近欧式距离特征点分别为p2nt和p2snt，用d1nt和d1snt来描述正向的最近距离和次最近距离，其中从参考图像目标集特征点在基准集中寻找的方向，定位正向，即确定目标集中一个特征点在基准集中进行遍历寻找满足条件的点；而从基准集在目标集中寻找的方向，定为反向；

(2)使用公式(16)计算最近距离和次最近距离的比率t1，即

如果t1<t，随之进行步骤(3)，否则，自动跳出；

(3)对于p2nt对应的特征点，用公式(15)在image1中寻得最近距离和次最近距离特征点分别为p1nt和p1snt，用d2nt和d2snt来描述反向的最近距离和次最近距离；

(4)使用公式(16)计算最近距离和次最近距离的比率t2

如果t2<t，同时p1nt所寻找的特征点和原始特征点p1是同一个特征点，则认为是正确匹配点对，否则返回步骤(1)对参考图像image1中的其它特征点进行新的特征点判断；

依次遍历参考图像目标集中所有特征点，最终实现双向匹配提纯。

步骤四、ransac算法求取变换矩阵参数

首先对变换矩阵求解，给定彩色遥感图像image1和image2上两点m'和m，根据透视变换矩阵h有：m'～hm，其中～表示成比例相等，设m和m'的坐标分别为(x1,y1)和(x2,y2)，写成齐次坐标形式为：(x1,y1,z1)和(x2',y2',z2)，其中x2'/z2＝x2，y2'/z2＝y2，那么则有公式

进而在齐次坐标的基础上得到(20)和(21):

不失一般性而且为了计算简便，令z1＝1，由(20)和(21)可得：

x1h11+y1h12+h13-x1x2h31-y1x2h32-x2h33＝0(22)

x1h21+y1h22+h23-x1y2h31-y1y2h32-y2h33＝0(23)

由此推导出：

由于透视变换模型最后一位数h33通常被我们设置为1，所以变换模型由八个参数组成，至少需要四对不共线的匹配点对，便可以求得透视变换模型各参数的值。但是不管是在实际计算过程中还是客观条件中，误差时一定存在的。首先，在实际计算中，由于噪声干扰、特征点检测的定位误差、选取的模型误差、错配等因素，只依靠四对不共线的匹配点对求得八个参数的话误差较大，无法保证配准精度，ransac算法恰好能够很好的解决上述问题。

ransac算法作为一种选择机制，它不但尽可能保证了提取匹配点对的个数和质量，而且可以对变换模型进行精确的估计，ransac算法实现变换矩阵参数求取步骤如下所示：

(1)给定的两幅彩色遥感图像image1和image2，寻找双向匹配提纯后的特征点对mi'←→mi，i＝1,2,......n；

(2)从匹配点对集合中，任意选取4对匹配点对，使用这4对匹配点对的坐标值，求得变换矩阵h；

(3)以欧氏距离为依据，在匹配点对集合中寻找符合d(hmi,mi')<t条件的点对，其中，hmi表示对image2中特征点进行矩阵变换映射到image1中的位置坐标，mi'表示mi在image1中对应特征点的位置坐标，d(hmi,mi')表示两个坐标向量的欧氏距离，将符合条件的点对作为最终的内点，并记录满足h约束的内点数量，t选取4；

(4)重复(2)和(3)两步n次(100≤n≤200)，记录每一次的内点数量；

(5)选取对应内点数量最大的hbest，寻找所有符合d(hbestmi,mi')<t条件的匹配点对，将它们作为最终的内点，即正确的匹配点对，不符合d(hbestmi,mi')<t条件的错误匹配点对，即为外点，予以剔除；

(6)根据随机抽样一致性算法求得n个匹配点对，我们对2n个包含匹配点对及对应变换矩阵h的方程所构成的矩阵记作a，对其进行svd奇异值分解得a＝udv^t，u和v中分别是a的奇异向量，d为对角线上的奇异值按降序排列，v的最后一列重构为3*3矩阵即为所求的透视变换矩阵参数估计，这个矩阵的参数满足最小二乘意义下误差最小；

步骤五、双三次插值实现配准

当获得变换矩阵后需要进行映射实现配准，在映射的过程中，由于像素都是以整数坐标值进行存储，为一个个非连续的值，这便会导致一种必然出现的情况发生：即原来在整数网格上的点(x1,y1)(坐标都是整数)，在映射之后(x2,y2)没有落在网格点上，因此插值的作用便显得尤为重要。

双三次插值法利用待插值点在待配准图像中的对应点4*4邻域内的十六个像素点的灰度信息来估计待插值点的灰度值，不仅考虑了四个直接相邻点灰度值的影响，还考虑了各相邻点灰度值变化率的影响。如果待插值点在待配准图像中的对应点附近存在4*4邻域内的十六个像素点，就可以采用此方法。

若参考图像上任意一点坐标为(x,y)，其灰度值为g(x,y)，该点在待配准图像上的对应点坐标为(u,v)，其灰度值为f(u,v)，[·]表示取整数，待配准图像中对应点4*4邻域内16个像素点组成的矩阵为b：

矩阵b中元素中函数f(u+i,v+j)(i＝-1,0,1,2；j＝-1,0,1,2)定义为在待配准图像上的对应点坐标为(u+i,v+j)时的灰度值。

则参考图像中待插值点的灰度值计算公式如下：

g(x,y)＝f(u,v)＝abc(26)

其中：

s(w)为双三次插值的基函数，是对sin(x*π)/x的逼近。通过双三次插值，将image2图像插值到参考图像image1中，实现最终的高精度配准。

双三次插值法采用了较大邻域内像素点的灰度信息，插值精度进一步提高，并且对图像的边缘细节也有很好的保持，视觉效果更加理想，基本是现有的插值算法中插值精度最高的算法。

以下实例对相同的两幅无人机彩色遥感图像image1和image2用原surf算法和一种鲁棒surf算法进行精度仿真实验验证，匹配效果良好，鲁棒surf算法没有出现交叉线(见附图2和附图3)。image1和image2之间实际上是一种经过2倍缩放，且有10度旋转，水平偏移239个像素点，纵向偏移28个像素点仿射变换。

实例：

首先根据前面所述的彩色遥感图像image1和image2之间几何变换关系，我们可以很容易的得到理论值hture。其次，我们用原surf算法对彩色遥感图像image1和image2进行配准，得到原始hsurf。最后用一种鲁棒surf算法对彩色遥感图像image1和image2进行配准，得到hrob。

不失一般性，随后选取八组无人机彩色遥感图像分别进行试验，这里进行8组实验求取平均值，对变换矩阵各参数进行统计如表1所示，表1是仿射变换矩阵参数对比表，发现surf算法其最大偏差平均值为0.218高于鲁棒的最大偏差平均值0.134，六参数接近程度，鲁棒算法明显优于surf算法。

表1仿射变换矩阵参数对照表