使用大气校正法对landsat-8tirs地表温度进行反演

软件平台

ENVI5.3软件

具体方法及原理

原理：首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相应的地表温度。

实现原理：卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值Lλ的表达式可写为（辐射传输方程）：
Lλ = [εB(TS) + (1-ε)L↓]τ + L↑
式中，ε为地表比辐射率，TS为地表真实温度(K)，B(TS)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度B(TS)为：
B(TS) = [Lλ - L↑- τ(1-ε)L↓]/τε
Ts可以用普朗克公式的函数获取。
TS = K2/ln(K1/ B(TS)+ 1)

对于TM，K1 =607.76 W/(m2µmsr)，K2 =1260.56K。
对于ETM+，K1=666.09 W/(m2µmsr)，K2 =1282.71K。
对于TIRS Band10，K1= 774.89 W/(m2µmsr)，K2 = 1321.08K。
这些走可以从MTL文件中获得
下图为例子是landsat7ETM+的数据

由以上公式我们可以得知如果要反演出该地区的温度情况，我们需要获得两个参数：其一是大气剖面参数，其二是地表比辐射率参数。其中大气剖面参数可以从NASA提供的网站上去查询(http://atmcorr.gsfc.nasa.gov/)。在此网站上输入影像成像时间和经纬度，接着选使用的卫星以及自己邮箱即可获得大气剖面参数

具体操作流程及步骤

辐射定标

此处以资阳的landsat8数据为例
首先打开数据，对第10波段进行辐射定标，获得辐射亮度图像
具体步骤：
在Toolbox工具箱中，选择Radiometric Correction/Radiometric Calibration。在File Selection对话框中，选择数据LC81290392014225LGN00_MTL_Thermal，单击Spectral Subset选择Thermal Infrared1（10.9），打开Radiometric Calibration面板。如图

计算NDVI并获得植被覆盖度数据

在ToolBox工具箱中选择spectral/vegetation/NDVI,对多光谱数据求得NDVI，然后使用Sobrino的NDVI阈值法计算地表比辐射率，具体公式：ε=0.004Pv+0.986
Pv为植被覆盖度，通常使用公式Pv = [(NDVI- NDVISoil)／(NDVIVeg - NDVISoil)]来计算，其中，NDVI为归一化植被指数，NDVISoil为完全是裸土或无植被覆盖区域的NDVI值，NDVIVeg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVIVeg = 0.70和NDVISoil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。
在Toolbox里面选择Band Algebra/bandmath输入公式：
(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，计算之后可以得到植被覆盖率图像
接着在bandmath中输入表达式：0.004b1 + 0.986

b1为刚刚获得的植被覆盖率，计算之后可以得到地表比辐射率数据

黑体辐射亮度计算

在上文提到的网站去查询大气剖面数据(http://atmcorr.gsfc.nasa.gov/)，输入相关参数可得到

点击Calculate即可得到大气剖面信息

大气在热红外波段的透过率τ：0.60
大气向上辐射亮度L↑：3.22 W/(m2·sr·μm)
大气向下辐射亮辐射亮度L↓：5.03W/(m2·sr·μm)
在bandmath中输入表达式：(b2-3.22-0.60*(1-b1)5.03)/(0.60b1)
其中b2指的是最开始求的第10波段的辐射亮度值
b1指的是刚求的地表比辐射率

计算之后可以得到黑体辐射亮度影像

地表摄氏温度计算

在bandmath中输入表达式：(1321.08)/alog(774.89/b1 + 1) - 273

其中b1指的是黑体辐射亮度
计算之后即可得到地表温度数据了