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

摘要

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

主权项

一种单状态变量Kalman滤波器时间尺度算法，其特征在于，包括以下步骤：S1建立原子钟模型$\left\{\begin{array}{c}{X}_1\left(t\right)=x_0+y_0·t+1/2·d·t^2+{\sigma }_1·W_1\left(t\right)+{\sigma }_2·{\int }_0^t{W}_2\left(s\right)ds\\ X_2\left(t\right)=y_0+d·t+{\sigma }_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方差的关系表示为：${\sigma }_y^2\left(\tau \right)={\sigma }_1^2/\tau +1/3·{\sigma }_2^2·\tau ---\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)={\sigma }_{i,1}\left(W_{i,1}\left(t_{k+1}\right)-W_{i,1}\left(t_k\right)\right)+{\sigma }_{i,2}{\int }_{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={\sigma }_{i,1}^2·T+1/3·{\sigma }_{i,2}^2·T^3---\left(5\right)$由式(5)和式(2)，得到：$Q_i={\sigma }_{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}=\phi ·{\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)=\left[1-K_{11}\left(k\right)-K_{12}\left(k\right)...-K_{1\left(N-1\right)}\left(k\right)\right]·\left[x_1\left(k\right)-{\hat{x}}_1\left(k-1\right)\right]\\ +K_{11}\left(k\right)·\left[x_2\left(k\right)-{\hat{x}}_2\left(k-1\right)\right]+K_{12}\left(k\right)·\left[x_3\left(k\right)-{\hat{x}}_3\left(k-1\right)\right]+...\\ +K_{1\left(N-1\right)}\left(k\right)·\left[x_N\left(k\right)-{\hat{x}}_N\left(k-1\right)\right]-\underset{i=1}{\overset{N-1}{\Sigma }}{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}{\Sigma }}{\omega }_{1,j}·\left(x_j\right(k)-{x_j}^{\prime }(k\left)\right)---\left(20\right)$其中，时间尺度的权重表示为：$\begin{array}{ccc}{\omega }_{1,1}=1-\underset{j=1}{\overset{N-1}{\Sigma }}{K}_{1j}\left(k\right)& {\omega }_{1,i}=K_{1\left(j-1\right)}\left(k\right)& (j\ne 1)\end{array}---\left(21\right)$时间尺度的预测值表示为：${x_j}^{\prime }\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}{\Sigma }}{\omega }_{i,j}·[x_j\left(k\right)-{\hat{x}}_j(k-1)]---\left(23\right)$由式(18)，当i≠1时，时间尺度的权重表示为：$\begin{array}{cccc}{\omega }_{i,1}=-\underset{j=1}{\overset{N-1}{\Sigma }}{K}_{ij}\left(k\right)& {\omega }_{i,i}=1+K_{i(i-1)}\left(k\right)& {\omega }_{i,j}=K_{i(j-1)}\left(k\right)& (j\ne 1,j\ne i)\end{array}---\left(24\right)$而时间尺度的预测值的表达式与式(22)相同。