地表水体是全球水循环的重要组分之一^[1]，是包括人类在内所有生物体的重要生存资源^[2]。受人类活动^[3-4]与自然过程^[5-6]的影响，地表水体处于持续变化当中^[7-8]。研究地表水体面积的变化趋势，有利于及时分辨和掌握当前水资源管理中的问题，优化水资源利用和保护^[9]。地表水体研究的关键在于如何从地表识别、提取水体^[10]。遥感影像具有大范围、便捷、长时间序列等特点，基于遥感的地表水体研究已成为主流^[11]。

水体指数法因形式简单及其易于使用的特点，是最常见的地表水体提取方法之一。1996年McFeatters^[12]在归一化植被指数(normalized difference vegetation index，NDVI)的基础上，基于水体与植被在绿波段反射率的差异，提出快速从遥感影像中提取水体的归一化差异水体指数(normalized difference water index, NDWI)。此后，基于不同研究目的，国内外学者构建出诸如：MNDWI(modified NDWI)^[13]、EWI(enhanced water index)^[14]、RNDWI(revised normalized different water index)^[15]、NWI(new water index)^[16]、NEW(new scheme extracted water body)^[17]、AWEI (automated water extraction index)^[18]、WI2015(water index 2015)^[19]、MBWI(multi-band water index)^[20]、WI2019(water index 2019)^[21]等一系列水体指数。

然而，水体指数的适用性存在差异，与其构造时间先后无明确关系^[22-23]，在使用时需要进行对比筛选。已有的水体指数对比研究，多局限于几种水体指数^[24-26]，对比结果差异较大，少有研究对十多种常用水体指数的分类精度进行综合对比。例如: 刘浩等^[27]对4种水体指数提取鄱阳湖水体信息进行对比，发现MNDWI效果最佳；王一帆和徐涵秋等^[25]对比MNDWI、AWEI、WI2015这3种水体指数，同样认为MNDWI整体更优；而邱煌奥等^[22]在鄱阳湖区域对比3种水体指数的结果表明，NWI要优于MNDWI。

水体指数法提取水体需要设置分类阈值。关于分类阈值的选择，一般使用大律法(OTSU)，基于类型间方差值划分阈值。但一些研究表明，大律法划定的阈值并不能最好地区分水体与非水体^[28-30]，人工调试虽然工作量大，但划定的阈值分类结果往往更好^[31]。

另一方面，验证分类结果的精度，是对比并筛选水体指数的重要依据。常见精度验证法可概括为抽样检验^[29]和面积检验^[32]2种。抽样检验是基于高分辨率影像或实测数据，随机选取一定规模水体样本点和非水体样本点^[33-34]，并结合混淆矩阵来验证精度，该方法具有易于操作的优点，但样本点数量随机性较大；面积检验是将分类结果与同时期高分辨影像目视解译结果做对比^[35-36]，反映不同水体指数提取出的水体面积与真实水体面积的差值^[37]，通过叠置分析在空间上展示错分、漏分的水体像元^[38]。土地利用转移矩阵可反映像元在空间上转化情况^[39]，将该方法用于水体指数提取结果精度验证中，能检测水体指数提取结果的漏分、错分误差。

基于上述分析，本研究选择石家庄市作为研究区，选用常见的11种水体指数，测试基于Landsat 8 OLI遥感数据源不同水体指数提取地表水体的差异，并使用同时期10 m分辨率Sentinel-2 MSI影像作为参考，使用结合转移矩阵的面积检验法，对比抽样检验法，筛选出最佳水体指数。根据选出的最佳水体指数，分析南水北调后，石家庄市永久地表水体面积的变化趋势，为区域水资源管理和利用提供参考。

1 研究区与研究方法 1.1 研究区介绍

石家庄市位于37°27′~38°47′N, 113°30′~115°20′E之间(图 1(a))，是河北省的省会城市(图 1(b))，区域面积约为14 069 km²。该区域西部地处太行山中段，东部为滹沱河冲积平原，地势西北高东南低，地貌组成复杂，包括：山地、丘陵和平原等。辖区内河流分属海河流域大清河水系和子牙河水系。属温带季风气候，降水多集中在夏季。多年平均地表水资源量约为7.8亿m³，是水资源极度短缺的地区之一。

南水北调中线是本地区重要的调水工程，用于缓解水资源短缺以及地下水超采问题，石家庄市引江水量分配定额为7.3亿m³，极大地增加了当地可用水资源量。该工程沿京广线西侧布设，总干渠由高邑县北渎村进入石家庄，经赞皇、元氏、鹿泉、新华区，从新乐市出境(图 1(a)红线部分)。工程自2014年开始向周边及北部的京、津地区供水。石家庄及周边区域增加了外调水利用，同时借助置换出水资源进行生态补水，大大增加了当地地表水体的面积。及时掌握南水北调后地表水体的增长态势，可以有效利用调水资源、防止地表水无序开发，保障南水北调工程的长期功能发挥和促进区域社会、经济、农业等可持续发展。

1.2 数据来源

在Google Earth Engine(GEE)云计算平台中，选择Landsat/LC0/C01/T1_SR数据集^[40]，该数据集空间分辨率为30 m，经过辐射校正、地形校正以及大气校正处理。从中选择2018年9月的遥感影像进行水体指数提取精度对比研究，影像Path/Row为124/33、124/34。作为验证样本的Sentinel-2 MSI影像，来自COPERNICUS/S2数据集，空间分辨率为10 m，影像获取时间同样为2018年9月。

用于提取永久地表水体的Landsat 8 OLI影像，时间范围为2013年5月—2020年12月，云量覆盖度≤5%，符合条件的共285景(图 1(c))。基于影像数据中QA波段进行云掩膜处理，基于DEM数据计算坡度去除山体阴影。数字高程DEM选择美国地质调查局雷达地形测绘SRTM数据，空间分辨率为30 m。

1.3 水体提取方法

选用常见的11种水体指数进行对比，具体水体指数计算方法见表 1。

1.4 转移矩阵

转移矩阵将土地利用变化转移面积以矩阵的形式表达，具体刻画出用地类型变化^[41]。反映在水体分类精度上，将作为验证的Sentinel-2 MSI影像目视解译为水体、非水体二分类标准图，与水体指数识别的Landsat 8 OLI影像中水体、非水体二分类图作对比，其中标准水体、非水体像元的转化结果，即水体指数提取水体结果的错分、漏分误差。转移矩阵的数学表达式如下

$ S_{\mathrm{iv}}=\left|\begin{array}{cc} w_{\mathrm{iV}} & N_{\mathrm{vi}} \\ w_{\mathrm{Vi}} & N_{\mathrm{iV}} \end{array}\right|, $

(1)

式中：S_iV为水体、非水体转移矩阵；w_iV为水体正确分类的面积；w_Vi为水体转换到非水体的水体面积；N_Vi为非水体转换到水体的面积；N_iV为非水体正确分类的面积。

1.5 精度评价体系

采用抽样精度检验法, 结合转移矩阵的面积精度检验法，建立基于混淆矩阵的评价体系，对各水体指数方法的水体提取结果进行精度评价。抽样精度检验法，在标准水体图中，结合水体和非水体面积，随机选出水体样本1 000个，非水体样本9 000个，样本总容量为10 000。将样本点与各水体指数提取结果叠置，得到样本点属性转化表，水体和非水体相转化的点，即错分或漏分的水体、非水体像元。

结合转移矩阵的面积检验法，利用各水体指数法提取出的二分类结果与标准水体图，计算转移矩阵，生成类型转化表。表中水体和非水体部分，即正确分类像元，相互转化部分为错分或漏分像元。根据转化表计算总体分类精度(O_A)、用户精度、Kappa系数、F₁-Score来检验水体指数的分类精度。

总体分类精度(O_A)表示在分类结果中，所有正确分类像元占总像元的比例, 公式为

$ O_{\mathrm{A}}=\frac{T_{\mathrm{w}}+T_{\mathrm{N}}}{T_{\mathrm{w}}+F_{\mathrm{w}}+T_{\mathrm{N}}+F_{\mathrm{N}}}, $

(2)

式中: O_A为总体分类精度，T_w为水体像元正确分类的个数，T_N为非水体像元被正确分类的个数，F_w为水体像元被错分为非水体的个数，F_N为非水体像元被分为水体的个数。

用户精度(U_A)表示在分类结果中，正确分类的水体像元占所有被分类为水体的像元比例, 公式为

$ U_{\mathrm{A}}=\frac{T_{\mathrm{w}}}{T_{\mathrm{w}}+F_{\mathrm{N}}} . $

(3)

生产者精度(P_A)表示在分类结果中，正确分类的水体像元占真实水体像元的比例，公式为

$ P_{\mathrm{A}}=\frac{T_{\mathrm{w}}}{T_{\mathrm{w}}+F_{\mathrm{w}}} . $

(4)

误判误差C_E和漏判误差O_E，分别反映分类结果中某一类别被错误分类、没有被分类出来的概率，公式分别为

$ C_{\mathrm{E}}=1-U_{\mathrm{A}}, $

(5)

$ O_{\mathrm{E}}=1-P_{\mathrm{A}} . $

(6)

F₁-Score综合了用户精度、生产者精度的评价结果，可综合判断模型整体分类性能。当生产者精度和用户精度都很大时，F₁-Score值会变大。F₁-Score的取值范围从0到100%，数值越高表示分类效果越好，计算公式如下

$ F_1 \text {-Score }=\frac{2 \times U_{\mathrm{A}} \times P_{\mathrm{A}}}{U_{\mathrm{A}}+P_{\mathrm{A}}} \text {. } $

(7)

Kappa系数是综合了生产者精度和总体分类精度参数用于一致性检验的综合性指标，可说明影像分类的精度，公式为

$ \text { Карра }=\frac{O_{\mathrm{A}}-P_{\mathrm{e}}}{1-P_{\mathrm{e}}}, $

(8)

$ P_{\mathrm{e}}=\frac{\left(T_{\mathrm{w}}+F_{\mathrm{w}}\right) \times\left(T_{\mathrm{w}}+F_{\mathrm{N}}\right)+\left(F_{\mathrm{w}}+T_{\mathrm{N}}\right) \times\left(T_{\mathrm{N}}+F_{\mathrm{N}}\right)}{\left(T_{\mathrm{w}}+F_{\mathrm{w}}+T_{\mathrm{N}}+F_{\mathrm{N}}\right)^2} . $

(9)

2 水体提取结果与讨论 2.1 不同水体指数水体识别结果对比

水体指数法分类精度能力对比是综合的，不仅要对比在提取宽阔水面上的精度，还应选择细小水体、城市水体作为对比对象。研究选用滹沱河上游岗南水库作为大型水面样地，岗南水库设计库容15.71亿m³，为大Ⅰ型水库；选择南水北调中线工程线路穿滹沱河处作为小型水体样地，该区域为滹沱河、南水北调、太平河、石津总干渠交汇处，且在石家庄主城区范围内。

使用水体指数法增强水体与非水体之间的差异是最重要目的，反映在遥感影像上，就是扩大水体像元与非水体像元的对比度。图 2是11种水体指数方法在岗南水库及其周围处运算后产生的图像，原始影像为红、绿、蓝3波段合成的真彩色图像，各水体指数法计算后的结果范围皆以由黑到白渐变色显示。

从图 2可以看出，11种水体指数均能从原始影像中识别出水体，但在水体与非水体对比度上，不同水体指数方法存在明显差异。NEW识别结果有明显的雾化，AWEInsh、MBWI、WI2015方法识别结果中水体与非水体对比较弱，非水体部分为灰色，对比度不如黑色明显。EWI、MNDWI、NDWI、NWI、RNDWI、WI2019方法水体为亮白色，非水体部分接近于黑色，水体与非水体部分对比度更大，差异更鲜明。

图 3是11种水体指数方法在南水北调工程穿滹沱河处(简称为穿滹沱河工程)的影像。该区域中南水北调中线工程线路水面宽约44 m，在影像中基本为1个像元，石津总干渠水面宽度29~30 m，在遥感影像中不足或接近1个像元大小。沿滹沱河两岸分别为石家庄市区与正定县县城，城中大量建筑、建筑阴影以及一些公园水体。

受太阳入射角、卫星侧摆、地物高度等因素影响，遥感影像中存在大量的阴影区。当建筑阴影与道路重合，光谱反射率随波长增大而降低^[42]，建筑阴影的光谱特征与水体光谱特征相似^[43]，易造成分类误差。从图 3可以看出，多数水体指数建筑阴影与水体的差异较小，不容易判别，影响了分类的准确率。特别是AWEInsh、AWEIsh、MBWI、NEW、NWI等指数，水体像元与建筑像元之间的差异不够明显，识别出的水体部分建筑物较多，误差较大。WI2015、WI2019这2种指数分辨率较高。

另一方面，受遥感影像的分辨率影响，细小河流及小块水体如：石—津总干渠、南水北调中线调水线路等的识别也存在较大误差。如图中左上角部分所示，EWI指数水体像元与非水体像元虽然对比明显，但水体判别缺失较多。仅AWEInsh、AWEIsh、MBWI、WI2015、WI2019方法对小块水体如人工鱼塘有比较好的识别效果。以上分析表明，WI2015、WI2019这2种指数无论在与建筑物阴影的区分方面，还是在小型水体的分辨方面，准确率都较好。

2.2 不同水体指数法水体提取结果对比

为起到筛选水体指数提取能力的目的，本研究使用人工调试划分阈值。各水体指数的分类阈值为AWEInsh>0;AWEIsh>0;EWI>-0.2;MBWI>-300;MNDWI>0.12;NDWI>0.05;NEW>0.2;NWI>-0.5;RNDWI < -0.12;WI2015>0;WI2019>0。

图 4为11种水体指数对岗南水库水面的提取结果，不同水体指数在大型水体、宽阔水面的提取能力基本相同，与同时期Sentinel-2 MSI 10 m分辨率遥感影像的目视解译结果相比较，主要差异分布在水库上、下游河道处(图中用圆圈标注)。AWEInsh、MBWI、NDWI和NEW方法在提取河道水面上有缺失，但11种水体指数对宽阔水面的提取能力差别不大，均能将宽阔水面的主体部分提取出来。故而大型水体提取结果不能反映出水体指数间的差异，需要选择细小水体做对比。

图 5为11种水体指数在南水北调中线工程穿滹沱河区域的提取结果。与同时期Sentinel-2 MSI 10 m分辨率遥感影像目视解译结果相比较，小块水体水面提取较差，如滹沱河与石—津调水总干渠之间的人工鱼塘(图 5中序号②与③之间)，所有水体指数均无法提取出来。其次，细小水系提取结果同样存在较大误差，除AWEIsh、NDWI、WI2015、WI2019方法外，其余指数均未能提取出南水北调工程线路。特别是EWI、MBWI分类精度最差，石—津输水总干渠、南水北调中线工程线路水体信息完全缺失。此外，AWEIsh、MNDWI、NDWI、NWI、RNDWI、WI2015方法的提取结果存在错分像元的现象，主要为建筑阴影像元被识别为水体像元。从提取结果来看，WI2019相较于其他水体指数在抑制阴影的影响上表现最佳。然而由于水体与阴影光谱特征相似的原因，WI2019并不能完全消除阴影，水体分类精度有所下降^[44]。已有研究^[23]结合SAR去消除多光谱影像中的阴影，水体分类精度有所提高。在SAR图像中，水体和阴影均为暗区域，提取水体时仍会提取出部分阴影^[45]，但水体与阴影覆盖区域地形特征存在差异。未来，如何结合SAR图像、水体指数以及地形特征来消除阴影，提高地表水体的分类精度，仍需进一步研究。

2.3 抽样精度检验法结果

利用水体样本点提取各水体指数各点位处的信息，建立混淆矩阵，计算总体精度、Kappa系数、F₁-Score综合评价指标，并基于各指标结果对水体指数的分类精度进行排序(表 2)。

由表 2可以看出，11种水体指数方法的总体分类精度均在97%以上，表明水体指数方法可以有效区分大多数水体和非水体像元。3个反映总体分类情况的指标，总体精度、Kappa系数以及F₁-Score的趋势一致，各水体指数法中分类精度从高到低排名前5位的为WI2015>AWEIsh>WI2019>MNDWI>NWI。

抽样精度检验结果中，错分误差均较小于2%，而漏分误差都在18%以上，MNDWI_S、WI2019_S分别为MNDWI、WI2019水体指数在10 m分辨率Sentinel-2 MSI影像上的抽样精度检验结果，用于解释抽样精度检验结果中漏分误差较大的原因。影像空间分辨率提高10 m时，水体指数能提取出部分30 m分辨率影像无法识别的水体信息，故而30 m分辨率影像水体漏分率较高。而这2个水体指数提取Sentinel-2 MSI影像中水体信息，错分误差高于提取Landsat 8 OLI影像，是由于其提取出的水体样点更多。以WI2019_S为例，将942个点识别为水体，其中分类正确915个，而在Landsat 8 OLI影像中WI2019法提取出811个水体点，其中正确分类799个。因此，同一水体指数的错分误差随分辨率变化产生差异。

2.4 结合转移矩阵的面积精度检验法结果

通过构造标准水体与各水体指数法提取水体结果的转移矩阵，结合混淆矩阵精度评价体系，计算总体精度、Kappa系数、F₁-Score综合评价指标，并基于各指标结果对水体指数的分类精度进行排序(表 3)。

结合转移矩阵的面积检验结果中，各水体指数总体分类精度更高，均在99%以上，这是由于研究区水体面积不足研究区面积的1%，水体和非水体的错分以及漏分面积占比较小所致。2种全局评价指标Kappa系数、F₁-Score结果表明，各水体指数分类精度从高到低排名前5位的为WI2019>MNDWI>AWEInsh>NEW>NDWI。

与抽样检验相比，结合转移矩阵的面积验证判定分类结果中错分误差更大，均在6%以上，主要是因为存在混合像元、山体阴影像元等。Sentinel-2 MSI影像分辨率为10 m, 水体边界信息相较于30 m分辨率的Landsat 8 OLI影像更加清晰，Landsat 8 OLI影像中水体、非水体像元的混合，导致错分误差高于抽样检验结果。此外，在10 m分辨率的影像中，南水北调中线工程线路、石—津总干渠以及城市内部的细小景观河、泄洪道可以被识别出来，而在30 m分辨率影像中无法识别出来。故而，漏分误差多来自于细小水系，属于机械误差，难以消除。

表 3的结果中，WI2019分类精度最高，在抽样检验中分类精度最佳的是WI2015，抽样检验结果存在偏差。相比于抽样检验法，面积检验法更能检测不同水体指数的分类性能，能够将错分像元、漏分像元在空间上进行展示，大大提高检测的可靠性。

上述分析表明，WI2019的水体分类精度更高，更适合用于提取石家庄市的地表水体。

2.5 基于WI2019法的2013—2020年石家庄市常年水体面积变化趋势

使用筛选出的最佳水体指数，基于GEE云计算平台，筛选石家庄区域2013年3月以后(Landsat 8发射)到2020年12月云量覆盖小于5%的遥感影像，共计285景。自然水系面积存在受降水变化影响的趋势，为明确研究区内较长期存在的固定水体，而非降水形成的短暂水体，Zou等^[46]学者依据水体频率指数来消除季节性降水对长期水体分类的影响，将一年中水体频率≥75%的水体认定为常年水体，这与Pekel等^[7]利用水体频率确定全球水体的方法基本一致。本研究也将一年中水体频率≥75%作为分辨永久水体的指标，通过WI2019提取每一年所有遥感影像中的水体像元，将结果进行叠置，获取水体出现频率大于75%的像元，确定为常年水体，从而得到2013—2020年石家庄市常年水体分布图及其变化趋势。

从永久水体空间分布(图 6(a))上看，石家庄市常年地表水体主要为滹沱河及滹沱河上游2个大型水库(黄壁庄水库、岗南水库)，水面变化集中在这些区域范围。除受人类调控影响最大的2个水库外(图 6(b))，南水北调中线工程通水以后，石家庄地表常年水体面积显著增加，从2014年42.0 km²增长到2020年62.5 km²，年际增幅约8.1%。其中重要自然水系滹沱河的有水河面显著增加(图 6(c))，从2013年9.6 km²增长到2020年16.9 km²。其余增加部分主要为城市中的景观河道、新建小型湿地，增长量约12.7 km²。

2017年开始，石家庄市作为河北省生态补水试点，通过引江水、黄壁庄水库、岗南水库向滹沱河全年持续补水，通过生态输水线路向太平河和环城水系补水，使得季节性断流、干涸的景观河道全线通水，为市民提供了良好的生态环境。然而据调查^[47]，石家庄城市明渠(如总退水渠、环城水系)、景观河道(如民心河)、输水渠道(如石-津输水总干渠、南水北调干渠)等多为人工浆砌的硬质沟槽，河道入渗较少、蒸发损失较多。除衬砌河道外，中、小型水库大多进行过防渗处理，此类景观水系、水库水面增长较快，虽然提高了城市的美观程度，但难以起到补充地下水的效果，无效蒸发较多。以相关学者^[48-49]监测白洋淀多年平均水面蒸发量约1 000.0~1 100.0 mm估算，相较于2013年，2020年石家庄市仅新增景观水面就造成超过1 270万m³的无效蒸发。固然景观河道能提升城市形象，给居民带来良好的休闲娱乐场所，但由于景观河道、城市水系等大多数经过衬砌，不能发挥补给地下水的功能，建议适当控制此类水面增长，减少外调水浪费。

3 结论

本研究通过Landsat 8 OLI影像提取常年水体，结合高分辨率影像Sentinel-2 MSI对比常见的11种水体指数的分类精度，筛选出最可靠的水体指数，并对南水北调通水前后2013—2020年间石家庄市常年地表水体的变化进行研究。结果表明：

1) 11种水体指数均能识别、提取大型水体及宽阔水面，但小型水体、河道、城区水系提取结果差异较大，WI2019无论在剔除建筑物阴影影响还是分辨小型水体方面都表现更好。

2) 南水北调通水以后，石家庄市地表水体面积显著增加，从2014年42.0 km²增长到2020年62.5 km²，年际增幅约8.1%。其中经过衬砌，不具备补给地下水功能的景观河道、中小型水库水面增加约12.7 km²，作为生态补水河道的滹沱河常年有水河道面积增加约5 km²。

但调查发现大部分新增水体为人工水系景观，存在人工浆砌的沟槽，地下水补给功能较差、无效蒸发较多，因此建议有关部门适当控制水体规模，以有效减少外调水资源的浪费。