張海濤,王偉林,張雅聲,王 浩,李 智*

(1. 航天工程大學(xué), 北京 101416; 2. 酒泉衛(wèi)星發(fā)射中心, 甘肅 酒泉 735000)

由于地球靜止軌道(geosynchronous orbit, GEO)的軌道優(yōu)勢(shì),許多高價(jià)值航天器都位于GEO區(qū)域。2021年,歐洲航天局空間碎片辦公室在《歐洲年度空間環(huán)境報(bào)告》中指出,截至2020年底,共有873個(gè)GEO空間目標(biāo)[1]。GEO空間目標(biāo)距離地面近36 000 km,給地面高清觀測(cè)[2-3]帶來(lái)巨大挑戰(zhàn)。

GEO區(qū)域的成像衛(wèi)星能夠在不受大氣干擾的情況下,拍攝更清晰的GEO空間目標(biāo)圖像,并能更好地監(jiān)測(cè)GEO航天器的工作狀態(tài)。因此,部署成像衛(wèi)星接近GEO航天器,實(shí)現(xiàn)對(duì)GEO航天器的持續(xù)觀測(cè)具有重要意義。利用GEO的軌道周期特點(diǎn),成像衛(wèi)星可以在整個(gè)自然繞飛任務(wù)中以良好的觀測(cè)角對(duì)GEO航天器進(jìn)行連續(xù)觀測(cè)。

根據(jù)用于表示相對(duì)運(yùn)動(dòng)的狀態(tài)量的類型,研究相對(duì)運(yùn)動(dòng)的模型可以分為兩類:第一類是用軌道參數(shù)作為狀態(tài)量的相對(duì)運(yùn)動(dòng)模型;第二類是用相對(duì)狀態(tài)矢量作為狀態(tài)量的相對(duì)運(yùn)動(dòng)模型,例如用直角坐標(biāo)系中的相對(duì)位置和速度[4]作為狀態(tài)矢量。第一類模型[5-8]的優(yōu)點(diǎn)是易于分析攝動(dòng)對(duì)相對(duì)運(yùn)動(dòng)的影響。因此,這種模型具有較高的精度,并能有效避免偏差的長(zhǎng)期累積。這種模型的缺點(diǎn)是不便于建立狀態(tài)微分方程,不利于研究連續(xù)小推力作用下的相對(duì)運(yùn)動(dòng)。在第二類模型[9-11]中,Clohessy-Wiltshire (CW)方程應(yīng)用最為廣泛,CW方程的優(yōu)點(diǎn)在于:建立狀態(tài)方程后,便于結(jié)合最優(yōu)控制理論開(kāi)展成像衛(wèi)星接近GEO航天器的研究,也便于對(duì)飛行編隊(duì)進(jìn)行研究。缺點(diǎn)是CW方程的偏差較大,且偏差隨時(shí)間不斷累積[12]。

在繞飛軌跡規(guī)劃問(wèn)題上,劉猛等[13]給出了航天器基于最小沖量進(jìn)入對(duì)空間目標(biāo)繞飛軌道的策略,在軌道機(jī)動(dòng)燃料消耗方面有借鑒價(jià)值。黃藝等[14]研究了對(duì)非合作目標(biāo)繞飛任務(wù)的姿軌耦合控制策略,給出能夠?qū)ν饨绺蓴_和不確定參數(shù)具有很好抑制效果的控制方案。譚天樂(lè)等[15-17]給出近圓軌道和大橢圓軌道的相對(duì)運(yùn)動(dòng)建模和解析方法,提供了對(duì)空間非合作目標(biāo)高精度的相對(duì)懸停和循跡繞飛控制方法。劉濤等[18]提出一種用于非合作目標(biāo)慣性指向軸位置捕獲的繞飛方法。常燕等[19]研究了對(duì)橢圓軌道上非合作目標(biāo)進(jìn)行長(zhǎng)期繞飛檢測(cè)的相對(duì)運(yùn)動(dòng)軌道構(gòu)型設(shè)計(jì)與構(gòu)型保持問(wèn)題。梁靜靜等[20]基于粒子群算法研究了雙脈沖繞飛的優(yōu)化問(wèn)題和安全繞飛軌跡設(shè)計(jì)問(wèn)題。王功波等[21]在連續(xù)小推力條件下,針對(duì)圓軌道推導(dǎo)了滿足快速繞飛條件的空間圓編隊(duì)動(dòng)力學(xué)模型。Dang[22-23]推導(dǎo)了在任意開(kāi)普勒軌道上相對(duì)運(yùn)動(dòng)的新?tīng)顟B(tài)轉(zhuǎn)換矩陣,給出了Tschauner-Hempel(TH)方程的解。

在CW方程的偏差問(wèn)題上,通過(guò)修正CW方程的非球形攝動(dòng)偏差和非線性二次項(xiàng)偏差改進(jìn)CW方程;在軌跡規(guī)劃問(wèn)題上,計(jì)算對(duì)GEO航天器繞飛的初始相位區(qū)間,成像衛(wèi)星在連續(xù)小推力作用下接近GEO航天器,以合適的初始相位對(duì)其開(kāi)始繞飛,實(shí)現(xiàn)對(duì)GEO航天器的自然繞飛并持續(xù)以有利的光照條件對(duì)其觀測(cè),如圖1所示。

圖1 成像衛(wèi)星P在整個(gè)繞飛任務(wù)中以良好的觀測(cè)角對(duì)GEO航天器E進(jìn)行拍照Fig.1 The optical satellite P can take pictures of the GEO spacecraft E with favorable observation angles throughout the entire fly-around mission

1 CW方程

CW方程采用當(dāng)?shù)卮怪碑?dāng)?shù)厮?local vertical local horizontal, LVLH)坐標(biāo)系,該坐標(biāo)系以航天器質(zhì)心為原點(diǎn),x軸沿地心到航天器質(zhì)心的方向;z軸垂直于軌道平面,與瞬時(shí)角動(dòng)量方向相同,y軸方向根據(jù)右手定則確定。

兩個(gè)非常接近的航天器分別記為M和S,M為主航天器,S為輔航天器,存在:

(1)

(2)

式中:

(3)

(4)

其中,ω為航天器M的軌道角速度。

記t0時(shí)刻,S的相對(duì)狀態(tài)矢量為X(t0),則t時(shí)刻,航天器S的相對(duì)狀態(tài)矢量為:

(5)

式中:

(6)

2 改進(jìn)的CW方程

非線性、參考軌道的偏心率和攝動(dòng)是CW方程偏差的三個(gè)來(lái)源。在成像衛(wèi)星對(duì)GEO航天器的接近和繞飛的任務(wù)中,以GEO航天器為原點(diǎn)建立LVLH坐標(biāo)系,則CW方程偏差來(lái)源主要為非線性項(xiàng)和攝動(dòng)項(xiàng)。記P為成像衛(wèi)星,E為被成像的GEO航天器。已知P和E的初始軌道參數(shù),改進(jìn)的CW(improved CW, ICW)方程研究相對(duì)運(yùn)動(dòng)的步驟如下:

步驟1：以航天器E的質(zhì)心為原點(diǎn)建立LVLH坐標(biāo)系,記為R;

步驟2：計(jì)算成像衛(wèi)星P的初始狀態(tài)矢量,即在LVLH坐標(biāo)系中成像衛(wèi)星P相對(duì)GEO航天器E的位置和速度矢量:

(7)

據(jù)文獻(xiàn)[22]的方法可求得:

(8)

步驟3：GEO航天器受到的主要攝動(dòng)為非球形、三體和太陽(yáng)光壓攝動(dòng)。然而,由于GEO軌道的軌道特性,非球形攝動(dòng)是導(dǎo)致兩航天器存在相對(duì)運(yùn)動(dòng)偏差的最主要因素?？紤]J2和J22項(xiàng)攝動(dòng),將CW方程中的萬(wàn)有引力常數(shù)μ修正為ICW方程中的萬(wàn)有引力常數(shù)μ′:

(9)

式中,μ=3.986 005×1014m3/s2為萬(wàn)有引力常數(shù),rE為地球的赤道半徑,ac為GEO軌道的半長(zhǎng)軸,λ為GEO航天器E的星下點(diǎn)地理經(jīng)度,J2=-1.082 627×10-3,J22=1.815 528×10-6。

步驟4：通過(guò)對(duì)CW方程初始狀態(tài)矢量的修正來(lái)補(bǔ)償非線性長(zhǎng)期項(xiàng)偏差。將非線性二次項(xiàng)假設(shè)為虛加的攝動(dòng),通過(guò)虛攝動(dòng)解可求解對(duì)CW方程初始狀態(tài)矢量的修正項(xiàng)。對(duì)CW方程初始狀態(tài)矢量的修正項(xiàng)為:

(10)

式中,nL為L(zhǎng)VLH坐標(biāo)系相對(duì)ECI坐標(biāo)系的旋轉(zhuǎn)角速度,ρ2為成像衛(wèi)星P與航天器E的相對(duì)距離的平方,β0為初始時(shí)刻y軸沿z軸反方向到成像衛(wèi)星P位置矢量的角度。

步驟5：將兩個(gè)獨(dú)立的修正結(jié)合起來(lái),以補(bǔ)償非線性偏差和非球形攝動(dòng)偏差。ICW方程的初始狀態(tài)矢量為:

(11)

在ICW方程中,通過(guò)步驟3補(bǔ)償非球形攝動(dòng)偏差,通過(guò)步驟4補(bǔ)償非線性偏差,2.1節(jié)和2.2節(jié)分別研究步驟3和步驟4補(bǔ)償攝動(dòng)偏差和非線性偏差的原理。

2.1 補(bǔ)償非球形攝動(dòng)偏差

定義航天器在離地球無(wú)限遠(yuǎn)處的重力勢(shì)能為0。r為ECI坐標(biāo)系中的位置矢量,r為位置矢量的大小。設(shè)r處單位質(zhì)量的重力勢(shì)能為:

(12)

(13)

化簡(jiǎn)后,引力勢(shì)為:

各項(xiàng)的權(quán)重為:

(15)

對(duì)于GEO航天器,權(quán)量最大的五項(xiàng)分別是J2,J22,J3,J31和J33,權(quán)重比為:

A2∶A22∶A3∶A31∶A33=11 053 ∶64 ∶3 ∶7 ∶6

(16)

因此,J2和J22是非球形攝動(dòng)中權(quán)重最大的項(xiàng)。將這兩項(xiàng)引入ICW方程中,包含J2和J22項(xiàng)的引力勢(shì)為:

(17)

式中,λ22=-14.929°。對(duì)于兩個(gè)相對(duì)接近的GEO航天器,其星下點(diǎn)的地理經(jīng)度近似相等。由于sin2φ≈0,cos2φ≈1,式(17)可簡(jiǎn)化為:

(18)

J2和J22對(duì)相對(duì)運(yùn)動(dòng)的影響體現(xiàn)在與CW方程相比,ICW方程中的萬(wàn)有引力常數(shù)修正為:

(19)

2.2 補(bǔ)償非線性偏差

在CW方程中,為了得到線性的狀態(tài)方程,重力加速度的泰勒級(jí)數(shù)展開(kāi)式中的非線性項(xiàng)在式(1)的條件下被舍棄,二次項(xiàng)是偏差的主要來(lái)源。被舍棄的二次項(xiàng)在位置分量上產(chǎn)生常數(shù)項(xiàng)、長(zhǎng)期項(xiàng)和短周期項(xiàng)三種類型的偏差。長(zhǎng)期項(xiàng)只出現(xiàn)在y軸方向上,是對(duì)精度影響最大的偏差。偏差的長(zhǎng)期項(xiàng)可通過(guò)對(duì)CW方程的初始狀態(tài)矢量的修正來(lái)補(bǔ)償,利用攝動(dòng)法可計(jì)算用于補(bǔ)償CW方程中非線性偏差的修正項(xiàng)。攝動(dòng)法中,將非線性二次項(xiàng)作為CW方程中的虛加攝動(dòng)力。通過(guò)虛攝動(dòng)解求解用于補(bǔ)償CW方程中非線性長(zhǎng)期項(xiàng)偏差的初始狀態(tài)矢量的修正項(xiàng)。

(20)

其中,R為在r處單位質(zhì)量的勢(shì)能。成像衛(wèi)星P相對(duì)GEO航天器E的位置矢量為ρ=rP-rc。根據(jù)式(12),對(duì)于航天器E和成像衛(wèi)星P,存在:

(21)

(22)

其中,fc和fP分別為航天器E和成像衛(wèi)星P受到的攝動(dòng)力,如太陽(yáng)光壓攝動(dòng)等。式(21)與式(22)相減得:

(23)

式中,Δf為成像衛(wèi)星P和航天器E受到的攝動(dòng)偏差。ICW方程在步驟3中對(duì)主要攝動(dòng)偏差進(jìn)行修正,在步驟4中對(duì)非線性偏差進(jìn)行修正。因此,在修正非線性偏差時(shí),令Δf=0。將式(23)從ECI坐標(biāo)系轉(zhuǎn)換到LVLH坐標(biāo)系:

(24)

(25)

在CW方程中,式(25)可化簡(jiǎn)為:

(26)

記CW方程中的狀態(tài)矢量為:

(27)

解為:

(28)

保留式(25)等號(hào)右邊的二次項(xiàng),舍棄高階項(xiàng)。將式(25)化簡(jiǎn)為:

(29)

根據(jù)式(29),將重力加速度泰勒級(jí)數(shù)展開(kāi)式中的非線性二次項(xiàng)作為CW方程中虛加的攝動(dòng)力,虛加攝動(dòng)力為:

(30)

將式(29)簡(jiǎn)記為:

(31)

式(31)為CW方程的非線性偏差傳播方程。記CW方程的非線性偏差為:

(32)

根據(jù)式(5)可得式(31)的解為:

(33)

CW方程中被舍棄的二次項(xiàng)導(dǎo)致y軸方向偏差的長(zhǎng)期項(xiàng)為:

(34)

為了修正y軸方向偏差的長(zhǎng)期項(xiàng),令ynl-long=0,可得非線性偏差的零累積準(zhǔn)則:

(35)

因此,根據(jù)虛攝動(dòng)解可得到對(duì)CW方程初始狀態(tài)矢量的修正項(xiàng)如式(10)所示。

2.3 偏差修正驗(yàn)證

為了研究ICW方程相對(duì)CW方程的改進(jìn)效果,選取成像衛(wèi)星P和GEO航天器E進(jìn)行仿真驗(yàn)證,其軌道半長(zhǎng)軸a、偏心率e、軌道傾角i、升交點(diǎn)赤經(jīng)Ω、近地點(diǎn)幅角w、平近地點(diǎn)角m如表1所示。

表1 成像衛(wèi)星P與GEO航天器E的軌道參數(shù)

LVLH坐標(biāo)系中成像衛(wèi)星P的初始狀態(tài)矢量如表2所示。初始時(shí)刻從y軸沿z軸反方向到成像衛(wèi)星P的位置矢量的角度為0°,X為CW方程的初始狀態(tài)矢量,Y為ICW方程的初始狀態(tài)矢量。

表2 成像衛(wèi)星P的初始狀態(tài)矢量

2.3.1 非球形攝動(dòng)偏差修正驗(yàn)證

為了驗(yàn)證ICW方程中步驟3的非球形攝動(dòng)解,將非球形攝動(dòng)解與高精度軌道傳播(high-precision orbit propagator, HPOP)模型作對(duì)比。HPOP模型中考慮前21階非球形攝動(dòng)、太陽(yáng)光壓攝動(dòng)和三體引力攝動(dòng)。利用二體模型、HPOP模型和步驟3中的非球形攝動(dòng)解,分別在ECI坐標(biāo)系中計(jì)算成像衛(wèi)星P和GEO航天器E的位置。將P在ECI坐標(biāo)系中的位置分別轉(zhuǎn)換到LVLH坐標(biāo)系R中,P在R中的位置即為P相對(duì)于E的位置,將HPOP模型得到的相對(duì)位置視為真實(shí)相對(duì)位置。在仿真場(chǎng)景中,航天器星下點(diǎn)的地理經(jīng)度為165″E。

相比二體模型,非球形攝動(dòng)解在相對(duì)位置的x、y和z方向?qū)鹊奶岣叻謩e如圖2、圖3和圖4所示。圖中紅色實(shí)線表示非球形攝動(dòng)解相對(duì)HPOP模型的偏差,藍(lán)色虛線表示二體模型相對(duì)HPOP模型的偏差。通過(guò)比較可得,非球形攝動(dòng)是攝動(dòng)偏差的最主要因素,非球形攝動(dòng)解提高了CW方程在LVLH坐標(biāo)系x、y軸方向的精度。

圖2 非球形攝動(dòng)解對(duì)x方向精度的提高Fig.2 Improvement of the non-spherical perturbation solution in x-coordinate

圖3 非球形攝動(dòng)解對(duì)y方向精度的提高Fig.3 Improvement of the non-spherical perturbation solution in y-coordinate

圖4 非球形攝動(dòng)解對(duì)z方向精度的提高Fig.4 Improvement of the non-spherical perturbation solution in z-coordinate

利用ICW方程和CW方程仿真成像衛(wèi)星P相對(duì)于GEO航天器E的運(yùn)動(dòng),比較三個(gè)軌道周期內(nèi)相對(duì)位置的x、y、z坐標(biāo)差異。由于CW方程存在非線性偏差,為了展現(xiàn)ICW方程中步驟3對(duì)CW方程的改進(jìn)效果,本節(jié)中不對(duì)ICW方程進(jìn)行非線性偏差修正。ICW方程與CW方程中相對(duì)位置的x、y、z坐標(biāo)差值分別如圖5、圖6和圖7所示。

圖5 ICW方程和CW方程在x軸方向的差值Fig.5 Deviation in x-coordinate between ICW equations and CW equations

圖6 ICW方程和CW方程在y軸方向的差值Fig.6 Deviation in y-coordinate between ICW equations and CW equations

圖7 ICW方程和CW方程在z軸方向的差值Fig.7 Deviation in z-coordinate between ICW equations and CW equations

從圖5～7中可以看出,ICW方程與CW方程中相對(duì)位置的x、y、z坐標(biāo)差值呈周期性變化,且差值幅值逐漸增大。x、y坐標(biāo)差值幅值的增長(zhǎng)率在同一個(gè)數(shù)量級(jí),z坐標(biāo)幅值的增長(zhǎng)率比x、y坐標(biāo)幅值的增長(zhǎng)率小一個(gè)數(shù)量級(jí)。

2.3.2 非線性偏差修正驗(yàn)證

CW方程和ICW方程計(jì)算出的成像衛(wèi)星P相對(duì)GEO航天器E的位置的y坐標(biāo)如圖8所示。圖中紅色實(shí)線代表ICW方程的y坐標(biāo),藍(lán)色虛線代表CW方程的y坐標(biāo)。ICW方程與CW方程的y坐標(biāo)差值如圖9所示。

圖8 ICW方程和CW方程中的y坐標(biāo)值Fig.8 y-coordinate in ICW equations and CW equations

圖9 ICW方程和CW方程在y軸方向的 非線性偏差值Fig.9 y-coordinate nonlinear deviation between ICW equations and CW equations

為了驗(yàn)證ICW方程步驟4的改進(jìn)效果,將CW方程和ICW方程計(jì)算的相對(duì)位置的y坐標(biāo)與真值進(jìn)行比較。由于CW方程存在攝動(dòng)偏差,為了展現(xiàn)ICW方程中步驟4對(duì)CW方程的改進(jìn)效果,本節(jié)不對(duì)ICW方程進(jìn)行攝動(dòng)偏差修正。如圖10所示,紅色實(shí)線表示ICW方程中相對(duì)位置的y坐標(biāo)與真值的差值,藍(lán)色虛線表示CW方程中相對(duì)位置的y坐標(biāo)與真值的差值。

圖10 ICW方程和CW方程在y軸方向與真值的差值Fig.10 y-coordinate deviation of ICW equations and CW equations from the truth value

成像衛(wèi)星P與GEO航天器E的初始距離為200 km,3個(gè)軌道周期后,ICW方程中相對(duì)位置的y坐標(biāo)相對(duì)真值出現(xiàn)36 m的長(zhǎng)期偏差,而CW方程的長(zhǎng)期累積偏差為33 563 m。ICW方程的長(zhǎng)期偏差僅約為CW方程的0.1%。通過(guò)比較,驗(yàn)證了ICW方程對(duì)非線性偏差長(zhǎng)期項(xiàng)補(bǔ)償?shù)挠行浴?/p>

ICW方程補(bǔ)償了相對(duì)位置y軸方向二次非線性偏差的長(zhǎng)期項(xiàng)和x、y軸方向上主要的攝動(dòng)偏差。ICW方程中,在建立狀態(tài)微分方程時(shí),不對(duì)非線性偏差的短期項(xiàng)進(jìn)行修正。在有更高的精度需求時(shí),可用短期項(xiàng)偏差修正結(jié)果狀態(tài)矢量。這樣處理的優(yōu)點(diǎn)在于:得到的控制系統(tǒng)是線性時(shí)不變系統(tǒng),同時(shí)在最終的狀態(tài)矢量中加入短周期項(xiàng)可保證結(jié)果的準(zhǔn)確性。

3 軌跡規(guī)劃

對(duì)GEO航天器繞飛持續(xù)觀測(cè)任務(wù)的軌跡規(guī)劃分為兩段:成像衛(wèi)星抵近GEO航天器和成像衛(wèi)星對(duì)GEO航天器繞飛。為了實(shí)現(xiàn)燃料消耗最少,成像衛(wèi)星最理想的軌跡規(guī)劃方案是:在連續(xù)小推力控制下接近GEO航天器,然后對(duì)GEO航天器自然繞飛完成持續(xù)觀測(cè)任務(wù)。

3.1 成像衛(wèi)星對(duì)GEO航天器繞飛

3.1.1 觀測(cè)角和太陽(yáng)的運(yùn)動(dòng)

觀測(cè)角定義如圖11所示。定義α為成像衛(wèi)星對(duì)GEO航天器的觀測(cè)角。0°觀測(cè)角表示:成像衛(wèi)星、GEO航天器和太陽(yáng)共線,成像衛(wèi)星處于中間位置。顯然,越接近0°,成像衛(wèi)星對(duì)GEO航天器的觀測(cè)越有利;越接近180°,成像衛(wèi)星對(duì)GEO航天器的觀測(cè)越不利。成像衛(wèi)星為了在整個(gè)繞飛階段都能對(duì)GEO航天器成像,需觀測(cè)角始終小于60°。成像衛(wèi)星能較好地對(duì)GEO航天器成像的觀測(cè)角的最大值取決于成像衛(wèi)星相機(jī)的成像能力和GEO航天器的材質(zhì),通常情況下,成像衛(wèi)星在觀測(cè)角小于60°的情況下,可實(shí)現(xiàn)對(duì)GEO航天器成像。

圖11 觀測(cè)角的定義Fig.11 Observation angle

在ECI坐標(biāo)系中,太陽(yáng)的軌道參數(shù)為:

(36)

其中,T為儒略世紀(jì),d為儒略日。據(jù)太陽(yáng)的軌道參數(shù),可得太陽(yáng)在ECI坐標(biāo)系中的位置矢量。由于太陽(yáng)到地心的距離遠(yuǎn)大于GEO航天器到地心的距離,假設(shè)地心到太陽(yáng)的方向與GEO航天器E到太陽(yáng)的方向重合。太陽(yáng)在以航天器E為原點(diǎn)的LVLH坐標(biāo)系R中的方向矢量為:

dsun=AECI→LVLH·kECI

(37)

式中,kECI為太陽(yáng)在ECI中的方向矢量。

3.1.2 繞飛階段的觀測(cè)角

設(shè)計(jì)成像衛(wèi)星P在GEO航天器E軌道平面內(nèi)的繞飛軌跡。設(shè)tp時(shí)刻太陽(yáng)方向矢量在LVLH坐標(biāo)系R的xoy平面上的投影與x軸重合,成像衛(wèi)星P在t0到tp時(shí)刻完成對(duì)航天器E的抵近,從tp時(shí)刻開(kāi)始繞飛。dP為航天器E到成像衛(wèi)星P的位置矢量,記θxy為dP在tp時(shí)刻的初始相位角,即θxy為tp時(shí)刻x軸沿z軸反方向到dP矢量的角度,如圖12所示。δ為太陽(yáng)的赤緯,即太陽(yáng)的星下點(diǎn)緯度,顯然δ∈(-23°26′,+23°26′)。t時(shí)刻,dP在LVLH坐標(biāo)系R中的x、y、z坐標(biāo)為:

圖12 太陽(yáng)和成像衛(wèi)星的相對(duì)位置Fig.12 Position of the sun and the optical satellite

dP=[x(t),y(t),z(t)]

(38)

(39)

(40)

當(dāng)θxy=0°時(shí),成像衛(wèi)星P在GEO航天器E的一個(gè)軌道周期內(nèi)對(duì)E的觀測(cè)角如圖13所示。

圖13 在θxy=0°時(shí)成像衛(wèi)星P對(duì)GEO航天器E的觀測(cè)角Fig.13 Observation angles of the optical satellite P to the GEO spacecraft E when θxy=0°

太陽(yáng)的赤緯在-23°26′到+23°26′間變化。記αm為一個(gè)軌道周期內(nèi)P對(duì)E的最大觀測(cè)角,αm隨太陽(yáng)赤緯而變化,如圖14所示。當(dāng)θxy=0°時(shí),αm在δ=0°時(shí)取得最小值19.47°,在δ=±23°26′時(shí)取得最大值26.35°。因此,成像衛(wèi)星P在tp時(shí)刻以θxy=0開(kāi)始對(duì)GEO航天器E繞飛,成像衛(wèi)星P對(duì)GEO航天器E的觀測(cè)角始終小于26.35°。

圖14 成像衛(wèi)星P對(duì)GEO航天器E的最大 觀測(cè)角隨太陽(yáng)赤緯的變化Fig.14 Maximum observation angle of the optical satellite P to the GEO spacecraft E varies with the declination of the sun δ

tp時(shí)刻,成像衛(wèi)星P以初始相位角θxy≠0°開(kāi)始對(duì)GEO航天器E繞飛。θxy取不同值時(shí),一個(gè)軌道周期內(nèi)的最大觀測(cè)角αm隨太陽(yáng)赤緯δ變化,如圖15所示。因此,無(wú)論太陽(yáng)的赤緯δ取何值,如果繞飛的初始相位θxy∈(-37.2°,37.2°),在整個(gè)繞飛任務(wù)中,P對(duì)E的觀測(cè)角始終小于60°,P可以以有利的觀測(cè)角對(duì)E拍照。因此,可對(duì)GEO航天器持續(xù)觀測(cè)的成像衛(wèi)星的繞飛初始相位角區(qū)間為θxy∈(-37.2°,37.2°)。

圖15 不同θxy時(shí)最大觀測(cè)角αm隨太陽(yáng)赤緯的變化Fig.15 The maximum observation angle in an orbital period αm varies with the declination of the sun δ at different θxy

3.1.3 繞飛軌道參數(shù)確定

設(shè)計(jì)成像衛(wèi)星P在GEO航天器E軌道平面內(nèi)的繞飛。在LVLH坐標(biāo)系R中,繞飛相對(duì)軌跡為橢圓,記rxy為相對(duì)軌跡的半長(zhǎng)軸,記yoff為橢圓相對(duì)軌道的中心與航天器E的距離。記Δe、Δi、ΔΩ、Δω和Δm為P與E軌道的偏心率、軌道傾角、升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角和平近地點(diǎn)角差值:

(41)

已知GEO航天器E的軌道參數(shù),根據(jù)任務(wù)要求的rxy、θxy和yoff,設(shè)計(jì)成像衛(wèi)星P的相對(duì)軌跡。根據(jù)式(41)可得Δa、Δe、Δi、ΔΩ、Δω、Δm,進(jìn)而可得成像衛(wèi)星P在tp時(shí)刻的軌道參數(shù)。

3.2 成像衛(wèi)星抵近GEO航天器

據(jù)3.1節(jié)的方法,可得成像衛(wèi)星P在tp時(shí)刻的軌道參數(shù),已知成像衛(wèi)星在t0時(shí)刻的軌道參數(shù)。以成像衛(wèi)星P的燃料消耗為優(yōu)化指標(biāo),計(jì)算成像衛(wèi)星接近GEO航天器的控制策略,從t0到tp的軌跡優(yōu)化是一個(gè)兩點(diǎn)邊界問(wèn)題。狀態(tài)方程為:

(42)

支付函數(shù)為:

(43)

根據(jù)式(43),哈密頓函數(shù)可寫(xiě)為:

(44)

式中,λ是協(xié)態(tài)變量。

根據(jù)最優(yōu)化的必要條件,協(xié)態(tài)方程為:

(45)

(46)

把式(46)代入式(42),得:

(47)

據(jù)式(47)和式(45),可得:

(48)

式(48)的解析解為:

(49)

式(49)可寫(xiě)為:

(50)

已知ICW方程中的初始狀態(tài)矢量Y(t0)和終端狀態(tài)矢量Y(tp),可得:

λ*(t0)=Π12(tp)-1[Y*(tp)-Π11(tp)Y(t0)]

(51)

根據(jù)式(46)、式(50)和式(51),可得最優(yōu)控制策略:

U*(t)=-BT[Π21(t)Y(t0)+Π22(t)λ*(t0)]

(52)

4 仿真

從t0到tp,成像衛(wèi)星P在連續(xù)小推力的作用下接近GEO航天器E;自tp時(shí)刻開(kāi)始,成像衛(wèi)星P對(duì)GEO航天器自然繞飛。成像衛(wèi)星P的相對(duì)繞飛軌跡在GEO航天器E的軌道平面內(nèi),橢圓相對(duì)軌跡的半長(zhǎng)軸為20 km,初始相位角θxy=0°。在該場(chǎng)景中,t0為GEO航天器星下點(diǎn)當(dāng)?shù)貢r(shí)間2021年6月10日00:00:00,tp為2021年6月11日12:00:00。成像衛(wèi)星P和GEO航天器E在t0時(shí)刻的軌道參數(shù)如表1所示。

根據(jù)CW方程和ICW方程,LVLH坐標(biāo)系R中成像衛(wèi)星P在t0和tp時(shí)的狀態(tài)矢量如表3所示,X為CW方程的狀態(tài)矢量,Y為ICW方程的狀態(tài)矢量。

表3 成像衛(wèi)星P在LVLH坐標(biāo)系中的狀態(tài)矢量

分別基于CW方程和ICW方程進(jìn)行仿真。計(jì)算成像衛(wèi)星P的最優(yōu)控制策略,將成像衛(wèi)星P在ECI坐標(biāo)系中的運(yùn)動(dòng)方程積分得到其位置和速度,再將ECI中得到的位置轉(zhuǎn)換到LVLH坐標(biāo)系中。積分采用龍格-庫(kù)塔法,積分得到的結(jié)果作為成像衛(wèi)星P的真實(shí)位置,進(jìn)而可得成像衛(wèi)星P相對(duì)GEO航天器E的位置以及P對(duì)E繞飛過(guò)程中的觀測(cè)角。將基于CW方程和ICW方程得到的相對(duì)位置和觀測(cè)角進(jìn)行對(duì)比。

4.1 基于CW方程的仿真

基于CW方程進(jìn)行仿真。t0到tp時(shí)間內(nèi),成像衛(wèi)星P在連續(xù)小推力作用下接近GEO航天器E,施加的總速度增量為4.70 m/s。圖16為成像衛(wèi)星P相對(duì)GEO航天器E的位置在LVLH坐標(biāo)系xoy平面中的投影。由于CW方程存在偏差,相對(duì)繞飛軌跡的中心在y軸方向上偏離GEO航天器E。圖17為P到E的距離,tp時(shí)刻,P到E的距離不是預(yù)先規(guī)劃的10 km,而是14 km。圖18為成像衛(wèi)星P在tp后的一個(gè)軌道周期內(nèi)對(duì)GEO航天器E的觀測(cè)角,最大觀測(cè)角為75.57°,成像觀測(cè)仿真任務(wù)失敗。

圖16 基于CW方程仿真的成像衛(wèi)星P相對(duì)GEO 航天器E的位置在LVLH坐標(biāo)系xoy平面中的投影Fig.16 Position of the optical satellite P relative to the GEO spacecraft E on the xoy plane based on CW equations simulation

圖17 基于CW方程仿真的成像衛(wèi)星P與 GEO航天器E的距離Fig.17 Distance between the optical satellite P and the GEO spacecraft E, simulation based on CW equations

圖18 基于CW方程仿真的成像衛(wèi)星P對(duì) GEO航天器E的觀測(cè)角Fig.18 Observation angles of the optical satellite P to the GEO spacecraft E based on CW equations simulation

4.2 基于ICW方程的仿真

基于ICW方程進(jìn)行仿真,圖19為成像衛(wèi)星P從t0到tp的控制推力加速度(工程上一般表述為控制推力),總速度增量為4.67 m/s。圖20～22分別表示成像衛(wèi)星P的推力控制的x、y、z坐標(biāo)分量。圖23為P相對(duì)于E的位置,圖24為P與E間的距離,tp時(shí)刻P與E的距離為10 km,P開(kāi)始對(duì)E自然繞飛。圖25為P對(duì)E在tp后一個(gè)軌道周期內(nèi)的觀測(cè)角,最大觀測(cè)角為22.6°。成像衛(wèi)星P能以良好的觀測(cè)角在整個(gè)繞飛任務(wù)中對(duì)GEO航天器E成像。

圖19 基于ICW方程仿真的成像衛(wèi)星P的控制推力Fig.19 Thrust control of the optical satellite P based on ICW equations simulation

圖20 基于ICW方程仿真的成像衛(wèi)星P在x方向的推力Fig.20 Thrust control of the optical satellite P in x-coordinate simulation based on ICW equations

圖21 基于ICW方程仿真的成像衛(wèi)星P在y方向的推力Fig.21 Thrust control of the optical satellite P in y-coordinate simulation based on ICW equations

圖22 基于ICW方程仿真的成像衛(wèi)星P在z方向的推力Fig.22 Thrust control of the optical satellite P in z-coordinate simulation based on ICW equations

圖23 基于ICW方程仿真的成像衛(wèi)星P相對(duì) GEO航天器E的位置Fig.23 Position of the optical satellite P relative to the GEO spacecraft E based on ICW equations simulation

圖24 基于ICW方程仿真的成像衛(wèi)星P與 GEO航天器E間的距離Fig.24 Distance between the optical satellite P and the GEO spacecraft E based on ICW equations simulation

圖25 基于ICW方程仿真的成像衛(wèi)星P對(duì) GEO航天器E的觀測(cè)角Fig.25 Observation angles of the optical satellite P to the GEO spacecraft E based on ICW equations simulation

5 結(jié)論

通過(guò)修正非球形攝動(dòng)偏差和重力加速度二次長(zhǎng)期項(xiàng)偏差改進(jìn)CW方程。將CW方程、ICW方程與真值作仿真對(duì)比,ICW方程補(bǔ)償了非線性偏差的長(zhǎng)期項(xiàng)和主要的攝動(dòng)偏差。在軌跡規(guī)劃問(wèn)題上,計(jì)算對(duì)GEO航天器繞飛的初始相位角區(qū)間,成像衛(wèi)星以該區(qū)間內(nèi)的初始相位角對(duì)GEO航天器開(kāi)始繞飛,可實(shí)現(xiàn)在整個(gè)繞飛任務(wù)中都能以良好的觀測(cè)角對(duì)GEO航天器進(jìn)行觀測(cè)。選擇合適的初始相位角,分別基于CW方程和ICW方程仿真接近和繞飛的全過(guò)程。在仿真中,成像衛(wèi)星在連續(xù)的小推力作用下接近GEO航天器,然后對(duì)GEO航天器自然繞飛。基于CW方程的仿真任務(wù)失敗,而基于ICW方程達(dá)到預(yù)期目標(biāo)。在基于ICW方程的仿真中,所需要的總速度增量?jī)H為4.67 m/s,工程上具有很強(qiáng)的可行性。ICW方程可滿足精度要求,可用于GEO空間態(tài)勢(shì)感知和在軌服務(wù)任務(wù)規(guī)劃。