基于 Google Earth Engine 的北京市中心区域地表温度反演

一、引言

地表温度（Land Surface Temperature, LST）是地球表面能量平衡和气候变化研究中的关键参数。利用卫星遥感数据结合相关算法，可以实现对大面积地表温度的快速获取。Google Earth Engine（GEE）提供了丰富的卫星影像资源和强大的计算能力，使得地表温度反演变得更加高效和便捷。本文将详细介绍如何使用 GEE 对北京市中心区域的地表温度进行反演，并将结果导出保存。

二、代码实现

2.1 定义区域和时间范围

// ------------------------- 步骤1：定义区域和时间范围 ------------------------- //
var roi = ee.Geometry.Point([116.3975, 39.9087]).buffer(10000); // 示例：北京市中心10公里缓冲区
Map.centerObject(roi, 10);  // 在地图中心显示区域
Map.addLayer(roi, {color:'red'}, 'ROI');  // 可视化ROI

var startDate = '2023-05-01';
var endDate = '2023-05-31';

这段代码定义了一个以北京市中心坐标为中心，半径为 10 公里的缓冲区作为研究区域（ROI），并将其显示在地图上。同时，设置了数据获取的时间范围为 2023 年 5 月 1 日至 5 月 31 日。

2.2 加载并预处理 Landsat 8 数据

// ------------------------- 步骤2：加载并预处理Landsat 8数据 ------------------------- //
var landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filterBounds(roi)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUD_COVER', 20))  // 过滤云量>20%的图像
  .first();  // 选择时间范围内的第一景影像

// 打印影像信息以确认波段名称
print('原始影像信息:', landsat);

这里从 GEE 的影像集合中加载 Landsat 8 数据，通过过滤条件筛选出位于研究区域内、在指定时间范围内且云量小于 20% 的影像，并选取其中的第一景影像进行后续处理。同时，打印出影像的信息以便确认波段名称。

2.3 计算 NDVI（归一化植被指数）

// ------------------------- 步骤3：计算NDVI（归一化植被指数） ------------------------- //
var ndvi = landsat.normalizedDifference(['SR_B5', 'SR_B4'])  // Landsat 8的SR_B5=近红外，SR_B4=红光
  .rename('NDVI')
  .clip(roi);  // 裁剪到ROI范围

// 可视化NDVI（范围0~1）
Map.addLayer(ndvi, {min: 0, max: 1, palette: ['red', 'yellow', 'green']}, 'NDVI');

使用 Landsat 8 影像的近红外波段（SR_B5）和红光波段（SR_B4）计算归一化植被指数（NDVI），将结果重命名为 “NDVI” 并裁剪到研究区域范围内，然后将其可视化显示在地图上。

2.4 计算地表比辐射率（Emissivity）

// ------------------------- 步骤4：计算地表比辐射率（Emissivity） ------------------------- //
// 计算植被覆盖度（FVC），假设NDVI在0.2（裸土）到0.5（茂密植被）之间
var fvc = ndvi.expression(
  '((NDVI - 0.2) / (0.5 - 0.2)) ** 2',  // 经验公式
  {'NDVI': ndvi}
).clamp(0, 1).rename('FVC');  // 限制范围在0~1

// 计算比辐射率（ε = FVC*ε_veg + (1-FVC)*ε_soil）
var emissivity = fvc.multiply(0.99)  // 植被比辐射率=0.99
  .add(fvc.multiply(-1).add(1).multiply(0.97))  // 裸土比辐射率=0.97
  .rename('Emissivity');

先根据 NDVI 值通过经验公式计算植被覆盖度（FVC），并将其限制在 0 到 1 的范围内。然后，根据植被覆盖度计算地表比辐射率（Emissivity），采用植被比辐射率为 0.99，裸土比辐射率为 0.97 的参数进行计算。

2.5 计算亮度温度（Brightness Temperature）

// ------------------------- 步骤5：计算亮度温度（Brightness Temperature） ------------------------- //
// Landsat 8 Band 10的辐射亮度转换为温度（单位：摄氏度）
var thermalBand = landsat.select('ST_B10');  // ST_B10是地表温度波段（单位：开尔文）
var bt = thermalBand.multiply(0.00341802).add(149.0)  // 转换为辐射亮度（原始单位为K，需校准）
  .subtract(273.15)  // 开尔文转摄氏度
  .rename('BT')
  .clip(roi);

// 可视化亮度温度（示例范围20~40°C）
Map.addLayer(bt, {min: 20, max: 40, palette: ['blue', 'yellow','red']}, 'Brightness Temp');

从 Landsat 8 影像中选择地表温度波段（ST_B10），通过一系列的转换公式将其从辐射亮度转换为摄氏度表示的亮度温度（BT），并裁剪到研究区域范围内，最后将亮度温度可视化显示在地图上。

2.6 反演地表温度（LST）

// ------------------------- 步骤6：反演地表温度（LST） ------------------------- //
// 使用单通道算法（Jiménez-Muñoz et al., 2014）参数
var a = 0.866;  // 大气透过率（需根据实际大气数据调整）
var b = 0.133;  // 经验参数

var lst = bt.expression(
  '(BT / (1 + (a * BT / 14388) * log(Emissivity))) - b',  // 公式来源：Jiménez-Muñoz (2014)
  {
    'BT': bt,  // 亮度温度
    'Emissivity': emissivity,  // 比辐射率
    'a': a,
    'b': b
  }
).rename('LST');  // 单位：摄氏度

// 可视化地表温度（可根据研究区域调整范围）
Map.addLayer(lst, {min: 25, max: 45, palette: ['blue', 'yellow','red']}, 'Land Surface Temp');

采用单通道算法（Jiménez-Muñoz et al., 2014），根据之前计算得到的亮度温度（BT）和地表比辐射率（Emissivity），结合设定的大气透过率（a）和经验参数（b），反演得到地表温度（LST），并将其可视化显示在地图上。

2.7 导出结果到 Google Drive

// ------------------------- 步骤7：导出结果到Google Drive ------------------------- //
Export.image.toDrive({
  image: lst,
  description: 'LST_Beijing_202305',  // 导出文件名
  folder: 'GEE_Exports',  // 指定Google Drive文件夹
  region: roi,
  scale: 30,  // Landsat分辨率（单位：米）
  maxPixels: 1e13,  // 允许的最大像素数量
  fileFormat: 'GeoTIFF'  // 导出格式
});

// 打印成功信息
print('代码执行完成！请在Google Drive的GEE_Exports文件夹中查看结果。');

将反演得到的地表温度影像导出到 Google Drive 的指定文件夹 “GEE_Exports” 中，设置好导出文件名、研究区域、分辨率、最大像素数量和文件格式等参数，并打印出代码执行完成的提示信息。

结果展示：