近50 a薩雷茲堰塞湖水域時序變化及其驅動因素研究

洪永欣,鄧晚倩,張 新,

(1. 河北工程大學 地球科學與工程學院, 河北 邯鄲 056038; 2. 中國科學院大學, 北京 100049;3. 中國科學院空天信息創新研究院 遙感科學國家重點實驗室, 北京 100101)

0 引言

薩雷茲堰塞湖位于亞洲中部的帕米爾高原上,是世界上最深的堰塞湖,于1911年因地震致使山體大面積崩塌而形成。形成至今,薩雷茲湖水位呈波動性顯著上漲趨勢[1],烏索伊大壩壩體一旦崩解將會給壩體下游中亞人民帶來巨大災難[2]。近年來,薩雷茲湖所在地區引起了國際社會的高度關注[3],加強災害的識別與評估,可有效提高災害風險防治的水平[4]。堰塞湖作為洪水災害的潛在因素,已有多學科學者對其形態[5]、穩定性[6]、崩潰預測[7]做相關研究;對于薩雷茲堰塞湖研究的主要集中在壩體監測預警[8-10]、水量監測[1]、風險評估[9]、防災減災[3,11]、地震對薩雷茲湖的影響[12-15],洪水演算模擬[2]等方面。因此,監測薩雷茲湖水資源動態變化,研究湖泊水位變化特征規律對于薩雷茲湖的綜合治理、水量調控等具有重要的指導作用以及科學價值。

國內外學者對薩雷茲湖水資源多從成因災害評估、潰壩模擬、影響等角度對其進行了闡述,而對薩雷茲湖水資源水位水量積蓄過程中的變化特點和發展趨勢、驅動因子分析等相關的內容研究還比較缺乏。然而了解薩雷茲堰塞湖積蓄過程,及影響該過程的相關因素,可以達到直接和間接的觀測效果,實現預測薩雷茲湖水位水量趨勢、潰壩風險,最大程度減少或者避免損失。歸一化水體指數NDWI(normalized water index, NDWI)是遙感水體最為有效的提取方法之一[16]。MCFEETERS等[17]發現利用NDWI方法在近紅外波段水體吸收最強,而植被和土壤在近紅外波段反射率較高,李飛等[18]在研究復雜水體信息提取方法研究中表示,NDWI能夠很大程度上去除植被信息,強調水體。已有的圖像分割方法有:基于區域[19]、邊緣[20]、學習分類[21]和閾值[22]的方法等,其中采用閾值分割對研究區內水體遙感影像信息進行提取的方法是最為廣泛,呂爭等[23]提取水體時,在初提取階段靈活運用全局閾值分割法中的最大類間方差算法,在精處理時采用局部閾值分割法,全局-部分相結合精度高,工程量小。文中遂以薩雷茲堰塞湖為研究對象,從區域、季節等多個角度對薩雷茲堰塞湖水域水量變化進行監測,并結合遙感數據以及氣溫、降水、蒸散發等氣象因素,對薩雷茲湖水資源變化的驅動因素進行探討。

1 研究區域及數據來源

1.1 研究區概況

薩雷茲湖位于塔吉克斯坦的帕米爾高原,位于72°32′30″～73°11′21″E,38°9′21″～38°20′20″N,薩雷茲湖湖區海拔高度約為3 263 m,湖集水面積16 506 km2,薩雷茲湖年平均長度約 55.8 km,平均寬度為 1.44 km;流域整體的平均海拔為3 263 m,堰塞湖總體積約17 km3,最深點在500 m,平均深度為202 m。巴爾坦格河為薩雷茲湖的下游,穆爾加布河為上游。薩雷茲湖年平均環境氣溫為1 ℃,近年來年平均降水量為108 mm。

1.2 數據資料

1.2.1 Landsat數據

Landsat是美國航天局和美國地質調查局共同管理維護的對地衛星,主要用于探究地下礦物資源、海洋能源和地下水資源,研究農作物生長規律、預測農作物產量,調查預測自然災害,制作專題地圖等。考慮時間節點、云層干擾、湖面結冰情況,最終選取了1972—2019年間的Landsat分辨率多光譜數據,用于分析薩雷茲湖面積變化長時間序列。

1.2.2 氣溫與降水數據

氣溫數據來源于全球歷史氣候學網絡數據庫(www.ncdc.noaa.gov/cdo-web)。文中選擇了距離薩雷茲湖最近、數據覆蓋最完整的霍羅格氣象臺站來研究氣溫數據。降水數據來自中亞生態與環境研究中心,采用了薩雷茲湖水文監測站點的1972—2018年薩雷茲湖日觀測水位數據以及2005—2019年薩雷茲湖入湖流量日觀測數據,用以分析薩雷茲湖水資源變化及其變化的驅動因素。

1.2.3 潛在蒸散

潛在蒸散數據來自東英吉利大學氣候研究格網單元與英國國家大氣研究中心聯合發布的1901—2019年期間全球陸面月平均數據集(CRU-TS-4.04)。其空間分辨率為0.5°×0.5°(https://catalogue.ceda.ac.uk/)。文中采用的是薩雷茲湖湖面1972—2019年間潛在蒸散量的月平均值。

2 研究方法

綜合考慮影像的時間節點、云層覆蓋以及湖面結冰狀況,以及薩雷茲湖周圍山體眾多,存在影像被干擾的情況,選用了歸一化水體指數與全局-局部閾值分割技術相結合的自適應水體提取算法[24]對薩雷茲湖的湖岸廓線進行提取。

2.1 歸一化水體指數提取

本研究采用的是歸一化水體指數NDWI,可以實現有效地增強水體光譜特征并抑制植被等環境信息。MCFEETERS等[17]給出的NDWI公式如式(1)所示:

(1)

式中:ρGreen為遙感影像波段中大氣反射率的綠波段,對應Landsat的MSS中的波段1;ρNIR為遙感影像中表觀反射率的近紅外波段,對應于Landsat中的波段3,在TM/ETM+中對應的波段為波段2與波段4,因此可以靈活地獲取長時間序列的提取。由于不同傳感器的標定系數和太陽光譜輻照度不同,NDWI的計算是基于大氣表觀反射率。

2.2 全局-局部分割閾值選擇

通過對NDWI的物理特征的了解,從中選取全局閾值,其閾值的范圍通常為[-0.1,0.1]。第1步,全局分割整幅圖像。首先通過設定初始全局閾值(如T0=-0.1)對影像進行分割,對分割結果向外擴充若干單元設立相應的緩沖區,其中緩沖區的像元總數與湖泊像元總數相同,得到初分割湖泊。第2步,針對薩雷茲湖的特征作局部分割,其周圍的背景像素作為局部分析單元。將上述包含緩沖區的局部單元進行直方圖統計,若頻率分布直方圖呈雙峰分布,則通過式(2)閾值將湖泊與緩沖區進一步分割,得到新湖岸廓線;對該湖泊矢量建立新的緩沖區,按照下述閾值進行再次分割迭代,至2次湖泊矢量的像元個數差值小于給定閾值。

(2)

式中:μwater為水體的平均值;μland為陸地像元的平均值;σwater、σland分別為水體與陸地像元的標準差。

第3步,混合像素。由于山體陰影與湖泊水體具有相似的光譜特征,因此僅通過影像中的光譜信息難以對其進行有效的區分。通過DEM數據生成的坡度圖進行分割,將坡度圖與湖泊矢量進行疊加,剔除山體陰影的像素點。最后,為保證提取矢量圖層的高精度,針對每期提取后的結果進行目視判讀,迭代檢驗矯正,最終可以得出薩雷茲湖矢量數據。

2.3 精度檢驗

對于遙感數據預處理后所得薩雷茲湖岸線精度,本研究聚焦于薩雷茲湖的動態變化,因此采用目視判讀的方法。憑借實踐經驗和專業知識在薩雷茲湖DEM遙感影像上獲取地物特征,實踐表明判讀的結果滿足薩雷茲湖水體研究動態變化的條件。在分析水文參數趨勢變化檢驗方面,使用Mann-Kendall非參數統計法[25-26](M-K)來研究薩雷茲湖動態變化趨勢。文中統計了水域各季節的平均面積變化,運用M-K以四季為單位,統計構建得到薩雷茲堰塞湖水體面積的變化趨勢與突變時間拐點,具體計算方法如式(3),針對水體面積時間序列:x1,x2,x3,…,xn,構造一個秩序列ri,ri表示xi>xj(1≤j≤i)的樣本累積數。定義Sk:

(3)

式中,Sk的均值E(Sk)及方差var(Sk)定義如式(4):

(4)

水域面積時間序列是隨機獨立分布為假設的前提條件,那么定義統計量UF表示為:

(5)

式中,UFk為標準正態分布,其顯著性水平為0.05時的臨界值UFα=±1.96,當|UFk|>UFα時,測得的薩雷茲湖水域面積序列表現出顯著的增長趨勢或減少的趨勢;其中,把UFk繪構為曲線c1。

以上述的計算方法作用于反序列上,迭代計算,并將結果乘以-1,得到UBk;所有的UBk在圖中構成曲線c2。通過繪制c1、c2曲線圖,若UFk或UBk的值大于0,則表明薩雷茲湖水域面積時間序列中的該時段上為上升趨勢,反之則為下降;當曲線超過臨界值時,表明該時間段內存在劇烈的上升或下降幅度;若c1、c2兩曲線交點在臨界線區間內,表示在該交點處,為突變點的起始時間可能性大。

3 結果與分析

3.1 薩雷茲湖水域面積時序變化

經過對Landsat遙感影像數據進行NWDI和全局-局部分割閾值處理后的薩雷茲湖水域面積時序變化進行具體分析,如圖1所示。在1972—1998年間,水域面積呈緩慢上升趨勢,上升速率達到0.189 km2/a;1999—2004年間,水域面積呈緩慢下降趨勢,下降速率為-0.293 km2/a;2005年,薩雷茲湖水域面積急速擴張,上升速率高達到4.038 km2/a; 2006—2009年水域面積呈現出急速下降的趨勢,下降速率為-0.410 km2/a;2010—2019年水域面積始終保持周期性上升,上升速率為0.221 km2/a;將薩雷茲湖的面積變化時間序列按照逐年季節進行平均計算,春夏秋季節薩雷茲湖水域面積整體變化呈振蕩上升趨勢,年均的變化速率略有不同,春季水域面積的年平均增長率為0.157 km2/a,為一年中增長最快的季節;秋季在1972—1988年之間變化的水域面積增加幅度高達4.77 km2;四季的水域面積變化率,春夏秋按順序逐漸減少,冬季保持波動性平穩。

圖1 薩雷茲湖時間-面積變化序列圖Fig. 1 Time-area change sequences of Sarez Lake

將薩雷茲湖水域面積時序進行年平均結合Mann-Kendall非參數檢驗分析,如圖2。由于數據缺失的原因,暫不考慮1976年,近40 a 的UF統計量的值均為正值,發現在1991—1993年這3 a之間的UF統計量在臨界值之上,其顯著水平均大于0.05臨界值,則通過具體分析得到薩雷茲湖水表面積具有顯著的上升趨勢;具有同樣突出特征的年份有:1996—2001年以及2013年至今。然而,在2012—2013年,UF與UB兩統計量曲線有交點,則表明在該年份間薩雷茲湖湖面的年均面積變化劇烈。

圖2 薩雷茲湖年平均水域面積突變檢驗分析Fig. 2 Mutation test of annual average water areas of Sarez Lake

在春夏秋冬四季中水域面積的變化有所不同,經M-K突變檢測得春夏兩季表現規律相似,UF統計量的值始終為正值。2006年以后春季平均水域面積上升趨勢顯著,2012年以后水域面積具有顯著上升特征的季節為夏季,其UF統計量表現為在0.05臨界值線之上,此外春夏兩季的臨界線之間的交點分別位于2005年與2011年;秋季在1972—1983年UF統計量為負值,但UB統計量仍為正值,水體表面積始終處于上升趨勢,其余年份UF統計值均為正值,自2015年以來,其UF統計值均超過0.05顯著水平的臨界點,說明水域面積有明顯上升趨勢,水域面積變化的突變點與夏季相近,位于2011—2012年間;因此,1997—2000年及2003—2013年水域面積呈上升趨勢,且上升趨勢較平緩,在1972—2019年間上升下降趨勢變化幅度較大。因此,薩雷茲湖水域面積變化在2010—2013年間存在突變拐點,且變化的趨勢成劇烈狀態。

3.2 薩雷茲湖水域體量變化

本研究獲取了薩雷茲湖的實測站點水位數據,獲取了1972—1991年、2005—2011年、2014—2018年3個時段水位高程,分析水位在該3段時間序列里的變化,得出1971—2018年水位總體的趨勢為波動性上升。繼而,根據現有的觀測數據,主要包括5個時間段的結果:1972—1980年、1981—1985年、1986—1991年、2005—2011年、2014—2018年,表現出來的變化規律略有區別,對這5個階段的水位時間序列分別建立時間-水位變化函數,整體上1972—2018年間水位呈0.147 m/a的速率周期性上升,在1972—1980年間水位上升速率為0.560 m/a,而1981—1985年水位快速上漲速率達到0.784 m/a,在1986—1991年間水位呈0.356 m/a的波動型緩慢上漲;2005—2011年為唯一時間段水位呈周期性下降期速率達到-0.473 m/a,隨后在2014—2018年,水位再度呈現上漲趨勢速率達到0.444 m/a。為研究薩雷茲湖水域面積與同期水位之間的相關性,建立了兩者的關系曲線,通過時間匹配選取了133組面積水位數據點,二者顯示出了良好的相關性,并建立面積-水位關系模型:

S=-0.036 55×dh2+4.762 07×dh-68.717 94

(6)

式中:dh為相對水位(m),dh=h-3 200;S為水域面積(km2)。

根據湖泊水位時間序列與面積時間序列的計算,可以按梯形臺錐臺體積式(7)估算湖泊的體積變化,從而獲得薩雷茲湖體積變化:

(7)

式中:H1為最低水位高度;H2為最高水位高度;A1為高水位時湖泊的面積;A2為低水位時湖泊的面積。

由于水域面積與水位序列的時間匹配不佳的問題,文中通過面積-水位關系模型獲取與水位數據相匹配的面積序列,并結合水位時間序列來進行湖泊體積變化的計算。1972—2018年薩雷茲湖量總體表現為波動增長的趨勢,其中主要包括6個不同時間段的結果:1972—1980年、1980—1983年、1983—1985年、1985—1991年、2005—2011年、2014—2018年,對這6個階段的水量時間序列分別建立時間-水量變化函數,整體上1972—2018年間水量呈0.012 km3/a的速率周期性上升。

3.3 薩雷茲湖水域變化驅動力因素

湖泊水量變化的驅動因素眾多,比如有降水量、蒸散發、冰川補給等。為了深入了解薩雷茲堰塞湖水資源變化的驅動因素,文中結合氣象數據、估測的薩雷茲湖流域冰川變化以及薩雷茲湖體積變化,并進行了驅動力分析。

3.3.1 氣象因素

對薩雷茲湖水量變化的驅動力進行氣象因素方面研究,分析降水量的日數據變化發現,最大的降水量出現在春季(3—5月),但從7—9月幾乎沒有降雨。比較薩雷茲湖水量時序數據與降水量數據,發現兩者的月變化趨勢沒有較大關聯性,通過時間序列的匹配與篩選,將二者時間序列對齊進行相關系數的計算后結果為-0.265,這與常規的湖泊水量與總降水量呈正相關不符。

將薩雷茲湖月平均水量與氣溫變化比較,見圖3,發現兩者并無直接關系。采用全球陸面月平均地面數據集的潛在蒸散量對水量變化進行了進一步的相關分析,如圖4所示,1972—2019年間薩雷茲湖所在區域潛在蒸散量總體上呈相對平穩的趨勢,但水量與潛在蒸散相關系數為-0.103,整體上兩者之間表現出負相關性,并且相關性較弱。

圖3 薩雷茲湖月平均水量變化與月平均氣溫變化比較 圖4 薩雷茲湖月平均水量變化與潛在蒸散量變化比較Fig. 3 Comparison of the temperature data near Sarez Lake and the water volume change Fig. 4 Comparison of the potential evapotranspiration data near Sarez Lake and the water volume change

雖然月平均氣溫與薩雷茲湖月平均水量變化并無直接關系,但是薩雷茲湖水量的年際變化趨勢與年平均氣溫、年均最大氣溫(年內每月最大氣溫的平均值)、 年均最小氣溫(年內每月最小氣溫的平均值)的變化相比較,發現其趨勢一致(如圖5所示)。多年來,薩雷茲湖氣溫年際變化呈季節性波動, 并緩慢上升, 對齊二者時間序列相關分析結果分別為0.579、0.601、0.439,表明薩雷茲湖年際水量變化與區域年均溫度具有正相關性。推測其隨著氣溫的升高影響薩雷茲湖的水文特征變化狀況,使薩雷茲湖水資源得到了補給。

圖5 薩雷茲湖年平均水量變化與多年年均氣溫、多年年均最高氣溫、多年年均最低氣溫變化比較Fig. 5 Annual temperature data, the maximum temperature data, and the minimum temperature data near Sarez Lake and comparison with annual water volume change

3.3.2 入湖補給

薩雷茲湖水量近年來始終保持周期性上漲,為更好地顯示水資源變化的影響因素,本研究給出了薩雷茲湖水位變化與穆爾加布河入湖流量的變化關系,研究發現水位變化與湖泊入湖流量變化的趨勢在大部分時間內表現一致,兩者的相關系數為0.538,相關性良好,因此可以認為薩雷茲湖的水位上漲很大程度上源于穆爾加布河對其的補給。由于薩雷茲湖的水位變化呈現明顯的季節性周期變化,這表明穆爾加布河的補給作用可能受到其他因素的影響。研究表明,在1972—2000年間,由于氣候變暖導致帕米爾高原東部的冰川融溶體量增大[27],流域內的冰川凍土消融,融水一部分可能通過穆爾加布河補給匯入薩雷茲湖,導致穆爾加布河入流量的上升,最終給薩雷茲湖的水位變化特征造成影響。

3.3.3 流域冰川消融

分析流域冰川的變化對薩雷茲湖水量變化的影響,對入湖流量與流域冰川變化進行了相關分析。由于穆爾加布河上的產流區主要包括融雪產流與冰川產流,其中融雪產流主要集中在3—7月,冰川產流集中在6—9月,如圖6(a)為入湖流量的月平均變化,可以看出薩雷茲湖入湖流量較大的月份在7—9月,屬于穆爾加布河的冰川產流階段之內。此外,2003—2009年間,薩雷茲湖入湖流量呈季節性減少趨勢;在2010—2018年間,入湖流量呈增加趨勢。圖6(b)所示,顯示了薩雷茲湖流域冰川面積與年平均入湖流量變化的關系,可以看出冰川面積變化與薩雷茲湖流量變化趨勢呈負的相關性,年平均入湖流量與北部冰川面積、南部冰川面積的相關系數分別為-0.531、-0.462。圖6(c)中冰川體積變化的趨勢與流域內的河流流量的變化趨勢變現出顯著負相關,入湖流量與流域總體冰川體積變化,北部山脈冰川體積變化,南部山脈冰川體積變化之間的相關系數分別為-0.748、-0.781、-0.612,可以看出,無論是南北部山脈冰川變化還是整體流域山脈的冰川與流域內的水資源補給都有著密切的聯系。圖6(d)顯示了流域冰川體積變化對薩雷茲湖水量變化造成的影響,薩雷茲湖在2003—2009年間年平均水量變化呈減少趨勢,在2010—2018年呈緩慢上升趨勢,與流域冰川體積在2003—2009年增加,2010—2019年間減少的趨勢相反。

圖6 薩雷茲湖流域冰川變化與入湖流量變化比較Fig. 6 Comparison of the inflow discharges of Sarez Lake and glacier changes

可以看出,薩雷茲湖的水資源變化與流域冰川積雪消融密切相關,2003—2009年間受多種氣候因素影響流域冰川的消融過程減弱,此時,積雪消融的徑流為主要入湖補給,入湖流量與薩雷茲湖壩體滲透量相比可能相對較小,因此,總體上2003—2009年間薩雷茲湖水量呈下降趨勢;2010—2018年間流域冰川受到同期溫度的影響,冰川消融呈緩慢加劇趨勢,冰川融水導致穆爾加布河補給量增加,最終影響薩雷茲湖水位上升。假設在冰川系統和傳輸損耗(蒸發、滲透損失)過程中的蒸發和升華過程引起的質量損失可忽略不計,則根據入湖流量年平均變化與冰川體積年平均變化,計算2010—2018年間入湖流量在損耗的冰川體積中的年均占比可知:2010—2018年間薩雷茲湖流域冰川質量損失中有超過15.4%的比例被匯入穆爾加布河的河道產流,最終匯入薩雷茲湖,導致湖泊水量增多。

4 結論與討論

文中針對薩雷茲湖流域受地形因素和天氣要素的影響,傳統方法難以高效觀測的問題,采用光學影像、氣象數據等對薩雷茲湖水資源變化進行監測,并探討湖水變化的相關性因素,得到的主要結論如下:

1)結合長時間序列Landsat影像與薩雷茲湖實測水位數據,監測了薩雷茲湖1972—2019年間的湖岸廓線并計算了水域面積時間序列,結果表明薩雷茲湖水域面積變化在2010—2013年間存在突變拐點,變化趨勢加劇,薩雷茲湖呈現出一定的季節性周期變化規律,并且四季水域面積變化各有所區別,1972—2019年水域面積整體上呈現了逐漸增加的趨勢。

2)結合水文站點實測水位數據建模,同時進行了湖泊體積變化的計算,構建了薩雷茲湖1972—2018年的水量變化時間序列。表明1972—2018年薩雷茲湖水量整體以0.012 km3/a的變化速率波動上升。依據變化趨勢不同,分為了6個不同時間段的結果水量并呈現出“緩慢上升—緩慢下降—快速上升—緩慢上升—緩慢下降—緩慢上升”的變化規律。

3)結合氣溫、降水、潛在蒸散等氣象因素以及薩雷茲湖入湖流量、流域冰川變化,分析了薩雷茲湖水資源變化的驅動因素。分析結果表明薩雷茲湖水資源變化與降水補給、蒸散發在統計意義上均無顯著的相關關系,但與年平均溫度變化呈顯著正相關,還與穆爾加布河徑流補給量密切相關。綜合薩雷茲湖流域冰川的變化,表明,薩雷茲湖流域冰川的消融與其水文變化特征具有直接聯系,根據入湖流量年均變化與冰川體積年均變化可知2010—2019年間,薩雷茲湖流域冰川質量損失的15.4%多在穆爾加布河的河道產流中。

下一步研究需要考慮短時間內的降雪與融溶易對實驗造成影響以及對驅動因素研究需要增加大量驗證數據等問題,進而為全球變暖的大背景下薩雷茲堰塞湖潰壩的風險評估提供科學支撐。