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

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

1 引言

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

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

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

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

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

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

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

表1 錢德勒擺動周期和品質(zhì)因子的點估計和區(qū)間估計的現(xiàn)有結(jié)論①

在前人工作基礎(chǔ)上，本文創(chuàng)新性地提出了一種用于估計錢德勒擺動周期和品質(zhì)因子的數(shù)學(xué)模型和快速算法。新方法對模型誤差的假設(shè)更弱，估計結(jié)果的統(tǒng)計性質(zhì)更優(yōu)，新方法還采用了機制明確的優(yōu)化算法，計算效率更高。

2 極移理論

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

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

L是地球自轉(zhuǎn)角動量，在地固系下，其表達式是

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

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

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

3 估計算法

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

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

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

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

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

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

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

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

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

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

3.2 區(qū)間估計的方法——自助法

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

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

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

4 數(shù)據(jù)及預(yù)處理

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

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

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

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

圖1 預(yù)處理后的數(shù)據(jù)序列圖

5 結(jié)果

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

圖2 3 000組周期和品質(zhì)因子的參數(shù)估計值的統(tǒng)計直方圖

表2 3 000組周期和品質(zhì)因子的參數(shù)估計的均值、方差和90%置信區(qū)間

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

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

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

6 總結(jié)與討論

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

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

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

致謝

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