小浪底水庫水體時空變化特征與驅動因素分析

閆戈丁，景海濤，王 磊

（1. 河南理工大學測繪與國土信息工程學院，河南 焦作 454000）

水庫是調節水資源時空分布的重要措施，承擔著防洪、發電、航運、供水、環境保護等多方面的任務[1-2]。作為衡量水庫儲水量的一個重要指標，水體面積發生不同程度的擴張或萎縮可能會導致一系列生態環境問題，進而影響經濟社會發展[3]，因此監測水庫面積尤為重要。遙感技術可快速高效地提取水域和監測水體表面變化，具有宏觀、實時和低成本的特點，已被逐漸應用于水體識別與監測領域[4-5]。目前，地表水提取和動態監測的主要方法包括單波段法、水體指數法、譜間關系法、閾值法等[6]，曹榮龍[7]等利用Landsat TM 影像構造了修訂型歸一化水體指數，用于檢測密云水庫水面面積變化；周晗[8]等基于Sentinel-1/2衛星影像，對比分析了單波段法、水體指數法和監督分類法，得出歸一化水體指數（NDWI）法分類精度最高的結論；陸天啟[9]等研究松濤水庫面積的時序變化發現，水體面積呈先降后升的趨勢，并分析了氣候因素對水體面積變化的影響。然而，對水庫進行長時序以及月度細致時間序列的研究卻較少。Google Earth Engine（GEE）是Google提供的對全球尺度地球科學資料（衛星數據）進行在線可視化計算和分析的云平臺，可為長時序研究提供解決方案。用戶訪問并使用云平臺可快速、批量處理數量“巨大”的影像，極大地提高了效率[10-11]。本文通過GEE 云平臺獲取Sentinel-2 影像并提取小浪底水域面積，分析了水體面積變化的時空特征，進而討論了影響水體面積變化的驅動力。

1 研究區概況與數據來源

1.1 研究區概況

小浪底水庫位于黃河中游豫、晉兩省交界處，庫區范圍為三門峽水庫大壩至小浪底大壩之間的水域范圍。水庫集水區處于峽谷地段，地勢西北高東南低。庫區屬溫帶大陸性季風氣候，降水量年際變化較大，主要集中于夏、秋兩季，冬季雨量稀少。

1.2 數據來源

本文采用GEE平臺提供ALOS DSM:Global 30 m數據集和Sentinel 系列衛星。Sentinel-2 具有赤道5 d 重訪周期和中緯度地區2～3 d的高重復頻率。每顆Sentinel-2 衛星都有一個MSI 傳感器，包括13 個光譜波段，波長范圍為440～2 200 nm，覆蓋從可見光到短波紅外，具有10～60 m的3種高空間分辨率[12]。本文采用Sentinel-2 L1C 產品10 m、20 m 處的兩個波段（B3、B11），共選取2018-01-01—2021-12-31的574景遙感影像。為探究氣候因素對小浪底水庫水體面積變化的影響，本文采用中國1 km分辨率逐月近地表平均氣溫數據集與中國1 km分辨率逐月降水量數據集，數據來源于國家地球系統科學數據中心（http://www.geodata.cn）[13]。

2 研究方法

2.1 水體指數法

Mcfeeters S K[14]提出NDWI 用于水體提取。NDWI利用綠光波段最大化水體的反射率與近紅外波段最小化水體反射率的特點，有效提高了水體與非水體像元間的光譜差異；但當研究區域中有城鎮時會導致水體提取不精確，誤將建筑識別為水體。然而，MNDWI能很大程度地解決該問題[15]，降低建筑物對水體提取的影響，還能區分水體與陰影。MNDWI 是目前使用最廣泛的水體提取方法之一。本文采用MNDWI 從Sentinel-2 影像中提取小浪底水庫的水域面積，計算公式為：

式中，Green、MIR 分別為B3、B11波段。

使用水體指數方法的關鍵在于分割閾值的選擇，通常選擇0 作為閾值來分離水體但不一定是最佳分割。單張影像可反復調試選出最佳閾值，但對于大量時間序列影像的閾值分割，該方式費時費力，因此本文選用Otsu 算法計算提取水體的自適應閾值[16]。Otsu算法的基本原理是根據圖像的灰度特性將圖像分成背景和目標兩部分，再通過計算背景和目標之間的類間方差選取閾值，該閾值可將前景和背景最大程度地區分開[17]。

2.2 基于像元的水淹頻率

水淹頻率對于監測洪水和地表水體的變化趨勢具有重要意義，水淹頻率大的地區水陸交替頻繁；反之，水陸交替發生較少。水淹頻率是指每個像元在一定時間內被識別為水體的次數占總觀測次數的比例，可反映水體空間分布的變化特征[18]。本文以像元NDWI 值是否大于所計算的分割閾值為水體、非水體的判別標準。

2.3 基于GEE的庫區水域面積提取

基于GEE提取水體的具體步驟為：①篩選影像并去云處理，在GEE 云平臺上選取Sentinel-2 L1C 產品（Sentinel-2嵌入了具有云掩碼信息的位掩碼頻段，通過該頻段進行去云處理），利用Filter 方法篩選出2018—2021年小浪底庫區的遙感影像，對影像進行去云處理，獲取云量小于20%的像素點；②合成月度數據并計算MNDWI，逐月對遙感影像的每個像元取均值，得到月度合成影像，選擇B3 和B11 波段逐像元計算MNDWI；③計算分割閾值和水淹頻率，通過Otsu 算法計算MNDWI 影像的分割閾值，像元值大于閾值的為水體，否則為非水體，利用閾值逐一計算每個像元，像元值大于閾值的定義為1，反之定義為0；④提取水體并計算水體面積，利用分割閾值識別水體與非水體，但有少量山體陰影被識別為水體，利用ALOS衛星的DEM數據對提取的水體進行過濾，剔除坡度大于25°、高程大于300 m 的區域，得到較精準的水體分類；⑤統計水體像元總個數，計算水體面積。

3 研究結果與分析

3.1 小浪底水庫水體時間變化

2018—2021 年月尺度水域面積見圖1，可以看出，最大水域面積（234.51 km2）出現在2018年1月；最小水域面積（83.07 km2）出現在2021 年7 月；水庫水量和湖泊水量變化十分相似，受季節性影響存在枯水期和豐水期，枯水期出現在7—8月，豐水期出現在12 月—次年1 月；2018 年、2019 年水體面積呈“V”型變化，2020年水體面積變化趨于平穩，這與同年氣候變化平緩相符，2021年6—7月水體面積有一個劇烈下降，這是由于此時河南多地出現罕見的強降雨，導致水庫需要緊急減少庫容。

圖1 小浪底水庫月度水體面積變化

2018—2021年小浪底水庫水體面積具有明顯的季節波動和月際變化規律（圖2），水體面積整年呈先降后升的變化趨勢，最大水域面積出現在每年的1—3月，最小水域面積出現在7—8月，8—12月水體面積再次回升。小浪底水庫水域面積變化與自然湖泊水域變化呈負相關關系，這是由于夏季多雨水庫需提前降低庫容來應對防洪壓力，冬季少雨水庫需蓄水來保證生產生活用水和灌溉。

圖2 2018—2021年小浪底水庫水體面積變化

3.2 小浪底水庫水體空間變化

2018—2021 年小浪底水庫水淹頻率分布見圖3，可以看出，小浪底水庫水淹頻率總體呈內高外低的分布格局，水庫中心渠道水淹頻率最高，邊緣水淹頻率較低，表明水庫中心擁有永久性水體，邊緣則有少數期次為水體。本文統計了不同水淹頻率下的像元個數，計算得到不同水淹頻率的面積及其比例（表1），可以看出，水淹頻率較低的區域主要分布在石寺鎮、北冶村、西溝村、石井村、大路村、下馬新村、古城鎮、王茅鎮以及莘莊新村附近，其中石寺鎮、古城鎮、王茅鎮和莘莊新村尤為顯著，這4 個區域中水淹頻率為20%～80%的部分均有出現，表明該區域水陸更替頻繁，是防范旱澇災害的重點區域；小浪底水庫有34.84%的區域不存在水陸交替，面積約為85.29 km2（F=100%），該區域一直存在水體，說明水庫水體變化幅度較大，只有小部分區域為永久性水體，其余部分均存在水陸交替現象，其中約71.57 km2（占比29.23%）的區域屬于高水淹頻率區，存在居民在河道中進行耕種的現象，應勸導居民避免在該區域進行耕種活動，約66.49 km2（占比27.15%）的區域屬于中等水淹頻率區（20%＜F≤60%），由于水陸更替不是很頻繁這些區域最容易產生松懈心理，因此需要重點監控旱澇災害。

表1 2018—2021年水淹頻率統計表

圖3 2018—2021年小浪底水庫淹沒區域分布圖

3.3 精度驗證

為驗證水體分類的準確性，在研究區內隨機選取500個樣本點進行驗證，樣本點由Google Earth的高分辨率影像和Sentinel-2 影像目視解譯生成。由水體分類圖的樣本數據和目視解譯參考數據組成的混淆矩陣用于證明分類準確性[19]。本文選取4 張影像進行驗證，繪制混淆矩陣，并計算總體精度（OA）和Kappa系數，以定量表示提取精度[20]。結果見表2，可以看出，小浪底區域的平均OA為97.75%，Kappa 系數為0.932 5，驗證了水域動態檢測的可行性和有效性。

表2 水體分類精度

3.4 水體變化驅動力分析

3.4.1 氣候因素

氣溫和降水是可能影響水庫水體面積發生變化的因素，氣溫變化會影響水汽的蒸發循環，降雨則會直接提供水源補充水體。本文利用2018—2020年月平均氣溫與水域面積進行相關性檢驗，結果見圖4，可以看出，6—8月降水量較多，均超過80 mm；月平均氣溫也在6—8月達到最大值（約25℃）。雙變量Pearson檢驗結果顯示（圖5），平均氣溫與水域面積之間呈負相關關系（R=-0.824，P＜0.01）；降水量與水域面積之間也呈負相關關系（R=-0.865，P＜0.01），說明降水量與溫度對水域面積變化的影響顯著，隨著月平均氣溫和降水量的增加，水體面積明顯減少，這與水庫夏季排洪、冬季蓄水的特點吻合。

圖4 降水量與月平均氣溫的變化

圖5 2018—2021年小浪底水庫平均氣溫和降水量變化與水體面積的相關性

3.4.2 人類活動因素

通過查閱資料以及利用Google高分辨率影像調研水庫周邊用地類型發現，沿水庫邊緣絕大多數為耕地，耕地面積的擴張導致農忙時期需要用水量增加，與之對應的是河南省農忙季節冬小麥6 月下旬—7 月上旬；隨著庫區范圍內城鎮的發展，建設用地不斷擴張加劇了對區域水域面積的占用，同時生產生活需要用水量也在不斷加大。通過調查研究還發現，部分工業加工廠設立在庫區較近范圍內，存在水體污染的隱患。耕地和建設用地的不斷擴張導致了一系列問題，如陳村鄉、坡頭鄉、段村鄉、南村鄉等位于小浪底水庫周邊的村點出現了崩塌、滑坡、地面塌陷等現象，為今后庫區治理工作提供了參考。

4 結 語

本文基于GEE云平臺，采用水體指數法對小浪底水庫進行水體提取，再利用Sentinel 衛星的高重訪周期獲得了月尺度水體面積數據，結合氣候因素，分別從時間和水淹頻率方面分析了水體的變化趨勢。

1）利用GEE 云平臺可在線處理經過預處理的遙感數據產品，極大地提高了研究效率，節省了儲存空間，且在處理速度上具有量級的提升，可為后續進行大尺度海量數據處理提供有力支撐。

2）小浪底水庫水域面積月際波動較大，最大水域面積約為210～230 km2，出現在每年冬季；最小水域面積約為90～100 km2，出現在每年夏季，具有明顯的季節性變化特點。水庫中央區域水體穩定、邊緣區域變化幅度較大，僅有34.84%的水體不存在水陸交替變化。研究庫區水淹頻率并繪制水淹頻率圖可為不同區域旱澇防控提供合理參考。

3）水體面積的變化與月平均氣溫和平均降雨量有顯著關系，月平均氣溫和降水量增加，水體面積減少。小浪底水庫區域的降水主要集中在6—9月，形成了7—8月夏汛和9—10月秋汛，水庫為保持警戒水位將加大排洪；汛期結束后，為應對冬季枯水期的生產生活和灌溉用水，水庫進入蓄水期，水體面積明顯增大。