非合作目標彈道系數解算研究及應用

劉舒蒔，李勰，*，滿海鈞，陳光明，曹建峰，鞠冰

(1.北京航天飛行控制中心，北京 100094；2.航天飛行動力學技術重點實驗室，北京 100094)

非合作目標泛指不能獲得有效信息的空間目標，如來襲導彈、敵機、失效或故障航天器、空間碎片等。自1957 年第一顆人造衛星上天以來，已有近24 000 個目標隕落墜入地球大氣層，其中大部分是非合作目標。再入過程中，雖然大多數航天器在與大氣相互作用過程中燒蝕、熔融和解體，但仍有10%～40%的器件落到地面，對人類安全和財產帶來潛在威脅[1-3]。隨著中國航天事業的不斷發展，服役期滿航天器再入大氣的頻率越來越高，對這些目標開展隕落預報，體現了中國作為航天大國負責任的態度。另外，掌握他國在軌空間目標的運動特性，尤其是非合作目標，對于國土防空早期預警及增強空間態勢感知能力，也具有十分重要的意義。

影響航天器再入隕落預報精度的主要因素包括軌道初始狀態、彈道系數及大氣密度。其中，彈道系數是計算大氣阻力的重要參數，包含了航天器質量、外形和氣動特性的組合參數，對于非合作目標，需要利用其軌道信息來解算彈道系數。

有學者利用卡爾曼濾波算法將航天器運動狀態和彈道系數進行聯合濾波估計，并對傳統卡爾曼濾波進行了改進。例如，采用不敏卡爾曼濾波算法對非線性函數的概率分布進行近似[4-5]，基于期望最大化迭代優化框架；文獻[6]利用精密定軌過程中的殘差信息，并結合空間環境變化對大氣阻力系數修正。這些方法需要豐富的軌道觀測數據，不適用于觀測條件不理想的非合作目標。

美國空間監視網對大部分在軌目標進行監視、編目和維護，北美防空司令部將目標的軌道信息以雙行根數（two line element，TLE）形式發布在Space-Track 網站。TLE 由于其完備性、實時性、精確性及開放性而備受關注，為非合作目標提供了有效的軌道信息[7-8]，是解算非合作目標彈道系數的重要數據源。Gupta[9]和Shoemaker[10]等根據ti時刻TLE根數推算t時刻近地點hpi和 遠地點高度hai，再根據t時刻TLE 根數計算t時刻近地點高度hpt和遠地點高度hat，令 |hpi?hpt|和 |hai?hat|最小來獲得最優彈道系數。據報道，該方法在航天器隕落的末期較準確，但是過于依賴TLE 中的偏心率，而偏心率的精度往往較低。為了彌補TLE 中偏心率精度較低的缺陷，Sharma 等[11]提出了一種同時估算彈道系數和初始偏心率的方法，使用響應曲面法擬合遠地點高度來估計偏心率和彈道系數。由于TLE 中的半長軸精度優于偏心率，Sang 等[12]利用不同TLE 中半長軸的衰減估算彈道系數，Saunders 等[13-14]對其進行了改進，通過使用多天TLE，多次迭代提高彈道系數估算精度。但對于非合作目標，由于不知道其軌道機動情況，在使用多天數據解算時，如果發生軌道機動，則會降低彈道系數估算精度。

針對上述問題，本文在Saunders 方法的基礎上進行了改進。首先，闡述了利用TLE 進行彈道系數估計的理論方法，分析了構造軌道半長軸衰減觀測量的依據；然后，針對非合作目標特點，對TLE 中的半長軸做了平滑檢測和二次判別，識別出野值、軌道機動和地磁暴3 種情況，比較了不同觀測弧長對彈道系數解算的影響；最后，利用本文算法對天宮一號正常在軌和隕落期間解算的彈道系數進行了精度評估，并在天宮一號隕落預報中進行了驗證。

1 數據與方法

1.1 基于TLE 的彈道系數估計

彈道系數B表達式[15-16]為

式中：m為航天器質量；CD為大氣阻力系數；s為航天器迎風面積。

大氣阻力對半長軸的影響主要體現為半長軸的長期衰減效應，即

式中：da/dt為半長軸的衰減率；v為航天器瞬時速度；μ為引力常數；r¨drag為阻力攝動加速度。

由TLE根數計算的半長軸記為atle，由式（3）計算：

式中：nk為TLE 根數中對衛星平均運動的描述量。

對于任意2 個連續時間間隔范圍的半長軸衰減，可記為

將?atle作為半長軸衰減觀測量，構建條件方程。半長軸的衰減主要受大氣阻力的影響，因此在給定時間間隔內，對式（2）進行積分，即可獲取半長軸衰減的理論解算結果?acomp，即

式中：ρ為大氣密度。

積分區間即為2 個TLE 的時間間隔。對于式（6），可利用常規求積算法給出。需要說明的是，在積分求和過程中，給定了積分步長，則積分區間內每個步點上的 ρ、a和v，可利用數值積分計算。數值積分的初值包含初始歷元時刻、初始位置速度和彈道系數，可通過SGP4 模型對TLE 進行單點轉換得到。將半長軸的衰減展開，并略去高階項，得到

式中：?B為彈道系數改進量。

由于式（7）為超定方程，采用最小二乘迭代求解，設定相鄰2 次迭代彈道系數改進量小于給定門限，方程收斂，從而得到最終的彈道系數。

計算流程如圖1 所示。

圖1 彈道系數解算流程Fig.1 Flowchart of ballistic coefficient solution

1.2 軌道半長軸衰減分析

TLE 是用特定方法平滑掉了周期擾動項的平均軌道根數，對應的軌道模型是SGP4/SDP4 解析模型[17]。TLE 數據沒有包含軌道精度信息，因此，利用天宮一號正常工作時段的精密星歷對TLE 的預報星歷進行了精度評估，從而確定構造半長軸衰減觀測量的依據。

選取天宮一號2012 年1 月1 日至1 月10 日期間的TLE 預報星歷和精密星歷進行比對，由于大氣阻力主要體現在高度變化率上，對TLE 軌道和精密軌道的高度變化率進行了比較分析。首先，去除瞬時精密軌道根數中的周期性變化（如地球重力場非球形的影響），得到每日的平根數；其次，根據TLE中的軌道傾角、偏心率和平運動周期計算半長軸；然后，利用樣條插值將TLE 的半長軸插值到與精密軌道根數相同的歷元；最后，依次計算精密軌道平根數半長軸和TLE 半長軸的每日衰減，并標示出兩者高度變化率的偏差，圖2 和圖3 分別展示了兩者每日和每2 日高度衰減情況。結果表明，TLE 高度變化率偏差的標準差為24.95 m/d 和30.6 m/(2d)，軌道高度衰減平均為211 m/d 和422 m/(2d)，TLE 反演出的每日軌道高度衰減誤差為11.85%，每2 日軌道高度衰減誤差為7.25%。由于常用大氣模式精度約10%，利用TLE 的半長軸每2 日衰減量構造觀測量，其精度能夠較好地反映大氣阻力的變化。

圖2 TLE 星歷與精密星歷每日高度衰減對比Fig.2 Comparison of daily height attenuation between TLE and precision ephemeris

圖3 TLE 星歷與精密星歷每2 日高度衰減對比Fig.3 Comparison of every two days height attenuation between TLE and precision ephemeris

1.3 TLE 數據預處理

由于飛行器在軌飛行過程中可能會發生軌道機動，造成半長軸變化，需要在預處理中檢測軌道機動，對軌道機動前后的TLE 分別處理，即式（5）積分區間不能包含軌道機動。此外，為提高彈道系數計算精度，應剔除精度較低的TLE，在這個過程中需結合地磁指數變化進行識別，避免地磁暴引起的半長軸快速衰減被誤判為低精度TLE。采用二次多項式擬合的方法對TLE 的半長軸進行初步檢測和二次判別，具體方法如下：首先，利用SGP4 模型計算出所有TLE 的半長軸；其次，對連續的第1,2,3,5,6,7 個歷元的半長軸進行二次多項式擬合，根據擬合的多項式計算第4 個歷元的半長軸的擬合值a4p；然后，計算第4 個歷元的擬合值與觀測值的差異ε=a4p?a4t，并根據3 倍方差準則進行判別；最后，對于超過門限的 ε進行二次判別，結合地磁指數，識別出3 類情況，即野值、軌道機動和地磁暴。

對2012 年1 月至3 月的天宮一號TLE 半長軸進行初步平滑檢測，門限值3 倍方差取20 m，如圖4所示，發現有5 個時段（紅框標識）的半長軸超過門限，圖5 給出了同時期地磁指數ap（每3 h 均值），并用紅框標出了超出門限的時段?？梢钥闯?，P1、P2、P3 和P4 時段都發生了地磁暴，導致P1、P2 和P4 時段航天器半長軸比平靜期衰減得多，即a4t比擬合值a4p小，使得 ε超出門限，因此P1、P2 和P4 時段的TLE 不是野值，不應剔除；P3 時段雖然有微小地磁擾動，但a4t反 而大于擬合值a4p，可以認定P3 是野值；P5 時段平滑值與實測值差異達到100 m 以上，即使這個時段發生了地磁擾動，也應認定為發生了軌道機動。

圖5 地磁指數變化Fig.5 Variations of geomagnetic index

2 結果與分析

2.1 彈道系數精度比較與評估

為了評估TLE 解算彈道系數的精度，本節比較了不同方法解算彈道系數的差異。在精密定軌中，往往把大氣阻力系數CD作為未知量與衛星的運動狀態矢量一起解算，再利用式（1）求解彈道系數。定軌中盡可能采用高精度的動力學模型，包括JGM64×64 地球重力場模型、MSISE2000 大氣密度模型，全面考慮了日月引力、太陽光壓和地球固體潮汐等多種攝動力。選擇天宮一號2011 年11 月19 日至2012 年12 月31 日期間的軌道資料進行比對，這段時期軌道資料豐富，飛行器性能穩定，地磁活動和太陽輻射處于中等活躍水平。利用每天的軌道觀測數據進行軌道確定，并解算當日CD，結果如圖6 所示，CD平均值為2.221 9，彈道系數平均值為0.009 4。同時，圖6（b）給出了期間的太陽輻射指數F10.7 和地磁指數AP（日均值），可以看出，彈道系數與F10.7 有較強的相關性，尤其在27 d 周期變化方面有很強的一致性。飛行器質量是精確已知的，且隨著發動機每次開機略有降低，飛行器姿態穩定的工況下每圈次的平均迎風面積也是穩定的，因此彈道系數的周期性變化是由CD變化引起的。大氣阻力系數CD是描述大氣粒子與衛星表面相互作用的無量綱系數，反映了粒子與衛星表面碰撞過程中的動量傳遞和能量耗散效率，其大小取決于衛星表面的材質、大氣的成分和溫度，以及入射大氣粒子與衛星表面的夾角等因素[18-21]。在這些影響因素中，衛星表面材質在正常飛行過程中是穩定的，大氣粒子與衛星表面的入射夾角平均值在每個圈次也是穩定的，因此可以推斷，CD的周期性變化是由大氣成分、溫度的變化引起的，而大氣成分、溫度受太陽活動影響與F10.7 有很強的相關性。

圖6 精密定軌解算彈道系數和空間環境指數Fig.6 Ballistic coefficients solved by orbit determination and space environment index

圖6 精密定軌解算的彈道系數展示了與F10.7的相關性，說明解算結果是準確可信的，但是在定軌過程中，彈道系數不可避免地吸收了大氣模型、軌道測量數據和飛行器迎風面積等因素引入的誤差，為此，需要對每日解算的彈道系數進行平滑處理。由于天宮一號彈道系數中的質量和迎風面積是精確已知的，從彈道系數中提取了物理意義更加明確的大氣阻力系數。如圖7 所示，分別采用3 d、5 d 和7 d 均值法對阻力系數平滑?？梢钥闯?，平滑處理后，不僅去除了隨機誤差，還保留了原來的周期性變化。因此，選取7 d 平滑處理后的阻力系數與TLE 解算的阻力系數比對。

圖7 精密定軌解算的大氣阻力系數及其平滑值（起始日期2012-01-01）Fig.7 Drag coefficients solved by orbit determination and their smoothed values (start date Jan.1, 2012)

利用TLE 解算天宮一號2012 年全年的大氣阻力系數CD和彈道系數，解算過程中采用了不同時間間隔的TLE 構建觀測量 ?atle，即式（6）中的積分區間，分別設置為2 d 和3 d；觀測弧段分別選取7 d和10 d，即分別采用7 d 和10 d 的TLE 解算一個彈道系數。4 組解算結果如圖8 所示，同時給出了精密定軌的解算結果?？梢钥闯?，TLE 解算的4 組結果與精密定軌解算結果的變化趨勢十分吻合，都能反映出太陽活動的周期性變化；不同時間間隔的觀測量對解算結果影響不大；觀測弧段為10 d 的解算結果更穩定，曲線更加平滑。

圖8 TLE 解算的大氣阻力系數（起始日期2012-01-01）Fig.8 Drag coefficients solved by TLE (start date Jan.1, 2012)

2.2 天宮一號隕落分析

2.2.1 隕落階段彈道系數解算分析

天宮一號是中國第一個非受控再入大氣層的超大型航天器，于2016 年3 月16 日正式終止數據服務。在大氣阻力下影響下，天宮一號軌道不斷衰減，最終于2018 年4 月2 日8 時15 分落入南太平洋。天宮一號最后30 d 平均軌道高度如圖9 所示。隕落之前數月難以得到精確的姿態信息，從而難以估算式（1）中的迎風面積，因此適宜采用1.1 節的方法，對包含了CD和迎風面積的彈道系數聯合估算。

圖9 天宮一號隕落前軌道高度變化Fig.9 Height variation of Tiangong-1 before reentry

使用1.2 節的方法對隕落前40 d 的TLE 的半長軸進行平滑檢測，結果如圖10 所示。隕落前40 d至隕落前23 d，平滑值與TLE 半長軸差異20 m 之內，說明這期間飛行器高度衰減變化穩定，TLE 精度較高。隕落前22 d 至隕落，平滑值與TLE 半長軸差異逐漸增大，隕落前2 d 的差值甚至超過200 m，初步判斷可能是以下原因造成的：①隕落前15 d 至隕落前14 d 有地磁擾動；②飛行器姿態不穩定導致迎風面積變化；③200 km 以下大氣密度更加復雜多變；④TLE 本身的精度較差。前3 個原因造成的差異不影響彈道系數解算，不需要剔除，只需剔除精度較差的TLE。經過二次判別，隕落前3 d 至隕落的數據大部分不可用，因此后續軌道預報中隕落前2 d 和隕落前1 d 的彈道系數使用隕落前3 d 的結果。

圖10 天宮一號TLE 半長軸平滑檢測結果Fig.10 Semi-major axis smoothing test results for Tiangong-1 TLE

基于2 d 間隔的TLE 構建觀測量 ?atle，分別采用每5 d、7 d 和10 d 的觀測弧段解算一個彈道系數，結果如圖11 所示。彈道系數隨著觀測弧段的增長，彈道系數的變化趨勢趨于平緩。文獻[22]利用Sentman 模型計算了150～300 km 高度范圍不同構型衛星的大氣阻力系數，認為圓柱體衛星（與天宮一號構型較接近）的阻力系數在2.3～2.7 之間變化，最大值是最小值的1.17 倍。圖11 給出的類似高 度 范 圍 的 彈 道 系 數（TLE-10d-2d）在4.7×10?3～5.4×10?3之間變化，最大值是最小值的1.15 倍，變化幅度與Sentman 模型計算結果接近。

圖11 隕落前30 d 不同觀測弧段解算的彈道系數Fig.11 Ballistic coefficients solved with different observation arcs before 30 days of reentry

2.2.2 天宮一號隕落時間預報

從2018 年3 月3 日起，使用2 種方法基于美國空間監視網每日發布的天宮一號TLE 對隕落時間進行預報。第1 種方法：利用SGP4 模型計算TLE歷元時刻的衛星位置速度 (r,r˙)，以此作為初始軌道，利用自主定軌軟件做軌道預報，計算衛星高度降至125 km 的時刻Tr，預報過程中使用64×64 階重力場模型和MSISE00 大氣密度模型，空間環境參數太陽流量F10.7 和地磁指數AP 使用隕落前30 d 的均值，彈道系數使用圖11 中TLE-5d-2d 解算值；第2 種方法：利用SGP4 模型和TLE 進行軌道預報，計算衛星高度降至125 km 的時刻Tr′。用上述2 種方法逐日預報天宮一號隕落時間，如圖12 所示?？梢钥闯?，利用第1 種方法解算的彈道系數預報天宮一號隕落，最大誤差15 d，平均誤差7.03 d，精度顯著優于第2 種方法（最大誤差25 d，平均誤差15.6 d）。

圖12 使用不同彈道系數的預報誤差Fig.12 Prediction errors using different ballistic coefficients

天宮一號隕落事件在國際上備受關注，世界主要航天大國都發布了各自預報的再入窗口，如圖13所示。統計了俄羅斯國家航天集團、美國國家航空航天局、歐洲航天局、日本宇宙航空研究開發機構和印度空間研究組織?4 d～?1 d 發布的再入窗口，其中缺少美國國家航空航天局（?3 d）、歐洲航天局（?2 d）、俄羅斯國家航天集團（?1 d）和日本宇宙航空研究開發機構（?1 d）的數據。需要說明的是，在軌道預報中使用的太陽輻射指數F10.7 和地磁指數AP 是隕落前30 d 的實測平均值，因此可直接計算出航天器再入大氣層（125 km）的時刻，而各國是基于F10.7 和AP 的一個變化范圍做出的軌道預報，因此發布的再入窗口也是一個時間范圍，為方便比較，取各國發布的時間窗口的中間值。從圖13 的統計結果可看出，本文第?4 d、?1 d 的預報精度最高，第?3 d、?2 d 的預報精度處于中間水平，總的來看，本文的隕落預報精度在各航天大國中處于領先水平，說明彈道系數解算方法是正確可行的。

圖13 各國與本文發布的天宮一號隕落預報誤差比較Fig.13 Comparison of Tiangong-1 reentry forecast errors issued by various countries and this paper

3 結 論

本文研究了利用TLE 的半長軸衰減反演非合作目標彈道系數的方法，總結如下：

1）對TLE 的半長軸精度進行了評估，驗證了半長軸衰減信息用于反演彈道系數的可行性，確定了使用2 d 間隔的TLE 半長軸構造軌道衰減觀測量。

2）研究了TLE 預處理方法，通過多項式平滑檢測和二次判別識別出野值、軌道機動和地磁暴3 種情況。

3）利用天宮一號2012 年TLE 解算了全年的彈道系數，并與精密定軌解算的彈道系數進行了比對，兩者變化趨勢十分吻合，且都能反映出太陽活動的周期性變化。

4）解算了天宮一號隕落前30 d 的每日彈道系數，并利用解算結果進行了隕落預報，與世界其他國家發布的預報比較，解算的彈道系數能夠獲得優于其他國家的預報精度。

5）利用TLE 的半長軸衰減信息解算非合作目標彈道系數在航天工程中是切實可行的方法，對于非合作目標的監視，碰撞預警和規避提供了更精準的目標特性信息，對提升中國空間安全和態勢感知能力具有重要意義。

本文提出的彈道系數解算算法僅適用軌道高度600 km 以下的目標，對于軌道高度400～600 km的航天器，其半長軸衰減較慢，在半長軸衰減觀測量時應擴大時間間隔。

在后續工作中，計劃對不同形狀、不同高度及不同太陽活動條件下的航天器彈道系數進行解算，分析該方法在不同條件下的適用性。