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

目录TOC\o"1-5"\h\z\o"CurrentDocument"1星历计算的时间和坐标系统 2\o"CurrentDocument"1.1有关的时间系统与坐标系统 2\o"CurrentDocument"1.1.1时间系统及其换算 2\o"CurrentDocument"1.1.2坐标系统及其换算 3\o"CurrentDocument"1.2计算单位和有关常数 6\o"CurrentDocument"2轨道动力学计算的基本数学模型 12\o"CurrentDocument"2.1二体问题 12\o"CurrentDocument"2.2地球非球形引力摄动 13\o"CurrentDocument"2.3日、月摄动 15\o"CurrentDocument"2.4太阳直接辐射压摄动 16\o"CurrentDocument"2.5地球固体潮摄动 18\o"CurrentDocument"2.6大气阻力摄动 19\o"CurrentDocument"2.7Y轴偏差加速度摄动 19\o"CurrentDocument"2.8巡航姿态控制动力摄动 20\o"CurrentDocument"2.9其它摄动影响 20\o"CurrentDocument"附录：日月位置计算 21\o"CurrentDocument"3轨道计算方法 25\o"CurrentDocument"RungeKutta积分法 26\o"CurrentDocument"AdamsCowell积分 263.3轨道计算 283.4星历的快速插值 28\o"CurrentDocument"4轨道根数与位置矢量、速度矢量的关系 32\o"CurrentDocument"4.1由位置矢量和速度矢量计算轨道根数 32\o"CurrentDocument"4.2由轨道根数计算位置矢量和速度矢量 331星历计算的时间和坐标系统1.1有关的时间系统与坐标系统轨道计算过程重要涉及到不同的时间系统和坐标系统，下面将空间战场环境系统中所涉及到的时间系统和坐标系统进行定义，并说明各系统之间的相互关系。一般情况下，仿真系统采用的是TDT时间系统和J2000地心惯性坐标系。1.1.1时间系统及其换算在轨道计算中，时间是独立变量。但是，在计算不同的物理量时，却使用不同的时间系统。例如：在计算恒星时用世界时UT1；定位解算时采用GPS时GPST；岁差和章动量的计算采用TDB时等。所以必须清楚各时间系统的定义和各时间系统之间的转换，下面给出各种时间系统的定义及它们之间的转换公式。格林尼治恒星时格林尼治恒星时为春分点对格林尼治平天文子午面的时角。由于岁差、章动原因，它由格林尼治真恒星时（GAST）和平恒星时（GMST）之分。两者的关系是：其中：Avcos8为赤经章动GMST=67310，.54841+（8640184，.812866+876600h）T+0s.093104T2-0..62x10-5T3T为自J2000.0（JD2451545.0）起算至观测UT1时刻的儒略世纪数，即世界时UT1UT1是以平北极（国际习惯用原点）为统一标准的观测世界时，是反映地球实际自转的时间，恒星时计算与此有关。国际原子时TAITAI时以铯原子C133基态两 能级间跃迁辐射的9192631770周所经历的时间作为1秒长的均 2匀时间，起点在1958年1月1日0hUT。国际协调时UTCUTC是经跳秒修改后的国际原子时，它与世界时UT1的差0s.9，观测纪录时间是以此为准的。质心动力学时tdb(BarycentricDynamicalTime)TDB为相对于太阳质心的运动方程给出的历表、引数等所用的时间尺度，岁差及章动量的计算是以此为依据的。地球动力学时tdt(TerrestrialDynamicalTime)T^D^T为视地心历表所用的时间尺度，它具有均匀连续的特性，卫星运动方程就是以此为独立的时间变量。GPS时间GPSTGPST是由系统定义和应用的一种时间尺度，于1980年1月6日0hGPST与UTC相等，在此以后由系统主控站密切跟踪UTC以保持高度统一。但GPST不作跳秒修正，因此它与UTC具有整秒的差异(1997年1月至6月相差为iis)。在计算GPS卫星轨道的初值时将涉及到GPST，GPS精密星历的参考时间为GPST。以上各时间尺度的相互关系如下：其中：AUT1可从地球自转参数文件中获得；AAT=19s+AGPST；ATD=0s.001658sin(v+0.0167sinv),v=6.240040768+628.3019501T(rad)。上式中的T为自j2000.0年起算至观测TDB时刻的儒略世纪数，即：不同时间系统间的关系如下图所示：图1几种时间系统之间的关系1.1.2坐标系统及其换算卫星轨道计算和实际定位解算分别是在J2000.0惯性坐标系与WGS-84地固系算中还涉及到星固坐标系。下面的定义及相互转换关系。中进行的，此外，卫星加速度计给出与本课题有关的主要坐标系算中还涉及到星固坐标系。下面的定义及相互转换关系。J2000.0地心惯性系原点：地球质心Z轴：向北指向J2000.0年平赤道面（基面）的极点X轴：指向J2000.0平春分点Y轴：符合右手系法则位置矢量：r星固坐标系原点：卫星质心Z轴：指向卫星的天线方向，即指向地心X轴：在轴与太阳构成的平面内，完成右手系法则Y轴：沿卫星太阳能翼的支轴位置矢量：ra星固坐标系坐标轴伐，a,Z）在惯性系中的方向余弦分别为：（r,r分别为太阳aaa s和卫星在J2000.0地心惯性系中的位置矢量）WGS-84坐标系WGS-84为1984年世界大地坐标系（WGS为WorldGeodeticSystem的简称），WGS-84的坐标定义及其采用的椭球参数为：原点：地球质心Z轴：指向BIH1984.0定义的协议地球极（CTP）方向X轴：指向BIH1984.0的零子午面和CTP赤道的交点Y轴：与X、Z轴成右手系地球椭球长半径：a=6378137m地球引力常数（含大气层）:GM=3986005x108m3/s2正常化二阶带球谐系数：F20=-484.16685x10-6地球自转角速度：①=7292115x10-11rad/s2地球椭球扁率：f= 1/298.257223563地固坐标系与惯性坐标系的 转换模型惯性坐标系转换到＞地固坐标系模型地固坐标系转换到＞惯性坐标系模型式中：［A］为极移矩阵；［B］ 为自转矩阵；［C］ 为章动矩阵；［D］ 为岁差矩阵。上述各矩阵的意义及具体定义如下：极移：由于地球不是刚体及其它一些地球物理因素的影响，地球自转轴相对于地球的位置随时间而变化从而引起观察者的天顶在天球上的位置发生变化，称为极移，矩阵为［A］：式中：x,y为地极坐标，可从地球自转参数文件中给出的极移值插值得到。自转：即地球公转的同时也在绕自转轴旋转。矩阵［B］:式中：9G为格林尼治恒星时章动：是指外力作用下，地球自转轴在空间运动的短周期摆动部分，即同一瞬间真天极相对平天极的运动，月球对地球引力的变化是形成章动现象的主要外力作用，其次是太阳。矩阵［C］：式中：屁=如dcos(£nA) 交角章动J=1 k=1△中=如csin也nA) 黄经章动J=1 k=1其中：\dj,n狄都为常数，可自章动系数表1中查出。而月亮平近点角：太阳平近点角：月亮平升交角：日月平角矩：月亮轨道对黄道平均升交点 5的黄经：

其中：1r=360T0为自J2000.0起算至t的儒略世纪数岁差：地球在太阳、月球和行星的引力作用下，地球的自转轴在空间不断发生变化，其长期运动称为岁差，矩阵［D］：式中：T0的意义同上。旋转矩阵的求法分别为1.2计算单位和有关常数轨道计算采用人卫单位系统，具体定义为:长度单位:质量单位:地球赤道半径a长度单位:质量单位:e地球的总质量Me1时间单位:T=时间单位:T=as=8.068111241279087x102在以上人卫单位系统中，引力常数G=1。为完整起见，给出以下常数:地球赤道半径：6378137m地球扁率：f=1298.257地球总质量：M=5.974x1024kg地球自转角速度：o=7.292115x10-5rad/s地心引力常数：GM=3.986005x1014m3/s2日心引力常数：GM=1.32712438x1020m3/s2天文单位长度：a=1.4959787x1011m月球地球质量比：GM；GM=0.01230002光速：c=299792458m/s引 力 常 数： G=6.672x10-11m3/(kg-s2)6太 阳 吊 数： P=4.5605x10-6kg/(m•s2)地球引力位系数采用WGS-84EGM的规格化值。参见表2。表1黄经和倾角章动序列表引数黄经章动倾角章动LLFD0".00010”.0001T0".00010”.0001T100001-171996-174.2920258.920000220620.2-8950.53-20201460-2404202-00110005-20202-301061-0-10-3000702-2-21-20108202-0110009002-22-13187-1.65736-3.110010001426-3.454-0.111012-22-5171.2224-0.6120-2-22217-0.5-950.313002-211290.1-70014200-204801015002-20-22000160200017-0.1001701001-1509018022-22-160.17019020021-6030210-2-21-5030

22200-2140-2023012-2140-2024100-10-400025210-20100026002-2110002701220-10002801002100029-10011100030012-20-10003100202-227-0.2977-0.532100007120.1-703300201-386-0.420003410202-3010129-0.135100-20-1580-1036-102021230-5303700020630-203810001630.1-33039-10001-58-0.132040-10222-5902604110201-5102704200222-3801604320000290-1044102-22290-1204520202-3101304600200260-1047-10201210-10048-10021160-8049100-21-1307050-10221-1005051110-20-7000520120270-30530-202-70305410222-80305510020600056202-2260-305700021-60305800221-703059102-2160-3060000-21-5030611-00050006220201-503063010-20-400064102-0040006500010-40006611000-300067102003000681-202-301069-1-222-301070-20001-20107130202-3010720-222-3010731120220-1074-102-21-2010752000120-107610002-201077300002000780021220-1079-1000210-1080100-40-100081-2022210-1082-10242-201083200-40-100084112-2210-108510221-101086-20242-1010

87-104021000881-0-20100089202-2110-109020222-10009110021-100092004-22100093302-22100094102-20-10009501201100096-1-021100097002-01-100098002-12-10009901020-1000100102--20-10001010-201-1000102110-21-1000103102-20-100010420020100010500242-1000106010101000表2计算地球引力位加速度的引力位系数（GEM-T3）20-0.48416510e-30.021-0.17e-91.19e-9220.24390658e-5-0.14000946e-5300.95720109e-60.0310.20277142e-50.24921712e-6320.90447073e-6-0.61944767e-6330.72034249e-60.14138845e-5400.53952118e-60.0

41-0.53615108e-6-0.47343598e-6420.35021814e-60.66301523e-6430.99093372e-6-0.20092742e-644-0.18877065e-60.30942370e-6500.68343345e-70.051-0.58280231e-7-0.96083937e-7520.65271099e-6-0.32386369e-653-0.45233006e-6-0.21529578e-654-0.29558409e-60.49690346e-7550.17376349e-6-0.66890704e-660-0.14951352e-60.061-0.76894161e-70.26998384e-7620.48734482e-7-0.37401311e-6630.57203194e-70.93727543e-864-0.86826474e-7-0.47130637e-665-0.26733038e-6-0.53678021e-6660.96846267e-7-0.23713482e-6700.91300885e-70.0710.27486873e-60.97465920e-7720.32779498e-60.93246712e-7730.25122012e-6-0.21529269e-674-0.27556101e-6-0.12376718e-6750.13261698e-80.18620005e-776-0.35883140e-60.15173875e-6770.97027653e-90.24083642e-7800.48883181e-70.0810.23628239e-70.58847236e-7820.77598534e-70.66008696e-783-0.17785247e-7-0.86346983e-784-0.24633984e-60.70179584e-785-0.25041147e-70.89426831e-786-0.64923712e-70.30912257e-6870.67462214e-70.75094831e-788-0.12419836e-60.12017218e-62轨道动力学计算的基本数学模型卫星在轨道上运行要受到各种力因素的影响，产生的摄动是多方面的。国内外一些学者对卫星轨道的受摄问题作了详细的研究与分析，尤以澳大利亚的C.Rizos、A.Stolz和美国的H.F.Fliegel等人为代表。统筹考虑精度的需要和时间耗费，通过大量试算，本课题考虑了卫星所受的以下作用力来进行轨道计算（注：时间系统采用TDT时间系统、坐标系统采用J2000.0惯性坐标系），主要包括：地球质心引力F。、除质心外的地球引力七、太阳和月球引力廿、太阳辐射压力F「大气阻力（低卫星轨卫星）、Y轴偏差Fy、地球潮汐附加力Ft。F=F0+Fe+Fn+FA+Ft+七 （2.1）其中地球质心引力F0是最主要的，其次是地球的非质心引力Fe，称为地球非球形摄动力。如果将地球质心引力视为1，地球非球形摄动力可达10-3量级，而其它摄动力则大多在10-6以下。2.1二体问题在惯性系中，卫星运动方程被描述为(2.2)——GM———(2.2)+r(t,r,r)=一 er+r(t,r,r)片3其中.——和—分别表示时刻卫星在惯性系中的位置速度和加速度矢量.其中：r,r和r分力〃表示、时刻卫生在1贝性系中的位置、速度和加速度大量；G和M为分别为引力常数和地球总质量。（3.2）式的第一项为地球质心引力项，称为二体运动，是力模型的主项；第二项为摄动力引起的总摄动项，是t,r,—的函数矢量。

2.2地球非球形引力摄动在地固坐标系中，地球引力位函数作为拉普拉斯方程的解，其非球形部分为:(2.3)(2.4)n- —— — —(2.3)(2.4)U—^^ (CUm+SVm)n—2m—0式中：U_GManPm(sin©)cosmk

n rn+1-GManPm(sin©)sinmkVm_ ean' n rn+1其中：入和6表示单位质点在地固坐标系中的地心经、纬度;ae表示地球平均半径；Pm(sin©)为规格化的缔合勒让德多项式；nC和S为规格化的地球引力位系数；n和m为多项式的阶和次，N为取的最高阶数。非球形引力摄动的求解，按如下步骤进行：a.首先求解Um和Vm。用下列公式逐阶次推算得到，即:如VmLn+1刀a—―arIm如VmLn+1刀a—―arIm{扣'2n+n+11sin©UiV"In-Jmn-1Um-n-1Vm-n-1」Un+1n+1a=TKn+1{Unncos©coskVnncos©Vn+1Ln+1-rn+1Vnn[Unn}sink}(2.5)式中：k=arctg(y),©—arctg(.七, ■ 2n+3Im—1n+1 \(n+m+1)(n-m+1).(n+m)(n-m)aI Jmn-1(2.6)递推方法如下所示:递推初值为:Kn+1=v

n+12n-1rv;3,n=02n+3 nJ-——-,n>0\2n+213U0=竺e,U0=有土U0=竺e,U0=有土sineU0,V0=0,V0=00r i r 0 0 i递推方法采用先求对角线，再按行递推（m不变，n增加）进行。U0♦U0♦U0♦U0♦U0…U0 0 1 2 3 40Ui♦Ui♦Ui♦Ui…Un0U2—U2—U2…U20、U3—U3…U30U44...U41U,11VmLn+1刀aer(n—m+1){(2n+1)sin©[U：1Vmn―J(n+m)rI"}VmUn—1刀Un+1aUnVn递推初值为：n+1Vn+1Ln+1-=伽+1)TrV：］c。帅c°s入不nncos©sinUnJn••一 *..•...........如果采用非规格化的系数U和V进行计算，递推公式为:(2.8)(2.9)(2.10)递推方法同上。此方法完成递推后，还需将U和V规格化，其公式为:U=＜：2n+IU亓 「、(〃—刃"TOC\o"1-5"\h\zUm=-2(2n+I) UmnY (n+m)!nb.计算Um和Vm的偏导数VUm和VVmi.首先计算k■- — —— ——— —Am=结(n—m+i)(n一m+2)(2n+1)/(2n+3)n2Ak. — —Bm=〜(n+m+i)(n+m+2)(2n+1)/(2n+3)n2'Em=、.(n—m+i)(n+m+i)(2n+1)/(2n+3)n其中乘数因子为(2.11)VUmnVVmnaeTFT. FT.AmUm-1—BmUm+1(2.11)VUmnVVmnaeTFT. FT.AmUm-1—BmUm+1—AmVm—1—BmVm+1nn+1 nn+1—EmUmn n+1AmVm—1—BmVm+1

nn+1nn+1AmUm—1+BmUm+1n n+1 n n+1—EmVmnn+1(2.12)frr /i\ftU-m=(—1)mUm

/n n\V—m=(—1)m+1Vmc.求非球形引力摄动引力位而:^N、n._. — - —.VU=^^ (CVUm+SVVm)n=2m=0(2.13)则摄动加速度为:,,—►「=[W]TVU式中：矩阵[W]为地固坐标系转换至惯性坐标系的旋转矩阵。(2.14)2.3日、月摄动考虑卫星的N体影响，只需顾及太阳和月球的引力作用已满足精度要求。N体摄动模型的建立是基于牛顿第二运动定律和万有引力定律。由图2所示，分析卫星的受力可得卫星三体的几何关系的日、月摄动加速度矢量为:其中：A•=r—r—=—zgmN

j=S,LA.为摄A r(y^+j)jA3 r3jj(2.15)动体至卫星的中心距，即矢量K=^/2/2,•巨，ii.计算VUm和VVmn nA的模;jr,r分别为摄动体m和卫星在惯性系中的位置矢量，r,r的计算参见本节附i i i录。这里S和L分别代表太阳和月球。太阳和月球对卫星的摄动影响主要呈长周期变化，且与卫星轨道对太阳和月球的定向有关。2.4太阳直接辐射压摄动照射在卫星上的太阳光，一部分被其吸收转化为热能，另一部分被反射向太空。因此，卫星会受到照射时的辐射力和反射时的反射力的作用，这里统称为直接辐射压摄动。直接辐射压摄动与光压强度、卫星表面的反射系数和光照面积有关。由光压和牛顿第二运动定律建立的直接辐射压对卫星产生的运动加速度矢量为：r=yr=ykA其中：fa2T2)

―sun~as-"a3A2J

esQGsin(a+0.015sin2a)-Gcos(a+0.025sin2a)) (2.16)k=1.013x10-7-2.4x10-9cos2以+1.3x10-9cos4以+4.0x10-10sin6a；电,G对应卫星表面反射系数项，是为补偿模型不足而引入的拟合参数;fa=a=arccosa=1.4959787x1011m为天文单位长度（地球轨道长半径）；T=8.068111241279087x102为人卫时间单位；a为卫星的日心距，即矢量A=r-r的模，r，r分别为卫星和太阳在惯性系中s ss S的位置矢量；y为蚀因子，具体定义为:卫星在地影和月影之外卫星在地球或月球的本影之内A' — , ，■八一 - .1、.1-三卫星在地球或月球的伪本影、半影之内

、A

s其中：4,A分别为太阳的被蚀视面积和视面积。要确定蚀因子，需计算下列诸里：(2.17)A=8：5气=口；A=na2以s<am=sin-i(土)As=sin-i(^m)Ama'、a=sin-i(f)气=气Ji-fa:=(a+20000)/1-fcos2帝八 r•a6广cos-1(M)< rr6 =cos-1(—m△s)、ms A•A式中：A,A为月球和地球的视面积；f,f分别代表太阳和地球的扁率，f=3.352813178x10-3；】为考虑地球大气衰减以及地球扁率效应的有效地球赤道半径;ea=1738000m为月球半径；J=6.96x108m为考虑太阳扁率的太阳有效半径；s4为地心纬度；(2.18)(2.19)(2.20)(2.21)6是地球一卫星一太阳张角；6是月球一卫星一太阳张角。地影和月影判断过程如图3所示。当卫星在地球半影中时：62+a2—a2其中：6= s ee 26当卫星在地球伪本影中时：A-*2当卫星在月球半影中时：其中：6=°m；+a2-amm 26当卫星在月球伪本影中17时：A'=^a2则饨因子为. 4—1max(A，人■e—兀a2s图3地影和月影判断流程人UU-A|J/.丫=1——卜■e—兀a2s图3地影和月影判断流程AS太阳直接辐射摄动对卫星轨道的影响是十分复杂的，它与卫星表面的反射特性、卫星轨道相对太阳的定向以及太阳活动等有关。卫星是由各种不同折射性质的原材料构成的不规则形体，在其运行过程中，太阳对它照射的面积也在不断地改变(它的太阳能翼始终是面向太阳的)。此外，由于太阳活动的变化，所谓太阳常数七也并非常数。因此，给出卫星受太阳直接辐射压摄动的精确模型是很困难的。所以采用较简单的平面模型计算太阳辐射加速度影响。2.5地球固体潮摄动地球并非刚体，它受日月引力作用会产生弹性形变，称为潮汐现象。这种形变使得地球内部物质发生小的变化，随之导致引力位函数产生小的形变位差一一潮汐位。卫星运动的地球固体潮摄动就是潮汐位效应的结果。已知日(或月)对地面点的引力位球谐展开式为：式中：七为日(或月)的质量，,次'分别为引力体至地心距和与地面点的地心角，P(cosV)为n阶勒让德多项式。从上式中排除不产生形变位差的0或1阶项，且只取n=2阶项，可得日月潮汐形变对卫星产生的摄动位为.其中：k为二阶Love数(取k=0.3)。顾及关系式P(cosv)—»3(命「-1］，最后得到卫星的固体潮摄动加速度矢量j2j为：•• (dUVk▽GMa5…<_ , 、r=*I=—Z^-e{［3一15("『)2］『+6(rTr)r} (2.22)t{dr) 2r3r4 j jjj-S,Lj式中：r,r(j-s,l)分别为•,•的单位矢量。jj2.6大气阻力摄动高层大气对卫星的运行将产生阻力，这种阻力对低轨道卫星是主要摄动力之一，但大气密度变化机制非常复杂，不但随高度变化，也与太阳活动、时间、季节、纬度和地磁活动的变化有关。本文采用了静止球型大气密度模型（HP模型）。若只考虑大气分子对卫星表面的法向作用力而忽略其切向作用力，大气阻力使卫星产生的摄动加速度为：••••——=r••••——=r+r1 (A、-2cdp\-\vv—jcRR2DPZ %cos甲V2ARs(32.23)其中：rB为卫星星体部分的大气阻力摄动加速度;r—p为卫星太阳帆板的大气阻力摄动加速度；cd为大气阻力系数；P是大气密度，由HP模型计算得到；V—是卫星相对大气的速度矢量V—=r-oTxr；R RA是卫星参考面积与卫星质量之比；mCDP为太阳帆板的大气阻力系数；Ap为太阳帆板的面积；中为太阳帆板的法向与卫星相对于大气的速度方向的夹角;As为As向量的单位向量。2.7Y轴偏差加速度摄动Y轴偏差加速度主要是GPS卫星的结构失调和卫星体的热辐射而产生的一个附加加速度在星固坐标系Y轴方向上的分量。在设计上，为使卫星的太阳翼以最大面积朝向太阳，两翼的支轴应保持在一条直线上，并要求垂直于卫星和太阳方向的连线，用于控制两翼俯仰的太阳传感器也应完全垂直于卫星翼支轴，卫星的偏航高度控19制应保持正确，但事实上并不完全这样；另一方面，由卫星本身产生的超高温要从Y轴方向的百叶孔排出，这样使处在不稳定状态中的卫星体也会产生结构失调现象。由此导致了Y轴偏差加速度的存在。H.F.Fliegel等人把结构失调的影响表示为［14］T=Gy (2.24)Yya其中：Gy为考虑模型剩差而引入的待估参数；ya为星固坐标系的Y轴在J2000.0惯性坐标系中的方向余弦，由下式决定：2.8巡航姿态控制动力摄动有些卫星在巡航过程中需要保持三轴稳定姿态，需通过姿态控制实现，有的卫星其姿态控制的动力来源于高压气瓶的喷气。这样，在姿态控制的同时也影响了卫星质心的运动。该摄动加速度矢量可以表示为：二二二=.rA=\WKC+CT+Ccosu+SsinuJ (2.25)其中："=o+f；C为卫星在RTN坐标系中姿控动力摄动加速度的常数分量；oC］为卫星在RTN坐标系中姿控动力摄动加速度的时间变化率；C和U为周期项的系数；对全局参数，为由历元时刻起算的相对时间；对弧段相关参数，为观测时刻t在相应弧段内的相对时间。2.9其它摄动影响在轨道运行的卫星除受到上述摄动作用外，还受其它一些摄动的影响，如反照辐射摄动、地球自转形变摄动、广义相对论摄动、海潮摄动、大气潮摄动等。这些摄动影响对卫星轨道摄动非常小，但其计算却要耗费大量的时间片。考虑到2.1节所述的摄动已能满足课题的精度需求和时间消耗的限制，因此，本课题中忽略了其它摄动的影响。附录：日月位置计算太阳位置矢量r计算其中：[D]t为计算历元时刻的平春分点到J2000.0平春分点的转换矩阵，即岁差矩阵。计算当时平春分点的几何平黄经为：平近点角为：平黄赤交角为：以上式中月球位置矢量r”计算其中：^=Q+0.47x10-3sinM'+乎Qsinaj=1 .■系数K、P、Q.和J.的值以及幅角a.的计算列入表5中，表内D为日月平角距，F=L-Q，其计算式为：

表5计算月球位置的系数表j0.02122300.0100200.0008200.0000002D-M'0.01216100.008250-0.002300-0.0001902D0.00332400.0001200.0028300.000000M+180°0.00400000.000000-0.002130-0.0001902F+180°0.00510300.0000900.0004500.0000002M'-2D+180°0.00610000.0004200.0000000.0000002D-M'-M0.00709300.000900-0.0000800.0000002D+M'0.00808000.000560-0.0001200.0000002D-M0.00907200.0003400.0000000.000000M'-M01 0.0006100.0002900.0000000.000000D+PI11 0.0005300.0002800.0000000.000000M'+M+180°210.0002700.000000-0.027240-0.0024102F-2D+180°10.000.00-0.00.002F+M'+300000000000800000180°10.000.000.000.002F-M'+40410020007000000180°10.000.000.000.004D-M'5018001800000000010.000.000.000.004D-2M'6015001200000000010.000.000.000.002D+M-M70140007000000000'+180°10.000.000.000.002D+M+18012000900000000080°10.000.000.000.00M'-D9009000000000000020.000.000.000.00D+M0009000000000000020.000.000.000.002D-M+M10070007000000000'20.000.000.000.002D+2M'2007000900000000020.000.000.000.004D30070008000000000

420.0000000.0000000.0001900.0001802F-2M'20.000.000.000.002F-2D+50000000011000100M20.000.000.000.002D-2F+60000000004700000M20.000.000.000.002F-2D-70000000002300000M'20.000.000.000.002D-2F-80000000000230000M'3轨道计算方法卫星运动方程的解有分析法、数值法和半分析半数值方法等。分析法是将力模型展开取有限项，给出任意时刻解的表达式，通常称之为一般摄动法。其优点是便于定性地给出轨道的特征，如长期项、长周期项影响等。但对摄动力模型处理限制较大，使精度受到影响。数值法即常微分方程的初值解问题，在轨道计算中称之为特别摄动法。其优点是能完整地顾及所选择的力模型、处理简单、计算精度高，是高精度卫星轨道计算最实际有效的方法。经研究，Runge_Kutta型和Cowell型数值积分方法都可用于产生高精度的卫星轨道，尤其是后者通过预报一一校准公式进行迭代计算，收敛快（一般只需迭代2~3次），累积误差影响要比前者小的多。所以这一方法是轨道计算最常用的方法之一。不足的是它为多步型积分器，必须用其它方法为它准备起步值。本课题提供的积分模式有Runge_Kutta4、Runge_Kutta8和高精度的Adams_Cowell方法，高精度计算时采用Runge_Kutta单步积分起步，用Adams_Cowell多步积分求解卫星运动方程，确定卫星在J2000.0惯性坐标系中的位置和速度矢量。3・1Runge_Kutta积分法Runge_Kutta法是一种单步积分方法，公式形式简单，应用方便，将它用于解算卫星运动方程和变分方程也具有很好图4Runge_Kutta积分框图的稳定性。但该方法随着积分式阶数的增 -高，计算右函数的次数也会随之增多，无疑会导致累积误差增加，同时也会降低计算效率。因此，在大规模的轨道计算中，通常只将其作为多步型积分器的起步方法。流程图如图4所示。设有初值问题：!》(t)2((，^) (3.1)[y(七)=七则求解该问题的Runge_Kutta积分公式为：y=y+汹Cf(3.2)-- i=0- -式中：h=t1-1，为积分步长；k为积分式所取的阶数；C.,a,b均为已知的系数，具体参见表6所示。由运动方程和变分方程形成的矩阵常微分方程的初值问题为:上式中：左端y(t)=右函数F=初值y（1）=0[w]

r（t右函数F=初值y（1）=0y，V,r,wb,r,y；3003.2Adams_Cowell积分Adams_Cowell积分适用显始值需要位置项和速度项。积含Adams_Cowell积分适用显始值需要位置项和速度项。积含1阶导的二阶微分方程，初分的预报公式为：V〃广h[七+&.七-']j=0七+1=h2[〃疽工%.七j]I j=0其中：a.,p.为Adams_Cowell积分公式的预报系数，按下列公式计算，见表4.2。校正公式为：Vn+1=h[In+&j孔-j"j=0七+1=h2[II疽£pjynj+1]I j=0式中：a：,p:为Adams_Cowell积分公式的校正系数，按下列公式计算，见表4.2。(3.3)、(3.4)式中的『，⑴和『〃⑴按以下公式递推计算：Im和II皿分别为yn的一次和分和二次和分，按下列公式递推计算：I和II的初值定义为：.、广万一¥snsn 〃“广y空pj吭j=0(3.3)结果参(3.4)结果参(3.3)结果参(3.4)结果参3.3轨道计算图3.3轨道计算图5Adams_Cowell积分流程图就是由一组正确的初始外推出未来的卫星轨型较完备的情况下，初响外推轨道的精度。且值误差传播影响，推算就越大。轨道计算步骤程图值预测卫星轨道实际上值数值积分卫星运动方程，道。显然，在卫星运动力模始值的精确与否将直接影在推算过程中，由于受初始区间越长，积分累积的误差为：图6卫星轨道计算流3.4星历的快速插通过分布式计算提供的星历数据间隔太大，而在实际定位计算中，有时需要间隔为1秒甚至步长更小的精确星历坐标。因此，高精度、快速的精密星历插值就成为一项重要工作。以精度和计算耗费小为出发点，本课题在研究的过程中，分析和比较了以下几种多项式插值方法：拉格朗日多项式插值、牛顿多项式插值、三角函数插值、样条函数插值、切比雪夫多项式插值等。拉格朗日多项式插值的多项式构造和插值时间耗费非常大。其时间耗费可由下式表示：其中：〃表示已知点数；A表示加法运算、M和D分别表示乘、除运算。构造插值多项式还需额外的(n-1)A+nM运算。因此，总运算耗费为：比拉格朗日多项式插值更有效的方法是利用数据均差的牛顿多项式插值[W2]假设有n+1个控制点r,t,,t， 风定义在t处的零阶均差为01n 28 i

f［门=f(t)，在t,t处的1阶均差定义为a=fl,tL区'fllk阶均差定义为：/ / ^'j 0ijt-1.则牛顿插值多项式可表示为：p(t)=a+(t-1){a+(t-1){a+,・・+(t-1){a+(t-1)a},・・}}(4.12)n 0 01 12 n-2n-1 n-1n上述表达式在每级中都由2个加法和1个乘法运算组成，其1时间耗费为：n(2A+M)。与经典的拉格朗日多项式插值相比，牛顿多项式插值经济得多。给定n+1个不同的点t,t,…,t和相应的均差系数a,a，…,a，对任意时刻tet,t］插值多项式

0 1n 0 1 n 0n的结果由下列迭代产生(Horner算法)。上述算法的精度和拉格朗日的精度相同，但时间耗费得到了明显的改善，而且各阶均差可以预先运算好储存起来，需要时直接调用即可。插值点等间隔条件下的插值是牛顿多项式插值的一种特例，可以导出更简洁的形式，而不降低精度。IGS提供的精密星历数据是等间隔的，步长为15分钟。在实际应用中，我们采用牛顿前向差分形式。牛顿前向均差的定义为：同理可推出各高阶前向均差为:(k)(k)(-1)k-j \f(t)Ij) 1+jj=0 3E(t)=A(Ak-1f(t))=Ak-1f(t)-Ak-1f(t)=£

i i i+1 i假设插值多项式的阶数为n，则插值多项式为：(t—t) (t—t)(t—t),・♦(t—t)P(t)=f(t)+I JAf(t)+...+——0 1 n-d-Anf(t) (4.14)n0h0 n!hn 0令，=(t-t"h，则上式可表示为：p(T)=£，一， 、『—，，，，J ，一, j=0\Jj),,,,, 、、r, -用(4.13)式构造插值多项式，再用(4.14)式用作插值，该方法在时间耗费上:Vjf(t0)，(4.13)比Horner算法有一定的改善。与经典的拉格朗日插值方法相比，三次样条多项式插值速度快，但精度低，三角函数多项式插值必须用DFT求解系数，时间耗费大，且这两种方法都将星历数据近似为周期函数，插值精度低，不能满足高精度插值的需要。切比雪夫多项式插值方法在端点处抑制了误差的扩大，但计算量较大。牛顿多项式插值方法的精度高，速度快，充