王 琰，張倩倩，車通宇，劉 永

（1.解放軍信息工程大學(xué) 地理空間信息學(xué)院，鄭州 450052；2.北京衛(wèi)星導(dǎo)航中心，北京 100094）

時(shí)間序列異常值探測的Bayes方法及其GNSS數(shù)據(jù)質(zhì)量控制中的應(yīng)用

王 琰1,2，張倩倩1，車通宇1，劉 永1

（1.解放軍信息工程大學(xué) 地理空間信息學(xué)院，鄭州 450052；2.北京衛(wèi)星導(dǎo)航中心，北京 100094）

GNSS衛(wèi)星精密定軌定位的觀測量驗(yàn)后殘差分析是數(shù)據(jù)質(zhì)量控制的重要組成部分。目前的方法都是根據(jù)經(jīng)驗(yàn)設(shè)定異常值判斷的閾值，沒有利用驗(yàn)后殘差序列的時(shí)間序列特性。引入基于識(shí)別變量的AR模型異常值探測的Bayes方法對(duì)驗(yàn)后殘差序列中的異常值進(jìn)行探測；將異常值探測問題轉(zhuǎn)化為識(shí)別變量后驗(yàn)概率值的計(jì)算問題，并給出明確的判別閾值，通過后驗(yàn)概率值與事先給定的閾值進(jìn)行比較判別出異常值的位置；運(yùn)用Gibbs抽樣設(shè)計(jì)了識(shí)別變量后驗(yàn)概率值的計(jì)算方法和異常值的估算方法。實(shí)測數(shù)據(jù)對(duì)該算法的驗(yàn)證表明該算法能夠準(zhǔn)確地探測出異常值的位置，將異常值剔除之后的觀測值序列應(yīng)用到精密單點(diǎn)定位中，結(jié)果表明該算法的使用提高了精密單點(diǎn)定位的定位精度。

驗(yàn)后殘差；異常值探測；Bayes方法；AR模型；精密單點(diǎn)定位

GNSS數(shù)據(jù)處理中，質(zhì)量控制是一個(gè)很重要的部分，質(zhì)量控制的優(yōu)劣影響最終定位、定軌的精度。完整的質(zhì)量控制主要包括兩部分：數(shù)據(jù)的預(yù)處理與驗(yàn)后殘差分析[1-3]。數(shù)據(jù)預(yù)處理不能完全探測出GNSS觀測數(shù)據(jù)中的異常值，殘余的異常值必然對(duì)最終參數(shù)估計(jì)的結(jié)果產(chǎn)生影響，因此就需要從觀測量的驗(yàn)后殘差序列中發(fā)現(xiàn)殘余的異常點(diǎn)，確保參與估計(jì)的觀測量都是干凈的，從而加速衛(wèi)星精密定軌定位參數(shù)估計(jì)的收斂。與觀測量的預(yù)處理相比，觀測量的驗(yàn)后殘差中不包含一些系統(tǒng)誤差，比如對(duì)流層誤差，因此呈現(xiàn)的是偶然誤差的特性。參數(shù)估計(jì)一次迭代之后，觀測量的驗(yàn)后殘差是一組隨著歷元規(guī)律而平滑變化的序列，驗(yàn)后殘差分析也就是對(duì)這組GNSS時(shí)間序列進(jìn)行分析。

觀測量的驗(yàn)后殘差分析是綜合數(shù)據(jù)處理軟件中一個(gè)很重要的模塊，比如GAMIT、PANDA、BERNESE。從本質(zhì)上講，基于驗(yàn)后殘差的異常值探測問題就是如何從這組GNSS時(shí)間序列觀測值中尋找出可能存在的異常值并進(jìn)行有效的處理。但是目前對(duì)于GNSS時(shí)間序列觀測值的異常值探測大多數(shù)是進(jìn)行簡單地統(tǒng)計(jì)分析，憑經(jīng)驗(yàn)設(shè)定閾值，未能借助于當(dāng)代時(shí)間序列建模理論與分析技術(shù)對(duì)GNSS時(shí)間序列異常值進(jìn)行有效的探測。因此，尋求有效地抵制GNSS時(shí)間序列異常值影響的策略，以便獲得更精確、更穩(wěn)健、更可靠的GNSS導(dǎo)航定軌定位成果以適應(yīng)社會(huì)發(fā)展的需要，顯得越來越重要。

目前處理GNSS時(shí)間序列異常值的方法大致可分為兩大類：一類是基于抗差估計(jì)思想的異常值處理方法。例如：楊元喜將抗差估計(jì)用于處理GNSS時(shí)間序列中的異常值[4]；Nikolaidis基于抗差估計(jì)的思想設(shè)計(jì)了處理 GPS時(shí)間序列中異常值的 IQR（Interquartile Range）準(zhǔn)則[5]；Cannavo et.al提出了GPS時(shí)間序列異常值的抗差自適應(yīng)實(shí)時(shí)處理方法[6]；張恒璟和程鵬飛給出了 GPS高程時(shí)間序列的抗差探測和插補(bǔ)研究方法[7]；Khodabandeh et.al基于漸近M估計(jì)對(duì)GPS定位時(shí)間序列進(jìn)行了分析[8]。另一類是以統(tǒng)計(jì)假設(shè)檢驗(yàn)為基礎(chǔ)的GNSS時(shí)間序列異常值定位、定值和修復(fù)的探測方法。例如：Perfetti[9]等學(xué)者利用 DIA（Detection, Identification and Adaptation）方法設(shè)計(jì)的GNSS時(shí)間序列異常值探測方法，并應(yīng)用于GPS參考框架時(shí)間序列或 GPS永久性監(jiān)測站坐標(biāo)時(shí)間序列中異常值的探測；de Lacy[10]、Borghi[12]、Benciolini[13]等學(xué)者基于Gauss-Markov模型提出的GNSS時(shí)間序列異常值探測的Bayes方法，并應(yīng)用于GNSS載波相位觀測序列中的周跳探測、意大利某地區(qū)GPS監(jiān)測站坐標(biāo)時(shí)間序列中異常值的探測以及地震帶中板塊漂移現(xiàn)象的分析，以及在其他應(yīng)用領(lǐng)域提出的一些方法[14-16]。

上述這些GNSS時(shí)間序列異常值探測方法各有其特點(diǎn)和相應(yīng)的應(yīng)用場合，但是應(yīng)用以上方法對(duì) GNSS觀測量驗(yàn)后殘差的時(shí)間序列異常值進(jìn)行分析存在以下兩個(gè)問題：未充分考慮或利用GNSS時(shí)間序列的先驗(yàn)信息；往往憑經(jīng)驗(yàn)設(shè)定野值剔除的閾值。因此，充分利用先驗(yàn)信息，研究基于Bayes統(tǒng)計(jì)理論的時(shí)間序列異常值探測方法是必要的，也是可行的。

鑒于此，本文將時(shí)間序列異常值探測的Bayes方法運(yùn)用到GNSS觀測量驗(yàn)后殘差序列的異常值處理中來。首先，采用Box-Jenkins建模方法對(duì)驗(yàn)后殘差的時(shí)間序列進(jìn)行建模；再次，運(yùn)用Bayes統(tǒng)計(jì)理論，提出時(shí)間序列異常值探測的Bayes方法以期對(duì)驗(yàn)后殘差序列進(jìn)行異常值處理；最后，設(shè)計(jì)了兩個(gè)算例來驗(yàn)證本文算法的有效性。首先，通過精密單點(diǎn)定位一次迭代之后得到的觀測量的驗(yàn)后殘差序列應(yīng)用本文算法探測異常值，實(shí)測數(shù)據(jù)表明本文的方法能夠準(zhǔn)確的探測出異常值的位置；其次，采用剔除異常值的觀測值序列進(jìn)行第二次迭代，結(jié)果表明本文算法提高了精密單點(diǎn)定位最終的定位精度，也證明了本文算法是有效的。

1 觀測值驗(yàn)后殘差的時(shí)間序列建模

時(shí)間序列分析的建模方法有很多種，本文采用Box-Jenkins建模方法，具體的建模流程如圖1所示。

圖 1 觀測值驗(yàn)后殘差序列建模流程圖Fig.1 Flowchart of post-fit residuals analysis on observations

從圖1可知，Box-Jenkins方法的時(shí)間序列建模主要包括原始觀測序列的平穩(wěn)性檢驗(yàn)、模型類型的識(shí)別和定階、參數(shù)估計(jì)[15-16]。

由于AR(p)模型較之ARMA(p,q)模型在實(shí)際處理中更為方便，且我們可以對(duì)ARMA(p,q)序列通過一個(gè)高階的AR模型來近似逼近，繼而再做相關(guān)處理。所以，我們不妨假設(shè)，經(jīng)建模處理后的驗(yàn)后殘差序列符合如下的AR(p)模型：

2 時(shí)間序列異常值探測的Bayes方法

2.1 異常值探測的Bayes模型和準(zhǔn)則

式中： 和2σ為未知參數(shù)，假定其先驗(yàn)分布為其中λ為超參數(shù)，ta為白噪聲，且當(dāng)s t＞ 時(shí)，中的異常值，對(duì)每個(gè)觀測值tx引入一個(gè)識(shí)別變量：

0,,,V φ ν

為基于上述AR模型探測和識(shí)別觀測值

并假設(shè)每個(gè)觀測值tx為異常值的先驗(yàn)概率都為α ，即由此構(gòu)造下列探測模型：

t式中：tz代表正常觀測值，tw代表異常擾動(dòng)的大小，并假設(shè)這里2μ ξ、 為超參數(shù)，

t = 1,… ,n。

為判別觀測值是否為異常值以及確定判別閾值，構(gòu)造如下Bayes假設(shè)檢驗(yàn)問題[15]：

H ：為正常觀測值；1H：為異常值

根據(jù)Bayes假設(shè)檢驗(yàn)的思想和原理，當(dāng)備選假設(shè)1H對(duì)應(yīng)的后驗(yàn)概率

0 大于原假設(shè)H0對(duì)應(yīng)的

2.2 基于 Gibbs抽樣的識(shí)別變量后驗(yàn)概率和異常值的估計(jì)

如上所述，異常值探測問題就歸結(jié)為計(jì)算識(shí)別變

式中，

有了以上的完全條件分布，我們就可以進(jìn)行Gibbs抽樣。假設(shè)()2()()()()r R= … )為用Gibbs抽樣算法從上述完全條件分布中抽取的樣本，則由式(9)可得識(shí)別變量后驗(yàn)概率值的計(jì)算公式：

φ σ δ

r

、 、 、 ( 1,,

rr r

W

其中，他參數(shù)先驗(yàn)分布按照類似的方法確定。超參數(shù)文本根據(jù)經(jīng)驗(yàn)Bayesian估計(jì)確定的，例如參數(shù)α取0.05代表的含義是觀測數(shù)據(jù)中含有粗差的比例為 5%。本文在算例中給出這些超參數(shù)的一組具體取值如下：

3 試驗(yàn)與分析

為驗(yàn)證本文算法的效果，采用 GPS實(shí)測數(shù)據(jù)驗(yàn)證，收集了2014-10-27日(年積日300)全球分布的10個(gè)IGS監(jiān)測站的數(shù)據(jù)進(jìn)行靜態(tài)精密單點(diǎn)定位實(shí)驗(yàn)，采樣間隔30 s，測站分布如圖3。

進(jìn)而，根據(jù) Bayes點(diǎn)估計(jì)原理，由式(10)可得異常擾動(dòng)的估值公式：

2.3 異常值探測Bayes方法的實(shí)施過程

異常值探測Bayes方法的實(shí)施流程如圖2所示。

該算法實(shí)施過程中最為關(guān)鍵的是確定先驗(yàn)分布中的超參數(shù)，本文參數(shù)的先驗(yàn)分布是根據(jù)共軛準(zhǔn)則確定的：以參數(shù)φ為例，首先由十幾個(gè)歷元的觀測值*X得出參數(shù)φ的似然函數(shù)，然后選取與似然函數(shù)具有相同核的分布作為參數(shù)φ的先驗(yàn)分布；其

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

圖2 算法的實(shí)施流程圖

Fig.2 Process of the algorithm

3.1 時(shí)間序列異常值探測的Bayes方法的性能試驗(yàn)

以ASPA站的結(jié)果為例，其他站的結(jié)果與該站的結(jié)果類似。圖4為ASPA站對(duì)PRN01星一次迭代之后獲得的一天的隨高度角變化的殘差序列。從圖4可知，該星在一天中只有一個(gè)模糊度，時(shí)段是134~870歷元，采用本文的算法分析這個(gè)模糊度時(shí)段的觀測值中是否含有異常值。

為了驗(yàn)證算法的有效性，取高度角大于80°的數(shù)據(jù)進(jìn)行模擬試驗(yàn)（第369~437歷元），選擇這段時(shí)間的數(shù)據(jù)原因是該段數(shù)據(jù)期間高度角大，觀測數(shù)據(jù)受多路徑誤差的影響較小，觀測量的驗(yàn)后殘差序列在這段時(shí)間內(nèi)是平緩的。在第400歷元人為加入0.1 m的粗差，采用本文提出的算法進(jìn)行異常值探測，結(jié)果如圖5?？梢钥闯觯?00歷元的識(shí)別變量后驗(yàn)概率值為1.0，根據(jù)判別規(guī)則，驗(yàn)后概率值大于0.5判斷為異常值，可知第400歷元的觀測量為異常值，說明算法能夠識(shí)別驗(yàn)后殘差序列的中的異常值，因此本文算法是有效的。

圖4 ASPA對(duì)PRN01星的驗(yàn)后殘差序列Fig.4 Post-fit residuals series of ASPA to PRN01

圖5 加入粗差后的驗(yàn)后殘差序列和識(shí)別變量后驗(yàn)概率值Fig.5 Post-fit residuals series after adding outlier and posterior probability values of identification variable

3.2 定位結(jié)果比較

將3.1節(jié)的ASPA對(duì)PRN01星的原始?xì)埐钚蛄邪凑毡疚奶岢龅乃惴ㄟM(jìn)行異常值探測，探測結(jié)果和原始的殘差序列如圖6所示。

從圖6可以看出，有36個(gè)歷元對(duì)應(yīng)的識(shí)別變量后驗(yàn)概率值大于0.5，如圖6中紅色星號(hào)點(diǎn)標(biāo)示的位置，其它位置的識(shí)別變量后驗(yàn)概率值都遠(yuǎn)遠(yuǎn)小于0.5，根據(jù)判別規(guī)則可知這36個(gè)歷元的觀測值為異常值。由于在衛(wèi)星高度角較低（衛(wèi)星升起和降落）時(shí)，觀測數(shù)據(jù)受多路徑影響較大，觀測量的噪聲較大，此時(shí)觀測量的驗(yàn)后殘差較大，異常值剔除的數(shù)量會(huì)相應(yīng)的增多，說明算法是有效的。

為了檢驗(yàn)該算法能否對(duì)最終的定位結(jié)果有所提高，對(duì)上述異常值剔除之后，再進(jìn)行一次迭代。與不經(jīng)過異常值處理的精密單點(diǎn)定位結(jié)果進(jìn)行比較，兩個(gè)方法的結(jié)果都分別與IGS周解進(jìn)行比較，10個(gè)測站的結(jié)果如表1。

從表1的結(jié)果可以看出，應(yīng)用本文算法后的精密單點(diǎn)定位結(jié)果都有一定改善，大部分的改善都在mm量級(jí)上，而且是對(duì)于 U方向改善明顯，這也說明了本文算法的使用能夠提高精密單點(diǎn)定位最終的精度和可靠性。

表1 兩種算法的精密單點(diǎn)定位結(jié)果對(duì)比Tab.1 Comparison on the two algorithms for precise point positioning

4 結(jié) 論

本文引入時(shí)間序列異常值探測的 Bayes方法對(duì)驗(yàn)后殘差數(shù)據(jù)的質(zhì)量控制問題進(jìn)行研究和探討，并利用全球分布的10個(gè)IGS監(jiān)測站的數(shù)據(jù)進(jìn)行靜態(tài)精密單點(diǎn)定位實(shí)驗(yàn)，驗(yàn)證算法的有效性。主要工作和創(chuàng)新點(diǎn)如下：

1） 首次將基于識(shí)別變量的AR模型異常值探測的Bayes方法成功地應(yīng)用于驗(yàn)后殘差數(shù)據(jù)的優(yōu)化處理中，并驗(yàn)證了新方法的可行性和有效性；

2） 根據(jù)Bayes假設(shè)檢驗(yàn)原理確定出了異常值的判別閾值，通過比較這些識(shí)別變量后驗(yàn)概率值與事先給定的閾值來進(jìn)行異常值探測，有效地克服了以往探測方法的模糊性及探測標(biāo)準(zhǔn)選擇困難的問題；

3） 在正態(tài)—Gamma先驗(yàn)分布下，基于Gibbs抽樣算法，提出了識(shí)別變量后驗(yàn)概率值的計(jì)算方法和異常值的估算方法；

4） 利用全球分布的10個(gè)IGS監(jiān)測站的數(shù)據(jù)進(jìn)行靜態(tài)精密單點(diǎn)定位實(shí)驗(yàn)，并對(duì)異常值剔除前后的定位結(jié)果進(jìn)行了分析。

試驗(yàn)表明本文算法能夠較好地保障定位結(jié)果的可靠性。

（References）：

[1] 李敏. 多模GNSS融合精密定軌理論及其應(yīng)用研究[D].武漢大學(xué), 2014.

Li Min. Research on multi-GNSS precise orbit determination theory and application[D]. Wuhan University, 2014. [2] 蔡華, 趙齊樂, 孫漢榮, 等. GNSS實(shí)時(shí)數(shù)據(jù)質(zhì)量控制[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 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精密單點(diǎn)定位中的實(shí)時(shí)質(zhì)量控制[J]. 武漢大學(xué)學(xué)報(bào)(信息科學(xué)版), 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] 張恒璟, 程鵬飛. 基于高程時(shí)間序列粗差的抗差探測與插補(bǔ)研究[J]. 大地測量學(xué)與測量工程, 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

國家自然科學(xué)基金（41104047，41174026）

王琰（1990—），男，博士研究生，從事測量數(shù)據(jù)處理理論與方法研究。E-mail: wang1yan.hi@163.com