空間物體解體碎片云的長期演化建模與分析

張斌斌，王兆魁，張育林，2

1.國防科技大學(xué)　航天科學(xué)與工程學(xué)院，長沙 410073 2.清華大學(xué)　航天航空學(xué)院，北京 100084

?

空間物體解體碎片云的長期演化建模與分析

張斌斌1，王兆魁2，*，張育林1，2

1.國防科技大學(xué)航天科學(xué)與工程學(xué)院，長沙 410073 2.清華大學(xué)航天航空學(xué)院，北京 100084

空間碎片云由空間物體解體產(chǎn)生的大量空間碎片組成，由于其相對集中地分布在有限的空間內(nèi)，將會對臨近航天器產(chǎn)生較大的碰撞威脅。為了分析解體碎片云長期分布特點，文章首先利用數(shù)值積分方法對空間碎片云短期分布規(guī)律進(jìn)行了研究；在此基礎(chǔ)上，針對處于環(huán)狀分布的碎片云，根據(jù)碎片所在的軌道高度和具有的面質(zhì)比值，將碎片劃分到不同分組，以每個組作為研究對象，建立了描述碎片云在大氣阻力作用下的解析演化模型。模型避免了對單個解體碎片的運動狀態(tài)進(jìn)行積分，可大大降低對計算資源和計算時間的需求?？紤]在高度為1 422 km圓軌道上運行的物體，解體產(chǎn)生了1 780個碎片，利用解析演化模型得到碎片云未來50年內(nèi)的演化分布狀態(tài)。數(shù)值結(jié)果表明，碎片云的峰值密度在解體物體軌道高度附近，并在大氣阻力作用下向更大高度區(qū)間內(nèi)擴散；較低高度區(qū)間內(nèi)碎片密度具有先增加，然后在大氣阻力作用下不斷減少的特點。

空間碎片云；空間密度；解體模型；碎片分布；演化模型

空間碎片云由大量解體碎片構(gòu)成，這些解體碎片由火箭上面級、在軌運行航天器等人造物體解體產(chǎn)生。隨著人類航天活動的不斷增加，由空間物體相互碰撞而解體產(chǎn)生的空間碎片數(shù)量日益增多，截止到2016年6月，已經(jīng)記錄的空間物體超高速災(zāi)難性碰撞解體事件有4次，分別是1991年COSMOS 1934航天器和編目的13475號空間碎片碰撞，1996年法國“櫻桃”號航天器和編目的18208號空間碎片碰撞，2005年THOR BURNER 2A火箭體與編目的26207號空間碎片碰撞，以及2009年COSMOS 2251航天器IRIDIUM 33航天器發(fā)生碰撞[1]。每次航天器的超高速劇烈碰撞事件，都會引入大量的解體碎片，如2009年COSMOS 2251和IRIDIUM 33碰撞解體產(chǎn)生2 201個編目碎片，占在軌總編目碎片數(shù)量的1/8。可以看出，空間物體之間相互碰撞形成的解體碎片，會對在軌工作航天器產(chǎn)生巨大的碰撞威脅，將是未來空間碎片增長主要原因[2-4]。對空間物體解體產(chǎn)生的碎片云的長期分布演化狀態(tài)進(jìn)行研究，對于分析預(yù)測空間碎片環(huán)境的未來演化趨勢、分析航天器受碎片云的碰撞威脅等具有重要的意義。

空間碎片云的分布演化不僅與解體物體的質(zhì)量、尺寸以及運動狀態(tài)相關(guān)，還受大氣阻力、非球形引力等復(fù)雜攝動力作用的影響。通常分別建立空間物體解體模型和長期演化計算模型，來分別描述空間物體解體產(chǎn)生碎片情況和碎片云的長期分布演化情況。文獻(xiàn)中一般利用美國國家航空航天局(NASA)標(biāo)準(zhǔn)解體模型來模擬產(chǎn)生物體解體碎片[5]，但該模型是有局限的，未能考慮解體過程中解體物體與碎片的質(zhì)量和動量守恒。針對空間碎片分布的長期演化分析問題，文獻(xiàn)[6]建立了空間碎片碰撞風(fēng)險分析工具(SDIRAT)；文獻(xiàn)[7]設(shè)計了空間碎片仿真系統(tǒng)(SDS)，利用該仿真系統(tǒng)得到，若星座衛(wèi)星附近發(fā)生碰撞解體事件，解體碎片云將會使星座衛(wèi)星受到的碰撞威脅急劇增加，但該仿真系統(tǒng)僅適合在幾個小時或數(shù)月內(nèi)對解體碎片云的分布狀態(tài)進(jìn)行仿真分析。上述空間碎片分布演化模型，以及文獻(xiàn)[8]中建立的LEGEND模型、文獻(xiàn)[9]中建立的SDM模型，均以單個解體碎片的運動狀態(tài)為演化計算對象，當(dāng)進(jìn)行長期分析時，不僅需要較長計算時間，消耗大量計算資源，而且不能保證碎片云分布狀態(tài)長期演化計算結(jié)果的精度[10]。

本文首先在NASA標(biāo)準(zhǔn)模型的基礎(chǔ)上，通過數(shù)值模擬產(chǎn)生解體碎片；然后通過對解體碎片的運動狀態(tài)進(jìn)行數(shù)值積分，演化計算得到碎片云短期分布狀態(tài)，進(jìn)而分析得到碎片云的分布特點；最后，基于解體碎片云的分布特點，建立了大氣阻力作用下空間碎片云的解析演化模型，并利用該解體模型分析了碎片云的長期分布演化狀態(tài)。

1　解體模型和碎片云短期分布特點的數(shù)值分析

1.1解體模型及其數(shù)值仿真實現(xiàn)

(1)

(2)

利用式(2)可以確定每個解體碎片的尺寸。進(jìn)一步可以通過下述雙正態(tài)分布方程確定解體碎片的面質(zhì)比：

(3)

在解體時刻，解體碎片獲得相對于解體物體的速度增量Δv滿足正態(tài)分布函數(shù)，即

(4)

綜合空間物體解體產(chǎn)生碎片過程中質(zhì)量守恒和動量守恒約束，利用解體模型模擬生成解體碎片的流程如圖1所示。計算時，首先確定要研究解體碎片的最小尺寸lc，然后根據(jù)關(guān)系式(1)和計算解體產(chǎn)生大于尺寸碎片的總數(shù)目，再根據(jù)分布律式(2)、(3)和(4)確定解體碎片具有的尺寸、面質(zhì)比和速度增量。應(yīng)用分布律式(2)、(3)和(4)時，需要對隨機變量抽樣，可采用反函數(shù)的方法進(jìn)行隨機變量的抽樣[11]。

以一顆Globalstar星座衛(wèi)星受空間碎片碰撞解體為例，說明上述解體模型數(shù)值實現(xiàn)。Globalstar星座衛(wèi)星質(zhì)量約為800 kg，運行在高度為1 422 km的圓軌道上，解體時刻在J2000坐標(biāo)系中的位置和速度矢量分別為r=[7 784.4 00 0 -0.001 305] km，v=[0 4.311 5.721] km/s。星座衛(wèi)星受空間碎片碰撞后，完全解體，產(chǎn)生1 780個尺寸大于5 cm的碎片。利用解體模型模擬產(chǎn)生的碎片分布情況如圖2和圖3所示。

圖1　NASA標(biāo)準(zhǔn)解體模型的實現(xiàn)流程Fig.1　Computation flow diagram of the breakup model

圖2　解體碎片面質(zhì)比分布的對比分析Fig.2　Breakup fragments′ area-mass ratio distribution

圖3　解體碎片速度增量分布的對比分析Fig.3　Breakup fragments′ velocity increment distribution

圖2(a)和圖3(a)分別是利用解體模型模擬得到的解體碎片面質(zhì)比分布和速度增量分布，圖2(b)和圖3(b)則是利用1 780個在軌解體碎片的長期觀測數(shù)據(jù)，分析得到的解體碎片面質(zhì)比分布和速度增量分布情況[5]。對比圖2(a)和圖2(b)，圖3(a)和圖3(b)可以看出，利用解體模型得到的解體碎片的面質(zhì)比、速度增量分布與實際觀測數(shù)據(jù)一致，說明了解體模型和實現(xiàn)方法的正確性。

1.2碎片云分布特點分析

為了分析碎片云的分布特點，綜合考慮大氣阻力、地球非球形引力、太陽光壓以及日、月三體引力等攝動力的作用，利用數(shù)值積分的方法對解體碎片的運動狀態(tài)進(jìn)行計算更新。數(shù)值計算中，大氣阻力計算采用Harris-Priester大氣模型，該模型能夠反映太陽輻射對大氣密度的長期影響，且計算復(fù)雜度較低；地球非球形引力中考慮了J2、J3和J4等帶諧項，以及J2，2田諧項的影響；數(shù)值積分采用4階Adams-Bashforth/Adams-Moulton預(yù)測修正方法，利用4階Runge-Kutta方法進(jìn)行積分的初始化[12]。

對第1.1節(jié)模擬生成的1780個解體碎片的運動狀態(tài)進(jìn)行數(shù)值積分計算，計算結(jié)果如圖4～圖6所示。

圖4(a)、圖5(a)、圖6(a)給出了解體碎片云的空間密度，即單位空間體積內(nèi)解體碎片的數(shù)量，在不同時刻、不同高度上的分布。圖4(b)、圖5(b)、圖6(b)給出了解體碎片云在3維空間中的分布情況。對比圖4～圖6可以看出，解體碎片云的分布狀態(tài)隨時間演化大致可以分為3個階段：集中分布階段、帶狀分布階段及環(huán)狀分布階段。空間物體解體后很短一段時間內(nèi)，解體碎片處于集中分布階段，如圖4(a)所示，此階段內(nèi)解體碎片呈現(xiàn)密集團狀分布，該分布特點主要由解體時刻碎片獲取的速度增量決定。如圖5(b)所示，由于解體碎片獲取的速度增量不同，導(dǎo)致解體碎片的軌道周期和相位出現(xiàn)差異，經(jīng)過1個月左右解體碎片云在近似一個軌道面上形成均勻的圓帶狀，進(jìn)入了帶狀分布階段。對比圖4(a)和圖5(a)可以看出，相比于集中分布階段，帶狀階段內(nèi)的解體碎片不僅相位上呈現(xiàn)差異，也在逐漸向不同高度上擴散。在帶狀分布的基礎(chǔ)上，受主要攝動力地球非球形J2項引力的作用：

圖4　物體解體時刻碎片云的空間分布Fig.4　Space cloud′s spatial density distribution when the breakup event accured

圖5　物體解體后1個月尺寸大于5cm解體碎片的空間分布Fig.5　Space cloud′s spatial density distribution after the breakup event

(5)

式中：RE為地球半徑；μ為地球引力常數(shù)；a、e和i分別為碎片軌道半長軸、偏心率和傾角。

不同解體碎片軌道的升交點赤經(jīng)和近地點幅角的變化率呈現(xiàn)差異，并在一段時間后趨于均勻分布，形成環(huán)狀分布，如圖6(b)所示。碎片云進(jìn)入環(huán)狀分布階段后，將主要在大氣阻力作用下，沿軌道高度方向變化。

圖6　物體解體后5個月尺寸大于5cm解體碎片的空間分布Fig.6　Space cloud′s spatial density distribution after the breakup event

從上述分析過程可以看出，利用數(shù)值積分方法，可以得到解體碎片在未來一段時間內(nèi)的運動狀態(tài)，進(jìn)而分析解體碎片云的分布特點。但數(shù)值積分方法，以每個解體碎片為計算對象，對計算資源要求高，演化計算時間也較長。如在上述分析中，利用配置為32核(單核主頻為2.13 GHz)的計算服務(wù)器，需要15 h的計算時間才能得到如圖6(b)所示的碎片分布狀態(tài)。

2　大氣阻力作用下空間碎片云的解析演化模型

考慮到利用數(shù)值積分方法對解體碎片的運動狀態(tài)進(jìn)行長期演化分析，將需要大量計算資源和時間。本節(jié)將在一定假設(shè)的基礎(chǔ)上，基于空間碎片云環(huán)狀分布特點，將解體碎片云按高度進(jìn)行分層離散化，以每個高度層內(nèi)解體碎片作為一個整體進(jìn)行研究，建立碎片云的解析演化模型，實現(xiàn)對碎片云分布狀態(tài)的長期分析。

當(dāng)解體碎片云進(jìn)入如圖6所示環(huán)狀分布狀態(tài)后，在地球非球形攝動力作用下，解體碎片軌道的緯度幅角和升交點赤經(jīng)不斷均勻變化，同一軌道高度上的解體碎片的分布將是均勻和相對穩(wěn)定的。此后，解體碎片云將會在大氣阻力的作用下，沿軌道高度方向變化。將解體碎片云按軌道高度劃分到Nh個高度區(qū)間內(nèi)。同時考慮到不同面質(zhì)比物體，受到大氣阻力攝動力不同，進(jìn)一步將同一高度區(qū)間內(nèi)的解體碎片劃分到Na個面質(zhì)比區(qū)間內(nèi)。同一高度和面質(zhì)比區(qū)間內(nèi)的解體碎片可視為一個獨立分組，以每個分組為研究對象，建立解體碎片云演化的微分控制方程，實現(xiàn)對碎片云分布的長期分布演化。

取第i高度區(qū)間為控制體，該控制體內(nèi)物體軌道高度差為Δr，如圖7所示。

圖7　碎片云分布空間分層與控制體選取示意Fig.7　Grouping and control volume setting of space debris cloud

假設(shè)解體碎片均運行在圓軌道上，則在大氣阻力攝動作用下，控制體內(nèi)位于第j個面質(zhì)比區(qū)間的解體碎片的狀態(tài)變化：

(6)

式中：si,j(r,t)是第i高度區(qū)間和j個面質(zhì)比區(qū)間內(nèi)解體碎片的空間密度；vri,j(r)是在大氣阻力作用下，第i高度區(qū)間和j個面質(zhì)比區(qū)間內(nèi)解體碎片軌道高度衰減的平均速度。從式(6)中可以看出，由于大氣阻力作用進(jìn)入和離開第i高度區(qū)的解體碎片數(shù)量的變化。

解體碎片的空間密度si,j(r,t)：

(7)

近地點地心距；分母中ak為碎片軌道半長軸。在大氣阻力攝動作用下，解體碎片的平均衰減速度vri,j(r)利用高斯型一般攝動力方程組[13]求解：

(8)

式中：fu、fh和fn分別表示大氣阻力沿解體碎片飛行速度方向、軌道面內(nèi)垂直速度方向以及軌道面法向上的分量。基于指數(shù)大氣密度模型，大氣阻力各分量取值為

(9)

考慮解體碎片運行在圓軌道上，綜合方程組式(8)和方程式(9)可以得到：

(10)

式中：CD是阻力系數(shù)；ρi0是參考橢球面r=ri0處的大氣密度；Hi是參考密度標(biāo)高。

令Δr→0，即研究極限情況下，同一高度上解體碎片的變化，得到描述解體碎片分布狀態(tài)變化的連續(xù)偏微分模型為

(11)

式(11)是一階擬線性偏微分方程。利用微分方程的特征線法，可以得到方程式(11)的解為

(12)

式中：εi,j定義為

(13)

式(12)是碎片云分布演化的解析表達(dá)式，描述了第i高度區(qū)間和j個面質(zhì)比區(qū)間內(nèi)的解體碎片，在不同高度上隨時間的演化分布情況?？梢钥闯?，解體碎片的分布主要受到碎片初始分布狀態(tài)和大氣阻力作用的影響。Nh個高度區(qū)間和Na個面質(zhì)比區(qū)間內(nèi)解體碎片分布的累加，即為解體碎片云的分布：

(14)

3　空間碎片云分布狀態(tài)的長期演化分析

在第1.2節(jié)中，利用數(shù)值積分計算方法，將1 780個解體碎片演化到環(huán)狀分布狀態(tài)。本節(jié)將利用第2節(jié)建立的解析演化模型，對處于環(huán)狀分布的碎片云的長期分布特點進(jìn)行分析。根據(jù)解體碎片的軌道高度和面質(zhì)比，將碎片云劃分為(Nh=210個高度區(qū)間)×(Na=10面質(zhì)比區(qū)間)=2 100個分組。表1給出了碎片云的初始分布狀態(tài)si,j(r,0)。

表1　2 100個碎片分組的初始空間密度si,j(r,0)

應(yīng)用解析分布演化模型式(14)，在配置為4核(單核主頻為2.93 GHz)的個人計算機上，0.2 s內(nèi)可以將碎片云的分布狀態(tài)演化到未來的50年，結(jié)果如圖8所示。分析長期分布演化結(jié)果可以得到以下結(jié)論：

1)碎片云的峰值密度在解體物體軌道高度附近，解體物體軌道高度附近的航天器將持續(xù)受到解體碎片云的碰撞威脅；

2)受大氣阻力的作用，碎片云峰值空間密度隨時間演化而減少，碎片云向更大高度區(qū)間內(nèi)擴散；

3)由于峰值密度處解體碎片不斷向低軌衰減，較低高度區(qū)間內(nèi)碎片在一段時期內(nèi)不斷增加，如演化10年時刻的碎片云分布曲線所示，隨后在大氣阻力衰減作用下較低高度區(qū)間內(nèi)碎片的密度也將不斷減少。

4　結(jié)束語

隨著空間碎片的不斷增加，航天器受到碎片的碰撞威脅將會日益嚴(yán)重。在軌工作航天器、失效載荷及火箭上面級等空間大物體，由于爆炸或相互碰撞，解體產(chǎn)生的大量碎片將成為未來空間碎片不斷增加的主要原因。本文提出的碎片云解析演化模型，描述了碎片云整體分布演化特點，能夠反映大氣阻力等主要攝動力對解體碎片云分布的影響。利用該模型，能夠快速地獲取碎片云的長期演化狀態(tài)，從而為評估碎片云對航天器產(chǎn)生的碰撞威脅、分析空間物體解體作用對近地軌道碎片環(huán)境的長期影響提供長期的碎片分布狀態(tài)數(shù)據(jù)輸入。

References)

[1]PARDINI C，ANSELMO L. Review of past on-orbit collisions among cataloged objects and examination of the catastrophic fragmentation concept [J]. Acta Astronaut.，2014，100: 30-39.

[2]LIOU J C. Collision activities in the future orbital debris environment[J]. Advances in Space Research，2006，38(9): 2102-2106.

[3]WANG T. Analysis of debris from the collision of the Cosmos 2251 and the Iridium 33 satellites[J]. Science & Global Security，2010，18: 87-118.

[4]ZHANG Y L，WANG Z K. Space traffic safety management and control[J]. Intelligent Transportation Systems，IEEE Transactions on，2015，PP99: 1-4.

[5]JOHNSON N L，KRISKO P H，LIOU J C，et al. NASA′s new breakup model of evolve 4.0[J]. Advances in Space Research，2001，28(9): 1377-1384.

[6]PARDINI C，ANSELMO L. Assessing the risk of orbital debris impact[J]. Adv. Space Res.，1999，1(1): 59-80.

[7]SWINERD G G，BARROWS S P，CROWTHER R. Short-term debris risk to large satellite constellations[J]. J. Guid. Control Dynam.，1999，22(2)：291-295.

[8]LIOU J C，HALL D T，KRISKO P H，et al. LEGEND-a three-dimensional LEO-to-GEO debris evolutionary model[J]. Advances in Space Research，2004，34(5): 981-986.

[9]ROSSI A，ANSELMO L，PARDINI C，et al. The new space debris mitigation(SDM 4.0)long-term evolution code[C]∥Proceedings of the 5th European Conference on Space Debris，Noordwijk，The Netherlands，2009.

[10]NAZARENKO A I. Prediction of the space debris spatial distribution on the basis of the evolution equations[J]. Acta Astronaut.，2014，100(1)：47-56.

[11]高惠璇. 統(tǒng)計計算[M]. 北京：北京大學(xué)出版社，1995.

GAO H X. Statistical computation[M]. Beijing： Peking University Press，1995.

[12]MONTENBRUCK O，GILL E. Satellite orbits: models，methods and applications[M]. Berlin:Springer Science & Business Media，2012.

[13]郗曉寧，王威. 近地航天器軌道基礎(chǔ)[M]. 長沙：國防科技大學(xué)出版社，2003.

XI X N，WANG W. Fundamentals of near-Earth spacecraft orbit[M]. Changsha: National University of Defense Technology Press，2003.

(編輯：車曉玲)

Modeling and analysis on the long-term evolution of the space debris cloud

ZHANG Binbin1，WANG Zhaokui2，*，ZHANG Yulin1，2

1.College of Aerospace Science and Engineering，National University of Defense Technology，Changsha 410073，China 2.School of Aerospace Engineering，Tsinghua University，Beijing 100084，China

A space debris cloud consists of numerous breakup fragments which are originated from an instinct space object. As the breakup fragments are in a concentrated distribution，those spacecraft in the neighborhood of the breakup objects will face a serious threat of collision. In order to evaluate the long-term evolution of space debris，the propagation of space debris cloud for a short-term period was analyzed numerically.Based on the stable distribution characteristic of the debris cloud，fragments were divided into several groups according to their orbital heights and area-mass ratios. The long-term distribution under the perturbation of atmosphere drag was described by an analytic evolution model. Comparing with all the fragments propagating individually，the analytic model was more efficient and less time-consuming. The evolution of a breakup debris cloud consists 1 780 breakup fragments，and was obtaired for the next fifty years. Results show that the peak spatial density of the debris cloud is around the height of the breakup object and will spread into a broader height interval continuously,and the fragments in the lower height interval will firstly increase and then decrease under the perturbation of air drag.

space debris cloud；spatial density; breakup model；debris distribution；evolution model

10.16708/j.cnki.1000-758X.2016.0046

2016-03-31；

2016-07-01；錄用日期：2016-06-30；

時間：2016-08-0214:20:06

http:∥www.cnki.net/kcms/detail/11.1859.V.20160802.1420.003.html

國家自然科學(xué)基金(11572168)

張斌斌(1987-)，男，博士研究生，zhangbinbin10@126.com，研究方向為航天器動力學(xué)與控制、空間碎片演化建模與分析

王兆魁(1978-)，男，教授，wangzk@tsinghua.edu.cn，主要研究方向為航天器動力學(xué)與控制

V528

A

http:∥zgkj.cast.cn

引用格式：張斌斌，王兆魁，張育林. 空間物體解體碎片云的長期演化建模與分析[J].中國空間科學(xué)技術(shù)，2016，36(4)：1-8.

ZHANGBB，WANGZK，ZHANGYL.Modelingandanalysisonthelong-termevolutionofthespacedebriscloud[J].ChineseSpaceScienceandTechnology，2016，36(4)：1-8 (inChinese).