崔祜濤，張振江，余 萌

(哈爾濱工業(yè)大學(xué) 深空探測基礎(chǔ)研究中心，150080 哈爾濱，river18202@gmail.com)

形狀不規(guī)則是小行星與類球形大行星的主要區(qū)別之一，不規(guī)則形狀使得小行星的引力場非常復(fù)雜.這不但給研究衛(wèi)星在小行星引力場中的運(yùn)動規(guī)律帶來了非常大的困難，也使得繞小行星衛(wèi)星軌道的設(shè)計(jì)問題變得異常復(fù)雜.因此，建立小行星引力場模型不但是小行星探測中所要解決的關(guān)鍵問題，更是研究小行星衛(wèi)星軌道動力學(xué)和設(shè)計(jì)小行星衛(wèi)星軌道的基礎(chǔ).傳統(tǒng)方法很難對形狀如此復(fù)雜的小行星引力場建模，隨著計(jì)算機(jī)科學(xué)的發(fā)展，特別是近20 年來計(jì)算機(jī)3D 建模與仿真技術(shù)的提高，使得精確地建立小行星引力場成為可能［1］.到目前為止，常用的小行星引力場建模方法按其本質(zhì)可分為兩大類:第一類是采用無窮級數(shù)逼近小行星的引力勢函數(shù)，常用的有球諧函數(shù)法［2］和橢球諧函數(shù)法［3］;第二類方法是用可計(jì)算引力勢能的模型來逼近小行星的形狀，然后求得近似模型的引力勢函數(shù).這類方法主要有三軸橢球體法［4］和多面體模型法.地球引力場建模常用第一類方法，但對于小行星而言，應(yīng)用第一類方法存在兩個(gè)問題，一是缺乏環(huán)繞軌道數(shù)據(jù)來確定(橢)球諧系數(shù)，二是在接近小行星表面(布里淵(橢)球［3］以內(nèi))時(shí)，諧函數(shù)方法的計(jì)算結(jié)果會存在較大誤差甚至發(fā)散.三軸橢球體方法的缺點(diǎn)是橢球體與小行星形狀相似度差，精度不高.

為解決上述方法存在的問題，R.A.Werner［5］和P.S.RYAN［6］提出用多面體來逼近小行星的形狀進(jìn)而求其引力勢能的方法，2002 年Miller 等［7］研究了應(yīng)用地面天文觀測確定小行星的形狀、自旋情況、類別和密度等參數(shù)的方法.2006 年E.G.Fahnestock 等［8］應(yīng)用多面體模型方法對不規(guī)則形狀雙星系統(tǒng)的繞飛軌道進(jìn)行仿真計(jì)算.2009 年胡維多［9］系統(tǒng)地分析了近小行星區(qū)域的動力學(xué)環(huán)境及其對環(huán)繞小行星飛行器軌道的影響情況.為解決小行星探測任務(wù)前無法獲得高精度引力場球諧系數(shù)的問題，筆者在2010 年提出一種基于多面體模型的不規(guī)則形狀小行星引力場球諧系數(shù)的求取方法［10］.Shao Wei［11］提出一種利用三維散亂點(diǎn)對小天體表面進(jìn)行三角剖分，并建立小行星多面體模型的方法.

1 小行星形狀與多面體模型

圖1 為美國近地小行星交會任務(wù)(Near Earth Asteroid Rendezvous Mission，NEAR)的蘇梅克探測器(Shoemaker spacecraft)在2000 年2 月12 日從距離Eros433 小行星1 800 km 之外所拍攝的一系列照片中的6 幅.由圖1 可知，Eros433 小行星的形狀極不規(guī)則，事實(shí)上，太陽系中絕大多數(shù)小行星的形狀都不規(guī)則.為了盡可能精確地描述小行星的形狀，可以使用1 個(gè)表面由一系列三角形構(gòu)成的多面體來逼近小行星，這個(gè)多面體即為小行星的多面體模型.仿真證明，當(dāng)多面體由5 184 個(gè)面構(gòu)成時(shí)，已經(jīng)可以較好地逼近Eros433 小行星的形狀.在本文的分析計(jì)算中，應(yīng)用的多面體模型包含49 152 個(gè)面.讀者可以在NASA 的Planetary Data System 中獲得更多的小行星多面體模型.

多面體模型可以通過分析小行星照片獲得［12］，這些照片可以是地面天文臺拍攝的［13］，也可以是掠飛過程中衛(wèi)星拍攝的.對于已經(jīng)探測的小行星，探測器的激光雷達(dá)可以提供更為精確的多面體模型［14］.多面體模型可以充分利用地面天文觀測或掠飛任務(wù)所拍攝的小行星圖像信息，因而在任務(wù)設(shè)計(jì)階段就可以取得較高的精度.

圖1 Eros433 全貌照片

2 多面體模型引力勢能的計(jì)算

小行星引力場建模實(shí)質(zhì)為求取引力場中任意檢驗(yàn)點(diǎn)P(x，y，z)處的引力加速度函數(shù)F(x，y，z)的過程.而引力加速度可由引力勢能V(x，y，z)得到，二者的關(guān)系為

式(1)表明小行星引力場中某一點(diǎn)的引力加速度F(x，y，z)等于該點(diǎn)引力勢函數(shù)對位置的偏導(dǎo)數(shù)，即引力勢函數(shù)的梯度.這樣，引力場建模問題就轉(zhuǎn)化為求小行星引力場中任意檢驗(yàn)點(diǎn)引力勢函數(shù)的問題.

圖2 為引力勢能函數(shù)計(jì)算的矢量圖.

圖2 引力勢能函數(shù)的計(jì)算

如圖2 所示，檢驗(yàn)點(diǎn)P 在小行星固連坐標(biāo)系下的位置矢量為R' =［x y z］T，大小用R'表示.欲求此點(diǎn)的引力勢能，可以將中心天體分解成若干小體積微元，則每個(gè)體積微元可以看成1 個(gè)質(zhì)點(diǎn)，設(shè)1 個(gè)質(zhì)量記為dm 的體積微元S 在小行星固連坐標(biāo)系內(nèi)的位置矢量為R=［ξ η ζ］T，大小用R 表示.則由檢驗(yàn)點(diǎn)到該體積微元的矢量為r=［x-ξ y-η z-ζ］T，大小用r 表示，檢驗(yàn)點(diǎn)處的引力勢能可由下面三重積分定義［15］:

應(yīng)用高斯公式，式(4)中的三重積分可以轉(zhuǎn)化為曲面積分，即

對于多面體模型而言，式(5)中的曲面積分可以寫成如下形式:

式中eedge表示邊，re為多面體邊e 上任意一點(diǎn)到檢驗(yàn)點(diǎn)的矢量，且

其中Ee為3×3 矩陣;^nA為面A 的外法線方向矢量為在面A 內(nèi)的邊e 的外法線方向矢量;re1、re2為檢驗(yàn)點(diǎn)到邊的兩個(gè)端點(diǎn)的距離;e12為邊e 的長度;^nf為面f 的外法線方向矢量;Ff為3×3矩陣.

圖5(a)～圖5(c)分別是軸承3種運(yùn)行狀態(tài)下某一采樣樣本的雙譜時(shí)頻圖.圖5中可以看出,雖然軸承在不同故障運(yùn)行狀態(tài)下雙譜圖存在很大差異,可以為后續(xù)SVDD診斷所需的特征提取提供很好的素材.但由于滾動軸承復(fù)合故障運(yùn)行狀態(tài)下，各個(gè)單個(gè)故障信號之間的相互干擾極有可能出現(xiàn)耦合現(xiàn)象,很難直接在其雙譜時(shí)頻圖上提取有效的故障特征信息.

上面就是多面體模型計(jì)算引力勢能的全部公式，下面介紹一下引力加速度和其梯度的計(jì)算公式.

式(9)是用來計(jì)算引力加速度的，式(10)為引力加速度的梯度矩陣.式(10)中的▽2V(R')稱為引力場拉普拉斯算子，它等于引力梯度矩陣的跡.根據(jù)拉普拉斯算子的值可以方便的判斷檢驗(yàn)點(diǎn)是否位于多面體模型的外側(cè).判斷準(zhǔn)則如下式所示.

3 Eros433 引力場計(jì)算與分析

`本文計(jì)算中使用的Eros433 小行星多面體模型來源于NASA 的Planetary Data System.數(shù)據(jù)共包含4 個(gè)不同精度的多面體模型，面的個(gè)數(shù)分別為49 152，196 608，786 432，3 145 728.為減小計(jì)算量，節(jié)約時(shí)間，本文應(yīng)用49 152 個(gè)面的多面體模型進(jìn)行計(jì)算.

3.1 引力場計(jì)算

NASA 提供的數(shù)據(jù)包中的內(nèi)容共分兩部分，第一部分為多面體表面上點(diǎn)的位置列表，如數(shù)據(jù)包中第1 個(gè)點(diǎn)為“1 8.332 72E+00-4.809 02E+00-4.784 54E+00”，其中1 表示該點(diǎn)的編號，后面3 個(gè)數(shù)為該點(diǎn)在小行星固連坐標(biāo)系中的位置分量，單位為km;第二部分為組成多面體的三角形索引，如其中1 個(gè)三角形表示為“1 1 105 101”其中第1 個(gè)1 表示該三角形編號，后面3 個(gè)數(shù)表示三角形的3 個(gè)頂點(diǎn)的編號，查詢第一部分即可確定該三角形在小行星固連坐標(biāo)系中的位置.

根據(jù)此數(shù)據(jù)包提供的數(shù)據(jù)，結(jié)合上文給出的公式，即可計(jì)算多面體模型外任意一點(diǎn)的引力勢能和引力加速度.計(jì)算中首先應(yīng)用式(11)測試檢驗(yàn)點(diǎn)是否位于多面體外側(cè)，對于位于外側(cè)的檢驗(yàn)點(diǎn)，應(yīng)用式(7)和式(8)求其引力勢能和引力加速度.小行星的密度取為2 670 kg/m3，計(jì)算結(jié)果如圖3 ～圖10 所示.在引力加速度的等高線圖中，顏色用來表示引力加速度的大小，箭頭表示引力加速度的方向.

圖3 ～4 分別為小行星固連坐標(biāo)系的XOY，YOZ，ZOX 三個(gè)坐標(biāo)平面內(nèi)的引力勢能與引力加速度的分布情況，由這6 幅圖像可知，Eros433 小行星引力場中引力勢能和引力加速度的分布都不規(guī)則.小行星表面的引力加速度分布也極不均衡，最大加速度與最小加速度相差近三倍，而且引力加速度的方向也較復(fù)雜，南半球的表面平坦區(qū)域，引力大體上與小行星表面垂直，而在北半球的巨大凹槽處，引力的方向與小行星表面的夾角最大.

圖3 Eros433 小行星在坐標(biāo)平面內(nèi)引力勢能分布

圖5 ～6 為距離質(zhì)心18 km 的球面上Eros433小行星引力勢能和引力加速度的分布情況.由這兩幅圖像可知，在半徑為18 km 處，引力加速度存在3 個(gè)極大值(圖6 中點(diǎn)A，B，C 周圍)區(qū)域和兩個(gè)極小值(圖6 中點(diǎn)D，E 周圍)區(qū)域.引力加速度最大的區(qū)域大約位于西經(jīng)165°的赤道(點(diǎn)A)附近.兩個(gè)極小值區(qū)域位于西經(jīng)75°的南半球和東經(jīng)100°的北半球附近.在這些重力異常區(qū)，衛(wèi)星的軌道會受到很大的影響，軌道根數(shù)可能在短時(shí)間內(nèi)發(fā)生較大的變化，這與繞大行星的衛(wèi)星軌道有很大區(qū)別［16］.

圖4 Eros433 小行星在坐標(biāo)平面內(nèi)引力加速度分布

圖5 R=18 km 的引力勢能分布

3.2 對比分析

上文給出了應(yīng)用多面體模型方法計(jì)算得到的Eros433 小行星引力場的分布情況，為檢驗(yàn)算法的準(zhǔn)確性，將上述計(jì)算結(jié)果與JPL 實(shí)驗(yàn)室公布的Eros433小行星16 階引力場模型進(jìn)行對比.該引力場模型根據(jù)蘇梅克探測器對小行星近一年的環(huán)繞飛行的軌道數(shù)據(jù)和探測器上攜帶的激光高度計(jì)的觀測值綜合計(jì)算得到.由于在小行星表面附近(布里淵球內(nèi))球諧函數(shù)發(fā)散，因此無法計(jì)算貼近小行星表面區(qū)域的引力勢能，所以這里只比較兩種方法計(jì)算的距小行星質(zhì)心18 km 處的引力勢能值之間的差.其結(jié)果如圖7 所示.

圖6 R=18 km 的引力加速度分布

圖7 R=18 km 時(shí)兩種方法計(jì)算結(jié)果的偏差

由圖7 可知兩種方法計(jì)算得到的引力勢能誤差非常小，經(jīng)計(jì)算表明，最大相對誤差出現(xiàn)在小行星的長軸一端附近，對應(yīng)于圖6 的A 點(diǎn)附近區(qū)域，其最大誤差為1.52%，其余區(qū)域的最大誤差均不超過0.5%.這充分說明本文方法具有較高的計(jì)算精度.為清楚的描述計(jì)算誤差與小行星形狀間的關(guān)系，將誤差曲面投影到Eros433 表面，可以得到小行星表面的誤差分布圖，如圖8 所示.

圖8 計(jì)算誤差在小行星表面的分布

3.3 偽勢能與平衡點(diǎn)

除了不規(guī)則的形狀外，小行星的自旋對繞其運(yùn)動的衛(wèi)星軌道亦有很大的影響，這兩個(gè)因素相互疊加，使得小行星衛(wèi)星的軌道動力學(xué)變得異常復(fù)雜.在小行星固連坐標(biāo)系中，衛(wèi)星的運(yùn)動方程為［21］

其中［x y z］T為衛(wèi)星在小行星固連坐標(biāo)系中的位置矢量;ωT為小行星的自旋角速度，Eros433的自旋角速度為ωT=3.311 7×10-4rad/s;V 為前文建模的引力勢函數(shù).運(yùn)動方程(12)存在如下形式積分［17］:

其中J 稱為Jacobi 積分常數(shù);0.5(˙x2+˙y2+˙z2)為飛行器的動能，由式(13)右端的前兩項(xiàng)的和可以定義飛行器的偽勢能，其形式為

偽勢能描述了旋轉(zhuǎn)坐標(biāo)系中的引力場分布情況，因?yàn)榭紤]了小行星自身旋轉(zhuǎn)對引力場的影響，偽勢能可以比引力勢能更好的描述小行星的引力場狀況.圖9 為Eros433 小行星赤道平面內(nèi)的偽勢能曲面，圖10 為偽勢能的等高線圖.由式(13)可知偽勢能U 相當(dāng)于動能為零時(shí)的jacobi 積分常數(shù)J，因此，圖10 又稱為零速度曲線圖.

圖9 Eros433 小行星偽勢能曲面

圖9 ～11 反映了赤道平面內(nèi)的偽勢能分布情況，由圖9 可知，在小行星外側(cè)，隨著與小行星距離的增大，偽勢能先減小后增大，在Y 軸附近，存在兩個(gè)偽勢能的極小值點(diǎn)(對應(yīng)于圖10 中的E3，E4).而在X 軸附近，偽勢能曲面存在兩個(gè)鞍狀平衡點(diǎn)(對應(yīng)于圖10 中的E1，E2).這4 個(gè)點(diǎn)稱為平衡點(diǎn)，由圖11 可知，在平衡點(diǎn)附近，偽勢能分布情況復(fù)雜，因此在此區(qū)域內(nèi)運(yùn)動的衛(wèi)星軌道受其影響表現(xiàn)出強(qiáng)烈的非線性，這些軌道往往是不穩(wěn)定的［18］.E1～E4四個(gè)平衡點(diǎn)的性質(zhì)與圓型限制性三體問題中的動平衡點(diǎn)性質(zhì)相似.故在其周圍也存在周期軌道和擬周期軌道，筆者在文獻(xiàn)［19］中研究了E1、E2兩個(gè)平衡點(diǎn)附近的軌道特性和軌道控制問題.Scheeres 在文獻(xiàn)［20］研究了平衡點(diǎn)附近周期軌道的確定與計(jì)算問題.

圖10 Eros443 零速度曲線和平衡點(diǎn)

圖11 平衡點(diǎn)E2 點(diǎn)附近偽勢能的分布

4 結(jié) 論

本文介紹了基于形狀逼近的小行星引力場建模方法——多面體模型法.給出了Eros433 小行星的多面體模型，并推導(dǎo)了多面體模型方法中引力勢能、引力加速度和拉普拉斯算子的計(jì)算公式.應(yīng)用 NASA 的 Planetary Data System 提供的Eros433小行星多面體模型分析了其引力場分布情況，繪制了引力勢能和引力加速度的分布圖.考慮到小行星旋轉(zhuǎn)對引力場的影響，本文還分析了小行星赤道平面內(nèi)偽勢能曲面與零速度曲線分布情況.由以上分析可知:

1)Eros433 小行星表面的引力加速度分布極不均衡，最大加速度與最小加速度相差近三倍，而且引力加速度的方向也較復(fù)雜，南半球的表面平坦區(qū)域，引力大體上與小行星表面垂直，北半球的巨大凹槽處，引力的方向與小行星表面存在很大的夾角.

2)在距離小行星中心18 km 處，引力加速度存在3 個(gè)極大值和2 個(gè)極小值.其中引力加速度最大的區(qū)約位于西經(jīng)165°的赤道附近，2 個(gè)極小值分別位于南半球的西經(jīng)75°和北半球的東經(jīng)100°區(qū)域.

3)在自旋和不規(guī)則形狀的共同影響下，Eros433小行星引力場中存在4 個(gè)平衡點(diǎn)，其中2 個(gè)平衡點(diǎn)位于Y 軸附近，對應(yīng)于偽勢能的極小值;另外2 個(gè)位于X 軸附近，與之對應(yīng)的偽勢能曲

面呈馬鞍面狀分布.

［1］ZUBER M T，SMITH D E，CHENG A F.The shape of 433 Eros from the NEAR-shoemaker laser rangefinder［J］.Science，2000，289(5487):2097-2100.

［2］STEFANO C，SUSANNA M.Methods for computing the potential of an irregular，homogeneous，solid body and its gradient［C］//AIAA/AAS Astrodynamics Specialist Conference.Denver，CO:AIAA，2000:4-17.

［3］GARMIER R，BARRIOT J P.Ellipsoidal harmonic expressions of the gravitational potential:theory and applications［J］.Celestial Mechanics and Dynamical Astronomy，2001，79(4):235-275.

［4］SCHEERES D J.Dynamics about uniformly rotating triaxial ellipsoids.Applications to Asteroids［J］.Icarus，1994，110(2):225-238.

［5］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(3):313-44.

［6］RYAN S P，ROBERT A W，BHASKARANZ S.Estimating small-body gravity field from shape model and navigation data［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Honolulu，Hawaii:AIAA，2008，AIAA-2008-6603.

［7］MILLER J K，KONOPLIV A S，ANTREASIAN P G.Determination of shape，gravity，and rotational state of asteroid 433 Eros［J］.Icarus，2002，155(1):3-17.

［8］FAHNESTOCK E G，TAEYOUNG L，MELVIN L.Polyhedral potential and variational integrator computation of the full two body problem［C］//AIAA/AAS Astrodynamics Specialist Conference and Exhibit.Keystone，Colorado:AIAA，2006，AIAA 2006-6289.

［9］胡維多，SCHEERES D J，向開恒.飛行器近小行星軌道動力學(xué)的特點(diǎn)及研究意義［J］.天文學(xué)進(jìn)展，2009，27(6):152-166.

［10］ZHANG Zhenjiang，CUI Hutao，REN Gaofeng.Modeling for the gravitation potential environment of an irregular-shaped asteroid and the spherical harmonic coefficient estimation［J］.Spacecraft Environment Engineering，2010，27(3):383-388.

［11］SHAO Wei，CUI Pingyuan，CUI Hutao.Physical properties calculation of small body using points triangulations［J］.Journal of Harbin Institute of Technology，2010，42(5):687-691.

［12］MITCHELLA D L，HUDSONB R S，OSTROA S J，et al.Shape of asteroid 433 Eros from inversion of goldstone radar doppler spectra［J］.Icarus，1998，131(1):4-14.

［13］SHEPARD M K，MARGOT Jean-Luc，CHRISTOPHER Magri，et al.Radar and infrared observations of binary near-Earth Asteroid 2002 CE26［J］.Icarus，2006，184(1):198-210.

［14］SHINSUKE A，TADASHI M，NARU H.Mass and local topography measurements of itokawa by hayabusa［J］.Science，2006，312(5778):1344-1347.

［15］劉林.航天器軌道理論［M］.北京，國防工業(yè)出版社，2000:104-106.

［16］LARA M，SCHEERES D J.Stability bounds for three dimensional motion close to asteroids［J］.Journal of the Astronautical Sciences，2002，50(4):389-409.

［17］SCHEERES D J，MILLER J K，YEOMANS D K.The orbital dynamics environment of 433 Eros:a case study for future asteroid missions［R］.Pasadena:Jet Propulsion Laboratory，2003，42(152):1-26.

［18］SCHEERES D J，HSIAO F Y，VINH N X.Stabilizing motion relative to an unstable orbit:applications to spacecraft formation flight［J］.Journal of Guidance，Control，and Dynamics，2003，26(1):62-73.

［19］ZHANG Zhenjiang，CUI Hutao.Orbit dynamics and control in the neighborhood of the asteroid＇s center point［C］//The 3rd International Symposium on Systems and Control in Aeronautics and Astronautics (ISSCAA 2010).Harbin:［s.n.］，2010.

［20］SCHEERES D J.Periodic orbits in rotating second degree and order gravity fields［J］.Chinese Journal of Astronomy and Astrophysics，2008，8(1)，108-118.