威脅小天體動(dòng)能撞擊防御誤差演化分析*

張高輦，李翔宇

（北京理工大學(xué)宇航學(xué)院，北京 100081）

1 引言

近地小天體的軌道與地球軌道接近，部分軌道甚至與地球相交，存在潛在的撞擊風(fēng)險(xiǎn)，是人類面臨的潛在重大威脅之一。地球上曾經(jīng)發(fā)生的22 次不同程度的生物滅絕事件至少有10 次以上是由小行星撞擊地球引起的［1］，如2013 年俄羅斯車?yán)镅刨e斯克地區(qū)隕石墜落事件等。因此研究小行星撞擊防御問(wèn)題具有重要意義。近年來(lái)，針對(duì)小行星防御的研究逐漸成為航天領(lǐng)域研究的熱點(diǎn)。目前現(xiàn)有提出的方法主要包括核爆［2］、動(dòng)能撞擊［3］、激光燒蝕［4］、引力拖曳［5-7］等。其中，動(dòng)能撞擊被認(rèn)為是針對(duì)直徑1 km 以下小行星防御的高效并且可行的方法［2］。2022 年，美國(guó)國(guó)家航空航天局（NASA）成功實(shí)施針對(duì)Didymos雙小行星系統(tǒng)第二顆小行星（65803）的撞擊實(shí)驗(yàn)DART，該實(shí)驗(yàn)將是采用動(dòng)能撞擊進(jìn)行小行星防御的第一次驗(yàn)證，撞擊結(jié)果將通過(guò)歐洲航天局（ESA）未來(lái)發(fā)射的LICIA 立方星觀測(cè)得到［8］。對(duì)于威脅小行星動(dòng)能撞擊防御任務(wù)而言，動(dòng)能撞擊效果的主要評(píng)價(jià)指標(biāo)是撞擊后小行星近地點(diǎn)距離的增加量，而這一指標(biāo)可近似等效于目標(biāo)小行星在接近地球的B 平面上的偏轉(zhuǎn)距離［9，10］。以此為基礎(chǔ)Izzo［11］提出了撞擊幾何（IG，Impact Geometry）的概念用于評(píng)估撞擊偏轉(zhuǎn)距離，2019年，Davide Farnocchia 等［12］針對(duì)考慮誤差下的行星接近分析的B 平面理論做了進(jìn)一步研究。但是以往的動(dòng)能撞擊研究大多針對(duì)對(duì)心碰撞，動(dòng)量傳遞效果僅通過(guò)一個(gè)系數(shù)K 決定［13］，或者僅考慮撞擊后沿撞擊點(diǎn)處濺射物對(duì)法向方向動(dòng)量傳遞的影響，而忽略了濺射物切向方向動(dòng)量的傳遞。清華大學(xué)的焦藝菲等［14］提出一種利用兩個(gè)動(dòng)量傳遞系數(shù)β和γ描述動(dòng)能撞擊的斜撞擊模型，該模型結(jié)合撞擊幾何可實(shí)現(xiàn)最優(yōu)斜撞擊方向的快速求解，并且可應(yīng)用于任意形狀小行星最優(yōu)撞擊方向的確定。此外，中科院的王藝睿等［15］針對(duì)Apophis 偏轉(zhuǎn)任務(wù)分別研究了考慮運(yùn)載能力限制下進(jìn)行撞擊器軌道的優(yōu)化設(shè)計(jì)，后續(xù)并在考慮地球借力和多脈沖機(jī)動(dòng)模式下進(jìn)行了改進(jìn)［16］。

由于動(dòng)能撞擊一般在深空進(jìn)行，而深空環(huán)境下撞擊器的導(dǎo)航制導(dǎo)控制誤差必然會(huì)帶來(lái)動(dòng)能撞擊任務(wù)的偏差，并且小行星動(dòng)能撞擊防御的周期一般在10 年以上，動(dòng)能撞擊時(shí)刻撞擊的微小誤差可能會(huì)在時(shí)間的累積下產(chǎn)生巨大的影響，因此有必要針對(duì)威脅小天體動(dòng)能撞擊防御進(jìn)行撞擊后的誤差演化分析，進(jìn)而評(píng)估撞擊任務(wù)設(shè)計(jì)的合理性與有效性。2010 年，Armellin R 等［17］利用微分代數(shù)方法分析了小行星Apophis 在2029 年飛掠地球時(shí)的近心點(diǎn)距離和近地點(diǎn)時(shí)刻偏差；2017年，F(xiàn)eldhacker 等［18］研究了動(dòng)能撞擊位置偏差和動(dòng)量傳遞系數(shù)偏差對(duì)不同形狀小行星撞擊效果的影響，但是該文獻(xiàn)沒(méi)有分析動(dòng)能撞擊后偏差對(duì)小行星軌道分布的影響。狀態(tài)轉(zhuǎn)移張量（STT）方法是2006 年首先由Scheeres 等［19］提出的一種計(jì)算航天器誤差演化的工具，該方法通過(guò)半解析地計(jì)算狀態(tài)轉(zhuǎn)移張量實(shí)現(xiàn)了誤差預(yù)報(bào)的快速高效分析。楊震等［20，21］針對(duì)帶有脈沖機(jī)動(dòng)過(guò)程的誤差演化方法進(jìn)行了研究，推導(dǎo)出可以連續(xù)預(yù)報(bào)的狀態(tài)轉(zhuǎn)移張量。

本文針對(duì)威脅小天體動(dòng)能撞擊防御過(guò)程的誤差演化問(wèn)題，考慮由撞擊器制導(dǎo)控制誤差帶來(lái)的撞擊位置誤差，根據(jù)斜撞擊模型下的動(dòng)量傳遞關(guān)系，計(jì)算出由撞擊位置誤差引起的小行星速度增量偏差，以及偏差的均值和方差的傳遞關(guān)系，然后，基于狀態(tài)轉(zhuǎn)移張量STT 方法，針對(duì)動(dòng)能撞擊帶來(lái)的速度增量偏差演化到小行星抵達(dá)近地點(diǎn)時(shí)刻，得到小行星考慮誤差情況下在地球附近的位置速度散布情況，結(jié)合蒙特卡洛方法進(jìn)行仿真驗(yàn)證，結(jié)果表明STT 方法可以有效預(yù)報(bào)動(dòng)能撞擊過(guò)程的誤差。

2 基于狀態(tài)轉(zhuǎn)移張量的小行星誤差演化方法

2.1 動(dòng)力學(xué)方程

下面給出航天器撞擊小行星后的動(dòng)力學(xué)模型，本文選擇考慮八大行星和月球作為第三體引力攝動(dòng)影響下的動(dòng)力學(xué)方程，該方程為：

式中：dpi表示第i個(gè)行星的日心矢量，ppi表示第i個(gè)行星到小行星的矢量，μS和μpi分別表示太陽(yáng)和第i個(gè)行星的引力常數(shù)，行星和月球的位置調(diào)用SPICE星歷［22］。

2.2 狀態(tài)轉(zhuǎn)移張量的推導(dǎo)

將航天器的動(dòng)力學(xué)寫為張量形式：

式中：x={xi|i=1， …， n}，n=6 是動(dòng)力學(xué)方程的維數(shù)。對(duì)于給定的初始狀態(tài)x0=x（t0），因此x（t）可以看作初始狀態(tài)x0的函數(shù)，進(jìn)而定義隱式解為：

注意到，該隱式解滿足：

可知，若獲得了航天器運(yùn)動(dòng)狀態(tài)的非線性映射關(guān)系φ，則從航天器初始狀態(tài)x0到任意時(shí)刻狀態(tài)x（t）的傳播關(guān)系就可以解析計(jì)算。給定航天器初始參考狀態(tài)及相對(duì)參考狀態(tài)的狀態(tài)偏差δx0，可將任意時(shí)刻航天器運(yùn)動(dòng)軌跡相對(duì)標(biāo)稱軌跡的偏差量δx（t）表示為

對(duì)其求導(dǎo)可得

則標(biāo)稱軌道與當(dāng)前狀態(tài)的偏差可以表示為：

式中：

同理，對(duì)狀態(tài)量的導(dǎo)數(shù)進(jìn)行Taylor展開(kāi)可得

式中：

式中：i表示向量函數(shù)的第i個(gè)分量，M為泰勒級(jí)數(shù)展開(kāi)的階次，為局部動(dòng)力學(xué)張量，沿標(biāo)稱狀態(tài)計(jì)算，是從t0到t時(shí)刻的狀態(tài)轉(zhuǎn)移張量，這里運(yùn)用了愛(ài)因斯坦求和約定。將絕對(duì)軌道偏差的一個(gè)具體樣本δx（t）看作航天器實(shí)際狀態(tài)x（t）相對(duì)標(biāo)稱狀態(tài)x(t)的一個(gè)相對(duì)運(yùn)動(dòng)狀態(tài)，式是一組非線性的相對(duì)運(yùn)動(dòng)方程，對(duì)公式求導(dǎo)可得

將公式展開(kāi)并對(duì)比其與公式相同項(xiàng)的系數(shù)，可以得到初始相對(duì)狀態(tài)到任意時(shí)刻相對(duì)狀態(tài)的狀態(tài)轉(zhuǎn)移張量所需要滿足的微分方程。狀態(tài)轉(zhuǎn)移張量可以通過(guò)積分微分方程得到，這樣δx（t）就可以由x0解析表示，這種方法稱為半解析方法。一般地，狀態(tài)轉(zhuǎn)移張量展開(kāi)的次數(shù)越高，預(yù)報(bào)的精度越高，但是隨之帶來(lái)計(jì)算量將指數(shù)增加。文獻(xiàn)［20］表明，對(duì)于航天器軌道的誤差演化問(wèn)題，采用2 階模型就可以達(dá)到較高的精度。本文針對(duì)小行星狀態(tài)的誤差演化選取2 階精度。

2.3 誤差統(tǒng)計(jì)矩的分析

根據(jù)概率論的基本理論，自由向量δx的均值和協(xié)方差矩陣定義為：

其中E［·］代表期望算符。將式代入方程可以得出相應(yīng)的均值和協(xié)方差演化結(jié)果，即為：

由公式可知，當(dāng)各階狀態(tài)轉(zhuǎn)移張量沿標(biāo)稱軌跡計(jì)算后，基于狀態(tài)轉(zhuǎn)移張量的各階統(tǒng)計(jì)矩預(yù)報(bào)僅為解析的數(shù)學(xué)運(yùn)算。當(dāng)狀態(tài)轉(zhuǎn)移張量展開(kāi)到M階，預(yù)報(bào)均值需要計(jì)算初始偏差的前M階中心距，預(yù)報(bào)協(xié)方差矩陣需要計(jì)算初始偏差的前2M階中心距。因此，Taylor展開(kāi)階次越高、預(yù)報(bào)的統(tǒng)計(jì)矩階次越高，計(jì)算公式越復(fù)雜，且計(jì)算量越大。

本文僅考慮前2 階矩的非線性傳播，預(yù)報(bào)協(xié)方差矩陣需要計(jì)算前4階中心距，終端均值m(t0)與協(xié)方差矩陣P(t0)可簡(jiǎn)化為：

對(duì)于高斯分布，前4階統(tǒng)計(jì)矩可以寫為均值m和協(xié)方差P的表達(dá)式：

為了簡(jiǎn)便起見(jiàn)，均值和協(xié)方差誤差演化結(jié)果可以簡(jiǎn)化為：

需要注意的是，在本文研究的問(wèn)題中，目標(biāo)威脅小行星在飛掠地球時(shí)可能處于地球影響球內(nèi)部，此時(shí)小行星受力主要為地球引力，而影響球外的大部分軌道主要受太陽(yáng)引力影響。因此在使用變步長(zhǎng)積分器計(jì)算狀態(tài)轉(zhuǎn)移張量STT 時(shí)，需要注意在影響球內(nèi)降低積分步長(zhǎng)。

3 基于斜撞擊的動(dòng)能撞擊過(guò)程分析

3.1 動(dòng)能撞擊任務(wù)設(shè)計(jì)原理

為了計(jì)算動(dòng)能撞擊對(duì)小行星防御的偏轉(zhuǎn)效果，通常需要積分計(jì)算小行星整個(gè)任務(wù)過(guò)程中的非線性軌道力學(xué)方程。文獻(xiàn)［11］表明，通過(guò)B 平面計(jì)算動(dòng)能撞擊造成小行星偏轉(zhuǎn)距離的可達(dá)范圍是有效的。B 平面是指過(guò)地心并且與小行星雙曲線軌道漸近線相垂直的平面。即該平面是與進(jìn)入雙曲超速方向U 相垂直的平面，如圖1 所示。相對(duì)速度U=|U|的大小定義如下［12］：

圖1 小行星接近地球及其B平面定義Fig.1 Asteroid approaching the Earth and its B-plane definition

為了利用B 平面來(lái)描述小行星飛越地球時(shí)的軌跡，可以建立行星中心坐標(biāo)系(ξ,η,ζ) ，如圖1所示。其中η軸沿著U 矢量方向并且ζ軸指向地球速度矢量vE在B 平面投影的反向。ξ軸滿足右手定則。

對(duì)于與地球存在碰撞可能的小行星，采用動(dòng)能撞擊方法需要最大化偏轉(zhuǎn)小行星到達(dá)地球的B平面距離。文獻(xiàn)［9］給出了近似的定量關(guān)系。

式中：K代表動(dòng)量傳遞系數(shù)，一般被認(rèn)為是給定值。aast表示小行星的半長(zhǎng)軸，vE代表地球的速度大小，θaE表示小行星雙曲線漸近線與地球速度方向夾角，μs表示太陽(yáng)的引力常數(shù)，mast表示小行星的質(zhì)量，ms/c表示航天器的質(zhì)量，Vsa表示撞擊航天器相對(duì)小行星的速度矢量，vast表示小行星的速度矢量。

通常將偏轉(zhuǎn)距離Δζ作為動(dòng)能撞擊的優(yōu)化指標(biāo)，Δζ定義為地球和漸近線之間的距離。從上式可以看出為了最大化偏轉(zhuǎn)距離|Δζ|，需要最大化內(nèi)積Vsa·vast的絕對(duì)值，該值常被稱為撞擊幾何。

3.2 斜撞擊模型下的動(dòng)量傳遞關(guān)系

當(dāng)撞擊器相對(duì)小行星的速度與小行星-撞擊器連線不重合時(shí)，稱為斜撞擊。撞擊速度vr一般表示為vsc-vast，即撞擊器相對(duì)于小行星的速度［14］，其中vsc表示航天器速度，vast表示小行星速度。

圖2 定義了小行星撞擊坐標(biāo)系Oxyz，該坐標(biāo)系可以用來(lái)描述和計(jì)算不同撞擊位置對(duì)小行星動(dòng)量變化的影響。其中k 為vr反方向的單位矢量，i為撞擊瞬時(shí)小行星角動(dòng)量方向的單位矢量，j滿足右手定則，i、j、k 分別為Ox、Oy 和Oz 軸的單位方向向量。由于撞擊過(guò)程極短，本文忽略了小行星自轉(zhuǎn)帶來(lái)的影響。圖3 定義了撞擊點(diǎn)位置的極坐標(biāo)角度參數(shù)θ和φ，其中θ為從i軸逆時(shí)針轉(zhuǎn)到撞擊位置矢量在xy平面投影所需角度，φ為撞擊點(diǎn)位置矢量與k軸夾角。其中n?表示撞擊點(diǎn)處切平面指向內(nèi)部的法向單位矢量，表示單位矢量k 與單位矢量n?作2次叉乘得到的單位矢量。根據(jù)幾何關(guān)系可以得到單位矢量n?和t?關(guān)于θ和φ的顯式表達(dá)式為

圖2 撞擊坐標(biāo)系定義Fig.2 Impact coordinate system definition

圖3 小行星撞擊坐標(biāo)系定義Fig.3 Asteroid impact coordinate system definition

式中：0≤θ≤2π，為防止斜撞擊角度過(guò)大而偏離小行星，一般限制φ∈ ［0， π /3］。

根據(jù)文獻(xiàn)［14］中相關(guān)的定義，引入法向動(dòng)量傳遞系數(shù)βf和切向動(dòng)量傳遞系數(shù)γf，用以表征斜撞擊帶來(lái)的撞擊點(diǎn)法向和切向動(dòng)量傳遞效果。此時(shí)小行星的速度改變量可以表示為：

代入公式可以得到速度改變量的表達(dá)式：

撞擊參數(shù)中，需要考慮的主要有撞擊角度θimp以及兩個(gè)動(dòng)量傳遞系數(shù)β和γ。一般動(dòng)量傳遞系數(shù)通過(guò)數(shù)值仿真軟件得到，以光滑粒子流體動(dòng)力學(xué)（SPH）方法為算法基礎(chǔ)的軟件可以作為求解手段［23］。

3.3 動(dòng)能撞擊后小行星的狀態(tài)誤差分析

由于動(dòng)能撞擊一般發(fā)生在深空，遠(yuǎn)離地球測(cè)控，動(dòng)能撞擊末制導(dǎo)段不可避免產(chǎn)生撞擊誤差，進(jìn)而使得小行星的速度改變量發(fā)生偏差。本節(jié)將分析動(dòng)能撞擊過(guò)程中小行星的狀態(tài)偏差傳遞過(guò)程，并給出動(dòng)能撞擊后小行星的狀態(tài)誤差模型。

本文主要關(guān)心動(dòng)能撞擊后小行星軌道的誤差演化情況，整個(gè)動(dòng)能撞擊后小行星的運(yùn)行軌跡流程如圖4 所示。根據(jù)動(dòng)能撞擊的實(shí)際物理過(guò)程，本文將小行星撞擊后的誤差δx+0分為2 部分處理，即動(dòng)能撞擊前小行星自身的測(cè)定軌誤差δx （t0），以及動(dòng)能撞擊位置偏差帶來(lái)的小行星速度改變量的誤差δxv0。

圖4 動(dòng)能撞擊小行星示意圖Fig.4 Diagram of kinetic impact on the asteroid

首先根據(jù)動(dòng)能撞擊的實(shí)際過(guò)程，分析動(dòng)能撞擊后小行星狀態(tài)的偏差。定義動(dòng)能撞擊時(shí)刻為t0，動(dòng)能撞擊前一時(shí)刻t-0小行星的狀態(tài)量誤差為δx0，撞擊帶來(lái)的小行星速度脈沖的誤差為δΔv0，則動(dòng)能撞擊后t+0時(shí)刻小行星的狀態(tài)量偏差可以表示為：

其中R= [ 03×3,I3×3]T，δΔv0是小行星速度增量與標(biāo)稱速度增量之間的偏差，δxv0是由動(dòng)能撞擊偏差引起的小行星狀態(tài)量的偏差，δx（ t0）是小行星在動(dòng)能撞擊前的定軌誤差。

其次考慮動(dòng)能撞擊帶來(lái)的小行星速度改變量的偏差。將由動(dòng)能撞擊帶來(lái)的小行星速度改變當(dāng)作脈沖增量，則可以記標(biāo)稱情況下小行星的速度增量為考慮小行星的真實(shí)脈沖誤差為δΔv0，那么小行星真實(shí)的速度改變量Δv0為：

對(duì)于速度改變量的偏差δΔv0，參考文獻(xiàn)［18］，本文僅考慮由航天器在垂直于相對(duì)小行星速度方向上的撞擊位置偏差帶來(lái)的小行星速度增量的偏差，如圖5 所示。航天器撞擊小行星不同位置時(shí)通過(guò)前文的斜撞擊模型求解小行星的速度增量。本文假定撞擊航天器在垂直于相對(duì)速度平面（即小行星撞擊坐標(biāo)系的x-y 平面）上撞擊位置的分布為2 維高斯分布。對(duì)于球形的小行星，通過(guò)在垂直于相對(duì)速度平面求撞擊速度增量Δv0相應(yīng)的數(shù)學(xué)期望即可得到小行星速度改變量的均值m（Δv）和方差P （Δv）。具體公式如下：

圖5 垂直相對(duì)速度方向撞擊位置誤差Fig.5 Impact position error perpendicular to relative velocity direction

其中，S表示球形小行星在x-y平面上的投影區(qū)域。注意到，當(dāng)假設(shè)小行星為球形時(shí)存在幾何關(guān)系：

其中（x，y）表示撞擊位置在垂直相對(duì)速度平面上的坐標(biāo)。小行星速度改變量Δv 的方差P （Δv）表達(dá)式為：

撞擊位置的二維高斯分布概率密度函數(shù)為：

式中：x = [x, y]T，表示撞擊位置在小行星撞擊坐標(biāo)系x-y 平面投影的坐標(biāo)，并且其中均值m=m（x），方差P=P （x）。此時(shí)，本文給出了動(dòng)能撞擊后小行星位置速度的誤差模型，后文將基于該模型進(jìn)行數(shù)值仿真。

4 仿真校驗(yàn)

4.1 算例配置

假定目標(biāo)小行星為2004年發(fā)現(xiàn)的威脅小行星Apophis，該小行星質(zhì)量估計(jì)為4.0×1010kg，忽略小行星原有形狀，假定其為等質(zhì)量的均質(zhì)圓球，假定小行星的密度為2600 kg/m3，可以計(jì)算得到小行星的等效半徑約為154.288 m。Apophis小行星將于2029年4月13日飛掠地球，近地點(diǎn)距離將低于同步軌道衛(wèi)星距離?？紤]撞擊器的質(zhì)量為4500 kg，撞擊處的相對(duì)速度大小為6 km/s，相對(duì)速度方向沿小行星速度方向，動(dòng)量傳遞系數(shù)β和γ分別為2.5和0.5。假定動(dòng)能撞擊時(shí)刻為2024年4月13日21：46：07.59（TDB），關(guān)于小行星和撞擊器狀態(tài)量的設(shè)置詳見(jiàn)表1，該狀態(tài)是通過(guò)積分文獻(xiàn)［16］中2019 年1 月1 日00：00：00（TDB）小行星Apophis的位置速度到動(dòng)能撞擊時(shí)刻得到。

表1 動(dòng)能撞擊時(shí)小行星與航天器日心赤道慣性系下?tīng)顟B(tài)量Table 1 States of the asteroid and spacecraft in the heliocentric equatorial inertial system at kinetic impact

仿真結(jié)束條件為小行星到達(dá)近地點(diǎn)，本文假定地球引力影響球的半徑為925000 km。

在航天器與小行星動(dòng)能撞擊時(shí)刻，設(shè)置標(biāo)稱撞擊點(diǎn)在小行星撞擊坐標(biāo)系下的角度θ=0°，角度φ=0°，從而得到日心赤道慣性系下小行星速度標(biāo)稱改變量為：［0.10678；1.57778；0.58899］mm/s。撞擊前后小行星的軌道以及對(duì)應(yīng)時(shí)刻地球軌道關(guān)系如圖6、圖7所示。

圖6 初始小行星Apophis、撞擊后小行星Apophis軌道和地球軌道Fig.6 Initial Apophis orbit， post-impact Apophis orbit and Earth orbit

圖7 小行星Apophis到達(dá)地球影響球邊界處的位置散布Fig.7 Distribution of the positions where Apophis reaching the boundary of the Earth influence sphere

4.2 仿真結(jié)果

針對(duì)小行星在與航天器動(dòng)能撞擊后的誤差演化問(wèn)題，本文分析了小行星在僅考慮撞擊位置偏差帶來(lái)的速度增量誤差情況下，到達(dá)地球影響球邊界t1時(shí)刻以及到達(dá)近地點(diǎn)t2時(shí)刻小行星位置速度的誤差分布情況。仿真結(jié)果如下。

4.2.1 僅考慮撞擊后小行星的速度增量偏差

對(duì)于僅考慮動(dòng)能撞擊造成小行星速度增量偏差的情況，如前所述，本文假定動(dòng)量傳遞系數(shù)和撞擊相對(duì)速度無(wú)偏差，僅考慮航天器撞擊小行星表面不同位置帶來(lái)的小行星速度增量偏差。設(shè)置撞擊器在x-y 平面投影位置（x， y）的均值為mxy=[0, 0]T，方差為Pxy=diag（[5m, 5m]2）。

小行星Apophis 到達(dá)影響球邊界時(shí)在日心赤道慣性系下的位置速度散布如圖8 和圖9所示。

圖8 小行星到達(dá)影響球邊界處的位置3σ橢球Fig.8 The 3σ ellipsoid of the positions where the asteroid reaching the boundary of the influence sphere

從上述圖中可以看出，在本算例中僅考慮撞擊位置偏差，演化5 年后，當(dāng)小行星到達(dá)地球影響球邊界時(shí)，其位置散布在10 km 量級(jí)，速度散布在0.1 cm/s量級(jí)。

如圖10 和圖11 所示，可以看出動(dòng)能撞擊前后小行星軌道在近心點(diǎn)附近的改變，從圖11可以清晰看出動(dòng)能撞擊后小行星位置誤差散布。對(duì)比圖11和圖7可以看出小行星在到達(dá)影響球邊界和到達(dá)近地點(diǎn)時(shí)的位置偏差3σ散布差異不大，均在10 km左右。

圖11 小行星到達(dá)近地點(diǎn)處的3σ橢球Fig.11 The 3σ ellipsoid at which the asteroid reaching its perigee

圖12 小行星到達(dá)近地點(diǎn)處的位置3σ橢球Fig.12 The 3σ ellipsoid of positions where the asteroid reaching its perigee

從上述圖中可以看出，在本算例中僅考慮撞擊位置偏差，當(dāng)小行星到達(dá)近地點(diǎn)時(shí)，其位置散布在10 km量級(jí)，速度散布在10 cm/s量級(jí)。但是對(duì)比圖13和圖9可以看出，小行星在到達(dá)近地點(diǎn)時(shí)的速度誤差散布較在地球影響球邊界處高2 個(gè)數(shù)量級(jí)。同時(shí)，上述結(jié)果表明，采用2 階STT 可以對(duì)小行星的終端位置速度誤差進(jìn)行有效預(yù)報(bào)，預(yù)報(bào)結(jié)果與蒙特卡洛仿真結(jié)果吻合。

圖13 小行星到達(dá)近地點(diǎn)處的速度3σ橢球Fig.13 The 3σ ellipsoid of velocities when the asteroid reaching its perigee

5 結(jié)論

本文基于STT 方法針對(duì)小行星動(dòng)能撞擊防御中的斜撞擊進(jìn)行誤差演化分析。首先，基于斜撞擊模型描述航天器撞擊小行星過(guò)程，通過(guò)兩個(gè)動(dòng)量傳遞參數(shù)β和γ得到撞擊后小行星的速度改變量；然后，針對(duì)考慮航天器撞擊小行星表面的位置偏差，在斜撞擊的模型基礎(chǔ)上推導(dǎo)了撞擊位置誤差到撞擊后速度增量誤差之間的均值和方差傳遞關(guān)系；最后，針對(duì)撞擊后小行星位置速度的偏差結(jié)果，通過(guò)狀態(tài)轉(zhuǎn)移張量STT 方法演化得到小行星進(jìn)入地球影響球邊界以及近地點(diǎn)處的位置速度誤差。結(jié)果表明，動(dòng)能撞擊位置誤差對(duì)小行星對(duì)動(dòng)能撞擊效果的影響較大，對(duì)于半徑150 m 左右的小行星、小行星速度增量1～2 mm/s量級(jí)的動(dòng)能撞擊而言，3σ邊界為15 m的撞擊位置偏差將在5 年后使得小行星在近地點(diǎn)處產(chǎn)生10 公里左右的位置偏差，0.1 m/s量級(jí)的速度偏差散布。如果考慮到小行星質(zhì)量估計(jì)的偏差、小行星材質(zhì)帶來(lái)的動(dòng)量傳遞系數(shù)偏差等因素的影響，真實(shí)的動(dòng)能撞擊帶來(lái)的位置速度散布將更大。因此在動(dòng)能撞擊任務(wù)設(shè)計(jì)的過(guò)程中需要考慮撞擊航天器的制導(dǎo)控制精度的影響，避免撞擊誤差對(duì)撞擊效果的不利影響。