羅永臻，董春

(中國測繪科學研究院，北京 100036)

0 引言

近50年當中青藏高原地區(qū)氣候發(fā)生顯著變化，融水增加，但由于青藏高原特殊的地理位置，該地區(qū)的水文資料等相關數(shù)據(jù)難以獲取，而數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)既保證了一定的地形精度，在數(shù)據(jù)量上又比其他遙感數(shù)據(jù)更為小巧且容易獲取，在盡可能小地耗費人力物力的情況下，可以得到較為滿意的研究結果。而在2011年前后，位于可可西里地區(qū)的卓乃湖發(fā)生潰決導致洪水流向下游庫塞湖區(qū)，受此影響，庫塞湖出現(xiàn)湖水外溢的情況，導致其下游的海丁諾爾出現(xiàn)類似情況，并最終匯入海丁諾爾東南側的鹽湖[1]。楊勇[2]認為，卓乃湖潰決的主要原因在于古河道季節(jié)河流的溯源侵蝕所導致，而姚曉軍等[3]認為可可西里地區(qū)降水增加和蒸發(fā)減少是該地區(qū)湖水潰決的主要原因。不管何種原因，卓乃湖、庫塞湖、海丁諾爾和鹽湖已經(jīng)建立出了自上而下影響的水文關系。

而下游的鹽湖作為該水文系統(tǒng)中的最終接收者，在2018年前后也出現(xiàn)了與上游相同的潰決溢流風險。不同于上游3個湖區(qū)的情況是，鹽湖東部有青藏鐵路和青藏公路高寒區(qū)國家工程，較大的湖水溢流可能會對這些工程的路基造成破壞[4]，所以需要進行鹽湖溢流的模擬分析，從而確定其大致影響范圍。本文就是要在快速反應迅速確定影響面積的條件下，基于數(shù)字高程模型來進行鹽湖區(qū)域的水文模擬分析。2002年，吳險峰等[5]利用DEM數(shù)據(jù)對黃河小花間流域進行了數(shù)字降水徑流模擬，得出了研究區(qū)域的流域模擬圖；2005年，黃金良等[6]利用DEM數(shù)據(jù)對九龍江流域進行水文模擬獲取到了該地區(qū)流量模擬；Tarekkegn等[7]2010年利用15 m的DEM數(shù)據(jù)作為水動力模型的輸入數(shù)據(jù)，精確地模擬了塔納湖水位變化影響；Gichamo等[8]于2012年利用全球數(shù)字高程模型進行垂直偏差校正后，對匈牙利Tisza河部分地區(qū)的洪水進行模擬，得到了較好結果，證明了DEM數(shù)據(jù)在該方面的應用潛力；在2016年，李林娜等[9]基于ArcGIS平臺與90 m DEM數(shù)據(jù)對元江縣內(nèi)的水域進行了河網(wǎng)提取，并通過水文模擬對居民地規(guī)劃提供了參考；2021年，劉鳳梅等[10]通過DEM數(shù)據(jù)水文分析對珠江三角洲流域進行了流域盆地范圍模擬，在當?shù)赝恋乩靡?guī)劃中發(fā)揮了重要的指導意義。

目前，基于DEM數(shù)據(jù)進行的水文分析都基本保持在河網(wǎng)提取和盆地分析層面，而本文在進行水文分析與模擬的基礎上又根據(jù)河網(wǎng)與研究區(qū)高程值進行了湖泊溢流影響面積模擬，便于確定鹽湖溢流后造成的影響面積并迅速對鹽湖溢流可能造成的影響范圍進行研究和探討。

1 研究區(qū)域與數(shù)據(jù)情況

1.1 研究區(qū)域概況

研究區(qū)域為卓乃湖、海丁諾爾的下游鹽湖地區(qū)，又名68道班鹽湖[11]，位于格爾木市西南方，屬于青海省玉樹藏族自治州治多縣，從海丁諾爾湖區(qū)至清水湖間海拔落差為4 462～4 497 m。湖區(qū)屬于青海南部高寒草原半干旱氣候，年平均氣溫為-4.0～-1.0 ℃，年降水量為150～200 mm，目前因為其位于可可西里自然保護區(qū)范圍內(nèi)，鹽業(yè)開采活動在20世紀80年代已經(jīng)停止。在鹽湖西側，自西向東分別為卓乃湖、庫塞湖和海丁諾爾湖。其中自2011年后，卓乃湖潰決后便迅速形成了與庫塞湖間的河道。

1.2 數(shù)據(jù)情況

受數(shù)據(jù)源限制，且為了實現(xiàn)快速預測湖泊溢流影響的目的，在保證一定精度的情況下，本研究用到的數(shù)字高程模型為航天飛機雷達地形測繪使命(shuttle radar topography mission，SRTM)數(shù)據(jù)，根據(jù)文獻[12-15]等的研究，選取分辨率為30 m的DEM初始數(shù)據(jù)。該數(shù)據(jù)由美國國家航空航天局(National Aeronautics and Space Administration，NASA)在2000年利用奮進號航天飛機上的雷達測觀測所得，是利用次數(shù)最多的高程數(shù)據(jù)，覆蓋了全球南北緯60°以內(nèi)的區(qū)域。由Wang等[16]研究可知，DEM數(shù)據(jù)可以計算水文資料并具有一定可靠性。2011年10月，海丁諾爾湖水還未進入到鹽湖當中，在這之后，海丁諾爾湖發(fā)生溢流現(xiàn)象，湖水溢流進入鹽湖[17]，如圖1所示。2012年開始，鹽湖庫容開始了快速增長的階段，從2010年10月至2012年10月間，鹽湖湖面積擴大了3倍左右，到2015年底，鹽湖面積增加到了將近150 km2，而在2015年之后，根據(jù)Landsat衛(wèi)星影像可以直接通過目視判別得出鹽湖的面積直到2020年9月份仍在持續(xù)不斷地增長(圖2)。由于其恰好與東方向的青藏鐵路與公路距離約12.5 km，所以對于鹽湖可能溢流后所造成的影響亟待進行提前預測與分析。

2 溢流淹沒分析方法

2.1 分析預測鹽湖溢流河道

1)DEM填洼處理。洼地是影響地表流水過程的重要因素。在自然條件下，水流從高處向低處流動，遇到洼地的時候肯定會先將其填滿后再從某一方向的最低出口處流出，但是洼地是局部的最低點，會導致無法確定該處的水流方向。因此，消除小洼地對確定水流的方向有重要的影響。在本研究中，研究區(qū)域為高寒高海拔凍土區(qū)域，且沒有較大的地表坑洞，一般都為小洼地，可以直接進行處理，對其高程進行賦值，將其填平(圖3)。

2)水流方向矩陣計算。根據(jù)物理模型，水流方向的矩陣計算有以下幾種方法：單流向算法(single flow direction algorithm，SFD)、多流向算法(multi-flow direction algorithms，MFD)以及其他算法。由于多流向算法認為水流分布具有分散性，即水流方向具有不確定性，可能同時流向相鄰8個鄰域中的幾個網(wǎng)格鄰域，這樣會造成計算量很大，算法復雜、工作效率低下。因此，本實驗采用單流向算法，根據(jù)DEM柵格單元與相鄰8個單元之間的最大坡降來確定水流方向。如圖4所示，將目標柵格的相鄰8個領域格網(wǎng)按二次冪值法進行順時針方向編碼，如值為2，就為東南向流向。

3)累計矩陣提取水流網(wǎng)絡。即用數(shù)學的方法來描述經(jīng)過流水的多少，通過流向格網(wǎng)進行流量累加然后得到所有格網(wǎng)數(shù)據(jù)的水流份數(shù)。而流域的匯流能力則是由確定流入每個柵格的累積上游柵格數(shù)目統(tǒng)計而成。將這種算法可以描述為：以每個柵格單元格為初始柵格，按順序?qū)α飨蚓仃囘M行遍歷計算，根據(jù)水流方向追蹤，直到DEM邊界，由此通過計算累積柵格數(shù)目即可確定水流方向和匯流量。需要注意的是，在實際水文預報研究模擬中，都要考慮地面的土壤滲透系數(shù)以及地表植被的影響、植物根系會吸收流水的系數(shù)等，而這些不是本研究的重點，只需根據(jù)研究區(qū)域高程與流向來計算水流量大小，并最終提取水流的路徑。

2.2 淹沒范圍模擬

在提取到鹽湖溢流河道后確定其上游溢出點與下游和清水河的匯水點的高程值后，用此范圍高程值從30 m的DEM數(shù)據(jù)中進行屬性值提取分析。由于高程值為一個范圍值，且提取次數(shù)較多，所以采用Python語言對DEM影像進行遍歷分析提取。

3 鹽湖湖水外溢模擬

3.1 提取溢出點

本研究通過ArcGIS平臺所提供的水文分析模塊，首先對研究區(qū)域的DEM數(shù)據(jù)進行洼地填充以消除小洼地對下一步流向計算的誤差影響。在獲得了洼地填充后的DEM數(shù)據(jù)后，利用D8算法對該數(shù)據(jù)進行流向計算獲得研究區(qū)域的流向圖(圖5)。根據(jù)其地形圖以及流向圖共同分析可以看出，鹽湖水的流向基本為正東至正南方向之間，在這之后，根據(jù)流向圖進行匯流累積量計算，并通過條件分析篩選出流入值大于2 000的像元值，由此獲得研究區(qū)域的河網(wǎng)圖。為了捕捉溢流點，利用Landsat 8全彩色影像疊加研究區(qū)域河網(wǎng)圖，根據(jù)流向圖在東南方向上目視捕捉到該河網(wǎng)匯流點即鹽湖溢流位置(圖6)，并根據(jù)相同方法提取到下游清水湖的入水點，各點地理位置及高程如表1所示，該地點位置與文獻[12]分析得到的入水點位置無較大差異，證明了該方法的有效性。

表1 出水點與入水點位置

3.2 淹沒模擬分析

當鹽湖湖面高程值超過其出水點高程時，則滿足其湖水外溢的條件。根據(jù)圖7位置示意圖發(fā)現(xiàn)溢流河道穿過青藏鐵路匯入清水湖，而鹽湖溢流點距離青藏鐵路直線距離約為12.5 km，鹽湖的擴容會極大威脅到青藏鐵路的安全。

當鹽湖發(fā)生溢流現(xiàn)象后，受地理因素影響會匯聚成新的河道并向外發(fā)展形成漫流，而本文根據(jù)所提取的溢流點和入流點高程值(表1)作為閾值，對原30 m DEM數(shù)據(jù)進行屬性值提取。通過屬性值提取，獲得不同高程值下的淹沒面積。由于研究區(qū)域水文數(shù)據(jù)難以獲取造成該研究無法利用水文模型進行時間序列上的分析，故本文將所模擬得到的淹沒影響范圍根據(jù)溢流湖面高程值疊加至2021年3月的Landsat 8真彩色影像上，并以鹽湖和清水湖為觀察對象劃定觀察范圍，從而進行被淹沒影響范圍的演示(圖8)，并得到鹽湖溢流后的湖區(qū)面積及所有受影響區(qū)域(圖9)。本文進行水文分析發(fā)現(xiàn)，如果鹽湖溢流造成下游清水湖也發(fā)生潰堤現(xiàn)象，最終將會導致溢流河道匯入楚瑪爾河形成長江支流。

4 結束語

本文基于D8算法，利用SRTM數(shù)據(jù)對研究區(qū)域進行水文分析與模擬，研究了鹽湖溢流后可能的淹沒區(qū)域，主要結論如下。

1)根據(jù)Landsat 8影像可知，可可西里鹽湖從2015—2020年間發(fā)生了明顯變化，由于上游卓乃湖潰堤導致下游庫塞湖、海丁諾爾發(fā)生了一系列的連鎖反應致使鹽湖出現(xiàn)了湖面積增加、湖區(qū)擴容。2015年，鹽湖湖面積為152.29 km2，至2020年12月，湖區(qū)面積已達到200.85 km2，在這5年間共增長了48.56 km2。

2)當鹽湖湖水發(fā)生溢流時，整個湖區(qū)湖泊面積將會達到219.94 km2。已知青藏鐵路路基單向?qū)挾葹?7 m，當湖泊發(fā)生溢流后，先從下游開始形成積水區(qū)，當下游清水湖湖面積突破5.84 km2時，青藏鐵路部分路基將會開始與水體接觸，在路基周圍形成小而淺的積水區(qū)從而影響路基的穩(wěn)定性。

3)若鹽湖發(fā)生溢流后不進行人為引流任由其自然發(fā)展，其溢流河道最終會匯入到下游楚瑪爾河中，成為一條長江支流。

本文結果證明，在缺失水文資料等數(shù)據(jù)的情況下，憑借遙感影像和DEM數(shù)據(jù)進行高原地區(qū)的湖泊溢流模擬分析是具有可行性的，由此可見DEM數(shù)據(jù)的廣泛應用性。但如果要對一些更加復雜的地形區(qū)域進行溢流模擬和精確估算，利用高精度地形和必要的水文數(shù)據(jù)會得到更加精確的結果。