有效可靠的卫星钟差数据是进行星载原子钟性能分析和卫星钟差建模与预报的前提和基础^[1]，也是卫星导航定位系统实现精准定位的关键。然而，由于原子钟本身性能、空间环境的影响和各种不确定因素的干扰以及有计划的卫星钟调相调频操作，所获取的卫星钟差数据中经常会出现异常扰动，从而降低了卫星钟性能分析的可靠性，破坏了钟差建模和预报的有效性，影响了导航定位结果的精准度。因此，需要对这些异常值进行处理。

近年来，已有不少学者对钟差异常值探测的方法开展了研究^[1-19]。目前常见的钟差异常值探测方法大致可以分为两类：一是通过比对直接探测观测数据中的异常值，如中位数法^[1-6]、Allan方差法及其扩展方法^[3-11]；二是针对观测数据首先建立相应模型，然后依托模型进行异常值探测，如基于二次多项式模型的递推遗忘因子最小二乘法^[12-13]、基于时间序列自回归（autoregressive, AR）模型的Bayes方法等^[14-19]。这些方法能够在一定程度上实现对异常值的探测，但也存在不足。直接法虽然简单方便，但对异常值的大小不敏感，易出现漏判现象；基于各种模型进行异常值探测的方法已逐渐为广大学者所认可，但基本上都是基于平稳序列进行研究的，且大部分方法都需要通过数值迭代来实现，不可避免地会受到迭代发散问题的困扰；一些方法还需要进行多重假设检验以及阈值确定等。这些问题的存在都增加了异常值探测的不确定性。

卫星钟差数据序列是一个非平稳时间序列，在长时间段的卫星钟差数据中存在着较为显著的数据异常^[1]，这些异常数据主要包括相位跳变、频率跳变、数据间断和数据粗差等^[1-6]。依据时间序列分析的观点，除频率跳变外，数据间断、数据粗差、相位跳变等形成的异常值都可归结为加性异常值（additive outlier，AO）^{[5, 15, 20-21]}。基于以上分析，本文基于求和自回归移动平均（autoregressive integrated move average, ARIMA）模型^[20-21]建立异常值探测模型，通过模型选择的方法对卫星钟差数据中的AO类异常值进行探测与处理。假定卫星钟差观测序列中含有多个异常值，针对异常值的不同情况建立相应的异常值探测模型；构造用于异常值探测的候选模型集合，从而将异常值的定位、定值问题转化为模型选择问题；基于Bayes统计原理确定各个候选模型的后验概率，导出模型选择的度量指标；选择指标值最大的候选模型为期望模型，以实现异常值的定位和定值。通过GPS和北斗导航卫星系统(BeiDou navigation satellite system, BDS)不同卫星钟差异常值的探测及钟差预报等试验，表明基于ARIMA模型的卫星钟差异常值探测的模型选择方法能同时解决异常值的定位和异常扰动的定值问题，验证了方法的正确性和有效性。

1 异常值探测模型

若有卫星钟差序列{x_t}符合如下形式，则称其满足ARIMA(p, d, q)模型^[20-21]：

$ \left\{ {\begin{array}{*{20}{l}} {\phi \left( B \right){\nabla ^d}{x_t} = \mathit{\Theta }\left( B \right){\varepsilon _t}}\\ {{\varepsilon _t}{\rm{i}}{\rm{.i}}{\rm{.d}}.N\left( {0, {\sigma ^2}} \right)} \end{array}} \right. $

(1)

式中，B为后移算子，如${B^k}{x_t} = {x_{t - k}}$；$\phi \left( B \right) = 1 - {\phi _1}B - {\phi _2}{B^2} - \cdots - {\phi _p}{B^p}$和$\mathit{\Theta }\left( B \right) = 1 - {\theta _1}B - {\theta _2}{B^2} - \cdots - {\theta _q}{B^q}$分别为B的p阶和q阶多项式，它们的零点都在单位圆外；$\nabla = 1 - B$为差分算子；d为差分次数；经过d次差分运算后，${\nabla ^d}{x_t}$为平稳时间序列；${{\varepsilon _t}{\rm{i}}{\rm{.i}}{\rm{.d}}.N\left( {0, {\sigma ^2}} \right)}$表示${\varepsilon _t}\left( {t = 1, 2 \cdots } \right)$为相互独立且同分布于均值为零、方差为${{\sigma ^2}}$的正态分布。

对于一组卫星钟差观测数据${y_1}, {y_2} \cdots {y_T}$，若其中含有m个AO类异常值，建立异常值探测模型${M_{{\tau _m}}}$：

$ \left\{ {\begin{array}{*{20}{l}} {\begin{array}{*{20}{l}} {{y_t} = {x_t} + {w_{{t_1}}}I_t^{\left( {{t_1}} \right)} + \cdots + {w_{{t_m}}}I_t^{\left( {{t_m}} \right)}}\\ {\varphi \left( B \right){\nabla ^d}{x_t} = \mathit{\Theta }\left( B \right){\varepsilon _t}} \end{array}}\\ {{\varepsilon _t}{\rm{i}}{\rm{.i}}{\rm{.d}}{\rm{.}}N\left( {0, {\sigma ^2}} \right)} \end{array}} \right. $

(2)

式中，x_t为不含任何异常值的纯净时间序列，且满足式（1）；${\mathit{\boldsymbol{\tau }}_m} = {[{t_1}\;{t_2} \cdots {t_m}]^{\rm{T}}}$为异常值位置指示向量，即对于${\mathit{\boldsymbol{\tau }}_m}$的任意一个分量h，如果t=h，则$I_t^{\left( h \right)} = 1$，否则$I_t^{\left( h \right)} = 0$；${w_{{t_1}}}, {w_{{t_2}}} \cdots {w_{{t_m}}}$表示相应位置处异常扰动的大小；${\mathit{\boldsymbol{\rho }}_{{\tau _m}}} = \left( {{\varphi _1}\;{\varphi _2} \cdots {\varphi _p}, {\theta _1}\;{\theta _2} \cdots {\theta _q}, {w_{{t_1}}}\;{w_{{t_2}}} \cdots {w_{{t_m}}}, \sigma } \right)$为模型参数；p_m = p + q + m + 1为参数个数。异常值探测就是找出观测序列中异常值的位置并计算其扰动大小，即确定${\mathit{\boldsymbol{\tau }}_m} = {[{t_1}\;{t_2} \cdots {t_m}]^{\rm{T}}}$和${w_{{t_1}}}, {w_{{t_2}}} \cdots {w_{{t_m}}}$。

2 异常值探测原理

常见的基于模型的异常值探测方法均需要通过多次迭代来实现。如基于AR模型的Bayes方法在计算时不得不考虑迭代的收敛问题，有时还需进行多重假设检验；基于二次多项式模型的方法存在临界值选取困难等问题。在此借鉴模型选择的思想^[22-23]，从整体上考虑异常值数量m、位置${\mathit{\boldsymbol{\tau }}_m} = {[{t_1}\;{t_2} \cdots {t_m}]^{\rm{T}}}$和扰动大小${w_{{t_1}}}, {w_{{t_2}}} \cdots {w_{{t_m}}}$同时确定的问题。

首先，设用于异常值探测的候选模型集合为$M = \left\{ {{M_{{\tau _0}}}, {M_{{\tau _1}}} \cdots {M_{{\tau _m}}}} \right\}$，其中${M_{{\tau _i}}}\left( {i = 0, 1, 2 \cdots m} \right)$表示观测序列中含有i个异常值的探测模型，m表示观测序列中可能存在异常值的最大个数。因此，候选模型的总数为：

$ S = \left( {\begin{array}{*{20}{c}} T\\ 0 \end{array}} \right) + \left( {\begin{array}{*{20}{c}} T\\ 1 \end{array}} \right) + \cdots + \left( {\begin{array}{*{20}{c}} T\\ m \end{array}} \right) $

(3)

式中，$\left( {\begin{array}{*{20}{c}} T\\ i \end{array}} \right)$表示从T个观测值中选择i个异常值的组合数，${i = 0, 1, 2 \cdots m}$。

其次，根据Bayes统计推断原理，计算含有i个异常值的候选模型的后验概率：

$ P({M_{{\tau _i}}}|\mathit{\boldsymbol{y}}) = \frac{{P\left( {{M_{{\tau _i}}}} \right)L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})}}{{f\left( \mathit{\boldsymbol{y}} \right)}} $

(4)

式中，${i = 0, 1, 2 \cdots m}$；${P\left( {{M_{{\tau _i}}}} \right)}$是候选模型${{M_{{\tau _i}}}}$的先验分布；${L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})}$为在给定候选模型${{M_{{\tau _i}}}}$条件下的边缘似然函数；$\mathit{\boldsymbol{y}} = {[{y_1}\;{y_2} \cdots {y_T}]^{\rm{T}}}$；$f\left( \mathit{\boldsymbol{y}} \right) = \sum\limits_{j = 0}^m {\sum\limits_{{\tau _j}} {P\left( {{M_{{\tau _j}}}} \right)L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})} } $是y的无条件似然函数。

最后，根据式(4)的结果选择后验概率最大的模型作为异常值探测模型；当模型确定后，异常值的个数、位置和异常扰动的大小同时确定。由此，将异常值的探测问题转化为模型的选择问题。

3 后验概率的计算方法

由式(4)可以看出，后验概率的计算主要是确定模型的先验分布$P\left( {{M_{{\tau _i}}}} \right)$和给定模型条件下的边缘似然函数$L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})$。在模型未确定之前，每个候选模型的地位相当，因此采用均匀分布作为模型的先验分布，即$P\left( {{M_{{\tau _i}}}} \right) = \frac{1}{m}\frac{1}{c}$，其中$c = \left( {\begin{array}{*{20}{c}} T\\ i \end{array}} \right)$。边缘似然函数$L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})$可以表示为：

$ L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}}) = \mathop \smallint \nolimits^ L({\mathit{\boldsymbol{\rho }}_{{\tau _i}}}|\mathit{\boldsymbol{y}}, {M_{{\tau _i}}})p({\mathit{\boldsymbol{\rho }}_{{\tau _i}}}|{M_{{\tau _i}}}){\rm{d}}{\mathit{\boldsymbol{\rho }}_{{\tau _i}}} $

(5)

式中，$L({\mathit{\boldsymbol{\rho }}_{{\tau _i}}}|\mathit{\boldsymbol{y}}, {M_{{\tau _i}}})$为观测序列的似然函数；$p({\mathit{\boldsymbol{\rho }}_{{\tau _i}}}|{M_{{\tau _i}}})$为给定模型${M_{{\tau _i}}}$条件下参数${\mathit{\boldsymbol{\rho }}_{{\tau _i}}}$的条件先验概率。鉴于边缘似然函数$L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}})$的解析表达式难以通过直接计算得到，故采取近似计算。对似然函数$L({\mathit{\boldsymbol{\rho }}_{{\tau _i}}}|\mathit{\boldsymbol{y}}, {M_{{\tau _i}}})$取对数，令：

在参数${{\mathit{\boldsymbol{\rho }}_{{\tau _i}}}}$的精确最大似然估计值${{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}$处，将对数似然函数${l_{{\tau _i}}}\left( {{\mathit{\boldsymbol{\rho }}_{{\tau _i}}}} \right)$进行二阶展开，从而导出式（5）的Laplace近似公式^[22]，故式（5）又可以表示如下：

$ L(\mathit{\boldsymbol{y}}|{M_{{\tau _i}}}) = {\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{T}} \right)^{\frac{{{p_i}}}{2}}}{\rm{exp}}\left( {{l_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right)} \right)p({{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}|{M_{{\tau _i}}})|{\mathit{\boldsymbol{H}}_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right){|^{ - \frac{1}{2}}} + O\left( {{T^{ - \frac{{{p_i}}}{2} - 1}}} \right) $

(6)

式中，p_i表示参数个数；$p({{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}|{M_{{\tau _i}}})$表示参数${{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}$的条件先验概率；${\mathit{\boldsymbol{H}}_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right)$表示${T^{ - 1}}{l_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right)$的Hessian矩阵；O(•)表示高阶无穷小。式（4）又可表示如下：

$ P({M_{{\tau _i}}}|\mathit{\boldsymbol{y}}) = \frac{{P\left( {{M_{{\tau _i}}}} \right)}}{{f\left( \mathit{\boldsymbol{y}} \right)}}{\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{T}} \right)^{\frac{{{p_i}}}{2}}}{\rm{exp}}\left( {{l_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right)} \right)p({{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}|{M_{{\tau _i}}})|{\mathit{\boldsymbol{H}}_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right){|^{ - \frac{1}{2}}} + O\left( {{T^{ - \frac{{{p_i}}}{2} - 1}}} \right) $

(7)

从理论研究的角度来看，采用后验概率$P({M_{{\tau _i}}}|\mathit{\boldsymbol{y}})$的数值最大的模型作为最终选择的模型是合适的；在实际应用中，顾及计算时间和计算量等问题，对后验概率$P({M_{{\tau _i}}}|\mathit{\boldsymbol{y}})$的形式进行进一步简化。对式(7)两端分别取对数，得到如下公式：

$ {\rm{ln}}(P({M_{{\tau _i}}}|\mathit{\boldsymbol{y}})) = {l_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right) + \frac{{{p_i}}}{2}{\rm{ln}}\left( {\frac{{2{\rm{ \mathsf{ π} }}}}{T}} \right) + {\rm{ln}}(p({{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}|{M_{{\tau _i}}})) - \frac{1}{2}{\rm{ln}}\left( {\left| {{\mathit{\boldsymbol{H}}_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right)} \right|} \right) - {\rm{ln}}\left( {f\left( \mathit{\boldsymbol{y}} \right)} \right) + O\left( {{T^{ - 1}}} \right) $

(8)

由贝叶斯信息准则的推导过程可知，式（8）的主要部分为其前两项，其中第1项反映的是模型拟合的好坏，第2项反映的是模型参数的多少^[22-23]。考虑到异常扰动项的加入会使参数增加，为限制参数过度化产生的影响，引入惩罚因子${\rm{ln}}\left( {\begin{array}{*{20}{c}} m\\ i \end{array}} \right)$，同时为了使表达式的系数为整数，式（8）的两端同时乘以2，得到可用于模型选择的度量指标（measure of detecting outliers, MDO）：

$ {\rm{MDO}} = 2{l_{{\tau _i}}}\left( {{{\mathit{\boldsymbol{\hat \rho }}}_{{\tau _i}}}} \right) - {p_i}{\rm{ln}}\left( T \right) - 2{\rm{ln}}\left( {\begin{array}{*{20}{c}} m\\ i \end{array}} \right) $

(9)

使MDO值最大的候选模型即为期望模型。

由式（3）易见，随着数据量T的增加，候选模型总数S将会迅速增长，使模型选择变得比较复杂。为降低计算复杂度，可采用两阶段法实施计算。第1阶段即假设数据中只有1个异常值，用逐步型法筛选出一部分MDO值较大的观测值标记为疑似异常值，通常疑似异常值数量小于数据总数的5%；第2阶段是在已选出的疑似异常值的基础上，再进一步计算其组合的MDO值，从而实现快速异常值探测。

4 算例与分析

为了验证本文提出的基于模型选择的卫星钟差异常值探测法的正确性和有效性，采用国际GNSS服务官网发布的GPS卫星和BDS卫星精密钟差数据进行试验和分析。

4.1 方案1

提取G31卫星钟为期30 d的钟差数据，时间为2017-09-16 00:00至2017-10-15 23:55，采样间隔为15 min，共计2 880个历元。从中截取一段为期3 d且质量好、无显著粗差的数据，其中前88个数据用于拟合钟差模型，后200个数据用于进行异常值探测试验。在后200个历元的第31、121和187个观测值上分别加上大小为3 μs、-3 μs和6 μs的单个异常扰动。采用本文提出的模型选择方法进行异常值探测，各种情形下相应的MDO值列于表 1。其中，第2列（i=0）最后一行代表数据中不含异常值时的MDO值；第3~5列（i=1）最后一行代表数据中含有1个异常值时最大的3个MDO值（异常值位于其他位置时，MDO值均小于-3.2），相应的τ_i为3个异常值出现的位置，w_ti为对应异常扰动值大小；第6~8列（i=2）最后一行表示数据中含有2个异常值时最大的3个MDO值（异常值位于其他位置时，MDO值均小于-2.8），相应的τ_i为3对异常值出现的位置，w_ti为对应位置处异常扰动值大小；第9列（i=3）最后一行表示数据中含有3个异常值时的最大MDO值（异常值位于其他位置时，MDO值均小于20）。

从表 1可以看出：（1）该钟差序列的第31、121和187个观测值为异常值，异常扰动大小分别为2.974 5 μs、-2.826 3 μs和6.040 1 μs，与所加的异常扰动比较接近。（2）对于固定的异常值数量i，候选模型的MDO值受到异常值大小和位置的影响。（3）对于不同的i，MDO值既与异常值的数量、大小和位置有关，也和惩罚因子有关。（4）最终探测模型的MDO值明显大于其他模型，说明所提方法能较为准确地探测到异常值的位置，并确定异常扰动的大小。

4.2 方案2

基于上述30 d的钟差数据进行异常值探测和钟差预报试验分析。首先对前100个历元的试验数据进行异常值探测并修正钟差序列预报模型，然后对后续2 780个历元的钟差进行预报。试验结果分析如下:（1）使用本文提出的基于模型选择的异常值探测方法对前100个数据进行探测时，发现其中存在4处异常，分别位于第5、6、7和38个历元，异常扰动分别为130 ns、120 ns、300 ns和-100 ns。采用基于AR模型的Bayes方法^[14-16]只能探测到第7个历元和第38个历元处的数据受到异常影响，异常扰动分别为180 ns和99 ns。这表明基于ARIMA模型的卫星钟差异常值探测的模型选择方法对于成片异常值的探测也有一定作用。（2）将数据修正后应用于后续钟差预报，并将预报结果与目前广泛使用的抗差二次多项式的预报结果和新近发展的基于AR模型的Bayes方法的预报结果分别进行比较，比较的指标为预测均方根误差（root mean square error of prediction, RMSEP）、平均误差（Mean）和最大偏差绝对值（maximum absolute biases，MAB），结果列于表 2。比较发现，在这3个指标下，基于AR模型的Bayes方法的预报精度优于抗差二次多项式模型的预报精度，而基于ARIMA模型的模型选择方法的预报精度又优于基于AR模型的Bayes方法的预报精度。

4.3 方案3

为进一步验证所提方法的有效性，采用北斗不同类型卫星的精密钟差数据进行试验分析。选用北斗地球同步轨道（geostationary earth orbit，GEO）卫星C02、倾斜地球同步轨道（inclined geosynchronous satellite orbit，IGSO）卫星C10和中地球轨道（medium earth orbit，MEO）卫星C14的钟差数据，时间从2018-09-09 00:00至2018-09-10 23:55，采样间隔为5 min。在3颗卫星的第一天精密钟差数据中分别随机抽取3个历元的数据加入粗差，以此为背景建立模型，对3颗卫星后一天的钟差分别进行预报。预报结果与实测值的误差分布见图 1。可以看出，本文提出的基于ARIMA模型的模型选择方法的预报精度比基于AR模型的Bayes方法的预报精度和采用抗差二次多项式模型的预报精度都高。

5 结语

卫星钟差异常值探测不仅有助于提高定位和时频传递精度，也对卫星完好性监测提供了有力支持。因此，寻求一种高效的钟差序列异常值探测方法显得尤为重要。本文构建了一种异常值探测新模型，提出了基于ARIMA模型的钟差异常值探测的模型选择法；以Bayes统计原理为理论基础，导出模型的后验分布，将异常值的探测问题转化为候选模型的选择问题；推导出可用于模型选择的度量指标MDO，解决了后验概率计算复杂的实际问题；针对数据量T的增加会使候选模型总数S迅速增长，从而导致模型选择变得比较复杂的新情况，提出了两阶段法实施计算，以降低计算复杂度；通过GPS和BDS卫星钟差数据异常值探测和钟差预报等仿真试验分析，验证了所提方法的正确性和有效性。