CN105572703A - 一种gps时间序列广义共模误差提取方法 - Google Patents

Abstract

本发明公开了一种GPS时间序列广义共模误差提取方法，本发明对共模误差的空间响应进行了分析，建立了大区域GPS网共模误差的提取机制，引入了相关系数ρ、距离及经纬度、本地效应、负荷效应、主分量贡献率及其空间响应等地理因素作为评价因子，通过聚类分析为大区域、大尺度下GPS网站点间共模误差，提供了可行的估计方法，获得更好的滤波效果。同时本发明考虑了传统GPS时间序列模型的局限性，减弱了时间序列周、半年年项拟合过程中引入的模型误差，且顾及了共模误差的周期性，即对残差序列进行共模误差剔除之前，保留了原始坐标序列的周年、半年项，有助于真实反映共模误差的空间变化及周期性变化，为进一步提高GPS坐标序列模型的精度提供依据。

Description

一种GPS时间序列广义共模误差提取方法

技术领域

本发明属于连续运行全球定位导航系统技术领域，涉及一种GPS时间序列广义共模误差提取方法。

背景技术

随着卫星定位系统测量定位精度的提高，全球分布的连续运行GPS(GlobalPositioningSystem)跟踪站积累了数十年的观测资料，得到了许多有价值的认识。

GPS坐标时间序列中不仅包含着构造运动信号，也包含着非构造运动信号、季节性信号等噪声，夹杂的噪声影响GPS解的可靠性，对一些地球物理现象甚至可能做出错误的解释。如何高效、快速的剔除这些噪声，分离非构造运动信号的影响，是GPS坐标时间序列中的关键问题之一。

在不同站点坐标时间序列(GPS网)中存在一种空间相关的误差，国外学者称之位共模误差(Commonmodeerror，CME)，被广泛认为是GPS时间序列数据误差的主要来源之一(Wdowinsk，1997；Dong，2006，2013；Jiang，2013；Shen，2014)，CME对GPS的影响较大，是影响GPS解精度及可靠性的主要因素之一。传统的CME滤波方法假定共模误差的空间响应存在一致性，即均匀分布，使得其在大区域GPS网中应用受限，甚至分离出错误的共模误差分量，对一些地球物理现象给出错误的解释。

目前对GPS网中共模误差的分离研究方法存在一些缺陷：1)什么是CME的物理起源，GPS网中的CME的影响因子包含哪些，缺乏深入的研究；2)目前大部分的CME研究局限在小区域网内，CME空间上相似性在多大的尺度上表现出一致性，直接影响CME分离的可靠性。

发明内容

针对现有技术存在的问题，本发明考虑到已有方法的局限性(区域性)，进一步对共模误差的空间分布进行探讨，将相关系数、距离、经纬度、本地效应(WRMS)、负荷效应作为评价因子，通过聚类分析，进一步探讨共模误差提取及分离方法，提供了一种GPS时间序列广义共模误差提取方法，高效、准确的分离出共模误差，提高GPS坐标时间序列的精度及可靠性。

本发明所采用的技术方案是：一种GPS时间序列广义共模误差提取方法，其特征在于，包括以下步骤：

步骤1：针对GPS观测值及相关文件，采用多种解算算法进行解算，分别获取GPS测站单日松弛解，通过公共基站进行不同解加权进行联合解算，获得GPS测站坐标时间序及速度参数；所述相关文件包括星历文件和表文件；

步骤2：对获取的GPS测站坐标时间序列进行滤波处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的测站坐标时间序列建立以下残差时间序列模型：

$y{\left(t\right)}_{E/N/U}=a+bt+c\sin \left(2\mathrm{\π}t\right)+d\cos \left(2\mathrm{\π}t\right)+e\sin \left(4\mathrm{\π}t\right)+f\cos \left(4\mathrm{\π}t\right)+\underset{j=1}{\overset{n_g}{\Σ}}{g}_{j}H(t-T_{gj})+v_i$

其中：y(t)为时刻t对应的GPS测站坐标观测值，包含E、N、U三个坐标分量；t_i(i＝1…,n)为GPS站点单日历元，以年为单位；a为GPS测站位置，为序列的平均值；b为线性速度，即趋势项；系数c、d、e、f为年周期和半年周期项的系数，c、d、e、f为待估计参数，经拟合获得；为跳变改正项，表示跳变振幅，g表示发生在历元T_g处的由于各种原因引起的阶跃式偏移量，表示跳变个数，j为跳变编号，这里假定发生偏移的时刻T_g已知，H为海维西特阶

梯函数，在跳变前H值为0，跳变后H值为1，v_i为时刻t的观测值残差；

步骤3：对获取的GPS测站坐标时间序列继续进行滤波处理，包括趋势项扣除和构造运动产生的跳跃项扣除；所述构造运动产生的跳跃项进行扣除，只包含周年、半年项，从而得到较高精度的坐标时间序列；

步骤4：计算GPS测站坐标时间序列E、N、U三坐标分量的各测站之间互相关系数及其加权均方根；所述相关系数包括互相关系数和自相关系数；

对步骤3中GPS测站坐标时间序列进行时间序列分析，得到各GPS测站站点的加权均方根；将加权均方根最大的点或者其值超过GPS网中所有点的加权均方根的中误差的2倍限差的点，认为是包含强烈的本地效应即噪声在进行滤波处理之前进行去除，防止将本地效应混叠到共模误差中，去除的这类点称之为本地效应点；

步骤5：站点空间分布图形绘制

通过采用地学绘图软件GMT根据站点的经纬度、站点信息计算站点之间的距离，并在图形中以中心圆方式进行站点图绘制；绘制得到的站点空间分布图，作为共模误差区域划分的评价因子之一；

步骤6：计算地表负荷对GPS测站位移影响；

采用mload程序分别计算大气压、非潮汐海洋、积雪和土壤水负荷引起的台站位移，通过负载改正，提高坐标序列的精度，消除部分非构造信号；同时通过分析地表负载的空间响应，为空间滤波的区域划分提供评价因子作为依据；

步骤7：将步骤4～6中得到的各站点的时空相关系数、本地效应点、经纬度、距离、地表负荷效应，作为评价因子，通过聚类分析，初步对GPS网进行子区域划分；

步骤8：采用主成分分析法对GPS测站坐标时间序列进行共模误差分离，提取共模误差。

作为优选，步骤4中所述互相关系数的公式如下：

${\mathrm{\ρ}}_{xy}=\frac{\underset{i=1}{\overset{N}{\Σ}}{x}_i{y}_i}{\sqrt{\underset{i=1}{\overset{N}{\Σ}}{x}_i{x}_i}\sqrt{\underset{i=1}{\overset{N}{\Σ}}{y}_i{y}_i}}$

其中x_i、y_i为时刻t对应的GPS测站坐标时间序列E、N、U三坐标分量的位移序列。

作为优选，步骤7中所述对GPS网进行子区域划分，主要依据为：

(1)经纬度进行初步区域划分，区域划分通过步骤5进行；

(2)根据步骤4获得的E、N、U三坐标分量的相关系及其加权均方根进行区域划分是否正确的检核；要求子网内站点的互相关系数大于0.15，将加权均方根最大的点或者其值超过GPS网中所有点的加权均方根的中误差的2倍限差的点，认为是包含强烈的本地效应即噪声予以进行剔除，降低本地噪声的影响；

(3)步骤6中计算地表负荷对GPS测站位移影响，把不同区域的负载效应作为区域划分的一个决策因子。

作为优选，步骤8中所述采用主成分分析法对GPS测站坐标时间序列进行共模误差分离，其具体实现过程为：

假定GPS台站获得的三维坐标观测值形成一个n×m的数据矩阵X，其中n>m，n为观测数或历元数，m为观测类型，其协方差阵为C_X，则C_X＝X^TX；数据矩阵如下：

$X=\left[\begin{array}{cccc}{x}_1\left(t_1\right)& x_1\left(t_2\right)& \mathrm{...}& x_1\left(t_m\right)\\ x_2\left(t_1\right)& x_2\left(t_2\right)& \mathrm{...}& x_2\left(t_m\right)\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \\ x_{1n-1}\left(t_1\right)& x_{n-1}\left(t_2\right)& \mathrm{...}& x_{n-1}\left(t_m\right)\\ x_n\left(t_1\right)& x_n\left(t_2\right)& \mathrm{...}& x_n\left(t_m\right)\end{array}\right];$

其中：m×1维列向量为其协方差阵的特征向量，λ_i为对应的特征值，令其中σ_i为正的奇异值，i＝1,2…r，则有：

$\left(X^{T}X\right)\stackrel{\→}{u_i}={\mathrm{\λ}}_i\stackrel{\→}{v_i};$

假定：

$\begin{array}{ccc}\stackrel{\→}{u_i}=\frac{1}{{\mathrm{\σ}}_i}X\stackrel{\→}{v_i},& \mathrm{\Λ}={\Σ}^T\Σ,& \Σ=diag({\mathrm{\σ}}_1,{\mathrm{\σ}}_2\mathrm{...}{\mathrm{\σ}}_r,\mathrm{0...0});\end{array}$

$\begin{array}{cc}V=\[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& \mathrm{...}& \stackrel{\→}{v_m}\end{array}\],& U=\left[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& \mathrm{...}& \stackrel{\→}{u_n}\end{array}\right];\end{array}$

其中是n×1列向量，U为n×n向量矩阵，V为m×m向量矩阵，则有：

X＝UΣV^T；

C_X＝VΛV^T；

即V构成X的正交基底，矩阵X展开可得：

$X(t_i,x_j)=\underset{k=1}{\overset{n}{\Σ}}{a}_k\left(t_i\right)v_k\left(x_j\right);$

a_k(t_i)可由下式求出：

$X=AV\⇒A={\mathrm{XV}}^{-1}\⇒A={\mathrm{XV}}^T;$

$a_k\left(t_i\right)=\underset{j=1}{\overset{n}{\Σ}}X(t_i,x_j)v_k\left(x_j\right);$

式中a_k(t)是第k个主成分，v_k(x)是对应主成分的响应特征矩阵，分别代表时间特征和空间响应，取前k个主分量计算得到的共模误差为：

${\mathrm{\ϵ}}_i\left(t_i\right)=\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(t_i\right){\mathrm{\υ}}_k\left(x_j\right).$

与现有的技术相比，本发明具有特点：

本发明的创新之处在于，一方面，对共模误差的空间响应进行了分析，建立了大区域GPS网共模误差的提取机制，突破了传统方法中大区域网下如何进行子网划分，引入了相关系数ρ、距离及经纬度、本地效应、负荷效应、主分量贡献率及其空间响应作为评价因子，通过聚类分析为大区域、大尺度下GPS网站点间共模误差，提供了可行的估计方法，获得更好的滤波效果。另一方面，顾及了共模误差的周期性，与传统方法相比，在进行主成分分析滤波之前，保留了周年、半年项，考虑了共模误差的周年性，同时降低了周年项拟合过程中引入的模型误差，有助于真实反映共模误差的空间变化及周期性变化，为进一步提高GPS坐标序列模型的精度提供依据。

附图说明

图1是本发明实施例的流程图；

图2是本发明实施例的GPS测站原始时间序列图和仅留周年、半年项的残差序列图；其中图2a为GPS测站原始时间序列图，图2b为仅留周年、半年项的残差序列图；

图3是本发明实施例的站点椭球面距离计算涉及的站点表述及选择依据；

图4本发明实施例的站点相关系统计结果；

图5是本发明实施例的测站加权均方根的全球分布图；

图6是本发明实施例的不同负载的位移效应的比较图；

图7是本发明实施例的PCA主分量的空间特征响应结果；

图8是本发明实施的广义共模误差提取方法与传统滤波方法的比较结果，其中图a是本发明方法结果，图b为传统方法的残差序列；

图9是本发明实施例经过共模误差提出后的功率频谱分析结果。

具体实施方式

为了使本发明的目的、技术方案和有益效果更加清楚明白，下面结合附图及具体实施方式，进一步说明本发明。应当理解，一下描述的集体实施方式仅用于解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种GPS时间序列广义共模误差提取方法，包括以下步骤：

步骤1：针对GPS观测值及相关文件(星历文件、表文件等)，采用多种解算算法进行解算，分别获取GPS测站单日松弛解，通过公共基站进行不同解加权进行联合解算，获得GPS测站坐标时间序列；

本步骤属于现有技术，具体可通过现有技术中成熟的GAMIT/GLOBK、Bernese、GIPSY等高精度数据处理软件或IGS分析中心获取数据。不同数据处理软件由于算法不的完善、模型系统偏差等往往会引入不可避免的解算误差，本发明的新颖之处在于才用了多种解算软件(GAMIT、GIPSY等)进行解算，通过公共基站进行不同解加权进行联合解算，能有效的消除单一软件解算的模型系统误差，进一步提高解的可靠性。

步骤2：对获取的GPS测站坐标时间序列进行预处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的测站坐标时间序列(E、N、U三分量)建立以下残差时间序列模型：

$y{\left(t\right)}_{E/N/U}=a+bt+csin\left(2\mathrm{\π}t\right)+dcos\left(2\mathrm{\π}t\right)+esin\left(4\mathrm{\π}t\right)+fcos\left(4\mathrm{\π}t\right)+\underset{j=1}{\overset{n_g}{\Σ}}{g}_{j}H(t-T_{gj})+v_i$

其中：y(t)为时刻t对应的GPS测站坐标观测值，包含E、N、U三个坐标分量；t_i(i＝1…,n)为GPS站点单日历元，以年为单位；a为GPS测站位置，为序列的平均值；b为线性速度，即趋势项；系数c、d、e、f为年周期和半年周期项的系数，c、d、e、f为待估计参数，经拟合获得；为跳变改正项，表示跳变振幅，g表示发生在历元T_g处的由于各种原因引起的阶跃式偏移量，表示跳变个数，j为跳变编号，这里假定发生偏移的时刻T_g已知，H为海维西特阶

梯函数，在跳变前H值为0，跳变后H值为1，v_i为时刻t的观测值残差；

步骤3：对获取的GPS测站坐标时间序列继续进行预处理，包括趋势项扣除和构造运动产生的跳跃(jump)项扣除；所述构造运动产生的跳跃项进行扣除，只包含周年、半年项，从而得到较高精度的坐标时间序列；

与传统方法不同的是，考虑到共模误差可能混叠在周年、半年项里面，在对GPS时间序列进行共模误差分离时，若事先扣除周年半年项，会是的共模误差的分离存在缺陷，因此本发明实施的新颖之处在于残差序列中保留了周年、半年项，顾及了共模误差的周期性，图2a为获得的原始坐标时间序列，图2b为经过去趋势等处理后的残差序列。

步骤4：计算GPS测站坐标时间序列E、N、U三坐标分量的时空相关系数及其加权均方根；所述时空相关系数包括互相关系数和自相关系数；

GPS测站坐标序列之间互相关系数的公式如下：

${\mathrm{\ρ}}_{xy}=\frac{\underset{i=1}{\overset{N}{\Σ}}{x}_i{y}_i}{\sqrt{\underset{i=1}{\overset{N}{\Σ}}{x}_i{x}_i}\sqrt{\underset{i=1}{\overset{N}{\Σ}}{y}_i{y}_i}}$

其中x_i、y_i为时刻t对应的GPS测站坐标时间序列E、N、U三坐标分量的位移序列。

另外通过MATLAB工具编写互相关计算程序，通过计算并绘制单变量随机时间序列的样本的ACF及置信区间，可以快速求得某两个时间序列的互相关系数，然后用polyfit、polyval这两个函数可以对画出的相关系数散点图进行拟合分析，得到这些散点的趋势图(请见图4)，在此基础上可以对其进行分析。

对步骤3中GPS测站坐标时间序列进行分析，得到各站点的加权均方根(WRMS)，对GPS坐标时间序列，其WRMS定义为：dsqrt[sum(resi*weight*resi)/sum(weight)]，单位是米(或毫米)，图5为GPS测站站点WRMS趋势图。

步骤5：站点空间分布图形绘制；

根据步骤3提供的先验信息，选取若干相关系数较好及WRMS较小的点作为中心，采用地学绘图软件GMT(TheGenericMappingTools)根据站点的经纬度、站点信息计算站点之间的距离，并在图形中以中心圆方式进行站点图绘制，初步将大尺度GPS网划分成若干子网，绘制得到的站点空间分布图，作为共模误差区域划分的评价因子之一；图3为站点空间分布实例，以A为中心，以固定半径为插值，通过等值线进行椭球面图形绘制，可以直观的表述站点之间的空间相对关系，为大区域子网划分提供判定依据。

步骤6：计算地表负荷对GPS测站位移影响；

采用mload程序分别计算大气压、非潮汐海洋、积雪和土壤水负荷引起的台站位移，通过负载改正，提高坐标序列的精度，消除部分非构造信号，避免将地表负载效应纳入共模误差；同时通过分析地表负载的空间响应(请见图6)，为空间滤波的区域划分是否合理提供决策因子。

步骤7：将步骤4～6中得到的各站点的时空相关系数、本地效应点、经纬度、距离、地表负荷效应，作为评价因子，通过聚类分析，初步对GPS网进行子区域划分；

对GPS网进行子区域划分，主要依据为：

(1)经纬度进行初步区域划分，区域划分通过步骤5进行；

(2)根据步骤4获得的E、N、U三坐标分量的相关系及其加权均方根进行区域划分是否正确的检核；要求子网内站点的互相关系数大于0.15(经验值，2000公里之后，站点互相关系数小于0.15)，将加权均方根较大的点作为本地效应进行剔除，降低本地噪声的影响；

(3)步骤6中计算地表负荷对GPS测站位移影响，把不同区域的负载效应作为区域划分的一个决策因子。

通过图6统计分析结果发现，本发明所划分的子网的负载效应存在明显的差异。

步骤8：采用主成分分析法(PCA方法)对GPS测站坐标时间序列进行共模误差分离，提取共模误差。

为了描述PCA方法在GPS时间序列中的具体实现，假定GPS台站获得的三维坐标观测值形成一个n×m(n>m，n为观测数或历元数，m为观测类型)的数据矩阵X，其协方差阵为C_X，则C_X＝X^TX。数据矩阵如下：

$X=\left[\begin{array}{cccc}{x}_1\left(t_1\right)& x_1\left(t_2\right)& \mathrm{...}& x_1\left(t_m\right)\\ x_2\left(t_1\right)& x_2\left(t_2\right)& \mathrm{...}& x_2\left(t_m\right)\\ \mathrm{...}& \mathrm{...}& \mathrm{...}& \\ x_{1n-1}\left(t_1\right)& x_{n-1}\left(t_2\right)& \mathrm{...}& x_{n-1}\left(t_m\right)\\ x_n\left(t_1\right)& x_n\left(t_2\right)& \mathrm{...}& x_n\left(t_m\right)\end{array}\right]$

其中：(m×1维列向量)为其协方差阵的特征向量，λ_i为对应的特征值，令其中σ_i为正的奇异值，i＝1,2…r。则有：

假定

$\begin{array}{ccc}\stackrel{\→}{u_i}=\frac{1}{{\mathrm{\σ}}_i}X\stackrel{\→}{v_i},& \mathrm{\Λ}={\Σ}^T\Σ,& \Σ=diag({\mathrm{\σ}}_1,{\mathrm{\σ}}_2\mathrm{...}{\mathrm{\σ}}_r,\mathrm{0...0})\end{array}$

$\begin{array}{cc}V=\[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& \mathrm{...}& \stackrel{\→}{v_m}\end{array}\],& U=\left[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& \mathrm{...}& \stackrel{\→}{u_n}\end{array}\right]\end{array}$

其中是n×1列向量，U为n×n向量矩阵，V为m×m向量矩阵，则有：

X＝UΣV^T

C_X＝VΛV^T

即V构成X的正交基底，矩阵X展开可得：

$X(t_i,x_j)=\underset{k=1}{\overset{n}{\Σ}}{a}_k\left(t_i\right)v_k\left(x_j\right)---\left(9\right)$

a_k(t_i)可由下式求出：

$X=AV\⇒A={\mathrm{XV}}^{-1}\⇒A={\mathrm{XV}}^T---\left(10\right)$

$a_k\left(t_i\right)=\underset{j=1}{\overset{n}{\Σ}}X(t_i,x_j)v_k\left(x_j\right)---\left(11\right)$

式中a_k(t)是第k个主成分，v_k(x)是对应主成分的响应特征矩阵，分别代表时间特征和空间响应，取前k个主分量计算得到的共模误差为：

${\mathrm{\ϵ}}_{i^{.}}\left(t_i\right)=\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(t_i\right){\mathrm{\υ}}_k\left(x_j\right).$

不同于以往方法，本发明估计了站点之间的相关系，考虑了站点的本地效应，将WRMS较大的测站进行去除，同时把站点的地理因素(经纬度及距离、外部负载变化)作为区域划分的因子，考虑的共模误差的空间不一致性，提出了依据地理环境因素的、多尺度评价因子的广义共模误差提取方法，经过主成分分析后，取前k个主分量(累计贡献了超过80％，若前k个主分量的累积贡献率小于80％，说明残差序列的相关性较差，或者不具有相关系，即不能较好的探测出共模误差，子区域划分不合理)计算得到的共模误差为：

${\mathrm{\ϵ}}_i\left(t_i\right)=\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(t_i\right){\mathrm{\υ}}_k\left(x_j\right)$

下面结合对比试验进一步说明本发明的有益效果。

为了对比传统空间滤波方法(Wdowinski，1997)与本发明方法处理共模误差情况，对大区域内的22个测站进行分析，表1为传统方法和本发明方法划归子区域后三坐标分量前4个主分量累积贡献率结果：

表1

从表1可知，本发明方法前4个主分量的贡献率均大于80％，即前几个主分量能较好的反映残差序列的主要性质，且前4个主分量的贡献率高于传统方法，E、N、U方向的累积贡献率约搞9.38％；此外，对主分量的空间响应分析，发现子网1、2的空间响应存在明显的差异，如子网1最多响应在E分量，而子网2出现在U分量，而子网内的空间响应一致性较好(请见图7)，说明共模误差在空间尺度上存在差异。

请见图8，为本发明方法与传统方法进行比较的结果，本发明所提出的方法能有效地提高GPS测站时间序列的精度，同时图9给出了滤波后残差序列的功率谱分析结果，结果表明经过滤波后，广泛存在的GPS交点年信号的振幅也得到了极大的抑制。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (4)

1.一种GPS时间序列广义共模误差提取方法，其特征在于，包括以下步骤：

步骤1：针对GPS观测值及相关文件，采用多种解算算法进行解算，分别获取GPS测站单日松弛解，通过公共基站进行不同解加权进行联合解算，获得GPS测站坐标时间序及速度参数；所述相关文件包括星历文件和表文件；

步骤2：对获取的GPS测站坐标时间序列进行滤波处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的测站坐标时间序列建立以下残差时间序列模型：

$y{\left(t\right)}_{E/N/U}=a+bt+c\sin \left(2\mathrm{\π}t\right)+d\cos \left(2\mathrm{\π}t\right)+e\sin \left(4\mathrm{\π}t\right)+f\cos \left(4\mathrm{\π}t\right)+\underset{j=1}{\overset{n_g}{\Σ}}{g}_{j}H\left(t-T_{gj}\right)+v_i;$

其中：y(t)为时刻t对应的GPS测站坐标观测值，包含E、N、U三个坐标分量；t_i(i＝1…,n)为GPS站点单日历元，以年为单位；a为GPS测站位置，为序列的平均值；b为线性速度，即趋势项；系数c、d、e、f为年周期和半年周期项的系数，c、d、e、f为待估计参数，经拟合获得；为跳变改正项，表示跳变振幅，g表示发生在历元T_g处的由于各种原因引起的阶跃式偏移量，表示跳变个数，j为跳变编号，这里假定发生偏移的时刻T_g已知，H为海维西特阶梯函数，在跳变前H值为0，跳变后H值为1，v_i为时刻t的观测值残差；

步骤3：对获取的GPS测站坐标时间序列继续进行滤波处理，包括趋势项扣除和构造运动产生的跳跃项扣除；所述构造运动产生的跳跃项进行扣除，只包含周年、半年项，从而得到较高精度的坐标时间序列；

步骤4：计算GPS测站坐标时间序列E、N、U三坐标分量的各测站之间互相关系数及其加权均方根；所述相关系数包括互相关系数和自相关系数；

对步骤3中GPS测站坐标时间序列进行时间序列分析，得到各GPS测站站点的加权均方根；将加权均方根最大的点或者其值超过GPS网中所有点的加权均方根的中误差的2倍限差的点，认为是包含强烈的本地效应即噪声在进行滤波处理之前进行去除，防止将本地效应混叠到共模误差中，去除的这类点称之为本地效应点；

步骤5：站点空间分布图形绘制；

通过采用地学绘图软件GMT根据站点的经纬度、站点信息计算站点之间的距离，并在图形中以中心圆方式进行站点图绘制；绘制得到的站点空间分布图，作为共模误差区域划分的评价因子之一；

步骤6：计算地表负荷对GPS测站位移影响；

采用mload程序分别计算大气压、非潮汐海洋、积雪和土壤水负荷引起的台站位移，通过负载改正，提高坐标序列的精度，消除部分非构造信号；同时通过分析地表负载的空间响应，为空间滤波的区域划分提供评价因子作为依据；

步骤7：将步骤4～6中得到的各站点的时空相关系数、本地效应点、经纬度、距离、地表负荷效应，作为评价因子，通过聚类分析，初步对GPS网进行子区域划分；

步骤8：采用主成分分析法对GPS测站坐标时间序列进行共模误差分离，提取共模误差。

2.根据权利要求1所述的GPS时间序列广义共模误差提取方法，其特征在于，步骤4中所述互相关系数的公式如下：

${\mathrm{\ρ}}_{xy}=\frac{\underset{i=1}{\overset{N}{\Σ}}{x}_i{y}_i}{\sqrt{\underset{i=1}{\overset{N}{\Σ}}{x}_i{x}_i}\sqrt{\underset{i=1}{\overset{N}{\Σ}}{y}_i{y}_i}}$

其中x_i、y_i为时刻t对应的GPS测站坐标时间序列E、N、U三坐标分量的位移序列。

3.根据权利要求1所述的GPS时间序列广义共模误差提取方法，其特征在于：步骤7中所述对GPS网进行子区域划分，主要依据为：

(1)经纬度进行初步区域划分，区域划分通过步骤5进行；

(2)根据步骤4获得的E、N、U三坐标分量的相关系及其加权均方根进行区域划分是否正确的检核；要求子网内站点的互相关系数大于0.15，将加权均方根最大的点或者其值超过GPS网中所有点的加权均方根的中误差的2倍限差的点，认为是包含强烈的本地效应即噪声予以进行剔除，降低本地噪声的影响；

(3)步骤6中计算地表负荷对GPS测站位移影响，把不同区域的负载效应作为区域划分的一个决策因子。

4.根据权利要求1所述的GPS时间序列广义共模误差提取方法，其特征在于：步骤8中所述采用主成分分析法对GPS测站坐标时间序列进行共模误差分离，其具体实现过程为：

假定GPS台站获得的三维坐标观测值形成一个n×m的数据矩阵X，其中n>m，n为观测数或历元数，m为观测类型，其协方差阵为C_X，则C_X＝X^TX；数据矩阵如下：

$X=\left[\begin{array}{cccc}{x}_1\left(t_1\right)& x_1\left(t_2\right)& ...& x_1\left(t_m\right)\\ x_2\left(t_1\right)& x_2\left(t_2\right)& ...& x_2\left(t_m\right)\\ ...& ...& ...& \\ x_{1n-1}\left(t_1\right)& x_{n-1}\left(t_2\right)& ...& x_{n-1}\left(t_m\right)\\ x_n\left(t_1\right)& x_n\left(t_2\right)& ...& x_n\left(t_m\right)\end{array}\right];$

其中：m×1维列向量为其协方差阵的特征向量，λ_i为对应的特征值，令其中σ_i为正的奇异值，i＝1,2…r，则有：

$\left(X^{T}X\right)\stackrel{\→}{u_i}={\mathrm{\λ}}_i\stackrel{\→}{v_i};$

假定：

$\stackrel{\→}{u_i}=\frac{1}{{\mathrm{\σ}}_i}X\stackrel{\→}{v_i},\mathrm{\Λ}={\mathrm{\Σ}}^T\mathrm{\Σ},\mathrm{\Σ}=diag({\mathrm{\σ}}_1,{\mathrm{\σ}}_2...{\mathrm{\σ}}_r,0...0);$

$V=\left[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& ...& \stackrel{\→}{v_m}\end{array}\right],U=\left[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& ...& \stackrel{\→}{u_n}\end{array}\right];$

其中是n×1列向量，U为n×n向量矩阵，V为m×m向量矩阵，则有：

X＝UΣV^T；

C_X＝VΛV^T；

即V构成X的正交基底，矩阵X展开可得：

$X(t_i,x_j)=\underset{k=1}{\overset{n}{\Σ}}{a}_k\left(t_i\right)v_k\left(x_j\right);$

a_k(t_i)可由下式求出：

$X=AV\⇒A={\mathrm{XV}}^{-1}\⇒A={\mathrm{XV}}^T;$

$a_k\left(t_i\right)=\underset{j=1}{\overset{n}{\Σ}}X(t_i,x_j)v_k\left(x_j\right);$

式中a_k(t)是第k个主成分，v_k(x)是对应主成分的响应特征矩阵，分别代表时间特征和空间响应，取前k个主分量计算得到的共模误差为：

${\mathrm{\ϵ}}_i\left(t_i\right)=\underset{k=1}{\overset{p}{\Σ}}{a}_k\left(t_i\right)v_k\left(x_j\right).$