祁连山托来南山6号冰川雷达测厚与冰储量分析

车 正， 王宁练，3， 梁 倩， 陈安安

（1.陕西省地表系统与环境承载力重点实验室，陕西 西安 710127；2.西北大学城市与环境学院地表系统与灾害研究院，陕西 西安 710127；3.中国科学院青藏高原地球科学卓越创新中心，北京 100101）

0 引言

青藏高原及周边地区的山地冰川作为“亚洲水塔”的重要组成部分，其变化关系着我国水资源供给与周边众多国家的水安全［1-2］。在我国西北干旱地区和中亚干旱区，山地冰川融水对河川径流起到“削峰填谷”的作用（在干旱年份，冰川消融量增大，对河川径流的补给增强，从而缓解下游地区的干旱；在湿润年份，冰川消融减弱，对河川径流的补给减少），是当前冰川补给型河流径流量变化相对平缓的重要保障［3］。随着全球变暖和冰川萎缩加剧，流域尺度的冰川融水峰值出现时间是目前冰川水文研究和流域水资源管理部门最为关注的重大科学问题之一。

冰储量作为冰川融水径流模拟的重要参数，其精度影响着冰川融水径流峰值出现时间的可信度［4］。现有的大空间尺度冰川储量多通过面积-体积公式或冰厚估算模型计算［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）。纵测线基本沿冰川中线，分布在海拔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，对雷达资料进行可视化后获取不同断面上测点的冰厚数据。由于本文的测厚点全部在主冰川上，支冰川上没有测点分布，因此本文的冰川厚度分布和冰床地形分析都集中在主冰川。基于ArcGIS平台，利用普通克里格法对主冰川的冰厚数据进行插值，插值时将冰川边界处的冰厚值设置为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）的雷达测量资料缺失。纵剖面的冰川厚度总体上呈现出比较均匀的分布状况，整个剖面的平均厚度为（78.61±1.67）m，厚度最大值出现在海拔4 820 m附近，约为（93.83±1.74）m，最小的测点深度也有（45.5±1.54）m。纵剖面冰下地形的起伏状况整体比较平缓，平均坡度仅为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形，且随着海拔升高冰川槽谷的宽度呈现出先变宽再变窄的趋势，冰床地形整体呈中间低两边高的典型槽谷地形，冰川上部的槽谷比较浅，谷底宽阔，中下部槽谷较深，谷底相对较窄。面积-体积公式是目前青藏高原及周边地区大范围冰储量估算的有效手段之一，本文初步探索结果显示出分类型进行拟合具有降低单条冰川储量结果误差的潜力。