一种恒定推力模式持续机动卫星的轨道预报方法

摘要

本发明涉及一种恒定推力模式持续机动卫星的轨道预报方法，包括：获取待轨道预报卫星的多组轨道平根数；建立半长轴摄动函数，根据多组轨道平根数，得到每一条轨道根数对应的历元时间差和半长轴差，根据历元时间差和半长轴差，计算得到半长轴摄动函数的系数；建立摄动方程，根据轨道平根数对摄动方程进行积分，得到每一条轨道根数对应的相位角差，构建新的平近点角摄动函数，根据每一条轨道根数对应的历元时间差和相位角差，计算得到新的平近点角摄动函数的系数；利用新的平近点角摄动函数对摄动方程进行修正，得到修正后的摄动方程，利用修正后的摄动方程对待轨道预报卫星进行轨道预报。本发明的方法提高了持续抬轨或降轨卫星的轨道预报精度。

权利要求

1.一种恒定推力模式持续机动卫星的轨道预报方法，其特征在于，包括：步骤1：获取待轨道预报卫星的多组轨道平根数；

步骤2：建立半长轴摄动函数，根据所述多组轨道平根数，得到每一条轨道根数对应的历元时间差和半长轴差，根据所述历元时间差和半长轴差，计算得到所述半长轴摄动函数的系数；所述步骤2包括：步骤2.1：建立半长轴摄动函数如下：

式中，a表示半长轴，t表示时间，A表示半长轴摄动函数的系数；

步骤2.2：依次计算得到第k条轨道根数与第N条轨道根数之间的历元时间差Δtk，其中，Δtk＝tk‑tN，k＝1,2,…,N，tk表示第k条轨道根数的历元时间，tN表示第N条轨道根数的历元时间，N表示轨道根数的个数；

步骤2.3：依次计算得到第k条轨道根数与第N条轨道根数之间的半长轴差Δak，其中，Δak＝ak‑aN，k＝1,2,…,N，ak表示第k条轨道根数的半长轴的平根数，aN表示第N条轨道根数的半长轴的平根数；

步骤2.4：根据N条轨道根数对应的历元时间差和半长轴差序列{(Δtk,ak)|k＝1,2,…,N}，利用最小二乘法求解得到所述半长轴摄动函数的系数A；

步骤3：建立摄动方程，根据所述轨道平根数对所述摄动方程进行积分，得到每一条轨道根数对应的相位角差，构建新的平近点角摄动函数，根据每一条轨道根数对应的历元时间差和相位角差，计算得到所述新的平近点角摄动函数的系数；所述步骤3包括：步骤3.1：建立摄动方程如下：

式中，e表示偏心率，i表示倾角，Ω表示升交点赤经，ω表示近地点幅角，M表示平近点角，rE表示地球半径，n表示平均运动速度，p表示轨道半通径，J2表示地球引力场二阶带谐系数，J3表示地球引力场三阶带谐系数；

步骤3.2：从第N条轨道根数开始，以预设步长对所述摄动方程进行积分，依次计算得到每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值；

步骤3.3：根据所述多组轨道平根数以及每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值，计算得到每一条轨道根数对应的相位角差，其中，第k条轨道根数的相位角差Δθk按照下式计算得到：Δθk＝ωk+Mk‑ωck‑Mck，k＝1,2,…,N；

式中，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，ωck表示第k条轨道根数的近地点幅角的积分值，Mck表示第k条轨道根数的平近点角的积分值；

步骤3.4：构建新的平近点角摄动函数如下：

式中，B表示新的平近点角摄动函数的第一系数，C表示新的平近点角摄动函数的第二系数；

步骤3.5：根据N条轨道根数对应的历元时间差和相位角差序列{(Δtk,Δθk)|k＝1,

2,…,N}，利用最小二乘法求解得到所述新的平近点角摄动函数的第一系数B和第二系数C；

步骤4：利用所述新的平近点角摄动函数对所述摄动方程进行修正，得到修正后的摄动方程，利用所述修正后的摄动方程对待轨道预报卫星进行轨道预报。

2.根据权利要求1所述的恒定推力模式持续机动卫星的轨道预报方法，其特征在于，在所述步骤1中，对所述待轨道预报卫星，按照历元时间顺序获取N条轨道根数对应的轨道平根数，其中，第k条轨道根数的轨道平根数表示为：σk(ak,ek,ik,Ωk,ωk,Mk)，k＝1,2,…,N；

式中，ak表示第k条轨道根数的半长轴的平根数，ek表示第k条轨道根数的偏心率的平根数，ik表示第k条轨道根数的倾角的平根数，Ωk表示第k条轨道根数的升交点赤经的平根数，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，N表示轨道根数的个数。

3.根据权利要求1所述的恒定推力模式持续机动卫星的轨道预报方法，其特征在于，所述步骤4包括：步骤4.1：利用所述新的平近点角摄动函数对所述摄动方程进行修正，得到修正后的摄动方程为：步骤4.2：对所述修正后的摄动方程进行积分，得到预报时刻的所述待轨道预报卫星的轨道平根数积分结果，根据所述轨道平根数积分结果得到所述待轨道预报卫星在该预报时刻的空间位置及速度。

说明书

一种恒定推力模式持续机动卫星的轨道预报方法

技术领域

[0001] 本发明属于航天测量与控制技术领域，具体涉及一种恒定推力模式持续机动卫星的轨道预报方法。

背景技术

[0002] 近年来，随着高度集成化和自动化航天技术快速发展、发射成本持续降低、市场需求不断扩大，巨型星座的研发和部署掀起前所未有的热潮，星座的发展逐渐从低轨卫星通信星座转向低轨巨型星座，巨型星座一般指卫星部署数量超过100颗的星座。

[0003] 为节约星座建造成本，这些卫星一般采用一箭多星的方式发射入轨，入轨高低一般较低，入轨后采用持续小推力的方式逐步抬升轨道至正常服务轨道。卫星持续抬轨或降轨时，可能采用恒定推力模式，即每天达到固定的轨道高度变化量，这种机动方式也可称为恒定半长轴变率模式的机动方式。

[0004] 虽然轨道高度日变化量相同，但是具体的轨控策略是未知的。对于小推力卫星而言，持续变轨期间每日点火时长超过万秒，推力大小可大致估算，但点火动时段分配未知，甚至可能是固定变轨量模式下的随机点火策略，导致轨道预报误差很大。

发明内容

[0005] 为了解决现有技术中存在的上述问题，本发明提供了一种恒定推力模式持续机动卫星的轨道预报方法。本发明要解决的技术问题通过以下技术方案实现：

[0006] 本发明提供了一种恒定推力模式持续机动卫星的轨道预报方法，包括：

[0007] 步骤1：获取待轨道预报卫星的多组轨道平根数；

[0008] 步骤2：建立半长轴摄动函数，根据所述多组轨道平根数，得到每一条轨道根数对应的历元时间差和半长轴差，根据所述历元时间差和半长轴差，计算得到所述半长轴摄动函数的系数；

[0009] 步骤3：建立摄动方程，根据所述轨道平根数对所述摄动方程进行积分，得到每一条轨道根数对应的相位角差，构建新的平近点角摄动函数，根据每一条轨道根数对应的历元时间差和相位角差，计算得到所述新的平近点角摄动函数的系数；

[0010] 步骤4：利用所述新的平近点角摄动函数对所述摄动方程进行修正，得到修正后的摄动方程，利用所述修正后的摄动方程对待轨道预报卫星进行轨道预报。

[0011] 在本发明的一个实施例中，在所述步骤1中，对所述待轨道预报卫星，按照历元时间顺序获取N条轨道根数对应的轨道平根数，其中，第k条轨道根数的轨道平根数表示为：

[0012] σk(ak,ek,ik,Ωk,ωk,Mk)，k＝1,2,…,N；

[0013] 式中，ak表示第k条轨道根数的半长轴的平根数，ek表示第k条轨道根数的偏心率的平根数，ik表示第k条轨道根数的倾角的平根数，Ωk表示第k条轨道根数的升交点赤经的平根数，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，N表示轨道根数的个数。

[0014] 在本发明的一个实施例中，所述步骤2包括：

[0015] 步骤2.1：建立半长轴摄动函数如下：

[0016]

[0017] 式中，a表示半长轴，t表示时间，A表示半长轴摄动函数的系数；

[0018] 步骤2.2：依次计算得到第k条轨道根数与第N条轨道根数之间的历元时间差Δtk，其中，Δtk＝tk‑tN，k＝1,2,…,N，tk表示第k条轨道根数的历元时间，tN表示第N条轨道根数的历元时间，N表示轨道根数的个数；

[0019] 步骤2.3：依次计算得到第k条轨道根数与第N条轨道根数之间的半长轴差Δak，其中，Δak＝ak‑aN，k＝1,2,…,N，ak表示第k条轨道根数的半长轴的平根数，aN表示第N条轨道根数的半长轴的平根数；

[0020] 步骤2.4：根据N条轨道根数对应的历元时间差和半长轴差序列{(Δtk,ak)|k＝1,2,…,N}，利用最小二乘法求解得到所述半长轴摄动函数的系数A。

[0021] 在本发明的一个实施例中，所述步骤3包括：

[0022] 步骤3.1：建立摄动方程如下：

[0023]

[0024] 式中，e表示偏心率，i表示倾角，Ω表示升交点赤经，ω表示近地点幅角，M表示平近点角，rE表示地球半径，n表示平均运动速度，p表示轨道半通径，J2表示地球引力场二阶带谐系数，J3表示地球引力场三阶带谐系数；

[0025] 步骤3.2：从第N条轨道根数开始，以预设步长对所述摄动方程进行积分，依次计算得到每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值；

[0026] 步骤3.3：根据所述多组轨道平根数以及每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值，计算得到每一条轨道根数对应的相位角差，其中，第k条轨道根数的相位角差Δθk按照下式计算得到：

[0027] Δθk＝ωk+Mk‑ωck‑Mck，k＝1,2,…,N；

[0028] 式中，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，ωck表示第k条轨道根数的近地点幅角的积分值，Mck表示第k条轨道根数的平近点角的积分值；

[0029] 步骤3.4：构建新的平近点角摄动函数如下：

[0030]

[0031] 式中，B表示新的平近点角摄动函数的第一系数，C表示新的平近点角摄动函数的第二系数；

[0032] 步骤3.5：根据N条轨道根数对应的历元时间差和相位角差序列{(Δtk,Δθk)|k＝1,2,…,N}，利用最小二乘法求解得到所述新的平近点角摄动函数的第一系数B和第二系数C。

[0033] 在本发明的一个实施例中，所述步骤4包括：

[0034] 步骤4.1：利用所述新的平近点角摄动函数对所述摄动方程进行修正，得到修正后的摄动方程为：

[0035]

[0036] 步骤4.2：对所述修正后的摄动方程进行积分，得到预报时刻的所述待轨道预报卫星的轨道平根数积分结果，根据所述轨道平根数积分结果得到所述待轨道预报卫星在该预报时刻的空间位置及速度。

[0037] 与现有技术相比，本发明的有益效果在于：

[0038] 本发明的一种恒定推力模式持续机动卫星的轨道预报方法，基于半长轴和相位角修正的摄动方程进行轨道预报，提高了持续抬轨或降轨卫星的轨道预报精度，克服了持续机动卫星控制策略未知导致轨道预报误差偏大的问题。

[0039] 上述说明仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，而可依照说明书的内容予以实施，并且为了让本发明的上述和其他目的、特征和优点能够更明显易懂，以下特举较佳实施例，并配合附图，详细说明如下。

附图说明

[0040] 图1是本发明实施例提供的一种恒定推力模式持续机动卫星的轨道预报方法的示意图。

具体实施方式

[0041] 为了进一步阐述本发明为达成预定发明目的所采取的技术手段及功效，以下结合附图及具体实施方式，对依据本发明提出的一种恒定推力模式持续机动卫星的轨道预报方法进行详细说明。

[0042] 有关本发明的前述及其他技术内容、特点及功效，在以下配合附图的具体实施方式详细说明中即可清楚地呈现。通过具体实施方式的说明，可对本发明为达成预定目的所采取的技术手段及功效进行更加深入且具体地了解，然而所附附图仅是提供参考与说明之用，并非用来对本发明的技术方案加以限制。

[0043] 实施例一

[0044] 请参见图1，图1是本发明实施例提供的一种恒定推力模式持续机动卫星的轨道预报方法的示意图，如图所示，本实施例的恒定推力模式持续机动卫星的轨道预报方法，包括：

[0045] 步骤1：获取待轨道预报卫星的多组轨道平根数；

[0046] 在本实施例中，对待轨道预报卫星，即需要进行轨道预报的卫星，按照历元时间顺序获取N条轨道根数对应的轨道平根数，其中，第k条轨道根数的轨道平根数表示为：

[0047] σk(ak,ek,ik,Ωk,ωk,Mk)，k＝1,2,…,N  (1)；

[0048] 式中，ak表示第k条轨道根数的半长轴的平根数，ek表示第k条轨道根数的偏心率的平根数，ik表示第k条轨道根数的倾角的平根数，Ωk表示第k条轨道根数的升交点赤经的平根数，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，N表示轨道根数的个数。

[0049] 其中，所有轨道根数对应的历元时间顺序为t1

[0050] 在一个可选地实施方式中，轨道平根数可以是由已知的TLE根数转换，也可基于实际测量数据进行短弧轨道确定得到。

[0051] 步骤2：建立半长轴摄动函数，根据多组轨道平根数，得到每一条轨道根数对应的历元时间差和半长轴差，根据历元时间差和半长轴差，计算得到半长轴摄动函数的系数；

[0052] 在一个可选地实施方式中，步骤2包括：

[0053] 步骤2.1：建立半长轴摄动函数如下：

[0054]

[0055] 式中，a表示半长轴，t表示时间，A表示半长轴摄动函数的系数；

[0056] 步骤2.2：依次计算得到第k条轨道根数与第N条轨道根数之间的历元时间差Δtk，其中，Δtk＝tk‑tN，k＝1,2,…,N，tk表示第k条轨道根数的历元时间，tN表示第N条轨道根数的历元时间，N表示轨道根数的个数；

[0057] 步骤2.3：依次计算得到第k条轨道根数与第N条轨道根数之间的半长轴差Δak，其中，Δak＝ak‑aN，k＝1,2,…,N，ak表示第k条轨道根数的半长轴的平根数，aN表示第N条轨道根数的半长轴的平根数；

[0058] 步骤2.4：根据N条轨道根数对应的历元时间差和半长轴差序列{(Δtk,ak)|k＝1,2,…,N}，利用最小二乘法求解得到半长轴摄动函数的系数A。

[0059] 步骤3：建立摄动方程，根据轨道平根数对摄动方程进行积分，得到每一条轨道根数对应的相位角差，构建新的平近点角摄动函数，根据每一条轨道根数对应的历元时间差和相位角差，计算得到新的平近点角摄动函数的系数；

[0060] 在一个可选地实施方式中，步骤3包括：

[0061] 步骤3.1：建立摄动方程如下：

[0062]

[0063] 式中，e表示偏心率，i表示倾角，Ω表示升交点赤经，ω表示近地点幅角，M表示平近点角，rE表示地球半径，n表示平均运动速度，p表示轨道半通径，J2表示地球引力场二阶带谐系数，J3表示地球引力场三阶带谐系数；

[0064] 在本实施例中，摄动方程中的 即为步骤2中得到的半长轴摄动函数。

[0065] 步骤3.2：从第N条轨道根数开始，以预设步长对摄动方程进行积分，依次计算得到每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值；

[0066] 在本实施例中，从第N条轨道根数开始，以‑60秒为步长对摄动方程进行积分，积分至每一条轨道根数对应的历元时间时，得到该条轨道根数对应的近地点幅角的积分值和平近点角的积分值。

[0067] 步骤3.3：根据多组轨道平根数以及每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值，计算得到每一条轨道根数对应的相位角差，其中，第k条轨道根数的相位角差Δθk按照下式计算得到：

[0068] Δθk＝ωk+Mk‑ωck‑Mck，k＝1,2,…,N  (4)；

[0069] 式中，ωk表示第k条轨道根数的近地点幅角的平根数，Mk表示第k条轨道根数的平近点角的平根数，ωck表示第k条轨道根数的近地点幅角的积分值，Mck表示第k条轨道根数的平近点角的积分值；步骤3.4：构建新的平近点角摄动函数如下：

[0070]

[0071] 式中，B表示新的平近点角摄动函数的第一系数，C表示新的平近点角摄动函数的第二系数；

[0072] 在本实施例中，新的平近点角摄动函数是通过对摄动方程中的原始平近点角摄动函数 引入关于时间t的函数Bt+C得到的。

[0073] 步骤3.5：根据N条轨道根数对应的历元时间差和相位角差序列{(Δtk,Δθk)|k＝1,2,…,N}，利用最小二乘法求解得到新的平近点角摄动函数的第一系数B和第二系数C。

[0074] 步骤4：利用新的平近点角摄动函数对摄动方程进行修正，得到修正后的摄动方程，利用修正后的摄动方程对待轨道预报卫星进行轨道预报。

[0075] 在一个可选地实施方式中，步骤4包括：

[0076] 步骤4.1：利用新的平近点角摄动函数对摄动方程进行修正，得到修正后的摄动方程为：

[0077]

[0078] 步骤4.2：对修正后的摄动方程进行积分，得到预报时刻的待轨道预报卫星的轨道平根数积分结果，根据轨道平根数积分结果得到待轨道预报卫星在该预报时刻的空间位置及速度。

[0079] 假设对待轨道预报卫星进行预报时刻为tp的轨道预报，在一个可选地实施方式中，首先，利用第N条轨道根数的轨道平根数以及修正后的摄动方程进行积分，则得到积分结束时刻，即预报时刻tp的轨道平根数积分结果为(ap,ep,ip,Ωp,ωp,Mp)；其次，将该预报(s) (s) (s) (s) (s) (s)时刻tp的轨道平根数积分结果转换为对应的瞬根数(ap ,ep ,ip ,Ωp ,ωp ,Mp )，其中，各轨道元素的瞬根数为对应平根数加上摄动力的短周期项；最后，将瞬根数转换为对应的空间位置及速度(xp,yp,zp,vxp,vyp,vzp)，得到待轨道预报卫星在预报时刻tp的轨道预报结果。

[0080] 本发明实施例的恒定推力模式持续机动卫星的轨道预报方法，基于半长轴和相位角修正的摄动方程进行轨道预报，提高了持续抬轨或降轨卫星的轨道预报精度，克服了持续机动卫星控制策略未知导致轨道预报误差偏大的问题。

[0081] 实施例二

[0082] 本实施例以具体星链‑1872为例，对实施例一的恒定推力模式持续机动卫星的轨道预报方法进行说明，星链‑1872卫星在2022年9月9日至9月17日期间持续降轨，轨道高度每日降低约4.5公里，每日降轨量恒定，可认为是恒定推力模型下的持续降轨。对其进行轨道预报的实施步骤如下：

[0083] 步骤一：获取星链‑1872卫星2022年9月10日至9月12日共6条TLE(Two Line Elements)根数，并按时间顺序排序t1

[0084] 表1.星链‑1872卫星9月10日至9月12日的TLE根数

[0085]

[0086] 将6条TLE根数转换为对应的轨道平根数，结果如表2所示。

[0087] 表2. 6条TLE根数对应的轨道平根数

[0088]

[0089] 步骤二：建立半长轴摄动函数 根据表2的轨道平根数得到6条轨道根数对应的历元时间差和半长轴差序列如表3所示，根据历元时间差和半长轴差序列计算得到半长轴摄动函数的系数A为‑4537.872，得到半长轴摄动函数为

[0090] 表3.历元时间差和半长轴差序列

[0091] Δtk 值(单位：天) Δak 值(单位：米)Δt1 ‑2.660 Δa1 12076.054
Δt2 ‑2.183 Δa2 9908.965
Δt3 ‑1.676 Δa3 7477.131
Δt4 ‑1.076 Δa4 4865.119
Δt5 ‑0.681 Δa5 3049.219
Δt6 0.000 Δa6 0.000

[0092] 步骤三：建立如公式(3)的摄动方程，其中，摄动方程中的 即为步骤二中得到的半长轴摄动函数 从第6条轨道根数开始，以‑60秒为步长对摄动方程进行积分，得到每一条轨道根数对应的历元时间的近地点幅角的积分值和平近点角的积分值，根据积分值和表2的平根数，得到6条轨道根数对应的历元时间差和相位角差序列如表4所示，构建新的平近点角摄动函数 根据历元时间差和相位角差序列
计算得到新的平近点角摄动函数的第一系数B为‑0.003，第二系数C为2.854，得到新的平近点角摄动函数为

[0093] 表4.历元时间差和相位角差序列

[0094]Δtk 值(单位：天) Δθk 值(单位：°)
Δt1 ‑2.660 Δθ1 ‑7.545
Δt2 ‑2.183 Δθ2 ‑6.180
Δt3 ‑1.676 Δθ3 ‑4.775
Δt4 ‑1.076 Δθ4 ‑3.052
Δt5 ‑0.681 Δθ5 ‑1.952
Δt6 0.000 Δθ6 0.000

[0095] 步骤四：根据新的平近点角摄动函数对摄动方程进行修正，得到修正后的摄动方程为：

[0096]

[0097] 假设预报24小时后，即2022‑9‑13 18:54:42时刻卫星的空间位置，利用第6条轨道平根数和修正后的摄动方程，积分得到2022‑9‑13 18:54:42时刻卫星的轨道平根数如表5所示。

[0098] 表5.2022‑9‑13 18:54:42时刻的平根数

[0099] ap(m) ep ip(°) Ωp(°) ωp(°) Mp(°)6888914.269 0.000992 52.929 275.110 87.711 244.800

[0100] 对各轨道元素的平根数加上摄动力的短周期项转换得到该时刻的瞬根数如表6所示。

[0101] 表6. 2022‑9‑13 18:54:42时刻的瞬根数

[0102]

[0103]

[0104] 将表6所示瞬根数转换为位置速度，得到卫星在2022‑9‑13 18:54:42时刻的轨道预报结果，如表7所示。

[0105] 表7.预报24小时后的惯性系位置速度

[0106] 历元时刻(UTC时) 2022‑9‑13 18:54:42X(m) ‑1373605.394
Y(m) ‑6254917.137
Z(m) ‑2546994.184
Vx(m/s) 4359.695052
Vy(m/s) ‑3142.974581
Vz(m/s) 5381.039858

[0107] 应当说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的物品或者设备中还存在另外的相同要素。“连接”或者“相连”等类似的词语并非限定于物理的或者机械的连接，而是可以包括电性的连接，不管是直接的还是间接的。

[0108] 以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。