随着高精度全球导航卫星系统(global navigation satellite system，GNSS)技术的发展，在连续运行的GPS参考站上已经积累了近20年的数据，形成了充足的GPS坐标时间序列。GPS台站的时间序列中包含着地壳、地表运动信息，主要分为构造信息和非构造信息^[1]。数据处理中需要对非构造运动信息进行分类剔除。非构造运行信息如海潮运动、天线相位中心、地表大气负荷、地表水迁移等，在GPS时序中还包含着其他一些信息，除纯净的信号之外的信息被认为是噪声。许多研究表明，GPS非构造信号残留有与时空相关的信号，国内外已有研究表明这种与时空相关的信号为共模误差，主要呈现出白噪声、有色噪声特性^[1]，它是引起GPS坐标时间序列中周期变化的主要因素。

本文将利用相关性分析的方法提取GPS时序中的共模噪声。通过主成分分析方法对GPS序列的残差进行主分量提取，并在时序分析中将其剔除，从而得到纯净的GPS时间序列。

1 GPS台站时间序列分析概述 1.1 GPS时间序列信息分类

GPS时间序列里所蕴涵的信息可以分为两大类：构造运动信息和非构造运动信息。

1) 构造运动信息主要包括：构造运动所引起的地球板块运动、地表沉降等，可以表现在时间序列的N、E、U方向体现出来。

2) 非构造运动信息划分为两类：①GPS观测技术类误差，包括GPS卫星轨道误差、基线解算模型、策略等。此类误差又可分为系统误差和随机误差两种^[2]。其中绝大部分的系统误差、随机误差可以通过建立相应的误差改正模型(理论模型、经验模型)，利用求差法或选择较好的硬件等进行消除或减弱，如对海潮运动、天线相位中心、对流层和电离层等通过建模进行改正。但是地表大气负荷、地表水迁移、非潮汐负荷、土壤水分负荷等改正仍包含在GPS时间序列中。②与时空相关的噪声(有色噪声等)即共模噪声(common mode error, CME) ^[3]。由于共模误差的存在，使得GPS台站坐标时间序列解的存在较大的不确定性，导致台站坐标、速度存在有偏估计^[1]。这个误差可以通过估计测试站点分布相关性趋势和周年周期信号来解算消除^[3]。

1.2 主成分分析法

主成分分析(principal component analysis, PCA)是一种非参数的正交分解数据处理方法。它把随时间变化的台站网时间序列分解成时间域的主分量和空间域的特征向量^[4]。

主成分分析的统计过程基于一个完整的观察矩阵，观察矩阵是由残差时间序列r(t)构造。首先，假设时间序列是规则和连续的，用n个站组成的GPS站网构造观测矩阵R(t_i, r_j)(i=1, 2, …, m和j=1, 2, …, n), 同时分别为每个地心分量(N、E、U方向)创建矩阵。

因此，基于观察矩阵的元素协方差的每个元素(n×n)对称矩阵B(t_i, r_j)定义为^[5]：

$ b_{i j}=\frac{1}{m-1} \sum\limits_{k=1}^{m} R\left(t_{k}, r_{i}\right) R\left(t_{k}, r_{j}\right) $

(1)

PCA通过对数据协方差矩阵B进行特征值分解为V^T(V^T是一个有正交的行的(n×n)矩阵)和具有k个非零对角特征值{λ_k}(k≥n)的Λ矩阵:

$ \boldsymbol{B}=\boldsymbol{V} \boldsymbol{\Lambda} \boldsymbol{V}^{\mathrm{T}} $

(2)

需要特征值和特征向量矩阵中对应列的降阶来定义主成分：

$ a_{k}\left(t_{i}\right)=\sum\limits_{j=1}^{n} R\left(t_{i}, r_{j}\right) v_{k}\left(r_{j}\right) $

(3)

式中，a_k(t_i)为R矩阵的k次主成分；v_k这是一个媒介。PCA是去除CME的有效算法，所得到的正交分量(主分量)按变化幅度从大到小排序。对于主分量在每个组成部分所占的方差方面是否显著，采用了文献[6]提出的CME定义来选择主分量，主分量超过一半的对应归一化响应的值大于25%。使用这种方法计算CME表达式为^[7]：

$ \mathrm{CME}_{j}\left(t_{i}\right)=\sum\limits_{k=1}^{p} a_{k}\left(t_{i}\right) v_{k}\left(r_{j}\right) $

(4)

PCA方法能够识别出共模误差，提高CME提取的可靠性^[8]，并顾及形式坐标误差。

2 试验与结果分析 2.1 试验数据与方案

本实验数据收集2010-01-01-2016-01-01跨度6年的台站X/Y/Z的时间序列。台站分布中国境内和周边的19个IGS站，包括：BJFS、URUM、WUHN、ULAB、TWTF、SUWN、SHAO、SELE、POL2、PIMO、YSSK、CHAN、GUAO、IISC、IRKT、KIT3、KHAJ、NTUS、LHAZ。每个站点的数据采集均为1 d，时间序列连续(对于少数缺失数据采用插值方法进行拟合)。

其他资料收集包括：全球大气资料、全球积雪资料、全球土壤湿度资料、非潮汐改正资料等，时间跨度为2010-2016年。

首先对时间序列进行去除阶跃项、趋势项的改正，只剩下周期项(年周期和半年周期)；然后计算残差，对坐标残差进行去除周期项分析，同时加入地壳质量负荷改正，将得到的相对纯净的时间序列数据进行主成分分析；最终使用经PCA滤波之后的纯净的GPS坐标序列来进行N、E、U方向的速度估计。

2.2 谱分析

对各台站时间序列进行频谱分析，可以看出各个台站的频谱图很相似，如图 1所示。从功率谱可明显看出在各站的时间序列中存在着周期性信息^[9]。最明显的就是周期为0.002 5左右的变化信息，初步分析为年变化周期。

2.3 计算各项地表负荷改正

分别计算每个GPS台站的大气质量负荷、非潮汐海洋负荷、积雪、土壤和水分负荷改正。负荷改正均为N、E、U 3个方向的改正量，以BJFS、SHAO站为例，负荷改正量如图 2、图 3所示。

2.4 PCA滤波分析

对原始数据进行PCA滤波分析之前，需要首先纠正原始坐标序列的趋势速度项、阶跃项、指数和对数衰减项等信号，在坐标时间序列中保留各类周期项。在进行各类纠正之后，得到GPS坐标的残差时间序列。

对GPS坐标的残余坐标残差进行时间序列分析，分析结果表明：坐标残差时间序列的3个坐标方向的振幅波动均比较大，说明了坐标方向的中误差相对较大；同时残差序列中存在明显的周期变化，其中U方向尤为明显，主要源自共模误差的影响^[10]。

为了对共模误差进行分离，在滤波过程中首先由PCA法分别计算得到前12个主成分对3个坐标方向(E、N、U)的空间响应，结果见表 1。

采用PCA模型得到的第1主成分(cpt1)对高程U方向的贡献率分别为69.29%、62.15%、70.92%，第2主成分(cpt2)对U方向的贡献率分别为11.59%、13.73%、9.53%，第3主成分(cpt3)对U方向的贡献率为6.28%、10.53%、4.41%。前3个主成分累计贡献率累计为87.16%、86.41%、84.86%，所有第1、2、3主成分的累计贡献率综合了大部分的空间响应(大于85%)。

表 1统计结果表明可以采用PCA方法得到的前3个主成分作为该区域的共性误差来计算，对坐标序列进行PCA滤波。通过PCA法进行滤波处理后，可以剔除共模误差信号，从而获得更为纯净的坐标时间序列。

2.5 PCA滤波结果分析

1) 前后振幅比较和相位大小比较：本实验对19个台站的GPS坐标序列采用地表负荷、PCA滤波前后的周年项的振幅大小进行比较和分析，结果表明：进行地表负荷改正和剔除共模误差实现空间滤波之后，序列周年运动的振幅有所减小。滤波前后年变化振幅减少情况如表 2所示。

从表 2可以看出平均年振幅变化减少最小为3.89%；最大减少了56.42%；平均减少了24%。

2) PCA滤波前后时间序列图比较，以BJFS、LHAZ为例其时序图如图 4、图 5所示。

2.6 坐标和速度场解算

时间序列最终目的通过长期时间序列数据分析得到准确的IGS站点在某个参考历元的坐标和速度场。本次试验将得到的2010.000历元的坐标和速度场与IGS网站公布的ITRF2010框架的坐标和速度进行了比较，结果表 3所示。

从表 3可以看出：采用滤波后的时间序列计算站点参考历元的坐标与IGS成果大部分符合精度在1 cm之内，解算的速度值跟IGS成果符合精度在2 mm之内。说明滤波效果很好。

3 结束语

在高精度GPS数据处理过程中，一般加入的改正模型不足以消除结果中的噪音和误差，所以需要用滤波的方法来探测噪音，从而得到纯净的坐标序列，使结果更加可靠。本文试验证明，针对复杂的时间序列，实施负荷改正和主成分分析等区域滤波后，台站的坐标时间序列趋于平稳，波动振幅有所减小，从而使得站点坐标计算精度提高，坐标序列的稳健性得到增强；同时PCA滤波后残差序列的周期性也有明显减弱，说明共模误差能解释GPS残余序列呈现周期信号的部分原因。