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

黃英志，泮斌峰，唐 碩

(西北工業(yè)大學(xué)航天學(xué)院，陜西西安 710072)

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

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

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

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

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

圖1 Lambert問題的幾何構(gòu)型Fig.1 Geometry of the Lambert problem

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

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

從而偏心率的大小可寫為

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

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

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

則由軌道力學(xué)基本關(guān)系可得軌道的半長軸為

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

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

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

由偏近點角與真近點角的關(guān)系可得到終端點處的偏近點角為

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

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

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

在得到以上公式及性質(zhì)以后，即可數(shù)值求解出給定Δt值后的對應(yīng)所有可能的解。

圖2 歸一條件下轉(zhuǎn)移時間隨eT的變化關(guān)系Fig.2 Transfer time as a function of eT

2 軌道性質(zhì)與橫向偏心率關(guān)系

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

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

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

令

式(16)化簡為

二次方程(18)的解為

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

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

圖3 長短徑分界點與轉(zhuǎn)移角的關(guān)系Fig.3 Locations of dividing point as a function of transfer angle

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

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

3 約束函數(shù)的特性分析

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

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

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

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

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

對于f函數(shù)，其極值點位于對eT導(dǎo)數(shù)為零處，對eT求導(dǎo)得:

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

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

代入式(17)并整理得

同理，對于F函數(shù)，其對eT求導(dǎo)得

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

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

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

式(31)等效為

整理為關(guān)于eT的二次方程:

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

圖5 約束函數(shù)極值點與轉(zhuǎn)移角的關(guān)系曲線Fig.5 Extremum points of constrain functions as a function of transfer angle

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

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

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

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

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

整理可得:

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

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

4 解算過程

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

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

如圖6所示，在式(23)的限制范圍內(nèi)，N－eT曲線與整數(shù)值的橫軸交點即為可行解，則在給定Lambert問題所有數(shù)學(xué)解中，可行解的運行圈數(shù)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」)。

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

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

圖6 軌道運行圈數(shù)Fig.6 Revolution number

5 算例與結(jié)果分析

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

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

現(xiàn)在考慮施加約束，令Rp=RE+350 km，Ra=RE+20 000 km，圖7給出f函數(shù)與F函數(shù)隨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)內(nèi)。軌道長、短徑分界點為=0.254 7。約束邊界處軌道圈數(shù)滿足「N(eTL)? =3，?N(eTU)」=5，故可行軌道滿足 N=3，4，5。分別對3條軌道以eTL= －0.532 9為初始猜測進行迭代計算，結(jié)果分別收斂至eT3，eT4和eT5，迭代次數(shù)分別為 4，5，5 次。至此，得到了全部11個解中滿足工程約束的3條可行軌道，且3條可行的軌道皆為短徑軌道。所求出的軌道解如圖8中上三角符號所示。

圖7 兩種約束函數(shù)隨eT的變化關(guān)系Fig.7 Constrain functions versus eT

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

6 結(jié) 論

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

[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.