Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments

Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2008, Article ID 372651, 11 pages doi:10.1155/2008/372651 Research Article Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments 1 2 2 Michael Lentmaier, Bernhard Krach, and Patrick Robertson Vodafone Chair Communications Systems, Dresden University of Technology, TU Dresden, 01062 Dresden, Germany German Aerospace Center DLR, Institute of Communications and Navigation, Oberpfaffenhofen, 82234 Wessling, Germany Correspondence should be addressed to Michael Lentmaier, michael.lentmaier@ifn.et.tu-dresden.de Received 1 August 2007; Revised 7 December 2007; Accepted 17 March 2008 Recommended by Letizia Presti 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. Copyright © 2008 Michael Lentmaier et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1. INTRODUCTION 2. A COMPARISON OF VARIOUS MULTIPATH MITIGATION APPROACHES A major error source within global navigation satellite systems (GNSSs) comes from multipath, the reception of Figure 1 gives an overview of the relationships between additional signal replica due to reflections, which introduce different multipath mitigation and estimation approaches. In a bias into the estimate of the delay lock loop (DLL) of fact, we have chosen to discriminate approaches according a conventional navigation receiver. For efficient removal of to their primary objective. The left column represents the this bias, it is possible to formulate advanced maximum class of techniques that attempt to mitigate the effect of likelihood (ML) estimators that incorporate the echoes into multipath in different ways. This can for example be achieved the signal model and are capable of achieving the theoretical by modifications of the antenna response, either by means of limits given by the Cramer ´ Rao bound. The drawback of hardware design or with signal processing techniques (e.g., ML estimator techniques is that the parameters are assumed beamforming). The majority of the remaining mitigation to be constant during the time of observation. Independent techniques are in some way aligning the more or less estimates are obtained for successive observation intervals, traditional receiver components (e.g., the early/late corre- whose length has to be adapted to the dynamics of the lator) to the signal received in the multipath environment. channel. The probably most simple technique is the adjustment of In this paper, we consider the important practical case the correlator spacing applied in the Narrow Correlator of dynamic channel scenarios and assess how the time- [4]. Other well-known examples of this category are the delay estimation can be improved if information is available Strobe Correlator [5], the Gated Correlator [6], or the Pulse about the temporal evolution of the channel parameters. Our Aperture Correlator [7]. To incorporate new signal forms approach is based on Bayesian filtering, the optimal and well- (such as BOC), these methods need “tuning” in order to known framework to address such dynamic state estimation suffer as little as possible from multipath. On the other problems [1]. Sequential Monte Carlo (SMC) methods [2, 3] hand, multipath estimation techniques (right column) treat are used for computing the posterior probability density multipath (in particular the delay of the paths) as something functions (PDFs) of the signal parameters. to be estimated from the channel observations, so that its 2 International Journal of Navigation and Observation Mitigation Estimation Static Dynamic Modification of standard DLL detector Maximum likelihood Sequential estimation • Narrow correlator • MEDLL • Bayesian filtering • Double-delta correlator • Vision correlator & ◦ Kalman filter • Strobe correlator MMT variants • Pulse aperture correlator • SAGE ◦ Sequential Monte • Newton-type Carlo methods • Antenna array signal . . processing techniques Antenna characteristics . Practical ML integration: • Choke ring • Signal compression • Beamforming • Loop-aided ML . • ML-in-the-loop MAP with static prior Figure 1: Classification of multipath mitigation approaches. effects can be trivially removed at a later processing stage. For e (t) is a binary function that controls the activity of the the estimation techniques, we have differentiated between i th path, and a (t)and τ (t) are their individual complex i i static and dynamic approaches, according to the underlying amplitudes and time delays, respectively. The signal is assumption of the channel dynamics. Examples for static disturbed by additive white Gaussian noise n(t). Grouping multipath estimation are those belonging to the family of blocks of L samples at times (m + kL)T , m = 0,... , L − 1, ML estimators, often using different efficient maximization together into vectors z , k = 0, 1,... , n whilst assuming strategies over the likelihood function [8–13]. For static the parameter functions e (t), a (t), and τ (t) being constant i i i channels without availability of prior information, the ML within the corresponding time interval and equal to e , a , i,k i,k approach is optimal and performs significantly better than and τ , this can be rewritten as i,k other techniques, especially if the echoes have short delay. An estimator based on sequential importance sampling (SIS) z = CG τ E a + n = s + n . (2) k k k k k k k methods (particle filtering) for static multipath scenarios has been considered in [14], which has the advantage that prior In the compact form on the right hand side, the samples of channel knowledge can be incorporated. the delayed pulses g(τ ) are stacked together as columns i,k As a first step towards addressing dynamic channels, one of the matrix G(τ ) = [g(τ ),... , g(τ )], C is a k 1,k N ,k can incorporate ML estimators in receiver loops or formulate matrix representing the convolution with the code, and the quasisequential estimators [15, 16]. Finally, dynamic esti- delays and amplitudes are collected in the vectors τ = mators that target the computation of the posterior PDF T T [τ ,... , τ ] and a = [a ,... , a ] ,respectively. 1,k N ,k k 1,k N ,k m m conditioned on the received channel output sequence at the Furthermore, for concise notation we use E = diag [e ] k k receiver can be applied on a per single range basis or operate whilst the elements of the vector e = [e ,... , e ] , k 1,k N ,k in the position domain. In this paper, we concentrate on m e ∈ [0, 1], determine whether the i th path is active or not i,k dynamic estimators applied per each range. The sequential by being either e = 1 corresponding to an active path or i,k Monte Carlo approach has also been suggested in the com- e = 0 for a path that is currently not active. The term s i,k k munications field for estimation of time-varying channel denotes the signal hypothesis and is completely determined responses in spread spectrum systems [17, 18]. by the channel parameters τ , a and e . Using (2), we can k k, k write the associated likelihood function as 3. SIGNAL MODEL Assume that the complex valued baseband-equivalent 1 1 p z | s = · exp − z − s z − s . k k k k k k received signal is equal to L 2 2L 2σ (2π ) σ N (3) z(t) = e (t)·a (t)· c(t)∗g t − τ (t) + n(t), (1) i i i i=1 The likelihood function will play a central role in the where c(t) is a delta-train code sequence that is modulated on algorithms discussed in this paper; its purpose is to quantify a pulse g (t), N is the maximum number of allowed paths the conditional probability of the received signal conditioned reaching the receiver (to restrict the modeling complexity), on the unknown signal (specifically the channel parameters). Michael Lentmaier et al. 3 3.1. Efficient likelihood computation into a canonical component decomposition,given by an L×N cc matrix Q ,and a principal component decomposition,given cc In [11], a general concept for the efficient representation by an N × N matrix Q .In[11], two choices for Q are cc pc pc cc of the likelihood (3) was presented, which is applicable to proposed many of the existing ML multipath mitigation methods. b −1 The key idea of this concept is to formulate (3) through ⎨ CG(τ )R Signal matched, cc Q = (8) avector z resulting from an orthonormal projection of cc c,k b −1 C(τ )R Code matched, cc the observed signal z onto a smaller vector space, so that z is a sufficient statistic according to the Neyman-Fisher c,k where the elements of the vector τ define the positions factorization [19] and hence suitable for estimating s .In of the individual correlators. The noise-free outputs of the other words, the reduced signal comprises the same infor- corresponding correlator banks are illustrated in Figure 2.To mation as the original signal itself. In practice, this concept H H b b decorrelate the bank outputs (CG(τ )) y and C(τ ) y,as becomesrelevantasthe projection canbeachievedby mentioned above, the whitening matrix R can be obtained cc processing the received signal (2) with a bank of correlators b b from a QR decomposition of CG(τ )and C(τ ), respectively. and a subsequent decorrelation of the correlator outputs. Apart from practical implementation issues, both correlation A variant of this very general concept, applied in [13], has methods given by (8) are equivalent from a conceptual point also been referred to as the Signal Compression Theorem of view. For details on the compression through Q the pc, in [20] for a set of special projections that do not require reader is referred to [11]. the step of decorrelation due to the structure of the used correlators. For instance, unlike the correlation technique 3.1.2. Interpolation used in [8], the one suggested in [13] already projects onto an orthogonal and thus uncorrelated subspace, similar to the In order to compute (4) independently of the sampling grid, code matched correlator technique proposed in [11]. Due to advantage can be made of interpolation techniques. Using complexity reasons, all practically relevant realizations of ML the discrete Fourier transformation (DFT), with Ψ being the estimators [8, 13] operate in a projected space, namely after −1 DFT matrix and Ψ being its inverse counterpart (IDFT), correlation. The corresponding mathematical background we get will be discussed below, including also interpolation of the −1 likelihood and elimination of complex amplitudes as further H s = Q CΨ diag Ψg(0) Ω τ E a c,k k k k methods for complexity reduction. (9) = M Ω τ E a , s k k k 3.1.1. Data compression with Ω(τ ) being a matrix of column-wise stacked vectors with Vandermonde structure [10, 11], such that the element As explained above, the large vector containing the received at row p and column q computes with signal samples z is linearly transformed into a vector z of k c,k much smaller size. Following this approach, the likelihood 2π (p − 1)τ q,k according to (2)can be rewrittenas R Ω τ = cos , p,q N T g s (10) z z p z | s = exp − 2π (p − 1)τ k k q,k 2L 2σ (2π ) σ I Ω τ =− sin . p,q N T g s H H H H R z Q Q s s Q Q s c k c k c c k k · exp − N is the length of the pulse g in samples. The advantage 2 2 σ 2σ of the interpolation is that it can take place in the reduced (4) 1 z z k space. The most costly computations in (9) can be carried out = exp − L 2 2L in precalculations as the matrix M , whose row dimension 2σ s (2π ) σ c corresponds to the dimension of the reduced space and H H R z s s s c,k c,k c,k c,k whose column dimension is N , is constant. · exp − , 2 2 σ 2σ 3.1.3. Amplitude elimination with the compressed received vector z and the compressed c,k signal hypothesis s : c,k In a further step, we reduce the number of parameters by H H z = Q z , s = Q s,(5) optimizing (4) for a given set of τ and e with respect to c,k k c,k k c c k k the complex amplitudes a , which can be achieved through a and the orthonormal compression matrix Q , which needs to closed-form solution. Using fulfill S = M Ω(τ )E (11) H H c,k s k k Q Q ≈ I, Q Q ≈ I,(6) c c c c and obtaining S by removing zero columns from S one c,k, to minimize the compression loss. According to [11], the c,k yields the corresponding amplitude values of the active paths: compression can be two-fold so that we can factorize −1 + +H + +H a = (S S ) S z . (12) Q = Q Q (7) c cc pc c,k k c,k c,k c,k 4 International Journal of Navigation and Observation 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Delay (μs) Delay (μs) (a) BPSK signal matched (b) BPSK code matched 1.2 1.5 0.8 0.6 0.5 0.4 0.2 0 −0.5 −0.2 −0.4 −1 −0.6 −0.8 −1.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Delay (μs) Delay (μs) (c) BOC(1,1) signal matched (d) BOC(1,1) code matched b b Figure 2: Output of the canonical component type correlator banks CG(τ )and C(τ ) for BPSK and BOC(1,1). As we have introduced a potential source of performance 3.2. Review of the ML Concept loss by eliminating the amplitudes and thus practically The concept of ML multipath estimation has drawn substan- are disregarding their temporal correlation, we propose to tial research interest since the first approach was proposed in optimize (4) using [8]. Despite being treated differently in various publications, the objective is the same for all ML approaches, namely to Q−1 find the signal parameters that maximize (3)for agiven z = · z (13) c,k c,k−l, observation z : l=0 s = arg max p z | s . (15) k k k with the adjustable averaging coefficient Q. When evaluating (4), we substitute s by c,k The signal parameters are thereby assumed as being constant throughout the observation period k.Different maximiza- s = S  a , (14) c,k c,k k tion strategies exist, which basically characterize the different approaches. Despite offering great advantages for theoretical where the elements of the vector  a that are indicated to analysis, the practical advantage of the generic ML concept is have an active path (a : i→e = 1) areset equalto k,i k,i questionable due to a number of serious drawbacks. the corresponding elements of  a . All other elements (a : k,i i→e = 0) can be set arbitrarily as their influence is masked k,i (i) The ML estimator assumes that the channel is static by the zero elements of e . for the observation period and is not able to exploit k Michael Lentmaier et al. 5 its temporal correlation throughout the sequence k = 1,... , n. Measured channel scenarios have shown Measurements (observed) significant temporal correlation [21]. (ii) Despite being of great interest in practice, the esti- z z z k−2 k−1 k mation of the number of received paths is often not addressed. The crucial problem here is to correctly Likelihood Time Time Time estimate the current number of paths to avoid over k − 2 k − 1 k p(z |x ) k k determination, since an overdetermined estimator will tend to use the additional degrees of freedom to match the noise by introducing erroneous paths. x x x k−2 k−1 k There exist various techniques based on model selec- p(x |x ) k k−1 tion that can be employed to estimate the number State transition PDF of paths [22] but they suffer from the problem of Hidden states having to dynamically adjust the decision thresholds. Typically, only a single hypothesis is tracked, which in practice causes error event propagation. Figure 3: Illustration of the hidden Markov estimation process for (iii) The ML estimator does only provide the most likely three time instances. Our channel measurements are the sequence parameter set for the given observation. No relia- {z , i = 1,... , k}, and the channel parameters to be estimated are {x , i = 1,... , k}. bility information about the estimates is provided. i Consequently, ambiguities and multiple modes of the likelihood are not preserved by the estimator. measurements that channel parameters are time varying but ML estimators require that the estimated parameters not independent from one time instance to the next; for remain constant during the observation period. Due to data example, an echo usually experiences a “life-cycle” from its modulation and phase variations in practice, this period, first occurrence, then a more or less gradual change in its which is often referred to as the coherent integration time, delay and phase over time, until it disappears [21]. These is limited to a range of 1 millisecond–20 milliseconds. To measurements also show that the common channel models reach a sufficient noise performance with a ML estimator considered in communication systems [17, 18]are not in practice, it is required to extend its observation interval adequately reflecting the properties that are crucial for high- approximately to the equivalent averaging time of a con- resolution signal delay estimation as required in navigation ventional tracking loop, which is commonly in the order of systems. several hundred coherent integration periods. To overcome Now that our major assumptions have been established, this problem, the observations are forced to be quasicoherent we may apply the concept of sequential Bayesian estimation. by aiding the ML estimator with a phase locked loop (PLL) The reader is referred to [23] which gives a derivation of and a data removal mechanism [8]. the general framework for optimal estimation of temporally evolving (Markovian) parameters by means of inference, 4. SEQUENTIAL ESTIMATION and we have chosen similar notation. The entire history of observations (over the temporal index k)can be writtenas 4.1. Optimal solution Z ={  z , k = 1,... , k}. (16) k k In Section 3, we have established the models of the under- lying time-variant processes. The problem of multipath Similarly, we denote the sequence of parameters of our mitigation now becomes one of sequential estimation of a hidden Markovian process by hidden Markov process. We want to estimate the unknown channel parameters based on an evolving sequence of X ={  x , k = 1,... , k}. (17) k k received noisy channel outputs z . The channel process for each range of a satellite navigation system can be modeled Here, x represents the characterization of the hidden as a first-order Markov process if future channel parameters channel state, including the parameters that specify the signal given the present state of the channel and all its past states hypothesis s  given in (2). Our goal is to determine the depend only on the present channel state (and not on any posterior PDF of every possible channel characterization past states). We also assume that the noise affecting successive given all channel observations: p(x | Z ) (see Figure 3). k k channel outputs is independent of the past noise values; so Once we have evaluated this posterior PDF, we can either each channel observation depends only on the present channel determine that channel configuration that maximizes it—the state. so-called maximum a posteriori (MAP) estimate; or we can Intuitively, we are exploiting not only the channel choose the expectation—equivalent to the minimum mean observations to estimate the hidden channel parameters (via square error (MMSE) estimate. In addition, the posterior the likelihood function), but we are also exploiting our distribution itself contains all uncertainty about the current prior knowledge about the statistical dependencies between range and is thus the optimal measure to perform sensor data successive sets of channel parameters. We know from channel fusion in an overall positioning system. 6 International Journal of Navigation and Observation It can be shown that the sequential estimation algorithm Movement model is recursive, as it uses the posterior PDF computed for time Prediction p(x |x ) stage k k−1 instance k − 1 to compute the posterior PDF for instance k (see Figure 4). For a given posterior PDF at time instance k− 1, p(x | Z ), the prior PDF p(x | Z ) is calculated k−1 k−1 k k−1 Prior p(x |Z ) in the so-called prediction step by applying the Chapman- k k−1 k = k +1 Kolmogorov equation: Measurements Likelihood Update p x | Z = p x | x p x | Z dx , (18) p(z |x ) k k−1 k k−1 k−1 k−1 k−1 k k stage with p(x | x ) being the state transition PDF of the k k−1 Posterior Markov process. In the update step, the new posterior PDF for p(x |Z ) k k step k is obtained by applying Bayes’ rule to p(x | z , Z ) k k k−1 yielding the normalized product of the likelihood p(z | x ) Figure 4: Illustration of the recursive Bayesian estimator. k k and the prior PDF: p x | Z = p x | z , Z k k k k k−1 where each particle with index j has a state x and has a p z | x , Z p x | Z weight w . The sum over all particles’ weights is one. In SIR- k k k−1 k k−1 PF, the weights are computed according to the principle of p z | Z (19) k k−1 importance sampling where the so-called proposal density p z | x p x | Z k k k k−1 = . is chosen to be p(x | x = x ), and with resampling k k−1 k−1 p z | Z k k−1 at every time step. For N →∞, the approximate posterior approaches the true PDF. The term p(z | x ) = p(z | s =  s ) follows from (4)and k k k c,k c,k The key step, in which the measurement for instance k represents the probability of the measured channel output is incorporated, is in the calculation of the weight w which (often referred to as the likelihood value), conditioned on a k for the SIR-PF can be shown to be the likelihood function: certain configuration of channel parameters at the same time p(z | x ). The characterization of the channel process enters step k. The denominator of (19)doesnot depend on x ,and k in the algorithm, when at each time instance k, the state so it can be computed by integrating the numerator of (19) of each particle x is drawn randomly from the proposal over the entire range of x (normalization). k To summarize so far, the entire process of prediction distribution, that is, from p(x | x ). k−1 and update can be carried out recursively to calculate the posterior PDF (19) sequentially, based on an initial value of 4.3. Exploiting linear substructures p(x | z ) = p(x ). The evaluation of the likelihood function 0 0 0 If there exist linear substructures in the model, it is possible p(z | x ) is the essence of the update step. Similarly, k k maximizing this likelihood function (i.e., ML estimation) to reduce the computational complexity of the filter by means of marginalization over the linear state variables [24], would be equivalent to maximizing p(x | Z ) only in the k k also known as Rao-Blackwellization [25]. In a marginalized case that the prior PDF p(x | Z )doesnot depend on k k−1 Z and when all values of x are a priori equally likely. filter, particles are still used to estimate the nonlinear states, k−1, k while for each of the particles the linear states can be Since these conditions are not met, evaluation of p(x | Z ) k k estimated analytically. In our case, since the measurement entails all the above steps. z is a linear function of the complex amplitudes a , the k k likelihood function can be factorized as 4.2. Sequential estimation using particle filters p z | s = p z |τ , e , a = f z , τ , e ·g z , τ , e , a , k c,k k k k k k k k k k k k The optimal estimation algorithm relies on evaluating the (21) integral (18), whichisusually averydifficult task, except for certain additional restrictions imposed on the model and the where the function g (z , τ , e , a ) is Gaussian with respect k k k k noise process. So, very often a suboptimal realization of a to a ,and f (z , τ , e ) is proportional to p(z |  s ). k k k k k c,k Bayesian estimator has to be chosen for implementation. In Marginalization of (21) over the linear state variables gives this paper, we use a sequential Monte Carlo (SMC) filter, in particular a sampling importance resampling particle filter p z | τ , e , a da = f z , τ , e ∝ p z |  s . k k k k k k k k k c,k (SIR-PF) according to [23]. In this algorithm, the posterior (22) density at step k is represented as a sum and is specified by a set of N particles: Assuming that the amplitudes are block-wise independent (Q = 1), it follows that the weights of the SIR particle filter j j j j are equal to p(z | x ) = p(z | s =  s ), whichisgiven by k k c,k p x | Z ≈ w ·δ x − x , (20) k c,k k k k k k (4). j=1 Michael Lentmaier et al. 7 If there exists a temporal correlation of the amplitudes, Note that the model implicitly represents the number of an optimal marginalized filter requires the implementation paths of a separate Kalman filter for each of the particles. As we do not want to increase the complexity per particle, we take the N = e (28) m,k i,k correlation into account by adjustment of Q in (13). i=1 as a time-variant parameter. 4.4. Choice of appropriate channel process When applied to our particle filtering algorithm, drawing To exploit the advantages of sequential estimation for from the proposal density is simple. Each particle stores our task of multipath mitigation/estimation, we must be the above-channel parameters of the model, and then the able to describe the actual channel characteristics (channel new state of each particle is drawn randomly from p(x | parameters) so that these are captured by p(x | x ). k k−1 x ) which corresponds to drawing values for n and n as α τ k−1 In other words, the model must be a first-order Markov well as propagating the “on”/“off ” Markov model, and then model, and all transition probabilities must be known. In our updating the channel parameters for k according to (23)– approach, we approximate the channel as follows. (26). (i) The channel is totally characterized by a direct path (index i = 1) and at most N − 1echoes, 4.5. Practical issues (ii) each path has complex amplitude a and relative i,k 4.5.1. Model matching delay and Δτ = τ − τ ; where echoes are i,k i,k 1,k constrained to have delay τ ≥ τ , that is, Δτ ≥ 0, i,k 1,k i,k It is important to point out that a sequential estimator is only as good as its state transition model matches the real (iii) the different path delays follow the process: world situation. The state model needs to capture all relevant τ = τ + α ·Δt + n , 1,k 1,k−1 1,k−1 τ hidden states with memory and needs to correctly model (23) their dependencies, while adhering to the first order Markov Δτ = Δτ + α ·Δt + n , i,k i,k−1 i,k−1 τ condition. Furthermore, any memory of the measurement noise affecting the likelihood function p(z | x )mustbe k k (iv) each parameter α that specifies the speed of the i,k explicitly contained as additional states of the model x,so change of the path delay follows its own process: that the measurement noise is i.i.d. The channel state model is motivated by channel α = 1 − ·α + n , (24) i,k i,k−1 α 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 (v) the magnitudes and phases of the individual paths, characterize the channel for signal level simulations and represented by the complex amplitudes a , are elimi- i,k receiver evaluation comes close to our task of building a nated according to (12)and (14) for the computation first-order Markov process for sequential estimation. For of the likelihood (4), particle filtering, the model needs to satisfy the condition (vi) each path is either “on” or “off,” as defined by channel that one can draw states with relatively low computational parameter e ∈{1 ≡ “on, 0 ≡ “off }, i,k complexity. Adapting the model structure and the model (vii) where e follows a simple two-state Markov process i,k parameters to the real channel environment is a task for with a-symmetric crossover and same-state probabil- current and future work. It may even be possible to envisage ities: hierarchical models in which the selection of the current model itself follows a process. In this case, a sequential p e = 0 | e = 1 = p , (25) i,k i,k−1 onoff estimator will automatically choose the correct weighting of these models according to their ability to fit the received p e = 1 | e = 0 = p . (26) i,k i,k−1 offon signal. The model implicitly incorporates three i.i.d noise sources: Gaussian n and n as well as the noise process driving the τ α 4.5.2. Integration into a receiver state changes for e . These sources provide the randomness i,k For receiver integration, the computational complexity of of the model. The parameter K defines how quickly the the filtering algorithm is crucial. From a theoretical point speed of path delays can change (for a given variance of n ). of view, it is desirable to run the sequential filter clocked Finally, Δt is the time between instances k − 1and k.We corresponding to the coherent integration period of the assume all model parameters (i.e., K , Δt, noise variances, and receiver and with a very large number of particles. From the the “on”/“off ” Markov model) to be independent of k (see practical point of view, however, it is desirable to reduce the Figure 5). The hidden channel state vector x can therefore sequential filter rate to the navigation rate and to minimize be represented as the number of particles. Existing ML approaches can help τ , Δτ ,... , Δτ , α ,... , α , e ,... , e . here to achieve a flexible complexity/performance tradeoff, 1,k 2,k Nm ,k 1,k Nm ,k i,k Nm ,k as strategies already developed to extend the observation (27) 8 International Journal of Navigation and Observation τ τ 1,k−2 1,k−1 1,k α α α Direct path 1,k−2 1,k−1 1,k e e e 1,k−2 1,k−1 1,k Δτ Δτ 2,k−1 Δτ 2,k−2 2,k First echo α α α 2,k−2 2,k−1 2,k e e e 2,k−2 2,k−1 2,k . . . . . . ··· ··· . . . State x State x State x k−2 k−1 k Δτ Δτ Δτ N ,k−2 N ,k−1 N ,k m m m α α α N ,k−2 N ,k−1 N ,k N − 1’th echo m m m m N ,k−2 e e m N ,k−1 N ,k m m Figure 5: The Markov process chosen in this paper to model the channel with N paths. Dotted arrows—shown only for a small subset of the transitions—indicate the constraint that Δτ ≥ 0. For an explanation of the terms, see the text leading up to (27). i,k 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,wehavecarried 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πΔτ f with i,k i,k c f = 1575.42 MHz being the frequency of the L1 carrier. The 0 0 10 20 30 40 50 60 70 Bayesian estimator uses a time interval of Δt = 1 millisecond Relative delay (m) corresponding to the duration of a codeword. The amplitude averaging coefficient is set to Q = 10, and signal compression SIR-PF, path activitiy tracking DLL narrow correlator is applied with N = 41 (code-matched correlators) N = SIR-PF, single path model DLL strobe correlator cc pc 2 2 SIR-PF, two path model 25. Thechannelparameters σ , σ , K,and p(e | e )are i,k i,k−1 i,τ i,α selected to fit the statistics of a real channel according to [21]. Figure 6: Static scenario: performance of SIR-PF for BPSK The SIR-PF uses the minimum mean square error (MMSE) modulation as function of relative multipath delay for different path criterion to estimate the parameters x from the posterior. models. 5.1. Static multipath scenario The capability of multipath mitigation techniques is com- and with Strobe Correlator is also shown. The simulated monly assessed by showing the systematic error due to signal corresponds to a GPS L1 signal with c(t) being a Gold a single multipath replica plotted as a function of the code of length 1023 that is modulated on a bandlimited relative multipath delay in a static channel scenario. In rectangular pulse. The chip rate is 1.023 Mchips/s so that Figure 6, the root mean square error (RMSE) is shown the duration of the codeword is 1millisecond. The one- for the proposed sequential estimator, implemented as a sided bandwidth of the resulting navigation signal is 5 MHz. SIR-PF with N = 2000 particles. For comparison, the Estimators with fixed two-path model or fixed single path performance of conventional DLLs with Narrow Correlator model are also shown for comparison with the implicit path RMSE (m) Michael Lentmaier et al. 9 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 20 40 60 80 100 120 0 10 20 30 40 50 60 70 Time (s) Relative delay (m) Figure 9: Considered dynamic multipath scenario: pseudoranges Figure 7: Static scenario: average probability of a two-path model over time of the direct path (continuous line) and the temporarily for the estimator with path activity tracking. present echoes. −5 20 40 60 80 100 120 500 1000 1500 2000 2500 3000 Time (s) Number of particles Figure 10: Dynamic scenario: performance of the DLL with narrow Figure 8: Static scenario: RMSE performance as function of correlator for BPSK modulation (upper/grey) and BOC(1,1) modu- number of particles N for BPSK modulation. lation (lower/black). The segment marked with a dashed box shows a situation where BOC(1,1) is clearly superior to BPSK. activity tracking. The performance of a single path estimator −8 is comparable to that of a DLL with infinitesimal correlator σ = 10 ,and p = p = 0.0001 were chosen nτ onoff offon spacing and shows a considerable bias over a large delay to resemble a typical urban satellite navigation channel range. The estimator with fixed two-path model successfully environment [26]. mitigates the multipath bias for delays greater than 30 m. Figure 10 shows for this multipath scenario at a carrier- However, for smaller delays it shows an increasing variance to-noise ratio of C/N = 45 dB-Hz the BPSK and the and is outperformed by the single path estimator. The BOC(1,1) performance of a DLL with a narrow correlator −7 estimator with path activity tracking is capable of combining spacing of 10 seconds, corresponding to a tenth of a the advantages of both models. From the posterior, it code chip. A DLL loop bandwidth of 1 Hz was selected, is possible to calculate the estimated average probability which appeared to deliver the best DLL performance for this P(N = 2 | Z ) of a two-path model, which is shown in m,k k channel. For a fair comparison of the different modulation Figure 7, and indicates the transition between the models: techniques, both signals were generated with the same C/A for small delays the two paths essentially merge to a single code sequence of length 1023, the same receiver bandwidth one. Note that in these simulations, the model parameters of 5 MHz (one-sided), a sampling rate of 1/T = 10.23 MHz, of the sequential estimator are still the ones designed for the and a coherent integration time of LT = 1 millisecond. dynamic channel and not optimal for this static scenario. The Although the results show an improvement for the BOC(1,1) RMSE as function of the number of particles N is shown modulation, there still remains a considerable error due to in Figure 8 for a relative multipath delay of 14 m, which the multipath. This behavior is confirmed by a comparison corresponds to the worst case in Figure 6. 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 5.2. Dynamic scenario by the dashed box in Figure 10), more sophisticated multi- The dynamic multipath channel scenario with up to N = 3 path mitigating techniques are required to reduce the errors paths, used in the following simulations and depicted in due to short range multipath. Figure 9, has been generated according to the movement The simulation results for the particle filter-based estima- −10 model, whereby the parameters K = 25000, σ = 10 , tor with N = 20000 particles are given in Figures 12(a) and n s Probability of two path model RMSE (m) Pseudorange (m) Error (m) 10 International Journal of Navigation and Observation −5 −10 0 50 100 150 200 250 300 350 0.5 1 1.5 2 2.5 3 ×10 Relative delay (m) Number of particles Figure 11: Multipath error envelope of the DLL with narrow cor- Figure 13: Dynamic scenario: RMSE performance as function of relator for BPSK modulation (dashed) and BOC(1,1) modulation number of particles N for BPSK modulation. (inner/solid). 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 effi- cient likelihood computation in combination with a particle −5 filter realization of the prediction and update recursion. The considered movement model has been adapted to dynamic −10 multipath scenarios and incorporates the number of echoes 0 20 40 60 80 100 120 as a time-variant hidden channel state variable that is tracked Time (s) together with the other parameters in a probabilistic fashion. (a) BPSK modulation 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 −5 as a simple form of Rao-Blackwellization to eliminate the complex amplitudes from the state space. Further work will 0 20 40 60 80 100 120 concentrate on additional complexity reduction techniques Time (s) such as more suitable proposal functions or particle filtering (b) BOC(1,1) modulation algorithms such as the auxiliary particle filter that are possibly more efficient with respect to the number of Figure 12: Dynamic scenario: Performance of the sequential particles when applied to our problem domain. estimator with particle filtering (black) compared to the DLL (grey). The sequential estimator is clearly superior for multipath mitigation. REFERENCES [1] A. H. Jazwinski, Stochastic Processes and Filtering Theory, Academic Press, New York, NY, USA, 1970. [2] A. Doucet, N. de Freitas, and N. Gordon, Eds., Sequential 12(b) for the same channel as in Figure 10. They show that Monte Carlo Methods in Practice, Springer, New York, NY, the performance could be improved substantially, resulting USA, 2001. in a reduction of the root mean square (RMS) error from [3] B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman 3.77 m to 0.769 m for BPSK modulation and from 2.61 m to Filter: Particle Filters for Tracking Applications,ArtechHouse, 0.694 m for BOC(1,1) modulation, respectively. Boston, Mass, USA, 2004. The RMSE as function of the number of particles N is s [4] A.J.Van Dierendonck, P. C. Fenton,and T. Ford,“Theory and shown in Figure 13. performance of narrow correlator spacing in a GPS receiver,” Error (m) Error envelope (m) Error (m) RMSE (m) Michael Lentmaier et al. 11 in Proceedings of the ION National Technical Meeting,San [18] T. Bertozzi, D. L. Ruyet, C. Panazio, and H. V. Thien, “Channel Diego, Calif, USA, January 1992. tracking using particle filtering in unresolvable multipath environments,” EURASIP Journal on Applied Signal Processing, [5] L. Garin, F. van Diggelen, and J.-M. Rousseau, “Strobe and vol. 2004, no. 15, pp. 2328–2338, 2004. edge correlator multipath mitigation for code,” in Proceedings of the 9th International Technical Meeting of the Satellite [19] S. M. Kay, Fundamentals of Statistical Signal Processing: Division of the Institute of Navigation (ION GPS ’96), pp. 657– Estimation Theory, Prentice Hall, Upper Saddle River, NJ, 664, Kansas City, Mo, USA, September 1996. USA, 1993. [6] G. MacGraw and M. Brasch, “GNSS multipath mitigation [20] L. R. Weill, “Achieving theoretical bounds for receiver-based using gated and high resolution correlator concepts,” in multipath mitigation using galileo OS signals,” in Proceedings Proceedings of the ION National Technical Meeting,San Diego, of the 19th International Technical Meeting of the Satellite Calif, USA, January 1999. Division of the Institute of Navigation (ION GNSS ’06),pp. 1035–1047, Fort Worth, Tex, USA, September 2006. [7] J.Jones,P.C.Fenton, andB.Smith,“Theory andperformance of the pulse aperture correlator,” NovAtel Technical Report, [21] A. Steingass and A. Lehner, “Measuring the navigation NovAtel, Calgary, Alberta, Canada, 2004. multipath channel—a statistical analysis,” in Proceedings of the 17th International Technical Meeting of the Satellite Division [8] R. D. J. van Nee, J. Siereveld, P. C. Fenton, and B. R. Townsend, of the Institute of Navigation (ION GNSS ’04), pp. 1157–1164, “The Multipath estimating delay lock loop: approaching Long Beach, Calif, USA, September 2004. theoretical accuracy limits,” in Proceedings of the IEEE Position Location and Navigation Symposium (PLANS ’94), pp. 246– [22] P. Stoica, Y. Selen, ´ and J. Li, “On information criteria and 251, Las Vegas, Nev, USA, April 1994. the generalized likelihood ratio test of model order selection,” IEEE Signal Processing Letters, vol. 11, no. 10, pp. 794–797, [9] L. R. Weill, “Achieving theoretical accuracy limits for pseudo- ranging in the presence of multipath,” in Proceedings of the 8th International Technical Meeting of the Satellite Division of the [23] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A Institute of Navigation (ION GPS ’95), pp. 1521–1530, Palm tutorial on particle filters for online nonlinear/non-Gaussian Springs, Calif, USA, September 1995. Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002. [10] J. Selva Vera, “Complexity reduction in the parametric estimation of superimposed signal replicas,” Signal Processing, [24] T. Schon, ¨ F. Gustafsson, and P.-J. Nordlund, “Marginalized vol. 84, no. 12, pp. 2325–2343, 2004. particle filters for mixed linear/nonlinear state-space models,” IEEE Transactions on Signal Processing, vol. 53, no. 7, pp. 2279– [11] J. Selva Vera, “Efficient multipath mitigation in naviga- 2289, 2005. tion systems,” Ph.D. dissertation, Universitat Politecnica de Catalunya, Barcelona, Spain, February 2004. [25] A. Doucet, N. de Freitas, K. Murphy, and S. Russell, “Rao- blackwellised particle filtering for dynamic Bayesian net- [12] F. Antreich, O. Esbr´ ı-Rodr´ ıguez, J. A. Nossek, and W. Utchick, works,” in Proceedings of the 16th Annual Conference on “Estimation of synchronization parameters using SAGE in Uncertainty in Artificial Intelligence (UAI ’00), pp. 176–183, a GNSS-receiver,” in Proceedings of the 18th International San Francisco, Calif, USA, June-July 2000. Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’05), pp. 2124–2131, Long Beach, [26] A. Steingass and A. Lehner, “Land mobile satellite na- Calif, USA, September 2005. vigation—characteristics of the multipath channel,” in Pro- ceedings of the 16th International Technical Meeting of the [13] P. C. Fenton and J. Jones, “The theory and performance of Satellite Division of the Institute of Navigation (ION GNSS ’03), NovAtel Inc.’s vision correlator,” in Proceedings of the 18th Portland, Ore, USA, September 2003. International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’05), pp. 2178–2186, Long Beach, Calif, USA, September 2005. [14] P. Closas, C. Fernandez-P ´ rades, J. A. Fernandez-R ´ ubio, and A. Ram´ ırez-Gonzalez, ´ “Multipath mitigation using particle filtering,” in Proceedings of the 19th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’06), pp. 1733–1740, Fort Worth, Tex, USA, September 2006. [15] M. Lentmaier and B. Krach, “Maximum likelihood multipath estimation in comparison with conventional delay lock loops,” in Proceedings of the 19th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’06), pp. 1741–1751, Fort Worth, Tex, USA, September 2006. [16] B. Krach and M. Lentmaier, “Efficient soft-output GNSS signal parameter estimation using signal compression techniques,” in Proceedings of the 3rd ESA Workshop on Satellite Navigation User Equipment Technologies (NAVITECH ’06), Noordwijk, The Netherlands, December 2006. [17] E. Punskaya, A. Doucet, and W. J. Fitzgerald, “Particle filtering for joint symbol and code delay estimation in DS spread spec- trum systems in multipath environment,” EURASIP Journal on Applied Signal Processing, vol. 2004, no. 15, pp. 2306–2314, 2004. International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Navigation and Observation Hindawi Publishing Corporation

Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments

Loading next page...
 
/lp/hindawi-publishing-corporation/bayesian-time-delay-estimation-of-gnss-signals-in-dynamic-multipath-WqcPU5Wfor
Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2008 Michael Lentmaier et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
ISSN
1687-5990
DOI
10.1155/2008/372651
Publisher site
See Article on Publisher Site

Abstract

Hindawi Publishing Corporation International Journal of Navigation and Observation Volume 2008, Article ID 372651, 11 pages doi:10.1155/2008/372651 Research Article Bayesian Time Delay Estimation of GNSS Signals in Dynamic Multipath Environments 1 2 2 Michael Lentmaier, Bernhard Krach, and Patrick Robertson Vodafone Chair Communications Systems, Dresden University of Technology, TU Dresden, 01062 Dresden, Germany German Aerospace Center DLR, Institute of Communications and Navigation, Oberpfaffenhofen, 82234 Wessling, Germany Correspondence should be addressed to Michael Lentmaier, michael.lentmaier@ifn.et.tu-dresden.de Received 1 August 2007; Revised 7 December 2007; Accepted 17 March 2008 Recommended by Letizia Presti 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. Copyright © 2008 Michael Lentmaier et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 1. INTRODUCTION 2. A COMPARISON OF VARIOUS MULTIPATH MITIGATION APPROACHES A major error source within global navigation satellite systems (GNSSs) comes from multipath, the reception of Figure 1 gives an overview of the relationships between additional signal replica due to reflections, which introduce different multipath mitigation and estimation approaches. In a bias into the estimate of the delay lock loop (DLL) of fact, we have chosen to discriminate approaches according a conventional navigation receiver. For efficient removal of to their primary objective. The left column represents the this bias, it is possible to formulate advanced maximum class of techniques that attempt to mitigate the effect of likelihood (ML) estimators that incorporate the echoes into multipath in different ways. This can for example be achieved the signal model and are capable of achieving the theoretical by modifications of the antenna response, either by means of limits given by the Cramer ´ Rao bound. The drawback of hardware design or with signal processing techniques (e.g., ML estimator techniques is that the parameters are assumed beamforming). The majority of the remaining mitigation to be constant during the time of observation. Independent techniques are in some way aligning the more or less estimates are obtained for successive observation intervals, traditional receiver components (e.g., the early/late corre- whose length has to be adapted to the dynamics of the lator) to the signal received in the multipath environment. channel. The probably most simple technique is the adjustment of In this paper, we consider the important practical case the correlator spacing applied in the Narrow Correlator of dynamic channel scenarios and assess how the time- [4]. Other well-known examples of this category are the delay estimation can be improved if information is available Strobe Correlator [5], the Gated Correlator [6], or the Pulse about the temporal evolution of the channel parameters. Our Aperture Correlator [7]. To incorporate new signal forms approach is based on Bayesian filtering, the optimal and well- (such as BOC), these methods need “tuning” in order to known framework to address such dynamic state estimation suffer as little as possible from multipath. On the other problems [1]. Sequential Monte Carlo (SMC) methods [2, 3] hand, multipath estimation techniques (right column) treat are used for computing the posterior probability density multipath (in particular the delay of the paths) as something functions (PDFs) of the signal parameters. to be estimated from the channel observations, so that its 2 International Journal of Navigation and Observation Mitigation Estimation Static Dynamic Modification of standard DLL detector Maximum likelihood Sequential estimation • Narrow correlator • MEDLL • Bayesian filtering • Double-delta correlator • Vision correlator & ◦ Kalman filter • Strobe correlator MMT variants • Pulse aperture correlator • SAGE ◦ Sequential Monte • Newton-type Carlo methods • Antenna array signal . . processing techniques Antenna characteristics . Practical ML integration: • Choke ring • Signal compression • Beamforming • Loop-aided ML . • ML-in-the-loop MAP with static prior Figure 1: Classification of multipath mitigation approaches. effects can be trivially removed at a later processing stage. For e (t) is a binary function that controls the activity of the the estimation techniques, we have differentiated between i th path, and a (t)and τ (t) are their individual complex i i static and dynamic approaches, according to the underlying amplitudes and time delays, respectively. The signal is assumption of the channel dynamics. Examples for static disturbed by additive white Gaussian noise n(t). Grouping multipath estimation are those belonging to the family of blocks of L samples at times (m + kL)T , m = 0,... , L − 1, ML estimators, often using different efficient maximization together into vectors z , k = 0, 1,... , n whilst assuming strategies over the likelihood function [8–13]. For static the parameter functions e (t), a (t), and τ (t) being constant i i i channels without availability of prior information, the ML within the corresponding time interval and equal to e , a , i,k i,k approach is optimal and performs significantly better than and τ , this can be rewritten as i,k other techniques, especially if the echoes have short delay. An estimator based on sequential importance sampling (SIS) z = CG τ E a + n = s + n . (2) k k k k k k k methods (particle filtering) for static multipath scenarios has been considered in [14], which has the advantage that prior In the compact form on the right hand side, the samples of channel knowledge can be incorporated. the delayed pulses g(τ ) are stacked together as columns i,k As a first step towards addressing dynamic channels, one of the matrix G(τ ) = [g(τ ),... , g(τ )], C is a k 1,k N ,k can incorporate ML estimators in receiver loops or formulate matrix representing the convolution with the code, and the quasisequential estimators [15, 16]. Finally, dynamic esti- delays and amplitudes are collected in the vectors τ = mators that target the computation of the posterior PDF T T [τ ,... , τ ] and a = [a ,... , a ] ,respectively. 1,k N ,k k 1,k N ,k m m conditioned on the received channel output sequence at the Furthermore, for concise notation we use E = diag [e ] k k receiver can be applied on a per single range basis or operate whilst the elements of the vector e = [e ,... , e ] , k 1,k N ,k in the position domain. In this paper, we concentrate on m e ∈ [0, 1], determine whether the i th path is active or not i,k dynamic estimators applied per each range. The sequential by being either e = 1 corresponding to an active path or i,k Monte Carlo approach has also been suggested in the com- e = 0 for a path that is currently not active. The term s i,k k munications field for estimation of time-varying channel denotes the signal hypothesis and is completely determined responses in spread spectrum systems [17, 18]. by the channel parameters τ , a and e . Using (2), we can k k, k write the associated likelihood function as 3. SIGNAL MODEL Assume that the complex valued baseband-equivalent 1 1 p z | s = · exp − z − s z − s . k k k k k k received signal is equal to L 2 2L 2σ (2π ) σ N (3) z(t) = e (t)·a (t)· c(t)∗g t − τ (t) + n(t), (1) i i i i=1 The likelihood function will play a central role in the where c(t) is a delta-train code sequence that is modulated on algorithms discussed in this paper; its purpose is to quantify a pulse g (t), N is the maximum number of allowed paths the conditional probability of the received signal conditioned reaching the receiver (to restrict the modeling complexity), on the unknown signal (specifically the channel parameters). Michael Lentmaier et al. 3 3.1. Efficient likelihood computation into a canonical component decomposition,given by an L×N cc matrix Q ,and a principal component decomposition,given cc In [11], a general concept for the efficient representation by an N × N matrix Q .In[11], two choices for Q are cc pc pc cc of the likelihood (3) was presented, which is applicable to proposed many of the existing ML multipath mitigation methods. b −1 The key idea of this concept is to formulate (3) through ⎨ CG(τ )R Signal matched, cc Q = (8) avector z resulting from an orthonormal projection of cc c,k b −1 C(τ )R Code matched, cc the observed signal z onto a smaller vector space, so that z is a sufficient statistic according to the Neyman-Fisher c,k where the elements of the vector τ define the positions factorization [19] and hence suitable for estimating s .In of the individual correlators. The noise-free outputs of the other words, the reduced signal comprises the same infor- corresponding correlator banks are illustrated in Figure 2.To mation as the original signal itself. In practice, this concept H H b b decorrelate the bank outputs (CG(τ )) y and C(τ ) y,as becomesrelevantasthe projection canbeachievedby mentioned above, the whitening matrix R can be obtained cc processing the received signal (2) with a bank of correlators b b from a QR decomposition of CG(τ )and C(τ ), respectively. and a subsequent decorrelation of the correlator outputs. Apart from practical implementation issues, both correlation A variant of this very general concept, applied in [13], has methods given by (8) are equivalent from a conceptual point also been referred to as the Signal Compression Theorem of view. For details on the compression through Q the pc, in [20] for a set of special projections that do not require reader is referred to [11]. the step of decorrelation due to the structure of the used correlators. For instance, unlike the correlation technique 3.1.2. Interpolation used in [8], the one suggested in [13] already projects onto an orthogonal and thus uncorrelated subspace, similar to the In order to compute (4) independently of the sampling grid, code matched correlator technique proposed in [11]. Due to advantage can be made of interpolation techniques. Using complexity reasons, all practically relevant realizations of ML the discrete Fourier transformation (DFT), with Ψ being the estimators [8, 13] operate in a projected space, namely after −1 DFT matrix and Ψ being its inverse counterpart (IDFT), correlation. The corresponding mathematical background we get will be discussed below, including also interpolation of the −1 likelihood and elimination of complex amplitudes as further H s = Q CΨ diag Ψg(0) Ω τ E a c,k k k k methods for complexity reduction. (9) = M Ω τ E a , s k k k 3.1.1. Data compression with Ω(τ ) being a matrix of column-wise stacked vectors with Vandermonde structure [10, 11], such that the element As explained above, the large vector containing the received at row p and column q computes with signal samples z is linearly transformed into a vector z of k c,k much smaller size. Following this approach, the likelihood 2π (p − 1)τ q,k according to (2)can be rewrittenas R Ω τ = cos , p,q N T g s (10) z z p z | s = exp − 2π (p − 1)τ k k q,k 2L 2σ (2π ) σ I Ω τ =− sin . p,q N T g s H H H H R z Q Q s s Q Q s c k c k c c k k · exp − N is the length of the pulse g in samples. The advantage 2 2 σ 2σ of the interpolation is that it can take place in the reduced (4) 1 z z k space. The most costly computations in (9) can be carried out = exp − L 2 2L in precalculations as the matrix M , whose row dimension 2σ s (2π ) σ c corresponds to the dimension of the reduced space and H H R z s s s c,k c,k c,k c,k whose column dimension is N , is constant. · exp − , 2 2 σ 2σ 3.1.3. Amplitude elimination with the compressed received vector z and the compressed c,k signal hypothesis s : c,k In a further step, we reduce the number of parameters by H H z = Q z , s = Q s,(5) optimizing (4) for a given set of τ and e with respect to c,k k c,k k c c k k the complex amplitudes a , which can be achieved through a and the orthonormal compression matrix Q , which needs to closed-form solution. Using fulfill S = M Ω(τ )E (11) H H c,k s k k Q Q ≈ I, Q Q ≈ I,(6) c c c c and obtaining S by removing zero columns from S one c,k, to minimize the compression loss. According to [11], the c,k yields the corresponding amplitude values of the active paths: compression can be two-fold so that we can factorize −1 + +H + +H a = (S S ) S z . (12) Q = Q Q (7) c cc pc c,k k c,k c,k c,k 4 International Journal of Navigation and Observation 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Delay (μs) Delay (μs) (a) BPSK signal matched (b) BPSK code matched 1.2 1.5 0.8 0.6 0.5 0.4 0.2 0 −0.5 −0.2 −0.4 −1 −0.6 −0.8 −1.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 Delay (μs) Delay (μs) (c) BOC(1,1) signal matched (d) BOC(1,1) code matched b b Figure 2: Output of the canonical component type correlator banks CG(τ )and C(τ ) for BPSK and BOC(1,1). As we have introduced a potential source of performance 3.2. Review of the ML Concept loss by eliminating the amplitudes and thus practically The concept of ML multipath estimation has drawn substan- are disregarding their temporal correlation, we propose to tial research interest since the first approach was proposed in optimize (4) using [8]. Despite being treated differently in various publications, the objective is the same for all ML approaches, namely to Q−1 find the signal parameters that maximize (3)for agiven z = · z (13) c,k c,k−l, observation z : l=0 s = arg max p z | s . (15) k k k with the adjustable averaging coefficient Q. When evaluating (4), we substitute s by c,k The signal parameters are thereby assumed as being constant throughout the observation period k.Different maximiza- s = S  a , (14) c,k c,k k tion strategies exist, which basically characterize the different approaches. Despite offering great advantages for theoretical where the elements of the vector  a that are indicated to analysis, the practical advantage of the generic ML concept is have an active path (a : i→e = 1) areset equalto k,i k,i questionable due to a number of serious drawbacks. the corresponding elements of  a . All other elements (a : k,i i→e = 0) can be set arbitrarily as their influence is masked k,i (i) The ML estimator assumes that the channel is static by the zero elements of e . for the observation period and is not able to exploit k Michael Lentmaier et al. 5 its temporal correlation throughout the sequence k = 1,... , n. Measured channel scenarios have shown Measurements (observed) significant temporal correlation [21]. (ii) Despite being of great interest in practice, the esti- z z z k−2 k−1 k mation of the number of received paths is often not addressed. The crucial problem here is to correctly Likelihood Time Time Time estimate the current number of paths to avoid over k − 2 k − 1 k p(z |x ) k k determination, since an overdetermined estimator will tend to use the additional degrees of freedom to match the noise by introducing erroneous paths. x x x k−2 k−1 k There exist various techniques based on model selec- p(x |x ) k k−1 tion that can be employed to estimate the number State transition PDF of paths [22] but they suffer from the problem of Hidden states having to dynamically adjust the decision thresholds. Typically, only a single hypothesis is tracked, which in practice causes error event propagation. Figure 3: Illustration of the hidden Markov estimation process for (iii) The ML estimator does only provide the most likely three time instances. Our channel measurements are the sequence parameter set for the given observation. No relia- {z , i = 1,... , k}, and the channel parameters to be estimated are {x , i = 1,... , k}. bility information about the estimates is provided. i Consequently, ambiguities and multiple modes of the likelihood are not preserved by the estimator. measurements that channel parameters are time varying but ML estimators require that the estimated parameters not independent from one time instance to the next; for remain constant during the observation period. Due to data example, an echo usually experiences a “life-cycle” from its modulation and phase variations in practice, this period, first occurrence, then a more or less gradual change in its which is often referred to as the coherent integration time, delay and phase over time, until it disappears [21]. These is limited to a range of 1 millisecond–20 milliseconds. To measurements also show that the common channel models reach a sufficient noise performance with a ML estimator considered in communication systems [17, 18]are not in practice, it is required to extend its observation interval adequately reflecting the properties that are crucial for high- approximately to the equivalent averaging time of a con- resolution signal delay estimation as required in navigation ventional tracking loop, which is commonly in the order of systems. several hundred coherent integration periods. To overcome Now that our major assumptions have been established, this problem, the observations are forced to be quasicoherent we may apply the concept of sequential Bayesian estimation. by aiding the ML estimator with a phase locked loop (PLL) The reader is referred to [23] which gives a derivation of and a data removal mechanism [8]. the general framework for optimal estimation of temporally evolving (Markovian) parameters by means of inference, 4. SEQUENTIAL ESTIMATION and we have chosen similar notation. The entire history of observations (over the temporal index k)can be writtenas 4.1. Optimal solution Z ={  z , k = 1,... , k}. (16) k k In Section 3, we have established the models of the under- lying time-variant processes. The problem of multipath Similarly, we denote the sequence of parameters of our mitigation now becomes one of sequential estimation of a hidden Markovian process by hidden Markov process. We want to estimate the unknown channel parameters based on an evolving sequence of X ={  x , k = 1,... , k}. (17) k k received noisy channel outputs z . The channel process for each range of a satellite navigation system can be modeled Here, x represents the characterization of the hidden as a first-order Markov process if future channel parameters channel state, including the parameters that specify the signal given the present state of the channel and all its past states hypothesis s  given in (2). Our goal is to determine the depend only on the present channel state (and not on any posterior PDF of every possible channel characterization past states). We also assume that the noise affecting successive given all channel observations: p(x | Z ) (see Figure 3). k k channel outputs is independent of the past noise values; so Once we have evaluated this posterior PDF, we can either each channel observation depends only on the present channel determine that channel configuration that maximizes it—the state. so-called maximum a posteriori (MAP) estimate; or we can Intuitively, we are exploiting not only the channel choose the expectation—equivalent to the minimum mean observations to estimate the hidden channel parameters (via square error (MMSE) estimate. In addition, the posterior the likelihood function), but we are also exploiting our distribution itself contains all uncertainty about the current prior knowledge about the statistical dependencies between range and is thus the optimal measure to perform sensor data successive sets of channel parameters. We know from channel fusion in an overall positioning system. 6 International Journal of Navigation and Observation It can be shown that the sequential estimation algorithm Movement model is recursive, as it uses the posterior PDF computed for time Prediction p(x |x ) stage k k−1 instance k − 1 to compute the posterior PDF for instance k (see Figure 4). For a given posterior PDF at time instance k− 1, p(x | Z ), the prior PDF p(x | Z ) is calculated k−1 k−1 k k−1 Prior p(x |Z ) in the so-called prediction step by applying the Chapman- k k−1 k = k +1 Kolmogorov equation: Measurements Likelihood Update p x | Z = p x | x p x | Z dx , (18) p(z |x ) k k−1 k k−1 k−1 k−1 k−1 k k stage with p(x | x ) being the state transition PDF of the k k−1 Posterior Markov process. In the update step, the new posterior PDF for p(x |Z ) k k step k is obtained by applying Bayes’ rule to p(x | z , Z ) k k k−1 yielding the normalized product of the likelihood p(z | x ) Figure 4: Illustration of the recursive Bayesian estimator. k k and the prior PDF: p x | Z = p x | z , Z k k k k k−1 where each particle with index j has a state x and has a p z | x , Z p x | Z weight w . The sum over all particles’ weights is one. In SIR- k k k−1 k k−1 PF, the weights are computed according to the principle of p z | Z (19) k k−1 importance sampling where the so-called proposal density p z | x p x | Z k k k k−1 = . is chosen to be p(x | x = x ), and with resampling k k−1 k−1 p z | Z k k−1 at every time step. For N →∞, the approximate posterior approaches the true PDF. The term p(z | x ) = p(z | s =  s ) follows from (4)and k k k c,k c,k The key step, in which the measurement for instance k represents the probability of the measured channel output is incorporated, is in the calculation of the weight w which (often referred to as the likelihood value), conditioned on a k for the SIR-PF can be shown to be the likelihood function: certain configuration of channel parameters at the same time p(z | x ). The characterization of the channel process enters step k. The denominator of (19)doesnot depend on x ,and k in the algorithm, when at each time instance k, the state so it can be computed by integrating the numerator of (19) of each particle x is drawn randomly from the proposal over the entire range of x (normalization). k To summarize so far, the entire process of prediction distribution, that is, from p(x | x ). k−1 and update can be carried out recursively to calculate the posterior PDF (19) sequentially, based on an initial value of 4.3. Exploiting linear substructures p(x | z ) = p(x ). The evaluation of the likelihood function 0 0 0 If there exist linear substructures in the model, it is possible p(z | x ) is the essence of the update step. Similarly, k k maximizing this likelihood function (i.e., ML estimation) to reduce the computational complexity of the filter by means of marginalization over the linear state variables [24], would be equivalent to maximizing p(x | Z ) only in the k k also known as Rao-Blackwellization [25]. In a marginalized case that the prior PDF p(x | Z )doesnot depend on k k−1 Z and when all values of x are a priori equally likely. filter, particles are still used to estimate the nonlinear states, k−1, k while for each of the particles the linear states can be Since these conditions are not met, evaluation of p(x | Z ) k k estimated analytically. In our case, since the measurement entails all the above steps. z is a linear function of the complex amplitudes a , the k k likelihood function can be factorized as 4.2. Sequential estimation using particle filters p z | s = p z |τ , e , a = f z , τ , e ·g z , τ , e , a , k c,k k k k k k k k k k k k The optimal estimation algorithm relies on evaluating the (21) integral (18), whichisusually averydifficult task, except for certain additional restrictions imposed on the model and the where the function g (z , τ , e , a ) is Gaussian with respect k k k k noise process. So, very often a suboptimal realization of a to a ,and f (z , τ , e ) is proportional to p(z |  s ). k k k k k c,k Bayesian estimator has to be chosen for implementation. In Marginalization of (21) over the linear state variables gives this paper, we use a sequential Monte Carlo (SMC) filter, in particular a sampling importance resampling particle filter p z | τ , e , a da = f z , τ , e ∝ p z |  s . k k k k k k k k k c,k (SIR-PF) according to [23]. In this algorithm, the posterior (22) density at step k is represented as a sum and is specified by a set of N particles: Assuming that the amplitudes are block-wise independent (Q = 1), it follows that the weights of the SIR particle filter j j j j are equal to p(z | x ) = p(z | s =  s ), whichisgiven by k k c,k p x | Z ≈ w ·δ x − x , (20) k c,k k k k k k (4). j=1 Michael Lentmaier et al. 7 If there exists a temporal correlation of the amplitudes, Note that the model implicitly represents the number of an optimal marginalized filter requires the implementation paths of a separate Kalman filter for each of the particles. As we do not want to increase the complexity per particle, we take the N = e (28) m,k i,k correlation into account by adjustment of Q in (13). i=1 as a time-variant parameter. 4.4. Choice of appropriate channel process When applied to our particle filtering algorithm, drawing To exploit the advantages of sequential estimation for from the proposal density is simple. Each particle stores our task of multipath mitigation/estimation, we must be the above-channel parameters of the model, and then the able to describe the actual channel characteristics (channel new state of each particle is drawn randomly from p(x | parameters) so that these are captured by p(x | x ). k k−1 x ) which corresponds to drawing values for n and n as α τ k−1 In other words, the model must be a first-order Markov well as propagating the “on”/“off ” Markov model, and then model, and all transition probabilities must be known. In our updating the channel parameters for k according to (23)– approach, we approximate the channel as follows. (26). (i) The channel is totally characterized by a direct path (index i = 1) and at most N − 1echoes, 4.5. Practical issues (ii) each path has complex amplitude a and relative i,k 4.5.1. Model matching delay and Δτ = τ − τ ; where echoes are i,k i,k 1,k constrained to have delay τ ≥ τ , that is, Δτ ≥ 0, i,k 1,k i,k It is important to point out that a sequential estimator is only as good as its state transition model matches the real (iii) the different path delays follow the process: world situation. The state model needs to capture all relevant τ = τ + α ·Δt + n , 1,k 1,k−1 1,k−1 τ hidden states with memory and needs to correctly model (23) their dependencies, while adhering to the first order Markov Δτ = Δτ + α ·Δt + n , i,k i,k−1 i,k−1 τ condition. Furthermore, any memory of the measurement noise affecting the likelihood function p(z | x )mustbe k k (iv) each parameter α that specifies the speed of the i,k explicitly contained as additional states of the model x,so change of the path delay follows its own process: that the measurement noise is i.i.d. The channel state model is motivated by channel α = 1 − ·α + n , (24) i,k i,k−1 α 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 (v) the magnitudes and phases of the individual paths, characterize the channel for signal level simulations and represented by the complex amplitudes a , are elimi- i,k receiver evaluation comes close to our task of building a nated according to (12)and (14) for the computation first-order Markov process for sequential estimation. For of the likelihood (4), particle filtering, the model needs to satisfy the condition (vi) each path is either “on” or “off,” as defined by channel that one can draw states with relatively low computational parameter e ∈{1 ≡ “on, 0 ≡ “off }, i,k complexity. Adapting the model structure and the model (vii) where e follows a simple two-state Markov process i,k parameters to the real channel environment is a task for with a-symmetric crossover and same-state probabil- current and future work. It may even be possible to envisage ities: hierarchical models in which the selection of the current model itself follows a process. In this case, a sequential p e = 0 | e = 1 = p , (25) i,k i,k−1 onoff estimator will automatically choose the correct weighting of these models according to their ability to fit the received p e = 1 | e = 0 = p . (26) i,k i,k−1 offon signal. The model implicitly incorporates three i.i.d noise sources: Gaussian n and n as well as the noise process driving the τ α 4.5.2. Integration into a receiver state changes for e . These sources provide the randomness i,k For receiver integration, the computational complexity of of the model. The parameter K defines how quickly the the filtering algorithm is crucial. From a theoretical point speed of path delays can change (for a given variance of n ). of view, it is desirable to run the sequential filter clocked Finally, Δt is the time between instances k − 1and k.We corresponding to the coherent integration period of the assume all model parameters (i.e., K , Δt, noise variances, and receiver and with a very large number of particles. From the the “on”/“off ” Markov model) to be independent of k (see practical point of view, however, it is desirable to reduce the Figure 5). The hidden channel state vector x can therefore sequential filter rate to the navigation rate and to minimize be represented as the number of particles. Existing ML approaches can help τ , Δτ ,... , Δτ , α ,... , α , e ,... , e . here to achieve a flexible complexity/performance tradeoff, 1,k 2,k Nm ,k 1,k Nm ,k i,k Nm ,k as strategies already developed to extend the observation (27) 8 International Journal of Navigation and Observation τ τ 1,k−2 1,k−1 1,k α α α Direct path 1,k−2 1,k−1 1,k e e e 1,k−2 1,k−1 1,k Δτ Δτ 2,k−1 Δτ 2,k−2 2,k First echo α α α 2,k−2 2,k−1 2,k e e e 2,k−2 2,k−1 2,k . . . . . . ··· ··· . . . State x State x State x k−2 k−1 k Δτ Δτ Δτ N ,k−2 N ,k−1 N ,k m m m α α α N ,k−2 N ,k−1 N ,k N − 1’th echo m m m m N ,k−2 e e m N ,k−1 N ,k m m Figure 5: The Markov process chosen in this paper to model the channel with N paths. Dotted arrows—shown only for a small subset of the transitions—indicate the constraint that Δτ ≥ 0. For an explanation of the terms, see the text leading up to (27). i,k 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,wehavecarried 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πΔτ f with i,k i,k c f = 1575.42 MHz being the frequency of the L1 carrier. The 0 0 10 20 30 40 50 60 70 Bayesian estimator uses a time interval of Δt = 1 millisecond Relative delay (m) corresponding to the duration of a codeword. The amplitude averaging coefficient is set to Q = 10, and signal compression SIR-PF, path activitiy tracking DLL narrow correlator is applied with N = 41 (code-matched correlators) N = SIR-PF, single path model DLL strobe correlator cc pc 2 2 SIR-PF, two path model 25. Thechannelparameters σ , σ , K,and p(e | e )are i,k i,k−1 i,τ i,α selected to fit the statistics of a real channel according to [21]. Figure 6: Static scenario: performance of SIR-PF for BPSK The SIR-PF uses the minimum mean square error (MMSE) modulation as function of relative multipath delay for different path criterion to estimate the parameters x from the posterior. models. 5.1. Static multipath scenario The capability of multipath mitigation techniques is com- and with Strobe Correlator is also shown. The simulated monly assessed by showing the systematic error due to signal corresponds to a GPS L1 signal with c(t) being a Gold a single multipath replica plotted as a function of the code of length 1023 that is modulated on a bandlimited relative multipath delay in a static channel scenario. In rectangular pulse. The chip rate is 1.023 Mchips/s so that Figure 6, the root mean square error (RMSE) is shown the duration of the codeword is 1millisecond. The one- for the proposed sequential estimator, implemented as a sided bandwidth of the resulting navigation signal is 5 MHz. SIR-PF with N = 2000 particles. For comparison, the Estimators with fixed two-path model or fixed single path performance of conventional DLLs with Narrow Correlator model are also shown for comparison with the implicit path RMSE (m) Michael Lentmaier et al. 9 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 20 40 60 80 100 120 0 10 20 30 40 50 60 70 Time (s) Relative delay (m) Figure 9: Considered dynamic multipath scenario: pseudoranges Figure 7: Static scenario: average probability of a two-path model over time of the direct path (continuous line) and the temporarily for the estimator with path activity tracking. present echoes. −5 20 40 60 80 100 120 500 1000 1500 2000 2500 3000 Time (s) Number of particles Figure 10: Dynamic scenario: performance of the DLL with narrow Figure 8: Static scenario: RMSE performance as function of correlator for BPSK modulation (upper/grey) and BOC(1,1) modu- number of particles N for BPSK modulation. lation (lower/black). The segment marked with a dashed box shows a situation where BOC(1,1) is clearly superior to BPSK. activity tracking. The performance of a single path estimator −8 is comparable to that of a DLL with infinitesimal correlator σ = 10 ,and p = p = 0.0001 were chosen nτ onoff offon spacing and shows a considerable bias over a large delay to resemble a typical urban satellite navigation channel range. The estimator with fixed two-path model successfully environment [26]. mitigates the multipath bias for delays greater than 30 m. Figure 10 shows for this multipath scenario at a carrier- However, for smaller delays it shows an increasing variance to-noise ratio of C/N = 45 dB-Hz the BPSK and the and is outperformed by the single path estimator. The BOC(1,1) performance of a DLL with a narrow correlator −7 estimator with path activity tracking is capable of combining spacing of 10 seconds, corresponding to a tenth of a the advantages of both models. From the posterior, it code chip. A DLL loop bandwidth of 1 Hz was selected, is possible to calculate the estimated average probability which appeared to deliver the best DLL performance for this P(N = 2 | Z ) of a two-path model, which is shown in m,k k channel. For a fair comparison of the different modulation Figure 7, and indicates the transition between the models: techniques, both signals were generated with the same C/A for small delays the two paths essentially merge to a single code sequence of length 1023, the same receiver bandwidth one. Note that in these simulations, the model parameters of 5 MHz (one-sided), a sampling rate of 1/T = 10.23 MHz, of the sequential estimator are still the ones designed for the and a coherent integration time of LT = 1 millisecond. dynamic channel and not optimal for this static scenario. The Although the results show an improvement for the BOC(1,1) RMSE as function of the number of particles N is shown modulation, there still remains a considerable error due to in Figure 8 for a relative multipath delay of 14 m, which the multipath. This behavior is confirmed by a comparison corresponds to the worst case in Figure 6. 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 5.2. Dynamic scenario by the dashed box in Figure 10), more sophisticated multi- The dynamic multipath channel scenario with up to N = 3 path mitigating techniques are required to reduce the errors paths, used in the following simulations and depicted in due to short range multipath. Figure 9, has been generated according to the movement The simulation results for the particle filter-based estima- −10 model, whereby the parameters K = 25000, σ = 10 , tor with N = 20000 particles are given in Figures 12(a) and n s Probability of two path model RMSE (m) Pseudorange (m) Error (m) 10 International Journal of Navigation and Observation −5 −10 0 50 100 150 200 250 300 350 0.5 1 1.5 2 2.5 3 ×10 Relative delay (m) Number of particles Figure 11: Multipath error envelope of the DLL with narrow cor- Figure 13: Dynamic scenario: RMSE performance as function of relator for BPSK modulation (dashed) and BOC(1,1) modulation number of particles N for BPSK modulation. (inner/solid). 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 effi- cient likelihood computation in combination with a particle −5 filter realization of the prediction and update recursion. The considered movement model has been adapted to dynamic −10 multipath scenarios and incorporates the number of echoes 0 20 40 60 80 100 120 as a time-variant hidden channel state variable that is tracked Time (s) together with the other parameters in a probabilistic fashion. (a) BPSK modulation 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 −5 as a simple form of Rao-Blackwellization to eliminate the complex amplitudes from the state space. Further work will 0 20 40 60 80 100 120 concentrate on additional complexity reduction techniques Time (s) such as more suitable proposal functions or particle filtering (b) BOC(1,1) modulation algorithms such as the auxiliary particle filter that are possibly more efficient with respect to the number of Figure 12: Dynamic scenario: Performance of the sequential particles when applied to our problem domain. estimator with particle filtering (black) compared to the DLL (grey). The sequential estimator is clearly superior for multipath mitigation. REFERENCES [1] A. H. Jazwinski, Stochastic Processes and Filtering Theory, Academic Press, New York, NY, USA, 1970. [2] A. Doucet, N. de Freitas, and N. Gordon, Eds., Sequential 12(b) for the same channel as in Figure 10. They show that Monte Carlo Methods in Practice, Springer, New York, NY, the performance could be improved substantially, resulting USA, 2001. in a reduction of the root mean square (RMS) error from [3] B. Ristic, S. Arulampalam, and N. Gordon, Beyond the Kalman 3.77 m to 0.769 m for BPSK modulation and from 2.61 m to Filter: Particle Filters for Tracking Applications,ArtechHouse, 0.694 m for BOC(1,1) modulation, respectively. Boston, Mass, USA, 2004. The RMSE as function of the number of particles N is s [4] A.J.Van Dierendonck, P. C. Fenton,and T. Ford,“Theory and shown in Figure 13. performance of narrow correlator spacing in a GPS receiver,” Error (m) Error envelope (m) Error (m) RMSE (m) Michael Lentmaier et al. 11 in Proceedings of the ION National Technical Meeting,San [18] T. Bertozzi, D. L. Ruyet, C. Panazio, and H. V. Thien, “Channel Diego, Calif, USA, January 1992. tracking using particle filtering in unresolvable multipath environments,” EURASIP Journal on Applied Signal Processing, [5] L. Garin, F. van Diggelen, and J.-M. Rousseau, “Strobe and vol. 2004, no. 15, pp. 2328–2338, 2004. edge correlator multipath mitigation for code,” in Proceedings of the 9th International Technical Meeting of the Satellite [19] S. M. Kay, Fundamentals of Statistical Signal Processing: Division of the Institute of Navigation (ION GPS ’96), pp. 657– Estimation Theory, Prentice Hall, Upper Saddle River, NJ, 664, Kansas City, Mo, USA, September 1996. USA, 1993. [6] G. MacGraw and M. Brasch, “GNSS multipath mitigation [20] L. R. Weill, “Achieving theoretical bounds for receiver-based using gated and high resolution correlator concepts,” in multipath mitigation using galileo OS signals,” in Proceedings Proceedings of the ION National Technical Meeting,San Diego, of the 19th International Technical Meeting of the Satellite Calif, USA, January 1999. Division of the Institute of Navigation (ION GNSS ’06),pp. 1035–1047, Fort Worth, Tex, USA, September 2006. [7] J.Jones,P.C.Fenton, andB.Smith,“Theory andperformance of the pulse aperture correlator,” NovAtel Technical Report, [21] A. Steingass and A. Lehner, “Measuring the navigation NovAtel, Calgary, Alberta, Canada, 2004. multipath channel—a statistical analysis,” in Proceedings of the 17th International Technical Meeting of the Satellite Division [8] R. D. J. van Nee, J. Siereveld, P. C. Fenton, and B. R. Townsend, of the Institute of Navigation (ION GNSS ’04), pp. 1157–1164, “The Multipath estimating delay lock loop: approaching Long Beach, Calif, USA, September 2004. theoretical accuracy limits,” in Proceedings of the IEEE Position Location and Navigation Symposium (PLANS ’94), pp. 246– [22] P. Stoica, Y. Selen, ´ and J. Li, “On information criteria and 251, Las Vegas, Nev, USA, April 1994. the generalized likelihood ratio test of model order selection,” IEEE Signal Processing Letters, vol. 11, no. 10, pp. 794–797, [9] L. R. Weill, “Achieving theoretical accuracy limits for pseudo- ranging in the presence of multipath,” in Proceedings of the 8th International Technical Meeting of the Satellite Division of the [23] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A Institute of Navigation (ION GPS ’95), pp. 1521–1530, Palm tutorial on particle filters for online nonlinear/non-Gaussian Springs, Calif, USA, September 1995. Bayesian tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 174–188, 2002. [10] J. Selva Vera, “Complexity reduction in the parametric estimation of superimposed signal replicas,” Signal Processing, [24] T. Schon, ¨ F. Gustafsson, and P.-J. Nordlund, “Marginalized vol. 84, no. 12, pp. 2325–2343, 2004. particle filters for mixed linear/nonlinear state-space models,” IEEE Transactions on Signal Processing, vol. 53, no. 7, pp. 2279– [11] J. Selva Vera, “Efficient multipath mitigation in naviga- 2289, 2005. tion systems,” Ph.D. dissertation, Universitat Politecnica de Catalunya, Barcelona, Spain, February 2004. [25] A. Doucet, N. de Freitas, K. Murphy, and S. Russell, “Rao- blackwellised particle filtering for dynamic Bayesian net- [12] F. Antreich, O. Esbr´ ı-Rodr´ ıguez, J. A. Nossek, and W. Utchick, works,” in Proceedings of the 16th Annual Conference on “Estimation of synchronization parameters using SAGE in Uncertainty in Artificial Intelligence (UAI ’00), pp. 176–183, a GNSS-receiver,” in Proceedings of the 18th International San Francisco, Calif, USA, June-July 2000. Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’05), pp. 2124–2131, Long Beach, [26] A. Steingass and A. Lehner, “Land mobile satellite na- Calif, USA, September 2005. vigation—characteristics of the multipath channel,” in Pro- ceedings of the 16th International Technical Meeting of the [13] P. C. Fenton and J. Jones, “The theory and performance of Satellite Division of the Institute of Navigation (ION GNSS ’03), NovAtel Inc.’s vision correlator,” in Proceedings of the 18th Portland, Ore, USA, September 2003. International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’05), pp. 2178–2186, Long Beach, Calif, USA, September 2005. [14] P. Closas, C. Fernandez-P ´ rades, J. A. Fernandez-R ´ ubio, and A. Ram´ ırez-Gonzalez, ´ “Multipath mitigation using particle filtering,” in Proceedings of the 19th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’06), pp. 1733–1740, Fort Worth, Tex, USA, September 2006. [15] M. Lentmaier and B. Krach, “Maximum likelihood multipath estimation in comparison with conventional delay lock loops,” in Proceedings of the 19th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS ’06), pp. 1741–1751, Fort Worth, Tex, USA, September 2006. [16] B. Krach and M. Lentmaier, “Efficient soft-output GNSS signal parameter estimation using signal compression techniques,” in Proceedings of the 3rd ESA Workshop on Satellite Navigation User Equipment Technologies (NAVITECH ’06), Noordwijk, The Netherlands, December 2006. [17] E. Punskaya, A. Doucet, and W. J. Fitzgerald, “Particle filtering for joint symbol and code delay estimation in DS spread spec- trum systems in multipath environment,” EURASIP Journal on Applied Signal Processing, vol. 2004, no. 15, pp. 2306–2314, 2004. International Journal of Rotating Machinery International Journal of Journal of The Scientific Journal of Distributed Engineering World Journal Sensors Sensor Networks Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Volume 2014 Journal of Control Science and Engineering Advances in Civil Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 Submit your manuscripts at http://www.hindawi.com Journal of Journal of Electrical and Computer Robotics Engineering Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 VLSI Design Advances in OptoElectronics International Journal of Modelling & Aerospace International Journal of Simulation Navigation and in Engineering Engineering Observation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2010 Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com http://www.hindawi.com Volume 2014 International Journal of Active and Passive International Journal of Antennas and Advances in Chemical Engineering Propagation Electronic Components Shock and Vibration Acoustics and Vibration Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation Hindawi Publishing Corporation http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014 http://www.hindawi.com Volume 2014

Journal

International Journal of Navigation and ObservationHindawi Publishing Corporation

Published: Apr 24, 2008

References