威脅小天體動能撞擊防御誤差演化分析*

張高輦，李翔宇

（北京理工大學宇航學院，北京 100081）

1 引言

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

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

本文針對威脅小天體動能撞擊防御過程的誤差演化問題，考慮由撞擊器制導控制誤差帶來的撞擊位置誤差，根據斜撞擊模型下的動量傳遞關系，計算出由撞擊位置誤差引起的小行星速度增量偏差，以及偏差的均值和方差的傳遞關系，然后，基于狀態轉移張量STT 方法，針對動能撞擊帶來的速度增量偏差演化到小行星抵達近地點時刻，得到小行星考慮誤差情況下在地球附近的位置速度散布情況，結合蒙特卡洛方法進行仿真驗證，結果表明STT 方法可以有效預報動能撞擊過程的誤差。

2 基于狀態轉移張量的小行星誤差演化方法

2.1 動力學方程

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

式中：dpi表示第i個行星的日心矢量，ppi表示第i個行星到小行星的矢量，μS和μpi分別表示太陽和第i個行星的引力常數，行星和月球的位置調用SPICE星歷［22］。

2.2 狀態轉移張量的推導

將航天器的動力學寫為張量形式：

式中：x={xi|i=1， …， n}，n=6 是動力學方程的維數。對于給定的初始狀態x0=x（t0），因此x（t）可以看作初始狀態x0的函數，進而定義隱式解為：

注意到，該隱式解滿足：

可知，若獲得了航天器運動狀態的非線性映射關系φ，則從航天器初始狀態x0到任意時刻狀態x（t）的傳播關系就可以解析計算。給定航天器初始參考狀態及相對參考狀態的狀態偏差δx0，可將任意時刻航天器運動軌跡相對標稱軌跡的偏差量δx（t）表示為

對其求導可得

則標稱軌道與當前狀態的偏差可以表示為：

式中：

同理，對狀態量的導數進行Taylor展開可得

式中：

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

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

2.3 誤差統計矩的分析

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

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

由公式可知，當各階狀態轉移張量沿標稱軌跡計算后，基于狀態轉移張量的各階統計矩預報僅為解析的數學運算。當狀態轉移張量展開到M階，預報均值需要計算初始偏差的前M階中心距，預報協方差矩陣需要計算初始偏差的前2M階中心距。因此，Taylor展開階次越高、預報的統計矩階次越高，計算公式越復雜，且計算量越大。

本文僅考慮前2 階矩的非線性傳播，預報協方差矩陣需要計算前4階中心距，終端均值m(t0)與協方差矩陣P(t0)可簡化為：

對于高斯分布，前4階統計矩可以寫為均值m和協方差P的表達式：

為了簡便起見，均值和協方差誤差演化結果可以簡化為：

需要注意的是，在本文研究的問題中，目標威脅小行星在飛掠地球時可能處于地球影響球內部，此時小行星受力主要為地球引力，而影響球外的大部分軌道主要受太陽引力影響。因此在使用變步長積分器計算狀態轉移張量STT 時，需要注意在影響球內降低積分步長。

3 基于斜撞擊的動能撞擊過程分析

3.1 動能撞擊任務設計原理

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

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

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

對于與地球存在碰撞可能的小行星，采用動能撞擊方法需要最大化偏轉小行星到達地球的B平面距離。文獻［9］給出了近似的定量關系。

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

通常將偏轉距離Δζ作為動能撞擊的優化指標，Δζ定義為地球和漸近線之間的距離。從上式可以看出為了最大化偏轉距離|Δζ|，需要最大化內積Vsa·vast的絕對值，該值常被稱為撞擊幾何。

3.2 斜撞擊模型下的動量傳遞關系

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

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

圖2 撞擊坐標系定義Fig.2 Impact coordinate system definition

圖3 小行星撞擊坐標系定義Fig.3 Asteroid impact coordinate system definition

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

根據文獻［14］中相關的定義，引入法向動量傳遞系數βf和切向動量傳遞系數γf，用以表征斜撞擊帶來的撞擊點法向和切向動量傳遞效果。此時小行星的速度改變量可以表示為：

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

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

3.3 動能撞擊后小行星的狀態誤差分析

由于動能撞擊一般發生在深空，遠離地球測控，動能撞擊末制導段不可避免產生撞擊誤差，進而使得小行星的速度改變量發生偏差。本節將分析動能撞擊過程中小行星的狀態偏差傳遞過程，并給出動能撞擊后小行星的狀態誤差模型。

本文主要關心動能撞擊后小行星軌道的誤差演化情況，整個動能撞擊后小行星的運行軌跡流程如圖4 所示。根據動能撞擊的實際物理過程，本文將小行星撞擊后的誤差δx+0分為2 部分處理，即動能撞擊前小行星自身的測定軌誤差δx （t0），以及動能撞擊位置偏差帶來的小行星速度改變量的誤差δxv0。

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

首先根據動能撞擊的實際過程，分析動能撞擊后小行星狀態的偏差。定義動能撞擊時刻為t0，動能撞擊前一時刻t-0小行星的狀態量誤差為δx0，撞擊帶來的小行星速度脈沖的誤差為δΔv0，則動能撞擊后t+0時刻小行星的狀態量偏差可以表示為：

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

其次考慮動能撞擊帶來的小行星速度改變量的偏差。將由動能撞擊帶來的小行星速度改變當作脈沖增量，則可以記標稱情況下小行星的速度增量為考慮小行星的真實脈沖誤差為δΔv0，那么小行星真實的速度改變量Δv0為：

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

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

其中，S表示球形小行星在x-y平面上的投影區域。注意到，當假設小行星為球形時存在幾何關系：

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

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

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

4 仿真校驗

4.1 算例配置

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

表1 動能撞擊時小行星與航天器日心赤道慣性系下狀態量Table 1 States of the asteroid and spacecraft in the heliocentric equatorial inertial system at kinetic impact

仿真結束條件為小行星到達近地點，本文假定地球引力影響球的半徑為925000 km。

在航天器與小行星動能撞擊時刻，設置標稱撞擊點在小行星撞擊坐標系下的角度θ=0°，角度φ=0°，從而得到日心赤道慣性系下小行星速度標稱改變量為：［0.10678；1.57778；0.58899］mm/s。撞擊前后小行星的軌道以及對應時刻地球軌道關系如圖6、圖7所示。

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

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

4.2 仿真結果

針對小行星在與航天器動能撞擊后的誤差演化問題，本文分析了小行星在僅考慮撞擊位置偏差帶來的速度增量誤差情況下，到達地球影響球邊界t1時刻以及到達近地點t2時刻小行星位置速度的誤差分布情況。仿真結果如下。

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

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

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

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

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

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

圖11 小行星到達近地點處的3σ橢球Fig.11 The 3σ ellipsoid at which the asteroid reaching its perigee

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

從上述圖中可以看出，在本算例中僅考慮撞擊位置偏差，當小行星到達近地點時，其位置散布在10 km量級，速度散布在10 cm/s量級。但是對比圖13和圖9可以看出，小行星在到達近地點時的速度誤差散布較在地球影響球邊界處高2 個數量級。同時，上述結果表明，采用2 階STT 可以對小行星的終端位置速度誤差進行有效預報，預報結果與蒙特卡洛仿真結果吻合。

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

5 結論

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