一种基于门限自回归模型的参考时间尺度产生方法

摘要

本发明提供了一种基于门限自回归模型的参考时间尺度产生方法，首先计算原子钟相位差，然后通过TAR建模预测出当前各原子钟的速率，计算原子钟速率与预测速率的差值，更新权重后再次计算原子钟速率，连续两次计算得到的原子钟速率相等时，利用相位差校准主钟。本发明运用门限自回归模型能够准确的预报出原子钟相对于参考时间尺度的速率，又根据速率的预报误差确定各原子钟权重的大小，对于性能不稳定的原子钟，能够很好的预测速率并在计算中加以扣除，充分利用了有限的原子钟资源，提高了小型守时实验室的TA的准确度和稳定度。

说明书

技术领域

本发明涉及一种守时系统中参考原子时的产生方法，属于时间频率技术领域。

背景技术

守时系统是由具有一定规模的钟组建立和保持的，采用时间尺度算法的目的是利 用统计的方法把钟组中的各原子钟综合成一个比其中的每一台原子钟更稳定、更准确、更 可靠的综合原子时间尺度，作为主钟频率驾驭的参考时间尺度。在时间尺度算法中，最关键 的两个因素是权重的确定和频率预报模型。当前多数守时实验室最常用的时间尺度算法是 经典Algos算法，比如国际权度局(BIPM)采用该方法得到UTC，日本国家信息与通信技术研 究室NICT以及国家授时中心所采用的方法都与Algos算法类似。Algos算法对铯原子钟和氢 原子钟的频率预报分别采用了一次多项式和二次多项式模型。

对于大型守时实验室的高性能原子钟而言，采用一次或二次模型预报铯原子钟和 氢原子钟的速率是较为合适的，对于庞大钟组中的个别性能差的原子钟，Algos算法会自动 将其权重降低甚至是取0权，几乎不会对综合时间尺度产生任何影响。但是在对守时需求日 益增多的今天，许多新建的小型实验室原子钟数量较少且性能偏差，采用一次或二次模型 已经不能满足对这些性能偏差的原子钟速率的预报要求，并且，降低权重甚至取0权的做法 不能充分利用有限的原子钟资源，因此，需要采用一种新的更适合于小型守时实验室的时 间尺度算法。

发明内容

为了克服现有技术的不足，本发明提供一种基于门限自回归模型的时间尺度产生 方法，能够更准确的预测出原子钟的频率变化趋势，并在计算中加以扣除，同时根据频率预 报误差的大小确定原子钟的权重，有效的利用了钟组中性能差的原子钟，可作为小型守时 实验室监控主钟的重要参考。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1，对于N台原子钟构成的钟组，采集主钟输出UTC(k)和第i个原子钟输出 Clock(i)，计算原子钟相位差数据UTC(k)-Clock(i)；将钟差数据进行预处理，检测和剔除 原子钟的异常数据；

步骤2，对每一计算区间各原子钟的输出数据进行线性拟合，得到原子钟的速率数 据，对速率数据进行TAR建模，预测出当前各原子钟的速率

步骤3，对于氢钟的趋势所做的预测对于铯钟的趋势所 做的预测其中，x_i(t₀)为t₀时刻的时间尺度与第i个原子钟的钟差；则第j个原子 钟与是本次计算得到的时间尺度TA′的相位差$x_j\left(t\right)=\underset{i=1}{\overset{N}{\Sigma }}{\omega }_i\left(t\right)\{h_i^{\prime }\left(t\right)-X_{\mathrm{ij}}\left(t\right)\}=\mathrm{UTC}\left(k\right)-{\mathrm{TA}}^{\prime },$其中，ω_i(t)为第i个原子钟的权重，X_ij(t)为第i个原子钟和第j个原子钟的相位差；

步骤4，根据UTC(k)-Clock(i)得到TA′-Clock(i)，对TA′-Clock(i)进行线性拟合 得到各原子钟的速率y(i,I_k)；

步骤5，计算y(i,I_k)与的绝对偏差计算各原子钟的 权重${\omega }_i=\frac{1/{\sigma }_i^2}{\underset{i=1}{\overset{N}{\Sigma }}1/{\sigma }_i^2},{\sigma }_i^2=\frac{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right){ϵ}_{i,I_k}^2}{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right)},$M_i代表采用的计算区间的个数；

步骤6，以y(i,I_k)代替重复步骤3，再次计算时间尺度，得到UTC(k)-TA″，TA″ 代表此次计算得到的时间尺度；

步骤7，再次根据UTC(k)-Clock(i)得到TA″-Clock(i)，对TA″-Clock(i)进行线性 拟合得到各原子钟的速率y′(i,I_k)；

步骤8，判断y′(i,I_k)与y(i,I_k)是否相等，如果相等则利用相位差UTC(k)-TA″校准 主钟；如果不相等，则重复计算步骤5至步骤8，直到连续两次计算得到的原子钟速率相等 时，利用相位差UTC(k)-TA″校准主钟。

所述的步骤2中，采用上一个计算区间的权重作为本计算区间第一次计算的权重， 对于第一个计算区间的第一次计算，每台原子钟的权重相等。

本发明的有益效果是：基于门限自回归模型的参考时间尺度算法中运用门限自回 归模型能够准确的预报出原子钟相对于参考时间尺度的速率，又根据速率的预报误差确定 各原子钟权重的大小，对于性能不稳定的原子钟，能够很好的预测速率并在计算中加以扣 除，充分利用了有限的原子钟资源，提高了小型守时实验室的TA的准确度和稳定度。

附图说明

图1是本发明涉及硬件的结构框图；

图2是基于门限自回归模型的参考时间尺度产生方法的实现流程图；

图3是门限自回归模型的建模流程图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施 例。

本发明解决了守时钟组中性能差的原子钟的有效利用问题，尤其适用于原子钟资 源少，且性能差异较大的小型守时实验室。

本发明用的具体技术方案是：

假设钟组中有N台原子钟，根据时间尺度基本方程的求解结果：其中，ω_i(t)为第i个原子钟的权重，h′_i(t)为对原子钟的趋势所做的预测，X_ij(t)为两台原子钟的相位差 数据，x_j(t)代表原子钟j与时间尺度TA(t)的相位差，计算时间尺度的目的即在于求得x_MC(t)，即 主钟与时间尺度的相位差，根据TA(t)与主钟的钟差驾驭主钟，最终得到实时的物理信号。求解x_j(t)的关键是h′_i(t)和ω_i(t)的确定，本发明中对于铯钟的趋势所做的预测 对于氢钟的趋势所做的预测这里，x_i(t_k)为t_k时刻TA与钟i的钟差，I_k表示计算时间间隔[t_k,t]，为I_k时间段原子 钟i相对于TA的频率预测值，z为氢原子钟的频率漂移的估计值。可见，的预测和z的 估计是得到时间修正量h′_i(t)的关键。

采用TAR模型实现原子钟相对于参考时间尺度TA的钟速预测。TAR模型即门限自回 归模型(ThresholdAutoregressivemodel)，该模型是通过采用门限值将非线性时间序列 按照状态空间逐步线性化来实现非线性系统的，在实际的应用中并不是门限域划分的越多 越好，一般采用一个门限值两个门限域。本发明中对钟速y(i,I_k)预测采用两个门限域，首 先根据各原子钟相对于TA的历史钟速数据进行TAR建模，假设某台原子钟的模型为TAR(2； k₁,k₂)，延迟量为d，门限值为r，其中2代表门限域的个数，k₁、k₂分别为两个门限域中AR模型 的阶数，具体的建模过程见下文的具体实施方式，则有

$\hat{y}(i,I_k)=\left(\begin{array}{cc}{\theta }_0^1+\underset{m=1}{\overset{k_1}{\Sigma }}{\theta }_m^{1}y(i,I_{k-m})& \begin{array}{cc}if& y(i,I_{k-d})\le r\end{array}\\ {\theta }_0^2+\underset{m=1}{\overset{k_2}{\Sigma }}{\theta }_m^{2}y(i,I_{k-m})& \begin{array}{cc}if& y(i,I_{k-d})>r\end{array}\end{array}\right)$

上式中，k₁、k₂为所划分的两个门限域中自回归模型的阶数，为两个门限域 中自回归模型的系数。对于氢原子钟还需考虑其频率漂移z，本算法中假设氢钟的频率漂移 在一个月内是保持不变的，采用上月原子钟相对于参考时间尺度TA的频率数据进行最小二 乘线性拟合，求得氢钟的频率漂移z。

新的权重算法取权的依据是强调“钟的可预测性”，认为“性能优良的原子钟具有 好的可预测性”，权重算法计算过程如下：

(1)采用上一个计算区间，最终选取的权重作为本计算区间第一次迭代的权重。其 中，对于第一个计算区间的第一次计算，将每台原子钟近似取等权。

(2)每个计算区间的实际频率值y(i,I_k)与预报值的绝对偏差如下：

${ϵ}_{i,I_k}=|y(i,I_k)-\hat{y}(i,I_k)|$

(3)计算得到每台钟的实际频率值与预报值绝对偏差的平方。

(4)利用1天时间的的数值来计算权重。考虑到新的测量数据具有更好的参考 性，因此采用一个滤波方法使新的测量数据对权重计算发挥更大的作用，而过去的测量数 据则发挥较小的作用：

${\sigma }_i^2=\frac{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right){ϵ}_{i,I_k}^2}{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right)}$

其中下表i是第i台钟，M_i代表采用的计算区间的个数。

(5)原子钟的相对权重的计算公式为：

${\omega }_i=\frac{1/{\sigma }_i^2}{\underset{i=1}{\overset{N}{\Sigma }}1/{\sigma }_i^2}$

当钟i的权重超过了最大权限制ω_max＝4/N时，则令ω_i＝ω_max。

本方法的实施所包括的硬件设备主要有原子钟、相位微调仪、比相仪、转换开关和 工控机，如图1所示，本发明的方法主要以软件的形式在工控机上运行，其实现流程见图2。

步骤1：首先采集比相仪输出的原子钟相位差数据UTC(k)-Clock(i)，UTC(k)代表 主钟输出，Clock(i)代表原子钟i的输出，将钟差数据进行预处理，检测和剔除原子钟的异 常数据。

步骤2：对每一计算区间(本实施例以一天为一个计算区间)的数据进行线性拟合 得到原子钟的速率数据，对得到的速率数据进行TAR建模，预测出当前各原子钟的速率 采用上一个计算区间的权重作为本计算区间第一次计算的权重，对于第一个计算 区间的第一次计算，将每台原子钟近似取等权。

步骤3：对于氢钟，将预测值代入铯钟则代入 根据时间尺度基本方程有$x_j\left(t\right)=\underset{i=1}{\overset{N}{\Sigma }}{\omega }_i\left(t\right)\{h_i^{\prime }\left(t\right)-X_{ij}\left(t\right)\}=UTC\left(k\right)-{\mathrm{TA}}^{\prime },$TA′是本次计算得到的时间尺度。

步骤4：根据UTC(k)-Clock(i)得到TA′-Clock(i)，对TA′-Clock(i)进行线性拟合 得到各原子钟的速率y(i,I_k)。

步骤5：计算y(i,I_k)与的绝对偏差根据如下公式计 算各原子钟的权重：

$\left(\begin{array}{cc}{\sigma }_i^2=\frac{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right){ϵ}_{i,I_k}^2}{\underset{m=1}{\overset{M_i}{\Sigma }}\left(\frac{M_i+1-m}{M_i}\right)}& {\omega }_i=\frac{1/{\sigma }_i^2}{\underset{i=1}{\overset{N}{\Sigma }}1/{\sigma }_i^2}\end{array}\right)$

M_i代表采用的计算区间的个数。

步骤6：以yi(,I)_k代替重复步骤3，再次计算时间尺度得到UTC(k)-TA″，TA″ 代表此次计算得到的时间尺度。

步骤7：再次根据UTC(k)-Clock(i)得到TA″-Clock(i)，对TA″-Clock(i)进行线性 拟合得到各原子钟的速率y′(i,I_k)。

步骤8：判断y′(i,I_k)与y(i,I_k)是否相等，如果相等则利用UTC(k)-TA″校准主钟。 如果不相等，则重复计算步骤5至步骤8，直到连续两次计算得到的原子钟速率相等时，利用 相位差UTC(k)-TA″校准主钟。

图3为步骤2中TAR模型的建模过程，具体方法为：

步骤1：将历史速率数据rate_i(t)进行从小到大的有序排序得到rate′_i(t)。

步骤2：给定参数D和K，D和K都是正整数。K表示分段线性AR模型阶数k_j的最大可能 值，D表示延迟数d的最大可能值，为了建模方便也可以设D＝K。比如设D＝K＝5，若计算结果 发现d或k_j的值达到5，则增大D和K的取值。

步骤3：令延迟参数取初值d＝1，采用模型TAR(2；k₁,k₂)，因此仅有一个tq₁作为一 个门限r₁的候选值，则有R₁∈(-∞,tq₁]和R₂∈(tq₁,+∞]。若取tq₁＝rate′_i(q₁)，那么，q₁表 示rate′_i(t)的总资料数的百分比。首先取q₁＝0.3N，利用门限值tq₁将速率数据rate_i(t)划 分为R₁和R₂两个门限域。

步骤4：取n₀＝max(d,K)，若速率数据rate_i(t)中,有rate_i(t)_0-d+1，rate_i(t)_0-d+4， rate_i(t)_0-d+5，…，的一些样本值小于或等于tq₁,而另一些样本值则大于tq₁，那么就构成下 述的两个线性方程：

用矩阵可以表示为

$\left(\begin{array}{c}{Z}_1=A_1{\theta }_1\\ Z_2=A_2\theta \end{array}\right)$

式中，Z₁和Z₂分别表示序列rate_i(t)属于门限域R₁和R₂的序列向量，θ₁和θ₂分别是 待估计的两个线性AR模型系数向量。A₁和A₂分别是两个线性系统的系数矩阵。求解这两个方 程组，用最小二乘法可估计两个模型的自回归系数θ1和θ2。

步骤5：根据最小AIC估计值来选择两个分段线性AR模型阶k₁和k₂的最佳识别估计 和即有

$AIC\underset{0<k_1\le K}{\left({\hat{k}}_1\right)=min}[N_{1}ln\left({\mathrm{RSS}}_1\right(k_1)/N_1)+2(k_1+1)]$

$AIC\left(\underset{0<k_2\le K}{{\hat{k}}_2=min}\right[N_{2}ln\left({\mathrm{RSS}}_2\right(k_2)/N_2)+2(k_2+1)]$

上式中N₁和N₂分别是向量Z₁和Z₂的维数，RSS₁(k₁)和RSS₂(k₂)分别是两个线性系统 的剩余量平方和。所以，两个门限域内AR模型阶和的确定，以及相应的自回归系数和 的估计可以同时获得。

步骤6：对于某个固定的tq₁，有

$AIC\left({\mathrm{tq}}_1\right)=AIC\left({\hat{k}}_1\right)+AIC\left({\hat{k}}_2\right)$

再取tq₁＝0.4N，0.5N，……0.7N，按上述方法可以求得一系列的AIC(tq₁)，由

$AIC\left({\hat{r}}_1\right)=\underset{{\mathrm{tq}}_1}{\min }\left\{AIC\left(tq\right)\right\}$

可确定门限值r₁的估计值这样，就选出了取延迟参数为d＝1时的TAR(2；k₁,k₂) 模型。

步骤7：取延迟参数d＝2,3,……D，根据公式

$AIC\left(\hat{d}\right)=\underset{d}{\min }[AIC\left({\hat{r}}_1\right)/(N-\max (d,k\left)\right)]$

这样就建立了门限自回归模型TAR(2；k₁,k₂)的最佳模型。