一种GNSS位置时间序列周期特性挖掘方法

摘要

本发明公开了一种GNSS位置时间序列周期特性挖掘方法，本发明考虑了传统的谐函数存在的局限性，引入小波分析、FAMOUS时频分析方法、极大似然估计等方法，提供了一种准确的GNSS位置时间序列周期性信号分析及特性分析方法。该方法依据原始序列中不同信号的固有频率不相同，从而将原始时间序列分解成不同频段的新的时间序列，将信号进行分解和重构，对每一层进行处理，从而获得所需要的部分。本发明克服了传统GNSS时间序列模型的局限性，有助于真实反映坐标序列的周期性变化，进一步提高GNSS站坐标的精度与可靠性，获得高精度的位置与速度参数，为深入了解相关地球物理现象的影响机制和变化规律提供依据。

说明书

技术领域

本发明属于连续运行全球定位导航系统技术领域，涉及一种GNSS位置时间序列周期特性挖掘方法。

背景技术

近来的研究表明，GNSS坐标时间序列呈现出明显的周期性变化(Dong et al.,2002；Van Dam et al.,2010,2012；Jiang et al.,2013)。

传统的GNSS坐标序列模型认为GNSS坐标序列仅包含周年、半周年项，忽视了其他周期性变化，对一些地球物理现象甚至可能做出错误的解释。如何准确、高效的识别GNSS坐标序列的周期特性，是GNSS坐标时间序列中的关键问题之一。

目前对GNSS位置时间序列周期特性挖掘一些缺陷：1)什么是周期性的物理起源，GNSS坐标序列中包含哪些周期性信号，缺乏深入的研究；2)目前大部分的周期特性挖掘研究以经典的谐函数为基础，通过最小二乘方法进行估计，存在一定的局限性；3)大部分研究没有对坐标序列的周期特性进行深入分析，仅仅确定获得的周期新信号的周期(主频)，没有对其物理属性、影响因素、起源进行探讨。

发明内容

为了解决上述技术问题，本发明提供了一种GNSS位置时间序列周期特性挖掘方法。

本发明所采用的技术方案是：一种GNSS位置时间序列周期特性挖掘方法，其特征在于，包括以下步骤：

步骤1：针对原始GNSS观测值，解算获取GNSS测站单日松弛解；通过公共基站进行不同解加权进行联合解算，获得GNSS测站坐标时间序列及速度参数；

步骤2：对获取的GNSS测站坐标时间序列建立残差时间序列模型；

步骤3：对获取的GNSS测站坐标时间序列进行预处理；

步骤4：对步骤3获得的GNSS测站坐标时间序列进行地表负载引起的测站位移进行纠正；

步骤5：消除GNSS测站非构造运动引起的测站非线性变化及对应的周期性波动，对步骤4获得的GNSS测站坐标时间序列进行纠正，获得纠正后的测站坐标时间序列；

步骤6：确定小波子序列的主频率范围，采用小波分析的方法对纠正后的测站坐标时间序列进行周期特性分析；

步骤7：对步骤6中获得的各子序列进行周期特性分析，分析各子序列的频谱图，并根据频谱图确定各子序列ENU三坐标序列频谱图主峰值对应的频率，即确定子序列的主频，再与确定的主频与步骤6中确定的小波子序列主频率范围进行比对，最终确定各子序列的主峰的频率；

步骤8：对确定周期主频频段后的子序列进行处理，采用极大似然估计方法进一步确定子序列的噪声模型、周期性变化振幅，以获得GNSS坐标序列周期特性具体参数值。

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

本发明的创新之处在于，一方面，克服了传统GNSS利用时间序列谐函数模型进行周期性信号估计的局限性，并对GNSS坐标序列进行粗差、阶跃、地表负载及共模误差纠正，以降低测站的非线性运动及随机误差的影响；另一方面本发明对引入小波分析、时频分析、极大似然估计等方法，提供了一种时频分析的周期性信号分析方法。该方法依据原始序列中不同信号的固有频率不相同，从而将原始时间序列分解成不同频段的新的时间序列，将信号进行分解和重构，对每一层进行处理，从而获得所需要的部分；此外，本发明专利奖噪声模型估计方法引入到周期特性分析中，对经过小波分析分解后的信号采用频率分析结合极大似然估计，更加准确的识别各子序列的周期特性，包括主频、周年信号振幅、噪声模型等，有助于真实反映坐标序列的周期性变化，进一步提高GNSS站坐标的精度与可靠性，获得高精度的位置与速度参数，为深入了解相关地球物理现象的影响机制和变化规律提供依据。

附图说明

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

图2是本发明实施例的小波分解原理示意图。

具体实施方式

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

请见图1，本发明提供的一种GNSS位置时间序列周期特性挖掘方法，包括以下步骤：

步骤1：针对GNSS观测值及相关文件(星历文件、表文件等)，本实施例采用高精度GNSS数据后处理软件以及相应的解算模型解算，分别获取GNSS测站单日松弛解，通过公共基站进行不同解加权进行联合解算，获得GNSS测站坐标时间序列；

其中步骤1中GNSS测站坐标时间序列解算方法主要包括2种：1)双差观测值进行处理，根据观测值之间的相关性，可以消除或减弱GNSS观测手段本身的误差(如卫星接收机钟差、电离层延迟等)，从而达到提高精度的目的，以消除随机误差引入的周期性新号。基于双差观测值的解算软件主要有GAMIT。2)直接对载波非差观测量以单点定位(Precise Point Positioning,PPP)的方式进行处理分析。基于非差观测值的解算方法具有可用观测值多，不同测站的观测值不相关的优点，能更加准确反映测站的真实运动。基于非差观测值的解算软件主要有GIPSY、BERNESE、PANDA。3)对双差、非差的中间解(松弛解)进行联合解算，可以有效的消除(或减弱)不同软件及模型引入的随机及系统误差，如软件及其模型的不完善、模型存在的系统偏差等，提高坐标序列的精度。

其中步骤1中本发明专利采用GAMIT、GIPSY进行加权进行联合解算(采用SOPAC给出的GAMIT/GIPSY＝1:2.4经验权计算)，以获得最终的GNNS测站单日解位置时间序列(X、Y、Z以及协方差矩阵)。

其中步骤1解算过程中所提及解算模型见下表1。

表1

联合解解算策略采用的模型固体潮模型IERS 2003 convention海潮模型FES04 model(CMC)极潮模型IERS 2003 convention投影函数Global Mapping Function天线模型绝对IGS相位中心

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

步骤2：对获取的测站坐标时间序列建立以下残差时间序列模型：

其中：x(t)_E/N/U为时刻t对应的GNSS测站坐标观测值，包括ENU三坐标分量；x₀为测站初始位置，v_x为线性速度即趋势项，t代表时间；代表周期性项阶数，为跳变改正项，g_j表示跳变振幅，T_gj为跳变发生的时刻即历元，n_g表示跳变个数，j为跳变编号，这里假定发生偏移的时刻T_g已知，H为海维西特阶梯函数，在跳变前H值为0，跳变后H值为1；ε_x(t)为随机误差，O_x为粗差)。

步骤3：对获取的GNSS测站坐标时间原始序列继续进行预处理，包括去均值、粗差剔除(O_x)、阶跃改正去均值的目的是便于对数据进行分析，去粗差、去阶跃主要为了消除粗差对周期信号的干扰；

与传统方法不同的是，考虑到传统基于谐函数结合最小二乘拟合去除趋势项、周年、半周年项的局限性，在对GNSS时间序列进行周期性信号探测时，若事先扣除上述项，会使得探测出的周期性信号存在一定的偏差，因此本发明实施的新颖之处在顾及了传统的谐函数存在的局限性，依据原始序列中不同信号的固有频率不相同，从而将原始时间序列分解成不同频段的新的时间序列，将信号进行分解和重构，对每一层进行处理，从而获得所需要的部分，有助于真实反映坐标序列的周期性变化。

其中步骤3中粗差剔除采用四分位数间距法(interquartile range,IQR)，四分位数间距由P25、P50、P75将一组变量值等分为四部分，P25称下四分位数(Q1)，P75称上四分位数(Q3)，将P75与P25之差定义为四分位数间距(IQR)。分别计算a(Q1-3*IQR)、b(Q3+3*IQR)的值，原始序列中位于(a,b)区间之外的值，则为粗差，IQR方法具有较好的抗差性，能有效剔除GNSS坐标序列中存在的粗差。

对于步骤3中GNSS测站坐标时间序列中存在的阶跃，采用如下方法进行改正：1)对已知的Offset根据SOPAC发布的发生的时刻及影响进行改正，2)SOPAC未公布的offset采用Median Interannual Difference Adjusted for Skewness(MIDAS)method进行改正。

步骤4：对步骤3获得的GNSS测站坐标时间序列进行地表负载引起的测站位移进行纠正。

采用mload程序分别计算大气压、非潮汐海洋、积雪和土壤水负荷引起的台站位移，通过负载改正，提高坐标序列的精度，消除部分非构造信号。

在计算地表负载引起的测站位移中采用如下数据模型：大气负荷采用NCEP的全球表面大气压力数据，时间分辨率为6小时，空间分辨率为2.5°×2.5°；非海洋潮汐负荷采用ECCO海洋底压力模型，时间分辨率为12小时，空间分辨率为空间分辨率为1°×1°；积雪和土壤水负荷数据来自NCEP。

步骤5：进行共模误差剔除，以消除测站非构造运动引起的测站非线性变化及对应的周期性波动。

共模误差分离采用PCA方法进行剔除，其具体实现过程为：

假定GNSS台站获得的三维坐标观测值形成一个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个主成分(k＝1,2,3)，v_k(x)是对应主成分的响应特征矩阵，分别代表时间特征和空间响应，取前k个主分量(k≤4，累积贡献率大于80％)计算得到的共模误差为：

对原始序列进行纠正，获得纠正后坐标序列(S)。

步骤6：经步骤5后，获得纠正后的坐标序列(S)，对采用小波分析的方法其进行周期特性分析。首先确定小波子序列的主频率范围(附表2)，然后依据原始序列中不同信号的固有频率不相同，从而将原始时间序列进行小波分解(请见附图2)成不同频段的新的时间序列，将信号进行分解和重构，对每一层进行处理，从而获得各子序列(D₁，D₂，…D_n，A₁，A₂，…A_n)。；

表2

步骤7：对步骤6获得的各子序列进行租期特性分析，采用Frequency Analysis Mapping On Unusual Sampling(FAMOUS)方法分析各子序列的频谱图，确定其频谱主峰值对应的频率，确定各子序列的主频，并与步骤6中确定的主频范围进行比对，最终确定各子序列的主峰的频率。

步骤8：对确定周期属性后的子序列处理，采用极大似然估计方法对确定周期主频频段后的子序列进行处理，进一步确定其周期特性(噪声模型、周期性变化振幅等),，以获得GNSS坐标序列周期特性具体参数值。

针对现有技术存在的问题，本发明考虑到已有方法的局限性，进一步对GNSS位置时间序列周期特性进行探讨，考虑了传统的基于谐函数GNSS位置时间序列模型存在的局限性，引入小波分析、时频分析、极大似然估计等方法，提供了一种GNSS位置时间序列周期特性挖掘方法，高效、准确的分离出GNSS位置时间序列周期特性，有助于真实反映坐标序列的周期性变化，进一步提高GNSS站坐标的精度与可靠性，获得高精度的位置与速度参数，为深入了解相关地球物理现象的影响机制和变化规律提供依据。

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

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