小天体引力场建模技术进展

尚海滨，韦炳威，卢榉承

（北京理工大学 宇航学院，北京 100081）

引言

太阳系小天体通常是指太阳系中除行星、矮行星以及天然卫星以外的天体，它们一般包括小行星、彗星和流星体等。近年来，太阳系小天体已经成为人类深空探测活动的重要目标。对小天体实施探测主要有3方面的动机：探测太阳系小天体对揭开太阳系形成与演化之谜有重要意义[1-3]；小天体探测有助于应对近地小天体对地球的碰撞威胁，保障人类社会可持续发展[4-5]；通过光谱观测发现一些小天体蕴含着丰富资源[6-7]，对这些资源的开发与利用具有潜在经济价值，同时可为太空工业的发展提供重要支持。

早期小天体的探测方式主要以飞越、环绕为主，美国国家航空航天局（National Aeronautics and Space Administration，NASA）发射的“伽利略号”（Galileo）探测器依次飞越了小行星Gaspra 951 和Ida 243 及其卫星Dactyl[8]；“深空一号”（Deep Space 1）探测器飞越了小行星布雷尔（Braille 9969）和彗星布洛利（Borrelly 19P ）[9]；“黎明号”（Dawn）探测器对灶神星（Vesta）进行了长达1 a的环绕探测[10]；“嫦娥二号”（Chang’E-2，CE-2）探测器在其拓展任务中成功飞越了小行星图塔蒂斯（Toutatis 4179）[11-13]。这些任务获取了大量的小天体探测数据，加深了人类对太阳系小天体的认识。

为获得更为丰富的科学数据，例如岩石特征、矿物成分、地表热特性和磁特性等，需对小天体进行近距离和表面的深入探测。为此，各航天机构相继开展了附着、表面巡视及采样返回等复杂探测任务。日本宇宙航空研究开发机构（Japan Aerospace eXploration Agency，JAXA）开展了“隼鸟一号”（Hayabusa 1）任务[14]和“隼鸟二号”（Hayabusa 2）任务[15]、欧洲航天局（European Space Agency，ESA）开展了“罗塞塔”（Rosetta）任务[16]、NASA开展了“舒梅克探测器”（NEAR-Shoemaker）任务[17]和“欧西里斯号”（OSIRIS-REx）任务[18]。对于近距离探测任务而言，小天体引力场的作用变得不可忽视，特别对于无轨控能力的巡视器和着陆器，复杂引力场导致了更高的任务风险。Hayabusa探测器于2005年11月12日释放了着陆器“智慧女神”（Minerva），试图着陆于小行星糸川（Itokawa 25143 ）表面，由于释放位置高度过高导致Minerva最终逃逸了Itokawa[19]。近年来，各个国家也在积极开展双小行星系统探测的国际合作，由NASA和ESA合作开展的小行星防御与偏转评估任务（Asteroid Impact &Deflection Assessment，AIDA）已经完成了对目标双小行星系统Didymos 65803的初步探测[20]，并计划由ESA完成对该双小行星系统的近距离详细勘测。两个小行星复杂的引力场和双小行星系统独特的构型，为双小行星系统附近的引力场计算带来了新的挑战。

精确构建小天体附近及表面引力场模型是发展近距离探测技术，降低任务风险的重要基础。若能实现小天体引力场的高精度建模，将有助于实现探测器、巡视器和着陆器的任务可靠性设计，包括环绕轨道、悬停轨道、转移轨道、受控附着和弹道附着轨道的设计和稳定性分析，以及表面巡视轨迹的鲁棒规划等。除了实际任务需求，小天体精确引力场建模还是一些动力学和行星科学研究的基础，例如周期轨道的搜索、表面溅射物的演化、双小行星系统的演化和小天体内部物质分布重构等。

由于小天体具有不规则形状和非均匀内部质量分布的重要特点，其引力场的引力势、引力加速度和引力梯度通常不具备解析闭合的表达形式。因此，从数学上描述小天体附近的复杂引力场具有很大的挑战性，多年来许多学者们致力于解决这一难题，获得了诸多重要成果。

本文梳理了两百多年来学者们对引力场建模问题的认识与发展，分别详细阐述了球谐函数展开法、椭球谐函数展开法、内外引力场法、球Bessel函数展开法、质点群法、多面体法、虚拟边界元法以及基于数据驱动的快速引力场计算方法的基本原理，评述了不同方法的优点与不足。除引力场表征问题，在此基础上还回顾了双小行星动力学系统建模技术的发展，按照从简单到复杂建模历程分别介绍了球形限制假设模型、级数展开模型和有限元模型，可为双小行星系统附近的引力场建模技术研究提供参考。

1 单小天体引力场建模技术

对于小天体近距离探测任务而言，探测器在探测目标附近的轨道运动主要受到小天体引力作用的影响，精确构建小天体附近的引力场有助于开展小天体探测任务设计。同时，小天体的引力场建模还是天体动力学分析的基础，能够支持行星科学的研究，具有较高的科学价值。

位于小天体外部空间一点r的探测器所感受到的小天体引力势定义为

其中：Ω为小天体体积区域；s为Ω内一点；G为引力常数；ρ 为小天体内部的密度分布函数。

探测器感受到的引力加速度则可表示为∇U，引力梯度张量可表示为∇∇U，具体表达式为

1.1 级数展开法

1.1.1 球谐函数展开法

无穷级数展开法是利用引力势函数的解析性质，将原求解不规则外形非均质分布小天体的引力势问题转化为计算有限项基函数的线性组合问题。球谐函数展开法是最为经典的引力场表征方法。球谐函数构成一组完全正交函数基，在球面上定义的连续函数都可写成球谐函数的线性组合。类似于定义在圆周上的周期函数，可通过Fourier级数表示为三角函数的线性组合。最初法国数学家Legendre于1782年在研究质点系引力场时发现观察点到小天体内微元距离的倒数能以级数形式表示为

其中：r=‖r‖；s=‖s‖；Pn(·)为n阶Legendre多项式；θ为矢量r与矢量s之间的夹角。

William 等[21]将这一工作进一步完善，在球坐标下任意天体的引力势可表示为

其中：M为天体质量；R为参考半径，通常选取天体最小包围球（称为Brillouin球）的半径；Cnm和Snm称为Stokes系数或球谐系数。

级数式（5）中缔合Legendre多项式与三角函数的乘积体现了球面函数的部分，称为球谐函数。级数式（4）形式仅当s

图1 Brillouin球与非球形小天体示意图Fig.1 Illustration of Brillouin sphere and non-spherical small body

1.1.2 椭球谐函数展开法

球谐函数展开法对于非球形天体引力场描述有局限性，为扩大收敛区域，Hobson等[22]总结出了椭球谐函数展开法。该方法首先引入椭球坐标系(λ1,λ2,λ3)描述空间一点。这3个坐标参数分别描述了3个共焦二次曲面：椭球曲面、单叶双曲面、双叶双曲面。3个共焦曲面在每个笛卡尔坐标象限中决定一个空间点。

已知引力势函数满足Laplace方程，即∇2U=0，则在椭球坐标下描述关于U的Laplace方程，再根据分离变量假设方程的解可写为U=E1(λ1)E2(λ2)E3(λ3)，代入椭球坐标系下的Laplace方程可得到对应的Lamé方程。Lamé方程的解称为Lamé函数或椭球谐函数，最终引力势函数可以表示为椭球谐函数的级数形式，即

其中：αnm为待定的椭球谐系数为参考椭球对应的坐标值；Enm(·)为第一类Lamé函数；Fnm(·)为第二类Lamé函数。

级数式（6）的收敛域为参考椭球外部，即λ1≥对于细长型小天体而言，可将参考椭球定义为包围小天体的最小椭球（即Brillouin椭球），椭球谐函数展开法能够描述更多靠近小天体表面的空间区域（如图 2所示）。

图2 Brillouin椭球与非球形小天体示意Fig.2 Illustration of Brillouin ellipsoid and non-spherical small body

Walter[23]就引力场建模中椭球谐系数和球谐系数相互转化的问题开展了研究，但其转换过程涉及大量椭圆积分的复杂数值计算，使得所提方法的应用受到了限制。Dechambre等[24]基于Dirichlet问题发展了另一套球谐系数到椭球谐系数的转换方法，简化了椭球谐系数的计算，之后学者致力于发展精确系数转换方法[25-27]。尽管椭球谐函数展开法具有更好的收敛特性，但其应用远不如球谐函数展开法，这是因为计算各项椭球谐函数相比各阶球谐函数具有更高的计算复杂度，甚至Lamé函数的递归公式还没被发现，而缔合Legendre多项式具有成熟的递归公式体系，能够便于计算机的求解[28]。由于第一类和第二类Lamé函数存在零点和奇异点（当λi=h或λi=k），椭球谐级数在阶数较大时候容易出现数值不稳定的情况，导致计算精度退化[28]。

1.1.3 内外引力场法

为将级数展开法的收敛域扩大到小天体表面任意一点，Takahashi等[29]提出了内外球谐函数展开法。基本思想是以Brillouin球外部空间任意一点建立坐标系，并将式（4）改写为

在式（7）中，当r

其中：Rin为内Brillouin球半径；为内球谐系数，如图3所示，通过内球谐函数展开式可以将级数收敛域扩展到小天体表面一点。

图3 内外Brillouin球与非球形小天体示意Fig.3 Illustration of internal and external Brillouin sphere and non-spherical small body

虽然内外引力场法能通过内球谐函数展开式来描述原Brillouin球内的部分区域引力场，但该方法仍无法使用统一表达式实现引力场的全局表征。内Brillouin球仅能延伸至小天体表面一点，因此该方法可应用于探测器在小天体表面一点的附着模拟。若要描述小天体表面其他区域引力场，则还需建立更多内Brillouin球来覆盖更多小天体表面附近区域。

1.1.4 球Bessel函数展开法

为得到Brillouin球内部统一的引力场级数表征形式，Herrera-Sucarrat等[30-31]提出了球Bessel函数展开法。基本思想是将Brillouin球内视为一个具有一定质量分布的实物体区域，则引力势函数在该区域内满足Poisson方程，即

为求解式（9）方程，可假设Brillouin球内的质量被重新分布以满足下述Helmholtz方程。

其中：k为常数。

采用分离变量假设引力势函数具有的形式为：U(r,φ,λ)=(r)Φ(φ)Λ(λ)，其中径向势函数的波浪字符是为了与Brillouin球半径R做区分。将该形式代入Helmholtz方程并表示为Bessel方程（对应不同本征值k=km的本征方程）。

式（11）的解为第一类球Bessel函数jn(klnr)，于是Poisson方程（9）的解可表示为

其中：αln=klnR为归一化的本征值。

若将Brillouin球外部的引力势函数也写成球Bessel函数和球谐函数的形式可得

球Bessel函数展开法的难点是系数的求解。由于系数Alnm和Blnm的求解精度限制于原球谐系数精度以及人为引入的约束，通常只能获得低阶系数。

1.2 质点群法

质点群法的主要原理是利用有限个质点或均质小球填充小天体内部空间，则小天体外部的引力势可以表示为所有质点产生的引力势的叠加，即

该方法避免了复杂的体积积分计算，计算效率较高且建模方式简单易实现，广泛应用于小天体演化和附近轨道动力学的研究中，Richardson等[32-33]在N体问题开源求解器pkdgrav的框架下对小行星Itokawa的形成机理进行了研究；Geissler等[34]使用质点群模型研究了小行星Ida及其卫星Dactyl表面溅射物的再沉积过程；Shang等[35]以质点群法作为初始引力场模型，提出了小行星附近全局周期轨道的快速搜索方法，并研究了小行星Eros附近周期轨道的拓扑类型。

质点群法相比无穷级数展开法的主要优势：①能够较好拟合小天体的不规则形状；②通过给每个质点设置不同的质量可以表征小天体内部质量非均匀分布的特点；③不存在收敛域的限制，理论上属于小天体引力场的全局表征方法。然而，质点群法的缺点也较为明显，主要表现：①模型精度随质点个数增加提升速度较慢，导致当精度需求较高时，所需质点群总数将剧烈增加，计算效率显著降低；②虽然质点群能够近似表征小天体内部质量非均匀分布特征，但当质点数较大时会引起协方差矩阵奇异，导致各个质点的质量难以准确估计，Park等[36]证明了当质点数增大，通过轨道数据估计的各质点质量的不确定度迅速增加；③当观察点靠近小天体表面，引力场的精度会迅速退化，这是由观察点靠近质点群时所引起的近奇异现象所致，因此该方法并不适用于小天体表面的动力学研究。

1.3 多面体法

1.3.1 均质多面体法

Werner等[37]在前人研究的基础上首先完整推导了任意均质多面体引力势、引力加速度和引力梯度张量的闭合解，被誉为小天体引力场的优美解。该方法假设小天体内部的密度函数为常值，可寻找某个矢量场函数使得其关于小天体内质量微元坐标的散度恰好等于Newton核，利用Gauss散度定理将原计算引力势的体积分转化为关于小天体表面的面积分，式（1）可重写为

假设小天体形状可以由多面体近似表征，计算引力势的面积分可表示为有限个三角形面的积分之和。

其中：Ff表示多面体的第f三角形面积区域；rf表示空间点r到第f三角形面上任意一点的矢量；表示第f三角形面的单位外法线矢量。

通过Green第一恒等式的进一步表示以及面积分和线积分的转化，最终得到任意均值多面体的引力势、引力加速度和引力梯度张量

多面体法将小天体的引力场计算转化成了空间点r到多面体每个棱边和三角形面之间的几何关系计算。当使用较高分辨率的小天体多面体模型时，引力场的计算代价也通常较大。当空间点r位于多面体棱边上时多面体法出现奇异而无法计算。实际上棱边上的奇异点都是可去奇异点，多面体引力势在棱边上仍然是解析的，Tsoulis等[38]针对这一问题对引力势解析表达式进行了改进，消除了表达式的奇异部分，从而完善了多面体法。

1.3.2 非均质有限元多面体法

尽管多面体法是引力场的闭合解析形式，但在利用Gauss散度定理时引入了密度均匀假设，这导致多面体法无法反映小天体内部质量分布特征。实际上在小天体的形成过程中，物质的聚集往往是杂乱随机的，因此小天体内部的质量分布通常是非均匀的。为表征小天体内部质量分布情况，学者们采用有限元的思想，将原多面体分割成有限个更小的多面体，每个多面体赋予不同的密度值，最后通过叠加所有小多面体的解析引力场来得到真实引力场的近似值。Takahashi等[39]将小天体的多面体模型切割成有限个楔形体，通过球谐系数或者在轨观测数据反演每个楔形体的密度值。

一般来说，用于近似小天体外形的多面体具有较多的面数和边数，利用多面体法求解引力场时需要遍历计算所有面和边，这意味着该方法需要较大的计算成本。若采用有限元结合多面体法来求解非均质分布小天体的引力场，原二维流形问题变为三维问题，需要遍历的面和边的个数将剧烈增加，随之增加的计算成本通常难以承受。确定小天体内部每个有限单元的密度实际上是经典的引力场反问题，该问题并不存在唯一解，因此很难进行求解。为说明同一个引力场可对应多个质量分布，Takahashi等[39]描述了一种假想情况。学者们只能预先猜测小天体内部可能的质量分布规律，再进一步进行参数反演，但由于引力场反问题的不适定本质使得微小的测量误差能够导致很大的参数不确定性[36]。

1.4 虚拟边界元法

Wei等[40]发展了小天体引力场的虚拟边界元全局表征方法，有效解决了上述方法中不能同时满足：①表征的全局性；②表征的精确性；③表征参数的易求解性的问题。该方法从求解Laplace方程的边值问题的角度对不规则非均质小天体引力场进行建模，即小天体引力势函数是式（20）边值问题的解：

其中：Ω′表示小天体外部空间区域；n为小天体表面单位外法线矢量。

根据偏微分方程的基本理论，上述边值问题的解可表示为单层势。

其中：σ(x)为定义在边界上的待求面密度函数；G(x,r)为Laplace方程的基本解，或称为Green函数。

面密度函数σ(x)可以通过下述适定的第一类Fredholm积分方程求解。

为能够快速获得面积分(21)的具体数值，Wei和Shang[40]推导了单层势的半解析形式。首先假设小天体边界能够由多面体P 近似。

其中：N为多面体三角形面的总数；Γi是第i三角形面。

对于任意一三角形面 Γi，存在一个光滑映射 φi使得φi(T)=Γi，T ⊂R2是标准三角形区域：

根据上述多面体的定义，式（21）可被重写为

一种简单的光滑映射φi可以选为仿射变换。

其中：ξ3=1-ξ1-ξ2，Vi,1、Vi,2和Vi,3分别为第i三角形面上的三个顶点坐标矢量。

利用三角形单元上的线性插值近似估计面密度函数，即

其中：hi,1、hi,2和hi,3面密度函数分别在第i三角形面上三个顶点处的值。

故单层势可进一步表示为

函数ψi,k(ξ2,r)具有闭合的解析形式，这里不赘述其结果。根据ψi,k(ξ2,r)及其关于ξ2的导数在0和1处的性质，可以采用特殊的多项式近似逼近。

其中：Np为多项式的最高次数；系数an和bn可根据ψi,k(ξ2,r)的闭合解析形式构建线性方程组求解；多项式P(ξ2)用于近似ψi,1(ξ2,r)和ψi,3(ξ2,r)；多项式Q(ξ2)用于近似ψi,2(ξ2,r)。从而式（30）中的函数Φi,k(r)可写成如下半解析形式：

当空间点r靠近小天体边界时，由于Green函数G(x,r)的奇异性，式（34）和（35）的半解析形式精度会迅速退化。为了解决这一问题，Wei等[40]采用了虚拟边界法，其基本思想是利用函数的解析延拓概念，将原边值问题的解的定义域扩充到原定义域之外，在定义域外建立一个虚拟边界，从而空间观察点永远不会靠近边界，奇异性被根本消除。图 4给出了虚拟边界概念示意图，在小天体内部建立一个虚拟边界，使得在该边界定义的Laplace方程边值问题的解与原小天体的引力势函数相等，从而半解析形式（34）和（35）的全局精度能够被保证。

图4 虚拟边界概念示意图Fig.4 Illustration of concept of virtual boundary

这一方法被应用于电磁弹性材料问题和热传导问题的求解，然而通过虚拟边界构建的边界积分方程的适定性没有被严格证明，这限制了该方法在许多偏微分方程问题求解的应用。Wei等[40]首先利用调和函数的基本性质证明了引力势函数于小天体内部解析延拓的存在性和唯一性，利用Lax-Milgram定理证明了虚拟边界积分方程的适定性，即

虚拟边界元法与多面体法类似，需要遍历所有虚拟边界多面体的顶点和面的计算，当多面体分辨率较高时会带来较大的计算代价。

1.5 基于数据驱动的小天体引力场建模方法

针对模型计算精度和计算效率不能兼顾的问题，学者们提出了基于数据驱动的小天体引力场建模方法，旨在进一步提高引力场的计算效率。这些方法包括两大类：插值法和机器学习法。

1.5.1 插值法

空间局部插值方法试图通过扩大内存消耗来 提高引力场的计算效率，是将空间有限离散点处的引 力场信息（如引力势、引力加速度和引力梯度张量）事先存储在内存中，再通过插值计算得到空间其他任 意点的引力场信息。这一思想最早由Junkins[41]于 1976年在研究地球引力场时提出，后来由Engels等[42]一 同主张通过采用加权正交多项式作为插值基函数来提升插值精度。随着计算机硬件的发展，内存得到了迅速提升，基于空间局部插值方法被广泛应用于地球引力场的表征[43-45]。

Colombi等[46]提出了自适应八叉树空间网格结构，从而将这一思想应用到了小天体引力场的表征问题。自适应八叉树空间网格划分策略能够很好地抑制空间插值误差，使得全空间的引力场插值精度得到控制。Wei等[47]发展了结合复合引力逼近与快速网格跟踪的高精度高效率引力场建模方法。他们基于复合引力思想将引力场拆解为由椭球表征的主项与局部同胚网格插值表征的扰动项，针对扰动项，采用了15节点的楔形体单元进行插值近似。

为使得扰动项对全空间的“贡献”尽可能小，提出了引力最优椭球概念，证明了引力最优椭球求解问题等价于椭球与小天体多面体的布尔体积最小化问题，并推导了布尔体积的计算方法。基于星形多面体的特征发展了相容性空间网格划分方法，进一步，为了实现在轨道递推过程中当前点的快速准确网格定位，提出相容性网格数据结构和基于碰撞检测的网格搜索方法。

1.5.2 机器学习法

近些年来，随着计算机科学的发展与人工智能领域的兴起，基于数据驱动的机器学习方法在小天体的引力场计算中得到了广泛的应用。这类方法通过对小天体附近空间场点的数据采样，利用机器学习算法完成模型的训练，实现对小天体附近引力场的回归计算，克服了传统方法中计算精度和计算效率难以兼顾的问题。

Gao等[48]基于高斯过程回归，建立了场点位置坐标与引力场的直接映射关系，避免了复杂的建模计算过程。高斯过程的一个特点就是，对于每个变量x都有一个对应的高斯分布，而对于一组变量X={x1,x2,···,xn}，存在一个联合高斯分布，满足

对于一个新的样本xn+1，则新的高斯分布为

其中：

f是回归的目标函数；m是均值函数；k是高斯核函数。

在小天体的引力场建模计算中，回归的目标函数是小天体的引力加速度函数，变量x是小天体附近的场点位置。

人工神经网络得益于其灵活的网络结构，如图 5所示，在函数的拟合方面具有良好的性能，因此往往被应用于拟合小天体的引力场函数。Yu等[49]基于深度神经网络，拟合了小天体的引力加速度函数，同时，推导了深度神经网络引力加速度的偏导数，建立了引力梯度矩阵的近似形式，并应用于求解燃料最优的小天体着陆问题。Cheng等[50]同样地采用深度神经网络，实现了小天体引力场的快速计算，并解决了时间最优的小天体着陆问题。Furfaro等[51]采用单层反馈神经网络，建立了场点位置与小天体引力加速度的映射关系，同时，得益于单层反馈神经网络简洁的结构，避免了神经网络训练的迭代优化，有效地降低了训练神经网络的时间成本。

图5 人工神经网络示意图Fig.5 Illustration of artificial neural network

2 双小行星系统引力场建模技术

除了单小天体，太阳系还普遍存在大量双小行星系统，针对它们的探测与研究也同样对太阳系起源、行星形成与演化甚至行星防御具有重要意义。双小行星系统的引力场建模技术能为双星系统探测任务设计与稳定性分析提供动力学仿真环境，其中包含的双星系统动力学建模技术研究也是分析该系统形成与演化机制的基础。如图 6所示，位于双小行星系统附近空间一点r的探测器所感受到的小天体引力势定义为

图6 双小行星系统示意图Fig.6 Illustration of binary asteroid system

其中：si是小天体Bi内部一点相对于双星系统质心O的位置矢量；ρi是两个小天体内部的密度分布函数；r是探测器相对于双星系统质心O的位置矢量。

基于单个小天体引力场建模技术，探测器在双星系统附近的引力势是在两个小天体引力场中引力势的叠加

其中：Ui是两个小天体各自的引力势；ri是两个小天体各自的质心Oi相对于系统质心O的位置矢量；Ai是由小天体Bi固连坐标系转至惯性坐标系的姿态矩阵。

双小行星系统附近的引力场与两个小天体的指向和位置密切相关，研究两个小天体的姿态运动和轨道运动是双小行星系统引力场建模技术的关键。小天体的姿轨运动是由系统的动力学决定的，考虑两个小天体不规则的形状和非均匀的质量分布，这个系统的动力学往往呈现出姿轨耦合的特点，即全二体问题。由双小行星系统的运动学可知，双小行星系统的引力相互势是研究全二体问题的基础。然而，在一般情况下，两个不规则形状且非均匀质量分布的小天体之间的相互势没有闭合形式解。过去的几十年中，不同学者们都致力于发展能够计算双星系统引力相互势的方法，本章主要回顾了这些技术的发展历程。

2.1 双小行星系统的运动学

因为建立了相互势与两个小天体之间的引力、引力矩的偏导数关系，Maciejewski[52]的双刚体相对运动方程在双小行星系统建模中得到了广泛的应用。双小行星系统的相对运动示意图如图 7所示。双刚体相对运动方程利用两个小天体的相对轨道和姿态变量，减少了系统状态变量的数量，简化了问题的计算，在其中一个小天体B2固连坐标系中的形式为

图7 双小行星系统的相对运动Fig.7 Relative motion in binary asteroid system

其中：P是系统的线动量；R是小天体B1相对于小天体B2的位置矢量；Γ1和Γ2分别是小天体B1和B2的自转角动量；A是由小天体B2固连坐标系转至小天体B1固连坐标系的姿态矩阵；A2是由小天体B2固连坐标系转至惯性坐标系的姿态矩阵；Ω1和Ω2分别是小天体B1和B2的自转角速度。

引力相互势V的基本形式是两个天体的质量微元之间万有引力势在两个天体上的积分，如图8所示。

图8 引力相互势示意图Fig.8 Illustration of gravitational mutual potential

其中：G是万有引力常数，D1和D2分别是两个小天体内质量微元在各自小天体固连系下的位置矢量。

两个小天体之间的引力和引力矩可以通过相互势函数的偏导数表示，故建立准确的相互势模型是研究双星系统的基础。

2.2 球形限制假设模型

在双星系统引力和小天体结构强度的共同作用下，一部分双小行星系统的主星接近球形，次星接近椭球形，且与主星潮汐锁定。在早期的双小行星系统引力相互势建模中，通常将其中一颗小天体假设为匀质球体，其万有引力势等价于一个质点的引力势，这一假设往往称为球形限制假设（Sphere Restriction Assumption），假设小天体1是均质球体，相互势表示为

其中：M1是球形天体的质量。此时，双星系统的相互势就退化为一个质点在另一个不规则小天体引力场中的引力势。

若假设不规则小天体具有简单的几何形状，如三轴椭球，则该小天体的引力势函数具有简单的形式。进而，能够建立双星系统的相互势与系统模型参数（如形状尺寸参数、相对距离参数、相对指向参数、质量参数等）之间的定量关系。利用这种定量关系，分析了双星系统的平衡构型[53-54]、能量[55]，和接触双星系统的分裂[56]等。

2.3 级数展开模型

随着双星系统研究的发展，球形限制假设下对小天体形状的简化难以满足双星系统建模精度的要求。通过级数展开将相互势展开为无穷级数的形式，相互势可表示为双星系统相对位置参数、相对姿态参数、天体质量分布参数的函数。进而，通过相互势的偏导数计算双星系统引力和引力矩，基于级数展开的模型在双星系统的递推和演化分析中发挥了重要作用。

在全二体问题背景下，Borderies[57]将其中一个天体的引力势用球谐级数表示，将球谐级数在另一个天体的球坐标系中表示后，得到了基于两个天体的球谐系数的相互势表达形式。Compère 等[58]基于STF（symmetric trace free）张量，推导了一种与球谐方法等价的方法，并应用于双小行星系统的递推。Tricarico[59]基于幂级数推导了任意形状和质量分布的两个小天体之间的相互势，并证明了球谐系数与惯性积分作为天体质量分布参数的一致性。利用有限阶球谐系数和惯性积分的相互势、引力和引力矩的建模工作还可以参考文献[60～63]。

针对形状复杂的双小行星系统，学者们发展了一系列高阶级数展开方法，以提高建模精度。Werner 等[64]采用均质多面体描述两个小天体，通过叠加构成两个小天体任意四面体对之间的引力势解决总相互势的评估。利用勒让德级数来展开两个四面体之间微元的距离项，实现四面体对之间引力势的解析可积，得到基于各阶惯性张量的引力势

其中：a和b分别代表两个小天体中以多面体三角面为底面，质心为顶点的四面体；ρa和ρb分别代表两个小天体的密度；Ta和Tb分别代表四面体与标准四面体转换的雅克比行列式；Q是常系数张量；w和r是与两个小天体的位置、姿态相关的张量。

Fahnestock等[65]进一步地将Werner的工作推广至相互引力和引力矩的计算，并应用于1999KW4双小行星系统的动力学分析[66]。该方法对小天体形状描述的准确程度取决于多面体模型的分辨率，当多面体模型的分辨率很高时其计算代价是较大的。

为提高相互势评估的计算效率，Hou等[67]以两个天体边界为积分区间，利用勒让德级数对构成两个天体的任意微元间的距离进行逼近，从而将共同势表示为两个天体广义惯性积的函数

其中：t、a和b分别是满足递归关系的常系数；分别是小天体B2固连坐标系下，小天体B1和B2的广义惯性积分。由于两个小天体存在相对转动，是时变的，该方法通过坐标变换建立了两个小天体固连坐标系下广义惯性积分之间的关系，从而减少了重复计算，提高了计算效率。广义惯性积分的引入不仅避免了大量系统状态相关参数的重复计算，大幅度降低了共同势的计算代价，还方便地描述了小天体不规则的形状和不均匀的质量分布。此外，该模型具备成熟的递归公式体系，能够便于计算机的求解。

Shi等[68]将主星的引力场用多面体模型描述，通过在从星质心处利用泰勒级数对主星引力场进行展开，将共同势表示为从星惯性积和主星引力势及其导数的函数，给出了相互势展开至二阶项的形式

其中：U是主星的多面体引力势函数；∇∇U是主星的引力势能梯度阵；⊗ 表示张量计算；P是与从星的二阶惯性积分和两个小天体相对姿态相关的并矢

该模型能够对其中一颗星进行不规则引力场的精确建模，避免了同时舍去两颗星的高阶惯性积分。

基于级数展开的模型能够计算任意形状和质量分布的小天体之间的相互势，通过计算引力势偏导数的方法得到引力和引力矩。该模型的误差来源于级数的截断误差，随着两个小天体之间的距离减小，级数的截断误差增大，且此类模型的计算量随着级数阶数的提高显著增大。

2.4 有限元模型

针对基于级数展开的模型在两个小天体距离接近时，级数收敛变慢甚至发散。于洋等[69]提出了一种有限元方法来改善相互势、引力和力矩评估的收敛性和精度。这种方法并不依赖于相互势的偏导数计算引力和引力矩，而是根据直接的物理含义对两个小天体之间的引力相互作用进行建模。该方法通过对两个小天体进行四面体有限元划分，利用双线性插值对构成两个小天体任意两个微元之间的相互势、力和力矩进行计算。两个天体之间总的相互势、力和力矩可以通过微元间引力作用的简单叠加获得。

通过有限元的细致划分，基于有限元方法的模型对于小天体复杂的几何构型和局部的密度异常都是有效的。有限元模型的计算精度取决于有限元划分的精细程度，当模型的有限元划分精度较高时，该模型计算精度较高，但是大量的有限元单元也导致该方法的计算负担较大。高云峰等[70]在于洋工作的基础上利用CUDA并行计算框架进一步提升了这一方法的计算效率，在相同的精度下，与单核CPU相比，计算效率提升了两个量级。

3 未来发展趋势

小天体不规则形状和其内部不均匀的质量分布给小天体引力场计算带来了极大挑战。经历了两个世纪的漫长发展，已经就小天体的引力场建模问题有了较为全面的认识，但仍有一些问题需要进一步研究和解决。

3.1 引力场的解析建模

在天体引力场研究中，最先发展的是球谐函数展开法。该方法利用了Legendre多项式和三角函数构成的无穷级数，在球坐标下描述了任意天体的引力势函数、引力加速度和引力梯度张量等。由于Legendre展开式收敛条件的限制，球谐函数级数最终只能在天体的最小包围球（即Brillouin球）外部保证收敛，而球内部收敛无法保证。小天体普遍呈现不规则的形状，使得Brillouin球无法贴合小天体表面，导致小天体表面附近的空间引力场无法被描述。

学者们意识到球谐函数展开法的局限性后，利用椭球坐标系的描述将Brillouin球变换成Brillouin椭球，从而能够更好地贴合非球形小天体。对应的球谐函数被变换成了椭球谐函数，它由第一类和第二类Lamé函数构成。相比球谐函数，椭球谐函数具有更高的计算复杂度，且还未发现Lamé函数的递归公式，这限制了椭球谐函数的应用。另一方面，虽然Brillouin椭球能够更好地贴合非球形小天体的形状，但小天体表面绝大部分空间的引力场依然无法描述，因此椭球谐函数仍未解决无穷级数的非全局表征问题。

在过去的几十年中，学者们开始重点关注小天体表面的引力场全局表征问题，如无穷级数法中的内外引力场法和球Bessel函数展开法、质点群法、多面体法和有限元多面体法。这些方法均有侧重，例如内外引力场法侧重于将级数收敛域扩展到小天体表面一点；球Bessel函数展开法侧重于Brillouin球内部的低阶级数表达；质点群法的精度随观察点靠近小天体表面迅速退化，因此一般也用于较远距离的轨道动力学分析；多面体法虽然是全局的解析闭合引力场表征方法，但该方法不能表征非均质小天体引力场，一般仅用于理想均质多面体附近的动力学分析；有限元多面体法虽然能够近似表征非均质小天体的引力场，但每个有限单元的密度受到引力场反问题不适定本质的限制难以准确求解。总的来说这些方法均不能同时满足4方面要求：①表征的全局性；②能够表征非均质分布；③表征参数的易求解性；④表征的高效性。虚拟边界元法是有效解决前3个问题的新方法，但该方法依然存在有待解决的问题。对于任意给定小天体，虚拟边界的构建不唯一，不同虚拟边界最终影响着引力场表征精度，该方法尚未给出虚拟边界定义的标准；引力场计算精度只能通过增加虚拟边界的分辨率提高，由于该方法的计算复杂度为O(N2)，则计算代价会随着虚拟边界多面体的分辨率提升而迅速增加。表 1总结了本文梳理的目前主要小天体引力场模型特点，可见尚未有同时满足4方面要求的模型。

表1 不同小天体引力场模型特点比较Table 1 Comparison of the characteristics of gravitational field models of different small celestial bodies

3.2 基于数据驱动方法的计算效率

其他方法则是针对如何提升小天体引力场计算效率的问题提出的，这类方法包括机器学习法、八叉树插值法和复合引力场逼近法等。基于数据驱动的一类方法虽然计算效率得到大幅提升，但节点处的“真值”或训练集的采样数据依然依赖某个具体的引力场表征模型给出。因此，该类方法并没有从基本层面解决小天体引力场表征问题。未来小天体引力场表征方法的发展趋势是除保证方法的全局性、精确性和参数易求解性之外，还应保证方法的高效性。

3.3 双小行星系统动力学模型的高精度与高效率计算

双小行星系统对天体力学的科学探索和行星防御的工程实践都具有重要的价值，近些年来成为了深空探测领域的热点目标。双星系统的引力场建模技术能为双星探测任务设计与稳定性分析提供动力学仿真环境。其中，研究两个小天体的姿轨运动是双星系统引力场建模的基础。根据双星系统的运动学方程可知，引力相互势的建模是解决全二体问题的关键。然而，小行星不规则的形状、不均匀的质量分布和双星系统独特的构型，给双星系统引力相互势的计算带来了诸多挑战。

早期的研究一般基于单星的引力场模型，对另一颗星进行适当的简化，以得到闭合形式的引力相互势。这类方法建立了双星系统的相互势与系统模型参数（如形状尺寸参数、相对距离参数、相对指向参数、质量参数等）之间的定量关系。有助于进行双星系统能量、稳定构型方面的初步分析。然而，对小天体的形状或者引力场进行简化损失了对系统描述的精度。级数展开类的方法在保留对小天体不规则形状和不均匀质量分布描述的基础上，对相互势进行级数展开，相互势一般被表示为无穷级数的形式。此类方法的误差来源于级数的截断误差，随着两个小天体之间的距离减小，级数的截断误差增大，且此类模型的计算量随着级数阶数的提高显著增大。为了克服级数展开类方法收敛困难的问题，有限元方法在双星建模中的应用，为引力相互作用的计算提供了全新的思路。通过有限元划分，将两个小行星之间的相互势表示为有限元微元对之间相互势叠加的形式。通过细致的有限元划分，有限元方法能够准确地描述两个小行星的不规则形状和不均匀质量分布，并且截断误差可控。但是，该方法的计算负担会随着有限元的数量增加而增加。

基于单个小天体的引力场建模技术研究基础，未来双小行星系统附近的引力场建模技术的发展应主要关注于双小行星系统的动力学研究。在双小行星系统的动力学建模方面，已有的研究主要关注于双星系统的平衡构型以及姿轨运动演化。然而，在这些方法中，模型的计算精度和计算效率往往是矛盾的。在未来的双星系统引力场建模技术和双星系统动力学建模研究中，应主要关注以下几个方面。

1）对小天体不规则形状和不均匀质量分布的准确描述。根据已探明的小天体可知，小天体的形状不规则且质量分布不均匀，存在风化层等特殊结构。由于双小行星系统独特的动力学，双星系统的演化受到小天体质量分布和形状的影响较为显著，未来的模型应能够考虑小天体特殊的形态和物质组成。

2）兼顾计算精度和计算效率。在现有的方法中，模型的计算精度和计算效率往往是矛盾的，较高的计算精度常常会导致沉重的计算负担。发展兼顾计算精度和计算效率的方法，能够提升在双星系统递推方面的实用性，以实现对双星系统的长期演化递推和双星系统附近引力场建模的有效应用，例如，双星系统附近粒子的运动演化规律研究等。

4 结束语

本文介绍了单个小天体的引力场建模技术和双小行星系统附近的引力场建模技术。目前对目标小天体附近的引力场建模技术研究已经形成了比较成熟的理论体系，但是也仍有一些理论和技术问题尚待解决或完善。发展兼顾计算精度和计算效率、具备全局表征能力的建模技术仍是该领域的研究重点。未来将小天体探测的实际数据与现有理论方法进行对比，能够在验证现有理论的基础上，发现全新的科学、技术问题，推动小天体引力场建模方法和双小行星系统动力学建模方法的发展。