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

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

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

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

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

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

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

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

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

圖1 INPOP19a與INPOP17a歐拉角的差別

圖2 INPOP17a與DE430歐拉角的差別

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

圖3 INPOP19a與DE430歐拉角的差別

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

(4)

(5)

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

(6)

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

rdate=Rx(ε)BNP[ICRS]，

(7)

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

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

(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)

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

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

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

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

圖8 INPOP19a與DE430物理天平動(dòng)差別

3 物理天平動(dòng)的頻譜分析

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

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

圖7 INPOP19a與INPOP17a物理天平動(dòng)差別

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

圖9 INPOP17a與DE430物理天平動(dòng)的差別

圖10 INPOP19a與DE430物理天平動(dòng)的差別

圖11 INPOP19a與INPOP17a物理天平動(dòng)的差別

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

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

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

4 總 結(jié)

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

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