A High-Accuracy, High Anti-Noise, Unbiased Frequency Estimator Based on Three CZT Coefficients for Deep Space Exploration Mission

Deep space exploration navigation requires high accuracy of the Doppler measurement, which is equivalent to a frequency estimation problem. Because of the fence effect and spectrum leakage, the frequency estimation performances, which is based on the FFT spectrum methods, are significantly affected by the signal frequency. In this paper, we propose a novel method that utilizes the mathematical relation of the three Chirp-Z Transform (CZT) coefficients around the peak spectral line. The realization, unbiased performance, and algorithm parameter setting rule of the proposed method are described and analyzed in detail. The Monte Carlo simulation results show that the proposed method has a better anti-noise and unbiased performance compared with some traditional estimator methods. Furthermore, the proposed method is utilized to process the raw data of MEX and Tianwen-1 satellites received by Chinese Deep Space Stations (CDSS). The results show that the Doppler estimation accuracy of MEX and Tianwen-1 are both about 3 millihertz (mHz) in 1-s integration, which is consistent with that of ESA/EVN/CDSN and a little better than that of the Chinese VLBI network (CVN). Generally, this proposed method can be effectively utilized to support Chinese future deep space navigation missions and radio science experiments.


Introduction
In deep space exploration, the high-accuracy frequency of the downlink carrier or tone signal is a very important link to Doppler measurement, which is an essential spacecraft tracking technique that could support orbit determination of spacecraft and also provide important data for planetary radio science experiments, such as the measurement of the gravity field and planetary Radio Occultation (RO) [1]. The required accuracy of frequency estimation is typically several mHz, of which measurement accuracy is much higher than in other applications such as communications and the global navigation satellite system (GNSS). The phase-locked loop (PLL) receiver is used to estimate the frequency and phase at the deep space station, of which performance is negatively affected by the low Signal-to-Noise Ratio (SNR) and high dynamics of spacecraft [2,3]. In such situations, whether the accurate Doppler frequency can be derived or not is a critical requirement for spacecraft orbit determination. In addition to PLL, the open-loop Doppler measurement technique is also a preferable choice. In this case, the carrier or tone signal of the deep space probe is processed by modern signal theory to extract high-accuracy Doppler observables. Currently, the existing open-loop algorithms have the disadvantages of complicated algorithms or large computation [4,5] or special hardware platforms [6,7]. Chen et al. (2021) presented a Doppler frequency retrieving method via local correlation and segmented model, which was successfully applied to Chinese deep space exploration [8]. At the same time, there is a little complexity in setting the segment parameters for various dynamics.

Problem Description
An ideal single complex exponential waveform can be modeled as follows: x(n) = Ae j(2π nk 0 N fs +φ 0 ) , n = 0, 1, . . . , N − 1, k 0 = N f 0 / f s (1) where A and f 0 are the signal amplitude and frequency, respectively, and φ 0 is the initial phase, f s is the sampling frequency, and N is the sampling points number. Without the loss of generality, let A = 1. The FFT spectrum of x(n) can be calculated as following: )(k−k 0 )π] , k = 0, 1, . . . , N FFT − 1 (2) where N FFT is the length of FFT, so the frequency spectrum resolution ∆ f = f s /N FFT . Let k p be the peak position of the N FFT -point FFT spectrum, the real frequency, f 0, can be expressed as f 0 = k 0 ∆ f = (k p + δ FFT )∆ f , where δ FFT ∈ (−0.5, 0.5], meaning the frequency bias corresponding to the FFT spectrum bins. Once δ FFT is estimated, the final frequency is estimated asf 0 = (k p +δ FFT )∆ f , whereδ FFT is the estimation of δ FFT . So, our goal is to accurately estimate δ FFT . The classical estimators mentioned in Section 1 mainly use the peak FFT spectrum and its two neighbors because the three samples occupy most of the signal power [13]. Due to the fence effect and spectrum leakage of FFT, the estimation performance may be degraded when the signal frequency is located at the bin center or edge. On the other hand, the second and third peaks of the FFT samples are smaller than the maximum recurrent peak sample, which limits the estimator's ability to anti-noise. Moreover, because the frequency estimator is nonlinear, there is an SNR threshold for keeping the estimation performance [20].
Different from that, the FFT must be sampled for the entire unit circle of the z-plane uniformly, and the Chirp-Z transform (CZT) may provide a more centralized capability to perform a local z-transform that starts from an arbitrary point and samples with arbitrary uniform spacing, which is very suitable to overcome the low SNR in deep space situation [6]. The spectrum comparison of CZT and FFT is shown in Figure 1. As we can see, the fence effect and spectrum leakage of CZT are comparatively reduced; meanwhile, the spectrum resolution is higher than that of FFT. Therefore, in this paper, we propose a novel method by using the first three CZT samples.

Proposed Estimator
The Z-transform of x(n) is where M is the Z-transform length. Let z = AW −k , A= e jθ 0 , W = e jϕ 0 , we obtain the CZT of x(n), where M becomes the spectrum zooming factor, and ϕ 0 = 2πL N M , is the spectrum resolution of CZT, θ 0 = 2π k st N , is the start frequency of spectrum zooming, k st is the FFT sample index; and B is the spectrum analysis bandwidth of CZT, B = 2πL/N, L is the spectrum bandwidth factor.
We can obtain the expression of CZT spectrum by algebraic derivation, The magnitude of X(i) is , which stands for the real index of the signal frequency in the CZT samples. If i p stands for the position of peak of CZT spectrum, and k 0 = i p + δ CZT , |δ CZT | ≤ 0.5, then the magnitudes of the three largest CZT spectrum lines are, respectively: It is easy to know that L ≤ M, and N is usually set to a larger integer number for an accurate coarse frequency estimation. In Section 3, we set N = 1024 for the Monte Carlo Simulation. Therefore, it is easy to hold that L MN. In view of |δ CZT | ≤ 0.5, we can achieve the following approximation, as shown in Formula (8).
Therefore, substituting the expression in Formula (8) into Formula (7), we can obtain After the elimination of the public factors in Equation (10) and some mathematic deductions, we can find Finally, the frequency can be estimated by utilizing Equations (12)- (14).
The proposed method herein can be conducted in the following three steps: 1.
Calculate the FFT of x(n) and make the coarse frequency estimation by searching the peak position; 2.
Set the CZT parameters, including M, k st , and L, and calculate the CZT of x(n); 3.
Search the peak position of CZT and make the optimal frequency estimation by using Formulas (12)- (14).

Analysis of Unbiased Performance
Suppose that the single complex exponential waveform in Formula (1) is contaminated by Gaussian white noise, which is expressed as follows: where the noise terms w(n) is assumed to be zero mean, complex additive Gaussian white noise with a variance of σ 2 . Therefore, the SNR can be given as 1/σ 2 . The CZT of y(n) can be expressed as follows: where W(i) is the CZT of w(n). To acquire the unbiased performance of the proposed estimator, we firstly analyzed the statistical distribution of W(i), which can be expressed as Formula (17): where W R (i) and W I (i) are the real and imaginary part of W(i), respectively, denoted as: It is clear to see that W R (i) and W I (i) can be regarded as a linear combination of a series of Gaussian white noise sequences. Therefore, W R (i) and W I (i) are also the Gaussian distribution. The average value and variance of W R (i) are given by: where E [•] means the expectation operator. Then we can deduce that, In view of the Gaussian white noise and in the statistical sense, we can find Based on Formulas (7) and (8), Formulas (15) and (16) can obtained, Substituting Formula (22) into (12) and carrying out the necessary manipulations, we find that the mathematical expectation of frequency bias estimation equals its real value, Therefore, the proposed method of frequency estimation is unbiased under the Gaussian white noise condition.

Analysis of Parameter Setting
As mentioned in Section 2, B is the spectrum bandwidth of CZT, B = L∆f, and L is the number of FFT spectrum intervals. In order to reduce the calculation capacity and improve the spectrum resolution of CZT, L should be set to the minimum of reasonable values to cover the real frequency. The distribution of the FFT spectrum lines with different frequency biases, δ FFT , is displayed in Figure 2. When −0.5 ≤ δ FFT < 0, as shown in Figures 1a and 2c, the real frequency is located between the (k p − 1)th and the (k p )th spectrum line. When 0 < δ FFT ≤ 0.5, as shown in Figure 2b,d, the real frequency is located between the (k p )th and the (k p + 1)th spectrum line. Considering the random effect caused by noise and to improve the robustness of the algorithm, it is suitable to set L = 2, and the bandwidth of CZT is centered at k p .
Substituting Formula (22) into (12) and carrying out the necessary manipulations, we find that the mathematical expectation of frequency bias estimation equals its real value, Therefore, the proposed method of frequency estimation is unbiased under the Gaussian white noise condition.

Analysis of Parameter Setting
As mentioned in Section 2, B is the spectrum bandwidth of CZT, B = LΔf, and L is the number of FFT spectrum intervals. In order to reduce the calculation capacity and improve the spectrum resolution of CZT, L should be set to the minimum of reasonable values to cover the real frequency. The distribution of the FFT spectrum lines with different frequency biases, δFFT, is displayed in Figure 2. When −0.5 ≤ δFFT < 0, as shown in Figures  1a and 2c, the real frequency is located between the (kp − 1)th and the (kp)th spectrum line. When 0 < δFFT ≤ 0.5, as shown in Figure 2b,d, the real frequency is located between the (kp)th and the (kp + 1)th spectrum line. Considering the random effect caused by noise and to improve the robustness of the algorithm, it is suitable to set L = 2, and the bandwidth of CZT is centered at kp.
FFT spectrum lines with various δFFT (real frequency is marked by an oblique cross).

Numerical Simulations and Comparison
The algorithms above were implemented and simulated. The number of samples used in the simulation is N = 1024. The sampling frequency fs = 1024 Hz, so the spectrum resolution of FFT is 1 Hz. The base frequency of the signal was f0 = 120 Hz (the signal Figure 2. FFT spectrum lines with various δ FFT (real frequency is marked by an oblique cross).

Numerical Simulations and Comparison
The algorithms above were implemented and simulated. The number of samples used in the simulation is N = 1024. The sampling frequency f s = 1024 Hz, so the spectrum resolution of FFT is 1 Hz. The base frequency of the signal was f 0 = 120 Hz (the signal frequency changes from f 0 with a certain step as follows). With the CZT spectrum factor L = 2, the start and end frequencies of CZT are 119 Hz and 121 Hz, respectively. The spectrum zooming factor of CZT, M = 10, means the CZT spectrum resolution is 0.2 Hz. When the frequency bias zone is 0-0.5 Hz, and the frequency changing step is 0.025 Hz, there are 21 frequency values in total. The frequency estimation error σ f and bias ∆ f are calculated as follows: where N δ is the number of frequency values, herein N δ = 21, and N MC is the Monte Carlo simulation times, herein N MC = 10000. f i = f 0 + δ i means the true frequency value when the frequency bias changes, and δ i = (i − 1) × 0.025 Hz. The CRLB on each SNR condition can be calculated by using Formula (25) [10] (25) where N is the data length of the integration time and SNR is the Signal-to-Noise Ratio. Figure 3 displays the simulation results on the frequency estimation bias and error under three SNR conditions. As depicted in Figure 3a, the frequency bias randomly distributes around 0 Hz when the frequency bias changes and the mean bias of the three SNR conditions are −0.1918 mHz, 0.0327 mHz, and −0.0621 mHz, respectively, which shows that the proposed estimator here is unbiased. From Figure 3b, we can find that the frequency estimation error almost reaches the CRLB, especially when the SNR values are comparatively high, such as SNR = −10 dB and 0 dB here. The ratios of the frequency estimation error to the CRLB are 1.0905, 1.0173, and 1.0095 when SNR = −18 dB, −10 dB, and 0 dB, respectively.
where Nδ is the number of frequency values, herein Nδ = 21, and NMC is the Monte Carlo simulation times, herein NMC = 10000. fi = f0 + δi means the true frequency value when the frequency bias changes, and δi = (i − 1) × 0.025 Hz. The CRLB on each SNR condition can be calculated by using Formula (25) [10] ( ) where N is the data length of the integration time and SNR is the Signal-to-Noise Ratio. Figure 3 displays the simulation results on the frequency estimation bias and error under three SNR conditions. As depicted in Figure 3a, the frequency bias randomly distributes around 0 Hz when the frequency bias changes and the mean bias of the three SNR conditions are −0.1918 mHz, 0.0327 mHz, and −0.0621 mHz, respectively, which shows that the proposed estimator here is unbiased. From Figure 3b, we can find that the frequency estimation error almost reaches the CRLB, especially when the SNR values are comparatively high, such as SNR = −10 dB and 0 dB here. The ratios of the frequency estimation error to the CRLB are 1.0905, 1.0173, and 1.0095 when SNR = −18 dB, −10 dB, and 0 dB, respectively. Furthermore, the frequency estimation performances were compared with traditional algorithms, including Quinn (1994) [12], Rife (1974) [10], MacLeod (1998) [13], Aboutanios and Mulgrew (2005) [11], and Candan (2011) [15], as well as the CRLB [10]. The simulation conditions were set as mentioned above. There are 21 frequencies in total Furthermore, the frequency estimation performances were compared with traditional algorithms, including Quinn (1994) [12], Rife (1974) [10], MacLeod (1998) [13], Aboutanios and Mulgrew (2005) [11], and Candan (2011) [15], as well as the CRLB [10]. The simulation conditions were set as mentioned above. There are 21 frequencies in total for each SNR condition, and the final estimation bias and error are the mean value of all the 21 frequency conditions. Figure 4 shows the frequency estimation error and bias of the proposed and the five traditional estimators. The results in Figure 4 can be analyzed from the following three aspects.
(1) Anti-noise ability. As can be noted from this figure, there is a visible threshold effect except for the proposed method. Taking MacLeod's (1998) [13] method as an example, when SNR is higher than −13 dB, the estimation bias and error are significantly The results in Figure 4 can be analyzed from the following three aspects.
(1) Anti-noise ability. As can be noted from this figure, there is a visible threshold effect except for the proposed method. Taking MacLeod's (1998) [13] method as an example, when SNR is higher than −13 dB, the estimation bias and error are significantly decreased, and the error is very close to CRLB. While the estimation bias and error of the proposed estimator are more stable, even when the SNR is lower than −13 dB, showing that the proposed estimator has a high anti-noise ability. (2) Estimation bias. When the SNR is larger than the threshold, which is about −13 dB in this simulation, the estimation bias of MacLeod (1998) [13] and Candan (2011) [15] is significantly decreased, and the mean biases are 0.1 mHz and 1 mHz, respectively. There are obvious biases for Quinn (1994) [12], Rife (1974) [10], and Aboutanios and Mulgrew (2005) [11] under the same simulation conditions. However, the estimation bias of the proposed method is about 1 mHz when SNR = −20 dB, and the mean estimation bias of all the simulation SNR conditions is about 0.06 mHz, which means that the bias performance of the proposed method is comparatively better. The results in Figure 4b show that the proposed method is much closer to CRLB compared with other five methods.
From what has been discussed above, we may reasonably arrive at the conclusion that the proposed method has a better anti-noise ability, frequency estimation bias, and accuracy.

Results
In this section, the raw data obtained from the Mars Express and Tianwen-1 orbiter observation experiment were utilized to evaluate the frequency estimation performance when the two probes were both orbiting Mars. The observation experiments were simultaneously implemented by CDSN, which consists of three Chinese deep space stations, the Jiamusi (JM) station, the Kashi (KS) station, and the Argentina (AG) station [8].
Compared with the Monte Carlo simulation, the most significant difference in performing high-accuracy frequency estimation with digital raw data is the downlink signals of typical non-stationary due to the relative motion between the spacecraft and ground stations. The abovementioned frequency estimators, including the proposed method, are suitable for processing stationary signals with a constant frequency. Therefore, it is necessary to eliminate the Doppler Effect before frequency estimation [21].

The Elimination of Doppler Effect
Assume that the frequency of the deep spacecraft downlink signal satisfies the n-order polynomial model as follows, f (t) = a n t n + a n−1 t n−1 + . . . + a 1 t + a 0 (26) where {a i }, i = 0, 1, . . . , n are the frequency model coefficients. Considering phase is the time integral of frequency, we can achieve the corresponding phase model, ϕ(t) = 2π f (t)dt + ϕ 0 = 2π a n n + 1 t n + a n−1 n t n−1 + . . . + a 1 2 t 2 + a 0 t + ϕ 0 Now the signal model can be constructed as Formula (28), In order to obtain the model coefficients in Formula (26), we should process the measured raw data. Assume the observation time length is T, and the sampling interval is T s . The frequency of the raw data is coarsely estimated at the integration time T p , and the number of frequency estimation is K = [T/T p ], where [x] denotes the nearest integer number of x. The estimation results are remarked as: The time scale can be constructed at T p interval, t scale = (i + 0.5)T p , i = 0, 1, . . . , K − 1. Combine Formulas (26) and (29), and the frequency model can be obtained by using the Least Square Method, {â n ,â n−1 , . . . ,â 1 ,â 0 } Next, construct the time scale of the raw data at the sampling interval T s, t scale = (i + 0.5)T s , i = 0, 1, . . . , N − 1, where N = [T/T s ], denotes the total points of the sampled raw data. Letâ 0 = 0, and the phase model and signal model can be constructed as follows: Finally, we can find the residual data x res (t): where f bias (t) is the frequency bias caused by the inaccuracy of the frequency model, which can be nearly eliminated by iterative processing, when this is performed, the residual data x res (t) become a nearly stationary signal and can be processed by the proposed method, then a 0 in Formula (33) can be accurately estimated, which combines with the frequency model in Formula (30) to generate the frequency estimation of the spacecraft's downlink signal.

Mars Express Experiment
Before the first Chinese Mars exploration mission of Tianwen-1 was carried out, the China National Space Administration (CNSA) cooperated with the European Space Agency (ESA) to verify and confirm the Mars probe navigation ability of CDSN in 2020. MEX, which was launched on 2 June 2003, and has been orbiting Mars since December 2003 [22], was utilized to test and verify the feasibility of orbit measurement and orbit determination at the distance between Earth and Mars by CDSS. The data used herein were provided with the MEX observation experiments undertaken on 28 June 2020, from the JM and KS stations. The data were sampled and recorded with 0.5 MHz bandwidth and 8-bit quantification, whose format was the Delta-DOR Raw Data Exchange Format (RDEF) [23]. The sky frequency of the downlink signal is in the X band (about 8.4 GHz). The FFT spectrum of MEX at the JM station is shown in Figure 5, in which the sole peak strands for the carrier of the downlink signal. The carrier signal is utilized to estimate the high accuracy of the Doppler frequency using the proposed method.
After the elimination of the Doppler effect, the estimation of the residual frequency is shown in Figure 6a, which is the estimation of a 0 in Formula (33). In order to show the details of the stochastic characteristics, the estimation result minus its average value is shown in Figure 6b. We can see that the residual frequency results are mostly located in a region of ±10 mHz, indicating that the residual signal after Doppler effect elimination is comparatively stationary. The frequency estimation error in the presented situation is about 3.52 mHz.  The Doppler measurement results of the MEX obtained by the proposed method in this paper were compared with that of the MEX obtained by the Very Long Baseline Array (VLBA), European VLBI Network (EVN), and the Chinese VLBI network (CVN) as well as CDSN [3,6,8], respectively, as shown in Table 1. For quantitative comparison, the Dop- The Doppler measurement results of the MEX obtained by the proposed method in this paper were compared with that of the MEX obtained by the Very Long Baseline Array (VLBA), European VLBI Network (EVN), and the Chinese VLBI network (CVN) as well as CDSN [3,6,8], respectively, as shown in Table 1. For quantitative comparison, the Doppler frequency results of MEX by the proposed method were obtained within 1-s integration. The average accuracy Doppler frequency is 3.14 mHz in 1-s integration. The study by Rosenblatt et al. (2008) [3] showed that the Doppler accuracy of the MEX obtained by ESA and NASA was about 3.2 mHz in 1-s integration. These measurement results were obtained by the digital baseband receivers of ESA and NASA in the closed-loop mode. The study by Zhang et al. (2019) [6] showed that the average 7.0 mHz precision of MEX in 1-s integration was obtained by CVN. All of the above results of the MEX Doppler accuracy are displayed in Table 1; here, we can see that the accuracy of the proposed method in this paper is approximately consistent with EVN and VLBA, as well as with the previous work at CDSN, while is about two times better than CVN. Since the raw data processed in references [3,6,8] and in this paper were sampled and recorded by different ground station assemblies, the comparison results have proven to be feasible and effective for the frequency estimation of the proposed method.

Tianwen-1 Experiment
Tianwen-1 is the first Chinese deep probe to Mars in Martian science research, which was launched on 23 July 2020 [24]. The raw data were recorded when the CDSS supported the navigation mission of Tianwen-1. We selected the observation conducted by the JM and KS stations on 26 February 2021, to verify the proposed method further. At that time, Tianwen-1 was on an elliptical orbit with a periastron of~280 km, apastron of~57,815 km, and an orbital period of~49.0 h. The sampling frequency was 100 kHz, the quantification bit number was 8, and the data format was RDEF. The spectrum of Tianwen-1 is shown in Figure 7, which was observed by the KS station. The carrier signal in the spectrum is apparent. The Doppler measurement results of the MEX obtained by the proposed method in this paper were compared with that of the MEX obtained by the Very Long Baseline Array (VLBA), European VLBI Network (EVN), and the Chinese VLBI network (CVN) as well as CDSN [3,6,8], respectively, as shown in Table 1. For quantitative comparison, the Doppler frequency results of MEX by the proposed method were obtained within 1-s integration. The average accuracy Doppler frequency is 3.14 mHz in 1-s integration.
The study by Rosenblatt et al. (2008) [3] showed that the Doppler accuracy of the MEX obtained by ESA and NASA was about 3.2 mHz in 1-s integration. These measurement results were obtained by the digital baseband receivers of ESA and NASA in the closed-loop mode. The study by Zhang et al. (2019) [6] showed that the average 7.0 mHz precision of MEX in 1-s integration was obtained by CVN. All of the above results of the MEX Doppler accuracy are displayed in Table 1; here, we can see that the accuracy of the proposed method in this paper is approximately consistent with EVN and VLBA, as well as with the previous work at CDSN, while is about two times better than CVN. Since the raw data processed in references [3,6,8] and in this paper were sampled and recorded by different ground station assemblies, the comparison results have proven to be feasible and effective for the frequency estimation of the proposed method. Ref. [3] Ref. [6] Ref. [8] JM KS Average

Tianwen-1 Experiment
Tianwen-1 is the first Chinese deep probe to Mars in Martian science research, which was launched on 23 July 2020 [24]. The raw data were recorded when the CDSS supported the navigation mission of Tianwen-1. We selected the observation conducted by the JM and KS stations on 26 February 2021, to verify the proposed method further. At that time, Tianwen-1 was on an elliptical orbit with a periastron of ~280 km, apastron of ~57,815 km, and an orbital period of ~49.0 h. The sampling frequency was 100 kHz, the quantification bit number was 8, and the data format was RDEF. The spectrum of Tianwen-1 is shown in Figure 7, which was observed by the KS station. The carrier signal in the spectrum is apparent.     Table 2 displays the estimation results of both stations with different integration times, from which we can see that the Doppler estimation errors of the proposed method at the JM station are 2.97 mHz in 1-s integration, 1.86 mHz in 5-s integration, and 1.41 mHz in 10-s integration, respectively. Moreover, the Doppler estimation RMS at the KS station are 3.06 mHz in 1 s-integration, 1.85 mHz in 5 s-integration, and 1.55 mHz in 10 s-integration, respectively. The results are consistent with that of Chen et al. (2021) [8]. It is concluded that the frequency estimation results with an accuracy of about 3 mHz in 1-s integration can provide high accuracy orbit determination of the Mars probe and will be helpful for future Chinese deep space radio science experiments.  [8]. It is concluded that the frequency estimation results with an accuracy of about 3 mHz in 1-s integration can provide high accuracy orbit determination of the Mars probe and will be helpful for future Chinese deep space radio science experiments. Moreover, the Signal-to-Noise Ratios (SNR) of the received signal at both stations are also estimated, which are about 4.1 dB at JM and 2.3 dB at KS, respectively. The corresponding CRLB of frequency estimation was calculated and is displayed in Table 2. There are big gaps between the frequency estimation error and the CRLB.  Moreover, the Signal-to-Noise Ratios (SNR) of the received signal at both stations are also estimated, which are about 4.1 dB at JM and 2.3 dB at KS, respectively. The corresponding CRLB of frequency estimation was calculated and is displayed in Table 2. There are big gaps between the frequency estimation error and the CRLB.
The CRLB reflects the lower band of an unbiased estimation under the white noise condition. The thermal noise of the antenna's receiver is usually modeled as Gaussian white noise. Therefore, as the integration time progresses, the impact of thermal noise is reduced, and consequently, the frequency estimation error is lower. Considering the gap between estimation error and the CRLB, we deduce that the thermal noise is not the dominant error factor for the frequency estimation in Tianwen-1.

Error Sources Discussion
The main error sources of the Doppler estimation in deep space exploration include phase scintillation, thermal noise, and frequency stability of the oscillator at the station [25].
The total analyzed error of frequency estimation caused by the above three factors, σ f, can be calculated as follows: The phase scintillation is acquired by the downlink carrier when passing through the solar corona and introduces a random error to the Doppler measurement, which can be approximated by the following equation, θ SEP is the Sun-Earth-Probe angle (SEP) and T p is the estimation integration time. The constant parameter C band depends on the working mode and frequency band. If the spacecraft transmits a signal which is transmitted from the ground and the ground-based reference for the Doppler is the same one that drives the transmitter, the observation mode is "two-way". Three-way Doppler measurement is analogous to two-way mode, except that the downlink carrier is received at a different station than that from which the uplink carrier was transmitted. On the observation date of Tianwen-1, the KS station transmitted the uplink signal to Tianwen-1 in the X band and received the coherent frequency from Tianwen-1 in the X band, too. This means that the observation of the KS station utilized a two-way mode. At the same time, the JM station received the downlink signal in the X band by the three-way mode. In the two-way or three-way mode, the C band takes values of 5.5 × 10 −6 when both the uplink and downlink frequency are in the X band.
The SEP angle during the observation of Tianwen-1 is about 78 • , as shown in Figure 9. The corresponding frequency error caused by the phase scintillation is depicted in Figure 10. We can see that the errors are 1.754 mHz, 1.324 mHz, and 1.172 mHz for the 1 s, 5 s, and 10 s integration time, respectively.    The second error factor is thermal noise on both the uplink and downlink for the two-way and three-way mode in the Phase-locked-loop situation. As mentioned before, the thermal noise is always modeled as Gaussian white noise, and the error performance of the proposed method is no more than 1.0102 times that of CRLB when SNR is higher than 0 dB. Therefore, the white noise performance of the Tianwen-1 observation can be evaluated by the CRLB and has also been displayed in Table 2.  The second error factor is thermal noise on both the uplink and downlink for the twoway and three-way mode in the Phase-locked-loop situation. As mentioned before, the thermal noise is always modeled as Gaussian white noise, and the error performance of the proposed method is no more than 1.0102 times that of CRLB when SNR is higher than 0 dB. Therefore, the white noise performance of the Tianwen-1 observation can be evaluated by the CRLB and has also been displayed in Table 2.
The third error factor is the frequency stability of the oscillator at the station, which can be depicted as the Allan deviation of the carrier's frequency source [26]. where fsky is the sky frequency of downlink carrier, M equals to TLT/τ, and σA(τ) is the Allan deviation at time interval of τ. Therefore, the frequency estimation error caused by frequency source stability in the Tianwen-1 observation is about 0.03 mHz.  The third error factor is the frequency stability of the oscillator at the station, which can be depicted as the Allan deviation of the carrier's frequency source [26]. The Turnaround Light Time (TLT) for the Tianwen-1 observation is about 20 min. Moreover, the traditional Allan deviation of the hydrogen maser at 1000 s is about 2 × 10 −15 . The Doppler measurement error caused by frequency source stable, σ f,FS , can be calculated as the following formula, where f sky is the sky frequency of downlink carrier, M equals to TLT/τ, and σ A (τ) is the Allan deviation at time interval of τ. Therefore, the frequency estimation error caused by frequency source stability in the Tianwen-1 observation is about 0.03 mHz. Based on the above analysis, we can achieve the total analyzed error by using the Formula (34). For simplicity, only the results at the JM station are shown in Table 3. As we can see, the total analyzed errors are obviously smaller than the estimation errors. The ratios of the total analyzed error to estimation error are about 65%, 71%, and 83%, respectively. That is, the longer the integration time is, the closer the total analyzed error is to the estimation error. Since the thermal noise is related to the integration time, it is reasonable to deduce that the thermal noise performance of the proposed method is not comparable to but worse than CRLB when processing the raw data from the deep spacecraft. On the one hand, the thermal noise during the observation is possibly not ideal Gaussian white noise. On the other hand, the residual frequency of the raw data after the Doppler Effect elimination is still probably randomly changing, as depicted in Figures 6 and 8, which means that the residual signal is not ideally stationary. Therefore, CRLB could not perfectly reflect the thermal noise effect in this situation. In addition to the above three error factors, the terrestrial troposphere and ionosphere may also induce frequency estimation errors because of their temporal and spatial variety. It is known to all that tropospheric and ionospheric delay change with the observing elevation angle, changes with the motion of spacecraft. This means the tropospheric and ionospheric delay changes along with the motion of spacecraft during the observation. Figure 11 shows the tropospheric and ionospheric delay during the observation at the JM station, which is measured using the water vapor radiometer and the GNSS receiver with a sampling time of 1 s. It is easy to find that the tropospheric delay is time-varying, and the difference of tropospheric delay changes obviously with a random error of about 11.8 ps/s, which corresponds to a frequency error of about 95.6 mHz for the X band signal, much larger than that of frequency estimation error in Table 2. This is because, apart from the real variety of tropospheric delay, the water vapor radiometer has its own measuring accuracy and probably covers up the real variety of tropospheric delay. There is a similar phenomenon for the ionospheric delay but with a random error of about 7.2 mHz. Even so, it is reasonable to say that the tropospheric and ionospheric delay are also the error sources of the frequency estimation of the deep spacecraft downlink signal. In addition to the above three error factors, the terrestrial troposphere and ionosphere may also induce frequency estimation errors because of their temporal and spatial variety. It is known to all that tropospheric and ionospheric delay change with the observing elevation angle, changes with the motion of spacecraft. This means the tropospheric and ionospheric delay changes along with the motion of spacecraft during the observation. Figure 11 shows the tropospheric and ionospheric delay during the observation at the JM station, which is measured using the water vapor radiometer and the GNSS receiver with a sampling time of 1 s. It is easy to find that the tropospheric delay is timevarying, and the difference of tropospheric delay changes obviously with a random error of about 11.8 ps/s, which corresponds to a frequency error of about 95.6 mHz for the X band signal, much larger than that of frequency estimation error in Table 2. This is because, apart from the real variety of tropospheric delay, the water vapor radiometer has its own measuring accuracy and probably covers up the real variety of tropospheric delay. There is a similar phenomenon for the ionospheric delay but with a random error of about 7.2 mHz. Even so, it is reasonable to say that the tropospheric and ionospheric delay are also the error sources of the frequency estimation of the deep spacecraft downlink signal.

Conclusions
This paper presents a novel frequency estimator for Doppler measurement in deep space exploration. The proposed method is carried out with an FFT-based coarse frequency estimation and a fine estimation by utilizing the mathematic relation of the three CZT coefficients around the peak lobe. Firstly, the theoretical algorithm and signal processing procedures are described in detail. Monte Carlo simulations were implemented, and the results show that the unbiased frequency estimation error closely follows the CRLB in a lower SNR region in comparison to the previous estimators, including Rife (1974) [10], Macleod (1998) [13], Aboutanios and Mulgrew (2005) [11], and Candan (2011) [15], which indicate that the proposed frequency estimator has a better performance at

Conclusions
This paper presents a novel frequency estimator for Doppler measurement in deep space exploration. The proposed method is carried out with an FFT-based coarse frequency estimation and a fine estimation by utilizing the mathematic relation of the three CZT coefficients around the peak lobe. Firstly, the theoretical algorithm and signal processing procedures are described in detail. Monte Carlo simulations were implemented, and the results show that the unbiased frequency estimation error closely follows the CRLB in a lower SNR region in comparison to the previous estimators, including Rife (1974) [10], Macleod (1998) [13], Aboutanios and Mulgrew (2005) [11], and Candan (2011) [15], which indicate that the proposed frequency estimator has a better performance at anti-noise ability, frequency estimation bias, and accuracy. Then, the proposed method was utilized to process the received raw data of MEX and Tianwen-1 at the CDSS. The results show that the frequency estimation error of MEX and Tianwen-1 are both about 3 mHz in 1 s integration time. The accuracy of the Doppler frequency retrieving of MEX is consistent with ESA/EVN and is about two times better than CVN. Additionally, we evaluate the main error sources, including phase scintillation, frequency stability, and thermal noise, finding that phase scintillation is the dominant error source. However, there are some uncertain factors to be analyzed, such as the effects of tropospheric and ionospheric delay. Generally speaking, the proposed method herein can be effectively utilized to apply to future Chinese deep space navigation missions and can be a powerful support for radio science experiments in deep space exploration.