零基础实战用Landsat 8数据反演地表温度的完整流程解析第一次接触地表温度反演时看着满屏的公式和术语我盯着屏幕发呆了半小时——辐射定标、亮温转换、大气校正...每个词都认识连起来却像天书。直到导师扔给我一份批注过的代码才恍然大悟原来关键不是理解所有理论而是掌握正确的工具链和参数。本文将用最直白的语言带你走完从原始数据到温度图的完整流程无论你用ENVI、ArcGIS还是Python都能找到对应方案。1. 环境准备与数据获取1.1 Landsat 8数据下载与预处理登录USGS EarthExplorerhttps://earthexplorer.usgs.gov/搜索目标区域后选择Landsat 8 Collection 2 Level 1数据。关键点确保包含B1010.6-11.19μm和B1111.5-12.51μm这两个热红外波段下载时勾选MTL.txt元数据文件内含后续计算必需的辐射定标参数解压后文件结构应包含LC08_L1TP_123032_20220520_20220527_02_T1_MTL.txt LC08_L1TP_123032_20220520_20220527_02_T1_B10.TIF LC08_L1TP_123032_20220520_20220527_02_T1_B11.TIF ...其他可见光波段...1.2 工具选择对比根据你的技术栈选择合适工具工具优点缺点适用场景ENVI可视化操作内置遥感函数商业软件成本高快速验证、教学演示ArcGIS Pro空间分析工具链完整热红外处理功能有限需后续空间分析的项目Python灵活可定制免费需编程基础批量处理、算法研发提示初学者建议先用ENVI熟悉流程再用Python实现自动化。本文后续将提供三种工具的代码片段。2. 辐射定标与亮温计算2.1 辐射定标实战打开MTL文件找到以下关键参数示例值RADIANCE_MULT_BAND_10 3.3420E-04 RADIANCE_ADD_BAND_10 0.10000 RADIANCE_MULT_BAND_11 3.3420E-04 RADIANCE_ADD_BAND_11 0.10000辐射定标公式Lλ RADIANCE_MULT_BAND_X * DN RADIANCE_ADD_BAND_X其中DN是像元原始数字值。ENVI操作菜单栏选择Basic Tools Band Math输入表达式3.3420e-4 * b1 0.1b1对应B10波段Python代码import rasterio import numpy as np with rasterio.open(B10.TIF) as src: dn src.read(1) radiance 3.3420e-4 * dn 0.1 profile src.profile with rasterio.open(B10_radiance.TIF, w, **profile) as dst: dst.write(radiance, 1)2.2 亮温温度转换需要从MTL获取两个常量K1_CONSTANT_BAND_10 774.89 K2_CONSTANT_BAND_10 1321.08 K1_CONSTANT_BAND_11 480.89 K2_CONSTANT_BAND_11 1201.14亮温计算公式BT K2 / ln(K1/Lλ 1) - 273.15 # 转换为摄氏度常见错误排查单位混淆确认K1/K2值与波段对应无效值处理DN0的像元需设置为NaN温度范围异常正常地表温度应在-20℃~50℃之间3. 大气水汽含量估算3.1 简化估算方法对于Landsat 8可采用Jiménez-Muñoz(2014)提出的经验公式w 0.4 × (BT11 - BT10) 1.2其中w大气水汽含量g/cm²BT10/BT11B10/B11的亮温值注意该方法适用于中纬度地区夏季水汽含量2-3g/cm²冬季或极端气候需用MOD05数据辅助3.2 精度优化方案如需更高精度可结合ERA5再分析数据从Copernicus Climate Data Store下载同期数据提取total_column_water_vapour变量重采样至Landsat分辨率30mPython实现import xarray as xr era5 xr.open_dataset(era5.nc) w era5[tcwv].interp(latitudelat_range, longitudelon_range)4. 地表比辐射率计算4.1 NDVI与植被覆盖度首先计算归一化植被指数# B5: RED, B4: NIR ndvi (nir - red) / (nir red 1e-10) # 避免除零错误然后估算植被覆盖度FVCndvi_min, ndvi_max 0.2, 0.5 # 典型取值范围 fvc np.clip((ndvi - ndvi_min) / (ndvi_max - ndvi_min), 0, 1)4.2 比辐射率估算采用Sobrino(2008)的混合像元法ε10 0.973 * (1 - FVC) 0.987 * FVC # B10波段 ε11 0.961 * (1 - FVC) 0.980 * FVC # B11波段典型错误未做NDVI裁剪会导致FVC超出[0,1]范围波段混淆ε10对应B10ε11对应B11城市区域需额外考虑建筑材质影响5. 劈窗算法温度反演5.1 算法核心公式采用最成熟的Jiménez-Muñoz劈窗算法LST BT10 A0 A1*(BT10-BT11) A2*(BT10-BT11)^2 (A3 A4*w)*(1-ε) A5*Δε其中ε (ε10 ε11)/2Δε ε10 - ε11系数A0-A5随温度范围变化见下表温度区间A0A1A2A3A4A5-10~50℃-0.2681.3780.18354.300-2.238-129.205.2 全流程代码示例# 输入参数 bt10 ... # B10亮温 bt11 ... # B11亮温 w ... # 水汽含量 e10 ... # B10比辐射率 e11 ... # B11比辐射率 # 中间计算 e_mean (e10 e11) / 2 delta_e e10 - e11 # 劈窗算法 A [-0.268, 1.378, 0.183, 54.300, -2.238, -129.20] lst (bt10 A[0] A[1]*(bt10-bt11) A[2]*(bt10-bt11)**2 (A[3] A[4]*w)*(1-e_mean) A[5]*delta_e)5.3 结果验证技巧温度范围检查城市区域通常高于周边乡村3-5℃水体温度大型水体日间温度应稳定在15-25℃时间序列比对同区域不同时相结果应呈现合理变化6. 三种工具链完整实现6.1 ENVI操作流程辐射定标Basic Tools Band Math亮温计算同上输入温度转换公式波段运算依次计算NDVI、FVC、比辐射率最终温度用Band Math实现劈窗算法公式6.2 ArcGIS Pro方案使用Raster Calculator工具分步构建处理模型# 示例片段 - 亮温计算 arcpy.gp.RasterCalculator_sa( 1321.08 / Ln(774.89 / \%B10%\ 1) - 273.15, output_raster)6.3 Python自动化脚本推荐处理框架import rasterio, numpy as np def landsat_lst_pipeline(b10_path, b11_path, mtl_file): # 读取元数据 mtl parse_mtl(mtl_file) # 辐射定标 rad10 mtl[RADIANCE_MULT_B10] * b10 mtl[RADIANCE_ADD_B10] # 亮温转换 bt10 mtl[K2_B10] / np.log(mtl[K1_B10] / rad10 1) - 273.15 # ...完整处理流程... return lst7. 典型问题解决方案7.1 异常值处理条纹噪声使用ENVI Filter Destriping云污染结合QA波段掩膜边缘效应裁剪掉场景边缘10%区域7.2 精度提升技巧使用SRTM数据校正地形影响夜间数据可减少太阳辐射干扰城市区域引入建筑高度数据修正7.3 结果可视化优化import matplotlib.pyplot as plt plt.imshow(lst, cmapjet, vmin10, vmax40) plt.colorbar(labelTemperature (°C)) plt.savefig(lst_map.png, dpi300)记得第一次成功反演出温度图时我把结果和气象站数据对比——误差居然在1.5℃以内关键点在于严格检查每个中间步骤的输出范围比如辐射值不应出现负值NDVI必须在[-1,1]之间。当某个环节结果异常时就回溯到上一步排查。
手把手教你用Landsat 8数据反演地表温度：从辐射定标到劈窗算法的完整实操指南
零基础实战用Landsat 8数据反演地表温度的完整流程解析第一次接触地表温度反演时看着满屏的公式和术语我盯着屏幕发呆了半小时——辐射定标、亮温转换、大气校正...每个词都认识连起来却像天书。直到导师扔给我一份批注过的代码才恍然大悟原来关键不是理解所有理论而是掌握正确的工具链和参数。本文将用最直白的语言带你走完从原始数据到温度图的完整流程无论你用ENVI、ArcGIS还是Python都能找到对应方案。1. 环境准备与数据获取1.1 Landsat 8数据下载与预处理登录USGS EarthExplorerhttps://earthexplorer.usgs.gov/搜索目标区域后选择Landsat 8 Collection 2 Level 1数据。关键点确保包含B1010.6-11.19μm和B1111.5-12.51μm这两个热红外波段下载时勾选MTL.txt元数据文件内含后续计算必需的辐射定标参数解压后文件结构应包含LC08_L1TP_123032_20220520_20220527_02_T1_MTL.txt LC08_L1TP_123032_20220520_20220527_02_T1_B10.TIF LC08_L1TP_123032_20220520_20220527_02_T1_B11.TIF ...其他可见光波段...1.2 工具选择对比根据你的技术栈选择合适工具工具优点缺点适用场景ENVI可视化操作内置遥感函数商业软件成本高快速验证、教学演示ArcGIS Pro空间分析工具链完整热红外处理功能有限需后续空间分析的项目Python灵活可定制免费需编程基础批量处理、算法研发提示初学者建议先用ENVI熟悉流程再用Python实现自动化。本文后续将提供三种工具的代码片段。2. 辐射定标与亮温计算2.1 辐射定标实战打开MTL文件找到以下关键参数示例值RADIANCE_MULT_BAND_10 3.3420E-04 RADIANCE_ADD_BAND_10 0.10000 RADIANCE_MULT_BAND_11 3.3420E-04 RADIANCE_ADD_BAND_11 0.10000辐射定标公式Lλ RADIANCE_MULT_BAND_X * DN RADIANCE_ADD_BAND_X其中DN是像元原始数字值。ENVI操作菜单栏选择Basic Tools Band Math输入表达式3.3420e-4 * b1 0.1b1对应B10波段Python代码import rasterio import numpy as np with rasterio.open(B10.TIF) as src: dn src.read(1) radiance 3.3420e-4 * dn 0.1 profile src.profile with rasterio.open(B10_radiance.TIF, w, **profile) as dst: dst.write(radiance, 1)2.2 亮温温度转换需要从MTL获取两个常量K1_CONSTANT_BAND_10 774.89 K2_CONSTANT_BAND_10 1321.08 K1_CONSTANT_BAND_11 480.89 K2_CONSTANT_BAND_11 1201.14亮温计算公式BT K2 / ln(K1/Lλ 1) - 273.15 # 转换为摄氏度常见错误排查单位混淆确认K1/K2值与波段对应无效值处理DN0的像元需设置为NaN温度范围异常正常地表温度应在-20℃~50℃之间3. 大气水汽含量估算3.1 简化估算方法对于Landsat 8可采用Jiménez-Muñoz(2014)提出的经验公式w 0.4 × (BT11 - BT10) 1.2其中w大气水汽含量g/cm²BT10/BT11B10/B11的亮温值注意该方法适用于中纬度地区夏季水汽含量2-3g/cm²冬季或极端气候需用MOD05数据辅助3.2 精度优化方案如需更高精度可结合ERA5再分析数据从Copernicus Climate Data Store下载同期数据提取total_column_water_vapour变量重采样至Landsat分辨率30mPython实现import xarray as xr era5 xr.open_dataset(era5.nc) w era5[tcwv].interp(latitudelat_range, longitudelon_range)4. 地表比辐射率计算4.1 NDVI与植被覆盖度首先计算归一化植被指数# B5: RED, B4: NIR ndvi (nir - red) / (nir red 1e-10) # 避免除零错误然后估算植被覆盖度FVCndvi_min, ndvi_max 0.2, 0.5 # 典型取值范围 fvc np.clip((ndvi - ndvi_min) / (ndvi_max - ndvi_min), 0, 1)4.2 比辐射率估算采用Sobrino(2008)的混合像元法ε10 0.973 * (1 - FVC) 0.987 * FVC # B10波段 ε11 0.961 * (1 - FVC) 0.980 * FVC # B11波段典型错误未做NDVI裁剪会导致FVC超出[0,1]范围波段混淆ε10对应B10ε11对应B11城市区域需额外考虑建筑材质影响5. 劈窗算法温度反演5.1 算法核心公式采用最成熟的Jiménez-Muñoz劈窗算法LST BT10 A0 A1*(BT10-BT11) A2*(BT10-BT11)^2 (A3 A4*w)*(1-ε) A5*Δε其中ε (ε10 ε11)/2Δε ε10 - ε11系数A0-A5随温度范围变化见下表温度区间A0A1A2A3A4A5-10~50℃-0.2681.3780.18354.300-2.238-129.205.2 全流程代码示例# 输入参数 bt10 ... # B10亮温 bt11 ... # B11亮温 w ... # 水汽含量 e10 ... # B10比辐射率 e11 ... # B11比辐射率 # 中间计算 e_mean (e10 e11) / 2 delta_e e10 - e11 # 劈窗算法 A [-0.268, 1.378, 0.183, 54.300, -2.238, -129.20] lst (bt10 A[0] A[1]*(bt10-bt11) A[2]*(bt10-bt11)**2 (A[3] A[4]*w)*(1-e_mean) A[5]*delta_e)5.3 结果验证技巧温度范围检查城市区域通常高于周边乡村3-5℃水体温度大型水体日间温度应稳定在15-25℃时间序列比对同区域不同时相结果应呈现合理变化6. 三种工具链完整实现6.1 ENVI操作流程辐射定标Basic Tools Band Math亮温计算同上输入温度转换公式波段运算依次计算NDVI、FVC、比辐射率最终温度用Band Math实现劈窗算法公式6.2 ArcGIS Pro方案使用Raster Calculator工具分步构建处理模型# 示例片段 - 亮温计算 arcpy.gp.RasterCalculator_sa( 1321.08 / Ln(774.89 / \%B10%\ 1) - 273.15, output_raster)6.3 Python自动化脚本推荐处理框架import rasterio, numpy as np def landsat_lst_pipeline(b10_path, b11_path, mtl_file): # 读取元数据 mtl parse_mtl(mtl_file) # 辐射定标 rad10 mtl[RADIANCE_MULT_B10] * b10 mtl[RADIANCE_ADD_B10] # 亮温转换 bt10 mtl[K2_B10] / np.log(mtl[K1_B10] / rad10 1) - 273.15 # ...完整处理流程... return lst7. 典型问题解决方案7.1 异常值处理条纹噪声使用ENVI Filter Destriping云污染结合QA波段掩膜边缘效应裁剪掉场景边缘10%区域7.2 精度提升技巧使用SRTM数据校正地形影响夜间数据可减少太阳辐射干扰城市区域引入建筑高度数据修正7.3 结果可视化优化import matplotlib.pyplot as plt plt.imshow(lst, cmapjet, vmin10, vmax40) plt.colorbar(labelTemperature (°C)) plt.savefig(lst_map.png, dpi300)记得第一次成功反演出温度图时我把结果和气象站数据对比——误差居然在1.5℃以内关键点在于严格检查每个中间步骤的输出范围比如辐射值不应出现负值NDVI必须在[-1,1]之间。当某个环节结果异常时就回溯到上一步排查。
