羅永臻,董春
(中國測繪科學研究院,北京 100036)
近50年當中青藏高原地區氣候發生顯著變化,融水增加,但由于青藏高原特殊的地理位置,該地區的水文資料等相關數據難以獲取,而數字高程模型(digital elevation model,DEM)數據既保證了一定的地形精度,在數據量上又比其他遙感數據更為小巧且容易獲取,在盡可能小地耗費人力物力的情況下,可以得到較為滿意的研究結果。而在2011年前后,位于可可西里地區的卓乃湖發生潰決導致洪水流向下游庫塞湖區,受此影響,庫塞湖出現湖水外溢的情況,導致其下游的海丁諾爾出現類似情況,并最終匯入海丁諾爾東南側的鹽湖[1]。楊勇[2]認為,卓乃湖潰決的主要原因在于古河道季節河流的溯源侵蝕所導致,而姚曉軍等[3]認為可可西里地區降水增加和蒸發減少是該地區湖水潰決的主要原因。不管何種原因,卓乃湖、庫塞湖、海丁諾爾和鹽湖已經建立出了自上而下影響的水文關系。
而下游的鹽湖作為該水文系統中的最終接收者,在2018年前后也出現了與上游相同的潰決溢流風險。不同于上游3個湖區的情況是,鹽湖東部有青藏鐵路和青藏公路高寒區國家工程,較大的湖水溢流可能會對這些工程的路基造成破壞[4],所以需要進行鹽湖溢流的模擬分析,從而確定其大致影響范圍。本文就是要在快速反應迅速確定影響面積的條件下,基于數字高程模型來進行鹽湖區域的水文模擬分析。2002年,吳險峰等[5]利用DEM數據對黃河小花間流域進行了數字降水徑流模擬,得出了研究區域的流域模擬圖;2005年,黃金良等[6]利用DEM數據對九龍江流域進行水文模擬獲取到了該地區流量模擬;Tarekkegn等[7]2010年利用15 m的DEM數據作為水動力模型的輸入數據,精確地模擬了塔納湖水位變化影響;Gichamo等[8]于2012年利用全球數字高程模型進行垂直偏差校正后,對匈牙利Tisza河部分地區的洪水進行模擬,得到了較好結果,證明了DEM數據在該方面的應用潛力;在2016年,李林娜等[9]基于ArcGIS平臺與90 m DEM數據對元江縣內的水域進行了河網提取,并通過水文模擬對居民地規劃提供了參考;2021年,劉鳳梅等[10]通過DEM數據水文分析對珠江三角洲流域進行了流域盆地范圍模擬,在當地土地利用規劃中發揮了重要的指導意義。
目前,基于DEM數據進行的水文分析都基本保持在河網提取和盆地分析層面,而本文在進行水文分析與模擬的基礎上又根據河網與研究區高程值進行了湖泊溢流影響面積模擬,便于確定鹽湖溢流后造成的影響面積并迅速對鹽湖溢流可能造成的影響范圍進行研究和探討。
研究區域為卓乃湖、海丁諾爾的下游鹽湖地區,又名68道班鹽湖[11],位于格爾木市西南方,屬于青海省玉樹藏族自治州治多縣,從海丁諾爾湖區至清水湖間海拔落差為4 462~4 497 m。湖區屬于青海南部高寒草原半干旱氣候,年平均氣溫為-4.0~-1.0 ℃,年降水量為150~200 mm,目前因為其位于可可西里自然保護區范圍內,鹽業開采活動在20世紀80年代已經停止。在鹽湖西側,自西向東分別為卓乃湖、庫塞湖和海丁諾爾湖。其中自2011年后,卓乃湖潰決后便迅速形成了與庫塞湖間的河道。
受數據源限制,且為了實現快速預測湖泊溢流影響的目的,在保證一定精度的情況下,本研究用到的數字高程模型為航天飛機雷達地形測繪使命(shuttle radar topography mission,SRTM)數據,根據文獻[12-15]等的研究,選取分辨率為30 m的DEM初始數據。該數據由美國國家航空航天局(National Aeronautics and Space Administration,NASA)在2000年利用奮進號航天飛機上的雷達測觀測所得,是利用次數最多的高程數據,覆蓋了全球南北緯60°以內的區域。由Wang等[16]研究可知,DEM數據可以計算水文資料并具有一定可靠性。2011年10月,海丁諾爾湖水還未進入到鹽湖當中,在這之后,海丁諾爾湖發生溢流現象,湖水溢流進入鹽湖[17],如圖1所示。2012年開始,鹽湖庫容開始了快速增長的階段,從2010年10月至2012年10月間,鹽湖湖面積擴大了3倍左右,到2015年底,鹽湖面積增加到了將近150 km2,而在2015年之后,根據Landsat衛星影像可以直接通過目視判別得出鹽湖的面積直到2020年9月份仍在持續不斷地增長(圖2)。由于其恰好與東方向的青藏鐵路與公路距離約12.5 km,所以對于鹽湖可能溢流后所造成的影響亟待進行提前預測與分析。
1)DEM填洼處理。洼地是影響地表流水過程的重要因素。在自然條件下,水流從高處向低處流動,遇到洼地的時候肯定會先將其填滿后再從某一方向的最低出口處流出,但是洼地是局部的最低點,會導致無法確定該處的水流方向。因此,消除小洼地對確定水流的方向有重要的影響。在本研究中,研究區域為高寒高海拔凍土區域,且沒有較大的地表坑洞,一般都為小洼地,可以直接進行處理,對其高程進行賦值,將其填平(圖3)。
2)水流方向矩陣計算。根據物理模型,水流方向的矩陣計算有以下幾種方法:單流向算法(single flow direction algorithm,SFD)、多流向算法(multi-flow direction algorithms,MFD)以及其他算法。由于多流向算法認為水流分布具有分散性,即水流方向具有不確定性,可能同時流向相鄰8個鄰域中的幾個網格鄰域,這樣會造成計算量很大,算法復雜、工作效率低下。因此,本實驗采用單流向算法,根據DEM柵格單元與相鄰8個單元之間的最大坡降來確定水流方向。如圖4所示,將目標柵格的相鄰8個領域格網按二次冪值法進行順時針方向編碼,如值為2,就為東南向流向。
3)累計矩陣提取水流網絡。即用數學的方法來描述經過流水的多少,通過流向格網進行流量累加然后得到所有格網數據的水流份數。而流域的匯流能力則是由確定流入每個柵格的累積上游柵格數目統計而成。將這種算法可以描述為:以每個柵格單元格為初始柵格,按順序對流向矩陣進行遍歷計算,根據水流方向追蹤,直到DEM邊界,由此通過計算累積柵格數目即可確定水流方向和匯流量。需要注意的是,在實際水文預報研究模擬中,都要考慮地面的土壤滲透系數以及地表植被的影響、植物根系會吸收流水的系數等,而這些不是本研究的重點,只需根據研究區域高程與流向來計算水流量大小,并最終提取水流的路徑。
在提取到鹽湖溢流河道后確定其上游溢出點與下游和清水河的匯水點的高程值后,用此范圍高程值從30 m的DEM數據中進行屬性值提取分析。由于高程值為一個范圍值,且提取次數較多,所以采用Python語言對DEM影像進行遍歷分析提取。
本研究通過ArcGIS平臺所提供的水文分析模塊,首先對研究區域的DEM數據進行洼地填充以消除小洼地對下一步流向計算的誤差影響。在獲得了洼地填充后的DEM數據后,利用D8算法對該數據進行流向計算獲得研究區域的流向圖(圖5)。根據其地形圖以及流向圖共同分析可以看出,鹽湖水的流向基本為正東至正南方向之間,在這之后,根據流向圖進行匯流累積量計算,并通過條件分析篩選出流入值大于2 000的像元值,由此獲得研究區域的河網圖。為了捕捉溢流點,利用Landsat 8全彩色影像疊加研究區域河網圖,根據流向圖在東南方向上目視捕捉到該河網匯流點即鹽湖溢流位置(圖6),并根據相同方法提取到下游清水湖的入水點,各點地理位置及高程如表1所示,該地點位置與文獻[12]分析得到的入水點位置無較大差異,證明了該方法的有效性。

表1 出水點與入水點位置
當鹽湖湖面高程值超過其出水點高程時,則滿足其湖水外溢的條件。根據圖7位置示意圖發現溢流河道穿過青藏鐵路匯入清水湖,而鹽湖溢流點距離青藏鐵路直線距離約為12.5 km,鹽湖的擴容會極大威脅到青藏鐵路的安全。
當鹽湖發生溢流現象后,受地理因素影響會匯聚成新的河道并向外發展形成漫流,而本文根據所提取的溢流點和入流點高程值(表1)作為閾值,對原30 m DEM數據進行屬性值提取。通過屬性值提取,獲得不同高程值下的淹沒面積。由于研究區域水文數據難以獲取造成該研究無法利用水文模型進行時間序列上的分析,故本文將所模擬得到的淹沒影響范圍根據溢流湖面高程值疊加至2021年3月的Landsat 8真彩色影像上,并以鹽湖和清水湖為觀察對象劃定觀察范圍,從而進行被淹沒影響范圍的演示(圖8),并得到鹽湖溢流后的湖區面積及所有受影響區域(圖9)。本文進行水文分析發現,如果鹽湖溢流造成下游清水湖也發生潰堤現象,最終將會導致溢流河道匯入楚瑪爾河形成長江支流。
本文基于D8算法,利用SRTM數據對研究區域進行水文分析與模擬,研究了鹽湖溢流后可能的淹沒區域,主要結論如下。
1)根據Landsat 8影像可知,可可西里鹽湖從2015—2020年間發生了明顯變化,由于上游卓乃湖潰堤導致下游庫塞湖、海丁諾爾發生了一系列的連鎖反應致使鹽湖出現了湖面積增加、湖區擴容。2015年,鹽湖湖面積為152.29 km2,至2020年12月,湖區面積已達到200.85 km2,在這5年間共增長了48.56 km2。
2)當鹽湖湖水發生溢流時,整個湖區湖泊面積將會達到219.94 km2。已知青藏鐵路路基單向寬度為17 m,當湖泊發生溢流后,先從下游開始形成積水區,當下游清水湖湖面積突破5.84 km2時,青藏鐵路部分路基將會開始與水體接觸,在路基周圍形成小而淺的積水區從而影響路基的穩定性。
3)若鹽湖發生溢流后不進行人為引流任由其自然發展,其溢流河道最終會匯入到下游楚瑪爾河中,成為一條長江支流。
本文結果證明,在缺失水文資料等數據的情況下,憑借遙感影像和DEM數據進行高原地區的湖泊溢流模擬分析是具有可行性的,由此可見DEM數據的廣泛應用性。但如果要對一些更加復雜的地形區域進行溢流模擬和精確估算,利用高精度地形和必要的水文數據會得到更加精確的結果。