基于地形起伏度的冰湖潰決隱患研究
——以希夏邦馬峰東部為例

賀鵬， 童立強(qiáng)， 郭兆成， 涂杰楠， 王根厚

(1.中國地質(zhì)大學(xué)(北京)地球科學(xué)與資源學(xué)院，北京 100083；2.中國自然資源航空物探遙感中心，北京 100083)

0 引言

在全球氣候變暖背景下，近年來青藏高原暖濕化特征顯著，冰川活動性增強(qiáng)，冰湖潰決事件頻發(fā)[1]。西藏地區(qū)是我國冰湖潰決災(zāi)害高發(fā)區(qū)，以喜馬拉雅山區(qū)尤甚[2]。自20世紀(jì)30年代以來，該區(qū)有文獻(xiàn)記載的潰決冰湖33處，冰湖潰決事件37次，對我國西藏南亞大通道和川藏鐵路等重要交通干線建設(shè)構(gòu)成嚴(yán)重威脅[3-6]。冰湖潰決的形成機(jī)制復(fù)雜，誘發(fā)成因多樣，其中冰崩、冰滑坡及相關(guān)組合誘因占70%[7]。冰崩、冰滑坡等災(zāi)害多發(fā)于斜坡處，由冰川斷裂崩解形成[8]。目前，冰湖潰決危險(xiǎn)性判別及評價流程已形成較為完整的體系[9]，其中斜坡地質(zhì)災(zāi)害敏感性是定量評價過程中的重要指標(biāo)之一[10]，而地形起伏度又是實(shí)現(xiàn)準(zhǔn)確評價的前提[11]。

作為描述區(qū)域地形特征的宏觀指標(biāo)，地形起伏度是指某一確定面積中最高點(diǎn)與最低點(diǎn)海拔的差值[12]，可以真實(shí)反映地表起伏特征變化，與斜坡地質(zhì)災(zāi)害，特別是冰崩災(zāi)害的發(fā)生顯著相關(guān)[13]。近年來，基于地貌發(fā)育的基本理論，前人在不同地區(qū)和尺度下展開研究： 涂漢明等[12]利用國家數(shù)字地形數(shù)據(jù)庫，基于最大高差法得出中國地形起伏度最佳統(tǒng)計(jì)單元為21 km2； 劉愛利等[14]基于中國1∶100萬數(shù)字高程模型(digital elevation model,DEM)，探討DEM定量化、自動化的地貌分類應(yīng)用，得出在1∶100萬比例尺下最佳統(tǒng)計(jì)單元為6.5 km2； 韓海輝等[15]應(yīng)用均值變點(diǎn)法提取青藏高原地區(qū)地形起伏特征并分級分析，認(rèn)為最佳統(tǒng)計(jì)窗口為1.17 km2； 張軍等[16]以新疆1∶25萬DEM數(shù)據(jù)為基礎(chǔ)，利用鄰域分析及均值變點(diǎn)法得出2.56 km2的窗口尺度為最佳統(tǒng)計(jì)單元； 畢曉玲等[17]基于人工判讀法得出四川省地形起伏度的最佳統(tǒng)計(jì)單元為2.25 km2，并進(jìn)一步探討了其在區(qū)域地質(zhì)災(zāi)害評價中的應(yīng)用； 千琳勃等[18]利用均值變點(diǎn)法獲取了六盤山區(qū)地形起伏度，確認(rèn)最佳統(tǒng)計(jì)單元為1.96 km2，并分析了研究區(qū)地形空間變化特征； 張競等[19]針對京津冀地區(qū)不同地貌樣本計(jì)算了最佳分析窗口，認(rèn)為地貌起伏在400 m內(nèi)以4.64 km2為佳，而起伏度大于400 m的以5.35 km2更為適合； 王巖等[20]運(yùn)用均值變點(diǎn)法，在青海德貴地區(qū)開展1∶5萬比例尺下地形起伏度研究，論證了最佳分析窗口為0.16 km2。以往研究中以中低尺度下全國及西部地區(qū)的居多，以均值變點(diǎn)法的應(yīng)用頻率最高、效果最好[21]。然而，在我國藏南冰湖集中區(qū)1∶5萬比例尺下的研究較少，將地形起伏度應(yīng)用于冰湖潰決隱患識別方面的分析案例還較少。在當(dāng)前冰湖潰決災(zāi)害頻發(fā)的嚴(yán)峻形勢下，亟須開展研究，以豐富并完善我國冰湖潰決隱患分析評價方法。本文以希夏邦馬峰東部冰湖集中區(qū)為研究區(qū)，以國產(chǎn)高分二號(GF-2)衛(wèi)星影像及ASTER GDEM V3為數(shù)據(jù)源，利用均值變點(diǎn)法確定地形起伏度最佳提取單元，結(jié)合冰湖分布特征開展?jié)Q隱患相關(guān)性分析。

1 研究區(qū)概況

研究區(qū)位于聶拉木縣以北，地處喜馬拉雅山中段希夏邦馬峰東南部區(qū)域，在行政區(qū)劃上隸屬西藏自治區(qū)日喀則市聶拉木縣，范圍涵蓋波絨鄉(xiāng)、門布鄉(xiāng)、亞來鄉(xiāng)及聶拉木鎮(zhèn)等4個鄉(xiāng)鎮(zhèn)，總面積約730 km2(圖1)。地理坐標(biāo)介于E28.17°～28.5°，N85.75°～86.00°之間。

圖1 研究區(qū)遙感影像

受青藏高原新構(gòu)造運(yùn)動影響，形成整體以高落差基巖山地夾持縱深河谷的地貌特征，區(qū)內(nèi)冰湖分布廣泛、重力地貌發(fā)育、各類地質(zhì)災(zāi)害頻發(fā)[22]。受近年暖濕化氣候影響，喜馬拉雅山區(qū)冰川退縮明顯，冰湖顯著擴(kuò)張，潰決風(fēng)險(xiǎn)加劇，潰后引發(fā)鏈?zhǔn)酱紊鸀?zāi)害事件也隨之增多[23]。據(jù)前人研究，希夏邦馬峰東部及周邊地區(qū)存在多個高潰決風(fēng)險(xiǎn)冰湖，其中嘉龍錯、次仁瑪錯和帕曲錯等已發(fā)生潰決，對下游居民的生命財(cái)產(chǎn)安全造成巨大損失[24-25]。

2 研究數(shù)據(jù)與方法

2.1 冰湖信息提取

采用人機(jī)交互解譯方法，基于國產(chǎn)GF-2衛(wèi)星正射級產(chǎn)品數(shù)據(jù)，并配合圖像增強(qiáng)提取處理等方法，實(shí)現(xiàn)對區(qū)內(nèi)冰湖信息的準(zhǔn)確獲取。數(shù)據(jù)由中國自然資源航空物探遙感中心提供。

2.2 地形起伏度計(jì)算

2.2.1 地形數(shù)據(jù)

本文采用ASTER GDEM V3版地形數(shù)據(jù)。該數(shù)據(jù)于2019年8月發(fā)布，覆蓋了從N83°～S83°的所有陸地表面，其全球空間分辨率為30 m，垂直分辨率約為10 m，較V2版主要減少了高程值空白區(qū)域及水域數(shù)值異常[26]。數(shù)據(jù)來源于美國國家航空航天局(https: //earthdata.nasa.gov/)。

2.2.2 均值變點(diǎn)法

地形起伏度一般基于DEM數(shù)據(jù)通過窗口分析法求取，其關(guān)鍵是確定最佳統(tǒng)計(jì)窗口大小。以某一起伏地形上任意點(diǎn)為中心計(jì)算地形起伏度，隨著統(tǒng)計(jì)窗口尺度的擴(kuò)大，相對高程差與地形起伏度的變化率呈對數(shù)曲線特征，其中必定存在一個變化率由大變小(即曲線由陡轉(zhuǎn)緩)的“拐點(diǎn)”，該點(diǎn)有且唯一，對應(yīng)的統(tǒng)計(jì)窗口大小即適宜計(jì)算尺度。本文通過均值變點(diǎn)法確定研究區(qū)最佳統(tǒng)計(jì)單元，具有客觀性，相關(guān)計(jì)算原理參見相關(guān)文獻(xiàn)[15-20]，在此不再贅述。

3 適宜計(jì)算尺度確定

3.1 主要步驟

以n×n的矩形作為分析窗口，對希夏邦馬峰東部冰湖集中區(qū)進(jìn)行窗口尺度遞增的起伏度計(jì)算，規(guī)定起止矩形窗口尺度分別為3×3個像元和65×65個像元，單次增幅為2個像元(表1)。

表1 地形起伏度分析窗口尺度設(shè)置

基于ArcGIS軟件的Spatial Analyst模塊下neighborhood statistics工具對區(qū)內(nèi)DEM柵格數(shù)據(jù)進(jìn)行不同尺度窗口下的遍歷計(jì)算，得到各窗口的柵格最大及最小值； 利用Raster Calculator工具計(jì)算起伏度值，并將分析窗口大小與該尺度下起伏度值進(jìn)行關(guān)聯(lián)(表1)； 運(yùn)用Band Collection Statistics工具進(jìn)行柵格數(shù)據(jù)統(tǒng)計(jì)分析，并將結(jié)果屬性表輸出為表格； 利用VBA編程采用均值變點(diǎn)法計(jì)算統(tǒng)計(jì)量，其計(jì)算公式為：

(1)

(2)

將表1中各項(xiàng)數(shù)值代入式(1)—(2)中，實(shí)現(xiàn)對窗口面積與起伏度的曲線擬合與均值變點(diǎn)法中S-Sj的計(jì)算。變點(diǎn)的存在會使S和Sj的差距增大，S-Sj的最大值對應(yīng)的分析窗口大小即為最佳分析窗口。

3.2 地形起伏度擬合曲線

前人研究表明地勢起伏度與分析窗口大小的變化曲線為邏輯斯蒂曲線[12]，因此，對地形起伏度和統(tǒng)計(jì)單元面積(分析窗口)進(jìn)行曲線擬合，如圖2所示。

圖2 窗口大小與地形起伏度擬合曲線

在圖2的擬合曲線中，決定系數(shù)R2為0.950 5，擬合度良好，通過統(tǒng)計(jì)學(xué)檢驗(yàn)。地形起伏度隨統(tǒng)計(jì)單元尺度的擴(kuò)大而增加，變化曲線表現(xiàn)為對數(shù)函數(shù)或冪函數(shù)特征。統(tǒng)計(jì)單元面積小于32.49×104m2(窗口尺度19×19)時，其增長趨勢迅速； 統(tǒng)計(jì)單元面積在32.49×104～47.61×104m2區(qū)間(窗口尺度在19×19～23×23)時，地形起伏度的增速由急轉(zhuǎn)緩； 而統(tǒng)計(jì)單元面積超過47.61×104m2(窗口尺度23×23)之后，地形起伏度的增速放緩并趨于平穩(wěn)。

3.3 最佳窗口尺度確定

如圖2所示，擬合曲線內(nèi)存在一個增速逐漸趨緩的點(diǎn)(非數(shù)學(xué)拐點(diǎn))，即所謂“變點(diǎn)”，該點(diǎn)對應(yīng)的統(tǒng)計(jì)面積即為最佳分析窗口尺度。將起伏度計(jì)算結(jié)果逐一代入統(tǒng)計(jì)量公式，得到非線性系統(tǒng)數(shù)列樣本j的S與S-Sj值，其中S-Sj最大值所對應(yīng)的分析窗口大小，即所求的適宜計(jì)算尺度，計(jì)算結(jié)果見表2。

表2 均值變點(diǎn)法分析統(tǒng)計(jì)結(jié)果

根據(jù)表2數(shù)據(jù)，構(gòu)建S-Sj與非線性系統(tǒng)數(shù)列樣本j的變化擬合曲線(圖3)，可以看到該曲線為近似拋物線型，在j=10時出現(xiàn)S-Sj的最大值，可知此處為所求的變點(diǎn)，即最佳分析窗口的尺度，其窗口尺度為21×21，統(tǒng)計(jì)單元面積為39.69×104m2。

圖3 S-Sj值的變化趨勢

4 地形起伏度對冰湖潰決的影響分析

4.1 基于遙感的冰湖提取結(jié)果

基于GF-2衛(wèi)星數(shù)據(jù)的冰湖信息提取結(jié)果顯示(圖4)，區(qū)內(nèi)現(xiàn)存冰湖多達(dá)1 020個，總面積17.67 km2。從規(guī)模上來看，可分為大型(≥1 km2)、中型(0.1～1 km2)及小型(≤0.1 km2)3類； 就類型而言，包括冰蝕湖、冰磧阻塞湖、冰斗湖及冰川槽谷湖4類。區(qū)內(nèi)冰蝕湖僅3處，分別是位于希夏邦馬峰東部地區(qū)的崗西錯、貢錯及嘎龍錯，其總面積達(dá)12.21 km2，是全區(qū)僅有的3個大型冰湖； 冰磧阻塞湖18處，主要是分布在希夏邦馬峰東南部地區(qū)的嘉龍錯和達(dá)熱措等冰湖，其中中型冰湖8處、小型冰湖10處，總面積為3.12 km2； 冰斗湖30處，均為小型冰湖，以中南部地區(qū)分布居多，總面積0.41 km2； 冰川槽谷湖969處，規(guī)模上除2處中型冰湖外其余均為小型，多分布于冰川侵蝕作用形成的谷地內(nèi)。

圖4 冰湖提取結(jié)果及分類

4.2 地形起伏度分級

依據(jù)上述統(tǒng)計(jì)分析結(jié)果，采用21×21的網(wǎng)格大小(39.69×104m2)對研究區(qū)內(nèi)地形起伏度進(jìn)行計(jì)算。為便于后續(xù)分析，進(jìn)一步對區(qū)內(nèi)起伏度進(jìn)行分級，具體分為6個等級： 起伏度在[0,30) m為平原，[30,70) m為臺地，[70,200) m為丘陵，[200,500) m為小起伏山地,[500,1 000) m為中起伏山地； [1 000,2 500) m為大起伏山地。

結(jié)果顯示(圖5)，區(qū)內(nèi)平原分布相對較少，整體以小起伏山地為主，在拉布日峰至俄熱村之間、希夏邦馬峰東南側(cè)多有分布； 丘陵及中起伏山地所占比例相近，前者多見于北部納朵熱至大來村、崗西錯東部等地，后者主要分布于中部及南部冰湖集中分布區(qū)域； 大起伏山地相對較少，主要位于中西部的希夏邦馬峰及西南部的多爾雷山附近。

圖5 研究區(qū)地形起伏度分級

4.3 地形起伏度對冰湖潰決的影響分析

如前文所述，地形起伏度是定量化分析斜坡地質(zhì)災(zāi)害敏感性的重要因子之一，更是冰湖潰決危險(xiǎn)性評價體系中不可或缺的指標(biāo)，對研究全區(qū)冰湖潰決隱患具有重要的指示意義。本文以冰湖周邊地形起伏組合特征為依據(jù)，疊加冰湖空間分布特征，以分析區(qū)內(nèi)地形起伏度對冰湖潰決的影響。

4.3.1 地形起伏度與冰湖潰決的關(guān)系

通過地形起伏度與冰湖空間分布特征的疊加分析發(fā)現(xiàn)就冰湖后緣及兩側(cè)起伏度特征而言，主要可分為單一起伏等級型、多起伏等級組合型及跨起伏等級型3類冰湖(以下簡稱單一型、組合型和跨級型)。

單一型冰湖周緣起伏度類型一致，多發(fā)育于冰川槽谷區(qū)，湖泊類型上以小型冰磧湖為主。區(qū)內(nèi)地形起伏度較小、地勢變化緩、穩(wěn)定性高，不具備潰決條件。從影像特征及實(shí)地驗(yàn)證結(jié)果看，沖堆普流域(希夏邦馬峰與多爾雷山之間)的典型冰湖(圖6(a))特征與分析結(jié)果一致，作為冰川退縮消融過程中的殘留湖泊，多處于同一地形起伏等級(圖6(b))，周邊地勢平緩(圖6(c))，冰湖間多由常年地表徑流相連，由于湖面較淺、流動性好，整體潰決隱患極低。

(a) 單一型冰湖影像 (b) 單一型冰湖起伏特征 (c) 單一型冰湖實(shí)地調(diào)查照片

(d) 組合型冰湖影像 (e) 組合型冰湖起伏特征 (f) 組合型冰湖實(shí)地調(diào)查照片

(g) 跨級型冰湖影像 (h) 跨級型冰湖起伏特征 (i) 跨級型冰湖實(shí)地調(diào)查照片

圖6 典型冰湖與地形起伏度疊加分析圖及實(shí)地照片

組合型冰湖周邊起伏度類型具顯著差異，湖泊后緣多與冰川相接，湖泊類型以中-小型冰斗湖和冰磧阻塞湖為主。此類湖泊所處地質(zhì)環(huán)境條件相對復(fù)雜，周邊各地勢起伏等級單元均有分布，后緣冰川易崩落入湖，形成涌浪可能誘發(fā)潰決，潰決概率較高，以嘎龍錯較為典型。遙感影像及實(shí)地驗(yàn)證結(jié)果顯示，嘎龍錯周邊地形特征(圖6(d))與分析結(jié)果相符，屬組合型(圖6(e))，冰湖前緣終磧堤兩側(cè)及后緣部分區(qū)域地形高陡(圖6(f))，局部具備冰滑坡、冰崩等成災(zāi)條件，可能引發(fā)冰湖潰決災(zāi)害。

跨級型冰湖周邊起伏度的等級變化不連續(xù)，多為冰蝕湖、冰磧湖，規(guī)模均在中型以下，以多爾雷山東北部的嘉龍錯最為典型。從遙感影像(圖6(g))及實(shí)地調(diào)查結(jié)果看，該湖東側(cè)地勢較緩，而西北側(cè)起伏度變化明顯(圖6(h))，與分析結(jié)果一致。冰湖側(cè)磧壟與后緣相接處可見顯著地形起伏度跨級界線(圖6(i)紅線標(biāo)識位置)，湖面與后緣冰川具有一定高差，上部冰川易失穩(wěn)發(fā)生大規(guī)?；?，進(jìn)而誘發(fā)冰湖潰決災(zāi)害，因此這類冰湖發(fā)生潰決的可能性極高。

4.3.2 對冰湖潰決隱患的評價

基于上述分析，可進(jìn)一步依據(jù)地形起伏度對冰湖潰決隱患進(jìn)行分級評價(表3),其中單一型冰湖規(guī)模不一、以小型為主，由于周邊地形平緩，不具備成災(zāi)條件，故整體潰決隱患低； 組合型冰湖規(guī)模以小型居多，具一定成災(zāi)條件，潰決隱患中等； 跨級型冰湖規(guī)模適中，后緣冰川所處地形落差較大，具有顯著的成災(zāi)條件，潰決隱患高。區(qū)內(nèi)屬低潰決隱患的冰湖993處，占比97.35%，在規(guī)模上涵蓋小、中、大3種類型冰湖，較典型的有大型冰湖貢錯及中型冰湖達(dá)熱措等； 中潰決隱患的冰湖18處，占比1.77%，規(guī)模上包括小型及大型2類，以大型冰湖崗西錯、嘎龍錯等最為典型； 高潰決隱患的冰湖9處，占比0.88%，包括小型及中型2種規(guī)模的冰湖，其中屬中型冰湖嘉龍錯潰決隱患最高(圖7)。

表3 冰湖潰決隱患評價結(jié)果

圖7 地形起伏度對冰湖潰決隱患的影響分析

從評價結(jié)果中不同隱患等級冰湖的數(shù)量來看，全區(qū)近98%的冰湖潰決隱患極低，僅有不足2%的冰湖存在潰決的可能； 但就湖泊規(guī)模而言，中-高隱患的冰湖面積達(dá)11.37 km2，占到全區(qū)冰湖總面積的64.31%，其對下游河道兩側(cè)牧民以及聶拉木縣城構(gòu)成了極大的威脅。

值得注意的是，冰湖潰決造成的損失與其規(guī)模大小并非絕對正相關(guān)。例如2014年那隆藏布支溝發(fā)生的冰湖潰決事件[6]，該冰湖潰前面積僅0.24 km2，按本文湖泊規(guī)模劃分方案屬中型，由于冰湖后緣冰川與常年積雪區(qū)邊坡高陡，在連續(xù)降雨條件下冰川發(fā)生崩塌后沖入湖中，進(jìn)而誘發(fā)潰決； 在2016年7月聶拉木縣樟木鎮(zhèn)上游發(fā)生的次仁瑪錯右側(cè)支溝冰湖潰決事件中[6]，冰湖潰決前面積僅0.01 km2，屬小型規(guī)模，但潰后形成洪水導(dǎo)致樟木口岸被淹，造成較大經(jīng)濟(jì)損失。

因此，本次依據(jù)地形起伏度評價出的9處高潰決隱患冰湖雖規(guī)模偏小，但在當(dāng)前全球氣候暖濕化的背景下，由冰川退縮消融、冰崩/冰滑坡、埋藏冰融化、強(qiáng)降水等誘發(fā)潰決的可能性極高，潛在威脅較大。特別是位于多爾雷山東北處沖堆普流域內(nèi)的嘉龍錯曾發(fā)生過2次潰決事件[27]，再次發(fā)生災(zāi)害的概率較高，需重點(diǎn)關(guān)注。

5 結(jié)論

本文基于GIS技術(shù)對西藏聶拉木縣境內(nèi)希夏邦馬峰東南部冰湖集中分布區(qū)域開展地形起伏度研究，從冰湖地形起伏度組合特征的角度出發(fā)，對區(qū)內(nèi)冰湖潰決隱患進(jìn)行分析評價，得出如下結(jié)論：

1)利用均值變點(diǎn)法實(shí)現(xiàn)區(qū)內(nèi)地形起伏度的提取，經(jīng)有效計(jì)算及準(zhǔn)確性檢驗(yàn)，該區(qū)在1∶5萬比例尺下最佳采樣單元為21×21(采樣單元面積為0.39 km2)，提取結(jié)果可準(zhǔn)確反映地形變化特征。

2)研究區(qū)整體地勢以小起伏山地為主，中起伏山地及丘陵次之，大起伏山地、臺地及平原分布較少，多見高落差基巖山地夾持縱深河谷的地貌特征。

3)區(qū)內(nèi)冰湖周邊地形起伏度特征可分為單一起伏等級型、多起伏等級組合型、跨起伏等級型3種組合類型，對應(yīng)潰決隱患等級依次升高，以跨起伏等級型冰湖潰決隱患最大。

4)全區(qū)1 020個冰湖中，潰決隱患為低、中、高級別的數(shù)量分別占總體的97.35%，1.77%和0.88%，中-高隱患的冰湖面積占到全區(qū)總體的64.31%，高潰決隱患冰湖中以嘉龍錯特征最為典型，對下游河道兩側(cè)牧民以及聶拉木縣城構(gòu)成威脅。

5)結(jié)果顯示，以地形起伏度作為冰湖潰決隱患評價分析指標(biāo)，對完善現(xiàn)有冰湖潰決危險(xiǎn)性判別方法和體系具有實(shí)際意義，可有效提高評價結(jié)果精度，為同類地區(qū)開展相關(guān)研究提供理論參考。