GEE高分辨率LST反演实战：从源码到可复现的完整链路

简介这份源码资源面向遥感、地理信息与生态环境方向的学习者和研究者聚焦于在Google Earth Engine平台上实现高分辨率地表温度反演。它解决的是单一卫星热红外数据空间分辨率不足的问题通过融合Landsat 8与Sentinel-2多源数据结合随机森林回归建模与残差融合技术将30米LST降尺度至10米可服务于城市热岛分析与生态环境监测等场景。资源包共4个文件以js脚本为核心辅以html页面、inscode配置与gitignore文件整体约9KB结构轻量便于直接导入GEE运行与二次修改。目前已有170人学习。读者可从中获取从数据加载、预处理、指数计算、模型训练到残差处理与结果验证的完整代码流程理解多源数据融合的建模思路与精度优化方法并据此迁移到自己的研究区域或课题中。1. GEE高分辨率LST反演从源码到可复现的完整链路城市热岛、矿区热异常、农业旱情监测这些场景都绕不开一个核心指标——地表温度LST。传统做法是下载 Landsat 影像到本地用 ENVI 或 ArcGIS 跑单窗算法一张影像折腾半小时一个研究区几十景影像就是几天的工作量。GEEGoogle Earth Engine把这个流程搬到了云端LST 反演从「本地跑单张」变成「云端批量出图」源码写一次后面换研究区只改几何边界就行。但问题也在这里网上流传的 GEE LST 反演源码质量参差不齐有的用错热红外波段有的把亮温当 LST 直接输出有的没做大气校正就敢标「高分辨率」。这篇笔记不讲概念科普直接拆一套能跑通的 GEE 高分辨率 LST 反演源码从波段选择、亮温计算、比辐射率估算到最终出图每一步的参数含义和踩坑点都写清楚。适合已经会用 GEE 基础 API、想把手上的热红外影像真正变成可用 LST 产品的从业者。2. 高分辨率LST反演的技术底座为什么选Landsat 8/9 TIRS2.1 热红外波段的选择逻辑与分辨率权衡Landsat 8/9 搭载的 TIRSThermal Infrared Sensor有两个热红外波段Band 1010.60–11.19 μm和 Band 1111.50–12.51 μm。很多人第一反应是用 Band 11因为它的波长更长理论上对低温目标更敏感。但实际操作中Band 11 从 2013 年起就存在杂散光stray light问题USGS 官方建议在大多数定量反演中优先使用 Band 10。我一般直接锁定 Band 10除非研究区有特殊需求。空间分辨率方面TIRS 原始分辨率是 100 米但 GEE 上的 Landsat 8/9 Collection 2 产品已经把它重采样到 30 米和 OLI 的多光谱波段对齐。这意味着你拿到的 LST 产品名义上是 30 米分辨率但热红外本身的物理分辨率还是 100 米——重采样不会凭空增加信息量。这一点在写论文或做报告时要诚实说明别把 30 米重采样当成真正的 30 米热红外。提示如果研究区需要更高分辨率的热红外数据可以考虑 ECOSTRESS约 70 米或 Sentinel-3 SLSTR约 1 公里但前者时间分辨率低后者空间分辨率更差。Landsat 8/9 在「免费 30 米重采样 16 天重访」这个组合下仍然是大多数陆地应用的最优解。2.2 亮温计算从DN值到辐射亮温的完整公式GEE 上的 Landsat Collection 2 Level-1 产品已经提供了 TOA大气顶层辐射亮度但要做 LST 反演第一步是把辐射亮度转成亮温Brightness Temperature。公式如下// Landsat 8/9 Band 10 亮温计算 // ST_B10 是 Collection 2 Level-2 的地表温度产品但我们要自己反演 // 先用 Level-1 的 TOA 辐射亮度 var radiance image.select(B10).multiply(0.0003342).add(0.1); // 亮温计算常数Landsat 8/9 TIRS Band 10 var K1 774.8853; // 单位W/(m²·sr·μm) var K2 1321.0789; // 单位K // 亮温单位开尔文 var bt radiance.divide(K1).add(1).log().multiply(K2).rename(BT);这段代码的逻辑是先把 DN 值乘以辐射定标系数0.0003342再加上偏移量0.1得到 TOA 辐射亮度然后用 Planck 函数的反函数计算亮温。K1 和 K2 是 Landsat 8/9 TIRS Band 10 的定标常数直接从元数据里取。注意这里的log()是自然对数GEE 的ee.Image.log()默认就是自然对数不需要换底。参数说明0.0003342和0.1是 Collection 2 Level-1 的辐射定标系数不同传感器版本可能不同建议从影像元数据里动态读取。K1 和 K2 的值在 USGS 官方文档里有明确给出不要凭记忆写。2.3 比辐射率估算NDVI阈值法的三个分支亮温只是「传感器看到的温度」要得到真正的 LST还需要除以比辐射率emissivity。目前最常用的方法是 NDVI 阈值法Sobrino 等人 2008 年提出的方案根据 NDVI 值分三种情况// NDVI 计算 var ndvi image.normalizedDifference([B5, B4]).rename(NDVI); // 比辐射率估算NDVI 阈值法 var em ee.Image(0).expression( (ndvi 0.157) ? 0.923 : (ndvi 0.727) ? 0.994 : (0.923 (ndvi - 0.157) * (0.994 - 0.923) / (0.727 - 0.157)), {ndvi: ndvi} ).rename(EM);这段表达式的三个分支分别对应NDVI ≤ 0.157 时地表视为裸土比辐射率取 0.923NDVI ≥ 0.727 时地表视为全植被覆盖比辐射率取 0.994中间区域用线性插值。这两个阈值来自 Sobrino 的论文适用于大多数陆地场景。注意水体区域的比辐射率通常取 0.985 左右但 NDVI 阈值法对水体的处理不够准确。如果研究区有大面积水体建议单独用 NDWI 掩膜水体再赋一个固定比辐射率值。2.4 单窗算法从亮温和比辐射率到LST有了亮温和比辐射率最后一步是单窗算法Single-Channel Algorithm。最简化的形式是// 单窗算法反演 LST var lst bt.divide( ee.Image(1).add( (bt.multiply(0.00115)).divide(em.multiply(14380)) ) ).rename(LST); // 转成摄氏度 var lst_c lst.subtract(273.15).rename(LST_C);这里的0.00115是波长相关的常数对于 Band 10 的中心波长 10.9 μm14380是 Planck 常数与玻尔兹曼常数的比值。这个公式是简化版没有考虑大气水汽含量。如果研究区大气水汽变化大需要引入大气透过率参数那就得用 MODIS 水汽产品或者 NCEP 再分析数据做大气校正。我一般会在代码里加一个开关让用户选择是否做大气校正。对于大多数城市热岛研究简化版已经够用如果要做精度验证建议至少用 MODIS 的 LST 产品做交叉对比。3. 在GEE上跑通LST反演从影像筛选到批量导出的完整代码3.1 研究区定义与影像筛选日期、云量、边界第一步是定义研究区和时间范围。GEE 的ee.Geometry支持手动输入坐标、导入 Shapefile 或直接用ee.Geometry.Rectangle。我一般用ee.Geometry.Rectangle快速测试正式跑的时候换成导入的矢量边界。// 定义研究区以北京市为例 var roi ee.Geometry.Rectangle([116.0, 39.6, 116.8, 40.2]); // 时间范围 var startDate 2023-06-01; var endDate 2023-08-31; // 筛选 Landsat 8/9 Collection 2 Level-1 影像 var collection ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .merge(ee.ImageCollection(LANDSAT/LC09/C02/T1_L2)) .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt(CLOUD_COVER, 20)) .sort(CLOUD_COVER);这里有几个关键点merge把 Landsat 8 和 9 的影像合并增加可用影像数量filterBounds筛选覆盖研究区的影像filterDate限定时间范围filter用云量过滤我一般设 20% 以下sort按云量升序排列方便后续取最少云的影像。提示Collection 2 Level-2 产品已经提供了ST_B10波段地表温度但那是 USGS 官方反演的结果精度和参数不可控。如果你想自己控制反演过程用 Level-1 产品从头算如果只是要一个快速结果直接用 Level-2 的ST_B10也行。3.2 逐影像反演封装一个可复用的LST计算函数把前面几节的步骤封装成一个函数输入是单景影像输出是带 LST 波段的影像。这样后面用map批量处理时代码更干净。// 封装 LST 反演函数 function calculateLST(image) { // 辐射定标 var radiance image.select(B10).multiply(0.0003342).add(0.1); // 亮温计算 var K1 774.8853; var K2 1321.0789; var bt radiance.divide(K1).add(1).log().multiply(K2).rename(BT); // NDVI 计算 var ndvi image.normalizedDifference([B5, B4]).rename(NDVI); // 比辐射率估算 var em ee.Image(0).expression( (ndvi 0.157) ? 0.923 : (ndvi 0.727) ? 0.994 : (0.923 (ndvi - 0.157) * (0.994 - 0.923) / (0.727 - 0.157)), {ndvi: ndvi} ).rename(EM); // 单窗算法 var lst bt.divide( ee.Image(1).add( (bt.multiply(0.00115)).divide(em.multiply(14380)) ) ).rename(LST); // 转摄氏度 var lst_c lst.subtract(273.15).rename(LST_C); return image.addBands([bt, ndvi, em, lst, lst_c]); } // 批量处理 var lstCollection collection.map(calculateLST);这个函数的逻辑很清晰先算亮温再算 NDVI 和比辐射率最后用单窗算法得到 LST。返回的影像里包含了中间变量BT、NDVI、EM方便后续做精度分析。参数说明0.0003342和0.1是辐射定标系数774.8853和1321.0789是亮温常数0.00115和14380是单窗算法常数。这些值都是针对 Landsat 8/9 TIRS Band 10 的换成其他传感器需要重新查。3.3 时间合成与导出中值合成、裁剪、Export.image批量反演完下一步是做时间合成。我一般用中值合成median因为 LST 受云影响大中值比均值更稳健。// 时间合成 var lstMedian lstCollection.select(LST_C).median().clip(roi); // 可视化参数 var visParams { min: 20, max: 45, palette: [blue, green, yellow, orange, red] }; // 添加到地图 Map.centerObject(roi, 10); Map.addLayer(lstMedian, visParams, LST Median); // 导出到 Drive Export.image.toDrive({ image: lstMedian, description: LST_Beijing_2023_Summer, folder: GEE_LST, scale: 30, region: roi, maxPixels: 1e13, crs: EPSG:4326 });导出时注意几个参数scale设 30 米和 Landsat 重采样后的分辨率一致region是研究区边界maxPixels设大一点避免导出失败crs我一般用 EPSG:4326方便后续在 QGIS 或 ArcGIS 里打开。注意如果研究区很大导出时间可能很长。建议先用小区域测试确认代码没问题再跑全区域。另外GEE 的导出任务在 Tasks 面板里手动启动不会自动跑。3.4 精度验证用MODIS LST做交叉对比自己反演的 LST 到底准不准最直接的方法是和 MODIS 的 LST 产品做交叉对比。MODIS 的 MOD11A1 产品提供 1 公里分辨率的 LST虽然空间分辨率低但精度经过验证。// 加载 MODIS LST 产品 var modisLST ee.ImageCollection(MODIS/061/MOD11A1) .filterBounds(roi) .filterDate(startDate, endDate) .select(LST_Day_1km) .median() .multiply(0.02) .subtract(273.15) .clip(roi); // 重采样到 30 米方便对比 var modisLST_30m modisLST.resample(bilinear).reproject({ crs: EPSG:4326, scale: 30 }); // 计算差值 var diff lstMedian.subtract(modisLST_30m).rename(LST_Diff); // 统计差值 var stats diff.reduceRegion({ reducer: ee.Reducer.mean().combine({ reducer2: ee.Reducer.stdDev(), sharedInputs: true }), geometry: roi, scale: 30, maxPixels: 1e13 }); print(LST Difference Stats:, stats);这段代码的逻辑是加载 MODIS LST转成摄氏度重采样到 30 米然后和自研 LST 做差值统计。如果平均差值在 ±2°C 以内说明反演结果基本可靠如果差值超过 5°C需要检查比辐射率估算或大气校正环节。参数说明0.02是 MODIS LST 的缩放因子273.15是开尔文到摄氏度的转换。resample(bilinear)用双线性插值重采样reproject指定目标分辨率和坐标系。4. GEE LST反演避坑指南5个血泪教训4.1 坑一Band 11 杂散光导致LST异常偏高现象用 Band 11 反演的 LST 比 Band 10 高出 5–10°C尤其在冬季或高纬度地区更明显。原因Landsat 8 TIRS Band 11 从 2013 年起存在杂散光问题USGS 在 2017 年才公开承认。杂散光来自 Band 10 的辐射泄漏导致 Band 11 的亮温被高估。解决直接用 Band 10别碰 Band 11。如果非要用双波段算法至少做杂散光校正但 GEE 上没有现成的校正工具得自己写。4.2 坑二NDVI阈值法在水体区域失效现象水体区域的 LST 反演结果异常要么偏高要么偏低和实测数据对不上。原因NDVI 阈值法的两个阈值0.157 和 0.727是针对陆地表面的水体 NDVI 通常为负值落在第一个分支里被当成裸土处理比辐射率取 0.923但实际水体的比辐射率在 0.985 左右。解决用 NDWI 掩膜水体单独赋比辐射率值。代码里加一个where判断var ndwi image.normalizedDifference([B3, B5]).rename(NDWI); var em em.where(ndwi.gt(0), 0.985);4.3 坑三云掩膜不彻底导致LST被云污染现象LST 图上出现大片异常低温区域形状和云层一致。原因Landsat 的CLOUD_COVER属性是整景影像的云量百分比不代表研究区内的云量。如果研究区刚好被云覆盖但整景云量低于 20%这景影像还是会被选进来。解决用QA_PIXEL波段做逐像元云掩膜而不是只靠CLOUD_COVER属性。代码var qa image.select(QA_PIXEL); var cloudMask qa.bitwiseAnd(1 3).eq(0); var imageMasked image.updateMask(cloudMask);4.4 坑四导出时maxPixels不够导致任务失败现象Export 任务在 Tasks 面板里显示「Failed」错误信息提到「maxPixels」。原因GEE 默认的maxPixels是 1e8对于大区域或高分辨率导出不够用。解决在Export.image.toDrive里显式设置maxPixels: 1e13。如果还失败检查scale是不是设得太小或者研究区是不是太大。4.5 坑五时间合成用均值导致LST被异常值拉偏现象中值合成的 LST 图看起来正常但均值合成的图上有明显的异常高值或低值区域。原因LST 受云、阴影、气溶胶影响大均值对异常值敏感。中值更稳健但也不是万能的——如果研究区在时间范围内大部分被云覆盖中值也可能不可靠。解决先用qualityMosaic或者自定义的云掩膜做预处理再用中值合成。如果影像数量太少少于 5 景考虑扩大时间范围或换传感器。5. 进阶技巧用GEE Python API批量跑LST并做时序分析5.1 从JavaScript到PythonAPI迁移的注意事项GEE 的 JavaScript API 适合快速原型但要做批量处理或集成到本地工作流Python API 更方便。迁移时注意几个差异JavaScript 的var在 Python 里直接赋值ee.Image.expression的语法基本一致Map.addLayer在 Python 里用geemap库替代。import ee import geemap ee.Initialize() roi ee.Geometry.Rectangle([116.0, 39.6, 116.8, 40.2]) def calculate_lst(image): radiance image.select(B10).multiply(0.0003342).add(0.1) bt radiance.divide(774.8853).add(1).log().multiply(1321.0789).rename(BT) ndvi image.normalizedDifference([B5, B4]).rename(NDVI) em ee.Image(0).expression( (ndvi 0.157) ? 0.923 : (ndvi 0.727) ? 0.994 : (0.923 (ndvi - 0.157) * (0.994 - 0.923) / (0.727 - 0.157)), {ndvi: ndvi} ).rename(EM) lst bt.divide( ee.Image(1).add( (bt.multiply(0.00115)).divide(em.multiply(14380)) ) ).rename(LST) return image.addBands(lst.subtract(273.15).rename(LST_C)) collection ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) \ .merge(ee.ImageCollection(LANDSAT/LC09/C02/T1_L2)) \ .filterBounds(roi) \ .filterDate(2023-06-01, 2023-08-31) \ .filter(ee.Filter.lt(CLOUD_COVER, 20)) lst_collection collection.map(calculate_lst) lst_median lst_collection.select(LST_C).median().clip(roi) # 导出 geemap.ee_export_image( lst_median, filenameLST_Beijing_2023.tif, scale30, regionroi, file_per_bandFalse )Python 版的逻辑和 JavaScript 完全一致只是语法不同。geemap.ee_export_image直接下载到本地不用经过 Drive。5.2 时序LST分析按月合成与变化趋势如果要分析 LST 的季节变化可以按月做合成然后计算趋势。# 按月合成 months ee.List.sequence(1, 12) def monthly_lst(month): start ee.Date.fromYMD(2023, month, 1) end start.advance(1, month) return lst_collection.filterDate(start, end).select(LST_C).median() \ .set(month, month) monthly_collection ee.ImageCollection(months.map(monthly_lst)) # 计算趋势简单线性回归 def add_trend(image): month image.get(month) return image.addBands(ee.Image(month).rename(month)) trend_collection monthly_collection.map(add_trend) trend trend_collection.select([month, LST_C]).reduce(ee.Reducer.linearFit()) # 导出趋势图 geemap.ee_export_image( trend.select(scale), filenameLST_Trend_2023.tif, scale30, regionroi )这段代码的逻辑是按月合成 LST然后对每个像元做线性回归scale波段就是趋势斜率单位°C/月。正值表示升温趋势负值表示降温趋势。5.3 一个容易忽略的细节投影与重采样GEE 默认用 EPSG:4326 投影但 Landsat 原始投影是 UTM。如果直接导出GEE 会自动重投影可能引入误差。我一般会在导出时指定crs: EPSG:32650北京地区的 UTM Zone 50N保持和原始影像一致。另外重采样方法也很关键。对于 LST 这种连续变量双线性插值bilinear比最近邻nearest更合适。在 GEE 里可以用resample(bilinear)显式指定。var lst_resampled lstMedian.resample(bilinear).reproject({ crs: EPSG:32650, scale: 30 });这个细节看起来小但在做多时相对比或和实测数据验证时投影不一致会导致明显的空间偏移。我在这上面翻过车后来养成了导出前先检查投影的习惯。希望帮到你。本文还有配套的精品资源点击获取