劉媛琪，周珊羽，鞏子嘉

(1.山東大學(xué)(威海) 空間科學(xué)與物理學(xué)院/空間科學(xué)研究院，山東 威海 264209；2.山東大學(xué)(威海) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東 威海 264209)

?

N體問題下日地三角拉格朗日點(diǎn)的數(shù)值研究

劉媛琪1，周珊羽1，鞏子嘉2

(1.山東大學(xué)(威海) 空間科學(xué)與物理學(xué)院/空間科學(xué)研究院，山東 威海 264209；2.山東大學(xué)(威海) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，山東 威海 264209)

針對日地系統(tǒng)的拉格朗日點(diǎn)相關(guān)的N體問題(N>3)目前尚無可用的解析理論解的問題，試圖通過數(shù)值方法，運(yùn)用攝動理論的思想，得出限制性三體問題下的日地三角拉格朗日點(diǎn)(L4、L5)在N體問題下的引力值(非零引力合力偏差)。通過對拉格朗日點(diǎn)的攝動規(guī)律進(jìn)行定性和定量分析，從來源、周期性等方面總結(jié)分析，為深空探測相關(guān)研究提供參考。

拉格朗日點(diǎn)；N體問題；數(shù)值方法

0 引言

1772年，意大利科學(xué)家拉格朗日推導(dǎo)出2個(gè)質(zhì)量相差懸殊的天體在同一平面上存在5個(gè)特殊的點(diǎn)，在這5個(gè)點(diǎn)上物體所受到的引力與物體運(yùn)動產(chǎn)生的向心力達(dá)到平衡，這5個(gè)點(diǎn)也被稱為“拉格朗日點(diǎn)”(記為L1～L5)(如圖1所示)。

從圖中可知，L4、L5點(diǎn)在距離太陽和地球相等距離的2個(gè)等邊三角形頂點(diǎn)處，所以也稱為三角拉格朗日點(diǎn)。眾所周知，三角形結(jié)構(gòu)具有穩(wěn)定性，科學(xué)家在1906年發(fā)現(xiàn)第617號小行星出現(xiàn)在木星軌道落后60°處；20世紀(jì)80年代，科學(xué)家又發(fā)現(xiàn)在土星和它的大衛(wèi)星所構(gòu)成的運(yùn)動系統(tǒng)中也有類似的等邊三角形。數(shù)學(xué)的推導(dǎo)和天文觀測結(jié)論證明了拉格朗日點(diǎn)的特殊性確實(shí)存在。

圖1 日地三體系統(tǒng)的拉格朗日點(diǎn)

日地拉格朗日點(diǎn)因與太陽和地球相對位置的不變性而成為國際空間探測的熱點(diǎn)。在拉格朗日點(diǎn)上，航天器可以基本保持穩(wěn)定的軌道而不消耗太多燃料來維持；因此拉格朗日點(diǎn)也是公認(rèn)的望遠(yuǎn)鏡和探測器的絕佳位置以及未來太空城的最佳位置。目前，威爾金森微波背景各向異性探測器已經(jīng)在日-地系統(tǒng)的L2點(diǎn)上運(yùn)行，詹姆斯韋伯太空望遠(yuǎn)鏡將要被放置在日-地系統(tǒng)的L2點(diǎn)上。我國的嫦娥二號在各項(xiàng)科學(xué)目標(biāo)都取得圓滿成功后開始了新的征程，奔向150萬km外的拉格朗日L2點(diǎn)。2011-08-25 T 23：27：00，經(jīng)過77 d的飛行，嫦娥二號在世界上首次實(shí)現(xiàn)從月球軌道出發(fā)，受控準(zhǔn)確進(jìn)入日-地系統(tǒng)的L2點(diǎn)環(huán)繞軌道，第一次實(shí)現(xiàn)中國對月球以外的太空進(jìn)行探測，是中國第一次開展拉格朗日點(diǎn)轉(zhuǎn)移軌道和使命軌道的設(shè)計(jì)和控制，并實(shí)現(xiàn)了150萬km遠(yuǎn)距離測控通信。

在天體力學(xué)中，行星系統(tǒng)的起源和演化問題一直是大家關(guān)注的重點(diǎn)。對這類問題，需要解N體運(yùn)動的微分方程；但這些方程一般都相當(dāng)復(fù)雜，除了二體問題等幾種少數(shù)情況外，都不能得到嚴(yán)格的分析解。目前常用的研究方法是攝動理論，先用二體問題近似，然后再進(jìn)一步考慮攝動因素影響。如為了得到精確的行星軌道，除了考慮太陽的影響，還要考慮太陽系中其他行星的引力作用，尤其是木星這樣的大質(zhì)量行星對其他行星產(chǎn)生的影響?，F(xiàn)代精密的天體測量技術(shù)已經(jīng)可以測定出對于二體運(yùn)動非常微小的偏離。

至今為止對日地系統(tǒng)三體限制性拉格朗日點(diǎn)的研究做了近似，尤其是忽略了其他天體的影響、將所有運(yùn)動軌道都近似看作在黃道面上。本文考慮了黃道坐標(biāo)系三維空間內(nèi)8大行星在拉格朗日點(diǎn)的引力作用，以期得到限制性三體問題下日地拉格朗日點(diǎn)在太陽系內(nèi)受行星影響的情況。

本文利用開普勒根數(shù)計(jì)算給定時(shí)間黃道坐標(biāo)系內(nèi)8大行星的精確位置。假設(shè)太陽和地球?qū)窭嗜拯c(diǎn)的引力和與其向心力完全相同，其他行星在拉格朗日點(diǎn)的加速度通過最基礎(chǔ)的加速度矢量疊加，即可求解拉格朗日點(diǎn)受到影響(數(shù)值)的規(guī)律。

通過研究發(fā)現(xiàn)，在黃道坐標(biāo)系內(nèi)，金星對拉格朗日點(diǎn)z方向上有一定周期性的擾動，與金星距離地球較近且有3°的軌道傾角情況相符。同時(shí)，在黃道面內(nèi)，拉格朗日點(diǎn)受木星擾動最大，且有一定的周期性規(guī)律，與木星距離地球近且質(zhì)量較大的情況相符。

通過攝動理論，可以解決許多天文學(xué)與空間科學(xué)中的重要問題，例如證明太陽系存在的第九大行星、深空探測衛(wèi)星軌道設(shè)計(jì)等。如能夠通過數(shù)值方法將攝動理論充分應(yīng)用，將對攝動理論的發(fā)展產(chǎn)生重要的影響。

1 計(jì)算過程

1.1 確定8大行星在給定日期的狀態(tài)矢量

其步驟如下：

1)求儒略日J(rèn)D[2]203-208。

根據(jù)儒略日的定義，用下式推算某日世界時(shí)0時(shí)刻的儒略日

(1)

式中：y為年份；m為月份；d為天數(shù)；int表示取整運(yùn)算符。

2)求J2000至給定日期的儒略世紀(jì)T0。

J2000表示儒略紀(jì)元法，本文中所用軌道根數(shù)均在儒略紀(jì)元2000時(shí)為起始值，并有世紀(jì)變化率。因此在計(jì)算時(shí)需得出J2000至給定日期的儒略世紀(jì)

(2)

3)求得軌道根數(shù)在T0時(shí)的值

(3)

4)求JD時(shí)的近日點(diǎn)幅角ω、真近點(diǎn)角θ[2]155-161。

ω=?-Ω

M=L-?，

(4)

式中：?為近日點(diǎn)緯度；Ω為升交點(diǎn)赤經(jīng)；M為平近點(diǎn)角。

開普勒方程為

M=E-e sinE。

(5)

式中：E為偏近點(diǎn)角；e為橢圓軌道偏心率。

可以通過迭代逼近求得偏近點(diǎn)角E為

(6)

即可求得真近點(diǎn)角θ。

5)由軌道根數(shù)求解狀態(tài)矢量R、v[1]224-228，[2]203-208。

①軌道面坐標(biāo)系下(以行星各自軌道面為z軸，軌道焦點(diǎn)為原點(diǎn))有

(7)

式中a為橢圓軌道半長軸。

②轉(zhuǎn)換為日心黃道坐標(biāo)系的變換矩陣為

(8)

式中：Ω為升交點(diǎn)赤經(jīng)；ω為近日點(diǎn)幅角；i為軌道傾角。則

(9)

1.2 L4、L5拉格朗日點(diǎn)的精確軌道

1)拉格朗日點(diǎn)L4：L4點(diǎn)在以限制性3體問題的2個(gè)主天體(本文中均在日地系統(tǒng)中進(jìn)行研究，分別指太陽和地球)連線為底邊的等邊三角形的第3個(gè)頂點(diǎn)上，在較小天體圍繞2天體系統(tǒng)質(zhì)心運(yùn)行軌道的前方，則

(10)

式中：R代表日心黃道坐標(biāo)系內(nèi)位置參量；下標(biāo)L4代表拉格朗日4點(diǎn)上測試質(zhì)量的位置參量；下標(biāo)3代表地球。由于L4點(diǎn)的軌道也在地球公轉(zhuǎn)平面上，因此在黃道坐標(biāo)系的xOy平面上可通過地球狀態(tài)參量旋轉(zhuǎn)得到其狀態(tài)參量。

2)拉格朗日點(diǎn)L5：L5點(diǎn)在以2個(gè)主天體連線為底邊的等邊三角形的第3個(gè)頂點(diǎn)上，在較小天體圍繞2天體系統(tǒng)質(zhì)心運(yùn)行軌道的后方。則有

(11)

式中下標(biāo)L5代表拉格朗日5點(diǎn)上測試質(zhì)量的位置參量。由于L5點(diǎn)的軌道也在地球公轉(zhuǎn)平面上，因此在黃道坐標(biāo)系的xOy平面上可通過地球狀態(tài)參量旋轉(zhuǎn)得到其狀態(tài)參量。

1.3 計(jì)算拉格朗日點(diǎn)受其他行星運(yùn)行的影響

1)目標(biāo)函數(shù)(向量)

(12)

式中：F為拉格朗日點(diǎn)軌道上實(shí)驗(yàn)質(zhì)量m受到的萬有引力合力矢量；下標(biāo)N-body代表N體問題條件下，3-body代表三體問題條件下。

2)三體問題下

(13)

式中:G為萬有引力常數(shù)G=6.67×10-11·m3·kg-1·s-2;M為天體質(zhì)量；下標(biāo)0為太陽；1～8依次代表由近及遠(yuǎn)的太陽系8大行星。方程右邊由萬有引力定律得出。

3)N體問題下

(14)

4)目標(biāo)函數(shù)

(15)

目標(biāo)函數(shù)是N體問題和三體問題萬有引力合力的差值，也就是三體問題下拉格朗日點(diǎn)受到行星的攝動力影響。

2 數(shù)值分析

2.1 L4點(diǎn)的受力分析

圖2所示為L4點(diǎn)所受太陽及8大行星加速度之和大小(標(biāo)量)的變化，橫軸單位為d(指經(jīng)過2015-01-01 T 00∶00∶00的天數(shù)，下文如無特別指出均為此含義)，縱軸單位為10-5m·s-2。取若干點(diǎn)對日地合加速度和其余行星合加速度數(shù)量級進(jìn)行計(jì)算，均有近10個(gè)數(shù)量級之差。由此可見，拉格朗日點(diǎn)受8大行星攝動影響較為有限，基本處于十分穩(wěn)定的狀態(tài)；但合加速度變化幅度有一定的周期性，可以進(jìn)行進(jìn)一步研究。

圖2 太陽及太陽系8大行星在L4的合加速度

圖3分別為除金星、木星和地球外5個(gè)行星(小圖題目所示)在拉格朗日L4點(diǎn)的攝動加速度大小的三維分布(深色)與L4點(diǎn)所受太陽及8大行星合加速度大小的三維分布(淺色)的對比圖。如圖3所示，以上5個(gè)行星的加速度比合加速度小至少1個(gè)數(shù)量級，如此微小的攝動可以忽略不計(jì)。

圖3 行星在L4引力加速度與合加速度比較圖

圖4中深色線為金星在L4點(diǎn)的加速度在10 000 d內(nèi)的軌跡圖。不難看出金星在L4點(diǎn)的引力作用產(chǎn)生的加速度與合加速度(淺色)的z分量周期性變化規(guī)律及極值大小幾乎相同。金星與地球距離較小，且具有3.4°的軌道傾角，相比其他行星，認(rèn)為其傾角足夠大而足以產(chǎn)生較為明顯的影響；所以不難理解L4點(diǎn)在z方向的攝動主要是受金星位置的影響。而且由于金星周期比地球周期小;所以z方向加速度的周期性變化十分明顯，峰值往往出現(xiàn)在金星近地點(diǎn)。接下來對其變化周期與極值點(diǎn)出現(xiàn)時(shí)金星與L4點(diǎn)的相對位置進(jìn)行分析討論。

圖5為木星對L4點(diǎn)引力作用產(chǎn)生加速度與合加速度在10 000 d內(nèi)的變化軌道圖。深色線為木星在L4點(diǎn)的加速度，與合加速度(淺色)在xOy平面分量非常接近。木星是距離地球最大的巨行星，軌道傾角1.3°，其對拉格朗日點(diǎn)在黃道面上的影響最大也是可以合理解釋的。雖然木星周期較大，但拉格朗日點(diǎn)的軌道周期與地球相同;所以其共線的周期在10 000 d內(nèi)有多次的體現(xiàn)。

圖4 金星在L4引力加速度與合加速度比較圖

圖5 木星在L4引力加速度與合加速度比較圖

圖6為L4點(diǎn)處金星萬有引力在10 000 d內(nèi)隨時(shí)間的變化。圖7為L4點(diǎn)處木星萬有引力隨時(shí)間的變化。與金星在L4點(diǎn)引力的極大值變化規(guī)律略有不同，木星在L4點(diǎn)引力作用產(chǎn)生的加速度極大值有明顯的浮動，這與木星軌道與地球軌道相比離心率較大(是地球的2.8倍)，同時(shí)半長軸較大(是地球的5倍)有關(guān)，與金星更近似為正圓(離心率是木星的1/8)的軌道差別較大。

圖6 金星在L4引力加速度標(biāo)量變化圖

圖7 木星在L4引力加速度標(biāo)量變化圖

根據(jù)以上分析可以看出，拉格朗日點(diǎn)攝動影響有明顯的周期性規(guī)律，既有長周期，又有短周期，且與三角函數(shù)等形式較為接近；因此通過時(shí)間序列譜分析可以提取其中隱含的部分信息并直觀理解其原因。本文選取木星在L4點(diǎn)引力加速度隨時(shí)間的變化作為示例。

通過離散傅里葉變換(discrete Fourier transform，DFT)把原始序列信號從時(shí)間域變換到頻率域，即把原始序列x(n)映射成新的序列

k=0,1,2,…,N-1。

(16)

設(shè)第k個(gè)點(diǎn)對應(yīng)頻率是fk，則

(17)

式中NΔt是原始點(diǎn)列x(n)中相鄰2個(gè)點(diǎn)橫坐標(biāo)的間隔，即采樣周期。

在對木星在L4點(diǎn)引力加速度隨時(shí)間的變化進(jìn)行分析時(shí)，取Δt=10，點(diǎn)數(shù)N=10 000，經(jīng)過離散傅里葉變換，得到頻譜圖如圖8所示。

圖8 木星在L4引力加速度變化頻譜圖

由頻譜圖可以看出，振幅出現(xiàn)峰值的3個(gè)頻率從小到大依次是

(18)

其對應(yīng)周期分別是

(19)

式中:T1與木星的公轉(zhuǎn)周期接近；T3與木星、地球的相遇周期接近，與將天體運(yùn)動看作勻速圓周運(yùn)動近似計(jì)算的兩行星相遇周期400 d幾乎一致。

根據(jù)直觀分析我們認(rèn)為所分析周期可以看作是三角函數(shù)形式的周期(將離散的點(diǎn)看作連續(xù)函數(shù)a=a(t)，所以假設(shè)函數(shù)是由3個(gè)余弦函數(shù)相加而成，即

a(t)=a0+∑Ricos(2πfit) i=1,2,3。

(20)

式中3個(gè)頻率是通過頻譜中的峰值頻率得到的。

利用最小二乘法擬合對點(diǎn)列x(n)做線性回歸，求出系數(shù)a=a0、Ri，其結(jié)果如圖9所示。

圖9 木星在L4引力加速度周期性變化與線性擬合函數(shù)對比圖

深色線表示的擬合曲線與淺色的原曲線相比圖像擬合十分理想，長短周期都基本吻合，在相位上也無較大誤差。因此，數(shù)值分析可以大致認(rèn)為攝動的周期性基本由拉格朗日點(diǎn)公轉(zhuǎn)周期和引起攝動的行星公轉(zhuǎn)周期決定。

2.2 L5點(diǎn)的受力分析圖示

L5點(diǎn)在受力方面與L4點(diǎn)的結(jié)果非常相似，互相驗(yàn)證。這里只將結(jié)果展示在圖10～圖14中以供參考。

圖10 行星在L5引力加速度與合加速度比較圖

圖11 金星在L5引力加速度與加速度比較圖

圖12 木星在L5引力加速度與合加速度比較圖

圖13 金星在L5引力加速度標(biāo)量變化圖

圖14 木星在L5引力加速度標(biāo)量變化圖

3 結(jié)果檢驗(yàn)

根據(jù)L4、L5點(diǎn)的分析結(jié)果，選取幾個(gè)特殊的點(diǎn)分析了極值點(diǎn)出現(xiàn)時(shí)有特別特征的天體金星與拉格朗日點(diǎn)的相對位置。圖15為1 050 d內(nèi)金星與三角拉格朗日點(diǎn)的軌道圖。

圖15 引力作用極值時(shí)位置示意圖

圖中深色線是金星的運(yùn)行軌道，淺色線為三角拉格朗日點(diǎn)的軌道(與地球軌道完全相同)；選取1 495 dL4點(diǎn)與金星的相對位置(深色圓點(diǎn))與717 dL5點(diǎn)與金星的相對位置(淺色圓點(diǎn))作為示例，2個(gè)時(shí)間點(diǎn)天體與拉格朗日點(diǎn)的位置是在周期內(nèi)達(dá)到最近：以此驗(yàn)證我們以上的分析是符合科學(xué)原理與常識的。

4 結(jié)束語

本文對N體問題的下拉格朗日點(diǎn)進(jìn)行定性和定量分析，從時(shí)間域、空間域、頻率域中找尋規(guī)律。通過研究發(fā)現(xiàn)，三角拉格朗日點(diǎn)在N體問題中(考慮8大行星)受到的攝動非常有限，與向心加速度相差在10個(gè)數(shù)量級左右。同時(shí)可以得出在黃道坐標(biāo)系內(nèi)，金星對拉格朗日點(diǎn)z方向上有一定周期性的擾動，與金星距離地球較近且有3°的軌道傾角情況相符。在黃道面內(nèi)拉格朗日點(diǎn)受木星擾動最大，且有一定的周期性規(guī)律，與木星質(zhì)量較大的情況相符。攝動加速度的周期性與行星公轉(zhuǎn)周期有密切關(guān)系，尤其是拉格朗日點(diǎn)與行星相遇周期有一定的復(fù)合規(guī)律。

致謝：感謝許國昌教授、山東大學(xué)空間科學(xué)研究院導(dǎo)航與遙感實(shí)驗(yàn)室杜玉軍等各位老師、師兄師姐和朋友們對本文的指導(dǎo)。

[1] XU G C，XU J.Orbits[M].Berlin：Springer Berlin，2013：103-120.

[2] CURTIS H D.Orbital mechanics for engineering students[M].3rd ed.Oxford：Butterworth-Heinemann Elsevier Ltd，2014.

[3] MURALIDHAR K，SARATHY R.A theoretical basis for perturbation methods[J].Statistics & Computing，2003，13(4)：329-335.

[4] 李廣宇.天體測量和天體力學(xué)基礎(chǔ)[M].北京：科學(xué)出版社，2015.

[5] 孫義燧.現(xiàn)代天體力學(xué)導(dǎo)論[M].高等教育出版社，2008.

[6] 李廣俠，姜勇，孫玉華.拉格朗日點(diǎn)在深空通信中的應(yīng)用[J].數(shù)字通信世界，2011(1)：84-87.

[7] 張寧博，喬棟.拉格朗日點(diǎn)觀測與探測任務(wù)發(fā)展現(xiàn)狀與分析[C]//中國宇航學(xué)會深空探測技術(shù)專業(yè)委員會.中國宇航學(xué)會深空探測技術(shù)專業(yè)委員會第八屆學(xué)術(shù)年會論文集(下篇).上海：中國宇航學(xué)會深空探測技術(shù)專業(yè)委員，2011：497-505.

[8] 林輝慶.拉格朗日L-4點(diǎn)的理論驗(yàn)算[J].物理教師，2012，33(4)：42-43.

[9] 胡先進(jìn)，李力.簡明推導(dǎo)“三角拉格朗日點(diǎn)”的新方法[J].物理教師，2013，34(5)：54-55.

[10]張文昭，劉文忠，孫艷春，等.N體問題中心構(gòu)型的數(shù)學(xué)定義[J].北京師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2011，47(4)：359-361.

[11]劉文中，張同杰.N體問題共線解的簡明數(shù)值方法[J].北京師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2005，41(1)：54-57.

[12]楊遠(yuǎn)玲，聶清香，吳曉梅，等.N體問題的幾種數(shù)值算法比較[J].計(jì)算物理，2006，23(5)：599-603.

Numerical study on Sun-Earth Lagrangian points of N-body problem

LIU Yuanqi1，ZHOU Shanyu1，GONG Zijia2

(1.School of Space Science and Physics/Institute of Space Science，Shandong University，Weihai，Shandong 264209，China；2.School of Mathematics and Statistics，Shandong University(Weihai)，Weihai，Shandong 264209，China )

Aiming at the problem that theN-body problem (N>3) has no analytical solution in the related existed study on Lagrangian points in Sun-Earth system，the paper tried to establish a 3D model to describe the solar system of gravitational fields based on the planet orbits，by calculating the precise gravitational value (the differences caused by theN-body gravity) of triangular Lagrangian pointsL4andL5using perturbation theory.Through the qualitative and quantitative analysis on the perturbation law of Lagrangian points，the conclusion about their sources，periodicity and so on could provide a reference for related research.

Lagrangian points；N-body problem；numerical study

2016-04-21

劉媛琪(1995—)，女，山東泰安人，大學(xué)本科生，主持國家級大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練項(xiàng)目，研究方向?yàn)橐岳窭嗜拯c(diǎn)相關(guān)問題為主的天體力學(xué)。

劉媛琪，周珊羽，鞏子嘉.N體問題下日地三角拉格朗日點(diǎn)的數(shù)值研究[J].導(dǎo)航定位學(xué)報(bào)，2016，4(4)：5-11.(LIU Yuanqi，ZHOU Shanyu，GONG Zijia.Numerical study on Sun-Earth Lagrangian points ofN-body problem[J].Journal of Navigation and Positioning，2016，4(4)：5-11.)

10.16547/j.cnki.10-1096.20160402.

P128.1

A

2095-4999(2016)04-0005-07