基于多种水体指数的城市碧道水体遥感提取方法对比研究

2022-03-01刘海洋胡泓达刘旭拢刘万侠

广东水利水电 2022年2期

刘海洋,胡泓达,刘平,刘旭拢,刘万侠

(1.广东省水利水电技术中心，广州 510635；2.广东省科学院广州地理研究所,广州 510070；3.南方海洋科学与工程广东省实验室，广州 511458)

万里碧道是以水为主线，兼顾生态、安全、文化、景观、经济等功能的生态廊道，快速准确获取河湖水体分布对碧道规划建设管理及碧道生态功能研究和评价具有重要意义[1-2]。卫星遥感技术在水体提取应用中发挥的作用日益凸显[3]。水体遥感提取研究从最初的人工目视解译，到基于光谱特征的半自动提取与分类，再发展到结合光谱特征与空间信息的方法[4]。其中，基于光谱水体指数的方法由于其易于理解、计算简便等特征，应用范围最为广泛。

遥感水体指数使得快速实现大范围碧道水体提取成为可能[5]。波段的不同组合和计算公式产生了多种水体指数。MCFEETERS[6]于1996年首次提出了归一化差异水体指数(NDWI)，充分抑制了植被等非水体信息。徐涵秋[7]提出了一种改进的归一化差异水体指数(MNDWI)，通过将近红外波段替换为短波红外波段，有助于消除来自土壤和建筑物的干扰。相关学者利用Landsat7 ETM+影像各波段的自然对数来反映反射系数和相互作用条件，创建了水体指数WI2006。WI2015水体指数[8]则以WI2006为基础，利用线性判别分析分类确定最佳分割训练区类别的系数，进一步提高了水体提取精度。

哨兵2号影像具有20 m分辨率的远红外(SWIR)波段和10 m分辨率的可见光(VIS)和近红外(NIR)波段，将SWIR波段的空间分辨率提高到10 m即可计算得到10 m分辨率的水体指数[9]。亟待解决的关键问题是如何准确提高SWIR波段的空间分辨率，即遥感影像空间降尺度问题。一般来说，空间插值和全色锐化是两种最流行的用于遥感影像空间降尺度方法。传统的双线性插值方法没有引入其它高分辨率波段信息，而是直接应用于低分辨率波段，计算简便但效果有限。相比之下，全色锐化方法是以遥感影像的高分辨率全色(PAN)波段的可用性为前提，广泛应用于具有低分辨率多光谱波段和高分辨率PAN波段的遥感影像[10]。

1 研究区域与方法

1.1 研究区域与流程

研究区域位于广州市区，总面积超过100 km2(如图1所示)，研究区域内主要分布珠江航道、海珠湿地等都市型碧道水体，另外还有一部分细小水体,对碧道水体提取造成干扰的主要是城区分布的大片建筑物及其阴影。

图1 研究区域示意

本文的研究流程见图2。首先利用哨兵2号遥感影像绿和近红外波段计算NDWI水体指数；然后分别通过双线性插值与GS波段锐化将远红外波段空间分辨率提高至10 m，与可见光波段计算得到10 m分辨率的MNDWI和WI2015水体指数；针对各水体指数影像进行掩膜阈值训练，利用分水岭算法实现水体对象分割与提取；最后利用研究区域地表水体真值数据，对不同波段降尺度方法与水体指数组合的碧道水体提取精度开展评价对比。

图2 研究流程示意

1.2 水体指数计算

NDWI是基于遥感影像绿波段和近红外波段组合而成，计算公式如下：

(1)

式中：

ρGreen——绿波段的地表反射率；

ρNIR——近红外波段的地表反射率。

NDWI主要利用近红外波段水体吸收能力强，而植被反射性强的特点，达到抑制植被和突出水体信息的目的。

将NDWI中的近红外波段替换为远红外波段，计算公式如下：

(2)

式中：

ρSWIR——远红外波段的地表反射率。

MNDWI能够显著增强水体与建筑物的对比度，降低两者之间的混淆度，有利于准确提取城镇地区的水体信息[11]。

WI2015则是基于各波段自然对数创建，计算公式如下：

WI2015=1.7204+171ρGreen+3ρRed-70ρNIR-45ρSWIR1-71ρSWIR2

(3)

1.3 远红外波段降尺度方法

空间插值和全色锐化算法通常被用来实现低分辨率波段的空间降尺度。Gram-Schmidt (GS)是应用最广泛的全色锐化算法之一，具有最优的光谱保真度，可保持全色锐化前后影像光谱特征的一致性，即降尺度得到的高分辨率波段数据保留了原始低分辨率波段的光谱特征[12]。

哨兵2号影像不具备全色波段，而GS全色锐化需要融入类似全色波段的高分辨率波段信息。VIS/NIR多光谱波段分辨率均为10 m，但波段之间不同程度的相关性会导致数据冗余。主成分分析方法PCA被广泛用于遥感多光谱影像波段信息降维[13]。利用PCA将VIS/NIR多光谱波段压缩成一个类全色波段。通过对VIS/NIR四波段进行线性变换，生成相互正交的光谱空间，其中第一主成分包含的信息最丰富。据此可将第一主成分视为10 m类全色波段，通过GS方法融合20 m远红外波段，将远红外波段空间分辨率提高至10 m。

1.4 标记控制的分水岭水体分割算法

标记控制的分水岭算法尤为适合水体分割，与最大类间方差等使用单阈值分割水体/非水体的算法相比，其在水体边缘表现更好[14]。

使用该算法进行水体提取包括3个步骤：

① 标记水体/非水体区域：针对各水体指数影像，分别标记可信度高的水体和非水体区域；

② 梯度影像生成：将Sobel算子应用到各水体指数影像，计算生成相应的梯度影像，用于确定水体和非水体区域标记之间的分界线；

③ 执行水体分割：基于水体/非水体区域标记和梯度影像，分水岭算法迭代扩展每个标记，直到所有未标记的像元均被标记为水体或非水体。

1.5 水体提取精度评价指标

采用用户精度(User’s Accuracy,UA)，生产者精度(Producer’s Accuracy,PA)，总体精度(Overall Accuracy,OA)和Kappa系数等4项评价指标进行水体提取结果评价。计算方法分别如下：

(4)

(5)

(6)

(7)

式中：

TP——正确提取的水体像元数目；

FP——错误提取的水体像元数目；

FN——未被提取的水体像元数目；

TN——正确提取的非水体像元数目；

T——整幅影像总像元数目,∑=(TP+FP)×(TP+FN)+(FN+TN)×(FP+TN)。

2 实验结果分析

2.1 水体指数计算结果

图3分别展示了NDWI水体指数、BIL和GS方法对应的MNDWI和WI水体指数空间分布状况。不难看出，3种水体指数均能够较好地突出水体与陆地之间的反差，水域边界已经较为明显。相比较而言，依靠双线性插值方法降尺度的远红外波段计算得到的水体指数值域范围比GS方法分布更广。

2.2 水体遥感提取精度评价

利用分水岭算法在各水体指数影像上实现水体/非水体对象分割，从而实现水体区域的精准提取。将水体提取结果与真实水体数据进行叠加对比，即可得到水体误提区域与漏提区域，分别如图4a和图4b所示。除了对水体误提和漏提结果进行定性直观表达，另外利用精度评价指标对不同降尺度方法与水体指数组合的水体提取结果开展定量评价，评价结果见表1。图4a、b的误提和漏提结果分别对应表1的生产者精度UA和用户精度OA。

a 水体误提结果

表1 不同降尺度方法与水体指数组合的水体提取结果精度评价

在城市区域易误分为水体的主要是建筑物阴影。重点关注建筑物密集分布区域，在原始影像上呈现出明显的阴影现象，与水体同属低反射率地物。可以看到各水体指数影像提取水体结果一定程度均会受到建筑物阴影混淆影响。相比较而言，NDWI指数和GS方法对应的WI指数水体提取结果均表现较好，只有少数建筑物阴影被误提为水体。特别是NDWI指数具有最优的UA精度(95.02%)，说明NDWI对应的水体误提情形最少，这与已有研究结果不太相符。一般认为MNDWI、WI等水体指数在处理建筑物及其阴影干扰时表现比NDWI更优，因此出现水体误提概率相对更低。然而，对于哨兵2号影像，由于计算MNDWI和WI指数需要的10 m分辨率远红外波段是通过空间降尺度而得，该过程必然会引入降尺度误差，直接降低水体指数计算精度，进而一定程度上会影响水体提取结果。

BIL方法对应的MNDWI和WI两种指数均易将建筑物阴影误分为水体，水体提取精度基本全面落后于NDWI指数。这是由于BIL插值方法虽然一定程度上保持了低分辨率远红外波段的原有光谱信息，但是由于没有引入其它高分辨率波段信息，对空间细节提升作用有限。结果表明高分辨率波段信息融入对于改进低分辨率波段降尺度结果尤为重要，可显著降低水体误提概率，仅依靠简单的空间插值方法难以达到相同的效果。

整体来说，每一种组合都取得了不错的水体提取结果，UA等各项精度评价指标均在80%以上。而GS方法对着的WI指数则提供了最优的PA精度(89.10%)、OA精度(96.79%)与Kappa系数(0.897)。综合来看，利用GS降尺度波段信息计算得到的WI指数取得了最佳的水体提取效果。这也说明选择可利用高分辨率波段信息的空间降尺度方法对于提升水体提取精度尤为重要。

研究区域内分布着众多细小水体，3种水体指数均出现了漏提现象。这是由于GS全色锐化和BIL空间插值均是将参与水体指数计算的所有波段空间分辨率统一至10 m，自然无法提取宽度低于10 m的亚像元级别细小水体。