NDWI/MNDWI 水体提取全流程：指数计算、阈值确定到导出矢量

一、引言：为什么水体提取值得反复打磨

水体提取看似是遥感入门级的操作，但真正做过业务化生产的人都知道，它远不止“算个指数、卡个阈值”这么简单。同一景影像，换一个阈值，湖泊面积可能相差百分之几到百分之十几；同一片水域，换一种后处理策略，矢量边界可能从锯齿状变成平滑曲线，也可能把细小的坑塘全部吃掉。这些差异在科研中影响结论的稳健性，在工程中直接影响面积统计、淹没范围制图和变化检测的可信度。

NDWI（Normalized Difference Water Index）由 McFeeters 于 1996 年提出，利用绿波段与近红外波段的归一化差异来增强水体信号。随后 Xu 于 2006 年提出 MNDWI（Modified NDWI），用短波红外替换近红外，显著改善了对建筑物阴影和土壤背景的抑制能力。这两个指数至今仍是水体提取的“基础设施”。本文评述：指数法的生命力恰恰来自它的“笨”——波段运算透明、可解释、可复现，不依赖训练样本，在数据匮乏或算力受限的场景下依然是首选。

但“笨”不等于“简单”。本文试图确立一条贯穿全文的分析主线：水体提取的精度瓶颈，往往不在指数公式本身，而在阈值决策与后处理链路的工程细节上。围绕这条主线，下文将从物理基础讲到矢量导出，每一步都给出可操作的路径和可检验的判断依据。

二、水体提取的物理基础与指数家族谱系

2.1 水体在光学波段的响应特征

清洁水体在可见光蓝绿波段反射率较高，随波长增加反射率迅速下降，在近红外（NIR）和短波红外（SWIR）波段几乎全吸收，反射率接近零。这一“绿高红低、近红外趋零”的光谱形态，是所有水体指数赖以成立的物理基础。浑浊水体因悬浮泥沙散射，在红光和近红外波段的反射率会抬升，导致指数值下降；富营养化水体因藻类叶绿素在近红外的强反射，甚至可能出现“水华”像元被误判为植被的情况。本文评述：任何水体指数都不是“水体探测器”，而是“绿-近红外反差放大器”，理解这一点才能解释为什么阈值需要随水体类型调整。

2.2 指数家族谱系

从 NDWI 出发，学界衍生出一系列变体，核心思路都是“选一对波段，让水体在一端高、另一端低”。下表梳理了主要成员及其设计意图。

需要特别提醒的是，Gao 提出的 NDWI 与 McFeeters 的 NDWI 同名但物理含义完全不同，前者面向植被冠层含水量，后者面向开放水体。文献检索时若不加区分，极易张冠李戴。笔者认为，命名冲突是遥感指数领域的历史遗留问题，工程实践中应始终以公式和波段定义为准，而非以缩写为准。

三、NDWI 与 MNDWI 的公式推导与波段选择

3.1 归一化差异的数学本质

归一化差异指数的一般形式为 (A−B)/(A+B)，其值域为 [−1, 1]。当 A 远大于 B 时趋近 +1，反之趋近 −1。这种结构的优势在于：其一，比值形式对乘性光照差异不敏感，能部分抵消地形阴影和太阳高度角变化的影响；其二，值域有界，便于跨影像比较和阈值设定。本文评述：归一化差异本质是一种“对比度拉伸”，它把两个波段的绝对差异压缩到有界区间，代价是丢失了绝对反射率信息——这也是为什么单靠指数无法区分“深水”和“浅水”。

3.2 波段选择对照表

不同传感器的波段编号和波长范围不同，直接套用“第几波段”极易出错。下表给出主流传感器的对应关系。

Sentinel-2 的 B8 带宽较宽（0.78–0.90 μm），包含了部分水汽吸收带，在湿润地区可能引入轻微噪声；若追求更纯净的水体信号，可考虑 B8A（0.86–0.88 μm，窄近红外）。本文评述：波段选择不是“越窄越好”，窄波段信噪比通常更低，实际项目中应以信噪比和空间分辨率的平衡为准。

四、数据准备：传感器、波段与预处理

4.1 数据源选择

当前免费可用的中分辨率光学数据主要有 Sentinel-2（10–20 m，5 天重访）和 Landsat 8/9（30 m，16 天重访，双星联合 8 天）。Sentinel-2 在空间分辨率和重访周期上占优，适合中小水体和动态监测；Landsat 系列时间跨度长（可追溯至 1984 年），适合长时序分析。若需更高分辨率，可考虑 GF-1/2、PlanetScope 等商业或半商业数据，但需注意波段配置是否包含 SWIR。

4.2 预处理链路

标准预处理包括辐射定标、大气校正、云掩膜和地形校正。对于水体指数计算，大气校正的影响相对温和，因为归一化差异在一定程度上抵消了加性大气程辐射。但在浑浊水体或气溶胶光学厚度较大的区域，未校正的影像会导致指数值系统性偏移。本文评述：如果只是做相对变化检测，TOA 反射率往往够用；如果要做绝对阈值分割或跨时相对比，建议使用地表反射率产品（如 Sentinel-2 L2A、Landsat C2 L2）。

云掩膜是容易被忽视的环节。Sentinel-2 L2A 自带 Scene Classification Layer（SCL），可直接提取云、云影、雪等类别；Landsat C2 提供 QA_PIXEL 波段。工程实践中，建议将云影也纳入掩膜——云影在绿波段偏暗、在近红外也偏暗，其 NDWI 值可能落入水体区间，造成“影子水体”误判。

五、指数计算的工程实现

5.1 分块计算与内存管理

一景 Sentinel-2 L2A 的 10 m 波段约 1 亿像元，若同时读入多个波段做浮点运算，内存占用可达数 GB。工程上推荐两种策略：一是使用 rasterio 的窗口读取（windowed reading），按块处理；二是使用 xarray + dask 构建惰性计算图，最后一次性写出。本文评述：对于单景影像，窗口读取足够；对于时序分析，dask 的惰性求值能显著降低峰值内存。

import rasterio
import numpy as np

def compute_mndwi(green_path, swir_path, out_path, block_size=1024):
    with rasterio.open(green_path) as src_g, rasterio.open(swir_path) as src_s:
        profile = src_g.profile.copy()
        profile.update(dtype='float32', count=1, compress='lzw')
        with rasterio.open(out_path, 'w', **profile) as dst:
            for ji, window in src_g.block_windows(1):
                green = src_g.read(1, window=window).astype('float32')
                swir = src_s.read(1, window=window).astype('float32')
                denom = green + swir
                # 避免除零
                mndwi = np.where(denom == 0, 0, (green - swir) / denom)
                dst.write(mndwi, 1, window=window)

5.2 无效值处理

原始影像中常存在 0 值填充、云掩膜后的 NoData 区域。若直接参与运算，会出现 0/0 或异常大值。建议在计算前统一将 NoData 设为 NaN，并在输出时保留 NoData 标记。本文评述：很多“指数异常”问题，根源其实是 NoData 处理不当，而非公式错误。

六、阈值确定：从经验值到自适应算法

阈值确定是整条链路中主观性最强、对结果影响最大的环节。本节按“由简到繁”的顺序梳理四类策略。

6.1 固定阈值法

最经典的做法是取 NDWI > 0 或 MNDWI > 0 作为水体判据。这一经验值在多数场景下能给出合理结果，但在以下情况会失效：浑浊水体指数值偏低，可能被漏分；山地阴影、深色屋顶、沥青路面指数值偏高，可能被误分。本文评述：固定阈值适合快速预览和教学演示，不建议直接用于业务化面积统计。

6.2 Otsu 自适应阈值

Otsu 法通过最大化类间方差自动寻找分割阈值，前提是影像直方图呈双峰分布。水体指数影像通常满足这一条件：水体像元聚集在高值端，非水体聚集在低值端。但若水体占比极小（如干旱区的小水库），直方图双峰不明显，Otsu 可能失效。工程上可先做直方图统计，确认双峰性再决定是否使用。

from skimage.filters import threshold_otsu
import numpy as np

def otsu_water_mask(mndwi, valid_mask=None):
    data = mndwi[valid_mask] if valid_mask is not None else mndwi.ravel()
    data = data[np.isfinite(data)]
    thresh = threshold_otsu(data)
    return mndwi > thresh, thresh

6.3 分位数法与双阈值法

分位数法假设水体在指数直方图中占据固定比例，取第 90 或 95 分位数作为阈值。这一假设在时序监测中尤其危险——丰水期和枯水期的水体占比差异巨大，固定分位数会系统性高估或低估。双阈值法（如 MNDWI > 0.1 为确定水体，0 < MNDWI < 0.1 为候选水体）则通过引入“不确定带”来降低误判，但需要额外的辅助数据（如坡度、纹理）来消解候选区。

6.4 阈值策略对比

本文评述：没有“万能阈值”。工程上更务实的做法是——先用固定阈值快速出图，再根据目视检查结果调整；若需批量处理，则用 Otsu 配合人工抽检；若精度要求极高，则引入双阈值加辅助特征。阈值的选择应当被记录在元数据中，作为结果可复现性的一部分。

七、后处理：形态学、连通域与边界优化

7.1 形态学操作

二值水体掩膜常含有椒盐噪声（孤立像元）和孔洞（水中的小岛或误判）。开运算（先腐蚀后膨胀）可去除孤立小斑块，闭运算（先膨胀后腐蚀）可填补小孔洞。结构元素的大小需与目标水体尺度匹配：对于坑塘，3×3 或 5×5 足够；对于大型湖泊，可适当增大以平滑边界。

7.2 连通域过滤

按面积过滤连通域是最有效的去噪手段。设定最小面积阈值（如 0.5 公顷或 10 个像元），剔除小于该阈值的斑块。阈值的设定需结合研究目标：若关注大型湖泊，可设得较大；若关注零星坑塘，则需设得较小，但会保留更多噪声。本文评述：连通域过滤的阈值本质上是一个“尺度先验”，它隐含了研究者对“什么算水体”的定义，应当在方法部分明确说明。

7.3 边界规整

栅格转矢量后，边界常呈锯齿状。可通过两种方式规整：一是对栅格先做平滑（如高斯滤波后再二值化），二是对矢量做简化（如 Douglas-Peucker 算法）。前者改变像元归属，后者只改变几何表达。本文评述：若面积统计是核心目标，建议保留原始栅格边界；若制图美观是核心目标，可对矢量做适度简化，但需记录简化容差。

八、栅格转矢量：参数、陷阱与坐标系

8.1 转换工具与参数

GDAL 的 polygonize 是最常用的栅格转矢量工具，QGIS 中的“栅格转矢量”和 ArcGIS 的 Raster to Polygon 底层均调用它。关键参数包括：是否简化边界（simplify）、是否只保留指定值（mask）、是否连通八邻域（8-connected）。八邻域连通会把对角相邻的像元视为同一区域，通常更符合视觉直觉，但可能产生自相交多边形。

# 使用 GDAL 命令行
gdal_polygonize.py water_mask.tif -f "ESRI Shapefile" water.shp water_id

# 使用 Python (rasterio + shapely)
import rasterio.features
import geopandas as gpd
from shapely.geometry import shape

with rasterio.open('water_mask.tif') as src:
    mask = src.read(1)
    transform = src.transform
    shapes = rasterio.features.shapes(mask, mask=(mask==1), transform=transform)
    geoms = [shape(s) for s, v in shapes]
    gdf = gpd.GeoDataFrame({'geometry': geoms}, crs=src.crs)
    gdf.to_file('water.shp')

8.2 坐标系陷阱

栅格数据通常以投影坐标系存储（如 UTM），转矢量后若未正确设置 CRS，会导致面积计算错误。面积统计应在等面积投影下进行（如 Albers 等积投影），而非地理坐标系（WGS84）或 UTM（UTM 在带边缘有变形）。本文评述：很多“面积对不上”的问题，根源是投影选择不当，而非提取算法有误。

8.3 拓扑检查

转出的矢量可能存在自相交、重叠、缝隙等拓扑错误。建议在导出前用 shapely 或 PostGIS 做一次拓扑修复（make_valid），并检查多边形是否闭合。对于时序分析，还需确保不同时相的水体矢量使用同一坐标系和同一边界范围。

九、精度评估与不确定性分析

9.1 精度指标

水体提取的精度评估通常基于混淆矩阵，核心指标包括总体精度（OA）、用户精度（UA）、生产者精度（PA）和 Kappa 系数。对于水体这一类不平衡问题（水体占比通常远小于非水体），Kappa 系数可能产生误导，建议同时报告 F1 分数和交并比（IoU）。

9.2 验证样本的获取

验证样本可通过高分辨率影像目视解译、无人机影像或实地测量获取。样本应覆盖不同水体类型（清洁水、浑浊水、小水体）和不同背景（城市、农田、山地）。本文评述：验证样本的独立性至关重要——若用训练阈值时看过的区域做验证，精度会被高估。建议采用分层随机抽样，并保证每类样本量不少于 50 个。

9.3 不确定性来源

水体提取的不确定性主要来自四个方面：传感器噪声与大气校正误差、混合像元（水陆过渡带）、阈值选择的主观性、后处理参数的任意性。其中混合像元是最难消除的——一个 10 m 像元若一半是水一半是岸，其光谱是两者的加权平均，指数值介于水陆之间，无论阈值如何设定都无法完美分类。本文评述：在精度报告中，应当明确区分“可消除误差”和“不可消除误差”，后者是数据分辨率的物理极限。

十、典型场景实战：从平原湖泊到山区河流

10.1 平原大型湖泊

以洞庭湖、鄱阳湖为例，水体面积大、边界相对清晰，MNDWI > 0 即可获得较好结果。主要挑战在于：丰枯水期面积差异巨大，固定阈值在枯水期可能把湿生植被误判为水体；湖中洲滩的光谱与陆地接近，容易被漏分。建议结合 NDVI 做辅助判断——NDVI 高且 MNDWI 低的区域判为植被，NDVI 低且 MNDWI 高的判为水体。

10.2 山区河流

山区河流窄（常不足一个像元宽）、比降大、阴影多。NDWI 在阴影区容易误判，MNDWI 因引入 SWIR 而有所改善，但仍有残余。工程上可引入坡度掩膜——坡度大于一定阈值（如 15°）的区域排除水体可能。此外，山区河流的混合像元问题尤为突出，提取出的矢量往往断断续续，需结合流向数据做连通性修复。

10.3 城市水体

城市水体（公园湖泊、河道、景观水池）面积小、形状规则，但背景复杂——深色屋顶、沥青路面、建筑阴影都可能被误判。MNDWI 相比 NDWI 在此场景下优势明显，因为 SWIR 波段对建筑材料的反射特性与水体差异更大。本文评述：城市水体提取是检验指数鲁棒性的“试金石”，若某方法在城市场景下表现良好，通常在其他场景也不会太差。

十一、前沿进展与混合范式展望

11.1 深度学习语义分割

近年来，U-Net、DeepLab、SegFormer 等语义分割网络被广泛用于水体提取，在复杂场景下精度显著优于指数法。但深度学习依赖大量标注样本，且模型泛化能力受训练区域限制。本文评述：深度学习不是指数法的替代，而是补充——指数法可用于生成初始标签，经人工修正后训练网络，形成“指数预标注—人工精修—模型迭代”的闭环。

11.2 多源数据融合

光学与 SAR 融合是水体提取的重要方向。SAR 对水体表面粗糙度敏感，平静水面呈镜面反射，后向散射极低，与光学指数形成互补。在云雨频繁地区，SAR 可填补光学数据的空缺。本文评述：融合的关键不是简单叠加，而是建立光学指数与 SAR 后向散射之间的物理关联，避免“数据堆砌式”融合。

11.3 时序水体与动态监测

随着 Sentinel-2 和 Landsat 免费数据的积累，时序水体监测成为可能。Google Earth Engine 等云平台提供了便捷的时序分析能力。但时序分析对阈值一致性要求极高——若不同时相使用不同阈值，面积变化可能来自阈值漂移而非真实变化。本文评述：时序分析中，建议固定阈值策略，或使用相对变化检测（如与基准期比较）来规避绝对阈值的漂移问题。

十二、结论与操作清单

回到本文的主线：水体提取的精度瓶颈不在指数公式，而在阈值决策与后处理链路。NDWI 和 MNDWI 的公式早已写定，但如何选阈值、如何去噪、如何转矢量、如何评估，才是决定最终成果可用性的关键。以下是一份可落地的操作清单，供工程实践参考。

水体提取是一项“细节决定成败”的工作。指数法看似简单，但每一个环节都有值得推敲的地方。希望本文的梳理能为读者提供一条清晰的操作路径，也为后续引入更复杂的方法打下坚实基础。