張煒，崔文，張育衛(wèi)，劉興，李菊清

(西安衛(wèi)星測控中心，西安710600)

1 引言

大氣阻力是即將再入空間目標所受的最主要非保守攝動力，精確的空間目標面質(zhì)比及合理的大氣阻力特性建模是準確計算大氣阻力加速度、預報再入時間的關(guān)鍵。大氣阻力系數(shù)與空間目標外形、表面材料、大氣組成及溫度等密切相關(guān)，同一個目標的大氣阻力系數(shù)在不同軌道高度、不同太陽活動水平等情況下會有較大差異[1]。通常，空間目標的準確外形、質(zhì)量、姿態(tài)和表面材料等都是未知的，分別確定大氣阻力系數(shù)、迎風面積和質(zhì)量難度很大，因此引入彈道系數(shù)B進行統(tǒng)一處理[2]。彈道系數(shù)的定義為:

式中，CD為大氣阻力系數(shù)；A為迎風面積；m為質(zhì)量。

在大氣阻力的影響下，空間目標的軌道不斷衰減，軌道形狀越來越圓，最終墜入稠密大氣層。按初始遠地點高度Ha_ini(本文的 “初始”指空間目標再入前10天)對再入目標軌道進行分類:

(1)近圓軌道:Ha_ini＜500km

(2)小橢圓軌道:500km≤Ha_ini＜5000km

(3)大橢圓軌道:Ha_ini≥5000km

依此統(tǒng)計，2012—2018年非受控再入的大型空間目標數(shù)量，結(jié)果如圖1所示?？梢钥闯?，每年再入的大型空間目標中小橢圓軌道目標不在少數(shù)，近7年平均占1/6以上，2018年更是超過40%。

圖1 2012年至2018年非受控再入的大型空間目標數(shù)量Fig.1 Number of large uncontrolled reentry objects from 2012 to 2018

再入預報的主要難點在于軌道確定和大氣阻力建模[3]，一般在足夠測量數(shù)據(jù)的基礎上確定精確軌道并解算彈道系數(shù)，但是只有少數(shù)國家和組織具有獲取再入目標測量數(shù)據(jù)的能力。即使可以，相鄰圈次的間隔可能太長導致難以聯(lián)合使用確定精密軌道，容易出現(xiàn)定軌殘差太大、不收斂等問題，這也是本文 “稀疏數(shù)據(jù)”這一基本背景?！暗厍蛑亓龊秃Ｑ蟓h(huán)流探測衛(wèi)星”(GOCE)再入前20天起，布設于德國Wachtberg的跟蹤及成像雷達 (TIRA)共獲取12圈次測量數(shù)據(jù)，不同數(shù)據(jù)組合的定軌結(jié)果明顯差異[4]，在沒有基準星歷的情況下難以選優(yōu)。目前大部分再入預報研究只能依靠美國戰(zhàn)略司令部于Space-track網(wǎng)站公開發(fā)布的兩行根數(shù) (TLE)，基于TLE的再入預報方法研究主要聚焦于TLE的預處理[5，6]，或利用TLE解算彈道系數(shù)、太陽光壓系數(shù)、狀態(tài)矢量等[7-11]。幾乎所有的研究均假設再入目標的彈道系數(shù)是一個常數(shù)，顯然這是不夠準確的，尤其對于橢圓軌道再入目標而言，容易忽略彈道系數(shù)變化帶來的問題。

相較于近圓軌道目標，小橢圓軌道目標的再入預報難度更大。主要原因是小橢圓軌道目標再入過程中幾乎每圈都要穿過整個稠密大氣層，高度變化導致同一軌道周期內(nèi)彈道系數(shù)變化明顯。實際應用中發(fā)現(xiàn)，若采用整個數(shù)據(jù)弧段內(nèi)解算單個彈道系數(shù)的策略，各測量元素的殘差可能較大，難以判斷軌道確定結(jié)果的準確性，而采用分段解算多個彈道系數(shù)的策略，雖然可以獲得更加精確的目標軌道，但是進行再入預報時彈道系數(shù)的初值難以選擇。

本文針對稀疏數(shù)據(jù)情況和小橢圓軌道的特點，提出使用平均彈道系數(shù)進行再入預報。利用半數(shù)值法計算平均彈道系數(shù)、數(shù)值法進行短弧數(shù)據(jù)軌道確定和軌道外推，并使用實測數(shù)據(jù)進行效果驗證。結(jié)果表明，本文提出的方法可有效避免稀疏數(shù)據(jù)情況下小橢圓軌道目標軌道確定收斂難或殘差過大、再入預報時彈道系數(shù)難確定等問題，預報精度達到與近圓軌道目標同等水平。

2 方法

本文針對稀疏數(shù)據(jù)情況和小橢圓軌道的特點，提出半數(shù)值法和數(shù)值法相結(jié)合預報再入時刻。 使用經(jīng)典軌道根數(shù)σ(a，e，i，Ω，ω，M) 描述再入目標軌道， (a，e，i，Ω，ω，M) 分別為半長軸、 偏心率、傾角、升交點赤經(jīng)和平近點角。

2.1 軌道確定

稀疏數(shù)據(jù)情況是指測量數(shù)據(jù)圈次太少，例如24h內(nèi)僅有1～2圈測量數(shù)據(jù)，多圈數(shù)據(jù)聯(lián)合確定軌道時容易出現(xiàn)殘差太大、不收斂等問題。本文使用短弧數(shù)據(jù)定軌方法對單圈數(shù)據(jù)分別確定軌道，考慮的攝動項包括:地球非球形引力、大氣阻力和日月引力。短弧數(shù)據(jù)定軌的問題主要在于軌道外推時誤差容易迅速發(fā)散，但是本文僅需數(shù)據(jù)弧段內(nèi)的6個參數(shù)，可有效避免不收斂和誤差發(fā)散等問題。軌道結(jié)果精度主要受限于測量精度，與基于大量測量數(shù)據(jù)的精密軌道相比，短弧數(shù)據(jù)軌道位置誤差一般為幾百米，且主要集中在對再入預報影響最小的沿跡方向[12]。

基于單圈測量數(shù)據(jù)的軌道確定結(jié)果為密切根數(shù) (σ1，σ2，…，σN) ， 本文還用到對應的平根數(shù)。去除密切根數(shù)中的短周期項即可得到平根數(shù)，即:

式中，為平根數(shù)；σ為密切根數(shù)；Δσs為短周期項。平根數(shù)的計算方法及短周期項的詳細討論可參考文獻 [13]—文獻 [16]。沒有可用測量數(shù)據(jù)的情況下可以使用TLE根數(shù)作為偽測量數(shù)據(jù)，使用SGP4/SDP4模型解析TLE歷元時刻的位置速度，并轉(zhuǎn)換成軌道根數(shù)。

2.2 平均彈道系數(shù)

如何計算再入目標的彈道系數(shù)是許多再入預報研究的焦點。在大多數(shù)研究中再入目標的彈道系數(shù)被視為常數(shù)，采用將彈道系數(shù)作為待估參數(shù)的策略，在軌道確定時與狀態(tài)參量(r，r)一同求解，再基于彈道系數(shù)和狀態(tài)參量的解算結(jié)果進行再入預報。但是這種方法對小橢圓軌道目標并不完全適用，由于空間目標的彈道系數(shù) (主要是大氣阻力系數(shù))隨高度變化，如果整個數(shù)據(jù)弧段內(nèi)解算單個彈道系數(shù)，各測量元素的殘差可能較大導致難以判斷軌道確定結(jié)果的準確性；如果分段解算多個彈道系數(shù)，進行再入預報時難以選擇彈道系數(shù)的初值。為避免這些問題，提出基于空間目標軌道周期內(nèi)的彈道系數(shù)平均值進行再入預報。平均彈道系數(shù)反映了一軌道周期內(nèi)大氣阻力的平均效果，定義為:

式中，ρ為大氣密度值。為書寫方便，下文統(tǒng)一寫為B。

大氣阻力僅有長期作用效果，因此進行彈道系數(shù)擬合時可以采用開普勒平根數(shù)作為根數(shù)系統(tǒng)，并使用半數(shù)值法進行軌道積分，積分考慮的攝動項包括地球非球形J2項、J3項、大氣阻力。使用半數(shù)值法時高階攝動力建模較為困難和復雜，且精度改進效果不明顯[17]。積分模型如下[18]:

式中，v為再入目標相對大氣的運動速度；f為真近點角；ωE為地球運動速率；μ為地球引力常數(shù)；rE為地球半徑；n為平運動；p為軌道半通徑；J2、J3分別為地球引力場二階、三階帶諧系數(shù)。本文使用MSIS-90模型計算大氣密度。由于根數(shù)系統(tǒng)不包含短周期項，因此積分步長可以使用軌道周期的整數(shù)倍。

假設已有空間目標連續(xù)N組軌道平根數(shù)(，)，平均彈道系數(shù)的計算過程如下:

(1)以及初始平均彈道系數(shù)B0為初值，一軌道周期為步長對公式 (4)進行反向積分。反向積分的目的是防止彈道系數(shù)初值偏差太大導致積分過程 “鉆地”。

(2)若積分時刻與歷元之差小于半周期，則計算當前時間半長軸計算值與實測值之差Δak及半長軸對彈道系數(shù)的偏導數(shù)根據(jù)公式 (5)計算:

(3)第1條根數(shù)計算結(jié)束后，使用公式 (6)計算彈道系數(shù)修正值:

(4)若ΔB小于收斂閾值，則計算結(jié)束，否則對彈道系數(shù)進行修正:Bj+1=Bj+ΔB，返回 (1)重新計算。

2.3 再入時間預報

以第N圈數(shù)據(jù)確定的軌道密切根數(shù)σN及彈道系數(shù)計算結(jié)果B為初值，用數(shù)值積分的方法進行軌道外推，直到再入目標的平均軌道高度低于80km，積分結(jié)束時間后約5min即為空間目標的再入時間。數(shù)值積分使用10階定步長KSG積分器，用到的攝動力及其模型如表1所示。

表1 再入預報使用的攝動力及模型Table 1 Dynamical and geometrical models employed in the reentry prediction

由于不同大氣模型計算得到的大氣密度值存在差異[19]，因此彈道系數(shù)擬合和再入預報時應使用同一大氣模型，本文均使用MSIS-90模型。

3 數(shù)據(jù)與結(jié)果

使用誤差百分比δ評估再入預報的精度，誤差百分比δ用下式計算:

式中，treal為空間目標的實際再入時間；tpred為預報的再入時間；t0為進行再入預報所用初軌的歷元。從公開資料看，近圓軌道空間目標的再入預報誤差基本在10%～20%之間，部分預報結(jié)果誤差甚至可能超過30%，這主要取決于再入過程中再入目標姿態(tài)和大氣環(huán)境的變化[20，21]。

選擇美國銥星52(NORAD編號為25169)作為算例目標進行驗證，其初始近地點高度約為154km、遠地點高度約為608km。根據(jù)Space-track網(wǎng)站發(fā)布的信息，銥星52的再入時間為2018-11-05 21∶20(UTC)，與本文最后一次預報結(jié)果基本一致，因此treal取為 2018-11-05 21∶20(UTC)。從2018年10月27日起，選擇間隔10h以上的連續(xù)3圈測量數(shù)據(jù)進行再入預報，銥星52的平均彈道系數(shù)及再入時間預報結(jié)果如表2、圖2所示。

為了更直觀地理解平均彈道系數(shù)的計算效果，預報時大氣環(huán)境參數(shù)使用的均為實際觀測值，期間大氣環(huán)境平靜，大氣模型誤差對平均彈道系數(shù)結(jié)果影響較小。從計算結(jié)果看，銥星52的彈道系數(shù)在0.0251～0.0256m2/kg之間變化，均值為0.02548m2/kg，上下浮動1%左右，遠小于大氣阻力系數(shù)軌道周期內(nèi)約10%的變化幅度，驗證了平均彈道系數(shù)計算方法的準確性和有效性。10次再入時間預報結(jié)果中的最大誤差僅為3.25%，其中8次計算的誤差小于1%。

按照本文提出的再入預報方法，選擇TLE作為數(shù)據(jù)源對2017年再入的11個小橢圓軌道目標進行再入預報。再入預報時將5天內(nèi)每天1組TLE組成根數(shù)序列，并將根數(shù)序列視為再入目標的實測軌道值。再入前10天起每天進行一次預報，預報時大氣環(huán)境參數(shù)使用預報值，具有不確定性，各目標10次預報中的實際最大預報誤差如圖3所示?？梢钥闯?，11個目標中9個目標的再入預報誤差小于15%，預報誤差與初始偏心率 (或初始遠地點高度)無明顯相關(guān)性，說明本文方法具有較好的適用性。另外兩個分別為35.6% (目標NORAD編號為41726)和42.5%(目標NORAD編號為09980)，由于TLE使用前未作異常值篩選，分析認為異常TLE是造成這兩個目標最大誤差偏大的主要原因。從圖3也可以看到，使用TLE進行再入預報也能獲得較高的預報精度。

表2 銥星52再入預報結(jié)果Table 2 Reentry prediction results of Iridium-52

圖2 銥星52的再入預報誤差百分比及彈道系數(shù)Fig.2 Relative error and mean ballistic coefficient of Iridium-52

圖3 2017年再入的11個小橢圓軌道目標的最大預報誤差Fig.3 Maximum relative error of 11 reentry objects with low-eccentricity orbits in 2017

4 結(jié)論

由于大氣阻力系數(shù)隨高度變化，小橢圓軌道目標的再入預報比近圓軌道目標難度更大。本文針對稀疏數(shù)據(jù)情況和小橢圓軌道的特點，提出使用平均彈道系數(shù)進行再入預報，并使用再入目標的實測數(shù)據(jù)進行驗證。結(jié)果表明，基于平均彈道系數(shù)進行小橢圓軌道目標再入預報，可有效降低再入預報誤差，預報精度達到與近圓軌道目標同等水平。