小行星重力場全球及局部球諧系數計算與仿真

郭延寧，李曉宇，馬廣富，崔祜濤

(1.哈爾濱工業大學控制科學與工程系，150001哈爾濱;2.哈爾濱工業大學深空探測基礎研究中心，150001哈爾濱)

近年來，作為深空探測重要組成部分，小行星探測活動從早期的近距離飛越、低空繞軌探測［1］，發展至目前的懸停［2］、表面軟著陸和采樣返回［3］，對于推動科學技術進步和幫助人類了解宇宙發展等有重要科學意義［4］.受質量小、形狀不規則、自旋復雜等因素影響，小行星重力場與傳統近球形大行星或天然衛星重力場存在極大差異.因此高精度的重力場建模不僅是研究設計小行星衛星軌道所需要解決的首要問題，也是小行星探測任務的主要科學目標之一.

小行星重力場建模方法的研究已經有一百多年歷史，方法眾多［5］，目前最為廣泛使用的是球諧函數模型［6］與多面體模型［7］.球諧函數模型計算效率高、易于解析表達，可根據各階次球諧系數計算出小行星非球形引力攝動，便于軌道設計.但在布里淵球域內發散，且存在截斷誤差.多面體模型可進行重力場的高精度建模，無發散區域，且能提供檢驗點相對小行星的位置信息，但因其是數值算法，不能給出受攝軌道解析式，計算量大，存在一定的局限性.實際使用的球諧系數大都是由軌道數據解算出來的［8］，而對于還未進行實際探測任務的小行星無法得到精確值.

針對小行星球諧系數的獲取問題，大多采用三軸橢球體模型法，計算由小行星簡化的橢球體的球諧系數，以此近似代替小行星的球諧系數［9］.可見小行星的形狀越不規則，誤差將越大.Scheeres［9］和胡維多等［10］以小行星二階二次重力場為出發點對小行星附近的動力學情況進行了研究，并利用小行星的三個主軸轉動慣量唯一確定球諧系數的二階項，但未給出高階項結果.張振江等［12］提出一種基于多面體模型的不規則形狀小行星重力場球諧系數求取方法，不僅提高了重力場建模的精度，同時可求其高階系數.但就計算過程中檢驗點的選取規則、多面體模型的劃分形式以及局部重力場建模等問題并未進行討論.較小行星重力場問題而言，地球的重力場研究已很成熟.眾多優秀學者［13-14］先后研究了適于描述地球重力場全球特性的全球模型和能夠精確構建局部地區重力場的局部模型.

鑒于此，本文在已有研究基礎上進一步研究了小行星重力場球諧系數的計算及優化問題.首先，分析說明了檢驗點的均勻分布形式和等面劃分的均勻多面體模型可以提高全球球諧系數的計算精度.隨后，為獲得高精度的局部重力場模型，提出局部球諧系數概念，并通過建立新的性能指標檢驗基于此組系數的球諧函數法在布里淵球域內的可行性，有效拓展了球諧函數法的應用范圍.最后，以小行星Eros433為實際算例，計算了全球及局部球諧系數，驗證本文方案的有效性和準確性.

1 重力場模型與球諧系數的求取

1.1 重力場模型

1)球諧函數模型.利用球諧函數描述的重力勢能函數為［5］

其中:G是萬有引力常數;M是小行星質量;R0是布里淵球半徑;是標準締合勒讓德多項式;是標準球諧系數;n、m是多項式階數、次數;R、φ、λ是檢驗點到中心天體質心的距離、緯度和經度.

將U對位置矢量求一階偏導數，可得重力加速度為

其中:

2)多面體模型.利用多面體模型描述的重力勢能函數為［7］

其中Eedges和Ffaces分別表示多面體模型的邊集合與面集合.

對應重力加速度函數為

其中:σ為天體密度;e表示邊;f表示面，向量是由檢驗點指向表面任一點的位置矢量;re是由檢驗點指向多面體邊上任一點的位置矢量.

1.2 球諧系數計算方法

對于未探測的星體，在小行星固連坐標系下，通過引入勢能值的虛擬觀測量，求解由球諧函數描述的勢能的拉普拉斯方程即可得到球諧系數.

任給一檢驗點Ri=(λi，φi，ri)，對于不規則形狀小行星利用多面體模型算得該檢驗點處的勢能作為勢能虛擬觀測量Ui，對于形狀規則的、勢能可解析表達的小行星，Ui用解析值代替，根據式(1)，則可構造一個以和為變量的方程［12］:

其中方程未知量的個數為l(l+2)(0無意義，為0，00為1).理論上只要在小行星勢能場中選取l(l+2)個不同的檢驗點建立線性方程組，構成Ax=b形式，就可以求解出共l(l+2)個的變量和但由于方程系數存在因子(R0/ri)nnm(sinφi)，使得檢驗點選取時無法確保方程線性無關.所以在實際計算過程中，通過增加方程個數組成超定方程，再求出此超定方程的最小二乘解來解決上述問題.

2 全球球諧系數

本文將基于小行星全球信息得到的球諧系數稱為全球球諧系數，即傳統意義上的球諧系數，其流程如圖1.因其滿足球諧函數的基函數具有整個空間域上的正交性，因此用它來描述的函數刻畫的都是全空間上的性質，適用于小行星布里淵球域外任意位置.利用精確的全球球諧系數表示的重力場全球模型可以比較直觀地反映出小行星某些全球性特征.

圖1 全球球諧系數研究流程

2.1 多面體模型面數的選取

鑒于勻質三軸橢球體重力場具有解析表達且存在高階球諧系數，以橢球體為例進行分析和驗證.

橢球體的球諧系數解析表達式如下［12］:

1)當n、m為所有正整數時，Snm=0;

2)當n或m為奇數時，Cnm=0;

3)當n、m為其他情況時

式中克羅內克符號

取橢球體半長軸a、b、c分別為16、8、6 km，密度為2.7 g/cm3，則布里淵球半徑R0為16 km.本文將星體表面用面積大小相當的三角形剖分得到的多面體模型稱為均勻多面體模型，上述橢球體多面體模型如圖2.分別取面數N為760、20 000和54 000，計算橢球體球諧系數，結果見表1.

圖2 勻質橢球體均勻多面體模型

表1 橢球體不同面數球諧系數計算結果

由表1結果可知，隨著多面體模型面數的增加，模型越逼近真實形狀，計算得到的球諧系數結果也越精確，但以犧牲計算效率為代價.

2.2 檢驗點選取方式

大量計算表明，檢驗點選取方式將直接影響球諧系數的計算精度.圖3給出了檢驗點的四種選取方式，固定相同的檢驗點個數(總數均為700個)和包絡體形狀(球體，且半徑為20 km)，圖中網格交點即為檢驗點:1)檢驗點均勻分布在球面上;2)檢驗點集中z軸兩端，對應橢球體南北兩極附近;3)檢驗點集中在x軸兩端，對應橢球體長軸兩端;4)檢驗點集中在橢球體北極上空.

采用上述檢驗點選取方式求解橢球體球諧系數，其中多面體模型面數N為54 000，所有選取方式對應的計算結果見表2.

圖3 檢驗點選取方式

表2 橢球體不同檢驗選取方式求解橢球體球諧系數

由表2計算結果易知，選取均勻分布的檢驗點(1)進行計算，可以得到相對其他方式更加精確的結果，各項相對誤差更小，全局取點方式優于局部取點方式.

2.3 多面體模型劃分形式

下面將橢球體表面由經緯度信息構成的四邊形網格結構分割成三角形，由此得到的多面體模型表面各三角形大小不一，這里簡稱為非均勻多面體模型，如圖4所示.為與均勻多面體模型進行比較分析，計算N=54 000的非均勻多面體模型球諧系數，結果見表3.

圖4 非均勻多面體模型

由表3可知，均勻多面體模型要優于非均勻多面體模型，前者20、22相對誤差僅為0.062 3%、0.020 7%，而 后 者 結 果 分 別 為 0.203 1% 和0.084 3%.非均勻多面體的三角形面積大小不均，位于南北兩極的三角形面積極小且相當密集，而在赤道附近，尤其長軸兩端處，三角形面積較大，分布稀疏，這就使得在相同面數條件下均勻模型更加精準的逼近小行星不規則形狀，進而得到更加精確的重力場模型.不過，理想橢球體僅為一特殊情況，對于實際的小行星而言，如何逼近其真實形狀的劃分方法并不絕對，如在平坦地區可以采用大三角塊，而在崎嶇地區建議細化分割.

表3 橢球體不同多面體模型球諧系數計算結果

3 局部球諧系數

針對探測器活動在某一特定空間范圍時，特別是在布里淵球域內，全局球諧系數往往不能滿足精度和效率的計算需求.因此本文借鑒地球重力場局部模型的概念，通過在特定活動區域集中選取檢驗點，求取局部球諧系數.引入局部球諧系數不旨在得到更加精確的全球重力場，而是以提高該局部區域重力場計算精度、增強軌道定軌和預報精度為目的.傳統的球諧函數法只適用于布里淵球域之外，而局部球諧系數模型可用于任意局部空間，不存在任何局限性.

為驗證局部球諧系數模型的性能優勢，在探測器活動區域選取軌道，抽取測試點，利用局部球諧系數模型求其勢能值，并與真實值進行分析比較.其研究流程如圖5所示.

為衡量局部球諧系數對局部活動區域勢能場的擬合程度，定義性能指標如下:

其中:N為軌道上測試點的抽取數目是局部球諧系數模型計算得到的測試點勢能值;U是測試點的真實勢能值，若小行星形狀規則，使用其解析解，若小行星形狀不規則，利用多面體模型結果近似代替真實值.PI越接近于1，說明結果擬合程度越高.

下面，仍以橢球體為例進行討論.

圖5 局部球諧系數研究流程

3.1 布里淵球域外局部球諧系數

以探測器活動區域集中在橢球體北極上空為例.因橢球體勢能場是解析的，此處，以其解析值［5］為式(5)左側的勢能虛擬觀測值，并利用圖3的(a)、(b)、(c)三種檢驗點方式，算得球諧系數見表4.為考察局部球諧系數性能，選取橢球體北緯85°，高度25 km處的盤旋軌道(見圖6)，分別抽取不同數目的測試點，計算PI，結果見表5.

圖6 北極上空盤旋軌道

表4 球諧系數計算結果

表5 橢球體北極盤旋軌道PI計算結果

由表 4、5結果知，檢驗點方式(1)、(2)到(4)，檢驗點逐步集中于橢球體北極，球諧系數與真實值的偏差增大，但對該活動區域的重力場建模精度顯著提高.隨著測試點的增多，局部球諧系數模型建模偏差積累效應明顯比全局球諧模型小的多.針對局部活動區域得到的局部球諧系數模型對該區域的重力場建模精度可以大大提高.

3.2 布里淵球域內局部球諧系數

假設探測器需要在圖7所示的近小行星區域進行低軌繞飛以獲得到高精度觀測數據，高度為10 km，處于布里淵球域內.

圖7 探測器局部活動區域

集中在圖7所示的局部活動區域內選取均勻分布形式的檢驗點，共計700個.求取球諧系數，見表6.選用圖8所示的繞飛軌道進行驗證.分別利用上述局部球諧系數和由檢驗點方式(1)得到的全球球諧系數對軌道抽取的測試點計算PI，結果見表7.由結果可知，局部球諧系數模型仍有極高的重力場擬合程度，而全球球諧系數已經發散，無法在該區域使用.

表6 局部球諧系數計算結果

圖8 探測器繞飛軌道

表7 橢球體統飛軌道PⅠ計算結果

4 小行星Eros球諧系數的求取

Eros433小行星于1989年發現，美國國家航空航天局(NASA)的NEAR-Shoemaker探測器對其進行了探測.據NASA公布，其質量(6.690 4±0.003)×1015kg，體積 2 503 ±25 km3，平均密度2.67 ±0.000 3 g·cm-3，大小34.4 ×11.2 ×11.2 km，均勻多面體模型如圖9所示［15］.

圖9 Eros多面體模型

4.1 全球球諧系數結果

利用本文方法求解其前15階全球球諧系數，結果見表8.將結果與傳統的三軸橢球體法和NEAR-Shoemaker實際環繞探測解算出的球諧系數相比較.計算過程中采用面數N分別為49 152和196 608的兩種Eros均勻多面體模型.

表8 Eros球皆系數求取結果

4.2 局部球諧系數結果

以探測器活動在高度為25 km的Eros北極上空為例，利用圖3檢驗點方式(1)、(2)、(4)，基于Eros多面體模型求解球諧系數，并選取圖6所示的軌道進行驗證，計算PI，結果見表9.可見，對不規則形狀小行星，由探測器局部活動區域得到的局部球諧系數相比全球球諧系數對該區域的重力場擬合程度更高.

以探測器活動在圖7所示的局部區域為例，集中在此區域選取檢驗點，基于Eros多面體模型計算局部球諧系數，利用圖3檢驗點方式(1)求取其全球球諧系數，并選擇圖8所示的繞飛軌道進行驗證，計算PI，結果見表10.可見，基于局部球諧系數的球諧函數法在布里淵球域內針對該局部活動區域仍能保證較高的重力場建模精度.

表9 Eros北極盤旋軌道PI計算結果

表10 Eros繞飛軌道PI計算結果

5 結語

本文討論了小行星球重力場球諧系數計算相關問題.通過定性討論和定量分析，結果表明均勻多面體模型和均勻分布的檢驗點保證了高精度的全球球諧系數的獲得.針對探測器在某一特定區域活動情況，提出局部球諧系數概念，提高了局部區域重力場建模精度，擴展了球諧函數法的應用范圍.

重力場模型提供分析、描述和設計地球表面及其外空間一切物體力學行為的先驗重力場約束.精細的重力場信息為小行星探測任務的成功提供保障.而到目前為止，小行星重力場描述雖然眾多，但都存在一定缺陷.因此，尋找一種以動態觀、整體論的方法描述的具有普適意義重力場模型，更加直觀地解釋分析小行星重力場特性，具有一定的必要性和緊迫性，也將成為未來小行星探測的一大挑戰.

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

［2］ BROSCHART S B，SCHEERES D J.Control of hovering spacecraft near small bodies:application to asteroid 25143 Itokawa ［J］.Journal of Guidance，Control，and Dynamics，2005，28(2):343-354.

［3］YANO H，KUBOTA T，MIYAMOTO H，et al.Touchdown of the hayabusa spacecraft at the Muses Sea on Itokawa［J］.Science，2006，312(5778):1350-1353.

［4］徐偉彪，趙海斌.小行星深空探測的科學意義和展望［J］.地球科學進展，2005，20(11):1183-1190.

［5］CASTOTTO S，MUSOTTO S.Methods for computing the potential of an irregular，homogeneous，solid body and its gradient［J］.AIAA Paper，2000，4023:82-96.

［6］KAULA W M.Theory of Satellite Geodesy:Applications ofSatellitesto Geodesy［M］. Mineola:Dover Publications，2000:4-8.

［7］WERNER R A，SCHEERES D J.Exterior gravitation of a polyhedron derived and compared with harmonic and mascon gravitation representationsofasteroid 4769 castalia［J］. CelestialMechanicsand Dynamical Astronomy，1997，65(3):313-344.

［8］ PARK R S， WERNER R A， BHASKARAN S.Estimating small-body gravity field from shape model and navigation data［J］.Journal of Guidance，Control，and Dynamics，2010，33(1):212-221.

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

［10］HU W，SCHEERES D J.Numerical determination of stability regions for orbital motion in uniformly rotating second degree and order gravity fields［J］.Planetary and Space Science，2004，52(8):685-692.

［11］SCHEERES D J.Orbit mechanics about asteroids and comets［J］. JournalofGuidance Controland Dynamics，2012，35(3):987.

［12］張振江，崔祜濤，任高峰.不規則形狀小行星引力環境建模及球諧系數求取方法［J］.航天器環境工程，2010，27(3):383-388.

［13］寧津生.地球重力場模型及其應用［J］.冶金測繪，1994，3(2):2-6.

［14］LI J C，CHAO D B，NING J S.Spherical cap harmonic expansion for local gravity-field representation ［J］.Manuscripta geodaetica，1995，20(4):265-277.

［15］National Aeronautics and Space Administration.http://sbn.psi.edu/pds/resource/rshape.html.2010-4-20.