-
摘要: 针对多种分布形式混合的观测数据，建立了p范混合模型，考虑到模型中混合数属于不完全数据，引入期望最大化(expectation-maximum, EM)算法，对该混合模型的参数进行估计，详细推导了p范混合模型参数估计的迭代公式，并给出了相应的迭代步骤。采用混合高斯分布数据、拉普拉斯分布与高斯分布混合数据及实测GPS观测值残差数据，验证了公式的正确性和适应性。算例结果表明，与单一概率分布相比，p范混合模型能够准确反映数据分布的实际情况，同时利用EM算法估计的模型参数具有较高的精度。Abstract:Objectives Aiming at the mixed observation data of multiple distribution forms, a expectation-maximum (EM) combined p-norm distributed model(EM_p) is established.Methods Considering that the mixed number in the mixture model belongs to incomplete data, the EM algorithm is introduced to estimate the parameters of the mixture model and the p-model mixture model parameters are derived in detail. The estimated iteration formula and the corresponding iteration steps are given.The mixture Gaussian distribution data, Laplace distribution and Gaussian distribution mixture data, and the residual data of measured global positioning system(GPS) observations are used to verify the correctness and adaptability of the formula in this paper.Results and Conclusions The results of the calculation examples show that, compared with the single probability distribution, the p-norm mixture model can accurately reflect the actual situation of the data distribution, and the model parameters estimated by the EM algorithm have higher accuracy.
-
守时实验室需要建立和维持一个准确、稳定、可靠的时间尺度作为时间基准[1-3]，其核心算法包括时间尺度算法[4]、钟差预测算法[5-6]、驾驭算法[7-9]等。原子钟是守时系统的核心部件。时间基准的性能不仅和算法设计有关，还和参与计算的各台原子钟的性能有关。每种算法需要针对各参与计算的原子钟的模型和频率稳定度来设计。因此分析原子钟的模型和稳定度具有重要意义。
目前的研究中，文献[10-13]只分析了指定平滑时间的稳定度，或只对时差观测量中的某一具体分量展开分析，没有具体分析时差观测量中各分量对Allan偏差的贡献，并通过估计值外推得到所有平滑时间的Allan偏差估计值；文献[12-14]对Kalman滤波器估计原子钟状态原理描述不清晰；当观测噪声过大、存在周期性波动时，无法使用斜率法准确直接估计原子钟噪声强度[15]；当周期性波动在时差中不明显时，目前估计方法较为复杂[12, 16-17]。
本文针对上述问题展开研究，探索原子钟模型和频率稳定度分析方法。大量实验表明，典型氢钟和铯钟的观测模型可以表示为[2, 18-20]：①确定性部分，用二次多项式（时差、频差和线性频漂）加周期性波动项表示；②随机性部分，即原子钟噪声，主要为频率白噪声（white frequency modulation noise，WFM）和频率随机游走噪声（walk random frequency modulation noise，RWFM）；③观测噪声，噪声类型为相位白噪声（white phase modulation noise，WPM）。本文详细分析了该模型各分量的Allan偏差表达式。
在此基础上，本文从最优估计和低通滤波器[7-9]两个角度描述Kalman滤波器（Kalman filter，KF）估计原子钟状态的原理；提出了综合KF状态估计的结果和Allan偏差图估计原子钟噪声和观测噪声强度的方法；提出了3种不同的估计线性频漂幅度的方法；结合原子钟随机微分方程模型，提出了综合Kalman滤波器状态估计的结果和对数Allan偏差图估计周期性波动周期和幅度的方法。
1 原子钟模型及频率稳定度
1.1 原子钟时差和时差量模型
典型的氢钟和铷钟，其确定性分量用二次多项式来表示；典型的铯钟，其确定性分量用一次多项式来表示。二次多项式模型[2]表示为：
$$ x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + {\varepsilon _x}\left( t \right) $$ (1) 式中，等式右边前3项对应确定度分量，其中x0为初始时差，y0为初始频差，d为线性频漂；εx(t)对应时差的随机性部分，即噪声。一次多项式模型即式（1）中线性频漂d=0的情况。
某些时候，时差的确定性分量中还包含周期性波动分量。为简化分析，假设瞬时频差的周期性波动为标准正弦波或余弦波形式[19]，即：
$$ {y_s}\left( t \right) = A\cos \left( {2{\rm{ \mathit{ π} }}{f_0}t + \varphi } \right) $$ (2) 式中，A为幅度；f0为周期性波动的频率；φ为初始相位。
这时，频差周期性波动在时差上表现为：
$$ \begin{array}{*{20}{c}} {{x_s}\left( t \right) = \int_0^t {{y_s}\left( s \right){\rm{d}}s} = A\int_0^t {\cos \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right){\rm{d}}s} = }\\ {\frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right.} \end{array} $$ (3) 于是，二次多项式叠加周期性波动项的时差模型表示为：
$$ \begin{array}{*{20}{c}} {x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + }\\ {\frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + {\varepsilon _x}\left( t \right)} \end{array} $$ (4) 通过大量观测发现，氢钟、铯钟时差中起主导作用的噪声是RWFM[2, 18, 20-21]和WFM。这两种噪声在时差上可以分别用维纳过程和积分维纳过程来建模[2, 18, 20-21]，即：
$$ {\varepsilon _x}\left( t \right) = {\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} $$ (5) 式中，εx(t)为时差的随机性部分，即原子钟时差的噪声；W1(t)和W2(t)分别代表两个独立的维纳过程，并且有W(t)~N(0, t)，即每个维纳过程服从均值为0、方差为时间t的正态分布；σ1和σ2分别是这两个维纳过程的扩散系数，用于表明噪声的强度。
把式（5）代入式（4），时差可表示为：
$$ \begin{array}{*{20}{c}} {x\left( t \right) = {x_0} + {y_0}t + \frac{1}{2}d{t^2} + \frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + }\\ {{\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} } \end{array} $$ (6) 时差的观测量表示为：
$$ \begin{array}{*{20}{c}} {z\left( t \right) = x\left( t \right) + \sigma \varepsilon \left( t \right) = {x_0} + {y_0}t + }\\ {\frac{1}{2}d{t^2} + \frac{A}{{2{\rm{ \mathit{ π} }}{f_0}}}\sin \left( {2{\rm{ \mathit{ π} }}{f_0}s + \varphi } \right)\left| {_0^t} \right. + }\\ {{\sigma _1}{W_1}\left( t \right) + {\sigma _2}\int_0^t {{W_2}\left( s \right){\rm{d}}s} + \sigma \varepsilon \left( t \right)} \end{array} $$ (7) 式中，ε(t)为WPM；σ用于表明观测噪声的强度。
1.2 各分量Allan方差的表达式
时域上通常用Allan方差来表征频率稳定度。Allan方差的定义如下[2, 19]：
$$ \sigma _y^2\left( \tau \right) = \frac{1}{2}E\left[ {{{\left( {\bar y\left( {{t_k} + \tau } \right) - \bar y\left( {{t_k}} \right)} \right)}^2}} \right] $$ (8) 式中，τ为平滑时间；y为平均频差，定义为：
$$ \bar y\left( {{t_k}} \right) = \frac{1}{\tau }\int_{{t_k} - \tau }^{{t_k}} {y\left( t \right){\rm{d}}t} = \frac{{x\left( {{t_k}} \right) - x\left( {{t_k} - \tau } \right)}}{\tau } $$ (9) 式中，y为瞬时频差；x即为式（6）中定义的时差。
实际上常用Allan方差的平方根，即Allan偏差来表征稳定度。文献[22]详细推导了扩散系数与Allan方差的关系，即：
$$ \sigma _y^2\left( \tau \right) = \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau $$ (10) 式中，σy2(τ)代表平滑时间为τ时的Allan方差。式（10）中等式右边第一项斜率为-1，第二项斜率为1，分别对应WFM和RWFM。这说明在对数Allan方差图中，WFM的斜率为-1，RWFM的斜率为1。在实际应用中，很容易通过对Allan方差拟合得到扩散系数的值。
文献[20]推导了WPM与Allan方差的关系：
$$ \sigma _y^2\left( \tau \right) = 3{\sigma ^2}/{\tau ^2} $$ (11) 式（11）表明，WPM在对数Allan方差图中斜率为-2。
把式（1）和式（9）代入式（8），得到线性频漂d与Allan方差的关系：
$$ \sigma _y^2\left( \tau \right) = \frac{1}{2}{d^2}{\tau ^2} $$ (12) 式（12）表明，d在对数Allan方差图中斜率为2。式（2）所示的周期性波动项和Allan方差的关系为[19]：
$$ \sigma _y^2\left( \tau \right) = {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}} $$ (13) 综上，式（6）所示的时差的Allan方差表示为：
$$ \sigma _y^2\left( \tau \right) = \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau + \frac{1}{2}{d^2}{\tau ^2} + {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}} $$ (14) 式（7）所示的时差观测量的Allan方差表示为：
$$ \begin{array}{*{20}{c}} {\sigma _y^2\left( \tau \right) = 3{\sigma ^2}/{\tau ^2} + \sigma _1^2/\tau + \frac{1}{3}\sigma _2^2\tau + }\\ {\frac{1}{2}{d^2}{\tau ^2} + {A^2}\frac{{{{\sin }^4}\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}}{{{{\left( {{\rm{ \mathit{ π} }}{f_0}\tau } \right)}^2}}}} \end{array} $$ (15) 式中，等式右边第1项为观测噪声；第2、3项为原子钟噪声；第4项为频漂；第5项为周期性波动在Allan方差中的分量。
由式（15）看出，WPM、WFM、RWFM和线性频漂在对数Allan方差图中的斜率分别为-2、-1、1和2，在对数Allan偏差图中的斜率分别为-1、-1/2、1/2和1。假如时差中存在较大幅度的周期性波动，那么Allan方差也存在明显的周期性波动，即Allan方差在中间某一平滑时间段内会凸起来。
但是，并不是每台钟的模型中都含有上述所有分量。例如大部分铯钟几乎没有线性频漂，大部分型号的地面钟都没有周期性波动。
2 时差模型及频率稳定度分析方法
本节首先从最优估计和低通滤波器两个角度描述了KF估计原子钟状态的原理；然后展示分析WPM、WFM、RWFM强度、线性频漂幅度、周期性波动的幅度和频率的方法。
2.1 使用Kalman滤波器估计原子钟状态
式（1）所示的原子钟模型用随机微分方程（stochastic differential equations，SDEs）闭合解的离散形式表示为[2, 18, 20-21]：
$$ \left\{ \begin{array}{l} x\left( {{t_{k + 1}}} \right) = x\left( {{t_k}} \right) + {x_2}\left( {{t_k}} \right)T + \frac{1}{2}d{T^2} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\sigma _1}\left( {{W_1}\left( {{t_{k + 1}}} \right) - {W_1}\left( {{t_k}} \right)} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;{\sigma _2}\int_{{t_k}}^{{t_{k + 1}}} {\left( {{W_2}\left( s \right) - {W_2}\left( {{t_k}} \right)} \right){\rm{d}}s} \\ {x_2}\left( {{t_{k + 1}}} \right) = {x_2}\left( {{t_k}} \right) + dT + {\sigma _2}\left( {{W_2}\left( {{t_{k + 1}}} \right) - {W_2}\left( {{t_k}} \right)} \right) \end{array} \right. $$ (16) 式中，x和x2分别代表了两个状态分量，x和式（1）完全相同，代表时差，x2代表频差的一部分，即不含WFM的频差；T为时间间隔，T=tk+1-tk；其他符号的含义和§1.1相同。其中，噪声分量为：
$$ {J_k} = \left[ {\begin{array}{*{20}{c}} {{\sigma _1}\left( {{W_1}\left( {{t_{k + 1}}} \right) - {W_1}\left( {{t_k}} \right)} \right) + {\sigma _2}\int_{{t_k}}^{{t_{k + 1}}} {\left( {{W_2}\left( s \right) - {W_2}\left( {{t_k}} \right)} \right){\rm{d}}s} }\\ {{\sigma _2}\left( {{W_2}\left( {{t_{k + 1}}} \right) - {W_2}\left( {{t_k}} \right)} \right)} \end{array}} \right] $$ (17) 使用KF对上述两个状态分量进行估计，得到估计值记为$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$，其步骤用以下5个方程表示：
$$ {\mathit{\boldsymbol{\hat s}}_{k, k - 1}} = \phi {\mathit{\boldsymbol{\hat s}}_{k - 1, k - 1}} $$ (18) $$ {P_{k, k - 1}} = \phi {P_{k - 1, k - 1}}{\phi ^{\rm{T}}} + \mathit{\boldsymbol{Q}} $$ (19) $$ {\mathit{\boldsymbol{K}}_k} = {P_{k, k - 1}}{H^{\rm{T}}}{\left( {H{P_{k, k - 1}}{H^{\rm{T}}} + R} \right)^{ - 1}} $$ (20) $$ {\hat s_{k, k}} = {\hat s_{k, k - 1}} + {K_k}\left( {{z_k} - H{{\hat s}_{k, k - 1}}} \right) $$ (21) $$ {\mathit{\boldsymbol{P}}_{k, k}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_k}\mathit{\boldsymbol{H}}} \right){\mathit{\boldsymbol{P}}_{k, k - 1}} $$ (22) 式中各符号意义见文献[12, 21]。其中，H=[1 0]；$\phi $=$\left[ {\begin{array}{*{20}{c}} 1&T\\ 0&1 \end{array}} \right]$；过程噪声方差矩阵Q和观测噪声方差矩阵R分别表示为：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Q}} = {\rm{E}}\left[ {\left( {{J_k} - 0} \right){{\left( {{J_k} - 0} \right)}^{\rm{T}}}} \right] = }\\ {\left[ {\begin{array}{*{20}{c}} {\sigma _1^2T + \frac{1}{3}\sigma _2^2{T^3}}&{\frac{1}{2}\sigma _2^2{T^2}}\\ {\frac{1}{3}\sigma _2^2{T^3}}&{\sigma _2^2T} \end{array}} \right]} \end{array} $$ (23) $$ \mathit{\boldsymbol{R}} = {\sigma ^2} $$ (24) 当原子钟时差符合模型（16），Q和R的值符合${\sigma ^2}$、${\sigma _1}^2$和${\sigma _2}^2$时，KF估计得到的$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$是最小均方意义下的最优估计。从频域角度，确定Q和R的值相当于确定KF的带宽[8-9]；KF作为低通滤波器，其作用就是滤除时差观测量中高频的WPM，$\hat x\left({{t_k}} \right)$中保留了低频的WFM和RWFM，${\hat x_2}\left({{t_k}} \right)$中只保留了最低频的RWFM[8-9]。
文献[12-14]将x2理解为频差，但笔者认为，频差是由时差的差分得到的，其中依然包含了WFM和RWFM，而x2中只包含RWFM。
由于$\hat x\left({{t_k}} \right)$中包含WFM和RWFM，而${\hat x_2}\left({{t_k}} \right)$中只包含RWFM，噪声分量更少，所以，从${\hat x_2}\left({{t_k}} \right)$中可以更有效地估计频漂、周期性分量等。当含有周期性波动项时，尽管原子钟模型不完全符合式（16），但从KF作为低通滤波器的角度[8-9]出发，KF依然可以有效滤除观测噪声，估计出$\hat x\left({{t_k}} \right)$和${\hat x_2}\left({{t_k}} \right)$。
2.2 使用斜率法和Kalman滤波器分析WPM、WFM和RWFM的强度
WPM、WFM和RWFM在对数Allan偏差图中的斜率分别为-1、-1/2、1/2，理论上可以通过斜率拟合出σ2、σ12和σ22的估计值，记为${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$。但是，实际上很多情况，尤其是远程比对的情况下，观测噪声很大，造成对数Allan偏差图中WFM分量淹没在WPM分量中。此外，周期性波动幅度很大时，Allan偏差图中间部分会向上凸起。这些都可能造成无法观察到WFM分量。针对这种情况，这里提出一种在${\hat \sigma _1}^2$无法准确获知的情况下，使用KF结合时差观测量z(t)和时差估计值$\hat x\left(t \right)$的Allan偏差，获取${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$的方法。步骤如下：
1) 计算式（7）所示z(t)的Allan偏差，斜率拟合${\hat \sigma _2}^2$和${\hat \sigma ^2}$；
2) 设置一个σ12的大致估计值，代入式（23）确定Q的值，使用KF估计得到Q的估计值，并画出$\hat x\left(t \right)$的Allan偏差图；理论上讲，$\hat x\left(t \right)$的噪声分量应该只包含WFM和RWFM；
3) 根据$\hat x\left(t \right)$的Allan偏差图，反复调整σ12的估计值，重复步骤2），观察$\hat x\left(t \right)$的Allan偏差图，确保噪声分量斜率为-1/2和1/2，即只含有WFM和RWFM，最终获取${\hat \sigma _1}^2$。
2.3 分析线性频漂的幅度
多种方法可以估计线性频漂$\hat d$。如直接对时差观测量z(t)进行最小二乘拟合（方法1），或对Allan偏差的斜率进行拟合（方法2），或对KF估计得到的${\hat x_2}\left(t \right)$的斜率进行最小二乘拟合（方法3）。
参照文献[5-6]，可以从理论上分析这些方法的估计不确定度，本文不展开分析。按照§2.1的分析，由于KF得到的${\hat x_2}\left(t \right)$中滤除了WFM，只含有RWFM，此外典型氢钟或铯钟WFM强度远大于RWFM强度，即σ12远大于σ22，所以从直观上理解，方法3的估计不确定度最小。
2.4 分析周期性波动分量
1）分析f0。按照§2.1方法，尽管式（16）没有对周期性波动建模，KF依然可以估计得到$\hat x\left(t \right)$和${\hat x_2}\left(t \right)$，由于${\hat x_2}\left(t \right)$中只含有RWFM，周期性波动在z(t)和$\hat x\left(t \right)$中不明显，在${\hat x_2}\left(t \right)$中很明显。可以通过观测较长时间段内${\hat x_2}\left(t \right)$波峰波谷的位置，得到波动频率的估计值，记为${\hat f_0}$。
2）分析A。由于噪声的存在，观察${\hat x_2}\left(t \right)$无法准确获取A的估计值，记为$\hat A$。本文采用如下方法：在$\hat x\left(t \right)$的Allan偏差图中，反复调整$\hat A$的值，观察式（13）平方根所示的周期性波动的Allan偏差与$\hat x\left(t \right)$的Allan偏差的吻合程度，最终得到$\hat A$。
3 实验分析
本文采用两台国产氢钟（分别记为Hm1和Hm2）相对于参考时间基准（记为Ref）的实测钟差（Hm1-Ref，Hm2-Ref）进行分析，以验证提出的原子钟模型和频率稳定度分析方法。其中Ref和国际协调世界时同步，稳定度远高于单台国产氢钟。
3.1 使用斜率法和KF分析WPM、WFM和RWFM的强度和线性频漂的幅度
取某一段长度约75 d的Hm1-Ref实测数据，作为时差观测量z(t)。按照§1.2的方法估计${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$，得到${\hat \sigma ^2}$=1×10-22 s2，${\hat \sigma _1}^2$=3×10-26 s和${\hat \sigma _2}^2$=1.2×10-33 s-1。其中，${\hat \sigma ^2}$=1×10-22 s2意味着观测噪声的标准差为0.01 ns或10 ps，符合测量设备的指标。
图 1展示了KF前后z(t)和$\hat x\left(t \right)$的Allan偏差，同时也画出了根据上述估计值，按照式（10）、（11）的平方根计算得到的WPM、WFM、RWFM部分的Allan偏差。综上，本实验表明§2的方法可以在观测噪声和周期性波动幅度很大时，有效估计出${\hat \sigma ^2}$、${\hat \sigma _1}^2$和${\hat \sigma _2}^2$。
![]()
图 1 z(t)和$\hat x\left(t \right)$的Allan偏差Figure 1. Allan Deviations of z(t) and $\hat x\left(t \right)$3.2 采用多种方法分析和验证线性频漂的幅度
本节分别对§2.3的3种方法进行验证。
方法1：对$\hat x\left(t \right)$分别用一阶多项式和二阶多项式进行最小二乘拟合，得到线性频漂的估计值为$\hat d$=-3.891×10-20 s/s2。图 2(a)和2(b)分别画出了拟合残差。对比图 2(a)和2(b)，可以明显看出Hm1中存在线性频漂。
![]()
图 2 $\hat x\left(t \right)$的一阶残差和二阶残差Figure 2. First Order and Second Order Residuals of $\hat x\left(t \right)$方法2：验证方法1。将方法1得到的$\hat d$值代入式（12）再平方根，计算得到Allan偏差的频漂部分。图 1中画出了Allan偏差的频漂部分，可知两种方法得到的$\hat d$值基本吻合。
方法3：图 3画出了KF估计得到的${\hat x_2}\left(t \right)$，可知${\hat x_2}\left(t \right)$中存在明显的线性频漂。对${\hat x_2}\left(t \right)$的斜率进行最小二乘拟合，得到$\hat d$=-3.86×10-20 s/s2，和方法1的结果基本吻合。
综上，§2.3的3种方法估计结果基本吻合。实际上，直观上分析，方法3估计不确定度较小，原因在于KF估计得到的${\hat x_2}\left(t \right)$中滤除了WFM，只含有RWFM。假如直接对时差一次差分或者二次差分，得到的序列中同样含有WFM和RWFM差分得到的噪声，由于WFM强度远大于RWFM，将导致无法直接在这些差分序列中观测到图 3中所示的周期性波动，周期性波动已经淹没于噪声中。
3.3 分析周期性波动项
采用§2.4的方法，观察图 3波峰波谷的位置，得到${\hat f_0}$=1/86 400 Hz，$\hat A$大约在1×10-14~2×10-14 s/s左右。反复比较式（13）开平方后所示的周期性波动和$\hat x\left(t \right)$的Allan偏差，最终得到$\hat A$=1.6×10-14 s/s。图 4画出了周期性波动项和$\hat x\left(t \right)$的Allan偏差，可见两者比较吻合。本实验表明了§2的方法可以有效估计周期性波动的周期和幅度。
![]()
图 4 周期性波动项和$\hat x\left(t \right)$的Allan偏差Figure 4. Allan Deviations of the Sinusoidal Component and $\hat x\left(t \right)$3.4 实例2分析
采用和§2相同的方法分析第2台国产氢钟Hm2的性能，得到该氢钟的参数估计值分别为：${\hat \sigma ^2}$=1×10-22 s2，${\hat \sigma _1}^2$=3×10-26 s，${\hat \sigma _2}^2$=8×10-33 s-1，$\hat d$=-3.8×10-20 s/s2，${\hat f_0}$=1/86 400 Hz和$\hat A$=1.0×10-14 s/s。图 5画出了KF前后z(t)和$\hat x\left(t \right)$的Allan偏差，以及通过参数估计值计算得到的WPM、WFM、RWFM、线性频漂、周期性波动项的Allan偏差。
![]()
图 5 z(t)和$\hat x\left(t \right)$的Allan偏差Figure 5. Allan Deviations of z(t) and $\hat x\left(t \right)$把§3.1~§3.4得到的参数估计值代入式（14）和式（15），得出τ < 1 d的实验结果和国产氢钟的说明书相符；当τ > 10 000 s时，观测噪声对Allan偏差的影响很小。
综上，本节实验结果验证了本文方法可以有效分析WPM、WFM、RWFM、线性频漂、周期性波动项各自的Allan偏差，以及总的Allan偏差，并通过这些估计值，拟合出任意平滑时间的Allan偏差估计值。
4 结语
本文展示了原子钟模型和频率稳定度分析方法，详细分析了原子钟时差观测量中的各分量，包括确定性部分（时差、频差、线性频漂和周期性波动项）、随机性部分（WFM、RWFM）和观测噪声（WPM）；分析了WPM、WFM、RWFM、线性频漂、周期性波动项在Allan偏差中的表达式，描述了KF用于估计原子钟状态的原理；提出了当在对数Allan偏差图中，WFM完全淹没于WPM时，使用KF估计WPM、WFM、RWFM强度的方法；提出了3种估计线性频漂幅度的方法和估计周期性波动周期和幅度的方法。通过两台国产氢钟的实测数据验证了本文方法的实用性。实际上，可以通过这些估计值拟合出任意平滑时间的Allan偏差估计值。本文提出的方法物理原理清晰，操作简便易行。该研究对于时间尺度、钟差预测、原子钟驾驭等算法具有重要意义。
-
![]()
图 2 直方图及真实和估计的分布曲线(算例1)
Figure 2. Histogram and True and Estimated Distribution Curves (Example 1)
![]()
图 3 模型参数估计的均方根误差(算例1)
Figure 3. Root Mean Square Error of Model Parameter Estimates(Example 1)
![]()
图 5 直方图及真实和估计的分布曲线(算例2)
Figure 5. Histogram and True and Estimated Distribution Curves (Example 2)
![]()
图 6 模型参数估计的均方根误差(算例2）
Figure 6. Root Mean Square Error of Model Parameter Estimates (Example 2)
表 1 高斯分布混合下的EM算法估计结果
Table 1 Estimation Results of EM Algorithm Under Gaussian Distribution Mixing
次数 EM_p EM_G EM_p EM_G EEM_p EM_G EM_p EM_G EM_p EM_G EM_p EM_G 1 3.019 1 2.904 8 3.946 3 4.072 0 0.992 8 1.171 2 2.903 1 2.708 0 0.601 7 0.585 7 0.398 3 0.414 3 2 3.019 1 2.881 8 3.946 8 4.099 2 0.993 1 1.126 7 2.903 5 2.717 4 0.601 9 0.583 9 0.398 1 0.416 1 3 3.019 1 2.904 7 3.946 0 4.070 0 0.992 7 1.170 2 2.902 9 2.707 2 0.601 6 0.584 9 0.398 4 0.415 1 4 3.019 1 2.880 4 3.945 8 4.089 7 0.992 6 1.122 2 2.902 8 2.712 8 0.601 5 0.580 0 0.398 5 0.420 0 5 3.019 0 2.878 7 3.946 4 4.079 2 0.992 9 1.117 0 2.903 0 2.707 5 0.601 7 0.575 5 0.398 3 0.424 5 6 3.019 1 2.879 4 3.945 8 4.083 2 0.992 6 1.119 1 2.902 8 2.709 5 0.601 5 0.577 2 0.398 5 0.422 8 7 3.019 0 2.885 7 3.946 4 4.125 3 0.992 9 1.138 9 2.903 0 2.729 8 0.601 7 0.594 5 0.398 3 0.405 5 8 3.019 1 2.888 4 3.946 0 4.143 6 0.992 7 1.147 3 2.902 9 2.738 3 0.601 6 0.601 7 0.398 4 0.398 3 9 3.019 0 2.886 3 3.946 4 4.129 4 0.992 9 1.140 7 2.903 1 2.731 7 0.601 7 0.596 1 0.398 3 0.403 9 10 3.019 1 2.883 4 3.946 6 4.110 0 0.993 0 1.131 7 2.903 4 2.722 5 0.601 8 0.588 3 0.398 2 0.411 7 11 3.019 0 2.885 3 3.946 7 4.122 6 0.993 0 1.137 5 2.903 3 2.728 5 0.601 8 0.593 4 0.398 2 0.406 6 12 3.019 1 2.887 7 3.946 5 4.138 2 0.993 0 1.144 9 2.903 3 2.735 7 0.601 8 0.599 6 0.398 2 0.400 4 均值 3.019 1 2.887 2 3.945 8 4.105 2 0.992 8 1.139 0 2.903 1 2.720 7 0.601 7 0.588 4 0.398 3 0.411 6 
下载: 导出CSV
表 2 拉普拉斯分布与高斯分布混合下的EM算法估计结果
Table 2 Estimation Results of EM Algorithm Under the Mixture of Laplace and Gaussian Distributions
次数 EM_p EM_G EM_p EM_G EM_p EM_G EM_p EM_G EM_p EM_G EM_p EM_G 1 0.003 8 0.118 4 0.929 9 0.883 8 0.889 5 0.644 7 4.054 1 3.744 0 0.448 8 0.406 1 0.551 2 0.593 9 2 0.003 8 0.118 5 0.929 8 0.883 9 0.889 4 0.645 2 4.053 9 3.744 5 0.448 8 0.406 3 0.551 2 0.593 7 3 0.003 8 0.118 3 0.929 7 0.883 7 0.889 3 0.644 7 4.053 8 3.743 9 0.448 7 0.406 1 0.551 3 0.593 9 4 0.003 9 0.118 4 0.930 2 0.883 8 0.889 8 0.644 9 4.054 7 3.744 2 0.449 0 0.406 2 0.551 0 0.593 8 5 0.003 9 0.118 3 0.930 0 0.883 7 0.889 6 0.644 6 4.054 3 3.743 8 0.448 9 0.406 0 0.551 1 0.594 0 6 0.003 8 0.118 3 0.929 9 0.883 7 0.889 5 0.644 6 4.054 1 3.743 8 0.448 8 0.406 0 0.551 2 0.594 0 7 0.003 9 0.118 5 0.930 1 0.883 9 0.889 8 0.645 1 4.054 6 3.744 4 0.449 0 0.406 2 0.551 0 0.593 8 8 0.003 9 0.118 4 0.930 0 0.883 8 0.889 7 0.644 9 4.054 4 3.744 1 0.448 9 0.406 2 0.551 1 0.593 8 9 0.003 9 0.118 5 0.930 2 0.884 0 0.889 9 0.645 4 4.054 8 3.744 7 0.449 0 0.406 4 0.551 0 0.593 6 10 0.003 9 0.118 5 0.930 2 0.883 9 0.889 9 0.645 3 4.054 8 3.744 6 0.449 0 0.406 3 0.551 0 0.593 7 11 0.003 8 0.118 4 0.929 7 0.883 8 0.889 3 0.645 0 4.053 7 3.744 3 0.448 7 0.406 2 0.551 3 0.593 8 12 0.003 8 0.118 3 0.929 7 0.883 7 0.889 3 0.644 6 4.053 8 3.743 7 0.448 7 0.406 0 0.551 3 0.594 0 均值 0.003 9 0.118 4 0.930 0 0.883 8 0.889 6 0.644 9 4.054 2 3.744 2 0.448 9 0.406 2 0.551 1 0.593 8
下载: 导出CSV
表 3 观测值残差的EM算法估计结果
Table 3 Estimation Results of EM Algorithm of Observed Value Residuals
次数 1 -0.066 -0.028 0.543 0.425 0.488 0.511 2 -0.075 -0.024 0.556 0.436 0.433 0.566 3 -0.065 -0.029 0.543 0.428 0.485 0.515 4 -0.071 -0.026 0.549 0.431 0.464 0.536 5 -0.075 -0.025 0.556 0.437 0.434 0.566 6 -0.073 -0.025 0.552 0.434 0.449 0.551 7 -0.067 -0.028 0.545 0.430 0.478 0.523 8 -0.070 -0.026 0.549 0.431 0.467 0.534 9 -0.076 -0.025 0.557 0.437 0.433 0.567 10 -0.068 -0.027 0.546 0.427 0.480 0.520 11 -0.072 -0.025 0.552 0.433 0.452 0.548 12 -0.071 -0.026 0.551 0.432 0.458 0.542 均值 -0.071 -0.026 0.550 0.432 0.460 0.540
下载: 导出CSV
表 4 伪距单点定位结果/m
Table 4 Pseudorange Single Point Positioning Results/m
坐标 LS Lp 估值 中误差 估值 中误差 X -1 304 154.105 2.582 -1 304 153.996 2.472 Y -4 831 834.409 4.474 -4 831 833.913 4.154 Z 3 943 236.838 5.344 3 943 236.658 5.212 精度 1.138 0.920
下载: 导出CSV
-
[1] 潘雄. 半参数模型的估计理论及其应用[D]. 武汉: 武汉大学, 2005 Pan Xiong. The Estimation Theory and Application Research in Semi-Parametric Model[D]. Wuhan: Wuhan University, 2005
[2] 潘雄, 程少杰, 赵春茹. 一元p范分布的参数快速估计方法[J]. 武汉大学学报·信息科学版, 2010, 35(2): 189-192 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201002017.htm Pan Xiong, Cheng Shaojie, Zhao Chunru. A Fast Parameter Estimation in p-Norm Distribution[J]. Geomatics and Information Science of Wuhan University, 2010, 35(2): 189-192 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201002017.htm
[3] 孙海燕. p范分布理论及其在现代测量数据处理中的应用[D]. 武汉: 武汉测绘科技大学, 1995 Sun Haiyan. p-Distribution Theory and Its Application in Modern Survey Data Processing[D]. Wuhan: Wuhan University of Surveying and Mapping, 1995
[4] Booth J G, Hobert J P. Maximizing Generalized Linear Mixed Model Likelihoods with an Automated Monte Carlo EM Algorithm[J]. Journal of the Royal Statistical Society, 1999, 61(1): 265-285 doi: 10.1111/1467-9868.00176
[5] 连军艳. EM算法及其改进在混合模型参数估计中的应用研究[D]. 西安: 长安大学, 2006 Lian Junyan. The Application Research of EM Algorithm and Its Improvement in Mixed Model Parameter Estimation[D]. Xi'an: Chang'an University, 2006
[6] Tuaç Y, Güney Y, Arslan O. Parameter Estimation of Regression Model with AR(p)Error Terms Based on Skew Distributions with EM Algorithm [J]. Soft Computing, 2020, 24(5): 3309-3330 doi: 10.1007/s00500-019-04089-x
[7] 吴柯, 何坦, 杨叶涛. 基于混合像元分解与EM算法的中低分辨率遥感影像变化检测[J]. 武汉大学学报·信息科学版, 2019, 44(4): 555-562 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201904012.htm Wu Ke, He Tan, Yang Yetao. Change Detection Method Based on Pixel Unmixing and EM Algorithm for Low and Medium Resolution Remote Sensing Imagery[J]. Geomatics and Information Science of Wuhan University, 2019, 44(4): 555-562 https://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201904012.htm
[8] 肖琴琴, 宋迎春, 杜琨. EM算法在广播星历计算卫星位置中的应用[J]. 测绘工程, 2013, 22 (6): 73-76 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201306020.htm Xiao Qinqin, Song Yingchun, Du Kun. Application of EM Algorithm to the Calculation of the Satellite Position Based on Broadcast Ephemeris[J]. Engineering of Surveying and Mapping, 2013, 22(6): 73-76 https://www.cnki.com.cn/Article/CJFDTOTAL-CHGC201306020.htm
[9] 鲁纳纳, 余旌胡. EM算法的参数分辨率[J]. 数学物理学报, 2019, 39(3): 638-648 doi: 10.3969/j.issn.1003-3998.2019.03.021 Lu Nana, Yu Jinghu. Research on Resolution Based on EM Algorithm[J]. Acta Mathematica Scientia, 2019, 39(3): 638-648 doi: 10.3969/j.issn.1003-3998.2019.03.021
[10] 赵杨璐, 段丹丹, 胡饶敏, 等. 基于EM算法的混合模型中子总体个数的研究[J]. 数理统计与管理, 2020, 39 (1): 35-50 https://www.cnki.com.cn/Article/CJFDTOTAL-SLTJ202001005.htm Zhao Yanglu, Duan Dandan, Hu Raomin, et al. On the Number of Components in Mixture Model Based on EM Algorithm[J]. Journal of Applied Statistics and Management, 2020, 39(1): 35-50 https://www.cnki.com.cn/Article/CJFDTOTAL-SLTJ202001005.htm
[11] 李仁忠, 张缓缓, 景军锋, 等. 基于EM算法的高斯混合型的织物疵点检测研究[J]. 计算机工程与应用, 2014, 50(10): 184-187 https://www.cnki.com.cn/Article/CJFDTOTAL-JSGG201410040.htm Li Renzhong, Zhang Huanhuan, Jing Junfeng, et al. Fabric Defect Detection Based on Gaussian Mixture Models of EM Algorithm[J]. Computer Engineering and Applications, 2014, 50(10): 184-187 https://www.cnki.com.cn/Article/CJFDTOTAL-JSGG201410040.htm
[12] 冯杭, 王胜兵. 基于EM算法的离散-连续型混合分布参数估计[J]. 统计与决策, 2019, 35(3): 85-88 https://www.cnki.com.cn/Article/CJFDTOTAL-TJJC201903020.htm Feng Hang, Wang Shengbing. Discrete-Continuous Mixed Distribution Parameter Estimation Based on EM Algorithm[J]. Statistics & Decision, 2019, 35 (3): 85-88 https://www.cnki.com.cn/Article/CJFDTOTAL-TJJC201903020.htm
[13] Guo X, Li Q Y, Xu W L. Acceleration of the EM Algorithm Using the Vector Aitken Method and Its Steffensen Form[J]. Acta Mathematicae Applicatae Sinica, English Series, 2017, 33(1): 175-182
[14] 潘雄, 赵启龙, 王俊雷, 等. 一元非对称p范分布的极大似然平差[J]. 测绘学报, 2011, 40(1): 33-36 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201101009.htm Pan Xiong, Zhao Qilong, Wang Junlei, et al. Maximum Likelihood Adjustment of the Monadic Unsym metrical P-Norm Distribution[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(1): 33-36 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201101009.htm
[15] 潘雄, 罗静, 汪耀. p范分布的实数阶与对数矩估计法[J]. 测绘学报, 2016, 45(3): 302-309 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201603009.htm Pan Xiong, Luo Jing, Wang Yao. Real Order and Logarithmic Moment Estimation Method of p-Norm Distribution[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(3): 302-309 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201603009.htm
-
期刊类型引用(21)
1. 李国俊，史丰丰，徐金锋，赵润，付桂涛，杨玉婷，李兆南，孙文龙，程梦飞，林勇昕. 基于钙离子光钟的160 d连续守时性能分析. 光学学报. 2025(04): 126-133 . 
百度学术
2. 崔海波，王伟，武建锋，李帅辰，李冲，王康. 机动状态下国产铯原子钟稳定性分析. 时间频率学报. 2024(01): 26-33 . 百度学术
3. 宋鑫，张秀荣，李卫东. 关联噪声对维纳过程和拉比振荡的影响. 山西大学学报(自然科学版). 2024(03): 578-582 . 百度学术
4. 李方能，梁益丰，许江宁，吴苗，朱兵. 基于鲁棒双参数指数平滑法的BDS卫星钟差预报. 中国惯性技术学报. 2024(07): 645-653 . 百度学术
5. 王芳敏，李汶林，戴鸿飞，李春怡，周建华，薛申辉，王波. A real-time performance improvement method for composite time scale. Chinese Physics B. 2024(09): 354-361 . 百度学术
6. 李宗源，袁海波，张虹，张继海，刘素芳，王一桁. 基于低通滤波器的原子钟时差测量数据降噪处理及其性能分析. 时间频率学报. 2024(03): 171-179 . 百度学术
7. 李昂，周铁中，薛潇博，易航，陈德好. 一种智能整定氢原子频标PID控制系统参数的方法. 宇航计测技术. 2024(05): 39-44 . 百度学术
8. 曹建峰，孔静，满海钧，鞠冰，张宇，刘荟萃. 天问一号星载USO的长期漂移标定. 武汉大学学报(信息科学版). 2024(12): 2181-2186 . 百度学术
9. 赵军华，张玉军，刘安斐，展昕，陈亮. 附加钟特性的原子钟预测不确定度分析. 测绘科学. 2024(11): 30-39 . 百度学术
10. 阚昊宇，胡志刚，吕逸飞，谢新，周仁宇，赵齐乐. 利用不同时间同步体制钟差评估北斗三号星载原子钟性能. 武汉大学学报(信息科学版). 2023(04): 604-610 . 百度学术
11. 李方能，梁益丰，许江宁，吴苗. BDS/GPS新型铷原子钟长期特性分析. 中国惯性技术学报. 2023(05): 452-461 . 百度学术
12. 李帅辰，武建锋，崔海波，方婧. 国产氢原子钟移动守时性能测试与分析. 导航定位学报. 2023(03): 38-44 . 百度学术
13. 郭文飞，朱萌萌，辜声峰，左鸿铭，陈金鑫. GNSS精密时频接收机时钟调控模型与参数设计方法. 武汉大学学报(信息科学版). 2023(07): 1126-1133 . 百度学术
14. 安景浩，王怀智，梁益丰，许江宁. BDS快速钟差预报方法对比分析. 舰船电子工程. 2023(03): 55-60 . 百度学术
15. 卢鋆，武建峰，袁海波，申建华，孟轶男，宿晨庚，陈颖. 北斗三号系统时频体系设计与实现. 武汉大学学报(信息科学版). 2023(08): 1340-1348 . 百度学术
16. 秦晓伟，杜二旺，赵春勃，王国永，杨涛. 基于频率源稳定度的钟差序列生成方法研究. 空间电子技术. 2023(04): 66-71 . 百度学术
17. 陈汗龙，董哲，周真帆，郑晓雪. 国产原子钟频率稳定度评估分析. 导航定位学报. 2023(06): 28-33 . 百度学术
18. 梁益丰，许江宁，吴苗，何泓洋. 北斗卫星钟差的CEEMDAN分解与周期项提取方法. 中国惯性技术学报. 2022(04): 476-484 . 百度学术
19. 张杰梁，陈鲤文，陈静，张煌辉，姜立斌. 守时系统钟组规模技术验证. 电子质量. 2022(10): 151-157 . 百度学术
20. 伍贻威，王世超. 不通过自由纸面时建立时间基准的方法与性能分析. 测绘学报. 2021(03): 343-355 . 百度学术
21. 惠恬，赵高长，苏军，龚莹莹. 原子钟数据的改进经验模态分解降噪. 河南科技大学学报(自然科学版). 2021(05): 45-50+56+7 . 百度学术
其他类型引用(11)
计量
- 文章访问数: 1023
- HTML全文浏览量: 230
- PDF下载量: 65
- 被引次数: 32
