小天體引力場(chǎng)建模技術(shù)進(jìn)展

尚海濱，韋炳威，盧櫸承

（北京理工大學(xué) 宇航學(xué)院，北京 100081）

引言

太陽(yáng)系小天體通常是指太陽(yáng)系中除行星、矮行星以及天然衛(wèi)星以外的天體，它們一般包括小行星、彗星和流星體等。近年來(lái)，太陽(yáng)系小天體已經(jīng)成為人類深空探測(cè)活動(dòng)的重要目標(biāo)。對(duì)小天體實(shí)施探測(cè)主要有3方面的動(dòng)機(jī)：探測(cè)太陽(yáng)系小天體對(duì)揭開太陽(yáng)系形成與演化之謎有重要意義[1-3]；小天體探測(cè)有助于應(yīng)對(duì)近地小天體對(duì)地球的碰撞威脅，保障人類社會(huì)可持續(xù)發(fā)展[4-5]；通過(guò)光譜觀測(cè)發(fā)現(xiàn)一些小天體蘊(yùn)含著豐富資源[6-7]，對(duì)這些資源的開發(fā)與利用具有潛在經(jīng)濟(jì)價(jià)值，同時(shí)可為太空工業(yè)的發(fā)展提供重要支持。

早期小天體的探測(cè)方式主要以飛越、環(huán)繞為主，美國(guó)國(guó)家航空航天局（National Aeronautics and Space Administration，NASA）發(fā)射的“伽利略號(hào)”（Galileo）探測(cè)器依次飛越了小行星Gaspra 951 和Ida 243 及其衛(wèi)星Dactyl[8]；“深空一號(hào)”（Deep Space 1）探測(cè)器飛越了小行星布雷爾（Braille 9969）和彗星布洛利（Borrelly 19P ）[9]；“黎明號(hào)”（Dawn）探測(cè)器對(duì)灶神星（Vesta）進(jìn)行了長(zhǎng)達(dá)1 a的環(huán)繞探測(cè)[10]；“嫦娥二號(hào)”（Chang’E-2，CE-2）探測(cè)器在其拓展任務(wù)中成功飛越了小行星圖塔蒂斯（Toutatis 4179）[11-13]。這些任務(wù)獲取了大量的小天體探測(cè)數(shù)據(jù)，加深了人類對(duì)太陽(yáng)系小天體的認(rèn)識(shí)。

為獲得更為豐富的科學(xué)數(shù)據(jù)，例如巖石特征、礦物成分、地表熱特性和磁特性等，需對(duì)小天體進(jìn)行近距離和表面的深入探測(cè)。為此，各航天機(jī)構(gòu)相繼開展了附著、表面巡視及采樣返回等復(fù)雜探測(cè)任務(wù)。日本宇宙航空研究開發(fā)機(jī)構(gòu)（Japan Aerospace eXploration Agency，JAXA）開展了“隼鳥一號(hào)”（Hayabusa 1）任務(wù)[14]和“隼鳥二號(hào)”（Hayabusa 2）任務(wù)[15]、歐洲航天局（European Space Agency，ESA）開展了“羅塞塔”（Rosetta）任務(wù)[16]、NASA開展了“舒梅克探測(cè)器”（NEAR-Shoemaker）任務(wù)[17]和“歐西里斯號(hào)”（OSIRIS-REx）任務(wù)[18]。對(duì)于近距離探測(cè)任務(wù)而言，小天體引力場(chǎng)的作用變得不可忽視，特別對(duì)于無(wú)軌控能力的巡視器和著陸器，復(fù)雜引力場(chǎng)導(dǎo)致了更高的任務(wù)風(fēng)險(xiǎn)。Hayabusa探測(cè)器于2005年11月12日釋放了著陸器“智慧女神”（Minerva），試圖著陸于小行星糸川（Itokawa 25143 ）表面，由于釋放位置高度過(guò)高導(dǎo)致Minerva最終逃逸了Itokawa[19]。近年來(lái)，各個(gè)國(guó)家也在積極開展雙小行星系統(tǒng)探測(cè)的國(guó)際合作，由NASA和ESA合作開展的小行星防御與偏轉(zhuǎn)評(píng)估任務(wù)（Asteroid Impact &Deflection Assessment，AIDA）已經(jīng)完成了對(duì)目標(biāo)雙小行星系統(tǒng)Didymos 65803的初步探測(cè)[20]，并計(jì)劃由ESA完成對(duì)該雙小行星系統(tǒng)的近距離詳細(xì)勘測(cè)。兩個(gè)小行星復(fù)雜的引力場(chǎng)和雙小行星系統(tǒng)獨(dú)特的構(gòu)型，為雙小行星系統(tǒng)附近的引力場(chǎng)計(jì)算帶來(lái)了新的挑戰(zhàn)。

精確構(gòu)建小天體附近及表面引力場(chǎng)模型是發(fā)展近距離探測(cè)技術(shù)，降低任務(wù)風(fēng)險(xiǎn)的重要基礎(chǔ)。若能實(shí)現(xiàn)小天體引力場(chǎng)的高精度建模，將有助于實(shí)現(xiàn)探測(cè)器、巡視器和著陸器的任務(wù)可靠性設(shè)計(jì)，包括環(huán)繞軌道、懸停軌道、轉(zhuǎn)移軌道、受控附著和彈道附著軌道的設(shè)計(jì)和穩(wěn)定性分析，以及表面巡視軌跡的魯棒規(guī)劃等。除了實(shí)際任務(wù)需求，小天體精確引力場(chǎng)建模還是一些動(dòng)力學(xué)和行星科學(xué)研究的基礎(chǔ)，例如周期軌道的搜索、表面濺射物的演化、雙小行星系統(tǒng)的演化和小天體內(nèi)部物質(zhì)分布重構(gòu)等。

由于小天體具有不規(guī)則形狀和非均勻內(nèi)部質(zhì)量分布的重要特點(diǎn)，其引力場(chǎng)的引力勢(shì)、引力加速度和引力梯度通常不具備解析閉合的表達(dá)形式。因此，從數(shù)學(xué)上描述小天體附近的復(fù)雜引力場(chǎng)具有很大的挑戰(zhàn)性，多年來(lái)許多學(xué)者們致力于解決這一難題，獲得了諸多重要成果。

本文梳理了兩百多年來(lái)學(xué)者們對(duì)引力場(chǎng)建模問(wèn)題的認(rèn)識(shí)與發(fā)展，分別詳細(xì)闡述了球諧函數(shù)展開法、橢球諧函數(shù)展開法、內(nèi)外引力場(chǎng)法、球Bessel函數(shù)展開法、質(zhì)點(diǎn)群法、多面體法、虛擬邊界元法以及基于數(shù)據(jù)驅(qū)動(dòng)的快速引力場(chǎng)計(jì)算方法的基本原理，評(píng)述了不同方法的優(yōu)點(diǎn)與不足。除引力場(chǎng)表征問(wèn)題，在此基礎(chǔ)上還回顧了雙小行星動(dòng)力學(xué)系統(tǒng)建模技術(shù)的發(fā)展，按照從簡(jiǎn)單到復(fù)雜建模歷程分別介紹了球形限制假設(shè)模型、級(jí)數(shù)展開模型和有限元模型，可為雙小行星系統(tǒng)附近的引力場(chǎng)建模技術(shù)研究提供參考。

1 單小天體引力場(chǎng)建模技術(shù)

對(duì)于小天體近距離探測(cè)任務(wù)而言，探測(cè)器在探測(cè)目標(biāo)附近的軌道運(yùn)動(dòng)主要受到小天體引力作用的影響，精確構(gòu)建小天體附近的引力場(chǎng)有助于開展小天體探測(cè)任務(wù)設(shè)計(jì)。同時(shí)，小天體的引力場(chǎng)建模還是天體動(dòng)力學(xué)分析的基礎(chǔ)，能夠支持行星科學(xué)的研究，具有較高的科學(xué)價(jià)值。

位于小天體外部空間一點(diǎn)r的探測(cè)器所感受到的小天體引力勢(shì)定義為

其中：Ω為小天體體積區(qū)域；s為Ω內(nèi)一點(diǎn)；G為引力常數(shù)；ρ 為小天體內(nèi)部的密度分布函數(shù)。

探測(cè)器感受到的引力加速度則可表示為?U，引力梯度張量可表示為??U，具體表達(dá)式為

1.1 級(jí)數(shù)展開法

1.1.1 球諧函數(shù)展開法

無(wú)窮級(jí)數(shù)展開法是利用引力勢(shì)函數(shù)的解析性質(zhì)，將原求解不規(guī)則外形非均質(zhì)分布小天體的引力勢(shì)問(wèn)題轉(zhuǎn)化為計(jì)算有限項(xiàng)基函數(shù)的線性組合問(wèn)題。球諧函數(shù)展開法是最為經(jīng)典的引力場(chǎng)表征方法。球諧函數(shù)構(gòu)成一組完全正交函數(shù)基，在球面上定義的連續(xù)函數(shù)都可寫成球諧函數(shù)的線性組合。類似于定義在圓周上的周期函數(shù)，可通過(guò)Fourier級(jí)數(shù)表示為三角函數(shù)的線性組合。最初法國(guó)數(shù)學(xué)家Legendre于1782年在研究質(zhì)點(diǎn)系引力場(chǎng)時(shí)發(fā)現(xiàn)觀察點(diǎn)到小天體內(nèi)微元距離的倒數(shù)能以級(jí)數(shù)形式表示為

其中：r=‖r‖；s=‖s‖；Pn(·)為n階Legendre多項(xiàng)式；θ為矢量r與矢量s之間的夾角。

William 等[21]將這一工作進(jìn)一步完善，在球坐標(biāo)下任意天體的引力勢(shì)可表示為

其中：M為天體質(zhì)量；R為參考半徑，通常選取天體最小包圍球（稱為Brillouin球）的半徑；Cnm和Snm稱為Stokes系數(shù)或球諧系數(shù)。

級(jí)數(shù)式（5）中締合Legendre多項(xiàng)式與三角函數(shù)的乘積體現(xiàn)了球面函數(shù)的部分，稱為球諧函數(shù)。級(jí)數(shù)式（4）形式僅當(dāng)s

圖1 Brillouin球與非球形小天體示意圖Fig.1 Illustration of Brillouin sphere and non-spherical small body

1.1.2 橢球諧函數(shù)展開法

球諧函數(shù)展開法對(duì)于非球形天體引力場(chǎng)描述有局限性，為擴(kuò)大收斂區(qū)域，Hobson等[22]總結(jié)出了橢球諧函數(shù)展開法。該方法首先引入橢球坐標(biāo)系(λ1,λ2,λ3)描述空間一點(diǎn)。這3個(gè)坐標(biāo)參數(shù)分別描述了3個(gè)共焦二次曲面：橢球曲面、單葉雙曲面、雙葉雙曲面。3個(gè)共焦曲面在每個(gè)笛卡爾坐標(biāo)象限中決定一個(gè)空間點(diǎn)。

已知引力勢(shì)函數(shù)滿足Laplace方程，即?2U=0，則在橢球坐標(biāo)下描述關(guān)于U的Laplace方程，再根據(jù)分離變量假設(shè)方程的解可寫為U=E1(λ1)E2(λ2)E3(λ3)，代入橢球坐標(biāo)系下的Laplace方程可得到對(duì)應(yīng)的Lamé方程。Lamé方程的解稱為L(zhǎng)amé函數(shù)或橢球諧函數(shù)，最終引力勢(shì)函數(shù)可以表示為橢球諧函數(shù)的級(jí)數(shù)形式，即

其中：αnm為待定的橢球諧系數(shù)為參考橢球?qū)?yīng)的坐標(biāo)值；Enm(·)為第一類Lamé函數(shù)；Fnm(·)為第二類Lamé函數(shù)。

級(jí)數(shù)式（6）的收斂域?yàn)閰⒖紮E球外部，即λ1≥對(duì)于細(xì)長(zhǎng)型小天體而言，可將參考橢球定義為包圍小天體的最小橢球（即Brillouin橢球），橢球諧函數(shù)展開法能夠描述更多靠近小天體表面的空間區(qū)域（如圖 2所示）。

圖2 Brillouin橢球與非球形小天體示意Fig.2 Illustration of Brillouin ellipsoid and non-spherical small body

Walter[23]就引力場(chǎng)建模中橢球諧系數(shù)和球諧系數(shù)相互轉(zhuǎn)化的問(wèn)題開展了研究，但其轉(zhuǎn)換過(guò)程涉及大量橢圓積分的復(fù)雜數(shù)值計(jì)算，使得所提方法的應(yīng)用受到了限制。Dechambre等[24]基于Dirichlet問(wèn)題發(fā)展了另一套球諧系數(shù)到橢球諧系數(shù)的轉(zhuǎn)換方法，簡(jiǎn)化了橢球諧系數(shù)的計(jì)算，之后學(xué)者致力于發(fā)展精確系數(shù)轉(zhuǎn)換方法[25-27]。盡管橢球諧函數(shù)展開法具有更好的收斂特性，但其應(yīng)用遠(yuǎn)不如球諧函數(shù)展開法，這是因?yàn)橛?jì)算各項(xiàng)橢球諧函數(shù)相比各階球諧函數(shù)具有更高的計(jì)算復(fù)雜度，甚至Lamé函數(shù)的遞歸公式還沒(méi)被發(fā)現(xiàn)，而締合Legendre多項(xiàng)式具有成熟的遞歸公式體系，能夠便于計(jì)算機(jī)的求解[28]。由于第一類和第二類Lamé函數(shù)存在零點(diǎn)和奇異點(diǎn)（當(dāng)λi=h或λi=k），橢球諧級(jí)數(shù)在階數(shù)較大時(shí)候容易出現(xiàn)數(shù)值不穩(wěn)定的情況，導(dǎo)致計(jì)算精度退化[28]。

1.1.3 內(nèi)外引力場(chǎng)法

為將級(jí)數(shù)展開法的收斂域擴(kuò)大到小天體表面任意一點(diǎn)，Takahashi等[29]提出了內(nèi)外球諧函數(shù)展開法?；舅枷胧且訠rillouin球外部空間任意一點(diǎn)建立坐標(biāo)系，并將式（4）改寫為

在式（7）中，當(dāng)r

其中：Rin為內(nèi)Brillouin球半徑；為內(nèi)球諧系數(shù)，如圖3所示，通過(guò)內(nèi)球諧函數(shù)展開式可以將級(jí)數(shù)收斂域擴(kuò)展到小天體表面一點(diǎn)。

圖3 內(nèi)外Brillouin球與非球形小天體示意Fig.3 Illustration of internal and external Brillouin sphere and non-spherical small body

雖然內(nèi)外引力場(chǎng)法能通過(guò)內(nèi)球諧函數(shù)展開式來(lái)描述原Brillouin球內(nèi)的部分區(qū)域引力場(chǎng)，但該方法仍無(wú)法使用統(tǒng)一表達(dá)式實(shí)現(xiàn)引力場(chǎng)的全局表征。內(nèi)Brillouin球僅能延伸至小天體表面一點(diǎn)，因此該方法可應(yīng)用于探測(cè)器在小天體表面一點(diǎn)的附著模擬。若要描述小天體表面其他區(qū)域引力場(chǎng)，則還需建立更多內(nèi)Brillouin球來(lái)覆蓋更多小天體表面附近區(qū)域。

1.1.4 球Bessel函數(shù)展開法

為得到Brillouin球內(nèi)部統(tǒng)一的引力場(chǎng)級(jí)數(shù)表征形式，Herrera-Sucarrat等[30-31]提出了球Bessel函數(shù)展開法?；舅枷胧菍rillouin球內(nèi)視為一個(gè)具有一定質(zhì)量分布的實(shí)物體區(qū)域，則引力勢(shì)函數(shù)在該區(qū)域內(nèi)滿足Poisson方程，即

為求解式（9）方程，可假設(shè)Brillouin球內(nèi)的質(zhì)量被重新分布以滿足下述Helmholtz方程。

其中：k為常數(shù)。

采用分離變量假設(shè)引力勢(shì)函數(shù)具有的形式為：U(r,φ,λ)=(r)Φ(φ)Λ(λ)，其中徑向勢(shì)函數(shù)的波浪字符是為了與Brillouin球半徑R做區(qū)分。將該形式代入Helmholtz方程并表示為Bessel方程（對(duì)應(yīng)不同本征值k=km的本征方程）。

式（11）的解為第一類球Bessel函數(shù)jn(klnr)，于是Poisson方程（9）的解可表示為

其中：αln=klnR為歸一化的本征值。

若將Brillouin球外部的引力勢(shì)函數(shù)也寫成球Bessel函數(shù)和球諧函數(shù)的形式可得

球Bessel函數(shù)展開法的難點(diǎn)是系數(shù)的求解。由于系數(shù)Alnm和Blnm的求解精度限制于原球諧系數(shù)精度以及人為引入的約束，通常只能獲得低階系數(shù)。

1.2 質(zhì)點(diǎn)群法

質(zhì)點(diǎn)群法的主要原理是利用有限個(gè)質(zhì)點(diǎn)或均質(zhì)小球填充小天體內(nèi)部空間，則小天體外部的引力勢(shì)可以表示為所有質(zhì)點(diǎn)產(chǎn)生的引力勢(shì)的疊加，即

該方法避免了復(fù)雜的體積積分計(jì)算，計(jì)算效率較高且建模方式簡(jiǎn)單易實(shí)現(xiàn)，廣泛應(yīng)用于小天體演化和附近軌道動(dòng)力學(xué)的研究中，Richardson等[32-33]在N體問(wèn)題開源求解器pkdgrav的框架下對(duì)小行星Itokawa的形成機(jī)理進(jìn)行了研究；Geissler等[34]使用質(zhì)點(diǎn)群模型研究了小行星Ida及其衛(wèi)星Dactyl表面濺射物的再沉積過(guò)程；Shang等[35]以質(zhì)點(diǎn)群法作為初始引力場(chǎng)模型，提出了小行星附近全局周期軌道的快速搜索方法，并研究了小行星Eros附近周期軌道的拓?fù)漕愋汀?/p>

質(zhì)點(diǎn)群法相比無(wú)窮級(jí)數(shù)展開法的主要優(yōu)勢(shì)：①能夠較好擬合小天體的不規(guī)則形狀；②通過(guò)給每個(gè)質(zhì)點(diǎn)設(shè)置不同的質(zhì)量可以表征小天體內(nèi)部質(zhì)量非均勻分布的特點(diǎn)；③不存在收斂域的限制，理論上屬于小天體引力場(chǎng)的全局表征方法。然而，質(zhì)點(diǎn)群法的缺點(diǎn)也較為明顯，主要表現(xiàn)：①模型精度隨質(zhì)點(diǎn)個(gè)數(shù)增加提升速度較慢，導(dǎo)致當(dāng)精度需求較高時(shí)，所需質(zhì)點(diǎn)群總數(shù)將劇烈增加，計(jì)算效率顯著降低；②雖然質(zhì)點(diǎn)群能夠近似表征小天體內(nèi)部質(zhì)量非均勻分布特征，但當(dāng)質(zhì)點(diǎn)數(shù)較大時(shí)會(huì)引起協(xié)方差矩陣奇異，導(dǎo)致各個(gè)質(zhì)點(diǎn)的質(zhì)量難以準(zhǔn)確估計(jì)，Park等[36]證明了當(dāng)質(zhì)點(diǎn)數(shù)增大，通過(guò)軌道數(shù)據(jù)估計(jì)的各質(zhì)點(diǎn)質(zhì)量的不確定度迅速增加；③當(dāng)觀察點(diǎn)靠近小天體表面，引力場(chǎng)的精度會(huì)迅速退化，這是由觀察點(diǎn)靠近質(zhì)點(diǎn)群時(shí)所引起的近奇異現(xiàn)象所致，因此該方法并不適用于小天體表面的動(dòng)力學(xué)研究。

1.3 多面體法

1.3.1 均質(zhì)多面體法

Werner等[37]在前人研究的基礎(chǔ)上首先完整推導(dǎo)了任意均質(zhì)多面體引力勢(shì)、引力加速度和引力梯度張量的閉合解，被譽(yù)為小天體引力場(chǎng)的優(yōu)美解。該方法假設(shè)小天體內(nèi)部的密度函數(shù)為常值，可尋找某個(gè)矢量場(chǎng)函數(shù)使得其關(guān)于小天體內(nèi)質(zhì)量微元坐標(biāo)的散度恰好等于Newton核，利用Gauss散度定理將原計(jì)算引力勢(shì)的體積分轉(zhuǎn)化為關(guān)于小天體表面的面積分，式（1）可重寫為

假設(shè)小天體形狀可以由多面體近似表征，計(jì)算引力勢(shì)的面積分可表示為有限個(gè)三角形面的積分之和。

其中：Ff表示多面體的第f三角形面積區(qū)域；rf表示空間點(diǎn)r到第f三角形面上任意一點(diǎn)的矢量；表示第f三角形面的單位外法線矢量。

通過(guò)Green第一恒等式的進(jìn)一步表示以及面積分和線積分的轉(zhuǎn)化，最終得到任意均值多面體的引力勢(shì)、引力加速度和引力梯度張量

多面體法將小天體的引力場(chǎng)計(jì)算轉(zhuǎn)化成了空間點(diǎn)r到多面體每個(gè)棱邊和三角形面之間的幾何關(guān)系計(jì)算。當(dāng)使用較高分辨率的小天體多面體模型時(shí)，引力場(chǎng)的計(jì)算代價(jià)也通常較大。當(dāng)空間點(diǎn)r位于多面體棱邊上時(shí)多面體法出現(xiàn)奇異而無(wú)法計(jì)算。實(shí)際上棱邊上的奇異點(diǎn)都是可去奇異點(diǎn)，多面體引力勢(shì)在棱邊上仍然是解析的，Tsoulis等[38]針對(duì)這一問(wèn)題對(duì)引力勢(shì)解析表達(dá)式進(jìn)行了改進(jìn)，消除了表達(dá)式的奇異部分，從而完善了多面體法。

1.3.2 非均質(zhì)有限元多面體法

盡管多面體法是引力場(chǎng)的閉合解析形式，但在利用Gauss散度定理時(shí)引入了密度均勻假設(shè)，這導(dǎo)致多面體法無(wú)法反映小天體內(nèi)部質(zhì)量分布特征。實(shí)際上在小天體的形成過(guò)程中，物質(zhì)的聚集往往是雜亂隨機(jī)的，因此小天體內(nèi)部的質(zhì)量分布通常是非均勻的。為表征小天體內(nèi)部質(zhì)量分布情況，學(xué)者們采用有限元的思想，將原多面體分割成有限個(gè)更小的多面體，每個(gè)多面體賦予不同的密度值，最后通過(guò)疊加所有小多面體的解析引力場(chǎng)來(lái)得到真實(shí)引力場(chǎng)的近似值。Takahashi等[39]將小天體的多面體模型切割成有限個(gè)楔形體，通過(guò)球諧系數(shù)或者在軌觀測(cè)數(shù)據(jù)反演每個(gè)楔形體的密度值。

一般來(lái)說(shuō)，用于近似小天體外形的多面體具有較多的面數(shù)和邊數(shù)，利用多面體法求解引力場(chǎng)時(shí)需要遍歷計(jì)算所有面和邊，這意味著該方法需要較大的計(jì)算成本。若采用有限元結(jié)合多面體法來(lái)求解非均質(zhì)分布小天體的引力場(chǎng)，原二維流形問(wèn)題變?yōu)槿S問(wèn)題，需要遍歷的面和邊的個(gè)數(shù)將劇烈增加，隨之增加的計(jì)算成本通常難以承受。確定小天體內(nèi)部每個(gè)有限單元的密度實(shí)際上是經(jīng)典的引力場(chǎng)反問(wèn)題，該問(wèn)題并不存在唯一解，因此很難進(jìn)行求解。為說(shuō)明同一個(gè)引力場(chǎng)可對(duì)應(yīng)多個(gè)質(zhì)量分布，Takahashi等[39]描述了一種假想情況。學(xué)者們只能預(yù)先猜測(cè)小天體內(nèi)部可能的質(zhì)量分布規(guī)律，再進(jìn)一步進(jìn)行參數(shù)反演，但由于引力場(chǎng)反問(wèn)題的不適定本質(zhì)使得微小的測(cè)量誤差能夠?qū)е潞艽蟮膮?shù)不確定性[36]。

1.4 虛擬邊界元法

Wei等[40]發(fā)展了小天體引力場(chǎng)的虛擬邊界元全局表征方法，有效解決了上述方法中不能同時(shí)滿足：①表征的全局性；②表征的精確性；③表征參數(shù)的易求解性的問(wèn)題。該方法從求解Laplace方程的邊值問(wèn)題的角度對(duì)不規(guī)則非均質(zhì)小天體引力場(chǎng)進(jìn)行建模，即小天體引力勢(shì)函數(shù)是式（20）邊值問(wèn)題的解：

其中：Ω′表示小天體外部空間區(qū)域；n為小天體表面單位外法線矢量。

根據(jù)偏微分方程的基本理論，上述邊值問(wèn)題的解可表示為單層勢(shì)。

其中：σ(x)為定義在邊界上的待求面密度函數(shù)；G(x,r)為L(zhǎng)aplace方程的基本解，或稱為Green函數(shù)。

面密度函數(shù)σ(x)可以通過(guò)下述適定的第一類Fredholm積分方程求解。

為能夠快速獲得面積分(21)的具體數(shù)值，Wei和Shang[40]推導(dǎo)了單層勢(shì)的半解析形式。首先假設(shè)小天體邊界能夠由多面體P 近似。

其中：N為多面體三角形面的總數(shù)；Γi是第i三角形面。

對(duì)于任意一三角形面 Γi，存在一個(gè)光滑映射 φi使得φi(T)=Γi，T ?R2是標(biāo)準(zhǔn)三角形區(qū)域：

根據(jù)上述多面體的定義，式（21）可被重寫為

一種簡(jiǎn)單的光滑映射φi可以選為仿射變換。

其中：ξ3=1-ξ1-ξ2，Vi,1、Vi,2和Vi,3分別為第i三角形面上的三個(gè)頂點(diǎn)坐標(biāo)矢量。

利用三角形單元上的線性插值近似估計(jì)面密度函數(shù)，即

其中：hi,1、hi,2和hi,3面密度函數(shù)分別在第i三角形面上三個(gè)頂點(diǎn)處的值。

故單層勢(shì)可進(jìn)一步表示為

函數(shù)ψi,k(ξ2,r)具有閉合的解析形式，這里不贅述其結(jié)果。根據(jù)ψi,k(ξ2,r)及其關(guān)于ξ2的導(dǎo)數(shù)在0和1處的性質(zhì)，可以采用特殊的多項(xiàng)式近似逼近。

其中：Np為多項(xiàng)式的最高次數(shù)；系數(shù)an和bn可根據(jù)ψi,k(ξ2,r)的閉合解析形式構(gòu)建線性方程組求解；多項(xiàng)式P(ξ2)用于近似ψi,1(ξ2,r)和ψi,3(ξ2,r)；多項(xiàng)式Q(ξ2)用于近似ψi,2(ξ2,r)。從而式（30）中的函數(shù)Φi,k(r)可寫成如下半解析形式：

當(dāng)空間點(diǎn)r靠近小天體邊界時(shí)，由于Green函數(shù)G(x,r)的奇異性，式（34）和（35）的半解析形式精度會(huì)迅速退化。為了解決這一問(wèn)題，Wei等[40]采用了虛擬邊界法，其基本思想是利用函數(shù)的解析延拓概念，將原邊值問(wèn)題的解的定義域擴(kuò)充到原定義域之外，在定義域外建立一個(gè)虛擬邊界，從而空間觀察點(diǎn)永遠(yuǎn)不會(huì)靠近邊界，奇異性被根本消除。圖 4給出了虛擬邊界概念示意圖，在小天體內(nèi)部建立一個(gè)虛擬邊界，使得在該邊界定義的Laplace方程邊值問(wèn)題的解與原小天體的引力勢(shì)函數(shù)相等，從而半解析形式（34）和（35）的全局精度能夠被保證。

圖4 虛擬邊界概念示意圖Fig.4 Illustration of concept of virtual boundary

這一方法被應(yīng)用于電磁彈性材料問(wèn)題和熱傳導(dǎo)問(wèn)題的求解，然而通過(guò)虛擬邊界構(gòu)建的邊界積分方程的適定性沒(méi)有被嚴(yán)格證明，這限制了該方法在許多偏微分方程問(wèn)題求解的應(yīng)用。Wei等[40]首先利用調(diào)和函數(shù)的基本性質(zhì)證明了引力勢(shì)函數(shù)于小天體內(nèi)部解析延拓的存在性和唯一性，利用Lax-Milgram定理證明了虛擬邊界積分方程的適定性，即

虛擬邊界元法與多面體法類似，需要遍歷所有虛擬邊界多面體的頂點(diǎn)和面的計(jì)算，當(dāng)多面體分辨率較高時(shí)會(huì)帶來(lái)較大的計(jì)算代價(jià)。

1.5 基于數(shù)據(jù)驅(qū)動(dòng)的小天體引力場(chǎng)建模方法

針對(duì)模型計(jì)算精度和計(jì)算效率不能兼顧的問(wèn)題，學(xué)者們提出了基于數(shù)據(jù)驅(qū)動(dòng)的小天體引力場(chǎng)建模方法，旨在進(jìn)一步提高引力場(chǎng)的計(jì)算效率。這些方法包括兩大類：插值法和機(jī)器學(xué)習(xí)法。

1.5.1 插值法

空間局部插值方法試圖通過(guò)擴(kuò)大內(nèi)存消耗來(lái) 提高引力場(chǎng)的計(jì)算效率，是將空間有限離散點(diǎn)處的引 力場(chǎng)信息（如引力勢(shì)、引力加速度和引力梯度張量）事先存儲(chǔ)在內(nèi)存中，再通過(guò)插值計(jì)算得到空間其他任 意點(diǎn)的引力場(chǎng)信息。這一思想最早由Junkins[41]于 1976年在研究地球引力場(chǎng)時(shí)提出，后來(lái)由Engels等[42]一 同主張通過(guò)采用加權(quán)正交多項(xiàng)式作為插值基函數(shù)來(lái)提升插值精度。隨著計(jì)算機(jī)硬件的發(fā)展，內(nèi)存得到了迅速提升，基于空間局部插值方法被廣泛應(yīng)用于地球引力場(chǎng)的表征[43-45]。

Colombi等[46]提出了自適應(yīng)八叉樹空間網(wǎng)格結(jié)構(gòu)，從而將這一思想應(yīng)用到了小天體引力場(chǎng)的表征問(wèn)題。自適應(yīng)八叉樹空間網(wǎng)格劃分策略能夠很好地抑制空間插值誤差，使得全空間的引力場(chǎng)插值精度得到控制。Wei等[47]發(fā)展了結(jié)合復(fù)合引力逼近與快速網(wǎng)格跟蹤的高精度高效率引力場(chǎng)建模方法。他們基于復(fù)合引力思想將引力場(chǎng)拆解為由橢球表征的主項(xiàng)與局部同胚網(wǎng)格插值表征的擾動(dòng)項(xiàng)，針對(duì)擾動(dòng)項(xiàng)，采用了15節(jié)點(diǎn)的楔形體單元進(jìn)行插值近似。

為使得擾動(dòng)項(xiàng)對(duì)全空間的“貢獻(xiàn)”盡可能小，提出了引力最優(yōu)橢球概念，證明了引力最優(yōu)橢球求解問(wèn)題等價(jià)于橢球與小天體多面體的布爾體積最小化問(wèn)題，并推導(dǎo)了布爾體積的計(jì)算方法?；谛切味嗝骟w的特征發(fā)展了相容性空間網(wǎng)格劃分方法，進(jìn)一步，為了實(shí)現(xiàn)在軌道遞推過(guò)程中當(dāng)前點(diǎn)的快速準(zhǔn)確網(wǎng)格定位，提出相容性網(wǎng)格數(shù)據(jù)結(jié)構(gòu)和基于碰撞檢測(cè)的網(wǎng)格搜索方法。

1.5.2 機(jī)器學(xué)習(xí)法

近些年來(lái)，隨著計(jì)算機(jī)科學(xué)的發(fā)展與人工智能領(lǐng)域的興起，基于數(shù)據(jù)驅(qū)動(dòng)的機(jī)器學(xué)習(xí)方法在小天體的引力場(chǎng)計(jì)算中得到了廣泛的應(yīng)用。這類方法通過(guò)對(duì)小天體附近空間場(chǎng)點(diǎn)的數(shù)據(jù)采樣，利用機(jī)器學(xué)習(xí)算法完成模型的訓(xùn)練，實(shí)現(xiàn)對(duì)小天體附近引力場(chǎng)的回歸計(jì)算，克服了傳統(tǒng)方法中計(jì)算精度和計(jì)算效率難以兼顧的問(wèn)題。

Gao等[48]基于高斯過(guò)程回歸，建立了場(chǎng)點(diǎn)位置坐標(biāo)與引力場(chǎng)的直接映射關(guān)系，避免了復(fù)雜的建模計(jì)算過(guò)程。高斯過(guò)程的一個(gè)特點(diǎn)就是，對(duì)于每個(gè)變量x都有一個(gè)對(duì)應(yīng)的高斯分布，而對(duì)于一組變量X={x1,x2,···,xn}，存在一個(gè)聯(lián)合高斯分布，滿足

對(duì)于一個(gè)新的樣本xn+1，則新的高斯分布為

其中：

f是回歸的目標(biāo)函數(shù)；m是均值函數(shù)；k是高斯核函數(shù)。

在小天體的引力場(chǎng)建模計(jì)算中，回歸的目標(biāo)函數(shù)是小天體的引力加速度函數(shù)，變量x是小天體附近的場(chǎng)點(diǎn)位置。

人工神經(jīng)網(wǎng)絡(luò)得益于其靈活的網(wǎng)絡(luò)結(jié)構(gòu)，如圖 5所示，在函數(shù)的擬合方面具有良好的性能，因此往往被應(yīng)用于擬合小天體的引力場(chǎng)函數(shù)。Yu等[49]基于深度神經(jīng)網(wǎng)絡(luò)，擬合了小天體的引力加速度函數(shù)，同時(shí)，推導(dǎo)了深度神經(jīng)網(wǎng)絡(luò)引力加速度的偏導(dǎo)數(shù)，建立了引力梯度矩陣的近似形式，并應(yīng)用于求解燃料最優(yōu)的小天體著陸問(wèn)題。Cheng等[50]同樣地采用深度神經(jīng)網(wǎng)絡(luò)，實(shí)現(xiàn)了小天體引力場(chǎng)的快速計(jì)算，并解決了時(shí)間最優(yōu)的小天體著陸問(wèn)題。Furfaro等[51]采用單層反饋神經(jīng)網(wǎng)絡(luò)，建立了場(chǎng)點(diǎn)位置與小天體引力加速度的映射關(guān)系，同時(shí)，得益于單層反饋神經(jīng)網(wǎng)絡(luò)簡(jiǎn)潔的結(jié)構(gòu)，避免了神經(jīng)網(wǎng)絡(luò)訓(xùn)練的迭代優(yōu)化，有效地降低了訓(xùn)練神經(jīng)網(wǎng)絡(luò)的時(shí)間成本。

圖5 人工神經(jīng)網(wǎng)絡(luò)示意圖Fig.5 Illustration of artificial neural network

2 雙小行星系統(tǒng)引力場(chǎng)建模技術(shù)

除了單小天體，太陽(yáng)系還普遍存在大量雙小行星系統(tǒng)，針對(duì)它們的探測(cè)與研究也同樣對(duì)太陽(yáng)系起源、行星形成與演化甚至行星防御具有重要意義。雙小行星系統(tǒng)的引力場(chǎng)建模技術(shù)能為雙星系統(tǒng)探測(cè)任務(wù)設(shè)計(jì)與穩(wěn)定性分析提供動(dòng)力學(xué)仿真環(huán)境，其中包含的雙星系統(tǒng)動(dòng)力學(xué)建模技術(shù)研究也是分析該系統(tǒng)形成與演化機(jī)制的基礎(chǔ)。如圖 6所示，位于雙小行星系統(tǒng)附近空間一點(diǎn)r的探測(cè)器所感受到的小天體引力勢(shì)定義為

圖6 雙小行星系統(tǒng)示意圖Fig.6 Illustration of binary asteroid system

其中：si是小天體Bi內(nèi)部一點(diǎn)相對(duì)于雙星系統(tǒng)質(zhì)心O的位置矢量；ρi是兩個(gè)小天體內(nèi)部的密度分布函數(shù)；r是探測(cè)器相對(duì)于雙星系統(tǒng)質(zhì)心O的位置矢量。

基于單個(gè)小天體引力場(chǎng)建模技術(shù)，探測(cè)器在雙星系統(tǒng)附近的引力勢(shì)是在兩個(gè)小天體引力場(chǎng)中引力勢(shì)的疊加

其中：Ui是兩個(gè)小天體各自的引力勢(shì)；ri是兩個(gè)小天體各自的質(zhì)心Oi相對(duì)于系統(tǒng)質(zhì)心O的位置矢量；Ai是由小天體Bi固連坐標(biāo)系轉(zhuǎn)至慣性坐標(biāo)系的姿態(tài)矩陣。

雙小行星系統(tǒng)附近的引力場(chǎng)與兩個(gè)小天體的指向和位置密切相關(guān)，研究?jī)蓚€(gè)小天體的姿態(tài)運(yùn)動(dòng)和軌道運(yùn)動(dòng)是雙小行星系統(tǒng)引力場(chǎng)建模技術(shù)的關(guān)鍵。小天體的姿軌運(yùn)動(dòng)是由系統(tǒng)的動(dòng)力學(xué)決定的，考慮兩個(gè)小天體不規(guī)則的形狀和非均勻的質(zhì)量分布，這個(gè)系統(tǒng)的動(dòng)力學(xué)往往呈現(xiàn)出姿軌耦合的特點(diǎn)，即全二體問(wèn)題。由雙小行星系統(tǒng)的運(yùn)動(dòng)學(xué)可知，雙小行星系統(tǒng)的引力相互勢(shì)是研究全二體問(wèn)題的基礎(chǔ)。然而，在一般情況下，兩個(gè)不規(guī)則形狀且非均勻質(zhì)量分布的小天體之間的相互勢(shì)沒(méi)有閉合形式解。過(guò)去的幾十年中，不同學(xué)者們都致力于發(fā)展能夠計(jì)算雙星系統(tǒng)引力相互勢(shì)的方法，本章主要回顧了這些技術(shù)的發(fā)展歷程。

2.1 雙小行星系統(tǒng)的運(yùn)動(dòng)學(xué)

因?yàn)榻⒘讼嗷?shì)與兩個(gè)小天體之間的引力、引力矩的偏導(dǎo)數(shù)關(guān)系，Maciejewski[52]的雙剛體相對(duì)運(yùn)動(dòng)方程在雙小行星系統(tǒng)建模中得到了廣泛的應(yīng)用。雙小行星系統(tǒng)的相對(duì)運(yùn)動(dòng)示意圖如圖 7所示。雙剛體相對(duì)運(yùn)動(dòng)方程利用兩個(gè)小天體的相對(duì)軌道和姿態(tài)變量，減少了系統(tǒng)狀態(tài)變量的數(shù)量，簡(jiǎn)化了問(wèn)題的計(jì)算，在其中一個(gè)小天體B2固連坐標(biāo)系中的形式為

圖7 雙小行星系統(tǒng)的相對(duì)運(yùn)動(dòng)Fig.7 Relative motion in binary asteroid system

其中：P是系統(tǒng)的線動(dòng)量；R是小天體B1相對(duì)于小天體B2的位置矢量；Γ1和Γ2分別是小天體B1和B2的自轉(zhuǎn)角動(dòng)量；A是由小天體B2固連坐標(biāo)系轉(zhuǎn)至小天體B1固連坐標(biāo)系的姿態(tài)矩陣；A2是由小天體B2固連坐標(biāo)系轉(zhuǎn)至慣性坐標(biāo)系的姿態(tài)矩陣；Ω1和Ω2分別是小天體B1和B2的自轉(zhuǎn)角速度。

引力相互勢(shì)V的基本形式是兩個(gè)天體的質(zhì)量微元之間萬(wàn)有引力勢(shì)在兩個(gè)天體上的積分，如圖8所示。

圖8 引力相互勢(shì)示意圖Fig.8 Illustration of gravitational mutual potential

其中：G是萬(wàn)有引力常數(shù)，D1和D2分別是兩個(gè)小天體內(nèi)質(zhì)量微元在各自小天體固連系下的位置矢量。

兩個(gè)小天體之間的引力和引力矩可以通過(guò)相互勢(shì)函數(shù)的偏導(dǎo)數(shù)表示，故建立準(zhǔn)確的相互勢(shì)模型是研究雙星系統(tǒng)的基礎(chǔ)。

2.2 球形限制假設(shè)模型

在雙星系統(tǒng)引力和小天體結(jié)構(gòu)強(qiáng)度的共同作用下，一部分雙小行星系統(tǒng)的主星接近球形，次星接近橢球形，且與主星潮汐鎖定。在早期的雙小行星系統(tǒng)引力相互勢(shì)建模中，通常將其中一顆小天體假設(shè)為勻質(zhì)球體，其萬(wàn)有引力勢(shì)等價(jià)于一個(gè)質(zhì)點(diǎn)的引力勢(shì)，這一假設(shè)往往稱為球形限制假設(shè)（Sphere Restriction Assumption），假設(shè)小天體1是均質(zhì)球體，相互勢(shì)表示為

其中：M1是球形天體的質(zhì)量。此時(shí)，雙星系統(tǒng)的相互勢(shì)就退化為一個(gè)質(zhì)點(diǎn)在另一個(gè)不規(guī)則小天體引力場(chǎng)中的引力勢(shì)。

若假設(shè)不規(guī)則小天體具有簡(jiǎn)單的幾何形狀，如三軸橢球，則該小天體的引力勢(shì)函數(shù)具有簡(jiǎn)單的形式。進(jìn)而，能夠建立雙星系統(tǒng)的相互勢(shì)與系統(tǒng)模型參數(shù)（如形狀尺寸參數(shù)、相對(duì)距離參數(shù)、相對(duì)指向參數(shù)、質(zhì)量參數(shù)等）之間的定量關(guān)系。利用這種定量關(guān)系，分析了雙星系統(tǒng)的平衡構(gòu)型[53-54]、能量[55]，和接觸雙星系統(tǒng)的分裂[56]等。

2.3 級(jí)數(shù)展開模型

隨著雙星系統(tǒng)研究的發(fā)展，球形限制假設(shè)下對(duì)小天體形狀的簡(jiǎn)化難以滿足雙星系統(tǒng)建模精度的要求。通過(guò)級(jí)數(shù)展開將相互勢(shì)展開為無(wú)窮級(jí)數(shù)的形式，相互勢(shì)可表示為雙星系統(tǒng)相對(duì)位置參數(shù)、相對(duì)姿態(tài)參數(shù)、天體質(zhì)量分布參數(shù)的函數(shù)。進(jìn)而，通過(guò)相互勢(shì)的偏導(dǎo)數(shù)計(jì)算雙星系統(tǒng)引力和引力矩，基于級(jí)數(shù)展開的模型在雙星系統(tǒng)的遞推和演化分析中發(fā)揮了重要作用。

在全二體問(wèn)題背景下，Borderies[57]將其中一個(gè)天體的引力勢(shì)用球諧級(jí)數(shù)表示，將球諧級(jí)數(shù)在另一個(gè)天體的球坐標(biāo)系中表示后，得到了基于兩個(gè)天體的球諧系數(shù)的相互勢(shì)表達(dá)形式。Compère 等[58]基于STF（symmetric trace free）張量，推導(dǎo)了一種與球諧方法等價(jià)的方法，并應(yīng)用于雙小行星系統(tǒng)的遞推。Tricarico[59]基于冪級(jí)數(shù)推導(dǎo)了任意形狀和質(zhì)量分布的兩個(gè)小天體之間的相互勢(shì)，并證明了球諧系數(shù)與慣性積分作為天體質(zhì)量分布參數(shù)的一致性。利用有限階球諧系數(shù)和慣性積分的相互勢(shì)、引力和引力矩的建模工作還可以參考文獻(xiàn)[60～63]。

針對(duì)形狀復(fù)雜的雙小行星系統(tǒng)，學(xué)者們發(fā)展了一系列高階級(jí)數(shù)展開方法，以提高建模精度。Werner 等[64]采用均質(zhì)多面體描述兩個(gè)小天體，通過(guò)疊加構(gòu)成兩個(gè)小天體任意四面體對(duì)之間的引力勢(shì)解決總相互勢(shì)的評(píng)估。利用勒讓德級(jí)數(shù)來(lái)展開兩個(gè)四面體之間微元的距離項(xiàng)，實(shí)現(xiàn)四面體對(duì)之間引力勢(shì)的解析可積，得到基于各階慣性張量的引力勢(shì)

其中：a和b分別代表兩個(gè)小天體中以多面體三角面為底面，質(zhì)心為頂點(diǎn)的四面體；ρa(bǔ)和ρb分別代表兩個(gè)小天體的密度；Ta和Tb分別代表四面體與標(biāo)準(zhǔn)四面體轉(zhuǎn)換的雅克比行列式；Q是常系數(shù)張量；w和r是與兩個(gè)小天體的位置、姿態(tài)相關(guān)的張量。

Fahnestock等[65]進(jìn)一步地將Werner的工作推廣至相互引力和引力矩的計(jì)算，并應(yīng)用于1999KW4雙小行星系統(tǒng)的動(dòng)力學(xué)分析[66]。該方法對(duì)小天體形狀描述的準(zhǔn)確程度取決于多面體模型的分辨率，當(dāng)多面體模型的分辨率很高時(shí)其計(jì)算代價(jià)是較大的。

為提高相互勢(shì)評(píng)估的計(jì)算效率，Hou等[67]以兩個(gè)天體邊界為積分區(qū)間，利用勒讓德級(jí)數(shù)對(duì)構(gòu)成兩個(gè)天體的任意微元間的距離進(jìn)行逼近，從而將共同勢(shì)表示為兩個(gè)天體廣義慣性積的函數(shù)

其中：t、a和b分別是滿足遞歸關(guān)系的常系數(shù)；分別是小天體B2固連坐標(biāo)系下，小天體B1和B2的廣義慣性積分。由于兩個(gè)小天體存在相對(duì)轉(zhuǎn)動(dòng)，是時(shí)變的，該方法通過(guò)坐標(biāo)變換建立了兩個(gè)小天體固連坐標(biāo)系下廣義慣性積分之間的關(guān)系，從而減少了重復(fù)計(jì)算，提高了計(jì)算效率。廣義慣性積分的引入不僅避免了大量系統(tǒng)狀態(tài)相關(guān)參數(shù)的重復(fù)計(jì)算，大幅度降低了共同勢(shì)的計(jì)算代價(jià)，還方便地描述了小天體不規(guī)則的形狀和不均勻的質(zhì)量分布。此外，該模型具備成熟的遞歸公式體系，能夠便于計(jì)算機(jī)的求解。

Shi等[68]將主星的引力場(chǎng)用多面體模型描述，通過(guò)在從星質(zhì)心處利用泰勒級(jí)數(shù)對(duì)主星引力場(chǎng)進(jìn)行展開，將共同勢(shì)表示為從星慣性積和主星引力勢(shì)及其導(dǎo)數(shù)的函數(shù)，給出了相互勢(shì)展開至二階項(xiàng)的形式

其中：U是主星的多面體引力勢(shì)函數(shù)；??U是主星的引力勢(shì)能梯度陣；? 表示張量計(jì)算；P是與從星的二階慣性積分和兩個(gè)小天體相對(duì)姿態(tài)相關(guān)的并矢

該模型能夠?qū)ζ渲幸活w星進(jìn)行不規(guī)則引力場(chǎng)的精確建模，避免了同時(shí)舍去兩顆星的高階慣性積分。

基于級(jí)數(shù)展開的模型能夠計(jì)算任意形狀和質(zhì)量分布的小天體之間的相互勢(shì)，通過(guò)計(jì)算引力勢(shì)偏導(dǎo)數(shù)的方法得到引力和引力矩。該模型的誤差來(lái)源于級(jí)數(shù)的截?cái)嗾`差，隨著兩個(gè)小天體之間的距離減小，級(jí)數(shù)的截?cái)嗾`差增大，且此類模型的計(jì)算量隨著級(jí)數(shù)階數(shù)的提高顯著增大。

2.4 有限元模型

針對(duì)基于級(jí)數(shù)展開的模型在兩個(gè)小天體距離接近時(shí)，級(jí)數(shù)收斂變慢甚至發(fā)散。于洋等[69]提出了一種有限元方法來(lái)改善相互勢(shì)、引力和力矩評(píng)估的收斂性和精度。這種方法并不依賴于相互勢(shì)的偏導(dǎo)數(shù)計(jì)算引力和引力矩，而是根據(jù)直接的物理含義對(duì)兩個(gè)小天體之間的引力相互作用進(jìn)行建模。該方法通過(guò)對(duì)兩個(gè)小天體進(jìn)行四面體有限元?jiǎng)澐郑秒p線性插值對(duì)構(gòu)成兩個(gè)小天體任意兩個(gè)微元之間的相互勢(shì)、力和力矩進(jìn)行計(jì)算。兩個(gè)天體之間總的相互勢(shì)、力和力矩可以通過(guò)微元間引力作用的簡(jiǎn)單疊加獲得。

通過(guò)有限元的細(xì)致劃分，基于有限元方法的模型對(duì)于小天體復(fù)雜的幾何構(gòu)型和局部的密度異常都是有效的。有限元模型的計(jì)算精度取決于有限元?jiǎng)澐值木?xì)程度，當(dāng)模型的有限元?jiǎng)澐志容^高時(shí)，該模型計(jì)算精度較高，但是大量的有限元單元也導(dǎo)致該方法的計(jì)算負(fù)擔(dān)較大。高云峰等[70]在于洋工作的基礎(chǔ)上利用CUDA并行計(jì)算框架進(jìn)一步提升了這一方法的計(jì)算效率，在相同的精度下，與單核CPU相比，計(jì)算效率提升了兩個(gè)量級(jí)。

3 未來(lái)發(fā)展趨勢(shì)

小天體不規(guī)則形狀和其內(nèi)部不均勻的質(zhì)量分布給小天體引力場(chǎng)計(jì)算帶來(lái)了極大挑戰(zhàn)。經(jīng)歷了兩個(gè)世紀(jì)的漫長(zhǎng)發(fā)展，已經(jīng)就小天體的引力場(chǎng)建模問(wèn)題有了較為全面的認(rèn)識(shí)，但仍有一些問(wèn)題需要進(jìn)一步研究和解決。

3.1 引力場(chǎng)的解析建模

在天體引力場(chǎng)研究中，最先發(fā)展的是球諧函數(shù)展開法。該方法利用了Legendre多項(xiàng)式和三角函數(shù)構(gòu)成的無(wú)窮級(jí)數(shù)，在球坐標(biāo)下描述了任意天體的引力勢(shì)函數(shù)、引力加速度和引力梯度張量等。由于Legendre展開式收斂條件的限制，球諧函數(shù)級(jí)數(shù)最終只能在天體的最小包圍球（即Brillouin球）外部保證收斂，而球內(nèi)部收斂無(wú)法保證。小天體普遍呈現(xiàn)不規(guī)則的形狀，使得Brillouin球無(wú)法貼合小天體表面，導(dǎo)致小天體表面附近的空間引力場(chǎng)無(wú)法被描述。

學(xué)者們意識(shí)到球諧函數(shù)展開法的局限性后，利用橢球坐標(biāo)系的描述將Brillouin球變換成Brillouin橢球，從而能夠更好地貼合非球形小天體。對(duì)應(yīng)的球諧函數(shù)被變換成了橢球諧函數(shù)，它由第一類和第二類Lamé函數(shù)構(gòu)成。相比球諧函數(shù)，橢球諧函數(shù)具有更高的計(jì)算復(fù)雜度，且還未發(fā)現(xiàn)Lamé函數(shù)的遞歸公式，這限制了橢球諧函數(shù)的應(yīng)用。另一方面，雖然Brillouin橢球能夠更好地貼合非球形小天體的形狀，但小天體表面絕大部分空間的引力場(chǎng)依然無(wú)法描述，因此橢球諧函數(shù)仍未解決無(wú)窮級(jí)數(shù)的非全局表征問(wèn)題。

在過(guò)去的幾十年中，學(xué)者們開始重點(diǎn)關(guān)注小天體表面的引力場(chǎng)全局表征問(wèn)題，如無(wú)窮級(jí)數(shù)法中的內(nèi)外引力場(chǎng)法和球Bessel函數(shù)展開法、質(zhì)點(diǎn)群法、多面體法和有限元多面體法。這些方法均有側(cè)重，例如內(nèi)外引力場(chǎng)法側(cè)重于將級(jí)數(shù)收斂域擴(kuò)展到小天體表面一點(diǎn)；球Bessel函數(shù)展開法側(cè)重于Brillouin球內(nèi)部的低階級(jí)數(shù)表達(dá)；質(zhì)點(diǎn)群法的精度隨觀察點(diǎn)靠近小天體表面迅速退化，因此一般也用于較遠(yuǎn)距離的軌道動(dòng)力學(xué)分析；多面體法雖然是全局的解析閉合引力場(chǎng)表征方法，但該方法不能表征非均質(zhì)小天體引力場(chǎng)，一般僅用于理想均質(zhì)多面體附近的動(dòng)力學(xué)分析；有限元多面體法雖然能夠近似表征非均質(zhì)小天體的引力場(chǎng)，但每個(gè)有限單元的密度受到引力場(chǎng)反問(wèn)題不適定本質(zhì)的限制難以準(zhǔn)確求解?？偟膩?lái)說(shuō)這些方法均不能同時(shí)滿足4方面要求：①表征的全局性；②能夠表征非均質(zhì)分布；③表征參數(shù)的易求解性；④表征的高效性。虛擬邊界元法是有效解決前3個(gè)問(wèn)題的新方法，但該方法依然存在有待解決的問(wèn)題。對(duì)于任意給定小天體，虛擬邊界的構(gòu)建不唯一，不同虛擬邊界最終影響著引力場(chǎng)表征精度，該方法尚未給出虛擬邊界定義的標(biāo)準(zhǔn)；引力場(chǎng)計(jì)算精度只能通過(guò)增加虛擬邊界的分辨率提高，由于該方法的計(jì)算復(fù)雜度為O(N2)，則計(jì)算代價(jià)會(huì)隨著虛擬邊界多面體的分辨率提升而迅速增加。表 1總結(jié)了本文梳理的目前主要小天體引力場(chǎng)模型特點(diǎn)，可見尚未有同時(shí)滿足4方面要求的模型。

表1 不同小天體引力場(chǎng)模型特點(diǎn)比較Table 1 Comparison of the characteristics of gravitational field models of different small celestial bodies

3.2 基于數(shù)據(jù)驅(qū)動(dòng)方法的計(jì)算效率

其他方法則是針對(duì)如何提升小天體引力場(chǎng)計(jì)算效率的問(wèn)題提出的，這類方法包括機(jī)器學(xué)習(xí)法、八叉樹插值法和復(fù)合引力場(chǎng)逼近法等?；跀?shù)據(jù)驅(qū)動(dòng)的一類方法雖然計(jì)算效率得到大幅提升，但節(jié)點(diǎn)處的“真值”或訓(xùn)練集的采樣數(shù)據(jù)依然依賴某個(gè)具體的引力場(chǎng)表征模型給出。因此，該類方法并沒(méi)有從基本層面解決小天體引力場(chǎng)表征問(wèn)題。未來(lái)小天體引力場(chǎng)表征方法的發(fā)展趨勢(shì)是除保證方法的全局性、精確性和參數(shù)易求解性之外，還應(yīng)保證方法的高效性。

3.3 雙小行星系統(tǒng)動(dòng)力學(xué)模型的高精度與高效率計(jì)算

雙小行星系統(tǒng)對(duì)天體力學(xué)的科學(xué)探索和行星防御的工程實(shí)踐都具有重要的價(jià)值，近些年來(lái)成為了深空探測(cè)領(lǐng)域的熱點(diǎn)目標(biāo)。雙星系統(tǒng)的引力場(chǎng)建模技術(shù)能為雙星探測(cè)任務(wù)設(shè)計(jì)與穩(wěn)定性分析提供動(dòng)力學(xué)仿真環(huán)境。其中，研究?jī)蓚€(gè)小天體的姿軌運(yùn)動(dòng)是雙星系統(tǒng)引力場(chǎng)建模的基礎(chǔ)。根據(jù)雙星系統(tǒng)的運(yùn)動(dòng)學(xué)方程可知，引力相互勢(shì)的建模是解決全二體問(wèn)題的關(guān)鍵。然而，小行星不規(guī)則的形狀、不均勻的質(zhì)量分布和雙星系統(tǒng)獨(dú)特的構(gòu)型，給雙星系統(tǒng)引力相互勢(shì)的計(jì)算帶來(lái)了諸多挑戰(zhàn)。

早期的研究一般基于單星的引力場(chǎng)模型，對(duì)另一顆星進(jìn)行適當(dāng)?shù)暮?jiǎn)化，以得到閉合形式的引力相互勢(shì)。這類方法建立了雙星系統(tǒng)的相互勢(shì)與系統(tǒng)模型參數(shù)（如形狀尺寸參數(shù)、相對(duì)距離參數(shù)、相對(duì)指向參數(shù)、質(zhì)量參數(shù)等）之間的定量關(guān)系。有助于進(jìn)行雙星系統(tǒng)能量、穩(wěn)定構(gòu)型方面的初步分析。然而，對(duì)小天體的形狀或者引力場(chǎng)進(jìn)行簡(jiǎn)化損失了對(duì)系統(tǒng)描述的精度。級(jí)數(shù)展開類的方法在保留對(duì)小天體不規(guī)則形狀和不均勻質(zhì)量分布描述的基礎(chǔ)上，對(duì)相互勢(shì)進(jìn)行級(jí)數(shù)展開，相互勢(shì)一般被表示為無(wú)窮級(jí)數(shù)的形式。此類方法的誤差來(lái)源于級(jí)數(shù)的截?cái)嗾`差，隨著兩個(gè)小天體之間的距離減小，級(jí)數(shù)的截?cái)嗾`差增大，且此類模型的計(jì)算量隨著級(jí)數(shù)階數(shù)的提高顯著增大。為了克服級(jí)數(shù)展開類方法收斂困難的問(wèn)題，有限元方法在雙星建模中的應(yīng)用，為引力相互作用的計(jì)算提供了全新的思路。通過(guò)有限元?jiǎng)澐郑瑢蓚€(gè)小行星之間的相互勢(shì)表示為有限元微元對(duì)之間相互勢(shì)疊加的形式。通過(guò)細(xì)致的有限元?jiǎng)澐?，有限元方法能夠?zhǔn)確地描述兩個(gè)小行星的不規(guī)則形狀和不均勻質(zhì)量分布，并且截?cái)嗾`差可控。但是，該方法的計(jì)算負(fù)擔(dān)會(huì)隨著有限元的數(shù)量增加而增加。

基于單個(gè)小天體的引力場(chǎng)建模技術(shù)研究基礎(chǔ)，未來(lái)雙小行星系統(tǒng)附近的引力場(chǎng)建模技術(shù)的發(fā)展應(yīng)主要關(guān)注于雙小行星系統(tǒng)的動(dòng)力學(xué)研究。在雙小行星系統(tǒng)的動(dòng)力學(xué)建模方面，已有的研究主要關(guān)注于雙星系統(tǒng)的平衡構(gòu)型以及姿軌運(yùn)動(dòng)演化。然而，在這些方法中，模型的計(jì)算精度和計(jì)算效率往往是矛盾的。在未來(lái)的雙星系統(tǒng)引力場(chǎng)建模技術(shù)和雙星系統(tǒng)動(dòng)力學(xué)建模研究中，應(yīng)主要關(guān)注以下幾個(gè)方面。

1）對(duì)小天體不規(guī)則形狀和不均勻質(zhì)量分布的準(zhǔn)確描述。根據(jù)已探明的小天體可知，小天體的形狀不規(guī)則且質(zhì)量分布不均勻，存在風(fēng)化層等特殊結(jié)構(gòu)。由于雙小行星系統(tǒng)獨(dú)特的動(dòng)力學(xué)，雙星系統(tǒng)的演化受到小天體質(zhì)量分布和形狀的影響較為顯著，未來(lái)的模型應(yīng)能夠考慮小天體特殊的形態(tài)和物質(zhì)組成。

2）兼顧計(jì)算精度和計(jì)算效率。在現(xiàn)有的方法中，模型的計(jì)算精度和計(jì)算效率往往是矛盾的，較高的計(jì)算精度常常會(huì)導(dǎo)致沉重的計(jì)算負(fù)擔(dān)。發(fā)展兼顧計(jì)算精度和計(jì)算效率的方法，能夠提升在雙星系統(tǒng)遞推方面的實(shí)用性，以實(shí)現(xiàn)對(duì)雙星系統(tǒng)的長(zhǎng)期演化遞推和雙星系統(tǒng)附近引力場(chǎng)建模的有效應(yīng)用，例如，雙星系統(tǒng)附近粒子的運(yùn)動(dòng)演化規(guī)律研究等。

4 結(jié)束語(yǔ)

本文介紹了單個(gè)小天體的引力場(chǎng)建模技術(shù)和雙小行星系統(tǒng)附近的引力場(chǎng)建模技術(shù)。目前對(duì)目標(biāo)小天體附近的引力場(chǎng)建模技術(shù)研究已經(jīng)形成了比較成熟的理論體系，但是也仍有一些理論和技術(shù)問(wèn)題尚待解決或完善。發(fā)展兼顧計(jì)算精度和計(jì)算效率、具備全局表征能力的建模技術(shù)仍是該領(lǐng)域的研究重點(diǎn)。未來(lái)將小天體探測(cè)的實(shí)際數(shù)據(jù)與現(xiàn)有理論方法進(jìn)行對(duì)比，能夠在驗(yàn)證現(xiàn)有理論的基礎(chǔ)上，發(fā)現(xiàn)全新的科學(xué)、技術(shù)問(wèn)題，推動(dòng)小天體引力場(chǎng)建模方法和雙小行星系統(tǒng)動(dòng)力學(xué)建模方法的發(fā)展。