简介一份基于Google Earth EngineGEE的高分辨率地表温度LST反演源码适合遥感、GIS及相关领域的研究人员与学生。该方案整合Landsat 8热红外数据与Sentinel-2多光谱数据通过随机森林回归与残差融合技术将30米分辨率降尺度至10米有效弥补高分辨率热红外数据稀缺的不足可服务于城市热岛效应分析、精细化生态监测等应用。压缩包共4个文件总大小仅9KB包含核心JavaScript源码.js、HTML辅助页面便于查看说明或结果、开发环境配置.inscode及版本管理文件.gitignore文件结构简洁js主程序可在GEE在线环境直接加载运行。目前已有170人学习。借助这份源码读者可完整掌握数据预处理、遥感指数计算、模型训练、残差校正到精度验证的全流程包括大气校正、云影去除等关键细节同时可针对不同研究区调整参数进行二次开发显著减少GEE编码与LST降尺度方法实现的时间成本。1. GEE高分辨率LST反演30 米地表温度为什么值得自己写源码GEE 高分辨率 LST 反演说白了就是在地球引擎上把 Landsat 热红外波段的辐射数字换算成能直接支撑业务判断的物理温度。我第一次在成都做城市热岛分析时官方 LST 产品够用但想换大气参数、想和 NDVI 同步出图、想批量处理三年影像时官方产品就成了黑匣子于是自己写反演源码这件事就绕不过去了。本文面向遥感、生态、农林和规划方向的工程师目标是让新手能在半天内跑通单景反演让熟手直接拿走避坑清单。下面从算法选型开始逐步落到 GEE 代码、参数调整、验证和批量生产。2. 反演算法选型单窗算法的三个输入与为什么绕开劈窗2.1 单窗算法为什么是 GEE 上 LST 反演的默认解地表温度反演主流有三条路辐射传输方程法RTE、单窗算法、劈窗算法。RTE 需要实时大气剖面在 GEE 里虽然能用 MOD07 或再分析数据凑但逐像元订正流程重单景调试一次要反复等任务队劈窗算法理论上只需要两个热红外通道做强大气校正可 Landsat 8/9 的 TIRS 虽然给了 Band 10 和 Band 11Band 11 的定标不确定性偏大官方都不建议把它当成高精度输入拿它做劈窗等于把噪声喂给模型这个方向天然不踏实。剩下的单窗算法覃志豪等提出只需要一个热红外通道加上地表比辐射率、大气透射率和大气平均作用温度三个输入就能把星上亮温修正到地表温度。这特别适合 GEE 的逐像元运算不需要把大气剖面展开成体数据两三个近似参数就能把误差压在 23 K 以内。公式长这样Ts {a×(1-C-D) [b×(1-C-D)CD]×T - D×Ta} / C其中 C ε×τD (1-τ)×[1(1-ε)×τ]T 是亮温Ta 是大气平均作用温度a、b 是随传感器变化的回归系数。这套结构在 GEE 里拆成 Image 运算很顺手而且跨传感器通用Landsat 5 用 Band 6Landsat 8/9 用 Band 10代码逻辑完全一致只是定标常数不同。所以我个人做 GEE 上的高分辨率 LST默认就是单窗。2.2 数据选型L2 算 NDVI、L1 算亮温别混用 Collection 1确定算法后先解决数据源。GEE 上的 Landsat Collection 2 分两级L1 是原始几何校正产品保留热红外 DN 值和辐射定标系数L2 是表面反射率产品自带 QA_PIXEL 云掩膜波段还直接给了一个官方 LST 波段 SR_B10。很多人直接取 SR_B10 当反演结果这没问题但它不是“自己反演”你拿不到中间过程也没法针对局部传感器换比辐射率。我常用的搭配是 L2 提供云掩膜、NDVI、FVC 和比辐射率计算L1 提供热红外 DN 转辐射亮度和亮温。原因很简单L2 的表面反射率比 TOA 反射率干净NDVI 端元更稳定而 L1 的 TIR 波段是未参与表面反射率反演的原始信号单窗算法要的就是它。两者通过 system:index 配对同一景影像的 L1/L2 产品索引一致空间坐标系也都是 CEP90不需要额外做几何配准。这里特别提醒Collection 2 和 Collection 1 的波段名、QA 位设计、定标系数已经完全不一样。网上老脚本里的“CLOUD_MASK”“RADIANCE_MULT_BAND_103.3420E-04”这类写死参数在 C2 上大概率翻车。我的习惯是任何系数都从影像属性动态读取绝不写死。下面这个表是我每次新建反演工程时的固定对照。用途数据集合波段关键注意NDVI 与 FVCLANDSAT/LC09/C02/T1_L2SR_B4、SR_B5用表面反射率别用 TOA云掩膜同上QA_PIXEL按 bit 位运算不能读 CFMask亮温 BTLANDSAT/LC09/C02/T1B10定标系数从属性读官方 LST 对比LANDSAT/LC09/C02/T1_L2SR_B10验证用单位是 K2.3 大气参数固定值、经验公式与 MOD07 水汽单窗算法的输入里最影响结果的是大气透射率 τ其次是大气平均作用温度 Ta。后者由近地面气温折算前者由大气水汽含量 w 近似得到。GEE 上常见做法分三档快速预览用常数档例如夏季中纬度取 w1.5 g/cm² 代入认真做单景用 MOD07 水汽产品上业务化流程则用再分析数据按影像时间插值。经验公式是这么一组近似τ 0.9743 - 0.0801×wTa 17.326 0.92621×Tair其中 Tair 是近地面气温、单位 K。这套系数在多篇单窗算法扩展论文里都能对上我用下来在 30°N45°N 之间比较稳定到了沿海高湿或青藏高原这种极端区就得换数据源。用 MOD07 取水汽的代码很短// 按影像日期关联 MOD07 水汽产品取研究区均值 var wvImg ee.ImageCollection(MODIS/061/MOD07_L2) .filterBounds(roi) .filterDate(2023-07-01, 2023-07-02) .first(); var wv wvImg.select(Water_Vapor).rename(WV); var wvMean wv.reduceRegion({ reducer: ee.Reducer.mean(), geometry: roi, scale: 5000, maxPixels: 1e9 }).get(WV); var wvVal ee.Number(wvMean).divide(10); // 按产品实际单位换算到 g/cm2这段代码的逻辑是先按研究区和日期取到一景 MOD07 影像再把水汽波段在研究区内平均最后换算成单窗算法需要的单位。需要注意的是不同时期 MOD07 产品的水汽单位有过调整脚本里打印一下 wvMean 的值如果量级在几百到一千说明原始单位是 kg/m² 或 10⁻² cm要按比例折算如果量级在 14 之间才是 g/cm²。单位是这个环节最常见的黑匣子多少人反演温度系统性偏差就是栽在这里。3. 在 GEE 把反演跑通核心代码与参数逐段拆解3.1 筛选影像与云掩膜QA_PIXEL 位运算反演第一步是把影像筛干净。GEE 上的 Landsat Collection 2 云掩膜必须用 QA_PIXEL 波段做位运算不能再翻 Collection 1 时代的 CFMask 老皇历。QA_PIXEL 是一个 16 bit 整型第 2 位是卷云、第 3 位是云、第 4 位是云影、第 5 位是雪想保留干净像元就得把这几类全部排除。var roi ee.FeatureCollection(projects/your_asset/roi); var start 2023-07-01; var end 2023-08-01; // 加载研究区夏季的 Landsat 9 L2 表面反射率产品 var l9sr ee.ImageCollection(LANDSAT/LC09/C02/T1_L2) .filterBounds(roi) .filterDate(start, end) .filter(ee.Filter.lt(CLOUD_COVER, 20)); var maskClouds function(img) { var qa img.select(QA_PIXEL); var cirrusBit 1 2; var cloudBit 1 3; var shadowBit 1 4; var snowBit 1 5; var mask qa.bitwiseAnd(cirrusBit).eq(0) .and(qa.bitwiseAnd(cloudBit).eq(0)) .and(qa.bitwiseAnd(shadowBit).eq(0)) .and(qa.bitwiseAnd(snowBit).eq(0)); return img.updateMask(mask); }; var imgSR l9sr.map(maskClouds).first();这里 CLOUD_COVER 只是一个快速预筛参数它描述整景影像的云量占比并不能逐像元判断所以真正干活的是后面的位运算。bitwiseAnd 的写法是在按位检查每个 flag 是否为 0四个条件同时满足才保留。我建议把 cirrus 位也加进去否则夏季高频出现的薄卷云会在 LST 上制造大量虚假低温斑块。3.2 亮温计算与 L1/L2 配对单窗算法需要的亮温 T只能从 L1 产品的热红外波段来。L2 的 SR_B10 已经是表面温度直接拿去算等于用成品回推原料。同一景影像的 L1 和 L2 在 GEE 里 system:index 相同手动配对很简单不需要 ee.Join代码反而更直白。// 加载同一时间范围的 Landsat 9 L1 TOA 产品 var l9toa ee.ImageCollection(LANDSAT/LC09/C02/T1) .filterBounds(roi) .filterDate(start, end) .filter(ee.Filter.lt(CLOUD_COVER, 20)); // 按 system:index 从 L1 集合里取出对应影像 var toaImg l9toa .filter(ee.Filter.eq(system:index, imgSR.get(system:index))) .first(); // 动态读取定标系数而不是写死常数 var radMult ee.Number(toaImg.get(RADIANCE_MULT_BAND_10)); var radAdd ee.Number(toaImg.get(RADIANCE_ADD_BAND_10)); var k1 ee.Number(toaImg.get(K1_CONSTANT_BAND_10)); var k2 ee.Number(toaImg.get(K2_CONSTANT_BAND_10)); // DN - 辐射亮度 - 亮温公式来自 Landsat 官方定标 var rad toaImg.select(B10).multiply(radMult).add(radAdd); var bt rad.expression( K2 / log(K1 / rad 1), { K1: k1, K2: k2, rad: rad } ).rename(BT);这段代码的核心是“动态读取”四字。Landsat 8 和 Landsat 9 的热红外定标常数不完全相同Landsat 7 更是换了波段号一旦把 L8 的 K1/K2 写死到 L9 影像上整景 LST 会系统性偏移 23 K这是很难排查的玄学偏差。第二个关键点是 rad.expression 的用法rad 本身是一个单波段 ImageK1/K2 是 ee.Numberexpression 会自动把数字提升为常量图参与逐像元运算输出仍然是一张图。3.3 FVC 与地表比辐射率NDVI 阈值混合像元地表比辐射率 ε 是把星上亮温修正到地表温度的关键它不能直接测只能用 NDVI 反演的植被覆盖度 FVC 做线性混合。GEE 上常见的近似公式是 ε 0.004×FVC 0.986其中 FVC 由 NDVI 阈值法给出。这一步我强烈建议用 L2 表面反射率算 NDVI因为 TOA NDVI 在山地和城市阴影区会偏低进而把 FVC 压小最后 ε 整体被拉低。// NDVI 和 NDWINDWI 用于识别水面 var ndvi imgSR.normalizedDifference([SR_B5, SR_B4]).rename(NDVI); var ndwi imgSR.normalizedDifference([SR_B3, SR_B5]).rename(NDWI); // 植被覆盖度NDVI 0.2 到 0.5 之间做非线性插值 var fvc ndvi.expression( ((NDVI - 0.2) / (0.5 - 0.2)) ** 2, {} ).rename(FVC); fvc fvc.where(ndvi.lt(0.2), 0).where(ndvi.gt(0.5), 1); // 混合像元比辐射率水体单独给 0.991 var em fvc.expression(0.004 * FVC 0.986).rename(EM); var em em.where(ndwi.gt(0), 0.991);FVC 端点的选择有讲究。NDVI 低于 0.2 通常视为裸土或建筑区高于 0.5 视为完全植被覆盖这是阈值法的通用做法但到了干旱区或城市核心区裸土端还要再细分。水体识别用 NDWI 而不是 NDVI是因为水体在近红外波段吸收很强NDWI 为正的像元基本可以判定为水面比 NDVI 阈值稳得多。这里把 GEE 上常见的 NDWI 用法嵌进了比辐射率流程城市热岛研究里这一步做不好湖面温度能反演出比周围还低的“凉岛”看着合理其实是错的。3.4 单窗公式落成 Image 运算并导出三样输入齐了亮温 BT、比辐射率 EM、大气参数 τ 和 Ta。下面这段把单窗公式完整落到 GEE 的 Image 运算里。为了便于复现我先用常数档大气参数替换成 MOD07 的做法在第 2.3 节已经给了。// 常数档大气参数夏季中纬度 w1.5Tair300 var tau ee.Image.constant(0.9743 - 0.0801 * 1.5); // 约 0.854 var ta ee.Image.constant(17.326 0.92621 * 300); // 约 295 K var a -62.8065; var b 0.4338; var c em.multiply(tau); var d ee.Image.constant(1).subtract(tau) .multiply(ee.Image.constant(1) .add(ee.Image.constant(1).subtract(em).multiply(tau))); var part1 ee.Image.constant(a).multiply(ee.Image.constant(1).subtract(c).subtract(d)); var part2 ee.Image.constant(b).multiply(ee.Image.constant(1).subtract(c).subtract(d)) .add(c).add(d); var lst part1.add(part2.multiply(bt)).subtract(d.multiply(ta)).divide(c) .clip(roi).rename(LST); // 导出为 30 米分辨率 GeoTIFF Export.image.toDrive({ image: lst.toFloat(), description: LST_ imgSR.get(system:index), folder: GEE_LST, scale: 30, crs: EPSG:4326, maxPixels: 1e13 });参数说明在这里要一次讲透。a、b 是单窗算法针对 Landsat 8/9 TIRS 的回归系数参考覃志豪单窗算法扩展的常见取值不同论文里略有差异但量级一致C 是比辐射率与透射率的乘积D 是大气下行辐射的等效项这两项决定了大气在整条辐射路径上的损耗。导出时 scale 设为 30 并不是把 100 米热红外原生分辨率变成真的 30 米而是以 Landsat 多光谱的 30 米网格重采样输出这也是标题里“高分辨率”的工程含义。若用 panchromatic 锐化去硬凑 15 米 LST只会引入大量伪纹理不建议。4. 反演结果验证与避坑一条温度曲线是怎么偏出 5℃ 的4.1 三级对照官方 ST 产品、MODIS、气象站自己写的反演代码出结果快但没人能保证第一次就对。我的验证习惯是三级对照第一级用官方 L2 的 SR_B10 产品做逐像元差值和 RMSE 统计第二级用 MODIS MYD11A1 的 1 公里 LST 做尺度对比第三级用气象站或实测点位做地面验证。第一级代码最短也最直接// 官方 LST 与自反演 LST 做差值统计 var diff lst.subtract(imgSR.select(SR_B10)).rename(DIFF); var stats diff.reduceRegion({ reducer: ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true), geometry: roi, scale: 30, maxPixels: 1e10 }); print(mean diff:, stats.get(DIFF_mean)); print(std diff:, stats.get(DIFF_stdDev));差值均值如果落在 ±1 K 以内说明大气参数和比辐射率基本合理如果达到 2 K 以上优先查水汽常数如果空间分布上呈板块状偏移大概率是云掩膜还有漏网。第二级对照要注意过境时刻Landsat 过境约上午 10:30MODIS Terra 约上午 10:30 相近Aqua 则接近下午 13:30选 Terra 结果做比较更合理。第三级不推荐直接拿百叶箱气温对比因为气温和地表温度本来就不是一个物理量最好找辐射温度计或通量站数据。4.2 避坑清单坑一Collection 1 老脚本直接套 C2。现象是反演温度整体偏移 24 K且影像边缘出现条带。原因是 C2 改了 QA 波段设计、辐射定标系数也从属性里换了一套写死的参数全部失效。解决方法是统一走属性读取并彻底放弃旧的 CFMask 写法。坑二薄卷云没掩掉出现孤立低温像元。现象是 LST 图上出现大量 -10℃ 左右的异常斑块恰好和影像上的薄云位置重合。原因是 QA_PIXEL 只掩了 bit3 云没掩 bit2 卷云。解决办法是把 bit2 也加进位运算里代价是会丢少量干净像元但对温度产品来说值得。坑三沿海高湿区反演偏低。现象是同一套代码在成都正常放到广州就整体低 34 K。原因是固定 w1.5 g/cm² 在沿海潮湿大气下严重低估水汽含量导致 τ 偏大。解决方法是改用第 2.3 节的 MOD07 水汽并按天匹配影像过境时间。坑四城市核心区白天 LST 比官方产品高 8℃。现象是在沥青、屋顶这类像元上温度异常高。原因是 NDVI 接近 0 的像元被当成了裸土采用了植被端混合公式而实际不透水面的比辐射率只有 0.920.96 左右。解决方法是引入不透水面比例数据集把混合公式的裸土端换成 0.958 左右。坑五官方 ST 和自反演 LST 完全一致反而要警惕。现象是两者偏差正好为零看起来完美。原因是官方 ST 用了更完善的 RTE 和实时大气剖面近似算法在个别像元上不可能分毫不差如果逐像元完全吻合说明代码里很可能错把官方 ST 波段当成了输入。正确的做法是把官方产品当成“可用锚点”允许存在均值为零附近、空间随机分布的 1 K 级差异。5. 进阶批量反演与月合成把源码变成业务产品单景跑通只是开始真正让 LST 源码产生业务价值的是批量化和时间序列化。我会把所有计算步骤封装成一个纯函数输入 L2 和 L1 两景配对影像输出 LST 单波段然后对整个影像集合做 map。这一步可以顺手把 function 保存成 GEE Module换研究区时只改 roi 和日期范围。封装后的批量月合成大概长这样// 把单景反演封装成函数输入 SR 影像和 TOA 影像输出 LST var calcLST function(srImg, toaImg) { // 依次调用云掩膜、NDVI/FVC/epsilon、亮温、单窗公式 return lst; }; // 对集合逐景反演并保留时间标记 var lstCol l9sr.map(function(sr) { var toa l9toa .filter(ee.Filter.eq(system:index, sr.get(system:index))) .first(); return calcLST(sr, toa) .set(system:time_start, sr.get(system:time_start)) .set(system:index, sr.get(system:index)); });批量跑之前务必先做一件事拿一个月的影像和官方 ST 产品逐景计算差值确认每景的 mean diff 都落在合理范围。我最早的版本把水面像元错判成裸土湖面温度比实测高了 7℃就是因为只看了单景漂亮结果就上了批量后来养成了“先看差值图和任务队列再谈业务合成”的习惯。批量任务在 GEE 的 Admin 后台会排长队导出前把 scale、crs、maxPixels 写死避免几十个任务因为参数错误集体失败。月合成时不要直接 mean 所有时相白天和夜间 LST 分开处理并优先用同一过境时段的数据否则月均值会被过境时刻差异污染。这套流程跑通后你手里的就不只是几景温度图而是一套可回归、可校准、可追溯的地表温度时间序列产品。希望帮到你。本文还有配套的精品资源点击获取