Landsat数据在Google Earth Engine中的实战1986-2023年中国植被覆盖度分析指南植被覆盖度FVC作为衡量地表植被健康状况的核心指标其动态监测对生态研究、农业评估和气候变化分析具有不可替代的价值。Google Earth EngineGEE平台的出现让研究者无需下载海量遥感数据即可完成跨年度的植被覆盖分析。本文将手把手带您用JavaScript代码实现从Landsat数据预处理到年度FVC产品生成的全流程特别针对Landsat 5/7/8/9不同传感器的波段差异提供兼容性解决方案。1. 环境准备与数据加载在GEE代码编辑器https://code.earthengine.google.com/中新建脚本首先需要加载中国行政边界和研究时段。推荐使用GEE官方提供的FAO国家边界数据集var china ee.FeatureCollection(FAO/GAUL/2015/level0) .filter(ee.Filter.eq(ADM0_NAME, China)); var startYear 1986; var endYear 2023;对于Landsat数据GEE已预置了经过辐射校正和大气校正的地表反射率产品。我们需要创建按传感器分类的数据集集合// Landsat5 TM (1984-2012) var l5 ee.ImageCollection(LANDSAT/LT05/C02/T1_L2) .filterDate(startYear-01-01, 2012-12-31) .filterBounds(china); // Landsat7 ETM (1999-2023) var l7 ee.ImageCollection(LANDSAT/LE07/C02/T1_L2) .filterDate(1999-01-01, endYear-12-31) .filterBounds(china); // Landsat8 OLI (2013-2023) var l8 ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterDate(2013-01-01, endYear-12-31) .filterBounds(china);注意Landsat9数据目前仍在积累期如需使用可加载LANDSAT/LC09/C02/T1_L2集合其波段结构与Landsat8完全一致。2. 云掩膜与影像预处理不同Landsat卫星的QA波段结构略有差异需要分别构建云掩膜函数。这里展示一个支持多传感器的通用解决方案var cloudMask function(image) { // 通用QA波段解译 var qa image.select(QA_PIXEL); var cloudBitMask 1 3; var shadowBitMask 1 4; // Landsat7特有的条带修复 if (image.get(SPACECRAFT_ID) LANDSAT_7) { var dilatedCloud qa.bitwiseAnd(cloudBitMask).neq(0); var cloudShadow qa.bitwiseAnd(shadowBitMask).neq(0); var mask dilatedCloud.or(cloudShadow).not(); return image.updateMask(mask); } else { // Landsat5/8/9的标准处理 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)); return image.updateMask(mask); } };波段反射率需要按官方文档进行缩放处理var applyScaleFactors function(image) { var opticalBands image.select([SR_B.]).multiply(0.0000275).add(-0.2); return image.addBands(opticalBands, null, true); };3. 多传感器NDVI计算Landsat系列传感器的红光和近红外波段编号存在差异需要构建自适应计算函数var getNDVI function(image) { var spacecraft image.get(SPACECRAFT_ID); var nir, red; if (spacecraft LANDSAT_5) { nir SR_B4; red SR_B3; } else if (spacecraft LANDSAT_7) { nir SR_B4; red SR_B3; } else { // Landsat8/9 nir SR_B5; red SR_B4; } var ndvi image.normalizedDifference([nir, red]).rename(NDVI); return image.addBands(ndvi); };4. 年度最佳像元合成采用最佳NDVI策略从每年所有可用影像中选择代表值var createAnnualComposites function(collection, year) { var start ee.Date.fromYMD(year, 1, 1); var end start.advance(1, year); return collection .filterDate(start, end) .map(cloudMask) .map(applyScaleFactors) .map(getNDVI) .qualityMosaic(NDVI) // 选择NDVI最高的像元 .set(year, year); };5. 植被覆盖度FVC计算基于NDVI的混合像元分解模型计算FVCvar calculateFVC function(annualComposite, region) { var ndvi annualComposite.select(NDVI); // 动态计算NDVI极值5%和95%分位数 var percentiles ndvi.reduceRegion({ reducer: ee.Reducer.percentile([5, 95]), geometry: region, scale: 30, maxPixels: 1e13 }); var min ee.Number(percentiles.get(NDVI_p5)); var max ee.Number(percentiles.get(NDVI_p95)); // FVC计算公式 var fvc ndvi.subtract(min).divide(max.subtract(min)) .clamp(0, 1) // 限制在0-1范围内 .rename(FVC); return annualComposite.addBands(fvc); };6. 批量处理与结果导出构建完整处理流程并导出年度FVC结果// 生成年份序列 var years ee.List.sequence(startYear, endYear); // 合并多源Landsat数据 var mergedCollection l5.merge(l7).merge(l8); var annualFvc years.map(function(year) { year ee.Number(year).toInt(); var composite createAnnualComposites(mergedCollection, year); return calculateFVC(composite, china); }); // 转换为ImageCollection var fvcCollection ee.ImageCollection.fromImages(annualFvc); // 导出示例2020年FVC数据 Export.image.toDrive({ image: fvcCollection.filter(ee.Filter.eq(year, 2020)).first(), description: China_FVC_2020, scale: 30, region: china, maxPixels: 1e13, fileFormat: GeoTIFF, formatOptions: { cloudOptimized: true } });7. 可视化与验证在GEE中快速验证结果质量// 加载调色板 var fvcPalette [ FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400, 3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301 ]; // 创建时间序列动画 var animation ui.Thumbnail({ image: fvcCollection.select(FVC), params: { min: 0, max: 1, palette: fvcPalette, dimensions: 800x600, region: china, framesPerSecond: 5 } }); // 添加到地图 Map.addLayer(fvcCollection.filter(ee.Filter.eq(year, 2020)).first(), {min:0, max:1, palette: fvcPalette}, 2020 FVC); Map.centerObject(china, 4);8. 进阶技巧与优化建议多传感器一致性处理Landsat7 ETM的条带补偿使用focal_mean()进行局部平滑跨传感器归一化引入伪不变特征(PIF)校正计算效率优化// 使用指数移动窗口减少计算量 var optimizedReducer ee.Reducer.percentile({ percentiles: [5, 95], maxRaw: 1e5, outputNames: [p5, p95] }).unweighted().forEachBand([NDVI]);质量评估指标指标计算方法理想值云覆盖率QA波段统计10%NDVI动态范围(p95-p5)0.3有效像元比例非空像元计数80%常见问题排查结果中出现异常高值检查云掩膜是否完全生效验证NDVI极值计算区域是否合理空间不连续确认研究区域边界是否准确检查不同Landsat场景间的重叠处理年度数据缺失尝试放宽云量过滤阈值考虑使用相邻年份插值
Landsat数据在Google Earth Engine中的实战：1986-2023年中国植被覆盖度分析指南
发布时间：2026/5/22 18:28:12
Landsat数据在Google Earth Engine中的实战1986-2023年中国植被覆盖度分析指南植被覆盖度FVC作为衡量地表植被健康状况的核心指标其动态监测对生态研究、农业评估和气候变化分析具有不可替代的价值。Google Earth EngineGEE平台的出现让研究者无需下载海量遥感数据即可完成跨年度的植被覆盖分析。本文将手把手带您用JavaScript代码实现从Landsat数据预处理到年度FVC产品生成的全流程特别针对Landsat 5/7/8/9不同传感器的波段差异提供兼容性解决方案。1. 环境准备与数据加载在GEE代码编辑器https://code.earthengine.google.com/中新建脚本首先需要加载中国行政边界和研究时段。推荐使用GEE官方提供的FAO国家边界数据集var china ee.FeatureCollection(FAO/GAUL/2015/level0) .filter(ee.Filter.eq(ADM0_NAME, China)); var startYear 1986; var endYear 2023;对于Landsat数据GEE已预置了经过辐射校正和大气校正的地表反射率产品。我们需要创建按传感器分类的数据集集合// Landsat5 TM (1984-2012) var l5 ee.ImageCollection(LANDSAT/LT05/C02/T1_L2) .filterDate(startYear-01-01, 2012-12-31) .filterBounds(china); // Landsat7 ETM (1999-2023) var l7 ee.ImageCollection(LANDSAT/LE07/C02/T1_L2) .filterDate(1999-01-01, endYear-12-31) .filterBounds(china); // Landsat8 OLI (2013-2023) var l8 ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .filterDate(2013-01-01, endYear-12-31) .filterBounds(china);注意Landsat9数据目前仍在积累期如需使用可加载LANDSAT/LC09/C02/T1_L2集合其波段结构与Landsat8完全一致。2. 云掩膜与影像预处理不同Landsat卫星的QA波段结构略有差异需要分别构建云掩膜函数。这里展示一个支持多传感器的通用解决方案var cloudMask function(image) { // 通用QA波段解译 var qa image.select(QA_PIXEL); var cloudBitMask 1 3; var shadowBitMask 1 4; // Landsat7特有的条带修复 if (image.get(SPACECRAFT_ID) LANDSAT_7) { var dilatedCloud qa.bitwiseAnd(cloudBitMask).neq(0); var cloudShadow qa.bitwiseAnd(shadowBitMask).neq(0); var mask dilatedCloud.or(cloudShadow).not(); return image.updateMask(mask); } else { // Landsat5/8/9的标准处理 var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(shadowBitMask).eq(0)); return image.updateMask(mask); } };波段反射率需要按官方文档进行缩放处理var applyScaleFactors function(image) { var opticalBands image.select([SR_B.]).multiply(0.0000275).add(-0.2); return image.addBands(opticalBands, null, true); };3. 多传感器NDVI计算Landsat系列传感器的红光和近红外波段编号存在差异需要构建自适应计算函数var getNDVI function(image) { var spacecraft image.get(SPACECRAFT_ID); var nir, red; if (spacecraft LANDSAT_5) { nir SR_B4; red SR_B3; } else if (spacecraft LANDSAT_7) { nir SR_B4; red SR_B3; } else { // Landsat8/9 nir SR_B5; red SR_B4; } var ndvi image.normalizedDifference([nir, red]).rename(NDVI); return image.addBands(ndvi); };4. 年度最佳像元合成采用最佳NDVI策略从每年所有可用影像中选择代表值var createAnnualComposites function(collection, year) { var start ee.Date.fromYMD(year, 1, 1); var end start.advance(1, year); return collection .filterDate(start, end) .map(cloudMask) .map(applyScaleFactors) .map(getNDVI) .qualityMosaic(NDVI) // 选择NDVI最高的像元 .set(year, year); };5. 植被覆盖度FVC计算基于NDVI的混合像元分解模型计算FVCvar calculateFVC function(annualComposite, region) { var ndvi annualComposite.select(NDVI); // 动态计算NDVI极值5%和95%分位数 var percentiles ndvi.reduceRegion({ reducer: ee.Reducer.percentile([5, 95]), geometry: region, scale: 30, maxPixels: 1e13 }); var min ee.Number(percentiles.get(NDVI_p5)); var max ee.Number(percentiles.get(NDVI_p95)); // FVC计算公式 var fvc ndvi.subtract(min).divide(max.subtract(min)) .clamp(0, 1) // 限制在0-1范围内 .rename(FVC); return annualComposite.addBands(fvc); };6. 批量处理与结果导出构建完整处理流程并导出年度FVC结果// 生成年份序列 var years ee.List.sequence(startYear, endYear); // 合并多源Landsat数据 var mergedCollection l5.merge(l7).merge(l8); var annualFvc years.map(function(year) { year ee.Number(year).toInt(); var composite createAnnualComposites(mergedCollection, year); return calculateFVC(composite, china); }); // 转换为ImageCollection var fvcCollection ee.ImageCollection.fromImages(annualFvc); // 导出示例2020年FVC数据 Export.image.toDrive({ image: fvcCollection.filter(ee.Filter.eq(year, 2020)).first(), description: China_FVC_2020, scale: 30, region: china, maxPixels: 1e13, fileFormat: GeoTIFF, formatOptions: { cloudOptimized: true } });7. 可视化与验证在GEE中快速验证结果质量// 加载调色板 var fvcPalette [ FFFFFF, CE7E45, DF923D, F1B555, FCD163, 99B718, 74A901, 66A000, 529400, 3E8601, 207401, 056201, 004C00, 023B01, 012E01, 011D01, 011301 ]; // 创建时间序列动画 var animation ui.Thumbnail({ image: fvcCollection.select(FVC), params: { min: 0, max: 1, palette: fvcPalette, dimensions: 800x600, region: china, framesPerSecond: 5 } }); // 添加到地图 Map.addLayer(fvcCollection.filter(ee.Filter.eq(year, 2020)).first(), {min:0, max:1, palette: fvcPalette}, 2020 FVC); Map.centerObject(china, 4);8. 进阶技巧与优化建议多传感器一致性处理Landsat7 ETM的条带补偿使用focal_mean()进行局部平滑跨传感器归一化引入伪不变特征(PIF)校正计算效率优化// 使用指数移动窗口减少计算量 var optimizedReducer ee.Reducer.percentile({ percentiles: [5, 95], maxRaw: 1e5, outputNames: [p5, p95] }).unweighted().forEachBand([NDVI]);质量评估指标指标计算方法理想值云覆盖率QA波段统计10%NDVI动态范围(p95-p5)0.3有效像元比例非空像元计数80%常见问题排查结果中出现异常高值检查云掩膜是否完全生效验证NDVI极值计算区域是否合理空间不连续确认研究区域边界是否准确检查不同Landsat场景间的重叠处理年度数据缺失尝试放宽云量过滤阈值考虑使用相邻年份插值
