单窗算法反演地表温度教程1.1 算法原理1.1.1 单窗算法单窗算法( MW 算法)是覃志豪于2001 年提出的针对TM 数据只有一个热红外波段的地面温度反演算法。

经过众多学者验证，单窗算法具有很高的反演精度，且同样适用于ETM+和landsat 8数据。

公式如下：T s (a6(1 C6 D6) (b6(1 C6 D6) C6 D6 )T sensor D6T a)/C6 式中，LST为地表温度(K)，T sensor是传感器上的亮度温度( K)，T a 是大气平均温度( K)；a、b为参考系数，当地表温度为0-70℃时，a = -67.355351，b =0.458606；C 、D 为中间变量，计算公式为：式中，为地表比辐射率，为地面到传感器的大气总透射率。

因此单窗算法反演地表温度的关键是计算得到亮度温度T senso、地表比辐射率、大气透射率和大气平均作用温度T a。

1.1.2 参数计算1.1.2.1 辐射亮温计算利用Planck 公式将图像像元对应传感器辐射强度值转换为对应的亮度温度值。

公式如下式中，T senso 为亮度温度值； L 影像预处理后得到的光谱辐射值，2单位为w/(m2 sr m)，K1 、K2 为常量，可由数据头文件获取。

计算图像辐射亮温之前，需采用辐射定标参数将像元灰度值DN 转换为热辐射强度值，公式如下：式中，M L 为增益参数，A L 为偏移参数，该参数可直接在影像通文件数据中获取，且ENVI 软件中已经集成，不需要自己在查找。

1.1.2.2 地表比辐射率计算根据覃志豪针对TM 影像提出的混合像元分解法来确定区域地表福辐射率。

对于城市区域，我们简单的将其分为水体、自然表面和建筑表面三因此针对混合像元尺度上的地表比辐射率通过下式来种，估算：式中，为混合像元的地表比辐射率；P V为植被覆盖率；R V为植被的温度比率；R M为建筑表面的温度比率；V 表示植被法地表比辐射率，m 表示建筑表面的地表比辐射率；d 表示辐射校正项。

根据覃志豪经验公式，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 -(P1v)当Pv=0.5 时, d ε =0.0019RV 和Rm 的确定采用覃志豪等提出的经验公式计算:RV=0.9332+0.0585*PVRm=0.9886+0.1287*PVPV 为植被覆盖率可用下式计算:PV=(NDVI-NDVI S)/(NDVI V-NDVI S)其中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)大气透射率估计方程相关系数平方（R2）标准误差（SEE）高气温0.4-1.6 τ6=0.97429-0.08007w0.99611 0.0023681.6-3.0 τ6=1.03141-0.99827 0.00253低气温0.4-1.6 τ6=0.98200-70.09611w0.99463 0.003341.6-3.0 τ6=1.05371-0.99899 0.00237接下来只需要计算出大气水分含量便可计算得到大气透射率。

学者杨景梅根据地面、高空气象要输资料，提出利用地面湿度计算整层大气可降水量，由此计算得出整层大气水汽含量的方法（杨景梅邱金恒，用地面湿度参量计算我国整层大气可降水量及有效水汽含量的方法的研究，大气科学，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 ，双击此工具，打开辐射定标的数据，进行相关的参数设置进行大气校正。

注意：如果在多光谱数据辐射定标时Scalefactor 值忘记设置，可在本步骤中打开辐射定标数时设置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.tx”t 中找到，具体名称：DATE_ACQUIRED = 2013-05-12 ；SCENE_CENTER_TIME = 02:55:26.6336980Z；8)大气模型参数选择：Sub-Arctic Summer(根据成像时间和纬度信息选择)；9)气溶胶模型Aerosol Model ：Urban，气溶胶反演方法AerosolRetrieval：2-band（K-T ）；10) 其他参数按照默认设置即可。

11) 多光谱参数设置中，K-T 反演选择默认模式：Defaults->Over-Land Retrieval standard(600:2100)波谱响应函数：默认指向\ProgramFiles\Exelis\ENVI51\classic\filt_func\landsat8_oli.sli把它重新指向：..\Progr amFiles\Exelis\ENVI51\resource\filterfuncs\landsat8_oli. sli注：这是因为ENVI5.1 版本的一个小bug，即Classic中的L8 的波谱响应函数不正确，另外一个一劳永逸的方法是：将“ filterfuncs 中”的Landsat8_oli.sli 和Landsat8_oli.hdr 两个文件拷贝覆盖：“... filt_func 中” 的两个文件。

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

见下图经过 FLAASH 校正的影像基本去除了空气中水汽颗粒等因子的影响，植被的波谱曲线趋于正常12） 高级参数设置：根据内存大小设置 Tile Size （Mb ）：100（8g 物理内存），其他参数默认即可，详细见下图(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 gt0.5)*0.0019 b1：植被覆盖度Pv第二步，计算R V 与R m，在波段运算器中分别输入下列式中：R V : 0.9332+0.0585*b1R m : 0.9886+0.1287*b1b1: 植被覆盖度Pv第三步，计算地表比辐射率利用B10 波段计算地表比辐射率: b1*b2*0.98672+(1-b1)*b3*0.96767+b4 利用B11 波段计算地表比辐射率: b1*b2*0.9899+(1-b1)*b3*0.97515b1: 植被覆盖度Pv; b2:植被表面温度比率RV; b3:建筑表面温度比率R m; b4：指的是dεK1 K2 需要根据不同影像查表得到：(1321.0789)/alog(774.8853/b1 +1)b1：波段10 或者波段11 经过辐射校正后的数据 L1.2.3 地表温度计算第一步，计算C: b1*b2b1：地表比辐射率；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。