基于 Google Earth Engine （GEE） 的区域植被景观指标计算与分析

一、引言
本文将详细介绍如何使用 GEE 对南京江宁区（以特定多边形区域为例）2010 - 2015 年的植被景观进行分析，计算一系列景观指标，并对结果进行可视化和导出。

2.1 研究区域

var ROI = ee.Geometry.Polygon([
[[118.60, 31.80], [118.60, 32.00], [118.90, 32.00], [118.90, 31.80]]
]);

2.2 数据与时间范围

var startYear = 2010;
var endYear   = 2015;

三、数据预处理

3.1 云掩膜处理

由于卫星影像可能受到云层和云影的影响，需要对数据进行云掩膜处理。通过 maskLandsat 函数，基于 QA_PIXEL 波段去除云及云影覆盖的区域。

function maskLandsat(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 3).eq(0)  // 云
             .and(qa.bitwiseAnd(1 << 4).eq(0)); // 云影
  return image.updateMask(mask);
}

3.2 波段缩放与重命名

// Landsat 5 缩放波段并重命名
function scaleL5(image) {
  return image.select(
      ['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B7']
    )
    .multiply(0.0000275).add(-0.2)
    .rename(['blue','green','red','nir','swir1','swir2'])
    .copyProperties(image, ['system:time_start']);
}

// Landsat 8 缩放波段并重命名
function scaleL8(image) {
  return image.select(
      ['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7']
    )
    .multiply(0.0000275).add(-0.2)
    .rename(['blue','green','red','nir','swir1','swir2'])
    .copyProperties(image, ['system:time_start']);
}

分别对 Landsat 5 和 Landsat 8 数据的波段进行缩放和重命名，以统一波段名称和数据范围。

3.3 年度合成

function getAnnualComposite(year) {
  year = ee.Number(year);
  var start = ee.Date.fromYMD(year,   1, 1);
  var end   = ee.Date.fromYMD(year.add(1), 1, 1);

  var col5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
               .filterDate(start, end)
               .filterBounds(ROI)
               .map(maskLandsat).map(scaleL5);
  var col8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
               .filterDate(start, end)
               .filterBounds(ROI)
               .map(maskLandsat).map(scaleL8);

  var merged = col5.merge(col8);

  var fallback = ee.Image.constant([0.1,0.1,0.1,0.1,0.1,0.1])
                    .rename(['blue','green','red','nir','swir1','swir2'])
                    .reproject('EPSG:4326', null, 30);

  return merged
    .merge(ee.ImageCollection([fallback]))
    .median()
    .clip(ROI);
}

为了得到每年的影像数据，将 Landsat 5 和 Landsat 8 数据合并，并添加一个 “回退” 影像，以确保 median 操作不会因空集而失败。

四、景观指标计算

4.1 计算 NDVI 与植被掩膜

根据年度合成影像计算归一化植被指数（NDVI），并通过设定的 NDVI 阈值提取植被区域。

var ndviThreshold = 0.3;
function computeMetrics(year) {
  year = ee.Number(year);
  var comp = getAnnualComposite(year);
  var ndvi = comp.normalizedDifference(['nir','red']).rename('NDVI');
  var veg  = ndvi.gt(ndviThreshold).selfMask();

4.2 总植被面积计算

通过 ee.Image.pixelArea() 函数计算每个像素的面积，并根据植被掩膜求和得到总植被面积。

  var vegArea = ee.Number(
    ee.Image.pixelArea().updateMask(veg)
      .reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: ROI,
        scale: 30, maxPixels: 1e13
      }).get('area')
  );

4.3 斑块提取与景观指标计算

如果总植被面积大于 0，则进行斑块提取和景观指标计算，包括斑块数量、斑块密度、香农多样性指数，以及新增的平均斑块面积、最大斑块面积和斑块面积标准差。

  return ee.Algorithms.If(
    vegArea.lte(0),
    ee.Feature(null, {
      year:               year,
      patchNumber:        0,
      patchIndex:         0,
      shannonIndex:       0,
      totalVegArea:       0,
      avgPatchArea:       0,
      maxPatchArea:       0,
      patchAreaStdDev:    0
    }),
    (function() {
      // 连通斑块
      var patches = veg.connectedComponents({
        connectedness: ee.Kernel.square(1),
        maxSize: 256
      }).select('labels').rename('patch');

      // 斑块数量、密度
      var patchCount = ee.Number(
        patches.reduceRegion({
          reducer: ee.Reducer.countDistinct(),
          geometry: ROI,
          scale: 30, maxPixels: 1e13
        }).get('patch')
      );
      var regionHa = ROI.area().divide(10000);
      var patchIndex = patchCount.divide(regionHa);

      // 斑块面积分布（列表 of {patch, sum}）
      var groups = ee.List(
        ee.Dictionary(
          ee.Image.pixelArea().addBands(patches)
            .reduceRegion({
              reducer: ee.Reducer.sum().group({
                groupField: 1,
                groupName: 'patch'
              }),
              geometry: ROI,
              scale: 30, maxPixels: 1e13
            })
        ).get('groups')
      );

      // 香农多样性
      var shannon = ee.Number(
        groups.map(function(g) {
          var d = ee.Dictionary(g);
          var area = ee.Number(d.get('sum'));
          var p = area.divide(vegArea).max(1e-6);
          return p.multiply(p.log());
        }).iterate(function(val, acc) {
          return ee.Number(acc).add(ee.Number(val));
        }, 0)
      ).multiply(-1);

      // 计算新指标：平均斑块面积、最大斑块面积、斑块面积标准差
      var areaList = groups.map(function(g) {
        return ee.Number(ee.Dictionary(g).get('sum'));
      });
      var avgArea = vegArea.divide(patchCount);
      var maxArea = ee.Number(areaList.reduce(ee.Reducer.max()));
      var variance = ee.Number(
        areaList.iterate(function(a, acc) {
          a = ee.Number(a);
          return ee.Number(acc)
                   .add(a.subtract(avgArea).pow(2));
        }, 0)
      ).divide(patchCount);
      var stdDev = variance.sqrt();

      return ee.Feature(null, {
        year:            year,
        patchNumber:     patchCount,
        patchIndex:      patchIndex,
        shannonIndex:    shannon,
        totalVegArea:    vegArea,
        avgPatchArea:    avgArea,
        maxPatchArea:    maxArea,
        patchAreaStdDev: stdDev
      });
    })()
  );
}

五、批量处理与结果展示

5.1 批量计算景观指标

使用 ee.FeatureCollection(years.map(computeMetrics)) 对 2010 - 2015 年的景观指标进行批量计算。

var years = ee.List.sequence(startYear, endYear);
var metricsFC = ee.FeatureCollection(years.map(computeMetrics));
print('2010–2015 年度景观指标（含新增指标）', metricsFC);

5.2 可视化示例

以 2025 年为例，展示 NDVI、植被掩膜和植被斑块的可视化结果。

var visYear = 2025;
var comp2025 = getAnnualComposite(visYear);
var ndvi2025 = comp2025.normalizedDifference(['nir','red']).rename('NDVI');
var veg2025  = ndvi2025.gt(ndviThreshold).selfMask();
var patch2025= veg2025.connectedComponents({
  connectedness: ee.Kernel.square(1), maxSize: 256
}).randomVisualizer();

Map.centerObject(ROI, 10);
Map.addLayer(ndvi2025, {min:0, max:1, palette:['white','green']}, 'NDVI ' + visYear);
Map.addLayer(veg2025,  {palette:['blue']},               'Vegetation Mask ' + visYear);
Map.addLayer(patch2025,{ },                              '植被斑块 ' + visYear);

5.3 结果导出

将计算得到的景观指标数据导出为 CSV 文件到 Google Drive。

Export.table.toDrive({
  collection:  metricsFC,
  description: 'LandscapeMetrics_ROI_2010_2015_Expanded',
  fileFormat:  'CSV'
});

六、结果展示：
完整代码：

// 1. 区域、时间与全局参数
var ROI = ee.Geometry.Polygon([
[[118.60, 31.80], [118.60, 32.00], [118.90, 32.00], [118.90, 31.80]]
]);
var startYear = 2010;
var endYear   = 2015;
var ndviThreshold = 0.3;
var years = ee.List.sequence(startYear, endYear);

// 2. 云掩膜：基于 QA_PIXEL 波段
function maskLandsat(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 3).eq(0)  // 云
             .and(qa.bitwiseAnd(1 << 4).eq(0)); // 云影
  return image.updateMask(mask);
}

// 3a. Landsat 5 缩放波段并重命名
function scaleL5(image) {
  return image.select(
      ['SR_B1','SR_B2','SR_B3','SR_B4','SR_B5','SR_B7']
    )
    .multiply(0.0000275).add(-0.2)
    .rename(['blue','green','red','nir','swir1','swir2'])
    .copyProperties(image, ['system:time_start']);
}

// 3b. Landsat 8 缩放波段并重命名
function scaleL8(image) {
  return image.select(
      ['SR_B2','SR_B3','SR_B4','SR_B5','SR_B6','SR_B7']
    )
    .multiply(0.0000275).add(-0.2)
    .rename(['blue','green','red','nir','swir1','swir2'])
    .copyProperties(image, ['system:time_start']);
}

// 4. 年度合成：合并 Landsat 5/8，并追加“回退”影像，保证 median 不空集
function getAnnualComposite(year) {
  year = ee.Number(year);
  var start = ee.Date.fromYMD(year,   1, 1);
  var end   = ee.Date.fromYMD(year.add(1), 1, 1);

  var col5 = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
               .filterDate(start, end)
               .filterBounds(ROI)
               .map(maskLandsat).map(scaleL5);
  var col8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
               .filterDate(start, end)
               .filterBounds(ROI)
               .map(maskLandsat).map(scaleL8);

  var merged = col5.merge(col8);

  var fallback = ee.Image.constant([0.1,0.1,0.1,0.1,0.1,0.1])
                    .rename(['blue','green','red','nir','swir1','swir2'])
                    .reproject('EPSG:4326', null, 30);

  return merged
    .merge(ee.ImageCollection([fallback]))
    .median()
    .clip(ROI);
}

// 5. 斑块提取与景观指标计算（新增三个指标）
function computeMetrics(year) {
  year = ee.Number(year);
  var comp = getAnnualComposite(year);
  var ndvi = comp.normalizedDifference(['nir','red']).rename('NDVI');
  var veg  = ndvi.gt(ndviThreshold).selfMask();

  // 总植被面积（m2）
  var vegArea = ee.Number(
    ee.Image.pixelArea().updateMask(veg)
      .reduceRegion({
        reducer: ee.Reducer.sum(),
        geometry: ROI,
        scale: 30, maxPixels: 1e13
      }).get('area')
  );

  // 若无植被，直接返回所有指标为 0
  return ee.Algorithms.If(
    vegArea.lte(0),
    ee.Feature(null, {
      year:               year,
      patchNumber:        0,
      patchIndex:         0,
      shannonIndex:       0,
      totalVegArea:       0,
      avgPatchArea:       0,
      maxPatchArea:       0,
      patchAreaStdDev:    0
    }),
    (function() {
      // 连通斑块
      var patches = veg.connectedComponents({
        connectedness: ee.Kernel.square(1),
        maxSize: 256
      }).select('labels').rename('patch');

      // 斑块数量、密度
      var patchCount = ee.Number(
        patches.reduceRegion({
          reducer: ee.Reducer.countDistinct(),
          geometry: ROI,
          scale: 30, maxPixels: 1e13
        }).get('patch')
      );
      var regionHa = ROI.area().divide(10000);
      var patchIndex = patchCount.divide(regionHa);

      // 斑块面积分布（列表 of {patch, sum}）
      var groups = ee.List(
        ee.Dictionary(
          ee.Image.pixelArea().addBands(patches)
            .reduceRegion({
              reducer: ee.Reducer.sum().group({
                groupField: 1,
                groupName: 'patch'
              }),
              geometry: ROI,
              scale: 30, maxPixels: 1e13
            })
        ).get('groups')
      );

      // 香农多样性
      var shannon = ee.Number(
        groups.map(function(g) {
          var d = ee.Dictionary(g);
          var area = ee.Number(d.get('sum'));
          var p = area.divide(vegArea).max(1e-6);
          return p.multiply(p.log());
        }).iterate(function(val, acc) {
          return ee.Number(acc).add(ee.Number(val));
        }, 0)
      ).multiply(-1);

      // 计算新指标：平均斑块面积、最大斑块面积、斑块面积标准差
      var areaList = groups.map(function(g) {
        return ee.Number(ee.Dictionary(g).get('sum'));
      });
      var avgArea = vegArea.divide(patchCount);
      var maxArea = ee.Number(areaList.reduce(ee.Reducer.max()));
      var variance = ee.Number(
        areaList.iterate(function(a, acc) {
          a = ee.Number(a);
          return ee.Number(acc)
                   .add(a.subtract(avgArea).pow(2));
        }, 0)
      ).divide(patchCount);
      var stdDev = variance.sqrt();

      return ee.Feature(null, {
        year:            year,
        patchNumber:     patchCount,
        patchIndex:      patchIndex,
        shannonIndex:    shannon,
        totalVegArea:    vegArea,
        avgPatchArea:    avgArea,
        maxPatchArea:    maxArea,
        patchAreaStdDev: stdDev
      });
    })()
  );
}

// 6. 批量处理并打印
var metricsFC = ee.FeatureCollection(years.map(computeMetrics));
print('2010–2025 年度景观指标（含新增指标）', metricsFC);

// 7. 2025 年可视化示例
var visYear = 2025;
var comp2025 = getAnnualComposite(visYear);
var ndvi2025 = comp2025.normalizedDifference(['nir','red']).rename('NDVI');
var veg2025  = ndvi2025.gt(ndviThreshold).selfMask();
var patch2025= veg2025.connectedComponents({
  connectedness: ee.Kernel.square(1), maxSize: 256
}).randomVisualizer();

Map.centerObject(ROI, 10);
Map.addLayer(ndvi2025, {min:0, max:1, palette:['white','green']}, 'NDVI ' + visYear);
Map.addLayer(veg2025,  {palette:['blue']},               'Vegetation Mask ' + visYear);
Map.addLayer(patch2025,{ },                              '植被斑块 ' + visYear);

// 8. 导出至 Google Drive（CSV）
Export.table.toDrive({
  collection:  metricsFC,
  description: 'LandscapeMetrics_ROI_2010_2025_Expanded',
  fileFormat:  'CSV'
});