一種LSTM模型預(yù)測BC值的空間碎片無控再入預(yù)報方法

蒯家偉，趙柯昕，孫立剛，廖名傳

(1. 南京航空航天大學(xué)航天學(xué)院，南京 211106； 2. 北京電子工程總體研究所，北京 100854；3. 北京遙測技術(shù)研究所，北京 100094； 4. 中國科學(xué)院微小衛(wèi)星創(chuàng)新研究院，上海 201304)

0 引 言

低地球軌道(LEO)空間碎片在大氣阻力耗散作用下，軌道能量降低并再入大氣層，進(jìn)而解體隕落。空間碎片再入有利于降低駐空物體(RSO)碰撞風(fēng)險，維持空間環(huán)境物體數(shù)量穩(wěn)定，由于其中存在大型碎片如運(yùn)載火箭上面級、燃料儲箱等，可能因再入大氣層時燃燒不充分而撞擊地面，對地面人類財(cái)產(chǎn)等帶來安全隱患[1]，因而對空間碎片的再入時刻預(yù)報對空間環(huán)境的治理與地面安全的保障有著重要意義，為了降低這類事故風(fēng)險需要對空間碎片再入時刻進(jìn)行高精度的預(yù)測。現(xiàn)階段空間碎片再入時刻預(yù)報誤差主要來源于大氣阻力的估計(jì)誤差。大氣阻力攝動fdrag主要與碎片的大氣阻力系數(shù)Cd,面積質(zhì)量比A/m,大氣密度ρ和相對于大氣速度矢量v有關(guān)，其表達(dá)式如式(1)所示：

(1)

式中：v=|v|，其中大氣密度受太陽和地磁活動影響，存在較強(qiáng)的短期時變特征，較長時間預(yù)報導(dǎo)致軌道誤差較大，因此空間碎片再入時刻精確預(yù)測一直是空間碎片研究的重點(diǎn)和難點(diǎn)[2]。

國際機(jī)構(gòu)間空間碎片協(xié)調(diào)委員會(IADC)自1998年以來為了提高空間碎片再入預(yù)報水平，組織開展了27次再入隕落預(yù)報聯(lián)測活動，從歷史結(jié)果看總體統(tǒng)計(jì)精度水平在30%～10%之間[3]。在預(yù)測方法方面，通過實(shí)時衛(wèi)星數(shù)據(jù)矯正密度模型提高大氣密度預(yù)測值[4]，由高精度大氣模型得到較精確的相對大氣速度[5]等方法可一定程度上提高短期預(yù)報精度，但因空間碎片外型、材質(zhì)、運(yùn)動復(fù)雜多樣，這類方法難以適用。空間碎片大氣阻力系數(shù)和面質(zhì)比參數(shù)主要受大氣密度、材質(zhì)、碎片姿態(tài)的影響，直接測量阻力系數(shù)和面質(zhì)比參數(shù)十分困難。因此在難以確定面積質(zhì)量的情形下，通常軌道預(yù)報將阻力系數(shù)和面質(zhì)比耦合成一個系數(shù)，稱為彈道系數(shù)(Ballistic coefficient，BC)，本文在公式中以β表示。

(2)

BC的值可以在具有雷達(dá)、光學(xué)等觀測條件下通過精密定軌進(jìn)行估算，從而提高再入預(yù)測精度。而當(dāng)信息不充分時，要提升空間碎片再入彈道系數(shù)估算精度是很困難的。對于無觀測設(shè)備、或者觀測能力不足的情況下，BC的估算還可以通過機(jī)構(gòu)(space-track.org)發(fā)布的兩行根數(shù)進(jìn)行擬合預(yù)測，這是當(dāng)前IADC國際聯(lián)測中各機(jī)構(gòu)基于space-track.org數(shù)據(jù)的主要方法。以Pardini和Anselmo提出基于BC序貫修正預(yù)測BC值的再入時刻預(yù)測方法為例，其BC值的不確定度為±20%[6]，再入時刻預(yù)測偏差接近20%。

近年來隨著人工智能(AI)技術(shù)的發(fā)展，深度學(xué)習(xí)(DL)如深度神經(jīng)網(wǎng)絡(luò)(DNN)、卷積神經(jīng)網(wǎng)絡(luò)(CNNs)、長短周期記憶神經(jīng)網(wǎng)絡(luò)(Long short term memory neural networks，LSTM)等方法被廣泛應(yīng)用于參數(shù)預(yù)測、建模。航天任務(wù)與空間科學(xué)中人工智能算法研究近年來也是熱門課題，涉及地磁活動性指數(shù)、熱層大氣密度、電離層電子含量等在內(nèi)的多種參量都使用深度學(xué)習(xí)法進(jìn)行了較高精度地預(yù)測或建模[7]，包括與本文密切相關(guān)的對大氣密度存在重要影響的參數(shù)地磁指數(shù)的預(yù)測，Kp指數(shù)[8]和Dst指數(shù)[9]都可以通過神經(jīng)網(wǎng)絡(luò)算法進(jìn)行預(yù)報；Chen等[10]開發(fā)出R-DCGAN模型實(shí)現(xiàn)了TEC Map缺失數(shù)據(jù)的補(bǔ)全，以及全球預(yù)測和未來預(yù)測，受此啟發(fā)，Pan等[11]設(shè)計(jì)開發(fā)了DCGAN-PB；隨后，Liu等[12]和Tang等[13]也分別使用長短周期記憶神經(jīng)網(wǎng)絡(luò)預(yù)測了常態(tài)及磁暴期間電離層TEC。張君彪等[14]提出了一種基于集合經(jīng)驗(yàn)?zāi)B(tài)分解和注意力長短時記憶神經(jīng)網(wǎng)絡(luò)的高超聲速滑行飛行器軌跡智能預(yù)測方法，對氣動加速度進(jìn)行實(shí)時估計(jì)，進(jìn)而智能預(yù)測飛行軌跡。此外在空間物體有控再入以及路徑規(guī)劃方面，深度學(xué)習(xí)的應(yīng)用也取得了較多的研究成果，周大鵬等[15]設(shè)計(jì)了一種智能神經(jīng)網(wǎng)絡(luò)自適應(yīng)復(fù)合控制方法，基于誤差反饋學(xué)習(xí)準(zhǔn)則在線更新神經(jīng)網(wǎng)絡(luò)權(quán)重以自適應(yīng)提升空天飛機(jī)再入精度。相比于傳統(tǒng)方法，深度學(xué)習(xí)可以從大量觀測樣本中提取有價值的特征，并利用多神經(jīng)網(wǎng)絡(luò)層重現(xiàn)相關(guān)特征，從而實(shí)現(xiàn)狀態(tài)預(yù)測。

本工作提出了一種利用兩行軌道根數(shù)(Two line element, TLE)數(shù)據(jù)獲得精確彈道系數(shù)BC的深度學(xué)習(xí)方法，采用長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM模型對其進(jìn)行預(yù)測，繼而利用高精度動力學(xué)模型計(jì)算再入時刻。本文首先介紹了利用簡化通用攝動模型(SGP4)從TLE數(shù)據(jù)中獲取精度較高BC值的算法，并實(shí)現(xiàn)其計(jì)算過程與結(jié)果，建立BC值訓(xùn)練集，然后詳細(xì)闡述了LSTM預(yù)測模型結(jié)構(gòu)的再入時刻預(yù)報模型的建立過程，最后對訓(xùn)練好的模型性能進(jìn)行驗(yàn)證。使用LSTM模型剔除了BC值中存在的誤差因素，可以在空間碎片再入前中長期獲得比較準(zhǔn)確的再入時刻。不僅如此，這種方法在一次學(xué)習(xí)之后可以快速地對大批量空間碎片再入時刻進(jìn)行動態(tài)跟蹤預(yù)測，提高預(yù)測精度，并可大幅減少積分計(jì)算量，顯著提升運(yùn)算效率。

1 彈道系數(shù)的計(jì)算

1.1 TLE數(shù)據(jù)的預(yù)處理

TLE數(shù)據(jù)中包含了空間碎片的軌道根數(shù)、平運(yùn)動以及和大氣阻力相關(guān)的B*數(shù)據(jù)。在大批量使用TLE數(shù)據(jù)時，大概率會遇到空間碎片TLE數(shù)據(jù)異常情況，包括重復(fù)歷元、數(shù)據(jù)格式錯亂、編號異常等。因此需要對整個TLE數(shù)據(jù)進(jìn)行預(yù)處理，以去除異常和錯誤的數(shù)據(jù)[16]。預(yù)處理過程主要包括去除重復(fù)TLE數(shù)據(jù)、平均運(yùn)動項(xiàng)異常數(shù)據(jù)、傾角異常、B*為負(fù)值的數(shù)據(jù)等，由于TLE格式可讀性問題，預(yù)處理過后的TLE數(shù)據(jù)依然會存在某些異常情況，會明顯反應(yīng)在后續(xù)計(jì)算結(jié)果中。

1.2 彈道系數(shù)BC估計(jì)方法

以Saunders等[17]和Gondelach等[18]提出利用TLE數(shù)據(jù)修正BC值的方法作為本文研究的比對基準(zhǔn)。其基本原理是通過利用TLE配套的SGP4模型從TLE數(shù)據(jù)中獲得碎片的軌道狀態(tài)、物理參數(shù)(特別是B*數(shù)據(jù))等信息，建立連續(xù)TLE數(shù)據(jù)間的平運(yùn)動參數(shù)變率迭代估計(jì)任意時刻的BC值。

由于TLE數(shù)據(jù)中不包含短周期項(xiàng)的變化，數(shù)據(jù)中半長軸的變化可以假設(shè)為僅受大氣阻力影響下的長期項(xiàng)變化。考慮到太陽光壓導(dǎo)致的半長軸長周期變化可能包括在TLE數(shù)據(jù)中，但由于太陽光壓與低軌衛(wèi)星再入過程中大氣阻力攝動相比，其影響要小得多，因此本文忽略光壓的影響。由此可假設(shè)基于TLE數(shù)據(jù)平運(yùn)動參數(shù)計(jì)算得出的結(jié)果與由軌道外推計(jì)算結(jié)果之間的任何差異都是因?yàn)榇髿庾枇椀老禂?shù)BC值的不正確猜測引起的。利用連續(xù)TLE數(shù)據(jù)計(jì)算BC值的過程如下：

1)利用TLE數(shù)據(jù)中的平運(yùn)動參數(shù)計(jì)算兩個TLE數(shù)據(jù)的半長軸的變化量Δa。

(3)

Δa=a2-a1

(4)

式中：a為軌道半長徑；μ為地球引力常數(shù)；n0為空間物體軌道平均軌道運(yùn)動速度，也就是TLE的平運(yùn)動項(xiàng)。a1與a2分布表示前后兩次TLE數(shù)據(jù)的半長徑。

2)使用高精度的軌道動力學(xué)模型進(jìn)行外推，計(jì)算兩個TLE時刻(tTLE1,tTLE2)內(nèi)僅由大氣阻力引起的軌道半長軸的變化。其中p為軌道半通徑，θ為真近點(diǎn)角，e為偏心率，p=a(1-e2),fRdrag,fTdrag分別是由阻力在徑向和橫向分量產(chǎn)生的加速度。

(5)

(6)

3)利用TLE數(shù)據(jù)中的B*參數(shù)確定BC的初值。B*參數(shù)是SGP4模型中的近似反應(yīng)阻力的系數(shù)，可以通過(7)式計(jì)算得到BC的初值β1。第二個TLE根數(shù)的BC值β2可以通過第一個BC的值β1和Δaprop之間的線性關(guān)系計(jì)算得到。當(dāng)Δaprop小于10-4千米時，BC值收斂。如果在后續(xù)計(jì)算過程中出現(xiàn)了過早再入的情況，可以將β1的值設(shè)置為原來值的90%，但是需要更多次的迭代方能收斂。

(7)

4)對BC值進(jìn)行迭代計(jì)算，直到迭代收斂。βn是迭代過程中第n個BC值。

(8)

式中：Δadiff=Δa-Δaprop。

在默認(rèn)情況下，使用正向傳播來進(jìn)行BC值的估算，取最早TLE時刻的狀態(tài)，將它外推到最新TLE數(shù)據(jù)時刻的狀態(tài)。也可以使用反向傳播，從最新的TLE時刻的狀態(tài)向后傳播到前一個TLE時刻的狀態(tài)，反向傳播可以防止預(yù)測時空間碎片過早的再入。在臨近再入時刻估計(jì)時十分有用，因?yàn)椴粶?zhǔn)確的BC估計(jì)值可能導(dǎo)致再入預(yù)報異常。以上估計(jì)過程其本質(zhì)原理類似Newton-Raphson方法。如果TLE數(shù)據(jù)充分，且TLE數(shù)據(jù)在時序上連續(xù)且間隔不大，不失為一種較優(yōu)的方法，然后對于中短期預(yù)報來說，由于其數(shù)據(jù)使用率較高，如果其中出現(xiàn)異常點(diǎn)，其預(yù)測值會突然拉偏，導(dǎo)致預(yù)報誤差突然增大。空間物體的軌道是攝動情況下連續(xù)可微的擬二次曲線，因此有必要在預(yù)測過程中考慮短期的數(shù)據(jù)情況的變化，同時又能兼顧長期的動力學(xué)演化趨勢，因此本文參考了長短期記憶網(wǎng)絡(luò)預(yù)測模型LSTM在空間任務(wù)中的研究結(jié)論，將其應(yīng)用于彈道系數(shù)BC值的預(yù)測。

2 長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM模型原理

為了解決原有的循環(huán)神經(jīng)網(wǎng)絡(luò)(RNN)存在梯度爆炸或梯度消失的問題[19]，Hochreiter等[20]于1997年提出了長短期記憶神經(jīng)網(wǎng)絡(luò)LSTM，并被Graves等[21]進(jìn)行了改進(jìn)和推廣。LSTM神經(jīng)網(wǎng)絡(luò)經(jīng)過了特殊的設(shè)計(jì)，使得其適合在重要事件之間存在很長時間滯后時從經(jīng)驗(yàn)中學(xué)習(xí)，以對時間序列進(jìn)行分類、處理和預(yù)測，目前在多個領(lǐng)域得到了廣泛的應(yīng)用，其單元內(nèi)部包含遺忘門、輸入門和輸出門，具體結(jié)構(gòu)如圖1所示。

圖1 LSTM單元內(nèi)部結(jié)構(gòu)Fig.1 Internal structure of LSTM unit

遺忘門用來確定狀態(tài)傳遞的過程中上一時刻的單元狀態(tài)ct-1有多少保留到當(dāng)前時刻的狀態(tài)ct。讀取上一個時刻輸出值ht-1和當(dāng)前時刻輸入的BC值數(shù)據(jù)xt，并利用設(shè)定的連接權(quán)重矩陣Wf和sigmoid函數(shù)σ計(jì)算得到經(jīng)過處理的需要保留下來的歷史信息狀態(tài)，如式(9)所示，其中bf為遺忘門的偏移值。將會輸出一個0至1之間的數(shù)值給歷史狀態(tài)ct-1，1表示完全保留，0表示完全遺棄。

ft=σ(Wf[ht-1,xt]+bf)

(9)

輸入門利用連接權(quán)重參數(shù)Wi、偏移值bi等參數(shù)來確定當(dāng)前時刻輸入的BC值數(shù)據(jù)xt有多少可以作為輸入門的輸出it，如式(10)所示。再利用計(jì)算新狀態(tài)的連接權(quán)重矩陣WC、偏移值bC確定臨時的當(dāng)前時刻新狀態(tài)，如式(11)。再利用遺忘門和輸入門的數(shù)據(jù)更新狀態(tài)，將上一時刻狀態(tài)ct-1更新到當(dāng)前時刻狀態(tài)ct，如式(12)。

it=σ(Wi·[ht-1,xt]+bi)

(10)

(11)

(12)

輸出門利用連接權(quán)重矩陣Wo、偏移值bo確定當(dāng)前時刻的輸出值ht，如式(13)、式(14)所示。

ot=σ(Wo[ht-1,xt]+bo)

(13)

ht=ot·tanh(Ct)

(14)

3 再入預(yù)報流程與性能評估

3.1 再入時刻預(yù)測流程

首先建立預(yù)測模型的訓(xùn)練數(shù)據(jù)集，別選取近三年再入大氣層的80個和20個空間碎片的TLE數(shù)據(jù)作為預(yù)測模型的訓(xùn)練集和測試集。TLE數(shù)據(jù)和再入時刻都是從北美防空司令部空間目標(biāo)監(jiān)視編目軌道根數(shù)發(fā)布網(wǎng)站(https://www.space-track.org)上獲取的，該網(wǎng)站發(fā)布了空間物體歷史TLE數(shù)據(jù)以及再入日期，確切的再入時刻未給出，因此最終再入時刻僅能精確到天。

LSTM模型預(yù)測彈道系數(shù)BC值的具體流程如圖2所示。首先對所有TLE數(shù)據(jù)進(jìn)行預(yù)處理，剔除其中誤差較大及錯誤的數(shù)據(jù)，利用1.2節(jié)的方法計(jì)算出每一TLE時刻的BC值，隨后將得到的BC值隨時間變化數(shù)據(jù)代入LSTM模型，確定模型中的特征參數(shù)，由此得到了BC值的精確預(yù)測模型，再將測試集代入預(yù)測模型中，得到對未來BC值的預(yù)測，即可計(jì)算出精確的再入時刻，最后對再入時刻進(jìn)行評估。

圖2 空間碎片再入時刻預(yù)測流程Fig.2 The reentry time prediction scheme of space debris

3.2 彈道系數(shù)BC的估計(jì)

圖3為編號為32712的空間碎片再入大氣層前200天分別利用1.2節(jié)所述的迭代過程和由TLE數(shù)據(jù)中B*參數(shù)兩種方法計(jì)算得到的BC值隨時間的變化曲線。圖中橫軸為再入前的時間，縱軸為BC值的大小。圖示TLE數(shù)據(jù)是經(jīng)過預(yù)處理的，可以看出兩種方法計(jì)算的BC值具有相同的變化趨勢，但是有一個偏移量，同時BC值中仍然存在一些異常點(diǎn)，如果進(jìn)一步利用近地點(diǎn)半徑進(jìn)行濾波，這些異常值會被剔除。

圖3 編號32712空間碎片BC值與再入時刻關(guān)系Fig.3 The relationship between BC value and reentry time of the space debris No. 32712

利用1.2節(jié)的迭代方法，選取了100個已再入隕落的空間碎片作為目標(biāo)集，計(jì)算了目標(biāo)集對應(yīng)的從再入前200天到再入時刻彈道系數(shù)BC的最小值、平均值和最大值。計(jì)算結(jié)果表明基于迭代修正的BC值預(yù)測再入時刻精度在10%的置信度為60%～70%，預(yù)測精度在20%的置信度為95%。根據(jù)歐洲航天局的DISCOS數(shù)據(jù)庫提供的信息，彈道系數(shù)BC估計(jì)值精度20%的置信度為85%。考慮到大氣密度10～20%誤差因素，表明BC可以利用1.2節(jié)的流程建立訓(xùn)練集。因此以此為參考對象建立訓(xùn)練集樣本，利用LSTM模型預(yù)測BC值，進(jìn)而預(yù)報再入時刻，并評估預(yù)報精度與模型性能。

3.3 建立預(yù)測BC值的LSTM模型

LSTM模型的建立主要分3個步驟：首先需要對BC值數(shù)據(jù)集進(jìn)行整理；其次是完成對模型的構(gòu)建，包括確定用于訓(xùn)練和預(yù)測所需要的特征參數(shù)和確定LSTM神經(jīng)網(wǎng)絡(luò)模型的具體參數(shù)；最后是利用測試集數(shù)據(jù)代入得到預(yù)測結(jié)果，并對預(yù)測模型準(zhǔn)確性進(jìn)行驗(yàn)證,如圖4所示。

圖4 LSTM預(yù)測模型框架Fig.4 The framework of LSTM prediction model

采用多時刻輸入單一時刻輸出的模式對彈道系數(shù)BC值進(jìn)行預(yù)測，利用40天的TLE數(shù)據(jù)作為預(yù)測模型的輸入，預(yù)測結(jié)果作為輸出。設(shè)置兩層LSTM預(yù)測模型，并使用全連接層對LSTM層進(jìn)行線性化。訓(xùn)練過程中，反向傳播的誤差項(xiàng)是損失函數(shù)對輸出值的倒數(shù)，損失函數(shù)是均方根誤差函數(shù)，給定各層神經(jīng)元個數(shù)、學(xué)習(xí)率以及訓(xùn)練迭代步數(shù)，利用梯度下降算法更新網(wǎng)絡(luò)權(quán)重并獲得最終的隱藏層。隨后將測試集進(jìn)行迭代預(yù)測，并輸出預(yù)測結(jié)果。

建立彈道系數(shù)預(yù)測模型時需要考慮5個參數(shù)，即訓(xùn)練步長、輸出層維數(shù)、隱藏層層數(shù)、隱藏層維數(shù)和輸出變量維數(shù)。訓(xùn)練步長的設(shè)置要考慮該步長包含的信息量是否可以滿足模型訓(xùn)練，若設(shè)置的過小，則不足以使模型在訓(xùn)練中學(xué)習(xí)BC數(shù)據(jù)的變化規(guī)律，同時考慮到模型訓(xùn)練難易程度，設(shè)置過大可能使得數(shù)據(jù)序列中包含過多無用數(shù)據(jù)，導(dǎo)致預(yù)測性能下降并增加訓(xùn)練時間。本工作測模型中訓(xùn)練步長設(shè)置為20步。

針對于本文的隕落時間預(yù)測場景，輸入層維數(shù)的設(shè)置根據(jù)雙行軌道根數(shù)中軌道和動力學(xué)特征數(shù)據(jù)來確定,在本模型中使用TLE數(shù)據(jù),所以輸入層維數(shù)設(shè)置為1。多隱藏層的網(wǎng)絡(luò)從理論上來說有足夠的表達(dá)能力可以擬合任何復(fù)雜的函數(shù),但網(wǎng)絡(luò)結(jié)構(gòu)越復(fù)雜越難以訓(xùn)練,并且可能出現(xiàn)過擬合的情況,即在訓(xùn)練數(shù)據(jù)上表現(xiàn)很好,但泛化能力差。隱藏層維數(shù)即為隱藏層神經(jīng)元個數(shù),對于較復(fù)雜、非線性程度大的序列,設(shè)置的神經(jīng)元個數(shù)越多擬合效果越好,但花費(fèi)的時間也越多,同時可能出現(xiàn)過擬合的現(xiàn)象。本模型中隱藏層由兩層LSTM層和全連接層構(gòu)成,第二層的LSTM層神經(jīng)元個數(shù)是第一層個數(shù)的兩倍，第三層的全連接層神經(jīng)元個數(shù)和第一層LSTM層個數(shù)相等，由于本文主要是利用算法，算法具體細(xì)節(jié)不作贅述。

3.4 LSTM預(yù)測模型的性能評估

使用平均絕對百分比誤差(MAPE)作為評估參數(shù)與預(yù)測模型輸出值之間的預(yù)測準(zhǔn)確性度量，用以量化訓(xùn)練后預(yù)測模型的性能。MAPE定義為測試數(shù)據(jù)集中每個數(shù)據(jù)點(diǎn)上的BC預(yù)測值βpre與計(jì)算值βcal之間的真實(shí)誤差之和與計(jì)算值βcal之和的比值，下式中以ξ標(biāo)志，其數(shù)學(xué)表達(dá)式為：

(15)

式中：n是參與訓(xùn)練的有效數(shù)據(jù)數(shù)量。利用預(yù)測模型預(yù)測出結(jié)果越精確，式(15)將越逼近其下邊界0，MAPE的數(shù)值越低說明預(yù)測模型的性能越好。這里選擇的度量標(biāo)準(zhǔn)與預(yù)測模型的目標(biāo)一致，該目標(biāo)是通過最大程度減少訓(xùn)練數(shù)據(jù)中的總誤差，使所有數(shù)據(jù)點(diǎn)盡可能平滑。

3.5 LSTM模型預(yù)測BC值效能分析

采用編號20230的空間碎片，利用其再入前200天至160天之間40天的TLE數(shù)據(jù)作為預(yù)測模型的輸入，預(yù)測再入前160天的BC值隨時間的變化關(guān)系。

圖5橫軸表示再入預(yù)報所用TLE根數(shù)時間與www.space-track.org公布的再入時間的間隔，即當(dāng)前歷元距離再入日期的天數(shù)，從再入前160天至再入時刻為止，縱軸表示了BC值的大小。其中橙色折線的中心點(diǎn)表示了10天內(nèi)的BC計(jì)算值的平均值，藍(lán)色折線的中心點(diǎn)表示了10天內(nèi)的BC預(yù)測值的平均值，誤差線表示了10天內(nèi)的最大值和最小值。總數(shù)據(jù)統(tǒng)計(jì)預(yù)測模型的性能指標(biāo)MAPE達(dá)到13.9%，在某些時間段，MAPE更低，其中有一些未能預(yù)處理掉的異常數(shù)據(jù)點(diǎn)也清晰地反映在圖上,如果剔除一些異常點(diǎn)，預(yù)測結(jié)果的MAPE值會更低。數(shù)據(jù)結(jié)果表明訓(xùn)練好的預(yù)測模型可以較好地學(xué)習(xí)BC值變化趨勢，表明本文建立的LSTM預(yù)測模型是有效的。

圖5 編號20230的空間碎片再入前160天的BC預(yù)測值和計(jì)算值與預(yù)測時間的關(guān)系Fig.5 The relationship between BC predicted/calculated values and predicted times for 160 days prior to reentry of the space debris No.20230

3.6 再入預(yù)報能分析

BC值預(yù)測的目標(biāo)是為了更好地預(yù)報空間碎片再入時刻，或者是在軌壽命。在獲得BC預(yù)測值之后，基于最新的TLE根數(shù)，考慮地球非球形攝動、大氣阻力攝動、太陽光壓攝動和日月三體引力攝動，利用RKF7/8階數(shù)值積分進(jìn)行軌道外推，外推截止軌道高度設(shè)定80 km即認(rèn)為目標(biāo)已經(jīng)發(fā)生了再入，采集此時的歷元作為再入時刻。

為了評估再入時刻的預(yù)報精度，使用時間的相對誤差估計(jì)的精度，再入預(yù)報精度ε定義為實(shí)際預(yù)測時間誤差與最后所用軌道根數(shù)歷元距離隕落日期的時間之比值。

(16)

式中：tpre是預(yù)測的再入日期,tact是space-track.org網(wǎng)站上公布的再入日期,tTLElast是最后用于預(yù)測再入時刻的TLE軌道根數(shù)歷元。

以編號42727目標(biāo)為例，取再入前160天內(nèi)所有有效的TLE數(shù)據(jù)，利用當(dāng)天的TLE數(shù)據(jù)計(jì)算出衛(wèi)星的狀態(tài)，并預(yù)測出從TLE數(shù)據(jù)時刻到再入時刻的BC變化趨勢。利用衛(wèi)星軌道參數(shù)和BC變化計(jì)算出衛(wèi)星再入時刻并與space-track.org網(wǎng)站給出的真實(shí)時刻計(jì)算再入時刻誤差，如圖6所示。可以發(fā)現(xiàn)，利用本文方法得到平均預(yù)報誤差為7.3%，再入前160天時預(yù)測誤差在4%以下，再入前30天時預(yù)測誤差在10%以下。但是再入前20天預(yù)測的效果沒有長期預(yù)測的理想，最大預(yù)測誤差可達(dá)到51%，這是因?yàn)閟pace-track.org給出的目標(biāo)再入時刻只精確到天，依據(jù)式(16)的標(biāo)準(zhǔn)，越到末端誤差會急劇增大，導(dǎo)致即使預(yù)測是準(zhǔn)確的，再入前10天預(yù)測時每差24小時就帶有近10%的誤差，再入最后2天每一至兩個小時就接近10%誤差。同時，短期預(yù)測誤差率大也可能是因?yàn)榻咏偃氪髿鈱訒r氣動力快速變化導(dǎo)致無法準(zhǔn)確預(yù)測與計(jì)算，因此后續(xù)評估主要對大于30天臨近再入時間的再入預(yù)報有較好的指導(dǎo)意義，但本文同時也給出了小于30天的預(yù)測結(jié)果。

圖6 編號42727空間碎片再入時刻預(yù)報結(jié)果Fig.6 Re-entry time forecast results of the space debris No.42727

選取5個已再入物體作為測試目標(biāo)，利用本文的方法預(yù)測其再入前160天至再入時刻的BC變化，基于再入前160天，90天和30天的TLE數(shù)據(jù)計(jì)算再入時刻，并計(jì)算再入時刻誤差。計(jì)算結(jié)果如表1所示。

表1 測試目標(biāo)再入時刻預(yù)報結(jié)果Table 1 Re-entry time forecast results of the test targets

對選取的90個目標(biāo)于再入前10天、20天、30天、60天、90天和160天時預(yù)測再入時刻并得到累計(jì)誤差。所得到的再入預(yù)報時間累積誤差和95%置信區(qū)間結(jié)果如圖7所示。圖7表明，再入前160天時的預(yù)測最準(zhǔn)確，總誤差低于5%，再入前10、20天預(yù)測的結(jié)果不如整體結(jié)果準(zhǔn)確是由于公布的再入時刻精確到天導(dǎo)致。總體來看，再入前160天，90天，30天分別得到的預(yù)測誤差在5%以下，8%以下和12%以下，算例表明利用長短期記憶神經(jīng)網(wǎng)絡(luò)先預(yù)測出物體BC變化規(guī)律再進(jìn)行再入時刻預(yù)報可以得到較好的空間碎片無控再入預(yù)報時刻。

圖7 90個目標(biāo)的再入前10、20、30、60、90和160天的再入時刻預(yù)報誤差的累積分布Fig.7 Cumulative distributions of reentry time prediction errors of 90 objects for 10, 20, 30, 60, 90, and 180 days prior to reentry

4 結(jié) 論

空間碎片無控再入時刻預(yù)報是一項(xiàng)國際上非常關(guān)注的重要工作，其受到大氣密度、軌道數(shù)據(jù)精度、彈道系數(shù)估計(jì)等諸多因素的制約，以至于準(zhǔn)確預(yù)測空間碎片的再入時刻是十分困難的。在已有大氣阻力彈道系數(shù)BC值預(yù)測方法的基礎(chǔ)上，建立了基于長短期記憶神經(jīng)網(wǎng)絡(luò)深度學(xué)習(xí)LSTM模型預(yù)測BC值，進(jìn)而預(yù)測空間碎片再入時刻的方法，獲得了較好的研究結(jié)果。從選取的樣本訓(xùn)練模型以及測試結(jié)果看，在95%的置信度內(nèi)，90天預(yù)測精度達(dá)到8%，30天預(yù)測精度達(dá)到12%，考慮到大氣密度模型帶來的不確定影響，結(jié)果表明LSTM模型預(yù)測方法在空間碎片再入預(yù)報具有很好應(yīng)用價值。后續(xù)我們將進(jìn)一步精細(xì)化考慮空間碎片類型，將相近類型的空間碎片進(jìn)行聚類，建立更精確數(shù)據(jù)樣本，提高再入或者在軌壽命預(yù)測精度。