基于偏轉距離近似模型的動能撞擊小行星防御任務脈沖軌道優化研究1)

王藝睿 李明濤 周炳紅

(中國科學院國家空間科學中心,北京100190)(中國科學院大學,北京 100049)

引言

太陽系內的小行星主要分布于火星軌道與木星軌道之間,也有一部分小行星會穿越地球軌道威脅到地球上的生命.諸多嚴重的小行星撞擊地球事件(如6500 萬年前的K-T 事件[1]、1908 年通古斯大爆炸事件[2]、2013 年車里雅賓斯克事件[3]等)已引發人類對于行星防御的關注.在能夠對近地小行星提前預警的前提下,摧毀小行星結構或者偏轉小行星軌道是避免其撞擊地球的兩種基本方式[4].目前已提出的防御手段主要有核爆[5]、動能撞擊、激光燒蝕[6]、拖船[7]、引力拖船[8]、離子束牽引[9]等,其中動能撞擊技術是目前最易實現且成熟度較高的手段.

動能撞擊技術是指撞擊器以一定的速度和角度撞擊小行星,使小行星軌道發生改變.2005 年美國實施了深度撞擊(deep impact)任務[10],任務釋放一顆約370 kg 的撞擊器,以約10.3 km/s 的相對速度撞擊“Tempel 1”彗星的核,撞擊后的3 年內該彗星位置變化了約10 km;雙小行星重定向測試任務(double asteroid redirection test,DART)將于2021 年7 月執行,由一個約650 kg 的撞擊器以約6.65 km/s 的相對速度撞擊小行星Didymos 的衛星,預計產生0.8 ～2 mm/s 的速度改變量(取決于動量傳遞因子的大小,動量傳遞因子用于描述濺射物對動量傳遞的影響[11]).Atchison 等[12]總結了DART 軌道設計需求,包括撞擊日期、撞擊角度以及為了滿足光學需求的太陽相位角需求;Ozimek 和Atchison[13]研究了DART 軌道中逃逸段的月球借力方法;Atchison 等[14]介紹了DART的終端制導、載荷等信息.

早期對于動能撞擊處置小行星的研究,多關注動能撞擊對小行星的軌道偏轉機理,一般不考慮運載能力約束.Ahrens 和Harris[15]的研究表明切向偏轉速度增量對偏轉小行星軌道效率最高,并估計了偏轉小行星所需的速度增量,但如果小行星在撞擊地球前會近距離飛掠大行星,實際所需的偏轉速度增量會有所降低[16];Park 和Ross[17]推導了小行星與地球距離的一、二階導數表達式,發現只有在撞擊前某些特殊時刻最優偏轉速度才平行于小行星速度方向;Ross 等[18]基于圓錐曲線拼接法考慮了地球引力對最優偏轉速度增量的影響.Carusi 等[16]提出一種估算在多體動力學模型下將小行星偏轉1 倍地球半徑所需速度增量的方法,但僅研究了切向速度增量;Kahle 等[19]對其方法進行了改進,在多體問題下同時優化了切向、徑向的速度增量,發現Carusi 等[16]的“切向速度增量”結論僅在二體階段下是合理的,當小行星近距離飛掠行星時,最優速度增量方向并不是切向的了,近距離飛掠行星時受行星引力影響可能會降低偏轉所需的速度增量.Vasile 和Colombo[20]發現當偏轉時間小于小行星公轉周期時,最優速度增量的徑向分量占主導,當偏轉時間更大時,切向分量占主導.Yamaguchi 和Yamakawa[21]提出了IGMs(impactgeometry maps),用于可視化偏轉距離與IG(impact geometry)(小行星速度與航天器撞擊速度的內積)之間的關系.Greenstreet 等[22]基于統計的方法,發現與行星發生近距離飛掠的小行星,所需的偏轉速度增量可以非常小.

動能撞擊任務設計的一種優化指標是最大化小行星被偏轉前后的近地點距離改變量(偏轉距離),換言之是使小行星被偏轉后的近地點盡可能的離地球更遠[23].采用多體型計算小行星被偏轉前后的近地點距離,可以提供關于問題最完整的解,但軌道演化和近地點搜索計算過程中會產生冗余計算量,在任務設計初期的快速搜索研究(quick scoping studies)中,多體數值積分是不必要的[24].采用二體模型替代多體模型的方法可以顯著降低軌道演化過程的計算量,但在搜索近地點的過程中仍存在冗余計算量.為同時降低軌道演化和近地點搜索的計算量,可采用解析偏轉模型對偏轉距離進行估算.

除了Ahrens 和Harris[15]提出的偏轉距離線性解析模型外,目前使用較為廣泛的的偏轉距離解析模型還有:2001 年Conway[25]提出了基于拉格朗日系數解析遞推偏轉距離的方法(后文簡稱Conway 模型);2007 年Izzo[26]基于幾何法推導了偏轉距離的解析模型(后文簡稱Izzo 模型);2008 年Vasile 和Colombo[20]基于高斯攝動方程推導了偏轉距離的解析模型(后文簡稱Vasile 模型),其相比于Conway 模型具有更高的計算效率.Vasile 模型是目前使用較為廣泛的解析偏轉模型,文獻[27-28]均基于Vasile 模型,以偏轉距離為指標優化了動能撞擊任務.但是目前尚無研究對偏轉距離解析模型的最優性進行分析.

小行星被偏轉前后,雖然軌道根數發生了改變,但到達地球近地點的時刻可能變化不大.本文首先研究了近地點時刻受偏轉影響的變化規律,提出在動能撞擊任務的設計初期,可以假定小行星被偏轉前后近地點時刻不變(近似偏轉模型),由此可以省略掉近地點時刻搜索的計算量.相比于解析模型,這種近似模型可以在提升最優性的同時、降低求解復雜性.考慮運載約束,將化學推進變軌簡化為脈沖推力變軌,建立了直接轉移(兩脈沖及三脈沖)和行星借力飛行轉移(單次及兩次借力)的動能撞擊軌道優化模型.以偏轉小行星Apophis 為例,設計了動能撞擊任務,研究了利用兩脈沖/三脈沖直接轉移、單次/兩次行星借力飛行技術是否可以提升偏轉距離.通過對比二體動力學模型與高精度動力學模型對計算偏轉距離的誤差,進一步探討了在動能撞擊任務設計過程中采用二體動力學模型的合理性.

1 偏轉距離的現有計算方法

本文中偏轉距離的定義為:偏轉前后小行星近地點(closest-approach,CA)距離的改變量(偏轉距離).偏轉距離δρCA的數學表達式為

其中,ρ(t)為t時刻小行星的地心位置矢量,a和b為時間搜索窗口的下界和上界.關于偏轉距離的計算方法,本節將介紹數值模型、兩種經典解析模型(Vasile 模型和Izzo 模型).

1.1 數值模型

對于小行星的高精度軌道演化動力學模型,本文考慮的攝動力有八大行星引力、月球引力、四大主帶小行星引力(Ceres,Vesta,Pallas,Hygiea)、相對論效應(general relativity,GR)、地球非球形引力J2、太陽光壓(solar radiation pressure,SRP)以及Yarkovsky效應(yarkovsky effect,YE).高精度動力學模型的數學表達式如下所示

其中,dpi表示第i個行星地日心矢量,ρpi表示第i個行星到小行星的矢量.行星和月球的星歷調用JPL DE430 歷表[29],小行星星歷調用JPL Horizons On-Line Ephemeris System[30]公布的真實星歷.動力學模型中考慮非引力項(太陽光壓、Yarkovsky 效應)的原因為:太陽光壓對一些尺寸較小的小行星影響較大,對于尺寸較大的小行星,Yarkovsky 效應不能忽略[31-33].利用Runge-Kutta-Fehlberg7(8)積分器[34-35]對軌道運動方程的一階形式進行數值積分,并采用一種一維最小值算法[36]求解近地點距離,因為這種算法不需要地心距離的一階導數信息.

為驗證本文高精度軌道預報方法的準確性,本節對小行星Apophis 在2029 年的近距離飛掠地球事件進行了預報,并與JPL 公布的結果[30]進行了對比.小行星Apophis 是目前行星防御領域高度關注的小行星,根據JPL[30]公布的最新信息(表1),Apophis 將在2029 年4 月13 日與地球近距離飛掠,近地點距離為0.00025AU,其質量約為6.1×107t.將Apophis 在2019 年1 月1 日00:00:00(TDB)的位置速度(坐標系為ICRF)作為積分起點.

表1 小行星Apophis 初始位置速度[30]Table 1 Initial state of Apophis

表2 給出了JPL 公布的結果與本文預報的結果.仿真結果表明,利用本文方法計算的Apophis 近地事件,時間誤差約為0.02 s,位置誤差約為1 km.

表2 小行星Apophis 近地事件預報結果Table 2 Results of Apophis’closest approach

理論上,在動能撞擊軌道優化問題中,需要利用高精度動力學模型進行遞推以保證最優性,但是各種攝動力的引入會顯著降低求解效率.采用二體模型(僅考慮太陽引力)替代多體模型,可以顯著降低軌道演化過程的計算量,但仍需數值計算小行星的近地點.二體模型的數學表達式如下所示

假定小行星被偏轉前后日心位置r(td)不變、僅日心速度v(td)發生改變,即偏轉使小行星狀態由(r(td),v(td))改變為(r(td),v′(td)).將偏轉前后小行星的狀態利用二體模型由偏轉時刻(td)解析遞推(fkepler)至近地點時刻(tCA),可求得小行星近地點距離.小行星被偏轉前近地點位置的日心矢量r(tCA)及地心矢量ρ(tCA)表示為

其中,近地點時刻tCA采用一維最小值算法[36]對式(2)進行求解,可得偏轉距離為

1.2 解析模型

1.2.1 Vasile 模型

Vasile 模型可以解析計算小行星被偏轉前后MOID (minimum orbit intersection distance)的改變量.其基本思想是將MOID 點處位置矢量的改變量分解為徑向、切向及法向分量

其表示為軌道根數改變量的函數如下

其中,θMOID表示近地點時小行星的日心真近點角,.軌道根數改變量與偏轉速度增量之間的關系可由高斯攝動方程計算,由于撞擊改變了小行星的軌道周期,因此近地點處的平近點角改變量還需考慮由于軌道周期改變而產生的長期演化影響,公式如下

其中,θd表示撞擊時小行星的日心真近點角,σ 表示撞擊時刻距離近地點時刻的時間.

雖然MOID 概念并不等同于近地點距離[20,23],但通過觀察公式可以發現,如果將MOID 點真近點角θMOID替換為近地點真近點角θCA,則該公式可以用于估算近地距離改變量.δα(td)表示偏轉位置處小行星軌道根數的改變量,δv(td)表示偏轉位置處小行星速度的改變量,δr(tCA)表示近地點位置的偏移量,ACA和Gd表示狀態轉移矩陣

近地點距離的改變量δρCA可表示為

其中,ρCA表示偏轉前的近地點矢量.

1.2.2 Izzo 模型

Izzo 模型的基本思想是,通過計算偏轉前后近地點時間的差異來計算小行星在B 平面上的偏移量?ξ.其中,B 平面定義為垂直于小行星進入地球影響球雙曲線速度U且經過地球質心的平面[37],?pik[38]在B 平面上定義了(ξ,η,ζ)坐標系:η 軸沿U的方向,ζ 軸為地球日心速度在B 平面上投影的反方向,ξ軸的指向由右手定則給出,在該坐標系下ξ 坐標可代表小行星軌道與地球軌道的MOID 值,B 平面上的偏移量?ξ 可用來描述小行星的偏轉距離.在日心視角下,小行星被偏轉前后近地點位置的幾何關系如圖1所示.

圖1 Izzo 模型原理示意圖Fig.1 Schematic diagram of Izzo model

圖1 中,?σ 表示偏轉前后近地點時刻的改變量,vast表示小行星的近地點日心速度,U表示小行星進入地球影響球時的相對速度,ψ 表示小行星速度vast與U之間的夾角,θaE表示地球速度vE與U夾角的補角.?ξ 表示小行星被撞擊后的偏轉距離,按照圖中的幾何關系,可建立偏移量?ξ 與交會時刻改變量?σ 之間的等式如下

經過對?σ 表達式的推導(具體步驟參考文獻[26]),可將上式整理如下

其中,Usc表示撞擊時撞擊器相對于小行星的速度.

Izzo 模型可以比較直觀地給出偏轉距離的影響因素,比如小行星本身的物理屬性(小行星質量mast、動量傳遞因子β 等)、小行星與地球的軌道關系(小行星半長軸aast、近地點時地球的速度vE等)以及軌道優化結果(撞擊器質量msc、偏轉時間σ、撞擊速度Usc等).該公式的最后一項點乘表示在偏轉時間較長的情況下,在小行星速度增量的切向分量是影響偏轉距離的主要因素.另外需注意到,Izzo 模型只能計算小行星在B 平面內的偏移量,并不等于近地點距離的改變量,因此當采用該模型估算真實近地點距離改變量時會出現偏差.圖2 展示了一種利用Izzo 模型計算產生偏差的情況.

在該情況下,即使偏轉使小行星產生了較大的?ξ,但實際上小行星被偏轉后近地距離減小了.也就是說,Izzo 模型無法回答小行星被撞擊后近地點距離會更近還是更遠的問題.

圖2 Izzo 模型產生偏差的情況Fig.2 Deviation of Izzo model

2 近地點時刻預估近似模型

如1.1 節所述,求解小行星近地點距離改變量的本質在于求出偏轉前后的近地點時刻.如果能夠提前估算出小行星被撞擊后近地點的時刻,那么就可以利用軌道解析遞推模型,快速計算出被偏轉后的近地點矢量.下面通過公式推導分析近地點時刻的影響因素.

以Izzo 模型為基礎,在已知小行星軌道根數的情況下,撞擊時刻與近地點時刻之間的時間間隔σ可由下式計算

其中,td為偏轉時刻,?Mast為平近點角差值,nast為轉速.對上式取微分并忽略掉短期影響項后得

式中,aast為小行星的軌道半長軸,vast為小行星被撞擊時的日心速度,δvt為小行星被撞后產生的切向速度增量.對上式進行整理后,可得近地點時間改變量的解析表達式如下

經過推導可以發現,在小行星速度增量的各分量中,主要是切向分量影響近地點時刻的改變量.以偏轉Apophis 為例,下面研究了當撞擊產生不同的切向速度時,對近地點時刻的影響.以Apophis 在二體模型下求得的2029 年近地點時刻為參考(2029-04-13 21:40:09),當撞擊使Apophis 產生不同的切向速度增量(δvt=0.2 ～1 mm/s)時,Apophis 近地點時間改變量?σ 隨偏轉時間σ (2 ～10 a)的變化趨勢,如圖3所示.

圖3 近地點時刻改變量的變化曲線Fig.3 Change of perigee epoch

圖3 中坐標為(6,17.55)點的含義為:在Apophis到達近地點前6 年(2023/4/15 21:40:08)時進行撞擊,如果撞擊使Apophis 產生的切向速度增量為1 mm/s,那么將會使Apophis 在2029 年的近地點時刻晚17.55 s,也就是Apophis 被撞擊后的近地點時刻為2029/4/13 21:40:26.從這幅圖,可以獲得兩方面的信息:第一方面,在10 年預警時間的情況下,單顆動能撞擊器能使Apophis 近地點時刻改變不到30 s;第二方面,由Izzo 模型可知,?σ 與偏轉距離成正比,因此隨著偏轉時間σ 的增大,近地點距離改變量(偏轉距離)整體呈現波動上升的趨勢.

下面僅著眼于偏轉后的軌道,分析以不同時刻計算近地點時,近地近地距與真實近地距的差異.小行星進入地球影響球后,由于相對速度較高(數千米每秒),其運動軌跡近乎直線.小行星被偏轉后,在地球影響球內的運動軌跡如圖4 所示.其中,U表示小行星進入地球影響球時的相對速度,表示被偏轉后的真實近地點時刻,表示近似近地點時刻,.近似后的近地距滿足下面等式

圖4 小行星真實近地距與近似近地距Fig.4 Real perigee distance and approximate perigee distance

根據1.1 節的數值模型,可求得Apophis 進入地球影響球的相對速度∥U∥≈5.9 km/s,假設被偏轉后的真實近地距離為36 000 km,計算當近地點時間改變量?σ ∈[0,30]s 時,近似近地距與真實近地距之間的相對誤差,計算結果如圖5 所示.

圖5 近似近地距與真實近地距之間的相對誤差Fig.5 Relative error between real perigee distance and approximate perigee distance

由圖5 可知,在近地點改變時間小于30 s 時,近似近地距與真實近地距之間的相對誤差不足0.01%.例如,當近地點時刻相差30 s 時,近似近地距為36 000.4 km,與真實近地距差別小于1 km.近地距越大,二者之間的誤差越小.

因此,對于本文偏轉Apophis 算例,在任務設計初期,可以忽略偏轉前后小行星近地點時間的改變量,即.小行星的近地點距離改變量可由下面的近似模型進行估算

其中∥ρ′(tCA)∥表示小行星被偏轉后,在偏轉前近地點時刻tCA處的地心矢量.依據該模型,已知撞擊后小行星的軌道根數,可以解析地快速求得tCA時刻軌道根數,省略了用于搜索近地點時刻的計算量.

3 動能撞擊任務設計與優化建模

動能撞擊任務流程圖如圖6 所示.本章將首先建立航天器(撞擊器)軌道轉移模型,然后建立動能撞擊任務優化模型,最后介紹本文采用的優化方法.

圖6 動能撞擊任務流程Fig.6 Process of a kinetic impactor mission

3.1 轉移模型

本文采用圓錐曲線拼接法求解行星際轉移軌道.圓錐曲線拼接法的基本原理為將轉移軌道按照天體引力影響球進行分段,航天器在影響球內只考慮中心天體的引力,最后通過航天器在影響球邊界處的狀態對多段圓錐曲線進行拼接,天體影響球邊界的相對速度稱為逃逸速度或雙曲線超速.該簡化算法在深空軌道設計中是通用且高效的[39].下面對兩脈沖轉移模型、三脈沖轉移模型、借力飛行轉移模型進行建模.

3.1.1 兩脈沖轉移模型

首先,航天器在t0時刻從地球逃逸,經過?t1天后與撞擊危地小行星.假定轉移期間不施加任何速度增量,由運載火箭直接發射送入撞擊軌道(圖7).

圖7 兩脈沖直接轉移模型Fig.7 Two-impulse direct transfer model

通過求解Lambert 問題[40],可以計算出在地球影響球邊界處的逃逸速度v∞以及航天器撞擊小行星的相對速度Usc,對應的數學表達式為

其中,r0和rf分別表示航天器在t0和t0+t1時的日心位置矢量,v0和vf為對應的日心速度矢量,vE和vast分別表示地球和小行星的日心速度矢量.v∞可由運載火箭提供,對應的C3為

假設轉移過程不施加任何機動,則航天器在撞擊時的質量等于發射入軌時的質量

其中,撞擊時刻td=t0+t1,msc(t0)=f(C3).f(C3)為運載能力擬合公式,CZ-5 運載火箭的運載能力曲線如圖8 所示.

圖8 長征五號運載能力曲線Fig.8 CZ-5’s launch performance

當動能航天器撞擊目標小行星后,超高速撞擊會導致小行星表面產生濺射物,一部分濺射物最終會再次吸積在小行星表面,另一部分速度較高的濺射物將會從小行星的引力場中逃逸,這種現象類似于火箭產生動力的機理,將會放大動量傳遞效應,可采用動量傳遞因子β 描述濺射物對動量傳遞的影響[11].假定航天器與小行星的撞擊過程為完全非彈性碰撞,根據動量守恒原理,小行星被撞擊后產生的速度改變量可由下式估算[41-42]

其中,msc(td)和vsc(td)表示航天器在撞擊小行星時的質量和速度,mast和vast(td)表示小行星被偏轉時的質量和速度.本文不考慮濺射物的影響,因此取β=1.

綜上所述,給定節點(t0,?t1)后,可以確定航天器在撞擊小行星前的速度vsc(td)及質量msc(td).由式(25)計算出小行星被偏轉后的狀態,從而可計算出小行星被偏轉前后的近地點距離改變量.

3.1.2 三脈沖轉移模型

航天器在t0時刻從地球逃逸,經過?t1天后撞擊危地小行星,期間在t0+α1?t1時施加一次深空機動DS M1,α1的取值范圍為0 到1(圖9).

圖9 三脈沖直接轉移模型Fig.9 Three-impulse direct transfer model

對于三脈沖轉移模型,需要給定航天器從地球逃逸時的速度大小及方向,對應3 個參數(C3,DLA,RLA).由于航天器在轉移過程中施加了深空機動,所以航天器在撞擊時的質量應減去轉移過程的燃料消耗

其中,撞擊時刻td=t0+t1,轉移速度增量總和?vtot=DS M1,msc(t0)=f(C3),Isp表示發動機比沖.

綜上所述,給定入軌參數(C3,DLA,RLA)、時間節點(t0,?t1,α1)后,可以確定航天器在撞擊小行星前的速度及質量,由式(25)計算出小行星被偏轉后的狀態,從而可計算出小行星被偏轉前后的近地點距離改變量.

3.1.3 借力飛行轉移模型

為了提高動能撞擊的偏轉效率,本文引入借力飛行技術,探討借力飛行技術是否對偏轉效率有所提升.借力飛行的基本原理是利用自然天體的引力,對航天器的速度進行改變,如圖10 所示.

圖10 借力飛行示意圖Fig.10 Schematic diagram of gravity assist

受自然天體引力作用,在自然天體引力影響球邊界處航天器的進入與離開之間的夾角為

其中,μga和Rga分別表示借力天體的引力常數和半徑,hp表示借力高度,vga表示借力天體的速度.定義一個新的坐標系o-ijk

因此,借力后矢量可描述如下

其中,φ 表示在j-k的投影與k軸正向之間的夾角,取值范圍為[0,2π).

本文考慮的借力天體有金星、地球和火星.首先,航天器在t0時刻從地球逃逸,經過?t1到達借力天體,再經過?t2天后撞擊危地小行星.轉移期間,還需執行兩次深空機動,DS M1施加于t0+α1?t1時以及DS M2施加于t0+?t1+α2?t2時,α1與α2的取值范圍為0 到1(圖11).

圖11 借力飛行轉移模型Fig.11 Gravity assist transfer model

同樣的,利用式(26)可以計算出航天器在撞擊小行星時的質量,其中撞擊時刻td=t0+t1+t2,轉移速度增量總和?vtot=DS M1+DS M2.

綜上所述,給定入軌參數(C3,DLA,RLA)、時間節點(t0,?t1,?t2,α1,α2)、借力參數(hp1,φ1)后,可以確定航天器在撞擊小行星前的速度及質量,由式(25)計算出小行星被偏轉后的狀態,從而可計算出小行星被偏轉前后的近地點距離改變量.以此類推,每增加一個借力天體,需增加4 個新的優化變量(?t,α,hp,φ).

3.2 優化模型

對于兩脈沖直接轉移模型,共有2 個決策變量

對于三脈沖直接轉移模型,共有6 個決策變量

對于單次借力飛行轉移模型,共有10 個決策變量

對于兩次借力飛行轉移模型,共有14 個決策變量

分別利用數值模型、Vasile 模型、Izzo 模型和近似模型構建目標函數(偏轉距離),其表達式如下所示:

(1)利用數值模型構建目標函數

(2)利用Vasile 模型構建目標函數

(3)利用Izzo 模型構建目標函數

(4)利用近地點時刻預估模型構建目標函數

3.3 優化算法

通過下述優化算法,可計算出動能撞擊任務中直接轉移方案/借力飛行轉移方案的最優轉移軌跡.圖12 給出了優化算法的流程圖.

圖12 優化算法流程圖Fig.12 Optimization process

第一步:輸入決策變量范圍;

第二步:生成初始種群;

第三步:計算發射C3,并判斷是否在擬合區間內(本文CZ-5 擬合曲線需C3<60 km2/s2).如果不在區間內,返回上一步;如果符合區間范圍,進入下一步;

第四步:利用直接轉移模型或行星借力借力飛行轉移模型計算小行星被偏轉后的狀態;

第五步:按照指定的目標函數式(34)～式(37),計算每個個體的適值函數;

第六步:判斷是否達到終止條件.如果未達到條件,進入第七步;如果達到終止條件,直接進入第八步;

第七步:選擇種群中的優質個體,進行遺傳運算,更新種群,返回第三步;

第八步:輸出適值度最高的個體,即輸出最優轉移軌跡.

4 仿真結果與討論

本文以偏轉Apophis 為例,在10 年預警時間的條件下,參考CZ-5 的運載能力,設計了直接轉移和行星借力飛行轉移的動能撞擊偏轉任務.

4.1 偏轉距離計算方法比較

本節以直接轉移方案為例,對比分析了利用不同模型進行優化的結果,對優化模型中的目標函數進行了選擇.

航天器逃逸時間t0的范圍為2020 年1 月1 日至2021 年1 月1 日,轉移時間?t1∈[100,1400]d.分別以數值模型、Vasile 模型、Izzo 模型、近似模型作為目標函數,利用遺傳算法對直接轉移方案進行了全局優化.該優化問題中,種群個體數設置為1000,種群代數設置為10.將優化出的發射窗口代入數值模型中計算數值偏轉距離,表3 對比了采用不同目標函數時求得的最優解.

通過表3 優化結果的對比,可以發現在優化問題中,以Vasile 模型和近似模型作為目標函數時,可以較好地替代數值模型.尤其對于近似模型,利用其優化出的發射窗口與數值模型的基本一致,同時計算時間相比于數值模型降低了約90%;Vasile 模型也可以較好還原優化結果,但優化出的發射窗口與數值模型的存在偏差,會損失一定的最優性;利用Izzo 模型作為目標函數,即使其目標函數值達到398 km,但優化出的解不具備參考性,這是由于小行星被偏轉前后出現了圖2 所示的位置關系.這也驗證了1.2.2 節中關于Izzo 模型初步分析的結論,不可將其作為優化的目標函數.從Izzo 模型的優化結果也可以得知,不管撞擊器對小行星的偏轉能力有多強,一旦選擇了錯誤的偏轉時機,則可能會給地球帶來更大的威脅.

為了直觀對比數值模型、解析模型、近似模型的發射窗口差異,下面在同一張圖中繪制出利用不同模型求解出的偏轉距離等于180 km 的等高線(圖13).

表3 基于不同偏轉距離模型的優化結果Table 3 Optimization results based on different deflection models

圖13 不同模型的180 km 等高線圖Fig.13 180 km contour maps of different models

圖13 中,黃色粗線表示利用數值模型求解的180km 等高線,紅色和藍色細線分表表示用近似模型和Vasile 模型求解的結果,“*”標識出的位置表示利用數值模型優化出的發射窗口.紅線與黃線以及藍線與黃線的交叉點,表示與數值模型誤差為0%的點.可以看出,紅線與紅線的擬合性較好,說明近似模型可以較好的替代數值模型;而藍線與黃線交叉的點較少,說明Vasile 模型的解與數值解存在一定偏差.同時,可以明顯的看到“*”標識落在藍線上,這也說明如果以Vasile 模型的解作為優化目標,優化出的窗口會與實際最優窗口有一定偏差.

4.2 直接轉移方案

4.2.1 兩脈沖轉移首先利用Pork-Chop 圖分析發射窗口,具體計算的方法為:將航天器逃逸時間t0的范圍選為2020年1 月1 日至2021 年1 月1 日,轉移時間?t1∈[600,1400]d,以1 d 為步長,利用近似模型計算所有情況產生的偏轉距離,最終繪制得偏轉距離等高線如圖14 所示.

圖14 Apophis 偏轉距離等高線圖(2020 年)Fig.14 Apophis deflection distance contour map(2020)

由圖14 中可知,在2020 年中,利用兩脈沖轉移方案偏轉Apophis 共有兩個發射窗口,一個位于2020年5 月2 日前后,另一個位于2020 年5 月13 日前后.這2020 年中兩個發射窗口的具體信息總結如表4.

可以看到,2 號窗口雖然擁有更大的撞擊器質量和Apophis ?v,但最終與1 號窗口的偏轉距離相差不大.這是由于2 號窗口所需的轉移時間更長,偏轉相位次于1 號窗口,從而導致Apophis 被偏轉后到達2029 年近地點的偏轉時間降低.因此,不能直接用撞擊體質量、速度增量等因素單方面來評估偏轉距離,而是需要綜合考慮多項影響因素,如Izzo 模型所示,包括小行星本身的物理屬性(小行星質量mast,β 等)、小行星與地球的軌道關系(小行星半長軸aast、近地點時地球的速度vE等)以及軌道優化結果(撞擊器質量msc、偏轉時間σ、撞擊速度Usc等).

表4 2020 年內兩個發射窗口的優化結果Table 4 Optimization results of two launch windows in 2020

以總任務時間較少的1 號窗口為例,直接轉移方案的轉移軌跡如圖15 所示.航天器自2020 年5 月2日從地球逃逸,經過670 d 后,在2022 年3 月4 日與Apophis 相撞.經過撞擊,Apophis 在2029 年的近地距可增加約192 km.

圖15 2020 年1 號發射窗口的兩脈沖直接轉移軌跡Fig.15 Two-impulse direct transfer trajectory of the 1st launch window in 2020

4.2.2 三脈沖轉移

將航天器逃逸時間t0的范圍選為2020 年1 月1日至2021 年1 月1 日,轉移時間?t1∈[100,1400]d.優化結果如表5 所示,轉移軌跡如圖16 所示.

表5 三脈沖直接轉移方案優化結果Table 5 Optimization results of three-impulse direct transfer trajectory

圖16 三脈沖直接轉移軌跡Fig.16 Three-impulse direct transfer trajectory

仿真結果表明,三脈沖直接轉移解與兩脈沖直接轉移的1 號窗口非常接近,三脈沖解雖然具有更高的發射入軌質量,但是由于施加了中途機動,同時又會損失一部分撞擊器質量,利用三脈沖直接轉移方案對偏轉距離不會有明顯提升.

4.3 借力飛行轉移方案

航天器逃逸時間t0范圍為2020 年1 月1 日至2021 年1 月1 日,轉移時間?t∈[100,1000]d,hp∈[100,3000]km,φ ∈[0,2π]rad,C3∈[0,50]km2/s2,DLA∈[0,π]rad,RLA∈[0,2π]rad.分別設計了如下的借力飛行偏轉方案:金星、地球、火星、金星?金星、金星?地球、金星?火星、地球?金星、地球?地球、地球?火星、火星?金星、火星?地球、火星?火星.遺傳算法的參數設置:單次借力的種群個體數設置為5000,種群代數設置為30;兩次借力的種群個體數設置為8000,種群代數設置為30,設計結果見表6.

表6 借力飛行轉移方案優化結果Table 6 Gravity assist transfer trajectory optimization results

仿真結果表明,只有單次地球借力的轉移方案對偏轉距離有所提升,相比直接轉移方案提升了約20 km.關于單次地球借力轉移方案的細節如表7所示.

表7 地球借力轉移方案優化結果Table 7 Optimization results of Earth gravity assist transfer trajectory

地球借力轉移方案的轉移軌跡如圖17 所示.航天器自2020 年2 月19 日從地球逃逸,2021 年5 月17 日進行地球借力,2023 年1 月25 日與Apophis 相撞,航天器轉移時間共計1071 d.經過撞擊,Apophis在2029 年的近地距可增加約210 km.

圖17 地球借力轉移軌跡Fig.17 Earth gravity assist transfer trajectory

由借力飛行轉移方案的設計結果可知,利用借力飛行未必能夠提升偏轉距離.這是因為借力飛行概念的初衷是降低深空探測過程的燃料消耗,而對于動能撞擊這種特定的深空探測任務來說,其性能指標是偏轉距離,這是一個需要綜合考慮燃料消耗、轉移時間、撞擊相位等因素的指標.如果在航天器轉移的過程中執行不當的借力飛行,航天器為了與借力天體相位匹配,可能會錯過最佳發射窗口,從而降低偏轉距離.

5 二體模型對防御窗口的影響

以上仿真均是在二體模型的假設下進行的,本章將主要討論基于二體模型與基于高精度模型優化結果的差別.

實際上,在高精度模型中某些攝動力對于偏轉距離的影響非常微弱,在優化過程中是可以忽略的因素.因此,首先對小行星Apophis 在運動過程中所受的攝動力進行分析,篩選出對偏轉距離影響最大的攝動力,對高精度模型進行簡化,從而減少優化過程中不必要的計算量.以小行星Apophis 為例,下圖給出了Apophis 在2019 年1 月1 日至2030 年1 月1 日之間各項攝動加速度的變化趨勢.其中,圖18為內太陽系4 顆行星、廣義相對論、太陽光壓以及Yarkovsky 效應的攝動加速度變化曲線;圖19 為外太陽系5 顆行星、3 顆矮行星、4 顆主帶小行星、地球非球形引力J2 的攝動加速度變化曲線.

圖18 Apophis 攝動力量級(1)Fig.18 Orders of magnitude of Apophis perturbation(1)

由圖18 可知,Apophis 在2029 年到達近地點之前,金星、地球、火星和月球的攝動加速度會發生量級上的顯著變化.不同攝動源按照攝動加速度量級的分布區間總結于表8 中.

圖19 Apophis 攝動力量級(2)Fig.19 Orders of magnitude of Apophis perturbation(2)

表8 Apophis 攝動力量級Table 8 Orders of magnitude of Apophis perturbation

基于3.1.1 節所述的兩脈沖轉移模型,以第4 章兩脈沖轉移仿真結果中的1 號窗口(2022 年5 月2日)為例,在2022 年3 月2 日對小行星Apophis 施加0.38 mm/s 的速度增量,分別利用高精度模型和二體模型對偏轉后的小行星進行遞推,計算二體模型與高精度模型對計算偏轉距離的誤差.仿真發現,利用二體模型計算偏轉距離為192.03 km,利用高精度模型計算的偏轉距離176.27 km.進一步,在二體模型上每次增加一項攝動力,用不同力模型遞推求得的偏轉距離,結果如表9 所示中,發現對偏轉距離影響最大的攝動力為地球和金星引力.在考慮金星+地球攝動模型基礎上,表10 每次增加一項攝動力,分析除金星和地球外,對偏轉距離影響最大的因素.由表9和表10 可知,對小行星Apophis 偏轉距離影響最大的因素為金星、地球、木星.在二體模型基礎上引入地球、金星、金星的引力攝動后,計算出的偏轉距離與高精度解的誤差在1%以內.

在保證計算精度的前提下,為了提高計算速度,本文引入金星、地球、木星引力攝動對二體模型修正,利用修正模型替代高精度模型.修正模型的數學表達式如下

表9 在二體模型基礎上每次增加一項攝動力Table 9 Add one perturbation source at a time based on the two-body model

表10 在金星+地球攝動模型基礎上每次增加一項攝動Table 10 Add one perturbation source at a time based on the Venus+Earth model

以直接轉移方案為例,分別基于二體模型和修正模型對偏轉Apophis 的動能撞擊任務進行優化.航天器的逃逸時間t0的范圍選為2020 年1 月1 日至2021 年1 月1 日,轉移時間為?t1∈[600,1400]d.種群個體數設置為200,種群代數設置為10.利用二體模型及修正模型優化的結果如表11 所示.提取出二者優化出的防御窗口,代入高精度模型計算出對應的高精度結果.

表11 二體模型與修正模型優化解對比Table 11 Comparison of optimization results between two-body model and modified two-body model

仿真結果表明,利用二體模型進行優化,在不影響防御窗口的情況下,計算時間會比精度更高的修正模型更少.在動能撞擊任務前期設計過程中,可以利用二體模型快速評估防御效果.

6 結論

本文對動能撞擊小行星防御任務的脈沖軌道優化方法開展了研究.為提升動能撞擊任務的求解效率,將兩種經典的偏轉距離解析模型引入優化模型,同時提出一種基于近地點時刻預估的偏轉距離近似模型;考慮運載約束,建立了直接轉移和行星借力飛行轉移的動能撞擊優化模型.

仿真結果表明,相比于經典的偏轉距離解析模型,本文提出的基于近地點時刻預估的近似模型更加簡潔直觀,在提升計算效率的同時,也能保證優化問題的最優性.三脈沖直接轉移方案與兩脈沖直接轉移方案的最優偏轉效果基本一致,借力飛行轉移方案相比于直接轉移方案對偏轉距離的提升效果并不明顯.在動能撞擊任務的前期設計中,可以基于二體模型進行防御效果的快速評估,雖然對偏轉距離有一定誤差,但對防御窗口的優化結果影響不大.