地球物理学报  2017, Vol. 60 Issue (1): 37-49   PDF    
利用运动学轨道提高GRACE时变重力场解算
杨帆1,2 , 王长青2 , 许厚泽1,2 , 钟敏2 , 周泽兵1     
1. 华中科技大学物理学院地球物理研究所, 武汉 430074;
2. 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室, 武汉 430077
摘要: 基于变分方程法,本文利用GARCE高精度K波段星间测速数据KBRR,结合德国格拉茨大学发布的运动学轨道和GFZ发布的简动力学轨道作为两种伪观测值,分别解算了2005-2010年60阶全球时变重力场模型Hust-IGG01与Hust-IGG02.通过与GRACE官方机构发布的模型和其他国际主流权威模型进行对比,发现基于运动学轨道结合KBRR解算的模型Hust-IGG01优于基于简动力学轨道结合KBRR解算的模型Hust-IGG02:在重力场系数C20时间序列的统计数据上,Hust-IGG01比Hust-IGG02更接近SLR结果,在如C60C70C80以及C90等重力场低阶项上的数学统计均更接近CSR RL05;Hust-IGG01的重力场系数误差分布和GFZ RL05在同一水平,而Hust-IGG02的误差估计过于乐观;Hust-IGG02在主要质量变化区域上存在5%~10%信号低估,而Hust-IGG01能完全达到国际主流机构利用GPS观测数据的解算水平,Hust-IGG01与官方机构CSR、JPL和GFZ最新模型在格陵兰岛的冰川消融年际趋势分别是-125.4、-125.4、-127.3、-124.3 Gt·a-1,在亚马逊流域的平均等效水高周年振幅分别是17.56、17.40、17.46、17.22 cm,在撒哈拉沙漠的平均等效水高均方差分别是0.87、0.77、1.10、0.87 cm;另外在Hust-IGG01的实际应用上,本文分析了全球32个主要流域质量变化的年际趋势、周年振幅和半周年振幅三种信号模式,统计结果显示Hust-IGG01与CSR RL05结果基本吻合.
关键词: GRACE      时变重力场      运动学轨道      简动力学轨道      等效水高     
Towards a more accurate temporal gravity model from GRACE observations through the kinematic orbits
YANG Fan1,2, WANG Chang-Qing2, HSU Hou-Tse1,2, ZHONG Min2, ZHOU Ze-Bing1     
1. Institute of Geophysics, Huazhong University of Science and Technology, Wuhan 430074, China;
2. State Key Laboratory of Geodesy and Earth's Geodynamics Chinese Academy of Sciences, Wuhan 430077, China
Abstract: Based on the GRACE Level 1b raw data from 2005 to 2010, we have successfully produced an unconstrained monthly gravity field model (Hust-IGG01) up to d/o 60. Unlike the official data-processing centers, we employ the kinematic orbits instead of the GPS measurements as the pseudo observations. Meanwhile, an alternative model (Hust-IGG02) using the reduced-dynamic orbits as the pseudo observations is provided as well. We aim to understand the impacts of orbit pseudo observations on the accuracy of the ultimate gravity products. To this end, Hust-IGG01 and Hust-IGG02 are fully compared to each other, such that we are able to identify which type of orbit pseudo observations is more desired for the gravity inversion. Experiments demonstrate that Hust-IGG01 performs better in terms of signal to noise level than Hust-IGG02 in the following aspects:(a) Hust-IGG01 improves the estimation of geopotential coefficients at low degrees, e.g., C20 C60 C70 C80 C90, which are closer to SLR or CSR RL05 results; (b) the induced formal error of Hust-IGG01 is appropriate and comparable to that of GFZ RL05a, while Hust-IGG02 poses too optimistic formal error; (c) over three typical regions of interest (Amazon, Greenland and Sahara), the mass variation derived from Hust-IGG02 has been under-estimated by about 5%~10% with respect to those from the official products, while Hust-IGG01 has achieved a fairly comparable accuracy. The latter is supported by the numerical results such as:the yearly trend of glacial melting over Greenland derived from Hust-IGG01, CSR RL05, GFZ RL05a and JPL RL05 are -125.4 Gt·a-1, -125.4 Gt·a-1, -127.3 Gt·a-1 and -124.3 Gt·a-1, respectively; the annual amplitude of mass change in terms of EWH (equivalent water height) over Amazon is 17.56 cm, 17.40 cm, 17.46 cm and 17.22 cm, respectively; the RMS of mass change over Sahara desert is 0.87 cm, 0.77 cm, 1.10 cm and 0.87 cm, respectively. An additional validation is undertaken as well, to investigate the performance of Hust-IGG01 on the scale of basins, and the results demonstrate that the induced annual amplitude, semi-annual amplitude and yearly trend of mass variations agree with those of CSR RL05 over the 32 selected major river basins. In summary, our comparisons above suggest that an appropriate kinematic orbit is more beneficial than the reduced-dynamic orbits, for the accurate gravity recovery from the GRACE observations..
Key words: GRACE      Temporal gravity field      Kinematic orbit      Reduced-dynamic orbit      Equivalent water height     
1 引言

GRACE(Gravity Recovery And Climate Experiment)重力卫星探测任务由美国NASA和德国DLR发起(Tapley et al.,2004),于2002年发射,迄今为止,GRACE累积的十多年的观测数据为我们提供了前所未有的地球质量变化信息.尤其由三大官方机构美国德克萨斯大学空间中心CSR(Center for Space Research)、美国宇航局喷气推进实验室JPL(Jet Propulsion Laboratory)和德国地学研究中心GFZ(GeoFoschungs Zentrum)根据GRACE星载加速度计、行星姿态摄像仪、K波段测距系统以及GPS接收机等(Touboul et al.,1999; Dunn et al.,2003)所采集的数据解算出的Level 2地球月解时变重力场产品(Bettadpur,2007),通过更新原始数据质量或者提高数据处理流程,逐渐发展到目前广泛使用的第五版,已经达到了相当高的精度,并且被广泛运用在海洋学(Chambers et al.,2004)、水文学(Richey et al.,2015)、地震学(Han and Simons,2008)等等研究领域(Kusche et al.,2009).除了三大官方机构之外,国内外还有其他单位发布了自己的产品,根据建立观测方程的方式分类:德国格拉茨大学和波恩大学用短弧长法分别解算了ITSG和ITG模型(Mayer-Gürr,2006; Mayer-Gürr et al.,2014),荷兰代尔夫特理工大学利用加速度法得到了DMT模型(Liu,2008),瑞士伯尔尼大学利用天体力学法得到了AIUB模型(Jäggi et al.,2012),法国空间中心变分方程法解算得到的正则化重力场模型GRGS(Bruinsma et al.,2010),以及近期美国俄亥俄大学用能量法也成功解算到时变重力场信号(Shang,2015),此外,在国内,同济大学根据改进短弧长法和加速度法解算得到了Tongji系列模型(Chen et al.,2015a2015b),中国科学院测量与地球物理研究所利用短弧长法和变分方程方法分别解算得到了IGG-CAS和IGG系列模型(冉将军等,2014; 王长青等,2015),武汉大学也利用变分方程法解算得到了时变重力场模型WHU-Grace01s(Zhou et al.,2015).

然而参考国内公开发表的论文计算结果,在发生较大质量变化的区域,国内模型和国外权威模型相比存在一定程度的信号低估:Tongji模型、CSR RL05和JPL RL05在亚马逊流域以及其周边一定范围([21° S,5° N],[45° W,80° W]),从2003—2011年的周年振幅以等效水高表示分别是26.0 cm、27.6、27.3 cm(Chen et al.,2015);中国科学院测量与地球物理研究所IGG-CAS模型、CSR RL05和JPL RL05在亚马逊流域,从2004—2010年周年振幅是14、16、16 cm(冉将军,2014);中国科学院测量与地球物理研究所IGG模型、CSR RL05和JPL RL05在亚马逊流域,从2005—2010年周年振幅是14.7 cm、17.1、16.9 cm(Wang et al.,2015).可以发现,官方单位CSR RL05、JPL RL05的结果相当一致,而相比之下国内反演出的时变信号存在一定程度低估,这种“信号低估”可能表明采用的数据处理方案并未充分完整的挖掘出GRACE的潜力.

在以上国内外各机构发布的模型中,三大官方模型CSR RL05、GFZ RL05a、JPL RL05以及法空局GRGS模型均采用GPS数据联合K波段星间测距数据反演时变重力场,瑞士伯尔尼大学AIUB模型、德国波恩大学ITG模型以及德国格拉茨大学ITSG模型均采用运动学轨道联合K波段星间测距数据反演时变重力场,而目前国内同济大学Tongji模型、中国科学院测量与地球物理研究所IGG-CAS和IGG模型均采用简动力学轨道联合K波段星间测距数据反演时变重力场.因此作者认为目前国内解算信号偏弱的一个可能原因是使用简动力学轨道替代了GPS观测数据作为伪观测值.尽管目前GRACE任务所得到时变重力场中高阶信号基本由K波段观测数据体现,但是GPS或者轨道数据对重力场低阶信号的约束作用仍然不可忽视(Schwintzer and Reigber,2002).事实上,简动力学轨道在GPS定轨过程中引入了较强的先验信息和一定程度的随机信息,对于轨道参数的估计约束比较强(Jäggi et al.,2007),并且这些定轨过程中引入的先验信息还会传递到重力场反演过程中,最终影响到重力场参数的解算(Meyer et al.,2016),比如会造成解算信号相比直接利用GPS作为观测值的解算结果偏小(部分信号被轨道参数吸收).然而,同样利用GPS定轨得到的运动学轨道尽管整体误差比简动力学轨道大,但是在定轨过程中却并未引入动力学模型,不包含先验信息(švehla and Rothacher,2005),理论上在时变重力场解算过程中用运动学轨道来代替GPS数据作为观测值,可以减少先验信息对解算信号的压制,从而得到更加真实的时变重力场信号.

为了比较不同类型轨道数据替代GPS作为伪观测值反演重力场的效果,本文首先针对简动力学轨道、运动学轨道、K波段观测数据的数据特性进行分析,然后基于变分方程方法,采用自主开发的重力场解算软件,分别利用运动学轨道与简动力学轨道做为伪观测值,并结合KBRR(K-band range rate)观测数据解算了60阶时变重力场Hust-IGG01与Hust-IGG02模型;将以上两模型与以GPS数据为观测值的CSR RL05、以运动学轨道为观测值的AIUB、以简动力学轨道为观测值的Tongji等时变重力场模型进行对比. 通过以上实验我们将验证运动学轨道对GRACE时变重力场反演精度的提高,并且验证国内普遍采取的GPS定轨和重力场解算分开施行的做法是合理可行的(Han and Xu,2007).

2 GRACE数据处理方案

JPL,GFZ与CSR处理GRACE数据的一般流程:第一步,单独利用GPS观测数据对GRACE双星进行精密定轨,仅以GRACE轨道相关参数为未知数建立观测方程,以获得比如双星起始位置、加速度计矫正参数、GPS相关参数等等参数的初值;第二步,联合GPS和K波段测距测量值作为观测值,以重力场Stokes参数和上述部分轨道相关参数同时做为未知数再次建立观测方程,同时解算所有参数.在上述的处理过程中,不同机构处理数据的流程略有不同,比如:GFZ RL05版本模型在解算的第二步固定了所有轨道参数,仅仅解算了重力场Stokes参数,造成了隐式正则化(Meyer et al.,2016),压制了重力场信号,于是在后续发布的RL05a中第二步同时解算了所有参数,获得了更真实的时变信号;Luthcke等(2006)在第一步联合GPS和K波段观测数据共同校正轨道参数,然后在第二步固定了大部分轨道参数,仅仅利用K波段观测值成功解算了重力场Stokes参数;波恩大学、格拉茨大学、伯尔尼大学各自发布的模型则在第二步中,分别用各自计算的运动学轨道取代GPS观测值,然后联合K波段测距观测值来解算部分轨道参数和Stokes参数;国内同济大学和中国科学院测量与地球物理研究所的策略则是在第一步和第二步均采用GFZ发布的简动力学轨道代替GPS数据,从而避免了第一步GPS定轨过程中的复杂参数估计.本文参照国内机构做法,分别采取两种不同类型的轨道数据代替GPS观测值反演重力场,探讨轨道数据对反演结果造成的影响.

2.1 背景模型和参数估计方案

本文利用变分方程法,以3 h为一个弧段,建立观测方程,观测方程的建立参考相关文献(Wang et al.,2015),不再赘述,建立法方程的过程中所利用的背景力模型参见表 1.

表 1 背景模型介绍 Table 1 Background model

利用上述力模型在每个卫星历元上建立观测方程后,对于每个弧段,局部参数包括:加速度计偏差常数项和偏差线性趋势项(Bettadpur,2009),每个弧段起始点的位置和速度矢量.特别指出,本文局部参数的估计周期和弧长(3 h)一致,另外,加速度计尺度因子直接采用了Bettadpur(2009)发布的初值,而并未作为待估参数进行解算.最终建立每个弧段的观测方程后,消除局部参数,构建成求解全局60阶Stokes参数的法方程.

2.2 观测数据质量分析

一般的,通过线性化建立观测方程后,解算地球重力场的函数模型和随机模型可看作是高斯-马尔科夫过程:

(1)

其中,y为GRACE观测值,e为观测值误差或者线性化误差等,且满足E{e}=0;D{y}为协方差矩阵;σ2是先验方差因子;A代表观测方程的设计矩阵;x是待求重力场参数,P为不同类型观测数据的权重比,在本文使用过程中被认为是一个常数.

根据公式(1),通常的无偏最小二乘估计:

(2)

GRACE原始观测值y是GPS和K波段观测数据,本文采用轨道数据替换GPS数据,特别是简动力学轨道会引入先验信息如下:

(3)

此时,$\tilde{x}$为待求重力场参数的有偏估计;y*为简动力学轨道引入的先验信息.通过这样的方式,对于KBRR和轨道数据两类观测值,可以分别建立法方程,并通过一定的权重比最终联合构建总法方程如下:

(4)

其中,NkbrrNorbit为两类观测值各自建立的法方程左阵,LkbrrLorbit为其对应的右阵;σkbrrσorbit为两类观测值的先验均方差;$\tilde{x}$为待解算重力场参数.(4)式可变换为

(5)

其中,一般认为GPS观测数据定轨精度约为1~2 cm(Ries et al.,2011),而KBRR观测值的精度约为0.1 μm·s-1,故可以赋值P=1×1010~4×1010.两类观测数据通过该权重比融合在一起组合成新的法方程,共同解算出重力场参数.可以看到,两类观测值的质量和性质将直接决定解算重力场的精度.

KBRR残差的数据特征见图 1a.KBRR残差主要包含了时变重力场中高阶信号,是GRACE反演高精度时变重力场的保证,因此其残差序列可以作为衡量重力场解算质量的指标之一,通过图 1a可以发现本文结果与GFZ的KBRR残差分布、量级和走势基本相似(Flechtner et al.,2013).除此之外,轨道或者GPS观测数据对重力场的影响也不可忽略,本文主要对两类轨道观测数据进行讨论:(1)德国格拉茨大学发布的GRACE运动学轨道,其历元非均匀分布,间隔不一,并且存在较大的粗差和较长时间的间断,需要一定的处理后才能使用;(2)而GFZ发布的GNV1b简动力学轨道历元均匀分布,每5 s一个观测值,基本无间断和较大粗差.为了观察这两类轨道数据的特征,本文利用校正后的轨道参数积分某一个月的轨道,然后分别与简动力学轨道和运动学轨道做差求得残差序列,将求得的残差序列进行对比,下面以2010年1月的数据为例.

图 1 2005—2010年KBRR残差每日RMS序列;(b)在2010年1月,运动学轨道观测值与积分轨道的残差在X,Y,Z三个方向的日均方差Ki-int X,Ki-int Y,Ki-int Z; 简动力学轨道观测值与积分轨道的残差在X,Y,Z三个方向的日均方差Rd-int X,Rd-int Y,Rd-int Z Fig. 1 Daily root mean square of K band range rate residuals span from 2005 to 2010;(b)The Root Mean Square of residuals between Kinematic orbit and integrated orbit on X,Y,Z direction respectively,denoted as Ki-int X,Ki-int Y,Ki-int Z; The Root Mean Square of residuals between Reduced Dynamic orbit and integrated orbit on X,Y,Z direction respectively,denoted as Ki-int X,Ki-int Y,Ki-int Z

根据图 1b,以运动学轨道作为观测值的三个方向的残差RMS在2 cm左右,与星载GPS观测数据定轨精度类似,而以简动力学轨道作为观测值的三个方向的残差RMS远小于1 cm,与星载GPS观测数据定轨精度不匹配.而在法方程定权的过程中,我们常用星载GPS观测数据定轨精度替代轨道残差精度,而星载GPS观测数据定轨精度是1~2 cm,从这个意义上讲,运动学轨道更接近GPS观测数据.然而简动力学轨道残差相比之下偏小,减少了随机误差的同时也可能减弱了重力场信号,而且该误差减弱并不会对解算有太大的帮助,这点可以通过残差频域图看出.图 2为2010年1月24日两类轨道数据残差和K波段星间测速(KBRR)残差的频谱分析.需要注意的是,由于KBRR往往被赋予一个相对于轨道数据1×1010~4×1010的权重比,故在能量谱中其能量被放大了1×1010~4×1010倍,最终和轨道数据的能量谱在同等水平.在图 2中发现如下:

图 2 两类轨道数据残差和K波段星间测速(KBRR)残差的频谱分析 蓝色实心线为KBRR残差的能量谱密度,红色实心线为运动学轨道作为观测值时的残差能量谱密度,黑色实心线为简动力学轨道作为观测值时的残差能量谱密度;绿色垂直线(1)代表 1 cpr频率;绿色垂直线(2)代表 2 cpr频率处;绿色垂直线(3)代表高频噪声开始占据主导的位置;绿色垂直线(4)代表运动学轨道的截止频率. Fig. 2 PSD analysis of the residuals due to the orbit pseudo observations and KBRR observations Blue solid curve denotes the power spectral density(PSD)of the K-band range rate residuals; Red solid curve denotes the PSD of orbit residuals when the reduced-dynamic orbit is taken as the observations; Black solid curve denotes the PSD of orbit residuals when the kinematic orbit is taken as the observations; Green vertical line(1)is the frequency of one cycle per revolution(CPR); Green vertical line(2)is the frequency of two CPR; Green vertical line(3)is the frequency where the high frequent noise dominates the signal; Green vertical line(4)is the truncated frequency of the kinematic orbit.

·绿色线(1)为1 CPR(cycle per revolution)效应对应频率,KBRR数据在此处往低频方向迅速下降,说明在KBRR在校正阶段低频噪声被很好的压制,但是同时低阶重力场信号也在一定程度上被压制 .

·绿色线(2)为2 CPR对应频率,如果认为2 CPR之前为重力场低阶信号,可以发现在低阶处运动学轨道信号占主导地位,简动力学轨道低阶信号偏小,而GPS数据对重力场的贡献一般认为主要在低阶.

·绿色线(3)代表信噪比开始急剧下降的频率,可以发现在绿色线(2)和(3)之间,代表了重力场信号的中高阶信号,这时KBRR对反演的贡献大于两种轨道数据.

·绿色线(3)之后代表了重力场高阶信号,可以发现,虽然简动力学轨道的高频噪声远小于运动学轨道,但是此时占据主导地位的仍是KBRR的高频噪声,因而运动学轨道的高频噪声并不会明显污染重力场高频段 .

·绿色线(4)代表了运动学轨道的截止频率(由采样率决定),运动学轨道的低采样率,恰恰使得其高频能量消失,从而减少解算重力场的误差水平,而简动力学轨道的高采样率(每5s一个观测值),不仅增加了计算负荷,并且增加了高频误差,对解算结果并无实质帮助,因为以目前GPS的精度水平来看,轨道数据没有能力反演超高阶信号.

3 利用不同类型轨道反演重力场模型结果比较

本文基于德国格拉茨大学发布的运动学轨道(Zehentner and Mayer,2013)计算了2005—2010年时变重力场模型Hust-IGG01,同时基于GFZ简动力学轨道(GNV)计算了2005—2010年时变重力场模型Hust-IGG02,并与其他国际主流模型进行了一系列对比分析.

3.1 2005—2010年72个月的平均重力场系数比较

本文首先计算每个模型在2005—2010年共72个月的平均重力场系数,然后减去背景模型GIF48后,再导出其大地水准面阶方差.在图 3中对比以下四个模型的结果:Hust-IGG01,Hust-IGG02,GFZ RL05a和CSR RL05.可以发现,在低阶项15之前,Hust-IGG01、GFZ RL05a、CSR RL05结果相当一致,而Hust-IGG02模型在低阶的信号相对偏低,这和之前简动力学轨道低阶信号偏低的特征吻合;在中高阶处,Hust-IGG02与GFZ RL05a走势较为相似,说明其高阶误差并没有因为简动力学轨道的引入而显著下降;特别注意的是,Hust-IGG01模型与CSR RL05模型的线性相关度高达0.994,而GFZ RL05a和CSR RL05的相关度是0.895.分析导致该现象的可能原因是:GFZ RL05a是90阶模型,而CSR RL05和Hust-IGG01都是60阶模型,截断引起的隐性正则化造成了60阶和90阶模型在高阶处出现明显差异.

图 3 扣除背景模型GIF48后,2005—2010年共72个月的平均重力场系数导出的每阶大地水准面高 Fig. 3 The degree geoid heights derived from the average geopotential coefficients within 72 months covering the period from 2005 to 2010,with respect to the background model GIF48
3.2 轨道数据对低阶重力场系数的影响

根据本文在第二节中对轨道数据特征的探讨可知,不同轨道数据对低阶重力场系数可能会造成影响,因此本节进一步对各模型低阶重力场系数的时间序列进行分析.首先,我们对比CSR RL05,JPL RL05,SLR,Hust-IGG01和Hust-IGG01重力场系数C20的时间序列(2005—2010年共72个月,见图 4).由于SLR观测技术估算的C20被公认为是比较准确的结果,因此这里以SLR结果为参照标准.可以看到,Hust-IGG01和Hust-IGG02模型在时间序列上更接近SLR观测结果,甚至优于CSR RL05和JPL RL05;同时根据数学统计(表 2)也发现,相比Hust-IGG02来看,Hust-IGG01与SLR的残差值的均值和标准差均微弱下降,其时间序列与SLR的线性相关度也得到一定程度提高,该结果显示Hust-IGG01的C20更加接近SLR.

图 4 2005—2010年共72个月各模型的C20时间序列 Fig. 4 Time series of C20 span from 2005 to 2010 derived from SLR,CSR RL05,JPL RL05,Hust-IGG01 and Hust-IGG02
表 2 各模型C20时间序列和SLR做差后的残差数学统计 Table 2 Statistics of differences between the following model and SLR at C20

尽管公认GRACE对C20的估计不够准确,但是对于其他低阶系数的估计是比较准确的,因此接下来以C60C70C80 以及C90为例来再次研究轨道数据对重力场低阶恢复的影响(无法对全体低阶项逐一分析,仅以上述为例,对于其他项结论类似).另外,本文引入了瑞士伯尔尼大学的AIUB模型参加对比,因为该模型的数据处理流程和本文Hust-IGG01相似(利用Bernese软件精密定轨得到运动学轨道,然后结合K波段星间测速数据反演重力场),可以作为衡量Hust-IGG01解算水平的参照.同样的,我们从CSR RL05,AIUB,Hust-IGG01与 Hust-IGG02时变重力场模型分别导出各模型2005—2010年C60C70C80C90的时间序列,并且以CSR RL05的重力场系数时间序列为参照,将其余三个模型(AIUB、Hust-IGG01与Hust-IGG02)的时间序列与之进行做差,通过对差异值进行数学统计判断各模型与CSR RL05的接近程度.结果如表 3所示,Hust-IGG01与AIUB统计结果在各种数学指标上十分接近,并且均一定程度优于Hust-IGG02.这不仅初步表明Hust-IGG01的数据处理流程正确,并且也说明运动学轨道确实对低阶项估计存在影响.

表 3 各模型C60C70C80C90时间序列与CSR RL05相应项差值的数学统计 Table 3 Statistics of differences between the following model and CSR RL05 at C60C70C80C90
3.3 重力场系数误差比较

如前所述,简动力学观测值在定轨过程引入了一定程度的约束,该约束可能会传递到重力场解算过程中,并造成解算的重力场系数的formal error估计偏小.因此,本文对Hust-IGG01和Hust-IGG02模型的formal error进行分析.图 5给出了2009年5月的结果,可以发现Hust-IGG01和Hust-IGG02各系数的formal error分布特征类似(带谐项到田谐项,误差逐渐增大;随着阶数增加,误差也逐渐增大).同时可以发现,Hust-IGG02整体的formal error比Hust-IGG01小,但是这并不意味着Hust-IGG02的误差水平更低,因为仅从数学意义上得到的formal error是一种较乐观的估计,更为真实的误差估计是校正误差(calibrated error).由于校正误差综合考虑了更多实际因素,因此一般比仅从最小二乘理论出发得到的formal error更大(Wagner et al.,2012).本文在图 6中画出了GFZ RL05a在该月的calibrated error和formal error作为参考,并与Hust-IGG01,Hust-IGG02,AIUB模型的formal error进行对比.可以发现,基于GPS数据反演的GFZ RL05a、基于运动学轨道反演的Hust-IGG01和AIUB模型的formal error在同一水平,而基于简动力学轨道Hust-IGG02的formal error过于乐观.因此,我们认为运动学轨道和GPS数据在重力场系数的误差控制上有类似效果,而简动力学轨道显然带来了一定程度的约束.

图 5 2009年5月重力场系数CnmSnm的formal error(色标以数学对数表示)分布(左Hust-IGG01,右Hust-IGG02) Fig. 5 The formal error of each geopotential coefficient Cnm and Snm(Left panel is Hust-IGG01 and the right is Hust-IGG02)
图 6 2009年5月以每阶大地水准面高表示的各个模型的标准误差或校正误差 Fig. 6 Formal errors or calibrated errors in terms of geoid height per degree on May 2009
3.4 格陵兰岛、亚马逊流域、撒哈拉沙漠典型区域信号比较

在重力场全球周年振幅模式中,以亚马逊流域的信号最为典型(信号季节性变化幅度最大);在重力场全球长期趋势模式中,以格陵兰地区的信号最为典型(或南极区域,冰川消融趋势信号最大);而在撒哈拉沙漠地区,信号十分平稳,几乎没有质量变化,故常常以撒哈拉沙漠地区信号的RMS来衡量重力场解算的误差水平.因此,本文以这三个典型区域在2005—2010年的质量变化信号来反映各重力场模型的信噪比.需要注意的是,获取上述地区的质量时变信号均采用以下后处理过程:SLR的C20结果代替GRACE解算重力场的C20(Cheng et al.,2004);重力场一阶地心改正项采用(Swenson et al.,2008)建议;500 km高斯平滑滤波;所有球谐系数转化为等效水高(Wahr et al.,1998);所有区域信号的提取均在频域中完成(Swenson et al.,2002).从图 7结果可以看出,Hust-IGG01模型和CSR、GFZ、JPL在亚马逊流域和格陵兰岛的对应质量变化曲线基本重合,而Hust-IGG02的曲线稍微偏离,总体呈现变化幅度偏小;而在撒哈拉沙漠地区([10° E,20° N],[20° N,30° N]),这5个模型所反映出的误差水平基本在同等水平,难以区分.

图 7 基于5个时变模型导出的区域质量变化(以平均等效水高表示(EWH)) (a,b,c)分别表示亚马逊、 格陵兰、 撒哈拉沙漠地区. Fig. 7 Mass variations derived from five different gravity products (a,b,c)represent the mass changes over Amazon,Greenland and Sahara area in terms of EWT,respectively.

为了进一步量化时变信号水平,本文给出了各模型在格陵兰岛、亚马逊流域和撒哈拉沙漠3个典型区域不同信号模式的详细数据统计,并与国际主流权威重力场模型进行对比.这里参与对比的模型主要包含ICGEM(International Centre for Global Earth Models)自2014年以来最新发布的7个无约束模型GFZ RL05a,CSR RL05,JPL RL05,Tongji-GRACE01,Tongji-GRACE02,ITSG2014 以及 AIUB RL02.以上模型所得到的具体结果参见表 4图 8,需要注意的是,在图 8中,为了直观的显示各模型之间的差异程度,所有模型结果均除以CSR RL05的结果来计算其比率.由表 4图 9可以看出,Hust-IGG01模型相比Hust-IGG02模型有很大改善:在周年振幅以及长期趋势两种模式上,与GFZ RL05a、CSR RL05、JPL RL05、ITSG、AIUB相比,Hust-IGG01模型完全达到同等水平,几家结果的互差完全在精度水平之内,而基于简动力学轨道的Hust-IGG02,Tongji01和Tongji02模型总体来说相对于其他模型存在大概5%~10%的信号低估;在撒哈拉沙漠地区,可以看出,以上模型的误差基本都在同等水平,虽然基于简动力学轨道的Hust-IGG02的RMS相比基于运动学轨道的Hust-IGG01的RMS下降了9%,但根据三家官方机构的RMS计算值来看(存在较大的波动),作者认为本文所计算的两个模型仍在同一误差水平.

表 4 各模型在典型区域质量变化的数学统计 Table 4 Statistics of the induced mass variations over typical regions from the following models
图 8表 4的数据与CSR模型的计算比率后的统计图 Fig. 8 The corresponding bar graph after a unit normalization applied to data in Table 4
图 9 2005—2010年全球质量异常模式分析 (a,b)分别表示CSR RL05和Hust-IGG01质量异常长期趋势(cm·a-1);(c,d)分别表示CSR RL05和Hust-IGG01质量异常周年振幅(cm) Fig. 9 This figure demonstrates the different modes of global mass anomaly from 2005 to 2010 in terms of EWH (a,b)stand for the yearly trends of mass anomaly derived from CSR and Hust-IGG01,respectively;(c,d)stand for the annual amplitudes of mass anomaly derived from CSR and Hust-IGG01,respectively.

特别注意的是,由于AIUB模型的计算策略和本文Hust-IGG01基本相同(用运动学轨道代替GPS观测值参与解算),因此所得到的时变信号与Hust-IGG01也是几乎完全相同(三个典型区域的时变信号不符度小于1%),这也反映了本文建立Hust-IGG01模型的数据处理流程的正确性.同时,基于简动力学轨道的模型Hust-IGG02在这些典型区域信号水平的偏弱,说明了简动力学轨道会在一定程度造成重力场时变信号的低估,与运动学轨道或GPS数据反演结果存在差异.

4 Hust-IGG01模型在全球32个流域的实际应用

经过以上讨论,可以看出由运动学轨道结合KBRR数据获取的模型Hust-IGG01精度得到了明显提升.为了进一步验证Hust-IGG01可靠性,在图 9中给出了其在全球尺度上质量变化的长期趋势模式和周年振幅模式,以等效水高表示,并与CSR RL05结果进行对比.可以发现,无论在空间分布上,还是在信号幅度上,Hust-IGG01与CSR RL05几乎没有差别,达到了国际主流机构的处理精度水平.

同时,GRACE时变重力场在实际应用中常常用于分析局部流域的质量变化信号,如长江流域、黄河流域等等,对这些区域的分析可以反映解算模型的精细程度.为了验证Hust-IGG01模型在局部细节上的实际应用价值,本文对全球32个主要流域进行了分析(见图 10,包含所有待分析流域的名称以及其边界分布).在图 11的统计结果中给出了这32个主要流域的质量变化的三种模式: 年际趋势、周年振幅和半周年振幅(均以区域平均等效水高表示).通过比较不同流域不同的模式可以看出,Hust-IGG01得到的结果和CSR RL05基本一致.以上分析表明,Hust-IGG01模型在全球尺度和流域尺度的重力场精度水平与国际主流机构相当.

图 10 全球32个主要流域 Fig. 10 (1 )Amur (2 )Aral (3 )Brahmaputra (4 )Caspian (5 )Colorado (6 )Congo (7 )Danube (8 )Dnieper (9 )Euphrates (10 )Eyre (11 )Ganges (12 )Indus (13 )Lena (14 )Mackenzie (15 )Mekong (16 )Mississippi (17 )Murray (18 )Nelson (19 )Niger (20 )Nile (21 )Ob (22 )Okavango (23 )Orange (24 )Orinoco (25 )Parana (26 )St.Lawrence (27 )Tocantins (28 )Yangtze (29 )Yellow (30 )Yenisey (31 )Yukon (32 )Zambeze
图 11 全球32个流域的平均等效水水高模式分析((a)年际趋势;(b)周年振幅;(c)半周年振幅) Fig. 11 Yearly trends,annual amplitudes and semi-annual amplitudes of the mass variations in terms of EWH over 32 major basins
5 结论

本文首先分析了运动学轨道和简动力学轨道数据的时间序列特征和频谱特征,初步阐释了运动学轨道在重力场信号反演的优势;然后基于GRACE实测数据,我们分别用简动力学轨道和运动学轨道代替GPS数据作为伪观测值,并联合K波段星间测距观测值,解算了60阶Hust-IGG01和Hust-IGG02月解时变重力场模型.初步对比发现,基于运动学轨道解算的模型Hust-IGG01在2005—2010年平均时变重力场系数与CSR RL05结果的相关性高达0.994,明显优于Hust-IGG02;为了分析不同轨道数据对重力场低阶项的影响,本文对各模型C20C60C70C80以及C90时间序列进行分析,发现在各个数学统计指标上Hust-IGG01均略微优于Hust-IGG02,更加接近SLR、AIUB和CSR RL05等权威模型结果;同时为了验证轨道数据控制误差的效果,本文计算了两个模型的formal error,结果发现Hust-IGG01与GFZ RL05、AIUB在同等水平,而Hust-IGG02 的formal error估计过于乐观.

Hust-IGG01推导出的格陵兰岛的长期冰川消融趋势、亚马逊流域的质量变化周年振幅以及撒哈拉沙漠的质量变化RMS,与目前各大权威机构发布的模型的差距均在精度范围之内.除了这些典型区域之外,本文还研究了全球32个主要流域的质量变化模式:年际趋势、周年振幅和半周年振幅,统计结果均显示Hust-IGG01和CSR RL05在同等水平.Hust-IGG02模型在格陵兰地区和亚马逊流域的质量变化信号,虽然与同样采用简动力学轨道作为伪观测值的Tongji01、Tongji02模型结果类似,但相对其他主流模型(基于运动学轨道或者GPS数据)仍存在一定程度的“信号低估”.

尽管造成“信号低估”可能有多方面因素决定(如国内各机构处理数据和解算方法均有差异),但是结合本文试验,我们认为简动力学轨道代替GPS作为伪观测值可能是造成国内解算“信号低估”的一个潜在原因.其根本在于简动力学轨道在GPS定轨阶段引入了一定程度的约束信息,使得解算结果更加偏向先验背景模型,从而压制了实际信号而造成“信号低估”,然而运动学轨道可以一定程度避免这个问题而得到更加真实的重力场时变信号.通过全面的对比,本文利用运动学轨道联合K波段所解算的模型Hust-IGG01完全达到了官方机构的信号水平,因此本文建议在当前GRACE观测值的精度水平下,可以采用运动学轨道来代替GPS数据作为伪观测值进行解算,不仅可以降低计算负荷,也可以达到与直接采取GPS数据类似的效果而没有“信号低估”.本文试验结果同时也说明,在采取合理的数据处理策略的前提下,将GPS定轨过程和重力场反演过程分开并不会明显影响到当前GRACE时变重力场的解算精度,而且比起一步法(直接利用GPS观测数据)同时解算轨道和重力场,能有效的降低计算难度,增加计算效率.

致谢

感谢JPL和GFZ提供的GRACE L1B数据,以及德国格拉茨大学提供的运动学轨道数据,感谢国家留学基金委CSC对本人在德国所有研究工作的资助,特别感谢匿名评审专家对本文文字、图件以及其他相关内容提出的宝贵的修改意见,特别感谢中国科学院测量与地球物理研究所冯伟博士对文章结构组织提出的宝贵意见.

参考文献
Bettadpur S. 2007. CSR Level-2 processing standards document for product release 04 GRACE. The GRACE Project. Center for Space Research, University of Texas at Austin, 327-742.
Bettadpur S. 2009. Recommendation for a-priori bias and scale parameters for Level-1B ACC data (version 2). GRACE tn-02.
Bruinsma S, Lemoine J, Biancale R, et al. 2010. CNES/GRGS 10-day gravity field models (release 2) and their evaluation. Adv. Space Res., 45(4): 587-601. DOI:10.1016/j.asr.2009.10.012
Case K, Kruizinga G, Wu S. 2010. GRACE level 1B data product user handbook. JPL Publication D-22027.
Chambers D P, Wahr J, Nerem R S. 2004. Preliminary observations of global ocean mass variations with GRACE. Geophysical Research Letters, 31(13): L13310.
Chen Q J, Shen Y Z, Chen W, et al. 2015a. A modified acceleration-based monthly gravity field solution from GRACE data. Geophysical Journal International, 202(2): 1190-1206. DOI:10.1093/gji/ggv220
Chen Q J, Shen Y Z, Zhang X F, et al. 2015b. Monthly gravity field models derived from GRACE Level 1B data using a modified short-arc approach. Journal of Geophysical Research:Solid Earth, 120(3): 1804-1819. DOI:10.1002/2014JB011470
Cheng M K, Tapley B D. 2004. Variations in the Earth's oblateness during the past 28 years. Journal of Geophysical Research:Solid Earth, 109(B9): B09402.
Dunn C, Bertiger W, Bar-Sever Y, et al. 2003. Instrument of GRACE:GPS augments gravity measurements. GPS World, 14(2): 16-28.
Flechtner F, Dahle C, Gruber C, et al. 2013. Status GFZ RL05 and RL05a GRACE L2 products.//Paper presented at the GRACE Science Team Meeting. Austin, TX.
Flechtner F, Dobslaw H. 2013. AOD1B product description document for product release 05. GFZ German Research Centre for Geosciences.
Han B M, Xu H Z. 2007. GPS-based reduced-dynamic orbit determination for gravity field recovery. Progress in Geophysics, 22(1): 73-79.
Han S C, Simons F J. 2008. Spatiospectral localization of global geopotential fields from the gravity recovery and climate experiment (GRACE) reveals the coseismic gravity change owing to the 2004 Sumatra-Andaman earthquake. Journal of Geophysical Research:Solid Earth, 113(B1).
Jäggi A, Beutler G, Bock H, et al. 2007. Kinematic and highly reduced-dynamic LEO orbit determination for gravity field estimation.//Dynamic Planet. Berlin Heidelberg:Springer, 354-361.
Jäggi A, Beutler G, Meyer U, et al. 2012. AIUB-GRACE02S:status of GRACE gravity field recovery using the celestial mechanics approach.//Geodesy for Planet Earth. Berlin Heidelberg:Springer, 161-169.
Kusche J, Schmidt R, Petrovic S, et al. 2009. Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model. Journal of Geodesy, 83(10): 903-913. DOI:10.1007/s00190-009-0308-3
Liu X L. 2008. Global gravity field recovery from satellite-to-satellite tracking data with the acceleration approach[Ph. D. thesis]. Dutch:Delft University of Technology.
Luthcke S B, Rowlands D D, Lemoine F G, et al. 2006. Monthly spherical harmonic gravity field solutions determined from GRACE inter-satellite range-rate data alone. Geophysical Research Letters, 33(2).
Mayer-Gürr T. 2006. Gravitationsfeldbestimmung aus der Analyse kurzer Bahnbögen am Beispiel der Satellitenmissionen CHAMP und GRACE[Ph. D. thesis]. Bonn:Universitäts und Landesbibliothek Bonn.
Mayer-Gürr T, Zehentner N, Klinger B, et al. 2014. ITSG-Grace 2014:a new GRACE gravity field release computed in Graz.//Proceedings of GRACE Science Team Meeting. Potsdam.
Meyer U, Dahle C, Sneeuw N, et al. 2016. The effect of pseudo-stochastic orbit parameters on GRACE monthly gravity fields:Insights from lumped coefficients.//VIII Hotine-Marussi Symposium on Mathematical Geodesy. Berlin, Heidelberg:Springer, 177-183.
Petit G, Luzum B. 2010. IERS conventions (2010). DTIC Document.
Ran J J, Xu H Z, Zhong M, et al. 2014. Global temporal gravity filed recovery using GRACE data. Chinese Journal of Geophysics (in Chinese), 57(4): 1032-1040. DOI:10.6038/cjg20140402
Richey A S, Thomas B F, Lo M H, et al. 2015. Quantifying renewable groundwater stress with GRACE. Water Resources Research, 51(7): 5217-5238. DOI:10.1002/2015WR017349
Ries J C, Bettadpur S, Poole S, et al. 2011. Mean background gravity processing mean gravity field from space and ground.//Paper presented at the GRACE Science Team Meeting.
Rieser D, Mayer-Gürr T, Savcenko R, et al. 2012. The ocean tide model EOT11a in spherical harmonics representation. Technical Note.
Schwintzer P, Reigber C. 2002. The contribution of GPS flight receivers to global gravity field recovery. Journal of Global Positioning Systems, 1(1): 61-63. DOI:10.5081/jgps
Shang K. 2015. GRACE time-variable gravity field recovery using an improved energy balance formalism[Ph. D. thesis]. Columbus:The Ohio State University.
švehla D, Rothacher M. 2005. Kinematic precise orbit determination for gravity field determination.//A Window on the Future of Geodesy. Berlin Heidelberg:Springer, 181-188.
Swenson S, Wahr J. 2002. Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment (GRACE) measurements of time-variable gravity. Journal of Geophysical Research:Solid Earth (1978-2012), 107(B9): ETG 3-1-ETG 3-13. DOI:10.1029/2001JB000576
Swenson S, Chambers D, Wahr J. 2008. Estimating geocenter variations from a combination of GRACE and ocean model output. Journal of Geophysical Research:Solid Earth, 113(B8): B08410.
Tapley B D, Bettadpur S, Watkins M, et al. 2004. The gravity recovery and climate experiment:Mission overview and early results. Geophysical Research Letters, 31(9).
Touboul P, Willemenot E, Foulon B, et al. 1999. Accelerometers for CHAMP, GRACE and GOCE space missions:synergy and evolution. Boll. Geof. Teor. Appl., 40(3-4): 321-327.
Wagner C A, McAdoo D C. 2012. Error calibration of geopotential harmonics in recent and past gravitational fields. Journal of Geodesy, 86(2): 99-108. DOI:10.1007/s00190-011-0494-7
Wahr J, Molenaar M, Bryan F. 1998. Time variability of the Earth's gravity field:Hydrological and oceanic effects and their possible detection using GRACE. Journal of Geophysical Research, 103(B12): 30205-30229. DOI:10.1029/98JB02844
Wang C Q, Xu H Z, Zhong M, et al. 2015. Monthly gravity field recovery from GRACE orbits and K-band measurements using variational equations approach. Geodesy and Geodynamics, 6(4): 253-260. DOI:10.1016/j.geog.2015.05.010
Wang C Q, Xu H Z, Zhong M, et al. 2015. An investigation on GRACE temporal gravity field recovery using the dynamic approach. Chinese J. Geophys. (in Chinese), 58(3): 756-766. DOI:10.6038/cjg20150306
Zehentner N, Mayer-Guerr T. 2013. Kinematic orbits for GRACE and GOCE based on raw GPS observations.//Poster presented at the IAG Scientific Assembly. Potsdam, Germany.
Zhou H, Luo Z C, Zhong B. 2015. WHU-Grace01s:A new temporal gravity field model recovered from GRACE KBRR data alone. Geodesy and Geodynamics, 6(5): 316-323. DOI:10.1016/j.geog.2015.07.004
冉将军, 许厚泽, 钟敏, 等. 2014. 利用GRACE重力卫星观测数据反演全球时变地球重力场模型. 地球物理学报, 57(4): 1032–1040. DOI:10.6038/cjg20140402
王长青, 许厚泽, 钟敏, 等. 2015. 利用动力学方法解算GRACE时变重力场研究. 地球物理学报, 58(3): 756–766. DOI:10.6038/cjg20150306