基于偏轉(zhuǎn)距離近似模型的動能撞擊小行星防御任務(wù)脈沖軌道優(yōu)化研究1)

王藝睿 李明濤 周炳紅

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

引言

太陽系內(nèi)的小行星主要分布于火星軌道與木星軌道之間,也有一部分小行星會穿越地球軌道威脅到地球上的生命.諸多嚴(yán)重的小行星撞擊地球事件(如6500 萬年前的K-T 事件[1]、1908 年通古斯大爆炸事件[2]、2013 年車?yán)镅刨e斯克事件[3]等)已引發(fā)人類對于行星防御的關(guān)注.在能夠?qū)匦⌒行翘崆邦A(yù)警的前提下,摧毀小行星結(jié)構(gòu)或者偏轉(zhuǎn)小行星軌道是避免其撞擊地球的兩種基本方式[4].目前已提出的防御手段主要有核爆[5]、動能撞擊、激光燒蝕[6]、拖船[7]、引力拖船[8]、離子束牽引[9]等,其中動能撞擊技術(shù)是目前最易實現(xiàn)且成熟度較高的手段.

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

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

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

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

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

1 偏轉(zhuǎn)距離的現(xiàn)有計算方法

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

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

1.1 數(shù)值模型

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

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

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

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

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

表2 小行星Apophis 近地事件預(yù)報結(jié)果Table 2 Results of Apophis’closest approach

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

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

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

1.2 解析模型

1.2.1 Vasile 模型

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

其表示為軌道根數(shù)改變量的函數(shù)如下

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

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

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

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

其中,ρCA表示偏轉(zhuǎn)前的近地點矢量.

1.2.2 Izzo 模型

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

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

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

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

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

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

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

圖2 Izzo 模型產(chǎn)生偏差的情況Fig.2 Deviation of Izzo model

2 近地點時刻預(yù)估近似模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 動能撞擊任務(wù)設(shè)計與優(yōu)化建模

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

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

3.1 轉(zhuǎn)移模型

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

3.1.1 兩脈沖轉(zhuǎn)移模型

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

圖7 兩脈沖直接轉(zhuǎn)移模型Fig.7 Two-impulse direct transfer model

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

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

假設(shè)轉(zhuǎn)移過程不施加任何機動,則航天器在撞擊時的質(zhì)量等于發(fā)射入軌時的質(zhì)量

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

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

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

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

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

3.1.2 三脈沖轉(zhuǎn)移模型

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

圖9 三脈沖直接轉(zhuǎn)移模型Fig.9 Three-impulse direct transfer model

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

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

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

3.1.3 借力飛行轉(zhuǎn)移模型

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

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

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

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

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

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

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

圖11 借力飛行轉(zhuǎn)移模型Fig.11 Gravity assist transfer model

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

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

3.2 優(yōu)化模型

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

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

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

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

分別利用數(shù)值模型、Vasile 模型、Izzo 模型和近似模型構(gòu)建目標(biāo)函數(shù)(偏轉(zhuǎn)距離),其表達式如下所示:

(1)利用數(shù)值模型構(gòu)建目標(biāo)函數(shù)

(2)利用Vasile 模型構(gòu)建目標(biāo)函數(shù)

(3)利用Izzo 模型構(gòu)建目標(biāo)函數(shù)

(4)利用近地點時刻預(yù)估模型構(gòu)建目標(biāo)函數(shù)

3.3 優(yōu)化算法

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

圖12 優(yōu)化算法流程圖Fig.12 Optimization process

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

第二步:生成初始種群;

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

第四步:利用直接轉(zhuǎn)移模型或行星借力借力飛行轉(zhuǎn)移模型計算小行星被偏轉(zhuǎn)后的狀態(tài);

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

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

第七步:選擇種群中的優(yōu)質(zhì)個體,進行遺傳運算,更新種群,返回第三步;

第八步:輸出適值度最高的個體,即輸出最優(yōu)轉(zhuǎn)移軌跡.

4 仿真結(jié)果與討論

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

4.1 偏轉(zhuǎn)距離計算方法比較

本節(jié)以直接轉(zhuǎn)移方案為例,對比分析了利用不同模型進行優(yōu)化的結(jié)果,對優(yōu)化模型中的目標(biāo)函數(shù)進行了選擇.

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

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

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

表3 基于不同偏轉(zhuǎn)距離模型的優(yōu)化結(jié)果Table 3 Optimization results based on different deflection models

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

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

4.2 直接轉(zhuǎn)移方案

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

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

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

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

表4 2020 年內(nèi)兩個發(fā)射窗口的優(yōu)化結(jié)果Table 4 Optimization results of two launch windows in 2020

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

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

4.2.2 三脈沖轉(zhuǎn)移

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

表5 三脈沖直接轉(zhuǎn)移方案優(yōu)化結(jié)果Table 5 Optimization results of three-impulse direct transfer trajectory

圖16 三脈沖直接轉(zhuǎn)移軌跡Fig.16 Three-impulse direct transfer trajectory

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

4.3 借力飛行轉(zhuǎn)移方案

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

表6 借力飛行轉(zhuǎn)移方案優(yōu)化結(jié)果Table 6 Gravity assist transfer trajectory optimization results

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

表7 地球借力轉(zhuǎn)移方案優(yōu)化結(jié)果Table 7 Optimization results of Earth gravity assist transfer trajectory

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

圖17 地球借力轉(zhuǎn)移軌跡Fig.17 Earth gravity assist transfer trajectory

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

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

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

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

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

由圖18 可知,Apophis 在2029 年到達近地點之前,金星、地球、火星和月球的攝動加速度會發(fā)生量級上的顯著變化.不同攝動源按照攝動加速度量級的分布區(qū)間總結(jié)于表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 節(jié)所述的兩脈沖轉(zhuǎn)移模型,以第4 章兩脈沖轉(zhuǎn)移仿真結(jié)果中的1 號窗口(2022 年5 月2日)為例,在2022 年3 月2 日對小行星Apophis 施加0.38 mm/s 的速度增量,分別利用高精度模型和二體模型對偏轉(zhuǎn)后的小行星進行遞推,計算二體模型與高精度模型對計算偏轉(zhuǎn)距離的誤差.仿真發(fā)現(xiàn),利用二體模型計算偏轉(zhuǎn)距離為192.03 km,利用高精度模型計算的偏轉(zhuǎn)距離176.27 km.進一步,在二體模型上每次增加一項攝動力,用不同力模型遞推求得的偏轉(zhuǎn)距離,結(jié)果如表9 所示中,發(fā)現(xiàn)對偏轉(zhuǎn)距離影響最大的攝動力為地球和金星引力.在考慮金星+地球攝動模型基礎(chǔ)上,表10 每次增加一項攝動力,分析除金星和地球外,對偏轉(zhuǎn)距離影響最大的因素.由表9和表10 可知,對小行星Apophis 偏轉(zhuǎn)距離影響最大的因素為金星、地球、木星.在二體模型基礎(chǔ)上引入地球、金星、金星的引力攝動后,計算出的偏轉(zhuǎn)距離與高精度解的誤差在1%以內(nèi).

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

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

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

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

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

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

6 結(jié)論

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

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