A novel Doppler frequency measurement method based on the closed-loop signal correlation for deep space exploration

Signal correlations have been used for Doppler shift measurements since long time ago. For example, in the application of measuring blood flow velocity, the phase increment of the latter signal relative to the previous signal can be calculated through the autocorrelation of signals at adjacent moments, and the blood flow velocity can be obtained (Bonnefous & Pesqué 1986). Signal correlations are also widely used in other fluid velocity measurements (Cadel & Lowe 2015). Signal correlations also have some important applications in other fields, including high precision measurement of the speed of moving objects (Hirata & Kurosawa 2012), and microwave Doppler sensors (Khablov 2017).

The cross-correlation method is often used in astronomical research to measure the relative velocity of celestial bodies. By calculating the cross-correlation between the signal of celestial bodies and the signal of the template, the velocity of celestial bodies can be measured (Allende Prieto 2007). In the technical field of astronomy, the key principle of very long baseline interferometry (VLBI) is also the cross-correlation method. The delay difference between the signals of the same detector received by different VLBI stations can be measured through the cross correlation operation, and then the delay difference can be used for the orbit determination (Schuh & Behrend 2012).

Recently, a cross correlation has been used in carrier tracking for deep space exploration, and an open-loop cross-correlation Doppler frequency shift estimation algorithm has been obtained. The frequency measurement results of this method are compared with that of the actual Radio Science Receiver, and the tracking accuracy of Doppler frequency is slightly improved (He et al. 2020). However, the shortcoming of the open-loop crosscorrelation method is that when the sampling rate of the signal is high, the calculation amount of the algorithm is too large, and the measurement accuracy of the algorithm will be significantly reduced at high dynamics.

Phase-locked loop (PLL) is a classical high precision carrier tracking method, which is widely used in radar, deep space exploration and other fields. In PRIDE (Duev et al. 2016) (Planetary Radio Interferometry and Doppler Experiment) of Joint Institute for VLBI in Europe (JIVE), phase-locked loop is used in carrier tracking. Shanghai Astronomical Observatory also developed a Doppler measurement method by a second order phaselocked loop (Ma et al. 2016). In order to realize high precision Doppler frequency and phase estimation, and to meet the needs of further distance and richer forms of deep space exploration, it is necessary to improve tracking performance of PLL.

In our study, a phase detector is taken as the key module to improve the performance of PLL. At present, an arctangent phase detector is most widely used in PLL due to its largest convergence range and the smallest phase error (Di et al. 2014). As is seen in Figure 1, the function of the traditional phase detector is to extract the phase difference between the input signal X₁(t) and the reconstructed signal X₂(t). Because the coherence of adjacent signals' spectrums can be used to identify phases, therefore, a novel phase detector can be designed to replace the arctangent phase detector in PLL. We refer to the newly designed detector as "coherent phase detector". Compared with the open-loop cross-correlation algorithm, the new closed-loop signal correlation algorithm proposed in this paper can greatly reduce the computation especially when sampling rate is high.

Zoom In Zoom Out Reset image size

The first section of this paper is introduction. Second, the realization of the coherent phase detector is detailed presented. Third, Cramer-Rao lower bound (CRLB) and optimal tracking accuracy of frequency are calculated. Fourth, simulation results are given to verify this proposed algorithm. Fifth, the processing results of Longjiang-2 are shown. Finally, we give a summary.

The corresponding structure of PLL based on a coherent phase detector is shown in Figure 2. The input signal and the reconstructed signal of the two branches are respectively mixed, and the mixed signal then pass through integrate/dump filter. The difference of phase difference between input signal and reconstructed signal at adjacent moment could be obtained by coherent phase detector. The noise in the output of phase detector is filtered out by a loop filter. The reconstructed phase is obtained by integrator and can be used to reconstruct signal. Frequency initialization is performed by using the open-loop frequency maximum likelihood estimation method (Satorius et al. 2003). In most cases, this method can make the estimation accuracy of the initial frequency within 1 Hz when Fourier transform (FFT) resolution is set to 1 Hz, which is conducive to frequency acquisition.

Zoom In Zoom Out Reset image size

Figure 3 shows the basic principle of a coherent phase detector. The time lengths of d₁(t) and d₂(t) are equal. In Figure 3, it assumes that the phase increment of the latter signal d₂(t) relative to the preceding signal d₁(t) is Δφ and Δφ can be identified. So a new type of phase detector is obtained. This kind of phase detector can also be used as a key module of phase locked loop. The derivation process of the coherent phase detector is shown in Figure 3. Let the form of the input signal be

$$\begin{eqnarray}&&X(t)=A\cos (2\pi ft),\,t\gt 0.\end{eqnarray} \tag{ 1 }$$

The actual expression of the input signal is X(t) = Acos(2πf_c t + ϕ), f_c is carrier frequency, ϕ is initial phase. In the expression of Equation (1), we take the initial phase ϕ as part of 2πft. Over time, the value of f will approach the real carrier frequency f_c and both become nearly the same. The form of Equation (1) is convenient for subsequent derivation. The reconstructed in-phase (I) component of input signal is ${S}_{{\rm{I}}}(t)=\cos (2\pi \hat{f}t)$, while the reconstructed quadrature (Q) component is ${S}_{{\rm{Q}}}(t)=-\sin (2\pi \hat{f}t)$. After mixing the input signal and the in-phase component, the output is

$$\begin{eqnarray}\begin{array}{ll}{D}_{{\rm{I}}}(t) & ={S}_{{\rm{I}}}(t)\times X(t)\\ & =\displaystyle \frac{A}{2}\cos (2\pi (f-\hat{f})t)+\displaystyle \frac{A}{2}\cos (2\pi (f+\hat{f})t).\end{array}\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

After mixing the input signal and quadrature component, the output is

$$\begin{eqnarray}\begin{array}{ll}{D}_{{\rm{Q}}}(t) & ={S}_{{\rm{Q}}}(t)\times X(t)\\ & =\displaystyle \frac{A}{2}\sin (2\pi (f-\hat{f})t)-\displaystyle \frac{A}{2}\sin (2\pi (f+\hat{f})t).\end{array}\end{eqnarray} \tag{ 3 }$$

We assume the period used in integrate/dump filter in Figure 2 is T₀, and the function is summing up the data for each continuous period of time (T₀). F_s is the sampling rate of the data at the output of integrate/dump filter. If T₀ is set to 25 μs, it means that the 25 μs data is accumulated to get an output data. The effect of integrate/dump filter is equivalent to a low-pass filter, which can eliminate interference signals of high frequency. Therefore, the signal form after integrating and dumping operation can be expressed as follows. The in-phase component

$$\begin{eqnarray}&&{d}_{{\rm{I}}}(t)=\displaystyle \frac{A}{2}\cos (2\pi (f-\hat{f})t).\end{eqnarray} \tag{ 4 }$$

The quadrature component

$$\begin{eqnarray}&&{d}_{{\rm{Q}}}(t)=\displaystyle \frac{A}{2}\sin (2\pi (f-\hat{f})t).\end{eqnarray} \tag{ 5 }$$

Construct a variable

$$\begin{eqnarray}&&Z(t)={d}_{{\rm{I}}}(t)+j{d}_{{\rm{Q}}}(t)=\displaystyle \frac{A}{2}{e}^{j(2\pi (f-\hat{f})t)}.\end{eqnarray} \tag{ 6 }$$

Discrete form of Z(t) is

$$\begin{eqnarray}&&{Z}_{{\rm{k}}}=Z(k{T}_{0})=\displaystyle \frac{A}{2}{e}^{j(2\pi (f-\hat{f})k{T}_{0})}.\end{eqnarray} \tag{ 7 }$$

Phase compensation is then performed for Z_k. Compensation formula is

$$\begin{eqnarray}&&{Z}_{{\rm{k}}}^{{\rm{c}}}={Z}_{{\rm{k}}}{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}k{T}_{0})}.\end{eqnarray} \tag{ 8 }$$

The spectral resolution of the Fourier transform of M₁ is ${\delta }_{{\rm{f}}}=\displaystyle \frac{{F}_{{\rm{s}}}}{N}$ and N is the length of M₁, M₁ is defined below. Phase compensation can shift the spectrum of the sequence. The reason for the phase compensation is that the position of the spectral line with the maximum amplitude can be at the non-edge of the obtained coherent spectrum, which is conducive to further signal processing.

Define ${M}_{1}=\{{Z}_{{\rm{m}}}^{{\rm{c}}},{Z}_{{\rm{m}}+1}^{{\rm{c}}},\cdots,{Z}_{{\rm{m}}+{\rm{N}}-1}^{{\rm{c}}}\}$. First, wavelet soft threshold denoising is used to preprocess M₁ to eliminate part of the noise and the discrete Meyer wavelet (Dmey) is selected here. Then apply Fourier transformation to M₁, the detail expression is

$$\begin{eqnarray}\begin{array}{ll}{F}_{1}(n) & ={\rm{FFT}}({M}_{1})={\rm{FFT}}({Z}_{{\rm{m}}}^{{\rm{c}}},{Z}_{{\rm{m}}+1}^{{\rm{c}}},\cdots,{Z}_{{\rm{m}}+{\rm{N}}-1}^{{\rm{c}}})\\ & =\displaystyle \frac{A}{2}\displaystyle \sum _{k=0}^{N-1}{e}^{j2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}(m+k){T}_{0}}{Z}_{{\rm{m}}+{\rm{k}}}{e}^{-j(\displaystyle \frac{2\pi n(k+m)}{N})}\\ & =\displaystyle \frac{A}{2}\displaystyle \sum _{k=0}^{N-1}{e}^{j2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}(m+k){T}_{0}}{Z}_{{\rm{m}}+{\rm{k}}}{e}^{-j(\displaystyle \frac{2\pi nk}{N})}\\ & =\displaystyle \frac{A}{2}{e}^{j(2\pi \Delta fm{T}_{0}+2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}m{T}_{0})}\\ & \,\,\,\times \displaystyle \sum _{k=0}^{N-1}{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})k-j\displaystyle \frac{2\pi nk}{N}}\\ & =\displaystyle \frac{A}{2}{e}^{j(2\pi \Delta fm{T}_{0}+2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}m{T}_{0})}\\ & \,\,\,\times \displaystyle \frac{1-{e}^{jN(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j2\pi n}}{1-{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j\displaystyle \frac{2\pi n}{N}}},\\ \Delta f & =f-\hat{f}.\end{array}\end{eqnarray} \tag{ 9 }$$

Define ${M}_{2}=\{{Z}_{{\rm{m}}+{\rm{N}}}^{{\rm{c}}},{Z}_{{\rm{m}}+{\rm{N}}+1}^{{\rm{c}}},\cdots,{Z}_{{\rm{m}}+2{\rm{N}}-1}^{{\rm{c}}}\},\displaystyle \frac{m}{N}\in {Z}^{+}$, M₂ is next to M₁ and both have equal length. Assume that M₂ = e^jΔφ M₁. Wavelet soft threshold denoising is also used to preprocess M₂ to eliminate part of the noise. Then Apply FFT transformation to M₂, the expression is

$$\begin{eqnarray}\begin{array}{ll}{F}_{2}(n) & ={\rm{FFT}}({M}_{2})={e}^{j\Delta \varphi }{\rm{FFT}}({M}_{1})\\ & =\displaystyle \frac{A}{2}{e}^{j\Delta \varphi }{e}^{j(2\pi \Delta fm{T}_{0}+2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}m{T}_{0})}\\ & \times \displaystyle \frac{1-{e}^{jN(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j2\pi n}}{1-{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j\displaystyle \frac{2\pi n}{N}}},\\ \Delta f & =f-\hat{f}.\end{array}\end{eqnarray} \tag{ 10 }$$

Multiply F₁(n) by F₂(n), then a coherent function of two frequency spectrums is obtained.

$$\begin{eqnarray}\begin{array}{ll}C(n) & ={F}_{1}(n){F}_{2}(n)\\ & =\displaystyle \frac{{A}^{2}}{4}{e}^{j\Delta \varphi }{e}^{j(4\pi \Delta fm{T}_{0}+\pi N{\delta }_{{\rm{f}}}m{T}_{0})}\\ & \times \displaystyle \frac{{(1-{e}^{jN(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j2\pi n})}^{2}}{{(1-{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j\displaystyle \frac{2\pi n}{N}})}^{2}}.\end{array}\end{eqnarray} \tag{ 11 }$$

It can be obtained from the analysis of Equation (11) that, when $(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-\displaystyle \frac{2\pi n}{N}=2q\pi,q\in Z$, the amplitude of C(n) is largest. Then we obtain that ${e}^{j(4\pi \Delta fm{T}_{0}+\pi N{\delta }_{{\rm{f}}}m{T}_{0})}=1$ and $\displaystyle \frac{{(1-{e}^{jN(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j2\pi n})}^{2}}{{(1-{e}^{j(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-j\displaystyle \frac{2\pi n}{N}})}^{2}}={N}^{2}$. In this way, a relationship can be further obtained from Equation (11).

$$\begin{eqnarray}&&C({n}_{0})=\displaystyle \frac{{N}^{2}{A}^{2}}{4}{e}^{j\Delta \varphi }.\end{eqnarray} \tag{ 12 }$$

n₀ is the frequency point corresponding to the maximum amplitude of C(n). In this case, the relative phase increment between M₂ and M₁ is argument of C(n₀). In fact, we know that the frequency resolution ${\delta }_{{\rm{f}}}=\displaystyle \frac{{F}_{s}}{N}$ is large. At this time, there is still some small deviation between the obtained value n₀ and the actual value n, and n satisfies the equation $(2\pi \displaystyle \frac{N}{4}{\delta }_{{\rm{f}}}{T}_{0}+2\pi \Delta f{T}_{0})-\displaystyle \frac{2\pi n}{N}=2q\pi,q\in Z$. For this reason, two points adjacent to n₀ are also used to calculate Δφ. For the complex Fourier transform of N points, $0\sim \displaystyle \frac{N}{2}$ represents positive frequencies, $\displaystyle \frac{N}{2}+1\sim N$ represents negative frequencies. The actual negative frequency has to be subtracted by N + 1. Therefore, the Δφ calculation strategy is as follows:

If ${n}_{0}\in (\displaystyle \frac{N}{2},\displaystyle \frac{N}{2}+1,N,1)$

$$\begin{eqnarray}&&\Delta \varphi =\Delta {\varphi }_{{{\rm{n}}}_{0}}=\arg (C({n}_{0})).\end{eqnarray} \tag{ 13 }$$

If ${n}_{0}\notin (\displaystyle \frac{N}{2},\displaystyle \frac{N}{2}+1,N,1)$, then

$$\begin{eqnarray}&&\Delta \varphi =\displaystyle \frac{\displaystyle \sum _{i=-1}^{1}\Delta {\varphi }_{{{\rm{n}}}_{0}+{\rm{i}}}|C({n}_{0}+i)|}{\displaystyle \sum _{i=-1}^{1}|C({n}_{0}+i)|}.\end{eqnarray} \tag{ 14 }$$

Δφ can be used as the input of the loop filter. On the other hand, in order to prevent integer cycle ambiguity to appear in phase discrimination, coherence time T, that is, the corresponding time length of M₁ or M₂ used in coherent computation should satisfies inequality |2πΔf₁ T| < 2π and Δf₁ is the relative frequency increment between M₁ and M₂. Therefore, it is necessary to ensure that the first derivative of the Doppler frequency f_d satisfies the relation

$$\begin{eqnarray}&&|{f}_{{\rm{d}}}|=\left|\displaystyle \frac{\Delta {f}_{1}}{T}\right|\lt \displaystyle \frac{1}{{T}^{2}}.\end{eqnarray} \tag{ 15 }$$

As is seen in Figure 4, θ₁(t) is the phase of the input signal, θ₂(t) is the phase of the reconstructed signal. Compared to the traditional phase detectors shown in Figure 1, the phase difference δφ between input signal and output signal is not directly obtained by coherent phase detector. Instead, the difference of phase difference δφ at adjacent moments is obtained. The discrete linear structure of PLL in this paper is shown in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The measured Doppler frequency is

$$\begin{eqnarray}&&{f}_{2}=\displaystyle \frac{{w}_{1}}{2\pi T}.\end{eqnarray} \tag{ 16 }$$

The transfer function of phase error

$$\begin{eqnarray}&&{\phi }_{{\rm{e}}}=\displaystyle \frac{T{s}^{3}+{s}^{2}}{T{s}^{3}+(1+2{w}_{{\rm{n}}}T){s}^{2}+2{w}_{{\rm{n}}}^{2}Ts+{w}_{{\rm{n}}}^{3}T}.\end{eqnarray} \tag{ 17 }$$

When the input phase is θ(t) = 0.5at², in which a is a constant, the steady-state error is

$$\begin{eqnarray}&&{e}_{2}(\infty )=\mathop{\mathrm{lim}}\limits_{s\to 0}s\cdot \displaystyle \frac{a}{{s}^{3}}{\phi }_{{\rm{e}}}=\displaystyle \frac{a}{{w}_{{\rm{n}}}^{3}T}.\end{eqnarray} \tag{ 18 }$$

Although the JR (Jaffe Rechtin) filter is a second-order filter, the PLL is still a second-order PLL due to the use of coherent phase detector. The closed-loop transfer function of the phase-locked loop is

$$\begin{eqnarray}\begin{array}{ll}H(s) & =1-{\phi }_{{\rm{e}}}\\ & =\displaystyle \frac{2{w}_{{\rm{n}}}T{s}^{2}+2{w}_{{\rm{n}}}^{2}Ts+{w}_{{\rm{n}}}^{3}T}{T{s}^{3}+(1+2{w}_{{\rm{n}}}T){s}^{2}+2{w}_{{\rm{n}}}^{2}Ts+{w}_{{\rm{n}}}^{3}T}.\end{array}\end{eqnarray} \tag{ 19 }$$

The calculation formula of loop noise bandwidth (Wang 2009) is

$$\begin{eqnarray}&&BW=\displaystyle {\int }_{0}^{\infty }|H(j2\pi f){|}^{2}df.\end{eqnarray} \tag{ 20 }$$

Write the input signal in another form

$$\begin{eqnarray}\begin{array}{ll}X[n] & =A\cos (2\pi {f}_{{\rm{a}}}n+\phi )+\omega [n],\\ n & =0,1,\cdots,K-1;\\ {f}_{{\rm{a}}} & =\displaystyle \frac{{f}_{{\rm{c}}}}{{f}_{{\rm{s}}}},\end{array}\end{eqnarray} \tag{ 21 }$$

where A is amplitude, K is the number of points that are sampled independently, f_s is sampling rate of original signal, f_c is carrier frequency, ϕ is initial phase, and ω[n] is a Gaussian noise which mean is 0 and variance is σ². In PLL, the accuracy of frequency estimation is almost unaffected by amplitude and initial phase. Then when calculating the CRLB of the estimation of f_c, we assume that the amplitude and the initial phase are known. The derivation process of CRLB is given below.

The probability density function of X[n] is

$$\begin{eqnarray}\begin{array}{ll} & p(X;{f}_{{\rm{a}}})=\displaystyle \frac{1}{{(2\pi {\sigma }^{2})}^{\displaystyle \frac{N}{2}}}\\ & \times \exp {\{-\displaystyle \frac{1}{2{\sigma }^{2}}\displaystyle \sum _{n=0}^{K-1}[X[n]-A\cos (2\pi {f}_{{\rm{a}}}n+\phi )]}^{2}\}.\end{array}\end{eqnarray} \tag{ 22 }$$

Take the derivative of the likelihood function

$$\begin{eqnarray}\begin{array}{ll} & \displaystyle \frac{\partial lnp(X;{f}_{{\rm{a}}})}{\partial {f}_{{\rm{a}}}}=-\displaystyle \frac{2\pi A}{{\sigma }^{2}}\\ & \times \displaystyle \sum _{n=0}^{K-1}[X[n]\sin (2\pi {f}_{{\rm{a}}}n+\phi )-\displaystyle \frac{A}{2}\sin (4\pi {f}_{{\rm{a}}}n+2\phi )]n.\end{array}\end{eqnarray} \tag{ 23 }$$

And

$$\begin{eqnarray}\begin{array}{ll} & \displaystyle \frac{{\partial }^{2}lnp(X;{f}_{{\rm{a}}})}{\partial {f}_{{\rm{a}}}^{2}}=-\displaystyle \frac{4{\pi }^{2}A}{{\sigma }^{2}}\\ & \times \displaystyle \sum _{n=0}^{K-1}[{n}^{2}X[n]\cos (2\pi {f}_{{\rm{a}}}n+\phi )-{n}^{2}A\cos (4\pi {f}_{{\rm{a}}}n+2\phi )].\end{array}\end{eqnarray} \tag{ 24 }$$

Fisher information is

$$\begin{eqnarray}\begin{array}{ll}I({f}_{{\rm{a}}}) & =-E\left[\displaystyle \frac{{\partial }^{2}lnp(X;{f}_{{\rm{a}}})}{\partial {f}_{{\rm{a}}}^{2}}\right]\\ & =\displaystyle \frac{2{\pi }^{2}{A}^{2}}{{\sigma }^{2}}\displaystyle \sum _{n=0}^{K-1}[{n}^{2}-{n}^{2}\cos (4\pi {f}_{{\rm{a}}}n+2\phi )]\\ & \approx \displaystyle \frac{2{\pi }^{2}{A}^{2}}{{\sigma }^{2}}\displaystyle \sum _{n=0}^{K-1}{n}^{2}\\ & =\displaystyle \frac{{A}^{2}}{{\sigma }^{2}}\displaystyle \frac{{\pi }^{2}K(K-1)(2K-1)}{3}.\end{array}\end{eqnarray} \tag{ 25 }$$

Then

$$\begin{eqnarray}\begin{array}{ll}var({f}_{{\rm{c}}}-\hat{{f}_{{\rm{c}}}}) & =var({f}_{{\rm{a}}}-\hat{{f}_{{\rm{a}}}}){f}_{{\rm{s}}}^{2}\ge {f}_{{\rm{s}}}^{2}/I({f}_{{\rm{a}}})\\ & =\displaystyle \frac{{\sigma }^{2}}{{A}^{2}}\displaystyle \frac{3{f}_{{\rm{s}}}^{2}}{{\pi }^{2}K(K-1)(2K-1)}.\end{array}\end{eqnarray} \tag{ 26 }$$

$\hat{{f}_{{\rm{c}}}}$ represents an estimate of f_c and $\hat{{f}_{{\rm{a}}}}$ represents an estimate of f_a. Finally, the optimal root mean square (RMS) of frequency is obtained

$$\begin{eqnarray}\begin{array}{ll}rms({f}_{{\rm{c}}}-\hat{{f}_{{\rm{c}}}}) & ={f}_{{\rm{s}}}\sqrt{var({f}_{{\rm{a}}}-\hat{{f}_{{\rm{a}}}})}\\ & \ge {f}_{{\rm{s}}}\sqrt{\displaystyle \frac{{\sigma }^{2}}{{A}^{2}}\displaystyle \frac{3}{{\pi }^{2}K(K-1)(2K-1)}}.\end{array}\end{eqnarray} \tag{ 27 }$$

The square root of CRLB will be used as an absolute criterion to evaluate the tracking accuracy of proposed algorithm.

4.1. Tracking Accuracy and Robustness Analysis

Simulation parameters setting: sampling rate f_s = 4 MHz, the initial frequency of the signal f₀ = 1 MHz. The definition formula of the carrier to noise ratio (CNR) generated in simulation is

$$\begin{eqnarray}&&{\rm{CNR}}=10{\mathrm{log}}_{10}(\displaystyle \frac{{P}_{{\rm{s}}}}{{N}_{0}})=10{\mathrm{log}}_{10}(\displaystyle \frac{{A}^{2}}{2{\sigma }^{2}{T}_{{\rm{s}}}}).\end{eqnarray} \tag{ 28 }$$

Here P_s is signal power, and N₀ is the noise spectral density, and A is amplitude of signal, and σ² is the variance of the noise and T_s is sampling period of input signal. For phase locked loop based on coherent phase detector, the period of integrate/dump filter T₀ is set to 250 μs and coherence time T, which is loop updating time, is set to 5 ms.

The phase model is defined as

$$\begin{eqnarray}&&\theta (t)=\phi +2\pi {f}_{0}t+2\pi \displaystyle \frac{{f}_{1}}{2}{t}^{2}+2\pi \displaystyle \frac{{f}_{2}}{6}{t}^{3}.\end{eqnarray} \tag{ 29 }$$

The units of f₀, f₁, f₂ are Hz, Hz s⁻¹, Hz s⁻², respectively. ϕ is the initial phase of input signal.

As is seen in Figure 6, the phase-locked loop algorithm based on coherent phase detector can accurately track Doppler frequency of the signal, which proves the correctness of the algorithm. Figure 6 shows the tracking accuracy of the signal frequency under different CNRs. In the simulation experiments of Figure 6, the loop noise bandwidth was set to decrease gradually from 4.35 Hz to 0.13 Hz (ω_n decreases from 20 to 1) in the first 5 seconds. The data of 65 seconds were tracked each time, and when caculating tracking accuracy of frequency, that is root mean square of residual frequency, the measured frequency of the first five seconds were eliminated.

Zoom In Zoom Out Reset image size

In Figure 6, when setting the loop noise bandwidth to 0.13 Hz, even though the third derivative of the phase cannot be ignored, the estimation precision of the frequency of 1 second integral still approaches the square root of CRLB. At the same time, with the bandwidth set to 0.13 Hz, the frequency of relatively high-dynamic signals can still be measured by this algorithm. The results illustrate that the measured frequency here is a nearly minimum variance unbiased estimate for the measured frequency of 1 second integration. Although the CRLB derivation assumes that the frequency is completely unknown, and in this new PLL, the frequency estimate at the previous time is a priori information of the frequency for the current second. This leads to the assumption in the CRLB derivation that the frequency is completely unknown is not entirely correct. However, in general, when set proper parameters for applications, the possible highest accuracy of the actual estimated frequency is only slightly higher than CRLB due to the long integration time of the frequency as long as 1 second. Therefore, the CRLB can still be used as an absolute benchmark to evaluate the accuracy of the measured frequency of 1 second integration. For the estimated frequency of 5 s and 60 s integration, the optimal estimation accuracy of frequency which is square root of CRLB cannot be achieved by using the current algorithm proposed by this paper.

As is seen in Figure 6, when the CNR of signal is larger than 30 dB-Hz, the signal frequency can be tracked stably under the current parameter settings. Actually, the CNR of actual data is often larger than 30 dB-Hz in deep space exploration. In addition, the dynamic of the simulation signal is large enough to basically cover the application scenarios of carrier tracking of X-band (8.4 GHz) and S-band (2.2 GHz) signals in most cases.

4.2. Computing Time

In order to further verify the reliability of the new proposed Doppler measurement algorithm, the observation data of Longjiang-2 on 2018 August 24 with 2 bits sampling were processed in this paper. Longjiang-2 is a minor satellite launched alongside Queqiao relay satellite of Chang'e-4. The sampling rate and bandwidth of the signals received by the telescopes after down-conversion are 4 MHz and 2 MHz, respectively. During the tracking process, set the loop noise bandwidth to 0.13 Hz when the loop becomes stable. In observation, the VLBI station alternate observations of Longjiang-2 and radio source. These data were observed by VLBI station in Kunming.

Figure 8 shows the results of estimated CNR of the signal. The process of estimating CNR of actual signals is as follows: First, the phase of the signal is measured using the algorithm proposed in this paper. Second, the measured phase is used to compensate the raw data so that the signal energy is concentrated within 1 Hz. Third, the power spectrum of the compensated data is calculated. Fourth, the maximum value of the power spectrum is used as an estimate of the total power of the useful signal. Fifth, take some amplitude values except the maximum value in the power spectrum, and calculate its mean value as the estimate value of noise power. Finally, Equation (28) is used to calculate the estimated value of CNR of actual signal.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The auto power spectrums of data before and after phase compensation are both shown in Figure 9. It can be seen that the energy of signal is obviously more concentrated after phase compensation.

Zoom In Zoom Out Reset image size

For these observation data, Longjiang-2 operated in coherent mode while generating signal, which means that frequency standard of Longjiang-2 is hydrogen atomic clock on the ground. Because the frequency of the atomic clock is highly stable, the error introduced by atomic clock can be ignored for the measured frequency of 1 second integration. Then it can be considered that the frequency error is mainly caused by the thermal noise, and the noise caused by the ionosphere and atmosphere. Then the phase noise can be approximated as the Gaussian process. The main steps to analyze the accuracy of frequency are as follows: First, estimate the signal frequency of the 1 second integration. Second, because the loop bandwidth slowly decreases when the frequency of first 5 seconds is estimated in each scan, the estimated frequency of the first 5 seconds are removed to ensure the correctness of subsequent RMS estimation. Third, the frequency residuals are obtained by fitting the estimated frequency values with the sixth order polynomial in each scan. Fourth, the RMS of frequency residuals is calculated and this RMS can be used as estimate of accuracy of frequency.

In addition, the mean value of estimated CNRs of each scan is calculated. Then the mean value is substituted into Equation (27) to obtain the estimate of the optimal estimation accuracy of frequency which is square root of CRLB. The performance of the current algorithm is evaluated by comparing the optimal estimation accuracy with the actual estimation accuracy of frequency.

As seen in Figure 10, the RMS of frequency obtained by novel method is very close to the square root of CRLB and the relative error between them is less than 22%. This result illustrates that the estimated frequency obtained by this method is almost the minimum variance unbiased estimator when processing actual data. In addition, Shanghai Astronomical Observatory has developed an engineering software for carrier tracking (Zheng et al. 2013) and the software has been used in recent deep space exploration missions, including the first Chinese Mars exploration and Chang'e-5 lunar exploration mission. We compare the tracking accuracy of measured frequency obtained by the software with that of the new method presented in this paper. It is clear that the new method can achieve nearly 6 times higher tracking accuracy of frequency.

Zoom In Zoom Out Reset image size

Data from Chang'e-5 on 2020 November 28 with 2 bits sampling were recorded in VLBI data processing center. The sampling rate and bandwidth of these data are 4 MHz and 2 MHz, respectively. The differences between the velocity of precision ephemeris obtained by orbit determination and the measured velocity are shown Figure 11. The velocity residuals further illustrate that the algorithm proposed in this paper is feasible for frequency measurement. At about 14:30 Coordinated Universal Time, an orbit transfer of Chang'e-5 occurred, which resulted in the velocity residual is relatively large.

Zoom In Zoom Out Reset image size

In this paper, a coherent phase detector is designed based on the phase relationship between adjacent signals in the time domain. Simulation and actual data processing results show that the coherent phase detector can obviously improve the estimation accuracy of Doppler frequency and the frequency estimate is close to the minimum variance unbiased estimate when setting proper parameters.

The weak signal tracking capability of the phase-locked loop can be further improved by extending the coherence time T, but dynamic tracking performance of the loop filter will deteriorate due to the increasing of time interval, and the assumption that the phase increment of the latter segment relative to the former segment become not necessarily true for the two adjacent signals. Then the dynamic tracking capability of PLL becomes poor and the frequency may not be locked. Therefore, the limitation of dynamic tracing capability should be taken into account when increasing the coherence time.

The research presented in this paper has been supported by the National Natural Science Foundation of China (11773060, 11973074, U1831137, 11703070 and 11803069), the National Key Basic Research and Development Program (2018YFA0404702), Shanghai Key Laboratory of Space Navigation and Positioning (3912DZ227330001) and the Key Laboratory for Radio Astronomy of CAS.