## 2 Methodology  
The frequency accuracy of the atomic clock is defined as:  
\[\sigma=\frac{f-f_{0}}{f_{0}} \tag{1}\]  
where \(\sigma\) denotes the frequency accuracy, \(f\) the actual frequency, and \(f_{0}\) the nominal frequency. Equation (1) describes the absolute concept of the frequency accuracy, which is not derived directly in actual data processing. Generally, the frequency accuracy is calculated by the time difference method, that is, the phase difference between two consecutive epochs is obtained:  
\[y_{i}=\frac{x(t_{1}+v)-x(t_{1})}{v} \tag{2}\]  
where \(x\) is the clock (phase) data, \(v\) the data sampling interval, and \(y_{i}\) the frequency data at epoch \(t_{i}\). With the phase data in a long time series, the frequency accuracy can be obtained by calculating the mean value of the differenced data series. It should be noted here that Equation. (2) is also the formula of converting the phase data to frequency data. The frequency drift rate can be written as:  
\[D=\frac{\sum_{i=1}^{N}(y_{i}-\overline{y})\big{(}t_{i}-\overline{t}\big{)}}{ \sum_{i=1}^{N}\big{(}t_{i}-\overline{t}\big{)}^{2}} \tag{3}\]  
where \(D\) denotes the frequency drift rate and \(N\) the number of \(y_{i}\). \(\overline{y}\) and \(\overline{t}\) denote the mean value of \(y_{i}\) and \(t_{i}\), respectively. The evaluation of frequency stability usually uses the ADEV or HDEV method. The Allan deviation is the most common time domain measure of frequency stability [28], but it is sensitive to the linear drift of the frequency data. Hence, we choose the HDEV method to evaluate the frequency stability. The HDEV examines the second difference of the frequency data and is insensitive to the linear frequency drift. The formula reads:  
\[\sigma_{y}^{2}(\tau)=\frac{1}{6(M-2)}\sum_{i=1}^{M-2}(\gamma_{i+2}-2\gamma_{i+ 1}+\gamma_{i})^{2} \tag{4}\]where \(\gamma_{i}\) denotes the ith of \(M\) fractional frequency values with the averaging time \(\tau\)[28]. With the above formulas and the ODTS clock products, the frequency accuracy, drift rate, and stability of the atomic clock can be derived by a data processing procedure. For one clock data series, the steps are as follows:  
1. The clock differenced data between the clock data of one station or one satellite and reference clock data (AHM station) are derived, named as phase data. The sampling interval is 300 s.
2. The preprocessing of data editing is completed on a daily basis. If the data loss rate on a given day (i.e., the percentage between the lost epochs and the sum of epochs (288)) is more than 20%, this day is flagged as an unavailable day. If the phase data on a given day are available, the linear interpolation method is used to fill the phase data for those lost data on a daily basis. 3. The phase data of one day are converted to the frequency data by the first-order time difference. Then, with the frequency data, the outliers, the phase jumps, and the frequency jumps are detected and removed by using the median absolute deviation (MAD) method: \[\theta=\ \mathrm{Median}\Bigg{(}\frac{\left|y_{i}-\mathrm{Median}(y_{i}) \right|}{0.6745}\Bigg{)}\] (5) where \(y_{i}\) is the frequency value and the \(\mathrm{Median}(y_{i})\) denotes the median value of the time series of \(y_{i}\). \(\theta\) denotes the sigma of the median deviation under the condition of normal distribution. Five times \(\theta\) is regarded as the threshold to reject the bad data (i.e., if \(\left|y_{i}\right|>\mathrm{Median}(y_{i})+5\cdot\theta\), then \(y_{i}\) is treated as an outlier) [28]. 4. The lost days are counted and flagged, and then the data are fragmented. If the number of days in a continuous segment is less than seven, these data are deleted. There exists the day boundary value between two consecutive days. Thus, the connection point value is deleted and retrieved by using the linear interpolation method. Then, the clean frequency data can be derived. 5.