郭忠臣 ，鄒 慧

1.宿州學(xué)院環(huán)境與測繪工程學(xué)院，宿州,234000；2.河南省地圖院，鄭州,450003

地球自轉(zhuǎn)參數(shù)(Earth Rotation Parameters,ERP)包括極移(PMX、PMY)和日長變化(LOD)，是地球參考框架(ITRF)和天球參考框架(ICRF)之間相互轉(zhuǎn)換必不可少的參數(shù)[1-2]。當(dāng)前有多種空間大地測量技術(shù)(VLBI、GNSS、SLR等)可精確測定ERP，但由于數(shù)據(jù)處理過程復(fù)雜，且不能實(shí)時得到高精度的ERP測定結(jié)果[3-5]?？紤]到ERP在現(xiàn)代空間導(dǎo)航中的作用日漸顯著，高精度ERP預(yù)報方法也成為當(dāng)今亟待解決的問題之一[6]。對ERP預(yù)報需對其時間序列的有效信息進(jìn)行提取，傅里葉變換因有諸多優(yōu)勢，在信號處理、物理學(xué)等領(lǐng)域已得到廣泛應(yīng)用，本文采用快速傅里葉變換(Fast Fourier Transform,FFT)法對ERP序列的周期性進(jìn)行分析，為后續(xù)ERP預(yù)報提供基礎(chǔ)。

1 基2時域抽取FFT算法原理

FFT是在離散傅里葉變換(Discrete Fourier Transform,DFT)基礎(chǔ)上依據(jù)DFT的奇、偶、虛、實(shí)等特性發(fā)展而來的一種快速算法，其工作原理是將長序列的DFT依據(jù)其內(nèi)在的周期性和對稱性分解成許多較短序列的DFT之和[7-8]。

一組長度為N的序列x(n)，其DFT變換為：

(1)

因DFT的計算量與N2成正比，故當(dāng)序列較長時，為了減少DFT的計算量，可考慮將長序列x(n)的DFT分解成若干短序列的DFT之和?？紤]到式(2)中WNnk的周期性、對稱性和可約性，并將序列x(n)按照式(3)分成兩組，則可推出由兩個序列長度為N/2的DFT來表示的式(1)，結(jié)果見式(4)。

(2)

(3)

其中,u=0,1,…,N/2-1

X(k)=DFT[x(n)]

=DFT[x(2u)]+DFT[x(2u+1)]

(4)

當(dāng)N比較大時，對一次分解后的兩個序列作DFT,較為復(fù)雜，故可對該序列進(jìn)行多次分解,以簡化運(yùn)算。由于每步均按照輸入序列屬于偶數(shù)還是奇數(shù)進(jìn)行，故稱該方法為按時間抽取的FFT算法。

2 實(shí)驗(yàn)分析

本文采用的數(shù)據(jù)為IERS提供的EOP C04序列(http://hpiers.obspm.fr/eoppc/eop/)，采樣間隔為1 d，時間跨度取1962/01/01至2015/01/01，圖1為IERS發(fā)布的原始數(shù)據(jù)。

由圖1可判斷出極移序列(PMX、PMY)和LOD序列具有明顯的周期性和趨勢項，不考慮除周期項與趨勢項外的較小殘差項，使用Matlab工具箱對極移序列的趨勢項進(jìn)行大致擬合，其余部分均作周期項以進(jìn)行頻譜分析，擬合結(jié)果見圖2。

圖1 ERP時間數(shù)據(jù)

圖2 極移周期項分離結(jié)果

考慮到日長變化受固體地球潮汐項影響較大，對LOD預(yù)報時首先要對固體地球潮汐項進(jìn)行改正。1981年，Yoder給出了計算固體潮的方法[9]，后被IERS采用，本文所用改正模型見IERS Convention 2010，改正后的LOD記為LODR。圖3給出了固體地球潮汐項對LOD的影響值，從圖中可以看出，固體地球潮汐項對LOD的影響不能忽略。

圖3 固體地球潮汐改正項對LOD的影響值 圖4 扣除固體地球潮汐改正項后的LOD

圖5 LODR周期項分離結(jié)果

與極移類似，不考慮除周期項與趨勢項外的較小殘差項，使用Matlab工具箱對LODR序列的趨勢項進(jìn)行大致擬合，其余部分均作周期項以進(jìn)行頻譜分析，擬合結(jié)果見圖5。

利用Matlab編寫FFT程序?qū)鄢厔蓓椇蟮臉O移序列和LOD序列進(jìn)行處理，處理結(jié)果見圖6和圖7，通過分析可知，極移受1.2周年項(Chandler項)和周年項的影響最大，LODR受周年項和半周年項的影響最大。

為進(jìn)一步分析ERP序列的特性，依據(jù)其性質(zhì)，可分別構(gòu)建最小二乘模型對極移和LODR的周期項和趨勢項進(jìn)行分離，對極移和LODR構(gòu)建的模型見公式(5)和(6)。

圖6 極移功率譜圖

圖7 LOD功率譜圖

極移最小二乘模型：

其中,ai、bi為趨勢項的擬合系數(shù)，ci、di為各周期項的振幅，P1、P2分別為周年項和1.2周年項的周期，ωi為各周期項的相位，t為時間。

LODR最小二乘模型：

L(t)=a+bt+ccos(2×pi×t/P1+ω1)

+dcos(2×pi×t/P2+ω2)

(6)

其中,a、b為趨勢項的擬合系數(shù)，c、d為各周期項的振幅，P1、P2分別為周年項和半周年項的周期，ωi為各周期項的相位，t為時間。

選擇2005/01/01至2015/01/01年間10年數(shù)據(jù)代入極移和LODR最小二乘模型中迭代解算趨勢項和周期項系數(shù)，并利用解算結(jié)果對原始序列的趨勢項和周期項進(jìn)行擬合，擬合結(jié)果見圖8。由圖8可知，對于極移，殘差項的影響在-0.05 as～+0.05 as之間；對于LODR，殘差項的影響在-0.5 ms～+0.5 ms之間。與原始序列相比，殘差項的影響不容忽視，后續(xù)預(yù)報時需再對殘差項進(jìn)行建模處理。

3 結(jié)束語

本文利用Matlab編程實(shí)現(xiàn)了ERP序列周期性的提取，實(shí)驗(yàn)結(jié)果表明：

(1)極移主要受到1.2周年項和周年項的影響比較大，LOD主要受到周年項和半周年項的影響比較大。

圖8 地球自轉(zhuǎn)參數(shù)擬合結(jié)果

(2)ERP序列主要受到周期項和線性趨勢項的影響，但扣除趨勢項和周期項后，殘差項影響依然較大，需構(gòu)建其他模型對殘差項進(jìn)行分析。

本文僅對ERP序列的周期性進(jìn)行了分析，并依據(jù)其特性構(gòu)造最小二乘模型對周期項和趨勢項進(jìn)行提取，但通過實(shí)驗(yàn)可知，ERP序列中同樣含有影響較大的殘差項，下一步會構(gòu)建函數(shù)模型，對剩余的殘差序列進(jìn)行分析，并對ERP序列進(jìn)行預(yù)報。

參考文獻(xiàn):

[1]徐天河,王潛心,于素梅,等.利用區(qū)域網(wǎng)GPS/BDS數(shù)據(jù)確定地球自轉(zhuǎn)參數(shù)[J].導(dǎo)航定位學(xué)報,2015,3(3):13-17

[2]Yao Y B,Yue S Q,Chen P.A new LS+AR model with additional error correction for polar motion forecast[J].Science China Earth Science,2013,56(5):818-828

[3]姚宜斌,岳順強(qiáng),陳鵬.一種適用于極移預(yù)報的附加誤差修正的LS+AR新模型[J].中國科學(xué):地球科學(xué),2013,43(4):665-676

[4]魏二虎,常亮,劉經(jīng)南.我國進(jìn)行激光測月的研究[J].測繪信息與工程,2006,31(3):1-3

[5]張志斌,王廣利,劉祥,等.中國VLBI網(wǎng)觀測地球定向參數(shù)能力分析[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2013,38(8):911-915

[6]魏二虎,萬麗華,金雙根,等.聯(lián)合GNSS和SLR觀測對地球自轉(zhuǎn)參數(shù)的解算與分析[J].武漢大學(xué)學(xué)報:信息科學(xué)版,2014,39(5):581-585

[7]張登奇,李宏民,李丹.按時間抽取的基2 FFT算法分析及MATLAB實(shí)現(xiàn)[J].電子技術(shù),2011,38(2):75-77

[8]胡林林,鄧超兵,付龍.基于C的FFT算法在波長掃描干涉中的應(yīng)用[J].工業(yè)控制計算機(jī),2016,29(12):35-36

[9]Yoder F C,Williams J G,Parke E M.Short period tidal variations of earth rotation[J].Journal of Geophysical Research Solid Earth,1981,86(B2):881-891