一种6-DOF小行星着陆轨迹序列凸优化方法.pdf

第 45 卷第 3 期 2024 年 3 月宇 航 学 报Journal of AstronauticsNo.32024MarchVol.45一种6-DOF小行星着陆轨迹序列凸优化方法初彦峰1，穆荣军1，梁浩2，崔乃刚1（1.哈尔滨工业大学航天学院，哈尔滨 150001；2.北京宇航系统工程研究所，北京 100076）摘要：针对6自由度的小行星动力着陆多约束轨迹优化问题，提出了一种基于序列凸优化的小行星着陆轨迹优化算法。首先采用多面体引力场模型计算小行星的引力场，其可用于描述任意形状的小行星引力场，且比质点群法和球谐函数展开法等方法计算精度更高。为了在凸约束下解决小行星着陆过程燃料消耗最优问题，通过对探测器动力学模型及其约束进行线性化和离散化，将原来的非凸连续时间优化问题转化为凸化的子问题，即二阶锥规划问题（SOCP）；然后引入了虚拟控制项和信赖域，以增强算法的鲁棒性。通过在形状不规则小行星上的着陆模拟，验证了所提出算法的有效性，仿真结果满足各项约束条件，可实现高精度着陆和燃料消耗最优的目标。关键词：小行星着陆；多面体引力场模型；6自由度；序列凸优化中图分类号：V448.2 文献标识码：A 文章编号：1000-1328（2024）03-0341-11 DOI：10.3873/j.issn.1000-1328.2024.03.002A Sequential Convex Optimization Method for 6-DOF Asteroid Landing TrajectoryCHU Yanfeng1，MU Rongjun1，LIANG Hao2，CUI Naigang1（1.School of Astronautics，Harbin Institute of Technology，Harbin 150001，China；2.Beijing Institute of Aerospace Systems Engineering，Beijing 100076，China）Abstract:Aiming at the multi-constraint trajectory optimization problem of 6-degrees-of-freedom（6-DOF）dynamic asteroid landing，an optimization algorithm of asteroid landing trajectory based on sequence convex optimization is proposed.Firstly，the polyhedron gravity field，a suitable approach for describing the gravitational field of irregular-shaped asteroid，is adopted to calculate the gravitational field of the asteroid，which can be used to describe the gravitational field of the asteroid of any shape，and the calculation accuracy is higher than that of the spherical harmonic gravity field and the particle group gravity field.In order to solve the problem of optimal fuel consumption during asteroid landing under convex constraints，the initial non-convex continuous-time optimization problem is transformed into convex sub-problem，namely the second-order cone programming problem（SOCP），by linearizing and discretization of the dynamics model and constraints.Virtual control terms and trust domains are then introduced to enhance the robustness of the algorithm.The effectiveness of the proposed algorithm is demonstrated by the simulation of landing on an irregular-shaped asteroid，and the results meet all the constraints，achieving the goals of high-precision landing and optimal fuel consumption.Key words:Asteroid landing；Polyhedron gravitation model；Six-degree-of-freedom（6-DOF）；Sequential convex optimization收稿日期：2023-09-26；修回日期：2024-01-31基金项目：国家载人航天专项第四批预研项目（18123060201）宇航学报第 45 卷0引言作为太阳系起源的原始物质，小行星蕴含了行星形成及动力学过程的信息，其为研究太阳系演化、深空资源开发、行星防御等空间科学问题提供了重要途径1。小行星探测已成为当代太空争夺的新前沿，也是一个国家综合国力的体现。最初人类只能利用天文望远镜对小行星进行观测，而随着科学技术的发展，自20世纪80年代以来，人类已经实现了对小行星的环绕飞行及采样返回等探测任务2。NASA 发射的伽利略木星探测器于 1991 年对加斯普拉小行星进行了环绕飞行，这是人类第一次实现对小行星近距离探测，并于1993年飞越了编号为243的Ida小行星3；2018年黎明号探测器完成对小行星灶神星（Vesta）和谷神星（Ceres）的探测4；2020年奥西里斯-雷克斯探测器实现了在Bennu小行星上着陆及采样，之后其携带约 250 g 样本于2023年在美国犹他州安全着陆，这是NASA首个小行星采样返回任务5。日本宇宙航空研究开发机构（JAXA）于2014年发射了国际上首颗采样返回式小行星探测器隼鸟2号，并于2020年实现了对小行星龙宫（1999JU3）采集样本及返回的探测任务6。ESA和NASA联合推出了小行星撞击与偏转评 估 计 划，于 2022 年 利 用 DART 航 天 器 撞 击Didymos小行星，以迫使其偏离运行轨道7。我国嫦娥二号月球探测器于2012年对Toutatis小行星进行了交会探测，拍摄了多张高清的小行星光学图像8。小行星引力场建模是探测器实现精准着陆的基础。由于形状不规则，小行星附近的引力场呈现出不规则性。球谐系数引力场模型使用球形势的谐波级数展开，所以该方法对于非球形天体引力场的描述有局限性；内球谐引力场模型的收敛域可以扩展到小天体表面任意一点，但该方法仍无法使用统一表达式实现引力场的全局表征9。质点群法采用有限个质点填充小天体内部空间的方式进行建模，当观察点靠近小天体表面时引力场的精度会迅速退化，因此该方法并不适用于小天体表面的动力学研究10。多面体引力场模型可以精准地描述任意均质多面体的引力势、引力加速度和引力梯度张量，被誉为小天体引力场的优美解10-11。多面体引力场被视作高精准的模型，NASA已将该模型用于433Eros小行星探测任务的规划中12。着陆轨迹对探测器能否安全着陆至关重要。当探测器接近小行星表面时将处于潜在的碰撞危险中，并且探测器将着陆到地形复杂区域以采集到高价值的样本。这都要求探测器的着陆轨迹满足一定的约束条件，同时还需考虑如燃料消耗和着陆时间等性能指标。小行星着陆轨迹优化方法已经成为一个热门的研究问题。文献 13 采用内球谐引力场模型对小行星引力场进行建模，通过变量松弛和离散化处理将非凸的 3-DOF着陆问题转化为更容易求解的等效凸优化问题。文献 14 针对分段状态约束的轨迹优化问题，采用归一化方法，并引入了时间膨胀系数将其转化成为固定时间的轨迹优化问题。文献 15 提出了基于凸优化的轨迹优化方法，通过建立最小着陆误差目标函数，实现了不同飞行时间下太阳辐射压力-控制的轨迹优化，然后采用信赖域约束得到一系列 SOCP 问题。文献 16 在Lambert求解器的基础上提出了一种变飞行时间轨迹优化算法，用于解决小行星着陆过程中的燃料消耗问题。上述方法都是基于3-DOF的制导算法，即只涉及探测器质心的平移运动，而不考虑探测器姿态变化。3-DOF动力着陆制导起源于二十世纪的阿波罗登月计划17，其实现相对容易。对于平动和转动紧密耦合的未来小行星探测任务，3-DOF制导方案可能无法满足探测器的任务需求，实现位置和姿态同步控制的6-DOF制导算法更具优势。轨迹优化方法包括直接法和间接法。间接法采用 Pontryagin 原理求解，并高度依赖于初始估计的准确性，所以使用间接法求解比较困难18。直接法将轨迹优化问题转化为参数优化问题，进而使用非线性规划对参数优化问题进行求解。直接法主要包括伪谱法、凸优化方法和打靶法，其中凸优化方法具有计算速度快和实时性高的特点，被广泛地应用于飞行器轨迹优化问题19。本文提出一种 6-DOF 小行星着陆轨迹序列凸优化方法。首先使用多面体引力场模型建立不规则小行星的引力场；通过对6-DOF动力学模型及其约束进行线性化和离散处理，将原着陆轨迹优化问题转化为二阶锥规划问题，然后通过虚拟控制和信赖域增强算法的鲁棒性；最后通过模拟着陆433Eros342第 3 期初彦峰等：一种6-DOF小行星着陆轨迹序列凸优化方法小行星验证了算法的可行性和有效性。1小行星着陆动力学模型1.1多面体引力场模型多面体引力场模型适合于描述形状不规则小天体的引力场，通过一系列三维实体来逼近小行星形状。Werner等20利用Gauss散度定理推导出了多面体的引力场模型，由面引力和边引力两部分组成。不规则小行星附近空间任意一点的引力加速度可表示为：U(r)=Gf=1nf()Ff rf f-Ge=1ne()Ee re Le（1）式中：U(r)是引力加速度；G为引力常数；为小行星密度；nf是多面体的面数；ne是多面体的边数；rf和re分别是探测器到第f个面（三角形）和第e条边的向量；Ff和Ee是和多面体形状相关的并矢对称矩阵；f和Le是描述计算点与多面体面和边相对关系的无量纲因子10，16。1.2小行星着陆动力学模型本文涉及2种与探测器运动相关的坐标系，如图1所示。小行星固联坐标系OI XIYIZI（简写为OI系）：坐标原点位于小行星的质心，OIZI轴为小行星自转轴，OIXI轴为小行星最大惯性主轴，OIYI与OIXI，OIZI构成右手坐标系。探测器本体坐标系OB XBYBZB（简写为OB系）：坐标原点位于探测器的质心，OBXB，OBYB和OBZB分别为探测器的惯性主轴。探测器的质量消耗计算如下：m(t)=TB()t（2）式中：m(t)是探测器质量；TB(t)是探测器在OB系下的推力；=-1()Ispg为质量流量系数，Isp和g分别为推力器比冲和地球标准重力常数。探测器的动力学方程可以表示为：rI()t=vI()tvI()t=-2evI()t-()e2rI()t+TI()tm()t+U()rqB I()t=12B()t-CBI()qB I()teqB I()tJ B()t=()t-B()t JB()t（3）式中：rI是OI系下探测器位置向量；vI是速度向量；qB I=q0 q1 q2 q3T是探测器姿态四元数；J为探测器的惯性张量；e=ex ey ezT是小行星自转角速度；e=0 -ez ey ez 0 -ex -ey ex 0 是e的反对称矩阵；B和分别为探测器的角速度和力矩；TI是探测器的推力，TB和TI满足TI=CIB(qB I)TB，其中CIB是描述OB系到OI系的方向余弦矩阵。方向余弦矩阵CIB和矩阵()定义分别如下：CBI()qB I=1-2()q22+q232(q1q2+q0q3)2(q1q3-q0q2)2(q1q2-q0q3)1-2()q21+q232(q2q3+q0q1)2(q1q3+q0q2)2(q2q3-q0q1)1-2()q21+q22（4）()=0 -1 -2 -31 0 3 -22 -3 0 13 2 -1 0（5）式中：=1 2 3T。四元数和欧拉角（偏航角、俯仰角和滚转角）的转换关系定义如下：=arctan2()q1q2+q0q3q20+q21-q22-q23arcsin()2()q0q2-q1q3arctan2()q0q1+q2q3q20-q21-q22+q23（6）XIOIYIYBZBOBXBrIZI图1坐标系的定义Fig.1Definition of coordinate systems343宇航学报第 45 卷1.3状态约束和控制约束首先探测器的质量m(t)应始终大于它的干重mdry，质量约束如下：mdry m(t)（7）输出推力T(t)应介于最大推力Tmax和最小推力Tmin之间，推力约束如下：TminT()t Tmax（8）输出力矩M(t)小于最大力矩Mmax即可，约束如下：M()t Mmax（9）考虑到小行星形状不规则，并且存在着自转，探测器下降时可能会发生硬着陆现象。为了防止探测器撞击到小行星，需要对探测器的飞行路径进行约束。从距离小行星较远的位置到着陆点上方，探测器主要在小行星的外部空间飞行，因此形成包含小行星的椭球，也就是说椭球体积应大于小行星的实际尺寸。如图2所示，最外层黄色的网格代表椭球，最里面灰色的是小行星（本文绘制的小行星是433Eros，其 半 长 轴 为16 km 6.5 km 6.5 km）。椭球以外的点满足x2a2+y2b2+z2c2 1，整理后得到的椭球约束数学表达式如下：1 rTIRrI=rTI 1a2 0 0 0 1b2 0 0 0 1c2rI，t 0，tc（10）式中：a，b，c是椭球的半长轴，由其构成的椭球要想将 433Eros 小行星包围，需要满足a 16 km，b 6.5 km，c 6.5 km；rI=x y zT是探测器位置；R是椭球半长轴形成的对角矩阵；tc为动力下降段到最终着陆段的时间。当探测器进入到最终着陆段时，探测器与着陆点的距离足够近，则椭球约束不再适用。要实现高精度着陆，在最终着陆段需要确保探测器配备的相机能始终对着陆点进行监测。相机视野范围应该包括着陆点，也就是说相机视野范围大于相机光轴和视线向量之间的夹角，所形成的视场角约束如图3所示，定义如下：cos-CBI()rI-rf-BTdB-CBI()rI-rf-BdB，t(tc，tf（11）式中：是相机视场角；B是相机在探测器上的安装位置；dB是相机光轴向量；tf为最终着陆时间；rI和rf分别为OI系下探测器当前位置和最终着陆位置。2小行星着陆序列凸优化算法2.1状态方程的凸化系统状态方程式（3）及式（10）、式（11）等都是非凸的，可对其进行线性化，进而获得凸化的表达式。定义探测器的状态向量为x=rT vT T qTT R13，控制变量为u=TT MTT R6。系统状态方程可表示为：f(x(t)，u(t)=dx()tdt=r v qT（12）式（12）的线性化形式为：x(t)=A(x，u)x(t)+B(x，u)u(t)+c(x，u)（13）式中相关的系数矩阵如下：图2椭球约束Fig.2The ellipsoid constraintXIYIrIrIrfZIrfdBOBYBZBXBB CIB(rIrf)B,dBOI图3视场角约束Fig.3The field-of-view constraint344第 3 期初彦峰等：一种6-DOF小行星着陆轨迹序列凸优化方法|A()x，u=f()x，uxx，u R13 13|B()x，u=f()x，uux，u R13 6c()x，u=f()x，u-A()x，u x-B()x，u u R13（14）将探测器的状态向量分为平动状态xp=r vT和转动状态xr=qT，则式（13）可以分解为2部分，如下：xp()t=Ap()xp，u pxp()t+Bp()xp，u pup()t+cp()xp，u pxr()t=Ar()xr，u rxr()t+Br()xr，u rur()t+cr()xr，u r（15）其中相关的系数矩阵如式（16）（17）所示Ap()xp，u p=03 3 E3-()e2-E3r3 -2eBp()xp，u p=03 3CIBmcp()xp，u p=03 1U()r+E3r3r（16）式中：E3代表三阶单位对角矩阵。Ar()xr，u r=Cf()x，u B 03 4Cqf()x，u B Cqf()x，u qBr()xr，u r=J-104 3cr()xr，u r=CCqf()x，u-Ar()xr，u rxr-Br()xr，u ru r（17）式中：C=03 6，E3，03 4，Cq=04 9，E4。根据式（15）（17）可得式（13）中相关矩阵的具体形式如下：A()x，u=Ap()xp，u p 06 7 07 6 Ar()xr，u rB()x，u=Bp()xp，u p 06 3 07 3 Br()xr，u rc()x，u=cp()xp，u pcr()xr，u r（18）2.2约束的凸化引入松弛因子=x y zT代替控制力T=x y zT，两者关系可以定义为j=Tj|TjTj，j=x，y，z（Tj代表估计值），则有：Tmin|j Tmax（19）对椭球约束（10）进行一阶泰勒展开，则有：1+xT(t)Rx(t)-2xT(t)Rx(t)0（20）式中：x(t)是状态x(t)的估计值。同理，对视场角约束（11）进行一阶泰勒展开，则有：h(x)x(t)-h(x)x(t)+h(x)0（21）式中：h(x)计算如下：h(x)=CBI()x()1 3-rf+B+dTBcosdB(CBI(x()1 3-rf)+B)（22）2.3离散化将飞行时间0，tf离散化为N个子区间，设每个区间长度为t，则有：tk=kt，tf=Nt，k=0，1，2，N（23）系统的状态方程离散化后可以表示为：xi(k+1)=Ai(k)xi(k)+Bi(k)ui(k)+c i(k)（24）式中：k K，K=0，1，2，N-1，i代表迭代次数，相关系数矩阵定义如下：Ai()k=A()tk+1，tk=etA()xi-1，ui-1Bi()k=tktk+1A()，tkB()xi-1，ui-1dc i()k=tktk+1A()，tkc()xi-1，ui-1d（25）当k Kc，Kc=0，1，2，Nc，Nc=tc t，椭球约束离散化后：1+()xi-1()kTRxi-1(k)-2()xi-1()kTRx(k)0（26）345宇航学报第 45 卷当k Kf，Kf=Nc+1，N，视场角约束离散化后：h(xi-1(k)(xi(k)-xi-1(k)+h(xi-1(k)0（27）控制力约束离散化后可以表示为：ij(k)=Ti-1j()k|Ti-1j()kTij(k)，j=x，y，z，k K（28）2.4虚拟控制和信赖域由于对系统方程和约束进行了线性化，这样的做法将存在人为不可行性，对于飞行时间较短的轨迹优化问题可能存在不可行的解21。通过在线性化方程中引入虚拟控制项Vx可以解决该类问题，则式（24）可表示为：xi(k+1)=Ai(k)xi(k)+Bi(k)ui(k)+c i(k)+Vix(k)（29）椭球约束和视场角约束同样引入Ve和Vf，则式（26）和式（27）改写为：1+xT(k)Rx(k)-2xT(k)Rx(k)Vie(k)（30）h(xi-1(k)(x(k)-xi-1(k)+h(xi-1(k)Vif(k)（31）式中：Vie(k)0，Vif(k)0。将Vx，Ve和Vf加入到最小燃料消耗的目标函数中，则有：min Ji=k=0N-1i()k1t+Wxk=0N-1Vix()k1+Wek=0NcVie()k1+Wfk=Nc+1NVif()k1（32）式中：Wx 0，We 0和Wf 0为相应的权重系数；在求解过程中，虚拟控制项可视作强制惩罚项以保证算法的收敛性。为确保线性化后能够捕获原始问题的非线性，可以引入信赖域约束使前后两步迭代之间没有明显的偏离22。在燃料最优的目标函数中，状态方程和约束的非线性主要是和xi(k)相关的，因此只引入状态信赖域约束即可。设第i次和i-1次的状态偏差为：xik=xi(k)-xi-1(k)（33）与状态相关的信赖域约束可以表示为：xikTxik ik，k K（34）上述约束引入到目标函数后，式（32）可进一步改写为：min Ji=i1t+WxVix1+WeVie1+WfVif1+Wi1（35）式中：W是信赖域约束的权重系数。2.5序列凸优化算法经过上述凸化、离散化以及引入虚拟控制和信赖域后，小行星着陆轨迹凸优化子问题，即二阶锥规划问题（SOCP）如下所示：min Ji=i1t+WxVix1+WeVie1+WfVif1+Wi1s.t.xi()k+1=Ai()k xi()k+Bi()k ui()k+c i()k+Vi()k，k Kmi()k+1=mi()k+i()kt，k K1+xT()k Rx()k-2xT()k Rx()k Vie()k0 Vie()k，k Kch()xi-1()k()x()k-xi-1()k+h()xi-1()k Vif()k0 Vif()k，k KfxikTxik ik，k KTmin|j Tmax，j=x，y，z，M Mmax，k Kmi()0=mwet，mi()k mdry，xi()0=x0|xin()N-xf，n ef，n，n=1，2，13（36）式中：x0和xf分别表示初始状态和预计着陆状态；ef R13表示允许着陆误差；mwet和mdry分别表示探测器初始质量和探测器的净质量。设最小状态偏差为，最大迭代次数为Niter，初始参考着陆参数表示为x0(k)，u0(k)，m0(k)，则小行星着陆轨迹序列凸优化算法流程如图4所示。首先利用初始状态x0和预计着陆状态xf计算初始参考着陆参数，然后利用其计算系数矩阵。在i小于最大迭代次数时，利用第i-1次迭代的结果进行凸优化求解，该过程需满足状态约束、控制约束及边界条件，并且根据探测器位置矢量ri-1(k)实时更新引力加速度U(ri-1(k)；优化后可得到第i次迭代结果xi(k)，ui(k)，mi(k)，将第i-1次结果和第i次结果进行比较，两者小于一定范围（最小状态偏差），可认为已经收敛并输出优化轨迹；否则以此次计算结果重新执行凸优化。以此类推，直到状态偏差满足最小状态偏差或者达到最大迭代次数。346第 3 期初彦峰等：一种6-DOF小行星着陆轨迹序列凸优化方法3仿真验证为了快速进行序列求解，初始轨迹可利用初始状态x0和预计着陆状态xf均匀插值获取，如下：x0k=N-kNx0+kNxf，k K（37）燃料消耗会导致惯性张量发生变化，定义惯性张量更新方式，如下：Ji(k)=2.10001.970001.41mi(k)（38）式中：Ji(k)代表第i次迭代k时刻的惯性张量；mi(k)是第i次迭代k时刻探测器质量。为了校验所提出算法的有效性，本文以433Eros小行星为典型分析对象，它的自转角速度和密度分别为3.3110-4 rad/s和2.67103 kg/m3。仿真采用了856个顶点和1 708个面的小行星多面体模型，相关数据可以在PDS网站上获取23。时间参数设定为动力下降段至最终着陆段历时tc=900 s及最终着陆时间tf=1 200 s，允许状态偏差为=10-3，其他参数如表1和表2所示。本文数值模拟使用了 MATLAB 软件下的 CVX工具箱进行解算。根据上述模型和算法，小行星着陆轨迹序列凸优化仿真结果如图512所示；为了直观地表示探测器姿态变化，本文将四元数转换为欧拉角（图8）。如图58所示，所提出的算法可以引导探测器到达目标着陆点，探测器的最终着陆状态误差小于允许着陆误差ef，利用其优化出来的轨迹是有效的。开始初始化x0(k),u0(k),m0(k)计算系数矩阵A0(k),B0(k),c0(k)结束i=1iNiter否是是否利用第i1次迭代的结果解决小行星着陆轨迹SOCP问题多面体引力场面引力计算边引力计算i=i+1状态约束椭圆约束视场约束边界条件xi(0)=x0mi(0)=mwetmi(k)=mdry|xin(N)xf,n|ef,n控制约束力约束力矩约束计算得到xi(k),ui(k),mi(k)计算状态偏差xik=xi(k)xi1(k)|xik|xi(k),ui(k),mi(k)Ai(k),Bi(k),ci(k)subject to图4小行星着陆轨迹序列凸优化算法流程Fig.4Flow chart of sequence convex optimization algorithm for asteroid landing trajectory表1仿真参数Table 1Simulation parameters参数探测器初始质量/kg探测器净质量/kg最大推力/N最小推力/N最大力矩/（Nm）相机位置/m相机光轴向量/m探测器视场角/（）推力器比冲/s椭球约束的半长轴/km标称地球重力加速度/（ms-2）万有引力常量/（Nm2kg-2）初始惯性张量/（kgm2）离散区间长度/s最大迭代次数数值1 4001 0002550.5B=0.9，0，-1TdB=0，0，-1T2522522，10.5，7.5T9.806 656.6710-112 940 0 0；0 2 758 0；0 0 1 9741015表2初始状态、预计着陆状态及容许着陆误差Table 2Initial state，final landing state and tolerating landing errors参数r0/mv0/（ms-1）0/（rads-1）q0rf/mvf/（ms-1）f/（rads-1）qfef数值7 143.78，-6 020.65，-8 475.25T1.22，1.43，-0.42T0，0，0T-0.292 504 233 760 79，0.715 419 718 821 09，0.607 314 365 739 738，0.183 807 400 068 947T6 825.68，-4 665.87，-4 533.93T0，0，0T0，0，0T0.175 915 217 312 454，0.730 948 899 081 422，0.600 657 493 003 202，-0.271 989 189 764 371T1，1，1，0.02，0.02，0.02，0.01，0.01，0.01，0.01，0.01，0.01，0.005，0.005，0.005，0.005T347宇航学报第 45 卷从图910可知，探测器的推力和力矩均满足控制约束。探测器的质量变化如图11所示，最终剩余质量为1 394.8 kg，也就是说消耗燃料5.2 kg，满足质量约束。式（10）表明椭球约束与探测器的位置相关，进而影响探测器着陆过程中的速度等参数；经计算，探测器的位置在动力下降段0，tc满足椭球约束。式（11）表明视场角约束的作用是限制探测器姿态运动，即探测器要持续调整姿态才可以追踪着陆点；如图12所示，探测器的视场角在最终着陆段(tc，tf满足视场角约束。为了进一步验证轨迹优化算法对初始状态偏差的鲁棒性，本文进行了500次蒙特卡洛打靶仿真。9640103200400600时间/s(c)rz曲线rz/m8001 0001 2000864103200400600时间/s(b)ry曲线ry/m8001 0001 2000678103200400600时间/s(a)rx曲线rx/m8001 0001 200图5位置变化曲线Fig.5Curves of position0202200400600时间/s(b)vy曲线8001 0001 200vy/(ms1)0202200400600时间/s(a)vx曲线vx/(ms1)8001 0001 2002020200400600时间/s(c)vz曲线8001 0001 200vz/(ms1)图6速度变化曲线Fig.6Curves of velocity02.502.5104200400600时间/s(a)x曲线x/(rads1)8001 0001 2000202103200400600时间/s(b)y曲线y/(rads1)8001 0001 2000101104200400600时间/s(c)z曲线z/(rads1)8001 0001 200图7角速度变化曲线Fig.7Curves of angular velocity01062200400600时间/s(b)曲线8001 0001 200/()01021009896200400600时间/s(a)曲线/()8001 0001 20018001800200400600时间/s(c)曲线8001 0001 200/()图8姿态变化曲线Fig.8Curves of attitude348第 3 期初彦峰等：一种6-DOF小行星着陆轨迹序列凸优化方法初始位置和速度参数如表3所示，其他未列举的参数与表12中相同。着陆误差统计如表4所示，各个参数的单位与表2一致，其中最终着陆误差都达到了边界值（允许着陆误差ef）。由仿真结果可知，即使存在初始状态偏差，在制导算法作用下，探测器仍能以较小的状态误差到达着陆点，表明所提出的制导算法有很强的鲁棒性。0200400600时间/s视场角/()8001 0001 200705060403020100图12视场角变化曲线Fig.12Curve of angle of field of view表3蒙特卡洛打靶仿真初始位置和速度Table 3Initial positions and velocities for Monte Carlo simulation参数rx/mry/mrz/mvx/（ms-1）vy/（ms-1）vz/（ms-1）最小值6 500-6 500-9 000-2-2-2最大值7 700-5 500-8 000222015015200400600时间/s(b)Fy曲线8001 0001 200Fy/N150150200400600时间/s(c)Fz曲线8001 0001 200Fz/N015015200400600时间/s(a)Fx曲线Fx/N8001 0001 200图9推力变化曲线Fig.9Curves of thrust1010103200400600时间/s(c)Mz曲线Mz/(Nm)8001 0001 2000202102200400600时间/s(b)My曲线My/(Nm)8001 0001 2000202103200400600时间/s(a)Mx曲线Mx/(Nm)8001 0001 200图10力矩变化曲线Fig.10Curves of moment0200400600时间/s质量/kg8001 0001 2001 4001 3991 3981 3971 3961 3951 394图11质量变化曲线Fig.11Curve of mass表4着陆状态误差统计Table 4Error statistics for landing state参数均值最大值rx0.621ry0.711rz0.551vx0.020.02vy0.020.02vz0.020.02x0.0050.01y0.0050.01z0.0050.01q00.0010.005q10.0010.005q20.0010.005q30.0010.005349宇航学报第 45 卷4结论本文提出了一种 6-DOF 小行星着陆轨迹序列凸优化算法。通过线性化和离散化将小行星着陆轨迹优化问题转化为一个可以进行序列凸优化求解的SOCP问题；并引入虚拟控制项和信赖域，以保证线性化后采用凸优化方法是可行的。以形状不规则的 433Eros小行星为着陆目标，采用多面体引力场模型对该小行星进行引力场建模；仿真结果表明探测器能在节省燃料的同时进行高精度着陆，并满足各项约束。所提出小行星着陆轨迹优化算法可以扩展到6自由度耦合的其他行星着陆问题，如火星或月球探测；同时，该方法精度高、鲁棒性强等优势对机载应用非常有吸引力。参 考 文 献1 BIERHAUS E B，CLARK B C，HARRIS J W，et al.The OSIRIS-REx spacecraft and the touch-and-go sample acquisition mechanism（TAGSAM）J.Space Science Reviews，2018，214（7）：107.2 TAKAO Y，MORI O，MATSUMOTO J，et al.Sample return system of OKEANOS：The solar power sail for Jupiter Trojan exploration J.Acta Astronautica，2023，213：121-137.3 VEVERKA J，THOMAS P C，HELFENSTEIN P，et al.Dactyl：Galileo observations of Idas satellite J.Icarus，1996，120（1）：200-211.4 苏炳志.小行星探测视觉辅助导航系统滤波方法研究 D.哈尔滨：哈尔滨工业大学，2020.SU Bingzhi.Research on filtering method for vision-aided navigation system in asteroid explorationD.Harbin：Harbin Institute of Technology，2020.5 OBA Y，TAKANO Y，DWORKIN J P，et al.Ryugu asteroid sample return provides a natural laboratory for primordial chemical evolution J.Nature Communications，2023，14（1）：3107.6 KIKUCHI S，MIMASU Y，TAKEI Y，et al.Preliminary design of the Hayabusa2 extended mission to the fast-rotating asteroid 1998 KY26 J.Acta Astronautica，2023，211：295-315.7 TRGOLO N，BAGATIN A C，MORENO F，et al.Lifted particles from the fast spinning primary of the near-earth asteroid（65803）Didymos J.Icarus，2023，397：115521.8 张莹莹.小行星自主着陆轨迹规划和控制方法研究 D.哈尔滨：哈尔滨工业大学，2021.ZHANG Yingying.Trajectory planning and control method for autonomous asteroid landingD.Harbin：Harbin Institute of Technology，2021.9 肖尧，阮晓钢，魏若岩.基于433 Eros的多面体引力模型精度与运行时间研究 J.深空探测学报，2016，3（1）：41-46.XIAO Yao，RUAN Xiaogang，WEI Ruoyan.Research on t