简介：面向遥感与地理信息分析学习者，提供基于Google Earth Engine（GEE）的Landsat系列指数归一化处理教程，涵盖1985—2024年NDVI、EVI、SAVI、NDMI等常用植被与水分指数的计算流程。教程针对Landsat 5/7/8集合分别整合原始波段并统一波段命名，优化了云掩膜与预处理环节，使长时序影像的批量归一化操作更简洁高效。资源为单份PDF文档，共1个文件，大小约900KB，适合需要系统掌握GEE指数计算与数据预处理方法的入门及进阶用户。目前已有977人学习下载。通过该PDF可获取去云函数、波段缩放与重命名、指数计算等核心代码片段，并了解如何将不同Landsat传感器的数据纳入同一处理框架，便于后续开展植被动态监测、生态环境评估等长时序分析。

1. 用 GEE 把 1985—2024 年 Landsat 指数一次性归一化：先搞懂 C02 再动手

遥感数据处理里最烦的一件事，就是不同传感器之间的指数没法直接比较。Landsat 5、7、8 的波段设置不一样，大气校正后的数据格式也一直在变，加上云遮挡和条带丢失，你从 USGS 拉下来的 SR 数据如果直接算 NDVI、EVI，不同年份之间数值漂移很严重。这份 GEE 教程的价值在于：它把 Landsat 5/7/8 的 C02 集合统一做了去云、缩放、波段重命名，再通过一个 band-wise 的 mean±3σ 归一化函数，把 1985 到 2024 年任意时相的影像压到 0—1 区间，方便后续做变化检测、分类或者时序分析。适合正在做长时序生态环境监测、植被覆盖度估算，或者被多传感器数据整合搞得头疼的研究生和工程师。我拆完这份代码后发现，它解决的不只是「归一化」这个动作，更关键的是把「不同卫星的数据怎么对齐」这件事给了个可复现的模板。

2. 为什么必须处理 C02 数据：从缩放因子到 QA 波段

2.1 C02 和 C01 在数值上的本质差异

如果你之前用的是 COLE/RADIOMETRIC 那套老数据，切到 C02 后第一个要改的就是缩放方式。C02 的 SR 波段已经是反射率乘以 10000 的整数存储，而 C01 的 surface reflectance 需要乘以 0.0001 才能得到 0—1 的反射率。这份代码里用的是：

var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);

为什么不是 0.0001？因为 C02 Collection 2 的 Surface Reflectance 产品实际采用的缩放系数是 0.0000275，偏移量是 -0.2。也就是说，DN 值乘以 0.0000275 再减去 0.2，才能还原成真实反射率。这个参数如果你还用 C01 的习惯，算出来的 NDVI 整体偏低，因为红光和近红外的值域被压窄了。

热红外波段更需要注意，Landsat 5/7 的 ST_B6 用的是 0.00341802 和 149.0 的缩放，而 Landsat 8 的 ST_B10、ST_B11 同样是这两个系数，但波段名不同。代码里对 maskL8sr 用的是image.select('ST_B.*')，也就是用正则把 ST_B10 和 ST_B11 都选出来，再统一乘系数。这个细节很关键，因为 Landsat 8 有两个热红外波段，你如果只选了 ST_B10，地表温度产品就少了一个维度。

2.2 QA_PIXEL 位掩码的真正含义

C02 产品的质量评估波段从 C01 的 QA_BAND 换成了 QA_PIXEL，而且 C02 的 QA 波段用的不是简单的整数标记，而是位掩码。代码里这样写：

var qaMask = image.select('QA_PIXEL').bitwiseAnd(parseInt('11111', 2)).eq(0);

parseInt('11111', 2)是二进制 31，也就是把低 5 位全部置 1。这 5 位分别对应 Fill、Dilated Cloud、Cirrus、Cloud、Cloud Shadow。bitwiseAnd的结果如果等于 0，说明这五个条件都不满足，也就是这个像元既不是填充值、也不是云、也不是云影，可以放心使用。

这里有一个隐藏的坑：Landsat 5/7 的 QA_PIXEL 没有 Cirrus 位（因为 TM 和 ETM+ 没有卷云波段），所以代码里 maskL457sr 的注释中 Bit 2 写的是 Unused，而 maskL8sr 的 Bit 2 是 Cirrus。但两个函数用的parseInt('11111', 2)是一样的，也就是低 5 位都做了掩膜。对于 Landsat 5/7 来说，Bit 2 永远是 0，不影响结果；对于 Landsat 8 来说，Bit 2 的 Cirrus 检测如果被误判，会把一些薄云区域的像元剔除掉。如果你发现某些区域的有效像元比预期少，可以检查一下 Cirrus 位的影响。

2.3 为什么要把 QA_RADSAT 也用上

var saturationMask = image.select('QA_RADSAT').eq(0);

QA_RADSAT 是辐射饱和度标记，如果某个波段的 DN 值达到传感器饱和，这个像元在这一个波段上就是无效的。在极端亮目标（雪、云、盐碱地）上，饱和像元在 C02 产品里仍然会给出一个数值，但那个数值已经不可信了。eq(0)表示所有波段都没有饱和。这里要注意，QA_RADSAT 有 8 个位，分别对应该影像的 8 个波段，如果只是局部波段饱和，整个像元会被剔除，这在有些场景下会丢掉不少边缘像元，但为了时间序列的一致性，这个代价是值得的。

3. 把 Landsat 5/7/8 对齐到同一套波段名：rename 函数与集合整合

3.1 波段名不一致是所有多源遥感分析的第一个坎

Landsat 5/7 的 SR 波段叫 SR_B1、SR_B2、SR_B3、SR_B4、SR_B5、SR_B7（注意没有 SR_B6，因为那个位置是热红外），而 Landsat 8 的 SR 波段是 SR_B2、SR_B3、SR_B4、SR_B5、SR_B6、SR_B7（因为多了海岸气溶胶波段 SR_B1）。如果你不统一命名，后面写指数函数时就得为每颗卫星写一套，代码冗余不说，还容易张冠李戴。这份教程的做法是把它们分别重命名为 blue、green、red、nir、swir1、swir2：

function rename57 (image) { return image.select(['SR_B1','SR_B2','SR_B3','SR_B4',"SR_B5","SR_B7"], ['blue', 'green', 'red','nir',"swir1","swir2"]); } function rename89 (image) { return image.select(['SR_B2','SR_B3','SR_B4',"SR_B5","SR_B6","SR_B7"], ['blue', 'green', 'red','nir',"swir1","swir2"]); }

这两段代码的逻辑是：先按位置选波段，再给新名字。select的第一个数组是原波段名，第二个数组是新波段名，一一对应。注意 Landsat 5/7 选的是 SR_B1 到 SR_B7（跳过 SR_B6），Landsat 8 选的是 SR_B2 到 SR_B7（跳过 SR_B1，因为那是海岸波段）。重命名之后，Landsat 8 的 SR_B6（短波红外 1）会被命名为 swir1，SR_B7 命名为 swir2，这样就和 Landsat 5/7 的 SR_B5、SR_B7 对齐了。

一个常见的疑惑是：为什么不用image.select(['SR_B.*'])把可见光和近红外波段全选出来？因为那样会把 Landsat 8 的海岸波段也包含进来，而且波段顺序是 B1、B2、B3、B4、B5、B6、B7，对应到 blue、green、red 时会错位。显式列出波段名虽然啰嗦，但最可靠。

3.2 三个集合构建与 merge 的细节

var collection1 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2') .filterBounds(geometry) .filterDate('1985-01-01', '2012-01-01') .map(maskL457sr).map(rename57) .map(NDVI).map(NDWI).map(NDBI).map(EVI); var collection2 = ee.ImageCollection('LANDSAT/LE07/C02/T1_L2') .filterBounds(geometry) .filterDate('1985-01-01', '2012-01-01') .map(maskL457sr).map(rename57) .map(NDVI).map(NDWI).map(NDBI).map(EVI); var collection3 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterBounds(geometry) .filterDate('2013-01-01', '2024-01-01') .map(maskL8sr).map(rename89) .map(NDVI).map(NDWI).map(NDBI).map(EVI); var col = collection1.merge(collection2).merge(collection3);

这里有几层逻辑值得拆开说。第一，Landsat 5 和 7 的日期范围都写了 1985—2012，这意味着 LT05 在 2012 年之后其实已经没有数据了，LE07 在 2003 年之后出现了 SLC 故障（条带丢失），但代码没有对 Landsat 7 做额外的条带掩膜。如果你直接把 2003 年之后的 Landsat 7 影像纳入时序，条带区域会在归一化后出现明显的空洞。实际上 ETM+ 的 SLC-off 影像在 GEE 里可以通过image.updateMask(image.select('SR_B3').mask().and(image.select('SR_B4').mask()))之类的方式做处理，但这份代码没有做，你如果要用 2003 年后的 Landsat 7，需要自己补一步。

第二，merge 的顺序是 collection1、collection2、collection3，也就是 Landsat 5、7、8 拼接。优先级上，如果三个集合里同时存在同一时相的影像（比如 2012 年 Landsat 5 和 Landsat 7 都有数据），merge 会保留先出现的那个，也就是 Landsat 5 优先。这个顺序不是绝对的，你可以按质量或时间偏好调整。

第三，每个集合的map顺序是先去云、再重命名、再算指数。这一步一个常见失误是：如果你把 rename 放在 mask 之前，QA_PIXEL 波段还在，但 SR_B 波段已经被重命名，后面 maskL457sr 里image.select('SR_B.')就会失效。这份代码把 mask 放在 rename 前，顺序是对的。我自己写的时候习惯把 mask 和 rename 合成一个函数，减少 map 次数，但分开写的好处是方便单独调试每一步。

3.3 指数函数的两种写法对比

代码里给了两组指数计算的写法：NDVI、NDWI、NDBI 用的是normalizedDifference，EVI 用的是expression：

function NDVI(image) { return image.addBands( image.normalizedDifference(["nir", "red"]).rename("NDVI")); } function EVI(image) { var evi = image.expression( '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', { 'NIR': image.select('nir'), 'RED': image.select('red'), 'BLUE': image.select('blue') }); return image.addBands(evi.rename('EVI')); }

normalizedDifference(["nir", "red"])计算的是(nir - red) / (nir + red)，这是 NDVI 的标准公式。NDWI 在这份代码里用的是(red - swir1) / (red + swir1)，这是 Gao 提出的基于短波红外的水体指数，和 McFeeters 的 NDWI（绿波段和近红外）是两码事。你如果是做水体提取，要确认自己用的是哪个版本；这份代码的 NDWI 公式对植被含水量更敏感，对开阔水体的响应反而一般。

EVI 的expression写法把公式直接写在字符串里，参数通过字典传入，好处是公式一目了然，坏处是expression的执行效率比normalizedDifference低一些。在一个长时序的 ImageCollection 上，如果每一景都要跑一遍 expression，计算量会明显增大。一个优化方式是把 EVI 公式里的常数提出来用乘法和减法表达，但可读性会下降。对于 1985 到 2024 年的长时序，我建议先把集合按年份合成（比如用median()），再做指数计算，效率会高很多。

4. band-wise 归一化：mean±3σ 的截断逻辑与 reduceRegion 参数

4.1 为什么用 mean±3σ 而不是 min-max

代码里的归一化函数是这份教程最核心的部分，它没有用传统的 min-max 拉伸，而是基于每个波段的均值和标准差做截断：

function normalization(image, region, scale) { var mean_std = image.reduceRegion({ reducer: ee.Reducer.mean() .combine(ee.Reducer.stdDev(), null, true), geometry: region, scale: scale, maxPixels: 10e9 }); var unitScale = ee.ImageCollection.fromImages( image.bandNames().map(function(name) { name = ee.String(name); var band = image.select(name); var mean = ee.Number(mean_std.get(name.cat('_mean'))); var std = ee.Number(mean_std.get(name.cat('_stdDev'))); var max = mean.add(std.multiply(3)); var min = mean.subtract(std.multiply(3)); var band1 = ee.Image(min).multiply(band.lt(min)) .add(ee.Image(max).multiply(band.gt(max))) .add(band.multiply(ee.Image(1) .subtract(band.lt(min)).subtract(band.gt(max)))); var result_band = band1.subtract(min).divide(max.subtract(min)); return result_band; })).toBands().rename(image.bandNames()); return unitScale; }

这段代码的骨架逻辑是：先对影像的每个波段算均值和标准差，然后用mean ± 3 * std作为上下界，把超出这个范围的像元值截断到边界，最后用(value - min) / (max - min)线性拉伸到 0—1。

为什么用 3σ 而不是 min-max？因为遥感影像里如果有一小块极亮目标（比如云边缘、建筑物屋顶），min-max 归一化会把整个波段的动态范围压缩得很厉害，导致正常地表像元在归一化后数值集中在 0.9 以上，区分度极差。3σ 截断相当于把极端值剔除后再做线性拉伸，鲁棒性要好得多。这个思路在遥感反演里很常见，但多数人直接用normalizedDifference之后就不管量纲了，很少去做这种 band-wise 的归一化。

4.2 代码里容易出错的三个地方

第一个是reduceRegion的 reducer 组合方式。ee.Reducer.mean().combine(ee.Reducer.stdDev(), null, true)的第三个参数true表示两个 reducer 共享输入，输出字段名会以_mean和_stdDev后缀区分。注意是_stdDev，不是_std，下面mean_std.get(name.cat('_stdDev'))必须和这个后缀完全一致，差一个字母就会报 null 引用错误。

第二个是image.bandNames().map()这个写法。bandNames()返回的是一个 ee.List，对它做map时，回调函数里的name是一个 ee.String 对象，必须用ee.String(name)包一层才能调用cat方法。这个细节新手经常漏，导致回调里name.cat('_mean')报错。

第三个是toBands()的用法。ee.ImageCollection.fromImages(...)把每个波段的结果组装成一个 ImageCollection，然后再用toBands()把集合变成多波段影像。如果不做rename(image.bandNames())，生成的波段名会是类似"0_NDVI"、"1_EVI"这样的前缀格式，后续用select("NDVI")就找不到了。这一步的 rename 是必须的。

4.3 maxPixels 和 scale 的选择逻辑

maxPixels: 10e9表示这个 reducer 最多处理 100 亿个像元。如果你的研究区很大（比如整个省份），且 scale 设成 30 米，那么像元总量很可能超过这个数。一个更稳妥的做法是把maxPixels设成 1e13，或者用bestEffort: true让它自动降低采样分辨率。默认的scale参数如果从影像本身的尺度来算，Landsat SR 是 30 米，但如果你的 geometry 边界范围非常大，一次 reduceRegion 的计算量会非常大，很容易超时。实际使用中我会先print(geometry.area())估算一下面积，再决定 scale 是否放宽到 100 米或者 250 米。

还有一点，这段代码对tileScale做了注释（// tileScale: 16），说明作者遇到过计算超时的问题。tileScale是 GEE 用来控制每个计算任务拆分粒度的参数，默认是 1，调大到 16 可以让任务拆得更碎，从而避免单个 worker 内存溢出。但代价是计算时间变长。如果你的区域特别大，建议先默认跑一次，如果报Computation timed out再把这个参数打开。

5. 避坑指南：C02 迁移里最常见的五个翻车现场

5.1 报了Band 'SR_B1' does not exist错误

现象：对 Landsat 8 影像执行 maskL457sr 时报错，提示找不到 SR_B1。

原因：Landsat 8 的波段从 SR_B2 开始，没有 SR_B1。你如果对 L8 集合用了 L5/7 的 mask 函数，select('SR_B.')倒不会报错，但如果你在 rename 时用了['SR_B1','SR_B2',...]这样的显式波段名，就会直接出错。这份代码里 maskL457sr 和 maskL8sr 已经区分了这种情况，但你自己扩展其他传感器（比如 Landsat 9）时容易踩同一个坑。

解决：每新建一个集合，先print(collection.first().bandNames())确认实际波段名，再决定用哪个 mask 和 rename 函数。Landsat 9 的波段结构和 Landsat 8 完全一致，可以直接复用 maskL8sr 和 rename89。

5.2 归一化结果出现大片纯黑色区域

现象：归一化后的影像在某个区域全是 0 值，直方图在 0 处有个巨大的峰。

原因：mean ± 3σ的下界把该波段的低值全部截断为 0。如果某个波段的直方图本身是严重左偏的（比如水体区域的 NIR 值普遍很低），3σ 下界可能远大于实际的最小值，导致大量像元被截断为 0，看起来就是一片黑色。这不是代码 bug，是统计假设的问题：3σ 截断的前提是数据近似正态分布，对双峰分布（比如水体加植被的混合区域）并不适用。

解决：在归一化前先看直方图，如果某个波段的分布不是单峰的，可以改用分位数截断（比如 1%—99%），或者对不同的土地覆盖类型分别归一化。另一个办法是把 σ 倍数从 3 降到 2，牺牲一些动态范围来保留更多像元。

5.3 Landsat 7 条带区的 NDVI 异常跳变

现象：时间序列里 2003 年之后的 NDVI 曲线上，某些年份出现规则的空洞或者数值骤降。

原因：Landsat 7 的 SLC 故障导致 2003 年 5 月之后所有影像都存在约 22% 的数据缺失（条带）。这份代码的 QA_PIXEL 掩膜只处理了云和云影，没有处理 SLC 条带，所以条带区在计算 NDVI 时会出现黑色条带，时间序列分析时会把这些空洞当成 0 值参与统计。

解决：对 Landsat 7 集合增加一步updateMask(image.select('SR_B3').mask().and(image.select('SR_B4').mask()))，把缺失条带直接掩掉。更精细的做法是用landsat7.simpleCloudScore()做云评分，但那个是 TOA 数据的方法，对 SR 数据不适用。最实用的是在时间序列合成时用median()合成器，条带区在多时相取中位数时会被其他时相的有效像元填补。

5.4 merge 之后影像数量比预期多出好几倍

现象：打印 collection 的 size，发现某个年份有几十景影像，远超该年的实际覆盖次数。

原因：Landsat 5 和 Landsat 7 的时间范围有重叠（2012 年前后），同一区域可能在同一天被两颗卫星都拍到。merge 只是把两个集合按顺序拼接，不会按日期去重。如果你不关心具体是哪颗卫星，直接col.first()取第一景倒没影响，但如果做时间序列，重叠日期会造成重复采样。

解决：在 merge 之后加一步ee.ImageCollection(col.distinct(['system:time_start']))按时间戳去重。更稳妥的是按日期区间分开处理，比如 1985—1999 只用 Landsat 5，2000—2012 用 Landsat 5/7 合并，2013 之后用 Landsat 8/9，这样每颗卫星负责一个时期，衔接处用平均值过渡。

5.5 用col.first()取出来的影像不是你想要的那一景

现象：在 Map 上显示的影像，时间戳和你的研究时段完全对不上。

原因：ee.ImageCollection默认按system:time_start排序，first()返回的是时间上最早的一景，而不是最近的一景，也不是云量最少的那一景。如果你忘了做日期筛选，first()可能取到 1985 年的第一景影像。

解决：取影像前先按日期过滤，再按云量排序。常见做法是：

var col = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2') .filterBounds(geometry) .filterDate('2023-06-01', '2023-09-30') .sort('CLOUD_COVER'); var image = col.first();

按CLOUD_COVER升序排序后取第一景，就是研究时段内云量最少的那一景。如果你想取中位数合成影像，直接col.median()就能得到一个无缝的合成结果，但median()会改变像素值的统计特性，做定量反演时不推荐。

6. 把归一化函数封装成通用模块：批量处理与导出

这份教程的归一化函数在单景影像上表现不错，但要真正投入生产，还得把它封装成一个可复用的模块，支持批量导出和跨研究区复用。我会在 GEE 的代码编辑器里把它保存为一个函数文件，然后通过require的方式引入到其他脚本里。

一个完整的批处理框架大致是这样的：

// 封装：从影像集合生成归一化后的影像集合 function normalizeCollection(collection, region, scale) { var normalized = collection.map(function(img) { var clipped = img.clip(region); return normalization(clipped, region, scale); }); return normalized; } // 导出：按年合成并导出为 GeoTIFF var yearly = ee.ImageCollection.fromImages( ee.List.sequence(1985, 2024).map(function(year) { year = ee.Number(year); var start = ee.Date.fromYMD(year, 1, 1); var end = start.advance(1, 'year'); var annual = col.filterDate(start, end).select(['NDVI']); var composite = annual.reduce(ee.Reducer.median()); return composite.set('system:time_start', start).set('year', year); }) ); Export.table.toDrive({ collection: yearly, description: 'NDVI_annual_median', folder: 'GEE_exports', scale: 30, region: geometry, maxPixels: 1e13 });

这个导出模板里的关键点：select(['NDVI'])在合成前只保留目标指数，减少内存占用；reduce(ee.Reducer.median())每年合成一景中值影像，消除单景的云残留和噪声；set('system:time_start', start)让每景年度影像带上时间戳，方便后续做时序分析。

实际用的时候还有几个参数值得调：

• scale保持 30 米会导出大量数据，如果是省级范围，建议改成 100 米甚至 250 米导出，处理速度会快一个量级。
• maxPixels在导出时需要给足，建议 1e13，否则大区域导出会报错。
• Export.table.toDrive的description会作为文件名前缀，最好带上年份和指数名，避免覆盖。

我现在自己做长时序分析，已经习惯把这套流程固定下来：先print集合的 size 和波段名，再抽出一景影像跑一次归一化，把归一化前后的直方图对比打出来确认截断效果，最后再批量合成和导出。这套模板最大的价值不是那几行归一化代码，而是它把「多源数据对齐」这个步骤标准化了——从那以后我每次拿到新的合成产品，都强制走一遍去云、重命名、指数计算、归一化、导出这个流程，省掉了大量来回调试的时间。希望帮到你。

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