CN114394260B - 地球同步卫星小推力南北保持控制解析计算方法 - Google Patents

Abstract

本发明公开的地球同步卫星小推力南北保持控制解析计算方法，具体计算过程包括如下步骤：步骤1、解析计算地球同步卫星倾角矢量的长期摄动变化量；步骤2、计算电推力器点火序列，解析计算出整个规划周期内每一次电推点火的开关机时刻。本发明实现地球同步卫星小推力南北保持点火序列的计算，通过解析方程精确计算出同步卫星倾角矢量的长期项摄动变化量，点火序列计算简便精准，且点火时长接近理论最优，具有计算速度快、适用性强的优点。

Description

地球同步卫星小推力南北保持控制解析计算方法

技术领域

本发明属于航天导航技术领域，具体涉及地球同步卫星小推力南北保持控制解析计算方法。

背景技术

通常来说，为满足星上的通信、光学等载荷的地面指向性或者覆盖性需求，地球同步卫星的倾角矢量要始终处于指定的较小范围内。在地球非球形引力和日月引力摄动的影响下，卫星的倾角矢量会逐渐漂移，为了克服这种自然摄动产生的倾角矢量漂移，而定期进行的倾角矢量控制，就是工程上俗称的“南北保持”。

在长期管理阶段，地球同步卫星的南北保持是最耗燃料的卫星控制工作，仅抵消一年的轨道倾角长期项摄动就需要给卫星提供约43m/s的速度增量，消耗数十千克燃料质量。随着小推力技术逐渐成熟，越来越多的地球静止卫星采用小推力推力器进行南北控制，小推力推力器最重要的优点是比冲大，提供相同速度增量，仅需要消耗化学推进器10％左右的燃料质量。然而，目前的小推力推力器推力一般只有几十到几百毫牛，实现南北保持需求的倾角控制量需要长时间开机(每天数小时)。此外，小推力推力器工质的动能来源是星上电能(化学推力器是燃料化学能)，过长的开机时间对卫星能源供应是个挑战。因此，根据卫星应用约束合理设计电推力器南北控制策略，使卫星以尽可能少的开机时间，对于完成南北保持控制是很有意义的。

传统的化学推进南北控制是在倾角即将超出限制区域时，选取合适的时刻利用法，向速度增量将倾角矢量向年平均摄动方向的反向控制，只在倾角矢量i_x方向到达边界限制时，对控制方向做一些调整。此策略的优点是，基本控制方向一年才做一次调整，便于操作，且控制方向也近似是燃料最优控制方向。年平均摄动方向的计算通常有两种方法——解析公式或轨道外推，解析公式计算会不可避免的存在一定误差，而轨道外推计算不能提出长周期项的影响。因此，如果在小推力条件下长期固定倾角矢量控制方向，方向误差造成的累积效应会使倾角矢量与控制目标显著偏差。

发明内容

本发明的目的是提供地球同步卫星小推力南北保持控制解析计算方法，解决了现有装配了电推力器的地球同步卫星南北保持任务，解析公式计算存在误差的问题。

本发明所采用的技术方案是，地球同步卫星小推力南北保持控制解析计算方法，具体计算过程包括如下步骤：

步骤1、解析计算地球同步卫星倾角矢量的长期摄动变化量；

步骤2、计算电推力器点火序列，解析计算出整个规划周期内每一次电推点火的开关机时刻。

本发明的其他特点还在于，

步骤1的具体步骤为：确定倾角矢量为(i_xi_y)，其中i_x＝icosΩ，i_y＝isinΩ，i和Ω分别为地球同步卫星轨道要素的倾角和升交点赤经；计算规划时段t₀到t_f内，倾角矢量的长期摄动变化量(Δi_x)_长期和(Δi_y)_长期。

步骤1具体按照以下步骤：

步骤1.1、计算地球非球形引力对倾角矢量的摄动量；

步骤1.2、计算太阳引力对倾角矢量的摄动量；

步骤1.3、计算月球引力对倾角矢量的摄动量；

步骤1.4、将步骤1.1至步骤1.3分别求得的地球非球形引力、太阳引力、月球引力对倾角矢量的摄动量求和，即为地球同步卫星倾角矢量的长期项摄动变化量。

步骤1.1中地球非球形引力对倾角矢量的摄动量计算公式为：

其中，n为卫星轨道角速度，a为卫星轨道半长轴，μ_e为地球引力常数，R_e为地球半径，J₂、J₂₂、J₃₁、J₃₃、J₄、J₄₃、λ₂₂、λ₃₁、λ₃₃、λ₄₂、λ₄₃是地球非球形四阶主项摄动系数，i_x、i_y为倾角矢量初值。

步骤1.2中太阳引力对倾角矢量的摄动量计算公式为：

其中，μ_s为太阳引力常数，r_s为日地距离，Ω_s为地心惯性系下的太阳轨道升交点赤经，i_s为太阳轨道倾角。

步骤1.3中月球引力对倾角矢量的摄动量计算公式为：

其中，μ_l为月球引力常数，r_l为地月距离，Ω_l为月球轨道升交点赤经，i_l为月球轨道倾角。

步骤1.4中地球同步卫星倾角矢量的长期项摄动变化量由上述三个部分组成，即

(Δi_x)_长期＝(Δi_x)_{非球形，长期}+(Δi_x)_{太阳，长期}+(Δi_x)_{月球，长期} (7)

(Δi_y)_长期＝(Δi_y)_{非球形，长期}+(Δi_y)_{太阳，长期}+(Δi_y)_{月球，长期} (8)。

步骤2的具体过程如下：

步骤2.1、根据步骤1中地球同步卫星倾角矢量的长期项摄动变化量得到一个周期南北保持的倾角矢量目标控制量；一个周期南北保持的倾角矢量目标控制量能够抵消周期内的长期项摄动量，即

(Δi_y)_control＝-(Δi_y)_长期 (9)

(Δi_x)_control＝-(Δi_x)_长期 (10)；

步骤2.2、计算每次小推力点火时长，得到规划周期内的倾角控制量；

根据步骤2.1中一个周期南北保持的倾角矢量目标控制量得到控制方向角，控制方向角α_i＝arctan2[(Δi_y)_control,(Δi_x)_control]，南侧电推力器的点火中心赤经为InnerTwoPi(α_i)，北侧电推力器的点火中心赤经为InnerTwoPi(α_i+π)；

点火序列规划周期为N天，其中的M天进行电推点火(M≤N，为轨道确定、载荷控制预留空间)，每天有2次南北控制的时机，电推点火总次数为2M，计算每次小推力点火时长为

T_fire＝2arcsin[(Δi)_controlV_sω_e/acc/M/4]/ω_e (11)

其中，(Δi)_control为规划周期内的倾角控制量，acc为电推力器法向合加速度，V_s为地球同步卫星标称轨道速度，ω_e为地球自转角速度；

电推力器一天内的两次点火时刻t_1,on、t_1,off和t_2,on、t_2,off分别由南侧和北侧点火中心赤经和点火时长计算出来，设当天00:00:00.00时刻的格林尼治恒星时角为Θ₀，卫星平经度为λ；

t_1,on＝InnerTwoPi(α_i-λ-Θ₀)/ω_e-T_fire/2 (13)

t_1,off＝t_1,on+T_fire (14)

t_2,on＝InnerTwoPi(α_i+π-λ-Θ₀)/ω_e-T_fire/2 (15)

t_2,off＝t_2,on+T_fire (16)

计算出来的t_1,on、t_1,off和t_2,on、t_2,off是从当天00:00:00.00时刻开始的积秒，以此类推，解析计算出整个规划周期内每一次电推点火的开关机时刻。

本发明的有益效果是，本发明地球同步卫星小推力南北保持控制计算方法，适应于小推力条件下的地球同步卫星南北位置保持控制机动规划，针对装配了电推力器的地球同步卫星南北保持任务，提出一种解析计算电推力器点火序列的方法，控制卫星抵消长期项倾角矢量摄动，使卫星倾角矢量在小范围内做周期性运动，从而实现地球同步卫星小推力南北保持点火序列的计算，通过解析方程精确计算出同步卫星倾角矢量的长期项摄动变化量，点火序列计算简便精准，且点火时长接近理论最优，提高了计算速度，满足快速规划计算的需求，具备广泛的适用性。

附图说明

图1是本发明一种实施例连续进行一年小推力南北控制的卫星倾角矢量变化情况图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明的地球同步卫星小推力南北保持控制解析计算方法，具体包括如下步骤：

步骤1、解析计算地球同步卫星倾角矢量的长期摄动变化量；

步骤1的具体步骤为：确定倾角矢量为(i_xi_y)，其中i_x＝icosΩ，i_y＝isinΩ，i和Ω分别为地球同步卫星轨道要素的倾角和升交点赤经；计算规划时段t₀到t_f内，倾角矢量的长期摄动变化量(Δi_x)_长期和(Δi_y)_长期；

步骤1.1、计算地球非球形引力对倾角矢量的摄动量；

地球非球形引力对倾角矢量的摄动量计算公式为：

其中，n为卫星轨道角速度，a为卫星轨道半长轴，μ_e为地球引力常数，R_e为地球半径，J₂、J₂₂、J₃₁、J₃₃、J₄、J₄₃、λ₂₂、λ₃₁、λ₃₃、λ₄₂、λ₄₃是地球非球形四阶主项摄动系数，i_x、i_y为倾角矢量初值。

步骤1.2、计算太阳引力对倾角矢量的摄动量；

太阳引力对倾角矢量的摄动量计算公式为：

其中，μ_s为太阳引力常数，r_s为日地距离，Ω_s为地心惯性系下的太阳轨道升交点赤经，i_s为太阳轨道倾角。

步骤1.3、计算月球引力对倾角矢量的摄动量；

月球引力对倾角矢量的摄动量计算公式为：

其中，μ_l为月球引力常数，r_l为地月距离，Ω_l为月球轨道升交点赤经，i_l为月球轨道倾角。

步骤1.4、将步骤1.1至步骤1.3分别求得的地球非球形引力、太阳引力、月球引力对倾角矢量的摄动量求和，即为地球同步卫星倾角矢量的长期项摄动变化量。

地球同步卫星倾角矢量的长期项摄动变化量由上述三个部分组成，即

(Δi_x)_长期＝(Δi_x)_{非球形，长期}+(Δi_x)_{太阳，长期}+(Δi_x)_{月球，长期} (7)

(Δi_y)_长期＝(Δi_y)_{非球形，长期}+(Δi_y)_{太阳，长期}+(Δi_y)_{月球，长期} (8)

步骤2、计算电推力器点火序列，解析计算出整个规划周期内每一次电推点火的开关机时刻。

电推力器能够提供的推力只有几十/几百毫牛，完成南北保持任务需要平均每天开机数小时，所以工程上通常采取一次性计算包含多天多次控制的电推力器点火序列的方式进行。

步骤2.1、根据步骤1中地球同步卫星倾角矢量的长期项摄动变化量得到一个周期南北保持的倾角矢量目标控制量，一个周期南北保持的倾角矢量目标控制量能够抵消周期内的长期项摄动量，即

(Δi_y)_control＝-(Δi_y)_长期 (9)

(Δi_x)_control＝-(Δi_x)_长期 (10)步骤2.2、计算每次小推力点火时长，得到规划周期内的倾角控制量；

根据步骤2.1中一个周期南北保持的倾角矢量目标控制量得到控制方向角，控制方向角α_i＝arctan2[(Δi_y)_control,(Δi_x)_control]，南侧电推力器的点火中心赤经为InnerTwoPi(α_i)，北侧电推力器的点火中心赤经为InnerTwoPi(α_i+π)。

点火序列规划周期为N天，其中的M天进行电推点火(M≤N，为轨道确定、载荷控制预留空间)，每天有2次南北控制的时机，电推点火总次数为2M，计算每次小推力点火时长为

T_fire＝2arcsin[(Δi)_controlV_sω_e/acc/M/4]/ω_e (11)

其中，(Δi)_control为规划周期内的倾角控制量，acc为电推力器法向合加速度，V_s为地球同步卫星标称轨道速度，ω_e为地球自转角速度。

电推力器一天内的两次点火时刻t_1,on、t_1,off和t_2,on、t_2,off分别由南侧和北侧点火中心赤经和点火时长计算出来。设当天00:00:00.00时刻的格林尼治恒星时角为Θ₀，卫星平经度为λ。

t_1,on＝InnerTwoPi(α_i-λ-Θ₀)/ω_e-T_fire/2 (13)

t_1,off＝t_1,on+T_fire (14)

t_2,on＝InnerTwoPi(α_i+π-λ-Θ₀)/ω_e-T_fire/2 (15)

t_2,off＝t_2,on+T_fire (16)

计算出来的t_1,on、t_1,off和t_2,on、t_2,off是从当天00:00:00.00时刻开始的积秒，以此类推，解析计算出整个规划周期内每一次电推点火的开关机时刻。

通过上述方法，即可实现地球同步卫星小推力南北保持点火序列的计算，其特点是通过解析方程精确计算出同步卫星倾角矢量的长期项摄动变化量，点火序列计算简便精准，且点火时长接近理论最优。

实施例

本实施例以亚太6D卫星一个周期的小推力南北控制序列计算为例，如图1所示，南北控制周期为7天(2021年6月7日至2021年6月13日)，其中前5天进行电推点火。

利用式(1)至式(8)计算7天的倾角矢量长期摄动变化量。

带入常数λ₂₂＝-14.929°、λ₃₁＝6.968°、λ₃₃＝20.994°、λ₄₄＝30.280°、J₂₂＝1.811528e^-6、J₃₁＝2.2091169e^-6、J₃₃＝0.2213602e^-6、J₄₄＝7.63937899e^-9、R_e＝6378.1363、μ＝398600.4415和a＝42165.700、n＝7.292116e^-5、i_x＝0.02°、i_y＝0.03°、t_f-t₀＝86400.0×7。

带入常数μ_s＝1.3271222e¹¹、r_s＝1.496e⁸、Ω_s＝0.0°、i_s＝23.4437°；

带入常数μ_l＝4.90280107e³、r_l＝3.85e⁵(月球倾角i_l和升交点赤经Ω_l是时变量，由当前日期计算)；

(Δi_x)_长期＝(Δi_x)_{非球形，长期}+(Δi_x)_{太阳，长期}+(Δi_x)_{月球，长期} (7)

(Δi_y)_长期＝(Δi_y)_{非球形，长期}+(Δi_y)_{太阳，长期}+(Δi_y)_{月球，长期} (8)代数式(1)至式(6)计算结果

(Δi_x)_长期＝-0.0023°

(Δi_y)_长期＝0.0169°

由式(9)和(10)计算倾角矢量目标控制量

(Δi_y)_control＝-(Δi_y)_长期 (9)

(Δi_x)_control＝-(Δi_x)_长期 (10)

(Δi_y)_control＝0.0023°

(Δi_x)_control＝-0.0169°

控制方向角α_i＝arctan2[(Δi_y)_control,(Δi_x)_control]＝277.8198°，南侧电推力器的点火中心赤经为InnerTwoPi(α_i)，北侧电推力器的点火中心赤经为

InnerTwoPi(α_i+π)。

N＝7,M＝5

带入卫星电推加速度acc，计算每次小推力点火时长为

T_fire＝2arcsin[(Δi)_controlV_sω_e/acc/M/4]/ω_e (11)

T_fire＝10785(s)

规划周期内的倾角控制量

(Δi)_control＝0.171°

由式(13)至式(16)计算周期内每次小推力点火的开关机时刻

t_1,on＝InnerTwoPi(α_i-λ-Θ₀)/ω_e-T_fire/2 (13)

t_1,off＝t_1,on+T_fire (14)

t_2,on＝InnerTwoPi(α_i+π-λ-Θ₀)/ω_e-T_fire/2 (15)

t_2,off＝t_2,on+T_fire (16)

例如，对2021年6月7日的两次点火有：

Θ₀＝261.5293°，亚太6D卫星的定点经度为134.0°E，带入式(13)至式(16)得到t_1,on＝26787.5(s)，t_1,off＝37572.5(s)，t_2,on＝69987.5(s)，t_2,off＝80772.5(s)，对应时刻如下表1所示。

表1亚太6D卫星对应时刻数据表

开机时刻t1,on	开机时刻t1,off	开机时刻t2,on	开机时刻t2,off
				07:26:26.5	09:57:01.5	19:26:26.5	21:57:01.5

Claims (2)

1.地球同步卫星小推力南北保持控制解析计算方法，其特征在于，具体计算过程包括如下步骤：

步骤1、解析计算地球同步卫星倾角矢量的长期摄动变化量；

步骤2、计算电推力器点火序列，解析计算出整个规划周期内每一次电推点火的开关机时刻；

步骤1的具体步骤为：确定倾角矢量为(i_x i_y)，其中i_x＝icosΩ，i_y＝isinΩ，i和Ω分别为地球同步卫星轨道要素的倾角和升交点赤经；计算规划时段t₀到t_f内，倾角矢量的长期摄动变化量(Δi_x)_长期和(Δi_y)_长期；

所述步骤1具体按照以下步骤：

步骤1.1、计算地球非球形引力对倾角矢量的摄动量；

步骤1.2、计算太阳引力对倾角矢量的摄动量；

步骤1.3、计算月球引力对倾角矢量的摄动量；

步骤1.4、将步骤1.1至步骤1.3分别求得的地球非球形引力、太阳引力、月球引力对倾角矢量的摄动量求和，即为地球同步卫星倾角矢量的长期项摄动变化量；

所述步骤1.1中地球非球形引力对倾角矢量的摄动量计算公式为：

其中，n为卫星轨道角速度，a为卫星轨道半长轴，μ_e为地球引力常数，R_e为地球半径，J₂、J₂₂、J₃₁、J₃₃、J₄、J₄₃、λ₂₂、λ₃₁、λ₃₃、λ₄₃是地球非球形四阶主项摄动系数，i_x、i_y为倾角矢量初值；

所述步骤1.2中太阳引力对倾角矢量的摄动量计算公式为：

其中，μ_s为太阳引力常数，r_s为日地距离，Ω_s为地心惯性系下的太阳轨道升交点赤经，i_s为太阳轨道倾角；

所述步骤1.3中月球引力对倾角矢量的摄动量计算公式为：

其中，μ_l为月球引力常数，r_l为地月距离，Ω_l为月球轨道升交点赤经，i_l为月球轨道倾角；

所述步骤1.4中地球同步卫星倾角矢量的长期项摄动变化量由上述三个部分组成，即

(Δi_x)_长期＝(Δi_x)_{非球形，长期}+(Δi_x)_{太阳，长期}+(Δi_x)_{月球，长期} (7)

(Δi_y)_长期＝(Δi_y)_{非球形，长期}+(Δi_y)_{太阳，长期}+(Δi_y)_{月球，长期} (8)。

2.如权利要求1所述的地球同步卫星小推力南北保持控制解析计算方法，其特征在于，所述步骤2的具体过程如下：

步骤2.1、根据步骤1中地球同步卫星倾角矢量的长期项摄动变化量得到一个周期南北保持的倾角矢量目标控制量；一个周期南北保持的倾角矢量目标控制量能够抵消周期内的长期项摄动量，即

(Δi_y)_control＝-(Δi_y)_长期 (9)

(Δi_x)_control＝-(Δi_x)_长期 (10)；

步骤2.2、计算每次小推力点火时长，得到规划周期内的倾角控制量；

根据步骤2.1中一个周期南北保持的倾角矢量目标控制量得到控制方向角，控制方向角α_i＝arctan2[(Δi_y)_control,(Δi_x)_control]，南侧电推力器的点火中心赤经为InnerTwoPi(α_i)，北侧电推力器的点火中心赤经为InnerTwoPi(α_i+π)；

点火序列规划周期为N天，其中的M天进行电推点火，M≤N，为轨道确定、载荷控制预留空间，每天有2次南北控制的时机，电推点火总次数为2M，计算每次小推力点火时长为

T_fire＝2arcsin[(Δi)_controlV_sω_e/acc/M/4]/ω_e (11)

其中，(Δi)_control为规划周期内的倾角控制量，acc为电推力器法向合加速度，V_s为地球同步卫星标称轨道速度，ω_e为地球自转角速度；

电推力器一天内的两次点火时刻t_1,on、t_1,off和t_2,on、t_2,off分别由南侧和北侧点火中心赤经和点火时长计算出来，设当天00:00:00.00时刻的格林尼治恒星时角为Θ₀，卫星平经度为λ；

t_1,on＝InnerTwoPi(α_i-λ-Θ₀)/ω_e-T_fire/2 (13)

t_1,off＝t_1,on+T_fire (14)

t_2,on＝InnerTwoPi(α_i+π-λ-Θ₀)/ω_e-T_fire/2 (15)

t_2,off＝t_2,on+T_fire (16)

计算出来的t_1,on、t_1,off和t_2,on、t_2,off是从当天00:00:00.00时刻开始的积秒，以此类推，解析计算出整个规划周期内每一次电推点火的开关机时刻。