紅崖山水庫近20年面積變化遙感調查及驅動力分析

郝固狀，甘甫平，閆柏琨，李賢慶，胡輝東

(1.中國礦業大學(北京)煤炭資源與安全開采國家重點實驗室，北京 100083；2.中國礦業大學(北京)地球科學與測繪工程學院，北京 100083；3.中國自然資源航空物探遙感中心，北京 100083)

0 引言

水資源是人類生存和社會發展不可缺少的生命之源和物質基礎，而水庫作為水資源中不可或缺的淡水資源，為人類的生產、生活提供了重要的保障。水庫的分布狀況在一定程度上反映了水資源的區域變化和時空差異。水庫水量變化在工業、農業、城市以及防汛抗旱等方面發揮著重要的調節作用[1]，而水庫的水域面積變化直接反映著水庫水量的增減狀況。紅崖山水庫是位于石羊河流域的大型水庫，是上游冰川融水、降水水源補給的主要匯集地與下游工、農、生活用水的主要補給與調節地，其水量、水面分布變化受氣候、人為因素的雙重作用與影響。

隨著遙感技術的不斷成熟，充分發揮遙感技術宏觀性、動態性和實時性的優勢，通過遙感技術可以提供大范圍、長時間序列的動態監測，可以節省大量的人力、物力、財力，并快速準確地獲取湖泊和水庫面積信息。李均力等[2]通過對柴窩湖面積近50 a的時序變化進行分析，得出面積變化分為3個階段，面積呈現先緩慢增加后快速退縮的趨勢；魏善蓉等[3]利用MODIS09衛星數據通過對比6種水體指數法的湖泊提取結果，對柴達木盆地連續14 a的湖泊面積進行解譯，湖泊呈現先增加后減少的趨勢變化；張文等[4]利用Landsat和GF影像對鄱陽湖進行了長時間序列監測，研究結果表明在空間上，鄱陽湖水體面積呈緩慢收縮的趨勢，在時間上，除冬季面積有明顯下降趨勢外，其他季節整體趨勢變化不大，不同季節鄱陽湖的湖面范圍有明顯的差異，湖泊面積具有顯著的季節性；張兆鵬等[5]以95景多時相Landsat衛星影像為主，還原了1987—2016年年際變化(以8月水面分布為對比基準)及部分年內月度變化，并定性分析了水量變化的自然、人為因素影響，認為年際水面變化呈“W”形，總體為增加趨勢，人工增雨增雪、農業灌溉、水庫建設是主要的影響因素。該研究對紅崖山水庫的水面面積變化及驅動力得出了基本與初步結論。

本文利用谷歌地球引擎(Google Earth Engine，GEE)平臺進行水域面積的提取，這是一種更為便捷快速的研究方法；利用更高分辨率的遙感影像來驗證真實的水庫水域面積變化，并進行精度對比，提取出更加真實的水域面積；利用氣象數據，綜合氣溫、降水、蒸散發的數據以及多角度影響因素對變化原因進行分析。在全球氣候顯著變化的背景下，還原監測其水量、水面歷史分布變化，分析解剖變化的驅動力，對于研究我國西北內陸河流域和水庫水資源變化、保障、氣候變化影響與適應對策等具有重要的代表意義。

1 研究區概況與數據源

1.1 研究區概況

亞洲最大的人工沙漠水庫——紅崖山水庫位于西北地區甘肅省民勤縣，被騰格里和巴丹吉林兩大沙漠包圍，屬于典型沙漠平原水庫,因位于石羊河下游，其水源補給主要來自石羊河流域(圖1)，紅崖山水庫是沙漠地區的一座中型洼地蓄水工程，也是唯一的水利調蓄工程，西面依紅崖山而建，其他三面都是人工所筑，而且又修建在沙漠中，這在全國甚至全世界都是罕見的[6-7]。

圖1 研究區地理位置Fig.1 Location of the study area

1.2 數據源

1.2.1 遙感數據及預處理

本研究主要采用的是Landsat系列影像數據(數據來源：https://developers.google.com/earth-engine/datasets/catalog/landsat)和高分二號(GF-2)衛星影像(數據來源：http://geocloud.cgs.gov.cn)，其中Landsat包括Landsat5 TM，Landsat7 ETM+和Landsat8 OLI數據，影像獲取的時間是2000—2019年，空間分辨率為30 m，其中Landsat5有254景，Landsat7有34景，Landsat8有171景(表1)，GF-2包括多光譜(PMS)和全色(PAN)，分辨率分別為4 m和1 m，共12景(表2)。本文對Landsat5/7/8的表面反射率(surface reflectance，SR)數據進行了幾何糾正、云雪去除、影像拼接、影像融合等預處理，為后續水體的提取奠定了基礎；同時對GF-2數據進行了輻射定標、大氣校正、正射校正、影像配準、影像融合和影像裁剪等預處理。本文采用Landsat影像數據與采集時間基本一致的6景GF-2衛星遙感影像作水體提取精度對比。

表1 Landsat遙感影像數據Tab.1 Landsat remote sensing image data

表2 高分二號(GF-2)遙感影像數據Tab.2 GF-2 remote sensing image data

1.2.2 氣象數據

獲取研究區周邊5個氣象站點的2000—2018年的氣溫、降水數據(數據來源：中國氣象數據網http://data.cma.cn)。本次研究的氣象站點的數量有限，為獲取區域內各點的氣溫降水值，采用克里金插值法得到氣象要素柵格圖[8]，并提取紅崖山水庫質心處氣溫降水數據。

2 研究方法

本文紅崖山水庫水域面積的提取方法，精度驗證以及氣象因子、植被覆蓋度和人工自然等其他因素對水域面積變化驅動力的分析流程如圖2所示。

圖2 總體技術路線Fig.2 Overall technical route

2.1 水體指數法提取水體

基于遙感影像數據的水體提取方法主要包括光譜特性提取和紋理特征提取2類。基于光譜信息的提取方法有單波段閾值法[9]、光譜指數法、分類器法；基于紋理信息的提取方法主要是水體特殊的紋理。本文Landsat遙感數據主要對改進歸一化差異水體指數(modified normalized difference water index，MNDWI)、歸一化差異水體指數(normalized difference water index，NDWI)、增強型水體指數(enhance water index，EWI)、新水體指數(new water index，NWI)4種光譜指數法進行對比分析。

運用水體指數方法的關鍵在于閾值的選擇，本文選用Ostu算法[10]進行閾值的自動分割。Ostu是一種確定圖像分割閾值的算法，按其求得的閾值進行圖像二值化分割后，前景與背景圖像的類間方差最大；其被認為是圖像分割中閾值選取的最佳算法，計算簡單，不受圖像亮度和對比度的影響，因此在數字圖像處理上得到了廣泛的應用。

1)McFeeters[9]提出的NDWI主要是根據水體的光譜特性利用水體在綠波段(Green)和近紅外波段(NIR)的反射特性，使得影像中的水體和其他地物區分開來，同時抑制植被和土壤信息。其計算公式如下：

(1)

式中：Green為TM/ETM+/OLI遙感數據中綠光波段數據；NIR為近紅外波段數據。

2)MNDWI[11]水體指數是對NDWI指數的進一步修正，削弱了提取水體時由于山體陰影的影響，可有效抑制建筑物信息，能夠更好地提取城市范圍內的水體信息。其計算公式如下：

(2)

式中：MIR為TM/ETM/OLI遙感數據中紅外波段數據。

3)在GEE平臺使用NWI[12]指數時，Landsat數據的表面反射率(surface reflectance,SR)和大氣表觀反射率(top of atmosphere,TOA)產品的水體提取效果是不同的，經過試驗，NWI(TOA)只要大于0就可以很好地提取出水體，而NWI(SR)則要自己調整閾值，根據普適性，NWI(SR)的組合更能滿足多地區的水體提取。其計算公式如下：

(3)

式中：B1,B4,B5,B7分別對應TM/ETM+遙感數據中的第1,4,5,7波段數據(OLI對應遙感數據中的第2,5,6,7波段)；C為常數，本文中取200，也可取100或其他正值。

4)EWI[13]主要是將MNDWI和NDWI相結合，利用綠波段，紅外波段和短波紅外波段進行水體提取，其計算公式如下：

(4)

2.2 植被覆蓋度的計算

歸一化植被指數(normalized difference vegetation index,NDVI)計算公式如下：

NDVI=(NIR-R)/(NIR+R)。

(5)

NDVI與植被覆蓋度之間有緊密相關性，綜合大量研究，把NDVI與像元二分模型結合起來，建立了基于NDVI的植被覆蓋度估算模型[14]，計算公式如下：

(6)

式中：NDVI為混合像元的值；NDVIveg為純植被像元的值；NDVIsoil為純非植被像元的值，理論上NDVIveg=1，NDVIsoil=0。

由于在現實情況中，植被覆蓋受到多種氣象因子的影響，一般不可能取到理論值，又因為缺乏大面積地表實測數據，經大量研究取定置信區間內的NDVI最大值(NDVImax)和最小值(NDVImin)來代替[15]，因此公式變為：

(6)

本文通過此種方法，將2000—2019年中偶數年份的NDVI值進行了統計，根據累積百分比選擇累積百分數2%～98%為置信區間，累積百分比小于2%為近似純非植被覆蓋，小于98%為純植被覆蓋，分別將對應的NDVI值代入NDVImin和NDVImax進行計算[16]。

3 結果分析

3.1 GF-2/Landsat影像水體面積提取的精度驗證

由于Landsat數據為低分辨率遙感影像，為保證提取的水體精度的可信度，故采用空間分辨率為1 m的GF-2影像數據作為精度驗證數據。通過4種水體指數與GF-2影像目視解譯法提取的水域面積進行對比分析(表3)，可知MDNWI水體指數法提取的水體面積精度最高，絕對誤差僅為0.006 km2，相對誤差0.029%，結合4種水體指數提取的效果圖，MNDWI指數(圖3(b))提取的水體更加符合原始的遙感影像(圖3(a))中的水體，邊界上的提取更加全面真實。NDWI指數(圖3(c))和NWI指數(圖3(d))在石羊河下游與水庫入口的連接處的水體未提取出來，EWI指數(圖3(e))提取效果總體上很符合真實情況，但是在水庫邊界的提取上有部分水體無法提取，因此綜合對比分析4種指數(圖3(f))的提取效果，本文將采用MNDWI指數法作為水體提取的方法。

(a)Landsat8原始影像 (b)MNDWI水體指數法 (c)NDWI水體指數法

表3 4種水體指數法與GF-2目視解譯法水體面積提取精度對比Tab.3 Comparison of water area extraction accuracy between four water index methods and GF-2 visual interpretation method

由于GF-2衛星影像數據缺少中紅外波段，采用GF-2 PMS傳感器獲取的4 m多光譜數據和1 m全色數據進行圖像融合，得到1 m空間分辨率的多光譜數據，由于波段限制，GF-2影像采用NDWI水體指數法進行水體的提取(圖4)。根據表4和圖5的精度對比，Landsat影像數據和GF-2影像數據的相關系數高到0.99，相對誤差最大為6.63%，最小僅為0.33%，絕對誤差最大不超過1 km2，明顯看出不同分辨率的遙感影像數據對紅崖山水庫的水域面積的提取影響不大，因此Landsat影像數據對水域面積的提取可以反映紅崖山水庫的真實性。

(a)GF-2原始影像 (b)GF-2二值影像 (c)GF-2水體和非水體

表4 GF-2/Landsat數據水域面積提取精度對比Tab.4 Comparison of precision of water area extraction in GF-2/Landsat data

圖5 GF-2/Landsat水域面積相關性Fig.5 GF-2/Landsat water area correlation

3.2 水庫水域面積變化特征

3.2.1 月際面積變化特征

水庫的水量和湖泊的水量變化十分相似，受季節性影響也存在枯水期和豐水期，本文選取每4年間隔(2000年、2004年、2008年、2012年、2016年、2019年)的年份進行月際之間的對比分析。綜合分析6 a的柱狀圖和趨勢線(圖6)，結果表明：2000年、2004年、2008年、2016年和2019年5 a的水域面積月際趨勢呈倒“正態曲線”分布，且每年出現兩次峰值，主要集中在3月份和9—10月份，最低值出現在6月份，根據豐水期和枯水期的劃分，2000年、2004年、2008年、2016年和2019年的豐水期為春秋季節的3月份和9—10月份，枯水期為夏季6月份；但2012年較為特殊，主要表現在枯水期為4月份和7月份，豐水期仍然是春秋季節的3月份和9—10月份。總體上分析2000—2019年月均水域面積變化(圖7)，月際面積變化呈倒“正態曲線”分布，且豐水期主要集中在春秋季節的3月份和9—10月份，枯水期主要集中在夏季6月份。

(a)2000年 (b)2004年 (c)2008年

圖7 2000—2019年紅崖山水庫月均水域面積Fig.7 Average monthly water area of Hongyashan Reservoir from 2000 to 2019

3.2.2 年際面積變化特征

由于水庫不同季節的水量有很大變化，本研究選取不同季節年際之間的水域面積變化進行分析。春季選取4月份，夏季選取7月份，秋季選取10月份，冬季選取12月份進行年際之間水域面積的比較。根據不同季節年際變化情況的研究(圖8)，結果表明：春夏秋冬四季水域面積都是逐年上升的，春秋季節平均年際增長率為5.03%和5.22%，秋季平均年際增長率最低僅為2.42%，夏季水域面積平均年際增長率為22.19%，年際變化幅度最大，秋季變化幅度最小，因此夏季水域面積變化最不穩定，秋季水域面積變化最穩定。夏季水域面積呈“V”字波動上升且變化幅度很大。

(a)春季2000—2019年面積變化 (b)夏季2000—2019年面積變化

3.2.3 水庫水域面積變化總體特征

通過上述年際變化得出夏季的增長趨勢與2000—2019年的總體上(圖9)的增長趨勢極為相似，也呈現“V”形波動上升。根據圖9的研究結果表明，總體上水域面積年際增長的相關系數達到0.78，相關性良好；2000—2019年水域總面積增加8.98 km2，面積變化率高達42.6%，水域面積經歷了先增加再減少，后持續波動上升的過程，具體表現為：2000—2001年水域面積減少，且在2001年達到了最低值10.7 km2。2001—2003年水域面積增長迅速，面積變化速率為27.7%，為近20 a水域面積增長幅度最大。2003—2004年水域面積驟降4.33 km2，為2000—2019年降幅最大。2004—2017年水域面積呈“V”字波動上升，總體面積增加6.97 km2，增長率為56.6%，其中2004—2005年,2014—2015年為小幅增長；2005—2006年水域面積變化速率為27.4%，接近最大變化速率；2006—2010年，2010—2012年，2012年—2014年,2015—2017年水域面積分別為3個“V”形增長；2017—2019年水域面積近乎直線增長，最接近趨勢線走勢。根據年際變化率(圖10)可以看出年際之間的變化幅度在逐年下降且趨向于0，說明紅崖山水庫的水域面積趨向于穩定增長趨勢。

圖9 紅崖山水庫年均水域面積變化Fig.9 Annual average water area change of Hongyashan Reservoir

圖10 紅崖山水庫水域面積年際變化率Fig.10 Interannual change rate of water area of Hongyashan Reservoir

3.3 水域面積變化驅動力

3.3.1 氣象因子的影響

紅崖山水庫位于西北干旱地區甘肅省民勤縣境內，是最大的人工沙漠水庫，選取民勤縣周邊5個氣象站點的氣溫和降水數據進行研究，分析影響水域面積變化的驅動因子(圖11)。

1)氣溫。紅崖山水庫的主要水源是祁連山東部的冰雪融水[7]，在降水量不變的情況下，氣溫的升高會導致冰雪融化，入庫水量增加，水域面積增大。根據對近20 a紅崖山水庫的年均氣溫(圖11(a))研究，結果表明：2000—2005年均低于距平氣溫6.72 ℃，2006年氣溫驟升至7.17 ℃，水域面積增加了3.43 km2；2006—2008年氣溫連續下降了0.89 ℃，此時水域面積也在持續下降，2008—2009年氣溫上升0.8 ℃，2009—2012年氣溫降至最低氣溫6.01 ℃，2013年氣溫達到最高7.26 ℃，2014—2018氣溫波動上升，年均氣溫變化總體來說，研究區氣溫呈上升趨勢，并通過氣溫和水域面積相關性分析(圖11(b))得出，氣溫對水域面積的增加起到正向作用，因此氣溫上升是水域面積增加的重要的間接驅動因素。

(a)年平均氣溫距平變化 (b)水域面積-年均氣溫相關性

2)降水。西北干旱區降水稀少且分布不均，但卻是水庫水量增加的最直接的補給水源。根據近20年降水量的年際變化(圖11(c))表明：2000—2006年(除2003年)降水量低于距平降水量，2003年、2007年和2018年降水量明顯高于相鄰年份，其對應的年份水域面積也明顯增加，2008—2012年降水量波動增加，2012—2015年出現2次明顯下降，一次上升，同年的水域面積也發生同樣的變化，2015—2018年降水量連續4年上升，總體上看，降水量整體呈明顯增加態勢，通過水域面積與降水量進行的相關性分析，相關系數高達0.93(圖11(d))，進一步說明二者有很好的相關性，因此降水量的增加是影響水域面積增加的主要驅動因素。

3.3.2 植被蓋度的影響

植被不僅對區域氣候變化有著重要的影響，植被面積的增加可以減少地表對太陽短波輻射的反射率，緩和地層增溫，進而減少地氣湍流交換和水分蒸發，同時植被恢復也有利于地表水源涵養功能的增加，減少降水的無效下滲和散失，有助于徑流的產生，增加水量。通過對研究區周邊的植被指數提取，根據年際之間植被覆蓋度(圖12)顯示，總體上植被覆蓋度在逐漸增加，2000—2002年植被蓋度增加了10%，但是2002—2008年出現持續下降15%，2008年達到最低值，2008—2018年植被蓋度逐年上升，增加26%。通過水域面積和植被覆蓋度的相關性分析(圖13)，表明植被覆蓋度對水域面積的增加起了正向作用，植被蓋度的增加是導致水域面積增加的重要的間接驅動力。

圖12 植被覆蓋度的年際變化Fig.12 Interannual changes in vegetation coverage

圖13 水域面積-植被覆蓋度相關性Fig.13 Water area-vegetation coverage correlation

3.3.3 其他因素的影響

1)加高擴建工程。為緩解紅崖山水庫防洪與興利，來水與需水的主要矛盾，實施了紅崖山水庫加高擴建工程，庫區水域總面積約22.0 km2，和通過遙感影像獲取的最大面積基本吻合，同時使得總庫容從0.993億m3增加到1.48億m3[17-18]，不僅使得西北干旱地區的水源得到充足使用，而且提高了水庫蓄水的調節能力，加強了農業灌溉的力度，這是紅崖山水庫水域面積逐漸增大的主要驅動力。

2)生態環境修復。為了石羊河流域的綜合治理，恢復青土湖生態，紅崖山水庫向青土湖的下泄生態水量逐年增大[19]，使得青土湖水量從無到有，水位有所上升，地表植被覆蓋逐漸增加，這是紅崖山水庫水量逐年增加的重要原因。每年夏天農業灌溉結束之后，紅崖山水庫將在8—11月份有計劃地向青土湖生態補水3 180 m3[20]，這是導致每年11月份紅崖山水庫水域面積變小的重要原因，因需求量逐漸增加，紅崖山水庫的水量也在不斷增加，這是導致紅崖山水庫逐漸上升的間接驅動力。

3)入庫徑流變化。1962—2010年均入庫徑流最小值為2001年的0.70億m3，據金彥兆等[21]研究，紅崖山水庫站徑流量年際變化曲線表明，2001—2010年入庫徑流量逐漸增加，根據紅崖山水庫站季節徑流量變化圖發現，春夏秋冬四季的徑流量自2000年以后都出現回升，近20 a紅崖山水庫水域面積不斷增加的直接原因是入庫徑流量自2000年逐漸增加。

4)工、農、生活用水。紅崖山水庫的最重要的職能是供工業園區、農業灌溉和百姓生活用水。民勤縣紅沙崗工業園區的需水主要來自紅崖山水庫，以滿足工業用水，保證工業生產順利進行，為民勤縣的經濟發展提供便利[22]，為當地的國民經濟打下堅實基礎。每年的夏季6—8月份為農業、林業灌溉期，氣溫升高，蒸散發增加，導致該時期紅崖山水庫出現枯水期，紅崖山擔負了紅崖山灌區13個鄉鎮，2個國營農林場農田及生態經濟林草灌溉任務,實際灌溉面積440.4 km2，是主要的灌溉水來源。由于工農業用水的增加導致水域面積處于波動上升的趨勢。2010年全縣總人口比2000年減少近2.8萬人，人口總量與用水呈正相關，人口的減少意味著用水量的減少，這是水域面積增加的間接驅動力。

4 結論

本文通過對西北地區紅崖山水庫2000—2019年的水面面積變化遙感調查及驅動力分析，得出以下結論：

1)與傳統的遙感監測方法相比，利用GEE云計算平臺可以在線處理多種類型的遙感影像，節省存儲和運行空間，不受時空的制約，處理速度得到數量級提升，效率得到指數提高，為后續的大尺度區域的遙感影像處理提供了有力支撐。

2)紅崖山水庫的水域總面積增加8.98 km2，面積變化率高達42.6%，呈波動上升的趨勢，年際變化率逐年減小，水域面積整體趨于穩定增加狀態。2000—2019年水庫水域面積月際變化呈倒“正態曲線”，豐水期主要集中在每年春秋季節的3月份和9—10月份，枯水期主要集中在每年夏季的6月份；春夏秋冬四季年際水域面積都是逐年上升，春秋季節平均年際增長率為5.03%和5.22%，秋季平均年際增長率最低僅為2.42%，夏季水域面積平均年際增長率為22.19%，年際變化幅度最大，秋季變化幅度最小，因此夏季水域面積變化最不穩定，秋季水域面積變化最穩定，夏季水域面積呈“V”字波動上升且變化幅度很大。

3)根據對紅崖山水庫水域面積增加的驅動力分析，直接驅動因素為降水量的變化，加高擴建工程，以及入庫徑流變化，其中降水量的增加和加高擴建工程是水域面積增加的主要驅動力；間接驅動因素為氣溫的變化，植被覆蓋度的變化，生態環境修復以及工、農、生活用水，都是十分重要的驅動力。