CN113282097B - 一种geo博弈航天器相对位置非球形摄动误差的控制方法 - Google Patents

Abstract

本发明公开了一种GEO博弈航天器相对位置非球形摄动误差的控制方法，包括：以参考航天器的质心为原点建立LVLH坐标系；分别计算追踪航天器和逃逸航天器相对参考航天器的初始状态矢量；将追踪航天器和逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程，可得追踪航天器的状态矢量和逃逸航天器的状态矢量；逃逸航天器的状态矢量减去追踪航天器的状态矢量可得控制非球形摄动误差后的GEO博弈航天器相对位置。本发明针对地球静止轨道博弈航天器相对运动过程中，用CW方程得到的博弈航天器间相对位置存在非球形摄动误差的不足，添加对相对运动影响较大的摄动项来提高CW方程精度，同时又没有改变CW方程的线性特性，便于用变分方法研究航天器的博弈策略。

Description

一种GEO博弈航天器相对位置非球形摄动误差的控制方法

技术领域

本发明涉及一种地球静止轨道博弈航天器相对位置非球形摄动误差的控制方法，属于航天器轨道动力学领域。

背景技术

研究博弈相对运动的方法可以分为两类：基于运动学的相对运动，基于动力学的相对运动。基于运动学的相对运动一般的做法是：第一、建立主从航天器绝对轨道根数的差值与相对轨道根数的转换关系，一般可以得到精确的转换矩阵；第二、建立对描述编队问题比较有利的坐标系，研究相对轨道根数与该坐标系下状态量的转换矩阵；第三，由于根据高斯摄动方程可以得到三个方向的脉冲控制力与绝对轨道根数差值间的关系，进而建立控制力、绝对轨道根数差值、相对轨道根数、状态量间的转换关系，第四、对控制力施加的时机和大小进行优化，可用于分析编队的构建、重构最优控制问题。

其中比较有代表性的是殷建丰、韩潮等人基于相对轨道要素推导的椭圆相对运动方程。殷建丰、韩潮等人从运动学的方式提出了一套可用于近圆、椭圆编队设计、重构、优化控制的方法：在文献([1]韩潮,殷建丰.基于相对轨道要素的椭圆轨道卫星相对运动研究[J].航空学报,2011,32(12):2244-2258.Han Chao,Yin Jianfeng.Study of satelliterelative motion in elliptical orbit using relative orbit elements[j].ActaAeronautica et Astronautica sinica,2011,32(12):2244-2258.)中严格定义了相对轨道要素，并推导了其与绝对轨道要素间的转换矩阵，用于分析编队相对运动关系；在文献([1]Yin J,RaoY,Han C.Inverse Transformation of Elliptical Relative StateTransition Matrix[J].International Journal ofAstronomy&Astrophysics,2014,04(3):419-428.)中推导了相对轨道要素到质心轨道系状态量的转换矩阵和转换逆矩阵，建立了脉冲速度增量与轨道系中状态量的转化方程；在文献([1]Yin J,Han C.Ellipticalformation control based on relative orbit elements[J].中国航空学报(英文版),2013,26(006):1554-1567.)中，基于相对轨道根数，将轨道面内的控制与轨道法向的控制解耦，通过四次轨道机动，实现状态控制，并用数值方法对控制的时机以及控制量进行优化。殷建丰、韩潮的工作在博弈的单方的优化控制中具有借鉴意义，不足在于：基于相对轨道要素的运动学编队相对运动方程没法建立状态量对时间微分的状态方程，不便于航天器双方交互式的博弈优化策略的研究。

基于运动学研究相对运动的优点在于：模型的精度高，便于引入摄动对相对运动的影响，可以有效控制模型误差随时间的累积。缺点在于：无法通过微分方程对***进行描述，在微分博弈问题中，不便于研究双方交互式的博弈优化对策。

基于动力学方法对相对运动进行研究，Clohessy-Wiltshire方程，简称：CW方程，就是其中应用十分广泛的方法。基于动力学对相对运动进行研究的一般思路：第一、选取合适的状态量，确定绝对轨道根数的差值与状态量的转换关系；第二，建立状态量对时间的微分方程，即状态方程；第三、对微分方程进行分析，用于编队的构建、重构，优化控制。

所述CW方程为：

其中，X(t)表示t时刻的状态矢量，X(t₀)表示初始状态矢量，U为施加的控制力，B是常数矩阵，

s是积分自变量符号，U(s)是s时刻的控制力，φ(t,s)是s时刻到t时刻的状态转移矩阵，φ(t,t₀)是t₀时刻到t时刻的状态转移矩阵；

其中：τ＝ω(t-t₀)表示参考航天器从t₀时刻到t时刻转过的角度，

a_c为参考航天器的轨道半长轴，μ＝3.986005×10¹⁴m³/s²为地球引力常数。

不足在于：CW方程误差来源于线性化误差、摄动误差，CW方程还有小偏心率假设误差(比如cw方程)，误差较大，特别是不容易控制随时间累积的误差。优点在于：建立状态方程后，便于运用现代控制、最优控制的理论对相对运动开展研究，特别方便解决双方交互式的微分博弈问题。

针对GEO航天器的相对运动模型，比较有代表性的是STM模型。STM模型从最初建立，历经十多年的时间逐步完善。D’Amico([1]Montenbruck,O.,Kirschner,M.,D’Amico,S.,and Bettadpur,S.,“E/I-Vector Separation for Safe Switching ofthe GRACEFormation,”Aerospace Science and Technology,Vol.10,No.7,2006,pp.628–635.doi:10.1016/j.ast.2006.04.001；[2]D’Amico,S.,and Montenbruck,O.,“ProximityOperations ofFormation-Flying Spacecraft Using an Eccentricity/InclinationVector Separation,”Journal of Guidance,Control,and Dynamics,Vol.29,No.3,2006,pp.554–563.doi:10.2514/1.15114；[3]D’Amico,S.,“Autonomous Formation Flying inLow Earth Orbit,”Ph.D.)研究了基于相对E/I矢量的相对轨道根数编队动力学模型STM模型，在STM模型中考虑了J₂项摄动，并研究了编队安全路径问题；Gaias,G(Gaias,G.,Ardaens,J.-S.,and Montenbruck,O.,“Model ofJ2 Perturbed Satellite RelativeMotion with Time-Varying Differential Drag,”Celestial Mechanics and DynamicalAstronomy,Vol.123,No.4,2015,pp.411–433.doi:10.1007/s10569-015-9643-2)引入了大气阻力摄动，分析了大气阻力摄动对轨道半长轴的长期影响，改进了STM模型；Spiridonova,S(Spiridonova,S.,“Formation Dynamics in Geostationary Ring,”Celestial Mechanics and Dynamical Astronomy,Vol.125,No.4,2016,pp.485–500.doi:10.1007/s10569-016-9693-0)在STM模型中加入了太阳光压和三体摄动，消除了相对轨道根数漂移的所有长期和长周期项，完善了STM模型，仿真表明模型在10天之后也具有较高的精度，使STM模型适用于包含近距离空间操作等在轨服务任务的GEO编队问题。虽然具有较高的精度，但Spiridonova,S完善后的STM模型复杂，包含非线性部分，无法建立线性的状态控制方程，不便于用于GEO航天器博弈问题。

发明内容

本发明提供了一种GEO博弈航天器相对位置非球形摄动误差的控制方法，地球静止轨道记为GEO，是一种地球静止轨道博弈航天器相对位置非球形摄动误差的控制方法，本发明针对CW方程中存在的摄动误差,造成航天器相对位置算的不准，从而导致航天任务无法正常开展的技术问题，通过修正提高模型的精度，在研究GEO航天器博弈问题时，使用该方法计算航天器间的相对位置和速度，控制了地球静止轨道博弈航天器摄动误差的传播，提高了相对位置的精度。

为了达到上述目的，本发明提供了一种GEO博弈航天器相对位置非球形摄动误差的控制方法，该方法包括：

步骤一，选择参考航天器，参考航天器的轨道半长轴为GEO标称轨道半长轴，偏心率为0，轨道倾角为0，升交点赤经、近地点俯角和平近地点角与追踪航天器或逃逸航天器相同，以参考航天器的质心为原点建立LVLH坐标系；

步骤二，在LVLH坐标系中，分别计算追踪航天器和逃逸航天器相对参考航天器的初始状态矢量；

步骤三，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程，可得追踪航天器的状态矢量；将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程，可得逃逸航天器的状态矢量；逃逸航天器的状态矢量减去追踪航天器的状态矢量可得控制非球形摄动误差后的GEO博弈航天器相对位置；

所述修正CW方程为：

其中，X(t)表示t时刻的状态矢量，X(t₀)表示初始状态矢量，U为施加的控制力，B是常数矩阵，

ε是积分自变量符号，U(ε)是ε时刻的控制力，φ(t,ε)是ε时刻到t时刻的状态转移矩阵，φ(t,t₀)是t₀时刻到t时刻的状态转移矩阵；

其中：τ＝ω(t-t₀)表示参考航天器从t₀时刻到t时刻转过的角度，

a_c为参考航天器的轨道半长轴；

式中，r_E为地球半径，J₂＝-1082.627×10^-6、J₂₂＝1.815528×10^-6，a_c为参考航天器的轨道半长轴，λ为参考航天器的星下点地理经度，μ＝3.986005×10¹⁴m³/s²为地球引力常数。

进一步的，所述LVLH坐标系，为当地水平当地垂直坐标系，是指以参考航天器的质心为原点，以地心到参考航天器的矢径为x轴的正向，z轴指向参考航天器轨道的正法向，y轴的方向根据右手定则确定。

所述步骤二中，在LVLH坐标系中，计算追踪航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_P(t₀),y_P(t₀),z_P(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在初始时刻的值；T表示转置；

在LVLH坐标系中，计算逃逸航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_E(t₀),y_E(t₀),z_E(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标，

为x_E,y_E,z_E对时间导函数在初始时刻的值；T表示转置。

所述步骤三中，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程

中，可得t时刻追踪航天器的状态矢量

其中，(x_P(t),y_P(t),z_P(t))为t时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在t时刻的值，T表示转置；

将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程

可得t时刻逃逸航天器的状态矢量

其中，(x_E(t),y_E(t),z_E(t))为t时刻逃逸航天器在LVLH坐标系中的位置坐标；

为x_E,y_E,z_E对时间导函数在t时刻的值，T表示转置；

(x_E(t)-x_P(t),y_E(t)-y_P(t),z_E(t)-z_P(t))为控制非球形摄动误差后的GEO博弈航天器的相对位置。

其中，所述追踪航天器可用于对空间目标实施抓捕任务；

逃逸航天器躲避追踪航天器的抓捕；

所述参考航天器为虚拟的不做轨道机动的航天器。

本发明有益效果如下：

相比CW方程计算的地球静止轨道博弈航天器的相对位置，控制了摄动误差的传播，体现在：

a.针对传统动力学相对运动CW方程中存在的摄动误差问题，考虑GEO航天器的轨道共振特性，引入对相对运动影响最大的J₂和J₂₂，避免了非球形摄动的影响；

b.针对传统动力学相对运动CW方程中存在的偏心率误差问题，步骤一在选取参考航天器时，令参考航天器的偏心率为0，避免偏心率误差的传播。

附图说明：

图1是μ′/μ随GEO航天器定点经度的变化示意图。

图2是本方法与CW方程中x状态量的差值随时间的变化示意图。

图3是本方法与CW方程中y状态量的差值随时间的变化示意图。

图4是本方法与CW方程中z状态量的差值随时间的变化示意图。

图5是本方法与CW方程中x状态量的差值随时间的变化示意图。

图6是本方法与CW方程中y状态量的差值随时间的变化示意图。

图7是本方法与CW方程中z状态量的差值随时间的变化示意图。

图8是本方法与CW方程计算的x状态量10T后差值的幅值示意图。

图9是本方法与CW方程计算的y状态量10T后差值的幅值示意图。

图10是本方法与CW方程计算的z状态量10T后差值的幅值示意图。

具体实施方式

一种地球静止轨道博弈航天器相对位置摄动误差的控制方法，针对CW方程中存在的摄动误差问题，通过修正提高CW方程的精度，在研究GEO航天器博弈问题时，使用该方法计算航天器间的相对位置和速度，控制了地球静止轨道博弈航天器摄动误差的传播，提高了相对位置的精度。

为了达到上述目的，本发明提供了一种GEO博弈航天器相对位置非球形摄动误差的控制方法，该方法包括：

步骤一，选择参考航天器，参考航天器的轨道半长轴为GEO标称轨道半长轴，偏心率为0，轨道倾角为0，升交点赤经、近地点俯角和平近地点角与追踪航天器或逃逸航天器相同，以参考航天器的质心为原点建立LVLH坐标系；

步骤二，在LVLH坐标系中，分别计算追踪航天器和逃逸航天器相对参考航天器的初始状态矢量；

步骤三，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程，可得追踪航天器的状态矢量；将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程，可得逃逸航天器的状态矢量；逃逸航天器的状态矢量减去追踪航天器的状态矢量可得控制非球形摄动误差后的GEO博弈航天器相对位置；

所述修正CW方程为：

其中，X(t)表示t时刻的状态矢量，X(t₀)表示初始状态矢量，U为施加的控制力，B是常数矩阵，

ε是积分自变量符号，U(ε)是ε时刻的控制力，φ(t,ε)是ε时刻到t时刻的状态转移矩阵，φ(t,t₀)是t₀时刻到t时刻的状态转移矩阵；

其中：τ＝ω(t-t₀)表示参考航天器从t₀时刻到t时刻转过的角度，

a_c为参考航天器的轨道半长轴；

式中，r_E为地球半径，J₂＝-1082.627×10^-6、J₂₂＝1.815528×10^-6，a_c为参考航天器的轨道半长轴，λ为参考航天器的星下点地理经度，μ＝3.986005×10¹⁴m³/s²为地球引力常数。

进一步的，所述LVLH坐标系，为当地水平当地垂直坐标系，是指以参考航天器的质心为原点，以地心到参考航天器的矢径为x轴的正向，z轴指向参考航天器轨道的正法向，y轴的方向根据右手定则确定。

具体的，地球静止轨道记为GEO。当地垂直当地水平坐标系(Local VerticalLocal Horizontal,LVLH),该坐标系以航天器的质心为坐标原点，以地心到航天器的质心的方向e_x为x轴的方向，以航天器轨道面的正法向e_z为z轴的方向，y轴的方向e_y通过下式确定：

e_y＝e_z×e_x

记P为追踪航天器，E为逃逸航天器。选择距离i∈{P,E}都不太远的虚拟的不做轨道机动的航天器作为参考航天器，参考航天器记为c，建立以c的质心为原点的LVLH坐标系，记为：

所述步骤二中，在LVLH坐标系中，计算追踪航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_P(t₀),y_P(t₀),z_P(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在初始时刻的值；T表示转置；

在LVLH坐标系中，计算逃逸航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_E(t₀),y_E(t₀),z_E(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标，

为x_E,y_E,z_E对时间导函数在初始时刻的值；T表示转置。

具体的，记r_LVLH为追踪航天器P在LVLH坐标系中的位置矢量，r_ECI为追踪航天器P在地心惯性坐标系ECI(Earth centered-inertial frame)中的位置矢量，A_ECI→LVLH为ECI坐标系到LVLH坐标系的坐标变换矩阵，A_LVLH→ECI为LVLH坐标系到ECI坐标系的坐标变换矩阵，故：

r_ECI＝A_LVLH→ECI·r_LVLH (1)

r_LVLH＝A_ECI→LVLH·r_ECI (2)

其中：

对于列向量a＝[a₁ a₂ a₃]^T和b＝[b₁ b₂ b₃]^T，存在

其中：

A_ECI→LVLH和A_LVLH→ECI对时间的变化率为：

式中n为LVLH坐标系相对ECI坐标系转动的角速度，且：

式中

为追踪航天器P与参考航天器c的位置速度矢量在ECI坐标系中的差值，

即为追踪航天器P相对参考航天器c的状态量，初始时刻的状态量即为所需的初始状态量

相同的方法可得逃逸航天器E相对参考航天器c的状态量，初始时刻的状态量即为所需的初始状态量

所述步骤三中，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程

中，可得t时刻追踪航天器的状态矢量

其中，(x_P(t),y_P(t),z_P(t))为t时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在t时刻的值，T表示转置；

将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程

可得t时刻逃逸航天器的状态矢量

其中，(x_E(t),y_E(t),z_E(t))为t时刻逃逸航天器在LVLH坐标系中的位置坐标；

为x_E,y_E,z_E对时间导函数在t时刻的值，T表示转置；

(x_E(t)-x_P(t),y_E(t)-y_P(t),z_E(t)-z_P(t))为控制非球形摄动误差后的GEO博弈航天器的相对位置。

其中，所述追踪航天器可用于对空间目标实施抓捕任务；

逃逸航天器躲避追踪航天器的抓捕；

所述参考航天器为虚拟的不做轨道机动的航天器。

具体的，航天器近距离相对运动，特别是对于博弈的两个航天器，相对距离较近，当：

式中：ρ表示两个航天器间的相对距离，a表示其中一个航天器的轨道半长轴。

以参考航天器的质心为原点的LVLH坐标系

中，记从航天器的相对状态量为

推力加速度U＝(u_x,u_y,u_z)^T，则：

式中：

t₀时刻，从航天器的相对状态量为：

则，t时刻，从航天器的状态为：

式中：

φ_rr(t,t₀)、φ_rv(t,t₀)、φ_vr(t,t₀)和φ_vv(t,t₀)都是3×3的矩阵，且：

其中：τ＝ω(t-t₀)，s＝sinτ，c＝cosτ。

以航天器在距离地球无穷远处的势为0，地球引力场的势函数为：

其中，μ＝3.986005×10¹⁴m³/s²为地球引力常数，r为航天器距离地心的距离，r_E为地球半径，J₂＝-1082.627×10^-6、J₂₂＝1.815528×10^-6、J₃＝2.532435×10^-6、J₃₁＝2.2091169×10^-6、J₃₃＝0.2213602×10^-6、

2λ₂₂＝-29.8568°、λ₃₁＝6.9683°、3λ₃₃＝62.9815°，将

归一化后，记为：

归一化系数为：

整理后，地球引力场的势函数为：

各项的权重为：

归一化后的谐系数，除J₂的量级较大外，其余的量级一般都较小。针对地球静止轨道附近的航天器，计算得权重最大的五项为：J₂、J₂₂、J₃、J₃₁、J₃₃，其权重比值为：

A₂:A₂₂:A₃:A₃₁:A₃₃＝11053:64:3:7:6

因此，在本方法中引入地球非球形摄动中权重最大的两项J₂和J₂₂项。考虑J₂和J₂₂项摄动的地球引力场势函数：

式中，

为航天器的星下点地理纬度，由于参与博弈的GEO航天器相对距离较近，且博弈过程中星下点地理经度变化较小，即：

则：

式中，λ为航天器的星下点地理经度，a_s为地球静止轨道航天器标称半长轴。

J₂和J₂₂项对CW方程的影响体现在：本方法中修正的地球引力常数μ′不再是步骤二中的地球引力常数μ。其中：

相比CW方程，引入J₂和J₂₂项摄动后，GEO航天器地球引力势函数的绝对值变大，且随定点经度的变化，μ′与μ的比值如图1所示。J₂₂项表明地球静止轨道存在引力势的两个极小值点：75.0716°E和104.9284°W。J₂₂项导致地球静止轨道发生共振特性，在不加控制的情形下，定点在14.9284°W～165.0716°E的地球静止轨道航天器在J₂₂的作用下，受到向75.0716°E漂移的摄动；定点在165.0716°E～14.9284°W的地球静止轨道航天器在J₂₂的作用下，受到向104.9284°W漂移的摄动。

一优选实例

为了分析不同初始条件下，本方法的误差，设置1号航天器和2号航天器的轨道根数如表1所示。

表1

取参考航天器c的轨道根数如表2所示。

表2

1号航天器和2号航天器在

中的初始状态量如表3所示。

表3

分别运用本方法和CW方程仿真表1中的1号和2号航天器与参考航天器的相对运动关系，对比本方法与CW方程的x,y,z状态量的差值随时间的变化情况。仿真案例中，航天器的星下点地理经度取165.0716°E。

a.案例1：ρ＝200km,α₀＝180°

本仿真案例下，图2，图3，图4分别表示本方法与CW方程中x,y,z状态量的差值随时间的变化。从图中可以看出：本方法与CW方程中x,y,z状态量的差值呈现周期性变化，差值的幅值逐渐扩大，x状态量差值变化的周期与案例2相差90°相位，z状态量差值变化的周期与前述案例相同。

b.案例2：ρ＝200km,α₀＝270°

本仿真案例下，图5，图6，图7分别表示本方法与CW方程中x,y,z状态量的差值随时间的变化。从图中可以看出：本方法与CW方程中x,z状态量的差值呈现周期性变化，差值的幅值逐渐扩大，x状态量差值变化的周期与案例3相差90°相位，z状态量差值变化的周期与前述案例相同；y状态量的差值包含随时间的累积项和短周期项，差值幅值逐渐扩大，10个轨道周期后相差约300m。

c.不同定点经度对差值增长幅值的影响

从案例1和2中，α₀＝180°,270°两种初始状态下，本方法与CW方程中x状态量的差值随时间都呈周期性变化，且变化的幅值逐渐增大；z状态量的差值随时间也都呈周期性变化，且变化的幅值逐渐增大，但幅值增大的速率比x状态量差值幅值增大的速率小一个量级。

α₀＝180°两种初始状态下，本方法与CW方程中y状态量的差值随时间呈周期性变化，且变化的幅值逐渐增大，幅值增大的速率与x状态量差值的幅值增大的速率相当。α₀＝270°两种初始状态下，y状态量的差值包含随时间的累积项和短周期项，差值幅值逐渐扩大。取ρ＝200km,α₀＝270°，分析不同定点经度处的博弈航天器用本方法与CW方程计算的x,y,z状态量10T后差值的幅值如图8，图9，图10所示。

从图中不同定点经度处的博弈航天器用本方法与CW方程计算的x,y,z状态量10T后差值的幅值与J₂₂项摄动势函数呈现相同的变化规律。非球形摄动会造成x状态量出现约40.5m的偏差；非球形摄动会造成y状态量出现约311.5m的偏差；非球形摄动会造成z状态量出现约2.87m的偏差。非球形摄动对x,y状态量的影响大于对z状态量的影响。

本发明针对地球静止轨道博弈航天器相对运动过程中，用CW方程得到的博弈航天器间相对位置存在非球形摄动误差的不足，添加对相对运动影响较大的摄动项来提高CW方程精度，同时又没有改变CW方程的线性特性，便于用变分方法研究航天器的博弈策略。

本发明实施例有益效果如下：

相比CW方程，避免了CW方程误差的传播，体现在：

a.针对传统动力学相对运动CW方程中存在的摄动误差问题，考虑GEO航天器的轨道共振特性，本方法中加入了对相对运动影响最大的J₂和J₂₂避免了非球形摄动的影响；

b.针对传统动力学相对运动CW方程中存在的偏心率误差问题，在选取c时，令c的偏心率为0，避免偏心率误差的传播。

Claims (5)

1.一种GEO博弈航天器相对位置非球形摄动误差的控制方法，其特征在于，该方法包括：

步骤一，选择参考航天器，参考航天器的轨道半长轴为GEO标称轨道半长轴，偏心率为0，轨道倾角为0，升交点赤经、近地点俯角和平近地点角与追踪航天器或逃逸航天器相同，以参考航天器的质心为原点建立LVLH坐标系；

步骤二，在LVLH坐标系中，分别计算追踪航天器和逃逸航天器相对参考航天器的初始状态矢量；

步骤三，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程，可得追踪航天器的状态矢量；将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程，可得逃逸航天器的状态矢量；逃逸航天器的状态矢量减去追踪航天器的状态矢量可得控制非球形摄动误差后的GEO博弈航天器相对位置；

所述修正CW方程为：

其中，X(t)表示t时刻的状态矢量，X(t₀)表示初始状态矢量，U为施加的控制力，B是常数矩阵，

ε是积分自变量符号，U(ε)是ε时刻的控制力，φ(t,ε)是ε时刻到t时刻的状态转移矩阵，φ(t,t₀)是t₀时刻到t时刻的状态转移矩阵；

其中：τ＝ω(t-t₀)表示参考航天器从t₀时刻到t时刻转过的角度，

a_c为参考航天器的轨道半长轴；

式中，r_E为地球半径，J₂＝-1082.627×10^-6、J₂₂＝1.815528×10^-6，λ为参考航天器的星下点地理经度，μ＝3.986005×10¹⁴m³/s²为地球引力常数。

2.如权利要求1所述的方法，其特征在于，所述LVLH坐标系，为当地水平当地垂直坐标系，是指以参考航天器的质心为原点，以地心到参考航天器的矢径为x轴的正向，z轴指向参考航天器轨道的正法向，y轴的方向根据右手定则确定。

3.如权利要求1或2所述的方法，其特征在于，所述步骤二中，在LVLH坐标系中，计算追踪航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_P(t₀),y_P(t₀),z_P(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在初始时刻的值；T表示转置；

在LVLH坐标系中，计算逃逸航天器相对参考航天器的初始状态矢量的方法包括：

其中，(x_E(t₀),y_E(t₀),z_E(t₀))为初始时刻追踪航天器在LVLH坐标系中的位置坐标，

为x_E,y_E,z_E对时间导函数在初始时刻的值；T表示转置。

4.如权利要求1所述的方法，其特征在于，所述步骤三中，将追踪航天器相对参考航天器的初始状态矢量代入修正CW方程

中，可得t时刻追踪航天器的状态矢量

其中，(x_P(t),y_P(t),z_P(t))为t时刻追踪航天器在LVLH坐标系中的位置坐标；

为x_P,y_P,z_P对时间导函数在t时刻的值，T表示转置；

将逃逸航天器相对参考航天器的初始状态矢量代入修正CW方程

可得t时刻逃逸航天器的状态矢量

其中，(x_E(t),y_E(t),z_E(t))为t时刻逃逸航天器在LVLH坐标系中的位置坐标；

为x_E,y_E,z_E对时间导函数在t时刻的值，T表示转置；

(x_E(t)-x_P(t),y_E(t)-y_P(t),z_E(t)-z_P(t))为控制非球形摄动误差后的GEO博弈航天器的相对位置。

5.如权利要求1所述的方法，其特征在于，所述追踪航天器可用于对空间目标实施抓捕任务；

逃逸航天器躲避追踪航天器的抓捕；

所述参考航天器为虚拟的不做轨道机动的航天器。

Images