潛在威脅小行星碰撞防御的計算與分析

姜宇，程彬，寶音賀西，李恒年

（1. 西安衛(wèi)星測控中心宇航動力學國家重點實驗室，西安 710043；2. 清華大學 航天航空學院，北京 100084）

潛在威脅小行星碰撞防御的計算與分析

姜宇1,2*，程彬2，寶音賀西2，李恒年1

（1. 西安衛(wèi)星測控中心宇航動力學國家重點實驗室，西安 710043；2. 清華大學 航天航空學院，北京 100084）

研究了采用碰撞的方式進行小行星防御的動力學問題。采用多面體模型來建立小行星的外形模型，以碎石堆模型來建立小行星的結構模型，計算了小行星受到與其密度和材質相同的球體高速碰撞過程和碰撞后的碎石分布。計算過程在考慮了小行星與碰撞球體的接觸形變以及小行星內部組成碎石堆的接觸形變條件下，計算了碎石堆內部的相互引力、法向接觸力、切向靜摩擦力、切向動摩擦力和滾動摩擦力矩。以小行星101955 Bennu（中文名貝努）為對象計算了潛在威脅小行星的碰撞防御過程的動力學行為。結果顯示：采用高速碰撞的方法進行小行星防御可以有效地將小行星撞成大量碎小的石塊，且該方法具有核爆的方法不可比擬的優(yōu)勢，即對空間環(huán)境無污染。

潛在威脅小行星；小行星防御；小行星碰撞；碰撞防御；碎石堆小行星

0 引 言

有的小行星屬于潛在威脅[1-6]的小行星，這類小行星可能與地球相撞。對于已確定將要撞擊地球的小行星[3]，可以采用對小行星進行變軌控制[1]、對小行星進行碰撞[4-6]或使用核武器將其引爆[7]等方式將其炸毀。Quarta和Mengali[1]研究了采用電動太陽帆進行潛在威脅小行星的交會任務，并與采用平坦的太陽帆進行潛在威脅小行星的交會進行了比較。Vardaxis等[4]研究了針對潛在威脅小行星進行碰撞的從地球至小行星的轉移軌跡設計。Michel等[5]研究了采用沖擊動能技術偏轉雙小行星系統(tǒng)65803 Didymos過程中的動力學參數，該研究的背景是小行星碰撞與偏轉評估（Asteroid Impact & Deflection Assessment，AIDA）任務[5-6]，AIDA是第一個針對雙小行星系統(tǒng)表面及內部進行觀測與探索的任務。Lomov等[7]研究了采用核爆的方式摧毀小行星來防止?jié)撛谕{小行星碰撞地球，研究表明該方式可以提供足夠的能量來摧毀小行星內部結構的完整性。此外，他們發(fā)現(xiàn)無孔結構的固體巖石狀小行星在受到核爆時的損傷程度比含大量孔隙結構的小行星在受到核爆后損傷程度大，因為無孔結構無法像多孔介質那樣衰減核爆形成的沖擊波，盡管多孔介質小行星的結構強度不如前者。Koenig和Chyba[8]從動量定理出發(fā)研究了采用動能沖擊來偏轉小行星軌道的方法，考慮了小行星的質量和相對速度。

本文以潛在威脅的小行星101955 Bennu（中文名貝努）[2,9-10]為例，研究了采用碰撞來摧毀小行星的方法并計算了碰撞過程的相關參數。美國于2016年9月8日發(fā)射的OSIRIS-REx任務探測器預計于2018年抵達小行星101955 Bennu，采樣后將于2023年返回地球[2]。2013年俄羅斯車里雅賓斯克州隕石事件的隕石在進入大氣層以前的直徑約為15 m，造成了超過700人受傷和大量建筑物受損[11-12]，而小行星101955 Bennu的直徑約為492 m[2]，如果撞擊地球，災難將是毀滅性的。

本文研究通過高速碰撞將小行星撞擊成大量碎片。對于將要撞擊地球的小行星來說，在其撞擊地球以前，人類將該小行星撞擊成大量碎片，進入大氣層的部分碎片會燒毀，不至于造成嚴重災難。研究表明，尺度在小行星五分之一大小量級的球體高速碰撞小行星可以把小行星完全撞碎。與采用核爆的方式撞碎小行星相比，這種方法無輻射污染，在保護環(huán)境方面具有明顯的優(yōu)勢。

1 小行星101955 Bennu的形狀與物理特征

使用多面體模型來建立小行星的外形與引力場模型，小行星101955 Bennu的引力勢[13]

其中：G = 6.67 × 10–11m3kg–1s–2表示牛頓引力常數，σ表示小行星101955 Bennu的體密度，re和rf分別為小行星本體固連坐標系中表示的從場點分別到邊e和面f的向量，Ee和Ff分別為多面體模型的邊和面的參數，Le為多面體模型中場點和邊的積分因子，ωf為相對多面體模型中場點的帶符號的角度。小行星101955 Bennu[14-15]的自旋周期為4.288 h，平均直徑為492 ± 20 m，體密度為0.95 g/cm3。

圖 1給出了小行星101955 Bennu的有效勢[9-10,14-16]的整體結構情況和在各平面的投影結構。圖 2給出了小行星101955 Bennu不規(guī)則幾何形狀及平衡點的位置關系以及平衡點和有效勢能分布的位置關系。計算過程采用的小行星101955 Bennu物理模型含有1 348個頂點和2 692個面。由圖 1（b）和圖 1（c）可見，小行星101955 Bennu的有效勢在xy平面的分布和在yz平面及zx平面的分布完全不同，有效勢在yz平面及zx平面的結構基本類似。圖 1（a）給出了小行星101955 Bennu的有效勢在xy平面的大小數值，由此可見，小行星101955 Bennu的xy平面有效勢存在臨界點。對比圖 1（a）和圖 2（b）、圖 2（e）可知該小行星有效勢和平衡點的位置關系。小行星101955 Bennu共有9個平衡點，其中1個位于小行星體內，8個位于小行星體外。從穩(wěn)定性和拓撲類型的角度來看，位于小行星體內的平衡點是穩(wěn)定的，另外位于小行星體外的8個平衡點都是不穩(wěn)定的。由圖 2（c）和圖 2（d）可見，小行星101955 Bennu的平衡點并未落在一個平面內。事實上，這些平衡點都是非赤道面的平衡點，都落在小行星101955 Bennu自旋的赤道面以外。

由圖 1和圖 2可見，小行星101955 Bennu的形狀與物理特征較為復雜，外部有8個相對平衡點，也就是不規(guī)則外形與引力分布引起的有效勢能臨界點有8個。臨界點偏離赤道面，說明該小行星南北并非對稱。因此，在研究小行星101955 Bennu作為潛在威脅的小行星進行碰撞時，不能簡單地將該小行星作為圓球或橢球來考慮，否則將不能真實地模擬該小行星的不規(guī)則幾何外形及其引起的引力場環(huán)境。此外，雖然101955 Bennu的物理模型與外形數據是通過多面體模型給出的，包含數千個頂點和面，但多面體模型作為一個整體，僅在計算小行星的幾何外形和引力場乃至計算其引力場中的動力學行為方面存在諸多優(yōu)勢，而無法計算小行星受到碰撞發(fā)生斷裂或碎石堆小行星解體成大量碎石的過程。因此我們需要尋找一種方法，既能良好地保留多面體模型關于小行星不規(guī)則外形的模擬優(yōu)點，能較為精確地還原小行星的真實外形，又能良好地模擬碎石堆小行星的內部結構。

圖 1 小行星101955 Bennu的有效勢Fig. 1 Effective potential of the asteroid 101955 Bennu

圖 2 小行星101955 Bennu的形狀及平衡點Fig. 2 The shape and equilibrium points of the asteroid 101955 Bennu

2 小行星101955 Bennu碰撞的數值計算

根據已有的小行星近距離探測及熱紅外成像數據，我們發(fā)現(xiàn)大部分直徑在100 m～100 km之間的小行星（包括101955 Bennu）是由直徑幾毫米至幾十米的碎石顆粒在引力作用下聚合形成的，即“碎石堆”結構。這些顆粒之間具有接觸力和引力的作用，相鄰顆粒之間存在碰撞、摩擦和滾動。在力的作用下，相鄰顆粒之間并非剛性接觸，而是存在著接觸形變。采用軟球離散元耦合引力N體模型可以很好地模擬這樣的力學行為[17]。這里我們采用軟球離散元方法來建立小行星的碎石堆模型。通過足夠多的球來填充小行星的通過多面體模型數據建立的幾何外形的內部，從而得到采用碎石堆模型表示的小行星的內部結構和外形。在小行星碎石堆模型中，相鄰碎石顆粒之間考慮法向的彈性力和阻尼力，而在切向則考慮彈性力、阻尼力、滑動力、滾動力及合力矩。

第p個小球受到的合力Fp與合力矩Mp采用下式計算[18]

其中：kn為法向接觸剛度；kt為切向接觸剛度；Cn為法向阻尼系數；CT為切向阻尼系數；FN為法向接觸力；S為相對于平衡接觸點的切向位移；為S的單位矢量；為相鄰顆粒間的法向單位矢量；為切向單位矢量；μs為靜摩擦系數；μr為滾動摩擦系數；un為總的相對速度的法向分量；ut為總的相對速度的切向分量；A為中間變量；vrot為滾動摩擦系數計算的中間過程量；x為相鄰顆粒之間的重疊量。

法向與切向接觸剛度需保證碰撞過程中顆粒最大重疊量不超過最小顆粒半徑的1%，即xmax約等于1%Rmin，一般通過下式確定

其中：m為顆粒質量，vmax為顆粒系統(tǒng)最大碰撞速度，Rmin為最小顆粒半徑，εn為法向彈性恢復系數，μ為特征碰撞對約化質量。

以碎石堆模型來模擬小行星101955 Bennu的外形與內部結構，采用4 875個顆粒來填充多面體模型。設正向彈性恢復系數為0.5，滑動摩擦系數為0.5，滾動摩擦系數為0.55。計算的時間步長為微秒，即1 μs = 1.0 × 10–6s。計算得到正向剛度系數為2.094 395 × 1011，正向阻尼系數為4.512 453 × 108，切向剛度系數為5.983 986 × 1010，切向阻尼系數為7.931 700 × 107。

圖 3給出了小行星101955 Bennu受到碰撞過程的數值計算結果。為了便于觀察，采用碎石顆粒的顏色來表示小行星表面及內部碎石顆粒的相對速度。

圖 3 小行星101955 Bennu碰撞過程的顆粒受力分布Fig. 3 Calculation of the impact of the asteroid 101955 Bennu

由圖 3（a）可見，小行星101955 Bennu受到其五分之一直徑大小的球體碰撞初始碰撞時，小行星的星體上顆粒受力大小分布與顆粒和碰撞球體之間的距離有關，小行星表面相鄰顆粒的受力大小比較接近，同時緊鄰撞擊物附近顆粒受力遠遠高于背離撞擊區(qū)域。從圖 3（b）和圖 3（c）可見，隨著碰撞的進一步深入進行，雖然球體已經逐漸進入小行星內部，但小行星表面顆粒的受力大小分布仍然和表面顆粒與初始碰撞接觸時刻小行星表面碰撞點的距離有關，只是出現(xiàn)了一些受力大小與相鄰顆粒受力大小相差較大的顆粒，但從小行星整體來看，小行星表面顆粒速度大小呈現(xiàn)出明顯的分層特點，即初始撞擊點附近顆粒由于被拋射出小行星本體，不再受到由于撞擊而產生的沖擊波的影響，受力迅速下降；而受力最大區(qū)域（白色區(qū)域）隨著撞擊物的深入逐漸沿撞擊方向移動，最終完全貫穿小行星本體。將圖 3（d）和圖 3（a）對比可見，小行星受到碰撞后期的表面顆粒受力分布和碰撞早期的表面顆粒速度分布幾乎完全相反，同時小行星表層顆粒在碰撞過程中由表面向內部逐層被完全拋射出小行星。

圖 4給出了小行星101955 Bennu受到碰撞以后的碎石顆?？臻g瞬時分布情況。碰撞以后，坐標軸x方向的碎石顆粒最大速度分量為449.191 929 2 m/s，最小速度分量為–368.694 478 2 m/s；坐標軸y方向的碎石顆粒最大速度分量為385.1557403 m s–1，最小速度分量為–328.650 189 7 m/s；坐標軸z方向的碎石顆粒最大速度分量為384.660 715 9 m/s，最小速度分量為–348.184 552 8 m/s。碰撞后碎石顆粒的速度的模的最大值為610.867 m/s，碎石顆粒的速度的模的最小值為21.775 49 m/s。

從圖 4的碰撞以后碎石顆?？臻g瞬時分布情況來看，以小行星五分之一大小的球體高速碰撞小行星，可以把小行星完全摧毀，撞擊成大量碎石。以碰撞前小行星質心為坐標原點，最終碎石顆粒的速度大小約在數百米每秒量級，在該速度下，小行星碎片云無法在引力作用下重聚起來，因此可以認為該小行星已被完全摧毀瓦解，對地球的撞擊威脅已基本降低到安全水平以下。這里需要特別指出，在實際的撞擊過程中，碎石顆粒之間的碰撞破碎將會產生大量碎屑和粉塵，但由于計算機計算能力的限制，不能完全模擬出撞擊產生的所有碎屑和粉塵的動力學行為，只能近似地計算出撞擊后特定尺寸大小以上的碎石的動力學行為。

圖 4 小行星101955 Bennu碰撞以后的碎石顆粒分布Fig. 4 Calculation of the impact of the asteroid 101955 Bennu

Lomov等[7]在2013年的研究中采用核爆的方式來炸毀潛在威脅的小行星，對于個頭較大的小行星，例如數百千米量級的小行星來說，目前核爆的方法是摧毀這類小行星的唯一有效方法[7,19]。此外，對于事先并未發(fā)現(xiàn)的突然出現(xiàn)的將要撞擊地球的幾十千米以下量級的小行星，核爆的方法也能快速有效地將其摧毀。然而，核爆的方法摧毀潛在威脅的小行星也會產生大量碎屑和粉塵。同本文的結果不同的是，核爆的方式炸毀小行星產生的碎屑和粉塵具有強烈的放射性。這些具有強烈放射性的碎屑和粉塵以及個頭稍大的碎石會向四面八方散開，一部分可能會進入地球大氣層，對人類的生存環(huán)境造成難以估量的影響。

因此，Lomov等[7]研究的采用核爆的方式來炸毀將要撞擊地球的小行星的辦法，在針對數百千米大小的小行星或者事先未發(fā)現(xiàn)的突然出現(xiàn)的小行星方面，是萬不得已的唯一辦法。然而通常情況下，本文所采用的方法既能有效防御潛在威脅的小行星，又不對太空環(huán)境產生放射性污染，是兩全其美的辦法。

3 結 論

研究了針對潛在威脅的小行星的碰撞攔截問題。采用高精度多面體模型數據來建立小行星的形狀與物理模型。采用碎石堆填充多面體模型來建立小行星的內部結構模型。研究以有較為精確的外形數據的小行星101955 Bennu為例計算了針對小行星的碰撞過程。使用一個大小為小行星五分之一的球狀物碰撞小行星。研究表明：這種方式可以有效地將小行星撞成大量碎小的石塊，可以用來對潛在威脅的小行星進行防御。

[1]Quarta A A，Mengali G. Electric sail missions to potentially hazardous asteroids[J]. Acta Astronautica，2010，66（9）：1506-1519.

[2]Chesley S R，F(xiàn)arnocchia D，Nolan M C，et al. Orbit and bulk density of the OSIRIS-REx target Asteroid（101955） Bennu[J]. Icarus，2014，235（1）：5-22.

[3]Saito T，Kaiho K，Abe A，et al. Hypervelocity impact of asteroid/comet on the oceanic crust of the earth [J]. International Journal of Impact Engineering，2008，35：1770-1777.

[4]Vardaxis G，Sherman P，Wie B. Impact risk assessment and planetary defense mission planning for asteroid 2015 PDC[J]. Acta Astronautica，2016，122：307-323.

[5]Michel P，Cheng A，Küppers M，et al. Science case for the Asteroid Impact Mission（AIM）：a component of the Asteroid Impact & Deflection Assessment（AIDA） mission[J]. Advances in Space Research，2016，57（12）：2529-2547.

[6]Cheng F A，Atchison J，Kantsiper B，et al. Asteroid impact and deflection assessment mission [J]. Acta Astronautica，2015，115：262-269.

[7]Lomov I，Herbold E B，Antoun T H，et al. Influence of mechanical properties relevant to standoff deflection of hazardous asteroids [J]. Procedia Engineering，2013，58（4）：251-259.

[8]Koenig J D，Chyba C F. Impact deflection of potentially hazardous asteroids using current launch vehicles[J]. Science & Global Security，2007，15（1）：57-83.

[9]Jiang Y，Baoyin H X，Li H N. Periodic motion near the surface of asteroids [J]. Astrophysics and Space Science，2015，360（2）：63.

[10]Wang X Y，Jiang Y，Gong S P. Analysis of the potential field and equilibrium points of irregular-shaped minor celestial bodies[J]. Astrophysics and Space Science，2014，353（1）：105-121.

[11]Popova O P，Jenniskens P，Emelyanenko V，et al. Chelyabinsk airburst，damage assessment，meteorite recovery，and characterization [J]. Science，2013，342（6162）：1069-1073.

[12]Brumfiel G.Russian meteor largest in a century [N]. Nature News. England，Macmillan Publishers Limited，part of Springer Nature，2013. doi：10.1038/nature.12438.

[13]Werner R A，Scheeres D J. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia[J]. Celestial Mechanics and Dynamical Astronomical，1997，65（3）：313-344.

[14]Nolan M C，Magri C，Howell E S，et al. Shape model and surface properties of the OSIRIS-REx target Asteroid（101955） Bennu from radar and lightcurve observations [J]. Icarus，2013，226（1）：629-640.

[15]Nolan M C，Magri C，Howell E S，et al. Asteroid（101955） Bennu Shape Model V1.0 [Z]. NASA Planetary Data System 211. Los Angeles：University of California，2013.

[16]Jiang Y，Baoyin H X，Li J F，et al. Orbits and manifolds near the equilibrium points around a rotating asteroid [J]. Astrophysics and Space Science，2014，349（1）：83-106.

[17]Jiang Y，Zhang Y，Baoyin H X. Surface motion relative to the irregular celestial bodies [J]. Planetary and Space Science，2016，127：33-43.

[18]Schwartz S R，Richardson D C，Michel P. An implementation of the soft-sphere discrete dlement method in a high-performance parallel gravity tree-code [J]. Granular Matter，2012，14：363-380.

[19]馬鵬斌，寶音賀西. 近地小行星威脅與防御研究現(xiàn)狀[J]. 深空探測學報，2016，3（1）：10-17. Ma P B，Baoyin H X. Research status of the near-Earth asteroids' hazard and mitigation[J].Journal of Deep Space exploration，2016，3（1）：10-17.

通信地址：清華大學蒙民偉科技大樓北樓N904（100084）

E-mail：jiangyu_xian_china@163.com

Calculation and Analysis of the Impact Defense to the Potentially Hazardous Asteroids

JIANG Yu1,2*，CHEN Bin2，BAOYIN Hexi2，LI Hengnian1
（1. State Key Laboratory of Astronautic Dynamics，Xi’an Satellite Control Center, Xi’an 710043，China；2. School of Aerospace Engineering，Tsinghua University，Beijing 100084，China）

The dynamical problems of the asteroid defense are investigated by using the impact method. The polyhedron model isused to build the shape model of the asteroids, and the rubble-pile model is used to build the structure model of the asteroids. A sphere which has the same dense and material quality with the asteroid is assumed to impact to the asteroid, and the high velocity impact and the distribution of the rocks and rubbles after the impact are computed. Considering the contact deformation between the asteroid and the impact sphere, and the contact deformation inside the body of the asteroid among the rocks and rubbles which comprise the structure of the asteroid, the mutual gravitational force, the normal contact force, the tangential static frictional force, the tangential kinetic friction force, and the moment of rolling friction are calculated. Asteroid 101955 Bennu is used to calculate the dynamical behaviors of the impact defense of the potentially hazardous asteroids. The results show that the method of high velocity impact is useful for the defense of potentially hazardous asteroids, the asteroids can be crashed into a large number of rocks and rubble. The presented methodhas no pollution to the space environment. In this sense, the method is much better than the nuclear explosion method.

potentially hazardous asteroids；asteroid defense；asteroid impact；impact defense；rubble-pile asteroids

P185.7

A

2095-7777（2017）02-0190-06

10.15982/j.issn.2095-7777.2017.02.014

姜宇（1983– ），男，博士后，工程師。主要研究方向：太空目標碰撞與沖擊動力學，多小行星系統(tǒng)動力學與探測，衛(wèi)星星座與編隊控制。

[責任編輯：宋宏，英文審校：朱魯青]

姜宇，程彬，寶音賀西，等. 潛在威脅小行星碰撞防御的計算與分析[J]. 深空探測學報，2017，4（2）：190-195.

Reference format: Jiang Y，Cheng B，Baoyin H X，et al. Calculation and analysis of the impact defense to the potentially hazardous asteroids [J]. Journal of Deep Space Exploration，2017，4（2）：190-195.

2016-10-01

2017-04-09

國家自然科學基金資助項目（11372150）