2014—2019年北京密云和官厅水库时空变化分析

王梦琦 张 文 孟令奎

1武汉大学遥感信息工程学院，湖北 武汉，430079

随着遥感技术不断发展，遥感影像的数据质量和分辨率不断提升，遥感技术也逐渐成为水库时空变化监测的重要手段。孙建芸等［1］利用GF-1影像提取丹江口水库面积进行分析；马艳敏等［2］基于MODIS数据对2000—2017年月亮泡水库、向海水库和查干湖面积进行时空变化分析；龚文峰等［3］利用Landsat8 OLI数据提取磨盘山水库面积进行分析；郝小翠等［4］利用MODIS数据分析甘肃省主要水库的面积变化特征。但目前对北京密云水库和官厅水库时间序列的面积变化分析较少，本文以北京地区密云水库和官厅水库为研究对象，采用GF-1数据对水库面积进行提取和变化分析，结合气温降雨等因素分析南水北调工程对密云水库和官厅水库的变化影响。

1 研究区概况及数据预处理

1.1 研究区概况

密云水库位于北京密云区，是北京最大的唯一水源供应地，位于潮河、白河中游偏下，主要水流来源为潮白河［5］，负责供水、防洪和灌溉。官厅水库横跨北京和河北，位于河北省张家口市怀来县和北京市延庆县界内，官厅水库是我国建国后所建设的第一座大型水库，主要水流来源为河北怀来永定河，在防洪、灌溉、发电等方面发挥了巨大的作用。密云水库和官厅水库在北京地区的分布如图1所示。

图1 密云、官厅水库分布图Fig.1 The Distribution of Miyun Reservoir and Guanting Reservoir

1.2 数据选取及预处理

本文所采用的数据源是GF-1卫星的宽幅相机WFV所获取的1A级多光谱遥感影像（影像来源：中国资源卫星应用中心http：//www.cresda.com/cn/）。选择2014年1月至2019年12月每月15日左右水库处无云或云量较少的数据，最终选取了130景GF-1卫星的WFV数据。其中密云水库2018年6月影像缺失，官厅水库2014年9月影像缺失。对获取的遥感影像进行波段合成、辐射定标、大气校正和正射校正，本文以2017年1月影像为基准对其余时间序列的影像进行配准。

2 研究方法

2.1 水体提取方法

由于密云水库是丘陵型水库，官厅水库为山谷型水库，一年中3月至10月水库周围被植被覆盖。为了更好的区分水库和植被，基于归一化水体指数（normalized difference water index，NDWI）进行水体提取，NDWI可以有效消除植被信息，突出水体，NDWI计算公式为［6］：

式中，Green表示绿光波段；NIR表示近红外波段。

一般的，对于NDWI影像进行水体和非水体分割的阈值选取为0，但由于本文所选的水库分布位置不同且时间跨度较大，选择NDWI=0作为水体和非水体的分割阈值得到的水库提取效果较差，因此对不同影像选取不同阈值进行水库提取，本文选取的水库分布于山脉间，水库周围植被较多，在NDWI影像上，水体的值较大主要集中在NDWI≥0处，而背景值较小主要集中在NDWI＜0处，对NDWI影像进行像频值统计常常会出现双峰现象，在NDWI值较大处集中的峰可以认为是水体，而在NDWI值较小处集中的峰可以认为是背景，一般两个峰值间的像频值最小的位置可以认为是最佳分割阈值，用来区分水体和非水体，为了实现自动化阈值选取，需要找到影像的两个像频值最大的峰值，两个峰值之间的像频值最小处的NDWI值即是影像最佳分割阈值［7］，实现水体自动化提取。

2.2 结冰水面提取方法

京津冀地区冬季平均气温仅为0～8℃，在1月、2月和12月会出现部分水库结冰现象，结冰水面和水体的反射特性差异较大，使用NDWI无法提取出结冰水面。结冰水面和积雪具有相似的反射特性，但传统的积雪指数需要短波红外波段，GF-1影像不具备短波红外波段，无法使用积雪指数进行结冰水面的提取。在遥感影像中结冰水面在各个波段的反射率均较高，本文采用四波段相加之和自适应阈值法进行影像分割，利用结冰水面在4个波段均有较高反射率的特性区分冰面和背景，实现对水库中结冰处反射率较高、在影像上颜色较浅的冰面进行提取［8，9］。

2.3 水库面积变化指标

对水库面积变化采用定量分析和定性分析，定量分析采用水库面积变化幅度指标和水库面积变化强度指标。水库变化幅度指标可以反映水库在一定时间范围内的面积变化幅度，有正变化和负变化之分，水库面积变化幅度指标的表达式为［10］：

式中，L1表示在研究时间范围内水库面积变化总幅度；Ua表示研究期初水库面积；Ub表示研究期末水库面积；T表示研究时间间隔年数；L2表示在研究时间范围内水库面积年变化幅度。当L1＞0时，水库面积增加，水库面积为正变化，L1＜0时，水库面积减小，水库面积为负变化。

水库面积变化强度指标可以描述水库在不同时期的相对面积变化，使用该指标可以定量描述水库的面积变化强弱［11，12］，其表达式为：

式中，C表示水库面积变化强度；ΔA ba表示研究期a～b时间内水库的面积变化；A a表示a年水库的面积；Δt表示研究时间间隔年数。

3 结果与分析

3.1 密云水库面积变化分析

密云水库在2014年2月、2017年1月和2019年2月影像出现结冰现象，利用NDWI自适应阈值方法无法提取出结冰水面，以2014年2月影像为例，如图2所示，而采用四波段相加之和的自适应阈值法可以很好地提取出结冰水面范围，结合NDWI的提取结果实现对结冰水库的完整提取。

图2 密云水库2014年2月结冰水面提取结果Fig.2 Extraction Results of Miyun Reservoir Frozen Water in Feb.2014

对2014—2019年各年内的水库面积变化进行统计得到的结果如图3所示，2014—2019年密云水库面积总体变化趋势如图4所示。

图3 密云水库2014—2019年各年面积变化Fig.3 Area Change Each Year of Miyun Reservoir in 2014 to 2019

图4 密云水库2014—2019年总体面积变化Fig.4 Area Change of Miyun Reservoir in 2014 to 2019

1）密云水库2014年年内面积变化总体呈现逐月减小的趋势；2015年年内面积变化总体呈现先减小后增大的趋势，变化波动范围较小；2016年、2017年和2018年年内面积变化总体呈现逐月增加的趋势，其中2016年水库面积变化波动范围较大；2019年年内面积变化呈现波动增加的趋势，变化波动范围很小。

2）2014—2019年间密云水库面积总体呈现先减小后增大的趋势，如表1所示，2019年面积较2014年极大增加。在2014—2019年间，2019年6月密云水库面积达到最大值为145.65 km2，2015年7月面积达到最小值为64.39 km2，2015年平均面积最小，2015年之后面积逐渐增加。2014—2015年水库面积呈现负变化，其余年份间均呈现正变化，2015—2016年面积变化强度最大。

表1 2014—2019年密云水库年变化幅度、年变化强度Tab.1 Annual Change Range and Intensity of Miyun Reservoir in 2014 to 2019

3）密云水库在夏季出现暴雨、大暴雨时，水库面积会出现一个增加的峰值，如2014年7月和2015年8月，其他时间密云水库面积变化无明显规律，受人为因素影响较大。

4）密云水库是北京重要的供水水库，承担北京的生活、灌溉用水，自南水北调工程通水后，极大缓解了密云水库的供水压力，由于近年来南水北调工程输水、周围水库向密云水库输水、密云水库流域降雨增加、供水减少的“三增一减”，有效解决了密云水库干涸问题［13］，近年来密云水库面积逐年增加，2019年面积达到近年来最大值。

3.2 官厅水库面积变化分析

对2014—2019年各年年内的官厅水库面积变化进行统计得到的结果如图5所示，2014—2019年官厅水库面积总体变化趋势如图6所示。

图5 官厅水库2014—2019年各年面积变化Fig.5 Area change each year of Guanting Reservoir in 2014 to 2019

图6 官厅水库2014—2019年总体面积变化Fig.6 Area Change of Guanting Resevoir in 2014 to 2019

1）2014—2018年年内官厅水库面积变化呈现先增大后减小再增大的趋势，2019年年内面积变化总体呈现先减小后增大的趋势，除2019年外，官厅水库各年年内面积变化波动范围均不大，如表2所示。官厅水库面积一般1月至4月开始增加，随着温度升高，冰雪消融，河流结冰处融化，水量增加，水库得到补充，加上水库水量蒸发小，灌溉用水需求量不大，水库面积增加；水库面积一般5月至8月逐渐减小，由于温度继续升高，蒸发量增加，且夏季灌溉用水需求量增加，水库面积减小。由于京津冀地区降雨主要集中在夏季，会发生多次中到大雨，甚至暴雨和大暴雨，会在降雨时间结束后出现一个短暂的面积增加，但由于夏季蒸发量大于降雨量，降雨结束后一段时间水库面积会继续缩减；水库面积9月至12月逐渐增加，由于气温逐渐降低，蒸发量减小，官厅水库面积有所增加。

表2 2014—2019年官厅水库年变化幅度、年变化强度Tab.2 Annual Change Range and Intensity of Guanting Reservoir in 2014 to 2019

2）在2014—2019年间官厅水库面积总体呈现波动增加趋势，最大面积出现在2019年1月，为101.26 km2，最小面积出现在2015年9月，面积为55.93 km2，在2014—2019年间，官厅水库面积最大值一般出现在1月至4月或9月至12月，面积最小值一般出现在6月至9月。2014—2018年官厅水库平均面积逐年增加，2019年平均面积较2018年略微减小，2018—2019年官厅水库面积为负变化，其余均为正变化，2015—2016年年变化强度最大。

3）官厅水库面积变化季节性规律较强，受气温和降水影响较大，官厅水库自1997年不再向城市供应生活用水，主要负责灌溉和防洪，因此面积变化受人为因素影响较小。

4）随着近年来官厅水库流域降雨量增加和南水北调工程的供水，缓解了官厅水库对工农业用水的供水压力，使得官厅水库面积增加，2019年面积达到近年来最大值。

4 结束语

本文利用NDWI自适应阈值法实现对水体的有效提取，对于结冰水面，采用4波段相加之和的自适应阈值法进行提取，结合提取结果对水库面积进行统计及分析，得到的结果如下：

2014—2019年间密云水库流域冬季气温较官厅水库流域气温更低，密云水库2014年2月、2017年1月和2019年2月均出现了结冰现象，官厅水库冬季没有出现结冰现象；由于密云水库具有供应城市生活用水、防洪灌溉等作用，面积变化受人为因素影响较大，因此密云水库各年内面积变化季节性规律不强；由于周围水库输水、南水北调输水，密云水库流域降雨量增加和密云水库本身供水减少，密云水库近年来面积逐渐增加。由于1997年之后官厅水库不再供应生活用水，官厅水库面积变化受气候和降水影响较大，受人为因素影响较小，因此官厅水库各年内面积变化季节性规律较强，近年来，官厅流域降雨量增加，南水北调工程供水，缓解了官厅水库的工农业用水压力，官厅水库面积增加。

本文所选取的遥感影像为月度数据，为了实现精确的水库变化监测，后续可以选取时间间隔更小的遥感影像进行水库变化分析。