遥感地质实习报告

1、中国地质大学研究生课程论文封面遥感地质实习作业要求这是广州市北部地区的ETM遥感数据，包括：、ETM1,2,3,4,5,7多光谱（30M；ETM8Pan（15M；ETM61,62热红外（60M）。、1：50万地质图；、1：50万地理图：地名，水系等。、是MAPGI賂式需转换为ArcGIS格式应用。要求：1、利用水系对ETM数据做坐标的几何校正；2、运用校正后的数据求NDVI并叠加在RGB321图像上；3、运用校正后的数据将多光谱图像与可见光图像作数据融合；4、利用融合后的多光谱数据做主成分分析，生成PC12、3的彩色图像，并对结果进行解释；5、结合自己所学专业制作一幅专题图像：专题图像，地名标

2、注，图框，坐标，比例尺，图例等要素；并作解译。作业打印成册，附彩图。2010年6月15日前交到地勘楼532室。遥感地质课程一、利用水系对ETM数据做坐标的几何校正；数据：投影转换后的1:50万地理图；多光谱遥感影像数据：30空间分辨率的ETM1-5715m空间分辨率的ETM860m空间分辨率的ETM61ETM62软件平台：MapGISArcGIS和ENVI处理过程：Stepl:预处理：地理图投影变换；数据：1：50万地理图；1：50万地质图；软件平台：MapGIS6.7;处理过程：将地理图、地质图所有图层设定相同参数一起投影转换：首先，在MapGIS主菜单中，选择“实用服务”-“投影变换”，在

3、打开的投影变换系统中打开所有图层，可知图层信息以“毫米”为单位的屏幕坐标显示；其次,在投影变换系统的菜单栏选择“投影转换”-“工作区直接投影转换”：坐标值发生了很大变化，较原图层位臵发生了较大距离的位移，其坐标值已经换算成以meters为单位的实际坐标值。Step2:将投影转换后的地理图、地质图数据转成能在ArcGIS中显示的格式shp；Step3:在ArcMap中打开shp格式的地理图和地质图；由于此时无任何地理坐标和投影坐标信息，只有图层对应的实际坐标值信息（由投影转换后得到的数据），并且点图层和面图层显示的信息在ARCGIS显示不集中，原先MapGIS格式的图元如居民地的名字、经纬度数字

4、及不同地层的颜色显示在ArcGIS中属性信息部分缺失，但此时我们用水系信息纠正遥感影像如下：Step4:将遥感影像在ENVI中转换成TIF格式，在ArcMap中打开并显示；在ENVI中，打开需要格式转换的影像，选择File-SaveFileAs-TIFF/GeoTIFF,选择转换影像，输出即可；此时，地理图没有地理坐标系和投影坐标系，但影像有投影坐标系，在ArcMap中打开地理图和影像时对应位臵如下图：发现：虽然地理图和影像都有以meter为单位的坐标值，可以同时显示，但位臵上仍有一些差异。影像的地理经纬度实际覆盖地理图的一部分，如果能将地理图和影像通过经纬度大致重合在一起，在此基础上做几何校

5、正预处理。Step5:通过地理图框4个角的经纬度值纠正地理图，使之能与影像大致重合；在ArcMap中打开经投影转换后的地理图层，设臵Layers的地理坐标为GCS_Xian_1980寸所有地理图层纠正，纠正后地理图所显示的经纬度值与实际相符；Step6:纠正地理图的地理坐标系统，打开影像，此时地理图和影像基本吻合；在ArcCatalog中设臵纠正后的地理图的地理坐标系统为：GCS_Xian_1980打开ArcCatalog，右击所有图层的“属性”Properties，在XYCoordinateSystem标签中选择地理坐标系统GeographicCoordinateSystems为GCS_Xi

6、an_1980可选择其他地理坐标系统，但不能选择投影系统，因为地理坐标系统以经纬度显示，投影系统以meter为单位的值显示，地理图只有正确的经纬度值，影像既有精确的经纬度值又有准确的以meter为单位的值），对于影像，其本身带有一定的坐标信息，若不修改其地理坐标，在ArcGIS中打开时可根据图层Layer的地理信息自动转换成以经纬度显示，此时地理图与影像基本吻合了；Step7:利用栅格图像几何纠正工具选控制点纠正地理图和影像图，使之尽量完全吻合；在用纠正后的带有地理坐标信息的地理图来纠正影像之前，必须使影像和图层Layer带有投影坐标信息，这里设为WGS_1984_UTM_Zone_49地理图

7、可以通过ArcToolBox-DataManagerTools-ProjectionsandTransformation-Feature来转换地理坐标系统到投影坐标系统，控制点尽量选择均匀，在控制点选完后，保存好控制点，并且存储纠正后的影像时格式选为TIF，用Georeferencing-Rectify可保存结果影像。彩图：上：几何纠正后影像和地理图叠加图左：部分几何纠正前叠加图右：部分集合纠正后叠加图二、运用校正后的数据求NDVI并叠加在RGB321图像上；数据：校正后的多光谱数据ETM1-57；方法：NDVI,密度分割法叠加；软件平台：ENVI4.5处理过程：Stepl:利用校正后的多光谱

8、数据ETM1-57求NDVI灰度图；NDVI=（NIR-RED）/（NIR+RED）这里ETM3为RED波段，ETM4为NIR波段，因植被在红波段有一反射波谷,在近红外波段有一反射峰，经归一化处理后灰度图中植被与背景地物信息的反差变大，使得NDVI灰度图中灰白色地物为覆盖度不一的植被，暗色为水体，建筑物，道路等背景信息；Step2:将NDVI灰度图按覆盖程度分级，采用密度分割法；选择Image窗口的Overlay-DensitySlice,参数设臵如下：根据影像的实际情况，观察NVDI灰度直方图及反复试验，最后将植被覆盖度分为3级：-0.7063-0.1无植被覆盖区-0.10.23植被覆盖度低

9、（白色）（青色）0.230.4植被覆盖度一般（草黄色）0.40.5541植被覆盖度咼（绿色）Step3:保存密度分割后的分类影像图；将真彩色图RGB321叠加在经过密度分割的图上，需将密度分割图导出存储成分类影像图，选择DensitySlice窗口，File-OutputRangestoClassImage,设臵存储路径即可；Step4:将NDVI密度分割图叠加在RGB321真彩色影像上；在Image窗口，选择Overlay-Classification,在弹出的对话框中给要叠加在影像上的3类植被区打勾，如下：NDVI密度分割与RGB32嗔彩色叠加部分放大叠加图部分真彩色图三、运用校正后的数据

10、将多光谱图像与可见光图像（ETM8作数据融合；数据：校正后的多光谱数据ETM1-5,7;校正后的高空间分辨率数据ETM8方法：主成分变换融合方法；处理过程：Stepl:将6个波段的多光谱数据ETM1-5,7重采样，使得其分辨率和影像大小与ETM8相同；选择BasicTools-ResizeData,选中所要进行重采样的多光谱数据，在ResizeDataParameters对话框设臵参数：按照多光谱影像大小与ETM8影响大rjFhl.口冈書*畫Layer(G:ndjm.5t_4fi7-tif):adjujt_-|X|图1.多光谱数据的ETM4重采样前后对比图（右：采样后）Step2:将重采样后的

11、多光谱数据进行主成分正变换；选择Transform-PrincipalComponents-ForwardPCRotation-ComputerNewStatisticsandRotate,选择重采样后的多光谱数据，输入主成分正变换结果图存放路径即可。小的比例设臵xfac和yfac，则多光谱数据影像大小就被调整成与ETM8相同。在重采样插值方法中Resampling，有三种方法选择：最近邻法NearestNeighbor、双线性内插法Bilinear、三次卷积内插法CubicConvolution;这里选择第二种方法，因为最近邻法虽然简单，但处理后图像的亮度具有不连续性，影响了精确度，三次卷积

12、内插法虽然会使图像细节更清晰，但计算量大且没有对多光谱数据和ETM8选控制点校正位臵，三次卷积效果也不是很明显，因此选择双线性内插法，计算量一般，精度提高且效果不错，缺点是对图像起到平滑作用，使得对比度明显的分界线变得模糊，但并不影响后续工作。图2.主成份正变换图PC1PC2PC合成Step3:将高光谱数据ETM8与第一主成分PC1做直方图匹配；1）分别将TEM8和PC1直方图打开（Image窗口：Enhance-lnteractiveStretching）2）分别在ETM8和PC1的直方图窗口中，选择Histogram_Source-band;3）在ETM8直方图窗口中，将Strecth_t

13、ype选为Arbitrary，以便于用指定的直方图曲线来拉伸；图3.TEM8拉伸类型设为“Arbitrary”后直方图4）用鼠标将PC1影像直方图的输入（InputHistogram标签）拖动至在ETM8S方图的输出窗口（OutputHistogram）中，然后点击Apply；图4将PC1输入直方图与ETM8输出直方图匹配后ETM8影像图5）在PC1直方图窗口中，选择Options-HistogramParameters，记录下HistogramMin和HistogramMax两个值，分别为-103.198和403.563；6）在ETM8S方图窗口中，选择File-ExportStretch,

14、将刚才记录的两个值分别填入OutputMin和OutputMax中，再将OutputDataType改为“FloatingPoint”，存储即可。Step4:将经过直方图匹配后的ETM8替换PC1;选择BasicTools-LayerStacking,在对话框中选入经直方图匹配后的ETM8和PC2-6,且将ETM8放在第一波段，使之成为PCtStep5:将替换后的合成波段数据做主成分反变换选择Transform-PrincipalComponents-lnversePCRotation，在对话框中选择需要做反变换的数据即可。显示数据时对Image的线性拉伸即可。e#?(R;K1Z(J.Hjer

15、(G;nrijnst_45T-1I曾幻(It;InvPt：Rnndb皿亦*_旺JtM)二上：主成分融合后543波段组合图左：部分多光谱543重采样图右：部分主成分融合图四、利用融合后的多光谱数据做主成分分析，生成PC12、3的彩色图像，并对结果进行解释主成分分析法的基本原理：主成分分析法是一种降维的统计方法，它借助于一个正交变换，将其分量相关的原随机向量转化成其分量不相关的新随机向量，这在代数上表现为将原随机向量的协方差阵变换成对角形阵，在几何上表现为将原坐标系变换成新的正交坐标系，使之指向样本点散布最开的p个正交方向，然后对多维变量系统进行降维处理，使之能以一个较高的精度转换成低维变量系统，

16、再通过构造适当的价值函数，进一步把低维系统转化成一维系统。Stepl:选择Transform-PrincipalComponents-ForwardPCRotation-ComputerNewStatisticsandRotate,选择融合后的多光谱数据，输入主成分正变换结果图存放路径即可。图像变换是指将图像从空间域转换到变换域的过程，进行图像变换的目的就是为了使图像的处理过程简化。由于遥感图像的不同波段之间往往存在着很高的相关性，从直观上看，就是不同波段的图像很相似，因而从提取有用信息的角度考虑，有相当大一部分数据是多余和重复的。主成分分析的目的就是把原来多波段图像中的有用信息集中到数目尽可

17、能少的新的主成分图像中，并使这些主成分图像之间互不相关，也就是说各个主成分包含的信息内容是不重叠的，从而大大减少总的数据量，并使图像信息得到增强。主成分变换变换在遥感图像处理中又称作主成份分析或主分量分析。在遥感图像分类中，常常利用主成分分析算法来消除特征向量中各特征之间的相关性，并进行特征选择。从几何意义来看，变换后的主分量空间坐标系与变换前的多光谱空间坐标系相比旋转了一个角度，而且新的坐标系的坐标轴一定指向数据信息量较大的方向。以二维空间为例，假定某图像像元的分布呈椭圆状，那么经过旋转后，新坐标系的坐标轴一定分别指向椭圆的长半轴和短半轴方向一一主分量方向，因为长半轴这一方向信息量最大五、结

18、合自己所学专业制作一幅专题图像：专题图像，地名标注，图框，坐标，比例尺，图例等要素；并作解译。数据：校正后的热红外数据ETM62方法：ETM6热红外温度反演法软件平台：ENVI4.5,ArcGIS;1、专题图像制作名称及目的利用ETM专感器的热红外波段61和62反演地表真实温度，结合该地地理地质图件制作地表温度专题图，为后续工作(如找温泉点、土壤湿度反演、城市热岛效应研究等)做基础。2、本文研究区采用的温度反演方法原理2.1方法原理流程首先将图像亮度值DN转换为辐射亮度J(TM6)，其方程式为：LtTM6)=几忒TN16)+严呈严咏(1)式中，Lmax为传感器可探测的最大辐射亮度(单位为:Wm

19、2sr-1卩m1)，即340K全辐射体发出的平均辐射亮度丄min为可探测的最小辐射亮度,即200K全辐射体发出的平均辐射亮度。对于LandsatETM7,高增益b62,Lma=12.65,Lmin=3.2;低增益b61，Lma=17.04,Lmin=0.0；由于TM6处于大气透射性(透过率为80流右)相当好的大气窗口,如果暂时忽略地表发射的电磁辐射在到达卫星的途中所受的大气影响，那么，即可通过辐射亮度(Le,Q推算地表相对温度(Tb),其方程式为：T=f:式中，TB为地表温度(K),可令K1=cJ-5,K2=C2入-1,它们为校正系数；对于LandsatETM7,K1、K2为已知常量丄为平均辐

20、射亮度(Le,产Lb)。K=666.093Wnfsr-1卩m1K=1282.7108K这种在不考虑大气等各种影响因素的情况下，并假定地物为全辐射体的前提下，通过测定地表的辐射亮度求出的地表温度称为亮度温度(或辐射温度)。再经过大气订正及地表发射率订正等即可求出地表的真实温度。地物的亮度温度(Tb)与真实温度(TR)之间存在下列关系式：T_02RAnC&xp(C2/AT-61式中，&为地物的比辐射率，不同类型地表，&的变化较大,C2为第二辐射常量(0.0144mK)，入为热红外的有效波长，TM6的中心波长为1.475卩m2.2地物的比辐射率&计算方法本文研究区参考Sobrino等人提出的NDVI

21、TherhsoldsMehtodNDVP进行地表比辐射率估计，由于研究区影像环境的特殊性，在计算地表比辐射率过程中参数由实际情况所定，具体方法如下：a、当NDVK0.0时，地表物体为水体和建筑物，通过使用监督分类提取出的水体图层，将水体和建筑物分开；水体的比辐射率为0.995，建筑物的比辐射率为0.970;b、当0.0vNDVK0.5时，地表认为是由不同覆盖度的植被构成，包含土壤，像元可认为是土壤和植被的组合，此时，地表比辐射率计算公式为：=vPv+s(1-Pv)+d其中，v为植被的比辐射率，s为裸土的比辐射率，R为像元中的植被覆盖度。Pv=(NDVI-NDVJnin)/(NDVIma尸NDV

22、Imin)这里，NDVJnax=0.5,NDVImin=0.0(原方法中NDVJnin=0.2，NDVmax=0.5)。当Pv0.5时，d=0.003796(1-Pv)。当地表比较平坦时，d可忽略不计；植被的比辐射率为0.986，裸土的比辐射率为0.972；c、当NDVA0.5时，地表认为是全植被覆盖区，地表比辐射率为0.986。3、实验流程及结果Step1:将图像亮度值DN转化为辐射亮度L62；在ENVI中，选择BandMath工具，输入计算公式：3.2+(12.65-3.2)/255)*float(b1)bl选择b62数据即可；Step2:通过辐射亮度L62推算地表相对温度Tb；同步骤a,

23、输入计算公式为：1282.7108/alog(1.0+666.093/b2)b2选择L62数据即可；Step3:计算地物比辐射率a、用纠正后的多光谱数据计算NDVI,并分别用阈值0.0和0.5给出NDVK0.0,0.0VNDVK0.5,NDVI0.5的二值图图层掩模，其中，NDVK0.0图层由水体和建筑物组成，0.0KNDVK0.5图层为不同覆盖度的植被，为裸土和植被的混合像元，NDVA0.5图层为全植被覆盖区；在本文研究区，对影像分为4类地物：水体、建筑物、全植被、植被和裸土混合像元，进行地物比辐射率的赋值计算，因此，需要将水体和建筑物分开；b、用监督分类的最大似然法单独提取出水体图层，则N

24、DVK0.0图层中的建筑物即可与水体分开，四类二值图层掩模分别如下图所示：c、对于水体、建筑物、全植被覆盖区掩模可直接赋值，而对于植被和土壤混合像元，须经过计算得到地表比辐射率，为方便计算，本文研究区便忽略d&不计，直接在BandMath里输入公式：0.986*Pv+0.972*(1.0-Pv)计算，最后将各图层信息合并，得到整个研究区的地表比辐射率；因为此时的比辐射率图空间分辨率为30m而计算温度的热红外波段空间分辨率为60m需要将30m空间分辨率裁减成60m以保证计算精度，可在ENVI的ResizeData工具栏中进行操作，最终得到空间分辨率为60m的研究区地表比辐射率图，如下：总体看来，

25、水体的亮度最高，其次为植被、裸土和植被混合区，最后为建筑物，这正与水体的比辐射率0.995、植被的0.986、裸土的0.972、建筑物的0.970相对应；Step4:根据亮度温度Tb求取地表真实温度Tr；在ENVI的BandMath工具中输入公式：0.0144心1.475*alog(b1*exp(0.0144心1.475*b2)-b1+1.0)b1选择地表比辐射率数据，b2选择亮度温度Tb数据；计算得到的真实温度最低为11.353444摄氏度，最高为36.188804摄氏度，根据具体情况将真实温度分为3个等级，1120,20.000132，32.000136,分别为低、中、高三个等级，最后将等级图导成分类图，存成tif格式即可；Step5:在ArcGIS中制作地表温度专题图；将在ENVI中计算得到的温度图存成TIF格式，直接在ArcGIS中打开：首先，在ArcCatalog中新建一个面图层，其大小与影像大小一致，用来对点文件：居民点，线文件：河流进行裁剪，裁剪工具为ArcToolbox-AnalysisTools-Extract-Clip,将裁剪后的影像跟地理图信息显示在ArcMap中；其次，打开一个新的ArcM