-
摘要: 首先给出典型的原子钟时差观测量模型，包括确定性部分（时差、频差、线性频漂和周期性波动项）、随机性部分（即原子钟噪声）和观测噪声；分析了各分量对应的Allan偏差的表达式。针对部分文献对Kalman滤波器估计原子钟状态原理描述不清晰的问题，描述了原子钟随机微分方程模型和各物理量的含义，从最优估计和低通滤波器两个角度阐述其原理。针对观测噪声过大、存在周期性波动等原因造成无法准确估计原子钟噪声强度的情况，提出了综合Kalman滤波器状态估计结果和Allan偏差图，估计原子钟噪声和观测噪声强度的方法；提出了3种不同的估计线性频漂幅度的方法，并通过实测数据相互验证；针对周期性波动在时差中不明显的问题，结合原子钟随机微分方程模型，提出了综合Kalman滤波器状态估计的结果和对数Allan偏差图估计周期性波动周期和幅度的方法。对两台国产氢钟的实测数据进行了验证，证明该方法物理原理清晰，操作简便易行，具有实用性。通过该方法可以外推得到所有平滑时间的Allan偏差估计值。Abstract: This paper describes atomic clock models and frequency stability analyses methods. The time deviation observation model is comprised of the deterministic part (the time deviation, the frequency deviation, the linear frequency drift and the periodical part), the stochastic part which is the clock noise and the observation noise. This paper gives the Allan deviation expressions of these parts. The principle of utilizing Kalman filter to estimate the clock status is not illustrated in some paper clearly. Thus, this paper describes the atomic clock stochastic differential equations and parameters in detail. Then, we propose a method to estimate the clock noises and observation noise levels by means of the Kalman filter estimation results and Allan deviation pictures. Three methods of estimating the linear frequency drift level are proposed and validated by measurements. This paper also proposes the method of estimating the period and level of the periodical part by means of the atomic clock stochastic differential equations, the Kalman filter estimation results and Allan deviation pictures. The real measurements of two domestic hydrogen masers are used to validate these methods. The principles of these methods are distinct. They are practical and easily to be realized. The methods can be used to obtain the estimated Allan deviation values at any observation interval. This research is the basis of a series of successive researches such as time scale algorithms, clock prediction algorithms, steering algorithms.
-
Keywords:
- atomic clock /
- time deviation model /
- frequency stability /
- Allan deviation /
- Kalman filter
-
守时实验室需要建立和维持一个准确、稳定、可靠的时间尺度作为时间基准[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项对应确定度分量，其中x0为初始时差，y0为初始频差，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为幅度；f0为周期性波动的频率；φ为初始相位。
这时，频差周期性波动在时差上表现为：
$$ \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)为时差的随机性部分，即原子钟时差的噪声；W1(t)和W2(t)分别代表两个独立的维纳过程，并且有W(t)~N(0, t)，即每个维纳过程服从均值为0、方差为时间t的正态分布；σ1和σ2分别是这两个维纳过程的扩散系数，用于表明噪声的强度。
把式（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) 式中，σy2(τ)代表平滑时间为τ时的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和x2分别代表了两个状态分量，x和式（1）完全相同，代表时差，x2代表频差的一部分，即不含WFM的频差；T为时间间隔，T=tk+1-tk；其他符号的含义和§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]将x2理解为频差，但笔者认为，频差是由时差的差分得到的，其中依然包含了WFM和RWFM，而x2中只包含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，理论上可以通过斜率拟合出σ2、σ12和σ22的估计值，记为${\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) 设置一个σ12的大致估计值，代入式（23）确定Q的值，使用KF估计得到Q的估计值，并画出$\hat x\left(t \right)$的Allan偏差图；理论上讲，$\hat x\left(t \right)$的噪声分量应该只包含WFM和RWFM；
3) 根据$\hat x\left(t \right)$的Allan偏差图，反复调整σ12的估计值，重复步骤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强度，即σ12远大于σ22，所以从直观上理解，方法3的估计不确定度最小。
2.4 分析周期性波动分量
1）分析f0。按照§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 s2，${\hat \sigma _1}^2$=3×10-26 s和${\hat \sigma _2}^2$=1.2×10-33 s-1。其中，${\hat \sigma ^2}$=1×10-22 s2意味着观测噪声的标准差为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$。
![]()
图 1 z(t)和$\hat x\left(t \right)$的Allan偏差Figure 1. Allan Deviations of z(t) and $\hat x\left(t \right)$3.2 采用多种方法分析和验证线性频漂的幅度
本节分别对§2.3的3种方法进行验证。
方法1：对$\hat x\left(t \right)$分别用一阶多项式和二阶多项式进行最小二乘拟合，得到线性频漂的估计值为$\hat d$=-3.891×10-20 s/s2。图 2(a)和2(b)分别画出了拟合残差。对比图 2(a)和2(b)，可以明显看出Hm1中存在线性频漂。
![]()
图 2 $\hat x\left(t \right)$的一阶残差和二阶残差Figure 2. First Order and Second Order Residuals of $\hat x\left(t \right)$方法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/s2，和方法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的方法可以有效估计周期性波动的周期和幅度。
![]()
图 4 周期性波动项和$\hat x\left(t \right)$的Allan偏差Figure 4. Allan Deviations of the Sinusoidal Component and $\hat x\left(t \right)$3.4 实例2分析
采用和§2相同的方法分析第2台国产氢钟Hm2的性能，得到该氢钟的参数估计值分别为：${\hat \sigma ^2}$=1×10-22 s2，${\hat \sigma _1}^2$=3×10-26 s，${\hat \sigma _2}^2$=8×10-33 s-1，$\hat d$=-3.8×10-20 s/s2，${\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偏差。
![]()
图 5 z(t)和$\hat x\left(t \right)$的Allan偏差Figure 5. Allan Deviations of z(t) and $\hat x\left(t \right)$把§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偏差估计值。本文提出的方法物理原理清晰，操作简便易行。该研究对于时间尺度、钟差预测、原子钟驾驭等算法具有重要意义。
-
![]()
图 1 z(t)和$\hat x\left(t \right)$的Allan偏差
Figure 1. Allan Deviations of z(t) and $\hat x\left(t \right)$
![]()
图 2 $\hat x\left(t \right)$的一阶残差和二阶残差
Figure 2. First Order and Second Order Residuals of $\hat x\left(t \right)$
![]()
图 4 周期性波动项和$\hat x\left(t \right)$的Allan偏差
Figure 4. Allan Deviations of the Sinusoidal Component and $\hat x\left(t \right)$
-
[1] Whibberley P B, Davis J A, Shemar S L. Local Representations of UTC in National Laboratories[J]. Metrologia, 2011, 48(7):154-164 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=70ecf79baf14ffaa21430975a4719ee6
[2] Patrizia T. Statistical and Mathematical Tools for Atomic Clocks[J]. Metrologia, 2008, 45(6):183-192 doi: 10.1088/0026-1394/45/6/S24
[3] 伍贻威, 朱祥维, 龚航, 等.建立GNSS时间基准的构想和思考[J].电子学报, 2017, 45(8):1 818-1 826 http://d.old.wanfangdata.com.cn/Periodical/dianzixb201708003 Wu Yiwei, Zhu Xiangwei, Gong Hang, et al. Concepts and Thoughts of Forming a GNSS Time Scale[J]. Acta Electronica Sinica, 2017, 45(8):1 818-1 826 http://d.old.wanfangdata.com.cn/Periodical/dianzixb201708003
[4] 伍贻威, 龚航, 朱祥维, 等.单状态变量和两级Kalman滤波器时间尺度算法[J].中国科学:物理学, 力学, 天文学, 2016, 46(6):069502 http://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201606012.htm Wu Yiwei, Gong Hang, Zhu Xiangwei, et al. One-State and Twice Kalman Filter Time Scale Algorithms[J]. Scientia Sinica Physica, Mechanica & Astronomica, 2016, 46(6) :069502 http://www.cnki.com.cn/Article/CJFDTOTAL-JGXK201606012.htm
[5] Wu Y, Zhu X, Huang Y, et al. Uncertainty Derivation and Performance Analyses of Clock Prediction Based on Mathematical Model Method[J]. IEEE Transactions on Instrumentation and Measurement, 2015, 64(10):2 792-2 801 doi: 10.1109/TIM.2015.2418683
[6] Wu Y, Zhu X, Huang Y, et al. Optimal Observation Intervals for Clock Prediction Based on the Mathematical Model Method[J]. IEEE Transactions on Instrumentation and Measurement, 2015, 65(1):1-12 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=2c98dd1fb2ed4422a8bc688bdb1052f7
[7] 伍贻威, 龚航, 朱祥维, 等.原子钟两级驾驭算法及在建立GNSS时间基准中的应用[J].电子学报, 2016, 44(7):1 742-1 750 http://d.old.wanfangdata.com.cn/Periodical/dianzixb201607032 Wu Yiwei, Gong Hang, Zhu Xiangwei, et al. Twice Atomic Clock Steering Algorithm and Its Application in Forming a GNSS Time Reference[J]. Acta Electronica Sinica, 2016, 44(7):1 742-1 750 http://d.old.wanfangdata.com.cn/Periodical/dianzixb201607032
[8] Wu Y, Gong H, Zhu X, et al. A Clock Steering Method: Using a Third-order Type 3 DPLL Equivalent to a Kalman Filter with a Delay[J]. Metrologia, 2015, 52(6):864-877 doi: 10.1088/0026-1394/52/6/864
[9] Wu Y, Gong H, Zhu X, et al. A DPLL Method Applied to Clock Steering[J]. IEEE Transactions on Instrumentation & Measurement, 2016, 65(6):1 331-1 342 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=a5b12403d105bedb189794cdefdb550f
[10] 董绍武.国产守时型氢钟性能的初步测试[J].仪器仪表学报, 2008, 29(8):80-82 http://d.old.wanfangdata.com.cn/Conference/7290200 Dong Shaowu. Preliminary Test of Home-made Hydrogen Clock[J]. Chinese Journal of Scientific Instrument, 2008, 29(8):80-82 http://d.old.wanfangdata.com.cn/Conference/7290200
[11] 张敏, 董绍武, 张正.国产氢原子钟稳定度的测量与统计分析[C].第十七届全国测控计量仪器仪表学术年会, 厦门, 2007 Zhang Min, Dong Shaowu, Zhang Zheng. Measurement and Statistical Analysis of the Stability of Homemade Atomic Hydrogen Master[C]. MCMI'2007, Xiamen, 2007
[12] 张清华, 隋立芬, 贾小林.应用Jones-Tryon Kalman滤波器对在轨GPS Rb钟进行状态监测[J].武汉大学学报·信息科学版, 2012, 37(4):436-440 http://ch.whu.edu.cn/CN/abstract/abstract170.shtml Zhang Qinghua, Sui Lifen, Jia Xiaolin. Monitor State of GPS Rb Clock Using Jones-Tryon Kalman Filter[J]. Geomatics and Information Science of Wuhan University, 2012, 37(4):436-440 http://ch.whu.edu.cn/CN/abstract/abstract170.shtml
[13] 郭海荣, 杨元喜, 何海波, 等.导航卫星原子钟Kalman滤波中噪声方差-协方差的确定[J].测绘学报, 2010, 39(2): 147-149 http://d.old.wanfangdata.com.cn/Periodical/chxb201002007 Guo Hairong, Yang Yuanxi, He Haibo, et al. Determination of Covariance Matrix of Kalman Filter Used for Time Prediction of Atomic Clocks of Navigation Satellites[J].Acta Geodaetica et Cartographica Sinica, 2010, 39(2): 147-149 http://d.old.wanfangdata.com.cn/Periodical/chxb201002007
[14] Jones R H, Tryon P V. Estimating Time from Atomic Clocks[J]. Journal of Research of the National Bureau of the Standards, 1983, 88(1):112-115 http://d.old.wanfangdata.com.cn/NSTLHY/NSTL_HYCC025565343/
[15] 吴海涛, 李孝辉, 卢晓春, 等.卫星导航系统时间基础[M].北京:科学出版社, 2011 Wu Haitao, Li Xiaohui, Lu Xiaochun, et al. Time Basis of Satellite Navigation System[M]. Beijing:Science Press, 2011
[16] 黄观文, 张勤, 许国昌.基于频谱分析的IGS精密星历卫钟差度分析研究[J].武汉大学学报·信息科学版, 2008, 33(5): 496 -499 http://ch.whu.edu.cn/CN/abstract/abstract1585.shtml Huang Guanwen, Zhang Qin, Xu Guochang. IGS Precise Satellite Clock Model Fitting and Its Precision by Using Spectral Analysis Method[J]. Geomatics and Information Science of Wuhan University, 2008, 33(5): 496-499 http://ch.whu.edu.cn/CN/abstract/abstract1585.shtml
[17] Huang Guanwen, Zhang Qin, Xu Guochang. Real-Time Clock Offset Prediction with an Improved GPS Satellite Clocks Model[J]. GPS Solution, 2014, 18:95-104 doi: 10.1007/s10291-013-0313-0
[18] Panfilo G, Tavella P. Atomic Clock Prediction Based on Stochastic Differential Equations[J]. Metrologia, 2008, 45(6):108-116 doi: 10.1088/0026-1394/45/6/S16
[19] Galleani L, Tavella P. The Dynamic Allan Variance IV: Characterization of Atomic Clock Anomalies[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2015, 62(5):791-801 doi: 10.1109/TUFFC.2014.006733
[20] Zucca C, Tavella P. The Clock Model and Its Relationship with the Allan and Related Variances[J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2005, 52(2):289-296 doi: 10.1109/TUFFC.2005.1406554
[21] Galleani L, Tavella P.Time and Kalman Filter[J]. IEEE Control System Magazine, 2010, 30(2):44-65 doi: 10.1109/MCS.2009.935568
[22] 郭海荣.导航卫星原子钟时频特性分析理论与方法研究[D].郑州: 信息工程大学, 2006 http://cdmd.cnki.com.cn/article/cdmd-90008-2007051629.htm Guo Hairong. Research on Theory and Method of Time-Frequency Characteristic Analysis of Navigation Satellite Atomic Clock[D]. Zhengzhou: Information Engineering University, 2006 http://cdmd.cnki.com.cn/article/cdmd-90008-2007051629.htm
-
期刊类型引用(21)
1. 李国俊，史丰丰，徐金锋，赵润，付桂涛，杨玉婷，李兆南，孙文龙，程梦飞，林勇昕. 基于钙离子光钟的160 d连续守时性能分析. 光学学报. 2025(04): 126-133 . 
百度学术
2. 崔海波，王伟，武建锋，李帅辰，李冲，王康. 机动状态下国产铯原子钟稳定性分析. 时间频率学报. 2024(01): 26-33 . 百度学术
3. 宋鑫，张秀荣，李卫东. 关联噪声对维纳过程和拉比振荡的影响. 山西大学学报(自然科学版). 2024(03): 578-582 . 百度学术
4. 李方能，梁益丰，许江宁，吴苗，朱兵. 基于鲁棒双参数指数平滑法的BDS卫星钟差预报. 中国惯性技术学报. 2024(07): 645-653 . 百度学术
5. 王芳敏，李汶林，戴鸿飞，李春怡，周建华，薛申辉，王波. A real-time performance improvement method for composite time scale. Chinese Physics B. 2024(09): 354-361 . 百度学术
6. 李宗源，袁海波，张虹，张继海，刘素芳，王一桁. 基于低通滤波器的原子钟时差测量数据降噪处理及其性能分析. 时间频率学报. 2024(03): 171-179 . 百度学术
7. 李昂，周铁中，薛潇博，易航，陈德好. 一种智能整定氢原子频标PID控制系统参数的方法. 宇航计测技术. 2024(05): 39-44 . 百度学术
8. 曹建峰，孔静，满海钧，鞠冰，张宇，刘荟萃. 天问一号星载USO的长期漂移标定. 武汉大学学报(信息科学版). 2024(12): 2181-2186 . 百度学术
9. 赵军华，张玉军，刘安斐，展昕，陈亮. 附加钟特性的原子钟预测不确定度分析. 测绘科学. 2024(11): 30-39 . 百度学术
10. 阚昊宇，胡志刚，吕逸飞，谢新，周仁宇，赵齐乐. 利用不同时间同步体制钟差评估北斗三号星载原子钟性能. 武汉大学学报(信息科学版). 2023(04): 604-610 . 百度学术
11. 李方能，梁益丰，许江宁，吴苗. BDS/GPS新型铷原子钟长期特性分析. 中国惯性技术学报. 2023(05): 452-461 . 百度学术
12. 李帅辰，武建锋，崔海波，方婧. 国产氢原子钟移动守时性能测试与分析. 导航定位学报. 2023(03): 38-44 . 百度学术
13. 郭文飞，朱萌萌，辜声峰，左鸿铭，陈金鑫. GNSS精密时频接收机时钟调控模型与参数设计方法. 武汉大学学报(信息科学版). 2023(07): 1126-1133 . 百度学术
14. 安景浩，王怀智，梁益丰，许江宁. BDS快速钟差预报方法对比分析. 舰船电子工程. 2023(03): 55-60 . 百度学术
15. 卢鋆，武建峰，袁海波，申建华，孟轶男，宿晨庚，陈颖. 北斗三号系统时频体系设计与实现. 武汉大学学报(信息科学版). 2023(08): 1340-1348 . 百度学术
16. 秦晓伟，杜二旺，赵春勃，王国永，杨涛. 基于频率源稳定度的钟差序列生成方法研究. 空间电子技术. 2023(04): 66-71 . 百度学术
17. 陈汗龙，董哲，周真帆，郑晓雪. 国产原子钟频率稳定度评估分析. 导航定位学报. 2023(06): 28-33 . 百度学术
18. 梁益丰，许江宁，吴苗，何泓洋. 北斗卫星钟差的CEEMDAN分解与周期项提取方法. 中国惯性技术学报. 2022(04): 476-484 . 百度学术
19. 张杰梁，陈鲤文，陈静，张煌辉，姜立斌. 守时系统钟组规模技术验证. 电子质量. 2022(10): 151-157 . 百度学术
20. 伍贻威，王世超. 不通过自由纸面时建立时间基准的方法与性能分析. 测绘学报. 2021(03): 343-355 . 百度学术
21. 惠恬，赵高长，苏军，龚莹莹. 原子钟数据的改进经验模态分解降噪. 河南科技大学学报(自然科学版). 2021(05): 45-50+56+7 . 百度学术
其他类型引用(11)
下载:
计量
- 文章访问数: 2195
- HTML全文浏览量: 359
- PDF下载量: 276
- 被引次数: 32
