CME的可能来源包括区域性的构造运动、残余的电离层和对流层延迟误差以及全球性的参考框架误差等^[1]。去除CME能显著提升GPS坐标时间序列的精度，有利于GPS地球动力学的应用研究。Wdowinski等^[2]提出叠加滤波，能在一定程度上减小噪声水平，但叠加滤波假定CME是均匀分布的，并且不能反映出各站CME的空间响应，存在较大的局限性。Nikolaidis^[3]在叠加滤波的基础上考虑GPS测站单天解的精度，并以此作为权重进行加权叠加滤波，但该方法忽略了GPS站点间的相关性，只适用于区域性且CME近似均匀分布的情况^[4]。Dong等^[5]提出主成分分析法与Karhunen-Loeve联合滤波的方法，取得了较好效果，但仍存在尺度局限性。对于大尺度的GPS观测网络，CME不再均一。为了提取美国达2 000 km尺度的PBO(plate boundary observatory)网络中的CME, SOPAC(Scripps Orbit and Permanent Array Center)将整个PBO网络分为若干个子区，在每个子区中进行叠加滤波。田云峰等^[1]指出该做法存在局限性：子区域的划分缺少依据，子区域交错带的CME差别较大，CME的起源和空间形态不明了，并且提出相关加权叠加滤波。该方法在叠加滤波的基础上考虑了各站残差序列间的相关系数大小，具有良好的空间响应，能显著降低大尺度GPS网络中的噪声水平。

本文在叠加滤波和相关叠加滤波的基础上，提出改进的相关加权叠加滤波方法。该方法顾及GPS站点位置残差序列总体噪声水平，从另一个角度使用相关系数进行加权，取得了比相关加权叠加滤波更好的实验结果。

1 提取共模误差的滤波方法 1.1 相关加权叠加滤波

假设GPS网络中有连续观测m天的n个测站，各测站某一坐标分量的位置残差时间序列组成大小为m×n的时序矩阵$\underset{m\times n}{\mathit{\boldsymbol{v}}}$根据叠加滤波计算CME的公式为：

$ {\mathit{\boldsymbol{\varepsilon }}_{i, j}} = \sum\limits_{k = 1}^n {\frac{{{\mathit{\boldsymbol{v}}_{i, k}}}}{{\sigma _{i, k}^2}}} /\sum\limits_{k = 1}^n {\frac{1}{{\sigma _{i, k}^2}}} $

(1)

式中，ε_{i, j}为第j个站点第i天的CME计算值，ν_{i, k}为第k个站点第i天的残差，σ_{i, k}为第k个站点第i天的单日解均方根误差。

叠加滤波法是一种简单的加权平均值算法，只适用于小尺度、CME分布均匀的GPS网，并且没有空间响应。相关加权叠加滤波引入相关系数作为权重因子，故具有良好的空间响应，克服了叠加滤波的距离局限性。其公式为^[1]：

$ {\mathit{\boldsymbol{\varepsilon }}_{i, j}} = \left( {\sum\limits_{k = 1}^n {\frac{{{\mathit{\boldsymbol{v}}_{i, k}}}}{{\sigma _{i, k}^2}} \times {r_{j, k}}} } \right)/\left( {\sum\limits_{k = 1}^n {\frac{1}{{\sigma _{i, k}^2}} \times {r_{j, k}}} } \right) $

(2)

式中，ε_{i, j}为第j个站点第i天的CME计算值，ν_{i, k}、σ_{i, k}分别为第k个站点第i天的残差和单日解均方根误差，r_{j, k}为第j站点和第k个站点的残差序列间的皮尔逊相关系数：

$ {r_{j, k}} = \frac{{\sum\limits_{i = 1}^m {({j_i} - \overline j )({k_i} - \overline k )} }}{{[\sum\limits_{i = 1}^m {{{({j_i} - \overline j )}^2}{{({k_i} - \overline k )}^2}{]^{1/2}}} }} $

(3)

式中，m为公共历元的个数，j和k为位置残差序列的均值。

相关加权叠加滤波要求参与计算CME的站点距离控制在一定范围内; 参与计算CME的站点与当前站点间的相关系数要高于某个阈值，从而剔除异常站点的影响; 同时，参与计算CME的站点数量也应大于某一阈值^[1]。本文改进的相关加权叠加滤波方法也有相应的要求。

相关加权叠加滤波在计算某一站点的CME时，按照参与计算的所有站点与其相关系数的大小，给予各站不同的权重，故具有良好的空间响应，不受站点间的距离影响。

1.2 改进的相关加权叠加滤波

假设时序矩阵$\underset{m\times n}{\mathit{\boldsymbol{v}}} $对应的单日解均方根误差矩阵为$\underset{m\times n}{\mathit{\boldsymbol{\sigma}}} $，求$\underset{m\times n}{\mathit{\boldsymbol{\sigma}}} $所有元素的平方，再求其倒数得到$\underset{m\times n}{\mathit{\boldsymbol{p}}} = [\underset{m\times 1}{\mathit{\boldsymbol{p}}_1}, \underset{m\times 1}{\mathit{\boldsymbol{p}}_2}, ..., \underset{m\times 1}{\mathit{\boldsymbol{p}}_n}]$计算$\underset{m\times n}{\mathit{\boldsymbol{v}}} $的列向量方差得到矩阵$\underset{1 \times n}{\mathit{\boldsymbol{f}}} $然后再求$\underset{1 \times n}{\mathit{\boldsymbol{f}}} $中各元素的倒数，得到$\underset{1 \times n}{\mathit{\boldsymbol{F}}} = [\underset{1 \times 1}{\mathit{\boldsymbol{F}}_1}, \underset{1 \times 1}{\mathit{\boldsymbol{F}}_2}, ..., \underset{1 \times 1}{\mathit{\boldsymbol{F}}_n}]$；$\underset{1 \times n}{\mathit{\boldsymbol{F}}} $中的元素乘以$\underset{m\times n}{\mathit{\boldsymbol{p}}}$中每列对应的向量，得到$\underset{m \times n}{\mathit{\boldsymbol{q}}} = [\underset{1 \times 1}{\mathit{\boldsymbol{F}}_1} \times \underset{m\times 1}{\mathit{\boldsymbol{p}}_1}, \underset{1 \times 1}{\mathit{\boldsymbol{F}}_2} \times \underset{m\times 1}{\mathit{\boldsymbol{p}}_2}, ..., \underset{1 \times 1}{\mathit{\boldsymbol{F}}_n} \times \underset{m \times 1}{\mathit{\boldsymbol{p}}_n}]$。

改进叠加滤波在叠加滤波的基础上顾及了各GPS站点残差序列本身的噪声水平(方差大小)，根据改进叠加滤波计算CME的公式为：

$ {\mathit{\boldsymbol{\varepsilon }}_{i, j}} = \sum\limits_{k = 1}^n {{\mathit{\boldsymbol{q}}_{i, k}} \times {\mathit{\boldsymbol{\nu }}_{i, k}}} /\sum\limits_{k = 1}^n {{\mathit{\boldsymbol{q}}_{i, k}}} $

(4)

式中，ε_{i, j}为第j个站点第i天的CME计算值，ν_{i, k}为第k个站点第i天的残差。

改进叠加滤波将残差序列总体噪声水平更低(方差更小)的站点赋予更大的权值，这样有2个优点：1)残差序列中的噪声由自身噪声+共性噪声(CME)组成，从大尺度GPS网上看，即使CME分布不均匀，各站CME有一定差异，但GPS网络中的CME大体相似，共性噪声水平相差不会很大，差距主要体现在各GPS站点的自身噪声上。如果某站点噪声水平较低，而各站中的共性噪声水平大致相当，也就意味着该站点自身噪声少，相对的共性噪声就更多，赋予此类站点更多的权重将有利于更彻底地提取CME。2)时间序列中难免会有残留的粗差，在计算CME时，给方差较大(可能残留粗差较多)的站点更小的权值，有利于减小粗差对计算结果的影响。

改进的相关加权叠加滤波在计算GPS观测网络内某一站点(代号为i)的共模误差ε_i时，首先和其余站点(代号为1, 2, …, n-1)进行两两之间的改进叠加滤波计算，解得n-1个CME向量$\left({{\mathit{\boldsymbol{\varepsilon }}_{\left({i, 1} \right)}}, {\mathit{\boldsymbol{\varepsilon }}_{\left({i, 2} \right)}}, {\mathit{\boldsymbol{\varepsilon }}_{\left({i, n - 1} \right)}}} \right)$，再根据它与这些站点的相关系数$\left({{r_{i, 1}}, {r_{i, 2}}, {r_{i, n - 1}}} \right)$将这n-1个不同的CME向量进行加权算术平均值计算，最终得到该站的共模误差ε_a。

用改进叠加滤波求2个GPS站点(代号为i和j)之间的CME时，式(4)可简化为：

$ {{\mathit{\boldsymbol{\varepsilon }}_{k, i}}} = {{\mathit{\boldsymbol{\varepsilon }}_{k, j}}} = \frac{{{\mathit{\boldsymbol{q}}_{k, i}} \times {\mathit{\boldsymbol{\nu }}_{k, i}} + {\mathit{\boldsymbol{q}}_{k, j}} \times {\mathit{\boldsymbol{\nu }}_{k, j}}}}{{{\mathit{\boldsymbol{q}}_{k, j}} + {\mathit{\boldsymbol{q}}_{k, j}}}} $

(5)

式中，ε_{k, i}、ε_{k, j}为第i个和第j个站点第k天的CME计算值，ν_{k, i}、ν_{k, j}为第i个和第j个站点第k天的残差。用ε_{(i, j)}来表示对站点i和j使用改进叠加滤波求得的CME向量，则有${\mathit{\boldsymbol{\varepsilon }}_{\left({i, j} \right)}} = \left[{{\mathit{\boldsymbol{\varepsilon }}_{1, i}}{\rm{; }}{\mathit{\boldsymbol{\varepsilon }}_{2, i}}; ...; {\mathit{\boldsymbol{\varepsilon }}_{m, i}}} \right] = \left[{{\mathit{\boldsymbol{\varepsilon }}_{1, j}}; {\mathit{\boldsymbol{\varepsilon }}_{2, j}}; ...; {\mathit{\boldsymbol{\varepsilon }}_{m, j}}} \right]$。

改进的相关加权叠加滤波方法求CME的最终表达式为：

$ {\mathit{\boldsymbol{\varepsilon }}_i} = \sum\limits_{j = 1}^{n - 1} {{\mathit{\boldsymbol{\varepsilon }}_{\left( {i, j} \right)}} \times {r_{i, j}}} /\sum\limits_{j = 1}^{n - 1} {{r_{i, j}}} $

(6)

式中，ε_i为第个i站点最终解得的CME向量，ε_{(i, j)}为对站点i和j使用改进叠加滤波得到的CME向量, r_{i, j}为站点i和j的残差序列间的相关系数。

2 算例分析 2.1 实验数据

本次实验选取中国陆态网8个GPS站点(DXIN、BJSH、JIXN、TAIN、YANC、XIAA、LUZH、XIAM)的坐标序列数据。对GPS坐标残差序列的CME进行滤波前还需进行数据预处理，本文采用四分位间距法去除粗差^[6]，使用GMIS软件作数据插值^[7]，然后进行阶跃修复、去线性趋势和去均值处理。规定参与计算CME的站点之间的相关系数不低于0.3，参与计算CME的GPS总站点数量不低于4个。以N方向为例，经预处理得到的最终残差时间序列如图 1。

2.2 实验结果

针对以上GPS站点，分别采用相关加权叠加滤波和本文改进的相关加权叠加滤波计算各站的CME，然后通过均方根(RMS)改进率和相关系数总改正量来评估2种方法的性能。

滤除2种方法提取到的CME后的N方向结果见图 2。可以直观地看出，2种方法滤波后对应GPS站点的残差序列形状基本一致，各站点残差序列间的相关性基本去除，噪声水平大大降低，但相对而言，本文方法的滤波结果更优。E、U方向结果与N方向类似。

表 1给出滤波后GPS站点平均RMS值的变化情况，其中的数值都是平均值，改进率即滤波后RMS的降低百分比。相比相关加权叠加滤波，本文方法的GPS站残差序列的平均RMS改进率要高20%左右，去噪性能更强。

在处理方法合理的情况下，去除CME后的GPS站点残差序列间的相关性越低，说明CME去除越彻底。表 2给出相关系数总量R的变化情况，R的计算式如下：

$ R = \sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^n {{r_{i, j}}} } $

(7)

式中，n为GPS站点残差序列编号，r_{i, j}为站i与站j的残差序列间的相关系数。可以看出，本文方法的相关系数总改进量略优于相关加权叠加滤波，说明CME去除更彻底。

3 结语

本文方法在计算GPS观测网络内某站点的CME时，假设该站点与其余站点之间都存在均匀的CME，利用改进叠加滤波逐一提取它与其余站点间的初步CME。而实际上该站点与其他站点间的空间相关性大小不一，网络中的CME分布也不均匀，该站点与相关性更大的站点解得的CME应具有更大的参考价值，与相关性过低的站点解得的CME则不应纳入计算之中。故引入相关系数作为权重因子给这些初步CME进行加权平均值计算，最终得到观测网络中该站点的CME。改进的相关加权叠加滤波具有良好的空间响应，能准确地提取GPS观测网络中的CME。通过分析我国陆态网的8个GPS站点的坐标数据得到，相对于相关加权叠加滤波，本文改进的相关加权叠加滤波拥有更好的性能，可进一步提高RMS改进率和信噪比，有利于GPS地球动力学应用研究，也为CME的起源分析提供了参考。