李佩鈺， 楊 震?， 羅亞中

（1.國防科技大學空天科學學院， 長沙 410073； 2.空天任務智能規(guī)劃與仿真湖南省重點實驗室， 長沙 410073）

1 引言

隨著航天技術的飛速發(fā)展，在軌空間目標的數量不斷增加，太空資產面臨的碰撞風險顯著提高。 截止2023 年6 月，人類成功送入軌道的航天器超過15 797 顆，其中約有8218 顆仍然在軌活動［1］。 且在地球軌道上存在約有50 萬顆1 ～10 cm 大小的碎片，1 mm 或更小的碎片估計超過1 億顆［2］。 此外，星鏈（Starlink）等大型商業(yè)星座的部署［3］，進一步加劇了軌道環(huán)境的擁擠程度，衛(wèi)星碰撞的可能性大幅增長。 由于空間目標在太空高速飛行，微小的碰撞也會產生嚴重的后果，可能導致航天器失效甚至解體［2］。 太空資產是國家戰(zhàn)略資產［4］，通過軌道機動規(guī)避碰撞風險是保護太空資產的重要手段。

碰撞規(guī)避的首要步驟一般是對目標進行接近分析，即通過定軌結果預報兩目標的接近距離、接近時刻等參數。 現有航天器接近分析方法主要有解析法和數值法2 類。 解析法是基于軌道根數之間幾何關系設計解析表達式，通過求封閉函數的解，獲取接近事件的信息。 Hoots 等［5］通過一系列篩選排除大量目標，而后考慮共面或異面情況，確定兩目標之間距離小于設定閾值時的最接近點和接近區(qū)間。 目前解析法相關研究大多是基于Hoots 方法開展［6-7］。 由于解析法對軌道類型十分敏感，在攝動影響下存在漏報的問題。 而隨著計算能力提升，數值法成為主要接近分析方法，即采用逼近、插值等近似計算方法，獲得目標接近時刻狀態(tài)信息。 基于部分已知軌道狀態(tài)，Alfano等［8］采用數值計算方法，通過計算定義距離函數的三次多項式方程得到接近時刻，并推導了誤差橢球區(qū)域進出點及進出時間，此方法即A-N（Alfano-Negron）算法。 李鑒等［9］針對A-N 算法中多項式易丟根情況，對其求根方法進行了改進。Denenberg 等［10］基于A-N 算法，對結果精度和速度進行折中計算，建立了更高效的接近分析方法。由于航天器可能存在與同一個目標在一定時間段內出現多次接近，或與多個空間目標同時存在接近的風險，因此為確保航天器的控后軌道跟任何已知目標在一定時間范圍內都不能發(fā)生接近風險，需要研究考慮多目標多圈次場景下的適用于碰撞規(guī)避的接近分析方法。

基于接近分析得到接近距離和接近時刻等參數后，需定義并計算碰撞風險評估指標來衡量碰撞規(guī)避過程中目標間的碰撞風險。 碰撞風險評估指標的計算方法目前主要有Box 區(qū)域判定和計算碰撞概率2 種。 Box 區(qū)域判定方法對位置關系進行考慮，是航天飛機碰撞預警的主要手段［11］。 李甲龍等［12］對碰撞風險的評估適用標準進行分析，在預報精度較差時，可利用Box 區(qū)域判定對風險進行綜合評定。 李翠蘭等［13］考慮軌道預報及軌控時的偏差傳播特性，提出了Box 動態(tài)閾值設定方法。 碰撞概率計算考慮了空間目標之間的相對幾何關系以及位置速度的誤差協方差，是目前使用較為廣泛的風險評估方法。 Chan［14］針對航天器長期碰撞建模問題，考慮橢球協方差誤差，利用3 個隨機變量替代笛卡爾坐標系下的位置速度，計算碰撞概率。 Patera［15－16］利用三維對稱協方差誤差矩陣，解決非線性相對運動下的碰撞概率問題，且為準確得出概率增量，將相遇段分成片段，將增量求和得到總概率值。 計算碰撞概率的方法在工程上易受測量精度及動力學模型的影響，當軌道精度不足時，難以保證誤差協方差的準確，使用起來并不方便。 對比之下采用Box 區(qū)域估計方法進行風險評估簡單效率較高，是更可靠更保守的方法［17－18］。

當航天器通過風險評估識別到碰撞風險時，需采取相應規(guī)避措施，即碰撞規(guī)避策略。 Patera等［19］基于開普勒傳播對空間目標狀態(tài)進行預報，并采用梯度法確定機動方向，進一步用一維求根方法得到機動大小，能夠將碰撞概率降低到預期。于大騰等［20］定義了潛在的碰撞威脅區(qū)，以弧長為優(yōu)化目標建立最優(yōu)規(guī)避模型。 基于航天器機動可達域，張賽等［21］定義了威脅評價指標，實現了最小燃料下的多脈沖威脅規(guī)避方法。 Bombardelli等［22－23］利用數值方法對施加的脈沖機動與B 平面的位移之間的相對關系進行優(yōu)化，得到使得碰撞脫靶量最大或碰撞概率最小時的機動規(guī)避策略。 以上理論研究針對碰撞規(guī)避機動優(yōu)化問題皆提出了相應的解決方法。 但為提高規(guī)避策略的高效且兼顧可靠性，需要從控后軌跡的建立方式和考慮控制工程約束的方向出發(fā)，對機動時間及機動量進行優(yōu)化，獲得燃料消耗相對較小的碰撞規(guī)避機動方案。 現有研究考慮工程任務需求和約束的較少，且大多關注單次接近情況，對多圈次接近考慮不多。 同時若采用數值的高精度軌道積分預報避撞后的控后軌道，將使得優(yōu)化環(huán)節(jié)計算量大，生成避撞機動方案的效率低。

本文基于航天器和空間目標軌道狀態(tài)數據，首先利用高精度非線性相對運動解析方程［24］，快速構建控后飛行軌跡；而后以步長推進，通過航天器與空間目標的長時間多次接近分析，有效識別分析時段的全部接近信息，同時進行碰撞風險評估；最后建立優(yōu)化模型，以半長軸為控制量，實現脈沖機動下的最優(yōu)碰撞規(guī)避控制量和變軌時間選擇，從而提高碰撞規(guī)避策略的安全性。

2 問題描述

碰撞規(guī)避背景下的機動規(guī)劃需要對目標之間的相互關系進行分析，識別與航天器碰撞風險較大的空間目標。 由于空間目標數量和分析時間不斷增加，如何快速預報風險且提高工程上應用的適配性，一直以來是軌控安全性分析關注的問題。本文基于航天器精密星歷及空間目標編目TLE（Two Line Elements）數據庫，基于解析非線性相對運動方程快速構建控后軌跡，建立碰撞規(guī)避機動規(guī)劃模型，流程圖如圖1 所示。

圖1 碰撞規(guī)避機動規(guī)劃方法流程Fig.1 Flowchart of collision avoidance maneuver planning method

首先，采用橢圓參考軌道，且包含J2 攝動項和中心引力梯度二階項的解析非線性相對運動方程，快速構建控后軌跡。

其次，針對航天器與空間目標庫大量數據，解決多圈次多目標下的接近分析問題。 將軌道信息粗篩選與多項式求根方法相結合，得到最接近距離對應的TCA（Time of Closest Approach）時刻、位置、速度等信息，同時進行基于等效距離量化指標的碰撞風險評估。

最后，由于航天器長期軌道維持一般只實施橫向控制，抬高或降低軌道半長軸，因此可用半長軸為控制量。 在允許的控制時間和控制量區(qū)間內，優(yōu)化規(guī)避機動時刻與機動量。

在不考慮控制偏差影響的情況下，再次調用接近分析模型，衡量控后軌跡的風險大小，獲得該機動量下的所有接近分析目標的軌道信息，排除再次引入風險的可能性。

3 碰撞風險分析基礎模型

3.1 坐標系定義與轉移矩陣

定義地心慣性坐標系（Earth Central Inertial，ECI）如圖2 所示。 其基準平面為天體平均赤道面，定義坐標原點O為地心，X軸位于平均赤道平面內并指向歷平春分點，Z軸沿著垂直于赤道平面并指向北極，Y軸符合右手笛卡爾坐標系。

圖2 航天器軌道坐標系Fig.2 Spacecraft orbital coordinate system

定義航天器當地軌道坐標系（Local Vertical，Local Horizontal， LVLH）如圖2 所示。 其坐標原點為航天器，基準平面為航天器或空間目標的軌道平面，x軸沿著航天器位置矢量方向，z軸垂直于軌道平面，為正法線方向，y軸根據右手定則確定。

本文使用非線性相對運動方程時，航天器LVLH 坐標系與ECI 坐標系相互轉換矩陣見式（1）［24］。

式中，r（t） 和v（t） 分別表示航天器位置矢量和速度矢量。

3.2 接近分析模型

由于空間目標庫數據量較大，進行接近分析時需首先對目標進行篩選，本文分別利用軌道歷元日期篩選、遠近地點篩選和最大速度及加速度篩選，將大部分不可能接近的空間目標排除。 針對篩選后的目標，結合多項式方程求解所有目標的最接近時刻對應的距離和角度等軌道信息，即得到碰撞規(guī)避規(guī)劃的風險量化參數。 接近分析的流程如圖3 所示。

圖3 接近分析過程Fig.3 Procedure of approach analysis

3.2.1 軌道接近事件篩選

已知碰撞分析的起始時刻t0和終止時刻tf，以及航天器精密星歷，將該星歷按時間步長Δt給出，時刻ti，i＝0，1，2，…，n－1 滿足ti＋1＝ti＋Δt，tn＝tf，使用七階拉格朗日插值法得到其位置和速度矢量Xch（ti） ＝[rch（ti），vch（ti）]T。 同時，根據空間目標編目數據庫中每個空間目標的兩行軌道根數TLE，采用SGP4 軌道預報模型，獲得任意時刻空間目標的位置和速度矢量。

首先，篩選軌道歷元日期。 提取空間目標TLE 數據中的歷元時刻tTLE信息，設空間目標編目數據庫中距離碰撞分析的起始時刻t0過長時（例如超過30 d），則不需要進行接近分析。 因此，軌道歷元日期的篩選條件為式（3）。

而后，考慮目標之間的相對幾何關系，分別計算分析開始時刻t0對應的航天器和空間目標的遠地點及近地點高度，使用近地點－遠地點篩選法，將不可能存在的碰撞點排除掉，篩選的標準［25］為式（4）。

式中，Q為遠地點高度中的航天器和空間目標中的小者，q為近地點高度中的大者，D1為設定的近地點－遠地點門限值。

對于通過篩選的空間目標，定義D2為危險距離，篩選掉大于該門限值的接近事件。 在碰撞分析的時間區(qū)間內，假設航天器在預報區(qū)間內以第二宇宙速度飛行時，依據步長分析航天器和空間目標之間的相對位置矢量rc＝［rc，x，rc，y，rc，z］T，設Rc，1為最大速度相對距離衡量指標［26］見式（5）。

式中，ve＝ 2μ／r。 若rc，x ＞Rc，1或rc，y ＞Rc，1或rc，z ＞Rc，1或rc＞Rc，1，那么在此時間區(qū)間內不可能發(fā)生接近事件。

最后，考慮地球引力加速度以及航天器和空間目標的相對速度，定義Rc，2為最大加速度相對距離衡量指標［26］見式（6）。

式中，D3為危險距離，若rc＞Rc，2，則通過安全距離篩選，該指標可進一步進行篩選。

3.2.2 接近時間識別

在ECI 坐標系中，根據航天器和空間目標的相對距離矢量為rc， 得到相對距離函數fc（t） 及其一二階導數見式（7）［8］。

式中，rh和rs分別為航天器和空間目標位置矢量，rc＝rh－rs。

而后基于距離變化函數Pc（τ） 的三次多項式方程式（8），得到相對距離函數fc（t） 達到極小值點時的根τTCA，從而得到其對應的時間tTCA見式（8）［8］。

3.3 碰撞風險評估指標

根據接近分析結果以及軌道預報誤差在橫向較大，徑向較小的理論基礎。 基于工程經驗建立的碰撞風險評估指標模型［17］。 將加權后的接近距離和徑向距離中的最大值等效為航天器與空間目標之間的最近距離，作為碰撞風險評估的等效距離指標見式（9）。

式中，l＝1，2，…，N，j＝1，2，…，pl，N為在時間區(qū)間［tm，tf］內與航天器接近距離小于設定值dtol的空間目標數目，pl為第l個空間目標與航天器在時間區(qū)間［tm，tf］內的總接近次數，djl（tTCA）為航天器與第l個空間目標的第j個最接近距離，δxjl（tTCA） 為航天器與第l個空間目標在第j個最接近時刻相對位置矢量的徑向分量。

等效距離指標sjl同時考慮最接近距離和徑向距離，能安全高效進行航天器與空間目標間的風險評估。 當sjl越小時，碰撞風險越大，sjl越大時，碰撞風險越小。

同時，定義碰撞預警的黃色門限為Syellow、紅色門限為Sred， 若sjl≤Syellow， 則為黃色警報，若則為紅色警報。 對于紅色警報的處置，一般需要通過執(zhí)行軌道機動來規(guī)避碰撞風險，同時復核該避撞機動，保證不引起與所有編目空間目標新的碰撞風險。

4 避撞機動優(yōu)化模型

4.1 控后軌跡快速生成

在碰撞規(guī)避的軌控策略規(guī)劃中，由于空間目標數據庫所包含的軌道信息量大，因此對于多圈次的碰撞安全分析，需要更高效的控后軌跡預報方式。 本文根據非線性相對運動方程，快速建立基于標稱軌道的控后軌道。

自初始時刻，將標稱狀態(tài)Xi加上該時刻相對狀態(tài)偏移量δxi，按照一定的步長高精度預報軌道，得到任意時刻的航天器相對運動狀態(tài)Xi， 控后軌道構建示意圖如圖4 所示。

圖4 軌控示意圖Fig.4 Diagram of orbit control

已知需要執(zhí)行碰撞規(guī)避的變量值為y＝［tm，Δam］，其中， Δam為半長軸控制量，根據高斯攝動方程，對應在航天器當地軌道坐標系下的橫向脈沖機動控制量Δvt為式（10）。

定義tm時刻后任意t（t ＞tm）時刻航天器的控后軌道，相對于控前軌道的相對運動狀態(tài)為式（11）。

以時間t為自變量，利用解析非線性相對運動方程，將δx(tm) 預報到tm時刻后的任意t時刻（t ＞tm），獲得預報控后軌道的相對運動狀態(tài)δx(t) ，見式（12）［24］。

式中，算子?表示張量積，即將2 個不同矢量空間中的矢量相乘得到一個更高維空間中的矢量。 狀態(tài)轉移矩陣Φ(t，tm) 與狀態(tài)轉移張量Ψ(t，tm) ，由航天器預報時間△t＝t－tm和基于精密星歷的非奇異修正初始軌道的軌道根數em＝［a，w＋f，i，ecosw，esinw］ 計算，e為偏心率。P～－1（tm，tm）＝D－1（tm）P－1（tm），D（t） 為線性變換矩陣，P(t)為6× 6 的雅克比矩陣，Q(t) 為6×6× 6 的Hessian 矩陣。

最后，基于以上非線性映射，以及航天器控前軌道精密星歷信息，得到控后軌道的任意時刻相對狀態(tài)X～ch（ti），如式（13） 所示。

4.2 避撞機動優(yōu)化建模

碰撞規(guī)避機動優(yōu)化變量一般包括機動大小、機動時間、機動方向等。 本文基于航天器控前精密星歷數據，通過七階拉格朗日插值法得到的航天器位置矢量rch（t） 和速度矢量vch（t）。

對于控制約束，考慮測控跟蹤計劃，根據最早可控制時刻tlow和最遲控制時刻tup， 確定規(guī)避機動控制時間窗口。 同時，基于航天器平臺機動能力限制和軌道維持約束，確定控制量門限，設定半長軸控制量上下限分別為Δaup和Δalow。

設碰撞規(guī)避機動執(zhí)行時刻為tm，半長軸控制量為Δam，碰撞規(guī)避機動后的控后飛行軌跡與空間目標間的最近距離門限值Stol。 碰撞規(guī)避優(yōu)化的設計變量為式（14）。

采用序列二次規(guī)劃算法求解優(yōu)化問題，優(yōu)化模型應用實施具體步驟如下：

1）根據工程需要和規(guī)避需求，基于第3.2 節(jié)接近分析得到的風險量參數，運用式（9）的指標結果，進行碰撞風險評估。

2）調用第4.2 節(jié)的序列二次規(guī)劃算法優(yōu)化軌控安全策略，優(yōu)化迭代總速度，并應用式（12）及式（13）的非線性相對運動方程，解算控后軌道。 當時間區(qū)間［tm，tf］內的所有接近等效距離均大于給定的門限值Stol， 得到最優(yōu)控制時刻tc和半長軸控制量Δac。

3）根據基于優(yōu)化算法模型得到最優(yōu)時刻tc和半長軸控制量Δac， 再次調用第4 節(jié)的控后軌道預報模型，當航天器的最小接近距離大于設定的碰撞預警門限值時滿足該條件，即通過軌控安全性復核，是可行的規(guī)避控制策略。

5 算例分析

對前文所述方法和模型進行仿真校驗，仿真分析時段為2022／08／02 23：00：00 至2022／08／08 12：00：00，在軌航天器（Satellite）與控前存在碰撞風險的目標航天器（Target）的飛行軌跡通過工業(yè)部門獲得，空間目標庫TLE 數據由國外公開網站獲得，如表1 所示。

表1 航天器與空間目標仿真數據源Table 1 Simulation data source of spacecraft and space targets

5.1 自然接近分析

項目組研發(fā)的自主航天任務設計軟件ATK（Aerospace Tool Kit）編寫實現了論文方法，采用ATK 仿真分析航天器與空間目標間的接近情況，設定接近距離門限100 km，分別取Syellow＝200 m，Sred＝50 m，得到接近距離較小的前3 個接近事件如表2 所示。

表2 航天器與空間目標接近事件列表Table 2 List of spacecraft and space target approach events

由表2 可知，Target 與空間目標國際編號53 316 的等效距離都小于200 m，都存在碰撞風險。 其中，接近時刻t?為8 月5 日17：25：54 的等效距離最近，此距離下的航天器（Satellite）與潛在接近航天器（Target）的接近過程，如圖5 所示。

圖5 ATK 仿真控前接近時刻Fig.5 ATK simulation of pre－control approach time

5.2 碰撞規(guī)避結果

根據航天器跟蹤測控計劃，確定規(guī)避機動控制時間窗口。 同時，由于Target 與空間目標國際編號53 316 同時存在碰撞風險，為規(guī)避所有目標的碰撞風險，時間區(qū)間的選擇要早于碰撞風險時刻2022／08／03 20：28：34。

由此設置最早可控制時刻tm為8 月3 日15：37：14，最遲控制時刻tf為8 月3 日19：15：24，采用橫向控制方式，抬高或降低軌道半長軸，控制量下限為Δalow＝－ 100 m，上限為Δaup＝ 100 m。調用序列二次優(yōu)化算法，規(guī)劃最優(yōu)規(guī)避機動策略。

選取各機動策略tci和Δai，在原來自然交會接近時刻t?的前后多個軌道周期內，分析規(guī)避機動控后軌道的多次接近最小等效距離，i＝1，2，3，…，得到如表3 所示的航天器與空間目標規(guī)避策略。 通過優(yōu)化模型得到最優(yōu)控制時刻和半長軸控制量，即第一個規(guī)避策略。 同時，依據等效距離排序，選取其他2 個滿足s?i≥200 m 且控制量較小的可行規(guī)避機動策略。

表3 航天器與空間目標規(guī)避策略Table 3 Avoidance strategy of spacecraft and space targets

5.3 機動復核結果

而后擬實施機動策略規(guī)避效果復核，分別選取3 個規(guī)避策略。 復核計算控制后在原來自然交會接近時刻t?的前后多次接近條件，如表4所示。

表4 3 個規(guī)避策略的控后規(guī)避效果Table 4 Control effects of three avoidance strategies

經過復核比對3 個可行的軌控策略的規(guī)避效果，得到最優(yōu)控制時間和控制量給出策略，可在最小機動量下，規(guī)避所有目標。 最優(yōu)控制時刻tc為2022 年8 月3 日17：25：09，最優(yōu)半長軸控制量為Δa＝－37.114 m。 利用ATK 軟件仿真航天器與空間目標的最優(yōu)機動下的規(guī)避結果，如圖6 所示。

圖6 ATK 仿真控后接近時刻Fig.6 ATK simulation of approach time after control

綜上所述，本文提出的碰撞規(guī)避機動規(guī)劃方法，對連續(xù)多個軌道周期出現多次接近的碰撞規(guī)避問題，能夠避免控制后短期內再次出現碰撞風險的問題。

6 結論

本文針對空間目標碰撞風險持續(xù)上升，且境內可軌控測控時間區(qū)間限制嚴格的碰撞規(guī)避機動優(yōu)化問題，提出了基于解析非線性方程的控后軌跡快速構建與碰撞規(guī)避機動優(yōu)化方法。

1） 基于解析非線性相對運動方程快速構建控后軌跡，結合接近分析方法，提高了多目標碰撞風險分析效率；

2） 構建避撞機動優(yōu)化模型，得到最優(yōu)脈沖機動下的控制變量和時間，并建立了基于等效距離指標評估和復核碰撞規(guī)避效果的模型。

3） 仿真結果表明，所提方法能夠準確識別在軌航天器的碰撞風險，可獲得最優(yōu)控制時間及半長軸控制量。