地球静止轨道航天器绕飞持续观测任务轨迹规划与控制

张海涛; 王伟林; 张雅声; 王浩; 李智; ZHANG Haitao; WANG Weilin; ZHANG Yasheng; WANG Hao; LI Zhi

doi: 10.11887/j.cn.202401008

张海涛¹ ，王伟林¹ ，张雅声¹ ，王浩² ，李智¹

1. 航天工程大学，北京 101416

2. 酒泉卫星发射中心，甘肃酒泉 735000

基金项目: 国家卓越青年基金资助项目(2017-JCJQ-ZQ-005) ；国家自然科学基金资助项目（61304228）

详细信息

作者简介

张海涛（1989—），男，山东淄博人，工程师，博士，E-mail: zhanghaitaoat@163.com

通讯作者

李智（1973—），男，北京怀柔人，教授，博士，博士生导师，E-mail: zhli_seu@163.com

中图分类号: V529

文献标识码: A

文章编号: 1001-2486(2024)01-074-13

Trajectory planning and control of continuous observation missions for geosynchronous orbit spacecraft fly-around

ZHANG Haitao¹ ， WANG Weilin¹ ， ZHANG Yasheng¹ ， WANG Hao² ， LI Zhi¹

1. Space Engineering University, Beijing 101416 , China

2. Jiuquan Satellite Launch Center, Jiuquan 735000 , China

摘要

针对地球静止轨道（geosynchronous orbit, GEO）航天器的高清观测任务，成像卫星在连续小推力作用下接近GEO航天器，对GEO航天器自然绕飞并以有利的光照条件对其持续观测。针对Clohessy-Wiltshire (CW)方程的偏差问题，通过修正非球形摄动和重力加速度二次长期项偏差对CW方程进行改进，补偿非线性偏差的长期项和主要的摄动项。在轨迹规划问题上，计算绕飞轨迹的初始相位角区间，以保证成像卫星在整个绕飞任务中都能够以良好的观测角观测GEO航天器。基于CW方程和改进的CW方程对成像卫星接近和绕飞GEO航天器全过程进行仿真，基于CW方程的仿真没有达到预期目标；基于改进的CW方程的仿真达到预期目标，全过程所需施加的总速度增量仅为4.67 m/s，工程上具有很强的可行性。

关键词

成像卫星 / 轨迹规划 / 绕飞任务 / 非线性偏差 / 摄动偏差 / 持续观测

Abstract

For high-definition observation of GEO(geosynchronous orbit) spacecraft, the optical satellite approaches the GEO spacecraft with continuous low-thrust and flies around the GEO spacecraft without thrust. On the deviation of Clohessy-Wiltshire (CW) equations, CW equations are improved by accommodating the non-spherical perturbation, and the quadratic terms of the nonlinearity in the differential gravitational acceleration. The secular growth of the nonlinear deviation and the most perturbation deviation have been accommodated in the ICW(improved CW） equations. On trajectory planning, we calculate the initial phase angle interval of the fly-around trajectory to ensure that the optical satellite can take pictures of the GEO spacecraft with favorable observation angles throughout the entire fly-around mission. Simulations are conducted based on CW equations and ICW equations, respectively. The simulation based on the CW equations fails, but the simulation based on the ICW equations succeeds. The total thrust required is only 4.67 m/s, which is highly feasible in engineering.

Keywords

optical satellite / trajectory planning / fly-around mission / nonlinear deviation / perturbation deviation / continuous observation

1 CW方程 2 改进的CW方程 2.1 补偿非球形摄动偏差 2.2 补偿非线性偏差 2.3 偏差修正验证 2.3.1 非球形摄动偏差修正验证 2.3.2 非线性偏差修正验证 3 轨迹规划 3.1 成像卫星对GEO航天器绕飞 3.1.1 观测角和太阳的运动 3.1.2 绕飞阶段的观测角 3.1.3 绕飞轨道参数确定 3.2 成像卫星抵近GEO航天器 4 仿真 4.1 基于CW方程的仿真 4.2 基于ICW方程的仿真 5 结论

由于地球静止轨道（geosynchronous orbit，GEO）的轨道优势，许多高价值航天器都位于GEO区域。2021年，欧洲航天局空间碎片办公室在《欧洲年度空间环境报告》中指出，截至2020年底，共有873个GEO空间目标^[1]。GEO空间目标距离地面近36 000 km，给地面高清观测^[2-3]带来巨大挑战。

GEO区域的成像卫星能够在不受大气干扰的情况下，拍摄更清晰的GEO空间目标图像，并能更好地监测GEO航天器的工作状态。因此，部署成像卫星接近GEO航天器，实现对GEO航天器的持续观测具有重要意义。利用GEO的轨道周期特点，成像卫星可以在整个自然绕飞任务中以良好的观测角对GEO航天器进行连续观测。

根据用于表示相对运动的状态量的类型，研究相对运动的模型可以分为两类：第一类是用轨道参数作为状态量的相对运动模型；第二类是用相对状态矢量作为状态量的相对运动模型，例如用直角坐标系中的相对位置和速度^[4]作为状态矢量。第一类模型^[5-8]的优点是易于分析摄动对相对运动的影响。因此，这种模型具有较高的精度，并能有效避免偏差的长期累积。这种模型的缺点是不便于建立状态微分方程，不利于研究连续小推力作用下的相对运动。在第二类模型^[9-11]中，Clohessy-Wiltshire（CW）方程应用最为广泛，CW方程的优点在于：建立状态方程后，便于结合最优控制理论开展成像卫星接近GEO航天器的研究，也便于对飞行编队进行研究。缺点是CW方程的偏差较大，且偏差随时间不断累积^[12]。

在绕飞轨迹规划问题上，刘猛等^[13]给出了航天器基于最小冲量进入对空间目标绕飞轨道的策略，在轨道机动燃料消耗方面有借鉴价值。黄艺等^[14]研究了对非合作目标绕飞任务的姿轨耦合控制策略，给出能够对外界干扰和不确定参数具有很好抑制效果的控制方案。谭天乐等^[15-17]给出近圆轨道和大椭圆轨道的相对运动建模和解析方法，提供了对空间非合作目标高精度的相对悬停和循迹绕飞控制方法。刘涛等^[18]提出一种用于非合作目标惯性指向轴位置捕获的绕飞方法。常燕等^[19]研究了对椭圆轨道上非合作目标进行长期绕飞检测的相对运动轨道构型设计与构型保持问题。梁静静等^[20]基于粒子群算法研究了双脉冲绕飞的优化问题和安全绕飞轨迹设计问题。王功波等^[21]在连续小推力条件下，针对圆轨道推导了满足快速绕飞条件的空间圆编队动力学模型。Dang^[22-23]推导了在任意开普勒轨道上相对运动的新状态转换矩阵，给出了Tschauner-Hempel（TH）方程的解。

在CW方程的偏差问题上，通过修正CW方程的非球形摄动偏差和非线性二次项偏差改进CW方程；在轨迹规划问题上，计算对GEO航天器绕飞的初始相位区间，成像卫星在连续小推力作用下接近GEO航天器，以合适的初始相位对其开始绕飞，实现对GEO航天器的自然绕飞并持续以有利的光照条件对其观测，如图1所示。

图1成像卫星P在整个绕飞任务中以良好的观测角对GEO航天器E进行拍照

Fig.1The optical satellite P can take pictures of the GEO spacecraft E with favorable observation angles throughout the entire fly-around mission

1 CW方程

CW方程采用当地垂直当地水平（local vertical local horizontal， LVLH）坐标系，该坐标系以航天器质心为原点，x轴沿地心到航天器质心的方向；z轴垂直于轨道平面，与瞬时角动量方向相同，y轴方向根据右手定则确定。

两个非常接近的航天器分别记为M和S，M为主航天器，S为辅航天器，存在：

\frac{ρ}{a} ≪ 1

(1)

式中， ρ为两航天器之间的距离，a表示航天器M的轨道半长轴。在以航天器M为原点的LVLH参考系中，定义S的状态矢量为

X = （ x ， y ， z ， \dot{x} ， \dot{y} ， \dot{z} ）^{T}

，连续推力控制U=（u_x，u_y，u_z）^T，则CW方程为：

\dot{X} = A X + B U

(2)

式中：

A = [\begin{matrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 3 ω^{2} & 0 & 0 & 0 & 2 ω & 0 \\ 0 & 0 & 0 & - 2 ω & 0 & 0 \\ 0 & 0 & - ω^{2} & 0 & 0 & 0 \end{matrix}]

(3)

(4)

其中，ω为航天器M的轨道角速度。

记t₀时刻，S的相对状态矢量为X（t₀），则t时刻，航天器S的相对状态矢量为：

X (t, t_{0}) = ϕ (t, t_{0}) X (t_{0}) + \int_{t_{0}}^{t} ϕ (t, ε) B U (ε) d ε

(5)

式中：

\begin{matrix} ϕ (t, t_{0}) = [\begin{matrix} ϕ_{r} (t, t_{0}) & ϕ_{v} (t, t_{0}) \end{matrix}] = [\begin{matrix} ϕ_{r r} (t, t_{0}) & ϕ_{r v} (t, t_{0}) \\ ϕ_{v r} (t, t_{0}) & ϕ_{v v} (t, t_{0}) \end{matrix}] \\ = [\begin{matrix} 4 - 3 c & 0 & 0 & s / ω & - 2 (c - 1) / ω & 0 \\ 6 (s - τ) & 1 & 0 & 2 (c - 1) / ω & (4 s - 3 τ) / ω & 0 \\ 0 & 0 & c & 0 & 0 & s / ω \\ 3 ω s & 0 & 0 & c & 2 s & 0 \\ 6 ω (c - 1) & 0 & 0 & - 2 s & - 3 + 4 c & 0 \\ 0 & 0 & - ω s & 0 & 0 & c \end{matrix}] \end{matrix}

(6)

其中，φ_rr（t，t₀）、φ_rv（t，t₀）、φ_vr（t，t₀）和φ_vv（t，t₀）都是3×3矩阵，

τ = \sqrt{\frac{μ}{a^{3}}} (t - t_{0})

，s=sinτ，c=cosτ。

2 改进的CW方程

非线性、参考轨道的偏心率和摄动是CW方程偏差的三个来源。在成像卫星对GEO航天器的接近和绕飞的任务中，以GEO航天器为原点建立LVLH坐标系，则CW方程偏差来源主要为非线性项和摄动项。记P为成像卫星，E为被成像的GEO航天器。已知P和E的初始轨道参数，改进的CW（improved CW，ICW）方程研究相对运动的步骤如下：

步骤1：以航天器E的质心为原点建立LVLH坐标系，记为

R

；

步骤2：计算成像卫星P的初始状态矢量，即在LVLH坐标系中成像卫星P相对GEO航天器E的位置和速度矢量：

X (t_{0}) = {(x (t_{0}), y (t_{0}), z (t_{0}), \dot{x} (t_{0}), \dot{y} (t_{0}), \dot{z} (t_{0}))}^{T}

(7)

据文献^[22]的方法可求得：

X (t_{0}) = [\begin{matrix} A_{E C I \to L V L H} & 0_{3 \times 3} \\ {\dot{A}}_{E C I \to L V L H} & A_{E C I \to L V L H} \end{matrix}] [\begin{matrix} δ r_{E C I} \\ δ v_{E C I} \end{matrix}]

(8)

式中，A_ECI→LVLH为从地心惯性（earth-centered inertial， ECI）坐标系到LVLH坐标系的坐标变换矩阵，

{\dot{A}}_{E C I \to L V L H}

为A_ECI→LVLH的时间导数，[δr^T_ECI δv^T_ECI]^T为ECI坐标系中P与E的位置和速度矢量之差。

步骤3：GEO航天器受到的主要摄动为非球形、三体和太阳光压摄动。然而，由于GEO轨道的轨道特性，非球形摄动是导致两航天器存在相对运动偏差的最主要因素。考虑J₂和J₂₂项摄动，将CW方程中的万有引力常数μ修正为ICW方程中的万有引力常数μ′：

μ^{'} = [1 - \frac{1}{2} J_{2} {(\frac{r_{E}}{a_{c}})}^{2} + 3 J_{22} {(\frac{r_{E}}{a_{c}})}^{2} c o s (2 λ + {29.8568}^{\circ})] μ

(9)

式中，μ=3.986 005×10¹⁴ m³/s²为万有引力常数，r_E为地球的赤道半径，a_c为GEO轨道的半长轴，λ为GEO航天器E的星下点地理经度，J₂=-1.082 627×10^-3，J₂₂=1.815 528×10^-6。

步骤4：通过对CW方程初始状态矢量的修正来补偿非线性长期项偏差。将非线性二次项假设为虚加的摄动，通过虚摄动解可求解对CW方程初始状态矢量的修正项。对CW方程初始状态矢量的修正项为：

X_{n l}^{c} (t_{0}) = {[\begin{matrix} 0 0 0 0 \frac{μ}{a_{c}^{4}} \cdot \frac{ρ^{2} [2 + 3 c o s (2 β_{0})]}{8 n_{L}} \end{matrix}]}^{T}

(10)

式中，n_L为LVLH坐标系相对ECI坐标系的旋转角速度，ρ²为成像卫星P与航天器E的相对距离的平方，β₀为初始时刻y轴沿z轴反方向到成像卫星P位置矢量的角度。

步骤5：将两个独立的修正结合起来，以补偿非线性偏差和非球形摄动偏差。ICW方程的初始状态矢量为：

Y (t_{0}) = [\begin{matrix} x (t_{0}) \\ y (t_{0}) \\ z (t_{0}) \\ \dot{x} (t_{0}) \\ \dot{y} (t_{0}) + \frac{μ^{'}}{a_{c}^{4}} \cdot \frac{ρ^{2} [2 + 3 c o s (2 β_{0})]}{8 n_{L}} \\ \dot{z} (t_{0}) \end{matrix}]

(11)

求解状态方程

\dot{Y} = A Y + B U

，得到Y（t），Y（t）为LVLH坐标系中成像卫星P相对于航天器E的位置和速度。

在ICW方程中，通过步骤3补偿非球形摄动偏差，通过步骤4补偿非线性偏差，2.1节和2.2节分别研究步骤3和步骤4补偿摄动偏差和非线性偏差的原理。

2.1 补偿非球形摄动偏差

定义航天器在离地球无限远处的重力势能为0。r为ECI坐标系中的位置矢量，r为位置矢量的大小。设r处单位质量的重力势能为：

\begin{matrix} R = - \frac{μ}{r} (1 + \sum_{n = 2}^{\infty} {(\frac{r_{E}}{r})}^{n} \{J_{n} P_{n} (s i n (φ)) + \sum_{m = 1}^{n} P_{n m} (s i n (φ)) \cdot \\ [C_{n m} \cdot c o s (m λ) + S_{n m} \cdot s i n (m λ)]\}) \end{matrix}

(12)

式中，J_n为地球非球形摄动的带谐项，C_nm和S_nm为田谐项，P_n和P_nm为勒让德球谐函数，φ是矢量r与赤道平面的夹角，λ是该单位质量的星下点地理经度。归一化P_n（sin（φ））和P_nm（sin（φ）），记为

{\bar{P}}_{n}

（sin（φ））和

{\bar{P}}_{n m}

（sin（φ）），归一化系数为：

\{\begin{matrix} ζ_{n} = \sqrt{\frac{1}{2 n + 1}} \\ ζ_{n m} = \sqrt{\frac{(n + m)!}{2 (2 n + 1) (n - m)!}} m ⩽ n \end{matrix}

(13)

化简后，引力势为：

\begin{matrix} R = - \frac{μ}{r} \{1 + \sum_{n = 2}^{\infty} {(\frac{r_{E}}{r})}^{n} [J_{n} ζ_{n} {\bar{P}}_{n} (s i n (φ)) + \\ \sum_{m = 1}^{n} J_{n m} ζ_{n m} {\bar{P}}_{n m} (s i n (φ)) c o s (m λ - m λ_{n m})]\} \end{matrix}

(14)

各项的权重为：

\{\begin{matrix} A_{n} = {(\frac{r_{E}}{r})}^{n} J_{n} ζ_{n} \\ A_{n m} = {(\frac{r_{E}}{r})}^{n} J_{n m} ζ_{n m} \end{matrix}

(15)

对于GEO航天器，权量最大的五项分别是J₂，J₂₂，J₃，J₃₁和J₃₃，权重比为：

A_{2} : A_{22} : A_{3} : A_{31} : A_{33} = 11053 : 64 : 3 : 7 : 6

(16)

因此，J₂和J₂₂是非球形摄动中权重最大的项。将这两项引入ICW方程中，包含J₂和J₂₂项的引力势为：

\begin{matrix} R = - \frac{μ}{r} [1 + \frac{1}{2} J_{2} {(\frac{r_{E}}{r})}^{2} (3 {s i n}^{2} φ - 1) + \\ 3 J_{22} {(\frac{r_{E}}{r})}^{2} {c o s}^{2} φ c o s (2 λ - 2 λ_{22})] \end{matrix}

(17)

式中，λ₂₂=-14.929°。对于两个相对接近的GEO航天器，其星下点的地理经度近似相等。由于sin²φ≈0，cos²φ≈1，式（17）可简化为：

R = - \frac{μ}{r} [1 - \frac{1}{2} J_{2} {(\frac{r_{E}}{a_{c}})}^{2} + 3 J_{22} {(\frac{r_{E}}{a_{c}})}^{2} c o s (2 λ + {29.8568}^{\circ})]

(18)

J₂和J₂₂对相对运动的影响体现在与CW方程相比，ICW方程中的万有引力常数修正为：

μ^{'} = [1 - \frac{1}{2} J_{2} {(\frac{r_{E}}{a_{c}})}^{2} + 3 J_{22} {(\frac{r_{E}}{a_{c}})}^{2} c o s (2 λ + {29.8568}^{\circ})] μ

(19)

2.2 补偿非线性偏差

在CW方程中，为了得到线性的状态方程，重力加速度的泰勒级数展开式中的非线性项在式（1）的条件下被舍弃，二次项是偏差的主要来源。被舍弃的二次项在位置分量上产生常数项、长期项和短周期项三种类型的偏差。长期项只出现在y轴方向上，是对精度影响最大的偏差。偏差的长期项可通过对CW方程的初始状态矢量的修正来补偿，利用摄动法可计算用于补偿CW方程中非线性偏差的修正项。摄动法中，将非线性二次项作为CW方程中的虚加摄动力。通过虚摄动解求解用于补偿CW方程中非线性长期项偏差的初始状态矢量的修正项。

航天器在ECI坐标系中的瞬时位置和速度矢量分别为r和

\dot{r}

，根据牛顿第二定律得：

\{\begin{matrix} \frac{d r}{d t} = \dot{r} \\ \frac{d \dot{r}}{d t} = - g r a d (R) \end{matrix}

(20)

其中，R为在r处单位质量的势能。成像卫星P相对GEO航天器E的位置矢量为ρ=r_P-r_c。根据式（12），对于航天器E和成像卫星P，存在：

{\ddot{r}}_{c} = - g r a d (R_{c}) = - \frac{μ}{r_{c}^{3}} r_{c} + f_{c}

(21)

{\ddot{r}}_{P} = - g r a d (R_{P}) = - \frac{μ}{r_{P}^{3}} r_{P} + f_{P}

(22)

其中， f_c和f_P分别为航天器E和成像卫星P受到的摄动力，如太阳光压摄动等。式（21）与式（22）相减得：

{\ddot{r}}_{P} - {\ddot{r}}_{c} = - (\frac{μ}{r_{P}^{3}} r_{P} - \frac{μ}{r_{c}^{3}} r_{c}) + Δ f

(23)

式中，Δf为成像卫星P和航天器E受到的摄动偏差。ICW方程在步骤3中对主要摄动偏差进行修正，在步骤4中对非线性偏差进行修正。因此，在修正非线性偏差时，令Δf=0。将式（23）从ECI坐标系转换到LVLH坐标系：

\ddot{ρ} + {\dot{ω}}_{c} \times ρ + 2 ω_{c} \times \dot{ρ} + ω_{c} \times (ω_{c} \times ρ) = - (\frac{μ}{r_{i}^{3}} r_{i} - \frac{μ}{r_{c}^{3}} r_{c})

(24)

式中，ω_c=[0 0 n]^T，

{\dot{ω}}_{c}

=[0 0 0]^T，r_c=[r_c 0 0]^T，

n = \sqrt{\frac{μ}{a_{c}^{3}}}

。式（24）可整理为：

\{\begin{matrix} \ddot{x} - 2 n \dot{y} - n^{2} x = - \frac{μ (r_{c} + x)}{{[{(r_{c} + x)}^{2} + y^{2} + z^{2}]}^{3 / 2}} + \frac{μ}{r_{c}^{2}} \\ \ddot{y} + 2 n \dot{x} - n^{2} y = - \frac{μ y}{{[{(r_{c} + x)}^{2} + y^{2} + z^{2}]}^{3 / 2}} \\ \ddot{z} = - \frac{μ z}{{[{(r_{c} + x)}^{2} + y^{2} + z^{2}]}^{3 / 2}} \end{matrix}

(25)

在CW方程中，式（25）可化简为：

\{\begin{matrix} \ddot{x} - 2 n \dot{y} - 3 n^{2} x = 0 \\ \ddot{y} + 2 n \dot{x} = 0 \\ \ddot{z} + n^{2} z = 0 \end{matrix}

(26)

记CW方程中的状态矢量为：

X_{c w} (t) = [\begin{matrix} x_{c w} (t) y_{c w} (t) z_{c w} (t) {\dot{x}}_{c w} (t) {\dot{y}}_{c w} (t) {\dot{z}}_{c w} (t) \end{matrix}]

(27)

解为：

\{\begin{matrix} x_{c w} = \frac{ρ}{2} s i n (n t + β_{0}) \\ y_{c w} = ρ c o s (n t + β_{0}) \\ z_{c w} = \frac{\sqrt{3} ρ}{2} s i n (n t + β_{0}) \\ {\dot{x}}_{c w} = \frac{ρ}{2} n c o s (n t + β_{0}) \\ {\dot{y}}_{c w} = - ρ n s i n (n t + β_{0}) \\ {\dot{z}}_{c w} = \frac{\sqrt{3} ρ}{2} n c o s (n t + β_{0}) \end{matrix}

(28)

保留式（25）等号右边的二次项，舍弃高阶项。将式（25）化简为：

\{\begin{matrix} \ddot{x} - 2 n \dot{y} - 3 n^{2} x = \frac{3 μ}{a_{c}^{4}} (\frac{y^{2}}{2} + \frac{z^{2}}{2} - x^{2}) \\ \ddot{y} + 2 n \dot{x} = \frac{3 μ}{a_{c}^{4}} x y \\ \ddot{z} + n^{2} z = \frac{3 μ}{a_{c}^{4}} x z \end{matrix}

(29)

根据式（29），将重力加速度泰勒级数展开式中的非线性二次项作为CW方程中虚加的摄动力，虚加摄动力为：

U_{n l} = {[\begin{matrix} \frac{3 μ}{a_{c}^{4}} (\frac{y^{2}}{2} + \frac{z^{2}}{2} - x^{2}) \frac{3 μ}{a_{c}^{4}} x y \frac{3 μ}{a_{c}^{4}} x z \end{matrix}]}^{T}

(30)

将式（29）简记为：

\dot{X} = A X + B U_{n l}

(31)

式（31）为CW方程的非线性偏差传播方程。记CW方程的非线性偏差为：

X_{n l} (t) = {[\begin{matrix} x_{n l} (t) y_{n l} (t) z_{n l} (t) {\dot{x}}_{n l} (t) {\dot{y}}_{n l} (t) \end{matrix} {\dot{z}}_{n l} (t)]}^{T}

(32)

根据式（5）可得式（31）的解为：

X_{n l} (t) = ϕ (t, t_{0}) X_{n l} (t_{0}) + \int_{t_{0}}^{t} ϕ (t, s) B U_{n l}^{c} (s) d s

(33)

CW方程中被舍弃的二次项导致y轴方向偏差的长期项为：

\begin{matrix} y_{n l - long} (t) = - \{6 n x_{n l} (t_{0}) + 3 {\dot{y}}_{n l} (t_{0}) - \\ \frac{3}{8} \frac{μ}{a_{c}^{4}} \frac{ρ^{2} [2 + 3 c o s (2 β_{0})]}{n}\} t \end{matrix}

(34)

为了修正y轴方向偏差的长期项，令y_nl_-long=0，可得非线性偏差的零累积准则：

\{\begin{matrix} x_{n l}^{c} (t_{0}) = 0 \\ {\dot{y}}_{n l}^{c} (t_{0}) = \frac{μ}{a_{c}^{4}} \frac{ρ^{2} [2 + 3 c o s (2 β_{0})]}{8 n} \end{matrix}

(35)

因此，根据虚摄动解可得到对CW方程初始状态矢量的修正项如式（10）所示。

通过对初始状态矢量X_cw（t₀）添加修正项X^c_nl（t₀），可以补偿CW方程中的非线性长期偏差，即在ICW方程中，初始状态矢量X_cw（t₀）被替换为X_cw（t₀）+X^c_nl（t₀）。

2.3 偏差修正验证

为了研究ICW方程相对CW方程的改进效果，选取成像卫星P和GEO航天器E进行仿真验证，其轨道半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点幅角w、平近地点角m如表1所示。

表1成像卫星P与GEO航天器E的轨道参数

Tab.1 Orbital elements of the optical satellite P and the GEO spacecraft E

LVLH坐标系中成像卫星P的初始状态矢量如表2所示。初始时刻从y轴沿z轴反方向到成像卫星P的位置矢量的角度为0°，X为CW方程的初始状态矢量，Y为ICW方程的初始状态矢量。

2.3.1 非球形摄动偏差修正验证

为了验证ICW方程中步骤3的非球形摄动解，将非球形摄动解与高精度轨道传播（high-precision orbit propagator， HPOP）模型作对比。HPOP模型中考虑前21阶非球形摄动、太阳光压摄动和三体引力摄动。利用二体模型、HPOP模型和步骤3中的非球形摄动解，分别在ECI坐标系中计算成像卫星P和GEO航天器E的位置。将P在ECI坐标系中的位置分别转换到LVLH坐标系

R

中，P在

R

中的位置即为P相对于E的位置，将HPOP模型得到的相对位置视为真实相对位置。在仿真场景中，航天器星下点的地理经度为165″E。

表2成像卫星P的初始状态矢量

Tab.2 The initial state vector of the optical satellite P

相比二体模型，非球形摄动解在相对位置的x、y和z方向对精度的提高分别如图2、图3和图4所示。图中红色实线表示非球形摄动解相对HPOP模型的偏差，蓝色虚线表示二体模型相对HPOP模型的偏差。通过比较可得，非球形摄动是摄动偏差的最主要因素，非球形摄动解提高了CW方程在LVLH坐标系x、y轴方向的精度。

图2非球形摄动解对x方向精度的提高

Fig.2Improvement of the non-spherical perturbation solution in x-coordinate

图3非球形摄动解对y方向精度的提高

Fig.3Improvement of the non-spherical perturbation solution in y-coordinate

图4非球形摄动解对z方向精度的提高

Fig.4Improvement of the non-spherical perturbation solution in z-coordinate

利用ICW方程和CW方程仿真成像卫星P相对于GEO航天器E的运动，比较三个轨道周期内相对位置的x、y、z坐标差异。由于CW方程存在非线性偏差，为了展现ICW方程中步骤3对CW方程的改进效果，本节中不对ICW方程进行非线性偏差修正。ICW方程与CW方程中相对位置的x、y、z坐标差值分别如图5、图6和图7所示。

图5ICW方程和CW方程在x轴方向的差值

Fig.5Deviation in x-coordinate between ICW equations and CW equations

从图5~7中可以看出，ICW方程与CW方程中相对位置的x、y、z坐标差值呈周期性变化，且差值幅值逐渐增大。x、y坐标差值幅值的增长率在同一个数量级，z坐标幅值的增长率比x、y坐标幅值的增长率小一个数量级。

图6ICW方程和CW方程在y轴方向的差值

Fig.6Deviation in y-coordinate between ICW equations and CW equations

图7ICW方程和CW方程在z轴方向的差值

Fig.7Deviation in z-coordinate between ICW equations and CW equations

2.3.2 非线性偏差修正验证

CW方程和ICW方程计算出的成像卫星P相对GEO航天器E的位置的y坐标如图8所示。图中红色实线代表ICW方程的y坐标，蓝色虚线代表CW方程的y坐标。ICW方程与CW方程的y坐标差值如图9所示。

为了验证ICW方程步骤4的改进效果，将CW方程和ICW方程计算的相对位置的y坐标与真值进行比较。由于CW方程存在摄动偏差，为了展现ICW方程中步骤4对CW方程的改进效果，本节不对ICW方程进行摄动偏差修正。如图10所示，红色实线表示ICW方程中相对位置的y坐标与真值的差值，蓝色虚线表示CW方程中相对位置的y坐标与真值的差值。

成像卫星P与GEO航天器E的初始距离为200 km，3个轨道周期后，ICW方程中相对位置的y坐标相对真值出现36 m的长期偏差，而CW方程的长期累积偏差为33 563 m。ICW方程的长期偏差仅约为CW方程的0.1%。通过比较，验证了ICW方程对非线性偏差长期项补偿的有效性。

图8ICW方程和CW方程中的y坐标值

Fig.8y-coordinate in ICW equations and CW equations

图9ICW方程和CW方程在y轴方向的非线性偏差值

Fig.9y-coordinate nonlinear deviation between ICW equations and CW equations

图10ICW方程和CW方程在y轴方向与真值的差值

Fig.10y-coordinate deviation of ICW equations and CW equations from the truth value

ICW方程补偿了相对位置y轴方向二次非线性偏差的长期项和x、y轴方向上主要的摄动偏差。ICW方程中，在建立状态微分方程时，不对非线性偏差的短期项进行修正。在有更高的精度需求时，可用短期项偏差修正结果状态矢量。这样处理的优点在于：得到的控制系统是线性时不变系统，同时在最终的状态矢量中加入短周期项可保证结果的准确性。

3 轨迹规划

对GEO航天器绕飞持续观测任务的轨迹规划分为两段：成像卫星抵近GEO航天器和成像卫星对GEO航天器绕飞。为了实现燃料消耗最少，成像卫星最理想的轨迹规划方案是：在连续小推力控制下接近GEO航天器，然后对GEO航天器自然绕飞完成持续观测任务。

3.1 成像卫星对GEO航天器绕飞

3.1.1 观测角和太阳的运动

观测角定义如图11所示。定义α为成像卫星对GEO航天器的观测角。0°观测角表示：成像卫星、GEO航天器和太阳共线，成像卫星处于中间位置。显然，越接近0°，成像卫星对GEO航天器的观测越有利；越接近180°，成像卫星对GEO航天器的观测越不利。成像卫星为了在整个绕飞阶段都能对GEO航天器成像，需观测角始终小于60°。成像卫星能较好地对GEO航天器成像的观测角的最大值取决于成像卫星相机的成像能力和GEO航天器的材质，通常情况下，成像卫星在观测角小于60°的情况下，可实现对GEO航天器成像。

图11观测角的定义

Fig.11Observation angle

在ECI坐标系中，太阳的轨道参数为：

\{\begin{matrix} a = 1.00000102 天 文 单 位 \\ e = 0.01670862 - 0.00004204 T - 0.00000124 T^{2} \\ i = {23.439291}^{\circ} - {0.01300417}^{''} T - {0.00000016}^{''} T^{2} \\ Ω = 0^{''} \\ ω = {282.937347}^{''} + {0.32256206}^{''} T - {0.00015757}^{''} T^{2} \\ m = {357.529100}^{\circ} + {0.98556200804}^{\circ} d - {0.0007734}^{\circ} d^{2} \end{matrix}

(36)

其中，T为儒略世纪，d为儒略日。据太阳的轨道参数，可得太阳在ECI坐标系中的位置矢量。由于太阳到地心的距离远大于GEO航天器到地心的距离，假设地心到太阳的方向与GEO航天器E到太阳的方向重合。太阳在以航天器E为原点的LVLH坐标系

R

中的方向矢量为：

d_{s u n} = A_{E C I \to L V L H} \cdot k_{E C I}

(37)

式中，k_ECI为太阳在ECI中的方向矢量。

3.1.2 绕飞阶段的观测角

设计成像卫星P在GEO航天器E轨道平面内的绕飞轨迹。设t_p时刻太阳方向矢量在LVLH坐标系

R

的xoy平面上的投影与x轴重合，成像卫星P在t₀到t_p时刻完成对航天器E的抵近，从t_p时刻开始绕飞。d_P为航天器E到成像卫星P的位置矢量，记θ_xy为d_P在t_p时刻的初始相位角，即θ_xy为t_p时刻x轴沿z轴反方向到d_P矢量的角度，如图12所示。δ为太阳的赤纬，即太阳的星下点纬度，显然δ∈（-23°26′，+23°26′）。t时刻，d_P在LVLH坐标系

R

中的x、y、z坐标为：

d_{P} = [x (t), y (t), z (t)]

(38)

\{\begin{matrix} x (t) = e_{A} a_{c} c o s (- n_{E} (t - t_{p}) - θ_{x y}) \\ y (t) = 2 e_{A} a_{c} s i n (- n_{E} (t - t_{p}) - θ_{x y}) \\ z (t) = 0 \end{matrix}

(39)

式中，

e_{A} = \sqrt{{(e_{c} + Δ e)}^{2} + e_{c}^{2} - 2 e_{c} (e_{c} + Δ e) c o s （ Δ m ）}

，n_E为ECI坐标系中航天器E的轨道角速度。成像卫星对GEO航天器的观测角为：

图12太阳和成像卫星的相对位置

Fig.12Position of the sun and the optical satellite

α = ⟨d_{P}, d_{sun}⟩ = a r c c o s \frac{d_{P} \cdot d_{sun}}{|d_{P}| \cdot |d_{sun}|}

(40)

当θ_xy=0°时，成像卫星P在GEO航天器E的一个轨道周期内对E的观测角如图13所示。

图13在θ_xy=0°时成像卫星P对GEO航天器E的观测角

Fig.13Observation angles of the optical satellite P to the GEO spacecraft E when θ_xy=0°

太阳的赤纬在-23°26′到+23°26′间变化。记α_m为一个轨道周期内P对E的最大观测角，α_m随太阳赤纬而变化，如图14所示。当θ_xy=0°时，α_m在δ=0°时取得最小值19.47°，在δ=±23°26′时取得最大值26.35°。因此，成像卫星P在t_p时刻以θ_xy=0开始对GEO航天器E绕飞，成像卫星P对GEO航天器E的观测角始终小于26.35°。

图14成像卫星P对GEO航天器E的最大观测角随太阳赤纬的变化

Fig.14Maximum observation angle of the optical satellite P to the GEO spacecraft E varies with the declination of the sun δ

t_p时刻，成像卫星P以初始相位角θ_xy≠0°开始对GEO航天器E绕飞。θ_xy取不同值时，一个轨道周期内的最大观测角α_m随太阳赤纬δ变化，如图15所示。因此，无论太阳的赤纬δ取何值，如果绕飞的初始相位θ_xy∈（-37.2°，37.2°），在整个绕飞任务中，P对E的观测角始终小于60°，P可以以有利的观测角对E拍照。因此，可对GEO航天器持续观测的成像卫星的绕飞初始相位角区间为θ_xy∈（-37.2°，37.2°）。

图15不同θ_xy时最大观测角α_m随太阳赤纬的变化

Fig.15The maximum observation angle in an orbital period α_m varies with the declination of the sun δ at different θ_xy

3.1.3 绕飞轨道参数确定

设计成像卫星P在GEO航天器E轨道平面内的绕飞。在LVLH坐标系

R

中，绕飞相对轨迹为椭圆，记r_xy为相对轨迹的半长轴，记y_off为椭圆相对轨道的中心与航天器E的距离。记Δe、Δi、ΔΩ、Δω和Δm为P与E轨道的偏心率、轨道倾角、升交点赤经、近地点幅角和平近地点角差值：

\{\begin{matrix} Δ a = 0 m \\ Δ e = \sqrt{{(\frac{r_{x y}}{2 a_{c}})}^{2} + \frac{e_{c} r_{x y}}{a_{c}} c o s θ_{x y} + e_{c}^{2}} - e_{c} \\ Δ i = 0^{\circ} \\ Δ Ω = 0^{\circ} \\ Δ ω = \frac{y_{o f f}}{a_{c}} - a r c t a n \frac{r_{x y} s i n θ_{x y}}{r_{x y} c o s θ_{x y} + 2 a_{c} e_{c}} \\ Δ m = a r c t a n \frac{r_{x y} s i n θ_{x y}}{r_{x y} c o s θ_{x y} + 2 a_{c} e_{c}} \end{matrix}

(41)

已知GEO航天器E的轨道参数，根据任务要求的r_xy、θ_xy和y_off，设计成像卫星P的相对轨迹。根据式（41）可得Δa、Δe、Δi、ΔΩ、Δω、Δm，进而可得成像卫星P在t_p时刻的轨道参数。

3.2 成像卫星抵近GEO航天器

据3.1节的方法，可得成像卫星P在t_p时刻的轨道参数，已知成像卫星在t₀时刻的轨道参数。以成像卫星P的燃料消耗为优化指标，计算成像卫星接近GEO航天器的控制策略，从t₀到t_p的轨迹优化是一个两点边界问题。状态方程为：

\dot{Y} = A Y + B U

(42)

支付函数为：

J = \frac{1}{2} \int_{t_{0}}^{t_{f}} (U^{T} U) d t

(43)

根据式（43），哈密顿函数可写为：

H = \frac{1}{2} U^{T} U + λ^{T} (A Y + B U)

(44)

式中，λ是协态变量。

根据最优化的必要条件，协态方程为：

{\dot{λ}}^{*} = - \frac{\partial H}{\partial Y} = - A^{T} λ^{*}

(45)

\frac{\partial H}{\partial U} = U^{*} + B^{T} λ^{*} = 0

(46)

把式（46）代入式（42），得：

{\dot{Y}}^{*} = A Y^{*} - B B^{T} λ^{*}

(47)

据式（47）和式（45），可得：

[\begin{matrix} {\dot{Y}}^{*} \\ {\dot{λ}}^{*} \end{matrix}] = [\begin{matrix} A & - B B^{T} \\ 0 & - A^{T} \end{matrix}] [\begin{matrix} Y^{*} \\ λ^{*} \end{matrix}]

(48)

式（48）的解析解为：

[\begin{matrix} Y^{*} (t) \\ λ^{*} (t) \end{matrix}] = e x p \{[\begin{matrix} A & - B B^{T} \\ 0 & - A^{T} \end{matrix}] (t - t_{0})\} [\begin{matrix} Y (t_{0}) \\ λ^{*} (t_{0}) \end{matrix}]

(49)

式（49）可写为：

[\begin{matrix} Y^{*} (t) \\ λ^{*} (t) \end{matrix}] = [\begin{matrix} Π_{11} (t) Π_{12} (t) \\ Π_{21} (t) Π_{22} (t) \end{matrix}] [\begin{matrix} Y (t_{0}) \\ λ^{*} (t_{0}) \end{matrix}]

(50)

已知ICW方程中的初始状态矢量Y（t₀）和终端状态矢量Y（t_p），可得：

λ^{*} (t_{0}) = Π_{12} {(t_{p})}^{- 1} [Y^{*} (t_{p}) - Π_{11} (t_{p}) Y (t_{0})]

(51)

根据式（46）、式（50）和式（51），可得最优控制策略：

U^{*} (t) = - B^{T} [Π_{21} (t) Y (t_{0}) + Π_{22} (t) λ^{*} (t_{0})]

(52)

4 仿真

从t₀到t_p，成像卫星P在连续小推力的作用下接近GEO航天器E；自t_p时刻开始，成像卫星P对GEO航天器自然绕飞。成像卫星P的相对绕飞轨迹在GEO航天器E的轨道平面内，椭圆相对轨迹的半长轴为20 km，初始相位角θ_xy=0°。在该场景中，t₀为GEO航天器星下点当地时间2021年6月10日00:00:00，t_p为2021年6月11日12:00:00。成像卫星P和GEO航天器E在t₀时刻的轨道参数如表1所示。

根据CW方程和ICW方程，LVLH坐标系

R

中成像卫星P在t₀和t_p时的状态矢量如表3所示，X为CW方程的状态矢量，Y为ICW方程的状态矢量。

表3成像卫星P在LVLH坐标系中的状态矢量

Tab.3 State vector of the optical satellite P in the LVLH frame

分别基于CW方程和ICW方程进行仿真。计算成像卫星P的最优控制策略，将成像卫星P在ECI坐标系中的运动方程积分得到其位置和速度，再将ECI中得到的位置转换到LVLH坐标系中。积分采用龙格-库塔法，积分得到的结果作为成像卫星P的真实位置，进而可得成像卫星P相对GEO航天器E的位置以及P对E绕飞过程中的观测角。将基于CW方程和ICW方程得到的相对位置和观测角进行对比。

4.1 基于CW方程的仿真

基于CW方程进行仿真。t₀到t_p时间内，成像卫星P在连续小推力作用下接近GEO航天器E，施加的总速度增量为4.70 m/s。图16为成像卫星P相对GEO航天器E的位置在LVLH坐标系xoy平面中的投影。由于CW方程存在偏差，相对绕飞轨迹的中心在y轴方向上偏离GEO航天器E。图17为P到E的距离，t_p时刻，P到E的距离不是预先规划的10 km，而是14 km。图18为成像卫星P在t_p后的一个轨道周期内对GEO航天器E的观测角，最大观测角为75.57°，成像观测仿真任务失败。

4.2 基于ICW方程的仿真

基于ICW方程进行仿真，图19为成像卫星P从t₀到t_p的控制推力加速度（工程上一般表述为控制推力），总速度增量为4.67 m/s。图20~22分别表示成像卫星P的推力控制的x、y、z坐标分量。图23为P相对于E的位置，图24为P与E间的距离，t_p时刻P与E的距离为10 km，P开始对E自然绕飞。图25为P对E在t_p后一个轨道周期内的观测角，最大观测角为22.6°。成像卫星P能以良好的观测角在整个绕飞任务中对GEO航天器E成像。

图16基于CW方程仿真的成像卫星P相对GEO 航天器E的位置在LVLH坐标系xoy平面中的投影

Fig.16Position of the optical satellite P relative to the GEO spacecraft E on the xoy plane based on CW equations simulation

图17基于CW方程仿真的成像卫星P与 GEO航天器E的距离

Fig.17Distance between the optical satellite P and the GEO spacecraft E, simulation based on CW equations

图18基于CW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.18Observation angles of the optical satellite P to the GEO spacecraft E based on CW equations simulation

图19基于ICW方程仿真的成像卫星P的控制推力

Fig.19Thrust control of the optical satellite P based on ICW equations simulation

图20基于ICW方程仿真的成像卫星P在x方向的推力

Fig.20Thrust control of the optical satellite P in x-coordinate simulation based on ICW equations

图21基于ICW方程仿真的成像卫星P在y方向的推力

Fig.21Thrust control of the optical satellite P in y-coordinate simulation based on ICW equations

图22基于ICW方程仿真的成像卫星P在z方向的推力

Fig.22Thrust control of the optical satellite P in z-coordinate simulation based on ICW equations

图23基于ICW方程仿真的成像卫星P相对 GEO航天器E的位置

Fig.23Position of the optical satellite P relative to the GEO spacecraft E based on ICW equations simulation

图24基于ICW方程仿真的成像卫星P与 GEO航天器E间的距离

Fig.24Distance between the optical satellite P and the GEO spacecraft E based on ICW equations simulation

图25基于ICW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.25Observation angles of the optical satellite P to the GEO spacecraft E based on ICW equations simulation

5 结论

通过修正非球形摄动偏差和重力加速度二次长期项偏差改进CW方程。将CW方程、ICW方程与真值作仿真对比，ICW方程补偿了非线性偏差的长期项和主要的摄动偏差。在轨迹规划问题上，计算对GEO航天器绕飞的初始相位角区间，成像卫星以该区间内的初始相位角对GEO航天器开始绕飞，可实现在整个绕飞任务中都能以良好的观测角对GEO航天器进行观测。选择合适的初始相位角，分别基于CW方程和ICW方程仿真接近和绕飞的全过程。在仿真中，成像卫星在连续的小推力作用下接近GEO航天器，然后对GEO航天器自然绕飞。基于CW方程的仿真任务失败，而基于ICW方程达到预期目标。在基于ICW方程的仿真中，所需要的总速度增量仅为4.67 m/s，工程上具有很强的可行性。ICW方程可满足精度要求，可用于GEO空间态势感知和在轨服务任务规划。

图1成像卫星P在整个绕飞任务中以良好的观测角对GEO航天器E进行拍照

Fig.1The optical satellite P can take pictures of the GEO spacecraft E with favorable observation angles throughout the entire fly-around mission

下载: 全尺寸图片

图2非球形摄动解对x方向精度的提高

Fig.2Improvement of the non-spherical perturbation solution in x-coordinate

下载: 全尺寸图片

图3非球形摄动解对y方向精度的提高

Fig.3Improvement of the non-spherical perturbation solution in y-coordinate

下载: 全尺寸图片

图4非球形摄动解对z方向精度的提高

Fig.4Improvement of the non-spherical perturbation solution in z-coordinate

下载: 全尺寸图片

图5ICW方程和CW方程在x轴方向的差值

Fig.5Deviation in x-coordinate between ICW equations and CW equations

下载: 全尺寸图片

图6ICW方程和CW方程在y轴方向的差值

Fig.6Deviation in y-coordinate between ICW equations and CW equations

下载: 全尺寸图片

图7ICW方程和CW方程在z轴方向的差值

Fig.7Deviation in z-coordinate between ICW equations and CW equations

下载: 全尺寸图片

图8ICW方程和CW方程中的y坐标值

Fig.8y-coordinate in ICW equations and CW equations

下载: 全尺寸图片

图9ICW方程和CW方程在y轴方向的非线性偏差值

Fig.9y-coordinate nonlinear deviation between ICW equations and CW equations

下载: 全尺寸图片

图10ICW方程和CW方程在y轴方向与真值的差值

Fig.10y-coordinate deviation of ICW equations and CW equations from the truth value

下载: 全尺寸图片

图11观测角的定义

Fig.11Observation angle

下载: 全尺寸图片

图12太阳和成像卫星的相对位置

Fig.12Position of the sun and the optical satellite

下载: 全尺寸图片

图13在θ_xy=0°时成像卫星P对GEO航天器E的观测角

Fig.13Observation angles of the optical satellite P to the GEO spacecraft E when θ_xy=0°

下载: 全尺寸图片

图14成像卫星P对GEO航天器E的最大观测角随太阳赤纬的变化

Fig.14Maximum observation angle of the optical satellite P to the GEO spacecraft E varies with the declination of the sun δ

下载: 全尺寸图片

图15不同θ_xy时最大观测角α_m随太阳赤纬的变化

Fig.15The maximum observation angle in an orbital period α_m varies with the declination of the sun δ at different θ_xy

下载: 全尺寸图片

图16基于CW方程仿真的成像卫星P相对GEO 航天器E的位置在LVLH坐标系xoy平面中的投影

Fig.16Position of the optical satellite P relative to the GEO spacecraft E on the xoy plane based on CW equations simulation

下载: 全尺寸图片

图17基于CW方程仿真的成像卫星P与 GEO航天器E的距离

Fig.17Distance between the optical satellite P and the GEO spacecraft E, simulation based on CW equations

下载: 全尺寸图片

图18基于CW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.18Observation angles of the optical satellite P to the GEO spacecraft E based on CW equations simulation

下载: 全尺寸图片

图19基于ICW方程仿真的成像卫星P的控制推力

Fig.19Thrust control of the optical satellite P based on ICW equations simulation

下载: 全尺寸图片

图20基于ICW方程仿真的成像卫星P在x方向的推力

Fig.20Thrust control of the optical satellite P in x-coordinate simulation based on ICW equations

下载: 全尺寸图片

图21基于ICW方程仿真的成像卫星P在y方向的推力

Fig.21Thrust control of the optical satellite P in y-coordinate simulation based on ICW equations

下载: 全尺寸图片

图22基于ICW方程仿真的成像卫星P在z方向的推力

Fig.22Thrust control of the optical satellite P in z-coordinate simulation based on ICW equations

下载: 全尺寸图片

图23基于ICW方程仿真的成像卫星P相对 GEO航天器E的位置

Fig.23Position of the optical satellite P relative to the GEO spacecraft E based on ICW equations simulation

下载: 全尺寸图片

图24基于ICW方程仿真的成像卫星P与 GEO航天器E间的距离

Fig.24Distance between the optical satellite P and the GEO spacecraft E based on ICW equations simulation

下载: 全尺寸图片

图25基于ICW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.25Observation angles of the optical satellite P to the GEO spacecraft E based on ICW equations simulation

下载: 全尺寸图片

表1成像卫星P与GEO航天器E的轨道参数

下载: 全尺寸图片

表2成像卫星P的初始状态矢量

下载: 全尺寸图片

表3成像卫星P在LVLH坐标系中的状态矢量

下载: 全尺寸图片

图1成像卫星P在整个绕飞任务中以良好的观测角对GEO航天器E进行拍照

Fig.1The optical satellite P can take pictures of the GEO spacecraft E with favorable observation angles throughout the entire fly-around mission

图2非球形摄动解对x方向精度的提高

Fig.2Improvement of the non-spherical perturbation solution in x-coordinate

图3非球形摄动解对y方向精度的提高

Fig.3Improvement of the non-spherical perturbation solution in y-coordinate

图4非球形摄动解对z方向精度的提高

Fig.4Improvement of the non-spherical perturbation solution in z-coordinate

图5ICW方程和CW方程在x轴方向的差值

Fig.5Deviation in x-coordinate between ICW equations and CW equations

图6ICW方程和CW方程在y轴方向的差值

Fig.6Deviation in y-coordinate between ICW equations and CW equations

图7ICW方程和CW方程在z轴方向的差值

Fig.7Deviation in z-coordinate between ICW equations and CW equations

图8ICW方程和CW方程中的y坐标值

Fig.8y-coordinate in ICW equations and CW equations

图9ICW方程和CW方程在y轴方向的非线性偏差值

Fig.9y-coordinate nonlinear deviation between ICW equations and CW equations

图10ICW方程和CW方程在y轴方向与真值的差值

Fig.10y-coordinate deviation of ICW equations and CW equations from the truth value

图11观测角的定义

Fig.11Observation angle

图12太阳和成像卫星的相对位置

Fig.12Position of the sun and the optical satellite

图13在θ_xy=0°时成像卫星P对GEO航天器E的观测角

Fig.13Observation angles of the optical satellite P to the GEO spacecraft E when θ_xy=0°

图14成像卫星P对GEO航天器E的最大观测角随太阳赤纬的变化

Fig.14Maximum observation angle of the optical satellite P to the GEO spacecraft E varies with the declination of the sun δ

图15不同θ_xy时最大观测角α_m随太阳赤纬的变化

Fig.15The maximum observation angle in an orbital period α_m varies with the declination of the sun δ at different θ_xy

图16基于CW方程仿真的成像卫星P相对GEO 航天器E的位置在LVLH坐标系xoy平面中的投影

Fig.16Position of the optical satellite P relative to the GEO spacecraft E on the xoy plane based on CW equations simulation

图17基于CW方程仿真的成像卫星P与 GEO航天器E的距离

Fig.17Distance between the optical satellite P and the GEO spacecraft E, simulation based on CW equations

图18基于CW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.18Observation angles of the optical satellite P to the GEO spacecraft E based on CW equations simulation

图19基于ICW方程仿真的成像卫星P的控制推力

Fig.19Thrust control of the optical satellite P based on ICW equations simulation

图20基于ICW方程仿真的成像卫星P在x方向的推力

Fig.20Thrust control of the optical satellite P in x-coordinate simulation based on ICW equations

图21基于ICW方程仿真的成像卫星P在y方向的推力

Fig.21Thrust control of the optical satellite P in y-coordinate simulation based on ICW equations

图22基于ICW方程仿真的成像卫星P在z方向的推力

Fig.22Thrust control of the optical satellite P in z-coordinate simulation based on ICW equations

图23基于ICW方程仿真的成像卫星P相对 GEO航天器E的位置

Fig.23Position of the optical satellite P relative to the GEO spacecraft E based on ICW equations simulation

图24基于ICW方程仿真的成像卫星P与 GEO航天器E间的距离

Fig.24Distance between the optical satellite P and the GEO spacecraft E based on ICW equations simulation

图25基于ICW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.25Observation angles of the optical satellite P to the GEO spacecraft E based on ICW equations simulation

表1成像卫星P与GEO航天器E的轨道参数

表2成像卫星P的初始状态矢量

表3成像卫星P在LVLH坐标系中的状态矢量

图1成像卫星P在整个绕飞任务中以良好的观测角对GEO航天器E进行拍照

Fig.1The optical satellite P can take pictures of the GEO spacecraft E with favorable observation angles throughout the entire fly-around mission

图2非球形摄动解对x方向精度的提高

Fig.2Improvement of the non-spherical perturbation solution in x-coordinate

图3非球形摄动解对y方向精度的提高

Fig.3Improvement of the non-spherical perturbation solution in y-coordinate

图4非球形摄动解对z方向精度的提高

Fig.4Improvement of the non-spherical perturbation solution in z-coordinate

图5ICW方程和CW方程在x轴方向的差值

Fig.5Deviation in x-coordinate between ICW equations and CW equations

图6ICW方程和CW方程在y轴方向的差值

Fig.6Deviation in y-coordinate between ICW equations and CW equations

图7ICW方程和CW方程在z轴方向的差值

Fig.7Deviation in z-coordinate between ICW equations and CW equations

图8ICW方程和CW方程中的y坐标值

Fig.8y-coordinate in ICW equations and CW equations

图9ICW方程和CW方程在y轴方向的非线性偏差值

Fig.9y-coordinate nonlinear deviation between ICW equations and CW equations

图10ICW方程和CW方程在y轴方向与真值的差值

Fig.10y-coordinate deviation of ICW equations and CW equations from the truth value

图11观测角的定义

Fig.11Observation angle

图12太阳和成像卫星的相对位置

Fig.12Position of the sun and the optical satellite

图13在θ_xy=0°时成像卫星P对GEO航天器E的观测角

Fig.13Observation angles of the optical satellite P to the GEO spacecraft E when θ_xy=0°

图14成像卫星P对GEO航天器E的最大观测角随太阳赤纬的变化

Fig.14Maximum observation angle of the optical satellite P to the GEO spacecraft E varies with the declination of the sun δ

图15不同θ_xy时最大观测角α_m随太阳赤纬的变化

Fig.15The maximum observation angle in an orbital period α_m varies with the declination of the sun δ at different θ_xy

图16基于CW方程仿真的成像卫星P相对GEO 航天器E的位置在LVLH坐标系xoy平面中的投影

Fig.16Position of the optical satellite P relative to the GEO spacecraft E on the xoy plane based on CW equations simulation

图17基于CW方程仿真的成像卫星P与 GEO航天器E的距离

Fig.17Distance between the optical satellite P and the GEO spacecraft E, simulation based on CW equations

图18基于CW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.18Observation angles of the optical satellite P to the GEO spacecraft E based on CW equations simulation

图19基于ICW方程仿真的成像卫星P的控制推力

Fig.19Thrust control of the optical satellite P based on ICW equations simulation

图20基于ICW方程仿真的成像卫星P在x方向的推力

Fig.20Thrust control of the optical satellite P in x-coordinate simulation based on ICW equations

图21基于ICW方程仿真的成像卫星P在y方向的推力

Fig.21Thrust control of the optical satellite P in y-coordinate simulation based on ICW equations

图22基于ICW方程仿真的成像卫星P在z方向的推力

Fig.22Thrust control of the optical satellite P in z-coordinate simulation based on ICW equations

图23基于ICW方程仿真的成像卫星P相对 GEO航天器E的位置

Fig.23Position of the optical satellite P relative to the GEO spacecraft E based on ICW equations simulation

图24基于ICW方程仿真的成像卫星P与 GEO航天器E间的距离

Fig.24Distance between the optical satellite P and the GEO spacecraft E based on ICW equations simulation

图25基于ICW方程仿真的成像卫星P对 GEO航天器E的观测角

Fig.25Observation angles of the optical satellite P to the GEO spacecraft E based on ICW equations simulation

表1成像卫星P与GEO航天器E的轨道参数

表2成像卫星P的初始状态矢量

表3成像卫星P在LVLH坐标系中的状态矢量

ESA Space Debris Office. ESA′s space environment report[R/OL].(2021-05-27)[2021-06-27].https://www.esa.int/Space_Safety/Space_Debris/ESA_s_Space_Environment_Report_2021.

HEILIGERS J, MCINNES C R, BIGGS J D,et al. Displaced geostationary orbits using hybrid low-thrust propulsion[J]. Acta Astronautica,2012,71:51-67.

DE ALMEIDA A K, Jr, PIÑEROS J O M, DE ALME PRADO A F B. On the use of a continuous thrust to find bounded planar trajectories at given altitudes in low earth orbits[J]. Scientific Reports,2020,10:8728.

SULLIVAN J, GRIMBERG S, D′AMICO S. Comprehensive survey and assessment of spacecraft relative motion dynamics models[J]. Journal of Guidance, Control,and Dynamics,2017,40(8):1837-1859.

YIN J F, RAO Y R, HAN C. Inverse transformation of elliptical relative state transition matrix[J]. International Journal of Astronomy and Astrophysics,2014,4(3):419-428.

YIN J F, HAN C. Elliptical formation control based on relative orbit elements[J]. Chinese Journal of Aeronautics,2013,26(6):1554-1567.

MONTENBRUCK O, KIRSCHNER M, D′AMICO S,et al. E/I-vector separation for safe switching of the GRACE formation[J]. Aerospace Science and Technology,2006,10(7):628-635.

D′AMICO S, MONTENBRUCK O. Proximity operations of formation-flying spacecraft using an eccentricity/inclination vector separation[J]. Journal of Guidance, Control,and Dynamics,2006,29(3):554-563.

D′AMICO S. Autonomous formation flying in low earth orbit[D]. The Netherlands: Ridderprint BV,2010.

GAIAS G, ARDAENS J S, MONTENBRUCK O. Model of J2 perturbed satellite relative motion with time-varying differential drag[J]. Celestial Mechanics and Dynamical Astronomy,2015,123(4):411-433.

SPIRIDONOVA S. Formation dynamics in geostationary ring[J]. Celestial Mechanics and Dynamical Astronomy,2016,125(4):485-500.

VADDI S S, VADALI S R, ALFRIEND K T. Formation flying:accommodating nonlinearity and eccentricity perturbations[J]. Journal of Guidance, Control,and Dynamics,2003,26(2):214-223.

刘猛, 李全. 航天器对空间目标追踪及绕飞研究[J]. 沈阳理工大学学报,2011,30(2):80-83. LIU M, LI Q. Research on achievement of tracking and flying around space targe[J]. Journal of Shenyang Ligong University,2011,30(2):80-83.(in Chinese)

黄艺, 贾英民. 非合作目标绕飞任务的航天器鲁棒姿轨耦合控制[J]. 控制理论与应用,2018,35(10):1405-1414. HUANG Y, JIA Y M. Robust relative position and attitude control for non-cooperative fly-around mission[J]. Control Theory & Applications,2018,35(10):1405-1414.(in Chinese)

谭天乐. 椭圆轨道交会、悬停与绕飞的全状态反馈控制[J]. 宇航学报,2016,37(7):811-818. TAN T L. Full state feedback control of rendezvous,hovering and fly-around in elliptical orbit[J]. Journal of Astronautics,2016,37(7):811-818.(in Chinese)

李彬, 张洪波, 郑伟. 航天器近距离相对悬停被动安全特性分析[J]. 国防科技大学学报,2019,41(6):1-11. LI B, ZHANG H B, ZHENG W. Passive safety analysis of close relative hovering for spacecraft[J]. Journal of National University of Defense Technology,2019,41(6):1-11.(in Chinese)

孙永军, 王钤, 刘伊威, 等. 空间非合作目标捕获方法综述[J]. 国防科技大学学报,2020,42(3):74-90. SUN Y J, WANG Q, LIU Y W,et al. A survey of non-cooperative target capturing methods[J]. Journal of National University of Defense Technology,2020,42(3):74-90.(in Chinese)

刘涛, 王勇, 解永春, 等. 一种用于非合作目标惯性指向轴位置捕获的绕飞方法[J]. 宇航学报,2018,39(5):524-531. LIU T, WANG Y, XIE Y C,et al. A fly-around method for capturing position on inertia axis of a non-cooperative target[J]. Journal of Astronautics,2018,39(5):524-531.(in Chinese)

常燕, 陈韵, 鲜勇, 等. 椭圆轨道上目标监测绕飞轨道构型设计与构型保持[J]. 系统工程与电子技术,2017,39(6):1317-1324. CHANG Y, CHEN Y, XIAN Y,et al. Configuration design and maintenance of flyaround trajectory for target monitoring in elliptical orbit[J]. Systems Engineering and Electronics,2017,39(6):1317-1324.(in Chinese)

梁静静, 解永春. 基于长方体禁飞区的安全绕飞轨迹设计[J]. 空间控制技术与应用,2018,44(1):45-50. LIANG J J, XIE Y C. Design method of safe fly-around trajectory based on rectangular keep-out-zone[J]. Aerospace Control and Application,2018,44(1):45-50.(in Chinese)

王功波, 孟云鹤, 郑伟, 等. 快速绕飞卫星空间圆编队设计方法[J]. 宇航学报,2010,31(11):2465-2470. WANG G B, MENG Y H, ZHENG W,et al. Fast fly around satellite space circle formation design[J]. Journal of Astronautics,2010,31(11):2465-2470.(in Chinese)

DANG Z H. New state transition matrix for relative motion on an arbitrary Keplerian orbit[J]. Journal of Guidance, Control,and Dynamics,2017,40(11):2917-2927.

DANG Z H. Solutions of Tschauner-Hempel equations[J]. Journal of Guidance, Control,and Dynamics,2017,40(11):2953-2957.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

期刊订阅

Rss

AI检索

English

1 CW方程

2 改进的CW方程

3 轨迹规划

4 仿真

5 结论