基于數(shù)值歷表的月球物理天平動研究

黃 凱，孫尚彪，楊永章，李祝蓮，李語強(qiáng)

(1. 中國科學(xué)院云南天文臺，云南 昆明 650216；2. 中國科學(xué)院大學(xué)，北京 100049；3. 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，吉林 長春 130026)

月球天平動是月球自轉(zhuǎn)、公轉(zhuǎn)軌道以及周圍天體對月球力的作用的綜合表現(xiàn)?；诩す鉁y月數(shù)據(jù)制作的月球歷表是進(jìn)行月球物理天平動研究的數(shù)據(jù)基礎(chǔ)[1]。

根據(jù)卡西尼定則，月球的公轉(zhuǎn)周期與自轉(zhuǎn)周期一致，但由于月球質(zhì)量的非球面分布，月球的自轉(zhuǎn)并不是均勻的。這種與均勻旋轉(zhuǎn)的偏離稱為物理天平動。物理天平動有兩種類型，受迫天平動和自由天平動。受迫天平動是由于地球、太陽和行星的引力引起的月球形狀上的時變力矩；自由天平動主要是地質(zhì)活動引起的，如沖擊，核-幔相互作用，或是與受迫天平動的共振[2]。理論上，受迫天平動的振幅、相位和周期可以根據(jù)月球形狀特征、彈性形變和旋轉(zhuǎn)耗散計(jì)算。對于自由天平動，只能計(jì)算周期，振幅和相位必須通過測量獲得[3]。

自由天平動有3種模式，一種在經(jīng)度，兩種在極位。經(jīng)度模式是圍繞月球極軸的周期為2.9年的擺狀振蕩[3]。極位模式描述了月球極位與周期為18.6年的受迫進(jìn)動以及受迫天平動所決定的極位置之間的變化。第一個極位模式稱為擺動模式，類似于地球的錢德勒擺動；另一個極位模式是自由進(jìn)動模式，對應(yīng)于空間中極軸周期為81年的運(yùn)動。隨著時間的推移，物理天平動由于能量耗散而逐漸減弱[4]。

月球物理天平動可以通過解析法和數(shù)值法兩種方式計(jì)算。早期因?yàn)橛^測數(shù)據(jù)缺乏以及數(shù)值算法的計(jì)算量大，數(shù)值法相對于解析法并沒有太大的優(yōu)勢，在工程上解析法的攝動理論完全可以達(dá)到精度要求[5]。但是隨著深空探測的不斷發(fā)展，對精度的要求逐漸提高，由于計(jì)算公式復(fù)雜，解析法開始顯出劣勢[6]。與此同時，在更多觀測數(shù)據(jù)和計(jì)算機(jī)技術(shù)的支持下，數(shù)值法能達(dá)到更高的精度。隨著地月數(shù)值歷表的建立，數(shù)值法逐漸成為月球天平動研究的主要方法。月球歷表包含月球在國際天球坐標(biāo)系中位置和速度狀態(tài)的信息[7-8]，本文對月球軌道及物理天平動的研究所使用的數(shù)據(jù)均來自地月歷表。

目前，我國月球探測已經(jīng)進(jìn)入月面著陸的快速發(fā)展期，特別是云南天文臺和中山大學(xué)相繼成功開展了月球激光測距試驗(yàn)，為我們加入國際月球激光測距分析社區(qū)提供了技術(shù)基礎(chǔ)[8-9]。隨著地面測月站數(shù)量的不斷增加，國內(nèi)激光測月向更高精度發(fā)展，未來會積累大量的高精度數(shù)據(jù)，這些數(shù)據(jù)為歷表的研制和物理天平動的研究提供支撐。進(jìn)入新世紀(jì)后，我國提出了月球和深空探測計(jì)劃，并成功實(shí)施了月球探測項(xiàng)目。文[4]基于快速傅里葉變換，提出了一種頻率提取和擬合的數(shù)值方法，從DE430歷表提取月球自由天平動周期頻率，與國際同行結(jié)果一致。在動力學(xué)模型方面，文[5-6]建立了基于國際上最新動力學(xué)模型的月球自轉(zhuǎn)雙層/單層模型及數(shù)值算法，結(jié)果精度與主流歷表相當(dāng)。本文使用最新的INPOP19a歷表實(shí)現(xiàn)對月球物理天平動數(shù)據(jù)的提取。月球物理天平動數(shù)值模型精度的提高，為進(jìn)一步高精度月球激光測距的研究打下基礎(chǔ)。

1 月球轉(zhuǎn)動歐拉角

月球歷表中包含各個天體在國際天球坐標(biāo)系中的位置和速度狀態(tài)信息，還包含月球的轉(zhuǎn)動歐拉角。利用INPOP19a，INPOP17a和DE430歷表中跨度600年的數(shù)據(jù)提取得到的歐拉角比較如圖1～圖3。

圖1 INPOP19a與INPOP17a歐拉角的差別

圖2 INPOP17a與DE430歐拉角的差別

由圖1～圖3可知，INPOP17a與DE430歐拉角差別較大，INPOP19a更為穩(wěn)定。其中，φ和θ在不同歷表中相差0.05″左右。但是ψ隨著時間逐漸增大，說明歷表中ψ的擬合模型不穩(wěn)定。對比歷表歐拉角的差別，計(jì)算的地心到月面反射器A15的距離最大誤差有30 cm，對月球激光測距的預(yù)報(bào)精度有較大的影響。

圖3 INPOP19a與DE430歐拉角的差別

以上完成了對天平動歐拉角的提取，利用快速傅里葉變換算法對提取的歐拉角之差進(jìn)行了頻譜分析，結(jié)果如表1。27.3天是極位置在慣性坐標(biāo)系下自由天平動的周期。

表1 對歐拉角進(jìn)行頻譜分析得到的周期(單位： 天)

2 月球天平動與月球轉(zhuǎn)動的力學(xué)模型

月球天平動是指地面觀測者所觀測的月球可見面上下左右小幅度的擺動。由于天平動效應(yīng)，觀測者能看到69%的月面，其中有19%時多時少。

月球相對于地心的月面經(jīng)度天平動和緯度天平動稱為光學(xué)天平動或幾何天平動。幾何天平動中的經(jīng)度項(xiàng)和緯度項(xiàng)除了與月球本身的物理天平動有關(guān)，還與月球的公轉(zhuǎn)軌道密切相關(guān)。月球?qū)嶋H擺動形成的天平動稱為物理天平動，物理天平動是月球自身相對于空間坐標(biāo)的轉(zhuǎn)動歐拉角，與月球的內(nèi)部結(jié)構(gòu)、自轉(zhuǎn)動量和轉(zhuǎn)動慣量有關(guān)。本文具體實(shí)現(xiàn)從歷表獲得的歐拉角到月球物理天平動的轉(zhuǎn)換過程。

月固坐標(biāo)系分為月球主軸坐標(biāo)系(Principal Axis Frame, PA)與平均地球指向坐標(biāo)系(Mean Earth, ME)。在描述月球轉(zhuǎn)動時，我們以月球質(zhì)心為原點(diǎn)，月球3個慣量主軸為3個坐標(biāo)軸，建立隨月球轉(zhuǎn)動的月球主軸坐標(biāo)系。同時，我們可以得到坐標(biāo)指向與國際天球參考系(International Celestial Reference System, ICRS)一致的月心天球參考系。兩個坐標(biāo)系通過3個歐拉角φ，θ和ψ進(jìn)行轉(zhuǎn)換：

[ICRS]=Rz(-φ)·Rx(-θ)·Rz(-ψ)·[PA].

(1)

這3個與月球轉(zhuǎn)動相關(guān)的歐拉角及變化率即月球天平動的數(shù)值表達(dá)，反映月球自轉(zhuǎn)的狀態(tài)。依據(jù)卡西尼定則：(1)月球自轉(zhuǎn)周期等于公轉(zhuǎn)周期；(2)月球赤道與黃道的夾角為1°32′；(3)月球自轉(zhuǎn)軸、月球軌道平面的法線以及黃道面的法線三者共面，且第三者處于前兩者之間。月球物理天平動的理論基礎(chǔ)主要是剛體轉(zhuǎn)動微分方程，由卡西尼第一定律(月球自轉(zhuǎn)周期等于公轉(zhuǎn)周期)簡化得到，對于剛體月球轉(zhuǎn)動微分方程為經(jīng)典的歐拉-劉維爾方程

(2)

其中，ω代表月球的轉(zhuǎn)動角速度；T為月球受到的攝動力。我們可以直接從DE歷表讀取月心天球參考系到月球主軸坐標(biāo)系的轉(zhuǎn)換矩陣。歐拉角從慣性參考系轉(zhuǎn)換到旋轉(zhuǎn)坐標(biāo)系時，

(3)

(4)

(5)

此轉(zhuǎn)換過程可用公式

(6)

表示。(6)式是月心天球參考系到月球主軸坐標(biāo)系轉(zhuǎn)換的理論基礎(chǔ)。將這一向量轉(zhuǎn)換到瞬時黃道坐標(biāo)系，有

rdate=Rx(ε)BNP[ICRS]，

(7)

其中，B,N和P分別為坐標(biāo)偏差矩陣、歲差矩陣和章動矩陣；ε為瞬時真黃道坐標(biāo)系的旋轉(zhuǎn)。

將月球主軸坐標(biāo)系中的單位向量R1=(1, 0, 0)和R2=(0, 0, 1)通過上述坐標(biāo)轉(zhuǎn)換至瞬時黃道坐標(biāo)系，得到Xdate和Zdate。i,j和k分別為在黃道坐標(biāo)系下沿著Ox，Oy和Oz的單位向量。n為坐標(biāo)原點(diǎn)指向月球赤道與黃道的升交點(diǎn)的單位向量，此時有

(8)

利用圖4與圖5中各量的關(guān)系，可以得到φc,θc和ψc：

cosφc=i·n，

(9)

sinφc=j·n，

(10)

cosθc=k·zdate，

(11)

cosψc=n·zdate，

(12)

sinψc=(zdate×n)·xdate.

(13)

根據(jù)Newhall[1]研究中使用的公式，我們可以得到

φc=Ω+σ，

(14)

θc=I+ρ，

(15)

ψc=τ-σ+F+180°.

(16)

在上述過程中，我們把圖4赤道坐標(biāo)系下的歐拉角ψ，θ和φ轉(zhuǎn)換為圖5黃道坐標(biāo)系下的φc，θc和ψc。Eckhardt分別使用Ισ，ρ和τ表示天平動中的周期項(xiàng)。其中，ρ是由θc分離出來，定義為解的周期部分，θc由(15)式定義，I是一個線性多項(xiàng)式。根據(jù)(14)式，Ω是軌道節(jié)點(diǎn)的經(jīng)度，σ是兩個節(jié)點(diǎn)之間的瞬時差，將(14)式變?yōu)棣│誧=ΙΩ+Ισ，Ισ定義為由一個常數(shù)和解的周期部分組成。同時(16)式可以變化為ψc+φc=τ+Ω+F+180°，定義τ由一個常數(shù)和結(jié)果的周期部分組成。我們分別使用Ισ，τ和ρ表示月球的物理天平動，根據(jù)上述過程計(jì)算得到結(jié)果最終如圖6。

圖4 月球主軸坐標(biāo)系到國際天球參考系中的旋轉(zhuǎn)歐拉角

圖5 月球主軸坐標(biāo)系到黃道坐標(biāo)系的歐拉角

至此完成了利用數(shù)值歷表中歐拉角對物理天平動的轉(zhuǎn)換。利用上述方法，從INPOP17a與DE430中分別提取物理天平動與INPOP19a作差比較，結(jié)果如圖7和圖8。結(jié)果顯示，INPOP19a與INPOP17a物理天平動差別更為穩(wěn)定。

圖8 INPOP19a與DE430物理天平動差別

3 物理天平動的頻譜分析

由圖6物理量τ繪制的圖與圖7中INPOP19a與INPOP17a之差對比，我們發(fā)現(xiàn)存在一個穩(wěn)定頻率的周期，且此頻率與τ本身的固有頻率不同，故在研究過程中重點(diǎn)對此頻率進(jìn)行分析。

圖6 利用INPOP19a歷表實(shí)現(xiàn)對月球物理天平動的提取

圖7 INPOP19a與INPOP17a物理天平動差別

首先，根據(jù)上述過程，利用INPOP19a，INPOP17a和DE430歷表中跨度600年的數(shù)據(jù)提取得到的物理天平動比較如圖9～圖11。

圖9 INPOP17a與DE430物理天平動的差別

圖10 INPOP19a與DE430物理天平動的差別

圖11 INPOP19a與INPOP17a物理天平動的差別

我們利用快速傅里葉變換算法對上面的圖像進(jìn)行頻譜分析，結(jié)果如表2。

表2 對物理天平動進(jìn)行頻譜分析得到的周期(單位： 天)

其中，1 056.53天的周期是經(jīng)度方向的自由天平動周期，27.21天和27.18天則是物理天平動中F模型和Wobble模型引起的。

4 總 結(jié)

本文利用已有的數(shù)值歷表數(shù)據(jù)開展月球物理天平動研究，得到INPOP19a歷表中3個歐拉角的相位時間圖像，與INPOP17a，DE430歷表進(jìn)行對比，發(fā)現(xiàn)其差值存在27.3天的周期，識別出了極位置在慣性坐標(biāo)系下的自由天平動。本文進(jìn)一步利用歷表數(shù)據(jù)對月球天平動進(jìn)行分析，根據(jù)卡西尼定則以及國際天球參考系到月球主軸坐標(biāo)系的轉(zhuǎn)換模型，提取得到月球相對于空間的擺動即月球的物理天平動。結(jié)果顯示，月球的物理天平動存在一定的周期，但振幅在200毫角秒內(nèi)，與我們能看到69%的月面相比，是一個極其微小的擺動量。即在地球上觀測月球的天平動時，月球的幾何天平動為主要影響因素。3個歷表提取得到的物理天平動作差比較，發(fā)現(xiàn)存在一個穩(wěn)定周期，且此周期與固有周期不同，頻譜分析以后將3個不同歷表得到的周期進(jìn)行了比較。

本文對提取得到的歐拉角進(jìn)行分析，計(jì)算得到新一代數(shù)值歷表INPOP19a與DE430，INPOP17a提取得到的歐拉角換算到月球激光測距上最多有30 cm的誤差。此誤差對月球激光測距的預(yù)報(bào)精度產(chǎn)生較大的影響。與DE430，INPOP17a相比，INPOP19a有較高的穩(wěn)定性，故在月球激光測距或者研究物理天平動時，推薦采用歷表INPOP19a。