保姆级教程：用ENVI+单窗算法反演Landsat8地表温度（附完整BandMath公式）
保姆级教程：ENVI+单窗算法实现Landsat8地表温度反演全流程解析
当我们需要从卫星影像中提取地表温度数据时，单窗算法因其相对简单的参数需求和较高的精度，成为许多研究者的首选方法。特别是对于Landsat8这样的开源数据，配合ENVI软件强大的BandMath功能，可以实现从原始数据到温度产品的完整处理流程。本文将手把手带你完成整个反演过程，重点解决实际操作中的公式输入和参数设置难题。
1. 数据准备与环境配置
在开始温度反演前，我们需要确保拥有完整的数据集和适当的软件环境。Landsat8数据可以从USGS EarthExplorer或地理空间数据云平台免费获取。下载时应选择云量低于10%的影像，并确保包含以下关键文件：
- MTL元数据文件：包含辐射定标参数
- B10波段：热红外波段（分辨率100米）
- B4-B5波段：用于NDVI计算（分辨率30米）
ENVI软件需要安装以下扩展模块：
- Radiometric Correction：用于辐射定标
- Atmospheric Correction：可选，用于反射率计算
- BandMath：核心计算工具
提示：建议在处理前创建新的ENVI工程目录，将原始数据、中间文件和最终结果分类存放，避免文件混乱。
2. 辐射与大气校正处理
2.1 辐射定标流程
辐射定标是将DN值转换为辐射亮度的关键步骤。在ENVI中操作如下：
- 打开
Radiometric Calibration工具 - 选择多光谱波段（B1-B7），设置参数：
定标类型：辐射亮度(Radiance) 缩放因子：0.1 输出数据类型：浮点型(Float) - 对热红外波段B10单独处理，参数设置：
定标类型：辐射亮度 缩放因子：0.0003342 偏移量：0.1
2.2 大气校正要点
虽然单窗算法对大气校正要求不高，但进行适当处理可提高NDVI计算精度：
# 示例：ENVI IDL大气校正脚本 pro atmospheric_correction input_file = dialog_pickfile(title='选择辐射定标后文件') output_file = input_file + '_atmcor' envi_doit, 'QUAC', $ input_file=input_file, $ output_file=output_file end3. 关键参数计算与实现
3.1 大气透射率(t)获取
大气透射率可通过以下方式确定：
| 方法 | 精度 | 数据需求 | 适用场景 |
|---|---|---|---|
| 在线计算工具 | 高 | 时间、位置、气象数据 | 精确研究 |
| 经验值 | 中 | 无 | 快速估算 |
| 历史数据拟合 | 中高 | 长期观测数据 | 区域研究 |
推荐使用NASA的AtmCorr工具在线计算，需要输入：
- 影像获取时间
- 经纬度坐标
- 当天天气数据（温度、湿度、气压）
3.2 植被指数与比辐射率计算
NDVI计算公式：
NDVI = \frac{B5 - B4}{B5 + B4}ENVI BandMath表达式：
(float(b5)-float(b4))/(float(b5)+float(b4))植被覆盖度(Fv)计算：
采用置信区间法确定NDVI极值：
NDVImin = 0.1 (5%置信区间) NDVImax = 0.57 (95%置信区间)BandMath表达式：
(b1 gt 0.57)*1 + (b1 lt 0.1)*0 + (b1 ge 0.1 and b1 le 0.57)*((b1-0.1)/(0.57-0.1))地表比辐射率(ε)分类计算：
- 水体：0.995
- 城镇：0.9589 + 0.086Fv - 0.0671Fv²
- 自然表面：0.9625 + 0.0614Fv - 0.0461Fv²
对应BandMath表达式：
(b1 le 0.1)*0.995 + (b1 gt 0.1 and b1 lt 0.57)*(0.9589+0.086*b2-0.0671*b2*b2) + (b1 ge 0.57)*(0.9625+0.0614*b2-0.0461*b2*b2)4. 单窗算法核心计算
4.1 中间变量计算
根据覃志豪单窗算法模型，需要先计算两个中间变量：
变量C：
C = ε × tBandMath实现：
0.73*b1 # 假设t=0.73，b1为比辐射率图变量D：
D = (1-t)[1+(1-ε)t]BandMath表达式：
(1-0.73)*(1+(1-b1)*0.73)4.2 亮温(T6)计算
从B10波段辐射亮度转换为亮温：
T6 = \frac{K2}{\ln(\frac{K1}{Lλ} + 1)}其中Landsat8常量：
K1 = 774.89 W/(m²·sr·μm) K2 = 1321.08 KBandMath表达式：
1321.08/alog(774.89/b1+1)4.3 大气平均作用温度(Ta)估算
根据研究区气候类型选择公式：
| 气候类型 | 公式 |
|---|---|
| 热带 | Ta = 17.9769 + 0.91715T0 |
| 中纬度夏季 | Ta = 16.0110 + 0.92621T0 |
| 中纬度冬季 | Ta = 19.2704 + 0.91118T0 |
注意：T0为当地气温（单位K），可从气象站获取并转换
4.4 最终地表温度计算
单窗算法完整公式：
Ts = \frac{a(1-C-D) + [b(1-C-D)+C+D]T6 - DT_a}{C}其中常量：
a = -67.355351 b = 0.458606对应BandMath表达式：
(-67.355351*(1-b1-b2)+(0.458606*(1-b1-b2)+b1+b2)*b3-b2*292.67455805)/b1变量对应关系：
b1 = C b2 = D b3 = T65. 结果后处理与验证
5.1 温度单位转换
计算结果为开尔文温度，转换为摄氏度：
温度°C = 温度K - 273.15BandMath表达式：
b1 - 273.155.2 结果验证方法
- 统计检验：检查温度值是否在合理范围内（通常地表温度在-20°C至50°C之间）
- 空间分布验证：比较水体、城市、植被区的温度差异是否符合预期
- 实地数据对比：如有气象站数据，可进行点对点验证
5.3 常见问题排查
异常高温/低温值：
- 检查辐射定标参数是否正确
- 确认大气透射率取值合理
- 验证比辐射率分类是否准确
空间分布异常：
- 检查NDVI计算是否正确
- 确认研究区边界裁剪准确
- 排除云污染像元
在处理宁夏中卫市的数据时，发现南部山区出现异常低温区，经检查是残留薄云影响。通过设置温度阈值（如<0°C视为无效值）可以有效过滤这类异常。