劉柔妮，陳 杰，孔祥龍，周世宏

(1.上海衛(wèi)星工程研究所，上海 201109； 2.上海航天技術(shù)研究院，上海 201109)

0 引言

航天器相對(duì)運(yùn)動(dòng)建模是編隊(duì)飛行等航天領(lǐng)域的研究熱點(diǎn)，研究空間自然條件下周期性相對(duì)繞飛軌道對(duì)長(zhǎng)期編隊(duì)飛行航天器有十分重要的意義。通過找到不消耗能量的自然繞飛軌道，可節(jié)省編隊(duì)航天器系統(tǒng)的能量損耗，延長(zhǎng)系統(tǒng)壽命。文獻(xiàn)[1]利用Hill方程和Clohessy-Willshire方程(簡(jiǎn)稱HCW方程)進(jìn)行編隊(duì)衛(wèi)星構(gòu)型設(shè)計(jì)，在圓參考軌道、無攝動(dòng)和線性化的假設(shè)條件下，找到了周期性相對(duì)運(yùn)動(dòng)軌道存在的初始條件。LAWDEN[2]研究了橢圓參考軌道的相對(duì)運(yùn)動(dòng)。TCHAUNER等[3]以真近點(diǎn)角為獨(dú)立變量，通過限制參考基準(zhǔn)軌道偏心率避免了Lawden方程存在的奇點(diǎn)問題。TILLERSON等[4]研究了橢圓軌道中編隊(duì)飛行衛(wèi)星的控制問題。ALFRIEND等[5]率先采用相對(duì)軌道要素描述相對(duì)運(yùn)動(dòng)。肖業(yè)倫等[6]也采用相對(duì)軌道要素描述相對(duì)運(yùn)動(dòng)。之后，韓潮等[7-8]在該概念的基礎(chǔ)上進(jìn)行了一定的研究工作。但以上研究成果均存在不足：Hill方程和Clohessy-Willshire方程使用條件苛刻，需要圓參考軌道、線性化的假設(shè)前提，工程應(yīng)用受到限制；真近點(diǎn)角、相對(duì)軌道要素等描述方式較難考慮周期性相對(duì)運(yùn)動(dòng)的初始條件，求解難度較大。

本文采用哈密爾頓-雅可比(HJ)方程和正則攝動(dòng)理論，引入適當(dāng)?shù)恼齽t變量描述相對(duì)軌道運(yùn)動(dòng)，建立適用于橢圓參考軌道，便于考慮各種攝動(dòng)的相對(duì)運(yùn)動(dòng)模型。經(jīng)推導(dǎo)發(fā)現(xiàn)，在各種攝動(dòng)下相對(duì)運(yùn)動(dòng)的哈密爾頓方程可分解為線性項(xiàng)和攝動(dòng)項(xiàng)，而相對(duì)運(yùn)動(dòng)在任意攝動(dòng)情況下的解可通過求解包含正則變量的微分方程和修正線性項(xiàng)得到，且模型滿足攝動(dòng)疊加原理，即如已得到一定攝動(dòng)下的相對(duì)運(yùn)動(dòng)狀態(tài)，則其他攝動(dòng)下的相對(duì)運(yùn)動(dòng)狀態(tài)只需通過攝動(dòng)項(xiàng)的疊加即可得到。得出周期性相對(duì)運(yùn)動(dòng)條件后，用時(shí)域配點(diǎn)法和列文伯格-馬夸爾特法(LM)對(duì)該條件的初值求解，可得到不消耗任何燃料的周期性繞飛軌道。

1 橢圓參考軌道相對(duì)運(yùn)動(dòng)建模

圖1 旋轉(zhuǎn)Euler-Hill參考坐標(biāo)系示意圖Fig.1 Schematic diagram of rotating Euler-Hill reference coordinate system

建立拉格朗日函數(shù)[9]，設(shè)主航天器在任意橢圓軌道上運(yùn)動(dòng)，考慮活動(dòng)坐標(biāo)系之間的相對(duì)導(dǎo)數(shù)關(guān)系，可得從航天器的速度

(1)

(2)

動(dòng)能表示為

(3)

勢(shì)能Ep(先只考慮正常引力勢(shì))為

Ep= -μ/‖r2‖=-μ/‖r1+ρ‖=

(4)

式中：μ為地球引力勢(shì)；Pk(cosα)為勒讓德多項(xiàng)式,其中,

(5)

從而得拉格朗日函數(shù)

(6)

只取正常引力勢(shì)中的第1，2階可得

(7)

哈密爾頓系統(tǒng)的正則動(dòng)量(px,py,pz)為

(8)

(9)

由式(9)可見，哈密爾頓函數(shù)不顯含時(shí)間t，則根據(jù)哈密爾頓動(dòng)力學(xué)原理可知，HJ方程有如下形式，即

(10)

考慮由(q,p)→(β,α)的正則變換，其中，q,p,β,α為正則變量。S為雅可比完全積分，其形式一般為

S=W(q1,q2,…,qf;α1,α2,…,αf)-

E(α1,α2,…,αf)t

(11)

式中：W(q1,q2,…,qf;α1,α2,…,αf)為哈密頓特征函數(shù)，(q1,q2,…,qf)為廣義坐標(biāo)；E(α1,α2,…,αf)為不含時(shí)間t的函數(shù)，(α1,α2,…,αf)為積分常數(shù)；t為時(shí)間；f為系統(tǒng)的自由度。則正則變換表示為

(12)

系統(tǒng)自由度f現(xiàn)為3，令式(11)中E(α1,α2,…,αf)=-α′1，則可得雅可比完全積分S為

(13)

從而可將HJ方程表示為W(x,y,z)的偏微分方程，即

(14)

將z方向單獨(dú)分離出來，W(x,y,z)=W′(x,y)+W3(z)，則HJ方程變?yōu)?/p>

(15)

(16)

將式(15)積分得

(17)

再令

(18)

同時(shí)令

(19)

由式(19)可得

(20)

(21)

積分得

W1(x)=

(22)

因此，可將雅可比完全積分S表示為

(23)

(24)

(25)

(26)

再根據(jù)正則變換式(12)可得βi，經(jīng)過一系列化簡(jiǎn)和小量忽略可得表達(dá)式

(27)

(28)

再聯(lián)立式(24)，(26)，(27)，(28)，解出x(t)，y(t)關(guān)于α1，β1，α3，β3的表達(dá)式，經(jīng)化簡(jiǎn)后可得

(29)

(30)

由此可得x，y方向上的運(yùn)動(dòng)形式與正則變量α1，β1，α3，β3的關(guān)系。如已知α1，β1，α3，β3的函數(shù)表達(dá)式，就可得到x(t)和y(t)的表達(dá)式，從而可知從航天器相對(duì)于主航天器在x-y平面內(nèi)的運(yùn)動(dòng)模式。

由上文推導(dǎo)過程可知，z方向的運(yùn)動(dòng)和x，y方向是完全解耦的?？上妊芯縳-y平面內(nèi)的運(yùn)動(dòng)特性，再通過分別控制x-y平面的運(yùn)動(dòng)和z方向的運(yùn)動(dòng)得到三維空間的運(yùn)動(dòng)模式。因此，本文主要針對(duì)x-y平面內(nèi)的運(yùn)動(dòng)特性進(jìn)行分析，未專門對(duì)z方向上的運(yùn)動(dòng)進(jìn)行研究。

以上推導(dǎo)都是基于正常引力勢(shì)的第1，2階項(xiàng)得出，將此時(shí)得到的哈密爾頓函數(shù)(9)記為H(0)，當(dāng)考慮引力勢(shì)的更高階項(xiàng)或J2項(xiàng)時(shí)，將哈密爾頓函數(shù)記作H=H(0)+H(1)+H(2)+…+H(n)。其中，H(1)，H(2)，…，H(n)是考慮攝動(dòng)或非線性引力勢(shì)附加的哈密爾頓函數(shù)ΔH，其被稱為攝動(dòng)哈密爾頓函數(shù)，即ΔH=H(1)+H(2)+…+H(n)。

當(dāng)考慮攝動(dòng)時(shí)，由哈密爾頓理論可能較難得出相對(duì)運(yùn)動(dòng)的顯示解析表達(dá)式，但無攝情況下的運(yùn)動(dòng)模型已精確得到，即式(29)，(30)。因此，根據(jù)正則攝動(dòng)理論，當(dāng)ΔH很小時(shí)，可對(duì)無攝條件下得到的動(dòng)力學(xué)模型作適當(dāng)修正，從而得到攝動(dòng)下動(dòng)力學(xué)模型的近似解。

將不考慮攝動(dòng)情況下精確的哈密爾頓函數(shù)記作H0(p,q,t)，則可將雅可比完全積分S(q,α,t)作為廣義函數(shù)，得到從(q,p)→(β,α)的正則變換，且變換后的新的哈密爾頓函數(shù)為零，則(β,α)是一組與時(shí)間無關(guān)的正則變量，由此可以得出相應(yīng)的在無攝情況下關(guān)于模型的精確解。當(dāng)考慮攝動(dòng)時(shí)，將哈密爾頓函數(shù)寫成

H(p,q,t)=H0(p,q,t)+ΔH(p,q,t)

(31)

考慮到對(duì)于給定坐標(biāo)變換的正則性與哈密爾頓函數(shù)的具體形式無關(guān)，因此，由S(q,α,t)作為廣義函數(shù)而產(chǎn)生的正則變換(q,p)→(β,α)在攝動(dòng)條件下還是一組正則變換，只是在攝動(dòng)下得到的新哈密爾頓函數(shù)將不再為零，新的正則變量(β,α)也不再是與時(shí)間無關(guān)的變量。新哈密爾頓函數(shù)的形式為

(32)

因此，正則變量滿足關(guān)系式

(33)

(34)

求解上述2f個(gè)方程(f為自由度，本文有4個(gè)自由度)，得到αi，βi關(guān)于時(shí)間的函數(shù)表達(dá)式，則由(q,p)→(β,α)的變換，可得相應(yīng)的qi，pi關(guān)于時(shí)間的表達(dá)式，從而得到動(dòng)力學(xué)模型。

將上述理論運(yùn)用到考慮非線性引力勢(shì)第3，4階，以及J2項(xiàng)的建模中。非線性引力勢(shì)第3階的攝動(dòng)哈密爾頓函數(shù)可表示為

(35)

非線性引力勢(shì)第4階的攝動(dòng)哈密爾頓函數(shù)可表示為

24μx2y2+24μx2z2-6μy2z2)

(36)

考慮J2攝動(dòng)項(xiàng)的攝動(dòng)哈密爾頓函數(shù)為

(37)

式中：Z=(r+x)sinisinθ+ysinicosθ+zcosi；引力常數(shù)μ=3.960 05×105km3/s2；地球赤道平均半徑Re=6 378.140 km；帶諧系數(shù)J2=0.001 082 63；i為參考軌道傾角；θ為真近點(diǎn)角。

將上述攝動(dòng)哈密爾頓函數(shù)代入求解后可得αi(t)，βi(t)關(guān)于時(shí)間t的表達(dá)式，再代入(q,p)→(β,α)變換關(guān)系式中，即式(29)，(30)中，可得x(t)和y(t)在攝動(dòng)情況下的表達(dá)式，即攝動(dòng)情況下的模型為

(38)

y(t)=-3α3(t)(t+β1(t))+β3(t)+

(39)

(40)

綜上，求解攝動(dòng)情況下的動(dòng)力學(xué)模型只需求解一組非線性方程組，即

(41)

(42)

則可得動(dòng)力學(xué)方程(38)，(39)。

2 周期性相對(duì)運(yùn)動(dòng)條件

橢圓基準(zhǔn)軌道的動(dòng)力學(xué)模型為式(38)，(39)，經(jīng)分析不難發(fā)現(xiàn)，獲得周期為T的周期性相對(duì)運(yùn)動(dòng)的條件為：?t滿足關(guān)系式

α1(t)=α1(t+T)

(43)

α3(t)=α3(t+T)

(44)

f1(t)=β3(t)-3α3(t)(β1(t)+t),

f1(t)=f1(t+T)

(45)

即

α1(t)=α1(t+T)

(46)

α3(t)=α3(t+T)

(47)

β1(t)=β1(t+T)

(48)

(49)

3 周期性相對(duì)運(yùn)動(dòng)初始化條件求解

將周期條件截?cái)酁楹蠳項(xiàng)的傅里葉級(jí)數(shù)的形式，即

(50)

(51)

(52)

(53)

(54)

則周期性相對(duì)運(yùn)動(dòng)的初值條件為

(55)

(56)

(57)

(58)

時(shí)域配點(diǎn)法是一種非線性系統(tǒng)的求解方法。它將常微分方程組轉(zhuǎn)化為非線性代數(shù)方程組，再用其他求解非線性代數(shù)方程組的方法對(duì)得到的方程組進(jìn)行求解，避免了直接求解常微分方程組，對(duì)于周期解的處理十分有效。該方法在文獻(xiàn)[10-11]中被提出，本文不對(duì)此作詳細(xì)說明。

(59)

利用時(shí)域配點(diǎn)法進(jìn)行相應(yīng)變換，式(59)左端有

(60)

因此，微分方程組轉(zhuǎn)化后的非線性代數(shù)方程組表示為

(61)

求解該非線性代數(shù)方程組，可得4×(2N+1)個(gè)配點(diǎn)的值，有

(62)

(63)

由此可求得4×(2N+1)個(gè)傅里葉系數(shù)，將其代入式(55)，(56)，(57)，(58)可得周期運(yùn)動(dòng)初值。

采用LM方法對(duì)式(61)的非線性代數(shù)方程組

(64)

(65)

進(jìn)行求解，可避免采用牛頓-拉斐森迭代法(NR)進(jìn)行求解時(shí)由于方程形式過于復(fù)雜，雅可比矩陣出現(xiàn)奇異或病態(tài)嚴(yán)重，無法繼續(xù)迭代的情況。

圖2 參考橢圓軌道偏心率為0.02仿真結(jié)果Fig.2 Simulation results of elliptical orbit with eccentricity of 0.02

采用一種全局收斂的LM算法[12-15]，考慮矢量形式的非線性方程組F(x)=0。利用線性方程組的解dk作為在點(diǎn)xk上的一個(gè)搜索方向：((Jk)T·Jk+μkI)dk=-(Jk)TFk。其中，I為單位矩陣，Jk為雅可比矩陣。通過引進(jìn)非負(fù)參數(shù)μk可保證矩陣((Jk)TJk+μkI)非奇異，同時(shí)也能避免出現(xiàn)過大的‖dk‖。同時(shí)，對(duì)于迭代步((Jk)TJk+μkI)dk=-(Jk)TFk，選取參數(shù)μk=αk(θ‖F(xiàn)k‖+(1-θ)‖(Jk)TFk‖),θ∈[0,1]。其中，對(duì)于參數(shù)αk的選取采用信賴域技巧來修正。定義價(jià)值函數(shù)φ(x)=‖F(xiàn)(x)‖2，第k步迭代的實(shí)際下降量Ak=‖F(xiàn)k‖2-‖F(xiàn)(xk+dk)‖2,預(yù)估下降量Pk=‖F(xiàn)k‖2-‖F(xiàn)k+Jkdk‖2，實(shí)際下降量與預(yù)估下降量的比值rk=Ak/Pk，由此來決定是否接受試探步dk和調(diào)整迭代過程中αk因子的大小。rk越大，表明目標(biāo)函數(shù)φ(x)=‖F(xiàn)(x)‖2下降越多，則接受dk，同時(shí)希望下一試探步dk+1更長(zhǎng)，故減小αk；反之，rk越小，則拒絕dk，同時(shí)增大αk。

4 仿真算例

對(duì)考慮正常引力勢(shì)的非線性3，4階項(xiàng)進(jìn)行仿真，取橢圓參考軌道長(zhǎng)半軸a=8 000 km，過近地點(diǎn)時(shí)刻tp=0，偏心率e=0.02，其余3個(gè)軌道根數(shù)均取為0，將傅氏級(jí)數(shù)截?cái)酁镹=3，在一個(gè)周期上均勻地取7個(gè)配點(diǎn)，仿真6個(gè)軌道周期，結(jié)果如圖2所示。

為定量描述每個(gè)周期的軌道漂移量的大小，以衡量周期運(yùn)動(dòng)的好壞，引入如下指標(biāo)來粗略描述經(jīng)過一個(gè)周期T的漂移量，即

(66)

經(jīng)計(jì)算可得ΔT=0.315 9 m。由此可見：每過一圈的漂移量是很小的。這充分說明周期條件初值在用于軌道計(jì)算時(shí)能保持較好的周期性，在軌道保持和控制等方面具有工程應(yīng)用價(jià)值。仿真結(jié)果驗(yàn)證了模型和求解方法的有效性。

5 結(jié)束語

本文主要基于哈密爾頓力學(xué)建立橢圓參考軌道的相對(duì)運(yùn)動(dòng)模型，利用正則攝動(dòng)理論將哈密爾頓函數(shù)分解為線性項(xiàng)和攝動(dòng)項(xiàng)。線性項(xiàng)的精確解很容易得到，攝動(dòng)項(xiàng)則可通過線性項(xiàng)的修正獲得，且各攝動(dòng)項(xiàng)滿足線性疊加的特點(diǎn)，各種軌道攝動(dòng)和人為控制易于加到運(yùn)動(dòng)模型中。同時(shí)，針對(duì)模型特點(diǎn)推導(dǎo)了周期性相對(duì)運(yùn)動(dòng)條件，利用時(shí)域配點(diǎn)法和LM方法對(duì)周期性相對(duì)運(yùn)動(dòng)的初始值進(jìn)行求解，得到了不消耗任何燃料的周期性繞飛軌道，為空間交會(huì)、編隊(duì)飛行、機(jī)動(dòng)伴飛等航天領(lǐng)域的研究提供了有益參考。