首页资源分类应用技术测试测量 > 卫星轨道动力学数值计算

卫星轨道动力学数值计算

已有 445110个资源

下载专区

文档信息举报收藏

标    签:卫星轨道数值计算

分    享:

文档简介

卫星轨道动力学的数值计算

文档预览

目录 1 星历计算的时间和坐标系统........................................................................................................2 1.1 有关的时间系统与坐标系统............................................................................................2 1.1.1 时间系统及其换算.................................................................................................2 1.1.2 坐标系统及其换算.................................................................................................4 1.2 计算单位和有关常数........................................................................................................7 2 轨道动力学计算的基本数学模型.............................................................................................12 2.1 二体问题..................................................................................................................12 2.2 地球非球形引力摄动..............................................................................................12 2.3 日、月摄动..............................................................................................................15 2.4 太阳直接辐射压摄动..............................................................................................16 2.5 地球固体潮摄动......................................................................................................19 2.6 大气阻力摄动..........................................................................................................19 2.7 Y 轴偏差加速度摄动...............................................................................................20 2.8 巡航姿态控制动力摄动..........................................................................................20 2.9 其它摄动影响..........................................................................................................21 附录:日月位置计算..................................................................................................... 21 3 轨道计算方法 ............................................................................................................................ 24 3.1 Runge_Kutta 积分法 ........................................................................................................24 3.2 Adams_Cowell 积分.........................................................................................................25 3.3 轨道计算.......................................................................................................................... 27 3.4 星历的快速插值..............................................................................................................28 4 轨道根数与位置矢量、速度矢量的关系.................................................................................32 4.1 由位置矢量和速度矢量计算轨道根数..........................................................................32 4.2 由轨道根数计算位置矢量和速度矢量..........................................................................33 1 1 星历计算的时间和坐标系统 1.1 有关的时间系统与坐标系统 轨道计算过程重要涉及到不同的时间系统和坐标系统,下面将空间战场环境系统中所涉 及到的时间系统和坐标系统进行定义,并说明各系统之间的相互关系。一般情况下,仿真系 统采用的是 TDT 时间系统和 J2000 地心惯性坐标系。 1.1.1 时间系统及其换算 在轨道计算中,时间是独立变量。但是,在计算不同的物理量时,却使用不同的时间系 统。例如:在计算恒星时用世界时 UT1;定位解算时采用 GPS 时 GPST;岁差和章动量的计算 采用 TDB 时等。所以必须清楚各时间系统的定义和各时间系统之间的转换,下面给出各种时 间系统的定义及它们之间的转换公式。 格林尼治恒星时 格林尼治恒星时为春分点对格林尼治平天文子午面的时角。由于岁差、章动原因,它由 格林尼治真恒星时(GAST)和平恒星时(GMST)之分。 两者的关系是: GAST  GMST   cos 其中:  cos 为赤经章动 GMST  67310s.54841 (8640184s.812866  876600h )Tu  0s.093104Tu2  0s.62 105Tu3 Tu 为自 J 2000.0(JD2451545.0) 起算至观测UT1时刻的儒略世纪数,即 Tu  JD(UT1)  2451545.0 36525.0 世界时UT1 UT1是以平北极(国际习惯用原点)为统一标准的观测世界时,是反映地球实际自转 的时间,恒星时计算与此有关。 国际原子时 TAI TAI 时以铯原子 C 133 s 基态两能级间跃迁辐射的 9192631770 周所经历的时间作为 1 秒长 的均匀时间,起点在 1958 年 1 月 1 日 0hUT 。 国际协调时UTC 2 UTC 是经跳秒修改后的国际原子时,它与世界时UT1的差 0s.9 ,观测纪录时间是以此 为准的。 质心动力学时TDB (Barycentric Dynamical Time) TDB 为相对于太阳质心的运动方程给出的历表、引数等所用的时间尺度,岁差及章动 量的计算是以此为依据的。 地球动力学时TDT (Terrestrial Dynamical Time) TDT 为视地心历表所用的时间尺度,它具有均匀连续的特性,卫星运动方程就是以此 为独立的时间变量。 GPS 时间 GPST GPST 是由系统定义和应用的一种时间尺度,于 1980 年 1 月 6 日 0h GPST 与 UTC 相等, 在此以后由系统主控站密切跟踪 UTC 以保持高度统一。但 GPST 不作跳秒修正,因此它与 UTC 具有整秒的差异(1997 年 1 月至 6 月相差为11s )。在计算 GPS 卫星轨道的初值时将涉及到 GPST,GPS 精密星历的参考时间为 GPST。 以上各时间尺度的相互关系如下: UT1  UTC  UT1 TAI  UTC  AT TDT  TAI  32s.184 TDB  TDT  TD GPST  UTC  GPST 其中: UT1可从地球自转参数文件中获得; AT  19s  GPST ; TD  0s.001658sin(v  0.0167sin v),v  6.240040768  628.3019501T (rad) 。 上式中的 T 为自 J 2000.0 年起算至观测 TDB 时刻的儒略世纪数,即: T  JD(TDB)  2451545.0 36525.0 不同时间系统间的关系如下图所示: 32.184s 动力学时(TDT) 原子时(TAI) 0s GPS 时(GPST) -19s 协调世界时(UTC) 世界时(UT1) 图 1 几种时间系统之间的关系 3 1.1.2 坐标系统及其换算 卫星轨道计算和实际定位解算分别是在 J2000.0 惯性坐标系与 WGS-84 地固系中进行的, 此外,卫星加速度计算中还涉及到星固坐标系。下面给出与本课题有关的主要坐标系的定义 及相互转换关系。 J 2000.0 地心惯性系 原点:地球质心 Z 轴:向北指向 J 2000.0 年平赤道面(基面)的极点 X 轴:指向 J 2000.0 平春分点 Y 轴:符合右手系法则  位置矢量: r 星固坐标系 原点:卫星质心 Z 轴:指向卫星的天线方向,即指向地心 X 轴:在轴与太阳构成的平面内,完成右手系法则 Y 轴:沿卫星太阳能翼的支轴 位置矢量: ra 星固坐标系坐标轴 ( Xˆ a , Yˆa , Zˆ a ) 在惯性系中的方向余弦分别为:(  rs ,  r 分别为太阳和卫 星在 J 2000.0 地心惯性系中的位置矢YXˆˆ量aZaˆa)rrssrrYYrraa  ,s   rs  r WGS-84 坐标系 WGS-84 为 1984 年世界大地坐标系(WGS 为 World Geodetic System 的简称),WGS-84 的坐标定义及其采用的椭球参数为: 原点:地球质心 Z 轴:指向 BIH1984.0 定义的协议地球极(CTP)方向 X 轴:指向 BIH1984.0 的零子午面和 CTP 赤道的交点 Y 轴:与 X、Z 轴成右手系 地球椭球长半径:ae=6378137 m 地球引力常数(含大气层): GM=3986005108 m3/s2 4 正常化二阶带球谐系数: c2.0 =-484.1668510-6 地球自转角速度:=729211510-11 rad/s2 地球椭球扁率:f = 1/298.257223563 地固坐标系与惯性坐标系的转换模型 转换到 a. 惯性坐标系 地固坐标系 模型 X  X  X    y    [ A][B][C][D]Y    [W ]Y   Z  D Z  J 2000.0 Z  J 2000.0 转换到 b. 地固坐标系 惯性坐标系 模型 X  X  X    y    [D]T [C]T [B]T [ A]T Y    [W ]T Y   Z  J 2000.0 Z  D Z  D 式中:[A]为极移矩阵; [B]为自转矩阵; [C]为章动矩阵; [D]为岁差矩阵。 上述各矩阵的意义及具体定义如下: 极移:由于地球不是刚体及其它一些地球物理因素的影响,地球自转轴相对于地球的位置随 时间而变化从而引起观察者的天顶在天球上的位置发生变化,称为极移,矩阵为[A]: [ A]  RY (x p )RX ( y p ) 式中: x p , y p 为地极坐标,可从地球自转参数文件中给出的极移值插值得到。 自转:即地球公转的同时也在绕自转轴旋转。矩阵[B]: [B]  RZ (G ) 式中:G 为格林尼治恒星时  G  67310s.54841 (8640184s.812866  876600h )Tu  0s.093104Tu2  0s.62 105Tu3   cos( M   ) JD(UT1)  2451545.0 Tu  36525.0 章动:是指外力作用下,地球自转轴在空间运动的短周期摆动部分,即同一瞬间真天极相对 平天极的运动,月球对地球引力的变化是形成章动现象的主要外力作用,其次是太阳。矩阵 [C]: [C]  RX ( M   )RZ ()RX ( M ) 式中:  M  2326'21''.448  46''.8150T0  0''.00059T02  0''.001813T03 5 106 5     d j cos( n jk Ak ) ——交角章动 j 1 k 1 106 5     c j sin( n jk Ak ) ——黄经章动 j 1 k 1 其中: c j , d j , n jk 都为常数,可自章动系数表 1 中查出。 而月亮平近点角: A1  485866''.733  (1325r  715922''.633)T0  31''.310T02  0''.064T03 太阳平近点角: A2  1287099''.804  (99r  1292581''.224)T0  0''.577T02  0''.012T03 月亮平升交角: A3  335778''.877  (1342r  295263'.137)T0 13''.257T02  0''.011T03 日月平角矩: A4  1072261''.307  (1236r 1105601''.328)T0  6''.891T02  0''.019T03 月亮轨道对黄道平均升交点的黄经: A5  450160''.280  (5r  482890''.539)T0  7''.455T02  0''.008T03 其中:1r  360 T0 为自 J2000.0 起算至 t 的儒略世纪数 T0  MJD(TDB)  51544.5 36525.0 岁差:地球在太阳、月球和行星的引力作用下,地球的自转轴在空间不断发生变化,其长期 运动称为岁差,矩阵[D]: [D]  RZ (Z p )RY ( p )RZ ( p ) 式中:  p  2306''.2181T0  0''.30188T02  0''.017998T03  p  2004''.3109T0  0''.42665T02  0''.041833T03 Z p  2306''.2181T0 1''.09468T02  0''.018203T03 T0 的意义同上。 旋转矩阵的求法分别为 1 RX ( )  0 0 0 cos  sin  0 sin    cos  cos 0  sin   RY ( )    0 1 0   sin  0 cos   cos sin  0 RZ ( )   sin  cos 0  0 0 1 6 1.2 计算单位和有关常数 轨道计算采用人卫单位系统,具体定义为: 长度单位:地球赤道半径 ae 质量单位:地球的总质量 M e 1 时间单位: Tas   ae3 GM e 2    8.068111241279087 102 在以上人卫单位系统中,引力常数 G=1。为完整起见,给出以下常数: 地球赤道半径:6378137m 地球扁率: f  1 298.257 地球总质量: M e  5.974 1024 kg 地球自转角速度:e  7.292115105 rad / s 地心引力常数: GM e  3.9860051014 m3 / s 2 日心引力常数: GM s  1.327124381020 m3 / s 2 天文单位长度: asun  1.4959787 1011 m 月球地球质量比: GM L GMe  0.01230002 光速: c  299792458m / s 引力常数: G  6.6721011m3 /(kg  s2 ) 太阳常数: Ps  4.5605106 kg /(m  s 2 ) 地球引力位系数采用 WGS-84 EGM 的规格化值。参见表 2。 表 1 黄经和倾角章动序列表 引数 i L L’ F D  1 00001 2 00002 3 -2 0 2 0 1 4 2 0 -2 0 0 5 -2 0 2 0 2 6 1 -1 0 -1 0 7 0 -2 2 -2 1 8 2 0 -2 0 1 9 0 0 2 -2 2 10 0 1 0 0 0 11 0 1 2 -2 2 黄经章动 0".0001 0".0001T -171996 -174.2 2062 0.2 46 0 11 0 -3 0 -3 0 -2 0 1 0 -13187 -1.6 1426 -3.4 -517 1.2 倾角章动 0".0001 0".0001T 92025 8.9 -895 0.5 -24 0 0 0 1 0 0 0 1 0 0 0 5736 -3.1 54 -0.1 224 -0.6 7 12 0 -1 2 -2 2 13 0 0 2 -2 1 14 2 0 0 -2 0 15 0 0 2 -2 0 16 0 2 0 0 0 17 0 1 0 0 1 18 0 2 2 -2 2 19 0 -1 0 0 1 20 -2 0 0 2 1 21 0 -1 2 -2 1 22 2 0 0 -2 1 23 0 1 2 -2 1 24 1 0 0 -1 0 25 2 1 0 -2 0 26 0 0 -2 2 1 27 0 1 -2 2 0 28 0 1 0 0 2 29 -1 0 0 1 1 30 0 1 2 -2 0 31 0 0 2 0 2 32 1 0 0 0 0 33 0 0 2 0 1 34 1 0 2 0 2 35 1 0 0 -2 0 36 -1 0 2 0 2 37 0 0 0 2 0 38 1 0 0 0 1 39 -1 0 0 0 1 40 -1 0 2 2 2 41 1 0 2 0 1 42 0 0 2 2 2 43 2 0 0 0 0 44 1 0 2 -2 2 45 2 0 2 0 2 46 0 0 2 0 0 47 -1 0 2 0 1 48 -1 0 0 2 1 49 1 0 0 -2 1 50 -1 0 2 2 1 51 1 1 0 -2 0 52 0 1 2 0 2 217 129 48 -22 17 -15 -16 -12 -6 -5 4 4 -4 1 1 -1 1 1 -1 -2274 712 -386 -301 -158 123 63 63 -58 -59 -51 -38 29 29 -31 26 21 16 -13 -10 -7 7 -0.5 -95 0.3 0.1 -70 0 0 1 0 0 0 0 -0.1 0 0 0 9 0 0.1 7 0 0 6 0 0 3 0 0 3 0 0 -2 0 0 -2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 977 -0.5 0.1 -7 0 -0.4 200 0 0 129 -0.1 0 -1 0 0 -53 0 0 -2 0 0.1 -33 0 -0.1 32 0 0 26 0 0 27 0 0 16 0 0 -1 0 0 -12 0 0 13 0 0 -1 0 0 -10 0 0 -8 0 0 7 0 0 5 0 0 0 0 0 -3 0 8 53 0 -1 2 0 2 54 1 0 2 2 2 55 1 0 0 2 0 56 2 0 2 -2 2 57 0 0 0 2 1 58 0 0 2 2 1 59 1 0 2 -2 1 60 0 0 0 -2 1 61 1 -1 0 0 0 62 2 0 2 0 1 63 0 1 0 -2 0 64 1 0 -2 0 0 65 0 0 0 1 0 66 1 1 0 0 0 67 1 0 2 0 0 68 1 -1 2 0 2 69 -1 -1 2 2 2 70 -2 0 0 0 1 71 3 0 2 0 2 72 0 -1 2 2 2 73 1 1 2 0 2 74 -1 0 2 -2 1 75 2 0 0 0 1 76 1 0 0 0 2 77 3 0 0 0 0 78 0 0 2 1 2 79 -1 0 0 0 2 80 1 0 0 -4 0 81 -2 0 2 2 2 82 -1 0 2 4 2 83 2 0 0 -4 0 84 1 1 2 -2 2 85 1 0 2 2 1 86 -2 0 2 4 2 87 -1 0 4 0 2 88 1 -1 0 -2 0 89 2 0 2 -2 1 90 2 0 2 2 2 91 1 0 0 2 1 92 0 0 4 -2 2 93 3 0 2 -2 2 -7 0 3 0 -8 0 3 0 6 0 0 0 6 0 -3 0 -6 0 3 0 -7 0 3 0 6 0 -3 0 -5 0 3 0 5 0 0 0 -5 0 3 0 -4 0 0 0 4 0 0 0 -4 0 0 0 -3 0 0 0 3 0 0 0 -3 0 1 0 -3 0 1 0 -2 0 1 0 -3 0 1 0 -3 0 1 0 2 0 -1 0 -2 0 1 0 2 0 -1 0 -2 0 1 0 2 0 0 0 2 0 -1 0 1 0 -1 0 -1 0 0 0 1 0 -1 0 -2 0 1 0 -1 0 0 0 1 0 -1 0 -1 0 1 0 -1 0 1 0 1 0 0 0 1 0 0 0 1 0 -1 0 -1 0 0 0 -1 0 0 0 1 0 0 0 1 0 0 0 9 94 1 0 2 -2 0 95 0 1 2 0 1 96 -1 -1 0 2 1 97 0 0 -2 0 1 98 0 0 2 -1 2 99 0 1 0 2 0 100 1 0 -2 -2 0 101 0 -1 2 0 1 102 1 1 0 -2 1 103 1 0 -2 2 0 104 2 0 0 2 0 105 0 0 2 4 2 106 0 1 0 1 0 -1 0 0 0 1 0 0 0 1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 1 0 0 0 -1 0 0 0 1 0 0 0 表 2 计算地球引力位加速度的引力位系数(GEM-T3) n m 2 0 2 1 2 2 3 0 3 1 3 2 3 3 4 0 4 1 4 2 4 3 4 4 5 0 5 1 5 2 5 3 5 4 5 5 6 0 6 1 6 2 6 3 6 4 6 5 6 6 7 0 7 1 7 2 7 3 7 4 7 5 7 6 7 7 8 0 8 1 Cnm -0.48416510e-3 -0.17e-9 0.24390658e-5 0.95720109e-6 0.20277142e-5 0.90447073e-6 0.72034249e-6 0.53952118e-6 -0.53615108e-6 0.35021814e-6 0.99093372e-6 -0.18877065e-6 0.68343345e-7 -0.58280231e-7 0.65271099e-6 -0.45233006e-6 -0.29558409e-6 0.17376349e-6 -0.14951352e-6 -0.76894161e-7 0.48734482e-7 0.57203194e-7 -0.86826474e-7 -0.26733038e-6 0.96846267e-7 0.91300885e-7 0.27486873e-6 0.32779498e-6 0.25122012e-6 -0.27556101e-6 0.13261698e-8 -0.35883140e-6 0.97027653e-9 0.48883181e-7 0.23628239e-7 S nm 0.0 1.19e-9 -0.14000946e-5 0.0 0.24921712e-6 -0.61944767e-6 0.14138845e-5 0.0 -0.47343598e-6 0.66301523e-6 -0.20092742e-6 0.30942370e-6 0.0 -0.96083937e-7 -0.32386369e-6 -0.21529578e-6 0.49690346e-7 -0.66890704e-6 0.0 0.26998384e-7 -0.37401311e-6 0.93727543e-8 -0.47130637e-6 -0.53678021e-6 -0.23713482e-6 0.0 0.97465920e-7 0.93246712e-7 -0.21529269e-6 -0.12376718e-6 0.18620005e-7 0.15173875e-6 0.24083642e-7 0.0 0.58847236e-7 10 8 2 8 3 8 4 8 5 8 6 8 7 8 8 0.77598534e-7 -0.17785247e-7 -0.24633984e-6 -0.25041147e-7 -0.64923712e-7 0.67462214e-7 -0.12419836e-6 0.66008696e-7 -0.86346983e-7 0.70179584e-7 0.89426831e-7 0.30912257e-6 0.75094831e-7 0.12017218e-6 11 2 轨道动力学计算的基本数学模型 卫星在轨道上运行要受到各种力因素的影响,产生的摄动是多方面的。国内外一些学者 对卫星轨道的受摄问题作了详细的研究与分析,尤以澳大利亚的 C.Rizos、A.Stolz 和美国 的 H.F.Fliegel 等人为代表。统筹考虑精度的需要和时间耗费,通过大量试算,本课题考虑 了卫星所受的以下作用力来进行轨道计算(注:时间系统采用 TDT 时间系统、坐标系统采用 J2000.0 惯性坐标系),主要包括:地球质心引力 F0、除质心外的地球引力 FE、太阳和月球 引力 FN、太阳辐射压力 FA、大气阻力(低卫星轨卫星)、Y 轴偏差 FY、地球潮汐附加力 FT。 F  F0  FE  FN  FA  FT  FY (2.1) 其中地球质心引力 F0 是最主要的,其次是地球的非质心引力 FE,称为地球非球形摄动力。 如果将地球质心引力视为 1,地球非球形摄动力可达 10-3 量级,而其它摄动力则大多在 10-6 以下。 2.1 二体问题 在惯性系中,卫星运动方程被描述为 r  r0  r(t,  r, r )   GM  r 3 e  r  r(t,  r, r ) 其中:  r, r 和 r 分别表示时刻卫星在惯性系中的位置、速度和加速度矢量; (2.2) G 和 M e 为分别为引力常数和地球总质量。 (3.2)式的第一项为地球质心引力项,称为二体运动,是力模型的主项;第二项为摄 动力引起的总摄动项,是 t ,  r, r 的函数矢量。 2.2 地球非球形引力摄动 在地固坐标系中,地球引力位函数作为拉普拉斯方程的解,其非球形部分为: Nn   U  (C nmU m n  S nmVnm ) n2 m0 式中: (2.3) 12 U m n  GM eaen Pnm (sin  ) cos m r n1 Vnm  GM eaen Pnm (sin  ) sin r n1 m 其中:λ和φ表示单位质点在地固坐标系中的地心经、纬度; ae 表示地球平均半径; Pnm (sin  ) 为规格化的缔合勒让德多项式; Cnm 和 S nm 为规格化的地球引力位系数; n 和 m 为多项式的阶和次,N 为取的最高阶数。 非球形引力摄动的求解,按如下步骤进行: a. 首先求解 U m n 和 Vnm 。用下列公式逐阶次推算得到,即: VUnnmm11     ae r I { m n1 2n  1 sin  VUnnmm     J m n1 VUnnmm11  }   VUnnnn1111     ae r K nn11{VUnnnn    cos  cos   UVnnnn    cos sin } 式中:   arctg( y ) ,   arctg( z ) x x2  y2 Im n 1  2n  3 (n  m  1)(n  m  1) Jm n1  (n  m)(n  m)  ae 2n 1 r  3, n0 K n1 n1      2n  3 , 2n  2 n0 递推方法如下所示: 递推初值为: U 0 0  GM e r , U 0 1  3 ae r sin U 0 0 , V00  0 , V10  0 递推方法采用先求对角线,再按行递推(m 不变,n 增加)进行。 (2.4) (2.5) (2.6) (2.7) 13 U 0 0  U 0 1  U 0 2  U 0 3  U 0 4  U 0 N  0 U11 U 1 2 U 1 3  U 1 4 U 1 N  0 U 2 2  U 2 3 U 2 4  U 2 N  0 U 3 3 U 3 4 U 3 N  0 U 4 4 U 4 N  U N N 如果采用非规格化的系数 U 和 V 进行计算,递推公式为: VUnmnm11     r(n ae m {(2n  1)  1) sin  VUnmnm     ae r (n  m)VUnmnm11  }   VUnnnn1111     (2n  1) ae r {VUnnnn    cos cos   UVnnnn    cos sin } 递推初值为: U 0 0  GM e r , U 0 1  ae r sin U 0 0 , V00  0 , V10  0 递推方法同上。此方法完成递推后,还需将 U 和 V 规格化,其公式为: U n  2n  1U n U m n  2(2n  1) (n (n   m)!U m)! m n b. 计算 U m n 和 Vnm 的偏导数 U m n 和 Vnm i. 首先计算 Anm  ka 2 (n  m  1)(n  m  2)(2n  1) /(2n  3) Bnm  kb 2 (n  m  1)(n  m  2)(2n  1) /(2n  3) Enm  (n  m  1)(n  m  1)(2n  1) /(2n  3) 其中乘数因子为 (2.8) (2.9) (2.10) 14 1 , m  1  K a   2 / 2, m  0   2, m 1 ii. 计算 U m n 和 Vnm 1 , m  0 Kb    2 / 2, m0   U m n    1 ae AnmAUnmVnmn1m111 BnmU m1 n1  BnmVnm11     EnmU m n1     Vnm    1 ae    AnmVnm11 AnmU m1 n1   BnmVnm11 BnmU m1 n1     EnmVnm1   而: U m n Vnm  ( 1) mU m n  (1) V m1 m n c. 求非球形引力摄动引力位 Nn   U  (C nm U m n  SnmVnm ) n2 m0 则摄动加速度为: rE  [W ]T U 式中:矩阵[W]为地固坐标系转换至惯性坐标系的旋转矩阵。 2.3 日、月摄动 (2.11) (2.12) (2.13) (2.14) 考虑卫星的 N 体影响,只需顾及太阳和月球的引力作用已满足精度要求。N 体摄动模型 的建立是基于牛顿第二运动定律和万有引力定律。由图 2 所示,分析卫星的受力可得卫星 m Mj Me 图 2 三体的几何关系 的日、月摄动加速度矢量为: 15     rN   jS,L GM j ( j 3j  rj ) r 3 j  (2.15) 其中:   j  r  rj , j 为摄动体至卫星的中心距,即矢量  j的模; rj , r 分别为摄动体 M j 和卫星在惯性系中的位置矢量, rj , r 的计算参见本节附录。 这里 S 和 L 分别代表太阳和月球。 太阳和月球对卫星的摄动影响主要呈长周期变化,且与卫星轨道对太阳和月球的定向有 关。 2.4 太阳直接辐射压摄动 照射在卫星上的太阳光,一部分被其吸收转化为热能,另一部分被反射向太空。因此, 卫星会受到照射时的辐射力和反射时的反射力的作用,这里统称为直接辐射压摄动。直接辐 射压摄动与光压强度、卫星表面的反射系数和光照面积有关。 由光压和牛顿第二运动定律建立的直接辐射压对卫星产生的运动加速度矢量为: rA   k  a 2 sun Ta2s a e3 2s  Gx  sin(  0.015sin 2 )  Gz cos(  0.025sin 2 ) (2.16) 其中: k  1.013107  2.4109 cos 2 1.3109 cos 4  4.0 1010 sin 6 ; Gx,Garzc对co应s卫星(rr表,面ss )反射;系数项,是为补偿模型不足而引入的拟合参数; asun  1.4959787 1011 m 为天文单位长度(地球轨道长半径); Tas  8.068111241279087 10 2 为人卫时间单位;   s 为卫星的日心距,即矢量 s  rs  r 的模,r , rs 分别为卫星和太阳在惯性系中的位 置矢量;  为蚀因子,具体定义为:  1    0 卫星在地影和月影之外 卫星在地球或月球的本影之内  1  As' 卫星在地球或月球的伪本影、半影之内  As (2.17) 其中: As' , As 分别为太阳的被蚀视面积和视面积。要确定蚀因子,需计算下列诸量:    As Am   2 s   2 m   Ae   2 e (2.18) 16   arctg(z / x2  y 2 ) cos2   (x2  y 2 ) / r 2   s   sin 1 ( a ' s ) s  m   sin 1 ( a ' m ) m    e   sin 1 ( ae' r ) aaem's' ses(aacscoeos1s121(0(f0rrsmm0 0)ss 1 ) s) s fe cos 2  式中: Am , Ae 为月球和地球的视面积; f s , fe 分别代表太阳和地球的扁率, fe  3.352813178103 ; a ' e 为考虑地球大气衰减以及地球扁率效应的有效地球赤道半径; am =1738000m 为月球半径; a ' s  6.96 108 m 为考虑太阳扁率的太阳有效半径;  为地心纬度;  es 是地球—卫星—太阳张角;  ms是月球—卫星—太阳张角。 地影和月影判断过程如图 3 所示。 当卫星在地球半影中时: Ae'   2 s cos 1 (e s )   2 e cos 1 ( es  e  e )   es 其中: e   2 es   2 s   2 e 2 es 当卫星在地球伪本影中时: Ae'   2 e  2 s   2 e 当卫星在月球半影中时: Am'   2 s cos 1  ( m s )   2 m cos 1  ( ms  m m )   ms 其中: m   2 ms   2 s   2 m 2 ms 当卫星在月球伪本影中时: Am'   2 m  2 s   2 m (2.19) (2.20) (2.21) 17 则蚀因子为:  1  As'  1  max( Ae' , Am' ) As  2 s 无地影,无月影 Y F=1 Y 结束 无地影 Y Ae=0 Y 无月影 Y Am=0 es<0 且ms<0 N es<0 N es(e+s) N es|e-s| N 在地球半影内 计算 Ae ms<0 N ms(m+s) N ms|m-s| N 在地球半影内 计算 Am 计算 F Y Y se N 在地球伪本影内 计算 Ae Y Y sm N 在地球伪本影内 计算 Am 在地球本影内 F=0 结束 在地球本影内 F=0 结束 结束 图 3 地影和月影判断流程 太阳直接辐射摄动对卫星轨道的影响是十分复杂的,它与卫星表面的反射特性、卫星轨 道相对太阳的定向以及太阳活动等有关。卫星是由各种不同折射性质的原材料构成的不规则 形体,在其运行过程中,太阳对它照射的面积也在不断地改变(它的太阳能翼始终是面向太 阳的)。此外,由于太阳活动的变化,所谓太阳常数 Ps 也并非常数。因此,给出卫星受太 阳直接辐射压摄动的精确模型是很困难的。所以采用较简单的平面模型计算太阳辐射加速度 影响。 18 2.5 地球固体潮摄动 地球并非刚体,它受日月引力作用会产生弹性形变,称为潮汐现象。这种形变使得地球 内部物质发生小的变化,随之导致引力位函数产生小的形变位差——潮汐位。卫星运动的地 球固体潮摄动就是潮汐位效应的结果。 已知日(或月)对地面点的引力位球谐展开式为:  v j  GM j aen r n j 1 Pn (cos j) 式中: M j 为日(或月)的质量, rj , j 分别为引力体至地心距和与地面点的地心角, Pn (cos j ) 为 n 阶勒让德多项式。 从上式中排除不产生形变位差的 0 或 1 阶项,且只取 n=2 阶项,可得日月潮汐形变对卫 星产生的摄动位为:  U et  GM j r3 jS,M j a 5 e r3 k2 P2 (cos j ) 其中: k 2 为二阶 Love 数(取 k2  0.3 )。 顾及关系式 P2 (cos j)  1 2 [3(rˆjT rˆ) 2 1] ,最后得到卫星的固体潮摄动加速度矢量为:  rT   Uet  r T    k2 2 GM j  jS ,L  r 3 j ae5 r4 {[3  15(rjT  r ) 2  ]r  6(rjT  r )rj } 式中: rˆ, rˆj ( j  S, L) 分别为 r , rj 的单位矢量。 (2.22) 2.6 大气阻力摄动 高层大气对卫星的运行将产生阻力,这种阻力对低轨道卫星是主要摄动力之一,但大气 密度变化机制非常复杂,不但随高度变化,也与太阳活动、时间、季节、纬度和地磁活动的 变化有关。本文采用了静止球型大气密度模型(HP 模型)。若只考虑大气分子对卫星表面的 法向作用力而忽略其切向作用力,大气阻力使卫星产生的摄动加速度为: rD  rDB  rDP   1 2 C D    A m VRVR  1 2 C DP    AP m   cos VR 2   s 其中: rrDDBP 为卫星星体部分的大气阻力摄动加速度; 为卫星太阳帆板的大气阻力摄动加速度; (32.23) C D 为大气阻力系数; VR是是大卫气星密相度对,大由气H的P速模度型矢计量算V得R 到 ;r      r ; 19 A 是卫星参考面积与卫星质量之比; m C DP 为太阳帆板的大气阻力系数; AP 为太阳帆板的面积;  为太阳帆板的法向与卫星相对于大气的速度方向的夹角;   s 为  s 向量的单位向量。 2.7 Y 轴偏差加速度摄动 Y 轴偏差加速度主要是 GPS 卫星的结构失调和卫星体的热辐射而产生的一个附加加速度 在星固坐标系 Y 轴方向上的分量。在设计上,为使卫星的太阳翼以最大面积朝向太阳,两翼 的支轴应保持在一条直线上,并要求垂直于卫星和太阳方向的连线,用于控制两翼俯仰的太 阳传感器也应完全垂直于卫星翼支轴,卫星的偏航高度控制应保持正确,但事实上并不完全 这样;另一方面,由卫星本身产生的超高温要从 Y 轴方向的百叶孔排出,这样使处在不稳定 状态中的卫星体也会产生结构失调现象。由此导致了 Y 轴偏差加速度的存在。H.F.Fliegel 等人把结构失调的影响表示为[14]: rY   Gy ya (2.24) 其中:yGa 为y 为星考固虑坐模标型系剩的差Y而轴引在入J的20待00估.0参惯数性;坐标系中的方向余弦,由下式决定:   ya s  rs  s s    r  r  r 2.8 巡航姿态控制动力摄动 有些卫星在巡航过程中需要保持三轴稳定姿态,需通过姿态控制实现,有的卫星其姿态 控制的动力来源于高压气瓶的喷气。这样,在姿态控制的同时也影响了卫星质心的运动。该   摄动加速度矢量可以表示为:rA  W T  Co   C1T    C cosu  S sin u (2.25) 其中: u    f ; Co 为卫星在 RTN 坐标系中姿控动力摄动加速度的常数分量; C1 为卫 星在 RTN 坐标系中姿控动力摄动加速度的时间变化率; C 和 S 为周期项的系数; 对全局参数,为由历元时刻起算的相对时间;对弧段相关参数,为观测时刻 t 在相应弧 段内的相对时间。 20 2.9 其它摄动影响 在轨道运行的卫星除受到上述摄动作用外,还受其它一些摄动的影响,如反照辐射摄动、 地球自转形变摄动、广义相对论摄动、海潮摄动、大气潮摄动等。这些摄动影响对卫星轨道 摄动非常小,但其计算却要耗费大量的时间片。考虑到 2.1 节所述的摄动已能满足课题的精 度需求和时间消耗的限制,因此,本课题中忽略了其它摄动的影响。 附录:日月位置计算 a. 太阳位置矢量 rs 计算  rs  [D]T rs cos s  R x ( 0 ) rs sin  s    0  其中:[D]T 为计算历元时刻的平春分点到 J2000.0 平春分点的转换矩阵,即岁差矩阵。 rs  asun [1  1 e2 2  e cos M  1 e2 2 cos 2M ] s  L  2e sin M  5 e2 4 sin 2M e  0.01675104  0.418 104 TB  0.126 106 TB2 计算当时平春分点的几何平黄经为: L  27941'48''.04 129602768''.13TB 1''.089TB2 平近点角为: M  35828'33''.04 129596579''.10TB  0''.54TB2 平黄赤交角为:  0  2327'08''.26  46''.845TB  0''.0059TB2 以上式中 TB  JD(TDB)  2415020 36525.0 b. 月球位置矢量 rm 计算  rm  [D]T R x ( 0 ) Rz (~ ) Rx ( J~)rrmm cos(~ sin(~   ~ ) ~ )     0    其中: ~    0.47 103 sin M '  n  j sin a j j 1 21  J~  J  n J j cos a j j 1  ~  L'  C1  n K j sin a j j 1  1 r  1 60.2664 [1  C 2  n Pj j 1 cos a j ] rm  r  ae J  5.1454 L' 27026'02".99 1732564379".31TB  4".08TB2   25910'59''.79  6962911''.23TB  7''.48TB2 M '  29606'16''.59 1717915856''.79TB  33''.0912TB2 C1  0.10976sin M '  0.00373sin 2M '  0.17 103 sin 3M ' C2  0.0545cos M '  0.00297cos 2M '  0.18103 cos3M ' 系数 K j 、 Pj 、  j 和 J j 的值以及幅角 a j 的计算列入表 5 中,表内 D 为日月平角距, F  L' ,其计算式为: D  35044'14''.95 160291611''.18TB  5''.17TB2 F  1115'03''.20 1739527290''.54TB 11''.56TB2 22 表 5 计算月球位置的系数表 j Kj Pj j Jj aj 1 0.022230 0.010020 0.000820 0.000000 2D-M' 2 0.011610 0.008250 -0.002300 -0.000190 2D 3 0.003240 0.000120 0.002830 0.000000 M+180° 4 0.000000 0.000000 -0.002130 -0.000190 2F+180° 5 0.001030 0.000090 0.000450 0.000000 2M'-2D+180° 6 0.001000 0.000420 0.000000 0.000000 2D-M'-M 7 0.000930 0.000900 -0.000080 0.000000 2D+M' 8 0.000800 0.000560 -0.000120 0.000000 2D-M 9 0.000720 0.000340 0.000000 0.000000 M'-M 10 0.000610 0.000290 0.000000 0.000000 D+PI 11 0.000530 0.000280 0.000000 0.000000 M'+M+180° 12 0.000270 0.000000 -0.027240 -0.002410 2F-2D+180° 13 0.000000 0.000000 -0.000080 0.000000 2F+M'+180° 14 0.000410 0.000200 0.000700 0.000000 2F-M'+180° 15 0.000180 0.000180 0.000000 0.000000 4D-M' 16 0.000150 0.000120 0.000000 0.000000 4D-2M' 17 0.000140 0.000070 0.000000 0.000000 2D+M-M'+180° 18 0.000120 0.000090 0.000000 0.000000 2D+M+180° 19 0.000090 0.000000 0.000000 0.000000 M'-D 20 0.000090 0.000000 0.000000 0.000000 D+M 21 0.000070 0.000070 0.000000 0.000000 2D-M+M' 22 0.000070 0.000090 0.000000 0.000000 2D+2M' 23 0.000070 0.000080 0.000000 0.000000 4D 24 0.000000 0.000000 0.000190 0.000180 2F-2M' 25 0.000000 0.000000 0.001100 0.000100 2F-2D+M 26 0.000000 0.000000 0.000470 0.000000 2D-2F+M 27 0.000000 0.000000 0.000230 0.000000 2F-2D-M' 28 0.000000 0.000000 0.000023 0.000000 2D-2F-M' 23 3 轨道计算方法 卫星运动方程的解有分析法、数值法和半分析半数值方法等。分析法是将力模型展开取 有限项,给出任意时刻解的表达式,通常称之为一般摄动法。其优点是便于定性地给出轨道 的特征,如长期项、长周期项影响等。但对摄动力模型处理限制较大,使精度受到影响。数 值法即常微分方程的初值解问题,在轨道计算中称之为特别摄动法。其优点是能完整地顾及 所选择的力模型、处理简单、计算精度高,是高精度卫星轨道计算最实际有效的方法。 经研究,Runge_Kutta 型和 Cowell 型数值积分方法都可用于产生高精度的卫星轨道, 尤其是后者通过预报——校准公式进行迭代计算,收敛快(一般只需迭代 2~3 次),累积误 差影响要比前者小的多。所以这一方法是轨道计算最常用的方法之一。不足的是它为多步型 积分器,必须用其它方法为它准备起步值。 本课题提供的积分模式有 Runge_Kutta 4、Runge_Kutta 8 和高精度的 Adams_Cowell 方法,高精度计算时采用 Runge_Kutta 单步积分起步,用 Adams_Cowell 多步积分求解卫星 运动方程,确定卫星在 J2000.0 惯性坐标系中的位置和速度矢量。 3.1 Runge_Kutta 积分法 Runge_Kutta 法是一种单步积分方法,公 式形式简单,应用方便,将它用于解算卫星 开始 运动方程和变分方程也具有很好的稳定性。 但该方法随着积分式阶数的增高,计算右函 数的次数也会随之增多,无疑会导致累积误 赋初值:t0,y0 求取右函数:F0 差增加,同时也会降低计算效率。因此,在 tk  t0  ak h 大规模的轨道计算中,通常只将其作为多步 k 1 型积分器的起步方法。流程图如图 4 所示。  yk  y0  ak h bkj Fj j0 k++ 设有初值问题: y(t)  f (t, y)   y(t n )  yn (3.1) 则求解该问题的 Runge_Kutta 积分公式为: k  yn1  yn  h Ci fi i0 (3.2) 求取右函数:Fk(tk,yk) k>m ? k  y  y0  h ci F i i0 求取右函数:f(t+h,y) 图 4 Runge_Kutta 积分框图 24  f0  f (tn , yn ) f1  f (tn  a1h, yn  a1b10 hf0 )  k    f k  f (tn  ak h, yn  ak h bki i0 fi ) 式中: h  tn1  tn ,为积分步长; k 为积分式所取的阶数; Ci , ai , bij 均为已知的系数,具体参见表 6 所示。 由运动方程和变分方程形成的矩阵常微分方程的初值问题为: y(t)  F(t, y) 上式中:左端 y(t)   rr    ; 右函数 F   r (t ) r(t,  r, r )   (t) (t,  r, r)  ; 初值 y(t0 )   rr00   0 0  3.2 Adams_Cowell 积分 Adams_Cowell 积分适用显含 1 阶导的二阶微分方程,初始值需要位置项和速度项。积 分的预报公式为:    y n1  h[I sn  k  j yn j ]  j0  k    y n1  h 2[II sn   j yn j ] j0 (3.3) 其中: j ,  j 为 Adams_Cowell 积分公式的预报系数,按下列公式计算,结果参见表 4.2。    j   (1) j k m j  m j  m 1 (1)    j  (1) j k m j  m j  m  2 (1) 校正公式为: 25    y n1  h[I sn  k  ' j yn j 1 ]  j0  k    y n 1  h2[II sn   ' j yn j 1 ] j0 (3.4) 式中:  ' j ,  ' j 为 Adams_Cowell 积分公式的校正系数,按下列公式计算,结果参见表 4.2。    ' j   1   " j , j    " j , j  0 0   " j   (1) j k m j  m j  m 1 (1)     ' j   (1) j k m j  m j  m  2 (1) (3.3)、(3.4)式中的 (1) 和 (1) 按以下公式递推计算:   0 (0)  1     i(0)   i 1 j0 i  1 j  1  j (0)   i(0)  i  j (0) i j (0)  j0  i   i(1)   j (0)  j0  i    i(1)   j(0) j0 I sn 和 II sn 分别为 yn 的一次和分和二次和分,按下列公式递推计算:    I sn II sn  I sn1  II sn  1 yn I sn  I sn 和 II sn 的初值定义为: I sn1    II sn1  y n h yn h2 k   ' j j0 k  j0 yn j ' j yn j 用 Adams_Cowell 积分解卫星运动方程积分流程图如图 5 所示: 26 求积分常数  输入积分初值r: r 求 I sn1 , II sn1 求 I sn , II sn 计算右函数r 预报 y 0 n1 ,i=0 计算右函数r i++ 校正 y i1 n 1 是否收敛? N yn  yYn1 N 积分完成? Y 结束 图 5 Adams_Cowell 积分流程图 3.3 轨道计算 预测卫星轨道实际上就是由一组正确的初始值数值积分卫星运动方程,外推出未来的卫 星轨道。显然,在卫星运动力模型较完备的情况下,初始值的精确与否将直接影响外推轨道 的精度。且在推算过程中,由于受初始值误差传播影响,推算区间越长,积分累积的误差就 越大。轨道计算步骤为: 27 GMST 计算 岁差计算 章动计算 极移计算 求坐标转换矩阵 W 求日月位置 求 r0 求 rE … 求 rY Runge-kuta 法计算积分初值 Adams_Cowell 积分 轨道预推 图 6 卫星轨道计算流程图 3.4 星历的快速插值 通过分布式计算提供的星历数据间隔太大,而在实际定位计算中,有时需要间隔为 1 秒甚至步长更小的精确星历坐标。因此,高精度、快速的精密星历插值就成为一项重要工作。 以精度和计算耗费小为出发点,本课题在研究的过程中,分析和比较了以下几种多项式插值 方法:拉格朗日多项式插值、牛顿多项式插值、三角函数插值、样条函数插值、切比雪夫多 项式插值等。 拉格朗日多项式插值的多项式构造和插值时间耗费非常大。其时间耗费可由下式表示: (2n  2)A  (n  2)M  (n 1)D 其中: n 表示已知点数; A 表示加法运算、M 和 D 分别表示乘、除运算。 构造插值多项式还需额外的 (n 1)A  nM 运算。因此,总运算耗费为: (3n  3)A  (2n  2)M  (n 1)D 比拉格朗日多项式插值更有效的方法是利用数据均差的牛顿多项式插值[1][2]。假设有 n 1个控制点 t0 ,t1,,tn ,定义在 ti 处的零阶均差为 f [ti ]  f (ti ) ,在 ti , t j 处的 1 阶均差       定义为 a0  f ti ,t j  f ti  f t j ti  t j ,k 阶均差定义为: 28 ak  f t0 ,t1,,tk   f   t0 ,t1,,tk1  f t1,t2 ,,tk  t0  tk k 0 则牛顿插值多项式可表示为: pn (t)  a0  (t  t0 ){a1  (t  t1){a2   (t  tn2 ){an1  (t  tn1)an}}}(4.12) n1 上述表达式在每级中都由 2 个加法和 1 个乘法运算组成,其时间耗费为:n(2A  M ) 。 与经典的拉格朗日多项式插值相比,牛顿多项式插值经济得多。给定 n 1个不同的点   t0 ,t1,,tn 和相应的均差系数 a0 , a1,, an ,对任意时刻 t  t0 ,tn 插值多项式的结果由下 列迭代产生(Horner 算法)。 bn  an for k  n 1 to 0 step by 1 bk  ak  (t  tk )bk1 end 上述算法的精度和拉格朗日的精度相同,但时间耗费得到了明显的改善,而且各阶均差 可以预先运算好储存起来,需要时直接调用即可。 插值点等间隔条件下的插值是牛顿多项式插值的一种特例,可以导出更简洁的形式,而 不降低精度。IGS 提供的精密星历数据是等间隔的,步长为 15 分钟。在实际应用中,我们 采用牛顿前向差分形式。 牛顿前向均差的定义为: f (ti )  f (ti1)  f (ti )  f (ti  h)  f (ti ),tk  t0  kh, k  0 同理可推出各高阶前向均差为:   k f (ti )  (k1 f (ti ))  k1 f (ti1 )  k1 f (ti )  k j0 (1) k  j  k j  f (ti j ) (4.13) 假设插值多项式的阶数为 n ,则插值多项式为: pn (t)  f (t0 )  (t  t0 h ) f (t0 )    (t  t0 )(t  t1 )(t n!h n  t n 1 ) n f (t0 )  令  (t  t0 ) h ,则上式可表示为: pn ( )  j n 0 j  j f (t0 ) , (4.14) 用(4.13)式构造插值多项式,再用(4.14)式用作插值,该方法在时间耗费上比 Horner 算法有一定的改善。 与经典的拉格朗日插值方法相比,三次样条多项式插值速度快,但精度低,三角函数多 项式插值必须用 DFT 求解系数,时间耗费大,且这两种方法都将星历数据近似为周期函数, 插值精度低,不能满足高精度插值的需要。切比雪夫多项式插值方法在端点处抑制了误差的 扩大,但计算量较大。牛顿多项式插值方法的精度高,速度快,充分利用了星历数据的特点, 且克服了拉格朗日方法的不足。在精密星历插值上,牛顿多项式插值方法比拉格朗日多项式 插值方法更为有效。由拉格朗日插值多项式和牛顿插值多项式的余项公式可推出截断误差在 29 端点处将迅速增长。为克服这种不足、提高插值的精度,在实现时提出并采用了“滑动式牛 顿多项式插值”方法[19],即采用滑动方式进行插值:每次取 13 个点,生成 12 阶多项式,始 终让被插值点在已知点的第 6 个点到第 7 个点之间。第 1 个点到第 13 个点为第一个插值区 间,仅用来插值第 6 个点到第 7 个点之间的时间段;第 2 个点到第 14 个点为第二个插值区 间,仅用来插值第 7 个点到第 8 个点之间的时间段;依次类推下去。这种滑动式牛顿多项式 插值精度非常高,可达到 cm 级精度,且易于实现,时间耗费小,确保了精密轨道解算时的 精度和实时要求。 表 6 8 阶 Runge_Kutta 积分公式系数表 ci 1 2 1 0 2 4 27 1 27 4 3 2 9 1 18 1 3 4 1 3 1 12 1 0 5 1 1 1 0 28 6 2 3 1 54 13 0 7 1 6 1 4320 389 0 8 1 1  -231 0 20 9 5 6 1 288 -127 0 10 1 1 820 1481 0 ai 1 840 41 0 bij 3 4 5 6 7 8 9 10 3 0 3 -27 42 8 -54 966 -824 243 81 -1164 656 -122 800 18 -678 456 -9 576 4 -81 7104 -3376 72 -5040 -60 720 0 27 272 27 216 0 216 41 表 7 12 阶 Adams_Cowell 积分公式系数表 i i 0 4.259634 1 -22.691026 i 0.877861 -4.800775  ' i 0.264351 0.823067  ' i 0.054794 0.165534 30 2 80.019378 3 -196.294318 4 349.412940 5 -462.480193 6 460.087276 7 -343.735149 8 190.394308 9 -75.976728 10 20.680772 11 -3.441245 12 0.264351 16.861107 -41.187179 73.086989 -96.516821 95.852539 -71.516859 39.571065 -15.777189 4.291462 -0.713663 0.054794 -2.071621 4.414893 -7.283104 9.192754 -8.853279 6.460362 -3.514964 1.383094 -0.372243 0.061367 -0.004677 -0.526813 1.189914 -2.009198 2.566623 -2.489666 1.825385 -0.996493 0.393084 -0.105997 0.017501 -0.001336 31 4 轨道根数与位置矢量、速度矢量的关系 4.1 由位置矢量和速度矢量计算轨道根数 已知某时刻 t 的位置矢量  r 和速度矢量 r ,计算轨道根数的步骤如下: 由  h   r  r  hX   hY     hZ  其中 h  h  hX2  hY2  hZ2 ,分别计算 i 和  而 计算  计算 a 类似可得 所以:   cosi   tan   hZ h hX   hY    e 1  r   h  r r eX     eY   eZ   e e  e 2 X  eY2  eZ2 tan   eY sin   eZ eX cos sin i   a   h2 1 e2 tan u  Y sin   Z X cos sin i f u 计算 ,由 tan E  1  e tan f 2 1e 2 得 (5.1) (5.2) (5.3) (5.4) (5.5) (5.6) (5.7) 32 nt    E  esin E (5.8) 4.2 由轨道根数计算位置矢量和速度矢量 若已知任何时刻 t 时天体的 r 和 f ,在轨道坐标系 O  xyz中有  x  r cos f   y   r sin f   z   0  根据惯性坐标系与个坐标系的转换关系可得位置矢量的表达式 (5.9)  r    X Y    R3  R1  iR3    r r cos sin f f     Z     0   将旋转矩阵展开可得     r  r cos f  P  r sin f  Q 其中 P 和 Q 分别表示 x 轴和 y 的单位矢量,即 (5.10)  P  R3  R1  iR3    1 0     cos  cos sin  cos  sin  sin  cos  sin   cos i cos i   0  sin  sin i  (5.11)    Q  R3  R1 iR3 1   cos  sin   sin  cos cosi    90o  0     sin  sin    cos  cos cosi  (5.12) 0  cos sin i  速度矢量表达式为: r   a2n sin E   P  a2n  1  e2 cos E  Q r r (5.13) 33

Top_arrow
回到顶部
EEWORLD下载中心所有资源均来自网友分享,如有侵权,请发送举报邮件到客服邮箱bbs_service@eeworld.com.cn 或通过站内短信息或QQ:273568022联系管理员 高进,我们会尽快处理。