CN105869168A - 基于多项式拟合的多源遥感图像形状配准方法 - Google Patents

Abstract

本发明提供一种新的基于多项式拟合的多源遥感图像形状配准方法。本发明方法首先提取遥感图像边缘轮廓特征，保留主轮廓边缘，并通过一种基于多项式拟合的特征点提取算法提取特征点并确定特征点的主方向；在此基础上，通过改进的形状内容描述子并最小化特征点之间的匹配代价来完成粗配准；最后去除误匹配，完成精配准。本发明方法能够实现对不同地形的多源遥感图像实现自动配准，具有配准精度高、计算速度快、鲁棒性和适用性好等优点。

Description

基于多项式拟合的多源遥感图像形状配准方法

技术领域

本发明涉及一种遥感图像处理技术，特别是一种基于多项式拟合的多源遥感图像形状配准方法。

背景技术

遥感技术经过半个多世纪的发展，在地学等领域得到了广泛地应用，为勘探资源、保护环境、监测灾害、分析全球变化等开辟了新的途径，成为人们观察、分析地球环境的有效手段。由于不同遥感卫星获得的遥感图像在分辨率、视角、尺度等方面都具有很大差异，故对遥感图像进行配准来尽可能地消除差异是后续分析应用的基础。根据配准过程中所利用的图像信息的不同，可将图像配准分为基于灰度的配准方法和基于特征的配准方法。基于灰度的配准方法主要包括互信息法和相位相关法等，互信息法因其配准精度高，常用于多模态医学图像的配准，近年来也常用于遥感图像的配准，但计算速度限制了它在遥感图像配准领域的发展和应用；相位相关法易于硬件实现.但它要求图像间有严格的线性关系.且对噪声非常敏感，局限性大。基于特征的配准方法不依赖于图像的灰度特性，将对整个图像的分析转化为对图像某种特征的分析，大大降低了计算量，是目前使用最多的遥感图像配准方法。但是在对不同类型的传感器获得多源遥感图像进行配准时，由于目前使用较多的基于特征的遥感图像配准算法都需要提取图像中明显的特征，如Harris、SIFT特征点或者直线特征等，此时因成像机理及像元表示形式不同使得图像间呈现很大的灰度和噪声差异，或使匹配精度很低，甚至匹配失败。

发明内容

本发明的目的在于提供一种基于多项式拟合的多源遥感图像形状配准方法，解决了传统方法对多源遥感图像进行配准时匹配精度低，运算耗时、适用局限性的问题。该方法包括以下步骤：

步骤1，输入两幅多源遥感图像作为参考图像和待配准图像；

步骤2，对两幅图像提取边缘轮廓特征，且消除密集区域无用的边缘特征；

步骤3，针对两幅图像提取到的轮廓上的像素点分别通过一种基于多项式拟合的特征点提取算法提取特征点并确定特征点的主方向；

步骤4，获取每一特征点的改进的形状内容描述子；

步骤5，基于像素点的改进的形状内容描述子，对参考图像中的每一个特征点，遍历待配准图像中特征点，获取与其匹配代价最小的待配准图像中特征点作为匹配点对；

步骤6，去除误匹配，完成精配准；

步骤7，根据匹配点对估计两图像间变换矩阵，将待配准图像经变换矩阵还原，完成配准；

其中步骤3所述的基于多项式拟合的特征点提取算法为：

步骤3.1，针对提取到的遥感图像轮廓边界上的每一个像素点，将其与周围的边界点进行三次多项式拟合，估计出拟合曲线的系数；

步骤3.2，计算拟合曲线的拟合误差以及该像素点处的曲率；

步骤3.3，抛弃拟合误差大的像素点，并在剩余的像素点中筛选出曲率大的点作为特征点；

步骤3.4，采用该像素点处曲线切线方向作为特征点的主方向；

其中步骤4中所述的改进的形状内容描述子为圆形模板，并赋予方向，该方向与特征点的主方向相同。

本发明与现有技术相比，具有以下优点：(1)本发明方法提出了一种新的基于多项式拟合的特征点提取算法，大大减少了计算量，使得经典的形状内容描述子在实际复杂的遥感图像中实现应用；(2)本发明方法利用多项式拟合的方法为每个特征点确定了主方向，改进了经典的形状内容描述子，使之具有旋转不变性，从而可以对不同地形的多源遥感图像进行配准，适用性得到了提高；(3)本发明方法通过由粗到精进行配准，在粗配准的基础上通过RANSAC算法去除误匹配，提高了配准精度。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明方法流程图。

图2是形状内容描述子圆形模板示意图，其中，图2(a)是传统的形状内容描述子圆形模板，图2(b)是本方法中改进的具有旋转不变性的形状内容描述子圆形模板。

图3是待配准的多源遥感图像及本发明方法获得的配准效果图。其中图3(a)(d)(g)为来自卫星探测器拍摄的卫星影像，分辨率为1024×1024，(b)(e)(h)为Google地图数据库中图像，分辨率为696×696。(c)(f)(i)为本发明方法获得的配准效果图。

图4是使用传统的SIFT算法得到的配准效果图。

具体实施方式

结合图1，本发明提出的基于多项式拟合的多源遥感图像形状配准方法，首先分别提取参考图像和待配准图像的边缘轮廓特征，并通过一种基于多项式拟合的特征点提取算法提取特征点并确定特征点的主方向，利用改进的形状内容描述子并最小化特征点之间的匹配代价来完成粗配准。最后去除误匹配，完成精配准。在此基础上根据匹配点对估计出参考图像和待配准图像之间的单应变换矩阵H，将待配准图像经变换矩阵还原，完成配准。

实现上述基于多项式拟合的多源遥感图像形状配准方法的具体实施步骤如下：

步骤1，输入两幅多源遥感图像作为参考图像和待配准图像；

步骤2，采用Canny算子对两幅图像提取边缘轮廓特征，并通过形态学变换消除密集区域无用的边缘特征，留下图像主轮廓；

步骤3，针对两幅图像提取到的轮廓上的像素点分别通过一种基于多项式拟合的特征点提取算法提取特征点并确定特征点的主方向；

步骤4，获取每一特征点的改进的形状内容描述子；

步骤5，基于像素点的改进的形状内容描述子，对参考图像中的每一个特征点，遍历待配准图像中特征点，获取与其匹配代价最小的待配准图像中特征点作为匹配点对；

步骤6，去除误匹配，完成精配准；

步骤7，根据匹配点对估计两图像间变换矩阵，将待配准图像经变换矩阵还原，完成配准。

具体地，在步骤3中，分别提取参考图像和待配准图像的主轮廓特征，针对提取到的遥感图像轮廓边界上的每一个像素点，将其与周围的边界点进行三次多项式拟合从而提取特征点并确定主方向。该操作共分为以下四个步骤：

步骤3.1，对轮廓上的像素点进行多项式拟合。

假定p_i是图像轮廓边界上的一个点，其坐标可以表示为（x_i,y_i)，提取图像轮廓边界上在点p_i左右的部分像素点，可以表示为p_i-m,…,p_i-1,p_i,p_i+1,…,p_i+n，像素点对应坐标可以表示为(x_i-m,y_i-m),…,(x_i,y_i),…,(x_i+n,y_i+n)。假定拟合曲线为三次多项式，故拟合曲线表示为：

y＝ax³+bx²+cx+d (1)

式中，a、b、c、d即表示拟合曲线的系数，x、y即为用于拟合的像素点坐标。

由于拟合的三次多项式曲线必须要保证经过当前点P_i。在此，本方法采用拉格朗日乘数法将其转化为多元函数求极值问题，从而来估计拟合曲线的系数。

以拟合曲线的系数a、b、c、d为自变量，估计拟合曲线系数问题可以转化为求函数z＝f(a,b,c,d)在附加条件下的极值问题。f(a,b,c,d)和具体公式为：

$f(a,b,c,d)=\underset{u=i-m}{\overset{i+n}{\Σ}}{({\mathrm{ax}}_u^3+{\mathrm{bx}}_u^2+{\mathrm{cx}}_u+d-y_u)}^2---\left(2\right)$

λ为系数，u从i-m到i+n，即选取的像素点拉格朗日函数公式为：

经计算，可以求解出系数a、b、c、d以及参数λ，表示如下：

[a,b,c,d,λ]＝(X^TX)^-1X^TY (5)

其中，

步骤3.2，进一步计算拟合曲线的拟合误差以及点P_i处的曲率来判断该点P_i能否作为特征点。拟合误差D(i)和点P_i处的曲率k(i)具体公式为：

$D\left(i\right)=\frac{1}{m+n+1}\underset{u=i-m}{\overset{i+n}{\Σ}}{({\mathrm{ax}}_u^3+{\mathrm{bx}}_u^2+{\mathrm{cx}}_u+d-y_u)}^2---\left(8\right)$

$k\left(i\right)=\frac{\left|y^{\′\′}\right|}{{\left(1+y^{\′2}\right)}^{\frac{3}{2}}}---\left(9\right)$

步骤3.3，拟合误差D(i)越小，即表示曲线与当前点P_i周围的真实轮廓边界一致性更高。对于轮廓形状特征来说，轮廓曲线弯曲度越大，即曲率越大，特征越明显。因此，为了提取特征更明显的特征点本方法中抛弃拟合误差D(i)大的点，并在剩余的拟合误差小的点中筛选出曲率大的点作为特征点。即设置一第一阈值，若拟合误差大于第一阈值，则抛弃该像素点。设置一第二阈值，若抛弃拟合误差大的点后，像素点处的曲率大于第二阈值，则该像素点为特征点。所述第一阈值、第二阈值根据想要获取的特征点个数来取值，特征点个数多，第一阈值大，第二阈值小。

步骤3.4，计算特征点主方向。本方法采用切线方向作为特征点的主方向，具体公式为：

$Q=3{\mathrm{ax}}_i^2+2{\mathrm{bx}}_i+c---\left(10\right)$

在步骤5中，对于提取出特征点计算改进的形状内容描述子并通过最小化匹配代价的方法实现粗配准。结合图2，图2(a)为经典的形状内容描述子圆形模板。假设图像轮廓边界上共有n个点，对于其中一点p_i，可以把它与边界上其它n-1采样点的结构关系用一个直方图h_i来表示：

h_i(k)＝#{q≠p_i；(q-p_i)∈bin(k)} (11)

其中，h_i(k)即定义为点p_i的形状内容描述子，为使得描述子更具鲁棒性，使用如图2(a)所示的对数—极坐标系，使形状内容描述子对距离参考点p_i越近的采样点变化越敏感。上式中bin即表示极坐标中的分段，k表示分段数，本方法中lgr方向上使用5bin，θ方向上使用12bin，形成一个60维特征直方图。

本方法中改进了原有的形状内容描述子圆形模板，在对特征点计算形状内容描述子的过程中，把特征点的主方向即切线方向作为圆形模板正轴，从而对圆形模板进行旋转，旋转的角度α具体公式如下：

α＝arctan Q (12)

为确定两幅待配准图像轮廓边界特征点之间的对应关系，在此采用C_ij≡C(p_i,q_j)表示这两特征点之间的匹配代价，具体公式为：

$C_{ij}=\frac{1}{2}\underset{k=1}{\overset{K}{\Σ}}\frac{{\[h_i\left(k\right)-h_j\left(k\right)\]}^2}{h_i\left(k\right)+h_j\left(k\right)}---\left(13\right)$

其中h_i(k)和h_j(k)分别表示点p_i和点q_j的K-bin归一化直方图。计算两幅图像中所有特征点对的匹配代价C_ij，通过最小化匹配代价来获得特征点一对一的对应关系。

步骤6中，对于粗配准得到的匹配点对，使用RANSAC算法去除误匹配，完成精配准。

本发明的效果可以通过以下仿真结果作进一步说明：

图3(a)(d)(g)为来自卫星探测器拍摄的卫星影像，分辨率为1024×1024，(b)(e)(h)为Google地图数据库中图像，分辨率为696×696。根据本发明提出基于多项式拟合的多源遥感图像形状配准方法得到的配准结果如(c)(f)(i)所示。图中地点从上之下分别为阿尔巴尼亚海岸，中心点纬度和经度分别为40°40'19.60"北，19°18'57.43"东；伊拉克Razazza湖，中心点纬度和经度分别为32°51'1.66"北，43°20'15.47"东和印度恒河(安拉阿巴德附近)，中心点纬度和经度分别为25°20'23.27"北，82°14'49.36"东。

为了说明本发明在配准精确度、适用性和计算速度上的优势，使用本发明方法和传统的SIFT配准方法分别对同样的多源遥感图像进行配准仿真实验。SIFT配准方法配准结果如图4所示。为进一步定量的评价遥感图像配准方法的性能，本文中分别通过互相关系数、均方根误差以及算法运算时间来综合评定配准效果。精度评价结果如下表所示，该表为本发明方法、SIFT配准方法进行仿真实验时，各方法的配准效果及运算时间对比。互相关系数在(0，1)范围内，值越大表明配准后的图像与目标图像之间相关性越大，配准效果越好。均方根误差通过来度量配准后图像与目标图像之间的距离离散度，均方根值越小，表明两者距离越小，即配准效果越好。本发明方法比SIFT算法适应性更广，得到的配准精度更高，计算速度比SIFT算法配准速度提高了2倍，能够应用于多源遥感图像的配准。

Claims (7)

1.一种基于多项式拟合的多源遥感图像形状配准方法，其特征在于，包括：

步骤1，输入两幅多源遥感图像作为参考图像和待配准图像，

步骤2，对两幅图像提取边缘轮廓特征，且消除密集区域无用的边缘特征，

步骤3，针对两幅图像提取到的轮廓上的像素点分别通过一种基于多项式拟合的特征点提取算法提取特征点并确定特征点的主方向，

步骤4，获取每一特征点的改进的形状内容描述子，

步骤5，基于像素点的改进的形状内容描述子，对参考图像中的每一个特征点，遍历待配准图像中特征点，获取与其匹配代价最小的待配准图像中特征点作为匹配点对，

步骤6，去除误匹配，完成精配准，

步骤7，根据匹配点对估计两图像间变换矩阵，将待配准图像经变换矩阵还原，完成配准；

其中步骤3所述的基于多项式拟合的特征点提取算法为：

步骤3.1，针对提取到的遥感图像轮廓边界上的每一个像素点，将其与周围的边界点进行三次多项式拟合，估计出拟合曲线的系数；

步骤3.2，计算拟合曲线的拟合误差以及该像素点处的曲率；

步骤3.3，抛弃拟合误差大的像素点，并在剩余的像素点中筛选出曲率大的点作为特征点；

步骤3.4，采用该像素点处曲线切线方向作为特征点的主方向；

其中步骤4中所述的改进的形状内容描述子为圆形模板，并赋予方向，该方向与特征点的主方向相同。

2.根据权利要求1所述的方法，其特征在于，所述步骤3.1中采用拉格朗日乘数法估计出拟合曲线的系数。

3.根据权利要求1所述的方法，其特征在于，所述步骤3.3中设置一第一阈值，若拟合误差大于第一阈值，则抛弃该像素点。

4.根据权利要求3所述的方法，其特征在于，所述步骤3.3中设置一第二阈值，若抛弃拟合误差大的点后，像素点处的曲率大于第二阈值，则该像素点为特征点。

5.根据权利要求1所述的方法，其特征在于，所述步骤4中的改进的形状内容描述子的旋转的角度α＝arctan Q，其中Q为特征点的主方向。

6.根据权利要求1所述的方法，其特征在于，所述步骤5中两特征点之间的匹配代价采用下式获得

$C_{ij}=\frac{1}{2}\underset{k=1}{\overset{K}{\Σ}}\frac{{\[h_i\left(k\right)-h_j\left(k\right)\]}^2}{h_i\left(k\right)+h_j\left(k\right)}$

其中i为参考图像中的特征点的索引值，j为待配准图像中的特征点的索引值，K为改进的形状内容描述子直方图中bin的个数，h_i(k)和h_j(k)分别表示特征点p_i和特征点q_j的K-bin归一化直方图。

7.根据权利要求1所述的方法，其特征在于，所述步骤6中采用RANSAC算法去除误匹配。