Harbin Institute of Technology
多时相遥感图像配准
实验报告
课程名称：遥感信息处理导论
院系：电子与信息工程学院
姓名：
学号：
授课教师：
哈尔滨工业大学
1. 实验目的
对于通过遥感平台获取的遥感数字影像，由于扫描过程中受地球曲率、地球自转、平台姿态、扫描方式等因素的影响，所获取的图像往往会产生不同程度的几何误差，这些误差如若不经处理，对图像的后续应用如分类、目标检测等会产生很大影响。
本实验的目的就是利用多项式映射和重采样等方法对两幅不同时间采集的遥感图像进行几何校正处理，从而方便从配准图像中找到不同时相内发生变化的地物目标信息。
2. 实验原理
图像的配准过程主要可以分为三个过程：控制点对选取（自动或手动）、坐标映射函数拟合、映射后像素重采样。
1) 控制点对选取
对于图像匹配过程，首先要通过控制点对建立两幅图像间的坐标联系。
控制点的选取有可分为自动和手动两种方式。
对于自动选取，可以采用SIFT 等算法实现；对于手动选取，需要实验人员通过目测观察的方法找到不同时相的遥感图像中的相同目标，这些控制点可以是建筑物的顶点、道路等的交叉点等或其他较容易分辨的位置。
选出控制点后将这些点的坐标信息记录下来，值得注意的是所选控制点的数量、分布情况以及精度会直接影响配准结果的精度和质量。
2) 图像几何校正
图像的配准实质上就是通过图像的几何校正，将产生几何失真的图像转换为标准的数字图像。
采用的方法就是通过一定的映射函数将原图像的像素坐标转换为标准图像中的坐标
(,)u f x y = (,)v g x y =
其中(,)x y 表示原图像中的像素坐标值，(,)u v 表示参考图像中的像素坐标值。
映射函数的选择可以有多种形式，一般较为简单常用的是多项式函数，以二次多项式函数为例，映射函数的具体形式为
22
01234522
012345u a a x a y a xy a x a y v b b x b y b xy b x b y
⎧=+++++⎨=+++++⎩ 其中，05,
,a a ，05,
,b b 分别表示二次多项式横纵坐标映射函数中对应项的系
数，通过选定控制点，可以将这些系数项求解出来。
假设在两幅图像中总共选择了n 对控制点，根据n 对控制点的坐标
(,)~(,)i i i i x y u v 可以得到如下两个方程组，
221011213114151222012223224252
22012345n n n n n n n u a a x a y a x y a x a y u a a x a y a x y a x a y u a a x a y a x y a x a y ⎧=+++++⎪=+++++⎪⎨
⎪
⎪=+++++⎩ 22
1011213114151222012223224252
22012345n n n n n n n
v b b x b y b x y b x b y v b b x b y b x y b x b y v b b x b y b x y b x b y ⎧=+++++⎪=+++++⎪⎨
⎪
⎪=+++++⎩ 如果另()1,
,T
n U u u =，()1,,T
n V v v =，05(,
,)T a a α=，05(,,)T b b β=，
2211
11112222
2222
22111n n
n n n n x y x y x y x y x y x y A x y x y x y ⎡⎤
⎢
⎥⎢
⎥=⎢⎥⎢⎥⎢⎥⎣⎦
那么上述两个方程组的向量和矩阵形式便可以写作
U A V A αβ
==
对于二次多项式来讲，每个方程组都有6个未知数，所以至少需要6对控制点才可以求解出各次项的系数，在实际实验中，为了使得到的模型更加准确，一般会选择多于6对的控制点进行计算，求解出最小二乘意义下最优系数α和β。
此处的误差根据映射后的待配准像素坐标和参考像素坐标间差的绝对值求平均来计算，横纵坐标的误差计算公式分别为
11
_error x A U n α=
- 11
_y error A V n
β=-
其中，1•表示向量或矩阵中所有元素的绝对值之和。
在计算出多项式系数后便得到了图像的几何校正多项式模型，根据待配准图像中每一像素的坐标值(,)i i x y 便可以计算出其在参考图像中的位置(,)i i u v 。
3) 图像重采样/栅格化
经过二次多项式映射后的各像素坐标并不能保证在参考图像的坐标系中处于整数点位置，它们在二维平面中的分布可能是散乱的无规律点，所以需要经过一个重采样/栅格化过程，根据这些无规律分布点的像素值将其对应的整数点位置的像素值计算出来。
一般常用的方法有最近邻插值、双线性插值、三次卷积插值。
其原理示意图如下所示：
g(u,v)f(x,y)
g(u,v)=f(x,y)
图1(a) 最近邻插值
图1(b) 双线性插值
图1(c) 三次多项式插值
3.实验步骤
本次遥感图像配准实验是在MATLAB平台下实现的，大致步骤如下：
1) 图像数据的导入和控制点选择
首先将两幅不同时间获取的遥感影像数据600X900_20121124.jpg和
600X900_20140705.jpg导入到MATLAB的工作区中，可以看到两幅图像分别存储在矩阵变量x600X900_20121124和x600X900_20140705中。
然后用cpslect函数实现控制点的手动选取，命令为
cpslect(‘x600X900_20121124’,’x600X900_20140705’);
其中前者为待配准图像，后者为参考图像。
在控制点选择窗口中，通过目测观察，寻找两幅图像中较为明显的相同的地物目标，并将其贴上相同的标签号码，本实验共选取了27对控制点，并使其在图像范围内分布较为均匀。
下图2即为控制点选择结果。
图2 控制点对选择结果图
在完成控制点对的选择后，在控制点选择窗口中选择file→Export Points to Workspace 命令，将选取的控制点对信息保存，然后在MATLAB 的工作区中就会出现两个变量，其中，movingPoints 中存储的是待配准图像中的控制点坐标，fixedPoints 中存储的是参考图像中的控制点坐标。
2) 几何校正模型参数计算及坐标映射
以选择二次多项式校正模型为例，根据控制点坐标信息(),i i x y 计算出待配准控制点信息矩阵A ，然后结合参考控制点坐标(),i i u v 计算出二次多项式系数的最小二乘解α、β，从而确定坐标映射函数的具体形式。
下面两个表格分别给出了一次和二次多项式映射函数的系数。
对于待配准图像中的每一像素坐标(,)i i x y ，通过多项式映射函数计算出其对应的参考图像坐标系中的坐标(,)i i u v ，
表1 一次多项式映射系数
i 0 1 2
α41.5969 0.9998 9.00e-04
β10.5613 4.22e-04 0.9987
表2 二次多项式映射系数
i 0 1 2 3 4 5
α41.3574 1.0000 0.0024 -1.79e-06 4.82e-07 -1.46e-06
β11.0188 -3.17e-04 0.9954 3.32e-06 -4.32e-07 2.05e-06
3) 图像重采样/栅格化
通过多项式映射后的像素坐标大多情况下不是整数，需要通过重采样过程得到映射后的栅格点像素值。
本实验运用griddata函数实现的重采样/栅格化功能，采用的是双线性插值的方法，得到的结果如下图
图3 映射结果经重采样得到的栅格图像
4.实验结果与分析
下表为采用一次和二次多项式映射后的平均坐标误差。
表3 映射坐标误差
多项式次数一次二次
error x0.2869 0.2699
_
error y0.4132 0.3899
_
下图是以二次多项式函数为映射函数采用双线性插值的方法的配准效果图
图4 二次多项式映射配准结果图
由表1和表2可以看出，横坐标方向的几何畸变主要是由由固定偏移量0a 和
x 坐标引起的，y 坐标的影响比较小；纵坐标方向的几何畸变主要是由由固定偏
移量0b 和y 坐标引起的，x 坐标的影响比较小。
对比表1和表2，并结合表3的坐标映射误差来看，多项式次数的提高对图像的几何校正并没有太大的改善效果，由此可以推测造成这两幅图的几何畸变的主要原因是线性畸变（平移、旋转）。
5. 程序源代码
%加载图像数据
%cpselect(x600X900_20121112,x600X900_20140705); %以前者为待配准图
参考图像
待配准图像
配准后图像
配准效果。