地球物理学报  2017, Vol. 60 Issue (2): 514-526   PDF    
CHAMP型卫星定轨顾及非线性改正的轨道扰动方程
于锦海1,2 , 朱永超1,2 , 孟祥超1,2     
1. 中国科学院计算地球动力学重点实验室, 北京 100049;
2. 中国科学院大学地球科学学院, 北京 100049
摘要: 本文针对CHAMP型卫星建立了顾及非线性改正的轨道扰动方程定轨理论与方法.首先从卫星运动的二阶微分方程出发,引入了正常引力位以及相应的参考轨道,然后分别推导了线性化轨道扰动方程与顾及非线性改正的轨道扰动方程,同时说明了建立的线性化轨道扰动方程与目前处理CHAMP卫星数据的动力学定轨方法是等价的.其次分别对线性化轨道扰动方程与顾及非线性改正的轨道扰动方程的精度进行了估计,在卫星定位精度为3 cm与非惯性力测量精度为3×10-10m·s-2的前提下证明了下列结论:当参考轨道与实际轨道之间的距离ρ≤4.7 m时线性化轨道扰动方程的精度能达到非惯性力的测量精度以及当ρ≤4.14×103 m时顾及非线性改正的轨道扰动方程能达到非惯性力的测量精度.由此便可得出结论:相对于线性化轨道扰动方程,顾及非线性改正的轨道扰动方程具有更高的精度,且适合在更长的时间弧段上建立关于引力场位系数的法方程组,特别是针对CHAMP卫星计划进行的模拟计算也完全验证了该结论.最后利用叠加原理,给出了顾及非线性改正的轨道扰动方程的求解方法.此外,还针对GRACE卫星计划利用顾及非线性改正的轨道扰动方程进行了恢复引力场的模拟计算,结果表明:分段建立位系数的法方程组时子弧段分别取值2 h、1 d、6 d对恢复引力场的结果几乎不产生影响,这表明在处理GRACE数据时能够以6 d的弧长来建立法方程组.
关键词: CHAMP型卫星      非线性改正      轨道扰动方程      参考轨道      法方程组     
The orbital perturbation differential equations with the non-linear corrections for CHAMP-like satellite
YU Jin-Hai1,2, ZHU Yong-Chao1,2, MENG Xiang-Chao1,2     
1. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China;
2. College of Earth Science, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Different from the existed derivation methods, the linearized orbital perturbation differential equations for CHAMP-like satellite are derived directly from second order differential equations of satellite motion after introducing the normal gravitational field as well as the reference orbit, and then introducing the omitted terms into the linearized orbital perturbation differential equations, the orbital perturbation differential equations with the nonlinear corrections are proposed for CHAMP-like satellite. Since the omitted terms in derivations can be expressed clearly, the accuracies of the linearized orbital perturbation differential equations and the orbital perturbation differential equations with the nonlinear corrections are estimated respectively. By complicated computations, it is concluded from the above estimations that if the error for satellite position is less than 3cm and the error for non-gravitational accelerations is less than 3×10-10m·s-2, the linearized orbital perturbation differential equations can retain the measurement accuracy of non-gravitational accelerations when the distance ρ between the reference orbit and actual one is less than 4.7 m and the orbital perturbation differential equations with the nonlinear corrections can retain the same accuracy when ρ≤4.14×103 m. Hence, compared with the linearized orbital perturbation differential equations, the orbital perturbation differential equations with the nonlinear corrections have much higher accuracy if the same reference orbit is used. This means that the orbital perturbation differential equations with the nonlinear corrections are suitable to establish the normal system of linear equations about spherical harmonic coefficients of the gravitational field for longer time interval. In addition, with the help of superposition principle, solving methods for the orbital perturbation differential equations with the nonlinear corrections are given. Therefore, the complete theory and method about the orbital perturbation differential equations with the nonlinear corrections are established. Three arithmetic examples for CHAMP mission are imitated. By these examples, it is examined that the orbital perturbation differential equations with the nonlinear corrections have higher accuracies than the linearized orbital perturbation differential equations given the same reference orbit. In addition, it is examined by numerical imitations for GRACE mission that the span of time's segment making the orbital perturbation differential equations with the nonlinear corrections valid can be extended to 6 days at least if the measurement accuracy of the non-gravitational accelerations reaches to 3×10-10m·s-2. Obviously, the span of 6 days is much longer than the recommended span in dealing with GRACE data which is less than one day usually. GRACE-follow-on mission will be launched in the near future. The measurement accuracy of the non-gravitational accelerations in GRACE-follow-on mission may be raised to 10-12 m·s-2, and this means that the orbital perturbation differential equations with higher accuracy are required. Obviously compared with the linearized orbital perturbation differential equations, the orbital perturbation differential equations with the nonlinear corrections are more suitable to deal with GRACE-follow-on data..
Key words: CHAMP-like satellite      Nonlinear correction      Orbital perturbation differential equations      Reference orbit     
1 引言

利用卫星轨道数据解算地球重力场的理论与方法很早就被提出了(例如:Kaula, 1966; Taff,1985Reigber, 1989Wolff, 1969),早期主要是利用激光测距等手段来求解重力场的低阶次位系数.随着测量技术的发展,先后实施了诸如CHAMP、GRACE以及GOCE专项重力卫星计划(ESA,SP-1233(1),1999;Reigber,2012;Rummel,2012),卫星重力才取得了丰硕的成果.目前利用重力卫星观测数据解算的重力场模型已经成为了许多学科的基础数据.正因为如此,研究或拓展利用重力卫星观测数据来解算重力场的方法无论在理论上、还是在应用上都具有很重要的意义.为了表述的方便,这里先给出CHAMP型卫星的定义.CHAMP卫星搭载了GPS接收机、姿态照相机以及加速度计,因此CHAMP卫星的位置与作用于卫星的非惯性加速度是可观测的,本文将称具有这样观测量的卫星为CHAMP型卫星.按照CHAMP型卫星这样的定义,无论是GOCE卫星,还是单个GRACE卫星都属于CHAMP型卫星.这意味着CHAMP型卫星是构成卫星重力计划的基础,因此研究CHAMP型卫星的定轨理论在处理重力卫星观测数据中起着基础性的作用.

众所周知,卫星的运动方程为

(1)

其中r是卫星在地球惯性系中的位置矢量,g (t, r)是地球引力,而a (t, r, )则是除了地球引力外作用于卫星的其他力的合力(包含了第三体引力、固体潮及海潮的附加引力、大气阻力、光压等).若用u表示引力位,则有g=gradu以及在地固坐标系下u可写成

(2)

这里GM是万有引力质量常数,a是赤道半径,λθ分别是经度和余纬,Pnm(t)是规范化的n阶m次连带Legendre函数,anmbnm便是u的球谐展开系数.虽然(2)式中和式对n求和理论上需要达到无穷大,但实际处理数据过程中仅取到某个特定的阶数,常用N表示,例如:N通常可取为N=60, 120, 180等.

若卫星运动方程(1)有解析解,那么在CHAMP型卫星的观测数据与系数anmbnm之间便能够建立起相应的函数关系(该函数关系通常称为法方程组),然后求解法方程组便可得到系数anm、bnm,从而得到地球重力场模型(或引力场模型).然而方程(1)是非线性极强的微分方程组,至少目前仍然无法对其进行解析求解,因而便产生了卫星轨道扰动理论.所谓卫星轨道扰动理论就是基于某种线性化途径来建立观测数据与系数anm、bnm之间的函数关系,这在处理卫星重力数据时也常被称为动力学定轨方法.

通常有两种途径来实现对观测数据的线性化处理.第一是(Montenbruck and Gill, 2000)引入状态转移矩阵(State transition matrix)和参数响应矩阵(Sensitivity matrix),然后通过求解变分方程求出这些矩阵,最后利用这些矩阵来建立初始状态、模型参数(包含了系数anm、bnm)与观测数据之间的线性观测方程组;第二(Reigber, 1989; Tapley, 1989; Tapley et al., 2004a, 2004b)是将卫星运动的二阶微分方程组化简为一阶微分方程组,引入参考轨迹后将该一阶微分方程组在参考轨迹处进行Taylor展开得到相应的线性化轨道扰动微分方程组,然后由此解得初始状态、模型参数与观测数据之间的函数关系.虽然形式上这两种途径存在着差异,但是本质上它们是一致的,其原因是在推导过程中只有线性扰动项被保留下来了.

除了上述讨论的动力学定轨方法外,还有其他的处理CHAMP型卫星观测数据的方法.例如:Douglas等(1980)Jekeli (1999)Visser等(2003)讨论过能量法;Ditmar等(2004, 2006, 2012)利用CHAMP卫星的加速度恢复了重力场模型.需要指出的是:在这些方法中使用了诸如卫星速度和加速度这样的量,这些量不是CHAMP型卫星的直接观测量,而是从位置数据派生出来的,因此这些方法是存在瑕疵的.为了能够直接处理直接观测量,Mayer-Gürr等(2005)IlK等(2005, 2008)、Rowlands等(2002)Wagner (2006)还将边值问题用于处理CHAMP或GRACE数据,并称这样的方法为短弧法.理论上看,短弧法与基于初始问题的动力学定轨方法一样都是直接处理观测数据,但要注意的是二阶常微分方程组边值问题在特定的边界条件下会出现解不唯一的情况,这也许正是短弧法对求解弧段需要进行约束的原因,即:弧段较长时会产生无解或多解的情况.另外还有基于Kaula解析方法的研究(Kaula,1966Cui and Lelgemann, 2000)以及Boutler等(2010)对常用方法进行了分析比较;Chen等(2004)Cheng (2002)对模型结果进行了验证与分析.

尽管在现有的处理卫星重力数据的方法中动力学定轨方法具有一定的优势,但是该方法仍存在着某些缺陷.例如:Xu (2008)曾指出对于动力学定轨方法而言线性化产生的误差随着弧段长度的增加会迅速地增大,因而在应用该方法时弧段的长度仍然受到限制.为了使得弧段的长度取的更长,Xu提出了利用卫星观测的位置矢量作为参考轨迹.由于这样的参考轨迹不再满足微分方程,所以Xu将卫星运动的微分方程转换成Volterra型积分方程,然后再进行线性化处理.由于利用位置观测数据得到的参考轨迹与实际轨道的偏差很小,所以线性化产生的误差就会相应减小,从而能够在更长的弧段内来针对重力场的位系数建立法方程组.理论上讲,Xu提出的方法确实可以减少线性化产生的误差,然而相应的计算复杂性也会有所增加,从而在应用上有一定的难度.

本文的目标是针对CHAMP型卫星的特点,在动力学定轨方法的基础上引入非线性改正从而提高定轨结果的精度,以便达到在更长的时间弧段内能够处理重力卫星的观测数据.为了便于理解与导出非线性改正项,本文以某种新的方式对动力学定轨方法进行了推导,该推导方法具有下列优点:第一、推导过程简洁,便于理解;第二、能够有效地估算线性化过程中略去的误差;第三、能够引入非线性改正项,从而提高定轨的精度;第四、相对于目前已有的关于动力学定轨方法的计算而言,解算的复杂性并没有增加.此外,本文还进行了系列的模拟计算,结果表明:无论是计算精度、还是时间弧段的延长,引入非线性项后的计算结果都比动力学定轨方法的计算结果有较大的提高.

为了使得后面的推导过程便于理解,这里事先对卫星轨道微分方程(1)的相关记号作出说明.例如:若用r来表示卫星的位置矢量,则用来表示r的观测值,而用ε r=r-ε r=|ε r|来表示的观测误差.同样地,对方程(1)中的a(t, r, )来说,也有相应的记号εaεa.对于方程(1)中的非地球引力项a(t, r, )而言,由于非惯性力是可以直接测量的,而第三体引力场、固体及海洋潮汐的附加引力等都可以进行模型化模拟,所以为了推导简单起见后面都假设a(t, r, )是可以直接观测的.也就是说,在方程(1)中ra(t, r, )都是观测量,其观测值分别用 (t)和(t)来表示之,而εrεa则分别是相应的观测误差.在本文中将应用GRACE计划的观测精度指标,即对于εrεa,分别假设maxεr≤3×10-2m和maxεa≤3×10-10m·s-2.

2 CHAMP型卫星的轨道扰动方程 2.1 参考轨道与轨道方程的变形

首先引入正常引力位U,该正常引力位应该是地球真实引力位u的足够好的近似,例如:U可以取为EGM2008等引力场模型.然后引入下列微分方程的初始问题:

(3)

这里初始位置r1(0)与速度v 1(0)都是已知的,并且取的足够接近CHAMP型卫星相应的真实值.由于(3)中的U(t)以及初始值r1(0)v 1(0)都是已知的,所以初始问题(3)有唯一的解,记为r=r1(t).以后称r=r1(t)表示的空间曲线为对应于U和(t)的参考轨道,或简称为参考轨道.令T=u-U,则T是扰动位,在地固坐标系下T有球谐展开式:

(4)

其中AnmBnmT的无量纲的球谐系数,而Rnm(θ, λ)与Snm(θ, λ)的表达式为

(5)

使用记号r=r2(t)表示卫星轨道方程(1)的解,并记r=r2(t)在t=0时刻的初始位置与速度分别是r2(0)v 2(0),那么则有

(6)

由于u=U+T,所以方程(6)可改写为

(7)

引入记号ρ=r2-r1,并将方程(7)与方程(3)对应相减,则有

(8)

其中ρ0=r2(0)-r1(0), η0=v 2(0)-v 1(0).在方程(8)中将gradU(t, r2)在r1处作Taylor展开,则有

(9)

这里运算都是在地固坐标系中进行的,而

(10)

是由U的二阶导数构成的3阶矩阵(通常称HU的Hessian矩阵),M(t, r1, ρ)则是O(ρ2)的高阶小量,其中ρ=|ρ|.事实上,M(t, r1, ρ)可表示为

(11)

这里Mj(t, r1)(ρ)是ρj次齐次多项式,例如:M2(t, r1)(ρ)和M3(t, r1)(ρ)的第i个分量可分别表示为

(12)

(13)

将(9)式代入(8)式,便有

(14)

这里E(t)是从惯性系到地固坐标系的转换矩阵.

在推导方程(14)时没有舍去任何项,因此方程(14)是严格成立的.事实上,方程(14)仅是在引入参考轨道并将求解函数转换成ρ后方程(1)的变形.

2.2 线性化的轨道扰动方程

如果在方程(14)中略去εa(t)和O(ρ2)、O(Anm·ρ)、O(Bnm·ρ)量级的小量,则方程(14)可写为

(15)

由于r1=r1(t)是已知的,所以H3×3(t, r1(t))仅是时间t的函数,这意味着方程(15)关于求解对象ρ是线性的.以后将方程(15)称为线性化轨道扰动方程.

注意:本文推导线性化轨道扰动方程的出发点是基于卫星运动的二阶微分方程(1),这与Tapley (1989)Montenbruck和Gill (2000)的推导方法不一样,然而其结果在本质上应该是等价的.

2.3 线性化轨道扰动方程的误差估计

根据从方程(14)到方程(15)的简化过程,容易看出被舍去的量分别是εa(t)、M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1).因为εa(t)是a(t, r, )的测量误差,其值是固有的,所以在舍去M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1)时必须参考εa(t)的大小.本文中已经假定误差εa(t)满足maxεa≤3×10-10 m·s-2,因此下面将δ=1.5×10-10m·s-2取为被舍去量的基准,即只有舍去量小于δ时,才认为舍去是合理的.事实上,取δ为舍去基准的目的就是为了保证在轨道扰动方程中舍去的误差不会超过测量误差,这样才能尽可能地在测量误差的范围内寻求出重力场的解.下面开始分别对M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1)进行估计.

2.3.1 M(t, r1, ρ)的估计

因为M2(t, r1)(ρ)分别是UM(t, r1, ρ)的主项,所以仅需去估计当U=U0M2(t, r1)(ρ)的值.根据(12)式,可得

(16)

对于CHAMP型卫星来说,r的变化范围几乎总在6700 km与6800 km之间,于是近似地有≈ 9.0 m·s-2.令|M2(t, r1)(ρ)|U=U0≤δ,即27×≤1.5×10-10,由此可得ρ≤14 m.这意味着只有当ρ≤14 m时,在方程(14)中舍去M(t, r1, ρ)才是合理的.

2.3.2 gradT(t, r2)-gradT(t, r1)的估计

根据(4)式,gradT(t, r2)-gradT(t, r1)可写成

(17)

因为分别等价于,所以从(17)式出发可得

(18)

注意恰好是扰动位T的阶方差.就目前对引力场模型精度分析的结论来说,若将正常引力位U取为精度较高的引力场模型,则其对应的扰动位的阶方差大致能够满足=O (10-10).另外还有下列近似表达式

(19)

因此便有

(20)

,注意到r1的变化范围大致在6700 km与6800 km之间,由此便可得到ρ≤4.7 m.这意味着当ρ≤4.7 m,舍去gradT(t, r2)-gradT(t, r1)是合理的.

总之,如果a(t, r, )的测量误差εa(t)能够达到精度maxεa≤3×10-10 m·s-2,而选择的正常引力位足够逼近真实的地球引力位以使得相应的扰动位阶方差具有量级O(10-10),那么只要参考轨道与真实轨道的差ρ≤4.7 m,此时线性化轨道扰动方程就能控制在测量误差εa(t)的范围内.事实上,轨道差ρ的大小也是由正常引力位U接近于真实引力位的逼近程度和参考轨道的初始条件决定的.因此正常引力位U与参考轨道初始条件的选择对于线性化轨道扰动方程的精度起着至关重要的作用.如果U取的较为简单,则轨道差ρ会较大,从而线性化轨道扰动方程的精度就会降低;如果将U取成较高阶次的引力场模型,虽然线性化轨道扰动方程的精度会有所提高,但是也会在计算方面带来更多的复杂性.

2.4 顾及非线性改正的轨道扰动方程

在推导线性化轨道扰动方程时诸如M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1)这样O(ρ2)以上的小量被略去了.对于CHAMP型卫星而言,卫星位置是可测量的,其观测值为2,因此在轨道变形方程(14)中O(ρ2)量级的项M(t, r1, ρ)和gradT(t, r2)能够通过观测值2进行计算,即它们可以分别被M(t, r1, )和gradT(t, 2)替代.在略去测量误差εa(t)后,轨道变形方程(14)便可写为

(21)

这里(t)=M(t, r1, E(t)),

(22)

以及

(23)

以后称方程(21)为顾及非线性改正的轨道扰动方程.与线性化轨道扰动方程(15)比较可见,方程(21)中不是简单地舍去O(ρ2)以上的小量M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1),而是使用观测值 2M(t, r1, ρ)和gradT(t, r2)-gradT(t, r1)进行了计算.事实上,M(t, r1, )和gradT(t, 2)-gradT(t, r1)就是非线性改正项.显然添加非线性改正项后,被舍去的项分别是εa(t)、M(t, r1, )-M(t, r1, ρ)和gradT(t, 2)-gradT(t, r2).

2.4.1 M(t, r1, )-M(t, r1, ρ)的估计

为了分析方程(21)的精度,需对| M(t, r1, )-M(t, r1, ρ) |和| gradT(t, 2)-gradT(t, r2) |进行估计.类似于(16)式的推导,可得

(24)

这里εr2 2的测量误差.根据本文对卫星位置的测量误差假设(GRACE误差指标),可知εr2≤0.03 m,因此便有

(25)

如果仍以δ作为舍去误差标准,则有2.428×10-7×≤1.5×10-10,从而得ρ≤6.178×10-4r1.对于CHAMP型卫星而言r1总是满足r1≥6700 km,于是可得结论:只要ρ≤4.14×103 m,使用M(t, r1, )或(t)去替代M(t, r1, ρ)时方程(14)的误差就会恒小于δ.

2.4.2 gradT(t, 2)-gradT(t, r2)的估计

类似于(17)式,可得

(26)

从而有

(27)

如果扰动位T的阶方差仍然满足,注意到εr2≤3×10-2m以及r的变化范围在6700 km与6800 km之间,那么沿用第2.3.2节的推导便得

(28)

这意味着使用gradT(t, 2)去替代gradT(t, r2)时方程(14)的误差总会小于δ.

事实上,此时即使扰动位的阶方差满足,仍然可以使用gradT(t, 2)去替代gradT(t, r2).这说明了如果顾及了非线性改正项,则正常引力位U的选择要求可以降低许多,而此时的轨道扰动方程仍能满足测量精度εa(t)的要求.

总之,对于顾及非线性改正的轨道扰动方程而言,若仍以εa(t)作为精度指标,则此时对ρ的限制要比线性化轨道扰动方程的限制弱得多.例如:在扰动位的阶方差满足=O(10-10)的情况下,线性化轨道扰动方程对ρ的限制是ρ≤4.7 m;而顾及非线性改正的轨道扰动方程对ρ的限制是ρ≤4.14×103 m.

由于r1(t)和 2(t)都是时间t的已知函数,所以将它们代入方程(21)后,顾及非线性改正的轨道扰动方程(21)便可写成下列形式:

(29)

这里3×3(t)=ET(t)H3×3(t, r1(t))E(t)、b1(t)=ET(t)(t)、而则是经过重新排列后的表示,Tk就是扰动位的位系数.

事实上,线性化轨道扰动方程(15)也可以写成方程(29)的形式,此时没有b1(t)项.这就是说,无论线性化轨道扰动方程、还是顾及非线性改正的轨道扰动方程,其最终的数学形式是一致的,因而它们在解算难度上也是一样的.但需要指出的是:在非惯性力精度maxεa≤3×10-10m·s-2下,顾及非线性改正的轨道扰动方程对ρ的要求是ρ≤4.14×103 m,而线性化轨道扰动方程对ρ的要求是ρ≤4.7 m,因此解算顾及非线性改正的轨道扰动方程的弧段要比线性化轨道扰动方程的弧段长得多.

注意在计算M(t, r1, )时需要计算U的高阶导数,为了避免U的高阶导数计算时带来的复杂性,这里利用(9)式给出关于M(t, r1, )的简单表达式如下:

(30)

3 轨道扰动方程的解法

本节讨论方程(29)的求解方法.注意到方程(29)关于ρ是线性的,因此可使用叠加原理对其进行求解.事实上,方程(29)可分解为

(31)

(32)

(33)

这里O3I3分别是三阶零矩阵和单位阵.因为初始问题(31)、(32)和(33)都是线性的,使用他们的解是唯一存在的,以Φ3×3(1)(t)、Φ3×3(2)(t)和Bk(t)分别表示它们的解,那么便有

定理1  方程(29)的解为

(34)

因为Φ3×3(1)(t)、Φ3×3(2)(t)和Bk(t)的导数3×3(1)(t)、3×3(2)(t)和k(t)在求解初始问题(31)、(32)和(33)时是可以与Φ3×3(1)(t)、Φ3×3(2)(t)和Bk(t)同步解出的,所以在得到(34)式时,还可以同步得到

(35)

至此,关于方程(29)的求解转换为初始问题(31)、(32)和(33)的求解.事实上,初始问题(31)、(32)和(33)在某种意义下等价于卫星轨道精密定轨理论中解算状态转移矩阵和参数响应矩阵的变分方程(Montenbruck and Gill, 2000).

下面介绍方程(29)的求解结果(34)与(35)式的潜在应用.

3.1 CHAMP卫星计划

对于CHAMP卫星而言, (t)是可测量的量,因此利用(34)式,在观测点上便有

(36)

这里tk是给出测量点的时刻.从(36)式极易看出,该式就是关于未知数ρ0η0Tk(k=2, …, K)的线性函数.因此只要观测数据足够多,便能构建关于这些未知数的法方程组,从而解得扰动位的位系数Tk.

3.2 GRACE卫星计划

GRACE卫星计划由两颗卫星GRACE-A和GRACE-B组成,其中每颗星都是CHAMP型卫星,而两颗卫星之间的距离变化率是主要的观测量.为了表述方便,下面分别使用上标(A)和(B)来表示对应于卫星GRACE-A和GRACE-B的观测量.

根据(35)式,对于卫星GRACE-A和GRACE-B来说,分别有

(37)

(38)

因为

(39)

所以将(37)与(38)式代入(39)式,便得

(40)

注意在(40)式中 1(A)- 1(B)3×3(A1)(t)、3×3(A2)(t)、3×3(B1)(t)、3×3(B2)(t)以及k(A)(t)-k(B)(t)都是已经解出的已知量,而 2(A)- 2(B)可是观测量,因此(40)式就给出了关于扰动位的位系数Tk(k=2, …, K)的法方程组.

4 算例验证

本节将做系列模拟CHAMP与GRACE卫星计划的计算,以便验证顾及非线性改正的轨道扰动方程比线性化的轨道扰动方程在恢复地球引力场方面具有更高的精度,而且更适合用于在更长弧段上建立法方程组.为了避免计算的复杂性对验证效果的影响,假设模拟计算中卫星轨道都是仅由引力场生成的.显然在上述假设下,方程(1)中的作用力a(t, r, )=0.

将前60完整阶次的EGM2008引力场模型作为真实的引力场,然后利用该引力场生成半长轴近似为6820 km、倾角近似为89°、偏心率近似为0.002的卫星轨道.在生成卫星轨道以及后面微分方程的数值求解时本文使用的解法都是8阶Runge-Kutta算法与8阶Adams算法的组合,即前10次迭代计算使用Runge-Kutta算法,以后的迭代计算则使用Adams算法,计算的步长均取为5 s.用r(t)表示生成的轨道的位置矢量,则r(t)是已知的(依次5 s间隔给出r(t)的值).下面就从已知的r(t)(即 2(t))的值分别利用线性化的轨道扰动方程与顾及非线性改正的轨道扰动方程来恢复真实的引力场.

例1  选择前完整60阶次的EGM96引力场模型作为正常引力场,并以与真实轨道相同的初始位置与速度来生成参考轨道(得到了r1(t)),然后分别利用线性化的轨道扰动方程与顾及非线性改正的轨道扰动方程(即(36)式)来恢复真实的引力场.在具体计算中对时间弧段分别选取了4组值:3、6、9、12天.恢复的引力场模型的误差阶方差由图 1给出.图中LOPE与NLOPE分别表示使用线性化与顾及非线性改正的轨道扰动方程恢复位系数的阶方差.

图 1 2(t)不存在观测误差、正常引力场是60阶次、参考轨道与真实轨道的初始条件相同时恢复的引力场的误差阶方差.(a)、(b)、(c)、(d)是时间弧段分别为3、6、9、12天的结果 Fig. 1 The error degree-variance of the recovered gravity field when there are no any errors in 2(t), the normal gravity field is chosen as the former 60 degree and the reference orbit has the same initial conditions as true orbit, where Fig. (a), (b), (c) and (d) are corresponding to time segments of 3, 6, 9 and 12 days respectively

图 1可见,利用顾及非线性改正的轨道扰动方程恢复的引力场的精度明显要优于线性化的轨道扰动方程的结果,特别是随着时间弧段的增长,利用线性化轨道扰动方程(例如弧段长度超过9天)不再能够完全恢复真实引力场了,而利用顾及非线性改正的轨道扰动方程恢复仍能至少将引力场恢复2个量级以上.

例2  在例1中仅变换初始条件如下:参考轨道的初始位置与速度和真实轨道的相应初始值分别偏差3 cm与0.05 cm·s-1,然后进行相同的计算,其恢复的引力场误差阶方差由图 2给出.

图 2 2(t)不存在观测误差、正常引力场是60阶次、参考轨道与真实轨道的初始条件存在偏差时恢复的引力场的误差阶方差.(a)、(b)、(c)、(d)是时间弧段分别为3、6、9、12天的结果 Fig. 2 The error degree-variance of the recovered gravity field when there are no any errors in 2(t), the normal gravity field is chosen as the former 60 degree and the initial conditions of the reference orbit and true orbit have small differences, where Fig. (a), (b), (c) and (d) are corresponding to time segments of 3, 6, 9 and 12 days respectively

综合图 1图 2可得出下列结论:第一、顾及非线性改正的轨道扰动方程恢复引力场的效果明显优于线性化轨道扰动方程;第二、利用顾及非线性改正的轨道扰动方程来恢复引力场时可以有效地增加时间弧段;第三、参考轨道初始条件对线性化轨道扰动方程恢复引力场的结果影响很大,而参考轨道初始条件的偏差对顾及非线性改正的轨道扰动方程的影响不大.在例2中初始位置与速度的偏差是依据CHAMP和GRACE的实际观测误差给出的,而这样的偏差能使得在长时间的弧段(例如:6天)上无法使用线性化轨道扰动方程来有效地恢复真实引力场.这也许正是Reigber等(2002)Tapley等(2004)建议使用线性化轨道扰动方程来构建引力场位系数的法方程组时时间弧长不要超过1天的原因.反观顾及非线性改正的轨道扰动方程,在初始条件存在偏差时即使时间弧长达到12天,仍然能够有效地恢复出真实的引力场.

另外需要说明的是,尽管可以通过不断迭代的方式来改进线性化轨道扰动方程的计算精度,但由于观测误差的原因,参考轨道初始条件的偏差总是存在的,因此如何控制线性化轨道扰动方程的迭代次数在应用上是很难把握的.与之相对应,由于参考轨道初始条件的偏差对于非线性改正的轨道扰动方程的计算结果影响很小,因此这也证明了非线性改正的轨道扰动方程的抗干扰性更强.

例3  将正常引力场选择为完整前3阶次EGM96模型,参考轨道的初始条件按例2给出,然后进行例2的计算,其结果由图 3给出.

图 3 2(t)不存在观测误差、正常引力场是3阶次、参考轨道与真实轨道的初始条件存在偏差时恢复的引力场的误差阶方差.(a)、(b)、(c)、(d)是时间弧段分别为3、6、9、12天的结果 Fig. 3 The error degree-variance of the recovered gravity field when there are no any errors in 2(t), the normal gravity field is chosen as the former 3 degree and the initial conditions of the reference orbit and true orbit have small differences, where Fig. (a), (b), (c) and (d) are corresponding to time segments of 3, 6, 9 and 12 days respectively

图 3可见,除了例1与例2得到的结论外,选择较为简单的正常引力场或参考轨道对顾及非线性改正的轨道扰动方程用于恢复真实引力场的影响仍然较小.例如:即使正常引力场仅取到3阶次,但时间弧长仍然可以达到了6天.从轨道扰动方程的建立来看,Hessian矩阵H3×3(t, r1(t))是由正常引力位U的二阶导数构成的,若U的表达式简单,则H3×3(t, r1(t))的计算量将会极大地减少.因此从计算的角度看,选择顾及非线性改正的轨道扰动方程无疑是非常有益的.

上述3个例子是针对CHAMP计划给出的,现在转向GRACE计划的模拟计算.依然用前60完整阶次EGM2008引力场模型作为真实引力场,GRACE两颗卫星的轨道设计仍与上面一致,此外两颗卫星之间的距离近似地等于220 km.这里将时间弧段的长度固定为6天,首先计算出两颗卫星的位置r2(A)(t)、r2(B)(t)与速度 2(A)(t)、 2(B)(t);然后计算出卫星间的距离变化率.下面的工作就是从这6天中已知的数据r2(A)(t)、r2(B)(t)和v2(AB)来恢复真实引力场.采用的是常用的分段建立法方程组的方法,即将整个6天的时间弧段分成若干个子弧段,在每个子弧段上使用顾及非线性改正的轨道扰动方程(即(40)式)来建立法方程组,然后合并所有子弧段上的法方程组对扰动位的位系数进行求解.本文具体的子弧段分段采用了3种,分别是2 h、1天、6天(或整个弧段).

例4  将6天的数据r2(A)(t)、r2(B)(t)和v2(AB)作为观测值,取完整前60阶次EGM96引力场模型作为正常引力场,在子弧段给出参考轨道的原则是参考轨道的初始条件与真实轨道的值相同.然后利用(40)式来恢复真实的引力场.恢复的引力场的误差阶方差由图 4给出.

图 4 观测数据r2(A)(t)、r2(B)(t)、v2(AB)没有误差、子弧段上参考轨道与真实轨道的初始条件相同时恢复的引力场的误差阶方差,这里“蓝线”表示子弧段为2 h的结果、“红线”表示子弧段为1天的结果、“黑线”表示弧段没有分段的结果 Fig. 4 The error degree-variance of the recovered gravity field when there are no any errors in r2(A)(t)、r2(B)(t)、v2(AB), where "blue", "red" and "black" are corresponding to lengths of time segment of two hours, 1 day and 6 days respectively

例5  设观测数据r2(A)(t)、r2(B)(t)、v2(AB)以及作用于卫星的非惯性力具有误差,其误差服从均值为零的正态分布,各项标准差如下:对于位置矢量r2(A)(t)和r2(B)(t),标准差为3 cm、对于距离变化率v2(AB),标准差为0.6×10-6m·s-1、对于非惯性力,标准差为3×10-10m·s-2.需要说明的是,这里标准差是依据GRACE计划对观测量的误差指标设定的.然后按照例4的方法来恢复真实引力场,恢复的引力场的误差阶方差由图 5给出.

图 5 观测数据r2(A)(t)、r2(B)(t)、v2(AB)以及非惯性力存在误差、子弧段上参考轨道与真实轨道的初始条件相同时恢复的引力场的误差阶方差,这里“蓝线”表示子弧段为2 h的结果、“红线”表示子弧段为1天的结果、“黑线”表示弧段没有分段的结果 Fig. 5 The error degree-variance of the recovered gravity field when there are some errors in r2(A)(t)、r2(B)(t)、v2(AB), where "blue", "red" and "black" are corresponding to lengths of time segment of two hours, 1 day and 6 days respectively

综合图 4图 5可知,无论观测数据r2(A)(t)、r2(B)(t)、v2(AB)以及作用于卫星的非惯性力是否存在误差,是否以分段的形式对顾及非线性改正的轨道扰动方程求解几乎没有任何影响.这说明了顾及非线性改正的轨道扰动方程适合在较长的时间弧段上建立关于引力场位系数的法方程组.例如:例5的结论说明可以在6天的弧段上利用顾及非线性改正的轨道扰动方程来恢复引力场,这显然比目前针对线性化轨道扰动方程仅使用最多1天的弧段要长的多.事实上,在使用线性化轨道扰动方程时有将弧段取得更短的趋势,例如:Arsov和Pail (2003)曾认为将弧段取为4 h是合适的.此外,对比图 4图 5可知,图 4显示的恢复真实引力场的精度要远高于图 5的结果,其主要原因是在例5中对直接用于建立法方程组的v2(AB)项增加了误差.事实上,这也同时说明了若v2(AB)的误差不超过0.6×10-6m·s-1,那么恢复的引力场的误差阶方差将好于10-10.5≈3×10-11(至少60阶次之前是这样).另外,顺便指出,针对GRACE计划的模拟计算,并没有使用线性化轨道扰动方程来恢复引力场,其原因是:第一、前面3个算例已经证明了线性化轨道扰动方程的精度远低于顾及非线性改正的轨道扰动方程,而且线性化轨道扰动方程不适合在较长弧段上建立引力场位系数的法方程组;第二、进行例4和例5计算的原因就是要验证顾及非线性改正的轨道扰动方程适合在长弧段上建立法方程组.

最后对为何要将时间弧段进行延长说明一下.首先是相对于短弧段来说,在长弧段上利用轨道扰动方程来建立法方程组在算法要简单一些,至少关于初始值的未知数要少很多.其次更重要的是卫星经过长时间的运行,某些作用于卫星的弱力对轨道的作用将表现出来,也就是说,这些弱力经过长时间对卫星的作用,其作用的效果将会被观测数据反映出来,因此建立适合在长时间弧段上的轨道扰动方程有助于进一步分析除引力场之外的弱力.另外,通过本文的研究发现,使用动力学方法处理CHAMP和GRACE数据时非惯性力的测量精度εa仅影响积分弧长,也就是说,提高εa的测量精度的本质就是延长积分弧长.

5 结论

本文基于卫星运动满足的二阶微分方程,在引入参考轨道后针对CHAMP型卫星直接推导出了相应卫星轨道的线性化轨道扰动方程以及顾及非线性改正的轨道扰动方程.该推导方法是首次给出,具有推导形式简洁,便于理解的特征.

对推导出的线性化轨道扰动方程与顾及非线性改正的轨道扰动方程,分别估计了相应的误差.在定位精度为3cm、非惯性力测量精度3×10-10 m·s-2的情况下证明了:只有当参考轨道与实际轨道的距离ρ≤4.7 m时线性化轨道扰动方程的精度才能达到非惯性力的测量精度以及当ρ≤4.14×103 m时顾及非线性改正的轨道扰动方程就能够达到非惯性力的测量精度.由此可知,在同样的参考轨道下顾及非线性改正的轨道扰动方程比线性化轨道扰动方程具有更高的精度,同时这也意味着顾及非线性改正的轨道扰动方程能够在更长的时间弧段上建立关于引力场位系数的法方程组.

利用叠加原理,给出了轨道扰动方程的求解方法,从而完整地建立起了利用顾及非线性改正的轨道扰动方程建立位系数的理论和方法.

针对CHAMP计划进行了模拟计算,其结果表明:正常引力场的选择以及参考轨道初始条件的选择对顾及非线性改正的轨道扰动方程的精度影响较小,而对线性化轨道扰动方程的精度影响较大.其原因在于线性化轨道扰动方程仅适合较小的ρ值,而顾及非线性改正轨道扰动方程适合较大的ρ值.这也验证了顾及非线性改正的轨道扰动方程能够在更长的时间弧段上建立关于引力场位系数的法方程组这样的结论.

利用顾及非线性改正的轨道扰动方程对GRACE卫星计划进行了模拟解算,解算过程中采用了3种子弧段方式(2天、1天、6天)来分段建立位系数的法方程组,计算结果表明:这3种不同的弧段方式几乎对恢复引力场不产生影响.这也就验证了对于GRACE卫星计划而言,顾及非线性改正的轨道扰动方程依然适合在长的时间弧段(例如:6天)上建立关于位系数的法方程组.

随着GRACE-follow-on卫星计划的即将实施(Loomis et al., 2012; Sheard et al., 2012),非惯性力的测量精度将会极大地提高,预计该精度将达到10-12 m·s-2,更高就需要适合这样精度的轨道扰动方程.相比于线性化轨道扰动方程,显然顾及非线性改正的轨道扰动方程将会更适合用于处理GRACE-follow-on数据.

参考文献
Arsov K, Pail R. 2003. Assessment of two methods for gravity field recovery from GOCE GPS-SST orbit solutions. Adv. Geosci, 1: 121-126. DOI:10.5194/adgeo-1-121-2003
Beutler G, Jäggi A, Mervart L, et al. 2010. The celestial mechanics approach:theoretical foundations. J. Geod., 84(10): 605-624. DOI:10.1007/s00190-010-0401-7
Chen J L, Wilson C R, Tapley B D, et al. 2004. Low degree gravitational changes from GRACE:Validation and interpretation. Geophy. Res. Lett., 31(22): L22607. DOI:10.1029/2004GL021670
Cheng M K. 2002. Gravitational perturbation theory for intersatellite tracking. J. Geod., 76(3): 169-185. DOI:10.1007/s00190-001-0233-6
Cui C, Lelgemann D. 2000. On the non-linear low-low SST observation equations for the determination of the geopotential based on an analytical solution. J. Geod., 74(5): 431-440. DOI:10.1007/s001900000103
Ditmar P, van Eck van der Sluijs A A. 2004. A technique for modelling the Earth's gravity field on the basis of satellite accelerations. J. Geod., 78(1): 12-33.
Ditmar P, Kuznetsov V, van Eck van der Sluijs A A, et al. 2006. 'DEOS CHAMP-01C 70':a model of the Earth's gravity field computed from accelerations of the CHAMP satellite. J. Geod., 79(10-11): 586-601. DOI:10.1007/s00190-005-0008-6
Ditmar P, da Encarnação J T, Farahani H H. 2012. Understanding data noise in gravity field recovery on the basis of inter-satellite ranging measurements acquired by the satellite gravimetry mission GRACE. J. Geod., 86(6): 441-465. DOI:10.1007/s00190-011-0531-6
Douglas B C, Goad C C, Morrison F F. 1980. Determination of the geopotential from satellite-to-satellite tracking data. J. Geophy. Res., 85(B10): 5471-5480. DOI:10.1029/JB085iB10p05471
ESA SP-1233(1). 1999. Gravity field and steady-state ocean circulation mission. ESA Reports for Mission Selection.
Ilk K H, Feuchtinger M, Mayer-Gürr T. 2005. Gravity field recovery and validation by analysis of short arcs of a satellite-to-satellite tracking experiment as CHAMP and GRACE.//Sanso F. A Window on the Future of Geodesy. Berlin:Springer, 189-194.
Ilk K H, Löcher A, Mayer-Gürr T. 2008. Do we need new gravity field recovery techniques for the new gravity field satellites?.//Xu P. L, Liu J. N, Dermanis A. VI Hotine-Marussi Symposium on Theoretical and Computational Geodesy, Berlin:Springer, 3-8.
Jekeli C. 1999. The determination of gravitational potential differences from satellite-to-satellite tracking. Celest. Mech. Dyn. Astron., 75(2): 85-101. DOI:10.1023/A:1008313405488
Kaula, W M. 1966. Theory of Satellite Geodesy. London:Blaisdell Publishing Company.
Loomis B D, Nerem R S, Luthcke S B. 2012. Simulation study of a follow-on gravity mission to GRACE. J. Geod., 86: 319-335. DOI:10.1007/s00190-011-0521-8
Mayer-Gürr T, Ilk K H, Eicker A, et al. 2005. ITG-CHAMP01:a CHAMP gravity field model from short kinematic arcs over a one-year observation period. J. Geod., 78(7-8): 462-480. DOI:10.1007/s00190-004-0413-2
Montenbruck O, Gill E. 2000. Satellite Orbits. Berlin, Heidelberg:Springer.
Pavlis N K, Holmes S A, Kenyon S C, et al. 2012. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). J. Geophy. Res., Solid Earth, 117: B04406. DOI:10.1029/2011JB008916
Reigber C. 1989. Gravity field recovery from satellite tracking data.//Sansò F, Rummel R. Theory of Satellite Geodesy and Gravity Field Determination. Lecture Notes in Earth Sciences. Berlin, Heidelberg:Springer, 197-234.
Reigber C, Lühr H, Schwintzer P. 2002. CHAMP mission status. Adv. Space Res., 30(2): 129-134. DOI:10.1016/S0273-1177(02)00276-4
Reigber C, Balmino G, Schwintzer P, et al. 2003. Global gravity field recovery using solely GPS tracking and accelerometer data from CHAMP. Space Sci Rev., 108(1-2): 55-66.
Rowlands D D, Ray R D, Chinn D S, et al. 2002. Short-arc analysis of inter-satellite tracking data in a gravity mapping mission. J. Geod., 76(6): 307-316.
Rummel R, Balmino G, Johannessen J, et al. 2002. Dedicated gravity field missions-principles and aims. J. Geodyn., 33(1-2): 3-20. DOI:10.1016/S0264-3707(01)00050-3
Sheard B S, Heinzel G, Danzmann K, et al. 2012. Intersatellite laser ranging instrument for the GRACE follow-on mission. J. Geod., 86(12): 1083-1095. DOI:10.1007/s00190-012-0566-3
Taff L G. 1985. Celestial mechanics:A Computational Guide for the Practitioner. New York:Wiley-Interscience.
Tapley B D. 1989. Fundamentals of orbit determination.//Sansò F, Rummel R. Theory of Satellite Geodesy and Gravity Field Determination. Lecture notes in Earth Sciences. Berlin, Heidelberg:Springer, 235-260.
Tapley B D, Bettadpur S, Watkins M, et al. 2004a. The gravity recovery and climate experiment:mission overview and early results. Geophy. Res. Lett., 31: L09607. DOI:10.1029/2004GL019920
Tapley B D, Schutz B E, Born G H. 2004b. Statistical Orbit Determination. Burlington, San Diego, London:Elsevier Academic Press.
Visser P N A M, Sneeuw N, Gerlach C. 2003. Energy integral method for gravity field determination from satellite orbit coordinates. J. Geod., 77(3-4): 207-216. DOI:10.1007/s00190-003-0315-8
Wagner C, McAdoo D, Klokocník J, et al. 2006. Degradation of geopotential recovery from short repeat-cycle orbits:application to GRACE monthly fields. J. Geod., 80(2): 4-103.
Wolff M. 1969. Direct measurements of the Earth's gravitational potential using a satellite pair. J. Geophy. Res., 74(22): 5295-5300. DOI:10.1029/JB074i022p05295
Xu P. 2008. Position and velocity perturbations for the determination of geopotential from space geodetic measurements. Celest. Mech. Dyn. Astron., 100(3): 231-249. DOI:10.1007/s10569-008-9117-x