A multipath mitigation algorithm for GNSS signals based on the steepest descent approach

Multipath interference seriously degrades the performance of Global Navigation Satellite System (GNSS) positioning in an urban canyon. Most current multipath mitigation algorithms suffer from heavy computational load or need external assistance. We propose a multipath mitigation algorithm based on the steepest descent approach, which has the merits of less computational load and no need for external aid. A new ranging code tracking loop is designed based on the steepest descent method, which can save an early branch or a late branch compared with the narrow-spacing correlation method. The power of the Non-Line-of-Sight (NLOS) signal is weaker than that of the Line-of-Sight (LOS) signal when the LOS signal is not obstructed and with a relatively high Carrier Noise Ratio (CNR). The peak position in the X-axis of the ranging code autocorrelation function does not move with the NLOS interference. Meanwhile, the cost function is designed according to this phenomenon. The results demonstrate that the proposed algorithm outperforms the narrow-spacing correlation and the Multipath Estimated Delay Locked Loop (MEDLL) in terms of the code multipath mitigation and computation time. The Standard Deviation (STD) of the tracking error with the proposed algorithm is less than 0.016 chips. Moreover, the computation time of the proposed algorithm in a software defined receiver is shortened by 24.21% compared with the narrow-spacing correlation.


Introduction
Multipath error is related to an observation environment (Sun et al., 2019;Xu et al., 2020) and difficult to predict. It is the primary error source in Global Navigation Satellite System (GNSS) positioning in an urban canyon. Unlike the errors in satellite clock, receiver clock, ionospheric propagation, tropospheric propagation, and satellite orbits, the differential approach cannot eliminate the multipath error (Jia et al., 2017).
Three ways to mitigate the multipath interference are antenna-based technology, measurement domain processing technology, and baseband signal processing technology. The antenna-based technology makes the antenna receive the indirect signals as little as possible.
Although the placement of a GNSS antenna in a welldesigned place (McGraw et al., 2004) is the most effective multipath mitigation way, it is impossible to always have such ideal environments in urban canyons. Another effective method to mitigate multipath interference is using well-designed antennas, including single antenna and antenna array. The multipath mitigation performance with the above antennas was compared in Maqsood (2013). A single antenna, such as a dual-polarization antenna (Egea-Roca et al., 2018;Jiang et al., 2012) and a choke ring antenna (Philippov et al., 1999), can suppress the short-delay multipath. The antenna array (Wu et al., 2019;Zhang et al., 2020) requires prior knowledge of the direct signal direction and the assistance of inertial navigation. Although these well-designed antennas mitigate multipath interference effectively, they are expensive and bulky (Wang, Jong, et al., 2015).
Measurement domain processing technologies are mainly based on a physical or empirical model to

Open Access
Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: zengqh@nuaa.edu.cn mitigate multipath interference, such as terrain assistance (Kbayer et al., 2018;Kumar et al., 2015;Ng et al., 2021;Zhang et al., 2021), vector tracking (Hsu et al., 2015), Bayesian receiver autonomous integrity monitoring (Pesonen et al., 2011), satellite selection (Blanco-Delgado et al., 2010), wavelet transform (Su et al., 2018), carrier phase multipath mitigation (Moradi et al., 2015), consistency checking (Jiang et al., 2011), Wireless Fidelity (WiFi) aided multipath mitigation (Nur et al., 2013), network-based multipath mitigation (Klimenko et al., 2021), Vondrak filtering (Zheng et al., 2005), empirical mode decomposition (Dai et al., 2014), stochastic state estimation (Zhang et al., 2018), Sparsity-Promoting Regularization (Chen et al., 2019), and Signal to Noise Ratio (SNR) based multipath error correction (Axelrad et al., 1996). However, these methods require external assistance or static long-term observation (Wen et al., 2019(Wen et al., , 2020. The multipath mitigation algorithm based on baseband signal processing technology can be divided into the nonparametric estimation and the parametric estimation (Qin et al., 2019). The methods based on the nonparametric estimation mainly include the narrowspacing correlation (Van Dierendonck et al., 1992), the High Resolution Correlator (HRC) (So et al., 2009), and the strobe correlator (Liu et al., 2016). The narrow-spacing correlation algorithm mitigates multipath interference by narrowing the correlation spacing. However, the infinite bandwidth assumption is invalid in practice. The multipath error cannot be further reduced by this method when it is less than the reciprocal of the channel bandwidth. The HRC and the strobe correlator are the variants of the double delta correlator. The HRC and the strobe correlator are inefficient in mitigating short-delay multipath interference. The nonparametric estimation method has a simple structure, but its mitigation capability is not as good as the parametric estimation method.
The parametric estimation method mainly adopts a Maximum Likelihood (ML) algorithm, and the correlator array structure is the hardware and software developed with this method (Blanco-Delgado et al., 2012). Multipath Estimated Delay Locked Loop (MEDLL) (Van Nee et al., 1994;Yan et al., 2017), Multipath Mitigation Technology (MMT) (Weill, 2002), vision correlator (Fenton et al., 2005), Coupled Amplitude Delay Lock Loop (CADLL), and Enhanced CADLL (ECADLL) (Chen et al., 2010(Chen et al., , 2013 are commonly leveraged in multipath estimation. MMT and vision correlator are the improved methods of the MEDLL. Sokhandan (2013) pointed out that the main idea of the MEDLL and its improved methods is to estimate the multipath components by a recursive process and then eliminate the influence of each element on the positioning result before the following estimation. Since the principle of the MEDLL is to use multiple correlators to determine the shape of the distorted correlation function accurately,  pointed out that although the MEDLL can reduce the multipath error significantly, the upper bound of estimation accuracy is determined by the correlator spacing. The computational load is directly proportional to the number of correlators. The CADLL adds two Amplitude Locked Loops (ALL) to track the Q and I signal amplitude branches and needs a batch of units to track the code phase. Jia (2017) pointed out that the CADLL is an algorithm with a heavy computational load. Parametric estimation methods require special digital signal processing, namely Application Specific Integrated Circuit (ASIC) hardware or enormous computational complexity, making them almost impossible to implement in real-time receivers.
A code multipath mitigation tracking loop based on the steepest descent is designed to reduce the computational load and dependence on external aids. The rest of this paper is organized as follows. The steepest descent method is first introduced, followed by the principle of the proposed multipath mitigation algorithm. Then a ranging code tracking loop based on the proposed algorithm is constructed. Finally, the performance of multipath mitigation is analyzed, and the feasibility of the tracking loop is verified by tracking the real-world signals.

Principle of code multipath mitigation based on the steepest descent
Non-Line-of-Sight (NLOS) is the phenomenon in which a signal arrives at a receiver through multiple reflection paths. Due to the reflection, the NLOS signal has two properties. One is that the NLOS signal always arrives at the antenna later than the Line-of-Sight (LOS) signal because it has a longer propagation path. The other is that the power of the NLOS signal is weaker than that of the LOS signal when the LOS signal is not obstructed and with a relatively high CNR (Tamazin et al., 2016;Ziedan, 2022). For coherent tracking, the distorted autocorrelation function of ranging code is shown as where α s is the amplitude of the NLOS signal relative to the LOS signal and within the range of (−1, 1). When the NLOS signal is in-phase with the LOS signal, α s is a positive number. While α s is a negative number if the NLOS signal is out-of-phase with the LOS signal. τ s is the code delay of the NLOS signal relative to the LOS signal and within the range of (0, 1) chips for two reasons. The first one is the NLOS signal always arrives at the antenna later than the LOS signal. The second one is the peak value and position in the X-axis of the autocorrelation function are not affected by the NLOS signal when the τ s exceeds one chip. R M (·) is the distorted autocorrelation function. Supposing the signal bandwidth is infinite and the correlation function has no sidelobe, the autocorrelation function R(·) of the ranging code can be expressed as The autocorrelation function R(·) within the domain of (− 1, 1) chips is considered in the tracking process. The phase difference between the ranging code and its duplicated code is compressed into one chip by the signal acquisition module before the raw signal enters the tracking module. Substituting (2) into (1) yields The shape of the distorted autocorrelation function shown in Fig. 1 can be obtained from (3).
As seen from Fig. 1, the monotonicity of function R(·) and its maximum position in the X-axis are not affected by the NLOS signal when the power of the NLOS signal is weaker than that of the LOS signal, which is the fundamental of the proposed code multipath mitigation algorithm. Although the position of its maximum in the X-axis is not shifted in the presence of NLOS interferences, there is significant rounding/smoothing/filtering with finite bandwidth. Thus the peak is not sharp and well-defined. This is the reason why we employ the steepest descent algorithm.
To make the function R M (·) and the function R(·) have the same value range, the function R M (·) needs to be normalized. The normalization process is to divide R M (·) by its maximum, and the normalized function is shown as The autocorrelation function R(·) is not suitable for the cost function of the steepest descent because it is non-differentiable at its maximum point. Therefore, defining a new cost function in the whole region of interest is necessary, which is differentiable and has only one extreme point. The quadratic function meets these requirements and is simple in control stability analysis, so the new cost function is defined as where P F (·) is the new cost function. Substituting the distorted autocorrelation function R M (·) into (5) yields The shape of the distorted cost function shown in Fig. 2 can be obtained from (6).
As seen in Fig. 2, the minimum point of the cost function P F (·) is also the stable point of the steepest descent method, and the location of the minimum does not move with the NLOS interference. Although the distorted cost function is non-differentiable at points τ = −1 + τ s , τ = τ s and τ = 1 , it has both left and right partial derivatives at these points. Whether in-phase or out-of-phase NLOS signals distort the cost function, they can meet the requirements of the steepest descent method.
The principle of the steepest descent (Zeng et al., 2021) is to minimize cost function P F (·) by adjusting τ . According to the state update process of the steepest descent, the iteration process of τ can be expressed as where ∂P F ∂τ | k is the derivation of P F (·) with respect to τ at the moment. µ is the control step, which is a dimensionless positive scalar. To analyze the stability of the iterative process τ in the absence of NLOS interferences, substituting (5) into (7), the iteration process of τ can be expressed as τ k can converge to 0 if and only if |1 − 2µ| < 1 . The iteration process of τ is in an over-damped state when µ ∈ (0, 0.5) , in an under-damped state when µ ∈ (0.5, 1) , and in a critical-damped state when µ = 0.5.

Ranging code tracking process based on the proposed algorithm
The principle of code multipath mitigation and the steepest descent method have been introduced previously. In this section, a new ranging code tracking loop is designed according to (7) and shown in Fig. 3.
As seen in Fig. 3, the proposed ranging code tracking loop is divided into four steps. The first step is to mix the received IF signal with two orthogonal local signals. The second step is to multiply the early code and the punctual code, or the late code and the punctual code generated locally with the signal entering the second step. The third step is an integral process, and the fourth step is to control the ranging code NCO based on the steepest descent. These four steps are described in detail below.
Step 1: Carrier mixing It is assumed that the structure of discrete IF signal input to the tracking loop is shown as where A is the amplitude of the signal x(n) , C(n) is the ranging code sequence, D(n) is the message code sequence, w I is the angular rate, and θ 0 is the initial phase. The carrier mixing process is to mix the IF signal with two orthogonal local signals and shown as where θ 1 is the initial phase of the local carrier signal.
Step 2: Strip the ranging code away If the left partial derivative of the cost function is applied in the iteration process of τ , the local late code and the local punctual code are utilized in this step. Similarly, if the right partial derivative of the cost function is adopted in the iteration process of τ , the local early code and the local punctual code are employed. Taking the iteration process based on the left partial derivative as an example, the stripping process of ranging code is shown as i(n) = (x(n))(cos(w I n + θ 1 )) = AC(n)D(n) cos(θ 0 − θ 1 ) + AC(n)D(n) cos(2w I n + θ 0 + θ 1 ) q(n) = (x(n))(sin(w I n + θ 1 )) = − AC(n)D(n) sin(θ 0 − θ 1 ) + AC(n)D(n) sin(2w I n + θ 0 + θ 1 ) Step 2 Step 3 Step 4 Fig. 3 Proposed ranging code tracking loop. If the iteration process of τ adopts the left partial derivative of the cost function, this method needs the late branch L and the punctual branch P . If the iteration process of τ employs the right partial derivative of the cost function, this method utilizes the early branch E and the punctual branch P . IF is the abbreviation of intermediate frequency. NCO is the abbreviation of numerical controlled oscillator where R(·) is the autocorrelation function of the ranging code, τ is the phase difference between the local ranging code and the modulation code in the received signal, and d is the punctual-late spacing.
Step 3: Integral process An integrator has the function of low-pass filtering, and the process is independent of navigation data since the integration time is set as 1 ms. So the integrated signal is shown as where T is the integration time, which is usually set as an integral multiple of the ranging code period, and f s is the sampling rate.
Step 4: The iteration process of τ.
To reduce the influence of the carrier on the ranging code tracking loop, the I / Q branch signal should be processed as follows The amplitude A of signal x(n) can be regarded as a constant in a short time, and the peak value of the autocorrelation function R(·) is 1. Therefore, the maximum value of S L or S P in a short time can be applied to normalize (13), and the normalization process is shown as where S max is the maximum value of S L or S P in a short time. According to the definition of the cost function P F (·) in (5), the new cost function P F (·) can be described as The iteration process of τ can be obtained by taking (15) into (7), which can be described as where P F (S P )| k is the cost value of the punctual correlator at the moment. P F (S L )| k is the cost value of the late correlator at the moment. ∂P F ∂τ | k − is the left derivation of P F (·) with respect to τ at the moment. Similarly, if the right partial derivative of the cost function is employed, the iteration process of τ can be expressed as where P F (S E )| k is the cost value of the early correlator at the moment. ∂P F ∂τ | k + is the right derivation of P F (·) with respect to τ at the moment.

Experiments
To verify the effectiveness of the proposed algorithm, a set of tests was conducted with the BeiDou Navigation Satellite System (BDS) B1I signal Zhao et al., 2022) (13) and Global Positioning System (GPS) (Wang, Ji, et al., 2015) L1 signals. The simulation test verifies the code multipath mitigation performance of the proposed algorithm. The real-world short-term test and long-term test confirm the practicability of code multipath mitigation. Considering the simple structure of the narrow-spacing correlation and the high tracking accuracy of the MEDLL, the narrow-spacing correlation and the MEDLL algorithm were adopted to evaluate the multipath mitigation performance. The punctual-late spacing employed for the narrow-spacing correlation is set to 0.05 chips, and the correlator numbers leveraged for the MEDLL algorithm are set to 65.

Simulation test
We designed a series of IF satellite signals interfered by NLOS signal or noise based on (1) to verify the performances in the following three situations.
• If no NLOS signal is combined with the LOS signal, the influence of the control step µ on the tracking performance under different CNRs is discussed. • If only one NLOS signal is added to the LOS signal, and the CNR of the LOS signal is set as 43 dB·Hz, the impact of the amplitude of the NLOS signal on the tracking error is analyzed. • If only one NLOS signal is combined with the LOS signal, and the noise and bandwidth are not consid-ered, the code multipath mitigation performance is studied with the multipath error envelope.
The simulation parameters are given in Table 1, and the influence of the control step on tracking error under different CNRs is illustrated in Fig. 4. Figure 4 demonstrates the influence of the control step on tracking error under different CNRs. The curvature of the cost function changes gently when the tracking error approaches zero. The tracking error is difficult to converge to zero if a small control step is adopted. The influence of noise on tracking error is the smallest when the control step is within the range of [0.6, 0.8] . When the control step is within the range of [0.2, 0.9] and the CNR is 43 dB·Hz or 53 dB·Hz, the tracking error of less than the set punctual-late spacing 0.1 chips has a probability of over 95%. The impact of the amplitude of the NLOS   Table 1, the NLOS signal delay is set as 0.2 chips, and the CNR is set as 43 dB·Hz. The results are shown in Fig. 5. Figure 5 shows the effect of the amplitude of the NLOS signal on the tracking error. The tracking error increases with an increase in the amplitude of the NLOS signal. When the system is in the under-damped state, the tracking error is less than that for the over-damped state and within 0.1 chips with a probability over 95%.
Since the multipath error envelope (Townsend et al., 1995) is always utilized to evaluate the code multipath mitigation performance, we adopt this curve to verify the advantages of the proposed algorithm. If only one NLOS signal is combined with the LOS signal and the noise and bandwidth are not considered, the multipath signal delay varies from 0 to 1.5 chips with an increment of 0.062 5 chips. The multipath error envelope is shown in Fig. 6. Figure 6 reveals the proposed algorithm outperforms other methods in code multipath mitigation. It is also remarkable that the MEDLL is better than the narrow-spacing correlation. The proposed algorithm is based on the phenomenon that the peak position in the X-axis of the ranging code autocorrelation function does not move with the NLOS interference. In theory, ignoring the noise and bandwidth, the tracking accuracy of the proposed algorithm is not affected by NLOS signals. The multipath error envelope of the proposed algorithm in Fig. 6 verifies this conclusion.   Fig. 9. The information on the satellites observed by antenna A and B is shown in Table 2.
According to the information in Table 2, satellites Pseudo Random Noise (PRN) number 3 and 17 can be detected only by antenna A, but in reality both antenna A and B simultaneously got the signals from satellites PRN3 and PRN5. This indicates that the signals of these satellites received by antenna C are interfered by NLOS signals. The four satellite signals collected by an IP-solutions sampler connected to antenna C were processed, and the results are shown in Fig. 10.
From Fig. 10, the MEDLL needs a longer convergence time than the narrow-spacing correlation and the proposed algorithm. The maximum tracking error of the three algorithms is less than 0.2 chips in the stable tracking state. The tracking accuracy of the proposed algorithm is better than that of the MEDLL, and that of the narrow-spacing correlation is the worst. To analyze the performance of the three algorithms, we listed their tracking errors, including Standard Deviation (STD) and mean error (Mean), in the stable tracking state in Table 3. Table 3 tells that the mean error of the three multipath mitigation algorithms is negligible, and the tracking accuracy of the proposed algorithm is better than those of the

Real-world long-term test
Cannon (2001) mentioned the oscillation of the longterm measurement data existed in multipath environments. If similar oscillations can be observed in two sets of satellite observation data with an interval of integer orbital periods, the data is likely to be interfered by NLOS signals. Therefore, two long-term measurements were carried out at antenna C in the real-world shortterm test. The first measurement was conducted from 03:40 to 06:10 Coordinated Universal Time (UTC) on July 14, 2021. The second measurement was carried out from 03:36 to 06:06 UTC, on July 15, 2021. Two measurements were separated by an interval of two orbital periods of 23 h and 56 min. The measurement CNR of PRN 6 is shown in Fig. 11. Figure 10 demonstrates the repeatability of the first and second measurement data and the oscillation from the start to the second 9 000, which indicates that the received signal of PRN 6 is interfered by NLOS signals. The IF signals of PRN 6 collected by the IP-solutions sampler connected to antenna C in the second long-term measurement were processed. The results are shown in Fig. 12. Figure 12 exhibits two basic phenomena. The first one is that the amplitude of the tracking error with the proposed algorithm is stable during the whole tracking stage, which is consistent with the fact that the proposed algorithm is not affected by NLOS signals in theoretical analysis. The second one is that the STDs of the tracking error with the narrow-spacing correlation, the MEDLL, and the proposed algorithm are 0.088 chips, 0.031 chips, and 0.016 chips, respectively. The proposed algorithm has a smaller tracking error than the other two methods. Figure 3 indicates that the proposed algorithm needs the late branch and the punctual branch, or the early branch and the punctual branch. The punctual branch needs to be adopted for carrier tracking. The punctual branch adopted in the proposed ranging code tracking loop is also employed for carrier tracking. Compared to the narrow-spacing correlation that early branch and late branch are not adopted by carrier tracking, our algorithm saves the late branch or the early branch for GNSS receivers. To verify the improvement in computation efficiency due to saving one branch, we recorded the time taken with the three multipath mitigation algorithms in offline processing the collected 12 s data in the short-term test. The three multipath suppression algorithms run on the Intel  (R) Core (TM) i5-9400F CPU@2.9 GHz processor. The computation time is shown in Table 4. Compared with the narrow-spacing correlation and the MEDLL, the proposed algorithm takes the least time. The narrow-spacing correlation method takes 13.708 s to offline process the 12 s data, while only 10.390 s to process the same data with the proposed algorithm. The proposed algorithm reduces 24.21% of the computation time compared with the narrow-spacing correlation.

Conclusions
A new code multipath mitigation algorithm is proposed based on the steepest descent method for GNSS signals without any external assistance. The proposed algorithm has a shorter computation time than the narrow-spacing correlation method and better multipath mitigation performance than the MEDLL. The experiments demonstrate that the proposed algorithm shortens the computation time by 24.21% compared with the narrow-spacing correlation and improves the tracking accuracy by 27.78% compared with the MEDLL. However, the proposed algorithm has poor anti-interference ability because it is based on the assumption that the LOS signal power is greater than the NLOS signal power.