2016, Vol. 3小行星探测一直是深空探测领域的研究热点，自20世纪90年代开始，美国、日本、欧盟已发射了多个探测器用于对小天体进行探测研究。截至目前，只有3个探测器成功在小天体上实现了着陆，分别是“会合- 舒梅克号(Shoemaker)”(2001，美国)、“隼鸟号(Muses-C)”(2005，日本)、“罗塞塔号(Rosetta)”(2014，欧盟)。美国于2007年发射了“黎明号(Dawn)”探测器，于2011年7月进入灶神星轨道，在经过一年多的绕飞探测之后，于2012年9月离开灶神星，驶向谷神星，并于2015年3月进入谷神星轨道，目前在轨运行。日本于2014年12月发射了“隼鸟2号(Hayabusa 2)”，预计于2018年7月到达目标小天体162173 Ryugu，并计划进行为期一年半的探索和着陆采样后，于2020年返回地球。即将要发射的探测器为OSIRIS-Rex，预计于2016年由NASA发射，任务计划7年，将对目标小天体101955 Bennu进行着陆与采样返回。
相对于地球、火星等近似于圆球状的行星，大多数小天体形状极不规则，所以其附近的引力场极其复杂，这使得小天体轨道动力学、绕飞、着陆导航与制导控制等相关研究变得异常困难。相关学者对不规则小天体的引力场建模已经做了大量研究，目前的方法主要有：球谐函数法[1]、质点群法[2]和多面体模型法[3, 4](为了防止与术语“多面体模型”混用，下面统一使用“多面体法”表示基于多面体模型的引力场建模方法)。球谐函数法通过无穷级数来逼近小行星的引力势函数，但其主要问题在于缺乏轨道数据来确定球谐系数Cnm和Snm，而且其收敛性依赖于位置，当所计算的位置接近于有效半径时，收敛非常缓慢，当所计算的位置位于有效半径(布里渊球)内时，其计算结果会存在很大的误差甚至完全发散[5]。因此该方法不适用于小天体表面着陆制导与导航的相关研究。质点群法使用足够多的体积元(小球体)逼近目标小天体的形状，从而计算出目标天体的引力场。但是，半径为r的球体大约只有边长为2r的立方体的体积的52%，当使用体积元逼近质量连续分布的目标小天体时，大约有48%的体积为空[4, 6]。文献[4]、文献[7]~文献[8]显示，随着计算点接近目标天体表面，质点群法计算误差随之增大。多面体法使用多面体模型来逼近不规则小行星的形状，进而通过积分变换对引力场体积分进行处理，求出任意多面体的引力势能，所以可以用于任意不规则匀质小天体的引力场建模，目前将该方法已广泛用于不规则小天体引力场特性分析及其轨道动力学相关研究[6, 7, 9]。相对于球谐函数法和质点群法等方法，多面体法可以对匀质不规则小天体表面附近引力进行更为精确的计算，因此更适用于不规则小天体着陆导航与制导控制相关研究。文献[10]~文献[12]将其用于小天体着陆燃耗最优控制研究仿真，但是该方法存在计算量大的问题[9, 10]。在小天体着陆导航控制仿真中，需要不断计算不同位置引力加速度，以此进行轨道数值积分计算，此过程将会耗费大量时间，部分文献为了减少计算量，将球谐函数法和多面体法进行了结合[7, 9, 10]，即在有效半径外，使用球谐函数法，当探测器位置小于有效半径时，使用多面体法进行仿真计算。然而，相关文献并没有对多面体法的时间复杂度进行分析，也没有分析不同精度的多面体模型对计算误差造成的影响。对这些问题的研究，有助于确定是否可以选择适合精度的小天体多面体模型，既满足计算精度的要求，又可以提高计算速度，减少相关研究仿真实验计算时间，甚至达到实时计算，使用于半物理仿真中。
因此，本文基于Planetary Data System(PDS)公布的433 Eros多面体模型[8]做了以下工作：使用多面体法计算了433 Eros表面的引力加速度分布情况；分析了不同精度多面体模型对多面体法计算时间和误差造成的影响；确定了以433 Eros为目标天体的着陆导航与制导控制相关仿真计算中，同时满足计算速度和精度要求所需要的多面体模型精度。
1 多面体模型的引力场建模多面体模型可用于描述任意物体的形状，目前广泛应用于计算机三维建模中。多面体模型的引力场建模[3, 4]为
| \[\begin{array}{l} U(r) = \frac{1}{2}G\rho \sum\limits_{e \notin edges} {r_e^{\rm{T}}} {E_e}{r_e} \cdot {L_e} - \\ \quad \quad \quad \frac{1}{2}G\rho \sum\limits_{f \notin facet{\rm{s}}} {r_f^{\rm{T}}} {F_f}{r_f} \cdot {\omega _f}n\\ \end{array}\] | (1) |
引力势能U对r进行一阶求导，得到引力加速度，其表达式为
| \[\begin{array}{l} g(r) = \nabla U(r) = - G\rho \sum\limits_{e \notin edges} {{E_e}{r_e} \cdot {L_e}} + \\ \quad \quad \quad G\rho \sum\limits_{f \notin facet{\rm{s}}} {{F_f}{r_f} \cdot {\omega _f}} \end{array}\] | (2) |
本文首先验证所实现的多面体法计算结果的正确性，并对该方法的误差进行初步分析。考虑到匀质球体的引力加速度有解析解，本文做了如下实验：首先使用3D Max建立了一个半径为16 km三维球模型，模型的面数为39 600，顶点个数为19 802个，将该多面体模型导出为obj格式文件，然后读取模型数据并使用多面体法计算各检验点的引力加速度，最后将计算结果与理论值进行比较，其中，检验点为在半径为16~ 48 km的三维球空间内随机选取的500个点；匀质球体的理论引力加速度计算方式为：$g(r) = - \frac{{GM}}{{{r^2}}}\widehat r$，式中，M为球体质量，r为位置向量，$\widehat r$为其对应的单位向量，r为r的幅值。本文所有计算中，万有引力常量 G = 6.67×10-11 N·m2/kg2。
本文所有实验使用的计算机为Inter® CoreTM 2 Quad CPU Q9400 @2.66GHZ，内存DDR2 800 4 G，操作系统Windows 7 32bit，计算所使用语言为Matlab，版本为R2013a。
定义误差为abs(gp-gt)，其中，gp表示多面体法计算的加速度大小，gt表示理论计算的加速度大小，同时定义相对误差为abs(gp-gt)/gt×100%。表 1为多面体法计算结果的误差统计，其中，gx、gy、gz分别表示x、y、z三个轴向的加速度，g为计算点的加速度矢量幅值。表 1中，最大值、最小值、平均值都是根据理论计算结果进行统计，可以看出，多面体法计算的误差量级在10-4 cm2·s-1以下，为了进一步比较，表 2和图 1分别给出了计算结果的相对误差统计和相对误差直方图，从中可以发现，多面体法计算结果相对与理论计算结果误差大约为0.04%，出现误差的原因是多面体模型并不能对理想球体的形状进行完美逼近。
|
表 1 半径为16～48 km三维球空间内随机500个点的计算结果统计 Table 1 Statistics of calculation results of 500 points random sampling from sphere space with radian from 16～48 km |
|
|
表 2 半径为16～48 km三维球空间内随机500个点的相对误差统计 Table 2 Statistics of relative errorsof 500 points random sampling from sphere space with radian from 16～48 km |
|
图 1 16～48 km三维球空间内随机500个点引力加速度相对误差直方图 Fig. 1 Histogram of gravitational acceleration relative error of 500 points random sampling from sphere space with radian from 16～48 km |
实验所使用的433 Eros多面体模型来源于PDS[8]。其中，Multispectral Imager(MSI)数据提供了6种不同精度的多面体模型，其面数分别为1 708、7 790、10 152、22 540、89 398、200 700，对应的数据文件为eros001708.tab、eros007790.tab、eros010152.tab、eros022540.tab、eros089398.tab、eros200700.tab，模型面数越高，其相应顶点数也越多，代表模型的精度越高，对小天体形状的描述越精确。本文使用多面体法计算了433 Eros表面的引力加速度分布情况，为了提高计算精度，使用面数为200 700的多面体模型，433 Eros物理参数为：平均密度2 670 kg/m3，自转周期为5.27 h。
图 2为使用多面体法计算得到的433 Eros表面引力加速度分布情况，其中，加速度的大小使用不同的颜色表示，单位为cm/s2。从图中可以看出，由于形状不规则，导致433 Eros表面的引力加速度分布极其复杂。根据计算结果，其表面引力加速度最小值约为0.43 cm/s2，最大值约为0.57 cm/s2。图 3为考虑天体自转引起离心加速度后，433 Eros表面重力加速度的分布情况，其最小值约为0.25 cm/s2，最大值约为0.56 cm/s2。
|
图 2 433 Eros表面引力加速度 Fig. 2 Surface gravitational acceleration of 433 Eros |
|
图 3 433 Eros表面重力加速度 Fig. 3 Surface gravitational acceleration of 433 Eros with centrifugal acceleration correction |
为了进一步研究在不同精度的多面体模型下，多面体法的计算精度以及耗费时间，本文做了如下两个实验：第一个实验使用了6个不同精度的多面体模型分别计算了eros001708.tab文件数据的856个顶点的引力加速度，并对计算时间、结果误差进行了统计分析，其中，误差计算以面数最高eros200700.tab文件数据计算结果为基准；第二个实验为使用两个不同精度的多面体模型对探测器的轨道进行数值仿真，并对比其仿真耗时和结果误差。
表 3为不同面数多面体模型计算统计结果，可以看出，随着计算所使用的模型精度不断提高，计算结果的平均误差和平均相对误差都逐渐减小，但是计算所花费的时间有所增大。当模型面数为1 708时，平均每点耗时仅为0.001 7 s，但是其相对误差较高，为 1.077 7%；当模型面数为200 700时，其计算每点所需耗时为0.318 1 s。
图 4为根据表 3的统计结果绘制的每个点计算耗时、相对误差和面数关系双纵轴坐标图，可以看到，每点耗费时间和面数近乎线性关系。通过对式(2)做进一步分析，定义多面体面数f，顶点数v，棱边数e，对于一个封闭的多面体，有e = v + f - 2。若组成多面体的面全部为三角形，则有e = 3f/2，故有v = f/2 + 2。假设加法和乘法各占一个时间单元，ln函数和arctan函数每计算一次时间固定，则式(2)的运行时间开销为 O(e) + O(f) = O(n)，即多面体法的时间开销为线性的，理论分析结果与实验结果一致。
|
图 4 不同面数多面体模型运行时间和相对误差 Fig. 4 The execution time and relative error of polyhedron shapemodelswith different facets |
|
|
表 3 不同面数多面体模型计算结果统计 Table 3 Statistics of calculation results of polyhedron shapemodel with different facets |
根据表 3和图 4可知，当多面体模型面数为89 398 和200 700时，其每点计算耗时分别为面数约为22 540 的模型耗时的4.5倍和10.5倍，但是平均误差下降趋势变小，当模型面数为89 398时，其计算结果相对误差为0.029 3 %，已经接近于面数200 700的模型的计算结果，即当多面体模型精度达到一定程度后，进一步提高模型精度，所取得的引力加速度精度提高有限，但是计算所需时间会继续呈线性增加。值得注意的是，200 700约为89 398的2.2倍，但是两个模型计算结果相对误差只有0.029 3 %，这说明，两个模型对433 Eros 的形状已经有很好的描述，使用面数为200 700模型计算结果作为基准做误差计算，并不影响结论的正确性。
第二个实验做了如下设定：假设探测器初始位置为x = 0 m，y = 30 000 m，z = 0 m，初始速度为vx = 0 m/s，vy = -3.5 m/s，vz = 0 m/s，不考虑天体自转，坐标系定义与模型中定义的坐标系一致，探测器所受的推力为0，在引力的作用下自由运动，轨道时间为105 s，使用Matlab的“ode45”函数(Dormand-Prince方法)，相对容差(relative tolerance)设置为10-7，绝对容差(absolute tolerance)设置为10-8，为了进行对比，使用了面数为 22 540和200 700的模型计算引力加速度，并分别进行两次仿真，积分所得轨道如图 5所示。从图 5可以看出，两条轨道近乎一致，其中，使用面数为22 540和200 700的模型进行仿真时，计算耗时约分别为29.5 s和370.9 s，均远小于轨道时间；但是，后者耗时大约是前者的13倍。仿真结束时，探测器的位置分别为：1) x = 20 077.567 m，y = -5 599.781,z = 104.504 m；2) x = 20 291.950 m，y = -4 463.600 m，z = 112.328 m；3个方向误差分别为214.383 m、1 136.180 m、7.823 m，这说明，随着积分时间的增长，引力计算误差会慢慢累积。根据文献[13]，“会合-舒梅克号(Shoemaker)”的最终着路段时间大约为45 min，对仿真时间为45 min后的轨道误差进行计算，3个方向误差分别为0.087 m、1.977 m、0.000 4 m，均小于2 m，满足≤5 m的精确着陆需求。实验过程中，CPU的负载率在50 %上下波动，考虑到计算所使用CPU为2008年上市，所以在最新计算机上，上述运行时间还可以进一步缩短，这说明：1)使用面数为22 540的Eros 433三维多面体模型进行探测器着陆导航与制导控制相关研究仿真，可以同时满足计算速度和精度的要求，同时仿真计算所需时间大大缩短，此时计算引力加速度每点耗时为0.030 4 s，接近实时要求；2)在实时控制系统中，一般控制周期为20～100 ms，假设每个控制周期只需要计算一次引力加速度，同时在此周期内，探测器所受的引力加速度恒为此值，那么计算每点≤0.030 4 s的耗时满足实时计算要求，所以在半物理仿真实验中使用多面体法是可行的。
|
图 5 使用面数为22 540和200 700的多面体模型对探测器轨道仿真计算 Fig. 5 Orbit simulations of spacecraft with 22 540-facets and 200 700-facets polyhedralmodels |
本文基于多面体法计算了433 Eros表面的引力加速度分布情况，研究了多面体法在不同精度多面体模型下的计算时间和精度问题，实验结果表明：1)多面体法的时间复杂度为O(n)，随着计算所使用的多面体模型精度不断提高，计算结果误差逐渐减小，但是计算所花费的时间呈线性增加；2)当多面体模型精度达到一定程度后，进一步提高模型精度，计算取得的引力加速度精度提高有限，所以在仿真实验中，选择一定精度的多面体模型，可以同时满足计算速度和精度的要求；3)在以433 Eros为目标天体的小天体着陆导航与制导控制相关研究中，使用面数为22 540的多面体模型进行引力加速度计算，可以满足≤5 m的精确着陆需求，同时计算速度接近实时，可以用于半物理仿真实验。
| [1] | Vallado D A, McClain W D. Fundamentals of astrodynamics and applications[M]. Springer Science & Business Media, 2001.( 1) |
| [2] | Geissler P, Petit J M, Durda D D, et al. Erosion and ejecta reaccretion on 243 Ida and its moon[J]. Icarus, 1996, 120(1): 140–157.(1) |
| [3] | Werner R A. The gravitational potential of a homogeneous polyhedron or don't cut corners[J]. Celestial Mechanics and Dynamical Astronomy, 1994, 59(3): 253–278.(2) |
| [4] | Werner R A, Scheeres D J. Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representations of asteroid 4769 Castalia[J]. Celestial Mechanics and Dynamical Astronomy, 1996, 65(3): 313–344.(3) |
| [5] | Hofmann-Wellenhof B, Moritz H. Physical geodesy[M]. Springer Science & Business Media, 2006.(1) |
| [6] | Chanut T G G, Winter O C, Tsuchida M. 3D stability orbits close to 433 Eros using an effective polyhedral model method[J]. Monthly Notices of the Royal Astronomical Society, 2014, 438(3): 56.(2) |
| [7] | Scheeres D J, Ostro S J, Hudson R S, et al. Dynamics of orbits close to asteroid 4179 Toutatis[J]. Icarus, 1998, 132(1): 53–79.(3) |
| [8] | Raugh A C. Near-a-5-collected-models-V1.0 [EB/OL].[2015-10-05]. http://sbn.psi.edu/pds/resource/nearbrowse.html.(3) |
| [9] | Rossi A, Marzari F, Farinella P. Orbital evolution around irregular bodies[J]. Earth, Planets and Space, 1999, 51(11): 1173–1180.(3) |
| [10] | Lantoine G, Braun R. Optimal trajectories for soft landing on asteroids[R]. AE8900 MS Special Problems Report, Space Systems Design Lab, Georgia Institute of Technology, Dec. 2006.(3) |
| [11] | Hawkins M, Guo Y, Wie B. ZEM/ZEV feedback guidance application to fuel-efficient orbital maneuvers around an irregular-shaped asteroid[C]//AIAA Guidance, Control, and Navigation Conference. Minneapolis, Minnesota: AIAA, 2012. |
| [12] | Yang H, Baoyin H. Fuel-optimal control for soft landing on an irregular asteroid[J]. Aerospace and Electronic Systems, IEEE Transactions on, 2015, 51(3): 1688–1697.(1) |
| [13] | Dunham D W, Farquhar R W, McAdams J V, et al. Implementation of the first asteroid landing[J]. Icarus, 2002, 159(2): 433–438.(1) |