CN104199058A - 基于Kalman滤波器实时预测值调整时间尺度算法 - Google Patents

Abstract

本发明属于时间频率、信号处理领域，提供一种基于Kalman滤波器的实时预测值调整时间尺度算法。时间尺度基本方程表明：时间尺度算法的核心是通过(N-1)组观测量(钟差)计算得到N台钟的权重和预测值。传统方法主要关注如何设置权重来提高时间尺度的稳定度。本发明提出通过实时调整预测值来建立时间尺度，并提高时间尺度的稳定度，同时保证时间尺度的连续性。本发明的技术方案包括：1)测量得到(N-1)组观测钟差；2)运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计；3)对时差序列进行重构；4)定义预测值；5)进行实时预测值调整，建立时间尺度。相比传统算法，本发明可以建立一个稳定度更高，同时也是连续的、实时的、可预测的时间尺度。

Description

基于Kalman滤波器实时预测值调整时间尺度算法

技术领域

本发明属于时间频率、信号处理领域，具体的说涉及一种基于Kalman滤波器的实时预测值调整时间尺度算法。

背景技术

在卫星导航系统的建设中，系统时间的建立和维持处于基础和核心的地位。卫星导航系统本质上是一个时间同步系统，卫星钟差的解算和预报要以高稳定度的时间基准为基础，而时间尺度的性能既取决于各台原子钟的性能，也取决于所设计的时间尺度算法。传统的时间尺度算法主要包括：BIPM采用的ALGOS算法、NIST采用的AT1算法、NIST在1980年代使用的标准Jones-Tyron型Kalman滤波器算法。

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

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

其中，TA(t)是建立的时间尺度，通常是纸面上的时间尺度；h_i(t)代表第i台(i＝1,2，…，N)原子钟的钟面读数，h_i′(t)代表第i台原子钟钟面读数的预测值；x_i(t)代表第i台原子钟的时差，x_i′(t)代表第i台原子钟时差的预测值；ω_i为第i台原子钟的权重。

时间尺度基本方程表明：时间尺度算法的核心是通过(N-1)组观测量(钟差)计算得到N台钟的权重和预测值。

大量文献都把研究重点放在如何设置权重来提高时间尺度的稳定度；对于预测值，考虑的仅仅是抑制时间尺度在两个相邻的时间段内由于权重改变引起的相位和频率跳变，以此来保证时间尺度的连续性。已经有文献证明：对于一个钟组，如果要通过设置合理的权重来提高时间尺度的稳定性，潜力是有限的，即如果不考虑预测值的情况下，加权平均算法用于提高组合钟时间尺度的稳定度是有一个理论极限的(卫国，原子钟时间尺度的稳定度分析，中国科学A辑，1992年1月，pp：80-86)。

发明内容

本发明的目的是针对现有时间尺度算法存在的缺点，提出一种基于Kalman滤波器的实时预测值调整时间尺度算法。该方法主要关注如何实时调整预测值来提高时间尺度的稳定度。

本发明的原理是：首先运用Kalman滤波器对观测钟差进行状态估计，Kalman滤波器每递推一次获得了一个新的预测值，通过实时调整预测值来建立时间尺度。通常的原子钟例如氢钟或铯钟主要包含频率白噪声(WFM)和频率随机游走噪声(RWFM)。通过本发明获取的预测值中主要只含有RWFM，所以建立的时间尺度中也主要只含有RWFM，所以该时间尺度的稳定度很高。

本发明采用的技术方案是：

一种基于Kalman滤波器的实时预测值调整时间尺度算法。算法主要包括以下步骤：

S1.对于一个含有N台钟的钟组，测量得到(N-1)组观测钟差；

S2.根据步骤S1得到的(N-1)组观测钟差，运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计，获取状态变量的估计值；

S3.根据步骤S2得到的状态变量的估计值，对时差序列进行重构；

S4.根据步骤S3得到的重构时差序列，定义预测值；

S5.根据步骤S4定义得到的预测值，进行实时预测值调整，建立时间尺度TA。

本发明的有益效果是：

1.观测钟差中含有WFM和RWFM，而使用本发明建立得到的时间尺度TA中不含WFM，其主要噪声为RWFM，所以稳定度明显提高；

2.该时间尺度在时差和频差上都保证连续；

3.该时间尺度是一个实时的、可预测的时间尺度。

附图说明

图1图解说明X1和

图2图解说明X₂和X₁的导数dX₁/dt(即频差)；

图3图解说明X₂和

图4图解说明观测钟差和时间尺度TA的时差；

图5图解说明观测钟差和和时间尺度TA的Allan偏差。

具体实施方式

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

本发明首先获取(N-1)组钟差；然后使用(N-1)个独立的Kalman滤波器对(N-1)组钟差进行状态估计；根据得到的估计量，对时差进行重构；由该重构时差序列，定义预测值；最后，实时调整预测值，建立时间尺度。

本发明所述方法包括以下步骤：

S1.对于一个含有N台钟的钟组，测量得到(N-1)组观测钟差；

S2.根据步骤S1得到的(N-1)组观测钟差，运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计，获取状态变量的估计值；

S2.1为了更清晰的描述步骤S2，首先描述根据步骤S1得到的1组观测钟差，运用1个相互独立的Kalman滤波器对这1组观测钟差分别进行状态估计，获取状态变量的估计值的情况；

对于任意一组观测钟差，可以用下面的模型来描述：

$\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(t\right)\mathrm{dt}\\ X_2\left(t\right)=y_0+d\·t+{\mathrm{\σ}}_2\·W_2\left(t\right)\end{array}\right.---\left(2\right)$

式(2)中，X₁(t)、X₂(t)分别为该观测钟差的两个状态变量。W₁(t)和W₂(t)分别代表两个独立的维纳过程(W(t))，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布。σ₁和σ₂分别是这两个维纳过程的扩散系数(diffusion coefficients)，用于表明噪声的强度。x₀、y₀、d分别代表初始时差、初始频差和频漂。式(2)表明：这两个维纳过程W₁(t)和W₂(t)在状态变量X₁(t)上分别表现为相位随机游走(对应于WFM)和积分相位随机游走(对应于RWFM)；维纳过程W₂(t)在状态变量X₂(t)上表现为频率随机游走(对应于RWFM)。

将式(2)按照间隔T进行离散化，得到：

$\left\{\begin{array}{c}{X}_1\left(t_{k+1}\right)=X_1\left(t_k\right)+X_2\left(t_k\right)\·T+d\·\frac{T^2}{2}+{\mathrm{\σ}}_1\·[W_1\left(t_{k+1}\right)-W_1\left(t_k\right)]+{\mathrm{\σ}}_2\·{\∫}_{t_k}^{t_{k+1}}{W}_2\left(s\right)\mathrm{ds}\\ X_2\left(t_{k+1}\right)=X_2\left(t_k\right)+d\·T+{\mathrm{\σ}}_2\·[W_2\left(t_{k+1}\right)-W_2\left(t_k\right)]\end{array}\right.---\left(3\right)$

式(3)就是观测钟差的状态方程。状态变量中的噪声分量J_k为：

$J_k=\left[\begin{array}{c}{\mathrm{\σ}}_1\·[W_1\left(t_{k+1}\right)-W_1\left(t_k\right)]+{\mathrm{\σ}}_2\·{\∫}_{t_k}^{t_{k+1}}{W}_2\left(s\right)\mathrm{ds}\\ {\mathrm{\σ}}_2\·[W_2\left(t_{k+1}\right)-W_2\left(t_k\right)]\end{array}\right]---\left(4\right)$

于是，状态变量噪声协方差Q可以表示为：

$Q=E[J_k\·{J_k}^T]=\left[\begin{array}{cc}{{\mathrm{\σ}}_1}^2\·T+1/3\·{{\mathrm{\σ}}_2}^2\·T^3& 1/2\·{{\mathrm{\σ}}_2}^2\·T^2\\ 1/2\·{{\mathrm{\σ}}_2}^2\·T^2& {{\mathrm{\σ}}_2}^2\·T\end{array}\right]---\left(5\right)$

符号E代表求均值，J_k ^T表示矩阵的转置。

观测钟差的离散化观测方程表述如下：

Z(t_k)＝X₁(t_k)+σ·ξ(t_k)              (6)

在式(6)中，Z(t_k)为观测量；σ·ξ(t_k)为观测噪声，其中ξ是标准高斯白噪声，被认为是相位白噪声(WPM)，σ用于表明观测噪声的强度。观测噪声协方差表示为：

R＝σ²                   (7)

使用Kalman滤波器估计原子钟的状态变量，可以通过下面5个方程进行描述。该5个方程中的变量的定义是相关领域内共知的，因此不再描述这些变量的定义。

1、状态预测

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

2、预测误差更新

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

3、Kalman增益更新

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

4、状态更新

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

5、估计误差更新

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

Kalman滤波器作为最小均方意义下的最优估计算法，可以有效地估计钟差的状态变量，最终得到这状态变量X₁和X₂的估计值，分别记为：和

S2.2在S2.1的基础上，将1组观测钟差扩展到(N-1)组观测钟差，根据步骤S1得到的(N-1)组观测钟差，运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计，获取状态变量的估计值；

为了对这(N-1)组钟差的状态变量进行区分，这些观测钟差的状态变量分别记为X₁ ^(1,n)、X₂ ^(1,n)，(n＝2，…，N)；其中上标(1,n)表示该组观测钟差是第1台钟与第n(2≤n≤N)台钟之间的观测钟差。对于每一组钟差，使用Kalman滤波器对进行状态估计，由式(8)-(12)所示的5个方程来描述。最终得到每组钟差的状态变量X₁ ^(1,n)和X₂ ^(1,n)的估计值，分别记为：和

S2.3验证Kalman滤波器用于状态估计的有效性

仿真生成两台钟组成钟组，来验证算法性能。因为此时只有钟组中两台钟，所以只能获取1组观测钟差，在这种情况下，把状态变量X₁ ^(1,n)、X₂ ^(1,n)分别简记为X₁、X₂；把状态变量估计值简记为使用Kalman滤波器对钟差进行状态估计，得到和

图1为X₁和从图1中可以看出：估计得到的和X₁的曲线基本重合，表明Kalman滤波器对于X₁的估计误差很小。

图2为X₂和X₁的导数dX₁/dt(即频差)。由式(2)可知，X₁的导数dX₁/dt和X₂分别表示为：

dX₁(t)/dt＝y₀+d·t+σ₁·dW₁(t)/dt+σ₂W₂(t)     (13)

和

X₂(t)＝y₀+d·t+σ₂W₂(t)                  (14)

对比式(13)和式(14)，状态变量X₂可以认为是不含WFM的频差。从图2可以看到：X₂的噪声分量明显小于频差dX₁/dt。

图3为X₂和图3表明：Kalman滤波器可以有效估计X₂，估计得到的中主要只含有RWFM。

S3.根据步骤S2得到的状态变量的估计值，对这(N-1)组钟差序列进行重构；

令：

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

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

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

显然 ${{\tilde{X}}_1}^{(1,n)}(n=2,...,N)$ 是通过估计量 ${{\hat{X}}_2}^{(1,n)}(n=2,...,N)$ 进行重构得到的钟差序列。

通过同样的方法可以得到(N-1)个重构的时差序列

S4.根据步骤S3得到的重构钟差序列，定义预测值；

把重构得到的钟差定义为每台钟的时差与时差预测值之间的偏差。在这种情况下，第1台钟的预测值表示为：

${x_1}^{\′}\left(t_k\right)=x_1\left(t_k\right)-{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(17\right)$

其中n可以是2，3，…，N中任意值。

第n(2≤n≤N)台钟的预测值表示为：

${x_n}^{\′}\left(t_k\right)=x_n\left(t_k\right)-{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(18\right)$

S5.根据步骤S4定义得到的预测值，进行实时预测值调整，建立时间尺度TA。

将预测值的表达式代入时间尺度基本方程，当钟组中含有2台钟时，得到时间尺度TA，即：

$\begin{array}{c}\mathrm{TA}\left(t_k\right)={\mathrm{\ω}}_1\·[x_1\left(t_k\right)-{x_1}^{\′}\left(t_k\right)]+{\mathrm{\ω}}_2\·[x_2\left(t_k\right)-{x_2}^{\′}\left(t_k\right)]\\ ={\tilde{X}}_1\left(t_k\right)={\tilde{X}}_1\left(t_{k-1}\right)+{\hat{X}}_2\left(t_{k-1}\right)\·T=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{\hat{X}}_2\left(t_m\right)\·T\end{array}---\left(19\right)$

当钟组中含有N(N≥3)台钟时，得到时间尺度TA，即：

$\mathrm{TA}\left(t_k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\·[x_i\left(t_k\right)-{x_i}^{\′}\left(t_k\right)]={\mathrm{\ω}}_1\·{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)+\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_n\·{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(20\right)$

至此，由S1～S5共5个步骤，建立了时间尺度TA。分析该5个步骤，可以得出如下结论：

1)本发明的核心思想是在Kalman滤波器每一次递推的过程中，调整一次预测值，通过每次实时调整预测值来建立时间尺度TA；

2)当钟组中只含有2台钟时，式(19)表明：权重对时间尺度TA没有任何影响；该TA最终可以表示为重构得到的时差序列(16)。所以，当使用本算法建立时间尺度时，可以抛开时间尺度的定义(即时间尺度基本方程(1))，直接通过估计量重构时差，得到时间尺度TA；

3)当钟组中含有N台钟时，可以抛开时间尺度标准方程，直接对这N-1个重构时差序列进行加权平均，最终可以得到如式(21)所示的时间尺度TA(tk)。

$\mathrm{TA}\left(t_k\right)=\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{1,n}\·{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)=\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_{1,n}\·[\underset{n=0}{\overset{k-1}{\mathrm{\Σ}}}{{\hat{X}}_2}^{(1,n)}\left(t_k\right)\·T]---\left(21\right)$

对比分析式(20)和式(21)，这两个表达式是等价的。

4)本算法建立了一个连续的、实时的、可预测的时间尺度。

下面从理论上分析该时间尺度的稳定度。

当钟组中含有2台钟时，式(19)可以改写为如下的形式：

$\begin{array}{c}\mathrm{TA}\left(t_k\right)={\tilde{X}}_1\left(t_k\right)=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{\hat{X}}_2\left(t_m\right)\·T=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}[X_2\left(t_m\right)+\mathrm{\ϵ}\left(t_m\right)]\·T\\ =y_0\·t_k+0.5\·d\·{t_k}^2+{\mathrm{\σ}}_2\·\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{W}_2\left(t_m\right)\·T+\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}\mathrm{\ϵ}\left(t_m\right)\·T\\ \≈y_0\·t_k+0.5\·d\·{t_k}^2+{\mathrm{\σ}}_2\·{\∫}_0^{t_k}{W}_2\left(t\right)\mathrm{dt}+\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}\mathrm{\ϵ}\left(t_m\right)\·T\end{array}---\left(22\right)$

其中，ε(t_m)为Kalman滤波器状态估计的误差。

式(22)表明，在不考虑误差项的情况下，该时间尺度TA只含有RWFM，所以稳定度很高。

当钟组中含有N台钟时，分析式(22)，时间尺度TA(t_k)是对主要含有RWFM的重构时差进行加权平均，因此同样主要含有RWFM，稳定度很高。随着钟组中原子钟的数量的增加，通过设计合理的权重，该时间尺度的长期稳定度也将得到提高。

图4描述了当钟组中只含有2台钟时，观测钟差和时间尺度TA的时差。图5描述了观测钟差和时间尺度TA的Allan偏差。

分析图4和图5可知：观测钟差的噪声主要表现为WFM和RWFM，而时间尺度TA中不含WFM，其主要噪声为RWFM；该TA在时差和频差上都保证连续。

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

Claims (1)

1.一种基于Kalman滤波器实时预测值调整时间尺度算法，其特征在于该算法主要包括以下步骤：

S1.对于一个含有N台钟的钟组，测量得到(N-1)组观测钟差；

S2.根据步骤S1得到的(N-1)组观测钟差，运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计，获取状态变量的估计值：

S2.1对于步骤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(t\right)\mathrm{dt}\\ X_2\left(t\right)=y_0+d\·t+{\mathrm{\σ}}_2\·W_2\left(t\right)\end{array}\right.---\left(2\right)$

式(2)中，X₁(t)、X₂(t)分别为该观测钟差的两个状态变量，W₁(t)和W₂(t)分别代表两个独立的维纳过程(W(t))，并且有W(t)～N(0,t)，即每个维纳过程服从均值为0，方差为时间t的正态分布，σ₁和σ₂分别是这两个维纳过程的扩散系数，用于表明噪声的强度，x₀、y₀、d分别代表初始时差、初始频差和频漂，将式错误！未找到引用源。按照间隔T进行离散化，得到：

$\left\{\begin{array}{c}{X}_1\left(t_{k+1}\right)=X_1\left(t_k\right)+X_2\left(t_k\right)\·T+d\·\frac{T^2}{2}+{\mathrm{\σ}}_1\·[W_1\left(t_{k+1}\right)-W_1\left(t_k\right)]+{\mathrm{\σ}}_2\·{\∫}_{t_k}^{t_{k+1}}{W}_2\left(s\right)\mathrm{ds}\\ X_2\left(t_{k+1}\right)=X_2\left(t_k\right)+d\·T+{\mathrm{\σ}}_2\·[W_2\left(t_{k+1}\right)-W_2\left(t_k\right)]\end{array}\right.---\left(1\right)$

式(1)就是观测钟差的状态方程，状态变量中的噪声分量J_k为：

$J_k=\left[\begin{array}{c}{\mathrm{\σ}}_1\·[W_1\left(t_{k+1}\right)-W_1\left(t_k\right)]+{\mathrm{\σ}}_2\·{\∫}_{t_k}^{t_{k+1}}{W}_2\left(s\right)\mathrm{da}\\ {\mathrm{\σ}}_2\·[W_2\left(t_{k+1}\right)-W_2\left(t_k\right)]\end{array}\right]---\left(2\right)$

状态变量噪声协方差Q可以表示为：

$Q=E[J_k\·{J_k}^T]=\left[\begin{array}{cc}{{\mathrm{\σ}}_1}^2\·T+1/3\·{{\mathrm{\σ}}_2}^2\·T^3& 1/2\·{{\mathrm{\σ}}_2}^2\·T^2\\ 1/2\·{{\mathrm{\σ}}_2}^2\·T^2& {{\mathrm{\σ}}_2}^2\·T\end{array}\right]---\left(3\right)$

符号E代表求均值，J_k ^T表示矩阵的转置；

观测钟差的离散化观测方程表述如下：

Z(t_k)＝X₁(t_k)+σ·ξ(t_k)     (4)

在式(4)中，Z(t_k)为观测量；σ·ξ(t_k)为观测噪声，其中ξ是标准高斯白噪声，被认为是相位白噪声(WPM)，σ用于表明观测噪声的强度，观测噪声协方差表示为：

R＝σ²           (5)

使用Kalman滤波器估计原子钟的状态变量，通过下面5个方程进行描述，

1、状态预测

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

2、预测误差更新

$P_{k,k-1}={\mathrm{\φP}}_{k-1,k-1}{\mathrm{\φ}}^T+Q---\left(7\right)$

3、Kalman增益更新

$K_k=P_{k,k-1}\·H^T{(H\·P_{k,k-1}\·H^T+R)}^{-1}---\left(8\right)$

4、状态更新 ${\hat{X}}_{k,k}={\hat{X}}_{k,k-1}+K_k\·(Z_k-H\·{\hat{X}}_{k,k-1})---\left(9\right)$

5、估计误差更新

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

Kalman滤波器作为最小均方意义下的最优估计算法，可以有效地估计钟差的状态变量，最终得到这状态变量X₁和X₂的估计值，分别记为：和

S2.2在S2.1的基础上，将1组观测钟差扩展到(N-1)组观测钟差，根据步骤S1得到的(N-1)组观测钟差，运用(N-1)个相互独立的Kalman滤波器对这(N-1)组观测钟差分别进行状态估计，获取状态变量的估计值；

为了对这(N-1)组钟差的状态变量进行区分，这些观测钟差的状态变量分别记为X₁ ^(1,n)、X₂ ^(1,n)，(n＝2，…，N)；其中上标(1,n)表示该组观测钟差是第1台钟与第n(2≤n≤N)台钟之间的观测钟差；对于每一组钟差，使用Kalman滤波器对进行状态估计，由式(6)-(10)所示的5个方程来描述；最终得到每组钟差的状态变量X₁ ^(1,n)和X₂ ^(1,n)的估计值，分别记为：和

S3.根据步骤S2得到的状态变量的估计值，对时差序列进行重构：

令：

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

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

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

显然(n＝2，…，N)是通过估计量(n＝2，…，N)进行重构得到的钟差序列；

通过同样的方法可以得到(N-1)个重构的时差序列(n＝2，…，N)；

S4.根据步骤S3得到的重构时差序列，定义预测值：

把重构得到的钟差定义为每台钟的时差与时差预测值之间的偏差，在这种情况下，第1台钟的预测值表示为：

${x_1}^{\′}\left(t_k\right)=x_1\left(t_k\right)-{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(13\right)$

其中n可以是2，3，…，N中任意值；

第n(2≤n≤N)台钟的预测值表示为：

${x_n}^{\′}\left(t_k\right)=x_n\left(t_k\right)-{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(14\right)$

S5.根据步骤S4定义得到的预测值，进行实时预测值调整，建立时间尺度TA：

将预测值的表达式代入时间尺度基本方程，当钟组中含有2台钟时，得到时间尺度TA，即：

$\begin{array}{c}\mathrm{TA}\left(t_k\right)={\mathrm{\ω}}_1\·[x_1\left(t_k\right)-{x_1}^{\′}\left(t_k\right)]+{\mathrm{\ω}}_2\·[x_2\left(t_k\right)-{x_2}^{\′}\left(t_k\right)]\\ ={\tilde{X}}_1\left(t_k\right)={\tilde{X}}_1\left(t_{k-1}\right)+{\hat{X}}_2\left(t_{k-1}\right)\·T=\underset{m=0}{\overset{k-1}{\mathrm{\Σ}}}{\hat{X}}_2\left(t_m\right)\·T\end{array}---\left(15\right)$

当钟组中含有N(N≥3)台钟时，得到时间尺度TA，即：

$\mathrm{TA}\left(t_k\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_i\·[x_i\left(t_k\right)-{x_i}^{\′}\left(t_k\right)]={\mathrm{\ω}}_1\·{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)+\underset{n=2}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ω}}_n\·{{\tilde{X}}_1}^{(1,n)}\left(t_k\right)---\left(16\right).$