IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO.

9, SEPTEMBER 2017 7663

Exploiting WSSUS Multipath for GNSS Ranging

Christoph Enneking and Felix Antreich

Abstract—We propose a low-dimensional parameter estimation obstructed [1]. Therefore, multipath conditions are often clas-
scheme for satellite-based ranging in the presence of an attenuated sified into the three categories diffuse/specular, dynamic/static,
line-of-sight (LOS) signal and a time-varying number of multipath and LOS/non-LOS (NLOS), the latter referring to the absence of
signals. This is of interest for landmobile users of global navigation
satellite systems in diverse scattering environments with limited a dominant direct signal due to signal blockage. While echoes
satellite visibility (e.g., a car moving through an urban canyon). and NLOS are in fact different phenomena, they often occur
Conventional ranging techniques, which rely solely on estimation simultaneously in challenging physical environments, such as
of LOS parameters, ignore or suppress multipath and thus sacrifice urban canyons. Numerous approaches have been designed for
available information; on the other hand, exhaustive approaches satellite-based ranging in the presence of multipath [2]–[17], but
involving detection and estimation of each time-variant echo are
too cumbersome in terms of parameter dimension. We suggest an a lack of apprioriate solutions can be clearly identified for the
unconditional maximum likelihood estimator (U-MLE) for joint diffuse/dynamic/NLOS scenario. Specifically, we will argue in
estimation of LOS synchronization parameters and a small set of the following that the joint occurrence of an attenuated LOS sig-
multipath distribution parameters. Our key assumptions are Gaus- nal and many rapidly appearing and disappearing echoes poses
sian wide-sense stationary uncorrelated scattering (WSSUS) and a a great difficulty to state-of-the-art processing techniques. Inter-
known structure of the power delay profile. The proposed estimator
uses as input data only a few complex-valued correlator outputs, estingly, these are precisely the propagation conditions encoun-
which are provided by conventional navigation receivers, and its tered by landmobile users in low-elevation, urban and indoor
complexity is independent of the actual number of present echoes. scenarios [18]–[22], where there is a growing desire for naviga-
Simulations demonstrate that the U-MLE improves accuracy and tion with GNSS.
robustness of LOS synchronization for various model-matched and Multipath suppression approaches are widespread in prac-
-mismatched multipath settings, including standardized and ap-
proved urban scenarios. tice, as they are computationally cheap and work regardless of
whether the echoes are dynamic or static, diffuse or specular.
Index Terms—Dense multipath, Global positioning system These techniques aim at filtering the incoming signal in such a
(GPS), high-sensitivity GNSS, non-line-of-sight (NLOS).
way that the echoes do not enter into the statistic that is used
for estimation of the LOS time-delay. This can be achieved,
I. INTRODUCTION at least for a part of the echoes, with directional antennas [3],
subspace or beamforming methods [4], [5], or narrow correlator
RECISE ranging in multipath environments is widely
P agreed to be one of the major challenges for users of
global navigation satellite systems (GNSS). To determine its
spacing combined with a large pre-correlation bandwidth [6].
If long coherent observations of the LOS signal are available,
the receiver can also simply take the mean over all observa-
distance to a satellite, a receiver must estimate the time-delay
tions, so that noise and non-coherent multipath are averaged out
of the line-of-sight (LOS) signal with high accuracy. This es-
[2]. However, suppression methods are often unable to provide
timation can be deteriorated in terms of bias and variance by
ranging in NLOS conditions, as they sacrifice the information
multipath signals (“echoes”) interfering with the LOS signal
that multipath alone can provide for ranging. The only way to
if standard synchronization techniques are used [1], [2]. The
obtain an accurate LOS time-delay estimate in this case is to
particular impact of multipath propagation is known to depend
use extremely long observation times of ten seconds or more
on the intensity and number of echoes, on the relative receiver-
to recover the weak LOS signal, which has become known as
scatterer movement, and on whether or not the LOS signal is
high-sensitivity GNSS [7], [8].
For multipath estimation, the echoes are formulated as a part
Manuscript received August 10, 2016; revised December 31, 2016; ac- of the signal model. Unknown parameters are introduced, which
cepted March 1, 2017. Date of publication March 28, 2017; date of current
version September 15, 2017. This work was supported by the Coordenação can include the time-variant amplitude, carrier-phase and time-
de Aperfeiçoamento de Pessoal de Nı́vel Superior under the PVE Grant delay of each multipath arrival. Once the correct estimates are
88881.030392/2013-01. The review of this paper was coordinated by Dr. K. obtained, they can be used to trivially subtract the multipath
Yu. (Corresponding author: Christoph Enneking.)
C. Enneking is with the Institute for Communications and Naviga- components from the received signal [9]–[15], the rationale be-
tion, German Aerospace Center (DLR), Wessling 82234, Germany (e-mail: ing that echoes are a mere nuisance that should be cancelled
christoph.enneking@dlr.de). at some stage of the receiver. A much more favorable option
F. Antreich is with the Institute for Communications and Navigation, German
Aerospace Center (DLR), Wessling 82234, Germany, and also with the Depart- in NLOS conditions is to actually benefit from the echo esti-
ment of Teleinformatics Engineering, Federal University of Ceará, Fortaleza mates for ranging or even positioning, which was demonstrated
60020-180, Brazil (e-mail: antreich@ieee.org). in an inspiring work by Gentner and Jost [16]. In any case, the
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org. joint estimation of several paths requires optimization or inte-
Digital Object Identifier 10.1109/TVT.2017.2688724 gration in many dimensions. To solve such problems optimally
0018-9545 © 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications standards/publications/rights/index.html for more information.

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7664 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

in a maximum likelihood (ML) or Bayesian sense, algorithms beforehand, but a parametric model can be used (e.g. exponen-
based on space-alternating generalized expectation maximiza- tial). WSSUS channels are commonly used in terrestrial com-
tion (SAGE) [12], Newton’s method [13] or Sequential Monte munication systems [26]–[30], and are mentioned as a model
Carlo [14]–[16] have been proposed. The parameter dimension candidate for the satellite-to-earth channel in [21].
grows linearly with the number of paths, which is acceptable In this work, we propose a U-MLE for CRLB-efficient joint
for very few and persistent echoes (e.g. a single reflection from estimation of deterministic LOS synchronization parameters
the ground during a low-altitude flight), but leads to impracti- and statistical WSSUS multipath distribution parameters,
cal estimators if many time-variant echoes are received. Two including the PDP’s delay spread and the average multipath
other major difficulties of exhaustive multipath estimation are power. We believe that this is a significant contribution, because
the detection problems before and after the actual estimation. in comparable works the multipath distribution is either known
First, the number of present echoes can vary with time and a priori [25] or the important estimation of the LOS time-delay
has to be determined for each observation period (using infor- “presents a noticable gap with respect to its CRLB” [24].
mation criteria [23]); this is computationally costly, and erro- Rather than exploiting distinct echoes and estimating path
neous detection will lead to over-/underdetermined estimators parameters [16], we estimate distribution parameters and
[15]. Second, once the parameters of all paths are estimated, exploit an ensemble of echoes. Especially, the fact that rapidly
the receiver has to decide which path actually is the LOS; in appearing and disappearing echoes can be exploited for precise
case the multiple arriving wavefronts are highly correlated and ranging in case the LOS signal is blocked, provides a very
have similar amplitude, outliers in the decision on the LOS interesting perspective for vehicular applications in demanding
occur [12]. environments. Attaining the CRLB is achieved in our work by
A strategy that has the potential to overcome the above de- assuming a mobile/vehicular receiver which records sufficiently
scribed difficulties is multipath distribution estimation: rather many realizations of an ergodic multipath channel as time is
than trying to resolve individual echoes explicitly, statistical passing while the receiver is moving through a diverse scatter-
multipath parameters are estimated after a suitable random ing environment. This approach has reasonable computational
model has been selected. A rank-one random multipath model complexity and thus can be envisaged for mobile applications.
was proposed for GNSS ranging in [17], which decomposes the However, we do make the assumption that the receiver is able
echoes into a signal which is fully correlated with the LOS signal to acquire the LOS signal and to remove its Doppler frequency
and a delay/temporally white Gaussian noise process with un- precisely. This can be achieved by using high-sensitivity GNSS
known power. For channel sounding of a mobile radio channel techniques on acquisition level or with external velocity aiding.
with many diffuse echoes, it was suggested to assume uncorre- The paper is organized as follows. In Section II, we in-
lated scattering (US) and estimate a parameterized power delay troduce the system model in pre- and post-correlation do-
profile (PDP) with unknown delay spread [24]. Furthermore, main. Section III states the WSSUS random channel model. In
for an a priori known PDP, Bialer, Raphaeli and Weiss [25] Section IV, we derive the multipath and noise statistics of the
propose to exploit this knowledge for ultra-wideband (UWB) post-correlation signal. Section V deals with parameter esti-
ranging in case of NLOS conditions, the key assumptions being mation, where we present the U-MLE as well as a simple
US and that the LOS arrival is the earliest arrival. The main echo-suppressing estimator. Simulation results are discussed in
advantage of approaches like [17], [24], [25] is that the actual Section VI.
number and realizations of random echoes are not of interest for
the parameter estimation, other than in [9]–[16]. This leads to II. SYSTEM MODEL
a multipath model of reduced and fixed order without the need
For the signal model, we distinguish between pre- and post-
for echo detection, but possibly improving LOS synchronization
correlation domain. Prior to correlation, the LOS and multipath
and ranging accuracy.
signals arriving at the receiver are delayed versions of a spread-
Clearly, the choice of the assumed multipath distribution
ing waveform, affected by different absolute delays and Doppler
should depend on the application and the required accuracy.
shifts. On the other side, the despread post-correlation sam-
A rank-one model as simple as [17] cannot exploit far echoes
ples are compensated up to unknown small-scale delay/Doppler
uncorrelated with the LOS; the US model in [24] does not ac-
residuals and will serve as data for parameter estimation. We
count for echo dynamics and leads to a synchronization estimate
will also use approximations for a slowly time-variant channel
unable to attain the Cramér-Rao Lower Bound (CRLB); finally,
and narrow bandwidth GNSS signals.
the US approach in [25] relies on measuring the PDP a priori
by transmitting very short UWB pulses, which is not possible
A. Pre- and Post-correlation Signal Model
when the receiver can only use GNSS signals. A reasonable
generalization of these models is a wide-sense stationary uncor- Let us assume that the satellite transmits the bandpass signal
related scattering (WSSUS) channel with Gaussian distribution.  
ssat (t) ∝ Re ejω c t s(t) t ∈ R, (1)
This channel includes the dimensions time and delay, allowing
for echo dynamics and arbitrary excess delays, but can still be where s(t) is a known dataless spreading waveform and ωc
characterized with low complexity by its mean, time-correlation is the angular carrier frequency. Moreover, s(t) is considered
function, and PDP. In an unknown environment, the PDP must strictly bandlimited to one-sided bandwidth Bs  ωc /2π and
be considered unknown as GNSS offers no way to determine it periodic with period T , i.e., s(t) ≡ s(t + kT ) for all k ∈ Z. Let

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7665

Fig. 1. Geometry of receiver R, scatterer S  and satellite T .

the autocorrelation function of the signal be denoted by


1
rs (τ ) = s∗ (t)s(t − τ ) dt. (2)
T T
Owing to the T -periodicity of s(t), this integral can be taken
over any interval of length T .
The geometry of satellite, receiver, and any scatterer  is
illustrated in Fig. 1. The scaterring environment is assumed
static in 2-D, while the receiver is moving with constant velocity
v = |v |.
The angle ε0 refers to satellite elevation, while ϑ0 and ϑ
denote the azimuth of satellite and scatterer , respectively, with
respect to receiver heading. Further, let r0 be the range between
Fig. 2. Correlator stage for compression of the received signal y(t) to a
satellite and receiver at time t = 0. vector-sequence of correlator outputs [Z (1) [n], . . . , Z (Q ) [n]]T .
The pre-correlation received signal in complex baseband no-
tation is given in terms of the th path’s received power P ,
carrier-phase φ , Doppler shift ω and time-delay τ as where ρ > 0 is the geometric excess path length of the th
    path with respect to the direct path at time t = 0. The Doppler
ω
y(t) = P ejφ  ejω  t s t − τ + t + η(t), (3) term Cω is the same as for the LOS, naturally assuming that the
ωc satellite is very far away.
∈L(t)
Presuming that a priori estimates τ̂0 and ω̂0 for time-delay
where L(t) ⊂ N0 is a finite index set containing all active paths and Doppler frequency are available at time t = 0 (e.g. from
at time t. The LOS path  = 0 is always considered active to the receiver’s acquisition stage or external aiding), the LOS
serve as a reference, but may have quasi-zero amplitude. The dynamics (e.g. Doppler frequencies of kHz-magnitude due to
process η(t) is complex zero-mean white Gaussian noise with high satellite velocity) can be removed and time-delay uncer-
power spectral density N0 . tainty can be reduced with some acquisition resolution TA > 0.
The LOS parameters τ0 , φ0 and ω0 are related to range, ele- A correlator stage as shown in Fig. 2 serves to despread the
vation and azimuth as follows signal y(t), incorporating the estimates τ̂0 and ω̂0 . The a priori
1 estimates are used to generate absolute time-varying references
τ0 = r0 + Cτ (4) τ̂0 (t), φ̂0 (t) for time-delay and phase for t ≥ 0. The time-delay
c
reference is generated as τ̂0 (t) = τ̂0 − ωω̂ c0 t, where the latter
φ0 = −kc r0 − Cφ (5)
term is provided by a numerically controlled oscillator (NCO)
ω0 = vkc cos(ε0 ) cos(ϑ0 ) + Cω (6) serving as an integrator. The time-delay reference is used to
drive the prompt replica s∗ (t − τ̂0 (t)). Likewise, the reference
where c denotes the speed of light and kc = ωc /c is the carrier for the carrier-phase is generated as φ̂0 (t) = ω̂0 t and drives the
wave number. The terms Cφ and Cτ represent instrumental
Doppler removal signal c∗ (t)  exp(− jφ̂0 (t)).
delays at receiver and satellite, or delays induced by atmospheric
Despreading of the waveform s(t) is accomplished with a set
effects, which are all considered (quasi-)constant and have to
of signal-matched correlators. We employ Q correlators con-
be determined by a later receiver stage for the actual ranging
trolled by delay hypotheses τ̂ (1) (t), . . . , τ̂ (Q ) (t) with
[1]. The Doppler shift Cω is due to satellite movement only.
The multipath signals’ time-delays, phases and Doppler shifts
for  > 0 are given by τ̂ (q ) (t) = τ̂0 (t) + (q )
, (10)

1 which cover at least the range of all possible LOS delays

τ = τ0 + ρ (7)
c (−TA , TA ) plus a margin TX ≥ 0 for excess delays; a reason-
φ = φ0 − kc ρ − π (8) able choice is −TA > (1) < (2) < . . . < (Q ) > TA + TX .
Selva [31] showed that the spacing between adjacent correlators
ω = vkc cos(ϑ ) + Cω (9) should not exceed 1/2Bs .

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7666 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

m
The discrete-time post-correlation signal obtained by the qth ejΔ ω  t at times tm = M Ti
correlator is given as the output of an integrate and dump filter
(q ) 1)  1 
nM

 x [n] ≈ P ejφ  I (tm )ejΔ ω  t m

1 nTi Ti
Z (q ) [n] = s∗ t − τ̂ (q ) (t) c∗ (t)y(t) dt, (11) m =(n −1)M +1
Ti (n −1)T i  tm    
ω̂0 ω
× s∗ t − τ̂ (q ) (0) + t s t − τ + t dt
t m −1 ωc ωc
with integration time Ti , from which N samples n = 1, . . . , N
shall be used for parameter estimation.  1 
nM

The signal (11) can be equivalently expressed as = P ejφ  I (tm )ejΔ ω  t m

M
m =(n −1)M +1
 mT  
∞
 1 Δω
∗
Z (q ) [n] =
(q )
x [n] + η (q ) [n] (12) × s t̃ − τ̂ (q )
(0) s t̃ − τ + t̃ dt̃
T (m −1)T ωc
=0
2) 
≈ P ejφ  rs Δτ − (q )
with the contribution of the th signal and noise, respectively,
1 
nM

   × I (tm )ejΔ ω  t m , (15)

1 nTi
ω M
(q ) ∗ m =(n −1)M +1
x [n] = s t − τ̂ (q )
(t) s t − τ + t
Ti (n −1)T i ωc
 where we used approximations 1) and 2) and the substitution
× I (t) P ejφ  ej(ω  −ω̂ 0 )t dt (13) t̃ = (1 + ωω̂ c0 )t. Defining η[n] = [η (1) [n], η (2) [n], . . . , η (Q ) [n]]T ,
√
 γ = P ejφ  and
1 nTi
(q )
η [n] = s∗ (t − τ̂ (q ) (t))c∗ (t) η(t) dt. (14) ⎡ ⎤
Ti (n −1)T i rs τ − (1)
⎢ ⎥
⎢ rs τ − (2) ⎥
In this representation of (11), we summarize over all  ≥ 0, indi- r(τ ) = ⎢
⎢ ..
⎥,
⎥ (16)
cating by I (t) ∈ {0, 1} whether  is included in L(t). However, ⎣ . ⎦
note that I0 (t) ≡ 1. rs τ − (Q )

the correlator output can be written in vector notation as

B. Approximations ∞

At this point, we define the residuals Δω  ω − ω̂0 and z[n] = γ r Δτ
=0
Δτ  τ − τ̂0 for  ≥ 0. Moreover, let the th path’s coherence
time be denoted by Δt . We make the following approximations 1 
nM

for a slowly time-varying channel, precise LOS Doppler removal × I (tm )ejΔ ω  t m + η[n]. (17)
M
and narrow bandwidth ranging signal, i.e., Δω Δt  2π  m =(n −1)M +1
ωc /Bs for all  ≥ 0.
1) For any t1 < t2 with |t2 − t1 | ≤ Δt , we may apply the III. WSSUS CHANNEL MODEL
piecewise approximations ejΔ ω  t ≈ ejΔ ω  t 2 and I (t) ≈
In the following, we present our assumptions on the randomly
I (t2 ) while t1 < t ≤ t2 . Typically, the LOS coherence
time-variant channel. Expressing (17) as a convolution
time is on the order of the total observation time, hence
Δt0 ≈ N Ti . The multipath coherence time for  > 0 can  
nM  ∞
 1
be much shorter; however, it can be bounded below by x[n]  z[n] = r(τ )h(τ, tm ) dτ,
Δt ≥ Ti /M ≈ T for some M ∈ N if receiver/scatterer η[n ]=0 M 0
m =(n −1)M +1
velocities are not excessive. (18)
2) Since s(t) has narrow bandwidth and its temporal resolu- we can identify the time-varying channel impulse response
tion is on the order of 1/Bs , we may apply the approxi- ∞
mation s(t + Δωωc  t) ≈ s(t) while t ≤ Δt [29]. 
h(τ, tm ) = γ I (tm )ejΔ ω  t m δ(τ − Δτ ). (19)
Note that with the above assumptions, the number of observa-
=0
tions N is essentially limited by the quality of the estimate ω̂0 .
Furthermore, we choose the integration time Ti = (1 + ωω̂ c0 )M T It is characterized by conditional (deterministic) LOS param-
equal to M code periods of duration T , multiplied by a factor to eters θ 0 = [Δτ0 , γ0 ]T and unconditional (random) multipath
correct for code expansion or compression due to the Doppler parameters λ (tm ) = [Δτ , Δω , γ , I (tm )]T for  > 0. Note
effect [8]. Now, we can approximate the integral in (13) by the that whenever we take expectations, θ 0 is in the condition-
integration over M subintervals with blockwise constant I (t), ing. Our central assumptions are that the channel mean and

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7667

exactly fulfilled if the receiver’s direction of travel is averaged

out [29]; for distinct directions of travel, it is an approximation
that is valid for chaotic propagation conditions in urban or in-
door environments, where the multipath angle of arrival ϑ can
be considered a uniform random variable [26]. These are the
relevant scenarios within the scope of this work.
5) Known multipath PDP structure: The PDP is defined as
E[|h(τ, tm )|2 ] = ψ(τ )(0) and indicates the expected multi-
path power density associated with delay τ at zero lag, i.e., for
Fig. 3. Power-delay channel statistics resulting from deterministic LOS and tm = tm . While the power term (0) must be considered un-
random multipath.
known, it is appropriate to make an assumption on the shape
of the PDP, hence on the delay profile ψ(τ ). Clearly, the PDP
autocorrelation have the following properties (cf. Fig. 3) has to be strictly zero for any τ < Δτ0 , since the LOS path is
  the shortest path. For urban terrestrial multipath channels, the
E h(τ, tm ) | θ 0 = const. (20) PDP is widely agreed to be a negative exponential, which is
  the immediate consequence of a uniform spatial distribution of
E h∗ (τa , tm )h(τb , tm ) | θ 0
scatterers determining the excess path lengths ρ in (7) and an
= δ(τa − τb ) |μ(τa )|2 + ψ(τa )(tm − tm ) (21) exponentially decaying path loss function [27]–[29]. Hansen
[30] showed analytically that the PDP inside a single room is
where the mean contribution is |μ(τa )|2 = P0 δ(τa − Δτ0 ), the indeed an exponential function, combining ray optical and sta-
sequence (tΔ ) expresses the (arbitrary) time-correlation of the tistical properties of wave propagation. Results in [18], [19] for
multipath channel, and ψ(τ  ∞) is an exponentially decaying mul- the satellite-to-earth channel also indicate an exponential PDP.
tipath delay profile with 0 ψ(τ ) dτ = 1. For convenience, we Therefore, we will always assume that the delay profile has
introduced tΔ  tm −m and τ  τa = τb . the form
Equations (20) and (21) express a series of assumptions on ⎧ τ −Δ τ
⎨ 1 e− τ ψ 0 if τ ≥ Δτ
which the estimator from Section V-A will rely. These are ex- ψ(τ ) = τψ
0
(22)
plained in the following from 1) to 5). While the time-correlation ⎩
0 otherwise,
(tΔ ) function is irrelevant for the estimator, the two additional
channel features 6) and 7) suggest one possible model, which where τψ is the unknown expected excess delay.
will be used in parts of our simulations. Note that we do not impose a general restriction on the time-
1) Constant mean :Since Δω0 ≈ 0 and I0 (tm ) ≡ 1, the LOS correlation function (tΔ ), but the following models for Δω
response γ0 I0 (tm )ejΔ ω 0 t m δ(τ − Δτ0 ) is quasi-constant. On the and I (tm ) will be used for simulations and are a plausible
other side, multipath contributions are modeled as having mean model for urban multipath channels.
zero, since the echo phases φ correspond to many thousands 6) Jakes distribution of multipath Doppler: Consider the ge-
of wavelengths. Even when conditioned on φ0 , they can be ometry shown in Fig. 1. Modeling chaotic receiver environment,
considered uniformly distributed in (−π, π) if the corresponding we assume iid uniform angles of arrival ϑ in (−π, π) following
scatterers are not located in the near field of the receive antenna [26]. Upon precise LOS Doppler removal (Δω0  2π/N Ti ),
[27]. each multipath’s relative Doppler frequency Δω is found to be
2) Uncorrelated scattering (US): A US channel is a direct iid according to a shifted Jakes distribution
consequence of assuming a spatially uncorrelated field of point 1
scatterers [29]. Real scattering objects such as buildings or trees pΔ ω  (ω) =  , |ω − λωv | ≤ ωv
π ωv2 − (ω − λωv )2
get more and more indistinguishable from point scatterers as the (23)
spatial resolution c/Bs of the ranging signal decreases. There- with λ = cos(ϑ0 ) cos(ε0 ) and ωv = vc ωc .
fore, while US yields an excellent approximation for the small 7) Multipath diversity (Ergodicity): The channel is assumed
bandwidths that are applied for mobile communications, it is to provide some kind of multipath diversity, so as to allow
still a reasonable assumption even for the example of a 10.23 estimation of multipath statistics after collecting echoes over
MHz GNSS receiver, for which c/Bs = 29.3 m. In terms of time. However, instead of the simple intuitive choice of very
our model, US means that λ1 (tm ), λ2 (tm ), . . . are uncorrelated many simultaneously active echoes |L(tm )| → ∞ with less and
conditioned on θ 0 . less power P → 0, we propose a more relaxed model with a
3) (Short-term) wide-sense stationarity: A WSS model can finite moderate number of simultaneously active paths chosen
only capture short-term channel variations, which implies the from a large set of possible paths
assumption of constant Doppler frequencies and WSS processes
L
I (tm ). Propagation conditions are not allowed to change sub-   m ax
 
stantially within the observation time N Ti . E |L(tm )| = lim E I (tm ) = L < ∞ (24)
L m a x →∞
=0
4) Temporal/delay separability: Assuming the separability
ψ(τ )(tΔ ) means that we consider time-delay Δτ and fre- and P = (0)/L. As the receiver is moving, echoes disappear
quency shift Δω of an echo as uncorrelated. This property is and new ones appear randomly following a birth-death model.

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7668 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

Fig. 4. Binary Markov chain model for appearing and disappearing multipath.

Fig. 6. Multipath time-correlation sequence for various receiver velocities v;

T = 1 ms, d = 1 m, λ = 0.

receiver velocity is illustrated in Fig. 6. Clearly, for v = 0,

|(tΔ )| → 0 as tΔ → ∞. Further, we define the power mix-
ing coefficients for averaging time Ti = (1 + ωω̂ c0 )M T as
Fig. 5. Simulation example for randomly evolving number of active multi-
paths over time; v = 5 m/s, d = 1 m, L = 15, L m a x = 5 × 103 .
1  
M M
P(k ) =  tk M +m −m
For finite Lm ax , this model can be used for Monte-Carlo simu- M2
m =1 m =1
lations as follows:

M −1  
For each multipath  = 1, . . . , Lm ax , the sequence I (tm ) ∈ 1 |Δ|
{0, 1} is a discrete-time first-order binary Markov chain with = 1−  tk M +Δ (28)
M M
Δ =1−M
birth probability α and death probability β, as illustrated in
Fig. 4. These processes are stationary [32]; additionally, we
show in Appendix B that for k = 0, . . . , N − 1.
  α Second, for v = 0, and if L, Lm ax are sufficiently large, the
E I (tm ) = (25) channel h(tm ) becomes a Gaussian random process: if L is
α+β
sufficiently large, then for any fixed tm , h(tm ) can be consid-
  αβ(1 − α − β)|m −m | + α2 ered a Gaussian random variable by the standard central limit
E I (tm )I (tm ) = . (26)
(α + β)2 theorem [32]; further, as Lm ax → ∞ and since the contribution
of each multipath is a stationary random process, the random
The death probability β = 1 − e−v |t m −m |/d depends on receiver process h(tm ) becomes ergodic and converges in distribution
speed and the average echo life span d (meaning the spatial dis- to a Gaussian process according to the central limit theorem for
tance across which a path is actually present); the birth proba- random processes [35]. In particular, the joint distribution of
bility can then be chosen as α = L mβa xL−L . For urban and indoor any collection of data z[1], . . . , z[N ] can be considered jointly
scenarios, a typical value for the life span is d = 1 m, while L Gaussian.
ranges mostly from 5 to 50 [18]–[21]. Similar birth-death pro-
cesses have been used to characterize the random occurrence
of multipaths [33], [34]. An exemplary realization of |L(tm )| is IV. RECEIVED SIGNAL STATISTICS
shown in Fig. 5.
We can make two observations from the properties 6) In the following, we derive the first- and second-order
and 7). First, time-variance and diversity of the channel are statistics of the data sequence z[n] across temporal and delay
a result from the receiver moving with speed v. For the dimension. Noise η[n] and useful signal x[n] are assumed
above model, the time-correlation function can be shown independent. While Section IV-A deals with the noise covari-
to be ance, Section IV-B provides a derivation of the useful signal’s
moments using the channel model from Section III.
L
αβ(1 − α − β)|Δ | + α2 jλω v t Δ m ax

(tΔ ) = e J0 ωv tΔ P ,
(α + β)2
=1
A. Covariance of Correlator Output Noise
(27)
where J0 ( · ) denotes the 0th-order Bessel function of the first When determining the noise covariance at the correlator out-
kind. The dependency of the time-correlation function on the puts, it is important to account for the correlations introduced

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7669

by the correlator bank. These are given by

 
E ηq∗ [n]ηr [n + k]
 (n +k )T i  n T i
1  
= 2 c(t1 )c∗ (t2 ) E n∗ (t1 )n(t2 )
Ti (n +k −1)T i (n −1)T i

× s t1 − τ̂ (q ) (t1 ) s∗ t2 − τ̂ (r ) (t2 ) dt1 dt2

  nTi
N0 (n +k )T i
= 2 c(t1 )c∗ (t2 )δ(t2 − t1 )
Ti (n +k −1)T i (n −1)T i

× s t1 − τ̂ (q ) (t1 ) dt1 s∗ t2 − τ̂ (r ) (t2 ) dt2


N0 (n +k )T i
= 2 Πn (t/Ti ) s t − τ̂ (q ) (t) s∗ t − τ̂ (r ) (t) dt
Ti (n +k −1)T i
 nM T (0) (1)
Fig. 7. Multipath power mixing coefficients P  (solid) and |P  | (dotted)
N0
= δ[k] s(t̃ − τ̂ (q ) (0))s∗ (t̃ − τ̂ (r ) (0)) dt̃ for d = 1 m as a function of the coherent averaging distance vT i , with various
(M T )2 (n −1)M T values of geometry coefficient λ.

N0
= δ[k] rs ( q − r ), (29) with the multipath delay correlation matrix
MT
 ∞ τ −Δ τ 0
where we used |c(t)|2 = 1 and Rψ =
1 −
r(τ )r H (τ )e τ ψ dτ. (35)
 τψ Δ τ 0
1 if n − 1 ≤ t
Ti <n
Πn (t/Ti ) = (30)
0 otherwise Note that Rψ does not only depend on τψ , but also on Δτ0 .
 Therefore, the received signal’s covariance is sensitive with re-
1 if k = 0 spect to the LOS time-delay, which reflects the fact that echoes
δ[k] = (31)
0 otherwise. do actually carry valuable information about the LOS. Even
Thus, we have noise that is uncorrelated over time but correlated with a (severely) attenuated LOS signal (γ0 → 0), it may still
among the correlator outputs with be possible to obtain an estimate for Δτ0 based on multipath
  signals.
E η[n + k]η H [n] (k )
The multipath power mixing coefficients P are a measure
N0  
of multipath correlation between samples x[n] and x[n + k].
= δ[k] r( 1 ) r( 2 ) . . . r( Q )  δ[k] Ση . (32) (0)
MT As can be observed in Fig. 7, we have always 0 < P ≤ (0).
The upper bound is obtained for static scenarios (vTi → 0)
B. Mean and Covariance of Useful Signal
and the lower bound is obtained for large integration distances
Using (18)-(22) and (28), the moments of x[n] are (0)
(vTi → ∞). In general, P decreases with vTi , which is the
 
μx|θ0  E x[n] | θ 0 = γ0 r(Δτ0 ) (33) well-known and often desired effect of ”averaging out” the mul-
tipath components. This is because the Doppler effect decorre-
and, for n = 1, . . . , N , k = 1 − n, . . . , N − n, lates multipath from the LOS as the receiver is traveling, if
(k )  
Σx|θ0  E x[n + k]xH [n] | θ 0 − μx|θ0 μH the receiver keeps track of ω0 . The decorrelation is stronger if
x|θ 0
|λ| → 0. As |λ| → 1 (e.g. in low elevation scenarios), more aver-

(n +k )M

nM  ∞ ∞ (0)
aging is necessary before P attains the asymptotic decay rate
1
= dτa dτb
of (vTi )−1 . The ratio P )/(0) can be interpreted as the degree
(0
M2 0 0
m =(n +k −1)M +1 m =(n −1)M +1
of multipath suppresion resulting from averaging. The integra-
 
× r(τb )r H (τa )E h∗ (τa , tm )h(τb , tm ) | θ 0 tion time Ti should be chosen depending on the specific scenario
(1) (0)
(λ, v) according to Fig. 7, such as to satisfy |P |  P . In
− P0 r(Δτ0 )r H (Δτ0 ) that case, successive samples z[n], z[n + 1] can be considered
(1)

(n +k )M

nM uncorrelated. For instance, the figure shows that |P | is 10 dB
1 (0)
=  tm −m below P already for moderate integration distances vTi (for
M2
m =(n +k −1)M +1 m =(n −1)M +1 instance, 15 ms × 30ms = 45cm). For k > 1, the correlations get
 ∞ ∞ even smaller.
× r(τb )r H (τa )ψ(τa )δ(τb − τa ) dτa dτb
0 0
V. PARAMETER ESTIMATION
 ∞
1  
M M
= 2  tk M +m −m r(τ )r H (τ )ψ(τ ) dτ From the N data samples z[1], . . . , z[N ] given by (17), the
M 0 estimation problem is to infer at least the LOS time-delay Δτ0 .
m =1 m =1
In Section V-A, we present an iterative scheme to produce
= P(k ) Rψ (34) the U-MLE that is optimal for a Gaussian WSSUS channel.
Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7670 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

 T 
Specifically, the U-MLE estimates θ 0 jointly with an addi-
(0) τ̂ψ = arg min Φ Δτ̂0 , γ̂0 , τψ , P̂(0) (43)
tional parameter θ X  [τψ , P ]T describing the multipath’s τ ψ ∈Tψ
power-delay statistics. On the other hand, the coherent averag-  T 
(0)
ing estimator (CAE) presented in Section V-B tries to suppress Δτ̂0 = arg min Φ Δτ0 , γ̂0 , τ̂ψ , P̂ (44)
Δ τ 0 ∈T0
multipath by taking the coherent average over all N blocks
before estimating θ 0 alone. until the sequence of estimates converges. In each of these steps,
parameters that are not subject to the optimization are kept fixed
A. Unconditional Maximum Likelihood Estimator (U-MLE) at the last available estimate.
We recommend line searches over finite intervals P , Tψ , T0
The proposed estimator assumes N independent and iden- for solving (42)-(44), which require the calculation of an inverse
tically distributed (iid) observations z[n] ∼ CN (μz , Σz ) fol- and a determinant of the Q × Q-matrix Σ for each cost function
lowing a Gaussian distribution with mean and covariance evaluation, Q being the number of correlators. The minimiza-
tion problem (41) should be solved analytically as follows. A
μz = μx|θ0 (36)
necessary condition for finding a minimum is that the Wirtinger
(0)
Σz = Σx|θ0 + Ση . (37) derivative [37] vanishes
∂Φ(θ)
This model is fully characterized by the unknown parameter = 0, (45)
∂γ0∗
⎡ ⎤
Δτ0 from which we obtain by a straightforward calculation
 
⎢ γ0 ⎥
θ=⎢ ⎥ = θ0 . (38) r H Δτ0 Σ−1
z z
⎣ τψ ⎦ θX γ̂0 = , (46)
(0) r H Δτ0 Σ−1
z r Δτ0
P
N
where z  N1 n =1 z[n]. The condition (45) is also sufficient
The U-MLE is given as the minimizer
for a minimum since Σz is positive definite; hence, for all
θ̂ = arg min Φ(θ) (39) parameter values
θ
∂ 2 Φ(θ)
of the nonlinear least-squares cost function = N r H (Δτ0 )Σ−1
z r(Δτ0 ) > 0. (47)
∂γ0 ∂γ0 ∗
 
N  Convergence of the above sequence of estimates to the U-MLE
Σ−1
H
Φ(θ) = tr z z[n] − μz z[n] − μz is not generally guaranteed, since it may converge to any local
n =1 minima of Φ(θ). Therefore, the initial guess θ̂ should already
+ N log det Σz (40) be located near the global minimum of Φ(θ). This issue is
discussed in Section VI-A.
Besides presuming a Gaussian model, the U-MLE neglects
(k )
the nonzero cross-covariances Σx|θ0 of the observations B. Coherent Averaging Estimator (CAE)
z[n], z[n + k] for k = 0. While no optimality claims can be The CAE is comparable with estimators used by high-
made for such an estimator in the strict sense, it is actually ML- sensitivity GNSS receivers with longer and longer coherent
(k ) (0)
optimal if |P |  P for k = 0 and if z[1], . . . , z[N ] can averaging times. It is ML-optimal for the multipath-free case
be considered jointly Gaussian. This was argued to be the case (0) (0)
(P → 0). Evaluating (46) for P = 0 leads to
for at least moderate averaging time Ti and moderate number
of simultaneously active echoes L in Sections III and IV. If r H (Δτ0 )Σ−1
η z
γ0 = , (48)
these assumptions can be made, the U-MLE can be considered r(Δτ0 )Σ−1
η r(Δτ0 )
asymptotically unbiased and asymptotically attaining the CRLB
which yields the estimation for the LOS time-delay
as N → ∞ [36]. Otherwise, the U-MLE may not be optimal but  T 
still have reasonable performance and low model complexity. In Δτ 0 = arg min Φ Δτ0 , γ 0 , 0, 0
any case, the CRLB for the iid Gaussian model provides an ul- Δτ0
timate performance benchmark for the proposed estimator, and  H 2
r (Δτ0 )Σ−1
η z

is derived in Appendix A. = arg max . (49)
To solve (39) iteratively with low complexity, we propose Δτ0 r H (Δτ0 )Σ−1
η r(Δτ0 )
to choose a starting point θ̂ and then to perform a grouped Obviously, if there is no multipath reception, the sample mean
coordinate descent by repeating the following steps z is a sufficient statistic for estimation of θ 0 and coherent aver-
 T  aging over all N is the optimal strategy.
γ̂0 = arg min Φ Δτ̂0 , γ0 , τ̂ψ , P̂(0) (41)
γ 0 ∈C
 VI. SIMULATIONS
T 
P̂(0) = arg min Φ Δτ̂0 , γ̂0 , τ̂ψ , P(0) (42) For the two estimators proposed in Section V, we deter-
( 0)
P  ∈P mine the root mean squared error (RMSE) with respect to Δτ0

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7671

TABLE I
GENERAL SIMULATION PARAMETERS

10 log 10 (C /N 0 ) v ε0 ϑ0 Ti

◦ ◦
43.2 dB-Hz 15 m/s 30 45 30 ms

Fig. 9. Scenario 1: RMSE as a function of the number of observations N .

TABLE II
PARAMETERS FOR SCENARIO 1

10 log 10 (P 0 /C ) 10 log 10 (P  /C ) τψ L
Fig. 8. Autocorrelation functions of ranging waveforms: 1.023 MHz binary
phase-shift keying and 7×1.023 MHz time-multiplexed binary offset carrier. −30.0 dB −30.0 dB 150 ns 18

numerically for various scenarios. To generate realizations yield the global minimum of Φ(θ) in some cases even in the
h(τ, tm ), we will use the WSSUS model discussed in asymptotic regime, i.e., for large N . This is because the auto-
Section IV-B as well as the landmobile multipath channel model correlation function rs (τ ) is generally multimodal, i.e., it can
(LMMCM) in urban environments [22] standardized by the In- have more than one peak. A prominent example is the L1C Pilot,
ternational Telecommunication Union (ITU). General parame- whose TMBOC(7) autocorrelation function has clear sidelobes
ters valid for all simulations are provided in Table I. The given at ±42 m. Furthermore, even for unimodal rs (τ ), a strong and
value for the nominal carrier-to-noise density C/N0 refers to persistent echo  with Δω ≈ Δω0 may lead to an extra peak of
unblocked conditions. the cost functions in (43) and (44). However, this is acceptable
For the ranging waveform s(t), we use the legacy GPS L1 C/A as the U-MLE is naturally mismatched and not intended for sce-
signal with binary phase-shift keying modulation as received by narios with quasi-static specular paths, which is also reflected
a mass-market receiver (Bs = 1.023 MHz), as well as the GPS in the simulation results for small L in Section VI-B. Locking
L1C (“New Civil”) Pilot signal with time-multiplexed binary on sidelobes introduced by multimodal rs (τ ) is avoided by let-
offset carrier modulation [38] as received by a professional re- ting TA ≤ 1/Bs . For the search spaces in (42)-(44), we suggest
ceiver (Bs = 7 × 1.023 MHz); we refer to them by BPSK(1) to use T0 = (−TA , TA ) and the (physically relevant) intervals
and TMBOC(7), respectively. Both are transmitted at the L1 P = (C × 10−3 , C), Tψ = (0, 500 ns), where C is the LOS
frequency ωc /(2π) = 1575.42 MHz. The respective autocorre- power in nominal unblocked conditions.
lation functions are shown in Fig. 8. While TMBOC(7) generally
leads to better performance at the cost of higher pre-correlation B. WSSUS Channel Model
sampling rate, it should be kept in mind that these kinds of sig-
nals also need a smaller acquisition interval to detect the narrow Using the WSSUS model from Section IV-B, we determine
autocorrelation peak and to prevent the fine estimation (track- the RMSE as a function of the characteristic parameters N ,
ing) from locking onto the autocorrelation sidelobes. We will P0 , τψ and L. For all simulations, we use a simulator size of
always assume an acquisition interval Δτ0 ∈ [−TA , +TA ] with Lm ax = 5 × 103 , an echo life span d = 1.0 m and 103 Monte-
TA = 1/Bs . This corresponds to 293 meters for BPSK(1) and Carlo runs. While 10 log10 (C/N0 ) = 43.2 dB-Hz in unblocked
to 42 meters for TMBOC(7). conditions, the actual LOS and multipath power used for the
respective shadowing/scattering scenarios will be given as the
A. Initial Guess, Search Space and Local Minima of Φ(θ) ratio P0 /C or P /C, respectively, in dB.
Note from Scenario 1 (Fig. 9, Table II) that relatively large
For the initital guess, we suggest θ̂ = [0, 0, 0, 0]T , i.e., to numbers of observations N are necessary to obtain acceptable
simply rely on the acquistion estimate τ̂0 , and then solving for range errors for all scenarios of interest, since LOS and multipath
(0)
γ0 , P , τψ , τ0 iteratively. Iterating the steps (42)-(44) may not power levels are generally far below the nominal carrier-to-noise

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7672 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

Fig. 10. Scenario 2: RMSE as a function of LOS obstruction P 0 /C .

TABLE III
PARAMETERS FOR SCENARIO 2 Fig. 11. Scenario 3: RMSE as a function of the expected excess delay τ ψ .

10 log 10 (P  /C ) τψ N L

−20.0 dB 150 ns 150 18

TABLE IV
PARAMETERS FOR SCENARIO 3

10 log 10 (P 0 /C ) 10 log 10 (P  /C ) N L

−36.0 dB −24.0 dB 250 18

density. For this scenario, we chose LOS power P0 and multipath

powers P1 , P2 , . . . to have the same value. In this case, using the
proposed U-MLE provides moderate gain when compared with
CAE. However, the benefit gets larger when LOS shadowing
gets stronger and/or echoes become more dominant. Fig. 12. Scenario 4: RMSE as a function of the expected number of simulta-
neously active paths L.
We consider the case of a LOS signal with different degrees
of obstruction in Scenario 2 (Fig. 10, Table III). The lower the TABLE V
LOS power P0 is, the more must estimators rely on echoes to PARAMETERS FOR SCENARIO 4
determine the LOS delay. Interestingly, for the U-MLE, using
TMBOC(7) instead of the simple BPSK(1) does not bring much 10 log 10 (P 0 /C ) 10 log 10 (L P  /C ) τψ N
benefit as P0 → 0. TMBOC was designed for good multipath
−36.0 dB −11.4 dB 250 ns 250
rejection in the presence of a strong LOS signal. However, for
urban scenarios with frequent shadowing, using a simple signal
such as BPSK in combination with multipath-exploiting esti-
mators may in fact outperform advanced wideband waveforms can be observed that the CRLB increases unexpectedly, whereas
in combination with long coherent averaging − presumed, of the RMSE decreases. We suppose that, as the echoes degener-
course, that there is scattering: in this scenario, there are (on ate to a LOS signal, the Fisher information measure becomes
average) 18 echoes with a power of −20 dB with respect to singular due to parameter identifiability problems [39].
nominal C/N0 . So far, we have assumed an average of L = 18 simultaneously
In Scenario 3 (Fig. 11, Table IV), RMSE is given as a func- active echoes, which would correspond to diffuse scattering.
tion of the average excess delay τψ . Besides the spatial distri- This number seems to be sufficiently large to assume a Gaussian
bution of scatterers, satellite elevation has been shown to be the channel, since the Gaussian CRLB has been attained. When the
main factor determining excess delay [18], [21]. The U-MLE same total power LP is concentrated in less active paths L
performs stable with respect to variations of τψ . For τψ → 0, it in Scenario 4 (Fig. 12, Table V), the CRLB is not attained.

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7673

Fig. 13. RMSE for landmobile multipath channel model “urban/car”.

Fig. 14. Percentage of estimates exceeding the acquisition interval for land-
mobile multipath channel model “urban/car”.

Especially for one or two specular echoes, direct estimation

of multipath parameters [9]–[15] is the better-suited solution;
nevertheless, the U-MLE performance is robust in this case of
model mismatch.

C. ITU: Landmobile Multipath Channel Model (LMMCM)

The LMMCM was derived from measurement data recorded
in Munich [18], [21] and was standardized by ITU-R [22]. Signal
propagation in urban and suburban areas was measured for car
and pedestrian applications with a 100 MHz channel sounding
signal near the L1 frequency. The model output is a complex
time-variant channel impulse response of the form


h(τ, tm ) = γ (tm )δ(τ − Δτ (tm )), (50)
∈L(t)
Fig. 15. Two exemplary realizations of urban LMMCM path amplitudes.
(a) Rich scattering. (b) Dead zone.
and allows for correlated scattering as well as delay drifts. As
user input to the model, we use the parameters given in Table I
to describe movement on a straight line, and select the scenario VII. CONCLUSION
“urban/car”. The number of observations is N = 100, with in- We introduced a methodology to improve accuracy and ro-
tegration time Ti = 30 ms per observation. Thus, for each of bustness of LOS synchronization and ranging for GNSS in mul-
the 2 × 103 Monte-Carlo runs, three seconds of the LMMCM tipath environments. Our solution is based on a low-complexity
output are used. The resulting RMSE is shown in Fig. 13. parametric estimation of a WSSUS echo distribution. In post-
In general, the U-MLE always outperforms the CAE, since correlation domain, this distribution can be characterized with
(0)
the channel outputs often resemble the “rich scattering” example two probabilistic parameters, namely, the expected power P
in Fig. 15(a) with many echoes. However, the performance of and excess delay τψ . In various dynamic scenarios with specu-
both estimators is significantly worse than with our WSSUS lar/dense scattering and LOS/NLOS conditions, we showed that
model in Section VI-B. This may be because the estimators the proposed U-MLE outperforms the high-sensitivity CAE,
often have to cope with channel realizations such as shown in since it is designed to exploit echoes with known PDP structure
Fig. 15(b), where there are only few weak echoes and no direct rather than to suppress them. Simulations with an ITU-approved
path for durations on the order of seconds. In such “dead zones”, non-WSSUS urban channel model proved that the method is ro-
estimates are severly affected by noise and even tend to leave bust to some degree of model mismatch. In the model-matched
the acquisition interval (−1/Bs , 1/Bs ). Whenever an estimate scenario, the U-MLE attains the CRLB if many uncorrelated
Δτ̂0 is equal to ±1/Bs , it is excluded from the statistic and multipath signals are received, which can be collected over time
considered as outlier (Fig. 14); this was observed in the case of in case of a mobile receiver; for static receivers or for very few
TMBOC. paths that lack diversity, this cannot be expected. The benefit of

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7674 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

using the proposed estimator instead of long coherent averaging APPENDIX B

is greatest when the echoes become more dominant and when AUTOCORRELATION OF MARKOVIAN MULTIPATH MODEL
the LOS is obstructed. While state-of-the-art GNSS receivers Consider a binary Markov sequence Xm ∈ {0, 1} with
often detect and discard satellites with NLOS propagation, the
m ∈ Z. To find the desired moments E[Xm ] and E[Xm +k Xm ],
proposed methodology provides a perspective on how these
we first calculate the joint probability mass function (PMF) of
satellites could still be useful for navigation. Xm +k and Xm , which is denoted as pX m + k ,X m (x, y), in steady
state. Let the transition probabilities as shown in Fig. 4 be given
APPENDIX A as follows:
CRAMÉR-RAO LOWER BOUND (CRLB)    
pX m + 1 |X m (0|0) pX m + 1 |X m (0|1) 1−α β
For a (real-valued) parameter θ, the CRLB is a lower bound P  =
pX m + 1 |X m (1|0) pX m + 1 |X m (1|1) α 1−β
for the RMSE of any unbiased estimator θ̂. The covariance (62)
matrix of such an estimator satisfies E[θ̂ θ̂ T ] − F −1 (θ)  0, for all m ∈ Z. First, we recall from [32] that the steady-state
where  0 means that the matrix is positive semidefinite [36]. PMF pX m (x) can be obtained from the condition
The Fisher information matrix F (θ), under the assumption of
N iid samples z[n] ∼ CN (μz , Σz ), is given as p(∞) = P p(∞) , (63)
 
∂Σz −1 ∂Σz where p(∞)  [pX m (0) pX m (1)]T , to yield
F (θ) = N tr Σ−1
z Σ ⎡ ⎤
ij ∂θ i z ∂θ j β
  ⎢α+β ⎥
∂μH ∂μz
+ N 2 Re z
Σ−1 (51) p(∞) = ⎢ ⎣ α ⎦
⎥ (64)
∂θ i z ∂θ j
α+β
and is evaluated at the true value of θ.
(0) and E[Xm ] = α/(α + β). Furthermore, we recall Bayes’ rule
For θ = [Δτ0 , |γ0 |, ∠(γ0 ), τψ , P ]T , the CRLB can be cal-
culated from the following partial derivatives pX m + k ,X m (x, y) = pX m + k |X m (y|x) pX m (x) (65)

∂μz ∂r(Δτ0 ) and the Markov property for positive k

= γ0 (52)
∂Δτ0 ∂Δτ0
pX m + k |X m + k −1 ,X m + k −2 ,... (x|y, z, . . .) = pX m + k |X m + k −1 (x|y).
∂μz (66)
= jγ0 r(Δτ0 ) (53)
∂∠(γ0 ) Marginalizing over all intermediate states, we can calculate
∂μz for k ≥ 0
= ej∠(γ 0 ) r(Δτ0 ) (54)
∂|γ0 | pX m + k ,X m (x, y)
∂μz
=0 (55) 
1
∂τψ = pX m + k ,X m + k −1 |X m (x, i1 |y) pX m (y)
∂μz i 1 =0
=0 (56)
∂P
(0)

1
= pX m + k |X m + k −1 (x|i1 ) pX m + k −1 |X m (i1 |y) pX m (y)
(0)
∂Σz P i 1 =0
= Rψ − r(Δτ0 )r H (Δτ0 ) (57)
∂Δτ0 τψ 
1 
1 
1

∂Σz = ... P x,i 1 P i 1 ,i 2 . . . P i k ,y p(∞)

y
=0 (58) i 1 =0 i 2 =0 i k =0
∂∠(γ0 )
∂Σz = eTx P k ey p(∞)
y , k ≥ 0, (67)
=0 (59)
∂|γ0 | (∞)
where P i,j refers to the (i + 1, j + 1)-th entry of P and py
(0)  ∞
∂Σz P −
τ −Δ τ 0
τ − Δτ0 − τψ refers to the (y + 1)-th entry of p(∞) . Moreover, e0 and e1
= r(τ )r H (τ )e τψ
dτ
∂τψ τψ Δτ0 τψ2 denote the unit vectors [1, 0]T and [0, 1]T . Using the fact that a
(60) Markov chain upon time axis inversion is still a Markov chain,
a similar result can be obtained for k < 0
∂Σz
= Rψ , (61) pX m + k ,X m (x, y) = eTx P −k ey p(∞)
∂P
(0) y , k <0 (68)

and merging the two results yields

given the autocorrelation function of the ranging signal r(τ ) as
sampled by the correlator bank. pX m + k ,X m (x, y) = eTx P |k | ey p(∞)
y , k ∈ Z. (69)

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
ENNEKING AND ANTREICH: EXPLOITING WSSUS MULTIPATH FOR GNSS RANGING 7675

The k-th power of P can be calculated by using the eigenvalue [11] N. Blanco-Delgado and F. Nunes, “Multipath estimation in multicorrelator
decomposition GNSS receivers using the maximum likelihood principle,” IEEE Trans.
Aerosp. Electron. Syst., vol. 48, no. 4, pp. 3222–3233, Oct. 2012.
P = V DV −1 , (70) [12] F. Antreich, J. Nossek, G. Seco-Granados, and A. Swindlehurst, “The ex-
tended invariance principle for signal parameter estimation in an unknown
spatial field,” IEEE Trans. Signal Process., vol. 59, no. 7, pp. 3213–3225,
where D is diagonal, its entries being the eigenvalues Jul. 2011.
[13] M. Sahmoudi and M. Amin, “Fast iterative maximum-likelihood algo-
d1 = 1 − α − β (71) rithm (FIMLA) for multipath mitigation in the next generation of GNSS
receivers,” IEEE Trans. Wireless Commun., vol. 7, no. 11, pp. 4362–4374,
d2 = 1, (72) Nov. 2008.
[14] P. Closas, C. Fernandez-Prades, and J. Fernandez-Rubio, “A Bayesian
while the matrices of left and right eigenvectors are given as approach to multipath mitigation in GNSS receivers,” IEEE J. Sel. Topics
Signal Process., vol. 3, no. 4, pp. 695–706, Aug. 2009.
⎡ 1 β ⎤ [15] M. Lentmaier, B. Krach, and P. Robertson, “Bayesian time delay estima-
√  tion of GNSS signals in dynamic multipath environments,” Int. J. Navig.
⎢ 2 α2 + β 2 ⎥
V =⎢ ⎥ Observ., vol. 2008, 2008, Art. no. 372651.
⎣ 1 ⎦ (73) [16] C. Gentner, R. Pöhlmann, T. Jost, and A. Dammann, “Multipath assisted
α
−√  positioning using a single antenna with angle of arrival estimations,” in
2 α2 + β 2 Proc. 27th Int. Tech. Meet. Satellite Division Inst. Navig., Tampa, FL,
⎡ √ √ ⎤ USA, Sep. 2014, pp. 1690–1697.
α 2 −β 2 [17] F. Wendler, F. Antreich, J. A. Nossek, and A. L. Swindlehurst, “Dual-
⎢ α+β α+β ⎥ polarization time delay estimation for multipath mitigation,” in Proc. 19th
⎢ ⎥
V −1 = ⎢   ⎥. (74) Int. ITG Workshop Smart Antennas, Mar. 2015, pp. 1–6.
⎣ α2 + β 2 α2 + β 2 ⎦ [18] A. Lehner and A. Steingaß, “Spatial dynamic wideband modeling of the
MIMO satellite-to-earth channel,” Int. J. Antennas Propag., vol. 2014,
α+β α+β 2014, Art. no. 102754.
[19] F. Fontan, M. Vazquez-Castro, C. Cabado, J. Garcia, and E. Kubista,
Finally, we obtain the desired autocorrelation from the PMF “Statistical modeling of the LMS channel,” IEEE Trans. Veh. Technol.,
vol. 50, no. 6, pp. 1549–1567, Nov. 2001.
  
1 
1 [20] T. Jost, W. Wang, U.-C. Fiebig, and F. Perez-Fontan, “A wideband satellite-
E Xm +k Xm = xy pX m + k ,X m (x, y) to-indoor channel model for navigation applications,” IEEE Trans. Anten-
x=0 y =0 nas Propag., vol. 62, no. 10, pp. 5307–5320, Oct. 2014.
[21] A. Lehner, “Multipath channel modelling for satellite navigation systems,”
= pX m + k ,X m (1, 1) Ph.D. dissertation, Universität Erlangen-Nürnberg, Erlangen, Germany,
2007.
[22] ITU-R, “Model parameters for an urban environment for the physical-
= eT1 V D |k | V −1 e1 p1
(∞)
statistical wideband LMSS model in recommendation ITU-R P.681-6,” In-
ternational Telecommunication Union, Geneva, Switzerland, Rep. P.2145,
αβ(1 − α − β)|k | + α2 2009.
= . (75) [23] J. Rissanen, Information and Complexity in Statistical Modeling. Berlin,
(α + β)2 Germany: Springer, 2007.
[24] C. Ribeiro, A. Richter, and V. Koivunen, “Joint angular- and delay-
REFERENCES domain MIMO propagation parameter estimation using approximate ML
method,” IEEE Trans. Signal Process., vol. 55, no. 10, pp. 4775–4790,
[1] E. Kaplan and C. Hegarty, Understanding GPS: Principles and Applica- Oct. 2007.
tions (Ser. Artech House Mobile Commun.). Norwood, MA, USA: Artech [25] O. Bialer, D. Raphaeli, and A. Weiss, “Efficient time of arrival estimation
House, 2005. algorithm achieving maximum likelihood performance in dense multi-
[2] R. van Nee, “Spread-spectrum code and carrier synchronization errors path,” IEEE Trans. Signal Process., vol. 60, no. 3, pp. 1241–1252, Mar.
caused by multipath and interference,” IEEE Trans. Aerosp. Electron. 2012.
Syst., vol. 29, no. 4, pp. 1359–1365, Oct. 1993. [26] W. Jakes, Microwave Mobile Communications (ser. IEEE Press Classic
[3] F. Scire-Scappuzzo and S. N. Makarov, “A low-multipath wideband GPS Reissue). Piscataway, NJ, USA: IEEE Press, 1974.
antenna with cutoff or non-cutoff corrugated ground plane,” IEEE Trans. [27] J. Parsons, The Mobile Radio Propagation Channel. Hoboken, NJ, USA:
Antennas Propag., vol. 57, no. 1, pp. 33–46, Jan. 2009. Wiley, 2000.
[4] M. Amin and W. Sun, “A novel interference suppression scheme for [28] C. of the European Communities, Digital Land Mobile Radio Communica-
global navigation satellite systems using antenna array,” IEEE J. Sel. tions. COST 207 (Ser. Information Technologies and Sciences) Brussels,
Areas Commun., vol. 23, no. 5, pp. 999–1012, May 2005. Belgium: Eur. Commission, 1989.
[5] S. Daneshmand, A. Broumandan, J. Nielsen, and G. Lachapelle, “Interfer- [29] J. Sadowsky and V. Kafedziski, “On the correlation and scattering func-
ence and multipath mitigation utilising a two-stage beamformer for global tions of the WSSUS channel for mobile communications,” IEEE Trans.
navigation satellite systems applications,” IET Radar, Sonar Navig, vol. 7, Veh. Technol., vol. 47, no. 1, pp. 270–282, Feb. 1998.
no. 1, pp. 55–66, Jan. 2013. [30] J. Hansen, “An analytical calculation of power delay profile and delay
[6] A. J. van Dierendonck, P. Fenton, and T. Ford, “Theory and performance spread with experimental verification,” IEEE Commun. Lett., vol. 7, no. 6,
of narrow correlator spacing in a GPS receiver,” J. Inst. Navig., vol. 39, pp. 257–259, Jun. 2003.
no. 3, pp. 265–284, 1992. [31] J. Selva, “Complexity reduction in the parametric estimation of superim-
[7] G. MacGougan et al., “Degraded GPS signal measurements with a stand- posed signal replicas,” Signal Process., vol. 84, no. 12, pp. 2325–2343,
alone high sensitivity receiver,” in Proc. 2002 Nat. Tech. Meet. Inst. Navig., 2004.
San Diego, CA, USA, Jan. 2002, pp. 191–204. [32] V. Krishnan and K. Chandra, Probability and Random Processes. Hobo-
[8] N. Ziedan, GNSS Receivers for Weak Signals (Ser. Artech House Space ken, NJ, USA: Wiley, 2015.
Technology and Applications). Norwood, MA, USA: Artech House, 2006. [33] M. Jakobsen, T. Pedersen, and B. Fleury, “Analysis of stochastic radio
[9] B. R. Townsend, P. C. Fenton, K. J. van Dierendonck, and D. J. R. van Nee, channels with temporal birth-death dynamics: A marked spatial point pro-
“Performance evaluation of the multipath estimating delay lock loop,” J. cess perspective,” IEEE Trans. Antennas Propag., vol. 62, no. 7, pp. 3761–
Inst. Navig., vol. 42, no. 3, pp. 503–514, 1995. 3775, Jul. 2014.
[10] J. Soubielle, I. Fijalkow, P. Duvaut, and A. Bibaut, “GPS positioning in [34] B. Krach and R. Weigel, “Markovian channel modeling for multipath mit-
a multipath environment,” IEEE Trans. Signal Process., vol. 50, no. 1, igation in navigation receivers,” in Proc. 3rd Eur. Conf. Antennas Propag.,
pp. 141–150, Jan. 2002. Mar. 2009, pp. 1441–1445.

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.
7676 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 66, NO. 9, SEPTEMBER 2017

[35] P. Hall and C. Heyde, Martingale Limit Theory and Its Application (Ser. Felix Antreich received the Diploma degree in elec-
Probability and Mathematical Statistics). New York, NY, USA: Academic, trical engineering and the Doktor-Ingenieur (Ph.D.)
1980. degree from the Munich University of Technology,
[36] S. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory Munich, Germany, in 2003 and 2011, respectively.
(Ser. Fundamentals of Statistical Signal Processing). Englewood Cliffs, From 2003 to 2016, he was an Associate Re-
NJ, USA: Prentice-Hall, 1993. searcher in the Department of Navigation, Institute of
[37] P. Schreier and L. Scharf, Statistical Signal Processing of Complex-Valued Communications and Navigation, German Aerospace
Data: The Theory of Improper and Noncircular Signals. Cambridge, U.K.: Center (DLR), Wessling, Germany. Since Septem-
Cambridge Univ. Press, 2010. ber 2016, he has been a Visiting Professor in the
[38] “NAVSTAR GPS Space Segment/User Segment L1C Interface,” Global Department of Teleinformatics Engineering, Federal
Positioning Systems Directorate, Los Angeles, CA, USA, Interface Spec- University of Ceará, Fortaleza, Brazil. From 2014
ification IS-GPS-800D, Sep. 2013. to 2016, he was a Special Visiting Researcher at the University of Brası́lia
[39] P. Stoica and T. Marzetta, “Parameter estimation problems with singular within the Programa Ciência sem Fronteiras funded by Coordenação de
information matrices,” IEEE Trans. Signal Process., vol. 49, no. 1, pp. 87– Aperfeiçoamento de Pessoal de Nı́vel Superior and Conselho Nacional de De-
90, Jan. 2001. senvolvimento Cientı́fico e Tecnológico. His research interests include sensor
array signal processing for GNSS and wireless communications, estimation
Christoph Enneking received the B.Sc. and M.Sc. theory, wireless sensor networks, and signal design for synchronization. He re-
degrees in electrical engineering from the Munich ceived the VDE Award 2011 for outstanding dissertation from Verband der Elek-
University of Technology, Munich, Germany, in 2012 trotechnik Elektronik Informationstechnik e.V., Germany in November 2011,
and 2014, respectively. as well as the Best Presentation Award at the ION GNSS 2012 and 2015 Con-
In September 2014, he joined the Institute of ference in September 2012 and 2015, respectively.
Communications and Navigation, German Aerospace
Center (DLR), Wessling, Germany. His research in-
terests include GNSS signal design, estimation the-
ory, and GNSS intra- and intersystem interference.

Authorized licensed use limited to: QSIO. Downloaded on December 17,2024 at 12:51:29 UTC from IEEE Xplore. Restrictions apply.