文章快速检索     高级检索
  大地测量与地球动力学  2019, Vol. 39 Issue (9): 924-927  DOI: 10.14075/j.jgg.2019.09.010

引用本文  

谢玖刚, 匡翠林, 周要宗. 一种改进的GPS坐标序列共模误差滤波方法[J]. 大地测量与地球动力学, 2019, 39(9): 924-927.
XIE Jiugang, KUANG Cuilin, ZHOU Yaozong. An Improved Method for Filtering Common-Mode Error of GPS Coordinate Time Series[J]. Journal of Geodesy and Geodynamics, 2019, 39(9): 924-927.

项目来源

国家自然科学基金(41774040)。

Foundation support

National Natural Science Foundation of China, No. 41774040.

第一作者简介

谢玖刚,硕士生,主要从事GPS时间序列分析研究,E-mail:xjgcsu@163.com

About the first author

XIE Jiugang, postgraduate, majors in GPS time series analysis, E-mail:xjgcsu@163.com.

文章历史

收稿日期:2018-09-01
一种改进的GPS坐标序列共模误差滤波方法
谢玖刚1     匡翠林1     周要宗1     
1. 中南大学地球科学与信息物理学院,长沙市麓山南路932号,410083
摘要:提出一种改进的相关加权叠加滤波方法,顾及GPS站点位置残差序列的总体噪声水平,有利于抑制站点自身噪声和残余粗差给共模误差(common-mode error,CME)计算带来的不利影响。该方法采用相关系数作为权重因子,拥有良好的空间响应。选取我国陆态网8个GPS站点的数据对该方法的适用性进行分析,结果表明,相对于相关加权叠加滤波,其拥有更好的性能,均方根(RMS)改进率提高20%左右,能更加准确地提取CME。
关键词GPS坐标时间序列共模误差相关加权叠加滤波

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天的残差和单日解均方根误差,rj, 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为公共历元的个数,jk为位置残差序列的均值。

相关加权叠加滤波要求参与计算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站点(代号为ij)之间的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)来表示对站点ij使用改进叠加滤波求得的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)为对站点ij使用改进叠加滤波得到的CME向量, ri, j为站点ij的残差序列间的相关系数。

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

图 1 GPS站点N方向坐标残差序列 Fig. 1 GPS coordinate residual series of N direction
2.2 实验结果

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

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

图 2 N方向残差序列滤波后对比 Fig. 2 The comparison of N direction residual series after filtering

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

表 1 滤波后GPS站点残差序列平均RMS变化 Tab. 1 The average RMS of GPS stations residual time series after filtering

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

表 2 滤波后GPS站残差序列相关系数总量变化 Tab. 2 Residual series correlation coefficient total change of GPS stations after filtering
$ R = \sum\limits_{i = 1}^{n - 1} {\sum\limits_{j = i + 1}^n {{r_{i, j}}} } $ (7)

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

3 结语

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

参考文献
[1]
田云锋, 沈正康. GPS观测网络中共模分量的相关加权叠加滤波[J]. 地震学报, 2011, 33(2): 198-208 (Tian Yunfeng, Shen Zhengkang. Correlation Weighted Stacking Filtering of Common-Mode Component in GPS Observation Networks[J]. Acta Seismologica Sinica, 2011, 33(2): 198-208 DOI:10.3969/j.issn.0253-3782.2011.02.007) (0)
[2]
Wdowinski S, Bock Y, Zhang J, et al. Southern California Permanent GPS Geodetic Array: Spatial Filtering of Daily Positions for Estimating Coseismic and Postseismic Displacements Induced by the 1992 Landers Earthquake[J]. Journal of Geophysical Research: Solid Earth, 1997, 102(B8): 18 057-18 070 DOI:10.1029/97JB01378 (0)
[3]
Nikolaidis R. Observation of Geodetic and Seismic Deformation with the Global Positioning System[D]. San Diego: University of California, 2004 (0)
[4]
贺小星, 姜卫平, 周晓慧, 等. GPS坐标时间序列广义共模误差分离方法[J]. 测绘科学, 2018, 43(10): 7-15 (He Xiaoxing, Jiang Weiping, Zhou Xiaohui, et al. A General Common Mode Error Filtering Method for GPS Time Series[J]. Science of Surveying and Mapping, 2018, 43(10): 7-15) (0)
[5]
Dong D, Fang P, Bock Y, et al. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve Expansion Approaches for Regional GPS Network Analysis[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B3) (0)
[6]
黄立人. GPS基准站坐标分量时间序列的噪声特性分析[J]. 大地测量与地球动力学, 2006, 26(2): 31-33 (Huang Liren. Noise Properties in Time Series of Coordinate Component at GPS Fiducial Stations[J]. Journal of Geodesy and Geodynamics, 2006, 26(2): 31-33) (0)
[7]
Liu N, Dai W J, Santerre R, et al. A MATLAB-Based Kriged Kalman Filter Software for Interpolating Missing Data in GNSS Coordinate Time Series[J]. GPS Solutions, 2018, 22(1): 25 DOI:10.1007/s10291-017-0689-3 (0)
An Improved Method for Filtering Common-Mode Error of GPS Coordinate Time Series
XIE Jiugang1     KUANG Cuilin1     ZHOU Yaozong1     
1. School of Geosciences and Info-Physics, Central South University, 932 South-Lushan Road, Changsha 410083, China
Abstract: This paper proposes an improved correlation weighted stacking filtering method. The method takes into account the overall noise level of the GPS station position residual series, which helps to suppress the adverse effects of the site's own noise and residual gross error on the CME calculation, and uses the correlation coefficient as the weighting factor to generate a good spatial response. Furthermore, the feasibility of the method is analyzed using the GPS time series data of 8 CMONOC stations. The analysis results show that the RMS improvement rate of the method is about 20% relative to correlation weighted stacking filtering, has better performance, and can extract CME more accurately.
Key words: GPS coordinate time series; common-mode error; correlation weighted stacking filtering