中亚高山冰川表面高程变化时序重建

都伟冰， 张世琼， 李均力， 包安明， 王双亭， 史宁可，许琳娟， 高 鑫， 马丹丹， 郑岩超

（1.河南理工大学测绘与国土信息工程学院，河南焦作 454003；2.中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室，新疆乌鲁木齐 830011；3.新疆遥感与地理信息系统应用重点实验室，新疆乌鲁木齐 830011；4.黄河水利委员会黄河水利科学研究院，河南郑州 450003）

在全球变暖背景下，气候变化加剧，中亚冰川变化表现出明显区域差异。青藏高原中部和西北部冰川相对稳定，而高原周边山区的冰川物质亏损严重［1］。喀喇昆仑地区冰川异常，冰川跃动现象延伸到昆仑西部和帕米尔，与中亚其他地区的负质量平衡形成鲜明对比［2］。目前，喀喇昆仑冰川质量增加暂缓，且开始消融［3］。冰川变化形式复杂多样，冰川积累或消融必须考虑表面高程变化因素。了解当前和未来冰川质量的变化，以此预测干旱区冰川融水径流的未来变化情况和自然灾害的发生，提前提出防治措施，对于淡水资源安全与经济可持续发展等有重要意义［4］。

目前，测定冰川变化主要有花杆测量法、卫星定位监测方法、激光测高等［5］。冰川区坡度大、冰碛覆盖和冰前湖发育特征使冰川变化的局部复杂化，增加了冰川高程变化监测的难度［6］。随着遥感技术的兴起，由于其多时相大范围覆盖的特点，被广泛应用于测定冰川高程变化。

本文采用2003—2009 年的ICESat 卫星激光测高数据，构建适合中亚高山区冰川表面高程数据的处理模型，获取中亚地区2003年各脚点每年的拟合高程，最后根据拟合高程结合降水、温度以及地形等因素分析研究区冰川表面高程变化特征及驱动因素。

1 原理与方法

1.1 点云去噪及其精度优化方法

ICESat卫星在中低纬度地区没有精密航天器指向控制，造成单一重复轨道在中低纬度地区并不完全匹配，在轨道上相隔几百米甚至几公里［7］。本文通过建立去噪算法和多项式模型等提高ICESat 激光测高数据在空间和时间上的一致性。

本文以SRTM 高程数据作为基准面，采用标准差法剔除由于受云层等因素影响的ICESat 高程异常值，从而保证实验数据的可靠性。所有激光脚点与基准DEM的差值服从正态分布时，设激光光束脚点与基准DEM 的高程差值为D=[d1,d2,d3,···，dN-1,dN]，其平均值为-d，标准差为e，di表示第i个高程脚点，N为激光脚点的总个数。则激光脚点与基准DEM的高程差值落在(-3e,+3e)的概率P如下式所示。

如果di服从正态分布，在3倍标准差的方法中，异常值就可定义为与平均值偏差超过3倍标准差的值；如果di不服从正态分布，可以采用远离平均值的多倍标准差来描述。后续在2.3.2节中，具体说明如何根据正态分布结果，采用远离1倍标准差，来剔除异常值。

1.2 大范围冰川高程变化时序重建

以SRTM高程为自变量，采用公式（1）剔除置信度低的高程脚点，并建立区域高程拟合多项式模型，确定每年激光雷达卫星ICESat 脚点高程与SRTM高程之间的关系。对每组给定的数据拟合点(xi,yi)(其中，i≤m，m为脚点数量；xi为中脚点i对应的重采样SRTM高程；yi为ICESat激光光束脚点i的高程)，为了获取xi与yi的n（n≤m）次多项式函数关系Pn(xi)。

式中：k为第k（k≤n）多项式；ak为拟合多项式的系数。当拟合值Pn(xi)与yi越接近，表示激光脚点i获取的高程值越精确。当所有激光脚点的拟合高程Pn(xi)与ICESat 脚点高程yi差值的平方和I值最小时，由系数ak构成的多项式模型也就越精确。

当拟合函数为n次多项式时，满足式（3）的Pn(xi)称为最小二乘拟合多项式；当n=1时称为线性拟合，I为系数a0,a1,···，ak(k=0,1,···，n)的多元函数。上述问题即为求I=I(a0,a1,···，ak)的极值问题。对公式（3）以aj为未知数求导数为0 时(j=0,1,···，n)，得多元函数求极值的必要条件公式如下：

将公式（4）按照关于aj(j=0,1,···，n)展开，用矩阵表示为式（5）的法方程组。

公式（5）的系数矩阵是对称正定矩阵，存在唯一解。当n=3时，求解得到的拟合高程Pn(xi)与基准高程xi相关系数R2≥0.95，选n=3 得到区域的拟合高程，即为激光测高某个时间激光脚点i的拟合高程值与基准DEM的拟合多项式如下：

式中：xi为基准高程值；Pn(xi)为待校正高程脚点数据；a3，a2，a1，a0分别为拟合多项式的系数。

2 实验

2.1 研究区概况

中亚山区位于亚洲中部海拔在3500 m以上，冰川面积约1.3×105km2。为研究中亚高山冰川区内部的冰川变化差异，本文根据地貌、地形等生态地理区域系统将研究区划分为3个区域（图1）。I区主要包括西藏及青海南部地区，该区东南部冰川大部分属于海洋性冰川，内部冰川大多属于大陆性冰川；II区主要包括新疆地区及青海北部和甘肃北部的部分区域，该区冰川以大陆性冰川为主，但天山、昆仑山东段以及喀喇昆仑山北坡等部分地区冰川属于亚大陆性冰川；III区在中国境外，与中国相邻，地势东高西低，大部分为温带大陆性气候。

图1 中亚冰川分区概况示意图Fig.1 Overview of glacier zoning in Central Asia

2.2 数据与预处理

本文使用数据主要包含以下3部分：（1）ICESat数据来源于冰、云和陆地测高卫星，卫星于2003年1月13 日成功发射，是首颗激光测高卫星，搭载地学激光测高系统GLAS［8］；ICESat数据产品在美国国家冰雪数据中心(http://www.nsidc.org)免费获取。根据研究目的，本文采用2003—2009 年的ICESat GLAH14数据，为便于时序重建，避免降雪和融化等季节间冰川高程变化的影响，选择每年的同一月份数据［9］；根据统计的数据量，选择数据量最多且最完整的每年3 月的数据进行处理分析。（2）本文选择能够免费获取的90 m 分辨率包含研究区范围的SRTM3数据作为参考DEM。（3）冰川编目数据采用2017年7月发布的RGI 6.0(Randolph Glacier Inventory 6.0)数据库，RGI 6.0 以2010 年左右的夏季无云Landsat TM/ETM+为主要影像来源，用于获取研究区冰川边界。

2.3 数据处理

首先从ICESat 数据中提取脚点的坐标、高程，进行坐标转换和地理配准；再剔除受地形、云层等因素影响产生的异常值，本文根据正态分布检验结果选定异常值并将其剔除。利用去除异常值后的ICESat 坐标数据，在SRTM DEM 数据中重采样，获取ICESat 脚点坐标对应的SRTM DEM 高程，建立每年的多项式拟合函数，以解决ICESat 数据重复轨道不完全匹配的问题，并计算相关系数R2，作为评估拟合效果参数。以2003 年的SRTM DEM 高程数据作为基准，获取2003年每个脚点数据每年的拟合高程，实现时序重建。

2.3.1 数据预处理 由于本文始终以SRTM DEM 为基准进行时序重建，因此需将Topex/Poseidon 高程转换为WGS84高程，如下式。

式中：H84为转换后点云WGS84高程；hglas为ICESat点云相对于Topex/Poseidon 参考椭球的高程（m）；N是大地水准面和参考椭球面的距离（m）；数值0.7的单位为m，为2 个参考椭球间的高程差异［10］。本文在30 km×30 km 尺度范围内选取最高点作为SRTM DEM 山顶点，对比参考椭球为WGS84 高分辨率四维地球影像筛选均匀分布可靠的山顶点，记录其作为控制点的150 个山顶点的经纬度坐标，对SRTM DEM 进行配准，提高不同坡度、坡向的ICESat 脚点精度［11］。

2.3.2 剔除异常值 根据获取的冰川边界筛选出研究区内ICESat 数据高程点。以Ⅰ区的2003 年数据为例，检验ICESat 脚点高程是否符合正态分布。如表1所示，显著性为0.00＜0.05，所以该区域高程值不服从正态分布，因此，选定远离平均值1倍标准差的值为异常值并将其剔除，确保数据的准确性。

表1 I 区的2003年ICESat脚点高程正态分布检验Tab.1 Normal distribution test of ICESat foot elevation in area I in 2003

2.3.3 多项式拟合函数 以SRTM DEM高程为自变量，与ICESat 高程进行多项式拟合，建立多项式拟合函数如式（6），并计算相关系数R2。如图2 所示，通过对比Ⅰ区内2003 年点云数据去除异常值前后拟合效果，发现去除异常值的R2更接近1，拟合函数的拟合效果更好，更符合冰川表面高程的实际情况。

图2 Ⅰ区2003年数据去除异常值前后对比Fig.2 Comparison of data in zone I before and after removing outliers in 2003

3 结果与分析

3.1 多项式拟合模型评估

由于卫星数据的处理过程带有保密性，难以给予物理模型对数据源进行校正。本文通过统计学的方法，采用确定系数R2评估多项式拟合效果，表示模型对现实数据拟合的程度，R2的计算公式如下。

式中：SSR 为预测数据与原始数据均值之差的平方和；SST 即原始数据和均值之差的平方和；SSE 为拟合数据和原始数据对应点的误差的平方和；R2的取值范围为［0,1］。当拟合数据和原始数据对应点的误差越小时，SSE越小，拟合数据越可靠，此时R2的取值越接近1。因此，当R2的取值越接近1，拟合效果越好。

由表2 所示，整体冰川区、I 区、II 区、III 区的R2都大于0.98，说明整体拟合效果很好。冰川区由物质平衡线分为消融区和积累区，在物质平衡线以下为消融区，物质平衡线以上为积累区，冰川表面高程变化可总体归为这2 种模式，因此利用三次多项式更能符合2 种模式的冰川表面高程变化现象，且符合冰川表面高程自身变化规律。

表2 研究区2003—2009年3月多项式拟合模型的R2Tab.2 R2of the polynomial fitting model from 2003 to March 2009 in the study area

3.2 冰川表面高程变化分析

在数据处理的过程中由于系统误差会产生部分异常值，本研究采用式（9），根据正态性检验结果，将数据中与变化平均值的偏差超过2 倍标准差的高程变化平均值作为异常值进行剔除。各区不同时间段内冰川表面高程平均变化量见表3。

表3 各区不同时间段内冰川表面高程平均变化量Tab.3 Average variation of glacier surface elevation in different periods of time in each region

式中：E为误差=高程变化值-高程变化均值。

由图3a可知，整体区域即整个中亚地区所有的冰川区域，总的冰川表面高程下降9.59±1.89 m，在2004年高程下降了5.89±3.76 m，次年回升到与2003年相当的高度，之后逐年缓慢下降；图3b中，Ⅰ区冰川表面高程在2003—2004 年间下降5.78±2.13 m，2005 年略有回升，之后高程逐年缓慢下降，变化斜率为-0.581 m·a-1；图3c 中，II 区其冰川表面平均高程在2003—2009 年间下降了7.87±5.03 m，2005 年冰川表面高程出现大幅上升，甚至超过2003年的高程值，但整体仍处于下降趋势，变化斜率为-1.226 m·a-1；图3d所示，与其他分区相比，Ⅲ区每年的冰川表面平均高程变化波动较大。

图3 各个研究区冰川表面高程变化时序重建Fig.3 Time series reconstruction of glacier surface elevation change in each study area

如图4所示，整个研究区内天山、帕米尔高原和青藏高原东南部地区的冰川高程下降剧烈，青藏高原南部地区冰川高程呈轻微下降趋势。中部的昆仑山西段地区部分冰川高程呈现升高状态，西段以及昆仑山东部的祁连山地区都呈下降趋势。北部的天山地区冰川高程变化呈现出明显的区域差异，东天山、北天山和西天山南部地区冰川高程下降明显，西天山北部冰川高程下降速率相对较慢。整个研究区内除昆仑山部分地区冰川高程增厚外，其他地区整体呈下降趋势，只是下降速率存在地区差异。与2019 年Treichler 等［2］的研究相比，各个区域的冰川表面高程变化趋势基本符合。但由于本文只采用了每年3 月的数据，数据量小，分布稀疏，部分区域缺少数据支持且表面高程变化特征表现性差，例如在Treichler等［2］的研究中，中天山东南部地区、昆仑山大部分地区以及青藏高原中部的部分地区冰川呈轻微增厚状态，但本文中在该部分地区数据稀少，无法确认冰川变化情况。且本文中颜色栅格表示小区域范围内激光点的变化率，因此，与Treichler 等［2］研究中大分区冰川表面高程变化率存在部分差异。

图4 各区冰川表面变化速率分布Fig.4 Distribution map of glacier surface change rate in each region

根据2019年Treichler等［2］对MERRA-2和ERAInterim气候再分析数据的处理，显示在2003—2009年间，气温呈上升趋势，除昆仑山以及周围少数地区外，亚洲冰川表面高程整体呈下降趋势。青藏高原夏季降水量，整体呈现自东南向西北递减的分布规律［12］，东北部干旱区雨季降水量也呈增加趋势［13］。图3 中，在2004—2005 年间各区的冰川表面高程变化剧烈。据相关研究［14］显示，2003 年12 月至2004年2月，研究区冬、春季大部分地区气温较常年同期明显偏高，降水偏少，造成冰川消融。2004 年新疆南部等地夏季高温日数普遍达到10～20 d，局部地区达到40 d［15］。2005 年西部地区伴有降水降雪的增加，造成各区的冰川表面高程在2005年都呈上升趋势［16］。其中，新疆东南部和西部以及青藏高原北部等地较常年同期偏高3倍［17］。

地形上，冰川在低海拔坡度小于30°时，冰川厚度减薄随着坡度的减小而加剧，冰川表面高程降低较快；中高海拔地区，温度逐渐降低，受降水的影响较大，易导致冰川表面高程增高；而在高海拔地区平均坡度＞35°，易发生冰崩、雪崩等现象，使得冰川物质减少［18］。结合地形因素，降水随海拔增高而增多，但有一个最大降水量高度，超过此高度，山地降水不再随高度递增，随气候干湿而异［19］。I区中，青藏高原平均海拔4000 m以上，海拔高度对气温的影响已超过纬度位置作用。因此，青藏高原除南部和东南部的冰川退缩较严重外，大部分区域冰川高程变化相对缓慢［20］。II 区南部有昆仑山和阿尔金山，这2个地区海拔高，干旱少雨，冰川高程变化相对稳定，但天山西北部受降水量影响较大的外山脉冰川融化速度较快。祁连山地表平均温度在海拔高的地区变暖趋势明显，因此，祁连山地区冰川表面高程缓慢下降［21］。Ⅲ区中西帕米尔高原山脉交错复杂，降水丰富，该地区冰川高程变化强烈。东帕米尔高原东南边缘高山阻隔了来自印度洋、太平洋的暖湿气流，气候干燥，冰川高程缓慢下降［22］。就下降趋势而言，研究区内境外区域冰川表面高程下降趋势最快，为-1.663 m·a-1，新疆地区次之，青藏高原地区最为缓慢。

4 结论

本文针对ICESat 卫星测高数据，考虑到重复轨道在中低纬度地区空间匹配性弱的问题，选取SRTM DEM 为高程基准面，建立冰川区点云去噪及其精度优化算法和多尺度冰川区表面高程时间序列分析模型，为点云类遥感监测冰川表面高程变化提供参考。主要结论如下:

（1）采用正态分布显著性检验，进行粗差剔除，在不服从正态分布情况下，选定远离线性回归1 倍标准差的高程值为异常值，能够有效剔除冰川区ICESat脚点高程中的异常值。

（2）利用三次多项式模型，对中亚山区整体和分区域的不同尺度冰川表面高程进行拟合，模型R2都在0.98 以上，表明ICESat 数据和SRTM 高程数据的三次多项式关系在该区具有较强的普适性。

（3）亚洲高山区的冰川表面高程2003—2009年整体下降幅度在10 m以内，表现出明显的区域差异，对于2004 年前后表现出强烈变化，由于环流和季风等因素的影响，造成的温度以及降水异常，从而导致冰川异常波动。

本研究方法对冰川区点云类高程脚点监测具有一定的通用性，但会使点云数据更稀疏，另对基准DEM的依赖度较高。后续，研究将进一步加入其他冰川区点云数据，研究时间序列建立所适用的空间尺度，建立更灵活多变的区域冰川表面高程时序监测模型，对影响冰川表面高程变化的地形、环境等因子进行更有效全面的分析。