CN114771873B

CN114771873B - 一种超低轨卫星轨道自主精确维持方法

Info

Publication number: CN114771873B
Application number: CN202210302709.1A
Authority: CN
Inventors: 李克行; 罗通; 刘洁; 龚立纲; 王淑一; 姚宁; 魏春岭; 何英姿; 丁建钊; 杨鸣; 黄盘兴; 严晗; 龚宇莲
Original assignee: Beijing Institute of Control Engineering
Current assignee: Beijing Institute of Control Engineering
Priority date: 2022-03-24
Filing date: 2022-03-24
Publication date: 2024-05-03
Anticipated expiration: 2042-03-24
Also published as: CN114771873A

Abstract

本发明提出了一种超低轨卫星轨道自主精确维持方法，采用卡尔曼滤波算法实现轨道平半长轴的精确获取，并针对地球高阶摄动引起的平半长轴波动项采用平均方法进行消除，获取米级精度平半长轴；在无加速度计进行大气阻力加速度测量情况下，通过平半长轴的变化确定大气阻力大小，同时根据确定的大气阻力实时修正补偿轨控推力，解决了超低轨卫星米级精度平半长轴精确获取、轨控推力实时修正的轨道自主精确维持控制问题。

Description

一种超低轨卫星轨道自主精确维持方法

技术领域

本发明属于航天器轨道控制领域，涉及一种超低轨卫星轨道自主精确维持方法。

背景技术

超低轨卫星轨道衰减影响主要来源于大气阻力，当前，低轨道以及超低轨道保持的控制方法主要有两种途径：一是假设参考卫星仅受到地球引力的作用，将超低轨道保持的控制问题转化为卫星编队飞行的队形保持问题进行优化；另一种无阻力轨道控制技术，也就是利用加速度计、大气成分探测器等星载仪器测量大气阻力加速度，通过施加控制力补偿大气阻力的轨道保持方法。

超低轨轨道维持的第一种方法将超低轨道保持的控制问题转化为卫星编队飞行的队形保持问题，引入一颗虚拟卫星作为参考卫星，其中虚拟卫星仅受到地球引力的作用，基于平均轨道根数建立真实卫星相对于参考卫星的相对运动模型，利用线性二次调节器(LQR)可求得编队飞行队形保持的控制律。此方法虚拟参考卫星选择尤为重要，由于超低轨任务只对轨道半长轴控制精度提出较高要求，未对其前后关系和沿飞行方向距离进行约束，在编队控制时由于相对运动方程约束需要对沿飞行方向距离进行约束，将带来额外燃料消耗；此外由于没有相对测量敏感器，平半长轴高阶摄动波动难以精确消除，难以实现几十米级控制精度。超低轨轨道维持的第二种方法利用加速度计、大气成分探测器等星载仪器测量大气阻力加速度，通过施加控制力补偿大气阻力，此种方法额外增加了载荷。

发明内容

本发明解决的技术问题是：基于GPS通道板数据确定的卫星平半长轴在测量噪声的影响下存在千米级误差，无法直接用于卫星轨道维持。采用卡尔曼滤波算法实现轨道平半长轴的精确获取，并针对地球高阶摄动引起的平半长轴波动项采用平均方法进行消除，获取米级精度平半长轴；在无加速度计进行大气阻力加速度测量情况下，通过平半长轴的变化确定大气阻力大小，同时根据确定的大气阻力实时修正补偿轨控推力，解决了超低轨卫星米级精度平半长轴精确获取、轨控推力实时修正的轨道自主精确维持控制问题。

本发明的技术解决方案是：

一种超低轨卫星轨道自主精确维持方法，包括如下步骤：

(1)将当前第i个采样的GPS数据由地固坐标系转换到瞬时真赤道坐标系下，获得在瞬时真赤道坐标系下GPS测量的卫星位置和速度

获得在瞬时真赤道坐标系下GPS测量的卫星位置和速度的方法，具体为：

x＝cos(S_g)·X-sin(S_g)·Y；

y＝sin(S_g)·X+cos(S_g)·Y；

z＝Z；

其中,X,Y,Z和为卫星在地固坐标系下的位置和速度，由星载GPS提供，单位为km和km/s；S_g为格林威治真恒星时角，单位为rad。

(2)若当前为首次采样，则跳转至步骤(4)；反之，则进入步骤(3)；

(3)将上一采样时刻更新后的卫星位置和速度X_fi-1外推到当前采样时刻，获得预估的卫星位置和速度然后进入步骤(4)；

步骤(3)中获得预估的卫星位置和速度的方法，具体为：

Y₃＝h×f(X₂)；X₃＝X_f+Y₃；

Y₄＝h×f(X₃)；

其中，h为当前采样时刻与上一采样时刻的时间间隔，即采样周期，取为GPS传递数据的时间间隔；f(X)为轨道动力学模型，f(X)的具体表达式为：

其中，X＝[X(1) X(2) X(3) X(4) X(5) X(6)]，R_e＝6378.14km为地球平均半径，μ＝3.986005*10⁵km³/s²为地球引力常数，J₂＝0.0010826266835531513为地球带谐项第二阶系数。

(4)若当前为首次采样，令更新后的卫星位置和速度X_f1＝X；否则，根据步骤(1)所述在瞬时真赤道坐标系下GPS测量的卫星位置和速度以及步骤(3)所述预估的卫星位置和速度/>采用卡尔曼滤波算法获得当前第i个采样时刻更新后的卫星位置和速度/>

步骤(4)中更新卫星位置和速度的卡尔曼滤波算法，具体为：

X_fi＝X₀+K(X-X₀)

其中，K为滤波增益矩阵。

K＝PH^T(HPH^T+R)^-1，

P＝Φ_6×6PΦ_6×6 ^T+Q，

其中，Q为模型噪声方差矩阵，R为测量噪声方差矩阵，P为估计均方差矩阵，Φ_6×6为状态转移矩阵，F为雅克比矩阵，h为采样周期，I_6×6为6维单位阵，f(X)为步骤(3)所述轨道动力学模型。

(5)根据步骤(4)所述更新后的卫星位置和速度X_f计算轨道平半长轴a；

步骤(5)中计算轨道平半长轴a可分为如下两步：

第一步，根据卫星状态计算卫星的密切(瞬时)轨道六根数；

第二步，通过消除短周期项，将密切轨道根数转化为轨道平根数。

以上两步的详细计算过程在很多轨道动力学方面的教科书里都有介绍，比如刘林《卫星轨道理论与应用》，第一步可参考2.3节，第二步可参考4.3节；由于其公式较为繁琐且不是本专利的重点，此处不再罗列。

(6)获得当前采样时刻以及前N-1个采样时刻的N个轨道平半长轴a；从而获得N个轨道平半长轴的平均值a₂；一个轨道周期内采样时刻的个数为N；

步骤(6)中一个轨道周期内采样时刻的个数N的确定方法，具体为：

其中，a_t为卫星轨道的目标半长轴，μ＝3.986005×10⁵km³/s²为地球引力常数，T为轨道周期，h为采样周期，表示向下取整。

(7)获得之前第2N-1个采样时刻到之前第N个采样时刻之间N个采样时刻的轨道平半长轴a，并获得平均值a₃；

(8)计算平均值a₂和平均值a₃的差值Δa_Real；

步骤(8)中差值Δa_Real的计算方法，具体为：

Δa_Real＝a₂-a₃。

(9)计算当前采样时刻以及前N-1个采样时刻的N个电推力引起的轨道平半长轴理论变化量的总和Δa_F；

步骤(9)中N个电推力引起的轨道平半长轴理论变化量的总和Δa_F的计算方法，具体为：

其中，Δa_F_j为之前第j个采样时刻的电推力引起的轨道平半长轴理论变化量，F_jet为之前第j个采样时刻的电推力大小(单位为N)，M_sat为之前第j个采样时刻的卫星总质量(单位为kg)，T为步骤(6)所述的轨道周期。

(10)根据步骤(8)差值Δa_Real和步骤(9)N个电推力引起的轨道平半长轴理论变化量的总和Δa_F，确定大气阻力大小F_drag；

步骤(9)中确定大气阻力大小F_drag的方法，具体为：

F_drag＝(Δa_F-Δa_Real×1000)×π×M_sat0/(T²)；

M_{sat 0}为当前采样时刻的卫星总质量。

(11)令下一采样周期内的电推力等于步骤(10)所述大气阻力大小F_drag，在新的采样周期来临后，返回步骤(1)；从而通过将每个采样周期内电推力大小设定为所求大气阻力大小实现轨道自主精确维持。

本发明与现有技术相比的优点在于：

1)无需采用加速度计测量大气阻力加速度，利用星载GPS实时定位数据即可确定大气阻力大小；

2)对卡尔曼滤波所的轨道平半长轴再次在一个轨道周期内求平均，获取了米级精度的轨道平半长轴；

3)令电推力实时等于大气阻尼大小，解决了动态不确定大气密度变化情况下，轨控推力大小自主调节问题。

附图说明

图1为本发明方法的流程图；

图2为步骤(5)计算所得轨道平半长轴a与真实动力学所得轨道平半长轴曲线图；

图3为步骤(10)确定的大气阻力结果曲线图；

图4为轨道半长轴的长期控制效果图；

图5为轨道平半长轴平均值与目标值a_t之间的偏差图。

具体实施方式

超低轨道卫星的大气阻力能够达到几十毫牛，长期积累，对卫星轨道影响较大；基于GPS通道板数据确定的卫星平半长轴在测量噪声的影响下存在千米级误差，无法直接用于卫星轨道维持。

本发明立足于“星上自主”轨道维持基本思想，在无高精度加速度计对大气阻力测量约束下，采用星载GPS数据滤波实时求解轨道，平均处理消除高阶摄动引起的平半长轴波动项，根据平半长轴的变化确定大气阻力大小并实时修正补偿推力，实现了仅依靠星载GPS数据完成超低轨卫星轨道高精度自主维持。

本发明采用图1所示流程完成一种超低轨卫星轨道自主精确维持方法，下面结合具体实例对本发明进行详细说明。

工况条件如下：

卫星阻力面积0.81m²，初始质量为550kg，大气阻力系数Cd＝2.2，初始轨道数据(2020-1-1 00:00:00，北京时间)为：

a(km)

e

i()

Ω()

ω()

M()

瞬时根数

6636.937

1E-6

96.7141

0

90

0

目标轨道半长轴为6646.14km；GPS数据与真实动力学模型的测量精度：位置误差10m(1σ)、速度误差0.2m/s(1σ)，每1秒传递一次数据没有时间延迟。

按照如下步骤可完成卫星轨道精确维持：

(1)将当前第i个采样时刻的GPS数据由地固坐标系转换到瞬时真赤道坐标系下，获得在瞬时真赤道坐标系下GPS测量的卫星位置和速度

(3)将上一采样时刻更新后的卫星位置和速度X_fi-1外推到当前采样时刻，获得预估的卫星位置和速度然后进入步骤4)；

(4)若当前为首次采样，令更新后的卫星位置和速度X_f1＝X；否则，根据步骤(1)所述在瞬时真赤道坐标系下GPS测量的卫星位置和速度以及步骤(3)所述预估的卫星位置和速度/>采用卡尔曼滤波算法获得更新后的卫星位置和速度/>

(5)根据步骤(4)所述更新后的卫星位置和速度计算轨道平半长轴a；

步骤(1)～(5)通过滤波消除GPS数据的噪声误差，滤波效果如图2所示：虚线表示真实动力学模型下的轨道平半长轴，实线表示利用存在测量误差的GPS数据滤波后获得的轨道平半长轴，由图可知，滤波结果可在5min以内收敛至真实动力学结果附近。

(6)获得当前采样时刻以及前N-1个采样时刻的N个轨道平半长轴a，并计算N个轨道平半长轴的平均值a₂；

一个轨道周期内采样时刻的个数

其中，目标半长轴a_t＝6646.14km，地球引力常数μ＝3.986005×10⁵km³/s²，采样周期h＝1s，T为轨道周期，表示向下取整；因此N计算为5381。

(8)计算平均值a_2i和平均值a_3i的差值Δa_Real＝a₂-a₃；

Δa_F的计算方法具体为：

其中，Δa_F_j为之前第j个采样时刻的电推力引起的轨道平半长轴理论变化量，F_{jet j}为之前第j个采样时刻的电推力大小(单位为N)，M_{sat j}为之前第j个采样时刻的卫星总质量(单位为kg)。

(10)确定大气阻力大小F_drag＝(Δa_F-Δa_Real×1000)×π×M_{sat 0}/(T^2)；

步骤(6)～(10)完成大气阻力F的实时确定，图3给出了大气阻力随时间变化的曲线，结果表明确定的大气阻力在一定范围内上下波动，与稀薄大气密度分布不均匀的特性一致。

(11)令下一采样周期内的电推力等于步骤(10)所述大气阻力大小F，在新的采样周期来临后，返回步骤(1)。

由上述步骤可完成轨道自主精确维持，该工况下轨道半长轴的长期控制效果如图4所示：卫星的初始半长轴略高于目标值a_t，经过一段自然衰减过程实现维轨构型建立，当轨道半长轴的平均值继续减小至目标值a_t后开启沿速度方向的电推进，并令推力值等于当前确定的大气阻力，结果表明卫星的轨道半长轴可长期维持在目标值附近。图5给出了某时段轨道半长轴的平均值与目标值a_t之间的偏差，其偏差量在米级范围内。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。

Claims

1.一种超低轨卫星轨道自主精确维持方法，其特征在于，包括如下步骤：

(5)根据步骤(4)所述更新后的卫星位置和速度，计算轨道平半长轴a；

(6)获得当前采样时刻以及前N-1个采样时刻的N个轨道平半长轴a；从而获得N个轨道平半长轴的平均值a₂；

(8)计算平均值a₂和平均值a₃的差值Δa_Real；

(11)令下一采样周期内的电推力等于步骤(10)所述大气阻力大小F_drag，在新的采样周期来临后，返回步骤(1)；将每个采样周期内电推力大小设定为所求大气阻力，完成轨道自主精确维持方法。

2.根据权利要求1所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(1)中，获得在瞬时真赤道坐标系下GPS测量的卫星位置和速度的方法，具体为：

x＝cos(S_g)·X-sin(S_g)·Y；

y＝sin(S_g)·X+cos(S_g)·Y；

z＝Z；

其中,X,Y,Z和为卫星在地固坐标系下的位置和速度，由星载GPS提供；S_g为格林威治真恒星时角。

3.根据权利要求1或2所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(3)中，获得预估的卫星位置和速度的方法，具体为：

Y₁＝h×f(X_fi)；

Y₂＝h×f(X₁)；

Y₃＝h×f(X₂)；X₃＝X_fi+Y₃；

Y₄＝h×f(X₃)；

其中，h为当前采样时刻与上一采样时刻的时间间隔，f(X)为轨道动力学模型，h等于GPS传递数据的时间间隔。

4.根据权利要求3所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述f(X)的具体表达式为：

其中，X＝[X(1) X(2) X(3) X(4) X(5) X(6)]，R_e＝为地球平均半径，μ为地球引力常数，J₂为地球带谐项第二阶系数。

5.根据权利要求4所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(4)中，更新卫星位置和速度的卡尔曼滤波算法，具体为：

X_fi＝X₀+K(X-X₀)，

其中，K为滤波增益矩阵。

6.根据权利要求5所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述滤波增益矩阵，具体为：

K＝PH^T(HPH^T+R)^-1，

P＝Φ_6×6PΦ_6×6 ^T+Q，

其中，Q为模型噪声方差矩阵，R为测量噪声方差矩阵，P为估计均方差矩阵，Φ_6×6为状态转移矩阵，F为雅克比矩阵，h为采样周期，I_6×6为6维单位阵。

7.根据权利要求1所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(6)中，N为一个轨道周期内采样时刻的个数。

8.根据权利要求7所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(6)中，个数N的确定方法，具体为：

其中，a_t为卫星轨道的目标半长轴，μ为地球引力常数，T为轨道周期，h为采样周期，表示向下取整。

9.根据权利要求1所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(8)中，差值Δa_Real的计算方法，具体为：

Δa_Real＝a₂-a₃。

10.根据权利要求1所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(9)中，N个电推力引起的轨道平半长轴理论变化量的总和Δa_F的计算方法，具体为：

其中，Δa_F_j为之前第j个采样时刻的电推力引起的轨道平半长轴理论变化量，F_{jet j}为之前第j个采样时刻的电推力大小，M_{sat j}为之前第j个采样时刻的卫星总质量，T为轨道周期，h为采样周期。

11.根据权利要求1所述的一种超低轨卫星轨道自主精确维持方法，其特征在于：所述步骤(10)中，确定大气阻力大小F_drag的方法，具体为：

F_drag＝(Δa_F-Δa_Real×1000)×π×M_{sat 0}/(T²)；

其中，T为轨道周期，M_{sat 0}为当前采样时刻的卫星总质量。