CN104309822B - 一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法 - Google Patents

Abstract

一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法，属于航天器轨道控制技术领域。本发明解决了现有的定点悬停方法要求控制量是连续的；现有的单脉冲水滴形绕飞方法来实现悬停，没有考虑由于悬停时间较长，悬停在目标航天器轨道平面的追踪航天器的燃料消耗的问题。技术方案为：目标航天器处于圆形高轨轨道，相对位置范围有上下边界，本发明采用带参数优化的单脉冲水滴形绕飞轨迹方案来实现，在基于hill方程的相对运动坐标系下考虑，只要在使整个水滴形轨迹都满足悬停的位置范围要求基础上，找到使性能指标值即燃料消耗最小的方案即可。本发明主要用于航天器的轨道控制。

Description

一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法

技术领域

本发明涉及航天器的近距离相对轨道运动控制方法，尤其涉及悬停控制方法，属于航天器轨道控制技术领域。

背景技术

当今航天领域的一个重要研究热点就是航天器的近距离相对轨道运动控制，航天器的相对轨道运动是研究一个航天器(追踪航天器)处于另一个航天器(目标航天器)周围的持续运动规律。常被应用到编队飞行、在轨维护、交会对接、跟踪监视等空间任务。目前最常用的相对轨道运动形式有悬停(追踪航天器与目标航天器保持相对位置始终不变)、伴随飞行(追踪航天器围绕目标航天器附近某点进行封闭轨迹飞行)和绕飞(伴随飞行的一种特殊情况，封闭轨迹的中心是目标航天器)等。

在实际工程中，由于航天器受到各种各样的摄动力的影响，运动规律变得非常复杂，真正意义上的“悬停”是不可能实现的。然而如果工程中只要求追踪航天器在一定时间内其位置限制在某个范围内，而相对速度不作过多的要求，这种广义上的悬停是可以实现的。

现有的定点悬停多用于航天器的交会对接中，或其他要求两航天器相对位置较严格地保持不变的场合。从hill方程或两体问题轨道动力学模型出发，可以获得实现定点悬停的开环控制方案，通过运用一些常用的控制方法，如PID控制、最优控制、变结构控制等，也可以获得相应的闭环控制方案，但无论是开环控制方案还是闭环控制方案，所施加的控制量都是连续的，连续控制要求发动机连续工作，造成燃料的持续消耗，而且由于姿态、轨道耦合控制问题的存在，影响航天器同时进行一些姿态指向的空间任务。

现有的单脉冲水滴形绕飞方法多用于局部限制轨迹，所谓局部限制轨迹，是指将相对运动轨迹限制在目标航天器附近的某个特定区域内。若用该方法来实现一定范围内的悬停，只需满足轨迹在空间范围约束内即可，不需要考虑燃料消耗的问题。

发明内容

本发明的目的是提出一种基于参数优化的单脉冲水滴形绕飞轨迹悬停控制方法，以解决针对现有的定点悬停方法要求控制量是连续的；现有的单脉冲水滴形绕飞方法来实现悬停，没有考虑由于悬停时间较长，悬停在目标航天器轨道平面的追踪航天器的燃料消耗的问题。

本发明为解决上述技术问题所采用的技术方案是：

本发明所述的一种基于参数优化的单脉冲水滴形绕飞轨迹悬停控制方法，是按照以下步骤实现的：

所述方法中的坐标系进行如下定义：在相对坐标系即轨道坐标系s-xyz中：坐标原点s与目标航天器的质心固连并随其沿轨道运动，x轴与目标航天器的地心矢量r_s重合，由地心指向s，y轴在目标航天器的轨道面内垂直于x轴，并指向运动方向，z轴与x轴，y轴组成直角右手坐标系；

其特征在于所述方法包括以下步骤：

步骤一、选取追踪航天器平均燃料消耗形式的性能指标，即如下目标函数，并使如下目标函数取得极小值

$J=\frac{\left|\mathrm{\Δ}v\right|}{\mathrm{\β}}---\left(5\right)$

其中，Δv——维持水滴形绕飞轨迹所需要施加的速度脉冲增量；

β——追踪航天器水滴形绕飞轨迹的周期T_w与目标航天器的轨道周期T的比值；

步骤二、选取使步骤一所述目标函数取得极小值的决策变量；

具体过程为：追踪航天器在水滴形绕飞轨迹中的x方向位置极限点A满足的绕飞条件如下：

$\{\begin{array}{c}{\stackrel{\·}{x}}_0=\frac{-3\mathrm{\π}\mathrm{\β}\sin \mathrm{\π}\mathrm{\β}}{4\sin \mathrm{\π}\mathrm{\β}-3\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}}{\mathrm{nx}}_0\\ {\stackrel{\·}{y}}_0=\frac{6(\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}-\sin \mathrm{\π}\mathrm{\β})}{4\sin \mathrm{\π}\mathrm{\β}-3\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}}{\mathrm{nx}}_0\end{array}---\left(3\right)$

其中，A点也为速度脉冲施加点，x₀——追踪航天器在A点的x轴方向的相对位置分量；

——追踪航天器在A点的相对线速度分量；n——目标航天器的平均轨道角速度；

每经过周期T_w，维持水滴形绕飞轨迹所需要施加的速度脉冲增量为

$\mathrm{\Δ}v=2{\stackrel{\·}{x}}_0---\left(4\right)$

结合式(3)、式(4)和式(5)，选取A点的相对位置分量x₀和β作为决策变量；

步骤三、求取使步骤一所述目标函数达到极小值需满足的约束条件；

具体过程为：当不施加主动控制时，xy平面内的状态转移方程为

$\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\end{array}\right]\left[\begin{array}{cccc}4-3\cos nt& 0& \frac{\sin nt}{n}& \frac{2(1-\cos nt)}{n}\\ 6(\sin nt-nt)& 1& \frac{2(\cos nt-1)}{n}& \frac{4\sin nt}{n}-3t\\ 3n\sin nt& 0& \cos nt& 2\sin nt\\ 6n(\cos nt-1)& 0& -2\sin nt& 4\cos nt-3\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ {\stackrel{\·}{x}}_0\\ {\stackrel{\·}{y}}_0\end{array}\right]---\left(39\right)$

其中，x,y——追踪航天器在相对坐标系中的位置分量；

——追踪航天器在相对坐标系中的线速度分量；

y₀——追踪航天器在A点的y轴方向的相对位置分量；

对于除A点之外的另一个x方向的位置极限点C，由式(39)可知

$\stackrel{\·}{x}=3{\mathrm{nx}}_0\sin nt+{\stackrel{\·}{x}}_0\cos nt+2{\stackrel{\·}{y}}_0\sin nt---\left(43\right)$

令式(43)的值为0，解得到达C点的时间t_C，将其代入式(39)，求得C点相对位置的分量x_C；同理，由(39)也可知

$\stackrel{\·}{y}=6{\mathrm{nx}}_0(\cos nt-1)-2{\stackrel{\·}{x}}_0\sin nt+{\stackrel{\·}{y}}_0(4\cos nt-3)---\left(44\right)$

令式(44)的值为0，解得到达y方向的位置极限点B,D的时间t_B,t_D，将其代入式(39)，求得B,D点相对位置的分量y_B,y_D；

y_B,x_C,y_D的表达式如下

$\{\begin{array}{c}{t}_C=\frac{1}{n}\arccos \left(\frac{3{\mathrm{nx}}_0+2{\stackrel{\·}{y}}_0}{\sqrt{{(3{\mathrm{nx}}_0+2{\stackrel{\·}{y}}_0)}^2+{\left({\stackrel{\·}{x}}_0\right)}^2}}\right)\\ x_C=(4-3\cos {\mathrm{nt}}_C)x_0+\frac{\sin {\mathrm{nt}}_C}{n}{\stackrel{\·}{x}}_0+\frac{2(1-\cos {\mathrm{nt}}_C)}{n}{\stackrel{\·}{y}}_0\\ t_B=\frac{1}{n}\[\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)-\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)\]\\ t_D=\frac{1}{n}\[\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)+\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)\]\\ y_B=6(\sin {\mathrm{nt}}_B-{\mathrm{nt}}_B)x_0+y_0+\frac{2(\cos {\mathrm{nt}}_B-1)}{n}{\stackrel{\·}{x}}_0+(\frac{4\sin {\mathrm{nt}}_B}{n}-3t_B){\stackrel{\·}{y}}_0\\ y_D=6(\sin {\mathrm{nt}}_D-{\mathrm{nt}}_D)x_0+y_0+\frac{2(\cos {\mathrm{nt}}_D-1)}{n}{\stackrel{\·}{x}}_0+(\frac{4\sin {\mathrm{nt}}_D}{n}-3t_D){\stackrel{\·}{y}}_0\end{array}---\left(45\right)$

综合考虑决策变量的约束条件，所有的约束条件如下：

$\left\{\begin{array}{c}{x}_{\min }\≤x_0\≤x_{max}\\ 0\≤\mathrm{\β}\≤1\\ x_{\min }\≤x_C\≤x_{max}\\ y_{\min }\≤y_B\≤y_{max}\\ y_{\min }\≤y_D\≤y_{max}\end{array}\right.---\left(7\right)$

其中，x_min,x_max,y_min,y_max——相对位置范围在x,y方向的上下边界分量；

步骤四、寻找使步骤一所述目标函数取极小值时的x₀和β；

步骤五、用步骤四得到的x₀和β，实施追踪航天器周期性地沿水滴形轨迹飞行的控制；

具体过程为：悬停任务开始时，当追踪航天器到达相对坐标系中的A点，施加速度脉冲，使速度满足绕飞条件，追踪航天器开始沿着水滴形轨迹飞行，经过一个周期T_w后回到A点，此时再沿x方向施加维持水滴形绕飞轨迹所需要的速度脉冲增量Δv，重新使速度满足绕飞条件，从而使追踪航天器周期性地沿水滴形轨迹飞行。

本发明的有益效果是：

1、与现有的定点悬停相比，现有的定点悬停要求控制量是连续的，本发明方法是不需要连续的，用的是脉冲控制，而且不存在姿态、轨道耦合控制问题，不影响航天器同时进行一些姿态指向的空间任务。

2、本发明将现有的单脉冲水滴形绕飞方法应用在悬停方面，并考虑到悬停在目标航天器轨道平面的追踪航天器的燃料消耗，因此本发明在其基础上进行了改进，通过对决定水滴形轨迹的两个因素——A点位置的x方向分量和比例系数β进行优化，并考虑了航天器结构和工程实际，从而在满足空间悬停范围约束的基础上实现燃料消耗的最优。

附图说明

图1为本发明的流程图；

图2为具体实施方式一和具体实施方式二中所述的水滴形绕飞轨迹示意图，其中方向箭头表示追踪航天器的运行方向，x轴和y轴表示追踪航天器在相对坐标系中的位置分量；

图3为仿真验证中目标航天器悬停的位置范围示意图；

图4为仿真验证中水滴形轨迹周期T_w与β的关系曲线图；

图5为仿真验证中A点x₀坐标与β的关系曲线图；

图6为仿真验证中维持水滴形轨迹所需单次速度脉冲增量Δv与β的关系曲线图；

图7为仿真验证中性能指标值J与β的关系曲线图；

图8为仿真验证中水滴形轨迹与允许的位置范围示意图；

图9为仿真验证中Simulink仿真模型图；

图10为图9的仿真结果中水滴形轨迹示意图；

图11为具体实施方式一中地心惯性系与轨道坐标系示意图，其中r表示位置矢量，h表示轨道平面法向量，N表示轨道平面与赤道平面交线的方向矢量；

图12为具体实施方式一中轨道平面示意图，其中r表示位置矢量，v表示速度矢量，e_r表示位置矢量方向的单位矢量，e_f表示轨道平面内，垂直于位置矢量方向的单位矢量；

图13为具体实施方式一中相对运动坐标系与地心惯性坐标系的关系示意图；

图14为本发明方法的示意图。

具体实施方式

结合附图进一步详细说明本发明的具体实施方式。

具体实施方式一：结合图1、图2、图14说明本实施方式。本实施方式所述的一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法，对所述方法中的坐标系进行如下定义：在相对坐标系即轨道坐标系s-xyz中：坐标原点s与目标航天器的质心固连并随其沿轨道运动，x轴与目标航天器的地心矢量r_s重合，由地心指向s，y轴在目标航天器的轨道面内垂直于x轴，并指向运动方向，z轴与x轴，y轴组成直角右手坐标系；

其特征在于所述方法包括以下步骤：

步骤一、为了考虑悬停时间内追踪航天器燃料消耗情况，选取追踪航天器平均燃料消耗形式的性能指标，即如下目标函数，并使如下目标函数取得极小值

$J=\frac{\left|\mathrm{\Δ}v\right|}{\mathrm{\β}}---\left(5\right)$

其中，Δv——维持水滴形绕飞轨迹所需要施加的速度脉冲增量；

β——追踪航天器水滴形绕飞轨迹的周期T_w与目标航天器的轨道周期T的比值；

步骤二、选取使步骤一所述目标函数取得极小值的决策变量；

具体过程为：因为A点在轨道平面xy内的位置[x₀,y₀]，决定了水滴形绕飞轨迹脉冲施加点的位置，而追踪航天器水滴形绕飞轨迹的周期T_w决定了整个轨迹的大小，两者结合则可以唯一地确定水滴形绕飞轨迹。追踪航天器在水滴形绕飞轨迹中的x方向位置极限点A满足的绕飞条件如下：

$\{\begin{array}{c}{\stackrel{\·}{x}}_0=\frac{-3\mathrm{\π}\mathrm{\β}\sin \mathrm{\π}\mathrm{\β}}{4\sin \mathrm{\π}\mathrm{\β}-3\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}}{\mathrm{nx}}_0\\ {\stackrel{\·}{y}}_0=\frac{6(\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}-\sin \mathrm{\π}\mathrm{\β})}{4\sin \mathrm{\π}\mathrm{\β}-3\mathrm{\π}\mathrm{\β}\cos \mathrm{\π}\mathrm{\β}}{\mathrm{nx}}_0\end{array}---\left(3\right)$

其中，A点也为速度脉冲施加点，x₀——追踪航天器在A点的x轴方向的相对位置分量；——追踪航天器在A点的相对线速度分量；n——目标航天器的平均轨道角速度；每经过周期T_w，维持水滴形绕飞轨迹所需要施加的速度脉冲增量为

$\mathrm{\Δ}v=2{\stackrel{\·}{x}}_0---\left(4\right)$

结合式(3)、式(4)和式(5)，得出平均燃料消耗的多少与A点的相对位置分量y₀的选取无关，因此y₀可取允许范围中的任意值。又因为追踪航天器水滴形绕飞轨迹的周期T_w＝βT，所以选取A点的相对位置分量x₀和β作为决策变量；

步骤三、求取使步骤一所述目标函数达到极小值需满足的约束条件；

具体过程为：采用带参数优化的单脉冲水滴形绕飞轨迹来实现脉冲控制下的悬停。当不施加主动控制时，xy平面内的状态转移方程为

$\left[\begin{array}{c}x\\ y\\ \stackrel{\·}{x}\\ \stackrel{\·}{y}\end{array}\right]\left[\begin{array}{cccc}4-3\cos nt& 0& \frac{\sin nt}{n}& \frac{2(1-\cos nt)}{n}\\ 6(\sin nt-nt)& 1& \frac{2(\cos nt-1)}{n}& \frac{4\sin nt}{n}-3t\\ 3n\sin nt& 0& \cos nt& 2\sin nt\\ 6n(\cos nt-1)& 0& -2\sin nt& 4\cos nt-3\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ {\stackrel{\·}{x}}_0\\ {\stackrel{\·}{y}}_0\end{array}\right]---\left(39\right)$

其中，x,y——追踪航天器在相对坐标系中的位置分量；

——追踪航天器在相对坐标系中的线速度分量；

y₀——追踪航天器在A点的y轴方向的相对位置分量；

要使水滴形绕飞轨迹能满足允许的悬停位置范围，只需使水滴形轨迹的4个位置极限点A,B,C,D在允许范围内即可。

对于除A点之外的另一个x方向的位置极限点C，由式(39)可知

$\stackrel{\·}{x}=3{\mathrm{nx}}_0\sin nt+{\stackrel{\·}{x}}_0\cos nt+2{\stackrel{\·}{y}}_0\sin nt---\left(43\right)$

令式(43)的值为0，解得到达C点的时间t_C，将其代入式(39)，求得C点相对位置的分量x_C；同理，由(39)也可知

$\stackrel{\·}{y}=6{\mathrm{nx}}_0(\cos nt-1)-2{\stackrel{\·}{x}}_0\sin nt+{\stackrel{\·}{y}}_0(4\cos nt-3)---\left(44\right)$

令式(44)的值为0，解得到达y方向的位置极限点B,D的时间t_B,t_D，将其代入式(39)，求得B,D点相对位置的分量y_B,y_D；

只要约束x_C,y_B,y_D满足允许的位置范围，则整个水滴形轨迹都在允许范围内,y_B,x_C,y_D的表达式如下

$\{\begin{array}{c}{t}_C=\frac{1}{n}\arccos \left(\frac{3{\mathrm{nx}}_0+2{\stackrel{\·}{y}}_0}{\sqrt{{(3{\mathrm{nx}}_0+2{\stackrel{\·}{y}}_0)}^2+{\left({\stackrel{\·}{x}}_0\right)}^2}}\right)\\ x_C=(4-3\cos {\mathrm{nt}}_C)x_0+\frac{\sin {\mathrm{nt}}_C}{n}{\stackrel{\·}{x}}_0+\frac{2(1-\cos {\mathrm{nt}}_C)}{n}{\stackrel{\·}{y}}_0\\ t_B=\frac{1}{n}\[\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)-\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)\]\\ t_D=\frac{1}{n}\[\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)+\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\·}{y}}_0}{\sqrt{{(6{\mathrm{nx}}_0+4{\stackrel{\·}{y}}_0)}^2+{\left(2{\stackrel{\·}{x}}_0\right)}^2}}\right)\]\\ y_B=6(\sin {\mathrm{nt}}_B-{\mathrm{nt}}_B)x_0+y_0+\frac{2(\cos {\mathrm{nt}}_B-1)}{n}{\stackrel{\·}{x}}_0+(\frac{4\sin {\mathrm{nt}}_B}{n}-3t_B){\stackrel{\·}{y}}_0\\ y_D=6(\sin {\mathrm{nt}}_D-{\mathrm{nt}}_D)x_0+y_0+\frac{2(\cos {\mathrm{nt}}_D-1)}{n}{\stackrel{\·}{x}}_0+(\frac{4\sin {\mathrm{nt}}_D}{n}-3t_D){\stackrel{\·}{y}}_0\end{array}---\left(45\right)$

综合考虑决策变量的约束条件，所有的约束条件如下：

$\left\{\begin{array}{c}{x}_{\min }\≤x_0\≤x_{max}\\ 0\≤\mathrm{\β}\≤1\\ x_{\min }\≤x_C\≤x_{max}\\ y_{\min }\≤y_B\≤y_{max}\\ y_{\min }\≤y_D\≤y_{max}\end{array}\right.---\left(7\right)$

其中，x_min,x_max,y_min,y_max——相对位置范围在x,y方向的上下边界分量；

步骤四、寻找使步骤一所述目标函数取极小值时的x₀和β；

步骤五、用步骤四得到的x₀和β，实施追踪航天器周期性地沿水滴形轨迹飞行的控制；

通过步骤四找到的x₀和β，得到A点的位置、水滴形轨迹的周期T_w以及维持水滴形绕飞轨迹所需要施加的速度脉冲增量Δv；

具体过程为：悬停任务开始时，当追踪航天器到达相对坐标系中的A点，施加速度脉冲，使速度满足绕飞条件，追踪航天器开始沿着水滴形轨迹飞行，经过一个周期T_w后回到A点，此时再沿x方向施加维持水滴形绕飞轨迹所需要的速度脉冲增量Δv，重新使速度满足绕飞条件，从而使追踪航天器周期性地沿水滴形轨迹飞行。

本实施方式中的相关技术介绍如下：

结合图11所示的地心惯性坐标系(O-XYZ)：坐标原点O在地球质心，Z轴向北指向平赤道面北极，X轴指向平春分点，Y轴与Z轴、X轴组成直角右手坐标系。

轨道坐标系(s-xyz)：坐标原点s与航天器的质心固连并随其沿轨道运动，x轴与航天器的地心矢量r_s重合，由地心指向s，y轴在航天器的轨道面内垂直于x轴，并指向运动方向，z轴与x轴，y轴组成直角右手坐标系。

开普勒方程：两体运动问题中，对于椭圆轨道，偏心率为e，偏近点角为E，平近点角为M，沿逆时针方向为正，开普勒方程可以表示为

E-esin(E)＝M(1)

Hill方程：假定两航天器仅受地球的引力作用，以目标航天器的轨道坐标系作为相对运动坐标系，通过一阶线性化，将相对运动动力学方程化为一组常系数线性微分方程

$\left\{\begin{array}{c}\stackrel{\·\·}{x}-2n\stackrel{\·}{y}-3n^{2}x=a_x\\ \stackrel{\·\·}{y}+2n\stackrel{\·}{x}=a_y\\ \stackrel{\·\·}{z}+n^{2}z=a_z\end{array}\right.---\left(2\right)$

脉冲控制：发动机短时工作情况下的轨道控制，由于发动机工作的时间比轨道飞行周期短得多，因而可以视为是脉冲作用。

单脉冲水滴形绕飞轨迹：假设追踪航天器初始时刻位于相对运动坐标系中的某点，当初始相对速度满足一定条件时，追踪航天器将沿着水滴形轨迹飞行，并最终返回到该点，施加一个速度脉冲增量Δv使其重新满足绕飞的条件，从而使追踪航天器周期性地沿水滴形轨迹飞行。

非线性规划问题：研究多元实函数在有非线性约束条件下的极值问题，包含决策变量、目标函数和约束条件三个部分，其中决策变量即待优化的参数，目标函数即问题所关心的目标与决策变量间的函数关系，约束条件包括线性或非线性的等式约束或不等式约束。

本实施方式中涉及的相关理论介绍如下：

一、两体问题轨道模型：

Kepler轨道及轨道要素——航天器轨道是指航天器在地心惯性系下的位置，描述的是航天器的平动运动。作为二体问题，两个天体均作为质点对待。考虑质量分别为m和M的两个天体间运动方程为：

$\stackrel{\·\·}{r}=-\frac{G(M+m)}{r^3}r=-\frac{\mathrm{\μ}}{r^3}r---\left(8\right)$

通过对二体问题的求解，天体(航天器)在惯性空间中的运动可以用六个经典的轨道要素(也称轨道根数，如图11所示)来描述：

a:轨道半长轴；

e:偏心率；

Ω:升交点赤经，在J2000地心惯性系中自Ox轴方向在xy平面内逆时针度量到升交点的角度，0≤Ω<2π；

i:轨道倾角，轨道正法向h与J2000地心惯性系Oz轴的夹角，0≤i≤π。若0≤i<π/2，为顺行轨道，航天器偏向东飞行；若π/2≤i<π，为逆行轨道，航天器偏向西飞行；若i＝π/2，为极轨道；

ω:近地点幅角，在轨道平面内自轨道升交点方向到偏心率矢量e之间的夹角，顺航天器方向度量，0≤ω<2π；

f:真近点角，在轨道平面内从e到r之间的夹角。

轨道计算与星历表计算——给定初始条件，从观测的航天器位置和速度数据出发，可以计算轨道根数，这一过程称为轨道计算。同时，轨道根数给定后，也可以计算出任意时刻航天器的位置矢量和速度矢量，该过程称为星历表计算。

以下简单介绍一下两种过程的计算方法：

1、轨道计算

求半长轴a：

$a=\frac{\mathrm{\&mu;}r}{2\mathrm{\&mu;}-{\mathrm{rv}}^2}---\left(9\right)$

其中v为航天器轨道速度。

求偏心率e：

$e=\sqrt{1-p/a}---\left(10\right)$

其中p为椭圆半通径。

轨道倾角i与升交点赤经Ω计算如下：

i＝arccos(h_z/h)(11)

Ω＝-arctan(2h_x/h_y)(12)

求真近点角f：

cosf＝(p/r-1)/e(13)

求近地点幅角ω：

令i_Ω为地心指向升交点的单位矢量，在地心惯性系中表示为i_Ω＝cosΩ·x_i+sinΩ·y_i。故可由下式确定轨道角u

$\cos u=\frac{1}{r}{r}^T{i}_{\mathrm{\&Omega;}}---\left(14\right)$

进而可得近地点幅角为

ω＝u-f(15)

注：当e＝0时，会出现奇异。因此一般考虑当e<10^-6时，轨道称为圆轨道，有ω＝0，f＝u。

2、星历表计算

在轨道平面内定义如图11所示的坐标系，有

$\begin{array}{c}r=a(\cos E-e)\&CenterDot;P+a\sqrt{1-e^2}\sin E\&CenterDot;Q\\ v=-\sqrt{\mathrm{\&mu;}/p}sinf\&CenterDot;P+\sqrt{\mathrm{\&mu;}/p}(e+\cos f)\&CenterDot;Q\end{array}---\left(16\right)$

其中，P和Q分别表示卫星近地点和半通径方向的单位矢量。通过3-1-3的欧拉转序旋转，P和Q的表达式如下所示：

$P=R_3(-\mathrm{\&Omega;})R_1(-i)R_3(-\mathrm{\&omega;})\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}cos\mathrm{\&Omega;}cos\mathrm{\&omega;}-sin\mathrm{\&Omega;}sin\mathrm{\&omega;}\cos i\\ sin\mathrm{\&Omega;}cos\mathrm{\&omega;}+cos\mathrm{\&Omega;}sin\mathrm{\&omega;}\cos i\\ sin\mathrm{\&omega;}\sin i\end{array}\right]---\left(17\right)$

$Q=R_3(-\mathrm{\&Omega;})R_1(-i)R_3(-\mathrm{\&omega;})\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]=\left[\begin{array}{c}-cos\mathrm{\&Omega;}sin\mathrm{\&omega;}-sin\mathrm{\&Omega;}cos\mathrm{\&omega;}\cos i\\ -sin\mathrm{\&Omega;}sin\mathrm{\&omega;}+cos\mathrm{\&Omega;}cos\mathrm{\&omega;}\cos i\\ cos\mathrm{\&omega;}\sin i\end{array}\right]---\left(18\right)$

二、相对轨道动力学模型

记目标航天器为s，追踪航天器为c。取目标航天器的轨道坐标系s-xyz作为相对运动坐标系，其原点与目标航天器的质心固连并随其沿轨道运动，x轴与目标航天器的地心矢量r_s重合，由地心指向s，y轴于目标航天器的轨道面内垂直于x轴，并指向运动方向，z轴由右手规则确定，亦即z轴与目标航天器轨道动量矩矢量的方向一致。轨道坐标系s-xyz与地心惯性坐标系O_E-XYZ的关系如图13所示。

在轨道坐标系中有

$r_s=\left[\begin{array}{c}{r}_s\\ 0\\ 0\end{array}\right]---\left(19\right)$

设追踪航天器的地心位置为r_c，则其对于目标航天器的位置矢量ρ为

$\mathrm{\&rho;}=r_c-r_S=\left[\begin{array}{c}x\\ y\\ z\end{array}\right]---\left(20\right)$

在地心惯性坐标中，目标航天器和追踪航天器的动力学方程如下

$\frac{d^2{r}_s}{{\mathrm{dt}}^2}=-\frac{{\mathrm{\&mu;r}}_s}{{r_s}^3}+a_s---\left(21\right)$

$\frac{d^2{r}_c}{{\mathrm{dt}}^2}=-\frac{{\mathrm{\&mu;r}}_c}{{r_s}^3}+a_c---\left(22\right)$

其中a_s和a_c分别为目标航天器和追踪航天器除地球中心重力以外的其他作用力的合力的加速度矢量，即对推力和摄动力(包括地球形状摄动，大气阻力摄动和太阳光压摄动)的加速度矢量。

由式(20)、式(21)和式(22)可得追踪航天器与目标航天器的绝对加速度之差为

$\frac{d^2\mathrm{\&rho;}}{{\mathrm{dt}}^2}=\frac{d^2{r}_c}{{\mathrm{dt}}^2}-\frac{d^2{r}_s}{{\mathrm{dt}}^2}=-\frac{{\mathrm{\&mu;r}}_c}{{r_c}^3}+\frac{{\mathrm{\&mu;r}}_s}{{r_s}^3}+a_c-a_s---\left(23\right)$

上式可进一步表示为下列等效的关系式

$\frac{d^2\mathrm{\&rho;}}{{\mathrm{dt}}^2}=\frac{\mathrm{\&mu;}}{{r_s}^3}\&lsqb;r_s-{\left(\frac{r_s}{r_c}\right)}^3{r}_c\&rsqb;+\mathrm{\&Delta;}a---\left(24\right)$

为了建立追踪航天器与目标航天器在动坐标系s-xyz中的相对运动方程，有

$\frac{d^2\mathrm{\&rho;}}{{\mathrm{dt}}^2}=\frac{{\mathrm{\&delta;}}^2\mathrm{\&rho;}}{{\mathrm{\&delta;t}}^2}+2\stackrel{\&CenterDot;}{\mathrm{\&theta;}}\&times;v+\stackrel{\&CenterDot;}{\mathrm{\&theta;}}\&times;(\stackrel{\&CenterDot;}{\mathrm{\&theta;}}\&times;\mathrm{\&rho;})+\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}\&times;\mathrm{\&rho;}---\left(25\right)$

上式中的和v分别为追踪航天器在目标航天器轨道坐标系中的相对加速度矢量和相对速度矢量，则

$\frac{{\mathrm{\&delta;}}^2\mathrm{\&rho;}}{{\mathrm{\&delta;t}}^2}=\left[\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{x}\\ \stackrel{\&CenterDot;\&CenterDot;}{y}\\ \stackrel{\&CenterDot;\&CenterDot;}{z}\end{array}\right]---\left(26\right)$

$v=\left[\begin{array}{c}\stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\\ \stackrel{\&CenterDot;}{z}\end{array}\right]---\left(27\right)$

式(25)中的分别为轨道参考坐标系旋转的角加速度矢量和角速度矢量。目标航天器的平均运动角速度为

由轨道要素描述的卫星位置和轨道角速度的形式如下

$\left|r_s\right|=\frac{a(1-e^2)}{1+ecos\mathrm{\&theta;}}---\left(28\right)$

$\stackrel{\&CenterDot;}{\mathrm{\&theta;}}=\frac{n{(1+ecos\mathrm{\&theta;})}^2}{{(1-e^2)}^{\frac{3}{2}}}---\left(29\right)$

针对圆轨道e＝0，由式(29)得并可得下列近似式

$\left\{\begin{array}{c}\stackrel{\&CenterDot;}{n}=0\\ n=\left[\begin{array}{c}0\\ 0\\ n\end{array}\right]\end{array}\right.---\left(30\right)$

对于追踪航天器和目标航天器的近距离的相对运动情况，卫星间距ρ是小量，特别是相对高度不大(即r_s/r_c接近1)时，从式(24)可以看出，其右端第一项为中心引力加速度差表达式，通过取一次近似(即线性化)可以进行简化。下面为具体的简化过程，可以看出简化实质是对中心引力取一次近似(即线性化)。

因为

$r_c={\&lsqb;{(x+r_s)}^2+y^2+z^2\&rsqb;}^{\frac{1}{2}}={({\mathrm{\&rho;}}^2+{r_s}^2+2{\mathrm{xr}}_s)}^{\frac{1}{2}}---\left(31\right)$

所以有

${\left(\frac{r_s}{r_c}\right)}^3={\&lsqb;1+{\left(\frac{\mathrm{\&rho;}}{r_s}\right)}^2+2\frac{x}{r_S}\&rsqb;}^{-\frac{3}{2}}---\left(32\right)$

在上式忽略以及更高的幂次项，则近似可得一下公式

${\left(\frac{r_s}{r_c}\right)}^3\&ap;1-3\frac{x}{r_s}---\left(33\right)$

将上式(33)和式(20)代入式(24)得

$\frac{d^2\mathrm{\&rho;}}{{\mathrm{dt}}^2}=\frac{\mathrm{\&mu;}}{{r_s}^3}\&lsqb;r_s-(1-3\frac{x}{r_c})(r_s+\mathrm{\&rho;})\&rsqb;+\mathrm{\&Delta;}a---\left(34\right)$

进一步简化上式，忽略小量xρ项，得

$\begin{array}{c}\frac{d^2\mathrm{\&rho;}}{{\mathrm{dt}}^2}=-\frac{\mathrm{\&mu;}}{{r_s}^3}\&lsqb;\mathrm{\&rho;}-\frac{3x}{r_s}{r}_s\&rsqb;+\mathrm{\&Delta;}a\\ =-n^2\left[\begin{array}{c}-2x\\ y\\ z\end{array}\right]+\mathrm{\&Delta;}a\end{array}---\left(35\right)$

将式(35)、式(26)、式(27)和式(30)代入式(25)可得

$\left\{\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{x}-2n\stackrel{\&CenterDot;}{y}-3n^{2}x=a_x\\ \stackrel{\&CenterDot;\&CenterDot;}{y}+2n\stackrel{\&CenterDot;}{x}=a_y\\ \stackrel{\&CenterDot;\&CenterDot;}{z}+n^{2}z=a_z\end{array}\right.---\left(36\right)$

通过简化，将相对运动动力学方程化为一组常系数线性微分方程。式(36)称为hill方程，也称Clohessey-Whiltshire方程(简称C-W方程)。

三、fmincon函数

Matlab优化工具箱中的fmincon函数可用来求解多元有非线性约束的极小值问题，其约束条件如下：

$\left\{\begin{array}{c}c\left(X\right)\&le;0\\ ceq\left(X\right)=0\\ A\&CenterDot;X\&le;b\\ Aeq\&CenterDot;X=beq\\ lb\&le;X\&le;ub\end{array}\right.---\left(37\right)$

常用表达式如下：

X＝fmincon(fun,X₀,A,b,Aeq,beq,lb,ub,nonlcon)(38)

其中：fun——表示目标函数；

X₀——表示选取的决策变量的初值；

A,b——表示约束条件中的线性不等式约束；

Aeq,beq——表示约束条件中的线性等式约束；

lb,ub——表示决策变量的下界和上界；

nonlcon——表示约束条件中的非线性约束。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤四的具体过程为:

不断地选取不同的x₀和β，采用非线性规划的方法，使用Matlab中的函数fmincon对A点位置的x方向分量和追踪航天器水滴形绕飞轨迹的周期T_w与目标航天器的轨道周期T的比值β进行优化，得到追踪航天器悬停所需燃料消耗最少的水滴形绕飞轨迹。其它步骤与具体实施方式一相同。

仿真验证

以目标航天器处于地球静止轨道(GEO)，悬停时间为一个轨道周期，即24小时，[x_min,x_max,y_min,y_max]＝[-60000,-40000,-10000,10000]m的正方形悬停位置范围为例，如图3所示。

选取y₀＝0，优化初值选取X₀＝[-50000,0.5]，优化后的结果如表1所示

表1

此时，水滴形轨迹的周期T_w≈3.7s，对于航天器结构和工程实际而言，这个周期过小。为了考虑工程实际中相邻两次速度脉冲的时间间隔约束，再进一步探求不同大小的水滴形轨迹(β不同)的优化结果，水滴形轨迹周期T_w与β的关系曲线如图4所示，A点x₀坐标与β的关系曲线如图5所示，维持水滴形轨迹所需单次速度脉冲增量Δv与β的关系曲线如图6所示，性能指标值J与β的关系曲线如图7所示。

由图5可以看出，当β约超过0.168后，水滴形轨迹将不能完全包含在允许的位置范围内，由图7可以看出，随着β的增大，性能指标值也逐渐增大，因此燃料消耗与水滴形轨迹周期直接存在矛盾，需根据实际情况对β值进行合理选择。

以β＝0.08为例，经过优化，脉冲施加点A的x方向分量x₀＝-4.393521106517433×10⁴m，理论上应维持的水滴形轨迹与允许的位置范围如图8所示。

在Simulink中进行24小时的仿真，如图9所示，仿真结果如图10所示，轨迹没有完全重叠是因为Simulink的仿真计算有数据位的限制，但基本可以证明绕飞条件式(40)及维持轨迹所需的速度脉冲增量式(41)的正确性。

Claims (2)

1.一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法，对所述方法中的坐标系进行如下定义：在相对坐标系即轨道坐标系s-xyz中：坐标原点s与目标航天器的质心固连并随其沿轨道运动，x轴与目标航天器的地心矢量r_s重合，由地心指向s，y轴在目标航天器的轨道面内垂直于x轴，并指向运动方向，z轴与x轴，y轴组成直角右手坐标系；

其特征在于所述方法包括以下步骤：

步骤一、选取追踪航天器平均燃料消耗形式的性能指标，即如下目标函数，并使如下目标函数取得极小值

$J=\frac{\left|\mathrm{\&Delta;}v\right|}{\mathrm{\&beta;}}---\left(5\right)$

其中，△v——维持水滴形绕飞轨迹所需要施加的速度脉冲增量；

β——追踪航天器水滴形绕飞轨迹的周期T_w与目标航天器的轨道周期T的比值；

步骤二、选取使步骤一所述目标函数取得极小值的决策变量；

具体过程为：追踪航天器在水滴形绕飞轨迹中的x方向位置极限点A满足的绕飞条件如下：

$\left\{\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_0=\frac{-3\mathrm{\&pi;}\mathrm{\&beta;}\sin \mathrm{\&pi;}\mathrm{\&beta;}}{4\sin \mathrm{\&pi;}\mathrm{\&beta;}-3\mathrm{\&pi;}\mathrm{\&beta;}\cos \mathrm{\&pi;}\mathrm{\&beta;}}{\mathrm{nx}}_0\\ {\stackrel{\&CenterDot;}{y}}_0=\frac{6\left(\mathrm{\&pi;}\mathrm{\&beta;}\cos \mathrm{\&pi;}\mathrm{\&beta;}-\sin \mathrm{\&pi;}\mathrm{\&beta;}\right)}{4\sin \mathrm{\&pi;}\mathrm{\&beta;}-3\mathrm{\&pi;}\mathrm{\&beta;}\cos \mathrm{\&pi;}\mathrm{\&beta;}}{\mathrm{nx}}_0\end{array}\right.---\left(3\right)$

其中，A点也为速度脉冲施加点，x₀——追踪航天器在A点的x轴方向的相对位置分量；——追踪航天器在A点的相对线速度分量；n——目标航天器的平均轨道角速度；每经过周期T_w，维持水滴形绕飞轨迹所需要施加的速度脉冲增量为

$\mathrm{\&Delta;}v=2{\stackrel{\&CenterDot;}{x}}_0---\left(4\right)$

结合式(3)、式(4)和式(5)，选取A点的相对位置分量x₀和β作为决策变量；

步骤三、求取使步骤一所述目标函数达到极小值需满足的约束条件；

具体过程为：当不施加主动控制时，xy平面内的状态转移方程为

$\left[\begin{array}{c}x\\ y\\ \stackrel{\&CenterDot;}{x}\\ \stackrel{\&CenterDot;}{y}\end{array}\right]=\left[\begin{array}{cccc}4-3\cos nt& 0& \frac{\sin nt}{n}& \frac{2\left(1-\cos nt\right)}{n}\\ 6\left(\sin nt-nt\right)& 1& \frac{2\left(\cos nt-1\right)}{n}& \frac{4\sin nt}{n}-3t\\ 3n\sin nt& 0& \cos nt& 2\sin nt\\ 6n\left(\cos nt-1\right)& 0& -2\sin nt& 4\cos nt-3\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\\ {\stackrel{\&CenterDot;}{x}}_0\\ {\stackrel{\&CenterDot;}{y}}_0\end{array}\right]---\left(39\right)$

其中，x,y——追踪航天器在相对坐标系中的位置分量；

——追踪航天器在相对坐标系中的线速度分量；

y₀——追踪航天器在A点的y轴方向的相对位置分量；

对于除A点之外的另一个x方向的位置极限点C，由式(39)可知

$\stackrel{\&CenterDot;}{x}=3{\mathrm{nx}}_0\sin nt+{\stackrel{\&CenterDot;}{x}}_0\cos nt+2{\stackrel{\&CenterDot;}{y}}_0\sin nt---\left(43\right)$

令式(43)的值为0，解得到达C点的时间t_C，将其代入式(39)，求得C点相对位置的分量x_C；同理，由(39)也可知

$\stackrel{\&CenterDot;}{y}=6{\mathrm{nx}}_0(\cos nt-1)-2{\stackrel{\&CenterDot;}{x}}_0\sin nt+{\stackrel{\&CenterDot;}{y}}_0(4\cos nt-3)---\left(44\right)$

令式(44)的值为0，解得到达y方向的位置极限点B,D的时间t_B,t_D，将其代入式(39)，求得B,D点相对位置的分量y_B,y_D；

y_B,x_C,y_D的表达式如下

$\left\{\begin{array}{c}{t}_C=\frac{1}{n}\arccos \left(\frac{3{\mathrm{nx}}_0+2{\stackrel{\&CenterDot;}{y}}_0}{\sqrt{{\left(3{\mathrm{nx}}_0+2{\stackrel{\&CenterDot;}{y}}_0\right)}^2+{\left({\stackrel{\&CenterDot;}{x}}_0\right)}^2}}\right)\\ x_C=\left(4-3\cos {\mathrm{nt}}_C\right)x_0+\frac{\sin {\mathrm{nt}}_C}{n}{\stackrel{\&CenterDot;}{x}}_0+\frac{2\left(1-\cos {\mathrm{nt}}_C\right)}{n}{\stackrel{\&CenterDot;}{y}}_0\\ t_B=\frac{1}{n}\&lsqb;\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0}{\sqrt{{\left(6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0\right)}^2+{\left(2{\stackrel{\&CenterDot;}{x}}_0\right)}^2}}\right)-\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\&CenterDot;}{y}}_0}{\sqrt{{\left(6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0\right)}^2+{\left(2{\stackrel{\&CenterDot;}{x}}_0\right)}^2}}\right)\&rsqb;\\ t_D=\frac{1}{n}\&lsqb;\arccos \left(\frac{6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0}{\sqrt{{\left(6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0\right)}^2+{\left(2{\stackrel{\&CenterDot;}{x}}_0\right)}^2}}\right)+\arccos \left(\frac{6{\mathrm{nx}}_0+3{\stackrel{\&CenterDot;}{y}}_0}{\sqrt{{\left(6{\mathrm{nx}}_0+4{\stackrel{\&CenterDot;}{y}}_0\right)}^2+{\left(2{\stackrel{\&CenterDot;}{x}}_0\right)}^2}}\right)\&rsqb;\\ y_B=6\left(\sin {\mathrm{nt}}_B-{\mathrm{nt}}_B\right)x_0+y_0+\frac{2\left(\cos {\mathrm{nt}}_B-1\right)}{n}{\stackrel{\&CenterDot;}{x}}_0+\left(\frac{4\sin {\mathrm{nt}}_B}{n}-3t_B\right){\stackrel{\&CenterDot;}{y}}_0\\ y_D=6\left(\sin {\mathrm{nt}}_D-{\mathrm{nt}}_D\right)x_0+y_0+\frac{2\left(\cos {\mathrm{nt}}_D-1\right)}{n}{\stackrel{\&CenterDot;}{x}}_0+\left(\frac{4\sin {\mathrm{nt}}_D}{n}-3t_D\right){\stackrel{\&CenterDot;}{y}}_0\end{array}\right.---\left(45\right)$

综合考虑决策变量的约束条件，所有的约束条件如下：

$\left\{\begin{array}{c}{x}_{min}\&le;x_0\&le;x_{max}\\ 0\&le;\mathrm{\&beta;}\&le;1\\ x_{min}\&le;x_C\&le;x_{max}\\ y_{min}\&le;y_B\&le;y_{max}\\ y_{\min }\&le;y_D\&le;y_{max}\end{array}\right.---\left(7\right)$

其中，x_min,x_max,y_min,y_max——相对位置范围在x,y方向的上下边界分量；

步骤四、寻找使步骤一所述目标函数取极小值时的x₀和β；

步骤五、用步骤四得到的x₀和β，实施追踪航天器周期性地沿水滴形轨迹飞行的控制；

具体过程为：悬停任务开始时，当追踪航天器到达相对坐标系中的A点，施加速度脉冲，使速度满足绕飞条件，追踪航天器开始沿着水滴形轨迹飞行，经过一个周期T_w后回到A点，此时再沿x方向施加维持水滴形绕飞轨迹所需要的速度脉冲增量△v，重新使速度满足绕飞条件，从而使追踪航天器周期性地沿水滴形轨迹飞行。

2.根据权利要求1所述的一种基于参数优化的航天器单脉冲水滴形绕飞轨迹悬停控制方法，其特征在于步骤四的具体过程为:

不断地选取不同的x₀和β，采用非线性规划的方法，使用Matlab中的函数fmincon对A点位置的x方向分量和追踪航天器水滴形绕飞轨迹的周期T_w与目标航天器的轨道周期T的比值β进行优化，得到追踪航天器悬停所需燃料消耗最少的水滴形绕飞轨迹。