CN106772498A - 一种gps位置时间序列噪声模型建立方法 - Google Patents

Abstract

本发明公开了一种GPS位置时间序列噪声模型建立方法，对GPS测站坐标位置时间序列的噪声模型特性进行分析，建立GPS测站坐标位置时间序列噪声模型分析机制，考虑了GPS位置序列解算策略、时间跨度、粗差、趋势速度项、地表环境负载效应、共模误差等因素对噪声模型建立的影响；本发明考虑了传统GPS时间序列噪声模型估计方法的局限性，顾及了噪声模型的多样性，考虑了多模型及其组合模型，且采用改进的AIC/BIC模型估计准则。通过对多因素进行分析，最终建立了一种稳健的GPS测站坐标位置时间序列噪声模型估计方法，以获得准确的噪声模型估计结果，最终获得准确的站速度及其不确定性参数，为进一步提高GPS坐标序列在高精度地球动力学中的应用提供可靠性数据基础。

Description

一种GPS位置时间序列噪声模型建立方法

技术领域

本发明属于卫星导航技术领域，涉及一种稳健的GPS位置时间序列噪声模型建立方法。

背景技术

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

GPS观测值噪声最初被假定为纯白噪声，即GPS连续观测值(坐标序列)是相互独立的，忽略了测站之间的相关性。通过对长周期的时间序列研究分析发现，该假设会使得数据处理结果产生偏差，使得参数估计带来偏差，表现为导致速度不确定度的过低估计，即得到过于乐观的误差估计结果。因此在利用GPS探索地球动力学等高精度问题时，有必要采用较为严格、准确的误差模型。

闪烁噪声加白噪声、幂律噪声及其组合被认为是描述GPS坐标时间序列的最佳随机特性的噪声模型，分析其原因可能有两方面：第一，已有的研究进行噪声模型估计时，采用的组合模型较少，即没有顾及噪声模型的多样性。大部分的学者仅仅采用FN、WN、FN+WN、PL+WN等噪声模型进行分析，并没有考虑更加广泛的噪声模型；第二，进行噪声估计时由于受到GPS坐标时间序列长度限制，未能采用长周期的时间序列进行分析，对GPS坐标时间序列中的长周期噪声分量，不能准确的进行估计。第三，外部环境因素对测站时间序列的潜在影响，如不同观测墩类型、时间序列数据缺失、不同跨度的时间序列、不同地理环境等可能导致噪声模型的改变；第四，噪声模型估计准则的可靠性，根据极大似然估计原理，不同的噪声模型组合将得到不同的极大似然对数值，MLE的数值越大，其结果越可靠。然而，当噪声模型包含的未知参数越多，其MLE值越大。因此不能简单选择MLE值较大的模型作为最优噪声模型，噪声模型估计算法有待进一步改进。因此，GPS坐标时间序列噪声模型的无偏估计，探讨噪声模型的稳健估计方法，对高精度GPS应用，尤其是速度场应用方面具有重要意义。

发明内容

针对现有技术存在的问题，本发明考虑到已有方法的局限性，对GPS测站坐标位置时间序列的噪声模型特性进行分析，建立GPS测站坐标位置时间序列噪声模型分析机制，考虑了GPS位置序列解算策略、时间跨度、粗差、趋势速度项(Trend)、地表环境负载效应、共模误差等因素对噪声模型建立的影响，提供了可行的噪声模型建立方法，以获得稳健的噪声模型估计结果，获取准确的站速度及其不确定性。

本发明所采用的技术方案是：一种GPS位置时间序列噪声模型建立方法，其特征在于，包括以下步骤：

步骤1：对原始观测文件进行解算，对不同松弛解进行加权解算，获取GPS测站坐标位置时间序列(E、N、U三坐标分离)；

步骤2：对获取的GPS测站坐标位置时间序列进行预处理；

步骤3：对步骤2中预处理后的GPS测站坐标位置时间序列继续进行地表环境负载修正，以改正地表环境负载引起的GPS测站非线性变化；

步骤4：对步骤3中修正后的GPS测站坐标位置时间序列继续进行共模误差改正，去除坐标序列中的共模误差，获得残差序列。为了准确的分离共模误差，本发明专利考虑的共模误差的时空特征，即对于大尺度下的共模误差分离，以地表环境负载为参考，对大区域进行子网划分；并顾及GPS坐标序列、共模误差的周期性，在采用主成分分析法去除共模误差之前，保留坐标序列的周年、半周年项，以避免周期性共模分量的残留；

步骤5：采用极大似然估计方法对不同组合的单一噪声模型、组合噪声模型分别对GPS测站坐标位置时间序列三坐标分量进行噪声估计，获得不同组合噪声模型在GPS测站坐标位置时间序列三坐标分量的MLE值；

步骤6：根据单一噪声模型、组合噪声模型GPS测站坐标位置时间序列三坐标分量的MLE值获得GPS测站坐标位置时间序列坐标分量的最优噪声模型。

作为优选，步骤2中所述预处理包括粗差分析、去均值、去趋势项，探测Offset并改正。

作为优选，步骤3中所述负载修正包括大气、地表水、非潮汐海洋、积雪负载改正。

作为优选，步骤6的具体实现包括以下子步骤：

步骤6.1：采用极大似然估计方法解算得到不同噪声模型的MLE值，并根据不同噪声模型的参数个数k，进一步采用AIC、BIC方法确定最佳噪声模型，并获得GPS站速度及其不确定性的估值；

步骤6.2：采用AIC、BIC模型估计准则分别计算不同噪声模型的AIC、BIC值，AIC/BIC的求解原理如下：

AIC＝-2lnL+2k；

BIC＝-2·lnL+k·ln(n)；

式中L为某一噪声模型下的似然函数值MLE，n为观测值个数，k为噪声模型的估计参数个数；

步骤6.3：当ACI、BIC估计出的GPS测站坐标位置时间序列噪声模型结果一致时，接受该估计结果；当ACI、BIC估计出的结果不一致时，采用功率谱分析的方法结合极大似然估计的结果，即通过功率谱分析方法分别对AIC、BIC估计后的结果以及原始数据进行频谱分析，并绘制频谱图，通过对不同方法估计出的频谱图与原始观测数据的频谱图的拟合趋近程度，进一步确定最优模型。

作为优选，步骤4中所述共模误差改正，是对于大尺度(500km尺度区域以上)下的共模误差分离，以地表环境负载为参考，对大区域进行子网划分；并顾及GPS坐标序列、共模误差的周期性，在采用主成分分析法去除共模误差之前，保留坐标序列的周年、半周年项，以避免周期性共模分量的残留。

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

本发明的创新之处在于，一方面，考虑了GPS位置序列解算策略(有效地减弱了解算软件及其算法模型对噪声模型的影响)、时间跨度(通过采用10年以上的时间跨度，提高了估计结果的无偏性)、粗差影响(减弱了粗差对噪声估计结果的影响)、趋势速度项(只对已知的过程进行trend去除，如冰川均衡调整、板块旋转引起的趋势变化，考虑了噪声的trend变化)、地表环境负载效应与共模误差改正(顾及了测站的非线性变化因素，并对其进行了改正)等因素对噪声模型建立的影响，通过对上述因素进行分析，建立了一个完整的、稳健的噪声模型估计方法；另一方面，本发明考虑了传统GPS时间序列噪声模型估计方法的局限性，顾及了噪声模型的多样性，考虑了多模型及其组合模型，且采用改进的AIC/BIC模型估计准则，最终建立了一种稳健的GPS测站坐标位置时间序列噪声模型建立方法，以获得准确的噪声模型估计结果，获取准确的站速度及其不确定性参数，为进一步提高GPS坐标序列在高精度地球动力学中的应用提供可靠性数据基础。

附图说明

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

图2是本发明实施例的噪声模型估计结果图。

具体实施方式

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

请见图1，本发明提供的一种GPS位置时间序列噪声模型建立方法，包括以下步骤：

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

获取GPS站坐标位置时间序列，为了保证最终用于噪声模型及站速度估计的位置序列的可靠性，降低不同解算模型、软件系统误差的影响，采用双差与精密单点定位方法分别解算原始观测值，再采用QOCA软件进行联合求解，联合解可以吸取不同软件及解算策略的优势，消除单一软件解算过程中引入的系统偏差，改善解的时空分布，提高坐标序列的精度；另外，为了保证噪声模型估计结果的准确性，原则上要求坐标序列尽量长于10年；

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

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

与传统方法不同的是，考虑了粗差、Offset、构造运动对噪声模型的影响。其中，粗差探测采用四分位数间距法(interquartile range,IQR)，四分位数间距由P25、P50、P75将一组变量值等分为四部分，P25称下四分位数(Q1)，P75称上四分位数(Q3)，将P75与P25之差定义为四分位数间距(IQR)。分别计算a(Q1-3*IQR)、b(Q3+3*IQR)的值，原始序列中位于(a,b)区间之外的值，则为粗差。

对于GPS测站坐标位置时间序列中存在的Offset，采用如下方法进行改正：1)对已知的Offset，根据IGS发布的相关测站Offset资料，即Offset发生的时刻及影响(见ftp://sopac-ftp.ucsd.edu/pub/gamit/setup/siteOffsets.txt)进行改正；2)对于未知的Offset，采用Sequential t test Analysis of Regime Shifts(STARS)方法进行探测并修正。另外，当站点位于地震活跃地带时，由于受地震运动的影响，当offset无法准确修正的情况下，并不适宜将整个跨度数据作为整体进行最优噪声模型分析并获取站速度及其不确定度等相关参数，故对其进行分段处理。

最后与传统方法不同的是，在进行噪声模型估计之前一般采用数学模型对坐标序列进行去趋势项(trend)，考虑到Trend项对噪声模型会产生一定的影响，且时间相关的噪声也可能呈现出趋势变化，因此在去trend的过程中，本发明专利只对已知的过程进行trend去除(如冰川均衡调整GI、板块旋转引起的趋势变化)，以保证后续结果分析的可靠性；

步骤3：对步骤2获取的GPS测站坐标位置时间序列继续进行环境负载修正，包括大气、地表水、非潮汐海洋、积雪负载改正，以改正地表环境负载引起的测站非线性变化；

步骤4：对步骤3获取的GPS测站坐标位置时间序列继续进行共模误差改正，采用PCA方法去除GPS站坐标序列中的共模误差ε_i(t_i)；

为了描述PCA方法在GPS时间序列(为去均值、趋势、负载改正后仅留周年、半周年项的残差序列)中的具体实现，假定GPS台站获得的三维坐标观测值形成一个n×m(n>m，n为观测数或历元数，m为观测类型)的数据矩阵X，其协方差阵为C_X，则C_X＝X^TX。数据矩阵如下：

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

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

X＝UΣV^T

C_X＝VΛV^T

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

a_k(t_i)可由下式求出：

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

经过步骤1～4，已获得“干净”的GPS测站坐标位置时间序列，考虑到GPS噪声模型的多样性，并结合噪声模型覆盖整个频率段的基本原则，选择单一噪声模型、组合噪声模型进行分析、解算；

其中备选的噪声模型包括：白噪声(White)、闪烁噪声(Flicker)、随机游走噪声(RandomWalk)、幂律噪声(Powerlaw)、高斯马尔科夫噪声(Generalized Gauss Markov，GGM)，ARMA(1)即FOGM、First order Gauss-Markov(FOGM)及其组合模型，充分考虑噪声模型的多样性。

步骤5：采用极大似然估计方法对不同组合的噪声模型分别对GPS测站坐标位置时间序列三坐标分量进行噪声估计，获得不同组合噪声模型在GPS测站坐标位置时间序列三坐标分量的MLE值；

步骤6：根据单一噪声模型、组合噪声模型GPS测站坐标位置时间序列三坐标分量的MLE值进一步求解GPS测站坐标位置时间序列坐标分量的最优噪声模型，其步骤如下：

步骤6.1：采用极大似然估计方法解算得到不同噪声模型的MLE值，并获取不同噪声模型(单一模型或组合模型)的待估参数个数k；

步骤6.2：采用Akaike Information Criteria(AIC)、Bayesian InformationCriteria(BIC)模型估计准则分别计算不同噪声模型的AIC、BIC值，AIC/BIC的求解原理如下：

AIC＝-2lnL+2k；

BIC＝-2·lnL+k·ln(n)；

式中L为某一噪声模型下的似然函数值MLE，n为观测值个数，k为模型的待估参数个数；根据AIC/BIC准则，AIC/BIC值越小，对应的模型越好(或者说趋近于真实模型)；

步骤6.3：当ACI、BIC估计出的GPS测站坐标位置时间序列噪声模型结果一致时，接受改估计结果；当ACI、BIC估计出的结果不一致时，采用功率谱分析的方法，通过频谱图的拟合趋近程度，对最优模型进行进一步确定。

请见图2，是本发明实施例的噪声模型估计结果图；经过本发明专利所描述方法进行数据处理后，通过AIC/BIC结合频谱图进行比较最终噪声模型的确定。图中给出了4个测站AIC、BIC估计结果不一致时处理方法，对于OHI2、MAS1站，FN+RW+WN模型对原始数据的拟合效果更佳，因此选该模型；同样对于CHUR站，FN+WN模型拟合效果更佳，对于KERG站，FN+WN、FN+RW+WN结果基本一致，对AIC、BIC估计出的RW分量进行分析，结果表明RW的比重小于0.1％，FN、WN模型的振幅、比重差异非常小，即AIC、BIC模型均可。

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

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

Claims (6)

1.一种GPS位置时间序列噪声模型建立方法，其特征在于，包括以下步骤：

步骤1：对原始观测文件进行解算，对不同松弛解进行加权解算，获取GPS测站坐标位置时间序列；

步骤2：对获取的GPS测站坐标位置时间序列进行数据预处理；

步骤3：对步骤2中预处理后的GPS测站坐标位置时间序列继续进行地表环境负载修正，以改正地表环境负载引起的GPS测站非线性变化；

步骤4：对步骤3中修正后的GPS测站坐标位置时间序列继续进行共模误差改正，去除坐标序列中的共模误差，获得残差序列；

步骤5：采用极大似然估计方法对不同组合的单一噪声模型、组合噪声模型分别对GPS测站坐标位置时间序列三坐标分量进行噪声估计，获得不同组合噪声模型在GPS测站坐标位置时间序列三坐标分量的MLE值；

步骤6：根据单一噪声模型、组合噪声模型GPS测站坐标位置时间序列三坐标分量的MLE值获得GPS测站坐标位置时间序列坐标分量的最优噪声模型，并获得站速度及其不确定性估计值。

2.根据权利要求1所述的GPS位置时间序列噪声模型建立方法，其特征在于：步骤2中所述预处理包括粗差分析、去均值、去趋势项，探测Offset并改正。

3.根据权利要求1所述的GPS位置时间序列噪声模型建立方法，其特征在于：步骤3中所述负载修正包括大气、地表水、非潮汐海洋、积雪负载改正。

4.根据权利要求1所述的GPS位置时间序列噪声模型建立方法，其特征在于：步骤4中所述共模误差改正，其具体实现过程是假定GPS台站获得的三维坐标观测值形成一个n×m的数据矩阵X，其协方差阵为C_X，其中n>m，n为观测数或历元数，m为观测类型；则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};$

假定：

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

5.根据权利要求1所述的GPS位置时间序列噪声模型建立方法，其特征在于，步骤6的具体实现包括以下子步骤：

步骤6.1：采用极大似然估计方法解算得到不同噪声模型的MLE值，并获根据不同噪声模型的待估参数个数k，进一步采用AIC、BIC结合功率谱分析的方法对最优噪声模型进行判定，并获得站速度及其不确定性的估值；

步骤6.2：采用AIC、BIC模型估计准则分别计算不同噪声模型的AIC、BIC值，AIC/BIC的求解原理如下：

AIC＝-2lnL+2k；

BIC＝-2·lnL+k·ln(n)；

式中L为某一噪声模型下的似然函数值MLE，n为观测值个数；

步骤6.3：当ACI、BIC估计出的GPS测站坐标位置时间序列噪声模型结果一致时，接受该估计结果；当ACI、BIC估计出的结果不一致时，采用功率谱分析的方法结合极大似然估计的结果，即通过功率谱分析方法分别对AIC、BIC估计后的结果以及原始数据进行频谱分析，并绘制频谱图，通过对不同方法估计出的频谱图与原始观测数据的频谱图的拟合趋近程度，进一步确定最优模型。

6.根据权利要求1-5任意一项所述的GPS位置时间序列噪声模型建立方法，其特征在于：步骤4中所述共模误差改正，是对于大尺度下的共模误差分离，以地表环境负载为参考，对大区域进行子网划分；并顾及GPS坐标序列、共模误差的周期性，在采用主成分分析法去除共模误差之前，保留坐标序列的周年、半周年项，以避免周期性共模分量的残留。