CN105930305A - 一种三脉冲交会接近制导方法 - Google Patents

Abstract

本发明公开了一种三脉冲交会接近制导方法，包括以下步骤：利用C‑W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程；利用相对轨道动力学方程建立带修正脉冲的三脉冲C‑W制导状态转移方程；根据状态转移方程确定燃料最省和制导精度的优化目标函数，确定脉冲大小、方向和时间为优化参数；根据优化目标函数得到最优结果。本发明提供的三脉冲交会接近制导方法，在交会接近的初始脉冲与最终脉冲之间施加一个修正脉冲以达到提高制导精度的目的。

Description

一种三脉冲交会接近制导方法

技术领域

本发明涉及航空器导航、制导与控制领域，尤其涉及一种三脉冲交会接近制导方法。

背景技术

近年来，利用空间机器人对地球静止轨道(Geostationary Orbit，GEO)卫星进行在轨服务的研究成为热点，空间机器人进行在轨服务的前提是实现对GEO卫星的自主交会接近，因而快速、准确、最优的交会接近制导方法至关重要。

交会接近制导方法的主要任务是利用航天器之间的相对轨道动力学，以交会接近时间、位置精度和燃料消耗等为约束条件，设计合理的、安全的交会接近运动轨迹。航天工程上，在近圆轨道上交会接近相对轨道动力学的描述方法主要是Clohessy-Wiltshire(C-W)方程。以C-W方程为基础，一些学者采用遗传优化算法，求得了以燃料和时间最优为指标的交会接近双脉冲速度解；一些学者研究了双脉冲最优交会时，考虑第一次脉冲位置变化的优化问题；一些学者推导出了多脉冲交会的一般算法，并以此为基础着重研究了双脉冲燃料消耗问题。上述研究成果多数为双脉冲C-W制导方法，然而，由于GEO轨道周期长、轨道引力弱，在外部轨道环境干扰下、有限精度相对导航和控制条件下，采用双脉冲C-W制导方法，空间机器人交会接近GEO卫星的实际运动轨迹与理想运动轨迹偏差较大。随着交会接近距离越远、时间越长，最终位置误差就越大。

随着GEO轨道在轨服务技术的不断发展，对交会接近技术的要求越来越高，需要进一步提高空间机器人对GEO卫星进行远距离交会接近的制导精度，同时提高优化算法的性能。

发明内容

为解决上述技术问题，本发明提供了一种三脉冲交会接近制导方法，在交会接近的初始脉冲与最终脉冲之间施加一个修正脉冲以达到提高制导精度的目的。

为达到上述目的，本发明采用以下技术方案：

本发明公开了一种三脉冲交会接近制导方法，包括以下步骤：

S1：利用C-W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程；

S2：利用步骤S1中建立的相对轨道动力学方程建立带修正脉冲的三脉冲C-W制导状态转移方程；

S3：根据步骤S2中的状态转移方程确定燃料最省和制导精度的优化目标函数，确定脉冲大小、方向和时间为优化参数；

S4：根据步骤S3中的优化目标函数得到最优结果。

优选地，步骤S4进一步包括：采用改进动态权重-因子的粒子群优化算法对步骤S3中的优化目标函数进行优化得到最优结果。

优选地，步骤S1中利用C-W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程为：

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

其中，r＝[x,y,z]^T，a_c＝[a_x,a_y,a_z]^T，为GEO卫星的轨道角速率，μ是地心引力常数，a_c为外界加速度矢量，r_t为GEO卫星在地心惯性坐标系下的位置矢量，r为空间机器人与GEO卫星之间的相对位置矢量；

将作为状态变量，以相对轨道动力学方程为系统状态方程的状态转移方程线性解为：

$X\left(t\right)=\mathrm{\Φ}(t,t_0)X\left(t_0\right)+{\∫}_{t_0}^t\mathrm{\Φ}(t,u)BU\left(u\right)du$

其中，

$\mathrm{\Φ}(t,t_0)=\left[\begin{array}{cccccc}1& 0& 6(\mathrm{\Δ}\mathrm{\τ}-S_0)& (4S_0-3\mathrm{\Δ}\mathrm{\τ})/n& 0& 2(1-C_0)/n\\ 0& C_0& 0& 0& S_0/n& 0\\ 0& 0& 4-3C_0& -2(1-C_0)/n& 0& S_0/n\\ 0& 0& 6n(1-C_0)& 4C_0-3& 0& 2S_0\\ 0& -{\mathrm{nS}}_0& 0& 0& C_0& 0\\ 0& 0& 3{\mathrm{nS}}_0& -2S_0& 0& C_0\end{array}\right]$

${\∫}_{t_0}^t\mathrm{\Φ}(t,u)BU\left(u\right)du=\left[\begin{array}{ccc}4\left(1-C_0\right)/n^2-1.5{\mathrm{\Δt}}^2& 0& 2\left(\mathrm{\Δ}t/n-S_0/n^2\right)\\ 0& \left(1-C_0\right)/n^2& 0\\ 2\left(S_0/n^2-\mathrm{\Δ}t/n\right)& 0& \left(1-C_0\right)/n^2\\ 4S_0/n-3\mathrm{\Δ}t& 0& 2\left(1-C_0\right)/n\\ 0& S_0/n& 0\\ -2\left(1-C_0\right)/n& 0& S_0/n\end{array}\right]U\left(t_0\right)$

其中，△t＝t-t₀，△τ＝△tn，S₀＝sin(△τ)，C₀＝cos(△τ)。

优选地，步骤S2中建立的带修正脉冲的三脉冲C-W制导状态转移方位为：

$X\left(t_2\right)=\mathrm{\Φ}(t_2-t_0)X\left(t_0\right)+\underset{j=0}{\overset{2}{\Σ}}\mathrm{\Φ}(t_2-t_j)\left[\begin{array}{c}0\\ {\mathrm{\Δv}}_j\end{array}\right]$

其中，t₀为第一个脉冲的施加时间，t₁为第二个脉冲的施加时间，t₂为第三个脉冲的施加时间，△v_j为t_j时刻的速度脉冲。

优选地，步骤S3中的优化目标函数为：

T(y,M)＝f(y)+MP(y)

其中，y为待优化参数，M为罚因子，MP(y)为罚函数，P(y)＝||r-r_f||，r_f为空间机器人对GEO卫星进行三脉冲制导的实际最终相对位置，在进行迭代计算时，若P(y)≤R时M＝0，反之M取极大值，R为空间机器人与GEO卫星的最终相对位置精度要求；

其中，待优化参数为：

y＝(△v₀,t₁)＝(△V₀cosα,△V₀sinα,t₁)

其中，△v₀为t₀时刻的速度脉冲，△V₀为施加第一个脉冲的大小，α为施加第一个脉冲的方向。

优选地，步骤S4具体包括：

S41：设定粒子群优化算法的粒子个数、搜索维数和最大迭代次数；

S42：初始化步骤S3中的优化参数的初始粒子的位置和速率；

S43：以步骤S3中的优化目标函数为粒子群适应度函数，计算每个粒子的适应度值，从中确定第i个粒子局部最优适应度值P_{i_best}和所有粒子的全局最优适应度值P_{g_best}；

S44：采用改进动态权重-因子的位置、速率更新方程对步骤S42中的初始粒子的位置、速率进行更新，以实现参数的优化；

△y_i(k+1)＝ξ(k){ω(k)△y_i(k)+c₁(k)Rand₁(k)(P_{i_best}-y_i(k))

+c₂(k)Rand₂(k)(P_{g_best}-y_i(k))}

y_i(k+1)＝y_i(k)+△y_i(k+1)

其中，k为迭代次数，y_i(k)为第i个粒子位置，△y_i(k)为第i个粒子速率，Rand₁(k)和Rand₂(k)为[0，1]之间的随机数，ω(k)为动态惯性权重因子，c₁(k)为动态局部学习因子，c₂(k)为动态全局学习因子，ξ(k)为动态收缩因子；

S45：重复步骤S43和步骤S44，直到达到最大迭代次数或者全局最优适应度值P_{g_best}满足精度要求，最终得到最优结果。

优选地，步骤S43中的第i个粒子局部最优适应度值P_{i_best}和所有粒子的全局最优适应度值P_{g_best}分别为：

$P_{i\_best}=\left\{\begin{array}{cc}{y}_i\left(k\right),& ifk=1\\ y_i\left(k\right),& elseif{F}_i\left(k\right)<F_{i\_best}\\ P_{i\_best},& else\end{array}\right.$

$P_{g\_best}=\left\{\begin{array}{cc}min\left(y_i\right(k\left)\right),& ifk=1\\ min\left(y_i\right(k\left)\right),& elseifmin\left(F_i\right(k\left)\right)<F_{g\_best}\\ P_{g\_best},& else\end{array}\right.,$

其中，F_i(k)为粒子群适应度函数，F_{i_best}为第i个粒子的粒子群适应度函数最优值，F_{g_best}为所有粒子的粒子群适应度函数最优值。

优选地，动态惯性权重因子ω(k)根据下式进行更新：

$\mathrm{\&omega;}\left(k\right)=\frac{N_{max}-k}{N_{max}}({\mathrm{\&omega;}}_{max}-{\mathrm{\&omega;}}_{min})+{\mathrm{\&omega;}}_{min}$

其中，ω_max为动态权重最大值，ω_min为动态权重最小值，步骤S41还包括设定ω_max和ω_min。

优选地，动态局部学习因子c₁(k)和动态全局学习因子c₂(k)根据下式进行更新：

$\left\{\begin{array}{c}{c}_1\left(k\right)=\frac{c_{1-max}+c_{1-min}}{2}-\frac{c_{1-max}-c_{1-min}}{2}cos\left(\mathrm{\&phi;}\right(k\left)\right)\\ c_2\left(k\right)=\frac{c_{2-max}+c_{2-min}}{2}+\frac{c_{2-max}-c_{2-min}}{2}cos\left(\mathrm{\&phi;}\right(k\left)\right)\end{array}\right.$

其中，c_1-max为动态局部学习因子最大值，c_1-min为动态局部学习因子最小值，c_2-max为动态全局学习因子最大值，c_2-min为动态全局学习因子最小值，步骤S41还包括设定c_1-max、c_1-min、c_2-max和c_2-min。

优选地，动态收缩因子ξ(k)根据下式进行更新：

$\mathrm{\&xi;}\left(k\right)=\frac{2}{|2-c\left(k\right)-\sqrt{c{\left(k\right)}^2-4c\left(k\right)}|}$

其中，

与现有技术相比，本发明的有益效果在于：本发明在交会接近的初始脉冲与最终脉冲之间施加一个修正脉冲，与现有的双脉冲C-W制导仿真结果进行对比，本发明提高了空间机器人对GEO卫星进行远距离交会接近制导的最终位置精度。

在进一步的方案中，本发明提出一种改进的动态权重-因子的粒子群优化算法，与时间线性递减的惯性权重粒子群优化算法、带收缩因子粒子群优化算法仿真对比，本发明的粒子群优化算法不仅提高了前期全局搜索能力以发现高质量的局部解空间，而且也提高了后期局部解空间搜索能力以发现最优值。

附图说明

图1是空间机器人与GEO卫星的相对运动关系图及坐标系；

图2是本发明优选实施例的基于改进粒子群算法的三脉冲交会接近制导优化方法的步骤图；

图3是本发明优选实施例的改进动态权重-因子的粒子群优化算法流程图；

图4是空间机器人的双脉冲C-W制导最终误差统计图；

图5是改进动态权重-因子的粒子群优化算法性能对比图；

图6是基于改进粒子群优化算法的速度脉冲适应度函数最优值变化图；

图7是空间机器人的三脉冲C-W制导运动轨迹图；

图8是空间机器人的三脉冲C-W制导最终误差统计图。

具体实施方式

下面对照附图并结合优选的实施方式对本发明作进一步说明。

为描述空间机器人与GEO卫星之间的相对运动关系，建立如图1所示的坐标系，包括：地心惯性坐标系∑I、空间机器人轨道坐标系∑c和GEO卫星轨道坐标系∑t。空间机器人和GEO卫星在远距离交会接近时可以视作质点。

地心惯性坐标系∑I原点位于地心，基准面XY是历元为J2000.0时的地球平赤道平面，+X轴指向此历元时的平春分点方向，+Z轴指向北极，Y轴与XZ轴构成右手坐标系。空间机器人轨道坐标系∑c原点为空间机器人的质心；X_cZ_c为轨道平面，+Z_c指向地心的方向；+X_c垂直于+Z_c指向飞行方向；+Y_c垂直于轨道面的方向，形成右手直角坐标系。GEO卫星轨道坐标系∑t-X_tY_tZ_t定义与空间机器人相同。r_c为空间机器人在地心惯性坐标系下的位置矢量；r_t为GEO卫星在地心惯性坐标系下的位置矢量；r为两者之间的相对位置矢量。

如图2所示，本发明优选实施例的基于改进粒子群算法的三脉冲交会接近制导优化方法包括以下步骤：

S1：利用C-W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程。

空间机器人与GEO卫星之间的相对轨道动力学表达式为：

$\stackrel{\&CenterDot;\&CenterDot;}{r}=-\frac{{\mathrm{\&mu;r}}_c}{|r_c{|}^3}+\frac{{\mathrm{\&mu;r}}_t}{|r_t{|}^3}+a_c---\left(1\right)$

上式(1)中，μ是地心引力常数，约为3.98603×10¹⁴m³/s²；a_c为外界加速度矢量。

对上式(1)在矢量r_t处进行一阶泰勒级数展开，得到空间机器人与GEO卫星之间的相对轨道动力学C-W方程为：

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

上式(2)中，r＝[x,y,z]^T，a_c＝[a_x,a_y,a_z]^T，为GEO卫星的轨道角速率。

将作为状态变量，以C-W方程为系统状态方程的状态转移方程线性解为：

$X\left(t\right)=\mathrm{\&Phi;}(t,t_0)X\left(t_0\right)+{\&Integral;}_{t_0}^t\mathrm{\&Phi;}(t,u)BU\left(u\right)du---\left(3\right)$

上式(3)中，

$\mathrm{\&Phi;}(t,t_0)=\left[\begin{array}{cccccc}1& 0& 6(\mathrm{\&Delta;}\mathrm{\&tau;}-S_0)& (4S_0-3\mathrm{\&Delta;}\mathrm{\&tau;})/n& 0& 2(1-C_0)/n\\ 0& C_0& 0& 0& S_0/n& 0\\ 0& 0& 4-3C_0& -2(1-C_0)/n& 0& S_0/n\\ 0& 0& 6n(1-C_0)& 4C_0-3& 0& 2S_0\\ 0& -{\mathrm{nS}}_0& 0& 0& C_0& 0\\ 0& 0& 3{\mathrm{nS}}_0& -2S_0& 0& C_0\end{array}\right]---\left(4\right)$

${\&Integral;}_{t_0}^t\mathrm{\&Phi;}(t,u)BU\left(u\right)du=\left[\begin{array}{ccc}4\left(1-C_0\right)/n^2-1.5{\mathrm{\&Delta;t}}^2& 0& 2\left(\mathrm{\&Delta;}t/n-S_0/n^2\right)\\ 0& \left(1-C_0\right)/n^2& 0\\ 2\left(S_0/n^2-\mathrm{\&Delta;}t/n\right)& 0& \left(1-C_0\right)/n^2\\ 4S_0/n-3\mathrm{\&Delta;}t& 0& 2\left(1-C_0\right)/n\\ 0& S_0/n& 0\\ -2\left(1-C_0\right)/n& 0& S_0/n\end{array}\right]U\left(t_0\right)---\left(5\right)$

其中，△t＝t-t₀，△τ＝△tn，S₀＝sin(△τ)，C₀＝cos(△τ)。

S2：利用步骤S1中建立的相对轨道动力学状态转移方程建立带修正脉冲的三脉冲C-W制导状态转移方程。

其中：双脉冲C-W制导只在转移初始时刻t₀和结束时刻t₁施加速度脉冲，转移过程并不进行控制。由于施加速度脉冲的时间很短，可以将上式(3)相对运动状态转移方程简化为：

$X\left(t_1\right)=\mathrm{\&Phi;}(t_1-t_0)\left(X\right(t_0)+\left[\begin{array}{c}0\\ {\mathrm{\&Delta;v}}_0\end{array}\right])+\left[\begin{array}{c}0\\ {\mathrm{\&Delta;v}}_1\end{array}\right]---\left(6\right)$

上式(6)中，△v₀为初始时刻t₀的速度脉冲，△v₁为结束时刻t₁的速度脉冲。状态转移矩阵拆分成4个3×3的子矩阵：

$\mathrm{\&Phi;}=\left[\begin{array}{cc}{\mathrm{\&Phi;}}_{rr}& {\mathrm{\&Phi;}}_{rv}\\ {\mathrm{\&Phi;}}_{vr}& {\mathrm{\&Phi;}}_{vv}\end{array}\right]---\left(7\right)$

利用状态转移矩阵的性质：

Φ(t_m+t_n)＝Φ(t_m)Φ(t_n) (8)

得到三脉冲C-W制导的状态转移方程为：

$X\left(t_2\right)=\mathrm{\&Phi;}(t_2-t_0)X\left(t_0\right)+\underset{j=0}{\overset{2}{\&Sigma;}}\mathrm{\&Phi;}(t_2-t_j)\left[\begin{array}{c}0\\ {\mathrm{\&Delta;v}}_j\end{array}\right]---\left(9\right)$

上式(9)中，△v_j为t_j时刻的速度脉冲。

S3：利用步骤S2中的状态转移方程确定燃料最省和制导精度的优化目标函数，确定脉冲大小、方向和时间为优化参数。

以固定时间燃料最省为优化问题，则三脉冲优化目标函数表示为：

$J=\min \left(\underset{j=0}{\overset{2}{\&Sigma;}}\right|{\mathrm{\&Delta;v}}_j\left|\right)=\min f\left(y\right)---\left(10\right)$

上式(10)中，y为待优化参数。

增加最终制导精度要求，作为固定时间燃料最省优化问题的约束函数：

P(y)＝||r-r_f||≤R (11)

上式(11)中，r为空间机器人对GEO卫星进行三脉冲制导的期望最终相对位置，r_f为实际最终相对位置，R为最终相对位置精度要求。

空间机器人的三脉冲制导优化问题转化为带约束的非线性最优问题，约束条件的处理是最优问题求解的重要环节，引入罚函数可以有效地将带约束的非线性优化问题转化为不带约束的非线性优化问题。因此，得到空间机器人三脉冲C-W制导新的优化目标函数为：

T(y,M)＝f(y)+MP(y) (12)

上式(12)中，MP(y)为罚函数，M为罚因子。当进行优化迭代计算时，若P(y)≤R满足相对位置精度要求时，M＝0；反之，则M取极大值。

三脉冲C-W制导待优化参数包括：第一个脉冲△v₀和t₀、第二个脉冲△v₁和t₁、第三个脉冲△v₂和t₂。由于初始时刻t₀、结束时刻t₂、初始位置r₀、结束位置r_f为已知；同时，速度脉冲△v₁和△v₂可以通过状态转移方程计算得到：

$\left\{\begin{array}{c}{r}_1={\mathrm{\&Phi;}}_{rr}\left(t_1-t_0\right)r_0+{\mathrm{\&Phi;}}_{rv}\left(t_1-t_0\right)\&lsqb;v_0+{\mathrm{\&Delta;v}}_0\&rsqb;\\ v_1={\mathrm{\&Phi;}}_{vr}\left(t_1-t_0\right)r_0+{\mathrm{\&Phi;}}_{vv}\left(t_1-t_0\right)\&lsqb;v_0+{\mathrm{\&Delta;v}}_0\&rsqb;\\ {\mathrm{\&Delta;v}}_1={\mathrm{\&Phi;}}_{rv}^{-1}\left(t_2-t_1\right)\&lsqb;r_f-{\mathrm{\&Phi;}}_{rr}\left(t_2-t_1\right)r_1\&rsqb;-v_1\\ {\mathrm{\&Delta;v}}_2=v_f-{\mathrm{\&Phi;}}_{vr}\left(t_2-t_1\right)r_1-{\mathrm{\&Phi;}}_{vv}\left(t_2-t_1\right)\&lsqb;v_1+{\mathrm{\&Delta;v}}_1\&rsqb;\end{array}\right.---\left(13\right)$

上式(13)中，v₀为第一个脉冲前的相对速度，v₁为第二个脉冲前的相对速度，v_f为第三个脉冲后的相对速度，r₁为第二个脉冲时刻的相对位置。

另外，空间机器人将进行轨道面内的远距离交会接近，Y轴方向的速度脉冲为零，那么空间机器人的三脉冲C-W制导待优化参数简化为：

y＝(△v₀,t₁)＝(△V₀cosα,△V₀sinα,t₁) (14)

上式(14)中，△V₀为施加第一个脉冲的大小，α为施加第一个脉冲的方向，t₁为第二个脉冲的施加时间。

S4：采用改进动态权重-因子的粒子群优化算法对步骤S3中的优化目标函数进行优化得到最优结果，具体流程如图3所示，包括以下步骤：

S41：设定粒子群优化算法的粒子个数n_p、搜索维数D、最大迭代次数N_max、罚因子M、动态权重最大值ω_max、动态权重最小值ω_min、动态局部学习因子最大值c_1-max、动态局部学习因子最小值c_1-min、动态全局学习因子最大值c_2-max、动态全局学习因子最小值c_2-min；

S42：初始化步骤S3中的优化参数的n_p个粒子的位置y₁(k)和速率△y₁(k)；

S43：以步骤S3中的式(12)的优化目标函数为粒子群适应度函数：

$F_i\left(k\right)=T\left(y_i\right(k),M(k\left)\right),\left\{\begin{array}{cc}M\left(k\right)=0,& ifP\left(y_i\right(k\left)\right)\&le;R\\ M\left(k\right)=M,& else\end{array}\right.---\left(15\right)$

计算每个粒子的适应度值，从中确定第i个粒子局部最优适应度值P_{i_best}：

$F_{i\_best}=\left\{\begin{array}{cc}{F}_i\left(k\right),& ifk=1\\ F_i\left(k\right),& elseif{F}_i\left(k\right)<F_{i\_best}\\ F_{i\_best},& else\end{array}\right.---\left(16\right)$

$P_{i\_best}=\left\{\begin{array}{cc}{y}_i\left(k\right),& ifk=1\\ y_i\left(k\right),& elseif{F}_i\left(k\right)<F_{i\_best}\\ P_{i\_best},& else\end{array}\right.---\left(17\right)$

从局部最优适应度中，确定所有粒子的全局最优适应度值P_{g_best}：

$F_{g\_best}=\left\{\begin{array}{cc}min\left(F_i\right(k\left)\right),& ifk=1\\ min\left(F_i\right(k\left)\right),& elseifmin\left(F_i\right(k\left)\right)<F_{g\_best}\\ F_{g\_best},& else\end{array}\right.---\left(18\right)$

$P_{g\_best}=\left\{\begin{array}{cc}min\left(y_i\right(k\left)\right),& ifk=1\\ min\left(y_i\right(k\left)\right),& elseifmin\left(F_i\right(k\left)\right)<F_{g\_best}\\ P_{g\_best},& else\end{array}\right.---\left(19\right)$

S44：利用改进动态权重-学习因子的位置、速率更新方程对步骤S42中粒子的初始位置、速率进行更新，以实现参数的优化；

$\begin{array}{c}{\mathrm{\&Delta;y}}_i\left(k+1\right)=\mathrm{\&xi;}\left(k\right)\{\mathrm{\&omega;}\left(k\right){\mathrm{\&Delta;y}}_i\left(k\right)+c_1\left(k\right){\mathrm{Rand}}_1\left(k\right)\left(P_{i\_best}-y_i\left(k\right)\right)\\ +c_2\left(k\right){\mathrm{Rand}}_2\left(k\right)\left(P_{g\_best}-y_i\left(k\right)\right)\}\end{array}---\left(20\right)$

y_i(k+1)＝y_i(k)+△y_i(k+1) (21)

上述式(20)和(21)中，y_i(k)为第i个粒子位置，△y_i(k)为第i个粒子速率，Rand₁(k)和Rand₂(k)为[0，1]之间的随机数，ω(k)为动态惯性权重因子，c₁(k)为动态局部学习因子，c₂(k)为动态全局学习因子，ξ(k)为动态收缩因子。

在进一步的实施例中，根据下式采用动态惯性权重：

$\mathrm{\&omega;}\left(k\right)=\frac{N_{max}-k}{N_{max}}({\mathrm{\&omega;}}_{max}-{\mathrm{\&omega;}}_{min})+{\mathrm{\&omega;}}_{min}---\left(22\right)$

上式(22)中，ω_max为动态权重最大值，ω_min动态权重最小值。通过采用动态惯性权重，使得在粒子群优化算法中，在运行初期具有更强的全局搜索能力以发现高质量的局部解空间，在运行后期则需要细致的局部解空间搜索能力。

在更进一步的实施例中提出一种随惯性权重动态变化的学习因子改进方法：

$\left\{\begin{array}{c}{c}_1\left(k\right)=\frac{c_{1-max}+c_{1-min}}{2}-\frac{c_{1-max}-c_{1-min}}{2}cos\left(\mathrm{\&phi;}\right(k\left)\right)\\ c_2\left(k\right)=\frac{c_{2-max}+c_{2-min}}{2}+\frac{c_{2-max}-c_{2-min}}{2}cos\left(\mathrm{\&phi;}\right(k\left)\right)\end{array}\right.---\left(23\right)$

上式(23)中，c_1-max为动态局部学习因子最大值，c_1-min为动态局部学习因子最小值，c_2-max为动态全局学习因子最大值，c_2-min为动态全局学习因子最小值。通过采用随惯性权重动态变化的学习因子，使得在粒子群优化算法中，在寻优的运行初期全局学习因子c₂较大，局部学习因子c₁较小，从而使得使粒子群的全局搜索能力增强；在寻优的运行后期局部学习因子c₁较大，全局学习因子c₂较小，粒子多向自身最优值学习，从而使得粒子群局部细致解的搜索能力增强。

进一步地，为了确保粒子群算法的收敛，在速度更新方程中增加动态收缩因子：

$\mathrm{\&xi;}\left(k\right)=\frac{2}{|2-c\left(k\right)-\sqrt{c{\left(k\right)}^2-4c\left(k\right)}|}---\left(24\right)$

上式(24)中，

S45：重复步骤S43和步骤S44直到达到最大迭代次数或者全局最优适应度值P_{g_best}满足精度要求，最终得到三脉冲C-W制导的最优结果。

下述对现有技术的双脉冲C-W制导和本发明的三脉冲C-W制导方法进行对比仿真。假定空间机器人对GEO卫星进行交会接近制导的初始位置r₀＝[-200.0,0,10.0]^Tkm、最终位置r_f＝[-8.0,0,0]^Tkm，初始速度为v₀＝[-1.0,0,0.5]^Tm/s，最终速度为v_f＝[0,0,0]^Tm/s，制导时间为18000s，相对位置确定误差为500m，相对速度确定误差为0.3m/s，速度脉冲控制误差为5％。要求交会接近制导的最终位置误差R<2km。

(1)双脉冲C-W制导仿真

通过上述条件，利用式(6)的双脉冲相对运动状态转移方程计算得到的速度脉冲为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_0={\left[\begin{array}{ccc}8.2632& 0& 8.2194\end{array}\right]}^{T}m/s\\ {\mathrm{\&Delta;v}}_1={\left[\begin{array}{ccc}-5.8088& 0& 9.6671\end{array}\right]}^{T}m/s\end{array}\right.---\left(25\right)$

基于以上两个速度脉冲，利用蒙特卡洛打靶法进行最终制导位置误差仿真，进行1000次仿真计算，误差统计结果如图4所示。由统计结果，位置误差在要求范围内(方框所示)的比例为23.5％，绝大多数不能满足最终位置误差要求，其中Z轴方向的最大偏差超过了6km，X轴方向的最大偏差则超过了10km。

(2)三脉冲C-W制导仿真

首先，利用优化评价函数Rastrigin函数，将本发明的改进动态权重-因子粒子群算法(Linear Decreasing Inertia Weight-PSO，DIWF-PSO)与时间线性递减的惯性权重粒子群算法(LDIW-PSO)、带收缩因子粒子群算法(ConstraintFactor-PSO，CF-PSO)进行优化性能比较，结果如图5所示。由结果可以看出，本发明的优化方法搜索时间略优于另外两种方法，且全局最优适应度值明显优于另外两种方法。

然后，利用本发明改进动态权重-因子粒子群算法的三脉冲交会接近制导优化方法进行远距离交会接近制导优化。将本发明的离子群优化算法的参数设定为：

$\left\{\begin{array}{c}{n}_p=50\\ D=3\\ N_{\max }=500\\ M={10}^5\\ {\mathrm{\&omega;}}_{\max }=0.9,{\mathrm{\&omega;}}_{\min }=0.5\\ c_{1-\max }=c_{2-\max }=2.5,c_{1-\min }=c_{2-\min }=1.6\end{array}\right.---\left(26\right)$

结合上式(26)的参数设置，根据本发明的方法得到的优化结果为：

$\left\{\begin{array}{c}{t}_1=15820s\\ {\mathrm{\&Delta;V}}_0=11.6551m/s\\ \mathrm{\&alpha;}=0.8169rad\\ {\mathrm{\&Delta;v}}_0={\&lsqb;7.9777,0,8.4969\&rsqb;}^{T}m/s\end{array}\right.---\left(27\right)$

在迭代过程中，空间机器人速度脉冲全局最优适应度值变化如图6所示。利用三脉冲C-W制导状态转移方程，得到空间机器人t₁中间时刻速度脉冲、t₂结束时刻速度脉冲以及总的速度脉冲增量为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;v}}_1={\&lsqb;-0.0898,0,1.2298\&rsqb;}^{T}m/s\\ {\mathrm{\&Delta;v}}_2={\&lsqb;-5.2885,0,8.8874\&rsqb;}^{T}m/s\\ J=\underset{j=0}{\overset{2}{\&Sigma;}}\left|{\mathrm{\&Delta;v}}_j\right|=23.2301m/s\end{array}\right.---\left(28\right)$

最后，利用上述式(27)和(28)三脉冲C-W制导优化结果进行交会接近制导，得到空间机器人的运动轨迹如图7所示；同样，利用蒙特卡洛打靶法进行最终制导位置误差仿真，进行1000次仿真计算，空间机器人的最终制导误差统计结果如图8所示。与双脉冲制导误差统计结果图4进行对比，本发明的三脉冲制导误差统计结果明显偏小，且全部满足最终位置误差要求。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。

Claims (10)

1.一种三脉冲交会接近制导方法，其特征在于，包括以下步骤：

S1：利用C-W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程；

S2：利用步骤S1中建立的相对轨道动力学方程建立带修正脉冲的三脉冲C-W制导状态转移方程；

S3：根据步骤S2中的状态转移方程确定燃料最省和制导精度的优化目标函数，确定脉冲大小、方向和时间为优化参数；

S4：根据步骤S3中的优化目标函数得到最优结果。

2.根据权利要求1所述的三脉冲交会接近制导方法，其特征在于，步骤S4进一步包括：采用改进动态权重-因子的粒子群优化算法对步骤S3中的优化目标函数进行优化得到最优结果。

3.根据权利要求1所述的三脉冲交会接近制导方法，其特征在于，步骤S1中利用C-W方程建立空间机器人与GEO卫星之间的相对轨道动力学方程为：

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

其中，r＝[x,y,z]^T，a_c＝[a_x,a_y,a_z]^T，为GEO卫星的轨道角速率，μ是地心引力常数，a_c为外界加速度矢量，r_t为GEO卫星在地心惯性坐标系下的位置矢量，r为空间机器人与GEO卫星之间的相对位置矢量；

将作为状态变量，以相对轨道动力学方程为系统状态方程的状态转移方程线性解为：

$X\left(t\right)=\mathrm{\&Phi;}(t,t_0)X\left(t_0\right)+{\&Integral;}_{t_0}^t\mathrm{\&Phi;}(t,u)BU\left(u\right)du$

其中，

$\mathrm{\&Phi;}(t,t_0)=\left[\begin{array}{cccccc}1& 0& 6(\mathrm{\&Delta;}\mathrm{\&tau;}-S_0)& (4S_0-3\mathrm{\&Delta;}\mathrm{\&tau;})/n& 0& 2(1-C_0)/n\\ 0& C_0& 0& 0& S_0/n& 0\\ 0& 0& 4-3C_0& -2(1-C_0)/n& 0& S_0/n\\ 0& 0& 6n(1-C_0)& 4C_0-3& 0& 2S_0\\ 0& -{\mathrm{nS}}_0& 0& 0& C_0& 0\\ 0& 0& 3{\mathrm{nS}}_0& -2S_0& 0& C_0\end{array}\right]$

${\&Integral;}_{t_0}^t\mathrm{\&Phi;}(t,u)BU\left(u\right)du=\left[\begin{array}{ccc}4(1-C_0)/n^2-1.5{\mathrm{\&Delta;t}}^2& 0& 2(\mathrm{\&Delta;}t/n-S_0/n^2)\\ 0& (1-C_0)/n^2& 0\\ 2(S_0/n^2-\mathrm{\&Delta;}t/n)& 0& (1-C_0)/n^2\\ 4S_0/n-3\mathrm{\&Delta;}t& 0& 2(1-C_0)/n\\ 0& S_0/n& 0\\ -2(1-C_0)/n& 0& S_0/n\end{array}\right]U\left(t_0\right)$

其中，Δt＝t-t₀，Δτ＝Δtn，S₀＝sin(Δτ)，C₀＝cos(Δτ)。

4.根据权利要求3所述的三脉冲交会接近制导方法，其特征在于，步骤S2中建立的带修正脉冲的三脉冲C-W制导状态转移方位为：

$X\left(t_2\right)=\mathrm{\&Phi;}(t_2-t_0)X\left(t_0\right)+\underset{j=0}{\overset{2}{\mathrm{\&Sigma;}}}\mathrm{\&Phi;}(t_2-t_j)\left[\begin{array}{c}0\\ \mathrm{\&Delta;}{v}_j\end{array}\right]$

其中，t₀为第一个脉冲的施加时间，t₁为第二个脉冲的施加时间，t₂为第三个脉冲的施加时间，Δv_j为t_j时刻的速度脉冲。

5.根据权利要求4所述的三脉冲交会接近制导方法，其特征在于，步骤S3中的优化目标函数为：

T(y,M)＝f(y)+MP(y)

其中，y为待优化参数，M为罚因子，MP(y)为罚函数，P(y)＝||r-r_f||，r_f为空间机器人对GEO卫星进行三脉冲制导的实际最终相对位置，在进行迭代计算时，若P(y)≤R时M＝0，反之M取极大值，R为空间机器人与GEO卫星的最终相对位置精度要求；

其中，待优化参数为：

y＝(Δv₀,t₁)＝(ΔV₀cosα,ΔV₀sinα,t₁)

其中，Δv₀为t₀时刻的速度脉冲，ΔV₀为施加第一个脉冲的大小，α为施加第一个脉冲的方向。

6.根据权利要求1至5任一项所述的三脉冲交会接近制导方法，其特征在于，步骤S4具体包括：

S41：设定粒子群优化算法的粒子个数、搜索维数和最大迭代次数；

S42：初始化步骤S3中的优化参数的初始粒子的位置和速率；

S43：以步骤S3中的优化目标函数为粒子群适应度函数，计算每个粒子的适应度值，从中确定第i个粒子局部最优适应度值P_{i_best}和所有粒子的全局最优适应度值P_{g_best}；

S44：采用改进动态权重-因子的位置、速率更新方程对步骤S42中的初始粒子的位置、速率进行更新，以实现参数的优化；

Δy_i(k+1)＝ξ(k){ω(k)Δy_i(k)+c₁(k)Rand₁(k)(P_{i_best}-y_i(k))

+c₂(k)Rand₂(k)(P_{g_best}-y_i(k))}

y_i(k+1)＝y_i(k)+Δy_i(k+1)

其中，k为迭代次数，y_i(k)为第i个粒子位置，Δy_i(k)为第i个粒子速率，Rand₁(k)和Rand₂(k)为[0，1]之间的随机数，ω(k)为动态惯性权重因子，c₁(k)为动态局部学习因子，c₂(k)为动态全局学习因子，ξ(k)为动态收缩因子；

S45：重复步骤S43和步骤S44，直到达到最大迭代次数或者全局最优适应度值P_{g_best}满足精度要求，最终得到最优结果。

7.根据权利要求6所述的三脉冲交会接近制导方法，其特征在于，步骤S43中的第i个粒子局部最优适应度值P_{i_best}和所有粒子的全局最优适应度值P_{g_best}分别为：

$P_{i\_best}=\left\{\begin{array}{cc}{y}_i\left(k\right),& ifk=1\\ y_i\left(k\right),& elseif{F}_i\left(k\right)<F_{i\_best}\\ P_{i\_best},& else\end{array}\right.$

$P_{g\_best}=\left\{\begin{array}{cc}\min \left(y_i\right(k\left)\right),& ifk=1\\ \min \left(y_i\right(k\left)\right),& elseif\min \left(F_i\right(k\left)\right)<F_{g\_best}\\ P_{g\_best},& else\end{array}\right.,$

其中，F_i(k)为粒子群适应度函数，F_{i_best}为第i个粒子的粒子群适应度函数最优值，F_{g_best}为所有粒子的粒子群适应度函数最优值。

8.根据权利要求6所述的三脉冲交会接近制导方法，其特征在于，动态惯性权重因子ω(k)根据下式进行更新：

$\mathrm{\&omega;}\left(k\right)=\frac{N_{max}-k}{N_{max}}({\mathrm{\&omega;}}_{max}-{\mathrm{\&omega;}}_{\min })+{\mathrm{\&omega;}}_{min}$

其中，ω_max为动态权重最大值，ω_min为动态权重最小值，步骤S41还包括设定ω_max和ω_min。

9.根据权利要求8所述的三脉冲交会接近制导方法，其特征在于，动态局部学习因子c₁(k)和动态全局学习因子c₂(k)根据下式进行更新：

$\left\{\begin{array}{c}{c}_1\left(k\right)=\frac{c_{1-max}+c_{1-min}}{2}-\frac{c_{1-max}-c_{1-min}}{2}cos\left(\mathrm{\&phi;}\left(k\right)\right)\\ c_2\left(k\right)=\frac{c_{2-max}+c_{2-\min }}{2}+\frac{c_{2-max}-c_{2-min}}{2}cos\left(\mathrm{\&phi;}\left(k\right)\right)\end{array}\right.$

其中，c_1-max为动态局部学习因子最大值，c_1-min为动态局部学习因子最小值，c_2-max为动态全局学习因子最大值，c_2-min为动态全局学习因子最小值，步骤S41还包括设定c_1-max、c_1-min、c_2-max和c_2-min。

10.根据权利要求9所述的三脉冲交会接近制导方法，其特征在于，动态收缩因子ξ(k)根据下式进行更新：

$\mathrm{\&xi;}\left(k\right)=\frac{2}{|2-c\left(k\right)-\sqrt{c{\left(k\right)}^2-4c\left(k\right)}|}$

其中，