祁連山七一冰川雷達測厚及冰儲量估算

魏文霞，李真，李亞楠

（1.中國科學院西北生態環境資源研究院冰凍圈科學國家重點實驗室，甘肅 蘭州 730000；2.中國科學院大學，北京 100049；3.中國科學院青藏高原研究所，北京 100101）

0 引言

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

1 研究區概況

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

2 數據與方法

2.1 數據來源

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

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

圖1 研究區位置示意和雷達測線分布圖Fig.1 The location of Qiyi Glacier and the distribution of the Ground Penetrating Radar（GPR）sounding lines

2.2 數據處理

提取七一冰川邊界時，首先把經過預處理的遙感影像的多光譜波段和全色波段數據進行了圖像融合處理，使影像分辨率從30 m提高到了15 m，然后裸冰區冰川邊界利用影像紅色波段與短波紅外波段的比值閾值來提取［30］，結合人工目視解譯進行修正，冰川分冰嶺處邊界則基于地形數據利用Arc-GIS軟件平臺中的水文分析模塊進行提取［31］，最后將提取的兩部分冰川邊界合并后進行線平滑處理。

雷達數據使用與雷達系統配套的EKKO-View Deluxe軟件處理。將雷達數據進行可視化處理后，再進行AGC（Automatic Gain Control）增益調節，使冰巖界面清晰（圖2）。假定雷達的波速為0.169 m·ns-1，將雷達信號的雙程走時轉化為冰厚度，逐個測點讀取冰川厚度數據。

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

冰川冰厚分布是基于ArcGIS軟件平臺，通過空間插值運算獲取。具體插值步驟為：（1）把測點的實測經度、緯度和冰厚數據導入軟件中，生成冰厚點圖層；（2）在冰川上部邊界選取若干點，冰厚值設為零，加入冰厚點圖層；（3）地形數據按冰川邊界裁剪并生成冰川表面坡度圖，獲取坡度數據；（4）以冰厚點圖層作為主要變量，坡度數據作為協變量［18］進行協同克里金空間插值運算；（5）插值結果按2015年七一冰川邊界裁剪，得到冰川冰厚分布圖。基于冰厚分布圖利用厚度積分法估算2015年七一冰川冰儲量。冰川冰下地形數據是利用冰川表面地形數據和冰厚分布柵格圖進行柵格運算獲得。

3 結果分析

3.1 冰川測線冰厚分布特征

七一冰川2015年實測冰厚最大是115 m，出現在海拔4 710 m處，冰厚最小值是6 m，出現在海拔4 768 m處。七一冰川縱測線共布設了5條，利用測點高程和冰厚數據，繪制縱測線剖面示意圖（圖3）。可以看出，縱測線冰川厚度隨海拔升高逐漸增大，至冰川中部粒雪盆區域達到實測冰厚最大值，在粒雪盆上方，冰川厚度隨海拔升高逐漸減小，但冰面地貌和冰床地形不盡一致。

圖3 七一冰川雷達測厚縱剖面Fig.3 Longitudinal GPR sounding profiles of Qiyi Glacier

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

圖4 七一冰川雷達測厚橫剖面Fig.4 Transverse GPR sounding profiles of Qiyi Glacier

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

表1 臨近冰川邊界測點的冰厚Table 1 The ice thickness that measured near the glacier boundary

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

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

冰川冰面地形數據和冰厚分布數據結合，經過柵格運算，獲取的七一冰川冰床地形數據如圖5（b）所示。七一冰川冰床海拔高度在4 246～5 144 m之間。冰川東部海拔4 550～4 850 m之間發育有西北方向的溝槽，冰川西側海拔4 550～4 700 m內發育有東北方向的溝槽。兩溝槽在海拔4 450～4 500 m處匯聚，因冰川對底部基巖侵蝕作用加強，導致該處形成較大范圍的圍椅狀洼地。冰床地形與冰川表面特征在海拔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 討論

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

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

表2 青藏高原部分冰川探地雷達和經驗公式所得冰儲量對比Table 2 Comparison of ice volumes obtained by GPR and empirical formula in some glaciers on the Tibetan Plateau

Farinotti等［34］通過集成多個冰厚模型的結果，發布了全球冰川冰厚數據集。在此數據集中，七一冰川平均冰厚為44.6 m，冰儲量為0.1129 km3，與本文實測結果幾乎一致。這說明該模型模擬結果精度高，可能適用于類似于七一冰川的其他冰川冰厚和冰儲量估算。圖6是該模型模擬的冰厚分布圖［34］，在細節方面與實測結果存在一些差異。模擬最大冰厚值為84 m，比實測最大冰厚值小31 m。模型模擬的七一冰川最大冰厚區出現在海拔4 550～4 570 m之間，與本文實測結果（圖5）分布不同。如果對該模型的輸入參數進行細微調整，可能會獲得與實際情況更符合的結果。

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

5 結論

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