檢測錢德勒擺動周期和Q值的一種新方法

胥燦燦，方 明，許雪晴，周永宏,4，段鵬碩

(1.中國科學院 上海天文臺，上海200030；2.麻省理工學院 地球、大氣和行星科學系，美國02139；3.中國科學院大學 天文與空間科學學院，北京100049；4.中國科學院 行星科學重點實驗室，上海200030)

1 引言

錢德勒擺動是地球自轉的一個本征模。自1891年首次被發現以來，引起錢德勒擺動的物理機制和描述它的兩個關鍵參數(周期和品質因子)一直是科學界關注的熱點問題?？臻g大地測量技術為檢測錢德勒擺動提供了豐富、高精度和高時空分辨率的觀測資料，人們利用這些高精度的觀測資料，在引起錢德勒擺動的物理機制有關領域的研究取得了很大進展，例如Smith和Dahlen[1]、朱耀仲[2?4]、Mathews等人[5]，以及Chen等人[6,7]在理論上闡述了各種地球物理因素(如海洋、地幔粘彈性、核幔耦合等)對錢德勒擺動的貢獻；Wahr[8]、Eubanks[9]、Barnes等人[10]、Gross[11,12]、Liao等人[13]，以及Fang等人[14]研究了大氣、海洋等地表流體激發錢德勒擺動的物理機制。盡管如此，由于理論假設的局限性[6]、部分觀測數據的不確定性[15]，以及觀測數據的不完備性[8]等原因，到目前為止，引起錢德勒擺動的物理機制，擺動周期和品質因子的取值也一直沒有明確的結論[16]。

由于錢德勒擺動與很多地球物理因素有關，因此，對其精確的檢測具有很重要的科學價值。目前，關于錢德勒擺動周期和品質因子的計算和估計方法主要分為如下四類：

(1)譜分析方法。直接對極移觀測序列進行譜分析，從而獲得錢德勒擺動周期和品質因子[17?19]。

(2)隨機激發假設方法。假設錢德勒擺動的激發機制是隨機的，采用統計方法來估計錢德勒擺動周期和品質因子[20?24]。

(3)地球物理激發理論方法。采用實際觀測數據，研究大氣、海洋以及陸地水等地球物理機制對錢德勒擺動的影響[16,25?28]。

(4)半解析方法?；谝欢ǖ睦碚撃Ｐ秃蛯嶋H觀測資料，例如海洋、地幔、核幔耦合和三軸橢球等等，給出錢德勒擺動周期和品質因子的解析表達式[1,2?5,29]。

本文將前三類方法獲得的周期和品質因子稱之為觀測值，將最后一類方法獲得的稱之為理論值。表1統計了現有方法對錢德勒擺動周期和品質因子的計算和估計結果。從表1可以看出，品質因子的參數估值具有極大的取值范圍(36,1 000)，盡管最新研究已經將其置信區間縮小至(56,255)[16]，但仍然比較大。

表1 錢德勒擺動周期和品質因子的點估計和區間估計的現有結論①

在前人工作基礎上，本文創新性地提出了一種用于估計錢德勒擺動周期和品質因子的數學模型和快速算法。新方法對模型誤差的假設更弱，估計結果的統計性質更優，新方法還采用了機制明確的優化算法，計算效率更高。

2 極移理論

假設地球不受外部天體的引力作用，在地固坐標系下，角動量定理可以寫成劉維爾方程，即：

其中，ω是地固系相對慣性系的轉動角速度。用Ω表示地球自轉的平均角速度，則ω可以進一步表示為：

L是地球自轉角動量，在地固系下，其表達式是

這里對式(7)做一點說明，傳統處理方法是先對式(5)化簡，再將其中的σcw替換為σc，即先化簡再替換，如文獻[31]中的式(13)和(14)，而本工作則是先替換再化簡。由于這兩種處理方法對m(t)的影響不超過1%[8]，并且觀測數據的測量精度有限，所以可以認為這兩種做法一致?？紤]極移m(t)與觀測極移p(t)的關系，可得[31]：

其中，p(t)=x(t)-iy(t)，x(t)和y(t)就是由國際地球自轉服務(International Earth Rotation and Reference Systems Service,IERS)或其他組織和機構提供的極移觀測序列，其中y(t)以指向90°W為正。將式(8)代入到式(7)中，可得

對這個微分方程進行求解，其中p(t)的自由項在經過很長一段時間后逐漸衰減消失，只剩下受迫激發項，因此p(t)的解析表達式為：

3 估計算法

基于式(10)，本章對誤差進行基本假設，提出用于估計Tc和Qc的數學模型；同時，基于這些誤差假設和新的數學模型，引入Tc和Qc的區間估計方法——自助法。

3.1 估計T c和Q c的快速算法

假設數據是等間隔采樣的，采樣間隔為δt，數據長度為N，初始歷元為t0。根據式(10)可以給出離散數據滿足如下一階自回歸形式：

考慮到Π(t)不完全是大氣和海洋的角動量，可能還包含了不可觀測的和未知的因素，如在式(6)的h(t)和ΔI(t)中，它們不僅包含了地表流體的作用，還包括了地球內部的作用(如地磁急變[19])，而后者不可直接觀測[8]。所以，可以將Π(t)分解為如下三個部分：

其中，Πobv(t)+Πerr(t)是可觀測部分的真值，Πobv(t)是與之對應的實際觀測數據，Πerr(t)是觀測數據的測量誤差(即觀測誤差)；Πres(t)是上面提到的其他不可觀測的或未知的激發源。這里將式(12)的后兩項合并，記為Πmer(t)，簡稱模型誤差。最后，假設模型誤差對應的離散序列Πmer(tn)是一個均值函數為0、一階自相關函數為0的隨機噪聲。

將Π(t)=Πobv(t)+Πmer(t)代入式(11)中，對關于Πobv(t)的積分用辛普森公式展開，而對關于Πmer(t)的積分用梯形公式展開，并引入以下表達式：

利用式(13)和(14)將式(11)展開并化簡為：

其中，n=0,1,2,···,N-2。引入矢量Y,X1,X2,ε和?ε，式(15)可以被改寫為回歸方程的形式：

考慮到序列Πmer(tn)的一階自相關函數為零，因此εH?ε=0成立(上標H代表矢量的共軛轉置)，可以獲得殘差平方和的解析表達式：

總之，該方法采用了一種機制明確的算法來直接獲得最優的?Tc和?Qc，而不用在某個特定區域內比較每個網格點上的結果，例如Nastula和Gross[16]采取的方法(下面將其簡稱為網格搜索法)，因此本方法具有更高的計算效率。

3.2 區間估計的方法——自助法

目前關于錢德勒擺動周期和品質因子的區間估計都是基于蒙特卡洛方法開展的[16,25]，這種方法需要事先獲取模型誤差Πmer(t)的理論分布。根據式(12)可知：Πres(t)是未知的或不可觀測的，假設它服從任何一種分布可能都不合理；此外，隨著觀測水平的提高，觀測數據的測量精度也在逐漸提高，因此觀測誤差Πerr(t)具有異方差性。綜合這兩個因素可以發現，理論模型難以描述Πmer(t)的分布函數，因此蒙特卡洛方法的區間估計結果可能不夠精確。

基于此，本工作采用自助法(Bootstrapping)做區間估計，該方法的優勢是不需要知道Πmer(t)的分布函數。自助法是一種現代的非參數統計方法，最早由Efron[32]提出，經過近幾十年的發展，它已經具備了扎實的理論背景并擁有廣泛的應用范圍[33]。自助法的操作流程簡捷，具體為：(1)將包含K個樣品的初始數據集當做總體，每次從中有放回地抽取K個樣品，如此重復B次，就得到了B個新樣本；(2)對每個再抽樣的樣本分別做參數估計，可以得到每個參數的B個估值；(3)基于這B組估值，可以給出基于初始數據集的區間估計，甚至研究參數估值的有效性。

本工作采用自助法來對Tc和Qc進行區間估計，按照以上步驟將式(15)中的N-1個方程當成總體(這一做法在參考文獻[33]中有詳細的描述)，每次從中有放回地抽取N-1個樣品，然后根據式(16)―(21)對再抽樣的樣本做參數估計，如此重復多次，可以獲取多組Tc和Qc的參數估值，據此給出Tc和Qc的區間估計結果。

4 數據及預處理

對于極移觀測序列，我們采用Ratcliff和Gross[34]解算的COMB2018數據集，該數據集在1990年后的時間段具有非常高的精度。大氣角動量(atmospheric angular momentum,AAM)數據來自美國國家環境預報中心/大氣研究中心(National Centers for Environmental Prediction/National Center for Atmospheric Research,NCEP/NCAR)數據集[35]，該數據集考慮了地形因素的影響，精度更高[36]。海洋角動量(oceanic angular momentum,OAM)數據來自IERS特殊海洋局(Special Bureau for Ocean)的ECCO-kf080i數據集。這三個數據集的采樣間隔都是1 d，選取的公共時間跨度為1993/01―2018/12。在參數估計之前，對這些數據做如下預處理[16,25]：

(1)扣除三個序列中的周年項、季節項和線性趨勢項。采用最小二乘法從原始序列中扣除頻率為1,2,3,4,5,6,7,8 cpy的諧波信號，兩個周期分別為13.661 d和13.633 d的潮汐項，一個一階多項式函數。

(2)扣除低頻項。除線性趨勢項外，極移序列中的低頻信號一般被認為來自地球內部，而這些激發源不能被直接觀測，據此設計一個高通濾波器，過濾掉三個序列中周期大于2 a的頻段。

預處理后的數據展示在圖1中，從該圖可以看出：錢德勒擺動的振幅一直處于調制中，總體在減小，而在2010―2019年期間，錢德勒擺動的振幅已經減小到非常顯著的程度(這也可以在Wang等人[37]的圖8中得到驗證)；對于預處理后的角動量序列，與平穩隨機噪聲非常相似，并且實部振幅比虛部振幅稍小。

圖1 預處理后的數據序列圖

5 結果

為了更全面地評估新方法的穩定性和可靠性，本工作選用三個不同的時間段對研究結果進行分析和討論，分別是1993―2010年、2000―2019年和1993―2019年。在每個時間段都使用自助法做1 000次再抽樣，每個時間段都得到1 000組?Tc和?Qc。這3 000組估計值的統計直方圖如圖2所示，它們的均值、方差和90%置信區間如表2所示。

圖2 3 000組周期和品質因子的參數估計值的統計直方圖

表2 3 000組周期和品質因子的參數估計的均值、方差和90%置信區間

從圖2和表2的統計結果可以看出，第一個時間段(1993―2010年)和第三個時間段(1993―2019年)給出的Tc和Qc的區間估計幾乎是重合的，并且?Tc的方差很??；對于Qc的估計，其對數據的選擇比較敏感，這也是現有估計方法共同的問題(如表1所示)。對比表1和表2可以看出，新方法對Qc的估計具有更小的方差，顯示了更優的統計有效性。

從圖2和表2的統計結果還可以發現，第二個時間段(2000―2019年)與其他兩個時間段的參數估計結果差距較大。這可能與近年來錢德勒擺動的振幅逐漸減小有關，該時間段的估計結果可能具有一定的系統性偏差，我們將在未來的研究中進一步探討。

綜合以上結果，可以確定第一個時間段(1993―2010年)的參數估計更加合理，即Tc和Qc的點估計為430.8 d和62.6，兩者的90%置信區間分別為(430.0,431.6)d和(43.5,75.7)。在該時間段內的參數置信區間更寬，具有更大的概率包含Qc的真值，是一個比較保守的估計。此外，這一結果與Mathews等人[5]及Nastula和Gross[16]的結果很接近，再一次驗證了新方法的可靠性。

6 總結與討論

本文從經典的極移理論出發，給出了一種用于估計錢德勒擺動周期和品質因子的新方法。新方法對于誤差的假設較弱，并采用更加合理的自助法進行區間估計，參數估計的結果具有更優的統計性質。此外，新的點估計算法采用了機制明確的雙點割線法來直接獲取最優解，與傳統的網格搜索法相比具有更高的計算效率。

對地球錢德勒擺動周期和品質因子的最優點估計結果分別為：430.8 d和62.6，兩者的標準差為0.50 d和9.63，最優的區間估計分別為(430.0,431.6)d和(43.5,75.7)。這一結果與Mathews等人[5]及Nastula和Gross[16]的結果很接近，且我們的參數估值具有更小的方差，再次驗證了新方法的可靠性和優越性。

由于目前只采用了大氣和海洋的角動量數據，對地球錢德勒擺動周期和品質因子的最優估計可能還存在一定偏差。未來我們將考慮陸地水或全球水文角動量的影響，進一步提高錢德勒擺動周期和品質因子的估計精度。此外，隨著對錢德勒擺動認識的逐步提高以及觀測數據的精度提升，本工作提出的新方法將對地球錢德勒擺動的精準預報提供可能。

致謝

感謝兩位審稿人對本文的肯定，他們提出的意見和建議也十分寶貴，這大大提高了文章的質量和可讀性。與上海天文臺廖新浩老師的交流討論讓我們受益匪淺，感謝廖老師在這個過程中提出的建設性意見，這促使新方法改進得更加完善。