不規(guī)則形狀小行星引力場組合建模方案

束磊，劉睿，張迎春

（1.哈爾濱工業(yè)大學衛(wèi)星技術(shù)研究所，黑龍江哈爾濱 150001；2.西北工業(yè)大學精確制導與控制研究所，陜西西安 710072）

引言

1996年，為探測近地 Eros 433小行星，美國NASA的NEAR項目發(fā)射了名為“鞋匠號”（Shoemaker）的探測器。2001年，經(jīng)過一年的在軌觀測，“鞋匠號”成功地在愛神星表面軟著陸。之后，歐洲和日本也各自發(fā)射了多個小行星探測器。

小行星的引力場模型對于探測器的動力學建模至關(guān)重要。由于多數(shù)小行星的形狀極不規(guī)則，其引力場與地球、月球這樣的近圓天體有很大差別。天體重力場分析通常采用球諧模型方法［1-4］。但是，球諧模型中假定了一個作為參考的限制球［2］。在限制球外空間，有限階的球諧展開式能夠較為準確地反映小行星引力場的強度和分布；但在限制球內(nèi)，球諧展開式會發(fā)散，無法準確描述真實引力場。為了解決這一問題，本文引入多面體模型法［5-6］，盡管該模型存在因為真實小行星密度非恒定引起的誤差，但由于該模型在限制球內(nèi)不存在發(fā)散問題，故在限制球內(nèi)其精度相對于球諧模型更高。文獻［7］提到了多面體模型，但也只是將其應(yīng)用于求解球諧模型的球諧系數(shù)，還沒有用多面體模型對不規(guī)則形狀小行星的引力場建模的報道。文獻［3］中引力場模型用的是球諧模型，但沒有考慮球諧模型的發(fā)散問題。

本文結(jié)合球諧模型和多面體模型各自的優(yōu)缺點，提出了一種基于軟著陸或超低軌道飛行過程中軌道高度不斷變化對引力場模型進行切換的思想，即在分界半徑外的引力場用球諧模型計算，分界半徑內(nèi)用多面體模型計算。

1 兩種引力場模型

1.1 球諧模型

小行星引力位的球諧展開式為:

球諧模型中，限制球定義為:球心星體為質(zhì)心；半徑為參考半徑Re，通常Re取星體的最大半徑。確定各階球諧系數(shù)后，可得到該引力位的解析式:

NASA官網(wǎng)上公布有Eros 433小行星的8階球諧系數(shù)，這些系數(shù)是根據(jù)探測器得到的數(shù)據(jù)計算而來的，本文利用這些系數(shù)建立球諧模型。但對于不規(guī)則形狀小行星，當R＜Re時，即檢驗點在限制球內(nèi)時，式（1）存在（Re/R）n，模型存在發(fā)散的可能。

1.2 多面體模型

多面體模型是由若干個小三角形組成的多面體來逼近真實小行星的形狀。在假設(shè)小行星內(nèi)部等密度情況下，可通過三重積分近似計算出整個小行星的外引力場強度和分布［6］。由于無法準確得到真實小行星的密度分布，通常將小行星的密度近似為其平均密度［6］。真實小行星密度非恒定會引起一定的誤差，但在限制球內(nèi)其精度相對于發(fā)散的球諧模型更高。將小行星等效為內(nèi)部密度恒定的多面體后，根據(jù)文獻［6］中的推導，其引力位為:

其中:

U對各坐標軸求導，可求得引力位梯度，其分量即為各坐標軸上的引力加速度。

以上各式中的符號定義如表1所示。

表1 多面體模型中相應(yīng)符號的定義

以愛神星Eros 433為例，參考NASA官網(wǎng)公布的探測數(shù)據(jù)，取其表面3 897個點，構(gòu)成7 790個三角形面，這樣構(gòu)成的多面體能夠高度模擬愛神星的真實形狀，如圖1和圖2所示。下文在這一模型的基礎(chǔ)上進行仿真計算。

圖1 愛神星多面體三維圖

圖2 愛神星多面體平面圖

2 軌道動力學方程

以小行星軟著陸為應(yīng)用目的，本文將探測器繞小行星運動的軌道動力學方程建立在小行星固定坐標系下。首先建立小行星的固定坐標系Oxyz:原點O在小行星質(zhì)心，z軸為小行星自轉(zhuǎn)軸（最大慣量軸），x軸為最小慣量軸，y軸滿足右手定則。

建立慣性坐標系Ox1y1z1:原點O在小行星質(zhì)心，x1軸平行于小行星赤道面與黃道面的交線，且指向春分點，z1軸為自轉(zhuǎn)軸，y1軸滿足右手定則。星體固定坐標系Oxyz和慣性坐標系Ox1y1z1之間的關(guān)系如圖3所示。

圖3 星體固定坐標系和慣性坐標系

設(shè)［xyz］T及［x1y1z1］T分別為物體在固定坐標系和慣性坐標系下的坐標。由坐標變換可得:

探測器在固定坐標系中的位置、速度和加速度矢量分別為r，r'和r″。

由于小行星自身旋轉(zhuǎn)，故星體固定坐標系為動坐標系，其相對于慣性坐標系的旋轉(zhuǎn)角速度為Ω。假設(shè)小行星只繞最大慣量軸旋轉(zhuǎn)，則其旋轉(zhuǎn)角速度在慣性空間的坐標可表示為Ω=［0 0ω］T。由力學知識可得:

式中，aa為絕對加速度；ae為牽連加速度；ar為相對加速度；ac為科氏加速度。其中:

則:

將各坐標代入上式得:

又因為=?U/?r，所以，探測器運行的軌道動力學方程為:

愛神星Eros 433繞其最大慣量軸自轉(zhuǎn)，自轉(zhuǎn)角速度約為每天1 639.4°，自轉(zhuǎn)不能忽略。

3 引力場模型仿真對比

3.1 引力加速度對比

以愛神星Eros為例，畫出兩種模型在不同軌道高度下引力加速度的分布曲線，結(jié)果如圖4～圖6所示。

圖4 xy平面初始高度50 km軌道處三種引力加速度對比

從圖4可以看出，在高軌道處（初始高度50 km），x，y軸上繞軌一周的加速度變化曲線基本一致，量級在10－7，z軸上有一定的差異，量級為10－10，可以忽略。

圖5 yz平面初始高度10 km軌道處三種引力加速度對比

圖6 yz平面初始高度10 km軌道處多面體模型與中心引力模型加速度對比

圖5表明，在低軌道處（初始高度10 km），應(yīng)用球諧模型得到的引力加速度變化劇烈，且與多面體模型和中心引力模型加速度差異很大，約為其105倍，說明球諧模型在此軌道處發(fā)散。

圖6為yz平面內(nèi)初始高度10 km軌道處多面體模型和中心引力模型的引力加速度繞軌一周的變化曲線。其變化趨勢基本一致，說明了多面體模型沒有發(fā)散問題，應(yīng)用于低軌是合理的。

3.2 非球形引力作用下探測器軌道變化仿真

設(shè)定探測器初始軌道位置及飛行速度，在兩種模型的引力場作用下，比較軌道高度的變化。

將慣性坐標系轉(zhuǎn)換為固定坐標系可得:

對方程兩邊同時求導，得:

假設(shè)探測器初始時刻運行在xy平面或yz平面內(nèi)，設(shè)定慣性坐標系下初始條件如表2所示。

表2 探測器在軌運行初始條件

結(jié)合軌道動力學方程，利用Simulink搭建仿真模塊，仿真得到不受控飛行軌道在兩種模型下隨時間變化的曲線，如圖7和圖8所示。

圖7為xy平面內(nèi)50 km軌道的變化曲線，仿真時間為105s。該軌道高度兩種模型的x和y坐標基本重合，而z坐標呈周期性變化。軌道半徑同樣呈周期性變化，且兩者變化趨勢一致。球諧模型軌道變化約為2.2 km，多面體模型的軌道變化約為1.8 km，差別不大。

圖7 xy平面初始高度50 km軌道變化

圖8為yz平面內(nèi)20 km軌道的變化曲線。兩種模型坐標初始基本重合，球諧模型的軌道變化很劇烈，當軌道半徑降為16 km左右后，呈明顯的發(fā)散趨勢，而此時多面體模型的軌道高度仍然在攝動的作用下緩慢地變化著。

圖8 yz平面初始高度20 km軌道變化

由此可見，球諧模型并不適用于小行星附近空間（限制球內(nèi)），而多面體模型在限制球內(nèi)還是限制球外都適用。多面體模型計算量較大，利用Matlab仿真耗時較多，但利用C語言等底層語言編程可以顯著提高速度，可以實時運算［8］，而球諧模型計算更加簡便快捷，故針對軟著陸或超低空飛行，考慮將兩種模型結(jié)合起來使用，一般情況下，分界半徑比限制球半徑稍大。

3.3 分界半徑確定

探測愛神星的探測器軌道一般是xy平面內(nèi)的逆行軌道，本文著重分析兩種重力場模型的偏差，通過計算兩種模型加速度，可得出兩種模型在平面內(nèi)各個點處的加速度偏差。加速度偏差分布見圖9。圖中，不同的灰度代表不同的偏差值。

圖9 球諧引力場模型與多面體引力場模型加速度偏差分布圖

由圖9可知，在20 km軌道半徑以內(nèi)，越接近星表，兩種加速度模型的偏差越大，故將分界半徑取為20 km。在分界半徑外，球諧模型能快速準確地計算引力場；在分界半徑內(nèi)，球諧模型是發(fā)散的，而多面體模型仍能準確描述小行星外引力場。

4 結(jié)論

本文提出了一種小行星引力場組合建模方案，以愛神星引力場模型為例，通過仿真可得以下結(jié)論:

（1）在高軌道處（初始高度50 km），多面體模型和球諧模型之間引力加速度誤差基本可以忽略，且對不受控探測器軌道攝動影響的差別不大，說明在高軌道處兩種模型都能較準確地描述引力場。

（2）在低軌道處（初始高度10 km），球諧模型引力加速度約為多面體模型的105倍，球諧模型在此軌道處發(fā)散，說明球諧模型不適用于限制球內(nèi)低軌衛(wèi)星軌道。

（3）考慮到球諧模型計算簡單快捷但存在發(fā)散問題，而多面體模型無發(fā)散問題但計算量大，故結(jié)合兩者的優(yōu)點，確定了xy平面內(nèi)的探測器軌道的分界半徑取20 km。對小行星引力場建模時，分界半徑外空間用球諧模型，分界半徑內(nèi)空間用多面體模型。

［1］Werner R A.Spherical harmonic coefficients for the potential of a constant-density polyhedron ［J］.Computers and Geosciences，1997，23（10）:1071-1077.

［2］Scheeres D J，Williams B G，Miller J K.Evaluation of the dynamic environment of an asteroid:Applications to 433 Eros［J］.Journal of Guidance，Control，and Dynamics，2000，23（3）:466-475.

［3］Huang X Y，Cui H T，Cui P Y.An autonomous optical navigation and guidance for soft landing on asteroids［J］.Acta Astronautica，2004，54（10）:763-771.

［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（1）:3-17.

［5］Cangahuala L A.Augmentations to the polyhedral gravity model to facilitate small body navigation［R］.AAS-05-4185，2005.

［6］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-344.

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

［8］Gregory Lantoine.Optimal trajectories for soft landing on asteroids［R］.AE8900 MS Special Problems Report，Space Systems Design Lab，2006.