遥感实验7教材.docx

1、遥感实验7教材实验7 利用ETM+数据反演地表温度1 实验目的与要求实习目的：利用成都市Landsat 7/8遥感影像中的热红外波段数据，计算两个时期内地表亮度温度，并分析地表温度的空间特征。实习要求： 掌握利用遥感影像进行地表温度反演的原理与方法。 几何校正。由于热红外波段影像丢失空间参考信息，故需要进行几何校正（配准）处理，实习中建议使用多项式法。 地表亮温空间特征分析。对计算结果进行分析，重点分析地表温度的空间分异特征，并分析其与地表覆被变化之间的关系。 采用ArcGIS、ENVI或Erdas软件进行地表温度。 详细记录操作过程，比较各种算法的优缺点并整理成实验报告上交。2 实验数据 成

2、都市Landsat 7/8遥感影像中的热红外波段数据。 矢量数据：成都市行政边界。3 实验分组与时间 本次实验采取分组方式进行，每组4-5人，设组长一名，组长负责协调并记录每人的操作内容（工作量）。 每组交一份实验报告。 时间为两次实验课，如果可能完成不了，每组利用课后时间完成。4 实验性质本实验属于软件操作练习。5 参考材料5.1 计算公式 辐射定标运算。即由像元DN值计算辐射亮度L： 亮度温度计算。即由辐射亮度L计算亮度温度：其中，K1的取值为666.09W/(m2srm)，K2的取值为1282.71K。5.2 ENVI下利用ETM+数据反演地表温度地表温度作为地球环境分析的重要指标，而遥

3、感技术作为现代重要的对地观测手段，使得基于遥感图像的地表温度反演的研究越来越多。主要的地表温度反演方法有：大气校正法，单窗算法，单通道法等等。本文介绍用辐射传输方程法对地表温度进行反演。技术流程：例子数据为2002年9月2日的襄樊市Landsat ETM+数据。根据数据的特点以及地表温度反演研究的技术要求，采用的技术路线为：先对Landsat ETM+数据进行预处理：数据读取、辐射定标、大气校正、襄樊区域裁剪，利用大气校正，即：辐射传输方程法对其影像热红外波段数据进行操作反演，实现襄樊市地区的地表真实温度的反演研究。具体的处理流程如下：具体的实现步骤如下：第一步：准备数据热红外数据使用的是La

4、ndsat的第六波段，已经做了传感器定标、几何校正、工程区裁剪，详细流程参考上面的流程图。文件为TM6-rad-subset-jz-xiangfan.img。由TM影像（已经过大气校正）生成的NDVI数据，已经利用主菜单-Basic Tools-Resize Data（SFatial/SFectral）重采样为60米分辨率，与TMi6数据保持一致，文件名为：TM-NDVI-60m.img。第二步：地表比辐射率计算物体的比辐射率是物体向外辐射电磁波的能力表征。它不仅依赖于地表物体的组成，而且与物体的表面状态(表面粗糙度等)及物理性质(介电常数、含水量等)有关，并随着所测定的波长和观测角度等因素有

5、关。在大尺度上对比辐射率精确测量的难度很大，目前只是基于某些假设获得比辐射率的相对值，本文主要根据可见光和近红外光谱信息来估计比辐射率。(一) 植被覆盖度计算计算植被覆盖度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主菜单-Basic Tools-Band Math，在公式输入栏中输入：（b

6、1 gt 0.7）*1+（b1 lt 0.）*0+（b1 ge 0 and b1 le 0.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和bui

7、lding分别代表自然表面像元和城镇像元的比辐射率。利用ENVI主菜单-Basic Tools-Band Math，在公式输入栏中输入：（b1 le 0）*0.995+（b1 gt 0 and b1 lt 0.7）*(0.9589 + 0.086*b2 - 0.0671*b22)+(b1 ge 0.7)*(0.9625 + 0.0614*b2 - 0.0461*b22)b1：NDVI值；b2：植被覆盖度值。得到地表比辐射率数据。第三步：计算相同温度下黑体的辐射亮度值辐射传输方程法，又称大气校正法，其基本思路为：首先利用与卫星过空时间同步的大气数据来估计大气对地表热辐射的影响。然后把这部分大气影

8、响从卫星高度上传感器所观测到的热辐射总量中减去。从而得到地表热辐射强度再把这一热辐射强度转化为相应的地表温度.卫星传感器接收到的热红外辐射亮度值L由三部分组成：大气向上辐射亮度L，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值的表达式可写为（辐射传输方程）：L = B(TS) + (1-) L + L （4）这里，为地表辐射率，TS为地表真实温度，B(TS)为普朗克定律推到得到的黑体在TS的热辐射亮度，为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：B(TS) = L- L-(1-) L

9、/ （5）在NASA官网(http:/atmcorr.gsfc.nasa.gov/)中输入成影时间以及中心经纬度，则会提供上式中所需要的参数。本专题输入的数据是襄樊市地区2002年9月2日北京时间10:30成像的Landsat7 ETM+影像，影像中心的经纬度为：32.51N，111.81 E。得到下图参数图：大气在热红外波段的透过率为0.6，大气向上辐射亮度L为3.39 W/(m2srm)，大气向下辐射亮辐射亮度L为5.12W/(m2srm)。图4. 41 2002年9月2日Landsat ETM+数据的大气辅助参数利用ENVI主菜单-Basic Tools-Band Math，在公式输入栏

10、中输入：（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/(m2srm)，K2 =1282.71K。利用ENVI主菜单-Basic Tools-Band Math，在公式输入栏中输入：（1282.71）/alog（666.09/b1 +1）-273b1：温度为T的黑体

11、在热红外波段的辐射亮度值。得到真实的地表温度值，单位是摄氏度。第五步：结果浏览与输出在DisFlay中显示温度值，是一个灰度的单波段图像。(1) 选择Tools-Color MaFFing-Density Slice，单击Clear Range按钮清除默认区间。(2) 选择OFions-Add New Ranges，增加以下四个区间： 39以上，红色 35至39，黄色 30至35，绿色 低于30，蓝色(3) 单击Apply。(4) 选择File-OutFut Range to Class Image，可以将反演结果输出。图4.42 地表温度反演结果地表温度反演结果可以用于很多地方，如城市热岛监

12、测、计算土壤湿度指数（NDVI/T）等。5.3 基于MODIS数据的海表温度反演劈窗算法最初是为反演海面温度开发的，具体地说是针对NOAA/AVHRR的4和5通道设计的，后来也被用来反演地表温度，这种算法较成熟，精度也高。劈窗算法以地表热辐射传导方程为基础，利用1013m 大气窗口内，两个相邻热红外通道(一般为10.511.5m、11.512.5m)对大气吸收作用的不同，通过两个通道测量值的各种组合来剔除大气的影响，进行大气和地表比辐射率的修正。表达式为：TS= T4+ A (T4- T5) + B其中：TS为地表真实温度，T4和 T5分别为AVHRR的4和5通道，A和B为常量。AVHRR 的

13、通道4(10.1511.13m) 和通道5 (11.1512.15m) 恰与MODIS 的第31 波段 (10.178 11.128m) 和32 波段 (11.177 12.127m ) 的中心波长相对应, 可将MODIS 的31、32 波段资料, 用于劈窗算法进行地表温度计算。很多学者对这个算法进行了推演，得到很多新的算法，如覃志豪、毛克彪等人。本文就是使用其他学者推演的算法。利用MODIS数据劈窗算法反演海表温度技术流程如下图：图：技术流程图注：（1）按照本流程反演出来的结果是SST。陆地上的值可以视为无效值，若要得到正确的陆表温度，需要加入海陆分离的步骤，以及城镇和自然表面的比辐射率计算

14、。（2）MODIS数据下载：nasa官网：http:/modis.gsfc.nasa.gov/.下面详细介绍处理流程。第一步：读取原始数据打开File-Open As-Generic formats-HDF4，选择hdf文件，选择发射率数据集，点击OK，选择数据存储方式BSQ，点击OK。打开的是热红外原始数据集，第20-36波段，共16个波段，分别是：20、21、22、23、24、25、27、28、29、30、31、32、33、34、35、36波段。第二步，辐射亮度定标打开/Raster Management/Data-Specific Utilities/View HDF Dataset A

15、ttributes，选择原始的hdf数据，点击OK，选择相应的热红外数据集：Earth View 1KM Emissive Bands Scaled Integers，点击OK，得到一个文件说明面板：radiance_scales，和radiance_offset这两项参数代表波段的增益和偏移量，是辐射定标的系数。比如要计算31波段的辐射亮度，读取到scales为0.00084002，offsets为1577.33972168，带入MODIS辐射定标的公式：Radiance=scales*（DN-offsets），即可以得到该波段的辐射亮度。打开/Band Algebra/Band Math工

16、具，输入公式：0.00084002*(B31-1577.33972168)，（B31是第31波段的DN值），点击OK，选择第31波段数据为B31，设置路径和文件名，点击OK。得到的结果就是31波段的辐射亮度。同样的方法得到32波段的辐射亮度，公式为：0.00072970*(B32-1658.22131348)第三步：几何校正MODIS1b数据是hdf格式，自带有经纬度坐标信息，可自动进行几何校正，先对反射率数据集进行自动几何校正。反射率数据几何校正：（1）打开反射率数据集：File-Open As-EOS-MODIS，选择hdf数据文件打开；（2）工具/Geometric Correction

17、/Reproject GLT with Bowtie Correction，选择反射率数据集文件，点击Spatial Subset，选择大亚湾的大概位置，点击Spectral Subet，选择2和19波段。点击OK；（3）打开工具/Geometric Correction/Georeference by Sensor/Georeference MODIS，选择反射率数据集，点击Spatial Subset，选择大亚湾的大概位置，点击Spectral Subet，选择2和19波段。点击OK。图：选择数据同时选择空间子集和光谱子集（4）在Georeference MODIS Parameter面板

18、设置为UTM， WGS84，49带，并输出GCP文件。图：Georeference MODIS Parameter面板（5）点击OK，在结果输出面板，设置路径和文件名输出。31、32波段的几何校正：(1) 双击/Geometric Correction/Registration/Warp from GCPs: Image to Map Registration，选择上一步保存出来的GCP文件，设置坐标系为：UTM WGS84 49带，设置输出分辨率均为1000米；(2) 点击OK，选择辐射定标后的31波段数据，点击Spatial Subset，起止行列号设置为和反射率子集一致，点击OK，(3)

19、 输出面板上，设置输出的方法为：Triangulation，重采样方法为：Bilinear，设置文件名输出；同样的方法对32波段辐射定标的结果进行几何校正。第四步：海表温度反演本文使用的是高懋芳，覃志豪等，MODIS数据反演地表温度的基本参数估计方法中分裂窗算法模型进行海表温度反演，旨在学习ENVI中的操作流程。算法为：Ts= A0+ A1*T31- A2* T32 (1)其中：Ts是 地 表 温 度 ( K ) , T31和T32分 别 是M OD I S 第 31 和 32 波段的亮度温度；A0、A1和A2是分裂窗算法的参数，分别定义如下：A0= D32(1-C31-D31) / (D32

20、C31-D31C32) a31- D31(1-C32-D32) /(D32C31- D31C32) a32 (2)A1= 1 + D31/ ( D32C31- D31C32) + D32( 1 -C31- D31) / ( D32C31- D31C32) b31 (3)A2= D31/ ( D32C31- D31C32) + D31( 1 - C32- D32) / ( D32C31- D31C32) b32 (4)式中，a31，b31，a32和 b32是常量, 根据 MODIS的波段特征确定, 在地表温度 0 5 e 范围内, 这些常量 分 别 可 取 a31=-64.60363 , b31

21、= 0.440817, a32=-68.72575, b32= 0.473453。上述公式的中间参数分别计算如下：Ci=ii() (6)Di= 1 -i() 1 + ( 1 -i) i() (7)式中：i 是指 MODIS的第31和32波段, 分别为 i= 31 或 32；i()是视角为的大气透过率；i是波段 i 的地表比辐射率。由以上公式可以看出, 该算法要求卫星遥感器的31和32波段数据来计算星上亮度温度, 同时还要求已知大气透过率和地表比辐射率, 才能进行地表温度的反演。下面是详细操作步骤：（1）大气透过率计算大气透过率i ()是计算地表温度的基本参数, 通常是通过大气水汽含量来估计。经

22、过前人研究，可以用MODIS第2和19波段来反演大气水分含量，然后再根据大气水分含量与大气透过率之间的关系来估计大气透过率。对于MODIS图像中的任何一个像元，其可能的大气水分含量用下式估计 （8）式中：是大气水分含量(g*cm-2)，和常量，分别= 0.02 和=0.6321；19和2分别是MODIS第19和2波段的地面反射率。使用bandmath工具计算大气水分含量：表达式：(0.02-alog(b19/b2)/0.6321)2B19：第19波段反射率B2：第2波段反射率大气透过率的计算中，水汽是最主要的考虑因素，毛克彪等将MODTRAN等大气模型模拟出来的两者的关系，应用到MODIS数据

23、中，提高了地表温度反演的精度和实时性，本文采用模拟效果较好的指数关系模拟方程，拟合度达到了0.99以上，公式为：31= 2.89798-1.88366*exp-/（-21.22704） （9）32= -3.59289+4.60414*exp-/（-32.70639） （10）式中，是水汽含量。使用bandmath工具计算大气透过率：31波段大气透过率表达式：2.89798-1.88366*exp(b1/21.22704)B1：大气水分含量。32波段大气透过率表达式：-3.59289+4.60414*exp(b1/32.70639)B1：大气水分含量。（2）地表比辐射率的估算地表比辐射率主要取决

24、于地表的物质结构，对MODIS来说，大致分为水面、城镇和自然表面。对于反演来说，利用混合像元分解的方法，根据植被覆盖率来计算自然表面和城镇的比辐射率，水体的可以用常量：31水体=0.996，32水体=0.992。有了这些参数，我们就可以计算C31、C32、D31、D32中间参数，BandMath表示式分别为：C31=0.996*b31C32=0.992*b32D31=(1-b31)*(1+(1-0.996)*b31)D32=(1-b32)*(1+(1-0.996)*b32)其中 B31：31波段大气透过率B32：32波段大气透过率（3）亮度温度的计算将图像DN值定标维热辐射强度之后, 可用Pl

25、anck函数求解出星上亮度温度, 计算公式如下：Ti= Ki 2/ l n ( 1 + Ki 1/Ii)式中，Ki 1和 Ki 2是常量，对于第 i = 31 波段, 分别为 K31 , 1= 729 .541636 Wm-2sr-1um-1，K31 , 2= 1304.413871K ; 对于第 i = 32 波段, 为 K32 , 1= 474 . 6 84780 Wm-2sr-1um-1，, K31 , 2= 1196 . 978785 K。使用bandmath工具计算31和32的亮温。31波段亮温：1304.413871/alog(1+729.541636/b31)B31: 31波段辐

26、射亮度值32波段亮温：1196.978785 /alog(1+474.684780/b32)B32：32波段辐射亮度值（4）A 0, A 1 和 A 2参数计算接下来我们计算A 0, A 1 和 A 2参数，Bandmath表达式分别为：A 0=b4*(1-b1-b3)/(b4*b1-b3*b2)*(-64.60363)-b3*(1-b2-b4)/(b4*b1-b3*b2)*(-68.72575)A 1=1+b3/( b4*b1-b3*b2)+b4*(1-b1-b3)/( b4*b1-b3*b2)* 0.440817A 2=b3/(b4*b1-b3*b2)+b3*(1-b2-b4)/(b4*b

27、1-b3*b2)* 0.473453其中，b1：C31B2：C32B3：D31B4：D32（5）温度计算把这些参数带入公式1中计算温度值，Bandmath表达式为：Ts=b0+b1*b31-b2*b32-273其中：b0：A0参数B1：A1参数B2：A2参数B31：B31亮温值B32：B32亮温值如下为得到的最终地表温度反演结果，单位为摄氏度。通过/Statistics/Compute Statistics工具统计看到最初的结果会有一些数量很少的异常值，几十个像素，这些异常值大多是处于影像的边缘。如下图为统计的结果，小于-40的像元占0.6%，大于32度像元占0.1%。可以将这些像元处理，如用以下bandmath处理：-40b132。图：初始结果统计图：最终反演结果