-
导航星座要长时间实现定位授时功能，星座时间必须与地面时间同步。然而由于星地时间观测受大气层和地球表面环境影响，时间观测序列容易出现粗差，原子钟性能也可能出现相应异常扰动，需要对粗差进行分析处理。
目前已有一些学者对钟差异常值的处理方法进行了研究。其中对于卫星钟差异常值的探测主要有中位数法[1,2]、Allan方差法[3-5]以及基于平稳时间序列的χ2检验法[6-8]等方法。但这些方法还存在一些问题，中位数法虽然简单易行，但是对粗差的大小不敏感；而且目前的处理方法未注意到区分卫星钟差数据中的异常值因数据间断、粗差、钟跳、相位跳变和频率跳变等不同原因而有多种类型，因而对多种类型的异常值处理不当；另外，探测出卫星钟差数据中的异常值以后，往往对异常位置处的数据赋值为0，或者是通过插值拟合的方法进行处理，未给出有效的异常值估值方法。鉴于此，本文引入时间序列异常值探测理论，分析钟差序列中多种类型的异常情况。其实除了频率跳变以外，数据间断、粗差、相位跳变等形成的异常值在频率数据中都为附加类异常值(additive outlier，AO)，可以将上述多种类型的异常形式都归结为AO类异常值的探测问题。
鉴于星地时间同步钟差应用问题的实际情况，本文利用基于自回归(auto-regrossive，AR)模型的AO类异常值探测的Bayesian方法对其中的异常值进行处理[9]。首先，引入基于识别变量的AR模型异常值探测的Bayesian方法对星地时间同步钟差序列中的异常值进行探测；将异常值探测问题转化为识别变量后验概率值的计算问题，并给出明确的判别阈值；运用Gibbs抽样设计了识别变量后验概率值的计算方法，通过后验概率值与事先给定的阈值进行比较，判别出异常值的位置。其次，基于迭代似然比检验法中的异常值描述模型，将异常值估值问题转化为简单的线性模型最小二乘估计问题，设计了异常值的估值方法；并将该估值方法用于钟差异常的估值中，以期对钟差序列中的异常值进行修复。实测数据对该算法的验证表明，本文方法能够准确地探测出异常值的位置并精确估计出异常值的大小；进一步从异常值剔除前后的钟差预报精度比较可知，异常处理后的预报精度得到了提高。
-
星地双向时间比对方法得到的是卫星钟差数据(为区别于时间序列，常称为相位数据)，但在实际数据处理中更常用星钟的频率数据，即瞬时相对频率偏差[10-13]。这两类数据可以通过一定的变换关系进行转换。一般来说，相位数据数值太大不利于进行粗差异常检测，而频率数据是描述原子钟频率源内部参数的基本观测量，便于异常检测。鉴于此，本文将卫星钟差数据转化为频率数据，然后基于频率数据进行异常值探测。
在原始相位数据和频率数据中经常出现的异常值有数据间断、粗差、相位跳变和频率跳变。除了频率跳变以外，其它的异常形式在频率数据中都表现为孤立的多个或者单个异常值[14]，从时间序列异常值探测理论来看，这些孤立的异常值就是时间序列中的AO类异常值[14,15]。因此考虑到钟差时间序列数据的随机特性，本文以AR模型作为钟差预报模型，进一步引入基于AR模型的AO类异常值探测Bayesian方法，将数据间断、粗差、相位跳变等异常类型都归结为AO类异常值的探测问题，并设计如下的探测规则对AO类异常值进行探测。
-
设有一组星地频率时间序列数据{x1,…,xn}符合如下的AR(p)模型，p为AR模型：
(1) 式中，at为模型拟合后的噪声；i.i.d.的含义为相互独立同分布；Φ=(φ1,…,φp)T和σ2为未知参数；N(0,σ2)代表均值为0、方差为σ2的一元正态分布；at i.i.d. N(0,σ2)表示at对于t=1,…,n相互独立且都服从一元正态分布N(0,σ2)。
对每个观测值xt分别引入异常值识别变量：
(2) 记w1、…、wn分别代表每个时刻观测值的异常值大小，并假设：
(1) 每个观测值xt受到异常扰动的先验概率为α[14]，即P(δt=1)=α；
(2) 根据共轭先验分布的选取准则[16,17]和实际应用的需要，取参数的先验分布为：
其中，μ、ξ、α、Φ0、V、υ和λ为超参数； wt i.i.d. N(μ,ξ2)表示wt对于对于t=1,…,n相互独立且都服从一元正态分布N(μ,ξ2)；b(·,·)代表两点分布；Np(·,·)代表p元正态分布；IG(·,·)代表倒伽玛分布。
根据以上假设，观测值xt可表示为：
(3) 其中,zt代表未受到异常扰动的干净数据。由此，可得基于AR模型的异常值探测的模型：
(4) 不妨设前p个观测值x1、…、xp为正常值[18]，为判定观测值中是否含有异常值以及确定判别阈值，构造如下Bayesian假设检验问题。
原假设H0：xt是正常观测值即δt=0；备选假设H1：xt为异常值，即δt=1。
根据Bayesian假设检验原理[16]，当备选假设H1对应的后验概率P(δt=1|X)大于原假设H0对应的后验概率P(δt=0|X)，即P(δt=1|X)＞0.5时，认为备选假设成立，从而认为观测值xt为异常值；反之，认为观测值xt是正常观测值，其中X=(xp+1,…,xn)T。这样，问题归结为计算每个观测值xt含有异常值的后验概率P(δt=1|X)。
由于后验概率P(δt=1|X)涉及的分布比较复杂，本文引入Gibbs抽样算法[19]来解决这些后验概率值的计算问题，进而判断钟差序列中异常值的位置。根据Bayesian定理[20]，可得下列参数的完全条件分布[9]：
(5) 
(6) 
(7) 
(8) 其中，W=(w1,…,wn)T，
(9) 其中，
(10) 
有了以上的完全条件分布，就可以进行Gibbs抽样。假设φ(r)、(σ2)(r)、δ(r)、W(r)(r=1,…,R)为用Gibbs抽样算法从上述完全条件分布中抽取的样本，则由式(9)可得识别变量后验概率值的计算公式：
(11) 其中，
-
步骤1 确定先验分布中的超参数。
步骤2 根据Bayesian估计方法以及超参数的取值，确定Gibbs抽样的初值。
步骤3 假定第s≥1次抽样开始时样本值向量为(φ(s-1),(σ2)(s-1),δ(s-1),W(s-1))，则第s次抽样按下列方式产生样本值向量(φ(s-1),(σ2)(s),δ(s),W(s))。
φ(s)从条件分布p(φ|X,(σ2)(s-1),δ(s-1),W(s-1))中抽取；(σ2)(s)从条件分布p(σ2|X,φ(s),δ(s-1),W(s-1))中抽取；(δj)(s)从条件分布p(δj|X,φ(s),(σ2)(s),δ(-j)(s-1,j+1),W(s-1,j))中抽取；(wj)(s)从条件分布p(wj|X,φ(s),(σ2)(s),δ(s-1,j+1),W(-j)(s-1,j+1))中抽取。
其中，向量的上角标(i,k)的含义为该向量的第1个分量至第k-1个分量是第i+1次抽样时抽取的样本，第k个分量至最后一个分量为第i次抽样时抽取的样本，例如：
重复上述抽样过程R次，得到Gibbs样本：
步骤4 异常值定位。按照式(11)计算识别变量的后验概率值。
-
当频率序列数据存在AO类异常扰动时，只要观测值xt受到异常扰动，则可将其表示为[21]：
(12) 式中，zt为未受异常扰动的干净数据；f(t)代表异常扰动，且f(t)=w0Φ-1(B)ξ(d)t，w0代表异常扰动大小Φ-1(B)作为一个乘数因子，是由异常值的类型决定的，
现在假设已经探测出了异常值的位置，即d是已知量，为了估计异常扰动大小w0，对式(12)两边同乘以Φ(B)，则可得出如下形式：
(13) 令lt=Φ(B)xt，mt=ξt(d)，at=Φ(B)zt则式(13)等价于：
(14) 从式(14)可以看出，复杂的异常扰动估值问题转化为简单的线性模型最小二乘估计问题：
(15) 即
(16) 由此即可估计出异常值的大小。
-
为验证本文提出的异常值处理方法的正确性和可行性，利用北斗卫星导航系统的实测数据进行试验和分析。考察的应用模式为钟差预报的快速恢复，适用于全部卫星在出现星钟调频/调相后的快速恢复，同时也适用于倾斜地球同步轨道(inclined geosynchronous satellite orbit，IGSO)和中地球轨道(medium earth orbit，MEO)入境后的快速恢复。快速恢复可设计为采用最新获取的20 min数据拟合线性模型并进行2 h预报，由于可利用时段数据较少，钟差模型对异常数据更为敏感。
数据源为某地球静止轨道(geosyn-chronous Earth orbit,GEO)卫星的3 d星地时间比对的卫星钟差数据，连续截取了72组实验数据，每组数据均为2 h 20 min，其中前20 min数据用于拟合钟差模型参数，后2 h用于钟差预报的检验，采样间隔为6 s。
-
选取一组质量较好、无显著粗差的数据段，首先取其中前20 min拟合数据，采用如下2种方案进行试验和计算。
方案1 不加入粗差；
方案2 在预报阶段的第22个、第162个和第962个观测值上分别加大小为8 ns、-16 ns和9 ns的单个异常值。
对应的识别变量后验概率值如图1所示。从图1(a)可以看出，当数据中不含异常值时，识别变量的后验概率值小于0.5;从图1(b)可知，当数据中加入3个孤立的异常值时，相应位置的识别变量的后验概率值均大于0.5，由异常值的判别规则可知，第22个、162个和962个位置上均含有3个孤立的AO类异常值。
为了提高数据的预测精度，下面利用§2的异常值估值方法对异常值的扰动大小进行估值，以期对数据进行修复。
图 1 识别变量的后验概率值
Figure 1. Posterior Probabilities of Classification Variables
表1给出了第22个、第162个和第962个数据点上含有的异常扰动估值，可以根据估值大小对这三个异常值进行修复。
表 1 异常值的扰动估值大小
Table 1. Magnitudes of Outliers
异常位置 扰动估值/10-9s 22 8.047 171 462 238 824 162 -15.675 900 824 393 180 962 8.975 020 263 531 422 -
利用本文提出的Bayesian异常值探测和估值方法对上述72组实验数据进行处理，图2给出了异常处理前后72组数据的2 h的最大预报偏差。从图2可以看出，异常值处理后，最大预报偏差明显降低。为了将本文方法与以往方法的效果进行比较，利用抗差最小二乘估计，采用IGG3等价权因子公式，对72组数据分别进行抗差估计，其中参数的取值范围为k0取(1,1.5),k1取(2.5,8)。图3给出了基于抗差最小二乘的预报2h的偏差最大值。从图3可知，抗差最小二乘估计对异常值的处理效果不佳。
图 2 异常值剔除前后的最大预报偏差
Figure 2. Maximum Prediction Biases Before and After Outliers Detection
图 3 2 h预报钟差最大值统计
Figure 3. Maximum Prediction Clock Offsets Based on the Observations of 2 h
-
本文提出了基于AR模型的星地钟差序列的异常值探测方法。该方法针对AR模型引入识别变量，通过比较这些识别变量的后验概率值与事先给定的阈值来进行异常值定位。 基于迭代似然比检验的异常值描述模型，设计了钟差大小的估值方法。该方法将异常值估值问题转化为了简单的线性最小二乘计算问题。
基于AR模型的星地钟差序列的异常值探测方法将多种类型的异常情况都转化为了AO类异常值的探测问题，适合于处理数据间断、粗差、相位跳变等多种类型的异常形式。 通过实测数据表明,本文的方法能够探测出钟差序列中的异常值而且能够以较高的精度估计出异常值的大小。
Bayesian Methods for Outliers Detection and Estimation in Clock Offset Measurements of Satellite-Ground Time Transfer
-
摘要: 由于星地时间观测受大气层和地球表面环境影响,时间观测序列容易出现粗差,原子钟性能也可能出现相应异常扰动,需要对粗差进行分析处理。对此,本文引入基于识别变量的自回归(auto-regressive,AR)模型异常值探测的Bayesian方法对星地时间同步钟差序列中的异常值进行探测,进一步基于迭代似然比检验法中的异常值描述模型,将异常值估值问题转化为简单的线性模型最小二乘估计问题,以期对钟差序列中的异常值进行修复。实验表明本文的方法能够准确的探测出异常值的位置并精确的估计出异常值的大小。Abstract: Clock offset measurements of satellite-ground time transfer are usually affected by outliers due to the impact of ionosphere errors, tropospheric errors, and multipath effects. Therefore, in this paper, we propose an autoregressive model based on Bayesian methods for detecting outliers in the clock offset measurements with the classification variables. Furthermore, the model for estimating the magnitude of outliers is given to correct the clock offset measurements, and solves the problem of outlier estimation by transforming it into a simple least square problem. Different schemes based on the real BDS data were designed to evaluate the performance of the new Bayesian method. We applied the new method ito the fast recovery of the clock offset prediction. Test examples illustrate that the Bayesian methods can detect the outliers effectively and estimate the magnitudes of outliers accurately.
-
Key words:
- satellite-ground clock offset /
- time transfer /
- outliers /
- auto-regressive model
-
图 2 异常值剔除前后的最大预报偏差
Figure 2. Maximum Prediction Biases Before and After Outliers Detection
图 3 2 h预报钟差最大值统计
Figure 3. Maximum Prediction Clock Offsets Based on the Observations of 2 h
表 1 异常值的扰动估值大小
Table 1. Magnitudes of Outliers
异常位置 扰动估值/10-9s 22 8.047 171 462 238 824 162 -15.675 900 824 393 180 962 8.975 020 263 531 422 
下载: 导出CSV
-
[1] Winkler G M R. Introduction to Robust Statistics and Data Filtering[J/OL]. Hamilton Technical Services, http://www.wriley.com/ROBSTAT.htm, 1993 [2] 郭海荣.导航卫星原子钟时频特性分析理论与方法研究[D].郑州:信息工程大学,2006 Guo Hairong. Study on the Analysis Theories and Algorithms of the Time and Frequeney Charaeterization for Atomie Clocks of Navigation Satellites[D]. Zhengzhou:Information Engineering University, 2006 [3] 杨剑,王泽民,王贵文,等.GPS接收机钟跳的研究[J].大地测量与地球动力学,2007,27(3):123-127 Yang Jian, Wang Zemin, Wang Guiwen, et al. Research on Clock Jumps of GPS Receiver[J].Journal of Geodesy and Geodynamics, 2007, 27(3):123-127. [4] 牛飞,韩春好,张义生,等.导航卫星星载原子钟异常监测分析[J].武汉大学学报·信息科学版,2009,34(5):585-588 Niu Fei, Han Chunhao, Zhang Yisheng, et al. Analysis and Detection on Atomic Clock Anomaly of Navigation Satellites[J]. Geomatics and Information Science of Wuhan University, 2009, 34(5):585-588 [5] Riley W J. Frequency Jump Detection and Analysis[C]. The 40th Annual Precise Time and Time Interval Meeting, Virginia, 2008 [6] 黄新明,龚航,朱祥维,等.基于Kalman滤波器的卫星钟实时异常监测算法[J].宇航计测技术,2011,31(5):6-11 Huang Xinming, Gong Hang, Zhu Xiangwei, et al. A Real Time Anomaly Monitoring Algorithm for Satellite Clock Based on Kalman Filter[J]. Journal of Astronautic Metrology and Measurement, 2011, 31(5):6-11 [7] 黄观文.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 [8] Riley W J. Outliers in Time and Frequency Measurements[R/OL]. Hamilton Technical Services, http://www.stable32.com/Outliers%20in%20Time%20and%20Frequency%20Measurements.pdf, 2013 [9] 张倩倩,归庆明,王延停.基于不同类型识别变量的AR模型异常值探测的Bayes方法[J].测绘学报,2012,41(3):378-384 Zhang Qianqian, Gui Qingming, Wang Yanting. Bayesian Methods for Outliers Detection in Autoregressive Model Based on Different Types of Classification Variables[J]. Acta Geodaeticaetica et Cartographica Sinica, 2012, 41(3):378-384 [10] 黄飞江.基于星间距离变化的动态双向时间同步算法[J].武汉大学学报·信息科学版,2010,35(1):13-16 Huang Feijiang. An Algorithm of Dynamic Two-Way Time Transfer Based on Inter-satellite Range Variation[J]. Geomatics and Information Science of Wuhan University, 2010, 35(1):13-16 [11] Koppang P, Wheeler P. Working Application of TWSTT for High Precision Remote Synchronization[C]. The Annual IEEE International Frequency Control Symposium, CA, USA, 1998 [12] 钟兴旺.卫星运动对星间双向法时间同步的影响分析与校正[J].中国空间科学技术,2007,27(6):54-58 Zhong Xingwang. Analysis and Correction Techniques of Movement Influence on Inter-satellite Two-Way Time Transfer[J]. Journal of Chinese Space and Technology, 2007,27(6):54-58 [13] 徐冬梅.利用条件平差实现导航星座自主守时方法[J].中国空间科学技术,2006,26(6):18-23 Xu Dongmei. Auto-Timing of Satellite Navigation System Using Conditional Estimation[J].Journal of Chinese Space and Technology, 2006, 26(6):18-23 [14] Abraham B, Box G E P. Bayesian Analysis of Some Outlier Problems in Time Series[J]. Biometrika, 1979, 66:229-236 [15] 韦博成.统计诊断引论[M].南京:东南大学出版社,1991 Wei Bucheng. The Introduction of Statistical Diagnosis[M]. Nanjing:Press of Southeast University, 1991 [16] Koch K R. Bayesian Inference with Geodetic Applications[M]. Berlin:Springer, 1990 [17] Gui Q, Li X, Gong Y. A Bayesian Unmasking Method for Locating Multiple Gross Errors Based on Posterior Probabilities of Classification Variables[J].Journal of Geodesy, 2011, 85:191-203 [18] McCulloch R E, Tsay R S. Bayesian Analysis of Autoregressive Time Series Via the Gibbs Sample[J]. Journal of Time Series Analysis, 1992, 2(15):235-250 [19] Christian P R. Monte Carlo Statistical Methods[M]. Berlin:Springer, 2004 [20] Box G E P, Tiao G C. Bayesian Inference in Statistical Analysis[M]. Massachusetts:Addison Wesley, 1973 [21] Tsay R S. Outliers, Level Shifts, and Variance Changes in Time Series[J]. Journal of Forecasting, 1988, 7:1-20 -
点击查看大图
图(3) / 表(1)
计量
- 文章访问数: 1098
- HTML全文浏览量: 41
- PDF下载量: 309
- 被引次数: 0
