基于多源数据的天山乌鲁木齐河源1号冰川变化研究

李宏亮， 王璞玉，， 李忠勤， 王盼盼， 徐春海， 刘爽爽，金 爽，张正勇，徐丽萍

（1.中国科学院西北生态环境资源研究院冰冻圈科学国家重点实验室/天山冰川观测试验站，甘肃兰州730000；2.中国科学院大学，北京100049；3.石河子大学理学院，新疆石河子832003；4.西北师范大学地理与环境科学学院，甘肃兰州730070）

0 引言

冰川是气候变化的敏感指示器，也是全球水循环的重要组成部分［1］。在气候变暖的影响下，全球大多数冰川呈现普遍变薄和退缩趋势，对水资源、水循环和生态环境等都产生了显著影响［2-3］，导致冰崩、冰湖溃决等灾变风险的增加［4-6］。过去50年间全球气温升高，我国天山冰川面积和冰储量均呈持续退缩趋势，其面积和冰储量年均退缩率分别为每年-0.7%和每年0.83%，且物质损失严重，减少了27%，同时1942—2014年间末端年均退缩速率为1.56 m·a-1［7-9］，短期内可使河川径流量增加，但在长时间尺度上这种影响可能使河川径流总量减少［10］，在冰川面积减少、厚度变薄及平衡线海拔升高等因素和机制的促进下，天山中部典型流域径流量自20世纪90年代中期后递减［11］，体现了气候变化驱动下以冰川为代表的固态水体对流域水资源的调控机制。因此，开展天山地区冰川变化与气候变化的联系至关重要，对提高天山冰川变化的认识和水资源利用与管理具有重要意义。

基于各冰川参数的变化可评估气候变化及管理水资源［12-13］。冰川物质平衡是表征冰川积累和消融量值的重要冰川参数之一，对气候变化响应敏感。物质平衡及其动态变化是引起冰川规模和径流变化的物质基础，是连结冰川与气候、冰川与水资源的重要纽带。面积和长度（末端）是冰川的重要几何形态参数，在冰川消融模拟、冰川几何形态变化模拟、冰川水文等各类研究中，均充当不可或缺的基础参数和边界限制条件。通过定期对参照冰川的面积、长度（末端）等参数进行重复观测，不仅可以精细研究冰川变化局部特征，捕捉变化过程，为冰川对气候变化响应及其模拟研究提供数据基础，而且是遥感资料不可或缺的验证数据。

目前国内外在青藏高原及周边地区开展了大量关于冰川变化时空格局及原因和对水资源及海平面变化的影响［14-20］等方面的研究。但仍存在一些问题。如：第一，大尺度冰川变化研究结果精度存在显著差异且缺少实测验证资料；第二，利用单条冰川观测结果可进行实测验证，但连续且详细的监测冰川在全球范围内非常少。因此，对现有定位观测的参照冰川持续观测尤为重要。乌鲁木齐河源1号冰川（简称1号冰川）作为全球定位观测的参照冰川之一，前人已经开展了大量研究［21-26］，随着近年来各种高新技术手段的涌现，使得冰川观测精度大幅提高，如三维激光扫描（Terrestrial Laser Scanner，TLS）灵活且更经济，可获得单条冰川年或季节内高时空分辨率冰面高程信息［27］；无人机航测技术（Unmanned Aerial Vehicle，UAV）具有机动性强、操作简单，低空飞行可获取冰面高分辨率正射影像（DOM）和数字高程模型（DEM），不仅能够清晰解读冰面微地貌特征，还在冰川地表模型构建方面具有明显的优势［28］；载波相位差分定位技术（Real Time Kinematic-Global Position System，RTK-GPS）是目前GPS测量中精度最高的一种定位方法，可实现密集点位覆盖，适合于冰川表面地形高程测量等优势［29］。前期观测1号冰川主要利用RTK-GPS和TLS技术，其中，物质平衡观测结果能够与花杆/雪坑法得到的物质平衡很好的对应，但观测结果仅由单一技术获取的冰面高程得到［24-25］，结合RTK-GPS和TLS技术，本文首次引入UAV技术对1号冰川继续观测，同时试验RTKGPS、TLS和UAV技术之间应用在冰川变化研究中的适用性。因此，本研究将基于不同时期RTKGPS、TLS和UAV技术获取的冰川高程数据，分析1号冰川近期变化，比较三种观测技术之间应用到冰川变化的适用性，同时提供近期准确的冰川变化信息，为下一步分布式冰川能量-物质模拟提供基础。

1 研究区概况

1号冰川（43°06′N，86°49′E）是我国监测时间最长、观测资料连续性最好的一条冰川，位于我国天山中部喀拉乌成山脉主脉北坡乌鲁木齐河河源上游［图1（a）］，属双支冰斗-山谷冰川，由东西两支组成，呈东北走向［图1（b）］。20世纪50年代至今，对该冰川进行连续的物质平衡、末端位置、表面高程等观测，2012年冰川厚度为44.5 m，1980—2012年间冰川厚度减薄速率为0.34 m·a-1，同期末端退缩速率及物质平衡分别为4.4 m·a-1和-0.46 m w.e.·a-1，东西支末端变化差异显著，在过去32年间呈加速退缩趋势［23］。

2 数据与方法

2.1 冰川表面高程、面积和末端变化

为了开展近期1号冰川表面高程、面积和末端变化，本研究结合不同时期数据源开展了对比研究（表1）。

表1 本研究使用数据列表Table 1 Data information used in this study

2012年9月1日利用RTK-GPS获得1号冰川表面高程和末端位置数据。GPS基站固定在临近末端的位置［图1（b）］，接收机在冰面进行同步观测，测量间距为20 m，除冰川上部陡峭部分外，实现了整个冰川的观测。

图1 乌鲁木齐河源1号冰川观测示意图Fig.1 The observation diagram on the Urumqi Glacier No.1：geographical location of the Urumqi Glacier No.1（a），observation sites（b）

2015年4月25日利用三维激光扫描技术获取1号冰川表面点云数据，其中用于确定目标对象三维坐标信息的纵向（ϴ）和横向（φ）扫描角度可在地面激光扫描仪中完成，最终，获得目标对象的三维坐标信息。同时，为了使三维激光扫描仪能够最大可能的扫描到1号冰川全貌，将脉冲频率分别设置为粗扫描（30 kHz）与精扫描（50 kHz），结合RTK-GPS确定的4个扫描站坐标信息［图1（b）］，实现了点云配准以及重复扫描区域比例大于30%的要求［21］，最后进行坐标系统配准和点云数据分类及滤波。

结合RTK-GPS及三维激光扫描仪获取的冰面高程数据，笔者于2018年4月24日使用大疆经纬M 200专业型四旋翼无人机对1号冰川进行了航空摄影测量以获取冰川正射影像和数字表面模型（DSM）。但由于电池供电不足原因，此次航测只获取到冰舌区冰面高程，具体航测范围如图1（b）所示，同时，在冰川基岩处和冰面花杆位置均匀的布设了8个地面控制点（GCPs）。解算软件使用Pix4D mapper，得到高精度的DEM及正射影像。GCPs精度利用平均误差（Emean）与均方根误差（RMSE）来评估［30］，结果表明，无人机航测平均水平误差为0.04 m，均方根误差（RMSE）为0.06 m。

将上述不同时期的冰川观测数据进行重采样（5 m×5 m），并进行坐标归一化处理，均采用统一的UTM投影和WGS84椭球体坐标系统。之后，通过不同时期的数据对比，开展冰川表面高程、面积和末端变化分析研究。

2.2 花杆/雪坑法冰川物质平衡

1号冰川物质平衡观测始于1959年，每年4月底和8月末各观测1次，研究期间均匀布设26～43个花杆，开展花杆/雪坑法冰川物质平衡观测，以期和不同时期冰面高程数据变化，即大地测量法冰川物质平衡结果进行对比。测量包括研究时期内花杆高度、雪坑深、雪类型及密度。单点物质平衡计算方法见式（1），之后利用插值手段，结合等值线法与等高线法，得到观测时段内整条冰川的物质平衡值。

式中：bs、bice、bsi分别为雪、冰川冰以及附加冰的物质平衡，具体可参见已有研究［31］。

花杆/雪坑法得到的冰川物质平衡误差主要来源于野外观测和计算整条冰川物质平衡时的插值，后者主要是将插值方法应用到冰川上未观测的区域中［32］，由于缺少实测数据，本文依据Andreassen等［33］假设的插值造成误差为±0.1 m w.e.·a-1作为本文插值造成的误差。考虑以上因素，整条冰川的物质平衡年误差（σglac）可以利用式（2）计算。

2.3 误差分析

2.3.1 DEM配准及表面高程误差

在计算多源DEM差值时，需要对原始高程数据进行配准和校正。本文采用Nuth等［34］提出的数据配准方法，以2015年DEM为参考，其余时期DEM为错位DEM，选择非冰川区作为感兴趣区，提取两期DEM非冰川区的高程差值dh、2015年DEM的坡度和坡向，通过余弦拟合求得平移矢量来校正DEM间水平和垂直误差。在完成DEM数据配准后，由于空间分辨率的差异，不同DEM数据间还存在高程偏差，且这种偏差在冰川区和非冰川区保持一致，利用地面最大曲率与非冰川区高程差之间的关系进行纠正［35］，本研究基于2015年DEM数据，在ENVI 5.1软件支持下计算像元最大曲率，建立了最大曲率与高程差的关系。

根据不同DEM数据间的高程残差误差满足高斯分布的假设，表面高程误差采用非冰川区高程残差的均方根误差或方差进行评估。DEM数据格网点的高程值存在较强的空间自相关，因此，对样本数据选择时，必须去除空间自相关，本文采用156 m去相关距离，基于Bolch等［36］研究DEM数据的误差可用高程残差平均值和标准误差计算：

式中：N为像元数；SD与SE分别为标准差和标准误差；MED为高程残差平均值；E为误差（表2）。

表2 感兴趣区DEM校正前后误差Table 2 The error of DEMs before and after correction in interesting area

2.3.2 冰川面积和末端变化误差

基于RTK-GPS、TLS及UAV数据的冰川边界，获得1号冰川末端和面积变化，利用式（5）和式（6）进行冰川面积误差评估［37-39］：

式中：UT与UA分别为末端及面积不确定误差；λ为空间分辨率；ε为矫正误差，保证在一个像元内。计算得到，2012—2015年和2015—2018年面积误差分别为1.4×10-4km2和7×10-6km2。

2.3.3 大地测量法物质平衡估算

冰川物质平衡估算采用大地测量法，大地测量法物质平衡（Bgeo）是通过将冰川体积变化与平均密度相乘，然后除以面积得出的，该面积是两期面积的平均值。Huss等［40］认为冰川密度是变化的且小于冰的密度，建议使用（850±60）kg·m-3较为合适。因此，本文采用（850±60）kg·m-3作为冰川体积-物质平衡转换参数。

3 结果与分析

3.1 冰川表面高程及物质平衡变化

1号冰川表面高程在低海拔区域呈现明显冰量损失，而在高海拔区域呈轻微减薄或变厚趋势（图2）。考虑到2018年无人机航测仅获得1号冰川冰舌区表面高程，所以，2012—2018年间表面高程下降（1.34±0.88）m·a-1，对应物质平衡为（-1.13±0.18）m w.e.·a-1，累积物质亏损量（-6.78±1.08）m w.e.，东支冰量损失较西支明显，而2015—2018年间冰川表面高程下降（2.03±0.96）m·a-1，物质平衡（-1.72±0.19）m w.e.·a-1，累 积 物 质亏 损 量（-5.16±0.57）m w.e.。2012—2015年间冰川表面高程下降（0.83±0.57）m·a-1，物质平衡（-0.71±0.17）m w.e.·a-1，累积物质亏损量（-2.13±0.51）m w.e.（表3）。研究期间表面高程减小区域主要位于冰舌区及东支右侧，且东西两支表面高程减薄速率随海拔高度的增加而减小。2012—2015年表面高程在-6～2 m间变化，冰川上部呈积累状态［图2（a）］；2015—2018年冰舌区表面高程变化区间为-13～1 m，东支减薄速率大于西支［图2（b）］，所以，2012—2018年冰舌区表面高程整体在-20～1 m间变化，东支末端减薄速率大于西支，海拔高度越高表面高程减薄越小，在西支表面高程变化上更为显著［图2（c）］。

图2 乌鲁木齐河源1号冰川表面高程变化Fig.2 Changes of surface elevation of the Urumqi Glacier No.1

表3 乌鲁木齐河源1号冰川不同时期表面高程及物质平衡Table 3 Surface elevation and mass balance of the Urumqi Glacier No.1 in different periods

3.2 大地测量法与花杆/雪坑法物质平衡对比

Wang等［24］计算得到1981—2009年间大地测量法与花杆/雪坑法所得物质平衡的差异小于10%（表4）；随着观测仪器的发展和误差评估方法的改进，Xu等［25］基于Zemp等［12］提出的关于大地测量法和花杆/雪坑法间的物质平衡减小差异（δ）的计算方法，得到了1981—2015年1号冰川物质平衡的减小差异 值（δ=0.53）通 过95%的置 信区 间检 验（|δ|＜1.96），大地测量法和花杆/雪坑法分别得到物质平衡值具有一致性。基于此，本文利用上述方法分析2012—2018年间两种方法的差异性，其中，2012—2015年间1号冰川物质平衡的差异［（Bgeo-Bgla）/Bgla］为2.9%（表4），两种方法的减小差异值（δ=0.78）通过95%的置信区间，这与前人研究结果一致。

表4 乌鲁木齐河源1号冰川不同时期花杆/雪坑法与大地测量法累积物质平衡比较Table 4 Comparison of glaciological and geodetic cumulative mass balances of the Urumqi Glacier No.1 for given time periods

在整条冰川上利用大地测量法与花杆/雪坑法所得物质平衡的差异（δ）进行评估，而2018年4月无人机航测区域仅为1号冰川冰舌区域，为了对2012—2018年及2015—2018年间大地测量法与花杆/雪坑法获取的物质平衡进行对比，提取对应花杆点的物质平衡。结果显示，大地测量法和花杆/雪坑法的物质平衡两者相关系数（R2）分别为0.93和0.91，均方根误差（RMSE）为1.19 m w.e.和0.81 m w.e.（图3），表明两种方法得到的物质平衡能够较好的对应。

图3 乌鲁木齐河源1号冰川冰舌区单点大地测量法与花杆/雪坑法物质平衡对比Fig.3 Geodetic versus glaciological mass balances for the ice tongue of Urumqi Glacier No.1

结合1号冰川物质平衡长期观测结果发现，自1959年以来1号冰川年物质平衡和累积物质平衡呈减少趋势，物质损失明显［26］，基于线性统计方法按照斜率不同得到1号冰川物质平衡2个时期的变化特征：1985—1996年物质平衡为-0.27 m w.e.·a-1，物质损失显著大于1959—1984年（-0.08 m w.e.·a-1）；从1997年开始，物质亏损更为强烈，使得1997—2018年间年物质平衡为-0.68 m w.e.·a-1，其中，2010年冰川物质平衡值低至-1.33 m w.e［42］。总的来说，2012—2018年间1号冰川物质损失（-0.64 m w.e.·a-1）大于1980—2012年（-0.47 m w.e.·a-1），2012—2014年物质平衡亏损减缓，而2015—2018年间物质平衡亏损增加，表明近期1号冰川物质平衡亏损仍在持续。将1号冰川年物质平衡与全球41条参照冰川物质平衡结果相比，其平均值变化趋势相一致［42］，说明1号冰川物质平衡在一定程度上能代表全球山地冰川平均物质平衡变化。

3.3 冰川面积及末端变化

自1959年以来，基于野外观测发现1号冰川退缩强烈［23］。考虑到2018年仅得到1号冰川3 785～4 185 m间冰川区，同时认为冰川上部面积保持不变的条件下，利用2015年4 185～4 485 m间冰川区得到2018年完整冰川区用于估算总面积。2012年，1号冰川总面积为1.59 km2，到2015年减少为1.56 km2，到2018年减少至1.52 km2，2012—2018年1号冰川面积整体减少了0.07 km2，年均变化率为-0.01 km2·a-1。其 中，2012—2015年 和2015—2018年均变化率分别为-0.01 km2·a-1和-0.02 km2·a-1，后者面积退缩强度大于前者，总的来说，2012—2018年1号冰川整体上呈现出持续退缩趋势。

2012—2018年1号冰川末端变化退缩速率为6.28 m·a-1，东、西支退缩速率分别为7.64 m·a-1和4.93 m·a-1，其中，2012—2015年末端变化基本处于稳定状态，在2016年东西支变化趋势相反，2017—2018年退缩幅度显著，但整体上东支退缩速率显著大于西支（图4）。在相同气候背景下，末端变化差异不仅取决于物质平衡的变化，而且与冰川地形和热力学参数有关的冰川动力学过程有关。首先，2012—2018年东支表面高程减薄速率大于西支［图2（c）］，这与东西支末端退缩趋势一致，体现了物质平衡对末端退缩的影响；其次，由于西支冰舌区坡向朝向为东南、东支受山体遮掩及西支面积和横截面比东支小等因素的影响，在假设消融速率相同的情景下，导致2000年代之前西支退缩速率大于东支。此外，已知冰川运动速度的年际变化作用于末端位置，已有研究表明1号冰川2012年后东支运动速度显著小于西支，低于1号冰川年均运动速度［41］，导致西支末端海拔较东支越来越高，且与西支相比东支末端受山体遮掩影响变小。最终，除2017年外西支退缩速率小于东支退缩速率。因此，通过分析物质平衡、面积及末端变化表明乌鲁木齐河源1号冰川近期呈加速退缩趋势。

图4 1980—2018年乌鲁木齐河源1号冰川末端变化Fig.4 Terminus changes of the Urumqi Glacier No.1 during 1980—2018

4 结论

本文基于RTK-GPS、TLS和UAV等资料分析了2012—2018年乌鲁木齐河源1号冰川面积、末端和物质平衡变化，对比了近期乌鲁木齐河源1号冰川变化，结果表明：

（1）2012—2018年期间，乌鲁木齐河源1号冰川面积减少0.07 km2，年平均面积退缩率为-0.01 km2·a-1；乌鲁木齐河源1号冰川末端在2012—2018年间呈退缩趋势，变化率为6.28 m·a-1，东支退缩速率显著大于西支，与2012—2015年相比，2015—2018年末端变化退缩幅度更为显著。

（2）乌鲁木齐河源1号冰川物质平衡在2012—2018年间表面高程下降（1.34±0.88）m·a-1，物质平衡为（-1.13±0.18）m w.e.·a-1，物质损失明显，其中，2012—2015年间冰川表面高程下降（0.83±0.57）m·a-1，物质平衡为（-0.71±0.17）m w.e.·a-1；2015—2018年间冰川表面高程下 降（2.03±0.96）m·a-1，物质平衡为（-1.72±0.19）m w.e.·a-1，物质损失主要发生在消融区及东支右侧，2012—2018年间物质平衡变化率大于1980—2012年间物质平衡变化率表明乌鲁木齐河源1号冰川近期呈加速消融。

致谢：感谢西北研究院刘宇硕工程师、西北师范大学沈思民硕士研究生和兰州大学郑续硕士研究生在无人机等数据处理过程中的大力帮助！