一种考虑长周期项的低轨星座系统高精度相位保持方法与流程

1.本发明涉及一种用于低轨星座系统的高精度相位保持方法。更具体地，涉及一种考虑长周期项的低轨星座系统高精度相位保持方法。


背景技术：

2.利用低轨星座提供覆盖全球的卫星互联网服务。星座作为一个庞大的空间系统，稳定的构型是其发挥正常功能的基础。由于各卫星的入轨偏差及所受摄动力存在差异，星座构型趋于发散，从而影响星座的覆盖特性及任务实现。此外，星座系统将面临卫星间以及卫星与空间碎片的碰撞威胁。为使低轨星座满足覆盖特性要求与空间轨道安全要求，低轨星座的卫星相位需要保持在一定的较小区间内(例如参考相位的
±
0.1
°
以内)。
3.在低轨星座高精度相位保持研究方面，现有的技术基于半长轴一阶变化率主要由大气阻力导致、近似为常值的假设，建立了平相位角偏差随时间呈二次函数关系的模型，计算了相位保持轨道高度改变量，实现了低轨星座系统相位保持在
±
0.1
°
以内。但该方法没有考虑日月三体引力等摄动因素导致的平相位角长周期项的影响，因此该方法虽能将相位保持在
±
0.1
°
的区间，但对相位保持区间的利用并不充分，相位保持轨道机动较为频繁。
4.因此，需要提供一种考虑长周期项的低轨星座系统高精度相位保持方法，将地球非球形、大气阻力、日月三体引力以及太阳光压摄动导致的平相位角长周期项作为平轨值的修正项，从平相位角中分离，再进行平相位角偏差系数的拟合，精确预测本周期以及下周期平相位角偏差，通过迭代计算平相位角偏差达到第二阈值的时刻，实现对平相位角偏差第二阈值的精准控制；通过迭代计算半长轴改变量大小，实现对下个相位保持周期平相位角偏差第一阈值的精准控制。本发明所提方法可以在包含地球非球形、大气阻力、日月三体引力、太阳光压的复杂摄动环境中，实现低轨星座系统的高精度相位保持。


技术实现要素：

5.本发明的目的在于提供一种考虑长周期项的低轨星座系统高精度相位保持方法。
6.为达到上述目的，本发明采用下述技术方案：
7.本发明提供一种考虑长周期项的低轨星座系统高精度相位保持方法，该方法的步骤包括：
8.s1、确定卫星参考平相位角随时间的变化关系；
9.s2、每隔固定时长采样一次卫星的位置和速度，并将其转化为实际平相位角，并计算相应时刻的平相位角长周期项；
10.s3、将所述实际平相位角与所述参考平相位角作差，计算平相位角偏差，并用所述平相位角长周期项做修正，获取平相位角偏差系数；
11.s4、预报下一采样时刻平相位角偏差，若超出第二阈值，计算所述平相位角偏差超出第二阈值时刻；
12.s5、根据所述超出第二阈值时刻与所述平相位角偏差系数确定相位保持轨道机动
半长轴改变量初始解；
13.s6、根据本相位保持周期平相位角偏差系数预测下一相位保持周期平相位角偏差，并用长周期项做修正，确定相位保持轨道机动半长轴改变量修正解；
14.s7、在平相位角偏差值达到第二阈值时刻施加速度增量，将轨道机动结束时刻作为下一控制周期的初始时刻，并重复s2至s7，直到完成卫星整个相位保持任务周期的相位保持。
15.优选地，所述第一阈值为平相位角偏差区间左端点的值，所述第二阈值为所述平相位角偏差区间右端点的值。
16.优选地，所述步骤s1包括：
17.根据地球中心引力场、地球非球形摄动和日月三体引力摄动确定参考平相位角随时间的变化率；
18.根据卫星在卫星相位保持任务周期的初始时刻的参考平相位角和所述参考平相位角随时间的变化率，确定卫星在某一时刻的参考平相位角。
19.优选地，所述步骤s2包括：
20.将平均轨道要素代入由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的长周期项公式，获得所述平相位角长周期项。
21.优选地，所述步骤s3包括：
22.用所述平相位角长周期项做修正，获取平相位角偏差修正值，根据平相位角偏差修正值随时间的变化趋势拟合计算平相位角偏差系数。
23.优选地，所述平相位角偏差系数为平相位角偏差修正值随时间变化的二次多项式系数，包括，初始时刻的平相位角偏差、初始时刻的半长轴偏差导致平相位角漂移一阶项、半长轴随时间线性变化导致平相位角漂移二阶项。
24.优选地，所述步骤s4还包括：
25.若平相位角偏差预报值未超出第二阈值则重新进行步骤s2至s3更新平相位角偏差系数，直至该时刻的平相位角偏差超出第二阈值；
26.若平相位角偏差预报值超出第二阈值，引进不足一个采样周期的自由漂移时长计算所述平相位角偏差超出第二阈值时刻。
27.优选地，所述步骤s5包括：
28.根据所述平相位角偏差系数，计算初始时刻实际半长轴与参考半长轴的偏差和初始时刻到所述平相位角偏差超出第二阈值时刻卫星的轨道半长轴随时间的一阶线性变化率，从而确定所述平相位角偏差超出第二阈值时刻的实际半长轴与参考半长轴的偏差；
29.根据平相位角偏差区间以及平相位角漂移二阶项，估算下一相位保持两次轨道高度抬升之间的时长；
30.根据两次轨道高度抬升之间的半长轴漂移量以及平相位角偏差超出第二阈值时刻的实际半长轴与参考半长轴的偏差，确定轨道机动的半长轴改变量初始解。
31.优选地，所述步骤s6包括：
32.预报下一相位保持周期日月及卫星轨道要素，计算对应的相位角偏差长周期项，获得所述下一相位保持周期平相位角偏差，确定使所述下一相位保持周期平相位角偏差最小值恰好达到第一阈值的相位保持轨道机动半长轴改变量修正解。
33.优选地，不同相位保持周期中日月及卫星轨道要素不同，由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的平相位角长周期项也不同。
34.本发明的有益效果如下：
35.1、本发明的考虑长周期项的低轨星座系统高精度相位保持方法，将平相位角长周期项从平相位角中分离，再进行平相位角偏差系数的拟合，精确预测本周期以及下周期平相位角偏差，从而精确控制平相位角偏差达到上下边界(阈值)，实现了在考虑地球非球形、大气阻力、日月三体引力、太阳光压的复杂摄动环境中低轨星座系统高精度相位保持。
36.2、本发明的一种低轨星座系统高精度相位保持方法，在卫星位置速度确定精度不高、而且定位数据采样率较低的情况下，可以实现高精度相位保持，降低了高精度相位保持对卫星定轨精度以及定轨数据采样存储要求，具有很大的工程应用价值。
附图说明
37.下面结合附图对本发明的具体实施方式作进一步详细的说明。
38.图1示出本发明提出的一种考虑长周期项的低轨星座系统高精度相位保持方法流程图。
39.图2示出各摄动因素导致的平相位角偏差短周期项时间历程。
40.图3示出各摄动因素导致的平相位角偏差长周期项时间历程。
41.图4示出相位保持任务周期内，相位角偏差瞬轨值、平轨值时间历程。
具体实施方式
42.为了更清楚地说明本发明，下面结合优选实施例和附图对本发明做进一步的说明。附图中相似的部件以相同的附图标记进行表示。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
43.本发明第一个实施例对本发明提出的一种考虑长周期项的低轨星座系统高精度相位保持方法进行详细说明。
44.本发明提出将地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的平相位角长周期项作为平轨值的修正项，从平相位角中分离，再进行平相位角偏差系数的拟合，精确预测本周期以及下周期平相位角偏差，通过迭代计算平相位角偏差达到第二阈值的时刻，实现对平相位角偏差第二阈值的精准控制；通过迭代计算半长轴改变量大小，实现对下个相位保持周期平相位角偏差第一阈值的精准控制。本发明所提方法可以在包含地球非球形、大气阻力、日月第三体引力、太阳光压的复杂摄动环境中，实现低轨星座系统的高精度相位保持。
45.采用绝对相位保持方法，星座内的每颗卫星，相位保持策略相同，仅参考相位不同，通过本发明的控制策略，使每颗卫星的实际平相位角跟踪其参考平相位角，即可实现整个星座的相位保持。由于每颗卫星的相位保持策略相同，后续不再针对星座展开描述，直接针对卫星展开描述。
46.如图1所示，该方法主要包括如下7个步骤：
47.s1、确定卫星参考平相位角随时间的变化关系；
48.卫星的任一速度位置状态对应着一组轨道根数，这组轨道根数称为瞬时轨道根
数。这组轨道根数处于时刻变化之中，用它作为一段时间内卫星的轨道根数，会产生较大的误差。采用平均化思路的解析法将瞬时轨道根数的变化分解为性质不同的三类：长期项，长周期项和短周期项。一般地，平均轨道根数(简称平根数、平轨值等)指的是将短周期项平均掉后，所对应的一组轨道根数，平均轨道根数能在一个更高层次上反应卫星的运动规律。通常在卫星轨道设计、轨道机动等任务中，都以平均轨道根数作为设计和控制的参数。
49.本发明中的平相位角λ，定义为近地点幅角平轨值ω与平近点角平轨值 m之和，即λ＝ω m。参考平相位角随时间的变化考虑地球中心引力场、地球非球形摄动与日月三体引力摄动。将参考平相位角λ随时间的变化率记为n
ref
，由地球中心引力场导致的平相位角变化率记为nr，由地球非球形摄动导致的平相位角变化率记为nz，由太阳引力导致的平相位角变化率记为ns，由月球引力导致的平相位角变化率记为nm，则：
50.n
ref
＝nr nz ns nmꢀꢀ
(1)
51.在一个可能的实现方式中，假设任务开始时刻，卫星的平相位角已经调整到保持区间内，将整个任务的开始时刻作为第一个相位保持周期的初始时刻。将卫星在初始时刻t0的参考平相位角记为λ
ref
(t0)，则卫星在时刻t的参考平相位角为：
52.λ
ref
(t)＝λ
ref
(t0) n
ref
·
(t-t0)
ꢀꢀ
(2)
53.s2、每隔固定时长采样一次卫星的位置和速度，并将其转化为实际平相位角，并计算相应时刻的平相位角长周期项；
54.从初始时刻开始，在同一相位保持周期内每隔固定时长采样一次卫星的位置和速度，并将其转化为实际平相位角，并计算相应时刻的平相位角长周期项。
55.在一个可能的实现方式中，将位置和速度采样间隔记为δt，从初始时刻开始获取l 1个采样点的位置矢量r和速度矢量v的值，其中l为某一相位保持周期的初始采样间隔数，在l
·
δt时间段内，无需对平相位角偏差是否超出第二阈值进行判断。通过位置、速度到轨道要素的转化方法，与瞬时轨道要素到平均轨道要素的转换方法，计算卫星在时刻t＝t0 l
·
δt(l∈[0，1，2，
…
，l])的实际平相位角λ
real
(t)。
[0056]
在一个可能的实现方式中，本发明要将平均轨道要素代入由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的长周期项公式，获得相应时刻由各摄动因素导致的平相位角长周期项之和λ1(t)。
[0057]
s3、将所述实际平相位角与所述参考平相位角作差，计算平相位角偏差，并用所述平相位角长周期项做修正，获取平相位角偏差系数；
[0058]
将相应时刻卫星的实际平相位角与其相应的参考平相位角作差，获得时刻t＝t0 l
·
δt的平相位角偏差：
[0059]
δλ(t)＝λ
real
(t)-λref(t)
ꢀꢀ
(3)
[0060]
在一个可能的实现方式中，本发明提出用平相位角长周期项做修正，获取平相位角偏差修正值再拟合计算平相位角偏差系数。计算以上相应时刻的平相位角长周期项，获得平相位角偏差修正值：
[0061][0062]
λ
l
(t0)是初始时刻t0的各摄动因素导致的平相位角长周期项之和。去除由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的长周期项影响后，在大气阻力的作用
下，半长轴随时间近似为一阶线性变化，平相位角偏差修正值随时间为二阶非线性变化，平相位角偏差修正值随时间t的变化趋势近似为：
[0063][0064]
其中：δλ0是初始时刻t0的平相位角偏差，δλ1是初始时刻t0的半长轴偏差导致平相位角漂移一阶项，δλ2是半长轴随时间线性变化导致平相位角漂移二阶项，这三项即为平相位角偏差随时间变化的二次多项式系数。
[0065]
根据时刻t＝t0 l
·
δt(l∈[0，1，2，
…
，l])，卫星的平相位角偏差δλ(t)数值，以及表达式(5)，利用最小二乘法拟合平相位角偏差系数δλ0、δλ1和δλ2。
[0066]
s4、预报下一采样时刻平相位角偏差，若超出第二阈值，计算所述平相位角偏差超出第二阈值时刻；
[0067]
确定平相位角偏差系数δλ0、δλ1和δλ2的值后，预报时刻t＝t0 (l h)
·
δt (h＝1，2，3，
…
)的平相位角偏差修正值同时预报卫星与日月在时刻 t＝t0 (l h)
·
δt轨道要素，获得平相位角长周期项的预报值，从而获得时刻t＝t0 (l h)
·
δt的平相位角偏差预报值δλ
p
(t)。
[0068][0069]
在一个可能的实现方式中，如果时刻t＝t0 (l h)
·
δt的平相位角偏差预报值未超出第二阈值δλ
max
，则在时间段[t0 (l h-1)
·
δt t0 (l h)
·
δt]内不施加相位保持轨道控制，采用步骤s2的方法获得时刻t0 (l h)
·
δt的实际平相位角与平相位角长周期项，并利用步骤s3的方法在考虑t0 (l h)
·
δt时刻平相位角偏差修正值后，对平相位角偏差随时间变化的二次多项式系数(偏差系数) 进行更新。
[0070]
如果时刻t＝t0 (l h)
·
δt的平相位角偏差预报值δλ
p
(t)超出第二阈值δλ
max
，则令：
[0071]
δλ
p
(t)＝δλ
max
ꢀꢀ
(7)
[0072]
结合式(6)解出平相位角偏差达到第二阈值δλ
max
的时刻tf[0073]
tf＝t0 (l h-1)
·
δt δtnꢀꢀ
(8)
[0074]
式中δtn＜δt。
[0075]
在一个可能的实现方式中，本发明通过引进不足一个采样周期的自由漂移时长δtn，使得对于平相位角偏差的第二阈值控制达到更加精准。由于计算长周期项[λ
l
(t)-λ
l
(t0)]涉及到对卫星轨道要素的预报，采用迭代法确定平相位角偏差超出第二阈值时刻tf。
[0076]
s5、根据所述超出第二阈值时刻与所述平相位角偏差系数确定相位保持轨道机动半长轴改变量初始解；
[0077]
确定了卫星从时刻t0到时刻tf的平相位角偏差修正值随时间变化的系数δλ0、δλ1和δλ2。通过下式可以计算初始时刻t0实际半长轴与参考半长轴的偏差δa(t0)。
[0078]
[0079]
式中，μ为地球引力常数，ac为参考半长轴。通过下式可以计算这段时间段内卫星的轨道半长轴随时间的一阶线性变化率
[0080][0081]
从而可以确定时刻tf的实际半长轴与参考半长轴的偏差δa(tf)。
[0082][0083]
根据平相位角偏差区间[δλ
min δλ
max
]，以及平相位角漂移二阶项δλ2，可以估算下一相位保持两次轨道高度抬升之间的时长为：
[0084][0085]
在两次轨道高度抬升之间的半长轴漂移量为：
[0086][0087]
所以相位保持轨道机动所需半长轴改变量初始解为：
[0088][0089]
s6、根据本相位保持周期平相位角偏差系数预测下一相位保持周期平相位角偏差，并用长周期项做修正，确定相位保持轨道机动半长轴改变量修正解；
[0090]
根据本相位保持周期平相位角偏差系数预测下一相位保持周期平相位角偏差，并用长周期项做修正。确定使下一相位保持周期平相位角偏差最小值恰好达到第一阈值的相位保持轨道机动半长轴改变量修正解。
[0091]
在一个可能的实现方式中，本发明考虑到不同相位保持周期日月及卫星轨道要素的不同，由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的平相位角长周期项也不同。因此，步骤s5获得的相位保持轨道机动半长轴改变量只是初始解，要精确地控制下一相位保持周期的平相位角偏差达到第一阈值，还应对以上获得的半长轴改变量进行修正。
[0092]
在一个可能的实现方式中，假设本相位保持周期为第k相位保持周期，在式中定义的平相位角偏差系数δλ0，δλ1，δλ2分别记为δλ
0，k
，δλ
1，k
，δλ
2，k
，则根据式计算出第k相位保持周期结束时刻半长轴偏差值δak(tf)，根据式计算出第k相位保持周期半长轴改变量初始解δa
′k，设半长轴修正量为δam，则第k相位保持周期半长轴改变量修正解为δak，满足：
[0093]
δak＝δa
′k δamꢀꢀ
(15)
[0094]
则第k 1相位保持周期初始时刻半长轴偏差值δa
k 1
(t0)，
[0095]
δa
k 1
(t0)＝δak(tf) δakꢀꢀ
(16)
[0096]
预测第k 1相位保持周期平相位角偏差系数δλ
0，k 1
，δλ
1，k 1
，δλ
2，k 1
为：
[0097][0098]
预报第k 1相位保持周期日月及卫星轨道要素，计算对应的相位角偏差长周期项[λ
l
(t)-λ
l
(t0)]，获得第k 1相位保持周期平相位角偏差δλ(k 1，t)
[0099]
δλ(k 1，t)＝δλ
0，k 1
δλ
1，k 1
(t-t0) δλ
2，k 1
(t-t0)2 [λ
l
(t)-λ
l
(t0)]
ꢀꢀ
(18)
[0100]
通过迭代法确定相位保持轨道机动半长轴改变量修正解δak，使下一相位保持周期平相位角偏差最小值min[δλ(k 1，t)]恰好达到第一阈值δλ
min
。
[0101]
s7、在平相位角偏差值达到第二阈值时刻施加速度增量，将轨道机动结束时刻作为下一控制周期的初始时刻，并重复步骤s2至s7，直到完成卫星整个相位保持任务周期的相位保持。
[0102]
在一个可能的实现方式中，对于近圆轨道，对应半长轴改变量δak需要的横向速度增量为：
[0103][0104]
相位保持轨道机动推进器工作时长为：
[0105][0106]
式中，m为卫星质量，f
p
为推进器推力值。在步骤s4获得的平相位角偏差达到第二阈值时刻tf＝t0 (l h-1)
·
δt δtn施加切向推力，获得横向速度增量δv
t
。相位保持周期的总时长定义为(l h-1)
·
δt δtn δtm。将轨道机动结束时刻t＝t0 (l h-1)
·
δt δtn δtm作为下一控制周期的初始时刻，并重复 s2至s7，直到完成卫星整个相位保持任务周期的相位保持。
[0107]
本发明另一个实施例结合附图2-4对本发明做进一步说明。
[0108]
如图2所示，给出了各摄动因素导致的平相位角偏差短周期项时间历程。由地球带谐主项导致的平相位角偏差瞬轨值振幅约为5.5
×
10-2
(
°
)，由地球田谐主项导致的平相位角偏差瞬轨值振幅约为5.3
×
10-3
(
°
)，由大气阻力导致的平相位角偏差瞬轨值振幅约为9.0
×
10-8
(
°
)，由太阳引力导致的平相位角偏差瞬轨值振幅约为5.8
×
10-6
(
°
)，由月球引力导致的平相位角偏差瞬轨值振幅约为 1.5
×
10-5
(
°
)，由太阳光压导致的平相位角偏差瞬轨值振幅约为3.6
×
10-6
(
°
)。
[0109]
假设任务要求将相位角偏差的瞬轨值保持在
±
0.1
°
的区间内，综合以上各摄动因素导致的短周期项以及瞬平转换中未考虑到的高阶短周期项，且保留一定控制裕度，将平相位角偏差的平轨值保持在区间[-0.03
°ꢀ
0.03
°
]内，即δλ
min
＝-0.03
°
，δλ
max
＝0.03
°
。
[0110]
本实施例中考虑了地球非球形摄动、大气阻力摄动、日月三体引力摄动与太阳光压摄动。地球非球形摄动采用wgs84-egm96模型，degree和order 设为21，大气阻力摄动采用us 1976 standard atmospheric density model。卫星参考轨道高度取为800km，初始轨
道要素见表1。卫星的面质比取为0.0175，阻力系数取为2.2，相位保持任务周期取为90天，定轨数据采样间隔取为2 小时，三轴位置确定误差取为5m(1σ)，三轴速度确定误差取为0.02m/s(1σ)。场景历元时间为1 jan 2022 00:00:00:000 utcg。
[0111]
表1初始平均轨道根数
[0112]
轨道根数值a7178.137kme0.001i87
°
ω20
°
ω90
°
f10
°
λ99.98
°
[0113]
确定卫星的参考平相位角随时间的变化率与参考平相位角随时间的变化关系。
[0114]
在卫星标称半长轴、偏心率、轨道倾角以及历元时间给定的情况下，参考平相位角随时间的变化率可以根据考虑的地球中心引力场和地球非球形摄动低阶项、以及日月三体引力摄动长期项直接解析获得：
[0115]nref
＝5.941
×
10-2
(
°
/s)
[0116]
假设任务开始时刻，卫星的平相位角偏差正好位于上边界δλ
max
处，将整个任务的开始时刻作为第一个相位保持周期的初始时刻。则卫星在初始时刻t0的参考平相位角λ
ref
(t0)应为99.95
°
。卫星在时刻t的参考平相位角为：
[0117]
λ
ref
(t)＝99.95
°
5.941
×
10-2
·
t(
°
)
[0118]
从初始时刻开始，在同一相位保持周期内每隔固定时长采样一次卫星的位置和速度，并将其转化为实际平相位角，并计算相应时刻的平相位角长周期项。
[0119]
取l＝30，从初始时刻开始，每隔δt＝2h采样一次卫星的位置矢量r和速度矢量v的值，并通过位置、速度到轨道要素的转化方法，与瞬时轨道要素到平均轨道要素的转换方法，计算卫星在时刻t＝l
·
7200sec(l∈[0，1，2，
…
，30])的实际平相位角λ
real
(l
·
7200sec)。将平均轨道要素代入由地球非球形、大气阻力、日月三体引力、太阳光压摄动导致的长周期项公式，获得相应时刻平相位角长周期项之和λ
l
(l
·
7200sec)。
[0120]
如图3所示，给出了由各种摄动因素导致的长周期项时间历程。由地球带谐主项导致的平相位角偏差瞬轨值振幅约为3.0
×
10-5
(
°
)，由大气阻力导致的平相位角偏差瞬轨值振幅约为2.0
×
10-16
(
°
)，由太阳引力导致的平相位角偏差瞬轨值振幅约为6.0
×
10-3
(
°
)，由月球引力导致的平相位角偏差瞬轨值振幅约为3.0
×
10-3
(
°
)，由太阳光压导致的平相位角偏差瞬轨值振幅约为 5.0
×
10-6
(
°
)。
[0121]
可以看出，由各摄动综合作用导致的平相位角偏差瞬轨值振幅约为 1.0
×
10-2
(
°
)，主要是由太阳引力与月球引力导致的长周期项叠加而成，因此日月三体引力导致的长周期项在各摄动因素导致的长周期项中占主导作用。也由于日月三体引力导致的长周期项已经与平相位角偏差的保持精度同等量级，因此在考虑日月影响的摄动模型中，用长周期项来修正平相位角偏差是非常有必要的。
[0122]
计算平相位角偏差，并用平相位角长周期项做修正，获取平相位角偏差系数。
[0123]
将以上相应时刻卫星的实际平相位角与其相应的参考平相位角作差，获得时刻t＝l
·
7200sec的平相位角偏差δλ(l
·
7200sec)，并用以上相应时刻的平相位角长周期项[λ
l
(l
·
7200sec)-λ
l
(0
·
7200sec)]做修正，获得平相位角偏差修正值表2给出了不同仿真时刻的平相位角偏差、长周期修正项与平相位角偏差修正项，其中δλ(l
·
δt)，λ
l
(t)-λ
l
(t0)，的单位为度 (
°
)。
[0124]
根据平相位角偏差修正值随时间t的变化规律，可以计算平相位角偏差系数为：
[0125][0126]
预报下一采样时刻平相位角偏差，并确定平相位角偏差超出第二阈值时刻。
[0127]
根据平相位角偏差随时间变化的因子参数，可以预测t＝31.7200s的平相位角为-0.0011
°
，没有达到边界值0.03
°
，在时间段[(31)
·
7200 (32)
·
7200]内不施加相位保持轨道控制，继续获得时刻(31)
·
7200s的实际平相位角偏差与平相位角长周期项，并更新平相位角偏差系数。
[0128]
采用上述步骤依次预报时刻(32)
·
7200s、(33)
·
7200s的平相位角偏差，直到相应时刻的平相位角偏差超出漂移边界(第二阈值)δλ
max
，确定从时刻t0＝0 开始的卫星无控制自由漂移时间为t
end
＝(30 159)
·
7200s 4479s，即h＝160，δtn＝4479s，以及该段时间内的平相位角偏差随时间变化的因子参数。
[0129][0130]
表2不同仿真时刻的平相位角偏差、长周期修正项与平相位角偏差修正项
[0131][0132]
确定相位保持轨道机动半长轴改变量初始解。
[0133]
计算初始时刻t0实际半长轴与参考半长轴的偏差δa(t0)。
[0134]
δa(0)≈14.886m
[0135]
计算这段时间段内卫星的轨道半长轴随时间的一阶线性变化率
[0136][0137]
从而可以确定时刻t
end
的实际半长轴与参考半长轴的偏差δa(tf)：
[0138]
δa(tf)≈-14.333m
[0139]
根据平相位角偏差区间[-0.03
°ꢀ
0.03
°
]以及平相位角漂移二阶项δλ2，可以估算下一相位保持两次轨道高度抬升之间的时长为：
[0140]
δt＝15.548day
[0141]
在两次轨道高度抬升之间的半长轴漂移量为：
[0142]
δad＝-28.748m
[0143]
所以相位保持轨道机动所需半长轴改变量为：
[0144]
δa
′
＝28.707m
[0145]
根据本相位保持周期平相位角偏差系数预测下一相位保持周期平相位角偏差，并用长周期项做修正。确定使下一相位保持周期平相位角偏差最小值恰好达到第一阈值的相位保持轨道机动半长轴改变量修正解。
[0146]
假设本相位保持周期为第k相位保持周期，平相位角偏差系数为δλ
0，k
，δλ
1，k
，δλ
2，k
，确定δλ
2，k
＝2.321
×
10-15
/s2，确定第k相位保持周期结束时刻半长轴偏差值δak(tf)＝-14.333m，以及第k相位保持周期半长轴改变量初始解δa
′k＝28.707m。采用迭代法计算半长轴修正解，首先假设半长轴修正量δam＝0，则半长轴改变量修正解为：
[0147]
δak＝28.707m
[0148]
则第k 1相位保持周期初始时刻半长轴偏差值：
[0149]
δa
k 1
(t0)＝14.374m
[0150]
预测第k 1相位保持周期平相位角偏差系数δλ
0，k 1
，δλ
1，k 1
，δλ
2，k 1
为：
[0151][0152]
预报第k 1相位保持周期日月及卫星轨道要素，计算对应的相位角偏差长周期项[λ
l
(t)-λ
l
(t0)]，获得第k 1相位保持周期平相位角偏差δλ(k 1，t)，其中最小值为：
[0153]
min[δλ(k 1，t)]＝-0.291
°
＞-0.03
°
[0154]
表明施加的半长轴改变量略小，应增大半长轴改变量，经过几步迭代，确定半长轴修正量δam＝0.106m，则半长轴改变量修正解为：
[0155]
δak＝28.813m
[0156]
第k 1相位保持周期初始时刻半长轴偏差值：
[0157]
δa
k 1
(t0)＝14.480m
[0158]
更新后的第k 1相位保持周期平相位角偏差系数δλ
0，k 1
，δλ
1，k 1
，δλ
2，k 1
为：
[0159][0160]
预测第k 1相位保持周期平相位角偏差最小值恰好为-0.03
°
。
[0161]
在计算获得的平相位角偏差达到上边界时刻施加速度增量。将轨道机动结束时刻作为下一控制周期的初始时刻，并重复上述步骤，直到完成卫星整个任务周期的相位保持。
[0162]
对于近圆轨道，对应半长轴改变量δak需要的横向速度增量为：
[0163]
δv
t
＝1.496
×
10-2
m/s
[0164]
取卫星质量为200kg，推进器推力值为80mn，则相位保持轨道机动推进器工作时长
为：
[0165]
δtm＝37.390s
[0166]
在获得的平相位角偏差达到上边界时刻tf＝(30 159)
·
7200s 4479s施加切向推力，获得横向速度增量δv
t
。相位保持周期的总时长定义为 t
end
＝(30 159)
·
7200s 4479s 37.390s。将轨道机动结束时刻t＝t
end
作为下一控制周期的初始时刻，并重复上述步骤，直到完成整个任务周期的相位保持。
[0167]
图4给出了6个相位保持周期的仿真结果，可以看出相位角偏差瞬轨值保持在
±
0.1
°
以内，相位角偏差平轨值基本控制在
±
0.03
°
以内，本发明所提方法实现了复杂摄动环境下低轨星座系统高精度相位保持。
[0168]
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。