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

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

1 引言

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

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

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

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

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

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

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

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

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

2 極移理論

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

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

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

這里對(duì)式(7)做一點(diǎn)說明，傳統(tǒng)處理方法是先對(duì)式(5)化簡，再將其中的σcw替換為σc，即先化簡再替換，如文獻(xiàn)[31]中的式(13)和(14)，而本工作則是先替換再化簡。由于這兩種處理方法對(duì)m(t)的影響不超過1%[8]，并且觀測(cè)數(shù)據(jù)的測(cè)量精度有限，所以可以認(rèn)為這兩種做法一致。考慮極移m(t)與觀測(cè)極移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)或其他組織和機(jī)構(gòu)提供的極移觀測(cè)序列，其中y(t)以指向90°W為正。將式(8)代入到式(7)中，可得

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

3 估計(jì)算法

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

3.1 估計(jì)T c和Q c的快速算法

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

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

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

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

利用式(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)置)，可以獲得殘差平方和的解析表達(dá)式：

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

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

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

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

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

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

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

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

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

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

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

5 結(jié)果

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

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

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

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

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

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

6 總結(jié)與討論

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

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

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

致謝

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