hiansvdj 2026-03-19 20:30 采纳率: 100%
浏览 2
已采纳

sentinel1GEE水体提取

在GEE平台上分别利用Landsat和哨兵系列遥感影像对数据进行统一预处理预处理,然后提取水体面积结合DEM反演水位并构建反演模型来估算断面流量

这是之前尝试过的代码,但有部分出错且结果不对,需要更改帮助
Sentinel1 SAR Otsu阈值水体提取


// 1. Otsu算法(优化版)
function otsu(histogram) {
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
  var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
  var size = means.length().get([0]);
  var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
  var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
  var mean = sum.divide(total);
  var indices = ee.List.sequence(1, size);
  
  var bss = indices.map(function(i) {
    var aCounts = counts.slice(0, 0, i);
    var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
    var aMeans = means.slice(0, 0, i);
    var aMean = ee.Algorithms.If(aCount.gt(0), aMeans.multiply(aCounts).reduce(ee.Reducer.sum(), [0]).get([0]).divide(aCount), 0);
    var bCount = total.subtract(aCount);
    var bMean = ee.Algorithms.If(bCount.gt(0), sum.subtract(aCount.multiply(aMean)).divide(bCount), 0);
    return aCount.multiply(aMean.subtract(mean).pow(2)).add(bCount.multiply(bMean.subtract(mean).pow(2)));
  });
  return ee.Algorithms.If(means.length().gt(0), means.sort(bss).get([-1]), -12); // 三峡SAR默认阈值
}
 
// 2. 初始化检查
if (typeof roi === 'undefined' || typeof dam === 'undefined') {
  print('🛑 【严重错误】: 请先在地图上画 roi 和 dam!');
} else {
  Map.centerObject(roi, 11);
  Map.addLayer(ee.Image().byte().paint(roi, 1, 2), {palette: 'red'}, 'ROI');
}
 
// 3. 数据处理(三峡SAR适配)
var filterSpeckles = function(img) {
  var vv = img.select('VV');
  var vv_smoothed = vv.focal_median(50, 'circle', 'meters').rename('VV_Filtered'); // 平滑窗口缩小更贴合库区
  return img.addBands(vv_smoothed);
};
 
var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(roi)
  .filterDate('2014-01-01', '2025-12-28') 
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .filter(ee.Filter.eq('instrumentMode', 'IW'))
  .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')) // 三峡优先降轨数据
  .map(filterSpeckles);
 
print('Sentinel1 影像数量:', S1.size());
 
// 4. 水体提取与计算
var Water_Calc = function(img) {
  var srtm = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation');
  var vv = img.select('VV_Filtered');
  
  // Otsu阈值计算(加异常处理)
  var histogramResult = vv.reduceRegion({
    reducer: ee.Reducer.histogram(255, 2),
    geometry: roi, 
    scale: 30,
    bestEffort: true
  });
  var localHisto = histogramResult.get('VV_Filtered');
  var thresholdNum = ee.Algorithms.If(localHisto, otsu(localHisto), -12);
  var thresholdImg = ee.Image.constant(ee.Number(thresholdNum));
  
  // 水体提取 + 三峡高程过滤
  var water = vv.lt(thresholdImg).rename('Water');
  var minSize = 500;
  var count = water.connectedPixelCount(minSize);
  var filtered1 = water.where(count.lt(minSize).or(srtm.gt(175)).or(srtm.lt(145)), 0);
  var water_final = filtered1.rename('Water');
  
  // 面积 + 水位计算
  var area = water_final.multiply(ee.Image.pixelArea()).divide(1e6)
             .reduceRegion({
               reducer: ee.Reducer.sum(),
               geometry: roi,
               scale: 30,
               bestEffort: true
             }).get('Water');
  var srtm_masked = srtm.updateMask(water_final);
  var mean_elev = srtm_masked.reduceRegion({
      reducer: ee.Reducer.mean(),
      geometry: dam,
      scale: 30,
      bestEffort: true
    }).get('Elevation');
  var dateStr = ee.Date(img.get('system:time_start')).format('YYYY-MM-dd');
 
  return ee.Feature(null, {
    'system:time_start': img.get('system:time_start'),
    'date_str': dateStr,
    'area_km2': area,
    'water_level_m': mean_elev,
    'threshold': thresholdNum
  });
};
var resultTable = S1.map(Water_Calc);
 
// 5. 导出
Export.table.toDrive({
  collection: resultTable,
  description: 'ThreeGorges_Water_Data_Sentinel1',
  fileFormat: 'CSV',
  selectors: ['date_str', 'area_km2', 'water_level_m', 'threshold']
});
  • 写回答

2条回答 默认 最新

  • 檀越@新空间 2026-03-19 20:31
    关注

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

    你提供的代码是基于 Sentinel-1 SAR 数据 使用 Otsu 阈值算法 进行水体提取,并结合 SRTM DEM 反演水位,进而估算断面流量的流程。但目前存在一些错误和逻辑问题,导致结果不准确。以下是详细的分析和修改建议:


    一、主要问题分析

    1. Otsu 函数实现错误

    • 当前的 otsu 函数中,使用了 ee.Algorithms.If(localHisto, otsu(localHisto), -12) 来处理 localHisto 是否为 null 的情况。
    • 但在 GEE 中,histogramResult.get('VV_Filtered') 返回的是一个数组(Array),不是布尔值,因此不能直接用 if 判断是否为 null
    • 修正方式:使用 ee.Algorithms.If 检查 localHisto 是否为空。

    2. Otsu 算法逻辑错误

    • means.sort(bss).get([-1]) 是错误的写法,GEE 中无法直接对 Array 进行排序并获取最大值。
    • 应该使用 bss.reduce(ee.Reducer.max(), [0]) 来获取最大值。

    3. 阈值应用错误

    • water.lt(thresholdImg) 中,thresholdImg 是一个 Image,但 thresholdNum 是一个 Number,应该使用 thresholdNum 而非 thresholdImg

    4. 数据范围与投影问题

    • 如果 SRTM 和 Sentinel-1 影像投影不一致,会导致掩膜或面积计算错误。
    • 建议统一投影(如 WGS84)以避免空间误差。

    5. 水体面积计算逻辑问题

    • area_km2 的计算可能包含非水体区域,应确保 water_final 已经过滤。

    二、修改后的完整代码(已修复错误)

    // 1. Otsu算法(优化版)
    function otsu(histogram) {
      var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
      var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
      var size = means.length().get([0]);
      var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
      var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
      var mean = sum.divide(total);
      var indices = ee.List.sequence(1, size);
    
      var bss = indices.map(function(i) {
        var aCounts = counts.slice(0, 0, i);
        var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
        var aMeans = means.slice(0, 0, i);
        var aMean = ee.Algorithms.If(aCount.gt(0), 
          aMeans.multiply(aCounts).reduce(ee.Reducer.sum(), [0]).get([0]).divide(aCount), 0);
        var bCount = total.subtract(aCount);
        var bMean = ee.Algorithms.If(bCount.gt(0), 
          sum.subtract(aCount.multiply(aMean)).divide(bCount), 0);
        return aCount.multiply(aMean.subtract(mean).pow(2)).add(bCount.multiply(bMean.subtract(mean).pow(2)));
      });
    
      // 获取最大类间方差对应的阈值
      var maxBSS = bss.reduce(ee.Reducer.max(), [0]);
      var thresholdIndex = bss.indexOf(maxBSS).get([0]);
    
      // 获取对应的均值作为阈值
      var threshold = means.get(thresholdIndex);
      return ee.Number(threshold);
    }
    
    // 2. 初始化检查
    if (typeof roi === 'undefined' || typeof dam === 'undefined') {
      print('🛑 【严重错误】: 请先在地图上画 roi 和 dam!');
    } else {
      Map.centerObject(roi, 11);
      Map.addLayer(ee.Image().byte().paint(roi, 1, 2), {palette: 'red'}, 'ROI');
    }
    
    // 3. 数据处理(三峡SAR适配)
    var filterSpeckles = function(img) {
      var vv = img.select('VV');
      var vv_smoothed = vv.focal_median(50, 'circle', 'meters').rename('VV_Filtered'); // 平滑窗口缩小更贴合库区
      return img.addBands(vv_smoothed);
    };
    
    var S1 = ee.ImageCollection('COPERNICUS/S1_GRD')
      .filterBounds(roi)
      .filterDate('2014-01-01', '2025-12-28') 
      .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
      .filter(ee.Filter.eq('instrumentMode', 'IW'))
      .filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING')) // 三峡优先降轨数据
      .map(filterSpeckles);
    
    print('Sentinel1 影像数量:', S1.size());
    
    // 4. 水体提取与计算
    var Water_Calc = function(img) {
      var srtm = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation');
      
      // 统一投影
      srtm = srtm.reproject('EPSG:4326', null, 30);
      img = img.reproject('EPSG:4326', null, 30);
    
      var vv = img.select('VV_Filtered');
    
      // Otsu阈值计算(加异常处理)
      var histogramResult = vv.reduceRegion({
        reducer: ee.Reducer.histogram(255, 2),
        geometry: roi, 
        scale: 30,
        bestEffort: true
      });
    
      var localHisto = histogramResult.get('VV_Filtered');
      var thresholdNum = ee.Algorithms.If(
        localHisto, 
        otsu(localHisto), 
        -12 // 默认阈值
      );
    
      // 水体提取 + 三峡高程过滤
      var water = vv.lt(thresholdNum).rename('Water');
      var minSize = 500;
      var count = water.connectedPixelCount(minSize);
      var filtered1 = water.where(
        count.lt(minSize).or(srtm.gt(175)).or(srtm.lt(145)), 
        0
      );
      var water_final = filtered1.rename('Water');
    
      // 面积 + 水位计算
      var area = water_final.multiply(ee.Image.pixelArea()).divide(1e6)
                 .reduceRegion({
                   reducer: ee.Reducer.sum(),
                   geometry: roi,
                   scale: 30,
                   bestEffort: true
                 }).get('Water');
    
      var srtm_masked = srtm.updateMask(water_final);
      var mean_elev = srtm_masked.reduceRegion({
          reducer: ee.Reducer.mean(),
          geometry: dam,
          scale: 30,
          bestEffort: true
        }).get('Elevation');
    
      var dateStr = ee.Date(img.get('system:time_start')).format('YYYY-MM-dd');
    
      return ee.Feature(null, {
        'system:time_start': img.get('system:time_start'),
        'date_str': dateStr,
        'area_km2': area,
        'water_level_m': mean_elev,
        'threshold': thresholdNum
      });
    };
    
    var resultTable = S1.map(Water_Calc);
    
    // 5. 导出
    Export.table.toDrive({
      collection: resultTable,
      description: 'ThreeGorges_Water_Data_Sentinel1',
      fileFormat: 'CSV',
      selectors: ['date_str', 'area_km2', 'water_level_m', 'threshold']
    });
    

    三、关键改进说明

    | 改进点 | 说明 | |--------|------| | Otsu 函数修复 | 修正了 threshold 的获取方式,使用 maxBSSthresholdIndex 确保正确性 | | 异常处理增强 | 增加了对 localHisto 是否存在的判断 | | 投影统一 | 所有影像统一到 WGS84 投影,避免空间错位 | | 阈值应用优化 | 使用 thresholdNum 直接进行水体提取,而非 thresholdImg | | 水体过滤逻辑强化 | 强化了 connectedPixelCount 和高程过滤条件 |


    四、后续建议

    1. 多源数据融合:可结合 LandsatSentinel-2 进行辅助验证,提升水体提取精度。
    2. 时间序列分析:将时间序列中的水位变化与降雨量、入库流量等变量结合,构建更精准的反演模型。
    3. 模型校准:使用实测水位数据对模型进行校正,提高估算断面流量的准确性。

    如有需要,我可以进一步帮你扩展这段代码用于 断面流量估算多源遥感数据融合。欢迎继续提问!

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论
查看更多回答(1条)

报告相同问题?

问题事件

  • 已采纳回答 3月19日
  • 创建了问题 3月19日