地球物理学报  2015, Vol. 58 Issue (3): 756-766   PDF    
利用动力学方法解算GRACE时变重力场研究
王长青1,2, 许厚泽1, 钟敏1, 冯伟1, 冉将军1, 杨帆3    
1. 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室, 武汉 430077;
2. 中国科学院大学, 北京 100049;
3. 华中科技大学物理学院, 武汉 430074
摘要:本文利用动力学方法建立GRACE (Gravity Recovery And Climate Experiment) K波段距离变率 (KBRR)观测、轨道观测与重力场系数的观测方程, 通过GRACE Level 1B观测数据,成功解算出全球月时变重力场模型——IGG时变重力场模型,并将2008—2009年的解算结果与GRACE三大数据处理机构美国德克萨斯大学空间中心CSR (Center for Space Research)、美国宇航局喷气推进实验室JPL(Jet Propulsion Laboratory)和德国地学研究中心GFZ (GeoForschungs Zentrum) 发布的最新全球时变重力场模型进行详细对比分析.结果表明:IGG结果在全球质量异常、中国及周边地区质量异常的趋势变化、全球质量异常均方差、2~60每阶位系数差值以及亚马逊流域和撒哈拉沙漠等典型区域平均质量异常等方面与CSR、JPL和GFZ解算的RL05结果较为一致.其中,IGG解算结果在2~20阶与CSR、GFZ和JPL最新解算结果基本一致,20~40阶IGG解算结果与GFZ、JPL单位最新解算结果较为接近,大于40阶IGG结果介于CSR与GFZ、JPL之间;亚马逊流域平均质量异常周年振幅IGG、CSR、GFZ和JPL获取到的结果分别为17.6±1.1 cm、18.9±1.2 cm、17.8±0.9 cm和18.9±1.0 cm等效水柱高.利用撒哈拉沙漠地区的平均质量异常做反演精度评定,IGG、CSR、GFZ和JPL的时变重力场获取到的平均质量异常均方差分别为1.1 cm、0.9 cm、0.8 cm和1.2 cm,表明IGG解算结果与CSR、GFZ和JPL最新发布的RL05结果在同一精度水平.
关键词GRACE     时变重力场     动力学方法     等效水柱高    
An investigation on GRACE temporal gravity field recovery using the dynamic approach
WANG Chang-Qing1,2, XU Hou-Ze1, ZHONG Min1, FENG Wei1, RAN Jiang-Jun1, YANG Fan3    
1. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophsics, Chinese Academy of Sciences, Wuhan 430077, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China;
3. School of Physics, Huazhong University of Science and Technology, Wuhan 430074, China
Abstract: The Gravity Recovery and Climate Experiment (GRACE) mission can significantly improve our knowledge of the temporal variability of the Earth's gravity field. We intend to obtain monthly gravity field solutions based on dynamic approach (variational equations approach) from GPS-derived positions of GRACE satellites and K-band range-rate measurements. Moreover, these solutions will validate through GRACE RL05 products published monthly gravity field solutions by the GRACE project, which including CSR (Center for Space Research), JPL (Jet Propulsion Laboratory) and GFZ (GeoForschungsZentrum).
Based on the theory of the dynamic approach, the processing began with the computation of purely dynamic orbits by fitting GNV1B orbits computed by the GRACE team using a reduced dynamic strategy. This step aimed to calibrate the GRACE A and B accelerometer data and the initial state vectors of these two satellites. Note that gravity field model parameters were not adjusted in this part. In the next step, the purely dynamic orbits were used as the reference orbit for the calculation of residual orbits, and were also applied to calculate the nominal value of K-band range-rate. In addition, the partial derivatives with respect to the initial state vectors, accelerometer parameters and geopentional coefficients were computed simultaneously with the integration of the purely dynamic orbits. Finally, normal equations of orbit positions and K-band range-rate measurements were synthesized using a fixed data weighting ratio to combine the both information matrices for the best combined solution.
The monthly gravity field solution obtained though above procedures was named as the IGG temporal gravity field model. IGG temporal gravity field models were compared with GRACE RL05 products in following aspects: (i) monthly gravity field solutions in February, 2008 and September, 2009; (ii) the root mean squares of the global mass anomaly during 2008 to 2009; (iii) the trend of the mass anomaly in China and its nearby regions within 2008 to 2009; (iv) the square root degree difference variance with respect to the mean gravity field GIF48 model from 2008 to 2009; (v) time-series changes in the mean water storage in the region of the Amazon Basin and the Sahara Desert between 2008 and 2009.
The results showed that IGG solutions were almost consistent with GRACE RL05 products in above aspects (i)—(iii). For aspect (iv), IGG solutions were in good agreement with latest published GRACE RL05 products at the degree from 2 to 20, and were at the same accuracy level with GFZ and JPL at the degree from 20 to 40, but were between CSR and GFZ (or JPL) at the degree exceeding 40. According to aspect (v), changes in the annual amplitude of mean water storage in the Amazon Basin were 17.6±1.1 cm for IGG, 18.9±1.2 cm for CSR, 17.8±0.9 cm for GFZ and 18.9±1.0 cm for JPL in terms of equivalent water height, respectively. The root mean squares of the mean mass anomaly in Sahara were 1.1 cm, 0.9 cm, 0.8 cm and 1.2 cm for temporal gravity field models of IGG, CSR, GFZ and JPL, respectively. Overall, it was noticeable that IGG temporal gravity field solutions were at the same accuracy level with the latest temporal gravity field solutions published by CSR, GFZ and JPL.
Key words: GRACE     Temporal gravity field     The dynamic approach     Equivalent water height    
 1 引言

GRACE(Gravity Recovery and Climate Experiment)是以空间分辨率400~40000 km,测量全球重力场为目标的卫星任务(Tapley et al., 2005),从GRACE时变重力场解可以获取月尺度至多年尺度 地球表面流体质量运移的时变特性.GRACE任务由NASA(National Aeronautics and Space Administration)和DLR(Deutsches Zentrum fur Luft-und Raumfahrt)共同实施,发射于2002年,包括两颗近圆轨道约为500 km、轨道倾角约89°相互分开约220 km的同轨卫星.GRACE装载高精度的K波段微波系统,其距离及距离变率测量精度优于10 μm和0.1 μm/s,有高精度GPS接收机、激光反射器、恒星敏感器、以及具有测量精 度1.0×10-10 m·s-2非保守力加速度计(Bruinsma et al., 2010).

基于GRACE观测数据,国内外的研究机构得到了一系列的时变重力场模型.国际上如GRACE项目团队,即美国德克萨斯大学奥斯汀分校空间研 究中心CSR(Center of Space Research)(Bettadpur,2012)、 德国地学研究中心GFZ(Geo Forschungs Zentrum)(Dahle et al., 2013)、美国宇航局喷气推 进实验室JPL(Jet Propulsion Laboratory)(Watkins,2012)、法国国家空间局等采用动力学方法(Bruinsma et al., 2010),瑞士伯尔尼大学采用天体力学方法(也可称为动力学方法)(Jggi, et al., 2012),荷兰代尔夫特理工大学采用加速度方法(Liu et al., 2010),德国波恩大学采用短弧长法(Mayer-Gürr,2006)等成功获取时变重力场模型;国内相关同行利用动力学方法、能量法、短弧长方法等在卫星重力反 演领域也做了大量的工作(罗佳,2003周旭华,2005王正涛,2005张兴福,2007邹贤才,2007郑伟,2007游为,2011Zhao et al., 2011冉将军,2014).其中,Zhao等(2011)利用动力学方法,冉将军等(2014)利用短弧长方法,Shen等(2013)利用改进的短弧长方法 成功获取了和国际上知名机构精度相当的时变重力场模型.

高精度时变重力场的获取是卫星重力技术中的难点,也是大地测量领域的研究热点.GRACE任务的实施使得时变重力场研究获取突破性的进展,我国在卫星重力场获取方面的研究尽管取得突破,但是由于相关研究的开展晚于国际同行,因此,相关的研究亦相对落后于国际同行.而近年来,我国的重力卫星任务逐渐提上日程,高精度时变重力场获取研究是目前我国卫星重力技术亟需解决的难题.在短弧长方法的研究方面,冉将军等(2014)Shen等(2013)利用最新发布的GRACE数据获取了与GRACE RL05精度相当的时变重力场模型;在动力学方法研究方面,继Zhao等(2011)之后,利用最新发布的GRACE数据解算时变重力场,并将解算结果与GRACE RL05时变重力场进行对比等仍缺乏后续研究.鉴于此,本文利用动力学两步法,采用最新发布的GRACE数据进行了时变重力场解算研究.GRACE观测数据中轨道信息对于低阶重力场较敏感,KBR测距(K-b and ranging,KBR)或距离变率(K-b and inter-satellite range-rate,KBRR)对于地球重力场中高阶较为敏感,因此,利用动力学方法采用GRACE观测数据获取地球时变重力场一般是结合轨道信息和KBR测距或距离变率信息的方式进行时变重力场的获取(Kim,2000Bruinsma et al., 2010Zhao et al., 2011Jggi, et al., 2012),如GRACE的RL05全球重力场解(CSR RL05、JPL RL05以及GFZ RL05).本文在重力场的解算过程中根据观测数据先验权,结合KBRR观测数据与轨道观测建立重力场系数的观测方程,成功地提取到全球时变重力场信号,同时将我们的解算结果与GRACE的RL05全球重力场解算结果进行了全面的对比分析. 2 数据处理

本文计算根据动力学两步法原理采用自主研制软件进行重力场的解算,关于轨道和星间距离变率反演重力场的基本观测模型具体参见(冉将军等,2012),在此不再赘述.数据处理步骤归结为:(1)通过拟合轨道观测计算动力学轨道;(2)使用轨道残差和星间距离变率残差求解无约束重力场解.轨道拟合之目的在于通过已有力学模型,调整加速计参数来逼近真实轨道,剩余的轨道残差看作是由重力场时变信号和观测噪声引起,同时已经调整好的轨道用作KBRR观测参考值的计算,以获取包含重力场时变信号以及观测噪声的KBRR观测残差.如前所述,轨道观测数据以及KBRR观测数据都包含重力场时变信号,因此,在利用GRACE观测数据进行时变重力场的解算时我们结合轨道观测和KBRR距离变率观测.轨道观测数据、KBRR观测数据属于不同类型、不同精度的观测数据,对于此问题假设已有残差观测方程如下:

公式(1)中,v i为观测残差,A i表示设计矩阵,d i观测值减参考值,w i表示观测权,σi表示观测精度(i=1,轨道观测;i=2,KBRR观测),I 为单位矩阵,δ x 为两类观测共同的求解参数向量,将上式整理后则有

此式即为我们在重力场解算中结合轨道数据和KBRR数据时采用的基本公式. 2.1 轨道拟合

通过JPL提供的精密的GNV1B轨道产品校正GRACE A 、GRACE B 加速度计数据偏差和初始 状态参数,采用数据包括GNV1B、SCA1B、ACC1B、 KBR1B、AOD1B等最新公布的GRACE Level 1B的产品(Case et al,2010Flechtner and Dobslaw, 2013).轨道运动方程采用以龙格库塔起步的8阶GAUSS-JACSON(Berry and Healy, 2004)积分器算法,以5 s为步长,轨道积分弧长选择6 h进行积分,积分过程中采用的力模型如表 1.需要说明的 是:GRACE加速度计数据被认为是1.0×10-10 m·s-2 的测量精度,重力场解算中海潮模型的18个主潮波是不够的,需要计算更多的潮波分量(Rieser et al., 2012);使用AOD1B数据时利用大汽潮模型(Biancale and Bode, 2006)扣除了大气潮S2的影响,大气潮S2的贡献在轨道积分过程中予以考虑,N体摄动中行星J2项效应也予以考虑,参见表 1.校正参数包括卫星初始状态参数(6个参数)、加速度计数据在SRF坐标框架下的每个方向偏差(3个参数),初始轨道状态参数的初始值由GNV1B获取,加速度计参数的初始值参照Bettadpur(2009)计算获得,利用校正后的参数重新进行轨道积分得到动力学轨道与GNV1B产品的残差RMS小于2 cm.以2008年2月为例给出该月份惯性坐标系下三个方向残差RMS序列(6 h拟合一次,总共有116个弧段),如图 1.

图 1 GRACE A 、B卫星拟合轨道与GNV1B轨道的残差的RMS,X、Y、Z分别表示在惯性系下的3个方向Fig. 1 Root mean squares of residuals between GRACE A,B satellite fitted orbit and GNV1B products. X,Y,Z represent three directions in inertial coordinate system respectively
2.2 残差计算以及重力场系数求解

将上一步中获取到的动力学轨道做参考轨道,以GNV1B和KBRR为观测值进行轨道残差、星间距离变率残差的计算以及粗差探测,同时计算观测值对初始状态参数、加速度计数据偏差以及重力场系数偏微分.轨道残差和星间距离变率残差获取后,求解参数(包括初始状态参数、加速度计偏差以及重力场系数)的观测方程可以通过以下式子建立:

Δ r i为卫星(i=A或B)轨道残差、为星间距离变率残差,为星间距离变率参考值,r i0,bi,β 表示卫星(i=A或B)的初始状态参数、加速度偏差参数以及重力场系数,在轨道运动方程积分的过程中求出,为星间距离变率对初始状态参数、加速度偏差参数以及重力场系数的偏导数,该偏导数由 线性组合获得(Kim,2000).

将轨道观测通过公式(3)整理后,GRACE A、GRACE B卫星的轨道观测可以形成法方程

xL、xG分别表示局部参数(GRACE A、B卫星的初始状态参数以及加速度偏差)和全局参数(重力场系数),N 11,N 22是对应于局部参数和全局参数的子矩阵,l xL,l xG为法方程右向量.类似地,根据公式(4)每个弧段KBRR观测可组成法方程

根据前述的结合轨道数据和KBRR数据的方法,有以下关系:

其中,

σkbrr为KBRR观测精度,σleo为轨道观测精度,KBRR残差RMS位于0.2~0.4 μm·s-1之间(见图 2,以2009年9月份计算得到的KBRR的RMS统计结果为例),在计算过程中我们取KBRR观测精度为0.2 μm·s-1,轨道观测精度为2 cm.公式(7)即为该弧段结合轨道数据和KBRR数据的法方程,消去法方程的局部参数可得到关于全局参数的法方程为

将上式记为

i表示为第i个观测弧段对应的法方程,将每个月的n个弧段法方程叠加则有

此式即为每个月重力场系数解算的法方程.
图 2 2009年9月KBRR残差的RMS(弧长6 h)Fig. 2 Root mean squares of KBRR residuals in Sep,2009(Arc length 6 hours)

为验证解算策略的可行性,利用已经发布的2008年1月份CSR RL05模型作为真实模型,以GIF48作为参考模型加入2 cm轨道误差、1.0×10-4 m·s-1的速度误差、1.0×10-7 m·s-1的KBRR观测误差,不考虑其他误差,对真实模型进行重建.图 3a为重建重力场模型和真实模型的每阶大地水准面高(实线)和累计大地水准面高(虚线)对比,图 3(b、c)为以等效水柱高表示的真实模型与重建模型的质量异常的空间分布.在上述误差设置下,真实模型与重建模型的累计大地水准面差距60阶时相差约为0.04 mm,详见图 3a,质量异常的全球分布和强度相当,如图 3(b、c);上述误差设置造成的质量异常的误差也即真实模型与重建模型质量异常差异的标准偏差约为1.92 cm,最大(小)误差约 5.97(-5.6)cm.通过图 3的综合对比可以看出在该误差设置下能够基本恢复出原来的信号.

图 3 重建重力场模型和真实模型的比较
(a)大地水准面高,TRUE 代表真实模型,MODEL代表重建的模型,实线表示每阶大地水准面高,虚线表示累计大地水准面高;(b)真实模型质量异常;(c)重建模型质量异常.单位cm,500 km平滑滤波.
Fig. 3 Comparison between constructed model and true model
(a)Geoid height,TRUE-true model,MODEL-constructed model,solid line represents geoid height per degree,dash line represents accumulative geoid height;(b)Mass anomaly from true model;(c)Mass anomaly from constructed model. Equivalent water height in cent-meter,filtering radius 500 km.

GRACE项目组发布的RL05时变重力场模型使用新的GRACE Level 1B产品(L1B_V2)、新的静态重力场模型、海潮模型以及去混频产品等,使之比时变重力场模型RL04具有更高的精度.为利用实测数据对处理策略进行验证,我们分别采用不同版本的GRACE L Level 1B产品(L1B_V1、L1B_V2)计算了2009年9月时变重力场模型,记为IGG_V01、IGG_V02,并将解算结果与CSR RL04、RL05时变重力场模型进行对比(IGG_V01解算时采用的力模型包括:海潮模型采用FES2004,静态重力场模型采用GIF22a,AOD1B数据采用RL04版本,IERS2003规范计算固体潮及固体极潮,其他如上述),如图 4 .IGG_V01、IGG_V02分别和CSR RL04、RL05处在同一精度水平,同时IGG_V02和CSR RL05相比IGG_V01、CSR RL04在30阶以上精度至少提高了1倍以上.通过图 4对比说明新的力学模型尤其是海潮模型、静态重力场模型以及新的GRACE Level 1B产品等能够很大程度上提高时变重力场解算精度(尤其体现在中高阶,如15阶以上).通过上述模拟数据和实测数据解算说明我们的处理策略具有恢复时变重力场信号的能力.

图 4 2009年9月CSR RL04、CSR RL05、 IGG_V01、IGG_V02 时变重力场模型每阶大地水准面差距Fig. 4 Geoid height per degree of CSR RL04,CSR RL05,IGG_V01,IGG_V02 temporal gravity field in Sep. 2009
3 解算结果分析

本文使用上述处理方法,利用最新的GRACE Level 1B产品(L1B_V2)以及力学模型(如表 1)解算了截至60阶的无约束时变重力场(以IGG代表 本文解算结果,下同).为评价解算精度,对2008— 2009年的解算结果,在单月时变重力场解、2008—2009年全球质量变化的均方差、每阶位系数变化以 及典型区域平均质量异常等方面与国际知名机构 CSR、GFZ和JPL等最新的解算结果RL05进行了详细对比分析.由于受GRACE卫星自身轨道设计与观测精度的制约,无约束时变重力场球谐系数存在较大的噪声,在空间分布上还存在南北方向的“条带”误差.我们对CSR、GFZ、JPL和IGG解算的GRACE结果采取了以下后处理策略:用卫星激光 测距(SLR)结果替换了GRACE的 C(2,0)项(Cheng and Tapley, 2004),一阶地心改正项采用Swenson等(2008)的结果,并进行了去条带和高斯平滑处理(Jekeli,1981; Swenson and Wahr, 2006),地壳均衡调整(GIA)导致的时变重力场长期变化效应也进行了扣除(Geru et al., 2013),球谐系数最终转化为全球分布的等效水柱高变化(Wahr et al., 1998).

表 1 重力场解算所用的力模型 Table 1 Force models used for gravity field determination

首先,本文以2008年2月和2009年9月解算的全球时变重力场模型为例,与CSR、GFZ和JPL等单位最新公布的RL05结果对比如下.一般认为30阶以下包含了较强的重力场信号,40阶以上主要反映噪声水平(Liu,2008),由图 5中每阶大地水准面可以看出,2008年2月以及2009年9月IGG解算结果与CSR最新发布的解算结果精度相当.图 6—7说明四单位2008年2和2009年9月全球质量异常空间分布以及信号强度基本上能够保持一致.

图 5 IGG、CSR、 GFZ和JPL每阶大地水准面差距
(a)2008年2月;(b)2009年9月.
Fig. 5 Geoid height per degree of IGG,CSR,GFZ and JPL
(a)Feb. 2008;(b)Sep. 2009.

图 6 2008年2月由GRACE时变重力场模型获取的质量异常
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. 等效水柱高单位为cm,滤波半径500 km.
Fig. 6 Mass anomaly from GRACE temporal gravity model in February 2008
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. Equivalent water height in centimeter,filtering radius 500 km.

图 7 2009年9月由GRACE时变重力场模型获取的质量异常
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. 等效水柱高单位cm,滤波半径500 km.
Fig. 7 Mass anomaly from GRACE temporal gravity model in February 2008
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. Equivalent water height in centimeter,filtering radius 500 km.

为进一步评价IGG时变重力场解的总体解算精度,本文利用四单位2008—2009年的解算结果计算了中国及周边地区质量异常的线性变化、全球质量异常均方差和这两年时间段内2~60阶的每阶位系数差值的平均值.中国及周边地区质量异常的2008—2009年线性变化信号如图 8所示,IGG、 CSR、GFZ和JPL等四单位的GRACE解算结果 均观测到印度北部和我国华北地区由于地下水超采导致的质量亏损信号(Feng et al., 2013; Rodell et al., 2009; Tiwari et al., 2009),以及横断山脉地区 山岳冰川融化导致的质量亏损信号(Yi and Sun, 2014).利用CSR、GFZ和JPL等最新发布的结果与IGG解算结果计算得到全球质量异常的均方差,如图 9所示.由图 9可知,不同机构的GRACE时变重力场模型计算得到的质量异常信号强度与空间分布有很好的一致性,在南美的亚马逊河流域、非洲的尼日尔河流域和东南亚的季风区等质量变化的强度较大.除空间域外,我们利用每阶位系数差值在频谱域进一步比较了不同机构的时变重力场模型的精度(Kim,2000游为,2011),我们以2008—2009年2~60阶每阶位系数差值的平均值作为四单位2008—2009年时变重力场模型的精度评价依据,如图 10所示.由图 10可知,IGG解算结果在2~20阶左右与其他三单位结果非常接近,20~40阶IGG解算结果与JPL、GFZ两单位的较为接近,在大于40阶部分,IGG解算结果介于CSR与GFZ、JPL最新发布解算结果RL05之间.

图 8 中国及周边地区 2008—2009年质量异常的变化趋势
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. 等效水柱高单位为cm,滤波半径500 km.
Fig. 8 Trend of mass anomaly in China and nearby regions over 2008—2009
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. Equivalent water height in centimeter,filtering radius 500 km.

图 9 2008—2009年24个月全球质量异常的均方差
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. 等效水柱高,单位cm,滤波半径500 km.
Fig. 9 Root mean square of global mass anomaly over 2008—2009
(a)IGG;(b)CSR;(c)GFZ;(d)JPL. Equivalent water height in centimeter,filtering radius 500 km.

图 10 2008—2009年24个月IGG、CSR、GFZ和JPL解算结果的每阶位系数差值的平均值比较 Fig. 10 Comparison of IGG,CSR,GFZ and JPL solution on degree difference variance mean over 2008—2009

此外,本文选取亚马逊整个流域地区和撒哈拉 沙漠地区(10°E—20°E,20°N—30°N),对四单位的时变重力场在上述两个地区的纬度加权平均质量变化信号进行对比.在亚马逊流域,四单位结果都表现出一致的平均质量周年变化,其主要反映了陆地水储量的变化.如图 11所示,IGG得到的2008—2009 年亚马逊流域平均水储量变化的周年振幅为17.6±1.1 cm 等效水柱高,略小于其他三机构的结果(CSR 为18.9±1.2 cm等效水柱高、GFZ为17.8±0.9 cm 等效水柱高和JPL为18.9±1.0 cm等效水柱高);撒哈拉沙漠地区时变重力场信号相对较小,不同机构反演的沙漠地区的时变结果可以用来评估其反演的 精度水平,CSR、GFZ、JPL和IGG四单位反演结果 RMS分别为0.9、0.8、1.2 cm和1.1 cm等效水柱高.IGG结果与其他三单位结果在同一精度水平.

图 11 基于IGG、CSR、JPL和GFZ时变重力场模型获取的2008—2009年亚马逊流域(a)和撒哈拉沙漠地区(b)的平均质量异常信号. 等效水柱高单位为cm,滤波半径500 kmFig. 11 Mean mass anomaly in Amazon(a) and Sahara desert(b)over 2008—2009 based on temporal gravity models from IGG,CSR,JPL, and GFZ. Equivalent water height in centimeter,filtering radius 500 km
4 结论

本文通过动力学两步法、利用模拟数据和实测数据验证解算策略的正确性,然后利用KBRR实测数据并结合轨道信息建立重力场系数的观测方程,解算了时变重力场模型IGG,并将IGG时变重力场模 型与CSR、JPL和GFZ三个国际权威机构的2008—2009年的最新RL05结果进行了详细的对比分析.

通过四单位2008—2009年解算的结果对比表明,IGG结果在全球质量异常空间分布,2008—2009年中国及周边地区质量异常线性变化,2008—2009年全球质量异常均方差空间分布及信号强度,2008—2009年2~60阶的每阶位系数差值的平均值,以及典型区域平均水储量等方面与CSR、JPL和GFZ等基本上保持一致.GRACE卫星重力场反演是复杂的数据处理过程,在精确的背景模型的保障下,还涉及修复间断数据、剔除观测数据粗差等数据预处理和重力场解算过程中参数处理策略的选择.观测数据中粗差对重力场的解算结果会有很大影响,不同的参数处理策略,如不同的尺度因子和加速度计参数的处理方式对于重力场解算结果也会有影响(冉将军,2014).由于不能详细了解CSR、JPL和GFZ等单位详细的数据处理过程,不能了解其参数处理策略,目前我们还无法推断造成IGG解算结果与其他三个单位解算结果差别的原因,这有待于今后进一步的研究.

综上所述,IGG全球时变重力场解算结果精度接近CSR、JPL、GFZ等单位最新发布的全球时变重力场解算结果. 致谢 中国科学院测量与地球物理研究所闫昊明研究员对文章结构以及结果比对方面提出了宝贵意见,德国GFZ提供GRACE Level 1B数据,本人和上海同济大学沈云中教授、中国科学院上海天文台周旭华研究员在重力场解算过程中的相关问题进行了讨论,两位匿名审稿专家对本文提出了修改意见,文中所涉图件由GMT软件绘制,在此一并表示感谢.

参考文献
[1] Berry M M, Healy L M. 2004. Implementation of gauss-jackson integration for orbit propagation. The Journal of the Astronautical Sciences, 52(3): 331-357.
[2] Bettadpur S. 2009. Recommendation for a-priori bias & scale parameters for Level-1B ACC data (version 2). GRACE TN-02 Version2.
[3] Bettadpur S. 2012. UTCSR Level-2 processing stands document for Level-2 product release 005. Center for Space Research, University of Texes at Austin.
[4] Biancale R, Bode A. 2006. Mean annual and seasonal atmospheric tide models based on 3-hourly and 6-hourly ECMWF surface pressure data. Potsdam: Deutsches GeoForschungsZentrum GFZ, Scientific Technical Report STR 06/01, doi: 10.2312/GFZ.b103-06011.
[5] 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.
[6] Case K, Kruizinga G, Wu S. 2010. GRACE Level 1B data product user handbook. Jet Propulsion Laboratory, California Institute of Technology.
[7] Cheng M K, Tapley B D. 2004. Variations in the earth's oblateness during the past 28 years. J. Geophys. Res., 109(B9): B09402, doi: 10.1029/2004JB003028.
[8] Dahle C, Flechtner F, Gruber C, et al. 2013. GFZ GRACE Level-2 processing standards document for Level-2 product release 0005: revised edition. Potsdam: Deutsches GeoForschungsZentrum GFZ, doi: 10.2312/GFZ.b103-1202-25.
[9] Desai S D. 2002. Observing the pole tide with satellite altimetry. J. Geophys. Res., 107(C11): 7-1-7-13, doi: 10.1029/2001JC001224.
[10] Feng W, Zhong M, Lemoine J M, et al. 2013. Evaluation of groundwater depletion in North China using the Gravity Recovery and Climate Experiment (GRACE) data and ground-based measurements. Water Resour. Res., 49(4): 2110-2118, doi: 10.1002/wrcr.20192.
[11] Flechtner F, Dobslaw H. 2013. GRACE AOD1B product description document for product release 05. GFZ German Research Centre for Geosciences.
[12] Geruo A, Wahr J, Zhong S J. 2013. Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada. Geophys. J. Int., 192(2): 557-572, doi: 10.1093/gji/ggs030.
[13] Jäggi A, Beutler G, Meyer U, et al. 2012. AIUB-GRACE02S: Status of GRACE gravity field recovery using the celestial mechanics approach. // Kenyon S, Pacino M C, Marti U eds. Geodesy for Planet Earth, International Association of Geodesy Symposia 136. Berlin Heidelberg: Springer, doi: 10.1007/978-3-642-20338-1_20.
[14] Jekeli C. 1981. Alternative methods to smooth the earth's gravity field. Report No. 327, Reports of the Department of Geodetic Science and Surveying. Columbus: The Ohio State University.
[15] Kim J. 2000. Simulation study of a low-low satellite-to-satellite tracking mission [Ph. D. thesis]. Texas: The University of Texas at Austin.
[16] Liu X. 2008. Global gravity field recovery from satellite-to-satellite tracking data with the acceleration approach [Ph. D. thesis]. Delft: Delft University of Technology.
[17] Liu X, Ditmar P, Siemes C, et al. 2010. DEOS Mass Transport model (DTM-1) based on GRACE satellite data: methodology and validation. Geophys. J. Int., 181(2): 769-788, doi: 10.1111/j.1365-246X.2010.04533.x.
[18] Luo J. 2003. Theory and methodology of earth gravity field determination using satellite-to-satellite tracking [Ph. D. thesis](in Chinese). Wuhan: Wuhan University.
[19] Mayer-Gürr T. 2006. Gravitationsfeldbestimmung aus der Analyse kurzer Bahnboegen am Beispiel der Satellitenmissioncn CHAMP und GRACE [Ph. D. thesis]. Bonn: Institute für Theoretische Geodäsie der Universität Bonn.
[20] Ran J. 2014. Theory, methodology and application of recovery using low-low tracking gravity satellite data[Ph. D. thesis](in Chinese). Wuhan: Institute of Geodesy and Geophysics. Chinese Academy Sciences.
[21] Ran J J, Xu H Z, Shen Y Z, et al. 2012. Expected accuracy of the global gravity field for next GRACE satellite gravity mission. Chinese J. Geophsys. (in Chinese), 55(9): 2898-2908, doi: 10.6038/j.issn.0001-5733.2012.09.009.
[22] Ran J J, Xu H Z, Zhong M, et al. 2014. Global temporal gravity filed recovery using GRACE data. Chinese J. Geophysics. (in Chinese), 57(4): 1032-1040, doi: 10.6038/cjg20140402.
[23] Rieser D, Mayer-Gürr T, Savcenko R, et al. 2012. The ocean tide model EOT11a in spherical harmonics representation. Technical Note.
[24] Rodell M, Velicogna I, Famiglietti J S. 2009. Satellite-based estimates of groundwater depletion in India. Nature, 460(7258): 999-102, doi: 10.1038/nature08238.
[25] Shen Y Z, Chen Q J, Hsu H, et al. 2013. A modified short arc approach for recovering gravity field model. Presented at the GRACE Science Team Meeting, University of Texas, Austin, Oct.22-25.
[26] Swenson S, Wahr J. 2006. Post-processing removal of correlated errors in GRACE data. Geophys. Res. Lett., 33(8): L08402, doi: 10.1029/2005GL025285.
[27] Swenson S, Chambers D, Wahr J. 2008. Estimating geocenter variations from a combination of GRACE and ocean model output. J. Geophys. Res., 113(B8): B08410, doi: 10.1029/2007JB005338.
[28] Tapley B, Ries J, Bettadpur S, et al. 2005. GGM02—An improved earth gravity field model from GRACE. J. Geod., 79(8): 467-478, doi: 10.1007/s00190-005-0480-z.
[29] Tiwari V M, Wahr J, Swenson S. 2009. Dwindling groundwater resources in northern India, from satellite gravity observations. Geophys. Res. Lett., 36(18): L18401, doi: 10.1029/2009GL039401.
[30] 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. J. Geophys. Res., 103(B12): 30205-30229.
[31] Wang Z T. 2005. Theory and methodology of earth gravity field recovery by satellite-to-satellite tracking data[Ph. D. thesis](in Chinese). Wuhan: Wuhan University.
[32] Watkins M. 2012. JPL Level-2 processing stands document for Level-2 product release 005, Jet Propulsion Laboratory.
[33] Yi S, Sun W K. 2014. Evaluation of glacier changes in high- mountain Asia based on 10 year GRACE RL05 models. J. Geophys. Res. Solid Earth, 119(3): 2504-2517, doi: 10.1002/2013JB010860.
[34] You W. 2011. Theory and methodology of earth's gravitational field model recovery by LEO data[Ph. D. thesis](in Chinese). Chengdu: Southwest Jiaotong University.
[35] Zhang X F. 2007. The earth's field model recovery on the basis of satellite-to-satellite tracking missions[Ph. D. thesis](in Chinese). Shanghai: Tongji University.
[36] Zhao Q, Guo J, Hu Z, et al. 2011. GRACE gravity field modeling with an investigation on correlation between nuisance parameters and gravity field coefficients. Adv. Space Res., 47(10): 1833-1850. doi: 10.1016/j.asr.2010.11.041.
[37] Zheng W. 2007. Theory and methodology of earth's gravitational field recovery based on satellite gravity measurement[Ph. D. thesis](in Chinese). Wuhan: Huazhong University of Science and Technology.
[38] Zhou X H. 2005. Study on satellite gravity and its application[Ph. D. thesis](in Chinese). Wuhan: Institute of Geodesy and Geophysics, Chinese Academy of Sciences.
[39] Zou X C. 2007. Theory of satellite orbit and earth gravity field determination[Ph. D. thesis](in Chinese). Wuhan: Wuhan University.
[40] 罗佳. 2003. 利用卫星跟踪卫星确定地球重力场的理论和方法[博士论文]. 武汉: 武汉大学.
[41] 冉将军. 2014. 低低跟踪模式重力卫星反演理论方法及应用[博士论文]. 武汉: 中国科学院测量与地球物理研究所.
[42] 冉将军, 许厚泽, 沈云中等. 2012. 新一代GRACE重力卫星反演地球重力场的预期精度. 地球物理学报, 55(9): 2898-2908, doi: 10.6038/j.issn.0001-5733.2012.09.009.
[43] 冉将军, 许厚泽, 钟敏等. 2014. 利用GRACE重力卫星观测数据反演全球时变地球重力场模型. 地球物理学报, 57(4): 1032-1040, doi: 10.6038/cjg20140402.
[44] 王正涛. 2005. 卫星跟踪卫星测量确定地球重力场的理论与方法[博士论文]. 武汉: 武汉大学.
[45] 游为. 2011. 应用低轨卫星数据反演地球重力场模型的理论和方法[博士论文]. 成都: 西南交通大学.
[46] 张兴福. 2007. 应用低轨卫星跟踪数据反演地球重力场模型[博士论文]. 上海: 同济大学.
[47] 郑伟. 2007. 基于卫星重力测量恢复地球重力场的理论和方法[博士论文]. 武汉: 华中科技大学.
[48] 周旭华. 2005. 卫星重力及其应用研究[博士论文]. 武汉: 中国科学院测量与地球物理研究所.
[49] 邹贤才. 2007. 卫星轨道理论与地球重力场模型的确定[博士论文]. 武汉: 武汉大学.