-
在卫星导航系统中, 系统时间的建立和维持处于基础和核心的地位。时间尺度的性能既取决于各台原子钟的性能, 也取决于所设计的时间尺度算法[1]。时间尺度算法用于综合钟组内所有的原子钟, 建立一个稳定度、可靠性更高的时间尺度。传统的时间尺度算法主要包括:国际权度局采用的加权平均算法(称之为ALGOS算法)[1-2], 美国国家标准与技术研究院(National Institute of Standards and Technolgoy, NIST)采用的加权平均算法(称之为AT1算法)[2-3], NIST在1980年代使用的标准Jones-Tyron型Kalman滤波器算法[4-6]。文献[6]研究表明NIST的Kalman滤波器算法实际上也是加权平均算法。任何由加权平均算法建立得到的时间尺度都可以由时间尺度基本方程[2](即本文中的式(2))来描述。时间尺度基本方程[2]表明时间尺度算法的核心是通过N-1组观测量(钟差)计算得到N台钟的权重和预测值[2]。目前对于时间尺度的研究重点为:(1)通过设置权重来提高时间尺度的稳定度。例如, 当权重反比于每台钟的中期稳定度时, 则优化的是时间尺度的中期稳定度; 而当权重反比于每台钟的长期稳定度时, 则优化的是时间尺度的长期稳定度。(2)通过选取预测值来抑制时间尺度在两个相邻的时间段内由于权重改变, 或是有原子钟加入或剔除出钟组时引起的时差和频差的跳变, 以此来保证时间尺度在时间和频率上的连续性[1-3]。
NIST的Kalman滤波器算法的不足[6-8]为:(1)该算法主要优化了时间尺度的长期稳定度, 无法有效提升时间尺度的中短期稳定度; (2)由于系统不完全可观造成Kalman滤波器发散。后人对该算法进行了一些改进[6-8]。同时, 文献[8]再次详细分析了NIST的Kalman滤波器算法。不同于这些传统Kalman滤波器算法, 本文提出了一种新的时间尺度算法, 主要关注如何实时调整预测值, 而不是权重来提高时间尺度的稳定度, 同时保证时间尺度的连续性。本文采用Kalman滤波器对观测钟差进行状态估计, Kalman滤波器每递推一次获得了一个新的预测值, 通过实时调整预测值来建立时间尺度。本文算法主要滤除了频率白噪声(white frequency modulation noise, WFM), 建立的时间尺度中主要只含有频率随机游走噪声(walk random frequency modulation noise, RWFM), 所以该时间尺度的中短期稳定度很高, 明显优于NIST的Kalman滤波器算法。
本文算法用于建立时间尺度的步骤和传统Kalman滤波器算法有很大的不同, 同时和Kalman滤波器用于状态估计和钟差预测[9-10]也不同。本文将首先描述时间尺度算法和Kalman滤波器的原理, 然后在此基础上论述本文算法的原理, 最后结合仿真实验进行分析。
-
原子钟的钟面读数定义为[5]:
(1) 式中, h(t)为原子钟的钟面读数; h0(t)=t为理想时间尺度; x(t)为原子钟的时差。
(2) 式中, TA(t)是建立的纸面时间尺度; hi(t)代表第i台原子钟的钟面读数; h'i(t)代表第i台原子钟的钟面读数的预测值; ωi为第i台原子钟的权重; N为钟组中原子钟的数量。
式(2)可以理解为传统的时间尺度算法中, 权重用于优化时间尺度的稳定度, 而预测值则用于消除每台钟确定性趋势项的影响。因此, 式(2)表明时间尺度TA是对去除确定性趋势项后只含有噪声分量的钟面读数的加权平均值。同时, 式(2)也暗示h'i(t)中确定性趋势项与hi(t)中确定性趋势项的符号相反。
在时间尺度算法中, 单台钟的时差xi(t)可以表示为[1-2]:
(3) 因此, 时间尺度基本方程也可以写为:
(4) 式中, xi(t)代表第i台原子钟的时差; x'i(t)代表第i台原子钟的时差的预测值。式(4)暗示x'i(t)中确定性趋势项与xi(t)中确定性趋势项的符号相同。
由式(2)和式(4)可知, 时间尺度算法的核心环节是通过N-1组观测量(钟差)计算得到N台钟的权重ωi和预测值h'i(t)或x'i(t)。
-
原子钟模型可以用随机微分方程[5, 11-12](stochastic differential equations, SDEs)来描述:
(5) 该SDEs的初值条件为:
(6) 求解该SDEs, 得到:
(7) 式中, X1、X2分别为原子钟的两个状态变量, X1代表时差, X2代表频差的随机游走部分, 即不含WFM的频差; W1(t)和W2(t)分别代表两个独立的维纳过程, 并且有W(t)~N(0, t), 即每个维纳过程服从均值为0、方差为时间t的正态分布, 该分布表明每个维纳过程的差分是白噪声, 即WFM的差分为相位白噪声(white phase modulation noise, WPM); σ1和σ2分别是这两个维纳过程的扩散系数, 用于表明噪声的强度; d为原子钟时差的频漂; s为(tk, tk+1)之间的时间。式(7)表明, 维纳过程W1(t)和维纳过程W2(t)对时间t的积分在状态变量X1上分别表现为WFM和RWFM。
将式(7)按照间隔T进行离散化, 得到:
(8) 式(8)就是原子钟的状态方程。状态变量中的噪声分量为:
(9) (10) 原子钟的观测方程[5]表述为:
(11) 式中, Z(tk)为观测量; σ·ξ(tk)为观测噪声, 其中ξ是标准高斯白噪声, 被认为是WPM, σ表示观测噪声的强度。观测噪声协方差表示为:
(12) 使用Kalman滤波器估计原子钟的状态变量, 其步骤可以通过下面5个方程[4]进行描述。
(13) (14) (15) (16) (17) 上述5个方程中各符号的含义是相关领域内所共知的, 在此不再描述各符号的含义。
Kalman滤波器作为最小均方意义下的最优估计算法, 可以有效地从观测量Z中估计状态X1和X2并得到它们的预测值。由状态方程(8)可知, X1的估计值和预测值中主要包含WFM和RWFM, 而X2的估计值和预测值中主要包含RWFM。另外, 按照文献[13]的方法进行验证, 该系统是一个完全可观测的系统。由系统的完全可观测性可以得知该Kalman滤波器最终会进入稳态, 即滤波器增益、估计误差矩阵、预测误差矩阵最终都是收敛的。于是, Kalman滤波器的初值设置不会影响Kalman滤波器进入稳态时的性能。但是, 为了使Kalman滤波器更快进入稳态, 本文的仿真实验中, 设置每台钟P矩阵的初值为P0=Q
, 矩阵的初值为[X1(1) (X1(2)-X1(1))/T]T。 -
本文算法可以分解为3个步骤。
1)重构钟差序列。对于一个含有N(N≥2)台钟的钟组, 可以获得N-1组钟差(不妨假设这些钟差为第n(n≠1)台钟与第1台钟之间的钟差)。分别使用N-1个相互独立的Kalman滤波器对这N-1组钟差进行估计, 得到N-1组估计值
1(1, n)(tk)和 2(1, n)(tk), 其中上标(1, n)表示第n台钟与第1台钟之间的钟差。对于每一组估计值, 令(18) 式中,
1(1, n)是重构的钟差序列。把
1(1, n)(tk)序列的初始值取为零, 式(18)可以改写为:(19) 显然
1(1, n)(tk)是通过估计量 2(1, n)(tk)进行重构得到的钟差序列。2)定义预测值。把重构得到的
1(1, n)(tk)定义为每台钟的时差与时差预测值之间的偏差。于是, 第1台钟的预测值表示为:(20) 其中, n可以是2, 3, …, N中任意值。
第n(2≤n≤N)台钟的预测值表示为:
(21) 3)建立时间尺度。将预测值的表达式(20)和(21)代入时间尺度基本方程(4), 可以得到时间尺度TA(tk):
(22) 式(22)最后一个等号右边的第一项中, 上标中的n可以是2, 3, …, N中任意值。
由式(22)可知, 时间尺度TA即为N-1组重构钟差序列
1(1, n)(tk)的加权平均。因此, 可以抛开时间尺度的标准定义式(4), 直接对这N-1个重构钟差序列进行加权平均, 最终可以得到如式(23)所示的时间尺度TA:(23) 对比分析式(22)和式(23), 这两个表达式是等价的。
需要说明一点的是, 对于权重的选取, 可以采用传统加权平均算法中权重的选取方法。大量的文献对此进行过研究, 所以在本文中不展开论述。
考虑一种特殊情况, 当钟组中只含有两台钟时, 将预测值的表达式(20)和(21)代入时间尺度基本方程(4), 得到时间尺度:
(24) 式(24)表明, 当钟组中只含有两台钟时, 权重对时间尺度TA没有任何影响, 时间尺度即为重构钟差序列。所以, 可以抛开时间尺度的定义(即时间尺度基本方程(4)), 直接把重构钟差序列作为时间尺度TA。
1)本文算法的核心思想是在Kalman滤波器每一次递推的过程中, 调整一次预测值, 通过每次实时调整预测值来建立时间尺度TA。
2)本文算法建立了一个实时、可预测的时间尺度。
3)钟组中只含有两台钟时, 可以抛开时间尺度的标准定义式(4), 直接把重构钟差序列作为时间尺度TA。当钟组中含有N(N≥3)台钟时, 也可以抛开时间尺度的标准定义式(4), 直接对这N-1组重构钟差序列进行加权平均, 得到时间尺度TA。
-
首先, 以钟组中只含有两台钟的情况为例进行说明。
式(24)可以改写为如下的形式:
(25) 其中, ε(tm)为Kalman滤波器的估计误差。
式(25)表明, 在不考虑误差项的情况下, 该时间尺度TA只含有RWFM, 所以中短期稳定度明显提高。
然后, 推广到钟组中含有N(N≥3)台钟的情况。分析式(23), 时间尺度TA(tk)是对N-1组主要含有RWFM的重构钟差进行加权平均, 所以同样主要含有RWFM, 中短期稳定度很高。随着钟组中钟数量的增加, 通过设计合理的权重, 该时间尺度的稳定度将得到进一步提高。
-
仿真生成三台铯钟组成钟组, 验证本文算法性能。
1)使用文献[14]的方法(文献[14]全文论述了生成原子钟5种不同类型噪声的方法, 本文使用该方法生成只含有WFM和RWFM两种噪声的噪声序列。)生成三台铯钟的时差, 每台铯钟共10 000个数据点。设置相邻数据点之间的时间间隔为T=14 400 s, 每台铯钟的噪声平方扩散系数都设置为σ12=1.50×10-23s, σ22=3.29×10-37s-1(该参数符合Symmetricom 5071A高性能铯束管铯钟的指标), 观测噪声方差为σ2=2×10-20s-2(即观测不确定度为0.14 ns左右, 目前的测量设备能达到该指标。)。为了简化分析, 设置确定性分量x0=y0=d=0。此时每台铯钟的状态变量真实值X1和X2都是已知量。
2)使用Kalman滤波器估计由这三台铯钟比对得到的两组钟差的状态。
3)运用本文提出的实时预测值调整算法建立时间尺度, 为了简化分析, 每组重构钟差的权重都设为0.5。
4)使用原始的NIST的Kalman滤波器建立时间尺度(记为T'A), 用于和本文建立的时间尺度TA进行相互比较。NIST的Kalman滤波器算法的基本原理和步骤见文献[4]。
-
使用文献[14]的方法生成一台铯钟(参数与上述描述相同), 使用Kalman滤波器对钟差进行状态估计, 得到
1和 2。图 1为X1和
1。从图 1中可以看出, 估计得到的 1和X1的曲线基本重合, 表明Kalman滤波器对于X1的估计误差很小。图 2为X2和X1的导数dX1/dt(即频差)。由式(7)可知, X1的导数dX1/dt和X2分别表示为:
图 2 X2和X1的导数dX1/dt(即频差)
Figure 2. X2 and dX1/dt Which is the Derivation of X1(the Frequency Deviation)
(26) (27) 对比式(26)和式(27), 状态变量X2可以认为是不含WFM的频差。从图 2可以看出, X2的噪声分量明显小于频差dX1/dt。
-
使用文献[14]的方法生成另外三台铯钟(参数与上述描述相同)组成钟组, 运用实时频差调整时间尺度算法, 建立时间尺度TA。同时使用NIST的Kalman滤波器算法, 建立时间尺度T'A。
图 4描述了两组观测钟差、TA和T'A的时差。图 5描述了其中一组观测钟差、TA和T'A的Allan偏差。分析图 4和图 5可知:观测钟差的噪声主要表现为WFM和RWFM, 而TA中不含WFM, 其主要噪声为RWFM; 相比于NIST的Kalman滤波器算法建立的T'A, 时间尺度TA的中短期稳定度很高。
图 5 其中一组观测钟差、TA和T'A的Allan偏差
Figure 5. Allan Deviations of One of the Observed Clock Differences, TA and T'A
综上, 本算法建立的时间尺度中短期稳定度很高。同时, 时间尺度TA的表达式(23)表明, 该TA是一个连续、实时、可预测的时间尺度。
-
本文提出了一种基于Kalman滤波器的实时预测值调整时间尺度算法。不同于传统算法主要通过调整权重来提高时间尺度的稳定度, 本文算法主要通过实时调整预测值来提高时间尺度的稳定度。本文算法的核心思想可以概括为用Kalman滤波器对观测钟差进行状态估计, 在Kalman滤波器每一次递推的过程中, 调整一次预测值, 通过每次实时调整预测值来建立时间尺度。理论分析和仿真实验表明, 本文算法滤除了WFM, 建立的时间尺度中主要只含有RWFM, 所以中短期稳定度很高。同时, 该时间尺度是一个连续、实时、可预测的时间尺度。
-
摘要: 时间尺度算法用于综合钟组内所有的原子钟,建立一个频率稳定度更高的时间尺度,其核心环节可以概括为通过N-1组观测量(钟差)计算得到N台钟的权重和预测值。传统算法主要关注如何调整权重来提高时间尺度的稳定度,本文算法通过调整预测值来提高时间尺度的稳定度。本文算法使用Kalman滤波器对观测钟差进行状态估计,在Kalman滤波器每一次递推的过程中,调整一次预测值,通过每次实时调整预测值来建立时间尺度。理论推导和仿真实验都表明,本文建立的时间尺度滤除了频率白噪声,主要只含有频率随机游走噪声,所以具有很高的中短期稳定度。该时间尺度是一个连续、实时、可预测的时间尺度。Abstract: The purpose of a time scale algorithm is to form a time scale with high frequency stability from an ensemble of clocks. The main component in a time scale algorithm generates weights and predictions for N clocks from the N-1 measured clock differences. Traditional algorithms mostly focus on how to set weights to improve the stability of the time scale. The algorithm we propose instead focuses on how to adjust predictions to improve stability. We use Kalman filters to estimate the states of measured clock differences. We adjust predictions with the every update of the Kalman filter, to form a time scale. Theoretical analyses and simulations both show that this algorithm filters out white frequency modulation noise, and the formed time scale mainly involves the random walk frequency modulation noise; hence, the time scale formed by our proposed algorithm has high short-term and middle-term stability. Itis also a continuous, real-time, and predictable time scale.
-
Key words:
- time scale algorithm /
- weight /
- prediction /
- frequency stability /
- Kalman filter
-
-
[1] 董绍武.守时中的若干重要技术问题研究[D].北京:中国科学院研究生院, 2007 http://cdmd.cnki.com.cn/article/cdmd-80024-2007178786.htm Dong Shaowu. Study on Several Important Technical Issues in Time-keeping[D]. Beijing: Graduate University of Chinese Academy of Sciences, 2007 http://cdmd.cnki.com.cn/article/cdmd-80024-2007178786.htm [2] Tavella P, Thomas C. Comparative Study of Time Scale Algorithms[J]. Metrologia, 1991, 28: 57-63 doi: 10.1088/0026-1394/28/2/001 [3] ParkerT E, Levine J.Impact of New High Stability Frequency Standards on the Performance of the NIST AT1 Time Scale[J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 1997, 44(6): 1 239-1 244 doi: 10.1109/58.656627 [4] JonesR H, Tryon P V. Estimating Time from Atomic Clocks[J]. Journal of Research of the National Bureau of the Standards, 1983, 88(1): 17-24 doi: 10.6028/jres.088.002 [5] Galleani L, Tavella P. Time and the Kalman Filter [J]. IEEE Control System Magazine, 2010: 44-65 http://cn.bing.com/academic/profile?id=2154685878&encoded=0&v=paper_preview&mkt=zh-cn [6] Greenhall C A. Forming Stable Timescales from the Jones-Tryon Kalman Filter[J]. Metrologia, 2003, 40: 335-341 doi: 10.1088/0026-1394/40/3/313 [7] Greenhall C A. A Kalman Filter Clock Ensemble Algorithm that Admits Measurement Noise[J]. Metrologia, 2006, 43: 311-321 doi: 10.1088/0026-1394/43/3/015 [8] Suess M, Greenhall C A. Congruence of Two Kalman Filter Composite Clocks[J]. Metrologia, 2013, 50: 499-508 doi: 10.1088/0026-1394/50/5/499 [9] 张清华, 隋立芬, 贾小林.应用Jones-Tryon Kalman滤波器对在轨GPS Rb钟进行状态监测[J].武汉大学学报·信息科学版, 2012, 37(4): 436-440 http://mall.cnki.net/magazine/article/whch201204013.htm 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://mall.cnki.net/magazine/article/whch201204013.htm [10] 黄观文, 张勤, 许国昌, 等, 基于频谱分析的IGS精密星历卫星钟差精度分析研究[J].武汉大学学报·信息科学版, 2008, 33(5): 496-499 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200805016.htm Huang Guanwen, Zhang Qin, Xu Guochang, et al. 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://www.cnki.com.cn/Article/CJFDTOTAL-WHCH200805016.htm [11] Zucca C, Tavella P. The Clock Model and Its Relationship with the Allan and Related Variances[J]. IEEE Transaction on Ultrasonics, Ferroelectrics, and Frequency Control, 2005, 52(2): 289-296 doi: 10.1109/TUFFC.2005.1406554 [12] Galleani L, Sacerdote L, Tavella P, et al. A Mathematical Model for the Atomic Clock Error [J]. Metrologia, 2003, 40: 257-264 doi: 10.1088/0026-1394/40/3/305 [13] 邓自立.卡尔曼滤波与维纳滤波:现代时间序列分析方法[M].哈尔滨:哈尔滨工业大学出版社, 2001 Deng Zili. Kalman Filter and Weiher Filter: Modern Time Sequence Analysis Method[M]. Harbin: Press of Harbin Institute of Technology, 2001 [14] Kasdin N J. Discrete Simulation of Colored Noise and Stochastic Processes and 1/f Power Law Noise Generation [J]. IEEE, 1995, 83(5): 802-827 doi: 10.1109/5.381848 -