时间序列异常值探测的Bayes方法及其GNSS数据质量控制中的应用

王 琰，张倩倩，车通宇，刘 永

（1.解放军信息工程大学 地理空间信息学院，郑州 450052；2.北京卫星导航中心，北京 100094）

时间序列异常值探测的Bayes方法及其GNSS数据质量控制中的应用

王 琰1,2，张倩倩1，车通宇1，刘 永1

（1.解放军信息工程大学 地理空间信息学院，郑州 450052；2.北京卫星导航中心，北京 100094）

GNSS卫星精密定轨定位的观测量验后残差分析是数据质量控制的重要组成部分。目前的方法都是根据经验设定异常值判断的阈值，没有利用验后残差序列的时间序列特性。引入基于识别变量的AR模型异常值探测的Bayes方法对验后残差序列中的异常值进行探测；将异常值探测问题转化为识别变量后验概率值的计算问题，并给出明确的判别阈值，通过后验概率值与事先给定的阈值进行比较判别出异常值的位置；运用Gibbs抽样设计了识别变量后验概率值的计算方法和异常值的估算方法。实测数据对该算法的验证表明该算法能够准确地探测出异常值的位置，将异常值剔除之后的观测值序列应用到精密单点定位中，结果表明该算法的使用提高了精密单点定位的定位精度。

验后残差；异常值探测；Bayes方法；AR模型；精密单点定位

GNSS数据处理中，质量控制是一个很重要的部分，质量控制的优劣影响最终定位、定轨的精度。完整的质量控制主要包括两部分：数据的预处理与验后残差分析[1-3]。数据预处理不能完全探测出GNSS观测数据中的异常值，残余的异常值必然对最终参数估计的结果产生影响，因此就需要从观测量的验后残差序列中发现残余的异常点，确保参与估计的观测量都是干净的，从而加速卫星精密定轨定位参数估计的收敛。与观测量的预处理相比，观测量的验后残差中不包含一些系统误差，比如对流层误差，因此呈现的是偶然误差的特性。参数估计一次迭代之后，观测量的验后残差是一组随着历元规律而平滑变化的序列，验后残差分析也就是对这组GNSS时间序列进行分析。

观测量的验后残差分析是综合数据处理软件中一个很重要的模块，比如GAMIT、PANDA、BERNESE。从本质上讲，基于验后残差的异常值探测问题就是如何从这组GNSS时间序列观测值中寻找出可能存在的异常值并进行有效的处理。但是目前对于GNSS时间序列观测值的异常值探测大多数是进行简单地统计分析，凭经验设定阈值，未能借助于当代时间序列建模理论与分析技术对GNSS时间序列异常值进行有效的探测。因此，寻求有效地抵制GNSS时间序列异常值影响的策略，以便获得更精确、更稳健、更可靠的GNSS导航定轨定位成果以适应社会发展的需要，显得越来越重要。

目前处理GNSS时间序列异常值的方法大致可分为两大类：一类是基于抗差估计思想的异常值处理方法。例如：杨元喜将抗差估计用于处理GNSS时间序列中的异常值[4]；Nikolaidis基于抗差估计的思想设计了处理 GPS时间序列中异常值的 IQR（Interquartile Range）准则[5]；Cannavo et.al提出了GPS时间序列异常值的抗差自适应实时处理方法[6]；张恒璟和程鹏飞给出了 GPS高程时间序列的抗差探测和插补研究方法[7]；Khodabandeh et.al基于渐近M估计对GPS定位时间序列进行了分析[8]。另一类是以统计假设检验为基础的GNSS时间序列异常值定位、定值和修复的探测方法。例如：Perfetti[9]等学者利用 DIA（Detection, Identification and Adaptation）方法设计的GNSS时间序列异常值探测方法，并应用于GPS参考框架时间序列或 GPS永久性监测站坐标时间序列中异常值的探测；de Lacy[10]、Borghi[12]、Benciolini[13]等学者基于Gauss-Markov模型提出的GNSS时间序列异常值探测的Bayes方法，并应用于GNSS载波相位观测序列中的周跳探测、意大利某地区GPS监测站坐标时间序列中异常值的探测以及地震带中板块漂移现象的分析，以及在其他应用领域提出的一些方法[14-16]。

上述这些GNSS时间序列异常值探测方法各有其特点和相应的应用场合，但是应用以上方法对 GNSS观测量验后残差的时间序列异常值进行分析存在以下两个问题：未充分考虑或利用GNSS时间序列的先验信息；往往凭经验设定野值剔除的阈值。因此，充分利用先验信息，研究基于Bayes统计理论的时间序列异常值探测方法是必要的，也是可行的。

鉴于此，本文将时间序列异常值探测的Bayes方法运用到GNSS观测量验后残差序列的异常值处理中来。首先，采用Box-Jenkins建模方法对验后残差的时间序列进行建模；再次，运用Bayes统计理论，提出时间序列异常值探测的Bayes方法以期对验后残差序列进行异常值处理；最后，设计了两个算例来验证本文算法的有效性。首先，通过精密单点定位一次迭代之后得到的观测量的验后残差序列应用本文算法探测异常值，实测数据表明本文的方法能够准确的探测出异常值的位置；其次，采用剔除异常值的观测值序列进行第二次迭代，结果表明本文算法提高了精密单点定位最终的定位精度，也证明了本文算法是有效的。

1 观测值验后残差的时间序列建模

时间序列分析的建模方法有很多种，本文采用Box-Jenkins建模方法，具体的建模流程如图1所示。

图 1 观测值验后残差序列建模流程图Fig.1 Flowchart of post-fit residuals analysis on observations

从图1可知，Box-Jenkins方法的时间序列建模主要包括原始观测序列的平稳性检验、模型类型的识别和定阶、参数估计[15-16]。

由于AR(p)模型较之ARMA(p,q)模型在实际处理中更为方便，且我们可以对ARMA(p,q)序列通过一个高阶的AR模型来近似逼近，继而再做相关处理。所以，我们不妨假设，经建模处理后的验后残差序列符合如下的AR(p)模型：

2 时间序列异常值探测的Bayes方法

2.1 异常值探测的Bayes模型和准则

式中： 和2σ为未知参数，假定其先验分布为其中λ为超参数，ta为白噪声，且当s t＞ 时，中的异常值，对每个观测值tx引入一个识别变量：

0,,,V φ ν

为基于上述AR模型探测和识别观测值

并假设每个观测值tx为异常值的先验概率都为α ，即由此构造下列探测模型：

t式中：tz代表正常观测值，tw代表异常扰动的大小，并假设这里2μ ξ、 为超参数，

t = 1,… ,n。

为判别观测值是否为异常值以及确定判别阈值，构造如下Bayes假设检验问题[15]：

H ：为正常观测值；1H：为异常值

根据Bayes假设检验的思想和原理，当备选假设1H对应的后验概率

0 大于原假设H0对应的

2.2 基于 Gibbs抽样的识别变量后验概率和异常值的估计

如上所述，异常值探测问题就归结为计算识别变

式中，

有了以上的完全条件分布，我们就可以进行Gibbs抽样。假设()2()()()()r R= … )为用Gibbs抽样算法从上述完全条件分布中抽取的样本，则由式(9)可得识别变量后验概率值的计算公式：

φ σ δ

r

、 、 、 ( 1,,

rr r

W

其中，他参数先验分布按照类似的方法确定。超参数文本根据经验Bayesian估计确定的，例如参数α取0.05代表的含义是观测数据中含有粗差的比例为 5%。本文在算例中给出这些超参数的一组具体取值如下：

3 试验与分析

为验证本文算法的效果，采用 GPS实测数据验证，收集了2014-10-27日(年积日300)全球分布的10个IGS监测站的数据进行静态精密单点定位实验，采样间隔30 s，测站分布如图3。

进而，根据 Bayes点估计原理，由式(10)可得异常扰动的估值公式：

2.3 异常值探测Bayes方法的实施过程

异常值探测Bayes方法的实施流程如图2所示。

该算法实施过程中最为关键的是确定先验分布中的超参数，本文参数的先验分布是根据共轭准则确定的：以参数φ为例，首先由十几个历元的观测值*X得出参数φ的似然函数，然后选取与似然函数具有相同核的分布作为参数φ的先验分布；其

图3 地面站分布图Fig.3 Distribution of receivers

图2 算法的实施流程图

Fig.2 Process of the algorithm

3.1 时间序列异常值探测的Bayes方法的性能试验

以ASPA站的结果为例，其他站的结果与该站的结果类似。图4为ASPA站对PRN01星一次迭代之后获得的一天的随高度角变化的残差序列。从图4可知，该星在一天中只有一个模糊度，时段是134~870历元，采用本文的算法分析这个模糊度时段的观测值中是否含有异常值。

为了验证算法的有效性，取高度角大于80°的数据进行模拟试验（第369~437历元），选择这段时间的数据原因是该段数据期间高度角大，观测数据受多路径误差的影响较小，观测量的验后残差序列在这段时间内是平缓的。在第400历元人为加入0.1 m的粗差，采用本文提出的算法进行异常值探测，结果如图5。可以看出，第400历元的识别变量后验概率值为1.0，根据判别规则，验后概率值大于0.5判断为异常值，可知第400历元的观测量为异常值，说明算法能够识别验后残差序列的中的异常值，因此本文算法是有效的。

图4 ASPA对PRN01星的验后残差序列Fig.4 Post-fit residuals series of ASPA to PRN01

图5 加入粗差后的验后残差序列和识别变量后验概率值Fig.5 Post-fit residuals series after adding outlier and posterior probability values of identification variable

3.2 定位结果比较

将3.1节的ASPA对PRN01星的原始残差序列按照本文提出的算法进行异常值探测，探测结果和原始的残差序列如图6所示。

从图6可以看出，有36个历元对应的识别变量后验概率值大于0.5，如图6中红色星号点标示的位置，其它位置的识别变量后验概率值都远远小于0.5，根据判别规则可知这36个历元的观测值为异常值。由于在卫星高度角较低（卫星升起和降落）时，观测数据受多路径影响较大，观测量的噪声较大，此时观测量的验后残差较大，异常值剔除的数量会相应的增多，说明算法是有效的。

为了检验该算法能否对最终的定位结果有所提高，对上述异常值剔除之后，再进行一次迭代。与不经过异常值处理的精密单点定位结果进行比较，两个方法的结果都分别与IGS周解进行比较，10个测站的结果如表1。

从表1的结果可以看出，应用本文算法后的精密单点定位结果都有一定改善，大部分的改善都在mm量级上，而且是对于 U方向改善明显，这也说明了本文算法的使用能够提高精密单点定位最终的精度和可靠性。

表1 两种算法的精密单点定位结果对比Tab.1 Comparison on the two algorithms for precise point positioning

4 结 论

本文引入时间序列异常值探测的 Bayes方法对验后残差数据的质量控制问题进行研究和探讨，并利用全球分布的10个IGS监测站的数据进行静态精密单点定位实验，验证算法的有效性。主要工作和创新点如下：

1） 首次将基于识别变量的AR模型异常值探测的Bayes方法成功地应用于验后残差数据的优化处理中，并验证了新方法的可行性和有效性；

2） 根据Bayes假设检验原理确定出了异常值的判别阈值，通过比较这些识别变量后验概率值与事先给定的阈值来进行异常值探测，有效地克服了以往探测方法的模糊性及探测标准选择困难的问题；

3） 在正态—Gamma先验分布下，基于Gibbs抽样算法，提出了识别变量后验概率值的计算方法和异常值的估算方法；

4） 利用全球分布的10个IGS监测站的数据进行静态精密单点定位实验，并对异常值剔除前后的定位结果进行了分析。

试验表明本文算法能够较好地保障定位结果的可靠性。

（References）：

[1] 李敏. 多模GNSS融合精密定轨理论及其应用研究[D].武汉大学, 2014.

Li Min. Research on multi-GNSS precise orbit determination theory and application[D]. Wuhan University, 2014. [2] 蔡华, 赵齐乐, 孙汉荣, 等. GNSS实时数据质量控制[J]. 武汉大学学报(信息科学版), 2011, 36(7): 1094-1097.

Cai Hua, Zhao Qi-le, Sun Han-rong, et al. GNSS realtime data quality control[J]. Geomatics and Information Science of Wuhan University, 2011, 36(7): 1094-1097. [3] 张小红, 郭斐, 李盼, 等. GNSS精密单点定位中的实时质量控制[J]. 武汉大学学报(信息科学版), 2012, 37(8): 940-944.

Zhang Xiao-hong, Guo Fei, LI Pan, et al. Real-time quality control procedure for GNSS precise point positioning[J]. Geomatics and Information Science of Wuhan University, 2012, 37(8): 940-944.

[4] Yang Y, Gao W, Zhang X. Robust Kalman fi ltering with constraints: a case study for integrated navigation[J]. Journal of Geodesy, 2010, 84: 373-381.

[5] Nikolaidis R. Observations of geodetic and seismic deformation with the global positioning system[D]. PhD Thesis, University of California, San Diego, 2002: 4470-4775.

[6] Cannavo F, Mattia M, Rossi M, et al. A new algorithm for automatic outlier detection in GPS time series[J]. Journal of Geophysical Research, 2010, 12: 5027.

[7] 张恒璟, 程鹏飞. 基于高程时间序列粗差的抗差探测与插补研究[J]. 大地测量学与测量工程, 2011, 31(4): 71-75. Zhang Heng-jing, Cheng Peng-fei. Study on robust detection and interpolation from gross errors of GPS height time series[J]. Journal of Geodesy and Geodynamics, 2011, 31(4): 71-75.

[8] Khodabandeh A, Amiri-Simkooei A R, Sharifi M A. GPS position time series analysis based on asymptotic normality of M-estimation[J]. Journal of Geodesy, 2012, 86: 15-33.

[9] Perfetti N. Detection of station coordinate discontinuities within the Italian GPS fiducial network[J]. Journal of Geodesy, 2006, 80: 381-396.

[10] de Lacy M C, Reguzzoni M, Sanso F, et al. The Bayesian detection of discontinuities in a polynomial regression and its application to the cycle slip problem[J]. Journal of Geodesy, 2008, 82: 527-542.

[11] Sabadini R, Aoudia A, Barzaghi R, et al. First evidences of fast creeping on a long lasting quiescent earthquake fault in Italy[J]. Geophysical Journal International, 2009, 179(2): 720-732.

[12] Borghi A, Cannizzaro L, Vitti A. Advanced techniques for discontinuity detection in GNSS coordinate time series: an Italian case study[C]//International Association of Geodesy Symposia. 2011, Vol. 136: 627-633.

[13] Benciolini B, Reguzzoni M, Venuti G, et al. Bayesian and variational methods for discontinuity detection: theory overview and performance comparison[M]//VII Hotine-Marussi Symposium on Mathematical Geodesy. International Association of Geodesy Symposia. 2012, Vol. 137: 147-152.

[14] Gui Q, Li X, Gong Y, Li B, et al. A Bayesian unmasking method for locating multiple gross errors based on posterior probabilities of classif i cation variables[J]. Journal of Geodesy, 2011, 85: 191-203.

[15] Zhang Qian-qian, Gui Qing-ming, Li Jian-wen, et al. Bayes method for cycle slips detection based on autoregressive model[C]//China Satellite Navigation Conference (CSNC) 2012 Proceedings: Lecture Notes in Electrical Engineering. 2012, Vol. 160: 317-335.

[16] Zhang Qian-qian, Gui Qing-ming. Bayesian methods for outliers detection in GNSS time series[J]. Journal of Geodesy, 2013, 87: 609-627.

Bayesian outlier-detection method based on autoregressive model for post-fit residuals analysis in GNSS

WANG Yan1,2, ZHANG Qian-qian1, CHE Tong-yu1, LIU Yong1
(1. Institute of Surveying and Mapping, Information Engineering University, Zhengzhou 450052, China; 2. China Beijing Satellite Navigation Center, Beijing 100094, China)

The post-fit residuals analyses on the measurements of GNSS positioning and orbit determination play an important role in the data quality control. In view that the present outlier-detection methods are based on the experience to set the threshold value, an Bayesian detection method for outliers based on autoregressive model is used for outlier detection of post-fit residuals, which takes into account the features of post-fit residuals time series. The outlier detection is changed to calculate the posterior probabilities of classification variables, and the threshold for outlier detection is explicitly given in this method. Gibbs sampler is used to compute the posterior probabilities of classification variables in outlier detection. The validation with the measured data shows that the proposed algorithm can accurately detect the location of the outliers, and clean observations can be used for precise point positioning. The validation result shows that the proposed algorithm can improve the accuracy of precise point positioning.

post-fit residuals; outlier detection; Bayesian detection method; autoregressive model; precise point positioning

P227

A

1005-6734(2016)01-0045-06

10.13695/j.cnki.12-1222/o3.2016.01.009

2015-11-18；

2016-01-25

国家自然科学基金（41104047，41174026）

王琰（1990—），男，博士研究生，从事测量数据处理理论与方法研究。E-mail: wang1yan.hi@163.com