汪宏波 趙長印 柳仲貴 張偉

(1中國科學院紫金山天文臺 南京 210008) (2中國科學院空間目標與碎片觀測重點實驗室 南京 210008) (3宇航動力學國家重點實驗室 西安 710043) (4北京跟蹤與通信技術研究所 北京 100094)

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

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

(1中國科學院紫金山天文臺 南京 210008) (2中國科學院空間目標與碎片觀測重點實驗室 南京 210008) (3宇航動力學國家重點實驗室 西安 710043) (4北京跟蹤與通信技術研究所 北京 100094)

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

天體力學:軌道計算和定軌,方法：數值,高層大氣阻力

1 引言

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

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

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

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

2 方法原理

2.1 低軌衛(wèi)星預報誤差的攝動源分析

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

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

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

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

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

圖1 GRACE-A衛(wèi)星的24 h軌道預報誤差(2010-02-03,平靜太陽活動)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軌道預報誤差(2014-10-23,劇烈太陽活動)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ā)散特征

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

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

其中λ=ω+M.

其中d t=t?t0,t0為預報初始時刻.

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

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

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

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

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

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

關于?a0,可以通過定軌段最后一個歷元與事后星歷直接比較給出.為了保證結果可靠,防止最后一個歷元的跳變,也可以用定軌段最后一個軌道周期內的?a0的加權平均值代替.關于?˙a,是真實的半長軸衰減變率減去模型計算的衰減變率,即?=true?model,主要來源是動力學模型誤差.因為

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

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

那么(9)式改寫為:

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

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

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

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

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

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

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

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

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

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

2.4 阻力系數計算步驟

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

圖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時刻之前24 h的資料進行第1次定軌,目的是給出T0時刻的軌道作為預報初值,同時給出?a0.求解的參數有:衛(wèi)星位置、速度狀態(tài)量,大氣阻力系數,沿跡方向經驗攝動加速度,不求解光壓系數,采用事后實測的太陽和地磁指數.說明:因為低軌衛(wèi)星光壓攝動量級比大氣阻力攝動低1–2個量級,是次要量,短期資料定軌往往不能較好地解算光壓系數,會出現負值等異常值.解算沿跡方向經驗加速度是為了盡可能好地擬合軌道,降低初值誤差?a0.

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

圖4 定軌預報問題示意圖Fig.4 The d iagram of the p rob lem of orb it determ ination and p red iction

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

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

表1 取不同定軌初始歷元時各種參數的計算結果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 算法設計中的幾個約定

為了方法更為嚴謹,在算法設計和編程上有如下兩個約定:

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

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

3 阻力系數應用于軌道預報的效果分析

3.1 分析方法

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

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

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

3.2 個例分析

為了直觀展示兩種方法的預報精度,我們按照不同太陽輻射等級選擇3組預報實例.

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

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

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

(2)太陽活動中等強度期:2012年5月23日

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

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

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

從以上3個計算實例可以看出,新方法對不同輻射水平的軌道預報都有效,精度也有顯著提升(60%–85%).同時,3個例子中都伴隨著中小程度的地磁擾動,可見該方法對于地磁擾動的資料也有一定的適用性.但實踐中發(fā)現,當預報段發(fā)生磁暴時(日平均地磁指數大于48),方法有很大可能失效,因為方法的基本假設(預報段大氣平均密度與定軌段基本接近)就不能成立.

圖5 兩種方法的預報誤差對比(2010-03-11,平靜太陽活動)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 兩種方法的預報誤差對比(2012-05-23,中等強度太陽活動)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 兩種方法的預報誤差對比(2014-10-28,劇烈太陽活動)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)兩種方法等效的實例:2014年4月3日

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

3.3 結果統(tǒng)計

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

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

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

圖8 兩種方法的預報誤差對比(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預報精度的統(tǒng)計Tab le 2 T he statistics of the 24-h orb it p red iction accu racy

表3 兩種方法24 h預報成功率和等效率的統(tǒng)計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 結論與討論

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

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

[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)星軌道力學.北京:高等教育出版社,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收到修改稿

?宇航動力學國家重點實驗室開放基金資助

?whb@pm o.ac.cn