面向Sentinel-2影像的亚像元级水体提取方法

熊龙海,何颖清,刘茉默,李 俊

(1. 水利部珠江河口治理与保护重点实验室,广东 广州 510611; 2.珠江水利委员会珠江水利科学研究院,广东 广州 510611)

准确获取地表水水体范围及面积对水资源保护和合理开发利用、海岸线变化监测、水环境监测、水资源调查、湿地调查、洪涝灾害评估等方面具有重要意义[1]。卫星遥感技术具有快速、大范围、周期性成像等特点,已经成为一种有效的地表水水体范围提取技术。不同空间分辨率、时间分辨率和光谱分辨率的卫星传感器已经被广泛应用于地表水水体范围的提取[2-8]。在地表水遥感的研究中,长序列、高频次、大范围对水面率、蓄水量、生态流量等水资源要素进行监测通常需要结合多源中高分辨率卫星数据实现。但是在结合高分辨率遥感影像与Sentinel-2的研究中,往往因为空间分辨率差异大导致水体提取结果存在空间尺度效应。因此,为解决空间尺度效应,研究面向Sentinel-2的亚像元级水体提取方法十分必要。

虽然目前有较多遥感水体提取方法,但绝大多数都是像元精度的,主要包括:①通过目视判读进行数字化;②单波段阈值法[9];③多波段光谱水指数法[10-12];④监督或非监督分类算法[13];⑤图像处理方法,如面向对象分析(object-oriented analysis, OOA)[14];⑥目标检测算法,如基于约束能量最小化 (constrained energy minimization, CEM) 的目标检测法[15];⑦上述方法的组合法[16-17]。上述方法将每个像元分为水体像元或陆地像元,即使是水陆混合像元也被分类为水体或陆地像元,无法提取水陆混合像元中水体丰度。然而,水陆混合像元普遍存在于卫星遥感影像上,现有像元级水体提取方法均不能准确获取混合像元中水体丰度,导致水体提取精度不足。

一些学者提出通过软分类、机器学习和光谱混合分析的方法解决水体提取过程中水陆混合像元的问题。软分类的方法并不计算混合像元中水体丰度,其水体提取结果还是像元级。机器学习方法需要大量的训练集,且依赖训练集的质量,限制了该方法的适用性。光谱混合分析方法是一种有效的、不需要大量训练数据的计算混合像元中水体丰度的方法,文献[18—19]基于Landsat 8 OLI数据,利用光谱混合分析方法提出了亚像元级水体提取技术,但该技术尚未在Sentinel-2数据上应用,Sentinel-2数据比Landsat 8 OLI数据波段更多,直接将上述技术应用在Sentinel-2数据上无法充分利用Sentinel-2数据的优势,难以取得最优效果。目前尚缺乏有效的面向Sentinel-2数据的亚像元级水体提取方法。

因此,本文在利用植被红边水体指数(vegetation red edge based water index,RWI)[20]提取纯水体像元的基础上,结合动态选择端元的多端元光谱混合分析方法构建面向Sentinel-2影像的亚像元级水体提取方法(subpixel water extraction method for Sentinel-2,SWES),以提升水体提取精度,获取更加精确的水域面积,将基于Sentinel-2数据的水体提取进行降尺度,解决多源中高分辨卫星数据提取水体时的空间尺度效应问题。

1 研究区及数据

1.1 研究区概况

广州市水系发达,包含多种水体类型,有不同的水陆混合像元,是亚像元级水体提取方法较为理想的研究区。本文在广州市分别选择了3个试验区,每个试验区都包含不同的水体类型、地表覆盖类型和不同数目的混合像元。这3个试验区包括河流、水库和坑塘养殖区,试验区内水域的水深、浑浊度及污染程度各异。选择这3个试验区是因为它们分别代表了不同光谱特征的水体。图1为试验区的地理位置示意图及假彩色合成图。其中,试验区1为水口水库,水体类型以清洁水体为主;试验区2为珠江干流,拥有城市背景下的复杂水生环境,主要污染来自于生活污水和工业废水;试验区3为坑塘养殖区,主要为水深较浅的桑基鱼塘。每个试验区都存在着不同数量的水陆混合像元,其中,珠江干流混合像元最少,水口水库有较多库湾包含较多混合像元,坑塘养殖区存在大量面积不一的坑塘,水陆混合像元数量最多。

图1 研究区位置

1.2 数据及处理

水体提取采用的Sentinel-2数据为从Google Earth Engine平台上获取的地表反射率产品。产品级别为L2A,该产品已经由欧空局预先进行了辐射定标、大气校正等预处理。覆盖研究区的Sentinel-2影像成像日期为2019年12月1日,并将所有空间分辨率不为10 m的波段重采样至10 m分辨率。

选取与Sentinel-2影像同一天成像的Google EarthTM高空间分辨率(<1 m)影像作为真实水体范围参考数据。统一两种影像的坐标系,通过选择道路交叉口、低层建筑角点等控制点对Sentinel-2影像与Google EarthTM高空间分辨率影像进行配准。将Sentinel-2数据叠置在真实的水体范围数据上,从而得到真实的纯水像元及真实的水体丰度图,并用于验证本文提出的SWES。

2 方法构建

卫星影像上的像元可分为3类:水体像元、陆地像元、水陆混合像元。对于地表水提取只需要提取水体像元及水陆混合像元,亚像元级地表水遥感提取结果由水体像元及水陆混合像元中水体丰度组成。本文构建的亚像元级水体遥感提取方法(SWES)共包括3个主要的步骤:①提取纯水像元;②提取水陆混合像元;③结合空间信息的多端元光谱混合分析提取水体丰度。

2.1 提取纯水像元

采用RWI提取纯水像元,RWI为面向Sentinel-2影像的细小水体提取指数,在一定程度上能够消除水陆边界混合像元的影响,在提取纯水像元上效果较好,该指数同样适用于大范围水体提取。RWI使用大气校正后的遥感反射率作为输入,其公式为

(1)

式中,Rrs(560)、Rrs(705)、Rrs(842)、Rrs(865)、Rrs(2190)分别为Sentinel-2第3、5、8、8A 和 12 波段的遥感反射率。

2.2 提取水陆混合像元

水陆混合像元的光谱非常复杂,常常表现出与其他地物相似的特征,因此单纯地利用光谱特征提取水陆混合像元是不可行的。而假设影像中所有像元都是水陆混合像元也是不切实际的,因此本文提出了利用空间信息提取水陆混合像元的方法。水陆混合像元往往在水体和陆地的过渡带上,即在水域的边界上,因此可利用这一空间特征提取水陆混合像元。具体步骤为在纯水像元的基础上搜寻8邻域所有的陆地像元,如图2所示,即在位于(i,j) 处纯水像元的8邻域内,所有的陆地像元被认为是水陆混合像元。

图2 低分辨率影像上水陆边界

2.3 结合空间信息的多端元光谱混合分析

2.3.1 光谱混合分析模型

光谱混合分析能够用来获取水陆混合像元中水体的丰度。常用的线性光谱混合模型假设混合像元的光谱是由构成该像元的地物端元的反射率以其所占像元面积比例为权重系数的线性组合,即

(2)

式中,Riλ为第i个混合像元的反射率;N为端元的数量;fji为第i个混合像元中第j个端元的丰度;Rjλ为第j个端元的反射率;eiλ为残差。在全约束条件的混合像元分解模型中,端元的丰度通常是非负而且满足和为1的约束,即

(3)

对于第i个混合像元,地物的丰度可通过最小二乘法得到,也即是求取使各个波段上残差的均方根误差(root mean squared error,RMSE)最小的解。可表示为

(4)

式中,T为波段数;eik为第i像元上第k个波段的残差。

2.3.2 结合空间信息的多端元动态选择

由于地表水体的光谱随水质变化而变化,浅水情况下还受水深及底质影响。一般情况下,邻近的水体像元的水质、水深、底质条件最为接近,因此,本文使用结合空间信息的水体端元动态提取方案。对于一个水陆混合像元,所有3×3的窗口8邻域内的纯水像元都被认为是这个混合像元的候选水体端元。如图2(d)所示,位于(k,l)处的水陆混合像元,其8邻域内所有纯水像元均被认为是该混合像元的候选水体端元。

3类陆地端元(植被、土壤、不透水面)通过利用纯净像元指数(pixel purity index,PPI)建立端元光谱库,并使用Count-based Endmember Selection(CoB),Endmember Average RMSE(EAR)和Minimum Average Spectral Angle(MASA)优化建立的端元光谱库。

2.3.3 端元组合

解混过程中考虑了多种陆地端元与水体端元的组合,主要分为2端元、3端元、4端元和5端元组合,具体的组合方式见表1。选择阴影端元的目的是消除光照的变化。根据前人的研究[21],分解结果同时满足下列条件的混合像元才被分解:①所有端元组分的丰度在-0.05到1.05之间;②阴影的丰度低于0.8;③各个波段上残差的均方根误差RMSE小于0.025。如果分解的结果无法同时满足上述条件则该混合像元不分解,并将水体丰度设置为0。对于一个水陆混合像元,可能有多个端元组合分解的结果满足上述条件,选择其中RMSE最小的组合作为最佳的端元组合,其分解结果作为最优结果。

表1 端元组合类型

3 结果与分析

3.1 目视分析

由于缺少专门面向Sentinel-2的亚像元级水体提取方法,本文选择可移植到Sentinel-2的自动亚像元水体制图方法(automatic subpixel water mapping method,ASWM)作为对比。在ASWM中,首先利用NDWI提取纯水像元和纯陆地像元,然后根据阈值提取水陆混合像元,最后利用固定端元组合的混合像元分解算法求解水陆混合像元中水体的丰度,纯水像元和水陆混合像元中水体丰度共同构成了最终的水体提取结果。在这一过程中纯水像元的丰度被设置为1,纯陆地像元的丰度被设置为0。

图3—图5分别为利用SWES和ASWM在3个研究区进行亚像元级水体提取的结果。从结果中可以看出:①整体而言,SWES在3个研究区的水体提取结果较好,明显优于ASWM。②对于水陆混合像元,SWES利用空间关系能够准确提取出来,并且获得相应的水体丰度,在结果图上可以明显地看出由水体到陆地的过渡,通过与真实水体丰度对比,可以看出SWES获取的水陆混合像元水体丰度也能与真实值保持较好的一致性。而ASWM根据光谱信息提取水陆混合像元带来较大误差,将较多陆地像元误提为水陆混合像元,而且难以准确提取3个试验区纯水体边界处的水陆混合像元。③对于纯水体像元,SWES能够较准确提取出来,漏提和误提较少,水体提取效果优于ASWM,尤其在水陆混合较严重的区域,如试验区3,SWES比ASWM取得更好的效果,在这些区域由于存在较多水陆混杂,水质、水深条件不一,水体光谱变化强烈,ASWM采用NDWI提取纯水体像元,很难取得较好效果。④通过图4可以看出,建筑及其阴影对ASWM的影响比SWES的大,SWES几乎没有将建筑及其阴影误提为水体,而ASWM将很多建筑阴影识别为水陆混合像元,也将图5中的道路识别为水陆混合像元。⑤两种方法在提取宽度小于3个像元的水体方面都存在一些误差。

图4 SWES和ASWM在试验区2估算水体丰度的结果

图5 SWES和ASWM在试验区3估算水体丰度的结果

3.2 精度评价

RMSE是亚像元分类中广泛使用的误差指标,量化了地物丰度分解的相对误差。本文选择RMSE对建立的亚像元级水体提取模型SWES进行精度评价。计算公式为

(5)

式中,fi(ref)和fi(esti)分别表示真实的和计算的水体丰度图中第i个像元的水体丰度;N为计算的水体丰度图上所有像元的数目。

表2为SWES和ASWM在3个试验区提取水体时的RMSE。整体上看,SWES在3个试验区的平均RMSE为0.147,而ASWM的平均RMSE为0.175,SWES在所有试验区都比ASWM的RMSE更低,特别是在有很多水陆混合像元的坑塘养殖区,相比于ASWM,本文提出的SWES的RMSE降低最多。

表2 两种方法水体提取结果RMSE

表3为SWES和ASWM在试验区提取的水域面积误差。通过表3可以看出,SWES在3个试验区提取的水体面积精度均较高,平均相对误差为8.03%,远低于ASWM的平均相对误差。两种方法提取的水体面积均小于真实面积,这是由于试验区有些区域水体宽度小于3个像元,两种方法皆难以提取。

表3 两种方法提取水体面积误差表

4 结 语

本文首先利用RWI提取纯水体像元;然后根据水陆混合与纯水像元的空间关系,提出用膨胀算法提取水陆边界混合像元;最后在对水陆边界混合像元进行分解时,为了解决地物的类内光谱变化问题,采取了多端元光谱混合分析算法,在提取水体端元时,基于地理相似定律提出了动态选择水体端元的方法,进而构建了面向Sentinel-2数据的亚像元级水体提取方法SWES。

SWES提取水体的精度达到亚像元级别,充分考虑了水体光谱的变化,结合空间信息提取水体端元,减少了解混计算成本。在3个不同水体类型的研究区开展了水体提取试验,结果均表明SWES取得较好效果,平均RMSE为0.147,在3个试验区上的水体提取效果均优于ASWM,尤其在水陆混合像元较多的坑塘养殖区。SWES在试验区获取的水体面积也有较高精度。由此可见,SWES既能有效提取纯水体像元又能有效提取水陆混合像元及水体丰度,准确获取水体面积,将基于Sentinel-2数据的水体提取进行降尺度,解决多源中高分辨卫星数据提取水体时的空间尺度效应问题。