Using cross correlation to estimate Doppler frequency

Journal Pre-proofs Using cross correlation to estimate Doppler frequency Qingbao He, Yongzhang Yang, Fei Li, Jianguo Yan, Yihao Chen PII: DOI: Reference:

S0273-1177(20)30025-9 https://doi.org/10.1016/j.asr.2020.01.007 JASR 14604

To appear in:

Advances in Space Research

Received Date: Revised Date: Accepted Date:

9 October 2019 16 December 2019 9 January 2020

Please cite this article as: He, Q., Yang, Y., Li, F., Yan, J., Chen, Y., Using cross correlation to estimate Doppler frequency, Advances in Space Research (2020), doi: https://doi.org/10.1016/j.asr.2020.01.007

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 COSPAR. Published by Elsevier Ltd. All rights reserved.

Using cross correlation to estimate Doppler frequency Qingbao He Yongzhang Yang Fei Li Jianguo Yan Yihao Chen State Key Laboratory of Information Engineering in Surveying, Mapping and remote sensing, Wuhan University, Wuhan 430070, China

[email protected]; [email protected]; [email protected]; [email protected]; [email protected]; Abstract: Doppler frequency estimation with high accuracy is very important for precise orbit determination and scientific studies in deep space exploration. In this paper, we propose a new method which implements cross correlation of only received signals from one single station to estimate Doppler frequency. The algorithm is relatively simple and it can achieve high accuracy even when signal to noise ratio (SNR) is very low. Moreover, the accuracy can be further improved by implementing frequency compensation, especially when the frequency dynamic range is high. Simulations were performed and results proved the validity of this method. X-band observations on MEX (Mars Express) and Juno were performed by a 13 m telescope of Wuhan University, which was equipped by two back ends, i.e. CDAS (Chinese Data Acquisition System) and RSR (Radio Science Receiver). Raw data recorded by the CDAS were processed by using the proposed method, and the estimated frequencies were in good agreement with the real time frequency estimation values by the RSR. What’s more, the accuracy by using cross correlation was higher than the one of RSR. Keywords: Frequency estimation; Doppler; Cross correlation; Deep space exploration

1. Introduction In deep space exploration, Doppler frequency is very important for determining spacecraft’s orbit and carrying out various scientific studies, such as gravity field measurement, atmospheric experiment and ring occultation study (Asmar et al., 2005; Showman et al., 2013). To obtain Doppler frequency, it is common to estimate received frequency at ground antennas and then subtract the transmitted frequency at the spacecraft. High accuracy of Doppler frequency estimation is required for precise orbit determination and various scientific studies. However, as spacecraft goes further away, signals received by ground antennas become weaker. In addition with the high dynamic range of frequency due to relative motion between the spacecraft and ground station, achieving high accuracy of frequency estimation becomes much more challenging. Using bigger antennas to observe spacecrafts can improve signal to noise ratio (SNR) and therefore improve the accuracy of frequency estimation. However, big antennas are not always available. In this paper, we present an improved data processing method for small antennas to achieve high accuracy of frequency estimation. Signals received by ground antennas are first down converted to intermediate frequency (IF) and digitized by analog to digital converter (ADC), then further converted to baseband and saved in storage unit, which is called raw data. Raw data contains all the information of frequency, phase, and amplitude. Here we discuss the methods of estimating the received frequency by using the raw data. Phase-locked loop (PLL) tracking measurement is a well known method for 1

frequency estimation. However, PLL would lose lock when the SNR is low and the frequency dynamic range is high (Paik et al., 2011). For improving the accuracy of frequency estimation, many algorithms have been proposed using open-loop measurement (Paik et al., 2011; Chen et al., 2015), but they usually have the disadvantages of complicated algorithms and huge amount of calculation. Moreover, time delay interferometry was proposed for precise Doppler frequency estimation, which required two tracking stations (Armstrong et al., 2008). Delay correlation of a local signal and a received signal was also proposed (Dai et al., 2010).In this paper, we propose a new method, which implements cross correlation of only received signals from one single station. The new method utilizes a simple algorithm and meanwhile achieves high accuracy without a predicted frequency model. The rest of this paper is organized as follows. Section 2 presents the details of the proposed new method. Section 3 shows simulation results. Section 4 gives frequency estimation results of Mars Express (MEX) and Juno using the new method. Section 5 concludes this paper.

2. Methodology 2.1. Primary principle The basic idea of the proposed method is to implement cross correlation of raw data and use both correlation amplitude and phase to estimate frequency. Figure 1 shows the primary steps in estimating the received frequency using the raw data. The first step is to divide the raw data into many small regions, which are equal in time length. In Figure 1, we set the time length as 50 ms, and of course, it can be changed depending on how strong or weak the signals are. The second step is to implement cross correlation of every two regions. By implementing cross correlation, we first implement FFT (Fast Fourier Transform) on each of the two regions, assuming the spectrums as Y1 and Y2, then obtain the cross correlation spectrum Y12 = conj(Y1)*Y2, where conj( ) is to get the conjugated value. The third step is to pick up the signal point (maximum in amplitude in cross correlation spectrum) and obtain its correlation phase. The signal point number is used to calculate primary frequency and the correlation phase is used to calculate residual frequency. The last step is to obtain the sum of the primary frequency and the residual frequency, which equals the received frequency. In Figure 1, f1, f2 and f3 are the received frequency at the time t1, t2 and t3, respectively. To simplify the discussion in this paper, we define the procedures in Figure 1 as 50 ms cross correlation for the time length of each region is 50 ms. Theoretical illustrations will be given in the following. The first and the second small regions of data in Figure 1 can be assumed as equations (1) and (2),

y1 (t )  cos[2 ( f  f d (t ))t  s ] ,

(1)

y2 (t )  cos[2 ( f  f d (t   ))(t   )  s ] ,

(2)

Where y1(t) and y2(t) represent the two regions of data, f is the frequency transmitted from spacecraft, fd(t) is the time-varying part, φs is the initial phase, andτis the time delay, which equals the time length of y1(t). So y2(t) = y1(t+τ). Assuming Y1(f) is the frequency spectrum of y1(t), frequency spectrum of y2(t) can be written as

Y2 ( f )  Y1 ( f )e j 2 f  .

(3) 2

The condition for equation (3) being correct is frequencies in the two regions are the same. Since Doppler frequency changes with time, frequencies in the two regions are obviously different. However, the two regions share the same frequency at the intersection time t1 (see Figure 1). Therefore, equation (3) is correct only when f equals f1, which is the frequency at the intersection time t1.

50 ms

50 ms

50 ms

t1

50 ms

cross correlation

cross correlation

Pick up signal point (fmain)；

Pick up signal point (fmain)；

Pick up signal point (fmain)；

Correlation phase (fres_1)；

Correlation phase (fres_2)；

Correlation phase (fres_3)；

f1 = fmain + fres _1

……

t3

t2

cross correlation

raw data

f2= fmain + fres_2

f3= fmain + fres_3

Figure 1. The primary steps in estimating received frequency using raw data After performing cross correlation, correlation phase can be written as

12  2 ( f main  f res )  2 (n  1)  2 f res

,

(4)

where n is the signal point number in cross correlation spectrum. In reality, actual correlation phase is within ±π, which equals φ12 – 2π (n–1). So the primary and the residual frequencies can be written as

Fs (n  1) , N   2 (n  1)  12 , 2

f main  f res

(5) (6)

where Fs is the sampling rate, N is the total samples in each region data. The sum of fmain and fres is the received frequency, which is the instantaneous frequency at the time t1. In actual calculation, for the convenience, we would connect the correlation phases by 2π ambiguities and choose one signal point number as a base point. Actual calculation procedures will be presented in the simulation part. In order to be able to connect the correlation phases, the time length of one region data (τ) and the frequency changing rate (k) should meet the following condition



1 2k

.

(7)

3

2.2. Frequency compensation The received frequency is changing with time because of Doppler shift, which leads to power decentralized in frequency domain and SNR declining. Therefore the accuracy of frequency estimation would decrease. However, frequency compensation can reduce the frequency changing rate and improve the SNR. Before implementing the frequency compensation, we first use non-compensated cross correlation (see section 2.1) to obtain relatively low-accuracy frequencies. The following illustrations are under the assumption that the low-accuracy frequencies have been obtained by using 50 ms cross correlation. The frequency compensation is implemented on every two regions of data (Figure 2). The first step is to centralize the dispersed frequency to the starting frequency. In Figure 2, the starting frequencies of part A and B are f0 and f1, respectively. The second step is to convert the centralized frequency to a fixed value. Usually, we would convert them to half of the bandwidth of the received channel. The function of the second step is to reduce the changing rate of the correlation phase, which is called fringe stopping. Theoretical illustrations are given in the following.

Part A t0

Frequency com pensation

50 m s

50 m s

t1

f0

t2

f1

50 m s

f2 Part B

Frequency com pensation

Figure 2. Outline of frequency compensation To simplify the discussion, we use complex exponent expressions and use part A in Figure 2 to be an example. Data series in part A can be represented by

y(t )  e j[2 ( f0  fv (t ))t s ]

,

(8)

where f0 is the frequency at starting time (t0) of the series, fv(t) is the time-varying part because of the relative motion between the spacecraft and the ground station, and t ranges from 0 to 0.1 s. Considering the time length of compensation data is very short, fv(t) can be represented by linear approximation, fv(t) ≈ mt, where m is the coefficient of the time-varying part. Frequency changing rate can be written as k = 2m. Using the low-accuracy frequencies, we can calculate a frequency changing rate k’ and a starting time frequency f0’ by using linear fitting within a few seconds results, and the calculated frequency changing rate k’ should be very close to the theoretical one k. Then the compensation equation can be written as

ycomp (t )  y (t )  e

 j 2

k' 2 t 2

e

 j 2 f LO t

e

k k' j [2 ( f 0  f LO ) t  2 (  ) t 2  s ] 2 2

,

(9)

where fLO equals (f0’ – B/2) with B of the bandwidth of the received channel. In equation (9),

e

 j 2

k' 2 t 2

is to centralize the dispersed frequency and e 4

 j 2 f LO t

is to convert frequency to the

middle of the received channel. After the frequency compensation, cross correlation is carried out between the two regions data in part A (see Figure 2), and the frequency we obtain is

f co  f 0  f LO  (k  k ') 

.

(10)

In the case of Figure 2, τ equals 50 ms. As we can see, fco is very close to half of the bandwidth (B/2). However, fco is not the received frequency at the time t1 yet. We have to plus the converted frequency value (fLO) and the reduced frequency changing rate (k’). Therefore the received frequency at the time t1 is

f re  f co  f LO  k ' .

(11)

3. Simulation 3.1. Validity of the methodology

100

(a)

5.2025

0

4

(b)

5.2015 5.201 5.2005

-50

0

2

4 6 FFT point

8

5.2

10 x 10

f / KHz

Phase / rad

-2

10 t/s

15

20

Theory Calculated value

1040.3

0

0.3

5

(d)

1040.4

2

-4

0

4

(c)

4

1040.2 1040.1

0

5

10 t/s

15

1040

20

0

0.2

0.005

0.1

0

0 -0.1 -0.2

5

10 t/s

15

20

(f)

0.01

(e) Error / mHz

Error / mHz

x 10

5.202

50

Max point

Cross correlation power / dB

Before processing the data collected by antennas observing spacecraft, it is important to know whether the calculation results using the methodology in section 2 are strictly consistent with the theoretical ones. Numerous simulations were performed under the condition of no random noise. Here we pick up one example to show the simulation results.

-0.005 -0.01 -0.015

-0.3 0

5

10 t/s

15

-0.02

20

0

5

10 t/s

15

20

Figure 3. Simulation results with no random noise. (a) Cross correlation spectrum; (b) Signal point 5

numbers; (c) Correlation phases; (d) Theoretical and calculated values; (e) Differential values between theoretical and calculated values; (f) Differential values after 1 s integration. The signal expression was just as equation (1) shows. The transmitted frequency f was set to 1040 KHz, the Doppler shift fd(t) was set to 10.11207*t, the initial phase φs was set to 0.2 rad, and the sampling rate Fs was set to 4 MHz. The total data series was as long as 20 s, and we implemented 50 ms cross correlation of the data series. The results are shown in Figure 3. Figure 3a shows the cross correlation spectrum between two regions. In simulation, we set only one spectral line in frequency domain, thus the signal point should be the maximum point in amplitude. Figure 3b shows the signal point numbers, which are stable for a short time and then change by a step of 1 unit. The signal point number at the starting time was set as a base point and the primary frequency was calculated using equation (5). We then extracted correlation phases of 5 points, in the middle of which was the signal point, and calculated the weighted mean value of the correlation phases. The weight is the cross correlation power of each point. Figure 3c shows the weighted mean values of correlation phases. The purpose of obtaining the weighted mean value is to balance the effect of changing power of the signal point. Afterwards we connected the correlation phases in Figure 3c by 2π ambiguities and used equation (6) to calculate the residual frequencies. Figure 3d shows the sum of the primary frequency and residual frequencies and the theoretical ones. The differential values between the two are shown in Figure 3e. As we can see, the differential values are very small (within ±0.2 mHz) but have small changing trends, which are caused by the changing power of the signal point. Figure 3f shows 1 s integration values of Figure 3e. The result shows the error after 1 s integration is within ±0.015 mHz and root mean square (RMS) is 0.009 mHz, which is very small and proves the validity of the methodology in section 2. 3.2. Accuracy before frequency compensation In the actual case, data collected by antennas always contains random noise. In order to know what accuracy the proposed algorithm can achieve in the case of containing random noise, we performed simulations in six different situations. In the six situations, the common characters were random noise impressed on a cosine-wave signal, noise power (σ2) was fixed at 1 W, and SNRs were set from –30 dB to 0 dB at a step of 0.5 dB. Signal power (Ps) can be calculated by SNR*σ2. The transmitted frequency (f) was set to 1040 KHz, and the sampling rate was set to 4 MHz. The different characters were the time length (τ) for performing cross correlation and the Doppler coefficient (m). Table 1 gives an inventory of the six situations. Table 1. Six different situations considered in simulations Situation No. Time length (τ) SNR

Signal expression

Doppler coefficient

Signal power

Noise power

1

10 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 1.11207

Ps unfixed

σ2 = 1

2

10 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 25.21207

Ps unfixed

σ2 = 1

3

20 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 1.11207

Ps unfixed

σ2 = 1

4

20 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 25.21207

Ps unfixed

σ2 = 1

5

50 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 1.11207

Ps unfixed

σ2 = 1

6

50 ms

–30 ~ 0 dB

2 Ps cos  2 ( f  mt )t  0.2

m = 25.21207

Ps unfixed

σ2 = 1

6

In the six simulations, t was set from 0 to 600 s, and frequencies were calculated for 1 s integration. The differential values between calculated frequencies and theoretical ones (f + 2mt) were obtained and RMS was calculated at every SNR value. The RMS values are the accuracy of frequency estimation using the proposed algorithm in section 2. The results are shown in Figure 4.

250 10 ms fast 10 ms slow 20 ms fast 20 ms slow 50 ms fast 50 ms slow

Error / mHz

200

150

100

50

0 -30

-25

-20

-15 -10 SNR / dB

-5

0

Figure 4. Frequency calculation errors at different SNR values In Figure 4, the legend of 10 ms fast implies 10 ms cross correlations and the frequency changing rate is fast (m = 25.21207). The legend of 10 ms slow implies the similar meaning, but the frequency changing rate is slow (m = 1.11207). The two correspond to the situation No. 2 and No. 1 in Table 1, respectively. Other legends follow the same logics. From Figure 4, we know that at a certain SNR value the error becomes smaller when the time length for cross correlation is longer. This becomes more obvious when the SNR becomes lower. The primary reason is the time length (τ) in equation (6) is in denominator, which means the error is inversely proportional to the time length. Moreover, in Figure 4, the errors are almost the same between the cases of situation No.1 and situation No.2, or situation No.3 and situation No.4. However, the error in the case of situation No.6 is bigger than the one of situation No.5. The primary reason is frequency changes much more within 50 ms than within 10 ms or 20 ms, which leads to frequency expanding and signal power decentralized. From the above information we can conclude that the frequency estimation error would become very small if we use longer length of data whose frequency changing rate is very slow to perform cross correlation. 3.3. Accuracy after frequency compensation Frequency compensation can slow down the frequency changing rate and help improve the accuracy. Here we discuss the simulations of frequency compensation in the situation No.5 and No.6 in Table 1. The simulation steps are as follows. First, at every SNR value, low-accuracy frequencies were obtained by using 50 ms cross correlation. Second, we set the length of each compensation part as 100 ms, just as Figure 2 shows. Third, the frequency changing rate and the starting frequency of each compensation part were calculated by using the low-accuracy frequencies. Fourth, equation (9) was implemented on each of the compensation part data and 7

afterwards 50 ms cross correlation was performed. Only correlation phase of the signal point was extracted, for there was almost no power leakage. Fifth, differential values between the calculated frequencies and the theoretical ones were obtained and RMS was calculated. We define these results as 50 ms compensation. Afterwards, we performed these procedures again, but set the compensation part as long as 2 s in the second step and performed 1 s cross correlation in the fourth step. The results are defined as 1 s compensation. Figures 5a and 5b show the results of frequency calculation errors before and after frequency compensation in the situation No.5 and No.6 in Table 1, respectively. The green dot line shows the results of 50 ms compensation, and the red dot line shows the results of 1 s compensation. As we can see, the errors after compensation are smaller than the uncompensated ones, and it is more obvious when the SNR is lower and the Doppler shift rate is higher. The errors after frequency compensation are almost the same in Figures 5a and 5b, and the errors of 1 s compensation is much smaller than the ones of 50 ms compensation. The above information shows that the frequency compensation can indeed improve the accuracy of frequency estimation and it is more distinctive when using long length of data, such as 1 s, to perform cross correlation after frequency compensation.

40

(a)

40

50 ms slow 50 ms compensation 1 s compensation

20

10

0 -30

50 ms fast 50 ms compensation 1 s compensation

30

Error / mHz

Error / mHz

30

(b)

20

10

-25

-20

-15 -10 SNR / dB

-5

0 -30

0

-25

-20

-15 -10 SNR / dB

-5

0

Figure 5. Frequency calculation errors before and after frequency compensation. (a) Situation No. 5; (b) Situation No.6.

4. Data processing results Though the simulation results have proved the validity of this algorithm, it would be more convincing if the frequency calculation results of actual spacecraft data could also achieve high accuracy. In this section, we would verify this by showing the frequency calculation results using data from observing spacecraft MEX (Mars Express) (Patzold et al., 2016) and Juno (Tommei et al., 2015). X-band observations on Juno and MEX were performed by a 13 m telescope of Wuhan University, China on January 10, 2018 and March 5, 2018, respectively. The 13 m telescope was equipped with two back ends, which were CDAS (Chinese Data Acquisition System) (Zhu et al., 2011) and RSR (Radio Science Receiver) (Zhang et al., 2018). The CDAS set baseband bandwidth to 2 MHz and recorded data by 4 bits. The raw data was saved in a storage unit. However, the RSR outputted estimated frequencies in real time, and saved no raw data. The method of estimating frequency in real time by the RSR can be seen in (Zhang et al., 2018). Figure 6 shows the auto 8

correlation power of MEX and Juno using 1 s raw data from the CDAS. As it shows, signals received by the 13 m telescope are very weak, especially from Juno. In Figure 6, signals in the middle of the baseband are phase calibration (Pcal) signals. We use the raw data recorded by the CDAS to estimate the received frequency using the algorithm in section 2, and compare the results with the one estimated by the RSR. Detailed information will be given below. Time in this section is in Coordinated Universal Time (UTC).

(a)

MEX

120

Pcal

Auto correlation power / dB

Auto correlation power / dB

120 100 80 60 40 20 0

0

0.5

1 FFT point

1.5

Juno

Pcal

100 80 60 40 20 0

2 x 10

(b)

0

0.5

1 FFT point

6

1.5

2 x 10

6

Figure 6. Auto correlation power of MEX (a) and Juno (b). Bandwidth is 2 MHz, half FFT size is 2x106.

20

-1.1

(a)

5

(b)

-1.2 fLO / Hz

Rate / (Hz/s)

18

x 10

16 14

-1.3 -1.4

12

0.8

1

1.2

-1.5

1.4

0.8

1

1

x 10

2

1.2

1.4

(d)

1 Phase / rad

Max point

1.4

6

(c)

1 1

0 -1

1 1

1.2 t/h

t/h

0.8

1

1.2

-2

1.4

0.8

1 t/h

t/h

Figure 7. Calculation results when implementing 1 s compensation on MEX raw data recorded by CDAS. (a) Frequency changing rates; (b) Frequencies compensated on the raw data; (c) Signal point numbers; (d) Correlation phases of signal points. We first obtained relatively low-accuracy frequencies by performing 50 ms cross correlation of the raw data of MEX, and then implemented 1 s compensation using the low-accuracy 9

frequencies when processing the raw data again. The definition of 1 s compensation and calculation procedures can be found in section 3. Figure 7 shows the results when implementing 1 s compensation on MEX raw data. Figures 7a and Figure 7b show the frequency changing rates and the values of fLO (see section 2.2), which are calculated by low-accuracy frequencies. Figures 7c and 7d show the signal point numbers and the correlation phases after frequency compensation, respectively. As we can see, the signal point numbers remain unchanged and the correlation phases are close to 0, which indicates they have reached the intended results. The frequencies estimated by 1 s compensation are shown in Figure 8a in red dot line. The blue dot line in Figure 8a shows the frequencies estimated by the RSR for 1 s integration. The values shown in Figure 8a are the frequencies within the baseband, the received frequencies are equal to the values in Figure 8a plus 8420 MHz. Figure 8b shows the differential values of the two in Figure 8a, which are close to 0 and distributed almost randomly. It indicates the results from two independent back ends are in good agreement. Figures 8c and 8d show the residual values of RSR and CDAS after polynomial fitting, respectively. The RMS for the residual values of RSR is 51 mHz. However, the RMS for the residual values of CDAS is 21 mHz, which is more than 2 times lower than the one of RSR.

200

(a) RSR CDAS

f / KHz

890

Difference / mHz

900

880 870 860 850

0.8

1

1.2

100

0

-100

-200

1.4

(b)

0.8

1

t/h 200

(c)

200

RSR RMS = 51 mHz

1.4

(d)

CDAS RMS = 21 mHz

100 Error / mHz

Error / mHz

100

0

-100

-200

1.2 t/h

0

-100

0.8

1

1.2

-200

1.4

t/h

0.8

1

1.2

1.4

t/h

Figure 8. Comparison between RSR and CDAS on MEX frequency estimation results. (a) Frequencies estimated by RSR and CDAS; (b) Differential values between RSR and CDAS; (c) Residual values of RSR after polynomial fitting; (d) Residual values of CDAS after polynomial fitting. We also calculated the received frequencies of Juno (Figure 9). For Juno, we first used 200 ms cross correlation to obtain relatively low-accuracy frequencies and then implemented 1 s 10

compensation. Figure 9a shows the estimated frequencies by the CDAS (blue dot line) and the RSR (red dot line). The results shown in Figure 9a are the frequencies within the baseband, the received frequencies are equal to the values in Figure 9a plus 8404 MHz. Figure 9b shows the differential values of the two in Figure 9a. Figures 9c and 9d show the residual values of RSR and CDAS after polynomial fitting, respectively. The RMS values of them are 96 mHz and 67 mHz, respectively. The random noise in the CDAS result is clearly smaller than the one in the RSR. However, there are some variations in both Figures 9c and 9d, which should be caused by the gravity field of the Jupiter. From the above results, we can conclude that the frequency estimation values using the proposed method are not only correct, but also with high accuracy.

800.75

400

(a)

RSR CDAS Difference / mHz

f / KHz

800.7

800.65

800.6

800.55 4.45

4.5

4.55

4.6

4.65

(b)

200

0

-200

-400 4.45

4.7

4.5

4.55

400

(c)

400

RSR RMS = 96 mHz

4.65

4.7

(d)

CDAS RMS = 67 mHz

200

Error / mHz

Error / mHz

200

0

-200

-400 4.45

4.6 t/h

t/h

0

-200

4.5

4.55

4.6

4.65

4.7

-400 4.45

4.5

4.55

4.6

4.65

4.7

t/h

t/h

Figure 9. Comparison between RSR and CDAS on Juno frequency estimation. (a) Frequencies estimated by RSR and CDAS; (b) Differential values between RSR and CDAS; (c) Residual values of RSR after polynomial fitting; (d) Residual values of CDAS after polynomial fitting.

5. Conclusion A method using cross correlation of only received signals from one single station to estimate Doppler frequency was proposed, and theoretical illustrations were presented. Simulation results proved the validity of this method and showed frequency compensation could improve the accuracy of frequency estimation. Actual spacecraft data from MEX and Juno were processed by using the proposed method. The results were in good agreement with the ones of RSR, and the accuracy was even higher than the RSR. This algorithm would be very useful for small antennas to track spacecraft in deep space missions.

11

Acknowledgements We are grateful for Xuan Yang, Weitong Jing, Qingyun Deng, Peng Feng and Xi Guo using the 13 m telescope to collect the data for this study, and we appreciate Dr. H. J. Rijks for giving useful suggestions for this paper. This work was supported by LIESMARS Special Research Founding.

Additional Information The authors declare that there is no conflict of interest regarding the publication of this paper.

References Asmar, S. W., Armstrong, J. W., Iess, L., et al, 2005. Spacecraft Doppler tracking: Noise budget and accuracy achievable in precision radio science observations. Radio Science 40, 1-9. doi:10.1029/2004RS003101 Armstrong, J. W., Estabrook, F. B., Asmar, S. W., et al, 2008. Reducing antenna mechanical noise in precision spacecraft tracking. Radio Science 43, 1-7. doi:10.1029/2007RS003766. Chen, W., Huang, L., 2015. Research on open-loop measurement technique for spacecraft lecture notes. Proc. 27th conference of spacecraft TT&C Technology in China 323, 185-197. https://doi.org/10.1007/978-3-662-44687-4_17 Dai, L., Wang, Z., Wang, J., et al, 2010. Joint code acquisition and Doppler frequency shift frequency shift estimation for GPS Signals. Vehicular Technology Conference Fall (VTC 2010-Fall), 2010 IEEE 72nd, Ottawa, ON, 1-5. doi:10.1109/VETECF.2010.5594093

Paik, M., Asmar, S. W., 2011. Detecting high dynamics signals from open-loop radio science investigations. Proc. IEEE 99, 881-888. doi:10.1109/JPROC.2010.2084550 Patzold, M., Hausler, B., Tyler, G.L., et al, 2016. Mars Express 10 years at Mars: Observations by the Mars Express Radio Science Experiment (MaRS). Planetary and Space Science 127, 44-90. http://dx.doi.org/10.1016/j.pss.2016.02.013 Showman, A. P., Fortney, J. J., Lewis, N. K., et al, 2013. Doppler signatures of the atmospheric circulation on hot Jupiters. Astrophysical Journal 762(1), 1-20. doi:10.1088/0004-637X/762/1/24 Tommei, G., Dimare, L., Serra, D., Milani, A., 2015. On the Juno radio science experiment: models, algorithms and sensitivity analysis. MNRNS 446, 3089-3099. doi:10.1093/mnras/stu2328 Zhu, R., Zhang, X., Wei, W., et al, 2011. The progress of modern Chinese data acquisition system (in Chinese). Progress in Astronomy 29 (2), 207-217. doi:10.1007/s11589-011-0776-4 Zhang, T., Meng, Q., Ping, J., et al, 2018. A real-time, high-accuracy, hardware-based integrated parameter estimator for deep space navigation and planetary radio science experiments. Measurement Science and Technology 30, 1-12. https://doi.org/10.1088/1361-6501/aaedec

12

13