GEE提取水体面积结果不对

Question

// ============================================================
// 三峡库区局部河段水体面积提取（巫峡—西陵峡段）
// 微调版：25°坡度 + 轻度形态学 + 浑水阈值 -0.05
// ============================================================

// ================= ROI（核心峡谷段） =================
var roi = ee.Geometry.Polygon([
  [[109.0, 30.0],
   [109.5, 30.0],
   [110.5, 30.5],
   [110.5, 30.75],
   [110.0, 30.8],
   [109.0, 30.3],
   [109.0, 30.0]]
]);

Map.centerObject(roi, 9);
Map.addLayer(roi, {color: 'red'}, 'ROI');

// ================= 参数 =================
var startDate = '2000-01-01';
var endDate   = '2025-12-31';

var MAX_CLOUD = 20;          // 云量阈值
var MIN_PIXELS = 300;        // 小区域最低有效像元数
var minSize = 50;            // 最小水体斑块（像元数），保留河道连续性
var MNDWI_THRESHOLD = -0.05; // 稍微放宽浑水阈值，适应汛期高含沙量
var morphRadius = 1;         // 圆形核半径（30m），轻量去噪

// ================= DEM 坡度掩膜（25°去山体阴影） =================
var dem = ee.Image('USGS/SRTMGL1_003');
var slopeMask = ee.Terrain.slope(dem).lt(25);

// ================= 云/阴影掩膜（所有传感器统一） =================
function maskLandsatSR(image) {
  var qa = image.select('QA_PIXEL');
  var mask = qa.bitwiseAnd(1 << 3).eq(0)   // 云
              .and(qa.bitwiseAnd(1 << 4).eq(0)) // 阴影
              .and(qa.bitwiseAnd(1 << 5).eq(0)); // 雪/冰
  return image.updateMask(mask);
}

// ================= MNDWI 计算 =================
function prep(img) {
  img = maskLandsatSR(img);   // 统一云掩膜

  var isL5L7 = img.bandNames().contains('SR_B2');
  var mndwi = ee.Image(
    ee.Algorithms.If(isL5L7,
      img.normalizedDifference(['SR_B2','SR_B5']),
      img.normalizedDifference(['SR_B3','SR_B6'])
    )
  ).rename('MNDWI');

  return mndwi.copyProperties(img, ['system:time_start']);
}

// ================= 数据源（L7 仅取 SLC-off 前：截止 2003-05-31） =================
var col = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
  .merge(ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
           .filterDate('2000-01-01', '2003-05-31'))
  .merge(ee.ImageCollection('LANDSAT/LC08/C02/T1_L2'))
  .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
  .filterBounds(roi)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUD_COVER', MAX_CLOUD))
  .map(prep);

print('影像数量:', col.size());

// ================= 形态学开运算（轻度） =================
function morph(img) {
  var k = ee.Kernel.circle(morphRadius);
  return img.focal_min({kernel: k, iterations: 1})
            .focal_max({kernel: k, iterations: 1});
}

// ================= 单景处理 =================
function process(img) {
  img = ee.Image(img).clip(roi);

  // 水体二值图（浑水阈值 -0.05）
  var water = img.gt(MNDWI_THRESHOLD).rename('water');

  // 轻度形态学去噪
  water = morph(water);

  // 连通域过滤（maxSize=1024，阈值50），剔除孤立噪点
  var conn = water.connectedPixelCount(1024, true);
  water = water.updateMask(conn.gte(minSize));

  // 坡度掩膜（25°）
  water = water.updateMask(slopeMask);

  // 面积 (km²)
  var area = water.multiply(ee.Image.pixelArea())
                  .reduceRegion({
                    reducer: ee.Reducer.sum(),
                    geometry: roi,
                    scale: 30,
                    maxPixels: 1e13,
                    tileScale: 4
                  }).get('water');
  var area_km2 = ee.Number(area).divide(1e6);

  // 水体比例
  var frac = water.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: roi,
    scale: 30,
    maxPixels: 1e13
  }).get('water');

  // 有效像元总数
  var pixels = img.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: roi,
    scale: 30,
    maxPixels: 1e13
  }).get('MNDWI');

  return ee.Feature(null, {
    'date': ee.Date(img.get('system:time_start')).format('YYYY-MM-dd'),
    'time': img.get('system:time_start'),
    'area': area_km2,
    'waterFrac': frac,
    'pixels': pixels
  });
}

// ================= 执行与过滤 =================
var fc = ee.FeatureCollection(col.map(process))
  .filter(ee.Filter.notNull(['area']))
  .filter(ee.Filter.gt('pixels', MIN_PIXELS))
  .filter(ee.Filter.gt('waterFrac', 0.005))
  .filter(ee.Filter.lt('waterFrac', 0.98));

print('有效结果数量:', fc.size());

// ================= 图表 =================
print(ui.Chart.feature.byFeature(fc, 'date', ['area'])
  .setOptions({
    title: 'Three Gorges Gorge Section Water Area (Tuned)',
    pointSize: 2,
    hAxis: {title: 'Date'},
    vAxis: {title: 'Area (km²)'}
  }));

// ================= 导出 CSV =================
Export.table.toDrive({
  collection: fc,
  description: 'ThreeGorges_LocalGorge_Area_Tuned',
  fileFormat: 'CSV'
});

Answer 1

晚上好🌙🌙🌙
本答案参考通义千问

你提供的 GEE 代码在提取水体面积时出现了问题，主要可能有以下几个原因：

---

一、问题分析

1. 水体二值化阈值设置不准确

• MNDWI_THRESHOLD = -0.05 是一个较宽松的阈值，但在某些情况下（如高含沙量、植被干扰）可能导致误判。
• 重点问题：如果该阈值未正确反映水体与非水体的区分边界，会导致水体面积过大或过小。

2. 形态学处理不当

• morph() 函数中使用了 focal_min + focal_max 的组合，这会起到“开运算”的作用，但若半径太小（如 morphRadius = 1），可能无法有效去除噪点。
• 重点问题：形态学处理不足可能导致孤立噪点被保留，从而影响面积计算。

3. 连通域过滤参数设置不合理

• minSize = 50 可能太小，导致过多的噪声被保留。
• 重点问题：连通域过滤不够严格，可能将非水体区域误认为是水体。

4. 坡度掩膜应用顺序问题

• 当前代码中 slopeMask 在最后才被应用，可能已经错过了对水体图的早期筛选。
• 重点问题：坡度掩膜应尽早应用于图像，以减少后续处理的复杂性。

5. 面积计算方式错误

• 使用 reduceRegion(...).get('water') 时，可能会因为 roi 外部区域的像素未被裁剪而造成误差。
• 重点问题：reduceRegion 应配合 clip(roi) 来确保只计算 ROI 内部的水体面积。

---

二、解决方案

1. 调整 MNDWI 阈值

• 尝试更精确的阈值，例如：

  var MNDWI_THRESHOLD = -0.1; // 更严格的水体识别
  

2. 优化形态学处理

• 增加形态学核的大小，提高去噪效果：

  var morphRadius = 3; // 增大到 3（相当于 90m）
  

3. 提高连通域过滤阈值

• 将 minSize 设置为更大值，避免误判：

  var minSize = 100; // 增大到 100 像元
  

4. 提前应用坡度掩膜

• 将 slopeMask 应用于 col.map(prep) 之前，提高效率和准确性：

  var col = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
    .merge(ee.ImageCollection('LANDSAT/LE07/C02/T1_L2').filterDate('2000-01-01', '2003-05-31'))
    .merge(ee.ImageCollection('LANDSAT/LC08/C02/T1_L2'))
    .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
    .filterBounds(roi)
    .filterDate(startDate, endDate)
    .filter(ee.Filter.lt('CLOUD_COVER', MAX_CLOUD))
    .map(function(img) {
      return img.updateMask(slopeMask); // 提前加入坡度掩膜
    })
    .map(prep);
  

5. 确保面积计算在 ROI 内部

• 在 process(img) 中增加 img.clip(roi)：

  function process(img) {
    img = ee.Image(img).clip(roi);
  
    // ... 其他处理 ...
  }
  

---

三、修改后的代码片段

// ================= 参数 =================
var startDate = '2000-01-01';
var endDate   = '2025-12-31';

var MAX_CLOUD = 20;          // 云量阈值
var MIN_PIXELS = 300;        // 小区域最低有效像元数
var minSize = 100;           // 最小水体斑块（像元数），保留河道连续性
var MNDWI_THRESHOLD = -0.1;  // 更严格的浑水阈值
var morphRadius = 3;         // 圆形核半径（90m），增强去噪

// ================= DEM 坡度掩膜（25°去山体阴影） =================
var dem = ee.Image('USGS/SRTMGL1_003');
var slopeMask = ee.Terrain.slope(dem).lt(25);

// ================= 数据源（L7 仅取 SLC-off 前：截止 2003-05-31） =================
var col = ee.ImageCollection('LANDSAT/LT05/C02/T1_L2')
  .merge(ee.ImageCollection('LANDSAT/LE07/C02/T1_L2')
           .filterDate('2000-01-01', '2003-05-31'))
  .merge(ee.ImageCollection('LANDSAT/LC08/C02/T1_L2'))
  .merge(ee.ImageCollection('LANDSAT/LC09/C02/T1_L2'))
  .filterBounds(roi)
  .filterDate(startDate, endDate)
  .filter(ee.Filter.lt('CLOUD_COVER', MAX_CLOUD))
  .map(function(img) {
    return img.updateMask(slopeMask); // 提前加入坡度掩膜
  })
  .map(prep);

// ================= 单景处理 =================
function process(img) {
  img = ee.Image(img).clip(roi);

  // 水体二值图（浑水阈值 -0.1）
  var water = img.gt(MNDWI_THRESHOLD).rename('water');

  // 轻度形态学去噪
  var k = ee.Kernel.circle(morphRadius);
  water = water.focal_min({kernel: k, iterations: 1})
            .focal_max({kernel: k, iterations: 1});

  // 连通域过滤（maxSize=1024，阈值100）
  var conn = water.connectedPixelCount(1024, true);
  water = water.updateMask(conn.gte(minSize));

  // 面积 (km²)
  var area = water.multiply(ee.Image.pixelArea())
                  .reduceRegion({
                    reducer: ee.Reducer.sum(),
                    geometry: roi,
                    scale: 30,
                    maxPixels: 1e13,
                    tileScale: 4
                  }).get('water');
  var area_km2 = ee.Number(area).divide(1e6);

  // 水体比例
  var frac = water.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: roi,
    scale: 30,
    maxPixels: 1e13
  }).get('water');

  // 有效像元总数
  var pixels = img.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: roi,
    scale: 30,
    maxPixels: 1e13
  }).get('MNDWI');

  return ee.Feature(null, {
    'date': ee.Date(img.get('system:time_start')).format('YYYY-MM-dd'),
    'time': img.get('system:time_start'),
    'area': area_km2,
    'waterFrac': frac,
    'pixels': pixels
  });
}

---

四、建议进一步优化的方向

| 项目 | 说明 | |------|------| | 多时相融合 | 可以考虑对多个时间点的水体结果进行平均或合成，提高稳定性。 | | 引入 NDVI 或其他指标 | 结合 NDVI 等指标，进一步排除植被干扰。 | | 使用机器学习模型 | 如随机森林或 U-Net，提升水体识别精度。 |

---

如果你需要我帮你生成最终的水体面积图表或导出数据，请告诉我！