CN114114332A - 一种有效探测gnss基准站坐标时间序列不连续点的方法 - Google Patents

Abstract

本发明提出了一种有效探测GNSS基准站坐标时间序列不连续点的方法。本发明以GNSS基准站坐标时间序列作为输入数据；对坐标时间序列中相邻时刻的观测值作差得到坐标差值序列，使用四分位距方法对差值序列进行探测得到初始不连续点；构造坐标时间序列的参数估计模型，通过优化求解得到未知参数的估计值及残差；构建检验统计量，对初始不连续点进行统计检验，求解获取新的未知参数估值；检查新的未知参数的估值与初始未知参数的估计值是否相等；输出参数估计值，即为GNSS基准站坐标时间序列中不连续点。本发明无需先验信息，即可实现不连续点的自动探测，可靠性和适应性强，可为基准站坐标时间序列的精细建模和分析提供准确的基础信息。

Description

一种有效探测GNSS基准站坐标时间序列不连续点的方法

技术领域

本发明属于时间序列分析领域，并具体涉及一种有效探测GNSS基准站坐标 时间序列不连续点的方法。

背景技术

GNSS基准站坐标时间序列中的阶跃(Offsets)现象，在时间序列分析中也 被称为不连续点(Discontinuities)，通常是由于接收机天线和天线罩更换、大地 震造成的同震位移、GNSS软件处理策略变化以及一些未知的人为因素和环境因 素引起的。一方面，GNSS时间序列中不连续点的存在会降低GNSS坐标时间序 列的精度，如果未知的不连续点没有被探测并正确修正，则会使速度产生偏差， 进而影响地球物理学估计量，降低GNSS成果的可靠性和其建立及维持地球参考 框架的精度。另一方面，由于时间序列中多个不连续点类似于一种随机游走过程， 不连续点的存在会改变坐标时间序列的噪声特性。因此，不连续点的探测分析已 经成为GNSS成果分析与应用中的重要问题。

GNSS基准站中的不连续点处理，通常依靠基准站事件记录的方式来设定阶 跃参数，例如：地震事件、接收机天线更换或人为事件等，实际处理时在事件发 生的时间点上增加一个待估参数。但事件记录往往存在疏漏，一些未知的偶然事 件也无法一一记录，这时只能从时间序列本身获取缺失的阶跃信息，因此有必要 研究基于数据本身而无需先验信息的阶跃探测方法。随着GNSS基准站网规模的 不断增加、观测数据的不断积累，实现基准站坐标时间序列中不连续性的快速可 靠探测势在必行。

发明内容

本发明主要是针对目前GNSS基准站坐标时间序列分析预处理时，难以准确 自动探测和识别序列中不连续点这一问题，提供一种无需先验信息实现自动探测 与识别序列中的不连续的方法，为后续坐标时间序列的精细建模和分析提供准确 的基础信息，进而获取更为可靠真实的基准站运动特征。

针对上述问题，本发明的主要技术解决方案为：一种有效探测GNSS基准站 坐标时间序列不连续点的方法。所述方法按以下步骤进行：

步骤1，输入GNSS基准站坐标时间序列；

步骤2，将GNSS基准站坐标时间序列中相邻时刻的观测值作差得到GNSS 基准站坐标差值序列，将GNSS基准站坐标差值序列使用四分位距方法进行探测 得到多个初始不连续点；

步骤3，构造GNSS基准站坐标时间序列的参数估计模型，以步骤1所述的 GNSS基准站坐标时间序列作为输入数据，以所述GNSS基准站坐标时间序列的 参数估计模型的截断项、所述GNSS基准站坐标时间序列的参数估计模型的线性 速度、所述GNSS基准站坐标时间序列的参数估计模型的年周期项正弦系数、所 述GNSS基准站坐标时间序列的参数估计模型的年周期项余弦系数、所述GNSS 基准站坐标时间序列的参数估计模型的半年周期项的正弦系数、所述GNSS基准 站坐标时间序列的参数估计模型的半年周期项的余弦系数、步骤2所述的多个初 始不连续点作为未知参数，通过优化求解得到未知参数的最优线性无偏估计值及 所述GNSS基准站坐标时间序列的残差序列；

步骤4，根据步骤3所述的残差序列构建检验统计量，对步骤2所述的多个 初始不连续点进行统计检验，每次从步骤3所述构造的GNSS基准站坐标时间序 列的参数估计模型中剔除不显著的初始不连续点，直到所有的初始不连续点经过 统计量检验都显著为止，以所述GNSS基准站坐标时间序列的参数估计模型的截 断项、所述GNSS基准站坐标时间序列的参数估计模型的线性速度、所述GNSS 基准站坐标时间序列的参数估计模型的年周期项正弦系数、所述GNSS基准站坐 标时间序列的参数估计模型的年周期项余弦系数、所述GNSS基准站坐标时间序 列的参数估计模型的半年周期项的正弦系数、所述GNSS基准站坐标时间序列的 参数估计模型的半年周期项的余弦系数、检验显著的不连续点作为新的未知参数， 通过步骤3所述优化求解得到新的未知参数的估计值及残差序列；

步骤5，检查步骤4获取的新的未知参数的最优线性无偏估计值X^new与步骤 3获取的未知参数的最优无偏估计值X^initial是否相等，若相等则进入步骤6，否则 用X^new替换X^initial重复步骤3～步骤4。

步骤6，输出新的未知参数的最优线性无偏估计值X^new，即为GNSS基准站 坐标时间序列中不连续点。

步骤1所述GNSS基准站坐标时间序列为：

S_i

i∈[1,m₀]

其中，m₀表示观测值的数量，S_i表示第i个时刻观测值；

步骤2所述GNSS基准站坐标差值序列为：

i∈[1,m₀],且

其中，

为第i个时刻观测值与第i-1个时刻观测值之差，m₀表示观测值的数 量，即时刻的数量，S_i表示第i个时刻观测值；

步骤2所述使用四分位距方法进行探测得到多个初始不连续点为：

将GNSS基准站坐标差值序列按从小到大顺序排列，得到排列后GNSS基准 站坐标差值序列；

将排列后GNSS基准站坐标差值序列中位于数据1/2位置处的数据为中位值， 所述中位值用于对排列后GNSS基准站坐标差值序列的集中程度进行度量；

排列后GNSS基准站坐标差值序列位于低四分位点的数值为Q1；

排列后GNSS基准站坐标差值序列位于高四分位点的数值为Q3；

IQR＝Q3-Q1，所述IQR用于对排列后GNSS基准站坐标差值序列的分散程 度进行度量；

NIQR＝IQR×0.7413，表示排列后GNSS基准站坐标差值序列的一倍中误差：

其中，median表示返回给定序列的中值，w为时间窗口；相应的探测条件 为：

满足以上条件的观测值所在序列中的位置即为不连续点出现的位置，也即完 成了初始不连续点的探测。

步骤2中多个初始不连续点为：

D_j

j∈[1,n₀]

其中，n₀表示探测出的不连续点总数，D_j表示第j个初始不连续点；

步骤3所述GNSS基准站坐标时间序列的参数估计模型如下：

式中，S_i为第i个时刻的观测值，t_i为第i个时刻的观测值对应的观测时间， 以年为单位，a为截断项，b为线性速度，c为年周期项正弦系数，d为年周期 项余弦系数，e为半年周期项的正弦系数，f为半年周期项的余弦系数；n₀表示 不连续点个数，D_j为第j个初始不连续点，T_j为第j个不连续点出现的观测时刻， H为海维西特阶梯函数，

v_i为第i个时刻的模型残差，i∈[1,T₀]， j＝1,2,…,n₀，T₀表示观测值的数量，n₀表示不连续点个数。

步骤3所述通过优化求解得到未知参数的估计值及残差序列，具体如下：

步骤3所述的未知参数为：

步骤3所述的GNSS基准站坐标时间序列的参数估计模型进一步转换为：

式中，A为系数矩阵，

为残差序列，v_i为第i个时刻的所述 GNSS基准站坐标时间序列对应的残差值，T₀表示观测值的数量，n₀表示不连续 点个数，i∈[1,T₀]；

所述A的定义如下：

其中，T_j，j＝1,2,…,n₀为第j个不连续点在所述GNSS基准站坐标时间序列 S_i中的位置，n₀表示不连续点个数，有n₀＜m₀，m₀表示GNSS基准站坐标时间 序列的数量；

步骤3所述未知参数的最优线性无偏估计值为：

X^*＝(A^TA)^-1A^TS

其中，X^*表示未知参数的最优线性无偏估计值，A表示系数矩阵，S为所 述GNSS基准站坐标时间序列，S＝S_i,i∈[1,m₀]

步骤3所述残差序列为：

V＝S-AX^*

其中，V表示残差序列，X^*表示未知参数的最优线性无偏估计值，A表示 系数矩阵，S为所述GNSS基准站坐标时间序列；

作为优选，步骤4构建检验统计量如下：

其中，T为统计量，V为步骤3获得的残差序列，U_rem为设置的阈值，一般 U_rem＝0.01，当T小于指定的阈值U_rem时，说明该参数不显著，参数有无对结果的 影响不大，可予以剔除；

本发明具有如下优点：

不依赖准确的先验信息输入，通过模型构建和参数显著性检验，能够自动探 测和剔除不显著参数，识别指定阈值标准以上的未知阶跃，具有较强的可靠性；

遵循从整体到局部，从大到小的原则，实现不连续性的自动探测，增强模型 的适应性。

附图说明

图1：方法流程图

图2：WUHN站坐标时间序列及探测结果；

具体实施方式

下面结合附图和相关实施案例，对本发明的技术方案做进一步的叙述。

实施案例

本发明的实施案例提出了一种有效探测GNSS基准站坐标时间序列不连续 点的方法。所述方法按以下步骤进行：

步骤1，输入GNSS基准站坐标时间序列；

步骤1所述GNSS基准站坐标时间序列为：

S_i

i∈[1,m₀]

其中，m₀＝6075表示观测值的数量，S_i表示第i个时刻观测值；

步骤2，将GNSS基准站坐标时间序列中相邻时刻的观测值作差得到GNSS 基准站坐标差值序列，将GNSS基准站坐标差值序列使用四分位距方法进行探测 得到多个初始不连续点；

步骤2所述GNSS基准站坐标差值序列为：

i∈[1,m₀],且

其中，

为第i个时刻观测值与第i-1个时刻观测值之差，m₀表示观测值的数 量，即时刻的数量，m₀＝6075表示观测值的数量，S_i表示第i个时刻观测值；

步骤2所述使用四分位距方法进行探测得到多个初始不连续点为：

将GNSS基准站坐标差值序列按从小到大顺序排列，得到排列后GNSS基准 站坐标差值序列；

将排列后GNSS基准站坐标差值序列中位于数据1/2位置处的数据为中位值， 所述中位值用于对排列后GNSS基准站坐标差值序列的集中程度进行度量；

排列后GNSS基准站坐标差值序列位于低四分位点的数值为Q1；

排列后GNSS基准站坐标差值序列位于高四分位点的数值为Q3；

IQR＝Q3-Q1，所述IQR用于对排列后GNSS基准站坐标差值序列的分散程 度进行度量；

NIQR＝IQR×0.7413，表示排列后GNSS基准站坐标差值序列的一倍中误差：

其中，median表示返回给定序列的中值，w为时间窗口；相应的探测条件 为：

满足以上条件的观测值所在序列中的位置即为不连续点出现的位置，也即完 成了初始不连续点的探测。

步骤2中多个初始不连续点为：

D_j

j∈[1,n₀]

其中，n₀表示探测出的不连续点总数，n₀＝6，D_j表示第j个初始不连续点；

步骤3，构造GNSS基准站坐标时间序列的参数估计模型，以步骤1所述的 GNSS基准站坐标时间序列作为输入数据，以所述GNSS基准站坐标时间序列的 参数估计模型的截断项、所述GNSS基准站坐标时间序列的参数估计模型的线性 速度、所述GNSS基准站坐标时间序列的参数估计模型的年周期项正弦系数、所 述GNSS基准站坐标时间序列的参数估计模型的年周期项余弦系数、所述GNSS 基准站坐标时间序列的参数估计模型的半年周期项的正弦系数、所述GNSS基准 站坐标时间序列的参数估计模型的半年周期项的余弦系数、步骤2所述的多个初 始不连续点作为未知参数，通过优化求解得到未知参数的最优线性无偏估计值及 所述GNSS基准站坐标时间序列的残差序列；

步骤3所述GNSS基准站坐标时间序列的参数估计模型如下：

式中，S_i为第i个时刻的观测值，t_i为第i个时刻的观测值对应的观测时间， 以年为单位，a为截断项，b为线性速度，c为年周期项正弦系数，d为年周期 项余弦系数，e为半年周期项的正弦系数，f为半年周期项的余弦系数；n₀表示 不连续点个数，D_j为第j个初始不连续点，T_j为第j个不连续点出现的观测时刻， H为海维西特阶梯函数，

v_i为第i个时刻的模型残差，i∈[1,T₀]， j＝1,2,…,n₀，T₀表示观测值的数量，n₀表示不连续点个数。

步骤3所述通过优化求解得到未知参数的估计值及残差序列，具体如下：

步骤3所述的未知参数为：

步骤3所述的GNSS基准站坐标时间序列的参数估计模型进一步转换为：

式中，A为系数矩阵，

为残差序列，v_i为第i个时刻的所述 GNSS基准站坐标时间序列对应的残差值，T₀表示观测值的数量，n₀表示不连续 点个数，i∈[1,T₀]；

所述A的定义如下：

其中，T_j，j＝1,2,…,n₀为第j个不连续点在所述GNSS基准站坐标时间序列 S_i中的位置，n₀表示不连续点个数，有n₀＜m₀，m₀表示GNSS基准站坐标时间 序列的数量；

步骤3所述未知参数的最优线性无偏估计值为：

X^*＝(A^TA)^-1A^TS

其中，X^*表示未知参数的最优线性无偏估计值，A表示系数矩阵，S为所 述GNSS基准站坐标时间序列，S＝S_i,i∈[1,m₀]

步骤3所述残差序列为：

V＝S-AX^*

其中，V表示残差序列，X^*表示未知参数的最优线性无偏估计值，A表示 系数矩阵，S为所述GNSS基准站坐标时间序列；

步骤4，根据步骤3所述的残差序列构建检验统计量，对步骤2所述的多个 初始不连续点进行统计检验，每次从步骤3所述构造的GNSS基准站坐标时间序 列的参数估计模型中剔除不显著的初始不连续点，直到所有的初始不连续点经过 统计量检验都显著为止，以所述GNSS基准站坐标时间序列的参数估计模型的截 断项、所述GNSS基准站坐标时间序列的参数估计模型的线性速度、所述GNSS 基准站坐标时间序列的参数估计模型的年周期项正弦系数、所述GNSS基准站坐 标时间序列的参数估计模型的年周期项余弦系数、所述GNSS基准站坐标时间序 列的参数估计模型的半年周期项的正弦系数、所述GNSS基准站坐标时间序列的 参数估计模型的半年周期项的余弦系数、检验显著的不连续点作为新的未知参数， 通过步骤3所述优化求解得到新的未知参数的估计值及残差序列，新的参数估计 值记为X^new；

步骤4构建检验统计量如下：

其中，T为统计量，V为步骤3获得的残差序列，U_rem为设置的阈值，一般 U_rem＝0.01，当T小于指定的阈值U_rem时，说明该参数不显著，参数有无对结果的 影响不大，可予以剔除；

步骤5，检查步骤4获取的新的未知参数的最优线性无偏估计值X^new与步骤 3获取的未知参数的最优无偏估计值X^initial是否相等，若相等则进入步骤6，否则 用X^new替换X^initial重复步骤3～步骤4。

步骤6，输出新的未知参数的最优线性无偏估计值X^new，即为GNSS基准站 坐标时间序列中不连续点。

WUHN站的坐标时间序列来自SOPAC分析中心，参数显著性水平保持控制 U_rem＝0.01，序列水平方向设置探测的最小阶跃为10mm，最小速率变化为 5mm/yr，序列垂直方向能够探测的最小阶跃为10mm，最小速率变化为10mm/yr， 周期信号的振幅最小值为1mm。利用本发明算法探测不连续点的结果如图2所 示，圆点表示原始时间序列，曲线为拟合曲线，垂直短虚线表示粗差点，垂直长 虚线表示阶跃及阶跃搜索函数曲线。探测的结果如表1所示，N、E方向序列不 存在明显的阶跃和速率变化，仅存在少部分粗差点。U方向存在两处阶跃，其中 历元2002.0699(即2002年1月26)发生明显的阶跃，该时间恰好更换了天线 罩，且天线高垂向改正数比原来大0.041m，分析认为可能是由于当时天线信息 未及时更新，分析中心解算时仍沿用之前的天线高计算，从而造成发生阶跃现象。 而另一处阶跃发生在历元2008.8511(即2008年11月7日)，设备更换表在当天 未有设备更换事件发生，此阶跃产生的原因尚未可知，但从图2中可见，该阶跃 是显著存在的。

表1：为本发明算法探测的不连续点信息及其估值。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属领 域的技术人员可以对所描述的具体实施例替换成其他基准站序列，做各种各样的 修改或补充，或采用相似方式替代，但并不会偏离本发明的精神或者超越所附权 利要求书所定义的范围。

Claims (4)

1.一种有效探测GNSS基准站坐标时间序列不连续点的方法，其特征在于，包括以下步骤：

步骤1，输入GNSS基准站坐标时间序列；

步骤2，将GNSS基准站坐标时间序列中相邻时刻的观测值作差得到GNSS基准站坐标差值序列，将GNSS基准站坐标差值序列使用四分位距方法进行探测得到多个初始不连续点；

步骤3，构造GNSS基准站坐标时间序列的参数估计模型，以步骤1所述的GNSS基准站坐标时间序列作为输入数据，以所述GNSS基准站坐标时间序列的参数估计模型的截断项、所述GNSS基准站坐标时间序列的参数估计模型的线性速度、所述GNSS基准站坐标时间序列的参数估计模型的年周期项正弦系数、所述GNSS基准站坐标时间序列的参数估计模型的年周期项余弦系数、所述GNSS基准站坐标时间序列的参数估计模型的半年周期项的正弦系数、所述GNSS基准站坐标时间序列的参数估计模型的半年周期项的余弦系数、步骤2所述的多个初始不连续点作为未知参数，通过优化求解得到未知参数的最优线性无偏估计值及所述GNSS基准站坐标时间序列的残差序列；

步骤4，根据步骤3所述的残差序列构建检验统计量，对步骤2所述的多个初始不连续点进行统计检验，每次从步骤3所述构造的GNSS基准站坐标时间序列的参数估计模型中剔除不显著的初始不连续点，直到所有的初始不连续点经过统计量检验都显著为止，以所述GNSS基准站坐标时间序列的参数估计模型的截断项、所述GNSS基准站坐标时间序列的参数估计模型的线性速度、所述GNSS基准站坐标时间序列的参数估计模型的年周期项正弦系数、所述GNSS基准站坐标时间序列的参数估计模型的年周期项余弦系数、所述GNSS基准站坐标时间序列的参数估计模型的半年周期项的正弦系数、所述GNSS基准站坐标时间序列的参数估计模型的半年周期项的余弦系数、检验显著的不连续点作为新的未知参数，通过步骤3所述优化求解得到新的未知参数的估计值及残差序列；

步骤5，检查步骤4获取的新的未知参数的最优线性无偏估计值X^new与步骤3获取的未知参数的最优无偏估计值X^initial是否相等，若相等则进入步骤6，否则用X^new替换X^initial重复步骤3～步骤4；

步骤6，输出新的未知参数的最优线性无偏估计值X^new，即为GNSS基准站坐标时间序列中不连续点。

2.根据权利要求1所述的有效探测GNSS基准站坐标时间序列不连续点的方法，其特征在于，步骤1所述GNSS基准站坐标时间序列为：

S_i

i∈[1,m₀]

其中，m₀表示观测值的数量，S_i表示第i个时刻观测值；

步骤2所述GNSS基准站坐标差值序列为：

i∈[1,m₀],且

其中，

为第i个时刻观测值与第i-1个时刻观测值之差，m₀表示观测值的数量，即时刻的数量，S_i表示第i个时刻观测值；

步骤2所述使用四分位距方法进行探测得到多个初始不连续点为：

将GNSS基准站坐标差值序列按从小到大顺序排列，得到排列后GNSS基准站坐标差值序列；

将排列后GNSS基准站坐标差值序列中位于数据1/2位置处的数据为中位值，所述中位值用于对排列后GNSS基准站坐标差值序列的集中程度进行度量；

排列后GNSS基准站坐标差值序列位于低四分位点的数值为Q1；

排列后GNSS基准站坐标差值序列位于高四分位点的数值为Q3；

IQR＝Q3-Q1，所述IQR用于对排列后GNSS基准站坐标差值序列的分散程度进行度量；

NIQR＝IQR×0.7413，表示排列后GNSS基准站坐标差值序列的一倍中误差：

其中，median表示返回给定序列的中值，w为时间窗口；相应的探测条件为：

满足以上条件的观测值所在序列中的位置即为不连续点出现的位置，也即完成了初始不连续点的探测；

步骤2中多个初始不连续点为：

D_j

j∈[1,n₀]

其中，n₀表示探测出的不连续点总数，D_j表示第j个初始不连续点。

3.根据权利要求1所述的有效探测GNSS基准站坐标时间序列不连续点的方法，其特征在于，步骤3所述GNSS基准站坐标时间序列的参数估计模型如下：

式中，S_i为第i个时刻的观测值，t_i为第i个时刻的观测值对应的观测时间，以年为单位，a为截断项，b为线性速度，c为年周期项正弦系数，d为年周期项余弦系数，e为半年周期项的正弦系数，f为半年周期项的余弦系数；n₀表示不连续点个数，D_j为第j个初始不连续点，T_j为第j个不连续点出现的观测时刻，H为海维西特阶梯函数，

v_i为第i个时刻的模型残差，i∈[1,T₀]，j＝1,2,…,n₀，T₀表示观测值的数量，n₀表示不连续点个数；

步骤3所述通过优化求解得到未知参数的估计值及残差序列，具体如下：

步骤3所述的未知参数为：

步骤3所述的GNSS基准站坐标时间序列的参数估计模型进一步转换为：

式中，A为系数矩阵，

为残差序列，v_i为第i个时刻的所述GNSS基准站坐标时间序列对应的残差值，T₀表示观测值的数量，n₀表示不连续点个数，i∈[1,T₀]；

所述A的定义如下：

其中，T_j，j＝1,2,…,n₀为第j个不连续点在所述GNSS基准站坐标时间序列S_i中的位置，n₀表示不连续点个数，有n₀＜m₀，m₀表示GNSS基准站坐标时间序列的数量；

步骤3所述未知参数的最优线性无偏估计值为：

X^*＝(A^TA)^-1A^TS

其中，X^*表示未知参数的最优线性无偏估计值，A表示系数矩阵，S为所述GNSS基准站坐标时间序列，S＝S_i,i∈[1,m₀]

步骤3所述残差序列为：

V＝S-AX^*

其中，V表示残差序列，X^*表示未知参数的最优线性无偏估计值，A表示系数矩阵，S为所述GNSS基准站坐标时间序列。

4.根据权利要求1所述的有效探测GNSS基准站坐标时间序列不连续点的方法，其特征在于，步骤4构建检验统计量如下：

其中，T为统计量，V为步骤3获得的残差序列，U_rem为设置的阈值，一般U_rem＝0.01，当T小于指定的阈值U_rem时，说明该参数不显著，参数有无对结果的影响不大，可予以剔除。

Images