A new approach to UTC calculation by means of the Kalman filter

Details concerning the clock frequency prediction indicated by term $h_{i}^{\prime}(t)$ and used in the calculation of UTC can be found in [5]. The prediction term is expressed in quadratic form to describe the frequency drift of the H-masers or the ageing of the caesium clocks. It includes three parameters: the phase, the frequency and the frequency drift:

• the estimation of the time correction relative to EAL of the clock H_i at the date t_k, added to avoid time jumps, is given by the last known difference between EAL and each clock
• the estimation of the frequency of clock H_i, relative to EAL, predicted for the calculation period, to avoid frequency jumps, from a theoretical point of view is the difference with respect to EAL between the first and last clock data points for the interval
  $$\begin{eqnarray}&&{{\hat{y}}_{i}}=\frac{{{x}_{i}}\left({{t}_{k+1}}\right)-{{x}_{i}}\left({{t}_{k}}\right)}{{{t}_{k+1}}-{{t}_{k}}}\end{eqnarray} \tag{ 3 }$$
• the estimation of the frequency drift of the clock H_i, relative to a frequency reference, predicted for the calculation period, under the hypothesis that this drift remains constant over a 1 month interval is given by using the frequency data ${{y}_{TT-{{h}_{i}}}}$ of the clock with respect to terrestrial time (TT) [9, 10]. The least-squares technique on 6 months of past data of ${{y}_{TT-{{h}_{i}}}}$ is used to evaluate the frequency drift for all kinds of clocks.

The weighting strategy is based on the principle that a good clock is a predictable clock [6, 11]. The main idea is to analyze the difference between the frequency of the clocks y_i and their prediction y_ip, obtained by using the relationship (3). In UTC calculation the upper limit of weight is set equal to $\text{EAL}{{\text{W}}_{\text{max}}}=4~{{\text{N}}^{-1}}$ to limit individual clock contributions. N is the number of the clocks used in the calculation. The weighting algorithm is a four-iteration process. The differences between the predicted and the real frequencies are evaluated for each 1 month interval over one year and these values are used to define the weight. Each iteration runs as follows:

• 1.  
  The values $\left[\text{EAL}-{{h}_{i}}\right]$ are found using a given set of relative weights. In the first iteration, the weights are those obtained in the previous computation interval after normalization. In the following iterations, the weights are obtained from the previous iteration;
• 2.  
  The absolute value of the difference between the frequency y_{i, j} as defined in (3) and the predicted frequency y_{ip, j} evaluated on the calculation interval j obtained by the relationship corresponds to:

  $$\begin{eqnarray}&&{{\epsilon}_{i,j}}\,=\,|\,{{y}_{i,j}}-{{y}_{ip,j}}|\end{eqnarray} \tag{ 4 }$$

  where the index i identifies the clock;
• 3.  
  One year of ${{\epsilon}_{i,j}}$ is considered to ensure long-term stability of EAL and UTC, j characterizes the calculation interval.
• 4.  
  A filter has been implemented to give a more predominant role to more recent measurements with respect to older ones, considering that new measurements have the most reliable statistics:

  $$\begin{eqnarray}&&\sigma _{i}^{2}=\frac{{\sum}_{j=1}^{M}\left(\frac{M+1-j}{M}\right)\epsilon _{i,j}^{2}}{{\sum}_{j=1}^{M}\left(\frac{M+1-j}{M}\right)}\end{eqnarray} \tag{ 5 }$$

  where i identifies the clock, j the calculation interval and M the number of available measurements.
• 5.  
  The relative weight of clock H_i is computed theoretically using a temporary value given by

  $$\begin{eqnarray}&&\text{EAL}{{\text{W}}_{i,\text{temp}}}=\frac{1/\sigma _{i}^{2}}{{\sum}_{i=1}^{N}1/\sigma _{i}^{2}}.\end{eqnarray} \tag{ 6 }$$

The new weight $\text{EAL}{{\text{W}}_{i}}$ of clock ${{\text{H}}_{i}}$ is equal to $\text{EAL}{{\text{W}}_{i,\text{temp}}}$ except in two cases:

• 1.  
  Clock H_i exceeds the limit value set for the limitation of weight as for the current algorithm.
• 2.  
  Clock H_i shows abnormal behaviour during the interval of computation so it cannot contribute. In this case the current value of the difference between the real frequency and the predicted frequency is checked. If the value is greater than 5 ns d⁻¹ the clock is temporarily excluded from the ensemble.

Once the necessary information has been obtained, the Kalman filter is constituted of two parts: an a priori estimation is made based solely on the knowledge we have of the system. Then, this estimation is updated to an a posteriori estimation using the additional information provided by the new set of measurement that are now available. Let us have a closer look at the procedure. For every time instant t_k, we proceed as follows:

• STEP 1: PredictionStarting from the given initial estimate $\mathbf{\hat{X}}\left({{t}_{0}}\right)$, for a set of time instants $\left\{{{t}_{k}}\right\}_{k=1}^{n}$ equally spaced at distance τ³, we compute
  $$\begin{eqnarray}&&\mathbf{\hat{X}}\left({{t}_{k}}|{{t}_{k-1}}\right)=\boldsymbol{}{{\mathbf{\Phi}}_{{{t}_{k}},{{t}_{k-1}}}}\mathbf{\hat{X}}\left({{t}_{k-1}}|{{t}_{k-1}}\right),\end{eqnarray} \tag{ 7 }$$
  where the notation means that the estimate is computed at time t_k using all the information up to (and including) the instant t_k−1. The covariance matrix for the estimation error is as well computed starting from $\boldsymbol{\Gamma}\left({{t}_{0}}\right)$ as
  $$\begin{eqnarray}&&\boldsymbol{\Gamma}\left({{t}_{k}}|{{t}_{k-1}}\right)=\boldsymbol{}{{\mathbf{\Phi}}_{{{t}_{k}},{{t}_{k-1}}}}\boldsymbol{\Gamma}\left({{t}_{k-1}}|{{t}_{k-1}}\right)\boldsymbol{\Phi}_{{{t}_{k}},{{t}_{k-1}}}^{T}+\mathbf{Q}\left({{t}_{k}},{{t}_{k-1}}\right),\end{eqnarray} \tag{ 8 }$$
  where Q is the covariance matrix of the noise vector $\mathbf{W}$.
• STEP 2: MeasureAt each time step t_k, the measurement vector can be expressed as $\mathbf{Z}\left({{t}_{k}}\right)=\mathbf{HX}\left({{t}_{k}}\right)+\mathbf{v}\left({{t}_{k}}\right)$, where $\mathbf{v}$ represents the measurement noise⁴ and $\mathbf{H}$ is the transformation matrix. We assume that the measurement noise $\mathbf{v}$ is a zero-mean Gaussian noise with covariance matrix $\mathbf{R}\left({{t}_{k}}\right)$.We then define a variable, called innovation, that computes the difference between the vector of true measurements $\mathbf{Z}$ and the estimated one.
  $$\begin{eqnarray}&&\boldsymbol{i}\left({{t}_{k}}\right)=\mathbf{Z}\left({{t}_{k}}\right)-\mathbf{H\hat{X}}\left({{t}_{k}}|{{t}_{k-1}}\right)\end{eqnarray} \tag{ 9 }$$
  whose covariance matrix is given by
  $$\begin{eqnarray}&&\mathbf{C}\left({{t}_{k}}\right)=\mathbf{H\Gamma}\left({{t}_{k}}|{{t}_{k-1}}\right){{\mathbf{H}}^{T}}+\mathbf{R}\left({{t}_{k}}\right).\end{eqnarray} \tag{ 10 }$$
  It is now possible to define the Kalman gain matrix as
  $$\begin{eqnarray}&&\mathbf{K}\left({{t}_{k}}\right)=\boldsymbol{\Gamma}\left({{t}_{k}}|{{t}_{k-1}}\right){{\mathbf{H}}^{T}}\mathbf{C}{{\left({{t}_{k}}\right)}^{-1}}.\end{eqnarray} \tag{ 11 }$$
• STEP 3: Estimation updateThe estimations are now updated using the innovation and the Kalman gain in the following way
  $$\begin{eqnarray}&&\mathbf{\hat{X}}\left({{t}_{k}}|{{t}_{k}}\right)=\mathbf{\hat{X}}\left({{t}_{k}}|{{t}_{k-1}}\right)+\mathbf{K}\left({{t}_{k}}\right)\mathbf{i}\left({{t}_{k}}\right)\end{eqnarray} \tag{ 12 }$$
  of covariance matrix
  $$\begin{eqnarray}&&\boldsymbol{\Gamma}\left({{t}_{k}}|{{t}_{k}}\right)=\left[\mathbf{I}-\mathbf{K}\left({{t}_{k}}\right)\mathbf{H}\right]\boldsymbol{\Gamma}\left({{t}_{k}}|{{t}_{k-1}}\right).\end{eqnarray} \tag{ 13 }$$
  Now the estimation is made using all the information up to the instant t_k.

In the case of an ensemble of N atomic clocks, the state vector will contain the information of phase, frequency and drift offsets for each clock in the ensemble. In this work we chose to order the state vector by first locating all the phase offsets, then all the frequencies and then the drifts, having

$$\begin{eqnarray}&&\mathbf{X}(t)={{\left[{{x}_{1}}(t)\;\cdots \;{{x}_{N}}(t)\;{{y}_{1}}(t)\;\cdots \;{{y}_{N}}(t)\;{{z}_{1}}(t)\;\cdots \;{{z}_{N}}(t)\right]}^{\prime}}.\end{eqnarray} \tag{ 14 }$$

The resulting $\boldsymbol{\Phi}$ matrix is a block matrix

$$\begin{eqnarray}\boldsymbol{\Phi}(\tau )=\left[\begin{array}{*{35}{l}} {{\mathbf{I}}_{N}} & \tau \centerdot {{\mathbf{I}}_{N}} & \frac{{{\tau}^{2}}}{2}\centerdot {{\mathbf{I}}_{N}} \\ \mathbf{O} & {{\mathbf{I}}_{N}} & \tau \centerdot {{\mathbf{I}}_{N}} \\ \mathbf{O} & \mathbf{O} & {{\mathbf{I}}_{N}} \end{array}\right]\end{eqnarray} \tag{ 15 }$$

where ${{\mathbf{I}}_{N}}$ is the $N\times N$ identity matrix, and $\mathbf{O}$ is a null matrix of the same dimension.

For each clock i, the noise vector ${{\mathbf{W}}^{i}}(t)$ is a three dimensional vector whose components are independent Brownian motion. Thus we have ${{\mathbf{W}}^{i}}(t)\sim N\left(\mathbf{0},{{\mathbf{Q}}^{i}}\right)$, where the noise covariance matrix for the ith clock is given by

$$\begin{eqnarray}{{\mathbf{Q}}^{i}}=\left[\begin{array}{*{35}{l}} \sigma _{1,i}^{2}\tau +\sigma _{2,i}^{2}\frac{{{\tau}^{3}}}{3}+\sigma _{3,i}^{2}\frac{{{\tau}^{5}}}{20} & \sigma _{2,i}^{2}\frac{{{\tau}^{2}}}{2}+\sigma _{3,i}^{2}\frac{{{\tau}^{4}}}{8} & \sigma _{3,i}^{2}\frac{{{\tau}^{3}}}{6} \\ \sigma _{2,i}^{2}\frac{{{\tau}^{2}}}{2}+\sigma _{3,i}^{2}\frac{{{\tau}^{4}}}{8} & \sigma _{2,i}^{2}\tau +\sigma _{3,i}^{2}\frac{{{\tau}^{3}}}{3} & \sigma _{3,i}^{2}\frac{{{\tau}^{2}}}{2} \\ \sigma _{3,i}^{2}\frac{{{\tau}^{3}}}{6} & \sigma _{3,i}^{2}\frac{{{\tau}^{2}}}{2} & \sigma _{3,i}^{3}\tau \end{array}\right]\,\\ \,\,\,\,\,\,\,\,=\left[\begin{array}{*{35}{l}} Q_{xx}^{i} & Q_{xy}^{i} & Q_{xz}^{i} \\ Q_{yx}^{i} & Q_{yy}^{i} & Q_{yz}^{i} \\ Q_{zx}^{i} & Q_{zy}^{i} & Q_{zz}^{i} \end{array}\right]\end{eqnarray} \tag{ 16 }$$

with ${{\sigma}_{j,i}}$, j  =  1, 2, 3, representing the diffusion coefficients of the three components of ${{\mathbf{W}}^{i}}(t)$⁵. We have the noise covariance matrix for the whole system is given by the $3N\times 3N$

$$\begin{eqnarray}\mathbf{Q}=\left[\begin{array}{*{35}{l}} \text{diag}\left({{Q}_{xx}}\right) & \text{diag}\left({{Q}_{xy}}\right) & \text{diag}\left({{Q}_{xz}}\right) \\ \text{diag}\left({{Q}_{yx}}\right) & \text{diag}\left({{Q}_{yy}}\right) & \text{diag}\left({{Q}_{yz}}\right) \\ \text{diag}\left({{Q}_{zx}}\right) & \text{diag}\left({{Q}_{zy}}\right) & \text{diag}\left({{Q}_{zz}}\right) \end{array}\right],\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}\text{diag}\left({{Q}_{ab}}\right)=\left[\begin{array}{*{35}{l}} Q_{ab}^{1} & 0 & \cdots & 0 \\ 0 & Q_{ab}^{2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & Q_{ab}^{N} \end{array}\right].\end{eqnarray}$$

If we assume that the reference clock of the ensemble is located in the first position, then the measurement matrix $\mathbf{H}$ has to transform $\mathbf{X}(t)$ from (14) into a column vector of differences $\mathbf{Z}(t)={{\left[{{x}_{2}}(t)-{{x}_{1}}(t),\cdots,{{x}_{N}}(t)-{{x}_{1}}(t)\right]}^{\prime}}$ of dimension N  −  1. $\mathbf{H}$ would then be the $(N-1)\times 3N$ matrix made of a column of  −1s, a $(N-1)\times (N-1)$ identity matrix, and zeros at the remaining entries.

Note that if we include measurement noise v(t) in the model, we assume it to be uncorrelated with all the components of the clock noise vector W(t).

One of the major problems we have to face when dealing with Kalman estimation is the uncontrolled growth of the estimation error. Let us recall that due to the ordering (14), the covariance matrix of the estimation error, $\rm{\hat \Gamma }(t)$, is a ($3N\times 3N$) block matrix where the first diagonal block (${{\rm{\hat \Gamma }}_{xx}}$) contains the estimation error on the clocks phase components, the second block (${{\rm{\hat \Gamma }}_{yy}}$) on the frequency and the third (${{\rm{\hat \Gamma }}_{zz}}$) on the drift components. The off-diagonal blocks contain the interaction terms. While the yz-submatrix of the error covariance matrix $\rm{\hat \Gamma }(t)$, relative to frequency and drift estimates, has been noted to be empirically well behaved [14], the xx-submatrix grows without limits. This is mainly due to the fact that the Kalman filter has to estimate the behaviour of N variables having as an input a vector of dimension N  −  1, the $\mathbf{Z}$ measurements. To overcome this problem, Greenhall proposed an approach, known as Covariance x-reduction, that can be proved to be transparent with respect to Kalman estimates of frequency and drift [14, 15]. This procedure is carried out in a different way for the two cases of noiseless or noisy measurements. We define

$$\begin{eqnarray}&&\tilde{x}(t)=x(t)-\hat{x}(t),\end{eqnarray} \tag{ 18 }$$

where $\hat{x}(t)$ represents the Kalman phase offset estimate and $\tilde{x}(t)$ is the correction applied to x(t) on the base of Kalman prediction. We define the ${{\tilde{x}}_{i}}$ centroid ${{\bar{x}}_{0}}$

$$\begin{eqnarray}&&{{\bar{x}}_{0}}=\underset{i=1}{\overset{N}{\sum}}\,\text{K}{{\text{W}}_{i}}{{\tilde{x}}_{i}}\quad \quad \text{and}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{{\tilde{x}}^{\prime}}=\tilde{x}-{{1}_{N}}{{\bar{x}}_{0}}=\left({{\mathbf{I}}_{N}}-{{\mathbf{1}}_{N}}\mathbf{K}{{\mathbf{W}}^{t}}\right)\tilde{x},\end{eqnarray} \tag{ 20 }$$

where ${{\mathbf{1}}_{N}}$ is a N-dimensional column vector of 1s, ${{\mathbf{I}}_{N}}$ the ($N\times N$) identity matrix and KW a vector of weighs to which we will refer as the Kalman weights. As (20) expresses the error of representation of $\tilde{x}$ by means of ${{\bar{x}}_{0}}$, we would like to determine the centroid in such a way that this error is minimized. This can be accomplished finding the weights $\text{K}{{\text{W}}_{i}}$, i  =  1,..., N, that minimize the quantity between brackets in (20). In order to obtain the required weights we solve the system

$$\begin{eqnarray}\left[\begin{array}{c} {{\widehat{\boldsymbol{\Gamma}}}_{xx}} & {{\mathbf{1}}_{N}} \\ \mathbf{1}_{N}^{t} & 0 \end{array}\right]\left[\begin{array}{*{35}{l}} \mathbf{K}{{\mathbf{W}}_{i}} \\ \theta \end{array}\right]=\left[\begin{array}{c} \mathbf{O} \\ 1 \end{array}\right].\end{eqnarray} \tag{ 21 }$$

If we define

$$\begin{eqnarray}{{\mathbf{S}}_{r}}=\left[\begin{array}{c} {{\mathbf{I}}_{N}}-{{\mathbf{1}}_{N}}\mathbf{K}{{\mathbf{W}}^{t}} & \mathbf{O} \\ \mathbf{O} & {{\mathbf{I}}_{y}} \end{array}\right],\end{eqnarray} \tag{ 22 }$$

the covariance x-reduction operation in the case of noisy measurements consists of the following transformation

$$\begin{eqnarray}&&{{\widehat{\boldsymbol{\Gamma}}}_{r}}={{\mathbf{S}}_{r}}\widehat{\boldsymbol{\Gamma}}\mathbf{S}_{r}^{t}.\end{eqnarray} \tag{ 23 }$$

For the derivation of (21) and further detail on methods and properties see [15].

We recall the definition of clock reading ${{h}_{i}}(t)={{h}_{0}}(t)-{{x}_{i}}(t)$, where h₀ denotes an ideal clock. Our goal would be to estimate h₀ by means of the measurements

$$\begin{eqnarray}&&{{x}_{ij}}(t)={{h}_{i}}(t)-{{h}_{j}}(t)={{x}_{j}}(t)-{{x}_{i}}(t).\quad i,j=1,...,N\end{eqnarray} \tag{ 24 }$$

Having at our disposal an estimate for the x_i, ${{\hat{x}}_{i}}$ (in our case the output of a Kalman filter), we define the corrected clocks as

$$\begin{eqnarray}&&{{\tilde{h}}_{i}}={{h}_{i}}+{{\hat{x}}_{i}}.\end{eqnarray} \tag{ 25 }$$

Using the definition of clock reading and (18) we have

$$\begin{eqnarray}&&{{\tilde{h}}_{i}}={{h}_{0}}-{{\tilde{x}}_{i}}.\end{eqnarray} \tag{ 26 }$$

In the case of noiseless measures all the ${{\tilde{x}}_{i}}$ are equal and therefore all the corrected clocks agree. In the presence of measurement noise, however, the corrected clocks do not agree, but they tend nevertheless to form a cluster. In the sight of (26), we choose to take the centroid of the cluster as the least square estimate of h₀ given the ${{\tilde{h}}_{i}}$. We define

$$\begin{eqnarray}&&{{\bar{h}}_{0}}(t)=\underset{i=1}{\overset{N}{\sum}}\,{{w}_{gi}}{{\tilde{h}}_{i}}(t)=\mathbf{w}_{g}^{t}\mathbf{\tilde{h}}.\end{eqnarray} \tag{ 27 }$$

Using (19) and once again the definition of clock reading we get

$$\begin{eqnarray}&&{{\bar{h}}_{0}}(t)={{h}_{0}}(t)-\underset{i=1}{\overset{N}{\sum}}\,{{w}_{gi}}{{\tilde{x}}_{i}}(t)={{h}_{0}}(t)-{{\bar{x}}_{0}}(t).\end{eqnarray} \tag{ 28 }$$

The quantity (27) defines a new time scale, the Kalman filter time scale (KF), where the vector ${{\mathbf{w}}_{g}}$ represents a generic weighting vector. If we define the representation error

$$\begin{eqnarray}&&\tilde{h}_{i}^{\prime}={{\tilde{h}}_{i}}-{{\bar{h}}_{0}},\quad \quad \text{for}~i=1,...,N,\end{eqnarray} \tag{ 29 }$$

that expresses how well each clock represents KF, we have the following result.

Theorem 1. The weight vector that, used to compute KF, minimizes the representation error is KW.

The reader can find the proof of this theorem in [15].

The first step we took in order to build a good simulation was to generate the clock data. This has been done according to the mathematical model presented in [13, 16]. The solution to the system of SDE can be expressed in the iterative form

$$\begin{eqnarray}&&\mathbf{X}\left({{t}_{k+1}}\right)=\boldsymbol{}{{\mathbf{\Phi}}_{{{t}_{k+1}},{{t}_{k}}}}\mathbf{X}\left({{t}_{k}}\right)+\mathbf{W}\left({{t}_{k}}\right).\end{eqnarray} \tag{ 30 }$$

where $\mathbf{X}\left(\mathbf{t}\right)=\left[{{X}_{1}}(t),{{X}_{2}}(t),{{X}_{3}}(t)\right]$. X₁(t) represents the phase offset of the modelled clock with respect to a reference clock, X₂(t) is part of the clock frequency deviation, while X₃(t) is part of the frequency drift. $\boldsymbol{}{{\mathbf{\Phi}}_{{{t}_{k+1}},{{t}_{k}}}}$ has the form (15) and $\mathbf{W}\left({{t}_{k}}\right)$ is the vector of clock noises whose covariance matrix, $\mathbf{Q}$, is the $3N\times 3N$ matrix defined in (17). The ${{\sigma}_{i}}$, the diffusion coefficients of the noise processes, express the amount of white noise on the frequency component, random walk noise on frequency and to a perturbation on the frequency drift.

Once generating the clock signals, we compute the differences (24) that will correspond to the measurements. All the differences will be taken with respect to a clock used as the reference. The measurement noise will then be inserted by means of the R covariance matrix introduced in section 3.1. The time scale itself is then computed applying (28).

We run different simulations using different parameter specifications for the clocks and measurement noises. At first we simulated a set of four clocks, two having characteristics resembling the ones of Caesium clocks (white frequency noise (WFN) at the level of about $3\times {{10}^{-14}}$ at 5 d and random walk frequency noise (RWFN) of $1\times {{10}^{-16}}$ at 5 d), and two modelling Maser clocks (WFN at the level of about $1\times {{10}^{-15}}$ at 5 d and RWFN of $1\times {{10}^{-15}}$ at 5 d). The resulting scale stability with respect to the clock's stability is shown in figure 1. We can observe that the scale is able to achieve a better stability than any clock of the ensemble both in the short and long term.

Zoom In Zoom Out Reset image size

Another simulation was run substituting the Maser-like clocks with Rubidium fountains (figure 2). Unlike the previous case, the KF Scale does not appear to have a better stability than the Rubidium fountains, although it does show a very good performance. This is because in presence of time transfer noise with atomic clocks characterized by very different stabilities (in the previous example caesium clocks and Rb fountains) the Kalman Filter solution does not seem very stable in cleaning the white phase noise. If we run the same simulation but without introducing measurement noise, KF is actually showing a better performance than any clock, including the Rb fountains. The amount of time transfer noise used in the simulations reported in figure 2 is equal to 15 ns. In conclusion when dealing with high precision devices, even a small amount of measurement noise spoils the received signal in such away that the stability of the corresponding time scale will be prevented from achieving stability results comparable to those of the devices themselves [17].

Zoom In Zoom Out Reset image size

Unobservability. The first thing to be underlined when passing from simulation to real data is that some of the quantities used to define the KF Scale are not observable outside the simulation. In particular, one does not have access to the single clock readings, but just to the differences between the clock data and the data from a reference (generally a better performing clock or a time scale). As a consequence, we will not be able to compute the KF Scale itself, but just the differences between the scale and the clocks of the ensemble. Therefore, instead of (27) we will compute

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{{\bar{h}}}_{0}}(t)-{{h}_{j}}(t) & =\underset{i=1}{\overset{N}{\sum}}\,{{w}_{gi}}\left({{{\tilde{h}}}_{i}}(t)-{{h}_{j}}(t)\right) \\ {} & =\underset{i=1}{\overset{N}{\sum}}\,{{w}_{gi}}\left[{{x}_{ij}}(t)+{{{\hat{x}}}_{i}}(t)\right]. \end{array}\end{eqnarray} \tag{ 31 }$$

Measurement noise correlation. A second remark is due regarding the covariance matrix of the measurement error R. Although formally defined as a full matrix, in literature, and particularly in simulation, this matrix is generally presented as diagonal. This implies the assumption that the noise components affecting two different measures are uncorrelated. In the case of UTC this condition does not hold. In fact, the uncertainties of the link between laboratories are not computed directly, but the system relies on a pivot laboratory, the PTB (Physikalisch-Technische Bundesanstalt), Germany. The link uncertainty between each laboratory and the PTB is estimated and then all the other inter-laboratory uncertainties are computed by triangulation. This procedure originates a correlation between measurement errors, and therefore the R matrix can no longer be considered diagonal.

4.2.1. A Kalman based UTC prototype.

In order to test the potentiality of the presented Kalman filter Scale in the UTC framework, we built a prototype of the scale. In order to deal with one problem at a time, we put ourselves in a favourable condition selecting a subset of the UTC participating clocks such that:

• every clock was present for the whole considered time interval (≃3 years from 2012 to 2014);
• the clocks were not affected by any kind of phase or frequency jump.

The resulting ensemble consisted of 139 clocks from 19 different laboratories. In this section we will present the computation of three different time scales:

• $\text{EA}{{\text{L}}_{139}}$: EAL computed on the restricted sample of clocks, with maximum weight set to 2.5 N⁻¹, suitable in the case of a smaller sample;
• $\text{K}{{\text{F}}_{\text{KW}}}$: the Kalman filter time scale with Kalman weights;
• $\text{K}{{\text{F}}_{\text{EALW}}}$: the Kalman filter time scale with the weights resulting from the procedure described in section 2.2 applied to the 139 clocks sample.

To assess the quality of $\text{EA}{{\text{L}}_{139}}$, the reduced EAL has been compared to ensemble of PFS used to steer the frequency of EAL. In figure 3 the frequency of EAL calculated with about 450 (red line) and EAL with 139 (blue line) atomic clocks are compared to the frequency of PFS. In figure 4 the Allan variance of the previous data is reported. From this result is clear that $\text{EA}{{\text{L}}_{139}}$ is less stable that the current time scale but the quality is sufficiently good to be used. We should observe that the parameters (frequency of the clocks, frequency drift etc) are optimized for having a good long term time scale in the case of 450 atomic clocks. Different choices should be done in the case of about 150 atomic clocks. However the results are good enough for the scope of our first analysis. The input data for our computation was therefore of the form $\text{EA}{{\text{L}}_{139}}-{{h}_{i}}$ with i  =  1,..., 139. From that, we estimated the clock parameters (initial states and noise coefficients) attributing the whole noise signature to the h_i. All these clock parameters, obtained by the calculation of $\text{EA}{{\text{L}}_{139}}$, are used to initialize the Kalman filter algorithm. This is the main idea of the paper, a mixed algorithm taking advantage of both algorithms. The noise coefficients are calculated by comparing the clocks to TT(BIPM) and by performing the Hadamard variance. Figure 5 shows the data used for the statistical analysis of the clock with respect to $\text{EA}{{\text{L}}_{139}}$ (phase offset of USNO(Rb)-$\text{EA}{{\text{L}}_{139}}$) and TT(BIPM) (normalized frequency and Hadamard variance of USNO(Rb)-TT(BIPM)). The phase data with respect to $\text{EA}{{\text{L}}_{139}}-{{h}_{i}}$ are used for the phase estimation, TT(BIPM) is the best frequency reference used for the estimation of the frequency drift and for the noise estimation. In the plot of the Hadamard variance are reported the typical behaviour of the WFN (blue line) and RWFN (green line) used for the noise estimation. In the case of data of Rubidium fountain given by the USNO (United States Naval Observatory) [18], as the extremely high stability (both long and short term) the three cornered hat is used in order to evaluate the order of magnitude of the noise. Considering its excellent metrological properties, the comparison between the stability of $\text{EAL}-{{h}_{i}}$ and $\text{KF}-{{h}_{i}}$ when h_i is a fountain is of particular interest. In this case, having the clock with a better stability than the time scales itself, the plot actually shows the comparison between the stabilities of the atomic time scales. The noise estimation is the critical aspect when a time scale is built by means the Kalman filter. The solution can be completely erroneous only due to bad estimations. This aspect is of particular relevance in the case of Rubidium fountains due to their excellent performance. An adaptive algorithm should be studied that automatically adjusts the parameters as required. Several tests have already been carried out in this direction by Greenhall and Davis [19].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before completing the time scales comparison, there is another aspect on which we need to focus our attention: the clock weights. The x-reduction procedure in the case of measurement noise delivers, together with Kalman's estimates, a weight $\text{K}{{\text{W}}_{i}}$ for each laboratory, based solely on the clock stability. However this kind of weight does not suit our purposes as it does not introduce a maximum for each clock weight (thus affecting the reliability and the traceability of the scale). In this case the weights are completely dominated by the Rubidium fountains one, as their stability is much better of those of the other clocks. Therefore we use the EAL weights (EALW) for our computation, based, and previously described, on the principle that a good clock is a predictable clock. Figure 6 shows the stability comparison between $\text{EA}{{\text{L}}_{139}}$, the $\text{K}{{\text{F}}_{\text{KW}}}$ and the $\text{K}{{\text{F}}_{\text{EALW}}}$, all with respect to a Rubidium fountain. We can observe that the KF Scales show an improved stability with respect to EAL, particularly in the case of $\text{K}{{\text{F}}_{\text{KW}}}$. This is not surprising as not having a maximum weight constraint implies that the scale can rely only on the better performing clocks. This fact, even it allows a better stability, makes the scale much less robust to possible problems affecting these high performance clocks.

Zoom In Zoom Out Reset image size