一种小天体附着探测下降轨迹优化方法与流程

本发明涉及一种小天体附着探测下降轨迹优化方法，属于航天航空领域。

背景技术：

小天体探测作为人类了解太阳系形成与演化、生命起源与进化以及防御外来天体撞击的重要途径，将是未来深空探测活动的主要内容之一，附着探测是未来一段时间内人类探索小天体的主要方式。下降阶段是着陆器附着小天体或完成采样返回探测的关键阶段，对能否安全、准确到达预设的具有科学探测价值的目标区域起着决定性的作用，这对下降阶段的轨迹设计、导航和制导控制都提出了很高的要求。小天体着陆器的下降过程可以转化为一个轨迹优化问题以及对标称轨迹的跟踪控制问题。设定的标称轨迹需要能够安全、准确地到达指定着陆点，满足始末状态约束、路径约束、控制约束等多重约束，同时使某项重要的性能指标最优化，如燃耗、时间等。轨迹优化方法主要包括基于参数化方法的直接法和基于庞特里亚金极小值原理的间接法。直接法无需推导优化问题的横截条件，避免了求解两点边值问题的协态初值敏感困难，因而得到了广泛的应用。此外，对于小天体附近的轨迹优化问题，还需要建立精确描述小天体附近引力加速度的引力场模型。

在已发展的小天体附着探测下降轨迹优化方法中，在先技术[1](参见Lantoine G,Braun R.Optimal trajectories for soft landing on asteroids.Space Systems Design Lab,Georgia Institute of Technology,Atlanta,GA,AE8900 MS Special Problems Report,Dec.2006.)采用多面体模型求解目标小天体附近的引力加速度，以能量作为优化指标，运用直接法进行大步长优化，从而估计得到优化问题的协态初值，然后基于庞特里亚金原理进行间接法轨迹优化计算。多面体模型求引力加速度效率低，优化算法繁琐，耗时长。

在先技术[2](参见Ren,Y.and Shan,J.,“Reliability-Based Soft Landing Trajectory Optimization near Asteroid with Uncertain Gravitational Field,”Journal of Guidance,Control,and Dynamics,Vol.38,No.9,2015,pp.1810-1820.)，首先构建收敛区域更明确的能量最优问题，然后通过调整同伦系数，序列求解，最终将能量最优问题转化为燃料最优问题，采用同伦算法仍然需解两点边值问题，优化过程耗时较长。

技术实现要素：

本发明目的是为了解决现有小天体附着探测下降轨迹优化方法因采用多面体引力场模型求引力加速度效率差；以及采用间接法进行轨迹优化解算，因协态初值难以估计求解困难的问题，提供一种小天体附着探测下降轨迹优化方法。

一种小天体附着探测下降轨迹优化方法，采用内球谐引力场模型估计目标小天体附近的引力加速度，采用凸优化算法解轨迹优化方程。

所述凸优化算法包括松弛、线性化、离散化和内点法。

一种小天体附着探测下降轨迹优化方法，包括以下步骤：

步骤一、由最小二乘法估计内球谐引力场模型的球谐系数：

在小天体外部构造内布里渊球，所构造内布里渊球应该与小天体表面的目标着陆点相切，球心选取应保证着陆器下降轨迹包含在内布里渊球内部。在内布里渊球内部任取N_data个点，由多面体模型计算N_data个点的引力加速度，并由所述引力加速度构造3N_data×1维矩阵由最小二乘法求取内球谐引力模型球谐系数。

其中，n和m分别表示勒让德多项式的阶次和幂次，为一个(n²+2n)×1维的矩阵，包含了所要求取的各阶内球谐系数

为引力加速度关于内球谐系数的导数构成的矩阵，维度为(n²+2n)×3N_data。W为一个3N_data维单位阵；

步骤二、构造小天体附着探测下降轨迹优化方程：

在小天体固连坐标系下，着陆器满足如下的动力学方程

其中，r＝[x,y,z]^T和分别为固连坐标系下着陆器的位置矢量和速度矢量；ω＝[0,0,ω]^T为小天体自旋角速度矢量；T＝[T_x,T_y,T_z]^T为着陆器推力矢量；m_e为着陆器的质量；I_sp为发动机比冲；g_e＝9.807为地球引力加速度常数；为引力加速度矢量，由以下公式计算得到：

其中，G为万有引力常数，M为目标小天体质量，R为步骤一中内布里渊球半径，δ_0,m为Kronecker delta(克劳内克)函数(当m＝0时，δ_0,m＝1)。为步骤一所求的内球谐系数；为步骤一中内球谐引力场模型基函数，满足如下递推关系。

着陆器满足如下的边界条件

其中，r₀，v₀和m₀分别为起始位置处着陆器的位置、速度和质量。t_f为终端时间，r_f为目标着陆点，v_f为终端速度约束，对于软着陆问题，v_f＝0。

着陆器推力矢量满足如下约束条件

0≤||T||≤T_max (7)

燃料最优性能指标表示如下

公式(3)、(6)、(7)、(8)构成了下降轨迹优化方程；

步骤三、向步骤二所得的下降轨迹优化方程中引入松弛变量Γ，对所述下降轨迹优化方程进行松弛：

引入松弛变量Γ替代下降轨迹优化方程中的||T||，则松弛后的轨迹优化方程为：

步骤四、对步骤三所得方程线性化处理：

定义新变量如下

将以上变量代入步骤三的优化方程中，得到线性化的优化方程为

其中，t表示时间变量。

步骤五、对步骤四所得的优化方程进行离散化处理：

将时间区间[0,t_f]等分成N份，得到对步骤四的优化方程进行离散化处理，经过离散化以后，轨迹优化方程转换为参数优化方程，参数优化方程的表达式如下：

其中，M＝[I₆,0_6×1]，t_k表示第k个时间节点，u_k和σ_k分别表示t_k时刻u和σ的取值。

步骤六、采用内点法对步骤五中的参数优化方程进行迭代求解，具体求解过程如下：

1)令引力加速度为一常值▽V₀，采用内点法求解步骤五中的参数优化方程，得到一条轨迹；

2)将步骤1)所得到的轨迹作为参考轨迹，计算参考轨迹中每个节点(即k＝0,…,N)处的引力加速度，并带入步骤五中的参数优化方程中，采用内点法进行求解，得到一条新的轨迹，将得到的新轨迹作为下一次迭代的参考轨迹；

3)当得到的轨迹收敛，则得到最优解，即得到探测最优下降轨迹。

步骤五所述离散化处理方法采用显式四阶龙格库塔积分公式法；

有益效果

本发明所给出的一种小天体附着探测下降轨迹优化方法，采用内球谐引力场模型估计小天体附近引力加速度，具有计算效率高的优点。采用凸优化算法求解燃耗最优控制问题，避免了间接法的推导复杂，以及协态变量无物理意义不易猜测的问题，推导步骤较为简便，节约了计算时间，是一种快速、高精度的估计优化方法，且所得结果符合初始和末端状态约束、动力学约束和控制约束。

附图说明

图1是本发明方法的流程图；

图2是仿真结果示意图，共经过四次迭代得到最有轨迹；其中(a)是小天体固连坐标系下着陆器沿x轴方向的轨迹，(b)是沿y轴方向轨迹示意图，(c)是沿z轴方向轨迹示意图，(d)是推力大小示意图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

实施例1

一种小天体附着探测下降轨迹优化方法，包括以下步骤：

步骤一、由最小二乘法估计内球谐引力场模型的球谐系数：

在小天体外部构造内布里渊球，所构造内布里渊球应该与小天体表面的目标着陆点相切，球心选取应保证着陆器下降轨迹包含在内布里渊球内部。在内布里渊球内部任取N_data个点，由多面体模型计算N_data个点的引力加速度，并由所述引力加速度构造3N_data×1维矩阵由最小二乘法求取内球谐引力模型球谐系数。

其中，n和m分别表示勒让德多项式的阶次和幂次，为一个(n²+2n)×1维的矩阵，包含了所要求取的各阶内球谐系数

为引力加速度关于内球谐系数的导数构成的矩阵，维度为(n²+2n)×3N_data。W为一个3N_data维单位阵；

步骤二、构造小天体附着探测下降轨迹优化方程：

在小天体固连坐标系下，着陆器满足如下的动力学方程

其中，r＝[x,y,z]^T和分别为固连坐标系下着陆器的位置矢量和速度矢量；ω＝[0,0,ω]^T为小天体自旋角速度矢量；T＝[T_x,T_y,T_z]^T为着陆器推力矢量；m_e为着陆器的质量；I_sp为发动机比冲；g_e＝9.807为地球引力加速度常数；为引力加速度矢量，由以下公式计算得到：

其中，G为万有引力常数，M为目标小天体质量，R为步骤一中内布里渊球半径，δ_0,m为(克劳内克)Kronecker delta函数(当m＝0时，δ_0,m＝1)。为步骤一所求的内球谐系数；为步骤一中内球谐引力场模型基函数，满足如下递推关系。

着陆器满足如下的边界条件

其中，r₀，v₀和m₀分别为起始位置处着陆器的位置、速度和质量。t_f为终端时间，r_f为目标着陆点，v_f为终端速度约束，对于软着陆问题，v_f＝0。

着陆器推力矢量满足如下约束条件

0≤||T||≤T_max (19)

燃料最优性能指标表示如下

公式(15)、(18)、(19)、(20)构成了下降轨迹优化方程；

步骤三、向步骤二所得的下降轨迹优化方程中引入松弛变量Γ，对所述下降轨迹优化方程进行松弛：

引入松弛变量Γ替代下降轨迹优化方程中的||T||，则松弛后的轨迹优化方程为：

步骤四、对步骤三所得方程线性化处理：

定义新变量如下

将以上变量代入步骤三的优化方程中，得到线性化的优化方程为

其中，t表示时间变量。

步骤五、对步骤四所得的优化方程进行离散化处理：

将时间区间[0,t_f]等分成N份，可以得到由显式四阶龙格库塔积分公式，对步骤四的优化方程进行离散化

其中，r_k，v_k和p_k分别表示r，v和p在时间节点t_k处的值，U_k和▽V_k分别表示状态变量、控制输入和引力加速度在时间节点t_k处的值。

u_k和σ_k分别表示t_k时刻u和σ的取值。经过离散化以后，轨迹优化方程转换为一个参数优化方程，其表达式如下：

其中，M＝[I₆,0_6×1]。

步骤六、采用内点法对步骤五中的参数优化方程进行迭代求解，具体求解过程如下：

1)令引力加速度为一常值▽V₀，采用内点法求解步骤五中的参数优化方程，得到一条轨迹；

2)将步骤1)所得到的轨迹作为参考轨迹，计算参考轨迹中每个节点(即k＝0,…,N)处的引力加速度，并带入步骤五中的参数优化方程中，采用内点法进行求解，得到一条新的轨迹，将得到的新轨迹作为下一次迭代的参考轨迹；

3)当得到的轨迹收敛，则得到最优解，即得到探测最优下降轨迹。

图2是在小天体216Kleopatra上进行的仿真示意图。着陆器的初始位置为[-150108,6010,-1034]m，目标位置为[-112600,6340,-12520]m。采用本发明所提出的小天体附着探测下降轨迹优化方法，共经过四次迭代，得到最终的燃耗最优轨迹。由图2可知，优化结果满足各项约束条件，着陆器以零速度软着陆于预设着陆点，整个计算过程耗时11.4秒，同时优化得到的控制力始终在设定的发动机推力上限之内。