守时实验室需要建立和维持一个准确、稳定、可靠的时间尺度作为时间基准^[1-3]，其核心算法包括时间尺度算法^[4]、钟差预测算法^[5-6]、驾驭算法^[7-9]等。原子钟是守时系统的核心部件。时间基准的性能不仅和算法设计有关，还和参与计算的各台原子钟的性能有关。每种算法需要针对各参与计算的原子钟的模型和频率稳定度来设计。因此分析原子钟的模型和稳定度具有重要意义。

目前的研究中，文献[10-13]只分析了指定平滑时间的稳定度，或只对时差观测量中的某一具体分量展开分析，没有具体分析时差观测量中各分量对Allan偏差的贡献，并通过估计值外推得到所有平滑时间的Allan偏差估计值；文献[12-14]对Kalman滤波器估计原子钟状态原理描述不清晰；当观测噪声过大、存在周期性波动时，无法使用斜率法准确直接估计原子钟噪声强度^[15]；当周期性波动在时差中不明显时，目前估计方法较为复杂^{[12, 16-17]}。

本文针对上述问题展开研究，探索原子钟模型和频率稳定度分析方法。大量实验表明，典型氢钟和铯钟的观测模型可以表示为^{[2, 18-20]}：①确定性部分，用二次多项式（时差、频差和线性频漂）加周期性波动项表示；②随机性部分，即原子钟噪声，主要为频率白噪声（white frequency modulation noise，WFM）和频率随机游走噪声（walk random frequency modulation noise，RWFM）；③观测噪声，噪声类型为相位白噪声（white phase modulation noise，WPM）。本文详细分析了该模型各分量的Allan偏差表达式。

在此基础上，本文从最优估计和低通滤波器^[7-9]两个角度描述Kalman滤波器（Kalman filter，KF）估计原子钟状态的原理；提出了综合KF状态估计的结果和Allan偏差图估计原子钟噪声和观测噪声强度的方法；提出了3种不同的估计线性频漂幅度的方法；结合原子钟随机微分方程模型，提出了综合Kalman滤波器状态估计的结果和对数Allan偏差图估计周期性波动周期和幅度的方法。

1 原子钟模型及频率稳定度 1.1 原子钟时差和时差量模型

典型的氢钟和铷钟，其确定性分量用二次多项式来表示；典型的铯钟，其确定性分量用一次多项式来表示。二次多项式模型^[2]表示为：

$ x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + {\varepsilon _x}\left( t \right) $

(1)

式中，等式右边前3项对应确定度分量，其中x₀为初始时差，y₀为初始频差，d为线性频漂；ε_x(t)对应时差的随机性部分，即噪声。一次多项式模型即式（1）中线性频漂d=0的情况。

某些时候，时差的确定性分量中还包含周期性波动分量。为简化分析，假设瞬时频差的周期性波动为标准正弦波或余弦波形式^[19]，即：

$ {y_s}\left( t \right) = A\cos \left( {2{\rm{ \mathit{ π} }}{f_0}t + \varphi } \right) $

(2)

式中，A为幅度；f₀为周期性波动的频率；φ为初始相位。

这时，频差周期性波动在时差上表现为：

$ \begin{array}{*{20}{c}} {{x_s}\left( t \right) = \int_0^t {{y_s}\left( s \right){\rm{d}}s} = A\int_0^t {\cos \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right){\rm{d}}s} = }\\ {\frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right.} \end{array} $

(3)

于是，二次多项式叠加周期性波动项的时差模型表示为：

$ \begin{array}{*{20}{c}} {x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + }\\ {\frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + {\varepsilon _x}\left( t \right)} \end{array} $

(4)

通过大量观测发现，氢钟、铯钟时差中起主导作用的噪声是RWFM^{[2, 18, 20-21]}和WFM。这两种噪声在时差上可以分别用维纳过程和积分维纳过程来建模^{[2, 18, 20-21]}，即：

$ {\varepsilon _x}\left( t \right) = {\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} $

(5)

式中，ε_x(t)为时差的随机性部分，即原子钟时差的噪声；W₁(t)和W₂(t)分别代表两个独立的维纳过程，并且有W(t)~N(0, t)，即每个维纳过程服从均值为0、方差为时间t的正态分布；σ₁和σ₂分别是这两个维纳过程的扩散系数，用于表明噪声的强度。

把式（5）代入式（4），时差可表示为：

$ \begin{array}{*{20}{c}} {x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + \frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + }\\ {{\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} } \end{array} $

(6)

时差的观测量表示为：

$ \begin{array}{*{20}{c}} {z\left( t \right) = x\left( t \right) + \sigma \varepsilon \left( t \right) = {x_0} + {y_0}t + }\\ {\frac{1}{2}d{t^2} + \frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + }\\ {{\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} + \sigma \varepsilon \left( t \right)} \end{array} $

(7)

式中，ε(t)为WPM；σ用于表明观测噪声的强度。

1.2 各分量Allan方差的表达式

时域上通常用Allan方差来表征频率稳定度。Allan方差的定义如下^{[2, 19]}：

$ \sigma _y^2\left( \tau \right) = \frac{1}{2}E\left[ {{{\left( {\bar y\left( {{t_k} + \tau } \right) - \bar y\left( {{t_k}} \right)} \right)}^2}} \right] $

(8)

式中，τ为平滑时间；y为平均频差，定义为：

$ \bar y\left( {{t_k}} \right) = \frac{1}{\tau }\int_{{t_k} - \tau }^{{t_k}} {y\left( t \right){\rm{d}}t} = \frac{{x\left( {{t_k}} \right) - x\left( {{t_k} - \tau } \right)}}{\tau } $

(9)

式中，y为瞬时频差；x即为式（6）中定义的时差。

实际上常用Allan方差的平方根，即Allan偏差来表征稳定度。文献[22]详细推导了扩散系数与Allan方差的关系，即：

$ \sigma _y^2\left( \tau \right) = \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau $

(10)

式中，σ_y²(τ)代表平滑时间为τ时的Allan方差。式（10）中等式右边第一项斜率为-1，第二项斜率为1，分别对应WFM和RWFM。这说明在对数Allan方差图中，WFM的斜率为-1，RWFM的斜率为1。在实际应用中，很容易通过对Allan方差拟合得到扩散系数的值。

文献[20]推导了WPM与Allan方差的关系：

$ \sigma _y^2\left( \tau \right) = 3{\sigma ^2}/{\tau ^2} $

(11)

式（11）表明，WPM在对数Allan方差图中斜率为-2。

把式（1）和式（9）代入式（8），得到线性频漂d与Allan方差的关系：

$ \sigma _y^2\left( \tau \right) = \frac{1}{2}{d^2}{\tau ^2} $

(12)

式（12）表明，d在对数Allan方差图中斜率为2。式（2）所示的周期性波动项和Allan方差的关系为^[19]：

$ \sigma _y^2\left( \tau \right) = {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}} $

(13)

综上，式（6）所示的时差的Allan方差表示为：

$ \sigma _y^2\left( \tau \right) = \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau + \frac{1}{2}{d^2}{\tau ^2} + {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}} $

(14)

式（7）所示的时差观测量的Allan方差表示为：

$ \begin{array}{*{20}{c}} {\sigma _y^2\left( \tau \right) = 3{\sigma ^2}/{\tau ^2} + \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau + }\\ {\frac{1}{2}{d^2}{\tau ^2} + {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}}} \end{array} $

(15)

式中，等式右边第1项为观测噪声；第2、3项为原子钟噪声；第4项为频漂；第5项为周期性波动在Allan方差中的分量。

由式（15）看出，WPM、WFM、RWFM和线性频漂在对数Allan方差图中的斜率分别为-2、-1、1和2，在对数Allan偏差图中的斜率分别为-1、-1/2、1/2和1。假如时差中存在较大幅度的周期性波动，那么Allan方差也存在明显的周期性波动，即Allan方差在中间某一平滑时间段内会凸起来。

但是，并不是每台钟的模型中都含有上述所有分量。例如大部分铯钟几乎没有线性频漂，大部分型号的地面钟都没有周期性波动。

2 时差模型及频率稳定度分析方法

本节首先从最优估计和低通滤波器两个角度描述了KF估计原子钟状态的原理；然后展示分析WPM、WFM、RWFM强度、线性频漂幅度、周期性波动的幅度和频率的方法。

2.1 使用Kalman滤波器估计原子钟状态

式（1）所示的原子钟模型用随机微分方程（stochastic differential equations，SDEs）闭合解的离散形式表示为^{[2, 18, 20-21]}：

$ \left\{ \begin{array}{l} x\left( {{t_{k + 1}}} \right) = x\left( {{t_k}} \right) + {x_2}\left( {{t_k}} \right)T + \frac{1}{2}d{T^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\sigma _1}\left( {{W_1}\left( {{t_{k + 1}}} \right) - {W_1}\left( {{t_k}} \right)} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\sigma _2}\int_{{t_k}}^{{t_{k + 1}}} {\left( {{W_2}\left( s \right) - {W_2}\left( {{t_k}} \right)} \right){\rm{d}}s} \\ {x_2}\left( {{t_{k + 1}}} \right) = {x_2}\left( {{t_k}} \right) + dT + {\sigma _2}\left( {{W_2}\left( {{t_{k + 1}}} \right) - {W_2}\left( {{t_k}} \right)} \right) \end{array} \right. $

(16)

式中，x和x₂分别代表了两个状态分量，x和式（1）完全相同，代表时差，x₂代表频差的一部分，即不含WFM的频差；T为时间间隔，T=t_k+1-t_k；其他符号的含义和§1.1相同。其中，噪声分量为：

$ {J_k} = \left[ {\begin{array}{*{20}{c}} {{\sigma _1}\left( {{W_1}\left( {{t_{k + 1}}} \right) - {W_1}\left( {{t_k}} \right)} \right) + {\sigma _2}\int_{{t_k}}^{{t_{k + 1}}} {\left( {{W_2}\left( s \right) - {W_2}\left( {{t_k}} \right)} \right){\rm{d}}s} }\\ {{\sigma _2}\left( {{W_2}\left( {{t_{k + 1}}} \right) - {W_2}\left( {{t_k}} \right)} \right)} \end{array}} \right] $

(17)

使用KF对上述两个状态分量进行估计，得到估计值记为$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$，其步骤用以下5个方程表示：

$ {\mathit{\boldsymbol{\hat s}}_{k, k - 1}} = \phi {\mathit{\boldsymbol{\hat s}}_{k - 1, k - 1}} $

(18)

$ {P_{k, k - 1}} = \phi {P_{k - 1, k - 1}}{\phi ^{\rm{T}}} + \mathit{\boldsymbol{Q}} $

(19)

$ {\mathit{\boldsymbol{K}}_k} = {P_{k, k - 1}}{H^{\rm{T}}}{\left( {H{P_{k, k - 1}}{H^{\rm{T}}} + R} \right)^{ - 1}} $

(20)

$ {\hat s_{k, k}} = {\hat s_{k, k - 1}} + {K_k}\left( {{z_k} - H{{\hat s}_{k, k - 1}}} \right) $

(21)

$ {\mathit{\boldsymbol{P}}_{k, k}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_k}\mathit{\boldsymbol{H}}} \right){\mathit{\boldsymbol{P}}_{k, k - 1}} $

(22)

式中各符号意义见文献[12, 21]。其中，H=[1 0]；$\phi $=$\left[ {\begin{array}{*{20}{c}} 1&T\\ 0&1 \end{array}} \right]$；过程噪声方差矩阵Q和观测噪声方差矩阵R分别表示为：

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}} = {\rm{E}}\left[ {\left( {{J_k} - 0} \right){{\left( {{J_k} - 0} \right)}^{\rm{T}}}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {\sigma _1^2T + \frac{1}{3}\sigma _2^2{T^3}}&{\frac{1}{2}\sigma _2^2{T^2}}\\ {\frac{1}{3}\sigma _2^2{T^3}}&{\sigma _2^2T} \end{array}} \right]} \end{array} $

(23)

$ \mathit{\boldsymbol{R}} = {\sigma ^2} $

(24)

当原子钟时差符合模型（16），Q和R的值符合${\sigma ^2}$、${\sigma _1}^2$和${\sigma _2}^2$时，KF估计得到的$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$是最小均方意义下的最优估计。从频域角度，确定Q和R的值相当于确定KF的带宽^[8-9]；KF作为低通滤波器，其作用就是滤除时差观测量中高频的WPM，$\hat x\left({{t_k}} \right)$中保留了低频的WFM和RWFM，${\hat x_2}\left({{t_k}} \right)$中只保留了最低频的RWFM^[8-9]。

文献[12-14]将x₂理解为频差，但笔者认为，频差是由时差的差分得到的，其中依然包含了WFM和RWFM，而x₂中只包含RWFM。

由于$\hat x\left({{t_k}} \right)$中包含WFM和RWFM，而${\hat x_2}\left({{t_k}} \right)$中只包含RWFM，噪声分量更少，所以，从${\hat x_2}\left({{t_k}} \right)$中可以更有效地估计频漂、周期性分量等。当含有周期性波动项时，尽管原子钟模型不完全符合式（16），但从KF作为低通滤波器的角度^[8-9]出发，KF依然可以有效滤除观测噪声，估计出$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$。

2.2 使用斜率法和Kalman滤波器分析WPM、WFM和RWFM的强度

WPM、WFM和RWFM在对数Allan偏差图中的斜率分别为-1、-1/2、1/2，理论上可以通过斜率拟合出σ²、σ₁²和σ₂²的估计值，记为${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$。但是，实际上很多情况，尤其是远程比对的情况下，观测噪声很大，造成对数Allan偏差图中WFM分量淹没在WPM分量中。此外，周期性波动幅度很大时，Allan偏差图中间部分会向上凸起。这些都可能造成无法观察到WFM分量。针对这种情况，这里提出一种在${\hat \sigma _1}^2$无法准确获知的情况下，使用KF结合时差观测量z(t)和时差估计值$\hat x\left(t \right)$的Allan偏差，获取${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$的方法。步骤如下：

1) 计算式（7）所示z(t)的Allan偏差，斜率拟合${\hat \sigma _2}^2$和${\hat \sigma ^2}$；

2) 设置一个σ₁²的大致估计值，代入式（23）确定Q的值，使用KF估计得到Q的估计值，并画出$\hat x\left(t \right)$的Allan偏差图；理论上讲，$\hat x\left(t \right)$的噪声分量应该只包含WFM和RWFM；

3) 根据$\hat x\left(t \right)$的Allan偏差图，反复调整σ₁²的估计值，重复步骤2），观察$\hat x\left(t \right)$的Allan偏差图，确保噪声分量斜率为-1/2和1/2，即只含有WFM和RWFM，最终获取${\hat \sigma _1}^2$。

2.3 分析线性频漂的幅度

多种方法可以估计线性频漂$\hat d$。如直接对时差观测量z(t)进行最小二乘拟合（方法1），或对Allan偏差的斜率进行拟合（方法2），或对KF估计得到的${\hat x_2}\left(t \right)$的斜率进行最小二乘拟合（方法3）。

参照文献[5-6]，可以从理论上分析这些方法的估计不确定度，本文不展开分析。按照§2.1的分析，由于KF得到的${\hat x_2}\left(t \right)$中滤除了WFM，只含有RWFM，此外典型氢钟或铯钟WFM强度远大于RWFM强度，即σ₁²远大于σ₂²，所以从直观上理解，方法3的估计不确定度最小。

2.4 分析周期性波动分量

1）分析f₀。按照§2.1方法，尽管式（16）没有对周期性波动建模，KF依然可以估计得到$\hat x\left(t \right)$和${\hat x_2}\left(t \right)$，由于${\hat x_2}\left(t \right)$中只含有RWFM，周期性波动在z(t)和$\hat x\left(t \right)$中不明显，在${\hat x_2}\left(t \right)$中很明显。可以通过观测较长时间段内${\hat x_2}\left(t \right)$波峰波谷的位置，得到波动频率的估计值，记为${\hat f_0}$。

2）分析A。由于噪声的存在，观察${\hat x_2}\left(t \right)$无法准确获取A的估计值，记为$\hat A$。本文采用如下方法：在$\hat x\left(t \right)$的Allan偏差图中，反复调整$\hat A$的值，观察式（13）平方根所示的周期性波动的Allan偏差与$\hat x\left(t \right)$的Allan偏差的吻合程度，最终得到$\hat A$。

3 实验分析

本文采用两台国产氢钟（分别记为Hm1和Hm2）相对于参考时间基准（记为Ref）的实测钟差（Hm1-Ref，Hm2-Ref）进行分析，以验证提出的原子钟模型和频率稳定度分析方法。其中Ref和国际协调世界时同步，稳定度远高于单台国产氢钟。

3.1 使用斜率法和KF分析WPM、WFM和RWFM的强度和线性频漂的幅度

取某一段长度约75 d的Hm1-Ref实测数据，作为时差观测量z(t)。按照§1.2的方法估计${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$，得到${\hat \sigma ^2}$=1×10^-22 s²，${\hat \sigma _1}^2$=3×10^-26 s和${\hat \sigma _2}^2$=1.2×10^-33 s^-1。其中，${\hat \sigma ^2}$=1×10^-22 s²意味着观测噪声的标准差为0.01 ns或10 ps，符合测量设备的指标。

图 1展示了KF前后z(t)和$\hat x\left(t \right)$的Allan偏差，同时也画出了根据上述估计值，按照式（10）、（11）的平方根计算得到的WPM、WFM、RWFM部分的Allan偏差。综上，本实验表明§2的方法可以在观测噪声和周期性波动幅度很大时，有效估计出${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$。

3.2 采用多种方法分析和验证线性频漂的幅度

本节分别对§2.3的3种方法进行验证。

方法1：对$\hat x\left(t \right)$分别用一阶多项式和二阶多项式进行最小二乘拟合，得到线性频漂的估计值为$\hat d$=-3.891×10^-20 s/s²。图 2(a)和2(b)分别画出了拟合残差。对比图 2(a)和2(b)，可以明显看出Hm1中存在线性频漂。

方法2：验证方法1。将方法1得到的$\hat d$值代入式（12）再平方根，计算得到Allan偏差的频漂部分。图 1中画出了Allan偏差的频漂部分，可知两种方法得到的$\hat d$值基本吻合。

方法3：图 3画出了KF估计得到的${\hat x_2}\left(t \right)$，可知${\hat x_2}\left(t \right)$中存在明显的线性频漂。对${\hat x_2}\left(t \right)$的斜率进行最小二乘拟合，得到$\hat d$=-3.86×10^-20 s/s²，和方法1的结果基本吻合。

综上，§2.3的3种方法估计结果基本吻合。实际上，直观上分析，方法3估计不确定度较小，原因在于KF估计得到的${\hat x_2}\left(t \right)$中滤除了WFM，只含有RWFM。假如直接对时差一次差分或者二次差分，得到的序列中同样含有WFM和RWFM差分得到的噪声，由于WFM强度远大于RWFM，将导致无法直接在这些差分序列中观测到图 3中所示的周期性波动，周期性波动已经淹没于噪声中。

3.3 分析周期性波动项

采用§2.4的方法，观察图 3波峰波谷的位置，得到${\hat f_0}$=1/86 400 Hz，$\hat A$大约在1×10^-14~2×10^-14 s/s左右。反复比较式（13）开平方后所示的周期性波动和$\hat x\left(t \right)$的Allan偏差，最终得到$\hat A$=1.6×10^-14 s/s。图 4画出了周期性波动项和$\hat x\left(t \right)$的Allan偏差，可见两者比较吻合。本实验表明了§2的方法可以有效估计周期性波动的周期和幅度。

3.4 实例2分析

采用和§2相同的方法分析第2台国产氢钟Hm2的性能，得到该氢钟的参数估计值分别为：${\hat \sigma ^2}$=1×10^-22 s²，${\hat \sigma _1}^2$=3×10^-26 s，${\hat \sigma _2}^2$=8×10^-33 s^-1，$\hat d$=-3.8×10^-20 s/s²，${\hat f_0}$=1/86 400 Hz和$\hat A$=1.0×10^-14 s/s。图 5画出了KF前后z(t)和$\hat x\left(t \right)$的Allan偏差，以及通过参数估计值计算得到的WPM、WFM、RWFM、线性频漂、周期性波动项的Allan偏差。

把§3.1~§3.4得到的参数估计值代入式（14）和式（15），得出τ < 1 d的实验结果和国产氢钟的说明书相符；当τ > 10 000 s时，观测噪声对Allan偏差的影响很小。

综上，本节实验结果验证了本文方法可以有效分析WPM、WFM、RWFM、线性频漂、周期性波动项各自的Allan偏差，以及总的Allan偏差，并通过这些估计值，拟合出任意平滑时间的Allan偏差估计值。

4 结语

本文展示了原子钟模型和频率稳定度分析方法，详细分析了原子钟时差观测量中的各分量，包括确定性部分（时差、频差、线性频漂和周期性波动项）、随机性部分（WFM、RWFM）和观测噪声（WPM）；分析了WPM、WFM、RWFM、线性频漂、周期性波动项在Allan偏差中的表达式，描述了KF用于估计原子钟状态的原理；提出了当在对数Allan偏差图中，WFM完全淹没于WPM时，使用KF估计WPM、WFM、RWFM强度的方法；提出了3种估计线性频漂幅度的方法和估计周期性波动周期和幅度的方法。通过两台国产氢钟的实测数据验证了本文方法的实用性。实际上，可以通过这些估计值拟合出任意平滑时间的Allan偏差估计值。本文提出的方法物理原理清晰，操作简便易行。该研究对于时间尺度、钟差预测、原子钟驾驭等算法具有重要意义。