《测绘学 》
构建与学术的桥梁 拉近与权威的距离
视频分享：国测一大队2020珠峰测量队员名单，复制以下链接打开抖音
国测一大队2020珠峰测量队员名单#测绘人 #珠峰
https://v.douyin.com/KNGsGq/ 复制此链接，打开【抖音短视频】，直接观看视频！
国际自由原子时原始权重算法初步分析
伍贻威
北京卫星导航中心, 北京 100094
收稿日期：2019-05-06；修回日期：2019-09-06
关键词：自由原子时(EAL) ALGOS算法 原子钟 加权平均 权重
Preliminary analyses of the original weighting algorithm of the echelle atomique libre
WU Yiwei
Beijing Satellite Navigation Center, Beijing 100094, China
First author: WU Yiwei(1987—), male, PhD, engineer, majors in key technologies of time references of satellite navigation systems.E-mail:Yiwei_Wu_sh@126.com.
Abstract: The echelle atomique libre (EAL) is the basis of international atomic time (TAI). It is a weighted average of several hundred atomic clocks of different laboratories spread worldwide. The weights of the original weighting algorithm of the weighted average algorithm (ALGOS) forming EAL is studied and discussed. The mathematical expectation values of the frequency deviation variances are derived. The mathematical distribution values of the frequency deviation variances are also derived with approximation. These expressions explain some experimental performances. Simulations validate the theoretical analyses.
Key words: echelle atomique libre (EAL) ALGOS algorithm atomic clock weighted average weight
协调世界时(coordinated universal time，UTC)是一个国际上的时间基准[1-5]，由国际计量局(Bureau International des Poids et Mesures，BIPM)计算并每月发布一次[2-3]。建立UTC，第一步是获取全世界几十个实验室的几百台原子钟的比对值，使用ALGOS算法建立自由原子时(echelle atomique libre，EAL)，目的是优化EAL的频率稳定度[6]。然后，使用基准频标对EAL进行驾驭，得到国际原子时(international atomic time，TAI)[6]。TAI的秒长与国际单位制(SI)中秒的定义一致。最后，在TAI的基础上引入闰秒，得到UTC[6]。EAL、TAI和UTC都是纸面时间。
ALGOS算法是一种加权平均时间尺度算法，经历了以下3个发展阶段。
(1) 在2011年前，权重反比于频差(也称为原子钟的速率)方差，所有氢钟和铯钟的预测值都按照一次多项式模型计算[6-7]。
(2) 从2011年开始，更改了预测算法，权重依然反比于频差方差，所有氢钟和铯钟预测值按照二次多项式模型计算[7-9]。
(3) 从2014年开始，更改了权重算法，权重按照频差预测误差的平方选取，预测算法依然采用2011年更改后的算法[9]。
目前国内外研究存在以下问题：①国内部分文献对ALGOS算法的原理描述还不够透彻；②没有推导给出权重的数学期望和数学分布的解析表达式，从而没有从理论上解释试验现象。
1 基本原理1.1 原子钟信号及其Allan方差
典型的氢钟和铯钟信号可以表示为[10-12]
(1)
式中，X1、X2分别为两个状态变量，X1代表时差，X2代表频差的一部分[10-12]；x0和y0分别为时差和频差的初值；d为频漂；W1(t)和W2(t)分别代表两个独立的维纳过程，并且有W(t)~N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布；σ1和σ2是扩散系数(diffusion coefficients)，用于表明噪声的强度。
W1(t)、W2(t)的积分在状态变量X1上分别表现为频率白噪声(white frequency modulation noise，WFM)和频率随机游走噪声(walk random frequency modulation noise，RWFM)[10]，即式(1)的第一个方程中，最后第二项代表WFM，最后一项代表RWFM。
对于WFM和RWFM这样的有色噪声[13-17]，一般用Allan方差来表征其频率稳定度。
Allan方差的定义为[18-21]
(2)
式中，E表示求数学期望；τ=tk+1–tk为平滑时间；x为时差，即式(1)中的X1；y为频差，定义为
(3)
由式(2)，Allan方差也可以写成
(4)
文献[11]经过详细推导得出(当y=d=0时)
(6)
将式(5)和式(6)代入式(4)，得到Allan方差的表达式为[10-11]
当y≠0时，式(7)依然成立。当d≠0时，同理推出Allan方差的表达式(上标添加波浪线，和无频漂时的Allan方差相区分)为[10-11]
1.2 时间尺度基本方程
任何加权平均时间尺度都可以化简，最终由时间尺度基本方程(basic time scale equations，BTSE)来表达。EAL通过BTSE表示为[6, 10]
(9)
式中，hi(t)代表第i台原子钟的钟面读数；h′i(t)代表第i台原子钟的钟面读数的预测值；ωi为第i台原子钟的权重，权重和等于1；N为钟组中原子钟的数量。
式(9)中，h′i(t)中确定性项和hi(t)中确定性项的符号相反。式(9)的物理含义是：EAL是对扣除确定性项后只含有噪声分量的各原子钟的钟面读数的加权平均值。由式(9)可见，算法的核心是确定权重ωi和预测值h′i(t)。
1.3 ALGOS算法基本步骤
由式(9)，得到[6, 10]
(10)
每台钟相对EAL的时差记为[6, 10]
把式(11)代入式(10)，得到
(12)
于是，得到如下方程组[6, 10]
(13)
该方程组共有N个方程。第1式为1个方程。第2式为N-1个方程，代表了其他N-1台钟相对于某一台钟的钟差。该方程组共有N个未知数，即每台钟相对于EAL的时差xj。权重和预测值都是已知量。求解上述方程组，结合式(11)，得到[6, 10]
(14)
综上，只要获知权重ωi和预测值h′i(t)，作为方程组的已知量，算法就可以自动计算出xj(t)，即计算得到了EAL。下节将描述计算这些已知量的方法。
1.4 ALGOS原始算法确定权重和预测值的方法
ALGOS原始算法中，每30 d或35 d为一个计算间隔。该间隔记为[tk,tk+1]，其中tk+1=tk+T，T=30 d或35 d。
2011年前，采用原始预测算法时，在本间隔[tk,tk+1]内，第i台钟的预测值表示为[6, 9]
式中，时差预测值表示为[6, 9]
频差预测值Bi(tk)或
(17)
式中，yi,k-1表示第i台钟在上一间隔[tk-1,tk]的频差。
2014年前，采用原始权重算法时，每个间隔的权重经过4次迭代计算方程组(13)得到[6, 9]。每次迭代的步骤如下[6, 9]:
(1) 在第1次迭代中，每台钟的权重等于上一间隔内该钟的权重；后续迭代中，权重等于本次迭代的结果。时差和频差的预测值分别根据式(16)和式(17)计算得到，并在本间隔内保持不变。将权重和预测值代入式(13)，计算得到每台钟的时差xj。
(2) 计算每台钟在本间隔的频差yi,k。
(3) 计算每台钟包括本间隔在内，最近12个频差的方差，记为σi2(12,T)。
(4) 每台钟的权重反比于其频差方差，即
(18)
(5) 若某台钟出现异常，则需要剔除出钟组，在后续迭代中权重为零。若某台钟达到权重上限2.5/N，则在后续迭代中权重值取上限值。除此之外，该钟在下次迭代中的权重等于式(18)的值。
最终，得到每台钟在该间隔[tk,tk+1]的权重ωi和频差yi,k(作为下一间隔的Bi(tk+1))。
1.5 新预测算法
2011年后，新预测算法采用二次多项式模型，即式(19)，预测所有原子钟[7-9]
(19)
式中，Ci(tk)为频漂预测值。
这时，Bi(tk)计算方法保持不变，采用每台钟相对于EAL的时差计算得到；Ci(tk)是采用每台钟相对于TT(BIPM)的频差计算得到。TT为地球时(terrestrial time)。TT(BIPM)是BIPM采用基准频标计算得到的TT，是独立于EAL的时间尺度，其秒长与SI一致。Ci(tk)的计算公式如下[7-8]
(20)
式中，yTT-hi为第i台钟相对于TT(BIPM)的频差，通过第i台钟相对于EAL的频差以及EAL相对于TT(BIPM)的频差换算得到；一个计算间隔的长度为T2，T2一般大于T，BIPM曾尝试过T2取3个月、4个月和6个月等。
改进的预测算法有效地消除所有原子钟频漂的影响，因此EAL的频漂幅度明显降低。
2 ALGOS原始权重算法的权重分析
本节将从理论上分析ALGOS原始权重算法的权重。为了聚焦研究内容，本节的推导是在理想情况下进行的：
(1) 假设原子钟模型严格符合式(1)，不存在时差跳变、频差跳变等异常情况。
(2) 假设每个间隔T都相等。
(3) 假设相对于每台钟，EAL为理想时间尺度，即式(11)中时差xi的确定性分量和噪声都为原子钟的确定性分量和噪声。
2.1 无频漂时算法权重分析
ALGOS原始权重算法中，12个频差组成了一个样本。其样本均值记为
(21)
频差方差σi2(12,T)本质上是一个样本方差，计算公式如下[22]
(22)
将式(21)代入式(22)，化简得到
(22)
设yi,k的起始时刻为tk=(k-1)T，由式(5)和式(6)，得到
(24)
(25)
把总体的方差记为σi2，以区别于样本方差σi2(12,T)。
按照《概率论与数理统计》，样本方差的数学期望等于总体的方差值[22]。把式(24)和式(25)代入式(23)，对比式(7)，得到
(26)
式(26)表明，频差方差的数学期望σi2略大于Allan方差的数学期望σi,y2(T)。当T足够小时，由于WFM起主导作用，对比式(26)和式(7)，这时σi2≈σi,y2(T)。
下面推导σi2(12,T)的数学分布。
令
(27)
由于
(28)
即在此近似条件下，σi2(12,T)服从χ2分布，其方差为2(σi2)2/11。
由上述推导过程，σi2(12,T)的数学期望是严格推导得到的，其数学分布和方差是在近似条件下推导得到的。
综上，得出结论：
(1) 无论T多大，σi2(12,T)的数学期望都等于σi2，即式(26)。
(2) σi2略大于σi,y2(T)。当T较小时，由于σi,y2(T)中WFM起主导作用，σi2≈σi,y2(T)。这时，按照频差方差σi2(12,T)取权时权重的数学期望约等于按照Allan方差σi,y2(T)取权时权重的数学期望。
(3) 当T越大时，式(27)越来越不近似于对角线元素值相等的对角矩阵，所以式(28)的近似不再精确。
(4) 由于σi2(12,T)存在一定的波动(由其数学分布和方差来表征)，每台钟的权重在不同间隔可能会略有改变。
2.2 存在频漂时算法权重分析
当存在频漂时，令
(29)
结合式(26)和式(8)，得出
(30)
根据式(30)，由于ALGOS算法中T=30 d，所以频差方差的数学期望
此时，频差协方差矩阵为
(31)
式中，
同理，当T和d都较小时，式(31)所示的协方差矩阵近似为对角矩阵，且对角线上元素的值近似相等。按照《概率论与数理统计》，这时有[22]
(32)
即在此近似条件下，
实际上，式(31)中，等号右边第1项的协方差表达式与式(27)完全相同，右边第2项的协方差表示为
(33)
将式(27)与式(33)代入式(31)，得到
(34)
对比观察式(27)、式(33)和式(34)，随着T和d的变大，式(34)中相应元素的值将会变得远大于式(27)中相应元素的值，所以相比式(27)，式(34)也就越来越不近似于对角矩阵，即12个频差越来越不近似于独立同分布。
这时，相比式(28)，式(32)的近似更不精确。
由上述推导过程，
综上，得出结论：
(1) 无论T和d多大，
(2) 由于ALGOS算法中T=30 d，所以
(3) 式(30)从理论上解释了ALGOS原始权重算法中氢钟权重很低的原因。可以推测：假如ALGOS原始权重算法中权重是按照反比于Allan方差来选取，氢钟的权重会有一定的提升。
(4) 随着T和d的变大，式(34)中协方差矩阵越来越不近似于对角线元素值相等的对角矩阵。因此，相比无频漂时的数学分布近似表达式(28)，式(32)的近似更不精确。
(5) 由于
表 1列出了采用ALGOS原始权重算法，无频漂和存在频漂时，作为取权依据的频差方差的数学期望、数学分布和方差的表达式。数学期望是严格推导得到的。数学分布和方差都是在近似条件下推导得到的，随着T的变大(铯钟)或T和d的变大(氢钟)，其近似都会越来越不精确。
表 1 频差方差的数学期望、数学分布和方差Tab. 1 Mathematical expectation, mathematical distribution and variances of frequency deviation variances
| 有无频漂 | 数学期望 |