王曉偉，劉 靜，崔雙星

(1. 中國(guó)科學(xué)院國(guó)家天文臺(tái)，北京 100101；2. 國(guó)家航天局空間碎片監(jiān)測(cè)與應(yīng)用中心，北京 100101； 3. 中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引 言

近年來(lái)，空間碎片數(shù)量激增使得在軌航天器的運(yùn)行安全受到極大威脅，Kessler雪崩效應(yīng)再次被提及[1]。空間碎片數(shù)量不斷增長(zhǎng)引發(fā)了國(guó)際上關(guān)于空間碎片環(huán)境長(zhǎng)期穩(wěn)定性的研究，空間碎片環(huán)境長(zhǎng)期演化模型成為該領(lǐng)域的國(guó)際研究熱點(diǎn)之一[2-9]。

空間碎片環(huán)境長(zhǎng)期演化模型可以預(yù)測(cè)未來(lái)幾十年至上百年的空間碎片環(huán)境演化，為調(diào)整太空發(fā)展戰(zhàn)略或制定相關(guān)空間政策提供技術(shù)支持，以保證未來(lái)太空活動(dòng)可持續(xù)發(fā)展?？臻g碎片演化模型能夠模擬空間碎片的主要增長(zhǎng)機(jī)制和減少機(jī)制，例如未來(lái)的航天器發(fā)射、在軌碰撞或爆炸解體、自然隕落、任務(wù)后處置等，通常由軌道預(yù)報(bào)模型、未來(lái)發(fā)射模型、碰撞概率評(píng)估模型、解體模型、任務(wù)后處置模型等幾個(gè)子模型組成。由于在軌碰撞解體是空間碎片在未來(lái)演化中的重要增長(zhǎng)來(lái)源，對(duì)碰撞概率的評(píng)估會(huì)影響空間碎片在軌碰撞事件的預(yù)測(cè)，因此碰撞概率評(píng)估模型是空間碎片演化模型中較為關(guān)鍵的環(huán)節(jié)之一。

傳統(tǒng)的瞬時(shí)碰撞概率計(jì)算方法多是基于空間物體的位置橢球誤差[10-11]，然而空間碎片演化模型預(yù)測(cè)的時(shí)間較長(zhǎng)(幾十年甚至幾百年)，空間物體軌道預(yù)報(bào)的精度有限，且高精度的位置誤差對(duì)于長(zhǎng)期演化并不具有實(shí)際意義，因此傳統(tǒng)的瞬時(shí)碰撞概率計(jì)算方法并不適用。要預(yù)測(cè)空間碎片在未來(lái)長(zhǎng)期演化中發(fā)生碰撞解體事件的可能性，需要一種適用于長(zhǎng)期軌道演化系統(tǒng)的碰撞概率估計(jì)方法。

Cube模型是由美國(guó)國(guó)家宇航局(National aeronautics and space administration, NASA)提出的一種快速成對(duì)算法，適用于任何軌道演化系統(tǒng)。Cube算法通過(guò)對(duì)整個(gè)演化系統(tǒng)進(jìn)行時(shí)間均勻采樣，能夠利用演化過(guò)程中不斷更新的軌道根數(shù)來(lái)評(píng)估空間物體之間的碰撞概率[12-13]。Cube模型算法具有原理簡(jiǎn)單、配對(duì)快速的優(yōu)點(diǎn)。Cube模型算法提出后，被多個(gè)國(guó)家的空間碎片環(huán)境長(zhǎng)期演化模型所采納。但在2014年第32屆機(jī)構(gòu)間空間碎片協(xié)調(diào)委員會(huì)(Inter-agency space debris coordination committee, IADC)全體會(huì)議上，法國(guó)宇航局(Centre national d’etudes spatiales, CNES)指出，利用其演化模型MEDEE模擬低地球軌道(Low earth orbit, LEO)空間環(huán)境的未來(lái)演化，當(dāng)Cube模型的立方體尺寸取值不同時(shí)，空間環(huán)境演化結(jié)果會(huì)有較大差異，如圖1所示，圖中粗線表示40次蒙特卡洛運(yùn)行的平均結(jié)果，細(xì)線表示1σ的標(biāo)準(zhǔn)差。表1為MEDEE模型預(yù)測(cè)的200年后的LEO空間碎片數(shù)量相比于初始碎片數(shù)量的變化百分比，及其1σ標(biāo)準(zhǔn)差[14-15]。但從物理規(guī)律而言，空間碎片在演化過(guò)程中的碰撞概率不會(huì)因碰撞概率評(píng)估算法不同而改變，因此若Cube模型因其參數(shù)設(shè)置不同而對(duì)未來(lái)空間環(huán)境演化結(jié)果帶來(lái)顯著影響，將直接導(dǎo)致空間碎片演化模型的可信度降低。

圖1 立方體大小不同時(shí)MEDEE模型的運(yùn)行結(jié)果Fig.1 Evolution results of MEDEE with different cube sizes

立方體尺寸/km200年后碎片數(shù)量相比于初始碎片數(shù)量的百分比[平均值+/-1σ標(biāo)準(zhǔn)差]h=5-13%+/-12%h=102%+/-16%h=15-5.3%+/-14%

本文對(duì)此問(wèn)題進(jìn)行了深入分析，研究了Cube模型算法及其在空間碎片演化模型中的應(yīng)用問(wèn)題，并在Cube模型算法的基礎(chǔ)上做了改進(jìn)，提出了I-Cube模型。經(jīng)過(guò)多次蒙特卡洛模擬運(yùn)行驗(yàn)證，使用改進(jìn)后的I-Cube算法，空間碎片演化模型的演化結(jié)果不再受碰撞概率算法參數(shù)的影響，顯著提高了演化模型的穩(wěn)定性與可信度。

本文所使用的SOLEM模型(Space objects long-term evolution model)是我國(guó)自主建立的空間碎片長(zhǎng)期演化模型[16-17]。SOLEM模型作為中國(guó)國(guó)家航天局(China national space administration, CNSA)的代表參與了IADC組織的多項(xiàng)國(guó)際聯(lián)合研究，并取得了與國(guó)際上其他演化模型較為一致的結(jié)果[18]。

本文第1部分對(duì)Cube模型算法及其在空間碎片演化模型中的應(yīng)用做了進(jìn)一步說(shuō)明與分析，第2部分介紹改進(jìn)的碰撞概率模型I-Cube，第3部分對(duì)比分析了采用Cube模型算法與I-Cube算法的演化結(jié)果，第4部分進(jìn)行了總結(jié)。

1 碰撞概率模型及其對(duì)演化的影響

碰撞概率的計(jì)算直接影響空間碎片碰撞解體次數(shù)的預(yù)測(cè)。碰撞解體是未來(lái)空間碎片數(shù)量增長(zhǎng)的重要來(lái)源，因而對(duì)空間碎片間的碰撞概率估計(jì)成為影響未來(lái)空間環(huán)境演化仿真結(jié)果的關(guān)鍵環(huán)節(jié)之一。

1.1 碰撞概率算法Cube模型

Cube模型是由NASA提出的適用于任何軌道演化系統(tǒng)的碰撞概率評(píng)估算法[12-13]。Cube算法通過(guò)對(duì)整個(gè)演化系統(tǒng)進(jìn)行時(shí)間均勻采樣的方法來(lái)計(jì)算空間碎片的碰撞概率。數(shù)學(xué)上，物體i和j之間在很長(zhǎng)一段時(shí)間內(nèi)(從tbegin到tend)的總碰撞次數(shù)可表達(dá)為：

(1)

式中：Pi,j是碰撞率;L是tbegin和tend之間的時(shí)間間隔數(shù);ts和ts+1代表第s個(gè)時(shí)間間隔的起止時(shí)刻。如果時(shí)間間隔ts+1-ts足夠短，兩個(gè)空間碎片之間的碰撞特征變化不大，那么可認(rèn)為Pi,j在這段時(shí)間內(nèi)是常數(shù)，則上述積分表達(dá)式可寫為

(2)

在每次采樣時(shí)刻，建立地心笛卡爾坐標(biāo)系，將三維近地空間劃分成多個(gè)邊長(zhǎng)為h的小立方體，每個(gè)空間物體的位置、速度根據(jù)該時(shí)刻的軌道根數(shù)計(jì)算。當(dāng)兩個(gè)空間物體處于同一個(gè)立方體時(shí)，碰撞率Pij由下式計(jì)算

Pij=sisjAcVimpdU

(3)

式中：si和sj表示目標(biāo)i和j在該立方體內(nèi)的空間密度;Ac表示兩目標(biāo)的碰撞截面;Vimp表示兩個(gè)目標(biāo)的碰撞速度;dU表示該立方體的體積。對(duì)Pij進(jìn)行時(shí)間積分后得到碰撞概率

pij=Pij·[ts+1-ts]=sisjAcVimpdUdt

(4)

計(jì)算出pij后，應(yīng)用蒙特卡羅方法生成一個(gè)隨機(jī)數(shù)Ri與pij進(jìn)行比較，以此來(lái)確定一次碰撞是否發(fā)生。若一個(gè)立方體內(nèi)同時(shí)有兩個(gè)或多個(gè)物體，計(jì)算其兩兩之間的碰撞概率，并確定碰撞是否發(fā)生；而若一個(gè)立方體內(nèi)只有一個(gè)空間物體，則不再考慮該物體與其他物體的碰撞可能性。對(duì)于N體系統(tǒng)，該方法能夠快速有效地找出碰撞配對(duì)物體，計(jì)算時(shí)間與N成正比，而非N2。

根據(jù)標(biāo)準(zhǔn)統(tǒng)計(jì)采樣方法，采樣次數(shù)越多，即采樣時(shí)間間隔越短越好，但在空間碎片長(zhǎng)期演化模型中，需要考慮整個(gè)模型的運(yùn)行速度，不能在碰撞概率計(jì)算上花費(fèi)過(guò)多的時(shí)間。在美國(guó)的空間碎片長(zhǎng)期演化模型LEGEND中，dt默認(rèn)值為5天[13]。

而關(guān)于采樣時(shí)刻所劃分的立方體大小，Cube模型認(rèn)為在平均軌道半長(zhǎng)軸的1%及以內(nèi)都是合理的[13]。將其應(yīng)用在空間碎片演化模型中時(shí)，主要研究對(duì)象是LEO軌道的碎片，軌道半長(zhǎng)軸在6578～8578 km之間，故依據(jù)Cube模型，其所劃分的立方體尺寸在65 km及以下時(shí)都是合理的。美國(guó)空間碎片長(zhǎng)期演化模型LEGEND中，h取值為10 km[13]。

1.2 Cube算法存在的問(wèn)題

依據(jù)Cube模型，上述碰撞概率算法式(4)相當(dāng)于在一個(gè)微觀尺度上(一個(gè)立方體內(nèi))應(yīng)用氣體動(dòng)力學(xué)理論。然而，根據(jù)式(4)計(jì)算出的pij并非真正的碰撞概率，而是兩個(gè)物體在所劃分的立方體dU內(nèi)在積分時(shí)間dt內(nèi)的平均碰撞次數(shù)。在積分時(shí)間并非趨于0的情況下，對(duì)于個(gè)別截面積超大的碎片，用該式計(jì)算出的碰撞概率存在pij>1的情況。實(shí)際上，在pij≤0.2的情況下，可將其近似看作碰撞概率，其誤差小于10%，而當(dāng)pij>0.2時(shí)，該近似不再成立，必須使用碰撞概率的嚴(yán)格表達(dá)式(見(jiàn)式(5))[19]。

圖2 立方體大小h對(duì)空間碎片演化預(yù)測(cè)結(jié)果的影響Fig.2 Evolution results of SOLEM model with the cube size h

2 改進(jìn)的碰撞概率算法I-Cube模型

圖3 考慮相鄰立方體中兩個(gè)物體的碰撞概率的二維示意圖Fig.3 Two-dimensional representation for considering possible collisions between debris residing in neighbouring cubes

I-Cube算法的具體實(shí)現(xiàn)步驟如下：

1)確定對(duì)演化系統(tǒng)進(jìn)行時(shí)間均勻采樣的間隔dt。在SOLEM模型中，時(shí)間采樣間隔dt為5天。

2)對(duì)近地空間劃分立方體。在采樣時(shí)刻，建立地心笛卡爾坐標(biāo)系，將地球周圍的近地空間劃分為一個(gè)個(gè)邊長(zhǎng)為h的立方體，并對(duì)每個(gè)立方體按h的倍數(shù)進(jìn)行編號(hào)。

3)計(jì)算空間物體所在的位置。在采樣時(shí)刻，根據(jù)每個(gè)空間物體(序號(hào)為i)更新后的軌道根數(shù)計(jì)算其所在的位置，并記錄其所處的立方體編號(hào)。

5)計(jì)算配對(duì)目標(biāo)之間的碰撞概率pij。

6)利用蒙特卡羅方法判斷碰撞是否實(shí)際發(fā)生。生成0～1之間的隨機(jī)數(shù)Ri，若Ri

步驟5)中碰撞概率pij的具體算法如下：

根據(jù)氣體動(dòng)力學(xué)，在體元dU內(nèi)，在dt時(shí)間內(nèi)，空間物體i和j的平均碰撞數(shù)為

c=SiSjVimpAcdUdt

(5)

pij=1-exp(-c)

(6)

與Cube算法相比，I-Cube算法采用了嚴(yán)格表達(dá)式來(lái)計(jì)算碰撞概率pij，此外增加了處于相鄰立方體內(nèi)的空間物體間碰撞可能性的計(jì)算，尋找碰撞配對(duì)物體的復(fù)雜度增加，因而計(jì)算所用的時(shí)間也相應(yīng)增加，基本相當(dāng)于原Cube算法計(jì)算時(shí)間的3倍。但從整個(gè)演化模型運(yùn)行上看，這個(gè)計(jì)算速度仍在可接受范圍內(nèi)。

3 仿真分析

據(jù)目前公開(kāi)的文獻(xiàn)看，只有CNES利用其MEDEE模型對(duì)原Cube算法參數(shù)的影響做過(guò)相關(guān)研究[14-15]，此外，UKSA的DAMAGE模型采用Cube算法并對(duì)其做了部分改進(jìn)[20]。但由于各國(guó)的空間碎片長(zhǎng)期演化模型在各個(gè)模塊算法的具體技術(shù)實(shí)現(xiàn)方式并未公開(kāi)，MEDEE模型、DAMAGE模型以及國(guó)際上其他長(zhǎng)期演化模型在Cube算法的應(yīng)用或改進(jìn)方面，以及太陽(yáng)活動(dòng)、大氣模型等關(guān)鍵影響因素的假設(shè)與參數(shù)設(shè)置方面均與SOLEM模型不盡相同，且在已公開(kāi)的文獻(xiàn)中，MEDEE和DAMAGE模型與SOLEM模型所使用的初始空間碎片數(shù)據(jù)也不同，故無(wú)法就Cube算法改進(jìn)前后其參數(shù)對(duì)空間環(huán)境演化仿真的影響與SOLEM模型形成橫向?qū)Ρ?。鑒于SOLEM模型在國(guó)際聯(lián)合研究中的出色表現(xiàn)，并出于避免不同演化模型因模型本身差異帶來(lái)影響的考慮，本文僅利用SOLEM模型在保證其他關(guān)鍵影響因素保持不變的前提下，對(duì)Cube算法改進(jìn)前后改變參數(shù)h對(duì)演化結(jié)果的影響進(jìn)行研究。

SOLEM演化模型中將太陽(yáng)輻射流量和地磁活動(dòng)設(shè)為平均值，即F10.7=130 sfu，Ap=9。軌道預(yù)報(bào)算法為簡(jiǎn)化的半分析法，考慮地心引力、地球非球形引力攝動(dòng)、大氣阻力攝動(dòng)、第三體引力攝動(dòng)和光壓攝動(dòng)，大氣密度模型采用NRLMSIS00模型。碰撞解體模型采用NASA標(biāo)準(zhǔn)解體模型。以2017年9月1日的空間碎片環(huán)境數(shù)據(jù)作為初始輸入，預(yù)報(bào)200年至2217年9月1日。以2009.09.01至2017.08.31之間8年的發(fā)射情況作為發(fā)射模型在未來(lái)200年內(nèi)不斷循環(huán)，航天器任務(wù)壽命假設(shè)為8年，不考慮軌道維持和人工避碰措施，任務(wù)后軌道處置率假設(shè)為30%。

SOLEM模型在上述假設(shè)條件均不變的情況下，碰撞概率評(píng)估模型分別采用Cube算法和I-Cube算法，僅改變參數(shù)h的大小，研究其對(duì)未來(lái)200年空間碎片環(huán)境演化的影響。本文算例中參數(shù)h取值分別為5 km、8 km、10 km、30 km、50 km。每個(gè)算例均運(yùn)行50次蒙特卡洛模擬并求平均。仿真結(jié)果如圖4～圖7所示，其結(jié)果量化如表2所示。

圖4 采用Cube算法，未來(lái)空間碎片數(shù)量的演化Fig.4 Predictions of future space environment evolution using Cube method

圖5 采用I-Cube算法，未來(lái)空間碎片數(shù)量的演化Fig.5 Predictions of future space environment evolution using I-Cube method

圖6 采用Cube算法，未來(lái)200年空間碎片碰撞次數(shù)的預(yù)測(cè)Fig.6 Predictions of cumulative number of space debris collisions in future 200 years using Cube method

圖4～圖7以及表2所示的仿真結(jié)果表明，演化模型在其他如太陽(yáng)輻射等關(guān)鍵影響因素保持不變的情況下，采用原Cube算法，改變所劃分立方體的大小h會(huì)嚴(yán)重影響未來(lái)空間碎片環(huán)境的演化結(jié)果。就算例中h的取值范圍(5～50 km)，不同的取值造成200年后空間碎片數(shù)量演化的預(yù)測(cè)值相差約17500個(gè)，為2017年初始碎片數(shù)量的134%；200年累計(jì)災(zāi)難性碰撞次數(shù)預(yù)測(cè)相差23次，非災(zāi)難性碰撞次數(shù)預(yù)測(cè)相差約27次。由于h的取值皆在“平均軌道半長(zhǎng)軸的1%及以內(nèi)”的合理范圍，因而無(wú)法判斷哪種演化結(jié)果可信。

而采用I-Cube算法后，改變參數(shù)h大小，未來(lái)空間碎片環(huán)境的演化結(jié)果不再受其顯著影響。就算例中h的取值范圍(5～50 km)，不同的取值造成200年后空間碎片數(shù)量演化的預(yù)測(cè)值相差僅為2017年初始碎片數(shù)量的18%，200年累計(jì)災(zāi)難性碰撞次數(shù)預(yù)測(cè)相差不到4次，非災(zāi)難性碰撞次數(shù)預(yù)測(cè)相差約6次。

圖7 采用I-Cube算法，未來(lái)200年空間碎片碰撞次數(shù)的預(yù)測(cè)Fig.7 Predictions of cumulative numbers of space debris collisions in future 200 years using I-Cube method

由此可見(jiàn)，采用I-Cube算法后，未來(lái)空間碎片的演化結(jié)果高度一致，碰撞解體事件的預(yù)測(cè)不再受參數(shù)h的顯著影響。相比于采用Cube算法，I-Cube算法對(duì)碰撞概率的計(jì)算更為準(zhǔn)確合理，進(jìn)而使空間碎片演化模型的穩(wěn)定性與可信度顯著提升。

4 結(jié) 論

I-Cube算法采用了嚴(yán)格表達(dá)式來(lái)計(jì)算空間物體間的碰撞概率，排除了原Cube算法中可能出現(xiàn)“碰撞概率”大于1的隱患。

I-Cube算法與Cube算法的模型假設(shè)不同。I-Cube算法假設(shè)只要兩個(gè)物體間的距離在給定閾值范圍內(nèi)都存在碰撞可能性，而Cube算法假設(shè)只有位于同一個(gè)立方體內(nèi)的物體之間存在碰撞可能性。相比之下，I-Cube算法的模型假設(shè)更符合物理真實(shí)。

I-Cube算法保留原Cube算法尋找碰撞配對(duì)物體的思路，仍通過(guò)在采樣時(shí)刻對(duì)近地空間劃分立方體的方法來(lái)快速尋找潛在的碰撞對(duì)，但尋找潛在碰撞對(duì)的范圍由原來(lái)的同一立方體擴(kuò)展至相鄰的立方體。由于增加了尋找碰撞配對(duì)物體的復(fù)雜度，I-Cube算法犧牲了一部分計(jì)算時(shí)間，但從空間碎片演化模型的整體運(yùn)行上看，計(jì)算速度仍可接受。

表2 采用Cube算法與I-Cube算法，改變參數(shù)h對(duì)未來(lái)空間環(huán)境演化的影響對(duì)比Table 2 Impacts of varying parameter values of h on the predictions of future space debris environment evolution using Cube algorithm compared with those using I-Cube algorithm

相比于Cube算法，I-Cube算法提升了碰撞概率計(jì)算方法的準(zhǔn)確度與合理性，使空間碎片環(huán)境長(zhǎng)期演化模型的演化結(jié)果不再受自身碰撞概率算法參數(shù)的影響，從而提高了空間碎片長(zhǎng)期演化模型的穩(wěn)定性與可信度。