随着全球导航卫星系统(global navigation satellite system, GNSS)连续运行观测站时序资料的不断积累，测站坐标时间序列的误差特性以及其组合噪声的定量估计引起了人们的不断重视。目前大量研究已经表明，测站空间相关的共模误差是坐标时间序列的主要误差源之一。自Bock等^[1]提出共模误差的概念，并基于其在研究区域均匀分布的假设，利用区域叠加滤波(即堆栈法)定量估计其大小，国内外学者开始进入深入研究^[1-17]，Marquez-Azua等^[2]考虑了共模误差空间分布的非均匀性, 在计算时加入了距离和时间序列长度作为权重因子；蒋志浩等^[3]统计了利用主成分分析法剔除共模误差前后坐标序列重复性以及周期项的变化，但是缺少估计速度的对比；高涵^[4]利用美国中南部seu1区域13个全球定位系统(global positioning system, GPS)参考站5年的坐标时间序列，仅分析了主成分分析法提取的共模误差对测站水平速度的影响，并认为在提取地壳微动态形变信息时不容忽视；虽然吴伟伟^[5]曾分析滤除共模误差前后坐标时间的速度估值及精度，但使用的是堆栈法滤除共模误差，且研究的测站坐标序列长度多为3年左右。

为了探究共模误差对估计中国陆态网(Crustal Movement Observation Network of China, CMONOC)连续运行观测站三维速度场的影响，本文基于精密单点定位(precise point positioning, PPP)技术解算出的62个GNSS连续运行观测站近5年坐标时间序列，利用谢树明等改进的，适用于大空间尺度GPS网的相关系数加权叠加滤波法提取共模误差^[9]，并使用时间序列分析软件(cointegration analysis of time series, CATS)估计测站各分量剔除共模误差前后坐标时间序列的速度及精度，对比对应估计值的变化，发现由剔除共模误差后的坐标时间序列而得的速度与陆态网发布速度的偏差均值更小，估计的精度更高。

1 测站数据处理及解算精度评定 1.1 测站数据处理

本文选取了中国境内分布均匀的52个陆态网连续运行观测站及周边10个国际GNSS服务(International GNSS Service, IGS)参考站(图 1)作为研究对象，时间跨度为2011-01-01~2015-12-31。

利用自己研发的GNSS数据分析处理工作站软件(PASSION)中的精密单点定位模块，对每个测站进行坐标解算，获取单日解，并将单天收敛序列的标准差作为单日解内符合精度。利用IGS提供的钟差文件以及事后精密星历来改正卫星的钟误差和轨道误差，测站坐标、接收机钟差和对流层延迟模型化后的残余误差均通过卡尔曼滤波进行估计。表 1给出了本文使用的定位程序中的部分处理策略和误差改正模型。

PPP解算坐标的框架与精密星历的框架一致，本次数据处理根据不同时段选择IGS05/IGS08/IGb08星历，解算结果是对应ITRF05/ITRF08框架下的坐标，然后利用国际地球参考框架(international terrestrial reference frame, ITRF)官网发布的转换参数，将所有坐标统一转换到ITRF08框架下，并分别选取各站地心坐标序列中2011年第一个有PPP解算坐标的历元的三维坐标作为站心，将对应的地心坐标转换成站心坐标序列^[22]，结果标记为S0。

最后对任一测站，处理其各年原始的各站心坐标分量，将满足式(1)的历元点剔除：

$ \left| {{X_{ij}} - {{\bar X}_i}} \right| > 3\sqrt {\frac{1}{{{n_i} - 1}}\sum\limits_{j = 1}^{{n_i}} {{{\left( {{X_{ij}} - {{\bar X}_i}} \right)}^2}} } $

(1)

式中，X_ij表示某一测站在任一分量第i年，第j个历元的坐标值；X_i表示该测站在对应分量上，第i年所有历元坐标的均值；n_i表示该分量上，测站第i年的观测历元总数。将剔除粗差后的各站站心坐标序列标记为S1。

1.2 坐标精度评定

在http://sopac-ftp.ucsd.edu/pub/timeseries/measures/中下载本文选取的10个IGS站在ITRF08框架下的地心坐标序列，它们分别是日本鹿儿岛站(AIRA)，中国房山站(BJFS)、长春站(CHAN)、拉萨站(LHAZ)、佘山站(SHAO)、乌鲁木齐站(URUM)、武汉站(WUHN)，韩国水原站(SUWN)，蒙古乌兰巴托站(ULAB)，并将之转换成站心坐标序列(转换方法与§2.1中对应测站相同)，然后与本文解算得到的10个IGS站的S1序列对比，将测站各分量上两个时间序列残差的均方根误差(root-mean-square error, RMSE)，作为评定PPP解算质量的依据。

$ {\rm{RMSE}} = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{X_{{\rm{sopac}}, i}} - {X_{{\rm{ppp}}, i}}} \right)}^2}} } $

(2)

式中，n表示测站在指定分量下的坐标总数; X_{sopac, i}表示斯克里普斯轨道和永久阵列中心(scripps orbit and permanent array center, SOPAC)发布的相同历元的站心坐标值; X_{ppp, i}表示PPP解算的第i个历元站心坐标值。各站东(E)、北(N)、高(U)方向的统计结果如表 2所示。

从表 2可以看出，对于所选的10个IGS站，本文利用PPP解算出的各坐标分量与SOPAC发布坐标差值的平均RMSE在15.5 mm内，其中各站N分量的RMSE最小。总体来说，PPP解算坐标与SOPAC发布的坐标差异较小，解算结果质量较好，可用于后续分析。

2 相关加权叠加滤波

在数据处理中，首先对所有测站各分量序列S1进行拟合，获取去除趋势项、年周期以及半年周期项的残差序列，然后在谢树明等^[9]提出的公式基础上，设置相关系数以及参与计算的测站数阈值，进行空间滤波，并比较剔除共模误差前后所有测站各分量残差序列的均方根误差的变化。

2.1 获取拟合残差

根据地震记录文件、测站日志文件以及各分量坐标序列图，整理S1中含有的已知阶跃以及肉眼判断明显的阶跃历元，并采用式(3)进行拟合，得到去除趋势项、年周期以及半年周期项的残差序列，标记为R：

$ \begin{array}{*{20}{c}} {y\left( {{t_i}} \right) = a + b{t_i} + c\sin \left( {2\pi {t_i}} \right) + d\cos \left( {2\pi {t_i}} \right) + }\\ {e\sin \left( {4\pi {t_i}} \right) + f\cos \left( {4\pi {t_i}} \right) + }\\ {\sum\limits_{i = 1}^{{n_k}} {{g_i}} H\left( {{t_i} - {T_{kj}}} \right) + {v_i}} \end{array} $

(3)

式中，y(t_i)表示t_i时间测站对应分量下，经预处理而得的站心坐标；a表示参考时刻的位置；b表示线性速率；c、d为年周期项系数；e、f表示半年周期项系数；t_i表示特定分量第i个历元的时间值；v_i表示残差。

2.2 计算相关系数

利用相关系数表征相同分量下测站间残差时间序列的线性相关程度，其计算公式为：

$ {r_{xy}} = \frac{{\sum\limits_{k = 1}^n {\left( {{x_k} - \bar X} \right)} \left( {{y_k} - \bar y} \right)}}{{{{\left[{\sum\limits_{k = 1}^n {{{\left( {{x_k}-\bar X} \right)}^2}} \sum\limits_{k = 1}^n {{{\left( {{y_k}-\bar y} \right)}^2}} } \right]}^{1/2}}}} $

(4)

式中，x和y分别表示相同分量下仅含所有公共历元的两个测站残差序列；x_k和y_k分别表示两站第k个历元的残差值；X、y分别表示两站所有公共历元残差的均值；n表示公共历元的个数。最后统计各分量测站间平均相关系数如表 3所示。

从表 3可以看出，所有测站E、N、U分量相关系数均值分别为0.32，0.17，0.26，N分量的均值较小，且各个分量中的均存在站间负相关的现象。

2.3 计算共模误差

基于测站各分量S1的残差序列R，利用相关系数和坐标单日解精度联合定权，即$\frac{{{r_{i, k}}}}{{\delta _{k, j}^2}}$，采用式(5)计算测站在特定分量下测站i第j天的共模误差ε_{i, j}：

$ {\varepsilon _{i, j}} = \frac{{\sum\limits_{k = 1}^{{N_{i, j}}} {\frac{{{v_{k, j}}}}{{\delta _{k, j}^2}}} \times {r_{i, k}}}}{{\sum\limits_{k = 1}^{{N_{i, j}}} {\frac{1}{{\delta _{k, j}^2}}} \times \left| {{r_{i, k}}} \right|}} $

(5)

式中，N_{i, j}表示在获取测站i第j天的共模误差时参与计算的测站数，设置其阈值为3；v_{k, j}，δ_{k, j}分别表示特定分量下测站k第j天的残差值和坐标精度；r_{i, k}表示测站i，k的在特定分量下的相关系数，其阈值为0。最后所有测站各分量滤波前后残差序列的RMS均值如表 4所示。

从表 4可以看出，相对于滤波前所有测站各分量原始的残差序列计算而得到的平均RMS，基于剔除共模误差后的残差序列计算得到的平均RMS有明显的下降，在E、N、U分量分别减小26.11%，20.49%，21.29%，这表明本文采用的相关加权叠加滤波法能够有效地提取出测站i第j天的共模误差，可用于后续的研究。

3 速度估计

由于GPS的坐标序列中含有与时空相关的有色噪声(colored noise, CN)已经成为共识，本文根据Nikolaidi^[23]、黄立人^[24]、吴伟伟等^[5]的研究成果，直接假设剔除共模误差前后测站各分量坐标时间序列的噪声组合均是白噪声(white noise, WH)与幂律噪声(power-law noise, PL)。将在S1基础上剔除共模误差后的序列标记为S2，并分别将S1、S2作为Williams^[25]编写的CATS软件的输入数据，基于趋势项、年周期、半年周期的函数模型，估计其速度及精度；以与本文处理数据的时间段最接近为原则，从陆态网官网收集其采用GIPSY、或GAMIT处理而得的相关测站站心速度，以便对比分析。整理由所有测站各分量序列S1、S2估计而得的幂律噪声和白噪声均值如表 5所示。

从表 5可以看出，相对于由所有测站各分量序列S1而得的幂率噪声和白噪声的平均估值，各分量序列S2估计而得的幂率噪声和白噪声均有较大幅度的减小。对于幂律噪声而言：E、N、U分量分别减小了30.44%，22.15%，28.07%；对于白噪声而言：E、N、U分量分别减小了10.46%，9.05%，16.26%。进一步验证出本文所采用的相关加权叠加滤波能够有效地估计出残差序列中含有的部分白噪声和幂律噪声。并且不管是在原始序列S1还是在剔除共模误差后的序列S2中，白噪声都不是主要成分，这与Williams^[26]、黄立人等^[24]的研究结果一致。

整理所有测站各分量序列S1与S2估计速度及其精度间偏差的均值发现，E、N、U分量分别为0.48±1.15 mm/a，0.06±0.07 mm/a，0.42±0.69 mm/a，这说明在剔除共模误差前后，测站估计速度具有较好的一致性，这与王伟^[27]的研究结果相同。同时整理由所有测站各分量序列S1、S2估计而得到的速度与陆态网发布速度偏差均值，以及平均的速度估计精度如表 6所示。

从表 6可以看出，对于E、N、U分量，基于剔除共模误差后的坐标序列整理而得的平均的速度偏差分别为2.12 mm/a，0.88 mm/a，3.27 mm/a，均小于由原始坐标序列估计而得的平均速度偏差，分别减小了7.32%，1.87%，7.84%；同样地，所有测站各分量估计的速度精度均值也表现出了相同情况，分别提高了33.43%，21.11%，59.26%。

4 结束语

本文利用独立研发的精密单点定位程序PASSION，解算出中国大陆及周边62个GNSS连续运行观测站站近5年的坐标时间序列。研究发现：相对于所有测站各分量原始的残差序列计算而得的平均RMS，剔除共模误差后残差序列的平均RMS至少减小了20%，即本文采用的相关系数加权的叠加滤波法，能够很好地估计出陆态网测站坐标序列中的存在的共模误差；利用CATS软件分别估计测站各分量剔除共模误差前后的坐标时间序列的组合噪声以及速度的结果显示，相对于原始坐标时间序列，不管是测站各分量组合噪声的平均估值，还是各分量估计速度与陆态网发布速度偏差的均值，由剔除共模误差的坐标时间序列整理而得值明显更小，同时各分量估计速度精度的均值至少提升了21%，这说明虽然由剔除共模误差前后的坐标时间序列估计而得的速度具有较好的一致性，但若不剔除原始坐标序列中的共模误差，陆态网测站的运动速度及其精度可能会被不恰当地估计^[28-32]，为了测站的速度估值更加准确可靠，建议在估计区域的多个测站速度场之前，先剔除原始坐标序列中的共模误差。