空間物體橢圓軌道假設(shè)下碰撞通量計(jì)算與分析

周威萍, 沈 鳴, 高鵬騏, 郭效忠, 趙 有

(1. 中國科學(xué)院國家天文臺(tái)， 北京 100012；2. 中國科學(xué)院大學(xué)， 北京 100049)

空間物體橢圓軌道假設(shè)下碰撞通量計(jì)算與分析

周威萍1,2, 沈 鳴1, 高鵬騏1, 郭效忠1, 趙 有1

(1. 中國科學(xué)院國家天文臺(tái)， 北京 100012；2. 中國科學(xué)院大學(xué)， 北京 100049)

利用陳氏模型，不依賴于氣體動(dòng)理論和泊松統(tǒng)計(jì)假設(shè)，對(duì)近地空間物體環(huán)境中環(huán)繞飛行航天器碰撞通量的計(jì)算問題進(jìn)行了研究。分析運(yùn)行在橢圓軌道上的空間物體通過指定區(qū)域的4種情況，推導(dǎo)得到碰撞通量關(guān)于物體穿越統(tǒng)計(jì)區(qū)域所用時(shí)間的表達(dá)式。文章亦對(duì)圓軌道和橢圓軌道假設(shè)作了比較，并以銥星86為例，討論了兩種假設(shè)情況下空間物體對(duì)它的碰撞通量。結(jié)果顯示，隨著統(tǒng)計(jì)區(qū)域厚度的增加，橢圓軌道假設(shè)下的碰撞通量從約為圓軌道假設(shè)的4倍減小到1.4倍；圓軌道假設(shè)下，計(jì)入統(tǒng)計(jì)的空間物體數(shù)量隨著統(tǒng)計(jì)區(qū)域的厚度變化產(chǎn)生的波動(dòng)較大，橢圓軌道假設(shè)下的統(tǒng)計(jì)目標(biāo)數(shù)量在區(qū)域厚度取5 km時(shí)可達(dá)圓軌道的60倍。

碰撞預(yù)警；碰撞通量；空間物體；航天器；衛(wèi)星軌道

0 引 言

空間碎片的存在給人造衛(wèi)星、空間站等航天器的安全運(yùn)行帶來了威脅，隨著人類航天事業(yè)的發(fā)展，空間物體的數(shù)量在不斷增加。航天器的撞擊解體、空間碎片的二次解體更是加速了碎片的增長(zhǎng)，給航天器的運(yùn)行帶來更具威脅的安全隱患[1-2]。為避免航天器因發(fā)生碰撞造成損失和空間災(zāi)難，須對(duì)航天器進(jìn)行預(yù)警分析，對(duì)將與空間物體發(fā)生碰撞的航天器進(jìn)行機(jī)動(dòng)規(guī)避[3-4]。目前，是否發(fā)布航天器碰撞預(yù)警的衡量標(biāo)準(zhǔn)多為碰撞概率閾值，F(xiàn)oster等[5]對(duì)國際空間站(International space station, ISS)的碰撞規(guī)避機(jī)動(dòng)進(jìn)行了分析，結(jié)果表明：若ISS的碰撞概率大于10-4，則應(yīng)該對(duì)其實(shí)施規(guī)避機(jī)動(dòng)策略。閾值一般根據(jù)三個(gè)量來確定：風(fēng)險(xiǎn)降低率、機(jī)動(dòng)次數(shù)、目標(biāo)的先驗(yàn)碰撞風(fēng)險(xiǎn)。作為碰撞判定基準(zhǔn)的先驗(yàn)碰撞風(fēng)險(xiǎn)理論上未知，可以用碰撞通量作為替代，碰撞通量根據(jù)背景空間物體的通量計(jì)算得到[6-7]。Kessler[8]在研究木星的八個(gè)衛(wèi)星之間的碰撞概率時(shí)，提出了Cube Approach模型[9],該模型將空間分割成一個(gè)個(gè)小空間體積元，分別計(jì)算空間物體在這些體積元內(nèi)的碰撞概率，相加得到這八個(gè)衛(wèi)星之間的總碰撞概率。根據(jù)Cube Approach模型，這些空間體積元的尺寸相比空間物體的軌道不確定度足夠小，當(dāng)空間物體穿越其中時(shí)，物體的位置可以被認(rèn)為是隨機(jī)的，這種情況下，空間物體的運(yùn)動(dòng)與氣體分子的運(yùn)動(dòng)相似，因而利用氣體動(dòng)理論模型計(jì)算空間物體之間的碰撞概率。Foster等[5]和Foster等[10]在確定航天器碰撞概率閾值的過程中，采用了Cube Approach模型計(jì)算空間物體對(duì)航天器的碰撞通量。為計(jì)算航天器在穿越空間物體云時(shí)的碰撞通量，NASA和ESA分別提出了ORDEM和MASTER模型[11-13]，同樣假設(shè)空間物體整體遵守氣體動(dòng)理論，分布上服從泊松統(tǒng)計(jì)。根據(jù)氣體動(dòng)理論，氣體分子作直線運(yùn)動(dòng)且方向隨機(jī)，它們的數(shù)密度是均勻的。然而，空間物體實(shí)際上遵循軌道運(yùn)動(dòng)規(guī)律，還受到地球引力攝動(dòng)、日月引力攝動(dòng)、大氣阻尼攝動(dòng)、太陽輻射壓力攝動(dòng)等因素的影響[14-15]。與氣體分子不同，空間物體的數(shù)密度在緯度空間內(nèi)并不均勻分布，它從赤道到碎片軌道所能到達(dá)的最大緯度不斷增大，在碎片的最高緯度，其數(shù)密度比低緯度大幾個(gè)數(shù)量級(jí)。基于這些因素，Chan[16]在1996年提出了不依賴于氣體動(dòng)理論模型和泊松統(tǒng)計(jì)假設(shè)的陳氏模型。

Dennis[17]假設(shè)空間物體的軌道為圓軌道，其升交點(diǎn)赤經(jīng)在赤道平面上均勻分布，軌道平面的相位均勻分布，從而給出了空間物體數(shù)量的統(tǒng)計(jì)學(xué)分布。Chan[16]在此基礎(chǔ)上對(duì)算法做了改進(jìn)，給出了空間物體的數(shù)密度表達(dá)式,并結(jié)合兩個(gè)物體之間的碰撞概率關(guān)系，得到不依賴于氣體動(dòng)理論模型和泊松統(tǒng)計(jì)假設(shè)的碰撞通量計(jì)算表達(dá)式，用以計(jì)算航天器在一段時(shí)間內(nèi)與多個(gè)空間物體發(fā)生交會(huì)時(shí)，可能發(fā)生的碰撞風(fēng)險(xiǎn)。然而，空間物體的運(yùn)行軌道實(shí)際上并不是絕對(duì)的圓，即使是偏心率為0.001的近圓軌道，也因?yàn)檐壍赖陌腴L(zhǎng)軸尺度而使得遠(yuǎn)、近地點(diǎn)的高度差量級(jí)為幾十千米。大橢圓軌道的區(qū)別更加明顯。在當(dāng)前設(shè)備探測(cè)精度較高的情況下，這一差距相對(duì)較大，可能影響航天器的碰撞風(fēng)險(xiǎn)分析，從而降低碰撞預(yù)警的精度。此外，對(duì)于近圓軌道物體，即使它們的遠(yuǎn)地點(diǎn)和近地點(diǎn)高度相差較小，但由于作為統(tǒng)計(jì)區(qū)域的圓環(huán)厚度和高度不同，空間物體在該區(qū)域內(nèi)的滯留時(shí)間(即有效統(tǒng)計(jì)時(shí)間)有可能遠(yuǎn)低于其軌道周期，從而導(dǎo)致該物體對(duì)航天器的碰撞概率遠(yuǎn)小于圓軌道假設(shè)下的結(jié)果。

文獻(xiàn)[18]從數(shù)據(jù)分析的角度，利用陳氏模型，對(duì)國際空間站在一定統(tǒng)計(jì)區(qū)域厚度內(nèi)來自空間物體的碰撞通量進(jìn)行了計(jì)算分析，并與MASTER模型和ORDEM模型的結(jié)果進(jìn)行對(duì)比。研究發(fā)現(xiàn)：隨著統(tǒng)計(jì)區(qū)域厚度與標(biāo)準(zhǔn)差比值的變化，ISS的碰撞通量先增后減，橢圓軌道假設(shè)下碰撞通量計(jì)算結(jié)果的變化范圍與ORDEM 3.0 模型以及MASTER模型的結(jié)果相接近；當(dāng)比值為2時(shí)，碰撞通量的值與ORDEM 3.0模型的平均碰撞通量值相近；利用陳氏模型計(jì)算碰撞通量的方法在時(shí)間上明顯優(yōu)于其它模型。具體內(nèi)容見參考文獻(xiàn)[18]。本文從理論出發(fā)，給出了空間物體橢圓軌道假設(shè)下航天器碰撞通量計(jì)算理論的詳細(xì)分析過程，以銥星86為例，計(jì)算了陳氏模型在圓軌道和橢圓軌道兩種假設(shè)條件下航天器的碰撞通量，并將計(jì)算結(jié)果進(jìn)行對(duì)比和分析。

1 圓軌道的碰撞通量

航天器的碰撞通量是指單位時(shí)間內(nèi)，航天器與多個(gè)空間物體發(fā)生碰撞的概率總和。當(dāng)航天器穿越多個(gè)空間物體時(shí)，航天器和空間物體各自有一個(gè)位置誤差橢球與它們的位置誤差相對(duì)應(yīng)，由于空間物體靠近航天器的方向未知，可以假設(shè)它們的聯(lián)合協(xié)方差各向同性，即誤差橢球?yàn)榍蛐?。令該球體的聯(lián)合協(xié)方差為σ2，則標(biāo)準(zhǔn)差為σ。

令航天器為主目標(biāo)，對(duì)應(yīng)的符號(hào)下標(biāo)用1表示，其它空間物體為從目標(biāo)，對(duì)應(yīng)的符號(hào)下標(biāo)用2表示。圖1[16]描述了緯度θ上主、從目標(biāo)的軌道，為便于說明，圖中只給出了主目標(biāo)作上升運(yùn)動(dòng)的情況；對(duì)應(yīng)的，從目標(biāo)可能作上升運(yùn)動(dòng)，也可能作下行運(yùn)動(dòng)。

為得到空間物體的數(shù)密度表達(dá)式，進(jìn)而得到統(tǒng)計(jì)區(qū)域內(nèi)空間物體對(duì)航天器的碰撞通量，將空間物體按軌道傾角進(jìn)行分類，首先考慮ntotal個(gè)同類別的空間物體在赤道平面上的分布情況，然后考慮其在各自軌道平面內(nèi)的統(tǒng)計(jì)分布，最后考慮一定厚度的高度空間，得到空間物體數(shù)密度關(guān)于統(tǒng)計(jì)區(qū)域厚度和軌道傾角的函數(shù)表達(dá)式。

兩個(gè)物體之間的碰撞概率是交會(huì)距離的函數(shù)，根據(jù)空間物體的數(shù)密度可以得到單位時(shí)間內(nèi)通過相遇平面的空間物體數(shù)量，乘以兩個(gè)物體之間的碰撞概率表達(dá)式，再對(duì)時(shí)間積分即可得到航天器穿越某空間區(qū)域時(shí)的碰撞通量。

航天器穿越空間物體云時(shí)，受到的碰撞通量如下所示[16]：

(1)

2 橢圓軌道的碰撞通量

2.1空間物體通過統(tǒng)計(jì)區(qū)域的情況分析

圖2給出了空間物體通過高度空間中某個(gè)統(tǒng)計(jì)區(qū)域時(shí)，其軌道相對(duì)該區(qū)域的位置關(guān)系。其中，兩個(gè)實(shí)線圓所圍起來的圓環(huán)區(qū)域?yàn)榻y(tǒng)計(jì)區(qū)域，圓環(huán)的內(nèi)、外半徑分別為r1和r2，橢圓軌道的近地點(diǎn)和遠(yuǎn)地點(diǎn)分別為d1和d2；虛線橢圓表示空間物體的運(yùn)行軌道；A、B是空間物體軌道和圓環(huán)的交點(diǎn)。

圓軌道假設(shè)下航天器的碰撞通量計(jì)算理論中，近地空間物體的軌道“幾乎為圓軌道”，它們的軌道完全處于厚度為h的統(tǒng)計(jì)區(qū)域內(nèi)，如圖2(a)所示。假設(shè)某個(gè)航天器的軌道半長(zhǎng)軸為R1，則該統(tǒng)計(jì)區(qū)域圓環(huán)的內(nèi)、外半徑分別為r1=(R1-h/2)，r2= (R1+h/2)。然而，空間物體的軌道并不一定完全被包含在所選定的某個(gè)統(tǒng)計(jì)區(qū)域中，根據(jù)空間物體的實(shí)際運(yùn)動(dòng)情況，可以將它們通過某個(gè)統(tǒng)計(jì)區(qū)域的狀態(tài)分成以下四種情況：

1)若某空間物體的軌道近地點(diǎn)離地心的距離大于等于該統(tǒng)計(jì)區(qū)域圓環(huán)的內(nèi)半徑，且遠(yuǎn)地點(diǎn)離地心的距離小于等于圓環(huán)的外半徑，則該空間物體的軌道完全包含在統(tǒng)計(jì)區(qū)域內(nèi)。

2)根據(jù)圖2(a)中對(duì)統(tǒng)計(jì)區(qū)域圓環(huán)內(nèi)外半徑的分析，當(dāng)通過該區(qū)域的空間物體的偏心率小于等于[h/(2R1)]時(shí)，該物體的軌道才有可能完全處于統(tǒng)計(jì)區(qū)域內(nèi)。然而，僅憑偏心率小這一條件并不能確定空間物體一定完全處在統(tǒng)計(jì)區(qū)域內(nèi)。有些空間物體的近地點(diǎn)處于區(qū)域內(nèi)，遠(yuǎn)地點(diǎn)在區(qū)域外，但仍然有很小的偏心率，如圖2(b)所示。根據(jù)軌道運(yùn)動(dòng)規(guī)律，在這種情況下，空間物體有較長(zhǎng)時(shí)間處在統(tǒng)計(jì)區(qū)域之外。

3)物體的近地點(diǎn)在區(qū)域外，遠(yuǎn)地點(diǎn)在區(qū)域內(nèi)，如圖2(c)所示。此種情況與情形2)類似，物體的近地點(diǎn)和遠(yuǎn)地點(diǎn)通常有且僅有一個(gè)處在統(tǒng)計(jì)區(qū)域內(nèi)，另一個(gè)在區(qū)域外。不同的是，物體由于遠(yuǎn)地點(diǎn)在區(qū)域內(nèi)部，使得其在統(tǒng)計(jì)區(qū)域內(nèi)的滯留時(shí)間相對(duì)較長(zhǎng)。但具有這種軌道的空間物體，其近地點(diǎn)軌道可能非常接近地球大氣層，受地球大氣阻力影響較大，導(dǎo)致其整個(gè)軌道運(yùn)行周期減短，同時(shí)導(dǎo)致遠(yuǎn)地點(diǎn)高度下降，只有小部分軌道處在統(tǒng)計(jì)區(qū)域內(nèi)，因而該空間物體滯留在統(tǒng)計(jì)區(qū)域內(nèi)的時(shí)間相對(duì)軌道周期仍然很短。

4)物體的近地點(diǎn)和遠(yuǎn)地點(diǎn)都不在區(qū)域內(nèi)，但軌道跨越了整個(gè)統(tǒng)計(jì)區(qū)域，如圖2(d)所示。此時(shí)空間物體的近地點(diǎn)距地心的距離小于r1，遠(yuǎn)地點(diǎn)距地點(diǎn)的距離大于r2。這種情況下，空間物體的有效統(tǒng)計(jì)時(shí)間是其通過統(tǒng)計(jì)區(qū)域時(shí)，中間兩部分軌道所花時(shí)間之和。

橢圓軌道假設(shè)下航天器的碰撞通量計(jì)算理論中，以空間物體在統(tǒng)計(jì)區(qū)域內(nèi)的滯留時(shí)間率作為其對(duì)碰撞通量的貢獻(xiàn)。所謂滯留時(shí)間率，是指對(duì)于某個(gè)空間物體，其代入計(jì)算的數(shù)量為該物體處在統(tǒng)計(jì)區(qū)域圓環(huán)內(nèi)的時(shí)間與其軌道周期的比值，即將某空間物體在統(tǒng)計(jì)區(qū)域內(nèi)滯留的時(shí)間記為Δt，軌道周期為T，則該物體的滯留時(shí)間率為Δt/T。圖2(a)中，物體完全處于區(qū)域內(nèi)，因而其數(shù)量統(tǒng)計(jì)值為1；圖2(b)～2(d)情況下的空間物體都只有部分軌道位于區(qū)域內(nèi)，因而這些物體的數(shù)量統(tǒng)計(jì)值應(yīng)小于1。

2.1.1 平近點(diǎn)角與軌道高度的關(guān)系表達(dá)式

橢圓軌道上某點(diǎn)對(duì)應(yīng)的時(shí)刻為物體的平均運(yùn)動(dòng)角速度與平近點(diǎn)角的函數(shù)：

M=n(t-τ)

(2)

式中:τ為物體過近地點(diǎn)的時(shí)刻。通過開普勒方程建立平近點(diǎn)角和高度之間的關(guān)系，根據(jù)高度得到平近點(diǎn)角，進(jìn)而可根據(jù)式(2)得到物體在對(duì)應(yīng)高度上的時(shí)刻。

空間物體的軌道方程為：

(3)

由真近點(diǎn)角f與偏近點(diǎn)角E的關(guān)系，結(jié)合開普勒方程

E-esinE=M

(4)

可得平近點(diǎn)角與軌道高度之間的關(guān)系，如下：

(5)

將式(5)代入式(2)得到橢圓軌道上的某點(diǎn)P與近地點(diǎn)之間的時(shí)間差如下：

(6)

根據(jù)式(3)可得對(duì)應(yīng)于高度r的真近點(diǎn)角f的值，但由于計(jì)算中涉及除以偏心率e的運(yùn)算，當(dāng)偏心率較小時(shí)，如小于0.01，如果直接利用式(3)求cosf進(jìn)而求真近點(diǎn)角f，容易造成結(jié)果的不穩(wěn)定，所以采用下面的方法進(jìn)行求解：

令

(7)

從而

(8)

由于0≤f≤π，式(8)不為負(fù)，根據(jù)泰勒展開：

當(dāng)0≤f≤π/2時(shí)，

(9)

當(dāng)π/2≤f≤π時(shí)，

(10)

將式(9)或式(10)代入式(6)計(jì)算軌道上P點(diǎn)與近地點(diǎn)之間的時(shí)間差。

當(dāng)偏心率較大時(shí)，可由式(3)直接計(jì)算真近點(diǎn)角f，代入式(6)得到軌道上P點(diǎn)與近地點(diǎn)之間的時(shí)間差。

2.1.2 空間物體統(tǒng)計(jì)數(shù)量的計(jì)算

兩行軌道根數(shù)(Two line element, TLE)是北美防空司令部(North America aerospace defend command, NORAD)發(fā)布的在軌編目目標(biāo)的軌道數(shù)據(jù)，它是目前最完備的地球軌道空間物體編目數(shù)據(jù)。根據(jù)TLE參數(shù)，容易得到航天器的軌道半長(zhǎng)軸，若給定統(tǒng)計(jì)區(qū)域圓環(huán)的厚度，則可以得到圓環(huán)的內(nèi)、外半徑r1和r2。分析圖2中的四種情況，可知空間物體通過該區(qū)域時(shí)，有效統(tǒng)計(jì)時(shí)間對(duì)應(yīng)的第一個(gè)高度為圓環(huán)的內(nèi)半徑或空間物體近地點(diǎn)離地心的距離，第二個(gè)高度為圓環(huán)的外半徑或物體遠(yuǎn)地點(diǎn)離地心的距離。針對(duì)不同的情況，分別將相應(yīng)的高度值代入式(3)、式(9)或式(10)，再根據(jù)式(6)計(jì)算得到空間物體與圓環(huán)的交點(diǎn)和近地點(diǎn)之間的時(shí)間差，以高度r的函數(shù)g(r)表示。

1)圖2(a)中，空間物體軌道完全處于統(tǒng)計(jì)區(qū)域內(nèi)，滯留時(shí)間為物體的軌道周期，即：Δt=T。

2)圖2(b)中，物體的近地點(diǎn)位于統(tǒng)計(jì)區(qū)域內(nèi)。其第一個(gè)高度為物體的近地點(diǎn)離地心的距離；第二個(gè)高度，即交點(diǎn)A2對(duì)應(yīng)的高度值，等于r2。令r等于r2，根據(jù)式(2)～(10)可得A2與近地點(diǎn)之間的時(shí)間差g(r2)，空間物體在統(tǒng)計(jì)區(qū)域內(nèi)的滯留時(shí)間為該時(shí)間差的2倍，即：Δt=2g(r2)。

3)圖2(c)中，物體的遠(yuǎn)地點(diǎn)位于統(tǒng)計(jì)區(qū)域內(nèi)。其第一個(gè)高度應(yīng)為交點(diǎn)A3對(duì)應(yīng)的高度值，等于r1；第二個(gè)高度為空間物體的遠(yuǎn)地點(diǎn)離地心的距離。令r為r1，得A3與近地點(diǎn)之間的時(shí)間差g(r1)，此時(shí)，空間物體的滯留時(shí)間為半軌道周期與該時(shí)間差的差值的2倍：Δt=2(T/2-g(r1))。

4)圖2(d)中，兩個(gè)交點(diǎn)高度對(duì)應(yīng)的值分別為圓環(huán)的內(nèi)、外半徑，分別將r1和r2代入計(jì)算，根據(jù)式(2)～(10)可得交點(diǎn)A4和B4與近地點(diǎn)之間的時(shí)間差，分別為g(r1)和g(r2)，物體的滯留時(shí)間為這兩者之間差值的2倍：Δt=2(g(r2)-g(r1))。

對(duì)于穿越統(tǒng)計(jì)區(qū)域圓環(huán)的空間物體，根據(jù)上述不同的軌道運(yùn)行情況，可以得到該空間物體相應(yīng)的滯留時(shí)間。計(jì)算航天器在統(tǒng)計(jì)區(qū)域內(nèi)的總碰撞通量時(shí)，本文采用依次計(jì)算其與單個(gè)物體的碰撞通量再累加的辦法。記物體j的滯留時(shí)間與其軌道周期的比值為nj，將nj代替式(1)中的ntotal計(jì)算K值，進(jìn)而得到該空間物體對(duì)航天器的碰撞通量。比值nj表示為：

nj=Δtj/Tj

(11)

2.2航天器基于橢圓軌道假設(shè)的碰撞通量

結(jié)合式(1)和式(11)，可以得到航天器實(shí)際運(yùn)行軌道中與空間物體發(fā)生碰撞的次數(shù)。航天器運(yùn)動(dòng)一個(gè)周期的總碰撞通量計(jì)算式如下所示：

(12)

式中：Kj、Ij分別表示物體j所代表的量，其余變量的含義同式(1)。

將航天器運(yùn)動(dòng)一周的碰撞通量乘以航天器一年的運(yùn)動(dòng)圈數(shù)，得到該航天器一年的總碰撞通量：

(13)

式中：T1為航天器運(yùn)行一圈所花的時(shí)間，單位為s，Y表示1年，計(jì)算時(shí)需將其轉(zhuǎn)化為s。

2.2.1 航天器碰撞通量計(jì)算

針對(duì)圓軌道假設(shè)和橢圓軌道假設(shè)兩種情況，空間物體計(jì)入統(tǒng)計(jì)的判斷條件有所不同。

圓軌道假設(shè)對(duì)應(yīng)的理論認(rèn)為，近地空間中，空間物體的軌道基本上為圓軌道，它們對(duì)碰撞概率的影響是主要的，大橢圓軌道對(duì)航天器的碰撞概率可以忽略不計(jì)。在這種情況下，只需要考慮軌道完全包括在統(tǒng)計(jì)區(qū)域圓環(huán)內(nèi)的空間物體。通過選取航天器平均軌道高度上下各一定空間范圍內(nèi)確定的統(tǒng)計(jì)區(qū)域，只要空間物體的遠(yuǎn)、近地點(diǎn)高度均位于統(tǒng)計(jì)區(qū)域內(nèi)，就令nj為1，再代入式(12)和式(13)計(jì)算該空間物體對(duì)航天器的年碰撞通量。

然而，在該假設(shè)中，要求統(tǒng)計(jì)區(qū)域足夠大，才可能得到可靠的計(jì)算結(jié)果。這是因?yàn)椋y(tǒng)計(jì)區(qū)域取得過小，符合條件的空間物體數(shù)量就很少，與實(shí)際空間環(huán)境相差較大，所得碰撞通量與實(shí)際值相差較大。并且統(tǒng)計(jì)區(qū)域取得過小，區(qū)域厚度微小的變化就可能導(dǎo)致碰撞通量值的較大波動(dòng)，這將導(dǎo)致計(jì)算結(jié)果的不穩(wěn)定。

橢圓軌道假設(shè)中，考慮了空間物體實(shí)際處在統(tǒng)計(jì)區(qū)域內(nèi)的有效軌道，包括了圖2中的四種情況。即空間物體的近、遠(yuǎn)地點(diǎn)均位于統(tǒng)計(jì)區(qū)域內(nèi)(圖2(a))；空間物體的近、遠(yuǎn)地點(diǎn)有一個(gè)處在統(tǒng)計(jì)區(qū)域內(nèi)，另一個(gè)在區(qū)域外(圖2(b)～2(c))；空間物體的近地點(diǎn)和遠(yuǎn)地點(diǎn)都不在統(tǒng)計(jì)區(qū)域內(nèi)，但中間的一部分軌道處在區(qū)域中(圖2(d))。將部分或全部軌道都處在統(tǒng)計(jì)區(qū)域內(nèi)的空間物體計(jì)入統(tǒng)計(jì)，計(jì)算它們?cè)诮y(tǒng)計(jì)區(qū)域內(nèi)的滯留時(shí)間，再利用式(12)和式(13)計(jì)算這些物體對(duì)航天器的碰撞通量。

2.2.2 參數(shù)取值分析

取航天器的高度為其軌道高度的平均值；取航天器高度的上、下各150 km作為統(tǒng)計(jì)區(qū)域，即厚度h=300 km；空間物體的大小不一，小的為幾毫米，大的可達(dá)幾十米， TLE編目目標(biāo)中，最小空間物體的尺寸約為10 cm，考慮到物體的尺寸越小，數(shù)量越多，本文將空間物體的半徑取為20 cm，將該值加上航天器的半徑值即得到碰撞半徑rA；積分上限取兩物體的傾角中較小正弦值對(duì)應(yīng)的軌道傾角或其補(bǔ)角。那么，若令航天器和空間物體的軌道傾角分別為i1，i2，令sini=min(sini1,sini2)，則θ″的取值為：

(14)

此外，為檢驗(yàn)碰撞通量對(duì)統(tǒng)計(jì)區(qū)域厚度的靈敏度，本文還計(jì)算了厚度h= 5, 10, 20, 30, 40, 50, 100, 200 km時(shí)航天器的總碰撞通量。

2.2.3 航天器碰撞通量算法實(shí)現(xiàn)

本文使用TLE編目數(shù)據(jù)計(jì)算航天器的碰撞通量。提取編目目標(biāo)的主要參數(shù)如軌道編號(hào)、軌道傾角、軌道半徑等，然后將編目目標(biāo)按軌道傾角進(jìn)行分組。分組間距為5°，通過分析，所有編目的空間物體的軌道傾角最大值為150°，因而將物體按(0,5]，(5,10]，…，(145,150]分成30組。對(duì)每一組中的物體逐一計(jì)算，并分類顯示其計(jì)算結(jié)果。

對(duì)TLE參數(shù)預(yù)處理后，提取統(tǒng)計(jì)區(qū)域內(nèi)符合條件的目標(biāo)，將其參數(shù)代入相應(yīng)的計(jì)算式，得到其對(duì)航天器的碰撞通量。根據(jù)圓軌道和橢圓軌道兩種假設(shè)對(duì)應(yīng)的理論，圓軌道假設(shè)下，當(dāng)空間物體的遠(yuǎn)、近地點(diǎn)軌道高度均在統(tǒng)計(jì)區(qū)域內(nèi)時(shí)，將該物體的參數(shù)代入計(jì)算；橢圓軌道假設(shè)下，只要空間物體有軌道在統(tǒng)計(jì)區(qū)域內(nèi)，就將該物體的參數(shù)代入計(jì)算，也就是說，除了考慮物體的近地點(diǎn)或遠(yuǎn)地點(diǎn)在統(tǒng)計(jì)區(qū)域內(nèi)的情形，還應(yīng)該考慮空間物體的軌道橫跨統(tǒng)計(jì)區(qū)域的情況。此外，圓軌道假設(shè)下，對(duì)計(jì)入統(tǒng)計(jì)的空間物體，令式(11)中的nj為1，再代入式(12)計(jì)算其對(duì)航天器的碰撞通量；橢圓軌道假設(shè)下，對(duì)滿足條件的空間物體，則應(yīng)該判斷物體軌道相對(duì)統(tǒng)計(jì)區(qū)域圓環(huán)的位置關(guān)系，進(jìn)而根據(jù)不同的情況計(jì)算物體在區(qū)域內(nèi)的滯留時(shí)間，再代入式(12)計(jì)算航天器的碰撞通量。圖3給出了圓軌道和橢圓軌道假設(shè)下的程序計(jì)算流程圖，圖中對(duì)兩種假設(shè)條件下空間物體的判斷依據(jù)進(jìn)行了區(qū)分。

3 計(jì)算結(jié)果對(duì)比與分析

以銥星86(NOARD ID: 25528)為例，計(jì)算了其在圓軌道和橢圓軌道兩種假設(shè)條件下的碰撞通量。取銥星的半徑為5 m，根據(jù)TLE數(shù)據(jù)，得到其主要參數(shù)如表1所示。根據(jù)式(12)得到銥星86運(yùn)行一圈的碰撞通量，再由式(13)計(jì)算得到一年的總碰撞通量。

表2給出了厚度h=5，10，20，30，40，50，100，200，300 km時(shí)，銥星86的年總碰撞通量以及計(jì)入統(tǒng)計(jì)的空間物體數(shù)量。

對(duì)比兩種情況下銥星86的數(shù)據(jù)統(tǒng)計(jì)結(jié)果發(fā)現(xiàn)，圓軌道假設(shè)下空間物體的統(tǒng)計(jì)數(shù)量隨著統(tǒng)計(jì)區(qū)域厚度的增加波動(dòng)幅度較大，由統(tǒng)計(jì)區(qū)域厚度為5 km時(shí)的72個(gè)增加到厚度為300 km時(shí)的4000多個(gè)。將圓軌道假設(shè)與橢圓軌道假設(shè)的統(tǒng)計(jì)值相比，后者對(duì)前者的比值由60倍降低為2倍。就空間物體對(duì)銥星86的碰撞通量來看，橢圓軌道假設(shè)與圓軌道假設(shè)的計(jì)算值之比由大約4倍降為1.4倍。

圖4給出了銥星86的統(tǒng)計(jì)區(qū)域內(nèi)圓軌道和橢圓軌道兩種假設(shè)情況下計(jì)入統(tǒng)計(jì)的空間物體軌道偏心率分布情況。限于篇幅，這里只給出h= 10，50，100 km的偏心率分布圖。其中，左列為圓軌道假設(shè)下空間物體的偏心率分布；右列為橢圓軌道假設(shè)下空間物體的偏心率分布。所有計(jì)入統(tǒng)計(jì)的空間物體按偏心率(0,0.01]，(0.01,0.02]，(0.02,0.03]，(0.03,0.04]，(0.04,0.05]，(0.05,0.1]，(0.1,0.2]，…，(0.8,0.9]進(jìn)行分組，統(tǒng)計(jì)每一組空間物體的數(shù)量。

由圖4可知，對(duì)于計(jì)入統(tǒng)計(jì)的空間物體，其軌道大部分屬于近圓軌道。但由于圓軌道假設(shè)將統(tǒng)計(jì)目標(biāo)的軌道嚴(yán)格限定在統(tǒng)計(jì)區(qū)域內(nèi)，使得當(dāng)統(tǒng)計(jì)區(qū)域厚度較小時(shí)，很大一部分沒有完全處于統(tǒng)計(jì)區(qū)域的近圓軌道空間物體被排除，從而造成了較大的計(jì)算誤差。

表1 銥星86的主要軌道參數(shù)Table 1 Orbital elements of the Iridium satellites 86

假設(shè)某空間物體的遠(yuǎn)、近地點(diǎn)軌道高度正好等于統(tǒng)計(jì)區(qū)域圓環(huán)的外半徑和內(nèi)半徑，則其對(duì)應(yīng)的偏心率最大。取銥星86統(tǒng)計(jì)區(qū)域厚度h=300 km，此時(shí)，內(nèi)圓環(huán)半徑為7005.8 km，外圓環(huán)半徑為7305.8 km，圓軌道假設(shè)中符合條件的空間物體的最大偏心率約為0.021。

表2 不同區(qū)域厚度下銥星的碰撞通量Table 2 Number of collisions of Iridium satellite in different layer thickness

根據(jù)圖4，當(dāng)統(tǒng)計(jì)區(qū)域厚度較大時(shí)，兩種假設(shè)條件下的碰撞通量計(jì)算結(jié)果比值較小，但仍然有很大一部分空間物體被圓軌道假設(shè)條件排除。例如，當(dāng)統(tǒng)計(jì)區(qū)域厚度為300 km時(shí)，偏心率大于0.021的空間物體將被排除；小部分偏心率小于0.021的空間物體因?yàn)闈M足圖2(b)～2(d)的情形被圓軌道假設(shè)條件剔除。這造成了統(tǒng)計(jì)結(jié)果的不準(zhǔn)確，進(jìn)而影響到碰撞通量的計(jì)算結(jié)果。

橢圓軌道假設(shè)中，部分計(jì)入統(tǒng)計(jì)的空間物體偏心率較大，但因?yàn)閷魰r(shí)間率作為空間物體的統(tǒng)計(jì)數(shù)量用于計(jì)算碰撞通量，所考慮的情形與實(shí)際空間環(huán)境相符，可以避免由于遺漏空間物體造成計(jì)算結(jié)果的不準(zhǔn)確。

4 結(jié) 論

考慮空間物體橢圓運(yùn)行軌道實(shí)際，本文根據(jù)開普勒方程和空間物體運(yùn)動(dòng)方程，采用陳氏模型，提出了空間物體橢圓軌道假設(shè)下航天器碰撞通量的計(jì)算方法。通過比較銥星86在兩種軌道假設(shè)下的計(jì)算結(jié)果發(fā)現(xiàn)，隨著統(tǒng)計(jì)區(qū)域厚度的增加，橢圓軌道與圓軌道的碰撞通量之比從大約4倍降為1.4倍。

由于在空間尺度上進(jìn)行計(jì)算，小偏心率的空間物體也可能具有較大的遠(yuǎn)、近地點(diǎn)高度差。這些空間物體由于穿越統(tǒng)計(jì)區(qū)域而對(duì)航天器的碰撞通量有貢獻(xiàn)。圓軌道假設(shè)只將軌道完全處于統(tǒng)計(jì)區(qū)域內(nèi)的空間物體計(jì)入統(tǒng)計(jì)，忽略了圖2(b)～2(d)的情況，與實(shí)際的空間環(huán)境相差較大；橢圓軌道考慮了空間物體的實(shí)際運(yùn)動(dòng)情況，可以避免遺漏統(tǒng)計(jì)目標(biāo)，提高計(jì)算結(jié)果的準(zhǔn)確度。

[1] 柳森, 蘭勝威, 李毅, 等. 航天器解體模型研究綜述[J]. 宇航學(xué)報(bào), 2010, 31(1):14-23.[Liu Sen, Lan Sheng-wei, Li Yi, et al. Review of spacecraft breakup model[J]. Journal of Astronautics, 2010, 31(1): 14-23.]

[2] 李怡勇, 沈懷榮, 李智, 等. 航天器撞擊解體碎片的短期危害評(píng)估[J]. 宇航學(xué)報(bào), 2010, 31(4):1231-1236.[Li Yi-yong, Shen Huai-rong, Li Zhi, et al. Assessing the short-term hazard of spacecraft collision breakup debris[J]. Journal of Astronautics, 2010, 31(4): 1231-1236.]

[3] NASA. Orbital debris quarterly news[EB/OL]. 2016[2017]. https://orbitaldebris.jsc.nasa.gov/quarterly- news/pdfs/odqnv20i4.pdf.

[4] 白顯宗, 陳磊, 張翼, 等. 空間物體碰撞預(yù)警技術(shù)研究綜述[J]. 宇航學(xué)報(bào), 2013, 34(8):1027-1039.[Bai Xian-zong, Chen Lei, Zhang Yi, et al. Survey on collision assessment and warning techniques for space objects[J]. Journal of Astronautics, 2013, 34(8): 1027-1039.]

[5] Foster J L, Wortham M B. ISS debris avoidance maneuver threshold analysis[R]. Washington, USA: National Aeronautics and Space Administration, May 2007.

[6] Carpenter J R, Markley F L. Wald sequential probability ratio test for space object conjunction assessment[J]. Journal of Guidance Control and Dynamics, 2014, 37(5): 1385-1396.

[7] Zhou W P, Guo X Z, Shen M, et al. Discussion on collision probability threshold[J]. Journal of Beijing Institute of Technology, 2015, 24: 61-65.

[8] Kessler D J. Derivation of the collision probability between orbiting objects: the lifetimes of Jupiter’s outer moons[J]. Icarus, 1981, 48(1): 39-48.

[9] Wang D F, Pang B J, Xiao W K. GEO space debris environment determination in the Earth fixed coordinate system[C]. 7th European Conference on Space Debris, Darmstadt, Germany, April 18-21, 2017.

[10] Foster J L, Frisbee J H. Comparison of the exclusion volume and probability threshold methods for debris avoidance for the STS orbiter and international space station[R]. Washington, USA: National Aeronautics and Space Administration, May 2007.

[11] Liou J C, Matney M J, Anz-meador P, et al. The new NASA orbital debris engineering model ORDEM2000[R]. Houston, USA: Johnson Space Center, May 2002.

[12] Matney M. NASA’s orbital debris environmental model ORDEM 3.0-implications for measurements and modeling[C]. 33rd Inter-Agency Space Debris Meeting, Houston, USA, March 30-April 2, 2015.

[13] Klinkrad H. Space debris models and risk analysis[M]. Chichester: Springer-Verlag, 2006.

[14] 劉舒蒔, 龔建村, 劉四清, 等. 中長(zhǎng)期軌道預(yù)報(bào)中大氣阻力系數(shù)補(bǔ)償算法的研究[J]. 宇航學(xué)報(bào), 2013, 34(2):157-162.[Liu Shu-shi, Gong Jian-cun, Liu Si-qing, et al. Atmospheric drag coefficient calibration in medium-term orbit prediction[J]. Journal of Astronautics, 2013, 34(2): 157-162.]

[15] 劉衛(wèi), 王榮蘭, 王四清, 等. 基于小波變換的衛(wèi)星阻力系數(shù)分析[J]. 宇航學(xué)報(bào), 2015, 36(2):142-150.[Liu Wei, Wang Rong-lan, Wang Si-qing, et al. Analysis of satellite drag coefficient based on wavelet transformation[J]. Journal of Astronautics, 2015, 36(2): 142-150.]

[16] Chan F K. Spacecraft collision probability[M]. California: The Aerospace Press, 2008.

[17] Dennis N G. Probabilistic theory and statistical distribution of Earth satellites[J]. Journal of the British Interplanetary Society, 1972, 25: 333-376.

[18] Zhou W P, Chan F K. Comparison of debris collision fluxes for the international space station[J]. Journal of Beijing Institute of Technology, 2016, 25: 41-49.

[19] 陳俊宇, 李彬, 章品,等. 低軌道空間碎片彈道系數(shù)及應(yīng)用[J]. 紅外與激光工程, 2016, 45(11).[Chen Jun-yu, Li Bin, Zhang Pin, et al. Low Earth orbit space debris ballistic coefficients and their applications[J]. Infrared and Laser Engineering, 2016, 45(11).]

CalculationsandAnalysisonImpactFluxesforSpaceObjectsunderAssumptionofEllipticalOrbits

ZHOU Wei-ping1,2, SHEN Ming1, GAO Peng-qi1, GUO Xiao-zhong1, ZHAO You1

(1. National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China; 2. University of Chinese Academy of Sciences, Beijing 100049, China)

This paper studies the impact fluxes calculation for the orbiting space objects via the Chan’s model which does not invoke either the kinetic theory of gases model or the Poisson distribution function. Four cases are discussed for the space objects on elliptical orbits passing through a thin spherical shell surrounding the spacecraft. The paper also compares the circular orbit assumption with the elliptical assumption. The satellite Iridium 86 is used to estimate the impact fluxes under both assumptions. It shows that the ratio of the impact flux between the elliptical orbits case and the circular case decreases from about 4 times to 1.4 times as the thickness increases. The statistical number of the objects gets a big fluctuation as the increased thickness in circular case. The elliptical case is 60 times more than the circular case when the shell’s thickness equals to 5 km.

Collision warning; Impact flux; Space objects; Spacecraft; Satellite orbits

V528

A

1000-1328(2017)11- 1234- 09

10.3873/j.issn.1000- 1328.2017.11.013

2017- 04- 18;

2017- 09- 05

國家重大科研裝備研制資助項(xiàng)目(ZDYZ2013-2)

周威萍(1990-)，女，博士生，主要從事空間物體碰撞概率閾值研究工作。

通信地址：北京市朝陽區(qū)大屯路甲20號(hào)國家天文臺(tái)(100012)

電話：18210938053

E-mail：wpzhou@nao.cas.cn

趙有(1964-)，男，博士，研究員，主要從事空間物體與碎片監(jiān)測(cè)預(yù)警的研究工作。本文通信作者。

通信地址：北京市朝陽區(qū)大屯路甲20號(hào)國家天文臺(tái)(100012)

電話：64807617

E-mail：youzhao@nao.cas.cn