空間目標(biāo)碰撞概率的等效矩形域快速算法

董文樸，任磊生，周 浩，陳 鴻，蘭勝威，李 毅

(中國空氣動力研究與發(fā)展中心，綿陽 621000)

0 引 言

近年來空間碎片問題日益嚴(yán)峻，給空間資源利用的可持續(xù)發(fā)展帶來了極大的阻礙。截止到2018年4月，地球外層空間編目的10 cm以上空間物體多達(dá)18900，其中空間碎片占到60%以上[1]。空間碎片不斷增加的主要因素有航天器碰撞或爆炸解體過程和空間碎片相互碰撞解體的反饋連鎖效應(yīng)等[2-3]。為了使空間碎片切實可控，已經(jīng)開展大量空間目標(biāo)的碰撞解體和碰撞預(yù)警研究，用于碎片減緩、防護設(shè)計、監(jiān)測預(yù)警等目的[4-6]。

依照解體模型[7]，空間目標(biāo)解體過程產(chǎn)生的碎片軌跡可以進行預(yù)報分析。如果能夠?qū)⑺槠壽E精確預(yù)測，則航天器可以提前采取規(guī)避措施，避免碰撞風(fēng)險。但是解體模型和軌跡預(yù)報會引入一定誤差，這使得碰撞變成概率事件。計算空間目標(biāo)碰撞概率成為風(fēng)險評估的一個基礎(chǔ)問題。在空間兩目標(biāo)相互靠近的過程中，根據(jù)軌道預(yù)報模型可以對交會時刻狀態(tài)進行預(yù)報，預(yù)報的精確度和初始誤差的傳播決定了交會時刻的誤差分布[8-9]。一般情況下空間目標(biāo)相交會的速度較大，交會瞬間可以近似為瞬時過程。在這種情況下，假設(shè)目標(biāo)在相遇期間速度保持不變，位置誤差保持不變，交會目標(biāo)的運動近似為線性相對運動。在線性相對運動假設(shè)下，交會目標(biāo)的碰撞概率計算可以轉(zhuǎn)換到交會平面上去求解，則把三維問題轉(zhuǎn)換成二維平面問題進行計算，使得計算問題得到簡化，提高了碰撞概率計算速度[10]。

在空間目標(biāo)碰撞概率預(yù)警工程中，要面對大量的碎片碰撞概率計算和預(yù)測，所以交會平面二維積分的計算速度直接影響了預(yù)警速率。針對空間目標(biāo)交會平面的碰撞概率計算問題，王華等[11]研究了三種簡化方法，一是將概率密度函數(shù)等效成均勻密度分布進行簡化計算；二是將二維積分問題轉(zhuǎn)化為一維曲線積分[12]；三是數(shù)值求解。這三種方法各有優(yōu)缺點，但是在精度和計算速度上不能同時滿足要求。安喜彬等[13]提出梯度數(shù)值積分方法，將橢圓區(qū)域積分轉(zhuǎn)化為適量的微元進行求和計算，減少了數(shù)值積分計算密度，但不利于解析分析問題。白顯宗等[15-16]在Chan[14]提出的等效圓域法基礎(chǔ)上，提出空間壓縮無窮級數(shù)方法，該方法將概率積分轉(zhuǎn)化為級數(shù)解，通常選取級數(shù)首項在保證計算精度的同時，提高了計算速度，同時該算法可以對交會目標(biāo)最大碰撞概率進行評估[17]。本文基于等效圓域法的思想，提出等效矩形域法對概率積分進行近似求解。在空間目標(biāo)交會算例中顯示，等效矩形域法比空間壓縮無窮級數(shù)法有更好的穩(wěn)定性，適用于碰撞概率預(yù)警評估。

1 空間目標(biāo)碰撞概率的計算

1.1 碰撞概率計算的提出

對于相互靠近的兩個空間目標(biāo)，可以通過軌道模型預(yù)測兩目標(biāo)交會時刻，給出最接近時刻和最接近距離。具體計算中將兩目標(biāo)狀態(tài)分布合并考慮，給出在交會時刻相對位置狀態(tài)概率分布，也就是聯(lián)合概率分布[10]。對于狀態(tài)誤差為正態(tài)分布的空間目標(biāo)，兩目標(biāo)的聯(lián)合概率分布仍然為正態(tài)分布，使得聯(lián)合誤差協(xié)方差矩陣等于兩目標(biāo)協(xié)方差矩陣之和。聯(lián)合誤差分布在坐標(biāo)空間中的等概率面為橢球形，對協(xié)方差矩陣進行對角化，可以得到誤差橢球的三個半軸σ1，σ2和σ3：

(1)

聯(lián)合分布給出了兩個空間目標(biāo)在交會時刻的相對位置誤差分布。如圖1所示，在線性相對運動假設(shè)下建立交會坐標(biāo)系，其中平面xy為目標(biāo)交會平面，z軸指向相對速度方向，最接近距離為ρh。誤差橢球在交會平面xy的投影半徑為σx和σy，其中x軸為半長軸方向。由于在軌空間目標(biāo)尺寸形狀不一，為了方便處理分析，將目標(biāo)等效成相應(yīng)尺寸的球形目標(biāo)。在交會時刻，如果兩目標(biāo)質(zhì)心相對距離小于兩目標(biāo)的等效半徑之和，則認(rèn)為碰撞發(fā)生。在交會坐標(biāo)系中，利用投影將兩目標(biāo)的碰撞概率，轉(zhuǎn)化為半徑為R的圓域內(nèi)的聯(lián)合誤差積分：

(2)

式中:積分區(qū)域S為x2+y2≤R2，μx和μy分別為ρh在坐標(biāo)軸上的分量。

圖1 聯(lián)合誤差橢球在交會平面的投影示意圖

1.2 等效矩形域法計算碰撞概率

圖2 等效矩形區(qū)域

(3)

其中，Φ(x)為標(biāo)準(zhǔn)正態(tài)分布函數(shù)。在等效矩形域下，二重積分公式近似成解析形式，使得二維積分計算得到簡化。同時根據(jù)積分域范圍給出等效矩形域法的誤差范圍，從圖2可以看出，當(dāng)矩形域外接積分圓時，可以得到積分PU的近似最大值Pmax，內(nèi)接時得到最小值Pmin，則有Pmin

1.3 基于等效矩形域法的最大碰撞概率計算

隨著空間目標(biāo)的誤差橢球形狀和大小的變化，碰撞概率存在最大值。當(dāng)目標(biāo)誤差不確定時或者被沖淡時[18]，最大碰撞概率可以為空間目標(biāo)接近狀態(tài)評估提供參考。利用等效矩形域法來近似求解誤差橢球形狀比例已知的碰撞概率最大值。對碰撞概率計算式(2)可將誤差橢圓半軸作等比例變換：

(4)

圖3 積分區(qū)域S′示意圖

在二維積分PU作等比例變換后，將其轉(zhuǎn)變到極坐標(biāo)系進行求解，令：

(5)

在坐標(biāo)變換下，式(4)可改寫為：

(6)

(7)

根據(jù)式(7)對σy求導(dǎo)，可以得到：

(8)

在橢圓長短軸比例為k的情況下，當(dāng)σy=σM時，根據(jù)式(3)得到最大碰撞概率PM。

2 算例分析

選取兩個算例進行驗證：算例1來自于空間目標(biāo)編號為27430和28064在交會時刻2005年3月8日10:11:56.599LT的交會狀態(tài)參數(shù)[19-20]；算例2為假設(shè)空間兩目標(biāo)接近相撞時的交會狀態(tài)參數(shù)。算

例1和算例2的交會狀態(tài)參數(shù)在表1中給出。

表1 交會時刻狀態(tài)參數(shù)Table 1 The position parameters at the encounter time m

圖4 算例1碰撞概率隨系數(shù)k的變化

對于算例1，假設(shè)聯(lián)合誤差橢球在交會平面的投影均方差為σy=20000 m，σx由第1.3節(jié)中的比例系數(shù)k給出。在給定狀態(tài)下，空間目標(biāo)碰撞概率隨系數(shù)k的變化結(jié)果如圖4所示。其中，碰撞概率積分式(2)的數(shù)值積分結(jié)果如PU所示，PE為利用等效矩形域方法近似求解結(jié)果。為了對比等效矩形域方法，利用基于空間壓縮無窮級數(shù)快速方法計算碰撞概率積分[15-16]，在該方法下，考慮計算效率以及收斂性，通常選取首項進行近似分析，根據(jù)首項得到的概率積分式(2)近似為：

(9)

利用式(9)計算得到的結(jié)果用PC表示。

圖4(a)中顯示，對于算例1中的目標(biāo)交會情況，空間目標(biāo)在交會時刻的相對距離較遠(yuǎn)，碰撞概率整體處在10-7量級，同時給出不同方法的估計范圍。隨著k的增加，碰撞概率繼續(xù)下降。兩種近似方法與數(shù)值積分法得到的曲線相一致，計算結(jié)果基本重合。為了進一步對比不同方法的近似程度，取近似值與數(shù)值積分值的相對偏差，將相對偏差與數(shù)值積分值作比，得到相對比例偏差曲線，如圖4(b)所示。PE和PC的近似結(jié)果與數(shù)值積分結(jié)果的相對偏差在同一個數(shù)量級。但是，當(dāng)系數(shù)k增大時，也就是誤差橢球被拉長時，PC的近似值偏差增大，而PE近似偏差相對較為平穩(wěn)。

在圖4(c)中，給出了最大碰撞概率隨均方差比例k的變化。在算例1交會情況下，兩種近似方法估計的最大碰撞概率結(jié)果都與數(shù)值方法相吻合。在圖4(d)中給出對應(yīng)的相對比例偏差結(jié)果，顯示出與碰撞概率相對比例偏差相類似的趨勢，估計偏差都在10-5量級，曲線PC的最大碰撞概率估計值隨k的增長偏差增大，而等效矩形域法給出的結(jié)果PE相對穩(wěn)定不變。

對于算例2，兩個空間目標(biāo)相互靠近，在交會時刻接近距離為100 m，目標(biāo)的聯(lián)合半徑為10 m，兩目標(biāo)接近碰撞，此時假設(shè)σy=10 m，碰撞概率隨系數(shù)k的變化情況如圖5(a)所示。對于算例2，可以看到空間壓縮無窮級數(shù)法PC給出的近似結(jié)果與數(shù)值方法PU偏差較大。在圖4(b)中給出了PC與PU對應(yīng)的相對比例偏差結(jié)果，可以看到在系數(shù)k接近30時，相對偏差達(dá)到9%。對于等效矩形域法，估計結(jié)果PE隨k的變化仍然與數(shù)值結(jié)果PU保持較小的偏差，并且能夠隨k的增長保持穩(wěn)定。

在圖5(c)和圖5(d)中給出了兩種方法對算例2的最大碰撞概率估計。在算例2情況下，最大碰撞概率PU隨k的增長持續(xù)增加。等效矩形域法和空間壓縮無窮級數(shù)法給出的最大碰撞概率估計結(jié)果PE和PC，顯示出相同的趨勢。系數(shù)k在5～25范圍內(nèi)時，PE與PC作對比分析，等效矩形域法PE與數(shù)值計算結(jié)果PU更為接近，等效矩形域法的估算基本在1%以內(nèi)。

圖5 算例2碰撞概率隨系數(shù)k的變化

結(jié)合算例1和算例2對不同碰撞概率近似方法進行計算效率分析。在計算碰撞概率過程中，等效矩形域法(式(3))中包含正態(tài)分布函數(shù)的積分運算。正態(tài)分布函數(shù)被廣泛應(yīng)用于工程計算中，在算法和應(yīng)用上較為成熟，其函數(shù)在Cmath庫和GSL等科學(xué)計算庫中都有調(diào)用。為了分析碰撞概率計算效率，在具體采用了誤差函數(shù)展開法[21]和調(diào)用GSL庫函數(shù)[22]兩種方法進行等效矩形域法算法實現(xiàn)，并與空間壓縮無窮級數(shù)法和數(shù)值積分法進行對比，其中算法實現(xiàn)語言為C++。表格2給出了算例1和算例2在比例系數(shù)k為5時不同方法的計算時長，其中計算結(jié)果保證4位有效數(shù)字的精度。

表2 計算時長Table 2 The calculation time s

表2給出了不同方法計算碰撞概率的時長，其中，PU為數(shù)值積分方法，PE-GSL和PE-erf分別為采用GSL庫函數(shù)和誤差函數(shù)展開方法求解的等效矩形域方法，PC為空間壓縮無窮級數(shù)法，計算時長為重復(fù)計算107次的平均CPU耗時，不同計算平臺和編譯環(huán)境的絕對時長可能會有差異，所以這里只作相對分析。計算結(jié)果顯示，PE和PC比PU有兩個數(shù)量級以上的優(yōu)勢。PC用時最短，PE-GSL和PE-erf時長為PC的1～4倍左右。計算結(jié)果表明，等效矩形域方法可以達(dá)到與空間壓縮無窮級數(shù)快速算法相同量級的計算速度，可以利用其實現(xiàn)對碰撞概率的快速預(yù)測。

3 結(jié) 論

本文基于誤差橢球交會平面投影算法，對接近空間目標(biāo)碰撞概率的計算方法進行探討，提出了等效矩形域快速計算方法。等效矩形域方法是將交會平面內(nèi)的碰撞概率積分域用等效矩形域代替，簡化了二維積分，得到了形式更簡單的碰撞概率解析近似式。在空間目標(biāo)誤差橢球形狀比例已知的情況下，利用等效矩形域方法推導(dǎo)了空間目標(biāo)的最大碰撞概率計算方法，可用于快速評估相互接近的空間目標(biāo)危險程度。

本文選取兩個算例進行計算分析，并與空間壓縮無窮級數(shù)近似法相對比，計算結(jié)果顯示了等效矩形域方法快速估計概率積分的可靠性。同時，等效矩形域方法可以適應(yīng)不同的聯(lián)合誤差橢球半軸比例情況，并且與數(shù)值結(jié)果的偏差隨半軸比例增長變化緩慢。特別在空間目標(biāo)相互靠近的情況下，本方法的估計精度要優(yōu)于空間壓縮無窮級數(shù)方法，在碰撞預(yù)警評估中，有助于對大量靠近的相遇實例進行碰撞概率快速對比評估。