0 引言

21世纪以来，全球气温上升已是不争的事实，其中受影响最大的就是北极地区.最近几十年来的监测数据表明，北极圈气温上升的速度几乎是同时期北半球平均升温速度的两倍(Bekryaev et al., 2010)，格陵兰冰盖消融的速度也逐渐加速.因此，对格陵兰冰盖的研究和保护刻不容缓.

ICESat(Ice, Cloud, and land Elevation Satellite)是世界上首颗绕地球飞行的激光测高卫星，其主要任务是确定两极冰盖质量的季节性和长期变化趋势.凭借激光地面光斑小、数据分辨率高等特点，ICESat测高任务在格陵兰冰盖变化的研究中取得了众多令人瞩目的成果.Slobbe等(2008)通过对ICESat重叠点的分析，得到2003至2007年格陵兰冰盖体积变化速度约为-149 km³·a^-1，其中高程2000 m以下区域冰盖表面正以-24 cm·a^-1的速度下降；Pritchard等(2009)采用三角内插的方法对ICESat数据进行处理，并分析了两极冰盖变化，结果表明格陵兰岛冰川高程下降速度与冰川流速有关；Sørensen等(2011)系统地研究了利用ICESat数据确定格陵兰冰盖质量变化的方法，得到格陵兰冰盖在2003至2008年间质量变化在-191~-240 Gt/a之间，但不同算法的结果存在一定差异；Schenk和Csathó(2012)则融合了ICESat及机载激光测高数据，并提出了一种称为交叉区域的新算法，得到格陵兰冰盖的体积变化速度为-259 km³·a^-1.国内也有不少学者利用ICESat数据对两极冰盖进行过相关研究，例如黄海兰(2011)对格陵兰冰盖的质量损失进行了估计，并与GRACE卫星重力观测结果进行了对比；史红岭等(2011)对南极冰盖的质量变化进行了计算，并分析了各冰川流域物质平衡变化尺度及空间分布特征；马跃等(2015)对ICESat重复点和交叉点数据确定冰盖变化的方法进行了讨论，并估计了格陵兰冰盖2000m以上区域的高程变化；李斐等(2016)在引入坡度改正后，利用ICESat重复轨道数据对南极冰盖的质量变化进行了分析.

在这些研究的基础上，本文针对坡度拟合函数的选择、平差法方程病态等问题，对ICESat重复轨道数据确定冰盖表面高程的算法进行了改进，并对2003年9月至2009年10月间格陵兰冰盖的高程和体积变化进行了估计.结果表明，改进后的算法可以有效地提高格陵兰冰盖边缘地区表面高程变化估值的空间分布率，从而更好地反映冰盖高程和体积变化.

1 数据介绍及粗差编辑 1.1 ICESat/GLAS(R34)数据介绍

ICESat卫星由美国国家航天局(NASA)于2003年1月在戈达德航天中心发射升空，轨道倾角为94°，纬度覆盖范围为±86°，飞行高度约600 km，搭载的主要观测仪器为地学激光测高系统GLAS(Geoscience Laser Altimeter System)，采样率40 Hz.ICESat/GLAS每年观测2~3次，每次连续观测约33天，轨道重复周期为91天，在6年半的工作时间内共传回18期观测数据.为使用方便，对其按观测顺序和使用的激光器号进行了编号，如表 1所示.其数据由美国冰雪数据中心(NSIDC)处理并分为15种产品进行发布，其中GLA12为Level 2两极冰盖测高数据，采用了专门针对冰盖表面的算法对脉冲波形进行处理，是本文的主要数据源.使用的数据版本为Release 34(R34)，是目前ICESat数据的最新版.与以往版本的ICESat数据相比，R34的最大变化是G-C(Gaussian-minus-centroid)误差的改正，更正了过去ICESat数据处理过程中波形拟合参考点错误的问题(Borsa et al., 2014).

1.2 粗差编辑与数据筛选

尽管NSIDC已经对GLA12数据进行了精密而富有针对性的处理，但由于观测条件复杂多变，数据集仍然存在一定数量的粗差观测值，因此，需要通过有效的数据编辑，剔除受系统误差或粗差影响的观测值.数据编辑的准则除了GLA12数据中提供的质量标记以外，主要通过对激光回波波形参数的筛选来实现，通过实验验证(陈国栋，2015)，本文所选的波形参数主要包括接收机增益、波形拟合标准差、回波反射率以及波形饱和改正，具体阈值如表 2所示.由于执行最后三个任务期时ICESat的激光器发射能量过低，接收机增益一直保持在较高的水平，导致其回波波形特征与前15个任务期有明显差别，因此，L2d—L2f三个任务期的编辑准则与其他任务期有显著区别.其他编辑准则包括：为避免普通陆地、海洋等无冰覆盖区域数据的干扰，根据格陵兰地区1 km分辨率地表类型格网数据(Bamber et al., 2013)进行筛选，仅保留冰盖及冰川表面的测高数据；为尽量避免地形坡度对计算的影响，剔除偏离ICESat设计轨迹500 m以上的观测值.经过筛选，每个任务期的格陵兰冰盖测高数据数量如表 1所示.其中，L1任务期的轨道为试运行轨道，与其余任务期所采用的正式轨道不重复，无法形成重复轨迹，故本文针对格陵兰冰盖变化的计算只采用L2a至L2f的17个任务期数据，时间跨度为6年.

2 确定冰盖表面高程变化的方法 2.1 基本方法

交叉点和重复轨迹是利用测高数据确定冰盖表面高程变化的两种常用方法，但在格陵兰地区，ICESat轨迹交叉点数量仅不到6000个，空间分辨率较低，特别是在轨迹分布较稀疏的南部地区，临近交叉点的间隔长达数十公里，难以体现冰盖变化细节.因此，本文采用重复轨迹的方法进行冰盖高程变化的估计.

由于卫星轨道具有周期性，因此，测高卫星在地球表面上的轨迹经过一段时间后会出现重复，称为重复轨迹，利用重复轨迹数据之间的高差和时间差，即可对地表高程变化进行估计.但由于卫星运行轨道、星下点指向等的变化，所谓的重复轨迹之间存在一定的偏差，并不能完全重合，因此，当地表存在坡度时，重复轨迹数据之间的高差除了地表随时间的变化外，还包含了由于重复轨迹位置差异所引起的高差.显然，为得到正确的高程变化速度，必须先将上述两种因素引起的高差进行有效的分离.然而，格陵兰地区现有的DEM分辨率和精度不足以对ICESat数据进行直接的坡度改正，因此，本文通过内插和参数估计的方法消除地面坡度的影响，具体做法如下：以ICESat设计轨迹为基准，将相邻的ICESat观测数据在沿轨方向上内插到与设计轨道垂直且通过最近的设计轨迹点的直线上，如图 1所示，其中黑点表示设计轨道的星下点位置，“+”表示卫星轨迹沿轨内插后的位置.对于每一排内插点(+号位置)，设内插点的高程为H₁，…，H_i，…，H_n，每个H_i可表示为

(1)

式中等号右侧第一项a为高程常数项，第二项为高程随时间的线性变化项，b和t_i分别表示高程年变化率及时间，三、四项是用三角函数表示的周年变化项，c、d为系数，H(D_i)则是用于拟合地面坡度的函数，自变量D_i表示内插点至设计轨迹的距离(规定位于设计轨迹东侧的距离为正，西侧为负).

2.2 坡度拟合函数选择

根据地表起伏的复杂程度，H(D_i)可以选择不同阶数的多项式进行拟合，显然，地面坡度变化越复杂，拟合函数的阶次就应更高.国外研究表明，格陵兰沿海地区高程起伏较大，因此，许多学者使用了较为复杂的曲面模型对格陵兰冰盖地面坡度进行拟合(Schenk and Csathó，2012；Ewert et al., 2012).然而，复杂的拟合模型会增加参数的数量，而激光难以穿透云雾的特点导致ICESat重复轨迹数据时常产生缺失，造成观测值数量不足而无法求解，从而使得局部区域数据分辨率下降，并引起对整个冰盖的高程和体积变化估计发生偏差.本文分别用一阶常数坡度和三次多项式(后文分别用d1和d3表示)对H(D_i)进行拟合，两者的差异如图 2所示.两种方法获得的冰盖体积变化约相差19 km³·a^-1(d1为-295.33 km³·a^-1，d3为-276.30 km³·a^-1)，差异主要集中在沿海区域.统计结果表明两者在高程2000 m以上的冰盖内陆仅相差1.06 km³·a^-1(d1—d3)，2000 m以下区域的差别则高达-20.09 km³·a^-1.进一步的分析发现，存在明显差异的区域都集中在三阶多项式拟合时由于参数过多而无法求解，造成数据空缺的区域，相比之下，采用一阶常数进行拟合时的结果则具有更好的空间分布.由上述结果可知，高阶的坡度拟合函数并没有对结果带来明显的改善，反而由于空间数据分布的缺失，在格网化时引起误差，得不偿失.因此，本文采用一阶函数对H(D_i)进行拟合，即认为图 1所示的同一排内插点之间的坡度变化一致，若以e表示坡度，则(1)式可表示为

(2)

2.3 方程病态问题解决

当重复轨迹在空间上的排列顺序恰好与观测时间顺序相同(或逆序)时(如图 3所示)，(2)式中的系数b和e有较强的相关性，从而造成方程中高程的时间变化趋势项bt_i和坡度变化项eD_i难以分离，甚至可能引起平差模型法方程病态而无法求解(Schenk and Csathó，2012).为解决这一问题，本文对算法进行了改进，加入虚拟观测值以改善平差模型结构，并通过权值调整尽量减小虚拟观测值对计算的影响.为此，首先利用ICESat交叉点数据获得了格陵兰冰盖高程变化格网(如图 4)，并将该格网结果作为虚拟观测值，即

(3)

其中b与(1)式相同，表示冰盖表面某点的高程年变化率，b′表示该处对应的虚拟观测值，即交叉点数据获得的高程变化格网值.联合上式与(2)式组成观测方程组.由于(3)式中参数b和e的系数分别为1和0，与(2)式组成系数阵后一定不相关，因此可以有效地改善法方程的病态性.同时，为了避免交叉点格网过多的影响解算结果，必须对虚拟观测值的权加以限制，通过实验(陈国栋，2015)，给出经验定权公式：

(4)

其中P_d表示虚拟观测值的权(其余观测值的权为1)，n为重复轨迹数量，N_p是本文定义的用于描述重复轨迹与观测时间之间相关程度的参数：对于观测时间为t_i的轨迹，判断所有在t_i时刻之后观测的重复轨迹与该轨迹的相对位置，规定在东侧为正、西侧为负，则对所有的n条轨迹共有N=n(n-1)/2次判断，其中N⁺个为正，N^-个为负，则N_p为

(5)

显然，N_p的取值范围为[0, 1]，并且N_p越接近1表示bt_i与eD_i相关程度越高.考虑到重复轨迹数量最多为17(ICESat正式轨道只有17个重复周期)，式(2)求解所需的最少观测量为5个，因此P_d的取值范围为[3/17, 1.6].当重复轨迹观测值数量较多或排列顺序较好时，P_d的取值都较小，因而可以保证在观测条件有利的情况下，虚拟观测值不会过多地影响平差结果.

用V₁、V₂分别表示加入虚拟观测值前后获得的高程变化，图 5所示为位于71.0°N，41.4°W的一段重复轨迹计算结果.从图中可以看到，V₁与N_p存在明显的相关性，在V₁发生跳变时，排列指数N_p的数值都明显偏高，而在加入虚拟观测值后，V₂的结果则比较平滑，没有明显的突变.由于该轨迹位于冰盖内陆地区，高程变化量级小且变化稳定(Zwally et al., 2012)，不可能出现V₁所显示的这种剧烈变化，因此可以判断加入虚拟观测值显著改善了估算结果.而对于较小的N_p值，V₁、V₂的差值则非常小，表明本文提出的经验定权方法也是有效的.

通过上述方法解决法方程病态问题后，可以有效地增加高程变化估值的数量，从而能够更好地反映格陵兰冰盖体积变化的情况.图 6所示为加入虚拟观测前后高程变化估值差异在1 m·a^-1以上数据的统计图及位置分布.尽管这些数据只占数据总量的1.2%，但从图 6的结果可以看出，其分布主要集中在冰盖边缘，这些区域由于坡度较大导致ICESat观测值数量较少，而冰盖边缘的冰川流失又是格陵兰冰盖体积变化的重要因素，因此，在加入虚拟观测值后得到改善的这些数据对于提高冰盖边缘数据密度、正确反应冰盖高程和体积变化具有重要意义.

实际计算冰盖高程变化时，只在N_p＞0.2时采用虚拟观测值参与计算，以尽量减少虚拟观测值对最终结果的影响.

3 结果 3.1 ICB改正的探讨

由于模型误差、卫星轨道变化、测高仪工作状态等众多已知和未知的因素，卫星测高数据本身可能存在一定的时变误差，ICESat各任务期的数据间也存在这种系统性的差异，称为任务期偏差(Inter-Campaign Bias，ICB)(Siegfried et al., 2011)，在对冰盖变化进行研究时必须予以考虑(Hofton et al., 2013；李斐等，2016).本文采用纬度±66°之间Release34的GLA15数据(ICESat海洋观测数据)对ICB进行了分析，首先以3.0 mm·a^-1的速率对观测数据中包含的海平面变化进行改正，然后通过交叉点不符值确定各任务期测高值之间的偏差.所得结果与国外相关研究获得的ICB变化进行对比，如图 7所示，其中，Shepherd等(2012)是根据±66°纬度之间ICESat测高数据与海平面模型的差值确定ICB的变化(Shepherd et al., 2012)，Ewert等(2012)是对南极沃斯托克湖表面测高数据分析获得的结果(Ewert et al., 2012)，Hofton等(2013)则通过对东南极地区数据的分析获得(Hofton et al., 2013).由于现有研究成果使用的均是未加入G-C误差改正的旧版ICESat数据，为增加与本文结果的可比性，对上述国外研究成果均进行了G-C误差变化趋势的改正(Borsa et al., 2014).由图 7可知，Shepherd等(2012)与本文所采用的数据分布范围一致，ICB分析结果也非常接近，而采用的数据分布区域不同时，Ewert等(2012)、Hofton等(2013)与前两组结果都有明显的差异，因此，本文认为，ICB分析结果与数据分布有关，上述几组结果并不能直接用于格陵兰地区的数据改正.另外，根据本文的分析结果拟合的ICB线性变化趋势不到0.1 mm·a^-1，表明R34版本数据中，由ICB造成的长期性变化趋势已基本消除，对高程时变计算的影响较小，这与Borsa等(2014)的结论一致.因此，本文没有对ICB进行改正.

3.2 地壳形变改正

根据作用原理的不同，地壳形变可以分为两类，弹性回弹(Elastic Rebound，ER)和冰后回弹(Postglacial Rebound，PGR)(Spada et al., 2012)，两者都不会引起冰盖体积的变化，需要在ICESat获得的冰盖表面高程变化中进行剔除.本文利用ICE-6G冰川均衡模型(Peltier et al., 2015)对PGR进行了改正，并利用Khan等(2010)提供的格陵兰岛弹性回弹模型进行ER改正.加入两项改正后，对格陵兰冰盖体积变化估值的影响合计约-6.04 km³·a^-1.

3.3 格陵兰冰盖体积变化

利用ICESat重复轨迹数据，采用第3节所述方法计算2003年10月至2009年10月间格陵兰冰盖表面高程变化，并按照5 km×5 km分辨率进行格网化，结果如图 8所示.参考国外相关研究(Sørensen et al., 2011)，采用自助抽样法对精度进行估计后，得到整个格陵兰地区冰盖、冰川的平均高程变化约为-16.79±0.84 cm·a^-1.将每个格网的高程变化乘以格网面积后得到整个冰盖的体积变化为-301.37±15.16 km³·a^-1，其中呈现高程上升趋势的面积约为866×10³ km²，占总面积的48.2%，主要集中在海拔较高的冰盖内陆地区；而体积损失则主要发生在冰盖边缘，特别是东北沿海、西南沿海以及雅各布港等冰川密集区域.

表 3所示为近年来不同学者利用ICESat数据获得的格陵兰冰盖体积变化与本文结果的对比(为与国外研究统一，表 3中结果均未进行PGR和ER改正)，尽管这些研究所采用的ICESat数据版本中尚不包含G-C误差改正，但Ewert等(2012)和Schenk & Csathó(2012)已对结果进行了ICB偏差的改正，可以在很大程度上消除G-C误差的影响(Borsa et al., 2014)，而未进行任何改正的Slobbe等(2008)和Sørensen等(2011)的结果则通过Borsa等(2014)提供的G-C误差变化趋势进行相应改正，以增强与本文结果的可比性.

表 3的结果表明，在以ICESat单一数据源进行的研究中，本文获得的格陵兰冰盖整体消融速度明显高于既往研究成果.虽然本文采用的观测数据时间跨度与其他研究略有不同，但当采用2003年10月至2008年3月间的数据并按照本文算法进行计算时，其结果与表 3所给出的本文结果之间的差异不到10 km³·a^-1，与Schenk和Csathó(2012)的研究结论一致，因此，时间差异并不是本文结果与国外结果产生差异的主要原因.

除ICESat数据以外，Schenk和Csathó(2012)的研究还采用了ATM机载测高数据，局部地区的数据空间分辨率更高，比表 3中的其他结果更具可靠性，而该结果与本文结果较为接近，说明本文的结果是准确的.对于冰盖局部区域结果的对比则表明本文并没有过高地估计格陵兰消融速度.例如，在高程2200 m以上的冰盖内陆区域(占格陵兰冰盖总面50%以上)，本文获得高程变化速度为+0.8 cm·a^-1，而Zwally等(2011)的结果(ICESat单一数据源)为+1.0 cm·a^-1，对应的体积变化差异仅为-2 km³·a^-1左右.在冰盖边缘区域，本文结果则与融合了多源数据的观测结果较为接近，例如Howat等(2005)联合ICESat测高和ASTER遥感影像数据得到的格陵兰冰盖东南沿海2000 m高程以下区域冰川体积变化为-80.4 km³·a^-1，Hurkmans等(2012)联合ICESat卫星测高数据、ATM机载测高数据对雅各布港冰川的研究表明该区域冰川流失速度为-20.5 km³·a^-1，Kjeldsen等(2013)同样利用ICESat和ATM数据对西北沿海冰川进行了研究，得到的体积变化为-55.8 km³·a^-1，而本文在这些地区的结果分别为-75.6、-20.3以及-51.9 km³·a^-1.

上述对比表明本文对格陵兰冰盖体积变化的估计是可靠的，而造成表 3中结果差异的原因主要是计算方法的不同，本文对算法的改进以及低阶坡度拟合函数的应用都提高了数据的空间分布率，特别是在冰盖边缘地区，数据量的增加能够更好地反映冰川流失情况，从而更准确地估计冰盖整体的体积变化速度.

4 冰川流域体积和高程变化分析

为更好地对格陵兰冰盖进行研究和分析，Zwally等(2012)根据冰川流向将格陵兰冰盖划分为8个冰川流域系统，每个流域又分割为若干个子流域(如图 8中灰色实线所示).在流域划分的基础上，本文以2000 m等高线(图 8中黑色虚线)为界进一步细分，分别计算了每个区域内冰盖的体积变化情况，如表 4所示.需要注意的是，表 4中结果只含格陵兰冰盖主体，与冰盖不相连的独立冰川、冰帽等都未统计在内，因此表 4中总和与3.3节的结果相差了约60 km³·a^-1，这几乎是格陵兰冰盖体积流失总量的四分之一，可见格陵兰岛周边独立冰盖和冰帽的流失情况也非常严重，甚至有研究认为这些独立冰川与冰帽的消融速度达到了格陵兰冰盖本体的2.5倍(Bolch et al., 2013).

根据表 4结果可知，从各流域系统来看，仅有DS2流域有体积增长的趋势，而DS3、DS4、DS7和DS8都呈现较严重的下降趋势.冰盖体积流失主要发生在2000 m以下的沿海区域，其中DS1.2-1.4以及DS2.2处在大体平衡的状况，其余区域由于冰川流失严重，则都呈现明显的体积减少趋势.国外早期研究表明，由于雪层的堆积，格陵兰冰盖2000 m以上的区域处于明显的体积积聚状态，在2007年之前的高程上升速度约为2~6 cm·a^-1(Johannessen et al., 2005；Zwally et al., 2005；Thomas et al., 2006；Slobbe et al., 2008).但表 4结果则显示该区域在2003—2009年间的体积变化仅为+1.90 km³·a^-1(平均高程变化约为+0.2 cm·a^-1)，相较于此前的结果，其增长趋势明显下降.

为更好地分析格陵兰冰盖高程/体积变化情况，以表 4中的每块小区域为单位，用三次函数拟合其平均高程变化情况，即

(6)

式中的a和e由(2)式的平差结果获得，b_i为高程随时间变化的各阶系数.根据(6)式得到各区域高程变化趋势及平均速度如图 9所示.

图 9结果表明，冰雪体积损失最严重的流域是DS4、DS5以及DS8，速度均超过-20 cm·a^-1，这些区域冰川流失严重，因此，即使是2000 m高程以上的区域也出现了不同程度的下降趋势，特别是DS5和DS8的各个子流域，从2007年开始，其高程下降速度(2000 m以下)有明显加速趋势.DS3流域平均高程下降速度虽不如上述几个流域严重，但其中的子流域DS33却是所有子流域中累积高程下降最多的，其2000 m以下的区域在6年间下降了近10 m，且下降速度仍在加剧，其中的重要原因是位于该子流域的Kangerd冰川的流失(Howat et al., 2005).类似的情况还有DS7流域中的子流域DS71，全球流速最快的雅各布港冰川位于此处(Hurkmans et al., 2012)，但不同与DS33的是，该区域的速度变化较为稳定，没有明显的加速趋势.

由于距离海岸较远，DS1、DS2、DS6流域以及DS3的DS31、DS32两个子流域体积损失并不严重甚至略有增加，多数区域2000 m以上都处在高程增长的状态，特别是位于气候更寒冷的格陵兰北部的DS22子流域，2000 m以下也处在高程增长的状态，这在格陵兰非常少见.但这些流域的变化情况并不能令人安心，其中的大部分区域在2007年之后都逐渐显现出下降趋势，例如DS1流域2000 m以上区域以及DS22子流域2000 m以下区域都由高程增加转为高程下降，DS12、DS13、DS14子流域2000 m以下区域下降速度明显加速，而DS21、DS32子流域以及DS6流域2000 m以下区域则是一直保持着下降趋势，速度较为稳定.

5 总结

本文基于ICESat卫星激光测高数据，在对算法进行改进的基础上，采用重复轨迹法对格陵兰冰盖在2003年9月至2009年10月间的高程和体积变化进行了分析.通过与国外研究结果的对比发现，由于在坡度拟合和平差模型方面的改进，本文的算法可以提供更多的高程变化估值，从而提高结果的空间分辨率，减小格网化时的内插误差，更好地体现冰盖边缘体积变化.根据本文的结果，2003—2009年间格陵兰冰盖年平均体积变化约为-301.37±15.16 km³·a^-1，体积损失主要发生在冰盖边缘，特别是西北沿海、东南沿海以及雅各布港等区域.而对格陵兰冰盖各流域高程变化的分析则表明，冰盖内陆地区高程2000 m以上的区域高程积聚的速度已经明显减缓，而且许多区域的高程变化趋势在2007年之后都发生了改变，表明格陵兰冰盖消融的情况仍在加剧.

致谢  感谢NSIDC提供ICESat数据产品、布里斯托大学Bamber教授提供格陵兰地区1 km分辨率地表类型格网数据、多伦多大学Peltier教授提供的ICE-6G冰后回弹模型、丹麦技术大学Khan教授提供的格陵兰岛弹性回弹模型以及美国国家航空航天局Zwally教授提供格陵兰冰盖冰川流域系统划分.感谢各位评审专家和编辑对本文的肯定以及提出的宝贵建议，为研究内容的完善提供了许多有益的帮助.