1 我们为什么要用EVI而不是NDVI在进行长时间序列的生态环境监测时很多同学首选 NDVI归一化植被指数。但在面对鄱阳湖流域这样水汽充沛、夏季植被茂密的区域时NDVI 往往会暴露两个短板高植被区饱和当植被覆盖度很高时NDVI 容易“封顶”Saturate很难区分茂密森林和非常茂密的森林。大气与土壤干扰气溶胶雾霾、水汽和土壤背景色会影响红波段的反射率。EVI增强植被指数就是为了解决这些问题而生的。它在公式中引入了**蓝波段Blue Band**来修正大气气溶胶的影响并加入了土壤调节参数L, C1, C2。简单来说EVI 比 NDVI 更抗干扰更能反映高覆盖度植被的真实生长状况。2 GEE代码分步讲解1研究区域初始化首先我们需要定义感兴趣的区域ROI。这里以江西省鄱阳湖流域为例。// 1. 导入研究区域 // 这里加载你上传的矢量边界如江西省或鄱阳湖流域 var roi ee.FeatureCollection(users/你的用户名/你的资产名称) .filterMetadata(NAME, equals, 江西); // 或者直接使用几何体绘制工具画框 Map.centerObject(roi, 7); Map.addLayer(roi, {color:red}, Study Area);2预处理核心——去云与辐射定标这是处理 Landsat Collection 2 (C02) 数据最关键的一步。Scale FactorsC02 数据为了压缩存储体积将反射率数值进行了缩放。我们需要按官方公式还原。QA_PIXEL利用质量波段进行位运算Bitwise精准剔除云Cloud和云阴影Cloud Shadow。// 2. 定义还原缩放与去云函数 // 还原缩放因子 (Scale Factors) // 公式来源USGS Landsat Collection 2 Product Guide function applyScaleFactors(image){ var opticalBands image.select(SR_B.).multiply(0.0000275).add(-0.2); var thermalBands image.select(ST_B.*).multiply(0.00341802).add(149.0); return image.addBands(opticalBands, null, true) .addBands(thermalBands, null, true); } // 去云掩膜函数 (Cloud Mask) // 筛选 QA_PIXEL 波段确保云和云阴影的位为0 function Mask(image){ var cloudMask (1 3); var cloudshadowMask (1 4); var QA image.select(QA_PIXEL); var mask QA.bitwiseAnd(cloudMask).eq(0) .and(QA.bitwiseAnd(cloudshadowMask).eq(0)); return image.updateMask(mask); }3多传感器解析这是最容易出错的地方Landsat 8 的波段号与 Landsat 5/7 不同。Landsat 5/7: 蓝波段是 B1。Landsat 8: 蓝波段是 B2。我们使用 image.expression 来书写 EVI 公式这样比连续调用 .add() .subtract() 更直观也方便对应波段。处理 Landsat 5 (2000-2011)// 3. 处理 Landsat 5 数据集 (2000-2011) var years_L5 ee.List.sequence(2000, 2011); var L5_COL ee.ImageCollection(years_L5.map(function(y) { var start ee.Date.fromYMD(y, 1, 1); var end start.advance(12, month); var evi_img ee.ImageCollection(LANDSAT/LT05/C02/T1_L2) .filterDate(start, end) .map(applyScaleFactors) // 辐射定标 .map(Mask) // 去云 .map(function(image){ // EVI 计算公式 // 注意L5 蓝波段为 B1 var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(SR_B4), RED: image.select(SR_B3), BLUE: image.select(SR_B1) }).rename(EVI); return image.addBands(evi); }); // 年合成使用中值合成去除极端值并裁剪到研究区 return evi_img.reduce(ee.Reducer.median()).clip(roi) .set(Year, y) .set(System:time_start, start.millis()); }));处理 Landsat 7 (2012 过渡年)注2012年虽然有L7但存在条带丢失问题通常作为补充数据。// 4. 处理 Landsat 7 数据集 var years_L7 ee.List.sequence(2012, 2012); // ... 代码逻辑同 L5因为 L7 波段定义与 L5 一致 ... // (此处省略重复代码逻辑同上仅更换 ImageCollection ID 为 LANDSAT/LE07/C02/T1_L2)处理 Landsat 8 (2013-2023)// 5. 处理 Landsat 8 数据集 (2013-2023) var years_L8 ee.List.sequence(2013, 2023); var L8_COL ee.ImageCollection(years_L8.map(function(y) { // ... 日期定义同上 ... var evi_img ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) // L8集合 .filterDate(start, end) .map(applyScaleFactors) .map(Mask) .map(function(image){ // 注意L8 波段发生漂移Blue变为 B2Red变为 B4NIR变为 B5 var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(SR_B5), RED: image.select(SR_B4), BLUE: image.select(SR_B2) }).rename(EVI); return image.addBands(evi); }); // ... reduce逻辑同上 ... }));4数据融合与时间序列分析将三个独立的 ImageCollection 合并为一个连续的时间序列并利用ui.Chart绘制趋势图。// 6. 数据合并与图表构建 var data ee.ImageCollection(L5_COL.merge(L7_COL).merge(L8_COL)); // 绘制时间序列图 var yearlychart ui.Chart.image.series({ imageCollection: data.select(EVI_median), region: roi, reducer: ee.Reducer.mean(), // 计算区域内的平均EVI scale: 500, xProperty: Year }).setOptions({ title: 鄱阳湖流域 EVI 年际变化趋势 (2000-2023), vAxis: {title: EVI Value}, hAxis: {title: Year, format: ####}, lineWidth: 2, pointSize: 3, trendlines: { 0: {color: red, visibleInLegend: true} } // 添加趋势线 }); print(yearlychart);5空间细节对比为了验证 EVI 对城市扩张和水体的敏感度我们选取南昌市核心区红谷滩-赣江段进行细节对比。这个对比图还是可以看到一些很明显的变化。// 7. 空间分布可视化对比 (南昌局部) // 定义南昌市中心及周边矩形区域 var roi1 ee.Geometry.Polygon([[[115.80, 28.60], [116.05, 28.60], [116.05, 28.75], [115.80, 28.75]]]); Map.centerObject(roi1, 10); Map.addLayer(roi1, {color:blue}, Nanchang ROI); // 提取首尾年份数据 var evi_2000 data.filterMetadata(Year, equals, 2000).first().clip(roi1); var evi_2023 data.filterMetadata(Year, equals, 2023).first().clip(roi1); // 可视化参数红-黄-绿 渐变 var visArgs {min: 0, max: 0.7, palette: [d7191c, fdae61, ffffbf, a6d96a, 1a9641]}; Map.addLayer(evi_2000.select(EVI_median), visArgs, 2000 EVI); Map.addLayer(evi_2023.select(EVI_median), visArgs, 2023 EVI);结果如下3 写在最后你可以更换你的研究区。运行代码后你会得到一张 EVI 变化趋势图和两张地图。在写论文或报告时可以关注以下几点趋势线斜率如果斜率为正说明区域整体变绿生态恢复如果为负可能涉及城市化扩张或退耕。空间差异对比研究区的 2000 和 2023 影像原来是绿色的区域农田/林地如果变红低值说明转为了不透水面城市建设。水体识别EVI 在水体处的值通常极低甚至为负这也是一种辅助提取水体范围的方法。今天的教程就到这里结束希望能帮到你我们下期见欢迎关注梧桐GIS我是加拿大一枝黄花再见
GEE基础教程：Landsat长时序EVI植被指数构建与分析
1 我们为什么要用EVI而不是NDVI在进行长时间序列的生态环境监测时很多同学首选 NDVI归一化植被指数。但在面对鄱阳湖流域这样水汽充沛、夏季植被茂密的区域时NDVI 往往会暴露两个短板高植被区饱和当植被覆盖度很高时NDVI 容易“封顶”Saturate很难区分茂密森林和非常茂密的森林。大气与土壤干扰气溶胶雾霾、水汽和土壤背景色会影响红波段的反射率。EVI增强植被指数就是为了解决这些问题而生的。它在公式中引入了**蓝波段Blue Band**来修正大气气溶胶的影响并加入了土壤调节参数L, C1, C2。简单来说EVI 比 NDVI 更抗干扰更能反映高覆盖度植被的真实生长状况。2 GEE代码分步讲解1研究区域初始化首先我们需要定义感兴趣的区域ROI。这里以江西省鄱阳湖流域为例。// 1. 导入研究区域 // 这里加载你上传的矢量边界如江西省或鄱阳湖流域 var roi ee.FeatureCollection(users/你的用户名/你的资产名称) .filterMetadata(NAME, equals, 江西); // 或者直接使用几何体绘制工具画框 Map.centerObject(roi, 7); Map.addLayer(roi, {color:red}, Study Area);2预处理核心——去云与辐射定标这是处理 Landsat Collection 2 (C02) 数据最关键的一步。Scale FactorsC02 数据为了压缩存储体积将反射率数值进行了缩放。我们需要按官方公式还原。QA_PIXEL利用质量波段进行位运算Bitwise精准剔除云Cloud和云阴影Cloud Shadow。// 2. 定义还原缩放与去云函数 // 还原缩放因子 (Scale Factors) // 公式来源USGS Landsat Collection 2 Product Guide function applyScaleFactors(image){ var opticalBands image.select(SR_B.).multiply(0.0000275).add(-0.2); var thermalBands image.select(ST_B.*).multiply(0.00341802).add(149.0); return image.addBands(opticalBands, null, true) .addBands(thermalBands, null, true); } // 去云掩膜函数 (Cloud Mask) // 筛选 QA_PIXEL 波段确保云和云阴影的位为0 function Mask(image){ var cloudMask (1 3); var cloudshadowMask (1 4); var QA image.select(QA_PIXEL); var mask QA.bitwiseAnd(cloudMask).eq(0) .and(QA.bitwiseAnd(cloudshadowMask).eq(0)); return image.updateMask(mask); }3多传感器解析这是最容易出错的地方Landsat 8 的波段号与 Landsat 5/7 不同。Landsat 5/7: 蓝波段是 B1。Landsat 8: 蓝波段是 B2。我们使用 image.expression 来书写 EVI 公式这样比连续调用 .add() .subtract() 更直观也方便对应波段。处理 Landsat 5 (2000-2011)// 3. 处理 Landsat 5 数据集 (2000-2011) var years_L5 ee.List.sequence(2000, 2011); var L5_COL ee.ImageCollection(years_L5.map(function(y) { var start ee.Date.fromYMD(y, 1, 1); var end start.advance(12, month); var evi_img ee.ImageCollection(LANDSAT/LT05/C02/T1_L2) .filterDate(start, end) .map(applyScaleFactors) // 辐射定标 .map(Mask) // 去云 .map(function(image){ // EVI 计算公式 // 注意L5 蓝波段为 B1 var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(SR_B4), RED: image.select(SR_B3), BLUE: image.select(SR_B1) }).rename(EVI); return image.addBands(evi); }); // 年合成使用中值合成去除极端值并裁剪到研究区 return evi_img.reduce(ee.Reducer.median()).clip(roi) .set(Year, y) .set(System:time_start, start.millis()); }));处理 Landsat 7 (2012 过渡年)注2012年虽然有L7但存在条带丢失问题通常作为补充数据。// 4. 处理 Landsat 7 数据集 var years_L7 ee.List.sequence(2012, 2012); // ... 代码逻辑同 L5因为 L7 波段定义与 L5 一致 ... // (此处省略重复代码逻辑同上仅更换 ImageCollection ID 为 LANDSAT/LE07/C02/T1_L2)处理 Landsat 8 (2013-2023)// 5. 处理 Landsat 8 数据集 (2013-2023) var years_L8 ee.List.sequence(2013, 2023); var L8_COL ee.ImageCollection(years_L8.map(function(y) { // ... 日期定义同上 ... var evi_img ee.ImageCollection(LANDSAT/LC08/C02/T1_L2) // L8集合 .filterDate(start, end) .map(applyScaleFactors) .map(Mask) .map(function(image){ // 注意L8 波段发生漂移Blue变为 B2Red变为 B4NIR变为 B5 var evi image.expression( 2.5 * ((NIR - RED) / (NIR 6 * RED - 7.5 * BLUE 1)), { NIR: image.select(SR_B5), RED: image.select(SR_B4), BLUE: image.select(SR_B2) }).rename(EVI); return image.addBands(evi); }); // ... reduce逻辑同上 ... }));4数据融合与时间序列分析将三个独立的 ImageCollection 合并为一个连续的时间序列并利用ui.Chart绘制趋势图。// 6. 数据合并与图表构建 var data ee.ImageCollection(L5_COL.merge(L7_COL).merge(L8_COL)); // 绘制时间序列图 var yearlychart ui.Chart.image.series({ imageCollection: data.select(EVI_median), region: roi, reducer: ee.Reducer.mean(), // 计算区域内的平均EVI scale: 500, xProperty: Year }).setOptions({ title: 鄱阳湖流域 EVI 年际变化趋势 (2000-2023), vAxis: {title: EVI Value}, hAxis: {title: Year, format: ####}, lineWidth: 2, pointSize: 3, trendlines: { 0: {color: red, visibleInLegend: true} } // 添加趋势线 }); print(yearlychart);5空间细节对比为了验证 EVI 对城市扩张和水体的敏感度我们选取南昌市核心区红谷滩-赣江段进行细节对比。这个对比图还是可以看到一些很明显的变化。// 7. 空间分布可视化对比 (南昌局部) // 定义南昌市中心及周边矩形区域 var roi1 ee.Geometry.Polygon([[[115.80, 28.60], [116.05, 28.60], [116.05, 28.75], [115.80, 28.75]]]); Map.centerObject(roi1, 10); Map.addLayer(roi1, {color:blue}, Nanchang ROI); // 提取首尾年份数据 var evi_2000 data.filterMetadata(Year, equals, 2000).first().clip(roi1); var evi_2023 data.filterMetadata(Year, equals, 2023).first().clip(roi1); // 可视化参数红-黄-绿 渐变 var visArgs {min: 0, max: 0.7, palette: [d7191c, fdae61, ffffbf, a6d96a, 1a9641]}; Map.addLayer(evi_2000.select(EVI_median), visArgs, 2000 EVI); Map.addLayer(evi_2023.select(EVI_median), visArgs, 2023 EVI);结果如下3 写在最后你可以更换你的研究区。运行代码后你会得到一张 EVI 变化趋势图和两张地图。在写论文或报告时可以关注以下几点趋势线斜率如果斜率为正说明区域整体变绿生态恢复如果为负可能涉及城市化扩张或退耕。空间差异对比研究区的 2000 和 2023 影像原来是绿色的区域农田/林地如果变红低值说明转为了不透水面城市建设。水体识别EVI 在水体处的值通常极低甚至为负这也是一种辅助提取水体范围的方法。今天的教程就到这里结束希望能帮到你我们下期见欢迎关注梧桐GIS我是加拿大一枝黄花再见
