衛(wèi)星軌道動力學(xué)數(shù)值計算

1、目 錄1星歷計算的時間和坐標(biāo)系統(tǒng)21.1 有關(guān)的時間系統(tǒng)與坐標(biāo)系統(tǒng)21.1.1 時間系統(tǒng)及其換算21.1.2 坐標(biāo)系統(tǒng)及其換算41.2 計算單位和有關(guān)常數(shù)72 軌道動力學(xué)計算的基本數(shù)學(xué)模型122.1 二體問題122.2 地球非球形引力攝動122.3 日、月攝動152.4 太陽直接輻射壓攝動162.5 地球固體潮攝動192.6 大氣阻力攝動192.7 y軸偏差加速度攝動202.8 巡航姿態(tài)控制動力攝動202.9 其它攝動影響21附錄：日月位置計算213 軌道計算方法243.1 runge_kutta積分法243.2 adams_cowell積分253.3 軌道計算273.4 星歷的快速插值284

2、 軌道根數(shù)與位置矢量、速度矢量的關(guān)系324.1 由位置矢量和速度矢量計算軌道根數(shù)324.2 由軌道根數(shù)計算位置矢量和速度矢量331星歷計算的時間和坐標(biāo)系統(tǒng)1.1 有關(guān)的時間系統(tǒng)與坐標(biāo)系統(tǒng)軌道計算過程重要涉及到不同的時間系統(tǒng)和坐標(biāo)系統(tǒng)，下面將空間戰(zhàn)場環(huán)境系統(tǒng)中所涉及到的時間系統(tǒng)和坐標(biāo)系統(tǒng)進行定義，并說明各系統(tǒng)之間的相互關(guān)系。一般情況下，仿真系統(tǒng)采用的是tdt時間系統(tǒng)和j2000地心慣性坐標(biāo)系。1.1.1 時間系統(tǒng)及其換算在軌道計算中，時間是獨立變量。但是，在計算不同的物理量時，卻使用不同的時間系統(tǒng)。例如：在計算恒星時用世界時ut1；定位解算時采用gps時gpst；歲差和章動量的計算采用tdb時等

3、。所以必須清楚各時間系統(tǒng)的定義和各時間系統(tǒng)之間的轉(zhuǎn)換，下面給出各種時間系統(tǒng)的定義及它們之間的轉(zhuǎn)換公式。格林尼治恒星時格林尼治恒星時為春分點對格林尼治平天文子午面的時角。由于歲差、章動原因，它由格林尼治真恒星時（gast）和平恒星時（gmst）之分。兩者的關(guān)系是：其中：為赤經(jīng)章動為自起算至觀測時刻的儒略世紀(jì)數(shù)，即世界時是以平北極（國際習(xí)慣用原點）為統(tǒng)一標(biāo)準(zhǔn)的觀測世界時，是反映地球?qū)嶋H自轉(zhuǎn)的時間，恒星時計算與此有關(guān)。國際原子時時以銫原子基態(tài)兩能級間躍遷輻射的9192631770周所經(jīng)歷的時間作為1秒長的均勻時間，起點在1958年1月1日。國際協(xié)調(diào)時是經(jīng)跳秒修改后的國際原子時，它與世界時的差，觀測紀(jì)

4、錄時間是以此為準(zhǔn)的。質(zhì)心動力學(xué)時（barycentric dynamical time）為相對于太陽質(zhì)心的運動方程給出的歷表、引數(shù)等所用的時間尺度，歲差及章動量的計算是以此為依據(jù)的。地球動力學(xué)時（terrestrial dynamical time）為視地心歷表所用的時間尺度，它具有均勻連續(xù)的特性，衛(wèi)星運動方程就是以此為獨立的時間變量。時間是由系統(tǒng)定義和應(yīng)用的一種時間尺度，于1980年1月6日gpst與utc相等，在此以后由系統(tǒng)主控站密切跟蹤utc以保持高度統(tǒng)一。但gpst不作跳秒修正，因此它與utc具有整秒的差異（1997年1月至6月相差為）。在計算gps衛(wèi)星軌道的初值時將涉及到gpst，g

5、ps精密星歷的參考時間為gpst。以上各時間尺度的相互關(guān)系如下：其中：可從地球自轉(zhuǎn)參數(shù)文件中獲得； ； 。 上式中的t為自年起算至觀測tdb時刻的儒略世紀(jì)數(shù)，即：不同時間系統(tǒng)間的關(guān)系如下圖所示：圖1 幾種時間系統(tǒng)之間的關(guān)系1.1.2 坐標(biāo)系統(tǒng)及其換算衛(wèi)星軌道計算和實際定位解算分別是在j2000.0慣性坐標(biāo)系與wgs-84地固系中進行的，此外，衛(wèi)星加速度計算中還涉及到星固坐標(biāo)系。下面給出與本課題有關(guān)的主要坐標(biāo)系的定義及相互轉(zhuǎn)換關(guān)系。地心慣性系原點：地球質(zhì)心z軸：向北指向年平赤道面（基面）的極點x軸：指向平春分點y軸：符合右手系法則位置矢量：星固坐標(biāo)系原點：衛(wèi)星質(zhì)心z軸：指向衛(wèi)星的天線方向，即指向

6、地心x軸：在軸與太陽構(gòu)成的平面內(nèi)，完成右手系法則y軸：沿衛(wèi)星太陽能翼的支軸位置矢量：星固坐標(biāo)系坐標(biāo)軸在慣性系中的方向余弦分別為：（分別為太陽和衛(wèi)星在地心慣性系中的位置矢量）wgs-84坐標(biāo)系wgs-84為1984年世界大地坐標(biāo)系（wgs為world geodetic system的簡稱），wgs-84的坐標(biāo)定義及其采用的橢球參數(shù)為：原點：地球質(zhì)心z軸：指向bih1984.0定義的協(xié)議地球極（ctp）方向x軸：指向bih1984.0的零子午面和ctp赤道的交點y軸：與x、z軸成右手系地球橢球長半徑：ae=6378137 m地球引力常數(shù)(含大氣層)： gm=3986005108 m3/s2正常化二

7、階帶球諧系數(shù)：=-484.1668510-6地球自轉(zhuǎn)角速度：w=729211510-11 rad/s2地球橢球扁率：f = 1/298.257223563地固坐標(biāo)系與慣性坐標(biāo)系的轉(zhuǎn)換模型a. 慣性坐標(biāo)系地固坐標(biāo)系 模型b. 地固坐標(biāo)系慣性坐標(biāo)系 模型式中：a為極移矩陣； b為自轉(zhuǎn)矩陣； c為章動矩陣； d為歲差矩陣。上述各矩陣的意義及具體定義如下：極移：由于地球不是剛體及其它一些地球物理因素的影響，地球自轉(zhuǎn)軸相對于地球的位置隨時間而變化從而引起觀察者的天頂在天球上的位置發(fā)生變化，稱為極移，矩陣為a：式中：為地極坐標(biāo)，可從地球自轉(zhuǎn)參數(shù)文件中給出的極移值插值得到。自轉(zhuǎn)：即地球公轉(zhuǎn)的同時也在繞自轉(zhuǎn)軸

8、旋轉(zhuǎn)。矩陣b：式中：為格林尼治恒星時 章動：是指外力作用下，地球自轉(zhuǎn)軸在空間運動的短周期擺動部分，即同一瞬間真天極相對平天極的運動，月球?qū)Φ厍蛞Φ淖兓切纬烧聞蝇F(xiàn)象的主要外力作用，其次是太陽。矩陣c：式中：交角章動黃經(jīng)章動其中：都為常數(shù)，可自章動系數(shù)表1中查出。而月亮平近點角：太陽平近點角：月亮平升交角：日月平角矩：月亮軌道對黃道平均升交點的黃經(jīng)：其中： t0為自j2000.0起算至t的儒略世紀(jì)數(shù) 歲差：地球在太陽、月球和行星的引力作用下，地球的自轉(zhuǎn)軸在空間不斷發(fā)生變化，其長期運動稱為歲差，矩陣d：式中：t0的意義同上。旋轉(zhuǎn)矩陣的求法分別為1.2 計算單位和有關(guān)常數(shù)軌道計算采用人衛(wèi)單位系統(tǒng)，

9、具體定義為：長度單位：地球赤道半徑質(zhì)量單位：地球的總質(zhì)量時間單位：在以上人衛(wèi)單位系統(tǒng)中，引力常數(shù)g=1。為完整起見，給出以下常數(shù)：地球赤道半徑：6378137m地球扁率：地球總質(zhì)量：地球自轉(zhuǎn)角速度：地心引力常數(shù)：日心引力常數(shù)：天文單位長度：月球地球質(zhì)量比：光速：引力常數(shù)：太陽常數(shù)：地球引力位系數(shù)采用wgs-84 egm的規(guī)格化值。參見表2。表1 黃經(jīng)和傾角章動序列表引 數(shù)黃經(jīng)章動傾角章動illfd0.00010.0001t0.00010.0001t100001-171996-174.2920258.920000220620.2-8950.53-20201460-240420-200110005

10、-20202-301061-10-10-300070-22-21-2010820-20110009002-22-13187-1.65736-3.110010001426-3.454-0.111012-22-5171.2224-0.6120-12-22217-0.5-950.313002-211290.1-70014200-204801015002-20-22000160200017-0.1001701001-1509018022-22-160.170190-1001-1206020-20021-6030210-12-21-503022200-2140-2023012-2140-2024100-

11、10-400025210-2010002600-22110002701-220-10002801002100029-10011100030012-20-10003100202-2274-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-22

12、290-1204520202-3101304600200260-1047-10201210-10048-10021160-8049100-21-1307050-10221-1005051110-20-7000520120270-30530-1202-70305410222-80305510020600056202-2260-305700021-60305800221-703059102-2160-3060000-21-5030611-100050006220201-503063010-20-40006410-20040006500010-40006611000-3000671020030006

13、81-1202-301069-1-1222-301070-20001-20107130202-3010720-1222-3010731120220-1074-102-21-2010752000120-107610002-201077300002000780021220-1079-1000210-1080100-40-100081-2022210-1082-10242-201083200-40-100084112-2210-108510221-101086-20242-101087-104021000881-10-20100089202-2110-109020222-10009110021-10

14、0092004-22100093302-22100094102-20-10009501201100096-1-102110009700-201-100098002-12-10009901020-100010010-2-20-10001010-1201-1000102110-21-100010310-220-100010420020100010500242-1000106010101000表2 計算地球引力位加速度的引力位系數(shù)（gem-t3）20-0.48416510e-30.021-0.17e-91.19e-9220.24390658e-5-0.14000946e-5300.95720109e

15、-60.0310.20277142e-50.24921712e-6320.90447073e-6-0.61944767e-6330.72034249e-60.14138845e-5400.53952118e-60.041-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-

16、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.

17、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.863469

18、83e-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 軌道動力學(xué)計算的基本數(shù)學(xué)模型衛(wèi)星在軌道上運行要受到各種力因素的影響，產(chǎn)生的攝動是多方面的。國內(nèi)外一些學(xué)者對衛(wèi)星軌道的受攝問題作了詳細的研究與分析，尤以澳大利亞的c.rizos、a.stolz和美國的h.f.fliegel等人為代表。統(tǒng)籌考慮精度的需要和時間耗費，通過大量試算，本課題考慮了衛(wèi)星所受的

19、以下作用力來進行軌道計算（注：時間系統(tǒng)采用tdt時間系統(tǒng)、坐標(biāo)系統(tǒng)采用j2000.0慣性坐標(biāo)系），主要包括：地球質(zhì)心引力f0、除質(zhì)心外的地球引力fe、太陽和月球引力fn、太陽輻射壓力fa、大氣阻力（低衛(wèi)星軌衛(wèi)星）、y軸偏差fy、地球潮汐附加力ft。 (2.1)其中地球質(zhì)心引力f0是最主要的，其次是地球的非質(zhì)心引力fe，稱為地球非球形攝動力。如果將地球質(zhì)心引力視為1，地球非球形攝動力可達10-3量級，而其它攝動力則大多在10-6以下。2.1 二體問題在慣性系中，衛(wèi)星運動方程被描述為 (2.2)其中：和分別表示時刻衛(wèi)星在慣性系中的位置、速度和加速度矢量； 和為分別為引力常數(shù)和地球總質(zhì)量。（3.2）

20、式的第一項為地球質(zhì)心引力項，稱為二體運動，是力模型的主項；第二項為攝動力引起的總攝動項，是的函數(shù)矢量。2.2 地球非球形引力攝動在地固坐標(biāo)系中，地球引力位函數(shù)作為拉普拉斯方程的解，其非球形部分為： (2.3)式中： (2.4)其中：和表示單位質(zhì)點在地固坐標(biāo)系中的地心經(jīng)、緯度； 表示地球平均半徑； 為規(guī)格化的締合勒讓德多項式； 和為規(guī)格化的地球引力位系數(shù)； n和m為多項式的階和次，n為取的最高階數(shù)。非球形引力攝動的求解，按如下步驟進行：a. 首先求解和。用下列公式逐階次推算得到，即： (2.5)式中： (2.6)遞推方法如下所示：遞推初值為： (2.7)遞推方法采用先求對角線，再按行遞推(m不變

21、，n增加)進行。l0 l0 l0 l0 lo m如果采用非規(guī)格化的系數(shù)u和v進行計算，遞推公式為： (2.8)遞推初值為：遞推方法同上。此方法完成遞推后，還需將u和v規(guī)格化，其公式為： (2.9)b. 計算和的偏導(dǎo)數(shù)和 i. 首先計算 (2.10)其中乘數(shù)因子為 (2.11) ii. 計算和 (2.12)而：c. 求非球形引力攝動引力位 (2.13)則攝動加速度為： (2.14)式中：矩陣w為地固坐標(biāo)系轉(zhuǎn)換至慣性坐標(biāo)系的旋轉(zhuǎn)矩陣。2.3 日、月攝動考慮衛(wèi)星的n體影響，只需顧及太陽和月球的引力作用已滿足精度要求。n體攝動模型的建立是基于牛頓第二運動定律和萬有引力定律。由圖2所示，分析衛(wèi)星的受力可

22、得衛(wèi)星圖2 三體的幾何關(guān)系的日、月攝動加速度矢量為： (2.15)其中：，為攝動體至衛(wèi)星的中心距，即矢量的模； 分別為攝動體和衛(wèi)星在慣性系中的位置矢量，的計算參見本節(jié)附錄。 這里s和l分別代表太陽和月球。太陽和月球?qū)πl(wèi)星的攝動影響主要呈長周期變化，且與衛(wèi)星軌道對太陽和月球的定向有關(guān)。2.4 太陽直接輻射壓攝動照射在衛(wèi)星上的太陽光，一部分被其吸收轉(zhuǎn)化為熱能，另一部分被反射向太空。因此，衛(wèi)星會受到照射時的輻射力和反射時的反射力的作用，這里統(tǒng)稱為直接輻射壓攝動。直接輻射壓攝動與光壓強度、衛(wèi)星表面的反射系數(shù)和光照面積有關(guān)。由光壓和牛頓第二運動定律建立的直接輻射壓對衛(wèi)星產(chǎn)生的運動加速度矢量為： (2.1

23、6)其中：； 對應(yīng)衛(wèi)星表面反射系數(shù)項，是為補償模型不足而引入的擬合參數(shù)； ； 為天文單位長度（地球軌道長半徑）； 為人衛(wèi)時間單位；為衛(wèi)星的日心距，即矢量的模，分別為衛(wèi)星和太陽在慣性系中的位置矢量；為蝕因子，具體定義為： (2.17)其中：分別為太陽的被蝕視面積和視面積。要確定蝕因子，需計算下列諸量： (2.18) (2.19) (2.20) (2.21)式中：為月球和地球的視面積； 分別代表太陽和地球的扁率，； 為考慮地球大氣衰減以及地球扁率效應(yīng)的有效地球赤道半徑； =1738000m為月球半徑； 為考慮太陽扁率的太陽有效半徑； 為地心緯度； 是地球衛(wèi)星太陽張角； 是月球衛(wèi)星太陽張角。地影和月

24、影判斷過程如圖3所示。當(dāng)衛(wèi)星在地球半影中時：其中：當(dāng)衛(wèi)星在地球偽本影中時：當(dāng)衛(wèi)星在月球半影中時：其中：當(dāng)衛(wèi)星在月球偽本影中時：則蝕因子為：圖3 地影和月影判斷流程太陽直接輻射攝動對衛(wèi)星軌道的影響是十分復(fù)雜的，它與衛(wèi)星表面的反射特性、衛(wèi)星軌道相對太陽的定向以及太陽活動等有關(guān)。衛(wèi)星是由各種不同折射性質(zhì)的原材料構(gòu)成的不規(guī)則形體，在其運行過程中，太陽對它照射的面積也在不斷地改變（它的太陽能翼始終是面向太陽的）。此外，由于太陽活動的變化，所謂太陽常數(shù)也并非常數(shù)。因此，給出衛(wèi)星受太陽直接輻射壓攝動的精確模型是很困難的。所以采用較簡單的平面模型計算太陽輻射加速度影響。2.5 地球固體潮攝動地球并非剛體，它受

25、日月引力作用會產(chǎn)生彈性形變，稱為潮汐現(xiàn)象。這種形變使得地球內(nèi)部物質(zhì)發(fā)生小的變化，隨之導(dǎo)致引力位函數(shù)產(chǎn)生小的形變位差潮汐位。衛(wèi)星運動的地球固體潮攝動就是潮汐位效應(yīng)的結(jié)果。已知日（或月）對地面點的引力位球諧展開式為：式中：為日（或月）的質(zhì)量，分別為引力體至地心距和與地面點的地心角， 為n階勒讓德多項式。從上式中排除不產(chǎn)生形變位差的0或1階項，且只取n=2階項，可得日月潮汐形變對衛(wèi)星產(chǎn)生的攝動位為：其中：為二階love數(shù)（取）。顧及關(guān)系式，最后得到衛(wèi)星的固體潮攝動加速度矢量為： (2.22)式中：分別為的單位矢量。2.6 大氣阻力攝動高層大氣對衛(wèi)星的運行將產(chǎn)生阻力，這種阻力對低軌道衛(wèi)星是主要攝動力之

26、一，但大氣密度變化機制非常復(fù)雜，不但隨高度變化，也與太陽活動、時間、季節(jié)、緯度和地磁活動的變化有關(guān)。本文采用了靜止球型大氣密度模型（hp模型）。若只考慮大氣分子對衛(wèi)星表面的法向作用力而忽略其切向作用力，大氣阻力使衛(wèi)星產(chǎn)生的攝動加速度為: (32.23)其中：為衛(wèi)星星體部分的大氣阻力攝動加速度； 為衛(wèi)星太陽帆板的大氣阻力攝動加速度；為大氣阻力系數(shù)；是大氣密度，由hp模型計算得到；是衛(wèi)星相對大氣的速度矢量；是衛(wèi)星參考面積與衛(wèi)星質(zhì)量之比； 為太陽帆板的大氣阻力系數(shù)；為太陽帆板的面積；為太陽帆板的法向與衛(wèi)星相對于大氣的速度方向的夾角；為 向量的單位向量。2.7 y軸偏差加速度攝動y軸偏差加速度主要是g

27、ps衛(wèi)星的結(jié)構(gòu)失調(diào)和衛(wèi)星體的熱輻射而產(chǎn)生的一個附加加速度在星固坐標(biāo)系y軸方向上的分量。在設(shè)計上，為使衛(wèi)星的太陽翼以最大面積朝向太陽，兩翼的支軸應(yīng)保持在一條直線上，并要求垂直于衛(wèi)星和太陽方向的連線，用于控制兩翼俯仰的太陽傳感器也應(yīng)完全垂直于衛(wèi)星翼支軸，衛(wèi)星的偏航高度控制應(yīng)保持正確，但事實上并不完全這樣；另一方面，由衛(wèi)星本身產(chǎn)生的超高溫要從y軸方向的百葉孔排出，這樣使處在不穩(wěn)定狀態(tài)中的衛(wèi)星體也會產(chǎn)生結(jié)構(gòu)失調(diào)現(xiàn)象。由此導(dǎo)致了y軸偏差加速度的存在。h.f.fliegel等人把結(jié)構(gòu)失調(diào)的影響表示為14： (2.24)其中：為考慮模型剩差而引入的待估參數(shù)； 為星固坐標(biāo)系的y軸在j2000.0慣性坐標(biāo)系中的

28、方向余弦，由下式?jīng)Q定：2.8 巡航姿態(tài)控制動力攝動有些衛(wèi)星在巡航過程中需要保持三軸穩(wěn)定姿態(tài)，需通過姿態(tài)控制實現(xiàn)，有的衛(wèi)星其姿態(tài)控制的動力來源于高壓氣瓶的噴氣。這樣，在姿態(tài)控制的同時也影響了衛(wèi)星質(zhì)心的運動。該攝動加速度矢量可以表示為: (2.25)其中：； 為衛(wèi)星在rtn坐標(biāo)系中姿控動力攝動加速度的常數(shù)分量； 為衛(wèi)星在rtn坐標(biāo)系中姿控動力攝動加速度的時間變化率； 和為周期項的系數(shù)；對全局參數(shù)，為由歷元時刻起算的相對時間；對弧段相關(guān)參數(shù)，為觀測時刻t在相應(yīng)弧段內(nèi)的相對時間。2.9 其它攝動影響在軌道運行的衛(wèi)星除受到上述攝動作用外，還受其它一些攝動的影響，如反照輻射攝動、地球自轉(zhuǎn)形變攝動、廣義相對

29、論攝動、海潮攝動、大氣潮攝動等。這些攝動影響對衛(wèi)星軌道攝動非常小，但其計算卻要耗費大量的時間片?？紤]到2.1節(jié)所述的攝動已能滿足課題的精度需求和時間消耗的限制，因此，本課題中忽略了其它攝動的影響。附錄：日月位置計算a. 太陽位置矢量計算其中：dt為計算歷元時刻的平春分點到j(luò)2000.0平春分點的轉(zhuǎn)換矩陣，即歲差矩陣。計算當(dāng)時平春分點的幾何平黃經(jīng)為：平近點角為：平黃赤交角為：以上式中b. 月球位置矢量計算其中： 系數(shù)、和的值以及幅角的計算列入表5中，表內(nèi)d為日月平角距，其計算式為： 表5 計算月球位置的系數(shù)表j10.0222300.0100200.0008200.0000002d-m20.011

30、6100.008250-0.002300-0.0001902d30.0032400.0001200.0028300.000000m+18040.0000000.000000-0.002130-0.0001902f+18050.0010300.0000900.0004500.0000002m-2d+18060.0010000.0004200.0000000.0000002d-m-m70.0009300.000900-0.0000800.0000002d+m80.0008000.000560-0.0001200.0000002d-m90.0007200.0003400.0000000.000000

31、m-m100.0006100.0002900.0000000.000000d+pi110.0005300.0002800.0000000.000000m+m+180120.0002700.000000-0.027240-0.0024102f-2d+180130.0000000.000000-0.0000800.0000002f+m+180140.0004100.0002000.0007000.0000002f-m+180150.0001800.0001800.0000000.0000004d-m160.0001500.0001200.0000000.0000004d-2m170.0001400

32、.0000700.0000000.0000002d+m-m+180180.0001200.0000900.0000000.0000002d+m+180190.0000900.0000000.0000000.000000m-d200.0000900.0000000.0000000.000000d+m210.0000700.0000700.0000000.0000002d-m+m220.0000700.0000900.0000000.0000002d+2m230.0000700.0000800.0000000.0000004d240.0000000.0000000.0001900.0001802f

33、-2m250.0000000.0000000.0011000.0001002f-2d+m260.0000000.0000000.0004700.0000002d-2f+m270.0000000.0000000.0002300.0000002f-2d-m280.0000000.0000000.0000230.0000002d-2f-m3 軌道計算方法衛(wèi)星運動方程的解有分析法、數(shù)值法和半分析半數(shù)值方法等。分析法是將力模型展開取有限項，給出任意時刻解的表達式，通常稱之為一般攝動法。其優(yōu)點是便于定性地給出軌道的特征，如長期項、長周期項影響等。但對攝動力模型處理限制較大，使精度受到影響。數(shù)值法即常微分方

34、程的初值解問題，在軌道計算中稱之為特別攝動法。其優(yōu)點是能完整地顧及所選擇的力模型、處理簡單、計算精度高，是高精度衛(wèi)星軌道計算最實際有效的方法。經(jīng)研究，runge_kutta型和cowell型數(shù)值積分方法都可用于產(chǎn)生高精度的衛(wèi)星軌道，尤其是后者通過預(yù)報校準(zhǔn)公式進行迭代計算，收斂快（一般只需迭代23次），累積誤差影響要比前者小的多。所以這一方法是軌道計算最常用的方法之一。不足的是它為多步型積分器，必須用其它方法為它準(zhǔn)備起步值。本課題提供的積分模式有runge_kutta 4、runge_kutta 8和高精度的adams_cowell方法，高精度計算時采用runge_kutta單步積分起步，用ad

35、ams_cowell多步積分求解衛(wèi)星運動方程，確定衛(wèi)星在j2000.0慣性坐標(biāo)系中的位置和速度矢量。 圖4 runge_kutta 積分框圖3.1 runge_kutta積分法runge_kutta法是一種單步積分方法，公式形式簡單，應(yīng)用方便，將它用于解算衛(wèi)星運動方程和變分方程也具有很好的穩(wěn)定性。但該方法隨著積分式階數(shù)的增高，計算右函數(shù)的次數(shù)也會隨之增多，無疑會導(dǎo)致累積誤差增加，同時也會降低計算效率。因此，在大規(guī)模的軌道計算中，通常只將其作為多步型積分器的起步方法。流程圖如圖4所示。設(shè)有初值問題： (3.1)則求解該問題的runge_kutta 積分公式為： (3.2)式中：，為積分步長； k

36、為積分式所取的階數(shù)； ,均為已知的系數(shù)，具體參見表6所示。由運動方程和變分方程形成的矩陣常微分方程的初值問題為：上式中：左端； 右函數(shù)； 初值3.2 adams_cowell積分adams_cowell積分適用顯含1階導(dǎo)的二階微分方程，初始值需要位置項和速度項。積分的預(yù)報公式為： (3.3)其中：為adams_cowell積分公式的預(yù)報系數(shù)，按下列公式計算，結(jié)果參見表4.2。校正公式為： (3.4)式中：為adams_cowell積分公式的校正系數(shù)，按下列公式計算，結(jié)果參見表4.2。(3.3)、(3.4)式中的和按以下公式遞推計算：和分別為的一次和分和二次和分，按下列公式遞推計算：和的初值定義

37、為：用adams_cowell積分解衛(wèi)星運動方程積分流程圖如圖5所示：圖5 adams_cowell積分流程圖3.3 軌道計算預(yù)測衛(wèi)星軌道實際上就是由一組正確的初始值數(shù)值積分衛(wèi)星運動方程，外推出未來的衛(wèi)星軌道。顯然，在衛(wèi)星運動力模型較完備的情況下，初始值的精確與否將直接影響外推軌道的精度。且在推算過程中，由于受初始值誤差傳播影響，推算區(qū)間越長，積分累積的誤差就越大。軌道計算步驟為：圖6 衛(wèi)星軌道計算流程圖3.4 星歷的快速插值通過分布式計算提供的星歷數(shù)據(jù)間隔太大，而在實際定位計算中，有時需要間隔為1秒甚至步長更小的精確星歷坐標(biāo)。因此，高精度、快速的精密星歷插值就成為一項重要工作。以精度和計算耗

38、費小為出發(fā)點，本課題在研究的過程中，分析和比較了以下幾種多項式插值方法：拉格朗日多項式插值、牛頓多項式插值、三角函數(shù)插值、樣條函數(shù)插值、切比雪夫多項式插值等。拉格朗日多項式插值的多項式構(gòu)造和插值時間耗費非常大。其時間耗費可由下式表示：其中：表示已知點數(shù)； a表示加法運算、m和d分別表示乘、除運算。構(gòu)造插值多項式還需額外的運算。因此，總運算耗費為：比拉格朗日多項式插值更有效的方法是利用數(shù)據(jù)均差的牛頓多項式插值12。假設(shè)有個控制點，定義在處的零階均差為，在處的1階均差定義為，k階均差定義為：則牛頓插值多項式可表示為：(4.12)上述表達式在每級中都由2個加法和1個乘法運算組成，其時間耗費為：。與經(jīng)

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