基于 GEE 的 Landsat-Sentinel 协同逐像元线性回归与残差校正的地表温度（LST）降尺度研究

目录

一、方法原理概述

（一）为什么需要降尺度

（二）技术路线

二、数据预处理与质量控制

（一）Landsat数据筛选与LST计算

（二）Landsat质量掩膜（QA_PIXEL）

（三）Sentinel-2数据筛选与SCL掩膜

（四）清晰像元比例（Clear Fraction）与影像评分

三、光谱指数计算

（一）Landsat光谱指数

（二）Sentinel-2光谱指数

四、多元线性回归建模

（一）自变量选取与影像构建

（二）线性回归系数求解

（三）30米拟合LST与10米预测LST

五、残差上采样与融合

（一）残差计算

（二）高斯卷积与上采样

（三）最终降尺度LST

六、精度评估指标

（一）回归拟合指标

（二）降尺度前后统计对比

七、结果导出与可视化

（一）地图可视化

（二）数据导出

八、总结

九、运行结果

若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！

高分辨率热红外数据稀缺是城市热岛精细分析、生态环境监测等领域长期面临的痛点——Landsat系列提供30米分辨率的热红外波段，但空间分辨率不足以支撑精细尺度分析；Sentinel-2虽拥有10米空间分辨率，却缺乏热红外波段。

本文基于Google Earth Engine（GEE）平台，系统解析一套将Landsat 8/9的30米地表温度（LST）降尺度至10米分辨率的技术方案。该方案通过多元线性回归建模与残差融合技术两大核心步骤，利用Sentinel-2的高分辨率光谱指数（NDVI、NDBI、NDWI）及DEM、坡度等辅助变量，实现LST空间分辨率的提升。文章将从数据预处理、指数计算、回归建模、残差上采样到精度评估逐层展开，兼顾理论原理解析与GEE实操避坑要点，适合遥感领域学习者和科研人员阅读参考。

一、方法原理概述

（一）为什么需要降尺度

Landsat 8/9的Level-2数据产品提供了经过大气校正的地表温度（LST）波段（ST_B10），空间分辨率为30米。然而，在城市热岛效应研究、微气候分析等场景中，30米分辨率往往不足以捕捉建筑物阴影、街谷、水体边界等细节尺度的温度差异。Sentinel-2多光谱仪（MSI）虽无热红外通道，但其可见光/近红外波段空间分辨率达10米，且重访周期短（5天），是理想的“辅助数据源”。

降尺度的核心逻辑是：建立LST与高分辨率辅助变量之间的统计关系，然后将该关系应用于更高分辨率的辅助变量上，从而“推断”出高分辨率的LST分布。

（二）技术路线

本方案采用的技术路线可概括为以下五个步骤：

二、数据预处理与质量控制

（一）Landsat数据筛选与LST计算

Landsat数据源选用Collection 2 Level-2产品（LANDSAT/LC08/C02/T1_L2和LANDSAT/LC09/C02/T1_L2），该级别产品已包含经大气校正的地表反射率和地表温度，大幅降低了预处理门槛。

var landsatCollection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
.filterBounds(region)
.filterDate(startDate, endDate)
.filter(ee.Filter.eq('PROCESSING_LEVEL', 'L2SP'))
.filter(ee.Filter.lte('CLOUD_COVER', landsatCloud));

• PROCESSING_LEVEL = 'L2SP'：仅保留经大气校正的Level-2科学产品；

• CLOUD_COVER ≤ 70%：云量阈值可根据研究区实际情况调整——云量过高会导致LST数据大面积缺失，过低则可能无可用影像；

• 时间范围：startDate至endDate，建议选择地表温度差异显著的季节（如北半球夏季）。

Landsat C2L2产品的ST_B10波段需通过官方缩放公式转换为摄氏温度：

其中DN为ST_B10波段的像元值，计算结果单位为℃。对应代码实现：

var lst = image.select('ST_B10')
.multiply(0.00341802)
.add(149)
.subtract(273.15)
.rename('LST');

光学波段（SR_B2~SR_B7）的反射率缩放公式为：

（二）Landsat质量掩膜（QA_PIXEL）

Landsat Level-2产品提供QA_PIXEL波段用于质量评估，需对以下条件进行掩膜：

• 填充像素（bit 0）、 dilated cloud（bit 1）、 cirrus（bit 2）、 cloud（bit 3）、 cloud shadow（bit 4）均需排除；

• QA_RADSAT波段需为0（无辐射饱和）。

var mask = qa.bitwiseAnd(1 << 0).eq(0)
.and(qa.bitwiseAnd(1 << 1).eq(0))
.and(qa.bitwiseAnd(1 << 2).eq(0))
.and(qa.bitwiseAnd(1 << 3).eq(0))
.and(qa.bitwiseAnd(1 << 4).eq(0))
.and(image.select('QA_RADSAT').eq(0));

（三）Sentinel-2数据筛选与SCL掩膜

Sentinel-2数据源选用COPERNICUS/S2_SR_HARMONIZED，该集合已对Sentinel-2A和2B数据进行辐射归一化处理。关键筛选逻辑是以Landsat影像日期为中心，向前后各扩展maxDifferenceDays天（默认30天），寻找最接近的Sentinel-2影像：

var sentinelCollection = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
.filterBounds(region)
.filterDate(
landsatDate.advance(-maxDifferenceDays, 'day'),
landsatDate.advance(maxDifferenceDays + 1, 'day')
)
.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', sentinelCloud));

Sentinel-2的云掩膜采用SCL（Scene Classification Layer） 波段，排除以下类别：

var mask = scl.neq(0).and(scl.neq(1)).and(scl.neq(3))
.and(scl.neq(8)).and(scl.neq(9))
.and(scl.neq(10)).and(scl.neq(11));

（四）清晰像元比例（Clear Fraction）与影像评分

为从影像集合中自动选出最优影像，代码引入清晰像元比例和影像评分两个指标：

• 清晰像元比例：有效像元（非掩膜像元）占研究区总面积的比例；

• 影像评分：综合云量、清晰像元比例、日期差异等因素的加权分数，分数越低表示影像质量越好。

var score = ee.Number(1)
.subtract(clearFraction)
.multiply(100)
.add(cloud.multiply(0.1));

最终通过.sort('imageScore').first()选出评分最优的Landsat和Sentinel-2影像对。

三、光谱指数计算

（一）Landsat光谱指数

基于Landsat地表反射率波段，计算以下三种常用光谱指数：

归一化植被指数（NDVI） ：

归一化水体指数（NDWI） ：

归一化建筑指数（NDBI） ：

代码中通过泛化函数实现指数计算，避免分母为零：

function getIndex(band1, band2, name) {
var sum = band1.add(band2);
return band1.subtract(band2)
.divide(sum)
.updateMask(sum.neq(0))
.rename(name);
}

Landsat波段对应关系：SR_B4（Red）、SR_B5（NIR）、SR_B3（Green）、SR_B6（SWIR1）。

（二）Sentinel-2光谱指数

Sentinel-2波段对应关系：B4（Red，10m）、B8（NIR，10m）、B3（Green，10m）、B11（SWIR，20m重采样至10m）。

var sentinelGreen = sentinel.select('B3')
.resample('bilinear')
.reproject({ crs: sentinelProjection, scale: 10 });

注意事项：Sentinel-2的B11（SWIR）原生分辨率为20米，需通过resample()和reproject()重采样至10米，重采样方法选择bilinear（双线性插值）可在保持空间连续性的同时避免过度平滑。

四、多元线性回归建模

（一）自变量选取与影像构建

回归模型的自变量包括：

• 光谱指数：NDVI、NDBI、NDWI（反映植被、建筑、水体等地表覆盖类型）；

• 地形因子：DEM（数字高程模型）、Slope（坡度）；

• 截距项：常数1。

因变量为Landsat LST（30米分辨率）。

回归方程可表示为：

其中为回归残差。

代码通过构建多波段影像，将6个自变量堆叠为一个多波段影像：

var regressionImage = constant
.addBands(landsatNDVI)
.addBands(landsatNDBI)
.addBands(landsatNDWI)
.addBands(dem30)
.addBands(slope30)
.addBands(lst30);

（二）线性回归系数求解

GEE的ee.Reducer.linearRegression()实现了普通最小二乘法（OLS） ，可同时处理多个自变量和一个因变量：

var regression = regressionImage.reduceRegion({
reducer: ee.Reducer.linearRegression({ numX: 6, numY: 1 }),
geometry: region,
scale: 30,
crs: landsatProjection,
maxPixels: 1e13,
tileScale: 4
});

回归结果为一个系数数组，依次对应：截距（b0）、NDVI、NDBI、NDWI、DEM、Slope。通过数组索引提取各系数：

var coefficientArray = ee.Array(regression.get('coefficients'));
var b0 = ee.Number(coefficientList.get(0));
var bNDVI = ee.Number(coefficientList.get(1));
// … 依次提取

（三）30米拟合LST与10米预测LST

利用回归系数，可计算30米尺度的拟合LST（用于残差计算）和10米尺度的预测LST（用于最终降尺度）：

// 30米拟合值
var fittedLST30 = landsatNDVI.multiply(bNDVI)
.add(landsatNDBI.multiply(bNDBI))
.add(landsatNDWI.multiply(bNDWI))
.add(dem30.multiply(bDEM))
.add(slope30.multiply(bSlope))
.add(b0);

// 10米预测值（基于Sentinel-2指数）
var predictedLST10 = sentinelNDVI.multiply(bNDVI)
.add(sentinelNDBI.multiply(bNDBI))
.add(sentinelNDWI.multiply(bNDWI))
.add(dem10.multiply(bDEM))
.add(slope10.multiply(bSlope))
.add(b0);

理论要点：回归系数在30米尺度上拟合，反映了LST与各辅助变量之间的统计关系。将同一组系数应用于10米尺度的辅助变量，本质上是假设这种统计关系在不同空间尺度上具有尺度不变性——这是统计降尺度方法的核心假设。

五、残差上采样与融合

（一）残差计算

30米尺度的回归残差定义为观测LST与拟合LST之差：

残差包含了回归模型未能解释的LST空间变异（如局地微气候效应、土壤湿度差异等）。

var residual30 = lst30.subtract(fittedLST30).rename('Residual_30m');

（二）高斯卷积与上采样

残差从30米上采样至10米的过程中，直接使用插值方法（如双三次插值）可能引入噪声。代码采用先卷积后插值的策略：

高斯核的定义如下：

var gaussianKernel = ee.Kernel.gaussian({
radius: 1.5,
sigma: 1,
units: 'pixels',
normalize: true
});

其中sigma=1控制平滑程度，radius=1.5为核半径（约覆盖3×3像元窗口）。

var residual10 = residual30
.unmask(0)
.convolve(gaussianKernel)
.resample('bicubic')
.reproject({ crs: sentinelProjection, scale: 10 })
.rename('Residual_10m');

（三）最终降尺度LST

将10米预测LST与10米残差相加，得到最终的10米降尺度LST：

var downscaledLST10 = predictedLST10
.add(residual10.unmask(0))
.updateMask(predictedLST10.mask())
.rename('LST_10m')
.clip(region);

残差融合技术的核心价值在于保证了降尺度结果在粗尺度上的总量一致性——即10米降尺度结果在30米聚合后的平均值与原始30米LST一致（或高度接近）。这避免了单纯依赖回归模型可能导致的系统性偏差。

六、精度评估指标

为量化回归拟合效果和降尺度结果质量，代码实现了一套完整的精度评估指标体系。

（一）回归拟合指标

基于30米尺度的观测LST（lst30）和拟合LST（fittedLST30），计算以下指标：

决定系数（R²） ：

均方根误差（RMSE） ：

平均绝对误差（MAE） ：

偏差（Bias） ：

function getMetrics(observed, predicted) {
// 计算SSE、SST、MSE、MAE、Bias等
// 返回包含R2、RMSE、MAE、Bias、N的字典
}

（二）降尺度前后统计对比

通过比较降尺度前后LST的统计特征（最小值、最大值、均值），可初步评估降尺度结果的合理性：

var beforeStatistics = lst30.reduceRegion({
reducer: ee.Reducer.minMax().combine(ee.Reducer.mean(), …),
geometry: region,
scale: 30
});

var afterStatistics = downscaledLST10.reduceRegion({
reducer: ee.Reducer.minMax().combine(ee.Reducer.mean(), …),
geometry: region,
scale: 10
});

七、结果导出与可视化

（一）地图可视化

代码提供降尺度前后的LST可视化对比：

var tempVis = { min: 20, max: 45, palette: ['blue', 'cyan', 'green', 'yellow', 'red'] };
Map.addLayer(lst30, tempVis, 'LST_30m (Before)');
Map.addLayer(downscaledLST10, tempVis, 'LST_10m (After)');

温度色阶的min/max需根据研究区实际LST范围调整——可通过print(lst30.reduceRegion({reducer: ee.Reducer.minMax()…}))获取。

（二）数据导出

最终结果通过Export.image.toDrive()导出为Cloud Optimized GeoTIFF：

Export.image.toDrive({
image: lst30.toFloat(),
description: 'LST_before_downscaling_30m',
folder: 'LST_downscaling',
scale: 30,
fileFormat: 'GeoTIFF',
formatOptions: { cloudOptimized: true }
});

建议同时导出降尺度前后的LST，便于后续在GIS或Python环境中进行对比分析。

八、总结

本文系统解析了一套基于GEE平台的Landsat-Sentinel-2 LST降尺度技术方案，核心要点可归纳如下：

方法框架：通过多元线性回归建立LST与高分辨率光谱指数（NDVI、NDBI、NDWI）及地形因子（DEM、坡度）之间的统计关系，结合残差融合技术实现从30米到10米的空间分辨率提升。

技术亮点：

• 利用Landsat C2L2产品的“开箱即用”LST波段，大幅降低预处理门槛；

• 引入清晰像元比例和影像评分机制，实现影像的自动优选；

• 高斯卷积+双三次插值的残差上采样策略，在保持空间细节的同时抑制噪声；

• 完整的精度评估体系（R²、RMSE、MAE、Bias），量化回归拟合效果。

适用场景：

• 城市热岛效应的精细尺度分析（如街区级热环境评估）；

• 生态环境监测（如自然保护区、湿地温度动态监测）；

• 农业干旱监测与作物热胁迫评估；

• 任何需要高分辨率LST数据但受限于热红外传感器空间分辨率的应用场景。

可复用经验：

• 多源遥感数据融合的“粗尺度建模 + 细尺度预测 + 残差校正”范式，可推广至其他地表参数（如土壤水分、蒸散发）的降尺度；

• GEE平台下ee.Reducer.linearRegression()配合reduceRegion()的回归建模流程，适用于各类像元级统计建模任务；

• 基于SCL波段的Sentinel-2像素级云掩膜策略，可复用于其他Sentinel-2分析流程。

本文展示的是基于多元线性回归的降尺度方法，其优势在于计算效率高、可解释性强。若研究区地表覆盖高度异质、LST与辅助变量之间存在非线性关系，可考虑将回归模型替换为随机森林或XGBoost等机器学习方法——GEE平台已原生支持这些算法，可实现类似的降尺度流程。

九、运行结果

原始的30m空间分辨率地表温度（LST）可视化结果

降尺度后的10m空间分辨率地表温度（LST）可视化结果

原始的30m空间分辨率地表温度（LST）局部细节展示（1）

降尺度后的10m空间分辨率地表温度（LST）局部细节展示（1）

原始的30m空间分辨率地表温度（LST）局部细节展示（2）

降尺度后的10m空间分辨率地表温度（LST）局部细节展示（2）

控制台输出的数据统计信息

点击RUN即可下载原始的30m空间分辨率和降尺度后的10m空间分辨率地表温度（LST）数据

若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！