不规则形状小行星引力环境建模及球谐系数求取方法

张振江，崔祜涛，任高峰

（哈尔滨工业大学 深空探测基础研究中心，哈尔滨 150080）

0 引言

近年来，小行星和彗星探测任务越来越多地受到人们的重视。这些已经发射或计划中的航天任务促使人们对于一个崭新的天体力学领域——绕不规则形状小天体的卫星轨道动力学的探索。这个极富挑战性的问题正吸引着越来越多的科学家对其进行研究[1-2]。小行星极不规则的形状是造成其卫星动力学特性与近球形大行星或天然卫星不同的主要原因之一。故对不规则形状小行星引力势能的建模就成为研究绕小行星卫星轨道动力学特性所需要解决的首要问题。小行星引力势能的建模方法很多[3]，这些方法可以分为数值法和解析法两大类。

数值法的主要思想是用一个多面体或多个特定形状质元（如球体或立方体）的组合来逼近小行星的形状，再通过积分变换将引力势能中的三重积分转化为可计算的曲线和曲面积分。此类方法充分考虑了小行星的不规则形状，利用了地面天文观测和掠飞任务所拍摄的小行星的图像信息，因而具有很高的精度。数值法的缺点在于它没有解析的形式，只能求出给定位置的引力势能和引力加速度值，无法分析各阶摄动项的大小及其对卫星轨道的影响。因此，数值法只能在数值仿真中使用而无法应用于轨道设计。

解析法的主要思想是用级数展开式直接逼近引力势能，这类方法具有解析的表达形式，算法简单，各阶摄动大小明显。因此解析法广泛应用于人造卫星轨道动力学和天体力学，尤其在分析摄动对轨道的影响、设计某些使命轨道时，只能使用此类方法。解析法的最大缺点是级数项的系数难于确定，例如在应用最广泛的球谐函数模型中，各阶次球谐系数Clm和Slm都是由卫星轨道数据通过导航算法“事后”解算出来的[4]。而对于还没有绕飞任务的小行星来说，因为没有卫星轨道数据作为观测量，故其准确的球谐系数是无法得到的。以往学者的做法是将小行星简化成三轴椭球体，再计算其各阶次的球谐系数，以此作为小行星引力场的球谐系数进行分析和计算[5-6]。显然，将形状十分复杂的小行星简化成球体或三轴椭球体这是不精确的。

针对在任务发射前小行星引力场建模过程中出现的三轴椭球体模型精度过低、高精度球谐系数由于没有观测量又无法获取的问题，本文提出了一种基于多面体模型的不规则形状小行星引力场球谐系数的求取方法。其思路是先由多面体模型方法重构出小行星外部引力场的引力势能情况，并以此作为求解引力场球谐系数的虚拟观测量，再根据球谐系数与引力势能间的关系，采用最小二乘方法解算出各阶次球谐系数。

1 多面体模型方法

1.1 模型的构建

1996年，R. A. Werner提出用多面体来逼近小行星的形状进而求其引力势能的方法[7-8]，即多面体模型方法。在小行星引力势能建模的数值方法中，多面体模型应用最为广泛。如图1所示，多面体模型就是用一个表面由一系列三角形构成的多面体来逼近小行星形状，这样只要求出多面体的引力势能分布，就可以知道小行星引力势能的分布情况了。

图1 小行星216 Kleopatra（左）和6489 Golevka（右）的多面体模型Fig. 1 Polyhedral model of asteroid 216 Kleopatra(left) and 6489 Golevka(right)

1.2 引力势能计算方法

任意形状中心天体的引力势能可由

式(1)为一个三重积分。对于形状不规则的物体，式(1)是不可积的；而对于匀质的多面体来说，可以采用积分变换的方法，将此三重积分转换成第二型曲线积分和曲面积分。而这些积分可以写成有限项之和的形式，如式(2)所示[9]。

以上即为多面体模型方法的全部公式，只要给定小行星固连坐标系中一个点的位置矢量，就可以由以上公式计算出该点的引力势能。下面介绍以引力势能为虚拟观测量求取球谐系数的方法。

2 球谐系数的计算

由多面体模型方法计算出小行星的引力势能后，可以将引力势能作为虚拟观测量，通过引力势能与球谐系数之间的关系采用类似导航算法中的最小二乘法来求解球谐系数。

2.1 引力势能与球谐系数间的关系

引力势能的球谐函数模型为

式中： GM为小行星的引力常量；R为检验点到中心天体质心的距离；a为小行星参考球体的半径；Plm(sin ϕ)为sinϕ的缔合勒让德多项式，当m=0时，即退化为一般的勒让德多项式；Clm和Slm为球谐系数。

略去式(3)中的高阶项，可得到下述形式的引力势能计算公式：

式(4)中n和l为所取球谐系数的阶数和次数，η为截断误差。由此可知，引力势能与各阶次球谐系数之间成线性关系，其系数项为 Plmsin ϕc os mλ和Plmsinϕ s inmλ。这样，在已知引力势能分布的情况下，我们就可以通过这个关系求得球谐系数的具体值。

2.2 构造线性方程组

式(5)中，方程未知量的个数为 N( N + 2)（Sl0无意义）。理论上只要在小行星的引力场中取N( N + 2)个不同的检验点就可以构造一个N( N + 2)维的线性方程组，其形式如式(6)所示：

上式为一个线性方程组，可以写成Ax=b的形式。解此方程组即可求出球谐系数Clm和Slm。

2.3 方程组的解法

理论上只要按式(6)构造一个 N( N + 2)维的线性方程组就可以唯一地解出球谐系数Clm和Slm，但实际发现，由于方程系数中存在这样的因子，使得不同检验点得到的方程可能是线性相关的。这使得原给定方程组变成欠定方程组，没有唯一解。有时虽然各个方程是线性无关的，但矩阵A的条件数很大，即矩阵A为病态，这时方程的解会因存在很大的误差而失去意义。

由线性方程理论可知，形如Ax=b的超定方程组一般有3种解法，分别是：

1）将原方程转化为正则方程 ATA x=ATb，并由此解得 x =(ATA)-1ATb，此方法虽然计算简单，但其缺点是矩阵 (ATA)的条件数过大，因而结果精度比较低。

2）由奇异值分解理论求得广义逆矩阵 A+，并由x =A+b求得方程的解，此方法可以获得可靠的方程解，即使矩阵A发生列秩亏损时，依然可以给出最小范数的最小二乘解。缺点是速度比较慢。

3）对原方程进行Householder变换，从而直接求出方程的解。此方法可靠性稍逊色于第二种方法，但速度较快，在矩阵A发生列秩亏损时，所给出的最小二乘解具有最少的非零元素。本文采用的就是这种方法。

关于方程数目的选择问题，理论上方程数目越多，解的结果就越理想，但同时计算速度也会越慢，因此我们需要在计算速度和精度之间做一个权衡。通过大量计算发现：一般方程数目取为未知数数目的2 3倍时即可获得较高的精度。当然，解的精度除了与方程数目有关外，还与中心小天体形状的复杂程度有关。若中心小天体的形状较复杂，应适当增加方程数目以保证解的精度。

3 对比验证

本部分包含两个算例：第一个算例应用本文方法计算一个给定三轴椭球体的球谐函数，计算结果与真实值相比较，以验证算法的正确性及算法精度与多面体面数之间的关系；第二个算例采用本文方法计算Eros 433小行星的球谐系数，其结果与传统的三轴椭球体方法和由NEAR探测器轨道数据得到的结果相比较，以验证本文所采用的方法其优点和不足。

3.1 三轴椭球体验证

匀质三轴椭球体球谐系数的真实值可由如下解析公式给出[10]：

1）当l , m 为所有正整数时，Slm=0；

2）当l或 m 为奇数时，Clm=0；

3）当l , m为其他情况时：

式(7)中：a为参考球体的半径；δ0m为克罗内克符号，其具体值为：

在本算例中，分别取三轴椭球体的3个半轴长为16 km、8 km和6 km。

图2为算例中使用的三轴椭球体的多面体模型，应用本文方法分别取多面体的外表面数目N为400和20 000进行计算。所得计算结果与由公式(7)计算得到的真实球谐系数相比较，结果如表1所示。

图2 三轴椭球体多面体模型Fig. 2 Polyhedral model of tri-axial ellipsoid

表1 三轴椭球体球谐系数Table 1 Spherical harmonic coefficients of tri-axial ellipsoid

续表1

由表1结果可知：

1）本文提出的方法可以准确地求出匀质三轴椭球体的球谐系数，且结果与真实值（解析方法计算的结果）之间的误差很小。如在表1的数据中，当多面体的面数取为20 000时，最大误差减小为0.0923%。

2）本方法的误差随多面体外表面数目的增加而减小，如N=400时，最大误差为2.5%；而当N增大到20 000时，最大误差减小为0.0923%。

3.2 Eros 433算例分析

1998年12月24日，NEAR探测器以1 km/s的速度，从距离Eros 433小行星4 100 km处飞过。在此前后对其的大小、形状、有无磁场和卫星等进行了观测。整个掠飞过程中，NEAR的多色照相机拍摄了222张Eros 433的照片，可分辨其表面小至500 km范围内的细节。天文学家根据这些照片信息建立了Eros 433小行星的三维模型。此模型形状数据可以由http://www.psi.edu/pds/resource/nearmod.html网站得到。本文采用的形状模型由64 800个数据点组成，每个数据点由小行星表面一点的经度、纬度和到小行星中心的距离构成。将这些数据点从球坐标系转化至笛卡尔坐标系下，经三角剖分后即可得到小行星的多面体模型。如图3所示。

图3 Eros 433小行星多面体模型Fig. 3 Polyhedral model of asteroid Eros 433

计算时，首先采用传统方法将其简化为三轴椭球体，简化得到的3个半轴长分别为16 km、8 km和6 km，将这3个参数代入公式(7)计算得出球谐系数。再应用本文方法计算出其球谐系数，最后将二者结果与由NEAR探测器环绕轨道数据解算的球谐系数相比较，结果见表2。

表2 Eros 433引力场球谐系数Table 2 Gravity spherical harmonic coefficients of Eros 433

由表2结果可知：

1）与将小行星近似为三轴椭球体的方法相比，本文提出的方法可以大幅度提高结果的精度，例如采用三轴椭球体方法计算的C20与飞行器轨道数据解算的结果相差17.4%，而采用本文的方法可以将这一误差减小到0.23%。

2）本方法计算结果与真实值尚存在一定偏差，例如前4阶次球谐系数中相差最大的S22项与实际值相差6%，其他阶次的球谐系数也存在着一定差异。这些差异主要来自两个方面：第一，多面体模型与小行星的真实形状并不完全相同，二者之间存在着一定的误差，我们可以通过增加多面体外表面数目的方法来减小这一误差；第二，多面体模型方法假定小行星是匀质的，但真实的小行星密度并不是处处相同的，这也会造成结果的误差。而这种误差是方法中没有建模的，所以无法消除。这也是本方法的局限所在。

4 结束语

本文提出的多面体模型方法可以在任务发射前得到比较精确的不规则形状小行星引力场的球谐系数，因为其充分利用了天文观测或掠飞任务所拍摄的图像信息，故其精度要比三轴椭球体模型方法有大幅度提高。本文方法计算结果的精度随多面体模型外表面的数目增加而提高。对于多面体模型与小行星形状不完全一致造成的误差，可通过增加多面体的面数减小误差；而对于小行星质量不均匀所引起的结果误差，因多面体模型中无法对密度分布建模，故无法消除或减小这种误差。

（

）

[1] Scheeres D J. The dynamical evolution of uniformly rotating asteroids subject to YORP[J]. Icarus, 2007, 188: 430-450

[2] Byram S M, Scheeres D J, Combi M R. Models for the comet dynamical environment[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(5)

[3] Casotto S, Musotto S. Methods for computing the potential of an irregular, homogeneous, solid body and its gradient, AIAA 2000-4023[R]

[4] 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: 3-17

[5] Scheeres D J. Dynamics about uniformly rotating tri-axial ellipsoids, applications to asteroids[J]. Icarus, 1994, 110: 225-238

[6] Washabaugh P D, Scheeres D J. Energy and stress distributions in ellipsoids[J]. Icarus, 2002 , 159(2): 314-321

[7] Werner R. On the gravity field of irregularly shaped celestial bodies[D]. The University of Texas at Austin, 1996

[8] 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, 1997, 65: 313-344

[9] Park R S, Werner R A, Bhaskaran S. Estimating small-body gravity field from shape model and navigation data, AIAA/AAS Astrodynamics Specialist Conference and Exhibit, 2008-08-18[R]

[10] Boyce W. Comment on a formula for the gravitational harmonic coefficients of a triaxial ellipsoid[J]. Celestial Mechanics and Dynamical Astronomy, 1997, 67: 107-110