0 引言

自GRACE(Gravity Recovery and Climate Experiment)卫星发射以来，在地球重力场的中长波信息以及时变特征研究上取得了丰硕的成果(Tapley et al., 2004a, b).基于GRACE数据解算重力场模型是一项非常复杂的工作，从事这方面研究的国外机构主要有：德克萨斯大学空间研究中心(CSR，Center for Space Research)(Bettadpur, 2012)、美国宇航局喷气推进实验室(JPL，Jet Propulsion Laboratory)(Watkins and Yuan, 2012)、德国地学研究中心(GFZ，GeoForshcungsZentrum Potsdam)(Dahle et al., 2012)、波恩大学理论大地测量研究所(ITG，Institute of Theoretical Geodesy)(Mayer-Gürr et al., 2010)、代尔夫特地球观测与空间系统研究所(DEOS, Delft institute of Earth Observation and Space System)(Liu et al., 2010)、伯尔尼大学天文研究所(AIUB，Astronomical Insitute)(Meyer et al., 2012)和法国空间研究中心空间大地测量研究所(GRGS，Groupe de Recherche de Geodesie Spatiale)(Bruinsma et al., 2010)等；国内的主要研究单位有同济大学(Chen et al., 2015; Chen et al., 2016)、中国科学院测量与地球物理研究所(冉将军等，2014；王长青等，2015)、武汉大学(罗志才等，2016)、华中科技大学(杨帆等，2017)等.

通过分析不同机构的处理过程，可以发现GRACE数据处理过程大体上可以分为一步法、两步法和三步法.一步法是指将GPS(Global Positioning System)相关观测变量(如整周模糊度、钟差等)、卫星初始状态、卫星加速度计校准参数等以及重力场位系数等变量一起求解；两步法是指首先利用GPS数据完成定轨获得几何轨道或简化动力学轨道，然后将轨道积分中卫星初始状态，加速度计参数等与重力场位系数一起求解；三步法是指首先利用GPS数据完成定轨获得几何轨道或者简化动力学轨道，然后利用已有的重力场模型标定卫星初始状态和加速度计变量，待标定完之后解算重力场位系数.此外数据处理过程中还需进行迭代处理、引入经验参数等.事实上，不论是一步法、两步法还是三步法获得时变重力场的结果基本上一致，但一步法计算最复杂需要解算更多的变量，为了平衡计算效率和计算复杂程度本文选择两步法处理GRACE数据.

卫星在运行的过程中受到保守力和非保守力的共同作用，其中保守力包括日月引力、固体潮、海潮、极潮等，这些可以通过模型计算给出；而非保守力(包括大气阻力、光压等)采用卫星上搭载的加速计得到.利用摄动力模型和加速计计算得到卫星所受到的力与卫星真实受到的力不一致即存在摄动力模型误差，因此对重力场位系数的恢复有着直接的影响.如果不对这些误差进行处理，最直观的感觉是南北条带状影响明显增强(Han et al., 2004；Ditmar et al., 2012；Seo et al., 2008).对于摄动力模型的误差以及仪器噪声的影响，Kim(2000)从希尔方程(Kaplan, 1976；Colombo, 1984)出发，导出星间距存在1CPR(Cycle Per Revolution)误差，并提出采用经验参数来吸收误差.随后国内外一些机构也都使用经验参数来吸收误差，例如使用五参数或七参数误差模型(Visser, 2005；Liu et al., 2010; Zhao et al., 2011; Zhou et al., 2016)，并且每1.5小时标定一次参数.此外，本文根据加速度计以及误差模型分段标定的特点，为了使得被标定的量仍然保持数学上的连续性，提出采用三次样条拟合来吸收误差，并通过模拟计算来说明三次样条拟合的可行性.

尽管动力学法在处理GRACE数据时具有一定的优势，但是该方法依旧存在某种缺陷，主要表现误差随着弧长的增加而快速累积.为了使得弧长变的更长，Xu(2008)提出利用卫星的观测位置矢量作为参考轨迹，将卫星运动方程转换为Volterra型积分方程，但该方法在实际应用中存在一定计算难度.于锦海等(2017)通过引入非线性改正来补偿扰动量的快速增大，并证明了若参考轨道与真实轨道之间的距离差超过4.7 m，基于线性化处理的动力学方法将无法保证相应微分方程组具有3×10^-10m·s^-2的非惯性力测量精度；若引入非线性改正，则可以将轨道间的距离差从4.7 m延长至4100 m；这就极大放宽了动力学方法应用限制.

本文的主要工作是将于锦海等(2017)在动力学方法中引入的非线性改正项用于实际的GRACE数据处理，解算出相应的时变重力场模型.此外，还将直接对轨道扰动量对应的状态转移方程以及变分方程进行简化、求解，导出剩余星间距变率的低频误差特征，再根据非惯性力分段标定的特点提出采用三次样条拟合吸收低频误差的方法.最后利用GRACE Level-1b数据解算重力场，并将解得的时变重力场模型与GRACE官方公布结果进行对比，从数值结果上来验证相应结论的正确性.

本文主要分为：第1节是直接给出于锦海等(2017)引入非线性改正后相应扰动量对应的主要微分方程组，导出星间距变率观测方程，其目的是保证论文的理论完整性；第2节是研究状态转移方程与变分方程的求解，以及剩余星间距变率的误差特征，并提出三次样条拟合用来吸收低频误差；第3节采用模拟计算说明三次样条拟合吸收低频误差是否可行；第4节则是将反演的时变重力场模型与国际上主要研究机构GFZ、CSR、JPL的相应结果进行比较；最后是总结与结论.

1 顾及非线性改正的动力学方法与KBRR观测方程

为便于叙述，在无特殊说明的情况下，数学记号采用以下约定：使用r表示卫星的位置矢量(r表示相应的模)，用v或表示速度，对应的下标“0”表示初始状态，对应的下标“1”表示卫星的真实值，而下标“2”则表示相应的参考轨道值，上标“(A)”与“(B)”则表示对应于GRACE-A星与GRACE-B星的值，例如：表示GRACE-A星的参考轨道的速度矢量.

首先，对于单颗卫星来说，顾及非线性改正的轨道扰动方程组如下(于锦海等，2017)：

(1)

这里β=r₁－r₂是卫星位置扰动矢量，H_3×3(t)是引入参考轨道时用到的正常重力场的Hessian矩阵在参考轨道上的值，b₁(t)是非线性改正项，T_k是扰动位重新排列后的位系数，b_k(t)(k≥2)是各个阶次球谐展开函数沿轨道的值，K则是所需求解重力场球谐级数的阶次数.在方程组(1)中，待求的量是β与位系数T_k.

为了求解方程组(1)，将其分解成下面三个初始问题：

(2)

(3)

与

(4)

这里O₃和I₃分别是三阶零矩阵和单位矩阵，特别当k=1时方程(4)的解即为非线性改正项.在卫星精密定轨理论中方程组(2)与(3)通常称为状态转移方程组，而方程组(4)称为变分方程组.

从方程组的结构来看，状态转移方程组与变分方程组都是线性的，而且状态转移方程组恰好是变分方程组的齐次形式.因为上述初始问题都是线性的，因此它们的解存在且唯一.若以Φ_3×3(t)、Ψ_3×3(t)和B_k(t)分别表示它们的解，那么轨道扰动方程组(1)的解可以写成

(5)

在使用数值积分求解方程(2)、(3)、(4)的过程中，可同步得到

(6)

以上是针对单颗卫星建立的扰动理论.下面推导星间距变率观测方程.事实上，除轨道外GRACE卫星还提供两颗卫星之间的星间距变率数据，具体形式可写为，其中是GRACE-A星与GRACE-B星之间的位置差，α₁是其模，而则是两颗卫星间的速度差.关于星间距变率，可构建相应的观测方程如下：

(7)

这里

(8)

是对时间的导数，这里不再给出公式.在方程(7)中，具有下标“2”的量来自于参考轨道，因而是已知的，是测量值，待求的量是初始状态以及扰动位系数T_k，它们都包含在ε(t)和中，来自于轨道扰动方程组(1)的解，即：(5)与(6)式，自此构建了星间距变率观测方程.

2 状态转移方程组与变分方程组的求解、低频误差

对于GRACE卫星来说，其轨道是近圆轨道，而生成参考轨道时重力场的零阶项是主项(其中GM为地球引力常数)，因此略去J₂项影响(约为10^-3)后，Hessian矩阵可写成

(9)

由于圆轨道在卫星轨道面坐标系下的方程可写成

(10)

这里a是轨道半径，是平均角速度.从方程(10)的形式可知，O－xz面是卫星的轨道面，y轴的方向是轨道面法向.将(9)与(10)式代入状态转移方程组(2)或(3)，便知：状态转移方程组具有下列简化的形式：

(11)

记β=(β₁, β₂, β₃)^T是坐标系O－xyz下的分量形式，则有

(12)

与

(13)

显然，(13)式的通解是

(14)

此解反映了轨道扰动量在轨道面法向的变化特征.关于(12)式，利用矩阵对角化方法，可得

(15)

令

(16)

这样(12)式化简为

(17)

进一步化简，可得

(18)

由此可知，对η₁而言，有特征根0，±ωi；对η₃而言，有特征根±ωi以及重根0，i为虚数单位，因此可得通解为

(19)

从变换公式(16)可见，η₁是轨道扰动量的径向方向分量(外法向r₀)，而η₃正是轨道扰动量的运动方向分量(切向τ)，因此结合(14)式，便知(η₃, β₂, η₁)正是扰动量在卫星切向坐标系的分量(由于轨道是近圆轨道，故切向坐标系与径向坐标系基本一致).以n表示轨道面法向，则{τ, n, r₀}构成了卫星的运动坐标系，而且.事实上，可进一步证明方程(17)中第一式和第二式正是轨道扰动方程在外法向和运动方向下表示形式.

至此，在略去J₂量级(约为10^-3)的前提下完成了状态转移方程组的求解，即：其通解由(14)与(19)式构成，并且(η₃, β₂, η₁)就是轨道扰动量在运动坐标系{τ, n, r₀}下的分量.

虽然Colombo(1984)基于拉格朗日转换给出了Hill方程的通解，但在求解的过程较难发现其实际的物理意义，本文给出的通解计算方法是将卫星所在的惯性坐标系转换到卫星的运动坐标系，具有直观实际的物理意义便于未来的应用，当然这已经超出了本文的研究范围，不做详细讨论.现在讨论变分方程组(4)的解的结构.假设摄动力在运动坐标系{τ, n, r₀}下的分量分别是{F₃, F₂, F₁}，由于方程组(4)对应的齐次方程组恰好是状态转移方程组，所以在该运动坐标系下的变分方程组是(Kaplan, 1976；Colombo, 1984；Kim, 2000)

(20)

根据线性微分方程解的结构可知，无论摄动力模型误差{F₁, F₂, F₃}是何种类型的函数，必定会产生对应齐次方程组的通解这样的项，也就是说，关于时间的线性函数的长期影响项以及1个轨道周期(1CPR)影响总是存在的.特别，GRACE卫星计划中的非惯性力是以线性函数的形式来标定的，而分量η₃的特征根具有2个重根0，因此沿卫星运动方向会产生t³这样的长期影响.若将此时变分方程组的解以数学形式表达出来，则其解在运动坐标系{τ, n, r₀}中必定会包含下列分量：

(21)

由此可以导出剩余星间距变率的五参数误差公式

(22)

此外，若在生成参考轨道时顾及的摄动力模型误差具有1CPR的项，则该项对轨道扰动量的影响将是(C₁+D₁t)cosωt+(C₂+D₂t)sinωt的形式.由此可导出剩余星间距变率的七参数误差公式

(23)

其中表示低频项误差的影响形式，a_k(k=1, 2, …, 6)是待解的系数.

综合上述讨论可得，非惯性力标定以及摄动力模型的误差都会对剩余星间距变率产生长期影响与1CPR的影响，这些影响来自于状态转移方程组通解的结构.相对于重力场而言，这些影响的频率明显是属于低频的.所以可以通过引入低频参数来减弱这些影响项.

目前在处理GRACE的星间距变率数据时，常用的方法就是引入五参数或七参数误差模型来减弱剩余星间距变率低频误差的影响.事实上，除了非惯性力标定以及摄动力中存在1CPR误差外，非惯性力通常都是分段标定的模型不精确，因此剩余星间距变率中存在的误差具有分段表示的形式.为此根据剩余星间距变率的误差具有分段的特点，本文引入了样条函数来拟合低频误差项.考虑到时间t的三次多项式已经出现在非惯性力标定的误差中，将采用三次样条函数(李岳生和黄友谦，1978)来拟合低频误差项.具体公式为

(24)

其中c_j是待求的未知系数，步长h取值为卫星的轨道周期，三次样条基函数具体表达式为

(25)

3 模拟计算

为了说明三次样条函数吸收低频误差的可行性，本文采用模拟计算进行验证.本文选取EMG08为真实重力场，GIF48为参考重力场.事实上，卫星在运行过程中受到的模型误差主要为海潮模型误差(Han et al., 2004)，为此选择两个海潮模型的差值作为模型误差如：EOT11A和FES2004.此外还假定卫星在运行过程中位置含有均值为0标准差为2 cm正态分布的随机误差，初始位置和初始速度的偏差分别为3 cm和0.05 cm·s^-1，星间距变率含有均值为0标准差为1.0×10^－7m·s^-1正态分布的随机误差，恢复重力场的弧长为24 h，卫星运行的时间为9天.分别采用三次样条拟合和七参数拟合来吸收低频误差恢复重力场位系数，并将结果绘成阶方差进行比较如图 1所示.从图中可以看出三次样条拟合和七参数拟合都可以用来处理低频误差，恢复出的重力场位系数的阶方差精度相当，并且在某些阶次上三次样条拟合的精度要略好于七参数拟合，如2阶项，15~20阶.说明本文提出的三次样条拟合吸收低频误差的方法是可行的.

4 时变重力场模型的比较

本文在求解微分方程时采用了8阶Runge-Kutta和12阶Adams相结合的方法，而线性方程组求解则采用了Openmp和MKL(Math Kernal Library)函数库.计算参考轨道时采用的力模型有：已知的重力场模型、三体摄动、固体潮、海潮、极潮、非潮汐大气海洋改正、相对论效应和加速度计测量的非保守力.地固系与非惯性系转换采用IERS(International Earth Rotation Service)2010协议，具体描述见表 1.GRACE数据由JPL提供发布，数据产品为Level-1b，包括GNVB1B(GPS Navigation Data Product)，SCA1B(Star Camera Data Product)，ACC1B(Accelerometer Data Product)，KBR1B(K-Band Ranging Data Product)，GPS1B(GPS Data Product)等多种数据(Case et al., 2010).根据本文研究的需要，本文采用的数据为GNVB1B，SCA1B，ACC1B，KBR1B.

目前国际上已有机构采用1天或1周的数据成功反演时变重力场模型，但多数机构解算时变重力场主要采用的是月观测数据.为了便于比较，本文选择以月观测数据计算时变重力场.此外，本文对反演时变重力场的弧长选择为6 h，反演时变重力场截断阶数为60阶，并命名为UCAS_Grace01.对于每个弧段根据公式(5)和(7)来构建轨道扰动、星间距变率与卫星初始位置、初始速度、加速度计参数、经验参数以及时变重力场位系数之间观测方程.加速度计采取每个小时标定一次偏差，对尺度参数不进行标定，尺度参数和偏差初始值按照Bettadpur(2009)统计给出.剩余星间距变率低频误差采用样条函数(24)来吸收.由于单独采用星间距变率观测方程解算时变重力场模型时会导致法方程奇异，因此本文采用组合轨道观测方程和星间距变率的观测方程共同解算时变重力场模型，对于轨道观测方程和星间距变率方程组合采用的权为1×10¹⁰.

本文随机选择两个月的数据如2006年2月与2007年12月采用上述方法反演了月重力场模型，并与GFZ、CSR、JPL公布的RL05的月重力场模型进行了比较，结果如图 2所示.从阶方差上看在60阶范围内与国际上的机构基本上一致，说明本文提出的方法是有效可行的.

为了分析时变重力场模型UCAS_Grace01的解算精度，本文选择了2006年1月至2009年12月共48个月的GRACE Level-1b数据反演了月时变重力场模型，并将其与GFZ、CSR和JPL公布的RL05版的GRACE数据产品进行比较分析.在进行比较之前，本文首先利用卫星激光测距(SLR)结果替换了二阶项(Cheng and Tapley, 2004)，引入一阶地心改正项，对南北条带的影响采用Swenson和Wahr(2006)提出的去相关滤波减弱南北条带误差，并使用350 km的高斯平滑降低GRACE时变重力场高阶位系数的影响(Jekeli, 1981)，最后将球谐系数转化为全球分布的等效水高变化进行表示(Wahr et al., 2004).

图 3与图 4分别给出了2006年2月与2007年12月的月重力场模型与GFZ、CSR、JPL公布的RL05结果画出等效水高(Equivalent Water Height，EWH)的比较结果.从图 3与图 4可以得出，UCAS_Grace01反演得到的等效水高与GFZ、CSR、JPL公布的月重力场模型在空间上符合得较好，如在亚马逊流域和刚果河流域等陆地水的季节性变化较大的地方，四家机构反演的质量变化的等效水高空间分布基本一致.

此外还利用2006—2009年四家机构反演的时变重力场计算了质量变化的周年振幅(图 5).从图 5可看出，不同机构计算的质量变化周年振幅在全球分布基本上较为一致，如振幅较大的亚马逊流域、非洲的尼日尼河流域、亚洲的印度北部和中国青藏高原地区.

最后选取特征流域进行检核，特征流域的范围按照Oki和Sud(1998)给出，选择的特征流域主要为陆地水信号变化最大的亚马逊流域，陆地水信号变化较大的长江流域和陆地水信号变化最小的撒哈拉沙漠地区.分别计算四家模型等效水高的周年变化和年际变化(图 6).在亚马逊流域(图 6a)，UCAS_Grace01模型得到周年振幅等效水高为16.5±1.2 cm，基本上与GFZ(17.6±1.1 cm)、CSR(17.7±1.2 cm)、JPL(17.5±1.1 cm)相一致，时间序列相关系数分别为0.9969、0.9976、0.9980，符合较好.在长江流域(图 6b)，时间序列的相关系数分别为0.9823、0.9857、0.9864，均超过0.98，基本一致.在撒哈拉沙漠地区(图 6c)计算出的等效水高标准差分别为1.69 cm、1.74 cm、1.75 cm和1.78 cm，说明陆地水质量变化较小，UCAS_Grace01模型与GFZ、CSR和JPL模型在撒哈拉沙漠地区基本一致.

5 结论

本文基于解的叠加原理求解了轨道扰动方程组，即：公式(2)、(3)、(4)，构建了位系数与轨道和星间距变率的观测方程.独立推导了状态转移方程组的解，即：公式(14)与(19)；根据线性微分方程解的结构，分析了剩余星间距变率数据误差的分布特征，导出了五参数和七参数吸收低频误差的公式，特别是在推导的过程中采用惯性系向运动系的坐标转换，比Colombo(1984)采用的拉格朗日转换而言更能显而易见其实际物理意义，希望未来能有更多应用.此外根据非惯性力误差模式分段标定的特点，提出三次样条函数来吸收剩余星间距变率数据中的低频误差，通过模拟计算说明三次样条拟合函数与七参数恢复重力场模型的精度基本一致，并且在某些阶次要略优于七参数.最后本文处理实际的GRACE Level-1b数据(2006-01—2009-12)，解算了时变重力场模型UCAS_Grace01，并将所得模型与国外主要机构的对应结果进行了比对，结论是：(1)从全球范围上看，时变信号的空间分布与官方机构基本一致；(2)在局部区域上统计了时变信号最大的亚马逊流域，时变信号较大的长江流域和时变信号最小的撒哈拉沙漠地区，通过分析水储量变化时间序列的周年振幅、相关系数以及标准差，结果显示UCAS_Grace01模型与国际官方机构相比精度基本一致.

综上所述，本文基于顾及非线性改正的动力学法并采用了三次样条拟合函数吸收低频误差恢复时变重力场模型UCAS_Grace01与国际官方机构发布的时变重力场模型精度基本一致.

致谢  感谢JPL提供所需的GRACE数据，感谢两位匿名审稿专家对本文提出的宝贵意见.