Landsat温度反演

定量遥感实验2 基于Landsat TM数据和大气校正法的地表温度反演一实验目的学习并掌握基于Landsat 5TM遥感影像，运用辐射传输方程法（大气校正法）对地表温度进行反演的方法。

二研究方法实验数据：2005年9月8日长春市Landsat 5TM 图像（景号：118-29和118-30）。

实验流程：完整流程涉及Landsat TM 的数据读取、辐射定标、反演模型、遥感反演过程、反演结果验证等内容。

软件功能：主要采用ENVI 主模块中的Landsat TM 数据读取、辐射定标、波段运算、结果统计分析等。

该方法需要进行大气校正，消除大气层对地表辐射能量的影响，这就需要从卫星观测得到的热辐射能量中扣除大气层的辐射分量，并利用热红外波段（Band 6）范围内的地表发射率作为参数，反演出地表的真实温度。

技术路线为：（1）Landsat TM数据预处理：数据读取、辐射定标。

（2）相关辅助数据的确定与查找：大气上行辐射以及下行辐射，采用数据当天的大气透过率信息等；（3）采用大气校正法（辐射传输方程法）利用Landsat TM band 6 进行地表温度反演：首先获取地表比辐射率值，其次，计算黑体在热红外波段的辐射亮度；最后，利用普朗克公式的反函数反演出整个研究区域的地表温度分布情况。

图1 流程图三实验步骤1、Landsat TM 数据读取和定标（1）数据读取在ENVI 主菜单中，选择File→ Open External File→ Landsat→Geo TIFF with Metadata ，打开。

其中包含两种数据：可见光波段数据（HRF）和热红外波段数据（HTM ）图2 Landsat TM 数据读取列表（2）数据辐射定标主菜单→Basic Tools→Preprocessing→Calibration Utilities→Landsat Calibration ，选择需要定标的数据，包括可见光波段数据和热红外波段数据，该过程功能可以从读取文件中直接获取实验数据的元数据（成像时间，定标参数等等）。

图3 辐射定标2、裁剪工程区（1）镶嵌影像主菜单→Map→Mosaicking→Georeferenced，在Map Based Mosaic窗口中Import File，输入两副已经辐射定标的热红外波段遥感影像。

影像导入完成后，可修改影像的投影信息，以便提高影像精度。

镶嵌模块中，主菜单上选择Options→Change Base Projection。

打开转变投影功能后，根据影像信息设置投影。

选择影像，单击右键选择Edit Entry，可以通过此功能，改变影像镶嵌参数。

完成点击应用File→Apply即可进入镶嵌过程，选择重采样方法等。

图4遥感影像镶嵌（2）裁剪工程区①打开要裁剪的区域范围的矢量数据：主菜单→File→Open Vector→选择文件类型Shapefile(*.shp)，注意在弹出的Import Vector Files Parameters 对话框选择正确的投影类型（与遥感影像的投影坐标一致）。

点击OK之后ENVI自动将矢量文件转为EVF格式。

图5 矢量数据读取②将矢量数据转为ROI：在Availabel Vectors List选择数据，通过File→Export Layers to ROI，在Select Data File to Associate with new ROIS 中选择镶嵌好的热红外影像，再在Export EVF Layers to ROI中选择Convert all records of an EVF layer to one ROI，点击OK，从而实现矢量数据转为ROI。

图6 Export EVF Layers to ROI对话框③裁剪栅格数据：在ENVI主菜单Basic Tools →Subset Data via ROIs，在Select Input File to Subset via ROI 中选择要裁剪的热红外波段影像（镶嵌好的），OK。

出现Spatial Subset via ROI Parameters对话框，在Select Input ROIs 中选择建立的ROI，选项Mask pixels outside of ROI ?选择Yes，在Mask Background Value默认为0，这样得到长春市的热红外波段遥感影像。

图7 Spatial Subset via ROI Parameters对话框对于可见光波段数据的裁剪方法与上述方法相同，由于可见光波段数据和热红外波段数据的分辨率不同，在利用该矢量数据裁剪可见光波段数据的时候需要重复②步骤，完成裁剪操作即可。

3、可见光和近红外波段大气校正（快速校正）主菜单→Basic Tools→Preprocessing→Calibration Utilities→QUick Atmospheric Correction，选择长春市裁剪的可见光波段，在QUick Atmospheric Correction Parameters对话框中，Sensor Type选择Landsat TM。

图8 QUick Atmospheric Correction Parameters对话框4、地表比辐射率值计算物体的比辐射率是物体向外辐射电磁波的能力表征。

它不仅依赖于地表物体的组成，而且与物体的表面状态(表面粗糙度等)及物理性质(介电常数、含水量等)有关，并随着所测定的波长和观测角度等因素有关。

在大尺度上对比辐射率精确测量的难度很大，目前只是基于某些假设获得比辐射率的相对值，在此主要根据可见光和近红外光谱信息来估计比辐射率。

（1）归一化植被指数NDVI利用TM3 、4波段的象元DN值利用公式（1）求得归一化植被指数NDVI：NDVI = (NIR-R)/(NIR+R) （1）NIR 和R 分别是TM 的近红外波段(波段4)和红光波段(波段3)的DN 值。

直接采用ENVI 主模块下的NDVI 模块获取：主菜单→Transform→NDVI：图9 NDVI值（2）植被覆盖度计算植被覆盖度F v 采用的是混合像元分解法，将整景影像的地类大致分为水体、植被和建筑，具体的计算公式如下：F V = [(NDVI-NDVI S)/ (NDVI V -NDVI S)]（2）其中，NDVI 为归一化差异植被指数，取NDVI V = 0.70 和NDVI S = 0.00 ，且有，当某个像元的NDVI 大于0.70 时，F V取值为1；当NDVI 小于，F V取值为0。

利用ENVI 主菜单→Basic Tools→Band Math ，在公式输入栏中输入：（）*1+（b1 lt 0.0）*0+（）*（(b1-0.0)/(0.7-0.0)）。

b1: 表示获取的NDVI 值。

图10 植被覆盖度Fv图11 植被覆盖度影像数据统计结果（3）地表比辐射率根据前人的研究，将遥感影像分为水体、城镇和自然表面3 种类型。

本实验采取以下方法计算研究区地表比辐射率：水体像元的比辐射率赋值为0.995 ，自然表面和城镇像元的比辐射率估算则分别根据下式（3）（4）进行计算：VV 2 （3）VV 2 （4）式中，εsurface 和εbuilding 分别代表自然表面像元和城镇像元的比辐射率。

利用ENVI 主菜单→Basic Tools→Band Math ，在公式输入栏中输入：（b1 le 0）*0.995+（）*(0.9589 + 0.086*b2 -0.0671*b2^2)+(b1 ge 0.7)*(0.9625 + 0.0614*b2 -0.0461*b2^2) 。

b1：NDVI 值b2：植被覆盖度值。

图12 30m 分辨率的地表比辐射率值上述公式计算的结果分辨率为30m ，由于热红外波段数据分辨率为60m，故而需要将其分辨率重采样为60m ，为反演地表温度提供基础数据。

具体的操作步骤：利用ENVI 主菜单→Basic Tools→Resize Data（Spatial/Spectral），出现对话框设置：图13 影像数据的重采样图14 60m 分辨率的地表比辐射率值5、大气参数获取卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，以及地面的真实辐射亮度经过大气层之后到达卫星传感器的能量。

地面的真实辐射亮度为同温度黑体的辐射亮度值L T与地物发射率ε的乘积ε·L T。

大气校正法的表达式可写为：Lλ= [ε·L T + (1-ε)L↓]·τ + L↑ （5）这里，T为地表真实温度，τ为大气在热红外波段的透过率.则温度为T的黑体在热红外波段的辐射亮度LT为：L T = [Lλ -L↑-τ·(1-ε)L↓]/τε （6）在NASA 官网(/) 中输入成像时间以及中心经纬度，则会提供上式中所需要的参数。

本实验输入的数据是2005 年9月8日北京时间2：9成像的Landsat5TM 影像，影像中心的经纬度为：N, E 。

得到下图参数图: 大气在热红外波段的透过率τ为0.67 ，大气向上辐射亮度L↑为W/(m2·sr·μm) ，大气向下辐射亮辐射亮度L↓为W/(m2·sr·μm) 。

图15 2005 年9 月8 日Landsat TM 数据的大气辅助参数利用ENVI 主菜单→Basic Tools→Band Math ，在公式输入栏中输入：（b2-7*（1-b1）*）/（*b1）b1：60m 分辨率的地表比辐射率值；b2：表示热红外波段辐射定标值。

获取的黑体在热红外波段的辐射亮度结果图如下：图16 热红外波段辐射亮度结果图17 热红外波段辐射亮度统计结果6、地表温度反演在获取热红外波段辐射亮度值以后根据普朗克公式的反函数，求得地表真实温度T：T = K2/ln(K1/L T + 1)对于TM，K1=6W/(m2·sr·μm)，K2 =1260.56K。

利用ENVI 主菜单→Basic Tools→Band Math ，在公式输入栏中输入：（12）/alog（6/b1 +1）b1 ：热红外辐射亮度值。

计算结果如下：图18 区域LST 统计结果7、反演结果分析图19地表的真实温度T 的分布图图20 长春市遥感影像假彩色合成利用ENVI 软件，采用单波段彩色变换方法对地表真实温度T 的灰度图进行密度分割，得到地表真实温度的分布情况（图19）。

其中：(1)红色代表的温度范围是：39℃以上，为异常点；(2)黄色代表的温度范围是：30℃至39℃；(3)绿色代表的温度范围是：25℃至30℃；(4)蓝色代表的温度范围是：低于25℃。