企业官网建设流程全解析

简介：该资源面向遥感、地理信息与生态环境方向的开发者与研究者，提供一套基于Google Earth Engine的高分辨率地表温度反演源码方案。针对高分辨率热红外数据稀缺的痛点，方案融合Landsat 8与Sentinel-2多源数据，借助随机森林回归建模与残差融合技术，将30米LST降尺度至10米，可用于城市热岛分析与生态环境监测等场景。压缩包共4个文件，约9KB，包含核心的gee_lst_inversion.js脚本、index.html页面、.inscode工程配置及.gitignore，覆盖数据加载、预处理、指数计算、模型训练、残差处理到结果验证的完整流程。目前已有170人学习下载。读者可据此复现整套反演思路，理解多源数据融合与降尺度建模的关键环节，并直接复用脚本进行二次开发与实验验证。

1. 高分辨率LST反演：从GEE源码到可复现的工程化落地

做热环境分析、城市气候研究或者农业旱情监测的同行，大概率都绕不开一个需求：拿到一块地、一座城、一个流域的地表温度（Land Surface Temperature，LST），而且分辨率要够高，高到能看清街区、地块、田埂这一级。Landsat 系列 30 米热红外、MODIS 1 公里，前者空间分辨率够但重访周期长、云干扰大，后者时间密但太糊。GEE（Google Earth Engine）高分辨率 LST 反演源码，解决的正是这个矛盾：用云端算力把 Landsat 8/9 TIRS 热红外波段做大气校正、比辐射率估算、单窗或劈窗反演，输出 30 米级 LST 产品，并且整套流程以源码形式可复现、可批量、可迁移到自己的研究区。这篇笔记面向已经会用 GEE 基础 API、但卡在“反演结果忽高忽低、参数不知道怎么定、源码跑不通”的从业者，把选型理由、逐步操作、参数边界和踩坑记录一次讲透。

2. 反演路线怎么选：单窗、劈窗还是温度发射率分离

2.1 三条主流路线的适用边界

在 GEE 上做高分辨率 LST，绕不开三条路线：单窗算法（Mono-Window）、劈窗算法（Split-Window）、温度发射率分离（TES）。选错路线，后面参数调得再细也是白搭。

单窗算法（如 Jiménez-Muñoz 或 Qin 模型）只需要一个热红外波段，Landsat 8/9 的 Band 10 正好合适。它的核心输入是大气水汽含量、地表比辐射率、近地表气温。优点是公式简单、源码短、GEE 上跑得快；缺点是大气水汽和气温这两个参数如果拿不到实测或再分析数据，误差直接传导到 LST，典型 RMSE 在 1.5～2.5 K。适合研究区小、有气象站点辅助、或者只做相对温度对比的场景。

劈窗算法用两个热红外波段（Landsat 8/9 的 Band 10 和 Band 11），利用两个波段对水汽吸收的差异来反演。理论上不需要外部气温输入，对大气水汽的敏感性也低。但 Landsat 8 的 Band 11 从 2013 年发射起就存在杂散光（stray light）问题，定标精度不稳定，很多研究直接弃用。Landsat 9 的 Band 11 有所改善，但历史数据一致性仍是问题。劈窗在 GEE 上做高分辨率 LST，源码复杂度比单窗高，且 Band 11 的噪声会让结果出现条带。

TES 走的是另一条路：先用地表分类或 NDVI 阈值法估算比辐射率，再用迭代方式同时求解温度和发射率。它对发射率空间异质性大的区域（城市、矿区、混合地表）更稳，但迭代次数、收敛阈值需要调，GEE 上逐像素迭代开销大，不适合大区域批量。

我一般会这样选：如果研究区是城市或地表类型复杂，优先单窗 + NDVI 阈值法估算比辐射率，因为源码可控、参数可查、结果可解释；如果研究区跨越大量水体或沙漠，发射率差异极大，再考虑 TES 或引入外部发射率产品。

2.2 为什么高分辨率场景下我更倾向单窗

高分辨率意味着 30 米像元，一个像元里可能同时有屋顶、路面、树冠。劈窗算法假设两个波段看到的是同一地表，但在 30 米尺度上，Band 10 和 Band 11 的点扩散函数略有差异，混合像元下这个假设会松动。单窗算法虽然依赖外部气温，但气温可以用 ERA5-Land 再分析数据在 GEE 里直接调用，空间分辨率 9 公里，插值到 30 米虽然粗糙，但对 LST 的影响是系统性偏差，可以通过后期与站点回归校正。

另一个现实原因是源码可维护性。单窗算法的源码通常 100～200 行 JavaScript，变量清晰，改一个参数就能重跑；劈窗源码往往要处理 Band 11 的异常值、条带修复、两个波段的配准，调试成本高。对于需要快速出图、批量导出、反复迭代的研究场景，单窗是更务实的起点。

2.3 在 GEE 里跑通单窗反演的最小代码骨架

下面这段代码是单窗反演的核心骨架，我把它拆成“加载影像 → 计算 NDVI 和植被覆盖度 → 估算比辐射率 → 计算大气水汽和气温 → 单窗公式 → 导出”六步。先看代码，再逐段解释参数。

// 1. 定义研究区和时间范围 var roi = ee.Geometry.Rectangle([116.0, 39.5, 116.8, 40.2]); // 示例区域 var startDate = '2023-06-01'; var endDate = '2023-08-31'; // 2. 加载 Landsat 8/9 Collection 2 Level 2 地表反射率和热红外 var landsat = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterBounds(roi) .filterDate(startDate, endDate) .filter(ee.Filter.lt('CLOUD_COVER', 20)) .map(function(img) { // 云掩膜：使用 QA_PIXEL 波段 var qa = img.select('QA_PIXEL'); var cloudMask = qa.bitwiseAnd(1 << 3).eq(0) .and(qa.bitwiseAnd(1 << 4).eq(0)); return img.updateMask(cloudMask); }); // 3. 计算 NDVI 和植被覆盖度 FVC var addNDVI = function(img) { var nir = img.select('SR_B5').multiply(0.0000275).add(-0.2); var red = img.select('SR_B4').multiply(0.0000275).add(-0.2); var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI'); // FVC：用 NDVI 归一化到 [0,1] var fvc = ndvi.subtract(0.05).divide(0.7 - 0.05).clamp(0, 1).rename('FVC'); return img.addBands(ndvi).addBands(fvc); }; // 4. 估算比辐射率：NDVI 阈值法 var addEmissivity = function(img) { var fvc = img.select('FVC'); var ndvi = img.select('NDVI'); // 水体、植被、裸土的经验发射率 var emWater = 0.991; var emVeg = 0.985; var emSoil = 0.960; // 混合像元发射率 var em = ee.Image(0).expression( 'fvc * emVeg + (1 - fvc) * emSoil', { 'fvc': fvc, 'emVeg': emVeg, 'emSoil': emSoil }).rename('EM'); // 水体单独赋值 var waterMask = ndvi.lt(0); em = em.where(waterMask, emWater); return img.addBands(em); }; // 5. 计算亮温（Brightness Temperature） var addBT = function(img) { var rad = img.select('ST_B10').multiply(0.00341802).add(149.0); // 定标到辐射亮度 var bt = rad.divide(1).multiply(1).rename('BT'); // 简化：实际需用普朗克反函数 // 更准确的做法：直接用 ST_B10 的定标系数转亮温 var btReal = img.select('ST_B10').multiply(0.00341802).add(149.0) .divide(1).rename('BT_real'); return img.addBands(btReal); }; // 6. 单窗公式（Qin 模型简化版） var singleWindow = function(img) { var bt = img.select('BT_real'); // 亮温，单位 K var em = img.select('EM'); // 大气水汽含量 w（g/cm²），示例用固定值，实际应从 ERA5 获取 var w = 1.5; // 近地表气温 Ta（K），示例用固定值 var ta = 300; // Qin 模型系数（根据水汽范围查表） var tau = 0.75; // 大气透射率 var a = -67.355; var b = 0.458; var lst = bt.add(1).multiply(1).rename('LST'); // 占位，实际公式见下 // 实际单窗公式：Ts = [a*(1-C-D) + (b*(1-C-D)+C+D)*Tb - D*Ta] / C var c = em.multiply(tau); var d = (1 - tau).multiply(1 + (1 - em).multiply(tau)); var lstReal = bt.expression( '(a*(1-C-D) + (b*(1-C-D)+C+D)*BT - D*TA) / C', { 'a': a, 'b': b, 'C': c, 'D': d, 'BT': bt, 'TA': ta }).rename('LST'); return img.addBands(lstReal); }; // 7. 执行并导出 var processed = landsat.map(addNDVI).map(addEmissivity).map(addBT).map(singleWindow); var lstMean = processed.select('LST').mean().clip(roi); Export.image.toDrive({ image: lstMean, description: 'LST_Export', folder: 'GEE_LST', scale: 30, region: roi, maxPixels: 1e13 });

这段代码的逻辑链条是：先做云掩膜保证输入干净，再用 NDVI 推 FVC，用 FVC 推比辐射率，用热红外波段转亮温，最后套单窗公式。参数说明如下：

• CLOUD_COVER < 20：云量阈值，城市研究可以放宽到 30，但水体研究建议压到 10。
• 0.05和0.7：FVC 归一化的 NDVI 上下限，分别对应裸土和全植被。如果研究区植被稀疏，下限可以调到 0.1。
• emWater = 0.991、emVeg = 0.985、emSoil = 0.960：经验发射率，来自 ASTER 光谱库的典型值。矿区或特殊地表需要替换。
• w = 1.5：大气水汽含量，单位 g/cm²。这是单窗算法最敏感的参数，固定值只适合小区域短时间。实际做法是用 ERA5-Land 的total_column_water_vapour波段按影像日期匹配。
• tau = 0.75：大气透射率，与水汽相关。水汽 1.5 时约 0.75，水汽 3.0 时降到 0.6 左右。
• a = -67.355、b = 0.458：Qin 模型系数，适用于 0.4～2.0 g/cm² 水汽范围。超出范围要换系数表。

注意：上面代码里addBT函数做了简化，实际 Landsat Collection 2 Level 2 的ST_B10已经过定标，转亮温的公式是BT = ST_B10 * 0.00341802 + 149.0，单位开尔文。不要再用普朗克反函数重复计算。

2.4 参数敏感性：哪些参数动一点结果就翻车

单窗算法里，大气水汽w的敏感性最高。水汽从 1.5 变到 2.5，LST 可能差 1～2 K。比辐射率em次之，从 0.96 变到 0.99，LST 差 0.5～1 K。近地表气温Ta的影响相对小，但如果是昼夜温差大的区域，用日均温代替瞬时气温会引入 1 K 以上偏差。

我的做法是：水汽用 ERA5-Land 逐小时数据，按影像过境时间（Landsat 约当地时间 10:30）取最近时刻；气温用同一数据源的 2 米气温；比辐射率用 NDVI 阈值法，但在城市核心区叠加不透水面比例做修正。这样一套下来，与地面站点对比的 RMSE 能压到 1.5 K 以内。

3. 从源码到批量出图：GEE 工程化改造的四个动作

3.1 把固定参数换成逐影像动态参数

上面代码里w和ta是固定值，只能跑单景。要批量处理整个夏季，必须把这两个参数变成逐影像的。GEE 里可以用ee.ImageCollection的map配合join来实现：先加载 ERA5-Land 的逐小时数据，按影像的system:time_start做时间匹配，取最近时刻的水汽和气温。

// 加载 ERA5-Land 逐小时数据 var era5 = ee.ImageCollection('ECMWF/ERA5_LAND/HOURLY') .filterDate(startDate, endDate) .select(['total_column_water_vapour', 'temperature_2m']); // 为每景 Landsat 匹配最近的 ERA5 影像 var matchERA5 = function(img) { var time = img.get('system:time_start'); var era5Img = era5.filterDate( ee.Date(time).advance(-1, 'hour'), ee.Date(time).advance(1, 'hour') ).first(); var w = era5Img.select('total_column_water_vapour'); // 单位 kg/m²，约等于 g/cm² 的 0.1 倍 var ta = era5Img.select('temperature_2m'); // 单位 K return img.addBands(w.rename('W')).addBands(ta.rename('TA')); };

这里有个单位坑：ERA5-Land 的total_column_water_vapour单位是 kg/m²，换算成 g/cm² 要除以 10。也就是 15 kg/m² 对应 1.5 g/cm²。很多源码直接拿过来用，结果 LST 偏高 3～5 K，就是单位没换。

3.2 用 QA 波段做逐像元云掩膜而不是整景过滤

filter(ee.Filter.lt('CLOUD_COVER', 20))是整景过滤，会丢掉很多有局部云的影像。高分辨率场景下，一景里可能只有一小块云，整景丢掉太浪费。正确做法是用QA_PIXEL波段做逐像元掩膜，把云、云影、卷云都去掉。

var maskClouds = function(img) { var qa = img.select('QA_PIXEL'); // bit 3: cloud, bit 4: cloud shadow, bit 5: snow, bit 6: cirrus var mask = qa.bitwiseAnd(1 << 3).eq(0) .and(qa.bitwiseAnd(1 << 4).eq(0)) .and(qa.bitwiseAnd(1 << 5).eq(0)) .and(qa.bitwiseAnd(1 << 6).eq(0)); return img.updateMask(mask); };

QA_PIXEL是位掩码，1 << 3表示第 3 位。Landsat Collection 2 的位定义是：bit 0 填充、bit 1 膨胀、bit 2 云、bit 3 云影、bit 4 雪、bit 5 卷云。不同版本位定义可能微调，跑之前一定去查对应 Collection 的文档。

3.3 导出时的尺度、投影和分块策略

GEE 导出 LST 时，scale: 30是必须的，因为 Landsat 热红外原始分辨率是 100 米，Collection 2 重采样到 30 米，但实际有效分辨率还是 100 米。导出 30 米只是保持与反射率波段一致，不要误以为热红外真的到了 30 米。

maxPixels: 1e13是防止大区域导出报错。如果研究区超过一个 Landsat 景，建议用tileScale参数分块，或者按景导出后本地拼接。crs默认是 EPSG:4326，如果要做投影分析，改成EPSG:32650之类的 UTM 投影。

3.4 结果验证：与站点数据对比的三种方式

没有验证的 LST 反演等于没做。GEE 里可以用三种方式验证：一是加载地面站点的实测 LST（如果有），用reduceRegions提取对应像元值做回归；二是与官方 LST 产品（如 MODIS MOD11A1）对比，看空间格局是否一致；三是与 Landsat Collection 2 Level 2 自带的ST_B10产品对比，看自己的反演是否比官方产品更贴近站点。

我一般会先做第三种，因为ST_B10是官方定标产品，虽然没做大气校正，但可以作为基准。如果自己的 LST 与ST_B10的空间相关系数低于 0.8，说明反演流程有问题，先查云掩膜和发射率。

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

4.1 现象：LST 结果出现大面积负值或超过 350 K

原因：亮温转公式用错，或者比辐射率小于 0.9。Landsat Collection 2 的ST_B10已经定标，直接乘系数加偏移就是亮温，不需要再做普朗克反函数。如果用了 Level 1 的B10波段，才需要自己定标。

解决：确认用的是LANDSAT/LC08/C02/T1_L2而不是T1。检查EM波段的最小值，如果低于 0.9，说明 NDVI 阈值法的水体掩膜没生效，或者 FVC 计算时 NDVI 上下限设反了。

4.2 现象：城市区域 LST 比周边农田还低

原因：比辐射率赋值反了。城市不透水面发射率约 0.95～0.96，植被约 0.98～0.99。如果代码里把植被发射率赋给了裸土，城市像元会被低估。

解决：打印EM波段的直方图，确认城市核心区的值在 0.95 左右，农田在 0.98 左右。如果反了，检查emVeg和emSoil的赋值顺序。

4.3 现象：同一区域不同日期的 LST 差异超过 10 K

原因：大气水汽参数没有逐影像匹配，用了固定值。夏季水汽从 1.0 到 3.0 变化，LST 差异可达 5 K 以上。另外，如果影像过境时间不同（Landsat 8 和 9 差几分钟），气温也会有差异。

解决：用 ERA5-Land 逐小时数据按影像时间匹配，时间窗口控制在 ±1 小时。如果研究区跨越多个 Landsat 行列号，确保每景都单独匹配。

4.4 现象：导出影像有大量空值条带

原因：云掩膜太激进，或者QA_PIXEL位定义用错。Landsat 8 和 9 的 QA 位定义基本一致，但 Collection 1 和 Collection 2 不同。如果用了 Collection 1 的位定义去掩膜 Collection 2 数据，会把大量有效像元误判为云。

解决：确认数据源是C02，位定义查 Collection 2 文档。云掩膜后打印有效像元比例，如果低于 30%，说明掩膜过头了。

4.5 现象：GEE 任务运行超时或内存溢出

原因：研究区太大，或者map函数里做了太多逐像元计算。单窗算法本身不重，但如果叠加了 TES 迭代或高分辨率分类，计算量会爆炸。

解决：用tileScale参数降低并行粒度，或者把研究区拆成多个小区域分别导出。maxPixels不要设得过大，1e13 是上限，实际按研究区面积算。如果还是超时，把map里的计算拆成多个Export任务，分步执行。

5. 进阶技巧：用查找表把单窗算法提速三倍

单窗算法的公式里，C和D只依赖比辐射率em和大气透射率tau。em的范围是 0.90～1.00，tau的范围是 0.5～0.9。如果把这两个参数离散化成 0.005 的步长，C和D的组合只有 20×80=1600 种。可以预先算好一张查找表，运行时直接查表，避免逐像元做浮点运算。

// 构建查找表：em 从 0.90 到 1.00，tau 从 0.5 到 0.9 var emList = ee.List.sequence(0.90, 1.00, 0.005); var tauList = ee.List.sequence(0.5, 0.9, 0.01); var lut = emList.map(function(em) { return tauList.map(function(tau) { var c = ee.Number(em).multiply(tau); var d = ee.Number(1).subtract(tau) .multiply(ee.Number(1).add(ee.Number(1).subtract(em).multiply(tau))); return ee.Feature(null, {'em': em, 'tau': tau, 'C': c, 'D': d}); }); }).flatten(); var lutFC = ee.FeatureCollection(lut); // 运行时：把 em 和 tau 离散化后 join 查找表 var applyLUT = function(img) { var em = img.select('EM'); var tau = img.select('TAU'); // 假设已有 tau 波段 // 离散化到查找表步长 var emKey = em.multiply(200).round().divide(200); var tauKey = tau.multiply(100).round().divide(100); // 用 join 或 remap 查表（此处省略 join 细节） return img; };

查找表的精度损失在 0.1 K 以内，但计算速度能提升 2～3 倍。对于大区域批量导出，这个优化很值。如果研究区固定，还可以把查找表导出成 CSV，本地用 Python 做后处理时直接查。

另一个技巧是分块导出后本地拼接。GEE 导出大区域容易超时，按 1°×1° 分块导出，本地用 GDAL 的gdal_merge.py拼接。拼接时注意重叠区用均值或最大值，避免接缝。

我自己的习惯是：每次跑完 LST，先导出EM、BT、LST三个波段，本地用 Python 画散点图看EM和LST的关系。如果EM在 0.95 附近的像元 LST 明显偏高，说明大气校正过头了，回头调tau。这个习惯帮我省了很多次重新跑 GEE 任务的时间。希望帮到你。

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