CN102819266B - 一种拟周期j2不变相对轨道编队飞行控制方法 - Google Patents

Abstract

一种拟周期J2不变相对轨道编队飞行控制方法，以编队飞行卫星中飞行在标称轨道上的卫星为基准星，对编队飞行卫星中未飞行在标称轨道上的卫星施加控制力矩Tc进行调整。在给定相对位置初值后，对于每个时刻，计算具有Hamiltonian结构非受控J2型C-W方程的双曲特征值及稳定和不稳定流形，选择适当的控制增益，最终计算得到反馈控制力矩Tc。本发明方法基于更精确的模型得到的控制律，且只需要两颗卫星之间的相对位置反馈信息，对相对位置初值没有约束，计算量小，燃料消耗小，易于工程实现。

Description

一种拟周期J2不变相对轨道编队飞行控制方法

技术领域

本发明涉及一种编队飞行航天器的控制方法。

背景技术

编队飞行时，受限于星间链路的有效范围，编队卫星之间的相对位置需要保持在一定范围之内。同时，位置保持也保证了编队卫星群的相对构型，避免卫星之间的碰撞。此外，如干涉测量等特殊探测任务，还需要保持卫星间的相对方位，这对编队卫星的位置保持提出了较高的要求。

由C-W方程设计的相对周期运动，在地球J2摄动的影响下将不再闭合，这主要是由于主星和从星的子午面转动速率不同导致。(当编队卫星多于2颗时，称运行于既定轨道的为主星，除主星外的其他编队卫星为从星。)现有编队飞行卫星的维持控制主要存在以下几点问题：

1)大都是基于经典C-W方程设计；而C-W方程没有考虑J2摄动的影响，故控制器在消除J2项引起迹向漂移时，将耗费较多的燃料。

2)受C-W方程的影响，长期以来人们都将编队轨迹限定在周期轨道；实际上，编队飞行的有界运动，也可以拟周期轨道的形式实现。

3)经典C-W方程对编队卫星的初始状态有一定约束，编队构型难以快速改变；然而，在许多空间任务中卫星群快速改变相对构型能力非常重要。

4)现有编队飞行控制方法，计算量较大、不易于在工程上实现。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种拟周期J2不变相对轨道编队卫星的控制方法，燃料消耗少，以拟周期轨道运行，可以快速改变相对构型，且易于工程实现。

本发明的技术解决方案是：一种拟周期J2不变相对轨道编队飞行控制方法，以编队飞行卫星中飞行在标称轨道上的卫星为基准星，对编队飞行卫星中未飞行在标称轨道上的卫星施加控制力矩Tc进行调整，T_c＝B^-Δr-σ²G[u₊u₊ ^H+u_-u_- ^H]Δr，其中B^-＝(B-B^T)/2，B⁺＝(B+B^T)/2，Δr为基准星轨道坐标系下待调整卫星与基准星的相对位置矢量，B_λ＝(b₁₁+b₂₂+a₁₁a₂₂+a₁₂ ²)，C_λ＝(b₁₁b₂₂-b₁₂ ²)，a_ij表示矩阵A的第i行第j个元素，b_ij表示矩阵B⁺的第i行第j个元素，u₊、u_-分别满足[σ²I-σA+B]u₊＝0，[σ²I+σA+B]u_-＝0，G为控制增益，

$A=2\{\left[\begin{array}{ccc}0& -\cos \left(i\right)& \cos \left(u\right)\·\sin \left(i\right)\\ \cos \left(i\right)& 0& -\sin \left(i\right)\·\sin \left(u\right)\\ -\cos \left(u\right)\·\sin \left(i\right)& \sin \left(i\right)\·\sin \left(u\right)& 0\end{array}\right]\stackrel{\·}{\mathrm{\Ω}}+\left[\begin{array}{ccc}0& 0& -\sin \left(u\right)\\ 0& 0& -\cos \left(u\right)\\ \sin \left(u\right)& \cos \left(u\right)& 0\end{array}\right]\stackrel{\·}{i}+\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·}{u}\}$

$B=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·\·}{u}+\left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 0\end{array}\right]{\stackrel{\·}{u}}^2+\left[\begin{array}{ccc}\frac{1}{2}(-3+\cos \left(2u\right))& -\cos \left(u\right)\·\sin \left(u\right)& 0\\ -\cos \left(u\right)\·\sin \left(u\right)& \frac{1}{2}(-3-\cos \left(2u\right))& 0\\ 2\·\cos \left(u\right)& -2\·\sin \left(u\right)& -1\end{array}\right]\stackrel{\·}{u}\·\stackrel{\·}{i}$

$+\left[\begin{array}{ccc}\frac{1}{2}(1+\cos \left(i\right))(-3-\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos \left(2u\right)& \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& \cos \left(i\right)\·\sin \left(i\right)\·\sin \left(u\right)\\ \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& -\frac{1}{2}(1+\cos \left(i\right))(-3+\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos \left(2u\right))& \cos \left(i\right)\·\sin \left(i\right)\·\cos \left(u\right)\\ (2+\cos \left(i\right))\·\sin \left(i\right)\·\sin \left(u\right)& (2+\cos \left(i\right))\·\sin \left(i\right)\·\cos \left(u\right)& {-\sin }^2\left(i\right)\end{array}\right]\stackrel{\·}{u}\·\stackrel{\·}{\mathrm{\Ω}}$

$+\frac{\mathrm{\μ}}{r^3}\left[\begin{array}{ccc}-2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]+\frac{\mathrm{\μ}{J}_2{\mathrm{Re}}^2}{r^5}\left[\begin{array}{ccc}6(3\·{\sin }^2\left(i\right)\·{\sin }^2\left(u\right)-1)& -6\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& -6\·\sin \left(2i\right)\·\sin \left(u\right)\\ -\frac{9}{2}\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& 3\·{\sin }^2\left(i\right)\·\cos \left(2u\right)& \frac{3}{2}\·\sin \left(2i\right)\·\cos \left(u\right)\\ -\frac{9}{2}\·\sin \left(2i\right)\·\sin \left(u\right)& 3\·\sin \left(2i\right)\·\cos \left(u\right)& 3\·\cos \left(2i\right)\end{array}\right]$

u为主星纬度幅角，i为主星轨道倾角，Ω为主星升交点赤经，μ为地球的引力势能常数，Re为地球半径，J₂＝-108264×10^-3为常数，r为主星的矢径的模，I为两阶单位阵，角标T为矩阵的转置，字母上方的一点和两点分别表示一阶微分和二阶微分。

所述的控制增益G的选取方法为：令G从0开始逐渐增大，直至采用T_c＝B^-Δr-σ²G[u₊u₊ ^H+u_-u_- ^H]Δr描述的系统的所有Floquet乘子均分布在复域的单位圆上，则此时对应的G可以作为取值范围的下限，选取G时应大于此值。

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

1)本发明方法基于较C-W方程更精确的模型得到的控制律，在消除J2摄动影响时燃料消耗量相对于C-W方程得到的控制率要小；

2)采用本发明方法得到的相对轨道为拟周期的，便于实现更灵活的编队飞行卫星相对构型；

3)本发明方法在对卫星进行控制时，对相对位置初值没有约束，且编队构型可以快速改变；

4)本发明控制方法计算量小，且只需要相对位反馈信息，工程上易于实现。

附图说明

图1为本发明方法的控制原理图；

图2为J₂型C-W方程与精确J₂相对动力学方程所得到的运动轨迹的误差的时间历程；

图3(a)为受控系统的Floquet乘子的模随控制增益G变化的曲线，图3(b)为复平面上Floquet乘子分布情况；

图4为荡模式频率的时间历程；

图5为拟周期J₂不变轨道相对位置的时间历程；

图6控制加速度的时间历程。

具体实施方式

假设编队卫星中的一颗卫星无偏差地飞行在其预定地球轨道或者日地系统的Halo轨道上，为简便起见，称该星为主星；本发明拟周期J₂不变相对轨道编队飞行控制方法，仅对编队中除此星之外的其他卫星分别进行控制，称这些卫星为从星。对每一颗从星分别控制，控制方式相同，如图1所示，通过下述步骤来完成：

1、任意给定主星轨道坐标系下从星和主星之间的相对位置矢量初值Δr₀。

定义主星(Leader)轨道坐标系S_L(x y z)：原点取为主星质心，x轴为矢径(地心至质心)方向，z轴为动量矩方向，y轴由右手法则确定。在主星轨道坐标系下，从星与主星的相对位置矢量可表示Δr＝[x y z]^T。在实际编队任务中，相对位置矢量在主星轨道坐标系下的投影可以由相对导航系统直接得到。

因而，这里给定在主星轨道坐标系下，主星和从星之间的相对位置矢量Δr的初值为Δr₀＝[x₀ y₀ z₀]^T。

2、计算具有Hamiltonian结构非受控J₂型C-W方程的双曲特征值±σ及稳定和不稳定流形u±。

首先对精确J₂相对动力学方程进行一阶展开，得到线性时变的J2型C-W方程，然后对其系数矩阵进行修改以得到具有Hamiltonian结构的非受控J₂型C-W方程，随后计算该非受控J2型C-W方程的双曲特征值±σ及稳定和不稳定流形u_±。具体过程如下：

定义惯性坐标系：原点取为地心、x_i轴由地心指向春分点、z_i轴为地球旋转的动量矩方向、y_i轴由右手法则确定。显然，x_i轴和y_i轴张成地球赤道平面。在惯性坐标系下，从星和主星的相对位置矢量为Δr_I＝r_Follower-r_Leader，其中r_Follower和r_Leader分别为从星和主星在惯性坐标系下的位置矢量；

惯性坐标系到主星轨道坐标系的转换关系为：

R(Ω，i，u)＝R_z(u)R_x(i)R_z(Ω)        (1)

其中，u为主星纬度幅角，i为主星轨道倾角，Ω为主星升交点赤经，表示绕*轴转动角的转换矩阵。

转换矩阵(1)的逆变换为

F(Ω，i，u)＝R^-1＝R_z(-Ω)R_x(-i)R_z(-u)   (2)

则相对位置矢量在两坐标系的转换为

Δr_I＝R^-1·Δr＝F·Δr    (3)

对式(3)进行2次微分运算，可以得到

$\mathrm{\Δ}{\stackrel{\·\·}{r}}_I=F\·\mathrm{\Δ}\stackrel{\·\·}{r}+2\stackrel{\·}{F}\·\mathrm{\Δ}\stackrel{\·}{r}+\stackrel{\·\·}{F}\·\mathrm{\Δr}---\left(4\right)$

其中，为惯性坐标系下相对位置矢量的二阶微分；和分别表示主星轨道坐标系下相对位置矢量的二阶微分和一阶微分；和表示F的一阶微分和二阶微分。

考虑J2摄动时，从星(Follower)的引力势能为

$U_F=-\frac{\mathrm{\μ}}{r_F}+\frac{\mathrm{\μ}{J}_2{\mathrm{Re}}^2}{{2r}_F^3}(3\frac{z_F^2}{r_F^2}-1)---\left(5\right)$

式中，μ为地球的引力势能常数，Re为地球半径，J₂＝-1.08264×10^-3为常数，从星在惯性坐标系中的位置矢量的模以及其位置矢量在惯性坐标系中z_i轴上的分量分别为和z_F＝(r+x)sinusini+ycosusini+zcosi，其中r为主星的矢径的模。

上述引力势能也适用于主星，即

U_L＝U_F|_{(x，y，z)＝(0，0，0)}         (6)

因为当(x y z)＝(0 0 0)时，即为主星轨道坐标系的原点，根据主星轨道坐标系的定义，原点取在主星质心，因此该点的引力势能即为主星的引力势能。

设主星轨道径向、迹向、法向的单位矢量在惯性坐标系中的分量列阵分别为e_r、r_u、e_i，并记由于从主星轨道坐标系到惯性坐标系的转换矩阵可由轨道坐标系的单位矢量在惯性坐标系上的分量列阵表示，因而则在球坐标系下的梯度算子为

其中(与该球坐标系相对应的直角坐标系为地心节点坐标系S_n，即原点在地球中心，轴x_n由地心指向升交点；轴z_n垂直于轨道平面，与动量矩同方向；轴y_n按右手法则决定)。

则主星的J₂绝对动力学方程可表示为

${\stackrel{\·\·}{r}}_L=-{\▿}_L{U}_L---\left(8\right)$

其中为惯性坐标系下主星位置矢量的二阶微分。

从星轨道坐标系定义类似主星轨道坐标系：原点取为从星质心，x_F轴为从星矢径方向，z_F轴为动量矩方向，y_F轴由右手法则确定。设从星轨道径向、迹向、法向的单位矢量在惯性坐标系中的分量列阵分别为e_rF、e_uF、e_iF，并记类似地，以表示惯性坐标系下从星位置矢量的二阶微分，则从星的J₂绝对动力学方程为

${\stackrel{\·\·}{r}}_F=-{\▿}_F{U}_F---\left(9\right)$

其中 为从星的轨道坐标系矢量。

将式(9)和(8)相减，则得到J₂相对动力学方程

$\mathrm{\Δ}{\stackrel{\·\·}{r}}_I=-{\▿}_F{U}_F+{\▿}_L{U}_L---\left(10\right)$

显然，与具有不同的基底，不能同时参与运算。

因此，为在同一个基底下展开J2相对动力学方程(10)，定义线性算子 ${\▿}_{\mathrm{\Δr}}=\left(\begin{array}{ccc}\frac{\∂}{\∂x}& \frac{\∂}{\∂y}& \frac{\∂}{\∂z}\end{array}\right),$ 可得

其中，为在Δr＝0处的梯度值。联立式(4)和(11)并消去可得

$\mathrm{\Δ}\stackrel{\·\·}{r}+2F^{-1}\stackrel{\·}{F}\mathrm{\Δ}\stackrel{\·}{r}+F^{-1}\stackrel{\·\·}{F}\mathrm{\Δr}=-(\▿\·{\▿}_{\mathrm{\Δr}})U_F{|}_{\mathrm{\Δr}=0}\·\mathrm{\Δr}---\left(12\right)$

交换微分线性算子和的作用顺序，可得

$\mathrm{\Δ}\stackrel{\·\·}{r}+2F^{-1}\stackrel{\·}{F}\mathrm{\Δ}\stackrel{\·}{r}+F^{-1}\stackrel{\·\·}{F}\mathrm{\Δr}=-({\▿}_{\mathrm{\Δr}}\·\▿)U_F{|}_{\mathrm{\Δr}=0}\·\mathrm{\Δr}---\left(13\right)$

上式即为J₂型C-W方程的一般形式。

近一步展开式(13)，可以得到线性时变的J2型C-W方程的矩阵形式，即

$\mathrm{\Δ}\stackrel{\·\·}{r}+\mathrm{A\Δ}\stackrel{\·}{r}+\mathrm{B\Δr}=0---\left(14\right)$

其中，

$A=F^{-1}\stackrel{\·}{F}=2\{\left[\begin{array}{ccc}0& -\cos \left(i\right)& \cos \left(u\right)\·\sin \left(i\right)\\ \cos \left(i\right)& 0& -\sin \left(i\right)\·\sin \left(u\right)\\ -\cos \left(u\right)\·\sin \left(i\right)& \sin \left(i\right)\·\sin \left(u\right)& 0\end{array}\right]\stackrel{\·}{\mathrm{\Ω}}+\left[\begin{array}{ccc}0& 0& -\sin \left(u\right)\\ 0& 0& -\cos \left(u\right)\\ \sin \left(u\right)& \cos \left(u\right)& 0\end{array}\right]\stackrel{\·}{i}+\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·}{u}\}$

$B=F^{-1}\stackrel{\·\·}{F}+({\▿}_{\mathrm{\Δr}}\·\▿)U_F{|}_{\mathrm{\Δr}=0}$

$=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·\·}{u}+\left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 0\end{array}\right]{\stackrel{\·}{u}}^2+\left[\begin{array}{ccc}\frac{1}{2}(-3+\cos \left(2u\right))& -\cos \left(u\right)\·\sin \left(u\right)& 0\\ -\cos \left(u\right)\·\sin \left(u\right)& \frac{1}{2}(-3-\cos \left(2u\right))& 0\\ 2\·\cos \left(u\right)& -2\·\sin \left(u\right)& -1\end{array}\right]\stackrel{\·}{u}\·\stackrel{\·}{i}$

$+\left[\begin{array}{ccc}\frac{1}{2}(1+\cos \left(i\right))(-3-\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos \left({2u}_{-}\right)& \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& \cos \left(i\right)\·\sin \left(i\right)\·\sin \left(u\right)\\ \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& -\frac{1}{2}(1+\cos \left(i\right))(-3+\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos \left(2u\right))& \cos \left(i\right)\·\sin \left(i\right)\·\cos \left(u\right)\\ (2+\cos \left(i\right))\·\sin \left(i\right)\·\sin \left(u\right)& (2+\cos \left(i\right))\·\sin \left(i\right)\·\cos \left(u\right)& {-\sin }^2\left(i\right)\end{array}\right]\stackrel{\·}{u}\·\stackrel{\·}{\mathrm{\Ω}}$

$+\frac{\mathrm{\μ}}{r^3}\left[\begin{array}{ccc}-2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]+\frac{\mathrm{\μ}{J}_2{\mathrm{Re}}^2}{r^5}\left[\begin{array}{ccc}6(3\·{\sin }^2\left(i\right)\·{\sin }^2\left(u\right)-1)& -6\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& -6\·\sin \left(2i\right)\·\sin \left(u\right)\\ -\frac{9}{2}\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& 3\·{\sin }^2\left(i\right)\·\cos \left(2u\right)& \frac{3}{2}\·\sin \left(2i\right)\·\cos \left(u\right)\\ -\frac{9}{2}\·\sin \left(2i\right)\·\sin \left(u\right)& 3\·\sin \left(2i\right)\·\cos \left(u\right)& 3\·\cos \left(2i\right)\end{array}\right]$

其中角标T为矩阵的转置，I为两阶单位阵。

J₂型C-W方程的矩阵形式(14)包含经典C-W方程：令J₂＝0、ü＝0、和则J₂型C-W方程完全退化为经典C-W方程。其中， 分别为主星升交点赤经、轨道倾角和纬度幅角的一阶微分；ü为主星纬度幅角的二阶微分；为天龙角速度， 为一个轨道周期内平均轨道半长轴， 为一个轨道周期内平均轨道倾角。

下面以表1的主星轨道根数为算例，验证所得到J₂型C-W方程的有效性。相对位置初始为x₀＝y₀＝z₀＝1km，相对初始速度由经典C-W方程得到，即其它在系数矩阵A和B中出现的主星瞬时轨道根数u、i、Ω及其导数需要由外部输入或由近似理论给出。这里采用圆轨道根数的2阶近似解析解计算，参见文献Born G.H.，Motion of a Satellite under theInfluence of an Oblate Earth[R].Technical Report ASEN 3200，Univrsity ofCorado，Boulder，2001。

表1主星轨道根数

图2所示为J₂型C-W方程与精确J₂相对动力学方程所得到的运动轨迹的误差的时间历程。对于相对距离为1km量级的编队，其在6个小时内的误差小于50m；误差大小与Schweighart所得到J₂线性模型的误差一致。这表明所推导J₂型C-W方程的有效性；同时也说明了J₂相对动力学方程1阶近似的精度，即为该误差精度。J₂型C-W方程与精确J₂相对动力学方程得到的运动轨迹基本吻合。

由以上分析可知，J₂型C-W方程具有极高的精度；故以该方程作为设计模型得到的控制器，可完全适用于精确J₂相对模型甚至全摄动模型。

在J₂型C-W方程(14)右端加入控制力矩T_c，则可得到受控J₂型C-W模型方程，即控制器设计模型：

$\mathrm{\Δ}\stackrel{\·\·}{r}+\mathrm{A\Δ}\stackrel{\·}{r}+\mathrm{B\Δr}=T_c---\left(15\right)$

其中，A^T＝-A，B不具有对称性，因此该线性系统不具有Hamiltonian结构。为使方程具有Hamiltonian形式，对矩阵B作恰当的分解

$\mathrm{\Δ}\stackrel{\·\·}{r}+\mathrm{A\Δ}\stackrel{\·}{r}+B^{+}\mathrm{\Δr}=T_c-B^{-}\mathrm{\Δr}---\left(16\right)$

其中，B⁺＝(B+B^T)/2，B^-＝(B-B^T)/2，则B⁺和B^-分别为对称和斜对称矩阵。令T_c-B^-Δr＝0则得到修改后的非受控系统为

$\mathrm{\Δ}\stackrel{\·\·}{r}+\mathrm{A\Δ}\stackrel{\·}{r}+B^{+}\mathrm{\Δr}=0---\left(17\right)$

显然，修改后的非受控系统具有Hamiltonian结构，且原点为任意时刻的平衡点。

对于时刻t ∈[0，T_d]，T_d为天龙周期，式(17)

所示系统的原点具有双曲特征值_±σ(t)，称为双曲型平衡点；双曲型平衡点存在稳定流形u₊和不稳定流形u_-(u_±(t)实际上是式(17)所示系统通过双曲平衡点的解曲线，u₊(t)随时间增大接近原点，u_-(t)随时间增大远离原点)，且满足

其中I为2阶单位阵。考虑式(17)所示非受控系统，其特征方程为

det(Iλ²+Aλ+B⁺)＝λ⁴+(b₁₁+b₂₂+a₁₁a₂₂+a₁₂ ²)λ²+(b₁₁b₂₂-b₁₂ ²)＝0   (19)

式中，λ为特征方程的特征根，a_ij表示矩阵A的第i行第j个元素，b_ij表示矩阵B⁺的第i行第j个元素，记B_λ＝(b₁₁+b₂₂+a₁₁a₂₂+a₁₂ ²)，C_λ＝(b₁₁b₂₂-b₁₂ ²)，则可计算其双曲特征值

${\mathrm{\σ}}^2=\frac{-B_{\mathrm{\λ}}+\sqrt{{B_{\mathrm{\λ}}}^2-{4C}_{\mathrm{\λ}}}}{2}---\left(20\right)$

于是可以计算得到满足式(18)的特征向量u₊和u_-，即：

$u_{\±}=\frac{1}{\sqrt{1+{{\mathrm{\β}}_{\±}}^2}}\left[\begin{array}{c}1\\ {\mathrm{\β}}_{\±}\end{array}\right]---\left(21\right)$

其中， ${\mathrm{\β}}_{+}=-\frac{{\mathrm{\σ}}^2+b_{11}}{b_{12}-\mathrm{\σ}{a}_{12}}=-\frac{b_{12}+\mathrm{\σ}{a}_{12}}{{\mathrm{\σ}}^2+b_{22}},$ ${\mathrm{\β}}_{-}=-\frac{{\mathrm{\σ}}^2+b_{11}}{b_{12}+\mathrm{\σ}{a}_{12}}=-\frac{b_{12}-\mathrm{\σ}{a}_{12}}{{\mathrm{\σ}}^2+b_{22}}.$

3、根据双曲特征值±σ(t)及稳定和不稳定流形u_±(t)，选择适当的控制增益G，计算反馈控制力矩T_c。

根据式(18)，系统存在稳定和不稳定流形u_±，它们的存在会使得系统不稳定，因此，为了满足稳定性要求，需要在控制力矩T_c添加稳定和不稳定流形，从而将这两部分的作用抵消。同时需要注意，只反馈不稳定流形不足以使系统稳定。

因此，保Hamiltonian结构控制器(即反馈控制力矩T_c)可构造为

T_c＝B^-Δr-σ(t)²G[u₊u₊ ^H+u_-u_- ^H]Δr        (22)

其中，上标“H”表示Hermite转置，G为控制增益；σ、u₊和u_-分别为当前时刻的双曲特征值及稳定和不稳定流形。需要说明的是：若某时刻t，原点为椭圆型平衡点(即不具有双曲特征值)，则控制器中参数σ＝0。一般情况下，原点的拓扑性质交替变化，这是由J₂相对动力学迹向漂移近线性增长的属性决定的。

式(22)描述的控制器可以将平衡点的拓扑性质，即控制器通过在虚轴上配置极点，将双曲型平衡点转变为椭圆型平衡点。对于一般的系统(非Hamiltonian系统)，极点是不允许配置在虚轴上，因为被忽略的高阶非线性项可能致使系统不稳定。但对于Hamiltonian系统，Morse引理仍然可以保证受控非线性系统的Lyapunov稳定性。可以证明，由保Hamiltonian结构控制器生成的运动将是拟周期的且具有Lyapunov稳定。

具体地，将控制力矩T_c带入式(16)描述的受控Hamiltonian系统得到

$\mathrm{\Δ}\stackrel{\·\·}{r}+\mathrm{A\Δ}\stackrel{\·}{r}+\{B^{+}+{\mathrm{\σ}}^{2}G[u_{+}{u_{+}}^H+u_{-}{u_{-}}^H\left]\right\}\·\mathrm{\Δr}=0---\left(23\right)$

仍旧为周期系统。由线性周期Hamiltonian系统的稳定性定理可知，如果所有Floquet乘子均分布在复域的单位圆上(或Floquet乘子的模为1)，该Hamiltonian系统是Lyapunov稳定的。

满足稳定性要求的受控系统，其原点在任意时刻都为椭圆型平衡点，故当前时刻的运动可由简约的Hamiltonian函数描述

H＝H₀+ω₁τ₁+ω₂τ₂+ω₃τ₃+O(3)       (24)

其中，H₀为Hamiltonian函数在原点的值，τ_i，i＝1，2，3为作用量，ω_i，i＝1，2，3为三种振荡模式的频率，O(3)为高阶小量。三种振荡模式只是稳定Lissajous轨道的线性分解，且其x、y、z三轴的运动不可能完全解耦。数值仿真可以证明，通过调节式(22)所描述控制器的控制增益G，可以生成任意周期的拟周期相对运动。

因此，为了确定能够使得系统稳定的控制增益G，应该找出Floquet乘子随控制增益G的变化规律，从而确定使得所有Floquet乘子均分布在复域的单位圆上的G的取值范围，从而使得式(22)描述的控制器可以保证J₂拟周期不变轨道的长期稳定性。可以令G从0开始逐渐增大，直至系统的所有Floquet乘子均分布在复域的单位圆上(模为1)，则此时对应的G可以作为G取值范围的下限，选取G时应大于此值。

下面结合具体算例进一步说明：

根据

表1轨道根数建立的式(23)描述受控系统的Floquet乘子随控制增益G的变化规律如图3所示：当G＞0.27时，所有Floquet乘子均分布在复域的单位圆上，即式(22)描述的控制器可以保证J₂拟周期不变轨道的长期稳定性。图4给出了模式频率的时间历程：

4、对从星施加反馈控制力矩T_c，控制器可以由连续推进的高比冲小推力发动机实现。

式(22)描述的控制器根据当前相对位置进行实时反馈，增益选择越大，控制器输出也越大；编队半径越大，输出也越大。

若控制增益取G＝0.5；相对初始位置为x₀＝y₀＝z₀＝1km，则生成图5的稳定Lissajous轨道，所需要的控制加速度如图6所示：最大输出加速度在10^-5m/s²量级，5天内所耗费的燃料仅为ΔV＝2.4077m/s。以当前小推力发动机的工艺水平，燃气喷射速度可达到16.434km/s(SMART-1星载PPS-1350型Hall离子发动机)，则控制器在1年内耗费的燃料仅占总质量的1.06％。

可见，在从星的±x轴和±y轴方向分别安装连续推进的高比冲小推力发动机，本发明所用控制器即可以实现。

5、重复步骤(2)～(4)，可以使得编队飞行卫星以拟周期J₂不变相对轨道进行编队飞行。

根据步骤3中对受控系统稳定性的分析，式(22)描述的保Hamiltonian结构控制器可以生成稳定有界的相对轨道(称为稳定Lissajous轨道)，且生成条件适用于任意的初始状态。因此，本步骤可以使得编队飞行卫星以拟周期J₂不变相对轨道进行编队行。另外，式(22)描述的控制器仅需要相对位置Δr的反馈信息，因此十分适于无速度测量编队任务。

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

Claims (2)

1.一种拟周期J₂不变相对轨道编队飞行控制方法，其特征在于：以编队飞行卫星中飞行在标称轨道上的卫星为基准星，对编队飞行卫星中未飞行在标称轨道上的卫星施加控制力矩Tc进行调整，T_c＝B^-Δr-σ²G[u₊u₊ ^H+u_-u_- ^H]Δr，其中B^-＝(B-B^T)/2，Δr为基准星轨道坐标系下待调整卫星与基准星的相对位置矢量，B_λ＝(b₁₁+b₂₂+a₁₁a₂₂+a₁₂ ²)，C_λ＝(b₁₁b₂₂-b₁₂ ²)，a_ij表示矩阵A的第i行第j个元素，b_ij表示矩阵B⁺的第i行第j个元素，B⁺＝(B+B^T)/2，u₊、u_-分别满足[σ²I-σA+B]u₊＝0，[σ²I+σA+B]u_-＝0，G为控制增益，

$A=2\{\left[\begin{array}{ccc}0& -\cos \left(i\right)& \cos \left(u\right)\·\sin \left(i\right)\\ \cos \left(i\right)& 0& -\sin \left(i\right)\·\sin \left(u\right)\\ -\cos \left(u\right)\·\sin \left(i\right)& \sin \left(i\right)\·\sin \left(u\right)& 0\end{array}\right]\stackrel{\·}{\mathrm{\Ω}}+\left[\begin{array}{ccc}0& 0& -\sin \left(u\right)\\ 0& 0& -\cos \left(u\right)\\ \sin \left(u\right)& \cos \left(u\right)& 0\end{array}\right]\stackrel{\·}{i}+\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·}{u}\}$

$\begin{array}{c}B=\left[\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]\stackrel{\·\·}{u}+\left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& 0\end{array}\right]{\stackrel{\·}{u}}^2+\left[\begin{array}{ccc}\frac{1}{2}(-3+\cos \left(2u\right))& -\cos \left(u\right)\·\sin \left(u\right)& 0\\ -\cos \left(u\right)\·\sin \left(u\right)& \frac{1}{2}(-3-\cos \left(2u\right))& 0\\ 2\·\cos \left(u\right)& -2\·\sin \left(u\right)& -1\end{array}\right]\stackrel{\·}{u}\·\stackrel{\·}{i}\\ +\left[\begin{array}{ccc}\frac{1}{2}(1+\cos \left(i\right))(-3-\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos \left(2u\right)& \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& \cos \left(i\right)\·\sin \left(i\right)\·\sin \left(u\right)\\ \cos \left(u\right)\·{\sin }^2\left(i\right)\·\sin \left(u\right)& -\frac{1}{2}(1+\cos \left(i\right))(-3+\cos \left(i\right)-\cos \left(2u\right)+\cos \left(i\right)\·\cos (2u)& \cos \left(i\right)\·\sin \left(i\right)\·\cos \left(u\right)\\ (2+\cos \left(i\right))\·\sin \left(i\right)\·\sin \left(u\right)& (2+\cos \left(i\right))\·\sin \left(i\right)\·\cos \left(u\right)& -{\sin }^2\left(i\right)\end{array}\right]\stackrel{\·}{u}\\ +\frac{\mathrm{\μ}}{r^3}\left[\begin{array}{ccc}-2& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]+\frac{\mathrm{\μ}{J}_2{\mathrm{Re}}^2}{r^5}\left[\begin{array}{ccc}6(3\·{\sin }^2\left(i\right)\·{\sin }^2\left(u\right)-1)& -6\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& -6\·\sin \left(2i\right)\·\sin \left(u\right)\\ -\frac{9}{2}\·{\sin }^2\left(i\right)\·\sin \left(2u\right)& 3\·{\sin }^2\left(i\right)\·\cos \left(2u\right)& \frac{3}{2}\·\sin \left(2i\right)\·\cos \left(u\right)\\ -\frac{9}{2}\·\sin \left(2i\right)\·\sin \left(u\right)& 3\·\sin \left(2i\right)\·\cos \left(u\right)& 3\·\cos \left(2i\right)\end{array}\right]\end{array}\·\stackrel{\·}{\mathrm{\Ω}}$

u为主星纬度幅角，i为主星轨道倾角，Ω为主星升交点赤经，μ为地球的引力势能常数，Re为地球半径，J₂＝-1.08264×10^-3为常数，r为主星的矢径的模，I为两阶单位阵，角标T为矩阵的转置，字母上方的一点和两点分别表示一阶微分和二阶微分。

2.根据权利要求1所述的一种拟周期J₂不变相对轨道编队飞行控制方法，其特征在于：所述的控制增益G的选取方法为：令G从0开始逐渐增大，直至采用T_c＝B^-Δr-σ²G[u₊u₊ ^H+u_-u_- ^H]Δr描述的系统的所有Floquet乘子均分布在复域的单位圆上，则此时对应的G可以作为取值范围的下限，选取G时应大于此值。