简介这份PDF教程面向需要在Google Earth EngineGEE中处理长时序Landsat遥感影像的科研人员和GIS学习者聚焦1985—2024年间NDVI、EVI、SAVI、NDMI等常用植被指数的归一化流程并已更新适配Landsat Collection 2数据集。教程针对Landsat 5/7/8不同传感器进行了数据整合完成原始波段的统一处理与波段名称重写优化了去云、辐射定标等预处理函数使多源影像能以一致口径参与指数计算。作者进一步简化了不同集合中的指数调用逻辑研究者可按需选取1985年至今任意时期影像快速完成NDVI、EVI等指数的批量归一化。资源为1个PDF文件大小约900KB内部包含核心代码、函数模块和操作说明方便离线阅读与代码复用。目前已有977人学习下载适合植被遥感、生态监测与长时序地表变化分析方向的入门及进阶用户。通过该文档读者可以系统掌握QA_PIXEL去云、SR波段缩放、波段重命名以及NDVI、EVI、NDWI、NDBI等指数在GEE中的实现方法大幅提升Landsat数据的预处理效率与批量计算稳定性也为后续时间序列分析和趋势研究提供可复用的代码基础。1. 用GEE归一化NDVI/EVI/SAVI/NDMILandsat C02 数据预处理与计算用GEE做长时序NDVI的人多半有同一个体会数据源源不断但每次换一颗卫星就要重写一遍预处理做完NDVI还想补EVI、SAVI、NDMI却发现Landsat 5/7和8在C02版本里的波段名、缩放因子、去云位完全对不上。这份GEE在线归一化教程把1985-2024年的Landsat 5/7/8数据整合成同一套流程统一完成去云、缩放、波段重命名、指数计算和均值±3倍标准差归一化解决的就是“想复现却要重复造轮子”的问题。我把代码里的每个函数都跑了一遍确认它可以按需调用任一时期的影像做归一化输出适合做长时序植被覆盖、土壤调节植被指数和水体指数分析的人直接复用。2. 先把数据拉齐C02 数据的去云、缩放与波段重命名Landsat C02 相比 C01 最大的变化是 Level-2 产品直接提供地表反射率但波段命名同时改了Landsat 5/7 的 SR 波段为 SR_B1 到 SR_B7Landsat 8/9 为 SR_B2 到 SR_B7热红外波段也从 ST_B6 变成 ST_B10/ST_B11。如果不统一命名每新增一个指数就要给不同卫星写不同的波段选择逻辑代码量直接翻倍。这份资源的关键做法就是把数据准备阶段切成“去云掩膜、缩放因子、波段重命名”三个函数先让所有影像在统一波段名上工作再进入指数计算。2.1 掩膜函数用 QA_PIXEL 位运算代替老的 CFmaskC02 数据的云掩膜不再推荐用 CFmask 的 quality 波段而是用 QA_PIXEL 和 QA_RADSAT。QA_PIXEL 是 16 位整型低 5 位的含义如下表BitLandsat 5/7Landsat 8/90FillFill1Dilated CloudDilated Cloud2UnusedCirrus3CloudCloud4Cloud ShadowCloud Shadow教程里给 L5/L7 写的掩膜函数是整个流程的基础所有 L5/L7 影像在进入指数计算之前都必须经过它function maskL457sr(image) { // QA_PIXEL 低5位Fill/Dilated Cloud/Unused/Cloud/Cloud Shadow var qaMask image.select(QA_PIXEL).bitwiseAnd(parseInt(11111, 2)).eq(0); var saturationMask image.select(QA_RADSAT).eq(0); // 光学波段官方缩放系数 0.0000275 与偏移 -0.2 var opticalBands image.select(SR_B.).multiply(0.0000275).add(-0.2); var thermalBand image.select(ST_B6).multiply(0.00341802).add(149.0); return image.addBands(opticalBands, null, true) .addBands(thermalBand, null, true) .updateMask(qaMask) .updateMask(saturationMask); }代码逻辑分三层。第一层qaMask 读取 QA_PIXEL用 bitwiseAnd 与 31 做位与运算只要低 5 位任何一个位为 1就说明该像元是云、云影、填充或膨胀云eq(0) 表示只保留低 5 位全为 0 的干净像元。第二层saturationMask 读取 QA_RADSAT要求它等于 0即所有波段都没有进入饱和这在高反射裸土和城市区域尤其关键。第三层光学波段用 0.0000275 与 -0.2 两个系数做线性缩放热红外波段用另一组系数 0.00341802 与 149.0然后 addBands 加回原影像。addBands 的第三个参数 true 表示覆盖原始波段名称后面不需要再对缩放后的波段做重命名。maskL8sr 与 maskL457sr 的区别只有两处一是热红外波段正则表达式从 ST_B6 改为 ST_B.*二是 Landsat 8 多了 Cirrus 波段。这份资源里统一掩掉低 5 位Cirrus 也会一并剔除。有人会担心这样误伤薄云覆盖的植被区实际测试中大部分研究区仍能保留足够像元真正影响 NDVI 统计的是厚云和云影所以这套保守策略适合作为默认。2.2 波段重命名把 L5/L7 与 L8 的 SR 波段对齐掩膜和缩放完成后影像保留的是原始波段名。L5 的 SR_B1 对应蓝色L8 的 SR_B1 却是海岸波段两者完全不是一回事。如果直接拿原始波段名算指数L5 和 L8 的 NDVI 公式内部选的波段就不一致最终结果根本没法合并比较。因此下一步必须用重命名函数统一成 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] ); }这里有一个很容易踩的坑Landsat 8 的 SWIR1 是 SR_B6SWIR2 是 SR_B7Landsat 5/7 的 SWIR1 是 SR_B5SWIR2 是 SR_B7。两代卫星的 SR_B7 都叫 SWIR2但 L5 的 SR_B5 和 L8 的 SR_B6 才是 SWIR1。如果照抄 L5 的映射到 L8后面算 NDBI 时会把两个短波红外波段位置调换建筑物和水体的区分度会完全错乱。我每次改完这套函数都会先执行一次 print(collection1.first().bandNames()) 确认波段顺序无误再继续往下的指数计算。2.3 为什么用 T1_L2 而不是 T2 或 Level-1教程里三个数据集合分别是 LANDSAT/LT05/C02/T1_L2、LANDSAT/LE07/C02/T1_L2、LANDSAT/LC08/C02/T1_L2。T1 是经过精密几何校正的地形校正产品适合做跨年度时间序列T2 精度低一级在拼接或趋势分析中会出现明显的轨道接缝噪声。L2 表示已经是地表反射率省去了大气校正环节。如果换成 Level-1 的 TOA 数据缩放系数要从反射率定标那套重新算而且没有 QA_PIXEL 这么干净的位掩码等于把简单事情复杂化。三颗卫星的时间边界也要留意LT05 覆盖 1984-2012LE07 覆盖 1999-至今LC08 覆盖 2013-至今。教程里把 L5 和 L7 放在同一个 1985-2012 的日期段L8 单独放 2013-2024最大程度避免同一日期出现两景不同传感器的影像。真正可用的覆盖是 1985 年到 2012 年以 L5/L7 为主2013 年以后以 L8 为主这个分段方式本身就值得直接复用。3. 指数计算与数据集合合并同一个函数打通NDVI、EVI、SAVI、NDMI3.1 指数函数库的写法教程里把 NDVI、NDWI、NDBI、EVI 都封装成接收 Image 输出 Image 的函数每个函数用 addBands 把指数加到影像上保留原始光学波段。多函数串成 map 链的好处是后续每处理一景影像都同时得到所有指数波段不需要为每个指数单独写一遍集合遍历。这样写还有一个隐藏优势新增指数时只需要加一个函数再在 map 链里补一个 .map(函数名)完全不用改动前面的掩膜和重命名逻辑。function NDVI(image) { return image.addBands( image.normalizedDifference([nir, red]).rename(NDVI)); } function NDWI(image) { return image.addBands( image.normalizedDifference([red, swir1]).rename(NDWI)); } function NDBI(image) { return image.addBands( image.normalizedDifference([swir2, swir1]).rename(NDBI)); } 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)); }NDVI、NDWI、NDBI 用的都是 normalizedDifference只需在数组里调换波段名。教程里 NDWI 用的是 red 和 swir1 的组合Gao 版本它突出的是植被水分信息想要 McFeeters 版本的水体指数把数组改成 [green, nir] 即可。EVI 没有用 normalizedDifference是因为公式里有 6、7.5、2.5 这些系数expression 比手工乘除更不容易写错。标题里提到的 SAVI 和 NDMI 可以沿同一套模式扩展。SAVI 的 L 值一般取 0.5适用于中等植被覆盖区域如果研究区是稀疏灌丛L 可以降到 0.25。NDMI 就是归一化差异水分指数用 nir 和 swir1在干旱半干旱区NDMI 比 NDVI 更容易看出季节水分变化function SAVI(image) { var savi image.expression( (1 0.5) * ((NIR - RED) / (NIR RED 0.5)), { NIR: image.select(nir), RED: image.select(red) }); return image.addBands(savi.rename(SAVI)); } function NDMI(image) { return image.addBands( image.normalizedDifference([nir, swir1]).rename(NDMI)); }这两个函数可以直接插入原来的 map 链不需要改动任何预处理函数。SAVI 在裸地占比高的区域比 NDVI 稳定NDMI 在含水植被监测里比 NDVI 敏感两者互补性很强。3.2 三个集合的过滤、map 与 merge数据集合的构建是整个流程的骨架。L5 与 L7 的处理链完全一致L8 需要换掩膜函数和重命名函数三者都经过 mask、rename、指数计算后才合并。这一步的顺序不能错具体原因会在避坑章节展开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) .map(SAVI).map(NDMI); 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) .map(SAVI).map(NDMI); 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) .map(SAVI).map(NDMI); var col collection1.merge(collection2).merge(collection3); var image col.first().clip(geometry);filterDate 的日期范围决定了你能覆盖的时间长度。要处理 2024 年的数据把 collection3 的 end 改成 2024-12-31要把 L7 的早期数据也纳入可以单独给 collection2 设一个 1999-01-01 的起点。merge 要求三个集合所有影像的波段名完全一致这一步依赖前面 rename57/rename89 的成功执行。如果某一步忽略重命名merge 后取 first() 不会马上报错但 export 或 reduceRegion 时会出现波段缺失。map 链的执行顺序不能乱mask 必须在 rename 之前因为缩放时要用原始波段名rename 必须在指数计算之前因为指数函数里 select 的是统一的 blue、red、nir。虽然 GEE 的 map 是延迟执行但函数逻辑上依然是严格序列化的。我合并完三个集合后会马上跑一次 col.size().getInfo() 和 col.first().bandNames().getInfo()第一时间发现空集合或波段错位避免在后续大型计算里反复出 bug。4. 归一化的边界在哪均值±3倍标准差逐波段归一化4.1 归一化到底处理了什么在 GEE 里算完 NDVI 后直接输出数值范围并不是严格的 [0,1]。EVI 在部分地表会超过 1不同年份不同传感器的绝对数值也不可比。为了把多个时间点、多个指数拉到同一尺度教程用了均值±3倍标准差的统计归一化方法。具体做法是对区域内每个波段求均值和标准差把均值减 3 倍标准差作为下限均值加 3 倍标准差作为上限超过上限的裁成 1、低于下限的裁成 0中间部分线性拉伸到 [0,1]。这种方法的优点是对离群值不敏感因为 3σ 之外的数据被直接压平缺点是被裁掉的像元失去了内部差异。如果研究区是大范围水体加城市混合统计均值容易被极端值拖偏后面避坑章节会专门讲怎么处理这种情况。4.2 逐波段归一化代码拆解normalization 函数是这份资源的另一个核心它一次性完成所有波段的统计与映射不需要手动指定波段名称。函数接收三个参数image 是要处理的单景影像region 是统计范围scale 是采样分辨率function normalization(image, region, scale) { // 同时计算均值与标准差输出键为“波段名_mean”和“波段名_stdDev” 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)); // 小于下限的置为min大于上限的置为max中间保留原值 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)))); // 线性拉伸到 [0,1] var result_band band1.subtract(min).divide(max.subtract(min)); return result_band; }) ).toBands().rename(image.bandNames()); return unitScale; }代码有五个关键点。第一reduceRegion 的 reducer 是 mean 和 stdDev 的组合combine 的第三个参数传 true 表示两个统计量共用同一份像元采样避免重复计算。第二maxPixels 设成 10e9也就是 100 亿像元覆盖大区域统计时不至于算出部分结果就报 Too many pixels。第三band1 表达式做的是上下限裁剪band.lt(min) 返回 0/1 掩膜乘以 min 后把低于下限的像元置为 minband.gt(max) 同理中间项保留原值。第四result_band 用 (x-min)/(max-min) 把值域映射到 [0,1]。第五ImageCollection.fromImages 把每个波段的结果包成集合再 toBands最后 rename 回原波段名保证输出影像的波段名与输入一致。实际调用示例var normal_image normalization(image, geometry, 1000); Map.addLayer(normal_image, {bands: [EVI, NDVI, NDWI], min: 0, max: 1}, normalized);这里的 scale 是 reduceRegion 的采样尺度。Landsat 原始分辨率 30m我一般传 30 或 100。传 1000 会显著加快统计速度但统计结果用的是金字塔重采样数据均值标准差会被平滑归一化结果轻微偏暗。小区域直接用 30中等区域用 100省级以上大范围才考虑 1000并且要结合 tileScale 一起调。注意bandNames() 的顺序决定了最终影像的波段顺序如果某个影像缺少某个波段整条 toBands 链会因此断掉所以前面一定确保每个集合的波段完整。4.3 归一化前后效果对比教程用 ui.Chart.image.histogram 对比了归一化前后的 nir 和 EVI。归一化前直方图通常没有铺满横轴大部分像元挤在某个狭窄区间归一化后直方图占满 [0,1] 的整个范围。影像显示上归一化前整体偏灰或偏暗归一化后层次分明。这个对比对实际项目很重要。我在做多年 NDVI 合成图时如果直接输出原始 NDVI1985 年的数据和 2015 年的数据因为传感器差异整体数值会系统性偏移 0.05 到 0.1做差值分析时会出现虚假的“趋势”。归一化之后分析尺度一致差值才能解释为真实的植被变化。5. 避坑指南复现这份GEE代码时踩过的四个坑5.1 现象merge 之后 first() 影像波段对不上我最早复现时把 rename57 和 rename89 的顺序写反L8 集合没有重命名就直接 map 指数函数。结果是 col.first() 打印出来正常但导出时 GEE 提示 Band nir not found。原因L8 原始波段名是 SR_B2 到 SR_B7而 NDVI 函数里 select(nir) 在这个影像上根本选不到因为影像里没有这个名。merge 其实成功执行了三个集合的波段列表不一致导出时才暴露出缺失。解决把 map 链的顺序固定为 mask - rename - index并在 merge 后用 print(col.first().bandNames()) 确认六个光学波段和六个指数波段都在。5.2 现象QA_PIXEL 位掩码把边缘像元全部抹掉有一阵子我的影像在研究区边界出现大量空洞山区阴影区域全是 0。原因是照搬了 parseInt(11111, 2)把 Dilated Cloud 也当成云处理。Dilated Cloud 是云边缘的缓冲区域很多是薄云或半透明云气对 NDVI 均值影响有限全部掩掉会让有效像元数量减少 10% 到 20%。解决改成 parseInt(11000, 2) 只掩 Cloud 和 Cloud Shadow。如果研究区云量本来就高还可以进一步保留 Cirrus只掩 Bit3 和 Bit4。这个参数没有绝对最优取决于研究区的云覆盖类型。5.3 现象归一化结果大量堆积在 0 和 1我用包含大面积水体的区域测试时输出影像出现大范围纯黑和纯白直方图两端堆积严重。原因是水体在 EVI 上接近 0城市建筑区域又很大均值被两个极端拉偏3σ 区间没有覆盖中间的植被分布相当一部分像元被裁到边界。解决归一化前先做水体掩膜把 NDWI 0 的像元排除出统计范围或者干脆在 reduceRegion 之前 clip 到去掉水域的矢量区域。如果想要更稳健的边界用值域两端 2% 的分位数代替 3σ效果比固定标准差好很多。5.4 现象scale 参数随意改导致归一化结果整体偏暗我把 scale 从 30 改到 1000 加速统计结果归一化影像亮度和对比度明显变化直方图形状完全不一样。原因是尺度放大后reduceRegion 采样的是金字塔聚合层的数据均值和标准差被空间平滑3σ 边界随之变化。解决reduceRegion 的 scale 应该与最终产品的使用尺度一致。如果导出 30m 产品就把 scale 设成 30想控制内存用 tileScale: 16 而不是加大 scale。这是资源代码里没有写明的参数取舍。5.5 现象L5 和 L7 重叠期出现同一日期两景影像我把 collection1 的时间范围设成 1985-2013collection3 从 2013 开始结果 2013 年里 L7 和 L8 各有一景同一天的数据做逐日时间序列时影像出现跳变。原因Landsat 7 和 Landsat 8 在 2013 年有将近一年的共同覆盖直接用 merge 合并时同一日期会保留两景影像。解决时间边界切清楚LE07 截止到 2013-01-01LC08 从 2013-01-01 开始或者 merge 后用日期分组取中值。我倾向后者因为 L7 的 SLC-off 数据在中纬度地区仍有较好覆盖率保留它可以增加样本量。6. 把归一化接入长时序趋势分析一段能直接换指数的结尾normalization 函数不止能做单景输出。实际项目中我把它封装成通用函数对每年的合成影像做归一化再对归一化后的 NDVI 序列做线性趋势拟合得到 1985-2024 年的植被变化斜率。这样做能避开单年云覆盖的随机噪声也把三颗传感器之间绝对数值差异的影响降到最低。function annualNDVITrend(col, geometry, startYear, endYear) { var years ee.List.sequence(startYear, endYear); var annual ee.ImageCollection(years.map(function(y) { var year ee.Number(y); var start ee.Date.fromYMD(year, 1, 1); var end ee.Date.fromYMD(year, 12, 31); var yearCol col.filterDate(start, end).map(function(img) { return img.clip(geometry); }); return yearCol.median().select([NDVI, EVI, SAVI, NDMI]).set(year, year); })); return annual; }这个函数把集合按年切分每年返回中值合成影像。注意我显式 select 了四个指数而不是用全部波段因为中值合成对每个波段独立计算波段越多越耗时。之后用 linearFit 对年份和 NDVI 做线性回归var annual annualNDVITrend(col, geometry, 1985, 2024); var withYear annual.map(function(img) { return img.addBands(ee.Image(img.get(year)).rename(year)); }); var slope withYear.select([year, NDVI]).reduce(ee.Reducer.linearFit());linearFit 输出的 scale 波段就是趋势斜率正值表示植被改善负值表示退化。这里有一个容易忽略的点趋势分析必须用归一化后的 NDVI原始 NDVI 在三颗传感器之间存在系统性偏移斜率会被假趋势污染。我做长时序项目还有个习惯最初手动 export 时只导出 slope 和 rms 两个波段在本地打开后叠加研究区边界确认没有异常的条带噪声再继续。这个检查能滤掉很多因云掩膜不彻底造成的虚假趋势。从那以后我每次写 GEE 脚本都先 print 波段、再预览单景、最后跑归一化直方图比对这一条固定流程帮我拦住过无数次波段错位和无意义趋势希望这个习惯也能帮到你。本文还有配套的精品资源点击获取