地球同步軌道衛(wèi)星結構爆炸解體碎片的軌道仿真

周宵燈，崔村燕，趙蓓蕾，詹 翔，王 巖

(1. 航天工程大學 研究生管理大隊,北京 101416； 2. 航天工程大學 宇航技術系，北京 101416)

0 引言

當前，空間碎片已成為全球共同面臨的空間環(huán)境污染問題。為加強空間碎片管理，美國航空航天局(NASA)依據(jù)空間監(jiān)視網(wǎng)(SSN)的觀測數(shù)據(jù)和軌道測量結果對碎片進行編目，并形成與之配套的“兩行軌道根數(shù)”(TLE)數(shù)據(jù)。但美國SSN所用的雷達和望遠鏡的觀測能力有限，在低軌(LEO)范圍內可觀測和編目的空間目標的最小尺寸約為10 cm，在地球靜止軌道(GEO)上約為1 m。此外，盡管國外已建立了EVOLVE、ORDEM、MASTER[1-5]等多種空間碎片環(huán)境模型，但建模的數(shù)據(jù)主要源于SSN編目、雷達探測、望遠鏡觀測、航天器回收表面分析等10類探測數(shù)據(jù)，缺乏cm級碎片的直接探測數(shù)據(jù)，只能通過μm級(由長期暴軌裝置LDEF獲得)和cm級以上碎片數(shù)據(jù)插值獲得[6]。國內尚無公開的具有國際影響力的空間碎片環(huán)境模型，相關研究基本圍繞碎片探測、空間目標碰撞、軌道優(yōu)化等方面展開[7-10]。

衛(wèi)星爆炸解體產(chǎn)生的碎片一直是尺寸從幾毫米至幾分米的空間碎片的主要來源，這一尺寸范圍內的碎片包括直徑為1～10 cm的特別危險的碎片，現(xiàn)有的在軌技術無法對其進行防御，監(jiān)視網(wǎng)也無法跟蹤。因此，開展對衛(wèi)星爆炸碎片軌道的研究，快速掌握爆炸碎片軌道分布特性，特別是難以探測的cm級碎片的運行軌跡極具現(xiàn)實意義。

本文基于ANSYS/AUTODYN軟件對模擬衛(wèi)星結構進行爆炸數(shù)值仿真，將仿真得到的碎片參數(shù)作為初始數(shù)據(jù)，利用MATLAB將大量數(shù)據(jù)導入STK的SGP4模型中，得到所有碎片的軌道信息，并從碎片的軌道分布、速度增量、軌道演化、壽命估計等方面展開分析，以期為快速捕獲衛(wèi)星爆炸碎片，及時規(guī)避航天器碰撞提供參考。

1 衛(wèi)星爆炸仿真

真實衛(wèi)星形式多樣，結構復雜，爆炸解體過程和碎片分布特性受諸多因素影響。為簡化分析計算過程，本文以在軌柱形衛(wèi)星為研究對象。根據(jù)能量守恒定律，將衛(wèi)星推進劑等效為TNT炸藥，推進劑液態(tài)肼爆炸的TNT當量系數(shù)Y為0.93[11]，轉化后的TNT炸藥質量為164 g；根據(jù)結構強度等效方法將衛(wèi)星殼體等效為厚度為3 mm，殼體尺寸Φ為100 mm×150 mm，殼體材料為鋁合金的薄殼圓柱；選取位于圓柱中心線上距離端面10 mm處的衛(wèi)星推進劑的質心為炸點。

AUTODYN是一個顯式有限元分析程序，由ANSYS子公司研發(fā)，在國際軍工市場占比超過80%，為爆炸的應用案例求解提供了更高的精確度[12]。本文采用AUTODYN 15.0進行數(shù)值仿真，數(shù)值模型由TNT炸藥和鋁合金殼體2部分組成，在對爆炸碎片形成過程進行模擬時采用SPH算法，選擇粒子間距為1 mm。

柱形衛(wèi)星因內部推進劑爆炸而膨脹與破碎的過程如圖1所示。衛(wèi)星爆炸仿真最終生成了533塊質量大于10 mg的碎片，其中質量小于100 mg的碎片占大部分，具體分布如圖2所示。

圖1 衛(wèi)星殼體膨脹與破碎過程Fig.1 Expansion and fragmentation of satellite shell

圖2 碎片質量分布Fig.2 Mass distribution of debris

2 碎片軌道仿真

SGP4模型是在美國廣泛使用的標準外推模型，與TLE一起使用。它考慮了由地球非球形引起的長期和周期變化、日月引力、引力場共振效應，但并未充分考慮大氣阻力、太陽風等空間環(huán)境因素，只是簡單地使用一個阻力模型來計算產(chǎn)生的軌道衰變，導致計算結果精度不高，在進行軌道演化和壽命分析時，只能用作大致參考。SGP4模型為解析模型，優(yōu)點是計算速度快，本文選擇該模型的目的是在第一時間獲得爆炸碎片的軌道概況。SGP4模型是一種高效率的通用算法，但只有當輸入為TLE軌道參數(shù)時，計算結果才有效。TLE文件由美國戰(zhàn)略空間司令部開發(fā)，是以SGP4模型對空間目標進行軌道預報的數(shù)據(jù)文件格式，文獻[13]對其格式進行了詳細說明。

本文以地球同步軌道衛(wèi)星的爆炸仿真數(shù)據(jù)作為爆炸碎片的初始參數(shù)，其中包括碎片位置、速度等關鍵信息。AUTODYN軟件可輸出碎片在x、y、z方向上的速度值，通過與地面爆炸破片速度計算公式(即Gurney公式)進行對比可驗證仿真結果的準確性。相比于碎片與地球之間的距離，爆炸能量使碎片產(chǎn)生的初始位移可忽略不計。通過假設碎片的初始位置不變，確定碎片的直角坐標，再通過公式將坐標轉化為軌道六要素，得到碎片的初始TLE軌道參數(shù)。整個分析過程如圖3所示。首先利用拉格朗日法對某簡化柱形衛(wèi)星進行爆炸仿真，通過后處理得到衛(wèi)星爆炸后碎片的初始位置和速度參數(shù)，然后通過與經(jīng)驗公式作對比驗證爆炸仿真的準確性，如果仿真結果與理論計算相一致，則繼續(xù)下一步，否則重新設置參數(shù)進行衛(wèi)星爆炸仿真。通過MATLAB編程，將爆炸碎片初始位置和速度參數(shù)轉化為TLE數(shù)據(jù)文件格式，將其作為輸入導入SGP4模型，得到該衛(wèi)星爆炸后碎片的初始軌道。

圖3 衛(wèi)星爆炸碎片分析流程Fig.3 Analysis process of satellite explosion debris

2.1 爆炸碎片軌道分布

將TLE文件導入到SGP4模型后得到的碎片軌道分布情況如圖4所示。其中：紅色線條為目前正在運行的GPS導航衛(wèi)星的軌道，黑色線條為爆炸碎片的運行軌道。在衛(wèi)星和碎片運動軌道參數(shù)已知的條件下，借助STK軟件，就可實現(xiàn)實時的軌道預警。在可能發(fā)生碰撞的位置，衛(wèi)星采取軌道機動的方式，避免與碎片發(fā)生碰撞。

圖4 碎片軌道分布圖Fig.4 Debris orbit distribution

為進一步分析爆炸碎片的軌道分布情況，選取爆炸碎片中的4個典型單元進行模擬研究。因爆炸能量使碎片產(chǎn)生的初始位移可忽略不計，故碎片的位置參數(shù)可選用衛(wèi)星發(fā)生爆炸時的位置(42 154.8，-982.731，0)，單位為km。衛(wèi)星和各碎片單元在該位置的速度參數(shù)見表1。

表1 衛(wèi)星和各單元速度參數(shù)

將經(jīng)轉化后的數(shù)據(jù)代入模型，得到衛(wèi)星與各單元的運行軌道，如圖5所示。結合軌道仿真結果和力學分析，可將爆炸產(chǎn)生的碎片運行軌跡分為4種類型：當碎片速度小于衛(wèi)星運行速度且小于某一臨界值V1時，碎片就會像A單元一樣落入大氣層，不會變成太空垃圾污染空間環(huán)境；當碎片速度小于衛(wèi)星速度但大于臨界值V1時，碎片的運行軌跡如B單元，碎片繞著地球做橢圓運動，但運行范圍小于原來衛(wèi)星的范圍；當碎片速度大于衛(wèi)星運行速度且大于某一臨界值V2時，碎片就會像D單元一樣脫離地球引力飛向更遠的地方；當碎片速度大于衛(wèi)星運行速度且小于某一臨界值V2時，碎片的運行軌跡如C單元，碎片繞著地球做橢圓運動，但運行范圍大于原來衛(wèi)星的范圍。

圖5 衛(wèi)星與各單元運行軌道Fig.5 Orbits of satellite and each unit

通過活力公式可求得臨界值V1和V2的理論值，將其與仿真結果進行對比，對比結果見表2。

由表可見，數(shù)值仿真與經(jīng)驗估算的結果非常接近，其相對誤差在允許范圍之內。數(shù)值仿真結果之所以略大于公式推導結果，是因為在公式推導中，未考慮軌道攝動力的影響。

表2 臨界值V1、V2數(shù)值仿真與公式推導結果對照

2.2 爆炸碎片速度增量分析

衛(wèi)星發(fā)生爆炸時，爆炸沖量使碎片獲得的速度增量為Δv，碎片速度的改變引起了軌道根數(shù)的變化。爆炸碎片的軌道根數(shù)變化與速度增量的關系為

式中：Δvx、Δvy、Δvz為軌道坐標系O-xyz各軸上的速度增量；Δa、Δe、Δω、ΔΩ、Δi、ΔM為軌道根數(shù)的變化，其中，a為半長軸，e為偏心率，ω為近地點幅角，Ω為升交點赤經(jīng)，i為軌道傾角，M為平近地點角；μ=3.986 005×1014m3/s2，為地球引力常數(shù)；r=(x2+y2+z2)1/2，為衛(wèi)星到地心的距離；u=1/r；p為近地點。由上式可知，爆炸碎片的速度增量Δvx和Δvy決定了平近地點角、偏心率和半長軸的變化，而Δvz則決定了升交點赤經(jīng)和軌道傾角的變化。

圖6 爆炸碎片在x，y，z三個方向上的速度增量頻次分布Fig.6 Velocity increment frequency distribution of explosion debris in three directions of x, y and z

以衛(wèi)星爆炸仿真得到的碎片參數(shù)為依據(jù)，分析爆炸碎片在x、y、z三個方向上的速度增量，結果如圖6所示。由圖可知，在x方向上的速度增量分布范圍為-200～300 m/s，其中在-25～25 m/s之間碎片分布最為集中，爆炸碎片在x方向上的速度增量均值為7.9 m/s，y方向上的速度增量均值為9.2 m/s，z方向上的速度增量均值為-9.6 m/s。

2.3 爆炸碎片軌道演化

軌道演化是一個動態(tài)過程，想要清楚地看到軌道的變化情況，需要一個較大的時間差。圖7是將STK中的場景時長從原來的1個月，增加到3年后的B單元的軌道情況。由圖可見，B單元的運行軌跡并非是一條線，而是一個區(qū)域，這是因為隨著運行時間的增加，原來的計算模型的誤差將會不斷積累，軌道會發(fā)生很大變化。

圖7 B單元運行軌道Fig.7 Operating orbit of unit B

3 結束語

本文利用MATLAB將AUTODYN爆炸仿真軟件與航天分析軟件STK聯(lián)系到一起。仿真結果表明：SGP4模型運算速度快，滿足大量爆炸碎片的軌道仿真的要求。通過應用軌道動力學和天體力學的相關知識分析碎片受力情況，可以發(fā)現(xiàn)，攝動分析的精度嚴重受到攝動力模型的影響，而不同仿真模型采用的內部攝動力模型也各不相同。因此，從理論與仿真兩方面進行分析驗證極為必要。本文對爆炸碎片的運行軌道進行了分類。在具體分析碎片運行對航天器的危害時，可去除落入大氣層和飛入僵尸軌道的碎片，以減少任務量，提高計算效率；可對碎片進行威脅評估、跟蹤和編目。

如何快速確定與預報衛(wèi)星在軌爆炸解體產(chǎn)生的大量空間碎片，以提供及時的空間碰撞預警是目前空間態(tài)勢感知的一個重要研究方向。本文提出的方法受SGP4模型限制，對攝動力模型進行了簡化，未充分考慮高低溫、太陽風等空間環(huán)境的影響，導致仿真結果精度不夠。后續(xù)研究可從模型本身出發(fā)，進一步考慮空間環(huán)境因素，應用更精確的大氣密度模型，利用地面觀測數(shù)據(jù)對分析軌道進行修正，以提高計算精度。