祁連山托來南山6號冰川雷達測厚與冰儲量分析

車 正， 王寧練，3， 梁 倩， 陳安安

（1.陜西省地表系統與環境承載力重點實驗室，陜西 西安 710127；2.西北大學城市與環境學院地表系統與災害研究院，陜西 西安 710127；3.中國科學院青藏高原地球科學卓越創新中心，北京 100101）

0 引言

青藏高原及周邊地區的山地冰川作為“亞洲水塔”的重要組成部分，其變化關系著我國水資源供給與周邊眾多國家的水安全［1-2］。在我國西北干旱地區和中亞干旱區，山地冰川融水對河川徑流起到“削峰填谷”的作用（在干旱年份，冰川消融量增大，對河川徑流的補給增強，從而緩解下游地區的干旱；在濕潤年份，冰川消融減弱，對河川徑流的補給減少），是當前冰川補給型河流徑流量變化相對平緩的重要保障［3］。隨著全球變暖和冰川萎縮加劇，流域尺度的冰川融水峰值出現時間是目前冰川水文研究和流域水資源管理部門最為關注的重大科學問題之一。

冰儲量作為冰川融水徑流模擬的重要參數，其精度影響著冰川融水徑流峰值出現時間的可信度［4］?，F有的大空間尺度冰川儲量多通過面積-體積公式或冰厚估算模型計算［5-11］。面積-體積公式是基于冰儲量與冰川面積間存在良好關系而構建的冪指數模型，目前中國第二次冰川編目數據［12］中冰儲量估算公式就是參考了Radi?等［6］和Grinsted［7］分別在WGI和RGI冰川編目數據基礎上統計分析而得出的經驗公式，國內Liu等［13］根據早期天山地區和祁連山地區的冰川測厚數據提出了適合于相應地區的冰川面積-體積公式。近年來，基于一定物理過程的冰厚估算模型被應用到冰儲量估算研究當中。例如利用淺冰近似理論估算出全球727×103km2冰川的冰儲量約為（140.8±40.4）×103km3［9］，比Farinotti等［10］在2019年基于五種模型平均結果估算的全球705×103km2冰川的冰儲量［（158±41）×103km3］少了約10%；相較Huss等［11］在2012年基于冰川流動定律估算的全球734×103km2冰川冰儲量［（170±21）×103km3］少了約20%。為準確估算冰儲量，Farinotti等［14］結合實測數據對現有17個冰川厚度估算模型進行了對比實驗，結果表明盡管多種模型估算結果的均值能有效提高冰川厚度估算結果的精度，但冰厚估算模型未來改進的關鍵仍是基于實測數據提升模型輸入項的質量。冰川多位于高海拔山區，目前全球約3 000條冰川（占全球冰川總數約1%）進行了冰川測厚［15］，使得現有冰儲量估算模型仍缺乏足夠的實測數據進行模型參數率定，尤其是青藏高原及周邊地區。

冰川測厚可追溯到20世紀20年代，60年代英國科學家Bailey等［16］利用無線電在冰川中具有良好的穿透能力對極地冰蓋進行了探測，隨后Bogorodsky等［17］將雷達無線電回波探測方法用于冰川測厚，使得利用無線電在冰川測厚研究領域得到了系統應用。我國首次冰川測厚是1979年5月用國產KDL-A型礦井地質雷達在祁連山羊龍河冰川開展的［18］，隨后原蘭州冰川凍土研究所研制了專門用于冰川測厚的B-1型雷達，并在天山1號冰川［19］和南極半島柯斯冰帽［20］進行了冰厚探測。21世紀以來隨著探地雷達設備的快速發展，冰川測厚工作已經在全球不同區域的冰川展開［21-36］，積累了一定的冰川厚度實測數據。

本文基于2019年7月利用pulse EKKO PRO探地雷達在托來南山6號冰川的測厚數據，估算了該冰川的冰儲量，分析該冰川的厚度分布和冰床地形狀況。本次工作是托來南山地區第一次開展冰川厚度測量工作，為建立該地區的冰儲量工作估算提供了基礎數據，也為未來在該地區開展的冰川融水徑流模擬等工作奠定基礎。

1 研究區概況

托來南山6號冰川（38.639°N、98.283°E，冰川編碼：G098282E38637N）位于祁連山中段地區的托來南山北坡（圖1），冰川融水最終匯入疏勒河，該冰川是一條小型冰斗-山谷冰川，無表磧物覆蓋，冰面坡度較小，比較平整。依據2019年Landsat-8影像通過目視解譯獲取當年該冰川的總面積為（1.34±0.05）km2，末端海拔為4 586 m，頂端為5 231 m。其中主冰川面積為（1.25±0.04）km2，支冰川面積為（0.09±0.01）km2。第二次中國冰川編目數據［12］顯示該冰川2007年的面積（1.39±0.07）km2，表明該冰川過去十余年間僅萎縮了約3.6%，萎縮速率僅為0.3%·a-1。

2 數據與方法

2.1 數據獲取

為獲取托來南山6號冰川的儲量信息，2019年7月利用加拿大SSI公司生產的pulse EKKO PRO探地雷達（ground penetrating radar，GPR）在該冰川共獲取1條縱測線和8條橫測線，共計216個測點（圖1）?？v測線基本沿冰川中線，分布在海拔4 540～4 940 m，八條橫測線依次分布在海拔約4 620 m（AA′）、4 660 m（BB′）、4 700 m（CC′）、4 720 m（DD′）、4 770 m（EE′）和4 790 m（FF′）、4 850 m（GG′）和4 930 m（HH′）附近。測量時雷達天線頻率設定為100 MHz，發射和接收天線間距為1.5 m，縱剖面測點間隔為15 m，橫剖面測點間隔為20 m。該套雷達系統及設定參數在八一冰川［24］及煤礦冰川［27］都取得了良好的結果，其中八一冰川的探測結果與冰芯長度的誤差僅為1%［24］。測點空間平面坐標用北京合眾思壯公司生產的MG868型手持式GPS進行單點測量，其單點平面定位精度為1.2 m。數字高程模型為ASTER GDEM V3，其空間分辨率為30 m。冰川邊界基于2019年8月14日的Landsat-8 OLI影像經人工目視解譯獲取，計算冰川面積時將投影設置為Albers等面積投影（與第二次中國冰川編目一致）。

圖1 托來南山6號冰川位置及雷達測點分布Fig.1 Location of Tuolainanshan Glacier No.6 and distribution of the ground penetrating radar sounding points

2.2 數據處理

2.2.1 主冰川雷達測厚數據處理

雷達測點資料在EKKO-View Deluxe軟件中處理，將電磁波在冰川中的傳播速度設定為169 m·μs-1，對雷達資料進行可視化后獲取不同斷面上測點的冰厚數據。由于本文的測厚點全部在主冰川上，支冰川上沒有測點分布，因此本文的冰川厚度分布和冰床地形分析都集中在主冰川?；贏rcGIS平臺，利用普通克里格法對主冰川的冰厚數據進行插值，插值時將冰川邊界處的冰厚值設置為0，零值點盡可能均勻分布在主冰川邊界上，主冰川與支冰川交匯處的冰厚設定為支冰川平均厚度的估計值。對插值后的冰厚數據通過柵格化獲取主冰川區內的冰厚分布及冰儲量信息，冰厚數據的空間分辨率保持為30 m（與ASTER GDEM數據一致），進一步利用冰厚柵格數據繪制主冰川區10 m等間距的冰川厚度等值線圖以分析主冰川冰厚分布特征；之后將冰面高程數據（ASTER GDEM）與冰厚數據相減獲得測厚區冰床地形圖。

2.2.2 支冰川平均厚度估算

支冰川因為沒有測厚資料，無法通過插值得到該區域的冰儲量，需尋找適合的方法進行估算。有學者基于半球物理公式提出了估算山地冰川平均厚度的方法［8，37-38］，其表達式為

式中：HF和Hf分別為整條冰川和冰川中線的平均厚度；f為形態因子（山谷冰川通常取值0.8）；ρ為冰密度（900 kg·m-3）；g為重力加速度（9.81 m·s-2）；α為冰川中線的表面坡度；底部剪應力（τ）與冰川作用區的高程差（?H）相關，其計算公式為

2.3 誤差評價

2.3.1 冰川厚度誤差評估

主冰川平均厚度誤差一方面來源于利用GPR進行冰川測厚時產生的測量誤差（εHdata），另一方面來源于利用普通克里格插值進行冰川厚度預測時，預測結果和實測結果之間存在的誤差（εHRMSE）。因此，冰川平均厚度的估算誤差εHˉ可以表示為

Lapazaran等［39］指出雷達脈沖測量誤差和平面定位誤差是基于探地雷達獲取冰川厚度數據誤差的主要來源。本文在托來南山6號冰川共計獲取了216個測點，用所有測點誤差的均值為作為估算冰川平均厚度誤差中的雷達脈沖測量誤差［式（5）］，各個測點誤差（εHdatai）如式（6）所示。

式中：εHGPRi和εHxyi分別為測點i的雷達脈沖測量誤差和平面定位誤差。其中，εHxyi可分為兩部分，一是GPS自身定位誤差，二是當GPS與GPR連接時，GPS的位置信息更新時間和GPR測厚信息獲取時間存在差異而導致的位置偏差。本文未將GPS和GPR連接，二者因信息獲取時間偏差而導致的位置誤差可以忽略，因此平面定位誤差主要來源于GPS儀器自身定位誤差。雷達脈沖測量誤差εHGPR可以表示為

式中：εHc和εHt分別為雷達信號在冰川內傳播速度的誤差和雷達信號在冰川內傳播時間的誤差，二者的計算公式為

式中：t和c分別為雷達信號在冰川中的傳播時間和傳播速度（取中值為169 m·μs-1）；εc和εt分別為二者的誤差，取為雷達天線的頻率（本文為100 MHz）。

利用普通克里格法對測厚點數據進行插值，獲取冰川范圍內空間分辨率為30 m的冰川厚度預測結果（y^i），結合對應位置的實測冰厚（yi）計算二者的均方根誤差作為利用克里格插值法進行冰川厚度估算結果的誤差（εHRMSE）。

2.3.2 主冰川冰儲量誤差評估

本文考慮通過遙感影像獲取冰川面積的誤差和平均冰川厚度的估算誤差，主冰川冰儲量（V）計算公式為

本文利用2019年的Landsat影像通過目視解譯獲取了托來南山6號冰川的面積信息，冰川面積的誤差（εArea）可通過冰川邊緣的像元數計算［12］。

式中：N為冰川邊界經過的像元數；S為Landsat-8經全色波段融合后的像元面積（225 m2）?？紤]誤差傳播，可以將主冰川冰儲量誤差（εV）近似表示為

3 結果與分析

3.1 冰川縱剖面冰川厚度分布特征

托來南山6號冰川的縱測線基本是沿冰川中線布設的，共計有99個測點。圖2為冰川厚度與冰下基巖在縱剖面的變化圖，其中在海拔4 670～4 720 m范圍（測線距離為210～375 m）的雷達測量資料缺失?？v剖面的冰川厚度總體上呈現出比較均勻的分布狀況，整個剖面的平均厚度為（78.61±1.67）m，厚度最大值出現在海拔4 820 m附近，約為（93.83±1.74）m，最小的測點深度也有（45.5±1.54）m?？v剖面冰下地形的起伏狀況整體比較平緩，平均坡度僅為15°，相較于平整的冰川表面，基巖的縱剖面僅在局部區域有微弱起伏。例如，在海拔4 725 m附近有一個微弱的凸起區域，該區域的冰川厚度小于65 m，海拔4 820 m和4 770 m附近有兩個凹陷區域，兩處冰川厚度均在85 m以上。一般來說，山地冰川的運動速度和冰通量均會在冰川中值高度附近達到最大［40］，此處往往冰川厚度也最大。第二次中國冰川編目數據顯示托來南山6號冰川的中值高度約為海拔4 886 m，略高于本文觀測到冰川厚度最大的海拔，其原因可能是該冰川西側山坡的支冰川在海拔4 850 m左右匯入了主冰川，此處山谷對冰川的側向約束變小，支冰川的冰開始向主冰川流動，主冰川受力表現為橫向擠壓和縱向拉伸，使得冰川最厚的地方出現在匯入點下方。

圖2 托來南山6號冰川縱剖面Fig.2 Longitudinal GPR sounding profile of Tuolainanshan Glacier No.6

3.2 冰川橫剖面形態特征

本文通過在托來南山6號冰川共布設的8條橫測線，共計117個有效測點（圖1）。結合GPS實測的高程數據獲取了不同海拔區間的冰下橫剖面（圖3）：冰舌區的幾條橫剖面（AA′、BB′、CC′和DD′）呈現出明顯的U形冰川槽谷形態，且隨著海拔的上升，冰川槽谷底部逐漸寬闊，谷壁逐漸陡峭，橫剖面形狀向梯形靠近；冰川中部EE′和FF′兩個剖面的冰川槽谷盡管仍呈U形谷形態，但谷底部的寬度卻變窄，其中在EE′剖面測得該冰川的最大冰厚［（100.78±1.78）m］，該位置與冰川縱剖面海拔4 770 m處附近的凹陷盆地接近；GG′剖面冰川厚度自西向東呈現出逐漸減小趨勢，其原因可能是該剖面西段位于西側支冰川往主冰川匯入點附近，使得該剖面西段幾個測點的冰厚較大；HH′剖面測點主要位于冰川上部比較平緩區域，因此該剖面冰厚整體呈現出比較平緩，僅在偏西段呈現比較微弱的V形。托來南山6號冰川不同海拔的橫剖面結果整體表明，該冰川的槽谷地形整體呈現出明顯的U形，冰川谷槽的寬度隨海拔升高呈現出先變寬再變窄的趨勢，冰川谷槽寬度在DD′剖面處最大。冰川槽谷是冰川長期作用山谷的結果，冰川的侵蝕能力決定了槽谷的形態，通常冰川規模越大、侵蝕能力越強，冰川槽谷越寬、越深，谷壁越陡［41］。按照施雅風［42］對我國冰川的類型分區，托來南山6號冰川屬亞大陸型冰川，亞大陸型冰川的溫度相對較低，往往與冰床凍結在一起，對冰床的侵蝕較弱，對冰川兩側山坡的侵蝕較強，此時冰川谷槽易形成谷底較寬，谷壁較陡的U形谷，這與本文得到的剖面結果基本一致。

3.3 冰川厚度分布及冰床地形

基于已有的測厚點，利用普通克里格法在主冰川范圍內進行插值，得到主冰川的平均厚度為（39.61±5.32）m，估算結果最大厚度為（100.33±6.03）m（與實測最大厚度相近）。從主冰川的冰川厚度等值線圖（圖4）可以看出：冰厚分布在整體上呈現出自邊緣向中間逐漸增厚的特征；在縱測線附近形成了數個封閉的等值線區域，其中冰舌區中下部形成了2個厚度值為80 m的區域，冰川中部形成了2個厚度為90 m的閉合區；結合冰川形態分析，2個90 m的閉合區域均位于冰川西部支冰川與主冰川交匯處的附近。

圖4 主冰川冰厚等值線Fig.4 Ice thickness contours of the trunk glacier

基于冰厚數據，結合DEM數據獲取主冰川的冰床地形圖（圖5），從等高線的分布可看出，冰床地形整體上為中間低兩邊高的典型槽谷地形，且以4 800 m等高線為分界線，上部冰床槽谷的寬深比相較于下部冰舌區域冰床的更大，也就是上部槽谷比較淺，谷底寬闊，中下部槽谷較深，谷底相對較窄。另外，在冰川東側的山坡區域，冰床地形與冰川表面地形基本一致，且冰床等高線近似平行，說明此處受到的冰川侵蝕作用小，這與該區域冰川厚度普遍較小相吻合。通過冰面與冰床地形的對比可看出，冰川表面地形與冰床地形的對應關系比較明顯，這說明冰床地形相對比較規整，僅在冰川末端形成了閉合等值線，這里隨著后期冰川進一步退縮有形成冰湖的潛力。

圖5 主冰川冰床地形Fig.5 Topography of ice bed of the trunk glacier

3.4 冰儲量估算

結合冰川矢量邊界，對插值后的雷達測厚區的冰厚柵格數據進行二重積分，獲得測厚區冰儲量為（0.0495±0.0082）km3。支冰川作用區的高程差為0.17 km，用式（3）計算出支冰川的底部剪應力（τ）為26.99 kPa，支冰川中線處的平均坡度為16.9°，用式（2）計算出支冰川中線的平均厚度為13.14 m，支冰川的平均厚度僅有10.32 m，支冰川儲量的估算結果為（0.0009±0.0001）km3。將雷達測厚區的主冰川冰儲量與支冰川的冰儲量相加，得到托來南山6號冰川的總冰儲量為（0.0504±0.0082）km3。利用第二次中國冰川編目兩組面積-體積公式［12］估算出托來南山6號冰川主冰川2019年的儲量分別0.0496 km3和0.0577 km3，其結果表明其中一個公式的估算結果與基于實測數據估算結果比較接近，另外一組則有16.7%的差異。

4 討論

冰儲量是描述冰川水資源狀況最直接的參數，青藏高原及其周邊地區冰儲量的準確估算對評估該地區冰川融水資源具有重要意義。因此，冰儲量估算一直是青藏高原及周邊地區冰川變化研究的核心問題，該區域現有的冰儲量主要依靠經驗公式估算。例如我國第一次冰川編目用于估算冰川平均厚度的經驗公式就是依據27條冰川測厚結果建立的［42］，而第二次中國冰川編目中的冰儲量估算則分別引用了Radi?等［6］總結的適用于WGI冰川編目中山地冰川的式（4）和Grinsted［7］在RGI冰川編目數據基礎上經統計分析得出的估算全球冰儲量的式（5）。此外，Liu等［13］根據早期天山地區和祁連山地區的冰川測厚數據提出了適合于中國地區的冰川面積-體積公式，為

式中：V為冰川冰儲量（km3）；A為冰川面積（km2）。然而，Bahr等［5］指出面積-體積公式通常適用于大區域的冰儲量估算，當被應用于單條冰川計算時其準確度僅能與真實結果保持在同一數量級，且不同類型冰川的面積-體積公式參數也存在一定差異。Grinsted［7］在進行全球冰儲量估算時分別對冰川和冰蓋，以及不同規模的冰川面積-體積公式進行了構建。為了嘗試提升冰川面積-體積經驗公式在青藏高原地區的適用性，收集了青藏高原及周邊地區12條已發表的測厚冰川資料（表1），嘗試分析冰川類型對利用面積-體積公式進行冰儲量估算的影響。

利用已有的面積-體積公式對12條已測厚冰川進行冰儲量估算，與實測數據對比結果表明盡管二者處于同一數量級，但單條冰川的誤差存在較大差異。例如在八一冰川、七一冰川、烏魯木齊河源1號冰川和羌塘1號冰川等冰川的估算誤差僅在10%左右，而在抗物熱冰川和四工河4號冰川的估算誤差約為1倍。冰川類型通常是由冰川規模、運動及其所處地形共同決定。以下嘗試依據冰川類型來進行面積和體積關系的分析。表1的整體結果顯示，已有的3組經驗公式對冰帽型（平頂）冰川（八一冰川）的估算誤差均較小，其后依次為冰斗-山谷型冰川和冰斗冰川，其中懸冰川（抗物熱冰川）的估算誤差相對較大。據此，本文分別對冰斗型和冰斗-山谷型兩類實測數據相對較多的冰川的面積-體積關系進行了擬合（懸冰川和冰帽型冰川各自僅有1條，未擬合）。結果如圖6所示，兩種形態冰川的面積-體積擬合曲線有明顯的差異，即冰斗型冰川面積與體積之間的關聯程度要明顯要小于冰斗-山谷型冰川，且5條冰斗-山谷型冰川的擬合關系明顯要優于冰斗型冰川。利用兩種形態的擬合結果，分別對表1中的冰斗型冰川和冰斗-山谷型冰川的冰儲量進行了估算。結果表明利用分形態擬合公式估算的冰儲量誤差整體比之前的幾種方法都有明顯的降低（表1），這說明分類型進行冰川面積-體積公式在提升單條冰川冰儲量估算精度領域存在一定潛力。但由于樣本量太少，使得驗證的結果缺乏足夠的說服力，后續還應該通過增加測厚冰川的數量來提升分類型面積-體積公式的適用性。

表1 青藏高原及周邊地區部分冰川雷達實測冰儲量與經驗公式估算值對比Table 1 Ice volumes，sounded by GPR，of some glaciers on the Tibetan Plateau and surrounding areas，compared with the ice volumes estimated from volume-area scaling of glaciers

圖6 基于實測數據的冰斗型和冰斗-山谷型冰川面積-體積擬合關系Fig.6 Fitting relationship between area and volume for the cirque glaciers and the cirque-valley glaciers based on surveyed data

5 結論

冰川測厚是目前獲取精確冰川儲量的關鍵，本文利用探地雷達對托來南山6號主冰川進行了系統的冰厚測量，在此基礎上結合普通克里格法和半球物理公式，揭示了該冰川的最大厚度為（100.78±1.78）m，主冰川的平均厚度為（39.61±5.32）m，支冰川的平均厚度為10.32 m，進而估算出整條冰川的冰儲量為（0.0504±0.0082）km3。主冰川繪的厚度等值線圖顯示該冰川冰厚呈自邊緣向中間逐漸增厚的分布特征，冰川中部有兩個閉合區域的冰厚超過了90 m。主冰川冰床縱剖面地形整體比較平緩，僅在海拔4 820 m和4 770 m附近有兩個凹陷區域，橫剖面形態呈典型的U形，且隨著海拔升高冰川槽谷的寬度呈現出先變寬再變窄的趨勢，冰床地形整體呈中間低兩邊高的典型槽谷地形，冰川上部的槽谷比較淺，谷底寬闊，中下部槽谷較深，谷底相對較窄。面積-體積公式是目前青藏高原及周邊地區大范圍冰儲量估算的有效手段之一，本文初步探索結果顯示出分類型進行擬合具有降低單條冰川儲量結果誤差的潛力。