基于单窗算法反演地表温度的ENVI操作教程

单窗算法反演地表温度教程

1.1 算法原理 1.1.1 单窗算法

单窗算法（MW算法）是覃志豪于2001年提出的针对TM数据只有一个热红外波段的地面温度反演算法。经过众多学者验证，单窗算法具有很高的反演精度，且同样适用于ETM+和landsat 8数据。公式如下：

Ts?(a6(1?C6?D6)?(b6(1?C6?D6)?C6?D6)Tsensor?D6Ta)/C6

式中，LST为地表温度（K），Tsensor是传感器上的亮度温度（K），Ta是大气平均温度（K）；a、b为参考系数，当地表温度为0-70℃时，a = -67.355351，b = 0.458606；C 、D为中间变量，计算公式为：

式中，为地表比辐射率，为地面到传感器的大气总透射率。 因此单窗算法反演地表温度的关键是计算得到亮度温度Tsenso、地表比辐射率、大气透射率和大气平均作用温度Ta。 1.1.2 参数计算 1.1.2.1 辐射亮温计算

利用Planck公式将图像像元对应传感器辐射强度值转换为对应的亮度温度值。公式如下

式中，Tsenso为亮度温度值；L?影像预处理后得到的光谱辐射值，

2w/(m?sr??m)，K1 、K2为常量，可由数据头文件获取。 单位为

计算图像辐射亮温之前，需采用辐射定标参数将像元灰度值DN转换为热辐射强度值，公式如下：

式中，ML为增益参数，AL为偏移参数，该参数可直接在影像通文件数据中获取，且ENVI软件中已经集成，不需要自己在查找。 1.1.2.2 地表比辐射率计算

根据覃志豪针对TM影像提出的混合像元分解法来确定区域地表福辐射率。对于城市区域，我们简单的将其分为水体、自然表面和建筑表面三种，因此针对混合像元尺度上的地表比辐射率通过下式来估算：

式中，为混合像元的地表比辐射率；PV为植被覆盖率；RV为植被的温度比率；RM为建筑表面的温度比率；射率，

m表示建筑表面的地表比辐射率；d

V表示植被法地表比辐

表示辐射校正项。

根据覃志豪经验公式，V=0.986；m=0.972。

根据Aster提供的常用地物比辐射率光谱库Landsat 8数据特点，对于B10波段，V=0.98672；m=0.96767。对于B11波段，V=0.98990；

m=0.977515。

在地表相对较平整的区域,一般可取dε=0;在地表起伏较大的区域,dε可根据植被的构成比例估计。研究表明热辐射相互作用在植被与裸土分别占一半时达到最大,所以可以根据如下经验公式来估计dε：

当Pv≤0.5时,dε=0.003 8Pv 当Pv>0.5时,dε=0.003 8 (1-Pv) 当Pv=0.5时, dε=0.0019

RV和Rm的确定采用覃志豪等提出的经验公式计算:

RV=0.9332+0.0585*PV Rm=0.9886+0.1287*PV PV为植被覆盖率可用下式计算: PV=(NDVI-NDVIS)/(NDVIV-NDVIS)

其中NDVI是归一化植被指数,NDVIV、NDVIS分别是植被和裸土的NDVI值,由于没有详细的区域植被和土壤光谱资料,采用覃志豪经验值以NDVIV=0.70和NDVIS=0.05来进行植被覆盖度的近似估计。

1.1.2.3 大气平均作用温度的计算

对于大气平均温度，覃志豪等也提供了一个经验公式，该公式是根据Modtran所提供的标准大气推导出来的，如下表所示，其中T0为近地表温度(K)。

一般研究区在重庆的，选择中纬度夏季平均大气。 1.1.2.4 大气透射率的计算

由于单窗算法对输入的参数较敏感，因此需要较精确的大气透射率参数来反演地表温度。影响大气透射率的因素较多,如气压、气温、气溶胶含量、大气水分含量、O3、CO2、CO、NH4等对热辐射传导均有不同程度的作用,从而使地表的热辐射在大气中的传导产生衰减。因此,准确的大气透射率求算比较复杂,需要较详细的大气剖面数据，但在实际应用中很难获取详细的大气剖面数据。研究表明大气水分含量对大气透射率的变化影响较大,而其它因素由于很难在短时间内发生较大变化而对大气透射率的影响较小。因此,水分含量就成为大气

透射率估计的主要考虑因素。

覃志豪等根据大气中的水汽含量动态变化对大气透射率的影响程度进行了研究，并对大气透射率进行了估计，如表3.1。其中:w是大气水汽含量。以重庆为研究区的研究一般选用是气温较高时的估算公式。

大气透射率估计方程

水分含量大气剖面 大气透射率估计方程

w/(g·cm-2) 高气温 低气温

0.4-1.6 1.6-3.0 0.4-1.6 1.6-3.0

τ6=0.974290-0.08007w τ6=1.031412-0.11536w τ6=0.982007-0.09611w τ6=1.053710-0.14142w

相关系数平方(R2) 0.99611 0.99827 0.99463 0.99899

标准误差(SEE) 0.002368 0.002539 0.003340 0.002375

接下来只需要计算出大气水分含量便可计算得到大气透射率。学者杨景梅根据地面、高空气象要输资料，提出利用地面湿度计算整层大气可降水量，由此计算得出整层大气水汽含量的方法（杨景梅 邱金恒，用地面湿度参量计算我国整层大气可降水量及有效水汽含量的方法的研究，大气科学，2002.01）公式如下：

式中，ω为大气水汽含量（g/cm2），W为整层大气可降水量（cm）；C0 C1为经验系数，以重庆为研究对象，一般可选取成都的经验值来计算。C0= 0.1274；C1=0.6878.计算大气可降水量公式如下：

式中，W为整层大气可降水量（cm），e为地面水汽压（hPa），a0 a1为经验系数，可从以下公式计算得出。

式中，φ为地理纬度（°），H为海拔高度。重庆市主城区纬度取29.5，平均高程取0.35 km，则有。

ω=0.19604e+0.04691

其中：e是绝对水汽压（hPa）。

采用WMO在1996年公布的《气象仪器和观测方法指南》中公布的饱和水汽压E公式。对于纯水面用公式计算，适用于温度-45~60℃。 E=6.112exp[17.62*T0/(243.12+T0)] RH=e/E*100 推算出：

e=6.112exp[17.62*T0/(243.12+T0)]*RH

式中，RH为相对湿度，可以通过气象观测数据获得，T0为气温(℃)。 1.2 地表温度反演步骤

1.2.1 影像预处理 （1）辐射定标

选择Radiometric Correction/Radiometric Calibration。在File Selection对话框中，选择数据LC81230322013132LGN02_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（10.9），打开Radiometric Calibration面板。

Scale factor 不能改变，否则后续计算会报错。保持默认1即可。

（2）大气校正

本教程选择Flaash 校正法。FLAASH Atmospheric Correction，双击此工具，打开辐射定标的数据，进行相关的参数设置进行大气校正。

注意：如果在多光谱数据辐射定标时Scale factor值忘记设置，可在本步骤中打开辐射定标数时设置single scale faceor 值为0.1，若已设置，则默认值为1即可。

1) Input Radiance Image：打开辐射定标结果数据； 2) 设置输出反射率的路径，由于定标时候；

3) 设置输出FLAASH校正文件的路径，最优状态：路径所在磁盘空间足够大；

4) 中心点经纬度Scene Center Location：自动获取；

5) 选择传感器类型：Landsat-8 OLI；其对应的传感器高度以及影像数据的分辨率自动读取；

6) 设置研究区域的地面高程数据；

7) 影像生成时的飞行过境时间：在layer manager中的Lc8数据图层右键选择View Metadata，浏览time字段获取成像时间；

注：也可以从元文件“LC81230322013132LGN02_MTL.txt”中找到，具体名称：DATE_ACQUIRED = 2013-05-12；SCENE_CENTER_TIME = 02:55:26.6336980Z；

8) 大气模型参数选择：Sub-Arctic Summer（根据成像时间和纬度信息选择）；

9) 气溶胶模型Aerosol Model：Urban，气溶胶反演方法Aerosol Retrieval：2-band（K-T）；

10) 其他参数按照默认设置即可。

11) 多光谱参数设置中，

K-T反演选择默认模式：Defaults->Over-Land Retrieval standard（600:2100）

波谱响应函数：默认指向..

\\Program Files\\Exelis\\ENVI51\\classic\\filt_func\\landsat8_oli.sli 把

它

重

新

指

向

：

..\\Program

Files\\Exelis\\ENVI51\\resource\\filterfuncs\\landsat8_oli.sli

注：这是因为ENVI5.1版本的一个小bug，即Classic中的L8的波谱响应函数不正确，另外一个一劳永逸的方法是：将“\\Exelis\\ENVI51\\resource\\filterfuncs”中的

Landsat8_oli.sli

和

Landsat8_oli.hdr两个文件拷贝覆盖：“...\\ENVI51\\classic\\filt_func” 中的两个文件。

否则SWIR1波段大气校正后的结果全为0。见下图

12) 高级参数设置：根据内存大小设置Tile Size（Mb）：100(8g物理内存)，其他参数默认即可，详细见下图

经过FLAASH校正的影像基本去除了空气中水汽颗粒等因子的影响，植被的波谱曲线趋于正常。

1.2.1.1 地表比辐射率计算 （1）植被覆盖度计算

第一步，利用多光谱遥感影像，计算研究区域归一化植被指数NDVI，ENVI软件中有工具集成，使用即可。

第二步，利用波段运算工具，根据覃志豪经验公式，计算研究区植被覆盖度。

(b1 gt 0.7)*1+(b1 lt 0.05)*0+(b1 ge 0.05 and b1 le 0.7)*((b1-0.05)/(0.7-0.05))

b1：NDVI （2）地表比辐射率计算 第一步，计算dε

(b1 le 0.5)*0.0038*b1+(b1 eq 0.5)*0.0038*(1-b1)+(b1 gt 0.5)*0.0019 b1：植被覆盖度Pv

第二步，计算RV与Rm，在波段运算器中分别输入下列式中：

RV : 0.9332+0.0585*b1 Rm : 0.9886+0.1287*b1 b1: 植被覆盖度Pv

第三步，计算地表比辐射率

利用B10波段计算地表比辐射率: b1*b2*0.98672+(1-b1)*b3*0.96767+b4 利用B11波段计算地表比辐射率: b1*b2*0.9899+(1-b1)*b3*0.97515

b1: 植被覆盖度Pv; b2:植被表面温度比率RV; b3:建筑表面温度比率Rm; b4：指的是dε

1.2.2 辐射亮温计算

K1 K2需要根据不同影像查表得到：

（1321.0789）/alog（774.8853/b1 +1）

b1：波段10 或者波段11经过辐射校正后的数据L?

1.2.3 地表温度计算

第一步，计算C: b1*b2

b1：地表比辐射率； b2：大气透射率 第二步，计算D: (1-b2)*(1+(1-b1)*b2) b1：地表比辐射率； b2：大气透射率 第三步，计算LST:

(-67.355351*(1-b1-b2)+b3*(0.458606*(1-b1-b2)+b1+b2)-b2*295.4022465)/b1

b1:指的是C; b2：指的是D; b3：指的是辐射亮温；b4：指的是大气平均作用温度Ta

本文来源：https://www.bwwdw.com/article/smk8.html