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

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

(中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽(yáng) 621000)

0 引 言

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

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

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

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

1.1 碰撞概率計(jì)算的提出

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

(1)

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

1.2 等效矩形域法計(jì)算碰撞概率

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

(3)

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

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

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

(4)

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

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

(5)

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

(6)

(7)

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

(8)

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

2 算例分析

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

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

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

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

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

(9)

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

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

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

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

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

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

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

表2 計(jì)算時(shí)長(zhǎng)Table 2 The calculation time s

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

3 結(jié) 論

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

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