搜索

x

留言板

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

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

原子钟噪声变化时改进的Kalman滤波时间尺度算法

宋会杰 董绍武 王翔 章宇 王燕平

原子钟噪声变化时改进的Kalman滤波时间尺度算法

宋会杰, 董绍武, 王翔, 章宇, 王燕平
PDF
HTML
导出引用
  • Kalman滤波时间尺度算法是一种实时的原子钟状态估计方法, 在时间保持工作中具有重要的实用价值. 本文引入自适应因子改进Kalman滤波时间尺度算法, 对Kalman滤波时间尺度算法中状态预测协方差矩阵引入自适应因子, 利用统计量实时计算自适应因子的量值, 控制状态预测协方差矩阵的增长, 降低原子钟状态估计的扰动, 从而提高时间尺度的准确度和稳定度. 模拟数据和实测数据表明, 原子钟噪声强度变化时或噪声强度估计存在误差时改进的Kalman滤波时间尺度算法有效地提高了时间尺度的准确度和长期稳定度.
      通信作者: 宋会杰, songhuijie@ntsc.ac.cn
    [1]

    陈卫东, 刘要龙, 朱奇光, 陈颖 2013 物理学报 62 170506

    Chen W D, Liu Y L, Zhu Q G, Chen Y 2013 Acta Phys. Sin. 62 170506

    [2]

    刘洋洋, 廉保旺, 赵宏伟, 刘亚擎 2014 物理学报 63 228402

    Liu Y Y, Lian B W, Zhao H W, Zhao H W, Liu Y Q 2014 Acta Phys. Sin. 63 228402

    [3]

    林旭, 罗志才 2015 物理学报 64 080201

    Lin X, Luo Z C 2015 Acta Phys. Sin. 64 080201

    [4]

    赵龙 2012 物理学报 61 104301

    Zhao L 2012 Acta Phys. Sin. 61 104301

    [5]

    郭海荣, 杨元喜, 何海波, 徐天河 2010 测绘学报 39 146

    Guo H R, Yang Y X, He H B, Xu T H 2010 Acta Geod. Cartogr. Sin. 39 146

    [6]

    Greenhall C A 2003 Metrologia 40 335

    [7]

    Davis J A, Greenhall C A, Stacey P W 2005 Metrologia 42 1

    [8]

    Greenhall C A 2006 Metrologia 43 311

    [9]

    Suess M, Greenhall C A 2012 Metrologia 49 588

    [10]

    Parisi F, Panfilo G 2016 Metrologia 53 1185

    [11]

    伍贻威, 牟卫华, 龚航, 朱祥维, 欧钢 2016 武汉大学学报 信息科学版 41 1253

    Wu Y W, Mu W H, Gong H, Zhu X W, Ou G 2016 Geomat. Inf. Sci. Wuhan Univ. 41 1253

    [12]

    宋会杰, 董绍武, 王燕平, 安卫, 侯娟 2019 武汉大学学报. 信息科学版 44 1205

    Song H J, Dong S W, Wang Y P, An W, Hou J 2019 Geomat. Inf. Sci. Wuhan Univ. 44 1205

    [13]

    宋会杰, 董绍武, 屈俐俐, 王翔, 广伟 2017 仪器仪表学报 38 1809

    Song H J, Dong S W, Qu L L, Wang X, Guang W 2017 Chin. J. Sci. Instrum. 38 1809

    [14]

    Chaffee J W 1987 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 34 655

    [15]

    Brown K R Jr 1991 Proceeding of the 4th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 1991) Albuquerque, NM, USA, September 11−13, 1991 p223

    [16]

    杨元喜 2006 自适应动态导航定位 (第1版) (北京: 测绘出版社) 第189页

    Yang Y X 2006 Adaptive Navigation and Kinematic Positioning (1st Ed.) (Beijing: Surveying and Mapping Press) p189 (in Chinese)

    [17]

    Yang Y, He H, Xu G 2001 J. Geodesy 72 109

    [18]

    Yang Y, Song L, Xu T 2002 J. Geodesy 76 353

    [19]

    Cui X, Yang Y 2006 Proc. Natl. Acad. Sci. U.S.A. 16 846

    [20]

    Panfilo G, Harmegnies A, Tisserand L 2014 Metrologia 51 285

    [21]

    Levine J 1999 Rev. Sci. Instrum. 70 2567

    [22]

    Zucca C, Tavella P 2005 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 289

    [23]

    Galleani L, Sacerdote L, Tavella P, Zucca C 2003 Metrologia 40 S257

    [24]

    Suess M, Matsakis D, Greenhall C A 2010 42nd Annual Precise Time and Time Interval Meeting Reston, Virginia, U.S, November 15−18, 2010 p481

  • 图 1  基于两种算法模拟的AHM1时间偏差改正值 (a) 基于Kalman滤波算法模拟的AHM1的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM1时间偏差改正值

    Fig. 1.  The corrected time deviations of modeling AHM1 based on two algorithms: (a) The corrected time deviations of modelling AHM1 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM1 based on modified Kalman filter algorithm.

    图 2  基于两种算法的模拟AHM2的时间偏差改正值 (a) 基于Kalman滤波算法模拟的AHM2的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM2的时间偏差改正值

    Fig. 2.  The corrected time deviations of modelling AHM2 based on two algorithms: (a) The corrected time deviations of modelling AHM2 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM2 based on modified Kalman filter algorithm.

    图 3  基于两种Kalman滤波算法的时间尺度

    Fig. 3.  Time scale based on two Kalman filter algorithms.

    图 4  基于两种Kalman滤波算法的时间尺度的重叠Allan偏差

    Fig. 4.  The overlapping Allan deviation of time scale based on two Kalman filter algorithms.

    图 5  基于两种Kalman滤波算法的时间尺度

    Fig. 5.  Time scales based on two Kalman filter algorithms.

    图 6  基于两种Kalman滤波算法的时间尺度的稳定度

    Fig. 6.  Time scale stabilities based on two Kalman filter algorithms.

    表 1  基于两种Kalman滤波算法的时间尺度的重叠Allan偏差

    Table 1.  The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

    取样时间/102 sKF scale/10–14R-KF scale/10–14
    39.8513.80
    66.738.75
    124.385.20
    302.502.70
    601.761.85
    1201.221.24
    3000.850.73
    6000.880.57
    12000.890.63
    下载: 导出CSV

    表 2  原子钟的相对权重

    Table 2.  Relative weights of atomic clocks.

    原子钟编号HM4926HM4967HM0296Cs3436Cs2098
    相对权重0.300.300.140.190.07
    下载: 导出CSV

    表 3  原子钟的噪声强度

    Table 3.  Noise intensity of atomic clocks.

    原子钟
    编号
    q1WFM/10–36 [s2/s]q2RWFM/10–50 [s2/s3]q3RWFM/[s2/s5]
    HM49264.644.994.27 × 10–60
    HM49675.884.603.87 × 10–60
    HM02962.014.106.80 × 10–49
    Cs343626.9000
    Cs20985.4300
    下载: 导出CSV

    表 4  基于两种Kalman滤波算法的统计值

    Table 4.  Statistical values based on two Kalman filter algorithms.

    计算方法最大偏
    差值/ns
    最小值偏
    差值/ns
    平均偏
    差值/ns
    KF scale4.36–3.130.98
    R-KF scale1.11–2.65–0.47
    下载: 导出CSV

    表 5  基于两种Kalman滤波算法的时间尺度的重叠阿伦偏差值

    Table 5.  The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

    Sample time/103 sKF scale/10–14R-KF scale/10–14
    3.601.941.64
    7.201.341.17
    14.401.030.92
    36.000.920.85
    72.000.740.69
    144.000.530.51
    360.000.270.25
    下载: 导出CSV
  • [1]

    陈卫东, 刘要龙, 朱奇光, 陈颖 2013 物理学报 62 170506

    Chen W D, Liu Y L, Zhu Q G, Chen Y 2013 Acta Phys. Sin. 62 170506

    [2]

    刘洋洋, 廉保旺, 赵宏伟, 刘亚擎 2014 物理学报 63 228402

    Liu Y Y, Lian B W, Zhao H W, Zhao H W, Liu Y Q 2014 Acta Phys. Sin. 63 228402

    [3]

    林旭, 罗志才 2015 物理学报 64 080201

    Lin X, Luo Z C 2015 Acta Phys. Sin. 64 080201

    [4]

    赵龙 2012 物理学报 61 104301

    Zhao L 2012 Acta Phys. Sin. 61 104301

    [5]

    郭海荣, 杨元喜, 何海波, 徐天河 2010 测绘学报 39 146

    Guo H R, Yang Y X, He H B, Xu T H 2010 Acta Geod. Cartogr. Sin. 39 146

    [6]

    Greenhall C A 2003 Metrologia 40 335

    [7]

    Davis J A, Greenhall C A, Stacey P W 2005 Metrologia 42 1

    [8]

    Greenhall C A 2006 Metrologia 43 311

    [9]

    Suess M, Greenhall C A 2012 Metrologia 49 588

    [10]

    Parisi F, Panfilo G 2016 Metrologia 53 1185

    [11]

    伍贻威, 牟卫华, 龚航, 朱祥维, 欧钢 2016 武汉大学学报 信息科学版 41 1253

    Wu Y W, Mu W H, Gong H, Zhu X W, Ou G 2016 Geomat. Inf. Sci. Wuhan Univ. 41 1253

    [12]

    宋会杰, 董绍武, 王燕平, 安卫, 侯娟 2019 武汉大学学报. 信息科学版 44 1205

    Song H J, Dong S W, Wang Y P, An W, Hou J 2019 Geomat. Inf. Sci. Wuhan Univ. 44 1205

    [13]

    宋会杰, 董绍武, 屈俐俐, 王翔, 广伟 2017 仪器仪表学报 38 1809

    Song H J, Dong S W, Qu L L, Wang X, Guang W 2017 Chin. J. Sci. Instrum. 38 1809

    [14]

    Chaffee J W 1987 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 34 655

    [15]

    Brown K R Jr 1991 Proceeding of the 4th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 1991) Albuquerque, NM, USA, September 11−13, 1991 p223

    [16]

    杨元喜 2006 自适应动态导航定位 (第1版) (北京: 测绘出版社) 第189页

    Yang Y X 2006 Adaptive Navigation and Kinematic Positioning (1st Ed.) (Beijing: Surveying and Mapping Press) p189 (in Chinese)

    [17]

    Yang Y, He H, Xu G 2001 J. Geodesy 72 109

    [18]

    Yang Y, Song L, Xu T 2002 J. Geodesy 76 353

    [19]

    Cui X, Yang Y 2006 Proc. Natl. Acad. Sci. U.S.A. 16 846

    [20]

    Panfilo G, Harmegnies A, Tisserand L 2014 Metrologia 51 285

    [21]

    Levine J 1999 Rev. Sci. Instrum. 70 2567

    [22]

    Zucca C, Tavella P 2005 IEEE Trans. Ultrason. Ferroelectr. Freq. Control 52 289

    [23]

    Galleani L, Sacerdote L, Tavella P, Zucca C 2003 Metrologia 40 S257

    [24]

    Suess M, Matsakis D, Greenhall C A 2010 42nd Annual Precise Time and Time Interval Meeting Reston, Virginia, U.S, November 15−18, 2010 p481

  • [1] 林旭, 罗志才. 一种新的卫星钟差Kalman滤波噪声协方差估计方法. 物理学报, 2015, 64(8): 080201. doi: 10.7498/aps.64.080201
    [2] 樊养余, 李慧敏, 贾蒙. 基于自适应因子轨道延拓法的不变流形计算. 物理学报, 2010, 59(11): 7686-7692. doi: 10.7498/aps.59.7686
    [3] 刘洋洋, 廉保旺, 赵宏伟, 刘亚擎. Kalman滤波辅助的室内伪卫星相对定位算法. 物理学报, 2014, 63(22): 228402. doi: 10.7498/aps.63.228402
    [4] 林敏, 方利民. 双稳系统演化的时间尺度与随机共振的加强. 物理学报, 2009, 58(4): 2136-2140. doi: 10.7498/aps.58.2136
    [5] 龚志强, 支蓉, 郑志海, 周磊. 基于矩阵理论的全球温度资料的尺度性研究. 物理学报, 2009, 58(3): 2113-2120. doi: 10.7498/aps.58.2113
    [6] 李 勇, 毕勤胜. 连续搅拌槽式反应器中自催化化学反应的延迟同步. 物理学报, 2008, 57(10): 6099-6102. doi: 10.7498/aps.57.6099
    [7] 贾梦源, 赵刚, 侯佳佳, 谭巍, 邱晓东, 马维光, 张雷, 董磊, 尹王保, 肖连团, 贾锁堂. 双重频率锁定的腔衰荡吸收光谱技术及信号处理. 物理学报, 2016, 65(12): 128701. doi: 10.7498/aps.65.128701
    [8] 张家树, 肖先赐. 混沌时间序列的自适应高阶非线性滤波预测. 物理学报, 2000, 49(7): 1221-1227. doi: 10.7498/aps.49.1221
    [9] 肖先赐, 甘建超. 基于相空间邻域的混沌时间序列自适应预测滤波器(Ⅰ)线性自适应滤波. 物理学报, 2003, 52(5): 1096-1101. doi: 10.7498/aps.52.1096
    [10] 肖先赐, 甘建超. 基于相空间邻域的混沌时间序列自适应预测滤波器(Ⅱ)非线性自适应滤波. 物理学报, 2003, 52(5): 1102-1107. doi: 10.7498/aps.52.1102
    [11] 吴振森, 柯熙政. 原子钟噪声中的混沌现象及其统计特性. 物理学报, 1998, 47(9): 1436-1449. doi: 10.7498/aps.47.1436
    [12] 张家树, 肖先赐. 用一种少参数非线性自适应滤波器自适应预测低维混沌时间序列. 物理学报, 2000, 49(12): 2333-2339. doi: 10.7498/aps.49.2333
    [13] 韦保林, 罗晓曙, 郭维, 傅金阶, 汪秉宏, 全宏俊. 一种基于三阶Volterra滤波器的混沌时间序列自适应预测方法. 物理学报, 2002, 51(10): 2205-2210. doi: 10.7498/aps.51.2205
    [14] 张家树, 肖先赐. 用于混沌时间序列自适应预测的一种少参数二阶Volterra滤波器. 物理学报, 2001, 50(7): 1248-1254. doi: 10.7498/aps.50.1248
    [15] 汝鸿羽, 齐 亮, 季敏标, 谢飞翔, 聂瑞娟, 马 平, 戴远东, 王福仁, 白 云, 刘新元, 何定武, 赵 巍. 在SQUID心磁测量中基于奇异值分解和自适应滤波的噪声消除法. 物理学报, 2006, 55(5): 2651-2656. doi: 10.7498/aps.55.2651
    [16] 张家树, 肖先赐. 混沌时间序列的Volterra自适应预测. 物理学报, 2000, 49(3): 403-408. doi: 10.7498/aps.49.403
    [17] 戴远东, 王福仁, 李壮志, 马 平, 谢飞翔, 杨 涛, 聂瑞娟, 刘新元, 谢柏青. 射频SQUID心磁图数据自适应滤波研究. 物理学报, 2005, 54(4): 1937-1942. doi: 10.7498/aps.54.1937
    [18] 王梦蛟, 周泽权, 李志军, 曾以成. 混沌信号自适应协同滤波去噪. 物理学报, 2018, 67(6): 060501. doi: 10.7498/aps.67.20172470
    [19] 李 磊, 沈永平, 张立杰. 基于多时间尺度分析的青藏铁路沿线土壤热流研究. 物理学报, 2005, 54(4): 1958-1964. doi: 10.7498/aps.54.1958
    [20] 盛峥. 电离层电子总含量不同时间尺度的预报模型研究. 物理学报, 2012, 61(21): 219401. doi: 10.7498/aps.61.219401
  • 引用本文:
    Citation:
计量
  • 文章访问数:  364
  • PDF下载量:  18
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-12-18
  • 修回日期:  2020-05-11
  • 上网日期:  2020-06-01
  • 刊出日期:  2020-09-05

原子钟噪声变化时改进的Kalman滤波时间尺度算法

  • 1. 中国科学院国家授时中心, 西安 710600
  • 2. 中国科学院时间频率基准重点实验室, 西安 710600
  • 3. 中国科学院大学天文与空间科学学院, 北京 100049
  • 通信作者: 宋会杰, songhuijie@ntsc.ac.cn

摘要: Kalman滤波时间尺度算法是一种实时的原子钟状态估计方法, 在时间保持工作中具有重要的实用价值. 本文引入自适应因子改进Kalman滤波时间尺度算法, 对Kalman滤波时间尺度算法中状态预测协方差矩阵引入自适应因子, 利用统计量实时计算自适应因子的量值, 控制状态预测协方差矩阵的增长, 降低原子钟状态估计的扰动, 从而提高时间尺度的准确度和稳定度. 模拟数据和实测数据表明, 原子钟噪声强度变化时或噪声强度估计存在误差时改进的Kalman滤波时间尺度算法有效地提高了时间尺度的准确度和长期稳定度.

English Abstract

    • Kalman滤波作为一种经典的估计技术, 广泛应用于控制理论、信号分析和卫星导航星载原子钟的异常检测等领域[1,2]. 作为最小二乘方法的动态形式, Kalman滤波是一种实时算法, 通过一组观测序列及系统的动力学模型来求解状态向量估值. 这个估计值是通过混合两种信息得到的, 分别为基于动力学模型的预测信息和基于测量方程的新息向量. 最后通过对这两部分信息进行加权平均得到Kalman滤波估计. 当动力学模型噪声和测量方程噪声服从高斯分布时, 那么在最小二乘意义下, 估计是线性最优的, 在其他情况下, Kalman滤波算法是线性次优的. 除此之外, 为了保证估计的无偏性, 每一种噪声均值为0. 在时间尺度计算中, 要对原子钟的噪声强度进行估计. 文献[3-5]研究了Kalman滤波算法中噪声模型及协方差矩阵的确定方法. 文献[6-11]研究了基于Kalman滤波的时间尺度算法, 其中包含原子钟噪声模型的建立. 文献[6]基于原子钟模型研究了三类Kalman滤波时间尺度, 并分析了降低协方差算法的最优性. 其中Kalman加权时间尺度算法需要考虑原子钟的基本时间尺度方程(BTSE), 考虑了原子钟的频率状态估计和频漂状态估计, 权重算法考虑了时间尺度的短期稳定度. 本文通过Kalman滤波估计原子钟相位状态, 并根据“可预测性”权重算法计算时间尺度, 保证了时间尺度长期稳定度前提下, 改善其短期稳定度. 文献[7]基于氢原子钟研究了噪声类型主要为调频闪变噪声的一种Kalman滤波时间尺度算法. Kalman滤波算法中调频闪变噪声模拟为马氏噪声过程的线性组合. 本文考虑原子钟噪声类型主要为调频白噪声、调频随机游走噪声和频漂随机游走噪声. 文献[7]通过对状态估计误差协方差矩阵进行改正, 提高了Kalman滤波时间尺度的短期稳定度. 本文是在保证时间尺度长期稳定度的基础上, 改进时间尺度的短期稳定度. 文献[8]基于Kalman滤波算法研究了远距离不同实验室比对噪声情况下的时间尺度算法, 是对同一实验室无比对噪声情况的一种扩展. 文献[9]研究了Kalman滤波时间尺度算法中组合Brown和Greenhall 的协方差降低运算, 进一步优化了时间尺度计算结果, 提高了时间尺度的短期稳定度. 文献[10]探讨了将Kalman滤波算法用于计算协调世界时(coordinated universal time, UTC), 并与目前原子时尺度算法做了比较, 表示Kalman时间尺度具有较高的短期稳定度. 文献[11]研究了使用Kalman滤波器调整预测值的时间尺度算法. 文中建立原子钟的两状态方程, 通过Kalman滤波器估计原子钟的状态, 根据状态预测值计算时间尺度. 本文基于原子钟的三状态模型, 通过Kalman滤波器估计原子钟状态, 并根据原子钟的“可预测性”权重算法计算时间尺度, 保证时间尺度长期稳定度, 同时根据Kalman滤波算法改进时间尺度的短期稳定度. 以上文献研究的前提是原子钟噪声强度不变并且噪声强度能够准确估计, 但是外界环境影响与原子钟长期运行导致其性能下降, 原子钟噪声可能发生变化或噪声强度不能够准确估计, 这时根据Kalman滤波时间尺度算法产生的时间尺度性能降低. 文献[12,13]分别是基于渐消因子的改进Kalman滤波时间尺度估计算法研究和Sage窗的自适应Kalman滤波用于钟差预报的研究. 文献[12]对状态预测协方差矩阵引入渐消因子, 降低原子钟状态估计的扰动, 通过原子钟的时差状态和时差状态估计之差定义时间尺度, 提高了时间尺度的短期稳定度. 而本文研究自适应因子改进Kalman滤波状态估计, 原子钟的状态估计作为原子钟的预测部分进行扣除, 基于“可预测性”取权方法计算时间尺度. 这样既保证了时间尺度的长期稳定度, 同时又改善了时间尺度的短期稳定度. 文献[13]研究了基于Sage窗的自适应Kalman滤波用于钟差预报, 基于Sage窗构造自适应因子, 取决于窗口内各历元状态参数与当前状态的适应程度, 若所开的窗口内的状态参数改正量正好反映载体当前状态扰动水平, 则求得的状态方程协方差矩阵能调节状态方程的噪声水平; 反之由于对历史信息的平滑, 反而使状态参数偏离. 本文时间尺度算法采用预测残差构造自适应因子, 能够较好地反映模型变化, 但对于稳定状态, 基于Sage窗的自适应Kalman滤波具有更好的效果. 鉴于此, 本文研究了噪声强度变化时改进的Kalman滤波时间尺度算法, 通过模拟数据和实测数据验证并分析了该算法的性能.

      Kalman 滤波时间尺度算法需要确定原子钟状态模型噪声, 通常情况下模型噪声由调频白噪声、调频随机游走噪声和频漂随机游走噪声组成, 并且假定原子钟噪声强度恒定, 但是受环境因素或原子钟老化影响, 原子钟噪声强度会发生改变. Kalman滤波时间尺度算法需要估计噪声强度, 估计方法的准确与否会影响时间尺度的性能, 通常根据原子钟测得的Allan 方差和Hadamard 方差来估计原子钟噪声强度, 但是受测量条件的影响, 原子钟各类方差估计存在误差, 因此原子钟噪声强度估计存在偏差, 最终影响时间尺度的性能. 文中改进Kalman滤波时间尺度算法的主要思想是原子钟噪声强度变化或噪声强度估计存在误差时, 计算自适应因子控制Kalman滤波预测协方差矩阵的增长, 降低状态模型信息不准确对状态估计的影响. 对于原子钟噪声强度恒定与噪声强度准确估计时, 文中改进的Kalman滤波时间尺度算法变为经典Kalman滤波时间尺度算法.

      文章组织如下: 第2节建立基于Kalman滤波算法的原子钟组模型; 第3节推导改进的Kalman滤波状态估计; 第4节介绍“可预测性”取权方法; 第5节为时间尺度的计算及分析; 第6节为本文的结论.

    • 在具有$N$台原子钟的钟组中, 状态向量包含每一台原子钟的相位偏差, 频率偏差和频率漂移偏差信息. 定义${t_k}$时刻状态向量的顺序依次为相位偏差、频率偏差和频率漂移偏差

      $\begin{split} & X(t_k) = \\ & [x_1(t_k)~y_1(t_k)~z_1(t_k), \cdots, x_N (t_k) y_N(t_k) z_N(t_k)], \end{split} $

      式中, ${x_i}\left( t \right)$, ${y_i}\left( t \right)$${z_i}\left( t \right)$分别表示钟$i$$t$时刻的相位偏差、频率偏差和频率漂移偏差. 第$i$台钟从时刻${t_{k - 1}}$到时刻${t_k}$的状态由下面方程描述[6]:

      $\begin{split} x_i(t_k)={}& x_i(t_{k -1}) + y_i(t_{k - 1}) \cdot \tau \\ &+ \frac{1}{2} \cdot {\tau ^2} \cdot {z_i} (t_{k - 1}) + w_{xi}(t_k), \end{split} $

      ${y_i}\left( {{t_k}} \right) = {y_i}\left( {{t_{k - 1}}} \right) + {z_i}\left( {{t_{k - 1}}} \right) \cdot \tau + {w_{yi}}\left( {{t_k}} \right), $

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

      (2)式和(3)式中$\tau = {t_k} - {t_{k - 1}}$, 表示相邻时刻的时间间隔; ${w_{xi}}\left( {{t_k}} \right)$, ${w_{yi}}\left( {{t_k}} \right)$${w_{zi}}\left( {{t_k}} \right)$分别表示${t_k}$时刻的调频白噪声、调频随机游走噪声和频漂随机游走噪声; ${{{W}}_i}({t_k}) = {\left[ {{w_{xi}}({t_k}), {w_{yi}}({t_k}), {w_{zi}}({t_k})} \right]^{\rm T}}$表示原子钟状态估计的过程噪声向量, ${{{W}}_i} \left( {{t_k}} \right)\sim N\left( {0, {{{Q}}_i}} \right)$, ${{{Q}}_i}$表示第$i$台钟的过程噪声协方差矩阵,

      $ \begin{split} & {{{Q}}_{i}} (\tau ) = \\ & \left[ {\begin{array}{*{20}{c}} {{q_{xi}}\tau + {q_{yi}}\dfrac{{{\tau ^3}}}{3} + {q_{zi}}\dfrac{{{\tau ^5}}}{{20}}}&{{q_{yi}}\dfrac{{{\tau ^2}}}{2} + {q_{zi}}\dfrac{{{\tau ^4}}}{8}}&{{q_{zi}}\dfrac{{{\tau ^3}}}{6}} \\ {{q_{yi}}\dfrac{{{\tau ^2}}}{2} + {q_{zi}}\dfrac{{{\tau ^4}}}{8}}&{{q_{yi}}\tau + {q_{zi}}\dfrac{{{\tau ^3}}}{3}}&{{q_{zi}}\dfrac{{{\tau ^2}}}{2}} \\ {{q_{zi}}\dfrac{{{\tau ^3}}}{6}}&{{q_{zi}}\dfrac{{{\tau ^2}}}{2}}&{{q_{zi}}\tau } \end{array}} \right] . \end{split} $

      (5)式中${q_{xi}}$, ${q_{yi}}$${q_{zi}}$分别表示钟$i$的调频白噪声、调频随机游走噪声和频漂随机游走噪声的强度, 与Hadamard方差的关系为[14]

      $H\sigma _{yi}^2\left( \tau \right) = \frac{{{q_{xi}}}}{\tau } + \frac{{{q_{yi}}\tau }}{6} + \frac{{11{q_{zi}}{\tau ^3}}}{{120}}.$

      假设以钟1为参考钟, ${t_k}$时刻钟$i$的钟差测量表示为

      ${z_{i1}}\left( {{t_k}} \right) = {x_i}\left( {{t_k}} \right) - {x_1}\left( {{t_k}} \right) + {v_{1i}}\left( {{t_k}} \right).$

      为了利用Kalman滤波算法, 需要建立原子钟组的矩阵向量方程, 其中状态向量方程表达为[15]

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

      文献[15]说明GPS是一个由原子钟组成的网络, 早期的GPS设计将其中的一台原子钟视为“主钟”, 当“主钟”出现故障时会影响系统的性能, 因此这种设计过多依赖于主钟. 提出了一种称为GPS组合时钟的解决方案, 本质上它是一个Kalman滤波应用到原始的轨道/时钟问题, 包括所有时钟, 没有做主钟的选择, 提高了系统的可靠性. (8)式中${{\varPhi}} \left( \tau \right)$为状态转移矩阵, 有对角块形式

      ${{\varPhi}} \left( \tau \right) = \left[ {\begin{array}{*{20}{c}} 1&\tau &{{{{\tau ^2}}}/{2}} \\ 0&1&\tau \\ 0&0&1 \end{array}} \right].$

      ${t_k}$时刻过程噪声向量${{W}}({{t_k}}) \!=\! [ {{{W}}_1}( {{t_k}}), \cdots\!, {{{W}}_N} (t_k) ]^{\rm{T}}$,其中${\rm{E}}{{W}}\left( {{{{t}}_k}} \right) \!=\! 0$, ${\rm{E}}{{W}}({{t_k}}){{W}}{({{t_j}})^{\rm{T}}} \!=\! {{Q}}( \tau ) {\delta _{kj}}$, 协方差矩阵${{Q}}\left( \tau \right)$是由${{{Q}}_{i}}\left( \tau \right)$组成的对角块矩阵.

      测量方程表达式为

      ${{Z}}\left( {{t_k}} \right) = {{H}}\left( {{t_k}} \right){{X}}\left( {{t_k}} \right) + {{V}}\left( {{t_k}} \right).$

      (9)式中${{Z}}\left( {{t_k}} \right) = {\left[ {{z_{21}}\left( {{t_k}} \right), L, {z_{N1}}\left( {{t_k}} \right)} \right]^{\rm T}}$, ${{H}}\left( {{t_k}} \right)$$\left( {N - 1} \right) \times 3 N$矩阵. 对于三台钟的情况:

      ${{H}}\left( {{t_k}} \right) = \left[ {\begin{array}{*{20}{c}} { - 1}&0&0&1&0&0&0&0&0 \\ { - 1}&0&0&0&0&0&1&0&0 \end{array}} \right].$

      测量噪声向量${{V}}\left( {{t_k}} \right) = {\left[ {{v_{12}}\left( {{t_k}} \right), \cdots, {v_{1 N}}\left( {{t_k}} \right)} \right]^{\rm T}}$, 其中${\rm{E}}{{V}}\left( {{t_k}} \right) = 0$, ${\rm{E}}{{V}}\left( {{t_k}} \right){{V}}{\left( {{t_j}} \right)^{\rm{T}}} = {\varSigma _k}{\delta _{kj}}$, ${\varSigma _k}$是测量噪声方差矩阵, 并且${\rm{E}}{{W}}\left( {{t_k}} \right){{V}}{\left( {{t_j}} \right)^{\rm{T}}} = 0$.

    • 改进Kalman滤波时间尺度算法与传统Kalman滤波时间尺度算法步骤基本相同, 分为两部分, 包括利用原子钟状态方程得到的状态预测估计和利用测量数据对状态预测估计的修正, 即状态验后估计, 对于时刻${t_k}$, 具体计算步骤如下.

    • 1)预测. 利用状态转移矩阵${{\varPhi }}\left( \tau \right)$和原子钟模型噪声${{W}}\left( {{t_k}} \right)$计算得到${t_k}$时刻原子钟状态的线性最小方差预测, 表示为

      $\bar {{X}}\left( {{t_k}} \right) = {{\varPhi}} \left( \tau \right) \cdot \hat {{X}}\left( {{t_{k - 1}}} \right), $

      式中, $\bar {{X}}\left( {{t_k}} \right)$表示${t_k}$时刻原子钟的预测状态估计, $\hat {{X}}\left( {{t_{k - 1}}} \right)$表示${t_{k - 1}}$时刻的原子钟的状态估计, 预测协方差矩阵表示为

      ${\varSigma _{{{{\bar X}}_k}}} = {{\varPhi}} \left( \tau \right) \cdot {\varSigma _{{{{\hat X}}_{k - 1}}}} \cdot {{{\varPhi}} ^{\rm T}}\left( \tau \right) +{{ Q}}\left( \tau \right), $

      式中, ${\varSigma _{{{{{\bar X}}}_{{k}}}}}$表示${t_k}$时刻原子钟的预测状态协方差矩阵, ${\varSigma _{{{{{\hat X}}}_{k - 1}}}}$表示${t_{k - 1}}$时刻原子钟的状态估计协方差矩阵.

      2)测量. 通过${t_k}$时刻原子钟测量数据${{Z}}\left( {{{{t}}_k}} \right)$估计新息向量, 新息向量表示为测量值与预测值之差, 表达式为

      ${{i}}\left( {{t_k}} \right) = {{Z}}\left( {{t_k}} \right) - {{H}}\left( {{t_k}} \right) \cdot \bar {{X}}\left( {{t_k}} \right).$

      新息向量的协方差矩阵表达式为

      ${{C}}\left( {{t_k}} \right) = {{H}}\left( {{t_k}} \right) \cdot {\varSigma _{{{{{\bar X}}}_{{k}}}}} \cdot {{H}}{\left( {{t_k}} \right)^{\rm{T}}} + {\varSigma _k}, $

      式中, ${\varSigma _k}$表示${t_k}$时刻测量值${{Z}}\left( {{t_k}} \right)$的误差协方差矩阵. 改进Kalman滤波时间尺度算法的增益矩阵表示为[16]

      $\begin{split}{\bar {{K}}_k} =& \dfrac{1}{{{\alpha _k}}} \cdot {\varSigma _{{{\bar X}_k}}} \cdot {{{H}}^{\rm T}}\left( {{t_k}} \right) \! \Big(\frac{1}{{{\alpha _k}}} \cdot {{H}}({{t_k}}) \\ & \times {\varSigma _{{{{{\bar X}}}_k}}} \cdot {{{H}}^{\rm T}}({{t_k}}) \!+\! {\varSigma _k} \!\Big)^{-1}, \end{split} $

      式中, ${\alpha _k}$为自适应因子, 用于原子钟状态预测信息的自适应估计. 文献[16]表明运动定位与导航一般应用Kalman滤波算法, 可靠的Kalman滤波算法要求有可靠的函数模型、随机模型以及合理的估计方法. 然而精确的函数模型的构造十分困难, 随机模型先验信息的获取一般都是基于验前统计信息, 而任何统计信息都有可能失真. 本文围绕如何利用当前观测信息和状态估计更新先验信息, 提高状态估计的准确性展开研究. 在Kalman滤波中, 这种由估计过程自适应地调整、更新先验信息的算法称为自适应Kalman滤波.

      根据以上信息得到${t_k}$时刻原子钟的状态估计, 表达式为

      $\hat {{X}}\!\left( {{t_k}} \right) \!=\! \bar {{X}}\!\left( {{t_k}} \right)\! +\! {\bar K_k} \cdot\! \left[ {{{Z}}\!\left( {{t_k}} \right) \!-\! {{H}}\left( {{t_k}} \right) \cdot \bar {{X}}\!\left( {{t_k}} \right)} \right].$

      ${t_k}$时刻原子钟状态估计的协方差矩阵表达式为

      ${\varSigma _{{{\hat X}}\left( {{t_k}} \right)}} = \left[ {{{I}} - {{\bar {{K}}}_k} \cdot {{H}}\left( {{t_k}} \right)} \right] \cdot {\varSigma _{{{{{\bar X}}}_{{k}}}}}, $

      式中, I表示单位矩阵, (14)式对预测状态协方差矩阵添加自适应因子${\alpha _k}$, 自适应因子用于控制模型误差对原子钟状态估计的影响, 通过构造模型误差判别统计量, 判断统计量与给定阈值的大小建立自适应因子函数, 当模型误差判别统计量小于给定阈值时, 自适应因子${\alpha _k}=1$, 改进的Kalman滤波算法变为传统Kalman滤波算法. 当模型误差判别统计量大于给定阈值时, 自适应因子增大, 模型预测状态协方差矩阵降低, 从而减小模型误差对原子钟状态估计的影响.

    • 改进Kalman滤波算法中自适应因子的确定, 对于原子钟状态估计, 由于测量信息不足以解算当前历元的状态参数, 故不能采用状态参数不符值构造统计量, 由于预测残差向量能够较好地反应模型变化, 考虑采用预测残差构造判别统计量, 利用三段函数构造自适应因子, 自适应因子表达式为[17-19]

      $ {\alpha _k}\! = \!\left\{ {\begin{array}{*{20}{l}} {1,} \\ {\dfrac{{{c_0}}}{{| {{{\tilde V}_k}}|}} \cdot {{\left( {\dfrac{{{c_1} - | {{{\tilde V}_k}}|}}{{{c_1} - {c_0}}}} \right)}^2}\!,} \\ {0,} \end{array}} \right.\;\; \begin{array}{*{20}{l}} {| {{{\tilde V}_k}}| \leqslant {c_0},} \\ {{c_0} < | {{{\tilde V}_k}}| \leqslant {c_1},} \\ {| {{{\tilde V}_k}} | > {c_1},} \end{array} $

      式中, ${c_0}$可取为1.0—1.5; ${c_1}$可取为3.0—8.5.

      ${\tilde {{V}}_k} = {{\left\| {{{\bar {{V}}}_k}} \right\|}\Big/ {\sqrt {{\rm{tr}}\left( {{\varSigma _{{{{{\bar V}}}_k}}}} \right)} }}, $

      式中, ${\bar {{V}}_k}$表示${t_k}$时刻预测残差向量, ${\rm{tr}}\left( {{\varSigma _{{{{{\bar V}}}_k}}}} \right)$表示预测残差向量协方差矩阵的迹.

    • 原子钟权重算法的原则为一台好的钟是一台“可预测”的钟[20,21], 主要思想是原子钟的预测频率与实际频率的差异在每个月的计算区间内进行估计, 使用这些值来定义权重, 并适当设定权重的上限. 算法进行4次迭代, 每次迭代运行如下.

      1) 对于值$\left[ {{{\bar h}_0} - {h_i}} \right]$是使用一组相对权重计算得到的, 式中${\bar h_0}$表示时间尺度, ${h_i}$表示第$i$台钟的读数. 第一次迭代的权重使用最近区间归一化后的权重值. 在接下来的迭代中, 权重根据上一次迭代计算得到.

      2) 频率值${y_{i, {I_k}}}$与预测频率${\hat y_{i, {I_k}}}$之差的绝对值表示为

      ${\varepsilon _{i,{I_k}}} = \left| {{y_{i,{I_k}}} - {{\hat y}_{i,{I_k}}}} \right|, $

      式中, $i$表示时钟的索引, ${I_k}$表示时间间隔, ${y_{i, {I_k}}}$表示频率实测值, ${\hat y_{i, {I_k}}}$表示频率预测值.

      3) 考虑到新的测量数据具有更可靠的统计性质, 通过滤波器使得较新的测量数据具有更重要的作用, 计算

      $\sigma _i^2 = \frac{{\displaystyle\sum\limits_{j = 1}^M {\left( {\frac{{M + 1 - j}}{M}} \right)\varepsilon _{i,j}^2} }}{{\displaystyle\sum\limits_{j = 1}^M {\left( {\frac{{M + 1 - j}}{M}} \right)} }}, $

      其中$i$表示钟; $j$表示计算区间; $M$为计算区间的数目.

      4) 计算时钟${H_i}$的权重为

      ${\omega _{i,{\rm{temp}}}} = \frac{{{1/ {\sigma _i^2}}}}{{\displaystyle\sum\nolimits_{i = 1}^N {{1/ {\sigma _i^2}}} }}.$

      除了时钟${H_i}$的权重大于设置的最大权和时钟${H_i}$在计算区间出现异常外, 时钟${H_i}$的相对权重${\omega _i}$$={\omega _{i, {\rm{temp}}}}$.

    • 为了检验噪声强度变化时改进Kalman滤波时间尺度算法的性能, 首先模拟原子钟的数据, 通过改变原子钟运行过程中的噪声强度, 模拟原子钟噪声强度变化的情况, 研究改进Kalman滤波时间尺度算法的性能, 最后根据实测数据验证噪声强度估计存在误差时算法的性能并分析.

    • 根据原子钟的数学模型模拟原子钟的数据[22,23], 原子钟系统随机微分方程的解能够表达为迭代形式:

      ${{X}}\left( {{t_{k + 1}}} \right) = {{\varPhi}} \left( \tau \right){{X}}\left( {{t_k}} \right) + {{W}}\left( {{t_{k + 1}}} \right).$

      对于单台原子钟, (22)式中${{X}}(t) = \left[ {x(t), y(t), z(t)} \right]$, $x(t)$表示模拟原子钟相对于参考原子钟的相位偏差, $y(t)$表示频率偏差的一部分, $z(t)$表示频率漂移偏差的一部分. ${{\varPhi}} \left( \tau \right)$表达式与(8)式相同, ${{W}}\left( {{t_k}} \right)$是原子钟的噪声向量; ${{Q}}\left( \tau \right)$为协方差矩阵, 与(8)式中的定义相同, ${q_{xi}}$, ${q_{yi}}$${q_{zi}}$与(5)式中表达的意义相同.

      产生钟的数据后, 计算钟差测量值, 钟差测量值为钟组中原子钟与某台参考钟的差值, 文中以基准钟为参考, 计算各台原子钟与参考钟的钟差, 通过测量协方差矩阵${\Sigma _k}$产生测量噪声加入到钟差, 计算方法如(9)式. 模拟中假设原子钟测量方差均为$1 \times {10^{ - 20}}$ s2, 且测量值互不相关. 根据原子钟系统模型, 钟的噪声强度和测量噪声对原子钟进行模拟, 文中模拟5台原子钟, 包括2台氢原子钟, 2台铯原子钟和1台基准钟[24].

      为了能够清楚的反映原子钟噪声强度变化的影响, 假设原子钟比对测量间隔为300 s, 即每间隔300 s测量一次, 测量周期为3000次. 当原子钟运行至1001次时改变原子钟的噪声强度, 取两台模拟氢原子钟为例, 改变模拟氢原子钟编号为AHM1和编号为AHM2的噪声强度, AHM1和AHM2的噪声强度分别为: ${q_{x1}}$= 5.8×10–26 (s2/s), ${q_{y2}}$= 5.1×10–34 (s2/s3), ${q_{z1}}$= 9.4×10–40 (s2/s5)和${q_{x2}}$= 5.9×10–26 (s2/s), ${q_{y2}}$= 6.2×10–34 (s2/s3), ${q_{z2}}$= 9.5×10–51 (s2/s5). 当原子钟运行至2001次时再次改变模拟氢原子钟的噪声强度, AHM1和AHM2的噪声强度分别为: ${q_{x1}}$= 5.8×10–26 (s2/s), ${q_{y1}}$= 9.1×10–34 (s2/s3), ${q_{z1}}$= 9.4×10–51 (s2/s5)和${q_{x2}}$= 5.9×10–26 (s2/s), ${q_{y2}}$= 9.2×10–34 (s2/s3), ${q_{z2}}$= 9.5×10–40 (s2/s5). 分别基于经典Kalman滤波时间尺度算法和改进Kalman滤波时间尺度算法估计原子钟的状态, 并计算各台原子钟的时间偏差改正值(corrected time deviations), 时间偏差改正值表达式为

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

      式中, ${\tilde x_i}({t_k})$表示${t_k}$时刻原子钟$i$的时间偏差改正值, ${\hat x_i}\left( {{t_k}} \right)$表示${t_k}$时刻原子钟$i$的时间偏差估计值. AHM1和AHM2的时间偏差改正值分别如图1图2所示. 计算表明, 基于经典Kalman滤波时间尺度算法, AHM1与AHM2的时间偏差改正值分别约有4 ns和3 ns的扰动, 而基于改进Kalman滤波时间尺度算法得到的时间偏差改正值不存在扰动现象, 说明基于经典Kalman滤波时间尺度算法得到的原子钟状态估计受原子钟噪声强度变化的影响, 而改进Kalman滤波时间尺度算法有效抵制了原子钟噪声强度变化对状态估计的影响, 提高了状态估计的准确度. 基于原子钟时间偏差改正值计算时间尺度, 计算公式为

      图  1  基于两种算法模拟的AHM1时间偏差改正值 (a) 基于Kalman滤波算法模拟的AHM1的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM1时间偏差改正值

      Figure 1.  The corrected time deviations of modeling AHM1 based on two algorithms: (a) The corrected time deviations of modelling AHM1 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM1 based on modified Kalman filter algorithm.

      图  2  基于两种算法的模拟AHM2的时间偏差改正值 (a) 基于Kalman滤波算法模拟的AHM2的时间偏差改正值; (b) 基于改进Kalman滤波算法模拟的AHM2的时间偏差改正值

      Figure 2.  The corrected time deviations of modelling AHM2 based on two algorithms: (a) The corrected time deviations of modelling AHM2 based on Kalman filter algorithm; (b) the corrected time deviations of modelling AHM2 based on modified Kalman filter algorithm.

      ${\bar x_0}\left( {{t_k}} \right) = \sum\limits_{i = 1}^N {{w_i}} {\tilde x_i}\left( {{t_k}} \right).$

      不妨假设${w_i}$取等权, 两种Kalman滤波算法的时间尺度如图3所示, 图3中经典Kalman滤波时间尺度用KF Scale表示, 改进Kalman滤波时间尺度用R-KF Scale表示, 图4表1中的表列与图3相同. 基于经典Kalman滤波算法的时间尺度受噪声强度变化的影响, 出现明显的扰动, 降低了时间尺度的准确度. 基于改进Kalman滤波算法的时间尺度有效抵制了噪声强度变化对时间尺度的影响, 提高了时间尺度的准确度.

      取样时间/102 sKF scale/10–14R-KF scale/10–14
      39.8513.80
      66.738.75
      124.385.20
      302.502.70
      601.761.85
      1201.221.24
      3000.850.73
      6000.880.57
      12000.890.63

      表 1  基于两种Kalman滤波算法的时间尺度的重叠Allan偏差

      Table 1.  The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

      图  3  基于两种Kalman滤波算法的时间尺度

      Figure 3.  Time scale based on two Kalman filter algorithms.

      图  4  基于两种Kalman滤波算法的时间尺度的重叠Allan偏差

      Figure 4.  The overlapping Allan deviation of time scale based on two Kalman filter algorithms.

      基于重叠阿伦偏差(overlapping Allan deviation)分析频率稳定度, 两种Kalman滤波算法的时间尺度的重叠阿伦偏差如图4表1所列, 分析得出相比于经典Kalman滤波算法, 基于改进Kalman滤波算法得到的时间尺度长期稳定度提高, 短期稳定度有所降低, 分析原因, 对Kalman滤波预测协方差矩阵引入自适应因子, 控制了噪声强度变化对时间尺度的影响, 降低了扰动, 提高了长期稳定度, 但同时控制作用影响了相位偏差改正值的短期稳定度, 进而影响了时间尺度的短期稳定度.

    • 不同于模拟钟情况, 对于原子钟实际测量数据定义的Kalman滤波时间尺度是非观测的. 比如, 不能得到单台原子钟的读数, 但是可以得到原子钟之间的钟差数据, 原子钟相对于参考钟或时间尺度的钟差数据. 因此, 只能通过计算得到Kalman滤波时间尺度与钟组中单台原子钟的钟差值. 原子钟实测数据的Kalman滤波时间尺度的计算式为

      $ \begin{split} {\bar h_0}\left( t \right) - {h_j}\left( t \right)& =\sum\limits_{i = 1}^N {{w_i}} \left( {{{\tilde h}_i}\left( t \right) - {h_j}\left( t \right)} \right)\\ &= \sum\limits_{i = 1}^N {{w_i}\left[ {{x_{ij}}(t) + {{\hat x}_i}(t)} \right]}, \end{split} $

      式中, ${\bar h_0}\left( t \right)$表示$t$时刻Kalman滤波时间尺度; ${h_j}\left( t \right)$表示$t$时刻原子钟$j$的读数; ${\tilde h_i}\left( t \right)$表示改正钟, 有关系式${\tilde h_i}\left( t \right) = {h_0}(t) - {\tilde x_i}\left( t \right)$, ${h_0}(t)$表示$t$时刻理想钟的读数; ${w_i}$表示第$i$台原子钟的权重; ${x_{ij}}(t)$表示$t$时刻钟$i$与钟$j$的钟差值.

      为了检验噪声强度估计存在误差时改进Kalman滤波时间尺度算法的性能, 随机选取中国科学院国家授时中心(national time service center, NTSC)时频基准实验室的5台原子钟, 在国际原子时(international atomic time, TAI)计算中的编号分别为HM4926, HM4967, HM0296, Cs3436和Cs2098. 钟差UTC(NTSC)-Clock(i)为用于时间尺度计算的数据形式, 数据长度为MJD58484 —MJD58726.96 (2019年 1月1日0时—8月31日23时), 数据采样间隔为1 h, 时间尺度根据第4节“可预测性”方法计算相对权重, 并且最大权设为0.3, 计算结果如表2所列, 得到HM4926与HM4967取得最大权重为0.3, 说明这2台氢原子钟的可预测性能优于其余3台原子钟, Cs2098取得最小权重0.07, 表明相对于钟组中其余原子钟, 其可预测性能低.

      原子钟编号HM4926HM4967HM0296Cs3436Cs2098
      相对权重0.300.300.140.190.07

      表 2  原子钟的相对权重

      Table 2.  Relative weights of atomic clocks.

      研究原子钟噪声强度估计存在误差时改进Kalman滤波时间尺度算法的性能. 随机改变原子钟的噪声强度, 如表3所列. 改变氢原子钟的噪声强度, 2台铯原子钟的调频随机游走噪声与频漂随机游走噪声强度为0. 分别基于经典Kalman滤波时间尺度算法和改进Kalman滤波算法计算时间尺度, 计算结果与UTC (NTSC) 进行比较, 如图5所示. 基于经典Kalman滤波算法的时间尺度明显偏离UTC (NTSC), 而基于改进Kalman滤波算法的时间尺度保持与UTC (NTSC)一致. 表4给出了通过两种Kalman滤波算法得到的时间尺度与UTC (NTSC)的最大偏差值, 最小偏差值与平均偏差值. 相比于经典Kalman滤波算法, 基于改进Kalman滤波算法得到的时间尺度与UTC (NTSC) 相比, 最大偏差值降低了3.25 ns, 最小偏差值改进了0.48 ns, 平均绝对偏差值改进了0.51 ns, 提高了时间尺度的准确度.

      原子钟
      编号
      q1WFM/10–36 [s2/s]q2RWFM/10–50 [s2/s3]q3RWFM/[s2/s5]
      HM49264.644.994.27 × 10–60
      HM49675.884.603.87 × 10–60
      HM02962.014.106.80 × 10–49
      Cs343626.9000
      Cs20985.4300

      表 3  原子钟的噪声强度

      Table 3.  Noise intensity of atomic clocks.

      计算方法最大偏
      差值/ns
      最小值偏
      差值/ns
      平均偏
      差值/ns
      KF scale4.36–3.130.98
      R-KF scale1.11–2.65–0.47

      表 4  基于两种Kalman滤波算法的统计值

      Table 4.  Statistical values based on two Kalman filter algorithms.

      图  5  基于两种Kalman滤波算法的时间尺度

      Figure 5.  Time scales based on two Kalman filter algorithms.

      基于重叠阿伦偏差估计时间尺度的频率稳定度, 分别计算基于两种Kalman滤波算法得到的时间尺度的重叠阿伦偏差. 因为系统中原子钟的采样时间为1 h, 因此重叠阿伦偏差的最小平滑时间也为1 h, 计算结果如图6表5所列. 相比于经典Kalman滤波算法, 基于改进Kalman滤波算法的时间尺度在相应平滑时间的重叠阿伦偏差有所降低, 说明稳定度有所提高. 分析其原因是受原子钟噪声强度估计误差的影响, 基于经典Kalman滤波的时间尺度的频率稳定度降低.

      图  6  基于两种Kalman滤波算法的时间尺度的稳定度

      Figure 6.  Time scale stabilities based on two Kalman filter algorithms.

      Sample time/103 sKF scale/10–14R-KF scale/10–14
      3.601.941.64
      7.201.341.17
      14.401.030.92
      36.000.920.85
      72.000.740.69
      144.000.530.51
      360.000.270.25

      表 5  基于两种Kalman滤波算法的时间尺度的重叠阿伦偏差值

      Table 5.  The overlapping Allan deviations of time scale based on two Kalman filter algorithms.

      原子钟噪声强度估计存在误差时, 基于经典Kalman滤波算法的时间尺度准确度降低, 噪声强度估计误差影响了时间尺度的准确度. 改进Kalman滤波算法对状态预测协方差矩阵引入自适应因子, 控制噪声强度估计误差引起的状态预测协方差矩阵的增长, 降低原子钟状态估计的扰动, 提高了时间尺度的准确度.

    • 针对原子钟运行过程中的噪声强度变化和噪声强度估计存在误差的情况, 本文研究了一种改进的Kalman滤波时间尺度算法, 通过原子钟模拟数据和实测数据验证了算法的性能并进行了分析. 该算法使用预测残差构造判别统计量, 并构造自适应因子, 通过对经典Kalman滤波预测协方差矩阵引入自适应因子, 实时控制原子钟噪声强度变化与噪声强度估计误差引起的状态估计扰动, 提高了时间尺度的准确度, 保证了时间尺度的长期稳定度. 同时结合模拟数据分析, 实时控制噪声强度变化引起的原子钟状态估计扰动过程也影响了时间尺度的短期稳定度. 进一步研究自适应因子的构造方法, 降低由于控制作用影响时间尺度的短期稳定性能是下一步的重点工作.

      模拟数据和实测数据的结果都计算了重叠Allan偏差, 如图4图6所示. 图4图6中的结果特征有一些差别, 分析原因, 与原子钟的取权算法有关, 文中实测数据采用“可预测性”取权算法是为了保证时间尺度的长期稳定度.

      根据时间尺度的计算过程可以得出, 基于经典Kalman滤波时间尺度的性能也与发生噪声强度变化原子钟的权重有关, 噪声强度变化的原子钟权重越大, 对时间尺度性能的影响越大.

参考文献 (24)

目录

    /

    返回文章
    返回