基于快速傅里叶变换法的地球自转参数周期性研究

郭忠臣 ，邹 慧

1.宿州学院环境与测绘工程学院，宿州,234000；2.河南省地图院，郑州,450003

地球自转参数(Earth Rotation Parameters,ERP)包括极移(PMX、PMY)和日长变化(LOD)，是地球参考框架(ITRF)和天球参考框架(ICRF)之间相互转换必不可少的参数[1-2]。当前有多种空间大地测量技术(VLBI、GNSS、SLR等)可精确测定ERP，但由于数据处理过程复杂，且不能实时得到高精度的ERP测定结果[3-5]。考虑到ERP在现代空间导航中的作用日渐显著，高精度ERP预报方法也成为当今亟待解决的问题之一[6]。对ERP预报需对其时间序列的有效信息进行提取，傅里叶变换因有诸多优势，在信号处理、物理学等领域已得到广泛应用，本文采用快速傅里叶变换(Fast Fourier Transform,FFT)法对ERP序列的周期性进行分析，为后续ERP预报提供基础。

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

FFT是在离散傅里叶变换(Discrete Fourier Transform,DFT)基础上依据DFT的奇、偶、虚、实等特性发展而来的一种快速算法，其工作原理是将长序列的DFT依据其内在的周期性和对称性分解成许多较短序列的DFT之和[7-8]。

一组长度为N的序列x(n)，其DFT变换为：

(1)

因DFT的计算量与N2成正比，故当序列较长时，为了减少DFT的计算量，可考虑将长序列x(n)的DFT分解成若干短序列的DFT之和。考虑到式(2)中WNnk的周期性、对称性和可约性，并将序列x(n)按照式(3)分成两组，则可推出由两个序列长度为N/2的DFT来表示的式(1)，结果见式(4)。

(2)

(3)

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

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

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

(4)

当N比较大时，对一次分解后的两个序列作DFT,较为复杂，故可对该序列进行多次分解,以简化运算。由于每步均按照输入序列属于偶数还是奇数进行，故称该方法为按时间抽取的FFT算法。

2 实验分析

本文采用的数据为IERS提供的EOP C04序列(http://hpiers.obspm.fr/eoppc/eop/)，采样间隔为1 d，时间跨度取1962/01/01至2015/01/01，图1为IERS发布的原始数据。

由图1可判断出极移序列(PMX、PMY)和LOD序列具有明显的周期性和趋势项，不考虑除周期项与趋势项外的较小残差项，使用Matlab工具箱对极移序列的趋势项进行大致拟合，其余部分均作周期项以进行频谱分析，拟合结果见图2。

图1 ERP时间数据

图2 极移周期项分离结果

考虑到日长变化受固体地球潮汐项影响较大，对LOD预报时首先要对固体地球潮汐项进行改正。1981年，Yoder给出了计算固体潮的方法[9]，后被IERS采用，本文所用改正模型见IERS Convention 2010，改正后的LOD记为LODR。图3给出了固体地球潮汐项对LOD的影响值，从图中可以看出，固体地球潮汐项对LOD的影响不能忽略。

图3 固体地球潮汐改正项对LOD的影响值 图4 扣除固体地球潮汐改正项后的LOD

图5 LODR周期项分离结果

与极移类似，不考虑除周期项与趋势项外的较小残差项，使用Matlab工具箱对LODR序列的趋势项进行大致拟合，其余部分均作周期项以进行频谱分析，拟合结果见图5。

利用Matlab编写FFT程序对扣除趋势项后的极移序列和LOD序列进行处理，处理结果见图6和图7，通过分析可知，极移受1.2周年项(Chandler项)和周年项的影响最大，LODR受周年项和半周年项的影响最大。

为进一步分析ERP序列的特性，依据其性质，可分别构建最小二乘模型对极移和LODR的周期项和趋势项进行分离，对极移和LODR构建的模型见公式(5)和(6)。

图6 极移功率谱图

图7 LOD功率谱图

极移最小二乘模型：

其中,ai、bi为趋势项的拟合系数，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为趋势项的拟合系数，c、d为各周期项的振幅，P1、P2分别为周年项和半周年项的周期，ωi为各周期项的相位，t为时间。

选择2005/01/01至2015/01/01年间10年数据代入极移和LODR最小二乘模型中迭代解算趋势项和周期项系数，并利用解算结果对原始序列的趋势项和周期项进行拟合，拟合结果见图8。由图8可知，对于极移，残差项的影响在-0.05 as～+0.05 as之间；对于LODR，残差项的影响在-0.5 ms～+0.5 ms之间。与原始序列相比，残差项的影响不容忽视，后续预报时需再对残差项进行建模处理。

3 结束语

本文利用Matlab编程实现了ERP序列周期性的提取，实验结果表明：

(1)极移主要受到1.2周年项和周年项的影响比较大，LOD主要受到周年项和半周年项的影响比较大。

图8 地球自转参数拟合结果

(2)ERP序列主要受到周期项和线性趋势项的影响，但扣除趋势项和周期项后，残差项影响依然较大，需构建其他模型对残差项进行分析。

本文仅对ERP序列的周期性进行了分析，并依据其特性构造最小二乘模型对周期项和趋势项进行提取，但通过实验可知，ERP序列中同样含有影响较大的残差项，下一步会构建函数模型，对剩余的残差序列进行分析，并对ERP序列进行预报。

参考文献:

[1]徐天河,王潜心,于素梅,等.利用区域网GPS/BDS数据确定地球自转参数[J].导航定位学报,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]姚宜斌,岳顺强,陈鹏.一种适用于极移预报的附加误差修正的LS+AR新模型[J].中国科学:地球科学,2013,43(4):665-676

[4]魏二虎,常亮,刘经南.我国进行激光测月的研究[J].测绘信息与工程,2006,31(3):1-3

[5]张志斌,王广利,刘祥,等.中国VLBI网观测地球定向参数能力分析[J].武汉大学学报:信息科学版,2013,38(8):911-915

[6]魏二虎,万丽华,金双根,等.联合GNSS和SLR观测对地球自转参数的解算与分析[J].武汉大学学报:信息科学版,2014,39(5):581-585

[7]张登奇,李宏民,李丹.按时间抽取的基2 FFT算法分析及MATLAB实现[J].电子技术,2011,38(2):75-77

[8]胡林林,邓超兵,付龙.基于C的FFT算法在波长扫描干涉中的应用[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