薩雷茲湖東向引水最低成本路徑及可行性分析

李 君，張 新，陳 燕，胡榮強(qiáng)，朱長明，王偉勝

(1.東華大學(xué) 環(huán)境科學(xué)與工程學(xué)院，上海 201620； 2.中國科學(xué)院空天信息創(chuàng)新研究院 遙感科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，北京 100101；3.天津飛眼無人機(jī)科技有限公司，天津 300459； 4.江蘇師范大學(xué) 地理測(cè)繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇 徐州 221116；5.中國科學(xué)院新疆生態(tài)與地理研究所 荒漠與綠洲生態(tài)國家重點(diǎn)實(shí)驗(yàn)室，新疆 烏魯木齊 830011；6.中國科學(xué)院中亞生態(tài)與環(huán)境研究中心，新疆 烏魯木齊 830011)

位于帕米爾高原中部的全球最大堰塞湖——薩雷茲湖，以其豐富的水資源而廣為人知，但近年來，由于湖水水位持續(xù)抬升、右岸滑坡等不穩(wěn)定因素卻時(shí)刻威脅著下游人民的安全，這些不安全因素持續(xù)為國內(nèi)外學(xué)者所擔(dān)憂。Schuster和Alford[1]考慮到薩雷茲湖蓄水以來，烏索伊(Usoi)大壩的穩(wěn)定性一直備受爭議；Risley等人[2]提出Usoi大壩一旦潰決，下游人民時(shí)時(shí)面臨洪水沖擊。該湖及其壩體的問題切實(shí)威脅著中亞五國的生態(tài)安全，雖也多次進(jìn)行中亞地區(qū)間合作，但遲遲懸而未決?，F(xiàn)有關(guān)于薩雷茲湖周圍地區(qū)研究，更多專注在防治其災(zāi)害、維持其穩(wěn)定性上，集中于周圍地質(zhì)斷層構(gòu)造[3]，附近地震活動(dòng)[4-5]以及西向Usoi大壩所面臨的滑坡隱患[6-8]。20世紀(jì)70年代中期有學(xué)者曾試想利用薩雷茲湖進(jìn)行農(nóng)業(yè)灌溉[9]后證實(shí)不可行,也有用薩雷茲湖和其他湖泊進(jìn)行水力發(fā)電的討論[10]。經(jīng)多方資料總結(jié),主要降低該地區(qū)風(fēng)險(xiǎn)性的措施有加固壩體[11]、修建排水口以降低湖水位[12]并利用水能資源[11]。薩雷茲湖作為世界上最大的堰塞湖，其位于高海拔地區(qū)，屬高風(fēng)險(xiǎn)、高危險(xiǎn)區(qū)，降低其水位的同時(shí)，可考慮在評(píng)估堰塞體穩(wěn)定[13-14]的基礎(chǔ)上同時(shí)進(jìn)行合理的水資源開發(fā)。綜合多方啟發(fā)，又發(fā)現(xiàn)薩雷茲湖距中國邊界直線距離約為80km，跨區(qū)域調(diào)配水資源至中國，也有一定的地理優(yōu)勢(shì)?？死锬镜萚15]提及新疆雖已建立一些跨流域調(diào)水工程和水利大壩，但水資源調(diào)配仍存在問題，缺水問題依舊嚴(yán)重。由此，可考慮將薩雷茲湖作為水源地，將其引致新疆缺水地區(qū)，一方面，降低壩體潰決的危險(xiǎn)性，另一方面可試圖解決新疆缺水地區(qū)的缺口問題。在維持烏索伊堰塞體穩(wěn)定的同時(shí)，利用遙感技術(shù)進(jìn)行模擬、預(yù)測(cè)并試圖建立相應(yīng)水利工程以增進(jìn)薩雷茲湖水力資源可利用性，對(duì)高原堰塞湖具有前瞻性、探索性、互惠性的科學(xué)意義。

引水工程選線為GIS中求解最低成本路徑的問題，Dijkstra算法為經(jīng)典算法之一[16]。GIS輔助選線、選址[17]將復(fù)雜工程問題進(jìn)行簡化處理，以技術(shù)化后得的結(jié)果進(jìn)行人工勘測(cè)、修正。引水工程涉及研究區(qū)水文、地質(zhì)、地形、地震帶[18-19]等多方面因素，利用DEM進(jìn)行水文分析可估測(cè)跨區(qū)域調(diào)水水文情況[20]。從地形分析入手，坡度、地勢(shì)起伏度、坡向則選為主要地形影響因子。其中地勢(shì)起伏度最佳統(tǒng)計(jì)單元的提取以均值變點(diǎn)法[21]較為主流。由此，引入水文、地形分析作為選線理論分析基礎(chǔ)。

本文以ASTER GDEM V3為數(shù)據(jù)源，參考相關(guān)最低成本路徑模型[22-23]，在此基礎(chǔ)上構(gòu)建薩雷茲湖引疆路線規(guī)劃模型，利用水文分析預(yù)估跨區(qū)域調(diào)水范圍，結(jié)合地形分析下的坡向、坡度、地勢(shì)起伏度作為主要空間分析指標(biāo)，借助Arcpy以均值變點(diǎn)法自動(dòng)提取地勢(shì)起伏度最佳單元，借助專家打分以及層次分析法分配各專題地圖權(quán)重，科學(xué)選址確定受水區(qū)，利用Dijkstra算法自動(dòng)生成基于最小成本路徑下的2條備選線路，借助Google Earth軟件輔助并解讀，從工程量、安全性、生態(tài)性角度比較兩方案優(yōu)劣及可能性，在維護(hù)高海拔地區(qū)堰塞湖生態(tài)安全性的同時(shí)，也可幫助缺水地區(qū)解決水源問題，不乏作為一舉雙得的新方向。

1 研究區(qū)及數(shù)據(jù)來源

1.1 研究區(qū)概況

1911年塔吉克斯坦東部的帕米爾高原上發(fā)生了強(qiáng)烈地震，烏索伊堰塞體阻塞了穆爾加布河而形成迄今全球最大的堰塞湖，即薩雷茲湖(圖1)。其地理坐標(biāo)范圍：72°35′7″E～73°13′25″E，38°9′2″N～38°18′43″N。其湖長55.8 km，湖面平均海拔3 263 m，蓄水量約17 km3，水位最大周期性波動(dòng)達(dá)12 m[24]。薩雷茲湖壩體很少有排泄口，水位波動(dòng)性較大，再加上薩雷茲湖位于高山區(qū)，科研人員僅在夏季才能對(duì)其進(jìn)行無人機(jī)勘測(cè)，這給薩雷茲湖的災(zāi)害防治增加了很多的工程難度和技術(shù)挑戰(zhàn)。

圖1 薩雷茲湖Fig.1 Sarez Lake

1.2 數(shù)據(jù)來源

文中的數(shù)字高程模型(DEM)數(shù)據(jù)選自2019年8月5日由NASA′s LP DAAC新發(fā)布的ASTER GDEM V3版本，分辨率為30m，下載地址為：https://e4ftl01.cr.usgs.gov/ASTT/ASTGTM.003/2000.03.01。相較于V2版本，新增了36萬光學(xué)立體像對(duì)數(shù)據(jù)，數(shù)據(jù)較新，另也減少高程值空白區(qū)域、水域數(shù)值異常。

2 研究方法及模型構(gòu)建

2.1 引水路線規(guī)劃模型

文中從距離、水資源合理開發(fā)等角度，提出薩雷茲湖東向引疆規(guī)劃的思路。路線總模型構(gòu)建如下(圖2)：(1)預(yù)選研究范圍，鑲嵌拼接，投影變換處理；(2)預(yù)處理后進(jìn)行水文分析提取流域范圍，即集水區(qū)范圍，掩膜提取研究樣區(qū)DEM圖，目的是為了框定合理研究范圍，作為空間研究范圍的依據(jù)；(3)坡度、地勢(shì)起伏度、坡向作為選線因子，主要是考慮到坡度、地勢(shì)起伏度越大則經(jīng)過該區(qū)域的成本就越高，反之則越?。黄矫?、陽坡較陰坡、半陽半陰坡受太陽輻射較多，積雪較少，對(duì)于引水施工、材料運(yùn)輸具有一定地形優(yōu)勢(shì)；(4)引入ArcGIS空間分析模塊，地形分析生成坡度、坡向圖，結(jié)合均值變點(diǎn)法思想并自編Arcpy自動(dòng)化提取地勢(shì)起伏度，確定最佳統(tǒng)計(jì)單元得地勢(shì)起伏圖，較精準(zhǔn)地呈現(xiàn)研究區(qū)復(fù)雜地貌特征，進(jìn)一步優(yōu)化選線模型；(5)邀請(qǐng)專家對(duì)地形因子打分，以層次分析法固定權(quán)重分配，得各專題分級(jí)圖；(6)將各專題分級(jí)圖進(jìn)行空間疊置分析，選取目的點(diǎn)和不同終點(diǎn)，借以Dijkstra算法思想[25]，進(jìn)行距離分析，得成本回溯鏈接、成本距離后生成備選最低成本引水路線圖；(7)結(jié)合備選路線及工程量等方面，簡析各方案優(yōu)缺點(diǎn)及可行性。

圖2 引水工程路線規(guī)劃模型Fig.2 The planning model of water diversion project route

2.2 薩雷茲湖引疆路徑

由于薩雷茲湖東向引水工程涉及范圍較廣，橫跨兩國，地形、地勢(shì)較為復(fù)雜，需進(jìn)行數(shù)據(jù)預(yù)處理。先預(yù)選取71°E～78°E，35°N～41°N范圍，對(duì)柵格圖層進(jìn)行鑲嵌、拼接；采用Albers投影，確定中央經(jīng)線95°E，第一條標(biāo)準(zhǔn)緯線15°N，第二條標(biāo)準(zhǔn)緯線65°N，基準(zhǔn)面選擇D_WGS_1984；在預(yù)選區(qū)域進(jìn)行填洼、流量、流向計(jì)算等水文分析步驟，得集水區(qū)域；掩膜提取薩雷茲湖集水區(qū)區(qū)域，進(jìn)一步確定研究樣區(qū)的范圍，以減免不必要的運(yùn)算步驟為后續(xù)數(shù)據(jù)處理提高精準(zhǔn)度并增快計(jì)算速率。圖3為經(jīng)數(shù)據(jù)預(yù)處理后的研究樣區(qū)DEM圖。

圖3 研究樣區(qū)DEM圖Fig.3 DEM of study sample area

經(jīng)上述預(yù)處理后，主要運(yùn)算步驟如下：引入Arcpy模塊編寫代碼得最佳窗口大小并統(tǒng)計(jì)研究樣區(qū)地勢(shì)起伏度導(dǎo)入Excel中，另利用均值變點(diǎn)法原理在Excel中對(duì)其進(jìn)行驗(yàn)證；利用ArcGIS空間分析模塊，得坡向、坡度、地勢(shì)起伏度圖，根據(jù)專家對(duì)地形要素三因子進(jìn)行兩兩打分，以層次分析法[26]思想建立如下權(quán)重計(jì)算：

(1)構(gòu)造判斷矩陣A，求最大特征根λmax及相對(duì)應(yīng)的特征向量W。

AW=λmaxW.

(1)

(2)運(yùn)用方根法，對(duì)特征向量進(jìn)行歸一化處理，得權(quán)重系數(shù)。

(2)

(3)

(3)驗(yàn)證其相容性，進(jìn)行一致性檢驗(yàn)。

經(jīng)上述計(jì)算并通過一致性檢驗(yàn)，確定合理權(quán)重值，并利用1～9標(biāo)度法打分如表1。

表1 地形要素打分賦值Table 1 Terrain factor score assignment

對(duì)坡度、地勢(shì)起伏度、坡向各圖層綜合賦值加權(quán)。利用重分類工具進(jìn)行打分值再分配，運(yùn)用柵格計(jì)算器對(duì)重分類后的新圖層進(jìn)行加權(quán)計(jì)算(式4)，疊加分析后生成路徑成本。路徑成本計(jì)算公式如下：

(4)

其中，Cost是路徑成本，g1，g2，g3為表1中坡度、地勢(shì)起伏度、坡向的分級(jí)打分，w1，w2，w3為表1中坡度、地勢(shì)起伏度、坡向的權(quán)重值。

最后，確定引水點(diǎn)和受水點(diǎn)，基于路徑成本下，融合Dijkstra算法思想，派生成本回溯鏈接圖，成本距離圖，最終生成最低成本引水路線。其中，本研究涉及的Dijkstra算法原理如圖4。

圖4 Dijkstra算法示例Fig.4 Example of Dijkstra algorithm

指定起始引水點(diǎn)A和最終受水點(diǎn)E，引入集合Q和集合U，分別包含已求出和未求出的最短路徑點(diǎn)及相應(yīng)最短路徑長度。初始化兩個(gè)集合，Q初始時(shí)只有起始節(jié)點(diǎn)，即A—>A=0；U初始時(shí)為A—>B=4，A—>C= ∞，A—>D=2，A—>E=∞。從U中找出路徑最短的點(diǎn)(A—>D=2)加入Q，然后判斷D到B、C、E的距離加上A到D的距離之和是否小于A到B、C、E的距離，是則更新U。迭代執(zhí)行該更新步驟，直至遍歷得到起始A到終止節(jié)點(diǎn)E的最短路徑。

3 自動(dòng)化提取地勢(shì)起伏度

地勢(shì)起伏度能較好的反映地勢(shì)起伏變化，由柵格鄰域最高高程值減去最低高程值得到[27]。它的提取一般采用窗口分析法，按矩形區(qū)域提取(窗口設(shè)置如表2)，統(tǒng)計(jì)在3×3，4×4，…,44×44窗口下的地勢(shì)起伏度。結(jié)合空間分析中的鄰域分析方法，并利用均值變點(diǎn)法[28]尋找拐點(diǎn)，綜合上述方法編寫的Arcpy模塊添加進(jìn)ArcGIS中，自動(dòng)提取矩形網(wǎng)格下的地勢(shì)起伏度并導(dǎo)入表格，運(yùn)算時(shí)長較短，一定程度上縮短了運(yùn)算時(shí)間。

表2 地勢(shì)起伏度分析窗口設(shè)置Table 2 Analysis window setting of terrain relief

利用均值變點(diǎn)法的原理進(jìn)行以下分步運(yùn)算：

1)計(jì)算單位地勢(shì)起伏度。

Ti=Qi/Pi.

(5)

式中：Ti為第i個(gè)分析窗口下單位起伏度，m；Qi為平均起伏度，m；Pi分析窗口面積，m2。

2)數(shù)列{Ti}取對(duì)數(shù)運(yùn)算，得數(shù)列{xi}；

3)數(shù)列{xi}，令i=1,2,3…,N,(N取44),將其分成兩A，B兩數(shù)列，A：{xi,x2,X3,…,xi-1}，B：{xi,xi+1,…,xN}統(tǒng)計(jì)量公式如下：

(6)

(7)

利用自編的Arcpy模塊自動(dòng)提取導(dǎo)入表格，整理后得圖5，對(duì)數(shù)據(jù)進(jìn)行擬合分析，舍去擬合度不高的模型，其中乘冪和對(duì)數(shù)函數(shù)為兩較好的擬合曲線，冪函數(shù)得y= 35.082x0.471，R2= 0.9861，對(duì)數(shù)函數(shù)得y= 63.728ln(x) - 6.0419，R2= 0.9324。結(jié)果表明，冪函數(shù)趨勢(shì)線對(duì)平均起伏度擬合度最高，對(duì)數(shù)函數(shù)次之。經(jīng)統(tǒng)計(jì)驗(yàn)證，可利用均值變點(diǎn)法尋找突變點(diǎn)。由表3發(fā)現(xiàn)，i從起點(diǎn)開始至16，S-Sk數(shù)據(jù)持續(xù)上升，當(dāng)i為16以后S-Sk數(shù)值明顯回落，由此可得，該點(diǎn)對(duì)應(yīng)的窗口面積(23.04×104m2)即為最佳分析窗口，對(duì)應(yīng)像元大小16×16為最佳統(tǒng)計(jì)單元。

圖5 平均起伏度與分析窗口面積Fig.5 Average fluctuation degree and analysis window

表3 S-Sk數(shù)值統(tǒng)計(jì)Table 3 S-Sk statistic

4 備選方案分析

預(yù)設(shè)引水項(xiàng)目的起始點(diǎn)為薩雷茲湖東部支流出水口，方案A目標(biāo)點(diǎn)選擇英吉沙國家濕地公園(地理坐標(biāo)：76°4′15″E～76°14′53″E，38°50′12″N～38°56′23″N)，屬噶爾河水系，具有較為優(yōu)質(zhì)的淡水資源，地理位置靠近喀什地區(qū)，方案B目標(biāo)點(diǎn)取布倫口水庫(地理坐標(biāo)：74°52′31″E～75°1′16″E，38°39′9″N～38°46′34″N)，現(xiàn)庫容4.4億m3，依托現(xiàn)有布倫口-公格爾梯級(jí)水電站改造升級(jí)，能有效降低一定經(jīng)濟(jì)成本。由此確定引水點(diǎn)和受水點(diǎn)，得兩條基于地形分析下最低成本路徑圖(見圖6)。

圖6 備選引水路線圖Fig.6 Alternative water diversion routes

初設(shè)規(guī)劃路線：路線A東起薩雷茲湖出水口，沿巴爾坦格河引入穆爾加布河，利用天然河道長約4.73 km后，朝東北向經(jīng)恰塔蓋山而后越過我國國境線到達(dá)我國新疆邊界阿克陶縣，朝北向引流臨近已建成公路，利用其天然水道途徑木吉鄉(xiāng)，而后沿拜斯闊特別勒山谷向東線延升約70 km，穿過國道G314公路，經(jīng)塔什克米鄉(xiāng)、龍甫鄉(xiāng)，最終到達(dá)英吉沙縣的英吉沙國家濕地公園，總線路長約334.81 km。路線B先從薩雷茲湖東向出水口，沿巴爾坦格河匯入普沙爾季加爾比河、紹克塞河，經(jīng)穆爾戈布區(qū)、切切克季地區(qū)抵達(dá)我國國界線，經(jīng)舒爾庫爾湖、蘭克庫爾湖，穿過依結(jié)布陶、科克迪克地爾塔格山朝東北向引入布依口水庫，總路線約191.78 km。

4.1 工程量及安全性

結(jié)合圖7所示，線路A、B主要規(guī)劃為，在薩雷茲湖東向出水口的下游新建一處水利大壩，以便調(diào)節(jié)蓄水量，遇高海拔地區(qū)(4 000～5 000m)考慮開挖引水隧洞。A、B兩線路中多處遇明顯高差地勢(shì)，可在其合理規(guī)劃下建造梯級(jí)型水電站，利用水位高差發(fā)電，提升經(jīng)濟(jì)效益，此外沿線多處利用自然湖泊河流，多為冰山積雪的水源，水源質(zhì)量有一定保障，這部分可作為引水明渠，在工程上也可采用避免因引流挖渠而增加工程量和技術(shù)成本。引水隧洞需考慮邊坡構(gòu)造、地質(zhì)、海拔等多種要素，A線路多段長距離引水則需開挖隧洞約100 km以上，B方案引水隧洞距離相對(duì)A較少，目前已有相關(guān)長距離輸水[29]論證，我國現(xiàn)有水利技術(shù)上可初步實(shí)現(xiàn)高海拔長距離引水管道鋪設(shè)。

圖7 規(guī)劃方案線路剖面圖Fig.7 Plan of line profile

結(jié)合有關(guān)引水技術(shù)資料及沿線圖發(fā)現(xiàn)，前期A線路走勢(shì)均為4 500 m以上的高寒地帶，盡管依托現(xiàn)有技術(shù)條件可以嘗試攻克，對(duì)于建設(shè)團(tuán)隊(duì)有一定的安全性考驗(yàn)，施工周期較長，A線路中后段可利用自由落差進(jìn)行梯級(jí)水電站的升級(jí)，最終可引流至英吉沙濕地公園，需另修蓄水水庫。B線路前50 km需建造引水爬坡的動(dòng)力渠道，需較大動(dòng)力能耗，此外還另挖引水隧洞數(shù)十千米以上，但相較于A海拔高度挑戰(zhàn)略低。

A、B線路各有其優(yōu)劣，但薩雷茲湖東向引疆的難點(diǎn)主要集中于高海拔地區(qū)引水隧洞的挖建、引流爬坡所費(fèi)能耗以及梯級(jí)型水電站的建造。

4.2 生態(tài)性

干旱區(qū)生態(tài)環(huán)境敏感、脆弱，薩雷茲湖引水進(jìn)疆要特別注重可能產(chǎn)生的生境破壞，尋求生態(tài)環(huán)境與開發(fā)利用水資源之間的平衡。經(jīng)Google Earth輔助觀測(cè)，備選沿線并無主要建筑樓宇，對(duì)于生態(tài)破壞主要集中于施工期挖建隧洞、排水污染等。備選線路受水區(qū)分別為濕地公園及水庫，提供了棲息水源，這樣既維持受水區(qū)的生物多樣性，生物量也會(huì)隨之增加。沿線的引水明渠還能涵養(yǎng)水土，改善土壤質(zhì)量，緩解土地荒漠化，調(diào)節(jié)部分環(huán)境的小氣候，這對(duì)沿線環(huán)境起到了生態(tài)修復(fù)的作用。

5 結(jié)束語

薩雷茲湖作為塔吉克斯坦國的一大堰塞湖，其潛在的風(fēng)險(xiǎn)性不言而喻，由于該國的政治局勢(shì)及其他中亞四國的經(jīng)濟(jì)壓力，導(dǎo)致薩雷茲湖的問題一直未能得到有效改善。本研究從地勢(shì)、距離方面著手，利用水文及空間分析重點(diǎn)結(jié)合地形要素，構(gòu)建適用的引水路線模型，自動(dòng)派生兩條最低成本下的備選路線，綜合實(shí)際引水路線所需工程量、安全性以及生態(tài)性角度簡要初探了薩雷茲湖東向引入我國境內(nèi)新疆自治區(qū)的兩條方案的優(yōu)缺點(diǎn)以及應(yīng)用。雖然高海拔地區(qū)施工有一定操作難度，但基于我國現(xiàn)有的水利水電工程水平，應(yīng)該是可以攻克的難題。通過上述研究為全球最大堰塞湖薩雷茲湖災(zāi)害轉(zhuǎn)移提供初步的規(guī)劃方案，同時(shí)也能增進(jìn)一帶一路中亞國家間互惠合作，對(duì)后續(xù)解決、防治、轉(zhuǎn)移高原堰塞湖的災(zāi)害問題有一定的啟發(fā)和探索意義。