祁連山七一冰川雷達(dá)測(cè)厚及冰儲(chǔ)量估算

魏文霞，李真，李亞楠

（1.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，甘肅 蘭州 730000；2.中國(guó)科學(xué)院大學(xué)，北京 100049；3.中國(guó)科學(xué)院青藏高原研究所，北京 100101）

0 引言

冰川厚度和冰儲(chǔ)量是冰川的重要屬性，是現(xiàn)代冰川學(xué)研究中的重要內(nèi)容，同時(shí)也是冰川水文模型［1］、冰川災(zāi)害評(píng)估［2］、冰川動(dòng)力模型［3］研究中的重要輸入?yún)?shù)。目前，冰川冰儲(chǔ)量的研究主要集中在冰量變化方面，但對(duì)冰川冰儲(chǔ)量的估算研究較少。冰川冰儲(chǔ)量的估算主要有經(jīng)驗(yàn)公式法［4-6］、冰厚模型估算法［7-8］和探地雷達(dá)法［9-11］，除經(jīng)驗(yàn)公式法外，冰厚模型法和探地雷達(dá)法都是通過(guò)獲取冰川厚度來(lái)估算冰川冰儲(chǔ)量。探地雷達(dá)法是目前獲取冰厚數(shù)據(jù)最為準(zhǔn)確的方法［12］，使用探地雷達(dá)對(duì)冰川厚度以及冰床地形進(jìn)行探測(cè)，可為冰川內(nèi)部沉積層位、冰川結(jié)構(gòu)、冰下河等方面研究提供豐富可靠的數(shù)據(jù)［13］。早在20世紀(jì)20年代，國(guó)外就已經(jīng)開(kāi)始使用探地雷達(dá)測(cè)量冰川厚度［14］。在20世紀(jì)60年代Bailey等［15］提出無(wú)線(xiàn)電回波探測(cè)方法后，探地雷達(dá)測(cè)量技術(shù)被更廣泛地應(yīng)用到冰厚測(cè)量中。我國(guó)探地雷達(dá)在冰川學(xué)中的應(yīng)用始于20世紀(jì)80年代，1980年中國(guó)科學(xué)院蘭州冰川凍土研究所研制了B-1型冰川專(zhuān)用測(cè)厚雷達(dá)，并在天山烏魯木齊河源1、3號(hào)冰川上試驗(yàn)成功［16］。近年來(lái)，隨著探地雷達(dá)技術(shù)的不斷發(fā)展，冰川探地雷達(dá)在天山、喜馬拉雅山、昆侖山等地區(qū)的冰川上得到廣泛應(yīng)用［17-24］。本文利用2015年8月七一冰川探地雷達(dá)測(cè)厚數(shù)據(jù)，對(duì)冰川的冰儲(chǔ)量進(jìn)行估算，并繪制冰川的冰厚分布圖和冰川冰下地形圖，為冰川動(dòng)力學(xué)模擬提供基礎(chǔ)數(shù)據(jù)。

1 研究區(qū)概況

七一冰川（GLIMS ID：G097755E39237N［25］；第一次冰川編目編碼5Y437C18［26］）位于青藏高原北部祁連山脈中段托賴(lài)山北坡，冰川融水注入北大河流域柳溝泉河（屬河西內(nèi)流水系）。七一冰川是我國(guó)開(kāi)展現(xiàn)代冰川綜合考察與研究的第一條冰川［27］，也是黑河流域內(nèi)唯一一條具有較長(zhǎng)時(shí)間觀測(cè)序列的冰川。根據(jù)中國(guó)第二次冰川編目［25］，2006年七一冰川面積為2.53 km2，冰川朝向?yàn)楸?。七一冰川按照冰川形態(tài)來(lái)劃分，屬于冰斗-山谷型冰川；按冰川物理性質(zhì)劃分，屬于亞大陸型冰川。冰川最低海拔為4 314 m，最高海拔為5 145 m，海拔跨度830 m，冰面平均坡度為20°，冰舌部分較為平坦，冰川后壁較陡峭，坡度跨度大。冰川積累區(qū)上部發(fā)育有東、中、西三個(gè)粒雪盆，其中東粒雪盆和西粒雪盆較為寬闊，中粒雪盆較小［28］。觀測(cè)結(jié)果顯示，七一冰川近年來(lái)冰川物質(zhì)收支以負(fù)平衡為主，末端退縮，冰川面積呈持續(xù)減小態(tài)勢(shì)［29］。

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

2.1 數(shù)據(jù)來(lái)源

本文使用兩景Landsat 8遙感影像數(shù)據(jù)和一幅DSM地形數(shù)據(jù)獲取七一冰川邊界和冰下地形。遙感影像數(shù)據(jù)下載自美國(guó)地質(zhì)調(diào)查局網(wǎng)站（USGS，https：//glovis.usgs.gov），影像的軌道號(hào)均為135/033，拍攝時(shí)間分別為2015年3月19日和2015年8月10日。提取冰川邊界時(shí)以2015年3月19日的遙感影像為主，該景影像在研究區(qū)內(nèi)無(wú)云量覆蓋且冰川積雪覆蓋范圍很小，僅冰川上部邊界處有少量積雪覆蓋，在提取此處冰川邊界時(shí)，使用2015年8月10日的遙感影像輔助。地形數(shù)據(jù)下載自日本宇宙航空研究開(kāi)發(fā)機(jī)構(gòu)網(wǎng)站（JAXA，https：//www.eorc.jaxa.jp/ALOS/en/aw3d30/index.htm），水平分辨率30 m。

冰厚數(shù)據(jù)是2015年8月利用pulse EKKO PRO型探地雷達(dá)獲取的實(shí)測(cè)資料。測(cè)量時(shí)雷達(dá)天線(xiàn)的中心頻率為100 MHz，天線(xiàn)間距為2 m，測(cè)點(diǎn)步長(zhǎng)為10 m，采樣時(shí)間間隔為0.8 ns，電磁波在冰川中的傳播速度設(shè)定為0.169 m·ns-1。此套雷達(dá)系統(tǒng)和相應(yīng)的工作參數(shù)曾在八一冰川［23］及煤礦冰川［24］的研究中使用，經(jīng)八一冰川兩根透底冰芯長(zhǎng)度和同位置雷達(dá)實(shí)測(cè)值對(duì)比驗(yàn)證，雷達(dá)測(cè)厚相對(duì)誤差小于1%［23］。在測(cè)量過(guò)程中，同步利用集思寶MG768W手持GPS定位記錄各測(cè)點(diǎn)的位置信息。七一冰川共獲取774個(gè)有效冰厚度數(shù)據(jù)，主要包含9條橫測(cè)線(xiàn)和5條縱測(cè)線(xiàn)（圖1）。

圖1 研究區(qū)位置示意和雷達(dá)測(cè)線(xiàn)分布圖Fig.1 The location of Qiyi Glacier and the distribution of the Ground Penetrating Radar（GPR）sounding lines

2.2 數(shù)據(jù)處理

提取七一冰川邊界時(shí)，首先把經(jīng)過(guò)預(yù)處理的遙感影像的多光譜波段和全色波段數(shù)據(jù)進(jìn)行了圖像融合處理，使影像分辨率從30 m提高到了15 m，然后裸冰區(qū)冰川邊界利用影像紅色波段與短波紅外波段的比值閾值來(lái)提取［30］，結(jié)合人工目視解譯進(jìn)行修正，冰川分冰嶺處邊界則基于地形數(shù)據(jù)利用Arc-GIS軟件平臺(tái)中的水文分析模塊進(jìn)行提?。?1］，最后將提取的兩部分冰川邊界合并后進(jìn)行線(xiàn)平滑處理。

雷達(dá)數(shù)據(jù)使用與雷達(dá)系統(tǒng)配套的EKKO-View Deluxe軟件處理。將雷達(dá)數(shù)據(jù)進(jìn)行可視化處理后，再進(jìn)行AGC（Automatic Gain Control）增益調(diào)節(jié)，使冰巖界面清晰（圖2）。假定雷達(dá)的波速為0.169 m·ns-1，將雷達(dá)信號(hào)的雙程走時(shí)轉(zhuǎn)化為冰厚度，逐個(gè)測(cè)點(diǎn)讀取冰川厚度數(shù)據(jù)。

圖2 七一冰川部分探地雷達(dá)圖像Fig.2 Example of GPR images sounded through Qiyi Glacier

冰川冰厚分布是基于ArcGIS軟件平臺(tái)，通過(guò)空間插值運(yùn)算獲取。具體插值步驟為：（1）把測(cè)點(diǎn)的實(shí)測(cè)經(jīng)度、緯度和冰厚數(shù)據(jù)導(dǎo)入軟件中，生成冰厚點(diǎn)圖層；（2）在冰川上部邊界選取若干點(diǎn)，冰厚值設(shè)為零，加入冰厚點(diǎn)圖層；（3）地形數(shù)據(jù)按冰川邊界裁剪并生成冰川表面坡度圖，獲取坡度數(shù)據(jù)；（4）以冰厚點(diǎn)圖層作為主要變量，坡度數(shù)據(jù)作為協(xié)變量［18］進(jìn)行協(xié)同克里金空間插值運(yùn)算；（5）插值結(jié)果按2015年七一冰川邊界裁剪，得到冰川冰厚分布圖?；诒穹植紙D利用厚度積分法估算2015年七一冰川冰儲(chǔ)量。冰川冰下地形數(shù)據(jù)是利用冰川表面地形數(shù)據(jù)和冰厚分布柵格圖進(jìn)行柵格運(yùn)算獲得。

3 結(jié)果分析

3.1 冰川測(cè)線(xiàn)冰厚分布特征

七一冰川2015年實(shí)測(cè)冰厚最大是115 m，出現(xiàn)在海拔4 710 m處，冰厚最小值是6 m，出現(xiàn)在海拔4 768 m處。七一冰川縱測(cè)線(xiàn)共布設(shè)了5條，利用測(cè)點(diǎn)高程和冰厚數(shù)據(jù)，繪制縱測(cè)線(xiàn)剖面示意圖（圖3）?？梢钥闯?，縱測(cè)線(xiàn)冰川厚度隨海拔升高逐漸增大，至冰川中部粒雪盆區(qū)域達(dá)到實(shí)測(cè)冰厚最大值，在粒雪盆上方，冰川厚度隨海拔升高逐漸減小，但冰面地貌和冰床地形不盡一致。

圖3 七一冰川雷達(dá)測(cè)厚縱剖面Fig.3 Longitudinal GPR sounding profiles of Qiyi Glacier

圖4是9條橫測(cè)線(xiàn)剖面示意圖，橫剖面上冰川表面平均坡度和地形起伏差異較大。其中，較為完整的測(cè)線(xiàn)BB′和CC′槽谷呈“U”形，DD′和EE′槽谷呈“V”形，說(shuō)明由上至下冰川槽谷谷底越來(lái)越寬闊，谷壁越來(lái)越平緩，逐漸由“V”形向“U”形轉(zhuǎn)變。據(jù)推測(cè)，其原因可能是在冰川向下運(yùn)動(dòng)時(shí)，運(yùn)動(dòng)的塑性冰川對(duì)底部巖塊進(jìn)行磨蝕與拔蝕，冰床兩壁上的巖石經(jīng)過(guò)凍融作用后也變得松散、易崩塌，冰川不斷下蝕與展寬，冰川槽谷兩側(cè)谷壁慢慢變得平直，槽谷形態(tài)逐漸轉(zhuǎn)變?yōu)椤癠”形形態(tài)。測(cè)線(xiàn)FF′和II′槽谷均出現(xiàn)兩個(gè)明顯凹槽，呈現(xiàn)復(fù)式槽谷特征。一般來(lái)說(shuō)，在冰床基巖軟硬特征確定的情況下，冰床如果遭受多個(gè)不同方向冰流的侵蝕，冰床截面就可能形成復(fù)式槽谷。測(cè)線(xiàn)FF′的槽谷主要因西粒雪盆和中粒雪盆兩股冰流侵蝕而成，測(cè)線(xiàn)II′槽谷則是中粒雪盆和東粒雪盆兩處冰流的侵蝕結(jié)果。

圖4 七一冰川雷達(dá)測(cè)厚橫剖面Fig.4 Transverse GPR sounding profiles of Qiyi Glacier

值得注意的是，在冰川消融區(qū)，臨近冰川邊界的測(cè)點(diǎn)，實(shí)測(cè)冰厚值都較大（表1）。測(cè)線(xiàn)AA′最東側(cè)測(cè)點(diǎn)距離冰川邊界10.3 m，實(shí)測(cè)冰厚值為60 m；測(cè)線(xiàn)DD′最東側(cè)測(cè)點(diǎn)距離邊界僅1.3 m，對(duì)應(yīng)的實(shí)測(cè)冰厚值為51 m。實(shí)際上，在消融期，七一冰川積累區(qū)冰川邊界處冰川冰和基巖相接，但在消融區(qū)，邊界冰體與基巖基本被水流分離，冰體明顯下切呈冰崖?tīng)?，邊界處冰厚并不為零。如果定義冰川消融區(qū)冰川邊界處冰厚值為零，顯然與實(shí)際不符。因此，本文在冰厚插值時(shí)，冰川消融區(qū)邊界處未設(shè)零值。大部分冰川厚度模型都假定冰川消融區(qū)邊界處的冰厚值為零［32］，若能在模型中對(duì)冰川消融區(qū)和積累區(qū)邊界處冰厚進(jìn)行不同定義，模擬結(jié)果可能更符合實(shí)際。

表1 臨近冰川邊界測(cè)點(diǎn)的冰厚Table 1 The ice thickness that measured near the glacier boundary

3.2 冰川冰厚分布情況和冰床地形

基于實(shí)測(cè)冰厚數(shù)據(jù)，利用空間插值運(yùn)算獲得七一冰川冰厚分布如圖5（a）所示。冰川厚度沿主流線(xiàn)向東西兩側(cè)逐漸減薄，東粒雪盆下方海拔4 650～4 700 m范圍是整條冰川冰厚值中心，平均厚度是96 m；在海拔4 470～4 560 m之間冰川消融區(qū)也存在一個(gè)冰厚較大區(qū)域，平均厚度是84 m；冰川海拔4 850 m以上的區(qū)域冰厚值較小，平均厚度僅有22 m。2015年七一冰川面積為2.517 km2，冰川平均厚度為44.9 m，冰儲(chǔ)量為0.1129 km3。

冰川冰面地形數(shù)據(jù)和冰厚分布數(shù)據(jù)結(jié)合，經(jīng)過(guò)柵格運(yùn)算，獲取的七一冰川冰床地形數(shù)據(jù)如圖5（b）所示。七一冰川冰床海拔高度在4 246～5 144 m之間。冰川東部海拔4 550～4 850 m之間發(fā)育有西北方向的溝槽，冰川西側(cè)海拔4 550～4 700 m內(nèi)發(fā)育有東北方向的溝槽。兩溝槽在海拔4 450～4 500 m處匯聚，因冰川對(duì)底部基巖侵蝕作用加強(qiáng)，導(dǎo)致該處形成較大范圍的圍椅狀洼地。冰床地形與冰川表面特征在海拔4 900 m以上幾乎一致，在海拔4 350～4 650 m之間差異較大。海拔下降到4 246～4 350 m，冰床和冰面地形特征又趨于相似，冰厚分布較均勻。

圖5 冰厚分布和冰床地形圖Fig.5 Ice thickness distribution and bed topography of Qiyi Glacier

4 討論

七一冰川冰儲(chǔ)量估算誤差可能來(lái)自以下幾個(gè)方面：（1）探地雷達(dá)系統(tǒng)設(shè)置的參數(shù)影響冰川厚度測(cè)量精度［33］。本文使用的探地雷達(dá)系統(tǒng)和參數(shù)實(shí)際測(cè)量的相對(duì)誤差小于1%［23］。（2）插值時(shí)，地形數(shù)據(jù)的精度對(duì)冰厚插值結(jié)果會(huì)產(chǎn)生一定的影響。JAXA發(fā)布的DSM數(shù)據(jù)集垂直精度±5 m，符合插值方法的要求。（3）測(cè)點(diǎn)空間分布情況和空間插值方法影響冰厚插值結(jié)果。七一冰川測(cè)點(diǎn)并未完全覆蓋整個(gè)冰川，尤其是冰川上部缺乏冰厚測(cè)量值，這可能導(dǎo)致冰厚插值結(jié)果較小，冰儲(chǔ)量估算值也會(huì)略小于真實(shí)值。本文實(shí)測(cè)點(diǎn)位置插值所得冰厚值和實(shí)測(cè)值相比，平均相對(duì)誤差僅為1.25%，說(shuō)明插值方法是可行的。本文的冰儲(chǔ)量估算誤差在誤差容許范圍內(nèi)，結(jié)果是可靠的。此外，本文使用的DSM數(shù)據(jù)集是JAXA發(fā)布的最新版本，該數(shù)據(jù)代表的為研究區(qū)2011年的地形，與探地雷達(dá)觀測(cè)時(shí)間相差4年。受冰川消融作用的影響，4年內(nèi)七一冰川表面高程會(huì)有所下降，用2011年地形數(shù)據(jù)進(jìn)行柵格運(yùn)算繪制七一冰川冰床地形，會(huì)導(dǎo)致冰床高程比實(shí)際值偏高。根據(jù)七一冰川物質(zhì)平衡觀測(cè)資料，2011—2015年期間，七一冰川累積物質(zhì)平衡為-1 840 mm w.e.［29］，相當(dāng)于冰川冰面減薄2.05 m，小于地形數(shù)據(jù)的垂直精度±5 m，說(shuō)明使用2011年地形數(shù)據(jù)引起的冰川冰床地形誤差在可接受范圍內(nèi)。

根據(jù)中國(guó)第二次冰川編目中使用的冰儲(chǔ)量估算公式計(jì)算七一冰川的冰儲(chǔ)量是0.130 km3［25］，與實(shí)測(cè)值相比，相對(duì)誤差為15.02%。表2匯集了當(dāng)前青藏高原具有完整雷達(dá)測(cè)厚資料的冰川，可以看出：大部分冰川冰儲(chǔ)量實(shí)測(cè)值和經(jīng)驗(yàn)公式計(jì)算結(jié)果有較大差異。這進(jìn)一步說(shuō)明適用于區(qū)域尺度的冰儲(chǔ)量-面積經(jīng)驗(yàn)公式應(yīng)用到單條冰川的冰儲(chǔ)量估算時(shí)有其局限性，方法帶來(lái)的誤差不可忽視。未來(lái)仍然需要積累大量冰川雷達(dá)厚度觀測(cè)資料，通過(guò)分析冰川規(guī)模，冰川形態(tài)和冰川所處地域等參數(shù)條件，進(jìn)一步優(yōu)化冰川冰儲(chǔ)量估算公式，提高單條冰川的冰儲(chǔ)量估算精度。

表2 青藏高原部分冰川探地雷達(dá)和經(jīng)驗(yàn)公式所得冰儲(chǔ)量對(duì)比Table 2 Comparison of ice volumes obtained by GPR and empirical formula in some glaciers on the Tibetan Plateau

Farinotti等［34］通過(guò)集成多個(gè)冰厚模型的結(jié)果，發(fā)布了全球冰川冰厚數(shù)據(jù)集。在此數(shù)據(jù)集中，七一冰川平均冰厚為44.6 m，冰儲(chǔ)量為0.1129 km3，與本文實(shí)測(cè)結(jié)果幾乎一致。這說(shuō)明該模型模擬結(jié)果精度高，可能適用于類(lèi)似于七一冰川的其他冰川冰厚和冰儲(chǔ)量估算。圖6是該模型模擬的冰厚分布圖［34］，在細(xì)節(jié)方面與實(shí)測(cè)結(jié)果存在一些差異。模擬最大冰厚值為84 m，比實(shí)測(cè)最大冰厚值小31 m。模型模擬的七一冰川最大冰厚區(qū)出現(xiàn)在海拔4 550～4 570 m之間，與本文實(shí)測(cè)結(jié)果（圖5）分布不同。如果對(duì)該模型的輸入?yún)?shù)進(jìn)行細(xì)微調(diào)整，可能會(huì)獲得與實(shí)際情況更符合的結(jié)果。

圖6 冰厚數(shù)據(jù)庫(kù)中七一冰川冰厚分布圖［34］Fig.6 Ice thickness distribution of Qiyi Glacier derived by inversion method［34］

5 結(jié)論

本文通過(guò)對(duì)2015年七一冰川實(shí)測(cè)GPR冰川厚度數(shù)據(jù)進(jìn)行協(xié)同克里金插值運(yùn)算，繪制出了冰川冰厚分布和冰床地形圖，得出其冰儲(chǔ)量。七一冰川2015年面積為2.517 km2，冰儲(chǔ)量為0.1129 km3，平均冰厚為44.9 m；海拔4 640～4 800 m之間與海拔4 640～4 800 m之間是冰厚值較大的區(qū)域。冰儲(chǔ)量實(shí)測(cè)結(jié)果與冰厚模型法估算結(jié)果一致，而與經(jīng)驗(yàn)公式法所得結(jié)果差異較大。冰厚模型估算法在冰川冰儲(chǔ)量估算方面具有很好的應(yīng)用前景，在未來(lái)需要更多冰川的雷達(dá)測(cè)厚資料，結(jié)合冰川規(guī)模和冰川流速等物理參數(shù)對(duì)模型進(jìn)行優(yōu)化，獲取與實(shí)際更符合的模擬結(jié)果。