閆戈丁，景海濤，王 磊

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

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

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來源

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

2 研究方法

2.1 水體指數(shù)法

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

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

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

2.2 基于像元的水淹頻率

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

2.3 基于GEE的庫區(qū)水域面積提取

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

3 研究結(jié)果與分析

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

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

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

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

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

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

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

表1 2018—2021年水淹頻率統(tǒng)計表

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

3.3 精度驗證

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

表2 水體分類精度

3.4 水體變化驅(qū)動力分析

3.4.1 氣候因素

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

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

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

3.4.2 人類活動因素

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

4 結(jié) 語

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

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

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

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