CN110632625B - 一种gnss时间序列阶跃探测与修复方法 - Google Patents

Abstract

本发明公开了一种GNSS时间序列阶跃探测与修复方法，针对GNSS时间序列阶跃问题，利用时间序列发生阶跃前后连续时间窗口内的趋势不变特性，通过移动窗口法和时间基函数保系数趋势项逼近方式构建了阶跃探测与修复模型，进而采用稳健最小二乘、广义最小二乘或其他最优估计理论进行阶跃参数估计和方差‑协方差估计，并通过统计假设检验进行阶跃参数显著性检验与阶跃修复。本发明适用于非等时间间隔或数据缺失下多种类型的GNSS时间序列探测与修复，阶跃探测与修复估计理论严密、易于编程，探测结果具有最优统计性质、修复结果具有统计可靠性。

Description

一种GNSS时间序列阶跃探测与修复方法

技术领域

本发明属于测绘科学与技术、GNSS测量数据处理与应用领域，具体涉及一种GNSS时间序列阶跃探测与修复方法。

背景技术

随着空间大地测量技术的发展，特别是北斗三代导航卫星系统的相继发射和组网，全球卫星导航定位系统(GNSS)广泛地应用于导航定位、变形监测、授时等领域。但是GNSS钟差时间序列、载波相位时间序列和GNSS站坐标时间序列(统称为GNSS时间序列)因各种原因相应地存在钟跳、周跳和阶跃(统称为阶跃项)问题，限制了GNSS应用的精度和可靠性。

目前国内外针对GNSS时间序列阶跃问题开展了一系列研究。如专利号为CN201711030935.4的中国专利“一种基于相位跳变的卫星钟差实时预报方法”利用钟差相位跳变的特性，基于多项式的拟合精度进行跳变的探测，并采用修复后的序列进行多项式钟差预报。专利号为CN201611108578.4的中国专利“一种顾及GNSS接收机钟跳的导出多普勒观测值构造方法”以钟差一次差为基础，以0.5毫秒为阈值进行钟跳的检验。专利号为CN201710469425.0的中国专利“一种GPS载波相位周跳探测与修复的方法”基于MW组合的差分值利用小波奇异值分解和自适应指数平滑进行载波相位观测值的周跳探测和修复。专利号为CN201710997434.7的中国专利“一种INS辅助的实时BDS单频周跳探测方法”利用真实双差相位观测值和INS信息计算双差相位预报值进行周跳探测。专利号为CN201710297705.8的中国专利“基于共模误差的GPS坐标时间序列不连续性的补足方法”提出了基于共模误差的数据插值方法。专利号为CN201710035054.5的中国专利“一种GNSS位置时间序列周期特性挖掘方法”对坐标时间序列中的周期项进行了探测。

综上，对于GNSS钟差时间序列，一般以时间序列一阶差分为基础采用固定阈值法进行钟跳探测，一阶差分无疑放大了随机误差而难以探测小阶跃；对于GNSS载波相位观测值中的小周跳，一般采用信号分解方法或其他辅助信息进行阶跃项探测，计算较为复杂、实时性差而不利于推广使用；GNSS站坐标时间序列的研究则鲜有针对阶跃项的探测和修复，即使进行探测和修复，大多采用目视确定和手动修复GNSS时间序列中的阶跃项，这种方法效率较低，且无法探测和修复较小的阶跃项。此外，现有专利技术重探测轻修复，而且未见到针对上述GNSS时间序列共性问题提出阶跃项探测与修复技术方案。

发明内容

发明目的：针对以上问题，本发明提出一种GNSS时间序列阶跃探测与修复方法，该方法简单可行、适用于不同类型的GNSS时间序列阶跃探测和修复，解决全球卫星导航定位系统(GNSS)时间序列阶跃探测与修复问题。

技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种GNSS时间序列阶跃探测与修复方法，基于GNSS时间序列在发生阶跃前、后连续时间窗口内的趋势不变特性，通过时间基函数保系数趋势项逼近方式对包含拟阶跃节点的移动窗口时间序列建立阶跃探测与修复模型，利用稳健最小二乘、广义最小二乘进行阶跃参数估计，并基于正态分布检验进行阶跃参数显著性判断，进而进行GNSS时间序列阶跃项修复；本发明方法包括以下步骤：

(1)获取GNSS基准站的时序文件，读取原始观测值；所述原始观测值指的是时序文件天顶方向的GNSS站坐标时间序列，GNSS组合钟差时间序列或GNSS组合相位观测时间序列；

(2)根据GNSS时间序列类别和移动窗口法，选取阶跃前、后连续时间窗口的移动窗口节点长度，确定时间基函数及其最高阶数；初始化移动窗口序号i＝1；

(3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内的GNSS时间序列进行时间基函数保系数趋势项逼近；

(4)构建基于最优估计理论的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；所述最优估计理论包括：稳健最小二乘估计、广义最小二乘估计或其他最优估计理论；

(5)根据参数估计及其方差-协方差估计，构建正态分布统计量，进行GNSS时间序列阶跃参数的显著性检验，若通过显著性检验，则修复GNSS时间序列阶跃项，并令修复后的时间序列替换原时间序列；若未通过显著性检验，则表示GNSS时间序列不存在阶跃项，无需修复；

(6)将移动窗口序号i增加为i+1，重复上述步骤(3)～步骤(5)，直至完成整个时间序列的阶跃探测和修复。

进一步，步骤(2)所述选取阶跃前、后连续时间窗口的移动窗口节点长度，确定时间基函数及其最高阶数；方法如下：

(2-1)设第i个移动窗口内的时间节点为t∈[t_i+1,t_i+2K]、窗口节点长度为2K，t_i+1为移动窗口起始时序节点，t_i+2K为移动窗口终止时序节点，相应的GNSS时间序列为y(t)；

(2-2)根据GNSS时间序列类别确定时间基函数及其最高阶数，时间基函数表示为{F_j(·)}_1≤j≤m，其中，F_j(·)表示j阶时间基函数，m表示时间基函数的最高阶数。

进一步，步骤(3)所述基于GNSS时间序列的趋势不变特性，将第i个移动窗口内的GNSS时间序列进行时间基函数保系数趋势项逼近；方法如下：

(3-1)将第i个移动时间窗口中的时间节点分为长度相等的两部分，前半部分和后半部分时间节点表示如下：

{t′_i,k＝t_i+k}_1≤k≤K

{t″_i,k＝t_i+K+k}_1≤k≤K

式中，t′_i,k、t″_i,k分别为移动窗口内前、后两个部分的时间节点；

(3-2)为防止出现时间节点值过大而导致数值计算不适定问题，分别对移动窗口内前、后两个部分的时间节点进行标准化处理；对第i个移动窗口内前、后两个部分的时间节点进行标准化的方法为：

τ′_i,k＝t′_i,k-t″_i,1 (1)

τ″_i,k＝t″_i,k-t″_i,1 (2)

式中，τ′_i,k、τ″_i,k分别为t′_i,k、t″_i,k的标准化时间节点；t″_i,1＝t_i+K+1表示拟发生阶跃时间节点；

(3-3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内前、后两个部分的GNSS时间序列进行保系数趋势项逼近，表达形式如下：

式中，y(τ′)表示第i个移动时间窗口内的前半部分时间序列；y(τ″)表示第i个移动时间窗口内的后半部分时间序列；F_j(·)表示j阶时间基函数，m表示时间基函数的最高阶数，τ′和τ″分别为移动窗口内前、后两部分时间序列的标准化时间节点；β_j表示j阶时间基函数的系数；ε′和ε″分别为移动窗口内前、后两部分时间序列的观测噪声，Step表示移动窗口内前、后两部分时间序列趋势项的阶跃参数。

进一步，步骤(4)所述最优估计理论选用稳健最小二乘估计；构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；方法如下：

(4-1-1)构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型；

设观测噪声服从偶然误差特性，根据时间基函数保系数趋势项逼近分析得到发生阶跃前、后两部分时间序列具有相同的拟合系数，建立GNSS时间序列阶跃探测与修复模型为：

式中，B_i为系数矩阵，L_i为常数向量，V_i表示残差向量，d_i＝[Step,β₀,…,β_m]^T表示第i个移动窗口的待估参数向量；

所述系数矩阵B_i，常数向量L_i，残差向量V_i表示如下：

(4-1-2)利用稳健最小二乘准则求解第i个移动窗口中的待估参数

所述稳健最小二乘准则表示为：

待估参数

计算公式为：

式中，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵；

(4-1-3)将参数估计结果

回代模型，即将

代入式(4)求得残差向量V_i，进而求得单位权方差

(4-1-4)根据方差-协方差传播定律计算模型参数估值

的方差-协方差阵D_i；所述方差-协方差阵计算公式如下：

式中，

为单位权方差，B_i为系数矩阵，

为第i个移动窗口内GNSS时间序列的稳健权矩阵。

进一步，步骤(4)所述最优估计理论选用广义最小二乘估计；构建基于广义最小二乘估计的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；方法如下：

(4-2-1)构建广义最小二乘的GNSS时间序列阶跃探测与修复模型；

设移动窗口节点长度为2K的第i个移动窗口的待估参数向量为：

d_i＝[Step,β₀,…,β_m]^T

将参数向量d_i进行扩展，重构待估参数向量为：

根据重构后的待估参数向量

建立GNSS时间序列阶跃探测与修复模型为：

式中，

表示系数矩阵，L_i表示常数向量，C表示约束矩阵，

表示残差向量，V_z表示约束误差向量；

所述系数矩阵

约束矩阵C的转置矩阵C^T表示如下：

式中，I_m+1表示m+1阶单位矩阵；

(4-2-2)利用广义最小二乘准则求解第i个移动窗口中的参数估值

所述广义最小二乘准则表示为：

待估参数

计算公式为：

式中，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵，P_z表示参数权矩阵；

(4-2-3)将阶跃参数估计结果

回代模型，即将

代入式(8)分别计算残差向量

和约束误差向量V_z，并估计单位权方差

(4-2-4)根据方差-协方差传播定律计算阶跃参数估值

的方差-协方差阵

所述方差-协方差阵计算公式如下：

式中，

为单位权方差，

为系数矩阵，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵，P_z表示参数权矩阵，C表示约束矩阵。

进一步，步骤(5)所述构建正态分布统计量，进行GNSS时间序列阶跃参数的显著性检验，若通过显著性检验，则修复GNSS时间序列阶跃，若未通过显著性检验，则无需修复GNSS时间序列；方法如下：

(5-1)根据参数估值及其方差-协方差估计，构建阶跃参数Step显著性检验的标准正态分布统计量U_i：

式中，Step表示阶跃参数，σ_Step表示阶跃中误差；当采用稳健最小二乘的时间序列阶跃探测与修复模型，

当采用广义最小二乘的时间序列阶跃探测与修复模型，

其中，

分别为向量

的第一个元素，D_i(1,1)、

分别为矩阵D_i、

的第一行第一列元素；

(5-2)根据预先选定的置信水平α，求得标准正态分布的接受域和拒绝域，所述接受域和拒绝域分别为式(13)和式(14)：

|U_i|≤u_α/2 (13)

|U_i|>u_α/2 (14)

式中，|·|表示取绝对值，u_α/2为正态分布的上α分位点；

(5-3)当正态分布统计量U_i落在接受域内，即满足式(13)时，表示阶跃参数Step不显著、未通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处没有发生阶跃，即无需修复；

反之，当正态分布统计量U_i落在拒绝域内，即满足式(14)时，表示阶跃参数Step显著、通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处发生阶跃，即需要修复；

(5-4)当所求阶跃参数通过显著性检验时，对时间节点t_i+K+1以后的时间序列进行修复，公式如下：

式中，y_Step(t)表示修复后的GNSS时间序列。

有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：

(1)本发明直接采用带时间标记的原始时间序列而非差分时序，适用于非等时间间隔或数据缺失情况下的GNSS时间序列阶跃探测与修复，有利于提高阶跃探测精度；

(2)本发明充分利用阶跃发生前、后连续时间窗口内的时序趋势不变特性，采用移动窗口法和时间基函数保系数趋势项进行逼近，从而实现了阶跃参数的直接分离和估计，适用于不同类型GNSS时间序列；

(3)本发明采用稳健最小二乘平差、广义最小二乘或者其他最优估计理论，并结合统计假设检验进行阶跃项修复，保证了GNSS时间序列阶跃探测的统计最优性质与阶跃修复的统计可靠性。

附图说明

图1是本发明方法的GNSS时序阶跃探测与修复流程图；

图2是本发明实施例的某GNSS基准站坐标时间序列探测与修复结果。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下文将结合附图对发明实施方式作进一步的描述。需说明的是，所提供附图仅仅是本发明的实施例示意图，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图，并获得其他的实施方式。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

选取中国大陆构造环境监测网络(CMONOC)中的某GNSS基准站2000.0年至2017.0年共17年天顶方向的GNSS站坐标时间序列，该坐标时间序列以ITRF2014为基准，时间分辨率为1天。采用本发明方法进行阶跃探测与修复，经过探测和修复后的坐标时间序列见图2，从该图可以看出在2000.59、2003.08、2007.68、2010.49共4处存在阶跃，探测修复结果与实际结果吻合。具体实施步骤如下：

(1)从CMONOC数据中心(www.cgps.ac.cn/cgs/index.action)下载GNSS基准站的时序文件，读取文件天顶方向的GNSS站坐标时间序列。

(2)根据GNSS时间序列类别和移动窗口法，选取移动窗口节点长度和时间基函数的最高阶数，初始化移动窗口序号i＝1。本实施例中，移动窗口节点长度取为60，时间基函数取一阶线性模型。

(3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内的GNSS时间序列进行时间基函数保系数趋势项逼近；具体如下：

(3-1)将第i个移动时间窗口中的时间节点分为长度相等的两部分，前半部分和后半部分时间节点表示如下：

{t′_i,k＝t_i+k}_1≤k≤K

{t″_i,k＝t_i+K+k}_1≤k≤K

式中，t′_i,k、t″_i,k分别为移动窗口内前、后两个部分的时间节点；

(3-2)分别对第i个移动窗口内前、后两个部分的时间节点进行标准化处理；标准化方法为：

τ′_i,k＝t′_i,k-t″_i,1 (1)

τ″_i,k＝t″_i,k-t″_i,1 (2)

式中，τ′_i,k、τ″_i,k分别为t′_i,k、t″_i,k的标准化时间节点；t″_i,1＝t_i+K+1表示拟发生阶跃时间节点；

(3-3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内前、后两个部分的GNSS时间序列进行保系数趋势项逼近，表达形式如下：

式中，y(τ′)表示第i个移动时间窗口内的前半部分时间序列；y(τ″)表示第i个移动时间窗口内的后半部分时间序列；F_j(·)表示j阶时间基函数，m表示时间基函数的最高阶数，τ′和τ″分别为移动窗口内前、后两部分时间序列的标准化时间节点；β_j表示j阶时间基函数的系数；ε′和ε″分别为移动窗口内前、后两部分时间序列的观测噪声，Step表示移动窗口内前、后两部分时间序列趋势项的阶跃参数。

(4)构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；具体如下：

(4-1-1)构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型；

设观测噪声服从偶然误差特性，根据时间基函数保系数趋势项逼近分析得到发生阶跃前、后两部分时间序列具有相同的拟合系数，建立GNSS时间序列阶跃探测与修复模型为：

式中，B_i为系数矩阵，L_i为常数向量，V_i表示残差向量，d_i＝[Step,β₀,…,β_m]^T表示第i个移动窗口的待估参数向量；

所述系数矩阵B_i，常数向量L_i，残差向量V_i表示如下：

(4-1-2)利用稳健最小二乘准则求解第i个移动窗口中的待估参数

所述稳健最小二乘准则表示为：

待估参数

计算公式为：

式中，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵；

(4-1-3)将参数估计结果

回代模型，即将

代入式(4)求得残差向量V_i，进而求得单位权方差

(4-1-4)根据方差-协方差传播定律计算模型参数估值

的方差-协方差阵D_i；所述方差-协方差阵计算公式如下：

式中，

为单位权方差，B_i为系数矩阵，

为第i个移动窗口内GNSS时间序列的稳健权矩阵。

(5)根据参数估计及其方差-协方差估计，构建正态分布统计量，进行GNSS时间序列阶跃参数的显著性检验，若通过显著性检验，则修复GNSS时间序列阶跃项，并令修复后的时间序列替换原时间序列；若未通过显著性检验，则表示GNSS时间序列不存在阶跃项，无需修复；具体如下：

(5-1)根据参数估值及其方差-协方差估计，构建阶跃参数Step显著性检验的标准正态分布统计量U_i：

式中，Step表示阶跃参数，σ_Step表示阶跃中误差；当采用稳健最小二乘的时间序列阶跃探测与修复模型，

当采用广义最小二乘的时间序列阶跃探测与修复模型，

其中，

分别为向量

的第一个元素，D_i(1,1)、

分别为矩阵D_i、

的第一行第一列元素；

(5-2)根据预先选定的置信水平α，求得标准正态分布的接受域和拒绝域，所述接受域和拒绝域分别为式(13)和式(14)：

|U_i|≤u_α/2 (13)

|U_i|>u_α/2 (14)

式中，|·|表示取绝对值，u_α/2为正态分布的上α分位点；

(5-3)当正态分布统计量U_i落在接受域内，即满足式(13)时，表示阶跃参数Step不显著、未通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处没有发生阶跃，即无需修复；

反之，当正态分布统计量U_i落在拒绝域内，即满足式(14)时，表示阶跃参数Step显著、通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处发生阶跃，即需要修复；

(5-4)当阶跃参数通过显著性检验时，对时间节点t_i+K+1以后的时间序列进行修复，公式如下：

式中，y_Step(t)表示修复后的GNSS时间序列。

(6)将移动窗口序号i增加为i+1，重复上述步骤(3)～步骤(5)，直至完成整个时间序列的阶跃探测和修复。

Claims (7)

1.一种GNSS时间序列阶跃探测与修复方法，其特征在于：该方法包括以下步骤：

(1)获取GNSS基准站的时序文件，读取原始观测值；所述原始观测值指的是时序文件天顶方向的GNSS站坐标时间序列、GNSS组合钟差时间序列或GNSS组合相位观测时间序列；

(2)根据GNSS时间序列类别和移动窗口法，选取阶跃前后连续时间窗口的移动窗口节点长度，确定时间基函数及其最高阶数；初始化移动窗口序号i＝1；

(3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内的GNSS时间序列进行时间基函数保系数趋势项逼近；

(4)构建基于最优估计理论的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；

(5)根据参数估计及其方差-协方差估计，构建正态分布统计量，进行GNSS时间序列阶跃参数的显著性检验，若通过显著性检验，则修复GNSS时间序列阶跃项，并令修复后的时间序列替换原时间序列；若未通过显著性检验，则表示GNSS时间序列不存在阶跃项，无需修复；

(6)将移动窗口序号i增加为i+1，重复上述步骤(3)～步骤(5)，直至完成整个时间序列的阶跃探测和修复。

2.根据权利要求1所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：步骤(2)所述选取阶跃前后连续时间窗口的移动窗口节点长度，确定时间基函数及其最高阶数；方法如下：

(2-1)设第i个移动窗口内的时间节点为t∈[t_i+1,t_i+2K]、窗口节点长度为2K，t_i+1为移动窗口起始时序节点，t_i+2K为移动窗口终止时序节点，相应的GNSS时间序列为y(t)；

(2-2)根据GNSS时间序列类别确定时间基函数及其最高阶数，时间基函数表示为{F_j(·)}_1≤j≤m，其中，F_j(·)表示j阶时间基函数，m表示时间基函数的最高阶数。

3.根据权利要求2所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：步骤(3)所述基于GNSS时间序列的趋势不变特性，将第i个移动窗口内的GNSS时间序列进行时间基函数保系数趋势项逼近；方法如下：

(3-1)将第i个移动时间窗口中的时间节点分为长度相等的两部分，前半部分和后半部分时间节点表示如下：

{t′_i,k＝t_i+k}_1≤k≤K

{t″_i,k＝t_i+K+k}_1≤k≤K

式中，t′_i,k、t″_i,k分别为移动窗口内前、后两个部分的时间节点；

(3-2)分别对第i个移动窗口内前、后两个部分的时间节点进行标准化处理；标准化方法为：

τ′_i,k＝t′_i,k-t″_i,1 (1)

τ″_i,k＝t″_i,′k-t″_i,1 (2)

式中，τ′_i,k、τ″_i,k分别为t′_i,k、t″_i,k的标准化时间节点；t″_i,1＝t_i+K+1表示拟发生阶跃时间节点；

(3-3)基于GNSS时间序列的趋势不变特性，将第i个移动窗口内前、后两个部分的GNSS时间序列进行保系数趋势项逼近，表达形式如下：

式中，y(τ′)表示第i个移动时间窗口内的前半部分时间序列；y(τ″)表示第i个移动时间窗口内的后半部分时间序列；F_j(·)表示j阶时间基函数，m表示时间基函数的最高阶数，τ′和τ″分别为移动窗口内前、后两部分时间序列的标准化时间节点；β_j表示j阶时间基函数的系数；ε′和ε″分别为移动窗口内前、后两部分时间序列的观测噪声，Step表示移动窗口内前、后两部分时间序列趋势项的阶跃参数。

4.根据权利要求3所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：步骤(4)所述最优估计理论选用稳健最小二乘估计；构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；方法如下：

(4-1-1)构建基于稳健最小二乘估计的GNSS时间序列阶跃探测与修复模型；

设观测噪声服从偶然误差特性，根据时间基函数保系数趋势项逼近分析得到发生阶跃前、后两部分时间序列具有相同的拟合系数，建立GNSS时间序列阶跃探测与修复模型为：

式中，B_i为系数矩阵，L_i为常数向量，V_i表示残差向量，d_i＝[Step,β₀,…,β_m]^T表示第i个移动窗口的待估参数向量；

所述系数矩阵B_i，常数向量L_i，残差向量V_i表示如下：

(4-1-2)利用稳健最小二乘准则求解第i个移动窗口中的待估参数

所述稳健最小二乘准则表示为：

待估参数

计算公式为：

式中，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵；

(4-1-3)将参数估计结果

回代模型，即将

代入式(4)求得残差向量V_i，进而求得单位权方差

(4-1-4)根据方差-协方差传播定律计算模型参数估值

的方差-协方差阵D_i；所述方差-协方差阵计算公式如下：

式中，

为单位权方差，B_i为系数矩阵，

为第i个移动窗口内GNSS时间序列的稳健权矩阵。

5.根据权利要求3所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：步骤(4)所述最优估计理论选用广义最小二乘估计；构建基于广义最小二乘估计的GNSS时间序列阶跃探测与修复模型，并进行阶跃参数估计与方差-协方差估计；方法如下：

(4-2-1)构建广义最小二乘的GNSS时间序列阶跃探测与修复模型；

设移动窗口节点长度为2K的第i个移动窗口的待估参数向量为：

d_i＝[Step,β₀,…,β_m]^T

将参数向量d_i进行扩展，重构待估参数向量为：

根据重构后的待估参数向量

建立GNSS时间序列阶跃探测与修复模型为：

式中，

表示系数矩阵，L_i表示常数向量，C表示约束矩阵，

表示残差向量，V_z表示约束误差向量；

所述系数矩阵

约束矩阵C的转置矩阵C^T表示如下：

式中，I_m+1表示m+1阶单位矩阵；

(4-2-2)利用广义最小二乘准则求解第i个移动窗口中的参数估值

所述广义最小二乘准则表示为：

待估参数

计算公式为：

式中，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵，P_z表示参数权矩阵；

(4-2-3)将阶跃参数估计结果

回代模型，即将

代入式(8)分别计算残差向量

和约束误差向量V_z，并估计单位权方差

(4-2-4)根据方差-协方差传播定律计算阶跃参数估值

的方差-协方差阵

所述方差-协方差阵计算公式如下：

式中，

为单位权方差，

为系数矩阵，

表示第i个移动窗口内GNSS时间序列的稳健权矩阵，P_z表示参数权矩阵，C表示约束矩阵。

6.根据权利要求4或5所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：步骤(5)所述构建正态分布统计量，进行GNSS时间序列阶跃参数的显著性检验，若通过显著性检验，则修复GNSS时间序列阶跃，若未通过显著性检验，则无需修复GNSS时间序列；方法如下：

(5-1)根据参数估值及其方差-协方差估计，构建阶跃参数Step显著性检验的标准正态分布统计量U_i：

式中，Step表示阶跃参数，σ_Step表示阶跃中误差；

(5-2)根据预先选定的置信水平α，求得标准正态分布的接受域和拒绝域，所述接受域和拒绝域分别为式(13)和式(14)：

|U_i|≤u_α/2 (13)

|U_i|>u_α/2 (14)

式中，|·|表示取绝对值，u_α/2为正态分布的上α分位点；

(5-3)当正态分布统计量U_i落在接受域内，即满足式(13)时，表示阶跃参数Step不显著、未通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处没有发生阶跃，即无需修复；

反之，当正态分布统计量U_i落在拒绝域内，即满足式(14)时，表示阶跃参数Step显著、通过显著性检验，则判定第i个移动窗口内时间序列的第K+1个时间节点t_i+K+1处发生阶跃，即需要修复；

(5-4)当阶跃参数通过显著性检验时，对时间节点t_i+K+1以后的时间序列进行修复，公式如下：

式中，y_Step(t)表示修复后的GNSS时间序列。

7.根据权利要求6所述的一种GNSS时间序列阶跃探测与修复方法，其特征在于：

若采用稳健最小二乘的时间序列阶跃探测与修复模型，则步骤(5-1)所述公式(12)中，

若采用广义最小二乘的时间序列阶跃探测与修复模型，则步骤(5-1)所述公式(12)中，

其中，

分别为向量

的第一个元素，D_i(1,1)、

分别为矩阵D_i、

的第一行第一列元素。

Images