*5.2. Performance Analysis of Detection Based on CCAF*

**Lemma 2.** *The probability distribution of* Ψ*i*(*u*, *f*) *under H*<sup>1</sup> *hypothesis is*

$$\Psi\_i \left( (u\_\prime f)|H\_1\right) \sim \text{CN}\left( \Psi\_{R\_i R\_i} (u\_\prime f)\_\prime \frac{\sigma\_\omega^2 \beta\_i^6}{N} \right). \tag{36}$$

**Proof.** See Appendix A.2.

Using the same analysis method to analyze the *H*<sup>0</sup> hypothesis, the distribution of the detection quantity Ψ*i*(*u*, *f*) of a single GPS satellite under the *H*<sup>0</sup> hypothesis is given by

$$\Psi\_i \left( (u\_\prime f) | H\_0 \right) \sim \text{CN} \left( 0, \frac{\sigma\_\omega^2 \beta\_i^6}{N} \right). \tag{37}$$

Using the data fusion technique of multiple GPS satellites, according to the cumulative nature of the Gaussian distribution, the distribution of the final detection statistic Ω (*Rr*, *V*) is

$$\left(\Omega\left(R\_{\prime}, V\right)|H\_{0}\right) \sim \text{CN}\left(0, \sum\_{i=1}^{M} \frac{\sigma\_{\omega}^{2} \beta\_{i}^{6}}{N}\right) \tag{38}$$

and

$$\left(\Omega\left(R\_{r\prime},V\right)|H\_1\right) \sim \text{CN}\left(\sum\_{i=1}^{M} \left\{\Psi\_{R,R\_i}(u\_{\prime}f)\right\}, \sum\_{i=1}^{M} \frac{\sigma\_{\omega\prime}^{2}\theta\_i^{6}}{N}\right). \tag{39}$$

According to Equation (38), the false alarm probability can finally obtain:

$$P\_{FA} = \int\_{\lambda}^{\infty} f(\Omega | H\_0) d\Omega = \exp\left(-\frac{\lambda}{\sum\_{i=1}^{M} \frac{\sigma\_{\omega}^2 \beta\_i^{\delta}}{N}}\right),\tag{40}$$

where *f* (Ω|*H*0) is the probability density function of (Ω|*H*0).

The detection threshold *λ* in the solution Equation (40) is

$$\lambda = -\ln\left(P\_{FA}\right) \cdot \sum\_{i=1}^{M} \frac{\sigma\_{\omega}^{2} \beta\_{i}^{6}}{N}.\tag{41}$$

According to Equations (A14), (39), and (41), the detection probability can be obtained by using the signal detection theory

$$P\_D = Q\_m \left( \sqrt{\frac{\left| \sum\_{i=1}^M \beta\_i^3 a\_i \right|^2}{\sum\_{i=1}^M \frac{\sigma\_\omega^2 \beta\_i^6}{N}}}, \sqrt{\frac{\lambda}{\sum\_{i=1}^M \frac{\sigma\_\omega^2 \beta\_i^6}{N}}} \right), \tag{42}$$

where *Qm*(·, ·) is the Marcum Q function.

It can be seen from Equation (42) that the theoretical detection probability of the GPS weak echo signal detection method based on multi-star data fusion is related to the parameters, such as the monitoring channel noise, the number of sampling points *N*, the number of satellites, and the false alarm probability. It can be seen that the detection probability is proportional to the number of sampling points *N* and the number of satellites, that is, the detection probability increases with the number of sampling points *N* and the number of satellites. It is theoretically proved that the GPS weak echo signal detection method based on multi-star data fusion can improve the detection performance from the number of sampling points and the number of satellites. In addition, by comparing Equation (32) and Equation (39), it can be found that the power of the noise of the detection statistic constructed based on the CCAF algorithm decreases as the number of sampling points increases. However, the detection statistic based on the CAF algorithm does not have this property, so the CCAF-based detection statistic construction method has good noise immunity in the proposed algorithm.

## **6. Numerical Results and Discussion**

In order to analyze and verify the effectiveness of the proposed detection algorithm and the influence of various factors on detection, this section presents several simulation experiments using MATLAB (9.5.0.944444 (R2018b), MathWorks Company, Natick, MA, USA) and simulation parameter setting according to [34,35].

Experiment 1: In order to compare the detection performance of CAF and CCAF, this experiment fixes the false alarm probability, the power difference between the reference signal and the echo signal, the number of sampling points, the number of satellites, etc. The delay and Doppler shift of the echo relative to the direct wave are set as 1 us and 500 Hz, respectively. The SNR of the echo is in the range from −90 dB to −30 dB. For different detection quantity construction methods, 2000 Monte Carlo simulation experiments are carried out on the GPS weak echo detection method based on multiple satellites data fusion. The parameters are substituted into the theoretical detection probability Equations (35) and (42) to compare them with the simulation. The specific simulation parameters are shown in Table 4.

**Table 4.** Parameter setting of experiment 1.


It can be seen from Figure 14 that, under the same conditions, the simulated detection probability could reach 99% at a SNR of −55 dB. When the SNR of the CCAF algorithm is −64 dB, the simulated detection probability also reaches 99%. Therefore, the CCAF algorithm is about 9 dB better than the CAF algorithm, and the difference between simulation results and theoretical results is 2 dB, which verifies the effectiveness of the method. The above simulation results show that CCAF has a certain noise suppression capability compared with CAF because the noise does not have cyclostationarity. Moreover, according to Equation (A19), it can be seen that the noise power in the CCAF-based detection decreases as the number of points increases, while CAF does not have this property. This experiment verifies the validity of the theory, showing the noise resistance of the CCAF.

Experiment 2: This verifies the detection performance of the GPS weak echo signal detection method based on multi-star data fusion for different satellite number conditions, the false alarm probability, GPS direct wave power, and sampling points. Assuming that the target is 10 km away from the receiver, the speed is 600 m/s, and the arrival angles *θ* of the three echoes are 45◦, 50◦, 60◦, respectively. It is concluded that the time delay and Doppler shift of the echo relative to the direct wave are 10 us, 3157 Hz, 11 us, 2669 Hz, 16 us, 2467 Hz, and the echo power is 70 dB different from the direct wave power. The SNR range of the echo is from −90 dB to −30 dB. For different satellite numbers, 2000 Monte Carlo experiments are carried out on the detection method of the GPS weak echo signal based on data fusion of multiple satellites. Finally, the parameters are substituted into Equations (35) and (42), and the theoretical detection probability is compared with the simulation. The specific simulation parameters are shown in Table 5.



From Figure 15, it can be seen that, under the same conditions, the detection probability of a satellite is 99% when the echo SNR is −65 dB, and the detection probability of two satellites is 99% when the echo SNR is −71 dB, and the detection probability of three satellites is 99% when the echo SNR is −74 dB, and the simulation and the theoretical detection performances are only about 2 dB in difference. The validity of the method is verified. The above simulation results show that, compared with the traditional single satellite detection, the fusion of multiple satellite detections can effectively enhance the detection performance. On the one hand, the algorithm can effectively convert multiple GPS delay–Doppler detection peak coordinates (10 us, 3157 Hz), (11 us, 2669 Hz), (16 us, 2467 Hz) to distance–speed detection peak (10 km, 600 m/s). Furthermore, a plurality of echo peaks can be superimposed, thereby enhancing the peak value of the echo peak. On the other hand, the method uses the CCAF algorithm to construct the detection amount, which suppresses the noise caused by the cross terms. Therefore, the method enhances the detection performance from two aspects: the cumulative effect of the multi-star detection fusion on the echo peak and the de-noise of the CCAF detection structure. However, it is necessary to consider that there are only 7–8 GPS satellites in the zenith at the same time, and detection performance improves as the number of satellites increases.

**Figure 14.** Comparison of detection performance between CCAF and CAF.

**Figure 15.** Comparison of detection performance between CCAF and CAF.

Experiment 3: In order to verify the influence of the accumulated time of the signal on the detection performance, we set the parameters of the false alarm probability, the direct wave power of GPS, the number of satellites, etc. The echo parameters are set according to Experiment 2. The SNR of the echo ranges from −90 to −30 dB. The Monte Carlo experiment was performed on the proposed algorithm under the condition that the number of sampling points changed from 10<sup>6</sup> to 109 points. Finally, the parameters are substituted into Equation (35) and Equation (42), and the theoretical

detection probability is compared with the actual simulation. The specific simulation parameters are shown in Table 6.


From the results in Figure 16, it can be seen that, under the same conditions, in the case of echo SNR of −60 dB, the detection probability of the cumulative time of 0.1 s reaches 99%; in the case of echo SNR of −71 dB, the detection probability of the cumulative time of 1 s reaches 99%, in the case of echo SNR of −78 dB, the detection probability of the cumulative time of 10 s reaches 99%, and the simulated and theoretical detection performances are only 2 dB in difference, which verifies the effectiveness of the method. From the simulation results, it can be found that, as the accumulation time increases, the detection probability of the algorithm will increase. This is because CCAF is two-dimensional correlation. Therefore, the longer the correlation accumulation time, the larger the energy accumulated by the target echo will be. At the same time, as the accumulation time increases, the noise floor of the detection amount is also suppressed. Therefore, the number of sampling points *N* is one of the key factors affecting the detection performance of the proposed method. It is possible to increase the detection probability by using a longer accumulation time as much as possible, but the accumulation time cannot be increased indefinitely due to the limited detection area.

**Figure 16.** Influence of different accumulation time on detection performance.

Experiment 4: When the detection system detects the target of the specified area, the detection performance is mainly affected by the radar cross section (RCS) area of the target due to the fixed detection distance, and the change of the RCS size causes the change of the direct wave and the echo power ratio (SDR). In order to verify the influence of the attenuation of GPS signal caused by RCS on the detection performance of the proposed algorithm, the experimental fixed false alarm probability, reference signal power, noise power, number of satellites, and sampling points. The echo parameters are set according to the second experiment. Under the condition that the SDR variation range is 105 dB to 60 dB, the RCS corresponding variation range is 0.01 m<sup>2</sup> to 96.3 m2, and 2000 Monte Carlo experiments are performed on the proposed algorithm. Finally, these parameters are substituted into Equations (35) and (42) to compare with the actual simulation. The specific simulation parameters are shown in Table 7.



From the results in Figure 17, it can be seen that the detection probability can reach 99% when the difference between the direct wave and echo power reaches 80 dB. At the same time, when the detection distance is about 10 km, the weak echo detection algorithm based on the multiple satellites data fusion proposed in this paper can effectively detect the weak GPS signal with an 80 dB attenuation. Both theoretical analysis and simulation verify that the proposed algorithm can effectively detect close-range targets.

**Figure 17.** Relationship between SDR and detection probability.

## **7. Conclusions**

In this paper, a GPS weak echo signal detection method has been proposed based on data fusion of multiple satellites, in which the SIC algorithm is utilized to separate and reconstruct multiple GPS reference signals on a reference channel. Then, the detection statistic of the weak echo of multiple GPS satellites has been constructed by using the anti-interference property of CCAF. Through the coordinate transformation algorithm, multiple detection statistics have been superimposed to obtain the final detection statistic. Finally, based on the theoretical analysis of the detection statistics, the relationships between the detection performance and a series of key parameters have been studied. Numerical results have shown that the proposed method is proved to be able to effectively detect the weak echo of the target in multiple GPS scenarios and significantly outperforms the existing methods.

**Author Contributions:** M.L. conceptualized and performed the algorithm and wrote the paper; Z.G. and Y.L. analyzed the experiment data; Y.C. helped modify the language; H.S. provided technical assistance to the research; F.G. is the research supervisor. The manuscript was discussed by all co-authors. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (61501348 and 61801363), the Shaanxi Provincial Key Research and Development Program (2019GY-043), the Joint Fund of Ministry of Education of the People's Republic of China (6141A02022338), the China Postdoctoral Science Foundation (2017M611912), the 111 Project (B08038), and the China Scholarship Council (201806965031).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

*Appendix A.1. Proof of Lemma 1*

Substituting Equations (26) and (27) into Equation (14), we can obtain

$$\begin{split} S\_i(\mathbf{r}, f) &= \sum\_{n=0}^{N} \mathbf{x}\_s^{IF} \left( nT\_s \right) \hat{\mathbf{x}}\_i^{IF} \left( nT\_s - \tau \right)^\* e^{j2\pi fnT\_s} \\ &= \sum\_{n=0}^{N} \text{Re} \left\{ \mathbf{x}\_s^{IF} \left( nT\_s \right) \hat{\mathbf{x}}\_i^{IF} \left( nT\_s - \tau \right)^\* e^{j2\pi fnT\_s} \right\} \\ &+ j \sum\_{n=0}^{N} \text{Im} \left\{ \mathbf{x}\_s^{IF} \left( nT\_s \right) \hat{\mathbf{x}}\_i^{IF} \left( nT\_s - \tau \right)^\* e^{j2\pi fnT\_s} \right\}. \end{split} \tag{A1}$$

Analyzing the probability distribution of the real part in Equation (A1), and the real part in Equation (A1) is given by

$$\begin{split} \mathrm{Re}\left(S\_{i}(\tau,f)\right) &= \sum\_{n=0}^{N} \mathrm{Re}\left\{x\_{s}^{IF}(nT\_{s})\dot{x}\_{i}^{IF}(nT\_{s}-\tau)^{\*}e^{j2\pi fnT\_{s}}\right\} \\ &= \sum\_{n=0}^{N} \mathrm{Re}\left\{a\_{s}^{\*}\,\beta\_{i}p\_{i}(nT\_{s}-\tau\_{i})e^{-j(2\pi f\_{d}nT\_{s})}p\_{i}^{\*}(nT\_{s}-\tau)e^{j2\pi fnT\_{s}}\right\} \\ &+ \sum\_{n=0}^{N} \mathrm{Re}\left\{\beta\_{i}p\_{i}(nT\_{s}-\tau\_{i})e^{-j(2\pi f\_{d}nT\_{s})}\left(\sum\_{\begin{subarray}{c}i=1\\i\neq i\neq j\end{subarray}}^{N}a\_{i}p\_{i}(nT\_{s}-\tau\_{i})\right)^{\*}e^{j2\pi fnT\_{s}}\right\} \\ &+ \sum\_{n=0}^{N} \mathrm{Re}\left\{\beta\_{i}p\_{i}^{\*}\,(nT\_{s}-\tau)e^{j2\pi fnT\_{s}}\omega(nT\_{s})\right\} \\ &\approx \sum\_{n=0}^{N} \mathrm{Re}\left\{a\_{s}^{\*}\,\beta\_{i}p\_{i}(nT\_{s}-\tau\_{i})e^{-j(2\pi f\_{d}nT\_{s})}p\_{i}^{\*}(nT\_{s}-\tau)e^{j2\pi fnT\_{s}}\right\} \\ &+ \sum\_{n=0}^{N} \mathrm{Re}\left\{\beta\_{i}p\_{i}^{\*}\,(nT\_{s}-\tau)e^{j2\pi fnT\_{s}}\omega(nT\_{s})\right\} \\ &= S\_{\mathbb{R}\mathbb{R}}(\tau\_{f},f) + S\_{\mathbb{R}\mathbb{R}}(\tau\_{f},f). \end{split} \tag{A2}$$

The reason for the sign being approximately equal in Equation (A2) is that the baseband signals of different GPS satellites are not correlated with each other, and the power of the echo signals is much smaller than the noise. The first term *Srisi* (*τ*, *f*) in Equation (A2) represents the real part of the mutual ambiguity function between the reference signal and its corresponding echo signal, and belongs to the determined detection amount, which can be expressed as

$$\begin{split} S\_{r\_i s\_i}(\boldsymbol{\tau}, f) &= \sum\_{n=0}^{N} \text{Re} \left\{ a\_i^\* \beta\_i p\_i \left( n \boldsymbol{T}\_s - \boldsymbol{\tau}\_i \right) e^{-j \left( 2 \pi f\_{d\_i} n \boldsymbol{T}\_s \right)} p\_i^\* \left( n \boldsymbol{T}\_s - \boldsymbol{\tau} \right) e^{j 2 \pi f\_{d\_i} n \boldsymbol{T}\_s} \right\} \\ &= \text{Re} \left\{ a^\* \beta \chi\_{pp} \left( \boldsymbol{\tau} - \boldsymbol{n}\_\boldsymbol{\tau} \boldsymbol{T}\_{s\epsilon} f - f\_d \right) \right\}, \end{split} \tag{A3}$$

where *n<sup>τ</sup>* stands for the time delay of the echo relative to the direct wave, *fd* is the frequency offset of the echo relative to the direct wave, and *χpp* (*τ*, *f*) is the self-fuzzy function of the amplitude-normalized baseband signal *pi* (*nTs*), which can be expressed as

$$\chi\_{pp}(\tau, f) = \sum\_{n=0}^{N-1} p\_i \left( n T\_s \right) p\_i^\* \left( n T\_s - \tau \right) e^{j2\pi f n T\_s} \,\tag{A4}$$

where |*pi* (*nTs*)| = 1, *E* [*pi* (*nTs*)] = 0. When *τ* = 0, *fd* = 0, the maximum value reached by Equation (A4) is max " |*χpp* (*τ*, *fd*)| # = *N*.

The second term *Sriω*(*τ*, *f*) in Equation (A2) represents the mutual blur function between the monitoring channel noise and the reference signal, which can be regarded as the product of the Gaussian noise and the amplitude scaling factor *pi* (*nTs*), and still obeys the Gaussian distribution. Therefore, the mean and variance of the term can represent the probability distribution of the term. The mean of *Sriω*(*τ*, *f*) is given by

$$E\left\{S\_{r\_i\omega}(\tau,f)\right\} = 0\tag{A5}$$

and the variance of *Sriω*(*τ*, *f*) is given by

*Var* [*Sriω*(*τ*, *f*)] = *E* % *<sup>N</sup>* ∑ *n*=0 Re *βiω*(*nTs*)*pi* <sup>∗</sup>(*nTs* <sup>−</sup> *<sup>τ</sup>*)*ej*2*<sup>π</sup> f nTs* <sup>2</sup> & − *E N* ∑ *n*=0 Re *βiω*(*nTs*)*pi* <sup>∗</sup>(*nTs* <sup>−</sup> *<sup>τ</sup>*)*ej*2*<sup>π</sup> f nTs* <sup>2</sup> = *E* % *<sup>N</sup>* ∑ *n*=0 Re *βiω*(*nTs*)*pi* <sup>∗</sup>(*nTs* <sup>−</sup> *<sup>τ</sup>*)*ej*2*<sup>π</sup> f nTs* <sup>2</sup> & <sup>=</sup> *<sup>N</sup>* ∑ *n*=0 *N* ∑ *m*=0 *E* Re *βiω*(*nTs*)*pi* <sup>∗</sup>(*nTs* <sup>−</sup> *<sup>τ</sup>*)*ej*2*<sup>π</sup> f nTs* · *<sup>β</sup>iω*(*mTs*) <sup>∗</sup> *pi*(*mTs* <sup>−</sup> *<sup>τ</sup>*)*e*−*j*2*<sup>π</sup> fmTs* = *β*<sup>2</sup> *i N* ∑ *n*=0 *N* ∑ *m*=0 *Rωω*(*nTs* − *mTs*) · *Rpi pi* (*nTs* <sup>−</sup> *mTs*)*ej*2*<sup>π</sup> <sup>f</sup>*(*n*−*m*)*Ts* = *β*<sup>2</sup> *i N* ∑ *n*=0 *σ*2 *ωδ*(0)*Rpp*(0) = *Nβ<sup>i</sup>* <sup>2</sup>*σ*<sup>2</sup> *ω*, (A6)

where *<sup>R</sup>ωω* (*nTs* − *mTs*) is the autocorrelation of the noise, and *<sup>R</sup>ωω* (*nTs* − *mTs*) = *<sup>σ</sup>*<sup>2</sup> *ωδ* (*nTs* − *mTs*), *δ* (*nTs*) is the sequence of unit samples, and *Rpi pi* (*nTs* − *mTs*) is the autocorrelation of the normalized signal. According to Equations (A4)–(A6), under the assumption of *H*1, the distribution of the real part of *Si*(*τ*, *f*) can be expressed as

$$\operatorname{Re}\left(S\_i(\tau,f)\right) \sim N\left(\operatorname{Re}\left\{a\_i^\*\beta\_i\chi\_{pp}\left(\tau - n\_{\tau\prime}f - f\_d\right)\right\}, N\beta\_i^2\sigma\_\omega^2\right). \tag{A7}$$

The same can be obtained. Under the *H*<sup>1</sup> hypothesis, the probability distribution of the imaginary part of *Si*(*τ*, *fd*) can be expressed as

$$\operatorname{Im}\left(S\_i(\tau,f)\right) \sim N\left(\operatorname{Im}\left\{a\_i^\*\beta\_i\chi\_{pp}\left(\tau - n\_{\tau\prime}f - f\_d\right)\right\}, N\beta\_i^2\sigma\_\omega^2\right). \tag{A8}$$

Therefore, the distribution of the detection statistic under the *H*<sup>1</sup> hypothesis is

$$\mathcal{R}\left(S\_i(\tau, f)|H\_1\right) \sim \text{CN}\left(\left\{a\_i^\* \beta\_i \chi\_{\mathcal{P}^p} \left(\tau - n\_{\tau \cdot} f - f\_d\right)\right\}, \quad 2N\beta\_i^2 \sigma\_\omega^2\right),\tag{A9}$$

where *CN*(.) represents the complex Gaussian process.
