保姆级教程：用ENVI+Erdas从Landsat数据反演地表温度（附完整模型与避坑指南）
从Landsat数据到地表温度：ENVI+Erdas全流程实战指南
当第一次接触遥感温度反演时，很多人会被复杂的公式和繁琐的步骤吓退。但事实上，只要掌握正确的工具链和清晰的流程，即使是初学者也能顺利完成从原始Landsat数据到地表温度产品的全流程处理。本文将带你一步步走过这个旅程，避开那些容易让人"掉坑"的陷阱。
1. 准备工作与环境配置
在开始温度反演之前，确保你的工作环境已经准备就绪。ENVI 5.3+和Erdas Imagine 2014+是最佳组合，它们对Landsat数据的支持最为完善。建议创建一个专门的项目文件夹，按照以下结构组织：
温度反演项目/ ├── 原始数据/ │ └── LC08_L1TP_123032_20200520_20200520_01_RT.tar ├── 中间结果/ │ ├── 辐射定标/ │ ├── 大气校正/ │ └── 植被指数/ └── 最终结果/关键检查点：
- 确认Landsat数据级别为L1TP（地形校正级）
- 检查影像覆盖区域是否有云层遮挡（云量<10%为佳）
- 准备足够磁盘空间（完整流程需要原始数据3-5倍的存储空间）
提示：从USGS官网下载数据时，优先选择Collection 2级别的数据，其辐射定标参数更为精确。
2. 数据预处理：从DN值到物理量
温度反演的精度很大程度上取决于预处理的质量。这个阶段需要完成三个关键转换：DN值→辐射亮度→反射率→地表反射率。
2.1 辐射定标：将数字信号转换为辐射亮度
在ENVI中操作步骤：
- 打开
Radiometric Correction→Radiometric Calibration - 选择热红外波段（B10或B11，取决于Landsat型号）
- 设置输出为
Radiance，单位W/(m²·sr·μm) - 关键参数核对：
- 太阳高度角（自动从元数据读取）
- 辐射定标系数（确保与数据版本匹配）
# 示例：Landsat 8 Band 10辐射定标公式 Lλ = ML * Qcal + AL 其中： ML = 3.3420E-04 (radiance multiplicative factor) AL = 0.10000 (radiance additive factor) Qcal = 原始DN值2.2 大气校正：消除大气干扰
推荐使用ENVI的FLAASH模块：
- 输入刚生成的辐射亮度图像
- 设置传感器类型为对应Landsat型号
- 大气模型根据影像经纬度和季节选择
- 气溶胶模型优先选择
Rural - 运行前务必检查：
- 地面高程（DEM数据可提升精度）
- 水汽反演选项（对温度反演至关重要）
常见问题排查：
- 如果结果出现异常值，尝试调整气溶胶光学厚度初值
- 边缘区域出现条带？检查影像镶嵌时的重叠区处理
3. 温度反演核心算法实现
3.1 亮度温度计算
在Erdas Modeler中构建第一个关键模型：
亮度温度模型： 输入：辐射定标后的热红外波段 处理： 1. 应用Planck公式反演： Tb = K2 / ln(K1/Lλ + 1) 其中： - Landsat 8 B10: K1=774.89, K2=1321.08 - Landsat 8 B11: K1=480.89, K2=1201.14 2. 设置输出为Float64类型 3. 勾选"Ignore Zero"选项注意：不同Landsat卫星的热红外波段参数差异很大，务必核对元数据中的具体参数。
3.2 植被覆盖度计算
基于NDVI的植被覆盖度(FVC)计算需要特别注意阈值处理：
当 NDVI < 0.05: FVC = 0 当 NDVI > 0.7: FVC = 1 否则: FVC = ((NDVI - 0.05) / (0.7 - 0.05))^2在Erdas中实现时，建议使用Conditional工具分步处理，避免公式一次性输入导致的语法错误。
3.3 地表比辐射率估算
根据地物类型采用不同计算方法：
| 地物类型 | 计算公式 | 备注 |
|---|---|---|
| 水体 | 0.991 | 需结合监督分类结果 |
| 植被 | 0.986 | 与FVC线性相关 |
| 城市 | 0.970 | 考虑建筑密度影响 |
实际操作技巧：
- 先完成监督分类（至少区分水体、植被、城镇三类）
- 对分类结果进行多数滤波（3×3窗口）消除噪声
- 使用
Reclassify工具为每类赋予相应比辐射率初值 - 对植被区域应用FVC加权公式
4. 完整温度反演模型构建
将所有组件集成到Erdas Modeler中，建议采用模块化设计：
主模型流程图： [输入热红外波段] → [亮度温度计算] ↓ [NDVI数据] → [植被覆盖度] → [比辐射率计算] ↓ ↓ [大气水汽数据] → [大气参数计算] → [地表温度反演]关键参数设置表：
| 参数项 | 推荐值 | 可调范围 | 影响分析 |
|---|---|---|---|
| 大气透射率τ | 0.8-0.9 | ±0.05 | 每0.1偏差导致1-2K误差 |
| 大气水汽含量w | 1.5-3.0 g/cm² | 根据气候区调整 | 湿润地区需现场测量 |
| 比辐射率ε | 0.96-0.99 | 地物依赖 | 主要误差来源之一 |
模型调试技巧：
- 分阶段验证：每个子模型单独运行检查
- 使用
View Histogram工具检查结果分布 - 对异常值区域进行局部采样分析
- 保存中间结果以便回溯问题
5. 结果验证与精度提升
5.1 交叉验证方法
- 空间一致性检查：
- 对比水体/植被/城市区域的温度差异合理性
- 检查山体阴/阳坡温度分布规律
- 时间序列分析：
- 同日不同过境时间的数据对比
- 同区域历史数据趋势验证
- 地面实测对比：
- 利用气象站数据校正系统偏差
- 移动测温设备现场验证（需同步观测）
5.2 常见问题解决方案
问题1：结果图像出现条带状异常
- 可能原因：辐射定标系数错误
- 解决方案：重新下载元数据文件核对参数
问题2：城市区域温度明显偏低
- 可能原因：比辐射率设置过高
- 调整建议：将城市区域ε值降低0.01-0.02
问题3：模型运行报"Divide by zero"错误
- 排查步骤：
- 检查所有除法运算是否有零值输入保护
- 确认"Ignore Zero"选项已勾选
- 对输入数据运行
Statistics查看极值
6. 高级技巧与效率优化
6.1 批处理实现
对于长时间序列分析，建议创建批处理脚本：
# 示例：Linux下使用GDAL批量处理 for file in LS08_*_B10.TIF; do gdal_calc.py -A $file --outfile=rad_${file} \ --calc="3.3420E-04*A+0.10000" --type=Float32 gdal_calc.py -A rad_${file} --outfile=BT_${file} \ --calc="1321.08/log(774.89/A+1)" --type=Float32 done6.2 精度提升方法
- DEM融合：
- 使用30m SRTM数据校正地形效应
- 实现高程相关的温度递减率调整
- 多算法集成：
- 单窗算法与分裂窗算法结果加权平均
- 不同大气校正结果交叉验证
- 机器学习后处理：
- 训练随机森林模型校正系统误差
- 应用CNN去除异常噪声
在实际项目中，我们发现最耗时的步骤往往是大气校正而非温度计算本身。一个实用的建议是：对静态区域可以预先建立大气参数查找表，将处理时间缩短60%以上。