基于拉普拉斯變換的空間目標碰撞概率計算方法

霍俞蓉， 李智， 韓蕾

(1. 裝備學院 航天裝備系， 北京 101416； 2. 裝備學院 航天指揮系， 北京 101416)

由于人類的航天活動日益增多，空間中在軌的空間目標數量不斷增大，因此航天器與空間碎片的碰撞風險越來越大。為了保證航天任務順利進行,保證航天員、航天器與空間環境的安全，空間碎片碰撞預警已成為航天任務設計與執行過程中必須進行的重要工作。

碰撞概率是目前國內外廣泛使用的判斷2個空間目標發生碰撞的可能性大小的重要判據。當前，國內外對碰撞概率方法的理論研究已經相對成熟。文獻[1]采用Monte Carlo方法得到了碰撞概率的計算結果；Foster和Estes[2]、Patera[3]和Alfano[4]提出了不同的碰撞概率數值積分模型；Chan[5-6]提出了一種將積分式進行近似的解析方法。

Monte Carlo方法以及Foster和Estes[2]提出的數值方法可以得到精度較高的碰撞概率結果但是計算速度較慢。Chan[5-6]解析方法計算速度快，但精度劣于數值方法。解析方法中比較經典的為Chan方法，根據文獻[7]中對該方法的精度分析，在概率密度函數(Probability Density Function，PDF)接近圓分布時，能夠較好地近似碰撞概率的真實值。但在PDF分布橢圓的長短軸之比較大時，誤差較大，并且由于該方法對積分域進行了近似，導致近似積分模型與原積分模型產生了無法準確衡量的偏差。

針對上述方法中的不足，本文在參考文獻[8]的基礎上，基于拉普拉斯變換的相關理論，推導了碰撞概率的冪級數表達式；研究了冪級數表達式的截斷誤差計算方法；討論了在不同精度需求下冪級數項數的取值；通過對2009年美俄衛星碰撞事件的仿真計算，將基于拉普拉斯變換方法的碰撞概率計算結果與Chan方法、Monte Carlo方法的計算結果進行對比，驗證了基于拉普拉斯變換方法在計算精度上的優勢。

1 碰撞概率計算的基本方法

1.1 碰撞概率計算假設

碰撞概率的計算不僅涉及到接近時刻(Time of Closet Approach，TCA)和接近距離(Miss Distance，MD)，還涉及到軌道的誤差協方差信息。由于2個空間目標在碰撞時所受的力比較復雜，誤差數據會根據目標運動狀態的改變發生變化，因此，為了簡化碰撞概率的計算過程，得到較準確的碰撞概率值，本文做出以下假設[9-11]：

1) 由于2個目標的相遇時間非常短，可將相遇期間內2個目標的相對運動視為速度恒定的線性運動。

2) 在碰撞時，目標的相對速度非常大，并且碰撞時間很短，因此，可以忽略2個目標速度的不確定性。

3) 2個目標的位置誤差互不關聯。

4) 在2個目標間的相對運動過程中，相遇時間非常短，因此認為誤差橢球在相遇期間內保持不變。

5) 2個空間目標均等效為半徑已知的球體。

6) 由于2個目標的誤差都是隨機的，可以認為2個目標的位置誤差都服從三維高斯分布。

當2個目標間的距離小于其等效半徑之和時，被判定為發生碰撞。

1.2 相遇坐標系與位置誤差的投影

1.2.1 相遇坐標系

令在TCA處，與2個目標的相對速度矢量垂直并且經過相對距離矢量的平面為相遇平面。在相遇平面內建立對應的相遇坐標系，即可消去速度方向的誤差并且進行碰撞概率計算。相遇坐標系以主目標O1為坐標原點，xe軸指向從目標O2在相遇平面內的投影點，ze軸指向2個目標間相對速度的方向，ye軸在相遇平面內垂直于xe軸,相遇坐標系如圖1所示。

3個軸的單位矢量可表示為

(1)

式中：ρTCA和vrel分別為TCA處目標間的相對距離矢量和相對速度矢量。

1.2.2 位置誤差的投影

通過誤差分析方法即可得到2個空間目標的誤差協方差矩陣。將2個目標的誤差橢球、相對運動狀態矢量分別投影到相遇平面內，由于2個目標的隨機位置矢量相互獨立且滿足三維高斯分布，2個目標在相遇平面內的位置矢量滿足三維高斯分布，相對位置矢量也是滿足高斯分布的隨機矢量。因此，在接近時刻處，令2個空間目標在相遇平面內的位置矢量分別為re1和re2，則相遇平面內的目標間相對位置矢量rerel的均值為

E(rerel)=E(re1-re2)=[ρTCA0]T

(2)

目標間相對位置矢量的誤差方差矩陣為

var(rerel)=var(re1-re2)=

var(re1)+var(re2)-2cov(re1,re2)=

var(re1)+var(re2)

(3)

通過式(2)和式(3)，即可將2個目標的位置誤差橢球投影到相遇平面內并形成聯合誤差橢圓域。若將2個目標看作是半徑分別為R1和R2的均勻球體，則可以根據2個目標的半徑大小形成聯合球體，將聯合球體投影到相遇平面內形成聯合圓域，圓域的有效半徑為R=R1+R2。圖2為相遇坐標系內的聯合圓域和聯合誤差橢圓域示意圖。

圖2 聯合圓域和聯合誤差橢圓域示意圖Fig.2 Schematic diagram of combined disk and combined error ellipse

1.3 碰撞概率計算的簡化

在定義了相遇坐標系后，碰撞概率的計算就被簡化為一個對PDF的二重積分過程，PDF表示為

(4)

式中:μx和μy分別為聯合誤差橢圓沿x和y方向的期望值；σx和σy分別為聯合誤差橢圓沿x和y方向的標準差。

碰撞概率Pc的計算式為

(5)

2 基于拉普拉斯變換的碰撞概率計算模型

2.1 計算流程

基于拉普拉斯變換的碰撞概率計算主要有3個步驟：

1) 將由式(5)表示的二維PDF積分公式重寫為g(λ)，其中λ=R2，函數g(λ)的拉普拉斯變換Lg在閉合域內進行。將Lg進行泰勒展開，然后使用拉普拉斯逆變換將原始的積分函數寫成冪級數的形式。

2) 計算得到冪級數的相關系數。

3) 得到碰撞概率計算結果。

從數值計算的角度考慮，直接計算g(λ)值是比較困難的。即使g(λ)的級數展開式是收斂的，但是由于每一項的大小接近、符號相反，會造成在一定精度下，級數的每一項加起來的和為0，結果為無效數字。在該情況下，對于較大數值的λ而言，級數展開是沒有意義的[12]。因此，在運算前應設置一個函數ψ，使用ψ·g替換g，以避免冪級數包含無用項，并且防止經過展開后的冪級數項在相加后出現高消去現象。

2.2 拉普拉斯變換與冪級數表達式

2.2.1 拉普拉斯變換

在進行拉普拉斯變換之前將碰撞概率改寫為Pc=g(R2)的形式，函數g:R+R+定義為

(6)

文獻[13]闡述了選擇ψ的方法，從拉普拉斯變換的角度考慮，本文使用的函數ψ的形式為ψ=epλ，其中p為待定因子。

2.1節因為把原二重積分公式重寫為了函數g(λ)，所以若要求解原積分公式的拉普拉斯變換，首先需求得函數g(λ)的拉普拉斯變換。由拉普拉斯變換定理可得函數g(λ)的拉普拉斯變換為

(7)

式中：s為復數變量。

根據富比尼定理[14]：

(8)

以及當f(x,y)=a(x)b(y)時有

(9)

可得Lg(s)的最終表達式為

(10)

定義函數h:R+R+為如下形式：

h(λ)=epλg(λ)

(11)

式中：p∈R+。則h的拉普拉斯變換Lh可表示為

(12)

式中：s>p，將式(12)轉化為

(13)

式中：η、Qx、Qy、a0各項被定義為

(14)

為了得到積分公式的冪級數表達式，可將Lh(s)進行泰勒展開，然后對該泰勒展開式的每一項進行拉普拉斯逆變換，再通過一系列計算得到積分公式的冪級數表達式。

2.2.2 冪級數表達式

令Lh(s-1)的值定義為

(15)

定義Hh(s)=s-2Lh(s-1)，則Hh(s)的導數為

(16)

由Hh(0)=a0以及式(16)可知，Hh(s)在復平面內是復可微的，根據冪級數的定義[14]，Hh(s)可展開為冪級數的形式。冪級數的求解復雜，因此使用了Maple數學工具(世界上最為通用的數學工具之一，如MATLAB、Mathematica等，可在www.maplesoft.com下載)對其進行系數的求值。根據Hh(s)=s-2Lh(s-1)，可知Lh(s)的冪級數形式為

(17)

式中:ak為系數。將式(17)中的每一項進行拉普拉斯逆變換，即可得到積分公式的冪級數形式。根據有理函數拉普拉斯變換與逆變換列表中的

(18)

式中：n=0,1,…。碰撞概率積分公式的冪級數形式可表示為

Pc=g(λ)=h(λ)e-pλ=

(19)

式中：

(20)

根據文獻[13]中級數展開的定義以及文獻[8]，函數ψ中的因子p可表示為

(21)

3 截斷誤差與冪級數項數的確定

3.1 截斷誤差

(22)

(23)

通過計算可以得到以下關系：

(24)

(25)

3.2 冪級數項數的確定

(26)

由于斯特林公式(Stirling’s Approximation)對n！的定義為

(27)

(28)

(29)

則可得到如下關系式：

(30)

令n+1=n1，由于n1≥1，則式(30)中的最后一個不等式可寫為

(31)

式中：A=η/2+(Qx+Qy)/p。通過計算，冪級數項數n的表達式為

(32)

若將Dd的值取為10-20～100，則n的計算結果如圖3所示。

若將2個目標的聯合有效半徑R的值取設為5、10、15、50、100、200、500和1000 m，則n在不同聯合有效半徑以及不同概率門限值下的計算結果如圖4所示。

從圖3和圖4結果可以看出，使用基于拉普拉斯變換的碰撞概率方法可以根據所需門限值(也作為計算精度的需求)的不同，得到不同的冪級數項數。從圖4可知，當R在1 000 m以內、計算精度需求為0～「-lg10-5?時，冪級數項數n不超過7，隨著計算精度的需求增大，項數逐漸增大。當R為1 000 m、計算精度需求為0～「-lg10-5?時，冪級數項數n不超過13。國際上通用的規避機動概率黃色門限值為10-5，碰撞概率的計算精度可設置為「-lg10-5?，因此，基于拉普拉斯的變換方法明顯減少了計算量，同時保證了碰撞概率的計算精度。

圖3 冪級數項數隨概率門限值的變化Fig.3 Variation of number of terms of power serieswith probability threshold

圖4 冪級數項數隨概率門限值以及聯合圓域有效半徑的變化Fig.4 Variation of number of terms of power series with probability threshold and radius of combined disk

4 碰撞實例仿真驗證

本文以2009年2月10日，美俄衛星發生的碰撞實例作為仿真算例。其中，美國衛星為銥星公司的Iridium-33(ID:24946)，俄羅斯失效衛星為Cosmos-2251(ID:22675)。通過接近分析得到，TCA[5]為2009年2月10日16:55:59.8 094 031，MD[5]為0.578 479 931 6 km。美俄衛星在TCA的運動狀態如表1和表2所示。

表1 主、從目標的位置參數Table 1 Position parameters of primary andsecondary object km

表2 主、從目標的速度參數Table 2 Velocity parameters of primary andsecondary object km/s

4.1 不同誤差橢球分布下的碰撞概率變化

令2個空間目標的半徑都為5 m，假設2個目標的位置誤差橢球大小相等，誤差橢球沿星基坐標系RSW的坐標軸R、S、W方向的比值σR:σS:σW不同時，Chan方法與基于拉普拉斯變換方法計算出的概率隨σS變化的情況如圖5所示。

從圖5可以看出，誤差值對碰撞概率的計算結果影響較大。使用基于拉普拉斯變換方法的計算結果曲線與Chan方法的計算結果曲線基本重合，初步證明了基于拉普拉斯變換的碰撞概率計算方法是正確的并且符合初步的碰撞預警的精度需求。

4.2 不同方法的碰撞概率計算結果比較

通過誤差分析，美俄衛星R、S、W方向的位置誤差標準差如表3所示。

令2個目標的半徑都為5 m，通過計算可得，在TCA處，2個目標在相遇平面內的3σ聯合誤差橢圓和碰撞圓域，如圖6所示。

對于Chan方法，取無窮級數的前10項以近似PDF積分式的值Pk，計算過程中，碰撞概率無窮級數形式的前10項結果如表4所示。

當使用基于拉普拉斯變換的方法對2個目標進行碰撞概率計算時，將冪級數的項數取值為k=1,2,…,10，則得到的碰撞概率值Pk和截斷誤差Sk如表5所示。

當無窮級數和冪級數項數分別取10時，Monte Carlo方法、Chan方法以及使用基于拉普拉斯變換方法的碰撞概率計算結果如表6所示。

從表6可以看出，基于拉普拉斯變換方法與Chan方法、Monte Carlo方法所計算得到的Pc都大于紅色門限值，說明在該接近點發生碰撞的可能性非常大，進一步證明了基于拉普拉斯變換的碰撞概率計算方法是正確有效的。表7列出了Chan方法、基于拉普拉斯變換方法所計算的碰撞概率與Monte Carlo方法所計算的碰撞概率的相對誤差。

圖5 不同球形分布下碰撞概率的變化趨勢Fig.5 Change trend of collision probability under different spherical distribution

目 標σR/kmσS/kmσW/km主目標0．2312070．20618850．071975從目標0．03632340．41020690．034114

圖6 3σ聯合誤差橢圓和碰撞圓域Fig.6 3σ combined error ellipse and collision disk

kPk/10-411．81743946141183421．82846608233051131．82848838960784341．82848841217507251．82848841218878661．82848841218880171．82848841218880681．82848841218880891．828488412188809101．828488412188809

表5 基于拉普拉斯變換方法下Pk和Sk的計算結果Table 5 Calculation results of Pk and Sk based onLaplace transformation method

基于拉普拉斯變換方法的碰撞概率與Chan方法、Monte Carlo方法的相對誤差分別為0.625 4%、0.078 9%，而Chan方法與Monte Carlo方法的相對誤差為0.549 9%。結果表明，基于拉普拉斯變換的碰撞概率計算精度不僅能夠滿足風險評估的要求，而且在項數取10時比Chan方法的計算精度高。

表6 Monte Carlo方法、Chan方法、基于拉普拉斯變換方法的碰撞概率計算結果Table 6 Collision probability calculation results ofMonte Carlo, Chan and based on Laplacetransformation methods

表7 Chan方法、基于拉普拉斯變換方法與Monte Carlo方法的碰撞概率相對誤差Table 7 Relative error of collision probability betweenMonte Carlo method and Chan, based on Laplacetransformation method

同時，隨著冪級數項數n的取值增大，基于拉普拉斯變換方法計算的碰撞概率的截斷誤差隨之減小。當n的取值小于5時，每一個項數所對應的碰撞概率結果相差較大，而在n大于5時，其結果基本沒有變化。因此，對于實例1，冪級數項數n取5，即可滿足初步的碰撞預警要求。需要注意的是，由于2種方法的模型不同，當Chan方法和拉普拉斯變換方法的級數同時取首項時，Chan方法的計算結果精度要比拉普拉斯方法高，但卻不滿足碰撞風險評估的要求，所以為了盡可能的提高計算精度，項數n的取值應選擇較大的值。

4.3 計算精度需求對冪級數項數的影響分析

為了分析計算精度需求對冪級數項數n及碰撞概率的影響，現將「-lgDd?取4、5、6三個值，計算對應的項數n、碰撞概率以及截斷誤差。表8列出了「-lgDd?取值不同時，對應的冪級數項數n的取值、相應的截斷誤差以及碰撞概率結果。

從表8可以看出，隨著計算精度要求不斷提高，冪級數項數n的值不斷增大，截斷誤差不斷減小，碰撞概率值也更加精確。可以預見，只需將冪級數項數設置為較大的值，碰撞概率的計算結果將會足夠精確，但是為了提高計算速度，可根據黃色門限值大小來確定精度需求，從而計算冪級數項數的值。

表8 精確位數-lg Dd不同時所需的冪級數項數、碰撞概率和截斷誤差Table 8 Number of terms of power series, collision probability and truncation error with different exact digits -lg Dd

4.4 基于拉普拉斯變換方法與Chan方法的計算時間比較

Monte Carlo方法是數值方法之一，由于數值方法沒有簡化，所以精度高，但是計算速度慢。因此本節只對解析方法中的Chan方法和拉普拉斯變換方法的計算時間進行實驗比較，以驗證拉普拉斯變換方法在計算速度上的優勢。

依舊以美俄衛星碰撞實例為算例，以表1～表3中的運動狀態與誤差數據作為2個目標的參數數據，在實驗環境相同的條件下，當冪級數項數n取不同的值時，Chan方法和拉普拉斯變換方法所花費的計算時間比較結果如圖7所示。

圖7 Chan方法和拉普拉斯變換方法的計算時間結果Fig.7 Results of computation time of Chan and Laplace transformation methods

通過計算，充分證明了基于拉普拉斯的變換方法與Chan方法相比，在速度上有所提高。基于拉普拉斯變換的方法是對原積分模型直接進行拉普拉斯變換和泰勒展開，從而獲得冪級數表達式的碰撞概率計算方法，雖然沒有對積分域進行近似，但卻是一種解析方法，因此兼具了精度和速度方面的優勢。

5 結 論

基于拉普拉斯變換的碰撞概率計算方法是一種解析方法，該方法避免了對二維積分近似過程中的誤差，因此，兼具計算精度和速度上的優勢。

1) 基于拉普拉斯變換的碰撞概率計算方法將概率密度函數的積分公式轉換到復平面內進行計算，再通過拉普拉斯逆變換并使用Maple數學工具將原積分公式表示為冪級數形式，這使得原本在實數域中計算較復雜的問題變得簡單和易理解。

2) 通過對冪級數項數的確定，可以在保證計算精度的條件下，提高碰撞概率的計算速度。

3) 針對美俄衛星碰撞實例，將基于拉普拉斯變換方法的碰撞概率計算結果與Chan方法和高精度的基于Monte Carlo數值方法的概率計算結果進行比較。結果表明，基于拉普拉斯變換的碰撞概率計算方法與Chan方法相比，能夠避免積分域的近似，在計算精度上有所提高，能夠滿足碰撞預警的精度要求。

參考文獻 (References)

[1] ALFANO S.Satellite conjunction Monte Carlo analysis[J].Advances in the Astronautical Sciences,2009,134:2007-2024.

[2] FOSTER J L,ESTES H S.A parametric analysis of orbital debris collision probability and maneuver rate for space vehicles:NASA/JSC-25898[R].Houston:NASA Johnson Space Flight Center,1992.

[3] PATERA R P.General method for calculating satellite collision probability[J].Journal of Guidance,Control,and Dynamics,2001,24(4):716-722.

[4] ALFANO S.Satellite collision probability enhancements[J].Journal of Guidance,Control,and Dynamics,2006,29(3):588-592.

[5] CHAN F K.Collision probability analysis for earth orbiting satellites[J].Advances in the Astronautically Sciences,1997(96):1033-1048.

[6] CHAN F K.Spacecraft collision probability[M].El Segundo,CA:Aerospace Press,2008.

[7] 陳磊,韓蕾,白顯宗,等.空間目標軌道力學與誤差分析[M].長沙:國防科技大學出版社,2010:178-180.

CHEN L,HAN L,BAI X Z,et al.Orbit target orbit mechanics and error analysis[J].Changsha：National University of Defense Technology Press,2010:178-180(in Chinese).

[8] SERRA R,ARZELIER D,JOLDES M,et al.Fast and accurate computation of orbital collision probability for short-term encounters[J].Journal of Guidance,Control,and Dynamics,2016,39(5):1-13.

[9] ALFRIEND K T.AKELLAM R,FRISBEE J,et al.Probability of collision error analysis[J].Space Debris,1999,1(1):21-35.

[10] AKELLA M R,ALFRIEND K T.Probability of collision between space objects[J].Journal of Guidance,Control,and Dynamics,2000,23(5):769-772.

[11] ALFANO S.A numerical implementation of spherical object collision probability[J].Journal of the Astronautical Sciences,2005,53(1):103-109.

[12] CHEVILLARD S,MEZZAROBBA M.Multiple-precision evaluation of the Airy Ai function with reduced cancellation[C]∥21st IEEE Symposium on Computer Arithmetic (ARITH).Piscataway,NJ:IEEE Press,2013:175-182.

[13] GAWRONSKI W,MüLLER J,REINHARD M.Reduced cancellation in the evaluation of entire functions and applications to the error function[J].SIAM Journal on Numerical Analysis,2007,45(6):2564-2576.

[14] HACKBUSCH W,SCHWARZ H R.Teubner-taschenbuch der mathematik[M].Berlin:Springer,2013:595.

[15] 李甲龍，熊建寧，許曉麗，等.碰撞風險評估標準適用性分析[J].天文學報，2014，55(5)：404-414.

LI J L,XIONG J N,XU X L,et al.A research on adaptability of collision criteria[J].Acta Astronomica Sinica,2014,55(5):404-414(in Chinese)