在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']
});