奇异谱分析在GPS坐标时间序列去噪中的应用

张拥军

（郑州市市政工程勘测设计研究院，河南 郑州 450046）

0 引言

针对海量的GPS 坐标时间序列进行分析是目前研究的热点之一，国内外学者开展了大量的研究。但GPS坐标时间序列中不可避免地会产生噪声，而噪声可能会与GPS 坐标时间序列中的信号产生混叠，淹没有用的信号，因此将GPS 坐标时间序列中的噪声进行剔除具有重要的意义。为了剔除GPS 坐标时间序列中的噪声，不同学者采用了多种不同的模型进行了研究。鲁铁定等人采用改进的经验模式分解对GPS 高程时间序列中的噪声进行剔除研究；黄声亨等人采用小波方法对坐标时间序列数据中的噪声进行剔除；马俊等人采用小波包系数信息熵法去除GPS 坐标时间序列中的有色噪声；焦雄风等人采用抗差卡尔曼滤波（Kalman 滤波）算法对GPS 监测数据去噪。上述方法各有优、缺点及适应的数据，该文采用奇异谱分析方法（Singular Spectrum Analysis，SSA）对GPS 坐标时间序列噪声进行剔除。SSA 方法是主成分分析方法的一种特殊场景，是对一维的时间序列进行分析的主成分分析方法。SSA 方法不受正弦波假定的约束，无须先验信息，具有稳定识别和强化周期信号的优点，能够较好地从时间序列中提取信号，目前已被广泛应用于大地测量领域。

1 SSA 方法与其他方法原理及精度比较指标

1.1 SSA 方法原理

设GPS 坐标时间序列为,,,… , x；x为GPS 时间序列中的观测值；为时间序列长度，且＞2；选用合适的窗口长度，为一个实数，且＜/2，从而构造出滞后矩阵，其矩阵形式如公式（1）所示。

式中：x为,,,… , x中的第个实数值，其他类同。根据公式（1）可知，滞后矩阵副对角线值相等，然后计算滞后矩阵，构造自协方差矩阵，自协方差矩阵=，其矩阵如公式（2）所示。

式中：c 为实数值，其他类同。

由于自协方差矩阵等于滞后矩阵与其转置矩阵相乘，即=，可知为×的矩阵，因此 c为实数值，其他类同。然后利用奇异值分解分别求出自协方差矩阵的特征值 λ与对应的特征向量E。再计算滞后矩阵在特征向量的投影 a，如公式（3）所示。

式中：为窗口长度，是一个正整数；为1～中的数，为滞后矩阵中的数； λ为自协方差矩阵的特征值，通过奇异值分解得到；E为自协方差矩阵特征值对应的特征向量；E为特征向量E中的一项；a 为滞后矩阵在特征向量上的投影。

奇异谱分析的重建则是由时间主成分和时间主分量来重建，重建后的值为ix，具体形式如公式（4）所示。

1.2 小波方法

小波方法的基本原理是利用多个小波基去逼近信号或可积函数。对任意信号或可积函数（），它的小波变换为该信号与小波函数, ()的内积，如公式（6）所示。

因此通过选择合适的小波函数，就能获取信号频域与时域的局部信息。小波去噪主要分为3 个步骤：1） 选择合适的小波函数，确定其分解层次，并分别计算其分解系数。2）对每个分解层次选择恰当的阈值，处理其高频系数，即剔除高频部分的噪声成分。3） 将处理后的高频系数与低频系数进行重构，重构后的信号就是去噪后信号。

1.3 经验模式分解法

经验模式分解（Empirical Mode Decomposition，EMD）是黄锷提出的自适应局部时频分析方法。经验模式分解是将任意一维时间序列分解为不同尺度的固有模态函数，固有模态函数必须同时满足以下2 个条件：1） 局部极值点和过零点的数目必须相等或最多相差1 个。2） 在任意时刻，局部最大值与最小值的包络平均为0。

经验模式分解过程如下：首先，确定时间序列的极大值点与极小值点，并用3 次样条插值法拟合极大值点与极小值点，分别得到其上、下包络线。其次，对上、下包络线进行均值处理，得到均值曲线。再次，用时间序列减去其均值曲线，若满足给定阈值要求，则得到本征函数。若不满足则重复上述过程。最后，用时间序列减去所有本征函数，当残差序列不超过2 个极值点时，经验模式分解完成。

1.4 评价指标

为更加精确地评价奇异谱分析去噪的有效性与优越性，该文采用了3 个评价指标对上述3 种方法，即SSA 方法、小波方法和经验模式分解法进行定量的比较分析。评价指标分别为噪声序列的均方根值；去噪后信号残差序列的均方根值；去噪后的信号与模拟信号的相关系数。

2 模拟试验分析

为验证SSA 去噪效果的有效性与优越性，模拟1 组具有不同周期的正弦波时间序列数据，然后加入噪声，其数据模拟模型y如公式（7）所示。设y由u与e构成，其中u由周期分别为1200 s、500 s 和200 s 的3 种正弦波信号构成，即u=sin（2π/1200）+sin（2π/500）+sin（2π/200），e为服从正态分布（0，0.5）的噪声。

设模拟数据的数量为3 000 个，时间间隔单位为s。为验证SSA 方法的稳健性，该文又分别将小波方法（db8 小波基）、经验模式分解法（EMD）与SSA 方法进行了比较，其结果如图1 所示。

图1中各行子图从上至下的含义如下：第一行为模拟数据；第二行为加入（0，0.52）噪声后的模拟数据；第三行是采用3 种方法去噪后的数据；第四行是采用3 种方法去噪后数据与模拟数据的差值；第五行是3 种方法去噪后模拟数据与加噪后模拟数据的差值。

图1 3 种滤波方法的模拟试验结果

通过对比可以看出，3 种方法均能有效地去除模拟数据中的噪声，其形态整体相似，说明SSA 方法具有良好的去噪能力。为了说明SSA 方法的优越性，该文又分别将值、值与相关系数这3 种指标进行比较，比较结果见表1。

表1 3 种方法去噪效果对比

从表1 中噪声部分的值与相关系数可以看出，其值大小基本相当，说明去除的噪声部分相同，模拟数据与去噪后的数据强相关，再次说明3 种方法都有效地剔除了噪声。但3 种方法的值有较大差异，其中SSA 方法模拟数据与去噪后的数据相差相对较小，而小波方法与EMD 方法则相差较大，说明后2 种方法在去噪过程中剔除了模拟的正弦波信号。因此通过综合对比可知，采用SSA 方法去噪更有优越性。

3 试例分析

该试验数据采用GPS 基准站ZIMM 站的GPS 坐标时间序列，该站有较长时间跨度的时间坐标序列，该文选用跨度为1999—2011 年共12 a 的数据，ZIMM 站坐标时间序列包括趋势项、周期信号以及噪声，较适于验证奇异谱分析方法去噪的效果。为了方便比较，该文对ZIMM 站坐标时间序列进行了预处理。首先，去除ZIMM 站坐标时间序列的趋势项及周期项。其次，采用3 倍中误差方法剔除ZIMM 站坐标时间序列中的中误差。再次，对剔除粗差后的坐标时间序列进行线性插值，得到处理后的坐标时间序列。最后，采用奇异谱分析方法（SSA 方法）进行去噪，其结果如图2 所示。

根据图2 可知，奇异谱分析方法能够有效地去除噪声，并能有效地保留GPS 坐标时间序列中的趋势项及周期项，且图2（b）的噪声序列比较平稳，说明奇异谱分析能够有效地剔除噪声。为了进一步说明奇异谱分析方法去噪的效果，该文又对其进行功率谱分析，其结果如图3 所示。

图2 GPS 坐标时间序列去噪前后的对比及残差序列

根据图3 可知，在低频域，GPS 坐标时间序列基本保持不变，其存在年周期信号与半年周期信号，且信号的波峰较为一致，说明奇异谱分析方法较好地保留了其趋势项与周期项。而在高频域，去噪后的坐标时间序列其功率谱密度明显低于未去噪的坐标时间序列，该现象说明奇异谱分析方法剔除的是GPS 坐标时间序列中的噪声。

图3 GPS 坐标时间序列去噪前后的功率谱对比

4 结论

奇异谱分析是一种高效的时间序列分析方法，该方法是对一维时间序列进行分析的方法，能够较好地从坐标时间序列中提取信息。该文通过模拟数据分析，验证了和小波方法与经验模式分解法相比，该方法去噪效果更好，能够在最大程度地保留坐标时间序列中有用的信号。该文还通过GPS 坐标时间序列实测数据分析，表明该方法能够有效地剔除GPS 坐标时间序列中的噪声，且去噪后的坐标时间序列更加稳定和平滑。