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個波段均有較高反射率的特性區分冰面和背景，實現對水庫中結冰處反射率較高、在影像上顏色較淺的冰面進行提?。?，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年之后官廳水庫不再供應生活用水，官廳水庫面積變化受氣候和降水影響較大，受人為因素影響較小，因此官廳水庫各年內面積變化季節性規律較強，近年來，官廳流域降雨量增加，南水北調工程供水，緩解了官廳水庫的工農業用水壓力，官廳水庫面積增加。

本文所選取的遙感影像為月度數據，為了實現精確的水庫變化監測，后續可以選取時間間隔更小的遙感影像進行水庫變化分析。