精密定軌中變分方程的數(shù)值解法及其檢驗(yàn)

洪 櫻,歐吉坤

(1.武漢紡織大學(xué)理學(xué)院,湖北武漢 430073;2.中國(guó)科學(xué)院測(cè)量與地球物理研究所,湖北武漢 430077)

精密定軌中變分方程的數(shù)值解法及其檢驗(yàn)

洪 櫻1,歐吉坤2

(1.武漢紡織大學(xué)理學(xué)院,湖北武漢 430073;2.中國(guó)科學(xué)院測(cè)量與地球物理研究所,湖北武漢 430077)

變分方程的求解是衛(wèi)星精密定軌中的關(guān)鍵步驟之一,通常都通過積分方法求解。首先從原理上詳細(xì)闡述一種新的數(shù)值方法求解變分方程,接著通過低軌衛(wèi)星精密定軌仿真軟件對(duì)其精度和效率進(jìn)行測(cè)試,且與常用數(shù)值積分方法的定軌結(jié)果進(jìn)行比較。結(jié)果表明,在實(shí)際軌道確定中,采用新的數(shù)值方法可達(dá)到較高的精度要求,且算法的效率更高。

變分方程;數(shù)值方法;精密定軌;積分方法

一、引 言

在精密軌道確定中,如果知道初始狀態(tài)矢量以及力矢量,則可通過積分衛(wèi)星的運(yùn)動(dòng)方程得到精密軌道。在求解初始狀態(tài)參數(shù)時(shí),需用目前狀態(tài)相對(duì)于初始狀態(tài)矢量的偏導(dǎo)數(shù),求解變分方程可得到構(gòu)成狀態(tài)轉(zhuǎn)移矩陣的這種偏導(dǎo)數(shù)。所以,在積分運(yùn)動(dòng)方程的同時(shí),必須求解變分方程[1-2]。對(duì)于變分方程的求解,許多學(xué)者都進(jìn)行了研究和探討,提出了許多有效的方法,如簡(jiǎn)化方法[3]、分析方法[3-4]、差分算法[5-6]、數(shù)值方法[2]等。其中數(shù)值方法是德國(guó)GFZ(GeoForschungs Zentrum)地學(xué)中心許國(guó)昌博士提出的一種新的一般化的求解變分方程的方法,其特點(diǎn)是通過差分代替微分使得變分方程轉(zhuǎn)化為線性代數(shù)方程組,進(jìn)而根據(jù)其初值條件推出嚴(yán)格的遞推解。雖然文獻(xiàn)[2]中給出了數(shù)值方法嚴(yán)格的推導(dǎo),但從未對(duì)該方法進(jìn)行過數(shù)值檢驗(yàn)和比較。使得該方法的應(yīng)用和推廣缺少必要的數(shù)值依據(jù)。

本文在此理論基礎(chǔ)之上,部分地應(yīng)理論原創(chuàng)者的要求,通過詳細(xì)的數(shù)值計(jì)算對(duì)這種新的求解變分方程的數(shù)值方法進(jìn)行了分析。首先從原理上詳細(xì)闡述了這種新的數(shù)值方法的求解過程,然后通過精密定軌仿真軟件對(duì)其進(jìn)行測(cè)試,且與常用的數(shù)值積分方法的結(jié)果進(jìn)行了比較,從而證明了新解法的可行性和有益性。

二、新的數(shù)值解法

在軌道確定和預(yù)報(bào)過程中,需解算運(yùn)動(dòng)方程[1-2]

式中,˙X(t)為衛(wèi)星的瞬時(shí)狀態(tài)向量;X(t0)為 t0時(shí)刻的初始狀態(tài)向量(記為 X0);F為時(shí)間 t、狀態(tài)向量 ˙X(t)以及動(dòng)力學(xué)參數(shù)向量 y的函數(shù),且

r和 ˙r分別為衛(wèi)星的位置和速度向量;f為衛(wèi)星所受力的向量和;m為衛(wèi)星質(zhì)量;Y是除位置、速度以外的動(dòng)力學(xué)參數(shù)向量。

在運(yùn)動(dòng)方程(1)兩邊對(duì) y求偏導(dǎo)數(shù)有

式中,＊號(hào)表示對(duì) F中顯含的參數(shù)矢量 y的偏導(dǎo),且記

式中,E為單位陣,0為零矩陣,下標(biāo)是其維數(shù)。則式(2)變?yōu)?/p>

其解構(gòu)成了狀態(tài)轉(zhuǎn)移矩陣,記為 Φ(t,t0)。則式(4)變?yōu)?/p>

式(5)稱為變分方程。

初值條件為

由于微分方程列與列之間是相互獨(dú)立的,因此只需考慮某列的方程的解。對(duì)于列 j,式(7)和初值條件式(8)可表示為

將式 (6)和式(3)代入式(5)得

若在時(shí)間間隔 [t0,t]中,步長(zhǎng) h=(t-t0)/m,則 tn=t0+nh,n=1,…,m,且

則式(9)變?yōu)?/p>

整理式(11)得

式中

則可得到 ψ(tn+1),速度向量 ψ˙(tn+1)可利用式(10)計(jì)算得到。

三、數(shù)值檢驗(yàn)算例

為了驗(yàn)證數(shù)值方法的可靠性,以CHAMP衛(wèi)星為例,在仿真過程中采用的動(dòng)力學(xué)模型包括：70×70階的 JG M3重力場(chǎng)模型、DT M78大氣阻力模型、Wahr潮汐模型、N體攝動(dòng)模型(行星位置由 JPL提供的行星歷表中精確地插值獲得)、光壓攝動(dòng)模型(假設(shè)衛(wèi)星帆板的法線方向始終與太陽垂直)以及廣義相對(duì)論攝動(dòng)模型[7-9]。分別采用積分方法和數(shù)值方法計(jì)算狀態(tài)轉(zhuǎn)移矩陣,歷元間隔取為 10 s,其中積分方法采用穩(wěn)定性較好、精度較高的 RKF7(8)法起步、10階的第 II種類型的 collocation-PECE算法[8]。由于低軌衛(wèi)星的偏心率較小,所以可采用固定步長(zhǎng)進(jìn)行計(jì)算[1]。取初始時(shí)刻為 2002年 1月 3日 0時(shí) 0分0秒,軌道初始狀態(tài)(慣性系下位置和速度)為x0=701 111.972 7 m, y0=-368 100.205 1 m, z0=6 723 471.518 4m, vx0=-4 695.158 9 m/s, vy0=6 017.559 8m/s, vz0=838.240 2m/s。為了避免比較數(shù)據(jù)過于龐大,這里僅列出狀態(tài)轉(zhuǎn)移矩陣的部分矩陣元,并且最少保留到數(shù)值方法與積分方法有差異的位數(shù)。將積分方法的計(jì)算結(jié)果記為Φ1,數(shù)值方法的計(jì)算結(jié)果記為Φ2,計(jì)算結(jié)果列于表 1。

表 1 積分方法和數(shù)值方法的計(jì)算結(jié)果比較

首先需要指出的是,隨著積分弧段的增加,積分方法所花的時(shí)間迅速增長(zhǎng)。在此算例計(jì)算中,積分方法和數(shù)值方法均積分了 1條軌道 (6個(gè)狀態(tài)方程),并且都需計(jì)算加速度的偏導(dǎo),而積分方法比數(shù)值方法多積分 6×6=36個(gè)變分方程。因此,積分方法共積分了 42個(gè)方程,而數(shù)值方法只積分了 6個(gè)方程。當(dāng)再增加一個(gè)動(dòng)力學(xué)參數(shù)時(shí),積分方法需積分6+7×7=55個(gè)方程,而數(shù)值方法仍僅需積分 6個(gè)方程,所以,待估動(dòng)力學(xué)參數(shù)越多,數(shù)值方法節(jié)省時(shí)間的優(yōu)越性就越明顯。當(dāng)然其前提條件是要保證計(jì)算的精度。

表 1列出了積分方法和數(shù)值方法計(jì)算的狀態(tài)轉(zhuǎn)移矩陣中的部分矩陣元。從表中可以看出,當(dāng)弧段較短時(shí),兩種方法計(jì)算出的狀態(tài)轉(zhuǎn)移矩陣中矩陣元的值相差不大,對(duì)于數(shù)值較大的矩陣元仍有 3～4位數(shù)字相同。但隨著積分弧段的增加,兩種方法的計(jì)算結(jié)果差異變大。

由于狀態(tài)轉(zhuǎn)移矩陣描述的是初始狀態(tài)和模型參數(shù)的變化對(duì)動(dòng)力學(xué)系統(tǒng)演化的影響,即變分方程的解——狀態(tài)轉(zhuǎn)移矩陣是無量綱的,無法直接描述其精度。而在求解初始狀態(tài)參數(shù)時(shí),需要用到狀態(tài)轉(zhuǎn)移矩陣,因此狀態(tài)轉(zhuǎn)移矩陣的精度會(huì)影響到最終的定軌結(jié)果。所以,不妨以精密軌道確定的結(jié)果來考察變分方程解法的效果。同樣采用上述仿真軟件并模擬 100min(略大于低軌衛(wèi)星的運(yùn)行周期)的P3碼偽距觀測(cè)值用于軌道確定,歷元間隔為 10 s。分別采用數(shù)值法和積分法求解變分方程,利用批處理算法得到的軌道初始狀態(tài)的改正量列于表 2中。

表 2 兩種方法求解變分方程得到的軌道初始狀態(tài)的改正量

表 2中方法A表示采用數(shù)值法求解變分方程,方法B表示采用積分法求解變分方程;OS0表示無誤差的 P3碼觀測(cè)值,OS1表示只包含接收機(jī)鐘差的P3碼觀測(cè)值,OS2表示包含接收機(jī)鐘差、GPS衛(wèi)星精密鐘差和觀測(cè)偶然誤差的 P3碼觀測(cè)值。

從表 2可以看出,在相同觀測(cè)條件下,用數(shù)值法和積分法求解變分方程得到的軌道初始狀態(tài)的改正量相差甚小,坐標(biāo)分量達(dá)到微米的數(shù)量級(jí),而速度分量達(dá)到 10-8～10-10m/s的數(shù)量級(jí)。

四、結(jié)束語

本文首次驗(yàn)證了變分方程的數(shù)值解算方法,從而確認(rèn)了定軌中一種全新的算法流程,即不再需要通過繁復(fù)的積分方法來求解變分方程。該方法并不像其他方法那樣受外部條件的限制,如軌道弧長(zhǎng)、動(dòng)力學(xué)模型的復(fù)雜程度等;且算法的推導(dǎo)嚴(yán)格,程序易于實(shí)現(xiàn)。從低軌衛(wèi)星軌道確定結(jié)果來看,在相同 GPS觀測(cè)條件下采用數(shù)值法求解變分方程的精度與積分法求解變分方程的精度相當(dāng)。由于積分方法的精度較高,可以推論這種新的數(shù)值方法的精度也是較高的。此外,從算法的效率來看,積分方法比數(shù)值方法多積分 n2(n為待估參數(shù)個(gè)數(shù))個(gè)變分方程。因此,這種新的數(shù)值方法在實(shí)際軌道計(jì)算中更有效。

[1] 李濟(jì)生.人造衛(wèi)星精密軌道確定[M].北京：解放軍出版社,1995：189-219.

[2] XU GC.GPS Theory,Algorithms andApplications[M]. Berlin：Springer-Verlag,2003：217-277.

[3] 劉林,張強(qiáng),廖新浩.人衛(wèi)精密定軌中的算法問題[J].中國(guó)科學(xué)：A輯,1998,28(9)：848-856.

[4] 張強(qiáng),劉林.軌道改進(jìn)中計(jì)算狀態(tài)轉(zhuǎn)移矩陣的分析方法[J].天文學(xué)報(bào),1999,40(2)：113-121.

[5] 張家祥,汪琦,楊捷興,等.彗木相撞預(yù)報(bào)[J].中國(guó)科學(xué)：A輯,1996,26(3)：280-288.

[6] 胡小工,黃輝,廖新浩.狀態(tài)轉(zhuǎn)移矩陣的差分算法及其應(yīng)用[J].天文學(xué)報(bào),2000,41(2)：113-122.

[7] 徐天河,楊元喜.基于能量守恒方法恢復(fù) CHAMP重力場(chǎng)模型[J].測(cè)繪學(xué)報(bào),2005,25(2)：75-81.

[8] MONTENBRUCK O,G ILL E.Satellite Orbits-Models, Methods andApplications[M].New York：Springer-Verlag,2000：53-154.

[9] BOCK H.EfficientMethods for Deter mining Precise Orbits of Low Earth Orbiters Using the Global Positioning System[D].Switzerland：Institute für Geod?sie und Photogrammetrie Eidg,2003：77-85.

Numerical Solution of Var iational Equation for Precise OrbitDeterm ination and Its Test

HONG Ying,OU Jikun

0494-0911(2010)12-0001-03

P228

B

2010-07-20

洪 櫻(1980—),女,湖北荊門人,碩士,講師,主要從事衛(wèi)星精密軌道確定與預(yù)報(bào)理論研究。