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

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

0 引言

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

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

1 橢圓參考軌道相對運動建模

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

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

(1)

(2)

動能表示為

(3)

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

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

(9)

由式(9)可見，哈密爾頓函數(shù)不顯含時間t，則根據(jù)哈密爾頓動力學(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)為不含時間t的函數(shù)，(α1,α2,…,αf)為積分常數(shù)；t為時間；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方向單獨分離出來，W(x,y,z)=W′(x,y)+W3(z)，則HJ方程變?yōu)?/p>

(15)

(16)

將式(15)積分得

(17)

再令

(18)

同時令

(19)

由式(19)可得

(20)

(21)

積分得

W1(x)=

(22)

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

(23)

(24)

(25)

(26)

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

(27)

(28)

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

(29)

(30)

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

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

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

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

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

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

(31)

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

(32)

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

(33)

(34)

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

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

(35)

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

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

(36)

考慮J2攝動項的攝動哈密爾頓函數(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為參考軌道傾角；θ為真近點角。

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

(38)

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

(39)

(40)

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

(41)

(42)

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

2 周期性相對運動條件

橢圓基準(zhǔn)軌道的動力學(xué)模型為式(38)，(39)，經(jīng)分析不難發(fā)現(xiàn)，獲得周期為T的周期性相對運動的條件為：?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 周期性相對運動初始化條件求解

將周期條件截斷為含有N項的傅里葉級數(shù)的形式，即

(50)

(51)

(52)

(53)

(54)

則周期性相對運動的初值條件為

(55)

(56)

(57)

(58)

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

(59)

利用時域配點法進行相應(yīng)變換，式(59)左端有

(60)

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

(61)

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

(62)

(63)

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

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

(64)

(65)

進行求解，可避免采用牛頓-拉斐森迭代法(NR)進行求解時由于方程形式過于復(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作為在點xk上的一個搜索方向：((Jk)T·Jk+μkI)dk=-(Jk)TFk。其中，I為單位矩陣，Jk為雅可比矩陣。通過引進非負參數(shù)μk可保證矩陣((Jk)TJk+μkI)非奇異，同時也能避免出現(xiàn)過大的‖dk‖。同時，對于迭代步((Jk)TJk+μkI)dk=-(Jk)TFk，選取參數(shù)μk=αk(θ‖F(xiàn)k‖+(1-θ)‖(Jk)TFk‖),θ∈[0,1]。其中，對于參數(shù)αk的選取采用信賴域技巧來修正。定義價值函數(shù)φ(x)=‖F(xiàn)(x)‖2，第k步迭代的實際下降量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，實際下降量與預(yù)估下降量的比值rk=Ak/Pk，由此來決定是否接受試探步dk和調(diào)整迭代過程中αk因子的大小。rk越大，表明目標(biāo)函數(shù)φ(x)=‖F(xiàn)(x)‖2下降越多，則接受dk，同時希望下一試探步dk+1更長，故減小αk；反之，rk越小，則拒絕dk，同時增大αk。

4 仿真算例

對考慮正常引力勢的非線性3，4階項進行仿真，取橢圓參考軌道長半軸a=8 000 km，過近地點時刻tp=0，偏心率e=0.02，其余3個軌道根數(shù)均取為0，將傅氏級數(shù)截斷為N=3，在一個周期上均勻地取7個配點，仿真6個軌道周期，結(jié)果如圖2所示。

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

(66)

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

5 結(jié)束語

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