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

第27卷第3期航天器环境工程2010年6月SPACECRAFTENVIRONMENTENGINEERING383不规则形状小行星引力环境建模及球谐系数求取方法张振江，崔祜涛，任高峰（哈尔滨工业大学深空探测基础研究中心，哈尔滨150080）摘要：在实施小行星探测任务之前，需要对不规则形状小行星的引力场有一个清楚的认识，以便于进行环绕或着陆小行星的轨道设计。文章提出了一种小行星引力场建模及球谐系数的求取方法。首先，由多面体模型方法重构出小行星外部的引力环境并以此作为虚拟观测量。在此基础上，根据球谐系数与引力势能间的关系，通过求解超定线性方程组得到小行星引力场的各阶次球谐系数。与传统的将小行星近似成三轴椭球体进而计算球谐系数的方法相比，该方法可大幅度提高引力场建模的精度。通过与由NEAR探测器轨道数据解算的Eros433小行星的球谐系数比较表明，其最大误差不超过6%。计算表明，该方法可以为小行星探测任务发射前的轨道设计提供更为精确的数据。关键词：小行星；多面体模型；引力势能；球谐系数中图分类号：V412.4文献标识码：A文章编号：1673-1379(2010)03-0383-06DOI：10.3969/j.issn.1673-1379.2010.03.0240引言近年来，小行星和彗星探测任务越来越多地受到人们的重视。这些已经发射或计划中的航天任务促使人们对于一个崭新的天体力学领域——绕不规则形状小天体的卫星轨道动力学的探索。这个极富挑战性的问题正吸引着越来越多的科学家对其进行研究[1-2]。小行星极不规则的形状是造成其卫星动力学特性与近球形大行星或天然卫星不同的主要原因之一。故对不规则形状小行星引力势能的建模就成为研究绕小行星卫星轨道动力学特性所需要解决的首要问题。小行星引力势能的建模方法很多[3]，这些方法可以分为数值法和解析法两大类。数值法的主要思想是用一个多面体或多个特定形状质元（如球体或立方体）的组合来逼近小行星的形状，再通过积分变换将引力势能中的三重积分转化为可计算的曲线和曲面积分。此类方法充分考虑了小行星的不规则形状，利用了地面天文观测和掠飞任务所拍摄的小行星的图像信息，因而具有很高的精度。数值法的缺点在于它没有解析的形式，只能求出给定位置的引力势能和引力加速度值，无法分析各阶摄动项的大小及其对卫星轨道的影响。因此，数值法只能在数值仿真中使用而无法应用于轨道设计。解析法的主要思想是用级数展开式直接逼近引力势能，这类方法具有解析的表达形式，算法简单，各阶摄动大小明显。因此解析法广泛应用于人造卫星轨道动力学和天体力学，尤其在分析摄动对轨道的影响、设计某些使命轨道时，只能使用此类方法。解析法的最大缺点是级数项的系数难于确定，例如在应用最广泛的球谐函数模型中，各阶次球谐系数lmC和lmS都是由卫星轨道数据通过导航算法“事后”解算出来的[4]。而对于还没有绕飞任务的小行星来说，因为没有卫星轨道数据作为观测量，故其准确的球谐系数是无法得到的。以往学者的做法是将小行星简化成三轴椭球体，再计算其各阶次的球谐系数，以此作为小行星引力场的球谐系数进行分析和计算[5-6]。显然，将形状十分复杂的————————————收稿日期：2009-09-09；修回日期：2009-10-09基金项目：国家自然科学基金（项目编号：10672044）作者简介：张振江（1981－），男，博士研究生，研究方向为小行星引力环境分析以及深空探测器轨道动力学与控制。E-mail：[e-mail]。384航天器环境工程第27卷小行星简化成球体或三轴椭球体这是不精确的。针对在任务发射前小行星引力场建模过程中出现的三轴椭球体模型精度过低、高精度球谐系数由于没有观测量又无法获取的问题，本文提出了一种基于多面体模型的不规则形状小行星引力场球谐系数的求取方法。其思路是先由多面体模型方法重构出小行星外部引力场的引力势能情况，并以此作为求解引力场球谐系数的虚拟观测量，再根据球谐系数与引力势能间的关系，采用最小二乘方法解算出各阶次球谐系数。1多面体模型方法1.1模型的构建1996年，R.A.Werner提出用多面体来逼近小行星的形状进而求其引力势能的方法[7-8]，即多面体模型方法。在小行星引力势能建模的数值方法中，多面体模型应用最为广泛。如图1所示，多面体模型就是用一个表面由一系列三角形构成的多面体来逼近小行星形状，这样只要求出多面体的引力势能分布，就可以知道小行星引力势能的分布情况了。图1小行星216Kleopatra（左）和6489Golevka（右）的多面体模型Fig.1Polyhedralmodelofasteroid216Kleopatra(left)and6489Golevka(right)1.2引力势能计算方法任意形状中心天体的引力势能可由1()d=−∫∫∫KMUGmrr(1)定义，式中：U表示检验点的引力势能；Kr表示检验点在小行星固连坐标系中的位置矢量；r为其幅值大小；G为引力常量。式(1)为一个三重积分。对于形状不规则的物体，式(1)是不可积的；而对于匀质的多面体来说，可以采用积分变换的方法，将此三重积分转换成第二型曲线积分和曲面积分。而这些积分可以写成有限项之和的形式，如式(2)所示[9]。TTedgeface111111ˆˆˆ()dddivdd2222∈∈=−=−=−=−⋅=−⋅+⋅∑∑∫∫∫∫∫∫∫∫∫∫∫KJKJKJJKJJKeeeeffffMVVSefUGmGVGVGSGELGErrρρρρωrrnrrrrr，(2)式中：Kr为检验点在小行星固连坐标系中的位置矢量；r为其幅值大小；ρ为小行星密度；ˆr为Kr的单位矢量，ˆdivr为ˆr的散度；ˆn为面积微元dS的外法线方向矢量；Ker为检验点到边e上任一点的矢量；其他参数分别表示如下。()()TTABA12B21ˆˆˆˆ=+eEnnnn，其中：Aˆn为面A的外法线方向矢量；Ee为3×3矩阵；12121212lneeeeerreLrre++=+−，其中：12,eerr为检验点到所求边的两个端点的距离；12e为边的长度；Kfr为检验点到面f上任一点的矢量；Tˆˆ=fffFnn，ˆfn为面f的外法线方向矢量，Ff为3×3矩阵；()()()()123123123221312arctan⋅×=+⋅+⋅+⋅KKKKKKKKKfrrrrrrωrrrrrrrrr。以上即为多面体模型方法的全部公式，只要给定小行星固连坐标系中一个点的位置矢量，就可以由以上公式计算出该点的引力势能。下面介绍以引力势能为虚拟观测量求取球谐系数的方法。2球谐系数的计算由多面体模型方法计算出小行星的引力势能后，可以将引力势能作为虚拟观测量，通过引力势能与球谐系数之间的关系采用类似导航算法中的最小二乘法来求解球谐系数。2.1引力势能与球谐系数间的关系引力势能的球谐函数模型为()[]M10()1sincossin∞==⎧⎫⎪⎪⎡⎤=+×+⎨⎬⎢⎥⎣⎦⎪⎪⎩⎭∑∑JKlllmlmlmlmGaUPCmSmRRϕλλR，(3)第3期张振江等：不规则形状小行星引力环境建模及球谐系数求取方法385式中：GM为小行星的引力常量；R为检验点到中心天体质心的距离；a为小行星参考球体的半径；(sin)lmPϕ为sinϕ的缔合勒让德多项式，当0m=时，即退化为一般的勒让德多项式；lmC和lmS为球谐系数。式(3)中，大括号中的第一项MGR表示一匀质球体的引力势能，后面的各项表示对匀质球体引力势能的修正。理论上这样一个无穷级数可以逼近任意形状中心天体的引力势能，但无穷级数本身是无法计算的，故在使用中通常按所需精度略去高阶项而取前有限项进行计算。略去式(3)中的高阶项，可得到下述形式的引力势能计算公式：()[]M10()1sincossin==⎧⎫⎪⎪⎡⎤=+×++⎨⎬⎢⎥⎣⎦⎪⎪⎩⎭∑∑JKlnllmlmlmlmGaUPCmSmRRRϕλλη。(4)式(4)中n和l为所取球谐系数的阶数和次数，η为截断误差。由此可知，引力势能与各阶次球谐系数之间成线性关系，其系数项为sincoslmPmλϕ和sinsinlmPmλϕ。这样，在已知引力势能分布的情况下，我们就可以通过这个关系求得球谐系数的具体值。2.2构造线性方程组任给一个检验点(),,=KiRRιιιλϕ，通过多面体模型方法算得此处的引力势函数iU，就可构造一个以lmC和lmS为变量的方程：()()()()1010111111MsincosSsinsincossin1⎡⎤⎡⎤⎡⎤++++=−⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦NiiiiiiNNNNNNiiiaaaRPCPCmmPCNSNURRRGϕλλϕλλ，(5)式(5)中，方程未知量的个数为(2)NN+（0lS无意义）。理论上只要在小行星的引力场中取(2)NN+个不同的检验点就可以构造一个(2)NN+维的线性方程组，其形式如式(6)所示：()()11111010111111McosSsincossin1;1111⎡⎤⎡⎤⎡⎤++++=−⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦NNNNNNNaaaRPCPCmmPCNSNURRRGλλλλ()()22221010111111McosSsincossin1;2222⎡⎤⎡⎤⎡⎤++++=−⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦NNNNNNNaaaRPCPCmmPCNSNURRRGλλλλ##()()1010111111McosSsincossin1⎡⎤⎡⎤⎡⎤++++=−⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦NiiiiiiiiNNNNNNaaaRPCPCmmPCNSNURRRGλλλλ##()()22222222222222221010101111McosSsincossin1⎡⎤⎡⎤⎡⎤++++=−⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦NN+N+N+N+N+N+N+N+NNNNNNaaaRPCPCmmPCNSNURRRGλλλλ。(6)上式为一个线性方程组，可以写成=xbA的形式。解此方程组即可求出球谐系数lmC和lmS。2.3方程组的解法理论上只要按式(6)构造一个(2)NN+维的线性方程组就可以唯一地解出球谐系数lmC和lmS，但实际发现，由于方程系数中存在(sin)nlmaPRϕ⎛⎞⎜⎟⎝⎠这样的因子，使得不同检验点得到的方程可能是线性相关的。这使得原给定方程组变成欠定方程组，没有唯一解。有时虽然各个方程是线性无关的，但矩阵A的条件数很大，即矩阵A为病态，这时方程的解会因存在很大的误差而失去意义。解决上述问题有两种方法：第一种是寻找合适的选择检验点(,,)iiiiRRλϕ=K，使得矩阵A的条件386航天器环境工程第27卷数约为1。第二种是增加方程的数目，使方程组变成超定方程组，再求出此超定方程组的最小二乘解。最小二乘解即为球谐系数lmC和lmS的最优估计值，本文即采用第二种方法进行计算。由线性方程理论可知，形如=xbA的超定方程组一般有3种解法，分别是：1）将原方程转化为正则方程TT=xbAAA，并由此解得T1T()−=xbAAA，此方法虽然计算简单，但其缺点是矩阵T()AA的条件数过大，因而结果精度比较低。2）由奇异值分解理论求得广义逆矩阵+A，并由+=xbA求得方程的解，此方法可以获得可靠的方程解，即使矩阵A发生列秩亏损时，依然可以给出最小范数的最小二乘解。缺点是速度比较慢。3）对原方程进行Householder变换，从而直接求出方程的解。此方法可靠性稍逊色于第二种方法，但速度较快，在矩阵A发生列秩亏损时，所给出的最小二乘解具有最少的非零元素。本文采用的就是这种方法。关于方程数目的选择问题，理论上方程数目越多，解的结果就越理想，但同时计算速度也会越慢，因此我们需要在计算速度和精度之间做一个权衡。通过大量计算发现：一般方程数目取为未知数数目的2～3倍时即可获得较高的精度。当然，解的精度除了与方程数目有关外，还与中心小天体形状的复杂程度有关。若中心小天体的形状较复杂，应适当增加方程数目以保证解的精度。3对比验证本部分包含两个算例：第一个算例应用本文方法计算一个给定三轴椭球体的球谐函数，计算结果与真实值相比较，以验证算法的正确性及算法精度与多面体面