南极冰盖作为地球上最大的冰体,如果全部融化,海平面将上升约60 m,这无疑会给人类带来巨大的环境灾难,因此精确评估南极冰盖质量平衡及其相关的气候响应具有显著的科学、社会和经济效益(Velicogna and Wahr, 2006).但由于南极冰盖地处偏远且地域广袤,探测其质量平衡实际上非常困难,近20年来随着卫星测高、卫星重力等空间大地测量技术的迅速发展,为南极冰盖质量平衡的监测提供了多种强有力的新型观测手段.空间大地测量技术具有受气候环境限制小、覆盖面积大、监测能力强、可连续大规模对冰盖进行监测并获取多源数据的技术优势,同时也是一种高效经济的观测手段和方法.其中2002年发射的重力卫星GRACE(Gravity Recovery and Climate Experiment)可通过测量地球重力场随时间的变化来获取南极冰盖质量的时空变化,GRACE探测南极冰盖质量平衡已成为近年来冰盖物质平衡研究领域的热点和重点,并涌现出一大批重大的研究成果(如：Tapley et al., 2004; Velicogna and Wahr, 2006; Jacob et al., 2012; King et al., 2012; Gardner et al., 2013; Shepherd et al., 2012).但这些研究结果彼此并不一致,甚至差异很大,最新的研究结果表明差异较大的原因在于南极GIA(Glacial Isostatic Adjustment)模型的不确定性太大(Shepherd et al., 2012; Velicogna and Wahr, 2013; 高春春等,2015).

GIA主要是指黏弹地球对末次最盛冰期以来冰退期地表冰负荷和海水负荷变化的响应,GIA能够引起地球内部物质的重新分布导致地球重力场的变化(Wang and Wu, 2006; 贾路路等,2011),因此为了估算冰盖质量变化,需要从GRACE重力变化中扣除GIA信号.但由于目前对冰盖历史变化和地幔黏滞度结构缺乏足够的认识,现有的GIA模型不可避免地存在着很大误差和不确定性(贾路路,2012).前期的研究结果表明,不同的GIA改正导致GRACE估算结果的差异可达80 Gt/a,量级甚至超过了冰盖质量变化本身(高春春等,2015),因此GIA模型本身的极大不确定性很大程度上降低了GRACE结果的可靠性,南极GIA信号的精确估算成为GRACE监测南极冰盖物质平衡研究中的重点和难点.很多研究已经提出采用其他观测数据对GIA信号进行约束和分离,如GPS(Global Positioning System)数据(Wang et al., 2013; Wu et al., 2010),ICESat(Ice,Cloud,and land Elevation Satellite)数据(Riva et al., 2009; Gunter et al., 2014).就南极地区而言,GPS观测站布设稀少很难形成密集的观测网,难以单独准确构建GIA,而2003年发射的激光测高卫星ICESat可直接观测南极冰盖表面高程变化,与 GRACE卫星有8年的重叠观测期(2003—2009年),因此联合ICESat和GRACE数据分离南极GIA信号是目前最佳的选择.

联合GRACE和ICESat数据分离南极GIA信号的想法最早由Wahr等(2000)提出,他们利用5年的GRACE和ICESat模拟数据对GIA进行了分离,指出两种卫星数据的联合能够有效减小GIA结果的不确定性(Wahr et al., 2000).随后Velicogna和Wahr(2002)提出利用卫星重力、卫星测高和GPS数据在频域中通过迭代计算实现南极 GIA和现今冰盖质量变化的分离,并做了数值模拟(Velicogna and Wahr, 2002).GRACE和ICESat卫星发射后,Riva等(2009)利用5年的GRACE和ICESat数据,采用粒雪和基岩的密度差结合质量守恒和体积守恒原理在空间域实现了南极 GIA和冰盖质量变化的分离,得到南极GIA信号在南极引起的质量变化为100±67 Gt/a(Riva et al., 2009).Gunter等(2014)利用2003—2009年的GRACE和ICESat数据,对Riva等的分离方法进行了改进,采用不同数据中心公布的GRACE数据得到南极GIA信号引起的质量变化为53～103 Gt/a(Gunter et al., 2014).

本文选取2003—2009年的GRACE和ICESat数据,并结合同一时期的南极粒雪密实化FDM(Firn Densification Model)模型和表面质量平衡SMB(Surface Mass Balance)模型对南极GIA信号进行了分离,在分离过程中提出冰流速度加权法和GPS球谐拟合改正法对分离的GIA信号进行进一步精化,同时对分离的GIA信号进行详细的精度评估和GPS验证,最后对南极GIA信号和冰盖质量变化及其时空分布特征进行了讨论和分析.

ICESat卫星于2003年1月13日在美国加利福尼亚范登堡空军基地成功发射,是人类历史上第一颗激光测高卫星,该卫星的主要任务之一就是测量冰盖的高程变化,卫星轨道倾角94°,覆盖南极大部分地区,卫星携带的激光测高仪GLAS(Geoscience Laser Altimeter System)对冰盖表面反射非常敏感,反射信号不会深入冰雪面,激光脚点小(约60 m),沿轨采样密集(水平间距约170 m),测量精度可达±10 cm(Schutz et al., 2005),ICESat激光测高仪的工作任务期见表 1.利用GLAS数据采用交叉点分析法可以获取ICESat任务期间的南极冰盖高程变化结果(见图 1)(史红岭,2010; 史红岭,2011).ICESat所得南极冰盖高程变化的空间分辨率约为100 km,时间跨度为2003年02月至2009年10月,高程变化负增长的区域主要集中在Amundsen海岸的Pine Island/Thwaites冰川周围、Queen Maud Land北部的Jutulstraumen等冰川、Wilkes Land北部的Totten冰川以及横贯南极山脉北部的Mercer冰流等地区,而高程变化正增长的区域主要分布在南极半岛的Palmer Land以及西南极的Kamb冰流和东南极广大内陆等地区.

表 1(Table 1)

表 1 激光测高仪的工作任务期汇总 Table 1 Laser operation campaign summary 序号 激光器任务期 起始时间 结束时间 工作天数 1 L1 2003-02-20 2003-03-29 38 2 L2A 2003-09-25 2003-11-18 54 3 L2B 2004-02-17 2004-03-20 33 4 L2C 2004-05-18 2004-06-20 34 5 L3A 2004-10-03 2004-11-08 37 6 L3B 2005-02-17 2005-03-24 36 7 L3C 2005-05-20 2005-06-22 34 8 L3D 2005-10-21 2005-11-23 34 9 L3E 2006-02-22 2006-03-27 34 10 L3F 2006-05-24 2006-06-25 32 11 L3G 2006-10-25 2006-11-27 34 12 L3H 2007-03-12 2007-04-14 34 13 L3I 2007-10-02 2007-11-04 36 14 L3J 2008-02-17 2008-03-21 34 15 L3K 2008-10-04 2008-10-18 14 16 L2D 2008-11-25 2008-12-17 23 17 L2E 2009-03-09 2009-04-11 34 18 L2F 2009-09-30 2009-10-11 12

表 1 激光测高仪的工作任务期汇总 Table 1 Laser operation campaign summary

图 1

Fig. 1

图 1 (a) ICESat数据估算的南极冰盖表面高程变化(2003-02—2009-10);(b)对应的误差 Fig. 1 (a) The height change trend for AIS estimates by ICESat data over the period between February of 2003 and October of 2009;(b)Corresponding uncertainties

为了更好地跟ICESat数据进行匹配,需要尽可能选择工作任务期跟ICESat相一致的GRACE月重力场模型,因此根据表 1所列的ICESat工作任务期将GRACE GSM RL05产品中2003-02和2003-03平均作为第1任务期,2003-09、2003-10和2003-11平均作为第2任务期,2004-02和2004-03平均作为第3任务期,2004-05和2004-06平均作为第4任务期,2004-10和2004-11作为第5任务期,2005-02和2005-03平均作为第6任务期,2005-05和2005-06平均作为第7任务期,2005-10和2005-11平均作为第8任务期,2006-02和2006-03平均作为第9任务期,2006-05和2006-06平均作为第10任务期,2006-10和2006-11平均作为第11任务期,2007-03和2007-04平均作为第12任务期,2007-10和2007-11平均作为第13任务期,2008-02和2008-03平均作为第14任务期,2008-10月作为第15任务期,2008-11和2008-12平均作为第16任务期,2009-03和2009-04平均作为第17任务期,2009-10作为第18任务期.这样将GSM数据产品归算到18个工作任务期后GRACE数据便与ICESat数据形成了良好的一致性,有助于后期的联合处理.

为了测试不同GRACE数据产品和后处理方法对GIA信号的影响情况,本文尽可能选择差别较大的GRACE结果,共选取了三种不同的GRACE数据产品并对应三种不同的后处理方法,分别为: CSR-60(Slepian法)(高春春,2015)、CSR-96(DDK5滤波)(Kusche et al., 2009)和GFZ(组合滤波法)(高春春等,2015).此外,由于GRACE对重力场的短波信息相对不敏感,而且球谐截断和后处理也会使重力场短波信号进一步丢失,从而导致GRACE结果出现泄漏误差,一方面,南极区域外的非冰盖质量变化信号会泄漏到南极冰盖内,形成内泄漏(leakage-in)误差; 另一方面,南极区域内的冰盖质量变化信号会削弱或泄漏到区域外,造成外泄漏(leakage-out)误差.内泄漏误差主要来源于冰盖以外的陆地水、非潮汐的海洋及冰盖附近冰帽等质量变化信号；陆地水文信号引起的泄漏,本文采用GLDAS(Global Land Data Assimilation System)模型(Rodell et al., 2004)进行模拟改正；海洋环流引起的泄漏,本文采用OGCM(Ocean General Circulation Model)数据(Dobslaw et al., 2013)进行估算扣除；对于海陆质量交换引起的泄漏,采用GLDAS模型计算的陆地水质量平衡将其平铺到海洋上进行估算扣除.而外泄漏误差本文采用一个新的信号恢复法——自适应边界选取法(高春春,2015)进行改正,以避免尺度因子恢复法(Velicogna and Wahr, 2006)引起的信号变形,泄漏误差改正的结果见表 2.

表 2(Table 2)

表 2 泄漏误差改正值(单位: Gt/a) Table 2 Correction values of leakage errors(Unit: Gt/a) 数据(方法) 内泄漏误差 外泄漏误差 陆地水文 海洋 海陆质量交换 CSR-60(Slepian法) +4.7 +9.1 +4.6 -6.2 CSR-96(DDK5滤波) +6.3 +15.0 +6.2 -12.4 GFZ(组合滤波法) +6.0 +15.3 +6.8 -12.7

表 2 泄漏误差改正值(单位: Gt/a) Table 2 Correction values of leakage errors(Unit: Gt/a)

GRACE估算的南极冰盖质量变化(见图 2a—2c)显示冰盖质量增加和减小的区域跟ICESat所得的冰盖高程上升和下降的区域基本一致,但在Queen Maud Land和横贯南极山脉北部区域,ICESat结果有明显的负增长,而GRACE没有体现.Queen Maud Land北部区域的不一致可能跟ICESat的不确定性或GRACE的空间分辨率不足有关,而横贯南极山脉北部的冰川区域,可能跟GIA信号在该区域对GRACE结果影响较大有关.

图 2

Fig. 2

图 2 GRACE数据估算的南极冰盖质量变化(2003-02—2009-10)及其不确定性 (a)采用CSR-60数据和Slepian法；(b)采用CSR-96数据和DDK5滤波法；(c)采用GFZ数据和组合滤波法；(d—f)图a—c对应的误差. Fig. 2 The mass change trends and corresponding uncertainties for AIS estimate from GRACE over the period between February of 2003 and October of 2009 (a)CSR-60 solutions and the Slepian method;(b)CSR-96 solutions and the DDK5; (c)GFZ solutions and the combined filter method;(d—f)the errors for Fig.a—c.

ICESat所观测到的南极冰盖表面高程变化dH/dt包含了很多因素的影响(Ligtenberg et al., 2011)：

(1)

v_acc为雪积累率,v_me为雪融化率,v_fc粒雪压实率,v_ice冰流速率,v_b基岩运动.其中粒雪压实率v_fc只影响ICESat所观测的冰盖高程变化,并不影响GRACE所观测的冰盖质量变化,如果不加以考虑,会降低分离GIA信号的可靠性,本文采用南极粒雪密实化FDM模型(Ligtenberg et al., 2011)来消除这种影响,FDM数据的时间分辨率为2天,依然参照ICESat工作任务期(表 1),将FDM的时间归算到与ICESat一致的工作任务期,然后进行时间序列分析,最后获得FDM冰盖表面高程变化趋势(见图 3a).

图 3

Fig. 3

图 3 (a) FDM数据估算的南极冰盖表面高程变化趋势(2003-02—2009-10)；(b)对应的误差 Fig. 3 Surface height change trend estimates from FDM for the period between February of 2003 and October of 2009;(b)corresponding uncertainties

ICESat在南极内陆冰盖观测到的高程变化主要包含两部分: 冰雪运动引起的冰盖高度变化ĥ_ice/snow和GIA引起的地壳隆升ĥ_GIA,

(2)

乘上密度转换为地表质量变化为

(3)

其中ρ_surf为冰盖表面的粒雪密度(见图 4a)(Ligtenberg et al., 2011),ρ_rock为岩石密度,在内陆冰盖为4000 kg·m^-3,在冰架区域为3400 kg·m^-3,因为在冰架区域岩石密度较小(Riva et al., 2009).由于GRACE监测到的质量变化为总质量变化$\hat{m}$_sum,即$\hat{m}$_GRACE=$\hat{m}$_sum,那么联合公式(2)和公式(3)就可以解得(Riva et al., 2009):

(4)

但是根据公式(1)知道,冰雪运动引起的冰盖高度变化ĥ_ice/snow是一个很复杂的过程,如果简单地乘上表面密度ρ_surf会造成结果不精确,例如在冰川区域ĥ_ice/snow的变化主要是由于冰排放引起,需要乘以冰的密度ρ_ice=917 kg·m^-3.因此这里引入FDM数据计算的冰盖表面高程变化ĥ_firn和对应的冰盖表面质量平衡SMB数据(见图 4b)(Lenaerts et al., 2012). 如果ĥ_ICESat－ĥ_firn<0且ĥ_ICESat－ĥ_firn|>时,则认为ĥ_ICESat－ĥ_firn的变化是由冰排放引起；如果ĥ_ICESat－ĥ_firn>0且|ĥ_ICESat－ĥ_firn|>2δ_h时,则认为是SMB数据被低估,需要加回；其他情况则认为ĥ_ICESat－ĥ_firn=0(图 5a展示了ĥ_ICESat－ĥ_firn,图 5b为2δ_h),那么公式(4)就可以改进为下式(Gunter et al., 2014)：

(5)

其中

(6)

图 4

Fig. 4

图 4 (a)南极冰盖表面密度;(b)南极冰盖表面质量平衡SMB(2003—2009),白线内区域为LPZ Fig. 4 (a)Surface density for AIS;(b)Surface mass balance(SMB)for AIS from 2003 to 2009,outline of the LPZ is white line

从图 5a可知ĥ_ICESat－ĥ_firn<0且|ĥ_ICESat－ĥ_firn|>2δ_h的区域集中在南极Amundsen海岸和横贯南极山脉北部冰川运动较为显著的地区,因此ĥ_ICESat－ĥ_firn<0的部分乘上冰密度是非常必要的. 而ĥ_ICESat－ĥ_firn>0且|ĥ_ICESat－ĥ_firn|>2δ_h的区域主要集中在南极半岛、西南极的Edith Ronne Land和东南极的Wilkes Land等降水较多的区域,SMB数据可能被低估.

图 5

Fig. 5

图 5 (a)ICESat冰盖高程变化与FDM之差ĥICESat－ĥfirn；(b)2δh Fig. 5 (a) Surface height change differences between the ICESat and FDM: ĥICESat－ĥfirn；(b)2δh

联合GRACE、ICESat、FDM、SMB和南极表面密度等数据,采用公式(5)即可分离出GIA信号ĥ_GIA,但值得注意的是,不同数据之间的空间分辨率并不一致,其中GRACE数据空间分辨率最低,因此联合处理时需要将空间分辨率较高的数据统一平滑到跟GRACE数据一致,这里采用高斯滤波进行空间平滑.虽然不同数据在空间分辨率上进行了统一,但不同类型的数据进行融合处理时,难免会在一些区域存在大的匹配误差,而GIA信号在南极地区主要表现为上升趋势,因此将结果中小于-2 mm·a^-1的值作为粗差剔除.为了便于区分,将三种不同GRACE结果获取的GIA信号分别命名为Ga、Gb和Gc,其空间分布见图 6a－6c.从图 6a－6c可以看出不同GRACE结果分离的GIA信号在空间分布的形态基本上一致,GIA垂直速率变化较大的区域主要集中在西南极的Ronne冰架、横贯南极山脉北部和Ross冰架,以及东南极Queen Maud Land、Enderby Land、Amery冰架、Princess Elizabeth Land和Kalser WilhelmⅡLand等地区.

图 6

Fig. 6

图 6 联合方法估算的GIA结果(a—i)以及GIA与GPS速率的比较,圆圈的半径为ĤGIA－ĥGPS,品红色圆圈代表ĤGIA－ĥGPS>0,蓝色圆圈代表ĤGIA－ĥGPS<0 Fig. 6 Antarctic GIA vertical rates(a—i)computed from the combination approach and comparison of the vertical rates from the GPS stations and the GIA uplift rates,the radius of the circles is ĤGIA－ĥGPS, magenta circles represent ĤGIA－ĥGPS>0, blue circles represent ĤGIA－ĥGPS<0

除了粗差,不同数据的联合处理不可避免会造成所得的结果出现一定系统偏差,本文采用Gunter等(2014)提出的LPZ(Low-Precipitation Zone)偏差改正法来消除系统偏差的影响(Gunter et al., 2014),LPZ即为图 4b中白线所圈的区域,由于GIA信号在LPZ区域信号非常小,理论上假设为0,那么Ga、Gb、Gc中LPZ区域的信号即为偏差改正值(见表 3),将Ga、Gb、Gc分别减去偏差改正值即可得到经过LPZ偏差改正后的GIA结果: Ga-LPZ、Gb-LPZ、Gc-LPZ(见图 6d—6f).

表 3(Table 3)

表 3 LPZ偏差改正值和GPS验证结果 Table 3 LPZ bias correction and the result of GPS validation GIA LPZ偏差改正(mm·a-1) GPS验证 WRMS(mm·a-1) RMS(mm·a-1) WA+AP EA 南极 WA+AP EA 南极 Ga - 2.4 3.4 3.2 2.6 3.8 3.3 Gb - 2.5 3.4 3.3 2.6 3.8 3.4 Gc - 1.7 3.1 2.9 1.9 3.7 3.1 Ga_LPZ 0.99 2.0 2.7 2.6 2.6 3.2 3.0 Gb_LPZ 0.88 2.2 2.8 2.7 2.6 3.3 3.0 Gc_LPZ 0.52 1.9 2.7 2.6 2.1 3.4 2.9 Ga_LPZ_w - 1.7 2.4 2.3 2.4 2.9 2.7 Gb_LPZ_w - 1.8 2.5 2.4 2.4 2.9 2.7 Gc_LPZ_w - 2.1 2.3 2.3 2.2 3.0 2.6 ICE5G - - - - 2.7 2.8 2.8 W12a - - - - 3.0 2.3 2.7 G14 - 3.0 1.8 2.1 3.0 2.3 2.6 Ga_LPZ_w_fit - 1.5 1.9 1.8 1.7 2.2 2.0 Gb_LPZ_w_fit - 1.7 1.9 1.8 1.7 2.2 2.0 Gc_LPZ_w_fit - 1.8 1.6 1.6 1.8 2.1 2.0 ICE5G_fit - - - - 1.7 1.9 1.8 W12a_fit - - - - 2.1 1.8 2.0 G14_fit - 2.5 1.5 1.7 2.4 2.0 2.2 注： WA+AP：西南极+南极半岛,EA：东南极.

表 3 LPZ偏差改正值和GPS验证结果 Table 3 LPZ bias correction and the result of GPS validation

为了评估所计算GIA结果的可靠性,还需对其进行精度评价,其中ICESat和FDM的不确定性见图 1b和图 3b,SMB数据的不确定性为模型本身的10%(Rignot et al., 2008),ρ_rock的不确定性为100 kg·m^-3(Wahr et al., 2000),ρ_α的不确定性为其本身的10%(Gunter et al., 2014),GRACE估算结果的误差(见图 2d—2f)需要综合考虑泄漏误差、大气误差、混频误差、测量误差等联合估算(高春春等,2015; 高春春,2015),各分项误差获取后,即可根据公式(5)按照误差传播定律计算Ga-LPZ、Gb-LPZ和Gc-LPZ的不确定性(见图 7a—7c).

图 7

Fig. 7

图 7 GIA结果的不确定性 Fig. 7 The uncertainties of GIA vertical rates

通过图 7a—7c可知,GIA信号误差较大的区域主要集中在南极沿海地区,因为这些区域是GRACE结果泄漏、混频等误差较为严重的区域,也是ICESat结果和FDM数据误差较大的区域,而利用这些数据联合估算的GIA信号在这些区域的误差自然就相对较大.我们研究发现误差较大的区域一般都集中在冰流速度较大的地区(见图 8a),因为冰流速度越大的区域冰盖高程变化越复杂,ICESat所监测的结果不确定性就越大,冰流速度大的区域质量变化信号也相对剧烈,GRACE结果中的泄漏误差也就越严重,而且冰流速度大的区域一般分布在沿海地区,气候变化多样,FDM和SMB数据误差也相对较大,GRACE的大气误差也较为明显,因此本文提出利用南极地区冰流速度分布数据(Rignot et al., 2011)对GIA结果进行加权处理,即冰流速度加权改正法.

图 8

Fig. 8

图 8 (a)南极冰流速度分布；(b)利用南极冰流速度计算的权值 Fig. 8 (a)Antarctic ice velocity;(b)Weight computed based on Antarctic ice velocity

具体做法：首先将冰流数据采用与ICESat等数据一样的高斯滤波进行平滑,这样得到一个与其他数据空间分辨率一致的冰流速度v_{ice flow},然后假定v_{ice flow}=100 m·a^-1时权值为0,v_{ice flow}=0 m·a^-1时权值为1,即

(7)

最后计算的权值w_{ice flow}见图 8b 那么将Ga-LPZ、Gb-LPZ、Gc-LPZ及其对应的误差都乘上w_{ice flow},即可得到经过加权改正后的GIA结果Ga-LPZ-w、Gb-LPZ-w、Gc-LPZ-w(见图 6g—6i)和相应的误差(见图 7d—7f).经过加权改正之后GIA信号误差(图 7d—7f)明显小了许多,特别是在Amundsen海岸冰盖质量变化最剧烈的地区误差得到了有效的控制.

为了验证分离GIA信号的可靠性,本文引入ICE5G、W12a两个GIA模型公布的地壳隆升速率,以及Gunter等利用GRACE CSR RL05 DDK数据和ICESat数据分离的南极GIA隆升速率(Gunter et al., 2014),这里简称之为G14,将这三种GIA信号作为外围的比较数据,然后利用南极35个GPS观测站(观测站位置分布为图 6圆圈中心)观测到的地壳垂直变化速率(Thomas et al., 2011)对本文估算的GIA信号和上述三种外围GIA信号分别进行GPS验证.GPS观测到的南极地壳变化速率包含GIA和地表质量变化引起的弹性形变两部分,在扣除弹性形变后就可以利用GPS数据作为地壳隆升的直接观测值对GIA信号Ĥ_GIA进行精度验证(Thomas et al., 2011),本文采用的验证方法为加权均方根WRMS(Whitehouse et al., 2012)和RMS,WRMS计算公式为

(8)

其中

(9)

GIA信号和ICE5G、W12a、G14的GPS验证结果列于表 3,由验证结果可知：①在东南极LPZ偏差改正法和冰流速度加权改正法能够有效减小WRMS和RMS值,也就是有效提升GIA信号的可靠性和精度,而在西南极和南极半岛只有Ga和Gb模型精度有所提高,而Gc精度则有所下降,就整个南极而言,两种方法总体改善了GIA信号的精度；②本文分离的GIA在西南极和南极半岛精度要明显优于ICE5G、W12a和G14这三个外围GIA模型,但是在东南极精度要明显差一些,在整个南极Ga_LPZ_w、Gb_LPZ_w和Gc_LPZ_w与外围GIA模型精度相当.

从图 6可以看出：本文分离的GIA信号与GPS速率的差值在东南极基本都大于0,在南极半岛都小于0,在西南极地区主要是大于0,也有少数站点小于0,这说明GIA结果存在着区域性的系统偏差.为了消除这种区域性系统偏差,本文采用GPS拟合法对GIA结果进行进一步精化,针对南极地区的区域特点,采用球谐函数进行曲面拟合是相对较优的选择,根据GPS的站点数量采用4阶球谐函数作为拟合模型进行GPS拟合效果最佳,其中球谐拟合的数学模型为

(10)

利用公式(10)根据最小二乘原理即可求得未知系数A_nm和B_nm,利用求得的球谐系数A_nm和B_nm即可以获取南极地区的拟合曲面ΔĤ,本文采用Ga_LPZ_w、Gb_LPZ_w、Gc_LPZ_w三个GIA结果和ICE5G、W12a、G14三个外围GIA信号一起进行GPS拟合,计算的拟合曲面见图 9.拟合曲面基本上表现出了GIA信号所存在的区域性系统偏差,将GIA信号加上ΔĤ即可得到GPS拟合改正后的GIA信号(见图 10),由拟合后的WRMS和RMS值(见表 3)可知GPS拟合改正使GIA精度得到了进一步改善,特别是东南极地区,三个外围GIA信号在GPS拟合改正后精度也明显得到了改善,但改善主要在西南极和南极半岛较为明显,经过GPS拟合改正后的所有GIA信号精度基本上趋于一致.从图 10可知经过拟合改正后的GIA结果Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit相比未改正前在东南极的上升趋势变的不明显,GIA信号的区域性偏差得到了很好的改正.

图 9

Fig. 9

图 9 GPS球谐拟合曲面 Fig. 9 GPS spherical fitting surfaces

图 10

Fig. 10

图 10 GPS球谐拟合改正后的GIA结果 Fig. 10 Antarctic GIA vertical rates corrected by GPS spherical fitting

经过联合方法分离、LPZ偏差改正、冰流速度加权改正以及GPS球谐拟合改正一系列处理后的GIA结果Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit(见图 10a－10c)显示: GIA信号上升速率较大的区域主要集中在西南极的Ronne冰架、横贯南极山脉北部和Ross冰架以及东南极Kalser Wilhelm Ⅱ Land等地区.东南极的广大内陆地区主要呈现微弱的下降趋势,这种下降趋势在W12a模型中有更为明显的表现(见图 10e).Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit相比同样采用GRACE和ICESat数据分离的G14最大的区别在于Amundsen海岸和东南极Kalser Wilhelm Ⅱ Land沿海区域,G14在Amundsen海岸表现出较强的隆升趋势,Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit并没有体现；Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit在Kalser Wilhelm Ⅱ Land沿海区域表现出较大的上升速率,而G14并没有那么强烈.将Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit所得的高程变化与岩石密度相乘即可转换为地表质量变化并从GRACE估算结果中扣除,可以得到2003年02月至2009年10月GRACE观测的南极冰盖质量变化趋势(见图 11a—11c).另外还可利用Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit的高程变化作为对ICESat观测的GIA改正,采用以下公式：

(11)

可以计算得到ICESat观测的南极冰盖质量变化趋势(见图 11d—11f),数值结果列于表 4.

图 11

Fig. 11

图 11 GRACE(a—c)和ICESat(d—f)估算的南极冰盖质量变化趋势(2003-02—2009-10) Fig. 11 Antarctic ice sheet mass change trends estimated from GRACE(a—c)and ICESat(d—f) for the period between February of 2003 and October of 2009

表 4(Table 4)

表 4 GRACE和ICESat估算的南极冰盖质量变化(2003-02—2009-10) Table 4 Antarctic ice sheet mass change estimated from GRACE and ICESat for the period between February of 2003 and October of 2009 GIA GIA(Gt/a) GRACE-GIA(Gt/a) ICESat-GIA(Gt/a) ICE5G 146.8±66.2 -170.4±69.6 -84.6±23.1 W12a 77.2±14.8 -97.2±26.2 -76.0±20.1 G14 75.6±47.7 -100±44.0 -81.5±21.6 Ga_LPZ_w 80.2±51.4 -104.0±56.6 -81.1±21.5 Gb_LPZ_w 86.9±49.3 -107.4±54.2 -81.4±21.6 Gc_LPZ_w 87.4±48.0 -110.7±52.6 -81.2±21.5 ICE5G_fit 106.5±66.2 -130.1±69.6 -81.2±23.1 W12a_fit 39.0±14.8 -59.0±14.8 -71.9±20.1 G14_fit 44.5±47.7 -68.9±44.0 -77.9±21.6 Ga_LPZ_w_fit 46.6±51.4 -70.4±56.6 -77.8±21.5 Gb_LPZ_w_fit 47.5±49.3 -68.1±54.2 -77.6±21.6 Gc_LPZ_w_fit 38.4±48.0 -61.7±52.6 -76.3±21.5

表 4 GRACE和ICESat估算的南极冰盖质量变化(2003-02—2009-10) Table 4 Antarctic ice sheet mass change estimated from GRACE and ICESat for the period between February of 2003 and October of 2009

从图 11可以看出2003年02月至2009年10月GRACE观测到的南极冰盖质量消融最大的区域集中在Amundsen海岸和南极半岛北部的冰川区域,而ICESat结果显示在Amundsen海岸和横贯南极山脉北部的冰川区域,两者最大的差别在于横贯南极山脉北部的Mercer冰流区域,这可能跟GIA信号在该区域较大有关,因为GIA对GRACE结果的影响相比对ICESat结果的影响要显著的多.从表 4可知：①未经GPS拟合改正的GIA结果Ga_LPZ_w、Gb_LPZ_w和Gc_LPZ_w所计算的地表质量变化分别为80.2±51.4 Gt/a、86.9±49.3 Gt/a和87.4±48.0 Gt/a,彼此相差很小,经过GPS拟合后的GIA所计算的地表质量变化分别为46.6±51.4 Gt/a、47.5±49.3 Gt/a和38.4±48.0 Gt/a,拟合改正后量级减小30～40 Gt/a；②利用Ga_LPZ_w、Gb_LPZ_w和Gc_LPZ_w分别进行GIA改正后的GRACE结果彼此相近,三者的平均质量变化趋势为-107.4±54.5 Gt/a,与G14的-100±44 Gt/a(Gunter et al., 2014)非常接近,而利用Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit进行GIA改正的GRACE结果平均质量变化趋势为-66.7±54.5 Gt/a,与G14的-68.9±44 Gt/a更为接近；③不同的GIA模型对ICESat的结果影响很小,最大差异不超过10 Gt/a,而对GRACE的结果影响很大,最大差异可达75 Gt/a；④采用Ga_LPZ_w、Gb_LPZ_w和Gc_LPZ_w进行GIA改正后的GRACE和ICESat结果彼此相差约25Gt/a,但是采用经过GPS拟合改正后的Ga_LPZ_w_fit、Gb_LPZ_w_fit和Gc_LPZ_w_fit,GRACE和ICESat结果差异缩小为10 Gt/a,可见GPS拟合改正法有效降低了GIA的不确定性,使得不同源的卫星数据估算的南极冰盖质量平衡趋于一致.

本文的GRACE结果平均为-66.7±54.5 Gt/a与Gunter等(2014)的结果-100±40 Gt/a有所差别,主要原因在于本文采用了GPS球谐拟合对GIA信号进行了改正,而未改正前本文的结果是-107.4±54.5 Gt/a与Gunter解算的结果基本一致,与King等(2012)的GRACE结果-69±18 Gt/a非常接近(King et al., 2012),与贾路路等(2011)的-82±18 Gt/a(贾路路等,2011)、罗志才等(2012)的-80 Gt/a(罗志才等,2012)、Luthcke等(2013)的-81±26 Gt/a(Luthcke et al., 2013)、Velicogna和Wahr(2013)的-83±49 Gt/a(Velicogna and Wahr, 2013)都相差不到20 Gt/a,GRACE结果彼此差异的原因主要在于选择GIA改正信号的不同(高春春等,2015).本文的ICESat结果平均为-77.2±21.5 Gt/a,低于Groh等(2014)的-126±20 Gt/a(Groh et al., 2014),高于李斐等(2016)的-44±21 Gt/a(李斐等,2016),主要原因在于ICESat数据解算的方法不同,本文采用的是交叉点分析法,而Groh和李斐采用的是重复轨迹法,此外ICESat数据偏差改正方法的不同也会导致解算结果的差异(李斐等,2016).

本文联合GRACE卫星重力和ICESat激光测高、南极FDM和SMB等数据分离了南极GIA信号,所分离的GIA信号不依赖于不确定性很大的地球模型和冰负荷模型,是基于观测数据直接估算的GIA信号,更具有可靠性.分离过程中,为了有效降低GIA信号的不确定性,本文在前人LPZ偏差改正法的基础上提出了冰流速度加权改正法和GPS球谐拟合改正法,实现了对GIA结果的进一步精化,GPS验证结果显示双改正法能够有效提高GIA结果的精度.最后本文对分离的GIA结果进行了详细的讨论和应用分析,利用联合数据所估算的GIA对GRACE和ICESat进行改正,得到南极冰盖质量变化的趋势分别为-66.7±54.5 Gt/a(GRACE)和-77.2±21.5 Gt/a(ICESat),相比采用其他的GIA模型,本文所得的GIA使GRACE和ICESat这两种不同测量方式所得的南极冰盖质量变化结果更加趋于一致.由于目前所能获得的ICESat数据只有18个工作任务期,时间段相对较短,工作任务期较少,GRACE数据的空间分辨率也相对较低,南极GPS站点也比较少,联合方法所分离的GIA结果不确定性还相对较大,但是随着未来ICESat-2和GRACE Follow-on任务的实施,相信随着更高质量更高精度的时变重力和激光测高数据用于南极GIA信号的分离中,南极GIA结果的精度将会得到极大的改善.

感谢匿名评审专家和编辑提供的宝贵建议和帮助,感谢CSR和GFZ提供GRACE数据产品,感谢荷兰乌得勒支大学Lenaerts博士提供南极SMB数据,感谢荷兰乌得勒支海洋大气研究所Ligtenberg博士提供南极FDM数据,感谢荷兰代尔夫特理工大学Gunter博士提供G14数据,感谢英国杜伦大学Whitehouse博士提供W12a模型,感谢美国科罗拉多大学Geruo博士提供ICE5G模型,感谢英国纽卡斯尔大学Thomas博士提供南极GPS数据,特别感谢中国科学院测量与地球物理研究所许厚泽院士、王勇研究员和汪汉胜研究员在本文完成中给予的指导和意见.