使用GDAL工具对卫星数据进行影像配准

一、        简介

本文将探讨使用GDAL来对卫星影像进行影像配准，依然以Orb-View3数据为例（选择北京市中心附近的影像为例）。其实按照文章中的方法，对任何影像都可以进行配准，不仅仅局限于卫星影像，只要能够提供图像行列号和真实的地理坐标就可以进行处理。

卫星图像来自免费的OrbView-3航天器，可以通过OrbView-3来了解更多信息。然而，在最原始的数据中，卫星图像是用没有地理位置的Tiff格式存储的（普通的Tiff格式）。关于Orb-View3的数据下载，可以参考博文http://blog.csdn.net/kupe87826/article/details/8078914和http://blog.csdn.net/kupe87826/article/details/8079001。

影像配准最常用的就是用来将多光谱影像和全色影像配准，然后进行融合得到高分辨率的多光谱影像。此外还有一个作用就是给一个没有投影的数据赋投影，群里面一直有童鞋问，怎么给一个没有投影的数据赋投影信息。看完这篇文章我觉得应该就知道了吧。

二、        软件和数据

准备要进行影像配准的卫星图像：我们选择的是北京市附近的一景Orb-View3的数据。参考影像使用GoogleMaps中的卫星影像。

所需要的软件：QGIS，GDAL工具集。选择控制点使用QGIS软件来进行辅助；之后使用GDAL的工具来对影像进行配准。

三、        准备正射校正必要的数据

下载的Orb-View3的数据的元数据信息如下。获取时间是2005年11月02日。现在基本上对于高分影像都需要进行正射校正才能达到使用的要求，但是对于北京市区，基本上都是平原，正射校正的意义不大，影像配准基本上也可以满足要求。（不过这里我们是要对配准进行一个说明，不管是哪里的数据，只是说明一下配准的步骤。）

表1Orb-View3卫星元数据信息表

Data Set Attribute	Attribute Value
Entity ID	87248
Acquisition Date	2005/12/17
Cloud Cover	0
Center Latitude	39°52'50.98"N
Center Longitude	116°22'54.31"E
NW Latitude	40°00'57.95"N
NW Longitude	116°22'14.62"E
NE Latitude	40°00'04.09"N
NE Longitude	116°27'42.31"E
SW Latitude	39°45'38.61"N
SW Longitude	116°18'03.99"E
SE Latitude	39°44'43.28"N
SE Longitude	116°23'36.33"E
Center Latitude dec	39.880828
Center Longitude dec	116.381753
NW Latitude dec	40.016097
NW Longitude dec	116.370728
NE Latitude dec	40.001136
NE Longitude dec	116.461754
SW Latitude dec	39.760724
SW Longitude dec	116.301108
SE Latitude dec	39.745355
SE Longitude dec	116.393424
Satellite	ORBVIEW-3
Sensor	Panchromatic
Sun Elevation	24.98
Sun Azimuth	164.14
Vendor	GeoEye Inc.
Vendor ID	3v051217p0001015651a520008201092m_001625916
Map Projection	GEOGRAPHIC
Processing Level	BASIC ENHANCED
Datum	WGS84
File Format	TIFF
Pixel Size X	1.0186931
Pixel Size Y	1.0101201
Date Entered	2011/10/11

 

下载的Orb-View3数据的范围如下图所示：

图1 Orb-View3数据的覆盖范围

接下来用QGIS打开Google的在线卫星地图，用将其缩放至北京地区。截图如下：

图2 QGIS打开的GoogleMaps数据

四、        影像配准

1、准备工作

要进行影像配准，首先要选取控制点，这里我们使用QGIS进行辅助选点工作。在上一小节中我们已经用QGIS打开了Google的卫星影像作为参考影像（Google的卫星数据的投影用的是Web墨卡托投影，不知道的童鞋请问问谷歌或者度娘）。然后使用QGIS中的Georeferencer插件来进行选取控制点。（Georeferencer插件是QGIS中用来选取控制点进行影像配准的工具，网址是：http://gis-lab.info/qa/qgis-georef-new-eng.html。）

点开插件之后，在插件中打开Orb-View3的数据作为待配准影像，如图3所示：

图3 选取控制点界面

2、选取控制点

接下来的主要工作就是选择控制点。选择控制点尽可能的均匀的分布在整幅图像上，然后我们将选择好的控制点保存成一个控制点文件，这个控制点文件是个文本格式，可以用记事本打开查看。下图4是选择的控制点列表，并且计算了中误差（这里的中误差有点大，可以移动控制点对其进行修改，我就不修改了）。

图4 选取的控制点列表

3、生成带有控制点的图像数据

由于GDALWarp工具不能直接加载控制点文件中的控制点，所以才有了这一步。这一步的目的就是把原来的图像和控制点信息都保存到一副新的图像中去。为接下来的第四步做准备工作。好了，废话不多说，这里我们用到的工具是gdal_translate（关于这个工具的作用和使用说明可以参考我之前的博文，或者直接看GDAL官方提供的帮助文档，地址是：http://www.gdal.org/gdal_translate.html。）。使用这个工具的命令行如下：

>gdal_translate -of GTiff-gcp 4280.35 11869 1.29566e+07 4.85182e+06 -gcp 3690.92 10956 1.29561e+074.85316e+06 -gcp 4462.46 10753.5 1.29571e+07 4.85321e+06 -gcp 4256.06 9726.071.29571e+07 4.85459e+06 -gcp 3390.84 9950.69 1.2956e+07 4.85454e+06 -gcp2338.91 2867.69 1.29565e+07 4.86395e+06 -gcp 5252.16 1409.01 1.29606e+074.86504e+06 -gcp 7577.26 199.105 1.29639e+07 4.86595e+06 -gcp 998.878 27260.61.29482e+07 4.83279e+06 -gcp 7211.95 28334.6 1.29559e+07 4.82966e+06 -gcp3907.23 24825.2 1.29526e+07 4.83515e+06 -gcp 165.189 21560 1.29487e+07 4.84042e+06-gcp 6352.08 23457.6 1.29561e+07 4.83624e+06 -gcp 4886.49 20770 1.2955e+074.84014e+06 -gcp 6701.2 17214.2 1.29582e+07 4.84425e+06 -gcp 2279.63 18098.31.29524e+07 4.84432e+06 -gcp 521.509 14313.2 1.29511e+07 4.8497e+06 -gcp4883.22 15076.1 1.29565e+07 4.84751e+06 -gcp 7595.89 11488.5 1.29609e+074.85141e+06 -gcp 203.457 8092.35 1.29524e+07 4.8578e+06 -gcp 7914.81 8033.131.29622e+07 4.85578e+06 -gcp 4791.71 6232.26 1.29587e+07 4.85895e+06 -gcp964.329 617.326 1.29554e+07 4.86722e+06 -gcp 3061.68 27991.1 1.29507e+074.83126e+06"F:/Data/Orb-View3/3V051217P0001015651A520008201092M_001625916/3v051217p0001015651a520008201092m_001625916.tif""F:/Data/Orb-View3/3v051217p0001015651a520008201092m_001625916_GCP.tif"
Input file size is 8016,28672
0...10...20...30...40...50...60...70...80...90...100- done.

图5是执行的截图，执行完成之后可以使用gdalinfo工具查看输出的图像中的信息，如果正常的话，应该会显示里面的控制点信息，这里我就不进行测试了，具体的可以参考之前使用GDAL正射的那篇博文。

图5 将控制点保存在文件中

4、使用GDALWarp进行图像配准

有了第三步的处理，下面可以直接使用GDALWarp工具来对图像进行配准了。配准的命令行如下：

> gdalwarp -r cubic -order2 -tr 1.000000 -1.000000"F:\Data\Orb-View3/3v051217p0001015651a520008201092m_001625916_GCP.tif""F:/Data/Orb-View3/3v051217p0001015651a520008201092m_001625916_georef.tif"

上面的命令中有个-order 2，意思是使用二次多项式进行校正，此外还可以用一次多项式和三次多项式，当然了，就是把3换成1或者3就行了，很简单吧。此外GDAL还提供了一种叫TPS（ThinPlane Spline）的校正方式，与Erdas的Rubber Sheeting算法比较类似。

处理的截图如图6所示。

图6 GDALWarp配准处理

五、        结论

至此，我们对影像的配准已经完成，接下来就是对其配准的精度进行验证。验证的方法是在QGIS中添加进来进行目视检查（当然这个是最直观的一种检查方式，缺点是不能定量分析，如果你不累的话，可以选些检查点进行定量分析一下）。下面是使用QGIS叠加显示的效果（图7），或者也可以用ArcMap或者Erdas的卷帘查看。

图7 图像配准完成后和基准影像叠加显示

      下面是把故宫部分放大进行比较，图8是googleMap的卫星图像，图9个校正后的图像，图10是把校正后的图像透明度设置为50%与GoogleMap的图像叠加显示的效果。

图8 GoogleMap故宫图像

图9 配准后的故宫图像

图10 使用50%透明度叠加显示

从目视效果来看，基本上大体上还是可以滴。由于选择的点的精度问题，导致在图10中的叠加显示有点模糊的效果，由于QGIS没有卷帘的功能，所以不能截一个卷帘显示的图，其实精度大体上还可以，如果调整控制点的精度的话，结果的精度会更高。如果对于山区的影像或者地形起伏以及卫星相机本身的畸变比较大的影像，使用常规的多项式可能达不到配准精度的要求，这时候就需要查找原因了。对于不同的算法原理，可以翻阅相关的参考书或者网络（基本上只要是本遥感图像处理的书肯定都会说几何多项式校正的）。

个人比较懒，没有使用Envi或者Erdas对这个数据进行处理对比。有兴趣的童鞋可以用同样的控制点使用Envi或者Erdas等对数据进行配准，然后与GDAL的结果比较一下，应该不会有太大的出入的。

再次说一句，GDAL提供的算法还是灰常强大滴！！^_^

六、        参考资料

[1]www.gdal.org/

[2] www.qgis.org

[3] http://pst.nst.pku.edu.cn/teaching/nuclear_medicine/DIPnotes/chapter10/chapter10.pdf

[4] http://www.mathworks.cn/products/image/description6.html

[5] http://blog.csdn.net/kupe87826/article/details/8078914

[6] http://www.mathworks.cn/cn/help/images/examples/registering-an-aerial-photo-to-an-orthophoto.html

[7] http://www.baike.com/wiki/%E5%9B%BE%E5%83%8F%E9%85%8D%E5%87%86