劉世興，史東華，陳菊，易中貴，王本亮，龔自正，郭永新*

(1.遼寧大學物理學院，沈陽 110036；2.遼寧大學空間科學與技術研究院，沈陽 110036；3.北京理工大學數(shù)學與統(tǒng)計學院，北京 100081；4.北京理工大學宇航學院，北京 100081；5.北京衛(wèi)星環(huán)境工程研究所，北京 100094)

1 引言

“杞人憂天”出自《列子·天瑞》：杞國有人憂天地崩墜，身亡所寄，廢寢食者。雖然這里的“杞人”是被人嘲諷的對象，但實際上，小行星撞擊地球的事件在歷史上時有發(fā)生。夜空劃過的美麗流星，每年11月份出現(xiàn)的獅子座流星雨等都是流星顆粒進入地球大氣層和空氣摩擦而形成的。我們將圍繞太陽運行、尺寸在1m～800km內，且不容易釋放出氣體和塵埃的天體，稱作小行星。太陽系中存在大量的小行星，主要分布在火星與木星軌道之間的小行星帶和海王星外的柯伊伯帶，如圖1所示。在天文學上，將距離地球最小距離在0.3AU(1AU=1.496×108km)范圍內的小行星稱為近地小行星。地球曾發(fā)生過22次不同程度的生物滅絕事件，其中至少10次以上是由小行星撞擊地球所引起[1]，如6500萬年前的恐龍滅絕事件，1908年的通古斯空爆事件，2013年的俄羅斯車里雅賓斯克地區(qū)空爆事件等[2-5]。僅2021年內全球已發(fā)生近地小行星飛掠地球事件1605次，觀測到29顆小行星進入地球大氣層發(fā)生火流星事件。我國境內也經常發(fā)生火流星事件，比如，2014年11月5日在內蒙古錫林郭勒，2017年10月4日在云南香格里拉，2018年在云南西雙版納，2019年10月11日在吉林松原地區(qū)，2020年在青海玉樹地區(qū)，2021年在河南駐馬店等發(fā)生了火流星事件。在2022年3月11日，編號為2022 EB5，直徑大約3m的小行星在冰島北部上空進入了地球大氣層。雖然上述事件均未造成災難，但已威脅到了人們的生產和生活安全。由此可見，近地小行星撞擊風險概率雖小，但危害極大，但近地小行星撞擊地球的危害可測可防。防范近地小行星撞擊風險，具有非常迫切的現(xiàn)實需求和深遠的戰(zhàn)略意義[6-10]。

圖1 太陽系中的八大行星和近地小行星帶分布圖

應對近地小行星撞擊風險，給國際天文界、航天界等多個領域帶來重大科學和技術挑戰(zhàn)。21世紀以來，國際社會高度重視小行星撞擊風險防御工作，聯(lián)合國成立多個機構專門負責協(xié)調應對近地小行星撞擊風險，各主要航天國家也積極致力于相關工作，并有力牽引、推動了空間科學研究、空間技術探索和空間資源利用。我國近地小行星撞擊風險應對工作起步較晚、基礎薄弱、國際貢獻率低、國際話語權小，這與我國航天大國地位不符[8]。近幾年我國已經開始開展近地小行星防御問題研究，受到國家高度重視，并已上升為國家戰(zhàn)略[7]，且在小行星的軌道預測、領航監(jiān)測器的軌道設計與控制、小行星的在軌處置和監(jiān)測預警等方面取得了一定的研究成果[11-24]。應對小行星撞擊風險主要涉及兩個方面：(1)監(jiān)測預警。搜索發(fā)現(xiàn)、跟蹤定軌和數(shù)據(jù)更新、物性測量、計算撞擊概率、預估撞擊風險走廊和落點，預報撞擊風險。即發(fā)現(xiàn)小行星和確定撞擊風險。(2)在軌處置。依據(jù)監(jiān)測預警信息，發(fā)射防御航天器，改變小行星軌道避免其撞擊地球(即偏轉，Deflection)，或將小行星分裂成碎片避免或降低其對地球的危害(即摧毀，Disruption)。這里無論是小行星軌道精確預報、在大天體攝動影響下的軌道不確定性，還是在人為外力作用下(如動能撞擊等瞬時作用，或拖曳等長期作用下)的小行星軌道偏移、防御航天器自身的深空自主導航與控制技術等，都離不開對小行星運行軌道的精確計算測定，離不開精確的軌道設計與動力學研究。小行星的軌道動力學問題是防御的基礎、前提和核心問題。精確的軌道設計與計算需要采用基于約束系統(tǒng)的幾何動力學方法進行建模，同時采用高精度的幾何數(shù)值積分方法進行數(shù)值模擬。而軌道動力學一直是分析力學發(fā)展進程中的核心問題之一，但它針對航天工程提出的更復雜的軌道動力學問題研究的不夠深入，特別是航天工程中需要考慮到小行星的非規(guī)則形狀和自轉、引力源非球對稱性對牛頓引力場的修正、太陽引力場的相對論修正、小行星運行環(huán)境的隨機性等問題，給小行星軌道動力學建模、軌道高精度計算以及在軌處置航天器(群)的控制提出了更高要求，也是分析力學面臨的新課題。

2 幾何力學與幾何數(shù)值積分的發(fā)展現(xiàn)狀

現(xiàn)代分析力學最具標志性的成果是幾何力學的發(fā)展，以及在此基礎上構造的保結構算法——幾何數(shù)值積分方法，將二者與控制結合起來所興起的計算幾何力學與控制使得若干動力學系統(tǒng)實現(xiàn)了長時間、高精度的數(shù)值建模分析，也將為小行星運行軌道的高精度數(shù)值模擬與測量、航天器(群)的精確軌道設計、動力學建模與控制提供一種理論基礎。同時，小行星軌道動力學建模和高精度計算也必將進一步完善分析力學的高精度算法，特別是Hamel場變分積分子算法。

幾何力學是基于流形、纖維叢、Lie群與Lie代數(shù)等現(xiàn)代微分幾何工具，將經典分析力學升級成無坐標語言的整體表述形式，形式更加簡潔優(yōu)美，為多體系統(tǒng)、流體、彈性力學、場論、量子系統(tǒng)及幾何控制理論等提供了統(tǒng)一的框架，它在空間探索的軌道設計、機器人運動控制、等離子體、超流、洋流動力學和氣候建模、計算解剖學等諸多領域已經得到成功的應用[25,26]。利用現(xiàn)代微分幾何語言表述Lagrange力學和Hamilton力學有以下幾點優(yōu)勢：首先，適用于位形流形上的一般力學系統(tǒng)，以及由纖維叢上聯(lián)絡所刻畫的約束力學系統(tǒng)[27,28]。聯(lián)絡可以自然地描述由非完整約束所決定的不可積分布，通過它將系統(tǒng)的基底變量和纖維變量聯(lián)系起來，并具有鮮明的物理意義，進而采用規(guī)范理論的觀點統(tǒng)一且有力地處理非完整力學和非線性控制理論中的諸多問題，包括力學系統(tǒng)的幾何表示和幾何結構、穩(wěn)定性和鎮(zhèn)定、最優(yōu)控制等，在車輛、機器人等工程領域有重要的應用；其次，方便用于描述系統(tǒng)的大范圍運動，也便于進行幾何類比而將有限維連續(xù)系統(tǒng)的結構推廣至無窮維和離散力學系統(tǒng)[29,30]。特別地，對于位形空間是Lie群情況，采用對稱性約化理論[31]可實現(xiàn)代數(shù)空間上的線性化，還可以得到判別穩(wěn)定性的能量—動量方法。最后，基于離散Hamilton原理的變分積分子構造保持力學系統(tǒng)幾何結構的數(shù)值格式——幾何數(shù)值積分算法，為獲得長時間高精度的動力系統(tǒng)模擬提供有效的方法[32-36]。將計算幾何力學同幾何控制理論結合，得到能長時間保持高精度的數(shù)值控制算法，這一研究將為精確、高效、實時地控制航天器、潛水器、各類機器人等現(xiàn)代工程系統(tǒng)提供合理的設計原理。

目前，在無窮維系統(tǒng)中的幾何力學和幾何數(shù)值積分算法是國內外相關領域的研究熱點，這既是分析力學研究的邏輯延伸，也是若干工程科學的實際需要，這其中Hamel形式及Hamel場變分積分子在動力學系統(tǒng)建模和高精度計算中發(fā)揮核心作用。采用無坐標的幾何語言，可以通過類比，將具有對稱性的有限維力學系統(tǒng)的成熟方法推廣到無窮維系統(tǒng)，特別是借助于幾何類比，可將分析力學中不依賴于坐標的活動標架法推廣到無窮維系統(tǒng)和場論情況，通過直接對非物質速度進行變分的方式構造Hamel場變分積分子[36,37]，這個算法的能量保持顯著優(yōu)于Shake和Rattle算法的數(shù)值特性[38]；對基于SE(3)的幾何精確梁進行動力學分析，驗證了該算法能夠長時間保持幾何精確梁的能量與動量[39,40]；從協(xié)變場論觀點與時空離散角度研究Hamel場變分積分子算法的保動量性質，在證明該算法保幾何結構性質的同時，間接說明該方法具有長時間數(shù)值穩(wěn)定的性質[41]。此外，Hamel場變分積分子在計算無窮維非完整力學系統(tǒng)時具有明顯優(yōu)勢[42]，該方法表述Chaplygin-Timoshenko雪橇模型時表述簡單，便于分析該動力學方程解的適定性和解釋其特有的非線性行為。Hamel場變分積分子已用在工業(yè)機械臂、球面上移動機器人等的最優(yōu)控制，展現(xiàn)出了較好的數(shù)值特性[39,43,44]。

盡管Hamel場變分積分子算法在無窮維系統(tǒng)的動力學建模和高精度計算中得到了驗證，并且有理由相信不依賴于坐標選取的Hamel形式適于描述不規(guī)則形狀和自轉的小行星軌道動力學和姿-軌耦合問題，但針對小行星的精確軌道計算和預測，還需要考慮各種實際非理想因素，如小行星不規(guī)則形狀、后牛頓效應的攝動等對構造Hamel變分積分子的影響。此外，小行星的運行環(huán)境是動態(tài)變化的，甚至是隨機的，盡管在有關隨機微分方程的處理上，幾何力學和幾何數(shù)值積分方法具有明顯的優(yōu)勢[45]，但是Hamel場變分積分子對隨機問題和基于學習的離散場論保結構算法的拓展還缺乏研究。

3 近地小行星軌道動力學高精度建模研究進展

小行星在太陽系中的運動非常復雜，存在各種攝動因素的影響，如各大行星和小行星帶上質量較大小行星的質點引力攝動、太陽和大行星扁率產生的引力攝動、太陽引力場的相對論效應、小行星自轉和形狀不規(guī)則以及表面熱物理環(huán)境的差異導致的亞爾科夫斯基(Yarkovsky)效應所產生的非引力攝動項等。因此，小行星的運動不是簡單的Kepler問題，要精確計算和預測小行星運行軌道，一方面需要考慮諸多影響小行星運行的攝動因素，包括小行星在中心天體(太陽)和第三天體(大行星)的不規(guī)則形狀造成的引力攝動[46-49]。為精確計算小行星軌道，并考慮大天體(八大行星和月球)與主帶67顆質量較大的小行星的質點引力攝動、地球扁率項的引力攝動、太陽質心的后牛頓效應、太陽輻射壓對小行星軌道的攝動，通過引入測量隨機誤差、利用運動方程的積分計算、觀測量的理論歸算等求解小天體運行軌道[12, 14]。文獻[14]考慮各大行星的引力攝動、后牛頓效應等因素，采用改進的顯式辛算法和嵌套的RKF7(8)積分器對43顆已命名的近地小行星軌道演化進行數(shù)值研究，在一定程度上彌補了對于大偏心率問題采用定步長的不足。除此之外，由于小行星自轉、形狀不規(guī)則以及表面熱物理環(huán)境的差異導致了Yarkovsky效應也要產生非引力攝動項。另一方面在計算方法上，構造包含非保守力項的高階Galerkin變分積分子能夠處理小天體姿-軌耦合、太陽輻射力等非保守力項的情況[50]，而使用Chebyshev控制點與Lagrange插值方法構造的高階譜變分積分子比常規(guī)的Galerkin變分積分子具有更好的穩(wěn)定性[51]。采用變步長的數(shù)值算法，構造精細的力學模型并減少計算機截斷誤差，以達到小行星軌道計算精度和效率要求[11, 13]。雖然以上成果對影響小行星運行軌道的各種攝動因素有所研究，但是針對不同的小行星，對這些引力攝動和非引力效應如何對比觀測數(shù)據(jù)進行系統(tǒng)考慮，以提高計算精度和效率，還缺乏深入研究?，F(xiàn)有的小行星運行軌道的計算方法難以在動態(tài)環(huán)境下實現(xiàn)長時間高精度的效果，缺乏對復雜攝動因素的系統(tǒng)分析。此外，小行星運行環(huán)境瞬息萬變，各種攝動要素的影響交替變化，具有一定的隨機性。為此需要采用分析力學的幾何數(shù)值積分方法，特別是Hamel場變分積分子算法，確保在引力和非引力攝動效應下的長時間高精度計算需要加強幾何數(shù)值積分方法與隨機問題的結合研究。

4 近地小行星探測/在軌處置航天器(群)軌道動力學與控制研究進展

近地小行星具有弱引力場、形狀不規(guī)則的特點，從而使其外部引力場環(huán)境極其復雜，從而導致其運動形態(tài)復雜多樣。這些復雜性導致軌道優(yōu)化算法對初始值極其敏感，為探測/在軌處置航天器(群)的軌道設計和控制帶來了極大困難和挑戰(zhàn)[52-54]。

4.1 小行星處置航天器軌道設計與最優(yōu)控研究進展

目前，近地小行星處置航天器的軌道設計通常利用Lagrange點L1和L2 周圍的周期軌道的穩(wěn)定和不穩(wěn)定流形，揭示相空間中的管道網絡，這解釋了圓型限制性三體問題(CR3BP)中的遷移現(xiàn)象[55-59]。從不變流形理論出發(fā)，利用行星際高速公路來設計返回的轉移軌道[15]，基于三體問題的不變流形設計低成本登月軌道[16,17]。最近發(fā)現(xiàn)了太陽系中的一個從小行星帶延伸到天王星及更遠地方的流形管道結構[60]；利用Lagrange擬序結構(LCS)理論[61]計算橢圓型限制性三體問題(ER3BP)中的遷移勢壘，證明了ER3BP相空間中存在周期性脈動形式的Lagrange擬序結構，同時限制性四體問題中也存在依賴時間的不變流形[19-20]。將不變流形理論和離散力學與最優(yōu)控制方法(DMOC)結合，成功計算出從地球到月球的具有脈沖推力的優(yōu)化軌跡[62]。文獻[18]綜述了利用不變流形理論和弱穩(wěn)定邊界理論研究三體系統(tǒng)中低能量轉移與捕獲軌道設計的研究進展。在航天器群重構制導方面，基于笛卡爾坐標的重構技術，采用隨機樹搜索、遺傳算法、非線性優(yōu)化技術、序列二次規(guī)劃、相關狀態(tài)Riccati方程(SDRE)技術和腦風暴優(yōu)化算法等，計算了跟隨航天器從給定初態(tài)到給定末態(tài)的最省能量控制軌跡[63-68]。

4.2 伴飛航天器軌道動力學的綜合建模研究進展

考慮伴飛航天器受太陽光壓、其它星體引力攝動、隨機擾動等綜合影響的高精度軌道動力學模型具有重要意義。已有工作主要考慮類地J2攝動而忽略其它的引力項對小行星弱引力場中的攝動軌道建立引力攝動模型[69]，也有僅考慮太陽光壓效應而不考慮其它因素，建立繞小行星的軌道攝動模型，用于小行星航天器群編隊控制開發(fā)[70]。為實現(xiàn)從密切狀態(tài)到平均狀態(tài)的轉換，現(xiàn)有的轉換公式往往也是需要類地的關于引力場J2項的假設, 利用一個基于J2冪展開的Hamilton量來轉換密切根數(shù)和平均根數(shù)，需假設J2比余緯向勢大幾個數(shù)量級，有的工作甚至忽略了除J2之外的任何攝動項[71]；對于中心天體的旋轉速率等于或大于軌道平均運動，使用頻率分解從相關信號中剔除密切(短周期)頻率的做法計算代價大，且需頻繁測量來檢測短周期頻率[72]。另一種思路是設計消除密切效應的濾波器，利用Kalman濾波器，考慮二階引力勢從平均絕對根數(shù)中重構密切絕對根數(shù)[73]。

4.3 探測/在軌處置航天器群編隊制導與控制研究進展

近地小行星防御需要航天器群在深空軌道高精度編隊飛行，與繞地軌道的編隊場景(星鏈(圖2)、GPS(圖3))相比，具有更大的挑戰(zhàn)性，描述其相對運動的方程更為復雜，必須考慮攝動影響，編隊高精度控制的連續(xù)小推力往往控制在納牛到微牛的量級[21]。要完成航天器群小行星防御任務需要解決兩個基本問題：編隊制導和編隊控制，前者指集群軌道設計問題，后者為集群在軌保持或重構隊形。在編隊制導方面，常見的方法通常使用線性規(guī)劃、迭代學習控制、線性二次調節(jié)器(LQR)、相關狀態(tài)Riccati方程(SDRE)、Lagrange乘子、模型預測控制和感知框架等技術，通過跟蹤預定義的參考點實現(xiàn)[74-81]。另外，基于相對軌道根數(shù)的動力學模型是將編隊制導描述為跟蹤特定相對軌道根數(shù)的形式[82-84]。在編隊控制方面，可利用人工勢場來實現(xiàn)重構[85,86]，也有使用多種生物優(yōu)化技術來確定最優(yōu)的重構方案[87,88]，還有提出一種基于相對軌道根數(shù)的解析式燃料最優(yōu)脈沖編隊重構策略[89]，以及基于相對軌道根數(shù)的ΔV最優(yōu)方案和相關軌道動力學的精確線性化方法[90,91]。繞地球運行的集群編隊控制方法，需要考慮大氣阻力和J2對相對運動的影響[82-92]。為避免集群中個體的相對漂移、內部的碰撞，并減少整個任務中使用的推進劑，可以采用幾種新的開環(huán)控制方法[91]。另外，編隊保持方法也可針對特定軌道根數(shù)來設置[93,94]。

圖2 美國星鏈計劃[95]

圖3 GPS星座[96]

5 問題與結論

小行星防御涉及數(shù)學、物理、化學、力學和天體力學等多個學科領域，是一個復雜的系統(tǒng)工程，但其防御離不開四步：一來、二去、三繞、四控。小行星防御過程如圖4所示。

圖4 小行星防御過程示意圖

首先，所謂的“一來”，即如何對近地小行星軌道進行精確計算與預測問題。盡管國內外對影響小行星運行軌道的各種攝動因素有所研究，包括各大行星和小行星帶上質量較大小行星的質點引力攝動、太陽和大行星扁率產生的引力攝動、太陽引力場的相對論效應、小行星自轉和形狀不規(guī)則以及表面熱物理環(huán)境的差異導致的Yarkovsky效應所產生的非引力攝動項等，但是針對不同的小行星，對這些引力攝動和非引力效應如何對比觀測數(shù)據(jù)進行系統(tǒng)考慮，以提高計算精度和效率，還缺乏深入研究。同時，現(xiàn)有的小行星運行軌道的計算方法難以在動態(tài)環(huán)境下實現(xiàn)長時間高精度的效果，影響了復雜攝動因素的系統(tǒng)分析，為此需要采用分析力學的幾何數(shù)值積分方法，特別是Hamel場變分積分子算法，確保在引力和非引力攝動效應下的長時間高精度計算。此外，小行星運行環(huán)境瞬息萬變，各種攝動要素的影響交替變化，具有一定的隨機性，需要加強幾何數(shù)值積分方法與隨機問題的結合研究。

其次，所謂的“二去”，即如何設計發(fā)射近地小行星處置航天器的軌道與實現(xiàn)最優(yōu)控制問題。根據(jù)時變系統(tǒng)的Lagrange擬序結構理論，通過有限時間Lyapunov指數(shù)(FTLE)可以找到時變系統(tǒng)中的分界線，以此可以區(qū)分流場中有不同動力學特征的區(qū)域。然而，關于小行星防御航天器的軌道設計方面，多采用不變流形理論來設計軌道，缺乏采用Lagrange擬序結構為最優(yōu)控制算法提供初始值的研究。另一方面，航天器群制導未利用領航器的優(yōu)化軌道，且基于笛卡爾坐標的重構技術計算量很大，此外，采用DMOC方法計算編隊飛行衛(wèi)星群重構的最優(yōu)開環(huán)控制律，利用Halo軌道，限制了對小行星軌道的應用。

再次，所謂的“三繞”，即如何對近地小行星伴飛航天器進行綜合建模問題。在構建繞小行星的受擾Kepler軌道模型時，主要考慮類地球小行星的J2攝動建立引力攝動模型，無法達到任意小行星所需精度水平，也缺乏考慮相對運動中太陽光壓、其它第三星體引力攝動、隨機擾動等綜合影響，因而這樣的軌道動力學模型精度不高。此外，天體力學N體問題中存在一種規(guī)范對稱性，忽略其規(guī)范自由度可能造成軌道演化動力學在數(shù)值上的不穩(wěn)定性[72]。由于現(xiàn)有模型主要考慮類地球系統(tǒng)的引力場J2項，密切狀態(tài)到平均狀態(tài)轉換方法計算效率和計算精度較低，不適用于任意非地球小行星和其它行星體。有的剔除密切(短周期)頻率的頻率分解法具有計算代價大且需頻繁檢測的缺點，也不適于旋轉速度慢的小行星。即使采用Kalman濾波器，僅考慮低階引力場并不能準確得到小行星中心軌道的密切或平均軌道根數(shù)。

最后，所謂的“四控”，即如何對探測/處置航天器群實現(xiàn)有效的編隊制導與控制，從而對具有潛在威脅的小行星實行干預或清除問題。目前，關于小行星防御航天器群隊形重構和保持研究的多數(shù)技術是簡單地推廣近地編隊飛行的方法和成果，不能滿足集群任務關鍵要求。而且，航天器群制導未利用領航器的優(yōu)化軌道，且基于笛卡爾坐標的重構技術計算量很大，缺乏計算量少、能有效利用ΔV、基于相對軌道根數(shù)的大規(guī)模編隊重構控制方案。當前的一些方法只能考慮一個軌道周期的短時間跨度，動力學的線性化只在極短的時間內是準確的，而頻繁的線性化又會增加計算成本。基于相對軌道根數(shù)的編隊重構控制方案通常局限于兩個航天器，并不能解決大群體的避碰問題。另外，缺乏解決在小行星特有的動態(tài)環(huán)境下航天器群全姿態(tài)最優(yōu)控制方法，忽略航天器姿態(tài)可能會影響任務的執(zhí)行。

綜上所述，目前針對小行星防御問題的軌道預測、偏移與設計、航天器群的編隊最優(yōu)制導與控制的研究還不完善，尚缺乏綜合考慮各種攝動效應的精確的數(shù)值計算方法和合理的小行星動力學模型，以及帶隨機擾動的在軌航天器(群)全位姿動力學模型；同時，圍繞地球以外小行星特有動態(tài)環(huán)境下的航天器群編隊最優(yōu)制導與控制方法仍需進一步發(fā)展。上述問題的解決可為小行星的軌道偏轉、天基監(jiān)測系統(tǒng)的建立、引力拖車等編隊任務提供理論指導和技術支持。鑒于此，采用基于幾何力學和幾何數(shù)值積分建立動態(tài)環(huán)境下的高精度Hamel變分積分子，通過高精度積分子研究長期作用力下小行星軌道的演化規(guī)律，對小行星的軌道偏轉問題進行動力學建模，實現(xiàn)移步易景，實時計算小行星的運動狀態(tài)，做到對小行星偏轉軌道的精準預測，給小行星的災害防御提供精準的動態(tài)數(shù)據(jù)。在此基礎上，構建處置航天器基于Lagrange擬序結構的最優(yōu)控制算法；建立考慮各種攝動的小行星在軌航天器全位姿動力學模型；針對航天器群編隊的構建、保持和重構等任務，研究其場論框架下柔性集群編隊控制理論和算法，為小行星防御的在軌處置提供理論支撐。