N體問題下日地三角拉格朗日點的數值研究

劉媛琪，周珊羽，鞏子嘉

(1.山東大學(威海) 空間科學與物理學院/空間科學研究院，山東 威海 264209；2.山東大學(威海) 數學與統計學院，山東 威海 264209)

?

N體問題下日地三角拉格朗日點的數值研究

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

(1.山東大學(威海) 空間科學與物理學院/空間科學研究院，山東 威海 264209；2.山東大學(威海) 數學與統計學院，山東 威海 264209)

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

拉格朗日點；N體問題；數值方法

0 引言

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

從圖中可知，L4、L5點在距離太陽和地球相等距離的2個等邊三角形頂點處，所以也稱為三角拉格朗日點。眾所周知，三角形結構具有穩定性，科學家在1906年發現第617號小行星出現在木星軌道落后60°處；20世紀80年代，科學家又發現在土星和它的大衛星所構成的運動系統中也有類似的等邊三角形。數學的推導和天文觀測結論證明了拉格朗日點的特殊性確實存在。

圖1 日地三體系統的拉格朗日點

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

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

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

本文利用開普勒根數計算給定時間黃道坐標系內8大行星的精確位置。假設太陽和地球對拉格朗日點的引力和與其向心力完全相同，其他行星在拉格朗日點的加速度通過最基礎的加速度矢量疊加，即可求解拉格朗日點受到影響(數值)的規律。

通過研究發現，在黃道坐標系內，金星對拉格朗日點z方向上有一定周期性的擾動，與金星距離地球較近且有3°的軌道傾角情況相符。同時，在黃道面內，拉格朗日點受木星擾動最大，且有一定的周期性規律，與木星距離地球近且質量較大的情況相符。

通過攝動理論，可以解決許多天文學與空間科學中的重要問題，例如證明太陽系存在的第九大行星、深空探測衛星軌道設計等。如能夠通過數值方法將攝動理論充分應用，將對攝動理論的發展產生重要的影響。

1 計算過程

1.1 確定8大行星在給定日期的狀態矢量

其步驟如下：

1)求儒略日JD[2]203-208。

根據儒略日的定義，用下式推算某日世界時0時刻的儒略日

(1)

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

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

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

(2)

3)求得軌道根數在T0時的值

(3)

4)求JD時的近日點幅角ω、真近點角θ[2]155-161。

ω=?-Ω

M=L-?，

(4)

式中：?為近日點緯度；Ω為升交點赤經；M為平近點角。

開普勒方程為

M=E-e sinE。

(5)

式中：E為偏近點角；e為橢圓軌道偏心率。

可以通過迭代逼近求得偏近點角E為

(6)

即可求得真近點角θ。

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

①軌道面坐標系下(以行星各自軌道面為z軸，軌道焦點為原點)有

(7)

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

②轉換為日心黃道坐標系的變換矩陣為

(8)

式中：Ω為升交點赤經；ω為近日點幅角；i為軌道傾角。則

(9)

1.2 L4、L5拉格朗日點的精確軌道

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

(10)

式中：R代表日心黃道坐標系內位置參量；下標L4代表拉格朗日4點上測試質量的位置參量；下標3代表地球。由于L4點的軌道也在地球公轉平面上，因此在黃道坐標系的xOy平面上可通過地球狀態參量旋轉得到其狀態參量。

2)拉格朗日點L5：L5點在以2個主天體連線為底邊的等邊三角形的第3個頂點上，在較小天體圍繞2天體系統質心運行軌道的后方。則有

(11)

式中下標L5代表拉格朗日5點上測試質量的位置參量。由于L5點的軌道也在地球公轉平面上，因此在黃道坐標系的xOy平面上可通過地球狀態參量旋轉得到其狀態參量。

1.3 計算拉格朗日點受其他行星運行的影響

1)目標函數(向量)

(12)

式中：F為拉格朗日點軌道上實驗質量m受到的萬有引力合力矢量；下標N-body代表N體問題條件下，3-body代表三體問題條件下。

2)三體問題下

(13)

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

3)N體問題下

(14)

4)目標函數

(15)

目標函數是N體問題和三體問題萬有引力合力的差值，也就是三體問題下拉格朗日點受到行星的攝動力影響。

2 數值分析

2.1 L4點的受力分析

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

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

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

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

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

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

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

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

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

圖6 金星在L4引力加速度標量變化圖

圖7 木星在L4引力加速度標量變化圖

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

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

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

(16)

設第k個點對應頻率是fk，則

(17)

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

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

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

由頻譜圖可以看出，振幅出現峰值的3個頻率從小到大依次是

(18)

其對應周期分別是

(19)

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

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

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

(20)

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

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

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

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

2.2 L5點的受力分析圖示

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

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

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

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

圖13 金星在L5引力加速度標量變化圖

圖14 木星在L5引力加速度標量變化圖

3 結果檢驗

根據L4、L5點的分析結果，選取幾個特殊的點分析了極值點出現時有特別特征的天體金星與拉格朗日點的相對位置。圖15為1 050 d內金星與三角拉格朗日點的軌道圖。

圖15 引力作用極值時位置示意圖

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

4 結束語

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

致謝：感謝許國昌教授、山東大學空間科學研究院導航與遙感實驗室杜玉軍等各位老師、師兄師姐和朋友們對本文的指導。

[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] 李廣宇.天體測量和天體力學基礎[M].北京：科學出版社，2015.

[5] 孫義燧.現代天體力學導論[M].高等教育出版社，2008.

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

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

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

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

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

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

[12]楊遠玲，聶清香，吳曉梅，等.N體問題的幾種數值算法比較[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—)，女，山東泰安人，大學本科生，主持國家級大學生創新創業訓練項目，研究方向為以拉格朗日點相關問題為主的天體力學。

劉媛琪，周珊羽，鞏子嘉.N體問題下日地三角拉格朗日點的數值研究[J].導航定位學報，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