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

目錄TOC\o"1-5"\h\z\o"CurrentDocument"1星歷計算的時間和坐標(biāo)系統(tǒng) 2\o"CurrentDocument"1.1有關(guān)的時間系統(tǒng)與坐標(biāo)系統(tǒng) 2\o"CurrentDocument"1.1.1時間系統(tǒng)及其換算 2\o"CurrentDocument"1.1.2坐標(biāo)系統(tǒng)及其換算 3\o"CurrentDocument"1.2計算單位和有關(guān)常數(shù) 6\o"CurrentDocument"2軌道動力學(xué)計算的基本數(shù)學(xué)模型 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巡航姿態(tài)控制動力攝動 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軌道根數(shù)與位置矢量、速度矢量的關(guān)系 32\o"CurrentDocument"4.1由位置矢量和速度矢量計算軌道根數(shù) 32\o"CurrentDocument"4.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)進(jìn)行定義，并說明各系統(tǒng)之間的相互關(guān)系。一般情況下，仿真系統(tǒng)采用的是TDT時間系統(tǒng)和J2000地心慣性坐標(biāo)系。1.1.1時間系統(tǒng)及其換算在軌道計算中，時間是獨立變量。但是，在計算不同的物理量時，卻使用不同的時間系統(tǒng)。例如：在計算恒星時用世界時UT1；定位解算時采用GPS時GPST；歲差和章動量的計算采用TDB時等。所以必須清楚各時間系統(tǒng)的定義和各時間系統(tǒng)之間的轉(zhuǎn)換，下面給出各種時間系統(tǒng)的定義及它們之間的轉(zhuǎn)換公式。格林尼治恒星時格林尼治恒星時為春分點對格林尼治平天文子午面的時角。由于歲差、章動原因，它由格林尼治真恒星時（GAST）和平恒星時（GMST）之分。兩者的關(guān)系是：其中：Avcos8為赤經(jīng)章動GMST=67310，.54841+（8640184，.812866+876600h）T+0s.093104T2-0..62x10-5T3T為自J2000.0（JD2451545.0）起算至觀測UT1時刻的儒略世紀(jì)數(shù)，即世界時UT1UT1是以平北極（國際習(xí)慣用原點）為統(tǒng)一標(biāo)準(zhǔn)的觀測世界時，是反映地球?qū)嶋H自轉(zhuǎn)的時間，恒星時計算與此有關(guān)。國際原子時TAITAI時以銫原子C133基態(tài)兩 能級間躍遷輻射的9192631770周所經(jīng)歷的時間作為1秒長的均 2勻時間，起點在1958年1月1日0hUT。國際協(xié)調(diào)時UTCUTC是經(jīng)跳秒修改后的國際原子時，它與世界時UT1的差0s.9，觀測紀(jì)錄時間是以此為準(zhǔn)的。質(zhì)心動力學(xué)時tdb(BarycentricDynamicalTime)TDB為相對于太陽質(zhì)心的運(yùn)動方程給出的歷表、引數(shù)等所用的時間尺度，歲差及章動量的計算是以此為依據(jù)的。地球動力學(xué)時tdt(TerrestrialDynamicalTime)T^D^T為視地心歷表所用的時間尺度，它具有均勻連續(xù)的特性，衛(wèi)星運(yùn)動方程就是以此為獨立的時間變量。GPS時間GPSTGPST是由系統(tǒng)定義和應(yīng)用的一種時間尺度，于1980年1月6日0hGPST與UTC相等，在此以后由系統(tǒng)主控站密切跟蹤UTC以保持高度統(tǒng)一。但GPST不作跳秒修正，因此它與UTC具有整秒的差異(1997年1月至6月相差為iis)。在計算GPS衛(wèi)星軌道的初值時將涉及到GPST，GPS精密星歷的參考時間為GPST。以上各時間尺度的相互關(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時刻的儒略世紀(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地固系算中還涉及到星固坐標(biāo)系。下面的定義及相互轉(zhuǎn)換關(guān)系。中進(jìn)行的，此外，衛(wèi)星加速度計給出與本課題有關(guān)的主要坐標(biāo)系算中還涉及到星固坐標(biāo)系。下面的定義及相互轉(zhuǎn)換關(guān)系。J2000.0地心慣性系原點：地球質(zhì)心Z軸：向北指向J2000.0年平赤道面（基面）的極點X軸：指向J2000.0平春分點Y軸：符合右手系法則位置矢量：r星固坐標(biāo)系原點：衛(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ù)為：原點：地球質(zhì)心Z軸：指向BIH1984.0定義的協(xié)議地球極（CTP）方向X軸：指向BIH1984.0的零子午面和CTP赤道的交點Y軸：與X、Z軸成右手系地球橢球長半徑：a=6378137m地球引力常數(shù)（含大氣層）:GM=3986005x108m3/s2正常化二階帶球諧系數(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］ 為歲差矩陣。上述各矩陣的意義及具體定義如下：極移：由于地球不是剛體及其它一些地球物理因素的影響，地球自轉(zhuǎn)軸相對于地球的位置隨時間而變化從而引起觀察者的天頂在天球上的位置發(fā)生變化，稱為極移，矩陣為［A］：式中：x,y為地極坐標(biāo)，可從地球自轉(zhuǎn)參數(shù)文件中給出的極移值插值得到。自轉(zhuǎn)：即地球公轉(zhuǎn)的同時也在繞自轉(zhuǎn)軸旋轉(zhuǎn)。矩陣［B］:式中：9G為格林尼治恒星時章動：是指外力作用下，地球自轉(zhuǎn)軸在空間運(yùn)動的短周期擺動部分，即同一瞬間真天極相對平天極的運(yùn)動，月球?qū)Φ厍蛞Φ淖兓切纬烧聞蝇F(xiàn)象的主要外力作用，其次是太陽。矩陣［C］：式中：屁=如dcos(￡nA) 交角章動J=1 k=1△中=如csin也nA) 黃經(jīng)章動J=1 k=1其中：\dj,n狄都為常數(shù)，可自章動系數(shù)表1中查出。而月亮平近點角：太陽平近點角：月亮平升交角：日月平角矩：月亮軌道對黃道平均升交點 5的黃經(jīng)：

其中：1r=360T0為自J2000.0起算至t的儒略世紀(jì)數(shù)歲差：地球在太陽、月球和行星的引力作用下，地球的自轉(zhuǎn)軸在空間不斷發(fā)生變化，其長期運(yùn)動稱為歲差，矩陣［D］：式中：T0的意義同上。旋轉(zhuǎn)矩陣的求法分別為1.2計算單位和有關(guān)常數(shù)軌道計算采用人衛(wèi)單位系統(tǒng)，具體定義為:長度單位:質(zhì)量單位:地球赤道半徑a長度單位:質(zhì)量單位:e地球的總質(zhì)量Me1時間單位:T=時間單位: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)和傾角章動序列表引數(shù)黃經(jī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計算地球引力位加速度的引力位系數(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軌道動力學(xué)計算的基本數(shù)學(xué)模型衛(wèi)星在軌道上運(yùn)行要受到各種力因素的影響，產(chǎn)生的攝動是多方面的。國內(nèi)外一些學(xué)者對衛(wèi)星軌道的受攝問題作了詳細(xì)的研究與分析，尤以澳大利亞的C.Rizos、A.Stolz和美國的H.F.Fliegel等人為代表。統(tǒng)籌考慮精度的需要和時間耗費，通過大量試算，本課題考慮了衛(wèi)星所受的以下作用力來進(jìn)行軌道計算（注：時間系統(tǒng)采用TDT時間系統(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，稱為地球非球形攝動力。如果將地球質(zhì)心引力視為1，地球非球形攝動力可達(dá)10-3量級，而其它攝動力則大多在10-6以下。2.1二體問題在慣性系中，衛(wèi)星運(yùn)動方程被描述為(2.2)——GM———(2.2)+r(t,r,r)=一 er+r(t,r,r)片3其中.——和—分別表示時刻衛(wèi)星在慣性系中的位置速度和加速度矢量.其中：r,r和r分力〃表示、時刻衛(wèi)生在1貝性系中的位置、速度和加速度大量；G和M為分別為引力常數(shù)和地球總質(zhì)量。（3.2）式的第一項為地球質(zhì)心引力項，稱為二體運(yùn)動，是力模型的主項；第二項為攝動力引起的總攝動項，是t,r,—的函數(shù)矢量。

2.2地球非球形引力攝動在地固坐標(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ì)點在地固坐標(biāo)系中的地心經(jīng)、緯度;ae表示地球平均半徑；Pm(sin?)為規(guī)格化的締合勒讓德多項式；nC和S為規(guī)格化的地球引力位系數(shù)；n和m為多項式的階和次，N為取的最高階數(shù)。非球形引力攝動的求解，按如下步驟進(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)行計算，遞推公式為:(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.計算Um和Vm的偏導(dǎo)數(shù)VUm和VVmi.首先計算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.求非球形引力攝動引力位而:^N、n._. — - —.VU=^^ (CVUm+SVVm)n=2m=0(2.13)則攝動加速度為:,,—?「=[W]TVU式中：矩陣[W]為地固坐標(biāo)系轉(zhuǎn)換至慣性坐標(biāo)系的旋轉(zhuǎn)矩陣。(2.14)2.3日、月攝動考慮衛(wèi)星的N體影響，只需顧及太陽和月球的引力作用已滿足精度要求。N體攝動模型的建立是基于牛頓第二運(yùn)動定律和萬有引力定律。由圖2所示，分析衛(wèi)星的受力可得衛(wèi)星三體的幾何關(guān)系的日、月攝動加速度矢量為:其中：A?=r—r—=—zgmN

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

表5計算月球位置的系數(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軌道計算方法衛(wèi)星運(yùn)動方程的解有分析法、數(shù)值法和半分析半數(shù)值方法等。分析法是將力模型展開取有限項，給出任意時刻解的表達(dá)式，通常稱之為一般攝動法。其優(yōu)點是便于定性地給出軌道的特征，如長期項、長周期項影響等。但對攝動力模型處理限制較大，使精度受到影響。數(shù)值法即常微分方程的初值解問題，在軌道計算中稱之為特別攝動法。其優(yōu)點是能完整地顧及所選擇的力模型、處理簡單、計算精度高，是高精度衛(wèi)星軌道計算最實際有效的方法。經(jīng)研究，Runge_Kutta型和Cowell型數(shù)值積分方法都可用于產(chǎn)生高精度的衛(wèi)星軌道，尤其是后者通過預(yù)報一一校準(zhǔn)公式進(jìn)行迭代計算，收斂快（一般只需迭代2~3次），累積誤差影響要比前者小的多。所以這一方法是軌道計算最常用的方法之一。不足的是它為多步型積分器，必須用其它方法為它準(zhǔn)備起步值。本課題提供的積分模式有Runge_Kutta4、Runge_Kutta8和高精度的Adams_Cowell方法，高精度計算時采用Runge_Kutta單步積分起步，用Adams_Cowell多步積分求解衛(wèi)星運(yùn)動方程，確定衛(wèi)星在J2000.0慣性坐標(biāo)系中的位置和速度矢量。3?1Runge_Kutta積分法Runge_Kutta法是一種單步積分方法，公式形式簡單，應(yīng)用方便，將它用于解算衛(wèi)星運(yùn)動方程和變分方程也具有很好圖4Runge_Kutta積分框圖的穩(wěn)定性。但該方法隨著積分式階數(shù)的增 -高，計算右函數(shù)的次數(shù)也會隨之增多，無疑會導(dǎo)致累積誤差增加，同時也會降低計算效率。因此，在大規(guī)模的軌道計算中，通常只將其作為多步型積分器的起步方法。流程圖如圖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)動方程和變分方程形成的矩陣常微分方程的初值問題為:上式中：左端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積分適用顯始值需要位置項和速度項。積含Adams_Cowell積分適用顯始值需要位置項和速度項。積含1階導(dǎo)的二階微分方程，初分的預(yù)報公式為：V〃廣h[七+&.七-']j=0七+1=h2[〃疽工%.七j]I j=0其中：a.,p.為Adams_Cowell積分公式的預(yù)報系數(shù)，按下列公式計算，見表4.2。校正公式為：Vn+1=h[In+&j孔-j"j=0七+1=h2[II疽￡pjynj+1]I j=0式中：a：,p:為Adams_Cowell積分公式的校正系數(shù)，按下列公式計算，見表4.2。(3.3)、(3.4)式中的『，⑴和『〃⑴按以下公式遞推計算：Im和II皿分別為yn的一次和分和二次和分，按下列公式遞推計算：I和II的初值定義為：.、廣萬一￥snsn 〃“廣y空pj吭j=0(3.3)結(jié)果參(3.4)結(jié)果參(3.3)結(jié)果參(3.4)結(jié)果參3.3軌道計算圖3.3軌道計算圖5Adams_Cowell積分流程圖就是由一組正確的初始外推出未來的衛(wèi)星軌型較完備的情況下，初響外推軌道的精度。且值誤差傳播影響，推算就越大。軌道計算步驟程圖值預(yù)測衛(wèi)星軌道實際上值數(shù)值積分衛(wèi)星運(yùn)動方程，道。顯然，在衛(wèi)星運(yùn)動力模始值的精確與否將直接影在推算過程中，由于受初始區(qū)間越長，積分累積的誤差為：圖6衛(wèi)星軌道計算流3.4星歷的快速插通過分布式計算提供的星歷數(shù)據(jù)間隔太大，而在實際定位計算中，有時需要間隔為1秒甚至步長更小的精確星歷坐標(biāo)。因此，高精度、快速的精密星歷插值就成為一項重要工作。以精度和計算耗費小為出發(fā)點，本課題在研究的過程中，分析和比較了以下幾種多項式插值方法：拉格朗日多項式插值、牛頓多項式插值、三角函數(shù)插值、樣條函數(shù)插值、切比雪夫多項式插值等。拉格朗日多項式插值的多項式構(gòu)造和插值時間耗費非常大。其時間耗費可由下式表示：其中：〃表示已知點數(shù)；A表示加法運(yùn)算、M和D分別表示乘、除運(yùn)算。構(gòu)造插值多項式還需額外的(n-1)A+nM運(yùn)算。因此，總運(yùn)算耗費為：比拉格朗日多項式插值更有效的方法是利用數(shù)據(jù)均差的牛頓多項式插值[W2]假設(shè)有n+1個控制點r,t,,t， 風(fēng)定義在t處的零階均差為01n 28 i

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

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