CN104076819B - 一种圆参考轨道下卫星有界伴飞的边界控制方法 - Google Patents

Abstract

本发明提供一种圆参考轨道下卫星有界伴飞的边界控制方法，包括：当需要终端时刻实现跟踪卫星相对参考卫星的在有界边界间的有界伴飞构型时，预建立圆参考轨道下卫星有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式，求解得到与给定边界伴飞任务对应的期望轨道根数差；使用预给定的定量解析控制模型计算，得到初始时刻和终端时刻需要施加给跟踪卫星的脉冲速度增量。以边界解析模型求解期望轨道根数差，拓展了伴飞构型设计的空间，给编队或集群飞行任务构型的选取带来一定的灵活性。以定量解析控制模型精确计算得到具有给定边界的脉冲速度增量需求，具有控制精度高、控制过程简单的优点。

Description

一种圆参考轨道下卫星有界伴飞的边界控制方法

技术领域

本发明属于航天器轨道动力学和控制技术领域，具体涉及一种圆参考轨道下卫星有界伴飞的边界控制方法。

背景技术

多卫星在轨协同工作的空间任务，如卫星编队重力场测量、在轨服务、非合作目标监视等任务，需要使编队或集群飞行的卫星彼此保持在有限的距离范围内。伴飞，是一种以参考卫星为中心的轨迹封闭的周期性相对运动。当初始相对状态给定后，伴飞构型将具有确定的星间距离边界，可作为带有边界约束的多星协同空间任务的首选方案。传统上，由于缺失伴飞构型的边界理论模型，为实现具有期望星间距离边界的空间任务，轨道设计人员常将伴飞构型限定在少数几种简单构型上，例如：空间圆构型、轨道面投影圆构型、同轨道跟飞构型、Cartwheel构型、Pendulum构型等。由于可选的具有明确边界的伴飞构型数目较少，这大大限制了空间任务设计的灵活性。

为实现带有边界约束的伴飞构型，传统上采用多次脉冲推力控制或连续推力控制的方式，使编队或集群飞行系统里的各个卫星均保持在设定好的站位上。该种方法存在的主要问题为：由于传统上已知的具有解析边界表达式的伴飞构型类型较少，所设定的卫星控制目标较为保守，使得维持期望边界的相对运动具有控制消耗较大的缺点；另外，由于需要多次尝试进行脉冲推力控制或连续推力控制，也就是说，需要向被控制的卫星施加多次推力，具有控制过程复杂的问题。

发明内容

针对现有技术存在的缺陷，本发明提供一种圆参考轨道下卫星有界伴飞的边界控制方法，用以解决上述问题。

本发明采用的技术方案如下：

本发明提供一种圆参考轨道下卫星有界伴飞的边界控制方法，包括以下步骤：

S1，以自由飞行的圆轨道卫星作为参考卫星，以另一个在参考卫星附近运动的卫星为跟踪卫星；

S2，当下达以下任务要求时，执行S3：终端时刻t₁，实现跟踪卫星相对参考卫星的在d_max～d_min边界间的有界伴飞构型；其中，d_max指跟踪卫星和参考卫星的星间距离的最大值；d_min指跟踪卫星和参考卫星的星间距离的最小值；

S3，在当前的初始时刻t₀，测量得到跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数；

S4，预建立圆参考轨道下卫星有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式；

将d_max、d_min、跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数代入所述边界定量解析关系式，求解该边界定量解析关系式，得到与给定边界伴飞任务d_max和d_min对应的跟踪卫星与参考卫星之间的期望轨道根数差；

S5，以所述初始轨道根数差、所述期望轨道根数差、所述初始时刻t₀和所述终端时刻t₁作为已知的输入参数，使用预给定的定量解析控制模型进行计算，得到初始时刻t₀需要施加给跟踪卫星的初始脉冲速度增量，以及，还得到终端时刻t₁需要施加给跟踪卫星的终端脉冲速度增量；

S6，在初始时刻t₀，将S5计算得到的初始脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹；当达到终端时刻t₁时，将S5计算得到的终端脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹，进而形成满足给定边界d_max和d_min的有界伴飞构型。

优选的，S3中，所述参考卫星的初始轨道根数包括：轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω和初始轨道平近点角M₀；其中，轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω为常数；轨道平近点角M为变量，其在初始时刻的轨道平近点角记为初始轨道平近点角M₀；

所述初始轨道根数差包括：初始轨道半长轴差δa₀、初始轨道偏心率差δe₀、初始轨道倾角差δi₀、初始轨道升交点赤经差δΩ₀、初始轨道近地点幅角差δω₀和初始轨道平近点角差δM₀。

优选的，S4中，所建立的边界定量解析关系式为：

$d_{\max }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\δ}}^{2}e+({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)\\ +\sqrt{9{\mathrm{\δ}}^{4}e+{6\mathrm{\δ}}^{2}e({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)\cos 2\mathrm{\α}+{({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)}^2}\end{array}}---\left(1\right)$

$d_{\min }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\δ}}^{2}e+({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)\\ -\sqrt{9{\mathrm{\δ}}^{4}e+{6\mathrm{\δ}}^{2}e({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)\cos 2\mathrm{\α}+{({\mathrm{\δ}}^{2}i+{\mathrm{\δ}}^2\mathrm{\Ω}{\sin }^{2}i)}^2}\end{array}}---\left(2\right)$

其中：

$\mathrm{\α}=\mathrm{\ω}-\arctan \left(\frac{\mathrm{\δ\Ω}\sin i}{\mathrm{\δi}}\right)---\left(3\right)$

另外：

δa＝0，δM₀+δω+cosiδΩ＝0 (4)

上述公式(4)为公式(1)(2)的有效条件。

优选的，以所述边界定量解析关系式为基础，建立给定边界下的跟踪卫星与参考卫星之间的期望轨道根数差的解析关系式为：

$\left\{\begin{array}{c}{\mathrm{\δa}}_1=0\\ {\mathrm{\δe}}_1=\±\sqrt{x_1}\\ {\mathrm{\δi}}_1=\±\sqrt{\frac{x_2}{1+{\tan }^2(\mathrm{\ω}-\mathrm{\α})}}\\ {\mathrm{\δ\Ω}}_1=\frac{\tan (\mathrm{\ω}-\mathrm{\α})}{\sin i}{\mathrm{\δi}}_1\\ {\mathrm{\δM}}_1=-{\mathrm{\δ\ω}}_1-\cos i{\mathrm{\Ω}}_1\end{array}\right.---\left(5\right)$

其中，δa₁为期望轨道半长轴差、δe₁为期望轨道偏心率差、δi₁为期望轨道倾角差、δΩ₁为期望轨道升交点赤经差、δM₁为期望轨道平近点角差；δω₁为期望轨道近地点幅角差；

其中：δω₁可在0到2π之间任意取值；而x₁、x₂、x₃由以下公式计算：

$0\&le;{\mathrm{\&eta;}}_2<\frac{16}{25}:\left\{\begin{array}{c}-1\&le;x_3\&le;1\\ x_1=\frac{-({3x}_3-5)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)-\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\\ x_2=\frac{(9-{15x}_3)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)+5\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\end{array}\right.---\left(6\right)$

上述公式(6)～(7)中的参数x₃可在公式给定的范围内自由取值，参数η₂、分别为：

${\mathrm{\&eta;}}_2=\frac{{4d}_{\max }^2{d}_{\min }^2}{{(d_{\max }^2+d_{\min }^2)}^2}---\left(8\right)$

${\stackrel{\&OverBar;}{d}}_{\max }=\sqrt{2}{d}_{\max }/a---\left(9\right)$

${\stackrel{\&OverBar;}{d}}_{\min }=\sqrt{2}{d}_{\min }/a---\left(10\right)$

将d_max、d_min、参考卫星的轨道半长轴a、参考卫星的轨道近地点幅角ω、参考卫星的轨道倾角i各已知参数代入上述公式(5)～(10)，得到期望轨道根数差。

优选的，S5具体为：

S51，求和的值：

：跟踪卫星在t₀时刻施加脉冲速度增量控制后相对参考星的速度；

：跟踪卫星在t₁时刻施加脉冲速度增量控制前相对参考星的速度；

和的解析表达式为：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}={\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}^{-1}({\mathrm{\&rho;}}_1-{\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}{\mathrm{\&rho;}}_0)\\ {\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}={\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}{\mathrm{\&rho;}}_0+{\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}\end{array}\right.---\left(11\right)$

其中：

${\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}=\left[\begin{array}{ccc}4-3c& 0& 0\\ 6(s-\mathrm{\&psi;})& 1& 0\\ 0& 0& c\end{array}\right]---\left(12\right)$

${\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}s/n& 2(1-c)/n& 0\\ -2(1-c)/n& (4s-3\mathrm{\&psi;})/n& 0\\ 0& 0& s/n\end{array}\right]---\left(13\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}=\left[\begin{array}{ccc}3\mathrm{ns}& 0& 0\\ 6n(c-1)& 0& 0\\ 0& 0& -\mathrm{ns}\end{array}\right]---\left(14\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}c& 2s& 0\\ -2s& 4c-3& 0\\ 0& 0& c\end{array}\right]---\left(15\right)$

其中，公式(12-15)中出现的几个参数，其计算方式如下：

Δt＝t₁-t₀，ψ＝nΔt，s＝sin(ψ)，c＝cos(ψ)，μ＝3.986×10¹⁴

另外：初始时刻相对位置矢量ρ₀通过以下公式计算：

${\mathrm{\&rho;}}_0={[x_0,y_0{,z}_0]}^T=\left\{\begin{array}{c}{x}_0={\mathrm{\&delta;a}}_0-a\cos f_0{\mathrm{\&delta;e}}_0\\ y_0=a({\mathrm{\&delta;M}}_0+{\mathrm{\&delta;\&omega;}}_0+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0)+2a\sin f_0{\mathrm{\&delta;e}}_0\\ z_0=a\sin (f_0+\mathrm{\&omega;}){\mathrm{\&delta;i}}_0-a\cos (f_0+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0\end{array}\right.---\left(16\right)$

f₀＝M₀ (17)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δM₀、δω₀、δi₀、δΩ₀代入公式(16)～(17)，即可求得ρ₀值；

终端时刻相对位置矢量ρ₁通过以下公式计算：

${\mathrm{\&rho;}}_1={[x_1,y_1{,z}_1]}^T=\left\{\begin{array}{c}{x}_1={\mathrm{\&delta;a}}_1-a\cos f_1{\mathrm{\&delta;e}}_1\\ y_1=a({\mathrm{\&delta;M}}_1+{\mathrm{\&delta;\&omega;}}_1+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1)+2a\sin f_1{\mathrm{\&delta;e}}_1\\ z_1=a\sin (f_1+\mathrm{\&omega;}){\mathrm{\&delta;i}}_1-a\cos (f_1+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1\end{array}\right.---\left(18\right)$

f₁＝M₁＝M₀+n(t₁-t₀) (19)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δM₁、δω₁、δi₁、δΩ₁代入公式(18)～(19)，即可求得ρ₁值；

t₁和t₀的值作为已知值，通过公式(11)～(19)即可求得和的值；

S52，求解和值：

：跟踪卫星在t₀时刻施加脉冲速度增量控制前相对参考星的速度；

：跟踪卫星在t₁时刻施加脉冲速度增量控制后相对参考星的速度；

(1)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0={[{\stackrel{\&CenterDot;}{x}}_0,{\stackrel{\&CenterDot;}{y}}_0,{\stackrel{\&CenterDot;}{z}}_0]}^T---\left(20\right)$

其中：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_0=\mathrm{na}\sin f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{y}}_0=-\frac{3}{2}n{\mathrm{\&delta;a}}_0+2\mathrm{na}\cos f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{z}}_0=\mathrm{na}[\cos (\mathrm{\&omega;}+f_0){\mathrm{\&delta;i}}_0+\sin (\mathrm{\&omega;}+f_0)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0]\end{array}\right.---\left(21\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δi₀、δΩ₀、以及通过公式17计算得到的f₀代入公式(20)～(21)，即可求得值；

(2)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1={[{\stackrel{\&CenterDot;}{x}}_1,{\stackrel{\&CenterDot;}{y}}_1,{\stackrel{\&CenterDot;}{z}}_1]}^T---\left(22\right)$

其中，

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_1=\mathrm{na}\sin f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{y}}_1=-\frac{3}{2}n{\mathrm{\&delta;a}}_1+2\mathrm{na}\cos f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{z}}_1=\mathrm{na}[\cos (\mathrm{\&omega;}+f_1){\mathrm{\&delta;i}}_1+\sin (\mathrm{\&omega;}+f_1)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1]\end{array}\right.---\left(23\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δi₁、δΩ₁、以及通过公式(19)计算得到的f₁代入公式(22)～(23)，即可求得值；

其中，S51和S52不分先后顺序；

S53，将和代入公式(24)，求得Δv₀和Δv₁：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_0={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}\\ {\mathrm{\&Delta;v}}_1={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}\end{array}\right.---\left(24\right)$

其中：

Δv₀代表初始时刻t₀需要施加给跟踪卫星的脉冲速度增量；

Δv₁代表终端时刻t₁需要施加给跟踪卫星的脉冲速度增量。

优选的，终端时刻t₁通过以下方法确定：

如果伴飞控制任务要求形成伴飞构型的时间不超过T，则有界伴飞控制的终端时刻满足0＜t₁≤T，具体数值可通过优化方法搜索确定；

如果伴飞控制任务要求形成伴飞构型的时刻为T，则有界伴飞控制的终端时刻t₁＝T。

本发明提供的圆参考轨道下卫星有界伴飞的边界控制方法，具有以下优点：

本发明给出了圆参考轨道下伴飞构型的边界解析模型，能够定量给出具有给定星间距离边界的期望轨道根数差，从而大大拓展了有界伴飞构型的设计空间，给编队或集群飞行任务设计带来显著的灵活性。本发明以给定星间距离边界为控制目标，通过求解定量解析控制模型能够精确计算得到在初始时刻和终端时刻需要施加给跟踪卫星的脉冲速度增量，因此，只需要向跟踪卫星施加两次脉冲速度增量，即能形成满足给定边界的有界伴飞构型，具有控制精度高、控制过程简单的优点。

附图说明

图1为本发明提供的圆参考轨道下有界伴飞边界控制的过程原理图；

图2为两次脉冲速度增量施加前后相对运动的变化情况图；

图3为两脉冲速度增量施加前后的相对运动随时间的变化情况图；

图4为两脉冲速度增量实施结束后星间距离随参考星真近点角的变化情况图。

具体实施方式

以下结合附图对本发明进行详细说明：

本发明提供一种圆参考轨道下卫星有界伴飞的边界控制方法，只需要向跟踪卫星施加两次脉冲速度增量，即能形成满足给定边界的有界伴飞构型。如图1所示，为圆参考轨道下有界伴飞边界控制的过程原理图，其中，O_E代表地心，O₀为参考卫星质心，O₁为跟踪卫星质心，x为参考卫星轨道径向，z为参考卫星轨道面法向，y垂直于x和z方向并与x和z构成直角坐标系，r₀(t₀)、r₀(t₁)分别为初始时刻和终端时刻参考卫星的惯性位置矢量；r₁(t₀)、r₁(t₁)分别为初始时刻和终端时刻跟踪卫星的惯性位置矢量，Δv₀表示初始时刻施加的脉冲速度增量，Δv₁表示终端时刻施加的脉冲速度增量。

具体的，本发明提供一种圆参考轨道下卫星有界伴飞的边界控制方法，包括以下步骤：

S1，以自由飞行的圆轨道卫星作为参考卫星，以另一个在参考卫星附近运动的卫星为跟踪卫星；

S2，当下达以下任务要求时，执行S3：终端时刻t₁，实现跟踪卫星相对参考卫星的在d_max～d_min边界间的有界伴飞构型；其中，d_max指跟踪卫星和参考卫星的星间距离的最大值；d_min指跟踪卫星和参考卫星的星间距离的最小值；

具体的，终端时刻t₁通过以下方法确定：

如果伴飞控制任务要求形成伴飞构型的时间不超过T，则有界伴飞控制的终端时刻满足0＜t₁≤T，具体数值可通过优化方法搜索确定；

如果伴飞控制任务要求形成伴飞构型的时刻为T，则有界伴飞控制的终端时刻t₁＝T。

S3，在当前的初始时刻t₀，测量得到跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数；

其中，参考卫星的初始轨道根数包括：轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω和初始轨道平近点角M₀；其中，轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω为常数；轨道平近点角M为变量，其在初始时刻的轨道平近点角记为初始轨道平近点角M₀；

所述初始轨道根数差包括：初始轨道半长轴差δa₀、初始轨道偏心率差δe₀、初始轨道倾角差δi₀、初始轨道升交点赤经差δΩ₀、初始轨道近地点幅角差δω₀和初始轨道平近点角差δM₀。

S4，预建立圆参考轨道下卫星有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式；

将d_max、d_min、跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数代入所述边界定量解析关系式，求解该边界定量解析关系式，得到与给定边界伴飞任务d_max和d_min对应的跟踪卫星与参考卫星之间的期望轨道根数差；

具体的，建立的边界定量解析关系式为：

$d_{\max }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\\ +\sqrt{9{\mathrm{\&delta;}}^{4}e+{6\mathrm{\&delta;}}^{2}e({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\cos 2\mathrm{\&alpha;}+{({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)}^2}\end{array}}---\left(1\right)$

$d_{\min }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\\ -\sqrt{9{\mathrm{\&delta;}}^{4}e+{6\mathrm{\&delta;}}^{2}e({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\cos 2\mathrm{\&alpha;}+{({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)}^2}\end{array}}---\left(2\right)$

其中：

$\mathrm{\&alpha;}=\mathrm{\&omega;}-\arctan \left(\frac{\mathrm{\&delta;\&Omega;}\sin i}{\mathrm{\&delta;i}}\right)---\left(3\right)$

另外：

δa＝0，δM₀+δω+cosiδΩ＝0 (4)

上述公式(4)为公式(1)(2)的有效条件。

另外，以所述边界定量解析关系式为基础，建立给定边界下的跟踪卫星与参考卫星之间的期望轨道根数差的解析关系式为：

$\left\{\begin{array}{c}{\mathrm{\&delta;a}}_1=0\\ {\mathrm{\&delta;e}}_1=\&PlusMinus;\sqrt{x_1}\\ {\mathrm{\&delta;i}}_1=\&PlusMinus;\sqrt{\frac{x_2}{1+{\tan }^2(\mathrm{\&omega;}-\mathrm{\&alpha;})}}\\ {\mathrm{\&delta;\&Omega;}}_1=\frac{\tan (\mathrm{\&omega;}-\mathrm{\&alpha;})}{\sin i}{\mathrm{\&delta;i}}_1\\ {\mathrm{\&delta;M}}_1=-{\mathrm{\&delta;\&omega;}}_1-\cos i{\mathrm{\&Omega;}}_1\end{array}\right.---\left(5\right)$

其中，δa₁为期望轨道半长轴差、δe₁为期望轨道偏心率差、δi₁为期望轨道倾角差、δΩ₁为期望轨道升交点赤经差、δM₁为期望轨道平近点角差；δω₁为期望轨道近地点幅角差；

其中：δω₁可在0到2π之间任意取值；而x₁、x₂、x₃由以下公式计算：

$0\&le;{\mathrm{\&eta;}}_2<\frac{16}{25}:\left\{\begin{array}{c}-1\&le;x_3\&le;1\\ x_1=\frac{-({3x}_3-5)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)-\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\\ x_2=\frac{(9-{15x}_3)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)+5\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\end{array}\right.---\left(6\right)$

上述公式(6)～(7)中的参数x₃可在公式给定的范围内自由取值，参数η₂、分别为：

${\mathrm{\&eta;}}_2=\frac{{4d}_{\max }^2{d}_{\min }^2}{{(d_{\max }^2+d_{\min }^2)}^2}---\left(8\right)$

${\stackrel{\&OverBar;}{d}}_{\max }=\sqrt{2}{d}_{\max }/a---\left(9\right)$

${\stackrel{\&OverBar;}{d}}_{\min }=\sqrt{2}{d}_{\min }/a---\left(10\right)$

将d_max、d_min、参考卫星的轨道半长轴a、参考卫星的轨道近地点幅角ω、参考卫星的轨道倾角i各已知参数代入上述公式(5)～(10)，得到期望轨道根数差。

S5，以所述初始轨道根数差、所述期望轨道根数差、所述初始时刻t₀和所述终端时刻t₁作为已知的输入参数，使用预给定的定量解析控制模型进行计算，得到初始时刻t₀需要施加给跟踪卫星的初始脉冲速度增量，以及，还得到终端时刻t₁需要施加给跟踪卫星的终端脉冲速度增量；

本步骤具体包括以下步骤：

S51，求和的值：

：跟踪卫星在t₀时刻施加脉冲速度增量控制后相对参考星的速度；

：跟踪卫星在t₁时刻施加脉冲速度增量控制前相对参考星的速度；

和的解析表达式为：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}={\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}^{-1}({\mathrm{\&rho;}}_1-{\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}{\mathrm{\&rho;}}_0)\\ {\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}={\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}{\mathrm{\&rho;}}_0+{\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}\end{array}\right.---\left(11\right)$

其中：

${\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}=\left[\begin{array}{ccc}4-3c& 0& 0\\ 6(s-\mathrm{\&psi;})& 1& 0\\ 0& 0& c\end{array}\right]---\left(12\right)$

${\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}s/n& 2(1-c)/n& 0\\ -2(1-c)/n& (4s-3\mathrm{\&psi;})/n& 0\\ 0& 0& s/n\end{array}\right]---\left(13\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}=\left[\begin{array}{ccc}3\mathrm{ns}& 0& 0\\ 6n(c-1)& 0& 0\\ 0& 0& -\mathrm{ns}\end{array}\right]---\left(14\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}c& 2s& 0\\ -2s& 4c-3& 0\\ 0& 0& c\end{array}\right]---\left(15\right)$

其中，公式(12-15)中出现的几个参数，其计算方式如下：

Δt＝t₁-t₀，ψ＝nΔt，s＝sin(ψ)，c＝cos(ψ)，μ＝3.986×10¹⁴

另外：初始时刻相对位置矢量ρ₀通过以下公式计算：

${\mathrm{\&rho;}}_0={[x_0,y_0{,z}_0]}^T=\left\{\begin{array}{c}{x}_0={\mathrm{\&delta;a}}_0-a\cos f_0{\mathrm{\&delta;e}}_0\\ y_0=a({\mathrm{\&delta;M}}_0+{\mathrm{\&delta;\&omega;}}_0+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0)+2a\sin f_0{\mathrm{\&delta;e}}_0\\ z_0=a\sin (f_0+\mathrm{\&omega;}){\mathrm{\&delta;i}}_0-a\cos (f_0+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0\end{array}\right.---\left(16\right)$

f₀＝M₀ (17)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δM₀、δω₀、δi₀、δΩ₀代入公式(16)～(17)，即可求得ρ₀值；

终端时刻相对位置矢量ρ₁通过以下公式计算：

${\mathrm{\&rho;}}_1={[x_1,y_1{,z}_1]}^T=\left\{\begin{array}{c}{x}_1={\mathrm{\&delta;a}}_1-a\cos f_1{\mathrm{\&delta;e}}_1\\ y_1=a({\mathrm{\&delta;M}}_1+{\mathrm{\&delta;\&omega;}}_1+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1)+2a\sin f_1{\mathrm{\&delta;e}}_1\\ z_1=a\sin (f_1+\mathrm{\&omega;}){\mathrm{\&delta;i}}_1-a\cos (f_1+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1\end{array}\right.---\left(18\right)$

f₁＝M₁＝M₀+n(t₁-t₀) (19)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δM₁、δω₁、δi₁、δΩ₁代入公式(18)～(19)，即可求得ρ₁值；

t₁和t₀的值作为已知值，通过公式(11)～(19)即可求得和的值；

S52，求解和值：

：跟踪卫星在t₀时刻施加脉冲速度增量控制前相对参考星的速度；

：跟踪卫星在t₁时刻施加脉冲速度增量控制后相对参考星的速度；

(1)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0={[{\stackrel{\&CenterDot;}{x}}_0,{\stackrel{\&CenterDot;}{y}}_0,{\stackrel{\&CenterDot;}{z}}_0]}^T---\left(20\right)$

其中：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_0=\mathrm{na}\sin f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{y}}_0=-\frac{3}{2}n{\mathrm{\&delta;a}}_0+2\mathrm{na}\cos f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{z}}_0=\mathrm{na}[\cos (\mathrm{\&omega;}+f_0){\mathrm{\&delta;i}}_0+\sin (\mathrm{\&omega;}+f_0)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0]\end{array}\right.---\left(21\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δi₀、δΩ₀、以及通过公式17计算得到的f₀代入公式(20)～(21)，即可求得值；

(2)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1={[{\stackrel{\&CenterDot;}{x}}_1,{\stackrel{\&CenterDot;}{y}}_1,{\stackrel{\&CenterDot;}{z}}_1]}^T---\left(22\right)$

其中，

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_1=\mathrm{na}\sin f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{y}}_1=-\frac{3}{2}n{\mathrm{\&delta;a}}_1+2\mathrm{na}\cos f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{z}}_1=\mathrm{na}[\cos (\mathrm{\&omega;}+f_1){\mathrm{\&delta;i}}_1+\sin (\mathrm{\&omega;}+f_1)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1]\end{array}\right.---\left(23\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δi₁、δΩ₁、以及通过公式(19)计算得到的f₁代入公式(22)～(23)，即可求得值；

其中，S51和S52不分先后顺序；

S53，将和代入公式(24)，求得Δv₀和Δv₁：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_0={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}\\ {\mathrm{\&Delta;v}}_1={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}\end{array}\right.---\left(24\right)$

其中：

Δv₀代表初始时刻t₀需要施加给跟踪卫星的脉冲速度增量；

Δv₁代表终端时刻t₁需要施加给跟踪卫星的脉冲速度增量。

S6，在初始时刻t₀，将S5计算得到的初始脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹；当达到终端时刻t₁时，将S5计算得到的终端脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹，进而形成满足给定边界d_max和d_min的有界伴飞构型。

本发明提供的圆参考轨道下卫星有界伴飞的边界控制方法，具有以下优点：

(1)显式地表示出一般条件下伴飞构型边界与卫星初始轨道根数差的解析关系，能够定量给出具有给定星间距离边界的期望轨道根数差，从而大大拓展了有界伴飞构型的设计空间，给编队或集群飞行任务设计带来显著的灵活性。

(2)显示地给出确定边界约束下的两脉冲速度增量控制规律，以给定星间距离边界为控制目标，通过求解定量解析控制模型能够精确计算得到在初始时刻和终端时刻需要施加给跟踪卫星的脉冲速度增量，因此，只需要向跟踪卫星施加两次脉冲速度增量，即能形成满足给定边界的有界伴飞构型，具有控制精度高、控制过程简单的优点。

以上两个优点为有界编队或集群飞行的构型设计及控制系统研发奠定了理论基础。

下面举例对本发明进一步详细说明：

一颗近地圆轨道卫星作为参考卫星自由飞行，另一颗卫星在参考卫星附近运动，其为跟踪卫星；在初始时刻t₀＝0s，参考卫星的初始轨道根数、跟踪卫星的初始轨道根数以及跟踪卫星与参考卫星之间的初始轨道根数差如表1所示。根据某次伴飞任务的需求，要求跟踪卫星在t₁＝4023s后实现相对参考卫星星间距离上界50km、下界10km的伴飞，即：d_max＝50km；d_min＝10km。

表1

为实现50km上界、10km下界的有界伴飞运动，首先推导有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式，即公式(1)～(2)：

(一)有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式

根据卫星编队动力学的一般理论，跟踪卫星相对参考卫星的相对位置和相对速度在参考卫星当地轨道坐标系下的坐标分量形式可表示为：

$\left\{\begin{array}{c}x=\mathrm{\&delta;a}-a\cos \mathrm{f\&delta;e}\\ y=-\frac{3}{2}n\&CenterDot;(t-t_0)\mathrm{\&delta;a}+a({\mathrm{\&delta;M}}_0+\mathrm{\&delta;\&omega;}+\cos \mathrm{i\&delta;\&Omega;})+2a\sin \mathrm{f\&delta;e}\\ z=a\sin (f+\mathrm{\&omega;})\mathrm{\&delta;i}-a\cos (f+\mathrm{\&omega;})\sin \mathrm{i\&delta;\&Omega;}\end{array}\right.---\left(100\right)$

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{x}=\mathrm{na}\sin \mathrm{f\&delta;e}\\ \stackrel{\&CenterDot;}{y}=-\frac{3}{2}\mathrm{n\&delta;a}+2\mathrm{na}\cos \mathrm{f\&delta;e}\\ \stackrel{\&CenterDot;}{z}=\mathrm{na}[\cos (\mathrm{\&omega;}+f)\mathrm{\&delta;i}+\sin (\mathrm{\&omega;}+f)\sin \mathrm{i\&delta;\&Omega;}]\end{array}\right.---\left(101\right)$

f＝M＝M₀+n(t-t₀) (102)

其中，a、e、i、Ω、ω、M₀为参考卫星的初始轨道根数；δa、δe、δi、δΩ、δω、δM₀为跟踪卫星和参考卫星之间的初始轨道根数差。

由公式(100)中y方向运动表达式可知，当两星轨道半长轴偏差不为零时，在y方向存在长期发散的趋势。在相对运动有界性问题中，需增加δa＝0的约束条件，此时，相对运动公式简化为：

$\left\{\begin{array}{c}x=-a\cos \mathrm{f\&delta;e}\\ y=a(\mathrm{\&delta;}{M}_0+\mathrm{\&delta;\&omega;}+\cos \mathrm{i\&delta;\&Omega;})+2a\sin \mathrm{f\&delta;e}\\ z=a\sin (f+\mathrm{\&omega;})\mathrm{\&delta;i}-a\cos (f+\mathrm{\&omega;})\sin \mathrm{i\&delta;\&Omega;}\end{array}\right.---\left(103\right)$

对于伴飞运动，运动的中心应该是参考卫星质心，为此，需使y方向的偏移项a(δM₀+δω+cosiδΩ)为零，此时伴飞相对运动的表达式化简为：

$\left\{\begin{array}{c}{x}^{\&prime;}=-a\cos \mathrm{f\&delta;e}\\ y^{\&prime;}=2a\sin \mathrm{f\&delta;e}\\ z^{\&prime;}=a\sin (f+\mathrm{\&omega;})\mathrm{\&delta;i}-a\cos (f+\mathrm{\&omega;})\sin \mathrm{i\&delta;\&Omega;}\end{array}\right.---\left(104\right)$

跟踪卫星相对参考卫星的星间距离可由下式表示：

$d=\sqrt{x^{\&prime;2}+y^{\&prime;2}+z^{\&prime;2}}---\left(105\right)$

将(104)式带入(105)式，可得到星间距离的具体形式：

$d=a\sqrt{c_1{\sin }^{2}f+c_2\sin f\cos f+c_3}---\left(106\right)$

其中：

c₁＝3δ²e+(δ²i+δ²Ωsin²i)cos2α (107)

c₂＝(δ²i+δ²Ωsin²i)sin2α (108)

c₃＝δ²e+(δ²i+δ²Ωsin²i)sin²α (109)

$\mathrm{\&alpha;}=\mathrm{\&omega;}-\arctan \left(\frac{\mathrm{\&delta;\&Omega;}\sin i}{\mathrm{\&delta;i}}\right)---\left(3\right)$

上述星间距离表达式(106)是变量f的函数，为寻找星间距离的极值，记s＝sinf，记则星间距为：

$d=a\sqrt{g+c_3}---\left(110\right)$

容易证明函数g的最大值和最小值分别为

$g_{\max }=\frac{1}{2}(c_1+\sqrt{c_1^2+c_2^2})---\left(111\right)$

$g_{\min }=\frac{1}{2}(c_1-\sqrt{c_1^2+c_2^2})---\left(112\right)$

因此d的最大值和最小值分别为：

$d_{\max }=a\sqrt{\frac{1}{2}(c_1+\sqrt{c_1^2+c_2^2})+c_3}---\left(113\right)$

$d_{\min }=a\sqrt{\frac{1}{2}(c_1-\sqrt{c_1^2+c_2^2})+c_3}---\left(114\right)$

带入c₁、c₂、c₃的表达式，可得到星间距离边界具体表达式：

$d_{\max }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\\ +\sqrt{9{\mathrm{\&delta;}}^{4}e+{6\mathrm{\&delta;}}^{2}e({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\cos 2\mathrm{\&alpha;}+{({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)}^2}\end{array}}---\left(1\right)$

$d_{\min }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}{5\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\\ -\sqrt{9{\mathrm{\&delta;}}^{4}e+{6\mathrm{\&delta;}}^{2}e({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)\cos 2\mathrm{\&alpha;}+{({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2\mathrm{\&Omega;}{\sin }^{2}i)}^2}\end{array}}---\left(2\right)$

其中：

a—参考卫星的轨道半长轴

e—参考卫星的偏心率

i—参考卫星的轨道倾角

Ω—参考卫星的升交点赤经

ω—参考卫星的近地点幅角

M—参考卫星的平近点角

δa为轨道半长轴差、δe为轨道偏心率差、δi为轨道倾角差、δΩ为轨道升交点赤经差、δM为轨道平近点角差；δω为轨道近地点幅角差；

另外：

δa＝0，δM₀+δω+cosiδΩ＝0 (4)

公式(4)为公式(1)(2)的有效条件。

(二)求解得到期望轨道根数差

以边界定量解析关系式为基础，建立给定边界下的跟踪卫星与参考卫星之间的期望轨道根数差的解析关系式为：

$\left\{\begin{array}{c}{\mathrm{\&delta;a}}_1=0\\ {\mathrm{\&delta;e}}_1=\&PlusMinus;\sqrt{x_1}\\ {\mathrm{\&delta;i}}_1=\&PlusMinus;\sqrt{\frac{x_2}{1+{\tan }^2(\mathrm{\&omega;}-\mathrm{\&alpha;})}}\\ {\mathrm{\&delta;\&Omega;}}_1=\frac{\tan (\mathrm{\&omega;}-\mathrm{\&alpha;})}{\sin i}{\mathrm{\&delta;i}}_1\\ {\mathrm{\&delta;M}}_1=-{\mathrm{\&delta;\&omega;}}_1-\cos i{\mathrm{\&Omega;}}_1\end{array}\right.---\left(5\right)$

其中，δa₁为期望轨道半长轴差、δe₁为期望轨道偏心率差、δi₁为期望轨道倾角差、δΩ₁为期望轨道升交点赤经差、δM₁为期望轨道平近点角差；δω₁为期望轨道近地点幅角差；

其中：δω₁可在0到2π之间任意取值；而x₁、x₂、x₃由以下公式计算：

$0\&le;{\mathrm{\&eta;}}_2<\frac{16}{25}:\left\{\begin{array}{c}-1\&le;x_3\&le;1\\ x_1=\frac{-({3x}_3-5)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)-\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\\ x_2=\frac{(9-{15x}_3)({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)+5\sqrt{{({3x}_3-5)}^2{({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2)}^2-4(34-{30x}_3){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2(34-{30x}_3)}\end{array}\right.---\left(6\right)$

上述公式(6)～(7)中的参数x₃可在公式给定的范围内自由取值，参数η₂、分别为：

${\mathrm{\&eta;}}_2=\frac{{4d}_{\max }^2{d}_{\min }^2}{{(d_{\max }^2+d_{\min }^2)}^2}---\left(8\right)$

${\stackrel{\&OverBar;}{d}}_{\max }=\sqrt{2}{d}_{\max }/a---\left(9\right)$

${\stackrel{\&OverBar;}{d}}_{\min }=\sqrt{2}{d}_{\min }/a---\left(10\right)$

将d_max、d_min、参考卫星的轨道半长轴a、参考卫星的轨道近地点幅角ω、参考卫星的轨道倾角i各已知参数代入上述公式(5)～(10)，得到期望轨道根数差。

根据公式(8)和给定的边界d_max＝50km、d_min＝10km，计算得到η₂＝0.148，可知选取公式(6)计算期望轨道根数。由于x₃是个自由变量，本实施例取值为1，则计算得到相应的x₁＝2.06×10^-6，x₂＝4.33×10^-5。进而通过公式(5)，可计算得到期望轨道根数差的具体数值为：

δa₁＝0，δe₁＝0.0014，δi₁＝0.0275°，δΩ₁＝0.379°，δω₁＝0.1°，δM₁＝-0.049°。

(三)初始时刻脉冲速度增量和终端时刻脉冲速度增量的计算

在初始时刻t₀，由表1可知，跟踪卫星相对参考卫星的轨道根数差为：

δσ(t₀)＝[10000,0.001,0.1°,0.1°,0.1°,0.1°]^T (115)

边界控制的要求为：在时刻t₁，跟踪卫星相对参考卫星的轨道根数差满足下述(116)式约束，即：

δσ(t₁)＝[0,0.0014,0.0275°,0.379°,0.1°,-0.049°]^T (116)

参考卫星是圆轨道，则初始时刻t₀＝0和终端时刻t₁＝4023s其轨道根数分别为：

σ₀(t₀)＝[6961181,0,97.622°,262.154°,265.810°,323.848°]^T (117)

σ₀(t₁)＝[6961181,0,97.622°,262.154°,265.810°,214.4°]^T (118)

初始时刻和终端时刻跟踪卫星的轨道根数分别为：

σ₁(t₀)＝σ₀(t₀)+δσ(t₀) (119)

σ₁(t₁)＝σ₀(t₁)+δσ(t₁) (120)

为使跟踪卫星从t₀时刻的轨道根数σ₁(t₀)变为t₁时刻的轨道根数σ₁(t₁)，可以通过相对运动的C-W方程描述经过两次脉冲速度增量后的运动转移情况。为此，首先需要将跟踪卫星在初始时刻和终端时刻的轨道根数转化为Hill系下的相对状态，具体的，通过下列方法计算得到t₀时刻和t₁时刻的相对位置矢量和相对速度矢量：

初始时刻相对位置矢量ρ₀通过以下公式计算：

${\mathrm{\&rho;}}_0={[x_0,y_0{,z}_0]}^T=\left\{\begin{array}{c}{x}_0={\mathrm{\&delta;a}}_0-a\cos f_0{\mathrm{\&delta;e}}_0\\ y_0=a({\mathrm{\&delta;M}}_0+{\mathrm{\&delta;\&omega;}}_0+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0)+2a\sin f_0{\mathrm{\&delta;e}}_0\\ z_0=a\sin (f_0+\mathrm{\&omega;}){\mathrm{\&delta;i}}_0-a\cos (f_0+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0\end{array}\right.---\left(16\right)$

f₀＝M₀ (17)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δM₀、δω₀、δi₀、δΩ₀代入公式(16)～(17)，即可求得ρ₀值；

终端时刻相对位置矢量ρ₁通过以下公式计算：

${\mathrm{\&rho;}}_1={[x_1,y_1{,z}_1]}^T=\left\{\begin{array}{c}{x}_1={\mathrm{\&delta;a}}_1-a\cos f_1{\mathrm{\&delta;e}}_1\\ y_1=a({\mathrm{\&delta;M}}_1+{\mathrm{\&delta;\&omega;}}_1+\cos \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1)+2a\sin f_1{\mathrm{\&delta;e}}_1\\ z_1=a\sin (f_1+\mathrm{\&omega;}){\mathrm{\&delta;i}}_1-a\cos (f_1+\mathrm{\&omega;})\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1\end{array}\right.---\left(18\right)$

f₁＝M₁＝M₀+n(t₁-t₀) (19)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δM₁、δω₁、δi₁、δΩ₁代入公式(18)～(19)，即可求得ρ₁值；

初始时刻相对速度矢量通过以下公式计算：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0{=[{\stackrel{\&CenterDot;}{x}}_0,{\stackrel{\&CenterDot;}{y}}_0,{\stackrel{\&CenterDot;}{z}}_0]}^T$

其中：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_0=\mathrm{na}\sin f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{y}}_0=-\frac{3}{2}n{\mathrm{\&delta;a}}_0+2\mathrm{na}\cos f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{z}}_0=\mathrm{na}[\cos (\mathrm{\&omega;}+f_0){\mathrm{\&delta;i}}_0+\sin (\mathrm{\&omega;}+f_0)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_0]\end{array}\right.---\left(21\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δi₀、δΩ₀、以及通过公式17计算得到的f₀代入上式，即可求得值；

终端时刻相对速度矢量通过以下公式计算：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1={[{\stackrel{\&CenterDot;}{x}}_1,{\stackrel{\&CenterDot;}{y}}_1,{\stackrel{\&CenterDot;}{z}}_1]}^T$

其中，

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_1=\mathrm{na}\sin f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{y}}_1=-\frac{3}{2}n{\mathrm{\&delta;a}}_1+2\mathrm{na}\cos f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{z}}_1=\mathrm{na}[\cos (\mathrm{\&omega;}+f_1){\mathrm{\&delta;i}}_1+\sin (\mathrm{\&omega;}+f_1)\sin \mathrm{i\&delta;}{\mathrm{\&Omega;}}_1]\end{array}\right.---\left(23\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δi₁、δΩ₁、以及通过公式(19)计算得到的f₁代入上式，即可求得值；

因此，将表1的轨道根数分别代入上式，可计算得到初始时刻和终端时刻跟踪卫星的相对位置和相对速度，具体计算结果见表2。

表2

参数	x(m)	y(m)	z(m)	vx(m/s)	vy(m/s)	vz(m/s)
							t0时刻	-5001.9	766.2	-1464.6	-11.9	16.3	-18.5
t1时刻	8250	-11302	25897	-6	-18	41

令Δt＝t₁-t₀，ψ＝nΔt，s＝sin(ψ)，c＝cos(ψ)，则相对运动的状态转移矩阵可表示为：

${\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}=\left[\begin{array}{ccc}4-3c& 0& 0\\ 6(s-\mathrm{\&psi;})& 1& 0\\ 0& 0& c\end{array}\right]---\left(12\right)$

${\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}s/n& 2(1-c)/n& 0\\ -2(1-c)/n& (4s-3\mathrm{\&psi;})/n& 0\\ 0& 0& s/n\end{array}\right]---\left(13\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}=\left[\begin{array}{ccc}3\mathrm{ns}& 0& 0\\ 6n(c-1)& 0& 0\\ 0& 0& -\mathrm{ns}\end{array}\right]---\left(14\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}c& 2s& 0\\ -2s& 4c-3& 0\\ 0& 0& c\end{array}\right]---\left(15\right)$

则有

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}={\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}^{-1}({\mathrm{\&rho;}}_1-{\mathrm{\&phi;}}_{\mathrm{\&rho;\&rho;}}{\mathrm{\&rho;}}_0)\\ {\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}={\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}{\mathrm{\&rho;}}_0+{\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}\end{array}\right.---\left(11\right)$

最后，可得到在t₀时刻和t₁时刻分别施加的脉冲速度增量分别为

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_0={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}\\ {\mathrm{\&Delta;v}}_1={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}\end{array}\right.---\left(24\right)$

带入相对位置和相对速度的具体数值，计算结果见表3：

表3

参数	vx(m/s)	vy(m/s)	vz(m/s)
				t0时刻	6.99	-4.48	-10.76
t1时刻	-0.86	-0.94	32.85

如图2所示，为两次脉冲速度增量施加前后相对运动的变化情况图。其中，曲线C代表未施加脉冲速度增量控制前，初始相对运动轨迹；曲线B代表施加初始脉冲速度增量控制后形成的转移轨道；曲线A代表施加两次脉冲速度增量控制后形成的伴飞轨道；其中，点D为第一次脉冲施加点；点E为第二次脉冲施加点。由图可知，在未施加脉冲速度增量控制前，相对运动轨迹是发散的、不封闭的，即图2中的曲线C；而在完成两次脉冲速度增量后，相对运动进入封闭的周期性伴飞轨道，即图2中的曲线A；因此，所施加的两次脉冲速度增量控制是有效的。如图3所示，为两脉冲速度增量施加前后的相对运动随时间的变化情况图，其中，前三幅子图表示三个相对位置坐标分量随时间的变化，最后一个子图表示星间距离的随时间的变化。由图3可知，施加脉冲速度增量控制后，相对运动在三个坐标分量方向及星间距离上均是周期的、有界的。如图4所示，为两脉冲速度增量实施结束后星间距离随参考星真近点角的变化情况图，由图4可知，星间距离在10km～50km之间周期性震荡，说明形成的伴飞运动是有界的，且边界与预期相同，从而进一步说明了本发明提供的圆参考轨道下卫星有界伴飞的边界控制方法的有效性。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视本发明的保护范围。

Claims (2)

1.一种圆参考轨道下卫星有界伴飞的边界控制方法，其特征在于，包括以下步骤：

S1，以自由飞行的圆轨道卫星作为参考卫星，以另一个在参考卫星附近运动的卫星为跟踪卫星；

S2，当下达以下任务要求时，执行S3：终端时刻t₁，实现跟踪卫星相对参考卫星的在d_max～d_min边界间的有界伴飞构型；其中，d_max指跟踪卫星和参考卫星的星间距离的最大值；d_min指跟踪卫星和参考卫星的星间距离的最小值；

S3，在当前的初始时刻t₀，测量得到跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数；

其中，S3中，所述参考卫星的初始轨道根数包括：轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω和初始轨道平近点角M₀；其中，轨道半长轴a、轨道偏心率e、轨道倾角i、轨道升交点赤经Ω、轨道近地点幅角ω为常数；轨道平近点角M为变量，其在初始时刻的轨道平近点角记为初始轨道平近点角M₀；

所述初始轨道根数差包括：初始轨道半长轴差δa₀、初始轨道偏心率差δe₀、初始轨道倾角差δi₀、初始轨道升交点赤经差δΩ₀、初始轨道近地点幅角差δω₀和初始轨道平近点角差δM₀；

S4，预建立圆参考轨道下卫星有界伴飞的边界与卫星轨道根数差之间的边界定量解析关系式；

将d_max、d_min、跟踪卫星与参考卫星之间的初始轨道根数差和参考卫星的初始轨道根数代入所述边界定量解析关系式，求解该边界定量解析关系式，得到与给定边界伴飞任务d_max和d_min对应的跟踪卫星与参考卫星之间的期望轨道根数差；

S4中，所建立的边界定量解析关系式为：

$d_{max}=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}5{\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i)\\ +\sqrt{9{\mathrm{\&delta;}}^{4}e+6{\mathrm{\&delta;}}^{2}e\left({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i\right)\cos 2\mathrm{\&alpha;}+{\left({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i\right)}^2}\end{array}}---\left(1\right)$

$d_{\min }=\frac{\sqrt{2}}{2}a\sqrt{\begin{array}{c}5{\mathrm{\&delta;}}^{2}e+({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i)\\ -\sqrt{9{\mathrm{\&delta;}}^{4}e+6{\mathrm{\&delta;}}^{2}e\left({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i\right)\cos 2\mathrm{\&alpha;}+{\left({\mathrm{\&delta;}}^{2}i+{\mathrm{\&delta;}}^2{\mathrm{\&Omega;sin}}^{2}i\right)}^2}\end{array}}---\left(2\right)$

其中：

$\mathrm{\&alpha;}=\mathrm{\&omega;}-\arctan \left(\frac{\mathrm{\&delta;}\mathrm{\&Omega;}\sin i}{\mathrm{\&delta;}i}\right)---\left(3\right)$

另外：

δa＝0，δM₀+δω+cosiδΩ＝0 (4)

上述公式(4)为公式(1)(2)的有效条件；

以所述边界定量解析关系式为基础，建立给定边界下的跟踪卫星与参考卫星之间的期望轨道根数差的解析关系式为：

$\left\{\begin{array}{c}{\mathrm{\&delta;a}}_1=0\\ {\mathrm{\&delta;e}}_1=\&PlusMinus;\sqrt{x_1}\\ {\mathrm{\&delta;i}}_1=\&PlusMinus;\sqrt{\frac{x_2}{1+{\tan }^2(\mathrm{\&omega;}-\mathrm{\&alpha;})}}\\ {\mathrm{\&delta;\&Omega;}}_1=\frac{tan(\mathrm{\&omega;}-\mathrm{\&alpha;})}{\sin i}{\mathrm{\&delta;i}}_1\\ {\mathrm{\&delta;M}}_1=-{\mathrm{\&delta;\&omega;}}_1-\cos {\mathrm{i\&delta;\&Omega;}}_1\end{array}\right.---\left(5\right)$

其中，δa₁为期望轨道半长轴差、δe₁为期望轨道偏心率差、δi₁为期望轨道倾角差、δΩ₁为期望轨道升交点赤经差、δM₁为期望轨道平近点角差；δω₁为期望轨道近地点幅角差；

其中：δω₁可在0到2π之间任意取值；而x₁、x₂、x₃由以下公式计算：

$0\&le;{\mathrm{\&eta;}}_2<\frac{16}{25}:\left\{\begin{array}{c}-1\&le;x_3\&le;1\\ x_1=\frac{-\left(3x_3-5\right)\left({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2\right)-\sqrt{{\left(3x_3-5\right)}^2{\left({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2\right)}^2-4\left(34-30x_3\right){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2\left(34-30x_3\right)}\\ x_2=\frac{\left(9-15x_3\right)\left({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2\right)+5\sqrt{{\left(3x_3-5\right)}^2{\left({\stackrel{\&OverBar;}{d}}_{\max }^2+{\stackrel{\&OverBar;}{d}}_{\min }^2\right)}^2-4\left(34-30x_3\right){\stackrel{\&OverBar;}{d}}_{\max }^2{\stackrel{\&OverBar;}{d}}_{\min }^2}}{2\left(34-30x_3\right)}\end{array}\right.---\left(6\right)$

上述公式(6)～(7)中的参数x₃可在公式给定的范围内自由取值，参数η₂、分别为：

${\mathrm{\&eta;}}_2=\frac{4d_{max}^2{d}_{\min }^2}{{(d_{\max }^2+d_{\min }^2)}^2}---\left(8\right)$

${\stackrel{\&OverBar;}{d}}_{\max }=\sqrt{2}{d}_{\max }/a---\left(9\right)$

${\stackrel{\&OverBar;}{d}}_{\min }=\sqrt{2}{d}_{\min }/a---\left(10\right)$

将d_max、d_min、参考卫星的轨道半长轴a、参考卫星的轨道近地点幅角ω、参考卫星的轨道倾角i各已知参数代入上述公式(5)～(10)，得到期望轨道根数差；

S5，以所述初始轨道根数差、所述期望轨道根数差、所述初始时刻t₀和所述终端时刻t₁作为已知的输入参数，使用预给定的定量解析控制模型进行计算，得到初始时刻t₀需要施加给跟踪卫星的初始脉冲速度增量，以及，还得到终端时刻t₁需要施加给跟踪卫星的终端脉冲速度增量；

S5具体为：

S51，求和的值：

跟踪卫星在t₀时刻施加脉冲速度增量控制后相对参考星的速度；

跟踪卫星在t₁时刻施加脉冲速度增量控制前相对参考星的速度；

和的解析表达式为：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}={\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}^{-1}({\mathrm{\&rho;}}_1-{\mathrm{\&phi;}}_{\mathrm{\&rho;}\mathrm{\&rho;}}{\mathrm{\&rho;}}_0)\\ {\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}={\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}{\mathrm{\&rho;}}_0+{\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}\end{array}\right.---\left(11\right)$

其中：

${\mathrm{\&phi;}}_{\mathrm{\&rho;}\mathrm{\&rho;}}=\left[\begin{array}{ccc}4-3c& 0& 0\\ 6(s-\mathrm{\&psi;})& 1& 0\\ 0& 0& c\end{array}\right]---\left(12\right)$

${\mathrm{\&phi;}}_{\mathrm{\&rho;}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}s/n& 2(1-c)/n& 0\\ -2(1-c)/n& (4s-3\mathrm{\&psi;})/n& 0\\ 0& 0& s/n\end{array}\right]---\left(13\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\mathrm{\&rho;}}=\left[\begin{array}{ccc}3ns& 0& 0\\ 6n(c-1)& 0& 0\\ 0& 0& -ns\end{array}\right]---\left(14\right)$

${\mathrm{\&phi;}}_{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}=\left[\begin{array}{ccc}c& 2s& 0\\ -2s& 4c-3& 0\\ 0& 0& c\end{array}\right]---\left(15\right)$

其中，公式(12-15)中出现的几个参数，其计算方式如下：

$\mathrm{\&Delta;}t=t_1-t_0,\mathrm{\&psi;}=n\mathrm{\&Delta;}t,s=sin\left(\mathrm{\&psi;}\right),c=cos\left(\mathrm{\&psi;}\right),n=\sqrt{\mathrm{\&mu;}/a^3},\mathrm{\&mu;}=3.986\&times;{10}^{14}$

另外：初始时刻相对位置矢量ρ₀通过以下公式计算：

${\mathrm{\&rho;}}_0={\&lsqb;x_0,y_0,z_0\&rsqb;}^T=\left\{\begin{array}{c}{x}_0={\mathrm{\&delta;a}}_0-a\cos f_0{\mathrm{\&delta;e}}_0\\ y_0=a\left({\mathrm{\&delta;M}}_0+{\mathrm{\&delta;\&omega;}}_0+\cos {\mathrm{i\&delta;\&Omega;}}_0\right)+2a\sin f_0{\mathrm{\&delta;e}}_0\\ z_0=a\sin \left(f_0+\mathrm{\&omega;}\right){\mathrm{\&delta;i}}_i-a\cos \left(f_0+\mathrm{\&omega;}\right)\sin {\mathrm{i\&delta;\&Omega;}}_0\end{array}\right.---\left(16\right)$

f₀＝M₀ (17)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δM₀、δω₀、δi₀、δΩ₀代入公式(16)～(17)，即可求得ρ₀值；

终端时刻相对位置矢量ρ₁通过以下公式计算：

${\mathrm{\&rho;}}_1={\&lsqb;x_1,y_1,z_1\&rsqb;}^T=\left\{\begin{array}{c}{x}_1={\mathrm{\&delta;a}}_1-a\cos f_1{\mathrm{\&delta;e}}_1\\ y_1=a\left({\mathrm{\&delta;M}}_1+{\mathrm{\&delta;\&omega;}}_1+\cos {\mathrm{i\&delta;\&Omega;}}_1\right)+2a\sin f_1{\mathrm{\&delta;e}}_1\\ z_1=a\sin \left(f_1+\mathrm{\&omega;}\right){\mathrm{\&delta;i}}_1-a\cos \left(f_1+\mathrm{\&omega;}\right)\sin {\mathrm{i\&delta;\&Omega;}}_1\end{array}\right.---\left(18\right)$

f₁＝M₁＝M₀+n(t₁-t₀) (19)

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δM₁、δω₁、δi₁、δΩ₁代入公式(18)～(19)，即可求得ρ₁值；

t₁和t₀的值作为已知值，通过公式(11)～(19)即可求得和的值；

S52，求解和值：

跟踪卫星在t₀时刻施加脉冲速度增量控制前相对参考星的速度；

跟踪卫星在t₁时刻施加脉冲速度增量控制后相对参考星的速度；

(1)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0={\&lsqb;{\stackrel{\&CenterDot;}{x}}_0,{\stackrel{\&CenterDot;}{y}}_0,{\stackrel{\&CenterDot;}{z}}_0\&rsqb;}^T---\left(20\right)$

其中：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_0=na\sin f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{y}}_0=-\frac{3}{2}{\mathrm{n\&delta;a}}_0+2na\cos f_0{\mathrm{\&delta;e}}_0\\ {\stackrel{\&CenterDot;}{z}}_0=na\&lsqb;\cos \left(\mathrm{\&omega;}+f_0\right){\mathrm{\&delta;i}}_0+\sin \left(\mathrm{\&omega;}+f_0\right)\sin {\mathrm{i\&delta;\&Omega;}}_0\&rsqb;\end{array}\right.---\left(21\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及初始轨道根数差δe₀、δa₀、δi₀、δΩ₀、以及通过公式17计算得到的f₀代入公式(20)～(21)，即可求得值；

(2)求解的值：

由于：

${\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1={\&lsqb;{\stackrel{\&CenterDot;}{x}}_1,{\stackrel{\&CenterDot;}{y}}_1,{\stackrel{\&CenterDot;}{z}}_1\&rsqb;}^T---\left(22\right)$

其中，

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_1=na\sin f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{y}}_1=-\frac{3}{2}{\mathrm{n\&delta;a}}_1+2na\cos f_1{\mathrm{\&delta;e}}_1\\ {\stackrel{\&CenterDot;}{z}}_1=na\&lsqb;\cos \left(\mathrm{\&omega;}+f_1\right){\mathrm{\&delta;i}}_1+\sin \left(\mathrm{\&omega;}+f_1\right)\sin {\mathrm{i\&delta;\&Omega;}}_1\&rsqb;\end{array}\right.---\left(23\right)$

将参考卫星的初始轨道根数n、a、i、ω和M₀，以及通过公式(5)求解得到的期望轨道根数差δe₁、δa₁、δi₁、δΩ₁、以及通过公式(19)计算得到的f₁代入公式(22)～(23)，即可求得值；

其中，S51和S52不分先后顺序；

S53，将和代入公式(24)，求得Δv₀和Δv₁：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_0={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_0^{-}\\ {\mathrm{\&Delta;v}}_1={\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{+}-{\stackrel{\&CenterDot;}{\mathrm{\&rho;}}}_1^{-}\end{array}\right.---\left(24\right)$

其中：

Δv₀代表初始时刻t₀需要施加给跟踪卫星的脉冲速度增量；

Δv₁代表终端时刻t₁需要施加给跟踪卫星的脉冲速度增量；

上述公式(24)即是所述的预给定的定量解析控制模型，该模型中所需要的中间参数由公式(11)～(23)计算得到；

S6，在初始时刻t₀，将S5计算得到的初始脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹；当达到终端时刻t₁时，将S5计算得到的终端脉冲速度增量施加给跟踪卫星，改变跟踪卫星相对参考卫星的运动轨迹，进而形成满足给定边界d_max和d_min的有界伴飞构型。

2.根据权利要求1所述的圆参考轨道下卫星有界伴飞的边界控制方法，其特征在于，终端时刻t₁通过以下方法确定：

如果伴飞控制任务要求形成伴飞构型的时间不超过T，则有界伴飞控制的终端时刻满足0＜t₁≤T，具体数值可通过优化方法搜索确定；

如果伴飞控制任务要求形成伴飞构型的时刻为T，则有界伴飞控制的终端时刻t₁＝T。