祁连山七一冰川雷达测厚及冰储量估算

魏文霞，李真，李亚楠

（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 km²，冰川朝向为北。七一冰川按照冰川形态来划分，属于冰斗-山谷型冰川；按冰川物理性质划分，属于亚大陆型冰川。冰川最低海拔为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 km²，冰川平均厚度为44.9 m，冰储量为0.1129 km³。

冰川冰面地形数据和冰厚分布数据结合，经过栅格运算，获取的七一冰川冰床地形数据如图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 km³［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 km³，与本文实测结果几乎一致。这说明该模型模拟结果精度高，可能适用于类似于七一冰川的其他冰川冰厚和冰储量估算。图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 km²，冰储量为0.1129 km³，平均冰厚为44.9 m；海拔4 640～4 800 m之间与海拔4 640～4 800 m之间是冰厚值较大的区域。冰储量实测结果与冰厚模型法估算结果一致，而与经验公式法所得结果差异较大。冰厚模型估算法在冰川冰储量估算方面具有很好的应用前景，在未来需要更多冰川的雷达测厚资料，结合冰川规模和冰川流速等物理参数对模型进行优化，获取与实际更符合的模拟结果。