衛(wèi)星軌道動(dòng)力學(xué)數(shù)值計(jì)算1

目錄TOC\o"1-5"\h\z\o"CurrentDocument"1星歷計(jì)算的時(shí)間和坐標(biāo)系統(tǒng) 2\o"CurrentDocument"1.1有關(guān)的時(shí)間系統(tǒng)與坐標(biāo)系統(tǒng) 2\o"CurrentDocument"1.1.1時(shí)間系統(tǒng)及其換算 2\o"CurrentDocument"1.1.2坐標(biāo)系統(tǒng)及其換算 3\o"CurrentDocument"1.2計(jì)算單位和有關(guān)常數(shù) 6\o"CurrentDocument"2軌道動(dòng)力學(xué)計(jì)算的基本數(shù)學(xué)模型 12\o"CurrentDocument"2.1二體問題 12\o"CurrentDocument"2.2地球非球形引力攝動(dòng) 13\o"CurrentDocument"2.3日、月攝動(dòng) 15\o"CurrentDocument"2.4太陽直接輻射壓攝動(dòng) 16\o"CurrentDocument"2.5地球固體潮攝動(dòng) 18\o"CurrentDocument"2.6大氣阻力攝動(dòng) 19\o"CurrentDocument"2.7Y軸偏差加速度攝動(dòng) 19\o"CurrentDocument"2.8巡航姿態(tài)控制動(dòng)力攝動(dòng) 20\o"CurrentDocument"2.9其它攝動(dòng)影響 20\o"CurrentDocument"附錄：日月位置計(jì)算 21\o"CurrentDocument"3軌道計(jì)算方法 25\o"CurrentDocument"RungeKutta積分法 26\o"CurrentDocument"AdamsCowell積分 263.3軌道計(jì)算 283.4星歷的快速插值 28\o"CurrentDocument"4軌道根數(shù)與位置矢量、速度矢量的關(guān)系 32\o"CurrentDocument"4.1由位置矢量和速度矢量計(jì)算軌道根數(shù) 32\o"CurrentDocument"4.2由軌道根數(shù)計(jì)算位置矢量和速度矢量 331星歷計(jì)算的時(shí)間和坐標(biāo)系統(tǒng)1.1有關(guān)的時(shí)間系統(tǒng)與坐標(biāo)系統(tǒng)軌道計(jì)算過程重要涉及到不同的時(shí)間系統(tǒng)和坐標(biāo)系統(tǒng)，下面將空間戰(zhàn)場環(huán)境系統(tǒng)中所涉及到的時(shí)間系統(tǒng)和坐標(biāo)系統(tǒng)進(jìn)行定義，并說明各系統(tǒng)之間的相互關(guān)系。一般情況下，仿真系統(tǒng)采用的是TDT時(shí)間系統(tǒng)和J2000地心慣性坐標(biāo)系。1.1.1時(shí)間系統(tǒng)及其換算在軌道計(jì)算中，時(shí)間是獨(dú)立變量。但是，在計(jì)算不同的物理量時(shí)，卻使用不同的時(shí)間系統(tǒng)。例如：在計(jì)算恒星時(shí)用世界時(shí)UT1；定位解算時(shí)采用GPS時(shí)GPST；歲差和章動(dòng)量的計(jì)算采用TDB時(shí)等。所以必須清楚各時(shí)間系統(tǒng)的定義和各時(shí)間系統(tǒng)之間的轉(zhuǎn)換，下面給出各種時(shí)間系統(tǒng)的定義及它們之間的轉(zhuǎn)換公式。格林尼治恒星時(shí)格林尼治恒星時(shí)為春分點(diǎn)對格林尼治平天文子午面的時(shí)角。由于歲差、章動(dòng)原因，它由格林尼治真恒星時(shí)（GAST）和平恒星時(shí)（GMST）之分。兩者的關(guān)系是：其中：Avcos8為赤經(jīng)章動(dòng)GMST=67310，.54841+（8640184，.812866+876600h）T+0s.093104T2-0..62x10-5T3T為自J2000.0（JD2451545.0）起算至觀測UT1時(shí)刻的儒略世紀(jì)數(shù)，即世界時(shí)UT1UT1是以平北極（國際習(xí)慣用原點(diǎn)）為統(tǒng)一標(biāo)準(zhǔn)的觀測世界時(shí)，是反映地球?qū)嶋H自轉(zhuǎn)的時(shí)間，恒星時(shí)計(jì)算與此有關(guān)。國際原子時(shí)TAITAI時(shí)以銫原子C133基態(tài)兩 能級間躍遷輻射的9192631770周所經(jīng)歷的時(shí)間作為1秒長的均 2勻時(shí)間，起點(diǎn)在1958年1月1日0hUT。國際協(xié)調(diào)時(shí)UTCUTC是經(jīng)跳秒修改后的國際原子時(shí)，它與世界時(shí)UT1的差0s.9，觀測紀(jì)錄時(shí)間是以此為準(zhǔn)的。質(zhì)心動(dòng)力學(xué)時(shí)tdb(BarycentricDynamicalTime)TDB為相對于太陽質(zhì)心的運(yùn)動(dòng)方程給出的歷表、引數(shù)等所用的時(shí)間尺度，歲差及章動(dòng)量的計(jì)算是以此為依據(jù)的。地球動(dòng)力學(xué)時(shí)tdt(TerrestrialDynamicalTime)T^D^T為視地心歷表所用的時(shí)間尺度，它具有均勻連續(xù)的特性，衛(wèi)星運(yùn)動(dòng)方程就是以此為獨(dú)立的時(shí)間變量。GPS時(shí)間GPSTGPST是由系統(tǒng)定義和應(yīng)用的一種時(shí)間尺度，于1980年1月6日0hGPST與UTC相等，在此以后由系統(tǒng)主控站密切跟蹤UTC以保持高度統(tǒng)一。但GPST不作跳秒修正，因此它與UTC具有整秒的差異(1997年1月至6月相差為iis)。在計(jì)算GPS衛(wèi)星軌道的初值時(shí)將涉及到GPST，GPS精密星歷的參考時(shí)間為GPST。以上各時(shí)間尺度的相互關(guān)系如下：其中：AUT1可從地球自轉(zhuǎn)參數(shù)文件中獲得；AAT=19s+AGPST；ATD=0s.001658sin(v+0.0167sinv),v=6.240040768+628.3019501T(rad)。上式中的T為自j2000.0年起算至觀測TDB時(shí)刻的儒略世紀(jì)數(shù)，即：不同時(shí)間系統(tǒng)間的關(guān)系如下圖所示：圖1幾種時(shí)間系統(tǒng)之間的關(guān)系1.1.2坐標(biāo)系統(tǒng)及其換算衛(wèi)星軌道計(jì)算和實(shí)際定位解算分別是在J2000.0慣性坐標(biāo)系與WGS-84地固系算中還涉及到星固坐標(biāo)系。下面的定義及相互轉(zhuǎn)換關(guān)系。中進(jìn)行的，此外，衛(wèi)星加速度計(jì)給出與本課題有關(guān)的主要坐標(biāo)系算中還涉及到星固坐標(biāo)系。下面的定義及相互轉(zhuǎn)換關(guān)系。J2000.0地心慣性系原點(diǎn)：地球質(zhì)心Z軸：向北指向J2000.0年平赤道面（基面）的極點(diǎn)X軸：指向J2000.0平春分點(diǎn)Y軸：符合右手系法則位置矢量：r星固坐標(biāo)系原點(diǎn)：衛(wèi)星質(zhì)心Z軸：指向衛(wèi)星的天線方向，即指向地心X軸：在軸與太陽構(gòu)成的平面內(nèi)，完成右手系法則Y軸：沿衛(wèi)星太陽能翼的支軸位置矢量：ra星固坐標(biāo)系坐標(biāo)軸伐，a,Z）在慣性系中的方向余弦分別為：（r,r分別為太陽aaa s和衛(wèi)星在J2000.0地心慣性系中的位置矢量）WGS-84坐標(biāo)系WGS-84為1984年世界大地坐標(biāo)系（WGS為WorldGeodeticSystem的簡稱），WGS-84的坐標(biāo)定義及其采用的橢球參數(shù)為：原點(diǎn)：地球質(zhì)心Z軸：指向BIH1984.0定義的協(xié)議地球極（CTP）方向X軸：指向BIH1984.0的零子午面和CTP赤道的交點(diǎn)Y軸：與X、Z軸成右手系地球橢球長半徑：a=6378137m地球引力常數(shù)（含大氣層）:GM=3986005x108m3/s2正?；A帶球諧系數(shù)：F20=-484.16685x10-6地球自轉(zhuǎn)角速度：①=7292115x10-11rad/s2地球橢球扁率：f= 1/298.257223563地固坐標(biāo)系與慣性坐標(biāo)系的 轉(zhuǎn)換模型慣性坐標(biāo)系轉(zhuǎn)換到＞地固坐標(biāo)系模型地固坐標(biāo)系轉(zhuǎn)換到＞慣性坐標(biāo)系模型式中：［A］為極移矩陣；［B］ 為自轉(zhuǎn)矩陣；［C］ 為章動(dòng)矩陣；［D］ 為歲差矩陣。上述各矩陣的意義及具體定義如下：極移：由于地球不是剛體及其它一些地球物理因素的影響，地球自轉(zhuǎn)軸相對于地球的位置隨時(shí)間而變化從而引起觀察者的天頂在天球上的位置發(fā)生變化，稱為極移，矩陣為［A］：式中：x,y為地極坐標(biāo)，可從地球自轉(zhuǎn)參數(shù)文件中給出的極移值插值得到。自轉(zhuǎn)：即地球公轉(zhuǎn)的同時(shí)也在繞自轉(zhuǎn)軸旋轉(zhuǎn)。矩陣［B］:式中：9G為格林尼治恒星時(shí)章動(dòng)：是指外力作用下，地球自轉(zhuǎn)軸在空間運(yùn)動(dòng)的短周期擺動(dòng)部分，即同一瞬間真天極相對平天極的運(yùn)動(dòng)，月球?qū)Φ厍蛞Φ淖兓切纬烧聞?dòng)現(xiàn)象的主要外力作用，其次是太陽。矩陣［C］：式中：屁=如dcos(￡nA) 交角章動(dòng)J=1 k=1△中=如csin也nA) 黃經(jīng)章動(dòng)J=1 k=1其中：\dj,n狄都為常數(shù)，可自章動(dòng)系數(shù)表1中查出。而月亮平近點(diǎn)角：太陽平近點(diǎn)角：月亮平升交角：日月平角矩：月亮軌道對黃道平均升交點(diǎn) 5的黃經(jīng)：

其中：1r=360T0為自J2000.0起算至t的儒略世紀(jì)數(shù)歲差：地球在太陽、月球和行星的引力作用下，地球的自轉(zhuǎn)軸在空間不斷發(fā)生變化，其長期運(yùn)動(dòng)稱為歲差，矩陣［D］：式中：T0的意義同上。旋轉(zhuǎn)矩陣的求法分別為1.2計(jì)算單位和有關(guān)常數(shù)軌道計(jì)算采用人衛(wèi)單位系統(tǒng)，具體定義為:長度單位:質(zhì)量單位:地球赤道半徑a長度單位:質(zhì)量單位:e地球的總質(zhì)量Me1時(shí)間單位:T=時(shí)間單位:T=as=8.068111241279087x102在以上人衛(wèi)單位系統(tǒng)中，引力常數(shù)G=1。為完整起見，給出以下常數(shù):地球赤道半徑：6378137m地球扁率：f=1298.257地球總質(zhì)量：M=5.974x1024kg地球自轉(zhuǎn)角速度：o=7.292115x10-5rad/s地心引力常數(shù)：GM=3.986005x1014m3/s2日心引力常數(shù)：GM=1.32712438x1020m3/s2天文單位長度：a=1.4959787x1011m月球地球質(zhì)量比：GM；GM=0.01230002光速：c=299792458m/s引 力 常 數(shù)： G=6.672x10-11m3/(kg-s2)6太 陽 吊 數(shù)： P=4.5605x10-6kg/(m?s2)地球引力位系數(shù)采用WGS-84EGM的規(guī)格化值。參見表2。表1黃經(jīng)和傾角章動(dòng)序列表引數(shù)黃經(jīng)章動(dòng)傾角章動(dòng)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計(jì)算地球引力位加速度的引力位系數(shù)（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軌道動(dòng)力學(xué)計(jì)算的基本數(shù)學(xué)模型衛(wèi)星在軌道上運(yùn)行要受到各種力因素的影響，產(chǎn)生的攝動(dòng)是多方面的。國內(nèi)外一些學(xué)者對衛(wèi)星軌道的受攝問題作了詳細(xì)的研究與分析，尤以澳大利亞的C.Rizos、A.Stolz和美國的H.F.Fliegel等人為代表。統(tǒng)籌考慮精度的需要和時(shí)間耗費(fèi)，通過大量試算，本課題考慮了衛(wèi)星所受的以下作用力來進(jìn)行軌道計(jì)算（注：時(shí)間系統(tǒng)采用TDT時(shí)間系統(tǒng)、坐標(biāo)系統(tǒng)采用J2000.0慣性坐標(biāo)系），主要包括：地球質(zhì)心引力F。、除質(zhì)心外的地球引力七、太陽和月球引力廿、太陽輻射壓力F「大氣阻力（低衛(wèi)星軌衛(wèi)星）、Y軸偏差Fy、地球潮汐附加力Ft。F=F0+Fe+Fn+FA+Ft+七 （2.1）其中地球質(zhì)心引力F0是最主要的，其次是地球的非質(zhì)心引力Fe，稱為地球非球形攝動(dòng)力。如果將地球質(zhì)心引力視為1，地球非球形攝動(dòng)力可達(dá)10-3量級，而其它攝動(dòng)力則大多在10-6以下。2.1二體問題在慣性系中，衛(wèi)星運(yùn)動(dòng)方程被描述為(2.2)——GM———(2.2)+r(t,r,r)=一 er+r(t,r,r)片3其中.——和—分別表示時(shí)刻衛(wèi)星在慣性系中的位置速度和加速度矢量.其中：r,r和r分力〃表示、時(shí)刻衛(wèi)生在1貝性系中的位置、速度和加速度大量；G和M為分別為引力常數(shù)和地球總質(zhì)量。（3.2）式的第一項(xiàng)為地球質(zhì)心引力項(xiàng)，稱為二體運(yùn)動(dòng)，是力模型的主項(xiàng)；第二項(xiàng)為攝動(dòng)力引起的總攝動(dòng)項(xiàng)，是t,r,—的函數(shù)矢量。

2.2地球非球形引力攝動(dòng)在地固坐標(biāo)系中，地球引力位函數(shù)作為拉普拉斯方程的解，其非球形部分為:(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表示單位質(zhì)點(diǎn)在地固坐標(biāo)系中的地心經(jīng)、緯度;ae表示地球平均半徑；Pm(sin?)為規(guī)格化的締合勒讓德多項(xiàng)式；nC和S為規(guī)格化的地球引力位系數(shù)；n和m為多項(xiàng)式的階和次，N為取的最高階數(shù)。非球形引力攝動(dòng)的求解，按如下步驟進(jì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增加）進(jì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??一 *..?...........如果采用非規(guī)格化的系數(shù)U和V進(jìn)行計(jì)算，遞推公式為:(2.8)(2.9)(2.10)遞推方法同上。此方法完成遞推后，還需將U和V規(guī)格化，其公式為:U=＜：2n+IU亓 「、(〃—刃"TOC\o"1-5"\h\zUm=-2(2n+I) UmnY (n+m)!nb.計(jì)算Um和Vm的偏導(dǎo)數(shù)VUm和VVmi.首先計(jì)算k■- — —— ——— —Am=結(jié)(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其中乘數(shù)因子為(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.求非球形引力攝動(dòng)引力位而:^N、n._. — - —.VU=^^ (CVUm+SVVm)n=2m=0(2.13)則攝動(dòng)加速度為:,,—?「=[W]TVU式中：矩陣[W]為地固坐標(biāo)系轉(zhuǎn)換至慣性坐標(biāo)系的旋轉(zhuǎn)矩陣。(2.14)2.3日、月攝動(dòng)考慮衛(wèi)星的N體影響，只需顧及太陽和月球的引力作用已滿足精度要求。N體攝動(dòng)模型的建立是基于牛頓第二運(yùn)動(dòng)定律和萬有引力定律。由圖2所示，分析衛(wèi)星的受力可得衛(wèi)星三體的幾何關(guān)系的日、月攝動(dòng)加速度矢量為:其中：A?=r—r—=—zgmN

j=S,LA.為攝A r(y^+j)jA3 r3jj(2.15)動(dòng)體至衛(wèi)星的中心距，即矢量K=^/2/2,?巨，ii.計(jì)算VUm和VVmn nA的模;jr,r分別為攝動(dòng)體m和衛(wèi)星在慣性系中的位置矢量，r,r的計(jì)算參見本節(jié)附i i i錄。這里S和L分別代表太陽和月球。太陽和月球?qū)πl(wèi)星的攝動(dòng)影響主要呈長周期變化，且與衛(wèi)星軌道對太陽和月球的定向有關(guān)。2.4太陽直接輻射壓攝動(dòng)照射在衛(wèi)星上的太陽光，一部分被其吸收轉(zhuǎn)化為熱能，另一部分被反射向太空。因此，衛(wèi)星會(huì)受到照射時(shí)的輻射力和反射時(shí)的反射力的作用，這里統(tǒng)稱為直接輻射壓攝動(dòng)。直接輻射壓攝動(dòng)與光壓強(qiáng)度、衛(wèi)星表面的反射系數(shù)和光照面積有關(guān)。由光壓和牛頓第二運(yùn)動(dòng)定律建立的直接輻射壓對衛(wèi)星產(chǎn)生的運(yùn)動(dòng)加速度矢量為：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對應(yīng)衛(wèi)星表面反射系數(shù)項(xiàng)，是為補(bǔ)償模型不足而引入的擬合參數(shù);fa=a=arccosa=1.4959787x1011m為天文單位長度（地球軌道長半徑）；T=8.068111241279087x102為人衛(wèi)時(shí)間單位；a為衛(wèi)星的日心距，即矢量A=r-r的模，r，r分別為衛(wèi)星和太陽在慣性系中s ss S的位置矢量；y為蝕因子，具體定義為:衛(wèi)星在地影和月影之外衛(wèi)星在地球或月球的本影之內(nèi)A' — , ，■八一 - .1、.1-三衛(wèi)星在地球或月球的偽本影、半影之內(nèi)

、A

s其中：4,A分別為太陽的被蝕視面積和視面積。要確定蝕因子，需計(jì)算下列諸里：(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；】為考慮地球大氣衰減以及地球扁率效應(yīng)的有效地球赤道半徑;ea=1738000m為月球半徑；J=6.96x108m為考慮太陽扁率的太陽有效半徑；s4為地心緯度；(2.18)(2.19)(2.20)(2.21)6是地球一衛(wèi)星一太陽張角；6是月球一衛(wèi)星一太陽張角。地影和月影判斷過程如圖3所示。當(dāng)衛(wèi)星在地球半影中時(shí)：62+a2—a2其中：6= s ee 26當(dāng)衛(wèi)星在地球偽本影中時(shí)：A-*2當(dāng)衛(wèi)星在月球半影中時(shí)：其中：6=°m；+a2-amm 26當(dāng)衛(wèi)星在月球偽本影中17時(shí)：A'=^a2則飩因子為. 4—1max(A，人■e—兀a2s圖3地影和月影判斷流程人UU-A|J/.丫=1——卜■e—兀a2s圖3地影和月影判斷流程AS太陽直接輻射攝動(dòng)對衛(wèi)星軌道的影響是十分復(fù)雜的，它與衛(wèi)星表面的反射特性、衛(wèi)星軌道相對太陽的定向以及太陽活動(dòng)等有關(guān)。衛(wèi)星是由各種不同折射性質(zhì)的原材料構(gòu)成的不規(guī)則形體，在其運(yùn)行過程中，太陽對它照射的面積也在不斷地改變(它的太陽能翼始終是面向太陽的)。此外，由于太陽活動(dòng)的變化，所謂太陽常數(shù)七也并非常數(shù)。因此，給出衛(wèi)星受太陽直接輻射壓攝動(dòng)的精確模型是很困難的。所以采用較簡單的平面模型計(jì)算太陽輻射加速度影響。2.5地球固體潮攝動(dòng)地球并非剛體，它受日月引力作用會(huì)產(chǎn)生彈性形變，稱為潮汐現(xiàn)象。這種形變使得地球內(nèi)部物質(zhì)發(fā)生小的變化，隨之導(dǎo)致引力位函數(shù)產(chǎn)生小的形變位差一一潮汐位。衛(wèi)星運(yùn)動(dòng)的地球固體潮攝動(dòng)就是潮汐位效應(yīng)的結(jié)果。已知日(或月)對地面點(diǎn)的引力位球諧展開式為：式中：七為日(或月)的質(zhì)量，,次'分別為引力體至地心距和與地面點(diǎn)的地心角，P(cosV)為n階勒讓德多項(xiàng)式。從上式中排除不產(chǎn)生形變位差的0或1階項(xiàng)，且只取n=2階項(xiàng)，可得日月潮汐形變對衛(wèi)星產(chǎn)生的攝動(dòng)位為.其中：k為二階Love數(shù)(取k=0.3)。顧及關(guān)系式P(cosv)—?3(命「-1］，最后得到衛(wèi)星的固體潮攝動(dòng)加速度矢量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大氣阻力攝動(dòng)高層大氣對衛(wèi)星的運(yùn)行將產(chǎn)生阻力，這種阻力對低軌道衛(wèi)星是主要攝動(dòng)力之一，但大氣密度變化機(jī)制非常復(fù)雜，不但隨高度變化，也與太陽活動(dòng)、時(shí)間、季節(jié)、緯度和地磁活動(dòng)的變化有關(guān)。本文采用了靜止球型大氣密度模型（HP模型）。若只考慮大氣分子對衛(wèi)星表面的法向作用力而忽略其切向作用力，大氣阻力使衛(wèi)星產(chǎn)生的攝動(dòng)加速度為：????——=r????——=r+r1 (A、-2cdp\-\vv—jcRR2DPZ %cos甲V2ARs(32.23)其中：rB為衛(wèi)星星體部分的大氣阻力攝動(dòng)加速度;r—p為衛(wèi)星太陽帆板的大氣阻力攝動(dòng)加速度；cd為大氣阻力系數(shù)；P是大氣密度，由HP模型計(jì)算得到；V—是衛(wèi)星相對大氣的速度矢量V—=r-oTxr；R RA是衛(wèi)星參考面積與衛(wèi)星質(zhì)量之比；mCDP為太陽帆板的大氣阻力系數(shù)；Ap為太陽帆板的面積；中為太陽帆板的法向與衛(wèi)星相對于大氣的速度方向的夾角;As為As向量的單位向量。2.7Y軸偏差加速度攝動(dòng)Y軸偏差加速度主要是GPS衛(wèi)星的結(jié)構(gòu)失調(diào)和衛(wèi)星體的熱輻射而產(chǎn)生的一個(gè)附加加速度在星固坐標(biāo)系Y軸方向上的分量。在設(shè)計(jì)上，為使衛(wèi)星的太陽翼以最大面積朝向太陽，兩翼的支軸應(yīng)保持在一條直線上，并要求垂直于衛(wèi)星和太陽方向的連線，用于控制兩翼俯仰的太陽傳感器也應(yīng)完全垂直于衛(wèi)星翼支軸，衛(wèi)星的偏航高度控19制應(yīng)保持正確，但事實(shí)上并不完全這樣；另一方面，由衛(wèi)星本身產(chǎn)生的超高溫要從Y軸方向的百葉孔排出，這樣使處在不穩(wěn)定狀態(tài)中的衛(wèi)星體也會(huì)產(chǎn)生結(jié)構(gòu)失調(diào)現(xiàn)象。由此導(dǎo)致了Y軸偏差加速度的存在。H.F.Fliegel等人把結(jié)構(gòu)失調(diào)的影響表示為［14］T=Gy (2.24)Yya其中：Gy為考慮模型剩差而引入的待估參數(shù)；ya為星固坐標(biāo)系的Y軸在J2000.0慣性坐標(biāo)系中的方向余弦，由下式?jīng)Q定：2.8巡航姿態(tài)控制動(dòng)力攝動(dòng)有些衛(wèi)星在巡航過程中需要保持三軸穩(wěn)定姿態(tài)，需通過姿態(tài)控制實(shí)現(xiàn)，有的衛(wèi)星其姿態(tài)控制的動(dòng)力來源于高壓氣瓶的噴氣。這樣，在姿態(tài)控制的同時(shí)也影響了衛(wèi)星質(zhì)心的運(yùn)動(dòng)。該攝動(dòng)加速度矢量可以表示為：二二二=.rA=\WKC+CT+Ccosu+SsinuJ (2.25)其中："=o+f；C為衛(wèi)星在RTN坐標(biāo)系中姿控動(dòng)力攝動(dòng)加速度的常數(shù)分量；oC］為衛(wèi)星在RTN坐標(biāo)系中姿控動(dòng)力攝動(dòng)加速度的時(shí)間變化率；C和U為周期項(xiàng)的系數(shù)；對全局參數(shù)，為由歷元時(shí)刻起算的相對時(shí)間；對弧段相關(guān)參數(shù)，為觀測時(shí)刻t在相應(yīng)弧段內(nèi)的相對時(shí)間。2.9其它攝動(dòng)影響在軌道運(yùn)行的衛(wèi)星除受到上述攝動(dòng)作用外，還受其它一些攝動(dòng)的影響，如反照輻射攝動(dòng)、地球自轉(zhuǎn)形變攝動(dòng)、廣義相對論攝動(dòng)、海潮攝動(dòng)、大氣潮攝動(dòng)等。這些攝動(dòng)影響對衛(wèi)星軌道攝動(dòng)非常小，但其計(jì)算卻要耗費(fèi)大量的時(shí)間片?？紤]到2.1節(jié)所述的攝動(dòng)已能滿足課題的精度需求和時(shí)間消耗的限制，因此，本課題中忽略了其它攝動(dòng)的影響。附錄：日月位置計(jì)算太陽位置矢量r計(jì)算其中：[D]t為計(jì)算歷元時(shí)刻的平春分點(diǎn)到J2000.0平春分點(diǎn)的轉(zhuǎn)換矩陣，即歲差矩陣。計(jì)算當(dāng)時(shí)平春分點(diǎn)的幾何平黃經(jīng)為：平近點(diǎn)角為：平黃赤交角為：以上式中月球位置矢量r”計(jì)算其中：^=Q+0.47x10-3sinM'+乎Qsinaj=1 .■系數(shù)K、P、Q.和J.的值以及幅角a.的計(jì)算列入表5中，表內(nèi)D為日月平角距，F(xiàn)=L-Q，其計(jì)算式為：

表5計(jì)算月球位置的系數(shù)表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軌道計(jì)算方法衛(wèi)星運(yùn)動(dòng)方程的解有分析法、數(shù)值法和半分析半數(shù)值方法等。分析法是將力模型展開取有限項(xiàng)，給出任意時(shí)刻解的表達(dá)式，通常稱之為一般攝動(dòng)法。其優(yōu)點(diǎn)是便于定性地給出軌道的特征，如長期項(xiàng)、長周期項(xiàng)影響等。但對攝動(dòng)力模型處理限制較大，使精度受到影響。數(shù)值法即常微分方程的初值解問題，在軌道計(jì)算中稱之為特別攝動(dòng)法。其優(yōu)點(diǎn)是能完整地顧及所選擇的力模型、處理簡單、計(jì)算精度高，是高精度衛(wèi)星軌道計(jì)算最實(shí)際有效的方法。經(jīng)研究，Runge_Kutta型和Cowell型數(shù)值積分方法都可用于產(chǎn)生高精度的衛(wèi)星軌道，尤其是后者通過預(yù)報(bào)一一校準(zhǔn)公式進(jìn)行迭代計(jì)算，收斂快（一般只需迭代2~3次），累積誤差影響要比前者小的多。所以這一方法是軌道計(jì)算最常用的方法之一。不足的是它為多步型積分器，必須用其它方法為它準(zhǔn)備起步值。本課題提供的積分模式有Runge_Kutta4、Runge_Kutta8和高精度的Adams_Cowell方法，高精度計(jì)算時(shí)采用Runge_Kutta單步積分起步，用Adams_Cowell多步積分求解衛(wèi)星運(yùn)動(dòng)方程，確定衛(wèi)星在J2000.0慣性坐標(biāo)系中的位置和速度矢量。3?1Runge_Kutta積分法Runge_Kutta法是一種單步積分方法，公式形式簡單，應(yīng)用方便，將它用于解算衛(wèi)星運(yùn)動(dòng)方程和變分方程也具有很好圖4Runge_Kutta積分框圖的穩(wěn)定性。但該方法隨著積分式階數(shù)的增 -高，計(jì)算右函數(shù)的次數(shù)也會(huì)隨之增多，無疑會(huì)導(dǎo)致累積誤差增加，同時(shí)也會(huì)降低計(jì)算效率。因此，在大規(guī)模的軌道計(jì)算中，通常只將其作為多步型積分器的起步方法。流程圖如圖4所示。設(shè)有初值問題：!》(t)2((，^) (3.1)[y(七)=七則求解該問題的Runge_Kutta積分公式為：y=y+洶Cf(3.2)-- i=0- -式中：h=t1-1，為積分步長；k為積分式所取的階數(shù)；C.,a,b均為已知的系數(shù)，具體參見表6所示。由運(yùn)動(dòng)方程和變分方程形成的矩陣常微分方程的初值問題為:上式中：左端y(t)=右函數(shù)F=初值y（1）=0[w]

r（t右函數(shù)F=初值y（1）=0y，V,r,wb,r,y；3003.2Adams_Cowell積分Adams_Cowell積分適用顯始值需要位置項(xiàng)和速度項(xiàng)。積含Adams_Cowell積分適用顯始值需要位置項(xiàng)和速度項(xiàng)。積含1階導(dǎo)的二階微分方程，初分的預(yù)報(bào)公式為：V〃廣h[七+&.七-']j=0七+1=h2[〃疽工%.七j]I j=0其中：a.,p.為Adams_Cowell積分公式的預(yù)報(bào)系數(shù)，按下列公式計(jì)算，見表4.2。校正公式為：Vn+1=h[In+&j孔-j"j=0七+1=h2[II疽￡pjynj+1]I j=0式中：a：,p:為Adams_Cowell積分公式的校正系數(shù)，按下列公式計(jì)算，見表4.2。(3.3)、(3.4)式中的『，⑴和『〃⑴按以下公式遞推計(jì)算：Im和II皿分別為yn的一次和分和二次和分，按下列公式遞推計(jì)算：I和II的初值定義為：.、廣萬一￥snsn 〃“廣y空pj吭j=0(3.3)結(jié)果參(3.4)結(jié)果參(3.3)結(jié)果參(3.4)結(jié)果參3.3軌道計(jì)算圖3.3軌道計(jì)算圖5Adams_Cowell積分流程圖就是由一組正確的初始外推出未來的衛(wèi)星軌型較完備的情況下，初響外推軌道的精度。且值誤差傳播影響，推算就越大。軌道計(jì)算步驟程圖值預(yù)測衛(wèi)星軌道實(shí)際上值數(shù)值積分衛(wèi)星運(yùn)動(dòng)方程，道。顯然，在衛(wèi)星運(yùn)動(dòng)力模始值的精確與否將直接影在推算過程中，由于受初始區(qū)間越長，積分累積的誤差為：圖6衛(wèi)星軌道計(jì)算流3.4星歷的快速插通過分布式計(jì)算提供的星歷數(shù)據(jù)間隔太大，而在實(shí)際定位計(jì)算中，有時(shí)需要間隔為1秒甚至步長更小的精確星歷坐標(biāo)。因此，高精度、快速的精密星歷插值就成為一項(xiàng)重要工作。以精度和計(jì)算耗費(fèi)小為出發(fā)點(diǎn)，本課題在研究的過程中，分析和比較了以下幾種多項(xiàng)式插值方法：拉格朗日多項(xiàng)式插值、牛頓多項(xiàng)式插值、三角函數(shù)插值、樣條函數(shù)插值、切比雪夫多項(xiàng)式插值等。拉格朗日多項(xiàng)式插值的多項(xiàng)式構(gòu)造和插值時(shí)間耗費(fèi)非常大。其時(shí)間耗費(fèi)可由下式表示：其中：〃表示已知點(diǎn)數(shù)；A表示加法運(yùn)算、M和D分別表示乘、除運(yùn)算。構(gòu)造插值多項(xiàng)式還需額外的(n-1)A+nM運(yùn)算。因此，總運(yùn)算耗費(fèi)為：比拉格朗日多項(xiàng)式插值更有效的方法是利用數(shù)據(jù)均差的牛頓多項(xiàng)式插值[W2]假設(shè)有n+1個(gè)控制點(diǎn)r,t,,t， 風(fēng)定義在t處的零階均差為01n 28 i

f［門=f(t)，在t,t處的1階均差定義為a=fl,tL區(qū)'fllk階均差定義為：/ / ^'j 0ijt-1.則牛頓插值多項(xiàng)式可表示為：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上述表達(dá)式在每級中都由2個(gè)加法和1個(gè)乘法運(yùn)算組成，其1時(shí)間耗費(fèi)為：n(2A+M)。與經(jīng)典的拉格朗日多項(xiàng)式插值相比，牛頓多項(xiàng)式插值經(jīng)濟(jì)得多。給定n+1個(gè)不同的點(diǎn)t,t,…,t和相應(yīng)的均差系數(shù)a,a，…,a，對任意時(shí)刻tet,t］插值多項(xiàng)式

0 1n 0 1 n 0n的結(jié)果由下列迭代產(chǎn)生(Horner算法)。上述算法的精度和拉格朗日的精度相同，但時(shí)間耗費(fèi)得到了明顯的改善，而且各階均差可以預(yù)先運(yùn)算好儲(chǔ)存起來，需要時(shí)直接調(diào)用即可。插值點(diǎn)等間隔條件下的插值是牛頓多項(xiàng)式插值的一種特例，可以導(dǎo)出更簡潔的形式，而不降低精度。IGS提供的精密星歷數(shù)據(jù)是等間隔的，步長為15分鐘。在實(shí)際應(yīng)用中，我們采用牛頓前向差分形式。牛頓前向均差的定義為：同理可推出各高階前向均差為:(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假設(shè)插值多項(xiàng)式的階數(shù)為n，則插值多項(xiàng)式為：(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)式構(gòu)造插值多項(xiàng)式，再用(4.14)式用作插值，該方法在時(shí)間耗費(fèi)上:Vjf(t0)，(4.13)比Horner算法有一定的改善。與經(jīng)典的拉格朗日插值方法相比，三次樣條多項(xiàng)式插值速度快，但精度低，三角函數(shù)多項(xiàng)式插值必須用DFT求解系數(shù)，時(shí)間耗費(fèi)大，且這兩種方法都將星歷數(shù)據(jù)近似為周期函數(shù)，插值精度低，不能滿足高精度插值的需要。切比雪夫多項(xiàng)式插值方法在端點(diǎn)處抑制了誤差的擴(kuò)大，但計(jì)算量較大。牛頓多項(xiàng)式插值方法的精度高，速度快，充