基于 GEE 开发的一种利用 OTSU 算法实现水体提取的便捷工具
目录
一、OTSU 算法（大津法）
二、五种水体指数
三、水体指数计算函数（原理 + 代码细节）
（一）归一化水体指数（NDWI）
（二）修改后的归一化水体指数（MNDWI）
（三）无阴影水体指数（AWEI_nsh）
（四）阴影水体指数（AWEI_sh）
（五）2015 版水体指数（WI_2015）
四、影像数据加载与预处理（流程 + 参数意义）
（一）ROI 初始化与地图配置
（二）Landsat 8 影像集合筛选（多条件过滤）
（三）影像合成与裁剪（数据融合 + 范围聚焦）
五、OTSU 阈值算法（数学原理 + 代码实现）
（一）输入数据结构
（二）OTSU 阈值算法核心代码
（三）关键数学公式
六、交互界面设计（UI 组件 + 事件绑定）
（一）面板容器初始化关键代码
（二）核心 UI 组件实现
（三）事件绑定逻辑
七、交互功能函数（业务逻辑 + 异常处理）
（一）水体提取函数（extractWater）—— 核心业务函数
（二）导出函数（exportWaterIndex + exportThresholdWater）
八、代码运行流程与注意事项
九、运行结果
若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！
本代码基于Google Earth Engine（GEE） 平台开发，依赖 GEE 的云原生地理空间计算能力和影像库支持。GEE 提供了 Landsat、Sentinel 等海量卫星影像数据，以及ee.Image、ee.ImageCollection等核心数据结构，支持通过 JavaScript API 实现遥感数据处理与分析。本代码的核心目标是解决水体信息自动化提取问题，适用于水资源调查、洪涝灾害监测、生态环境评估等场景。
一、OTSU 算法（大津法）
OTSU 算法（大津法）是一种确定图像二值化分割阈值的算法，由日本学者大津于1979年提出。从大津法的原理上来讲，该方法又称作最大类间方差法，因为按照大津法求得的阈值进行图像二值化分割后，前景与背景图像的类间方差最大。
它被认为是图像分割中阈值选取的最佳算法，计算简单，不受图像亮度和对比度的影响，因此在数字图像处理上得到了广泛的应用。它是按图像的灰度特性，将图像分成背景和前景两部分。因方差是灰度分布均匀性的一种度量，背景和前景之间的类间方差越大，说明构成图像的两部分的差别越大，当部分前景错分为背景或部分背景错分为前景都会导致两部分差别变小。因此，使类间方差最大的分割意味着错分概率最小。
- 应用：是求图像全局阈值的最佳方法，应用不言而喻，适用于大部分需要求图像全局阈值的场合。
- 优点：计算简单快速，不受图像亮度和对比度的影响。
- 缺点：对图像噪声敏感；只能针对单一目标分割；当目标和背景大小比例悬殊、类间方差函数可能呈现双峰或者多峰，这个时候效果不好。
二、五种水体指数
- NDWI（归一化水体指数）
- 计算公式：(Green - NIR) / (Green + NIR)
- 使用波段：Green = B3，NIR = B5
- 原理：利用水体在绿光波段反射率高、近红外波段反射率低的特性
- MNDWI（修改后的归一化水体指数）
- 计算公式：(Green - SWIR) / (Green + SWIR)
- 使用波段：Green = B3，SWIR = B6
- 优势：相比NDWI，能更好地区分水体与建筑物
- AWEI_nsh（无阴影水体指数）
- 计算公式：4*(Green - SWIR1) - (0.25NIR + 2.75SWIR2)
- 使用波段：Green=B3，NIR=B5，SWIR1=B6，SWIR2=B7
- 特点：减少阴影对水体提取的干扰
- AWEI_sh（阴影水体指数）
- 计算公式：Blue + 2.5Green - 1.5(NIR + SWIR1) - 0.25*SWIR2
- 使用波段：Blue=B2，Green=B3，NIR=B5，SWIR1=B6，SWIR2=B7
- 适用场景：适合提取存在阴影的水体区域
- WI_2015（2015版水体指数）
- 计算公式：1.7204 + 171Green + 3Red - 70NIR - 45SWIR1 - 71*SWIR2
- 使用波段：Red=B4，Green=B3，NIR=B5，SWIR1=B6，SWIR2=B7
- 优势：经过优化的经验公式，在多种场景下表现稳定
三、水体指数计算函数（原理 + 代码细节）
水体指数的本质是通过不同波段的辐射特性差异，增强水体与非水体（植被、土壤、建筑、云等）的灰度对比。代码中 5 种指数的设计逻辑和代码细节如下：
（一）归一化水体指数（NDWI）
- 设计原理：水体在绿色波段（B3，波长～0.56μm）反射较强，在近红外波段（B5，波长～0.86μm）因强吸收而反射较弱；而植被在近红外波段反射极强。通过
(Green - NIR)/(Green + NIR)的归一化计算，可使水体像元值为正，植被 / 土壤像元值为负，实现水体初步分离。 - 关键代码：
function NDWI(image) {// normalizedDifference：GEE内置函数，直接计算两波段归一化差值// addBands：将计算得到的NDWI作为新波段添加到原影像中，不覆盖原始波段return image.addBands(image.normalizedDifference(["B3", "B5"]).rename("NDWI") // rename：为新波段命名为"NDWI"，便于后续筛选调用); } - 局限性：在植被密集区域，若植被覆盖的水体边缘易被误判为非水体，需结合其他指数辅助验证。
（二）修改后的归一化水体指数（MNDWI）
- 优化逻辑：针对 NDWI 对植被干扰敏感的问题，将近红外波段（B5）替换为短波红外波段（B6，波长～1.61μm）。土壤和植被在短波红外波段反射较强，水体在该波段吸收更强，通过
(Green - SWIR)/(Green + SWIR)可进一步抑制植被 / 土壤干扰，更适合城镇水体、小面积水体提取。 - 代码差异：仅将
normalizedDifference的输入波段从["B3","B5"]改为["B3","B6"]，其余逻辑与 NDWI 一致，保证函数接口统一性。
（三）无阴影水体指数（AWEI_nsh）
- 设计背景：阴影（如建筑物阴影、山体阴影）的灰度值与水体接近，易导致 NDWI/MNDWI 误判。AWEI_nsh 通过多波段加权组合，强化水体与阴影的差异，公式为
4*(G - SWIR1) - (0.25*NIR + 2.75*SWIR2)。 - 关键代码：
function AWEI_nsh(image) {var awei_nsh = image.expression("4 * (G - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)", // 数学表达式{ // 波段映射：将表达式中的变量与影像实际波段关联"G": image.select("B3"), // 绿色波段"NIR": image.select("B5"), // 近红外波段"SWIR1": image.select("B6"),// 短波红外1波段"SWIR2": image.select("B7") // 短波红外2波段});return image.addBands(awei_nsh.rename("AWEI_nsh")); } - 权重意义：4 倍放大绿色与 SWIR1 的差异，同时用 0.25 和 2.75 的权重平衡 NIR 与 SWIR2 的抑制效果，确保无阴影水体的灰度值显著高于其他地物。
（四）阴影水体指数（AWEI_sh）
- 适用场景：针对山区、高层建筑区等存在大量阴影的区域，公式引入蓝色波段（B2，波长～0.48μm），通过
B + 2.5*G - 1.5*(NIR + SWIR1) - 0.25*SWIR2增强阴影区域水体的识别能力。 - 蓝色波段作用：阴影在蓝色波段的反射率低于水体，加入蓝色波段可进一步拉开阴影与水体的灰度差距，减少误判。
（五）2015 版水体指数（WI_2015）
- 技术特点：由学者提出的高精度指数，引入红色波段（B4，波长～0.65μm） 和固定常数项（1.7204），公式为
1.7204 + 171*G + 3*R - 70*NIR - 45*SWIR1 - 71*SWIR2。 - 优势：通过大权重（如 171 倍 G、71 倍 SWIR2）强化关键波段的区分能力，适用于复杂下垫面（如混合植被、浅水区、浑浊水体） 的提取，精度高于传统指数。
四、影像数据加载与预处理（流程 + 参数意义）
该模块是水体提取的 “数据准备阶段”，需确保输入影像的空间准确性、时间有效性和质量可靠性，具体流程如下：
（一）ROI 初始化与地图配置
- 关键代码：
var roi = geometry; // 需提前定义ROI（如ee.Geometry.Polygon） Map.centerObject(roi, 11); // 将地图中心定位到ROI，缩放级别11 Map.addLayer(roi, {color: "red"}, "roi", false); // 添加ROI边界图层，颜色红色，图层名"roi"，默认关闭（false） - 缩放级别意义：GEE 地图缩放级别 1-20，级别 11 对应地面分辨率约 100 米（因投影而异），既能完整显示 ROI，又不会因级别过高导致加载缓慢。
（二）Landsat 8 影像集合筛选（多条件过滤）
- 数据来源：
ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA")，其中 “C02” 表示 Landsat 8 的第二版数据（经过辐射定标、大气校正等预处理），“T1_TOA” 表示大气顶部反射率产品，可直接用于指数计算。 - 关键代码：
var l8Col = ee.ImageCollection("LANDSAT/LC08/C02/T1_TOA").filterBounds(roi) // 空间筛选：仅保留与ROI相交的影像.filterDate("2017-1-1", "2017-12-30") // 时间筛选：2017全年.map(ee.Algorithms.Landsat.simpleCloudScore) // 计算云评分.map(function(image) {// 云去除：保留云评分≤10的像元（0=无云，100=全云）return image.updateMask(image.select("cloud").lte(10));}).map(NDWI) // 依次添加5种水体指数波段.map(MNDWI).map(AWEI_nsh).map(AWEI_sh).map(WI_2015); - 云去除关键：
simpleCloudScore为 Landsat 影像添加 “cloud” 波段，updateMask函数根据云评分掩膜掉高云量像元（掩膜后的像元为null，不参与后续计算），避免云对水体指数的干扰。
（三）影像合成与裁剪（数据融合 + 范围聚焦）
- 中值合成：
l8Col.median()将 2017 年多幅影像按像元取中值，优势在于：- 抑制异常值（如单幅影像的噪声、云残留）；
- 生成全年 “无云” 效果的合成影像，代表年度平均下垫面状态。
- 裁剪操作：
clip(roi)将合成影像的范围限制为 ROI，减少无效区域计算量，同时确保后续提取结果仅聚焦目标区域。
五、OTSU 阈值算法（数学原理 + 代码实现）
Otsu 算法是自适应阈值分割的经典方法，无需人工设定阈值，核心是找到使 “前景（水体）与背景（非水体）类间方差最大” 的阈值，此时两类的区分度最高。
（一）输入数据结构
算法输入为reduceRegion计算得到的直方图字典，包含两个关键属性：
histogram：数组，每个元素代表对应灰度区间的像素数量；bucketMeans：数组，每个元素代表对应灰度区间的均值（即该区间的 “代表灰度值”）。
（二）OTSU 阈值算法核心代码
var otsu = function(histogram) {// 1. 将直方图属性转换为ee.Array，便于后续计算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); // 全局灰度均值// 2. 生成分割点索引（从1到size-1，避免全选/全不选）var indices = ee.List.sequence(1, size);// 3. 遍历每个分割点，计算类间方差（BSS）var bss = indices.map(function(i) {// 前景（假设分割点左侧为前景）：前i个区间的像素计数和均值var aCounts = counts.slice(0, 0, i); // 前i个区间的计数var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]); // 前景总像素数var aMeans = means.slice(0, 0, i); // 前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); // 背景均值// 计算类间方差：BSS = (aCount*(aMean-mean)²) + (bCount*(bMean-mean)²)return aCount.multiply(aMean.subtract(mean).pow(2)).add(bCount.multiply(bMean.subtract(mean).pow(2)));});// 4. 找到类间方差最大的分割点，返回对应区间均值（即最优阈值）return means.sort(bss).get([-1]); // 按BSS排序，取最后一个（最大BSS）对应的均值
};
（三）关键数学公式
- 全局均值：
mean = 总灰度和 / 总像素数 - 类间方差：
BSS = aCount*(aMean-mean)² + bCount*(bMean-mean)²（aCount= 前景像素数，aMean= 前景均值，bCount= 背景像素数，bMean= 背景均值） - 最优阈值：对应
BSS最大值的bucketMeans元素。
六、交互界面设计（UI 组件 + 事件绑定）
GEE 的ui模块支持自定义交互界面，代码设计的面板位于地图左上角，包含 “选择 - 提取 - 导出” 完整流程的控件，降低用户操作门槛。
（一）面板容器初始化关键代码
var waterPanel = ui.Panel({style: {position: 'top-left', padding: '8px', width: '200px'} // 位置：左上角，内边距8px，宽度200px（适配常用屏幕）
});
waterPanel.add(ui.Label("水体提取工具", {fontWeight: 'bold', fontSize: '16px'}));
// 添加标题，加粗+16号字体，提升辨识度
（二）核心 UI 组件实现
| 组件类型 | 功能描述 | 关键代码片段 |
|---|---|---|
| 下拉选择框 | 选择待使用的水体指数，默认 MNDWI | var waterIndexSelect = ui.Select({items: ["NDWI",...], value: "MNDWI", ...}) |
| 提取水体按钮 | 触发水体提取逻辑，调用extractWater函数 | var extractButton = ui.Button({label: "提取水体", onClick: extractWater}) |
| 导出按钮（2 个） | 分别导出水体指数连续影像和二值影像，调用对应导出函数 | var exportIndexButton = ui.Button({label: "导出水体指数图像", ...}) |
| 阈值显示标签 | 显示 Otsu 算法计算的最优阈值，初始为空，提取后更新 | var thresholdLabel = ui.Label("阈值：", {fontWeight: 'bold', ...}) |
（三）事件绑定逻辑
- 下拉框
onChange事件：选择指数后，打印当前选择，并更新阈值标签的提示文本（如 “阈值：MNDWI”），告知用户当前阈值对应的指数； - 按钮
onClick事件：直接绑定对应功能函数（如extractWater），点击即执行，无需额外触发条件。
七、交互功能函数（业务逻辑 + 异常处理）
（一）水体提取函数（extractWater）—— 核心业务函数
function extractWater() {// 1. 获取用户选择的指数，提取对应波段var indexName = waterIndexSelect.getValue();var indexImage = l8Image.select(indexName); // 筛选出选定的水体指数波段// 2. 计算指数波段在ROI内的直方图（用于Otsu算法）var histoDict = indexImage.reduceRegion({reducer: ee.Reducer.histogram(), // 直方图 reducergeometry: roi, // 计算范围：ROIscale: 30, // 分辨率：30米（Landsat 8原始分辨率）maxPixels: 1e13 // 最大像素数：避免因ROI过大导致计算失败});var histo = ee.Dictionary(histoDict).get(indexName); // 提取当前指数的直方图// 3. 调用Otsu算法计算最优阈值var threshold = otsu(histo);print("计算得到的 " + indexName + " 阈值：", threshold); // 打印阈值到控制台// 4. 移除旧的水体图层（避免图层堆积）var layers = Map.layers();for (var i = layers.length() - 1; i >= 0; i--) {var layer = layers.get(i);// 匹配图层名以"水体 ("开头的图层（如"水体 (MNDWI)"）if (layer.getName().indexOf("水体 (") === 0) {Map.layers().remove(layer);}}// 5. 生成水体二值图像（1=水体，掩膜=非水体）var waterMask = ee.Image(1) // 创建全1图像.updateMask(indexImage.gte(threshold)); // 掩膜：仅保留指数≥阈值的像元var waterLayerName = "水体 (" + indexName + ")";Map.addLayer(waterMask, {palette: "ff0000"}, waterLayerName, true); // 添加到地图，红色显示，默认开启// 6. 更新阈值标签（需用evaluate转换为客户端值，因GEE计算在服务端）threshold.evaluate(function(t) {thresholdLabel.setValue("阈值：" + t.toFixed(4)); // 保留4位小数，提升可读性});// 7. 保存当前二值图像，供导出使用currentWaterMask = waterMask;
}
reduceRegion的scale:30：与 Landsat 8 的原始分辨率一致，确保直方图计算的准确性；updateMask的作用：掩膜后的非水体像元不显示（透明），仅保留红色的水体像元，便于可视化；evaluate函数：GEE 中服务端数据（如threshold）需通过evaluate转换为客户端 JavaScript 值，才能更新 UI 标签。
（二）导出函数（exportWaterIndex + exportThresholdWater）
-
共性配置：
- 导出目标：Google Drive；
- 文件夹：
folder: "EarthEngineOutputs"（需用户提前在 Drive 创建该文件夹）； - 投影：
crs: "EPSG:4326"（WGS84 坐标系，全球通用）； - 分辨率：
scale:30（与原始影像一致）； - 最大像素数：
maxPixels:1e13（支持大区域导出，避免像素数超限报错）。
-
差异点：
exportWaterIndex：导出连续值影像（如 MNDWI 的灰度值影像），文件名格式为WaterIndex_MNDWI；exportThresholdWater：导出二值影像（1 = 水体，0 = 非水体，需先执行extractWater生成currentWaterMask），文件名格式为Threshold_WaterMask_MNDWI，并添加前置判断：if (currentWaterMask === null) {print("请先点击【提取水体】按钮获取二值图像！");return; // 未提取时终止导出，避免报错 }
八、代码运行流程与注意事项
完整运行流程：
- 准备工作：在 GEE 代码编辑器中定义
roi（如通过绘图工具绘制多边形）； - 加载数据：运行代码后，GEE 自动加载 2017 年 ROI 内的 Landsat 8 影像并预处理；
- 用户操作：
- 从下拉框选择水体指数（如 MNDWI）；
- 点击 “提取水体”，等待 Otsu 算法计算阈值并显示红色水体图层；
- 点击对应导出按钮，在 GEE 任务面板中启动导出，完成后在 Google Drive 查看结果。
关键注意事项：
- ROI 定义：
roi必须是ee.Geometry类型（如ee.Geometry.Polygon），否则会导致filterBounds和clip操作失败； - 数据权限：确保用户 GEE 账号有访问 Landsat 8 C02 数据的权限（默认开放）；
- 导出等待：大区域导出可能耗时较长（几分钟到几十分钟），需在 GEE “Tasks” 面板查看进度；
- 指数选择建议：
- 无阴影区域：优先选 MNDWI；
- 阴影密集区域：优先选 AWEI_sh；
- 复杂下垫面：优先选 WI_2015。
九、运行结果
若觉得代码对您的研究 / 项目有帮助，欢迎点击打赏支持！需要完整代码的朋友，打赏后可在后台私信（复制文章标题发给我），我会尽快发您完整可运行代码，感谢支持！