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

Abstract

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

Description

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

技术领域

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

背景技术

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

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

在不同站点坐标时间序列(GPS网)中存在一种空间相关的误差，国外学者称之位共模误差(Common mode error，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测站坐标时间序列进行滤波处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的GPS测站坐标时间序列建立以下残差时间序列模型：

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

其中：y(t_i)为时刻t_i对应的GPS测站坐标观测值，包含E、N、U三个坐标分量，y(t_i)_E/N/U表示ENU不同方向的坐标序列；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{...}& \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};$

假定：

Λ＝Σ^TΣ，Σ＝diag(σ₁,σ₂…σ_r,0…0)；

$V=\left[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& \mathrm{...}& \stackrel{\→}{v_m}\end{array}\right],U=\left[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& \mathrm{...}& \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){\mathrm{\υ}}_k\left(x_j\right);$

其中p为主分量的个数。

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

本发明的创新之处在于，一方面，对共模误差的空间响应进行了分析，建立了大区域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测站坐标时间序列进行预处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的GPS测站坐标时间序列(E、N、U三分量)建立以下残差时间序列模型：

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

其中：y(t_i)为时刻t_i对应的GPS测站坐标观测值，包含E、N、U三个坐标分量，y(t_i)_E/N/U表示ENU不同方向的坐标序列；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(The Generic Mapping Tools)根据站点的经纬度、站点信息计算站点之间的距离，并在图形中以中心圆方式进行站点图绘制，初步将大尺度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{...}& \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。则有：假定

Λ＝Σ^TΣ，Σ＝diag(σ₁,σ₂…σ_r,0…0)

$V=\left[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& \mathrm{...}& \stackrel{\→}{v_m}\end{array}\right],U\left[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& \mathrm{...}& \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)---\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);$

其中p为主分量的个数。

不同于以往方法，本发明估计了站点之间的相关系，考虑了站点的本地效应，将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测站坐标时间序列进行滤波处理，包括粗差探测、阶跃探测、粗差修复、阶跃修复；并对获取的GPS测站坐标时间序列建立以下残差时间序列模型：

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

其中：y(t_i)为时刻t_i对应的GPS测站坐标观测值，包含E、N、U三个坐标分量，y(t_i)_E/N/U表示ENU不同方向的坐标序列；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)& \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};$

假定：

Λ＝Σ^TΣ，Σ＝diag(σ₁,σ₂…σ_r,0…0)；

$V=\[\begin{array}{cccc}\stackrel{\→}{v_1}& \stackrel{\→}{v_2}& \mathrm{...}& \stackrel{\→}{v_m}\end{array}\],U=\[\begin{array}{cccc}\stackrel{\→}{u_1}& \stackrel{\→}{u_2}& \mathrm{...}& \stackrel{\→}{u_n}\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);$

其中p为主分量的个数。