基于偶然误差特性的时间序列粗差探测方法及装置与流程

1.本发明涉及粗差探测方法技术领域，具体涉及一种基于偶然误差特性的时间序列粗差探测方法及装置。


背景技术：

2.在时间序列预处理方面，kleijer等以及pernetti等利用dia(detection identification adaptation)程序对gps坐标时间序列的异常值和中断进行了探测，发现gps坐标时间序列中大约有6％的数据为粗差。目前，应用较广泛的粗差探测方法分为两类，一类是统计量粗差探测法，一类是非统计量粗差探测法。
3.统计量粗差探测法是基于观测数据的误差分布特征，构造与之相适应的粗差判别统计量，采用假设检验对粗差进行探测。其中常用的一个是经典3倍中误差(3σ)准则，另外一个是相对稳健的四分位距(interquartile range，iqr)准则，3σ准则的原理是当某一观测值与时间序列的均值的偏差大于3倍标准差时，则其判定为粗差，该方法当先验模型准确、所含粗差量较少的情况下效果较好，但其抗差能力较差，受粗差大小影响较大；iqr准则是一种稳健估计法，其样本中位数和样本iqr受粗差的影响较小，用它们分别代替统计方法中的平均值和标准偏差来进行总体估计，再通过稳健z比分数统计量来探测数据中的极端异常值(粗差)。在3σ准则和iqr准则的基础上衍生出不少新的算法，wang x等在2015年提出了数据驱动的奇异谱分析技术结合iqr准则的方法(ssa-iqr)，明锋等在2016年提出了一种将l1范数(l1-norm)估计与四分位距统计量iqr(inter quartile range)组合的移动开窗粗差探测算法，杨凯钧等在2015年提出了一种顾及偏态的iqr法则，吴浩等在2019年提出了一种基于小波分析的改进型3σ粗差探测方法，这些方法对传统的3σ准则及iqr准则做了改进，取得了一定效果。
4.非统计量粗差探测法主要以小波分析为主，它从小波变换后小波系数具有的模量极大值特性出发，寻找出观测数据的突变点来实现粗差探测，其在数据统计规律不明确的情况下，对粗差含量少、偏离度较大的数据序列具有一定的探测效果。但小波信号的奇异性分为两个类型：一是信号幅值发生突变，这是信号中的粗差；二是信号虽然光滑但其一阶微分不连续也引起信号小波变换后发生突变，因此难以判断序列奇异点是否是粗差点。
5.目前，对于粗差的处理常用方法有3σ准则和四分位间距(inter quartile range，iqr)准则，这两种方法都属于统计量粗差探测法。3σ准则作为一种经典的粗差探测方法，其利用最小二乘法获得残差序列，再进行粗差探测。当先验模型准确、所含粗差量较少的情况下效果较好。但该方法很大程度上依赖于先验模型的准确度，且最小二乘估计的抗差干扰能力较差，其计算的残差容易受粗差污染，从而导致使用3倍中误差作为限差来剔除粗差的准则失去了功效。iqr准则是一种稳健估计法，其具有较好的抗差性，但对于离散度较大的观测数据，其四分位间距随之变大，使得稳健z比分数统计量对于偏离度较小的粗差判定不够敏感，存在一定程度上的“误判”和“漏判”。因此，标准的iqr准则在gnss时间序列探测上的准确性和可靠性受到质疑。


技术实现要素：

6.为此，本发明实施例提供一种基于偶然误差特性的时间序列粗差探测方法及装置，以解决现有技术中存在的基于3σ准则和四分位间距准则的缺陷和不足的问题。
7.为了实现上述目的，本发明实施例提供如下技术方案：
8.第一方面，一种基于偶然误差特性的时间序列粗差探测方法，包括：
9.根据时间序列拟合模型对原始数据时间序列进行拟合，并得到时间时序的残差序列r；
10.计算残差序列r的绝对值，并将结果进行排序；
11.选取rbin值对残差序列r进行分段；
12.根据残差序列r《rbin分别进行正态分布拟合，并分别求出不同rbin分段的σi；
13.判断是否存在粗差；
14.若不存在粗差，则σi服从同一正态分布，并随着rbin范围所包含的数据增多，σi逐渐趋于固定值；
15.若σi显现出不同趋势，则存在粗差，并对相邻σi进行差分计算；
16.对σi进行插值，并对其去趋势，得到不含趋势项的插值σi；
17.通过局部峰值查找函数，确定查找差分结果表现出不同运动趋势时刻的σ
x
；
18.根据偶然误差分布特性进行粗差去除。
19.进一步的，所述时间序列拟合模型公式为：
20.设任意观测序列y(ti)，则：
[0021][0022]
其中，i为坐标历元时刻标识；y(ti)为测站某一分量ti时刻的坐标；b为截距；v0为线性速度；m0为谐波个数；am和bm是频率为fm时周期项的振幅；为随机噪声，即
[0023]
进一步的，所述计算残差序列r的绝对值，并将结果进行排序时是按照从小到大的顺序进行排序的。
[0024]
进一步的，所述选取rbin值对残差序列r进行分段时，rbin值在gnss高精度解算中，水平方向可取至1mm。
[0025]
进一步的，所述对σi进行插值时可进行不同类型插值。
[0026]
第二方面，一种基于偶然误差特性的时间序列粗差探测装置，包括：
[0027]
拟合模块，用于根据时间序列模型对原始数据时间序列进行拟合，并得到时间时序的残差序列r；
[0028]
排序模块，用于计算残差序列r的绝对值，并将结果进行排序；
[0029]
分段模块，用于选取rbin值对残差序列r进行分段；
[0030]
计算模块，用于根据残差序列r《rbin分别进行正态分布拟合，并分别求出不同rbin分段的σi；
[0031]
判断模块，用于判断是否存在粗差；
[0032]
若σi服从同一正态分布，并随着rbin范围所包含的数据增多，逐渐趋于固定值，则
不存在粗差；
[0033]
若σi显现出不同趋势，则存在粗差并对相邻σi进行差分计算；
[0034]
去趋势模块，用于对σi进行插值，并对其去趋势，得到不含趋势项的插值σi；
[0035]
查找模块，用于通过局部峰值查找函数，确定查找差分结果表现出不同运动趋势时刻的σ
x
；
[0036]
粗差去除模块，用于根据偶然误差分布特性进行粗差去除。
[0037]
进一步的，所述时间序列拟合模型公式为：
[0038]
设任意观测序列y(ti)，则：
[0039][0040]
其中，i为坐标历元时刻标识；y(ti)为测站某一分量ti时刻的坐标；b为截距；v0为线性速度；m0为谐波个数；am和bm是频率为fm时周期项的振幅；为随机噪声，即
[0041]
进一步的，所述计算残差序列r的绝对值，并将结果进行排序时是按照从小到大的顺序进行排序的。
[0042]
第三方面，一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现基于偶然误差特性的时间序列粗差探测方法的步骤。
[0043]
第四方面，一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现基于偶然误差特性的时间序列粗差探测方法的步骤。
[0044]
本发明至少具有以下有益效果：本发明提供了一种基于偶然误差特性的时间序列粗差探测方法及装置，包括：根据时间序列拟合模型对原始数据时间序列进行拟合，得到残差序列r；计算r的绝对值，并将结果进行排序；选取rbin值对r进行分段；根据r《rbin分别进行正态分布拟合，并求出不同rbin分段的σi；若不存在粗差，则σi服从同一正态分布，并随着rbin范围所包含的数据增多，σi逐渐趋于固定值；若σi显现出不同趋势，则存在粗差并对相邻σi进行差分计算；对σi差分后数据进行插值，并对其去趋势；通过波峰查找函数，确定查找差分结果表现出不同运动趋势时刻的σ
x
；根据偶然误差概率分布特性进行粗差去除；本发明基于偶然误差特性改善了由于原始时间序列中粗差污染造成无法精确估算数据粗差的影响，提高了粗差探测的准确性。
附图说明
[0045]
为了更清楚地说明现有技术以及本发明，下面将对现有技术以及本发明实施例描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引申获得其它的附图。
[0046]
本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均
应仍落在本发明所揭示的技术内容能涵盖的范围内。
[0047]
图1为本发明实施例提供的基于偶然误差特性的时间序列粗差探测方法流程图；
[0048]
图2为本发明实施例提供的rbin范围σi分布曲线图；
[0049]
图3为本发明实施例提供的gnss模拟时间序列分布图；其中，a为理想先验模型分布图，b为非理想先验模型；
[0050]
图4为3σ准则粗差探测结果，其中，a为理想先验模型分布图，b为非理想先验模型；
[0051]
图5为iqr准则粗差探测结果，其中，a为理想先验模型分布图，b为非理想先验模型；
[0052]
图6为基于偶然误差特性的时间序列粗差探测结果，其中，a为理想先验模型分布图，b为非理想先验模型；
[0053]
图7为ynmh站n方向时间序列粗差探测结果；
[0054]
图8为ynmh站e方向时间序列粗差探测结果；
[0055]
图9为ynmh站u方向时间序列粗差探测结果；
[0056]
图10为xiag站n方向时间序列粗差探测结果；
[0057]
图11为xiag站e方向时间序列粗差探测结果；
[0058]
图12为xiag站u方向时间序列粗差探测结果。
具体实施方式
[0059]
为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本技术，并不用于限定本技术。
[0060]
在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”、“第三”“第四”等(如果存在)旨在区别指代的对象。对于具有时序流程的方案，这种术语表述方式不必理解为描述特定的顺序或先后次序，对于装置结构的方案，这种术语表述方式也不存在对重要程度、位置关系的区分等。
[0061]
此外，术语“包括”、“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包括了一系列步骤或单元的过程、方法、系统、产品或设备不必限于已明确列出的那些步骤或单元，而是还可包含虽然并未明确列出的但对于这些过程、方法、产品或设备固有的其它步骤或单元，或者基于本发明构思进一步的优化方案所增加的步骤或单元。
[0062]
偶然误差∈n(0,σ2)，即偶然误差服从均值为0，标准差为σ的正态分布。而经统计r出现超过三倍σ的概率为1-99.7％，可以通过超出概率范围的方法去除粗差。但同时，通过偶然误差分布特性得知，如何更加准确的求出σ，是有效准确去除粗差的关键。现在普遍的方法是通过对数据序列求取标准差，并将它作为σ。但因数据中存在粗差，则这样求取的σ往往由于包含粗差，而导致求取的σ结果偏大无法准确的描述原有的n(0,σ2)分布，进而导致粗差探测的失败。
[0063]
请参阅图1，本发明的实施例提供一种基于偶然误差特性的时间序列粗差探测方法，即基于偶然误差特性的正态分布拟合gnss时序粗差探测方法(nls-ndf)，包括：
[0064]
s1：根据时间序列拟合模型对原始数据时间序列进行拟合，并得到时间时序的残
差序列r；
[0065]
具体的，假设任意观测序列y(ti)，对于该序列可以选取一个有效函数对它进行回归，如gnss时间序列中常采用上述时间序列提出的拟合模型，该模型能够有效的对台站的线性运动趋势、周年和半年项、阶跃、震后形变、共模误差噪声等进行拟合。该模型公式如所示：
[0066][0067]
其中，i为坐标历元时刻标识；y(ti)为测站某一分量ti时刻的坐标；b为截距；v0为线性速度；m0为谐波个数；am和bm是频率为fm时周期项的振幅；为随机噪声，即
[0068]
s2：计算残差序列r的绝对值，并将结果进行排序；
[0069]
具体的，计算残差序列r的绝对值，并将结果按从小到大的顺序进行排序。
[0070]
s3：选取rbin值对残差序列r进行分段；
[0071]
具体的，选取合适的rbin值对r进行分段，rbin的选取一般取决于数据序列的理想精度，如gnss高精度解算中，水平方向rbin可取至1mm。
[0072]
s4：根据残差序列r《rbin分别进行正态分布拟合，并分别求出不同rbin分段的σi；
[0073]
具体的，请参阅图2。
[0074]
s5：判断是否存在粗差；
[0075]
若σi服从同一正态分布，并随着rbin范围所包含的数据增多，逐渐趋于固定值，则不存在粗差，如图2中的σi＝3；
[0076]
若σi显现出不同趋势，则存在粗差并对相邻σi进行差分计算；
[0077]
具体的，因粗差的存在，σi必然显现出不同趋势。本方法对相邻σi进行差分计算，且从偶然误差的“绝对值较小的误差，出现的机率较大；绝对值较大的误差，出现的机率较小”统计特性得知，rbin范围较小时求得的σ比rbin大求得的σ更加可靠。
[0078]
s6：对σi进行插值，并对其去趋势，得到不含趋势项的插值σi；
[0079]
具体的，因σi只在固定的n*rnbin点位有值，无法较为精确的获取其变化位置，因此下一步对σi进行插值，并对其去趋势，得到不含趋势项的插值σi。
[0080]
s7：通过局部峰值查找函数，确定查找差分结果表现出不同运动趋势时刻的σ
x
；
[0081]
具体的，此时σ
x
为不含粗差的准确值，后面趋势的变化是因为粗差的加入导致了σi的跳跃。
[0082]
s8：根据偶然误差分布特性进行粗差去除。
[0083]
具体的，得到上面求得的精确σ
x
后，即可按照偶然误差分布特性进行粗差去除。
[0084]
本发明基于偶然误差“绝对值较小的误差，出现的机率较大。绝对值较大的误差，出现的机率较小”的偶然特性来进行粗差探测。通过对偶然误差特性的数据描述及变化过程，更加准确的求出数据中的标准差σ(有效准确去除粗差的关键：因数据中存在粗差，造成其它方法求取的σ往往由于包含粗差，而导致求取的σ结果偏大无法准确的描述原有的n(0,σ2)分布，进而导致粗差探测的失败)。
[0085]
为了更加清楚的说明本发明提供的基于偶然误差特性的时间序列粗差探测方法
的有效性，首先基于时间序列函数模型，生成模拟gnss时间序列。针对是否加入高频周期项，将模拟的数据分为理想先验模型与非理想先验模型。
[0086]
理想先验模型生成模拟数据的参数为：b＝5，v＝0，m＝2，1＝5，b1＝5，a2＝3，b2＝3，一年频率f1＝1/365和半年频率f2＝1/182；而非理想先验模型生成的模拟数据中，加入了高频周期项：a3＝6，b3＝6，f3＝1/30。
[0087]
请参阅图3，模拟过程如下：首先采用标准差为∈n(3σ,6σ)的正态分布得到一组随机误差序列，然后将小于3σ的数据剔除，再将剩余序列任意的加入到原始时间序列中，最终得到被粗差污染的时间序列。模拟的时间序列中粗差总数为219，占总观测历元的10％。
[0088]
请参阅图4至图6，分别利用3σ准则、iqr准则，以及本发明提供的nls-ndf三种方法对模拟的数据进行粗差探测，结果如图4至图6所示。可以看到，3σ准则探测粗差效果最差，只能探测到少部分的粗差，粗差探测率仅为24.7％，iqr准则能探测到大部分粗差，粗差探测率为82.7％，而本发明提供的nls-ndf几乎能探测到所有粗差点，粗差探测率达99.7％。
[0089]
为进一步验证三种粗差探测方法对于不同数值σ粗差的敏感性，本发明分别统计了3σ准则、iqr准则，以及nls-ndf对于粗差1倍至9倍的粗差探测有效率和非粗差误报个数，结果如表一所示。从表一中可以发现：当面临不同数值大小的粗差时，3σ粗差探测方法的平均有效率只有24.5％，iqr和nls-ndf方法的平均有效率分别为84.0％和99.6％，3σ粗差探测方法与iqr和nls-ndf方法有效率要低效许多。三种粗差探测方法的粗差探测效率均较为稳定，伴随σ的增加，并没有较大的波动。误报个数方面，模拟数据中添加的总误差个数为219，可以发现3σ准则、iqr准则和nls-nd粗差探测方法的误报个数均《＝4个，即误报率小于2％。
[0090]
表一：
[0091][0092]
同时，为进一步验证三种粗差探测方法对于原始时序中粗差占比的敏感性，本发明分别统计了3σ准则、iqr准则以及新方法nls-ndf对于粗差占比5％至30％的粗差探测有效率和非粗差误报个数，结果如表二所示。从表二中可以发现：当面临不同占比的粗差时，3σ粗差探测方法和iqr粗差探测方法相比于本发明提出的nls-ndf方法的伴随粗差占比的提升，有效率出现明显的下降。当粗差占原始序列比例为30％时，3σ只能探测到2.44％的误差，iqr只能探测38.97％的误差，而本发明提出的nls-ndf方法伴随粗差占比的升高，有效率并未出现明显的下降。误报个数方面，伴随模拟数据中我们添加的总误差个数增多，可以发现3σ准则、iqr准则和nls-nd粗差探测方法的误报个数均《＝4个，即误报率小于2％，并不存在明显差异。
[0093]
表二：
[0094][0095]
请参阅图7至图9，为了验证本发明提出的nls-ndf的粗差探测效果，本发明采用云南区域的2个itrf2014框架下的区域cors站(ynmh站和xiag站)时间序列作为实测数据。数据时间跨度为2011.0014～2019.9603。图7至图9分别为ynmh站n、e、u方向的3σ、iqr、nls-ndf算法粗差探测实验对比。从上至下依次为ynmh站原始时间序列、3σ粗差探测结果、iqr粗差探测结果、nls-ndf粗差探测结果，其中符号
“×”
为算法探测得到的粗差。
[0096]
从ynmh站粗差探测结果得知，三种方法当面临较大数值粗差时，均能表现出较好的粗差探测能力。但从ynmh站的n方向时间序列结果，可以发现3σ算法并未探测到2014～2015时间跨度内的粗差，这表明iqr和nls-ndf算法在大数值粗差探测的有效性要明显优于3σ算法。
[0097]
请参阅图10至图12，分别为xiag站n、e、u方向的3σ、iqr、nls-ndf算法粗差探测实验对比。
[0098]
从xiag站粗差探测结果得知，nls-ndf算法与3σ算法、iqr算法的粗差探测结果具有较大差异。传统的3σ算法、iqr算法由于其只是简单数学计算算法，无法根据识别时间序列中年变周期信号，而错误的把部分年变信号识别为粗差。而本发明提出的nls-ndf算法在由于其核心思想是偶然误差分布特性，在存在复杂年变周期性变化时间序列的粗差探测过程中有效性和准确性要明显优于3σ算法、iqr算法。
[0099]
通过对模拟数据和实测数据的分析：结果表明相对于传统3σ算法、iqr算法，nls-ndf算法法能够有效地探测到gnss坐标时间序列中的粗差。模拟数据验证了三种粗差探测方法对于不同数值σ粗差的敏感性，对于粗差1倍至9倍的粗差探测有效率和非粗差误报个数，在误报率结果小于2％的情况下，nls-ndf方法的平均有效率99.6％要明显优于3σ粗差探测方法的24.5％和iqr方法的84.0％。同时本发明进一步验证三种粗差探测方法对于原始时序中粗差占比的敏感性问题，当面临不同占比的粗差时，3σ粗差探测方法和iqr粗差探测方法相比于本发明提出的nls-ndf方法的伴随粗差占比的提升，有效率出现明显的下降。实测数据表明，nls-ndf方法由于其核心思想是偶然误差分布特性，在存在复杂年变周期性变化时间序列的粗差探测过程中有效性和准确性要明显优于传统数学粗差算法。nls-ndf算法，既解决了gnss时间序列在大规模粗差影响下无法准确提取粗差的难题，也改善了3σ法等传统粗差算法因残差受到粗差污染影响而无法准确预估中误差所带来的影响，提高了粗差算法对gnss时间序列粗差探测的适应性。
[0100]
虽然本发明所采用的数据为gnss时间序列，但本方法并不只适用于gnss时间序列，本发明提供的方法适用于所有存在偶然误差特性的时间序列粗差探测处理。
[0101]
本发明的实施例一种基于偶然误差特性的时间序列粗差探测装置，包括：
[0102]
拟合模块，用于根据时间序列拟合模型对原始数据时间序列进行拟合，并得到时间时序的残差序列r；
[0103]
排序模块，用于计算残差序列r的绝对值，并将结果从小到大进行排序；
[0104]
分段模块，用于选取rbin值对残差序列r进行分段；
[0105]
计算模块，用于根据残差序列r《rbin分别进行正态分布拟合，并分别求出不同rbin分段的σi；
[0106]
判断模块，用于判断是否存在粗差；
[0107]
若σi服从同一正态分布，并随着rbin范围所包含的数据增多，逐渐趋于固定值，则不存在粗差；
[0108]
若σi显现出不同趋势，则存在粗差并对相邻σi进行差分计算；
[0109]
去趋势模块，用于对σi进行插值，并对其去趋势，得到不含趋势项的插值σi；
[0110]
查找模块，用于通过局部峰值查找函数，确定查找差分结果表现出不同运动趋势时刻的σ
x
；
[0111]
粗差去除模块，用于根据偶然误差分布特性进行粗差去除。
[0112]
具体的，时间序列拟合模型公式为：
[0113]
设任意观测序列y(ti)，则：
[0114][0115]
其中，i为坐标历元时刻标识；y(ti)为测站某一分量ti时刻的坐标；b为截距；v0为线性速度；m0为谐波个数；am和bm是频率为fm时周期项的振幅；为随机噪声，即
[0116]
关于基于偶然误差特性的时间序列粗差探测装置的具体限定可以参见上文中对于基于偶然误差特性的时间序列粗差探测方法的限定，在此不再赘述。
[0117]
本发明的实施例提供一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现基于偶然误差特性的时间序列粗差探测方法的步骤。
[0118]
本发明的实施例提供一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现基于偶然误差特性的时间序列粗差探测方法的步骤。
[0119]
以上几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。
[0120]
以上实施例的各技术特征可以进行任意的组合(只要这些技术特征的组合不存在矛盾)，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述；这些未明确写出的实施例，也都应当认为是本说明书记载的范围。
[0121]
上文中通过一般性说明及具体实施例对本发明作了较为具体和详细的描述。应当指出的是，在不脱离本发明构思的前提下，显然还可以对这些具体实施例作出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。