(1.航天工程大學 研究生院，北京 101416； 2.航天工程大學 太空安全研究中心，北京 101416；3.中國人民解放軍66350部隊，河北 保定 071000)

0 引言

地球靜止軌道(GEO)因其特殊軌道位置特點，已成為國際上眾多大型衛(wèi)星尤其是通信衛(wèi)星的集中部署區(qū)，衛(wèi)星和空間目標十分密集。截至2018年6月1日，美國空間監(jiān)視網公布地球靜止軌道航天器850個，還有許多空間目標以大橢圓軌道穿越該區(qū)域，另外大量無法觀測到的空間碎片也分布在該區(qū)域[1]。這些空間目標密集地運行在海拔相差不到400公里，赤道兩側維度15°以內的球面約束的空間中，一旦發(fā)生碰撞將嚴重惡化日益嚴重的地球靜止軌道環(huán)境[2]，對航天器構成嚴重威脅，甚至給地球靜止軌道帶來災難性的影響[3]。

國外在航天器解體特性方面的研究比較成熟，通過對在軌撞擊試驗與地面撞擊實驗的統(tǒng)計分析，Reynolds等人于20世紀90年代末，擬合建立了“NASA標準解體模型”，該模型應用于NASA的EVOLVE4.0和LEGEND等空間碎片環(huán)境模型[4]。后來，改進后用在ESA的MASTER2009空間碎片環(huán)境模型[5]。2013年，德國Ernst-Mach-Institute的F.Schafer等開展了簡單立方體衛(wèi)星撞擊解體的試驗和仿真[6]。

國內許多學者也對空間環(huán)境、碰撞風險分析等方面有深入研究。中國空氣動力研究與發(fā)展中心開展了多次超高速撞擊解體試驗，蘭勝威和柳森等人對比了國外航天器解體模型的發(fā)展歷程，并基于CSBM模型開發(fā)了航天器碰撞解體碎片分析軟件SFA2.0[7-9]，哈爾濱工業(yè)大學的龐寶君團隊對空間碎片環(huán)境預測算法等進行了研究[10-11]，國防科技大學的張斌斌和白顯宗等人研究了碰撞預警與態(tài)勢的相關問題[12-13]，航天工程大學李怡勇、王衛(wèi)杰等人對空間目標的解體碎片的擴散特性進行了相關研究并評估了航天器撞擊解體碎片的短期危害[14-17]。

本文基于NASA標準解體模型，對生成的碎片信息進行檢驗，選取符合要求的解體碎片數(shù)據(jù)，分析其演化規(guī)律。一旦GEO衛(wèi)星發(fā)生碰撞解體，為分析新產生的碎片對GEO區(qū)域其他航天器的碰撞風險分析提供支撐。

1 碰撞相對速度分析

靜止軌道空間碎片的軌道傾角在攝動作用下會在±15°[18-19]內發(fā)生周期變化，運行在GEO附近的衛(wèi)星如果與軌道傾角為15°的空間碎片發(fā)生碰撞，碰撞相對速度約為：802.6 m/s。俄羅斯的“閃電”通信衛(wèi)星，運行軌道400 km×40000 km，傾角63.4°，近地點幅角270°。對于來自該軌道的空間碎片，如果近地點幅角發(fā)生漂移，當近地點幅角超過330°時，該碎片會穿越GEO區(qū)域，當近地點幅角為345.037°時，該碎片軌道會與GEO交會，此時，與GEO區(qū)域的衛(wèi)星碰撞概率最大，假設在交會點處存在一個GEO衛(wèi)星，則碎片相對衛(wèi)星的速度為2980.6 m/s。

碰撞相對速的計算如圖1所示。碰撞相對速度確定的關鍵是根據(jù)開普勒軌道六根數(shù)得到交會位置的速度矢量，即將圖中的4個空間目標開普勒根數(shù)轉換成空間目標的位置和速度。

圖1 碰撞相對速度示意圖

第一，根據(jù)開普勒公式:

M0=E0-esinE0

將空間目標交會處的平近地點角轉化成偏近地點角，然后根據(jù):

計算真近地點角f；

第二，令：幅角u=ω+f，據(jù):

計算地心距；

第三，根據(jù)球面三角形公式計算空間目標在J2000坐標系下的位置(x,y,z)：

第四，據(jù)活力公式：

計算空間目標速度的大小，據(jù)空間目標角動量守恒：

h=r×v

計算位置矢量和速度矢量的夾角θ，故速度矢量與由地心指向升交點方向的夾角為：

uv=u+θ

第五，同第三步，據(jù)球面三角形公式，得速度矢量在J2000坐標系下三個坐標方向的分量：

對于第一步中的超越方程M0=E0-esinE0求解問題，有多種解法，常用的方法有簡單迭代法和牛頓法。簡單迭代法的迭代公式為：

Ek+1=M+esinEk,k∈N

牛頓法求解超越方程。定義函數(shù)：

f(E)=E-esinE-M

超越方程的根即為函數(shù)f(E)的零點。將f(E)泰勒展開，

因此，取迭代公式為：

無論是簡單迭代法還是牛頓法，迭代初值可以簡單選取，不妨取E0=M0，迭代終止條件為：|Ek+1-Ek|<ε，不妨取ε=1×10-12。經過驗證，牛頓法比簡單迭代法迭代次數(shù)要少，效率相對更高。故，選用牛頓法通過平近地點角求偏近地點角。

來自“閃電”軌道的空間碎片，當近地點幅角為345.037°時，該碎片軌道會與GEO交會，計算方法如下：

設空間碎片的近地點俯角取ω0時，其軌道與GEO存在交匯，則：當空間碎片的近地點角滿足條件ω0+f0=π時，空間碎片的星下點維度為0，即：此時空間碎片處在其軌道與GEO的交點處。據(jù)此可計算得到ω0=345.037°。

2 解體碎片生成建模

空間目標高速碰撞會產生大量的空間碎片，根據(jù)空間目標碰撞前的運動狀態(tài)和質量等參數(shù)，模擬生成的空間碎片的尺寸、面質比、速度增量等，建立碎片生成模型[20-21]。

2.1 NASA標準解體模型

撞擊產生大于等于特征尺寸的碎片的數(shù)量為：

N(Lc) = 0.1mtot0.75Lc-1.71

式中，特征尺寸Lc為x、y、z三者的均值，其中，碎片最大長度記為x，與最大長度垂直的最大的長度記為y，與x-y平面垂直的最大長度記為z，解體部分的質量：

對于某一特征尺寸的碎片其面質比A/M的對數(shù)χ=lg(A/M)與特征尺寸的對數(shù)λc=lg(Lc)滿足：

DA/M=α(λc)N(μ1(λc),σ1(λc),χ)+

(1-α(λc))N(μ2(λc),σ2(λc),χ)

其中:N為正太分布，其概率密度函數(shù)為：

對于航天器，相關系數(shù)為：

對于特征尺寸小于1.67 mm的碎片，假設其形狀為正方體；對于特征尺寸大于1.67 mm的碎片，假設其為正方形。碎片平均截面積與特征尺寸之間的函數(shù)關系為：

撞擊產生碎片的分離速度方向為全向均勻分布，分離速度增量的對數(shù)δ=lg(△v)滿足：

D△v(χ,δ)=N(μ(χ),σ(χ),δ)

其中：

2.2 解體碎片仿真生成

碰撞解體事件遵循質量守恒定律，即解體后系統(tǒng)質量和等于解體前的質量和。碰撞解體事件遵循動量守恒定律，即解體后系統(tǒng)動量和等于解體前的動量和。

在對某一碰撞解體事件進行仿真時，仿真生成的數(shù)據(jù)具有隨機性，為了增加仿真數(shù)據(jù)的可靠性，可對仿真生成的數(shù)據(jù)運用質量守恒定律進行檢驗，如果仿真生成的數(shù)據(jù)不滿足質量守恒定律，則重新調用該模型，生成碎片數(shù)據(jù)，直到滿足質量守恒定律。運用質量守恒定律對NASA標準解體模型檢驗的流程如圖2所示。

圖2 運用質量守恒對NASA標準解體模型驗證

將碰撞解體碎片的速度分解為原始速度和由于碰撞獲得的速度增量兩者之和。碰撞解體事件滿足動量守恒一方面體現(xiàn)在碰撞解體碎片的原始速度v0的確定上，另一方面體現(xiàn)在碰撞解體后的系統(tǒng)的動量之和等于碰撞解體前系統(tǒng)的動量之和。該模型中，設定解體產生碎片的速度增量的方向滿足全向均勻分布，對于完全解體事件，解體后的初始速度可通過動量守恒定律確定，即：

(mt+mp)v0=mtvt+mpvp

對于非完全解體事件，假定解體產生碎片的質量來源按照碰撞解體前兩空間目標質量的加權獲得，即解體部分來自兩空間目標的質量分別為：

未解體部分的速度保持不變，解體產生的碎片的初始速度參照完全解體情形確定：

針對生成的仿真數(shù)據(jù)，無論是否完全解體，分別計算解體事件前后的動量和，如果不符合動量守恒定律，則舍棄重新運用該模型生成仿真數(shù)據(jù)，直到滿足動量守恒。

3 解體碎片演變

空間目標碰撞產生的碎片會迅速向外擴張，形成碎片云，以碎片云的密度和形狀作為參數(shù)描述演變過程，可將其分為球形、橢球形、繩形、螺旋線形、全方位擴散形和球殼形6個階段。[22]

1)球形階段。該階段是解體初期碎片急劇向外擴張的過程。空間目標碰撞解體產生碎片獲得的速度增量的方向服從各項均勻分布，在分離速度的作用下，空間碎片釋放初期呈現(xiàn)球形。解體點為初始時刻空間碎片的球心，球心速度可通過動量守恒獲得，碎片空間密度期望的峰值出現(xiàn)在一個球面上，該球面以解體速度增量的期望沿球徑向擴散。該階段，影響碎片云形狀的主要因素為碎片的解體速度增量。

2)橢球形階段。這一階段碎片云在空間目標間的相對運動規(guī)律作用下碎片云的形狀由球形變成橢球型。定義球形和橢球形階段的分界點為0.05個軌道周期，橢球形階段的起始時刻為：

3)繩形階段。該階段是指碎片云的形狀沿橢球的半長軸拉伸，形成兩頭尖中間粗的繩形的過程。當分離速度的作用不能增大橢球在軌道法向和徑向上的軸距，該軸距完全由相對運動規(guī)律決定時，橢球形階段結束，繩形階段開始。碎片在0.25個周期內沿各個方向的位移發(fā)生一次相位轉換，比如初始沿軌道切向的位移最大的碎片，在經過0.25個周期后，位移變?yōu)?，此即為繩形階段的起始時刻：

在此階段，由于繩頭和繩尾處碎片的軌道半長軸不同，他們之間的相對距離會越來越大。

4)螺旋線階段。繩形階段中，繩頭與繩尾間的距離逐漸增大，當繩頭追上繩尾，即繩頭處的碎片軌道運行超過繩尾處碎片一圈，此時繩形階段結束，開始螺旋線階段。

設解體前空間目標的軌道半長軸a0，角速度n0，繩頭處的碎片軌道半長軸a1，角速度n1，繩尾處的碎片軌道半長軸a2，角速度n2，解體碎片獲得的最大速度增量△vmax，有：

a0-a1=a2-a0=°(a0)

則螺旋線階段開始的時刻為：

故：

5)全方位擴散階段。該階段是指碎片向空間各個方向不斷擴散的過程，攝動是主要作用力，一般認為碎片升交點赤經相差10°時進入全方位擴散階段。當升交點赤經漂移率最大的碎片追上升交點赤經漂移率最小的碎片，全方位擴散階段結束。全方位擴散結束時刻[10]：

圖3 情形一仿真實驗

6)球殼形階段。全方位擴散階段結束，碎片在球殼內分布相對均勻，空間密度隨時間變化不會發(fā)生較大變化，此即為球殼階段。

4 算例仿真

解體碎片云的形狀在演變過程中會經歷六個階段。全方位擴散階段和球殼階段的主要作用力是攝動，這兩個階段持續(xù)時間較長。螺旋線階段的開始標志著平經度漂移率最大的碎片追上平經度漂移率最小的碎片，即：空間碎片蔓延到整個GEO區(qū)域，在地球赤道上空形成空間碎片的“環(huán)帶”。螺旋線階段的開始標志著空間碎片對整個GEO區(qū)域的航天器產生潛在的威脅。

引入發(fā)生碰撞兩空間目標的質量比值k，k等于質量較小的空間目標與質量較大空間目標的質量之比：

顯然，k∈(0,1]。

根據(jù)第1節(jié)確定的碰撞相對速度，運用NASA標準解體模型進行分析，若以802.6 m/s的相對速度發(fā)生碰撞，當k小于0.1242時，空間目標碰撞發(fā)生非完全解體；則對任意質量比值k大于等于0.1242時，空間目標碰撞發(fā)生完全解體。若以2980.6 m/s的相對速度發(fā)生碰撞，當k小于0.0090時，空間目標碰撞發(fā)生非完全解體；則對任意質量比值k大于等于0.0090時，空間目標碰撞發(fā)生完全解體。

選取SSN編號為36106的美國商業(yè)通信衛(wèi)星與空間碎片發(fā)生碰撞進行分析。該衛(wèi)星于2009年11月30日(UTCG)發(fā)射，運行在地球靜止軌道，據(jù)2018年05月31日的TLE數(shù)據(jù)，2018年05月31日22:47.870星下點經度85°E，質量2550 kg。

4.1 完全解體仿真

4.1.1 情形一

假設該衛(wèi)星與來自“閃電”軌道的碎片發(fā)生碰撞，碰撞相對速度2980.6 m/s，當k大于等于0.0090時，發(fā)生完全解體，計算知，當碎片的質量大于等于22.90 kg時，發(fā)生完全解體。當前情形下，不妨取碎片質量25 kg進行仿真，2018年6月1日0時，空間碎片軌道根數(shù)如表1所示。

表1 “閃電”軌道空間碎片軌道根數(shù)

仿真4次解體事件，計算解體產生的每一個碎片的平經度漂移率，橫坐標表示碎片的編號，縱坐標表示碎片平經度漂移率，正值表示方向由西向東，負值表示方向由東向西，單位為：°/h，如圖3所示。

從圖中可以看出，該情形下的四次仿真，碎片的平經度漂移率普遍在5°/h以內，最大的接近15°/h，解體碎片云分別經過13.58 h、14.4 h、12.41 h和13.33 h后進入螺旋線階段，即碎片會在赤道上空附近繞地球形成環(huán)帶，標志著空間碎片對整個GEO區(qū)域的航天器產生潛在的威脅。

圖4 情形二仿真實驗

4.1.2 情形二

假設該衛(wèi)星與運行在GEO附近的碎片發(fā)生碰撞，碰撞相對速度802.6 m/s，當k大于等于0.1242時，發(fā)生完全解體，計算知，當碎片的質量大于等于316.45 kg時，發(fā)生完全解體。當前情形下，取碎片質量320 kg進行仿真，2018年6月1日0時，空間碎片軌道根數(shù)如表2所示。

a/mei/(°)Ω/(°)ω/(°)M/(°)265560000.677863.4336.600

仿真4次解體事件，計算解體產生的每一個碎片的平經度漂移率，橫坐標表示碎片的編號，縱坐標表示碎片平經度漂移率，正值表示方向由西向東，負值表示方向由東向西，單位為：°/h，如圖4所示。

從圖4可以看出，該情形下的四次仿真，碎片的平經度漂移率普遍在5°/h以內，最大的達到25°/h，解體碎片云分別經過8.57 h、14.12 h、14.4 h和15.65 h后進入螺旋線階段，即碎片會在赤道上空附近繞地球形成環(huán)帶，標志著空間碎片對整個GEO區(qū)域的航天器產生潛在的威脅。

4.2 非完全解體仿真

4.2.1 情形三

假設該衛(wèi)星與來自“閃電”軌道的碎片發(fā)生碰撞，碰撞相對速度2 980.6 m/s，由2.6.1.1的分析知，當碎片的質量小于22.90 kg時，發(fā)生非完全解體。當前情形下，取碎片質量10 kg進行仿真，2018年6月1日0時，空間碎片軌道根數(shù)如表所示。

仿真4次解體事件，計算解體產生的每一個碎片的平經度漂移率，橫坐標表示碎片的編號，縱坐標表示碎片平經度漂移率，正值表示方向由西向東，負值表示方向由東向西，單位為：°/h，如圖5所示。

從圖6可以看出，該情形下的4次仿真，碎片的平經度漂移率普遍在2°/h以內，較完全解體的兩種情形漂移率小，解體碎片云分別經過60 h、12 h、25.71 h和32.73 h后進入螺旋線階段，即碎片會在赤道上空附近繞地球形成環(huán)帶，標志著空間碎片對整個GEO區(qū)域的航天器產生潛在的威脅，不過該情形產生的空間碎片的數(shù)量較完全解體的情形要少的多。

4.2.2 情形四

假設該衛(wèi)星與運行在GEO附近的碎片發(fā)生碰撞，碰撞相對速度802.6 m/s，由2.6.1.2的分析知，當碎片的質量小于316.45 kg時，發(fā)生非完全解體。當前情形下，取碎片質量10 kg進行仿真，2018年6月1日0時，空間碎片軌道根數(shù)如表1所示。

仿真4次解體事件，計算解體產生的每一個碎片的平經度漂移率，橫坐標表示碎片的編號，縱坐標表示碎片平經度漂移率，正值表示方向由西向東，負值表示方向由東向西，單位為：°/h，如圖6所示。

從圖6可以看出，碎片的平經度漂移率普遍在1.5°/h以內，4種情形中該情形下碎片的平經度漂移率最小。該情形下的4次仿真，解體碎片云分別經過51.43 h、83.72 h、72 h和102.86 h后進入螺旋線階段，即碎片會在赤道上空附近繞地球形成環(huán)帶，標志著空間碎片對整個GEO區(qū)域的航天器產生潛在的威脅，但該情形產生的空間碎片的數(shù)量較前三種情形要少的多，對GEO區(qū)域航天器的威脅要小的多。

圖6 情形四的仿真實驗

5 結論

通過第4節(jié)4種情形16次仿真實驗的結果來看，有14次實驗新產生的碎片在3天內擴散到整個GEO區(qū)域，另外兩次實驗，新產生的空間碎片也近乎擴散到整個GEO區(qū)域。對于來自“閃電”軌道的空間碎片與GEO衛(wèi)星發(fā)生碰撞的情形，兩者相對速度約2.98 km/s，相對速度較大，通過仿真實驗可知，碰撞極易發(fā)生完全解體，產生大量空間碎片，新產生的空間碎片迅速漂移到整個GEO區(qū)域。

對于來自GEO區(qū)域的空間碎片與GEO衛(wèi)星發(fā)生碰撞的情形，兩者相對速度較小，最大碰撞相對速度約0.8 km/s，只有較大質量的空間目標與GEO衛(wèi)星發(fā)生碰撞，才會導致衛(wèi)星完全解體，比如情形二中，只有當空間碎片質量達到320 kg的時候，與編號36106的航天器發(fā)生碰撞才會產生大量空間碎片。通過仿真實驗可知，較小質量的空間碎片與GEO衛(wèi)星發(fā)生碰撞產生的空間碎片對GEO區(qū)域的短期影響相對較小，不過產生的少量空間碎片對GEO區(qū)域的長期影響有待進一步研究。