# Research on Distributed Autonomous Timekeeping Algorithm for Low-Earth-Orbit Constellation
## 3 Methodology
### Parameter Estimation of Filter Noise Matrix  
When estimating the clock error using the Kalman filter, the observation noise is typically modeled as white noise, while the internal noise of the atomic clock is treated as colored noise. The autocorrelation function of white noise is an impulse function, indicating no correlation between epochs. In contrast, colored noise shows correlations between epochs, which can be linearly estimated based on the space formed by historical data and described by an \(AR(M)\) model [24] as follows:  
\[\begin{array}{l}x(k)+a_{1}x(k-1)+\ldots+a_{M}x(k-M)=\varepsilon(k)\\ H(z)=\frac{1}{1+a_{1}z^{-1}+\ldots+a_{M}z^{-M}}\end{array} \tag{5}\]  
where \(M\) represents the order of the filter, \(a_{1},\ldots,a_{M}\) represent the filter coefficients, \(x(k)\) is the observed data at epoch k, \(\varepsilon(k)\) is the white noise sequence at epoch k, and \(H(z)\) is the frequency domain expression of the filter.  
Therefore, an M-order prediction error filter is designed based on the noise model to filter the clock error. The variance of the prediction error derived from the filtered data is used as the estimated value of the observation noise covariance matrix in the Kalman filter. The canonical equation for the M-order prediction error filter is as follows:  
\[\begin{bmatrix}r(0)&r(1)&\cdots&r(M)\\ r(1)&r(0)&\cdots&r(M-1)\\ \vdots&\vdots&\ddots&\vdots\\ r(M)&r(M-1)&\cdots&r(0)\end{bmatrix}\begin{bmatrix}1\\ a_{1}\\ \vdots\\ a_{M}\end{bmatrix}=\begin{bmatrix}\rho\\ 0\\ \vdots\\ 0\end{bmatrix} \tag{6}\]  
or  
\[R_{x}H=\varepsilon \tag{7}\]  
where \(R_{x}\) is the M-order autocorrelation matrix, \(r(m)\) is the autocorrelation value with a difference of \(m\) epochs, \(H\) is the filter vector, \(\rho\) is the minimum value of error energy, and \(\varepsilon\) is the error vector. The observed data are passed through the filter, which is expressed as follows:  
\[y=X^{T}H \tag{8}\]where \(X\) is the input data vector, \(y\) is the prediction error, and \(\sigma^{2}\) is the variance of the prediction error series. The estimated observation noise covariance matrix \(R\) is as follows:  
\[R=\sigma^{2} \tag{9}\]  
The discrete state equation of the clock can be characterized as follows [25]:  
\[\begin{bmatrix}x(k)\\ f(k)\end{bmatrix}=\begin{bmatrix}1&T_{s}\\ 0&1\end{bmatrix}\begin{bmatrix}x(k-1)\\ f(k-1)\end{bmatrix}+w_{k} \tag{10}\]  
where \(k\) is the number of epochs, \(T_{s}\) is the sampling interval between epochs, \(w_{k}\) is the process noise vector, \(x(k)\) is the clock error at epoch k, and \(f(k)\) is the frequency deviation. The covariance matrix of the process noise is as follows:  
\[Q=E(w_{k}{w_{k}}^{T})=\begin{bmatrix}S_{t}T_{s}+S_{f}\frac{T_{s}^{3}}{3}&S_{f} \frac{T_{s}^{2}}{2}\\ S_{f}\frac{T_{s}^{2}}{2}&S_{f}T_{s}\end{bmatrix} \tag{11}\]  
where \(S_{t}\) and \(S_{f}\) are the noise power spectral density of the clock error and frequency offset, respectively, which can be obtained by fitting the Allan variance. The Allan variance is calculated as follows:  
\[\sigma^{2}(\tau)=\frac{1}{(N-2)\tau^{2}}\sum_{i=1}^{N-2}\frac{(x(k+2)-2x(k+1)+ x(k))^{2}}{2} \tag{12}\]  
where \(k\) is the number of epochs, and \(\sigma^{2}(\tau)\) is the calculated Allan variance. The formula for estimating \(S_{t}\) and \(S_{f}\) by the least square method is as follows:  
\[\sigma^{2}(\tau)=S_{t}\tau^{-1}+S_{f}\tau/3 \tag{13}\]