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

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

地表温度是研究城市热岛效应、农业干旱监测、气候变化等领域的重要参数。Landsat 8卫星搭载的热红外传感器TIRS（Thermal 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^-1）
• Qcal是影像的原始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_add

3.2 亮温计算

将辐射亮度转换为亮温（Brightness Temperature）需要使用Planck公式的反函数：

T = K2 / ln(K1 / Lλ + 1)

对于Landsat 8的B10和B11波段：

• B10: K1 = 774.89, K2 = 1321.08
• B11: K1 = 480.89, K2 = 1201.14

Python实现代码：

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）是另一个关键参数。我们可以通过以下步骤计算：

1. 计算NDVI（归一化植被指数）：

   NDVI = (B5 - B4) / (B5 + B4)

2. 计算植被覆盖度（FVC）：

   FVC = ((NDVI - NDVI_soil) / (NDVI_veg - NDVI_soil))^2

   通常取NDVI_soil=0.2，NDVI_veg=0.5

3. 计算比辐射率：

   ε = 0.004 * FVC + 0.986

Python实现代码：

def calculate_ndvi(band4, band5): """计算NDVI""" return (band5 - band4) / (band5 + band4 + 1e-10) # 加小量避免除零 def calculate_fvc(ndvi, ndvi_soil=0.2, ndvi_veg=0.5): """计算植被覆盖度""" return ((ndvi - ndvi_soil) / (ndvi_veg - ndvi_soil + 1e-10)) ** 2 def calculate_lse(fvc): """计算地表比辐射率""" return 0.004 * fvc + 0.986

5. 劈窗算法应用

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 - ε11

Python实现代码：

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 lst

6. 结果验证与优化

6.1 结果验证

获得地表温度后，我们需要验证结果的合理性。常用的方法包括：

1. 检查温度范围是否合理（通常在250-330K之间）
2. 比较水体温度（应该比较稳定，约300K左右）
3. 检查城市热岛效应是否明显（城市区域温度应高于周边）

6.2 常见问题处理

在实际项目中，我遇到过几个常见问题：

1. 异常值处理：有时会出现极端高温或低温值，可以使用中值滤波或设置合理的温度范围来过滤。

2. 云污染：云会严重影响温度反演结果，建议先进行云检测和掩膜处理。

3. 比辐射率估算误差：在裸土或冰雪覆盖区域，比辐射率估算可能不准确，需要特殊处理。

Python代码示例（异常值处理）：

def filter_outliers(lst, min_temp=250, max_temp=330): """ 过滤异常温度值 :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_filtered

7. 完整处理流程示例

下面是一个完整的处理流程示例代码：

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(dtype=rasterio.float32, nodata=np.nan) with rasterio.open('lst_result.tif', 'w', **profile) as dst: dst.write(lst_filtered.astype(np.float32), 1)

在实际项目中，你可能需要根据具体需求调整参数和处理流程。比如在干旱地区，可能需要调整比辐射率计算公式；在城市区域，可能需要考虑建筑材料的发射率特性。