張 志，廖 瑛，余 越

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

極移是表征地球自轉(zhuǎn)的一個重要物理量，它與歲差、章動和日長變化一起構(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ù)據(jù)，因此，有必要對以極移為代表的EOP參數(shù)進(jìn)行預(yù)報。

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

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

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)、世界時UT1-UTC和日長LOD數(shù)據(jù)等，采樣間隔為1天。

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

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

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

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

2 極移模型

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

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

采用插值基礎(chǔ)序列觀測值對式(2)進(jìn)行最小二乘擬合［11］，可知極移趨勢項(xiàng)x(t)包括線性趨勢項(xiàng)A+Bt，錢德勒項(xiàng)(含 pc項(xiàng))和周年項(xiàng)(含 pa項(xiàng))三個部分?？鄢厔蓓?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ù)報是因其具有逼近非線性函數(shù)能力強(qiáng)和網(wǎng)絡(luò)收斂速度快的優(yōu)點(diǎn)。

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

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

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

圖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ù)據(jù)，x(t-1)為最近的過去時刻的數(shù)據(jù)，兩者時間間隔為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ù)長度的選擇對輸出有很大的影響，r的取值范圍為5≤r≤15。隱含層節(jié)點(diǎn)數(shù)s的取值范圍為1≤s≤round(s max)，round是四舍五入取整符號，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 誤差精度評定

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

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

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

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

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

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

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

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

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

由表1可以看出，對于10～120d的預(yù)報跨度，本文方法在極移分量x，y的預(yù)報中MAE小于22mas，RMSE小于27mas。綜合圖6～圖9可知，10～120d的預(yù)報跨度，較文獻(xiàn)［9］方法，本文方法具有明顯優(yōu)勢。對于10～360d的跨度，采用本文方法的極移預(yù)報RMSE值優(yōu)于文獻(xiàn)［9］中的結(jié)果，這說明本文極移預(yù)報方法穩(wěn)定性更好。極移y分量預(yù)報的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ù)報極移得到了很好的精度，這主要得益于插值基礎(chǔ)序列采用了10d間隔的極移數(shù)據(jù)，在不損失數(shù)據(jù)信息的情況下可以選取更多的插值基礎(chǔ)序列數(shù)據(jù)預(yù)報下一時刻數(shù)據(jù)。BP神經(jīng)網(wǎng)絡(luò)模型在中短期(跨度在120d以內(nèi))預(yù)報精度較高，隨著預(yù)報跨度增加，其預(yù)報精度略有下降，但總體效果較好。

6 結(jié)論

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

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

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ù)報［D］．上海:中國科學(xué)院研究生院上海天文臺，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］．測繪通報，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é)報，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．