基于TLE的彈道系數計算方法及應用分析

張煒，崔文，張育衛，劉興，朱俊

1. 西安衛星測控中心，西安 710600 2. 宇航動力學國家重點實驗室，西安 710043

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

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

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

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

1 彈道系數計算方法

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

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

(1)

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

(2)

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

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

因此，彈道系數計算：

(3)

其中：

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

1.2 基于多組TLE的彈道系數計算方法

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

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

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

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

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

(4)

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

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

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

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

2 結果與分析

2.1 彈道系數結果

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

表1 兩個算例目標的基本信息

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

圖1 00694的彈道系數結果Fig.1 Ballistic coefficient results of 00694

圖2 42821的彈道系數結果Fig.2 Ballistic coefficient results of 42821

m2/kg

2.2 用于再入預報的比較

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

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

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

以42821目標為例，取最后一次計算結果為實際再入時間，利用兩種方法的彈道系數結果進行再入預報的精度情況如圖3所示。可以看到，基于兩組TLE的最大預報誤差為53.2%，平均誤差約12.2%，而基于多組TLE的最大預報誤差為39.8%，平均誤差約10.6%。總體而言，無論是預報精度還是預報結果的穩定性，基于多組TLE的再入預報方法均有較大優勢，原因就是基于多組TLE的彈道系數結果具有更好的穩定性。

圖3 42821目標的再入預報結果Fig.3 Re-entry prediction results of 42821

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

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

綜合以上分析，在中長期再入預報中，建議采用基于多組TLE的彈道系數計算方法，以獲得更穩定的再入預報結果；當距離目標再入不足1天，建議采用基于兩組TLE的彈道系數計算方法，以獲得更高精度的再入預報結果。

2.3 TLE歷元間隔對結果的影響

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

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

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

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

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

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

3 結束語

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

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