# Bayesian Time Delay Estimation of GNSS Signals in  Dynamic Multipath Environments

Authors: Michael Lentmaier; Bernhard Krach; Patrick Robertson
Journal: International Journal of Navigation and Observation
 (2008)

---

## Abstract

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


---

## Body

## 1. Introduction

A  major error source within global
navigation satellite systems (GNSSs) comes from
multipath, the reception of additional signal replica due to reflections, which
introduce a bias into the estimate of the delay lock loop (DLL) of a
conventional navigation receiver. For efficient removal of this bias, it is
possible to formulate advanced maximum likelihood (ML) estimators that
incorporate the echoes into the signal model and are capable of achieving the
theoretical limits given by the Cramér Rao bound. The drawback of ML estimator
techniques is that the parameters are assumed to be constant during the time of
observation. Independent estimates are obtained for successive observation
intervals, whose length has to be adapted to the dynamics of the channel.In this paper, we consider the important practical
case of dynamic channel scenarios and assess how the time-delay estimation can
be improved if information is available about the temporal evolution of the
channel parameters. Our approach is based on Bayesian filtering, the optimal
and well-known framework to address such dynamic state estimation problems
[1]. Sequential Monte
Carlo (SMC) methods [2, 3] are used for computing the posterior probability
density functions (PDFs) of the signal parameters.


## 2. A Comparison of Various Multipath Mitigation Approaches

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


## 3. Signal Model

Assume that the  complex valued baseband-equivalent received signal is equal to(1)z(t)=∑i=1Nmei(t)⋅ai(t)⋅[c(t)∗g(t−τi(t))]+n(t), where c(t) is a
delta-train code sequence that is modulated on a pulse g(t), Nm is the maximum
number of allowed paths reaching the receiver (to restrict the modeling
complexity), ei(t) is a binary
function that controls the activity of the i′th path, and ai(t) and τi(t) are their
individual complex amplitudes and time delays, respectively. The signal is
disturbed by additive white Gaussian noise n(t). Grouping blocks of L samples at
times (m+kL)Ts, m=0,…,L−1, together into vectors zk, k=0,1,…,n whilst assuming
the parameter functions ei(t), ai(t), and τi(t) being constant
within the corresponding time interval and equal to ei,k, ai,k, and τi,k, this can be rewritten as(2)zk=CG(τk)Ekak+nk=^sk+nk.In the compact form on the right
hand side, the samples of the delayed pulses g(τi,k) are stacked
together as columns of the matrix G(τk)=[g(τ1,k),…,g(τNm,k)], C is a matrix
representing the convolution with the code, and the delays and amplitudes are
collected in the vectors τk=[τ1,k,…,τNm,k]T and ak=[a1,k,…,aNm,k]T, respectively.
Furthermore, for concise notation we use Ek=diag[ek] whilst the
elements of the vector ek=[e1,k,…,eNm,k]T, ei,k∈[0,1], determine whether the i′th path is
active or not by being either ei,k=1 corresponding
to an active path or ei,k=0 for a path that
is currently not active. The term sk denotes the
signal hypothesis and is completely determined by the channel parameters τk, ak, and ek. Using
(2), we can write the associated likelihood
function as(3)p(zk∣sk)=1(2π)Lσ2L⋅exp[−12σ2(zk−sk)H(zk−sk)].The likelihood function will
play a central role in the algorithms discussed in this paper; its purpose is
to quantify the conditional probability of the received signal conditioned on
the unknown signal (specifically the channel parameters).


## 3.1. Efficient Likelihood Computation

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


## 3.1.1. Data Compression

As explained above, the large vector containing the
received signal sampleszk is linearly
transformed into a vector zc,k of much smaller
size. Following this approach, the likelihood according to (2) can be rewritten
as(4)p(zk∣sk)=1(2π)Lσ2Lexp[−zkHzk2σ2]⋅exp[ℜ{zkHQcQcHsk}σ2−skHQcQcHsk2σ2]=1(2π)Lσ2Lexp[−zkHzk2σ2]⋅exp[ℜ{zc,kHsc,k}σ2−sc,kHsc,k2σ2],with the compressed received
vector zc,k and the
compressed signal hypothesis sc,k:(5)zc,k=QcHzk,sc,k=QcHsk,and the orthonormal compression
matrix Qc, which needs to fulfill(6)QcQcH≈I,QcHQc≈I,to minimize the compression
loss. According to [11], the compression can be two-fold so that we can
factorize(7)Qc=QccQpcinto a canonical component
decomposition, given by an L×Ncc matrix Qcc, and a principal component decomposition,
given by an Ncc×Npc matrix Qpc. In [11], two choices for Qcc are
proposed(8)Qcc={CG(τb)Rcc−1Signal matched,C(τb)Rcc−1Code matched,where the elements of the vector τb define the
positions of the individual correlators. The noise-free outputs of the
corresponding correlator banks are illustrated in Figure 2.  To decorrelate the bank
outputs (CG(τb))Hy and C(τb)Hy, as mentioned above, the whitening matrix Rcc can be obtained
from a QR decomposition of CG(τb) and C(τb), respectively.
Apart from practical implementation issues, both correlation methods given by
(8) are equivalent from a conceptual point of view. For details on the
compression through Qpc, the reader is
referred to [11].Output of the
canonical component type correlator banksCG(τb) and C(τb) for BPSK and
BOC(1,1).
(a)
BPSK signal matched(b)
BPSK code matched(c)
BOC(1,1) signal matched(d)
BOC(1,1) code matched


## 3.1.2. Interpolation

In order to
compute (4) independently of the sampling grid,
advantage can be made of interpolation techniques. Using the discrete Fourier
transformation (DFT), with Ψ being the DFT
matrix and Ψ−1 being its
inverse counterpart (IDFT), we get(9)sc,k=QcHCΨ−1diag[Ψg(0)]Ω(τk)Ekak=^MscΩ(τk)Ekak,with Ω(τk) being a matrix
of column-wise stacked vectors with Vandermonde structure [10, 11], such that the element at
row p and column q computes
with(10)ℜ{[Ω(τk)]p,q}=cos(2π(p−1)τq,kNgTs),ℑ{[Ω(τk)]p,q}=−sin(2π(p−1)τq,kNgTs).Ng is the length
of the pulse g in samples. The
advantage of the interpolation is that it can take place in the reduced space.
The most costly computations in (9) can be carried out in precalculations as
the matrix Msc, whose row dimension corresponds to the dimension of
the reduced space and whose column dimension is Ng, is constant.


## 3.1.3. Amplitude Elimination

In a further
step, we reduce the number of parameters by optimizing (4) for a given set of τk and ek with respect to
the complex amplitudes ak, which can be achieved through a closed-form
solution. Using(11)Sc,k=MscΩ(τk)Ekand obtaining Sc,k+ by removing
zero columns from Sc,k, one yields the
corresponding amplitude values of the active paths:(12)a^k+=(Sc,k+HSc,k+)−1Sc,k+Hz¯c,k.As we have introduced a
potential source of performance loss by eliminating the amplitudes and thus
practically are disregarding their temporal correlation, we propose to optimize
(4) using(13)z¯c,k=1Q⋅∑l=0Q−1zc,k−l,with the adjustable averaging
coefficient Q. When evaluating (4), we substitute sc,k by(14)s˜c,k=Sc,ka^k,where the elements of the vector a^k that are
indicated to have an active path (ak,i:i→ek,i=1) are set equal
to the corresponding elements of a^k+. All other elements (ak,i:i→ek,i=0) can be set
arbitrarily as their influence is masked by the zero elements of ek.


## 3.2. Review of the ML Concept

The concept of
ML multipath estimation has drawn substantial research interest since the first
approach was proposed in [8]. Despite being treated differently in various
publications, the objective is the same for all ML approaches, namely to find
the signal parameters that maximize (3) for a given observation zk:(15)s^k=argmaxsk{p(zk∣sk)}.The signal parameters are
thereby assumed as being constant throughout the observation period k. Different maximization strategies exist, which
basically characterize the different approaches. Despite offering great
advantages for theoretical analysis, the practical advantage of the generic ML
concept is questionable due to a number of serious drawbacks.(i)
The ML
estimator assumes that the channel is static for the observation period and is
not able to exploit its temporal correlation throughout the sequencek=1,…,n. Measured channel scenarios have shown significant
temporal correlation [21].(ii)
Despite being
of great interest in practice, the estimation of the number of received paths
is often not addressed. The crucial problem here is to correctly estimate the
current number of paths to avoid over determination, since an overdetermined
estimator will tend to use the additional degrees of freedom to match the noise
by introducing erroneous paths. There exist various techniques based on model
selection that can be employed to estimate the number of paths [22] but they suffer from the problem of having to
dynamically adjust the decision thresholds. Typically, only a single hypothesis
is tracked, which in practice causes error event propagation.(iii)
The ML
estimator does only provide the most likely parameter set for the given
observation. No reliability information about the estimates is provided.
Consequently, ambiguities and multiple modes of the likelihood are not
preserved by the estimator.ML estimators require that the estimated parameters
remain constant during the observation period. Due to data modulation and phase
variations in practice, this period, which is often referred to as the coherent
integration time, is limited to a range of 1 millisecond–20 milliseconds. To
reach a sufficient noise performance with a ML estimator in practice, it is
required to extend its observation interval approximately to the equivalent
averaging time of a conventional tracking loop, which is commonly in the order
of several hundred coherent integration periods. To overcome this problem, the
observations are forced to be quasicoherent by aiding the ML estimator with a
phase locked loop (PLL) and a data removal mechanism [8].


## 4. Sequential Estimation

## 4.1. Optimal Solution

In Section3,
we have established the models of the underlying time variant processes. The
problem of multipath mitigation now becomes one of sequential estimation of
a hidden Markov process. We want to estimate the unknown channel parameters
based on an evolving sequence of received noisy channel outputs zk. The channel process for each range of a satellite
navigation system can be modeled as a first-order Markov process if future
channel parameters given the present state of the channel and all its past
states depend only on the present channel state (and not on any past states).
We also assume that the noise affecting successive channel outputs is
independent of the past noise values; so each channel observation depends
only on the present channel state.Intuitively, we are exploiting not only the channel
observations to estimate the hidden channel parameters (via the likelihood
function), but we are also exploiting our prior knowledge about the statistical
dependencies between successive sets of channel parameters. We know from
channel measurements that channel parameters are time varying but not
independent from one time instance to the next; for example, an echo usually experiences
a “life-cycle” from its first occurrence, then a more or less gradual
change in its delay and phase over time, until it disappears [21]. These measurements also
show that the common channel models considered in communication systems
[17, 18] are not adequately
reflecting the properties that are crucial for high-resolution signal delay
estimation as required in navigation systems.Now that our major assumptions have been established,
we may apply the concept ofsequential Bayesian estimation. The reader
is referred to [23]
which gives a derivation of the general framework for optimal estimation of
temporally evolving (Markovian) parameters by means of inference, and we have
chosen similar notation. The entire history of observations (over the temporal
index k) can be
written as(16)Zk=^{zk′,k′=1,…,k}.Similarly, we denote the
sequence of parameters of our hidden Markovian process by(17)Xk=^{xk′,k′=1,…,k}.Here, xk′ represents the
characterization of the hidden channel state, including the parameters that
specify the signal hypothesis sk′ given in
(2).
Our goal is to determine the posterior PDF of every possible channel
characterization given all channel observations: p(xk∣Zk) (see Figure  3).
Once we have evaluated this posterior PDF, we can either determine that channel
configuration that maximizes it—the so-called maximum a posteriori (MAP) estimate;
or we can choose the expectation—equivalent to the minimum mean square error
(MMSE) estimate. In addition, the posterior distribution itself contains all
uncertainty about the current range and is thus the optimal measure to perform
sensor data fusion in an overall positioning system.Figure 3
Illustration of the hidden Markov estimation process
for three time instances. Our channel measurements are the sequence{zi,i=1,…,k}, and the channel parameters to be estimated are {xi,i=1,…,k}.It can be shown that the sequential estimation
algorithm is recursive, as it uses the posterior PDF computed for time instancek−1 to compute the posterior
PDF for instance k (see Figure  4).
For a given posterior PDF at time instance k−1, p(xk−1∣Zk−1), the prior PDF p(xk∣Zk−1) is calculated
in the so-called prediction step by applying the Chapman-Kolmogorov
equation:(18)p(xk∣Zk−1)=∫p(xk∣xk−1)p(xk−1∣Zk−1)dxk−1,with p(xk∣xk−1) being the state
transition PDF of the Markov process. In the update step, the new
posterior PDF for step k is obtained by
applying Bayes' rule to p(xk∣zk,Zk−1) yielding the
normalized product of the likelihood p(zk∣xk) and the prior
PDF:(19)p(xk∣Zk)=p(xk∣zk,Zk−1)=p(zk∣xk,Zk−1)p(xk∣Zk−1)p(zk∣Zk−1)=p(zk∣xk)p(xk∣Zk−1)p(zk∣Zk−1).The term p(zk∣xk)=p(zk∣sc,k=s˜c,k) follows from
(4) and represents the probability of the measured channel output (often
referred to as the likelihood value), conditioned on a certain configuration of
channel parameters at the same time step k. The denominator of (19) does not depend on xk, and so it can
be computed by integrating the numerator of (19) over the entire range of xk (normalization).Figure 4
Illustration of the recursive Bayesian estimator.To summarize so far, the entire process of prediction
and update can be carried out recursively to calculate the posterior PDF (19)
sequentially, based on an initial value of p(x0∣z0)=p(x0). The evaluation of the likelihood function p(zk∣xk) is the essence
of the update step. Similarly, maximizing this likelihood function (i.e., ML
estimation) would be equivalent to maximizing p(xk∣Zk) only in the
case that the prior PDF p(xk∣Zk−1) does not depend
on Zk−1, and when all
values of xk are a priori
equally likely. Since these conditions are not met, evaluation of p(xk∣Zk) entails all the
above steps.


## 4.2. Sequential Estimation Using Particle Filters

The optimal estimation algorithm relies on evaluating
the integral (18), which is usually a very difficult task, except for certain
additional restrictions imposed on the model and the noise process. So, very
often a suboptimal realization of a Bayesian estimator has to be chosen for
implementation. In this paper, we use a sequential Monte Carlo (SMC) filter, in
particular a sampling importance resampling particle filter (SIR-PF) according
to [23]. In this
algorithm, the posterior density at step k is represented
as a sum and is specified by a set of Np particles:(20)p(xk∣Zk)≈∑j=1Npwkj⋅δ(xk−xkj),where each particle with index j has a state xkj and has a weightwkj. The sum over all particles' weights is one. In
SIR-PF, the weights are computed according to the principle of importance
sampling where the so-called proposal density is chosen to be p(xk∣xk−1=xk−1j), and with resampling at every time step. For Np→∞, the approximate
posterior approaches the true PDF.The key step, in which themeasurement for
instance k is
incorporated, is in the calculation of the weight wkj which for the
SIR-PF can be shown to be the likelihood function: p(zk∣xkj). The characterization of the channel process enters in the algorithm, when at each time instance k, the state of each particle xki is drawn
randomly from the proposal distribution, that is, from p(xk∣xk−1j).


## 4.3. Exploiting Linear Substructures

If there exist
linear substructures in the model, it is possible to reduce the computational
complexity of the filter by means of marginalization over the linear state
variables [24], also
known as Rao-Blackwellization [25]. In a marginalized filter, particles are still used
to estimate the nonlinear states, while for each of the particles the linear
states can be estimated analytically. In our case, since the measurement zk is a linear
function of the complex amplitudes ak, the likelihood function can be factorized
as(21)p(zk∣sc,k)=p(zk∣τk,ek,ak)=f(zk,τk,ek)⋅g(zk,τk,ek,ak),where the function g(zk,τk,ek,ak) is Gaussian
with respect to ak, and f(zk,τk,ek) is proportional
to p(zk∣s˜c,k). Marginalization of
(21) over the linear state
variables gives(22)∫akp(zk∣τk,ek,ak)dak=f(zk,τk,ek)∝p(zk∣s˜c,k).Assuming that the amplitudes are
block-wise independent (Q=1), it follows
that the weights of the SIR particle filter are equal to p(zk∣xkj)=p(zk∣sc,k=s˜c,kj), which is given by
(4).If there exists a temporal correlation of the
amplitudes, an optimal marginalized filter requires the implementation of a
separate Kalman filter for each of the particles. As we do not want to increase
the complexity per particle, we take the correlation into account by adjustment
ofQ in (13).


## 4.4. Choice of Appropriate Channel Process

To exploit the advantages of sequential estimation for
our task of multipath mitigation/estimation, we must be able to describe the
actual channel characteristics (channel parameters) so that these are captured
byp(xk∣xk−1). In other words, the model must be a first-order
Markov model, and all transition probabilities must be known. In our approach,
we approximate the channel as follows.(i)
The channel is
totally characterized by a direct path (indexi=1) and at most Nm−1 echoes,(ii)
each path has
complex amplitudeai,k and relative
delay and >Δτi,k=τi,k−τ1,k; where echoes are constrained to have delay τi,k≥τ1,k, that is, >Δτi,k≥0,(iii)
the different
path delays follow the process:(23)τ1,k=τ1,k−1+α1,k−1⋅Δt+nτ,Δτi,k=Δτi,k−1+αi,k−1⋅Δt+nτ,(iv)
each parameterαi,k that specifies
the speed of the change of the path delay follows its own
process:(24)αi,k=(1−1K)⋅αi,k−1+nα,(v)
the magnitudes
and phases of the individual paths, represented by the complex amplitudesai,k, are eliminated according to (12) and (14) for the
computation of the likelihood (4),(vi)
each path is
either “on” or “off,” as defined by channel parameterei,k∈{1≡ “on,” 0≡ “off” },(vii)
whereei,k follows a
simple two-state Markov process with a-symmetric crossover and same-state
probabilities:(25)p(ei,k=0∣ei,k−1=1)=ponoff,(26)p(ei,k=1∣ei,k−1=0)=poffon. The model
implicitly incorporates three i.i.d noise sources: Gaussian nτ and nα as well as the
noise process driving the state changes for ei,k. These sources provide the randomness of the model.
The parameter K defines how
quickly the speed of path delays can change (for a given variance of nα). Finally, >Δt is the time
between instances k−1 and k. We assume all model parameters (i.e., K, >Δt, noise variances, and the “on”/“off”
Markov model) to be independent of k (see Figure 5).
The hidden channel state vector xk can therefore
be represented as(27)[τ1,k,>Δτ2,k,…,>ΔτNm,k,α1,k,…,αNm,k,ei,k,…,eNm,k]T.Note that the model implicitly
represents the number of pathsNm,k=^∑i=1Nmei,kas a time-variant parameter.Figure 5
The Markov
process chosen in this paper to model the channel withNm paths. Dotted
arrows—shown only for a small subset of the transitions—indicate the
constraint that >Δτi,k≥0. For an explanation of the terms, see the text
leading up to (27).When applied to our particle filtering algorithm,
drawing from the proposal density is simple. Each particle stores the
above-channel parameters of the model, and then the new state of each particle
is drawn randomly fromp(xk∣xk−1j) which
corresponds to drawing values for nα and nτ as well as
propagating the “on”/“off” Markov model, and then updating the
channel parameters for k according to
(23)–(26).


## 4.5. Practical Issues

## 4.5.1. Model Matching

It is important
to point out that a sequential estimator is only as good as its state
transition model matches the real world situation. The state model needs to
captureall relevant hidden states with memory and needs to correctly
model their dependencies, while adhering to the first order Markov condition.
Furthermore, any memory of the measurement noise affecting the likelihood
function p(zk∣xk) must be
explicitly contained as additional states of the model x, so that the measurement noise is i.i.d.The channel state model is motivated by channel
modeling work for multipath prone environments such as the urban satellite
navigation channel [21, 26]. In fact, the process of constructing a channel model
in order to characterize the channel for signal level simulations and receiver
evaluation comes close to our task of building a first-order Markov process for
sequential estimation. For particle filtering, the model needs to satisfy the
condition that one can draw states with relatively low computational
complexity. Adapting the model structure and the model parameters to the real
channel environment is a task for current and future work. It may even be
possible to envisage hierarchical models in which the selection of the current
model itself follows a process. In this case, a sequential estimator will
automatically choose the correct weighting of these models according to their
ability to fit the received signal.


## 4.5.2. Integration into a Receiver

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


## 5. Performance Evaluation

To demonstrate the capabilities of the SMC-based
Bayesian time delay estimator proposed in Section4, we have carried out
computer simulations for both static and dynamic channel environments. The
signal-to-multipath ratio (SMR) was chosen constant and equal to 6 dB, while the
relative phases are changing according to >Δφi,k=2π>Δτi,kfc with fc=1575.42 MHz being the
frequency of the L1 carrier. The Bayesian estimator uses a time interval of >Δt=1 millisecond
corresponding to the duration of a codeword. The amplitude averaging
coefficient is set to Q=10, and signal
compression is applied with Ncc=41 (code-matched
correlators) Npc=25. The channel parameters σi,τ2, σi,α2, K, and p(ei,k∣ei,k−1) are selected to
fit the statistics of a real channel according to [21]. The SIR-PF uses the
minimum mean square error (MMSE) criterion to estimate the parameters xk from the
posterior.


## 5.1. Static Multipath Scenario

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


## 5.2. Dynamic Scenario

The dynamic multipath channel scenario with up toNm=3 paths, used in
the following simulations and depicted in Figure 9, has been generated according
to the movement model, whereby the parameters K=25000, σnα=10−10, σnτ=10−8, and ponoff=poffon=0.0001 were chosen to
resemble a typical urban satellite navigation channel environment [26].Figure 9
Considered dynamic multipath scenario: pseudoranges
over time of the direct path (continuous line) and the temporarily present
echoes.Figure10 shows for this multipath scenario at a
carrier-to-noise ratio of C/N0=45 dB-Hz the BPSK
and the BOC-(1,1) performance of a DLL with a narrow correlator spacing of 10−7 seconds,
corresponding to a tenth of a code chip. A DLL loop bandwidth of 1 Hz was
selected, which appeared to deliver the best DLL performance for this channel.
For a fair comparison of the different modulation techniques, both signals were
generated with the same C/A code sequence of length 1023, the same receiver
bandwidth of 5 MHz (one-sided), a sampling rate of 1/Ts=10.23 MHz, and a
coherent integration time of LTs=1 millisecond.
Although the results show an improvement for the BOC-(1,1) modulation, there
still remains a considerable error due to the multipath. This behavior is
confirmed by a comparison with the multipath error envelope, depicted in Figure 11. While in some multipath delay regions, the error is reduced by the BOC-(1,1)
modulation (see, e.g, the segment marked by the dashed box in Figure 10), more sophisticated
multipath mitigating techniques are required to reduce the errors due to short
range multipath.Figure 10
Dynamic scenario: performance of the DLL with narrow
correlator for BPSK modulation (upper/grey) and BOC-(1,1) modulation (lower/black). The segment marked with a dashed box shows a situation where BOC(1,1)
is clearly superior to BPSK.Figure 11
Multipath error envelope of the DLL with narrow
correlator for BPSK modulation (dashed) and BOC(1,1) modulation (inner/solid).The simulation results for the particle filter-based
estimator withNs=20000 particles are
given in Figures 12(a)  and 12(b) for the same channel as in Figure 10. They show that the performance could be improved substantially, resulting
in a reduction of the root mean square (RMS) error from 3.77 m to 0.769 m for
BPSK modulation and from 2.61 m to 0.694 m for BOC(1,1) modulation,
respectively.Dynamic scenario: Performance of the sequential
estimator with particle filtering (black) compared to the DLL (grey). The
sequential estimator is clearly superior for multipath mitigation.
(a)
BPSK modulation(b)
BOC(1,1) modulationThe RMSE as function of the number of particlesNs is shown in
Figure 13.Figure 13
Dynamic scenario: RMSE performance as function of
number of particlesNs for BPSK
modulation.


## 6. Conclusions

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



---

*Source: 372651-2008-04-24.xml*