时频卡尔曼滤波技术在原子钟的应用和前景展望

郭衍莹

[摘要]时频卡尔曼滤波是指这种滤波技术在校正、监测、控制原子钟时间频率方面的应用。它是在多学科（原子钟，卡尔曼滤波和精密时频计量技术）基础上发展的新技术。本文介绍它的发展历程，及其前景展望。并介绍时频卡尔曼滤波算法的基本概念，以及它和常规卡尔曼滤波的不同之处。随着原子钟小型化和芯片化的迅速发展，随着卡尔曼滤波算法芯片化，以及精密时频测量技术的推广，时频卡尔曼滤波将在卫星导航，武器制导，时统等方面获得越来越广泛的应用

[关键词]原子钟  卡尔曼滤波  时频卡尔曼滤波  铷原子钟  芯片原子钟  GPS导航  GPS制导

一、时频卡尔曼滤波技术的发展历程

近年来，对高频率稳定度的原子钟信号（无论信号来自GPS星种或地面守时系统，来自一般的原子钟甚至来自高稳定度石英钟)进行时间频率卡尔曼滤波处理，获得越来越广泛的应用。它既可对原子钟时差和频率进行监测（包括长期慢漂移），实时估计和预报；又可通过测试和对比进一步提高星载或地面原子钟的频率稳定度（包括短期，长期稳定度和频率漂移）。有些文献资料称这种滤波技术是对原子钟的“驯服”或“最优估计”等。但从本质上讲它是一种实时数据处理技术（从相位噪声中提取相位信息），所以本文遵循阿伦（D.W.Allan)论文的提法，称为“时间频率卡尔曼滤波”（Time & Frequency Kalman Filters.）或简称“时频卡尔曼滤波”。

众所周知，卡尔曼滤波是1960年匈牙利科学家卡尔曼通过他的博士论文提出一种新的滤波算法，通称卡尔曼滤波或卡尔曼算法。这是一种实时算法，它一方面利用系统的数学模型不断外推（预估）系统状态的估值，另方面不断测量（利用精密测量设备或仪器）状态的实时数据。最后通过对这两部分数据进行加权平均，得到卡尔曼滤波最优估值。 当系统模型噪声和测量噪声都服从高斯分布时（高斯白噪声）,在最小二乘意义下, 估计是线性最优的。由于它能对现场采集数据进行实时运算，更新和预测，所以从问世后就作为一种经典的数据处理技术, 广泛 应用于控制理论、信号分析和卫星导航星载原子钟的监测等领域。据说目前高校信息类专业，都将《卡尔曼滤波》作为必修课程。

将卡尔曼滤波应用于原子钟预测，则是上世纪80年代后之事。当时美国NBS的科学家阿伦等人于82年首先发表有关这方面论文。阿伦是时频学术界人士熟悉的科学家。他在60年代就和巴纳斯等科学家提出“振荡器相位噪声”的理论基础。他们指出振荡器的频率所以不稳定，主要是由于有相位噪声，而不是幅度噪声引起的。这是电子学一个重大的理论进展。所以卡尔曼滤波要处理的状态信息并非信号的幅度，而是它的相位信息。这就和传统的卡尔曼滤波的概念有所不同。阿伦还指出，振荡器的相位噪声基本上可用5类噪声来表征，它们的谱密度数学模型可用“幂律”来表征，即：白噪声调相（谱密度为h₂),闪烁噪声调相（h₁f^-1）,白噪声调频（h₀f^-2）,闪烁噪声调频（h_-1f^-3）,频率随机游走（h_-2f^-4）。1966年阿伦提出“阿伦方差”，用时域表证各种振荡源的相位和频率稳定度，直至今天它仍是相位噪声一种最广泛的计量方法。

1982年阿伦等人对GPS星载铯钟和铷钟的时差数据进行卡尔曼滤波（以地面标准时间为基准）试验，目的是进一步提高GPS星种时间预报的精度（GPS通过导航电文不间断地向地面发送时差数据；并定时更新）。他们在引用卡尔曼算法时遇到第一个难题是原子钟的相位噪声不是严格的高斯白噪声。不能直接利用卡尔曼运算公式。解决的办法是状态矩阵只包括频率（相位的微分，记作）和频率变化率（频率的微分，记作）二项。阿伦认为，前者的随机项就是白噪声调频，后者就是频率随机游走。它们的谱密度系数h₀和h_-2，以及它们方差值都是可以通过仪器测量和运算而得。然后就可以利用标准的卡尔曼公式求出最佳增益k,以及滤波后最佳估值。从阿伦发表的论文可见，GPS星钟的阿伦方差经滤波后，在测量时间段τ较小时有明显的改善；在τ较大时则改善很小。此项研究公布后（公布的主要是计算机仿真结果）引起世界各国科研人员的注意和兴趣；也遭到个别专家的质疑（认为他的一些近似公式缺乏理论根据）。1985年作者曾指导研究生在航天部铷原子钟专家熊承锡同志的支持和帮助下进行计算机仿真（并首次用BASIC语言编制卡尔曼滤波算法流程）。仿真结果请查阅文献{4][13]。可惜熊承锡研究员不幸英年早逝，此项研究未能继续进行下去。

从1985至1990年间，这一技术在世界范围内获得迅速发展。1988年在意大利召开的国际时频讨论会上，将近三分之一的论文与时频卡尔曼滤波有关。其中有关预测GPS星种数据的论文又占一大部分。美国NBS还在会上宣布，他的三个标准时间尺度[UTS(NIST),TA(NIST)和TA1]中的后二项都是经卡尔曼滤波处理的。并称二者频率稳定度都在1×10^-14量级（τ为一个月左右）。至于在理论上的进展有：一是状态矩阵的参数从二个增至三个。即相位（θ），频率（）和频率变化率（）。三个方程式中的噪声项可分别视为白噪声调相，白噪声调频和频率随机游走。三种噪声的方差值都可从仪器测量和运算中求得。这就大大提高卡尔曼滤波的效果。二是对星种一些频率慢漂移，国内外有些学者开始运用自适应卡尔曼滤波的办法来处理。所谓自适应，主要指在利用测量数据完成滤波的同时，还可对模型参数和噪声统计特性进行估计、修正和自动更新，起到修改滤波设计和减少滤波误差的作用。

这一阶段无论理论或是应用，都取得不小的成就。不过仍有一些主要问题没能得到解决或完全解决，有的还遗留至今。这些问题主要有：1）虽然状态矩阵已增至3项，但仍有二项相位噪声项，即闪烁噪声调相和闪烁噪声调频没能得到考虑和处理；2）卡尔曼滤波在处理原子钟频率慢漂移和频率异变时的效果较差。3）卡尔曼算法的芯片化尚没得到普及；成本、可靠性等方面亟需改进。

二、卡尔曼滤波和时频卡尔曼滤波的基本概念

1）卡尔曼滤波的概念

卡尔曼滤波有其严密的理论依据和繁复的数学推导。有兴趣的同学和读者可以参阅有关教科书和文献。本文只是试图从科普的角度，结合具体例子从概念上而不是从数学上向非本专业读者和同学介绍卡尔曼滤波的基本概念。当然如果读者是信息专业，完全可以忽略这一节，直接阅读下一节。

现以“反导系统是如何跟踪弹道导弹的”为例，介绍卡尔曼滤波算法的基本概念。滤波过程首先就是进行预估和测量。预估的第一步是建立状态方程式。假定地面某反导系统已在时刻k-1时得到目标（弹道导弹）的坐标参数（多参数时可写成状态矩阵，例如经度、纬度、高度）的最优估值为X_k-1，那么在时刻k时系统对目标坐标的估值应为

X_k=ΦX_k-1+W_k               (1)

式中Φ状态转移矩阵。它的物理概念是弹道导弹的轨迹基本上是抛物线或准抛物线，Φ是表征目标在抛物线上前后位置的几何关系。W_k代表目标偏离抛物线的随机误差；一般认为高斯白噪声。其协方差矩阵为Q。

第二步是建立测量方程式。反导系统在对目标的位置（坐标）不断进行预估的同时，系统中的相控阵雷达也在不间断地实测目标的坐标值。测量方程式为：

Y_K=HX_k+V_k                   （2）

式中Y为测量数据矩阵；H为观测矩阵，它代表预估值和观测值之间固定转换关系；V为测量噪声矩阵（协方差矩阵为R）。请读者注意，（2）式只是为了说清概念，不可直接使用（一些教科书上也是如此处理）。

现在我们手头就有二组重要数据：一组来自雷达测量数据Y_K，即式（2）；其误差为V。一组来自预估值X_k，即式（1）。其误差为W。将此预估值乘上因子H(如式（2）所示)，得到的HΦX_k-1，就可以直接和数据Y_K进行比较了。可以想象最优估值一定在二者之间。而且哪个误差小最优估值就向哪个靠拢。数学上讲就是使得最优估值误差的方差为最小。于是卡尔曼就定义了一个参数，称为卡尔曼增益k，代表最优估值向哪一组数据倾斜程度，使最后得到的最优估值的方差为最小。图1表示一个一维系统，它的Y_k（测量数据）和HΦX_k-1（预估值折算成测量数据）的相互关系，以及增益k的含义。弄清这些概念后接下来的工作就是代公式了。最重要的公式有二，一是增益k；二是最后得到的最优估值X_k。增益k的数学推导过程很复杂，这里从略；只将最后得到的最佳X_kOPT公式抄录如下：

X_kOPT=ΦX_k-1+k(Y_k-HΦX_k-1)               （3）

（3）式有明确的物理意义，可从图1来理解。

图1  一维系统的预估值、测量值、最佳估值和增益之间关系

2）时频卡尔曼滤波的基本概念

这一节试图从概念上介绍和分析时频卡尔曼滤波和一般传统卡尔曼滤波的不同之处，尽量避免繁复的数学推导。

时频卡尔曼滤波主要不同处，一是研究的状态变量是相位而不是幅度；二是原子钟是频率高稳定器件，状态方程式中的转移因子Φ可认为等于1；即认为原子钟频率保持不变；而将变化部分归结到“噪声”中。目前大多数学者研究的是三维系统，即相位、频率和频率变化率。今假定原子钟的相位变化为θ(t）；则频率变化y（t）=dθ(t）/dt；频率变化率a（t）=dy（t）/dt.于是我们定义：

θ(t）=x(t)+n₁（t）

dθ(t）/dt=y（t）+n₂(t)

dθ²(t）/dt²=a（t）+n₃(t)

式中n₁（t）、n₂（t）和n₃（t）分别为阿伦定义的相位白噪声，频率白噪声，和频率随机游走。它们的噪声方差分别为（可用仪器测得，τ为其采样周期）：

σ₁²=3h₂f_c/4π²τ²

σ₂²=h₀/2τ

σ₃²=2π²h_-2τ/3                         （4）

因此可得一组离散的状态方程：

x（k+1）=x(k)+τy（k)+τ²a(k)/2+n₁（k）

y(k+1)=y（k）+τa（k)+n₂(k)

a(k+1)=a(k)+n₃(k)

写成矩阵形式，有

卡尔曼滤波的过程就是一边进行预估，一边用精密的“相位噪声测试仪”或“频率稳定度测试仪”等仪器测量原子钟的数据。因此状态方程和测试方程仍可分别写成式（1）和（2）的形式。并进行卡尔曼滤波运算。式中Φ为状态转移矩阵:

当然以上所述仅限于三维系统。

当然以上所述仅限于三维系统。如果维数增加，运算就复杂得多，有的至今仍是没有解决的难题。

三、美国NBS学者早期所做工作简介

上世纪80年代初，美国NBS的阿伦等科学家对GPS星载铯钟和铷钟的时差数据进行卡尔曼滤波（以地面标准时间为基准）处理。他们假定状态参数只需要二维就可以了，所以从今天看来比较简单，但它具有启蒙性，它的设计思路有助于初学者建立清晰的概念。因此下面做一简要介绍。

阿伦等人为了简化运算，认为只考虑了频率白噪声，和频率随机游走二项误差已经足够。因此可得离散状态方程：

y(k+1)=y（k）+τa（k)+u₁(k)

a(k+1)=a(k)+u₂(k)

因此有：

u1和u2分别为频率白噪声和频率随机游走，他们的噪声方差公式已在前面给出。接下来的工作就是代公式，求得协方差Q、R以及增益K等参数。详细运算过程读者可参阅阿伦等人的原著。

不过在阿伦等人论文中公布的主要是他们的计算机仿真结果。并以一台地面钟代替星钟。该钟的σ₂=10ns;σ₃=3ns;τ=1天。然后经卡尔曼滤波处理，信号的频率稳定度得到改善。如图2，分三种情况：图a用共视法，由图可见阿伦方差大约改善半个数量级。图b是测单颗卫星的短期频率稳定度，τ取15s至900s，由图可见τ=15s时阿伦方差改善二个数量级。图c是对GPS-5号导航星钟（可能是铷钟）处理的实际数据，图示在11天时，阿伦方差改善11db。

图2  美国NBS科学家对原子钟的时频卡尔曼滤波结果

四、新世纪以来时频卡尔曼滤波在小型（芯片）原子钟控制中的应用及其发展前景

新世纪以来，一项引世人瞩目的成就就是中国科学院国家授时中心时间频率基准重点实验室的科学家们利用自适应卡尔曼滤波对国家授时的原子钟的频率异常进行预测和估计。这一成果对我国时间基准保持工作具有重要价值。成果的情况已发表在时频领域国际期刊Metrologia 上，它的发表标志着授时中心在守时理论与方法领域的研究工作获国际认可。

从科学院的公报来看，这一成果主要得益于采用自适应卡尔曼滤波。其实新世纪以来国内外已发表过多篇重要的自适应卡尔曼滤波的论文。不过可能限于专利和知识产权，详细细节都没有公布。

新世纪以来，利用卡尔曼滤波监测和改善原子钟稳定度的方法已在国民经济一些领域，如通信，电力控制中心等的原子钟得到应用，并引起人们兴趣和关注。以下二个方面既是原子钟当前研发的热点和未来若干年重要发展方向，也是时频卡尔曼滤波应用的热点和发展方向。值得我们注意：

1）小型铷钟的发展和卡尔曼滤波的应用

铷钟是三大经典原子钟之一（另二个为铯钟和氢钟）。它在三种钟中的优势是短期稳定度好；但长期稳定度差，频率漂移严重。在应用方面它的突出优点是造价最低，体积也最小，如果进一步小型化，再利用时频卡尔曼滤波提高其长期稳定度，则肯定会有广阔的发展前途。

在铷钟小型化方面，国内航天科工二院203所自主研发的拥有自主知识产权的超薄铷原子钟，本体尺寸只有76mm×76mm×17mm，是世界上厚度最薄的高可靠铷原子钟，可在强振动条件下保持可靠锁定，技术指标达到国际顶尖水平，目前已经完成研制并已实现量产。图3是中国航天科工集团官网于2019年8月23日公开发布的，关于203所研发超薄铷钟的信息和图片。

图3

小型铷钟在军事上有很好的发展前景。目前很多战术导弹、炸弹、炮弹等多采用GPS制导和定位（作者曾在文献[12]中指出，如GPS用户钟也改用原子钟，GPS制导的精度将获得明显提高）。据资料报道，导弹中GPS定位系统的用户钟过去大都用高精密，且能耐苛刻军用条件的石英钟（晶振）。近年来有的已改用小型超薄铷原子钟。因为无论从GPS定位精度，从体积重量，从环境条件以及成本来考虑，无论是导弹、炮弹还是炸弹，小型原子钟都远胜于石英钟，是未来战术武器弹载钟的重要发展方向。

小型铷钟如果再配上卡尔曼滤波，优势更明显。事实上已有学者在做这个工作。如文献[6]提出的卡尔曼滤波控制铷原子钟算法，可使铷钟获得了较好的准确性和长期稳定性,且对其短期稳定性影响最小,是一种可靠稳定的铷钟控制方法。该法已应用于实际，它以GPS接收机输出的1pps信号为参考信号,采用卡尔曼滤波算法对铷原子钟的参数进行预估,计算铷原子钟的频率调整量,对铷原子钟进行调整,使其和UTC时间保持同步.实验结果表明,受控铷原子钟输出1pps与UTC(NTSC)钟差的标准差优于3.5 ns,钟差峰值优于15 ns,100 s采样的阿伦方差为1.83×10-12,10000 s采样的阿伦方差为6.1×10-13.实验证明采用卡尔曼算法的思路是完全正确的。

图4  卡尔曼滤波控制某铷钟的原理框图

2）芯片原子钟（CPT）的发展和时频卡尔曼滤波的应用

CPT原子钟是一种新生事物，它是利用原子的相干布局囚禁（CPT）原理而实现的一种新型原子钟，由于不需要微波谐振腔，因此可以做到真正的微型化。CPT原子钟由于可以集成到一个芯片上，因此也被称为芯片尺度原子钟；国内则称之为芯片级原子钟。它是迄今为止能够用电池供电长时间工作的唯一的一种原子钟。美国商品化的最小CPT芯片原子钟为手表尺寸，约16毫米，用纽扣电池供电；稳定度达2×10^-10

一般商品化CPT钟，比起小型化铷钟，他们的稳定度指标接近，但前者体积和功耗都要小1至2倍，甚至1到2个数量级。在一些不能承受铷钟体积、功耗的应用场合，CPT钟无疑是个理想的选择。

目前国内宣称研发CPT钟的单位有航天系统203所，中国科学院武汉物理数学研究所，北京波时科技公司等几家。据武汉所声称，他们研制出体积200cm³，功耗4W，稳定度优于6×10^-11t^--1/2的小型CPT钟。波时公司研发的B-AC01型芯片级原子钟样机，体积为40×35×11mm³，已于今年8月经中国计量院正式计量校准，结果为：秒级稳定度（阿伦方差）为7.1×10^-11；日漂移率为5.4×10^-11。有人预言，CPT钟的发展和推广，会对GPS推广应用（尤其是军用）带来革命性的变化。

图5  国外芯片级原子钟产品

利用卡尔曼滤波对芯片原子钟进行控制，国内外均已有一些成功的例子，并已在文献上公开发表。例如文献[9][10]，虽是国内专利，但对一些细节报道得比较清楚，值得我们参考和借鉴。

（感谢）航天二院23所的秦中宇、蔡庆宇、刘兴术三位研究员（也是三位卡尔曼滤波方面的专家）对本文提出很多宝贵的、关键性的修改意见，谨此表示深切感谢!

参考资料：

[1]Jorgensen，P.S.,“Kalman Filter Estimation of the NAVSTAR

Satellite Clock Parameter, Proc.15th PTTI Application Planning Meeting”1983，P.531

[2]J.A.Barnes,D.W.Allan,Time Scale Stabilities Based on Time & Frequency Kalman Filters. Proc. of 39th Annual Frequency Control Symposium(1985)

[3]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

[4]郭衍瑩：信号处理技术在计量测试中的应用。宇航计测技术。4期连载；（1988年1期p.66;1988年3期p.64;1988年5期p.1;1989年3期p.69)

[5]黄新明等：基于Kalman滤波器的原子钟实时异常监测算法。2011年全国时频学术会议论文集，p.255

[6]樊多盛等：Kalman滤波的铷原子钟控制算法，宇航学报，2015年1期，p.90

[7]《基于自适应卡尔曼滤波的原子钟频率异常检测算法》

https://www.xianjichina.com/news/details_7022

[8]李孝辉等：用卡尔曼滤波器对原子钟进行控制。控制理论与应用，2003年4期p.551

[9]北大陈敬标等：基于卡尔曼滤波和/或DFB的激光芯片原子钟及实现方法。国知局专利（26049510，2021/7/27）

[10]孟红玲等：一种芯片原子钟驯服方法、装置和系统与流程。国知局专利（33324796，2023/3/3）

[11]Song hui-jie“An improved Kalman filter time scale algorithm for atomic clock noise variation”Acta Phys.2020 69(17)

[12]郭衍瑩：星载弹载原子钟的发展概况和展望，中国航天，2022年4期p.37

[13]郭衍莹：现代电子设备的频率稳定度。宇航出版社，1989年8月