Abstract
Multipath propagation is a major source of error in global navigation satellite systems (GNSSs), especially in urban environments. This is because conventional GNSS receivers can provide biased range estimates that lead to positioning errors. In this paper, an Extended Kalman Filter (EKF)-based solution that relies on a multi-correlator structure is proposed to replace the conventional delay locked loop (DLL). The underlying signal model incorporates the radio propagation channel between the satellite and receiver and thus inherently accounts for reflected signal replicas. The algorithm was evaluated with simulations and hardware emulations and a measurement campaign was conducted in road traffic with different multipath environments. Our results revealed that the proposed EKF was very effective against multipath-associated error.
1 INTRODUCTION
Global navigation satellite systems (GNSSs) are critical for a variety of applications, including surveying, time synchronization of mobile telecommunication networks, and power grids. The use of GNSSs in aeronautics, maritime, and automobile navigation has been expanding. This technology is also critical for the navigation of autonomous vehicles. Current GNSSs already provide reliable and precise positioning under routine conditions. However, conventional GNSS receivers are vulnerable to multipath propagation that can degrade their performance or even lead to failure. As a result, the use of GNSS technology in safety-critical applications (e.g., autonomous cars or unmanned aerial vehicles) remains challenging and requires additional countermeasures.
Multipath propagation is the term used to describe a radio channel in which a signal reaches the user via more than one path. This can lead to the superposition of a direct line-of-sight (LOS) satellite signal and multiple delayed replicas of itself at the receiver. Under these circumstances, the correlation of the received signal with the local replica recorded by the GNSS receiver will lead to distorted results. A conventional delay locked loop (DLL) cannot handle these types of signals and might provide erroneous code-tracking results. This may lead to position errors of up to multiple meters (Breivik et al., 1997; Kos et al., 2010; van Nee, 1992).
Many solutions have been proposed to reduce multipath errors, ranging from pure detection to active multipath mitigation. For example, signal quality monitoring (SQM) techniques, such as Delta or Ratio metrics (Phelts, 2001), are commonly used to detect simple multipath events. These methods can be used to detect signal anomalies and to reveal the distorting effects of multipath propagation. Satellites contributing to these anomalies can then be weighted less in the position, velocity, and time (PVT) solution or simply excluded altogether.
Several more advanced receiver tracking loops are available that extend beyond pure detection. One well-known approach is known as Narrow Correlator™ that was first described by Van Dierendonck et al. (1992). Using this method, the Early-Late spacing of the correlators for the DLL is reduced and the influence of multipaths is thus lowered. Other common techniques include those of the double-delta correlator (DDC) class, e.g., the high-resolution correlator (HRC) described by McGraw & Braasch (1999) or the strobe correlator described by Garin & Rousseau (1997). These approaches add little additional computational complexity, but exhibit poor noise performance (McGraw & Braasch 1999).
Active mitigation of multipath requires more complex solutions. The underlying concept featured in the following approaches involves an estimate of the signal parameters of each of the individual multipaths and the use this information to reduce their impact on the PVT solution (Pany, 2010). One of the most widely used approaches is the multipath estimating delay locked loop (MEDLL) (van Nee et al., 1994). This algorithm relies on a multipath-incorporating signal model that samples the correlation function with a correlator bank and matches it to the expected outcome from the signal model using maximum likelihood (ML) criteria. The resulting ML estimator determines the delays and complex amplitudes of the individual multipaths. Additional replicas are then subtracted from the received signal. The remaining signal (ideally, without multipath errors) is then processed by a conventional DLL-based GNSS receiver. A considerable amount of work has been done to reduce the complexity involved in solving the ML problem. Among others, these efforts led to the development of multipath mitigation technology (MMT) (Weill, 2002). Fenton & Jones (2005) used this technique to develop the vision correlator, which is a method that decomposes LOS and multipath in the domain of the chip transitions. Furthermore, Antreich et al. (2005) proposed a method that featured the use of the space alternating generalized expectation-maximization (SAGE) algorithm for iterative approximations of the ML estimator. Bayesian solutions have also been proposed (Lentmaier et al., 2008). One drawback shared by of these advanced algorithms is that they all require an accurate replica of the received signal. Differences between the signal forms that are actually received, for example, and those that are due to unknown radio frequency (RF) filter characteristics at the receiver front-end, may be mistaken by the algorithms as multipath. The computational complexity also remains high, even with the simplifications and approximations discussed in the text above. Of note, Chen et al. (2010) proposed the coupled amplitude delay locked loop (CADLL). This method features consecutive tracking units in which each unit tracks one signal propagation path and subtracts it from the signal for subsequent units. However, instabilities can arise when more units than propagation paths are used.
Here, we present a solution to this problem that involves only moderate complexity. Iltis (1990) proposed an extended Kalman filter (EKF)-based multi-correlator structure for multipath suppression that samples the correlation function with a correlator bank. The outputs are then processed by an EKF to replace the conventional code tracking structure with DLL and early and late correlators. The EKF includes a signal model that incorporates the radio propagation channel with the help of a tapped-delay line channel. This inherently facilitates the consideration of multipath and other signal distortions, including RF front-end characteristics. When using this method, it is less crucial to have an accurate local replica of the received signal. Moreover, because this procedure provides an estimate of both the code delay and the channel impulse response (CIR), classification and characterization applications are enabled, e.g., for the detection of multipath propagation or repeater signals. This concept was explored further by Iliopoulos et al. (2017). Likewise, Siebert et al. (2021a) re-examined and refined the approach and demonstrated its potential with more sophisticated simulation and testing. Now that Galileo has reached full operational capability (FOC) and the deployment of the new L1 Civil (L1C) navigation signal has been initiated as part of the modernization of the Global Positioning System (GPS), the use case of multi-signal GNSS receivers has increased. To address this issue, Siebert et al. (2021b) proposed an extension to multiple signals accompanied by testing with actual measurement data.
In this work, we refine and verify the algorithm presented by Siebert et al. (2021b). We propose an additional scaling of the measurement noise matrix to account for the effects of approximations made in the signal model. Multipath error envelopes have been determined and used to quantify the resulting multipath mitigation capability. The solution was then tested in a dynamic environment with actual measurements. For this purpose, we evaluated the results of a measurement campaign that was conducted near Munich, Germany involving a test vehicle that was driven through a suburban area.
The remaining sections of this paper are organized as follows. The derivation of the signal model is presented in Section 2, followed by the description of the algorithm in Section 3. Simulation and hardware emulation results, as well as the evaluation of the measurement campaign that demonstrated the capabilities of the algorithm are presented in Section 4. Our conclusions are presented in Section 5.
2 SIGNAL MODEL
In this section, we describe the derivation of a multi-signal GNSS model. This model was originally described by Siebert et al. (2021b) in a publication that was an extension of Siebert et al. (2021a) and Iliopoulos et al. (2017). For these calculations, it is assumed that a navigation satellite of a GNSS constellation can broadcast Nsig navigation signals on a single carrier frequency fc. The number of navigation signals per band would be, for example, Nsig = 3 when using the GPS coarse/ acquisition (C/A) code and the data and pilot component of the L1C signal. The i-th navigation data stream di(t) is spread in frequency before transmission by a pseudorandom noise (PRN) sequence as shown in Equation (1):
1
where Tcode, i is the code period of the i-th signal. Each PRN code consists of individual chips with duration Tc. As all currently deployed civil GNSSs have equal chip duration in one band, Tc is assumed to be equal for all Nsig. Before the raw PRN code was used to spread the navigation message, it modulated itself with the sequence mi(t). The legacy GPS L1 C/A code uses a binary phase-shift keying (BPSK) for this purpose to generate a modulating sequence that is simply 
. Because the Galileo Open Service (OS) relies on binary offset carrier (BOC) modulation schemes, mi(t) is the corresponding sub-carrier signal. Furthermore, the individual chip sequences 
undergo a chip-shaping that is modeled with the impulse response g(t). The resulting modulated PRN signal is referred to as si(t) as shown in Equation (2):
2
where * denotes the convolution operator. The modulated navigation signals are then transmitted on their carrier frequency fc and correspondingly down-converted on the receiver side. The relative velocity between the satellite and the user introduces a Doppler shift that leads to a frequency shift as well as a shortening or stretching of the received PRN signals. To maintain readability, this effect was omitted in the derivation. However, we will revisit this topic in Subsection 3.3.3. Furthermore, we assumed that the navigation data bits were di, (t) = 1, 
for the same reason. The signal propagation path between the satellite and the user introduces a temporal delay τ(0)(t) that is equal for all Nsig. This is also referred to as the code delay. Additional distortions are introduced by channel characteristics (e.g., multipath) or hardware imperfections. Assuming that these characteristics can be modeled by a filtering operation, the complex received baseband signal can be expressed as indicated in Equation (3):
3
where, in Equation (4):
4
with the zero mean complex Gaussian noise 
. The channel characteristics were approximated by a time-varying tapped delay line channel model as shown in Equation (5):
5
where 
, l = −Lh,…,Lh are the complex channel coefficients and Th is the channel tap spacing with in total Ntap = 2Lh + 1 channel taps. Correspondingly, the CIR spans over ±WCIR with WCIR = LhTh. This value is equal for all Nsig as all signals propagate with the same carrier frequency fc along the same path. In contrast to the methods used by Iliopoulos et al. (2017); Iltis (1990), we allowed channel coefficients with l < 0, despite the fact that GNSSs are causal systems. Imperfect estimates of the code delay τ(0) result in a sinc-shaped CIR, which can only be accurately reflected with an acausal value for CIR. The received baseband signal was then expressed as indicated in Equation (6):
6
Next, the baseband signal was brought to the discrete signal space by sampling it at fs = 1/Ts ≥ 2B, with B representing the one-sided bandwidth of y(t), thereby fulfilling the Nyquist-Shannon sampling theorem. Therefore, we approximate the blockwise constant channel coefficients 
and code delay 
with 
and the integration time Tint. By doing so, we assumed that the channel coherence time was larger than Tint; this is a valid assumption for most applications. Moreover, without loss of generality, we limit Tint to common multiples of all Tcode,i. Thus, the periodicity of all si(t) will be preserved over the integration time Tint, thereby facilitating simpler mathematical notation in the following equations. Considering the aforementioned periodicity, this led us to the following sampled received baseband signal as shown in Equation (7):
7
where n = 0, …, N– 1 and 
. Without loss of generality, the sampling frequency fs was assumed to be chosen such that N will be an integer. For each time step k, N samples were stacked to a column vector before further processing. The signal in matrix notation is then given as shown in Equation (8):
8
with [·]⊤ denoting the transposition of a matrix or a vector. The noise vector and the vector of PRN signals were defined as shown in Equations (9), (10), and (11):
9
10
11
while GNSS receivers are designed to operate simultaneously with multiple signals from the same carrier, the phase locked loop (PLL) can hold only one signal in-phase. Because signal components are not necessarily received in-phase, other signals may be left in quadrature as a result of this process. To account for this in the signal model, the PRN signals si(t) were defined as complex. This completes the signal model. In Section 3, these principles will be used to explain the proposed algorithm.
3 PROPOSED ALGORITHM
This section describes a multi-correlator-based EKF with improved multipath suppression that was designed to replace conventional DLL for code tracking. This algorithm is an extension of the one first described by Siebert et al. (2021b), which itself presented a further development of the work described by Siebert et al. (2021a) and Iliopoulos et al. (2017) based on an approach that was initially proposed by Iltis (1990). In the following sections, we first define the correlator bank. The EKF was then described in detail in Subsection 3.2. Implementation-related information will be provided in Subsection 3.3.
3.1 Correlator Bank
Before transmission, the navigation data streams were spread in frequency with the PRN signals. As a first step, this process must be reverted in the receiver by correlating the received signal yk with the corresponding PRN signal si(t). For most conventional receivers, it is sufficient to sample the correlation function at three points, i.e., early, prompt, and late. Given that our objective was multipath suppression, we needed to collect more information. Accordingly, we introduced a multi-correlator approach with a correlator bank that sampled the correlation function more frequently. We included one bank for each of the Nsig superimposed signals per satellite. Each bank, as described in Equation (12):
12
consisted of Ncorr = 2Lc + 1 signal-matched correlators placed symmetrically and equidistantly within ±Wbank around the prompt with Wbank = LcΔτ. The hat notation indicates here and in the following the estimate of a parameter. The Nyquist condition determines the maximal correlator spacing in the bank with Δτ ≤ 1/(2B). Falling below this threshold, i.e. conducting an oversampling, is not associated with an information gain as the correlators become increasingly correlated. The correlation with the received signal results then in Equation (13):
13
The conjugate transposition of a matrix is denoted by [·]H. Because the different PRN codes are largely orthogonal to each other, cross-correlation components were neglected. The correlation colors the originally white noise ηk and led to 
with covariance matrix (Misra & Enge, 2006, Chapter 10) as indicated in Equation (14):
14
where Equation (15):
15
is the autocorrelation function of the i-th PRN code. The asterisk indicates the complex conjugate. Now that the correlation step has been fully defined, the precise nature of the EKF can be described.
3.2 Kalman Filtering Algorithm
This section explains the multipath rejection function of the EKF in detail. The EKF will be used to replace the DLL and perform the code tracking with active multipath mitigation. To begin, we defined the state vector by Equation (16):
16
with Equation (17):
17
holding all Nst = 4(Lh + 1) unknown parameters. The operators Re (·) and Im (·) are the real and imaginary parts of a complex number, respectively. In contrast to findings described in earlier publications, for example, Iliopoulos et al. (2017); Iltis (1990), the state vector was augmented by the code Doppler 
as described by Siebert et al. (2021b). This facilitated the operation of the proposed EKF without carrier-aiding, if desired. Moreover, the complex channel coefficients 
were divided into real and imaginary parts to maintain the state vector as real-valued, as initially proposed by Siebert et al. (2021a). This was required to prevent the substitution of complex values for real-valued parameters for code delay 
and Doppler 
. To maintain the observability of the system, the channel tap spacing must equal the correlator spacing, i.e., Th = Δτ, and the CIR width must fulfill WCIR ≤ Wbank. We set WCIR = Wbank so that there would be no artificial limitations to the degrees of freedom of the Kalman Filter (KF). It thus followed that Ntap = Ncorr. The process model, which determines the temporal evolution of the state parameters, is set as shown in Equation (18):
18
where Equation (19):
19
is the block diagonal process matrix. 
denotes an identity matrix. The process model predicts the code delay and Doppler according to a constant velocity model. Furthermore, we assumed that the channel coefficients follow, independent of that, a random walk model. The process noise 
was given accordingly as described by (Bar-Shalom et al., 2001) and shown here in Equation (20):
20
The interdependence of CIR and code delay presents an ambiguity for the EKF. A code delay that is too small can be compensated by a CIR where the LOS path is represented in a later channel tap and vice versa. A constraint must be added to the EKF to ensure that the LOS path is always on the central channel tap. Siebert et al. (2021a) proposed an additional constraining measurement as shown in Equation (21):
21
This constraint holds the LOS in place by maximizing the power of the central channel tap 
while minimizing the remaining non-line-of-sight (NLOS) taps. The measurement variance 
determines the tightness of the constraint. The choice of 
represents a trade-off between distorting the CIR estimate in favor of the LOS amplitude and, with a too loose constraint, letting the LOS drift off the central channel tap 
. The latter may cause the EKF to diverge. The constraining measurement is suitable for multipath environments, as the reflected signals are generally weaker than direct ones, given that the satellite has an unobstructed LOS path. Additionally we note that the constraint is also resilient against cycle slips and loss of lock from the PLL. Carrier phase errors reflect as an overall phase offset in the estimated CIR taps 
. Because the constraint considers only the absolute values of the channel taps, phase offsets have no impact on this process. Moreover, sudden carrier phase changes that appear with cycle slips may induce a transient process in which the EKF adjusts its CIR estimate, which might also affect channel tap amplitudes. However, this does not have a lasting effect on code delay estimates. The complete measurement model of the EKF, which includes the newly-introduced constraining measurement as well as the post-correlation signal model from Equation (13), can then be summarized as indicated in Equation (22):
22
with Equation (23):
23
where Nmeas = 2NcorrNsig +1. Similar to the state vector, the complex measurements were also decomposed into their real and imaginary parts. This was necessary to maintain the entire EKF and the state vector xk as real-valued. The function f(xk) has been introduced to permit cleaner notation in the following. The signal model also accounts for the radio propagation channel with the help of a CIR. Because WCIR ≪ Tcode, i /2 applies in general, this is only an approximation; multipaths with delays > WCIR are not considered, see Equation (6). In particular the outer correlators highlight this shortcoming because they sample correlation triangles of multipaths with delays Wbank < τ < Wbank + Tc that still partially extend into the correlator bank, where Tc is the PRN code chip duration. The correlation result obtained cannot be accurately reproduced with the signal model. This modeling mismatch effectively reflects in an increase in the observed measurement noise level of affected correlators. To account for this issue, we propose to scale the measurement noise covariance matrix of the correlator outputs from Equation (14) as shown in Equation (24):
24
where ○ depicts the Hadamard product. The weighting function w is defined as indicated in Equation (25):
25
which relies on the shape of the well-established Tukey window (Prabhu, 2014) as indicated in Equation (26):
26
We set W = Wbank + ΔT to ensure that fTukey(τ) ≠0, ∀τ ∈ [−LcΔτ, LcΔτ]. A Tukey window was chosen because its tuning parameter 0 < α ≤ 1 permits a convenient selection between no scaling, i.e., α = 0, and a Hann window shape with α = 1. The entire measurement noise covariance matrix of the EKF is then given by the following block diagonal matrix as shown in Equation (27):
27
where in Equation (28):
28
Now that the process and measurement model has been fully defined, the EKF can be applied. This estimates the state vector xk based on the measurements z0, …, zk in the minimum mean square error (MMSE) sense. This is done in two steps. First, the state vector of the next time step k is predicted. Then, when the new measurement zk became available, the state prediction is updated. The two steps are explained in the sections to follow.
3.2.1 Prediction Step
In the prediction step, the state vector of the upcoming time step k was predicted based on the previous state estimate 
and the process model A. Because this is a linear process model, one can use the standard KF Equations (29) and (30):
29
30
where P is the covariance matrix of the estimated state vector 
. The indexing notation k|k−1 denotes that the estimate at time step k is based on all measurements up to and including k−1. The prediction of these measurements is then as shown in Equation (31):
31
As the EKF is an iterative filter, an initial state estimate 
with an initial state covariance must be provided as shown in Equation (32):
32
For the initial code delay estimate, this value must approximately match the actual LOS delay.
3.2.2 Update Step
The update step can be performed once the next measurement zk becomes available. It corrects the state prediction 
with the new measurement. Because the measurement model in Equation (22) is non-linear, we needed to use the full EKF equations (Bar-Shalom et al., 2001) as shown in Equations (33) and (34):
33
34
where Kk and Sk are the Kalman gain and the residual covariance matrix, respectively. The measurement model was linearized at the current state estimate 
with the Jacobian matrix shown in Equation (35):
35
The individual partial derivatives can be found in Appendix A. With the assumption that the code delay estimate 
is close to the actual code delay 
the following approximation to the measurement model from Equation (22) can be made as shown in Equation (36):
36
The first transformation is possible because of the circularity of the PRN signals si and accordingly, of the correlator bank Ci. We then approximated 
and 
as equal. The resulting vector of autocorrelation functions in the final expression is time-invariant and thus, can be pre-computed in an initialization phase. This will reduce the load during run-time.
3.3 Implementation Details
Implementation-relevant aspects of the proposed algorithm are discussed in the sections to follow.
3.3.1 Receiver Architecture
Figure 1 is a schematic illustration of the proposed multipath mitigation algorithm integrated into a GNSS receiver architecture. All required components needed to track all Nsig signals from a single satellite and provide the tracking results to the PVT solution are included. The RF signal from the front-end is introduced into the architecture from the left and is first down-converted with a mixer. The carrier replica required for this process is generated using a numerically-controlled oscillator (NCO). The baseband signal yk, which was introduced in Equation (8), is then passed to the correlator banks Ci from Equation (12), i.e., one for each signal component of the tracked satellite. Correspondingly, Nsig correlator output vectors 
are obtained, as defined in Equation (13). The data bits are estimated and eliminated from the correlator outputs for the code and carrier tracking. This is possible at sufficiently high carrier-to-noise density ratios (C/N0s), which is assumed to be the case in our application. The EKF processes all correlator outputs. The estimated code delay and Doppler 
and 
, respectively, are then used to drive the correlator banks, ideally free of multipath errors. Furthermore, 
will be required for the PVT solution. The carrier phase tracking was implemented with a PLL which relies on all central correlator outputs, i.e., the prompts. The individual prompts were summed up coherently for the PLL to reduce noise after the data bits were eliminated. If the navigation signals are broadcast in different phase positions, all prompts would need to be rotated accordingly to match the primary signal i = 1.
Proposed receiver architecture with multipath estimating and mitigating EKF
3.3.2 Local PRN Replicas
The choice of the local replicas for the measurement model is very critical when using MEDLL-like techniques. These replicas must precisely reflect all signal distortions the signal experiences, with the exception of multipath propagation, including transmitter chain imperfections of the satellite and RF filter characteristics of the receiver. Discrepancies might be misinterpreted as multipath by the algorithms and thus may potentially lead to biases. In van Nee et al. (1994), for example, the correlation result obtained in a low multipath environment was averaged over time to obtain an authentic reference correlation function for the measurement model with a negligible noise level. This is not necessary for our proposed solution. Instead, deviations between the local replica and the signal shape that was actually received will be absorbed by the estimated CIR. This feature also allows to adjust to time-varying characteristics. The effect of these deviations on tracking performance can be neglected for most applications. Therefore, we simply built the local replicas si(t) for the correlator banks and the EKF measurement model analytically from the ground up. A pulse shape with a finite spectrum was used to compose the PRN signal out of the superposition of adjacent sampled pulses. As band-limited pulses are infinite in time, contributions to neighboring pulses must be neglected inevitably at some point. For our model, we assumed that the contribution can be neglected after ±10Tc. The residual modeling error was insignificant.
3.3.3 Effect of Doppler
The relative velocity between the satellite and the user introduces a Doppler shift to the received signal. The signal model in Section 2 and correspondingly, the algorithm described in Section 3 focused on improved readability and did not consider this effect. Nevertheless, the local PRN replicas would need to undergo constant shortening or stretching according to the current Doppler frequency, not only for the correlator banks Ci but also for the local replicas in the EKF measurement model as described in Equation (22). This would eliminate the opportunity to pre-compute the autocorrelation functions for the measurement prediction and the Jacobian as shown in Equations (36) and (35), respectively. However, for most applications, the observed Doppler frequency is in the range of ±10 kHz. The resulting shortening or stretching of si(t) over one code period Tcode,i will be found in the region of Tc/100 for the GPS C/A code. Correspondingly, the effect on the autocorrelation functions Φsisi and the correlation results 
can be neglected. Therefore, we chose to use local replicas with zero Doppler at all times so that we would not need to make adjustments to maintain a lower level of computational complexity.
4 PERFORMANCE EVALUATION
The performance of the proposed multipath mitigating algorithm was demonstrated in the following three steps. First, a simulation with synthetic signals was performed to evaluate the basic functionality of the system. This was followed by a hardware emulation using a GNSS constellation simulator to test the solution in a more complex scenario. Finally, a measurement campaign was carried out to evaluate its real world performance. In all cases, the IQ rate was set to fs = 20 MHz. An NI PXIe-51740 4-channel oscilloscope card was used as a digitizer for the hardware emulation and measurement campaign. Figure 2 is a schematic that illustrates the different setups used in this scenario. The efficacy of the proposed solution was compared to that of a conventional state-of-the-art GNSS receiver. The non-coherent DLL has a second-order loop filter with a noise bandwidth of 0.5 Hz and an early-late spacing of 0.1Tc. For the PLL, the correlator outputs were integrated over 4 ms and subsequently pass through a two-quadrant ATAN Costas discriminator; this was followed by a third-order loop filter with 9 Hz bandwidth. When multiple signal components per satellite are available, as is the case, for example, with the Galileo OS with its data and pilot components, the DLL and PLL sum up the correlation results to reduce noise. Moreover, the correlator outputs were integrated over Tint = 20 ms for the DLL and the proposed solution. As mentioned earlier, data bits were eliminated before integration. Furthermore, the EKF uses a correlator and channel tap spacing of Δτ = Th = 0.05Tc with the bank and CIR width set to Wbank = WCIR = Tc, resulting in Nst = 84 states. Unless otherwise specified, the tuning parameter for the scaling of the correlator noise covariances was set to α = 1, see Subsection 3.2. To obtain accurate initialization values for the EKF, we rely on the DLL until the tracking reaches a steady state before switching over. The initial code delay and Doppler 
and 
, respectively, are then set according to the latest DLL tracking results. Moreover, for the initial CIR, we set the LOS tap to 
and the NLOS taps to 
, l ≠ 0. The standard deviation of the constraint was set to σconstr = 10−3. The measurement noise variance 
was determined by correlating the received baseband signal yk with the PRN code of a non-existing satellite. The variance of the obtained correlator output corresponds to the diagonal elements of Rηci. Thus, 
can be determined. Additional EKF parameters are presented in Table 1.
Schematic illustrations of the different setups (a) Synthetic signals (b) Hardware emulation (c) Measurement campaign
EKF Parameters for Performance Evaluation
4.1 Synthetic Signals
The basic functionality of the algorithm was demonstrated using synthetic signals. GNSS signals with a one-sided bandwidth of B = 10 MHz were generated and stored for subsequent offline processing, as illustrated in the block diagram in Figure 2(a). In all scenarios, a single geostationary GPS satellite was considered that broadcasted the legacy C/A code and reached the user with 45 dB Hz C/N0. The received signals have zero Doppler and were generated in the same way as the local replicas as described in Subsection 3.3.2.
As a first step, we showed that the proposed technique and the conventional DLL exhibited comparable convergence times. To demonstrate this point, we simulated a sudden 5 m change in the received pseudorange. The results of this simulation are shown in Figure 3(a). The plot of the proposed solution begins a few seconds after the one that represents the DLL results. This is because the multipath mitigating EKF is only enabled after initial conventional code tracking with a DLL. This initial tracking phase was not plotted because it would be identical to that of the conventional DLL-based code tracking (shown in red). Once a steady tracking state is reached, the proposed algorithm switches over to EKF-based tracking. After 15 s, a 5 m jump occurred. Both code-tracking algorithms show very similar step responses. We note that the convergence behavior of the proposed solution depends on the selection of the process as well as on measurement noise covariance matrices; this is similar to the DLL, in which different step responses can be achieved by adjusting the loop filter. We then explored the response time to a change in the propagation path. A multipath with a 50 m path delay and a 3 dB LOS-to-multipath power ratio appeared after 15 s seconds of LOS conditions (Figure 3(b)). The DLL responds to this multipath as expected, and reached a steady-state tracking error of 10 m after approximately 10 s. By contrast, the proposed solution with its multipath mitigation settled at an average error of 0.5 m and adjusted to the new conditions in under 3 s.
Step response and convergence time of proposed solution compared against a conventional DLL (a) Step response: at 15 s, the pseudorange of the satellite performs a unit jump of 5 m. (b) Convergence time: at 15 s, a multipath with 50 m path delay and 3 dB LOS to multipath to power ratio appears.
The multipath error envelopes of the proposed solution were then determined and compared to the envelopes of a DLL. A LOS signal was superimposed with an additional delayed, attenuated, and phase-shifted replica that simulated a two-path multipath environment. In contrast to the analysis above, the GNSS signals were generated completely free of noise, and we are now also considering the Galileo E1 OS. The multipath delays were increased stepwise from 0 to 650 m. The envelopes show the range of code tracking errors obtained over an entire multipath phase rotation. Each delay and phase combination was maintained until a steady tracking state was reached, with a LOS-to-multipath power ratio set at 3 dB. Figure 4(a) and Figure 4(b) present the multipath error envelopes for the GPS C/A code and for the Galileo E1 OS C ranging code with a simplified BOC(1,1) modulation, respectively. The superposition of the three different phenomena led to the multipath error envelope shapes obtained for the proposed multipath mitigation algorithm. These are discussed in detail in the following sections.
Multipath error envelopes. The signals considered in this analysis were limited to a one-sided bandwidth of B = 10 MHz with a LOS to multipath power ratio of 3 dB. For the DLL, the early-late spacing was set to 0.1Tc. Setting the tuning parameter to α = 1 accounts for the approximation made in the signal model, see Subsection 3.2. (a) GPS C/A code (b) Galileo E1 OS C with BOC(1,1) modulation
First, as expected and similar to the DLL, the multipath error increased with increasing multipath delay up to errors of 5.4 m. For delays between 5 and 15 m, depending on the multipath phase, the algorithm begins to differentiate between the LOS and the multipath. The error then begins to decrease and reaches zero once again for delays of approximately 24 m. From this point forward, an oscillating error was observed. This fluctuation is easily explained by the shape of the estimated CIR, see Figure 5. Side lobes arise whenever a multipath does not align perfectly with one of the CIR taps. This can be seen, for example, in Figure 5(a) for multipath delays <300 m with small ripples between the main LOS and multipath signal components. These secondary peaks were superimposed on the LOS. As a result, the global maximum of the CIR deviated slightly from the LOS delay, alternately to the left and the right. Correspondingly, the constraining measurement ensured that the code delay was adjusted until the central channel tap held the maximal power once again, leading to oscillating errors. This effect can be observed equally for both signals under consideration.
CIRs for different multipath delays. In-phase multipath sweeps from 0 to 650 m with a 3 dB LOS-to-multipath power ratio. Only the real parts of the CIRs were plotted because the imaginary parts are zero in this in-phase multipath. (a) GPS C/A, α = 0 (b) GPS C/A, α = 1 (c) Galileo E1 C BOC(1,1), α = 0 (d) Galileo E1 C BOC(1,1), α = 1
We also identified an increasing error for both GPS and Galileo signals with delays of approximately 
and 
, respectively. This was because of the constraining measurement that distorted the amplitudes of the channel taps in favor of the central LOS tap 
. Correspondingly, the estimated LOS amplitude tended to be too high and the multipath amplitudes too low. While this satisfies the constraint, it was at the cost of the other optimization criteria, including matching of the correlator outputs. In an attempt to achieve a better match of the correlator outputs, side peaks arose in the CIR. This effect was clearly observed in the estimated CIRs shown in Figure 5. The LOS and the sweeping multipath were followed by side peaks every Tc for the GPS C/A code. For the Galileo signal with the BOC(1,1) modulation these peaks occurred every Tc/2 because of the doubled chip rate after modulation. The side peaks of the multipath were superimposed with the LOS; this led to a shift of the global maximum in the CIR and thus to tracking errors. This effect becomes even more pronounced with tighter constraints. However, constraints that are too loose may result in divergence because the LOS is not held reliably in place.
A third phenomenon was observed for multipaths with delays between Tc and 2Tc and no correlator output covariance scaling, i.e., α = 0. These multipaths lay outside the correlator bank width Wbank but their correlation triangle still extends into the correlation result. Correspondingly, this cannot be properly represented by the CIR. However, the EKF attempts to approximate this situation with available CIR taps, thus leading to additional side peaks. These peaks are also superimposed on the LOS and lead to tracking errors, see Figures 5(a) and 5(c). These phenomena can be resolved by extending the correlator bank and CIR width to Tcode,i/2, which will ensure that the observed correlator outputs remain periodic. However, this would add an enormous computational complexity. Alternatively, one can apply a scaling of the correlator output covariance matrix to account for the model mismatch, as we proposed in Subsection 3.2. When setting α = 1, a Hann window-shaped scaling is used that largely resolves this problem. The effect of scaling can also be observed in Figures 5(b) and 5(d), where the side peaks from 300 m onward are largely suppressed.
4.2 Hardware Emulation
A hardware emulation was performed to examine the proposed method in a more sophisticated environment at the German Aerospace Center (DLR) multi-output advanced signal test environment for receivers (MASTER). This facility has several GSS9000 series Spirent GNSS constellation simulators (Konovaltsev et al., 2019) which have been used to generate authentic GNSS signals of a full GPS constellation in real time and to broadcast the legacy C/A code and the new L1C signals. The signals were then digitized and stored for subsequent offline processing, as illustrated in Figure 2(b). The constellation simulator provides a controlled environment with known ground truths. With the user position set to 48°N and 11°E at an altitude of 550 m, 9 satellites were in sight; 4 of these satellites were selected for the processing. In the simulated scenario, 3 of the 4 satellites were disturbed after approximately 53 s, each by a single static multipath. Table 2 presents more details on this finding.
Satellite Parameters of Hardware Emulation
The resulting position root-mean-square errors (RMSEs) were compared to one another as shown in Figure 6(a). This representation shows the results with a conventional DLL and the proposed solution. Each was evaluated twice, once using only the C/A code and once with the joint tracking of C/A and L1C. We note that the two cases that included the L1C achieved a first PVT solution earlier than those that used the C/A alone. This is because, in this particular example, the information required for the PVT solution was decoded earlier by the receiver from L1C navigation messages than from those that used the C/A alone. In the absence of multipaths, RMSEs below 2 m were obtained in all cases; errors developed once the multipaths began to occur. Using the DLL, the RMSE settled at 13 m with the use of the C/A code. Adding L1C reduced the errors to 8.1 m because of the superior multipath error envelope of the new signals (Hein et al., 2006). This new approach achieves a steady-state error on average of 0.6 m with single-signal tracking. Multi-signal tracking does not lead to reductions in the absolute error, but does successfully reduce the estimate variance because of the higher overall signal power, similar to what was observed for the DLL. This can also be observed in Figure 6(b), in which the standard deviation of the position estimates are presented. As shown, adding L1C to the code tracking reduces the estimate variance. Note, in this representation, values between 50 and 55 s were omitted to exclude the jump that occurs in response to the multipath that arises.
Results of hardware emulation (a) Position RMSE (b) Standard deviation of position RMSE
4.3 Measurements
A measurement campaign was conducted to evaluate the performance of the proposed EKF-based multipath estimation and mitigation solution in a real world environment. In this experiment, a measurement vehicle drove along a route with different multipath environments while an Antcom G8 antenna on the roof of the car received GNSS signals. The antenna signal was divided so that it provided information to the NovAtel PwrPak7 reference receiver and the digitizer to generate a reference track and to store the raw data, respectively, as illustrated in Figure 2(c). The system relies on inertial measurement units (IMUs) that help to bridge short GNSS outages. The log files were processed by NovAtel’s software, GrafNav, to determine afterwards position estimates with forward/backward post-processing taking into account precise point positioning (PPP) correction data. The car was driven from the Oberpfaffenhofen DLR site near Munich to the neighboring town Germering and back. The route included highways, forests, open fields, and suburban environments as shown in Figure 7(a). Most sections of the 42 min run featured minimal multipath propagation. Stronger multipath was observed in the suburban environment. The satellite constellation is shown in Figure 7(b). The GPS C/A code and the Galileo E1 OS with data and pilot components were used for processing with an elevation mask of 10°. Position estimates were rejected if the sum of the squared residual pseudoranges exceeded 3 m, and the estimated time scale offset between GPS and Galileo deviated from the broadcast Galileo-GPS Time Offset (GGTO) by more than 20 ns or the geometric dilution of precision (GDOP) was greater than 20. Moreover, the proposed approach and the conventional solution, both rely on a snapshot-based least-squared algorithm for the PVT solution. Thus, no user motion models were involved in this evaluation.
Route and satellite constellation of measurement campaign (a) Route of the measurement campaign projected onto a satellite image from Google Earth: the arrow indicates the direction of travel. Start and endpoint were at the DLR site. (b) Satellite constellation: the integers correspond to the respective PRN code numbers. Crosses and circles indicate the final positions of GPS and Galileo satellites, respectively.
Figure 8(a) documents the position RMSE of the entire measurement campaign over time. As indicated in the plot, we observed increased multipath propagation most notably in the suburban environment. Here, the differences between the conventional DLL and the proposed solution became apparent. On route sections with less severe multipath propagation, for example, while driving on highways or in open fields, both solutions perform for the most part identically with smoother estimates presented occasionally by the new approach. Altitude estimates are shown in Figure 8(b). Here, we detected a constant offset between the ground truth and estimates of approximately 3 m. A comparison of the slant total electron contents (STECs) from the broadcast Klobuchar model coefficients as well as from ionospheric maps revealed deviations of up to 25 TECU per satellite. This suggests that the altitude bias may be caused by erroneous ionospheric corrections that were generated by the Klobuchar model. After approximately 24 min, a larger PVT outage was observed. This was because the measurement vehicle stopped at an unfortunate location, in which buildings on both sides of the street blocked the sky view and caused most of the low-elevation satellites to disappear. However, the GPS satellites 2 and 4 continued to be tracked. As the direct LOS was no longer accessible, the code tracking became locked in with the reflected signal. Accordingly, the resulting pseudorange estimates were significantly in error. This led to time scale offsets that exceeded the specified threshold of 20 ns and thus a rejection of the position estimate. Both the conventional DLL and the proposed multipath rejecting EKF had difficulties and were not prepared for this situation. The constraint assumes that the LOS will always be the strongest signal path. When it disappears, the EKF adjusts to accommodate the second strongest path. This drawback could be addressed, for example, by using receiver autonomous integrity monitoring (RAIM) (van Graas & Farrell, 1993), which is an algorithm that monitors the consistency of the available pseudoranges. Under appropriate conditions, this facilitates the detection and exclusion of satellites that are only received only via a NLOS path. Alternatively one could add C/N0 (Groves & Jiang, 2013) or elevation-based weighting (Jin & de Jong, 1996) to the PVT solution or rely directly on vector tracking loops Siebert et al. (2022). Both would result in improved positioning by relying more on signals from the remaining satellites. Another PVT outage occurred after 33 min when the measurement vehicle was traveling across the second highway section. This outage occurred because of a tunnel that resulted in a complete loss of signal. The fusion with other sensors, such as IMUs or, for vehicles, wheel speed sensors, would be needed to overcome the blocked GNSS.
Results of the measurement campaign. The color-coded time axes indicate the different route sections, as introduced in Figure 7(a). (a) RMSE of position estimates (b) Altitude estimates
One section of the suburban route was examined in detail, as shown in Figure 9. Here, residential buildings result in multipaths from alternating sides of the street. As a result of the additional signal replicas, the position estimates of the DLL drift from the ground truth and include altitude errors of up to 15 m. Furthermore, the DLL revealed repeated PVT outages. This was due to pseudorange estimate biases generated by strong multipaths. The resulting PVT solutions have residual errors that exceed the specified threshold that were ultimately discarded. At the same time, the proposed solution suppressed multipath effects and managed to keep its track closer to the ground truth with minimal PVT outages most of the time. However, large outliers appeared at approximately 22.2 min when the measurement vehicle was driving around a corner. These outliers were caused by the Galileo satellite 31. While this satellite disappeared at this point because of its low elevation and its position behind a building, its signal continued to be received via its reflection from a building on the other side of the street until its eventual loss. As discussed earlier for the PVT outage after 24 min, the EKF itself was not designed for NLOS-only reception and must be combined with other approaches to overcome this gap.
Close up of position estimates and altitude errors from 21.24 to 22.57 min in the suburban route section. The arrow indicates the direction of travel. Residential buildings cause multipaths from alternating sides of the street that were mostly mitigated by the proposed solution, while the DLL deviates repeatedly from the ground truth and looses its PVT solution several times. Errors around min 22.2 are caused by NLOS-only reception. (a) Position estimates (b) Altitude errors
5 CONCLUSION
In this paper, a multipath mitigation algorithm was proposed. This algorithm was based on a multi-correlator approach using an EKF to replace the conventional DLL and its early and late correlators. Because the measurement model of the EKF incorporates a radio propagation channel between the satellite and the user, multipath effects are naturally taken into account. This results in a more robust code tracking, as well as an estimate of the CIR. In contrast to other multipath estimating and mitigating approaches, for example, MEDLL, precise knowledge of the shape of the received signal is less crucial for the calculation. In this algorithm, the estimated CIR will absorb deviations from the actual signal shape, for example, those resulting from unknown receiver front-end filter characteristics. Moreover, the proposed solution is of only moderate complexity due to the applied KF-based filtering.
Synthetically-generated GNSS signals were used to determine multipath error envelopes. The results revealed the multipath mitigation characteristics of the proposed EKF-based solution. Step responses were also used to demonstrate that the proposed solution had a similar convergence time as that provided by conventional DLL-based code tracking. Hardware emulation using a GNSS constellation simulator was performed to determine the positioning performance in a model multipath scenario that included both for single- and multi-signal tracking. The multipath mitigation capabilities of the proposed solution were demonstrated. Subsequently, a measurement campaign was conducted that provided an extensive comparison between the results from the proposed solution and a conventional receiver under real world conditions. For this purpose, a measurement vehicle was driven along a diversified route that included suburban sections and highways as well as short forest passages. The proposed algorithm outperformed a conventional DLL-based GNSS receiver in these settings. As expected, the differences were particularly large in environments with stronger multipaths. Moreover, the limitations of the proposed solution were explored.
Future work in this field might focus on simplifying the presented approach. For example, one might attempt to increase the correlator spacing and determine the performance of this algorithm in response to sub-sampling-induced information loss.
HOW TO CITE THIS ARTICLE
Siebert, C., Konovaltsev, A., & Meurer, M. (2023). Development and validation of a multipath mitigation technique using multi-correlator structures. NAVIGATION, 70(4). https://doi.org/10.33012/navi.609
A JACOBIAN OF MEASUREMENT MODEL
In the following, the individual partial derivatives of the Jacobian matrix Jf,k from Equation (35) of the EKF measurement model were derived under consideration of the approximation made in Equation (36). We begin with the code delay. The derivative with respect to 
is given as shown in Equation (A1):
A1
where 
(τ). Furthermore, the partial derivatives with respect to the channel coefficients can be expressed as Equation (A2):
A2
For the constraining measurement, the following derivations exist, as indicated in Equations (A3) and (A4):
A3
A4
All other partial derivatives not explicitly stated above are zero. Similar to the pre-computation of the autocorrelation functions for the measurement prediction, one can also reduce the load during run-time by pre-computing the time-invariant components of Equations (A1) and (A2).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.