張振江，崔祜濤，任高峰

（哈爾濱工業(yè)大學(xué) 深空探測基礎(chǔ)研究中心，哈爾濱 150080）

0 引言

近年來，小行星和彗星探測任務(wù)越來越多地受到人們的重視。這些已經(jīng)發(fā)射或計劃中的航天任務(wù)促使人們對于一個嶄新的天體力學(xué)領(lǐng)域——繞不規(guī)則形狀小天體的衛(wèi)星軌道動力學(xué)的探索。這個極富挑戰(zhàn)性的問題正吸引著越來越多的科學(xué)家對其進行研究[1-2]。小行星極不規(guī)則的形狀是造成其衛(wèi)星動力學(xué)特性與近球形大行星或天然衛(wèi)星不同的主要原因之一。故對不規(guī)則形狀小行星引力勢能的建模就成為研究繞小行星衛(wèi)星軌道動力學(xué)特性所需要解決的首要問題。小行星引力勢能的建模方法很多[3]，這些方法可以分為數(shù)值法和解析法兩大類。

數(shù)值法的主要思想是用一個多面體或多個特定形狀質(zhì)元（如球體或立方體）的組合來逼近小行星的形狀，再通過積分變換將引力勢能中的三重積分轉(zhuǎn)化為可計算的曲線和曲面積分。此類方法充分考慮了小行星的不規(guī)則形狀，利用了地面天文觀測和掠飛任務(wù)所拍攝的小行星的圖像信息，因而具有很高的精度。數(shù)值法的缺點在于它沒有解析的形式，只能求出給定位置的引力勢能和引力加速度值，無法分析各階攝動項的大小及其對衛(wèi)星軌道的影響。因此，數(shù)值法只能在數(shù)值仿真中使用而無法應(yīng)用于軌道設(shè)計。

解析法的主要思想是用級數(shù)展開式直接逼近引力勢能，這類方法具有解析的表達形式，算法簡單，各階攝動大小明顯。因此解析法廣泛應(yīng)用于人造衛(wèi)星軌道動力學(xué)和天體力學(xué)，尤其在分析攝動對軌道的影響、設(shè)計某些使命軌道時，只能使用此類方法。解析法的最大缺點是級數(shù)項的系數(shù)難于確定，例如在應(yīng)用最廣泛的球諧函數(shù)模型中，各階次球諧系數(shù)Clm和Slm都是由衛(wèi)星軌道數(shù)據(jù)通過導(dǎo)航算法“事后”解算出來的[4]。而對于還沒有繞飛任務(wù)的小行星來說，因為沒有衛(wèi)星軌道數(shù)據(jù)作為觀測量，故其準(zhǔn)確的球諧系數(shù)是無法得到的。以往學(xué)者的做法是將小行星簡化成三軸橢球體，再計算其各階次的球諧系數(shù)，以此作為小行星引力場的球諧系數(shù)進行分析和計算[5-6]。顯然，將形狀十分復(fù)雜的小行星簡化成球體或三軸橢球體這是不精確的。

針對在任務(wù)發(fā)射前小行星引力場建模過程中出現(xiàn)的三軸橢球體模型精度過低、高精度球諧系數(shù)由于沒有觀測量又無法獲取的問題，本文提出了一種基于多面體模型的不規(guī)則形狀小行星引力場球諧系數(shù)的求取方法。其思路是先由多面體模型方法重構(gòu)出小行星外部引力場的引力勢能情況，并以此作為求解引力場球諧系數(shù)的虛擬觀測量，再根據(jù)球諧系數(shù)與引力勢能間的關(guān)系，采用最小二乘方法解算出各階次球諧系數(shù)。

1 多面體模型方法

1.1 模型的構(gòu)建

1996年，R. A. Werner提出用多面體來逼近小行星的形狀進而求其引力勢能的方法[7-8]，即多面體模型方法。在小行星引力勢能建模的數(shù)值方法中，多面體模型應(yīng)用最為廣泛。如圖1所示，多面體模型就是用一個表面由一系列三角形構(gòu)成的多面體來逼近小行星形狀，這樣只要求出多面體的引力勢能分布，就可以知道小行星引力勢能的分布情況了。

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

1.2 引力勢能計算方法

任意形狀中心天體的引力勢能可由

式(1)為一個三重積分。對于形狀不規(guī)則的物體，式(1)是不可積的；而對于勻質(zhì)的多面體來說，可以采用積分變換的方法，將此三重積分轉(zhuǎn)換成第二型曲線積分和曲面積分。而這些積分可以寫成有限項之和的形式，如式(2)所示[9]。

以上即為多面體模型方法的全部公式，只要給定小行星固連坐標(biāo)系中一個點的位置矢量，就可以由以上公式計算出該點的引力勢能。下面介紹以引力勢能為虛擬觀測量求取球諧系數(shù)的方法。

2 球諧系數(shù)的計算

由多面體模型方法計算出小行星的引力勢能后，可以將引力勢能作為虛擬觀測量，通過引力勢能與球諧系數(shù)之間的關(guān)系采用類似導(dǎo)航算法中的最小二乘法來求解球諧系數(shù)。

2.1 引力勢能與球諧系數(shù)間的關(guān)系

引力勢能的球諧函數(shù)模型為

式中： GM為小行星的引力常量；R為檢驗點到中心天體質(zhì)心的距離；a為小行星參考球體的半徑；Plm(sin ?)為sin?的締合勒讓德多項式，當(dāng)m=0時，即退化為一般的勒讓德多項式；Clm和Slm為球諧系數(shù)。

略去式(3)中的高階項，可得到下述形式的引力勢能計算公式：

式(4)中n和l為所取球諧系數(shù)的階數(shù)和次數(shù)，η為截斷誤差。由此可知，引力勢能與各階次球諧系數(shù)之間成線性關(guān)系，其系數(shù)項為 Plmsin ?c os mλ和Plmsin? s inmλ。這樣，在已知引力勢能分布的情況下，我們就可以通過這個關(guān)系求得球諧系數(shù)的具體值。

2.2 構(gòu)造線性方程組

式(5)中，方程未知量的個數(shù)為 N( N + 2)（Sl0無意義）。理論上只要在小行星的引力場中取N( N + 2)個不同的檢驗點就可以構(gòu)造一個N( N + 2)維的線性方程組，其形式如式(6)所示：

上式為一個線性方程組，可以寫成Ax=b的形式。解此方程組即可求出球諧系數(shù)Clm和Slm。

2.3 方程組的解法

理論上只要按式(6)構(gòu)造一個 N( N + 2)維的線性方程組就可以唯一地解出球諧系數(shù)Clm和Slm，但實際發(fā)現(xiàn)，由于方程系數(shù)中存在這樣的因子，使得不同檢驗點得到的方程可能是線性相關(guān)的。這使得原給定方程組變成欠定方程組，沒有唯一解。有時雖然各個方程是線性無關(guān)的，但矩陣A的條件數(shù)很大，即矩陣A為病態(tài)，這時方程的解會因存在很大的誤差而失去意義。

由線性方程理論可知，形如Ax=b的超定方程組一般有3種解法，分別是：

1）將原方程轉(zhuǎn)化為正則方程 ATA x=ATb，并由此解得 x =(ATA)-1ATb，此方法雖然計算簡單，但其缺點是矩陣 (ATA)的條件數(shù)過大，因而結(jié)果精度比較低。

2）由奇異值分解理論求得廣義逆矩陣 A+，并由x =A+b求得方程的解，此方法可以獲得可靠的方程解，即使矩陣A發(fā)生列秩虧損時，依然可以給出最小范數(shù)的最小二乘解。缺點是速度比較慢。

3）對原方程進行Householder變換，從而直接求出方程的解。此方法可靠性稍遜色于第二種方法，但速度較快，在矩陣A發(fā)生列秩虧損時，所給出的最小二乘解具有最少的非零元素。本文采用的就是這種方法。

關(guān)于方程數(shù)目的選擇問題，理論上方程數(shù)目越多，解的結(jié)果就越理想，但同時計算速度也會越慢，因此我們需要在計算速度和精度之間做一個權(quán)衡。通過大量計算發(fā)現(xiàn)：一般方程數(shù)目取為未知數(shù)數(shù)目的2 3倍時即可獲得較高的精度。當(dāng)然，解的精度除了與方程數(shù)目有關(guān)外，還與中心小天體形狀的復(fù)雜程度有關(guān)。若中心小天體的形狀較復(fù)雜，應(yīng)適當(dāng)增加方程數(shù)目以保證解的精度。

3 對比驗證

本部分包含兩個算例：第一個算例應(yīng)用本文方法計算一個給定三軸橢球體的球諧函數(shù)，計算結(jié)果與真實值相比較，以驗證算法的正確性及算法精度與多面體面數(shù)之間的關(guān)系；第二個算例采用本文方法計算Eros 433小行星的球諧系數(shù)，其結(jié)果與傳統(tǒng)的三軸橢球體方法和由NEAR探測器軌道數(shù)據(jù)得到的結(jié)果相比較，以驗證本文所采用的方法其優(yōu)點和不足。

3.1 三軸橢球體驗證

勻質(zhì)三軸橢球體球諧系數(shù)的真實值可由如下解析公式給出[10]：

1）當(dāng)l , m 為所有正整數(shù)時，Slm=0；

2）當(dāng)l或 m 為奇數(shù)時，Clm=0；

3）當(dāng)l , m為其他情況時：

式(7)中：a為參考球體的半徑；δ0m為克羅內(nèi)克符號，其具體值為：

在本算例中，分別取三軸橢球體的3個半軸長為16 km、8 km和6 km。

圖2為算例中使用的三軸橢球體的多面體模型，應(yīng)用本文方法分別取多面體的外表面數(shù)目N為400和20 000進行計算。所得計算結(jié)果與由公式(7)計算得到的真實球諧系數(shù)相比較，結(jié)果如表1所示。

圖2 三軸橢球體多面體模型Fig. 2 Polyhedral model of tri-axial ellipsoid

表1 三軸橢球體球諧系數(shù)Table 1 Spherical harmonic coefficients of tri-axial ellipsoid

續(xù)表1

由表1結(jié)果可知：

1）本文提出的方法可以準(zhǔn)確地求出勻質(zhì)三軸橢球體的球諧系數(shù)，且結(jié)果與真實值（解析方法計算的結(jié)果）之間的誤差很小。如在表1的數(shù)據(jù)中，當(dāng)多面體的面數(shù)取為20 000時，最大誤差減小為0.0923%。

2）本方法的誤差隨多面體外表面數(shù)目的增加而減小，如N=400時，最大誤差為2.5%；而當(dāng)N增大到20 000時，最大誤差減小為0.0923%。

3.2 Eros 433算例分析

1998年12月24日，NEAR探測器以1 km/s的速度，從距離Eros 433小行星4 100 km處飛過。在此前后對其的大小、形狀、有無磁場和衛(wèi)星等進行了觀測。整個掠飛過程中，NEAR的多色照相機拍攝了222張Eros 433的照片，可分辨其表面小至500 km范圍內(nèi)的細節(jié)。天文學(xué)家根據(jù)這些照片信息建立了Eros 433小行星的三維模型。此模型形狀數(shù)據(jù)可以由http://www.psi.edu/pds/resource/nearmod.html網(wǎng)站得到。本文采用的形狀模型由64 800個數(shù)據(jù)點組成，每個數(shù)據(jù)點由小行星表面一點的經(jīng)度、緯度和到小行星中心的距離構(gòu)成。將這些數(shù)據(jù)點從球坐標(biāo)系轉(zhuǎn)化至笛卡爾坐標(biāo)系下，經(jīng)三角剖分后即可得到小行星的多面體模型。如圖3所示。

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

計算時，首先采用傳統(tǒng)方法將其簡化為三軸橢球體，簡化得到的3個半軸長分別為16 km、8 km和6 km，將這3個參數(shù)代入公式(7)計算得出球諧系數(shù)。再應(yīng)用本文方法計算出其球諧系數(shù)，最后將二者結(jié)果與由NEAR探測器環(huán)繞軌道數(shù)據(jù)解算的球諧系數(shù)相比較，結(jié)果見表2。

表2 Eros 433引力場球諧系數(shù)Table 2 Gravity spherical harmonic coefficients of Eros 433

由表2結(jié)果可知：

1）與將小行星近似為三軸橢球體的方法相比，本文提出的方法可以大幅度提高結(jié)果的精度，例如采用三軸橢球體方法計算的C20與飛行器軌道數(shù)據(jù)解算的結(jié)果相差17.4%，而采用本文的方法可以將這一誤差減小到0.23%。

2）本方法計算結(jié)果與真實值尚存在一定偏差，例如前4階次球諧系數(shù)中相差最大的S22項與實際值相差6%，其他階次的球諧系數(shù)也存在著一定差異。這些差異主要來自兩個方面：第一，多面體模型與小行星的真實形狀并不完全相同，二者之間存在著一定的誤差，我們可以通過增加多面體外表面數(shù)目的方法來減小這一誤差；第二，多面體模型方法假定小行星是勻質(zhì)的，但真實的小行星密度并不是處處相同的，這也會造成結(jié)果的誤差。而這種誤差是方法中沒有建模的，所以無法消除。這也是本方法的局限所在。

4 結(jié)束語

本文提出的多面體模型方法可以在任務(wù)發(fā)射前得到比較精確的不規(guī)則形狀小行星引力場的球諧系數(shù)，因為其充分利用了天文觀測或掠飛任務(wù)所拍攝的圖像信息，故其精度要比三軸橢球體模型方法有大幅度提高。本文方法計算結(jié)果的精度隨多面體模型外表面的數(shù)目增加而提高。對于多面體模型與小行星形狀不完全一致造成的誤差，可通過增加多面體的面數(shù)減小誤差；而對于小行星質(zhì)量不均勻所引起的結(jié)果誤差，因多面體模型中無法對密度分布建模，故無法消除或減小這種誤差。

（

）

[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