作者：CSDN @ _养乐多_
本文介绍了通过Google Earth Engine平台，并使用哨兵数据提取水体掩膜的方法和代码。通过裁剪和去除云等处理步骤，最终得到具有水体掩膜的影像，并进行可视化和导出。这种方法基于归一化水体指数（NDWI）和OTSU阈值计算技术，无需复杂的图像处理算法，适用于快速获取水体信息的需求。
 图1 彩色合成
|
 图2 水体掩膜
|
一、OTSU阈值方法
1.1 原理
OTSU阈值方法是一种常用的图像分割算法，用于将图像分成背景和目标两个部分。在提取水体的应用中，可以使用OTSU阈值方法来自动分割水体和非水体区域。
OTSU阈值方法基于图像的灰度直方图，通过寻找最佳的阈值来实现分割。以下是OTSU阈值方法的原理步骤：
-
计算图像的灰度直方图：首先，将图像转换为灰度图像，并统计每个灰度级别的像素数量。这样可以得到图像的灰度直方图。
-
计算像素的总数：统计图像中的总像素数。
-
计算类间方差：选择一个阈值T，将图像分成两个类别：低于阈值T的像素类别为背景，高于阈值T的像素类别为目标（水体）。然后，计算两个类别的权重、均值和方差，并使用这些值来计算类间方差。类间方差是一个衡量分割质量的指标，当类间方差最大时，分割效果最佳。
-
寻找最大类间方差：通过遍历所有选择的阈值T，计算类间方差，并找到使类间方差最大的阈值T*。
-
应用阈值分割：将图像的每个像素与阈值T*进行比较，将像素分为水体和非水体两个部分。
-
OTSU阈值方法的原理是基于类间方差最大化的思想。通过最大化类间方差，可以实现最佳的分割效果，将水体和非水体区域有效地分开。这种方法不需要预先设定阈值，而是通过自动计算得到最优的阈值，因此具有较好的适应性和稳定性。
在应用中，可以将OTSU阈值方法应用于水体提取任务，通过对输入图像进行OTSU阈值分割，得到水体和非水体区域的二值图像，从而实现水体的提取和分割。
1.2 直方图
将像素直方图和类内方差最大的像素值打印出来，如下图所示。
用于演示阈值的示例链接（将像素直方图和类内方差最大的像素值打印了出来）：https://code.earthengine.google.com/4396f8dd90c07f824f132c568e350775?noload=true
二、代码详解
var roi = table
Map.centerObject(roi,10)
// 可视化参数
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
'74A901', '66A000', '529400', '3E8601', '207401', '056201',
'004C00', '023B01', '012E01', '011D01', '011301'];
//可视化参数，按843波段合成
var rgbVis = {
min: 0.0,
max: 0.35,
bands: ['B8', 'B4', 'B3'],
};
//按矢量边界裁剪
function roiClip(image){
return image.clip(roi)
}
//S2_SR去云
function remove_cloud(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(10000).select("B.*").copyProperties(image, ["system:time_start"])
}
//时间范围
var startDate = ee.Date('2020-01-01');
var endDate = ee.Date('2020-02-01');
//数据集选择
var S2 = ee.ImageCollection("COPERNICUS/S2_SR")
.filterDate(startDate, endDate)
.select('B3', 'B4', 'B8', 'B11', 'QA60')
.filterBounds(roi)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 5))
.map(remove_cloud)
.map(roiClip)
//运行同一天影像
var newImgS2 = S2.max();
print(newImgS2)
Map.addLayer(newImgS2, rgbVis, 'newImg_RGB');
//计算NDWI
function calNDWI(image)
{
var NDWI = image.normalizedDifference(['B3', 'B11']).rename('NDWI');
NDWI = NDWI.updateMask(NDWI.gt(-1).and(NDWI.lt(1)));
return image.addBands(NDWI);
}
//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 = aMeans.multiply(aCounts)
.reduce(ee.Reducer.sum(), [0]).get([0])
.divide(aCount);
var bCount = total.subtract(aCount);
var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
bCount.multiply(bMean.subtract(mean).pow(2)));
});
return means.sort(bss).get([-1]);
}
//计算水体掩膜
function calWaterMask(image){
var MNDWI_S2 = image.select("NDWI");
var histogram = MNDWI_S2.reduceRegion({
reducer: ee.Reducer.histogram(),
geometry: roi,
scale: 30,
maxPixels: 1e13,
tileScale: 8
});
var threshold_S2 = otsu(histogram.get("NDWI"));
var mask = MNDWI_S2.gte(threshold_S2);
return image.addBands(mask.rename('mask'));
}
//运行水体掩模函数
newImgS2 = calNDWI(newImgS2)
newImgS2 = calWaterMask(newImgS2).select('mask')
print(newImgS2)
Map.addLayer(newImgS2, {min:0, max:1, palette:['#DDDDDD', '#0099FF']}, 'S2_Water_body_mask');
//下载非水体概率影像
Export.image.toDrive({
image:newImgS2,
description: 'S2_water_mask',
scale:30,
region:roi,
fileFormat: 'GeoTIFF',
maxPixels:1e13,
});
这段代码的目的是对Sentinel-2遥感影像进行处理，提取水体掩膜，并将结果可视化和导出。
-
首先，定义了一个变量roi，表示感兴趣区域（Region of Interest），这是一个表格（table）对象。然后使用Map.centerObject(roi,10)将该区域设置为地图的中心，并指定缩放级别为10。
-
接下来，定义了一个调色板（palette）数组，用于可视化图像时的颜色映射。
-
然后，定义了一个可视化参数rgbVis，设置了最小值（min）和最大值（max），以及要显示的波段（bands）。
-
之后，定义了一个函数roiClip，用于将图像按照矢量边界进行裁剪。
-
接着是一个函数remove_cloud，用于去除图像中的云和雾。函数内部使用了位掩码（bitmask）的方式来选择需要保留的像素，然后将图像的像素值除以10000并选择特定波段，最后将处理后的图像的时间属性复制到结果中。
-
然后定义了起始日期和结束日期，用于选择要处理的影像时间范围。
-
接下来，使用ee.ImageCollection函数选择了名为"COPERNICUS/S2_SR"的Sentinel-2影像集合，并通过一系列操作对影像进行筛选、裁剪和去除云，最后将结果映射到ROI区域，并存储在变量S2中。
-
之后，使用S2.max()函数选取ROI区域内同一天的影像，并将结果存储在变量newImgS2中。
-
接着，定义了一个名为calNDWI的函数，用于计算归一化水体指数（NDWI）。
-
然后，定义了一个名为otsu的函数，用于计算阈值，该阈值用于将图像的像素分为水体和非水体。
-
接下来，定义了一个名为calWaterMask的函数，该函数通过计算NDWI的直方图，并使用OTSU方法计算阈值来生成水体掩膜。
-
然后，分别调用calNDWI和calWaterMask函数对newImgS2进行处理，并将结果可视化。
-
最后，使用Export.image.toDrive函数将非水体概率影像导出到Google Drive中。
三、代码链接
https://code.earthengine.google.com/f145e92abf2d0196aba61655745bad9c?noload=true
四、OTSU算法详解
// OTSU计算阈值
// 此函数接受一个直方图作为输入，并使用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);
// 创建从1到size的索引列表
var indices = ee.List.sequence(1, size);
// 计算类间方差（between-class variance）
var bss = indices.map(function(i) {
// 将counts数组切片，获取前i个桶的像素计数
var aCounts = counts.slice(0, 0, i);
// 计算前i个桶的像素计数之和
var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
// 将means数组切片，获取前i个桶的桶均值
var aMeans = means.slice(0, 0, i);
// 计算前i个桶的桶均值乘以像素计数的总和，然后除以像素计数之和
var aMean = aMeans.multiply(aCounts)
.reduce(ee.Reducer.sum(), [0]).get([0])
.divide(aCount);
// 计算后i个桶的像素计数之和
var bCount = total.subtract(aCount);
// 计算后i个桶的桶均值，公式为：(总和 - 前i个桶的像素计数之和 * 前i个桶的桶均值) / 后i个桶的像素计数之和
var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
// 计算类间方差
return aCount.multiply(aMean.subtract(mean).pow(2)).add(
bCount.multiply(bMean.subtract(mean).pow(2)));
});
// 将类间方差数组按照从小到大排序，获取最后一个（最大值），即最佳阈值
return means.sort(bss).get([-1]);
}
OTSU算法是一种自适应的图像阈值分割方法，它通过最大类间方差来选择最佳阈值。
具体步骤如下：
- 计算图像的灰度直方图。
- 计算总像素数和每个像素值的占比。
- 对每个像素值，计算类内方差，即以该像素值作为阈值进行分割时，前景和背景之间的方差。
- 选择使类内方差最大的像素值作为阈值。
- 利用该阈值将图像分割为前景和背景两部分。
这段代码的作用是实现OTSU算法的第4步，即计算类间方差并选择最佳阈值。
扫一扫
638
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
