張育林，張斌斌，王兆魁

(1. 國防科學(xué)技術(shù)大學(xué)空天科學(xué)學(xué)院，長沙 410073；2. 清華大學(xué)航天航空學(xué)院，北京 100084)

0 引 言

空間碎片環(huán)境是在人類航天事業(yè)發(fā)展過程中才逐漸形成的一類“人造”空間環(huán)境，是由大量空間碎片構(gòu)成的。根據(jù)美國空間監(jiān)視網(wǎng)的觀測數(shù)據(jù)，截至2016年仍在軌運行的空間目標超過17000個，而在軌運行的編目目標中，近90%是空間碎片[1]。由于空間碎片的運動是無控且很難精確預(yù)測的，碎片一旦與航天器相撞，將會引起航天器故障甚至?xí)购教炱鹘怏w，給航天器在軌運行帶來巨大威脅。歷史上已經(jīng)發(fā)生多次航天器被空間碎片撞擊的事件[2-3]，如1991年俄羅斯COSMOS 1934衛(wèi)星被一個碎片撞擊，導(dǎo)致衛(wèi)星部分解體；2009年美國Iridium 33衛(wèi)星和一個大的完整空間碎片——俄羅斯失效的通信衛(wèi)星COSMOS 2251相撞，直接導(dǎo)致兩個衛(wèi)星解體，產(chǎn)生近2000個編目解體碎片和大量未編目的碎片[4]。另外，還有10起以上已經(jīng)確認由空間碎片碰撞而引起的航天器失效或故障事件[5]?？臻g碎片環(huán)境已經(jīng)對人類航天活動已產(chǎn)生嚴重的影響。

隨著空間碎片規(guī)模的持續(xù)增加，空間碎片環(huán)境帶來的空間安全問題日益突出。對碎片環(huán)境的長期演化趨勢進行預(yù)測和分析，獲取碎片碎片環(huán)境演化的內(nèi)在機理和整體分布指標，則是開展碎片環(huán)境評估以及制定碎片環(huán)境減緩策略的重要基礎(chǔ)[6]。早在1978年，Kessler和Cour-Palais等人，通過分析空間目標的相互碰撞作用，建立描述空間碎片數(shù)量變化的數(shù)學(xué)模型[7]，對碎片的增長趨勢進行分析。近些年，碎片環(huán)境演化建模研究的工作得到了快速發(fā)展。根據(jù)演化模型選取的演化狀態(tài)量，可將演化模型可分為兩類[8]：

一類是以宏觀狀態(tài)量為變量的整體演化模型。利用宏觀量作為描述碎片演化的狀態(tài)量，如碎片在空間中的分布密度、碎片的數(shù)量等，建立碎片環(huán)境整體演化模型，利用解析或數(shù)值方法得到碎片環(huán)境的長期演化規(guī)律。Rossi[9]、Telant[10]、Lewis[11]、Nazarenko[12]以及Kebschull[13]等人在整體演化模型方面開展了研究。當(dāng)選取空間碎片數(shù)量作為狀態(tài)量時，利用常微分方程(組)建立演化模型，能夠得到碎片數(shù)量的長期增長趨勢；當(dāng)選取空間密度作為狀態(tài)量時，利用偏微分方程(組)建立演化模型，能夠給出碎片在空間中的演化分布趨勢。由于選取宏觀量描述碎片的演化狀態(tài)，不跟蹤單個碎片的運動狀態(tài)，演化效率通常較高。另一方面，由于對碎片環(huán)境的演化進行了理想的近似假設(shè)，限制了整體演化模型結(jié)果的精度。如上述模型中，Rossi等人假設(shè)碎片運行在圓軌道上。整體演化模型中碎片受到的作用力模型一般較為簡單，如Rossi、Letizia等人在建立整體模型時僅考慮了大氣阻力作用。

另一類是以單個碎片的運動狀態(tài)為變量的演化模型。如低地球軌道到地球靜止軌道(GEO)環(huán)境的碎片模型(LEGEND)[14]、流星體和空間碎片環(huán)境參考模型(MASTER)[15]、空間碎片減緩長期分析模型(SDM)[16]、地球碎片環(huán)境演化模型(MEDEE)[17]、地球同步軌道碎片環(huán)境分析和監(jiān)視模型(DAMAGE)[18]以及長期碰撞分析工具(LUCA，LUCA-2)[19]等。上述模型中，首先利用確定性或隨機性方法，建立攝動力作用、航天發(fā)射活動、空間目標的相互碰撞解體以及在軌目標爆炸解體等影響因素的計算模型；再利用軌道更新計算方法，不斷對碎片的運動狀態(tài)進行推演，從而得到碎片環(huán)境的長期演化結(jié)果。通過跟蹤單個碎片的運動狀態(tài)，可以建立精確的力學(xué)模型來確定碎片的運動狀態(tài)。但跟蹤每一個空間碎片的運動狀態(tài)，當(dāng)碎片規(guī)模不斷增加時，需要消耗大量的計算資源才能得到碎片環(huán)境的長期演化結(jié)果。對于空間目標相互碰撞事件、目標在軌爆炸事件，一般通過蒙特卡羅隨機仿真的方法實現(xiàn)。基于上述模型運行得到的一次演化計算結(jié)果只是眾多隨機結(jié)果的一種可能，需要在多次計算結(jié)果的基礎(chǔ)上，統(tǒng)計得到期望的演化結(jié)果，這進一步提高了對演化計算資源的需求。因此，針對影響碎片環(huán)境演化的復(fù)雜因素，在建立精確可靠的計算模型的條件下，進一步設(shè)計高效的演化計算方法，是跟蹤單個碎片的運動狀態(tài)，構(gòu)建碎片環(huán)境演化計算模型的關(guān)鍵。

本文針對空間碎片環(huán)境的演化建模問題，從宏觀和微觀兩個方面，分別構(gòu)建了碎片環(huán)境的整體演化模型和數(shù)值演化計算模型。通過對比兩種模型的長期演化結(jié)果，分析了兩種模型的特點和適用條件。利用演化模型的結(jié)果，討論了碎片環(huán)境的長期演化趨勢和穩(wěn)定性。

1 空間碎片環(huán)境長期演化建模

1.1 分層離散化演化模型

針對低地球軌道碎片碎片環(huán)境，以碎片空間密度作為描述碎片環(huán)境分布狀態(tài)的宏觀量，考慮影響碎片環(huán)境演化的主要影響因素，建立空間碎片環(huán)境整體演化動力學(xué)方程。

首先，將碎片分布空間分層離散化。根據(jù)當(dāng)前編目目標的運動狀態(tài)[1]，近90%以上的編目目標運行在圓軌道上。且圓軌道目標解體產(chǎn)生的碎片中絕大部分運行在圓軌道上。為了方便討論，假設(shè)空間碎片均運行在圓軌道上，且碎片為勻質(zhì)球形?？蓪⒔剀壍郎系目臻g碎片劃分到Nh高度層內(nèi)。在攝動力作用下，認為空間碎片在同一個高度層內(nèi)是均勻分布的[9, 20]。為了體現(xiàn)大氣阻力對碎片環(huán)境的影響，進一步將同一個高度層內(nèi)的碎片，按照其尺寸大小劃分到Na個面質(zhì)比區(qū)間內(nèi)。綜上，我們可以將空間碎片分成Nh×Na個碎片組。

(1)

(2)

式中：(A/m)j是第j個面質(zhì)比區(qū)間的平均面質(zhì)比；ri0和Hi分別是第i個高度層的參考地心距和密度標高；ρi0是ri0處的參考大氣密度。對于不同高度層，ri0、ρi0和Hi取值不同，可以通過查表獲取，如參考文獻[21]。

為了簡化，進一步將式(2)中的系數(shù)定義為：

(3)

(4)

(5)

(6)

(7)

將方程(4)和(5)代入到約束控制方程(1)中，并令Δr足夠小，即令Δr→0，可以得到微分形式的分層離散化模型：

i=1,2,…,Nhj=1,2,…,Na

(8)

求解模型(8)中包含的Nh×Na個偏微分方程組，得到低地球軌道上空間碎片的長期演化結(jié)果。

1.2 空間碎片環(huán)境的長期演化計算模型

通過跟蹤逐個碎片的運動狀態(tài)，建立碎片演化過程中的攝動力作用、目標之間相互碰撞解體、航天發(fā)射活動等影響因素的計算模型，在高效數(shù)值計算方法的基礎(chǔ)上，構(gòu)建碎片環(huán)境的長期演化計算模型。長期演化計算模型的基本框架如圖2所示。

圖2中目標爆炸解體、目標相互碰撞解體以及人類的航天發(fā)射活動等會引起空間碎片規(guī)模不斷增加；而碎片在空間運動過程中的受力，尤其大氣阻力的作用，會引起碎片運行軌道高度降低，最終導(dǎo)致碎片再入大氣層銷毀。在多種影響因素綜合作用下，空間碎片數(shù)量將隨時間不斷動態(tài)變化。

長期演化計算模型中考慮的攝動力因素包括：大氣阻力，采用Harris-Priester模型計算大氣密度；基于4×4階重力場模型確定的地球非球形攝動；太陽、月球引力攝動；太陽光壓攝動。上述攝動力計算模型是在綜合計算精度和復(fù)雜度的基礎(chǔ)上確定的，確保能夠反映碎片環(huán)境長期分布規(guī)律，且具有較高的計算效率。

空間目標之間劇烈的碰撞作用會導(dǎo)致目標解體，產(chǎn)生大量的解體碎片，嚴重影響著空間碎片環(huán)境的長期演化。對空間碎片環(huán)境演化過程的碰撞事件進行建模研究，包含兩個方面，一方面是建立碰撞概率計算模型，確定演化過程可能發(fā)生碰撞的空間目標；另一方面是建立目標解體產(chǎn)生碎片的模擬方法，產(chǎn)生目標解體碎片。在確定目標之間發(fā)生碰撞的概率時，借鑒“立方體”碰撞概率模型[24]的離散化空間思想，即將空間碎片分布的空間劃分成離散的空間體積元，在每個體積元內(nèi)確定目標相互碰撞的概率。如圖3所示，沿著地心距、赤經(jīng)和赤緯方向，分別按照間隔Δh、Δλ和Δφ，將碎片分布空間劃分成若干空間體積元。空間體積元Ci,j,k定義為：

(9)

式中(ri,λj,φk)是體積元Ci,j,k中心的球坐標。

對于碎片環(huán)境演化過程的某一時刻，當(dāng)一個體積元目標數(shù)量大于1個時，通過對兩兩目標之間的相對運動狀態(tài)進行，確定碰撞概率。定義垂直于相對運動速度矢量的平面為相遇平面，在相遇平面內(nèi)將碰撞概率計算轉(zhuǎn)化為二維平面概率密度積分問題，進一步利用級數(shù)展開得到近似的解析解[25-26]。記同一個體積元內(nèi)兩目標的碰撞概率為PC，則

(10)

式中:u和v分別定義為

(11)

由于式(10)具有解析的表達式，提高了大規(guī)?？臻g目標之間的碰撞概率的計算效率。確定目標之間發(fā)生碰撞的概率后，可進一步利用NASA標準解體模型，通過隨機抽樣方法確定兩目標是否發(fā)生碰撞。

航天發(fā)射活動是導(dǎo)致空間大尺寸目標增加的唯一原因，是產(chǎn)生大量空間碎片的物質(zhì)基礎(chǔ)，對碎片環(huán)境的演化產(chǎn)生著重要影響。發(fā)射率L定義為每年發(fā)射入軌的航天數(shù)量，可用下式確定：

L=f(t,h)

(12)

式中:t是時間，單位是年；h是高度。對于給定的時刻和高度，L可取常值，也可依據(jù)歷史發(fā)射任務(wù)確定。本文選取2012～2016年的發(fā)射活動作為演化計算中連續(xù)5年的航天器發(fā)射活動，每隔5年重復(fù)一次該發(fā)射規(guī)模。即發(fā)射模型是時間的周期性函數(shù)，周期為5年。新發(fā)射入軌的航天器在不同高度的分布情況，由2012～2016年的航天器發(fā)射分布情況確定。2012～2016年間年發(fā)射規(guī)模在表1中給出。

表1 2012～2016年間的航天器發(fā)射情況Table 1 The spacecraft launch traffic over 2012～2016

進行碎片環(huán)境長期演化計算中，需要對碎片的運動狀態(tài)按照一定時間步長進行推演更新。由于空間碎片數(shù)量眾多，且在長期演化過程中不斷動態(tài)變化，對碎片的狀態(tài)推演更新方法的效率和精度均提出了較高的要求。在演化計算模型中，通過對碎片的運動狀態(tài)在一個軌道周期內(nèi)求平均，分離出短周期運動部分，得到碎片運行的平均軌道。建立攝動力作用下，描述平均軌道變化率的運動方程。由于分離了短周期運動部分，可以采用較大的積分步長，如可將積分步長設(shè)置為1天，能夠大大提高碎片運動狀態(tài)的演化計算效率[8]。

最后，基于MPI(Massage Passage Interface)并行標準[27]，設(shè)計了大規(guī)?？臻g碎片演化的并行計算框架，能夠高效利用“天河一號”等超級計算機平臺的計算資源進行大碎片環(huán)境的長期演化。

2 模型求解

2.1 簡化分層離散化模型的解析解

對分層離散化模型進行簡化，忽略碎片之間的相互碰撞作用，分析在大氣阻力作用下空間碎片環(huán)境的長期演化特點。將碰撞相關(guān)項從方程(8)中去除，整理可以得到僅有大氣阻力作用下，碎片環(huán)境演化的約束方程：

(13)

(14)

(15)

(16)

(17)

進一步求解方程(16)，得到簡化分層離散化模型的通解為：

(18)

式中:變量α利用(17)來確定。

i=1,2,…,Nhj=1,2,…,Na

(19)

2.2 分層離散化模型的數(shù)值求解

當(dāng)考慮空間目標之間的相互碰撞作用時，可通過求解方程(8)得到碎片環(huán)境的演化狀態(tài)。將方程(8)改寫成式(19)的形式。可以看出，考慮目標相互碰撞時的分層離散化模型由Nh×Na個高度耦合的一階線性雙曲型偏微分方程構(gòu)成，可通過逼近雙曲方程的差分格式進行數(shù)值求解。將改寫成標準的雙曲型形式為：

(20)

(21)

式中

(22)

式中:τ和h分別是差分逼近時的離散時間間隔和離散高度間隔，下標n表示第n個離散時間步長。容易證明，當(dāng)滿足條件約束條件(23)：

(23)

差分格式(21)可以對方程(19)實現(xiàn)穩(wěn)定逼近。需要說明的是，為了確保約束控制方程(19)成立，在選擇離散時間間隔和離散高度間隔時，應(yīng)保證每一步計算時，同一高度層內(nèi)的碎片最多衰減到下一個高度層內(nèi)。

2.3 長期演化計算模型的蒙特卡羅仿真實現(xiàn)方法

(24)

3 碎片環(huán)境的長期演化結(jié)果與分析

低地球軌道空間是當(dāng)前目標分布最為密集的區(qū)域。本節(jié)以低地球軌道碎片環(huán)境為研究對象，利用分層離散化模型和長期演化計算模型，開展碎片環(huán)境長期演化分析。這里選取低地球軌道碎片環(huán)境的空間分布范圍是距地面高度在200～3000 km的空間。僅考慮尺寸大于10 cm的空間碎片，并假設(shè)航天器在軌工作期間的軌道構(gòu)型不發(fā)生變化，不進行碰撞規(guī)避機動。

演化計算中，空間目標的初始狀態(tài)利用觀測數(shù)據(jù)確定。根據(jù)空間監(jiān)視網(wǎng)(SSN)的觀測數(shù)據(jù)，截至到2016年9月，低地球軌道上運行著12000個編目目標。將低地球軌道空間，按10 km間隔劃分為280個高度層；對于同一個高度層內(nèi)的碎片，按照等對數(shù)間隔的原則將碎片劃分到10個尺寸區(qū)間內(nèi)，即要求各尺寸區(qū)間邊界對數(shù)值的差值相等。根據(jù)1.1節(jié)中的分層離散化方案，在表2中，將低地球軌道上的空間碎片劃分到2800個碎片組，并給出了每個碎片組的初始空間密度。

表2 編目目標的初始空間密度分布(1/km3)Table 2 Initial spatial density of the catalogued objects (1/km3)

3.1 攝動力作用下碎片環(huán)境的自我凈化能力分析

僅考慮動力學(xué)約束，對碎片環(huán)境進行長期演化，可用來分析分層離散化模型和演化計算模型的力學(xué)條件是否一致，并討論碎片環(huán)境的自我凈化能力。

僅有攝動力的作用下，基于解析的分層離散化演化模型即方程來獲取碎片環(huán)境的長期演化結(jié)果。碎片環(huán)境的初始狀態(tài)si,j(r,0)在表2中給出。針對尺寸大于10 cm的碎片在未來200年內(nèi)的演化增長趨勢和分布狀態(tài)分別如圖4和圖5所示。

從圖4可以看出，僅有考慮攝動力的作用下，經(jīng)過200年的演化，碎片的規(guī)模從從初始12000個減少到5800多個。進一步從圖5可以看出，位于1000 km高度以下碎片的密度減少趨勢較為明顯，1000 km高度以上碎片的空間密度沒有明顯變化。

作為對比，利用長期演化計算模型，對低地球軌道碎片環(huán)境未來200年內(nèi)的狀態(tài)進行演化分析，如圖6和圖7所示的碎片環(huán)境演化結(jié)果。對比圖4和圖6、圖5和圖7可以看出，僅考慮攝動力作用下，利用兩種模型計算得到的碎片數(shù)量的減少趨勢、碎片在空間中的分布狀態(tài)的變化都是一致的。綜合兩種模型的演化結(jié)果可以看出，碎片數(shù)量在未來200年內(nèi)減少了約50%。這說明在攝動力的作用下空間碎片環(huán)境具有一定的自我凈化能力，但這種自然清除過程是緩慢的。

兩種模型均能給出碎片環(huán)境的長期演化結(jié)果，但兩種模型所需要的計算量相差較大。表3給出了兩種模型所需要計算量的統(tǒng)計，可以看出，長期演化計算模型需要的計算量大約是分層離散化模型所需計算量的104倍。

因此，分層離散化模型適合在計算資源有限的條件下，以較少的計算量，得到碎片規(guī)模的長期增長趨勢，從整體上把握碎片環(huán)境的演化特點。長期演化計算模型中包含的攝動力因素較為全面，能夠跟蹤空間目標之間發(fā)生的每一次碰撞事件，適合在具有大量計算資源的條件下，對碎片環(huán)境的增長規(guī)模、空間分布狀態(tài)進行精確、可靠的分析。

表3 僅考慮攝動力作用時兩種演化模型的計算量統(tǒng)計Table 3 Compute assumptions of the two evolution models

3.2 考慮目標相互碰撞作用時碎片環(huán)境的穩(wěn)定性分析

假設(shè)一種理想演化條件，即停止一切未來的航天發(fā)射活動，考慮碎片受到的作用力和空間目標之間的相互碰撞作用，分析當(dāng)前碎片環(huán)境的穩(wěn)定性。首先，定義碎片環(huán)境的兩種長期演化狀態(tài)：(1)穩(wěn)定演化狀態(tài)。處于穩(wěn)定演化狀態(tài)的空間碎片環(huán)境，空間碎片數(shù)量保持不變，甚至在不斷減少。(2)不穩(wěn)定的演化狀態(tài)。當(dāng)碎片環(huán)境處于不穩(wěn)定演化狀態(tài)時，空間碎片的數(shù)量在持續(xù)增加。

從理想條件下碎片數(shù)量的增長趨勢看以看出，近地軌道上總碎片數(shù)量隨演化時間在不斷增加，200年內(nèi)增加為初始碎片數(shù)量的近3倍，如圖8所示。其次，從碎片的來源看，空間碎片中碰撞解體碎片所占的比例在逐年增加，演化結(jié)束時80%以上的碎片均為碰撞解體碎片。因此，即使停止一切航天發(fā)射活動，大量解體碎片仍將促使總的碎片數(shù)量持續(xù)增加，表明當(dāng)前低地球軌道上的空間碎片處于不穩(wěn)定演化狀態(tài)。

圖10給出了，從長期演化模型計算結(jié)果中，統(tǒng)計得到的目標累積災(zāi)難性碰撞次數(shù)隨時間變化曲線?？梢钥闯觯?00年的演化時間內(nèi)，目標之間發(fā)生了近30次災(zāi)難性碰撞，累積碰撞次數(shù)隨演化年份呈現(xiàn)線性增長的趨勢，這與LEGEND模型的演化結(jié)果是一致的[29]。

表4給出了理想演化條件下不同尺寸目標之間發(fā)生碰撞的次數(shù)。統(tǒng)計時若包絡(luò)尺寸半徑超過1 m，則目標為空間大尺寸目標，包絡(luò)尺寸半徑小于1 m的為小尺寸目標?？梢钥闯?，近三分之二的碰撞解體發(fā)生在大尺寸目標和小尺寸目標之間，大尺寸目標和小尺寸目標發(fā)生碰撞次數(shù)是大尺寸目標之間發(fā)生碰撞次數(shù)的3倍多。因此，大尺寸目標解體會產(chǎn)生大量解體碎片，解體碎片會進一步與大尺寸目標發(fā)生碰撞，產(chǎn)生更多的解體碎片，這種“碰撞—解體碎片—碰撞”反饋連鎖碰撞效應(yīng)，使得解體碎片數(shù)量迅速增加。

3.3 航天發(fā)射活動對碎片環(huán)境演化的長期影響

針對低地球軌道碎片環(huán)境，利用航天發(fā)射活動模型，分析發(fā)射活動對空間碎片環(huán)境的影響。通過對碎片環(huán)境的20次演化結(jié)果進行統(tǒng)計平均，得到200年的演化結(jié)果如圖11和圖12所示。考慮航天發(fā)射活動后，空間碎片數(shù)量急劇增加，在200年內(nèi)尺寸大于10 cm碎片數(shù)量增加到約16萬。碰撞解體碎片占總碎片比例也在逐年增加，至演化結(jié)束時刻，95%以上空間碎片是由目標碰撞解體產(chǎn)生的，如圖11所示。

表4 理想演化條件下不同尺寸目標發(fā)生劇烈碰撞的次數(shù)Table 4 The catastrophic collision numbers between different size of space objects under the ideal condition

結(jié)合劇烈碰撞次數(shù)的統(tǒng)計結(jié)果，如圖13和圖14所示，可以看出，由于航天發(fā)射活動在不斷向空間中引入大尺寸空間目標，“碰撞-解體碎片-碰撞”的反饋連鎖碰撞效應(yīng)加強，導(dǎo)致近地軌道上空間碰撞解體碎片數(shù)量迅速增加，這在[500, 800] km高度范圍內(nèi)表現(xiàn)尤為明顯，如圖14所示。

表5給出了不同尺寸目標之間發(fā)生劇烈碰撞的次數(shù)，與理想演化條件下的演化結(jié)果相比，大尺寸目標之間發(fā)生碰撞的次數(shù)、大尺寸目標與小尺寸目標發(fā)生碰撞的次數(shù)都明顯增加了，進一步說明航天發(fā)射活動促進了碎片環(huán)境的連鎖碰撞效應(yīng)。

表5 考慮航天發(fā)射活動時不同尺寸目標發(fā)生碰撞的次數(shù)Table 5 The catastrophic collision numbers between different size of space objects with spacecraft launch traffic

4 結(jié) 論

隨著空間碎片數(shù)量的不斷增加，空間碎片在占據(jù)有限軌道資源的同時，也使在軌工作航天器面臨著日益嚴峻的碰撞風(fēng)險，空間碎片環(huán)境已嚴重影響了人類對空間資源的可持續(xù)開發(fā)利用。本文緊跟空間碎片環(huán)境研究這一熱點問題，對空間碎片環(huán)境的長期演化建模問題進行了研究，建立了空間碎片環(huán)境的長期演化計算模型和分層離散化演化模型，并討論了碎片環(huán)境的長期演化趨勢和穩(wěn)定性。結(jié)果表明，空間目標的相互碰撞解體，是空間碎片不斷增加的主要因素；即使停止一切航天發(fā)射活動，空間碎片的數(shù)量仍在不斷增加，表明低地球軌道空間碎片規(guī)模已經(jīng)超越穩(wěn)定臨界點；如不采取有效的碎片減緩策略，未來航天發(fā)射活動會進一步增強空間碎片環(huán)境演化的不穩(wěn)定性，加劇“碰撞-目標解體-碰撞”反饋連鎖碰撞效應(yīng)。論文的研究為掌握碎片環(huán)境的長期演化機理奠定了理論基礎(chǔ)，提供了有效的分析工具；將為航天器的安全運行和空間資源的可持續(xù)開發(fā)利用提供重要支撐。