云南GNSS 时间序列共模误差提取与站点形变分析

徐良叶，邵德盛，吴学群，牛 甜，韩 逍

（1. 云南省地震局，云南 昆明 650200；2. 昆明理工大学 国土资源工程学院，云南 昆明 650093）

0 引言

连续GNSS 时间坐标序列中存在一种空间相关的误差称之为共模误差（commno Mode Error，CME），共模误差被普遍认为是GNSS 坐标时间序列的主要来误差源之一，一些弱构造或瞬态构造信号会被共模误差所掩盖，它会影响区域网中站点形变特征获取的准确性[1-3]。关于共模误差的起源，还有待进一步研究，已有研究表明它可能是由以下几个方面引起的[4-5]：软件、算法及数据处理策略等不完善而引起的系统误差；GNSS 卫星轨道误差，海潮改正残差，地球旋转参数误差等模型误差；环境负载，小规模（区域性） 的地壳变形影响。

坐标时间序列中的共模误差可以通过PCA/KLE 相结合的滤波方法有效地提取出来[6-9]，为了更准确的提取共模误差，需要剔除本地效应强的站点，而被剔除的站点时间序列里仍然包含有用的构造运动信息，相关方面的研究较少，所以本文基于云南2015—2019 年GNSS 坐标时间序列，首先根据2 倍标准差的方法剔除本地效应强的站点，利用PCA/KLE 法来提取云南GNSS 时间序列的共模误差；接着对本地效应强的站点的时间序列进行去噪处理，结合地震信息，分析震前站点的形变信息。

1 PCA 和KLE 方法

主成份分析（Principal Component Analysis）PCA 也叫经验正交函数分析（Empirical Ortho—gonal function） 即EOF 或特征向量分析。它把区域网站点时间序列分解成时间域的主分量和空间域的特征向量，目的在于通过利用降维的思想把多指标降为少数几个综合指标[6]。

如果把所有站点的时间序列排列起来，1个站点时间序列为一列，可表为m×n的矩阵（通常m＞n）：

它是一个n×n的矩阵，这时X的方差矩阵可表为：

矩阵VT是由特征向量构成的n×n的正交矩阵，式中Λ 矩阵是由k个特征值构成的对角矩阵，这时V为X的正交基底，X可以由此正交基底表示为：

式中的展开系数ak（ti）的序列即为k阶主分量，它可以由下式求出：

式中的ak（t）表示第k个主成分，vk（x）为对应主成分的响应特征矩阵，ak（t）代表时间特征，vk（x）代表空间响应。将B 矩阵的特征值λ 按从大到小的排列，特征值越大代表对残余时间序列方差的贡献越大[10]。可以用特征值累积贡献率mk来代表主成分对矩阵的贡献率，mk可由下式表示：

相关矩阵C可由协方差矩阵B标准得到，即dij=bij/（σiσj），其中σ 的定义为：

这个过程就称为Karhunen-Loeve 展开（KLE）[11]。

PCA 和KLE 的唯一区别就是用来计算正交向量的矩阵不同。PCA 算法中使用的是协方差矩阵，KLE 算法中使用的是标准化后的相关系数矩阵。由于PCA 方法能得到精度较高的主成分值和空间响应，但如果本地噪声较大的话，则无法准确地分辨出CME，而KLE 方法正好弥补了这个缺点，它可以抑制本地噪声的影响，并提取出CME，它的缺点是得不到准确的空间相应[6，8]。

2 数据处理

本文数据来自于中国地震局第一监测中心解算并免费提供的GNSS 数据产品（http://www.eqdsc.com/)。选取云南及周边2015—2019 年数据比较完整的33 个GNSS 基准站的观测资料做分析（数据缺失率小于10%）。站点分布如图1所示。

图1 云南及周边站点分布图Fig.1 Distribution of stations in Yunnan and surrounding areas

2.1 获取GNSS 时间序列残差

GNSS 台站坐标每日观测时间序列拟合函数为[12]：

式中ti表示历元（单位：年），a为初始位置，b为速度，c、d为半周年系数，e、f为年数为Tj时刻的阶跃式偏移量，vi为残差序列。

使用QOCA 软件去除原始坐标时间序列中的阶跃、速度、周年、半周年项，得到坐标残差时间序列vi，限于篇幅本文以弥勒站（YNML）为例，如图2 所示。

图2 YNML 站原始坐标时间序列（左图） 及残差坐标时间序列（右图）Fig.2 Original coordinate time series of YNML Station（left） and residual coordinate time series（right）

2.2 异常站点剔除与共模误差提取

如果没有剔除本地效应较强的站点，会导致PCA 结果有一定程度的偏差，文献[6]研究表明，本地效应较强的站点具有显著的空间响应，且分布不均匀，此时PCA结果的PC1（第一主成份） 对应的空间响应与KLE 结果与大多数站相比有明显的差异。所以本文根据PCA结果与KLE 结果的PC1 对应的空间响应差值的2 倍标准差来剔除异常站点，如下式所示：

根据2.1 节获取的残差时间序列作为PCA/KLE 的输入数据，由（13）式来剔除异常站点，再重复进行PCA/KLE 计算来提取PC1，直到没有异常站点。本文循环了4 次，共剔除了9 个站点，剩余24 个站点，将剩余站点通过PCA方法提取共模误差。

PCA提取前3 个主成分在N、E、U 方向的贡献率分别为67.25%、4.07%、3.81%；43.04%、8.55%、5.90%；41.69%、6.74%、5.59%，在N、E、U 方向的累计贡献率分别为75.13%、57.49%、53.95%，可见主要的方差集中在PC1上。由PCA 得到的前3 个主成分空间响应如图3 所示，PC1 的空间响应分布均匀且都为正，在N、E、U 方向的均值为0.8867、0.8798、0.8622，PC2、PC3 的空间响应分布不均匀，有正有负，由于共模误差表示空间相关的误差，PC1 的空间响应一致性较好，所以本文选取PC1 作为共模误差。

图3 PCA 解算的前3 个主成分对应的空间响应Fig.3 Spatial responses of the first three principal components solved by PCA

PCA 提取的共模误差如图4 所示（以YNML 站为例），E、N 方向的共模误差大多数在-8～8mm 之间，U 方向的共模误差大多数在-14～14mm 之间，说明云南GNSS 坐标时间序列中存在比较明显的共模误差。

图4 PCA 提取的共模误差Fig.4 Common mode error extracted by PCA

对比PCA 去除共模误差前与去除共模误差后的坐标时间序列（图5） 可以看出，去除共模误差后，N、E、U 方向的坐标时间序列波动幅度减小了，并趋于平稳。

图5 YNML 站PCA 滤波前和滤波后的坐标时间序列Fig.5 Coordinate time series of YNML station before and after PCA filtering

如图6（蓝色小三角形表示滤波前的RMS，红色小三角形表示滤波后的RMS） 所示，PCA滤波前的坐标时间序列在在N、E、U 方向的均方差（RMS） 为2.27mm、2.25mm、6.14mm，滤波后为1.07mm、1.39mm、3.98mm，说明该方法很好地剔除了共模误差，提高了坐标精度。

图6 PCA 滤波前和滤波后的坐标时间序列均方差Fig.6 Mean square deviation of coordinate time series before and after PCA filtering

2.3 异常站点共模误差提取及形变分析

使用KLE 滤波能够有效地减小本地噪声的影响，并提取出CME，所以本文基于KLE 方法提取本地效应较强的站点（2.2 节中被剔除的异常站点） 坐标时间序列中的CME，如图7 所示，以新平站（YNXP）为例，图8 为去除共模误差后的坐标时间序列，N、E、U 方向的坐标时间序列比之前稳定。KLE 滤波前YNXP 站坐标时间序列在N、E、U 方向的均方差（RMS）为3.22mm、4.12mm、6.31mm，滤波后为2.10mm、2.76mm、4.73mm，说明剔除CME 后站点的坐标精度和可靠性得到提高。

图7 YNXP 站原始坐标时间序列（左图） 与残差坐标时间序列（右图）Fig.7 Original coordinate time series（left） and residual coordinate time series of YNXP Station（right）

研究表明[13]，去除共模误差后的坐标时间序列依旧含有噪声，包括白噪声、闪烁噪声、随机漫步噪声，其中的两种或三种混合噪声等。而自适应卡尔曼滤波（Kalman） 的前提要求是数据中的噪声为高斯分布的白噪声[14]，小波变换被誉为“数学显微镜”，它能从时频域的局部信号中有效地提取信息，具有较好的去噪特点[15]。

所以本文利用文献[16]中的小波变换和自适应Kalman 滤波算法，将2.3 节中去除共模误差后的坐标时间序列进行小波变换处理，再将处理后的结果作为自适应Kalman 滤波处理的输入数据，其目的是更科学、有效地降低坐标时间序列中的噪声，最终处理后的结果（图9）与没有进行去噪处理的坐标时间序列（图8 中右图）相比，可以看出，去噪后坐标时间序列曲线更光滑，更能直观地凸显地壳的形变信息。

图8 YNXP 站KLE 滤波前坐标时间序列（左图） 与滤波后坐标时间序列（右图）Fig.8 Coordinate time series before YNXP Station KLE filtering（left） and after filtering coordinate time series（right）

图9 YNXP 站去噪后的坐标时间序列Fig.9 Coordinate time series of YNXP station after noise remova

为了验证该方法的合理性，本文对滤除的噪声（即去噪之前的坐标时间序列与去噪后的时间序列的差值） 进行统计分析，如图10 所示，滤除的噪声符合正态分布，说明本文去除噪声方法的可靠性。

图10 去除的噪声统计直方图Fig.10 Histogram of removed noise statistics

本文结合2018 年8 月13 日通海5.0 级地震，分析地震前YNXP 站点的异常形变信息，如图11 所示，图中的红线表示坐标时间序列的均值加减2 倍标准差所对应的限差值，明显可以看出从2018 年5 月份开始，YNXP 站点在N、E 方向的运动趋势严重偏离了均值线，并超过了2 倍的标准差，N 方向的位移在2018 年5月30 日达到最大值-4.82mm，E 方向的位移在2018 年6 月17 日达到最大值-5.54mm；同样从图7 中也可以看出YNXP 站点在N、E 方向原始坐标时间序列和残差坐标时间序列均出现了凹槽，地震后站点恢复到原有的运动趋势。相关工作人员在2018 年5 月对站点进行异常核实，排除了数据质量问题和周边环境影响，认为YNXP 站点的变化为真实的构造运动引起。说明此方法能凸显地壳形变和构造运动信息，为提取GNSS 坐标时间序列中的信息提供了参考。

图11 YNXP 站去噪后的坐标时间序列Fig.11 Coordinate time series of YNXP station afternoise removal

3 结论

本文采用PCA 和KLE 方法提取云南GNSS坐标时间序列的共模误差，然后对坐标时间序列进行去噪处理，结合2018 年通海5.0 级地震，分析地震前站点异常形变信息。

结果表明，本文提出的定量去除本地效应较强的站点的方法能更准确地从GNSS 时间序列中去除共模误差，有效地提高了坐标精度；对减小本地效应误差后的时间序列进行去噪处理后，坐标时间序列曲线更加光滑，通海5.0级地震前YNXP 站点在N、E 方向的运动趋势严重偏离了均值线，加速向西南方向运动，并超过了2 倍的标准差，说明此方法能凸显地壳形变和构造运动信息，为提取GNSS 坐标时间序列中的信息提供了参考和新的思路。