面對接近式威脅的在軌飛行器安全性分析方法

于大騰

(北京跟蹤與通信技術研究所，北京 100094)

0 引 言

空間安全在國家安全中具有極其重要的戰略地位[1-2]，隨著在軌飛行器在未來起到的作用越來越重要，其可能受到的威脅也與日俱增，分析其在軌安全性并據此及時進行下一步規避決策，增強自身生存力顯得十分重要。目前，在軌飛行器面臨的威脅主要為非合作接近式威脅，不僅包括傳統的空間碎片等，更有部分潛在的主動接近式威脅。對于這類接近式威脅，通過規避機動進行防護是可行的。而進行規避防護的第一步便是安全性分析，所得結果可對后續自主規避決策提供支撐。本文研究飛行器面對接近式威脅時的安全性分析方法展開討論。

在安全性分析方面，諸多學者對空間碎片等傳統威脅進行了非常廣泛的研究[3-4]。對于空間目標接近分析問題，1984年Hoots等[5]提出了一種通過相位、高度和軌道幾何的解析方法，以此來確定接近目標的最小相對距離。文獻[6-10]利用軌道根數，對共焦點的開普勒軌道最小相對距離以及距離函數的臨界點解析計算方法進行了研究。Alfano[11]和Negron等[12]引入橢球函數和相對距離函數，提出了一種空間目標接近數值分析算法(即A-N算法)。Klinkrad等[13]和Alarcon-Rodriguez等[14]提出了一種多步篩選運動學分析算法，首先進行高度篩選，隨后進行粗略距離篩選和細致距離篩選，最后對距離一階導數求根來確定最小相對距離對應的時刻。Luo等[15]根據橢球距離、碰撞概率、預警門限等安全評價要素對GEO軌道軌跡安全性進行了分析與規劃。以上均是對接近問題進行分析，利用兩目標的軌道數據，設定接近門限值(給定的距離門限或者安全橢球)和分析起止時間，計算兩目標間相對距離到達距離門限或進入安全橢球范圍的時間。

在空間碎片碰撞風險評估方面，目前多采用碰撞概率作為評價指標，通過地面對空間目標的觀測和軌道預報計算得到評價指標，并據此評估碰撞風險。Chan[16]給出了碰撞概率的計算方法，證明如果飛行器和空間碎片誤差的協方差矩陣不相關，則它們的和即是相對位置誤差的協方差矩陣，故碰撞概率計算轉化為相遇平面二維概率問題。白顯宗等[17]推導了碰撞概率顯式表達式，對空間碎片的碰撞預警、碰撞危險等問題進行了研究。Patera[18]將二維積分問題轉化為一維曲線積分問題，簡化了不規則形狀飛行器之間的碰撞概率計算方法，并且提高了計算速度。Alfano[19]推導得到了一種級數表達式，該表達式由誤差函數和指數函數進行表示，同時還得到了最小級數項數。周威萍等[20]對比分析了圓軌道和橢圓軌道的碰撞通量。殷建豐等[21]基于相對軌道要素的方法，建立了計算碰撞概率的數學模型。王曉偉[22]研究了碰撞概率算法Cube模型參數對空間碎片演化模型的影響。以上方法中，Chan和白顯宗的方法是解析法，而Patera和Alfano的方法是數值法。

目前，航天領域的安全性分析多以空間碎片和失效飛行器為威脅目標進行研究，安全性指標主要為碰撞概率和相對距離。未來非合作主動接近飛行器等在軌目標的威脅日趨嚴重，現有方法在面對這些具有自主機動能力的飛行器時將難以進行有效判斷。如何對這些潛在威脅的機動軌跡進行分析，判斷其對在軌目標的威脅程度是十分必要的。本文將就這些問題開展研究，嘗試給出一種安全性分析解決方案。為便于統一，在后續研究中將我方飛行器稱為目標飛行器(簡稱目標器，Target)，具有自主交會能力的非合作接近飛行器稱為追蹤飛行器(簡稱追蹤器，Pursuer)。

1 安全性分析流程

追蹤器通常具有較強的機動變軌能力，一般的單次預測分析往往難以達到目的。本節將對追蹤器常見的交會過程進行分析，從而為安全性模型設計相應的算法分析流程。

本節安全性分析的基本思想為化繁為簡，在盡可能多的考慮重要安全要素的情況下兼顧計算效率。按照該思想對安全性分析流程進行設計，其基本架構與流程如圖1所示。

圖1 安全性分析流程圖

首先通過情報等渠道給出追蹤器當前脈沖機動能力，根據其剩余脈沖沖量計算軌道高度與軌道傾角的變化范圍，即威脅范圍。將我方目標器軌道參數與威脅范圍進行比較，處于威脅范圍之外的即可認為威脅度為零并終止后續計算，僅對處于威脅范圍內的航天器進行后續安全性分析。這里需要說明的是，為了提高計算效率，本文僅取軌道高度與傾角兩項作為威脅范圍參考指標。

當我方目標器處于追蹤器威脅范圍時，進入安全性分析階段。為了安全性分析的通用性以及合理性，本文選擇脈沖能力威脅、交會時間以及相對距離三個參數作為安全性特征指標進行分析。從追蹤器與目標器的相對軌道位置來分，可分為共面軌道和異面軌道兩種，分別建立相應的威脅評估模型，對更新的軌道數據進行計算，得到各指標的評估威脅度。

最后，在以上威脅度計算結果的基礎上，引入多指標加權綜合決策模型，將各個安全性特征指標下的威脅度進行歸一化，并利用專家權重系數將各項指標下的威脅度相加得到最終的多指標加權綜合威脅度，即可利用設置的相應門限值進行最終決策。

下文對已知剩余脈沖沖量情況下的威脅范圍快速估算方法進行介紹。將脈沖速度增量視為加速度與時間間隔的乘積，可由軌道攝動相關原理推導得到由脈沖沖量產生的軌道根數瞬時變化為[23]

(1)

式中:Δa和Δi為追蹤器與目標器的軌道根數差值，n為軌道角速度，eP,fP和EP分別為偏心率、真近點角以及偏近點角。ΔvPr和ΔvPt分別為所需的徑向和跡向脈沖沖量。r為當前位置地心距，u為軌道坐標系相對于赤道直角坐標系繞z軸的旋轉角，ΔvPh為所需法向的脈沖沖量。

(2)

則按照式(2)計算結果即可對目標器進行快速篩選，在僅考慮半長軸與傾角改變量時，對處于追蹤器威脅范圍內的目標器作后續的安全評估分析。

2 安全性分析指標選擇與計算模型

對于追蹤器，其變軌通常包括升降軌、軌道圓化、軌道傾角調整、以及升交點赤經調整等，其特點為變軌種類多、涉及參數廣且不可預測。若按照相應類別分別進行安全性分析，不僅對實際導航的處理要求過高，而且也很難提供一種可以統一處理的分析算法。

本節從能量約束以及交會安全性的角度出發，在涉及的眾多計算參數中選取交會脈沖沖量、交會時間以及最小相對距離三種參數作為安全性分析指標。下文對以上三種安全性分析指標構建相應的安全評估計算模型。

2.1 脈沖沖量計算模型

對于追蹤器來說，實施變軌機動對目標器進行接近，并最終完成接近是其固定的攻擊模式。因此，交會所需脈沖的大小是影響目標器安全性的重要因素。這里以共面和異面兩種不同情況來對交會所需脈沖進行快速估計。

當追蹤器與目標器共面時，若兩飛行器均為圓軌道，采用霍曼變軌進行計算。設追蹤器軌道半徑為rP，速度為vP，目標器軌道半徑為rE，速度為vE，橢圓轉移軌道近拱點和遠拱點的速度分別為vTP和vTE，則有[23]

(3)

故兩次轉移所需脈沖沖量為

(4)

式(4)為rPrE時，式(4)右側取負即可。總脈沖沖量為兩者之和，即

Δv=Δv1+Δv2

(5)

當追蹤器與目標器為共面非圓軌道時，兩者的軌道形狀有很多種，若根據每種特例均進行分析計算將十分復雜且計算時間成本過高。為了簡化分析，提高計算速度，本節根據追蹤器的交會特點，即其最終目的是要與目標器交會，故兩者軌道的半長軸及偏心率會逐漸接近并最終相等。因此，這里取軌道根數中半長軸和偏心率的變化作為所需脈沖計算依據。易知，沖量與半長軸和偏心率的變化關系為[23-24]

(6)

式中:Δa和Δe為追蹤器與目標器的軌道根數差值，n為軌道角速度，eP,fP和EP分別為偏心率、真近點角以及偏近點角。ΔvPr和ΔvPt分別為所需的徑向和跡向脈沖沖量。

式(6)中，Δa和Δe為已知量，聯立即可求得ΔvPr和ΔvPt，總脈沖沖量Δv為兩者之和。需要指出的是，為尋找最優脈沖，應對fP進行簡單遍歷，從而提高計算效率。

當追蹤器與目標器異面時，本節選用雙脈沖Lambert算法對變軌接近所需脈沖進行簡單估計。如圖2所示，rP和vP0為追蹤器某一時刻的初始位置和初始速度，vPf為轉移軌道的初始速度，rE和vEf為目標器在期望交會時刻的位置和速度，vE0為轉移軌道的終端速度。雙脈沖Lambert問題可以描述為給定追蹤器某一時刻的位置rP和速度vPf，目標器在期望交會時刻的位置rE和速度vEf，以及追蹤器沿轉移軌道從rP到rE的飛行時間ΔtTr，進而確定追蹤器的兩次脈沖沖量Δv1和Δv2。

圖2 異面交會示意

顯然，求解Lambert問題實質上就是求解高斯問題。而高斯問題的求解，研究者們早已推導給出一組經典方程，該方程對rP，rE，rPf和vE0之間的關系進行了詮釋[25]

(7)

為求解軌道轉移的最優脈沖，本節將構建雙重遍歷模型。首先，外層遍歷模型為對目標器的軌道參數進行遍歷，以一定步長計算其一個軌道周期內各時刻的位置速度。其次，內部遍歷模型為對轉移軌道的轉移時長進行遍歷，轉移時長遍歷范圍為零到目標器一個軌道周期。軌道轉移時長可代入Lambert問題進行求解得到所需脈沖總沖量，在遍歷的數據中取最小值即是Lambert交會最優脈沖沖量解。計算流程如圖3所示。需要指出的是，由于上節中威脅范圍的篩選為了計算的便捷性設計得較為簡單，因此這里計算得出的最優脈沖沖量解集中的解仍有可能超出追蹤器的脈沖沖量上限。這種情況將在下一節的多指標加權綜合評價中加以考慮，以使綜合評價更合理。

圖3 Lambert最優轉移脈沖計算流程

至此，便將共面和異面情況下的脈沖沖量計算模型構建完成，實際應用中可根據不同的情況選擇相應模型進行脈沖沖量計算。所得脈沖沖量可代入下節中的多指標加權綜合評價模型中，為最終的安全性評估提供支撐。

2.2 交會時間計算模型

對于追蹤器來說，執行接近任務時，雙方很可能已經處于交戰狀態，時間的重要性大大提升，因此交會時間是除了脈沖沖量外另一個重要的參考要素。交會時間計算模型與脈沖計算模型是相匹配的，兩者的計算結果將在下節中進行綜合評價。

需要指出的是，本節所作工作是在追蹤器與目標器相位角理想情況下進行的，即最有利于追蹤器的情況。實際情況中若存在相位差，且追蹤器不進行相位等待直接機動，則需要結合當時相位差進行相應計算，所需要脈沖消耗將大于本章中情況。

當追蹤器與目標器為共面軌道時，假設采用霍曼轉移，最短交會時間即轉移時長為轉移軌道的半個周期，則有

(8)

當追蹤器與目標器為共面非圓軌道時，為了簡化計算模型種類，這里將其歸到異面軌道情況下，統一采用Lambert算法進行計算。在采用Lambert算法進行計算時，遍歷的轉移軌道時長即是交會時間，其與所得到的脈沖沖量解集是一一對應的。在第3節，將利用多指標加權綜合評價模型對脈沖沖量及其對應的交會時間進行安全性評估，并選擇綜合性威脅較大者作為追蹤器的交會策略進行安全性評估。

2.3 最小相對距離計算模型

相對距離在傳統空間碎片安全性分析中常作為一項核心指標，由其來對目標器的安全性進行評估具有簡單客觀和表征性強等特點。無論追蹤器進行何種復雜機動，其最終目的都是接近目標器，因此兩者相對距離在其路徑規劃的終端必然呈逐漸減少的趨勢，其對目標器的威脅也將越來越大。

本文利用這一特點將相對距離引入安全性分析模型，以當前輸入的追蹤器與目標器兩者的軌道狀態為時間零點，以目標器軌道周期為時長分別對追蹤器和目標器進行軌道外推，計算兩者的相對距離，并取一個周期內的最小值作為最小相對距離計算模型的輸出。

3 多指標加權綜合評價與決策模型

在第2節中對安全性指標進行選擇與計算后，后續安全性綜合評價需要解決各項指標權重的確定以及如何綜合的問題。對多指標進行綜合的方法目前有很多種，但應用最廣、意義最直觀的是加權和法，該方法將各個指標評價值與相應的指標權重相乘后累加求和，即可得到整個方案綜合評價值，從而根據評價值的大小對方案進行排序或者決策。當然也可以采用其它非線性形式，但經驗判斷、理論推導和仿真計算都表明，雖然加權和法較非線性形式更為簡單，但兩者計算的結果非常接近，而前者具有簡單以及容易理解的優勢[27]。

設空間飛行器安全性分析的各評價指標最終權重向量為

σ=(σ1,σ2,…,σm)σj∈[0,1],

(9)

則采用加權和法的安全性分析綜合評價與決策模型可以表示為

(10)

式中:Zi即當前軌跡xi的綜合評價值，根據Zi的大小即可給出追蹤器對目標器i的威脅評價。qij為軌跡i的指標j標準化值。

本節接下來重點對指標結果標準化以及指標權重的確定進行介紹。

3.1 指標標準化

在應用加權和法之前，首先需要將第2節指標計算結果進行標準化。本文采用常見的向量規范化方法對所選擇的三個指標進行標準化。

對于脈沖沖量指標，其當前位置對目標器進行交會所需脈沖越小，交會越容易實現，則相應的威脅就越大。考慮極限情況，當Δv→0時，即兩者已趨近于交會狀態；當Δv>VPmax時，即目標器超出追蹤器可達域，暫無威脅，記指標值為0。因此以追蹤器當前剩余脈沖估計值為參考量進行標準化，則第i條軌跡的脈沖沖量指標為

(11)

由于前述威脅范圍計算較為簡單，故此處計算仍存在超出當前剩余燃料的情況，當Δv>VPmax時，計qi1=0。

對于交會時間指標，易知其當前位置對目標器進行交會所需時間越短則相應的威脅就越大，當交會時間為0時，即已經交會，則指標值為1；當交會所需時間超過當前軌道軌道周期時間時，認為暫時沒有威脅，指標值為0。因此以追蹤器當前軌道周期為參考量進行標準化，第i條軌跡的交會時間指標為

(12)

對于相對距離指標，考慮未來一個軌道周期內追蹤器與目標器外推值的相對距離越大則相應的威脅就越小，在相對距離過大時繼續進行威脅考慮將沒有意義，因此引入一個相對距離參考值進行標準化，小于該參考值時進行威脅計算，相對距離為0時指標值最大為1，相對距離大于該參考值時，則認為暫無威脅，指標值為0，則第i條軌跡的相對距離指標為

(13)

通過以上計算便完成了指標結果的標準化，隨后需要對各個指標的權重進行計算。

3.2 權重計算模型

權重體現了安全性分析評價指標的重要程度，一般來說最終權重是主觀權重與客觀權重的綜合，比較常用的方法是采用加權求和的方式。設σs與σo分別為指標主觀權重和客觀權重的向量，則組合權重可表示如下

σ=γσs+(1-γ)σo

(14)

式中：σ為組合權重向量，γ和1-γ分別為主觀權重和客觀權重的權系數，γ可根據實際需要確定。

主觀權重常由專家打分得到，下文主要介紹客觀權重的計算方法。通常情況下，某個指標計算結果的標準差越大，表明該指標值的差異程度越大，則提供的信息量就越多，在綜合評價中其起的作用越大，因此其權重也應該越大。相反，某個指標計算結果的標準差越小，表明該指標值的差異程度越小，則提供的信息量也越小，在綜合評價中其所起的作用越小，其權重也應越小。

本節采用指標標準差來計算相對應的客觀權重，設第i條軌跡的指標qij的標準差為δj，可由多組不同情況下的追蹤器軌跡實際計算得到。則其對應的客觀權重計算公式為

(15)

至此，再結合式(14)即可確定組合權重。需要指出的是，由于脈沖沖量與交會時間具有一定的相關性，因此需要先利用前幾節的相關內容對標準化指標進行計算，計算脈沖沖量與交會時間各自的權重并進行指標融合，第i條軌跡的計算公式如下

Ji1=σ′1qi1+σ′2qi2

(16)

式中:σ′1為脈沖沖量的組合權重，且σ′1=γσ′s1+(1-γ)σ′o1，σ′2為交會時間的組合權重，且σ′2=γσ′s2+(1-γ)σ′o2。

隨后將融合后的指標Ji1作為一個整體再與最小相對距離指標qi3一起按照前述步驟計算各自的組合權重σ1和σ2，其計算方式與前述權重計算方式相同，并由此進行融合，由式(10)可得到最終的安全性綜合指標值Zi，具體公式如下

Zi=σ1Ji1+σ2qi3

(17)

式中：σ1為脈沖沖量與交會時間指標融合后的組合權重，且σ1=γσs1+(1-γ)σo1，σ2為相對距離的組合權重，且σ2=γσs2+(1-γ)σo2。

根據實際情況，可對該綜合性指標設置一定的閾值，超過該閾值即認定為威脅不可避免，準備采取相應的規避機動方法進行規避。

4 仿真校驗

本節用對前面建立的模型進行仿真校驗，為保證仿真與實際情況相接近，仿真取美國2015年發射的通信衛星AprizeSat-7/8以及電視衛星AlSat-2A作為假想目標器，相關軌道參數由國外公開網站查詢獲得。同時假設追蹤器初始位于與AprizeSat-7共面的停泊軌道上，其余參數隨機設置以保證仿真的普適性。以上飛行器各自的初始軌道根數如表1所示，軌道歷元為UTC時間2019年4月17日04時0分0秒。

表1 各飛行器初始軌道根數

以軌道歷元為初始時刻，設追蹤器采用三脈沖變軌策略對AprizeSat-8進行接近。為檢驗安全性分析模型是否具有普適性，這里取的三脈沖變軌策略為非最優策略。三次脈沖在J2000坐標系下的矢量分別為

vm1=[248.584,634.62,364.293]Tm/s

vm2=[1084.415,1.221,628.763]Tm/s

vm3=[-541.598,302.413,-348.289]Tm/s

以初始時刻為零點，第二和第三次脈沖施加時間為tm1=1030 s和tm2=2660 s。設數據觀測頻率為每5 min一次，總時長80 min，共16個觀測數據，則該變軌軌跡與上述三個目標器軌道參數可按照上節中建立的模型計算相應的指標參數。設相對距離參考值Rref=800000 m，追蹤器剩余脈沖估計矩陣為

VPmax=[3000, 2200, 2200, 2200, 2200, 1800,1800, 1800, 1800, 1800, 1800, 1800,1800, 1800, 1100, 1100]m/s

現假設專家打分得到的脈沖沖量與交會時間的主觀權重系數為σs1=[0.7, 0.3]，主觀權重的權系數為γ=0.75，客觀權重系數可利用軌跡數據結合,上文公式計算得到，經計算可得σo1=[0.549, 0.451]，則脈沖沖量的組合權重σ′1=0.663，交會時間的組合權重σ′2=0.337。根據各自的組合權重可計算得到脈沖沖量與交會時間的融合指標，類似的，該融合指標與最小相對距離的主觀權重系數為σs2=[0.45, 0.55]，主觀權重的權系數不變，則上述融合指標與最小相對距離的組合權重分別為σ1=0.441和σ2=0.559。為便于比較觀察，在進行指標綜合評價時將各指標計算結果均乘以100。

根據以上數據可對追蹤器安全性指標進行計算，脈沖和交會時間的融合指標如圖4所示，最小相對距離安全指標如所圖5示。

圖4 脈沖沖量與交會時間的安全性融合指標

圖5 最小相對距離安全性指標

由圖4可知，初始時刻追蹤器對AlSat-2A的脈沖沖量與交會時間融合指標最大，在第一次脈沖施加后，其對AlSat-2A的脈沖沖量與交會時間融合指標快速減小，其對AprizeSat-8的指標值逐漸增大，可見安全性模型可以判斷出了脈沖意圖。由圖5可知，第一次脈沖使最小相對距離指標快速增大后又開始逐漸衰減，在第二次脈沖施加后該指標開始穩步增加，直至第三次脈沖后指標達到峰值。

圖6為多指標加權后的綜合評價結果，由圖可見，AprizeSat-7全程不受任何威脅，AlSat-2A初始時刻受到的威脅與AprizeSat-8相近，而在追蹤器施加脈沖后，AlSat-2A所受威脅逐漸減小，AprizeSat-8所受威脅呈不規則增大趨勢。之所以增大趨勢不規則，是由于本文為考慮普適性而選擇的交會策略并非最優，但本文算法仍對每次脈沖意圖進行了準確判斷。可見該模型對交會軌跡具有一定的魯棒性。

圖6 多指標加權安全性綜合評價結果

為進一步驗證本文模型的合理性，下文假設追蹤器在上述過程中不施加任何機動，即處于待命狀態，對上述指標進行重新計算，所得結果如圖7、圖8和圖9所示。

圖7 脈沖沖量與交會時間的安全性融合指標(待命狀態)

圖8 最小相對距離安全性指標(待命狀態)

圖9 多指標加權安全性綜合評價結果(待命狀態)

從待命狀態下的計算結果可知，多數指標均維持在一個較低的范圍內上下波動，且綜合評價最大值也不超過30。待命狀態下雖然評價值均很低，但通過比較仍可以發現追蹤器潛伏狀態下對我方哪些飛行器具有威脅，并可根據綜合評價值的大小對潛在的機動位置進行有效估計。

5 結 論

本文提出的分析方法可推動在軌飛行器面對各種主動接近威脅時的安全性分析研究，對未來衛星防護體系及相關技術的研究具有一定的參考意義。主要結論如下：

1)對于在軌飛行器的安全性分析問題，本文以主動接近式威脅為例，分析給出了一種安全性分析綜合評價模型，拓展了目前以空間碎片等非自主接近威脅的研究目標集，使飛行器進行安全分析時具有更多選擇。

2)本文所建立的模型可對在軌飛行器的安全性進行有效分析。若對綜合評價值設置安全閾值，受威脅的飛行器可據此對規避機動窗口進行自主決策，提升復雜態勢下的自身生存能力。

3)在本文分析模型下，追蹤器施加多脈沖機動變軌對目標器進行接近，其每次脈沖均會在安全性評價值中得到體現，算法可對脈沖真實意圖進行準確判斷。甚至在追蹤器未機動狀態下，仍可發現其對我方哪些飛行器具有威脅，并對潛在機動位置進行有效估計。