魏文霞,李真,李亞楠
(1.中國(guó)科學(xué)院西北生態(tài)環(huán)境資源研究院冰凍圈科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,甘肅 蘭州 730000;2.中國(guó)科學(xué)院大學(xué),北京 100049;3.中國(guó)科學(xué)院青藏高原研究所,北京 100101)
冰川厚度和冰儲(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ù)。
七一冰川(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]。
本文使用兩景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
提取七一冰川邊界時(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)算獲得。
七一冰川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
基于實(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
七一冰川冰儲(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]
本文通過(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é)果。