CN106100609A - 单状态变量和两级Kalman滤波器时间尺度算法 - Google Patents

Abstract

本发明提供了一种单状态变量和两级Kalman滤波器时间尺度算法，其中一种单状态变量Kalman滤波器时间尺度算法，包括以下步骤，首先建立原子钟模型；计算每台钟的过冲状态噪声方差；获取N‑1组观测钟差；用Kalman滤波器通过N‑1组观测钟差估计出N台钟的状态；计算时间尺度的权重和预测值；最后引入虚拟Kalman采样间隔，优化时间尺度TA。在该方法的基础上，本发明还提供了一种两级Kalman滤波器时间尺度算法，通过第一级Kalman滤波器和重构时差，获得N‑1组重构时差，然后，运行第二级的Kalman滤波器生成时间尺度TA。本发明可以通过选取虚拟Kalman采样间隔，使时间尺度的频率稳定度在任意某一个指定的平滑时间达到最优。

Description

单状态变量和两级Kalman滤波器时间尺度算法

技术领域

本发明涉及时间频率、信号处理领域，具体的说是设计了一种单状态变量Kalman滤波器时间尺度算法和一种两级Kalman滤波器时间尺度算法。

背景技术

时间尺度算法(time scale algorithms)在守时实验室和全球导航卫星系统(GNSS)中发挥着重要作用。时间尺度算法的目的是为了综合钟组中所有的原子钟，生成一个可靠性、频率稳定度更高的自由的纸面时间。单台钟和纸面时间都可以视为是一个时间尺度；没有被驾驭的时间尺度称为自由的时间尺度。所以说，使用时间尺度算法是为了生成一个时间尺度，该时间尺度一般通过它与某一台钟之间的时间差来表示；算法通常不考虑该时间尺度的频率准确度，将提升其频率准确度的任务留给驾驭算法来完成。

对于一个守时实验室，该自由的纸面时间通常记为TA(k)(其中k为实验室代号)；可以用于监测物理钟的性能；也可以在使用钟差预测算法使TA(k)与UTC保持同步后，对物理钟进行驾驭，得到UTC的本地物理实现，记为UTC(k)。

对于GNSS，该自由的纸面时间即为自由的“系统时间”，记为自由的“GNSST”；例如对于北斗系统，记为自由的“BDT”。之所以打引号，是因为这时的“BDT”是自由的时间尺度，还没有与UTC(BSNC)(其中BSNC为北京卫星导航中心的代号)同步，不能称之为真正的系统时间。通过UTC(BSNC)对其进行驾驭，得到BDT，此时BDT与UTC(BSNC)在时间和频率上保持了同步。

时间尺度算法通常可以分为两大类：1)一类是经典的加权平均算法，例如：国际计量局(BIPM)采用的加权平均算法(称为：ALGOS算法)、美国国家标准与技术研究院(NIST)采用的加权平均算法(称为：AT1算法)；对于任何加权平均算法，生成的时间尺度都可以用时间尺度基本方程来表示；2)另一类是采用Kalman滤波器实现的算法。由于AT1算法和Kalman滤波器算法所生成的时间尺度是实时的，比较适用于GNSS。

上世纪80年代NIST首次采用Kalman滤波器建立时间尺度。在一段时期内，该算法和AT1算法同时在NIST运行。该算法的不足在于：1)由于系统不完全可观测，造成估计误差矩阵发散，表明状态估计的误差不断变大；2)实验表明：对于铯钟和氢钟组成的钟组，时间尺度TA主要由长期频率稳定度最好的钟决定，从而没有利用到钟组中短期频率稳定度较好的钟的短期频率稳定度性能。

发明内容

为解决现有技术存在的问题，本发明提出一种单状态变量Kalman滤波器时间尺度算法；在单状态变量Kalman滤波器时间尺度算法的基础上，本发明还提供一种两级Kalman滤波器时间尺度算法。

本发明的技术方案是：

一种单状态变量Kalman滤波器时间尺度算法，包括以下步骤：

S1建立原子钟模型

$\left\{\begin{array}{c}{X}_1\left(t\right)=x_0+y_0\·t+1/2\·d\·t^2+{\mathrm{\σ}}_1\·W_1\left(t\right)+{\mathrm{\σ}}_2\·{\∫}_0^t{W}_2\left(s\right)ds\\ X_2\left(t\right)=y_0+d\·t+{\mathrm{\σ}}_2\·W_2\left(t\right)\end{array}\right.---\left(1\right)$

其中，X₁(t)、X₂(t)分别为两个钟差的状态变量；其中X₁(t)代表瞬时时差，X₂(t)代表不含频率白噪声(WFM)的瞬时频差；x₀和y₀分别代表瞬时时差和瞬时频差的初值；d代表频漂；W₁(t)和W₂(t)分别代表两个独立的维纳过程，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布；σ₁和σ₂分别是这两个维纳过程的扩散系数(diffusioncoefficients)，用于表明噪声的强度。

扩散系数与Allan方差的关系表示为：

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)={\mathrm{\σ}}_1^2/\mathrm{\τ}+1/3\·{\mathrm{\σ}}_2^2\·\mathrm{\τ}---\left(2\right)$

其中，代表平滑时间为τ时的Allan方差(是个约定俗成的表示Allan方差的方式，下标y的含义就是表明了Allan方差。)，τ为平滑时间。

S2计算每台钟的过冲状态噪声方差

钟组的状态方程为：

X(k+1)＝φ·X(k)+J(k) (3)

其中，X(k)＝[x₁(k),x₂(k),…,x_N(k)]^T是钟组在k时刻的时差；x_i(k)(i＝1,2,…,N)是钟组中第i台钟的时差；φ＝I_N×N是状态转移矩阵；J(k)是钟组的过程噪声，表示为J(k)＝[J₁(k),J₂(k),…,J_N(k)]^T，其中J_i(k)(i＝1,2,…,N)是钟组中第i台钟的过程噪声。

由式(1)，J_i(k)(i＝1,2,…,N)可以表示为：

$J_i\left(k\right)={\mathrm{\σ}}_{i,1}\left(W_{i,1}\right(t_{k+1})-W_{i,1}(t_k\left)\right)+{\mathrm{\σ}}_{i,2}{\∫}_{t_k}^{t_{k+1}}\left(W_{i,2}\right(t_s)-W_{i,2}(t_k\left)\right)ds,(i=1,2,...,N)---\left(4\right)$

其中，W₁(t)，W₂(t)，σ₁和σ₂的含义和式(1)中相应的符号相同，下标i代表第i台钟。

第i台钟的过程状态噪声方差即为式(4)的方差，表示为：

$Q_i={\mathrm{\σ}}_{i,1}^2\·T+1/3\·{\mathrm{\σ}}_{i,2}^2\·T^3---\left(5\right)$

由式(5)和式(2)，得到：

$Q_i={\mathrm{\σ}}_{i,y}^2\left(T\right)\·T^2---\left(6\right)$

其中，表示第i台钟的Allan方差。

于是，钟组的过程噪声方差表示为：

S3获取输入量即N-1组观测钟差

钟组的观测方程为

Z(k)＝H·X(k)+V(k) (8)

其中，H表示为：

V是测量噪声，表示为V(k)＝[V₁(k),V₂(k),…,V_N-1(k)]^T，其中V_i(k)(i＝1,2,…,N-1)是每一组钟差的观测噪声，它们之间相互独立；V_i(k)的噪声协方差为R_i，所以钟组的测量噪声V的协方差表示为：

由式(8)和式(9)，观测钟差可以表示为：

Z(k)＝[x₁(k)-x₂(k)+v₁(k)…x₁(k)-x_N(k)+v_N-1(k)]^T (11)

其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差。

S4使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态

${\hat{X}}_{k,k-1}=\mathrm{\φ}\·{\hat{X}}_{k-1,k-1}---\left(12\right)$

P_k,k-1＝φ·P_k-1,k-1·φ^T+Q (13)

K_k＝P_k,k-1·H^T(H·P_k,k-1·H^T+R)^-1 (14)

${\hat{X}}_{k,k}={\hat{X}}_{k,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k,k-1})---\left(15\right)$

P_k,k＝(I-K_k·H)·P_k,k-1 (16)

其中：代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、P_k-1,k-1是钟组时差的估计误差矩阵、P_k,k-1是钟组时差的预测误差矩阵、K_k是Kalman滤波器增益。

通过式(12)至至式(16)五个方程得到是第i台钟相对于TA的时差的估计值。

S5计算时间尺度的权重和预测值

将式(12)代入式(15)，得到：

${\hat{X}}_{k,k}={\hat{X}}_{k-1,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k-1,k-1})---\left(17\right)$

式(17)可以写成如下形式：

将每台钟相对于TA的时间差减去其时间差的估计值得到即为校准钟，由式(18)得到第1台校准钟表示为：

$\begin{array}{c}{x}_1\left(k\right)-{\hat{x}}_1\left(k\right)=\[1-K_{11}\left(k\right)-K_{12}\left(k\right)...-K_{1(N-1)}\left(k\right)\]\·\[x_1\left(k\right)-\widehat{x_1}(k-1)\]\\ +K_{11}\left(k\right)\·\[x_2\left(k\right)-{\hat{x}}_2(k-1)\]+K_{12}\left(k\right)\·\[x_3\left(k\right)-\widehat{x_3}(k-1)\]+...\\ +K_{1(N-1)}\left(k\right)\·\[x_N\left(k\right)-{\hat{x}}_N(k-1)\]-\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{1i}{v}_i\left(k\right)\end{array}---\left(19\right)$

当观测噪声为零时，式(19)满足时间尺度基本方程，即：

$TA\left(k\right)=x_1\left(k\right)-{\hat{x}}_1\left(k\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{1,j}\·\left(x_j\right(k)-{x_j}^{\′}(k\left)\right)---\left(20\right)$

其中，时间尺度的权重表示为：

$\begin{array}{ccc}{\mathrm{\ω}}_{1,1}=1-\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{1j}\left(k\right)& {\mathrm{\ω}}_{1,i}=K_{1(j-1)}\left(k\right)& (j\≠1)\end{array}---\left(21\right)$

时间尺度的预测值表示为：

${x_j}^{\′}\left(k\right)={\hat{x}}_j(k-1)---\left(22\right)$

当测量噪声为零时，由式(18)，任何一台“校准钟”都满足时间尺度基本方程，即：

$TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{i,j}\·\[x_j\left(k\right)-{\hat{x}}_j(k-1)\]---\left(23\right)$

由式(18)，当i≠1时，时间尺度的权重表示为：

$\begin{array}{cccc}{\mathrm{\ω}}_{i,1}=-\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{ij}\left(k\right)& {\mathrm{\ω}}_{i,i}=1+K_{i(i-1)}\left(k\right)& {\mathrm{\ω}}_{i,j}=K_{i(j-1)}\left(k\right)& (j\≠1,j\≠i)\end{array}---\left(24\right)$

而时间尺度的预测值的表达式与式(22)相同。

S6，引入虚拟Kalman采样间隔，优化时间尺度TA，方法为：

引入虚拟Kalman采样间隔T_v用于计算得到每台钟的虚拟状态噪声协方差，表示为：

${\mathrm{Qv}}_i={\mathrm{\σ}}_{i,1}^2\·T_v+1/3\·{\mathrm{\σ}}_{i,2}^2\·{T_v}^3,\∀i(i=1,...,N).---\left(25\right)$

按照步骤S2至S5的方法，TA将表示为：

$\begin{array}{c}TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right),(\∀i)\\ =\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{Qv}}_j}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{Qv}}_j}\·\[x_j\left(k\right)-{\hat{x}}_j(k-1)\]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}\·\[x_j\left(k\right)-\widehat{x_j}(k-1)\]\end{array}.---\left(26\right)$

本发明还提供一种两级Kalman滤波器时间尺度算法，包括以下步骤：

S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计，得到每组钟差的两个状态的估计量，记为和

N-1个Kalman滤波器组成第一级Kalman滤波器，第一级Kalman滤波器其输入量分别为N-1组观测钟差，把其中每组钟差分别记为Z^(1,n)(k)，其中上标(1,n)表示第1台钟与第n台钟之间的钟差，第i台钟的时差记为x_i(k)，通过这N-1个Kalman滤波器，估计得到N-1组钟差的两个状态的估计量；

S2、对时差进行重构

对于每一组钟差估计值，令

${{\tilde{X}}_1}^{(1,n)}(k+1)={{\tilde{X}}_1}^{(1,n)}\left(k\right)+{{\hat{X}}_2}^{(1,n)}\left(k\right)\·T---\left(27\right)$

把序列的初始值取为零，上式可以写为：

${{\tilde{X}}_1}^{(1,n)}\left(k\right)=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{{\hat{X}}_2}^{(1,n)}\left(m\right)\·T---\left(28\right)$

于是，式(28)就是通过估计量进行重构，得到的第1台钟和第n台钟之间的重构时差

把第i台钟的重构钟的时差记为于是有其中为重构时差；

S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成TA；

输入量为N-1组重构时差输出为N台重构钟的时差估计值，记为代表了第i台重构钟的时差与TA之间的偏差，于是得到：

$TA\left(k\right)={\tilde{x}}_i\left(k\right)-{\hat{\tilde{x}}}_i\left(k\right),\∀i.---\left(29\right)$

本发明提出一种单状态变量Kalman滤波器时间尺度算法。其中，每台钟只有1个状态变量，即每台钟的时差。在获得了权重的解析表达式的基础上，引入了虚拟Kalman采样间隔，用于计算过程噪声协方差；并证明在这种情况下，权重将反比于平滑时间为虚拟Kalman采样间隔的Allan方差；这样，可以通过选取虚拟Kalman采样间隔，使时间尺度的频率稳定度在任意某一个指定的平滑时间达到最优。

在此基础上，本发明还提出了一种两级Kalman滤波器时间尺度算法。第一级的Kalman滤波器用于估计钟差的两个状态。其中，第一个状态中包含了频率白噪声(WFM)和频率随机游走噪声(RWFM)两种噪声；而第二个状态中只包含了RWFM。第二级的Kalman滤波器实际上是一个单状态变量Kalman滤波器时间尺度算法，用于生成时间尺度；其输入为第一级Kalman滤波器的第二个状态估计值的“重构时差”；这样，算法相当于对只含有RWFM的时间序列进行加权，生成的时间尺度中也只含有RWFM这一种噪声，中短期频率稳定度更高。

附图说明

图1图解说明采用第一种算法时，在第一种情况下，每台钟、TA的时差(a)和Allan偏差(b)；

图2图解说明采用第一种算法时，在第二种情况下，每台钟、TA的时差(a)和Allan偏差(b)；

图3图解说明两级Kalman滤波器算法原理；

图4图解说明采用第二种算法时，每台“重构钟”、TA的时差(a)和Allan偏差(b)；

具体实施方式

下面结合附图对本发明的具体实施方式作出说明。

首先，时间尺度基本方程表示为：

$TA\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\left(x_i\right(t)-{x_i}^{\′}(t\left)\right)$

其中，TA(t)是生成的纸面上的时间尺度；N为钟组中原子钟的数量；ω_i为第i台原子钟的权重；x_i(t)代表第i台原子钟的时差；x_i’(t)代表第i台原子钟的时差的预测值。

时间尺度基本方程可以这样理解：对于加权平均时间尺度算法，权重用于优化时间尺度的频率稳定度，而预测值则消除了每台钟确定性趋势项的影响；所以TA的物理意义是对去除确定性趋势项后只含有噪声分量的原子钟时差的加权平均值。

事实上，不论是ALGOS算法还是AT1算法，其生成的时间尺度都可以进行化简，最终使用时间尺度基本方程来表示。

本发明提供了一种单状态变量Kalman滤波器时间尺度算法，包括以下步骤：

S1建立原子钟模型

$\left\{\begin{array}{c}{X}_1\left(t\right)=x_0+y_0\·t+1/2\·d\·t^2+{\mathrm{\σ}}_1\·W_1\left(t\right)+{\mathrm{\σ}}_2\·{\∫}_0^t{W}_2\left(s\right)ds\\ X_2\left(t\right)=y_0+d\·t+{\mathrm{\σ}}_2\·W_2\left(t\right)\end{array}\right.---\left(1\right)$

其中，X₁(t)、X₂(t)分别为两个钟差的状态变量；其中X₁(t)代表瞬时时差，X₂(t)代表不含频率白噪声(WFM)的瞬时频差；x₀和y₀分别代表瞬时时差和瞬时频差的初值；d代表频漂；W₁(t)和W₂(t)分别代表两个独立的维纳过程，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布；σ₁和σ₂分别是这两个维纳过程的扩散系数(diffusioncoefficients)，用于表明噪声的强度。

式(1)表明：W₁(t)和W₂(t)对时间t的积分在状态变量X₁上分别表现为WFM和RWFM；即式(1)的第一个方程中，最后第二项代表WFM，最后一项代表RWFM。

扩散系数与Allan方差的关系表示为：

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)={\mathrm{\σ}}_1^2/\mathrm{\τ}+1/3\·{\mathrm{\σ}}_2^2\·\mathrm{\τ}---\left(2\right)$

其中，代表平滑时间为τ时的Allan方差，τ为平滑时间.式(2)中等号后面第一项的斜率为-1，第二项斜率为1，这两项分别对应WFM和RWFM，这说明在对数Allan方差图中，WFM的斜率为-1，RWFM的斜率为1。

S2计算每台钟的过冲状态噪声方差

钟组的状态方程为：

X(k+1)＝φ·X(k)+J(k) (3)

其中，X(k)＝[x₁(k),x₂(k),…,x_N(k)]^T是钟组在k时刻的时差；x_i(k)(i＝1,2,…,N)是钟组中第i台钟的时差；φ＝I_N×N是状态转移矩阵；J(k)是钟组的过程噪声，表示为J(k)＝[J₁(k),J₂(k),…,J_N(k)]^T，其中J_i(k)(i＝1,2,…,N)是钟组中第i台钟的过程噪声；由式(1)，J_i(k)(i＝1,2,…,N)可以表示为：

$J_i\left(k\right)={\mathrm{\σ}}_{i,1}\left(W_{i,1}\right(t_{k+1})-W_{i,1}(t_k\left)\right)+{\mathrm{\σ}}_{i,2}{\∫}_{t_k}^{t_{k+1}}\left(W_{i,2}\right(t_s)-W_{i,2}(t_k\left)\right)ds,(i=1,2,...,N)---\left(4\right)$

第i台钟的过程状态噪声方差即为式(4)的方差，表示为：

$Q_i={\mathrm{\σ}}_{i,1}^2\·T+1/3\·{\mathrm{\σ}}_{i,2}^2\·T^3---\left(5\right)$

由式(5)和式(2)，得到：

$Q_i={\mathrm{\σ}}_{i,y}^2\left(T\right)\·T^2---\left(6\right)$

其中，表示第i台钟的Allan方差；

于是，钟组的过程噪声方差表示为：

S3获取输入量即N-1组观测钟差

钟组的观测方程为

Z(k)＝H·X(k)+V(k) (8)

其中，H表示为：

V是测量噪声，表示为V(k)＝[V₁(k),V₂(k),…,V_N-1(k)]^T，其中V_i(k)(i＝1,2,…,N-1)是每一组钟差的观测噪声，它们之间相互独立；V_i(k)的噪声协方差为R_i，所以钟组的测量噪声V的协方差表示为：

由式(8)和式(9)，观测钟差可以表示为：

Z(k)＝[x₁(k)-x₂(k)+v₁(k)…x₁(k)-x_N(k)+v_N-1(k)]^T (11)

其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差；

S4使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态

Kalman滤波器应用于时间尺度建立时，Kalman滤波器的作用是通过N-1组观测钟差估计出N台钟的状态，记为和所有Kalman滤波器一样，此时Kalman滤波器的步骤可以由下面5个方程来描述。

${\hat{X}}_{k,k-1}=\mathrm{\φ}\·{\hat{X}}_{k-1,k-1}---\left(12\right)$

P_k,k-1＝φ·P_k-1,k-1·φ^T+Q (13)

K_k＝P_k,k-1·H^T(H·P_k,k-1·H^T+R)^-1 (14)

${\hat{X}}_{k,k}={\hat{X}}_{k,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k,k-1})---\left(15\right)$

P_k,k＝(I-K_k·H)·P_k,k-1 (16)

其中：代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、P_k-1,k-1是钟组时差的估计误差矩阵、P_k,k-1是钟组时差的预测误差矩阵、K_k是Kalman滤波器增益。

通过式(12)至至式(16)五个方程得到是第i台钟相对于TA的时差的估计值。也就是说：Kalman估计得到的量是每台钟相对于TA的时间差。

S5计算时间尺度的权重和预测值

将式(12)代入式(15)，得到：

${\hat{X}}_{k,k}={\hat{X}}_{k-1,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k-1,k-1})---\left(17\right)$

式(17)可以写成如下形式：

将每台钟相对于TA的时间差减去其时间差的估计值得到即为校准钟，由式(18)得到第1台校准钟表示为：

$\begin{array}{c}{x}_1\left(k\right)-{\hat{x}}_1\left(k\right)=\[1-K_{11}\left(k\right)-K_{12}\left(k\right)...-K_{1(N-1)}\left(k\right)\]\·\[x_1\left(k\right)-\widehat{x_1}(k-1)\]\\ +K_{11}\left(k\right)\·\[x_2\left(k\right)-{\hat{x}}_2(k-1)\]+K_{12}\left(k\right)\·\[x_3\left(k\right)-\widehat{x_3}(k-1)\]+...\\ +K_{1(N-1)}\left(k\right)\·\[x_N\left(k\right)-{\hat{x}}_N(k-1)\]-\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{1i}{v}_i\left(k\right)\end{array}---\left(19\right)$

当观测噪声为零时，式(19)满足时间尺度基本方程，即：

$TA\left(k\right)=x_1\left(k\right)-{\hat{x}}_1\left(k\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{1,j}\·\left(x_j\right(k)-{x_j}^{\′}(k\left)\right)---\left(20\right)$

其中，时间尺度的权重表示为：

$\begin{array}{ccc}{\mathrm{\ω}}_{1,1}=1-\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{1j}\left(k\right)& {\mathrm{\ω}}_{1,i}=K_{1(j-1)}\left(k\right)& (j\≠1)\end{array}---\left(21\right)$

时间尺度的预测值表示为：

${x_j}^{\′}\left(k\right)={\hat{x}}_j(k-1)---\left(22\right)$

当测量噪声为零时，由式(18)，任何一台“校准钟”都满足时间尺度基本方程，即：

$TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right)=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{i,j}\·\[x_j\left(k\right)-{\hat{x}}_j(k-1)\]---\left(23\right)$

由式(18)，当i≠1时，时间尺度的权重表示为：

$\begin{array}{cccc}{\mathrm{\ω}}_{i,1}=-\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{ij}\left(k\right)& {\mathrm{\ω}}_{i,i}=1+K_{i(i-1)}\left(k\right)& {\mathrm{\ω}}_{i,j}=K_{i(j-1)}\left(k\right)& (j\≠1,j\≠i)\end{array}---\left(24\right)$

而时间尺度的预测值的表达式与式(22)相同。

S6，引入虚拟Kalman采样间隔，优化时间尺度TA，方法为：

引入虚拟Kalman采样间隔T_v用于计算得到每台钟的虚拟状态噪声协方差，表示为：

${\mathrm{Qv}}_i={\mathrm{\σ}}_{i,1}^2\·T_v+1/3\·{\mathrm{\σ}}_{i,2}^2\·{T_v}^3,\∀i(i=1,...,N).---\left(25\right)$

按照步骤S2至S5的方法，TA将表示为：

$\begin{array}{c}TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right),(\∀i)\\ =\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{Qv}}_j}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{Qv}}_j}\·\[x_j\left(k\right)-{\hat{x}}_j(k-1)\]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}\·\[x_j\left(k\right)-\widehat{x_j}(k-1)\]\end{array}.---\left(26\right)$

式(26)表明：在采用虚拟Kalman采样间隔T_v时，算法可以使TA在平滑时间为T_v的频率稳定度达到了最优。因此，在实际运用中，可以直接选取一个T_v的值，这样生成的TA就在指定的平滑时间T_v的频率稳定度达到了最优。

采用上述方法进行仿真验证：

首先仿真生成2台氢钟和1台铯钟组成钟组(每台钟共5000个数据点，相邻数据点之间的采样间隔为3600s)。每台氢钟如式(1)所示的平方扩散系数设为σ₁ ²＝5.0×10^-25s和σ₂ ²＝1.3×10^-35s^-1；铯钟的平方扩散系数设为σ₁ ²＝4.8×10^-23s和σ₂ ²＝1.9×10^-36s^-1。这些扩散系数符合典型氢钟和铯钟的指标。设置观测噪声为零。为了简化分析，设每台钟如式(1)所示的确定性分量为x₀＝y₀＝d＝0。

然后，按照上述方法生成TA。分为两种情况。第一种情况，设置虚拟Kalman采样间隔T_v＝432000s(5天)；第二种情况，设置T_v＝2592000s(30天)。当观测噪声为零时，所有“校准钟”都是相等的，所以不妨把第1台“校准钟”作为TA。

表1列出了在这两种情况下，由式(26)计算得到的每台钟权重的理论值。图1描绘了在第一种情况下，这三台钟、TA的时差(a)和Allan偏差(Allan方差的平方根)(b)。图2描绘了在第二种情况下，它们的时差(a)和Allan偏差(b)。

表1两种情况下每台钟的权重

由表1、图1和图2，可以看出：在T_v较小时，氢钟的权重较大，于是TA主要由两台氢钟决定；而在T_v较大时，铯钟的权重较大，于是TA主要由铯钟决定。

本发明提供的一种单状态变量Kalman滤波器时间尺度算法实质上是一种加权平均算法；当测量噪声为零时，权重反比于每台钟在平滑时间为T_v时的Allan偏差。所以，可以根据实际需求，灵活地设置T_v的值，使TA在该平滑间隔的频率稳定度达到最优。

本发明还提供了一种两级Kalman滤波器时间尺度算法，其原理如图3所示，包括以下步骤：

S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计，得到每组钟差的两个状态的估计量，记为和

N-1个Kalman滤波器组成第一级Kalman滤波器，第一级Kalman滤波器其输入量分别为N-1组钟差观测值(如式11所示)，把其中每组钟差分别记为Z^(1,n)(k)，其中上标(1,n)表示第1台钟与第n台钟之间的钟差，第i台钟的时差记为x_i(k)，；如果没有观测噪声，显然有Z^(1,n)(k)＝x₁(k)-x_n(k)。

通过这N-1个Kalman滤波器分别对这N-1组钟差进行状态估计，得到每组钟差的如式(1)所示的两个状态的估计量，记为和由式(1)可知，时间序列中的噪声为WFM和RWFM；而时间序列中将只含有RWFM。

S2、对时差进行重构

对于每一组钟差估计值，令

${{\tilde{X}}_1}^{(1,n)}(k+1)={{\tilde{X}}_1}^{(1,n)}\left(k\right)+{{\hat{X}}_2}^{(1,n)}\left(k\right)\·T---\left(27\right)$

把序列的初始值取为零，上式可以写为：

${{\tilde{X}}_1}^{(1,n)}\left(k\right)=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{{\hat{X}}_2}^{(1,n)}\left(m\right)\·T---\left(28\right)$

于是，式(28)就是通过估计量进行重构，得到的第1台钟和第n台钟之间的重构时差显然，“重构时差”中也只含有RWFM。

把第i台钟的重构钟的时差记为于是有其中为重构时差。

通过第一级Kalman滤波器和重构时差这两个环节，相当于滤除了输入时间序列(即观测钟差)中的WFM，得到了只含有RWFM的“重构时差”。能达到这样的效果，是利用了以下特性：Kalman滤波器用于状态估计时，按照原子钟模型即式(1)，估计得到的第1个状态变量中含有WFM和RWFM，而第2个状态变量中只含有RWFM。

S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成TA；

输入量为N-1组重构时差输出为N台重构钟的时差估计值，记为代表了第i台重构钟的时差与TA之间的偏差，于是得到：

$TA\left(k\right)={\tilde{x}}_i\left(k\right)-{\hat{\tilde{x}}}_i\left(k\right),\∀i.---\left(29\right)$

采用上述方法进行仿真验证：

实际应用中的步骤应该为：1)首先使用第一级的Kalman滤波器，并对时差进行重构，获得N-1组“重构时差”，即：H1-H2和H1-Cs的“重构时差”；2)然后，运行如第二级的Kalman滤波器，得到的估计量为每台“重构钟”与TA之间的时间差。

在本仿真实验中，假设单台钟的时差是已知的。

在本次仿真中，首先使用第一级的3个Kalman滤波器，分别对3台钟进行状态估计，对时差进行重构，获得3台“重构钟”的时差，显然也就得到了2组“重构时差”。然后，运行第二级的Kalman滤波器，其输入为两组“重构时差”，输出为估计得到的每台“重构钟”与TA之间的时间差，于是也就生成了TA。其中，设置虚拟Kalman采样间隔T_v＝432000s(即5天)。

表2每台“重构钟”的权重

Table 2the weights ofeach“reconstructed clock”

表2列出了每台“重构钟”的权重。从表2中可以看出：“重构铯钟”获得的权重更大，权重反比于σ₂ ²的值。图4描述了每台“重构钟”、TA的时差(a)和Allan偏差(b)。从图4中可以看出：TA中也只含有RWFM；它的中短期频率稳定度相比第3节生成的TA得到了明显的提高。

综上所述，虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，任何本领域普通技术人员，在不脱离本发明的精神和范围内，当可作各种更动与润饰，因此本发明的保护范围当视权利要求书界定的范围为准。

Claims (3)

1.一种单状态变量Kalman滤波器时间尺度算法，其特征在于，包括以下步骤：

S1建立原子钟模型

$\left\{\begin{array}{c}{X}_1\left(t\right)=x_0+y_0\·t+1/2\·d\·t^2+{\mathrm{\σ}}_1\·W_1\left(t\right)+{\mathrm{\σ}}_2\·{\∫}_0^t{W}_2\left(s\right)ds\\ X_2\left(t\right)=y_0+d\·t+{\mathrm{\σ}}_2\·W_2\left(t\right)\end{array}\right.---\left(1\right)$

其中，X₁(t)、X₂(t)分别为观测钟差的两个状态变量；其中X₁(t)代表瞬时时差，X₂(t)代表不含频率白噪声的瞬时频差；x₀和y₀分别代表瞬时时差和瞬时频差的初值；d代表频漂；W₁(t)和W₂(t)分别代表两个独立的维纳过程，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布；σ₁和σ₂分别是这两个维纳过程的扩散系数，用于表明噪声的强度；

扩散系数与Allan方差的关系表示为：

${\mathrm{\σ}}_y^2\left(\mathrm{\τ}\right)={\mathrm{\σ}}_1^2/\mathrm{\τ}+1/3\·{\mathrm{\σ}}_2^2\·\mathrm{\τ}---\left(2\right)$

其中，代表平滑时间为τ时的Allan方差，τ为平滑时间；

S2计算每台钟的过冲状态噪声方差

钟组的状态方程为：

X(k+1)＝φ·X(k)+J(k) (3)

其中，X(k)＝[x₁(k),x₂(k),…,x_N(k)]^T是钟组在k时刻的时差；x_i(k)(i＝1,2,…,N)是钟组中第i台钟的时差；φ＝I_N×N是状态转移矩阵；J(k)是钟组的过程噪声，表示为J(k)＝[J₁(k),J₂(k),…,J_N(k)]^T，其中J_i(k)(i＝1,2,…,N)是钟组中第i台钟的过程噪声；由式(1)，J_i(k)(i＝1,2,…,N)可以表示为：

$J_i\left(k\right)={\mathrm{\σ}}_{i,1}\left(W_{i,1}\left(t_{k+1}\right)-W_{i,1}\left(t_k\right)\right)+{\mathrm{\σ}}_{i,2}{\∫}_{t_k}^{t_{k+1}}\left(W_{i,2}\left(t_s\right)-W_{i,2}\left(t_k\right)\right)ds,\left(i=1,2,...,N\right)---\left(4\right)$

第i台钟的过程状态噪声方差即为式(4)的方差，表示为：

$Q_i={\mathrm{\σ}}_{i,1}^2\·T+1/3\·{\mathrm{\σ}}_{i,2}^2\·T^3---\left(5\right)$

由式(5)和式(2)，得到：

$Q_i={\mathrm{\σ}}_{i,y}^2\left(T\right)\·T^2---\left(6\right)$

其中，表示第i台钟的Allan方差；

于是，钟组的过程噪声方差表示为：

S3获取输入量即N-1组观测钟差

钟组的观测方程为

Z(k)＝H·X(k)+V(k) (8)

其中，H表示为：

V是测量噪声，表示为V(k)＝[V₁(k),V₂(k),…,V_N-1(k)]^T，其中V_i(k)(i＝1,2,…,N-1)是每一组钟差的观测噪声，它们之间相互独立；V_i(k)的噪声协方差为R_i，所以钟组的测量噪声V的协方差表示为：

由式(8)和式(9)，观测钟差可以表示为：

Z(k)＝[x₁(k)-x₂(k)+v₁(k)…x₁(k)-x_N(k)+v_N-1(k)]^T (11)

其中的每一个量代表在k时刻第1台钟与第i台钟(i≠1)之间的钟差；

S4使用Kalman滤波器通过N-1组观测钟差估计出N台钟的状态

${\hat{X}}_{k,k-1}=\mathrm{\φ}\·{\hat{X}}_{k-1,k-1}---\left(12\right)$

P_k,k-1＝φ·P_k-1,k-1·φ^T+Q (13)

K_k＝P_k,k-1·H^T(H·P_k,k-1·H^T+R)^-1 (14)

${\hat{X}}_{k,k}={\hat{X}}_{k,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k,k-1})---\left(15\right)$

P_k,k＝(I-K_k·H)·P_k,k-1 (16)

其中：代表了钟组在k-1时刻的时差的估计值、代表了钟组在k-1时刻对k时刻的时差的预测值、P_k-1,k-1是钟组时差的估计误差矩阵、P_k,k-1是钟组时差的预测误差矩阵、K_k是Kalman滤波器增益；

通过式(12)至至式(16)五个方程得到是第i台钟相对于TA的时差的估计值；

S5计算时间尺度的权重和预测值

将式(12)代入式(15)，得到：

${\hat{X}}_{k,k}={\hat{X}}_{k-1,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k-1,k-1})---\left(17\right)$

式(17)可以写成如下形式：

将每台钟相对于TA的时间差减去其时差的估计值得到即为校准钟，由式(18)得到第1台校准钟表示为：

$\begin{array}{c}{x}_1\left(k\right)-{\hat{x}}_1\left(k\right)=\[1-K_{11}\left(k\right)-K_{12}\left(k\right)...-K_{1\left(N-1\right)}\left(k\right)\]\·\[x_1\left(k\right)-{\hat{x}}_1\left(k-1\right)\]\\ +K_{11}\left(k\right)\·\[x_2\left(k\right)-{\hat{x}}_2\left(k-1\right)\]+K_{12}\left(k\right)\·\[x_3\left(k\right)-{\hat{x}}_3\left(k-1\right)\]+...\\ +K_{1\left(N-1\right)}\left(k\right)\·\[x_N\left(k\right)-{\hat{x}}_N\left(k-1\right)\]-\underset{i=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{1i}{v}_i\left(k\right)\end{array}---\left(19\right)$

当观测噪声为零时，式(19)满足时间尺度基本方程，即：

$TA\left(k\right)=x_1\left(k\right)-{\hat{x}}_1\left(k\right)=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{1,j}\·\left(x_j\right(k)-{x_j}^{\′}(k\left)\right)---\left(20\right)$

其中，时间尺度的权重表示为：

$\begin{array}{ccc}{\mathrm{\ω}}_{1,1}=1-\underset{j=1}{\overset{N-1}{\Σ}}{K}_{1j}\left(k\right)& {\mathrm{\ω}}_{1,i}=K_{1\left(j-1\right)}\left(k\right)& (j\≠1)\end{array}---\left(21\right)$

时间尺度的预测值表示为：

${x_j}^{\′}\left(k\right)={\hat{x}}_j(k-1)---\left(22\right)$

当测量噪声为零时，由式(18)，任何一台“校准钟”都满足时间尺度基本方程，即：

$TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right)=\underset{j=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_{i,j}\·\[x_j\left(k\right)-{\hat{x}}_j(k-1)\]---\left(23\right)$

由式(18)，当i≠1时，时间尺度的权重表示为：

$\begin{array}{cccc}{\mathrm{\ω}}_{i,1}=-\underset{j=1}{\overset{N-1}{\mathrm{\Σ}}}{K}_{ij}\left(k\right)& {\mathrm{\ω}}_{i,i}=1+K_{i(i-1)}\left(k\right)& {\mathrm{\ω}}_{i,j}=K_{i(j-1)}\left(k\right)& (j\≠1,j\≠i)\end{array}---\left(24\right)$

而时间尺度的预测值的表达式与式(22)相同。

2.根据权利要求1所述的单状态变量Kalman滤波器时间尺度算法，其特征在于，还包括步骤S6，引入虚拟Kalman采样间隔，优化时间尺度TA，方法为：

引入虚拟Kalman采样间隔T_v用于计算得到每台钟的虚拟状态噪声协方差，表示为：

${\mathrm{Qv}}_i={\mathrm{\σ}}_{i,1}^2\·T_v+1/3\·{\mathrm{\σ}}_{i,2}^2\·{T_v}^3,\∀i(i=1,...,N).---\left(25\right)$

按照步骤S2至S5的方法，TA将表示为：

$\begin{array}{c}TA\left(k\right)=x_i\left(k\right)-{\hat{x}}_i\left(k\right),\left(\∀i\right)\\ =\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{Qv}}_j}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{Qv}}_j}\·\[x_j\left(k\right)-{\hat{x}}_j\left(k-1\right)\]=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\frac{1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}{\underset{j=1}{\overset{N}{\mathrm{\Σ}}}1/{\mathrm{\σ}}_{j,y}^2\left(T_v\right)}\·\[x_j\left(k\right)-{\hat{x}}_j\left(k-1\right)\]\end{array}.---\left(26\right)$

3.一种两级Kalman滤波器时间尺度算法，其特征在于，包括以下步骤：

S1、通过N-1个Kalman滤波器分别对N-1组钟差进行状态估计，得到每组钟差的两个状态的估计量，记为和

N-1个Kalman滤波器组成第一级Kalman滤波器，第一级Kalman滤波器其输入量分别为N-1组观测钟差，把其中每组钟差分别记为Z^(1,n)(k)，其中上标(1,n)表示第1台钟与第n台钟之间的钟差，第i台钟的时差记为x_i(k)，通过这N-1个Kalman滤波器，估计得到N-1组钟差的两个状态的估计量；

S2、对时差进行重构

对于每一组钟差估计值，令

${{\tilde{X}}_1}^{(1,n)}(k+1)={{\tilde{X}}_1}^{(1,n)}\left(k\right)+{{\hat{X}}_2}^{(1,n)}\left(k\right)\·T---\left(27\right)$

把序列的初始值取为零，上式可以写为：

${{\tilde{X}}_1}^{\left(1,n\right)}\left(k\right)=\underset{m=0}{\overset{k-1}{\Σ}}{{\hat{X}}_2}^{\left(1,n\right)}\left(m\right)\·T---\left(28\right)$

于是，式(28)就是通过估计量进行重构，得到的第1台钟和第n台钟之间的重构时差

把第i台钟的重构钟的时差记为于是有其中为重构时差；

S3、第二级Kalman滤波器采用如权利要求1所述的单状态变量Kalman滤波器时间尺度算法生成TA；

输入量为N-1组重构时差输出为N台重构钟的时差估计值，记为代表了第i台重构钟的时差与TA之间的偏差，于是得到：

$TA\left(k\right)={\tilde{x}}_i\left(k\right)-{\hat{\tilde{x}}}_i\left(k\right),\∀i.---\left(29\right)$