鄭峰椿 湯靖師（ 南京大學(xué)天文與空間科學(xué)學(xué)院 南京 00）

（2 現(xiàn)代天文與天體物理教育部重點實驗室 南京 210023）

（3 南京大學(xué)空間環(huán)境與航天動力學(xué)研究所 南京 210023）

1 引言

分析法軌道預(yù)報由于具有較大的速度優(yōu)勢而被廣泛應(yīng)用于近地天體動力學(xué)研究中, 在構(gòu)造軌道分析解時除了計算二體引力還需要考慮地球非球形引力、日月第三體引力、大氣阻力、太陽光壓等攝動力的作用, 其中地球非球形引力攝動影響最大, 與二體引力相比量級約為O（10-3）. 地球非球形引力攝動函數(shù)一般以地固坐標(biāo)系中衛(wèi)星的球坐標(biāo)進行計算, 若將球坐標(biāo)轉(zhuǎn)換為衛(wèi)星軌道根數(shù), 那么函數(shù)表達式中可分離出只包含軌道傾角的傾角函數(shù), 而傾角函數(shù)是田諧項攝動解的重要組成因子.

在近地天體動力學(xué)領(lǐng)域, Kaula[1]最先提出傾角函數(shù)的概念并記為Flmp（I）, 表達式為關(guān)于sinI和cosI的三重求和多項式, 其中I是軌道傾角,l、m、p是傾角函數(shù)的3個指標(biāo), 其取值范圍為2 ≤l≤L,1 ≤m≤l, 0 ≤p≤l,L即地球引力場階次. 之后Izsak[2]給出更加簡潔的公式Kjnm（I）, 指標(biāo)n、j分別對應(yīng)了l、p, 只需要對sin和cos作一重求和.在此基礎(chǔ)上Allan[3]改寫了相應(yīng)符號而保持公式結(jié)構(gòu)不變, 得到迄今為止最泛用的傾角函數(shù)表達式,公式如下:其中v為求和指標(biāo).

盡管此公式已足夠簡潔,但在田諧項攝動解中,傾角函數(shù)需要多重求和, 在引力場階次較高時計算效率低, 且會出現(xiàn)損失精度的情況. Gooding[4–5]提出另一種遞推計算方法, 將歸一化傾角函數(shù)分解為兩個函數(shù)的乘積, 公式如下:

其中是關(guān)于sinI和cosI的多項式, 而則是關(guān)于sin和cos的多項式, 并可以通過三項遞推公式計算, 即Gooding遞推公式[4–5].與關(guān)系為:

其中Nlm為歸一化因子,k=l -2p. 值得注意的是傾角函數(shù)第3個指標(biāo)無論使用p還是k都可以, 本質(zhì)上沒有差別, 而Gooding等[5]指出雖然k取值不連續(xù), 但傾角函數(shù)關(guān)于k的正負具有對稱性, 在遞推計算時可借此減少計算量. 此外Gooding等[5]對遞推算法做出進一步優(yōu)化, 并將分離出不含傾角的因子以遞推方法計算, 解決了在高階次情況下出現(xiàn)上溢和下溢的問題. 同時Gooding等[5]還給出了傾角函數(shù)1階、2階偏導(dǎo)數(shù)的計算公式, 并證明了按照（2）式計算傾角函數(shù)時效率最高且具有高穩(wěn)定性, 因此Gooding遞推公式得到了廣泛應(yīng)用, 成為其他計算方法的衡量標(biāo)準.

在文獻[4–5]中Gooding等人并沒有詳細給出遞推公式的證明, 吳連大等[6–7]分別利用超幾何級數(shù)的遞推關(guān)系和Jacobi多項式的遞推關(guān)系, 完成了公式推導(dǎo), 并指出其實質(zhì)為Jacobi多項式的遞推. 之后吳連大等[8–9]在Gooding遞推方法的基礎(chǔ)上提出計算傾角函數(shù)的改進Gooding方法, 按照文中給出的標(biāo)準遞推過程計算, 而對于V klm則將Gooding[4]提出的表達式進行改寫并直接計算, 同時避免了出現(xiàn)上溢和下溢的情況. 測試結(jié)果表明改進Gooding方法計算速度更快、計算精度和穩(wěn)定性也有所提高, 另一方面傾角函數(shù)的計算程序也更加簡潔, 因此在分析法軌道預(yù)報過程中可以使用改進Gooding方法計算傾角函數(shù).

在分析法軌道預(yù)報過程中地球非球形引力攝動解的計算占據(jù)了較大比重, 而田諧項攝動分析解中需要大批量計算傾角函數(shù), 因此傾角函數(shù)的計算速度很大程度上影響了分析法預(yù)報耗時. 本文第2節(jié)以改進Gooding方法為基礎(chǔ), 對附加的Fortran程序進行改進, 進一步提升了大批量計算傾角函數(shù)的效率,為了方便,本文稱此方法為網(wǎng)格Gooding法.

一般情況下, 對地球衛(wèi)星作軌道預(yù)報時預(yù)報弧段不會太長, 預(yù)報過程中軌道傾角的變化量?I很小, 因此可以將傾角函數(shù)在初始傾角I0處作泰勒展開, 以1階或2階展開式計算傾角函數(shù). 傾角函數(shù)展開式為關(guān)于?I的線性或2次多項式, 與改進Gooding法相比,一方面計算效率得到大幅提升,但另一方面展開式存在一定的截斷誤差, 從而導(dǎo)致分析法預(yù)報精度受到影響. 如果展開式誤差帶來的軌道預(yù)報結(jié)果偏差在一定范圍內(nèi), 那么在分析法預(yù)報過程中可以使用展開式計算傾角函數(shù), 從而很大程度地提高預(yù)報速度, 對大量目標(biāo)作軌道預(yù)報時利用此方法可極大減輕工作量. 本文第3節(jié)分析了不同初始傾角、地球引力場階次、傾角變化量對傾角函數(shù)展開式誤差的影響, 比較了分別精確至1階、2階時展開式的計算速度以及精度. 第4節(jié)設(shè)計若干算例, 在分析法軌道預(yù)報中分別使用網(wǎng)格Gooding法、1階、2階展開式計算傾角函數(shù), 以改進Gooding方法為基準, 分析預(yù)報星歷偏差以及預(yù)報耗時,其中預(yù)報星歷為J2000平赤道坐標(biāo)系下的直角坐標(biāo).

2 傾角函數(shù)遞推算法介紹與程序優(yōu)化

由于吳連大等[8–9]提出的改進Gooding方法具有計算速度快、計算精度和穩(wěn)定性高、程序簡潔的優(yōu)點, 故可以直接應(yīng)用到分析法軌道預(yù)報中, 而在實際應(yīng)用中發(fā)現(xiàn)遞推程序有進一步提升速度的空間. 根據(jù)田諧項攝動解公式可知預(yù)報過程中需要計算大批量傾角函數(shù), 參考劉林[10]書中公式:

其中σ（2）s（t）表示軌道根數(shù)的2階短周期項, ?σlmpq中包含傾角函數(shù)Flmp及其1階偏導(dǎo)數(shù)F′lmp, 劉林[10]給出了田諧項攝動解的具體公式, 由于公式較長,本文不再羅列. 需要注意的是改進Gooding法計算結(jié)果為歸一化傾角函數(shù), 而常用的分析解中傾角函數(shù)為非歸一化的,所以實際計算時需對其進行轉(zhuǎn)換.

根據(jù)（4）式, 由l≤L可知傾角函數(shù)計算次數(shù)取決于地球引力場階次, 遞推次數(shù)N可由下式計算,

這表示至少需要循環(huán)調(diào)用N次遞推程序, 而根據(jù)改進Gooding法原理可知在遞推過程中可計算出若干傾角函數(shù), 如果能夠記錄下這些值, 那么在循環(huán)過程完全可以跳過這次計算,從而減少程序調(diào)用次數(shù).當(dāng)引力場階次較高時, 將節(jié)省不少傾角函數(shù)的計算時間.

同Gooding遞推方法一樣, 吳連大等[8–9]將傾角函數(shù)分解為（2）式的形式, 按照標(biāo)準過程遞推計算, 不同的是改寫了表達式并直接計算, 因此在批量計算時只需考慮的程序改進, 其遞推公式為:

其中x= cosI, 0 ≤J≤l, 公式沿著l遞推, 遞推初值為其中q= max（m,k）. 該遞推關(guān)系成立的條件為k≥0, 在k <0時可以利用傾角函數(shù)的對稱性直接計算, 如（7）式所示:

根據(jù)（6）式可知, 傾角函數(shù)在沿著l方向遞推過程中可得到中間值（q≤z≤l）, 由（2）式和（4）式可知這些中間值正是田諧項攝動分析解計算中所需要的, 而對遞推程序的改進原理為將其記錄下來以避免在循環(huán)過程中重復(fù)計算.

遞推程序改進的具體方法為將傾角函數(shù)以3維數(shù)組形式作為函數(shù)接口. 在吳連大等[8–9]給出的程序中, 函數(shù)輸入接口為軌道傾角和3個指標(biāo)l、m、k, 輸出接口為傾角函數(shù)及其偏導(dǎo)數(shù). 修改之后輸入接口為傾角、地球引力場階次L, 輸出接口為3維數(shù)組形式的傾角函數(shù)及其偏導(dǎo)數(shù), 在函數(shù)內(nèi)部同樣為3維數(shù)組.該算法的原理很簡單,規(guī)定2 ≤l≤L,1 ≤m≤l, 0 ≤p≤l, 在沿著l、m、p進行3重循環(huán)過程中將p轉(zhuǎn)換為k, 以改進Gooding法計算傾角函數(shù), 并將、分別賦給3維數(shù)組A[l][m][p]和F[l][m][p].每次計算之前對A[l][m][p]進行檢驗,如果已經(jīng)被賦值, 則跳過該過程, 并且在計算傾角函數(shù)時直接引用A[l][m][p]. 改進之后仍需計算N次, 而的計算次數(shù)減少, 因此傾角函數(shù)的計算耗時有所縮減, 且隨著階次的增加, 耗時減少比例逐漸增大. 由于該算法本質(zhì)為使用改進Gooding法計算傾角函數(shù)并在3維網(wǎng)格中賦值, 因此精度并未改變, 為了方便, 本文中稱其為網(wǎng)格Gooding法.

吳連大等[8–9]測試了1–200階傾角函數(shù)的計算時間, 與Gooding方法相比, 使用此方法時效率提升了41%. 本文設(shè)計若干算例, 對改進后程序的計算效率進行評估, 考慮到在分析法預(yù)報中地球引力場取50階即可獲得高精度星歷, 因此L選取10、50, 每組設(shè)置10個傾角, 選取0.1?、10?、20?、···、90?, 分別按照改進Gooding法和網(wǎng)格Gooding法批量計算傾角函數(shù)及其偏導(dǎo)數(shù), 且重復(fù)計算10次, 并記錄平均耗時, 其中兩種算法都以C++語言實現(xiàn).表1給出計算機配置以及編譯設(shè)置, 所有算例均以Release模式進行測試. 相比之下, 隨著L增大, 網(wǎng)格Gooding法的計算效率逐漸提高,L取10、50時計算耗時分別縮短了18.8%、23.8%（若將此方法應(yīng)用到分析法預(yù)報中, 預(yù)報耗時可參考本文第4節(jié)中的表5）.

3 傾角函數(shù)的展開式分析

3.1 傾角函數(shù)展開式

根據(jù)（1）式可知在3個指標(biāo)不變時傾角函數(shù)為關(guān)于sin和cos的多項式, 若已知初始傾角I0處的傾角函數(shù)及其偏導(dǎo)數(shù), 則可以在I0附近進行泰勒展開,以展開式計算傾角函數(shù)及其1階偏導(dǎo)數(shù), 如下式所示:

通常情況下, 分析法預(yù)報的弧段較短, 軌道平傾角變化幅度一般較小, 在3天弧段以內(nèi)基本不超過0.01?, 因此在分析法預(yù)報中能夠?qū)A角函數(shù)作泰勒展開, 以t0時刻的傾角函數(shù)（I0）及其1階、2階偏導(dǎo)數(shù)為初值, 在t時刻使用展開式計算傾角函數(shù). 傾角函數(shù)及其偏導(dǎo)數(shù)初值通過改進Gooding法計算, 由于只提供2階以內(nèi)的初值, 因此傾角函數(shù)展開式最高可以精確到2階,而偏導(dǎo)數(shù)只能展開至1階.

以泰勒展開式計算傾角函數(shù)可大幅提高計算速度,從而有效縮短分析法預(yù)報耗時.與改進Gooding法相比, 盡管使用網(wǎng)格Gooding法計算大批量傾角函數(shù)能在保持精度不變的同時提高計算速度, 但在軌道預(yù)報過程中傾角函數(shù)的計算只是其中一部分, 因此預(yù)報速度雖有所提升, 但并不明顯. 展開式進一步大幅減小了傾角函數(shù)的計算量, 經(jīng)測試在地球引力場取50階時,與網(wǎng)格Gooding法對比,1階、2階展開式的計算耗時均縮短了99%以上, 因此不難理解以展開式計算傾角函數(shù)能夠顯著提高分析法預(yù)報效率.

雖然泰勒展開方法可以有效提高計算效率, 但是計算精度卻受到一定影響. 由（8）式可知, 展開式存在截斷誤差, 傾角函數(shù)精度降低, 從而導(dǎo)致預(yù)報結(jié)果產(chǎn)生偏差. 只有當(dāng)截斷誤差較小、預(yù)報星歷偏差不足以影響預(yù)報精度時, 才可以將傾角函數(shù)展開式應(yīng)用到分析法軌道預(yù)報中.

3.2 展開式誤差分析

吳連大等[8–9]已經(jīng)驗證改進Gooding法的精度相比Gooding法更高一些, 而網(wǎng)格Gooding法計算結(jié)果與改進Gooding法一致, 因此可以將其作為計算傾角函數(shù)的精度對比標(biāo)準. 傾角函數(shù)展開式誤差按照下式計算:

其中ε1為傾角函數(shù)1階展開誤差,ε2為2階展開誤差,ε′為傾角函數(shù)偏導(dǎo)數(shù)1階展開誤差, 其中各項初值通過網(wǎng)格Gooding法計算得到.

3.2.1 算例

根據(jù)（9）式可知傾角函數(shù)展開式的誤差主要與初始傾角I0、傾角變化量?I以及3個傾角函數(shù)指標(biāo)有關(guān), 為了較為全面地分析傾角函數(shù)展開式的精度與這3種影響因素的關(guān)系, 我們設(shè)計了若干算例如表2所示. 其中地球引力場階次取為50, 指標(biāo)l取值范圍是2 ≤l≤50, 初始傾角I0從1?起每隔1?取值至179?, 此外也設(shè)置了0.01?及179.99?, 總共181組取值, 傾角變化量?I分別取為0.1?、0.01?、0.001?,傾角函數(shù)展開式精確至1階、2階, 分別使用泰勒展開方法和改進Gooding法計算傾角函數(shù), 通過這些算例分析其計算結(jié)果偏差. 在實際應(yīng)用場景中, 對于近圓軌道, 當(dāng)?shù)厍蛞鲭A次取50, 使用完整力模型時, 以分析方法預(yù)報3天后的平傾角變化量基本小于0.01?, 因此?I選擇0.01?為標(biāo)準. 傾角函數(shù)共有3個指標(biāo), 這里不再詳細討論m、p, 僅計算各階次l的傾角函數(shù)誤差最大值.

表2 傾角函數(shù)展開式誤差分析算例參數(shù)設(shè)置Table 2 Parameters setting for the error analysis test of the inclination function expansion

3.2.2 計算結(jié)果

由于傾角函數(shù)1階、2階展開計算速度差別很小, 而2階展開式精度更高, 因此主要分析不同初始傾角和傾角變化量對2階展開式誤差的影響, 并簡單對1階、2階展開式誤差進行比較. 2階展開式的誤差計算結(jié)果見圖1, 分別表示?I取不同角度時傾角函數(shù)2階展開式（左列）及其偏導(dǎo)數(shù)1階展開式誤差最大值（右列）的等高線圖, 其中誤差值取對數(shù), 橫軸為指標(biāo)l, 縱軸為初始傾角I0. 由圖1可知隨著l增加, 傾角函數(shù)2階展開誤差逐漸增大, 而根據(jù)等高線分布密度可知誤差增大速度逐漸降低. 另一方面在l <20時, 誤差大小幾乎不受初始傾角的影響,l≥ 20時, 誤差隨著初始傾角增加而減小, 當(dāng)超過90?之后又逐漸增大, 并且很明顯可以看出等高線關(guān)于I0=90?對稱. 傾角函數(shù)偏導(dǎo)數(shù)1階展開的誤差分布與上述規(guī)律基本一致, 只是誤差更大. 圖2為傾角函數(shù)1階展開誤差結(jié)果, 表示?I取為0.01?時展開式誤差最大值的等高線圖, 誤差取對數(shù), 橫軸為指標(biāo)l, 縱軸為初始傾角I0. 根據(jù)圖1 （d）與圖2結(jié)果,傾角函數(shù)偏導(dǎo)數(shù)的1階展開誤差約為傾角函數(shù)1階展開誤差的30倍.

圖2 傾角函數(shù)1階展開式的誤差最大值等高線圖Fig.2 Contour diagram for maximum error of 1st-order expansion of the inclination function

通過圖1中6個子圖的對比可知, 傾角函數(shù)及其偏導(dǎo)數(shù)誤差隨著?I增加而增大, 并且?I每增加10倍, 傾角函數(shù)誤差增大約1000倍, 而偏導(dǎo)數(shù)誤差增大約100倍, 這與展開式的截斷誤差是?I的n次冪相關(guān), 為便于理解, 將（9）式整理為如下形式:

據(jù)此可知傾角函數(shù)及其偏導(dǎo)數(shù)的展開式誤差分布和2階或3階偏導(dǎo)數(shù)分布是一致的.

根據(jù)計算結(jié)果顯示, 傾角函數(shù)展開誤差受指標(biāo)l和傾角變化量?I的影響較大,l越大誤差越大,而從（10）式可以直觀看出誤差與?I具有冪次關(guān)系.另外初始傾角對誤差影響較小, 尤其是l<20時, 隨著I0增加, 誤差幾乎沒有變化. 通過對比圖1 （c）和圖2可知傾角函數(shù)2階展開相比于1階展開誤差減小了約3個量級. 上述分析全部基于誤差最大值進行, 盡管該結(jié)果不能完全體現(xiàn)出傾角函數(shù)展開誤差在（l,m,k）空間的分布特征, 但由于在田諧項攝動分析解計算時傾角函數(shù)誤差最大值更能反映對計算精度的影響, 所以上述結(jié)果已經(jīng)足夠用來分析傾角函數(shù)展開式的適用場景了. 對于預(yù)報弧段較短、平傾角變化量較小的軌道, 當(dāng)?shù)厍蛞鋈≥^低的階次時, 在分析法預(yù)報中可以嘗試將傾角函數(shù)以2階展開式計算, 此時預(yù)報誤差相對較小, 不過實際應(yīng)用結(jié)果需要使用大量算例進行測試才能進一步分析.

4 分析法軌道預(yù)報誤差分析

分析法軌道預(yù)報中使用展開式計算傾角函數(shù)可有效提高預(yù)報速度, 但不可忽視展開式存在較大誤差, 而由（4）式可知在計算田諧項攝動分析解時誤差會在三重求和過程中逐漸累積, 并最終導(dǎo)致預(yù)報結(jié)果產(chǎn)生偏差. 雖然根據(jù)圖1與圖2可以初步分析傾角函數(shù)展開式的適用場景, 但為了進一步分析實際應(yīng)用中的預(yù)報精度, 需要使用大量算例作軌道預(yù)報測試. 一般情況下分析法預(yù)報誤差比較大,預(yù)報24 h的誤差為數(shù)十米至數(shù)百米不等, 而傾角函數(shù)展開導(dǎo)致的預(yù)報誤差并沒有這么大, 如果直接將展開后的預(yù)報結(jié)果與數(shù)值法預(yù)報結(jié)果比較, 則可能被分析法預(yù)報本身的誤差掩蓋. 由于展開式誤差不確定正負, 實際計算中可能出現(xiàn)展開式的預(yù)報精度比完整算法更高的“假象”, 因此以數(shù)值法預(yù)報星歷作為標(biāo)準無法準確評估展開式的預(yù)報精度. 因為整個預(yù)報過程中只有傾角函數(shù)的計算方法不同, 所以直接比較展開前后預(yù)報星歷的偏差即可. 本章以改進Gooding法為基準, 分別將傾角函數(shù)展開式精確到1階和2階, 對預(yù)報星歷偏差進行分析, 另一方面,為定量比較傾角函數(shù)展開前后分析法預(yù)報的速度,設(shè)計了若干算例, 并記錄預(yù)報耗時.

4.1 精度分析算例

傾角函數(shù)展開式誤差的影響因素為初始傾角I0、傾角變化量?I和階次l, 而階次l的取值范圍與地球引力場階次有關(guān), 因此選擇不同的初始軌道, 分析初始傾角I0和地球引力場階次L對預(yù)報星歷偏差的影響. 由于?I是在軌道預(yù)報過程中計算的, 無法通過給定初值分析, 所以不再考慮此因素.由于田諧項攝動分析解公式中存在因子（）l, 所以容易理解軌道半長徑a越大傾角函數(shù)展開誤差對預(yù)報結(jié)果的影響就越小, 因此主要選取低軌算例進行分析, 同時加入少量中高軌的算例, 算例的參數(shù)設(shè)置如表3所示. 分析法預(yù)報初始歷元按照簡約儒略日設(shè)為59000, 預(yù)報長度為3 d, 步長為30 min,L分別設(shè)5組取值, 地球參考半徑為6378.137 km, 軌道高度h通過將軌道半長徑a與地球參考半徑相減得到, 在300 km與40000 km之間取值, 共有46條低軌與14條中高軌, 軌道偏心率均設(shè)為0.001, 初始傾角則在90?以內(nèi)每隔5?選取了19組值, 軌道升交點經(jīng)度、近地點角和平近點角均設(shè)為50?.

表3 分析法預(yù)報精度分析算例參數(shù)設(shè)置Table 3 Parameters setting for test of calculation accuracy of analytical orbit prediction

分析法預(yù)報所用攝動模型為完整力模型, 包括帶諧項、田諧項、日月引力攝動、大氣阻力攝動和光壓攝動等. 一般情況下為了發(fā)揮分析法預(yù)報的速度優(yōu)勢, 步長設(shè)置為與預(yù)報時長相等, 而本次測試將步長設(shè)為30 min, 一方面是為了得到較密集的星歷, 便于分析星歷偏差, 另一方面可以更頻繁地更新參考星歷, 從而一定程度上提高預(yù)報精度,此外軌道預(yù)報是精密定軌過程中不可或缺的一部分, 而精密定軌中觀測數(shù)據(jù)的歷元間隔也決定了步長是較小的. 由于傾角函數(shù)展開誤差與偏心率無關(guān), 因此簡單地選取近圓軌道進行測試, 偏心率取0.001. 根據(jù)傾角函數(shù)展開式誤差分析結(jié)果可知誤差關(guān)于I0= 90?對稱, 因此不再考慮初始傾角大于90?的算例.

在計算田諧項攝動解時可能存在共振現(xiàn)象,即（4）式中某些項具有小分母問題（分母接近于0）,這主要取決于求和指標(biāo)l、m、p、q和衛(wèi)星平運動速率, 具體表達式可見文獻[10]. 一般情況下地球引力場階次越高, 滿足小分母關(guān)系的求和指標(biāo)越多,共振項也就越多, 另外低軌衛(wèi)星會更頻繁地遇見共振現(xiàn)象. 在本文所使用的分析法預(yù)報算法中, 計算田諧項時所有滿足小分母關(guān)系的求和指標(biāo)不再參與短周期項的計算.

測試中分別使用網(wǎng)格Gooding法、1階展開和2階展開這3種方法計算傾角函數(shù)及其偏導(dǎo)數(shù), 遍歷以上軌道預(yù)報算例, 得到3組預(yù)報星歷, 即各時刻的位置矢量:接著按照下式分別計算1階展開和2階展開的預(yù)報星歷偏差RMS:

其中Nt為星歷點數(shù)量, 根據(jù)RMS結(jié)果對傾角函數(shù)展開式在軌道預(yù)報中的適用性進行分析.

4.2 算例結(jié)果分析

為研究地球引力場階次和初始傾角與預(yù)報精度的關(guān)系, 挑選出相同軌道高度的算例, 對預(yù)報星歷偏差RMS作圖,并為了直觀分析,將RMS取對數(shù).另外挑選出若干初始傾角相同的算例, 比較不同軌道高度的預(yù)報星歷偏差RMS, 從而分析在不同軌道高度下, 傾角函數(shù)展開對預(yù)報結(jié)果的影響. 由于傾角函數(shù)2階展開精度更高, 且速度幾乎相同, 因此主要對2階展開的預(yù)報結(jié)果進行分析.

對于軌道高度在2000 km以內(nèi)的低軌算例, 一般情況下地球引力場階次越高, 預(yù)報星歷偏差RMS越大, 只是增加幅度很小.L≤50時, 預(yù)報星歷偏差RMS整體上隨著初始傾角的增加有逐漸減小的趨勢, 這和傾角函數(shù)誤差的變化規(guī)律相符, 而在I0=90?時RMS小幅增大.以500 km高度的軌道為例,傾角函數(shù)2階展開的預(yù)報星歷偏差如圖3所示, 橫坐標(biāo)為初始傾角I0, 縱坐標(biāo)為預(yù)報星歷偏差RMS, 并取對數(shù), 圖例表示地球引力場階次L. 由圖3可知地球引力場在50階以內(nèi)時RMS小于0.5 mm, 在I0=85?處RMS為極小值. 經(jīng)過進一步計算和分析發(fā)現(xiàn),由于各階次的引力場系數(shù)和傾角函數(shù)表達式不同,不同階次的RMS極小點存在差異, 例如L= 3時極小點在I0= 83.5?左右, 因為田諧項分析解為各階次計算結(jié)果的和, 所以圖3中極小點為綜合結(jié)果, 并不嚴格在90?處.

圖3 500 km高度軌道在取不同初始傾角時, 傾角函數(shù)2階展開的分析法預(yù)報星歷偏差Fig.3 The ephemeris deviations of analytical orbit prediction for 500 km altitude orbit that have different initial inclination with inclination function calculated by 2nd-order expansions

對于中高軌算例, 地球引力場階次對預(yù)報誤差影響很小, 而隨著初始傾角增加, 軌道高度不同時,預(yù)報星歷偏差RMS具有不同的變化規(guī)律, 但RMS整體很小. 比如軌道高度為10000 km算例的RMS隨初始傾角增加而逐漸增大, 而高度為36000 km的同步軌道算例的RMS卻有所浮動,如圖4所示,不過二者RMS均不超過10-7m.

另一方面為更清晰地分析預(yù)報誤差與軌道高度h的關(guān)系, 取L= 50, 將預(yù)報星歷偏差RMS關(guān)于h作圖. 軌道高度小于和大于2000 km, 結(jié)果分別如圖5、圖6所示, 其中橫坐標(biāo)為h, 縱坐標(biāo)為預(yù)報星歷偏差RMS的對數(shù)值, 圖例表示初始傾角I0. 可以發(fā)現(xiàn)當(dāng)軌道高度小于2000 km時, RMS隨軌道高度增加而逐漸減小, 但隨著初始傾角增加而減小, 且均小于1 mm. 而當(dāng)軌道高度大于2000 km時, RMS先隨軌道高度增加而減小, 然后逐漸增大并維持在10-8m至約10-7m之間.

圖5 不同高度的低軌傾角函數(shù)2階展開的分析法預(yù)報星歷偏差Fig.5 The ephemeris deviations of analytical orbit prediction for Low Earth Orbit at different altitude with inclination function calculated by 2nd-order expansions

圖6 不同高度的中高軌傾角函數(shù)2階展開的分析法預(yù)報星歷偏差Fig.6 The ephemeris deviations of analytical orbit prediction for Medium-High Earth Orbit at different altitude with inclination function calculated by 2nd-order expansions

通過上述算例結(jié)果可知,在分析法軌道預(yù)報中,當(dāng)?shù)厍蛞霾桓哂?0階時, 以2階泰勒展開式計算傾角函數(shù), 與改進Gooding法相比, 低軌衛(wèi)星的預(yù)報星歷偏差RMS小于1 mm, 而中高軌衛(wèi)星則更小, 以地球同步軌道衛(wèi)星為例, RMS約為10-7m.在實際應(yīng)用場景中, 分析法預(yù)報一般不需要使用高階次的地球引力場模型, 50階已經(jīng)足夠滿足精度要求, 由于分析法預(yù)報24 h的位置誤差一般大于10 m量級, 因此可以認為傾角函數(shù)2階展開不會影響分析法預(yù)報精度.

4.3 分析法預(yù)報速度分析

為定量比較傾角函數(shù)展開前后分析法預(yù)報的效率,設(shè)計若干算例進行分析,參數(shù)設(shè)置如表4所示.由于分析法預(yù)報耗時受軌道半長徑影響很小, 主要取決于地球引力場階次的大小以及預(yù)報節(jié)點的數(shù)量, 因此軌道半長徑不需考慮, 引力場階次分別取6個, 預(yù)報時長和預(yù)報步長分別設(shè)為3 d、30 min, 共有144個預(yù)報節(jié)點. 攝動模型依然為完整力模型, 軌道偏心率為0.001, 初始傾角在1?–90?之間間隔1?取值.

測試環(huán)境同樣如表1, 以Release模式編譯程序.分別使用改進Gooding法、網(wǎng)格Gooding法、1階展開和2階展開4種方法計算傾角函數(shù), 在L取不同值時, 分別對不同初始傾角的90組算例進行預(yù)報, 記錄預(yù)報耗時并求和,并以改進Gooding法為基準,計算其余3種方法的耗時減小比例. 統(tǒng)計結(jié)果如表5所示, 網(wǎng)格Gooding法使得分析法軌道預(yù)報速度有小幅提升, 在L較高時預(yù)報耗時減少13%左右, 而1階、2階展開對預(yù)報速度的提升效果十分顯著, 且二者差別很小, 隨著L的增加預(yù)報耗時減小比例逐漸增大,L= 50時達到48%左右. 因為分析法預(yù)報中除了傾角函數(shù)還有其他各種計算, 因此預(yù)報耗時的減小存在上限, 根據(jù)測試結(jié)果分析, 傾角函數(shù)的計算耗時約占預(yù)報總耗時的50%, 因此當(dāng)L增加時, 預(yù)報耗時最多可減少約50%.

5 結(jié)論

本文基于改進Gooding法調(diào)整了傾角函數(shù)的計算程序, 使其能夠一次性作大批量計算, 在計算50階以內(nèi)的傾角函數(shù)時, 相較于循環(huán)調(diào)用程序, 使用網(wǎng)格Gooding法耗時縮短了23.8%.

考慮到傾角函數(shù)計算公式為多項式, 可以進行泰勒展開, 而在短弧軌道預(yù)報中平傾角變化量一般很小, 所以本文對傾角函數(shù)展開式的計算精度進行分析, 發(fā)現(xiàn)2階展開式精度更高. 為了研究傾角函數(shù)展開式應(yīng)用到分析法軌道預(yù)報中的可行性, 設(shè)計了大量算例, 在預(yù)報過程中分別使用改進Gooding法和1階、2階展開式計算傾角函數(shù), 分析預(yù)報星歷偏差RMS. 根據(jù)測試結(jié)果, 在地球引力場取50階以內(nèi)時, 傾角函數(shù)2階展開后, 低軌衛(wèi)星的預(yù)報星歷偏差RMS小于1 mm, 而隨著軌道高度增加, RMS呈減小的趨勢. 一般在分析法預(yù)報中, 地球引力場階次取值較低, 基本不會超過50階, 因此將傾角函數(shù)作2階展開并不會影響分析法預(yù)報精度.

在計算耗時方面,與改進Gooding法相比,網(wǎng)格Gooding法的預(yù)報速度有小幅提升, 而1階、2階展開對預(yù)報速度的提升效果更加明顯, 且二者速度幾乎相同. 另一方面3種方法的預(yù)報耗時減小比例隨著地球引力場階次的增加而增大, 在50階時分別達到14%、49%、48%.