黃 凱，孫尚彪，楊永章，李祝蓮，李語強

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

月球天平動是月球自轉、公轉軌道以及周圍天體對月球力的作用的綜合表現?；诩す鉁y月數據制作的月球歷表是進行月球物理天平動研究的數據基礎[1]。

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

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

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

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

1 月球轉動歐拉角

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

圖1 INPOP19a與INPOP17a歐拉角的差別

圖2 INPOP17a與DE430歐拉角的差別

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

圖3 INPOP19a與DE430歐拉角的差別

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

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

2 月球天平動與月球轉動的力學模型

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

月球相對于地心的月面經度天平動和緯度天平動稱為光學天平動或幾何天平動。幾何天平動中的經度項和緯度項除了與月球本身的物理天平動有關，還與月球的公轉軌道密切相關。月球實際擺動形成的天平動稱為物理天平動，物理天平動是月球自身相對于空間坐標的轉動歐拉角，與月球的內部結構、自轉動量和轉動慣量有關。本文具體實現從歷表獲得的歐拉角到月球物理天平動的轉換過程。

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

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

(1)

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

(2)

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

(3)

(4)

(5)

此轉換過程可用公式

(6)

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

rdate=Rx(ε)BNP[ICRS]，

(7)

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

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

(8)

利用圖4與圖5中各量的關系，可以得到φ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)

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

φc=Ω+σ，

(14)

θc=I+ρ，

(15)

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

(16)

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

圖4 月球主軸坐標系到國際天球參考系中的旋轉歐拉角

圖5 月球主軸坐標系到黃道坐標系的歐拉角

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

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

3 物理天平動的頻譜分析

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

圖6 利用INPOP19a歷表實現對月球物理天平動的提取

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

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

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

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

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

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

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

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

4 總 結

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

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