-
星地双向时间同步[1-2]是目前精度最高的时间同步方法之一，可计算高精度的星地钟差，然而其直接计算得到的卫星钟差未经平滑，受到卫星钟、卫星收发设备、地面钟、地面收发设备、信号路径以及解算软件等众多因素影响，可能出现粗差或钟跳。针对卫星钟出现的异常可以通过监测的方式及时告警[3]，减少卫星钟差异常值数量或对异常情况进行记录。但产生异常值的原因很多，仅在物理层面对卫星钟的监测难以杜绝卫星钟差异常的产生。在卫星数据处理领域，针对粗差未完全剔除的情况，常采用抗差估计以减弱异常值对数据处理结果的影响，以期得到一些不是很好但仍然可用的观测值[4-5]，但抗差估计方法难以完全去除异常数据的影响。
星地双向时间同步方法常常用于星地时间同步试验中，得到的钟差数据未经后处理，卫星钟差观测值可能出现缺失数据段，同时有效数据段内各类异常也可能同时出现。当钟跳和粗差两类异常同时段出现时，将对对方的探测产生不利影响，必须同时解决钟跳和粗差的探测问题；在卫星测试阶段可能出现粗差量级远大于钟差数据本身的异常，将使最小二乘残差算法的异常值处理方法失效，钟跳和粗差难以分离。目前同时解决钟跳和粗差探测的研究较少。文献[6]针对GPS卫星钟差进行钟跳判定和粗差剔除，采用对钟差按历元单差，以单差序列的10倍中位数参数为阈值，利用单个粗差使单差序列相应历元形成正负两个粗大值，而单个钟跳使单差序列相应历元形成一个粗大值的特点，同时探测钟跳和剔除粗差，但该方法未考虑实际情况下卫星钟差多个粗差连续出现以及钟跳邻近历元出现粗差的情况，这些情况下钟跳和粗差无法被有效区分，且该算法难以探测小粗差。文献[7]利用频率序列识别钟跳，采用阈值法进行异常数据段剔除以及结合绝对中位数法(median absolute deviation，MAD)和Baarda数据探测法共同进行粗差探测与剔除，对北斗系统多种钟差异常同时出现的情况进行有效控制，但该方法采用恒定的经验阈值(文中取1×10-5 s)判定钟跳，对小于经验阈值的异常钟跳无法判定，同时也将使包含这些钟跳的数据段内的粗差探测失效。
本文对卫星钟差数据异常特点进行分析，顾及各类异常情况，据此探讨历元单差算法，使用中位数或中误差动态计算阈值，结合拟准检定法做进一步粗差探测，达到精准探测钟跳和粗差的目的，使用钟差模拟数据进行算法校验，利用卫星实测数据进行验证，分析了在复杂情况下剔除小粗差和探测钟跳的能力，另外还在算例中给出了探测成功后的进一步处理示例。
HTML
-
粗差是相对于观测精度而言的，在经典的最小二乘平差中，通常认为大于3倍中误差的真误差为粗差[8]。常使用3倍中误差或残差的绝对中位数作为阈值。
3倍中误差准则，即3σ准则，是最常用也是最简单的判别粗差的准则，其以测量次数充分大为前提，在实际测量中是一个近似的准则[9]。3倍中误差准则计算全样本的单位权中误差，当粗差较大时, 单位权中误差数值会相应较大，需多次应用3倍中误差准则才能有效消除粗差。
残差的绝对中位数法在卫星数据处理中应用广泛，具有较强的抗差性，其判定准则为[10-11]：
式中，V为残差序列；Vj为第j次测量的残差；σMAD为中位数误差；n为放大因子，一般取3~5，有时为降低误判概率，可取10；函数median·表示求取中位数。若式(1)成立，则认为该次观测含有粗差，应予剔除。残差的绝对中位数受粗差的影响很小，一般使用1次就可有效消除粗差。
-
历元单差公式表示为：
式中，测量值序列记作L；相应的时间序列记作T；历元单差结果记作D；Lj、Tj和Dj分别表示各序列第j个元素。卫星钟差历元单差值Dj物理意义上可看成卫星钟速，在序列L没有异常值的情况下，D序列变化极小。
依据§1.1方法设置阈值，若某个历元单差值超过阈值，则认为测量值相应位置发生了值的突变，可判定发生了粗差或者钟跳。
典型的卫星钟差异常情况包括5种：单粗差情况、一般连续粗差情况、相近连续粗差值情况、钟跳情况以及钟跳相邻历元出现粗差的情况，分别如表 1-5所示。表中采用无量纲的数进行模拟分析，对历元单差真误差设置阈值0.1，真误差绝对值大于0.1的判定为异常。
项目 理论值 5 6 7 8 9 10 11 增量 0 0 0.5 0 0 0 0 观测值 5 6 7.5 8 9 10 11 历元单差 - 1 1.5 0.5 1 1 1 单差真误差 - 0 0.5 -0.5 0 0 0 Table 1. Single Gross Error
项目 理论值 5 6 7 8 9 10 11 增量 0 0 0.5 -0.3 0 0 0 观测值 5 6 7.5 7.7 9 10 11 历元单差 - 1 1.5 0.2 1.3 1 1 单差真误差 - 0 0.5 -0.8 0.3 0 0 Table 2. General Continuous Gross Errors
项目 理论值 5 6 7 8 9 10 11 增量 0 0 0.5 0.55 0.52 0 0 观测值 5 6 7.5 8.55 9.52 10 11 历元单差 - 1 1.5 1.05 0.97 0.48 1 单差真误差 - 0 0.5 0.05 -0.03 -0.52 0 Table 3. Continuous Gross Errors with Similar Value
项目 理论值 5 6 7 8 9 10 11 增量 0 0 1 1 1 1 1 观测值 5 6 8 9 10 11 12 历元单差 - 1 2 1 1 1 1 单差真误差 - 0 1 0 0 0 0 Table 4. Clock Jump
项目 理论值 5 6 7 8 9 10 11 增量 0 0 1 0.6 1 1 1 观测值 5 6 8 8.6 10 11 12 历元单差 - 1 2 0.6 1.4 1 1 单差真误差 - 0 1 -0.4 0.4 0 0 Table 5. Clock Jump and Gross Errors Happened in Adjacent Epochs
如表 1所示，若第N历元为单粗差情况，则第N-1和第N历元单差真误差判定为异常。
如表 2所示，若第N和第N+1历元出现连续粗差且粗差值相差较大，则第N-1至N+1这3个历元单差真误差判定为异常。
如表 3所示，第N至N+2历元出现连续粗差且粗差值相近，则第N-1和N+2这两个历元单差真误差判定为异常。
如表 4所示，第N个历元开始出现钟跳，则第N-1个历元单差真误差判定为异常。
如表 5所示，第N个历元出现钟跳，并且第N+1个历元出现粗差，则第N-1至N+1这3个历元单差真误差判定为异常。
探测粗差和钟跳应考虑这5种典型情况特点。利用历元单差法探测粗差和钟跳算法分为3步：
1) 历元单差序列阈值计算。利用原始观测值计算历元单差序列，利用§1.1阈值设定方法，以历元单差序列为观测量求取阈值。
2) 粗差剔除。若连续历元的单差值Dj和Dj+1均超过阈值，则判定历元Tj+1对应的原始观测值含粗差，予以剔除。
3) 疑似钟跳点判定。进行了第2)步粗差剔除之后，对剩余的数据计算历元单差序列，依第1)步所得的阈值Dth，若单差值Dj超过此阈值，则判定历元Tj和Tj+1可能发生钟跳，记疑似钟跳历元为Tj+1。
上述历元单差法探测粗差和钟跳算法成功剔除了表 1、表 2、表 5形式的粗差，准确判定表 4钟跳位置，能够判定表 5钟跳，但也存在缺陷：表 3粗差异常被判为疑似钟跳，表 5钟跳初始的正常观测值被剔除。值得注意的是，在有大量正常观测样本的情况下，有个别正常值被误剔除并不影响系统观测特性；相近连续粗差的出现一般有其系统性原因，对其作为疑似钟跳标记位置是有益的。
历元单差法能有效剔除粗差，其判定的疑似钟跳包含了真实钟跳情况。若在疑似钟跳观测值处对观测数据进行分段并做进一步处理，则可有效规避钟跳影响。
一般依次考虑以下3个因素对数据进行分段：①由历元单差法判定的疑似钟跳点；②数据中断点；③根据前两个因素分段后，若某段数据量过大，计算机难以处理，可对该段进行数据分割。应注意并非所有的数据中断都需要进行数据分段，只有当中断达到一定的时长时才进行数据分段，这个时长通常由经验给出，本文后面将介绍的实测数据算例中，1 s采样间隔的数据当中断时长达100 s时即进行分段。
-
对已经分段的钟差序列要进行更精准的粗差剔除，可分别对每段采用拟准检定法。
拟准检定法认为实际观测值大部分正确，粗差是少数。则可以选定r个拟准观测，此处直接给出真误差估值如下：
式中，${{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_r}$是拟准观测的真误差估值；${{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_1}$是非拟准观测的真误差估值；R表示平差因子阵；L是观测值；GQ=(0 ATr)取决于真误差估值对应的系数阵Ar；${\mathit{\boldsymbol{P}}_Q} = \left( {\begin{array}{*{20}{c}} 0&0 \\ 0&{{\mathit{\boldsymbol{P}}_r}} \end{array}} \right)$取决于真误差估值对应的权阵Pr。
拟准观测选取正确时，式(3)计算的真误差估值出现明显分群现象，即含粗差的真误差估值明显大于其他非正常的观测值的真误差估值。
拟准观测的选取非常关键，在没有特别先验信息的情况下，可采用初选和复选两步方案。
初选时，将观测值按一定规则分为4类：0类很可能含粗差；1类观测结构不好；3类含粗差可能性较小；其余为2类。实际操作中根据观测值计算指标对观测值分类，取第3类和标准化残差较小的2类观测值作为初选拟准观测值，参考文献[14-15]中对初选流程有详细阐述，此处不再赘述。
通过初选的拟准观测，可以计算真误差估值$\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}$，进而计算相应指标进行复选观测。
初次复选时, 计算指标为：
式中，${C_1} = 1.483{\text{median}}\left( {\left| {\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}} \right|} \right)$。
将$\left| {{{\hat \Delta }_i}} \right|$或|WiⅠ|各自按大小排序，排序后对相邻数据求差，若$\left| {{{\hat \Delta }_i}} \right|$和|WiⅠ|分群明显，则差值序列将有明显跳值, 且跳值点之下的$\left| {{{\hat \Delta }_i}} \right|$和|WiⅠ|值都比较小，此处分界，分界线以下复选为拟准观测。
初次复选后的分界线往往不够明显，需要迭代计算更为准确的复选指标，并且采用量化的分群准则，一般计算如下指标：
式中，${C_2} = {\left[{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}_r^{\text{T}}{{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_r}/\left( {r-1} \right)} \right]^{1/2}}$，其中${{{\mathit{\boldsymbol{ \boldsymbol{\hat \varDelta} }}}_r}}$为上一次复选得到的拟准观测的真误差估值。可采用拟准观测的3倍中误差为阈值进行分群，判定粗差，即将指标WiⅡ大于3倍中误差的观测判定为含有粗差。
也可对真误差估值按大小排序，计算以下指标：
式中，P为观测权。可对按大小排序的真误差估值差分后，对差分序列使用2倍的中误差分群，即从指标WiⅢ大于2倍中误差的位置进行分群。
多次复选后，若分群仍不明显，说明没有粗差或选择的拟准观测不正确，若拟准观测的单位权中误差比较大，说明可能有粗差，更换拟准观测后返回拟准复选过程，否则认为没粗差。
-
通过上述的处理过程，残余粗差的量值不大于拟准观测的3倍中误差，保证绝大多数的粗差被精确探测和剔除，此后可以对观测值的钟跳进行判定。判定方法有两种，第1种是对已经经过两步法剔除了粗差的数据再次计算历元单差序列，依据两步法中的第1步历元单差法计算的阈值对单差序列进行疑似钟跳的探测，探测结果即为真实钟跳点；第2种方法是分别利用疑似钟跳点之前和之后的数据拟合钟跳点钟差，并设定一个阈值，若拟合的两个值之差超过阈值，则判定此为钟跳。
1.1. 粗差阈值设定
1.2. 历元单差法探测粗差和钟跳
1.3. 拟准检定法探测粗差
1.4. 钟跳点判定
-
模拟数据能更好地模拟大粗差和各异常情况，更完整地体现卫星钟差异常探测过程，对卫星钟差实测数据的异常探测是算法的实践，以下分别进行算例分析。
-
卫星钟差一般可表示成下列二次多项式[16]：
式中，toc表示参考时刻；t为计算时刻；a0、a1、a2为系数。
给4个参数toc、a0、a1、a2分别依次取0、4×105 ns、0.86 ns/s、-6×10-6 ns/s2，每10 s得到一个钟差值，共计算200个数据。进行以下步骤, 产生模拟观测数据序列L。
1) 给所有历元加-0.5~0.5 ns的随机噪声，可将此随机噪声序列视为底噪，其中误差为0.283 ns；
2) 给历元9、13、16、49、51、99、100、109、123、133、134、148、154分别加入2.5、2.34、-5、0.77、-0.89、1.5、1.6、-0.66、0.84、0.68、0.54、1.6、-0.99 ns的较小误差；
3) 给历元4、5、20、21、68、71、121、128、141、142、145、150、157、181、184、185分别加入(-0.13, 0.23，0.06，-0.18，-0.17，-0.1，-0.39，-0.22，-0.46，0.16，0.4，-0.12，-0.24，0.37，0.24，0.240 1)×4×105 ns的较大量级误差；
4) 给历元60、61、62分别加入10、9、9.5 ns的误差；
5) 给第120及以后历元加-8×104 ns的偏差，模拟钟跳。
使用4种方案处理模拟数据：
方案1：使用历元单差法探测粗差和钟跳，使用绝对中位数法设置阈值，放大因子取10。
方案2：与方案1方法相同，但放大因子取3。
方案3：第1步，使用历元单差法探测粗差和钟跳，使用绝对中位数法设置阈值，放大因子取10；第2步，对各分段数据利用拟准检定法进行粗差探测，使用拟准观测的3倍中误差进行分群。
方案4：第1步，使用历元单差法探测粗差和钟跳，使用绝对中位数法进行判定，放大因子取10；第2步，对各分段数据利用拟准检定法进行粗差探测，对真误差估值按大小排序并作差分，对差分序列使用拟准观测的2倍中误差分群。
各方案处理结果如表 6所示。
粗差点位置 粗差模拟值/ns 方案1 方案2 方案3 方案4 4 -52 000 √ √ √ √ 5 92 000 √ √ √ √ 9 2.5 × √ √ √ 13 2.34 × √ √ √ 16 -5 × √ √ √ 20 24 000 √ √ √ √ 21 -72 000 √ √ √ √ 49 0.77 × × × √ 51 -0.89 × × √ √ 60 10 × × √ √ 61 9 × × √ √ 62 9.5 × × √ √ 68 -68 000 √ √ √ √ 71 -40 000 √ √ √ √ 99 1.5 × × √ √ 100 1.6 × × √ √ 109 -0.66 × × √ √ 121 -156 000 √ √ √ √ 123 0.84 × × × √ 128 -88 000 √ √ √ √ 133 0.68 × × √ √ 134 0.54 × × × × 141 -184 000 √ √ √ √ 142 64 000 √ √ √ √ 145 160 000 √ √ √ √ 148 1.6 × × √ √ 150 -48 000 √ √ √ √ 154 -0.99 × × √ √ 157 -96 000 √ √ √ √ 181 148 000 √ √ √ √ 184 96 000 √ √ √ √ 185 96 040 √ √ √ √ 注：“√”表示粗差被探测；“×”表示粗差未被探测 Table 6. Analysis Result of Simulation Data
方案1探测出第60、63和120历元的疑似钟跳点，以及使单差序列残差绝对值大于6.24 ns的钟差非连续粗差模拟值，本例中32个模拟粗差探测出16个。
方案2探测出第60、63、101和120历元的疑似钟跳点，以及大部分使单差序列残差绝对值大于1.87 ns的钟差粗差模拟值，本例中32个模拟粗差探测出19个。
方案3采用两步法，探测出疑似钟跳点后可探测小粗差，本例中32个模拟粗差探测出29个，探测出大于3倍底噪中误差的所有粗差模拟值和部分小于3倍底噪中误差但大于2倍底噪中误差的粗差模拟值。
方案4采用两步法，探测出疑似钟跳点后可探测小粗差，本例中32个模拟粗差探测出31个，探测出粗差模拟值中大于2倍底噪中误差的所有粗差模拟值。
粗差探测中难以避免地将部分正常观测值误判为粗差，本例误判统计：方案1为1个，方案2为1个，方案3为3个，方案4为4个。
由算例可见，方案1和方案2中，历元单差法能有效探测真实钟跳点，并消除大部分较大粗差，同时选取较小的放大因子能提高粗差探测能力，但疑似钟跳点数量也随之增加；采用两步法策略有良好的粗差探测能力，方案4探测能力最好。
粗差探测时不可避免地将某些正常观测值误判为粗差，本算例中方案4共判定粗差35个，包括31个模拟粗差和4个误判粗差。误判占总观测样本的2%，其减少了正常观测值个数，但在有大量正常观测样本的情况下，可认为少量的误判不影响正常观测特性。
剔除粗差后可以对钟跳点进行判定，本例方案4判定出唯一钟跳点在第122历元处。这和本例模拟的钟跳点第120历元有差异，原因在于第121历元为真实粗差被探测并剔除，第120历元被误判为粗差而被剔除，剩余的不含粗差的观测数据中钟跳将出现在第122历元。
-
本例针对某次星地双向时间同步钟差测试数据进行粗差剔除和钟跳探测，原始数据采样间隔1 s，数据时长约11 h，如图 1所示。
Figure 1. Primary Satellite Clock Bias Data
粗差剔除和钟跳探测使用§2.1中的方案4，并对钟跳点进行判定，处理结果如图 2所示，图 2中钟跳点局部特写如图 3所示。
Figure 2. Result of Method 4 of Practical Data
Figure 3. Features of Clock Jump Points
从图 2和图 3可以看出无明显粗差存在，钟跳被良好探测和标记。卫星钟的正常调钟引起的钟跳往往持续时间较长且跳值较大，而从图 2与图 3钟跳的分布特点和跳值可知这些钟跳并不是由正常的卫星调钟引起，对这些异常钟跳进行标记在实践中是有益的。
至此对粗差的剔除和钟跳的判定已经完成，对探测出的钟跳值的后续处理应根据实际需求进行处理。
2.1. 模拟数据校验分析算例
2.2. 实测数据验证分析
-
本文依据历元单差法和拟准检定法的算法特点，分析了数据异常的典型形式，探讨使用历元单差法探测钟跳和粗差的方法和步骤，总结拟准检定法探测粗差的过程，提出将二者结合应用于卫星钟差数据异常探测的两步法策略，即先利用历元单差法探测疑似钟跳并剔除部分粗差，再对按疑似钟跳分段的数据使用拟准检定法。分别通过模拟数据和实测数据进行算例分析，得出以下结论：
1) 历元单差算法简单，通过对历元单差序列使用残差的绝对中位数法检定，可有效探测粗差和疑似钟跳点，其探测的疑似钟跳点包含真实钟跳；
2) 使用两步法策略，按历元单差法确定的疑似钟跳点对数据分段，每个段内使用拟准检定法进行粗差探测，结果表明，该方法可精确探测小粗差，探测能力接近底噪的2倍中误差。
两步法策略步骤清晰，通过模拟数据和卫星钟差实测数据的分析，表明两步法策略能够精确剔除粗差和探测钟跳。
DownLoad:
