几何校正.docx

实习报告第10章 图像几何校正学生姓名：学生学号：学生年级：学生专业： 学生院系：任课教师：提交时间：资源与环境遥感模型实验实习报告图像几何校正1. 目的通过实习操作，掌握遥感图像几何校正的基本原理和方法，理解遥感图像几何畸变的来源、几何校正的意义以及坐标系在几何校正中的作用。2. 内容2.1 图像几何校正概述几何变形表现为影像上的像元相对于地面目标的实际位置发生的挤压、扭曲、拉伸和偏移等。其来源主要有：引起图像的几何变形一般分为两类：系统性和非系统性。系统性一般是由传感器本身引起的，有规律可循和课预测性，可以用传感器模型来校正，卫星地面接收站一般已完成这项校正；非系统性几何变形是不规律的，它可以是传感器平台本身的高度、姿态等不稳定性，也可以是地球曲率及空气折射的变化以及地形的变化等引起的。我们常说的几何校正就是校正这些非系统的几何变形。几何校正是利用地面控制点和几何校正数学模型来校正非系统性因素引起的误差，同时也是将图像投影到平面上，使其符合底图投影系统的过程；由于校正过程会将坐标系统赋予图像数据，所以此过程包括了地理编码。2.2 图像几何校正方法在ENVI4.8软件中，针对不同的数据源和辅助数据，有以下几种方法(1) 利用卫星自带的地理定位文件进行几何校正对于重返周期短，空间分辨率低的卫星数据，如AVHRR、MODIS等，地面控制点选择比较困难。这是就可以利用卫星传感器自带的定位文件进行几何校正，其校正的精度主要受定位文件的影响。选择主菜单MapGeoreference 传感器名称，启动这种校正方法。(2) Image to Image 几何校正以一幅已经经过几何校正过的的遥感影像为基准影像，在需要校正的影像和基准影像上选取相同的地物点，使相同的地物出现在校正后的图像上相同的位置处。这是大多数几何校正所采用的方法。选择主菜单MapRegistrationSelect GCPs：Image to Image，启动这种校正方法。(3) Image to Map 几何校正通过地面控制点对遥感影像的进行几何校正，控制点可以从键盘输入，矢量文件中获取或栅格文件中获取，如地图的校正。选择主菜单MapRegistrationSelect GCPs：Image to Map，启动这种校正方法。(4) Image to Image 自动图像配准根据像元灰度值或者地物特征自动寻找两幅图像的同名点，完成图像的几何校正。如，当同一地区的两幅影像由于各自的校正误差的影响，使得图上的相同地物的不重叠时，可以使用此方法进行调整。选择主菜单MapRegistrationAutomatic Registrations：Image to Image，启动这种校正方法。2.3 控制点选择方式ENVI提供以下控制点的选择方式（1）从栅格图像中选择如果有同一地区的已校正的图像，地形图等栅格数据，可以从中选择控制点，对应的控制点选择模式为Image to Image。（2）从矢量数据中选择如果有同一地区的已校正的矢量数据，可以从中选择控制点，对应的控制点选择模式为Image to Map。（3）从文本文件中导入如果有GPS测量数据，或者其它途径获得控制点坐标，可以导入到图像中，对应的控制点选择模式为Image to Map。（4）键盘输入如果只有控制点目标的坐标信息，则只好通过键盘输入坐标数据并在图像上找到对应点，对应的控制点选择模式为Image to Map。2.4 几何校正原理图像的几何校正：通过计算机对图像每个像元逐个地进行校正。其包括两项任务：一是像元坐标变换，二是像元灰度值重新计算（重采样）。如图 1图 1 几何校正原理（1） 坐标变换的两种方案像元坐标变换：确定原始图像和校正后的图像之间的坐标变换关系。如图 2（a） 直接法：从原始图像阵列出发，依次对其中每个像元分别计算其在输出（纠正后）图像的坐标。公式为：（b） 间接法：从空白（纠正后）图像阵列出发，依次计算空白图像每个像元在原始图像中的位置，然后把该点的灰度值计算后填充该空白像元。其公式为： 图 2 直接法与间接法其中：常用多项式来构建坐标变换的关系式， (二元)多项式可采用一次(阶)、二次、三次、n次。(x, y)：某像元的原始图像坐标；(X,Y)：纠正后同名点的地图坐标；ai, bi：多项式的系数 (i=0,1,2)。（2） 数字图像灰度值的重采样重新定位后的像元在图像中的分布是不均匀的，即输出图像像元点在输入图像中的行列号不是或不全是正数关系。因此，需要根据输出图像上的各像元在输入图像中的位置，在对原始图像按一定规则重新采样，进行亮度值的插值计算，建立新的图像矩阵。ENVI中的内插方法有：（a）最近邻法取与所计算点（x，y）周围相邻的四个点，比较它们与被计算点的距离，哪个点距离近，就取那个点的亮度值作为（x，y）点的亮度值。该方法的优点是输出图像仍然保持原来的像元值，简单、处理速度快。但这种方法会可能产生像元位置偏移，可能造成输出图像中某些地物的不连贯。（b）双线性内插法取（x，y）点周围的4个邻点，在y方向内插两次，再在x方向内插一次，得到（x，y）点的亮度值f(x，y)。双线内插法比最近邻法虽然计算量会增加，但精度会提高且对亮度不连续现象或线状地物特征的块状化现象有明显的改善。（c）三次卷积内插法取与计算点（x，y）周围16个点，先在一个方向内插，再根据计算结果在另一方向上内插，得到一个连续内插函数。这种方法对边缘有所增强，并具有均衡化和清晰化的效果，但会破坏原来的像元值，且计算量大。总的来说，最邻近法有利于保持原始图像中的灰阶，但对图像的几何结构会有破坏后两种方法虽然是对像元值的近似，但会在很大程度上保留了原图像的几何结构，如道路网、水系、地物边界。3. 结果3.1 Image to Iamge 几何校正以具有地理参考的SPOT 10米全色波段为基准影像（bldr_sp.img），对Landsat TM图像（bldr_tm.img）进行几何校正，其具体校正过程如下：3.1.1打开并显示文件选择主菜单FileOpen Image File，将SPOT影像（bldr_sp.img）和TM影像（bldr_tm.img）文件打开，分别在Display中显示。如图 3 图 3 TM影像与SPOT影像3.1.2启动几何校正模块（1）选择主菜单MapRegistrationSelect GCPs：Image to Image，打开几何校正。（2）选择显示SPOT影像的Display为基准影像（Base Image），显示TM影像的Display为待校正图像（Warp Image），点击OK，进入采集地面控制点。如图 4图 4 控制点采集3.1.3 采集地面控制点（1）在两个Display窗口中移动方框的位置，寻找明显的地物作为地面控制点（GCP）。在Zoom窗口中打开十字光标并分别定位到同名地物点上。如图 5图 5 定位同名地物（2）在Ground Control Points Selection窗口中，单击Add Point按钮，记录当前找到的地物同名点。如图 6图 6 添加控制点坐标（3）用同样的方法寻找其它的同名地物点，如果Ground Control Points Selection窗口中的Degree设为1，则选择的控制点数量达到3个时，可以使用该窗口中的Predict功能，这时在基准图像上定位一个特征点，单击Predict按钮，校正图像窗口上会自动预测相应的区域，适当调整一下位置，单击Add Point 按钮，记录同名地物点。控制点的数量越多预测点的精度越准确。（4）当控制点达到一定数量时，选择Ground Control Points Selection中的OptionsAuto Predict，打开自动预测功能，这时点击基准影像上的特征点，校正影像上会自动预测对应的点。在选择控制点时，数量要适当且要分布均匀。如图 7图 7 控制点分布（5）在Ground Control Points Selection上，点击Show List按钮，打开所有控制点的列表，如图 8图 8 控制点列表（5）在Image to Image GCP List上的OptionOrder Point by Error，按RMS值由高到低排列，对RMS值较高的点，可以直接删除或者重新定位到正确的位置，然后点击Image to Image GCP List上Updae按钮进行微调。根据实际的需要来控制RMS值的最大值。（6）在Ground Control Points Selection上，FileSave GCPs to ASCII，保存控制点。3.1.4 校正输出参数设置在ENVI中有两种校正输出方式：Warp File 和 Warp file as Image MapWarp File（1）在Ground Control Points Selection上，OptionWarp File，选择校正影像（TM数据）。（2）在校正参数对话框中，校正方法选择多项式（polynomial），多项式次数（degree）为2次，重采样选择Bilinear，背景值为0。Output Image Extent：默认是根据基准图像的大小计算，也可以进行调整，然后选择输出路径和文件名。如图 9图 9 校正参数及路径选择Warp file as Image Map（1）在Ground Control Points Selection上，OptionWarp File as Image Map，选择校正影像（TM数据）。（2）在校正参数对话框中，默认的投影参数和像元大小都与基准影像一致；这时输出像元的大小为30米；校正方法选多项式（2次）；重采样选Bilinear，背景值为0；输出影像范围默认根据基准影像的大小计算，也可以适当调整。选择输出路径和文件名，输出校正影像。如图 10图 10 投影机参数设置这两种方式的区别在于后者可以选择输出影像像元的大小。3.1.5 检验校正结果校正结果检验的基本方法是：在Display窗口中分别显示校正后的影像和基准影像，通过视窗链接（Link Display）或者地理链接（Geographic Link）进行关联显示来检验配准的效果。在视窗链接（Link Display）方法中设置一个透明度进行动态叠加显示，用中键画一个矩形框，点击左键移动该矩形框，可以实现卷帘检验的效果。从下图可以看出：校正前后线状地物基本重合。如图 11图 11 校正结果验证3.2 Image to Map 几何校正Image to Map 几何校正与Image to Iamge 几何校正基本类似，但控制点选择方式更加灵活，可从不同途径搜集或者直接从图像上读取然后键盘输入。3.2.1 打开影像选择主菜单FileOpen Image File，打开TM影像（bldr_tm.img）并显示在Display窗口中。如图 12图 12 TM影像3.2.2 启动校正模块（1）选择主菜单MapRegistrationSelect GCPs：Image to Map，打开几何校正，选择需校正图像的显示窗口，弹出校正参数设置对话框（Image to Map Registration）。（2）在Image to Map Registration对话框中，设置校正图像的投影参数和像元大小。如图 13图 13 Image to Map Registration对话框3.2.3 采集地面控制点地面控制点的可以通过以下几种方式进行采集（也可以结合使用）：（1）键盘输入在显示窗口中，移动方框选择明显它特征地物，在Zoom窗口中十字光标定位到特征地物上（可以利用键盘上的进行微调）。在Ground Control Points Selection对话框中输入该点的坐标（x，y），用同样的方法输入其它明显的特征地物点。如图 14图 14 Ground Control Points Selection对话框（2）从栅格文件中采集打开栅格文件，在Display中显示，在栅格图像上和校正的图像上寻找同名地物点并定位，在栅格定位点出右击打开快捷键，选择Pixel Locator，在Pixel Locator对话框中点击Export按钮，则该坐标会自动加载到Ground Control Points Selection对话框中对应坐标项。如图 15重复上述步骤采集其它同名地物点。图 15 从栅格数据获取控制点坐标（3）从矢量数据中采集在主菜单Open Vector File，在文件选择对话框中，文件类型选择USGS DLG（*.ddf,*.dlg），选择矢量数据，因不是ENVI的标准矢量类型，ENVI会自动转换成其标准类型，投影参数按默认类型设置。加载显示矢量数据，寻找明显的地物的特征点作为地面控制点，找到后，单击右键选择Export Map Location，会自动加载该点的坐标到Ground Control Points Selection对话框中对应坐标项。如图 16重复以上操作，选择其它的同名地物点作为控制点。图 16 从栅格数据获取控制点坐标在Ground Control Points Selection对话框中，点击按钮Show List查看以上方式选择的控制点误差的大小。如需要可以在该对话框的File菜单中保存控制点。3.2.4 校正输出参数设置在Ground Control Points Selection上，OptionWarp File，选择校正影像（TM数据）。参数设置与Image to Iamge 几何校正方法中的校正输出参数设置的“Warp file as Image Map”方式类似。如图 17图 17 校正参数设置3.2.5 校正结果比较校正好的TM影像与基准的SPOT影像分别在两个窗口中显示，通过地理联接随机选择一个特征地物，查看校正结果。如图 18图 18 Image to Map 校正结果3.3 具有已知几何信息的图像校正具有已知几何信息的图像校正方法是利用输入的几何文件生成一个地理位置查找表文件（Geographic Lookup Table，GLT），从该文件中可以了解到初始像元在输出结果中实际的地理位置。地理位置查找表文件是一个二维图像文件，文件中包含两个波段分别表示地理校正图像的行和列，文件对应的灰度值表示原始图像各像元对应的地理坐标信息，用符号整形存储，它的符号说明输出像元是对应真实的输入像元还是由邻近像元生成的填充像元。符号为正说明使用了输入像元的真实位置；符号为负说明使用了邻近像元的位置值；0值说明周围7个像元内没有邻近像元位置值。图像矫正的步骤如下：3.3.1打开文件主菜单FileOpen Image File，打开cup99hy_true.img影像和cup99hy_geo_igm文件，这两个文件被加载到波段列表中。如图 19图 19 影像显示及波段列表3.3.2 生成GLT文件(1)主菜单MapGeoreference from input GeometryBuild GLT(2)在Iput X Geometry Band对话框中，选择IGM Iput X Map作为X波段，选择IGM Iput Y Map作为Y波段(3)在Geometry Projection Information对话框中设置输出图像的投影参数。如图 20图 20Geometry Projection Information对话框（4）在Build Geometry Lookup File Parameters对话框中，选择默认输出像元大小，填写输出路径和文件名。图 21 图 21 Build Geometry Lookup File Parameters对话框3.3.3 利用GLT文件几何校正图像(1) 主菜单MapGeoreference from input GeometryGeometry from GLT(2) 在Iput Geometry Lookup File对话框中选择cup99hy_geo_glt文件和在Iput Data File对话框中选择待校正影像cup99hy_true.img。如图 22、图 23图 22 Iput Geometry Lookup File对话框图 23 Iput Data File对话框(3)在Georeference from GLT Parameters 对话框中填写背景值和输出路径。如图 24图 24 Georeference from GLT Parameters 对话框3.3.4 查看结果主菜单SpectralSPSRE ToolsGoogle Earth Bridge，将校正的结果叠加到Google Earth，可以看到校正结果基本和Google Earth基本重合。如图 25图 25 基于GLT文件的校正结果3.4 具有已知几何信息的图像校正（补充案例）待校正的数据是: 2010年11月的太阳辐射数据(W/m2)（solar_flux），在ENVI中打开并显示，如图 26图 26 太阳辐射数据（1）经纬度信息读取与显示数据每个像元的地理坐标有readme.txt文件中的说明可以计算出，如下：数据维数及地理坐标说明：The current settings for the grids aregrid: /X (degree_east) periodic (0.0) to (1.875W) by 1.875 N= 192 pts :gridgrid: /Y (degree_north) ordered (88.54196N) (86.65315N) (84.75323N) . (88.54196S) N= 94 pts :gridgrid: /T (months since 1960-01-01) ordered (Jan 1949) to (Nov 2010) by 1. N= 743 pts :gridGaussian grid with 192 longitude equally spaced and 94 latitude unequally spaced grid points由文件说明可知：经度的间隔周期是1.875，从东经0开始取192个点，经度超过180的要减去360转化为西经的度数。纬度的数据可以从已知的latitude.txt文件中获取，有94个数据。为此，可以在IDL中对经度和纬度数据分别生成一个192 samples，94 lines的数据，数据中每一位置的数值代表太阳辐射数据相应位置像元的经度或纬度信息，生成经度和纬度数据的程序如下：Pro calculate_lon_lat;根据每个像元的地理坐标建立IGM文件x=findgen(192)*1.875x97:*=x97:*-360.0temp=make_array(94, value=1.0)lon=temp#xfilename1=F:chapter10longitudeopenw,lun1,filename1,/get_lunwriteu,lun1,lonfree_lun,lun1y=ddread(F:chapter10datasolarlatitude.txt)y=reform(y)temp=make_array(192, value=1.0)lat=y#tempopenw,lun2,F:chapter10latitude,/get_lunwriteu,lun2,latfree_lun,lun2End生成的结果在ENVI中打开，填写头文件信息，注意数据类型选择Floating point，如图 27图 27 头文件信息设置打开后的纬度信息图像，如图 28图 28 纬度信息图像显示同样的方法打开，经度信息图像，如图 29图 29 经度信息图像显示（2）在EN