精密定轨中变分方程的数值解法及其检验

洪 樱,欧吉坤

(1.武汉纺织大学理学院,湖北武汉 430073;2.中国科学院测量与地球物理研究所,湖北武汉 430077)

精密定轨中变分方程的数值解法及其检验

洪 樱1,欧吉坤2

(1.武汉纺织大学理学院,湖北武汉 430073;2.中国科学院测量与地球物理研究所,湖北武汉 430077)

变分方程的求解是卫星精密定轨中的关键步骤之一,通常都通过积分方法求解。首先从原理上详细阐述一种新的数值方法求解变分方程,接着通过低轨卫星精密定轨仿真软件对其精度和效率进行测试,且与常用数值积分方法的定轨结果进行比较。结果表明,在实际轨道确定中,采用新的数值方法可达到较高的精度要求,且算法的效率更高。

变分方程;数值方法;精密定轨;积分方法

一、引 言

在精密轨道确定中,如果知道初始状态矢量以及力矢量,则可通过积分卫星的运动方程得到精密轨道。在求解初始状态参数时,需用目前状态相对于初始状态矢量的偏导数,求解变分方程可得到构成状态转移矩阵的这种偏导数。所以,在积分运动方程的同时,必须求解变分方程[1-2]。对于变分方程的求解,许多学者都进行了研究和探讨,提出了许多有效的方法,如简化方法[3]、分析方法[3-4]、差分算法[5-6]、数值方法[2]等。其中数值方法是德国GFZ(GeoForschungs Zentrum)地学中心许国昌博士提出的一种新的一般化的求解变分方程的方法,其特点是通过差分代替微分使得变分方程转化为线性代数方程组,进而根据其初值条件推出严格的递推解。虽然文献[2]中给出了数值方法严格的推导,但从未对该方法进行过数值检验和比较。使得该方法的应用和推广缺少必要的数值依据。

本文在此理论基础之上,部分地应理论原创者的要求,通过详细的数值计算对这种新的求解变分方程的数值方法进行了分析。首先从原理上详细阐述了这种新的数值方法的求解过程,然后通过精密定轨仿真软件对其进行测试,且与常用的数值积分方法的结果进行了比较,从而证明了新解法的可行性和有益性。

二、新的数值解法

在轨道确定和预报过程中,需解算运动方程[1-2]

式中,˙X(t)为卫星的瞬时状态向量;X(t0)为 t0时刻的初始状态向量(记为 X0);F为时间 t、状态向量 ˙X(t)以及动力学参数向量 y的函数,且

r和 ˙r分别为卫星的位置和速度向量;f为卫星所受力的向量和;m为卫星质量;Y是除位置、速度以外的动力学参数向量。

在运动方程(1)两边对 y求偏导数有

式中,＊号表示对 F中显含的参数矢量 y的偏导,且记

式中,E为单位阵,0为零矩阵,下标是其维数。则式(2)变为

其解构成了状态转移矩阵,记为 Φ(t,t0)。则式(4)变为

式(5)称为变分方程。

初值条件为

由于微分方程列与列之间是相互独立的,因此只需考虑某列的方程的解。对于列 j,式(7)和初值条件式(8)可表示为

将式 (6)和式(3)代入式(5)得

若在时间间隔 [t0,t]中,步长 h=(t-t0)/m,则 tn=t0+nh,n=1,…,m,且

则式(9)变为

整理式(11)得

式中

则可得到 ψ(tn+1),速度向量 ψ˙(tn+1)可利用式(10)计算得到。

三、数值检验算例

为了验证数值方法的可靠性,以CHAMP卫星为例,在仿真过程中采用的动力学模型包括：70×70阶的 JG M3重力场模型、DT M78大气阻力模型、Wahr潮汐模型、N体摄动模型(行星位置由 JPL提供的行星历表中精确地插值获得)、光压摄动模型(假设卫星帆板的法线方向始终与太阳垂直)以及广义相对论摄动模型[7-9]。分别采用积分方法和数值方法计算状态转移矩阵,历元间隔取为 10 s,其中积分方法采用稳定性较好、精度较高的 RKF7(8)法起步、10阶的第 II种类型的 collocation-PECE算法[8]。由于低轨卫星的偏心率较小,所以可采用固定步长进行计算[1]。取初始时刻为 2002年 1月 3日 0时 0分0秒,轨道初始状态(惯性系下位置和速度)为x0=701 111.972 7 m, y0=-368 100.205 1 m, z0=6 723 471.518 4m, vx0=-4 695.158 9 m/s, vy0=6 017.559 8m/s, vz0=838.240 2m/s。为了避免比较数据过于庞大,这里仅列出状态转移矩阵的部分矩阵元,并且最少保留到数值方法与积分方法有差异的位数。将积分方法的计算结果记为Φ1,数值方法的计算结果记为Φ2,计算结果列于表 1。

表 1 积分方法和数值方法的计算结果比较

首先需要指出的是,随着积分弧段的增加,积分方法所花的时间迅速增长。在此算例计算中,积分方法和数值方法均积分了 1条轨道 (6个状态方程),并且都需计算加速度的偏导,而积分方法比数值方法多积分 6×6=36个变分方程。因此,积分方法共积分了 42个方程,而数值方法只积分了 6个方程。当再增加一个动力学参数时,积分方法需积分6+7×7=55个方程,而数值方法仍仅需积分 6个方程,所以,待估动力学参数越多,数值方法节省时间的优越性就越明显。当然其前提条件是要保证计算的精度。

表 1列出了积分方法和数值方法计算的状态转移矩阵中的部分矩阵元。从表中可以看出,当弧段较短时,两种方法计算出的状态转移矩阵中矩阵元的值相差不大,对于数值较大的矩阵元仍有 3～4位数字相同。但随着积分弧段的增加,两种方法的计算结果差异变大。

由于状态转移矩阵描述的是初始状态和模型参数的变化对动力学系统演化的影响,即变分方程的解——状态转移矩阵是无量纲的,无法直接描述其精度。而在求解初始状态参数时,需要用到状态转移矩阵,因此状态转移矩阵的精度会影响到最终的定轨结果。所以,不妨以精密轨道确定的结果来考察变分方程解法的效果。同样采用上述仿真软件并模拟 100min(略大于低轨卫星的运行周期)的P3码伪距观测值用于轨道确定,历元间隔为 10 s。分别采用数值法和积分法求解变分方程,利用批处理算法得到的轨道初始状态的改正量列于表 2中。

表 2 两种方法求解变分方程得到的轨道初始状态的改正量

表 2中方法A表示采用数值法求解变分方程,方法B表示采用积分法求解变分方程;OS0表示无误差的 P3码观测值,OS1表示只包含接收机钟差的P3码观测值,OS2表示包含接收机钟差、GPS卫星精密钟差和观测偶然误差的 P3码观测值。

从表 2可以看出,在相同观测条件下,用数值法和积分法求解变分方程得到的轨道初始状态的改正量相差甚小,坐标分量达到微米的数量级,而速度分量达到 10-8～10-10m/s的数量级。

四、结束语

本文首次验证了变分方程的数值解算方法,从而确认了定轨中一种全新的算法流程,即不再需要通过繁复的积分方法来求解变分方程。该方法并不像其他方法那样受外部条件的限制,如轨道弧长、动力学模型的复杂程度等;且算法的推导严格,程序易于实现。从低轨卫星轨道确定结果来看,在相同 GPS观测条件下采用数值法求解变分方程的精度与积分法求解变分方程的精度相当。由于积分方法的精度较高,可以推论这种新的数值方法的精度也是较高的。此外,从算法的效率来看,积分方法比数值方法多积分 n2(n为待估参数个数)个变分方程。因此,这种新的数值方法在实际轨道计算中更有效。

[1] 李济生.人造卫星精密轨道确定[M].北京：解放军出版社,1995：189-219.

[2] XU GC.GPS Theory,Algorithms andApplications[M]. Berlin：Springer-Verlag,2003：217-277.

[3] 刘林,张强,廖新浩.人卫精密定轨中的算法问题[J].中国科学：A辑,1998,28(9)：848-856.

[4] 张强,刘林.轨道改进中计算状态转移矩阵的分析方法[J].天文学报,1999,40(2)：113-121.

[5] 张家祥,汪琦,杨捷兴,等.彗木相撞预报[J].中国科学：A辑,1996,26(3)：280-288.

[6] 胡小工,黄辉,廖新浩.状态转移矩阵的差分算法及其应用[J].天文学报,2000,41(2)：113-122.

[7] 徐天河,杨元喜.基于能量守恒方法恢复 CHAMP重力场模型[J].测绘学报,2005,25(2)：75-81.

[8] MONTENBRUCK O,G ILL E.Satellite Orbits-Models, Methods andApplications[M].New York：Springer-Verlag,2000：53-154.

[9] BOCK H.EfficientMethods for Deter mining Precise Orbits of Low Earth Orbiters Using the Global Positioning System[D].Switzerland：Institute für Geodäsie und Photogrammetrie Eidg,2003：77-85.

Numerical Solution of Var iational Equation for Precise OrbitDeterm ination and Its Test

HONG Ying,OU Jikun

0494-0911(2010)12-0001-03

P228

B

2010-07-20

洪 樱(1980—),女,湖北荆门人,硕士,讲师,主要从事卫星精密轨道确定与预报理论研究。