汪宏波 趙長印 柳仲貴 張偉

(1中國科學(xué)院紫金山天文臺(tái) 南京 210008) (2中國科學(xué)院空間目標(biāo)與碎片觀測重點(diǎn)實(shí)驗(yàn)室 南京 210008) (3宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室 西安 710043) (4北京跟蹤與通信技術(shù)研究所 北京 100094)

基于誤差發(fā)散規(guī)律的低軌衛(wèi)星大氣阻力系數(shù)計(jì)算方法?

汪宏波1,2,3?趙長印1,2柳仲貴4張偉1,2

(1中國科學(xué)院紫金山天文臺(tái) 南京 210008) (2中國科學(xué)院空間目標(biāo)與碎片觀測重點(diǎn)實(shí)驗(yàn)室 南京 210008) (3宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室 西安 710043) (4北京跟蹤與通信技術(shù)研究所 北京 100094)

低軌衛(wèi)星軌道預(yù)報(bào)精度受到大氣模型和大氣阻力系數(shù)精度的制約,給一些高精度的空間和航天任務(wù)帶來不利影響.提出了一種基于沿跡方向誤差發(fā)散規(guī)律的大氣阻力系數(shù)計(jì)算新方法.首先通過理論推導(dǎo)給出低軌衛(wèi)星軌道預(yù)報(bào)中沿跡誤差發(fā)散的分析表達(dá)式,定量描述初值誤差和模型誤差對(duì)沿跡誤差的綜合影響;提出利用定軌段的基本信息,優(yōu)選預(yù)報(bào)段所采用的阻力系數(shù),抑制沿跡誤差的發(fā)散速率,從而降低沿跡方向預(yù)報(bào)誤差的最大值,提高短期預(yù)報(bào)精度.以400 km附近的GRACE-A衛(wèi)星的全弧段星載GPS高精度資料為基礎(chǔ),檢驗(yàn)了方法的精度和成功率.結(jié)果表明:相對(duì)于傳統(tǒng)的定軌預(yù)報(bào)方法,新方法能提高24 h短期預(yù)報(bào)精度約45%,成功率約71%,總體有效率約86%;方法對(duì)低、中、高等3種太陽輻射水平均有效,對(duì)于中低等級(jí)的地磁擾動(dòng)也有效,具備較好的應(yīng)用價(jià)值.

天體力學(xué):軌道計(jì)算和定軌,方法：數(shù)值,高層大氣阻力

1 引言

一直以來,低軌衛(wèi)星大氣阻力攝動(dòng)加速度的計(jì)算存在較大誤差,是軌道預(yù)報(bào)誤差的最主要來源.大氣阻力難以精確計(jì)算的原因有兩個(gè)方面:一是大氣密度模型普遍不準(zhǔn),目前還沒有哪個(gè)模型能夠準(zhǔn)確描述和預(yù)測高層大氣在不同時(shí)空、不同環(huán)境下的大氣密度[1];二是大氣阻力系數(shù)存在較大偏差,從物理的角度看,它與大氣分子的成分和濃度、衛(wèi)星表面材料及溫度、撞擊角度和速度等因素有關(guān),很難被精確測定.

計(jì)算大氣阻力系數(shù)一般有數(shù)值模擬和軌道擬合兩種方式,其中數(shù)值模擬是基于大氣分子與衛(wèi)星表面撞擊作用的理論模型,輸入各種條件參數(shù)進(jìn)行大量的蒙特卡羅模擬,從而得到阻力系數(shù)[2].這種方法比較接近物理本質(zhì),結(jié)果也相對(duì)精確,但缺點(diǎn)是在實(shí)際任務(wù)中很多輸入條件是不確定的,而且計(jì)算耗時(shí)巨大,適用性和時(shí)效性較差.

通過軌道擬合方法計(jì)算大氣阻力系數(shù),是指在定軌階段將大氣阻力系數(shù)CD作為待估參數(shù),與衛(wèi)星狀態(tài)參數(shù)(→r,˙→r)共同求解,進(jìn)而應(yīng)用在軌道預(yù)報(bào)中.這種方法原理簡單,計(jì)算方便,能夠取得穩(wěn)健的預(yù)報(bào)精度,在實(shí)際任務(wù)中得到了普遍應(yīng)用.但是對(duì)于部分低軌衛(wèi)星的高精度空間任務(wù),例如:空間交會(huì)、碰撞預(yù)警、精密跟蹤等,目前的預(yù)報(bào)精度還不能滿足要求.這主要是因?yàn)槎ㄜ壎诬壍罃M合中,阻力系數(shù)與衛(wèi)星軌道狀態(tài)值、大氣模型形成了自洽系統(tǒng),阻力系數(shù)CD是定軌階段的最優(yōu)估計(jì),它吸收了觀測資料、大氣模型、衛(wèi)星面質(zhì)比等誤差,已經(jīng)成為一個(gè)統(tǒng)計(jì)意義上的綜合參數(shù),偏離了原先的物理本質(zhì).因此在預(yù)報(bào)階段系數(shù)CD的適用性就會(huì)受到挑戰(zhàn),尤其是在長期預(yù)報(bào)中準(zhǔn)確性會(huì)大大下降.

我們認(rèn)為,在大氣模型精度短時(shí)期難以提高的背景下,大氣阻力系數(shù)的計(jì)算方法是有可能突破的技術(shù)問題.本文將在大氣模型確定的情況下,尋找一種新的方法計(jì)算大氣阻力系數(shù),有效降低低軌衛(wèi)星24 h的軌道預(yù)報(bào)誤差,滿足高精度短期預(yù)報(bào)的需求.第2節(jié)將給出方法的原理,第3節(jié)給出實(shí)驗(yàn)數(shù)據(jù)的算例和統(tǒng)計(jì)結(jié)果,最后是討論.

2 方法原理

2.1 低軌衛(wèi)星預(yù)報(bào)誤差的攝動(dòng)源分析

對(duì)于低軌衛(wèi)星而言,軌道的短期預(yù)報(bào)主要受地球形狀、大氣阻力和太陽光壓攝動(dòng)的影響.地球形狀攝動(dòng)主要作用在徑向,大氣阻力攝動(dòng)主要作用在沿跡方向,表現(xiàn)為沿跡誤差的短周期變化以及接近二次多項(xiàng)式規(guī)律的發(fā)散.此外,大氣阻力攝動(dòng)還會(huì)引起能量的顯著耗散,使得軌道半長軸不斷衰減,在徑向誤差上表現(xiàn)出一定程度的線性發(fā)散.至于預(yù)報(bào)誤差的法向分量,主要由地球形狀和太陽光壓攝動(dòng)引起,總體量級(jí)很小,在短期預(yù)報(bào)中可以忽略不計(jì).

圖1–2分別給出了400 km左右的GRACE(Gravity Recovery and Climate Experiment)-A衛(wèi)星[3]在平靜、劇烈兩種太陽活動(dòng)水平下,24 h軌道預(yù)報(bào)誤差的典型實(shí)例.為了避免觀測資料的誤差以及觀測幾何的影響,凸顯出方法自身的誤差,我們選擇GRACEA衛(wèi)星全弧段星載GPS測量處理后的數(shù)據(jù)RSO(Rapid Science Orbit)為觀測量,精度在分米量級(jí).

圖1是2010年2月3日的預(yù)報(bào)結(jié)果,這一天太陽輻射指數(shù)為75 sfu(1 sfu=1× 10?22W·m?2·Hz?1),日平均地磁指數(shù)Ap為10.徑向、沿跡和軌道面法向的24 h最大誤差分別為1.9 m,102 m和0.6 m.沿跡方向的誤差最為顯著,整體呈線性遞增趨勢,同時(shí)具備時(shí)間的二次項(xiàng)特征.徑向誤差有明顯的短周期變化,伴隨著線性發(fā)散特征(大氣攝動(dòng)的能量耗散效應(yīng)).在軌道面法向,整體呈現(xiàn)不規(guī)則振幅的短周期變化,量級(jí)小于1 m.以上誤差特征均符合前文對(duì)攝動(dòng)源的分析.

圖2是2014年10月23日的預(yù)報(bào)結(jié)果,這一天太陽輻射指數(shù)為216 sfu,日平均地磁指數(shù)Ap為8.在徑向、沿跡和軌道面法向,24 h最大誤差分別為23m,855m和0.9m.誤差曲線的形態(tài)與圖1完全一致,只是徑向和沿跡方向誤差有了顯著提高,這主要是由于太陽活動(dòng)增強(qiáng),大氣模型誤差增大引起的.在軌道面法向誤差量級(jí)仍然小于1m,是短期預(yù)報(bào)誤差中的次要部分.

鑒于以上分析,我們以沿跡方向24 h最大誤差來表征低軌衛(wèi)星的預(yù)報(bào)精度,重點(diǎn)討論如何降低沿跡方向的最大預(yù)報(bào)誤差.

圖1 GRACE-A衛(wèi)星的24 h軌道預(yù)報(bào)誤差(2010-02-03,平靜太陽活動(dòng))Fig.1 The 24-h p red iction error of GRACE-A’s orb it(2010-02-03,low so lar activ ity)

圖2 GRACE-A衛(wèi)星的24 h軌道預(yù)報(bào)誤差(2014-10-23,劇烈太陽活動(dòng))Fig.2 T he 24-h p red iction error of GRACE-A’s orbit(2014-10-23,strong so lar activ ity)

2.2 沿跡方向誤差的發(fā)散特征

從理論上看,預(yù)報(bào)過程是一個(gè)對(duì)軌道初值的積分問題,預(yù)報(bào)結(jié)果與積分的初值有關(guān),同時(shí)又與微分方程中的力模型有關(guān).因此預(yù)報(bào)誤差可以分為兩個(gè)部分:初值誤差和模型誤差.根據(jù)衛(wèi)星動(dòng)力學(xué)理論和球面天文相關(guān)關(guān)系,能夠從理論上推導(dǎo)沿跡誤差的發(fā)散規(guī)律,給出誤差的分析表達(dá)式.對(duì)于相同時(shí)刻,用事后精密星歷減去軌道預(yù)報(bào)值,得到三維位置誤差,再按照徑向、沿跡和軌道面法向進(jìn)行分解,沿跡誤差為:

其中r是衛(wèi)星的地心距,a,e,i,?,ω,M是衛(wèi)星的開普勒根數(shù),f是真近點(diǎn)角,p=a(1?e2).在24 h短期預(yù)報(bào)中,e,i,?只有振幅很小的周期變化,這里可以視為常數(shù)(用平均值代替).對(duì)于低軌近圓軌道,經(jīng)過一系列簡化得到t時(shí)刻誤差?Vt的表達(dá)式:

其中λ=ω+M.

其中d t=t?t0,t0為預(yù)報(bào)初始時(shí)刻.

(2)式中第1部分為常數(shù).第2部分為軌道真近點(diǎn)角f的周期函數(shù),圖1和圖2中沿跡方向曲線在線性發(fā)散趨勢上疊加了一種短周期波動(dòng),就是來自這一項(xiàng)的貢獻(xiàn).這一項(xiàng)周期函數(shù)的振幅中含有?e0,而對(duì)于近圓軌道,偏心率e的量級(jí)已經(jīng)低于10?2,?e0一般都小于10?3.因此這部分周期變化的振幅會(huì)非常小,可以視為常數(shù).因此整個(gè)沿跡誤差最大的變化項(xiàng)即為(2)式中的第3部分,其詳細(xì)表達(dá)式為(3)式.那么沿跡誤差隨預(yù)報(bào)時(shí)間d t的發(fā)散規(guī)律可以表示為

鑒于篇幅限制,這里不再給出詳細(xì)推導(dǎo).實(shí)踐中我們常發(fā)現(xiàn),短期預(yù)報(bào)沿跡誤差會(huì)隨著預(yù)報(bào)時(shí)間d t遞增(正向或負(fù)向),呈線性或者二次多項(xiàng)式規(guī)律發(fā)散.那么(4)、(5)式恰恰就反映了這一點(diǎn):軌道初值的半長軸誤差?a0對(duì)線性發(fā)散有貢獻(xiàn),模型誤差?˙a對(duì)二次多項(xiàng)式發(fā)散有貢獻(xiàn).

2.3 抑制沿跡誤差發(fā)散的方法

既然沿跡誤差的發(fā)散服從(4)或(5)式的規(guī)律,因此就可以從這個(gè)關(guān)系式中尋找一種“特殊條件”來抑制沿跡誤差的發(fā)散,這是本文的思路.

從簡單入手,仍可以把短期預(yù)報(bào)的沿跡誤差看成按照一次線性規(guī)律發(fā)散,那么(5)式改寫為:

于是,誤差發(fā)散的“等效斜率”為:

關(guān)于?a0,可以通過定軌段最后一個(gè)歷元與事后星歷直接比較給出.為了保證結(jié)果可靠,防止最后一個(gè)歷元的跳變,也可以用定軌段最后一個(gè)軌道周期內(nèi)的?a0的加權(quán)平均值代替.關(guān)于?˙a,是真實(shí)的半長軸衰減變率減去模型計(jì)算的衰減變率,即?=true?model,主要來源是動(dòng)力學(xué)模型誤差.因?yàn)?/p>

式中S是衛(wèi)星飛行方向有效橫截面積,m是衛(wèi)星質(zhì)量,F=(1?cos i·rne/v)是與大氣旋轉(zhuǎn)速度ne有關(guān)的一個(gè)量,變化很小[4],因?yàn)樵跇?biāo)準(zhǔn)單位下ne=5.883×10?2.對(duì)于近圓軌道,上式在假設(shè)e=0后,能進(jìn)一步簡化為:

而且(9)式中的n和a可以看作常數(shù),為了描述方便,我們記常數(shù)項(xiàng)的組合為:

那么(9)式改寫為:

結(jié)合(7)式和(12)式,可以發(fā)現(xiàn)誤差發(fā)散的速度(等效斜率k)與6個(gè)參數(shù)有關(guān):初值誤差?a0,真實(shí)的阻力系數(shù)CDtrue,真實(shí)的大氣密度ρtrue,預(yù)報(bào)段阻力系數(shù)CDmodel,模型計(jì)算大氣密度ρmodel,以及預(yù)報(bào)時(shí)長d t.

其余的量(如n,a,η)都可視為常數(shù),用定軌段的平均值代替.要抑制沿跡誤差的發(fā)散,可以令“等效斜率”為0,即:

也就是說,可以通過調(diào)整預(yù)報(bào)段采用的CDmodel使得上式成立,從而控制誤差的發(fā)散速度,降低最大誤差.因此(13)式給出的方程可表示為

在(14)式的6個(gè)參數(shù)中,CDmodel是待求參數(shù),在剩下的5個(gè)參數(shù)中有2個(gè)是可以很快確定的:初值誤差?a0可以在定軌段給出;d t是預(yù)報(bào)時(shí)長,對(duì)于本文的研究取24 h.下面將分別討論剩余3個(gè)未知量(CDtrue,ρtrue,ρmodel)的計(jì)算方法.

真實(shí)的阻力系數(shù)CDtrue:大氣對(duì)衛(wèi)星的阻力系數(shù)是不斷變化的,如前文所述,它與很多物理參數(shù)有關(guān),在空間中是幾乎不可能實(shí)時(shí)求出的,我們可以按照通常的平均化處理,即認(rèn)為在高層大氣中阻力系數(shù)平均值近似為2.2.對(duì)于GRACE-A衛(wèi)星,考慮到吸收面質(zhì)比的誤差,我們?nèi)v史軌道擬合求解CD的長期平均值2.59作為近似.由于這個(gè)參數(shù)與ρtrue是耦合的,它們的乘積共同起作用,所以用定值近似是可行的.

真實(shí)的大氣密度ρtrue:這個(gè)量在預(yù)報(bào)階段會(huì)隨著時(shí)間、空間和空間環(huán)境產(chǎn)生變化,也是很難預(yù)知的.我們認(rèn)為大氣密度雖然變化很復(fù)雜,但是總體的平均密度在相鄰兩天仍然會(huì)保持連續(xù)緩變的特征,因此采用定軌段的真實(shí)大氣密度平均值true代替.

模型計(jì)算的大氣密度ρmodel:由于預(yù)報(bào)段的軌道本身就是未知的,就無法用模型逐點(diǎn)給出嚴(yán)格的模型密度,所以也采用定軌段的模型密度平均值model代替,只要定軌成功,將軌道代入大氣模型計(jì)算,然后求平均密度即可.

對(duì)于圖3的例子,定軌后可以求得以下信息

將這些值代入(15)式即可求出定軌段平均大氣密度

經(jīng)過以上分析,(14)式中6個(gè)參數(shù)有5個(gè)是可以給定的,只剩下1個(gè)待估參數(shù)CDmodel通過求解方程得到,此時(shí)的CDmodel就是應(yīng)用在預(yù)報(bào)段中的阻力系數(shù).這種新的阻力系數(shù)計(jì)算方法充分考慮初值誤差和模型誤差,用(13)式的規(guī)律尋求一種特殊條件,通過對(duì)預(yù)報(bào)段的CD進(jìn)行調(diào)節(jié),對(duì)初值和模型的誤差進(jìn)行“抑制”或者“抵消”,使得誤差發(fā)散斜率盡量趨近于0,從而保證預(yù)報(bào)精度.這就是新方法與傳統(tǒng)方法的區(qū)別.

2.4 阻力系數(shù)計(jì)算步驟

前文給出了計(jì)算阻力系數(shù)的基本原理,這里進(jìn)一步詳細(xì)闡明在定軌預(yù)報(bào)過程中阻力系數(shù)的計(jì)算步驟.問題可以描述為:以低軌衛(wèi)星全弧段高精度GPS位置測量資料為基礎(chǔ),進(jìn)行短期預(yù)報(bào),對(duì)于某一歷元T0,取前24 h的資料定軌,預(yù)報(bào)該歷元之后24 h的軌道,再與事后精密測量作比對(duì),給出沿跡誤差,表征預(yù)報(bào)精度.

圖3 GRACE-A衛(wèi)星定軌段沿跡方向的殘差曲線(2014-10-22)Fig.3 The residual error of GRACE-A in the along-track direction after orb it determ ination(2014-10-22)

STEP1:利用歷元T0時(shí)刻之前24 h的資料進(jìn)行第1次定軌,目的是給出T0時(shí)刻的軌道作為預(yù)報(bào)初值,同時(shí)給出?a0.求解的參數(shù)有:衛(wèi)星位置、速度狀態(tài)量,大氣阻力系數(shù),沿跡方向經(jīng)驗(yàn)攝動(dòng)加速度,不求解光壓系數(shù),采用事后實(shí)測的太陽和地磁指數(shù).說明:因?yàn)榈蛙壭l(wèi)星光壓攝動(dòng)量級(jí)比大氣阻力攝動(dòng)低1–2個(gè)量級(jí),是次要量,短期資料定軌往往不能較好地解算光壓系數(shù),會(huì)出現(xiàn)負(fù)值等異常值.解算沿跡方向經(jīng)驗(yàn)加速度是為了盡可能好地?cái)M合軌道,降低初值誤差?a0.

STEP2:在定軌段(T0?1,T0)中取多個(gè)歷元Ti作為定軌初始?xì)v元(如圖4中所示,可以以3 h為1個(gè)間隔),進(jìn)行二次定軌,目的是給出定軌段的模型平均密度odel.根據(jù)Ti時(shí)刻的?用(13)式求出rue,η.解算參數(shù)有:衛(wèi)星位置、速度狀態(tài)量,大氣阻力系數(shù),不解算光壓系數(shù)和經(jīng)驗(yàn)攝動(dòng)力.

圖4 定軌預(yù)報(bào)問題示意圖Fig.4 The d iagram of the p rob lem of orb it determ ination and p red iction

STEP4:利用STEP1中定軌后T0時(shí)刻的軌道作為預(yù)報(bào)初值,采用新方法計(jì)算的CDmodel進(jìn)行24 h預(yù)報(bào),再與事后高精度星歷作比對(duì),給出沿跡誤差.太陽和地磁指數(shù)使用預(yù)報(bào)值.

表1給出2014年10月22日資料選擇多組定軌初始?xì)v元計(jì)算得到的一系列參數(shù),以及最后預(yù)測23日的大氣阻力系數(shù).可以看出,取不同初始?xì)v元時(shí)?,ue變化較大,odel,η的值變化很小,定軌段變化和最后求解得到的預(yù)報(bào)段CDmodel變化也比較明顯,因此有必要選擇多組初值計(jì)算,得到更穩(wěn)健的結(jié)果.

表1 取不同定軌初始?xì)v元時(shí)各種參數(shù)的計(jì)算結(jié)果Tab le 1 The resu lts of various param eters with d ifferen t in itia l ep oches of orb it d eterm ination

2.5 算法設(shè)計(jì)中的幾個(gè)約定

為了方法更為嚴(yán)謹(jǐn),在算法設(shè)計(jì)和編程上有如下兩個(gè)約定:

(1)在STEP2中的二次定軌環(huán)節(jié),使用的太陽、地磁指數(shù)設(shè)定為定軌段實(shí)測值和預(yù)報(bào)段預(yù)報(bào)值二者的平均值,目的是使得rue,odel更加貼近預(yù)報(bào)段的實(shí)際空間環(huán)境條件,同時(shí)兼顧定軌段的資料.

(2)二次定軌中,在Ti–T0之間選擇定軌弧長為最大整數(shù)倍的軌道周期,因?yàn)榇藭r(shí)定軌的目的在于求rue,odel,其定義應(yīng)該為軌道周期內(nèi)的大氣密度平均值,而非全部時(shí)間段的平均密度,這樣的約定使得“平均密度”的物理意義更清晰,有利于推廣到預(yù)報(bào)段.

3 阻力系數(shù)應(yīng)用于軌道預(yù)報(bào)的效果分析

3.1 分析方法

為了驗(yàn)證新的阻力系數(shù)計(jì)算方法能否提高預(yù)報(bào)精度,需要與傳統(tǒng)方法進(jìn)行對(duì)比.對(duì)于任意歷元(一般設(shè)定為某天的世界時(shí)零時(shí)),取前一天的資料定軌預(yù)報(bào)該天的軌道.

傳統(tǒng)方法:在定軌中求解衛(wèi)星位置、速度狀態(tài)量,大氣阻力系數(shù),太陽光壓系數(shù),然后用定軌后最后一個(gè)歷元的軌道作為初值,并將阻力和光壓系數(shù)代入預(yù)報(bào)環(huán)節(jié),進(jìn)行24 h軌道預(yù)報(bào).

新方法:在定軌中求解衛(wèi)星位置、速度狀態(tài)量,大氣阻力系數(shù),沿跡方向經(jīng)驗(yàn)加速度,求得最后一個(gè)歷元的軌道作為初值,以及該初值所對(duì)應(yīng)的半長軸誤差?a0,然后按照第2節(jié)的方法計(jì)算得到適用于預(yù)報(bào)的大氣阻力系數(shù),進(jìn)行24 h軌道預(yù)報(bào).

3.2 個(gè)例分析

為了直觀展示兩種方法的預(yù)報(bào)精度,我們按照不同太陽輻射等級(jí)選擇3組預(yù)報(bào)實(shí)例.

(1)太陽活動(dòng)平靜期:2010年3月11日

利用3月10日的資料定軌,預(yù)報(bào)3月11日的軌道.3月10日的太陽輻射指數(shù)為80.3 sfu,日平均地磁指數(shù)Ap為10,3 h最大地磁指數(shù)ap為22.3月11日的太陽輻射指數(shù)為84.2 sfu,日平均地磁指數(shù)Ap為9,3 h最大ap為22.總體來看這兩天太陽輻射水平很低,地磁存在小擾動(dòng).

圖5是兩種方法預(yù)報(bào)24 h的沿跡誤差情況.可以看出傳統(tǒng)方法使用的預(yù)報(bào)段阻力系數(shù)為2.553,而新方法計(jì)算的阻力系數(shù)為3.883,兩者差別顯著.24 h預(yù)報(bào)最大誤差分別為288 m和84 m,新方法的預(yù)報(bào)誤差顯著降低了約200 m(約71%),可見新方法計(jì)算的阻力系數(shù)有效抑制了沿跡誤差的發(fā)散.雖然誤差發(fā)散斜率并沒有嚴(yán)格等于0,這是因?yàn)樵摲椒ㄊ褂玫念A(yù)報(bào)段大氣密度(真實(shí)和模型)是來自于定軌段的平均值,近似過程存在著誤差,所以(13)式不能嚴(yán)格成立,但它能顯著降低誤差發(fā)散速度,達(dá)到了降低預(yù)報(bào)誤差的目的.

(2)太陽活動(dòng)中等強(qiáng)度期:2012年5月23日

用5月22日的資料定軌,預(yù)報(bào)5月23日的軌道.5月22日太陽輻射指數(shù)為125 sfu,日平均地磁指數(shù)Ap為19,3 h最大地磁指數(shù)ap為32;5月23日太陽輻射指數(shù)為121 sfu,日平均地磁指數(shù)Ap為16,3 h最大地磁指數(shù)ap為39.總體來看太陽活動(dòng)在中等水平,地磁有顯著擾動(dòng).圖6是兩種方法的預(yù)報(bào)結(jié)果對(duì)比.傳統(tǒng)方法的阻力系數(shù)為3.357,最大預(yù)報(bào)誤差881 m,而新方法計(jì)算的阻力系數(shù)為4.656,最大預(yù)報(bào)誤差為134 m,兩者相比誤差減少了750 m(約85%).盡管新的阻力系數(shù)明顯偏離了典型值2.2,但它充分抑制了初值誤差和模型誤差,提高了預(yù)報(bào)精度.

(3)太陽活動(dòng)劇烈期:2014年10月28日

用10月27日的資料定軌,預(yù)報(bào)10月28日的軌道.10月27日的太陽輻射指數(shù)為216.6 sfu,日平均地磁指數(shù)為10,3 h最大地磁指數(shù)ap為15;10月28日的太陽輻射指數(shù)為187.8 sfu,日平均地磁指數(shù)為12,3 h最大地磁指數(shù)ap為22.總體來看太陽輻射較強(qiáng),地磁存在小擾動(dòng).圖7是兩種方法的軌道預(yù)報(bào)誤差.傳統(tǒng)方法的阻力系數(shù)為2.335,最大誤差2773 m;而新方法計(jì)算的阻力系數(shù)為2.930,最大誤差降為1072 m,精度顯著提升了1700m(約61%).

從以上3個(gè)計(jì)算實(shí)例可以看出,新方法對(duì)不同輻射水平的軌道預(yù)報(bào)都有效,精度也有顯著提升(60%–85%).同時(shí),3個(gè)例子中都伴隨著中小程度的地磁擾動(dòng),可見該方法對(duì)于地磁擾動(dòng)的資料也有一定的適用性.但實(shí)踐中發(fā)現(xiàn),當(dāng)預(yù)報(bào)段發(fā)生磁暴時(shí)(日平均地磁指數(shù)大于48),方法有很大可能失效,因?yàn)榉椒ǖ幕炯僭O(shè)(預(yù)報(bào)段大氣平均密度與定軌段基本接近)就不能成立.

圖5 兩種方法的預(yù)報(bào)誤差對(duì)比(2010-03-11,平靜太陽活動(dòng))Fig.5 The com parison of orbit p red iction errors betw een classicalm ethod and new m ethod(2010-03-11, qu iet so lar activ ity)

圖6 兩種方法的預(yù)報(bào)誤差對(duì)比(2012-05-23,中等強(qiáng)度太陽活動(dòng))Fig.6 The com parison of orbit p red iction errors betw een classicalm ethod and new m ethod(2012-05-23, m oderate solar activity)

圖7 兩種方法的預(yù)報(bào)誤差對(duì)比(2014-10-28,劇烈太陽活動(dòng))Fig.7 The com parison of orbit p red iction errors betw een classicalm ethod and new m ethod(2014-10-28, strong solar activ ity)

(4)兩種方法等效的實(shí)例:2014年4月3日

需要指出的是,在有些預(yù)報(bào)中會(huì)出現(xiàn)傳統(tǒng)方法和新方法基本等效的情況,也就是說傳統(tǒng)方法精度已經(jīng)很高了,新方法并不能繼續(xù)進(jìn)一步提高精度,與傳統(tǒng)方法的精度差異在50m以內(nèi),這時(shí)我們認(rèn)為兩種方法基本等效.例如2014年4月3日的預(yù)報(bào)(見圖8),用傳統(tǒng)方法預(yù)報(bào)時(shí)誤差為72 m,而用新方法誤差為68m,兩者非常接近,主要是因?yàn)閮煞N方法得到的阻力系數(shù)也接近(3.146和3.227).

3.3 結(jié)果統(tǒng)計(jì)

以上展示的僅是幾個(gè)實(shí)例,要觀察方法的應(yīng)用效果還需要更多的資料來驗(yàn)證.我們選擇GRACE-A衛(wèi)星2010—2014年共5 yr的資料進(jìn)行統(tǒng)計(jì)分析,對(duì)比傳統(tǒng)方法和新方法在5 yr內(nèi)的精度(剔除了日地磁指數(shù)超過48的資料).表2統(tǒng)計(jì)了兩種方法24 h預(yù)報(bào)精度的平均值,同時(shí)按照輻射水平進(jìn)行了分類.當(dāng)輻射指數(shù)低于100 sfu,設(shè)定為太陽活動(dòng)平靜期;輻射指數(shù)在100–150 sfu之間,設(shè)定為中等太陽活動(dòng)期;輻射指數(shù)大于150 sfu時(shí),設(shè)定為劇烈太陽活動(dòng)期.從表2中可以看出,新方法相對(duì)于傳統(tǒng)方法,絕對(duì)誤差的平均值下降了48–275m,相對(duì)誤差降低了約44%–51%,且對(duì)于不同太陽輻射等級(jí)的情況都有效.

需要指出的是,表2是先對(duì)同一方法的全部預(yù)報(bào)誤差進(jìn)行統(tǒng)計(jì)平均,然后再進(jìn)行不同方法之間的比較,這反映了一個(gè)方法整體的精度水平.還有一種統(tǒng)計(jì)方法是對(duì)每一天的預(yù)報(bào),一一對(duì)應(yīng)地比較兩種方法的差異,這樣就能統(tǒng)計(jì)出新方法的成功率.

為了描述清晰,我們定義:對(duì)于同一天的預(yù)報(bào),如果新方法的預(yù)報(bào)誤差小于傳統(tǒng)方法預(yù)報(bào)誤差,即認(rèn)為新方法成功;如果新方法誤差大于傳統(tǒng)方法,但是差異小于50m,且相對(duì)差異小于5%,則認(rèn)為兩種方法等效.表3給出了成功率和等效率的統(tǒng)計(jì).在3種輻射水平下,方法的成功率都超過70%,意味著對(duì)不同太陽活動(dòng)水平都有效.總體而言,新方法的成功率約為71%,等效率約為15%,兩者相加約為86%(有效率),在實(shí)際工作中具有較好的應(yīng)用價(jià)值.

圖8 兩種方法的預(yù)報(bào)誤差對(duì)比(2014-04-03,兩種方法等效)Fig.8 The com parison of orbit p red iction errors betw een classicalm ethod and new m ethod(2014-04-03, tw o m ethods are equiva len t)

表2 兩種方法24 h預(yù)報(bào)精度的統(tǒng)計(jì)Tab le 2 T he statistics of the 24-h orb it p red iction accu racy

表3 兩種方法24 h預(yù)報(bào)成功率和等效率的統(tǒng)計(jì)Tab le 3 T he statistics of the ratios of su ccessfu l and equ iva len t in stan ces for the 24-h orb it p red iction

4 結(jié)論與討論

本文在分析低軌衛(wèi)星軌道預(yù)報(bào)誤差發(fā)散規(guī)律的基礎(chǔ)上,提出了一種大氣阻力系數(shù)計(jì)算方法,充分考慮了軌道初值誤差和模型誤差,在定軌段已有信息的基礎(chǔ)上,通過計(jì)算合適的大氣阻力系數(shù)抑制預(yù)報(bào)階段沿跡方向的誤差發(fā)散.結(jié)果表明該方法能有效提高24 h短期預(yù)報(bào)精度:相對(duì)于傳統(tǒng)的定軌預(yù)報(bào)方法,新方法能提高預(yù)報(bào)精度約45%,成功率約71%,總體有效率約86%.數(shù)值試驗(yàn)還發(fā)現(xiàn),該方法對(duì)低、中、高3種輻射水平均有效,對(duì)于中低等級(jí)的地磁擾動(dòng)也有效,具備較好的應(yīng)用價(jià)值.

值得指出的是,該方法只適用于短期的軌道預(yù)報(bào)(比如24 h以內(nèi)),因?yàn)榇朔椒ㄓ幸粋€(gè)前提假設(shè),即認(rèn)為相鄰兩天的大氣密度環(huán)境是連續(xù)和緩變的,從定軌段提取的“平均大氣密度”(模型值和真實(shí)值)可以近似推廣到預(yù)報(bào)段.一旦預(yù)報(bào)時(shí)間過長,以上假設(shè)顯然不再成立,定軌段提取的信息與預(yù)報(bào)段不能匹配適用,所預(yù)測的大氣阻力系數(shù)就是不正確的,因此該方法不適用于中長期預(yù)報(bào).最后需要說明,本文所有試驗(yàn)是針對(duì)全弧段高精度GPS測量資料開展的,資料精度高、分辨率高，弧段覆蓋好,一定程度上保證了方法的可靠性,而對(duì)于光學(xué)、雷達(dá)等較稀疏的資料,仍需要進(jìn)一步的檢驗(yàn).

[1]Va llado D A,Fink lem an D.A cA au,2014,95:141

[2]M oe K,M oe M M.P&SS,2005,53:793

[3]Tap ley B D,Bettadpu r S,W atk ins M,et al.GeoRL,2004,31:L 09607

[4]劉林.人造地球衛(wèi)星軌道力學(xué).北京:高等教育出版社,1992:398

The M ethod for Calcu lating A tm ospheric D rag Coefficient Based on the Characteristics of A long-track E rror in LEO O rb it P red iction

WANG Hong-bo1,2,3ZHAO Chang-yin1,2LIU Zhong-gui4ZHANG Wei1,2

(1 Pu rp le M oun tain O bserva to ry,Chinese A cadem y of Scien ces,Nan jing 210008) (2 K ey Labo ra to ry of Space O bjec t and Debris O bserva tion,Pu rp le M oun tain O bserva to ry,Chinese A cadem y of Scien ces,Nan jing 210008) (3 Sta te K ey Labo ra to ry of A stronau tic D ynam ics,X i’an 710043) (4 Beijing In stitu te of T racking an d Telecomm un ication s Techno logy,Beijing 100094)

The errors of atmosphere densitym odeland d rag coefficient are them ajor factors to restrain the accuracy of orbit prediction for the LEO(Low Earth Orbit) ob jects,which would affect un favorably the spacem issions that need a high-precision orbit.This paper brings out a new method for calculating the drag coefficient based on the divergence laws of prediction error’s along-track com ponent.Firstly,we deduce the expression of along-track error in LEO’s orbit prediction,revealing the comprehensive effect of the initial orbit and model’s errors in the along-track direction.According to this expression,we work out a suitable drag coefficient adopted in prediction step on the basis of some certain information from orbit determ ination step,which w ill lim it the increasing rate of along-track error and reduce the largest error in this direction, then achieving the goal of im proving the accuracy of orbit prediction.In order to verify themethod’s accuracy and successful rate in the practice of orbit prediction,we use the full-arcs high precision position data from the GPS receiver on GRACE-A.The result shows that this new method can significantly improve the accuracy of prediction by about 45%,achieving a successful rate of about 71%and an effective rate of about 86%,with respect to classicalmethod which uses the fitted drag coefficient directly from orbit determ ination step.Furtherm ore,the new m ethod shows a p referable app lication value,because it is effective for low,moderate,and high solar radiation levels,as well as some quiet and moderate geomagnetic activity condition.

celestialm echanics:orbit calcu lation and determ ination,methods:numerical,upper-atmosphere drag

P135;

A

10.15940/j.cnki.0001-5245.2016.04.006

2015-11-14收到原稿,2015-12-03收到修改稿

?宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放基金資助

?whb@pm o.ac.cn