基于Kalman滤波的原子时算法研究

孙同川，王振岭，孙建设，刘铁强

(1.中国电子科技集团 第54研究所，石家庄 050051; 2.中国人民解放军61711部队，新疆 喀什 844000)

0 引言

时间对一个国家的经济、社会生活、军事等方面有着十分重要的作用，时间的稳定关系着国家和社会的安全稳定，时间频率体系建设作为国家的重大基础设施，具有非常重要的意义。卫星导航、通信互联、联合作战、反导拦截、金融、电力等对时间频率系统的准确度和可靠性的需求越来越高。只有一个高精度的时频体系，才能提供可靠的时间频率服务。随着当前国产原子钟技术水平的突破[1]，为了提高自主可控能力，采用国产原子钟组进行原子时算法的研究具有重要意义。通常情况下，构建时间频率体系的关键是建立和保持系统时间，时间尺度的性能一方面受时间频率系统内硬件(比如原子钟)的性能水平的影响，一方面受所采用的时间尺度算法的影响[2]。

较为著名的时候尺度方法主要包括:目前国外普遍使用的ALGOS算法[3]、AT1算法[4]等，20世纪80年代，针对传统加权平均算法采用Allan方差计算权重会忽视主要噪声过程之外的噪声过程的问题，美国学者Barnes将Kalman滤波用于钟组噪声的处理提出了Kalman算法[5-6],近年来，国家授时中心对原子钟数据异常情况下如何计算时间尺度进行了研究，并采用改进的Kalman滤波计算时间尺度取得了比较好的结果[14]，国家授时中心还对Vondrak-Cepek滤波进行研究，并研究分析这种方法在时间尺度计算中的应用，这种方法计算出的综合时间尺度稳定度进一步提高[18]。ALGOS算法和AT1算法的核心思想都是加权平均，具体做法是根据某个计算周期内钟组内各个原子钟的稳定度来确定每台原子钟在此计算周期内的权重，通关比较合理的权重分配，使钟组内的噪声最小，以此来提高时间尺度的稳定度。但是加权平均算法只考虑起主要影响的噪声过程，而不论其他噪声过程如何，在计算时间尺度时并没有消除或抑制钟组内的噪声，相反Kalman算法通过Kalman滤波器对原子钟噪声建模消除或抑制钟组内的噪声。从文献[7]可以知悉Kalman、AT1算法的短期稳定度较好，而ALGOS的长期稳定度较好，AT1算法和Kalman算法是根据当前计算周期的钟差数据来预测下个计算周期的钟差，具有实时性的特点，而ALGOS算法则是采用过去一段时间的钟差数据来计算时间尺度[8]，具有滞后性，造成滞后的原因主要是因为ALOGS算法的权重计算和速率预报需要采用数据计算前一个月的钟差预测值。Kalman算法的稳定度较好，但是存在发散性的问题[8]。

守时实验室的目标是产生一个稳定、准确、可靠的时间尺度，在二级守时节点中，守时钟组类型单一(大多为铯钟组)，针对此前我国大多数实验室采用的原子钟大多为进口原子钟的问题，为了提高自主可控水平，经典的ALGOS算法主要考虑原子时的长期稳定度，在二级守时节点中综合原子时应同时注重综合原子时的短期稳定度和灵活性，以便守时系统随时进行调整。

因此对加权平均算法进行了研究和改进，改进后的算法将针对传统加权平均算法中的噪声问题，加入Kalman滤波过程用于改善频率预测值并加入了频率跳变检测，可以进一步抑制原子钟的噪声，提高综合原子时的稳定度。首先研究分析加权平均算法的基本原理，然后对算法进行改进，最后通过搭建的试验平台对改进的算法性能进行验证分析。

1 原子钟噪声模型的研究

时间尺度算法目的就是研究原子钟间的噪声问题，通过设定符合实际钟组情况的参数使得钟组内的噪声最小，以此提高时间尺度的稳定度。而原子钟作为一种表征频率的设备，其纸面读数与绝对时间之间必然存在一定的偏差，同时原子钟的纸面读数还受到观测噪声的影响，因此为了分析原子钟的噪声，首先要对其建立合适的噪声模型。

原子钟的观测方程可以表示为[20]：

Z(tk)=X1(tk)+σ·ξ(k)

(1)

式中，Z(tk)为观测量；σ·ξ(k)为观测噪声，ξ是标准高斯白噪声，属于WPM，σ表示噪声强度。

原子钟模型可用随机微分方程[8-9]描述：

(2)

式中,X1(t)，X2(t)表示原子钟的两个状态变量，X1表示时差状态变量，X2表示频率偏差的随机游走噪声，W1(t)，W2(t)是两个独立的维纳过程。为了方便分析仿真实验时将此微分方程的初值设置为0。J.A.Barnes等人提出原子钟噪声的幂律谱模型即原子钟噪声是五种噪声分量线性叠加的结果，文献[11]不仅描述了各个噪声特性还给出了仿真方法，它的谱密度函数：

(3)

不同的α取值对应不同类型的噪声，具体为[17]：α=-2是随机游走调频噪声(RWFM)，α=-1是闪变调频噪声(FFM)，α=0是调频白噪声(WFM)，α=1是闪变调相噪声(FPM)，α=2是调相白噪声(WPM)，不同原子时之间的算法受到各种噪声不同程度的干扰，但这些噪声有时不能避免，有时可以采用加入滤波器的方式进行消除。AT1算法中一般都采用了指数滤波器，指数滤波器可以比较有效的控制原子钟的调频白噪声和频率随机游走噪声。本文算法采用的滤波器是卡尔曼形式的，卡尔曼滤波器是最小均方意义下的最优估计，可以有效避免观测量中频差随机游走噪声的影响。

加权平均算法在计算过程中用Allan方差表征一台原子钟的稳定度，Allan方差主要与原子钟的噪声有关，每台钟的权重大小主要与某个起主要作用的噪声过程有关，针对加权平均算法中的权重和噪声过程问题，Kalman算法从估值的角度出发，避免了权重的分配问题，对所选取的主钟和标准时间的差做最小均方差意义下的最优估计。这种方法通过对每台原子钟建立一个状态方程评估的噪声模型，来构成Kalman滤波器。

假设原子钟的噪声模型符合：

(4)

(5)

式中，xi(t)、yi(t)为第t次测量时的时间和频率与理想时间尺度的偏差；Qi是系统噪声矩阵；ηi为随机游走调频噪声，其方差是Hi；ξi为白色调频噪声，其方差是Ei。

由Kalman滤波原理，得到钟组的动态模型为：

X(t)=φX(t-1)+W(t)

(6)

式中，φ为转移矩阵，W(t)为系统噪声矩阵。

钟组的观测方程为：

Z(t)=HX(t)+V(t)

(7)

(8)

式中,H是测量矩阵，V(t)为测量噪声

(9)

根据卡尔曼滤波模型和噪声矩阵，通过卡尔曼滤波迭代计算得到每台钟的状态，具体迭代步骤如下：

(10)

式中,P为X的误差方差矩阵；Q为系统驱动噪声矩阵；K为卡尔曼增益，R为测量噪声矩阵。

2 综合原子时算法研究与改进

目前守时实验室常用的综合原子时算法主要有两大类：加权平均算法和各种滤波类算法[15]。时间尺度产生最典型的例子就是世界协调时(UTC)的产生，国际计量局(BIPM)采用ALGOS算法对全球的守时实验室数据加权平均，以此提高时间尺度的可靠性和稳定性。上文提到加权平均算法的基本原理都是通过给每台原子钟分配一个权重来调整钟组内的噪声关系[16]，以此来提高时间尺度的稳定性。但ALGOS算法采用连续一个月的数据来预测前一个月的钟差值，具有滞后性参数调整设置极其不方便，需要的钟差数据时间太长， 不适合在规模较小的守时实验室使用。

理想状态下，时间尺度算法的最主要步骤是根据钟组内N台原子钟，利用N-1组钟差对钟组内的原子钟进行权重调整，以及钟差预报等。而时间尺度算法就是通过计算每台钟的权重来调整钟组内的噪声关系使其对钟组产生的时间尺度影响最小化，由此生成的时间尺度比钟组中任何一个单独的原子钟具有更高的可靠性、稳定性和频率精度。钟组产生的时间尺度是采用数学方法计算得到的，并不是某台钟的时间也不是他们简单的集合，这就是综合原子时，它的物理实现是通过对钟组中任意一个原子钟进行适当的校正来实现的，如果没有测量噪声的影响，该值与所采用的原子钟无关。时间尺度算法的输入是每个原子钟与选取的主钟之间的钟差，原子钟作为一个物理系统产生一个频率，原子钟的时间是由频率经过人工推导得来的，频率才是真正的物理量，因此用来衡量原子钟性能的参数都包含了对频率的描述，原子时算法中还需要估计每个原子钟频率偏移参数。

类加权平均算法在计算时间尺度的过程上大致相同，ALOGOS采用之前的钟差数据进行计算，AT1算法则是根据当前的钟差值估计下一个时间的钟差值。在钟组内选定一台钟作为主钟，主钟的物理信号输入相位微跃器，微跃器可以调整主钟的物理信号相位，微跃器和主钟的组合称为组合钟，每台原子钟和组合钟的钟差值都可以通过与主钟的钟差间接计算得出。物理主钟信号与相位微跃器调整后的主钟信号的钟差值可通过钟i与主钟的钟差间接计算得出，每台钟的钟差估值为xri(t+τ)，对钟差估值xri(t+τ)加权平均计算出主钟与组合钟的钟差xr(t+τ)，对相位微跃器则参照这个值对主钟物理信号进行调整。

假设钟i在时刻t的钟差和速率分别为xi(t)，yi(t)，测量时间间隔为τ。则钟组内某台钟t+τ时刻相对于组合钟的钟差可由下式计算：

xi(t+τ)=xi(t)+yi(t+τ)·τ

(11)

将计数器在时刻t+τ测量得到的钟i和主钟的钟差记为ti(t+τ)，则主钟与组合钟的钟差可由下式计算：

xri(t+τ)=xi(t+τ)-ti(t+τ)

(12)

即：

xri(t+τ)=xi(t)+yi(t)·τ-ti(t+τ)

(13)

这样就通过钟i的钟差和速率计算了组合钟相对主钟的钟差值。

假设钟组内有N台中，那么可以得到N-1个钟差值，在对其进行加权平均可以得到：

(14) )

式中,xri(t+τ)为第i台钟间接计算出的钟差估算值，ωi为它的权重，xi(t+τ)是钟组内所有原子钟计算的得出的钟差值的加权平均，xi(t+τ)是相位微跃器调整主钟输出信号的相位的参照，由此便生成了综合原子时的物理信号。为了削弱噪声对时间尺度稳定性的影响，AT1算法需要引入一个指数滤波器，指数滤波器的计算方式为：

(15)

式中，k是一个重要参数，这里取k=30，第i台钟的频率估计值为fi，通过时间差分计算得到：

(16)

由上述计算过程可以看出：AT1算法钟每个原子钟有两个状态，原子钟相对主钟的时间偏移和频率偏移，频率偏移可以用来计算，但算法本身不能对频率偏移进行估计。改进的算法结合传统的AT1算法和Kalman滤波来估计表示频率随机游走噪声和原子钟的频率偏移状态和方差。对于一个给定的原子钟，这种状态不是一种物理状态，而是在白噪声调制频率的情况下对频率偏移的数学估计，假设这种状态称之为“Y”，那么这个状态的方差就表示对偏移的数学估计的置信度。

将AT1算法中的钟i的钟差写为：

(17)

式中,Di是频率漂移率。

如果钟j在t+τ时刻的刻度钟差测量值是原子钟xij(t+τ)，则钟j在t+τ时刻的钟差预报偏差为：

xj(t+τ)=

(18)

原子钟的权重计算方式通过时间残差计算：

(19)

频率偏移量可通过下式计算：

(20)

把频率偏移量合并到钟i当前平均频率偏移的指数滤波估计中：

(21)

mi是钟i的指数频率加权常数，mi由钟i的调频白噪声和随机闪烁调频白噪声决定，用来衡量时间预测值的稳定性。

其中：

(22)

τmin是σy(τ)到达最小值的时间，τ0是用于计算σy(τ)的最小时间。

接下来需要对状态Y进行频率偏移估计，文献[12]和[13]对原子钟的频率稳定度和钟差预测不确定度进行了估计，采用一个简单的卡尔曼形式对频率偏移量进行过滤[21],假设Q(T)是测量噪声，对状态Y建模，状态Y包含随机游走噪声和固定漂移，那么Y不是一个由原子钟直接产生的物理层面的频率，而是过滤了白色调频噪声的之后的频率随机游走分量加上频率漂移。

将式(4)和(5)代入卡尔曼方程[4]，系统模型为：

yi(t+τ)=yi(t)+Diτ+η(τ)

(23)

测量模型：

(24)

式中,yi是频率直接测量值，εi是白噪声。

残差估计：

(25)

其中：

(26)

这相当于一个自适应卡尔曼滤波器，Y的残差可以从初始值进行变化，Y的指数滤波器参数也可随时间变化。

方程的稳态形式中：

(27)

此原子钟的指数频率加权常数为：

(28)

频率预测为：

(29)

由此可以看出卡尔曼形式也得到了类似AT1中的指数滤波器形式[19]，此方法继承了AT1算法对闪变频率的建模能力。

频率跳变的检测需要观测检查一段时间内的数据，考虑到实时性的操作，可以将某一区间内的平均频率偏移估计值与该区间开始时的滤波估计值进行比较。将频率偏差的估计值作为离群值进行测试，就可以确定频率跳变是否在最近发生。

通过测量得到的是时钟之间的相位差，由于受到频率闪烁噪声和随机游走噪声的影响，频率跳变会在某一时刻发生，有些情况下时钟本身小的频率波动会被当成随机噪声，为了区分这种现象，需要对频率跳变的大小进行限定。

本文的频率跳变检测方法是通过迭代每个时钟的测量范围来检测频率跳变，从当前测量时间到之后的某个时间，此时间间隔为τmin，Lmax是测量次数的最大值，Lmax是由τmin确定的整数。在每次测量过程中，对每台时钟的频率跳变进行检测，定义一个时钟在测量时间间隔内的平均频率yavg：

(30)

x-L是在t-L时刻的估计值，x-1是在t-1时刻的估计值。通过将这个平均频率和当前测量时刻的估计值进行比较，实际就是把固定时间区间上的实时频率估计值与该区间的平均频率估计值进行比较，如果两者频率差大于阈值，就可以认为时钟发生了跳变,具体判断方法为：

|yavg-y-L|>

(31)

通过这种方式检测频率跳变的发生，当检测到某台时钟发生频率跳变时，此时算法会将发生频率跳变的时钟从这个计算周期中剔除，也就是把发生频率跳变的时钟的权重变为0，此后的计算周期中对该时钟的频率继续观测，直到算法自适应学习了此时钟的新频率特性，然后将此时钟重新加入到原子时的计算当中。这种方法的好处在于原子钟发生频率跳变时，此台钟对综合原子时的影响会降低，有利于提高综合原子时的长期稳定度。

3 守时系统硬件设计及原理

二级守时节点可分为硬件和软件两个部分，守时系统的硬件部分主要包括：原子钟组、时频信号产生与分配单元、时差测量单元、溯源比对单元、数据存储分析处理和监控单元。硬件部分主要考虑物理信号产生问题、与上级节点或同级节点间的数据交互问题和可拓展性等；软件部分主要包括：溯源比对软件、综合原子时计算软件、数据采集分析软件和监控软件。稳定可靠的守时系统由硬件和软件共同支撑，缺一不可。守时系统模块图如图1所示。

图1 硬件系统结构图

本次试验的守时系统采用4台国产高性能铯原子钟组成的钟组，铯原子钟独立工作，不对原子钟的数据在频率和相位上进行改变。4台铯原子钟的原始信号直接输出到频率信号切换器和多通道计数器。4台铯钟输出的1PPS信号直接进入多通道计数器，多通道计数器的作用就是对钟差进行测量，测量结果通过UDP协议端口传输至数据库。综合原子时软件部分每隔一定的周期会从数据库中提取需要的钟差数据进行综合原子时的计算，并将计算结果和每台钟相对于综合原子时的时差存入数据库。

在参与守时的钟组4台铯原子钟里面选用一个稳定度较好的原子钟作为主钟也就是参考钟。4台铯原子钟产生的10 M信号输入频率切换器，频率切换器输出主钟物理信号，随后进入相位微跃器，综合原子时软件对综合原子时进行计算，同时相位微跃器依据综合原子时与主钟的钟差对主钟物理信号的相位进行调整，同时对主钟产生的10 MHz信号进行驾驭，完成主钟向综合原子时的溯源。

4 仿真实验

采用本单位4台铯钟的数据对本文算法进行了仿真验证，铯钟数据得采样周期为60 s，采样时间为120天数据点172 800个,计算周期为30天，软件算法使用C语言完成，运行环境是国产麒麟操作系统。钟差测量设备一般使用时间间隔计数器(SR620),SR620的观测噪声方差为σ2=1×10-20。分别采用AT1算法和本文改进算法对综合原子时进行了计算，并对两种方法的频率偏差值进行了分析；对频率跳变检测进行仿真；对本文算法计算得出的原子时稳定度验证分析。

首先对正常状态的铯钟频率预测值进行仿真分析，改进的方法对频率预测值有明显改善，频率偏差曲线更加平滑,波动更小如图2～3所示。

图2 AT1算法频率偏差

图3 改进算法的频率偏差

其中一台铯钟在第76天铯钟频率发生跳变，图4是铯钟的频率跳变曲线，图5是由本文检测方法计算的平均频率估计值，通过平均频率计算的跳变值与阈值进行比较，可以判断出铯钟频率发生了跳变，因此将发生跳变的铯钟从当前原子时计算周期中剔除也即权重置0。

图4 跳变铯钟的频率

图5 跳变铯钟的平均频率

最后根据4台铯钟组成的钟组系统，分别采用AT1算法和改进后的算法对时间尺度进行了计算，原子时稳定度曲线如图5所示。因为本文算法改善了频率预测值且在一定程度上抑制了频率跳变带来的影响，所以理论上来讲综合原子时的中长期稳定度比AT1算法更好，从图6的曲线可以看出理论与仿真基本符合。

图6 综合原子时稳定度

5 结束语

本文提出了一种Kalman调整频率预测值和检测频率跳变的方法，和一般加权平均算法不同的是，一般的加权平均算法更加强调权重分配，以此提高时间尺度的稳定度。本文则是直接改善频率的预测值，本文算法核心思想是在每次频率预报的过程中，采用Kalman滤波器对预测值进行改善，然后再计算综合原子时，可以看到Kalman滤波器可以很好地抑制WFM和RWFM，随着钟组的老化，频率发生变化，本文算法还可以对其检测并重新适应原子钟的新频率特性，所以算法的中长期稳定度更高，但时间尺度的短期稳定度提升不大，后续还需要进一步研究。除此之外，与Kalman算法直接计算出的综合原子时不同，本文算法计算的时间尺度是连续、实时、可预测的，可以满足规模较小的守时实验室需求。