留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于ARMA模型的BDS卫星钟差异常值探测及其短期预报

韩松辉 宫轶松 李建文 马朝忠 李新娜 郭淑妹

韩松辉, 宫轶松, 李建文, 马朝忠, 李新娜, 郭淑妹. 基于ARMA模型的BDS卫星钟差异常值探测及其短期预报[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
引用本文: 韩松辉, 宫轶松, 李建文, 马朝忠, 李新娜, 郭淑妹. 基于ARMA模型的BDS卫星钟差异常值探测及其短期预报[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
HAN Songhui, GONG Yisong, LI Jianwen, MA Chaozhong, LI Xinna, GUO Shumei. Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
Citation: HAN Songhui, GONG Yisong, LI Jianwen, MA Chaozhong, LI Xinna, GUO Shumei. Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232

基于ARMA模型的BDS卫星钟差异常值探测及其短期预报

doi: 10.13203/j.whugis20190232
基金项目: 

国家自然科学基金 41774038

国家自然科学基金 41474009

详细信息
    作者简介:

    韩松辉, 博士, 副教授, 主要研究方向为测量导航数据处理。hansonghui@126.com

  • 中图分类号: P228

Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction

Funds: 

The National Natural Science Foundation of China 41774038

The National Natural Science Foundation of China 41474009

More Information
    Author Bio:

    HAN Songhui, PhD, associate professor, specializes in processing of surveying and navigation data. E-mail: hansonghui@126.com

  • 摘要: 星载原子钟在运行过程中会受到恶劣空间环境与设备老化等因素的影响, 使得卫星钟差数据中经常存在异常值, 其中AO(additive outlier)类异常值是钟差序列中常见的一类异常值。结合最大期望算法与自回归滑动平均(autoregressive moving average, ARMA)模型, 提出一种AO类异常值探测算法。该算法可以准确探测孤立AO类异常值与成片AO类异常值, 有效克服了其他算法经常出现的淹没与掩盖现象。在成功探测钟差序列AO类异常值的同时, 该算法可以估计得到精确的ARMA模型, 进而能准确地进行卫星钟差预报。利用仿真数据与北斗卫星钟差实测数据进行计算分析, 结果表明, 所提算法可以精确探测出钟差序列AO类异常值, 并且具有很好的卫星钟差预报效果。
  • 图  1  仿真时间序列的AO类异常值探测

    Figure  1.  Detection of AOs in Simulated Time Series

    图  2  两次差分后的钟差数据与wt的估计结果

    Figure  2.  Clock Error Data After Two Differences and Estimated Result of wt

    图  3  二次差分后的钟差数据与其预报结果的差

    Figure  3.  Difference Between Clock Error Data After Two Differences and Its Prediction Results

    图  4  原始钟差数据与其预报结果

    Figure  4.  Original Clock Error Data and Its Prediction Results

    图  5  原始钟差数据与本文算法预报结果的差

    Figure  5.  Difference Between Original Clock Errors and the Prediction Results of the Proposed Algorithm

    图  6  原始钟差数据与LP模型预报结果的差

    Figure  6.  Difference Between Original Clock Errors and the Prediction Results of LP Model

    图  7  原始钟差数据与QP模型预报结果的差

    Figure  7.  Difference Between Original Clock Errors and the Prediction Results of QP Model

    表  1  仿真时间序列的异常扰动估计值

    Table  1.   Estimations of Abnormal Perturbation in Simulated Time Series

    历元 异常扰动估计值
    100 -12.886
    200 10.826
    201 11.444
    202 13.651
    203 12.865
    204 9.937
    300 -12.324
    400 11.830
    下载: 导出CSV
  • [1] Yang Y X, Xu Y Y, Li J L, et al. Progress and Performance Evaluation of BeiDou Global Navigation Satellite System: Data Analysis Based on BDS-3 Demonstration System[J]. Science China Earth Sciences, 2018, 61(5): 614-624 doi:  10.1007/s11430-017-9186-9
    [2] Karunasinghea D, Liongb S. Chaotic Time Series Prediction with a Global Model: Artificial Neural Network[J]. Journal of Hydrology, 2006, 323(1-4): 92-105 doi:  10.1016/j.jhydrol.2005.07.048
    [3] Senior K L, Ray J R, Beard R L. Characterization of Periodic Variations in the GPS Clocks[J]. GPS Solution, 2008, 12(3): 211-225 doi:  10.1007/s10291-008-0089-9
    [4] Heo Y J, Cho J, Heo M B. Improving Prediction Accuracy of GPS Satellite Clocks with Periodic Variation Behavior[J]. Measurement Science and Technology, 2010, 21(7): 3 001-3 008
    [5] Li X, Dong X, Zheng K, et al. Research of Satellite Clock Error Prediction Based on RBF Neural Network and ARMA Model[C]. China Satellite Navigation Conference (CSNC), Wuhan, China, 2013
    [6] Kim J, Kim M.ARMA Prediction of SBAS Ephemeris and Clock Corrections for Low Earth Orbiting Satellites[J]. International Journal of Aerospace Engineering, 2015, DOI:  10.1155/2015/165178
    [7] 黄博华, 杨勃航, 李锡瑞, 等.顾及一次差分数据结构特征的钟差预报模型[J].武汉大学学报·信息科学版, 2020, DOI: 10.13203/j.whugis20190116

    Huang Bohua, Yang Bohang, Li Xirui, et al. Prediction of Navigation Satellite Clock Bias Considering Structure Characteristics of Single Difference Data[J]. Geomatics and Information Science of Wuhan University, 2020, DOI: 10.13203/j.whugis20190116
    [8] 赵琳琳, 刘万科, 李建龙, 等.Galileo/GPS星载原子钟在轨性能评估与分析[J].测绘地理信息, 2020, 45(1):13-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG202001003.htm

    Zhao Linlin, Liu Wanke, Li Jianlong, et al. Performance Assessment and Analysis of On-Board Atomic Clocks for Galileo and GPS Satellites[J]. Journal of Geomatics, 2020, 45(1):13-19 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXG202001003.htm
    [9] 崔先强, 焦文海.灰色系统模型在卫星钟差预报中的应用[J].武汉大学学报·信息科学版, 2005, 30(5): 447-450 灰色系统模型在卫星钟差预报中的应用

    Cui Xianqiang, Jiao Wenhai. Grey System Model for the Satellite Clock Error Prediction[J]. Geomatics and Information Science of Wuhan University, 2005, 30(5): 447-450 灰色系统模型在卫星钟差预报中的应用
    [10] 陆晓峰, 杨志强, 贾小林, 等.灰色系统理论的优化方法及其在卫星钟差预报中的应用[J].武汉大学学报·信息科学版, 2008, 33(5): 492-495 http://ch.whu.edu.cn/article/id/1584

    Lu Xiaofeng, Yang Zhiqiang, Jia Xiaolin, et al. Parameter Optimization Method of Gray System Theory for the Satellite Clock Error Predicating[J]. Geomatics and Information Science of Wuhan University, 2008, 33(5): 492-495 http://ch.whu.edu.cn/article/id/1584
    [11] 王利, 张勤, 黄观文, 等.基于指数平滑法的GPS卫星钟差预报[J].武汉大学学报·信息科学版, 2017, 42(7): 995-1 001 doi:  10.13203/j.whugis20150089

    Wang Li, Zhang Qin, Huang Guanwen, et al. GPS Satellite Clock Bias Prediction Based on Exponential Smoothing Method[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 995-1 001 doi:  10.13203/j.whugis20150089
    [12] 王宇谱, 吕志平, 王宁, 等.顾及卫星钟随机特性的抗差最小二乘配置钟差预报算法[J].测绘学报, 2016, 45(6): 646-655 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201606003.htm

    Wang Yupu, Lü Zhiping, Wang Ning, et al. Prediction of Navigation Satellite Clock Bias Considering Clock's Stochastic Variation Behavior with Robust Least Square Collocation[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(6): 646-655 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201606003.htm
    [13] 朱陵凤, 吴晓平, 李超, 等.灰色模型在卫星钟差预报中的缺陷分析[J].宇航计测技术, 2007, 27(4): 42-44 doi:  10.3969/j.issn.1000-7202.2007.04.010

    Zhu Lingfeng, Wu Xiaoping, Li Chao, et al. Vice Analysis of Gray System Model for the Satellite Clock Error Predicting[J]. Journal of Astronautic Metrology and Measurement, 2007, 27(4): 42-44 doi:  10.3969/j.issn.1000-7202.2007.04.010
    [14] 朱祥维, 肖华, 雍少为, 等.卫星钟差预报的Kalman算法及其性能分析[J].宇航学报, 2008, 29(3): 966-970, 1 052 doi:  10.3873/j.issn.1000-1328.2008.03.045

    Zhu Xiangwei, Xiao Hua, Yong Shaowei, et al. The Kalman Algorithm Used for Satellite Clock Offset Prediction and Its Performance Analysis[J]. Journal of Astronautics, 2008, 29(3): 966-970, 1 052 doi:  10.3873/j.issn.1000-1328.2008.03.045
    [15] Yang Y X. Some Numerical Prediction Methods for the Wind Speed in the Sea Based on ERS-1 Scattermeter wind Data[J]. Survey Review, 2001, 36(280): 121-131 doi:  10.1179/sre.2001.36.280.121
    [16] Kaloop M R, Hussan M, Kim D. Time-Series Analysis of GPS Measurements for Long-Span Bridge Movements Using Wavelet and Model Prediction Techniques[J]. Advances in Space Research, 2019, 63(11): 3 505-3 521 doi:  10.1016/j.asr.2019.02.027
    [17] Zhang G, Gui Q, Han S, et al. A Bayesian Method of GNSS Cycle Slips Detection Based on ARMA Model[C]. Forum on Cooperative Positioning & Service, Harbin, China, 2017
    [18] 归庆明, 李涛, 衡广辉.时间序列异常值探测的Bayes方法及其在电离层VTEC数据处理中的应用[J].武汉大学学报·信息科学版, 2011, 36(7): 802-806 http://ch.whu.edu.cn/article/id/599

    Gui Qingming, Li Tao, Heng Guanghui. Bayesian Method for Time Series Outliers Detection and Applications in Ionospheric VTEC Data Processing[J]. Geomatics and Information Science of Wuhan University, 2011, 36(7): 802-806 http://ch.whu.edu.cn/article/id/599
    [19] 黄观文. GNSS星载原子钟质量评价及精密钟差算法研究[D].西安: 长安大学, 2012

    Huang Guanwen. Research on Algorithms of Precise Clock Offset and Quality Evaluation of GNSS Satellite Clock[D]. Xi'an: Chang'an University, 2012
    [20] 郭海荣, 杨生, 杨元喜, 等.基于卫星双向时间频率传递进行钟差预报的方法研究[J].武汉大学学报·信息科学版, 2007, 32(1): 43-46 http://ch.whu.edu.cn/article/id/1809

    Guo Hairong, Yang Sheng, Yang Yuanxi, et al. Numerical Prediction Methods for Clock Difference Based on Two-Way Satellite Time and Frequency Transfer Data[J]. Geomatics and Information Science of Wuhan University, 2007, 32(1): 43-46 http://ch.whu.edu.cn/article/id/1809
    [21] Winkler G M R. Introduction to Robust Statistics and Data Filtering[EB/OL]. http://www.wriley.com/ROBSTAT.htm, 1993
    [22] Riley W J. Frequency Jump Detection and Analysis[C]. The 40th Annual Precise Time and Time Interval Systems and Applications Meeting, Reston, Virginia, 2008
    [23] 牛飞, 韩春好, 张义生.导航卫星星载原子钟异常监测分析[J].武汉大学学报·信息科学版, 2009, 34(5): 585-588 http://ch.whu.edu.cn/article/id/1261

    Niu Fei, Han Chunhao, Zhang Yisheng. Analysis and Detection on Atomic Clock Anomaly of Navigation Satellites[J]. Geomatics and Information Science of Wuhan University, 2009, 34(5): 585-588 http://ch.whu.edu.cn/article/id/1261
    [24] Huang G, Qin Z. Real-Time Estimation of Satellite Clock Offset Using Adaptively Robust Kalman Filter with Classified Adaptive Factors[J]. GPS Solutions, 2012, 16(4): 531-539 doi:  10.1007/s10291-012-0254-z
    [25] Abraham B, Box G E P. Bayesian Analysis of Some Outlier Problems in Time Series[J]. Biometrika, 1979, 66(2): 229-236 doi:  10.1093/biomet/66.2.229
    [26] Tsay R S. Time Series Model Specification in the Presence of Outliers[J]. Journal of American Statistical Association, 1986, 393(81): 132-141
    [27] Louni H. Outlier Detection in ARMA Models[J]. Journal of Time Series Analysis, 2005, 6(29): 1 057-1 065
    [28] Benkabou S E, Benabdeslem K, Canitia B. Unsupervised Outlier Detection for Time Series by Entropy and Dynamic Time Warping[J]. Knowledge & Information Systems, 2018, 54(2): 463-486
    [29] 牛飞, 韩春好, 张义生, 等.基于星间链路支持的导航卫星自主完好性检测设计仿真[J].测绘学报, 2011, 40(S1): 73-79 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB2011S1014.htm

    Niu Fei, Han Chuanhao, Zhang Yisheng, et al. Design and Simulation for Satellite Autonomous Integrity Monitoring Based on Inter-Satellite-Links[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(S1): 73-79 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB2011S1014.htm
    [30] Zhang Q, Gui Q. Bayesian Methods for Outliers Detection in GNSS Time Series[J]. Journal of Geodesy, 2013, 87(7): 609-627 doi:  10.1007/s00190-013-0640-5
    [31] 韩松辉, 张国超, 张宁, 等. AR模型中AO类异常值探测及其在GPS卫星钟差预报中的应用[J].测绘学报, 2019, 48(10): 1 225-1 235 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201910004.htm

    Han Songhui, Zhang Guochao, Zhang Ning, et al. New Algorithm for Detecting AO Outliers in AR Model and Its Application in the Prediction of GPS Satellite Clock Errors[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(10): 1 225-1 235 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201910004.htm
    [32] Hamilton J D. Time Series Analysis[M]. Princeton, New Jersey: Princeton University Press, 1994
    [33] 韩松辉, 杜兰, 归庆明, 等.诊断复共线性的特征分析法及其在GEO定轨中的应用[J].测绘学报, 2013, 42(1): 19-26 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201301006.htm

    Han Songhui, Du Lan, Gui Qingming, et al. Characteristics Analysis Approach for Multicollinearity Diagnosis and Its Applications in Orbit Determination of GEO Satellites[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(1): 19-26 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201301006.htm
    [34] Hoerl A E, Kennard R W. Ridge Regression: Biased Estimation for Non-orthogonal Problems[J]. Technometrics, 1970, 12(1): 55-58 doi:  10.1080/00401706.1970.10488634
    [35] Cryer J D, Chan K S. Time Series Analysis with Applications in R[M]. 2nd ed. New York: Springer-Verlag, 2008 Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction
  • [1] 马朝忠, 朱建青, 韩松辉.  基于ARIMA模型的卫星钟差异常值探测的模型选择方法 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 167-172. doi: 10.13203/j.whugis20180230
    [2] 汤俊, 毛文飞.  多尺度ARMA残差修正模型震前电离层TEC异常探测 . 武汉大学学报 ● 信息科学版, 2019, 44(6): 791-798. doi: 10.13203/j.whugis20170321
    [3] 刘立龙, 陈雨田, 黎峻宇, 田祥雨, 贺朝双.  活跃期区域电离层总电子短期预报及适用性分析 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1757-1764. doi: 10.13203/j.whugis20180145
    [4] 王利, 张勤, 黄观文, 田婕.  基于指数平滑法的GPS卫星钟差预报 . 武汉大学学报 ● 信息科学版, 2017, 42(7): 995-1001. doi: 10.13203/j.whugis20150089
    [5] 姜兆英, 刘国林, 于胜文.  病态方程基于Liu估计的一种迭代估计新方法 . 武汉大学学报 ● 信息科学版, 2017, 42(8): 1172-1178. doi: 10.13203/j.whugis20150218
    [6] 林东方, 朱建军.  附有奇异值修正限制的改进的岭估计方法 . 武汉大学学报 ● 信息科学版, 2017, 42(12): 1834-1839. doi: 10.13203/j.whugis20150581
    [7] 魏二虎, 李智强, 龚光裕, 张帅.  极移时间序列模型的拟合与预测 . 武汉大学学报 ● 信息科学版, 2013, 38(12): 1420-1424.
    [8] 刘斌, 龚健雅, 江万寿, 祝小勇.  基于岭参数的谱修正迭代法及其在有理多项式参数求解中的应用 . 武汉大学学报 ● 信息科学版, 2012, 37(4): 399-402.
    [9] 任锴, 宋小勇, 贾小林, 欧阳桂崇.  一种区域卫星导航系统导航策略 . 武汉大学学报 ● 信息科学版, 2012, 37(11): 1356-1359.
    [10] 归庆明, 李涛, 衡广辉.  时间序列异常值探测的Bayes方法及其在电离层VTEC数据处理中的应用 . 武汉大学学报 ● 信息科学版, 2011, 36(7): 802-806.
    [11] 袁修孝, 林先勇.  基于岭估计的有理多项式参数求解方法 . 武汉大学学报 ● 信息科学版, 2008, 33(11): 1130-1133.
    [12] 崔先强, 焦文海.  灰色系统模型在卫星钟差预报中的应用 . 武汉大学学报 ● 信息科学版, 2005, 30(5): 447-450.
    [13] 王振杰, 欧吉坤, 柳林涛.  一种解算病态问题的方法——两步解法 . 武汉大学学报 ● 信息科学版, 2005, 30(9): 821-824.
    [14] 王振杰, 欧吉坤.  L-曲线法确定岭估计中的岭参数 . 武汉大学学报 ● 信息科学版, 2004, 29(3): 235-238.
    [15] 范永弘.  SAR图像的几何校正 . 武汉大学学报 ● 信息科学版, 1997, 22(1): 39-41.
    [16] 王新洲.  在无偏估计类中改进最小二乘估计的方法 . 武汉大学学报 ● 信息科学版, 1995, 20(1): 46-50.
    [17] 邵建平, 沈鹤鸣, 宋光普.  抗SA影响的卡尔曼滤波技术 . 武汉大学学报 ● 信息科学版, 1994, 19(3): 274-279.
    [18] 张方仁, 汪晓庆.  平差参数的岭估计和压缩估计 . 武汉大学学报 ● 信息科学版, 1989, 14(3): 46-58.
    [19] 黄幼才.  岭估计及其应用 . 武汉大学学报 ● 信息科学版, 1987, 12(4): 64-73.
    [20] 李德仁.  自检校平差中过度参数化的克服 . 武汉大学学报 ● 信息科学版, 1986, 11(3): 95-104.
  • 加载中
图(7) / 表(1)
计量
  • 文章访问数:  27
  • HTML全文浏览量:  3
  • PDF下载量:  10
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-09-20
  • 刊出日期:  2021-02-05

基于ARMA模型的BDS卫星钟差异常值探测及其短期预报

doi: 10.13203/j.whugis20190232
    基金项目:

    国家自然科学基金 41774038

    国家自然科学基金 41474009

    作者简介:

    韩松辉, 博士, 副教授, 主要研究方向为测量导航数据处理。hansonghui@126.com

  • 中图分类号: P228

摘要: 星载原子钟在运行过程中会受到恶劣空间环境与设备老化等因素的影响, 使得卫星钟差数据中经常存在异常值, 其中AO(additive outlier)类异常值是钟差序列中常见的一类异常值。结合最大期望算法与自回归滑动平均(autoregressive moving average, ARMA)模型, 提出一种AO类异常值探测算法。该算法可以准确探测孤立AO类异常值与成片AO类异常值, 有效克服了其他算法经常出现的淹没与掩盖现象。在成功探测钟差序列AO类异常值的同时, 该算法可以估计得到精确的ARMA模型, 进而能准确地进行卫星钟差预报。利用仿真数据与北斗卫星钟差实测数据进行计算分析, 结果表明, 所提算法可以精确探测出钟差序列AO类异常值, 并且具有很好的卫星钟差预报效果。

English Abstract

韩松辉, 宫轶松, 李建文, 马朝忠, 李新娜, 郭淑妹. 基于ARMA模型的BDS卫星钟差异常值探测及其短期预报[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
引用本文: 韩松辉, 宫轶松, 李建文, 马朝忠, 李新娜, 郭淑妹. 基于ARMA模型的BDS卫星钟差异常值探测及其短期预报[J]. 武汉大学学报 ● 信息科学版, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
HAN Songhui, GONG Yisong, LI Jianwen, MA Chaozhong, LI Xinna, GUO Shumei. Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
Citation: HAN Songhui, GONG Yisong, LI Jianwen, MA Chaozhong, LI Xinna, GUO Shumei. Outliers Detection in BDS Satellite Clock Errors by Using ARMA Model and Corresponding Short-Term Prediction[J]. Geomatics and Information Science of Wuhan University, 2021, 46(2): 244-251. doi: 10.13203/j.whugis20190232
  • 北斗卫星导航系统(BeiDou navigation satellite system, BDS)目前正处于全面建设阶段, 随着其导航定位技术的发展, 高精度的钟差预报显得至关重要[1]。目前, 钟差预报的方法有多项式模型、谱分析模型、灰色模型、神经网络模型、Kalman滤波模型以及时间序列模型等[2-8]。多项式模型不适于卫星钟差的长期预报[9-11], 谱分析模型在钟差数据相对较少时存在周期确定不准确问题[12], 灰色模型的预报精度受模型指数系数的影响较大[13], 神经网络的训练过程非常耗时且训练结果容易陷入局部极小点[11-12], Kalman滤波模型在钟差数据较少时难以准确确定原子钟运行特性和随机先验信息[12, 14]。时间序列模型能反映卫星钟差数据间的相关性, 目前ARMA(p, q)模型已经在测绘导航中得到广泛的应用[15-16], 但ARMA(p, q)模型存在形式多样、定阶复杂的问题[17-19]。以上用于钟差预报的各种模型均有各自的优缺点与适用条件, 尚需要进一步研究与完善。

    由于星载原子钟所处的恶劣空间环境、传播路径和观测环境等因素的影响, 卫星钟差的历史观测序列可能包含AO(additive outlier)类异常值(AO类异常值是指只对出现异常值历元的观测值产生影响, 而不影响其他历元的观测值)[8, 20], 这严重影响了模型拟合与卫星钟差预报的精度。目前, 已有学者对钟差序列异常值的探测问题进行了研究并提出了一些解决方法, 如中位数法[19-21]、抗差估计法[20]、Allan方差法[22-23]和抗差自适应滤波法[24]等。但已有的方法均存在一定不足, 如需要差分、对异常值的定位存在偏差、无法有效地对异常扰动进行估计等[17]。时间序列模型充分考虑到卫星钟差的趋势性、周期性和随机性等特点, 一些学者已对时间序列异常值探测进行了研究分析[18, 25-28], 结果表明, 时间序列模型在卫星钟差数据AO类异常值探测与钟差预报方面具有明显的优势[5-6]。例如, 牛飞等[29]采用平稳时间序列的非稳态异常检测算法对星载原子钟相位抖动及频率的异常漂移进行了处理; Zhang等[30]将Bayes理论引入到时间序列异常值探测中。但是, 现有算法中, 有的未对异常扰动进行估计, 有的需要经过复杂的抽样过程, 有的需要计算后验概率。

    本文针对时间序列模型用于钟差预报的现状, 引入最大期望(expectation maximization, EM)算法与ARMA模型, 提出一种卫星钟差序列异常值探测与钟差预报的新算法。首先, 基于ARMA模型, 引入识别变量, 构建AO类异常值探测模型; 然后, 运用EM算法给出AO类异常值定位与异常扰动估计的方法, 并通过引入岭估计克服其中出现的复共线性问题, 该方法通过简单迭代计算得到ARMA模型系数、异常值位置和异常扰动估值; 最后, 利用仿真算例与BDS实测卫星钟差算例进行验证, 对单个AO类异常值和成片AO类异常值进行探测处理和卫星钟差预报分析。

    • 设卫星钟差序列数据$\left\{ {{y_t}} \right\}$满足如下ARMA (p, q)模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {{y_t} = {\phi _1}{y_{t - 1}} + {\phi _2}{y_{t - 2}} + \cdot \cdot \cdot + {\phi _p}{y_{t - p}} + }\\ \;\;\;\;{{a_t} + {\theta _1}{a_{t - 1}} + {\theta _2}{a_{t - 2}} + \cdots + {\theta _q}{a_{t - q}}}\\ {{a_t}{\rm{i}}.{\rm{i}}.{\rm{d}}.{\rm{}}N\left( {0, {\sigma ^2}} \right)} \end{array}} \right. $$ (1)

      式中, ${a_t}$为白噪声, 当$s > t$时, $E\left[ {{a_s}{y_t}} \right] = 0$; ${\rm{i}}.{\rm{i}}.{\rm{d}}.$的含义为相互独立同分布; $\mathit{\boldsymbol{\phi}} = {({\phi _1}, {\phi _2} \cdot \cdot \cdot {\phi _p})^{\rm{T}}} $、$\mathit{\boldsymbol{\theta }} = {({\theta _1}, {\theta _2} \cdot \cdot \cdot {\theta _q})^{\rm{T}}}$和${\sigma ^2}$为未知参数。

      引入识别变量标记异常值, 建立如下AO类异常值探测模型[17, 30]:

      $$ \left\{ {\begin{array}{*{20}{l}} {{x_t} = {y_t} + {\delta _t}{w_t}}\\ {{y_t} = {\phi _1}{y_{t - 1}} + {\phi _2}{y_{t - 2}} + \cdot \cdot \cdot + {\phi _p}{y_{t - p}} + {a_t} + }\\ \;\;\;\;\;\;\;{{\theta _1}{a_{t - 1}} + {\theta _2}{a_{t - 2}} + \cdots + {\theta _q}{a_{t - q}}}\\ {{a_t}{\rm{i}}.{\rm{i}}.{\rm{d}}.{\rm{}}N\left( {0, {\sigma ^2}} \right)} \end{array}} \right. $$ (2)

      式中, $\left\{ {{x_t}} \right\}$为观测得到的可能包含AO类异常值的钟差序列; $\left\{ {{y_t}} \right\}$为无法观测的正常钟差序列; ${\delta _t}$为AO类异常值的识别变量, 且${\delta _t}$服从两点分布, 取值为0或1。如果${\delta _t} = 1$, 则${x_t}$为AO类异常值, 相应异常扰动为${w_t}$; 如果${\delta _t} = 0$, 则${x_t}$不是AO类异常值。在上述构造的卫星钟差序列异常值探测模型中, 需要解决模型系数如何计算、异常值位置如何判定、异常扰动如何估计等问题。

    • 针对以上问题, 韩松辉等[31]基于AR(autoregressive)模型提出了简单直接的卫星钟差序列AO异常值探测算法, 并且对异常值的数量比例、探测成功率、估计精度等方面进行统计分析后发现, 该算法具有良好的稳定性和高效性。但是, 此算法无法处理ARMA(p, q)模型q > 0的情况。对于大部分卫星钟差序列而言, 当前历元的卫星钟差值不仅与以前历元的卫星钟差值有关, 还与以前历元的噪声存在一定的依存关系。此时, 采用AR模型处理卫星钟差序列已经不再合适, 需要采用ARMA(p, q)模型。因此, ARMA(p, q)模型在卫星钟差预报中有更重要的应用。

      考虑AO类异常值探测模型(2)的求解, 显然模型(2)严格的似然函数非常复杂, 使用起来非常困难。但是, 若ARMA(p, q)是可逆的, 则Hamilton[32]给出了一种近似的对数似然函数为:

      $$ {\rm{ln}}L = - \frac{{n - p}}{2}{\rm{ln}}\left( {2{\rm{ \mathit{ π} }}} \right) - \frac{{n - p}}{2}{\rm{ln}}{\sigma ^2} - \frac{1}{{2{\sigma ^2}}}\mathop \sum \limits_{t = p + 1}^n a_t^2 $$ (3)

      式中, ${a_t} = {x_t} - {\phi _1}{y_{t - 1}} - {\phi _2}{y_{t - 2}} - \cdots - {\phi _p}{y_{t - p}} - {\theta _1}{a_{t - 1}} - {\theta _2}{a_{t - 2}} - \cdots $$ - {\theta _q}{a_{t - q}} - {\delta _t}{w_t}$, 且${a_p} = {a_{p - 1}} = \cdots = {a_{p - q + 1}} = 0$。

      令$\mathit{\boldsymbol{X}} = {({x_1}, {x_2} \cdots {x_n})^{\rm{T}}}$, 则:

      $$ \begin{array}{*{20}{l}} {\xi = ({\phi _1}, {\phi _2} \cdots {\phi _p}, {\theta _1}, {\theta _2} \cdots {\theta _q}, {\sigma ^2}, {\delta _1}, {\delta _2} \cdots {\delta _n}, }\\ \;\;\;\;\;\;\;\;{{w_1}, {w_2} \cdots {w_n}{)^{\rm{T}}}} \end{array} $$

      式中, n为卫星钟差序列中的观测值个数。则可将文献[31]相关结论扩展到ARMA(p, q)模型, 推导得到AO类异常值探测与钟差预报的算法, 其步骤如下:

      1) 用最小二乘估计法或极大似然估计法得到ARMA(p, q)模型中相关参数${\phi _1}, {\phi _2} \cdots {\phi _p}, {\theta _1}, {\theta _2} \cdots {\theta _q}, {\sigma ^2}$的初始值$\hat \phi _1^0, \hat \phi _2^0 \cdots \hat \phi _p^0, \hat \theta _1^0, \hat \theta _2^0 \cdots \hat \theta _q^0, {({\hat \sigma ^2})^0}$, 并令$\hat \delta _1^0 = \hat \delta _2^0 = \cdots = \hat \delta _n^0 = 0$, $\hat w_1^0 = \hat w_2^0 = \cdots = \hat w_n^0 = 0$, $\hat a_t^0 = {x_t} - \hat \phi _1^0{x_{t - 1}} - \hat \phi _2^0{x_{t - 2}} - \cdots - \hat \phi _p^0{x_{t - p}} - \hat \theta _1^0\hat a_{t - 1}^0 - \hat \theta _2^0\hat a_{t - 2}^0 - \cdots - \hat \theta _q^0\hat a_{t - p}^0$, $\hat y_t^0 = {x_t}$。其中, $t = p + 1, p + 2 \cdots n$, 相关参数的上角标表示迭代次数。

      2) 已知$\delta _t^i$服从0-1分布, 令$P\left( {\delta _t^i = 1} \right) = \hat p_t^i$, $P\left( {\delta _t^i = 0} \right) = 1 - \hat p_t^i$。假设已得到${\xi ^i}$, 若$\hat w_t^i$为钟差序列中包含的随机误差, 则可以近似地得到$\hat w_t^i \sim N(0, \left( {{{\hat \sigma }^2}{)^i}} \right)$, 由正态分布的$k\sigma $原则(根据具体数据类型采用相应的k值)可判断得到的${\hat w_t}$是否是卫星钟差序列的AO类异常扰动。因此, 如果$\hat w_t^i$的绝对值超过$k{\hat \sigma ^i}$, 即:

      $$ \hat w_t^i < - k{\hat \sigma ^i}或\hat w_t^i > k{\hat \sigma ^i} $$ (4)

      则可认为$\hat w_t^i$为AO类异常扰动, ${x_t}$为AO类异常值, 设定$\hat p_t^i = 1$, 否则取$\hat p_t^i = 0$。同时, 得到消除当前历元AO类异常值影响的修正观测值:

      $$ \begin{array}{*{20}{l}} {\hat y_t^i = \hat \phi _1^i\hat y_{t - 1}^i + \hat \phi _2^i\hat y_{t - 2}^i + \cdot \cdot \cdot + \hat \phi _p^i\hat y_{t - p}^i + }\\ \;\;\;\;\;\;{\hat \theta _1^i\hat a_{t - 1}^i + \hat \theta _2^i\hat a_{t - 2}^i + \cdots + \hat \theta _q^i\hat a_{t - q}^i + \hat a_t^i} \end{array} $$ (5)

      3) 令$\hat z_t^i = {x_t} - \hat \phi _1^i{\hat y_{t - 1}} - \hat \phi _2^i{\hat y_{t - 2}} - \cdots - \hat \phi _p^i{\hat y_{t - p}}$, 则:

      $$ \hat w_t^{i + 1} = \hat z_t^i $$ (6)
      $$ \hat a_t^{i + 1} = \hat z_t^i - \hat p_t^i\hat w_t^i $$ (7)
      $$ {({\hat \sigma ^2})^{i + 1}} = \frac{{\sum\limits_{t = p + 1}^n {[{{(\hat z_t^i)}^2} + \hat p_t^i\left( {\hat w_t^i{)^2} - 2\hat p_t^i\hat w_t^i\hat z_t^i} \right]} }}{{n - p}} $$ (8)
      $$ \left( {\begin{array}{*{20}{c}} {\hat \phi _1^{i + 1}}\\ {\hat \phi _2^{i + 1}}\\ \vdots \\ {\hat \phi _p^{i + 1}}\\ {\hat \theta _1^{i + 1}}\\ {\hat \theta _2^{i + 1}}\\ \vdots \\ {\hat \theta _q^{i + 1}} \end{array}} \right) = {\left( {\begin{array}{*{20}{c}} {\mathop \sum \limits_{t = p + 1}^n \hat y_{t - 1}^2}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - p}}{{\hat y}_{t - 1}}}&{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - 1}}{{\hat y}_{t - 1}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - q}}{{\hat y}_{t - 1}}}\\ {\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - 1}}{{\hat y}_{t - 2}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - p}}{{\hat y}_{t - 2}}}&{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - 1}}{{\hat y}_{t - 2}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - q}}{{\hat y}_{t - 2}}}\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ {\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - 1}}{{\hat y}_{t - p}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n \hat y_{t - p}^2}&{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - 1}}{{\hat y}_{t - p}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - q}}{{\hat y}_{t - p}}}\\ {\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - 1}}{{\hat a}_{t - 1}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - p}}{{\hat a}_{t - 1}}}&{\mathop \sum \limits_{t = p + 1}^n \hat a_{t - 1}^2}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - q}}{{\hat a}_{t - 1}}}\\ {\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - 1}}{{\hat a}_{t - 2}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - p}}{{\hat a}_{t - 2}}}&{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - 1}}{{\hat a}_{t - 2}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - q}}{{\hat a}_{t - 2}}}\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ {\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - 1}}{{\hat a}_{t - q}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n {{\hat y}_{t - p}}{{\hat a}_{t - q}}}&{\mathop \sum \limits_{t = p + 1}^n {{\hat a}_{t - 1}}{{\hat a}_{t - q}}}& \cdots &{\mathop \sum \limits_{t = p + 1}^n \hat a_{t - q}^2} \end{array}} \right)^{ - 1}}\left( {\begin{array}{*{20}{c}} {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat y}_{t - 1}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat y}_{t - 1}}}\\ {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat y}_{t - 2}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat y}_{t - 2}}}\\ \vdots \\ {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat y}_{t - p}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat y}_{t - p}}}\\ {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat a}_{t - 1}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat a}_{t - 1}}}\\ {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat a}_{t - 2}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat a}_{t - 2}}}\\ \vdots \\ {\mathop \sum \limits_{t = p + 1}^n {x_t}{{\hat a}_{t - q}} - \mathop \sum \limits_{t = p + 1}^n \hat w_t^i\hat p_t^i{{\hat a}_{t - q}}} \end{array}} \right) $$ (9)

      4)由于式(9)中的系数矩阵在求逆时涉及到${\hat y_t}$和${\hat a_t}$两组数值差别较大的数据, 导致系数矩阵中存在相对于矩阵本身元素绝对值小2~3个数量级的特征值, 从而导致系数矩阵出现复共线性[33]。为此, 考虑采用岭估计克服复共线性对系数矩阵求逆的影响, 在式(9)系数矩阵的每一个对角线元素均加上岭参数k, 然后再求逆。这里的岭参数可利用Hoerl-Kennard公式[34]$k = \frac{{{{\hat \sigma }^2}}}{{{\rm{max}}\hat \alpha _j^2}}$得到, 其中$\mathit{\boldsymbol{\hat \alpha }} = {({\hat \alpha _1}, {\hat \alpha _2} \cdots {\hat \alpha _{p + q}})^{\rm{T}}}$$ = {Q^{\rm{T}}}{(\hat \phi _1^{i + 1}, \hat \phi _2^{i + 1} \cdots \hat \phi _p^{i + 1}, \hat \theta _1^{i + 1}, \hat \theta _2^{i + 1} \cdots \hat \theta _q^{i + 1})^{\rm{T}}}$是未知参数向量${(\hat \phi _1^{i + 1}, \hat \phi _2^{i + 1} \cdots \hat \phi _p^{i + 1}, \hat \theta _1^{i + 1}, \hat \theta _2^{i + 1} \cdots \hat \theta _q^{i + 1})^{\rm{T}}}$的典则参数向量, $Q = \left( {{\mathit{\boldsymbol{q}}_1}, {\mathit{\boldsymbol{q}}_2} \cdots {\mathit{\boldsymbol{q}}_t}} \right)$, ${\mathit{\boldsymbol{q}}_1}, {\mathit{\boldsymbol{q}}_2} \cdots {\mathit{\boldsymbol{q}}_t}$是系数矩阵的标准正交特征向量。

      5) 重复2)~4)步, 直到$\hat \phi _1^i, \hat \phi _2^i \cdots \hat \phi _p^i, \hat \theta _1^i, \hat \theta _2^i \cdots \hat \theta _q^i$与$\hat \phi _1^{i + 1}, \hat \phi _2^{i + 1} \cdots \hat \phi _p^{i + 1}, \hat \theta _1^{i + 1}, \hat \theta _2^{i + 1} \cdots \hat \theta _q^{i + 1}$之间的差别足够小。例如, 可以采用公式$\sqrt {\left( {\hat \phi _1^i - \hat \phi _1^{i + 1}} \right) + \left( {\hat \phi _2^i - \hat \phi _2^{i + 1}} \right) + \cdots + \left( {\hat \theta _q^i - \hat \theta _q^{i + 1}} \right)} \le \varepsilon $判断迭代是否收敛, $\varepsilon $是给定的一个比较小的数。

      6) 假设迭代m次后收敛, 则$\hat p_t^m$可作为判断第t个观测值是否为AO类异常值的依据, 而$\hat w_t^m$的大小可作为异常扰动的估值, $\hat \phi _1^m, \hat \phi _2^m \cdots \hat \phi _p^m, \hat \theta _1^m, \hat \theta _2^m \cdots \hat \theta _q^m$即为估计得到的ARMA(p, q)模型的相应系数。

      上述算法不但可以消除AO类异常值的影响, 而且可以同时得到准确的ARMA(p, q)模型, 并进行钟差预报。在此算法中, 模型系数、异常值识别变量和异常扰动大小表达式均为严格推导所得的结果, 具有理论上的严密性, 同时也不会额外引入不同迭代过程之间的计算误差。本算法估计得到的噪声方差与异常扰动比较准确, 利用式(4)可以准确探测出AO类异常值, 而不会出现淹没与掩盖现象(掩盖现象是指由于某些异常值的出现而导致其他异常值不能被探测出来, 从而无法发现全部异常值, 甚至根本无法发现异常值; 淹没现象是指由于某些异常值的出现而导致某些正常值被误判为异常值, 从而导致异常值定位结果与实际异常值位置不符)。准确判断AO类异常值之后, 基于准确的卫星钟差历史数据, 可得到消除AO类异常值影响的修正钟差数据。修正钟差数据与原始钟差数据一起对异常值进行探测, 使得计算结果既不偏离原始钟差数据, 又可消除AO类异常值的影响。

    • 考虑如下ARMA(2, 1)模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {{y_t} = 0.2{y_{t - 1}} + 0.5{y_{t - 2}} - 0.3{y_{t - 3}} + }\\ \;\;\;\;\;{{a_t} + 0.3{a_{t - 1}} - 0.1{a_{t - 2}} + 0.2{a_{t - 3}}}\\ {{a_t}\;{\rm{i}}.{\rm{i}}.{\rm{d}}.{\rm{}}N\left( {0, 1} \right)} \end{array}} \right. $$

      仿真产生500个观测值, 为了检验本文算法探测AO类异常扰动的效果, 在第200~204个历元的观测值中分别加上大小依次为11、12、13、11、10的成片AO类异常扰动, 并在第100、第300与第400个历元的观测值中分别加上大小为-13、-12与11的孤立AO类异常值。在式(4)中令k= 4, 则增加了AO类异常扰动的时间序列与异常值探测结果如图 1所示, 具体AO类异常扰动的估计值如表 1所示。

      图  1  仿真时间序列的AO类异常值探测

      Figure 1.  Detection of AOs in Simulated Time Series

      表 1  仿真时间序列的异常扰动估计值

      Table 1.  Estimations of Abnormal Perturbation in Simulated Time Series

      历元 异常扰动估计值
      100 -12.886
      200 10.826
      201 11.444
      202 13.651
      203 12.865
      204 9.937
      300 -12.324
      400 11.830

      图 1表 1可以看出, 无论是孤立AO类异常值还是成片AO类异常值, 本文算法均没有出现淹没与掩盖现象。对于这种同时包含孤立AO类异常值与成片AO类异常值的复杂序列, 该算法准确地探测出了AO类异常值的位置, 而且得到了准确的异常扰动估计值。令j=100, 200, 201, 202, 203, 204, 300, 400, 用${w_j}$表示这几个位置上的真实异常扰动, 则$\sqrt {\mathop \sum \limits_j {{({{\hat w}_j} - {w_j})}^2}} $可以表示异常扰动估计精度, 通过计算可得此次仿真的估计精度为2.248。

      为验证本文算法的有效性, 仿真计算100 000次, 异常值位置的成功探测率为96.631%。在异常值位置成功探测的前提下, 异常扰动估计的平均精度为2.988。计算结果说明了本文算法的稳健性与准确性。

    • 在iGMAS(International GNSS Monitoring & Assessment System)官网发布的BDS卫星精密钟差数据中, 提取为期1 d的C09卫星的最终钟差数据(单位为s), 时间为2019-01-22, 钟差数据的采样间隔为5 min, 共包括288个历元。利用本文算法探测BDS卫星钟差序列中的AO类异常值并进行钟差预报。

    • C09卫星的钟差数据经过两次差分后变为平稳序列, 计算发现, 其自相关函数、偏相关函数均不截尾, 综合运用扩展自相关法与赤池信息准则法计算分析[35], 可以认为两次差分后的序列服从ARMA(2, 3)模型。采用本文算法, 在式(4)中令k= 3, 可以探测出差分后钟差序列中的疑似AO类异常值, 计算得到的两次差分后的钟差数据和${w_t}$的计算结果如图 2所示。

      图  2  两次差分后的钟差数据与wt的估计结果

      Figure 2.  Clock Error Data After Two Differences and Estimated Result of wt

      分析图 2可以发现, 本文算法探测出一些疑似异常值, 并同时估计了疑似异常扰动的大小。总共探测出30个AO类疑似异常值, 在这些疑似异常值中, 一部分是成片的疑似异常值, 对于成片疑似AO类异常值, 本文算法没有出现淹没和掩盖现象。所以本文算法不但可以探测孤立疑似AO类异常值, 而且可以克服成片疑似AO类异常值探测时的淹没和掩盖现象。另外, 在计算过程中发现, 算法具有很好的迭代稳定性, 随着迭代次数的增加, 估计精度逐渐增高并趋于稳定。但是, ARMA模型本身比AR模型复杂, 而且还要克服系数矩阵的复共线性影响, 所以此处的计算比AR模型的相关计算要复杂一些。

      另外, 分析图 2发现, 异常扰动估计值接近于二次差分后的钟差值, 但是两者之间有一定的差别。这是合理的, 因为疑似异常扰动估计值只是疑似异常扰动的估计, 不是二次差分后钟差异常值的估计值。二次差分后的钟差值由真实钟差值与疑似异常扰动两部分组成。因为是实测数据, 此处不知道具体的疑似异常扰动大小, 所以只能得到接近于二次差分后钟差值的疑似异常扰动估计值。

    • 为分析卫星C09二次差分后的钟差数据的预报情况, 利用前面数据对后50个历元(约4.2 h)进行预报, 将二次差分后的钟差数据与相应的预报钟差数据作差, 结果如图 3所示。由图 3可以看出, 二次差分后的钟差数据的预报精度非常高, 最大偏差在0.4 ns左右, 充分验证了本文算法的有效性。此处的预报结果在最后几个历元有较大的偏差, 这和异常值探测结果是对应的, 因为在最后几个历元存在疑似AO类异常扰动。预报结果已经消除了疑似异常扰动的影响, 但二次差分后的钟差数据仍存在疑似异常扰动, 因此出现了图 3中最后几个位置预报偏差较大的情况, 出现此情况正说明了本文算法的有效性。为了验证卫星C09差分前初始钟差数据的预报精度, 给出初始钟差数据和反差分运算(差分运算的逆运算)后的预报钟差数据, 如图 4所示。

      图  3  二次差分后的钟差数据与其预报结果的差

      Figure 3.  Difference Between Clock Error Data After Two Differences and Its Prediction Results

      图  4  原始钟差数据与其预报结果

      Figure 4.  Original Clock Error Data and Its Prediction Results

      分析图 4发现, 仅采用238个历元的钟差数据, 预报数据的曲线和iGMAS提供的最终钟差数据的曲线几乎重合。计算结果验证了本文算法对ARMA模型估计的有效性, 即使在包含疑似AO类异常值的情况下, 本文算法不但可以探测出疑似AO类异常值, 而且可以估计得到准确的ARMA模型。为进一步检验预报精度, 将卫星C09初始钟差数据和相应反差分后的预报钟差数据作差, 其结果如图 5所示。同时为了对比分析, 采用一次多项式(linear polynomial, LP)与二次多项式(quadratic polynomial, QP)模型进行钟差预报, 其结果分别如图 6图 7所示。

      图  5  原始钟差数据与本文算法预报结果的差

      Figure 5.  Difference Between Original Clock Errors and the Prediction Results of the Proposed Algorithm

      图  6  原始钟差数据与LP模型预报结果的差

      Figure 6.  Difference Between Original Clock Errors and the Prediction Results of LP Model

      图  7  原始钟差数据与QP模型预报结果的差

      Figure 7.  Difference Between Original Clock Errors and the Prediction Results of QP Model

      图 5可以看出, 本文算法的钟差预报精度与iGMAS提供的钟差之间的差别较小, 利用历史卫星钟差数据进行钟差预报的精度比较高, 预报的最大偏差为0.4 ns。对比图 5图 6图 7发现, 本文算法在探测出疑似AO类异常值的同时, 相应的钟差预报精度高于LP模型与QP模型, 由此说明了本文算法的准确性与高效性。

    • 本文针对BDS卫星钟差数据中经常包含AO类异常值的问题, 提出了一种新的ARMA模型AO类异常值探测算法, 利用该算法进行卫星钟差预报。结果表明, 该算法不但可以处理孤立AO类异常值问题, 而且可以有效消除成片AO类异常值的影响。ARMA模型因为包含自回归过程与滑动平均过程, 模型比较复杂, 所以其计算过程比AR模型要复杂一些。总结起来, 该算法有以下优点:

      1) 本文算法可以同时估计ARMA模型系数、噪声方差、异常值位置和异常扰动大小, 所有待估参数的估值均在一个迭代过程中得到, 具有较好的估计相容性。

      2) 本文算法估计得到的噪声方差非常准确, 只要 $\sigma $的倍数选择得当, 算法可以准确找到卫星钟差序列中单个或成片异常值, 不会出现淹没和掩盖现象。

      3) 每一步迭代过程中均对原始卫星钟差观测进行修正得到纯净数据, 并与原始卫星钟差数据共同参与估计的迭代过程, 进而得到更准确的模型系数、噪声方差和异常扰动的估计值。

      4) 在异常扰动的判定过程中, 不仅考虑了估计的异常扰动大小, 而且考虑了估计的噪声方差, 充分利用了异常扰动的相关统计信息。

      5) 本文算法的迭代稳定性良好, 一般经过几次或几十次迭代即可得到最终的估计结果, 随着迭代次数的增加, 算法逐步收敛。

      6) 对于ARMA模型中$pq \ne 0$的情况, 由于涉及到观测值和噪声两类数据, 容易导致系数矩阵包含复共线性, 采用岭估计等有偏估计可克服复共线性的影响。

      基于本文算法处理的卫星钟差观测序列可以准确探测出观测序列中的疑似AO类异常值, 并同时拟合得到精确的ARMA模型, 从而实现对卫星钟差的高效、快速预报。

参考文献 (35)

目录

    /

    返回文章
    返回