3、水体指数法提取：三大指数的原理与计算
水体指数法，说白了就是利用不同波段的光谱差异，把水体和背景分开。我做了这么多年遥感，这招是最快、最稳的入门方法。今天咱们就聊聊三个经典指数：NDWI、MNDWI 和 AWEI。
3.1 归一化差异水体指数（NDWI）
NDWI 是 McFeeters 在 1996 年提出的。它的核心思路很简单：水体在绿光波段反射率高，在近红外波段反射率低。利用这个反差，就能把水凸显出来。
公式：
NDWI = (Green - NIR) / (Green + NIR)其中 Green 是绿光波段（比如 Landsat 的 Band 3），NIR 是近红外波段（比如 Band 5）。
取值逻辑：
- NDWI > 0：大概率是水体
- NDWI ≤ 0：非水体（植被、土壤、建筑等）
我个人习惯：NDWI 对纯净水体效果很好，但遇到建筑物阴影就容易翻车。阴影在近红外波段也很暗，NDWI 值会偏高，容易被误判为水。
避坑指南：我曾经在城区水体提取时，NDWI 把大片高楼阴影都识别成了水。后来加了坡度数据做掩膜，才把误判压下去。如果你做城市水体，记得留个心眼。
3.2 改进的归一化差异水体指数（MNDWI）
NDWI 的痛点在于对建筑和阴影敏感。徐涵秋老师在 2005 年做了改进——把近红外波段换成中红外波段（SWIR1）。
公式：
MNDWI = (Green - SWIR1) / (Green + SWIR1)SWIR1 对应 Landsat 的 Band 6（OLI 传感器）或 Band 5（TM/ETM+）。
为什么有效？
建筑物在 SWIR1 波段的反射率比 NIR 高得多。用 SWIR1 替换后，建筑和阴影的 MNDWI 值会明显降低，不容易和水体混淆。
| 地物类型 | NDWI 表现 | MNDWI 表现 |
|---|---|---|
| 纯净水体 | 高（>0） | 高（>0） |
| 建筑物 | 中等（可能>0） | 低（<0） |
| 植被 | 低（<0） | 低（<0） |
| 阴影 | 中等（可能>0） | 低（<0） |
你看，MNDWI 对建筑和阴影的抑制效果明显更好。我在项目中做过对比，MNDWI 的总体精度能比 NDWI 提高 5-10 个百分点。
3.3 自动化水体提取指数（AWEI）
AWEI 是 Feyisa 等人在 2014 年提出的。它有两个版本：AWEInsh（无阴影版）和 AWEIsh（阴影版）。
AWEInsh 公式：
AWEI_nsh = 4 × (Green - SWIR1) - (0.25 × NIR + 2.75 × SWIR2)AWEIsh 公式：
AWEI_sh = Blue + 2.5 × Green - 1.5 × (NIR + SWIR1) - 0.25 × SWIR2其中 SWIR2 对应 Landsat 的 Band 7。
核心思路：AWEI 通过多波段线性组合，把水体和背景的差异拉到最大。说白了就是给不同波段分配权重，让水体值更高、背景值更低。
怎么选版本？
- AWEInsh：适合阴影少的区域（比如开阔水域、农田区）
- AWEIsh：适合阴影多的区域（比如山区、城市）
注意：AWEI 的阈值不是固定的 0。你需要根据影像直方图手动调整。我一般先算一遍，然后看直方图的双峰分布，取谷值作为阈值。
3.4 三个指数的对比与选择
我整理了一张对比表，方便你快速决策：
| 指数 | 优势 | 劣势 | 适用场景 |
|---|---|---|---|
| NDWI | 计算简单，速度快 | 易受建筑和阴影干扰 | 纯净水体、大范围水域 |
| MNDWI | 抑制建筑和阴影效果好 | 对浑浊水体敏感度略降 | 城市水体、复杂背景 |
| AWEI | 精度高，适应性强 | 需要调阈值，计算稍复杂 | 高精度需求、阴影密集区 |
嗯，这里要注意：没有万能指数。我个人的工作流是先用 MNDWI 快速跑一遍，如果精度不够，再上 AWEI 做精细化处理。
3.5 核心逻辑框架图
下面这张图展示了三个指数的核心逻辑和适用场景：
3.6 代码实现（Python + GDAL）
下面是我常用的计算函数，直接复制就能用：
import numpy as np
from osgeo import gdal
def calc_ndwi(green_path, nir_path, output_path):
"""计算 NDWI"""
green_ds = gdal.Open(green_path)
nir_ds = gdal.Open(nir_path)
green = green_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
nir = nir_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
# 避免除零
ndwi = (green - nir) / (green + nir + 1e-10)
# 保存结果
driver = gdal.GetDriverByName('GTiff')
out_ds = driver.Create(output_path, green_ds.RasterXSize,
green_ds.RasterYSize, 1, gdal.GDT_Float32)
out_ds.SetGeoTransform(green_ds.GetGeoTransform())
out_ds.SetProjection(green_ds.GetProjection())
out_ds.GetRasterBand(1).WriteArray(ndwi)
out_ds.FlushCache()
print(f"NDWI 计算完成，结果保存至：{output_path}")
def calc_mndwi(green_path, swir1_path, output_path):
"""计算 MNDWI"""
green_ds = gdal.Open(green_path)
swir1_ds = gdal.Open(swir1_path)
green = green_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
swir1 = swir1_ds.GetRasterBand(1).ReadAsArray().astype(np.float32)
mndwi = (green - swir1) / (green + swir1 + 1e-10)
# 保存逻辑同上...
print(f"MNDWI 计算完成")
def calc_awei_nsh(green_path, nir_path, swir1_path, swir2_path, output_path):
"""计算 AWEI_nsh"""
# 读取各波段...
awei_nsh = 4 * (green - swir1) - (0.25 * nir + 2.75 * swir2)
# 保存逻辑同上...
print(f"AWEI_nsh 计算完成")小技巧：计算前记得检查波段顺序。我遇到过几次因为波段顺序搞反，算出来的指数全是负值。建议先打印一下波段信息确认。
好了，三个指数的原理和计算就聊到这儿。你想想看，从 NDWI 到 MNDWI 再到 AWEI，其实就是一个不断「打补丁」的过程——发现一个问题，改进一个版本。做工程就是这样，没有完美的算法，只有最适合当前场景的方案。
公众号:蓝海资料掘金营，微信deep3321