地表水变化对生态、社会和经济的影响是当前研究的热点^[1]。光学遥感水资源研究主要可以分为：气候和水资源^[2]、不同类型水体的提取^[3-5]、洪水监测^[6]、河流形态^[7]、水体长时间序列^[8]、海岸带变化^[9]。目前，Landsat影像是用于变化分析，识别和检测地表水的最广泛使用的数据源^[10]。在基于Landsat遥感影像水体提取的研究中，多数是针对大区域的河流、湖泊等水体的提取^[11-14]。近年来，一些学者采用Sentinel 2影像对大区域的地表水进行提取研究，取得了良好的效果^[15-16]。目视解译是最初的光学遥感影像提取水体的方法^[17]，其缺点是效率低、成本高，没有利用不同地物具有不同反射率的特性。遥感影像的水体信息提取是基于电磁波谱上的反映特性^[18]，文献[19]通过使用Landsat 4影像的6个波段提出缨帽变换湿度指数(tasseled cap wetness, TCW)。受到归一化差异植被指数^[20]的启发，文献[21]基于Landsat影像的绿波段和近红外波段提出了归一化差异水体指数(normalized difference water index, NDWI)，NDWI是简单的水体指数, 该指数能较好地抑制植被信息，但受薄云和山体阴影影响较大。文献[22]使用短波近红外波段代替近红外波段，改良了NDWI，提出了改进的归一化差异水体指数(modified normalized difference water index, MNDWI)。MNDWI指数能突出水体信息，但也容易把水体和阴影混淆。基于分类二叉树的思想，文献[23]提出了针对有阴影的场景自动水体提取指数AWEIsh和针对没有阴影的场景自动水体提取指数AWEInsh，二者都是多波段水体指数。AWEI指数在提出时，声明该指数不适用于高反射率的地物，比如冰雪、白色建筑物等。文献[24]在全球大范围多个地区进行测试了自身的稳定性和可行性, 提出了水体指数(water index 2015, WI2015)。针对复杂环境，文献[25]成功地构建了一个多波段水体指(multi-band water index, MBWI)，MBWI一共使用了5个波段并且在多个复杂环境的地区测试了MBWI的有效性。MBWI能较好地抑制非水体信息的能力，但其缺点之一是不适用提取细小水体。面向对象提取水体信息的方法^[26], 在抑制阴影信息的同时，能够准确提取细小水体。但是面向对象的方法的缺点之一是：不具有普遍性，难以在大范围下实现高精度的水体信息提取。

在地形复杂的环境下，如高山积雪较多的区域，对水体提取产生干扰的因素有：阴影、冰雪等^[27-28]。如何减少遥感影像的噪声干扰仍是水体信息提取的难点。本文分析多种地物的反射率，综合考虑多波段间的关系，通过建立六元不等式组，设计一个顾及多个地物噪声影响的光学遥感影像归一化多波段水体指数(normalized difference multi-band water index, NDMBWI)。

1 多光谱影像水体指数模型构建 1.1 已有的水体指数

选取7种典型的、常用的水体指数进行研究，7种常用的水体指数见表 1。这7种水体指数以绿波段和近红外波段、短波近红外波段为主要影响成分。TCW、MNDWI、AWEI、WI2015主要考虑绿波段和短波近红外波段的反射率，因为水体的绿波段反射率大于短波近红外波段的反射率，而其他非水体地物的绿波段小于短波近红外波段的反射率。正常情况下，这些指数都可以区分水体与非水体。但是，TCW、MNDWI、AWEIsh、WI2015会把可见光波段反射率大于短波近红外波段的非水体地物分类为水。文献[25]指出，雪就是可见光波段的反射率远远大于短波近红外波段的地物。

1.2 模型的机理

本文的技术路线如图 1所示。在不同的时段、不同的地区，从多幅Landsat 8、Sentinel 2、高分一号(GF-1)影像里选择不同地物不同波段的真实地表反射率光谱值进行拟合，具体采样点的信息见表 2。内陆水体、海洋水体、阴影、雪、植被、旱地、建筑物、云的平均反射率如图 2所示。图 2曲线是基于表 2中的采样点求出来的平均值。水体的反射主要在蓝绿光波段，其他波段的吸收率很强，尤其是在近红外波段、短波近红外波段有强烈的吸收带。研究已有的水体指数可以从各个光谱间的关系作为切入点，本文以不同波段与绿波段的比值，阐述已有水体指数模型的机理及其存在的缺陷，并建立本文水体指数NDMBWI。表 3为各个光谱与绿波段的关系，表 3中的数值可看成图 2曲线量化后的结果。图 1中，可以把TCW、NDWI、MNDWI、AWEIsh、AWEInsh、WI2015近似看成双波段水体指数。水体指数的模型机理在于不同波段的对不同地物具有不一样的敏感性和有效性。雪、云的反射特性与水体相似，表 3中，部分地物(内陆水体、海洋水体、雪、云)的Q₄、Q₅、Q₆的值都小于1，从机理上说明：主要依赖绿波段、近红外波段、短波近红外波段的双波段水体指数在有雪环境下不能有效地提取水体。

表 3中的数值是来源于Landsat 8、Sentinel 2、GF-1在各个地物、各个波段上的平均值，然后各自相除得到的比重。本文试验数据源为Landsat 8、Sentinel 2、GF-1，具体信息见表 4。

表 4中，Landsat 8、Sentinel 2、GF-1三者的中心波长不完全一致，这也导致了已有的水体指数模型在用于不同的数据源时，表现出阈值跳动的问题。经过研究，已有的水体指数的模型可用式(1)表示

式中，Y表示水体指数的值；x₁、x₂、x₃、x₄、x₅、x₆是6个波段前的系数；ρ₁至ρ₆是蓝波段、绿波段、红波段、近红外波段、短波近红外1波段、短波近红外2波段的反射率。

式(1)整体除以绿波段，可变形为

式中，Y表示水体指数的值；x₁、x₂、x₃、x₄、x₅、x₆是6个波段前的系数；ρ₁至ρ₆是蓝波段、绿波段、红波段、近红外波段、短波近红外1波段、短波近红外2波段的反射率。

根据水体的水体指数的值理论上要大于α，非水体的水体指数的值理论上要小于α。式(2)中，ρ₁、ρ₂、ρ₃、ρ₄、ρ₅、ρ₆的值用比值系数代替, 得到以下不等式组

式中，a、b、c、d、e、f、g、h分别代表求得的比值系数。

在本文把α设置为0，目的是希望本文指数能够以稳定的阈值提取水体信息，到达自动化提取水体的目的。结合式(3)与表 3的参数，组成一个六元不等式组

式中，x₁、x₂、x₃、x₄、x₅、x₆分别是蓝波段、绿波段、红波段、近红外波段、短波近红外1波段、短波近红外2波段的系数。

式(4)的系数源于Landsat 8、Sentinel 2、GF-1在各个地物、各个波段上的平均值，若是直接赋予式(4)每项等于定值，会因为推导不严谨，导致水体指数模型出现错漏现象，本文通过最小二乘法的方式进行消元、拟合，求出不同定值下的解。

由表 3中可以发现，已有水体指数Q₁、Q₂、Q₃、Q₄、Q₅、Q₆与1的大小关系并不能完全地实现水体的识别。例如，NDWI由绿波段与近红外波段组成，即是从Q₄的值可以解析NDWI对水体信息识别的能力，大于1的数值的地物包括阴影、植被、旱地、建筑物、云，小于1的值包括内陆水体、海洋水体、雪，此时导致NDWI无法以固定的阈值区分雪、水体。MNDWI、AWEIsh、AWEInsh、WI2015近似看成由绿波段与短波近红外波段组成的模型，可从Q₅的值可以解析此类型模型对水体信息识别的能力，大于1的数值的地物包括阴影、植被、旱地、建筑物，小于1的值包括内陆水体、海洋水体、雪、云。此时，以上指数只考虑Q₄或Q₅(双波段)，并没有综合考虑Q₁、Q₂、Q₃、Q₄、Q₅、Q₆的内在联系，所以本文建立六元不等式组，综合考虑多波段间的关系，设计一个顾及多个地物噪声影响的水体指数模型。

式(4)没有精确解，其原因是：式(4)是超定方程组的类型，即方程个数大于未知量个数。式(4)的解集是一个区域，解多元不等式(组)就是求它的各个不等式的解集的交集。为了使本文水体指数适用于国产遥感影像(没有短波近红外1波段、短波近红外2波段)，需要提前把x₅、x₆设置为0。本文通过逐步消元、拟合，求出其中的一种解：x₁、x₂、x₃、x₄、x₅、x₆分别是-1、3、2、-5、0、0，即初步构建的水体指数为

式中，ρ₁、ρ₂、ρ₃、ρ₄代表蓝波段、绿波段、红波段、近红外波段的反射率。

式(5)的系数是可行域中的简化解，目的是为了更好推广NDMBWI以及方便记忆。另外，可以通过使用式(5)的参数，反解式(4)进行验证。因为双波段水体指数难以抑制雪、云的影响，本文指数没有将绿波段作为唯一的正参数，额外引入了红波段作为正参数。为了使本文指数归一化，最终拟定的水体指数模型为

式中，ρ₁、ρ₂、ρ₃、ρ₄代表蓝波段、绿波段、红波段、近红外波段的反射率。

本文推算出来的系数，实则为式(6)的分子。在提取水体时，起决定性作用的是分子而不是分母。这样做的好处是：预处理时，在大气校正环节，某些算法为了减少存储空间会把结果扩大10 000倍。此时，本文使用了归一化的形式，则不会受以上的影响，无论大气校正有无扩大10 000倍，NDMBWI的范围始终在-1至1的范围。使用本文水体指数和表 1中列举的已有的水体指数进行仿真模拟，在不同的地区、不同的数据源选择水体、阴影、雪、植被、旱地、建筑物、云等7种地物作为对象，每种地物选择100个样本点，不同水体指数的稳定性测试结果如图 3所示。图 3表明，0可以作为NDMBWI的阈值，能够有效地区分水体与多种非水体；以上的水体指数均可区分水体与植被、旱地、建筑物；除了NDMBWI、MBWI，其他的水体指数均不能区分水体与冰雪；MBWI会对部分水体漏分；NDWI的阈值提高到0.2左右可以分区水体与冰雪。

2 水体提取实例分析 2.1 不同卫星影像在有雪环境下的水体提取

试验1选择的地区位于西藏林芝地区。数据源是Landsat 8影像和Sentinel 2影像，影像拍摄时间均为2020年2月13日。图 3中得知：8种指数以阈值为0均能将水体与部分非水体(植被、旱地、建筑物)区分。所以在本次试验中，NDMBWI、MBWI、WI2015、AWEInsh、AWEIsh、MNDWI、TCW统一以0作为阈值进行水体提取，而NDWI则是使用0.2作为阈值。提取雅鲁藏布江水体的结果如图 4所示。图 4(a)为Landsat 8影像的原始图。图 4(b)为Sentinel 2影像的原始图。对比图 4(c)、(d)可以发现，NDMBWI应用于Landsat 8影像、Sentinel 2影像时，均可以抑制了雪的影响，完成了提取水体信息。

在试验1中，选择100个水体、非水体样本组成混淆矩阵，计算出的Kappa系数、总体精度、错分误差、漏分误差、用户精度、生产者精度见表 5，使用Sentinel 2影像的效果均优于Landsat 8。为了进一步研究水指数的阈值，研究NDMBWI、MBWI、WI2015应用于Sentinel 2影像的数值分布，三者的直方图如图 5所示。在直方图中，最优阈值一般选定在“双峰”之间的低谷区间。图 5(a)验证了在试验1中，0可以作为NDMBWI的最优阈值。图 5(c)并无明显的低谷，说明WI2015本身难以抑制冰雪的干扰。

2.2 有云环境下的高分一号影像水体提取

云一直是光学遥感影像难以处理的问题之一，本文提出的NDMBWI适用于存在少量云的情景。云一般都有阴影，传统的水体指数无法处理阴影。由仿真模拟的结果可知：NDMBWI能够以0为阈值区分水体与云、阴影。GF-1是中国自主研制的第一颗民用高分辨率光学卫星，为我国水环境的遥感监测提供更高空间分辨率的光学遥感数据源。但是，GF-1只有可见光波段与近红外波段，使用了短波近红外波段的MBWI、WI2015、AWEI、MNDWI、TCW等指数不能应用在GF-1。试验2以中国香港的迪士尼公园为试验区域，影像的拍摄日期为2018年10月2日。使用NDMBWI提取水体的对比结果如图 6所示。

图 7作为NDMBWI在试验2的直方图，验证了0可以作为NDMBWI的最优阈值。云是光学影像无法处理的痛点，去云步骤一般分为两步：①识别云像素；②把云像素设置为背景值(无效值)。本文在设计指数模型时把云考虑到非水体地物的因子，所以通常情况下NDMBWI在会把云识别为非水体地物，云通过NDMBWI得到的值为负值。即NDMBWI无须做去云步骤。但是，值得说明的是NDMBWI适用于云量为20%以下的光学遥感影像。也就是说NDMBWI是无法适用于云量过多的场景。

2.3 多地区水体提取

选取表 2中的齐齐哈尔、巴彦淖尔、梧州、三亚作为试验3，从多个不同的地区，复杂的地形环境下提取不同类型的水体，验证本文水体指数模型的稳定性，图 8中的NDMBWI均以0作为阈值。

3 结论和讨论

如何减少噪声干扰对遥感影像的不利影响是水体提取的难点，本文设计一种归一化多波段水体指数NDMBWI，避免积雪、阴影等的影响。试验1使用相同时相、同一地区的Landsat 8影像、Sentinel 2影像，以西藏林芝地区为试验区域，试验结果说明NDMBWI能抑制雪的影响。试验2使用GF-1影像，在影像含有少量云的情况下完成水体提取，验证NDMBWI可以抑制云、阴影的影响。阈值，对于水体指数是至关重要的。从水体指数的构建开始，要注意阈值的选取问题。试验3选择多个不同地区进行试验，以说明本文水体指数模型的稳健性。今后，构建新的光学水体指数时，应该考虑多个波段间的敏感性和有效性。本文的创新点在于：充分利用可见光波段、近红外波段的信息，不需额外借助辅助数据，即可排除雪、云、阴影的影响，完成水体信息的提取。

从理论上，任何的水体指数都无法保证提取所有的水体，本文提出水体指数虽然在某些情况下能得到较高的精度，并在冰雪、阴影、建筑、云等区分具有一定的优势，然而因为自然水体的复杂性，并不是任意一种水体指数能区分所有水体。也就是说，NDMBWI的最优阈值不能一直保持在0，例如在钱塘江悬浮泥沙较多的出海口地区提取水体时，NDMBWI的阈值应该适当降低。目前主流的水体指数的研究大多数是固定阈值。本文的下一步工作是：基于原有的水体指数模型，构建动态阈值的模型。例如，以NDVI作为NDMBWI的阈值，使用NDMBWI与NDVI模型进行双通道并行处理，在NDMBWI的矩阵中，大于NDVI的像素点可以归为水体，反之为非水体。本文方法稳健，可以作为基于深度学习的水体检测分类的基础，解决难以区分雪、水的问题。