雷 雨 趙丹寧 喬?；?徐勁松 蔡宏兵

(1 西安郵電大學(xué)計(jì)算機(jī)學(xué)院 西安 710121)(2 寶雞文理學(xué)院電子電氣工程學(xué)院 寶雞 721016)(3 中國(guó)科學(xué)院國(guó)家授時(shí)中心 西安 710600)(4 江蘇師范大學(xué)江蘇圣理工學(xué)院 徐州 221116)

1 引言

地球自轉(zhuǎn)變化可以用一組地球定向參數(shù)(Earth Orientation Parameters,EOP)來(lái)表征,包括極移(polar motion,PM)、日長(zhǎng)變化、歲差與章動(dòng).其中極移是地球瞬時(shí)自轉(zhuǎn)軸相對(duì)于地球本體內(nèi)的運(yùn)動(dòng),是實(shí)現(xiàn)地球參考坐標(biāo)系和天球參考坐標(biāo)系相互轉(zhuǎn)換的必要參數(shù)之一,在衛(wèi)星導(dǎo)航、深空探測(cè)及跟蹤等空間科學(xué)工程應(yīng)用領(lǐng)域具有不可或缺的作用.此外,極移與地球上物質(zhì)分布等地球物理現(xiàn)象存在密切關(guān)系,因此,極移可以為全球性的地球物理現(xiàn)象提供重要信息.雖然甚長(zhǎng)基線干涉測(cè)量(Very Long Baseline Interferometry,VLBI)、全球衛(wèi)星導(dǎo)航系統(tǒng)(Global Navigation Satellite System,GNSS)以及衛(wèi)星激光測(cè)距(Satellite Laser Ranging,SLR)等空間測(cè)地技術(shù)能獲取高精度的極移數(shù)據(jù),但因復(fù)雜的資料處理過(guò)程,在極移觀測(cè)時(shí)效性上存在一定的滯后性,難以獲得實(shí)時(shí)的極移觀測(cè)量,所以只能通過(guò)預(yù)報(bào)手段實(shí)現(xiàn)極移參數(shù)的實(shí)時(shí)獲取.衛(wèi)星導(dǎo)航和深空探測(cè)器跟蹤、導(dǎo)航等空間科學(xué)工程應(yīng)用領(lǐng)域?qū)Ω呔葮O移預(yù)報(bào)數(shù)據(jù)有重大需求.其中深空探測(cè)及跟蹤通常需要1–7 d的極移預(yù)報(bào)數(shù)據(jù)以支持深空探測(cè)器跟蹤和導(dǎo)航,而北斗衛(wèi)星導(dǎo)航一般需要1–60 d的極移預(yù)報(bào)數(shù)據(jù)作為衛(wèi)星自主定軌的先驗(yàn)信息.因此,如何提高極移短期和中期預(yù)報(bào)精度成為天文地球動(dòng)力學(xué)領(lǐng)域的一個(gè)研究熱點(diǎn)[1–2].

國(guó)內(nèi)外學(xué)者提出了多種極移預(yù)報(bào)方法,例如最小二乘(Least Squares,LS)外推[3–5]、自回歸(Autoregressive,AR)[6–8]、神 經(jīng) 網(wǎng) 絡(luò)[9–11]、小 波 分析[12]和奇異譜分析模型[13]等.國(guó)際地球定向參數(shù)預(yù)報(bào)比較競(jìng)賽(Earth Orientation Parameters Prediction Comparison Campaign,EOP PCC)結(jié)果表明,最小二乘和自回歸模型的組合是預(yù)報(bào)極移的較優(yōu)模型,但該模型無(wú)法對(duì)任意跨度的極移預(yù)報(bào)均取得最優(yōu)[14].當(dāng)前,國(guó)際地球自轉(zhuǎn)服務(wù)局(International Earth Rotation and Reference Systems Service,IERS)A公報(bào)(Bulletin A)中發(fā)布的極移預(yù)報(bào)產(chǎn)品就是利用最小二乘和自回歸模型組合獲取的[15].

受大氣和海洋等多種激發(fā)因素的影響,極移序列不僅包含周年極移與錢德勒極移等確定周期性成分,還包含準(zhǔn)周期性甚至非線性不規(guī)則成分,其中,周期性成分預(yù)報(bào)精度較高,而準(zhǔn)周期性和非線性成分具有非線性、非平穩(wěn)特征,是影響極移預(yù)報(bào)精度的主要因素[16].混沌現(xiàn)象是確定性系統(tǒng)中因系統(tǒng)隨機(jī)性而產(chǎn)生的外在的、復(fù)雜及貌似無(wú)規(guī)律的行為,混沌系統(tǒng)對(duì)初始值敏感的特性使混沌系統(tǒng)輸入的變化可以快速地反映在系統(tǒng)輸出中,因此混沌理論提供一種更符合現(xiàn)實(shí)世界的非線性建模方法.混沌預(yù)報(bào)法作為一種新穎的非線性預(yù)報(bào)理論,通過(guò)對(duì)時(shí)間序列進(jìn)行相空間重構(gòu),能夠挖掘時(shí)間序列中隱含的規(guī)律,非常適合于總體呈現(xiàn)確定性但又具有隨機(jī)性的復(fù)雜系統(tǒng),所以混沌時(shí)間序列預(yù)報(bào)方法用于極移預(yù)報(bào)是合理可行的.混沌時(shí)間序列預(yù)報(bào)方法可以分為局部預(yù)報(bào)法[17]、全局預(yù)報(bào)法[18]和自適應(yīng)預(yù)報(bào)法[19].其中,局部預(yù)報(bào)法只利用部分樣本數(shù)據(jù)擬合非線性函數(shù),計(jì)算速度快但預(yù)報(bào)精度低;全局預(yù)報(bào)法利用全部樣本數(shù)據(jù)擬合函數(shù),當(dāng)樣本數(shù)據(jù)量大時(shí)計(jì)算速度慢且易過(guò)擬合;自適應(yīng)預(yù)報(bào)法能夠自適應(yīng)地跟蹤混沌運(yùn)動(dòng)軌跡,可以根據(jù)混沌時(shí)間序列的變化自動(dòng)調(diào)整模型參數(shù),具有較高的預(yù)報(bào)精度.而針對(duì)混沌時(shí)間序列的Volterra自適應(yīng)預(yù)報(bào)法能夠有效地預(yù)測(cè)低維混沌時(shí)間序列,已在電力、交通和天氣預(yù)報(bào)等領(lǐng)域取得成功應(yīng)用[20–22].

本文探索用新型的混沌動(dòng)力學(xué)理論來(lái)研究極移變化.首先,把只有一維的極移時(shí)間序列通過(guò)相空間重構(gòu)方法重構(gòu)出高維空間,恢復(fù)極移內(nèi)在的非線性特性;其次,對(duì)呈現(xiàn)非線性特性變化的極移序列應(yīng)用混沌理論分析與預(yù)報(bào),通過(guò)對(duì)極移序列的相空間重構(gòu)計(jì)算Lyapunov指數(shù),確定極移為混沌時(shí)間序列的前提下,結(jié)合Taken嵌入定理,確定Volterra自適應(yīng)預(yù)測(cè)濾波器的長(zhǎng)度,利用Volterra自適應(yīng)濾波對(duì)極移時(shí)序進(jìn)行預(yù)測(cè).文中對(duì)極移進(jìn)行了未來(lái)60 d的預(yù)報(bào),并將預(yù)報(bào)結(jié)果與EOP PCC和IERS A公報(bào)的預(yù)報(bào)結(jié)果進(jìn)行了對(duì)比,對(duì)比結(jié)果驗(yàn)證了本文方法的預(yù)報(bào)效果.

2 時(shí)間序列相空間重構(gòu)及混沌識(shí)別

常規(guī)的時(shí)間序列通常是在時(shí)間域中進(jìn)行時(shí)間序列的分析與處理,對(duì)于混沌時(shí)間序列的研究在相空間中進(jìn)行.所謂混沌時(shí)間序列,可以視作是考察混沌系統(tǒng)所得到的一組隨著時(shí)間而變化的觀察值.根據(jù)Taken嵌入定理,可以從一維混沌時(shí)間序列中重構(gòu)一個(gè)和原動(dòng)力系統(tǒng)在拓?fù)湟饬x下相同的相空間,混沌時(shí)間序列的識(shí)別、分析以及預(yù)測(cè)均在重構(gòu)的相空間中進(jìn)行,因此,混沌時(shí)間序列的相空間重構(gòu)是混沌時(shí)間序列分析與處理的關(guān)鍵.

2.1 相空間重構(gòu)

Packard等人于1980年提出了時(shí)間序列相空間重構(gòu)的兩種方法[23]:導(dǎo)數(shù)重構(gòu)法與坐標(biāo)延遲重構(gòu)法,其中,坐標(biāo)延遲重構(gòu)是通過(guò)一維時(shí)間序列的不同延遲時(shí)間τ來(lái)構(gòu)建D維的相空間矢量.設(shè)有長(zhǎng)度為N的一維時(shí)間序列x(t),t=1,2,···,N,嵌入維數(shù)為m,則時(shí)間序列x(t)重構(gòu)的相空間為

式中,M=N-(m-1)τ,X為M×m維矩陣,任一相點(diǎn)X(t)均包含m個(gè)元素,相空間的演化軌跡就是M個(gè)相點(diǎn)在m維空間的連線.相空間重構(gòu)的關(guān)鍵是選取適當(dāng)?shù)摩优cm.τ的選取方法有互信息法和自相關(guān)法,m的選取方法有最近鄰點(diǎn)法與飽和關(guān)聯(lián)維數(shù)法.近年研究表明,τ與m是相關(guān)的,聯(lián)合選取可以提高相空間重構(gòu)的質(zhì)量[24].本文利用C-C方法聯(lián)合求取τ與m,該方法使用關(guān)聯(lián)積分同時(shí)估計(jì)出時(shí)延與嵌入窗,具體過(guò)程見(jiàn)文獻(xiàn)[24].

2.2 混沌識(shí)別

時(shí)間序列的混沌識(shí)別主要是通過(guò)計(jì)算混沌特性參數(shù)來(lái)識(shí)別序列的混沌特性,混沌參數(shù)有K熵、關(guān)聯(lián)維數(shù)與Lyapunov指數(shù)等.本文選擇Lyapunov指數(shù)作為混沌參數(shù)來(lái)識(shí)別極移時(shí)間序列的混沌特性,如果最大Lyapunov指數(shù)(Largest Lyapunov Exponent,LLE)為正數(shù)表明時(shí)間序列具有混沌特性,且最大Lyapunov指數(shù)越大,混沌特性越強(qiáng).

基于小數(shù)量法計(jì)算最大Lyapunov指數(shù)的步驟如下[25]:

(1)對(duì)時(shí)間序列進(jìn)行快速傅里葉變換估計(jì)平均周期T;

(2)基于坐標(biāo)延遲重構(gòu)法重構(gòu)時(shí)間序列的相空間,對(duì)相空間中的每個(gè)點(diǎn)X(t),計(jì)算該相點(diǎn)和最鄰近點(diǎn)X(?t)的l個(gè)離散時(shí)間步長(zhǎng)后的距離dt(l),t為相點(diǎn)X(t)的歷元,?t為相點(diǎn)X(t)的最鄰近點(diǎn)的歷元,dt(l)可以表示為

式中,l=1,2,···,min(M-t,M-?t),?t=1,2,···,M,且t/=?t;

(3)取dt(l)的對(duì)數(shù)lndt(l),對(duì)每個(gè)l求所有非零lndt(l)的均值再除以時(shí)間序列的采樣周期Δt,即

式中,q為非零dt(l)的個(gè)數(shù),i為非零dt(l)的序號(hào);

3 Volterra自適應(yīng)預(yù)測(cè)濾波器

對(duì)于極移預(yù)報(bào),其本質(zhì)上是建立樣本“輸入–輸出”之間的映射關(guān)系,極移預(yù)報(bào)的準(zhǔn)確與否在一定程度上也是由樣本“輸入–輸出”之間的擬合程度決定的.在極移混沌時(shí)間序列預(yù)報(bào)中,設(shè)樣本輸入為重構(gòu)相空間中的相點(diǎn)X(n)=[x(n),x(n-τ),···,x(n-(m-1)τ)],輸出為x(n+1),n為歷元,則理論上存在一個(gè)映射模型F(·)使得:

由于Volterra級(jí)數(shù)能夠同時(shí)對(duì)時(shí)間序列中的線性行為和非線性行為進(jìn)行擬合,本文利用Volterra級(jí)數(shù)構(gòu)建樣本“輸入–輸出”映射模型來(lái)擬合F(·),Volterra級(jí)數(shù)的展開(kāi)表達(dá)式為[26]

式中,i1,i2,···,ip為滯后歷元,hp(i1,i2,···,ip)為p階Volterra核函數(shù).這種無(wú)窮級(jí)數(shù)展開(kāi)表達(dá)式在實(shí)際應(yīng)用中較難實(shí)現(xiàn),需要運(yùn)用有限次截?cái)嗯c有限次求和的形式對(duì)Volterra級(jí)數(shù)展開(kāi)式進(jìn)行簡(jiǎn)化,較為常用的是二階截?cái)嗪颓蠛偷男问?即

式中,N1和N2為濾波器長(zhǎng)度.對(duì)于混沌時(shí)間序列,根據(jù)Taken嵌入定理:一個(gè)混沌時(shí)間序列至少需要滿足m≥2D+1才能完全描述原系統(tǒng)的動(dòng)力學(xué)特征.因此,應(yīng)將N1和N2取為N1=N2=m≥2D+1,則用于極移預(yù)報(bào)的二階Volterra濾波預(yù)測(cè)器為

定義FIR(Finite Impulse Response)濾波器的輸入矢量U(n)=[1,x(n),x(n-τ),···,x(n-(m-1)τ),x2(n),x(n)x(n-τ),···,x2(n-(m-1)τ)]T,核函數(shù)向量H(n)=[h0,h1(0),h1(1),h1(m-1),h2(0,0),h2(0,1),···,h2(m-1,m-1)T],二 階Volterra濾波預(yù)測(cè)器可表示為

上式中的核函數(shù)向量可直接采用FIR濾波器的自適應(yīng)算法來(lái)確定,本文選用最小均方自適應(yīng)算法來(lái)確定核函數(shù)向量:

式中,e(n)為濾波器輸出值和期望值之差,核函數(shù)向量根據(jù)e(n)自適應(yīng)調(diào)整以獲得期望輸出;c為收斂控制參數(shù),可控制濾波器收斂于某個(gè)期望值.

4 預(yù)報(bào)方法及精度評(píng)定

4.1 預(yù)報(bào)方法

基于混沌時(shí)間序列的極移預(yù)報(bào)方法的基本流程如下:

(1)極移數(shù)據(jù)預(yù)處理:利用最小二乘擬合扣除極移時(shí)間序列中的線性趨勢(shì)項(xiàng)、周年項(xiàng)以及錢德勒項(xiàng),獲得趨勢(shì)項(xiàng)、周年項(xiàng)以及錢德勒項(xiàng)的外推值,最小二乘外推模型的數(shù)學(xué)表達(dá)式為

式中,ax和bx表示極移Xp分量的線性趨勢(shì)項(xiàng)參數(shù),和表示Xp的錢德勒項(xiàng)參數(shù),和表示Xp的周年項(xiàng)參數(shù),p1和p2分別表示錢德勒擺動(dòng)和周年擺動(dòng)周期,在擬合中取p1=1.183 yr、p2=1 yr.對(duì)于極移Yp分量,各對(duì)應(yīng)參數(shù)表示含義和Xp分量的模型相同.極移序列最小二乘擬合后的剩余殘差序列包含亞季節(jié)性準(zhǔn)周期性成分和非周期性成分,這些準(zhǔn)周期變化和非線性變化主要由大氣和海洋對(duì)極移的激發(fā)引起;

(2)極移殘差序列相空間重構(gòu):基于坐標(biāo)延遲重構(gòu)法重構(gòu)最小二乘擬合殘差時(shí)間序列的相空間,其中,相空間的τ與m通過(guò)C-C法確定;

(3)極移殘差序列混沌識(shí)別:采用小數(shù)據(jù)量方法計(jì)算最小二乘殘差時(shí)間序列的最大Lyapunov指數(shù),識(shí)別極移殘差的混沌特征;

(4)極移殘差序列預(yù)報(bào):利用二階Volterra自適應(yīng)濾波器對(duì)重構(gòu)相空間的極移殘差時(shí)間序列進(jìn)行預(yù)測(cè).在建模階段,濾波器輸入為殘差序列重構(gòu)相空間中的相點(diǎn)X(n)=[x(n),x(n-τ),···,x(n-(m-1)τ)],輸出為對(duì)應(yīng)輸入向量一步后的下一狀態(tài)x(n+1),利用最小均方自適應(yīng)算法建立輸入X(n)和輸出x(n+1)之間的映射關(guān)系;在預(yù)報(bào)階段,將濾波器預(yù)測(cè)值?x(n+1)加入到輸入向量X(n)中,并刪除向量中第1個(gè)元素,然后將其更新為X(n+1),以此逐步遞推,即可實(shí)現(xiàn)極移的多步預(yù)測(cè).圖1為基于混沌時(shí)間序列的極移預(yù)報(bào)流程.

圖1 基于混沌時(shí)間序列的極移預(yù)報(bào)流程Fig.1 Flowchart of PM prediction based on chaotic time series

4.2 精度評(píng)定

選取平均絕對(duì)誤差(mean absolute error,MAE)作為評(píng)價(jià)本文預(yù)報(bào)方法精度的指標(biāo)[14]:

式中,k表示預(yù)報(bào)跨度,L表示預(yù)報(bào)期數(shù),和分別表示第j期的第k天的極移預(yù)報(bào)值與觀測(cè)值.

5 極移預(yù)報(bào)分析

5.1 極移數(shù)據(jù)來(lái)源

采用國(guó)際地球自轉(zhuǎn)服務(wù)局發(fā)布的EOP C04序列和A公報(bào)資料進(jìn)行極移分析與預(yù)報(bào)研究,EOP C04序列包含1962年至今的極移觀測(cè)數(shù)據(jù),該觀測(cè)數(shù)據(jù)是通過(guò)VLBI、GNSS和SLR等多種空間測(cè)地觀測(cè)資料聯(lián)合解算得到的,極移觀測(cè)精度約為50 μas,每周更新兩次;A公報(bào)包含未來(lái)1 yr的極移預(yù)報(bào)數(shù)據(jù),該預(yù)報(bào)數(shù)據(jù)是通過(guò)最小二乘外推和自回歸模型組合外推獲得的,由IERS快速服務(wù)與產(chǎn)品中心每周更新一次.IERS EOP C04序列和A公報(bào)公布的極移數(shù)據(jù)的采樣間隔均為1 d,具有較高的精度和可靠性,是開(kāi)展極移分析與預(yù)報(bào)研究的常用資料來(lái)源.

5.2 極移混沌特性分析

為驗(yàn)證極移的混沌特性,選取EOP 14C04序列中1997年1月1日至2020年12月31日的極移數(shù)據(jù)進(jìn)行分析,將該時(shí)間段內(nèi)的極移數(shù)據(jù)每隔6 yr劃分為一個(gè)觀測(cè)時(shí)段,共分為4個(gè)觀測(cè)時(shí)段,以極移Xp分量為例.圖2給出了每個(gè)觀測(cè)時(shí)段內(nèi)的極移觀測(cè)序列、最小二乘擬合序列及其最小二乘擬合殘差序列,其中,實(shí)線表示極移觀測(cè)序列,點(diǎn)線表示極移最小二乘擬合序列,虛線表示極移最小二乘擬合殘差序列.從圖中可以看到,扣除線性趨勢(shì)項(xiàng)、周年項(xiàng)以及錢德勒項(xiàng)后的極移最小二乘擬合殘差呈現(xiàn)為明顯的非線性變化,主要包含亞季節(jié)性準(zhǔn)周期成分和剩余的不規(guī)則成分.

圖2 極移Xp分量觀測(cè)序列、最小二乘擬合序列及其殘差序列Fig.2 Observed,fitted and residual time series of the PM Xp component

根據(jù)小數(shù)據(jù)量方法分別計(jì)算上述4個(gè)時(shí)段的極移最小二乘擬合殘差序列的最大Lyapunov指數(shù),計(jì)算方法為對(duì)(3)式的計(jì)算結(jié)果ˉd(t)進(jìn)行最小二乘擬合,擬合直線的斜率即為最大Lyapunov指數(shù).限于篇幅,本文僅給出1997年1月1日至2002年12月31日的極移Xp和Yp序列的ˉd(t)及其擬合結(jié)果,如圖3所示,圖中的陰影部分為所選的線性擬合區(qū)間,擬合斜率分別為0.0275和0.0159,其余時(shí)段的ˉd(t)擬合斜率(最大Lyapunov指數(shù))見(jiàn)表1,表中LLE表示最大Lyapunov指數(shù).由表1可知,扣除線性趨勢(shì)項(xiàng)、周年項(xiàng)以及錢德勒項(xiàng)后的極移殘差序列具有低強(qiáng)度的混沌特性,說(shuō)明極移Xp、Yp分量時(shí)間序列可視為混沌序列.本文將扣除線性趨勢(shì)項(xiàng)、周年項(xiàng)以及錢德勒項(xiàng)后的極移最小二乘擬合殘差作為一組混沌序列,利用二階Volterra自適應(yīng)濾波器對(duì)混沌序列中的有序行為進(jìn)行外推預(yù)報(bào),然后將最小二乘外推結(jié)果和Volterra預(yù)報(bào)結(jié)果合并即為最終極移預(yù)測(cè)序列.為驗(yàn)證本文方法的極移預(yù)報(bào)效果,將本文 預(yù)報(bào)結(jié)果與EOP PCC和IERS A公報(bào)的預(yù)報(bào)結(jié)果進(jìn)行比較分析.

圖3 ˉd隨時(shí)間步長(zhǎng)的變化及其擬合區(qū)間Fig.3variations with the time step and the fitted interval

表1 不同觀測(cè)時(shí)段的極移序列的最大Lyapunov指數(shù)Table 1 Largest Lyapunov exponent of the PM in different observed periods

5.3 與EOP PCC比較

為評(píng)估各種方法的EOP預(yù)報(bào)效果,維也納理工大學(xué)于2005年10月1日至2008年2月28日組織了國(guó)際性的EOP預(yù)報(bào)比較競(jìng)賽,國(guó)際上共有12個(gè)小組提交EOP預(yù)報(bào)數(shù)據(jù),涉及20種預(yù)報(bào)方法.維也納理工大學(xué)對(duì)20種EOP預(yù)報(bào)方法進(jìn)行分析評(píng)估,通過(guò)分析評(píng)估發(fā)現(xiàn)極移預(yù)報(bào)精度較高的方法是最小二乘外推和自回歸組合模型,其中,超短期(1–10 d)預(yù)報(bào)精度最高的小組是Kosek 051小組和Kalarus 061小組,使用的是最小二乘外推與自回歸模型組合方法;短期(1–30 d)預(yù)報(bào)精度最高的小組是Kosek 051小組、Kalarus 061小組和Zotov 091小組,Zotov 091小組使用的是自回歸模型;中長(zhǎng)期(30 d以上)預(yù)報(bào)效果最好的小組是Kosek 051小組和Kalarus 061小組,其他小組編號(hào)及其使用的預(yù)報(bào)方法詳見(jiàn)文獻(xiàn)[10],此處不再贅述.

首先隨機(jī)選擇兩個(gè)時(shí)段進(jìn)行極移預(yù)報(bào),選擇的兩個(gè)時(shí)段分別是2006年11月25日至2007年1月23日和2007年5月19日至2007年7月17日,極移兩個(gè)分量的預(yù)報(bào)結(jié)果和預(yù)報(bào)誤差如圖4–5所示,圖中虛線表示本文極移預(yù)報(bào)結(jié)果,實(shí)線表示EOP 05C04序列極移觀測(cè)值,點(diǎn)線表示預(yù)報(bào)值和觀測(cè)值之差.從圖4–5中可以看出,本文方法1–30 d的極移預(yù)報(bào)值和觀測(cè)值符合得很好,絕對(duì)誤差小于5 mas,當(dāng)預(yù)報(bào)跨度超過(guò)30 d時(shí),預(yù)報(bào)值逐漸偏離觀測(cè)值,預(yù)報(bào)誤差表現(xiàn)出發(fā)散趨勢(shì),絕對(duì)誤差可以達(dá)到12 mas.

圖4 本文方法Xp預(yù)報(bào)值和EOP 05C04序列極移觀測(cè)值的比較Fig.4 Comparison between the Xp predictions by the proposed method and the EOP 05C04 observations

為了與EOP PCC預(yù)報(bào)結(jié)果進(jìn)行客觀的比較,本文預(yù)報(bào)了與EOP預(yù)報(bào)比較競(jìng)賽同時(shí)段的極移數(shù)據(jù),在進(jìn)行預(yù)報(bào)結(jié)果比較時(shí),采用EOP 05C04序列中的極移數(shù)據(jù)作為基礎(chǔ)數(shù)據(jù),基礎(chǔ)數(shù)據(jù)長(zhǎng)度取為周年極移和錢德勒極移的最小公倍數(shù)6 yr,選取平均絕對(duì)誤差作為精度評(píng)估指標(biāo).混沌時(shí)間序列預(yù)測(cè)方法的1–10 d、1–30 d和1–60 d預(yù)報(bào)結(jié)果和EOP PCC預(yù)報(bào)結(jié)果的平均絕對(duì)誤差如圖6所示,EOP PCC預(yù)報(bào)結(jié)果由維也納理工大學(xué)提供.由于各小組提交的EOP預(yù)報(bào)數(shù)據(jù)長(zhǎng)度不一致,如Zotov 091小組僅提交超短期和短期預(yù)報(bào)結(jié)果,而未提交中長(zhǎng)期預(yù)報(bào)結(jié)果,所以圖6中1–10 d、1–30 d和1–60 d預(yù)報(bào)結(jié)果比較中所涉及到平均絕對(duì)誤差曲線數(shù)量并不相同.圖6中,Kosek 051小組預(yù)報(bào)結(jié)果的平均絕對(duì)誤差用紅色實(shí)線表示,Kalarus 061小組的平均絕對(duì)誤差用黑色實(shí)線表示,Zotov 091小組用藍(lán)色實(shí)線表示,本文方法的平均絕對(duì)誤差在圖6中用紅色虛線表示,圖中所涉及的1–10 d平均絕對(duì)誤差根據(jù)約100期預(yù)報(bào)結(jié)果獲得,而1–30 d和1–60 d平均絕對(duì)誤差根據(jù)約50期預(yù)報(bào)結(jié)果獲得.表2給出了Kosek 051小組、Kalarus 061小組、Zotov 091小組和本文預(yù)報(bào)結(jié)果的平均絕對(duì)誤差數(shù)值.

從圖6中和表2可以發(fā)現(xiàn),本文方法的極移1–10 d超短期預(yù)報(bào)精度和Kosek 051小組、Kalarus 061小組相當(dāng),顯著優(yōu)于其他小組;對(duì)于10–30 d短期預(yù)報(bào),本文方法的平均絕對(duì)預(yù)報(bào)誤差與預(yù)報(bào)精度最高的Kosek 051小組、Kalarus 061小組和Zotov 091小組基本相當(dāng);在30–60 d中期預(yù)報(bào)中,本文方法的預(yù)報(bào)效果不如排在第1位Kalarus 061小組,極移Xp分量的平均絕對(duì)預(yù)報(bào)誤差與排在第2位Kosek 051小組的預(yù)報(bào)誤差處于同一水平,低于其他小組的平均預(yù)報(bào)誤差;對(duì)于極移Yp分量,30–50 d的平均絕對(duì)預(yù)報(bào)誤差與Kosek 051小組大體相當(dāng),當(dāng)預(yù)報(bào)跨度超過(guò)50 d時(shí),極移Yp分量的預(yù)報(bào)誤差發(fā)散趨勢(shì)明顯,預(yù)報(bào)精度不及Kosek小組.從表2還可以看出,沒(méi)有任何一種預(yù)報(bào)方法對(duì)于任意跨度和任何分量的預(yù)報(bào)精度都能夠達(dá)到最優(yōu),說(shuō)明各種方法均存在各自的局限性.

表2 本文預(yù)報(bào)結(jié)果和EOP PPC極移預(yù)報(bào)結(jié)果的平均絕對(duì)誤差Table 2 Mean absolute error of the polar motion predictions obtained by the proposed method and EOP PPC

圖5 本文方法Yp預(yù)報(bào)值和EOP 05C04序列極移觀測(cè)值的比較Fig.5 Comparison between the Yp predictions by the proposed method and the EOP 05C04 observations

圖6 本文方法和EOP PCC極移預(yù)報(bào)誤差的比較Fig.6 Comparison of the PM prediction errors between the proposed method and EOP PCC

5.4 與IERS A公報(bào)比較

自2017年1月起IERS利用EOP 14C05序列作為極移基礎(chǔ)數(shù)據(jù)發(fā)布常規(guī)極移預(yù)報(bào)值,本文選取EOP 14C05序列中的極移觀測(cè)資料作為建模數(shù)據(jù)進(jìn)行極移未來(lái)60 d預(yù)報(bào),建模數(shù)據(jù)長(zhǎng)度為6 yr,將預(yù)報(bào)結(jié)果與IERS A公報(bào)2021年第20周和25周發(fā)布的極移預(yù)報(bào)數(shù)據(jù)進(jìn)行對(duì)比,并將EOP 14C04序列極移觀測(cè)數(shù)據(jù)作為參考值.圖7給出了本文方法與IERS A公報(bào)極移預(yù)報(bào)值以及EOP 14C04序列極移觀測(cè)值的對(duì)比結(jié)果,圖中點(diǎn)線表示IERS A公報(bào)極移預(yù)報(bào)值,虛線表示本文方法的極移預(yù)報(bào)值,實(shí)線表示EOP 14C04序列極移觀測(cè)值.從圖7可以看到,本文方法在預(yù)報(bào)前期與EOP 14C04極移序列符合得很好,和IERS A公報(bào)預(yù)報(bào)效果相當(dāng)甚至更好,而在預(yù)報(bào)后期本文預(yù)測(cè)值逐漸偏離實(shí)際值,其中,本文方法對(duì)于極移Xp分量1–40 d預(yù)報(bào)精度較高,預(yù)報(bào)跨度超過(guò)40 d時(shí)誤差較大;本文方法對(duì)于極移Yp分量預(yù)報(bào)的改進(jìn)沒(méi)有Xp分量那樣明顯,1–20 d預(yù)報(bào)結(jié)果與IERS A公報(bào)比較接近,預(yù)報(bào)跨度超過(guò)20 d時(shí)本文方法的預(yù)報(bào)誤差開(kāi)始增大.整體而言,本文方法對(duì)于極移短期預(yù)報(bào)效果較好,當(dāng)預(yù)報(bào)跨度增加時(shí)預(yù)報(bào)結(jié)果不如IERS A公報(bào).

圖7 本文方法和IERS A公報(bào)極移預(yù)報(bào)值的比較Fig.7 Comparison of the predicted PM values by the proposed method and IERS Bulletin A

為進(jìn)一步對(duì)比分析本文方法與IERS A公報(bào)的極移預(yù)測(cè)結(jié)果,將本文預(yù)報(bào)結(jié)果與IERS A公報(bào)2021年第1周至第35周發(fā)布的極移預(yù)報(bào)結(jié)果進(jìn)行比較,將IERS EOP 14C04序列作為參考值,分別統(tǒng)計(jì)極移Xp分量和Yp分量的平均絕對(duì)預(yù)報(bào)誤差,統(tǒng)計(jì)結(jié)果見(jiàn)表3.從表3可以看到,本文方法和IERS A公報(bào)對(duì)于不同跨度的極移預(yù)報(bào)結(jié)果存在一定程度的差異,相對(duì)于IERS A公報(bào)發(fā)布的預(yù)報(bào)值,本文方法的預(yù)報(bào)結(jié)果在預(yù)報(bào)前期效果有所提高,其中,對(duì)于極移Xp分量,未來(lái)1 d預(yù)報(bào)結(jié)果提高最為明顯,預(yù)報(bào)精度比A公報(bào)提高15%.除第2 d、第3 d和第4 d外,本文方法對(duì)于未來(lái)30 d以內(nèi)其他跨度的預(yù)報(bào)有1%–12%程度的提高.對(duì)于極移Yp分量也是如此,只是本文方法對(duì)于Yp分量的預(yù)報(bào)改善程度沒(méi)有Xp分量明顯,改善程度在6%以內(nèi).無(wú)論是對(duì)于Xp分量還是Yp分量,在預(yù)報(bào)后期,即當(dāng)預(yù)報(bào)跨度超過(guò)一定的天數(shù)后,本文方法的極移預(yù)報(bào)結(jié)果相對(duì)于IERS A公報(bào)非但沒(méi)有提高,反而有所下降,其中,當(dāng)預(yù)報(bào)跨度超過(guò)30 d時(shí),本文方法對(duì)極移Xp分量的預(yù)報(bào)結(jié)果開(kāi)始惡化;當(dāng)預(yù)報(bào)跨度超過(guò)20 d時(shí),本文方法對(duì)極移Yp分量的預(yù)報(bào)效果逐漸降低,且Yp分量的預(yù)報(bào)效果降低得更為明顯,這說(shuō)明隨著預(yù)報(bào)跨度的增加,本文方法的極移預(yù)報(bào)誤差比IERS A公報(bào)增大得更快.根據(jù)第4.1節(jié)所述的極移預(yù)報(bào)過(guò)程,結(jié)合理論分析發(fā)現(xiàn),出現(xiàn)這種現(xiàn)象的可能原因在于極移預(yù)報(bào)值是由Volterra自適應(yīng)濾波器逐步遞推得到的,這說(shuō)明上一步濾波器輸出預(yù)報(bào)值的預(yù)報(bào)誤差會(huì)累積到下一步的預(yù)報(bào)值中,預(yù)報(bào)誤差會(huì)出現(xiàn)逐步傳遞的情況,換言之,預(yù)報(bào)結(jié)果存在誤差累積效應(yīng),且隨著預(yù)報(bào)跨度的增加,這種誤差累積效應(yīng)會(huì)愈發(fā)明顯.從表3還可以發(fā)現(xiàn),本文方法對(duì)于極移Xp分量預(yù)報(bào)的改善程度比Yp分量更顯著,這是由于Xp分量的混沌強(qiáng)度比Yp分量要強(qiáng)(見(jiàn)表1),因此基于混沌時(shí)間序列的預(yù)測(cè)方法對(duì)Xp分量具有更好的適用性.

表3 本文預(yù)報(bào)結(jié)果和IERS A公報(bào)極移預(yù)報(bào)結(jié)果的平均絕對(duì)誤差Table 3 Mean absolute error of the polar motion predictions obtained by the proposed method and IERS Bulletin A

6 總結(jié)與展望

極移由趨勢(shì)成分、周期成分與隨機(jī)成分組成,是多種激發(fā)源的物理機(jī)制共同作用的結(jié)果,考慮到極移這種復(fù)雜的非線性時(shí)變特征,提出基于混沌時(shí)間序列的極移預(yù)報(bào)方法.首先利用最小二乘擬合算法分離極移的確定項(xiàng)和殘差項(xiàng),隨后根據(jù)最大Lyapunov指數(shù)對(duì)極移殘差項(xiàng)的混沌特性進(jìn)行識(shí)別,計(jì)算結(jié)果表明:極移殘差項(xiàng)的最大Lyapunov指數(shù)為正數(shù),且數(shù)值較小,表明極移殘差項(xiàng)具有低強(qiáng)度的混沌特性.在此基礎(chǔ)上,構(gòu)建二階Volterra自適應(yīng)濾波器對(duì)極移殘差項(xiàng)進(jìn)行預(yù)測(cè),并疊加極移確定項(xiàng)的外推值得到極移預(yù)報(bào)值.將本文方法的極移預(yù)報(bào)結(jié)果和EOP PCC預(yù)報(bào)結(jié)果以及IERS A公報(bào)發(fā)布的極移預(yù)報(bào)產(chǎn)品分別進(jìn)行了對(duì)比,分析發(fā)現(xiàn):

(1)和EOP PPC極移預(yù)報(bào)結(jié)果相比,本文方法對(duì)于1–30 d的極移預(yù)報(bào)結(jié)果與EOP PPC最優(yōu)預(yù)報(bào)方法相當(dāng),當(dāng)預(yù)報(bào)跨度超過(guò)30 d時(shí),本文方法的預(yù)報(bào)精度低于EOP PCC最優(yōu)預(yù)報(bào)方法,但優(yōu)于參與EOP PCC的其他方法;

(2)和IERS A公報(bào)發(fā)布的極移預(yù)報(bào)產(chǎn)品相比,本文方法的前期預(yù)報(bào)結(jié)果相對(duì)于A公報(bào)有不同程度提高,當(dāng)達(dá)到一定的預(yù)報(bào)跨度以后,本文方法的預(yù)報(bào)精度出現(xiàn)明顯的降低,總體而言,本文方法更適合于極移的短期預(yù)報(bào).

通過(guò)分析發(fā)現(xiàn),本文方法在后期預(yù)報(bào)中誤差增大的主要原因可能是由Volterra自適應(yīng)濾波器的誤差累積效應(yīng)導(dǎo)致的,因此,如何設(shè)計(jì)更為優(yōu)良的Volterra自適應(yīng)濾波器,從而提高極移的預(yù)報(bào)精度是下一步研究的重點(diǎn).