基于奇异谱分析的GPS时间序列粗差探测及缺失数据插值 - 图文

基于奇异谱分析的GPS时间序列粗差探测及缺失数据插值

研究

王晓明1 成英燕1 蒋志浩2 王福丽1 陈 醒1

（1.中国测绘科学研究院，北京 100830；2.国家基础地理信息中心，北京 100048）

摘 要 由于传统方法在进行GPS时间序列粗差探测和数据插值的方法存在一定的缺陷，本文提出了基于奇异谱分析的SSA-IQR（Singular Spectrum Analysis-Inter Quartile Range）粗差探测方法，同时引入缺失数据插值方法SSA-M（Singular Spectrum Analysis for Missing Data）。两种方法都是基于改进的奇异谱分析可以通过自适应滤波的方法从不完整坐标序列中提取主要特性成分，并利用其建立模型完成粗差位置判断和缺失部分插值的基本原理。本文利用模拟数据对方法特性及效果进行验证，并利用该方法对国内部分CORS站完成了不完整坐标序列的粗差剔除和数据插值。 关键词 大地测量学；粗差探测；数据插值；奇异谱分析

Study on Gross Error Detection and Interpolation for GPS Time

Series Based on SSA

Xm Wang,Yingyan Cheng,Zhihao Jiang,Fuli Wang,Xing Chen

（1.Chinese Academy of surveying and Mapping, Beijing 100830；

2.National Geomantics Center of China,Beijing 100048）

Abstract：As the traditional method of gross error detection and data interpolation for GPS time series has some flaws and limitations, the SSA-IQR （Singular Spectrum Analysis-Inter Quartile Range） is used to detect gross error and SSA -M （Singular Spectrum Analysis for Missing Data） is used to fill the gaps of several data sets in this paper.SSA-M is a self-consistent method to infer missing data from GPS data series and extract the relevant empirical orthogonal functions.Since the filling and error detection method does not need a priori information about the data structure, the method is self-consistent and parameter free. Based on the theories of SSA and IQR criteria，the SSA-IQR method for gross error detection in observation data series is also studied. In this paper, the experiment is performed to identify the reliability of the methods,and the incomplete data of IGS stations in china is reconstructed.

Keywords: geodesy; gross error detection; interpolation; singular spectrum analysis

1 前言

近几年来，连续GPS观测在包括变形监测等多个领域中都开始发挥重要应用[1-3]。目前对于站点运动特性的研究主要是借助功率谱等方法分析站点的各类周期运动，而GPS连续站的粗差值和不均匀采样也使得这些方法的使用存在一定的障碍。本文认为传统数据插值方法存在一定的复杂性和不可靠性，因而采用SSA-M插值方法对模拟数据和GPS测站数据进行了方法验证，并利用模拟数据分析的方法对通常使用的粗差位置探测方法的缺陷进行了说明，提出了基于奇异谱分析的无先验模型粗差探测方法。

奇异谱分析（SSA，Singular Spectrum Analysis）是一种从序列的动力重构出发，并与经验正交函数相联系的一种统计技术[5-7]。对奇异谱分析方法的改进可以使其从不完整序列中

[4]

提取有效信息，继而建立适合数据结构的最佳模型，结合稳健估计粗差探测IQR准则进行粗差定位，将粗差剔除并同缺失数据一起最终完成插值工作。奇异谱分析的过程中涉及窗口嵌入维数M和重构阶次P两个参数，本文中采用交叉验证的方法完成上述两个参数的确定，在保证正确探测粗差和恢复缺失数据主要信息的前提下，最大限度的降低人为噪声的引入量。

2 方法原理介绍

本文提出的的SSA-IQR粗差探测及引入的SSA-M缺失数据插值方法都是在基于传统奇异谱分析的基础上，改变方法中的自协方差矩阵的求解方法，并通过迭代过程进行交叉验证的方法完成最佳数据结构模型的建立，继而根据模型完成粗差修复及缺失数据插值的数据分析方法。

2.1 传统奇异谱分析的基本原理

奇异谱分析是建立在相空间重构基础上的EOF分解，它首先计算原始数据的迟后自协方差矩阵的特征向量空间，然后将原始数据的迟后排列的矩阵向该正交空间进行投影，最后对重构相空间以主成分析的方法完成数据的内在结构分析。篇幅所限，本文不对该方法做详细介绍，具体可参考文献[5]。

2.2 SSA-M方法介绍

Schoellhamer[8]（2001）最早提出的了利用奇异谱分析进行数据插值，但方法存在诸多方面的缺陷。Kondrashov在参考该方法的基础上提出了新的奇异谱插值方法SSA-M[9]。该方法的具体步骤为：（i）对原始数据进行缺失数据位置的标记，并对缺失部分进行补零，对序列做中心化处理；（ii）对该序列采用传统奇异谱分析的方法进行分解，得到PCs及其对应的RCs，选取第一个RC（以下记为R1）中的值替换原序列中缺失位置处的数值，反复循环此过程直至收敛，该循环在整个数据插值过程中被称作内循环步骤；（iii）第一次内循环结束后，通过添加R重建缺失数据，即缺失数据由R和R通过线性叠加的方式得到，接下来的过程重复步骤ii，直至收敛，该过程在整个插值过程中被称作外循环步骤；（iv）在

212整个计算过程中，采用交叉验证的统计方法计算在选取不同嵌入维数和不同重构阶次时的测试误差，以此确定最佳的嵌入维数M和最佳的重构阶次P。

2.3 SSA-IQR方法介绍

目前对于坐标时间序列粗差探测的主要方法为，首先利用公式（假定线性趋势项、年周期和半年周期）通过最小二乘的方法获取站点的残差序列，然后利用三倍中误差等粗差准则进行异常值的剔除。这种方法对粗差的探测存在一定的缺陷：第一，已有学者的研究证明GPS坐标时间序列中含有部分有色噪声[10]，而通过最小二乘解算假定公式获取残差的过程是建立在残差为白噪声的假设之上的，这是不合理的；其次，对我国国内CORS站的分析表明，GPS坐标时间序列中除周年和半年周期项，还可能存在其他高频周期项[11]，同时对周期项做先验正弦波假设也不够严谨。由此本文认为先验假定公式同实际的站点运动形式未必完全一致，因此其解算的残差的性质还与先验模型与实际运动形式差异有关。

SSA在分析数据结构时的无先验假设特性及抗粗差特性使其在结合IQR准则进行粗差识别时更具有效性。该方法具体步骤为，先采用同SSA-M中一样的迭代-交叉验证方法完成最佳M值和P值的选择，以此为模型参数完成数据重构建模，然后从原始数据列中减去重构模型值即得到残差数据序列，最后利用IQR准则进行粗差探测，剔除残差并标记为数据缺失，利用SSA-M进行数据插值，然后重新检测，若仍存在粗差则重复SSA-M和SSA-IQR，直至不再检测到粗差为止。

3 模拟实验

本文结合GPS站点坐标序列的特性对多种模拟数据进行插值效果验证。首先构造时间序列x，序列中包含趋势项和三个周期项

tx(ti)?ati/yr?b1sin(2?ti/yr1)?b2cos(2?ti/yr1)?c1sin(2?ti/yr2?w1)?c2cos(2?ti/yr2?w1)?d1sin(2?ti/yr3?w2)?d2sin(2?ti/yr3?w2) （1）

本次试验中，依据对我国CORS站点的运行特性研究经验，将（1）式中的代表速率及振幅大小的参数a、b1、b2、c1、c2、d1、d2的值分别取为2、1、1.2、1.2、0.8、0.5、0.7，将表示周期的参数yr、yr1、yr2、yr3的值分别取为100、100、50和30，以代表年、半年和季节性周期，w和w随机取值。

利用AR模型构造有色噪声yt,其构造方程为式（2），其中f=0.1，f=0.2，wt为白噪

1212声，其方差为0.5，均值为0。最后合成得到模拟数据序列总长度为400。考虑到GPS在观测过程中会遇到地震的影响，因此在模拟数据中添加阶跃项，阶跃位置为序列号170处，纵坐标阶跃为2，粗差及缺失值位置不变。由于震后的站点的运动特性可能变化，因此阶跃后的模拟数据公式（1）中的参数a、b、b、c、c、d、d分别改为2.05、0.9、1、1、0.6、0.7、0.6，阶跃前后的噪声添加方式不变。

121212yt（i）=f1rt（i-1）+f2rt（i-2）+wt （2）

原始数据和最终结果如图1所示 ，其中图1.1为原始数据序列y，图1.2是表示在原始序列y的基础上增加模拟粗差值和数据缺失后的序列，本次试验中模拟粗差值10个，占总序列长度的2.5%，缺失数据为107个，占总序列长度的26.75%。插值结果见1.3，从图中可以看出利用SSA-M可以很好的还原缺失值，而且这种还原得到的是滤掉噪声后的平滑

tt值，另外从图1.3中也可以看出，SSA-IQR可以很好的识别序列中的粗差并结合SSA-M完成粗差的修复，最终得到干净完整的数据序列。图1.4中可以看出插值残差与模拟有色噪声有较好的一致性，且两者方差值接近（残差值方差为0.4，有色噪声方差值为0.5），所以我们可以认为SSA-M插值方法虽然不可以完全避免人为噪声的引入，但可以将这种影响尽量最小化。

图1 模拟地震影响数据的粗差修复及缺失数据插值结果

从图1.3中可以看出，当我们将地震时刻假定在序列号170处，并设置160~180区间存在数据缺失时，序列170~180处的数据插值效果并不理想，但除此之外其他各处的粗差修复和缺失数据插值都得到了理想的结果，从图1.4中也可以看出除170~180区间内残差较大外，其余各处的模型残差与序列中的有色噪声大致相等。

4 GPS连续运行站插值模型及粗差探测试验

本文采用的数据是由IGS提供的站点在ITRF05框架下国内6个IGS站的坐标序列，对

所有站点采用表一所示方法模拟缺失数据，以方便对方法结果进行验证。

表1 模拟缺失数据方式 起始时间 终止时间 （年 年积日） （年 年积日） 2002 47 2004 167 2005 301 2006 306 2002 147 2005 1 2005 331 2006 336 缺失天数 （天） 100 200 30 30 最大连续缺失天数 （天） 100 200 30 30 缺失段 1 2 3 4 5 缺失方式 连续 连续 连续 连续 随机抽取 2002 192 2003 232 84 15 因插值结果基本一致，本文以BJFS站为例进行说明。BJFS站采用的数据采集时间是自1999年294日至2008年319日之间，跨度为3313天，其中缺失170天的数据，缺失率为5.13%，数据连续缺失最多为35天。为进一步验证方法的精度和可靠性，在原数据序列中模拟5个缺失段，详细信息见表1。最终序列中的总缺失数据为613天，占总数据的18.5%。

图2~图4为利用奇异谱分析得到的数据插值模型与原始数据序列的比对图，横轴代表自1999年294日为起始历元的时间，单位为天，蓝色散点为站点的原始序列，红色圆圈表示利用SSA-IQR探测得到的站点粗差位置，红色线为利用SSA-M得到的插值模型。

图2 BJFS站E方向粗差修复及插值模型

图3 BJFS站N方向粗差修复及插值模型

图4 BJFS站U方向粗差修复及插值模型

图5为原始数据序列与插值模型之间的插值，横轴的意义与图2相同，绿色部分为模拟缺失数据处的残差，红色圆圈代表SSA-IQR识别的粗差点。从图中可以看出插值结果较稳定性及精度都较好。