0 引言

GNSS观测技术已成为目前主要的形变观测手段，被广泛应用于区域地壳形变研究(Mazzotti et al., 2003; Meade and Hager, 2005)，强地震位移场检测(Melbourne et al., 2002; Freed et al., 2006；Liu et al., 2012)，慢地震位移研究(Szeliga et al., 2004; Melbourne et al., 2005)，全球板块运动模型研究(Prawirodirdjo and Bock, 2004)，以及冰后期回弹分析(Sella et al., 2007)等领域，并发挥了重要的作用.然而，GPS定位受到众多因素的影响，某些在时空相关的公共噪声，如地壳运动、海潮及固体潮、地下水及气候变化、解算策略等诸多因素对区域GPS网的影响，通常称之为共模误差(common-mode error, CME).共模误差的有效提取与剔除对于提高区域连续GPS坐标时间序列精度具有重要的作用(Wdowinski et al., 1997; Nikolaidis, 2002; Dong et al., 2006).

近年来，有关学者针对共模误差展开了大量研究.早在1997年就有学者发现GPS解算的单日坐标残差有一定的相似性, 认为存在共模误差, 提出了CME校正算法(Wdowinski et al., 1997).对共模误差进一步深入研究，发现GNSS观测资料中除了含有白噪声外还有其他噪声，并进一步提出了日坐标残差的带权均值法(Nikolaidis, 2002；黄立人, 2006; 黄立人和符养，2007)、主成分分析法(Principal Component Analysis, PCA)、堆栈法(stacking)和Karhunen-Loeve展开联合滤波三种方法应用在GPS网单天坐标时间序列的滤波处理中，对三种方法的滤波效果进行比较，发现运用PCA进行共模误差的提取效果更优(Dong et al., 2006；胡守超等，2009；蒋志浩等，2010；Shen et al., 2014；马超等，2016).杨博等(2014)利用欧拉滤波法、均值滤波法和多核函数识别共模误差，发现共模误差的空间响应是不完全一致的.上述大部分对于共模误差的识别都是在小于500 km范围的较小空间尺度进行的.然而，Márquez-Azúa和DeMets(2003)根据测站间的距离提出空间权重叠加法进行大约3000 km范围空间尺度共模误差的识别.

多通道奇异谱分析(Multi-channel Singular Spectrum Analysis，MSSA)(Vautard et al., 1992)作为一种数字信号处理技术，在气候学(Wyatt et al., 2012)、测量学(Chen et al., 2013; Zabalza et al., 2014)、海洋学(Kondrashov and Berloff, 2015)等领域都有应用.MSSA可以从未知物理本质的数据中提取出尽可能多的可靠信息，并基于这些信息对数据进行趋势识别、周期判断、降噪处理和数据重构等(吴洪宝和吴蕾, 2005; Shen et al., 2018; Guo et al., 2018). MSSA从时间序列的动力重构出发，具备不受正弦波假定的约束，无需先验信息，具有稳定识别和强化周期信号等优点，特别适合于分析有周期震荡的时间序列数据.

为了有效的减弱或剔除共模误差，本文构建一种使用多通道奇异谱分析提取GNSS坐标时间序列中共模误差的新思路，通过对剔除共模误差前后的时间序列进行对比分析，验证利用MSSA能够有效的提取共模误差，提高GPS坐标时间序列的精度.

1 原理及方法

共模误差的存在使得GNSS测站坐标时间序列的局部变化特征被公共变化特征所掩盖，影响了GNSS坐标时间序列的分析，共模误差的剔除可以有效的提高GNSS坐标时间序列的可靠性.因此，本文将按照图 1给出的流程进行共模误差提取与分析.首先将GNSS坐标时间序列剔除某一历元的坐标后得到相对坐标时间序列；然后利用MSSA从相对坐标时间序列中提取残差序列，并对残差序列做预处理；然后利用MSSA从预处理后的残差时间序列中提取共模误差，并剔除共模误差；最后对剔除共模误差前后的测站坐标时间序列、相对坐标时间序列和残差序列进行对比分析.

1.1 数据预处理方法

本文中的数据预处理主要包括残差序列粗差的探测、剔除和对缺失数据的迭代插值.

文中使用的粗差探测方法为多通道奇异谱粗差探测方法，具体如下：首先将原始GNSS坐标时间序列的某个历元下的坐标扣除得到相对坐标时间序列；由于相对坐标时间序列中存在的已知周期有周年周期和半年周期且本文所用数据的时间分辨率为一天，因此，在提取残差序列时选择已知周期的最小公倍数365为窗口长度(吴洪宝和吴蕾，2005)；然后对相对坐标时间序列进行MSSA的分解，将序列中除了周年周期、半周年周期以及趋势项以外的信号进行重构作为残差；最后对残差序列进行分析并计算其标准差，依照3倍标准差(3σ)准则确定阈值，超过阈值的被认为是粗差并剔除.

对于因剔除粗差造成的数据缺失，或原时间序列中本就存在的数据缺失，利用MSSA迭代插值补全数据.MSSA迭代插值(Schoellhamer, 2001; Beckers and Rixen, 2003; Kondrashov and Ghil, 2006; Kondrashov et al., 2010; Shen et al., 2017, 2018；Guo et al., 2018)的具体方法是：首先对残差序列中的缺失数据用0做补全并标记，对补全后的时间序列做中心化处理；对于补全后的序列做MSSA分解；然后进行加权相关分析(w-correlation)，确定前n阶为残差序列的主成分，残差序列中缺失的数值用第1个重构成分中的值代替，反复重复此过程，根据经验(马超等，2016)，当前后两次插补的数据残差RMS小于0.01 mm时结束此次插补过程；然后继续对第2个重构成分重复上述插补过程，直至第n个重构成分插补完成.

1.2 共模误差的识别方法

GNSS相对坐标时间序列(X(t_i))为

(1)

其中，t_i表示时间；X_趋势(t_i)表示趋势项；X_周年(t_i)表示周年周期项；X_半周(t_i)表示半年周期项；Δ_残(t_i)为测站的残差.

MSSA方法不但可以将序列中的信号进行很好的分离，同时又能顾及不同通道之间的相干性，提取出时间序列的谱性质(Vautard et al., 1992).因此，本文中选用MSSA进行共模误差的提取.假设L维的GNSS相对坐标时间序列为X_lt(l=1, 2…, L, t=1, 2, …N)，其中t为时间序号，N为坐标时间序列的长度，l为空间点号，L为空间点的个数.则可以得到其时迟轨迹矩阵X为

(2)

X矩阵的第i列，也就是第i个状态向量为

(3)

i=0, 1, 2, …, N-M, 共N-M+1个状态向量. X矩阵有L×M行N-M+1列.类似于经验正交函数分解(EOF)，MSSA展开的分量形式为

(4)

这里的M为窗口长度.特征向量E^k是L×M维的，它的分量E_lj^k是第k个特征向量在l通道时间迟后j的分量. E^k称为空间-时间EOF(简称ST-EOF).

时间系数a_i^k是第i个状态向量X_i在E^k上的投影：

(5)

式中a_i^k为第k个空间-时间主成分(ST-PC)，0≤i≤N-M.

E^k是L×M阶分块Toeplitz矩阵T_X的特征向量，T_X是由L×L块矩阵组成：

(6)

其中的某个块矩阵T_ll′是第l通道x_li和第l′通道序列x_l′i之间迟后0，1，…，M-1的迟后协方差矩阵，T_ll′矩阵的第j行第j′列元素的无偏性最好的估计是

(7)

特征值λ^k等于第k个ST-PC序列的方差，为了选取重要的ST-EOF和ST-PC，特征值也是按降序排列.

重建成分是从被分析的对象中提取与某个或某些特征成分相联系的部分.注意矩阵X的结构(2)式和展开式(4)式，把(4)式右边的求和限制在第k个，就得到有第k个ST-EOF和ST-PC重建的x_li的部分，记为x_li^k，公式为

(8)

通过公式(2)至(8)就可以将预处理之后的相对坐标时间序列中的各类信号进行分离，利用公式(1)将分离结果中周年周期信号、半年周期信号和趋势项信号进行剔除得残差序列，按照1.1节的方法对残差序列进行预处理，对预处理后的残差序列利用公式(2)至(8)的方法进行分解，由于预处理后的残差序列中包含着共模误差和测站自身的噪声等，因此利用加权相关分析(w-correlation)方法(Hassani, 2007)对各重构成分(Reconstruction component, RC)之间的相关性进行分析，将具有相同信号特征的RC进行分组，将残差序列中相关性较高的部分进行重构.假设得到的时间序列为Y_i，则任意两个重建时间序列的w-correlation为

(9)

其中，为权重系数，w_k=min(k, M, N-k).ρ_{i, j}^w的绝对值越趋近于1，说明i, j对应部分之间的相关性越大；根据经验，相关性大于60%的即认为两者之间存在相关性.因此，我们将相邻两种重构成分相关性大于60%且与其他重构成分相关性之和小于60%的二者认为是一组信号，将其进行重构，认为其为共模误差.

2 实验结果与分析 2.1 实验数据及预处理

使用的GPS坐标时间序列是由The University of La Rochelle(ULR)利用GAMIT/GLOBK软件解算，将坐标时间序列扣除掉2004.4973历元坐标值得到相对坐标时间序列，相对坐标时间序列可由http://www.sonel.org获得.实验区域及测站分布如图 2所示.时间序列的时间分辨率为一天，各测站坐标时间序列的平均解算精度为7.45 mm，时间序列解算精度符合本文研究要求.本文仅针对测站U方向展开研究；由于MSSA要求序列长度一致，因此，均将序列长度截取为9年(2002—2010年).ULR坐标时间序列的解算策略(Santamaría-Gómez et al., 2012)为：截止高度角设置为10°；历元间隔设置为30 s；选用的对流层延迟模型和对流层映射函数分别为Saastamoinen模型和VMF1模型；对流层参数每一个小时估计一次；潮汐改正模型和太阳光压摄动模型分别选用FES2004模型和BERNE模型.GAMIT采用双差模式处理基线向量，消除了部分共模误差，但是实验GPS网中仍然存在共模误差的影响，如数据处理后残余的海潮、固体潮的影响，地壳运动的影响以及诸多未建模或模型不准确等因素对区域GPS网造成的影响.

时间序列中会存在一些粗差，粗差的存在会对时间序列的分析产生不良影响，所以有必要对粗差进行剔除.按照图 1所示的流程、第2节中粗差探测与剔除的方法对残差序列进行处理，处理结果如表 1所示.由表 1可得，平均剔除17个历元数据，平均剔除率为0.5%；剔除粗差最多的测站为VNDP，剔除33个数据，最大剔除率为0.9%.

由于时间序列中部分数据被剔除且原时间序列也有少数数据缺失，无法直接进行共模误差的提取，因此需要进行插值获得等时间间隔序列.按照如图 1所示的流程、第2节中奇异谱迭代插值的方法对残差序列进行插值.由于篇幅限制，仅给NANO测站残差序列插值后的时间序列，如图 3所示.

2.2 共模误差的提取

共模误差是区域连续GPS网中存在的一种与时空相关的主要误差源，MSSA方法能够顾及不同通道之间的相干性，因此，本文使用MSSA按照图 1所示的流程提取共模误差.通过对残差序列的分析发现，残差序列中存在60天和120天的周期，因此，在提取共模误差时选择的窗口长度为已知周期的最小公倍数120(吴洪宝和吴蕾，2005).首先将残差序列进行MSSA分解，然后将分解的结果通过w-correlation方法分析各重构成分之间的相关性，分析结果如图 4所示，根据1.3节提出的重构标准，发现重构成分的前10阶的相邻两种重构成分相关性大于60%且与其他重构成分相关性之和小于60%，因此将重构成分的前10阶提取出来作为共模误差.

Dong等(2006)使用主成分分析(PCA)对共模误差进行研究，利用共模模式给出了共模误差的判断标准，即某一主分量模式中超过50%的测站的标准化空间响应(每个特征向量除以其绝对值最大元素)大于25%，且该模式的特征值与所有特征值总和的比值超过1%，可认为是共有模式.根据这个准则，我们对利用MSSA提取出的共模误差进行分析，由于本文选用的方法为MSSA，因此，提取出的信号即为所有测站共有的模式.文中对利用MSSA提取出的前14阶分量的标准化空间响应和特征值占比的分析结果进行了统计，如表 2所示.由表 2可以看出，按照Dong等(2006)给出的共模误差的判断准则，前12阶符合要求；将MSSA提取共模误差的结果(表 2中前10阶)与按照Dong等提出的共模误差的判断标准的认定结果(表 2中前12阶)进行对比，前12阶特征值贡献率为84.36%，前10阶特征值贡献率为82.27%，两者的结果非常接近.

3 讨论 3.1 共模误差对时间序列的影响

为了分析共模误差对残差时间序列的影响，本文按照图 1所示的流程，对剔除共模误差前后的残差序列统计标准偏差，并统计剔除前后变化的百分比，统计结果如表 3所示.从表 3的统计结果可以看出，剔除共模误差前后残差的标准差均有明显的减小，减小的幅度均在20%以上，平均减少31.7%，最大的测站减少41.93%，效果比较显著.

通过对剔除共模误差前后的相对坐标时间序列进行分析，发现相对坐标时间序列的标准差最大改善36.89%，平均改善16.5%.由于篇幅限制，本文仅给出了2个测站(ALBH、HARV)剔除共模误差前后的相对坐标的时间序列图和测站振幅的统计结果，如图 5、图 6和表 4所示.从图 5、6和表 4可以看出，相对坐标时间序列在剔除共模误差以后明显收敛，序列的振幅明显变小，说明通过MSSA剔除共模误差后能够明显减小站坐标序列的不确定性，提高坐标的精度.

3.2 共模误差对测站噪声的影响

连续GPS站坐标时间序列中不仅存在白噪声(white noise，WN)，还存在有色噪声，例如闪烁噪声(flicker noise，FN)、随机游走噪声(random walk noise，RWN)和非整数谱指数幂律噪声(power law noise，PL)等(Nikolaidis, 2002).

为了进一步分析共模误差对坐标时间序列精度的影响，本文运用极大似然估计方法对剔除共模误差前后站坐标序列的噪声特性进行了研究，分析了测站的最佳噪声模型，并估计了各噪声分量的量级，如表 5所示.

从表 5可以看出，实验区域GPS坐标时间序列的共模误差同时具有白噪声、闪烁噪声、随机游走噪声以及非整数谱指数幂律噪声的特性.通过MSSA剔除共模误差后，各测站的噪声分量的量级有了明显的降低，白噪声和有色噪声量级平均分别下降30.68%和45.78%，这说明实验区域GPS坐标时间序列的共模误差中的有色噪声占比重较多.通过以上分析可以表明利用MSSA剔除共模误差后能够有效的降低噪声的量级，提高坐标时间序列的可靠性.

通过最小二乘频谱分析法(Williams, 2008; Guo et al., 2012)分析了共模误差的剔除对相对坐标时间序列参数估计的影响，分析结果如表 6所示.从表 6中可以发现剔除共模误差站相对坐标时间序列参数估计的不确定度相比滤波前都有显著下降，速度、周年项振幅、周年项初相、半年项振幅和半年项初相的不确定度平均分别下降了51.48%、46.37%、17.94%、44.8%和21.2%.这说明剔除共模误差能够大幅减小线性项和周期项估计的不确定性；剔除共模误差前后大部分测站相对坐标时间序列的振幅与相位有了微小的变化，有少数测站的变化相对较大，如WRHS台站在去除共模误差前后，其周年和半年振幅出现了明显变化，至于变化较大的原因目前还不得而知，需日后进一步研究.

3.3 共模误差的周期分析

共模误差并不是特指某一种误差，而是某一区域大部分站点都存在的一种或多种误差的集合.为了对各种误差信号进行分离，我们使用MSSA对共模误差序列进行分解并使用快速傅里叶变换(Fast Fourier Transform, FFT)对其进行周期探测，分析结果呈现出比较高的相似性，由于篇幅限制，本文仅以VTIS站和NEAH站为例进行说明，如图 7-10所示.

图 7、图 9和图 8、图 10分别为VTIS测站和NEAH测站共模误差时间序列的MSSA和FFT的分析结果，每个图中包含的20个子图分别为分析结果的前20阶的振幅和周期探测结果，其中图 7和图 9中的子图、图 8和图 10中的子图分别为一一对应的关系.由图 7—10可以看出，共模误差随时间变化呈现出明显的周期特性，说明共模误差是一种与时间相关的公共误差，符合共模误差与时空相关的特性.各测站的共模误差存在的变化周期比较相似，但是不同测站不同周期对应的振幅的变化特性却不尽相同.以VTIS和NEAH两个测站为例来看，共模误差序列中存在诸多周期，且大多数周期为成对出现，其中比较明显的周期有：114天周期(图 9、图 10中f、g图)，此周期对应的振幅表现出不同的变化特征，VTIS站(图 7中f、g图)呈现出随时间逐渐减小的趋势，而NEAH站(图 8中f、g图)则表现出先变大后变小的趋势；85.8天周期(图 9、图 10中h、i图)，示例的两个测站的此周期所对应的振幅(图 7、图 8中h、i图)变化特性均呈现出随时间先增大后减小的趋势；95天周期(图 9、图 10中k、l图)，VTIS站此周期对应的振幅(图 7中k、l图)呈现出随时间逐渐减小的趋势，而NEAH站此周期对应的振幅(图 8中k、l图)的变化比较平稳；69天周期(图 9、图 10中m、n图)，VTIS站此周期对应的振幅(图 7中m、n图)总体上呈现出随时间先增大后逐渐保持平稳的趋势，而NEAH站此周期对应的振幅(图 8中m、n图)的变化比较平稳.以上周期可能是由于参考框架、大气质量负荷、未建模的卫星轨道、地球定向参数(EOP)和交点周期(Wdowinski et al., 1997)等不准确而引入的系统误差，其具体原因目前还没有合理的解释，还需要进一步的研究和探讨.另外，还有些信号的周期性并不明显，如图 9和图 10中的子图a—e，导致此现象的机制目前还不明确，需要进一步进行研究.

4 结论

本文提出了利用MSSA提取共模误差的新思路，并给出了利用MSSA提取GNSS坐标时间序列中共模误差的算法和流程.通过MSSA对实验区域GPS坐标时间序列的处理分析，有效地提取并剔除了GPS坐标时间序列中的共模误差.通过与Dong等(2006)利用PCA算法进行CME研究提出的判断标准进行比较发现，利用MSSA提取的共模误差与Dong等(2006)的标准基本吻合.通过分析剔除共模误差前后时间序列的噪声特性，参数估计的不确定度和时间序列的不确定度，发现剔除共模误差后时间序列的精度有了明显提升，测站噪声和参数估计不确定度有了明显的下降.同时还对共模误差序列进行了周期探测和分析，发现共模误差并非恒定不变，而是随着时间变化而改变；共模误差中有诸多周期信号，产生这些周期的物理机制目前还不得而知，需要进一步探讨.

致谢  感谢匿名审稿专家的意见和建议.感谢SONEL提供GPS坐标序列数据.