GNSS测站非线性运动建模方法与装置与流程

本发明涉及一种gnss测站非线性运动建模方法与装置，属于全球定位导航系统技术领域。

背景技术：

随着gnss(globalnavigationsatellitesystem，全球卫星导航系统)观测数据的不断丰富和数据处理技术的不断发展，分析测站非线性运动规律可以获得各种地球物理现象和季节性影响，如果对这部分残差规律进行建模拟合从而消除其影响，那么可以进一步提高gps台站坐标的精度。

尽管目前国际地球参考框架itrf(internationalterrestrialreferenceframe)中基准站的历元坐标和速度场已达到了毫米级，但几乎所有igs(internationalgpsservice)站坐标时间序列(尤其是垂直方向)都呈现显著的非线性运动趋势，且振幅可达1～2cm。因此，地面点相对于参考框架的运动无法用线性速度进行完整描述，仅基于线性速率的模型来维持参考框架会有一定的局限性，而且其精度只能达到厘米级。一般情况下，要精确表达站点坐标，通常采用两种方式描述非线性运动：①不考虑引起测站坐标非线性变化的各种物理机制，只根据坐标时间序列本身的运动趋势建模；②从产生形变的物理机理入手，分析各因素影响从而对测站坐标进行改正。目前，国内外对igs站坐标时间序列的分析已有不少相关研究，也取得较大进展。张鹏，蒋志浩等在2007年武汉大学学报(信息科学版)第32卷第3期第251-254页，公开了一篇名称为《我国gps跟踪站数据处理与时间序列特征分析》的文献，并在该文献中对我国igs站坐标时间序列进行频谱分析，发现高程方向周年特性表现明显；傅彦博，孙付平等在2018年测绘学报第47卷第10期第1337-1345页，公开了一篇名称为《全球gps测站垂向周年变化统计改正模型的建立》的文献，并在该文献中对全球461个gps测站垂向时间序列所包含的各类周期项进行了统计，发现周年项在全球具有普遍性，其次是半周年项。另有学者证实了大气负荷、水文负载、非潮汐海洋负载等负载变化是引起gps测站垂向位移的主要因素。

奇异谱分析(singularspectrumanalysis，ssa)是时间序列分析中一种常用的方法，在气候学、测量学和海洋学等领域应用广泛，它的优越性在于无需任何先验信息和假设条件，能稳定识别和强化周期信号，把最可预报的分量聚集到若干个时间序列中，以此选择若干有意义的分量重建序列，降低噪声影响，因此特别适合于分析有周期震荡的时间序列数据。然而，只根据贡献率较大的前几项来对时间序列进行重构，仅能有效地将半年及半年以上的周期项提取出来，而短周期项如季节周期项、月周期项则容易当作噪声被忽略掉。因此，需要一定方法削弱这种影响。

技术实现要素：

本发明的目的是提供一种gnss测站非线性运动建模方法，用以解决现有技术会忽略短周期项导致建模精度低的问题；本发明还提供一种gnss测站非线性运动建模装置，用以解决现有技术会忽略短周期项导致建模精度低的问题。

为实现上述目的，本发明提供了一种gnss测站非线性运动建模方法，包括以下步骤：

获取gnss测站的原始坐标时间序列；

对原始坐标时间序列进行小波分解和重构，得到低频部分和多层的高频部分；

对所述低频部分和每层高频部分分别进行奇异谱分解和重建；

将奇异谱分解和重建后的低频部分和每层高频部分进行合成得到原始坐标时间序列的拟合序列；

对所述拟合序列进行建模。

该方法的有益效果是：该方法先利用小波分解和重构将原始坐标时间序列分解成频率成分上相对单一、平滑的低频部分和各层高频部分，接着针对低频部分和各层高频部分分别进行奇异谱分解和重建，然后将奇异谱分解和重建后的各部分进行合成得到原始坐标时间序列的拟合序列，最后对合成后的拟合序列进行建模，在一定程度上降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率，从而提高了建模精度。

为了实现小波分解和重构，进一步地，采用haar小波、dbn小波、symn小波或bior小波模型对原始坐标时间序列进行小波分解和重构。

为了实现奇异谱分解和重建，进一步地，进行奇异谱分解和重建时选择的窗口长度为365。

本发明还提供了一种gnss测站非线性运动建模装置，包括存储器和处理器，所述处理器用于运行存储在所述存储器中的程序指令，以实现如下方法：

获取gnss测站的原始坐标时间序列；

对原始坐标时间序列进行小波分解和重构，得到低频部分和多层的高频部分；

对所述低频部分和每层高频部分分别进行奇异谱分解和重建；

将奇异谱分解和重建后的低频部分和每层高频部分进行合成得到原始坐标时间序列的拟合序列；

对所述拟合序列进行建模。

该装置的有益效果是：该装置先利用小波分解和重构将原始坐标时间序列分解成频率成分上相对单一、平滑的低频部分和各层高频部分，接着针对低频部分和各层高频部分分别进行奇异谱分解和重建，然后将奇异谱分解和重建后的各部分进行合成得到原始坐标时间序列的拟合序列，最后对合成后的拟合序列进行建模，在一定程度上降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率，从而提高了建模精度。

为了实现小波分解和重构，进一步地，采用haar小波、dbn小波、symn小波或bior小波模型对原始坐标时间序列进行小波分解和重构。

为了实现奇异谱分解和重建，进一步地，进行奇异谱分解和重建时选择的窗口长度为365。

附图说明

图1是本发明的gnss测站非线性运动建模方法原理图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施例对本发明进行进一步详细说明。

方法实施例：

如图1所示，本实施例的gnss测站非线性运动建模方法，包括以下步骤：

步骤1、获取gnss测站的原始坐标时间序列s。具体过程如下：

选取igs网站(该网站的网址为ftp://igs-rf.ensg.eu/pub/res/)上全球范围内某个纬度的gps站坐标时间序列数据，作为gnss测站的原始坐标时间序列s，设原始坐标时间序列s＝{s1,s2,s3…}。

步骤2、对原始坐标时间序列s进行小波分解和重构，得到低频部分和多层的高频部分。具体过程如下：

本实施例中，采用haar小波模型对原始坐标时间序列s＝{s1,s2,s3…}进行小波多尺度分解和重构，得到下式：

s＝ab+d1+d2+d3+…+dn(1)

式中，ab＝{ab1,ab2,ab3…}为原始坐标时间序列s中低频信号的重构结果(即低频部分)，d1＝{d11,d12,d13…},…,di＝{di1,di2,di3…},…,dn＝{dn1,dn2,dn3…}分别为原始坐标时间序列s中第1层到第n层的高频信号的重构结果(即多层的高频部分)。作为其他实施方式，还可以采用dbn小波、symn小波或bior小波模型对原始坐标时间序列s进行小波多尺度分解和重构。

那么，ti时刻的原始坐标时间序列s可以用下式表示：

si＝abi+d1i+d2i+…+dni(2)

式中，abi为低频信号在ti时刻的值，d1i,d2i,…,dni分别为各层高频信号在ti时刻的值，i表示变量。

步骤3、对低频部分和每层高频部分分别进行奇异谱分解和重建(即ssa建模)。

首先，利用ssa方法(即奇异谱分析方法)对低频部分ab＝{ab1,ab2,ab3…}进行奇异谱分解和重建，得到低频信号在ti时刻的拟合值为然后，利用ssa方法对每层高频部分di＝{di1,di2,di3…},i＝1,2,3,...,n分别进行奇异谱分解和重建，得到各层高频信号在ti时刻的拟合值分别为本实施例中，进行奇异谱分解和重建时，低频部分和每层高频部分对应的窗口长度、重构阶次是相同的。

下面以一维时间序列x＝(x1,x2,x3,...,xn)为例，详细介绍利用ssa方法对其进行奇异谱分解和重建的过程，对低频部分和每层高频部分分别进行奇异谱分解和重建的过程与此类似，主要由以下4个步骤构成：

第一步，构建时滞矩阵。

首先，选择适当的窗口长度l，满足一般情况下，l的选取不宜超过整个时间序列数据长度n的1/3。其中，窗口长度l并没有统一的选择标准，而是需要根据具体应用而定，一般说来，l过大，会导致奇异值分解得到不同的部分产生混叠，而l太小，则无法使用信号从弱到强渐进划分，从而使得部分信号无法获取。如果能根据经验大致确定时间序列中数据的周期特征，则l一般取周期的公倍数。本实施例中，由于选取的原始坐标时间序列s中存在的已知周期有周年周期和半年周期，且本实施例所用数据的时间分辨率为一天，因此，可以选择已知周期的最小公倍数365作为窗口长度l。

其次，根据一维时间序列x＝(x1,x2,x3,...,xn)构建时滞矩阵x，得到：

式(3)中，k＝n-l+1，x为l×k阶时滞矩阵，且矩阵x的每一条副对角线上的元素相等，即若将矩阵x中第i行、第j列的元素用xi,j表示，则有xi,j＝xi-1,j+1，因此矩阵x是一个hankel矩阵(即汉克尔矩阵)。

第二步，svd分解(即奇异值分解)。

对时滞矩阵x进行奇异值分解：

式(4)中，由矩阵xx^t和矩阵x^tx的前d(d＝min{l,k})个最大非零特征值构成对角矩阵为矩阵x的奇异值，等价于矩阵xx^t特征值的平方根；u为矩阵x的左奇异向量，等价于矩阵xx^t的特征向量；v为矩阵x的右奇异向量，等价于x^tx的特征向量。

由此推知：

式中，为矩阵x的奇异值，为奇异谱；ui通常由经验正交函数(empiricalorthogonalfunctions，eof)表示；vi为主成分(principlecomponents，pc)；为矩阵x的第i个三重特征向量。最大特征值对应的特征向量，代表了信号的最大变化趋势，而较小的特征值对应的特征向量一般被当作噪声。

因此，时滞矩阵x的奇异值分解可以表示为：

x＝x1+x2+…+xd(7)

第三步，分组。

将矩阵xi的下标{1,2,…,d}分成m个互不相交的集合i1,i2,…,im，设i＝[i1,i2,…,ip，那么与集合i相关的矩阵xi就可以表示为xi＝xi1+xi2+…+xip，时滞矩阵x就可以表示成：

x＝xi1+xi2+…+xim(8)

一般说来，在svd分解中，前r(r＜d)个特征向量对时滞矩阵x的贡献率ηr可以表示为：

第四步，对角平均化。

对角平均化的目的就是把第三步分解得到的矩阵xim重新转换成为长度为n的新时间序列，称之为重建成分(reconstructioncomponent，rc)，所有rc之和等于原始序列x。

假设z＝z1,z2,…,zn为z经过对角平均化得到的时间序列，则对角平均化的公式可表示为：

这样，所有重建成分rc叠加后的和与原始序列x相同，即：

截取前p个贡献较大的成分近似表示原始序列x，则有：

其中，重构阶次p的选择需要根据具体应用而定，一般来说，重构阶次p主要根据奇异值的贡献率来确定，若p太小，则表示后面部分信号将被当作噪声剔除掉，若p太大，则会使得部分噪声被当作信号而提取出来。实际应用时，可以将贡献率较大的前几项之和与设定的某一阈值进行比较，来确定重构阶次p。

本实施例中，由于选取的原始坐标时间序列s中的周期特征主要集中在前几个特征值重构的序列中，所以重构阶次p一般取6～10即可。

步骤4、将奇异谱分解和重建后的低频部分和每层高频部分进行合成得到原始坐标时间序列s的拟合序列

根据步骤3中得到的原始坐标时间序列s中低频信号和各层高频信号在ti时刻的拟合值，得到原始坐标时间序列s在ti时刻的拟合序列为：

式(13)中，为低频信号在ti时刻的拟合值，分别为各层高频信号在ti时刻的拟合值。

步骤5、对拟合序列进行建模。

利用步骤4得到原始坐标时间序列的拟合序列后，可以根据实际需要建立多种模型，例如：(1)利用拟合序列建立时间序列趋势模型，用于分离和提取拟合序列中的趋势成分，进而得到gnss测站三维方向存在的上升或移动趋势的量级；(2)利用拟合序列建立多路径效应改正模型，用于削弱原始坐标时间序列中多路径效应的影响；(3)利用拟合序列建立最优噪声模型，用于估计gnss测站速度场的运动趋势；(4)利用拟合序列建立共模误差估计模型，用于提取原始坐标时间序列中的共模误差并剔除，进而提高原始坐标时间序列的精度。

基于奇异谱分析的非线性运动建模方法能有效地将半年及半年以上的周期项提取出来，而短周期项如季节周期项、月周期项则容易被忽略。本发明有效融合了小波多尺度分解和奇异谱分析两种方法的优点，先利用小波分解和重构将原始坐标时间序列分解成频率成分上相对单一、平滑的低频部分和各层高频部分，接着针对低频部分和各层高频部分分别进行奇异谱分解和重建，然后将奇异谱分解和重建后的各部分进行合成得到原始坐标时间序列的拟合序列，最后对合成后的拟合序列进行建模，在一定程度上降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率，从而提高了建模精度。

装置实施例：

本实施例的gnss测站非线性运动建模装置，包括存储器和处理器，处理器用于运行存储在存储器中的程序指令，以实现如下方法：

获取gnss测站的原始坐标时间序列；

对原始坐标时间序列进行小波分解和重构，得到低频部分和多层的高频部分；

对低频部分和每层高频部分分别进行奇异谱分解和重建；

将奇异谱分解和重建后的低频部分和每层高频部分进行合成得到原始坐标时间序列的拟合序列；

对拟合序列进行建模。

其中，该方法的具体实施过程与方法实施例中gnss测站非线性运动建模方法的实施过程相同，此处不再赘述。