2016, Vol. 3
20世纪90年代以来，随着科学技术水平的不断提高，人类的目光已经不再局限于太阳系内的八大行星，深空探测迎来了空前的发展机遇。在科学研究方面，小行星研究作为21世纪深空探测重要主题之一，已然引起人们的广泛关注。据不完全统计，近10年来，Science杂志以asteroid为主题词搜索，有超过110篇原创论文(research)，Nature杂志更是超过有136篇原创论文，279篇新闻(news)，小行星的研究已经空前深入。在工程任务方面，主要有美国的“近地小行星探测”任务、“深度撞击”任务，日本的“隼鸟号”任务以及欧空局的“罗塞塔”任务等等。中国的“嫦娥2号”卫星也曾于700万千米深空飞越“图塔蒂斯”号小行星，开启了中国人的小行星探测记录。
相较于近球形的大行星和天然卫星，小行星的质量较小而且形状不规则，这使得它们的动力学特性也更加复杂，并且存在着巨大的个体差异。如何建立模型来真实全面地描述小行星的引力场，是研究小行星附近动力学环境的基础。常用的引力场描述方法主要分为两大类。第一类是级数逼近法，此类方法主要是通过无穷级数来逼近引力势函数，利用级数的规律来反映引力的规律，进而完成引力场建模。比较典型的有球谐函数法和椭球谐函数法。第二类方法是模型逼近法。该方法使用可以较为精确地求解势能的模型从形状上来逼近所研究的小行星，进而可以用模型的势函数来代替所研究小行星的势函数作为求解对象。这类方法主要有三轴椭球方法和多面体模型方法[1]。考虑到收敛性和准确性，基于多面体模型的引力场建模方法被认为是目前较为先进的小行星研究方法，自20世纪90年代以来得到长足发展。小行星的多面体模型划分如图 1所示。Werner于1993年提出了基于多面体模型的均质小行星引力场求解方法[2]。Miller等(2002)提出了应用地面天文观测数据确定小行星形状、密度和运动参数的方法[3]。Scheeres等于2004年利用雷达数据研究了小行星25143 Itokawa附近的动力学环境，并发现了其本体系下的4个平衡点[4]。Wang等(2014)利用均质多面体模型计算了一系列小行星的引力场，并研究了每颗小行星的平衡点、稳定性及其拓扑结构分类[5]。Fahnestock和Scheeres (2006)提出了基于双多面体模型的全二体运动仿真方法，第一次将多面体方法应用于多小行星系统的研究中[6]。张振江等(2010)利用均质多面体模型计算了小行星433 Eros附近的引力场数据作为虚拟测量量，以此反解出小行星433 Eros的球谐系数。通过对比该法求得的球谐系数与美国“近地小行星探测”任务中测得的球谐系数，证实了小行星433 Eros并非均质[7]。Scheeres等(2000)就利用多面体模型和球谐函数相结合的方法，给出了小行星433 Eros的一种可能的非均质模型[8]。
|
图 1 小行星的多面体模型划分 Fig. 1 Polyhedron models of the asteroids |
前文中曾提到，Fanhnestock和Scheeres(2006)提出基于双多面体模型的全二体运动仿真方法，该方法的基础是求取多面体间的相互势能及其偏导数。该方法先将多面体分别离散为若干四面体单体，进而把求解两个多面体间相互势能的问题转化为求解分别取自两个多面体的两个四面体单体间的相互势能并二重求和的问题，因此适用于非均质多面体的情况。不难想到，如果把其中一个多面体的体积缩小直至其成为一个质点，求解两个多面体间相互势能的方法就可以用来求解单个多面体的引力场。与双多面体模型类似，该方法可以考虑非均质多面体的情况，进而可以实现非均质小行星的引力场建模。但是该方法的存在截断误差，对于小行星形状极不规则或者距离小行星极近的情况，该方法精度有限。
本文的前两节分别从理论上给出了求解非均质小型附近引力场和平衡点的方法。在第3节中，验证了此方法的准确性，评估了小行星可能存在的非均质特性对小行星附近动力学环境造成的影响，希望借此引起读者对于非均质引力场模型的重视。最后，对后续工作做出了展望。
1 非均质小行星引力场的确定 1.1 小行星的多面体模型划分图 1给出了小行星多面体模型划分的样例。根据观测数据，可以将小行星表面划分为若干三角形单元。可以根据小行星形状的复杂程度和目标精度来决定三角形网格的细密程度。表面单元划分完成后，便可以使用多面体代替小行星作为研究对象。将多面体表面三角形的顶点分别与多面体质心相连，可以构成与表面三角形个数相当的若干四面体单体。求解整个多面体引力场的问题便可以转换为求解各个单体在该点的引力场并线性叠加的问题。
1.2 多面体引力场的表达式整个多面体对于多面体外部任意一点的引力势能为
| $ {U_A} = - G\iiint\limits_A {{\rho _A}\frac{1}{{{r}}}} {\rm{d}}a = - G\sum\limits_{a \in A} {{\rho _a}} \iiint\limits_a {\frac{1}{{{r}}}} {\rm{d}}a $ | (1) |
其中：a表示多面体A中的任意一个四面体单体；ρa代表其密度；r表示单体内任意一点到目标质点的矢量。不难看出这是一个比较复杂的三重积分问题，要解决这个问题，可以利用勒让德级数把被积函数展开，把任意形状的四面体单体变换为有三条棱单位正交的标准四面体单体，并作为积分区域进行三重积分。具体的过程可以参考文献[6]，在文献[9]中有更为详细的叙述。与上述两篇参考文献相比，本文有几点不同。其一，本文将上述文献中研究的双多面体模型简化为多面体-质点模型，所以与多面体B相关的尺寸变量都退化为0；其二，出于精度方面的考虑，本文保留了被积函数的0~4阶级数展开项；其三，本文假设每个四面体单体拥有独立的密度，考虑了非均质情况。综合上述3点，需要将文献中的公式做出一些调整，得到最终的势能表达式为
| $ U = - G\sum\limits_{a \in A}^{} {{\rho _a}\det ({T_a})\\ \left\{ \begin{array}{l} \left[{\frac{{{{{Q}}_0}}}{R}} \right] + \left[{ - \frac{{{{{Q}}_i}{{{w}}^i}}}{{{R^3}}}} \right] + \left[{ - \frac{{{{{Q}}_{ij}}{{{r}}^{ij}}}}{{2{R^3}}} + \frac{{3{{{Q}}_{ij}}{{{w}}^i}{{{w}}^j}}}{{2{R^5}}}} \right] + \left[{\frac{{3{{{Q}}_{ijk}}{{{r}}^{ij}}{{{w}}^k}}}{{2{R^5}}} - \frac{{5{{{Q}}_{ijk}}{{{w}}^i}{{{w}}^j}{{{w}}^k}}}{{2{R^7}}}} \right]\\ \quad \quad \quad \quad + \left[{\frac{{3{{{Q}}_{ijkl}}{{{r}}^{ij}}{{{r}}^{kl}}}}{{8{R^5}}} - \frac{{15{{{Q}}_{ijkl}}{{{r}}^{ij}}{{{w}}^k}{{{w}}^l}}}{{4{R^7}}} + \frac{{35{{{Q}}_{ijkl}}{{{w}}^i}{{{w}}^j}{{{w}}^k}{{{w}}^l}}}{{8{R^9}}}} \right] \cdots \end{array} \right\}} $ | (2) |
式(2)使用了张量符号，但是张量的上、下脚标仅被用来区分其内、外积运算，并无其他含义。式(2)中符号的定义式为
| $ {{{T}}_a} = \left[{\begin{array}{*{20}{c}} {\Delta x_1^a} & {\Delta x_2^a} & {\Delta x_3^a}\\ {\Delta y_1^a} & {\Delta y_2^a} & {\Delta y_3^a}\\ {\Delta z_1^a} & {\Delta z_2^a} & {\Delta z_3^a} \end{array}} \right] ,\quad {{v}}_j^i = \left[\begin{array}{l} {{{x}}^i}\\ {{{y}}^i}\\ {{{z}}^i} \end{array} \right] = \\ \quad\quad - {{{T}}_a} ,\quad {{{w}}^i} = {{{R}}^j}{{v}}_j^i ,\quad {{{r}}^{ij}} = {{v}}_p^i{{v}}_p^j = {{{x}}^i}{{{x}}^j} + {{{y}}^i}{{{y}}^j} + {{{z}}^i}{{{z}}^j} $ | (3) |
其中：$\Delta x_n^a$表示在多面体的本体系下表示的四面体单体a的第n个表面顶点的x坐标，Ta中的其他项同理；R表示在多面体本体系下表示的从多面体质心指向多面体外目标质点的矢量；Q0~Qijkl分别代表了一组0~4阶张量，对于多面体的所有单体而言，这些张量完全相同，且在积分过程中不会发生变化，这里直接给出其表达式
$ {{Q}_0} = \frac{1}{6} , $
$ {{Q}_i} = \frac{1}{{24}}\left[{\begin{array}{*{20}{c}} 1 \! & \! 1 \! & \! 1 \end{array}} \right] ,\quad {{Q}_{ij}} = \frac{1}{{120}}\left[{\begin{array}{*{20}{c}} 2 \! & \! 1 \! & \! 1\\ 1 \! & \! 2 \! & \! 1\\ 1 \! & \! 1 \! & \! 2 \end{array}} \right] , $
$ {{Q}_{ijk}} = \frac{1}{{720}}\left( {\left[{\begin{array}{*{20}{c}} 6 \! & \! 2 \! & \! 2\\ 2 \! & \! 2 \! & \! 1\\ 2 \! & \! 1 \! & \! 2 \end{array}} \right]\left[{\begin{array}{*{20}{c}} 2 \! & \! 2 \! & \! 1\\ 2 \! & \! 6 \! & \! 2\\ 1 \! & \! 2 \! & \! 2 \end{array}} \right]\left[{\begin{array}{*{20}{c}} 2 \! & \! 1 \! & \! 2\\ 1 \! & \! 2 \! & \! 2\\ 2 \! & \! 2 \! & \! 6 \end{array}} \right]} \right) , $
$ {{Q}_{ijkl}} \!=\! \frac{1}{{1260}} \!\! \left( \!\!\! \begin{array}{l} \left[\! {\begin{array}{*{20}{c}} 1 \!\! & \!\! 4 \!\! & \!\! 4\\ 4 \!\! & \!\! 6 \!\! & \!\! {12}\\ 4 \!\! & \!\! {12} \!\! & \!\! 6 \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} 4 \!\! & \!\! 6 \!\! & \!\! {12}\\ 6 \!\! & \!\! 4 \!\! & \!\! {12}\\ {12} \!\! & \!\! {12} \!\! & \!\! {12} \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} 4 \!\! & \!\! {12} \!\! & \!\! 6\\ {12} \!\! & \!\! {12} \!\! & \!\! {12}\\ 6 \!\! & \!\! {12} \!\! & \!\! 4 \end{array}} \! \right],\\ \left[\! {\begin{array}{*{20}{c}} 4 \!\! & \!\! 6 \!\! & \!\! {12}\\ 6 \!\! & \!\! 4 \!\! & \!\! {12}\\ {12} \!\! & \!\! {12} \!\! & \!\! {12} \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} 6 \!\! & \!\! 4 \!\! & \!\! {12}\\ 4 \!\! & \!\! 1 \!\! & \!\! 4\\ {12} \!\! & \!\! 4 \!\! & \!\! 6 \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} {12} \!\! & \!\! {12} \!\! & \!\! {12}\\ {12} \!\! & \!\! 4 \!\! & \!\! 6\\ {12} \!\! & \!\! 6 \!\! & \!\! 4 \end{array}} \! \right],\\ \left[\! {\begin{array}{*{20}{c}} 4 \!\! & \!\! {12} \!\! & \!\! 6\\ {12} \!\! & \!\! {12} \!\! & \!\! {12}\\ 6 \!\! & \!\! {12} \!\! & \!\! 4 \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} {12} \!\! & \!\! {12} \!\! & \!\! {12}\\ {12} \!\! & \!\! 4 \!\! & \!\! 6\\ {12} \!\! & \!\! 6 \!\! & \!\! 4 \end{array}} \! \right] \! \left[\! {\begin{array}{*{20}{c}} 6 \!\! & \!\! {12} \!\! & \!\! 4\\ {12} \!\! & \!\! 6 \!\! & \!\! 4\\ 4 \!\! & \!\! 4 \!\! & \!\! 1 \end{array}} \! \right] \end{array} \!\!\! \right) \!\! {。} $
综上，任意给定多面体外部的一个质点，利用公式即可求解该点处的引力势能，进而可以确定小行星附近的引力场。
2 非均质小行星的引力与平衡点 2.1 非均质多面体模型对周围质点的引力计算利用理论力学的相关知识可知，多面体附近质点所受的引力与其引力势能之间存在以下关系
| $ {{{F}}_{REL}} = - \frac{{\partial U}}{{\partial {{R}}}} $ | (3) |
求取上述偏导数的方法主要有两种。其一是数值方法，即利用导数的物理意义，使得R有一个微小的增量ΔR，并求解相应的ΔU，最后通过计算二者的比值来近似求解偏导数。这种方法需要多次计算多面体的引力势能，而且对于不同尺度的问题需要确定不同量级的ΔR，在实际应用中并不可取。其二是理论推导方法，该方法在公式的基础上直接对R求取偏导数而不需要求解多面体的引力势能。相比之下，第二种方法的计算量较小，且不会产生额外的数值误差，本文选用第二种方法。这里直接给出多面体引力的表达式
| $ \begin{array}{l} {{F}_{REL}} = G\sum\limits_{a \in A}^{} {{\rho _a}\det ({T_a})\left( {\frac{{\partial {{\hat U}_0}}}{{\partial {{F}_\theta }}} + \frac{{\partial {{\hat U}_1}}}{{\partial {{F}_\theta }}} + \frac{{\partial {{\hat U}_2}}}{{\partial {{F}_\theta }}} + } \right.} \\ \quad\quad\quad \left. {\frac{{\partial {{\hat U}_3}}}{{\partial {{F}_\theta }}} + \frac{{\partial {{\hat U}_4}}}{{\partial s{{F}_\theta }}} \cdots } \right) \end{array} $ | (4) |
其中
$ \begin{aligned} \frac{{\partial {{\hat U}_0}}}{{\partial {{{R}}_\theta }}} = & - \frac{{{{Q}}{{{R}}_\theta }}}{{{{{R}}^3}}}\\ \frac{{\partial {{\hat U}_1}}}{{\partial {{{R}}_\theta }}} = & \frac{{3{{{Q}}_i}{{{R}}_\theta }{{{w}}^i}}}{{{{{R}}^5}}} - \frac{{{{{Q}}_i}{{v}}_\theta ^i}}{{{{{R}}^3}}}\\ \frac{{\partial {{\hat U}_2}}}{{\partial {{{R}}_\theta }}} = & \frac{{3{{{Q}}_{ij}}{{{R}}^{ij}}{{{R}}_\theta }}}{{2{{{R}}^5}}} - \frac{{15{{{Q}}_{ij}}R_\theta ^{}{{{w}}^i}{{{w}}^j}}}{{2{{{R}}^7}}} + \frac{{3{{{Q}}_{ij}}{{{w}}^i}{{v}}_\theta ^j}}{{{{{R}}^5}}}\\ \frac{{\partial {{\hat U}_3}}}{{\partial {{{R}}_\theta }}} = & - \frac{{15{{{Q}}_{ijk}}{{{R}}^{ij}}{{{R}}_\theta }{{{w}}^k}}}{{2{{{R}}^7}}} + \frac{{3{{{Q}}_{ijk}}{{{R}}^{ij}}{{v}}_\theta ^k}}{{2{{{R}}^5}}} + \frac{{35{{{Q}}_{ijk}}{{{R}}_\theta }{{{w}}^i}{{{w}}^j}{{{w}}^k}}}{{2{{{R}}^9}}} - \frac{{15{{{Q}}_{ijk}}{{{w}}^i}{{{w}}^j}{{v}}_\theta ^k}}{{2{{{R}}^7}}} \\ \frac{{\partial {{\hat U}_4}}}{{\partial {{{R}}_\theta }}} = & - \frac{{15{{{Q}}_{ijkl}}{{{R}}^{ij}}{{{R}}^{kl}}{{{R}}_\theta }}}{{8{{{R}}^7}}} + \end{aligned} $
$ \begin{aligned} \frac{{105{{{Q}}_{ijkl}}{{{R}}^{ij}}{{{w}}^k}{{{w}}^l}{{{R}}_\theta }}}{{4{{{R}}^9}}} - \frac{{15{{{Q}}_{ijkl}}{{{R}}^{ij}}{{{w}}^k}{{v}}_\theta ^l}}{{2{{{R}}^7}}} - \frac{{315{{{Q}}_{ijkl}}{{{w}}^i}{{{w}}^j}{{{w}}^k}{{{w}}^l}{{{R}}_\theta }}}{{8{{{R}}^{11}}}} + \frac{{35{{{Q}}_{ijkl}}{{{w}}^i}{{{w}}^j}{{{w}}^k}{{v}}_\theta ^l}}{{2{{{R}}^9}}} \end{aligned} $
只要将多面体的形状参数、密度分布和目标点的位置矢量代入式(4)即可求得该点所受的引力。利用这种方法可以求取旋转小行星附近的平衡点位置。
2.2 旋转小行星附近的平衡点计算在小行星本体系中，任意质点(矢径为R)的有效势能为
| $ V = U - {T_0} = U - \frac{1}{2}(\Omega \times {{R}}) \cdot (\Omega \times {{R}}) $ | (5) |
其中：U可由式(4)得到。由于小行星绕其中心惯量主轴Z轴转动，因此式(5)可以简化为
| $ V = U - {T_0} = U - \frac{1}{2}{\omega ^2}\left( {{x^2} + {y^2}} \right) $ | (6) |
利用航天器动力学的相关知识可知，满足有效势对位移的一阶导数为零的点即为旋转小行星附近的平衡点。求解以下方程可以得到平衡点的位置，记为${L_i} = \left( {{x_i},{y_i},{z_i}} \right)$。
| $ \left\{ {\begin{array}{*{20}{c}} {{U_x} + {\omega ^2}x = 0}\\ {{U_y} + {\omega ^2}y = 0}\\ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\! {{U_z} = 0} \end{array}} \right. $ | (7) |
平衡点位置的确定是小行星附近轨道动力学研究中一项十分重要的工作，利用本文方法求解非均质小行星附近点的引力势，并代入式(7)，即可求解小行星附近平衡点位置。
3 准确性验证与算例展示 3.1 方法准确性检验前文提到，文献[2]开发出一套求解小行星附近引力场精确解的方法，文献[4]应用该方法求解了一系列小行星附近的引力场和平衡点分布，这里暂且将该方法称为传统方法。本文提出了求解非均质小行星附近引力场和平衡点的方法，为了检验算法和程序的准确性，令小行星各部分的密度相同，分别使用传统方法与本文方法求解小行星附近的引力场和平衡点，通过比较二者差异检验本文方法的准确性，结果如下。
选取小行星951 Gaspra，分别使用传统方法和本文方法求取该小行星引力场和平衡点。引力场对比图如图 2所示，平衡点位置如表 1所示。
|
|
图 1 以小行星951Gaspra的形心为中心，半径12 km的球面上的引力场分布图 Fig. 1 Contour plots of the gravitational potential of 951Gaspra on a sphere of 12 km radius |
|
表 1 小行星951 Gaspra平衡点相对于形心的坐标 Table 1 Positions of the equilibrium points for asteroid 951 Gaspra |
通过上述比较可以发现，使用本文方法求得的小行星附近引力场和平衡点位置与使用传统方法求得的结果相差较小(相对误差小于2%)，从而验证了本文方法的准确性。值得说明的是，本文利用勒让德级数求解引力势时只保留了其0~4阶项，故本方法存在一定的截断误差。在本例中，L2点、L4点的x坐标最大绝对误差达到近200 m，但是考虑到这两个点都位于小行星本体系的x轴附近，y坐标占主导，x坐标相对于y坐标的误差大约在1.5%左右，仍然认为误差较小。
3.2 均质、非均质小行星引力场对比在目前的研究中，大都把小行星视为均质体以降低问题的复杂程度。但是对于宇宙中真实存在的天体而言，均质体所占比例少之又少。如果考虑到小行星可能存在的非均质特性，应用传统方法求得的引力场和平衡点势必存在误差。为了衡量这种潜在的误差，我们假设我们观测到的小行星并非均质体，并研究它们相应的引力场分布与平衡点位置。通过比较非均质模型与均质模型之间的差别可以衡量上述误差，进而说明非均质引力场模型的研究价值。
依旧选取小行星951 Gaspra作为研究对象，利用本文方法分别求出均质情况下它的引力场分布和平衡点位置。在现有的研究数据中，小行星951 Gaspra被分割成了5 040个四面体单体。在考虑非均质引力场特性时，假设其靠近x轴正方向的1 040个单体的密度由原来的2 700 kg/m3增加到3 700 kg/m3，靠近x轴负方向的1 000个单体的密度减少为到1 700 kg/m3，并计算其引力场和平衡点位置。值得说明的是，对于951 Gaspra这类两端不对称的小行星而言，上述假想质量分布是完全可能存在的。两种情况下的对比结果如图 3和表 2所示。
|
图 3 以小行星951 Gaspra的形心为中心，半径12 km的球面上的引力场分布图 Fig. 3 Contour plots of the gravitational potential of 951 Gaspra on a sphere of 12 km radius |
|
|
表 2 均质、非均质情况下小行星951 Gaspra平衡点相对于形心的坐标 Table 2 Positions of the equilibrium points for asteroid 951 Gaspra |
本例仅仅考虑了小行星在本体系x轴方向上可能存在的非均质特性，观察图 3可以发现，其附近的引力场发生了明显的变化。通过表 2进一步得出，951 Gaspra的x轴上的平衡点相对于形心的位置发生了明显的变化，特别是L3点相对于原来位置的偏移量已经达到10%。选取其他小行星重复上述过程，结果类似，这里不再给出具体算例。
综上，小行星可能存在的非均质特性会在很大程度上影响小行星附近的动力学环境，需要予以足够的重视。
4 结束语本文以双多面体模型相互势能求解方法为基础，提出了一种新的基于多面体模型的非均质小行星的引力场研究方法。从理论上给出求解非均质小行星引力场、平衡点的公式，并通过与以往成熟工作的对比检验了本方法的准确性。考虑了小行星可能存在的非均质特性，评估了该特性对小行星附近动力学环境可能造成的影响，希望借此引起读者和广大研究人员对于非均质引力场模型的重视。
影响本文方法求解精度的主要因素是勒让德级数的截断误差。对于形状极不规则或者距离小行星极近的情况，本方法存在较大误差。在后续的非均质小行星研究过程中，应该尽可能地提高本方法的截断阶数，并着力开发其他更加精确的求解算法。
| [1] | 崔平远, 乔栋. 小天体附近轨道动力学与控制研究现状与展望[J]. 力学进展, 2013, 43(5): 526-539. Cui P Y,Qiao D.Research progress and prospect of orbital dynamics and control near small bodies[J].Advances in Mechanics,2013, 43(5): 526-539.( 1) |
| [2] | 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) |
| [3] | Miller J K, Konopliv A S, Antreasian P G, et al. Determination of shape, gravity, and rotational state of asteroid 433 Eros[J]. Icarus, 2002, 155(1): 3-17.(1) |
| [4] | Scheeres D J, Broschart S, Ostro S J, et al. The dynamical environment about Asteroid 25143 Itokawa: target of the Hayabusa Mission[C]//AIAA/AAS Astrodynamics Specialist Conference and Exhibit. [S. l]: AIAA, 2004: 16-19.(2) |
| [5] | Wang X, Jiang Y, Gong S. Analysis of the potential field and equilibrium points of irregular-shaped minor celestial bodies[J]. Astrophysics and Space Science, 2014, 353(1): 105-121.(1) |
| [6] | Fahnestock E G, Scheeres D J. Simulation of the full two rigid body problem using polyhedral mutual potential and potential derivatives approach[J]. Celestial Mechanics and Dynamical Astronomy, 2006, 96(3-4): 317-339.(2) |
| [7] | 张振江, 崔祜涛, 任高峰. 不规则形状小行星引力环境建模及球谐系数求取方法[J]. 航天器环境工程, 2012 (3): 383-388. Zhang Z J,Cui H T,Ren G F.Modeling for the gravitation potential environment of anirregular-shaped asteroid and the spherical harmonic coefficient estimation[J].Spacecraft Environment Engineering,2010(3): 383-388.(1) |
| [8] | Scheeres D J, Khushalani B, Werner R A. Estimating asteroid density distributions from shape and gravity information[J]. Planetary and Space Science, 2000, 48(10): 965-971.(1) |
| [9] | Werner R A, Scheeres D J. Mutual potential of homogeneous polyhedra[J]. Celestial Mechanics and Dynamical Astronomy, 2005, 91(3-4): 337-349.(1) |