基于TLE的彈道系數(shù)計算方法及應(yīng)用分析

張煒，崔文，張育衛(wèi)，劉興，朱俊

1. 西安衛(wèi)星測控中心，西安 710600 2. 宇航動力學(xué)國家重點實驗室，西安 710043

大氣阻力是低軌道空間目標(biāo)所受的最主要非保守攝動力，其大小主要由兩個因素決定：大氣密度和空間目標(biāo)的彈道系數(shù)。彈道系數(shù)B定義為大氣阻力系數(shù)CD與面質(zhì)比的乘積[1]，即B=CD×A/m，其中A為迎風(fēng)面積，m為質(zhì)量。

準(zhǔn)確計算低軌目標(biāo)的彈道系數(shù)對于開展大氣密度研究、低軌目標(biāo)軌道確定、再入預(yù)報等具有重要意義[2]。通常彈道系數(shù)被用作待估參數(shù)，與狀態(tài)矢量一起，在基于測軌數(shù)據(jù)的軌道確定過程中求解[3]。除實際測量數(shù)據(jù)外，美國戰(zhàn)略司令部于space-track網(wǎng)站發(fā)布的兩行要素(Two Line Element，TLE)也是重要數(shù)據(jù)源，具有易獲取、目標(biāo)數(shù)量大、時間及空間覆蓋廣等特點，有獨特的研究應(yīng)用價值。任廷領(lǐng)等利用TLE反演熱層密度并進行大氣模型誤差分析，反演過程中需要獲得低軌目標(biāo)的彈道系數(shù)“真值”[4-5]；Pardini等在機構(gòu)間空間碎片協(xié)調(diào)委員會(Inter-Agency Space Debris Coordination Committee，IADC)再入聯(lián)測中利用TLE進行再入預(yù)報，彈道系數(shù)計算也是其中重要一環(huán)[6-7]。

TLE中并沒有彈道系數(shù)的直接信息，較為相關(guān)的是B*項，B*與真實的彈道系數(shù)B可近似為：B=12.741621B*[8]，但是此結(jié)果仍不夠準(zhǔn)確，與真實值之間存在不小的偏差[9]，獲得較為準(zhǔn)確的彈道系數(shù)仍需要運用其他方法。基于TLE的彈道系數(shù)計算方法主要可以分為兩類：一類基于TLE生成偽測量數(shù)據(jù)，再利用偽測量數(shù)據(jù)重新確定軌道并解算彈道系數(shù)[10-12]，這類方法在實現(xiàn)上與基于實測數(shù)據(jù)的軌道確定過程并無本質(zhì)上的區(qū)別；第二類通過比較TLE之間的軌道衰減計算彈道系數(shù)，大部分研究基于兩組TLE進行計算[13-16]，彈道系數(shù)結(jié)果受TLE精度影響大，結(jié)果穩(wěn)定性較差，實際應(yīng)用中往往需要多次計算取結(jié)果均值。

本文主要研究第二類計算方法，即基于TLE軌道衰減計算彈道系數(shù)，考慮到TLE偏心率精度比半長軸精度低很多[17]，計算時僅考慮半長軸的衰減。介紹一種常見的基于兩組TLE的直接計算法，分析了不同TLE間隔對結(jié)果精度的影響；提出一種基于多組TLE的迭代計算法，并從彈道系數(shù)計算結(jié)果、在再入預(yù)報中的應(yīng)用等方面對這兩種方法進行了比較，為實際應(yīng)用提供依據(jù)。

1 彈道系數(shù)計算方法

1.1 基于兩組TLE的直接計算法

大氣阻力僅有長期作用效果，對于低軌目標(biāo)其他攝動力的長期項與大氣阻力相比基本可以忽略，因此可以通過平半長軸的衰減計算彈道系數(shù)。兩組TLE之間平半長軸的變化認(rèn)為由大氣阻力引起[18]，大氣阻力引起的半長軸變化可以表示為：

(1)

式中：a為半長軸；e為偏心率；n為平運動速度；p為軌道半通徑；f為真近點角；r為徑向距離；S和T分別為大氣阻力的徑向和橫向分量，

(2)

式中：v為空間目標(biāo)相對于大氣的運動速度；ωE為地球轉(zhuǎn)動速率；μ為地球引力常數(shù)；i為軌道傾角；ρ為大氣密度，本文使用MSIS-90模型計算。將式(2)代入式(1)，簡化后得到：

則兩組TLE之間半長軸的衰減量ΔaTLE可以表示為：

因此，彈道系數(shù)計算：

(3)

其中：

Sang等給出的計算公式略有差別[14]，但計算原理基本一致。

1.2 基于多組TLE的彈道系數(shù)計算方法

直接計算方法的問題是彈道系數(shù)結(jié)果受TLE精度影響大，結(jié)果穩(wěn)定性較差，實際應(yīng)用中往往需要多次計算取結(jié)果均值。因此，提出使用最小二乘法對N(N≥3)組TLE進行彈道系數(shù)擬合，使半長軸的計算值與實測值之間的平均誤差最小。假設(shè)已有N組連續(xù)的TLE，彈道系數(shù)計算過程如下：

2)設(shè)置彈道系數(shù)初值為B(0)=12.741 621B*，其中B*為第一組TLE的大氣阻力項。

3)從第一組TLE根數(shù)開始對Gauss型攝動方程進行軌道積分，考慮的攝動力包括地球非球形J2項、J3項和大氣阻力。

4)若積分結(jié)束時刻與第k組TLE歷元時刻之差小于積分步長的一半，則計算半長軸計算值與實測值之差Δak及半長軸對彈道系數(shù)的偏導(dǎo)數(shù)?ak/?B，否則繼續(xù)下一步積分。

5)第N條根數(shù)計算結(jié)束后，計算彈道系數(shù)修正值：

(4)

6)若ΔB(j)小于收斂閾值，則計算結(jié)束，否則對彈道系數(shù)進行修正：

B(j+1)=B(j)+ΔB(j)

返回步驟3)開始第j+1次迭代。

由于計算過程采用平根數(shù)作為軌道形式，不含周期項，因此軌道積分可采用半分析方法，步長取為軌道周期的整數(shù)倍，以提高計算效率。此外，計算過程中可以根據(jù)半長軸的殘差情況識別誤差較大的TLE并重新編輯數(shù)據(jù)資料，進而提高結(jié)果精度。

2 結(jié)果與分析

2.1 彈道系數(shù)結(jié)果

選擇NORAD編號為00694和42821的兩個目標(biāo)作為算例目標(biāo)進行驗證。其中，00694是一個外形為圓柱體的火箭箭體，近地點軌道高度約463 km，TLE選用的時間區(qū)間為2014年2月—2015年1月，期間F10.7指數(shù)90～200，太陽活動整體處于中等水平；42821是一個尺寸約為10 cm的立方體衛(wèi)星，質(zhì)量約1 kg，目前已經(jīng)隕落，TLE選用的時間區(qū)間為2018年5月—2019年5月，期間F10.7指數(shù)65～80，太陽活動水平極低。兩個目標(biāo)的基本信息如表1所示。

表1 兩個算例目標(biāo)的基本信息

圖1和圖2所示分別是00694和42821目標(biāo)的彈道系數(shù)序列，表2列出了兩個目標(biāo)彈道系數(shù)的統(tǒng)計結(jié)果。使用迭代法計算時，N取值為6，即基于6組TLE計算彈道系數(shù)，使用直接法時TLE選擇的是該TLE序列的第一組和最后一組，使兩種方法所用TLE的時間跨度一致。根據(jù)文獻[14]的分析結(jié)果，1980—2011年間00694的彈道系數(shù)均值為0.014 81 m2/kg，與本文結(jié)果僅差約2%；42821的理論彈道系數(shù)約為0.021 m2/kg，與本文結(jié)果也極為接近，證明了兩種方法的準(zhǔn)確性。對比兩種方法的彈道系數(shù)均值，發(fā)現(xiàn)差異極小。基于兩組TLE的彈道系數(shù)均值分別為0.014 48 m2/kg和0.019 81 m2/kg，基于多組TLE的彈道系數(shù)均值分別為0.014 49 m2/kg和0.019 82 m2/kg。認(rèn)為長期效果方面兩種方法并無差別，且兩種方法的彈道系數(shù)結(jié)果在變化趨勢上基本一致。但是基于多組TLE的彈道系數(shù)結(jié)果穩(wěn)定性略好，其原因是基于兩組TLE計算時TLE誤差直接被彈道系數(shù)吸收，彈道系數(shù)結(jié)果受TLE半長軸精度影響較為明顯。較為典型的是圖1和圖2中方框標(biāo)注的結(jié)果，明顯誤差偏大。需要說明的是，兩個算例目標(biāo)所用TLE均為根據(jù)平運動平滑過濾后的結(jié)果，已經(jīng)剔除大部分誤差較大TLE，若使用原始數(shù)據(jù)進行計算，基于多組TLE迭代計算的優(yōu)勢將更明顯。

圖1 00694的彈道系數(shù)結(jié)果Fig.1 Ballistic coefficient results of 00694

圖2 42821的彈道系數(shù)結(jié)果Fig.2 Ballistic coefficient results of 42821

m2/kg

2.2 用于再入預(yù)報的比較

確定空間目標(biāo)的彈道系數(shù)是再入預(yù)報的關(guān)鍵步驟，彈道系數(shù)結(jié)果的準(zhǔn)確性直接影響再入預(yù)報結(jié)果的精度。本節(jié)將討論兩種彈道系數(shù)計算方法用于再入預(yù)報時的效果。進行再入預(yù)報時，基于最新兩組或多組(取值為6)TLE計算彈道系數(shù)，并使用SGP4/SDP4模型直接計算最后一組TLE在其歷元時刻的狀態(tài)矢量，作為再入預(yù)報的初始軌道。最后基于初始軌道和彈道系數(shù)結(jié)果，采用數(shù)值積分的方法進行軌道外推直到軌道高度低于80 km，積分結(jié)束時刻即為再入時刻?？紤]的攝動項包括地球非球形引力、大氣阻力攝動、太陽光壓攝動和日月引力攝動。

使用再入預(yù)報相對誤差δ評估再入預(yù)報的精度，再入預(yù)報相對誤差δ為:

式中：treal為再入目標(biāo)的實際再入時間；tpred為預(yù)報的再入時間；t0為進行再入預(yù)報時使用軌道的歷元，本文中即為兩組或多組TLE中最后一組TLE的歷元。

以42821目標(biāo)為例，取最后一次計算結(jié)果為實際再入時間，利用兩種方法的彈道系數(shù)結(jié)果進行再入預(yù)報的精度情況如圖3所示?？梢钥吹?，基于兩組TLE的最大預(yù)報誤差為53.2%，平均誤差約12.2%，而基于多組TLE的最大預(yù)報誤差為39.8%，平均誤差約10.6%?？傮w而言，無論是預(yù)報精度還是預(yù)報結(jié)果的穩(wěn)定性，基于多組TLE的再入預(yù)報方法均有較大優(yōu)勢，原因就是基于多組TLE的彈道系數(shù)結(jié)果具有更好的穩(wěn)定性。

圖3 42821目標(biāo)的再入預(yù)報結(jié)果Fig.3 Re-entry prediction results of 42821

但是基于多組TLE的彈道系數(shù)結(jié)果用于再入預(yù)報時并非一直最優(yōu)。受方法本身特點的限制，使用多組TLE時數(shù)據(jù)區(qū)間不可能太小，而兩組TLE的選擇卻靈活許多?？疾炝穗S機28個再入目標(biāo)再入前10天起兩種方法的精度情況，結(jié)果如圖4所示。其中，兩組或多組TLE均為計算時獲取的最新TLE結(jié)果，即使用直接法時TLE選用的是迭代法所用TLE序列中的最后兩組。黑色部分為基于兩組TLE預(yù)報精度較優(yōu)的目標(biāo)個數(shù)，紅色部分為基于多組TLE預(yù)報精度較優(yōu)的目標(biāo)個數(shù)，比較時各再入目標(biāo)的實際再入時間均取自space-track網(wǎng)站發(fā)布的最終再入信息?？梢钥吹?，當(dāng)距離目標(biāo)再入不足1天，基于兩組TLE的彈道系數(shù)結(jié)果較優(yōu)的次數(shù)明顯多于基于多組TLE；當(dāng)距離目標(biāo)再入大于1天，兩種方法的優(yōu)劣次數(shù)基本相當(dāng)。其原因是使用直接法計算時，由于選擇的是最新的兩組TLE，數(shù)據(jù)區(qū)間較短，彈道數(shù)據(jù)結(jié)果與數(shù)據(jù)弧段內(nèi)再入目標(biāo)的軌道數(shù)據(jù)(即TLE)情況、空間環(huán)境情況等形成良好的耦合系統(tǒng)，可以更好地反映再入目標(biāo)短期軌道衰減特性。另外，在計算過程中發(fā)現(xiàn)，當(dāng)臨近再入，基于多組TLE擬合彈道系數(shù)時迭代次數(shù)明顯增加，甚至可能出現(xiàn)不收斂的情況。

圖4 28個再入目標(biāo)再入前10天起兩種方法的精度情況Fig.4 Re-entry prediction results for 28 test objects in the last 10 days

綜合以上分析，在中長期再入預(yù)報中，建議采用基于多組TLE的彈道系數(shù)計算方法，以獲得更穩(wěn)定的再入預(yù)報結(jié)果；當(dāng)距離目標(biāo)再入不足1天，建議采用基于兩組TLE的彈道系數(shù)計算方法，以獲得更高精度的再入預(yù)報結(jié)果。

2.3 TLE歷元間隔對結(jié)果的影響

假設(shè)兩組TLE之間半長軸衰減的誤差為Δaε，實際衰減值為Δat，根據(jù)式(3)有：

式中：Bc為彈道系數(shù)的計算值；Bt為真實值。顯然，TLE間隔越大或運行高度越低，大氣阻力的累計效應(yīng)λ越大，同等半長軸誤差引起的彈道系數(shù)誤差越小。

同樣以NORAD編號為00694和42821的兩個目標(biāo)為例，驗證使用基于兩組TLE直接計算彈道系數(shù)時TLE間隔對結(jié)果的影響。兩個目標(biāo)均選擇2018年的TLE作為數(shù)據(jù)源，使分析過程二者經(jīng)歷的空間環(huán)境一致。00694目標(biāo)2018年的彈道系數(shù)均值較第2.1小節(jié)分析結(jié)果有較大變化，僅為0.009 6 m2/kg。一方面是因為太陽活動水平極低情況下大氣密度模型的誤差增大[19]，另一方面空間目標(biāo)的大氣阻力系數(shù)較太陽活動高年也有所降低[20]。42821目標(biāo)的彈道系數(shù)均值與第2.1小節(jié)一致，取為0.019 8 m2/kg。TLE間隔分別取為0.25 d、0.5 d和1 d。兩個目標(biāo)彈道系數(shù)的相對誤差分別如圖5、圖6所示。

圖5 00694不同TLE間隔的彈道系數(shù)結(jié)果Fig.5 Ballistic coefficient results of 00694 based on different TLE intervals

圖6 42821不同TLE間隔的彈道系數(shù)結(jié)果Fig.6 Ballistic coefficient results of 42821 based on different TLE intervals

對比0.25 d、0.5 d和1 d三種間隔，00694目標(biāo)的最大誤差分別為290%，223%和123%，標(biāo)準(zhǔn)差分別為69%，49%和34%；42821目標(biāo)的最大誤差分別為85%，56%和56%，標(biāo)準(zhǔn)差分別為20%，16%和15%。可以看到：1)TLE間隔越小，彈道系數(shù)誤差越大，結(jié)果精度的穩(wěn)定性也更低；2)隨著軌道高度降低，相同間隔TLE的彈道系數(shù)誤差減小(00694的平均高度約為900 km，42 821的平均高度約為370 km)；3)TLE間隔的增大對42821彈道系數(shù)結(jié)果的改進效果明顯較弱，TLE間隔為0.5 d和1 d兩種情況的最大誤差和標(biāo)準(zhǔn)差已經(jīng)極為接近。其原因與前文分析一致，即在相同時間間隔內(nèi)，軌道高度越低的目標(biāo)所受的大氣阻力累積效應(yīng)越大，同等半長軸誤差引起的彈道系數(shù)誤差越小。因此，若使用基于兩組TLE直接計算彈道系數(shù)，應(yīng)根據(jù)空間目標(biāo)的軌道情況確定TLE間隔，當(dāng)空間目標(biāo)的軌道高度較高，兩組TLE之間的間隔應(yīng)偏大，當(dāng)空間目標(biāo)的軌道高度較低，TLE的間隔可以相應(yīng)縮小。

3 結(jié)束語

本文研究了兩種基于TLE的低軌空間目標(biāo)彈道系數(shù)計算方法。介紹了一種常見的基于兩組TLE的直接計算法，提出了一種基于多組TLE的迭代計算法，并對兩種方法的應(yīng)用效果進行了分析。得出以下結(jié)論：1)從長期效果看，兩種方法彈道系數(shù)結(jié)果的均值和變化趨勢基本一致，但是基于多組TLE的彈道系數(shù)計算方法具有更好的穩(wěn)定性。2)用于再入預(yù)報時，在中長期預(yù)報中建議采用基于多組TLE的彈道系數(shù)計算方法，以獲得更穩(wěn)定的再入預(yù)報結(jié)果；當(dāng)距離目標(biāo)再入不足1 d，建議采用基于兩組TLE的彈道系數(shù)計算方法，以獲得更高精度的再入預(yù)報結(jié)果。3)若基于兩組TLE計算彈道系數(shù)，當(dāng)空間目標(biāo)的軌道高度較高，兩組TLE之間的間隔應(yīng)偏大，當(dāng)空間目標(biāo)的軌道高度較低，TLE的間隔可以相應(yīng)縮小。

本文的研究結(jié)果可為實際應(yīng)用中彈道系數(shù)計算方法和相關(guān)參數(shù)的選擇提供依據(jù)。后續(xù)可在使用最小二乘法迭代計算彈道系數(shù)時如何改進軌道結(jié)果，及軌道改進后對軌道預(yù)報精度的影響等方面開展研究。