简介面向遥感与地理信息分析学习者提供基于Google Earth EngineGEE的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、swir2function 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_B6Landsat 8 选的是 SR_B2 到 SR_B7跳过 SR_B1因为那是海岸波段。重命名之后Landsat 8 的 SR_B6短波红外 1会被命名为 swir1SR_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 用的是normalizedDifferenceEVI 用的是expressionfunction 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 52000—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 和波段名再抽出一景影像跑一次归一化把归一化前后的直方图对比打出来确认截断效果最后再批量合成和导出。这套模板最大的价值不是那几行归一化代码而是它把「多源数据对齐」这个步骤标准化了——从那以后我每次拿到新的合成产品都强制走一遍去云、重命名、指数计算、归一化、导出这个流程省掉了大量来回调试的时间。希望帮到你。本文还有配套的精品资源点击获取