留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

原子钟模型和频率稳定度分析方法

伍贻威 杨斌 肖胜红 王茂磊

伍贻威, 杨斌, 肖胜红, 王茂磊. 原子钟模型和频率稳定度分析方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
引用本文: 伍贻威, 杨斌, 肖胜红, 王茂磊. 原子钟模型和频率稳定度分析方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
WU Yiwei, YANG Bin, XIAO Shenghong, WANG Maolei. Atomic Clock Models and Frequency Stability Analyses[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
Citation: WU Yiwei, YANG Bin, XIAO Shenghong, WANG Maolei. Atomic Clock Models and Frequency Stability Analyses[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058

原子钟模型和频率稳定度分析方法

doi: 10.13203/j.whugis20180058
详细信息
    作者简介:

    伍贻威, 博士, 主要从事时间基准与时间同步关键技术研究。Yiwei_Wu_sh@126.com

  • 中图分类号: P228

Atomic Clock Models and Frequency Stability Analyses

More Information
    Author Bio:

    WU Yiwei, PhD, specializes in key technologies of forming time references and time synchronizations. E-mail: Yiwei_Wu_sh@126.com

图(5)
计量
  • 文章访问数:  1293
  • HTML全文浏览量:  131
  • PDF下载量:  211
  • 被引次数: 0
出版历程
  • 收稿日期:  2018-12-16
  • 刊出日期:  2019-08-05

原子钟模型和频率稳定度分析方法

doi: 10.13203/j.whugis20180058
    作者简介:

    伍贻威, 博士, 主要从事时间基准与时间同步关键技术研究。Yiwei_Wu_sh@126.com

  • 中图分类号: P228

摘要: 首先给出典型的原子钟时差观测量模型,包括确定性部分(时差、频差、线性频漂和周期性波动项)、随机性部分(即原子钟噪声)和观测噪声;分析了各分量对应的Allan偏差的表达式。针对部分文献对Kalman滤波器估计原子钟状态原理描述不清晰的问题,描述了原子钟随机微分方程模型和各物理量的含义,从最优估计和低通滤波器两个角度阐述其原理。针对观测噪声过大、存在周期性波动等原因造成无法准确估计原子钟噪声强度的情况,提出了综合Kalman滤波器状态估计结果和Allan偏差图,估计原子钟噪声和观测噪声强度的方法;提出了3种不同的估计线性频漂幅度的方法,并通过实测数据相互验证;针对周期性波动在时差中不明显的问题,结合原子钟随机微分方程模型,提出了综合Kalman滤波器状态估计的结果和对数Allan偏差图估计周期性波动周期和幅度的方法。对两台国产氢钟的实测数据进行了验证,证明该方法物理原理清晰,操作简便易行,具有实用性。通过该方法可以外推得到所有平滑时间的Allan偏差估计值。

English Abstract

伍贻威, 杨斌, 肖胜红, 王茂磊. 原子钟模型和频率稳定度分析方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
引用本文: 伍贻威, 杨斌, 肖胜红, 王茂磊. 原子钟模型和频率稳定度分析方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
WU Yiwei, YANG Bin, XIAO Shenghong, WANG Maolei. Atomic Clock Models and Frequency Stability Analyses[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
Citation: WU Yiwei, YANG Bin, XIAO Shenghong, WANG Maolei. Atomic Clock Models and Frequency Stability Analyses[J]. Geomatics and Information Science of Wuhan University, 2019, 44(8): 1226-1232. doi: 10.13203/j.whugis20180058
  • 守时实验室需要建立和维持一个准确、稳定、可靠的时间尺度作为时间基准[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偏差图估计周期性波动周期和幅度的方法。

    • 典型的氢钟和铷钟,其确定性分量用二次多项式来表示;典型的铯钟,其确定性分量用一次多项式来表示。二次多项式模型[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;σ用于表明观测噪声的强度。

    • 时域上通常用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方差在中间某一平滑时间段内会凸起来。

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

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

    • 式(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)

      式中,xx2分别代表了两个状态分量,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),QR的值符合${\sigma ^2}$、${\sigma _1}^2$和${\sigma _2}^2$时,KF估计得到的$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$是最小均方意义下的最优估计。从频域角度,确定QR的值相当于确定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)$。

    • 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$。

    • 多种方法可以估计线性频漂$\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的估计不确定度最小。

    • 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$。

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

    • 取某一段长度约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)$

    • 本节分别对§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的结果基本吻合。

      图  3  ${\hat x_2}\left(t \right)$

      Figure 3.  ${\hat x_2}\left(t \right)$

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

    • 采用和§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偏差估计值。

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

参考文献 (22)

目录

    /

    返回文章
    返回