结合DEM与淹没频率的水库水体动态遥感提取优化方法

刘诗燕， 蔡晓斌*

(1.中国科学院大学， 北京 100049；2.中国科学院精密测量科学与技术创新研究院环境与灾害监测评估湖北省重点实验室， 武汉 430071)

地表水资源是制约区域经济发展、保障区域社会稳定、维护区域生态环境的关键[1].水库是地表水资源的重要组成部分，具有供水、防洪、抗旱、发电、灌溉、水产养殖、娱乐等多种生态系统服务价值[2].Chao等通过统计全球水库建设的历史资料重构了人工水库蓄水过程，计算发现由于水库蓄水使得海平面上升趋势减少了30 mm[3]，此外全球范围内大量的水库建设还会对水资源的时空分布格局产生直接的影响[4].因此，动态监测水库水资源量的时空变化对于维护生态环境稳定和保障人类用水安全具有十分重要的意义.

水库水体的水面范围提取是水库相关研究的基础.水体范围的变化能直观反映水量的增减[5]，同时对水循环领域相关的其他应用研究也至关重要.例如，水库蒸发的季节性/年际变化主要受水库表面积的变化驱动[6]，水环境承载力与生态系统的稳定性也与水面积的变化密切相关[7].随着对地观测技术的迅猛发展，遥感技术以其大面积同步观测、时效性强等优点，目前已经在水体范围提取方面有了广泛的应用[8].遥感技术根据传感器的不同通常分为光学、热红外和微波遥感三个大类，其中光学遥感因其历史数据资源丰富、提取精度高的特点，在长时序地表水识别中应用最为广泛[9].

目前常用的光学遥感影像提取水体信息的方法包括监督与非监督分类[10]、水体指数阈值分割法[11]、以及决策树法[12]等，其中阈值分割法应用范围最为广泛.张超等利用人工目视确定阈值对柴达木盆地湖泊水面变化进行了遥感提取[13]，吴川等、孙爱民等利用OTSU最大类间方差法确定阈值，分别对丹江口水库和博斯腾湖水域变化进行了监测研究[14-15].

然而，现有研究中所采用的水体信息提取方法往往运行效率较低，大范围、长时序的水体监测实现难度较大，此外，常规的光学遥感监测阈值方法在长时序、高精度监测上还存在一些明显不足.首先，云及云阴影的影响难以消除[16].一些遥感影像上被云遮挡的地物信息无法获取，在长时序水体提取中通常被舍弃，而未被遮挡区域所包含的有用信息没有得到充分利用，因此水体监测结果也会存在大量的数据缺失.其次，大量水库位于山区，山体阴影造成的水体错分问题难以避免[17]，山体阴影会导致提取结果存在与水库水体的混分现象，这些明显的错分增加了水体动态监测的误差.此外，在水体的长时序监测中，海量的遥感影像下载及处理工作耗时费力，传统的提取方法常常难以满足监测快速响应的需求.

近年来遥感云计算技术的发展为海量遥感数据处理和分析提供了前所未有的机遇[18].谷歌地球引擎(google earth engine，GEE)是一个专门用于卫星遥感数据运算处理的云平台，通过网页浏览器编程运算，大大加快了算法的测试分析与大尺度的快速分析和应用的进程[19].针对全球尺度的水库信息遥感提取，已有很多学者基于该平台进行了相关研究工作.在全球地表水体制图(joint research centre monthly water history， JRC MWH)数据集[20]的基础上，Zhao等提出了针对水库污染像元信息自动纠正的算法，生成了全球6 000多个水库的面积时间序列[21]，该算法能充分提高遥感影像在水体提取中的利用率而被广泛使用[22].然而，全球水体数据集往往更新周期长，且提取结果中存在着一些明显的错误[23]，以JRC全球水体数据集为例，直至2021年初才更新前五年的数据，因此要想实现水库水体范围的动态监测，仍需探索一种简单高效且准确的遥感监测方法.

鉴于此，本研究基于GEE平台，结合Landsat影像和NASA DEM数据，提出了一种结合DEM和淹没频率的水库水体遥感提取DF优化方法.针对全球九个不同的水库实现了长时序水体动态提取，并结合精度验证及测高水位数据集，与Zhao等基于JRC MWH数据集的水库水体动态监测结果进行了对比，分析了新方法的精度与有效性.

1研究区及数据来源

1.1研究区域

全球水库和大坝数据集(GRanD)[24]包含了全球范围内7 250座水库的范围矢量数据及建成时间、大坝经纬度等基本信息.研究从中选取了九个不同地形和水文地质条件的水库作为研究案例，用于测试本论文提出的研究方法.九个水库分布在五大洲不同纬度区，高程分布从230 m到1 740 m.水库具体分布如图1所示，基本信息如表1所示.

图1 九个研究水库分布图Fig.1 Distribution maps of nine reservoir

表1 水库基本信息

1.2遥感数据

1.2.1Landsat影像 所采用的Landsat影像为USGS Landsat 5/7/8 Level 2，Collection 2，Tier 1影像数据集.该数据集包含了经过大气校正的正射地表反射率(Landsat5/7包含三个可见光波段，一个近红外波段以及两个短波红外波段；Landsat8包含四个可见光波段，一个近红外波段以及两个短波红外波段)以及通过CFMASK算法[25]生成的云、雪卷标信息.数据空间分辨率均为30 m，重访周期为16 d，获取时间为1984年1月1日至2020年12月31日.

1.2.2DEM数据 选用的高程数据为NASA-DEM全球数据集，其空间分辨率也为30 m.该数据通过合并ASTER GDEM、ICESat GLAS和PRISM数据集的辅助对NASA SRTM DEM数据进行了再处理，显著提高了其精度[26].其中最重要的改进处理在于使用了改进的相位解缠技术，通过ICESat GLAS数据进行控制来减少原有SRTM DEM数据中的数据空值区.

1.2.3JRC Monthly Water History地表水月历史数据 Pekel等利用1984年3月16日至2020年12月31日期间采集的Landsat 5、7和8的400多万幅影像，通过专家系统对每个像素点分别进行了水/非水分类，并将结果整理成月度历史，用于动态监测全球水体变化[20].该数据集包含1984年3月至2020年12月间442个月的全球地表水的位置和时间分布的地图，并提供了有关水体的范围及其变化的统计数据.

2研究方法

利用Landsat Level2，Collection2，Tier1产品中的云、雪卷标去除影像中的受云、雪影响的无效地表像元，通过OTSU大津法对单一时相影像进行水库水体信息的初步提取，在时序提取结果的基础上统计每个像素的淹没频率.考虑到光学影像中无法回避的无效地表像元问题，根据单景影像提取结果中云、云影、雪等无效像元的长时序淹没频率以及同期正常像元的水体提取结果，修正异常像元的水体覆盖信息.在此基础上，结合修正后的水体淹没频率和DEM数据对每一景影像的提取结果进行优化.通过DEM和淹没频率约束剔除山体阴影及其他异常值的影响，从而最终获得水库水体的时序提取结果.工作流程如图2所示.

图2 处理流程图Fig.2 Processing flow chart

2.1基于MNDWI的大津法水体初步提取

为了突出水体特征，便于大津法提取.首先通过波段运算，计算改进的归一化差分水体指数(MNDWI)[27]，该指数计算公式为：

其中，XMNDWI为XNDWI计算值，Rgreen为绿光波段反射率(TM/ETM+的第2波段，OLI的第3波段)，RSWIR为短波红外波段反射率(TM/ETM+的第5波段，OLI的第6波段).

改进的归一化差分水体指数是在归一化差分水体指数(NDWI)的基础上利用短波红外代替近红外波段，该指数能更好的抑制建筑物和土壤的背景信息从而突出水体.

随后基于MNDWI数据运用大津法自动化确定分割阈值从而完成水体和非水体的区分.大津法即最大类间方差法，是由Otsu于1979年提出的一种无参数无监督的阈值分割算法[28].其实现算法如下：

假设使用阈值T将灰度图像分割为前景和背景；Size表示图像总像素个数，u为图像的平均灰度，W0为前景像素点占整幅图像大小的比例，U0为前景像素点的平均值，W1为背景像素点占整幅图像大小的比例，U1为背景像素点的平均值， 为类间方差，计算公式为：

u=W0×U0+W1×U1,

(1)

g=W0×(u-U0)2+W1×(u-U1)2.

(2)

将(1)式带入(2)式得到类间方差计算式：

g=W0×W1(U0-U1)2.

(3)

该算法根据图像的灰度特性，从图像的最小灰度值遍历到最大灰度值，找到使得分割的两类类间方差最大，即g最大时的灰度值T即为最佳分割阈值[29].

2.2高质量数据筛选及淹没频率计算

由于光学遥感影像存在云雪覆盖对地表地物的遮挡，以及Landsat7部分影像因传感器故障造成了影像中大量的无效像元，研究通过云、雪卷标剔除被影响的无效像元，按有效地表覆盖像元统计分析后将影像数据集分为三类：一类是缺失像元占比小于5%的遥感影像，将其视为高质量影像；第二类是缺失像元占比大于等于5%小于80%的遥感影像，这类影像存在有效信息可以利用，但需要对污染缺失像元进行信息的回填；第三类是缺失像元占比大于等于80%的影像，该类影像因有效信息参考太少，大津法分类误差过大而被舍弃.在三类影像数据分类的基础上，选择高质量影像通过大津法的水库水体提取结果计算水库淹没频率分布，并用于无效像元的信息回填以及提取结果的优化.

2.3污染像元信息回填

Zhao等提出了基于淹没频率的水体无效像元回填算法[30]，该算法的原理是：通过单景影像中有效水体像元淹没频率的频次直方图，计算不同淹没频率在该景影像中有效水体像元中的出现频次，当第一个频次大于当期所有有效水体像元淹没频率下的频次均值的0.17倍时，以该值所对应的淹没频率作为阈值.污染像元对应淹没频率高于该阈值的像元划分为水体，反之则划分为非水体.本研究基于此算法和计算的淹没频率，对缺失像元信息进行回填，以获得更完整、连续的水体信息.

2.4水体DF优化提取结果

水库水体提取修正基于以下假设：水库水体应为连续水体，远离主要水体的破碎斑块应视为误提.鉴于水库水体的特殊性，提取的水体高程应小于坝顶高程，因此DEM超出坝顶高程的水像元应视为误提；修正操作步骤为：

1) 以20%作为淹没频率阈值，得到淹没频率二值图；

2) 以坝顶高程作为DEM阈值，得到DEM二值图；

3) 提取及回填后同时满足淹没频率大于20%和DEM小于坝顶高程的水像元被分类为水，否则被分为非水；筛选后水体栅格转为矢量；

4) 选择水体矢量的最大面要素作为水库最终提取结果.

3结果与分析

3.1方法精度验证

为验证DF优化提取方法的精度，从九个水库的多时相提取结果中各自随机选择了一景影像的提取结果，通过随机生成的900个验证样点，验证了分类结果的准确性，并与Zhao等基于JRC MWH数据集的直接修正结果[30]进行了对比分析.对比结果如表2所示，采用DF方法优化提取的水库水面分类结果总体精度为95.67%，Kappa系数为0.912 6，而基于JRC MWH数据集的直接修正结果的总体精度为82.67%，Kappa系数为0.658 6.此外，图3还分别展示了两个数据的分类结果，可以看出基于JRC MWH数据集的直接修正结果中水库水体存在明显错分现象.

表2 两种方法水库水体分类精度对比表Tab.2 Comparison of classification accuracy of reservoir water body by two methods

图3 分类结果对比图Fig.3 Comparison of classification results

在此基础上，研究还计算了九个水库建成以来至2020年12月的水体面积时间序列结果，并按月合成得到各月份的水面积均值.基于美国农业部官方网站提供的水库水位卫星监测结果[31]，建立水面积与其同期水位监测结果的回归关系，如图4所示.从图中可以看出，优化后水面积和水位之间有较好的一致性，决定系数(R2)均在0.6以上.

通过水面积与水位的相关关系，计算均方根误差，结果如表3所示.实验的九个水库通过DF优化方法提取的水体面积，RMSE值均优于基于JRC MWH数据集的直接修正结果.总体RMSE均值从16.984 3减小到了11.796，其中五个水库减小幅度超过了30%.由此也可证明本论文提出的水库水体范围优化提取方法具有更高的准确性.

图4 水库水面积与水位的回归关系Fig.4 Regression relationship between water area and observed elevation of reservoirs

表3 各水库提取结果精度评价

3.2水库水体面积时间序列变化

利用上述算法，可以从遥感影像中准确地计算出水库的水体面积，从而构建水库水体面积的时序产品.然而，当水库范围内80%以上像元为无效像元时，DF优化方法无法使用.因此，使用该方法在部分月份无法获取到水库水体范围信息.利用其他年份多年该月面积均值替代的方法对缺失的月份面积估计值进行补全.研究所提取的水库面积时间序列结果如图5所示.由图5可以看出，从总体变化幅度上看，采用DF优化方法提取的水面积在1984年到2020年间波动幅度相较于直接修正的JRC MWH数据面积波动幅度更小，水库水面积极高或极低的异常值出现的频率也更低.

图5 九个水库月面积动态时间序列Fig.5 Dynamic time series of nine reservoirs per month

4讨论

水库水体范围的遥感提取受到山体阴影干扰[32]、云的遮挡、雪覆盖以及狭窄河道造成的水体不连续[33]等多种因素的制约，本文结合DEM和淹没频率，建立了DF优化方法消除部分山体阴影造成的水体错分误差，通过建立水体缓冲区消除断裂水体，选择最大水面要素作为水库范围，从而剔除远离水库的细小水体.DF方法在传统大津法水体提取的结果上进行了优化，通过对比JRC MWH数据集的直接修正结果，验证了DF优化提取方法的准确性.

当水陆像元占比接近时，大津法分割的结果最为理想[34]，但原始的影像边界来自于GRanD数据库shapefile矢量缓冲区，水陆占比与水库的形状及大小相关，难以控制确切的比例.因此，提取的最终结果也与大津法的分割效果有一定的关系.基于MNDWI大津法分割的水体提取对遥感影像的质量也有要求，当无效像元过多时，提取结果的可靠性降低，因此存在部分时段的影像被弃用的情况.这种影像弃用情况对于基于JRC MWH数据集的直接修正结果而言并不存在，因此其水库水面积监测数据的时相理论上应更为完整.当然，也必须认识到针对无效值过多的影像基于JRC MWH数据集的直接修正结果虽然不会缺失，但其结果的可靠性也存在较大问题.

表4显示了各水库时序覆盖数据以及月份缺失情况，其中污染像元占比超过80%的影像因无法准确进行大津法提取而被舍弃，水库建成后至2020年间，有效月份数据占比均值为61.98%.基于JRC MWH数据集的直接修正结果中，无效值占比超过90%的月份被舍弃，有效覆盖月份占比均值为69.48%，略高于本研究方法的提取结果.但是，由于JRC MWH数据集提取过程中算法的局部适宜性问题，也会出现部分水库有效覆盖数量明显少于DF优化方法(如Tahtakopru水库).

5结论

本论文提出了一种结合DEM和淹没频率的水库水面积遥感DF优化提取方法，并以全球九个不同类型的水库为研究区，对该方法进行了验证，计算其1984年至2020年间的水面积动态变化过程,得出以下主要结论.

表4 水库覆盖数据统计Tab.4 Reservoir coverage data statistics

1) 相较于基于JRC MWH数据集的水库地表水修正结果，利用DF优化提取方法得到的水体范围具有更高的提取精度，且水库水面积与水库水位之间有更强的一致性，均方根误差更小.

2) 受限于大津法分类原理，DF优化提取方法无法适用于污染像元占比过大的影像水体提取，但总体有效监测月数与已有基于JRC MWH数据集的水库地表水修正结果中的有效数据占比相差不大.而且，缺失时相的影像数据质量本身无效值较多，其结果的有效性也值得商榷.同时，该提取方法不用依托于JRC MWH数据集的更新，动态监测能力更强.

3) 与传统遥感监测方法相比，利用GEE云计算平台能更高效的进行数据处理和统计分析，为后续大流域甚至全球尺度的水库水体动态监测提供了有力支撑.