基于多项式逼近的行星际小推力转移轨道设计方法

技术领域

本发明涉及一种基于多项式的行星际小推力转移轨道快速设计方法，适用于行星际小推 力转移轨道的初始设计，属于深空探测转移轨道技术领域。

背景技术

在行星际探测任务中，探测器需要消耗多少燃料才能到达目标星体是任务设计首要关注 的问题。相比传统化学推进系统，小推力发动机具有比冲高、质量轻的特点，行星际探测器 利用其实现轨道转移时可以有效节省燃料消耗。由于小推力发动机需要长时间工作才能达到 改变轨道的目的，使得探测器轨道为典型的强非线性非开普勒轨道，传统脉冲轨道设计中的 许多理论与方法不再适用，寻求一种快速有效转移轨道设计方法成为目前研究的热点。基于 形状曲线逼近方法是目前实现转移轨道快速设计最为有效的途径，其是借鉴圆锥曲线描述脉 冲轨道的思想，利用合适的函数曲线描述小推力轨道形状，获得形状曲线与小推力转移轨道 间的参数关系，进而有效降低转移轨道设计难度。选取何种形状曲线对小推力轨道进行逼近， 在此基础上如何解算小推力轨道参数，决定了探测器转移轨道设计正确与否及设计效率。因 此基于形状曲线的小推力转移轨道设计方法是当前科技人员关注的重点问题之一。

在已发展的基于形状曲线的小推力转移轨道设计方法中，在先技术[1](Petropoulos A E， Longuski J M.Shape-based algorithm for automated design of low-thrust， gravity-assist trajectories[J].Journal of Spacecraft and Rockets，2004，41(5)： 87-796.)，采用正弦指数曲线对小推力轨道进行逼近，该方法首先通过假设推力方向始终沿 探测器速度方向或反方向，获得了转移过程中小推力发动机加速度等参数的解析表达式，然 后通过历遍搜索正弦指数函数参数，获得满足任务约束的最佳转移轨道。由于正弦指数函数 自由参数数量少，求得的转移轨道只能满足探测器始末端的位置约束，因此只能用于飞越型 轨道的设计；另外，该方法采用的是历遍搜索方式设计轨道，存在计算量大、设计效率低等 缺点。

在先技术[2](参见Wall B.J.and Conway B.A.Shape-Based Approach to Low-Thrust Rendezvous Trajectory Design[J].Journal of Guidance，Control，and Dynamics，2009， 32(1)：95-101.)，采用逆六次多项式对小推力轨道进行逼近，该方法同样也运用了推力方向 始终沿探测器速度方向或反方向的假设，求得了小推力轨道参数的解析解，然后采用遗传算 法获得最佳的多项式参数。由于逆六次多项式具有七个自由参数，可同时满足探测器始末端 位置速度以及转移时间约束，因此该方法可有效用于交会型转移轨道的设计。但是该方法仍 然无法突破推力方向假设，并且由于其是以转移相角为形状曲线自变量，在自由参数选择不 适的情况下，飞行时间约束很难满足，这直接影响了转移轨道设计的鲁棒性和效率。

发明内容

本发明针对目前基于形状曲线的行星际小推力转移轨道设计中存在推力方向假设、鲁棒 性差、效率低等问题，给出了一种基于多项式逼近的行星际小推力转移轨道设计方法，直接 避免了飞行时间约束，提高了鲁棒性。

该基于多项式逼近的行星际小推力转移轨道设计方法，包括以下步骤：

第一步：设计变量初值猜测，给定转移轨道设计变量的初值猜测；

第二步：计算探测器的始末端边界条件：通过读取行星星历文件，根据预设t₀时刻得到 出发星体的日心笛卡尔坐标系的位置矢量r_L和速度矢量v_L，根据预设t_f时刻得到目标星体的 位置矢量r_A和速度矢量v_A，得到探测器的始末端边界条件为$\left(\begin{array}{c}{r}_0=r_L\\ v_0=v_L+V_{L\infty }\end{array}\right),\left(\begin{array}{c}{r}_f=r_A\\ v_f=v_A+V_{A\infty }\end{array}\right);$

其中r₀为探测器发射时的日心位置矢量，v₀为探测器发射时的日心速度矢量，V_L∞为探测 器发射时相对地球的日心速度矢量；r_f为探测器到达目标星体时的日心位置矢量，v_f为到达 目标星体时的日心速度矢量，V_A∞为探测器到达目标星体时相对目标星体的日心速度矢量； 然后将笛卡尔坐标系中的边界条件转换到球坐标系中，并对相位角进行修正：

$Y_0^s=C_{\mathrm{sc}}^{-1}\left(t_0\right){[r_0^T,v_0^T]}^T$

$Y_f^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_f^T,v_f^T]}^T$

其中为探测器发射时在日心球坐标系中的轨道状态，为探测器到达目标星体时在日 心球坐标系中的轨道状态；

第三步：采用第二类切比雪夫多项式拟合探测器的转移轨道，多项式矩阵可以表示为

其中τ为归一化的时间变量，t₀时刻对应τ＝-1，t₀+t_f时刻对应τ＝1；然后利用球坐 标系下的探测器始末端边界条件计算切比雪夫多项式系数；

第四步：计算性能指标和约束条件：根据获得的切比雪夫多项式系数，计算转移轨道的 性能指标和小推力发动机推力约束；

第五步：根据计算的性能指标判断是否满足最优性条件，根据计算的推力约束判断是否 满足可行性条件，如果都满足，则优化成功，获得最佳转移轨道，如果有一项不满足，则调 整第一步中设计变量的初值猜测直至优化成功。

本发明方法利用切比雪夫多项式逼近小推力转移轨道形状，以时间为自变量避免了飞行 时间约束，增强了算法鲁棒性；然后通过探测器始末端轨道状态约束确定多项式系数，并直 接获得了转移过程中的加速度解析表达式，避免了推力方向假设，提高了计算效率；最后将 转移轨道设计问题转化为仅具有简单参数的混合整数非线性规划问题，有效降低了轨道设计 难度。该方法能够根据给定的始末端边界条件对不同任务类型的小推力转移轨道进行快速设 计。

附图说明

图1为行星际探测器转移轨道示意图。

具体实施方式

该行星际小推力转移轨道快速设计方法分采用迭代寻优计算方式求解，求解过程分为设 计变量初值猜测、始末端边界条件计算、多项式系数计算、性能指标和约束条件计算、设计 变量调整五个部分。

1)设计变量初值猜测

给定转移轨道设计变量的初值猜测Z₀。

2)始末端边界条件计算

通过读取行星星历文件，根据t₀得到出发星体的日心笛卡尔坐标系的位置速度矢量r_L和 v_L，根据t_f得到目标星体的位置速度矢量r_A和v_A，则可以得到探测器的始末端边界条件为

r₀＝r_L

v₀＝v_L+V_L∞

r_f＝r_A

v_f＝v_A+V_A∞

将笛卡尔坐标系中的边界条件转换到球坐标系中为

$Y_0^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_0^T,v_0^T]}^T$

$Y_f^s=C_{\mathrm{sc}}^{-1}\left(t_f\right){[r_f^T,v_f^T]}^T$

然后重新计算探测器在末端的方位角进行相位修正

θ_f＝θ_f+2πP

3)多项式系数计算

分别令τ＝-1和τ＝1，利用下式计算第二类切比雪夫多项式矩阵M(-1)和M(1)

然后利用球坐标系下的探测器始末端边界条件计算切比雪夫多项式系数

其中s为

4)性能指标和约束条件计算

根据获得的切比雪夫多项式系数，计算转移轨道的性能指标

$J={\int }_{t_0}^{t_f}\left|\right|A\left(t\right)\left|\right|\mathrm{dt}$

计算小推力发动机推力约束

Φ＝max(||A(t)||)-A_max

5)设计变量调整

寻优算法根据计算的性能指标判断是否满足最优性条件，根据计算的推力约束判断是否 满足可行性条件。如果都满足，则优化成功，获得最佳转移轨道，如果有一项不满足，则调 整设计变量，并重复步骤2)～5)，直至优化成功。

行星际探测器的转移轨道位置状态X可以在以黄道面为参考的日心球坐标系中描述， 其为飞行时间t的函数，可表示为

其中r，θ，分别为探测器在球坐标系中的日心距、方位角和仰角。

小推力发动机推力很小，需要长时间加速才能达到改变轨道的目的，这使得行星际转移 轨道一般为多圈螺旋线形状。切比雪夫多项式在曲线拟合方面具有很好的特性，选取飞行时 间为自变量的切比雪夫多项式分别对三个位置状态进行逼近得到

$r\left(t\right)=\underset{n=0}{\overset{N}{\Sigma }}{a}_n^r{U}_n\left(\tau \right)$

$\theta \left(t\right)=\underset{n=0}{\overset{N}{\Sigma }}{a}_n^{\theta }{U}_n\left(\tau \right)$

其中N为切比雪夫多项式次数，U_n(τ)为n阶第二类切比雪夫多项式，分别为 日心距、方位角和仰角对应的切比雪夫多项式系数，τ∈[-1，1]为归一化的转移时间，其与实 际飞行时间的关系为

$\tau =\frac{2(t-t_0)}{t_f-t_0}-1$

其中t₀为初始时刻，t_f为末端时刻。对上述三个切比雪夫逼近多项式求一阶导数可以得到探 测器速度各方向分量的表示为

其中n阶第二类切比雪夫多项式一阶导数的递归公式为

并且有

由探测器的位置和速度切比雪夫逼近式可以看出，如果多项式系数确定，则探 测器的轨道状态确定。探测器在行星间转移时，轨道必须满足行星的星历约束、飞行时间等 约束，既然切比雪夫多项式是以时间为自变量的，飞行时间约束已经满足，因此多项式系数 可仅通过星历约束就可获得。假定探测器在初始时刻的轨道状态约束为 在终端时刻的轨道状态约束为令

则M(τ)为2×(N+1)的矩阵，将始末端轨道状态约束分别带入切比雪夫逼近式中可以得到

由上述三式可以看出，当N＝3时方程组存在唯一解，这表明在满足探测器始末端轨道 状态约束的条件下，三次切比雪夫多项式就可实现小推力转移轨道的逼近。对方程进行求解 可以得到切比雪夫系数为

其中s为r，θ或

为了计算探测器在转移过程中的小推力加速度，将探测器球坐标系下轨道状态转换到笛 卡尔坐标系

Y^c(t)＝C_sc(t)Y^s(t)

其中Y^c＝[r_x，r_y，r_z，v_x，v_y，v_z]^T，C_sc为球坐标系到笛卡尔坐标系的转换矩阵，则小推力发动 机加速度可以表示为

其中μ为太阳的引力加速度常数。

由上式可以看出，小推力发动机在转移过程中的加速度大小和方向是由切比雪夫多项式 系数决定的，加速度方向无需提前进行假设限定，而加速度大小通常要受到最大推力的限制， 因此对转移轨道进行设计时必须考虑推力大小约束

max(||A(t)||)≤A_max

其中A_max为任务采用的小推力发动机允许的最大推力。

综上所述，影响切比雪夫多项式系数的变量为探测器的发射时间t₀和双曲线超速V_L∞， 到达目标行星的时间t_f和双曲线超速V_A∞，探测器小推力转移轨道的整圈数P，则轨道设计 问题可以归结为

设计参数

Z＝[t₀，t_f，V_L∞，V_A∞，P]^T

性能指标

$J\left(Z\right)={\int }_{t_0}^{t_f}\left|\right|A\left(t\right)\left|\right|\mathrm{dt}\to \min$

约束条件

Φ(Z)＝max(||A(t)||)-A_max≤0

以上归结的问题为简单的混合整数非线性规划问题，可以采用微分进化等算法直接进行 求解，最优解即对应满足任务约束的最佳小推力转移轨道。