遥感图像预处理

1.自定义坐标系1.1.坐标系概括常用坐标系分为地理坐标系和投影坐标系。1.1.1.地理坐标系地理坐标系中含椭球体(spheroid)和大地基准面(datum)两个重要部分。1)椭球体。地球并不是一个规则球体，因此我国常用以下椭球体模拟地球，见表1。表1我国常用椭球体椭球体名称年份长半轴/m短半轴/m扁率WGS841984年6378137.06356752.31:298.257克拉索夫斯基(Krasovsky)1940年6378245.06356863.01:298.3IAG0751975年6378140.06356755.31:298.2572)大地基准面大地基准面指参考椭球与WGS84参考椭球的相对位置关系（3个平移，3个旋转，1个缩放），可用3个、4个或者7个参数进行描述。每个椭球体都对应一个或多个大地基准面。1.1.2.投影坐标系投影坐标系是利用一定数学法则将地球表面上的经纬线网表示到平面上，属于平面坐标系。一个投影坐标系包含椭球体、投影方法和大地基准面。我国普遍采用的投影方法是高斯-克吕格投影，在欧美国家称为横轴莫卡托投影(TransverseMercator)。为保证地图精度，采用分带投影法，有3°分带和6°分带两种。对于每分带的投影，需指定中央经线以及坐标轴原点X、Y坐标。我国在北半球，为避免Y坐标出现负值，规定统一将各带Y值加500km。为避免各带坐标重复，可在不同带间Y值冠以带号，即通用坐标。商业软件坐标系分为标准坐标系和自定义坐标系，我国实际应用中往往需要自定义坐标系，以北京54坐标系和西安80坐标系为例讲解如何在Envi中自定义坐标系。1.2.Envi自定义坐标系1.2.1.定义椭球体使用记事本打开HOME\ITT\IDL80\products\envi48\map_proj\ellipse.txt，将“Krasovsky,6378245.0,6356863.3”(Beijing-54坐标系使用)和“IAG-75,6378140.0,6356755.3”(Xian-80坐标系使用)加入文本末端。1.2.2.定义基准面使用记事本打开HOME\ITT\IDL80\products\envi48\map_proj\datum.txt，将“D_BEIJING_1954,Krasovsky,-12,-113,-41”(Beijing-54坐标系使用)和“D_Xian_1980,IAG-75,0,0,0”(Xian-80坐标系使用)加入文本末端。注：此处的D_BEIJING_1954和D_Xian_1980在Envi4.7版本以上不可更改，否则会转换失败。1.2.3.定义投影系在Envi任何用到投影坐标的功能模块中，都可以新建坐标系（在任何地图投影选择对话框中，点击New按钮），也可通过主菜单→Map→CustomizeMapProjection，弹出CustomizedMapProjectionDefinition对话框，如图1。Beijing19546Degree111E图1定义投影系依序定义如下参数：1)ProjectionName：投影坐标系名称；2)ProjectionType：选择投影类型，如TransverseMercator；3)ProjectionDatum：选择基准面；4)Falseeasting：东偏距离500000，若需X轴含带号，添加相应带号于500000前；5)Falsenorthing：北偏距离，填写0；6)Latitude：中央纬线，填写0；7)Longitude：中央经线，填写某带的中央经线，如117；8)Scalefactor：中央经线长度比，填写0.9996，精度与高斯-克吕格等同。选择Projection→AddNewProjection，将投影添加到Envi使用的列表中。选择File→SaveProjection，存储新添或更改过的投影信息。打开HOME\ITT\IDL71\products\envi48\map_proj\map_proj.txt，即可看到新建的坐标信息，如图2。至此，投影系已添加完毕，可在Envi中使用。图2添加完成的投影坐标系1.3.图像投影转换描述一个栅格文件地理位置信息由两部分组成：坐标信息（map）和投影信息（projection）。坐标信息由起始点像素坐标以及对应的地理（投影）坐标和像素大小组成；投影信息就是坐标系信息。一般来说，如果坐标信息丢失，这个文件将会失去坐标；投影信息是用来描述坐标信息，如果投影信息丢失，可以重新设定。下面演示如何将北京54坐标系的栅格图像转化为西安80坐标系。1.3.1.定义初始投影系1)打开.\data\0-坐标定义\f49e011021.img文件；2)在AvailableBandsList中，MapInfo节点右击，单击EditMapInformation，弹出对话框；3)单击ChangeProjection按钮，在弹出对话框中选择目标投影系，如图3；图3设定影像投影系4)点击OK，AvailableBandsList中的MapInfo下拉节点坐标信息发生了相应改变，表明投影信息已成功赋予图像，如图4。图4设定影像投影系前后对比1.3.2.投影转换选择主菜单→Map→ConvertMapProjection，在弹出的对话框中设置如下参数，如图5：1)ChangeProjection：选择投影目标系名称；2)ConversionParameters：Method选用Polynomial，PolynomialDegree为2；3)OutputResultto：选择Memory。图5投影转换参数设置单击OK后，执行投影转换，转换结果如图6。图6投影转换后坐标信息2.图像几何校正几何校正是利用地面控制点和几何校正数学模型来矫正非系统因素产生的误差，同时也是将图像投影到平面上，使其符合地图投影系统的过程；由于校正过程中会将坐标系统赋予图像数据，所以此过程包括了地理编码(geo-coding)。2.1.基于自带定位信息的几何校正对于重返周期短、空间分辨率较低的卫星数据，如AVHRR、MODIS、SeaWiFS等，地面控制点的选择有相当难度。因此，可以用卫星传感器自带的地理定位文件进行几何校正。在Envi中，基本方法为主菜单→Map→Georeference传感器名称。本部分以MODISLevel1B级数据介绍基于自带定位信息的几何校正方法。2.1.1.打开数据文件打开.\data\1-Modis\MOD02HKM.A2002248.0345.005.2007348121959.hdf文件。2.1.2.选择校正模型选择主菜单→Map→GeoreferenceMODIS，弹出InputMODISFile对话框（图7），单击文件名，选择校正的文件单击OK，进入下一步。图7选择要处理的MODIS波段2.1.3.设置输出参数a)在GeoreferenceMODISParameters对话框中（图8），设置输出坐标系；b)在NumberWarpPoints中，键入X、Y方向校正点的数量。在X方向的校正点数量应该小于等于51，在Y方向的校正点数应该小于等于行数；c)EnterOutputGCPFilename：设置控制点输出文件；d)PerformBowTieCorrection：用来消除MODIS的“蝴蝶效应”，默认Yes。如图8，点击OK，进入下一步。图8选择投影坐标系2.1.4.设置校正参数a)系统自动计算起始点坐标值、像元大小、图像行列数据，也可根据要求更改；b)Background设置为0；c)设置输出影像路径。图9校正参数及输出路径如图9，点击OK，开始执行校正，经过错误!未找到引用源。两个步骤后，最终校正结果如图11所示。图10影像校正图11校正前(左)后(右)影像2.2.基于GLT的FY3数据校正GLT几何校正法利用输入的集合文件生成一个GLT(geographiclookuptable)地理位置查找表文件。它是一个二维图像文件，包含地理校正图像的行列两个波段，文件对应的灰度值表示原始图像每个像素对应的地理位置坐标信息，用符号整型存储，符号为正时，说明使用了真实的像元位置值，符号为负时，说明使用了邻近像元的位置值，值为0说明周围7个像元内没有邻近像元位置值。GLT文件包含初始图像每个像元的地理定位信息，它的校正精度是很高的，避免了通过地面控制点利用二次多项式几何校正法对低分辨率图像数据的处理。下面以我国(FY3)气象卫星可见光红外扫描辐射计(VisibleandInfraRedRadiometer,VIRR)数据为例，进行GLT几何校正。2.2.1.安装补丁安装HDF5补丁，将.\data\2-FY3中的open_hdf5_event.sav文件拷贝到HOME\ITT\IDL80\products\envi48\save_add文件夹中，重启ENVI。2.2.2.打开文件a)主菜单→File→OpenExternalFile→GenericFormats→HDF5，选择文件Z_SATE_C_BAWX_20090104070730_P_FY3A_VIRRX_GBAL_L1_20090104_0510_1000M_MS.HDF（文件类型选择*.*）b)如图12，选择EV_RefSB图像数据，点击ImporttoENVI。图12HDF5文件查看c)同上步，将定位经纬度文件打开(Latitude和Longitude)。2.2.3.生成GLT文件a)选择主菜单→Map→GeoreferencefromImputGeometry→BuildGLT。b)在InputXGeometryBand对话框中，选择精度Longitude作为X波段，如图13。图13选择X波段c)由于X波段左边边缘为0值，因此有必要对边缘进行掩膜处理，点击SpatialSubset按钮，在SelectSpatialSubset对话框中，设置Samples为3（去掉开始3个像素），如图14。图14掩膜去除0值像素d)在InputXGeometryBand对话框单击OK后，进入InputYGeometryBand对话框，选择纬度Latitude为Y波段。e)点击OK后，进入GeometryProjectionInformation对话框，设置输入输出投影参数。f)点击OK后，进入BuildGeometryLookupFile对话框，填写输出像元大小（默认），旋转角度(Rotation)为0（正上方为北），选择保存路径及文件名，如图15。图15输出GLT文件设置g)点击OK后，开始生成GLT文件，如图16。图16生成GLT文件h)生成的GLT文件含SampleLookup和LineLookup两个波段，如图17。图17生成的GLT文件2.2.4.利用GLT文件几何校正图像a)选择主菜单→Map→GeoreferencefromInputGeometry→GeoreferencefromGLT。b)在文件选择对话框中，依次选择GLT校正文件、待校正影像。c)在输出对话框中，选择输出路径及文件名，单击OK进行校正，校正后影像已包含MapInfo信息，如图18。图18校正后影像2.3.ImagetoImage几何校正本部分介绍以具有地理参考的SPOT410m全色波段为基础，对Landsat5TM30m影像进行校正的过程，流程如图19。图19几何校正一般流程显示图像文件采集地面控制点计算误差选择几何模型重采样输出检验校正结果结束误差太大2.3.1.打开影像文件选择主菜单→File→OpenImageFile，将.\data\3-几何校正中的bldr_sp.img和bldr_tm.img文件打开，并分别显示在Display中。2.3.2.启动几何校正模块a)选择主菜单→Map→Registration→SelectGCPs:ImagetoImage，打开几何校正模块。b)选择显示SPOT文件的Display为基准图像(BaseImage)，显示TM文件的Displ