時間序列異常值探測的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