余 彪，范東明，游 為，谷延超，蘇 勇

(1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都611756；2.西南石油大學(xué)土木工程與建筑學(xué)院，成都610500)

GRACE重力反演中的軌道數(shù)值積分方法分析

余 彪1，范東明1，游 為1，谷延超1，蘇 勇2

(1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，成都611756；2.西南石油大學(xué)土木工程與建筑學(xué)院，成都610500)

首先,對比分析Runge-Kutta積分法、Adams-Cowell聯(lián)合并行積分法、Gauss-Jackson積分法以及天體力學(xué)中常用的外插積分法用于GRACE衛(wèi)星軌道積分及變分方程解算時的優(yōu)缺點,得出Gauss-Jackson算法優(yōu)于其它算法并推薦了相對較優(yōu)的階數(shù)和步長，從算法的推導(dǎo)過程中分析出了Gauss-Jackson算法優(yōu)于其它算法的原因。其次,對比不同積分方法抵抗誤差的能力，結(jié)果表明這些方法都不能有效抵抗誤差的干擾；并分析初始狀態(tài)向量誤差和攝動力誤差對軌道積分的影響，結(jié)果表明,衛(wèi)星初始狀態(tài)向量中的速度誤差(0.1mm/s)對軌道積分的影響大于位置誤差(10mm)對軌道積分的影響，攝動力中的隨機誤差對軌道積分的影響較大且無規(guī)律可循。針對Gauss-Jackson算法中的不足，提出基于移動窗口的多項式內(nèi)插算法和Gauss-Jackson算法相結(jié)合的組合方法(改進Gauss-Jackson算法)，通過模擬和實測數(shù)據(jù)的計算表明該方法在保證積分精度的前提下不僅提高了積分效率而且可以得到任意時刻的積分值。

軌道積分；積分精度；積分效率；誤差影響

0 引 言

在利用重力恢復(fù)與氣候?qū)嶒?Gravity recovery and climate experiment, GRACE)衛(wèi)星數(shù)據(jù)反演靜態(tài)或時變地球重力場模型以及下一代重力衛(wèi)星(如NGGM: Next generation gravity measurements和GRACE-FO:GRACE follow-on)反演重力場的模擬計算中，需要用數(shù)值積分方法計算衛(wèi)星參考軌道和動力學(xué)反演方法中的變分方程等。而高精度、高效率和高穩(wěn)定性的數(shù)值積分方法對于地球重力場模型的計算至關(guān)重要。

目前,廣泛使用的數(shù)值積分方法主要有Runge-Kutta積分法(簡稱RK方法)、Adams-Cowell聯(lián)合并行積分法(簡稱AC方法)、Gauss-Jackson積分法(簡稱GJ方法)、天體力學(xué)中常用的外插積分法(簡稱EX方法)和Collocation積分法(配置法)[1-9]。不同的積分方法優(yōu)缺點不同[1-2]，單步法(如RK方法)可以自起步積分，但是計算效率較低；多步法(如AC方法，GJ方法)雖然計算效率高于單步法，但是無法自起步；Collocation積分法一般用于受力相對簡單的高軌衛(wèi)星軌道積分，如國際上著名的GNSS定軌定位軟件Bernese就是采用這種方法，但低軌GRACE衛(wèi)星，受力比較復(fù)雜，實際計算時并不能很好地模型化衛(wèi)星受力情況，采用此種方法進行軌道積分時積分精度不如其它方法，所以本文未對其進行對比分析。即使同一種積分方法，采用不同的階數(shù)和步長進行計算時也會得到不同精度的結(jié)果。在對衛(wèi)星進行軌道積分時，衛(wèi)星初始狀態(tài)向量及攝動力中含有各種誤差，為了得到高精度和高穩(wěn)定性的積分結(jié)果需要分析誤差對軌道積分的影響。低軌衛(wèi)星軌道積分往往效率比較低，而且一般只能得到歷元為步長的整數(shù)倍處的軌道狀態(tài)，若采用GPS解算的幾何軌道作為標(biāo)準(zhǔn)軌道進行地球重力場模型反演，則可能需要將歷元為非步長整數(shù)倍處的軌道歸算到歷元為步長整數(shù)倍處的軌道。

因此本文首先從軌道積分精度、積分效率及狀態(tài)轉(zhuǎn)移矩陣積分精度方面對比分析現(xiàn)有的常用數(shù)值積分方法，然后分析了初始狀態(tài)向量和攝動力誤差對軌道積分的影響。針對目前軌道積分算法的不足，基于本文的研究提出了改進GJ算法，該算法采用適當(dāng)階數(shù)和較大步長的GJ算法，保證積分精度并結(jié)合基于移動窗口的多項式內(nèi)插算法計算出步長范圍內(nèi)任意歷元處的軌道。由于采用了較大的積分步長，使得積分效率相比于采用5s步長的積分效率成倍增加，而積分精度并不會降低；由于采用了基于移動窗口的多項式內(nèi)插算法可以內(nèi)插出任意時刻的積分軌道，降低了不規(guī)則歷元處軌道處理的難度。

1 數(shù)值積分方法

RK方法是最常用的單步積分法且存在多種不同算法[1]，即便是同一階方法也存在不同的算法且算法之間計算精度相差較大，本文采用軌道積分精度較高的8階RK算法[5-6,8-9]，8階RK積分方法每積分一步都需計算10次右函數(shù)，對于GRACE衛(wèi)星實際觀測數(shù)據(jù)，計算一次右函數(shù)(主要指衛(wèi)星所受的所有攝動力)大約需要耗時0.1s(在配有Intel Core四核處理器，主頻為2.50GHz的PC機上采用Fortran語言計算)，這對于實際軌道積分計算效率低下。

EX積分法是單步積分法[1,4]，其計算過程分為兩步，首先以不同步長計算某一歷元處的近似值然后對近似值外插逐漸逼近軌道真值。通過計算發(fā)現(xiàn),對于GRACE衛(wèi)星軌道，三角EX法[1]比菱形EX法[4]具有更高的數(shù)值穩(wěn)定性，所以本文采用三角EX算法。

AC方法為多步積分法，計算效率高于單步法，為了提高計算精度，本文在計算中結(jié)合預(yù)估校正迭代算法[1-3,5]，并與單步法的計算效率進行對比，本文在計算中采用8階RK方法提供初值。

GJ算法也為多步積分法，為了起到更好的對比效果，GJ方法起步值的計算并未采用文獻[3]和[10]中介紹的Mid-Corrector公式，而采用與AC方法相同的起步算法(8階RK方法)，計算中同樣結(jié)合預(yù)估校正迭代算法。

2 不同方法軌道和狀態(tài)轉(zhuǎn)移矩陣積分對比

通過比較不同方法的軌道和狀態(tài)轉(zhuǎn)移矩陣積分精度，得出最優(yōu)方法；比較同一方法不同階數(shù)和步長的計算結(jié)果得出計算精度較高的階數(shù)和步長范圍。目前GRACE重力反演是分弧段進行處理的，弧段長度一般定為一天，但衛(wèi)星繞地球運行時受地球引力的作用是連續(xù)累積的，從理論上說采用連續(xù)的弧段而不進行分弧段處理有益于地球重力場的反演；并且由于目前GRACE時變地球重力場模型的時間尺度一般為一個月，本文在積分時直接采用一個月作為積分的時間跨度。

為了得到軌道的解析解，考慮衛(wèi)星只受地球中心引力的二體問題。根據(jù)開普勒軌道理論計算衛(wèi)星軌道解析解，將不同積分方法以不同的階數(shù)和步長計算積分軌道并與開普勒軌道比較得到積分誤差。采用表1中的GRACE衛(wèi)星軌道參數(shù)作為初始?xì)v元的開普勒軌道根數(shù)，并以不同步長計算得到一個月的開普勒軌道。

表1 GRACE衛(wèi)星初始軌道根數(shù)

2.1 衛(wèi)星軌道積分誤差對比

對于EX方法，外插計算一次算法階數(shù)增加2[1]，所以該方法只有偶數(shù)階。當(dāng)選取的階數(shù)很高時，積分軌道仍收斂，但是該算法的積分效率會隨著階數(shù)的增加迅速降低(第2.2節(jié)詳細(xì)說明)，考慮到實際計算效率本文只分析前20階的計算精度。圖1(b)表示不同階數(shù)和步長的EX方法積分軌道誤差。

對于AC方法，滿足積分精度要求的階數(shù)為8階到13階。圖1(c)表示不同階數(shù)和步長的AC方法積分軌道誤差。從圖1(c)可以看出，該方法在步長小于15s時，采用任何一種階數(shù)進行積分都無法滿足mm級的誤差要求(如常用的12階5s步長積分一個月的最大誤差為40.4mm)。

對于GJ方法，滿足積分精度要求的階數(shù)為3階到12階。圖1(d)表示不同階數(shù)和步長的GJ方法積分軌道誤差。從圖1(d)可以看出，GJ方法積分誤差普遍比上述幾種方法小(可達兩個量級)，而且該方法3到12階5s步長的積分誤差皆滿足要求。

對比分析圖1的結(jié)果可以得出以下兩點結(jié)論：1)對于衛(wèi)星軌道積分并不是積分步長越短越好、積分階數(shù)越高越好。其可能的原因為，當(dāng)積分步長較小時，積分計算的歷元個數(shù)增多，則軌道積分累計誤差或計算機的累計舍入誤差會增大；當(dāng)積分步長較大時，用于內(nèi)插右函數(shù)值的弧段較長，則內(nèi)插出的右函數(shù)值與實際右函數(shù)值的差值較大[11]，從而影響軌道積分的精度。當(dāng)積分的階數(shù)較大時，積分下一個歷元時所用到的前面歷元的個數(shù)增多，則用于內(nèi)插右函數(shù)的弧段增長，內(nèi)插出的右函數(shù)誤差可能增大[11]；而當(dāng)階數(shù)較小時，則用于內(nèi)插右函數(shù)的多項式的階數(shù)小，右函數(shù)的插值余項變大[11]，即右函數(shù)的近似誤差大。所以積分步長和積分階數(shù)應(yīng)根據(jù)精度、效率和觀測數(shù)據(jù)采樣要求合理選擇。2)從不同方法軌道積分誤差角度看，GJ方法優(yōu)于其它方法。通過對比分析各積分算法的推導(dǎo)過程，以GJ算法和AC算法為例來說明GJ算法優(yōu)于其它算法的原因。GJ算法與AC算法同為多步法，而且公式推導(dǎo)過程相似。兩套算法都是基于基本數(shù)學(xué)物理公式：

(1)

(2)

之前可認(rèn)為都是一樣的。AC算法將式(2)化為：

(3)

GJ算法將式(2)轉(zhuǎn)化為：

(rn+1-rn)-(rn-rn-1)=▽2rn+1=

(4)

2.2 衛(wèi)星軌道積分效率對比

在衛(wèi)星軌道積分時，除了要考慮積分誤差，還應(yīng)該考慮積分效率，影響軌道積分效率的主要因素是右函數(shù)的計算次數(shù)。一般情況下，單步法的計算效率低于多步法，但為了提高多步法的計算精度，常常結(jié)合預(yù)估校正迭代算法，本節(jié)主要比較結(jié)合預(yù)估校正迭代算法的多步法和單步法的計算效率。

對于8階RK方法，積分一步大約需要花費1s的時間。圖2(a)反應(yīng)了不同步長8階RK方法計算一個月的積分軌道時右函數(shù)計算次數(shù)。

對于EX方法，該方法需要以不同的步長計算某一歷元處衛(wèi)星的狀態(tài)得到EX方法的初值，在計算初值的過程中需要多次計算右函數(shù)。圖2(b)反應(yīng)的是不同步長和階數(shù)的EX方法計算一個月的積分軌道時右函數(shù)計算次數(shù)。從圖2(b)可以看出，EX方法右函數(shù)的計算次數(shù)會隨著階數(shù)的增高成倍的增加，而且相比于同階的RK方法，EX方法右函數(shù)的計算次數(shù)明顯要更多。

對于AC聯(lián)合并行預(yù)估校正迭代算法，雖然需要進行迭代計算，但在對GRACE衛(wèi)星軌道的積分計算中發(fā)現(xiàn)不需要迭代幾次即可滿足精度要求，所以并不會增加多少額外的計算時間。圖2(c)反應(yīng)的是不同步長的偶數(shù)階AC方法計算一個月的積分軌道時右函數(shù)的計算次數(shù)。通過與圖2(a)和圖2(b)對比發(fā)現(xiàn),該方法的右函數(shù)計算次數(shù)明顯少于RK方法和EX方法。

對于GJ方法，采用同樣的積分策略得到圖2(d)。從圖中可以看出，GJ方法也不需要迭代幾次就可滿足精度要求，且階數(shù)越高計算效率相對越高。對比圖2(c)與圖2(d)可知，AC方法與GJ方法的計算次數(shù)幾乎相等，且都小于RK方法和EX方法。

通過分析比較不同方法用于GRACE衛(wèi)星軌道積分時右函數(shù)計算次數(shù)，可以得出EX方法計算時右函數(shù)計算次數(shù)最多，所以積分效率最低，其次為RK方法，AC方法和GJ方法。

2.3 狀態(tài)轉(zhuǎn)移矩陣積分誤差對比

在用動力學(xué)法反演地球重力場模型時需要計算狀態(tài)轉(zhuǎn)移矩陣和參數(shù)敏感矩陣，一般通過對變分方程進行數(shù)值積分實現(xiàn)[1,5-6]。積分過程中，首先要積分每個歷元處衛(wèi)星的狀態(tài)向量，再由衛(wèi)星狀態(tài)向量計算狀態(tài)轉(zhuǎn)移矩陣或參數(shù)敏感矩陣的導(dǎo)數(shù)(右函數(shù))，最后由右函數(shù)計算下一歷元的狀態(tài)轉(zhuǎn)移矩陣或參數(shù)敏感矩陣。對于二體問題，可通過解析法得到狀態(tài)轉(zhuǎn)移矩陣，文獻[12]提供了二體問題狀態(tài)轉(zhuǎn)移矩陣解析解算法。本文在二體問題下比較不同方法用于狀態(tài)轉(zhuǎn)移矩陣積分時的積分誤差，仍然采用表1中的GRACE衛(wèi)星軌道參數(shù)作為初始?xì)v元的開普勒軌道根數(shù)。由于在積分狀態(tài)轉(zhuǎn)移矩陣時需要積衛(wèi)星軌道，EX方法積分衛(wèi)星軌道效率低下，所以此處并未對EX方法進行比較。狀態(tài)轉(zhuǎn)移矩陣積分效率主要受軌道積分的影響，所以軌道積分效率高的方法狀態(tài)轉(zhuǎn)移矩陣積分同樣效率高。表2為不同方法積分一個月的狀態(tài)轉(zhuǎn)移矩陣誤差(狀態(tài)轉(zhuǎn)移矩陣36個元素中積分誤差最大值)，表中符號*表示積分發(fā)散。從表2可以看出,GJ方法積分狀態(tài)轉(zhuǎn)移矩陣誤差比AC方法小一個量級，與RK積分方法的積分誤差相當(dāng)。

表2 不同參數(shù)下的各種方法狀態(tài)轉(zhuǎn)移矩陣積分誤差最大值

基于第2.1～2.3節(jié)的分析考慮到GJ算法用于GRACE衛(wèi)星軌道和狀態(tài)轉(zhuǎn)移矩陣積分相對于其它三種積分方法有較大優(yōu)勢，表3給出了GJ算法積分誤差相對較小的階數(shù)以及對應(yīng)的步長，其它三種積分方法不推薦用于長弧段的GRACE衛(wèi)星軌道及狀態(tài)轉(zhuǎn)移矩陣積分。

表3 GJ算法相對較優(yōu)階數(shù)及對應(yīng)步長

3 誤差對軌道積分影響分析

由于軌道積分需要初始時刻的衛(wèi)星狀態(tài)向量，而從GRACELevel1B的GNV1B文件或其它軌道文件中獲得的數(shù)據(jù)含有誤差，該誤差會影響積分軌道的計算。除了初始時刻的衛(wèi)星狀態(tài)含有誤差之外，積分計算時攝動力中含有的誤差也會對積分軌道產(chǎn)生影響。

3.1 初始狀態(tài)誤差對軌道積分的影響

為了分析衛(wèi)星初始時刻狀態(tài)向量誤差對軌道積分的影響，在軌道積分時分別僅對初始狀態(tài)向量中的位置和速度向量加入不同量級的誤差，分析該誤差對軌道積分的影響。通過采用不同的積分方法積分后發(fā)現(xiàn)初始狀態(tài)誤差對不同積分方法的影響幾乎相同，因此本文以GJ方法為例分析初始狀態(tài)誤差對軌道積分的影響。圖3(a)為僅對初始時刻衛(wèi)星位置向量加入不同量級誤差(100mm級誤差表示X、Y、Z方向各加100mm誤差，其它情形類似)后積分軌道誤差圖。圖3(b)為僅對衛(wèi)星初始時刻速度向量加入不同量級誤差(mm/s級誤差表示X、Y、Z方向各加1mm/s誤差，其它情形類似)后積分軌道誤差圖。為了分析初始狀態(tài)中位置和速度向量誤差對積分軌道的綜合影響，在初始狀態(tài)中共同加入不同量級的位置和速度誤差后得到積分軌道誤差圖3(c)。

通過比較分析圖3的計算結(jié)果可以得出以下結(jié)論：初始狀態(tài)誤差對積分軌道的影響為隨著積分時間增加積分軌道誤差呈現(xiàn)線性累加特性。衛(wèi)星初始狀態(tài)誤差對軌道積分的影響遠(yuǎn)遠(yuǎn)大于積分算法本身的誤差。對比圖3(a)、3(b)和3(c)發(fā)現(xiàn)，衛(wèi)星初始位置和速度向量中速度向量誤差(0.1mm/s)對軌道積分的影響大于位置向量誤差(10mm)對軌道積分的影響，若位置誤差為10mm量級，則當(dāng)速度誤差為0.01mm/s量級時，兩者對軌道積分的影響才能保持一致，而在當(dāng)前衛(wèi)星定軌精度下是很難達到的。其可能的原因為，軌道積分時速度為位置的右函數(shù)，而速度的右函數(shù)為衛(wèi)星所受的力。位置誤差所引起的攝動力變化相對較小，從而引起的速度變化也相對較小，所以在隨后的積分中位置誤差相對較?。欢绻跏紶顟B(tài)中速度存在誤差，由于速度為位置的右函數(shù)則在對位置進行數(shù)值積分時受該速度誤差影響較大，會使位置產(chǎn)生相對較大的偏差，具有較大偏差的位置對應(yīng)的攝動力偏差也相對較大，則在隨后的積分中位置誤差較大。

3.2 攝動力中的隨機誤差對軌道積分的影響

在對衛(wèi)星進行軌道積分時，由于先驗力模型及非保守力測量數(shù)據(jù)的不準(zhǔn)確會對軌道積分產(chǎn)生較大的影響。非保守力數(shù)據(jù)的誤差大致為10-10ms-2量級，保守力中主要是先驗重力場模型的誤差較大，雖然在重力反演中會對先驗重力場模型和非保守力數(shù)據(jù)進行校正，但如果誤差較大會影響參數(shù)的估計，本文取兩個不同的先驗重力場模型并取至不同的階次得到先驗重力場模型誤差引起的衛(wèi)星引力誤差量級。通過分別計算GGM05C模型(取至140階次)和EGM2008模型(取至70階次)對衛(wèi)星的引力并對引力的差做統(tǒng)計(時間為1d，間隔為5s)得到圖4(a)。根據(jù)圖4(a)的結(jié)果在二體問題下，對中心引力加上一個服從正態(tài)分布N(0, (5×10-7ms-2)2)的隨機誤差，圖4(b)為以GJ方法計算100次攝動力誤差對軌道積分結(jié)果影響圖。從圖4(b)可以看出，攝動力誤差對軌道積分的影響可超過20000m,大于初始狀態(tài)誤差對軌道積分的影響且無規(guī)律可循，影響也是隨機的?；诒疚牡姆治龅贸鋈糁苯硬捎?個月作為積分弧段進行重力反演，則必須考慮如何減小該誤差的影響。

4 改進Gauss-Jackson方法

雖然GJ算法具有較高的計算精度且有多種階數(shù)和步長可以選擇，但是該算法只能計算出歷元為步長整數(shù)倍處的積分軌道。由于GRACE LEVEL1B數(shù)據(jù)中給出的是5s間隔的標(biāo)準(zhǔn)軌道和星間距觀測數(shù)據(jù)，所以一般采用5s步長進行軌道積分，對于長弧段軌道積分積分效率較低?；谇拔牡姆治龅贸霾捎幂^大的步長進行積分也可以達到較高的積分精度，而且提高了積分效率。為了得到步長內(nèi)任意點的積分軌道可采用適當(dāng)?shù)膬?nèi)插算法[13]，只要步長的大小控制在一定的范圍內(nèi)，并選擇合適的內(nèi)插算法，內(nèi)插的精度可以很高。內(nèi)插計算不僅可以得到5s步長整數(shù)倍處的軌道值，而且可以得到任意歷元處的軌道值，當(dāng)以GPS解算的幾何軌道作為標(biāo)準(zhǔn)軌道時可簡化幾何軌道的預(yù)處理。

本文選擇8階30s作為積分階數(shù)和步長，在內(nèi)插算法上，采用10點9階基于移動窗口的多項式內(nèi)插算法[6,14]，在起始段和末尾段采用同階拉格朗日算法[11]，除了弧段兩端外，其它區(qū)域能保證為內(nèi)插窗口的中心區(qū)域。通過模擬數(shù)據(jù)計算得到該種內(nèi)插算法內(nèi)插出的軌道誤差在10-8～10-9m量級，起始段和末尾段內(nèi)插出的軌道誤差為10-6～10-7m量級，狀態(tài)轉(zhuǎn)移矩陣誤差在10-11量級，完全能滿足軌道積分和狀態(tài)轉(zhuǎn)移矩陣積分要求。內(nèi)插算法的計算時間相對于右函數(shù)的計算時間幾乎可以忽略不計，所以改進GJ算法不僅保證了軌道積分精度，而且提高了軌道積分效率。圖5為二體問題下采用改進GJ方法計算的積分軌道誤差圖，完全能滿足GRACE衛(wèi)星積分精度要求。

為了分析改進GJ方法在實際衛(wèi)星軌道積分中的性能，本文以GRACE Level1B數(shù)據(jù)以及相應(yīng)的力模型(見表4)進行軌道積分試驗。計算時分別采用8階5s步長的GJ方法和改進GJ方法進行積分，并從積分精度和積分效率方面進行比較。由于真實數(shù)據(jù)中存在數(shù)據(jù)間斷和粗差，需要對間斷數(shù)據(jù)進行填充并剔除粗差。為了計算方便，此處利用一天的實際數(shù)據(jù)進行軌道積分(未對初始狀態(tài)向量和力模型參數(shù)校正)并與JPL提供的軌道數(shù)據(jù)(GRACE GNV1B)對比。

對比發(fā)現(xiàn)，改進GJ方法的最大積分誤差比8階5s步長GJ方法略小(200mm左右)，并在計算效率上有了較大提高。采用8階5s步長GJ方法用時約57min，采用改進GJ方法用時約18min，積分效率提高了將近3倍。

表4 軌道積分計算所用的力模型

5 結(jié) 論

本文對比了多種軌道積分方法并分析了積分起始數(shù)據(jù)誤差對軌道積分結(jié)果的影響，提出了改進GJ算法?；趯Ρ确治龅贸隽艘韵掠幸娴慕Y(jié)論：

1)對于GRACE衛(wèi)星長弧段軌道積分，GJ算法由于在積分的過程中沒有直接采用前面歷元的積分軌道值使得其積分誤差比RK、EX和AC算法小(可達2個量級)。

2)EX算法在積分階數(shù)很高的情況下仍可以保持較高的積分精度，但是對于低軌衛(wèi)星計算效率低下；AC算法和GJ算法同為多步法，但對于GRACE衛(wèi)星軌道積分AC算法不及GJ算法，不推薦在長弧段軌道積分中使用EX算法和AC算法，RK算法只建議用于多步法的起步值計算。

3)對于GRACE衛(wèi)星軌道和狀態(tài)轉(zhuǎn)移矩陣積分，GJ算法較適宜的積分階數(shù)為8～10階，對應(yīng)的步長為5～40s和5～50s。

4)衛(wèi)星初始狀態(tài)誤差和力模型誤差對軌道積分的影響遠(yuǎn)大于算法本身的誤差；初始狀態(tài)中的速度誤差(0.1mm/s)對軌道積分的影響比位置誤差(10mm)大將近1個量級，且該類誤差對軌道積分的影響呈現(xiàn)線性累加特性；力模型誤差對軌道積分的影響較大且無規(guī)律可循。

5)GJ算法中取合適的較大積分步長(30s)并結(jié)合適宜階數(shù)(9階)的基于移動窗口的多項式內(nèi)插算法不僅可以大大提高積分效率而且可以得到任意時刻衛(wèi)星的積分軌道和狀態(tài)轉(zhuǎn)移矩陣，有利于GRACE衛(wèi)星不規(guī)則歷元的幾何軌道和星間觀測值聯(lián)合進行衛(wèi)星重力反演。

[1] Montenbruck O, Gill E. Satellite orbits: models, methods, and applications[M]. Berlin: Springer-Verlag, 2005.

[2] Beutler G. Methods of celestial mechanics[M]. Berlin: Springer-Verlag, 2005.

[3] Berry M M, Healy L M. Implementation of Gauss-Jackson integration for orbit propagation [J]. The Journal of the Astronautical Sciences, 2004, 52(3): 331-357.

[4] Bulirsch R, Stoer J. Numerical treatment of ordinary differential equations by extrapolation methods [J].Numerische Mathematik, 1966, 8(1): 1-13.

[5] 王正濤. 衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論與方法[D]. 武漢: 武漢大學(xué), 2005.[Wang Zheng-tao. Theory and methodology of earth gravity field recovery by satellite-to-satellite tracking data[D]. Wuhan: Wuhan University, 2005.]

[6] 游為. 應(yīng)用低軌衛(wèi)星數(shù)據(jù)反演地球重力場模型的理論和方法[D]. 成都: 西南交通大學(xué), 2011. [You Wei. Theory and methodology of earth′s gravitational field model recovery by LEO data[D]. Chengdu: Southwest Jiaotong University, 2011.]

[7] 龔曉穎, 蔡詩響, 李征航. 小步長低階Collocation方法軌道積分性能研究[J]. 宇航學(xué)報,2013, 34(5): 642-648. [Gong Xiao-ying, Cai Shi-xiang, Li Zheng-hang. Performance research of small-step length and low-order collocation method in orbit determination[J]. Journal of Astronautics, 2013, 34(5): 642-648.]

[8] Xu G C. GPS theory, algorithms and applications[M]. Berlin: Springer-Verlag, 2007.

[9] 劉大杰, 施一民, 過靜珺.全球定位系統(tǒng)(GPS)的理論與數(shù)據(jù)處理[M]. 上海: 同濟大學(xué)出版社, 1997: 257-263.

[10] 羅志才,周浩, 鐘波, 等. Gauss-Jackson積分算法分析與驗證[J]. 武漢大學(xué)學(xué)報·信息科學(xué)版, 2013, 38(11)：1364-1368.[Luo Zhi-cai, Zhou Hao, Zhong Bo, et al. Analysis and validation of Gauss-Jackson integral algorithm[J]. Geomatics and Information Science of Wuhan University, 2013, 38(11): 1364-1368.]

[11] 李大美, 李素貞, 朱方生. 計算方法[M]. 武漢: 武漢大學(xué)出版社, 2012: 94-130.

[12] 蔣虎, 黃珹, 胡小工. 二體問題狀態(tài)轉(zhuǎn)移矩陣的分析解和數(shù)值解比較[J]. 云南天文臺臺刊,1999, 4: 19-29. [Jiang Hu, Huang Cheng, Hu Xiao-gong. Comparison between analytic and numerical solutions of state transition matrices of two body problem[J]. Publications of the Yunnan Astronomical Observatory, 1999, 4: 19-29 .]

[13] 鄭偉, 許厚澤, 鐘敏, 等. 不同插值法對下一代衛(wèi)星重力反演精度的影響[J]. 宇航學(xué)報, 2014, 35(3): 269-276. [Zheng Wei, Xu Hou-ze, Zhong Min, et al. Influences of different interpolation methods on accuracy of next-generation satellite gravity recovery[J]. Journal of Astronautics, 2014, 35(3): 269-276.]

[14] 游為, 沈云中, 范東明, 等. 基于衛(wèi)星軌道擾動理論的重力反演方法[J]. 地球物理學(xué)報, 2010, 53(11): 2574-2581. [You Wei, Shen Yun-zhong, Fan Dong-ming, et al. The algorithm of earth′s gravitational field recovery based on satellite′s orbital perturbation theory[J]. Chinese Journal of Geophysics, 2010, 53(11): 2574-2581.]

通信地址：四川省成都市郫都區(qū)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院(611756)

電話：18382157390

E-mail:yubiao0221@qq.com

游 為(1985-)，男，講師，博士，主要從事衛(wèi)星重力測量及測量數(shù)據(jù)處理研究。本文通信作者。

通信地址：四川省成都市郫都區(qū)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院(611756)

電話：18123216821

E-mail: youwei1985@foxmail.com

(編輯：牛苗苗)

Analysis of Orbit Numerical Integration Methods in Earth′s Gravitational Field Recovery by GRACE

YU Biao1, FAN Dong-ming1, YOU Wei1, GU Yan-chao1, SU Yong2

(1.Faculty of Geosciences and Environmental Engineering, Southwest Jiao-tong University, Chengdu 611756, China; 2.School of Civil Engineering and Architecture, Southwest Petroleum University, Chengdu 610500, China)

This paper assesses the existing numerical integration methods including the Runge-Kutta methods, the Adams-Cowell methods, the Gauss-Jackson methods and the extrapolation methods in the computation of the GRACE satellite orbits and state transition matrix. According to the results, we recommend the Gauss-Jackson methods with optimal parameters to integral orbit and also analyze the advantages of the Gauss-Jackson methods relative to the other methods. Then the resistance capacity of the random errors which are included in the satellite initial state for those methods is analyzed, and the numerical results indicate that all methods show the similar little resistance capacity. The integral orbits are more sensitive to the satellite initial state velocity errors (0.1mm/s) than the initial state position errors(10mm). Compared to the errors in the satellite initial state, the perturbation force model errors have a significant impact on the orbit integral precision. Finally, we proposes a modified method which combines the Gauss-Jackson algorithm and the moving-window polynomial interpolation algorithm to overcome the shortage of the large step size of the Gauss-Jackson method. The simulated and actual results show that the modified method can give the satellite position and state transition matrix at any time and has a significant improvement in computing efficiency while possessing high accuracy.

Orbit integration; Integral precision; Integral efficiency; Error influence

2016-10-24;

2017-01-03

國家自然科學(xué)基金(41574018,41404018)

P228.1

A

1000-1328(2017)03-0253-09

10.3873/j.issn.1000-1328.2017.03.005

余 彪(1990-)，男，碩士生，主要從事衛(wèi)星重力測量研究。