从Landsat 8数据到地表温度：劈窗算法实战解析

1. 从Landsat 8数据到地表温度劈窗算法实战解析地表温度是研究城市热岛效应、农业干旱监测、气候变化等领域的重要参数。Landsat 8卫星搭载的热红外传感器TIRSThermal Infrared Sensor提供了两个热红外波段B10和B11非常适合用于地表温度反演。劈窗算法Split-Window Algorithm是一种经典的地表温度反演方法它利用两个相邻的热红外波段来消除大气水汽的影响从而获得更准确的地表温度。本文将带你一步步完成从Landsat 8原始数据到地表温度产品的完整流程。无论你是遥感专业的学生还是刚接触遥感数据处理的工程师都能跟着这个指南完成实际操作。我会分享我在项目中积累的经验包括常见问题的解决方法让你少走弯路。2. 数据准备与环境搭建2.1 获取Landsat 8数据首先你需要获取一景Landsat 8数据。推荐从USGS EarthExplorer网站下载选择包含TIRS波段B10和B11的数据产品。下载后你会得到一个压缩包解压后包含多个文件其中我们需要重点关注的是以_B10.TIF和_B11.TIF结尾的热红外波段数据以_MTL.txt结尾的元数据文件元数据文件非常重要它包含了辐射定标所需的参数比如RADIANCE_MULT_BAND_x和RADIANCE_ADD_BAND_x等。2.2 搭建Python处理环境我推荐使用Python来处理Landsat 8数据因为Python有丰富的遥感数据处理库。你需要安装以下库pip install numpy gdal matplotlib rasterio如果你处理大量数据建议使用conda环境管理工具。我在实际项目中发现使用conda可以避免很多依赖冲突问题。3. 辐射定标与亮温计算3.1 辐射定标辐射定标是将DN值Digital Number转换为辐射亮度的过程。Landsat 8的辐射定标公式为Lλ RADIANCE_MULT_BAND_x * Qcal RADIANCE_ADD_BAND_x其中Lλ是辐射亮度W·m^-2·sr^-1·μm^-1Qcal是影像的原始DN值RADIANCE_MULT_BAND_x和RADIANCE_ADD_BAND_x可以从MTL文件中获取以下是Python实现代码import numpy as np import rasterio def radiance_calibration(dn_band, rad_mult, rad_add): 辐射定标 :param dn_band: 输入DN值数组 :param rad_mult: RADIANCE_MULT_BAND参数 :param rad_add: RADIANCE_ADD_BAND参数 :return: 辐射亮度数组 return rad_mult * dn_band rad_add3.2 亮温计算将辐射亮度转换为亮温Brightness Temperature需要使用Planck公式的反函数T K2 / ln(K1 / Lλ 1)对于Landsat 8的B10和B11波段B10: K1 774.89, K2 1321.08B11: K1 480.89, K2 1201.14Python实现代码def brightness_temperature(radiance, k1, k2): 计算亮温 :param radiance: 辐射亮度 :param k1: Planck常数K1 :param k2: Planck常数K2 :return: 亮温数组开尔文 return k2 / np.log(k1 / radiance 1)4. 关键参数估算4.1 大气水汽含量估算大气水汽含量Column Water Vapor, CWV是劈窗算法中的重要参数。对于Landsat 8数据我们可以使用以下经验公式估算CWV 0.02 0.631 * (BT11 - BT10) - 0.113 * (BT11 - BT10)^2其中BT10和BT11分别是B10和B11波段的亮温。在实际应用中我发现这个经验公式在大多数情况下都能给出合理的结果但在极端湿润或干燥地区可能需要调整系数。4.2 地表比辐射率估算地表比辐射率Land Surface Emissivity, LSE是另一个关键参数。我们可以通过以下步骤计算计算NDVI归一化植被指数NDVI (B5 - B4) / (B5 B4)计算植被覆盖度FVCFVC ((NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil))^2通常取NDVI_soil0.2NDVI_veg0.5计算比辐射率ε 0.004 * FVC 0.986Python实现代码def calculate_ndvi(band4, band5): 计算NDVI return (band5 - band4) / (band5 band4 1e-10) # 加小量避免除零 def calculate_fvc(ndvi, ndvi_soil0.2, ndvi_veg0.5): 计算植被覆盖度 return ((ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil 1e-10)) ** 2 def calculate_lse(fvc): 计算地表比辐射率 return 0.004 * fvc 0.9865. 劈窗算法应用5.1 劈窗算法原理劈窗算法利用两个相邻热红外波段对大气水汽吸收的差异来消除大气影响。基本形式为LST BT10 A * (BT10 - BT11) B其中A和B是与大气水汽含量和地表比辐射率相关的系数。5.2 算法实现我推荐使用Sobrino等人提出的劈窗算法它在各种地表条件下都表现良好LST BT10 1.378 * (BT10 - BT11) 0.183 * (BT10 - BT11)^2 - 0.268 (54.30 - 2.238 * CWV) * (1 - ε) (-129.20 16.40 * CWV) * Δε其中ε (ε10 ε11)/2Δε ε10 - ε11Python实现代码def split_window_algorithm(bt10, bt11, epsilon10, epsilon11, cwv): 劈窗算法计算地表温度 :param bt10: B10波段亮温 :param bt11: B11波段亮温 :param epsilon10: B10波段比辐射率 :param epsilon11: B11波段比辐射率 :param cwv: 大气水汽含量 :return: 地表温度开尔文 delta_bt bt10 - bt11 epsilon_mean (epsilon10 epsilon11) / 2 delta_epsilon epsilon10 - epsilon11 lst (bt10 1.378 * delta_bt 0.183 * delta_bt**2 - 0.268 (54.30 - 2.238 * cwv) * (1 - epsilon_mean) (-129.20 16.40 * cwv) * delta_epsilon) return lst6. 结果验证与优化6.1 结果验证获得地表温度后我们需要验证结果的合理性。常用的方法包括检查温度范围是否合理通常在250-330K之间比较水体温度应该比较稳定约300K左右检查城市热岛效应是否明显城市区域温度应高于周边6.2 常见问题处理在实际项目中我遇到过几个常见问题异常值处理有时会出现极端高温或低温值可以使用中值滤波或设置合理的温度范围来过滤。云污染云会严重影响温度反演结果建议先进行云检测和掩膜处理。比辐射率估算误差在裸土或冰雪覆盖区域比辐射率估算可能不准确需要特殊处理。Python代码示例异常值处理def filter_outliers(lst, min_temp250, max_temp330): 过滤异常温度值 :param lst: 地表温度数组 :param min_temp: 最小合理温度 :param max_temp: 最大合理温度 :return: 过滤后的温度数组 lst_filtered np.copy(lst) lst_filtered[lst min_temp] np.nan lst_filtered[lst max_temp] np.nan return lst_filtered7. 完整处理流程示例下面是一个完整的处理流程示例代码import numpy as np import rasterio # 1. 读取数据 with rasterio.open(LC08_L1TP_123045_20200520_20200520_01_RT_B10.TIF) as src: b10 src.read(1) with rasterio.open(LC08_L1TP_123045_20200520_20200520_01_RT_B11.TIF) as src: b11 src.read(1) with rasterio.open(LC08_L1TP_123045_20200520_20200520_01_RT_B4.TIF) as src: b4 src.read(1) with rasterio.open(LC08_L1TP_123045_20200520_20200520_01_RT_B5.TIF) as src: b5 src.read(1) # 2. 辐射定标 rad_b10 radiance_calibration(b10, 3.3420E-04, 0.10000) rad_b11 radiance_calibration(b11, 3.3420E-04, 0.10000) # 3. 亮温计算 bt10 brightness_temperature(rad_b10, 774.89, 1321.08) bt11 brightness_temperature(rad_b11, 480.89, 1201.14) # 4. 计算NDVI和比辐射率 ndvi calculate_ndvi(b4, b5) fvc calculate_fvc(ndvi) epsilon10 calculate_lse(fvc) epsilon11 calculate_lse(fvc) # 简化处理实际可能不同 # 5. 计算大气水汽含量 cwv 0.02 0.631 * (bt11 - bt10) - 0.113 * (bt11 - bt10)**2 # 6. 应用劈窗算法 lst split_window_algorithm(bt10, bt11, epsilon10, epsilon11, cwv) # 7. 过滤异常值 lst_filtered filter_outliers(lst) # 8. 保存结果 profile src.profile profile.update(dtyperasterio.float32, nodatanp.nan) with rasterio.open(lst_result.tif, w, **profile) as dst: dst.write(lst_filtered.astype(np.float32), 1)在实际项目中你可能需要根据具体需求调整参数和处理流程。比如在干旱地区可能需要调整比辐射率计算公式在城市区域可能需要考虑建筑材料的发射率特性。