不規則形狀小行星引力環境建模及球諧系數求取方法

張振江，崔祜濤，任高峰

（哈爾濱工業大學 深空探測基礎研究中心，哈爾濱 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