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

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

0 引言

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

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

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

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

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

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

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

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

圖2 碎片質(zhì)量分布Fig.2 Mass distribution of debris

2 碎片軌道仿真

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

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

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

2.1 爆炸碎片軌道分布

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

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

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

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

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

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

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

由表可見，數(shù)值仿真與經(jīng)驗(yàn)估算的結(jié)果非常接近，其相對誤差在允許范圍之內(nèi)。數(shù)值仿真結(jié)果之所以略大于公式推導(dǎo)結(jié)果，是因?yàn)樵诠酵茖?dǎo)中，未考慮軌道攝動力的影響。

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

2.2 爆炸碎片速度增量分析

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

式中：Δvx、Δvy、Δvz為軌道坐標(biāo)系O-xyz各軸上的速度增量；Δa、Δe、Δω、ΔΩ、Δi、ΔM為軌道根數(shù)的變化，其中，a為半長軸，e為偏心率，ω為近地點(diǎn)幅角，Ω為升交點(diǎn)赤經(jīng)，i為軌道傾角，M為平近地點(diǎn)角；μ=3.986 005×1014m3/s2，為地球引力常數(shù)；r=(x2+y2+z2)1/2，為衛(wèi)星到地心的距離；u=1/r；p為近地點(diǎn)。由上式可知，爆炸碎片的速度增量Δvx和Δvy決定了平近地點(diǎn)角、偏心率和半長軸的變化，而Δvz則決定了升交點(diǎn)赤經(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三個方向上的速度增量，結(jié)果如圖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單元的運(yùn)行軌跡并非是一條線，而是一個區(qū)域，這是因?yàn)殡S著運(yùn)行時間的增加，原來的計算模型的誤差將會不斷積累，軌道會發(fā)生很大變化。

圖7 B單元運(yùn)行軌道Fig.7 Operating orbit of unit B

3 結(jié)束語

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

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