2006-2020年南四湖水體面積變化分析

周 石，李英冰

（武漢大學(xué)測繪學(xué)院，武漢430072）

南四湖位于山東省西南部，是華北地區(qū)最大的淡水湖泊，是南水北調(diào)東線工程主要調(diào)蓄樞紐之一，在維護經(jīng)濟發(fā)展和區(qū)域生態(tài)平衡方面，具有重要的戰(zhàn)略意義［1］。但是自20世紀(jì)80年代以來，由于各類污染型中小企業(yè)無節(jié)制發(fā)展，沿湖群眾對湖區(qū)濕地盲目圍墾及對湖區(qū)資源的掠奪性開發(fā)，致使南四湖湖區(qū)的生態(tài)環(huán)境明顯惡化，水域面積也不斷縮?。?］。目前國內(nèi)學(xué)者對南四湖的景觀格局［3］、生態(tài)［4，5］、氣候［6-8］、水質(zhì)［9-11］及水位變化［12］做了大量研究，為南四湖環(huán)境的改善做出了巨大貢獻(xiàn)。遙感監(jiān)測手段因其感測范圍大、更新時間快、實時動態(tài)且對被調(diào)查對象不產(chǎn)生破壞等一系列優(yōu)勢，被廣泛用于湖泊等水體的信息提取和水資源的監(jiān)測與保護，駱成鳳［13］利用長時間序列中分辨率遙感影像數(shù)據(jù)，通過人工提取湖岸水涯線信息，對青海湖水面面積進行監(jiān)測；劉曉雪［14］利用多年遙感影像，運用目視解譯和改進的歸一化差異水體指數(shù)方法提取可魯克湖和托素湖的面積，分析其變化特征及降水、氣溫、蒸發(fā)、徑流、農(nóng)業(yè)灌溉耗水等因素的影響；馬艷敏［15］等利用2000年以來的MODIS 數(shù)據(jù)，采用目視解譯的方法提取吉林省西部主要水體分布及面積數(shù)據(jù)，對年際、年內(nèi)變化規(guī)律及影響因素進行了分析；高耶［16］等則選擇了1979-2016年間多時期、多類型、多光譜遙感數(shù)據(jù)，分析評價洞庭湖區(qū)內(nèi)湖近40年的面積變化。

以往南四湖面積監(jiān)測受遙感影像分辨率和獲取時間的限制［17，18］，對不同季度的長時間序列監(jiān)測和年內(nèi)面積變化監(jiān)測相對缺乏。為了研究最新的南四湖水體面積的動態(tài)變化狀況，尋找其面積變化規(guī)律，為恢復(fù)其生態(tài)環(huán)境與湖域面積的相關(guān)措施做出決策參考，本文收集2006-2020年美國陸地衛(wèi)星（Land?sat5、7、8）獲取的95 景遙感影像，進行南四湖水體變化動態(tài)監(jiān)測。

1 研究區(qū)概況

南四湖是山東省最大的天然地表水資源，是南陽、昭陽、獨山和微山4 湖的總稱，位于山東省濟寧市南部，地跨濟寧市、棗莊市、徐州市3 個地級市（圖1）。南四湖具有大量水資源和豐富的濕地資源，屬于水庫型湖泊，濕地保護區(qū)總面積1 266.47 km2，包括南四湖整個湖泊及其灘地、周邊農(nóng)地、山林等。隨著經(jīng)濟社會的發(fā)展，南四湖水質(zhì)逐年惡化，生物群落負(fù)向演替速度加快，富營養(yǎng)化、沼澤化問題嚴(yán)重，水旱災(zāi)害頻發(fā)［19］，大大影響了當(dāng)?shù)厣鷳B(tài)環(huán)境。

圖1 研究區(qū)域位置Fig.1 The location of study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)源

本文從美國地質(zhì)勘探局官網(wǎng)（USGS）收集2006-2020年南四湖地區(qū)（東經(jīng)116°34'～117°24'，北緯34°27～35°20'）Landsat衛(wèi)星（Landsat5 TM、Landsat7 ETM、Landsat8 OLI）95 景影像作為數(shù)據(jù)源進行研究（表1）。

為消除大氣影響，獲得地物真實反射率，影像在ENVI 軟件中經(jīng)過了輻射定標(biāo)和FLAASH 大氣校正等預(yù)處理，ETM 影像需要進行額外的Gapfill處理。同時從歐洲航天局（ESA）官網(wǎng)下載2 景同區(qū)域的Sentinal-2 影像，并利用Sen2Cor 插件對其進行輻射定標(biāo)和大氣校正。該影像具有10 m空間分辨率，用于輔助目視解譯和精度評定。山東省年鑒數(shù)據(jù)來自山東省統(tǒng)計局官網(wǎng)，降雨數(shù)據(jù)來自美國國家海洋和大氣管理局（NOAA）氣候預(yù)測中心（CPC）的全球降雨插值數(shù)據(jù)集。

2.2 方法

2.2.1 水體指數(shù)

利用遙感影像進行水體提取的基本依據(jù)就是水體和其他地物在可見光和近紅外波段存在光譜反射差異［20］。水體指數(shù)是常用的水體自動提取方法，本文分別利用目前廣泛使用的NDVI［21］、NDWI［22］、MNDVI［23］、AWEI_sh［24］、NWI［25］對2017年12月10日的南四湖Landsat8 OLI 影像進行水體信息提取，提取的水體類型包括湖泊、河流、沼澤濕地、坑塘水田、臺田魚塘和人工水渠。各水體指數(shù)計算公式如表2。

表2 典型水體指數(shù)Tab.2 Typical water indexes

式中：ρBLUE、ρGREEN、ρRED、ρNIR、ρSWIR1、ρSWIR2分別表示遙感影像的藍(lán)光、綠光、紅光、近紅外、中紅外和遠(yuǎn)紅外波段。

使用統(tǒng)一的0 作為閾值，對水體指數(shù)圖像進行分割，得到Landsat8 OLI影像基于5種水體指數(shù)的水體分布，結(jié)合時間間隔相差1天的同期Sentinal-2影像輔助目視解譯，將各方法提取結(jié)果進行疊加，放大并比較水體信息差異較大的局部區(qū)域，對比選取最佳的水體指數(shù)用于本文實驗，結(jié)果如圖2。

圖2 不同水體指數(shù)提取結(jié)果Fig.2 Comparison of the extraction results with different water indexes

5 種方法基本可以提取南四湖整體的水體信息，但是在局部地區(qū)，結(jié)合Sentinal-2 影像疊加提取結(jié)果進行目視解譯可以看出，AWEI_sh的提取結(jié)果存在較大范圍的漏分情況，NDWI 和MNDWI 可以較大程度提取水體，但也將部分裸地和旱田錯分成了水體，NDVI 也存在部分漏分情況，并且不能很好地識別水田田埂等非水體信息，而NWI 可以較好的進行水體提取，很少出現(xiàn)漏分和錯分的情況，綜上所述，本文選取NWI 指數(shù)作為南四湖的水體提取方法。

2.2.2 水體面積計算

選取NWI需要的對應(yīng)影像波段進行波段計算，并對NWI提取結(jié)果按0閾值進行影像分割得到水體與非水體，利用ENVI統(tǒng)計水體的像元個數(shù)N，并利用下式計算水體面積：

Area=N R2

式中：R為影像空間分辨率。

春季水體提取結(jié)果如圖3所示。

圖3 南四湖2006年到2020年春季水體提取結(jié)果Fig.3 The water area in Spring in the Nansihu from 2006 to 2020

2.2.3 計算結(jié)果驗證

在能準(zhǔn)確區(qū)分水體和陸地的前提下，分辨率越高，岸線部分水與陸地混合像元中涉及的水體面積就越小，監(jiān)測得到的湖泊面積就約接近實際值［13］。Sentinal-2 MSI 影像空間分辨率為10 m，整體提取效果優(yōu)于Landsat8影像［26］，本文將2017/12/11和2019/11/11 的Sentinal-2 MSI 影像水體提取結(jié)果作為地面參考，對2017/12/10 和2019/11/14 的Landsat8 OLI 影像基于NWI 得到的水體分布進行精度評價，利用ENVI 的Confusion Matrix 工具建立混淆矩陣，計算總體精度和Kappa系數(shù)等定量評價指標(biāo)（表3）。從表3 中可以看出，兩期Landsat8 OLI 影像的總體精度均在90%以上，能夠滿足水體面積提取的精度需求。

表3 水體提取結(jié)果精度指標(biāo)Tab.3 Precision index of water extraction results

3 結(jié)果與分析

3.1 水體年際變化

3.1.1 年際變化

南四湖水面積受季節(jié)影響，為了更全面探究不同季節(jié)南四湖水體面積的年際變化情況，對南四湖水體面積分春（3-5月）、夏（6-8月）、秋（9-11月）、冬（12-2月）四季進行變化分析，結(jié)果如圖4。

從圖4 可以看出，總體上水體面積變化分為兩個階段，2006-2015年，水體面積呈下降趨勢，減少幅度28.768%，年均減少33.854 km2，2015-2020年，水體面積逐漸回升，回升幅度7.682%，年均增加14.488 km2。其中，夏季多年變化最不穩(wěn)定，多年最大面積（1 195.209 km2）和最小面積（503.256 km2）均出現(xiàn)在夏季，部分年份夏季水體面積可能高于當(dāng)年各季度水體面積，如2006年，也可能低于當(dāng)年各季度水體面積，如2014年。其他各季多年變化較為平穩(wěn)，但冬季水體變化幅度最大，2019年較2006年減少面積248.998 km2，變化幅度達(dá)24.284%。

圖4 南四湖15年間水體面積變化趨勢Fig.4 Changes of water surface area in the Nansihu during the past 15 years

各季度具體的年際變化情況如表4 的所示，各季年際變化差異明顯，以春季年際變?yōu)槔?，最大面積為2006年4月16日的1 130.627 km2，最小面積為2014年3月21日的775.345 km2，十五年間面積減少334.123 km2，減少幅度18.264%，年平均減少13.767 km2。

表4 南四湖水體面積年際變化Tab.4 The annual variation of the water area of the Nansihu

3.1.2 影響因素分析

湖泊面積變化跟自然因素和人為因素有關(guān)。對研究區(qū)年降雨總量和南四湖年平均面積進行相關(guān)性分析（圖5），兩者相關(guān)系數(shù)為0.583，并呈現(xiàn)0.05 水平的顯著性，說明南四湖水體面積在2006-2020年間整體上有所減少跟當(dāng)?shù)囟嗄陙斫涤曜兓哂幸欢P(guān)系。同時有研究表明南四湖年平均氣溫每年上升0.24 ℃，蒸發(fā)量增大，一定程度上也減少了湖泊的面積［17］；但水體面積下降至2015年后逐漸回升，這主要得益于南水北調(diào)東線一期工程的人為因素影響，該工程自2013年11月通水后，截至2020年度，分別向山東境內(nèi)調(diào)水1.7、3.28、6.02、8.89、10.88、8.44、7.03 億m3，累計向南四湖、東平湖及南水北調(diào)工程調(diào)蓄水庫進行生態(tài)補水3.74 億m3［27］，為維持和恢復(fù)南四湖水體面積提供了重要水源。

圖5 年降雨和年均水體面積變化Fig.5 Annual rainfall and annual average water area change

需要注意的是，2014年南四湖地區(qū)降雨持續(xù)偏少，特別是入汛以后降水偏少近五成，湖區(qū)蓄水不足歷年同期的二成，湖區(qū)漁業(yè)損失巨大，航運嚴(yán)重受阻生態(tài)瀕臨危機，有關(guān)部門于7月中旬先后開啟“引黃補湖”和“引江補湖”工作，歷時48 d，共向南四湖補水1.45 億m3［28］。所以圖4 中出現(xiàn)了2014年7月11 號水體面積達(dá)到研究期最低值后，當(dāng)年秋季又恢復(fù)到正常水平的現(xiàn)象。

南四湖作為山東省最大的天然地表水資源，是濟寧市地表水的主要來源［29］，濟寧市多年地表水資源的變化可以從側(cè)面反映南四湖水體面積的變化情況。從山東省多年年鑒中選取濟寧市地表水資源數(shù)據(jù)，分析其多年變化情況（圖6）。

圖6 濟寧市地表水資源變化Fig.6 Changes of surface water resources in Jining

從圖6 中可以看出，本文研究期間濟寧市地表水資源的變化情況與南四湖多年水體面積年際變化較為一致，均呈現(xiàn)先下降后回升的趨勢，且兩者的拐點時間（2014年和2015年）也相近，這在一定程度上驗證了南四湖水體面積年際變化規(guī)律的正確性。

3.2 水體年內(nèi)變化

3.2.1 年內(nèi)變化

本文選取數(shù)據(jù)較多的年份制作水體的年內(nèi)變化圖［圖7（a）］，同時統(tǒng)計研究期內(nèi)所有年份能采集的數(shù)據(jù)制作各月水體面積箱線圖［圖7（b）］。從圖7（a）中可以看出，南四湖水體面積表現(xiàn)出明顯的年內(nèi)周期變化規(guī)律：由冬入春時，水體面積逐漸增大，往往在5-6月達(dá)到最大值；入夏后水體面積開始減小，8-9月下降至最低水平；秋冬季水體面積有所回升，但波動面積不大，這也與年際變化分析中，春季的線性回歸線始終高于秋冬兩季的線性回歸線且后兩者基本重疊相吻合。圖7（b）反映的是多年各月水體面積的年內(nèi)分布情況，可以看出，各月中位數(shù)、上下四分位數(shù)的分布規(guī)律與圖7（a）基本保持一致；1、2、5月獲取的數(shù)據(jù)有限以及2014年7月南四湖地區(qū)旱情嚴(yán)重水體面體降至近年最低值，導(dǎo)致了圖中異常值的出現(xiàn)。

圖7 南四湖水體面積年內(nèi)變化Fig.7 Monthly variations of water area of the Nansihu

3.2.2 影響因素分析

湖泊面積反映蓄水量，年內(nèi)蓄水量變化和降雨、湖區(qū)工農(nóng)業(yè)用水量、水利工程建設(shè)等因素有關(guān)［12］。鑒于研究資料獲取的限制，本文僅對年內(nèi)降雨這一因素進行討論?？紤]到夏季能獲取的部分影像往往是在晴天，此時受高溫高蒸發(fā)以及大面積綠潮覆蓋影響（圖8），難以獲得夏季真實水面面積數(shù)據(jù)。在ArcMap中通過人工數(shù)字化綠潮覆蓋面要素并計算其面積，其值占比當(dāng)年平均面積達(dá)23.4%，故在進行相關(guān)分析時，去掉夏季月份的數(shù)據(jù)，以減少異常值的影響。選取可獲得影像較多的年份探究年內(nèi)降雨和水體面積變化的關(guān)系（圖9），計算年內(nèi)月平均降雨和水體面積的相關(guān)系數(shù)，結(jié)果如表5所示。

圖8 夏季綠潮覆蓋Fig.8 Green tide coverage in Summer

圖9 部分年份月均降雨和年內(nèi)水體面積變化Fig.9 Monthly rainfall and annual water area change in some years

從表5 可以看出：2013年以前，降雨量和該年內(nèi)水體面積呈正相關(guān)關(guān)系，均通過0.05的顯著性檢驗，且2010、2013年通過0.01 的顯著性檢驗；而2013年以后的年份，雖然部分年份的正相關(guān)系數(shù)超過0.5，但均未能通過0.05 的顯著性檢驗，說明降雨量和該年內(nèi)水體面積沒有明顯的相關(guān)關(guān)系。這再一次說明在南水北調(diào)東線一期工程2013年11月通水后，除了降雨之外，人工調(diào)水也在一定程度上影響著南四湖的水體面積變化。

表5 部分年份降雨和水體面積相關(guān)系數(shù)Tab.5 Correlation coefficient between rainfall and water area in some years

需要注意的是，在2014年之前，秋冬季的南四湖在夏季高溫蒸發(fā)量大的基礎(chǔ)上面積可能會進一步下降，如圖中2006年的變化情況。但在2013年11月南水北調(diào)東線第一期工程通水后，受人為調(diào)水影響，各年秋冬季的水體面積均為浮動上升變化。

4 結(jié)論

本文收集了2006-2020年美國陸地衛(wèi)星獲取的95 幅遙感影像，利用NWI 水體指數(shù)方法提取南四湖水體面積，分析其在部分年內(nèi)和十五年間的動態(tài)變化過程，主要研究工作和創(chuàng)新有。

（1）采用NWI進行水體面積提?。簩Ρ榷喾N水體指數(shù)，選取最適合南四湖水體提取的NWI水體指數(shù)進行水體提取，采用長時間序列Landsat遙感影像數(shù)據(jù)進行水體面積變化分析；

（2）分季度進行水體面積變化和影響因素的分析：南四湖水體面積呈現(xiàn)先下降后逐年回升的趨勢，2006-2015年間呈下降趨勢，年平均減少33.854 km2，此后轉(zhuǎn)為上升趨勢，年平均增長14.488 km2；各季節(jié)年際變化差異明顯，其中夏季多年變化最不穩(wěn)定，多年最大面積（1 195.209 km2）和最小面積（503.256 km2）均出現(xiàn)在夏季；年內(nèi)表現(xiàn)出明顯的季節(jié)性變化規(guī)律，即由冬入春，水體面積逐漸增大，入夏后開始減小，秋冬季則有所回升，但波動面積不大；2006-2013年間降雨量和該年內(nèi)水體面積呈正相關(guān)關(guān)系，但2013年以后，受南水北調(diào)的引水工程影響，南四湖水體面積變化與降雨關(guān)系減弱；水體面積在2015年開始回升，說明當(dāng)?shù)卣退块T加強了對南四湖的管理，水體面積逐漸恢復(fù)；

（3）采用Sentinal-2 影像和濟寧市地表水?dāng)?shù)據(jù)進行結(jié)果驗證：10 m 高分辨率的Sentinal-2 影像作為真實水體對水體提取結(jié)果進行對比驗證，總體精度在90%以上，能夠滿足水體面積提取的精度需求；濟寧年鑒水資源變化規(guī)律在一定程度上說明了南四湖水體面積變化規(guī)律的正確性。

本研究通過分析近年來南四湖水體面積變化規(guī)律以期為南四湖水資源保護和生態(tài)管理提供科學(xué)支持。但南四湖水體面積的變化受降雨、蒸發(fā)、出入湖水量已經(jīng)取用水量等眾多自然和人為因素影響，為更全面探究南四湖水體面積變化的原因，文章后續(xù)將在這方面進行深入研究。□