張 志，廖 瑛，余 越

(1．國防科技大學(xué)航天科學(xué)與工程學(xué)院，湖南長沙 410073;2．海軍蚌埠士官學(xué)校4系，安徽蚌埠 233012)

極移是表征地球自轉(zhuǎn)的一個(gè)重要物理量，它與歲差、章動(dòng)和日長變化一起構(gòu)成地球定向參數(shù)(Earth Orientation Parameters，EOP)，它在航天器軌道確定、衛(wèi)星導(dǎo)航定位等領(lǐng)域有重要的應(yīng)用［1］?，F(xiàn)代空間大地測量技術(shù)(衛(wèi)星激光測距、甚長基線干涉測量、全球定位系統(tǒng)等)是目前極移參數(shù)獲取的主要手段，但是由于復(fù)雜的分析處理過程，不能實(shí)時(shí)獲取測量數(shù)據(jù)，因此，有必要對(duì)以極移為代表的EOP參數(shù)進(jìn)行預(yù)報(bào)。

在極移預(yù)報(bào)方面，學(xué)者們提出了很多有意義的方法。Zhu［2］和 Chao［3］采用確定和不確定周期項(xiàng)構(gòu)建最小二乘諧波模型進(jìn)行極移預(yù)報(bào);Kosek［4-5］將自回歸模型應(yīng)用于最小二乘殘差序列中，采用LS模型與自回歸模型組合方法進(jìn)行預(yù)報(bào)。前述方法均為線性模型預(yù)報(bào)，隨著非線性科學(xué)的發(fā)展，以神經(jīng)網(wǎng)絡(luò)等為代表的非線性預(yù)報(bào)方法和理論研究取得了一些突破性的進(jìn)展。1992年，Egger［6］首先將神經(jīng)網(wǎng)絡(luò)應(yīng)用于地球自轉(zhuǎn)參數(shù)(Earth Rotation Parameters，ERP)的預(yù)報(bào)中;之后Egger和 Frehlich［7］又將該方法與 Fr?ehlich 的完全解析方法相比較，結(jié)果驗(yàn)證了神經(jīng)網(wǎng)絡(luò)在預(yù)報(bào)ERP上與完全解析方法的一致性，顯示了神經(jīng)網(wǎng)絡(luò)在預(yù)報(bào)不規(guī)則和準(zhǔn)周期過程中的巨大潛力;近年來，Schuh［8］、Liao［9］和王琪潔［10］均采用神經(jīng)網(wǎng)絡(luò)技術(shù)進(jìn)行EOP的預(yù)報(bào):Schuh和Liao分別采用斯圖加特神經(jīng)網(wǎng)絡(luò)模擬軟件和簡化的三層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)對(duì)EOP進(jìn)行預(yù)報(bào)，王琪潔在極移預(yù)報(bào)中采用“先—后—最佳拓?fù)洹钡姆绞竭M(jìn)行預(yù)報(bào)，三位學(xué)者均取得了比較好的預(yù)報(bào)效果。

文獻(xiàn)［9］采用簡化的三層網(wǎng)絡(luò)結(jié)構(gòu)對(duì)EOP作了短期與中長期預(yù)測，其中，輸入層選擇了6個(gè)神經(jīng)元，輸出層選擇了12個(gè)神經(jīng)元，即該神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)為多輸入-多輸出結(jié)構(gòu)。文獻(xiàn)［9］中每個(gè)輸入樣本的節(jié)點(diǎn)數(shù)為6個(gè)，而網(wǎng)絡(luò)的輸出預(yù)報(bào)點(diǎn)數(shù)比輸入節(jié)點(diǎn)還多，在輸出預(yù)報(bào)精度上不能很好保證。由此，考慮增加每個(gè)輸入樣本的節(jié)點(diǎn)數(shù)并減少神經(jīng)網(wǎng)絡(luò)的輸出預(yù)報(bào)點(diǎn)數(shù)是提高預(yù)報(bào)精度的有效方法之一。文中對(duì)極移數(shù)據(jù)的特點(diǎn)進(jìn)行了研究，根據(jù)極移數(shù)據(jù)的特點(diǎn)首先對(duì)極移序列重采樣處理，將得到的新序列進(jìn)行頻譜分析驗(yàn)證了重采樣間隔的合理性，提取新序列的趨勢項(xiàng)，然后采用多輸入-單輸出BP神經(jīng)網(wǎng)絡(luò)對(duì)新序列殘差進(jìn)行不同跨度預(yù)報(bào)，組合趨勢項(xiàng)和殘差預(yù)報(bào)得到最終的極移預(yù)報(bào)值。

1 頻譜分析

進(jìn)行頻譜分析用到的數(shù)據(jù)來源于國際地球自轉(zhuǎn)和參考系服務(wù)(International Earth Rotation and Reference Systems Service，IERS)發(fā)布的EOPC04序列①。該序列包含 1962年1月 1日(MJD 37665)至今的極移(包括x和y分量，分別簡寫為PMX和PMY)、世界時(shí)UT1-UTC和日長LOD數(shù)據(jù)等，采樣間隔為1天。

本文研究的是極移預(yù)報(bào)，所以要研究極移數(shù)據(jù)的周期性變化特性，也就是頻率域上的特征，選用傅里葉分析方法，離散傅里葉變換公式為:

其中，k=1，2，…，N;i是虛數(shù)單位;x 是觀測時(shí)間序列，此處為極移分量序列;N是時(shí)間序列x的長度。

首先對(duì)EOP C04序列數(shù)據(jù)(一般地，對(duì)采樣間隔為1d的極移數(shù)據(jù)，稱為基礎(chǔ)序列)進(jìn)行傅里葉變化，得到傅里葉變換幅值譜如圖1、圖2所示。然后對(duì)EOPC04序列進(jìn)行重采樣，選取采樣間隔為10d的極移序列(稱為插值基礎(chǔ)序列)進(jìn)行傅里葉分析，得到的傅里葉變化幅值譜如圖3，4所示。

由圖1～圖4可知，重采樣得到的10d間隔極移數(shù)據(jù)的幅值譜圖與基礎(chǔ)序列為1d間隔的幅值譜圖基本一致，均出現(xiàn)兩個(gè)明顯的周期-振幅峰值，說明極移是具有周期特性的，相應(yīng)的峰值橫坐標(biāo)對(duì)應(yīng)于極移的周期值，考慮到周期的大小(錢德勒周期大于周年周期)，兩個(gè)周期分別對(duì)應(yīng)周年周期和錢德勒周期。其中as表示角秒，重采樣得到極移的幅值譜分析很好地符合基礎(chǔ)序列特性，表現(xiàn)了極移的周期特性，因而，在對(duì)極移處理過程中可以采用采樣間隔為10d的極移數(shù)據(jù)。在極移預(yù)報(bào)中，這大大減少了極移數(shù)據(jù)點(diǎn)數(shù)目以及隨后的計(jì)算處理時(shí)間。

2 極移模型

根據(jù)前述頻譜分析結(jié)果，極移的三角函數(shù)模型形式表示如式(2)所示。

式中，x(t)是給定歷元t時(shí)的極移趨勢項(xiàng);A，B，C1，C2，D1，D2是待求參數(shù);pc和 pa分別為錢德勒擺動(dòng)和周年擺動(dòng)的周期。

采用插值基礎(chǔ)序列觀測值對(duì)式(2)進(jìn)行最小二乘擬合［11］，可知極移趨勢項(xiàng)x(t)包括線性趨勢項(xiàng)A+Bt，錢德勒項(xiàng)(含 pc項(xiàng))和周年項(xiàng)(含 pa項(xiàng))三個(gè)部分。扣除趨勢項(xiàng)的插值基礎(chǔ)序列觀測值稱為殘差序列數(shù)據(jù)，接下來將采用BP神經(jīng)網(wǎng)絡(luò)模型處理得到的殘差序列。

3 BP神經(jīng)網(wǎng)絡(luò)

由于提取趨勢項(xiàng)后的殘差序列數(shù)據(jù)具有非線性特性，采用BP神經(jīng)網(wǎng)絡(luò)建模預(yù)報(bào)是因其具有逼近非線性函數(shù)能力強(qiáng)和網(wǎng)絡(luò)收斂速度快的優(yōu)點(diǎn)。

BP神經(jīng)網(wǎng)絡(luò)的本質(zhì)是找出輸入和輸出之間未知但又存在的函數(shù)關(guān)系。從輸出層節(jié)點(diǎn)個(gè)數(shù)分類可以將BP網(wǎng)絡(luò)劃分為兩種結(jié)構(gòu):有多個(gè)輸出節(jié)點(diǎn)的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)［8-9］，這是當(dāng)前地球定向參數(shù)預(yù)報(bào)中比較常見的一種結(jié)構(gòu)，即多輸入-多輸出結(jié)構(gòu);另外一種結(jié)構(gòu)是多輸入-單輸出結(jié)構(gòu)，即只有一個(gè)輸出節(jié)點(diǎn)的網(wǎng)絡(luò)結(jié)構(gòu)，這是本文研究的內(nèi)容。

多輸入-單輸出BP神經(jīng)網(wǎng)絡(luò)的基本結(jié)構(gòu)如圖5所示，這是一種前向網(wǎng)絡(luò)，輸入層由信號(hào)源節(jié)點(diǎn)組成，第二層為隱含層，第三層為輸出層，它對(duì)輸入模式的作用做出響應(yīng)。

由圖5所示，要精確預(yù)報(bào)數(shù)據(jù)需要采用與輸出相關(guān)性好的參數(shù)作為輸入，分析極移數(shù)據(jù)發(fā)現(xiàn)，極移殘差序列的變化是一個(gè)漸變的過程，因此，根據(jù)歷史變化趨勢來預(yù)測后面的輸出可以得到較好的效果。文中構(gòu)建的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)是以若干連續(xù)實(shí)際值作為輸入值，預(yù)報(bào)下一時(shí)刻的輸出值。預(yù)報(bào)模型表示為:

圖1 基礎(chǔ)序列分量x的幅值譜Fig．1 Amplitude spectrum of PMX for basic series

圖2 基礎(chǔ)序列分量y的幅值譜Fig．2 Amplitude spectrum of PMY for basic series

圖3 插值基礎(chǔ)序列分量x的幅值譜Fig．3 Amplitude spectrum of PMX for interpolated basic series

式中，x(t)為給定歷元t時(shí)的殘差數(shù)據(jù)，x(t-1)為最近的過去時(shí)刻的數(shù)據(jù)，兩者時(shí)間間隔為10d;r為輸入節(jié)點(diǎn)數(shù);F為由神經(jīng)網(wǎng)絡(luò)確定的輸入-輸出映射關(guān)系。

圖4 插值基礎(chǔ)序列分量y的幅值譜Fig．4 Amplitude spectrum of PMY for interpolated basic series

圖5 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig．5 Structure of BP neural network

輸入節(jié)點(diǎn)數(shù)r，即歷史數(shù)據(jù)長度的選擇對(duì)輸出有很大的影響，r的取值范圍為5≤r≤15。隱含層節(jié)點(diǎn)數(shù)s的取值范圍為1≤s≤round(s max)，round是四舍五入取整符號(hào)，s max的取值［12］為:

式中，n為輸出層節(jié)點(diǎn)數(shù)，n的取值為1。輸入節(jié)點(diǎn)數(shù)r和隱含層節(jié)點(diǎn)數(shù)s根據(jù)神經(jīng)網(wǎng)絡(luò)的性能在各自取值范圍內(nèi)進(jìn)行調(diào)整，本文根據(jù)全局誤差準(zhǔn)則來訓(xùn)練［13］，以使BP網(wǎng)絡(luò)擬合訓(xùn)練數(shù)據(jù)最佳。

4 誤差精度評(píng)定

極移預(yù)報(bào)結(jié)果常用的誤差精度評(píng)定標(biāo)準(zhǔn)有平均絕對(duì)誤差(Mean Absolute Error，MAE)和均方根誤差(Root Mean Squared Error，RMSE)［14］，則預(yù)報(bào)跨度為i的MAE和RMSE的計(jì)算公式為:

其中，N為預(yù)報(bào)次數(shù)，o為極移觀測值，p為預(yù)報(bào)值，在本文研究中i=10，20，…，360，即預(yù)報(bào)跨度i的取值為10～360d。

5 預(yù)報(bào)結(jié)果及討論

運(yùn)用訓(xùn)練好的BP神經(jīng)網(wǎng)絡(luò)模型對(duì)殘差序列作不同跨度預(yù)報(bào)，預(yù)報(bào)后的數(shù)據(jù)加上極移趨勢項(xiàng)，得到最終的極移預(yù)報(bào)值。為更好地進(jìn)行試驗(yàn)精度比對(duì)，選用文獻(xiàn)［9］中的預(yù)報(bào)時(shí)間區(qū)間、預(yù)報(bào)結(jié)果用于比較分析。第一次預(yù)報(bào)時(shí)間為2001年4月6日(MJD 52005)到2002年3月31日(MJD 52364)，預(yù)報(bào)跨度為360d的極移;接著數(shù)據(jù)序列向前推移91d，繼續(xù)重復(fù)以上形式到下一個(gè)預(yù)報(bào)期間，最終共計(jì)預(yù)報(bào)37次，最后一次預(yù)報(bào)區(qū)間為2010年3月26日(MJD 55281)到2011年3月20日(MJD 55640)。限于篇幅，表1只列出了跨度為10～120d預(yù)報(bào)結(jié)果的MAE和RMSE。圖6～圖9繪出了本文方法與文獻(xiàn)［9］中方法的極移預(yù)報(bào)結(jié)果的MAE和RMSE值，其中mas表示毫角秒。

圖6 PMX預(yù)報(bào)的MAEFig．6 MAE of the predictions of PMX

圖7 PMX預(yù)報(bào)的RMSEFig．7 RMSE of the predictions of PMX

圖8 PMY預(yù)報(bào)的MAEFig．8 MAE of the predictions of PMY

圖9 PMY預(yù)報(bào)的RMSEFig．9 RMSE of the predictions of PMY

表1 極移預(yù)報(bào)的MAE和RMSETab．1 MAE and RMSE of the prediction for x and y components of polarmotion(PMX，PMY)

由表1可以看出，對(duì)于10～120d的預(yù)報(bào)跨度，本文方法在極移分量x，y的預(yù)報(bào)中MAE小于22mas，RMSE小于27mas。綜合圖6～圖9可知，10～120d的預(yù)報(bào)跨度，較文獻(xiàn)［9］方法，本文方法具有明顯優(yōu)勢。對(duì)于10～360d的跨度，采用本文方法的極移預(yù)報(bào)RMSE值優(yōu)于文獻(xiàn)［9］中的結(jié)果，這說明本文極移預(yù)報(bào)方法穩(wěn)定性更好。極移y分量預(yù)報(bào)的MAE值優(yōu)于文獻(xiàn)［9］中的結(jié)果，極移x分量的MAE在跨度150d之前和300d之后優(yōu)于文獻(xiàn)［9］，而在150～300d之間精度相當(dāng)。

結(jié)合圖表可以看出，運(yùn)用BP神經(jīng)網(wǎng)絡(luò)方法預(yù)報(bào)極移得到了很好的精度，這主要得益于插值基礎(chǔ)序列采用了10d間隔的極移數(shù)據(jù)，在不損失數(shù)據(jù)信息的情況下可以選取更多的插值基礎(chǔ)序列數(shù)據(jù)預(yù)報(bào)下一時(shí)刻數(shù)據(jù)。BP神經(jīng)網(wǎng)絡(luò)模型在中短期(跨度在120d以內(nèi))預(yù)報(bào)精度較高，隨著預(yù)報(bào)跨度增加，其預(yù)報(bào)精度略有下降，但總體效果較好。

6 結(jié)論

本文利用EOP C04的極移序列建模進(jìn)行中短期預(yù)報(bào)，其目的是能夠得到未來一段時(shí)間內(nèi)高精度的極移值。對(duì)EOP C04序列重采樣減少了數(shù)據(jù)點(diǎn)數(shù)目，而頻譜分析結(jié)果表明數(shù)據(jù)信息量并未減少，采樣間隔為10d的插值基礎(chǔ)序列可以作為極移預(yù)報(bào)的數(shù)據(jù)來源。對(duì)插值基礎(chǔ)序列提取趨勢項(xiàng)有利于建模預(yù)報(bào)，利用多輸入-單輸出的BP神經(jīng)網(wǎng)絡(luò)建模對(duì)不同跨度預(yù)報(bào)效果良好。

本文預(yù)報(bào)方法能夠得到較好的中短期極移預(yù)報(bào)值。然而，本文僅使用極移觀測序列的信息來建立模型實(shí)施預(yù)報(bào)，如果能結(jié)合地球自轉(zhuǎn)的物理特性，可使預(yù)報(bào)更加完善。

References)

［1］ Chin T M，Gross R S，Dickey JO．Modeling and forecast of the polarmotion excitation functions for short-term polarmotion prediction［J］．Journal of Geodesy，2004，78(6):343-353．

［2］ Zhu SY．Prediction of polarmotion［J］．Bulletin Géodésique，1982，56(3):258-273．

［3］ Chao B F．Predictability of the Earth's polar motion［J］．Bulletin Géodésique，1985，59(1):81-93．

［4］ Kosek W，Mccarthy D D，Johnson T J，et al．Comparison of polar motion prediction results supplied by the IERS subbureau for rapid service and predictions and results of other prediction methods［C］//Proceedings of Journees St．Petersburg:2003．

［5］ Kosek W，Kalarus M，Niedzielski T，et al．Forecasting of the Earth orientation parameters-comparison of different algorithms［C］．Paris，F(xiàn)rance，2008．

［6］ Egger D．Neuronales netz pr?dizierterdrotations parameter［J］．Allgemeine Vermessungsnachrichten(AVN)，1992:517-524．

［7］ Egger D，F(xiàn)r?hlich H．Pr?diktion von erdrotationsdaten-klassisch und neuronal［J］．Allgemeine Vermessungsnachrichten(AVN)，1993，10:366-375．

［8］ Schuh H，Ulrich M，Egger D，et al．Prediction of Earth orientation parametersby artificial neural networks［J］．Journal of Geodesy，2002，76(5):247-258．

［9］ Liao D C，Wang Q J，Zhou Y H，et al．Long-term prediction of the Earth orientation parameters by the artificial neural network technique［J］．Journal of Geodynamics，2012，62:87-92．

［10］ 王琪潔．基于神經(jīng)網(wǎng)絡(luò)技術(shù)的地球自轉(zhuǎn)變化預(yù)報(bào)［D］．上海:中國科學(xué)院研究生院上海天文臺(tái)，2007．WANG Qijie．Studies on the prediction of Earth's variable rotation by artificialneuralnetworks［D］．Shanghai:Shanghai Astronomical Observatory， Chinese Academy of Sciences，2007．(in Chinese)

［11］ 張志，廖瑛，文援蘭，等．基于基本多項(xiàng)式的GPS精密星歷插值方法研究［J］．測繪通報(bào)，2014(1):12-15．ZHANG Zhi，LIAO Ying，WEN Yuanlan，etal．Research on method of GPS precise ephemeris interpolation based on standard polynomial［J］．Bulletin of Surveying and Mapping，2014(1):12-15．(in Chinese)

［12］ 高大啟．有教師的線性基本函數(shù)前向三層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)研究［J］．計(jì)算機(jī)學(xué)報(bào)，1998，21(1):80-86．GAO Daqi．On structures of supervised linear basis function feed forward three-layered neural networks［J］．Chinese Journal of Computers，1998，21(1):80-86．(in Chinese)

［13］ Simon H．神經(jīng)網(wǎng)絡(luò)原理［M］．2版．葉世偉，史忠植，譯．北京:機(jī)械工業(yè)出版社，2004．Simon H．Neural networks［M］．2nd ed．Translated by YE Shiwei，SHIZhongzhi．Beijing:China Machine Press，2004．(in Chinese)

［14］ Niedzielski T，Kosek W．Prediction of UT1-UTC，LOD and AAMχ3 by combination of least-squares and multivariate stochastic methods［J］．Journal of Geodesy，2008，82(2):83-92．