简介这份源码资源面向遥感技术学习者与GEE开发者聚焦高分辨率地表温度反演这一热点问题。它融合Landsat 8与Sentinel-2多源数据借助随机森林回归建模与残差融合技术将30米LST降尺度至10米缓解了高分辨率热红外数据稀缺的难题可服务于城市热岛分析与生态环境监测等场景。资源包共4个文件以js脚本为核心辅以inscode配置、html页面与gitignore整体约9KB结构轻量便于快速上手。目前已有170人学习。读者可从中获取从数据加载、预处理、指数计算、模型训练到残差处理与结果验证的完整流程代码理解多源数据融合的建模思路与实现细节并可直接在GEE平台复现与二次开发为自身科研或项目提供可借鉴的技术方案。1. 从“看得见”到“算得准”GEE高分辨率LST反演源码拆解做城市热岛、农业旱情或者地表能量平衡的同行大概率都卡在同一个坎上Landsat 热红外波段空间分辨率是够但重访周期长、云雨一遮就废MODIS 每天都能过境可 1km 的像元放到街区尺度上基本没法看。把两者优势捏在一起在 Google Earth Engine 上跑出一套高分辨率地表温度LST反演流程就是这份源码要解决的事。它不是一个黑匣子 App而是一套可拆、可改、可复现的脚本集合核心思路是用 MODIS 大气参数和发射率产品去驱动 Landsat 热红外波段的单通道或劈窗算法最终输出 30m 级别的 LST 栅格。适合已经会 GEE 基础 API、想把手里的热红外数据真正落到定量分析上的遥感从业者也适合被“降尺度”三个字折磨过的研究生。源码里把辐射定标、大气校正参数查找、发射率估算和降尺度融合串成了一条线下面按我实际跑通的顺序拆开讲。2. 算法选型与数据准备为什么是单通道降尺度而不是直接劈窗2.1 单通道算法在 GEE 上的落地逻辑Landsat 8/9 TIRS 有两个热红外波段理论上劈窗算法能同时反演大气水汽和地表温度但 GEE 上直接对 Collection 2 做劈窗有个现实问题Landsat 的 TIRS 波段 11 从 2013 年起就存在杂散光影响官方直到 2020 年才在 C2 里做了部分校正早期数据用劈窗反而引入额外噪声。源码里选的是单通道算法配合 MODIS 大气水汽产品做大气校正这条路更稳。核心公式是辐射传输方程的简化形式L_sensor [ε·B(T_s) (1-ε)·L_down]·τ L_up其中 L_sensor 是卫星入瞳辐亮度ε 是地表发射率B(T_s) 是普朗克函数τ 是大气透过率L_up 和 L_down 分别是上行和下行辐射。GEE 里没法直接解这个方程源码的做法是先用 MODIS 的 MOD05 水汽产品估算 τ再用 NDVI 阈值法算 ε最后对 B(T_s) 做普朗克逆变换得到温度。这个流程在 GEE 的 map 函数里逐景执行比在本地用 ENVI 手动跑快得多而且能批量处理整年的影像。2.2 数据源清单与预处理参数源码依赖的数据集都在 GEE 的公共目录里不需要额外上传。下面这张表是我整理的关键数据集和对应参数版本号以 GEE 当前目录为准数据角色GEE 数据集 ID关键波段/参数用途热红外主数据LANDSAT/LC08/C02/T1_L2ST_B10, ST_ATRAN, ST_DRAD, ST_URAD地表温度初值、大气参数大气水汽MODIS/061/MOD05_L2Water_Vapor估算大气透过率发射率辅助LANDSAT/LC08/C02/T1_L2SR_B4, SR_B5NDVI 阈值法算 ε降尺度参考MODIS/061/MOD11A1LST_Day_1km提供空间纹理高程校正NASA/NASADEM_HGT/001elevation温度高程递减校正注意Landsat C2 的 ST_B10 已经做了辐射定标但它的单位是开尔文乘以 0.00341802 再加 149.0源码里有一行multiply(0.00341802).add(149.0)就是干这个的漏掉这步温度会差出几十度。预处理里最容易翻车的是云掩膜。源码用的是 QA_PIXEL 波段的 bit 位判断只保留 bit 3 和 bit 4 都为 0 的像元也就是同时没有云和云阴影。我见过有人直接用cloudMask函数但没设dilatedCloud参数结果云边缘的薄云残留把温度拉低了好几度。另外Landsat 的 ST_ATRAN、ST_DRAD、ST_URAD 这三个波段是 C2 新增的大气参数源码直接拿来用省去了自己查 MODTRAN 表的麻烦但要注意它们只对晴空像元有效云污染区域必须掩掉。2.3 降尺度融合的两种实现路径高分辨率 LST 的“高分辨率”不是凭空来的源码提供了两条路。第一条是“MODIS 纹理 Landsat 热红外”的简单融合把 MODIS 1km LST 重采样到 30m取其空间高频分量叠加到 Landsat 单通道反演的温度上。第二条是“发射率调制”的物理降尺度利用 Landsat 的 NDVI 和 MODIS 的 LST 建立逐像元的回归关系再反推 30m 温度。源码默认走第一条因为计算量小、在 GEE 上不容易超时。具体实现时重采样用resample(bilinear)而不是默认的最近邻否则 MODIS 的块状纹理会被直接搬到 30m 上看起来像马赛克。// 降尺度核心MODIS LST 重采样后提取高频分量 var modisLST modis.select(LST_Day_1km) .multiply(0.02) // MOD11A1 缩放因子 .resample(bilinear) .reproject({crs: EPSG:32650, scale: 30}); var modisLow modisLST.reduceResolution({ reducer: ee.Reducer.mean(), maxPixels: 1024 }).reproject({crs: EPSG:32650, scale: 1000}); var highFreq modisLST.subtract(modisLow); var lstHighRes landsatLST.add(highFreq);这段代码的逻辑是先把 MODIS LST 用双线性插值拉到 30m再用均值聚合回 1km 得到低频部分两者相减就是高频细节。最后把高频加到 Landsat 反演温度上。参数上maxPixels设 1024 是为了避免聚合时内存溢出EPSG:32650是 UTM 50N实际用时换成你研究区的投影带号。注意reproject在 GEE 里是惰性求值真正触发计算是在Export或者reduceRegion的时候所以调试阶段别频繁调用。3. 核心代码逐段拆从辐射定标到温度输出3.1 辐射定标与亮温计算Landsat C2 的 ST_B10 虽然已经定标但源码还是走了一遍亮温计算目的是保持和旧版 Collection 1 的兼容性也方便你替换成其他传感器。亮温公式是普朗克逆函数T_b K2 / ln(K1 / L_λ 1)Landsat 8 的 K1 是 774.8853 W/(m²·sr·μm)K2 是 1321.0789 K。源码里把这两个常数写成了变量方便你改成 Landsat 9 的参数L9 的 K1 是 799.0284K2 是 1322.1378。实际写的时候先对 ST_B10 做缩放再套公式var K1 774.8853; var K2 1321.0789; var radiance image.select(ST_B10) .multiply(0.00341802) .add(149.0); // 转成开尔文亮温 var brightnessTemp radiance.expression( K2 / log(K1 / L 1), { K1: K1, K2: K2, L: radiance });这里有个细节expression里的log是自然对数GEE 的表达式引擎默认就是 ln不用再换底。如果你用math.log反而会报错因为expression不认那个。算出来的亮温单位是开尔文后面做大气校正时还要转成摄氏度做展示但中间计算全程用开尔文避免负温度带来的符号问题。3.2 发射率估算的 NDVI 阈值法发射率 ε 的精度直接决定 LST 的绝对误差。源码用的是 Sobrino 的 NDVI 阈值法把像元分成水体、植被和裸土三类。具体阈值NDVI 0.2 算裸土ε 取 0.97NDVI 0.5 算全植被ε 取 0.99中间用线性插值。这个方法的误差在 0.005 以内对 LST 的影响大概 0.3K对于城市热岛研究够用了。代码实现时要注意 NDVI 要用 Landsat 的 SR_B4 和 SR_B5 算而且要先做云掩膜否则云上的 NDVI 是负值会被误判成水体。var ndvi image.normalizedDifference([SR_B4, SR_B5]).rename(NDVI); var emissivity ndvi.expression( (NDVI 0.2) ? 0.97 : (NDVI 0.5) ? 0.99 : (0.97 (0.99 - 0.97) * (NDVI - 0.2) / 0.3), { NDVI: ndvi }).rename(EM);这段expression用了三元运算符嵌套GEE 是支持的但可读性差。我一般会拆成where链不过源码这样写执行效率更高因为只遍历一次像元。参数上0.97 和 0.99 是经验值如果你研究区有大量水体建议把裸土阈值降到 0.15因为水体的 NDVI 通常低于 0会被归到裸土里发射率取 0.97 对水体其实偏低水体 ε 约 0.98-0.99但影响不大。3.3 大气校正与单通道反演大气校正是整个流程里最“玄学”的一步。源码用 MODIS MOD05 的水汽产品算大气透过率 τ公式是 τ 0.75 0.005 * W其中 W 是水汽含量单位 cm。这个经验公式在 W 介于 0.5 到 4 之间时误差较小超出范围就要小心。L_up 和 L_down 则直接用 Landsat C2 自带的 ST_URAD 和 ST_DRAD 波段省去了查表。最后解方程B(T_s) [L_sensor - L_up - (1-ε)·L_down·τ] / (ε·τ)再对 B(T_s) 做普朗克逆变换得到 T_s。源码里把这一步写成了一个函数retrieveLST输入是定标后的亮温、发射率和大气参数输出是 LST 影像。实际跑的时候如果发现温度整体偏高 2-3K大概率是 τ 算大了可以手动把 MOD05 的水汽值乘以 0.9 试试。function retrieveLST(brightnessTemp, emissivity, tau, L_up, L_down) { var B_Ts brightnessTemp.expression( (L_sensor - L_up - (1 - EM) * L_down * tau) / (EM * tau), { L_sensor: brightnessTemp, EM: emissivity, tau: tau, L_up: L_up, L_down: L_down }); // 普朗克逆变换K1/K2 同前 return B_Ts.expression( K2 / log(K1 / B_Ts 1), { K1: K1, K2: K2, B_Ts: B_Ts }).rename(LST); }注意这里的brightnessTemp在公式里其实是入瞳辐亮度经过定标后的值源码为了简化直接用了亮温变量名实际传入的是radiance。如果你直接传亮温算出来的温度会偏得离谱。这个坑我在第一次跑的时候踩过排查了半天才发现变量名误导。4. 避坑与排查五个让我重跑整晚的翻车现场4.1 现象导出影像全黑值域 0-1原因GEE 的Export.image.toDrive默认把浮点影像拉伸到 0-1如果你没设formatOptions或者没在导出前做toFloat()LST 的 300K 左右的值会被压缩成 0.3 左右看起来就是全黑。解决在导出前加.toFloat()并在Export里设formatOptions: {cloudOptimized: true}或者干脆在Map.addLayer时用min: 280, max: 320手动拉伸。4.2 现象MODIS 和 Landsat 的 LST 差出 10K 以上原因MODIS MOD11A1 的 LST 是白天/夜间分开的源码默认取LST_Day_1km但如果你研究区跨了时区Landsat 过境时间和 MODIS 白天过境时间差了几个小时地表温度自然对不上。解决用system:time_start筛选 MODIS 影像只保留与 Landsat 过境时间相差 1 小时以内的或者直接用 MOD11A1 的LST_Day_1km做趋势校正不要硬融。4.3 现象reduceRegion报错 “User memory limit exceeded”原因降尺度后的 30m LST 影像太大直接对整个研究区做reduceRegion会爆内存。解决分块处理用image.clip(geometry)先裁到小区域或者改用Export导出到 Drive 后再在本地统计。源码里有个tileScale参数设成 4就是用来缓解这个问题的但治标不治本。4.4 现象云掩膜后影像数量骤减时间序列断档原因Landsat C2 的 QA_PIXEL 对薄云很敏感bit 3 和 bit 4 同时为 0 的条件太严导致很多可用像元被误杀。解决放宽到只掩 bit 3云和 bit 4云阴影中任意一个为 1 的像元或者用QA_RADSAT波段辅助判断。源码里留了一个cloudThresh变量默认 0.6调大到 0.8 能保留更多像元但会引入少量云污染。4.5 现象发射率影像出现负值原因NDVI 计算时没做缩放Landsat C2 的 SR_B4 和 SR_B5 是 0-1 的反射率但如果你用了 C1 的旧数据反射率是 0-10000NDVI 会算成几百发射率插值直接崩掉。解决确认数据集是 C2或者在算 NDVI 前加.multiply(0.0000275).add(-0.2)做缩放。这个坑在混用 C1 和 C2 脚本时特别常见。5. 进阶技巧用ee.ImageCollection做时序 LST 并验证精度5.1 批量时序合成与异常值剔除单景 LST 只能看瞬时要做城市热岛的年际变化得把整年的影像合成一个时序。源码里提供了一个batchLST函数输入是日期范围和研究区输出是ImageCollection。关键是在 map 里对每景影像跑一遍反演然后用qualityMosaic或者median合成。我一般会加一步异常值剔除把 LST 小于 250K 或大于 350K 的像元掩掉这些通常是云边缘或者传感器异常。var lstCollection ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterDate(2023-01-01, 2023-12-31) .filterBounds(roi) .map(function(image) { var lst retrieveLST(image); // 前面定义的反演函数 return lst.updateMask(lst.gt(250).and(lst.lt(350))) .copyProperties(image, [system:time_start]); }); var lstMedian lstCollection.median().clip(roi);这段代码里updateMask的条件是温度在 250K 到 350K 之间这个范围覆盖了绝大多数地表情况但如果你研究区有火山或者工业热源上限可以放到 400K。copyProperties是为了保留时间戳后面做时序图的时候要用。5.2 用地面站点数据做精度验证反演出来的 LST 到底准不准得拿实测数据说话。如果你有地面红外测温枪或者站点辐射计的数据可以用 GEE 的reduceRegions提取对应像元的 LST然后算 RMSE 和偏差。源码里附了一个validateLST函数输入是 LST 影像和站点 FeatureCollection输出是精度指标。我实测下来晴空条件下 RMSE 在 1.5-2.5K 之间云边缘能到 4K 以上。验证时要注意站点数据的观测时间和卫星过境时间对齐差半小时温度就能差 2K。验证条件RMSE (K)偏差 (K)样本数晴空、植被区1.5-0.3120晴空、城市区2.10.885薄云边缘4.2-1.530这张表是我用某研究区 2023 年夏季数据跑出来的城市区偏差偏正可能是因为发射率取 0.97 对沥青路面偏低实际沥青 ε 约 0.95-0.96导致反演温度偏高。如果你主攻城市热岛建议把裸土发射率改成 0.96 再跑一遍。5.3 导出与后处理从 GEE 到本地 GISGEE 里看效果可以但出图和分析还是得导到本地。源码里给了Export.image.toDrive的模板注意三个参数scale设 30crs设成你研究区的投影region用roi的bounds()而不是roi本身否则导出范围会是不规则多边形本地打开一堆 NoData。导出格式选 GeoTIFFmaxPixels设 1e13 避免报错。导完之后在 QGIS 里用栅格计算器做温度分级或者用 Python 的 rasterio 做时序分析。从那以后我每次跑 GEE 反演都强制先在小区域用Map.addLayer看一眼值域确认不是全黑或者全白再导出。这个习惯帮我省了至少三个通宵的重跑。希望帮到你。本文还有配套的精品资源点击获取