基于偏心率向量含約束多圈Lambert問題

黃英志，泮斌峰，唐 碩

(西北工業大學航天學院，陜西西安 710072)

Lambert問題作為二體軌道的邊值問題，通常定義為[1]:“二體問題中，給定兩個終端位置向量和轉移時間，求轉移軌道”。作為天體力學中的經典問題，又由于其在空間技術應用方面和數學上的重要性，從高斯時期起，Lambert問題吸引了大量學者的關注，從而有多種求解方法被提出來。

Lambert問題最早的求解方法可追溯到高斯，他利用圓錐曲線的幾何性質解決了此問題。Lancaster[2]第一次提出使用無量綱量 x=作為迭代變量數值求解Lambert問題，統一了橢圓、拋物線和雙曲線全部轉移軌道。Battin[3]沿用 Lancaster提出的積分變量，利用二點邊值問題的不變量(Invariance)拓展改進了Gauss的方法，并彌補了Gauss方法在轉移角為180°時的奇異問題。Sun[4]同樣以 x為變量，詳細研究了軌道轉移時間的極值特性，給出了多圈Lambert問題中解與最大飛行圈數的關系。在以其他積分變量為基礎的研究中，Nelson[5]使用航跡角γ為迭代變量對問題進行求解，然而他未能給出轉移時間相對于γ的解析導數，只使用二分法進行搜索;Jaemyung[6]同樣使用γ為迭代變量，但給出了解析導數信息，使得數值解算更加容易收斂。特別值得注意的工作來自 Avanzini[7]，他首次提出使用偏心率向量垂直于轉移弦線的分量大小eT為迭代變量，建立了以其為變量的描述體系，并給出了轉移時間相對于eT的解析導數，便捷而直觀地解決了 Lambert問題;Quan He[8]將Avanzini的工作推廣到了多圈Lambert問題。以上對Lambert問題的解均是從純數學角度描述而并未考慮實際工程中的約束情況，Zhang[9]首次提出了將近地點與遠地點的約束加入Lambert問題的求解中，從而從多解中篩選出滿足工程約束的解。然而，Zhang 沿用 Prussing[10]的方法，使用軌道半長軸a為迭代變量。由于a并非轉移時間的單值函數，故而使得約束函數的描述復雜且難以分析其分類討論繁雜冗長。本文使用橫向偏心率向量大小eT為積分變量，詳細討論了以其為變量的軌道參數特性，給出了約束函數的特性分析，求解了引入軌道約束的多圈Lambert問題。從求解過程可看出，與Zhang提出的方法相比，本文提出的方法解算次數更少，求解過程更為直接且討論過程更為簡捷。

1 橫向偏心率向量描述方法

經典Lambert問題如圖1所示，其中F表示重力場的中心，P1和P2為初始與終端位置，其分別對應以F為中心的終端向量r1和r2;Δθ表示兩終端向量間的夾角，即轉移角;c表示連接P1與P2兩點的弦長，e為偏心律向量。

根據文獻[1]的推導，當給定了Lambert問題的幾何條件后，其偏心率向量e唯一決定了軌道的形狀，并且在P1P2的投影是一定值。這也就是說，偏心率向量可以分解成為垂直于P1P2的“橫向”可變分量eT和平行于P1P2的定常分量eF。這一性質使得我們可以使用偏心率向量的橫向分量eT來唯一描述轉移軌道的性質。

圖1 Lambert問題的幾何構型Fig.1 Geometry of the Lambert problem

縱向分量可表示為[1]:

其中ic表示沿P1P2方向的單位向量。橫向分量eT的大小用eT表示，定義與橫向單位向量ip同向時為正，ip與P1P2垂直。橫向分量eT與縱向分量eF的和即為偏心率向量:

從而偏心率的大小可寫為

當橫向分量eT等于零時，偏心率向量退化為eF且與弦線P1P2平行，此時的轉移軌道為Battin所定義的Lambert問題中的基本橢圓(Fundamental Ellipse)，該橢圓的長半軸與弦線P1P2平行，如圖1實線橢圓所示?；緳E圓的基本參數可以作為轉移軌道參數的中間量，其偏心率、半長軸以及半通徑由以下公式給出

式中角標F表示其為基本橢圓軌道下的量。

在基本橢圓軌道的基礎上，我們以半通徑為中間量建立一般轉移軌道與基本橢圓的聯系，從而以 eT為變量得到一般橢圓的參數。根據Avanzini[6]給出的公式有

則由軌道力學基本關系可得軌道的半長軸為

在軌道平面內定義r1及其垂線方向為參考直角坐標系兩軸方向，ωc為ic與r1的夾角，則ic可表示為(cos ωc，sin ωc)T，ip可表示為(－ sin ωc，cos ωc)T，從而對一般轉移軌道，其偏心率向量可表為

則e相對參考向量r1的角度ω可表為

式中的tan－1(x，y)為四象限反正切函數。根據上式由幾何關系可得到P1和P2兩點的真近點角為

由偏近點角與真近點角的關系可得到終端點處的偏近點角為

由式(1)－(11)可見，在給定了Lambert問題的幾何參數r1，r2以及Δθ以后，轉移軌道的基本參數如半長軸a，半通徑p，終端點的真近點角θ1θ2以及偏近點角E1E2都可表為橫向離心率eT的函數，從而建立了以eT為自變量的參數體系。Lambert問題的軌道轉移時間是以上參數的函數，從而成為了eT的單值函數。軌道轉移時間Δt可表示為[7]

式中N為轉移軌道的圈數。式(12)表明，當多圈Lambert問題幾何參數r1和r2以及圈數N給定以后，轉移時間Δt就成為了eT的函數，且由于函數曲線特性良好，最簡單的牛頓迭代法就足以很快收斂到給定Δt的解。牛頓迭代法所需的導數信息可由式(12)對eT直接求導得到，Quan He[7]給出了該式導數求解過程，本文不再贅述。值得注意的是，對于多圈Lambert問題，我們只考慮橢圓軌道，故eT應限制在一定范圍內。事實上，由于橢圓軌道偏心率滿足0＜e＜1，于是結合式(3)有:

由于eF的取值范圍被限制在(－1，1)內，故而數值迭代變得十分方便。圖2給出了多圈Lambert問題轉移時間Δt相對于eT的曲線，圈數N取0－5，其幾何參數采用歸一化參數，取r1=1，r2=2，Δθ=120°，μ =1。由圖2 可見，Δt是 eT的單值函數，且其導數在eT上下界兩點以外的地方皆為有界的。當N=0時，Δt隨eT單調遞增;當N ＞0時，同一Δt對應兩個eT的解，Δt在eT的上下界處趨于無窮大，并存在一個導數為零的極小值。正如Prussing[10]得出的結論，在給定Δt以后，多圈Lambert問題一共有2Nmax+1個解，其中Nmax是最大允許的轉移軌道圈數。

在得到以上公式及性質以后，即可數值求解出給定Δt值后的對應所有可能的解。

圖2 歸一條件下轉移時間隨eT的變化關系Fig.2 Transfer time as a function of eT

2 軌道性質與橫向偏心率關系

在以半長軸a為自變量描述Lambert問題時，對同一軌道轉移時間，有兩個不同的a值的解，分別對應所謂的長徑(long－path)和短徑(short－path)軌道，造成了分析的困難。由于Lambert問題中軌道與eT是一一對應的，故而同一個eT有兩個a值相對應。事實上，確定了軌道的半長軸以后，了解其屬于長徑軌道還是短徑軌道對軌道分析具有重要意義，因此本節討論半長軸與eT的關系。由式(6)和(7)知

軌道長徑與短徑的分界點位于a對eT的導數為零的點，對(14)式求導有:

對于橢圓軌道有e＜1，故而極值點有

令

式(16)化簡為

二次方程(18)的解為

注意到eT須滿足式(13)，故知其必處于(－1，1)內。對式(19)分析可知，當θ∈(0，π)時，eT1單調遞增且趨于正無窮(由于為正且趨于零);θ∈ (π，2π)時，eT2單調遞增且趨于正無窮，故其在對應的區間上都不滿足式(13)的約束，如圖3所示。由以上分析可知半長軸的極小值點為

如圖1所示，當eT為負時，e所指向的近地點位于轉移弧線上，轉移軌道屬于短徑，故當eT小于e*T時，轉移軌道為短徑，反之為長徑。圖4給出了 r1=1，r2=2，Δθ=120°時a與eT的關系，其中虛線表示短徑，實線表示長徑。由此關系分析圖2，同樣以虛實線表示了軌道的長徑短徑性質，可見其屬性由式(20)決定，即e*T完全由Lambert問題的幾何參數決定，故應注意在圖2中長短徑的分界點并不位于每條Δt－eT曲線的極值點，亦即同一N值Lambert問題的兩個解并不一定分別對應長徑和短徑軌道。

圖3 長短徑分界點與轉移角的關系Fig.3 Locations of dividing point as a function of transfer angle

圖4 軌道半長軸與eT的關系曲線Fig.4 Semimajor axis as a function of eT

在確定了e*T的值以后就可以很快地判斷所求出的Lambert解究竟屬于長徑軌道還是短徑軌道。

3 約束函數的特性分析

在采用以上方法可以求解出所有滿足Δt條件的Lambert轉移軌道并確定其軌道性質，然而這些軌道并不全都滿足物理和工程上的約束。Zhang[8]提出了兩種對軌道的約束:①軌道的近地點高度應大于某定值，以避免落入近地大氣中;②軌道的遠地點高度應小于某定值，以避免衛星可用觀測時間過短。然而，Zhang在文中使用半長軸a為變量進行分析，a的非單值性使得對兩種約束建立的函數無法精確分析，從而使得討論變得十分復雜。下面的討論使用eT為變量對兩種約束進行分析，并極大簡化討論的過程。

為了使得軌道滿足近地點及遠地點的約束，要求以下不等式成立:

其中Rp表示近地點高度下限，Ra表示遠地點高度的上限，其均為定值。由式(21)的約束條件，定義約束函數f函數和F函數，分別要求滿足

由式(22)的兩個不等式限制可以進一步縮小eT的范圍，從而減少迭代求解的次數，并排除不滿足實際需求的軌道。通過觀察f函數和F函數相對于eT的關系曲線形式，可見f函數隨eT變化為先增后減，F函數則為先減后增，二者分別存在一個極大值及一個極小值點。如果能證明兩個約束函數分別滿足這樣的增減關系，那么在式(22)限制下的eT取值范圍將滿足十分簡潔的形式，即:

其中eTp1，eTp2及eTa1，eTa2分別是f－eT曲線和F －eT曲線與y=0的交點，分別位于各自極值點左右兩邊。下文先證明這兩個約束函數在eT的范圍(－1，1)內分別只有一個極值點，并確定極值點的位置;然后給出式(23)中eT范圍的確定方法。

對于f函數，其極值點位于對eT導數為零處，對eT求導得:

令式(24)左端為零，則有等式右端括號內為零:

代入(3)，式(6)和(7)得

代入式(17)并整理得

同理，對于F函數，其對eT求導得

取式(28)右端括號內為零有:

代入式(3)，(6)，(7)及式(17)并整理得到:

綜合式(27)及式(30)，我們得到兩種約束函數f函數和F函數的極值條件:

式(31)等效為

整理為關于eT的二次方程:

為了確定根的分布情況，取eTr1和eTr2分別表示式(34)右端取“+”號和“－”的根，定義輔助函數來判斷所得根為哪個約束函數的極值，由式(31)可知，G=0時說明其為f函數的極值，G≠0時說明其為F函數的極值。圖5表示r1=1，r2=2時，eTr1和eTr2隨夾角Δθ變化的曲線以及對應的G函數值，其中實線表示方程的根值，虛線表示輔助函數G的值。

圖5 約束函數極值點與轉移角的關系曲線Fig.5 Extremum points of constrain functions as a function of transfer angle

對圖5分析可見，當θ∈(0，π)時，對eTr1有0 ＜ G ＜ 2，故eTr1為f函數的極值;當θ∈(θ*，π)時，對eTr2有G=0，故eTr2為F函數的極值點;當θ∈(0，θ*)時，對eTr2有G ＞2，說明此根無效應舍去。其他區間的分析類似，方程根分布情況綜合為表1所示。需注意的是θ*的值由式(34)右端的分母為零對應，即由解得。

表1 約束函數極值點分布Tab.1 Extremum points of constrain functions

綜上可得，在區間eT∈(－1，1)上可能有一個根或沒有根必有一個根，又由于f與F皆為連續函數，故其在討論區間內最多只有一次增減性的變化，式(23)成立。

在確定了兩約束函數的增減性質后，式(23)的求解本身并不困難。對于f函數，令f=0有

代入式(3)，(6)和(7)得:

整理可得:

由式(37)可見f=0的解是一個一元二次方程的兩個根，其解的分布情況符合以上對f函數增減性的分析。與上同理易求解出F=0的兩個解，同樣為一個二次方程的兩個根。值得注意的是，式(37)只是f=0的必要條件而非充分條件，對于超出eT取值范圍式(13)的根應去除而以相應范圍邊界代之。

將兩種約束條件下eT的取值范圍式(23)施加到式(12)的數值解算過程，即可得到滿足實際約束條件的轉移軌道。

4 解算過程

在確定了eT的范圍后，首先可確定可行解軌道運行的圈數。在Lambert問題中給定轉移時間后，轉移圈數成為eT的函數，有:

假設N(eT)在eTm處達到極大值，即

如圖6所示，在式(23)的限制范圍內，N－eT曲線與整數值的橫軸交點即為可行解，則在給定Lambert問題所有數學解中，可行解的運行圈數N滿足:

1)若eTm＞eTU，則

「N(eTL)?≤ N ≤?N(eTU)」;

2)若eTm＜eTL，則

「N(eTU)?≤ N ≤?N(eTL)」;

3)若eTL＜eTm＜eTU，則

N ∈ (min(「N(eTL)?，「N(eTU)?)，?Nmax」)。

應注意情況3)中相同的運行圈數N可能對應兩條軌道。

在求得滿足約束的可行圈數N之后，以求得的約束邊界eTL和eTU為迭代初始變量，經迭代計算即可求出所有可行解。須注意的是，通過這種迭代策略，有可能迭代到約束邊界以外的解，最極端的情況可能求出2n個解(n為可行解個數)，須用式(23)將邊界外的解濾除。在實際的解算過程中，迭代通常會順利收斂得到所求的n個解。

圖6 軌道運行圈數Fig.6 Revolution number

5 算例與結果分析

取終端向量r1=RE+1 000 km，r2=1.5 r1，轉移角為θ=2π/3，地球引力常數為μ=3.986×1014。圖8為轉移時間相對于eT的變化曲線。給定轉移時間Δt=5×104s，首先由牛頓迭代法求解出所有11個可能的eT解(Nmax=5)作為驗證的參考，如表2所示。

表2 所有多圈問題的解Tab.2 All solutions for multiple － revolution problem

現在考慮施加約束，令Rp=RE+350 km，Ra=RE+20 000 km，圖7給出f函數與F函數隨eT的變化曲線。由式(37)解出eTp1=－0.765 5，eTp2=0.083 5，eTa1= － 0.532 9，eTa2=0.762 2，由式(38)得滿足約束條件的軌道必須位于范圍eT∈(－0.532 9，0.083 5)內。軌道長、短徑分界點為=0.254 7。約束邊界處軌道圈數滿足「N(eTL)? =3，?N(eTU)」=5，故可行軌道滿足 N=3，4，5。分別對3條軌道以eTL= －0.532 9為初始猜測進行迭代計算，結果分別收斂至eT3，eT4和eT5，迭代次數分別為 4，5，5 次。至此，得到了全部11個解中滿足工程約束的3條可行軌道，且3條可行的軌道皆為短徑軌道。所求出的軌道解如圖8中上三角符號所示。

圖7 兩種約束函數隨eT的變化關系Fig.7 Constrain functions versus eT

圖8 轉移時間對比eT及可行解Fig.8 Transfer time versus eTand feasible solutions

6 結 論

本文研究了考慮工程約束的多圈Lambert問題，建立了使用eT描述Lambert問題的整個參數體系，證明了兩約束函數在eT的區間內最多只有一次增減性的變化，從而可以分別使用其與上下限的兩個交點作為eT的界限，縮小搜索范圍。使用約束邊界作為迭代計算的初始猜測，可直接解算出滿足約束條件的解，排除不可行解。與之前使用半長軸a作為積分變量分析約束函數相比，本文采用的方法所需解算次數更少，且更加簡單有效，避免了大量的討論過程。數值算例驗證了本文所述方法的簡便性與有效性。

[1]BATTIN R.An Introduction to the Mathematics and Methods of Astrodynamics[M].Revised Edition ed.Reston，VA:AIAA Education Series，AIAA，1999:295－297.

[2]LANCASTER R，BLANCHARD C.A Unified Form of Lambert's Theorem [R].NASA Technical Note，1969D－5368.

[3]BATTIN R H，VAUGHAN R M.An elegant lambert algorithm [J].Journal of Guidance，Control，and Dynamics，1984，7(6):662－670.

[4]SUN F T.On The Minimal Time Trajectory and Multiple Solutions of Lambert's Problem [R].Provincetown:AAS/AIAA Astrodynamics Specialists Conference，1979.

[5]NELSON S L，ZARCHAN P.Alternative approach to the solution of lambert's problem [J].Journal of Guidance，Control，and Dynamics，1992，15(4):1003－1009.

[6]AHN J，LEE S.Lambert algorithm using analytic gradients [J].Journal of Guidance，Control，and Dynamics，2013，36:1751－1761.

[7]AVANZINI G.A simple lambert algorithm [J].Journal of Guidance，Control，and Dynamics，2008，31(6):1587－1594.

[8]HE Q，LI J，HAN C.Multiple－Revolution solutions of the transverse－eccentricity－based lambert problem[J].Journal of Guidance，Control，and Dynamics，2010，33(1):265－268.

[9]ZHANG G，MORTARI D.Constrained multiple－revolution lambert's problem[J].Journal of Guidance，Control，and Dynamics，2010，33(6):1776－1784.

[10]PRUSSING J E.A class of optimal two－impulse rendezvous using multiple－revolution lambert solutions[J].The Journal of the Astronautical Sciences，2000，48(2 － 3):131－148.