在GEE平台上分别利用Landsat和哨兵系列遥感影像对数据进行统一预处理预处理，然后提取水体面积结合DEM反演水位并构建反演模型来估算断面流量
这是之前尝试过的代码，但有部分出错且结果不对，需要更改帮助
用于Landsat5789影像集的NDVI水体范围 提取
// 1. 初始化检查
// ============================================================
if (typeof roi === 'undefined' || typeof dam === 'undefined') {
print('🛑 【严重错误】: 请先在地图上画 roi 和 dam！');
print('👉 1. 画一个多边形包住整个水库，重命名为 roi');
print('👉 2. 画一个多边形在大坝深水区，重命名为 dam');
} else {
Map.centerObject(roi, 11);
Map.setOptions("TERRAIN");
}
// ============================================================
// 2. 参数设置
// ============================================================
var START_DATE = '1987-01-01';
var END_DATE = '2025-12-31';
var MAX_CLOUD_COVER = 30;
// ============================================================
// 3. 核心函数：去云与水体提取
// ============================================================
// Landsat 去云函数 (使用 QA_PIXEL 波段)
function maskLsr(image) {
var qa = image.select('QA_PIXEL');
// 位操作：云(Bit 3) 和 云影(Bit 4)
var mask = qa.bitwiseAnd(1 << 3).eq(0)
.and(qa.bitwiseAnd(1 << 4).eq(0));
return image.updateMask(mask)
.copyProperties(image, ["system:time_start"]);
}
// 核心水体提取函数 (NDWI)
var waterfunction = function(image){
// Landsat 8: B3=Green, B5=NIR；Landsat5/7:B2=green,B4=nir
var ndwi = image.normalizedDifference(['B3', 'B5']).rename('NDWI');
// 阈值判断：NDWI > 0 判定为水
var water = ndwi.gt(0.0).rename('Water');
return image.addBands(water).addBands(ndwi);
};
// ============================================================
// 4. 数据加载与处理
// ============================================================
// 合并 Landsat 8、7、5 数据集
var l5 = ee.ImageCollection('LANDSAT/LC05/C02/T1_TOA');
var l7 = ee.ImageCollection('LANDSAT/LC07/C02/T1_TOA');
var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_TOA');
var landsat = l8.merge(l7).merge(l5)
.filterBounds(roi)
.filterDate(START_DATE, END_DATE)
.filter(ee.Filter.lt('CLOUD_COVER', MAX_CLOUD_COVER))
.map(function(image){ return image.clip(roi); })
.map(maskLsr) // 1. 去云
.map(waterfunction); // 2. 算水体
print('筛选到的 Landsat 影像数量:', landsat.size());
// ============================================================
// 5. 数据清洗 (去除小斑块 & 错误高程)
// ============================================================
var Water_Clean = function(img) {
var water = img.select('Water');
// 使用 SRTM 高程
var alos = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation');
// 逻辑：去除小于500像素的小碎块
var minSize = 500;
var count = water.connectedPixelCount(minSize);
// 逻辑：去除高程异常区 (三峡水库水面通常 < 175m)
var filtered = water.where(count.lt(minSize).and(alos.gt(175)), 0);
return img.addBands(filtered.rename('Water_Clean'), null, true);
};
landsat = landsat.map(Water_Clean);
// ============================================================
// 6. 面积与水位计算
// ============================================================
var Water_Calc = function(img) {
var water = img.select('Water_Clean'); // 用清洗后的水体
var alos = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation');
// 1. 计算面积 (km²)
var area = water.multiply(ee.Image.pixelArea()).divide(1e6)
.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: roi,
scale: 30,
bestEffort: true,
maxPixels: 1e13
}).get('Water_Clean');
// 2. 计算水位 (m)
var alos_masked = alos.updateMask(water);
var max_elev = alos_masked.reduceRegion({
reducer: ee.Reducer.max(),
geometry: dam,
scale: 30,
bestEffort: true,
maxPixels: 1e13
}).get('Elevation');
// 防止无效值报错
var water_level = ee.Algorithms.If(max_elev, max_elev, 0);
// 3. 计算水深 (三峡库底约 90m)
var water_depth = ee.Number(water_level).subtract(90);
// 返回特征对象
return ee.Feature(null, {
'system:time_start': img.get('system:time_start'),
'area_km2': area,
'water_level_m': water_level,
'water_depth_m': water_depth
});
};
// 应用 Water_Calc 函数并收集结果
var waterStats = landsat.map(Water_Calc).toList(landsat.size());
// 打印结果
print('水体统计信息：', waterStats);
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('2006-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']
});
Sentinel2 NDWI+NDVI 水体提取
// 1. 初始化检查
if (typeof roi === 'undefined' || typeof dam === 'undefined') {
print('🛑 【严重错误】: 请先在地图上画 roi 和 dam！');
} else {
Map.centerObject(roi, 11);
Map.setOptions("TERRAIN");
}
// 2. 参数设置
var START_DATE = '2015-01-01';
var END_DATE = '2025-12-29';
var MAX_CLOUD_PERCENT = 20; // 改MIN为MAX，语义更对
// 3. 核心函数（优化阈值+去云）
function maskS2clouds(image) {
var qa = image.select('QA60');
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask).divide(1e4)
.select("B.*")
.copyProperties(image, ["system:time_start"]);
}
var addNDWI = function(image) {
return image.addBands(image.normalizedDifference(['B3', 'B8']).rename('NDWI'));
};
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']).rename('NDVI'));
};
// 优化组合掩膜：NDWI≥0.1 + NDVI≤0.2（更贴合三峡库区）
var CombinedMask = function(image) {
var ndwi = image.select('NDWI');
var ndvi = image.select('NDVI');
var waterMask = ndwi.gte(0.1).and(ndvi.lte(0.2)).rename('Water');
return image.addBands(waterMask);
};
// 4. 数据加载与按日合成
var s2 = ee.ImageCollection('COPERNICUS/S2_HARMONIZED')
.filterBounds(roi)
.filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', MAX_CLOUD_PERCENT))
.filterDate(START_DATE, END_DATE)
.map(function(image) { return image.clip(roi); });
function mosaicBy(imcol){
var imlist = imcol.toList(imcol.size());
var all_dates = imlist.map(function(im){ return ee.Image(im).date().format("YYYY-MM-dd"); });
var unique_dates = all_dates.distinct();
var mosaic_imlist = unique_dates.map(function(d){
var date1 = ee.Date(d);
var im = imcol
.filterDate(date1, date1.advance(1, "day"))
.mosaic();
return im.set(
"system:time_start", date1.millis(),
"date_str", d
);
});
return ee.ImageCollection(mosaic_imlist);
}
var s2day = mosaicBy(s2);
var S2_Processed = s2day
.map(maskS2clouds)
.map(addNDWI)
.map(addNDVI)
.map(CombinedMask)
.map(function(img) {
var count = img.select('Water').reduceRegion({
reducer: ee.Reducer.count(),
geometry: roi,
scale: 200,
maxPixels: 1e13
}).get('Water');
return img.set('valid_pixels', count);
})
.filter(ee.Filter.gt('valid_pixels', 100));
print('Sentinel2 有效影像数量:', S2_Processed.size());
// 5. 面积与水位计算
var Water_Calc = function(img) {
var water = img.select('Water');
var srtm = ee.Image("USGS/SRTMGL1_003").clip(roi).rename('Elevation');
var area = water.multiply(ee.Image.pixelArea()).divide(1e6)
.reduceRegion({
reducer: ee.Reducer.sum(),
geometry: roi,
scale: 30,
bestEffort: true,
maxPixels: 1e13
}).get('Water');
var srtm_masked = srtm.updateMask(water).updateMask(srtm.gte(145).and(srtm.lte(175)));
var mean_elev = srtm_masked.reduceRegion({
reducer: ee.Reducer.mean(),
geometry: dam,
scale: 30,
bestEffort: true,
maxPixels: 1e13
}).get('Elevation');
var water_level = ee.Algorithms.If(mean_elev, mean_elev, 0);
var dateStr = img.get('date_str');
return ee.Feature(null, {
'date': dateStr,
'area_km2': area,
'water_level_m': water_level
});
};
var resultTable = S2_Processed.map(Water_Calc);
// 6. 导出（修复语法错误）
Export.table.toDrive({
collection: resultTable,
description: 'ThreeGorges_Water_Data_Sentinel2',
fileFormat: 'CSV',
selectors: ['date', 'area_km2', 'water_level_m']
});
