基于有效角動(dòng)量函數(shù)的多參數(shù)極移預(yù)報(bào)*

趙 鑫 吳元偉 楊新宇 楊旭海 張首剛

(1中國(guó)科學(xué)院國(guó)家授時(shí)中心 西安 710600)

(2中國(guó)科學(xué)院大學(xué)電子電氣與通信工程學(xué)院 北京 100049)

(3中國(guó)科學(xué)院精密導(dǎo)航定位與定時(shí)技術(shù)重點(diǎn)實(shí)驗(yàn)室 西安 710600)

(4中國(guó)科學(xué)院時(shí)間頻率基準(zhǔn)重點(diǎn)實(shí)驗(yàn)室 西安 710600)

1 引言

地球在空間中的受力非常復(fù)雜，外部受到太陽(yáng)、月球及行星的萬(wàn)有引力作用，內(nèi)部受到大氣、海洋、地殼、地幔等圈層相互作用的影響[1]，導(dǎo)致了地球的自轉(zhuǎn)不恒定.一方面，自轉(zhuǎn)速率隨著時(shí)間變化;另一方面，自轉(zhuǎn)軸的方位也隨著時(shí)間變化.通常，地球的方位可以用地球定向參數(shù)(Earth Orientation Parameters，EOP)來(lái)描述，包括章動(dòng)參數(shù)、世界時(shí)(Universal Time，UT1)和極移參數(shù)(Polar Motion，PM)，其中，極移參數(shù)包括極移X和極移Y兩個(gè)分量，這些參數(shù)與國(guó)際天球參考架(International Celestial Reference Frame，ICRF)和國(guó)際地球參考架(International Terrestrial Reference Frame，ITRF)之間的轉(zhuǎn)換密切相關(guān)[2].地球定向參數(shù)對(duì)于許多科學(xué)研究和工程應(yīng)用具有重要的作用，例如，在大地測(cè)量學(xué)、天文學(xué)、航海以及時(shí)間保持等領(lǐng)域，不僅需要實(shí)時(shí)地球定向參數(shù)，也需要高精度的地球定向參數(shù)預(yù)報(bào)[3].國(guó)際地球自轉(zhuǎn)服務(wù)(International Earth Rotation Service，IERS)定期提供并維護(hù)地球定向參數(shù)的事后產(chǎn)品、快速產(chǎn)品和預(yù)報(bào)產(chǎn)品.

極移是地球自轉(zhuǎn)軸相對(duì)于地球表面的運(yùn)動(dòng).作為地球定向參數(shù)的組成部分，極移的預(yù)報(bào)同樣對(duì)于科學(xué)研究和工程應(yīng)用具有重要的意義[4].近些年來(lái)，很多學(xué)者參與了極移預(yù)報(bào)的研究，并提出了多種極移預(yù)報(bào)的方法.在這些預(yù)報(bào)方法中，絕大部分都是對(duì)事后極移數(shù)據(jù)的統(tǒng)計(jì)外推[5]，其中，最小二乘與自回歸組合(Least-Squares Extrapolation and Autoregressive Modeling，LS+AR)的方法很早就被應(yīng)用于極移的預(yù)報(bào)[6];于2005-2008年舉行的地球定向參數(shù)預(yù)報(bào)比賽(EOP Prediction Comparison Campaign)對(duì)極移預(yù)報(bào)方法的總結(jié)以及預(yù)報(bào)精度的提升具有里程碑的意義，除了最小二乘與自回歸組合的預(yù)報(bào)方法之外，其他預(yù)報(bào)方法包括:卡爾曼濾波[7]、自協(xié)方差預(yù)報(bào)[8]、模糊推理系統(tǒng)預(yù)報(bào)[9]、離散小波變換與自協(xié)方差的組合預(yù)報(bào)[10]以及人工神經(jīng)網(wǎng)絡(luò)等[11].極移與地球各個(gè)圈層之間的質(zhì)量重分布和相對(duì)運(yùn)動(dòng)有著緊密的聯(lián)系[12].大氣角動(dòng)量(Atmospheric Angular Momentum，AAM)的變化[13-14]、海洋角動(dòng)量(Oceanic Angular Momentum，OAM)的變化[15-16]、陸地水圈角動(dòng)量(Hydrospheric Angular Momentum，HAM)的變化以及海平面角動(dòng)量(Sea-Level Angular Momentum，SLAM)的變化對(duì)于極移有著不同程度的影響[17-18].由于傳統(tǒng)的統(tǒng)計(jì)預(yù)報(bào)方法沒(méi)有考慮地球流體的高頻質(zhì)量變化導(dǎo)致的隨機(jī)項(xiàng)，Dill等[19]在極移的短期預(yù)報(bào)中融合了有效角動(dòng)量函數(shù)(Effective Angular Momentum Functions，EAM)，從而顯著提高了極移短期預(yù)報(bào)的精度;最近，Dill等[20]在極移預(yù)報(bào)中組合了大氣、海洋和陸地水圈的有效角動(dòng)量函數(shù)，并聯(lián)合最小二乘與自回歸方法，進(jìn)一步提高了極移90 d內(nèi)的預(yù)報(bào)精度.本文以Dill等[20]的預(yù)報(bào)方法為基礎(chǔ)，通過(guò)使用分段最小二乘和改進(jìn)的自回歸模型，實(shí)現(xiàn)了極移更高精度的預(yù)報(bào).

本文通過(guò)聯(lián)合有效角動(dòng)量函數(shù)與LS+AR的數(shù)據(jù)外推方法，分別對(duì)短期(未來(lái)1-6 d)和中長(zhǎng)期(未來(lái)7-90 d)極移進(jìn)行了預(yù)報(bào).其中，短期極移預(yù)報(bào)外推了大地測(cè)量學(xué)的角動(dòng)量(Geodetic Angular Momentum，GAM)和EAM之間的殘差，并利用了EAM數(shù)據(jù)模型1-6 d的預(yù)報(bào)優(yōu)勢(shì);中長(zhǎng)期的極移預(yù)報(bào)直接外推了完整EAM(full EAM，GAM和短期預(yù)報(bào)中得到的有效角動(dòng)量之間的拼接);在兩個(gè)預(yù)報(bào)階段中，使用描述地球旋轉(zhuǎn)變化的劉維爾方程(Liouville Equation)實(shí)現(xiàn)了有效角動(dòng)量和極移之間的轉(zhuǎn)換.在第2節(jié)中，闡述了本文所采用的極移預(yù)報(bào)方法，對(duì)LS+AR方法的應(yīng)用進(jìn)行了適當(dāng)?shù)膬?yōu)化，一方面，采用了分段最小二乘以更好地?cái)M合及外推EAM數(shù)據(jù);另一方面，對(duì)于自回歸參數(shù)進(jìn)行了廣泛的測(cè)試和搜索，選取了一組能夠適應(yīng)不同預(yù)報(bào)階段的優(yōu)化參數(shù)，實(shí)現(xiàn)了極移的多參數(shù)并行預(yù)報(bào).本文第3節(jié)展示了多次預(yù)報(bào)實(shí)驗(yàn)的結(jié)果，在2019年9月至2021年2月的范圍內(nèi)進(jìn)行了441次極移預(yù)報(bào)，并將預(yù)報(bào)結(jié)果分別與EOP 14 C04(IAU2000A)(下文簡(jiǎn)稱(chēng)EOP C04)和IERS Bulletin A(Bulletin A)進(jìn)行了比較，完成了預(yù)報(bào)結(jié)果的評(píng)估.第4節(jié)對(duì)全文進(jìn)行了總結(jié)和展望.

2 方法與數(shù)據(jù)

2.1 有效角動(dòng)量及劉維爾方程

1-6 d的短期極移預(yù)報(bào)利用了EAM的預(yù)報(bào)值，并且對(duì)GAM和EAM之間的殘差(GAMresidual，大地測(cè)量學(xué)的有效角動(dòng)量殘差)進(jìn)行了外推.本文所采用的有效角動(dòng)量函數(shù)EAM包括以下幾個(gè)部分:AAM、OAM、HAM以及SLAM;EAM數(shù)據(jù)來(lái)源于德國(guó)地學(xué)研究中心地球系統(tǒng)模型組(Earth System Modelling group at Deutsches GeoForschungsZentrum，ESMGFZ)維護(hù)的ESMGFZ產(chǎn)品網(wǎng)站http://rz-vm115.gfz-potsdam.de:8080/repositor y.求解GAM residual是短期極移預(yù)報(bào)中的關(guān)鍵步驟，為此，需要首先獲得GAM，劉維爾方程示意極移和角動(dòng)量函數(shù)之間的微分關(guān)系[1]，如下所示:

其中，t為時(shí)間變量，極移的表達(dá)式為p(t)=pX(t)-ipY(t)，極移的X分量為pX(t)，極移X正變化量的方向指向本初子午線(xiàn);極移的Y分量為pY(t)，極移Y正變化量的方向指向西經(jīng)90°W.有效角動(dòng)量函數(shù)的表達(dá)式為χ(t)=χX(t)+iχY(t)，χX(t)和χY(t)分別是有效角動(dòng)量函數(shù)的X分量和Y分量，χX(t)正變化量的方向指向本初子午線(xiàn)，χY(t)正變化量的方向指向東經(jīng)90°E.可見(jiàn)，極移p(t)中Y分量pY(t)前的負(fù)號(hào)表達(dá)了極移與有效角動(dòng)量函數(shù)之間正變化量方向的不一致性.σc表示錢(qián)德勒復(fù)頻率[21]，并且σc=2π[1+i/(2Q)]/Tcw，本文中，錢(qián)德勒擺動(dòng)周期的取值為T(mén)cw=433 d，阻尼因子取值為Q=170[22-24].

Wilson進(jìn)一步推導(dǎo)出了下面的劉維爾方程的離散形式[25]:

以上兩式中，Δt為時(shí)間增量.

由于極移數(shù)據(jù)和角動(dòng)量數(shù)據(jù)的產(chǎn)品是離散的時(shí)間序列，因此，通過(guò)離散化的劉維爾方程即可實(shí)現(xiàn)極移和有效角動(dòng)量之間的轉(zhuǎn)換，于是GAM可以從極移數(shù)據(jù)中轉(zhuǎn)換得到.本文所使用的極移數(shù)據(jù)來(lái)自IERS提供的事后產(chǎn)品EOP C04，快速產(chǎn)品以及預(yù)報(bào)產(chǎn)品Bulletin A，詳見(jiàn)網(wǎng)站:https://www.iers.org/IERS/EN/DataProducts/data.html.

通過(guò)上文中有效角動(dòng)量函數(shù)與極移之間的轉(zhuǎn)換原理和相應(yīng)的數(shù)據(jù)集，可以得到最近4 yr的GAM數(shù)據(jù)和EAM數(shù)據(jù)如圖1所示.GAM由EOP C04的極移數(shù)據(jù)通過(guò)劉維爾方程得出.由圖可知，GAM residual表現(xiàn)較為穩(wěn)定.

圖1 上面兩幅子圖表示GAM數(shù)據(jù)和EAM數(shù)據(jù)，其中，黑色實(shí)線(xiàn)為GAM，黑色點(diǎn)線(xiàn)為EAM;下面兩幅子圖表示GAM和EAM之間的殘差;左圖表示X分量，右圖表示Y分量.此外，EAM=AAM+OAM+HAM+SLAM，GAM和EAM都是無(wú)量綱的.Fig.1 The top two sub-figures represent GAM data and EAM data，in which the black solid line is GAM，and the black dotted line is EAM;the bottom two sub-figures represent the GAM residual;the left figure represents the X component，and the right figure represents the Y component.In addition，EAM=AAM+OAM+HAM+SLAM，GAM and EAM are dimensionless.

此外，EAM的分量AAM、OAM、HAM、SLAM的數(shù)據(jù)特征如表1所示，可見(jiàn)各個(gè)分量數(shù)據(jù)的量級(jí)在10-9-10-7之間.

表1 2016—2021年EAM各分量的數(shù)據(jù)特征Table 1 Data characteristics of each component of the EAM from 2016 to 2021

2.2 最小二乘與自回歸組合方法

在極移的短期預(yù)報(bào)和中長(zhǎng)期預(yù)報(bào)的數(shù)據(jù)外推階段，皆使用了最小二乘和自回歸組合的方法LS+AR.為了更加準(zhǔn)確地?cái)M合極移短期預(yù)報(bào)中的GAM residual數(shù)據(jù)以及極移中長(zhǎng)期預(yù)報(bào)中的GAM數(shù)據(jù)，本文在擬合數(shù)據(jù)時(shí)采用了分段最小二乘的方法.最小二乘擬合公式包括線(xiàn)性項(xiàng)和周期項(xiàng)，對(duì)于不同的數(shù)據(jù)段采用了不同的線(xiàn)性擬合公式，如下所示:

圖2 灰色實(shí)線(xiàn)為GAM residual，黑色實(shí)線(xiàn)為L(zhǎng)S擬合曲線(xiàn)，黑色點(diǎn)線(xiàn)為L(zhǎng)S擬合誤差;左圖是X分量，右圖是Y分量.Fig.2 The gray solid line represents the GAM residual，the black solid line is LS fitting curve，and the black dotted line is LS fitting error;the left picture is the X component，and the right picture is the Y component.

圖3 灰色實(shí)線(xiàn)為full EAM，黑色實(shí)線(xiàn)為L(zhǎng)S擬合曲線(xiàn)，黑色點(diǎn)線(xiàn)為L(zhǎng)S擬合誤差;左圖是X分量，右圖是Y分量.Fig.3 The gray solid line represents the full EAM，the black solid line is LS fitting curve，and the black dotted line is LS fitting error;the left picture is the X component，and the right picture is the Y component.

由圖2-3可見(jiàn)，在使用分段最小二乘的方法對(duì)GAM residual以及full EAM進(jìn)行擬合時(shí)，兩個(gè)階段都產(chǎn)生了較為平穩(wěn)的擬合誤差.實(shí)際上，這些最小二乘擬合誤差接近于隨機(jī)噪聲，但由于其量級(jí)與擬合數(shù)據(jù)集(GAMresidual或者full EAM)相當(dāng)，因此，最小二乘擬合誤差的外推就顯得尤為重要，對(duì)于GAM residual和full EAM的外推有著顯著的影響.為了得到更好的最小二乘擬合誤差的外推值，本文在使用自回歸處理擬合誤差的基礎(chǔ)上，對(duì)于擬合誤差數(shù)據(jù)進(jìn)行了選擇性過(guò)濾，也就是從擬合誤差數(shù)據(jù)序列中間隔選取數(shù)據(jù)，之后再進(jìn)行自回歸操作，此方法實(shí)現(xiàn)了在數(shù)據(jù)選擇的同時(shí)保證了數(shù)據(jù)信息的完整度，具體方法如下所述:

在外推最小二乘擬合誤差時(shí)，本文使用了自回歸方法，并且對(duì)最小二乘擬合誤差時(shí)間序列進(jìn)行了濾波，使得在不同的極移預(yù)報(bào)階段以及不同分量的極移預(yù)報(bào)中，自回歸模型能夠間隔選取時(shí)間序列，以便于在更大范圍的參數(shù)空間中為不同的極移預(yù)報(bào)階段和不同分量的極移預(yù)報(bào)選擇更優(yōu)的參數(shù)，從而獲得更優(yōu)的預(yù)報(bào)結(jié)果.將實(shí)現(xiàn)時(shí)間序列濾波的參數(shù)稱(chēng)為間隔數(shù).因此，在如下式所示的自回歸模型中，可調(diào)節(jié)參數(shù)有階數(shù)和間隔數(shù)兩個(gè)，模型如下所示:

其中，sign(i)是符號(hào)函數(shù):

y(t)表示單變量時(shí)序數(shù)據(jù)，c是自回歸模型的常數(shù)項(xiàng)，p是自回歸階數(shù)，i表示數(shù)據(jù)的序號(hào)，φi是自回歸系數(shù)，?t為白噪聲;lag是自回歸間隔數(shù)，取值范圍是正整數(shù).sign(i)是與間隔數(shù)lag有關(guān)的符號(hào)函數(shù)，實(shí)現(xiàn)了在求解自回歸模型時(shí)，可以在階數(shù)p的范圍內(nèi)以lag為間隔選取時(shí)間序列，相當(dāng)于對(duì)時(shí)間序列進(jìn)行了濾波.自回歸模型的階數(shù)p和間隔數(shù)lag是本文極移預(yù)報(bào)中的重要參數(shù)，下文中將詳細(xì)論述.

2.3 極移預(yù)報(bào)流程

本文中，極移的預(yù)報(bào)總體上可以分為兩個(gè)階段，一是1-6 d的短期預(yù)報(bào);二是7-90 d的中長(zhǎng)期預(yù)報(bào).短期預(yù)報(bào)利用了EAM數(shù)據(jù)1-6 d的預(yù)報(bào)值，首先，劉維爾方程將極移數(shù)據(jù)轉(zhuǎn)換成GAM數(shù)據(jù)，與EAM數(shù)據(jù)做差，得到GAM residual.然后，通過(guò)LS+AR的方法將GAM residual外推至未來(lái)第6 d，EAM數(shù)據(jù)6 d的預(yù)報(bào)值與GAM residual數(shù)據(jù)6 d的外推值相加，獲得了未來(lái)6 d的有效角動(dòng)量預(yù)報(bào)值，再次利用劉維爾方程，即可將有效角動(dòng)量預(yù)報(bào)值轉(zhuǎn)換為極移的預(yù)報(bào).中長(zhǎng)期極移預(yù)報(bào)的關(guān)鍵在于外推GAM數(shù)據(jù)，將短期預(yù)報(bào)中得到的未來(lái)6 d的有效角動(dòng)量預(yù)報(bào)值與GAM數(shù)據(jù)的歷史值拼接起來(lái)，形成full EAM，利用LS+AR的方法對(duì)full EAM外推至未來(lái)第90 d，并通過(guò)劉維爾方程的轉(zhuǎn)換得到7-90 d的極移預(yù)報(bào).在極移的整個(gè)預(yù)報(bào)過(guò)程中，極移p(t)和有效角動(dòng)量χ(t)以復(fù)數(shù)的形式作為整體參與計(jì)算.在Dill等人研究的基礎(chǔ)上[20]，經(jīng)過(guò)多次實(shí)驗(yàn)測(cè)試，在短期預(yù)報(bào)中，最小二乘線(xiàn)性項(xiàng)的擬合數(shù)據(jù)段長(zhǎng)度選取為1 yr，周期項(xiàng)的擬合數(shù)據(jù)段長(zhǎng)度選取為4 yr，自回歸模型的擬合數(shù)據(jù)段長(zhǎng)度選取為4 yr;在中長(zhǎng)期預(yù)報(bào)中，最小二乘線(xiàn)性項(xiàng)的擬合數(shù)據(jù)段長(zhǎng)度選取為1 yr，周期項(xiàng)的擬合數(shù)據(jù)段長(zhǎng)度選取為4 yr，自回歸模型的擬合數(shù)據(jù)段長(zhǎng)度選取為1 yr.

3 結(jié)果與討論

3.1 自回歸模型參數(shù)

由上文可知，在整個(gè)極移預(yù)報(bào)過(guò)程中，自回歸模型有兩個(gè)可調(diào)節(jié)參數(shù)，分別是階數(shù)p和間隔數(shù)lag.在1-6 d的極移預(yù)報(bào)階段，將X分量和Y分量的自回歸階數(shù)固定為60階，使X分量和Y分量的間隔數(shù)在區(qū)間[1，20]的整數(shù)范圍內(nèi)變化，選擇符合條件的階數(shù)和間隔數(shù)，可以得到354組自回歸模型可調(diào)節(jié)參數(shù)，構(gòu)成了短期極移預(yù)報(bào)的自回歸參數(shù)空間.在7-90 d的極移預(yù)報(bào)階段，使X分量和Y分量的自回歸階數(shù)在區(qū)間[2，20]的整數(shù)范圍內(nèi)變化，同時(shí)，讓X分量和Y分量的自回歸間隔數(shù)在區(qū)間[1，20]的整數(shù)范圍內(nèi)變化，選擇符合條件的階數(shù)和間隔數(shù)，同樣可以得到354組自回歸模型可調(diào)節(jié)參數(shù)，構(gòu)成了中長(zhǎng)期極移預(yù)報(bào)的自回歸參數(shù)空間.兩個(gè)預(yù)報(bào)階段參數(shù)空間中的不同參數(shù)對(duì)于極移預(yù)報(bào)的影響如圖4所示.

圖4上方的兩個(gè)子圖表明，在1-6 d的短期極移預(yù)報(bào)中，自回歸模型的可調(diào)整參數(shù)p和lag的選擇對(duì)于短期極移預(yù)報(bào)的MAE(Mean Absolute Error)有較為明顯的影響，因此，短期極移預(yù)報(bào)的自回歸參數(shù)空間中存在更優(yōu)的參數(shù).同樣的，圖4下方的兩個(gè)子圖表明，在7-90 d的中長(zhǎng)期極移預(yù)報(bào)中，自回歸可調(diào)整參數(shù)p和lag的選擇對(duì)于中長(zhǎng)期極移預(yù)報(bào)的MAE有明顯的影響，所以，中長(zhǎng)期極移預(yù)報(bào)的自回歸參數(shù)空間中也存在更優(yōu)的參數(shù).于是，可以從參數(shù)空間中選擇出更優(yōu)的極移預(yù)報(bào)參數(shù).

圖4 不同的自回歸模型參數(shù)所對(duì)應(yīng)的極移預(yù)報(bào)值的MAEFig.4 The MAE of the polar motion prediction values corresponding to different AR model parameters

通過(guò)對(duì)上述的短期和中長(zhǎng)期極移預(yù)報(bào)自回歸參數(shù)空間中的參數(shù)進(jìn)行搜索，并且利用多次預(yù)報(bào)實(shí)驗(yàn)比較各組不同參數(shù)對(duì)于極移預(yù)報(bào)的影響，就可以選出兩組更優(yōu)的自回歸參數(shù)，如表2和表3所示.表2-3給出了兩組從參數(shù)空間中得到的優(yōu)化參數(shù)，其中，表2的參數(shù)用于極移X的預(yù)報(bào)，表3的參數(shù)用于極移Y的預(yù)報(bào)，使得對(duì)于不同分量的極移在不同階段的預(yù)報(bào)能夠采用更優(yōu)的參數(shù)，下文中的預(yù)報(bào)實(shí)驗(yàn)則使用了這兩組優(yōu)化參數(shù).

表2 對(duì)于極移X預(yù)報(bào)更優(yōu)的自回歸參數(shù)Table 2 Better autoregressive parameters for PMX prediction

表3 對(duì)于極移Y預(yù)報(bào)更優(yōu)的自回歸參數(shù)Table 3 Better autoregressive parameters for PMY prediction

由此，對(duì)于未來(lái)1-90 d中任意一天的極移預(yù)報(bào)，能夠?qū)?yīng)一組更優(yōu)的自回歸模型參數(shù)p和lag;例如，在預(yù)報(bào)未來(lái)第6 d的極移X時(shí)，自回歸參數(shù)選取為p=(60，60)、lag=(1，15);在預(yù)報(bào)未來(lái)第6 d的極移Y時(shí)，自回歸參數(shù)選取為p=(19，19)、lag=(1，1).可見(jiàn)，為了獲得更優(yōu)的結(jié)果，1-90 d的極移預(yù)報(bào)分別對(duì)應(yīng)多組不同的自回歸參數(shù)，這同時(shí)增加了預(yù)報(bào)過(guò)程中的計(jì)算次數(shù)，降低了極移預(yù)報(bào)的計(jì)算速度，為了解決此問(wèn)題，本文采用了不同參數(shù)分組并行的計(jì)算方法.

3.2 與IERS EOP C04比較

在時(shí)間跨度為2019年9月至2021年2月的范圍內(nèi)，采用上文中所述的相關(guān)參數(shù)，進(jìn)行了441次極移預(yù)報(bào)實(shí)驗(yàn)，將計(jì)算結(jié)果與EOP C04提供的極移數(shù)據(jù)相比較，如圖5所示.由圖可見(jiàn)，極移X在1-90 d的預(yù)報(bào)誤差分布在0-36 mas的范圍內(nèi)，并且10 d內(nèi)的預(yù)報(bào)誤差大約在4.5 mas以下;極移Y在1-90 d的預(yù)報(bào)誤差分布在0-16 mas的范圍內(nèi)，并且10 d內(nèi)的預(yù)報(bào)誤差大約在2 mas以下.此外，極移的預(yù)報(bào)效果與預(yù)報(bào)時(shí)間段有關(guān)，在有些月份預(yù)報(bào)得好(例如2019年9月)，在另一些月份預(yù)報(bào)得不好(例如2020年8月).

圖5 2019年9月至2021年2月之間的441次極移預(yù)報(bào)實(shí)驗(yàn)中1-90 d的預(yù)報(bào)誤差Fig.5 The prediction error from 1 to 90 days in 441 PM prediction experiments between September 2019 and February 2021

3.3 與IERS Bulletin A的比較

在時(shí)間跨度為2019年9月至2021年2月的范圍內(nèi)，進(jìn)行了441次極移預(yù)報(bào)實(shí)驗(yàn)，在實(shí)驗(yàn)過(guò)程中，采用了上文中所述的相關(guān)參數(shù)，例如，對(duì)于自回歸參數(shù)的選擇，極移X分量的預(yù)報(bào)中使用了表2所示的參數(shù)，極移Y分量的預(yù)報(bào)中則使用了表3所示的參數(shù).將多次預(yù)報(bào)結(jié)果與Bulletin A提供的極移預(yù)報(bào)數(shù)據(jù)相比較(以EOP C04為參考)，如圖6所示.

由圖6可知，整體上，相比于極移X，極移Y的預(yù)報(bào)優(yōu)于Bulletin A所占的比例更大.在極移X的預(yù)報(bào)中，短期有50%以上的預(yù)報(bào)結(jié)果優(yōu)于Bulletin A，這一比例在第6 d左右達(dá)到最大，隨著預(yù)報(bào)日期的增加而下降，30 d之后Bulletin A的預(yù)報(bào)占據(jù)了更多的優(yōu)勢(shì).在極移Y的預(yù)報(bào)中，短期有60%以上的預(yù)報(bào)結(jié)果優(yōu)于Bulletin A，同樣，這一比例在第6 d左右達(dá)到最大，雖然隨著預(yù)報(bào)日期的增加這一比例有所下降，但總體保持在50%以上.由此可見(jiàn)，與Bulletin A相比，極移X的預(yù)報(bào)在1-30 d更占優(yōu)勢(shì)，極移Y的預(yù)報(bào)則在1-90 d的所有預(yù)報(bào)階段均占優(yōu)勢(shì).更進(jìn)一步，將極移預(yù)報(bào)劃分為3個(gè)不同的階段:1-6 d、7-30 d、31-90 d，以EOP C04為參考，與Bulletin A的比較如圖7所示.

圖7 極移X預(yù)報(bào)與Bulletin A的分階段對(duì)比圖，白色紋理表示1-6 d的預(yù)報(bào)階段，灰色紋理表示7-30 d的預(yù)報(bào)階段，點(diǎn)狀紋理表示31-90 d的預(yù)報(bào)階段.Fig.7 The phase-by-phase comparison of PMX prediction and Bulletin A，the white texture represents the prediction phase from 1 to 6 days，the gray texture represents the forecast phase from 7 to 30 days，and the dotted texture represents the forecast phase from 31 to 90 days.

圖7中展示了不同階段的極移預(yù)報(bào)與Bulletin A之間的比較，左圖表示極移X，右圖表示極移Y.橫坐標(biāo)是與Bulletin A之間的差值，即本文相較于Bulletin A的預(yù)報(bào)精度提升量，其中，負(fù)值表示預(yù)報(bào)誤差小于Bulletin A，0值表示預(yù)報(bào)誤差與Bulletin A相差無(wú)幾，正值表示預(yù)報(bào)誤差大于Bulletin A.以左圖為例，在橫坐標(biāo)0值處的白色紋理柱子高度約為45%，那么可以認(rèn)為有45%的極移X預(yù)報(bào)與Bulletin A是一樣好的;同樣的，在橫坐標(biāo)-1處的白色紋理柱子高度約為20%，表明在1-6 d的極移X預(yù)報(bào)中，有20%的極移X預(yù)報(bào)誤差比Bulletin A小1 mas左右.仍然以白色紋理柱子為例，圖中左上角圖例的百分比表示在所有的1-6 d的極移X預(yù)報(bào)中，有56.9%的極移X預(yù)報(bào)誤差要小于Bulletin A.此外，灰色紋理、點(diǎn)狀紋理以及極移Y的含義同上.由此可見(jiàn)，在441次預(yù)報(bào)實(shí)驗(yàn)中，極移X的1-6 d預(yù)報(bào)、7-30 d預(yù)報(bào)和31-90 d預(yù)報(bào)比Bulletin A更優(yōu)的占比分別為56.9%、53.5%和40.5%，并且，大部分預(yù)報(bào)精度的提升分布在[-3，0]mas的范圍內(nèi);極移Y的1-6 d預(yù)報(bào)、7-30 d預(yù)報(bào)和31-90 d預(yù)報(bào)比Bulletin A更優(yōu)的占比分別為66.5%、59.7%和59.2%，并且，大部分預(yù)報(bào)精度的提升量分布在[-2，0]mas的范圍內(nèi).

4 總結(jié)與展望

本文通過(guò)分段最小二乘擬合及改進(jìn)的自回歸模型，使用ESMGFZ提供的數(shù)據(jù)，提升了1-90 d的極移預(yù)報(bào).將AAM、OAM、HAM和SLAM組成的EAM融合進(jìn)極移的預(yù)報(bào)，同時(shí)，選擇更優(yōu)的自回歸模型可調(diào)節(jié)參數(shù)，在多組更優(yōu)參數(shù)的條件下并行預(yù)報(bào)，從441次極移預(yù)報(bào)實(shí)驗(yàn)來(lái)看，極移的預(yù)報(bào)精度提升是比較明顯的.表4為本文中441次預(yù)報(bào)的MAE(以EOP C04為參考)和Bulletin A相應(yīng)預(yù)報(bào)的MAE(以EOP C04為參考)之間的比較.可見(jiàn)，整體上極移Y的預(yù)報(bào)效果更好;相比于IERS預(yù)報(bào)，極移X的預(yù)報(bào)精度在30 d內(nèi)均有提升，其中，在第1 d和第5 d的提升分別為2.62%和32.98%;極移Y的預(yù)報(bào)精度在90 d內(nèi)有明顯的提升，其中，在第1 d和第5 d的提升分別為20.77%和48.98%.綜上，采用本文中所述的方法，從極移X分量和極移Y分量的角度看，極移Y的預(yù)報(bào)更有優(yōu)勢(shì);從短期和長(zhǎng)期的角度看，短期預(yù)報(bào)更有優(yōu)勢(shì).

表4 不同階段的極移預(yù)報(bào)精度(以EOP C04為參考)Table 4 PM prediction accuracy at different stages(taking EOP C04 as reference)

本文中，極移Y的預(yù)報(bào)優(yōu)于極移X的預(yù)報(bào)，相比于極移Y，極移X的預(yù)報(bào)誤差更大.從有效角動(dòng)量數(shù)據(jù)的角度來(lái)看，X分量的數(shù)據(jù)規(guī)律性更弱，噪聲更大，造成這一現(xiàn)象的主要原因是地球表面的陸地相對(duì)于Y軸的分布比X軸更對(duì)稱(chēng).極移的預(yù)報(bào)結(jié)果與時(shí)間段也存在著一定的關(guān)系，在有些時(shí)間段預(yù)報(bào)得更好，在另外一些時(shí)間段預(yù)報(bào)得不好.在后續(xù)的研究中，或許可以通過(guò)進(jìn)一步研究極移的數(shù)據(jù)特征，尋找更好的最小二乘擬合模型，結(jié)合卡爾曼濾波或者人工智能等方法進(jìn)一步改善極移的預(yù)報(bào)精度.

致謝感謝Dill博士在本文的研究過(guò)程中對(duì)于若干問(wèn)題提供的悉心解答和建議.感謝ESMGFZ和IERS在數(shù)據(jù)上的支持.