基于Sentinel-1/2的水体提取方法对比研究——以斯里兰卡小型水体为例*

周 晗，叶虎平，魏显虎，张宗科，王法溧，齐永菊，马建如

(1 中国科学院遥感与数字地球研究所, 北京 100101； 2 中国科学院大学, 北京 100049； 3 中国科学院中-斯水技术联合研究中心，北京 100085； 4 中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室， 北京 100101)(2018年2月9日收稿； 2018年5月28日收修改稿)

作为维系地球生态系统稳定和健康的重要因素，水体可以由卫星遥感影像快速、准确地提取[1]。由于水体对电磁波吸收散射的特殊性，目前可用于水体提取的遥感手段很多，光学、红外、SAR等都有较成熟的提取方法，主要有单波段法、监督分类、非监督分类、水体指数等方法。

欧空局发射的哨兵系列卫星包括合成孔径雷达SAR卫星(Sentinel-1系列)和光学卫星(Sentinel-2系列)，具有分辨率高、重返周期短、数据易获取等优点[2-4]，广泛应用于水体提取等各类遥感应用研究。Binh等[2]利用Sentinel-1数据对湄公河流域进行基于神经网络的水体提取；孙亚勇等[3]利用Sentinel-1 A SAR数据快速提取缅甸伊洛瓦底江下游区洪水范围；Du等[4]对Sentinel-2的10 m数据和20 m数据波段融合后计算水体指数，得出6种10和20 m空间分辨率的水体指数，参与水体提取。常用的水体指数包括McFeeters[5]提出的归一化水体指数NDWI (normalized difference water index)，该指数对于城市建成区背景的水体具有局限性。Xu[6]在NDWI指数的基础上进行修改，提出改进的归一化差异水体指数MNDWI(modified NDWI)，能够揭示水体更微细的特征。监督分类和非监督分类方法对于细小水体的提取精度都有限，而单波段法和水体指数法的提取精度与当地的土地利用类型有关，如NDWI指数法不易区分城市建筑物与水体[7]。

斯里兰卡位于南亚次大陆南部，是中国“一带一路”倡议的重要节点国家。该国地处热带，降水量丰富，年均降雨量为800～5 000 mm，而降雨明显受热带季风的支配，中东部雨季集中于10月至次年2月，其他时间处于干季[1]，降水不均衡导致当地居民常常陷入用水困难的情境。斯里兰卡中部山区地势高，东部地区地势低洼。除人工修建的少数几个大型水库外，斯里兰卡东部地区分布的绝大多数是小型湖泊、河流等小面积水体，因此这类水体的提取对农村居民生活、农业灌溉、肾病防治等具有重要意义。

本文选取斯里兰卡中东部7月份干季的Sentinel-1与Sentinel-2影像，利用水体指数法、单波段法、决策树方法、监督分类等不同方法对所选研究区范围的小水体进行提取，分析比较不同方法的分类精度及针对不同情况水体的优缺，得出更适宜南亚小面积水体提取的方法，为哨兵系列数据在一带一路沿线国家的小型地表水体提取和应用提供参考。

1 数据与方法

1.1 研究区概况

研究区范围为81°0′E ～ 81°53′E，7°9′N ～ 8°9′N的斯里兰卡中东部地区，从行政区划上，这部分地区属于斯里兰卡的波隆纳鲁沃(Polonnaruwa)、拜蒂克洛(Batticoloa)、安帕赖(Ampara)、巴杜勒(Badulla)、内勒佳腊(Moneragala)省。斯里兰卡全境地处热带，属于热带季风气候，终年高温，年均气温约为28 ℃，降水丰富，只有湿季(wet season)和干季(dry season)的区别，其中中东部湿季为每年10月至次年2月。由于气候原因，研究区小型水体及周边都十分适宜植物生长，西部的山区常有茂密的森林覆盖，而整个研究区范围内都生长有多种热带植物。同时，研究区范围内西部地区为重要的水稻种植地，分布有大量水田。

研究区内涵盖斯里兰卡所有典型水体，包括面积最大的几个水库和湖泊，MaduruOya水库、Ulhitiya水库、Ratkinda水库、塞纳讷耶凯湖(Senanayake Samdra)等主要大面积水体，以及斯里兰卡中北部平原上数量众多的小型坑塘水体。此外，斯里兰卡最大的河流马哈威利河(Mahaweli River)流经研究区范围。研究区具体范围如图1。

图1 研究区哨兵(Sentinel)-2A真彩色影像(2017年7月27日)Fig.1 Sentinel-2A true color image of the study area(July 27, 2017)

1.2 数据介绍

本文使用哨兵(Sentinel)-1/2影像来自欧空局“哥白尼计划”的地球观测卫星系列的哨兵系列卫星，是由同轨的极轨卫星组成的星座。哨兵1号由A、B两颗合成孔径雷达星座，可以在所有天气条件下提供中高分辨率雷达影像，星座重访周期6 d。哨兵2卫星也是由A、B两颗光学卫星组成的陆地观测星座，星座重访周期为5 d，具有空间分辨率高(光学空间分辨率可达10 m(表1)，雷达空间分辨率可达5 m)和幅宽大的特征。所有数据均通过欧空局的数据共享网站(https:∥scihub.copernicus.eu/dhus/)下载。

由于斯里兰卡位于热带，多云雨天气，Sentinel-2过境拍摄的大部分影像云量较大，低云量的斯里兰卡地面影像数量不多。本文选取2017年7月27日研究区云量较少、质量较好的Sentinel-2 A一景少云光学影像以及时间相近的2017年7月23日Sentinel-1 A的SAR影像用于水面提取研究。二者成像时间相隔4 d，查阅当地天气记录，研究区内期间无大型降雨过程，可认为湖泊和水库水体面积不变或在研究可以接受的范围内。

表1 Sentinel-2 A的详细波段设置表Table 1 Detailed band settings of Sentinel-2 A

1.3 数据预处理

1.3.1 Sentinel-1数据预处理

本文使用雷达干涉宽刈幅IW(interferometric wide swath)模式的Level-1级别地距影像(GRD，ground range detected)，数据类型为单视复数数据(SLC)。数据的处理过程使用由Sarmap公司开发的SARscape软件，该软件架构于ENVI平台之上，处理过程主要包括多视处理、滤波、地理编码、辐射校正和图像镶嵌。

该时段研究区的Sentinel-1数据只有VV、VH两种极化方式的雷达影像，其中，VV极化方式的影像对水分信息更加敏感[8]。然而，VV极化方式的影像将水分含量较高的沼泽、水田等土壤覆盖类型的水分信息更明显地表现出来，使得江河、湖泊等水体与非水体在图像上的差异更小，对比不明显，因而选择VH极化方式的雷达数据用于接下来的研究处理。

1.3.2 Sentinel-2数据预处理

本文使用经过几何精校正的反射率产品Sentinel-2 L1C，对10 m分辨率的Sentinel-2影像使用ENVI5.4对数据进行辐射定标。遥感处理软件ENVI5.4对Sentinel-2的辐射校正及大气校正流程进行集成处理。软件中，对Sentinel-2数据进行大气矫正运用的是Flaash模型。该算法是由美国空军实验室(Air Force Research Labs，AFRL)发起，基于MODTRAN模型，可用于各种多光谱数据的大气校正处理，处理结果可以提供包括表面反照率、水汽柱信息、气溶胶和云的光学厚度、大气温度特性等的地表及大气属性信息[9]。根据影像的成像时间、中心点经纬度、成像区海陆类型、成像波段等参数选择适宜的气溶胶模型，得到Sentinel-2的大气反射率。在软件计算过程中，为降低结果存储空间，ENVI软件对反射率结果扩大10 000倍，因此得出的反射率结果是真实地表反射率的10 000倍。

1.4 各种不同水体提取方法

1.4.1 基于Sentinel-1的决策树方法

雷达图像的亮度值反映雷达的回波强度，根据微波在地表的反射特性，雷达的回波强度取决于其系统参数及后向散射系数。对于特定参数的雷达系统而言，其系统参数不改变，则回波强度与后向散射系数直接相关。雷达的后向散射系数由地物的自身特性决定，包括物理特性(介电常数、表面粗糙度、散射特点)及几何特性(雷达图像获得时的几何关系、坡度、形状等)[8]。研究区内的水体大部分属于平静表面的水体，微波信号在此发生镜面反射，回波信号极强，因此可以利用确定阈值建立决策树提取研究区内的水体。

1.4.2 基于Sentinel-2的单波段法

各类地物对不同波段辐射的反射能力各不相同，在近红外-短波红外波段，水体的吸收强度很大，如725 nm的纯水吸收aw(725)为1.489 m-1，而且随波长增大，aw(835)达3.702 2 m-1 [10]，高出其他地物的吸收好几个数量级。近红外-短波红外波段水体强吸收率而低反射，在影像上表现为水体相对其他地物具有更低的亮度值，与其他地物的弱吸收有明显区别[11]。Sentinel-2的B8近红外波段(835 nm)具有与可见光波段相同的10 m空间分辨率，B8单波段的灰度直方图也表现出十分明显的“双峰特征”，因此可以直接对B8波段的影像亮度值进行阈值分割。

1.4.3 水体指数法

现有的水体指数通常是利用水体在近红外波段、中红外波段的高吸收率和绿波段的高反射率，通过不同的数学运算强调水体而削弱非水体，以达到水体提取的目的[12]。常用的水体指数有归一化水体指数(NDWI)、改进的归一化水体指数(MNDWI)、增强型水体指数(EWI)[13]。其表达式分别为：

NDWI=(Green-NIR)/(Green+NIR)；

MNDVI=(Green-MIR)/(Green+MIR);

EWI=(Green-NIR-MIR)/(Green+NIR+MIR);

Sentinel-2卫星的B3波段为绿波段(Green)，B8、B8A波段为近红外波段(NIR)，未设置中红外波段。但是水体在B12波段(中心波长2 202.4 nm，半高宽242 nm)的光谱反射特性与在中红外波段的反射特性相似，因此用B12波段代替中红外波段(MIR)参与波段计算，得出水体指数。

对计算出的水体指数分别进行非监督分类，设置地物类型为8类，得到研究区内较为精细的地物分类结果。合并地物类别为水体、陆地两类，得出基于不同水体指数的研究区水体提取结果。

1.4.4 监督分类法

根据斯里兰卡东部的土地利用情况，设定水体、沼泽、旱地、森林、裸地、农田、居民地等7种类别，在进行监督分类操作时，选择基于统计学习理论的向量机分类方法(support vector machine)为监督分类的分类器。

地物的后向散射强度不同，其在微波波段的回波强度不同。将Sentinel-2数据与Sentinel-1数据进行融合，可以丰富影像的信息量，且能适当排除云对影像的影响。对Sentinel-1、Sentinel-2的10 m分辨率图像进行数据融合，得到重采样后的4波段融合数据，融合后的影像对于一般水体的目视解译效果高于Sentinel-2影像。融合后及真彩色影像上的水体显示结果如图2。

以融合后的影像为基础，参考斯里兰卡实地考察资料，选定适宜的训练样本，进行监督分类，对分类结果进行适当合并，最终得到“水体”与“非水体”两种类别的最终结果。

2 结果与讨论

2.1 水体阈值分割

1)Sentinel-1雷达影像

选取训练样本，分析训练样本的统计直方图，通过求拐点阈值分割的方法进行水陆二元分离，经过实验比对后，设定本景Sentinel-1影像的分割阈值为0.008 9，亮度值低于和高于所选取的阈值则分别为水体和非水体。

2)红外波段灰度影像

图2 融合前后水体的显示图像Fig.2 Images of water body before and after fusion

图3 影像光谱曲线Fig.3 Spectral profiles of various images

由于水体与非水体这个波段的灰度分布差异明显，分析水体训练样本及整副影像的灰度统计直方图，经过反复实验比对后，选取阈值为1 135，提取出研究区的水体。图3(a)、3(b)分别为影像的灰度分布图和水体训练样本在影像中的灰度分布，纵轴为像元数，由于前述ENVI软件大气校正处理原因，横轴为大气反射率(TOA)的10 000倍。

图4(a)、4(b)为Maduru Oya Reservior区域Sentinel-1影像及Sentinel-2的B8单波段影像。

图4 雷达影像与B8单波段影像Fig.4 Radar image and B8 single band image

由于水体附近土壤水分含量高，复介电系数ε高，在微波波段的辐射能力较强，Sentinel-1影像中水体附近区域非水体与水体的对比不如Sentinel-2的B8单波段影像明显，前者的目视效果显著低于后者。

2.2 全影像各水体提取方法提取结果

图5 不同方法提取结果图Fig.5 Extracting results using different methods

图5为研究区的不同方法水体提取结果，淡黄色部分为陆地，深蓝色部分为水体。分别为基于(a)NDWI指数、(b)MNDWI指数、(c)EWI指数、(d)Sentinel-1决策树、(e)NIR单波段的方法及(f)监督分类法的提取结果。比较不同方法的提取结果，可以得出不同方法在研究区的水体提取差异。以融合图像中选定总面积大致相同的水体和非水体感兴趣区域(ROI)为验证样本，计算不同方法的混淆矩阵，分析不同方法提取结果的精度，得到的精度指标如表2所示。

表2 各方法精度指标列表Table 2 A list of accuracy indexes of various methods %

在以上各种指标中，总体分类精度为正确分类的像元比例，错分误差(commission error)可以表现被错分为水体的其他土地类型，漏分误差(omission error)表现被错分为其他类型的水体占比。

分析结果显示，NDWI指数法和MNDWI指数法具有相对较高的精度，均达到94%以上的总体分类精度，MNDWI法略高于NDWI法约0.4%，达到94.53%，NIR单波段法与EWI指数分类精度也较高，约93%。具体到不同方法分类结果的差异，错分误差最高的为监督分类法，误差达13.88%，对分类精度影响较大，这是由于监督分类法容易将种植有水稻等农作物的水田误分为湖泊；由于只采用阈值进行二元分类，分类标准过于单一，NIR单波段的错分误差也较大；错分误差最小的方法是水体指数法，各类水体指数法的错分误差在5.0%～7.5%之间，被错分为水体的其他土地类型也较少，其中MNDWI指数的分类结果水田划分为水体，基于EWI指数的提取结果常把海边盐碱滩涂分类为水体，NDWI指数对于城市中的水体有较多局限，而研究区的城市化程度不高，在各种水体指数中错分误差最小，为7.22%。由于所用数据的空间分辨率较其他方法的空间分辨率更低，分类结果受山体阴影等因素的影响，基于Sentinel-1的决策树法精度较低。水体指数法的漏分误差位于4.0%～7.0%，由于参与指数计算的所有波段均为10 m空间分辨率，漏分误差最小的是NDWI水体指数法。

综合上述分析，MNDWI指数法和NDWI指数法的总体分类效果最佳。

2.3 典型水体提取效果分析

表3中“细小水体”列为采用不同方法，对马哈威利河部分河段及其支流的提取结果；“云覆盖区”列为光学数据有云区域及云阴影区域的水体提取结果对比；“水陆边界”列举不同方法对于水体与陆地的过渡地区的提取结果。

对于细小的河流、湖泊等规模较小不易从背景中分辨出来的水体，由于NDWI法和NIR单波段方法所采用的Sentinel-2波段10 m空间分辨率，较其他方法空间分辨率最高，且水体在B8近红外波段的强吸收特性明显区别于其他类型地物，而能够更加精细地将其从陆地背景中分离。在水体与陆地具有明显边界的区域(如堤坝)，这两种方法对水陆边界的准确度最高。相较而言，NIR单波段法的提取结果最为精细，不仅可以识别流经的马哈威利河主要河道，对于目视无法分辨的该流域的细小支流或被树木等部分遮挡的河流也能精确提取。NDWI指标同样对水体信息十分敏感，因此利用NDWI指标的方法也能识别出细小水体。

虽然Sentinel-2卫星B8波段数据对水体信息足够敏感，但由于NIR单波段信息提取只利用近红外波段的信息，分类标准过于单一，提取结果易受其他因素影响，容易将湖泊附近的湿润陆地分类为水体(表3)。这也是NIR方法具有比水体指数法更高的误分误差的原因之一。

南亚地区的四季常有丰富的云，光学影像中云的出现无法完全避免。基于不同水体指数的分类方法均无法避免云自身及云阴影带来的影响，水体的提取结果中混入了云阴影区。几种水体指数在不同程度上受到云对分类结果的影响，其中NDWI指数受到的影响最小，主要来自云覆盖区的干扰，而其他几种水体指数则同时受到云覆盖区和阴影区的影响。云阴影对地表信息有部分遮挡，降低了阴影区在各个波段的地表反射，使得阴影区在红外波段的反射强度与水体的反射强度十分相似，NIR单波段法会将云阴影误分为水体。雷达微波波段具有穿云透雨的能力，Sentinel-1决策树法和融合雷达影像的监督分类方法则不受云的影响。

虽然雷达数据不受天气和云的影响，单一的雷达影像用于水体提取时，由于影像噪声、空间分辨率、倾斜成像等原因，分类结果精度并不高，融合图像分类的结果受到雷达数据的影响，对地形起伏较为敏感，最终分类结果中受山体阴影的影响造成部分误差的存在。

表3 典型地区不同方法分类结果对比列表Table 3 Comparison of classification results in typical regions among different methods

3 结论

利用哨兵(Sentinel-)1/2数据，分析常用的水体提取算法，都能一定程度地提取各种湖泊、水库大型水体和线状细长河流水体以及坑洼池塘等小型水体。比较不同方法提取效果，在利用监督分类法的过程中，应当对水田的分类做进一步界定，或在后期根据图像纹理等特征对提取的水体再进行筛选或排除。NDWI法提取结果精度最高，被漏分的水体和被误分为水体的其他土地类型都较少，最适合广泛用于斯里兰卡的水体提取；NIR单波段法对于线状细长河流水体分离结果最好。在雷雨季节，云覆盖量少的光学影像较难获取时，Sentinel-1数据不失为一种补充，基于Sentinel-1的决策树方法适用于阴雨天气时水体的快速提取，可用于南亚地区的洪涝灾害灾情分析。基于小面积水体精确提取后，下一步开展水质和水量时空变化分析，以期对斯里兰卡民生有重大影响的饮用水安全和未知肾病(CKDu)[14]防治提供科学支撑。