CN102436652A - 一种多源遥感图像自动配准方法 - Google Patents

Abstract

本发明公开了一种多源遥感图像自动配准方法，包括对基准图像和待配准图像重叠区进行裁剪和对待配准图像进行小波尺度变换，使其分辨率与基准图像一致的步骤；对基准图像进行切片并进行特征点检测，提取均匀分布的特征点并获得初始匹配同名点对的步骤；利用初始匹配同名点对统计求得待配准图像经纬度方向整体偏移量的步骤；根据经纬度方向整体偏移量修正初始匹配同名点对，并对修正后特征点进行双向互相关特征点匹配获取亚像素定位精度的精匹配同名点，最终重采样得到配准图像的步骤。采用本发明有效克服了HJ星多源遥感图像大视场、大跨度、异特性等特点，可快速检测大量分布均匀、稳定、可靠的高精度匹配同名点对，并进行基于密集三角网的微分纠正。

Description

一种多源遥感图像自动配准方法

技术领域

本发明涉及一种多源遥感图像自动配准方法，属于卫星图像处理技术领域，可用于HJ星红外、CCD和高光谱同类传感器和异类传感器图像的自动配准。

背景技术

图像配准是图像处理的基本任务之一，是对取自不同时间、不同传感器或不同视角的同一景物的两幅图像或多幅图像进行匹配、叠加的过程。已经被广泛应用于图像分析、变化检测和信息融合等领域。设参考图像和待配准图像分别用I₁(x₁，y₁)和I₂(x₂，y₂)表示，则图像配准在数学上可以表示为

I₁(x₁，y₁)＝g(I₂(T(x₂，y₂)))

其中，g(·)是一维灰度变换，T(·)是二维空间位置变换。由于在大多数情况下，配准后的两幅图像一般不要求灰度一致，而且配准也并不能通过一个或多个变换来达到两幅图像的完全等同，而只能做到某种程度上的近似。

HJ星多源遥感图像配准具有其自身的特殊性，主要表现在HJ星CCD相机和红外相机分别通过双相机拼接(二台±30度视场)和双面镜扫描(±29度扫描范围)实现大视场观测；高光谱成像仪幅宽50km，侧向可视范围±30度。HJ星多源遥感图像空间分辨率从30m、100m、150m到300m不等，跨度大。将低分辨率图像向与之分辨率差异较大的较高分辨率转换，要利用插值算法，增加带宽等方式，填充大量虚假信息，会引起图像失真，影响配准精度。

典型的图像配准算法大致可以分为基于像素的、基于灰度的、基于特征的和基于物理模型的方法。目前，对于多源遥感图像的配准，主要采用的是基于特征的配准方法。在基于特征的图像配准方法中，特征的选择与提取无疑对图像配准的质量有着很大的影响，其将极大地影响图像配准成功率。传统的图像配准技术要达到较高的精度，通常需要引入一定的人工干预，手工方式不但降低了***的效率，而且影响了***的应用性。目前，还没有应用于业务化***的HJ数据快速自配准方法。

如何实现特征点的自动提取与自动匹配、如何提高宽幅遥感影像匹配效率、如何采用考虑成像几何的自动配准更好满足大视场、大跨度、异特性特点，是本发明重点解决的问题。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种多源遥感图像自动配准方法，该方法有效克服了HJ星多源遥感图像大视场、大跨度、异特性等特点，可自动快速检测大量分布均匀、稳定、可靠的高精度匹配同名点对，并进行基于密集三角网的微分纠正，支持HJ星红外、CCD和高光谱同类传感器和异类传感器的多源遥感图像的配准。本发明解决了传统多源影像配准过程中需要人为干预的问题，在处理初始阶段通过预设参数就可实现全流程自动配准。

本发明的技术解决方案是：一种多源遥感图像自动配准方法，包括下列步骤：

1.根据遥感数据自带的经纬度信息，对基准图I₁和待配准图像I₂的重叠区进行裁剪，得到基准图裁剪结果影像I₁ ^*和待配准图裁剪结果影像I₂ ^*。

2.根据基准图和待配准图像分辨率信息R₁和R₂，判断分辨率是否一致，一致则进行下一步；不一致则对待配准图裁剪结果影像I₂w进行尺度变换：

(1)对待配准图裁剪结果影像I₂ ^*进行Mallat小波分解；得到一个低频图像信号以及不同分辨率上的高频图像信号。

(2)图像经过小波分解后，待配准图裁剪结果影像I₂ ^*被分解成沿水平、垂直和对角三个方向的一系列高频信号图像与一个低频图像，其中同一方向的各高频图像信号之间具有相似性。小波尺度变换就是利用这种相似性，进行高频外推插值；

(3)利用小波反变换重构出比原始图像分辨率更高的插值图像WD，本方法在高频外推插值时采用的插值方法为双三次插值。

3.对基准图裁剪结果影像I₁ ^*并进行FAST特征点检测，提取均匀分布的特征点；依据2级图像的地理信息对应关系在待校正图像上定位匹配点，得到初始匹配同名点对：

(1)将基准图裁剪结果影像I₁ ^*按照预设大小进行切片，得到图像片I_1i(i＝1，2，...，M，M是图像片总数)。一个线程一次对一个切片的数据I_1i进行FAST特征点检测，实现特征点检测阶段的并行处理：

(2)记基准图裁剪结果影像I₁ ^*检测的特征点为PixI_1i(i＝1，2，3，...，N)。GDAL遥感图像库的OpenGIS坐标转换方法提供了用来描绘坐标***(投影和基准面)以及坐标***相互之间转换的服务，本发明通过这种坐标转换结合地理定位信息在待配准裁剪图像的尺度变换图像WD上获取初始匹配同名点PixW_i(i＝1，2，3，...，N)，N为特征点总数目。

4.选取初始匹配同名点对{PixI_1i，PixW_i}(i＝1，2，3，...，N)的一个子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)。基于金字塔和相关系数进行特征点匹配(匹配搜索范围优于HJ影像***几何校正的精度)，一个线程一次匹配一个点，达到并行计算目的，经粗差剔除后获得更准确的同名点对，统计求得待配准图像经纬度方向整体偏移量：

(1)将基准图裁剪结果影像I₁ ^*和待配准裁剪图像的尺度变换图像WD建立3层高斯金字塔，图像金字塔从下到上对应I₁ ^*、I₁₁ ^*、I₁₂ ^*和WD、WD₁、WD₂，从下往上每层图像缩小1/3。

(2)第1层：

(2.1)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(2.2)在顶层影像I₁₂ ^*上以PixI_1j/9为中心获取127*127图像块，WD₂上以PixW_j/9为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)和峰值处相关系数C_max，当C_max＞0.5时，进行下一步；否则，回到步骤(2.1)计算下1个匹配点。

(2.3)在WD₂上以(p，q)为中心获取127*127图像块，顶层影像I₁₂ ^*上以PixI_1j/9为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j/9对应的点位置相等时，认为PixI_1j/9与(p，q)为1对匹配同名点；否则，回到步骤(2.1)计算下1个匹配点。

(2.4)在顶层影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，进行RANSAC(Random Sample Consensus)剔除，获取可靠匹配同名点对{LvI_2t，LvW_2t}(t＝1，2，...，T)。

(3)第2层：

(3.1)根据顶层第1层匹配同名点{LvI_2t，LvW_2t}(t＝1，2，...，T)，各点位置乘以金子塔放大倍数3，计算第2层待配准图像到基准图像的仿射变换关系；

(3.2)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(3.3)在第2层影像I₁₁ ^*上，根据放射关系计算PixI_1j/3在WD₁上对应的位置PixW^* _j，计算WD₁上PixW_j/3与PixW^* _j在x和y方向上的距离，若均小于15则进行下一步；否则，回到步骤(3.2)计算下1个匹配点。

(3.4)在第2层影像I₁₁ ^*上，以PixI_1j/3为中心获取127*127图像块，WD₁上以PixW^* _j为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)。

(3.5)在WD₁上以(p，q)为中心获取127*127图像块，第2层影像I₁₁ ^*上，以PixI_1j/3为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j/3对应的点位置相等时，认为PixI_1j/3与(p，q)为1对匹配同名点；否则，回到步骤(3.2)计算下1个匹配点。

(3.6)在第二影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，获取匹配同名点对{LvI_1t，LvW_1t}(t＝1，2，...，T)。

(4)第3层(最底层)：

(4.1)根据顶层第2层匹配同名点{LvI_1t，LvW_1t}(t＝1，2，...，T)，各点位置乘以金子塔放大倍数3，计算第3层待配准图像到基准图像的仿射变换关系；

(4.2)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(4.3)在第3层影像I₁ ^*上，根据放射关系计算PixI_1j在WD上对应的位置PixW^* _j，计算WD₁上PixW_j与PixW^* _j在x和y方向上的距离，若均小于15则进行下一步；否则，回到步骤(4.2)计算下1个匹配点。

(4.4)在第3层影像I₁ ^*上，以PixI_1j中心获取127*127图像块，WD上以PixW^* _j为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)。

(4.5)在WD上以(p，q)为中心获取127*127图像块，在第3层影像I₁ ^*上，以PixI_1j为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j对应的点位置相等时，认为PixI_1j与(p，q)为1对匹配同名点；否则，回到步骤(3.2)计算下1个匹配点。

(4.6)在第3影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，获取匹配同名点对{LvI_t，LvW_t}(t＝1，2，...，T)。

(5)统计经纬度方向整体偏移量Δx，Δy：由匹配过程可知，最底层匹配同名点{LvI_t，LvW_t}(t＝1，2，...，T)中的LvI_t(t＝1，2，...，T)为始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)中PixI_1j的一个子集，可得到LvI_t(t＝1，2，...，T)对应的初始匹配点PixW_t(t＝1，2，...，T)，计算LvW_t(t＝1，2，...，T)与PixW_t(t＝1，2，...，T)在经纬度方向(对应图像横轴x和纵轴y方向)的均值偏移量Δx和Δy。

5.根据经纬度方向整体偏移量Δx，Δy修正初始匹配同名点对{PixI_1i，PixW_i}(i＝1，2，3，...，N)中待配准图对应匹配点PixW_i(i＝1，2，3，...，N)，其x坐标减去Δx，y坐标减去Δy；得到修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)。

6.对进行双向互相关特征点匹配，匹配峰值处采用二元三点插值获取亚像素定位精度的精匹配同名点；同时对精匹配同名点对进行多项式拟合验证，剔除显著误匹配点：

(1)提取修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)的一个匹配同名点{PixI_1i，PixW^* _i}

(2)在基准图裁剪结果影像I₁ ^*上以PixI_1j为中心获取127*127图像块，待配准裁剪图像的尺度变换图像WD上以PixW^* _i为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)和峰值处相关系数C_max，当C_max＞0.5时，进行下一步；否则，回到步骤(1)计算下1个匹配点。

(3)在WD上以(p，q)为中心获取127*127图像块，基准图裁剪结果影像I₁ ^*上以PixI_1j/9为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j对应的点位置相等时，认为PixI_1j与(p，q)为1对匹配同名点，进行下一步；否则，回到步骤(1)计算下1个匹配点。

(4)对相关系数矩阵C峰值对应位置(p，q)周围3*3邻域内采用二元三点插值算法对相关系数进行插值，插值间隔选取1/20，获取插值后3*3邻域相关系数矩阵Z，获取Z的对应的峰值位置(p^*，q^*)，此时PixI_1j与(p^*，q^*)为一对高精度匹配同名点。

(5)回到步骤(1)，遍历修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)的每一个同名点，得到亚像素匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，M)。

(6)对精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，M)进行RANSAC剔除，获取可靠精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，D)。

7.结合获取的大量精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，D)，构建密集Delaunay三角网，对各小面元分别进行纠正，得到高精度配准图像。

本发明与现有技术相比的有益效果在于：

(1)本发明对HJ星多源遥感图像特性进行了深入分析，设计了包括尺度变换、基于成像几何的初始匹配同名点对获取、经纬度整体偏移消除、同名点对精匹配、密集三角网纠正等步骤的全自动配准流程；在整个配准过程中考虑了HJ星配准数据量和计算性能，在特征点检测和匹配阶段均采用并行设计，提高了配准效率，满足业务化运行需求。

(2)HJ星多源遥感图像配准误差主要是经纬度方向误差，本发明结合成像几何对配准误差的影响，针对经纬度方向的整体误差采用了金字塔与互相关双向匹配的方法，提高了初始匹配精度，减小了后续匹配计算量。

(3)本发明针对同名点精匹配，采用互相关双向匹配与二元三点插值相结合的方法，将匹配精度提高到了子像元级。

(4)本发明结合HJ星几何纠正数据配准误差主要是平移误差的特点，具有同名点对分布均匀、稳定、可靠，定位精度高等优势，可检测大量匹配同名点，支持HJ星红外、CCD和高光谱同类传感器和异类传感器图像的配准，大量的试验结果表明，配准精度可达到平地优于2个像素，山地优于5个像素。

附图说明

图1为本发明的整体流程图；

图2为本发明基于金字塔和相关系数定位经纬度方向偏移量流程图；

图3为本发明特征点精匹配流程图。

具体实施方式

下面结合附图1、图2和图3对本发明的具体实施方式进行进一步的详细描述：

1.根据遥感数据自带的经纬度信息，对基准图I₁和待配准图像I₂的重叠区进行裁剪，得到基准图裁剪结果影像I₁ ^*和待配准图裁剪结果影像I₂ ^*。

2.根据基准图和待配准图像分辨率信息R₁和R₂，判断分辨率是否一致，一致则进行下一步；不一致则对待配准图裁剪结果影像I₂ ^*进行尺度变换：

(1)小波分解过程表述为：令待配准图裁剪结果影像I₂ ^*表示大小为M×N的原始图像，l(i)表示相对于分析小波的低通滤波器系数，i＝0，1，2，…N_l-1，N_l表示滤波器L的支撑长度；h(i)表示相对于分析小波的高通滤波器系数，i＝0，1，2，…N_h-1，N_h表示滤波器H的支撑长度，则小波分解为：

$I_L(x,y)=\frac{1}{N_l}\underset{i=0}{\overset{N_l-1}{\mathrm{\Σ}}}l\left(i\right)I_2^{*}((2x+i)\mathrm{mod}M,y)$

$I_H(x,y)=\frac{1}{N_h}\underset{j=0}{\overset{N_h-1}{\mathrm{\Σ}}}l\left(j\right)I_2^{*}((2x+j)\mathrm{mod}M,y)$

$x=\mathrm{0,1},\·\·\·,\frac{M}{2}-1;$ y＝0，1，…，N-1

$I_{\mathrm{LL}}(x,y)=\frac{1}{N_l}\underset{i=0}{\overset{N_l-1}{\mathrm{\Σ}}}l\left(i\right)I_L(x,(2x+i)\mathrm{mod}N)$

$I_{\mathrm{LH}}(x,y)=\frac{1}{N_h}\underset{j=0}{\overset{N_h-1}{\mathrm{\Σ}}}h\left(j\right)I_L(x,(2x+j)\mathrm{mod}N)$

$I_{\mathrm{HL}}(x,y)=\frac{1}{N_l}\underset{i=0}{\overset{N_l-1}{\mathrm{\Σ}}}l\left(i\right)I_H(x,(2x+i)\mathrm{mod}N)$

$I_{\mathrm{HH}}(x,y)=\frac{1}{N_h}\underset{j=0}{\overset{N_h-1}{\mathrm{\Σ}}}h\left(j\right)I_H(x,(2x+j)\mathrm{mod}N)$

$x=\mathrm{0,1},\·\·\·,\frac{M}{2}-1;$ $y=\mathrm{0,1},\·\·\·,\frac{N}{2}-1$

对逼近子图重复此过程，直到确定的分解水平，这样就能实现图像从高分辨率到低分辨率的小波分解。低频信号能反映出原始信号的基本波形，它是原始信号在低一级分辨率上的表示，而另一方面，高频信号则包含着由高分辨率变换到低分辨率所丢失的细节信息。如果输出的低频信号被当作新的信号输入到同一个滤波器组，再次进行离散小波变换，并且如此循环若干次，最终会得到一个低频信号以及不同分辨率上的高频细节信息。

(2)图像经过小波分解后，待配准图裁剪结果影像I₂ ^*被分解成沿水平、垂直和对角三个方向的一系列高频信号图像与一个低频图像，其中同一方向的各高频图像信号之间具有相似性。小波尺度变换就是利用这种相似性，进行高频外推插值，然后利用小波反变换重构出比原始图像分辨率更高的插值图像WD。

不妨设某图像为一低通滤波后得到的图像IL，先通过小波变换把IL分解成IH1，IV1，ID1以及低频图像IL1，然后采用相似变换外推得到高频细节IH，IV，ID，再用小波反变换得到比原始图像分辨率更高的图像I。

I＝DWT^-1(IL，IH，IV，ID)

本方法在高频外推插值时采用的插值方法为双三次插值。

3.对基准图裁剪结果影像I₁ ^*并进行FAST特征点检测，提取均匀分布的特征点；依据2级图像的地理信息对应关系在待校正图像上定位匹配点，得到初始匹配同名点对：

(1)将基准图裁剪结果影像I₁ ^*按照预设大小进行切片，得到图像片I_1i(i＝1，2，...，M，M是图像片总数)。一个线程一次对一个切片的数据I_1i进行FAST特征点检测，实现特征点检测阶段的并行处理：

对切片进行FAST特征点检测，表述为：检查待检测点c周围的圆，寻找其中最长的圆弧，如果圆弧中所有的点的灰度值都大于I(c)或都小于I(c)，则被判定为角点，其中I(c)为c的灰度值。在离散情况下，圆大小为3×3区域，弧长为离散点数目。灰度判定公式为：

|I(c′)-I(c)|≤t

其中I(c′)为3×3圆弧上的点的灰度值，t为阈值。只要满足上述公式的连续弧长的数目大于等于9，则该中心点为角点。

在对切片I_1i计算FAST特征点提取时，为了得到稀疏且分布均匀的特征点，对图像切片I_1i进行20*20个像素的细粒度划分，计算每个20*20像素粒度里面的特征点响应R_j(j＝1，2，...，400)。取角点响应的最大值R_max＝max(R_j)，如果R_max＞0，则标记该角点为特征点。通过上述计算，可以保证基准图上成像区每20*20个像素区域内包含一个特征点。

(2)记基准图裁剪结果影像I₁ ^*检测的特征点为PixI_1i(i＝1，2，3，...，N)。GDAL遥感图像库的OpenGIS坐标转换方法提供了用来描绘坐标***(投影和基准面)以及坐标***相互之间转换的服务，本发明通过这种坐标转换结合地理定位信息在待配准裁剪图像的尺度变换图像WD上获取初始匹配同名点PixW_i(i＝1，2，3，...，N)，N为特征点总数目。具体过程表述如下：

提取PixI_1i(i＝1，2，3，...，N)的一个点PixI_1i，根据其在图像上的坐标以及基准图的地理信息计算该点经纬度Geo_i；

根据经纬度Geo_i结合尺度变换图像WD的地理信息计算对应的图像坐标PixW_i；

遍历PixI_1i(i＝1，2，3，...，N)的每一个点，得到对应的初始匹配点PixW_i(i＝1，2，3，...，N)。

4.选取初始匹配同名点对{PixI_1i，PixW_i}(i＝1，2，3，...，N)的一个子集{PixI_1j，PixW_j}(j＝1，2，3，...，M＜N)，本发明中M＝50。基于金字塔和相关系数进行特征点匹配(匹配搜索范围优于HJ影像***几何校正的精度)，一个线程一次匹配一个点，达到并行计算目的，经粗差剔除后获得更准确的同名点对，统计求得待配准图像经纬度方向整体偏移量：

(1)将基准图裁剪结果影像I₁ ^*和待配准裁剪图像的尺度变换图像WD建立3层高斯金字塔，图像金字塔从下到上对应I₁ ^*、I₁₁ ^*、I₁₂ ^*和WD、WD₁、WD₂，从下往上每层图像缩小1/3。

(2)第1层：

(2.1)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(2.2)在顶层影像I₁₂ ^*上以PixI_1j/9为中心获取127*127图像块，WD₂上以PixW_j/9为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)和峰值处相关系数C_max，当C_max＞0.5时，进行下一步；否则，回到步骤(2.1)计算下1个匹配点。

相关系数匹配表述为：它是基于灰度相关匹配算法中最常用的算法，经过去均值和归一化处理后的相关系数对图像对比度、明暗变得不敏感，使得归一化互相关系数更可靠、适应性更强。归一化互相关系数的定义如下：

$C(A,B)=\frac{\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}\left[(A_{i,j}-\stackrel{\‾}{A})(B_{i,j}-\stackrel{\‾}{B})\right]}{\sqrt{\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{(A_{i,j}-\stackrel{\‾}{A})}^2*\underset{i}{\mathrm{\Σ}}\underset{j}{\mathrm{\Σ}}{(B_{i,j}-\stackrel{\‾}{B})}^2}}$

其中，A，B分别代表大小为N×N的参考子图和待配子图。

分别为二维矩阵A，B的均值。

由归一化互相关的定义可知，相关系数C(A，B)最大值时对应的位置即为两幅图像的匹配位置。

(2.3)在WD₂上以(p，q)为中心获取127*127图像块，顶层影像I₁₂ ^*上以PixI_1j/9为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j/9对应的点位置相等时，认为PixI_1j/9与(p，q)为1对匹配同名点；否则，回到步骤(2.1)计算下1个匹配点。

(2.4)在顶层影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，进行RANSAC(Random Sample Consensus)剔除，获取可靠匹配同名点对{LvI_2t，LvW_2t}(t＝1，2，...，T)。

RANSAC剔除表述为：迭代地在输入数据{PixI_1j，PixW_j}(j＝1，2，3，...，M)中采样所谓的最小点集(3对匹配同名点)。并利用每次采样所得到的最小点集(3对匹配同名点)估计出所要确定的参数(仿射变换参数)，同时根据阈值则来判别输入数据中哪些是与该参数相一致的，即内点，哪些是不一致的，即外点。如此迭代一定次数后，将对应输入数据中内点比例最高的所估计参数值以及所筛选出来的内点作为RANSAC最后解。将此解作为其它方法的初始值进一步优化计算，从而得到最后的估计参数值。抽样次数N的计算公式如下：

$N=\frac{\ln (1-s)}{\ln (1-(1-{\mathrm{\σ}}^R))}$

其中σ为预期的原始数据错误率，s代表至少有一个最小子集包含所有内点的概率。

(3)第2层：

(3.1)根据顶层第1层匹配同名点{LvI_2t，LvW_2t}(t＝1，2，...，T)，各点位置乘以金子塔放大倍数3，计算第2层待配准图像到基准图像的仿射变换关系；

(3.2)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(3.3)在第2层影像I₁₁ ^*上，根据放射关系计算PixI_1j/3在WD₁上对应的位置PixW^* _j，计算WD₁上PixW_j/3与PixW^* _j在x和y方向上的距离，若均小于15则进行下一步；否则，回到步骤(3.2)计算下1个匹配点。

(3.4)在第2层影像I₁₁ ^*上，以PixI_1j/3为中心获取127*127图像块，WD₁上以PixW^* _j为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)。

(3.5)在WD₁上以(p，q)为中心获取127*127图像块，第2层影像I₁₁ ^*上，以PixI_1j/3为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j/3对应的点位置相等时，认为PixI_1j/3与(p，q)为1对匹配同名点；否则，回到步骤(3.2)计算下1个匹配点。

(3.6)在第二影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，获取匹配同名点对{LvI_1t，LvW_1t}(t＝1，2，...，T)。

(4)第3层(最底层)：

(4.1)根据顶层第2层匹配同名点{LvI_1t，LvW_1t}(t＝1，2，...，T)，各点位置乘以金子塔放大倍数3，计算第3层待配准图像到基准图像的仿射变换关系；

(4.2)提取初始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的一个初始匹配同名点{PixI_1j，PixW_j}。

(4.3)在第3层影像I₁ ^*上，根据放射关系计算PixI_1j在WD上对应的位置PixW^* _j，计算WD₁上PixW_j与PixW^* _j在x和y方向上的距离，若均小于15则进行下一步；否则，回到步骤(4.2)计算下1个匹配点。

(4.4)在第3层影像I₁ ^*上，以PixI_1j中心获取127*127图像块，WD上以PixW^* _j为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)。

(4.5)在WD上以(p，q)为中心获取127*127图像块，在第3层影像I₁ ^*上，以PixI_1j为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j对应的点位置相等时，认为PixI_1j与(P，q)为1对匹配同名点；否则，回到步骤(3.2)计算下1个匹配点。

(4.6)在第3影像上遍历匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)的每一个点，获取匹配同名点对{LvI_t，LvW_t}(t＝1，2，...，T)。

(5)统计经纬度方向整体偏移量Δx，Δy：由匹配过程可知，最底层匹配同名点{LvI_t，LvW_t}(t＝1，2，...，T)中的LvI_t(t＝1，2，...，T)为始匹配同名点子集{PixI_1j，PixW_j}(j＝1，2，3，...，M)中PixI_1j的一个子集，可得到LvI_t(t＝1，2，...，T)对应的初始匹配点PixW_t(t＝1，2，...，T)，计算LvW_t(t＝1，2，...，T)与PixW_t(t＝1，2，...，T)在经纬度方向(对应图像横轴x和纵轴y方向)的均值偏移量Δx和Δy。

5.根据经纬度方向整体偏移量Δx，Δy修正初始匹配同名点对{PixI_1i，PixW_i}(i＝1，2，3，...，N)中待配准图对应匹配点PixW_i(i＝1，2，3，...，N)，其x坐标减去Δx，y坐标减去Δy；得到修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)。

6.对进行双向互相关特征点匹配，匹配峰值处采用二元三点插值获取亚像素定位精度的精匹配同名点；同时对精匹配同名点对进行多项式拟合验证，剔除显著误匹配点：

(1)提取修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)的一个匹配同名点{PixI_1i，PixW^* _i}

(2)在基准图裁剪结果影像I₁ ^*上以PixI_1j为中心获取127*127图像块，待配准裁剪图像的尺度变换图像WD上以PixW^* _i为中心在周围5个像素范围滑窗，进行相关匹配，获取相关系数矩阵C及其峰值对应位置(p，q)和峰值处相关系数C_max，当C_max＞0.5时，进行下一步；否则，回到步骤(1)计算下1个匹配点。

(3)在WD上以(p，q)为中心获取127*127图像块，基准图裁剪结果影像I₁ ^*上以PixI_1j/9为中心在其周围滑窗，进行相关匹配，获取相关系数矩阵及其峰值对应位置(m，n)，当(m，n)与PixI_1j对应的点位置相等时，认为PixI_1j与(p，q)为1对匹配同名点，进行下一步；否则，回到步骤(1)计算下1个匹配点。

(4)对相关系数矩阵C峰值对应位置(p，q)周围3*3邻域内采用二元三点插值算法对相关系数进行插值，插值间隔选取1/20，获取插值后3*3邻域相关系数矩阵Z，获取Z的对应的峰值位置(p^*，q^*)，此时PixI_1j与(P^*，q^*)为一对高精度匹配同名点。

(5)回到步骤(1)，遍历修正后匹配同名点对{PixI_1i，PixW^* _i}(i＝1，2，3，...，N)的每一个同名点，得到亚像素匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，M)。

(6)对精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，M)进行RANSAC剔除，获取可靠精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，D)。

7.结合获取的大量精匹配同名点对{PixID_1i，PixWD_i}(i＝1，2，3，...，D)，构建密集Delaunay三角网，对各小面元分别进行纠正，得到高精度配准图像。

本发明未详细阐述的部分属于本领域公知技术。

Claims (5)

1.一种多源遥感图像自动配准方法，其特征在于该方法包含以下步骤：

A.对基准图像和待配准图像的经纬度重叠区进行裁剪，若基准图裁剪结果影像和待配准裁剪结果影像的分辨率不一致，则对待配准裁剪结果影像进行Mallat小波尺度变换；

B.对基准图裁剪结果影像进行FAST特征点检测，提取均匀分布的特征点；依据基准图裁剪结果影像和待配准裁剪结果影像的地理信息对应关系在待配准裁剪结果影像上定位特征点的匹配点，得到初始匹配同名点对；

C.提取初始匹配同名点对的一个子集，在子集内基于金字塔和相关系数进行特征点匹配，剔除显著错误点，统计计算得到经纬度方向整体偏移量；

D.通过经纬度方向整体偏移量对初始匹配同名点对进行修正，剔除平移误差；对每个匹配同名点进行双向互相关匹配并在峰值处进行二元三点插值，获取亚像素级匹配精度的匹配同名点对，利用亚像素级匹配精度的匹配同名点对获得配准图像。

2.如权利要求1中所述的一种多源遥感图像自动配准方法，其特征在于：在所述步骤A中所述的对待配准裁剪结果影像进行Mallat小波尺度变换的方法为：

A1.对原始待配准图裁剪结果影像进行Mallat小波分解，得到一个低频图像信号以及不同分辨率上的高频图像信号；

A2.根据其中同一方向的各高频图像信号之间的相似性，进行高频外推插值；

A3.对插值后的低频图像信号以及不同分辨率上的高频图像信号利用小波反变换重构出比原始待配准图裁剪结果影像分辨率更高的待配准图裁剪结果影像。

3.如权利要求1中所述的一种多源遥感图像自动配准方法，其特征在于：在步骤B中获得初始匹配同名点对的步骤中：

利用GDAL遥感图像库的OpenGIS坐标转换方法以及坐标***相互之间转换的服务，通过地理定位信息确定待配准裁剪图像的尺度变换图像上对应的初始匹配同名点。

4.如权利要求1中所述的一种多源遥感图像自动配准方法，其特征在于：在所述步骤C中获得经纬度方向整体偏移量的步骤中：

对基准图像和待配准图像裁剪结果图像建立3层高斯金字塔，初始匹配同名点的子集进行变换后先在金子塔顶层进行双向互相关匹配，匹配关系逐步映射到下层金字塔。

5.如权利要求1中所述的一种多源遥感图像自动配准方法，其特征在于：在所述步骤D中获取亚像素级匹配精度的匹配同名点对的步骤为：

D1.通过经纬度方向整体偏移量对初始匹配同名点对进行修正，剔除平移误差后，对每个匹配同名点进行相关匹配；

D2.相关系数峰值处采用采用二元三点插值获取亚像素定位精度。所有修正后匹配同名点遍历完后，对精匹配同名点对进行RANSAC剔除显著误匹配点，获取匹配同名点对。

Images