束磊，劉睿，張迎春

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

引言

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

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

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

1 兩種引力場(chǎng)模型

1.1 球諧模型

小行星引力位的球諧展開(kāi)式為:

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

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

1.2 多面體模型

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

其中:

U對(duì)各坐標(biāo)軸求導(dǎo)，可求得引力位梯度，其分量即為各坐標(biāo)軸上的引力加速度。

以上各式中的符號(hào)定義如表1所示。

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

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

圖1 愛(ài)神星多面體三維圖

圖2 愛(ài)神星多面體平面圖

2 軌道動(dòng)力學(xué)方程

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

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

圖3 星體固定坐標(biāo)系和慣性坐標(biāo)系

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

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

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

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

則:

將各坐標(biāo)代入上式得:

又因?yàn)??U/?r，所以，探測(cè)器運(yùn)行的軌道動(dòng)力學(xué)方程為:

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

3 引力場(chǎng)模型仿真對(duì)比

3.1 引力加速度對(duì)比

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

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

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

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

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

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

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

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

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

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

對(duì)方程兩邊同時(shí)求導(dǎo)，得:

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

表2 探測(cè)器在軌運(yùn)行初始條件

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

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

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

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

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

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

3.3 分界半徑確定

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

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

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

4 結(jié)論

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

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

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

（3）考慮到球諧模型計(jì)算簡(jiǎn)單快捷但存在發(fā)散問(wèn)題，而多面體模型無(wú)發(fā)散問(wèn)題但計(jì)算量大，故結(jié)合兩者的優(yōu)點(diǎn)，確定了xy平面內(nèi)的探測(cè)器軌道的分界半徑取20 km。對(duì)小行星引力場(chǎng)建模時(shí)，分界半徑外空間用球諧模型，分界半徑內(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.