CN104142686A - 一种卫星自主编队飞行控制方法 - Google Patents

Abstract

本发明公开了一种卫星自主编队飞行控制方法，通过轨道平根数差进行编队飞行控制，由于轨道平根数差较准确的反映了卫星之间相对运动的长期趋势，这种方法可以较好的控制相对运动的长期变化。本发明通过设计轨道平面内的平半长轴差控制策略，采用分区间设置控制目标的方式，保证在控制区间内的漂移速度较小；在控制区间外时，能以较快的速度回到控制区间内。本发明通过轨控使用多次小脉冲喷气、姿控使用动量轮的方式，减少姿态喷气控制对轨道的影响，提高轨道控制执行精度。

Description

一种卫星自主编队飞行控制方法

技术领域

本发明涉及一种卫星自主编队飞行控制方法，适用于卫星长期编队飞行的控制。

背景技术

描述编队飞行的两颗卫星之间的相对运动时，一般使用线性化之后的C-W方程来描述。对编队进行控制时，也采用C-W方程来计算控制量。但近期的一些研究成果表明，C-W方程在描述相对运动时，对于编队的长期运动趋势的描述并不准确。近10年来发展起来的轨道根数差(又成为相对轨道根数)描述方式能更加好的描述长期相对运动趋势，可用于编队飞行的轨道控制。但现有文献中对于如何从测量数据计算轨道根数差的研究较少。若用两颗卫星的绝对轨道参数直接相减获得轨道根数差，由于绝对轨道参数的误差较大，两星的参数相减后的误差会更大，难以用于精确的编队飞行控制。

一般来说，在进行轨道控制时，轨控推力器会产生的干扰力矩影响卫星的姿态控制。卫星的姿态控制可以使用动量轮或姿控推力器两种方式。若轨控干扰力矩的持续时间较长，动量轮无法吸收轨控干扰力矩累积形成的角动量。因此，一般轨控时使用姿控推力器进行姿态控制。但姿控推力器同时也产生一定的速度增量，从而影响轨控执行的精度。

传统的轨道控制中，每次轨控的时间较长，轨控推力器产生干扰较大，轨控时的姿态控制一般使用推力器进行，而姿控推力器产生的速度增量对轨控的影响较小，可以作为误差考虑。但是在编队飞行中，两颗星处于稳定编队飞行状态下时，两星的轨道非常接近，轨道控制的速度增量非常小。若使用推力器进行姿控，姿控产生的速度增量会对严重影响轨控速度增量的执行精度。

发明内容

本发明所要解决的技术问题是：克服现有技术的不足，提供了一种卫星自主编队飞行控制方法，可以在燃料消耗非常少的情况下实现了卫星编队飞行的自主控制。

本发明包括如下技术方案：

一种卫星自主编队飞行控制方法，包括如下步骤：

(1)根据A星绝对轨道GPS数据获得A星绝对位置、速度；根据差分GPS数据获得A星和B星的相对位置、速度；根据A星绝对位置、速度与所述相对位置、速度得到B星的绝对位置、速度；

(2)根据A星绝对位置、速度计算A星瞬时轨道根数，根据B星绝对位置、速度计算B星瞬时轨道根数；A星瞬时轨道根数包括半长轴a_A,偏心率e_A，倾角i_A，升交点赤经Ω_A，近地点幅角ω_A，平近点角M_A，纬度幅角u_A以及幅角和λ_A，λ_A＝ω_A+M_A；

(3)根据A星瞬时轨道根数计算A星平根数，根据B星瞬时轨道根数计算B星平根数；A星平根数与B星平根数的计算方法相同；A星平根数包括平半长轴平偏心率矢量和平倾角和平升交点赤经B星平根数包括平半长轴平偏心率矢量和平倾角和平升交点赤经

(4)将A星平根数和B星平根数进行处理得到两星的平根数差，两星的平根数差包括平半长轴差平偏心率矢量差和平倾角矢量差和 $\mathrm{\Δ}{\tilde{i}}_y;$ $\mathrm{\Δ}\tilde{a}={\stackrel{\‾}{a}}_A-{\stackrel{\‾}{a}}_B,\mathrm{\Δ}{\tilde{e}}_x={\stackrel{\‾}{\mathrm{\ξ}}}_A-{\stackrel{\‾}{\mathrm{\ξ}}}_B,\mathrm{\Δ}{\tilde{e}}_y={\stackrel{\‾}{\mathrm{\η}}}_B-{\stackrel{\‾}{\mathrm{\η}}}_A,\mathrm{\Δ}{\tilde{i}}_x=\sin \left({\stackrel{\‾}{i}}_A\right)\sin ({\stackrel{\‾}{\mathrm{\Ω}}}_A-{\stackrel{\‾}{\mathrm{\Ω}}}_B),$ $\mathrm{\Δ}{\tilde{i}}_y=-\sin \left({\stackrel{\‾}{i}}_A\right)\cos \left({\stackrel{\‾}{i}}_B\right)\cos ({\stackrel{\‾}{\mathrm{\Ω}}}_A-{\stackrel{\‾}{\mathrm{\Ω}}}_B)+\sin \left({\stackrel{\‾}{i}}_B\right)\cos \left({\stackrel{\‾}{i}}_A\right);$

(5)对平根数差进行一阶滤波获得滤波后的平根数差，所述滤波后的平根数差包括平半长轴差Δa、偏心率矢量差Δe_x和Δe_y、倾角矢量差Δi_x和Δi_y；

(6)判断是否需要进行轨道平面外控制，当大于阈值Δi_max时，转入步骤(7)进行轨道平面外控制计算；当不大于阈值Δi_max时，不进行轨道平面外控制，转入步骤(10)；

同时判断是否需要进行轨道平面内控制，当x_r<X2或x_r>X3时，转入步骤(8)进行轨道平面内控制计算；当X2≤x_r≤X3时，不进行轨道平面内控制，转入步骤(10)；x_r为A星和B星在飞行方向上的相对位置，X2，X3为位置阈值；

(7)轨道平面外控制计算的步骤如下：

(7.1)根据公式计算平面外控制所需要的速度增量Δv_y；为平均轨道角速度，μ为地球引力常数；

(7.2)根据公式u_y＝arctan2(Δi_x,-Δi_y)计算平面外控制所需要的纬度幅角u_y；

(7.3)计算开机时间t_y和最终的轨道平面外控制所需要的速度增量 $\mathrm{\&Delta;}{\hat{v}}_y;$

令Δu₁＝u_y-u_A，

$\mathrm{\&Delta;}{\hat{u}}_1=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_1,2\mathrm{\&pi;}),$

${\mathrm{\&Delta;u}}_2=\mathrm{\&Delta;}{\hat{u}}_1+\mathrm{\&pi;},$

$\mathrm{\&Delta;}{\hat{u}}_2=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_2,2\mathrm{\&pi;})$

如果 $\mathrm{\&Delta;}{\hat{u}}_1<\mathrm{\&Delta;}{\hat{u}}_2,$ 则令 $t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_1}{{\stackrel{\&OverBar;}{n}}_B},$ t为当前时间，令 $\mathrm{\&Delta;}{\hat{v}}_y={\mathrm{\&Delta;v}}_y;$ 否则，令 $t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_2}{{\stackrel{\&OverBar;}{n}}_B},\mathrm{\&Delta;}{\hat{v}}_y=-{\mathrm{\&Delta;v}}_y;$

(7.4)转入步骤(9)；

(8)轨道平面内控制计算的步骤如下：

(8.1)根据如下公式计算两个脉冲速度增量Δv_x1和Δv_x2，以及两个脉冲对应的轨道纬度幅角和

${\mathrm{\&Delta;v}}_a=(1/4){\stackrel{\&OverBar;}{n}}_B({\mathrm{\&Delta;a}}_{\mathrm{new}}-\mathrm{\&Delta;a}),{\mathrm{\&Delta;a}}_{\mathrm{new}}=\left\{\begin{array}{cc}-{\mathrm{\&Delta;a}}_1& x_r<X1\\ -{\mathrm{\&Delta;a}}_2& X1\&le;x_r<X2\\ {\mathrm{\&Delta;a}}_2& X3<x_r\&le;X4\\ {\mathrm{\&Delta;a}}_1& x_r>X4\end{array}\right.,$ Δa₁为第一半长轴差控制目标、Δa₂为第二半长轴差控制目标，Δa₁>Δa₂，X1、X4为位置阈值，X1<X2<X3<X4，

${\mathrm{\&Delta;v}}_e=(1/4){\stackrel{\&OverBar;}{n}}_B{\stackrel{\&OverBar;}{a}}_B\left(\sqrt{{\mathrm{\&Delta;e}}_x^2+{\mathrm{\&Delta;e}}_y^2}\right),$ Δv_a、Δv_e为中间变量；

Δv_x1＝Δv_a+Δv_e，

Δv_x2＝Δv_a-Δv_e，

u_x1＝arctan2(-Δe_y,-Δe_x)，

${\hat{u}}_{x1}=\mathrm{mod}\_2\mathrm{PI}(u_{x1},2\mathrm{\&pi;}),$

$u_{x2}={\hat{u}}_{x1}+\mathrm{\&pi;},$ u_x1、u_x2为中间变量，

${\hat{u}}_{x2}=\mathrm{mod}\_2\mathrm{PI}(u_{x2},2\mathrm{\&pi;});$

(8.2)确定最终的第一脉冲的速度增量和第二脉冲的速度增量计算先执行的速度增量对应的轨道纬度幅角与当前轨道纬度幅角的差Δu；

如果( ${\hat{u}}_{x1}<\mathrm{\&pi;}$ 且( $u_A>{\hat{u}}_{x2}$ 或 $u_A\&le;{\hat{u}}_{x1}$ ))或( ${\hat{u}}_{x1}\&GreaterEqual;\mathrm{\&pi;}$ 且 $u_A>{\hat{u}}_{x2}$ 且 $u_A\&le;{\hat{u}}_{x1}$ )，则令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;u}={\hat{u}}_{x1}-u_A;$

否则，令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;u}={\hat{u}}_{x2}-u_A;$

(8.3)令获得最终的纬度幅角的差根据公式计算第一脉冲的开机时间t_x1，根据公式计算第二脉冲的开机时间t_x2；

(8.4)转入步骤(9)；

(9)速度增量脉冲执行，

根据每个方向计算的速度增量、轨控推力器的推力大小以及卫星质量确定相应方向的轨控推力器的开机时长；

根据开机时长确定轨道控制和姿态控制方式，当开机时长大于设定值时，姿态控制使用姿控推力器；在计算的开机时间，打开轨控推力器，轨控推力器连续喷气，使用加速度计测量速度增量，进行速度增量的累积，当速度增量累积到大于计算的速度增量时关闭轨控推力器；当开机时长小于等于设定值时，姿态控制使用动量轮，在计算的开机时间，打开轨控推力器，轨控推力器在开机后的每个控制周期内执行一个短时间的喷气，一直到喷气时间累积到大于开机时长时关闭轨控推力器；

(10)结束。

根据A星瞬时轨道根数计算A星平根数方法包括如下步骤：

a.平根数赋初值

令， $\stackrel{\&OverBar;}{a}=a_A$

$\stackrel{\&OverBar;}{\mathrm{\&xi;}}=e_A\cos {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&eta;}}=-e_A\sin {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{i}=i_A$

$\stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\mathrm{\&Omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\mathrm{\&lambda;}}_A,$

其中为平半长轴中间量，和为平偏心率矢量中间量，为平倾角中间量，为平升交点赤经中间量，为λ_A的平根数中间量；

b.短周期项计算，短周期项包括半长轴短周期项a_s，倾角短周期项i_s，升交点赤经短周期项Ω_s，偏心率矢量的短周期项ξ_s和η_s，幅角和的短周期项λ_s，

$\stackrel{\&OverBar;}{p}=\stackrel{\&OverBar;}{a}(1-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}^2-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}^2),$ 是轨道半通径

其中R_e＝6378.137km为地球赤道半径，J₂＝0.001623945为地球非球形J2项引力系数，是中间变量，

$a_s=\frac{{\tilde{A}}_2}{\stackrel{\&OverBar;}{a}}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

$i_s=\frac{{\tilde{A}}_2}{4{\stackrel{\&OverBar;}{p}}^2}\sin \left(2\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

${\mathrm{\&Omega;}}_s=\frac{{\tilde{A}}_2}{2{\stackrel{\&OverBar;}{p}}^2}\cos \left(\stackrel{\&OverBar;}{i}\right)\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

${\mathrm{\&xi;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\cos \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[1-1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]+\frac{7}{12}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\cos \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\},$

${\mathrm{\&eta;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\sin \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[-1+1.75{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]-\frac{7}{12}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\sin \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\},$

${\mathrm{\&lambda;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}[-0.5+1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right);$

c.平根数计算，得到一步迭代更新后的平半长轴偏心率矢量和平倾角平升交点赤经以及变量

${\stackrel{\&OverBar;}{a}}_n=a_A-a_s$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n=e_A\cos {\mathrm{\&omega;}}_A-{\mathrm{\&xi;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n=-e_A\sin {\mathrm{\&omega;}}_A-{\mathrm{\&eta;}}_s$

${\stackrel{\&OverBar;}{i}}_n=i_A-i_s$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n={\mathrm{\&Omega;}}_A-{\mathrm{\&Omega;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n={\mathrm{\&lambda;}}_A-{\mathrm{\&lambda;}}_s$

如果( $|{\stackrel{\&OverBar;}{a}}_n-\stackrel{\&OverBar;}{a}|<{\mathrm{\&epsiv;}}_a$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&xi;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&xi;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&eta;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&eta;}}$ 且 $|{\stackrel{\&OverBar;}{i}}_n-\stackrel{\&OverBar;}{i}|<{\mathrm{\&epsiv;}}_i$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&Omega;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&Omega;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&lambda;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}$ )，ε_a、ε_ξ、ε_η、ε_i、ε_Ω、ε_λ为阈值；转步骤e；

否则，令 $\stackrel{\&OverBar;}{a}={\stackrel{\&OverBar;}{a}}_n,\stackrel{\&OverBar;}{\mathrm{\&xi;}}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,\stackrel{\&OverBar;}{\mathrm{\&eta;}}={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n,\stackrel{\&OverBar;}{i}={\stackrel{\&OverBar;}{i}}_n,\stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n,\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n,$ 转步骤b；

e.输出平根数，

令 ${\stackrel{\&OverBar;}{a}}_A={\stackrel{\&OverBar;}{a}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n$

${\stackrel{\&OverBar;}{i}}_A={\stackrel{\&OverBar;}{i}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n$

本发明与现有技术相比的有益效果是：

(1)本发明通过轨道平根数差进行编队飞行控制，由于轨道平根数差较准确的反映了卫星之间相对运动的长期趋势，这种方法可以较好的控制相对运动的长期变化。

(2)本发明通过设计轨道平面内的平半长轴差控制策略，采用分区间设置控制目标的方式，保证在控制区间内的漂移速度较小；在控制区间外时，能以较快的速度回到控制区间内。

(3)本发明通过轨控使用多次小脉冲喷气、姿控使用动量轮的方式，减少姿态喷气控制对轨道的影响，提高轨道控制执行精度。

(4)本发明通过上述三点的措施，可以提高编队飞行控制精度，降低编队构型的漂移速度，更加节约燃料消耗。

附图说明

图1为本发明方法的实施步骤流程图。

图2为平根数迭代计算流程。

图3为轨道平面内半长轴差的控制策略示意图。

图4为轨道平面和轨道坐标系示意图。

具体实施方式

多颗卫星组成编队飞行能够完成单颗卫星不能完成功能。编队飞行需要卫星之间保持一定的距离范围。实际中，当两颗卫星的初始轨道有偏差、卫星受到的摄动作用的也略有不同，编队的构型会发生变化。此时，需要设计编队构型的保持控制策略，使卫星间的相对位置满足任务的要求。

卫星的瞬时轨道根数描述了在某一个瞬间卫星的轨道状态。轨道平根数理论将卫星轨道根数的变化按照长期变化、短周期变化、长周期变化等不同的特征进行分解。其中的长期变化部分和长周期变化部分表达了轨道在较长时间内的变化特征，这部分包含在轨道平根数中；短周期部分表述了轨道参数在一个轨道周期之内的变化特征。轨道平根数加上短周期部分等于瞬时轨道根数。对于编队飞行的两个卫星来说，它们轨道平根数差就反应了编队构型的长期变化。因此，为了长时间的维持编队飞行的构型，最好的方法是对轨道平根数进行控制。

利用星载GPS接收机能够得到A星绝对轨道GPS数据，即A星相对地心的位置速度(WGS84坐标系)。还能够得到差分GPS数据，即两颗星之间的相对位置速度。如图1所示，将A星绝对轨道GPS数据通过坐标系转换得到惯性系下的位置速度，再转换为A星瞬时轨道根数和A星平根数。将A星绝对位置速度与相对位置速度相“加”，得到B星的绝对位置速度。然后在转换为B星瞬时轨道根数和平根数。将A星和B星的平根数进行处理得到两星的平根数差。两个卫星的平根数差的变化非常缓慢，短时间内近似为常值，对得到的平根数差进行一阶滤波可以进一步减小随机噪声的影响。这种方法得到B星平根数可能和B星真实的平根数有一定的误差，但是最后得到的两星平根数差主要取决于相对测量的精度。由于两星相对位置速度较为精确，这种方法得到的平根数差也就较为精确。

平根数差包括平半长轴差、平偏心率矢量差与平倾角矢量差，其中平半长轴差与平偏心率矢量差反映了两星轨道平面内的相对运动形状，平倾角矢量差反映了轨道平面外的相对运动形状。根据相对运动规律，将编队构型的控制分为两部分：轨道平面内部分与轨道平面外部分。

轨道平面外部分控制策略为：对平倾角矢量差与期望的倾角矢量进行比较，若偏差超过设计的阈值，则进行一次平面外控制。

轨道平面内的控制策略为：使用间隔半个轨道周期的两次变轨进行平半长轴差和平偏心率矢量差的联合控制。平偏心率矢量差的控制目标为零。平半长轴差的控制目标根据沿飞行方向的相对距离的大小分成5个区间，如图3所示，编队飞行要求两星在飞行方向(X轴)上的相对位置保持在X1～X4范围内。当编队飞行x方向相对位置大于X4或小于X1时，设置较大的半长轴差控制目标Δa₁，使其较快回到编队飞行区域内。在X1～X2和X3～X4设置两个减速区，当处在这两个区时，设置较小的半长轴差控制目标Δa₂，减小漂移速度。在X2～X3之间时不进行X-Z平面内控制，使相对运动自由漂移，减少轨控的燃料消耗。

如图1所示，本发明的编队飞行控制方法，具体包括如下步骤：

(1)计算A星瞬时轨道根数

星上差分GPS接收机能够进行绝对定位与相对定位，提供WGS84坐标系下的绝对位置速度和两星相对位置速度。这里只利用A星的绝对位置速度计算A星的绝对轨道根数。B星的绝对位置和速度不使用。

根据接收机输出的A星绝对位置和速度，以及WGS84坐标系和J2000.0惯性系的转换关系，转换到惯性系下，并将位置速度转换为瞬时轨道根数，瞬时轨道根数包括半长轴a_A,偏心率e_A，倾角i_A，升交点赤经Ω_A，近地点幅角ω_A，平近点角M_A，纬度幅角u_A以及λ_A＝ω_A+M_A。将位置、速度转换为瞬时轨道根数的方法是现有技术。

(2)计算A星轨道平根数

为了避免近圆轨道小偏心率引起的奇异，定义三个新的变量ξ,η,λ来代替e,ω,M，三个新变量的定义为ξ＝ecosω，η＝-esinω，λ＝ω+M。轨道瞬时根数为表示为σ＝[a_A,ξ_A,η_A,i_A,Ω_A,λ_A]^T，轨道平根数记为根据轨道平根数理论，瞬时根数与平根数的关系为其中为短周期项，短周期项的值与平根数有关，这是一个复杂的非线性方程，无法写出平根数的解析公式。因此，本发明用迭代方法计算平根数，如图2所示，具体步骤中的计算公式如下：

a.平根数赋初值

$\stackrel{\&OverBar;}{a}=a_A$

$\stackrel{\&OverBar;}{\mathrm{\&xi;}}=e_A\cos {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&eta;}}={-e}_A\sin {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{i}=i_A$

$\stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\mathrm{\&Omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\mathrm{\&lambda;}}_A$

为平半长轴中间量，和为平偏心率矢量中间量，为平倾角中间量，为平升交点赤经中间量，为λ_A的平根中间量；

b.短周期项计算，包括半长轴短周期项a_s，倾角短周期项i_s，升交点赤经短周期项Ω_s，偏心率矢量的短周期项ξ_s和η_s，幅角和λ_A的短周期项λ_s

$\stackrel{\&OverBar;}{p}=\stackrel{\&OverBar;}{a}(1-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}^2-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}^2),$ 是轨道半通径

其中R_e＝6378.137km为地球赤道半径，J₂＝0.001623945为地球非球形J2项引力系数，是中间变量，

$a_s=\frac{{\tilde{A}}_2}{\stackrel{\&OverBar;}{a}}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)$

$i_s=\frac{{\tilde{A}}_2}{4{\stackrel{\&OverBar;}{p}}^2}\sin \left(2\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)$

${\mathrm{\&Omega;}}_s=\frac{{\tilde{A}}_2}{2{\stackrel{\&OverBar;}{p}}^2}\cos \left(\stackrel{\&OverBar;}{i}\right)\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)$

${\mathrm{\&xi;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\cos \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[1-1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]+\frac{7}{12}{\sin }^2\left(i_0\right)\cos \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\}$

${\mathrm{\&eta;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\sin \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[-1+1.75{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]-\frac{7}{12}{\sin }^2\left(i_0\right)\sin \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\}$

${\mathrm{\&lambda;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}[-0.5+1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)$

c.平根数计算

${\stackrel{\&OverBar;}{a}}_n=a_A-a_s$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n=e_A\cos {\mathrm{\&omega;}}_A-{\mathrm{\&xi;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n=-e_A\sin {\mathrm{\&omega;}}_A-{\mathrm{\&eta;}}_s$

${\stackrel{\&OverBar;}{i}}_n=i_A-i_s$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n={\mathrm{\&Omega;}}_A-{\mathrm{\&Omega;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n={\mathrm{\&lambda;}}_A-{\mathrm{\&lambda;}}_s$

d.迭代收敛判断

如果( $|{\stackrel{\&OverBar;}{a}}_n-\stackrel{\&OverBar;}{a}|<{\mathrm{\&epsiv;}}_a$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&xi;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&xi;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&eta;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&eta;}}$

且 $|{\stackrel{\&OverBar;}{i}}_n-\stackrel{\&OverBar;}{i}|<{\mathrm{\&epsiv;}}_i$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&Omega;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&Omega;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&lambda;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}$ )

转e

否则

$\begin{array}{ccc}\stackrel{\&OverBar;}{a}={\stackrel{\&OverBar;}{a}}_n& \stackrel{\&OverBar;}{\mathrm{\&xi;}}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n& \stackrel{\&OverBar;}{\mathrm{\&eta;}}={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n\end{array}$

$\begin{array}{ccc}\stackrel{\&OverBar;}{i}={\stackrel{\&OverBar;}{i}}_n& \stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n& \stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n\end{array}$

转b

其中收敛条件的设置与计算机字长、计算精度要求相关，典型的取值如下：

ε_a＝10^-8，ε_ξ＝10^-7，ε_η＝10^-7，ε_i＝10^-4，ε_Ω＝ε_λ＝10^-4

e.输出平根数

${\stackrel{\&OverBar;}{a}}_A={\stackrel{\&OverBar;}{a}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n$

${\stackrel{\&OverBar;}{i}}_A={\stackrel{\&OverBar;}{i}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n$

(3)计算B星轨道平根数

首先根据接收机输出的B星相对于A星的相对位置和速度，以及WGS84坐标系和J2000.0惯性系的转换关系，转换到惯性系下，得到惯性下B星相对与A星的相对位置[x_r y_r z_r]^T、相对速度 ${\left[\begin{array}{ccc}{\stackrel{\&CenterDot;}{x}}_r& {\stackrel{\&CenterDot;}{y}}_r& {\stackrel{\&CenterDot;}{z}}_r\end{array}\right]}^T.$

然后根据A星的绝对位置[x_A y_A z_A]^T、绝对速度 ${\left[\begin{array}{ccc}{\stackrel{\&CenterDot;}{x}}_A& {\stackrel{\&CenterDot;}{y}}_A& {\stackrel{\&CenterDot;}{z}}_A\end{array}\right]}^T.$ 以及B星相对与A星的相对位置[x_r y_r z_r]^T、相对速度 ${\left[\begin{array}{ccc}{\stackrel{\&CenterDot;}{x}}_r& {\stackrel{\&CenterDot;}{y}}_r& {\stackrel{\&CenterDot;}{z}}_r\end{array}\right]}^T.$ 计算B星的绝对位置、速度，计算公式如下：

$\left[\begin{array}{c}{x}_B\\ y_B\\ z_B\end{array}\right]=\left[\begin{array}{c}{x}_A\\ y_A\\ z_A\end{array}\right]+A_{\mathrm{oi}}^T\left[\begin{array}{c}{x}_r\\ y_r\\ z_r\end{array}\right]$

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_B\\ {\stackrel{\&CenterDot;}{y}}_B\\ {\stackrel{\&CenterDot;}{z}}_B\end{array}\right]=\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_A\\ {\stackrel{\&CenterDot;}{y}}_A\\ {\stackrel{\&CenterDot;}{z}}_A\end{array}\right]=A_{\mathrm{oi}}^T(\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_r\\ {\stackrel{\&CenterDot;}{y}}_r\\ {\stackrel{\&CenterDot;}{z}}_r\end{array}\right]+\left[\begin{array}{c}-{\mathrm{\&omega;}}_o{z}_r\\ 0\\ {\mathrm{\&omega;}}_0{x}_r\end{array}\right])$

其中：A_oi为惯性系到A星轨道坐标系(Z轴指向地心，Y轴指向轨道面负法向，X轴与Y轴和Z轴垂直且指向飞行方向)的转换矩阵，ω_o为A星的瞬时轨道角速度地球引力常数μ＝398600.5km³/s²。

然后，由B星绝对位置[x_B y_B z_B]^T、绝对速度 ${\left[\begin{array}{ccc}{\stackrel{\&CenterDot;}{x}}_B& {\stackrel{\&CenterDot;}{y}}_B& {\stackrel{\&CenterDot;}{z}}_B\end{array}\right]}^T$ 计算B星瞬时轨道根数，包括半长轴a_B、偏心率e_B、轨道倾角i_B、升交点赤经Ω_B、近地点幅角ω_B、平近点角M_B、真近点角f_B和纬度幅角u_B。

最后，由B星瞬根数计算B星平根数,算法与步骤二相同，只是将公式中a_A,e_A,i_A,Ω_A,ω_A,λ_A的下标A改为B即可，将计算结果记为并计算平均轨道角速度

(4)两星平根数差计算

两星平根数差包括平半长轴差偏心率矢量差和倾角矢量差和两星平根数差的计算如下：

$\mathrm{\&Delta;}\tilde{a}={\stackrel{\&OverBar;}{a}}_A-{\stackrel{\&OverBar;}{a}}_B$

$\mathrm{\&Delta;}{\tilde{e}}_x={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_B$

$\mathrm{\&Delta;}{\tilde{e}}_y={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_B-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_A$

$\mathrm{\&Delta;}{\tilde{i}}_x=\sin \left({\stackrel{\&OverBar;}{i}}_A\right)\sin ({\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_B)$

$\mathrm{\&Delta;}{\tilde{i}}_y=-\sin \left({\stackrel{\&OverBar;}{i}}_A\right)\cos \left({\stackrel{\&OverBar;}{i}}_B\right)\cos ({\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_B)+\sin \left({\stackrel{\&OverBar;}{i}}_B\right)\cos \left({\stackrel{\&OverBar;}{i}}_A\right)$

(5)平根数差滤波

对步骤(4)中得到的平根数进行一阶滤波以减小随机噪声，滤波计算公式如下：

$\mathrm{\&Delta;a}={\mathrm{\&Delta;a}}^{-}+0.01(\mathrm{\&Delta;}\tilde{a}-{\mathrm{\&Delta;a}}^{-})$

${\mathrm{\&Delta;e}}_x={\mathrm{\&Delta;e}}_x^{-}+0.01(\mathrm{\&Delta;}{\tilde{e}}_x-\mathrm{\&Delta;}{e}_x^{-})$

${\mathrm{\&Delta;e}}_y={\mathrm{\&Delta;e}}_y^{-}+0.01(\mathrm{\&Delta;}{\tilde{e}}_y-\mathrm{\&Delta;}{e}_y^{-})$

${\mathrm{\&Delta;i}}_x={\mathrm{\&Delta;i}}_x^{-}+0.01(\mathrm{\&Delta;}{\tilde{i}}_x-\mathrm{\&Delta;}{i}_x^{-})$

${\mathrm{\&Delta;i}}_y={\mathrm{\&Delta;i}}_y^{-}+0.01(\mathrm{\&Delta;}{\tilde{i}}_y-\mathrm{\&Delta;}{i}_y^{-})$

其中上标减号表示上一周期的值，等式左边均为滤波后的值。滤波状态量的初值均为零。公式中的0.001和0.01等数值为滤波系数。滤波系数可综合考虑噪声大小和滤波收敛速度进行选取。

(6)控制策略计算

由于编队飞行的两颗卫星的相对运动可以分为轨道平面内的相对运动以及垂直于轨道平面方向上的相对运动(轨道面外)两部分，如图4所示。这两部分相对运动是解耦的，可以分别进行控制。

a.轨道平面外控制策略计算

根据A星相对于B星在轨道面外的相对运动运动方程，在某一轨道纬度幅角处，施加轨道面法向的速度增量，会改变两星倾角矢量差。倾角矢量差改变量的大小与速度增量的大小和轨道纬度幅角有关。反之，若给定倾角矢量差改变量，则可以求出所需的速度增量以及施加控制所需的轨道纬度幅角。

本发明中，平面外控制的目标是将两星倾角矢量差控制为[0,0]^T，若当前相对倾角矢量差[Δi_x,Δi_y]^T，可得到速度增量Δv_y及轨控时的纬度幅角u_y有两组解：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_y={\stackrel{\&OverBar;}{a}}_B{\stackrel{\&OverBar;}{n}}_B\sqrt{{\mathrm{\&Delta;i}}_x^2+{\mathrm{\&Delta;i}}_y^2}\\ u_y=\arctan 2({\mathrm{\&Delta;i}}_x,-{\mathrm{\&Delta;i}}_y)\end{array}\right.$ 或 $\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_y=-{\stackrel{\&OverBar;}{a}}_B{\stackrel{\&OverBar;}{n}}_B\sqrt{{\mathrm{\&Delta;i}}_x^2+{\mathrm{\&Delta;i}}_y^2}\\ u_y=\mathrm{\&pi;}+\arctan 2({\mathrm{\&Delta;i}}_x,-{\mathrm{\&Delta;i}}_y)\end{array}\right.$

执行轨道平面外控制时，可以从上面两组解中按照轨道纬度幅角找一个距离当前纬度幅角最近处进行，以便尽早的施加控制。

根据以上所述的平面外控制原理，本发明具体做法如下：

首先根据测量误差、摄动影响和控制量大小等因素来设计倾角矢量偏差的阈值。编队飞行过程中，对倾角矢量差[Δi_x,Δi_y]^T与期望的倾角矢量差[0,0]^T进行比较若偏差Δi超过设计的阈值Δi_max(即Δi>Δi_max)，则进行一次平面外控制计算，得到开机时间t_y和速度增量Δv_y。具体计算过程如下。

${\mathrm{\&Delta;v}}_y={\stackrel{\&OverBar;}{a}}_B{\stackrel{\&OverBar;}{n}}_B\sqrt{{\mathrm{\&Delta;i}}_x^2+{\mathrm{\&Delta;i}}_y^2}$

u_y＝arctan2(Δi_x,-Δi_y)

Δu1＝u_y-uA

$\mathrm{\&Delta;}{\hat{u}}_1=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_1,2\mathrm{\&pi;}),$

${\mathrm{\&Delta;u}}_2=\mathrm{\&Delta;}{\hat{u}}_1+\mathrm{\&pi;},$

$\mathrm{\&Delta;}{\hat{u}}_2=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_2,2\mathrm{\&pi;}),$

如果 $\mathrm{\&Delta;}{\hat{u}}_1<\mathrm{\&Delta;}{\hat{u}}_2$ 则

$t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_1}{{\stackrel{\&OverBar;}{n}}_B},\mathrm{\&Delta;}{\hat{v}}_y={\mathrm{\&Delta;v}}_y;$

否则

$t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_2}{{\stackrel{\&OverBar;}{n}}_B}$

$\mathrm{\&Delta;}{\hat{v}}_y=-{\mathrm{\&Delta;v}}_y$

结束

其中arctan2(y,x)函数表示先计算arctan(y/x)，并根据x和y的正负取合适的象限，最后得到的角度值在[-π，π)区间内。mod_2PI(x,2π)表示将角度值x的取值归算到[0,2π)区间内。

b.轨道平面内控制策略计算

根据轨道平面内(X-Z平面内)的运动方程，在一定的轨道纬度幅角处时间的切向速度增量使得平半长轴差和偏心率矢量差均发生改变。使用相隔半个轨道周期(轨道纬度幅角相差π)的两个速度增量可以对平半长轴差和偏心率矢量差进行联合控制。即给定平半长轴差的改变量、偏心率矢量差的改变量，可以求解出第一个速度增量Δv_x1、第二个速度增量Δv_x2和分别施加两个速度增量的轨道纬度幅角u_x1和u_x2。

平半长轴差的控制目标根据沿飞行方向的相对距离的大小分成5个区间，如图3所示，编队飞行要求两星在飞行方向(X轴)上的相对位置保持在X1～X4范围内。当编队飞行x方向相对位置大于X4或小于X1时，设置较大的半长轴差控制目标Δa₁，使其较快回到编队飞行区域内。在X1～X2和X3～X4设置两个减速区，当处在这两个区时，设置较小的半长轴差控制目标Δa₂，减小漂移速度。在X2～X3之间时不进行X-Z平面内控制，使相对运动自由漂移，减少轨控的燃料消耗。

${\mathrm{\&Delta;a}}_{\mathrm{new}}=\left\{\begin{array}{cc}-{\mathrm{\&Delta;a}}_1& x_r<X1\\ -{\mathrm{\&Delta;a}}_2& X1\&le;x_r<X2\\ {\mathrm{\&Delta;a}}_2& X3<x_r\&le;X4\\ {\mathrm{\&Delta;a}}_1& x_r>X4\end{array}\right.$

其中x_r是两星在X轴方向(飞行方向)上的相对位置，5个区间之间的边界点满足X1<X2<X3<X4。

相对偏心率矢量的大小决定了相对运动轨迹椭圆的大小。本发明中，希望两星的相对位置在一个轨道周期内的变化尽量小，将相对偏心率矢量的控制目标为[0,0]^T。

按照上述轨道平半长轴差的控制策略，当需要控制时(x_r<X2或x_r>X3时)，进行一次平面内双脉冲控制计算,得到两个脉冲速度增量Δv_x1和Δv_x2，以及两个脉冲对应轨道纬度幅角u_x1和u_x2，计算公式如下：

$\mathrm{\&Delta;e}=\sqrt{{\mathrm{\&Delta;e}}_x^2+{\mathrm{\&Delta;e}}_y^2}$

dΔa＝Δa_new-Δa

${\mathrm{\&Delta;v}}_a=(1/4){\stackrel{\&OverBar;}{n}}_B\left(\mathrm{d\&Delta;a}\right)$

${\mathrm{\&Delta;v}}_e=(1/4){\stackrel{\&OverBar;}{n}}_B{\stackrel{\&OverBar;}{a}}_B\left(\mathrm{\&Delta;e}\right)$

Δv_x1＝Δv_a+Δv_e

Δv_x2＝Δv_a-Δv_e

u_x1＝arctan2(-Δe_y,-Δe_x)

${\hat{u}}_{x1}=\mathrm{mod}\_2\mathrm{PI}(u_{x1},2\mathrm{\&pi;})$

$u_{x2}={\hat{u}}_{x1}+\mathrm{\&pi;}$

${\hat{u}}_{x2}=\mathrm{mod}\_2\mathrm{PI}(u_{x2},2\mathrm{\&pi;})$

由于两个脉冲的执行不分先后顺序，为了尽早的实现控制，对两个轨道纬度幅角u_x1和u_x2进行比较，选择距离当前纬度幅角u_A较近的一个脉冲先执行。计算先执行的速度增量对应的轨道纬度幅角与当前轨道纬度幅角的差Δu。计算过程如下：

如果( ${\hat{u}}_{x1}<\mathrm{\&pi;}$ 且( $u_A>{\hat{u}}_{x2}$ 或 $u_A>{\hat{u}}_{x1}$ ))

或( ${\hat{u}}_{x1}\&GreaterEqual;\mathrm{\&pi;}$ 且 $u_A>{\hat{u}}_{x2}$ 且 $u_A\&le;{\hat{u}}_{x1}$ ))

令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;u}={\hat{u}}_{x1}-u_A,$

如果(( ${\hat{u}}_{x1}<\mathrm{\&pi;}$ 且 ${\hat{u}}_{x1}<u_A\&le;{\hat{u}}_{x2}$ ))

或( ${\hat{u}}_{x1}\&GreaterEqual;\mathrm{\&pi;}$ 且( $u_A>{\hat{u}}_{x1}$ 或 $u_A\&le;{\hat{u}}_{x2}$ )))

令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;u}={\hat{u}}_{x2}-u_A,$

经过上面判断后，第一个脉冲的速度增量总是为第二个脉冲的速度增量总是为最后计算脉冲执行的时刻，计算方法如下：

令获得最终的纬度幅角的差

根据公式计算第一脉冲的开机时间t_x1，根据公式计算第二脉冲的开机时间t_x2。

(7)速度增量脉冲执行

根据步骤(6)计算的相对倾角矢量控制所需的速度增量Δv_y、半长轴差和偏心率矢量差控制所需的速度增量Δv_x1和Δv_x2以及星上轨控推力器的推力大小、卫星质量分别计算每个速度增量所需的开机时间。

首先根据速度增量Δv计算开机时长Δt＝m×Δv/F，其中m为卫星质量，F为轨控推力器的推力大小。

然后，根据在轨控推力器开机时姿态控制方式的不同，分为两种发动机开机方式：

姿态控制使用动量轮。当开机时长较短(例如小于3秒)时使用这种方式。为了避免轨控推力器连续喷气所产生的干扰力矩影响姿态控制，轨控推力器使用小脉冲喷气控制，在每个控制周期内执行一个短时间的喷气，例如在250ms的控制周期内喷气50ms，一直到喷气时间累积到大于Δt时关闭轨控推力器。这种方式下可以避免姿态控制喷气对轨道产生干扰，影响轨道控制的精度。

姿态控制使用姿控推力器。当开机时长较长(例如大于3秒)时使用这种方式。这种方式下轨控推力器连续喷气，使用加速度计测量速度增量，进行速度增量的累积，当速度增量累积到大于Δv时关闭轨控推力器。这种方式下姿态控制喷气对轨道控制的影响相对较小。

对于轨道平面内控制，需要在计算的开机时间t_x1和t_x2分两次启动在X方向的轨控推力器，根据相应的速度增量和计算开机时长，并按照上述方式进行X方向的控制。

对于轨道平面外控制，需要在计算的开机时间t_y开启在Y方向的轨控推力器，根据计算的计算开机时长，并按照上述方式进行Y方向的控制。

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

Claims (2)

1.一种卫星自主编队飞行控制方法，其特征在于，包括如下步骤：

(1)根据A星绝对轨道GPS数据获得A星绝对位置、速度；根据差分GPS数据获得A星和B星的相对位置、速度；根据A星绝对位置、速度与所述相对位置、速度得到B星的绝对位置、速度；

(2)根据A星绝对位置、速度计算A星瞬时轨道根数，根据B星绝对位置、速度计算B星瞬时轨道根数；A星瞬时轨道根数包括半长轴a_A,偏心率e_A，倾角i_A，升交点赤经Ω_A，近地点幅角ω_A，平近点角M_A，纬度幅角u_A以及幅角和λ_A，λ_A＝ω_A+M_A；

(3)根据A星瞬时轨道根数计算A星平根数，根据B星瞬时轨道根数计算B星平根数；A星平根数与B星平根数的计算方法相同；A星平根数包括平半长轴平偏心率矢量和平倾角和平升交点赤经B星平根数包括平半长轴平偏心率矢量和平倾角和平升交点赤经

(4)将A星平根数和B星平根数进行处理得到两星的平根数差，两星的平根数差包括平半长轴差平偏心率矢量差和平倾角矢量差和 $\mathrm{\&Delta;}{\tilde{i}}_y;$ $\mathrm{\&Delta;}\tilde{a}={\stackrel{\&OverBar;}{a}}_A-{\stackrel{\&OverBar;}{a}}_B,\mathrm{\&Delta;}{\tilde{e}}_x={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_B,\mathrm{\&Delta;}{\tilde{e}}_y={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_B-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_A,\mathrm{\&Delta;}{\tilde{i}}_x=\sin \left({\stackrel{\&OverBar;}{i}}_A\right)\sin ({\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_B),$ $\mathrm{\&Delta;}{\tilde{i}}_y=-\sin \left({\stackrel{\&OverBar;}{i}}_A\right)\cos \left({\stackrel{\&OverBar;}{i}}_B\right)\cos ({\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A-{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_B)+\sin \left({\stackrel{\&OverBar;}{i}}_B\right)\cos \left({\stackrel{\&OverBar;}{i}}_A\right);$

(5)对平根数差进行一阶滤波获得滤波后的平根数差，所述滤波后的平根数差包括平半长轴差Δa、偏心率矢量差Δe_x和Δe_y、倾角矢量差Δi_x和Δi_y；

(6)判断是否需要进行轨道平面外控制，当大于阈值Δi_max时，转入步骤(7)进行轨道平面外控制计算；当不大于阈值Δi_max时，不进行轨道平面外控制，转入步骤(10)；

同时判断是否需要进行轨道平面内控制，当x_r<X2或x_r>X3时，转入步骤(8)进行轨道平面内控制计算；当X2≤x_r≤X3时，不进行轨道平面内控制，转入步骤(10)；x_r为A星和B星在飞行方向上的相对位置，X2，X3为位置阈值；

(7)轨道平面外控制计算的步骤如下：

(7.1)根据公式计算平面外控制所需要的速度增量Δv_y；为平均轨道角速度，μ为地球引力常数；

(7.2)根据公式u_y＝arctan2(Δi_x,-Δi_y)计算平面外控制所需要的纬度幅角u_y；

(7.3)计算开机时间t_y和最终的轨道平面外控制所需要的速度增量 $\mathrm{\&Delta;}{\hat{v}}_y;$

令Δu₁＝u_y-u_A，

$\mathrm{\&Delta;}{\hat{u}}_1=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_1,2\mathrm{\&pi;}),$

${\mathrm{\&Delta;u}}_2=\mathrm{\&Delta;}{\hat{u}}_1+\mathrm{\&pi;},$

$\mathrm{\&Delta;}{\hat{u}}_2=\mathrm{mod}\_2\mathrm{PI}({\mathrm{\&Delta;u}}_2,2\mathrm{\&pi;}),$

如果 $\mathrm{\&Delta;}{\hat{u}}_1<\mathrm{\&Delta;}{\hat{u}}_2,$ 则令 $t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_1}{{\stackrel{\&OverBar;}{n}}_B},$ t为当前时间，令 $\mathrm{\&Delta;}{\hat{v}}_y={\mathrm{\&Delta;v}}_y;$ 否则，令 $t_y=t+\frac{\mathrm{\&Delta;}{\hat{u}}_2}{{\stackrel{\&OverBar;}{n}}_B},\mathrm{\&Delta;}{\hat{v}}_y=-{\mathrm{\&Delta;v}}_y;$

(7.4)转入步骤(9)；

(8)轨道平面内控制计算的步骤如下：

(8.1)根据如下公式计算两个脉冲速度增量Δv_x1和Δv_x2，以及两个脉冲对应的轨道纬度幅角和

${\mathrm{\&Delta;v}}_a=(1/4){\stackrel{\&OverBar;}{n}}_B({\mathrm{\&Delta;a}}_{\mathrm{new}}-\mathrm{\&Delta;a}),{\mathrm{\&Delta;a}}_{\mathrm{new}}=\left\{\begin{array}{cc}-{\mathrm{\&Delta;a}}_1& x_r<X1\\ -{\mathrm{\&Delta;a}}_2& X1\&le;x_r<X2\\ {\mathrm{\&Delta;a}}_2& X3<x_r\&le;X4\\ {\mathrm{\&Delta;a}}_1& x_r>X4\end{array}\right.,$ Δa₁为第一半长轴差控制目标、Δa₂为第二半长轴差控制目标，Δa₁>Δa₂，X1、X4为位置阈值，X1<X2<X3<X4，

${\mathrm{\&Delta;v}}_e=(1/4){\stackrel{\&OverBar;}{n}}_B{\stackrel{\&OverBar;}{a}}_B\left(\sqrt{{\mathrm{\&Delta;e}}_x^2+{\mathrm{\&Delta;e}}_y^2}\right),$ Δv_a、Δv_e为中间变量；

Δv_x1＝Δv_a+Δv_e，

Δv_x2＝Δv_a-Δv_e，

u_x1＝arctan2(-Δe_y,-Δe_x)，

${\hat{u}}_{x1}=\mathrm{mod}\_2\mathrm{PI}(u_{x1},2\mathrm{\&pi;}),$

u_x1、u_x2为中间变量，

${\hat{u}}_{x2}=\mathrm{mod}\_2\mathrm{PI}(u_{x2},2\mathrm{\&pi;});$

(8.2)确定最终的第一脉冲的速度增量和第二脉冲的速度增量计算先执行的速度增量对应的轨道纬度幅角与当前轨道纬度幅角的差Δu；

如果( ${\hat{u}}_{x1}<\mathrm{\&pi;}$ 且( $u_A>{\hat{u}}_{x2}$ 或 $u_A\&le;{\hat{u}}_{x1}$ ))或( ${\hat{u}}_{x1}\&GreaterEqual;\mathrm{\&pi;}$ 且 $u_A>{\hat{u}}_{x2}$ 且 $u_A\&le;{\hat{u}}_{x1}$ )，则令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;u}={\hat{u}}_{x1}-u_A;$

否则，令 $\mathrm{\&Delta;}{\hat{v}}_{x1}=\mathrm{\&Delta;}{v}_{x2},\mathrm{\&Delta;}{\hat{v}}_{x2}=\mathrm{\&Delta;}{v}_{x1},\mathrm{\&Delta;u}={\hat{u}}_{x2}-u_A;$

(8.3)令获得最终的纬度幅角的差根据公式计算第一脉冲的开机时间t_x1，根据公式计算第二脉冲的开机时间t_x2；

(8.4)转入步骤(9)；

(9)速度增量脉冲执行，

根据每个方向计算的速度增量、轨控推力器的推力大小以及卫星质量确定相应方向的轨控推力器的开机时长；

根据开机时长确定轨道控制和姿态控制方式，当开机时长大于设定值时，姿态控制使用姿控推力器；在计算的开机时间，打开轨控推力器，轨控推力器连续喷气，使用加速度计测量速度增量，进行速度增量的累积，当速度增量累积到大于计算的速度增量时关闭轨控推力器；当开机时长小于等于设定值时，姿态控制使用动量轮，在计算的开机时间，打开轨控推力器，轨控推力器在开机后的每个控制周期内执行一个短时间的喷气，一直到喷气时间累积到大于开机时长时关闭轨控推力器；

(10)结束。

2.根据权利要求1所述的卫星自主编队飞行控制方法，其特征在于，根据A星瞬时轨道根数计算A星平根数方法包括如下步骤：

a.平根数赋初值

令， $\stackrel{\&OverBar;}{a}=a_A$

$\stackrel{\&OverBar;}{\mathrm{\&xi;}}=e_A\cos {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&eta;}}=-e_A\sin {\mathrm{\&omega;}}_A$

$\stackrel{\&OverBar;}{i}=i_A$

$\stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\mathrm{\&Omega;}}_A$

$\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\mathrm{\&lambda;}}_A,$

其中为平半长轴中间量，和为平偏心率矢量中间量，为平倾角中间量，为平升交点赤经中间量，为λ_A的平根数中间量；

b.短周期项计算，短周期项包括半长轴短周期项a_s，倾角短周期项i_s，升交点赤经短周期项Ω_s，偏心率矢量的短周期项ξ_s和η_s，幅角和的短周期项λ_s，

$\stackrel{\&OverBar;}{p}=\stackrel{\&OverBar;}{a}(1-{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}^2-{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}^2),$ 是轨道半通径

其中R_e＝6378.137km为地球赤道半径，J₂＝0.001623945为地球非球形J2项引力系数，是中间变量，

$a_s=\frac{{\tilde{A}}_2}{\stackrel{\&OverBar;}{a}}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

$i_s=\frac{{\tilde{A}}_2}{4{\stackrel{\&OverBar;}{p}}^2}\sin \left(2\stackrel{\&OverBar;}{i}\right)\cos \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

${\mathrm{\&Omega;}}_s=\frac{{\tilde{A}}_2}{2{\stackrel{\&OverBar;}{p}}^2}\cos \left(\stackrel{\&OverBar;}{i}\right)\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right),$

${\mathrm{\&xi;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\cos \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[1-1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]+\frac{7}{12}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\cos \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\},$

${\mathrm{\&eta;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}\left\{\sin \left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\right[-1+1.75{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]-\frac{7}{12}{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)\sin \left(3\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\},$

${\mathrm{\&lambda;}}_s=\frac{{\tilde{A}}_2}{{\stackrel{\&OverBar;}{p}}^2}[-0.5+1.25{\sin }^2\left(\stackrel{\&OverBar;}{i}\right)]\sin \left(2\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right);$

c.平根数计算，得到一步迭代更新后的平半长轴偏心率矢量和平倾角平升交点赤经以及变量

${\stackrel{\&OverBar;}{a}}_n=a_A-a_s$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n=e_A\cos {\mathrm{\&omega;}}_A-{\mathrm{\&xi;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n=-e_A\sin {\mathrm{\&omega;}}_A-{\mathrm{\&eta;}}_s$

${\stackrel{\&OverBar;}{i}}_n=i_A-i_s$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n={\mathrm{\&Omega;}}_A-{\mathrm{\&Omega;}}_s$

${\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n={\mathrm{\&lambda;}}_A-{\mathrm{\&lambda;}}_s$

d.迭代收敛判断

如果( $|{\stackrel{\&OverBar;}{a}}_n-\stackrel{\&OverBar;}{a}|<{\mathrm{\&epsiv;}}_a$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&xi;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&xi;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&eta;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&eta;}}$ 且 $|{\stackrel{\&OverBar;}{i}}_n-\stackrel{\&OverBar;}{i}|<{\mathrm{\&epsiv;}}_i$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&Omega;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&Omega;}}$ 且 $|{\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n-\stackrel{\&OverBar;}{\mathrm{\&lambda;}}|<{\mathrm{\&epsiv;}}_{\mathrm{\&lambda;}}$ )，ε_a、ε_ξ、ε_η、ε_i、ε_Ω、ε_λ为阈值；转步骤e；

否则，令 $\stackrel{\&OverBar;}{a}={\stackrel{\&OverBar;}{a}}_n,\stackrel{\&OverBar;}{\mathrm{\&xi;}}={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n,\stackrel{\&OverBar;}{\mathrm{\&eta;}}={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n,\stackrel{\&OverBar;}{i}={\stackrel{\&OverBar;}{i}}_n,\stackrel{\&OverBar;}{\mathrm{\&Omega;}}={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n,\stackrel{\&OverBar;}{\mathrm{\&lambda;}}={\stackrel{\&OverBar;}{\mathrm{\&lambda;}}}_n,$ 转步骤b；

e.输出平根数，

令 ${\stackrel{\&OverBar;}{a}}_A={\stackrel{\&OverBar;}{a}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&xi;}}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&eta;}}}_n$

${\stackrel{\&OverBar;}{i}}_A={\stackrel{\&OverBar;}{i}}_n$

${\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_A={\stackrel{\&OverBar;}{\mathrm{\&Omega;}}}_n.$