GRACE (Gravity Recovery And Climate Experiment)卫星自2002年发射至今，在恢复地球重力场的精度和空间分辨率方面均有巨大进展，一直是卫星重力领域的研究热点.GRACE卫星对地球重力场的巨大贡献得益于星上同时装载了星载GPS接收机、K波段测距系统(KBR)、加速度计以及星敏感器等主要载荷.其中，K波段测距系统是恢复地球重力场模型的核心载荷，用于测量两颗卫星之间的距离及距离变率；星载GPS接收机主要用于确定GRACE卫星精密轨道；加速度计用于高精度测量卫星的非保守力；星敏感器用于确定卫星的精确姿态.以上载荷提供高低卫卫跟踪和低低卫卫跟踪观测数据，成功应用于反演静态地球重力场和探测全球质量变化信息.基于GRACE观测资料的时变重力场模型在全球、区域水储量变化研究，极区冰盖质量变化研究等方面取得了丰硕的成果(Chambers et al., 2004; Chen et al., 2006, 2009)，对地球物理学、大地测量学、海洋学和冰川学等领域的发展具有重要的意义.随着我国低低卫卫跟踪重力卫星计划的提出，自主获取高精度的时变地球重力场模型的相关研究势在必行.

如何从GRACE卫星观测数据中提取出所需的地球质量变化信息，是卫星重力领域的研究重点.基于重力卫星观测数据恢复地球重力场模型主要有以下方法：能量法、加速度法、短弧法、天体力学法和动力学积分法.到目前为止，国际上主要有美国德克萨斯大学空间研究中心CSR、德国地学研究中心GFZ、美国宇航局喷气推进实验室JPL、法国国家空间局CNES、德国波恩大学理论大地测量研究所ITG、荷兰代尔夫特地球观测与空间系统研究所DEOS以及瑞士伯尔尼大学天文研究所AIUB等机构，分别采用不同方法成功反演时变地球重力场模型(Mayer-Gürr, 2006; Bruinsma et al., 2010; Liu et al., 2010; Meyer et al., 2012; Bettadpur et al., 2012; Watkins and Yuan, 2012; Dahle et al., 2013).国内同行也为GRACE模式下利用实测数据反演地球重力场模型做了大量的工作.目前采用短弧法和动力学法得到与国际上相同精度水平的时变重力场模型，如冉将军等(2014)、Chen等(2015)和Shen等(2015)采用短弧法，利用KBR星间观测数据结合精密轨道成功反演了时变重力场模型；王正涛(2005)、周旭华(2005)、肖云(2006)、张兴福(2007)、王庆宾(2009)等人联合GRACE精密轨道和KBR星间观测数据，利用动力学积分法反演地球重力场的工作，取得了显著的研究进展.Zhao等(2010)和Wang等(2015)采用动力学积分法，利用KBR星间测速数据与低轨卫星精密轨道数据成功恢复了时变地球重力场模型.

在恢复地球重力场模型的方法中，能量法、加速度法和短弧法均需要已知重力卫星的精密轨道.基于GRACE卫星观测数据的动力学积分法可分为同解法和两步法，两步法是利用GRACE卫星精密轨道作为伪观测值，同时结合KBR星间观测数据反演地球重力场模型，国内相关机构已采用两步法成功恢复时变地球重力场模型(Zhao et al., 2010; Wang et al., 2015)；同解法(Zou et al., 2010)即直接采用原始观测数据同时确定GRACE卫星精密轨道和地球重力场模型.同解法使观测数据和观测模型达到更大程度的同化，理论上更加严密，实现更优的参数估计.目前只有少数国际知名机构(如GFZ、CSR、CNES等)利用该方法成功反演地球重力场模型，因此基于星载GPS和KBR星间观测数据采用同解法恢复时变地球重力场模型的相关研究对我国发展自主卫星重力测量系统以及卫星重力模型的研制理论和实际应用都有积极的意义.

本文利用同解法反演月尺度时变重力场模型，即采用GRACE星载GPS观测数据和KBR星间测速观测数据同时求解GRACE卫星精密轨道和60阶月尺度地球重力场模型，成功提取地球重力场的时变信号，并详细分析了观测数据拟合残差、轨道确定精度、地球重力场时变信号以及特征区域水储量变化等计算结果.

GRACE卫星装载了星载GPS接收机和K波段测距系统，同时获取星载GPS观测数据和KBR星间测速观测数据.星载GPS观测数据用于确定低轨卫星精密轨道并恢复地球重力场模型，仅采用星载GPS观测数据在恢复中高阶的重力场精度具有一定的局限性，通过结合KBR星间测速观测数据可以提高中高阶部分的精度.本文采用动力学积分法结合两类原始观测数据实现同时确定GRACE卫星的精密轨道和地球重力场模型.

基于动力学积分法反演地球重力场模型和精密定轨，根据卫星的运动方程将卫星的瞬间状态矢量表示为初始状态矢量和重力位系数及其他力学模型参数的函数，具体表示如下：

(1)

式中：时刻的卫星j (j=1, 2) 的位置矢量和速度矢量，X_j0表示卫星j的初始状态矢量，包括表示重力场位系数，包括{C_nm, S_nm}，U_j表示其他待估力学模型参数，包括加速度计尺度因子和偏差等参数.

对(1) 式用泰勒级数展开，则其线性方程表示如下：

(2)

式中：δX_j0, δP, δU_j分别表示为各参数的改正量.

星载GPS观测数据包括双频伪距和相位观测值，由于伪距观测相对相位观测精度低，为了减小计算量，本文在恢复重力场时不采用伪距观测值.为了消除电离层影响形成无电离层组合相位观测值，将其作为进行精密轨道确定和反演地球重力场模型的基本观测量.建立基于双频无电离层相位观测值的观测方程，见(3) 式，具体公式推导过程参考郭南男(2015).

(3)

其中ρ_o^L₃为低轨卫星相对于GPS卫星的无电离层组合相位观测值；ρ_c为低轨卫星相对于GPS卫星的几何距离，，其中(x_r, y_r, z_r)和(x^s, y^s, z^s)分别为低轨卫星和GPS卫星的三维位置坐标；c为真空中光速；δt_r和δt^s分别为接收机钟差改正和GPS卫星钟差改正；δρ_rel为相对论延迟改正；δρ_wind-up为相位缠绕改正；δρ_pco^r为低轨卫星相位中心改正；δρ_pco^s为GPS卫星相位中心改正；b₃为相位观测值的模糊度；ε为观测噪声.结合(2)、(3) 式得到基于星载GPS相位观测数据的观测方程：

(4)

由(4) 式，采用最小二乘法即可得基于星载GPS相位观测值的法方程.星载GPS观测数据中涉及到的各类观测误差可通过模型进行改正，最后仅解算模糊度和接收机钟差参数，其中接收机钟差为历元参数.最终法方程中的解算参数包括：初始轨道状态、加速度尺度因子和偏差、重力场模型球谐系数、模糊度参数和接收机钟差.除地球重力场模型参数为全局参数外，其余参数均为局部参数，若用Q表示局部参数，P表示全局参数，可以将法方程简化表示为：

(5)

KBR星间测距系统的主要观测量是星间距离ρ及其一次变率.星间距离对重力场的变化及非保守力测定误差均很敏感，对非保守力的测定精度要求很高，因此一般不用星间距离而采用其一次变率作为恢复重力场的基本观测数据.

采用动力学积分法利用KBR星间测速观测值的具体观测方程如下：

(6)

星间测速观测值是卫星瞬时状态矢量的函数，表示为：

(7)

星间测速观测值对于两颗卫星瞬间状态矢量的偏导数如下：

式中

其中分别表示两颗卫星在空间坐标系中的位置和速度矢量，δP表示重力位系数的改正值，δU表示加速度计尺度和偏差参数的改正值，δZ表示KBRR经验参数的改正值.同样用Q表示局部参数，P表示全局参数，可以将法方程简化表示为：

(8)

联合星载GPS相位观测和KBR星间测速观测法方程，本文根据两类观测数据残差确定两类观测值的权，相位观测设置σ_PHA=1 cm，星间测速观测设置为σ_KBRR=0.2 μm·s^-1，叠加两类观测的法方程具体表示为：

(9)

对(9) 式约化局部参数向量，得到仅含地球重力场位系数的法方程：

(10)

最后的法方程由叠加每一弧段对应的约化法方程形成，求解法方程即可得到重力场位系数，将全局参数回带到每个弧段的法方程中可分别得到两颗低轨卫星的精密轨道.

在数据处理过程各类求解参数的具体设置以及力学模型的配置是决定解算精度的关键，不同的研究机构采用的模型配置和参数设置各不相同，本文具体的模型配置和参数设置详见表 1.

表 1(Table 1)

表 1 力学模型与解算参数 Table 1 Mechanical models and calculation parameters

表 1 力学模型与解算参数 Table 1 Mechanical models and calculation parameters

加速度尺度因子和偏差参数的设置对GRACE卫星轨道和地球重力场模型的精度有着重要影响.本文采用Bettadpur (2009)计算的加速度尺度因子参考值作为初始值，但由于加速度尺度因子与动力学模型中的CPR (once-per-rev)参数有较高的相关性，因此文中不设置该参数.仅每一小时解算一组加速度偏差参数.

KBR经验参数(也称之为运动学参数)，包括常数项、一阶项以及周期项参数，其中周期项参数可避免轨道周期相关的误差，参数设置具体表示为(11) 式，其中a, b每45 min求解一组，c, d每1.5 h求解一组.

(11)

式中a, b, c, d为求解参数，T为卫星轨道周期，在这里为1.5 h.

为了证明上述算法的可靠性，下面将从以下五个方面对解算结果进行分析，包括相位观测拟合残差、KBRR拟合残差、低轨卫星定轨精度、重力场模型系数阶方差以及反演的陆地水储量变化的空间分布.

(1) 相位数据拟合残差

基于星载GPS观测数据确定卫星精密轨道和地球重力场模型位系数的整体解算精度，可以通过相位观测值的拟合残差反映.以2007年为例，选取两颗卫星一年的相位观测值以天为单位进行统计.图 1为2007年4月1日一天的GRACE卫星相位残差，相位残差的RMS为5.8 mm，从图中可以看出相位残差分布不存在系统误差，呈现随机分布.图 2为以天为单位统计2007年的相位残差RMS计算结果，从图中可以看出每天的拟合残差的RMS均优于1 cm，计算一年的相位残差RMS的平均值为6.18 mm，与Bruinsma等(2010)计算的相位拟合残差RMS的平均值6.1 mm精度相当.

图 1

Fig. 1

图 1 2007年4月1日相位残差 Fig. 1 Post-fit phase residuals on 1 April 2007

图 2

Fig. 2

图 2 2007年相位残差的RMS Fig. 2 Root Mean Squares of phase residuals in 2007

(2) KBRR观测数据拟合残差

KBRR观测值作为恢复地球重力场模型的关键数据，其拟合残差可以反映解算重力场模型精度.图 3以2007年4月1日为例，给出一天KBRR拟合残差计算结果，拟合残差的RMS为0.18 μm·s^-1.

图 3

Fig. 3

图 3 2007年4月1日KBRR残差 Fig. 3 Residuals of KBRR on 1 April 2007

从图 3中可以看出KBRR拟合残差分布不存在系统误差，呈现随机分布.图 4为统计2007年KBRR拟合残差结果，计算一年KBRR拟合残差RMS的平均值为0.201 μm·s^-1，与Bruinsma等(2010)计算的KBRR拟合残差的RMS平均值0.199 μm·s^-1精度相当，且其变化趋势一致.

图 4

Fig. 4

图 4 2007年KBRR残差RMS Fig. 4 Root Mean Squares of KBRR residuals in 2007

轨道误差影响地球重力场模型恢复精度.通过对轨道质量评估，不仅可以验证算法的正确性，同时可以避免较大的轨道误差影响地球重力场模型反演精度.以JPL提供的高精度的GRACE约化动力学精密轨道为标准轨道，将2007年4月1日定轨结果与之比较.图 5表示GRACE-A和GRACE-B在地固坐标系下与JPL发布的约化动力学精密轨道在三个方向的位置之差.GRACE-A在X、Y、Z三个方向位置之差的RMS分别为1.81 cm、2.02 cm、1.40 cm，三维位置之差的RMS为3.05 cm；GRACE-B在X、Y、Z三个方向位置之差的RMS分别为2.02 cm、1.54 cm、1.82 cm，三维位置之差的RMS为3.13 cm，与目前GRACE卫星每个方向的位置精度为2 cm的定轨精度相当(Jaggi et al., 2007；Kang et al., 2009).

图 5

Fig. 5

图 5 与JPL比较，定轨结果位置残差 Fig. 5 Differences of position between our solution and JPL

通过对比国际上知名机构CSR、JPL、GFZ发布的RL05模型，分析本文解算的60阶月尺度时变重力场模型精度.为了方便表述，本文解算的结果在文中均命名为SHA，而其他三个机构的RL05模型均用机构代号表示，分别表示为CSR、JPL、GFZ.分别计算包括本文在内的4个机构月尺度时变重力场模型的时变信号的阶方差.由于GRACE卫星本身的限制，不能准确地测定C₂₀项，不同机构解算的C₂₀项存在一定的差异.为方便比较，图 6给出2007年4月和10月从3~60阶的时变信号的阶方差结果.已有研究结果表明地球重力场模型的时变信号主要集中在前30阶.从图 6中可以看出，在前30阶部分，包括本文在内的四组模型的计算结果基本一致.随着阶次的增加，噪声会逐步增大.由于不同机构采用的具体数据处理策略不同导致噪声引入程度不同，30阶以后四组模型的时变信号阶方差均出现微小差异.本文结果在高阶部分阶方差略小于其他三组模型的结果，但从总体上本文模型时变信号阶方差与其他三组模型具有很好的一致性.

图 6

Fig. 6

图 6 SHA与CSR、GFZ、JPL的阶方差比较 Fig. 6 Comparison of degree variances of CSR, GFZ, JPL and SHA

为了进一步分析本文解算的时变重力场模型对全球质量变化信号的恢复程度，将地球时变重力信号用等效水高表示，并对比CSR、JPL、GFZ发布的RL05模型.具体处理过程如下：首先用卫星激光测距解算的C₂₀项替换时变重力场模型结果，采用P4M6去相关滤波去除同一次的奇数阶和偶数阶位系数的相关性，减少南北条带误差影响(Swenson and Wahr, 2006)，并采用500 km滤波半径的Gauss滤波(Jekeli, 1981)对时变重力场模型进行处理，计算月尺度全球水储量信号.图 7和图 8分别表示2007年4月和10月四组模型的全球时变重力场信号空间分布情况.对比CSR、GFZ、JPL模型计算结果，本文模型与其他三个机构反演的陆地水时变信号在全球范围内的空间分布十分接近，具有较高的空间一致性.

图 7

Fig. 7

图 7 2007年4月GRACE时变重力场模型计算全球质量分布异常(等效水柱高，单位：cm) (a) SHA；(b) CSR；(c) GFZ；(d) JPL. Fig. 7 Mass anomalies from GRACE in April 2007 (equivalent water column height in cm)

图 8

Fig. 8

图 8 2007年10月GRACE时变重力场模型计算全球质量分布异常(等效水柱高，单位：cm) (a) SHA；(b) CSR；(c) GFZ；(d) JPL. Fig. 8 Mass anomalies from GRACE in October 2007 (equivalent water column height in cm)

为了量化比较本文恢复的月尺度时变地球重力场模型与CSR、GFZ和JPL三个机构发布的RL05的差异，本文针对亚马逊流域和长江流域，分别对包括本文在内的四个机构的时变重力场在上述两个地区质量变化信号进行对比.同样对四个机构的月尺度时变重力场模型进行相同的高斯滤波和去相关处理.亚马逊流域选取[21°S, 5°N]，[45°W, 80°W]区域范围，长江流域选取[24°N, 35°N]，[90°E, 122°E]，分别针对该区域范围的纬度进行加权平均得到该区域的质量变化对应的等效水高.

如图 9所示，在亚马逊流域，四个机构结果都表现出一致的周期性变化，其主要反映了陆地水储量的变化.计算SHA、CSR、GFZ和JPL四个机构2005—2007年亚马逊流域平均水储量变化的周年振幅分别为26.2±3.3 cm、27.2±3.9 cm、27.3±4.0 cm和27.0±3.9 cm等效水柱高，SHA与其他三组产品的相关系数均大于0.98，呈现出较强的相关性.

图 9

Fig. 9

图 9 基于SHA、CSR、GFZ和JPL反演的时变重力场计算亚马逊区域陆地水变化(等效水柱高，单位：cm) Fig. 9 Water storage changes of the Amazon basin calculated by GRACE monthly solutions from SHA, CSR, GFZ and JPL (equivalent water column height in cm)

如图 10所示，在长江流域，2005—2007年SHA得到的等效水柱高与其他三个机构的周年信号基本一致.时变信号对应的等效水高的周期和振幅与CSR、GFZ和JPL这三个机构符合得较好，计算的相关系数分别为0.96、0.94和0.95.证明本文模型与其他三个机构的陆地水时变信号在区域范围的等效水高十分接近，具有较高的一致性.

图 10

Fig. 10

图 10 基于SHA、CSR、GFZ和JPL反演的时变重力场计算长江区域陆地水变化(等效水柱高，单位：cm) Fig. 10 Water storage changes of the Yangtze basin calculated by GRACE monthly solutions from SHA, CSR, GFZ and JPL (equivalent water column height in cm)

反演高精度地球重力场模型是卫星重力测量中的核心内容.本文基于卫星动力学方法，采用2005—2007年GRACE卫星的Level 1B资料，在建立卫星重力场和卫星轨道参数联合解算基础上，给出一组完整的60阶月平均地球重力场模型和卫星精密轨道.从观测数据拟合残差、轨道确定精度、重力场模型精度三个方面详细分析解算结果，得到如下结论：

(1) 统计2007年卫星重力场和卫星轨道参数联合解算后的星载GPS相位数据和KBR星间测速数据残差，分别为5~8 mm、0.18~0.30 μm·s^-1，与国际上公布的一些高精度结果接近.

(2) 与JPL发布的约化动力学精密轨道相比，本文计算的GRACE卫星轨道的X、Y、Z三个方向精度约为2 cm，三维位置精度优于5 cm，表明本文确定的GRACE卫星轨道与JPL公布的精密轨道精度相当.

(3) 与CSR、GFZ和JPL发布的RL05模型相比，本文解算的60阶月平均地球重力场模型与RL05模型精度相当，且在反演陆地水储量变化的区域等效水高差异很小，并具有很好的空间分布一致性.如亚马逊和长江区域，本文计算2005—2007年水储量变化与CSR、GFZ、JPL的RL05模型结果差异很小.对于亚马逊流域，其相关系数均大于0.98；对于长江流域，其相关系数分别为0.96、0.94和0.95，且年变化周期明显.

综上所述，本文联合GRACE卫星的星载GPS和KBRR星间测速数据，反演的时变重力场模型达到国外最新的时变重力场模型精度，且具有同时获得高精度卫星精密轨道的能力，为我国未来卫星重力反演和轨道确定提供另一重要技术途径.

感谢CSR、GFZ、JPL以及IGS提供的相关数据.