地表温度反演操作总结

由于博主也是自己边做实验便总结，难免会有混乱出错的地方，欢迎大家指正~

整体实验流程：

1 数据预处理

1.1 数据下载

本实验使用的是Landsat8 OLI_TIRS数据，在地理空间数据云（http://www.gscloud.cn/）下载三幅广州地区的影像预备进行预处理。（由于这仅为一个练手实验，所以并未对云量有太高要求）

1.2 辐射定标

文件解压后使用File-Open as-Landsat-GeoTIFF with Metadata打开XML文件就可打开所需的多光谱和全色影像。
使用Radiometric Correction-Radiometric Calibration工具进行辐射定标，设定如下图（注意：辐射定标时只能对多光谱影像进行定标）

1.3 大气校正（温度反演可省去此步，对结果影响不大）

大气校正模型有两种，这里我们使用FLAASH大气校正，打开Radiometric Correction-Atmospheric Correcrion Module-FLAASH Atmospheric Correction工具，打开已定标的多光谱图像

1.4 图像融合

使用Image Sharpening-Gam-Schmidt Pan Sharpening工具进行图像融合，先选择低分辨率的多光谱影像（MultiSpectral）后选择高分辨率的全色影像（Panchromatic）

选择传感器类型和重采样方法，设置输出文件夹，点击OK输出融合图像

1.5 图像拼接

将三幅完成预处理的影像使用 Mosaicking-Seamless Mosaic工具进行拼接
点击加号添加三幅影像：

Color Matching Action可选择影像颜色的被参照图像和校正图像，选择一副图像为被参照图，另外两副的颜色就会校正成与被参照图相似的颜色

勾选Histogram Matching，这里主要是选择颜色改变是交界处改变还是全图改变

选择Auto Generate Seamlines自动生成缝纫线，能观察到如何裁剪图像

设置输出背景，输出文件夹和格式，点击Finish便可得到拼接文件

注：可勾选右上角Show Preview观察每一步操作后的图像

1.6 图像裁剪

现在我们得到三幅影像拼接的图像，我们要准确的切割出广州市的影像
打开广州市边界矢量图，使用Regions of Interest-Vector to ROI工具将矢量图转化为ROI方便切割（注意：若是两个坐标系不一致，还需要进行投影变换）

使用Subset Data from ROIs工具进行裁剪
先选择拼接好的图像

再选择广州市ROI，选择输出文件夹，点击OK输出裁剪图像（注意：Mask pixels output of ROI需要选YES才能输出精确广州市边界的图像，若选NO则输出一个矩形）

1.7 Band10辐射亮度图像获取

Band10波段后缀为Thermal，打开后照以上步骤拼接剪裁，然后进行辐射定标

打开Radiometric Calibration，选中Thermal图像，选择Thermal Infrared 1（10.9）光谱

然后开始进行辐射定标，设定如下

2 地表温度反演

地表温度反演算法主要有以下三种：大气校正法（也称为辐射传输方程：Radiative Transfer Equation——RTE）、单通道算法和分裂窗算法。
本实验采用大气校正法。
大气校正法基本原理：首先估计大气对地表热辐射的影响, 然后把这部分大气影响从卫星传感器所观测到的热辐射总量中减去, 从而得到地表热辐射强度, 再把这一热辐射强度转化为相应的地表温度。
卫星传感器接收到的热红外辐射亮度值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。

2.1 NDVI指数运算

使用Band Math计算NDVI指数(或者可以直接使用ENVI自带的计算NDVI指数工具)，指数公式:(近红外-红)/(近红外+红)
Band Math计算结果输出如下：

2.2 计算植被覆盖度

计算公式：(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输出结果
Band Math计算结果输出如下：

2.3 计算地表比辐射率

计算公式：0.004*B1+0.986
B1选择计算植被覆盖度后的结果
Band Math计算结果输出如下：

2.4计算同温度下的黑体辐射亮度图像

计算公式：(B2-0.75-0.9*(1-B1)1.29)/(0.9B1)
B1选择地表比辐射率图像
B2选Band10辐射亮度图像
Band Math计算结果输出如下：

2.5计算地表温度

计算公式：(1321.08)/alog(774.89/B1+1)-273
B1选择黑体辐射亮度图像
Band Math计算结果输出如下：

使用Cursor Value查看地表温度的数值，大多是在8-15℃的范围内，结合影像拍摄时间为凌晨两点，温度范围合理。