遥感实验7

实验7利用ETM+数据反演地表温度1实验目的与要求实习目的：利用成都市Landsat7/8遥感影像中的热红外波段数据，计算两个时期内地表亮度温度，并分析地表温度的空间特征。实习要求：⑴掌握利用遥感影像进行地表温度反演的原理与方法。⑵几何校正。由于热红外波段影像丢失空间参考信息，故需要进行几何校正（配准）处理，实习中建议使用多项式法。⑶地表亮温空间特征分析。对计算结果进行分析，重点分析地表温度的空间分异特征，并分析其与地表覆被变化之间的关系。⑷采用ArcGIS、ENVI或Erdas软件进行地表温度。⑸详细记录操作过程，比较各种算法的优缺点并整理成实验报告上交。2实验数据⑴成都市Landsat7/8遥感影像中的热红外波段数据。⑵矢量数据：成都市行政边界。3实验分组与时间⑴本次实验采取分组方式进行，每组4-5人，设组长一名，组长负责协调并记录每人的操作内容（工作量）。⑵每组交一份实验报告。⑶时间为两次实验课，如果可能完成不了，每组利用课后时间完成。4实验性质本实验属于软件操作练习。5参考材料5.1计算公式⑴辐射定标运算。即由像元DN值计算辐射亮度L：LDNgainbias⑵亮度温度计算。即由辐射亮度L计算亮度温度：)1(252,TkchTehcE)12ln(52LhckhcT)1ln(12LKKTB其中，K1的取值为666.09W/(m2srμm)，K2的取值为1282.71K。5.2ENVI下利用ETM+数据反演地表温度地表温度作为地球环境分析的重要指标，而遥感技术作为现代重要的对地观测手段，使得基于遥感图像的地表温度反演的研究越来越多。主要的地表温度反演方法有：大气校正法，单窗算法，单通道法等等。本文介绍用辐射传输方程法对地表温度进行反演。技术流程：例子数据为2002年9月2日的襄樊市LandsatETM+数据。根据数据的特点以及地表温度反演研究的技术要求，采用的技术路线为：先对LandsatETM+数据进行预处理：数据读取、辐射定标、大气校正、襄樊区域裁剪，利用大气校正，即：辐射传输方程法对其影像热红外波段数据进行操作反演，实现襄樊市地区的地表真实温度的反演研究。具体的处理流程如下：具体的实现步骤如下：第一步：准备数据热红外数据使用的是Landsat的第六波段，已经做了传感器定标、几何校正、工程区裁剪，详细流程参考上面的流程图。文件为TM6-rad-subset-jz-xiangfan.img。由TM影像（已经过大气校正）生成的NDVI数据，已经利用主菜单-BasicTools-ResizeData（SFatial/SFectral）重采样为60米分辨率，与TMi6数据保持一致，文件名为：TM-NDVI-60m.img。第二步：地表比辐射率计算物体的比辐射率是物体向外辐射电磁波的能力表征。它不仅依赖于地表物体的组成，而且与物体的表面状态(表面粗糙度等)及物理性质(介电常数、含水量等)有关，并随着所测定的波长和观测角度等因素有关。在大尺度上对比辐射率精确测量的难度很大，目前只是基于某些假设获得比辐射率的相对值，本文主要根据可见光和近红外光谱信息来估计比辐射率。(一)植被覆盖度计算计算植被覆盖度Fv采用的是混合像元分解法，将整景影像的地类大致分为水体、植被和建筑，具体的计算公式如下：FV=(NDVI-NDVIS)／(NDVIV-NDVIS)（2）其中，NDVI为归一化差异植被指数，取NDVIV=0.70和NDVIS=0.00，且有，当某个像元的NDVI大于0.70时，FV取值为1；当NDVI小于0.00，FV取值为0。利用ENVI主菜单-BasicTools-BandMath，在公式输入栏中输入：（b1gt0.7）*1+（b1lt0.）*0+（b1ge0andb1le0.7）*（(b1-0.0)/(0.7-0.0)）b1：选择NDVI图像得到植被盖度图像。(二)地表比辐射率计算根据前人的研究，将遥感影像分为水体、城镇和自然表面3种类型。本专题采取以下方法计算研究区地表比辐射率：水体像元的比辐射率赋值为0.995，自然表面和城镇像元的比辐射率估算则分别根据下式（3）（4）进行计算：εsurface=0.9625+0.0614FV-0.0461FV2（3）εbuilding=0.9589+0.086FV-0.0671FV2（4）式中，εsurface和εbuilding分别代表自然表面像元和城镇像元的比辐射率。利用ENVI主菜单-BasicTools-BandMath，在公式输入栏中输入：（b1le0）*0.995+（b1gt0andb1lt0.7）*(0.9589+0.086*b2-0.0671*b2^2)+(b1ge0.7)*(0.9625+0.0614*b2-0.0461*b2^2)b1：NDVI值；b2：植被覆盖度值。得到地表比辐射率数据。第三步：计算相同温度下黑体的辐射亮度值辐射传输方程法，又称大气校正法，其基本思路为：首先利用与卫星过空时间同步的大气数据来估计大气对地表热辐射的影响。然后把这部分大气影响从卫星高度上传感器所观测到的热辐射总量中减去。从而得到地表热辐射强度．再把这一热辐射强度转化为相应的地表温度.卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值的表达式可写为（辐射传输方程）：Lλ=[ε·B(TS)+(1-ε)L↓]·τ+L↑（4）这里，ε为地表辐射率，TS为地表真实温度，B(TS)为普朗克定律推到得到的黑体在TS的热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：B(TS)=[Lλ-L↑-τ·(1-ε)L↓]/τ·ε（5）在NASA官网()中输入成影时间以及中心经纬度，则会提供上式中所需要的参数。本专题输入的数据是襄樊市地区2002年9月2日北京时间10:30成像的Landsat7ETM+影像，影像中心的经纬度为：32.51N，111.81E。得到下图参数图：大气在热红外波段的透过率τ为0.6，大气向上辐射亮度L↑为3.39W/(m2·sr·μm)，大气向下辐射亮辐射亮度L↓为5.12W/(m2·sr·μm)。图4.412002年9月2日LandsatETM+数据的大气辅助参数利用ENVI主菜单-BasicTools-BandMath，在公式输入栏中输入：（b2-3.39-0.6*（1-b1）*5.12）/（0.6*b1）b1：60m分辨率的地表比辐射率值；b2：表示热红外波段的辐射定标值。得到了温度为T的黑体在热红外波段的辐射亮度值。第四步：反演地表温度在获取温度为TS的黑体在热红外波段的辐射亮度后，根据普朗克公式的反函数，求得地表真实温度TS：TS=K2/ln(K1/B(TS)+1)对于ETM+，K1=666.09W/(m2·sr·μm)，K2=1282.71K。利用ENVI主菜单-BasicTools-BandMath，在公式输入栏中输入：（1282.71）/alog（666.09/b1+1）-273b1：温度为T的黑体在热红外波段的辐射亮度值。得到真实的地表温度值，单位是摄氏度。第五步：结果浏览与输出在DisFlay中显示温度值，是一个灰度的单波段图像。(1)选择Tools-ColorMaFFing-DensitySlice，单击ClearRange按钮清除默认区间。(2)选择OFions-AddNewRanges，增加以下四个区间：39℃以上，红色35℃至39℃，黄色30℃至35℃，绿色低于30℃，蓝色(3)单击Apply。(4)选择File-OutFutRangetoClassImage，可以将反演结果输出。图4.42地表温度反演结果地表温度反演结果可以用于很多地方，如城市热岛监测、计算土壤湿度指数（NDVI/T）等。5.3基于MODIS数据的海表温度反演劈窗算法最初是为反演海面温度开发的，具体地说是针对NOAA/AVHRR的4和5通道设计的，后来也被用来反演地表温度，这种算法较成熟，精度也高。劈窗算法以地表热辐射传导方程为基础，利用10~13μm大气窗口内，两个相邻热红外通道(一般为10.5～11.5μm、11.5～12.5μm)对大气吸收作用的不同，通过两个通道测量值的各种组合来剔除大气的影响，进行大气和地表比辐射率的修正。表达式为：TS=T4+A(T4-T5)+B其中：TS为地表真实温度，T4和T5分别为AVHRR的4和5通道，A和B为常量。AVHRR的通道4(10.15~11.13μm)和通道5(11.15~12.15μm)恰与MODIS的第31波段(10.178~11.128μm)和32波段(11.177~12.127μm)的中心波长相对应,可将MODIS的31、32波段资料,用于劈窗算法进行地表温度计算。很多学者对这个算法进行了推演，得到很多新的算法，如覃志豪、毛克彪等人。本文就是使用其他学者推演的算法。利用MODIS数据劈窗算法反演海表温度技术流程如下图：图：技术流程图注：（1）按照本流程反演出来的结果是SST。陆地上的值可以视为无效值，若要得到正确的陆表温度，需要加入海陆分离的步骤，以及城镇和自然表面的比辐射率计算。（2）MODIS数据下载：nasa官网：下面详细介绍处理流程。第一步：读取原始数据打开File-OpenAs-Genericformats-HDF4，选择hdf文件，选择发射率数据集，点击OK，选择数据存储方式BSQ，点击OK。打开的是热红外原始数据集，第20-36波段，共16个波段，分别是：20、21、22、23、24、25、27、28、29、30、31、32、33、34、35、36波段。第二步，辐射亮度定标打开/RasterManagement/Data-SpecificUtilities/ViewHDFDatasetAttributes，选择原始的hdf数据，点击OK，选择相应的热红外数据集：EarthView1KMEmissiveBandsScaledIntegers，点击OK，得到一个文件说明面板：radiance_scales，和radiance_offset这两项参数代表波段的增益和偏移量，是辐射定标的系数。比如要计算31波段的辐射亮度，读取到scales为0.00084002，offsets为1577.33972168，带入MODIS辐射定标的公式：Radiance=scales*（DN-offsets），即可以得到该波段的辐射亮度。打开/BandAlgebra/BandMath工具，输入公式：0.00084002*(B31-1577.33972168)，（B31是第31波段的DN值），点击OK，选择第31波段数据为B31，设置路径和文件名，点击OK。得到的结果就是31波段的辐射亮度。同样的方法得到32波段的辐射亮度，公式为：0.00072970*(B32-1658.22131348)第三步：几何校正MODIS1b数据是hdf格式，自带有经纬度坐标信息，可自动进行几何校正，先对反射率数据集进行自动几何校正。反射率数据几何校正：（1）打开反射率数据集：File-OpenAs-EOS-MODIS，选择hdf数据文件打开；（2）工具/GeometricCorrection/ReprojectGLTwithBowtieCorrection，选择反射率数据集文件，点击SpatialSubset，选择大亚湾的大概位置，点击SpectralSubet，选择2和19波段。点击OK；（3）打开工具/GeometricCorrection/GeoreferencebySensor/GeoreferenceMODIS，选择反射率数据集，点击SpatialSubset，选择大亚湾的