在时间保持工作中，时间尺度相当于一台虚拟原子钟，这台虚拟原子钟是通过测量钟差实现的物理钟的加权平均。虚拟钟一般通过与某一台钟的偏差来表示。通常认为虚拟钟的短期频率稳定度与长期频率稳定度都优于任何单台原子钟^[1]。时间保持实验室通常使用的时间尺度算法主要包括ALGOS算法、AT1算法和Kalman算法。ALGOS算法是一种典型的加权平均算法，2011年和2014年，国际度量衡局分别对ALGOS算法中的频率预报模型和权重进行了改进^[2-3]。使用AT1算法的实验室主要有美国国家标准与技术研究院（National Institute of Standard and Technology，NIST）。AT1（NIST）算法是一种实时的原子时算法，也属于经典的加权平均算法。AT1（NIST）算法没有存储过去频率的真实值，只是考虑到频率的变化，用这种方法忽略了原子钟长期波动的信息^[4-6]。而Kalman算法是一种实时算法，通过最优估计理论计算时间尺度，放弃了权重的概念^[7]，对参考钟和理想钟之间的差作最优估计，并将此值作为修正量，计算时间尺度。该算法所采用的最优估值是通过一个Kalman滤波器来实现的^[8-10]。

作为最优的实时算法，Kalman滤波得到了深入研究^[11-15]。文献[11-12]通过研究模型误差的修复和噪声协方差阵的估计方法改进Kalman滤波算法，文献[8, 13]分析研究了Kalman滤波算法用于时间尺度计算的问题。文献[14]研究了自适应Kalman滤波解的性质。Kalman滤波时间尺度算法中状态模型随机误差噪声主要包括频率白噪声、频率随机游走噪声和频率随机奔跑噪声^[15-19]。模型误差噪声主要为频率闪变噪声，可以利用马氏过程模拟计算噪声参数^[20]。通常情况下，时域中误差噪声参数通过Allan方差或Hadamad方差估计。但是，根据“过时”原子钟的数据进行先验估计具有滞后性, 并且估计误差噪声本身具有一定误差。因此, 原子钟状态模型可能存在较大误差，导致Kalman滤波时间尺度算法的状态估计出现偏差，相应的时间尺度产生误差。鉴于此，本文基于渐消因子限制模型误差引起的状态估计偏差，即采用渐消因子限制Kalman滤波器的记忆长度，以便充分利用现实的观测数据增加时间尺度的准确性和稳定性。统计学界早在19世纪60年代和70年代就提出渐消滤波算法^[21-22]，渐消因子在模型预测与异常监测应用中取得了一定效果^[23-24]。

1 渐消Kalman滤波时间尺度算法 1.1 原子钟的Jones-Tryon模型

设$ t_k$表示测量时刻，当原子钟$ i\left( {i = 1 \cdots n} \right)$到达时刻$ t_k$时，第$ i$台钟的时间$ {h_i}\left( {{t_k}} \right)$产生一个脉冲。这个时间不能独立测量，只能通过测量两台钟脉冲的间隔来实现。假定无测量噪声，测量钟差可以表示为两台钟时间之差。钟$ i$在$ t_k$时刻的时差定义为$ {x_i}\left( {{t_k}} \right)$，钟$ i$的随机模型表示为：

$ \frac{{{\rm{d}}{x_i}\left( t \right)}}{{{\rm{d}}t}} = {y_i}\left( t \right) + {n_{{x_i}}}\left( t \right)\frac{{{\rm{d}}{y_i}\left( t \right)}}{{{\rm{d}}t}} = {n_{{y_i}}}\left( t \right) $

(1)

式中，$ {n_{{x_i}}}(t)$和$ {n_{{y_i}}}(t)$是白噪声；$ {{y_i}}(t)$是频率部分。在测量时刻$ {t_{k - 1}}$与$ {t_{k}}$求解式(1)，可以得到离散时间模型：

$ {x_i}\left( t \right) = {x_i}\left( {{t_{k - 1}}} \right) + {y_i}\left( {{t_{k - 1}}} \right){\tau _0} + {w_{{x_i}}}\left( {{t_k}} \right) $

(2)

$ {y_i}\left( {{t_k}} \right) = {y_i}\left( {{t_{k - 1}}} \right) + {w_{{y_i}}}\left( {{t_k}} \right) $

(3)

式中，$ {\tau _0} = {t_{k + 1}} - {t_k}$表示连续两次之间的测量间隔；$ {w_{{x_i}}}({t_k})$与$ {w_{{{\rm{y}}_i}}}({t_k})$分别表示$ {t_k}$时刻钟i的相位噪声和频率噪声。对于$ n$台钟的离散时间模型，可用矩阵表示为：

$ \mathit{\boldsymbol{X}}\left( {{t_k}} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\tau _0}} \right)\mathit{\boldsymbol{X}}\left( {{t_{k - 1}}} \right) + \mathit{\boldsymbol{W}}\left( {{t_k}} \right) $

(4)

式中，$ \mathit{\boldsymbol{X}}({t_k})$是系统状态；$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}({\tau _0})$是一个由对角块组成的$ 2n \times 2n$的矩阵; 过程噪声$ \mathit{\boldsymbol{W}}({t_k})$是一个均值为0向量、协方差阵为$ \mathit{\boldsymbol{Q}}({\tau _0})$的随机向量，$ \mathit{\boldsymbol{Q}}({\tau _0})$是由对角块$ \left[ {\begin{array}{*{20}{c}} {{q_{{x_i}}}{\tau _0} + {q_{{y_i}}}\tau _0^3/3}&{{q_{{y_i}}}\tau _0^2/2}\\ {{q_{{y_i}}}\tau _0^2/2}&{{q_{{y_i}}}{\tau _0}} \end{array}} \right]$组成的$ 2n \times 2n$的矩阵, $ 2{q_{{x_i}}}$和$ 2{q_{{y_i}}}$分别为白噪声$ {n_{{x_i}}}\left( t \right)$和$ {n_{{y_i}}}\left( t \right)$的谱密度。

1.2 渐消Kalman滤波时间尺度算法及原理

$ {t_k}$时刻测量钟差通过时差可以表示为：

$ {x_{ij}}\left( {{t_k}} \right) = {x_i}\left( {{t_k}} \right) - {x_j}\left( {{t_k}} \right) + {e_{ij}}\left( {{t_k}} \right) $

(5)

式中，$ {e_{ij}}({t_k})$表示$ t_k$时刻测量噪声。

原子钟比对数据测量方程的矩阵表示为：

$ \mathit{\boldsymbol{Z}}\left( {{t_k}} \right) = \mathit{\boldsymbol{HX}}\left( {{t_k}} \right) + \mathit{\boldsymbol{V}}\left( {{t_k}} \right) $

(6)

式中，$ \mathit{\boldsymbol{Z}}\left( {{t_k}} \right) = {\left[ {{x_{21}}\left( {{t_k}} \right) \cdots {x_{n1}}\left( {{t_k}} \right)} \right]^{\rm{T}}}$；H表示测量矩阵；$ \mathit{\boldsymbol{V}}\left( {{t_k}} \right)$表示$ t_k$时刻测量噪声向量。假定测量噪声向量均值为0，协方差矩阵为$ \mathit{\boldsymbol{R}}\left( {{t_k}} \right)$。假设第1台钟作为参考钟，对于3台钟的情况，$ \mathit{\boldsymbol{H}} = \left[ {\begin{array}{*{20}{c}} 1&0&{ - 1}&0&0&0\\ 1&0&0&0&{ - 1}&0 \end{array}} \right]$。Kalman滤波是高斯过程最优滤波的一种有效算法。当观测几何信息和动力学模型及统计信息可靠时，Kalman滤波计算性能较好。但当模型存在较大误差时，状态估计常常出现较大偏差^[25]。原子钟时间尺度的Kalman滤波算法中，由于模型误差估计存在偏差，估计$ n$台钟的状态，只有$ n-1$组测量值，随着时间的累积，估计将出现偏差。本文采用渐消因子限制Kalman滤波器的记忆长度，以便充分利用现实的观测数据，提高估值的准确度。

渐消Kalman滤波估计的状态预测模型及相应的预测协方差矩阵分别为：

$ \mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right) = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\tau _0}} \right)\mathit{\boldsymbol{\hat X}}\left( {{t_{k - 1}}} \right) $

(7)

$ {\mathit{\boldsymbol{ \boldsymbol{\bar \varSigma} }}_{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}} = {\lambda _k}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}\left( {{\tau _0}} \right){\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\hat X}}\left( {{t_{k - 1}}} \right)}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {{\tau _0}} \right) + \mathit{\boldsymbol{Q}}\left( {{\tau _0}} \right) $

(8)

式中，$ \mathit{\boldsymbol{\bar X}}({t_k})$为$ t_k$时刻状态参数预测向量；$ \mathit{\boldsymbol{\hat X}}({t_{k{\rm{ - 1}}}})$为$ {t_{k{\rm{ - 1}}}}$时刻状态参数向量的估值向量；$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\bar X}}({t_k})}}$为$ t_k$时刻$ \mathit{\boldsymbol{\bar X}}({t_k})$的协方差矩阵；$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\hat X}}({t_{k{\rm{ - 1}}}})}}$为$ {t_{k{\rm{ - 1}}}}$时刻$ \mathit{\boldsymbol{\hat X}}({t_{k{\rm{ - 1}}}})$的协方差矩阵；$ {\lambda _k}$为渐消因子。

渐消Kalman滤波状态估计为：

$ \mathit{\boldsymbol{\hat X}}\left( {{t_k}} \right) = \mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right) + {\mathit{\boldsymbol{\bar K}}_k}\left( {\mathit{\boldsymbol{Z}}\left( {{t_k}} \right) - \mathit{\boldsymbol{H\bar X}}\left( {{t_k}} \right)} \right) $

(9)

式中，$ {{\mathit{\boldsymbol{\bar K}}}_\mathit{k}}$为增益矩阵，

$ {\mathit{\boldsymbol{\bar K}}_k} = {\mathit{\boldsymbol{ \boldsymbol{\overline \varSigma} }} _{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\left( {\mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{ \boldsymbol{\overline \varSigma} }} }_{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}}{\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}}\left( {{t_k}} \right)} \right) $

(10)

$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}} = \left( {\mathit{\boldsymbol{I}} - {{\mathit{\boldsymbol{\bar K}}}_k}\mathit{\boldsymbol{H}}} \right){\mathit{\boldsymbol{ \boldsymbol{\overline \varSigma} }} _{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}} $

(11)

式中，$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_k}({t_k})}}$为$ t_k$时刻状态估计向量协方差矩阵。根据Kalman滤波算法时间尺度的定义，在同一实验室，原子钟钟差的测量噪声可以忽略，时间尺度表示为：

$ {x_e}\left( {{t_k}} \right) = {x_i}\left( {{t_k}} \right) - {\hat x_i}\left( {{t_k}} \right) $

(12)

式中，$ {x_e}({t_k})$表示$ t_k$时刻计算的Kalman时间尺度；$ {{\hat x}_i}({t_k})$表示$ t_k$时刻钟i的时差状态估计。

1.3 渐消因子的求解

渐消Kalman滤波算法的关键是渐消因子$ {\lambda _k}$的确定。理论上，若$ {t_{k{\rm{ - 1}}}}$时刻状态估计误差大，则$ {\lambda _k}$应该增大，即$ {{\mathit{\boldsymbol{\hat X}}}_{k{\rm{ - 1}}}}$的误差在$ {{{\mathit{\boldsymbol{\hat X}}}_k}}$中的影响应小。有多种确定渐消因子的方法^[26-27]，考虑到状态模型存在误差或状态受到干扰时，滤波算法仍然保持某种最优性，渐消因子的最优性原理描述如下。

设估计误差$ \mathit{\boldsymbol{r}}({t_k})$为：

$ \mathit{\boldsymbol{r}}\left( {{t_k}} \right) = Z\left( {{t_k}} \right) - \mathit{\boldsymbol{H}} \cdot \mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right) $

(13)

最优性要求不同时刻的估计误差是不相关的，表达式为：

$ {C_j}\left( {{t_k}} \right) = E\left[ {\mathit{\boldsymbol{r}}\left( {{t_{k + j}}} \right) \cdot {\mathit{\boldsymbol{r}}^{\rm{T}}}\left( {{t_k}} \right)} \right] \equiv 0,j \ne 0 $

(14)

经推导，式（14）成立的条件为：

$ {\mathit{\boldsymbol{ \boldsymbol{\overline \varSigma} }} _{\bar X\left( {{t_k}} \right)}} \cdot {\mathit{\boldsymbol{H}}^{\rm{T}}} - {\mathit{\boldsymbol{\bar K}}_k} \cdot {\mathit{\boldsymbol{C}}_0}\left( {{t_k}} \right) \equiv 0 $

(15)

式中，$ {\mathit{\boldsymbol{C}}_0}\left( {{t_k}} \right) = \mathit{\boldsymbol{H}} \cdot {{\mathit{\boldsymbol{ \boldsymbol{\bar \varSigma} }}}_{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}} \cdot {\mathit{\boldsymbol{H}}^{\rm{T}}} + \mathit{\boldsymbol{R}}\left( {{t_k}} \right)$。式(15)是构造最优自适应因子的基础，设函数$ f\left( {{t_k}} \right)$为：

$ f\left( {{t_k}} \right) = {\mathit{\boldsymbol{ \boldsymbol{\overline \varSigma} }} _{\mathit{\boldsymbol{\bar X}}\left( {{t_k}} \right)}} \cdot {\mathit{\boldsymbol{H}}^{\rm{T}}} - {\mathit{\boldsymbol{\bar K}}_k} \cdot {\mathit{\boldsymbol{C}}_0}\left( {{t_k}} \right) $

(16)

最优自适应因子求解要求$ f\left( {{t_k}} \right)$最小，同时考虑到滤波算法的实际误差大于理论误差，文中采用如下算法^[27]：

$ {\lambda _k} = \max \left\{ {1,{\rm{tr}}\left( {{\mathit{\boldsymbol{N}}_k}} \right)/{\rm{tr}}\left( {{\mathit{\boldsymbol{M}}_k}} \right)} \right\} $

(17)

式中，$ {\lambda _k}$表示$ {t_k}$时刻状态参数预测向量协方差矩阵的渐消因子；$ {\rm{tr}}\left( {\rm{ }} \right)$表示矩阵求迹符号；$ {\mathit{\boldsymbol{M}}_k}$和$ {\mathit{\boldsymbol{N}}_k}$为：

$ {\mathit{\boldsymbol{M}}_k} = \mathit{\boldsymbol{H \boldsymbol{\varPhi} }}\left( {{\tau _0}} \right){\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\hat X}}\left( {{t_{k - 1}}} \right)}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\rm{T}}}\left( {{\tau _0}} \right){\mathit{\boldsymbol{H}}^{\rm{T}}} $

(18)

$ {\mathit{\boldsymbol{N}}_k} = {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\bar V}}\left( {{t_k}} \right)}} - \mathit{\boldsymbol{HQ}}\left( {{\tau _0}} \right){\mathit{\boldsymbol{H}}^{\rm{T}}} - \mathit{\boldsymbol{R}}\left( {{t_k}} \right) $

(19)

式中，$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{\mathit{\boldsymbol{\bar V}}({t_k})}}$表示$ t_k$时刻预测残差向量的协方差矩阵；预测残差向量表达式为：

$ \mathit{\boldsymbol{\bar V}}\left( {{t_k}} \right) = \mathit{\boldsymbol{H\bar X}}\left( {{t_k}} \right) - \mathit{\boldsymbol{Z}}\left( {{t_k}} \right) $

(20)

为了敏感地反映当前观测历元动力学模型的误差状态，对预测残差向量的协方差矩阵进行估计^[28]：

$ {\mathit{\boldsymbol{ \boldsymbol{\hat \varSigma} }}_{\mathit{\boldsymbol{\bar V}}\left( {{t_k}} \right)}} = \left\{ {\begin{array}{*{20}{c}} {\frac{{{\lambda _{k - 1}}\mathit{\boldsymbol{\bar V}}\left( {{t_k}} \right){{\mathit{\boldsymbol{\bar V}}}^{\rm{T}}}\left( {{t_k}} \right)}}{{1 + {\lambda _{k - 1}}}},}&{k > 1}\\ {\frac{1}{2}{{\mathit{\boldsymbol{\bar V}}}_0}\mathit{\boldsymbol{\bar V}}_0^{\rm{T}},}&{k = 1} \end{array}} \right. $

(21)

式中，$ {\mathit{\boldsymbol{\bar V}}}_0$为$ k=0$时的预测残差向量。

2 实验与分析

为了验证本文方法的有效性，将本文方法与标准Kalman算法(以下简称Kalman算法)和基于预测残差构造自适应因子的Kalman算法(以下简称A-Kalman算法)作了对比。文中数据来源于国家授时中心时频基准实验室原子钟的比对数据，测量精度为0.01 ns，数据的采样间隔为1 h。取样时间段为2016-01-01—2016-04-29, MJD(modified Julian date)为57 388~57 507.96，总长度为120 d。考虑到计算量，随机选取Cs2959、Cs1011及Cs2573 3台铯原子钟进行实验。

采用Kalman算法，A-Kalman算法和基于渐消因子的改进Kalman算法（以下简称F-Kalman算法）进行分析。算例中3台铯原子钟的状态模型噪声误差取值方法如下：根据原子钟的测量数据计算不同取样时间的Allan方差估计值，频率白噪声参数$ q_1$和频率随机游走噪声参数$ q_2$与Allan方差估计值的关系为：

$ {\rm{ADE}}{{\rm{V}}^2}\left( \tau \right) = {q_1}{\tau ^{ - 1}} + \frac{1}{3}{q_2}\tau + \cdots $

(22)

式中, ADEV表示Allan偏差，对于原子钟可以忽略方程右侧的省略项，根据最小二乘原理估计参数$ q_1$和$ q_2$的值。对3台铯原子钟状态模型噪声误差参数估计频率白噪声分别为1.92×10^-15、1.62×10^-15、1.07×10^-15 ns²/ns，频率随机游走噪声为0 ns²/ns。

文中对4个连续的区间（I1~I4）用3种方法分别进行了时间尺度估计。I1、I2、I3及I4区间分别表示为MJD(57 388~57 417.96)、MJD(57 418~57 447.96)、MJD(57 448~57 477.96)和MJD(57 478~57 507.96)时间段。I1区间采用相对于世界协调时时间标准的前一时刻各台原子钟的相位偏差和前30 d数据估计的速率值作为3种Kalman滤波时间尺度估计的初始值。I2区间3种算法的Kalman滤波时间尺度估计采用I1区间相应算法得到的最后时刻的相位值估计和速率值估计作为初始值。I3和I4区间根据同样的方法计算初始值。测量噪声忽略，计算周期为30 d。图 1表示Kalman算法、A-Kalman算法与F-Kalman算法的计算结果比较。表 1与表 2给出了3种Kalman算法计算时间尺度的最大值与最小值。

由图 1、表 1和表 2可知，同一算法不同时间段计算的时间尺度准确度不同。这主要是因为在相应计算区间，原子钟的状态模型噪声参数不是常量，利用固定的噪声参数会给原子钟的状态估计引入一定误差。相比于Kalman算法和A-Kalman算法，F-Kalman算法能够进一步提高原子钟状态估计的准确度。可以看出，在I1、I2、I3及I4区间，时间尺度的最大值与最小值都有不同程度的改善，I4区间时间尺度的最大值降低了3.66 ns和0.62 ns，I3区间时间尺度的最小值改善了2.14 ns和0.02 ns。以上4个区间，最大值平均降低了1.65 ns和0.31 ns，最小值平均改善了1.26 ns和0.45 ns。这主要是因为模型误差较大时，F-Kalman尺度算法通过增大了渐消因子降低模型预测值的比重，相应增大了测量值的使用效率，提高了估计的准确度。

分析3种Kalman算法计算时间尺度的稳定度，由图 2和表 3可知，相比Kalman算法和A-Kalman算法，F-Kalman算法计算得到的时间尺度的稳定度有较明显改善。各区间稳定度提高情况分别为：I1区间平均稳定度提高2.42×10^-15和1.40×10^-15, 最大提高为7.40×10^-15和4.20×10^-15(取样时间为16 h)。I2区间平均稳定度提高2.00×10^-15和3.12×10^-16, 最大提高为4.90×10^-15(取样时间为128 h)和8.0×10^-15(取样时间为64 h)。I3区间平均稳定度提高1.03×10^-15和1.80×10^-16, 最大提高为2.80×10^-15(取样时间为16 h)和1.00×10^-15(取样时间为16 h)。I4区间平均稳定度提高3.50×10^-15和2.66×10^-15, 最大提高为6.25×10^-15(取样时间为64 h)和4.00×10^-15(取样时间为1 h)。分析得出，取样时间大于等于16 h，相比于Kalman算法，F-Kalman算法的尺度稳定度提高较明显，说明模型状态模型误差主要影响原子钟取样时间大于等于16 h时的稳定度。因此，相比Kalman算法，F-Kalman算法能够改善时间尺度的中长期稳定度。相比A-Kalman算法，F-Kalman算法的尺度稳定度也有提高。

根据标准Kalman滤波算法和F-Kalman算法计算时间尺度，设定两种算法的状态初始估计及相应的协方差矩阵相同，对于F-Kalman算法，根据式（8），将渐消因子引入到状态预测协方差矩阵进行计算。根据式（17）计算渐消因子的值，区间I1、I2、I3及I4的时间尺度计算过程中渐消因子如图 3所示。可以得出，4个区间的渐消因子都大于等于1，渐消因子等于1时，F-Kalman算法等同于标准Kalman滤波算法；渐消因子大于1时，改进Kalman算法降低了模型信息的权重，渐消因子在改进Kalman滤波时间尺度算法中通过膨胀状态预测协方差矩阵降低了模型信息的影响，增大测量信息对状态估计的贡献。根据实验结果可知，相比标准Kalman滤波算法，时间尺度的估计误差较小，准确度提高。对于时间尺度的稳定度，相比Kalman算法，F-Kalman算法估计误差降低，时间尺度波动变小，稳定度提高。本文算法主要是通过预测残差向量实时判断模型误差。如果预测残差向量较小，判断原子钟状态能够较好地符合模型；如果预测残差向量较大，说明原子钟运行受到异常扰动的影响，通过计算最优渐消因子降低了模型误差的影响。

对于渐消因子取值情况分析如下。文中3台铯原子钟运行于时频基准实验室稳定的环境中，原子钟分别通过双混频时差检测设备和SR620测量设备参加测量比对，未发现数据异常。一方面，通常对于铯原子钟利用两状态模型，但是由于部分铯原子钟长期运行出现老化现象，这时两状态模型已经不能很好地符合原子钟的实际运行状态，出现模型误差。另一方面，原子钟的随机噪声估计存在误差，原子钟的随机噪声主要通过原子钟历史测量数据进行估计。首先根据历史测量数据估计其Allan方差，然后根据估计的Allan方差计算噪声参数。这样当前时刻的噪声参数与历史时刻的噪声参数可能有所差异，通过Allan方差估计噪声参数的过程也可能存在误差。

3 结语

本文提出了基于渐消因子的改进Kalman滤波时间尺度算法，该方法通过对标准Kalman时间尺度算法引入渐消因子，能够对前一时刻的状态估值协方差矩阵进行膨胀，降低前一时刻状态估值的使用效率，即更多地依赖当前历元的观测信息。该方法有效抵制了由于原子钟状态模型误差引入的状态估计偏差。通过铯原子钟建立的时间尺度进行比较分析，相比于标准Kalman算法和A-Kalman算法，该方法提高了时间尺度的准确度和稳定度，适用于原子钟的实时频率驾驭。下一步将在本文算法的基础上，研究原子钟的频率驾驭算法，充分利用本文算法的性能，提高频率驾驭的准确性和稳定性。