GEE高分辨率LST反演实战：从亮温计算到劈窗算法的完整代码与避坑指南

简介：这份源码资源面向遥感、地理信息与生态环境方向的开发者及研究人员，聚焦于在Google Earth Engine平台上实现高分辨率地表温度反演。它针对单一卫星热红外数据分辨率不足的痛点，融合Landsat 8与Sentinel-2多源数据，通过随机森林回归建模与残差融合技术，将30米LST降尺度至10米，可用于城市热岛分析与生态环境监测等场景。资源包共4个文件，包含1个js核心脚本、1个inscode配置、1个html页面及1个gitignore，整体约9KB，结构精简，便于直接部署到GEE环境运行。目前已有170人学习下载。读者可获取从数据加载、预处理、指数计算、模型训练到残差处理与结果验证的完整代码链路，理解多源数据融合的建模思路与实现细节，为自身的热环境研究或遥感项目提供可复用的技术参考。

1. 从「GEE高分辨率LST反演」说起：为什么你拿到的地表温度总像蒙了层雾

做过城市热岛、蒸散发或者农业旱情监测的人，大概率都碰过同一个场景：在 GEE 里调出 Landsat 8/9 的 Collection 2 Level 2 产品，用ST_B10波段乘个系数、减个偏移，出图一看——整片城区像被糊了一层灰，水体温度偏高，裸地又偏低，跟地面实测差了三四度。这不是你代码写错了，而是「高分辨率 LST 反演」这件事本身在 GEE 上有一堆默认行为在跟你作对。

GEE 高分辨率 LST 反演，说白了就是在 Google Earth Engine 这个云端平台上，用 Landsat（30 米）、Sentinel-3 SLSTR（1 公里）甚至 MODIS 做辅助，把热红外波段的辐亮度一步步换算成地表温度。它解决的核心问题是：不用下载几十 GB 影像到本地、不用装 ENVI 或 ArcGIS，直接在浏览器里跑完辐射定标、大气校正、比辐射率估算、劈窗算法这一整条链路。适合谁？做城市热环境、矿区热异常、农田蒸散、火山监测的从业者，尤其是手头没有高性能工作站、又需要批量处理多年时序的人。

但「源码」两个字是关键词。GEE 上现成的 LST 产品（比如 MODIS MOD11A1）分辨率太粗，Landsat 官方又没有直接给 LST 产品，所以你必须自己写反演代码。而网上流传的很多脚本，要么是单窗算法参数写死、要么比辐射率用 NDVI 阈值法一刀切，换个研究区就翻车。这篇笔记就按我自己踩过的路，把这条链路从原理到代码、从参数到避坑，完整讲一遍。

2. 反演链路拆解：从热红外辐亮度到地表温度的四个环节

2.1 为什么不能直接用 ST_B10 乘系数了事

Landsat Collection 2 Level 2 的ST_B10波段，官方文档里写的是「Surface Temperature」，单位开尔文，缩放因子 0.00341802，偏移 149.0。很多人一看，直接ST_B10 * 0.00341802 + 149.0就出图了。这个做法在 2020 年之后的产品里确实能出一个「看起来像温度」的东西，但它有两个致命问题。

第一，这个ST_B10是 USGS 用他们自己的单窗算法反演出来的，算法内部用的大气参数是全球平均的，没有针对你的研究区、你的成像时间做校正。第二，它没有提供比辐射率波段，而比辐射率对 LST 的影响在 NDVI 低的区域（城市、裸地）能到 2-3 K。所以如果你要做定量分析，尤其是城市热岛这种地表类型复杂的场景，直接用ST_B10是不够的。

常见做法是：用 Level 1 的 TOA 辐亮度，自己走一遍辐射传输方程或者劈窗算法。GEE 上 Landsat 8/9 的 Level 1 数据是LANDSAT/LC08/C02/T1_TOA和LANDSAT/LC09/C02/T1_TOA，热红外波段是B10（10.6-11.19 μm）和B11（11.5-12.51 μm）。这两个波段就是劈窗算法的输入。

2.2 辐射定标与亮温计算：第一步就有人栽跟头

从 TOA 辐亮度到亮温，用的是普朗克函数的反函数。GEE 里B10波段的 TOA 辐亮度已经由 USGS 做了辐射定标，单位是 W/(m²·sr·μm)。亮温公式是：

T_b = K2 / ln(K1 / L + 1)

其中 L 是 TOA 辐亮度，K1 和 K2 是定标常数。Landsat 8 的 B10：K1 = 774.8853 W/(m²·sr·μm)，K2 = 1321.0789 K；B11：K1 = 480.8883，K2 = 1201.1442。Landsat 9 的常数略有不同，B10 的 K1 = 799.0284，K2 = 1329.2405，B11 的 K1 = 475.6581，K2 = 1198.817。这些常数在 GEE 的LANDSAT/LC08/C02/T1_TOA元数据里其实有，但很多人直接抄网上的旧值，Landsat 9 用了 Landsat 8 的常数，结果亮温差 0.5 K 左右。

代码上，我一般这样写：

// 计算亮温，以 Landsat 8 为例 var K1_B10 = 774.8853; var K2_B10 = 1321.0789; var K1_B11 = 480.8883; var K2_B11 = 1201.1442; var addBrightnessTemp = function(image) { var L_B10 = image.select('B10'); // TOA 辐亮度 var L_B11 = image.select('B11'); var BT_B10 = L_B10.multiply(0).add(K2_B10) .divide(L_B10.multiply(0).add(K1_B10).divide(L_B10).add(1).log()) .rename('BT_B10'); // 上面这种写法是为了规避 GEE 里除零和 log 的边界问题 var BT_B11 = L_B11.multiply(0).add(K2_B11) .divide(L_B11.multiply(0).add(K1_B11).divide(L_B11).add(1).log()) .rename('BT_B11'); return image.addBands([BT_B10, BT_B11]); };

这里multiply(0).add()的写法看起来绕，但它是为了处理 GEE 的惰性计算和无效值传播。直接写K2.divide(K1.divide(L).add(1).log())在 L 为 0 或负值时会产生 NaN，而 GEE 的log()对 0 返回 -Infinity，后续整幅影像就废了。用multiply(0).add()相当于给每个像素加了一个常数基底，避免除零。

2.3 比辐射率估算：NDVI 阈值法的三个版本和适用边界

比辐射率（ε）是 LST 反演里最容易被低估的变量。常见做法有三种：

第一种是 NDVI 阈值法（Sobrino 2008），把地表分成水体、裸地、植被三类。NDVI < 0.2 算裸地，ε = 0.97；NDVI > 0.5 算全植被，ε = 0.99；中间用公式 ε = 0.004 * Pv + 0.986，其中 Pv 是植被覆盖度，Pv = ((NDVI - NDVI_min) / (NDVI_max - NDVI_min))²。这个方法在 GEE 里实现简单，但 NDVI_min 和 NDVI_max 取全局还是局部，结果差很多。

第二种是 NDVI 阈值法的改进版，用 NDVI 和 NDVI 的方差来动态确定阈值。这个在 GEE 里需要先算研究区的 NDVI 统计，代码量翻倍，但城市区域效果明显好。

第三种是 TES（Temperature Emissivity Separation）算法，用多波段热红外数据同时反演温度和比辐射率。Landsat 8/9 只有两个热红外波段，TES 效果有限，一般还是用 NDVI 法。

我一般用第一种，但会把 NDVI_min 和 NDVI_max 从全局改成研究区内的 5% 和 95% 分位数。代码：

// 计算 NDVI 和植被覆盖度 var addNDVI = function(image) { var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI'); return image.addBands(ndvi); }; // 用研究区统计确定 NDVI_min 和 NDVI_max var ndviStats = ndviImage.reduceRegion({ reducer: ee.Reducer.percentile([5, 95]), geometry: roi, scale: 30, maxPixels: 1e13 }); var ndviMin = ee.Number(ndviStats.get('NDVI_p5')); var ndviMax = ee.Number(ndviStats.get('NDVI_p95')); // 植被覆盖度 var Pv = ndviImage.subtract(ndviMin).divide(ndviMax.subtract(ndviMin)).pow(2).clamp(0, 1); // 比辐射率 var emissivity = Pv.multiply(0.004).add(0.986).rename('EM');

这里clamp(0, 1)是必须的，因为 NDVI 可能超出 [ndviMin, ndviMax] 范围，Pv 会变成负数或大于 1，比辐射率就跑到 0.98 以下或 0.99 以上，温度反演直接偏 1-2 K。

2.4 劈窗算法：两个波段怎么组合出温度

劈窗算法的核心思想是：用两个相邻热红外波段对大气水汽的吸收差异，来校正大气影响。Landsat 8/9 的 B10 和 B11 正好满足这个条件。常用的劈窗公式是：

LST = BT_B10 + a * (BT_B10 - BT_B11) + b * (BT_B10 - BT_B11)² + c + d * (1 - ε) + e * Δε

其中 a, b, c, d, e 是系数，不同文献给的系数不一样。我一般用 Jiménez-Muñoz 2014 的系数：a = 1.02, b = 0.14, c = -0.15, d = 75, e = -120。这些系数是在全球范围内拟合的，对大部分区域够用，但如果你做的是极端干旱区或者高海拔区，最好用当地探空数据重新拟合。

代码实现：

// 劈窗算法 var a = 1.02, b = 0.14, c = -0.15, d = 75, e = -120; var deltaE = emissivity_B10.subtract(emissivity_B11); // 两个波段的比辐射率差 var meanE = emissivity_B10.add(emissivity_B11).divide(2); var LST = BT_B10.add(BT_B10.subtract(BT_B11).multiply(a)) .add(BT_B10.subtract(BT_B11).pow(2).multiply(b)) .add(c) .add(meanE.multiply(0).add(1).subtract(meanE).multiply(d)) .add(deltaE.multiply(e)) .rename('LST');

注意meanE.multiply(0).add(1).subtract(meanE)就是1 - meanE，写成这样是为了避免 GEE 里1 - image的语法问题。最后 LST 单位是开尔文，减 273.15 转摄氏度。

3. 在 GEE 里跑通一景 Landsat 9 的完整代码与参数表

3.1 数据筛选与云掩膜：别让云把温度带偏

Landsat 9 的 TOA 数据是LANDSAT/LC09/C02/T1_TOA，云掩膜用 QA_PIXEL 波段。QA_PIXEL 的 bit 3 是云，bit 4 是云阴影，bit 5 是雪。我一般把云和云阴影都掩掉，雪保留但标记出来。

var maskClouds = function(image) { var qa = image.select('QA_PIXEL'); var cloud = qa.bitwiseAnd(1 << 3).eq(0); var shadow = qa.bitwiseAnd(1 << 4).eq(0); var mask = cloud.and(shadow); return image.updateMask(mask); }; var collection = ee.ImageCollection('LANDSAT/LC09/C02/T1_TOA') .filterDate('2023-06-01', '2023-08-31') .filterBounds(roi) .filter(ee.Filter.lt('CLOUD_COVER', 20)) .map(maskClouds);

1 << 3是位运算，等于 8。QA_PIXEL 里云是 bit 3，所以bitwiseAnd(8)提取出来，eq(0)表示没有云。云阴影是 bit 4，1 << 4等于 16。两个条件都满足才保留。

3.2 逐景反演与合成：中值合成还是均值合成

单景 LST 反演完之后，如果研究区跨多景，需要合成。常见做法是取时间序列的中值。中值对异常值（比如残余薄云）更稳健，均值容易被带偏。

var lstCollection = collection.map(function(image) { var bt = addBrightnessTemp(image); var ndvi = addNDVI(image); var em = addEmissivity(ndvi, roi); var lst = addLST(bt, em); return lst; }); var lstMedian = lstCollection.median().clip(roi);

这里addEmissivity和addLST就是前面 2.3 和 2.4 的函数封装。注意median()是在 ImageCollection 上调用，返回一幅 Image。如果研究区只有一景，median()就是那一景本身。

3.3 关键参数表：K1/K2、比辐射率、劈窗系数

这张表建议直接抄进代码注释里，换传感器的时候改对应列就行。

3.4 导出与验证：用 MODIS LST 做交叉检验

GEE 里导出 LST 到 Drive 或者 Asset，用Export.image.toDrive。验证的话，可以调 MODIS 的MOD11A1产品，重采样到 30 米，算差值。注意 MODIS 是 1 公里，跟 Landsat 30 米比的时候，要先把 Landsat 聚合到 1 公里，或者把 MODIS 降尺度，不然尺度不匹配，差值会虚高。

// 导出 LST Export.image.toDrive({ image: lstMedian, description: 'LST_Export', scale: 30, region: roi, maxPixels: 1e13, fileFormat: 'GeoTIFF' }); // 交叉验证：MODIS LST var modis = ee.ImageCollection('MODIS/061/MOD11A1') .filterDate('2023-06-01', '2023-08-31') .filterBounds(roi) .select('LST_Day_1km') .map(function(img) { return img.multiply(0.02).rename('LST_MODIS'); // 缩放因子 }); var modisMedian = modis.median().clip(roi); var diff = lstMedian.subtract(modisMedian).rename('LST_Diff');

MODIS 的LST_Day_1km缩放因子是 0.02，单位开尔文。算完差值之后，看均值是不是在 ±2 K 以内，标准差是不是小于 3 K。如果偏差大，先检查比辐射率，再检查劈窗系数。

4. 避坑与排查：高分辨率 LST 反演里最容易翻车的五件事

4.1 亮温算出来是负数或者几千 K

现象：BT_B10的统计里出现 -200 K 或者 5000 K 这种离谱值。原因：TOA 辐亮度里有 0 或者负值，log()之后变成 -Infinity 或者 NaN，后续运算全乱。解决：在算亮温之前，先对辐亮度做max(0.001)，把 0 和负值抬到一个小正数。另外检查K1.divide(L).add(1)里 L 是不是影像，如果是影像，add(1)是对每个像素加 1，没问题；但如果 L 是常数，就要用ee.Image.constant()包一下。

4.2 水体温度比陆地还高

现象：湖泊、河流的 LST 比周边城区高 3-5 K。原因：NDVI 阈值法里，水体 NDVI < 0.2，被归到裸地类，ε 取 0.97，但水体的实际比辐射率在 0.98-0.99 之间，ε 偏低导致温度偏高。解决：单独把水体提出来，用 NDWI 或者 MNDWI 做掩膜，水体 ε 取 0.985。代码上就是在addEmissivity里加一个waterMask分支。

4.3 山区温度随高程变化不敏感

现象：海拔每升高 1000 米，LST 只降 2-3 K，实际应该降 6-7 K。原因：劈窗算法的大气校正系数是全球拟合的，没有考虑高程引起的大气水汽变化。解决：如果研究区高程差超过 1500 米，建议用当地探空数据重新拟合劈窗系数，或者直接用辐射传输模型（MODTRAN）算大气透过率和上行/下行辐射，走单窗算法。GEE 里做 MODTRAN 不现实，但可以用 NCEP/NCAR 再分析数据的水汽场做近似校正。

4.4 云掩膜之后影像还是花的

现象：QA_PIXEL 掩了云，但 LST 图上还是有块状异常。原因：Landsat 的 QA_PIXEL 对薄云和卷云识别不全，bit 3 只标了厚云。解决：加一个卷云波段（B9）的阈值，B9 > 0.005 的像素也掩掉。另外可以用QA_RADSAT波段检查饱和像元，饱和的热红外波段反演出来温度会偏低。

4.5 导出的时候内存爆了

现象：Export.image.toDrive报错 "User memory limit exceeded"。原因：研究区太大，或者scale设得太小，GEE 在导出前要做一次全图计算。解决：把研究区分块导出，或者先把 LST 聚合到 60 米、100 米再导出。如果一定要 30 米，用tileScale参数，设成 4 或者 8，让 GEE 分块计算。另外maxPixels设成 1e13 是上限，再大也没用。

5. 进阶技巧：用 GEE 的reduceRegions做时序 LST 统计与趋势检验

前面讲的都是单景或者短期合成。如果你要做 2013-2023 年的城市热岛趋势，需要把每年的 LST 中值提出来，再对每个像素做 Mann-Kendall 趋势检验。GEE 里没有现成的 MK 检验，但可以用ee.Reducer.linearFit做线性趋势，再用ee.Reducer.spearman做秩相关，近似替代。

具体做法：先按年合成 LST 中值，得到一个 ImageCollection，然后reduce成多波段影像，每个波段是一年。再用arrayReduce或者reduce沿时间维做线性拟合。

// 按年合成 var years = ee.List.sequence(2013, 2023); var annualLST = ee.ImageCollection(years.map(function(y) { var start = ee.Date.fromYMD(y, 1, 1); var end = start.advance(1, 'year'); return lstCollection.filterDate(start, end).median() .set('year', y); })); // 转成多波段影像 var annualStack = annualLST.toBands(); // 线性趋势：对每个像素，用年份做自变量，LST 做因变量 var trend = annualStack.reduce(ee.Reducer.linearFit()); // trend 的第一个波段是斜率，第二个是截距 var slope = trend.select('scale');

linearFit返回两个波段：scale是斜率，offset是截距。斜率单位是 K/年，正值表示升温趋势。如果要看显著性，可以用ee.Reducer.pearsonsCorrelation算 R²，R² > 0.5 且斜率绝对值 > 0.05 K/年，一般认为是显著趋势。

还有一个技巧：做时序统计的时候，把水体、云、雪掩膜后的无效像元用unmask填成 NaN，再reduce的时候用ee.Reducer.mean()会自动忽略 NaN。但linearFit对 NaN 的处理不一样，它会把 NaN 当 0 参与拟合，所以要先unmask成ee.Image.constant(-999)，再在拟合前用updateMask把 -999 掩掉。

我自己的习惯是：每次换研究区，先跑一景，把 LST 的直方图拉出来看，正常范围应该在 250-350 K 之间，超过这个范围先查比辐射率和亮温。另外，GEE 的median()合成对薄云残余不敏感，但对季节性差异大的区域，最好分季节合成，不然夏天和秋天的 LST 混在一起，趋势检验会失真。希望帮到你。

本文还有配套的精品资源，点击获取