Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments

Abstract

A sequential Bayesian estimation algorithm for multipath mitigation is presented, with an underlying movement model that is especially designed for dynamic channel scenarios. In order to facilitate efficient integration into receiver tracking loops, it builds upon complexity reduction concepts that previously have been applied within maximum likelihood (ML) estimators. To demonstrate its capabilities under different GNSS signal conditions, simulation results are presented for both BPSK-modulated and BOC-(1,1) modulated navigation signals.

1. Introduction

A major error source within global navigation satellite systems (GNSSs) comes from multipath, the reception of additional signal replica due to reflections, which introduce a bias into the estimate of the delay lock loop (DLL) of a conventional navigation receiver. For efficient removal of this bias, it is possible to formulate advanced maximum likelihood (ML) estimators that incorporate the echoes into the signal model and are capable of achieving the theoretical limits given by the Cramér Rao bound. The drawback of ML estimator techniques is that the parameters are assumed to be constant during the time of observation. Independent estimates are obtained for successive observation intervals, whose length has to be adapted to the dynamics of the channel.

In this paper, we consider the important practical case of dynamic channel scenarios and assess how the time-delay estimation can be improved if information is available about the temporal evolution of the channel parameters. Our approach is based on Bayesian filtering, the optimal and well-known framework to address such dynamic state estimation problems [1]. Sequential Monte Carlo (SMC) methods [2, 3] are used for computing the posterior probability density functions (PDFs) of the signal parameters.

2. A Comparison of Various Multipath Mitigation Approaches

Figure 1 gives an overview of the relationships between different multipath mitigation and estimation approaches. In fact, we have chosen to discriminate approaches according to their primary objective. The left column represents the class of techniques that attempt to mitigate the effect of multipath in different ways. This can for example be achieved by modifications of the antenna response, either by means of hardware design or with signal processing techniques (e.g., beamforming). The majority of the remaining mitigation techniques are in some way aligning the more or less traditional receiver components (e.g., the early/late correlator) to the signal received in the multipath environment. The probably most simple technique is the adjustment of the correlator spacing applied in the Narrow Correlator [4]. Other well-known examples of this category are the Strobe Correlator [5], the Gated Correlator [6], or the Pulse Aperture Correlator [7]. To incorporate new signal forms (such as BOC), these methods need “tuning” in order to suffer as little as possible from multipath. On the other hand, multipath estimation techniques (right column) treat multipath (in particular the delay of the paths) as something to be estimated from the channel observations, so that its effects can be trivially removed at a later processing stage. For the estimation techniques, we have differentiated between static and dynamic approaches, according to the underlying assumption of the channel dynamics. Examples for static multipath estimation are those belonging to the family of ML estimators, often using different efficient maximization strategies over the likelihood function [8–13]. For static channels without availability of prior information, the ML approach is optimal and performs significantly better than other techniques, especially if the echoes have short delay. An estimator based on sequential importance sampling (SIS) methods (particle filtering) for static multipath scenarios has been considered in [14], which has the advantage that prior channel knowledge can be incorporated.

Figure 1 

Classification of multipath mitigation approaches.

As a first step towards addressing dynamic channels, one can incorporate ML estimators in receiver loops or formulate quasisequential estimators [15, 16]. Finally, dynamic estimators that target the computation of the posterior PDF conditioned on the received channel output sequence at the receiver can be applied on a per single range basis or operate in the position domain. In this paper, we concentrate on dynamic estimators applied per each range. The sequential Monte Carlo approach has also been suggested in the communications field for estimation of time-varying channel responses in spread spectrum systems [17, 18].

3. Signal Model

Assume that the complex valued baseband-equivalent received signal is equal to𝑧(𝑡)=^𝑁𝑚_𝑖=1𝑒_𝑖(𝑡)⋅𝑎_𝑖(𝑡)⋅[𝑐(𝑡)∗𝑔(𝑡−𝜏_𝑖(𝑡))]+𝑛(𝑡), (1) where 𝑐(𝑡) is a delta-train code sequence that is modulated on a pulse 𝑔(𝑡), 𝑁_𝑚 is the maximum number of allowed paths reaching the receiver (to restrict the modeling complexity), 𝑒_𝑖(𝑡) is a binary function that controls the activity of the 𝑖^′th path, and 𝑎_𝑖(𝑡) and 𝜏_𝑖(𝑡) are their individual complex amplitudes and time delays, respectively. The signal is disturbed by additive white Gaussian noise 𝑛(𝑡). Grouping blocks of 𝐿 samples at times (𝑚+𝑘𝐿)𝑇_𝑠, 𝑚=0,…,𝐿−1, together into vectors 𝐳_𝑘, 𝑘=0,1,…,𝑛 whilst assuming the parameter functions 𝑒_𝑖(𝑡), 𝑎_𝑖(𝑡), and 𝜏_𝑖(𝑡) being constant within the corresponding time interval and equal to 𝑒_𝑖,𝑘, 𝑎_𝑖,𝑘, and 𝜏_𝑖,𝑘, this can be rewritten as𝐳_𝑘=𝐂𝐆(𝝉_𝑘)𝐄_𝑘𝐚_𝑘+𝐧_𝑘=𝐬_𝑘+𝐧_𝑘. (2)In the compact form on the right hand side, the samples of the delayed pulses 𝐠(𝜏_𝑖,𝑘) are stacked together as columns of the matrix 𝐆(𝝉_𝑘)=[𝐠(𝜏_1,𝑘),…,𝐠(𝜏_{𝑁𝑚,𝑘})], 𝐂 is a matrix representing the convolution with the code, and the delays and amplitudes are collected in the vectors 𝝉_𝑘=[𝜏_1,𝑘,…,𝜏_{𝑁𝑚,𝑘}]^𝑇 and 𝐚_𝑘=[𝑎_1,𝑘,…,𝑎_{𝑁𝑚,𝑘}]^𝑇, respectively. Furthermore, for concise notation we use 𝐄_𝑘=diag[𝐞_𝑘] whilst the elements of the vector 𝐞_𝑘=[𝑒_1,𝑘,…,𝑒_{𝑁𝑚,𝑘}]^𝑇, 𝑒_𝑖,𝑘∈[0,1], determine whether the 𝑖^′th path is active or not by being either 𝑒_𝑖,𝑘=1 corresponding to an active path or 𝑒_𝑖,𝑘=0 for a path that is currently not active. The term 𝐬_𝑘 denotes the signal hypothesis and is completely determined by the channel parameters 𝝉_𝑘, 𝐚_𝑘, and 𝐞_𝑘. Using (2), we can write the associated likelihood function as𝑝(𝐳_𝑘∣𝐬_𝑘1)=(2𝜋)^𝐿𝜎^2𝐿1⋅exp[−2𝜎²(𝐳_𝑘−𝐬_𝑘)^𝐻(𝐳_𝑘−𝐬_𝑘)]. (3)The likelihood function will play a central role in the algorithms discussed in this paper; its purpose is to quantify the conditional probability of the received signal conditioned on the unknown signal (specifically the channel parameters).

3.1. Efficient Likelihood Computation

In [11], a general concept for the efficient representation of the likelihood (3) was presented, which is applicable to many of the existing ML multipath mitigation methods. The key idea of this concept is to formulate (3) through a vector 𝐳_𝑐,𝑘 resulting from an orthonormal projection of the observed signal 𝐳_𝑘 onto a smaller vector space, so that 𝐳_𝑐,𝑘 is a sufficient statistic according to the Neyman-Fisher factorization [19] and hence suitable for estimating 𝐬_𝑘. In other words, the reduced signal comprises the same information as the original signal itself. In practice, this concept becomes relevant as the projection can be achieved by processing the received signal (2) with a bank of correlators and a subsequent decorrelation of the correlator outputs. A variant of this very general concept, applied in [13], has also been referred to as the Signal Compression Theorem in [20] for a set of special projections that do not require the step of decorrelation due to the structure of the used correlators. For instance, unlike the correlation technique used in [8], the one suggested in [13] already projects onto an orthogonal and thus uncorrelated subspace, similar to the code matched correlator technique proposed in [11]. Due to complexity reasons, all practically relevant realizations of ML estimators [8, 13] operate in a projected space, namely after correlation. The corresponding mathematical background will be discussed below, including also interpolation of the likelihood and elimination of complex amplitudes as further methods for complexity reduction.

3.1.1. Data Compression

As explained above, the large vector containing the received signal samples 𝐳_𝑘 is linearly transformed into a vector 𝐳_𝑐,𝑘 of much smaller size. Following this approach, the likelihood according to (2) can be rewritten as𝑝(𝐳_𝑘∣𝐬_𝑘1)=(2𝜋)^𝐿𝜎^2𝐿𝐳exp[−^𝐻𝑘𝐳_𝑘2𝜎²]⋅exp[ℜ{𝐳^𝐻_𝑘𝐐_𝑐𝐐^𝐻_𝑐𝐬_𝑘}𝜎²−𝐬^𝐻_𝑘𝐐_𝑐𝐐^𝐻_𝑐𝐬_𝑘2𝜎²]=1(2𝜋)^𝐿𝜎^2𝐿𝐳exp[−^𝐻𝑘𝐳_𝑘2𝜎²]⋅exp[ℜ{𝐳^𝐻_𝑐,𝑘𝐬_𝑐,𝑘}𝜎²−𝐬^𝐻_𝑐,𝑘𝐬_𝑐,𝑘2𝜎²], (4)with the compressed received vector 𝐳_𝑐,𝑘 and the compressed signal hypothesis 𝐬_𝑐,𝑘:𝐳_𝑐,𝑘=𝐐^𝐻_𝑐𝐳_𝑘,𝐬_𝑐,𝑘=𝐐^𝐻_𝑐𝐬_𝑘, (5)and the orthonormal compression matrix 𝐐_𝑐, which needs to fulfill𝐐_𝑐𝐐^𝐻_𝑐≈𝐈,𝐐^𝐻_𝑐𝐐_𝑐≈𝐈, (6)to minimize the compression loss. According to [11], the compression can be two-fold so that we can factorize𝐐_𝑐=𝐐_𝑐𝑐𝐐_𝑝𝑐 (7)into a canonical component decomposition, given by an 𝐿×𝑁_𝑐𝑐 matrix 𝐐_𝑐𝑐, and a principal component decomposition, given by an 𝑁_𝑐𝑐×𝑁_𝑝𝑐 matrix 𝐐_𝑝𝑐. In [11], two choices for 𝐐_𝑐𝑐 are proposed𝐐_𝑐𝑐=𝐂𝐆(𝝉^𝑏)𝐑⁻¹_𝑐𝑐Signalmatched,𝐂(𝝉^𝑏)𝐑⁻¹_𝑐𝑐Codematched, (8)where the elements of the vector 𝝉^𝑏 define the positions of the individual correlators. The noise-free outputs of the corresponding correlator banks are illustrated in Figure 2. To decorrelate the bank outputs (𝐂𝐆(𝝉^𝑏))^𝐻𝐲 and 𝐂(𝝉^𝑏)^𝐻𝐲, as mentioned above, the whitening matrix 𝐑_𝑐𝑐 can be obtained from a QR decomposition of 𝐂𝐆(𝝉^𝑏) and 𝐂(𝝉^𝑏), respectively. Apart from practical implementation issues, both correlation methods given by (8) are equivalent from a conceptual point of view. For details on the compression through 𝐐_𝑝𝑐, the reader is referred to [11].

(a) BPSK signal matched

(b) BPSK code matched

(c) BOC(1,1) signal matched

(d) BOC(1,1) code matched

Figure 2 

Output of the canonical component type correlator banks $CG({\mathit{\tau }}^b)$ and $C({\mathit{\tau }}^b)$ for BPSK and BOC(1,1).

3.1.2. Interpolation

In order to compute (4) independently of the sampling grid, advantage can be made of interpolation techniques. Using the discrete Fourier transformation (DFT), with 𝚿 being the DFT matrix and 𝚿⁻¹ being its inverse counterpart (IDFT), we get𝐬_𝑐,𝑘=𝐐^𝐻_𝑐𝐂𝚿⁻¹diag[𝚿𝐠(0)]𝛀(𝝉_𝑘)𝐄_𝑘𝐚_𝑘=𝐌_𝑠𝑐𝛀(𝝉_𝑘)𝐄_𝑘𝐚_𝑘, (9)with 𝛀(𝝉_𝑘) being a matrix of column-wise stacked vectors with Vandermonde structure [10, 11], such that the element at row 𝑝 and column 𝑞 computes withℜ{[𝛀(𝝉_𝑘)]_𝑝,𝑞}=cos(2𝜋(𝑝−1)𝜏_𝑞,𝑘𝑁_𝑔𝑇_𝑠),ℑ{[𝛀(𝝉_𝑘)]_𝑝,𝑞}=−sin(2𝜋(𝑝−1)𝜏_𝑞,𝑘𝑁_𝑔𝑇_𝑠). (10)𝑁_𝑔 is the length of the pulse 𝐠 in samples. The advantage of the interpolation is that it can take place in the reduced space. The most costly computations in (9) can be carried out in precalculations as the matrix 𝐌_𝑠𝑐, whose row dimension corresponds to the dimension of the reduced space and whose column dimension is 𝑁_𝑔, is constant.

3.1.3. Amplitude Elimination

In a further step, we reduce the number of parameters by optimizing (4) for a given set of 𝝉_𝑘 and 𝐞_𝑘 with respect to the complex amplitudes 𝐚_𝑘, which can be achieved through a closed-form solution. Using𝐒_𝑐,𝑘=𝐌_𝑠𝑐𝛀(𝝉_𝑘)𝐄_𝑘 (11)and obtaining 𝐒⁺_𝑐,𝑘 by removing zero columns from 𝐒_𝑐,𝑘, one yields the corresponding amplitude values of the active paths:𝐚⁺_𝑘=(𝐒^+𝐻_𝑐,𝑘𝐒⁺_𝑐,𝑘)⁻¹𝐒^+𝐻_𝑐,𝑘𝐳_𝑐,𝑘. (12)As we have introduced a potential source of performance loss by eliminating the amplitudes and thus practically are disregarding their temporal correlation, we propose to optimize (4) using𝐳_𝑐,𝑘=1𝑄⋅^𝑄−1_𝑙=0𝐳_{𝑐,𝑘−𝑙,} (13)with the adjustable averaging coefficient 𝑄. When evaluating (4), we substitute 𝐬_𝑐,𝑘 by^0𝑥0𝐬^02𝑑𝑐_𝑐,𝑘=𝐒_𝑐,𝑘𝐚_𝑘, (14)where the elements of the vector 𝐚_𝑘 that are indicated to have an active path (𝑎_𝑘,𝑖∶𝑖→𝑒_𝑘,𝑖=1) are set equal to the corresponding elements of 𝐚⁺_𝑘. All other elements (𝑎_𝑘,𝑖∶𝑖→𝑒_𝑘,𝑖=0) can be set arbitrarily as their influence is masked by the zero elements of 𝐞_𝑘.

3.2. Review of the ML Concept

The concept of ML multipath estimation has drawn substantial research interest since the first approach was proposed in [8]. Despite being treated differently in various publications, the objective is the same for all ML approaches, namely to find the signal parameters that maximize (3) for a given observation 𝐳_𝑘:𝐬_𝑘̂=argmax_𝐬𝑘{𝑝(𝐳_𝑘∣𝐬_𝑘)}. (15)The signal parameters are thereby assumed as being constant throughout the observation period 𝑘. Different maximization strategies exist, which basically characterize the different approaches. Despite offering great advantages for theoretical analysis, the practical advantage of the generic ML concept is questionable due to a number of serious drawbacks.(i) The ML estimator assumes that the channel is static for the observation period and is not able to exploit its temporal correlation throughout the sequence 𝑘=1,…,𝑛. Measured channel scenarios have shown significant temporal correlation [21].(ii) Despite being of great interest in practice, the estimation of the number of received paths is often not addressed. The crucial problem here is to correctly estimate the current number of paths to avoid over determination, since an overdetermined estimator will tend to use the additional degrees of freedom to match the noise by introducing erroneous paths. There exist various techniques based on model selection that can be employed to estimate the number of paths [22] but they suffer from the problem of having to dynamically adjust the decision thresholds. Typically, only a single hypothesis is tracked, which in practice causes error event propagation.(iii) The ML estimator does only provide the most likely parameter set for the given observation. No reliability information about the estimates is provided. Consequently, ambiguities and multiple modes of the likelihood are not preserved by the estimator.

ML estimators require that the estimated parameters remain constant during the observation period. Due to data modulation and phase variations in practice, this period, which is often referred to as the coherent integration time, is limited to a range of 1 millisecond–20 milliseconds. To reach a sufficient noise performance with a ML estimator in practice, it is required to extend its observation interval approximately to the equivalent averaging time of a conventional tracking loop, which is commonly in the order of several hundred coherent integration periods. To overcome this problem, the observations are forced to be quasicoherent by aiding the ML estimator with a phase locked loop (PLL) and a data removal mechanism [8].

4. Sequential Estimation

4.1. Optimal Solution

In Section 3, we have established the models of the underlying time variant processes. The problem of multipath mitigation now becomes one of sequential estimation of a hidden Markov process. We want to estimate the unknown channel parameters based on an evolving sequence of received noisy channel outputs 𝐳_𝑘. The channel process for each range of a satellite navigation system can be modeled as a first-order Markov process if future channel parameters given the present state of the channel and all its past states depend only on the present channel state (and not on any past states). We also assume that the noise affecting successive channel outputs is independent of the past noise values; so each channel observation depends only on the present channel state.

Intuitively, we are exploiting not only the channel observations to estimate the hidden channel parameters (via the likelihood function), but we are also exploiting our prior knowledge about the statistical dependencies between successive sets of channel parameters. We know from channel measurements that channel parameters are time varying but not independent from one time instance to the next; for example, an echo usually experiences a “life-cycle” from its first occurrence, then a more or less gradual change in its delay and phase over time, until it disappears [21]. These measurements also show that the common channel models considered in communication systems [17, 18] are not adequately reflecting the properties that are crucial for high-resolution signal delay estimation as required in navigation systems.

Now that our major assumptions have been established, we may apply the concept of sequential Bayesian estimation. The reader is referred to [23] which gives a derivation of the general framework for optimal estimation of temporally evolving (Markovian) parameters by means of inference, and we have chosen similar notation. The entire history of observations (over the temporal index 𝑘) can be written as𝐙_𝑘={𝐳_𝑘′,𝑘^′=1,…,𝑘}. (16)Similarly, we denote the sequence of parameters of our hidden Markovian process by𝐗_𝑘={𝐱_𝑘′,𝑘^′=1,…,𝑘}. (17)Here, 𝐱_𝑘′ represents the characterization of the hidden channel state, including the parameters that specify the signal hypothesis 𝐬_𝑘′ given in (2). Our goal is to determine the posterior PDF of every possible channel characterization given all channel observations: 𝑝(𝐱_𝑘∣𝐙_𝑘) (see Figure 3). Once we have evaluated this posterior PDF, we can either determine that channel configuration that maximizes it—the so-called maximum a posteriori (MAP) estimate; or we can choose the expectation—equivalent to the minimum mean square error (MMSE) estimate. In addition, the posterior distribution itself contains all uncertainty about the current range and is thus the optimal measure to perform sensor data fusion in an overall positioning system.

Figure 3 

Illustration of the hidden Markov estimation process for three time instances. Our channel measurements are the sequence $\{z_i\mathrm{,}\hspace{0.17em}\hspace{0.17em}i=1,\dots \mathrm{,}k\}$, and the channel parameters to be estimated are $\{x_i\mathrm{,}\hspace{0.17em}\hspace{0.17em}i=1,\dots \mathrm{,}k\}$.

It can be shown that the sequential estimation algorithm is recursive, as it uses the posterior PDF computed for time instance 𝑘−1 to compute the posterior PDF for instance 𝑘 (see Figure 4). For a given posterior PDF at time instance 𝑘−1, 𝑝(𝐱_𝑘−1∣𝐙_𝑘−1), the prior PDF 𝑝(𝐱_𝑘∣𝐙_𝑘−1) is calculated in the so-called prediction step by applying the Chapman-Kolmogorov equation:𝑝(𝐱_𝑘∣𝐙_𝑘−1)=𝑝(𝐱_𝑘∣𝐱_𝑘−1)𝑝(𝐱_𝑘−1∣𝐙_𝑘−1)d𝐱_𝑘−1, (18)with 𝑝(𝐱_𝑘∣𝐱_𝑘−1) being the state transition PDF of the Markov process. In the update step, the new posterior PDF for step 𝑘 is obtained by applying Bayes' rule to 𝑝(𝐱_𝑘∣𝐳_𝑘,𝐙_𝑘−1) yielding the normalized product of the likelihood 𝑝(𝐳_𝑘∣𝐱_𝑘) and the prior PDF:𝑝(𝐱_𝑘∣𝐙_𝑘)=𝑝(𝐱_𝑘∣𝐳_𝑘,𝐙_𝑘−1)=𝑝(𝐳_𝑘∣𝐱_𝑘,𝐙_𝑘−1)𝑝(𝐱_𝑘∣𝐙_𝑘−1)𝑝(𝐳_𝑘∣𝐙_𝑘−1)=𝑝(𝐳_𝑘∣𝐱_𝑘)𝑝(𝐱_𝑘∣𝐙_𝑘−1)𝑝(𝐳_𝑘∣𝐙_𝑘−1). (19)The term 𝑝(𝐳_𝑘∣𝐱_𝑘)=𝑝(𝐳_𝑘∣𝐬_𝑐,𝑘=𝐬^{0𝑥002𝑑𝑐}_𝑐,𝑘) follows from (4) and represents the probability of the measured channel output (often referred to as the likelihood value), conditioned on a certain configuration of channel parameters at the same time step 𝑘. The denominator of (19) does not depend on 𝐱_𝑘, and so it can be computed by integrating the numerator of (19) over the entire range of 𝐱_𝑘 (normalization).

Figure 4 

Illustration of the recursive Bayesian estimator.

To summarize so far, the entire process of prediction and update can be carried out recursively to calculate the posterior PDF (19) sequentially, based on an initial value of 𝑝(𝐱₀∣𝐳₀)=𝑝(𝐱₀). The evaluation of the likelihood function 𝑝(𝐳_𝑘∣𝐱_𝑘) is the essence of the update step. Similarly, maximizing this likelihood function (i.e., ML estimation) would be equivalent to maximizing 𝑝(𝐱_𝑘∣𝐙_𝑘) only in the case that the prior PDF 𝑝(𝐱_𝑘∣𝐙_𝑘−1) does not depend on 𝐙_𝑘−1, and when all values of 𝑥_𝑘 are a priori equally likely. Since these conditions are not met, evaluation of 𝑝(𝐱_𝑘∣𝐙_𝑘) entails all the above steps.

4.2. Sequential Estimation Using Particle Filters

The optimal estimation algorithm relies on evaluating the integral (18), which is usually a very difficult task, except for certain additional restrictions imposed on the model and the noise process. So, very often a suboptimal realization of a Bayesian estimator has to be chosen for implementation. In this paper, we use a sequential Monte Carlo (SMC) filter, in particular a sampling importance resampling particle filter (SIR-PF) according to [23]. In this algorithm, the posterior density at step 𝑘 is represented as a sum and is specified by a set of 𝑁_𝑝 particles:𝑝(𝐱_𝑘∣𝐙_𝑘)≈^𝑁𝑝_𝑗=1𝑤^𝑗_𝑘⋅𝛿(𝐱_𝑘−𝐱^𝑗_𝑘), (20)where each particle with index 𝑗 has a state 𝐱^𝑗_𝑘 and has a weight𝑤^𝑗_𝑘. The sum over all particles' weights is one. In SIR-PF, the weights are computed according to the principle of importance sampling where the so-called proposal density is chosen to be 𝑝(𝐱_𝑘∣𝐱_𝑘−1=𝐱^𝑗_𝑘−1), and with resampling at every time step. For 𝑁_𝑝→∞, the approximate posterior approaches the true PDF.

The key step, in which the measurement for instance 𝑘 is incorporated, is in the calculation of the weight 𝑤^𝑗_𝑘 which for the SIR-PF can be shown to be the likelihood function: 𝑝(𝐳_𝑘∣𝐱^𝑗_𝑘). The characterization of the channel process enters in the algorithm, when at each time instance 𝑘, the state of each particle 𝐱^𝑖_𝑘 is drawn randomly from the proposal distribution, that is, from 𝑝(𝐱_𝑘∣𝐱^𝑗_𝑘−1).

4.3. Exploiting Linear Substructures

If there exist linear substructures in the model, it is possible to reduce the computational complexity of the filter by means of marginalization over the linear state variables [24], also known as Rao-Blackwellization [25]. In a marginalized filter, particles are still used to estimate the nonlinear states, while for each of the particles the linear states can be estimated analytically. In our case, since the measurement 𝐳_𝑘 is a linear function of the complex amplitudes 𝐚_𝑘, the likelihood function can be factorized as𝑝(𝐳_𝑘∣𝐬_𝑐,𝑘)=𝑝(𝐳_𝑘∣𝝉_𝑘,𝐞_𝑘,𝐚_𝑘)=𝑓(𝐳_𝑘,𝝉_𝑘,𝐞_𝑘)⋅𝑔(𝐳_𝑘,𝝉_𝑘,𝐞_𝑘,𝐚_𝑘), (21)where the function 𝑔(𝐳_𝑘,𝝉_𝑘,𝐞_𝑘,𝐚_𝑘) is Gaussian with respect to 𝐚_𝑘, and 𝑓(𝐳_𝑘,𝝉_𝑘,𝐞_𝑘) is proportional to 𝑝(𝐳_𝑘∣𝐬^{0𝑥002𝑑𝑐}_𝑐,𝑘). Marginalization of (21) over the linear state variables gives_𝐚𝑘𝑝(𝐳_𝑘∣𝝉_𝑘,𝐞_𝑘,𝐚_𝑘)𝑑𝐚_𝑘=𝑓(𝐳_𝑘,𝝉_𝑘,𝐞_𝑘)∝𝑝(𝐳_𝑘∣𝐬^{0𝑥002𝑑𝑐}_𝑐,𝑘). (22)Assuming that the amplitudes are block-wise independent (𝑄=1), it follows that the weights of the SIR particle filter are equal to 𝑝(𝐳_𝑘∣𝐱^𝑗_𝑘)=𝑝(𝐳_𝑘∣𝐬_𝑐,𝑘=𝐬^{0𝑥002𝑑𝑐𝑗}_𝑐,𝑘), which is given by (4).

If there exists a temporal correlation of the amplitudes, an optimal marginalized filter requires the implementation of a separate Kalman filter for each of the particles. As we do not want to increase the complexity per particle, we take the correlation into account by adjustment of 𝑄 in (13).

4.4. Choice of Appropriate Channel Process

To exploit the advantages of sequential estimation for our task of multipath mitigation/estimation, we must be able to describe the actual channel characteristics (channel parameters) so that these are captured by 𝑝(𝐱_𝑘∣𝐱_𝑘−1). In other words, the model must be a first-order Markov model, and all transition probabilities must be known. In our approach, we approximate the channel as follows.(i) The channel is totally characterized by a direct path (index 𝑖=1) and at most 𝑁_𝑚−1 echoes,(ii)each path has complex amplitude 𝑎_𝑖,𝑘 and relative delay and >Δ𝜏_𝑖,𝑘=𝜏_𝑖,𝑘−𝜏_1,𝑘; where echoes are constrained to have delay 𝜏_𝑖,𝑘≥𝜏_1,𝑘, that is, >Δ𝜏_𝑖,𝑘≥0,(iii)the different path delays follow the process:𝜏_1,𝑘=𝜏_1,𝑘−1+𝛼_1,𝑘−1⋅Δ𝑡+𝑛_𝜏,Δ𝜏_𝑖,𝑘=Δ𝜏_{𝑖,𝑘−1}+𝛼_{𝑖,𝑘−1}⋅Δ𝑡+𝑛_𝜏, (23)(iv) each parameter 𝛼_𝑖,𝑘 that specifies the speed of the change of the path delay follows its own process:𝛼_𝑖,𝑘1=(1−𝐾)⋅𝛼_{𝑖,𝑘−1}+𝑛_𝛼, (24)(v)the magnitudes and phases of the individual paths, represented by the complex amplitudes 𝑎_𝑖,𝑘, are eliminated according to (12) and (14) for the computation of the likelihood (4),(vi) each path is either “on” or “off,” as defined by channel parameter 𝑒_𝑖,𝑘∈{1≡ “on,” 0≡ “off” },(vii) where 𝑒_𝑖,𝑘 follows a simple two-state Markov process with a-symmetric crossover and same-state probabilities:𝑝(𝑒_𝑖,𝑘=0∣𝑒_{𝑖,𝑘−1}=1)=𝑝_onoff,𝑝(𝑒_𝑖,𝑘=1∣𝑒_{𝑖,𝑘−1}=0)=𝑝_offon. (25)(26) The model implicitly incorporates three i.i.d noise sources: Gaussian 𝑛_𝜏 and 𝑛_𝛼 as well as the noise process driving the state changes for 𝑒_𝑖,𝑘. These sources provide the randomness of the model. The parameter 𝐾 defines how quickly the speed of path delays can change (for a given variance of 𝑛_𝛼). Finally, >Δ𝑡 is the time between instances 𝑘−1 and 𝑘. We assume all model parameters (i.e., 𝐾, >Δ𝑡, noise variances, and the “on”/“off” Markov model) to be independent of 𝑘 (see Figure 5). The hidden channel state vector 𝐱_𝑘 can therefore be represented as[𝜏_1,𝑘,>Δ𝜏_2,𝑘,…,>Δ𝜏_{𝑁𝑚,𝑘},𝛼_1,𝑘,…,𝛼_{𝑁𝑚,𝑘},𝑒_𝑖,𝑘,…,𝑒_{𝑁𝑚,𝑘}]^𝑇. (27)Note that the model implicitly represents the number of paths𝑁_𝑚,𝑘=^𝑁𝑚_𝑖=1𝑒_𝑖,𝑘 (1)as a time-variant parameter.

Figure 5 

The Markov process chosen in this paper to model the channel with $N_m$ paths. Dotted arrows—shown only for a small subset of the transitions—indicate the constraint that $\mathrm{>\Delta }{\tau }_{i\mathrm{,}k}\ge 0$. For an explanation of the terms, see the text leading up to (27).

When applied to our particle filtering algorithm, drawing from the proposal density is simple. Each particle stores the above-channel parameters of the model, and then the new state of each particle is drawn randomly from 𝑝(𝐱_𝑘∣𝐱^𝑗_𝑘−1) which corresponds to drawing values for 𝑛_𝛼 and 𝑛_𝜏 as well as propagating the “on”/“off” Markov model, and then updating the channel parameters for 𝑘 according to (23)–(26).

4.5. Practical Issues

4.5.1. Model Matching

It is important to point out that a sequential estimator is only as good as its state transition model matches the real world situation. The state model needs to capture all relevant hidden states with memory and needs to correctly model their dependencies, while adhering to the first order Markov condition. Furthermore, any memory of the measurement noise affecting the likelihood function 𝑝(𝐳_𝑘∣𝐱_𝑘) must be explicitly contained as additional states of the model 𝐱, so that the measurement noise is i.i.d.

The channel state model is motivated by channel modeling work for multipath prone environments such as the urban satellite navigation channel [21, 26]. In fact, the process of constructing a channel model in order to characterize the channel for signal level simulations and receiver evaluation comes close to our task of building a first-order Markov process for sequential estimation. For particle filtering, the model needs to satisfy the condition that one can draw states with relatively low computational complexity. Adapting the model structure and the model parameters to the real channel environment is a task for current and future work. It may even be possible to envisage hierarchical models in which the selection of the current model itself follows a process. In this case, a sequential estimator will automatically choose the correct weighting of these models according to their ability to fit the received signal.

4.5.2. Integration into a Receiver

For receiver integration, the computational complexity of the filtering algorithm is crucial. From a theoretical point of view, it is desirable to run the sequential filter clocked corresponding to the coherent integration period of the receiver and with a very large number of particles. From the practical point of view, however, it is desirable to reduce the sequential filter rate to the navigation rate and to minimize the number of particles. Existing ML approaches can help here to achieve a flexible complexity/performance tradeoff, as strategies already developed to extend the observation periods of ML estimators can be used directly to reduce the rate of the sequential filtering algorithm.

5. Performance Evaluation

To demonstrate the capabilities of the SMC-based Bayesian time delay estimator proposed in Section 4, we have carried out computer simulations for both static and dynamic channel environments. The signal-to-multipath ratio (SMR) was chosen constant and equal to 6 dB, while the relative phases are changing according to >Δ𝜑_𝑖,𝑘=2𝜋>Δ𝜏_𝑖,𝑘𝑓_𝑐 with 𝑓_𝑐=1575.42 MHz being the frequency of the L1 carrier. The Bayesian estimator uses a time interval of >Δ𝑡=1 millisecond corresponding to the duration of a codeword. The amplitude averaging coefficient is set to 𝑄=10, and signal compression is applied with 𝑁_𝑐𝑐=41 (code-matched correlators) 𝑁_𝑝𝑐=25. The channel parameters 𝜎²_𝑖,𝜏, 𝜎²_𝑖,𝛼, 𝐾, and 𝑝(𝑒_𝑖,𝑘∣𝑒_{𝑖,𝑘−1}) are selected to fit the statistics of a real channel according to [21]. The SIR-PF uses the minimum mean square error (MMSE) criterion to estimate the parameters 𝐱_𝑘 from the posterior.

5.1. Static Multipath Scenario

The capability of multipath mitigation techniques is commonly assessed by showing the systematic error due to a single multipath replica plotted as a function of the relative multipath delay in a static channel scenario. In Figure 6, the root mean square error (RMSE) is shown for the proposed sequential estimator, implemented as a SIR-PF with 𝑁_𝑠=2000 particles. For comparison, the performance of conventional DLLs with Narrow Correlator and with Strobe Correlator is also shown. The simulated signal corresponds to a GPS L1 signal with 𝑐(𝑡) being a Gold code of length 1023 that is modulated on a bandlimited rectangular pulse. The chip rate is 1.023 Mchips/s so that the duration of the codeword is 1millisecond. The one-sided bandwidth of the resulting navigation signal is 5 MHz. Estimators with fixed two-path model or fixed single path model are also shown for comparison with the implicit path activity tracking. The performance of a single path estimator is comparable to that of a DLL with infinitesimal correlator spacing and shows a considerable bias over a large delay range. The estimator with fixed two-path model successfully mitigates the multipath bias for delays greater than 30 m. However, for smaller delays it shows an increasing variance and is outperformed by the single path estimator. The estimator with path activity tracking is capable of combining the advantages of both models. From the posterior, it is possible to calculate the estimated average probability 𝑃(𝑁_𝑚,𝑘=2∣𝐙_𝑘) of a two-path model, which is shown in Figure 7, and indicates the transition between the models: for small delays the two paths essentially merge to a single one. Note that in these simulations, the model parameters of the sequential estimator are still the ones designed for the dynamic channel and not optimal for this static scenario. The RMSE as function of the number of particles 𝑁_𝑠 is shown in Figure 8 for a relative multipath delay of 14 m, which corresponds to the worst case in Figure 6.

Figure 6 

Static scenario: performance of SIR-PF for BPSK modulation as function of relative multipath delay for different path models.

Figure 7 

Static scenario: average probability of a two-path model for the estimator with path activity tracking.

Figure 8 

Static scenario: RMSE performance as function of number of particles $N_s$ for BPSK modulation.

5.2. Dynamic Scenario

The dynamic multipath channel scenario with up to 𝑁_𝑚=3 paths, used in the following simulations and depicted in Figure 9, has been generated according to the movement model, whereby the parameters 𝐾=25000, 𝜎_𝑛𝛼=10⁻¹⁰, 𝜎_𝑛𝜏=10⁻⁸, and 𝑝_onoff=𝑝_offon=0.0001 were chosen to resemble a typical urban satellite navigation channel environment [26].

Figure 9 

Considered dynamic multipath scenario: pseudoranges over time of the direct path (continuous line) and the temporarily present echoes.

Figure 10 shows for this multipath scenario at a carrier-to-noise ratio of 𝐶/𝑁₀=45 dB-Hz the BPSK and the BOC-(1,1) performance of a DLL with a narrow correlator spacing of 10⁻⁷ seconds, corresponding to a tenth of a code chip. A DLL loop bandwidth of 1 Hz was selected, which appeared to deliver the best DLL performance for this channel. For a fair comparison of the different modulation techniques, both signals were generated with the same C/A code sequence of length 1023, the same receiver bandwidth of 5 MHz (one-sided), a sampling rate of 1/𝑇_𝑠=10.23 MHz, and a coherent integration time of 𝐿𝑇_𝑠=1 millisecond. Although the results show an improvement for the BOC-(1,1) modulation, there still remains a considerable error due to the multipath. This behavior is confirmed by a comparison with the multipath error envelope, depicted in Figure 11. While in some multipath delay regions, the error is reduced by the BOC-(1,1) modulation (see, e.g, the segment marked by the dashed box in Figure 10), more sophisticated multipath mitigating techniques are required to reduce the errors due to short range multipath.

Figure 10 

Dynamic scenario: performance of the DLL with narrow correlator for BPSK modulation (upper/grey) and BOC-(1,1) modulation (lower/black). The segment marked with a dashed box shows a situation where BOC(1,1) is clearly superior to BPSK.

Figure 11 

Multipath error envelope of the DLL with narrow correlator for BPSK modulation (dashed) and BOC(1,1) modulation (inner/solid).

The simulation results for the particle filter-based estimator with 𝑁_𝑠=20000 particles are given in Figures 12(a) and 12(b) for the same channel as in Figure 10. They show that the performance could be improved substantially, resulting in a reduction of the root mean square (RMS) error from 3.77 m to 0.769 m for BPSK modulation and from 2.61 m to 0.694 m for BOC(1,1) modulation, respectively.

(a) BPSK modulation

(b) BOC(1,1) modulation

Figure 12 

Dynamic scenario: Performance of the sequential estimator with particle filtering (black) compared to the DLL (grey). The sequential estimator is clearly superior for multipath mitigation.

The RMSE as function of the number of particles 𝑁_𝑠 is shown in Figure 13.

Figure 13 

Dynamic scenario: RMSE performance as function of number of particles $N_s$ for BPSK modulation.

6. Conclusions

We have demonstrated how sequential Bayesian estimation techniques can be applied to the multipath mitigation problem in a navigation receiver. The proposed approach is characterized by code-matched, correlator-based signal compression together with interpolation techniques for efficient likelihood computation in combination with a particle filter realization of the prediction and update recursion. The considered movement model has been adapted to dynamic multipath scenarios and incorporates the number of echoes as a time variant hidden channel state variable that is tracked together with the other parameters in a probabilistic fashion. A further advantage compared to ML estimation is that the posterior PDF at the output of the estimator represents reliability information about the desired parameters and preserves the ambiguities and multiple modes that may occur within the likelihood function. Simulation results for BPSK- and BOC-(1,1) modulated signals show that in both cases significant improvements can be achieved compared to a DLL with narrow correlator. In this work, we have employed two methods to reduce complexity: the signal compression to facilitate the computation of the likelihood function as well as a simple form of Rao-Blackwellization to eliminate the complex amplitudes from the state space. Further work will concentrate on additional complexity reduction techniques such as more suitable proposal functions or particle filtering algorithms such as the auxiliary particle filter that are possibly more efficient with respect to the number of particles when applied to our problem domain.