GEE批量处理Landsat C02 1985-2024植被指数归一化实战

简介面向需要基于 Google Earth Engine 长时间序列遥感指数计算的科研人员与学习者这份教程以 Landsat C02 数据集为基准系统讲解 1985—2024 年 NDVI、EVI、SAVI、NDMI 等指数的归一化处理方法。教程针对 Landsat 5/7/8 分别整合数据优化了去云与预处理流程并统一重命名原始波段显著简化不同传感器之间指数计算的复杂度使读者可以直接调用任意时期影像完成批量归一化。资源为 PDF 文档共 1 个文件大小约 900 KB便于离线阅读与对照实操。文档从 QA_PIXEL 云掩膜、辐射定标缩放到波段重命名和 NDVI、NDWI、NDBI、EVI 等函数封装均有完整代码示例可帮助理解时序影像处理的关键环节。目前已有 977 人学习下载适合具备一定 GEE 基础、希望提升长时序植被指数处理效率的读者。1. 为什么2025年还要回头处理Landsat C02的1985—2024指数归一化如果你现在打开一个2021年之前写的GEE教程复制NDVI计算代码大概率会得到一堆数值在99左右浮动的鬼畜结果——这不是算法错了而是数据版本变了。Landsat Archive在2022年全面切到Collection 2后表面反射率产品从T1_SR变成了T1_L2波段名从B4/B5变成了SR_B4/SR_B5SR波段的自带Scale因子也从0.0001改成了0.0000275旧代码里凡是手写常量的地方全部失效。这篇文章要解决的就是用GEE在线把这套旧流程整体重写一遍从Landsat 5/7/8/9的C02 Level 2数据出发批量拉取1985—2024共40年年影像统一波段名算NDVI、EVI、SAVI、NDMI再做长时序逐像素归一化顺带给出FVC和MAD检验的进阶用法。全程在GEE Code Editor里跑不需要本地下载任何影像。适合被旧教程带偏、或者刚接手长时序植被指数任务的人直接照着复现。2. 从C02 Level 2批量拉取Landsat 1985—2024影像数据集选型与统一波段2.1 为什么只选Landsat C02而不是Collection 1或Sentinel-2Selection理由很直接1985年这个起点把Sentinel-2直接排除在外——哨兵2号2015年才发射覆盖不了前30年。MOD13Q1虽然从2000年开始且时间连续但它是250米分辨率的植被指数产品不是原始反射率后续没法自己算EVI/SAVI/NDMI。Collection 1T1_SR时间跨度倒是够但USGS在2022年后停止更新而且没有QA_PIXEL质量波段做云掩膜要靠CFMask属性猜精度差一截。C02 Level 2数据集ID后缀T1_L2是目前唯一能串起Landsat 5/7/8/9四颗卫星、原生30米分辨率、自带质量波段的表面反射率产品。它的关键变化有三点波段名带SR_前缀SR波段Scale因子是0.0000275新增QA_PIXEL位标志波段。这三个变化直接决定了旧教程为什么集体翻车也决定了下文所有代码的写法。对比项Collection 1 SRCollection 2 Level 2Sentinel-2 SR时间跨度1982—20211984—至今2015—至今波段命名B2/B3/B4…SR_B1/SR_B2/SR_B3…B2/B3/B4…SR Scale因子0.00010.00002750.0001QA波段CFMask属性QA_PIXELSCL/QA60适合长时序L5—L8可拼L5—L9可拼质量文件完整仅近10年2.2 统一四颗卫星的波段名L5/L7与L8/L9的错位是第一个大坑C02的SR_波段名在Landsat 5/7和8/9之间并不对齐。Landsat 5/7没有海岸气溶胶波段所以SR_B1是蓝波段、SR_B3是红、SR_B4是近红外而Landsat 8/9的SR_B1是海岸气溶胶、SR_B2才是蓝、SR_B4是红、SR_B5是近红外。如果所有卫星共用一套SR_B4当红波段L5/7算出来的是近红外NDVI公式直接分子分母用错。我一般会先把四颗卫星统一重映射成blue/red/nir/swir1/swir2五个逻辑波段后面所有指数公式只跟这五个名字打交道。以下是Landsat 5/7的预处理函数function prepL57(img) { var qa img.select(QA_PIXEL); var mask qa.bitwiseAnd(1).eq(0) // 非填充L7条带用这一位剔除 .and(qa.bitwiseAnd(8).eq(0)) // 非云 .and(qa.bitwiseAnd(16).eq(0)); // 非云阴影 return img .updateMask(mask) .select( [SR_B1, SR_B3, SR_B4, SR_B5, SR_B7], [blue, red, nir, swir1, swir2] ) .multiply(0.0000275) .copyProperties(img, [system:time_start]); }这段代码里做了三件事先用QA_PIXEL做掩膜再把波段序号映射成逻辑名最后乘Scale因子。bitwiseAnd(1)对应QA_PIXEL的第0位填充标志Landsat 7在2003年SLC传感器故障后条带空洞区域会被标记为填充这一位能同时解决条带问题。乘0.0000275是C02 L2的官方SR Scale因子注意是零点零零零零二七五不是C01的0.0001。Landsat 8/9的波段映射不一样函数单独写function prepL89(img) { var qa img.select(QA_PIXEL); var mask qa.bitwiseAnd(1).eq(0) .and(qa.bitwiseAnd(8).eq(0)) .and(qa.bitwiseAnd(16).eq(0)); return img .updateMask(mask) .select( [SR_B2, SR_B4, SR_B5, SR_B6, SR_B7], [blue, red, nir, swir1, swir2] ) .multiply(0.0000275) .copyProperties(img, [system:time_start]); }Landsat 8/9里SR_B1是海岸气溶胶不能当蓝波段用映射时直接从SR_B2开始。两个函数处理后输出波段名完全一致后续合并集合、算指数、归一化都不用再区分传感器。2.3 QA_PIXEL掩膜与合并一条语句带走云、阴影和L7条带QA_PIXEL是一个位标志波段每位代表一种质量状态第0位是填充、第1位是膨胀云、第2位是卷云、第3位是云、第4位是云阴影。按位bitwiseAnd做掩膜比C01里靠CLOUD_COVER属性过滤要精细得多因为它能作用到像元级而不是整景图一刀切。合并四颗卫星并构建全时序集合的代码var roi ee.FeatureCollection(users/your_username/roi).geometry(); var startDate 1985-01-01; var endDate 2024-12-31; var l7 ee.ImageCollection(LANDSAT/LE07/C02/T1_L2) .filterBounds(roi) .filterDate(startDate, endDate) .map(prepL57); var l89 ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) .merge(ee.ImageCollection(LANDSAT/LC09/C02/T1_L2)) .filterBounds(roi) .filterDate(startDate, endDate) .map(prepL89); var allImages l7.merge(l89); print(总景数, allImages.size());Landsat 5的数据在2013年就停了1985到1999年这段只有L5能覆盖所以还要单独把LANDSAT/LT05/C02/T1_L2加进集合。Landsat 5需要传到结束年份L7从1999年接力L8从2013年接力L9从2021年底接力四段拼起来正好覆盖1985到2024。合并前先filterBounds缩到ROI能少拉大量无关景导出和reduce都轻不少。3. NDVI/EVI/SAVI/NDMI统一计算与年合成公式、波段依赖与合成策略3.1 四个指数的公式、波段依赖和适用场景这四个指数的核心差异在于各自压制了不同的干扰因素。NDVI对土壤背景敏感EVI引入蓝波段压制大气气溶胶并降低土壤背景SAVI用土壤调节系数L来解决植被稀疏区的地表背景问题NDMI用的是近红外与短波红外用来监测植被含水量和干旱胁迫。指数公式依赖波段典型场景NDVI(NIR - Red) / (NIR Red)red, nir植被覆盖度、物候EVI2.5*(NIR - Red) / (NIR 6Red - 7.5Blue 1)blue, red, nir高植被区、大气校正后时序SAVI(NIR - Red) / (NIR Red L) * (1L)red, nir干旱半干旱区、稀疏植被NDMI(NIR - SWIR1) / (NIR SWIR1)nir, swir1植被含水量、干旱监测EVI里- 7.5*Blue这个系数不是玄学它是为了在蓝波段上抵消大气路径辐射的影响。SAVI的L系数在GEE社区最常见的取法是0.5对应中低植被密度如果你做的是极端荒漠区L取1更稳植被茂密区取0.25。NDMI的SWIR1在Landsat上对应1.55—1.65微米区间对水分极其敏感干旱年份这个值的下降速度比NDVI快得多。这套统一的blue/red/nir/swir1波段集合还能顺手扩展其他指数NDWI用green和nirNDBI用swir1和nirMarginal指数用swir1和red。标题里那个等指数本质就是同一个波段集合换组合公式的问题。3.2 逐景指数计算与年最大值合成指数计算直接在allImages集合上map输出多波段影像每一景都带system:time_start时间戳function calcIndices(img) { var ndvi img.normalizedDifference([nir, red]).rename(NDVI); var evi img.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: img.select(nir), RED: img.select(red), BLUE: img.select(blue) }).rename(EVI); var savi img.expression( (1 L) * (NIR - RED) / (NIR RED L), { NIR: img.select(nir), RED: img.select(red), L: 0.5 }).rename(SAVI); var ndmi img.normalizedDifference([nir, swir1]).rename(NDMI); return img.addBands([ndvi, evi, savi, ndmi]) .select([NDVI, EVI, SAVI, NDMI]) .copyProperties(img, [system:time_start]); } var indexed allImages.map(calcIndices);expression里的L: 0.5是作为常数传入的GEE的expression函数会把字典里的值当成数值常量参与运算这种方式比写死在公式字符串里更清晰。normalizedDifference([nir,red])就是(nir - red)/(nir red)省去手动写减法除法。年合成我倾向于用最大值合成而不是平均值合成。原因有两个一是最大值合成天然把云的残留影响压到最低——云区的指数通常显著低于真实地表取最大值时会被干净像元顶掉二是Landsat重访周期16天一年内生长期峰值本来就该用最大值来捕捉。按年合成代码如下var years ee.List.sequence(1985, 2024); var annual years.map(function(y) { var yearStr ee.Number(y).format(%04d); var yearImages indexed.filter(ee.Filter.calendarRange(y, y, year)); var maxImg yearImages.reduce(ee.Reducer.max()) .set(year, y) .set(system:time_start, ee.Date.fromYMD(y, 1, 1)); return maxImg; }); var annualCol ee.ImageCollection(annual);ee.Filter.calendarRange(y, y, year)按年过滤比字符串拼接日期稳得多尤其跨1月1日边界时不会漏景。reduce(ee.Reducer.max())输出波段名会带_max后缀比如NDVI_max后面归一化时select(NDVI_max)要对应上。这一步得到的是每个像元在当年所有有效观测里的最大指数值。3.3 按年导出到DriveExport的scale、crs与maxPixels参数导出前先做一步极低成本检查——打印一景缩略图确认数值范围正常这能避免导出半天后发现波段名写错// 缩略图预览 var testYear annualCol.filterMetadata(year, equals, 2020).first(); Map.centerObject(roi, 8); Map.addLayer(testYear.select([NDVI_max]).clip(roi), {min: 0, max: 1, palette: [brown, yellow, green]}, NDVI 2020);确认没问题后再批量导出。GEE导出是按任务走的一个Export对应一个任务40年就是40个任务可以在Tasks面板里批量提交annualCol.toList(annualCol.size()).forEach(function(img) { var imgEE ee.Image(img); var y imgEE.get(year); Export.image.toDrive({ image: imgEE.select([NDVI_max, EVI_max, SAVI_max, NDMI_max]).clip(roi), description: Landsat_C02_Index_ y, folder: GEE_Index, region: roi, scale: 30, crs: EPSG:32650, maxPixels: 1e13 }); });这里几个参数我踩过坑scale必须是30不要让GEE默认用影像原始分辨率干瞪眼crs建议直接用当地UTM分带比如中国东部用EPSG:32650用EPSG:4326导出会让栅格在每个纬度上的实际地面分辨率变形maxPixels设到1e13是因为ROI大时30米分辨率像元数很容易超过默认的1e8上限。四个指数一次性导出成四波段GeoTIFF比四个指数分别导出四次省时间。4. 长时序归一化min-max与分位数两种方案的参数设定4.1 归一化的数学含义与把归一化指数变成FVC的链路很多人以为归一化只是把数值压到0到1之间其实在植被指数语境里它做的是把绝对指数值变成相对于该像元历史区间的位置。同一个像元2000年NDVI是0.352023年NDVI是0.40绝对值只能说后者略高但如果这个像元40年NDVI范围是0.15到0.50归一化后0.40对应0.71的位置2000年只对应0.57这个差异才真正反映植被状态的变化。更重要的是归一化后的数值直接就是植被覆盖度FVC的近似FVC (NDVI - NDVI_min) / (NDVI_max - NDVI_min)。这就是这个标题里归一化最值得做的原因——它不只是数据预处理它把一个无量纲指数变成了一个有生态学含义的物理量。做法上最小值NDVI_min取裸土或水体在时序上的低值最大值NDVI_max取茂密植被的高值逐像元统计比用全球固定值精确得多。4.2 40年极值归一化min-max逐像素统计与分母保护先对40幅年合成影像逐像素求最小值、最大值然后逐年做归一化var ndviBands annualCol.select([NDVI_max]); var minImg ndviBands.reduce(ee.Reducer.min()); var maxImg ndviBands.reduce(ee.Reducer.max()); var ndviNorm annualCol.map(function(img) { var ndvi img.select(NDVI_max); var norm ndvi.subtract(minImg).divide(maxImg.subtract(minImg)) .rename(NDVI_norm); return norm.set(system:time_start, img.get(system:time_start)); });分母maxImg.subtract(minImg)就是NDVI的动态范围。问题在于这个分母在部分像元上会接近0——比如常年被水体覆盖的像元NDVI全年恒定在0附近分母趋近0归一化结果会变成极大值或NaN。我在做这个步骤时习惯给分母做一个下限保护var denom maxImg.subtract(minImg).max(0.01); var ndviNormFix annualCol.map(function(img) { var norm img.select(NDVI_max) .subtract(minImg) .divide(denom) .rename(NDVI_norm); return norm; });max(0.01)意思是分母最小取0.01当动态范围小于0.01时按0.01算。0.01不是拍脑袋定的NDVI本身理论范围就是-1到1动态范围0.01以下的像元要么是水体要么是永久裸岩归一化结果本身没有生态意义用下限保护防止导出数据里出现无穷值。如果ROI里水体占比很大可以把下限调到0.05但对常规植被研究0.01够用。4.3 分位数归一化5%—95%用median性抗异常min-max最大的弱点是扛不住异常值。40年里只要有几年Landsat 7条带没掩干净、或者某年大火把地表烧成极端低值全局最小值和最大值就会被这几景污染归一化结果整体被压缩。分位数归一化用5%分位替代最小、95%分位替代最大这两端的异常年份只影响各自5%的像元不会把全时序的范围都拉偏。var p5 ndviBands.reduce(ee.Reducer.percentile([5])); var p95 ndviBands.reduce(ee.Reducer.percentile([95])); var denomPct p95.subtract(p5).max(0.01); var ndviNormPct annualCol.map(function(img) { var norm img.select(NDVI_max) .subtract(p5) .divide(denomPct) .rename(NDVI_norm_pct); return norm; });分位数参数的选择湿润区、森林为主的研究区用5%—95%就行干旱半干旱区植被稀疏、NDVI整体偏低且年际波动大我会放宽到2%—98%否则稀疏植被的归一化结果会被有植被年份压缩在低值区间。ee.Reducer.percentile需要传入一个数字列表写[5]和[95]分别reduce两次比一次传[5, 95]得到双波段影像再拆更容易复用。用哪种方案取决于你的下游任务做FVC制图、物候分析min-max配合下限保护就够做长时序变化检测、干旱监测分位数归一化更抗噪。我建议两个都算出来一个用于生态参数标定一个用于变化检测输入。5. GEE处理Landsat C02的5个常见陷阱现象、原因与解决5.1 陷阱NDVI结果大面积出现99或负指数现象Map显示NDVI最大值超过99或者NDVI直方图集中在-30000附近。原因没有乘SR Scale因子直接在DN值上算(B4-B3)/(B4B3)。C02 L2的SR波段存储的是整数范围在0到32767之间真实反射率需要乘以0.0000275。DN值上千时比值会放大上百倍。解决在预处理函数里统一乘0.0000275不要在指数公式里乘。乘一次后续所有指数、合成、归一化都不用再管Scale问题。5.2 陷阱Landsat 5波段名在C02里跟旧教程对不上现象从LANDSAT/LT05/C02/T1_L2导出的影像select(B4)直接报错说波段不存在。原因C02把波段名从B1/B2/B3/B4改成了SR_B1/SR_B2/SR_B3/SR_B4旧教程里所有select都要改。更隐蔽的是SR_B4在L5/7上是近红外但在L8/9上是红波段两代卫星的波段序号错位。解决先print(image.bandNames())确认实际波段名再做统一波段映射。不要相信任何旧教程里写死的波段顺序。5.3 陷阱EVI计算结果整体偏低或出现负值现象同一研究区EVI均值比NDVI低了一半以上或者大范围负值。原因蓝波段选错。Landsat 8/9的SR_B1是海岸气溶胶波段不是蓝波段如果沿用Landsat 5时代的映射把SR_B1当蓝波段代入EVI公式-7.5*Blue这项会引入错误的气溶胶信号把EVI整体压低。解决L5/7的蓝波段是SR_B1L8/9的蓝波段是SR_B2在预处理阶段分别映射成统一的blueEVI公式只认blue不要再按传感器分开写。5.4 陷阱年合成影像布满细条纹现象2003年之后的L7影像年合成结果出现规则斜向条纹NDVI在条纹处断崖式降低。原因Landsat 7 SLC传感器2003年故障后影像边缘出现数据条带。C02的QA_PIXEL里这些条带区域被标记为填充bit 0如果掩膜只处理云和阴影没有处理填充位条带区域的DN值会以异常低值混入合成计算。解决QA_PIXEL掩膜里必须加上bitwiseAnd(1).eq(0)把填充像元全部剔除。我见过很多人只用云掩膜然后怪L7数据质量差其实一个bit位就能解决。5.5 陷阱Export导出任务跑几个小时然后报错现象Tasks面板里任务长期Running随后显示Error错误信息是Computed image is too large或User memory limit exceeded。原因直接导出整个annualCol集合或者ROI范围没裁剪GEE计算图里所有年份的像元数叠加超出内存限制。解决三个动作导出前clip(roi)每个Export只导出一年的影像不要一次导全部年份maxPixels显式设到1e13。如果还是超限把scale改到60先导出一次验证数据流确认没问题再回30做正式处理。6. 进阶把归一化结果做成FVC并用MAD做稳健性检验6.1 FVC落地直接用归一化NDVI剪裁FVC的公式和NDVI归一化公式完全一致所以上面算的NDVI_norm本身就是FVC的估计值只需要把超出0到1范围的值裁剪掉var fvc annualCol.map(function(img) { return img.select(NDVI_norm) .clamp(0, 1) .rename(FVC) .set(system:time_start, img.get(system:time_start)); });.clamp(0, 1)把低于0的像元置0、高于1的置1。这一步不能省因为min-max归一化在某些像元上会产出负值或大于1的值尤其是水体等非植被地类。6.2 用MAD检验归一化结果是否被异常年污染中位数绝对偏差MAD是比标准差更抗异常值的离散度指标。归一化后我一般会算每个像元40年时序的MADMAD过大的位置大概率是掩膜没掩干净或者地表发生真实突变var fvcCol fvc.select(FVC); var fvcMedian fvcCol.reduce(ee.Reducer.median()); var fvcMad fvcCol.map(function(img) { return img.subtract(fvcMedian).abs(); }).reduce(ee.Reducer.median()).rename(FVC_MAD); Map.addLayer(fvcMad, {min: 0, max: 0.3, palette: [white, red]}, FVC MAD);MAD超过0.2的像元建议抽查原始影像确认是云残留还是真实植被变化。这一步相当于给归一化结果做体检能提前发现5.4那种条带残留比跑到一半再开排查省心得多。实测经验是经过QA_PIXEL完整掩膜的数据常规植被区MAD大多在0.05以下超过0.15的像元基本都有问题。6.3 验证方法同一年份归一化前后叠加比对最后说一个我自己的验证习惯把2020年的原始NDVI和归一化后的FVC叠加在两个图层上打开混合模式目视检查。原始NDVI高值区0.6在FVC图里应接近1低值区0.2应接近0如果出现大面积错配说明min或max统计有误回头检查是不是合并时把未经掩膜的L7数据混了进来。这套流程跑完1985—2024年逐年FVC可以直接作为后续趋势分析和变化检测的输入。希望这些经验能帮你少走我走过的弯路。本文还有配套的精品资源点击获取