基于特征模型的Halo軌道維持控制

張斌，周敬

北京控制工程研究所，北京 100190

近年來，關于平動點任務的研究受到越來越多的關注。本文討論的平動點屬于三體問題范疇，是天體力學領域的基本問題之一。Euler和Lagrange分別在1765年和1772年研究限制性三體問題(Restricted Three Body Problem, RTBP)時，發現RTBP存在5個平衡點，稱之為平動點或拉格朗日點[1]，其中，位于兩個主天體連線上的3個平動點稱為共線平動點。由于平動點獨特的動力學特性，將航天器部署到平動點上時，航天器在引力作用下，與主天體的相對位置保持不變。然而在真實空間環境下，平動點位置無法被精確確定，而相比于平動點，其附近的周期軌道更適用于航天任務。因此實際任務中通常將航天器部署在平動點附近的軌道上。

平動點的動力學內涵十分豐富[2]，在其附近存在著豐富的軌道類型，包括周期軌道(如Lyapunov軌道、Halo軌道等)、擬周期軌道(如Lissajous軌道、擬Halo軌道等)以及不變流形等[3-5]。這些軌道為一些特殊的深空探測任務以及空間低能轉移提供了可能。例如，日地L2點因其獨特的空間位置和物理環境，為空間望遠鏡等深空觀測衛星提供了良好的部署條件[6]；日地L1點位于日地連線之間，是部署太陽活動監測衛星的絕佳地點[7]；地月L2點位于地月連線上月球的背面，選擇合適的軌道尺寸能夠使部署在上面的航天器有效規避月掩，是為地球與月球背面提供中繼通信的理想場所[8]；此外，平動點在太陽系中廣泛存在，這些平動點附近的不變流形相互連接，構成了一系列蜿蜒錯綜的管道，稱為行星際高速公路(InterPlanetary Superhighway, IPS)，利用IPS，可以實現行星際的低能轉移[9]。

平動點軌道在實際任務中的應用起步相對較晚。1978年，美國航空航天局(NASA)成功發射ISEE-3探測器，是平動點軌道第1次在實際任務中得到應用。以ISEE-3的成功實踐為起點，拉開了利用平動點進行深空探測的序幕[10]。在接下來的幾十年里，以NASA和歐洲航天局(ESA)為代表的航天機構，先后成功實施了SOHO、ACE、MAP、Genesis、ARTEMIS等多次平動點任務[11]，在太陽活動監測、深空觀測等領域取得了豐碩的成果。中國先后進行了嫦娥二號、嫦娥五號T1以及“鵲橋”衛星3次平動點任務，其中于2018年成功發射的“鵲橋”通信中繼衛星實現了人類歷史上首次月背中繼通信，為嫦娥四號進行人類首次月背軟著陸探測提供了有力支撐。

雖然共線平動點附近存在著豐富的軌道類型，但由于其本身是不穩定的，運行在其附近軌道上的航天器若不施加控制，將會很快偏離標稱軌道。因此，設計平動點任務時，進行平動點軌道維持研究十分必要。針對平動點軌道維持問題，許多學者進行了研究，并取得了一些成果。Simo等[12]利用Floquet模態理論針對擬Halo軌道和Halo軌道的維持問題進行了研究；Breakwell等[13]最早利用線性二次型調節器(LQR)方法實現了對Halo軌道的維持控制；Howell和Pernicka[14]提出了利用靶點法維持平動點軌道的策略，并在ARTEMIS任務中得到了應用；Lian等[15]利用離散時間滑模控制實現了對真實星歷下地月系平動點軌道的維持控制；Nazari等[16]分別利用時變LQR、退步控制和反饋線性化實現了地月L1附近Halo軌道的維持控制，并對3種方法的仿真結果進行了對比分析；徐明和徐世杰[17]設計了一種線性周期控制策略實現了Halo軌道的維持控制。

基于特征模型的自適應控制方法是由吳宏鑫院士在20世紀80年代初提出的一種從實際應用角度出發的控制方法。經過30多年的發展，不僅在理論上取得了重要進展，同時還成功應用于神舟飛船再入返回控制、交會對接等重大工程中[18]。所謂特征模型，即結合對象動力學特征、環境特征和控制性能要求而不是僅以對象精確動力學分析所建立的模型，其具有以下特點[19]：① 在同樣控制作用下，對象特征模型和實際對象在輸出上是等價的，在穩定情況下，輸出是相等的；② 特征模型的形式和階次除考慮對象特征外，主要取決于控制性能要求；③ 特征模型建立的形式比原對象的動力學方程簡單，易于控制器設計，工程實現容易、方便；④ 特征模型與高階系統的降階模型不同，它是把高階模型有關信息都壓縮到幾個特征參量中，并不丟失信息，一般情況下用慢時變差分方程描述。

由于三體問題相對于二體問題更為復雜，而基于特征模型的控制方法具有將復雜問題簡單化的優勢，因此本文基于特征模型理論，研究了地月系L2點附近Halo軌道的維持控制問題。本文結構安排如下：第1節介紹兩類動力學模型及標稱軌道的獲取；第2節介紹黃金分割控制器和PD控制器的設計過程；第3節介紹仿真結果及結果分析；第4節總結全文，給出結論。

1 動力學模型與標稱軌道

1.1 圓型限制性三體模型

圓型限制性三體問題(Circular Restricted Three Body Problem, CRTBP)是最簡化形式的三體問題，其可以描述為：一個質量無窮小的天體在兩個大質量天體的引力作用下的運動，其中兩個大質量天體(又稱為主天體)圍繞其公共質心做圓形周期運動。在航天應用領域，常涉及的主天體主要是地球-月球或者太陽-地球/月球。

通常在如圖1所示的會合坐標系下對CRTBP進行研究。圖中:m3為航天器;M1和M2為兩個主天體(M1≥M2);坐標原點O為它們的公共質心，Ox由M1指向M2，Oz指向M1和M2運動的角動量方向，xOy平面為M1和M2的運動平面，Oxyz構成右手坐標系；r、r1、r2分別為由O、M1和M2指向m3的矢徑。為了方便計算，對各物理量進行無量綱化處理。歸一化后的質量、長度、時間單位分別為

圖1 會合坐標系Fig.1 Syzygy frame

(1)

式中：等式右邊的各物理量均采用國際單位制，M1和M2分別為兩個主天體的質量；L12為主天體之間的距離；G為萬有引力常量，則可得到歸一化后的主天體M1和M2的質量分別為1-μ和μ，μ為質量參數，形式為

(2)

主天體在會合坐標系中的坐標分別為(-μ,0,0)和(1-μ,0,0)。在地月系中，μ=0.012 150 585 61。

(3)

式中：Ω為有效勢，其形式為

(4)

(5)

1.2 雙圓限制性四體模型

CRTBP是最理想情況下的簡化模型，在實際情況中，航天器在平動點軌道上的運動會受到第四體引力攝動、系統偏心率、太陽光壓等攝動因素的影響。對于地月三體系統，太陽引力攝動是影響航天器運動的主要攝動源，因此在CRTBP基礎上，引入太陽引力攝動，建立雙圓限制性四體模型(BiCircular Restricted Four Body Model, BCRFBM)能夠更加精確地模擬真實力學環境。

雙圓限制性四體坐標系如圖2所示，在BCRFBM假設下，太陽在地月運動平面內繞系統質心做勻速圓周運動。圖中:M4為太陽，坐標軸定義與會合坐標系相同；r、r1、r2、r4分別為由O、M1、M2和M4指向m3的矢徑；d4為O指向太陽的矢徑，θ4為d4與坐標系x軸的夾角。

圖2 雙圓限制性四體坐標系Fig.2 Coordinate system of BCRFBM

在BCRFBM假設下，航天器m3的動力學方程為

(6)

式中：

(7)

1.3 標稱Halo軌道

標稱軌道是Halo軌道維持研究的基礎。Halo軌道是共線平動點附近的一類周期軌道，是水平Lyapunov軌道分岔的結果。由于三體動力學系統具有強非線性特征，無法獲得Halo軌道的精確解析解，因此標稱軌道通常采用數值積分的方式求解。

數值積分的初值，通常通過修正由近似解析解提供的近似初值獲得。Richardson[20]在1980年基于Lindstedt-Poincaré方法推導了Halo軌道的三階近似解析解，其表達式為

(8)

式中：τ1為Halo軌道的線性化相位角，τ1=τ0+ωt，τ0為初始相位角；Ax、Az分別為Halo軌道在平動點旋轉坐標系中的線性化振幅，且Ax和Az之間存在限制性條件：

(9)

式(8)和式(9)中的其他各參數定義見文獻[20]。

利用式(8)獲得的積分初值只是近似初值，直接利用此初值代入動力學方程式(3)中進行數值積分，結果將很快發散。因此，還需對近似初值進行修正，以獲得滿足精度要求的積分初值。

對近似初值的修正可以采用微分修正的方法進行。由于Halo軌道具有關于xOz平面對稱的幾何性質，因此選取Halo軌道與xOz平面的一個交點作為軌道初值，記為

若初值為精確值，則積分半個軌道周期T/2時，軌道第一次穿越xOz平面，且有

可達到修正初值的效果。實際修正中，可固定z0，簡化修正過程，則微分修正量可通過如下形式獲得[21]：

(10)

設置Halo軌道z向振幅為Az=0.016 6，依據微分修正方法獲得了修正后的積分初值，然后將此初值代入動力學方程式(3)，數值積分可得如圖3所示的Halo軌道，其軌道周期在歸一化單位下為T=3.412 2，相當于約14.738 6天，圖中LEM表示地球與月球之間的距離。

圖3 微分修正初值之后的Halo軌道Fig.3 Halo orbit with initial value obtained by differential correction

值得注意的是，由于Halo軌道本身是不穩定的，經過微分修正初值后獲得的Halo軌道也不是完全閉合的軌道，經過多圈積分之后仍會發散。考慮到實際任務需求，僅依靠微分修正后的Halo軌道作為標稱軌道顯然無法完成較長時間周期內的軌道維持，因此還需對微分修正后的軌道做進一步的修正，獲得一條持續時間較長的標稱軌道。

為獲得滿足需求的標稱軌道，本文采用靶點法對Halo軌道進行了進一步的修正。靶點法的基本思想是通過最小化目標函數J，求解使軌道維持在標稱軌道附近所需的最小速度脈沖，其中目標函數J是機動速度和靶點狀態相對標稱軌道偏差的加權和，具有以下形式：

(11)

式中：Q、Rk和Tk為權值矩陣；pk和vk為航天器在第k個靶點處相對標稱軌道的位置和速度偏差，可利用狀態轉移矩陣Φ(tk,tc)獲得

(12)

式中：pc和vc為當前時刻tc的位置和速度偏差。將pk和vk代入式(11)中，令

(13)

即可求得Δvmin。

根據上述原理，以微分修正所獲得的Halo軌道第一圈為基礎，設計三靶點的修正策略，可以得到圖4所示的標稱軌道。修正過程中，修正間隔為0.01T，持續20個周期，圖5給出了速度修正量相對于速度的比例，可以看出除個別情況外，修正量均不超過速度的0.04%，因此可以認為標稱軌道是連續的。

圖4 標稱軌道(20周期)Fig.4 Reference orbit (20 periods)

2 控制器設計

(14)

式中：

其中：u稱為控制加速度。

2.1 速度跟蹤的多變量黃金分割控制器

根據文獻[19]，對多輸入多輸出系統建立特征模型，需輸入輸出同維。由于式(14)中控制量為3維，而為實現Halo軌道維持需要對全狀態(6維)進行跟蹤，因此無法對全狀態建立特征模型，故考慮僅對速度跟蹤特征建模。

在采樣周期Δt滿足一定條件下，速度跟蹤的特征模型可用如下輸出解耦型二階差分方程組描述：

X1(k+1)=F1(k)X1(k)+F2(k)X1(k-1)+G0(k)U1(k)+G1(k)U1(k-1)

(15)

式中：

且有當Δt→0時

式(15)中的特征參量F1(k)、F2(k)、G0(k)和G1(k)由參數估計得到。盡管是多變量，對于式(15)的特征模型仍可按每路獨立進行參數估計[19]。每一回路的特征方程可寫為

xj(k)=f1j(k)xj(k-1)+f2j(k)xj(k-2)+

(16)

式中：

參數估計采用如下最小二乘遞推算法：

(17)

(18)

2.2 位置跟蹤的PD控制器

u2=up+ud

(19)

式中：

(20)

為了便于調節，ud采用邏輯微分，其形式為

(21)

式中：cj、Nj為常數。

至此，可以得到最終的控制器為

u=u1+u2

(22)

由MIMO特征模型式(15)和黃金分割控制器式(18)組成的閉環系統的穩定性，一直是相關研究學者關注的問題，其中文獻[22-23]均給出了閉環系統穩定的充分條件，但該充分條件仍有很大的局限性，因此基于MIMO特征模型的黃金分割控制器的穩定性證明仍有待進一步研究和完善。由于控制器(22)為離散形式，而原系統(14)為連續模型，因此由控制器(22)和原系統組成的閉環系統是一個典型的混雜系統，目前關于混雜系統穩定性證明尚無有效的工具，仍有待進一步研究。基于以上分析，本文中暫不涉及控制器穩定性的證明。

3 仿真結果

為了驗證第2節中設計的控制器的有效性，選取了地月L2點附近的Halo軌道進行了仿真，并與LQR方法下的仿真結果進行了對比。設置軌道z向振幅為Az=0.016 6，相當于約6 391.5 km；周期為T=3.412 2，相當于14.738 6天；仿真總時長為20個周期，約300天。

3.1 CRTBP模型下的仿真結果

針對地月L2點附近的Halo軌道維持控制問題，利用第2節中設計的控制器，首先在CRTBP模型下進行了仿真。考慮到實際任務中，航天器進入Halo軌道時，不可避免地會產生入軌誤差，導致實際軌道偏離標稱軌道，因此，仿真時設置了初始位置誤差和速度誤差，其大小為

約相當于國際單位制下的：

設置系統采樣周期為Δt=0.001，仿真結果如圖6所示。

位置誤差曲線和速度誤差曲線如圖7(a)和7(b)所示。由圖7(a)可以看出，采用黃金分割控制器和PD控制器的情況下，位置和速度誤差均在很短時間內收斂到0；而采用LQR的情況下，速度誤差雖然仍能很快收斂，但位置誤差則表現出明顯的震蕩。

穩態后采用兩種控制器情況下的位置和速度誤差均值如表1所示。由表中數據可以看出，相比于LQR方法，采用本文設計的黃金分割控制器和PD控制器進行軌道維持，位置誤差均值和速度誤差均值均有顯著下降，由此可以說明本文設計的控制器相比于LQR方法具有更高的控制精度。

控制加速度曲線如圖7(c)所示。可以看出，采用黃金分割控制器和PD控制器情況下，控制加速度最大值出現在起始時刻，之后很快趨向于0；而采用LQR方法情況下，控制加速度則一直相對較小。對控制加速度進行時間積分可得到軌道維持過程中所需要的速度增量，如表2所示。

圖6 CRTBP模型下維持結果Fig.6 Station-keeping under CRTBP model

圖7 位置誤差、速度誤差和控制加速度(CRTBP)Fig.7 Curves of position errors, velocity errors and control acceleration (CRTBP)

表1 位置和速度誤差均值(CRTBP)
Table 1Averages of position and velocity errors(CRTBP)

誤差控制器x軸y軸z軸位置誤差均值/m黃金分割+PD10.34597.42340.8269LQR11007265571124速度誤差均值/(m·s-1)黃金分割+PD0.00150.00120.0002LQR0.07010.08610.0051

表2 速度增量(CRTBP)Table 2 Velocity increment( CRTBP)

從表2中數據可以看出，20個周期的仿真過程中，采用本文設計的控制器，所需要的總速度增量為95.513 0 m/s，采用LQR方法所需要的總速度增量為31.670 8 m/s。由于仿真初值加入了入軌誤差，起始階段為了克服初始誤差，需要比較大的控制加速度，同時，由于參數辨識未收斂，也會導致初始控制量較大，采用本文設計的控制器時，第一個周期所需要的速度增量為73.191 7 m/s，相比于LQR方法下的1.858 5 m/s，控制消耗較大。而達到穩態后，采用本文設計的控制器時，每個周期所需要平均速度增量則僅為1.322 4 m/s，相比于第1周期有顯著降低；而采用LQR方法，穩態后平均每個周期所需的速度增量為1.569 9 m/s，相比于第1周期也有所下降，但仍高于采用本文設計的控制器時的所需的速度增量。

由此可知，在CRTBP模型下，相比于LQR方法，利用第2節所設計的控制器進行軌道維持，可實現精度較高的控制效果，且進入穩態后，控制消耗較小。另外，提高入軌精度，減小入軌時的位置和速度誤差，可以顯著降低起始階段的控制消耗。

3.2 BCRFBM模型下的仿真結果

為驗證控制器在太陽引力攝動影響下的有效性，本文基于BCRFBM模型，利用第2節設計的控制器以及LQR方法分別進行了仿真。仿真過程中，仍然考慮了入軌誤差的影響，初始誤差大小與式(23)相同，采樣周期同3.1節，初始時刻太陽相位角取θ40=0。仿真結果如圖8所示。

圖8 BCRFBM模型下維持結果Fig.8 Station-keeping results under BCRFBM model

位置誤差曲線、速度誤差曲線以及控制加速度曲線如圖9所示。可以看出，采用第2節設計的控制器時，位置和速度誤差同CRTBP模型下的仿真結果一樣，均能在很短時間內收斂到0；而采用LQR方法時，位置誤差和速度誤差均呈現明顯的震蕩，不能收斂。

圖9 位置誤差、速度誤差以及控制加速度曲線(BCRFBM)Fig.9 Curves of position errors, velocity errors and control acceleration (BCRFBM)

穩態后采用兩種控制器情況下的位置和速度誤差均值如表3所示。由表中數據可以看出，采用第2節設計的黃金分割控制器和PD控制器進行軌道維持時，相比于初始誤差，位置誤差均值和速度誤差均值仍能維持在較低的水平；而采用LQR方法時，位置和速度誤差均值相比于初始誤差均有較大幅度上升。由此可以再次說明本文設計的控制器在控制精度方面的優勢。

表3 位置和速度誤差均值(BCRFBM)Table 3 Averages of position and velocity errors (BCRFBM)

從控制加速度曲線圖可以看出，與CRTBP模型下的仿真結果類似，采用第2節設計的控制器時，由于存在初始誤差，且參數辨識尚未收斂，控制加速度在起始時刻較大，之后很快趨向于0；而采用LQR方法，控制加速度一直比較平穩。

對控制加速度進行時間積分可得到軌道維持過程中所需要的速度增量，如表4所示。從表中數據可以看出，20個周期的仿真過程中，采用第2節設計的控制器時，所需要的總速度增量為745.024 6 m/s；由于仿真初值加入了入軌誤差，起始階段為了克服初始誤差，同時由于參數辨識未收斂，需要比較大的控制加速度，第1周期所需要的速度增量為130.037 4 m/s；達到穩態后，每個周期所需要平均速度增量則僅為33.166 3 m/s，相比于第1周期有顯著降低。而采用LQR方法時，所需的總速度增量、第1周期所需的速度增量以及穩態后每個周期所需要的平均速度增量，分別為610.551 6、30.294 2和30.544 7 m/s，相比于使用本文設計的控制器情況下均有所下降。

表4 速度增量(BCRFBM)Table 4 Velocity increment(BCRFBM)

與CRTBP模型下仿真結果對比可以看出，采用第2節設計的控制器時，在軌道維持過程中引入太陽引力攝動，所需要的速度增量提升了一個數量級；x軸和y軸的位置跟蹤誤差同樣也提升了一個數量級，分析原因，主要在于太陽引力攝動主要作用于航天器x軸和y軸運動，對z軸運動則影響不大，但考慮到Halo軌道本身尺度較大，這樣的位置誤差仍是可以接受的；z軸的位置跟蹤誤差以及三軸的速度跟蹤誤差則與CRTBP下的仿真結果相差不大。同時，與CRTBP模型下的仿真結果類似，提高入軌精度，起始階段的控制消耗有望進一步降低。而采用LQR方法時，雖然控制消耗有所下降，但其在消除跟蹤誤差方面表現很差，因此綜合來看，在達到穩態后，本文設計的控制器仍具有優勢。

從仿真位置和速度誤差來看，本文所設計的控制方法是一種跟蹤精度較高的策略。在Halo軌道的維持控制中，控制精度與控制消耗互相制約，更高的控制精度通常也意味著更高的控制消耗。此外，根據文獻[24]中的研究結論，標稱軌道的精確度與控制消耗密切相關。標稱軌道越精確，維持控制所需要的控制消耗則越低。本文所采用的基準Halo軌道是在CRTBP模型下修正的結果，對于考慮太陽引力的BCRFBM模型下的維持控制來說，精確度不夠。若能采用更精確的標稱軌道，維持控制所需的總的速度增量有望進一步降低。

4 結 論

本文首先基于Richardson三階近似解析解、微分修正及打靶法獲得了CRTBP模型下的地月L2點附近的標稱Halo軌道，然后設計了基于特征模型理論速度跟蹤的黃金分割控制器以及位置跟蹤的PD控制器，最后分別在CRTBP模型和BCRFBM下進行了地月L2點Halo軌道維持控制仿真，并與LQR方法下的仿真結果進行對比，驗證了控制器的有效性。通過分析仿真結果，可以得到以下結論：

1) 在CRTBP模型下，仿真結果跟蹤精度較高，控制消耗很小。

2) 在BCRFBM模型下，位置跟蹤精度有所下降，但相比于Halo軌道的尺度，仍然可以接受；但由于太陽引力作用以及標稱軌道精度影響，控制消耗顯著增大。

3) 提高入軌精度，可以顯著降低起始階段的控制消耗。