GLONASS卫星轨道积分算法分析

T ab. 2 Co mpar ison of Results Between Fix- Step
Integ ral and V ar-i Step Integr al/ m
两种积分方法坐 标分量差值
最大值
最小值
均值
x
0. 544
0. 076
0. 385
y
0. 864
0. 159
0. 521
z
0. 747
6 378 136 m; 重 力 位 第 二 带 谐 系 数 C20 =
- 0. 001 082 657; 地球自转速度 X= 7. 292 115 @ 10- 5 rad/ s; xdLS 、ydLS 、zdLS 为日 月摄动加速度之和, 在卫星星历中给出。
方程( 1) 为一个二阶常微分方程组, 且为隐函 数, 方程组中不含自变量 t。实际计算过程中, 可 采用定步长四阶龙格- 库塔法[ 5] 对其进行处理。
0. 135
0. 455
由表 2 的比较结果可以看到, 利用定步长积 分方法与变步长积分方法分别得到的相同时刻的 GL ONASS 卫星坐标 各分量 值的差 异均在 1 m 内, 可以认为这两种方法的精度是相当的。
90
30
x
7. 586 1. 866 5. 484
y
9. 622 2. 378 6. 369
z
11. 313 2. 745 6. 857
表 1 数据表明, 积分结果的误差随着积分时 间的增长而增大, 因此, 对于高精度定位而言, 最 好将积分时间控制在 60 min 以内。一般来讲, 可 以对星历中给出的每 个参考历元的 卫星状态向 前、向后分别积分 15 min, 则可得到整个观测时 段内的卫星位置和速度, 但在这种模式下, 由于定 步长积分得到的是各个均匀分布的节点上的卫星 状态, 而需要得到的卫星状态的时刻一般并不落 在某个节点上, 因此, 还需要利用多项式拟合。在 实际应用中, 可能会遇到星历中给出的个别卫星 的参考星历与测量时刻的时间间隔大于 30 min 的情况( 其他相近参考历元中没有给出此卫星星 历数据) , 因此, 轨道积分时, 还需单独对其增大积 分区间, 处理较为麻烦。
bc = f ( t , x , y , z , xc) ( 3)
cc = f ( t, x , y , z )
xc = a, yc = b, zc = c 则式( 2) 就化为一个一阶常微分方程组。
根据四阶龙格- 库塔格式[ 3] , 推导得该微分方
程组的积分公式为:
收稿日期: 2006-04-25。 项目来源: 国家基础测绘研究基金资助项目( 1469990324229) 。
hP 2 / 2, z ( ti ) + hQ 2 / 2, a( t i ) + hK 2 / 2) M3 = f ( ti + h/ 2, x ( t i ) + hN 2 / 2, y ( ti ) +
hP 2 / 2, z ( ti ) + hQ 2 / 2) N 3 = ai + hK 2 / 2 P3 = bi + hL 2 / 2 Q3 = ci + hM 2 / 2 K 4 = f ( ti + h, x ( ti ) + hN 3 , y( ti ) +
61 4
武 汉 大学 学报 # 信息 科 学版
wk.baidu.com
2006 年 7 月
ai+ 1 = ai + h( K 1 + 2K 2 + 2K 3 + K 4 ) / 6 bi+ 1 = bi + h( L 1 + 2L 2 + 2L 3 + L 4 ) / 6 ci+ 1 = ci + h( M1 + 2M2 + 2M3 + M 4 ) / 6 x i+ 1 = x i + h( N 1 + 2N 2 + 2N 3 + N 4 ) / 6 y i+ 1 = y i + h( P1 + 2P 2 + 2P 3 + P4 ) / 6 z i+ 1 = z i + h( Q1 + 2Q 2 + 2Q 3 + Q4 ) / 6
利用变步长法积分则可以直接得到测量时刻 的卫星位置和速度, 同时可以免去多项式拟合的 步骤, 且能达到与定步长方法相当的精度, 在编程
第 31 卷第 7 期
杨 剑等: G LO N ASS 卫星轨道积分算 法分析
61 5
实现 上更为简单方便, 有利 于实现自动化处理。 变步长积分的算法具体见文献[ 3] 。
hP 3 , z ( t i ) + hQ 3 , b( ti ) + hK 3 ) L 4 = f ( t i + h, x ( t i ) + hN 3 , y ( ti ) +
hP 3 , z ( t i ) + hQ 3 , a( t i ) + hK 3 ) M4 = f ( ti + h, x ( ti ) + hN 3 , y ( ti ) +
表 1 定步长积分法结果与 星历所给坐标之差的绝对值 T ab. 1 Absolute V alue of D ifference Between Fix- Step
Integ ral R esult and Coor dinates Ephemer is Offer ing
积分时 间/ min
30
zd = -
GM r3
z
+
3 2
C20
#
G
Ma r5
2 e
#z
#
1-
5z 2 r2
+ zdLS
( 1)
式中, x 、y 、z 表示地固系中的卫星坐标, 是时间 t 的函数; xc、yc、zc为坐标对时间的一阶导数, 表示 卫星的速度; xd、yd、zd表示卫星的加速度, 为坐标 对时 间 的 二 阶 导 数。 地 球 引 力 常 数 GM = 398 600 440 m 3 / s2 ; 参 考 椭 球 长 半 径 ae =
( 4) 其中,
K 1 = f ( ti , x ( t i ) , y ( t i ) , z ( ti ) , b( ti ) ) L 1 = f ( t i , x ( t i ) , y ( ti ) , z ( ti ) , a( ti ) ) M1 = f ( ti , x ( ti ) , y ( ti ) , z ( ti ) ) N 1 = ai P1 = bi Q1 = ci K 2 = f ( ti + h/ 2, x ( t i ) + hN 1 / 2, y( ti ) +
首先将方程( 1) 写为函数的形式为:
xd = f ( f , x , y , z , yc)
yd = f ( t , x , y, z , xc)
( 2)
zd = f ( t, x , y , z )
对式( 2) 降阶, 令 xc= a, yc= b, zc= c, 则有:
ac = f ( t, x , y , z , yc)
积分步 长/ s
30
坐标分量 差值
x y z
最大 值/ m
0. 698 1. 684 1. 188
最小 值/ m
0. 081 0. 682 0. 497
均值/ m
0. 395 1. 104 0. 849
60
30
x
3. 082 1. 529 2. 115
y
4. 840 0. 069 2. 550
z
4. 279 0. 576 2. 226
利用四阶龙格- 库塔方法进行轨道积分时, 可 采用定步长方法或变步长方法。定步长方法具有 积分公式简单、容易编程实现和运算速度快等优 点。但由于积分误差是累积的, 因而积分步长和积 分区间的有效性是影响积分精度的主要因素。实 验表明, 利用参考历元的星历分别向前、向后积分 15 min, 积分步长取 0. 1 s、1 s 和 30 s, 其精度水平 几乎相同[ 4] , 在 60 min 范围内精度相差也不大, 因 此, 步长取 30 s 可极大地减小计算量, 同时保证一 定的精度水平。对于积分区间的选择对积分精度 的影响, 本文利用 30 min 间隔的卫星星历分别积 分 30 min、60 min、90 m in, 得到的积分结果与星 历中给出的坐标的差值绝对值如表 1 所示。
hP 3 , z ( t i ) + hQ 3 ) N 4 = ai + hK 2 P4 = bi + hL 2 Q4 = ci + hM 2 式中, h 为积分步长。 根据以上公式, 则可由卫星星历中给出的参 考历元的各卫星位置、速度以及日月摄动加速度 进行轨道积分, 得到任意时刻的卫星状态。
2 积分步长、积分区间的选择
hP 1 / 2, z ( ti ) + hQ 1 / 2) N 2 = ai + hK 1 / 2 P2 = bi + hL 1 / 2 Q2 = ci + hM 1 / 2 K 3 = f ( ti + h/ 2, x ( t i ) + hN 2 / 2, y( ti ) +
hP 2 / 2, z ( ti ) + hQ 2 / 2, b( ti ) + hL 2 / 2) L 3 = f ( t i + h/ 2, x ( ti ) + hN 2 / 2, y ( ti ) +
文献标志码: A
GLONASS 卫星轨道积分算法分析
杨 剑1 王泽民1 孟 泱1 委民正2
( 1 武汉大学测绘学院, 武汉市珞喻路 129 号, 430079) ( 2 福建省地质测绘院, 福州市塔头路 2 号, 350011)
摘 要: 推导了利用四阶龙格- 库塔方法对 GL ON A SS 卫星运动方程进行轨道积分的计算公式, 根据实际数据 处理 结果对积分区间、积分步长作了一定分析, 在对定步长积 分方法和变 步长积分 方法分别进 行介绍的 基础 上, 对两者的结果作了比较。 关键词: GL O N ASS; 轨道 积分; 龙格- 库塔法; 积分步长 中图法分类号: P228. 41
第 31 卷 第 7 期 2006 年 7 月
武汉大 学学报 # 信息科学版 Geo matics and Informat ion Science of W uhan U niver sity
V ol. 31 N o. 7 July 2006
文章编号: 1671- 8860( 2006) 07-0613- 03
1 轨道积分的计算方法
地固系中卫星的运动方程为[ 1, 4] :
xd = -
GM r3
x
+
3 2
C20
#
G
M
a
2 e
r5
#x
#
1-
5z 2 r2
+ xdLS+ X2 x + 2X# yc
yd = -
GM r3
y
+
3 2
C 20
#
GM a r5
2 e
#
y
#
1-
5z 2 r2
+ ydL S + X2 y + 2 X# xc
本文分别利用了定步长积分法和变步长积分 法求解了某测量时刻的卫星状态, 并对积分结果 进行了比较, 其中定步长方法采用 30 s 积分步长 以及 8 阶多项式拟合, 积分区间为 30 m in; 变步 长方法在步长减半 4 到 5 次左右达到收敛。两种 方法积分得到的结果的比较如表 2 所示。
表 2 定步长积分和变步长积分的结果比较/ m
对 GLONASS 卫星轨道进行数值积分可用 的方法有欧拉( Euler) 法、龙格- 库塔( Rung- Kutt a) 法, 阿达姆斯 ( Adams) 法等, 其中以龙格- 库 塔法较为 常用。葛茂荣 等[ 1] 推导了 GL O N A SS 卫星在地固系中的运动方程, 并简要介绍了利用 龙格- 库塔法进行轨道积分的数值结果; 张永军[ 2] 和李建文[ 3] 分别利用定步长和变步长龙格- 库塔 法积分实现了 GL ONASS 卫星轨道求 解。但以 上研究均未给出 GL ONA SS 卫星运 动方程的具 体处理步骤及公式。本文在深入研究数值积分方 法的基础上, 详细推导了 GL ONA SS 卫星运动二 阶常微分方程组的处理公式, 并分别利用定步长 和变步长龙格- 库塔数值积分方法对 GL ONASS 轨道积分进行了比较分析。
hP 1 / 2, z ( ti ) + hQ 1 / 2, b( ti ) + hL 1 / 2) L 2 = f ( t i + h/ 2, x ( ti ) + hN 1 / 2, y ( ti ) +
hP 1 / 2, z ( ti ) + hQ 1 / 2, a( t i ) + hK 1 / 2) M2 = f ( ti + h/ 2, x ( t i ) + hN 1 / 2, y ( ti ) +