GPS坐标时间序列周期性分析

许焕盛

(广东省地质局第二地质大队，广东 汕头 515000)

山东省于2007年启动了山东省卫星定位连续运行综合应用服务系统(Shan Dong Continuously Operating Reference Station,SDCORS)项目，积累的观测数据不断丰富，为研究区域形变、地球板块运动和地球自转等提供重要的数据支撑[1]。研究GPS坐标时间序列的非线性变化，可获得基准站的位置变化及三维速度场，为工程应用基准站的选择、基准站的速度及数据分析提供参考[2-3]。相关研究表明：GPS坐标时间序列不仅具有线性变化趋势，还具有年、半年周期等的非线性变化特征，其中以U方向表征的最为明显[4-9]。本文以12个SDCORS站连续4年的观测序列为研究对象，计算坐标时间序列的功率谱，通过功率谱图分析其周期项，并将周期项特征进行统计，这对于分析基准站的稳定性及位置变化具有重要意义。

1 数据处理策略

1.1 数据预处理

选取的各基准站的基本信息如表1所示。连续均匀的坐标时间序列是数据分析的基础，因此在数据分析之前需对数据进行预处理，包括剔除孤立值与插值、阶跃项改正、趋势项拟合及去除趋势项等步骤。

表1 站点基本信息

1.1.1 孤立值的剔除与数据补全

由于外界观测条件、多路径效应及传输信号干扰等因素的影响，坐标时间序列中不可避免含有孤立值，也称粗差，这些粗差的存在会造成数据分析结果的偏差。因此，在数据进行周期分析和求站点速度之前，应将其剔除。GPS坐标时间序列由于基站的改造、观测条件不充分等干扰因素，导致在某天或某段时间出现部分数据缺失或不完整的现象，且在剔除粗差的步骤中，部分点被剔除，需要将这两部分的数据进行拟合补齐。

本文采用3倍中误差为判别原则，简称“3S”方法，探测坐标时间序列中的粗差，并将其剔除[9]。“3S”基本原理如下：

(1)

则认为Sa为粗差，需将其从观测序列中剔除。通常采用迭代的方法剔除坐标时间序列含有的粗差。

对于缺失和剔除粗差的数据，经过多次实验与分析，采用三次多项式插值的方法对数据进行补全，保证数据的连续性和完整性，为后续的研究和分析奠定基础。

1.1.2 阶跃项的修正

由于接收机故障、天线更换、地震等原因，GPS坐标时间序列在某时段出现阶跃或称中断，会导致测站速度估计有偏，需要对其进行修正。具体做法是从坐标突变的位置开始，往后的坐标时间序列统一加上阶跃值进行修正。由于sdrc基准站U方向的坐标时间序列中含有明显的坐标阶跃，本文就以sdrc为例对阶跃项的修正作以说明。图1所示为sdrc站U方向的原始坐标序列，从图可以看出存在两处坐标突变的地方，分别在第307天和第775天。图2所示为sdrc阶跃修正后的坐标时间序列，由图可知，经过阶跃修正后的坐标时间序列在66.42 m附近波动，说明具有良好的连续性与一致性。

图1 sdrc原始坐标时间序列

图2 sdrc阶跃修正后的坐标时间序列

1.1.3 趋势项拟合

在进行谱分析时要求坐标时间序列为零均值，即没有线性趋势项。为此，需要对坐标时间序列进行线性拟合，其原理如下：

yt=β0+β1t+xt

(2)

式中，β0、β1为线性趋势拟合的常数项和一次项系数。利用该模型计算得出各个时间点的拟合值，然后将对应时间点的坐标时间序列值减去拟合值，得到残差时间序列，即

(3)

图3所示为sdrc基准站三个方向经过粗差剔除、数据插值和阶跃修正后的坐标时间序列及线性趋势拟合。由图可知：sdrc站的N、E方向具有明显的线性变化，U方向不仅有线性变化，还表现出周期性变化，因此U方向的运动是十分复杂的；U方向数据的离散程度要比N、E方向大，这主要是因为U方向的精度要低于其它两个方向。在进行功率谱分析时要求数据零均值，因此需要将线性趋势项从坐标时间序列中去除，以满足功率谱分析的条件。图4所示为sdrc去除线性趋势项的坐标时间序列，由图可知：坐标时间序列经过上述步骤处理后，变化量在零的上下进行波动，说明具备了零均值的特性。同时图4的坐标时间序列表现出明显的季节性变化，说明存在一定的周期项，这就需要对坐标时间序列进行功率谱估计，探测潜在的周期。

图3 sdrc坐标时间序列及拟合直线

图4 sdrc残差序列

1.2 坐标时间序列的功率谱估计

本文采用传统功率谱估计——直接法，探测坐标时间序列中存在的潜在周期，基本思想为：直接对坐标时间序列作离散傅里叶变换，计算求得坐标时间序列的幅值，然后将幅值的平方除以时间长度作为功率谱密度，得到功率谱密度与周期之间的关系图，也称之为周期图法。周期图的峰值则为坐标时间序列可能存在的周期，且能量越高，对应的周期越明显[10]。其优点是计算简单，而且不依赖于相关函数。基本原理如下：

设x(n)为坐标时间序列，其长度为N的离散傅里叶变换为：

(4)

则功率谱密度定义为式(5)：

(5)

2 坐标时间序列周期性分析

坐标时间序列经过数据预处理得到部分基准站3个方向的残差序列及各方向的功率谱图，如图5所示。

图5 残差序列与功率谱图

图5所示(a)(c)(e)(j)为各站的残差序列图，横坐标为年份，纵坐标为各方向的振幅，单位为mm；(b)(d)(f)(h)为各站的功率谱图，横坐标为周期，单位为a，纵坐标为单位周期内的能量，单位为mm2/a。由图5可知：各站点在N方向的能量较小，均小于5×107，说明周期项不明显。其中sdrc站半年周期项能量最大，为4.963×107，说明sdrc站半年周期运动较明显，其他几个站半年周期项的能量差别不大。另外jimo站、sdrc站、rizh站还存在2年的周期项，cash站、jimo站、rizh站可能还存在更长的周期项。E方向上，各站半年周期项对应的能量突出，功率谱处于8～10×107之间，能量值相差不大，说明各站的半年周期项运动在E方向明显。U方向上，各基准站2年周期项的能量最高，其中sdrc站能量最小，为8.565×107，其他三个站的量级为108。

对其它基准站的功率谱图进行分析，统计各基准站的周期项及功率谱如表2所示，表中功率谱采用科学计数法，量级为107。

表2 基准站周期项及功率谱

由表2可知：

(1)基准站在N、E、U方向均表现出一定的周期性。N方向的周期运动不明显；E方向主要以0.5年周期项运动最为突出，且能量值差别不大；U方向存在2年周期项，绝大部分基准站存在1年周期项，且2年周期项的能量高于1年周期项的能量，说明U方向的2年周期运动较为明显。

(2)坐标时间序列中并不是严格存在半年周期项和年周期项。从功率谱可以看出，能量最强对应的周期并不是严格的半年周期项、年周期项等周期，而是在周期项的附近，称为类半年周期项、类年周期项等。

(3)坐标时间序列功率谱在E方向的0.5年周期项的能量普遍高于N方向的，说明E方向的0.5年周期性运动比N方向的要明显。

(4)N方向整体周期项不明显，且大多数周期项不明显的基准站具有更长周期项的趋势，具体原因有待分析。

3 结 论

本文对SDCORS站坐标时间序列的周期特征进行了分析与统计，得到：SDCORS站坐标时间序列存在周期项，且N、E方向年、半年周期项明显，U方向2年周期项明显，功率谱能量明显高于其它两个方向，这是由于U方向的影响因素较多。在后续的研究中，将对更长时间跨度的坐标时间序列进行周期特征的研究与分析，探讨周期项的影响因素，为得到精确的噪声模型提供数据基础。