Landsat8 TIRS 地表温度反演

热红外遥感（Infrared Remote Sensing）是指传感器工作波段限于红外波段范围之内的遥感。即利用星载或机载传感器收集、记录地物的热红外信息，并利用这种热红外信息来识别地物和反演地表参数如温度、湿度和热惯量等。目前有很多的卫星携带了热红外传感器，包括ASTER、AVHRR、MODIS、

TM/ETM+/ TIRS等。

目前，地表温度反演算法主要有以下三种：大气校正法（也称为辐射传输方程：Radiative Transfer Equation——RTE）、单通道算法和分裂窗算法。

本实例是基于大气校正法，利用Landsat8 TIRS反演地表温度。

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

具体实现为：卫星传感器接收到的热红外辐射亮度值Lλ由三部分组成：大气向上辐射亮度L↑，地面的真实辐射亮度经过大气层之后到达卫星传感器的能量；大气向下辐射到达地面后反射的能量。卫星传感器接收到的热红外辐射亮度值Lλ的表达式可写为（辐射传输方程）：

Lλ = [εB(T S) + (1-ε)L↓]τ+ L↑（1.1）

式中，ε为地表比辐射率，T S为地表真实温度(K)，B(T S)为黑体热辐射亮度，τ为大气在热红外波段的透过率。则温度为T的黑体在热红外波段的辐射亮度

B(T S)为：

B(T S) = [Lλ - L↑- τ(1-ε)L↓]/τε（1.2）

T s可以用普朗克公式的函数获取。

T S = K2/ln(K1/ B(T S)+ 1) （1.3）

对于TM，K1 =607.76 W/(m2*µm*sr)，K2 =1260.56K。

对于ETM+，K1=666.09 W/(m2*µm*sr)，K2 =1282.71K。

对于TIRS Band10，K1= 774.89 W/(m2*µm*sr)，K2 = 1321.08K。

从上可知此类算法需要2个参数：大气剖面参数和地表比辐射率。大气剖面参数在NASA提供的网站(/)中，输入成影时间以及中心经纬度可以获取大气剖面参数。适用于只有一个热红外波段的数据，如Landsat TM /ETM+/TIRS数据。

主要内容就是使用BandMath工具计算公式（1.2）和公式（1.3），处理流程如下图所示。

注：OLI大气校正是可选项，根据前人研究成果，大气校正处理对结果影响不是很大。

图1.1 基于大气校正法的TIRS反演流程图

1、图像辐射定标

(1)在主界面中，选择File→Open，在文件选择对话框中选择

"LC81230322013276LGN00_MTL.txt"文件，ENVI自动按照波长分为五个数据集：多光谱数据（1-7波段），全色波段数据（8波段），卷云波段数据（9波段），热红外数据（10，11波段）和质量波段数据（12波段）。

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

（3）在Radiometric Calibration面板中，设置以下参数：

•定标类型（Calibration Type）：辐射亮度值（radiance）。

•其他选择默认参数。

（4）选择输出路径和文件名1-LC81230322_band10_rad.dat，单击OK 按钮执行定标处理。得到Band10辐射亮度图像。

2、地表比辐射率计算

TIRS的Band10热红外波段与TM/ETM+6热红外波段具有近似的波谱范围，本例采用TM/ETM+6相同的地表比辐射率计算方法。使用Sobrino提出的NDVI阈值法计算地表比辐射率。

ε=0.004Pv+0.986 （1.4）

其中，Pv是植被覆盖度，用以下公式计算：

Pv = [(NDVI- NDVI Soil)／(NDVI Veg - NDVI Soil)] （1.5）

其中，NDVI为归一化植被指数，NDVI Soil为完全是裸土或无植被覆盖区域的NDVI值，NDVI Veg则代表完全被植被所覆盖的像元的NDVI值，即纯植被像元的NDVI值。取经验值NDVI Veg = 0.70和NDVI Soil = 0.05，即当某个像元的NDVI大于0.70时，Pv取值为1；当NDVI小于0.05，Pv取值为0。

注：这里采用简化的植被覆盖度计算模型，感兴趣的可以使用更加精确的植被覆盖度计算模型。

（1）在Toolbox工具箱中，双击Spectral/Vegetation/NDVI工具，在文件输入对话框中，选择Landsat8 OLI多光谱图像。

图1.2 NDVI文件输入对话框

提示：覃志豪提出使用原始的DN值图像计算NDVI对反演结果影响不大。

（2）在NDVI Calculaton parameters对话框中，自动识别NDVI计算波段：Red：4，Near IR：5。

图1.3 NDVI对话框

（3）选择输出文件名和路径。

（4）在Toobox中，选择Band Ratio/Band Math，输入表达式：(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。

计算得到植被覆盖度图像。

（5）在Toobox中，选择Band Ratio/Band Math，输入表达式：

0.004*b1+0.986

其中，b1：植被覆盖度图像。

计算得到地表比辐射率图像。