李志杰， 王寧練,3， 常佳雯

（1.陜西省地表系統(tǒng)與環(huán)境承載力重點實驗室，陜西 西安 710127；2.西北大學城市與環(huán)境學院地表系統(tǒng)與災害研究院，陜西 西安 710127；3.中國科學院青藏高原研究所，北京 100101）

喀喇昆侖公路（Karakoram Highway，KKH，以下簡稱公路）北起新疆喀什，南至巴基斯坦塔科特，是中巴兩國實現(xiàn)互聯(lián)互通，建設中巴經(jīng)濟走廊的戰(zhàn)略通道。公路穿越的東帕米爾高原、喀喇昆侖山是高亞洲冰川分布最集中的地區(qū)之一，冰川災害曾多次威脅和中斷公路通行，造成重大損失［1-3］。位于巴基斯坦北部罕薩河（洪扎河）中游的巴托拉段是整個公路沿線距離冰川最近、冰川災害最集中的路段，其中以巴托拉（Batura）、帕蘇（Pasu）、固爾金（Ghulkin）等3 條冰川對公路的威脅最大也最為直接［4-6］。

1973年，巴托拉冰川主排水道（巴托拉河）改道引發(fā)的冰川洪水沖毀了建設中的公路和公路橋。1974—1975 年，為提供合理的修復設計，中國考察組對巴托拉冰川開展了包括末端進退、流速、厚度、融水徑流等全方面、詳盡的科學考察［4］。最終考察組提出按照原線修復的科學建議，確保了公路的長期穩(wěn)定通行［7］。此后的1978 年、1980 年、1994 年等，考察組再次對巴托拉地區(qū)開展科學考察，積累了豐富的地質(zhì)地貌、冰川物理、水文、氣象等科學資料［7-8］。在此基礎上，學者們揭示了巴托拉、帕蘇、固爾金等冰川在20世紀近百年的進退變化，并預測巴托拉冰川將于2010s前后再次發(fā)生前進［4-5,7-8］。

進入21 世紀以來，隨著氣候波動的加劇，巴托拉、帕蘇、固爾金等冰川發(fā)生了許多新的變化，這些變化對公路的正常通行產(chǎn)生了新的威脅［9］。例如巴托拉冰川的主排水洞位置不斷北移，最終導致融水通道被冰舌擠壓切斷，時隔近50 a 后再次引發(fā)巴托拉河的大改道；帕蘇冰川末端的冰磧阻塞湖頻繁潰決，2008 年1 月和4 月2 次引發(fā)冰湖潰決洪水，中斷公路通行［10］；固爾金冰川末端的頻繁波動導致冰崖下的排水通道不斷遷移，僅在2008年就引發(fā)4次泥石流，致使公路橋底嚴重淤積［6,10-11］。因此，鑒于災害風險的日漸凸顯和公路無可取代的地位，有必要對巴托拉、帕蘇、固爾金冰川的歷史與新近變化進行揭示和分析，為公路的安全運行和災害防治提供參考。

1 研究區(qū)概況

巴托拉、帕蘇、固爾金冰川發(fā)源于喀喇昆侖山西段主脊——洪扎喀喇昆侖山北坡（圖1）。洪扎喀喇昆侖山平均海拔超過6000 m，最高峰巴托拉峰海拔高達7795 m，為冰川發(fā)育提供了極為有利的地形條件，使得巴托拉冰川群成為了公路沿線規(guī)模最大、末端海拔最低、進退最為頻繁的冰川群之一。巴托拉、帕蘇、固爾金冰川上游粒雪盆面積開闊，補給物質(zhì)充足，中下游則深入了干熱河谷，消融強烈，因此屬于典型的多溫型冰川，并表現(xiàn)出較高的物質(zhì)平衡水平［4］。本研究的解譯結(jié)果表明，2021 年巴托拉冰川長約56.5 km，面積307.12±4.32 km2，末端海拔約2620 m；帕蘇冰川長約26.4 km，面積64.81±0.72 km2，末端海拔約2700 m；固爾金冰川長約17.8 km，面積28.26±0.47 km2，末端海拔約2540 m。

圖1 巴托拉、帕蘇、固爾金冰川的地理位置Fig.1 Geographic location of the Batura,Pasu,and Ghulkin Glacier

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

2.1 數(shù)據(jù)源

2.1.1 Landsat 影像 本研究使用了獲取于1972—2022 年的25 景Landsat 系列影像，以揭示研究冰川的相關(guān)變化（表1）。所用影像均為L1T 級產(chǎn)品，由USGS（United States Geological Survey，https://earthex?plorer.usgs.gov/）進行了系統(tǒng)輻射校正和幾何校正，并結(jié)合DEM進行了地形校正。由于Landsat系列影像有相當高的正射校正精度，達到1/2 像元左右，因此本研究直接利用Landsat 影像開展冰川變化研究［12］。為提高冰川邊界、出水洞、融水徑流等地物的解譯或識別精度，本研究選取了消融季末期，冰川區(qū)無云、山體陰影覆蓋較小的影像［13］。

表1 遙感影像數(shù)據(jù)Tab.1 List of the Landsat images used in this study

2.1.2 冰面高程變化數(shù)據(jù)集 Brun等［14］利用ASTER立體像對提取DEM，基于大地測量學方法，通過DEM的空間匹配與誤差校正，重建了2000—2016年高亞洲冰川的表面高程變化與物質(zhì)平衡，并生成了高亞洲冰面高程變化數(shù)據(jù)集。數(shù)據(jù)集中的誤差主要源于DEM的生成精度及其空間匹配偏差，因此以非冰川區(qū)等地形穩(wěn)定區(qū)的高程殘差來校正并評價冰面高程變化的誤差［14］。該數(shù)據(jù)集的空間分辨率為30 m，可獲取于PANGAEA(https://doi.pangaea.de/10.1594/PANGAEA.876545)。

2.1.3 冰面流速數(shù)據(jù)集 NASA 的ITS_LIVE（Intermission time series of land ice velocity and elevation）產(chǎn)品包含了1985—2020 年高亞洲主要冰川的流速數(shù)據(jù)。該產(chǎn)品基于Landsat 影像，采用Gardner 等提出的自動裂縫特征跟蹤處理鏈方法，通過進行局部歸一化、過采樣、特征跟蹤等提取冰川流速［15-16］。影像的空間匹配偏差是該產(chǎn)品的主要誤差來源，盡管在生成過程中已基于基巖等穩(wěn)定地形對其進行校正，但無法徹底消除。因此，ITS_LIVE產(chǎn)品同樣基于非冰川區(qū)的流速殘差評價冰川流速的不確定性［15］。1999年之前（TM影像）冰川流速估算的不確定性超過10 m·a-1，1999年之后（ETM+影像）降低為約2 m·a-1，2013 年之后（OLI 影像）降低為約0.5 m·a-1［16-17］。本研究選用了其中1985—2018 年的冰川年平均流速合成數(shù)據(jù)，空間分辨率為120 m（https://its-live.jpl.nasa.gov/）。

2.2 方法

2.2.1 冰川邊界解譯 由于研究冰川表面存在大范圍表磧覆蓋，且進退變化頻繁，本研究采用目視解譯提取冰川邊界［13］。在中上游裸冰區(qū)，直接根據(jù)Landsat 假彩色影像（熱紅外、紅、綠波段組合）中顯著的色彩差異進行冰川邊界的數(shù)字化（圖1）［18］。在表磧覆蓋的冰川下游，則根據(jù)末端冰崖與出水洞的組合、高程變化等地形和水文特征識別冰川邊界［19］。假彩色影像中，裸冰區(qū)與表磧區(qū)、表磧區(qū)與基巖區(qū)、冰崖與水流等的色彩差異均得到了顯著增強，由此可準確識別冰舌、出水洞等的位置。野外考察資料表明，研究冰川冰舌消融強烈，因此末端冰崖與出水口的組合即代表了冰川運動的最前緣［4,6,20］。在實際解譯中，本研究首先完成2021年的冰川邊界提取。在此基礎上結(jié)合不同時期的影像，修改冰舌等發(fā)生變化的部位，以獲取其他年份的冰川邊界。2.2.2 解譯精度評價 基于光學影像解譯冰川邊界的精度主要受影像空間分辨率的影響，因此本研究采用統(tǒng)計冰川輪廓線經(jīng)過的像元數(shù)量的方法，評價目視解譯提取冰川邊界的不確定性［21］。公式如下：

E=N·λ/2

式中：E為冰川面積的不確定性；N為冰川邊界經(jīng)過的像元數(shù)量；λ為像元面積。

2.2.3 出水洞與融水通道識別 出水洞指位于冰川末端，冰內(nèi)和冰下水系的排水出口，本研究結(jié)合野外考察記錄和遙感影像識別冰川出水洞與融水通道的遷移變化。如圖2所示，在消融季末期的Land?sat 和Google Earth 影像中，得益于巴托拉冰川較大的排水量，冰川出水洞與融水通道的準確位置可獲得有效識別。此外，在1973—2021 年2 次大改道期間，巴托拉河下游保持穩(wěn)定，上游隨冰舌波動而延伸，因此結(jié)合不同時期的遙感影像即可準確識別出水洞與融水通道的變遷過程。

圖2 巴托拉冰川出水洞與融水徑流的識別Fig.2 Identification of drainage system of the Batura Glacier

3 結(jié)果與分析

3.1 冰川進退變化

近150 a 來，巴托拉冰川發(fā)生過多次進退波動。1873年巴托拉冰川發(fā)生前進，但此后近50 a間冰舌位置基本保持穩(wěn)定。1930s 冰川開始衰退，1944—1966年退縮827 m。1960s末轉(zhuǎn)入前進，1966—1980年前進約136 m。1980s初再次衰退，至1994年累計后退了150—200 m［4-5,7,22-23］。影像資料表明（圖3a），1994—2008 年巴托拉冰舌繼續(xù)后退約410 m，2008年冰舌再次前進，至2021 年7 月累計前進約230 m。綜上所述，近150 a 來巴托拉冰舌每隔50～60 a就會發(fā)生一輪前進，但每輪前進一般只可維持10～20 a，總體而言，近150 a 來巴托拉冰川以退縮為主導［23］。

圖3 1972—2021年巴托拉、帕蘇、固爾金冰川進退變化Fig.3 Terminus positions at different periods for the Batura,Pasu,and Ghulkin Glacier from 1972 to 2021

帕蘇冰川雖曾于1907—1913 年發(fā)生過大規(guī)模躍動，導致冰舌前進約1200 m［24］，但此后卻長期處于退縮狀態(tài)，僅在2000 年前后發(fā)生過一次前進。1972—1998年，帕蘇冰舌累計后退約180 m，1998—2002 年前進220 m，2002—2021 年再次后退約570 m。因此近百年來，帕蘇冰川總體處于加速退縮狀態(tài)（圖3b）。

近百年來，固爾金冰川進退頻繁，堪稱公路沿線最活躍的冰川之一。1885—1980年固爾金冰川發(fā)生了3 輪進退，期間還伴隨著多次小幅度波動［6］。1966—1978 年，固爾金南冰舌后退260 m，1978—1980年卻又前進260 m，此后前進速率減緩，1980—1994 年累計前進100 m［6］。之后再度衰退，至2021年累計退縮約390 m。北冰舌更加活躍（圖3c），1972—1990 年后退200 m，1990—1998 年又前進170 m。之后，冰舌開始在小范圍內(nèi)頻繁波動，最終在1998—2018 年累計前進120 m。2018 年起冰舌再次衰退，至2021年已后退約36 m。

3.2 融水通道變化

1973年6—7月，巴托拉河改道后冰川主出水洞由中流線附近轉(zhuǎn)移至冰舌南側(cè)，但此后主出水洞位置開始緩慢遷移。2000年之前遷移速率較慢，基本穩(wěn)定在冰舌南側(cè)。2000 年之后，遷移速率明顯加快，逐步越過中流線到達冰舌北側(cè)。至2021 年6月，主出水洞已處于冰舌西北側(cè)，排水道隨之延伸完全橫亙于冰舌前（圖4a）。2021年6—7月，冰舌的突然前進擠壓切斷了主排水道，引發(fā)了近50 a 來巴托拉河的又一次大改道（圖4b）。由于冰壩阻擋，融水徑流被迫在冰磧丘陵中侵蝕切割出一條新的通道，直接導致巴托拉河由2.9 km縮短為2.5 km，縱比降由65‰增大至76‰，因此流速顯著增快。

圖4 巴托拉、固爾金冰川融水徑流排水通道變化Fig.4 Variations of the meltwater runoff drainage channel of the Batura and Ghulkin Glacier

受地形影響，帕蘇冰川的排水通道十分穩(wěn)定。但在1981年之后，冰舌附近在融水補給下形成了一個小型終磧阻塞湖，2011年實測面積約0.12 km2［10］。在水量增加、冰舌崩塌、谷坡落石等的共同作用下，阻塞湖曾于2007年7月、2008年1月、4月等數(shù)次發(fā)生潰決，引發(fā)冰湖潰決洪水［3,10］。但經(jīng)過數(shù)次潰決以及融水徑流的長期侵蝕，近年來帕蘇冰湖終磧堤的完整性已遭到破壞，冰湖發(fā)育規(guī)模受到很大限制。

固爾金段是公路沿線受冰川融水影響最大的路段之一。由于冰舌分為南北兩支且融水徑流流程極短，不足1 km，因此冰川融水呈輻散狀，極不穩(wěn)定且頻繁發(fā)生改道。固爾金冰川主要排水通道約有6 條，小型通道則難以計數(shù)［6,10］。其中，北冰舌融水集中在M1通道，南冰舌的M2通道也已基本廢棄（圖4c）。長期以來，M5為南冰舌的主排水道，但近年來M4 的水量也在增大。M3 曾在2014 年前后短暫成為主排水道，但之后水量又迅速減少，M6 水量也較為有限。因此，近年來南冰舌融水通道以M5與M4為主，M3與M6為輔。

3.3 冰面高程與流速

近年來“喀喇昆侖異?！爆F(xiàn)象備受關(guān)注，但就處于喀喇昆侖山西北邊緣的巴托拉地區(qū)而言，2000年以來的冰川物質(zhì)虧損卻有所加速［20,25］?；贙H-9 DEM和SRTM DEM的研究發(fā)現(xiàn)，1973—1999年巴托拉和帕蘇冰川物質(zhì)平衡分別為0.00±0.10 m w.e.a-1和0.05±0.11 m w.e.a-1，表明2000年以前冰川物質(zhì)收支基本保持穩(wěn)定［20］。但此后的2000—2016年，巴托拉、帕蘇、固爾金冰川表面高程均發(fā)生了下降（圖5a），下降速率分別為0.09±0.29 m·a-1、0.15±0.17 m·a-1和0.08±0.26 m·a-1，尤以冰川末端表現(xiàn)最為明顯［14］。盡管下降速率較為微弱，但這依舊表明2000年以來巴托拉地區(qū)冰川發(fā)生了物質(zhì)虧損。

圖5 巴托拉、帕蘇、固爾金冰川表面高程變化（a）與年平均流速（b）Fig.5 Surface elevation change(a)and annual average flow velocity(b)of the Batura,Pasu,and Ghulkin Glacier

長期以來，3 條研究冰川均維持著較高的表面流速。1985—2018年，巴托拉、帕蘇、固爾金冰川的平均流速分別為59.60±0.28 m·a-1、106.44±0.44 m·a-1和19.50±0.35 m·a-1。就流速的空間分布而言（圖5b），高值區(qū)域均出現(xiàn)在冰川上游的粒雪線附近。巴托拉冰川的最大流速出現(xiàn)在粒雪盆出口，第一、第二冰流構(gòu)成的冰瀑布處，接近1200 m·a-1。帕蘇冰川的最大流速出現(xiàn)在相同的位置，約940 m·a-1。固爾金冰川流速相對平緩，最大流速也出現(xiàn)在了冰瀑布以下的附近位置。相較于巴托拉和固爾金冰川，帕蘇冰川流速較快，意味著其物質(zhì)平衡水平更高，這是帕蘇冰舌消融更強烈的主要原因。

4 討論

4.1 冰川變化的影響因素

山地冰川是氣候與地形共同作用的產(chǎn)物，但在十年或百年尺度上，冰川變化主要受氣候條件的影響［26］。冰川物質(zhì)平衡直接由當年氣候狀態(tài)決定，末端進退卻更多地受氣候長周期波動的影響，因為從氣候變化到冰舌波動需要一個能量積聚和動力反映的過程，即存在滯后性。滯后期的長短主要與冰川規(guī)模有關(guān)，小規(guī)模冰川為3～30 a，大規(guī)模冰川則會更長［4,27］。

就巴托拉、帕蘇、固爾金這種大型復式山地冰川而言，氣溫的長周期波動是影響冰舌進退的關(guān)鍵因素［8］。如圖5所示，2000年以來，冰川積累區(qū)基本保持穩(wěn)定，消融區(qū)高程則有所降低，表明近幾十年來罕薩河谷的升溫直接導致了冰川物質(zhì)虧損的加劇［28-29］。而巴托拉、帕蘇、固爾金冰舌在減薄的同時發(fā)生進退波動，則充分體現(xiàn)了冰舌進退對氣候響應的滯后性以及冰川規(guī)模對于滯后性的影響。有研究表明，近150 a來巴托拉冰舌進退相較于氣溫波動存在60 a 左右的滯后期，冰川的數(shù)次前進均與氣溫的階段性變化存在很好的對應關(guān)系［4,7-8］?；谶@一判斷，不僅成功預測了巴托拉冰舌在1970s—2000年的進退，還成功預測了冰舌在2010s 的最新一輪的前進。

躍動冰川指周期性發(fā)生快速運動的冰川，當躍動發(fā)生時其速度往往是平時的數(shù)倍至上百倍［9］。有研究將巴托拉第一、第二冰流等部位的高流速、大量裂隙與弧拱發(fā)育視為躍動的表現(xiàn)［24］，但實際上這對于大型山地冰川而言并不罕見。對巴托拉這種巨型冰川而言，支冰流躍動對末端進退的影響也十分有限，如克拉牙依拉克、費德琴科等冰川支冰流的躍動并未改變整體的退縮趨勢［30-31］。帕蘇冰川曾于1907—1913 年，固爾金冰川曾于1913—1925 年、1966—1978 年發(fā)生躍動［24］，但近50 a 來2 條冰川的積累區(qū)并未發(fā)生顯著增厚，因此在近期發(fā)生大規(guī)模躍動的可能性較低。

4.2 融水徑流的影響因素

較高的物質(zhì)平衡水平導致3條冰川的融水徑流量明顯高于其他地區(qū)的同等規(guī)模冰川，巴托拉冰川甚至高于近710 km2的費德琴科冰川［4］。巴托拉河多年平均流量為40.1 m3·s-1，瞬時最大流量可超過600 m3·s-1，其中近80%的徑流量集中在5—9 月，近60%集中在7—8 月［4,32］。巴托拉河90%的水量來源于冰川融水，因此徑流量主要受氣溫控制，隨氣溫上升而增加，在夏季降水條件下有時甚至會因氣溫下降而減少［4］。此外，即便在7—8月巴托拉河的徑流量波動也十分劇烈，常隨極端高溫天氣的出現(xiàn)而陡漲［32］。帕蘇、固爾金冰川融水徑流的年內(nèi)變化與巴托拉河類似，但變化過程遠不及前者劇烈。

巴托拉河徑流年際變化較為平緩，主要由于：（1）干旱年份降水偏少而融水偏多，濕潤年份降水偏多而融水偏少，起到了天然的調(diào)節(jié)作用［31］。（2）冰舌區(qū)年平均氣溫接近10 ℃，導致消融十分劇烈，存在熱融喀斯特現(xiàn)象，因此冰內(nèi)裂隙、孔穴，冰磧物的滯水能力也能起到一定的徑流調(diào)節(jié)作用［4,32］。（3）1970s 以來巴托拉流域冰川基本保持穩(wěn)定，即便2000年以來物質(zhì)虧損有所增加，但對融水徑流的影響依舊有限。實際上不僅是巴托拉河，1980s 以來整個罕薩河流域的年徑流量也無明顯變化［33］。

巴托拉河在1973 年和2021 年的2 次改道有一定的相似性，改道均發(fā)生在6—7月融水徑流快速上升的時期、改道前主排水洞均已發(fā)生長期緩慢的遷移、改道的發(fā)生均與冰舌突然前進有關(guān)，因此2次改道的原因也是基本相同的。如圖1 所示，巴托拉冰川的物質(zhì)來源極不對稱，冰舌北側(cè)主要由尤克蘇戈茲、北、西冰流等構(gòu)成、南側(cè)主要由第一、第二冰流構(gòu)成。由于南側(cè)冰流擁有更大的高差與更開闊的粒雪盆，因此物質(zhì)供給量更大，直接導致冰舌南側(cè)厚度與流速均高于北側(cè)［4,34］。冰川中下游，南側(cè)弧拱與裸冰的延伸范圍明顯大于北側(cè)就是二者存在流速差的直觀表現(xiàn)。在這種流速差的長期作用下，主排水洞受到推擠而緩慢北移，排水道隨之延長并橫亙于冰舌前，最終被突然前進的冰舌切斷而被迫改道。因此，冰舌南北兩側(cè)物質(zhì)供給量的差異是巴托拉河改道的根本原因，冰舌的突然前進則是直接原因［4］。

4.3 對公路的影響

本研究在公路沿線選取了4 個穩(wěn)定的參考點，分別為巴托拉橋、帕蘇橋、M1橋、M4橋，以參考點與冰舌的最近距離來衡量公路與冰川的距離。如圖6所示，近50 a來巴托拉冰舌與公路的距離總體上在增大，當前尚有1 km 以上，因此在短期內(nèi)不會對公路產(chǎn)生直接威脅。帕蘇冰舌自2000 年以來不斷加速后退，與公路的距離也在增大，因此短期內(nèi)也不會直接威脅公路。盡管固爾金南北冰舌僅在小范圍內(nèi)波動，但由于距離過近，且在將來不能排除發(fā)生躍動的可能，因此對公路的直接威脅將長期存在。

圖6 1972—2021年巴托拉、帕蘇、固爾金冰舌與公路的距離變化Fig.6 Variations of the distances from glacier tongues of the Batura,Pasu,and Ghulkin Glacier to the Karakoram highway from 1972 to 2021

巴托拉冰川對公路更現(xiàn)實的威脅在于巴托拉河新的大改道。2021年6—7月改道后，巴托拉河流程縮短，比降增大，下切侵蝕能力增強，近1 a年來，新河道明顯拓寬加深（圖7）。更關(guān)鍵的是新河道下游折向南流，與公路基本平行，因此水流將在轉(zhuǎn)彎處直接沖擊路基。甚至完全可以推測，正是由于路基堅固才迫使水流南流返回故道。巴托拉河具有流量大、徑流高度集中、陡漲陡落等特點［32］，因此在未來的強烈消融季中，隨著水量陡漲，超過600 m3·s-1的水流裹挾冰塊砂石奔騰下泄，將對路基產(chǎn)生強烈侵蝕，對公路構(gòu)成重大威脅。

圖7 2022年5月27日巴托拉冰舌區(qū)概況Fig.7 Overview of Batura Glacier tongue on May 27,2022

帕蘇冰川近年來消融強烈導致融水徑流增大，但由于冰湖發(fā)育規(guī)模受限，大規(guī)模冰湖潰決洪水對于下游公路的威脅趨于減弱。固爾金冰川融水輕則在路面積水結(jié)冰，影響通行，重則侵蝕路基、淤積橋底，甚至引發(fā)泥石流洪水。但在2008年采取工程措施后，固爾金段的融水災害已得到了有效治理［3］。然而，固爾金冰川融水徑流的遷移擺動問題仍不容忽視，一方面它可能會導致過水工程廢棄，但更重要的是徑流量可能會超過單一通道的排水能力，進而引發(fā)新的災害。

5 結(jié)論

中巴喀喇昆侖公路巴托拉段受冰川災害威脅最為嚴重，是影響公路暢通運行的關(guān)鍵路段。基于野外考察與遙感影像等資料，分析了巴托拉、帕蘇、固爾金冰川近百年來的變化、變化原因及對中巴公路的影響，主要結(jié)論如下：

（1）近百年來，特別是2000 年以來，巴托拉、帕蘇、固爾金冰川總體處于微弱退縮狀態(tài)，冰舌發(fā)生后退減薄或基本保持穩(wěn)定，對公路的直接威脅趨于減弱。

（2）在冰川消融加劇與冰舌波動的共同影響下，巴托拉段冰川融水徑流頻繁遷移改道，冰川融水洪水、冰湖潰決洪水、冰川泥石流等次生災害對公路的威脅正日漸凸顯。

（3）巴托拉段冰川變化對公路最直接的威脅在于2021 年6—7 月巴托拉河新的改道，改道后水流將直接沖擊侵蝕路基，影像公路安全通行。