地表温度反演（Land Surface Temperature, LST）是遥感应用中的一个重要任务，Landsat 8 卫星提供了热红外波段（TIRS）数据，可以用于计算地表温度。实现 Landsat 8 地表温度反演通常包括以下几个关键步骤：辐射定标、大气校正、比辐射率计算和温度反演算法（如单窗算法、分裂窗算法等）。以下是基于 Python 的 Landsat 8 数据进行地表温度反演的实现方法和代码示例。 ### 数据准备 Landsat 8 数据通常包括多个波段，其中用于地表温度反演的主要波段是 **热红外波段 10（TIRS Band 10）** 和 **近红外波段 5（NIR Band 5）**（用于计算 NDVI 和比辐射率）。 ### 1. 辐射定标 将数字数值（DN）转换为表观辐射亮度（Top-of-Atmosphere Radiance）： $$ L_{\lambda} = M_L \cdot Q_{cal} + A_L $$ 其中 $ M_L $ 和 $ A_L $ 是元数据中的增益和偏移值。 ```python import numpy as as np from osgeo import gdal # 读取热红外波段 10 thermal_band = "path_to_band10.tif" thermal_ds = gdal.Open(thermal_band) thermal_arr = thermal_ds.ReadAsArray().astype(np.float32) # 获取增益和偏移值（来自 MTL 文件） ML = 0.0003342 AL = 0.1 # 辐射定标 radiance = ML * thermal_arr + AL ``` ### 2. 亮温计算（使用普朗克函数） $$ T = \frac{K_2}{\ln\left(\frac{K_1}{L_{\lambda}} + 1\right)} $$ 其中 $ K_1 $ 和 $ K_2 $ 是传感器特定常数，对于 Band 10，分别为 774.89 和 1321.08。 ```python K1 = 774.89 K2 = 1321.08 # 计算亮温 brightness_temp = K2 / np.log((K1 / radiance) + 1) ``` ### 3. 比辐射率计算 使用 NDVI 计算比辐射率（Emissivity）： ```python # 读取近红外波段 5 和红光波段 4 nir_band = "path_to_band5.tif" red_band = "path_to_band4.tif" nir_ds = gdal.Open(nir_band) red_ds = gdal.Open(red_band) nir_arr = nir_ds.ReadAsArray().astype(np.float32) red_arr = red_ds.ReadAsArray().astype(np.float32) # 计算 NDVI ndvi = (nir_arr - red_arr) / (nir_arr + red_arr + 1e-10) # 比辐射率计算（经验公式） emissivity = np.where((ndvi < 0.2), 0.99, np.where((0.2 <= ndvi) & (ndvi <= 0.5), 0.98, 0.97)) ``` ### 4. 地表温度反演（单窗算法） 使用单窗算法对亮温进行修正： $$ LST = \frac{T_{b}}{1 + (\lambda \cdot T_{b} / \rho) \cdot \ln(\varepsilon)} $$ 其中 $ \lambda $ 是波长（约为 10.8 μm），$ \rho = 1.4388 \times 10^{-2} \, \text{m·K} $，$ \varepsilon $ 是比辐射率。 ```python # 常量 rho = 1.4388e-2 lambda_band = 10.8 # 热红外波段中心波长（μm） # 地表温度计算 lst = brightness_temp / (1 + (lambda_band * brightness_temp / rho) * np.log(emissivity)) ``` ### 5. 保存结果 ```python driver = gdal.GetDriverByName("GTiff") lst_ds = driver.Create("lst_output.tif", lst.shape[1], lst.shape[0], 1, gdal.GDT_Float32) lst_ds.SetGeoTransform(thermal_ds.GetGeoTransform()) lst_ds.SetProjection(thermal_ds.GetProjection()) band = lst_ds.GetRasterBand(1) band.WriteArray(lst) lst_ds.FlushCache() ``` ### 注意事项 - 实际应用中可能需要考虑大气校正（如使用大气辐射传输模型如 MODTRAN 或 6S）。 - 比辐射率的计算可根据具体研究区域的地表覆盖类型进行更精细的调整。 - 代码中路径和常量需根据实际数据和元数据进行修改。 ---