小浪底水庫(kù)水體時(shí)空變化特征與驅(qū)動(dòng)因素分析

閆戈丁，景海濤，王 磊

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

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

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

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)來(lái)源

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

2 研究方法

2.1 水體指數(shù)法

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

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

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

2.2 基于像元的水淹頻率

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

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

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

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

3.1 小浪底水庫(kù)水體時(shí)間變化

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

圖1 小浪底水庫(kù)月度水體面積變化

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

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

3.2 小浪底水庫(kù)水體空間變化

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

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

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

3.3 精度驗(yàn)證

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

表2 水體分類精度

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

3.4.1 氣候因素

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

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

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

3.4.2 人類活動(dòng)因素

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

4 結(jié) 語(yǔ)

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

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

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

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