an author's
https://oatao.univ-toulouse.fr/23553
http://doi.org/10.1109/TITS.2018.2848461
Lesouple, Julien and Robert, Thierry and Sahmoudi, Mohamed and Tourneret, Jean-Yves and Vigneau, Willy
Multipath Mitigation for GNSS Positioning in an Urban Environment Using Sparse Estimation. (2019) IEEE
Transactions on Intelligent Transportation Systems, 20 (4). 1316-1328. ISSN 1524-9050
Multipath Mitigation for GNSS Positioning in an
Urban Environment Using Sparse Estimation
Julien Lesouple , Student Member, IEEE, Thierry Robert, Mohamed Sahmoudi,
Jean-Yves Tourneret
, Senior Member, IEEE, and Willy Vigneau
Abstract—Multipath (MP) remains the main source of error
when using global navigation satellite systems (GNSS) in a
constrained environment, leading to biased measurements and
thus to inaccurate estimated positions. This paper formulates the
GNSS navigation problem as the resolution of an overdetermined
system whose unknowns are the receiver position and speed,
clock bias and clock drift, and the potential biases affecting
GNSS measurements. We assume that only a part of the satellites
are affected by MP, i.e., that the unknown bias vector has
several zero components, which allows sparse estimation theory
to be exploited. The natural way of enforcing this sparsity
is to introduce an
1
regularization associated with the bias
vector. This leads to a least absolute shrinkage and selection
operator problem that is solved using a reweighted-
1
algorithm.
The weighting matrix of this algorithm is designed carefully as
functions of the satellite carrier-to-noise density ratio (C/N
0
)and
the satellite elevations. Experimental validation conducted with
real GPS data show the effectiveness of the proposed method as
long as the sparsity assumption is respected.
Index Terms— GNSS, multipath mitigation, sparse, LASSO,
reweighted-l1 algorithm.
I. INTRODUCTION
M
ULTIPATH (MP) is one of the most difficult error
sources that needs to be tackled for GNSS posi-
tioning [2], [3]. Indeed, MP signals are generally due to
reflections on various obstacles, and thus strongly depend
on the geometric configuration of the scene in which the
receiver is located. More precisely, in the absence of obstacle,
the receiver is not affected by MP. Conversely, when the
receiver is located close to buildings, the received GNSS
measurements are very likely to be subjected to MP. The
problem of mitigating MP effects in GNSS measurements has
received a considerable attention in the literature. MP can
be mitigated at the antenna level, by exploiting the fact that
reflections change the polarization of the received signals. As a
This work was supported
in part by the CNES and in part M3 Systems.
The Associate Editor for this paper was C. F. Mecklenbräuker.
(Corresponding author: Julien Lesouple.)
J. Lesouple is with the TéSA Laboratory, 31500 Toulouse, France (e-mail:
julien.lesouple@tesa.prd.fr).
T. Robert is with the Centre national d’études spatiales, 31400 Toulouse,
France.
M. Sahmoudi is with the ISAE, University of Toulouse, 31400 Toulouse,
France (e-mail: mohamed.sahmoudi@isae.fr).
J.-Y. Tourneret is with the TéSA, IRIT, ENSEEIHT, University of
Toulouse, 31071 Toulouse, France (e-mail: jean-yves.tourneret@enseeiht.fr).
W. Vigneau is with M3 Systems, 31410 Lavernose-Lacasse, France.
Digital Object Identifier 10.1109/TITS.2018.2848461
consequence, antennas can be designed to be more sensitive to
the right polarization [4], [5]. Another important characteristic
of reflected signals is their low or negative elevation that can
be used at the antenna level to attenuate MP signals [4], [5].
Methods using antenna arrays also exist [6]–[8]. A recent
technique combines the two latest methods mentioned [9].
MP can also be mitigated at the receiver level, by modifying
the correlator, e.g., by using narrow correlators [10], double
delta correlators [11], early late slope [12] or vision corre-
lators [13]. Other techniques work at the discriminator level,
such as the Maximum likelihood techniques based on an MP
estimating delay lock loop (MEDLL) [14]–[16], the coupled
amplitude DLL (CADLL) [17], or the Multipath Insensitive
DLL (MIDLL) [18] have also been developed for MP sig-
nals. All the previously mentioned techniques need specific
and expensive hardware that cannot always be purchased.
Mitigating MP at a measurement or position level is thus an
interesting alternative. A first solution is to take advantage of a
3D model of the environment to predict MP signals [19]–[23],
and even to combine these techniques to other sensors, such
as cameras. However, this 3D model is not always available in
practical applications. A second option is to use the informa-
tion available at the receiver, such as pseudoranges, Doppler
shifts, satellite ephemeris and C/N
0
. A widespread technique
is to smooth the code measurements with phase measurements
that are more robust to MP [24]. Other techniques consist
in exploiting different measurements from the same satellite,
for instance code and phase measurements leading to the
code minus carrier (CmC) [4], or the difference between
the measurements from two receivers leading to differen-
tial GNSS [5] or even from two different users (collabo-
rative or cooperative positioning) [25], [26]. An interesting
family of MP mitigation methods rely on statistical tests trying
to exclude or correct the faulty measurements. The receiver
autonomous integrity monitoring (RAIM) method belongs to
this class of strategies [27], [28]. More recent technique uses
a-contrario modeling for discarding bad satellites [29]. Note
that these techniques require redundant measurements, that are
not always available in urban environment, and that the user
will only be able to detect/estimate up to two faulty measure-
ments. Other techniques consider non-gaussian error terms,
such as Gaussian mixtures, Markov process [30] or Dirichlet
process mixtures [31], [32].
The point of view considered in this work is to model
the effect of MP signals on GNSS measurements as additive
biases as in [33]. These biases have then to be estimated
and subtracted from the GNSS measurements to mitigate MP
effects. Sequential Monte Carlo methods also referred to as
particle filters were investigated in [33] for this estimation.
However, these methods are computationally intensive, making
a real time implementation very complicated in practical appli-
cations. The main contribution of this paper is to exploit sparse
estimation theory to estimate these biases with a significantly
reduced computational complexity. The main hypothesis mak-
ing this theory applicable is that many satellites are not
affected by MP making the biases sparse with respect to
number of received measurements. Note that sparse estimation
for mitigating multipath has already been considered in Radar
theory [34], [35], and that sparse assumption was also consid-
ered for GNSS applications in [36]. However, the proposed
approach is different. It results from the application of a
penalized least squares approach method taking advantage of
the recent developments in sparse estimation theory. Since the
measurement equation is linear with respect to the state vector,
we estimate it directly from the data and the biases term in
order to form a profile likelihood that is used to estimate
the MP biases. This sparse estimation formulation avoids to
consider an augmented state vector for bias estimation as
in [36]. Instead, the bias vector is directly obtained by the
minimization of a a penalized least-squares criterion resulting
from a sparse MP prior.
This paper is organized as follows: Section II summarizes
some basic principles on satellite navigation, describing how
measurements (code measurements and Doppler rates) are
related to the state vector (position, velocity) and the possible
MP biases. This section also recalls the Kalman filtering steps
that will be used to track the receiver position. Section III
summarizes the main ideas of sparse estimation theory and
the main estimation methods (LASSO, reweighted-
1
and
generalized LASSO) that can be used to estimate sparse
vectors. Section IV presents our contributions, i.e., how to
estimate MP biases using GNSS measurements with some
sparsity constraints. More precisely, we propose to consider a
linearized navigation equation and to assign some sparsity con-
straints on the biases possibly affecting GNSS measurements.
The positioning problem is then formulated as a penalized
least squares problem with an
1
regularization inspired by
the reweighted-
1
algorithm [37]. However, contrary to [37],
the weighting matrix used in this work is designed using
important information available at the receiver, based on the
value of the carrier-to-noise density ratio (C/N
0
)andthe
satellite positions. Experiments are presented in Section V
comparing the proposed algorithm with other navigation strate-
gies (Kalman filter, robust Kalman filter, classical LASSO,
coded filter [36]). Conclusions and future work are reported
in Section VI.
II. GNSS F
UNDAMENTALS
A. Observation Model
The navigation problem considered in GNSS consists in
estimating the position of a receiver from signals sent by
different satellites. More precisely, measuring the propagation
delay between the receiver and a given satellite, the receiver is
able to build a so-called pseudorange defined as follows [38]
ρ
i
=x
i
− x
u
2
+ b + ε
i
(1)
where
• ρ
i
denotes the pseudorange between the receiver and the
ith satellite, with i ∈{1,...,N}, N being the number of
in-view satellites,
• x
u
= (x, y, z)
T
is the receiver position to be estimated,
• x
i
= (x
i
, y
i
, z
i
)
T
is the known ith satellite position,
• x
i
− x
u
2
=
(x
i
− x)
2
+ (y
i
− y)
2
+ (z
i
− z)
2
is the
distance between the user and the ith satellite,
• b is the receiver clock bias, common to all measure-
ments (hence the name of pseudorange),
• ε
i
is the error term associated with the ith propagation
canal (modeling ionospheric delay, tropospheric delay,
satellite clock bias, satellite position uncertainty, MP and
receiver noise).
A classical way of estimating the receiver position from
N measurement equations as defined in (1) is to use an
iterative algorithm, which linearizes (1) around the previous
computed position. The resulting linearized problem for GNSS
navigation can be classically expressed as [39] and [40]
y
p
= Gx + m
p
+ n
p
(2)
with
• y
p
∈ R
N
the difference between the measured and
estimated pseudoranges (the subscript p is used for
pseudoranges),
• G ∈ R
N×4
the Jacobian matrix associated with the
linearized system,
• x ∈ R
4
the difference between the state (position and
receiver clock bias) estimated at the previous position and
the current state value,
• m
p
∈ R
N
an error term due to the possible presence of
MP affecting the pseudoranges (we assume that all errors
except MP have been corrected),
• n
p
∼ N (0, R
p
) ∈ R
N
a zero-mean Gaussian noise vector
with covariance matrix R
p
.
The expression of the matrix G can be found in many
textbooks such as [39] and is recalled here for completeness
G =
⎡
⎢
⎢
⎢
⎢
⎢
⎣
a
1
1
a
1
2
a
1
3
1
a
2
1
a
2
2
a
2
3
1
.
.
.
.
.
.
.
.
.
.
.
.
a
N
1
a
N
2
a
N
3
1
⎤
⎥
⎥
⎥
⎥
⎥
⎦
(3)
with
a
i
1
, a
i
2
, a
i
3
T
=
x
0
− x
i
x
0
− x
i
(4)
where x
0
is the point around which (1) has been linearized.
After differentiating (1), the following result can be
obtained [39]
˙ρ
i
=
a
i
1
, a
i
2
, a
i
3
( ˙x
u
−˙x
i
) +
˙
b +˙ε
i
(5)
which leads to the following linear equation that is associated
with pseudorange rates and is very similar to (2)
y
r
= G ˙x + m
r
+ n
r
(6)
where
• y
r
∈ R
N
is the difference between the measured and
estimated pseudorange rate (the subscript r is used to
indicate pseudorange rates),
• ˙x ∈ R
4
is the difference between the state (velocity and
receiver clock drift) estimated at the previous position and
the current state value,
• m
r
∈ R
N
is an error term due to the possible presence
of MP affecting the pseudorange rate,
• n
r
∼ N (0, R
r
) ∈ R
N
is a zero-mean Gaussian noise
vector with covariance matrix R
r
.
The main idea investigated in this paper is to exploit the
property that few satellites are suffering from MP such that
m = (m
p
, m
r
)
T
is a sparse vector, i.e., a vector containing a
lot of zeroes. As a consequence, sparse estimation theory can
be used to estimate m and x jointly. Note that the terms m
p
and m
r
in (2) and (6) are both resulting from the possible
presence of MP. However, as the receiver computes these
two measurements differently, there is a priori no relation
between these two terms. Before providing more details about
the proposed sparse estimation method, we recall some basic
elements about the extended Kalman filter that will be used
for state estimation.
B. Extended Kalman Filter for GNSS
The EKF considered in this work is very classical and
has been studied in many papers or textbooks including,
e.g. [41, p. 195] or [42]. The state vector at time instant k,
denoted as s
k
(for k = 1,...,K ,whereK is the sample size),
contains the receiver position and velocity expressed in the
ECEF frame, and the clock bias and drift (derivative of
clock bias), i.e.,
s
k
=
x
k
, ˙x
k
, y
k
, ˙y
k
, z
k
, ˙z
k
, b
k
,
˙
b
k
T
(7)
where
• the subscript k denotes the time instant,
• (x
k
, y
k
, z
k
)
T
is the receiver position at time instant k,
• ( ˙x
k
, ˙y
k
, ˙z
k
)
T
is the receiver velocity at time instant k,
• b
k
is the clock bias at time instant k,
•
˙
b
k
is the clock drift at time instant k.
The relationship between the state vectors s
k
and s
k−1
,
is defined by a propagation equation that is used to design
the Kalman filter. Let’s denote
C
k
=
1 t
k
01
(8)
and
F
k
=
⎡
⎢
⎢
⎣
C
k
0
2×2
0
2×2
0
2×2
0
2×2
C
k
0
2×2
0
2×2
0
2×2
0
2×2
C
k
0
2×2
0
2×2
0
2×2
0
2×2
C
k
⎤
⎥
⎥
⎦
. (9)
The propagation equation considered in this work is
s
k
= F
k−1
s
k−1
+ u
k−1
(10)
with
• F
k
the state transition matrix at time instant k,
• t
k
the duration between two consecutive measurements,
• u
k−1
∼ N (0, Q
k−1
) the process noise at time instant
k − 1, supposed to be zero-mean Gaussian with known
covariance matrix Q
k−1
(see [43] for a closed-form
expression).
The receiver can have access to various measurements depend-
ing on the application. In this paper, we focus on the code
pseudorange (ρ
i
) and the Doppler rate ( f
i
), even if the
proposed methodology is quite general and can be applied
to any kind of measurement affected by sparse additive
biases. Assuming that the receiver has access to N satellite
signals (N > 4), the measurement vector can be defined as
ρ
˙ρ
=
ρ
1
,...,ρ
N
, ˙ρ
1
,..., ˙ρ
N
T
(11)
where ˙ρ
i
=−λ
L
1
f
i
is the ith pseudorange rate, with
λ
L
1
the wavelength of the received signal. As explained in
Section II-A, the measurements are related to the state vector
thanks to (2) and (6) and can be concatenated into a single
equation defined as
y
k
= H
k
x
k
+m
k
+ n
k
(12)
where
• y
k
∈ R
2N
contains the difference between the actual and
predicted pseudoranges and pseudorange rates at time k,
• H
k
is the joint observation matrix for pseudoranges and
pseudorange rates at time k (see technical report [43] for
a closed-form expression),
• x
k
= s
k
−ˆs
k
∈ R
8
is the difference between the state
vector (receiver position, velocity, clock bias, and clock
drift) estimated at the previous position and the current
state vector at time k,
• m
k
∈ R
2N
is a bias term due to the possible presence of
MP at time k,
• n
k
∼ N (0, R
k
) ∈ R
2N
is a zero-mean Gaussian noise
vector at time k with covariance matrix R
k
(see technical
report [43] for a closed-form expression).
The classical EKF used to estimate the state vector from the
state equation (10) and the measurement equations (1) and (5)
can be summarized as follows
1) one step prediction
ˆs
k|k−1
= F
k−1
ˆs
k−1|k−1
2) one step prediction error covariance
P
k|k−1
= F
k−1
P
k−1|k−1
F
T
k−1
+ Q
k−1
3) Kalman gain
K
k
= P
k|k−1
H
T
k
H
k
P
k|k−1
H
T
k
+ R
k
−1
4) state estimation
ˆs
k|k
=ˆs
k|k−1
+ K
k
y
k
5) state error covariance
P
k|k
=
(
I − K
k
H
k
)
P
k|k−1
.
In this paper, we have adopted the square root implementation
of the Kalman filter discussed in [41, p. 181] for its known
robustness.
Remark: In order to compensate for relativistic effects,
time group delays and clock biases for the different satellites,
we followed [44]. Ionospheric delays and their derivatives
were compensated using the well known Klobuchar
model [45]. Zenith values were computed as in [46] and
mapped with the Niell global mapping function [47] to remove
the effects of tropospheric delays and their derivatives. Based
on these compensations, the residual error was considered as
Gaussian centered with variance σ
2
UERE
= σ
2
Ephemeris
+ σ
2
Iono
+
σ
2
Tropo
+ σ
2
Clock
for the pseudoranges and ˙σ
2
UERE
=
˙σ
2
Ephemeris
+˙σ
2
Iono
+˙σ
2
Tropo
+˙σ
2
Clock
for pseudorange rates,
see [39, p. 326] for justification and [39, p. 273] for typical
values. However, other additive biases due for instance to
errors in ionospheric or tropospheric corrections can also
affect the received measurements as additive biases. The
proposed method will include these errors in the vector m
k
and will correct them, providing the sparsity assumption
is satisfied.
III. S
PARSE ESTIMATION THEORY
This section recalls the principles of sparse estimation
theory and the main algorithms that can be used to estimate
sparse vectors. Assume that we have a vector of measurements
˜y ∈ R
2N
defined as ˜y =
˜
Hθ +˜n,where
˜
H ∈ R
2N×q
is
a known regression matrix, θ ∈ R
q
is an unknown sparse
vector (to be estimated) and ˜n ∈ R
2N
is an unknown error
term.
1
A classical way of estimating θ from the observed
measurement vector ˜y is to consider a data fidelity term
1
2
˜y−
˜
Hθ
2
2
penalized by an additive regularization promoting
the sparsity of θ. One can think of defining this additive
regularization as the
0
pseudo-norm of θ defined by
θ
0
= #{θ
i
= 0, i = 1,...,q}. (13)
This problem can be formulated in different ways [48] and we
choose the unconstrained one defined by
ˆ
θ = argmin
θ∈R
q
1
2
˜y −
˜
Hθ
2
2
+ λθ
0
(14)
where λ ∈ R is a fixed constant referred to as regularization
parameter. However, the problem (14) is NP-hard and non-
convex. Therefore, many relaxations have been proposed in
the literature to bypass this difficulty, as summarized below.
A. The LASSO Problem
A classical way of estimating a sparse vector from a linear
regression is to replace the
0
pseudo-norm in (14) by an
1
norm, i.e., to consider the so-called LASSO problem [49]
argmin
θ∈R
q
1
2
˜y −
˜
Hθ
2
2
+ λθ
1
. (15)
However, algorithms used to solve this problem can provide
solutions that are far from the solution of (14) [37]. This
has motivated the study of many different sparse estimation
strategies described in the next sections.
1
The meaning of the different vectors ˜y, θ , ˜n in the GNSS context will be
clarified in Section IV.
B. The Reweighted-
1
l1 and the Generalized LASSO
Candès et al. [37] investigated a so-called reweighted-
1
method defined as follows
argmin
θ∈R
q
1
2
˜y −
˜
Hθ
2
2
+ λWθ
1
(16)
where W ∈ R
q×q
is a diagonal weighting matrix. Ideally,
the weights contained in W should be inversely proportional
to the magnitude of the true unknown vector θ
0
, i.e., such that
w
i
=
1
|θ
0,i
|
,θ
0,i
= 0,
∞,θ
0,i
= 0.
(17)
However, this weight definition cannot be used in practice
since θ
0
is an unknown vector. One solution proposed by
Candès et al. [37] is summarized in the following iterative
algorithm
1) set = 0andw
(0)
i
= 1, i = 1,...,q,
2) solve the problem
θ
()
= argmin
θ∈R
q
1
2
˜y −
˜
Hθ
2
2
+ λW
()
θ
1
, (18)
3) update the weights as
w
(+1)
i
=
1
|θ
()
i
|+ε
, (19)
4) end when reaches a maximum value
max
.
This algorithm requires to adjust the two parameters ε and
max
. An adhoc choice was suggested in [37], namely ε = 0.1
and
max
= 2.
A generalized version of (15) referred to as “generalized
LASSO” was introduced in [50]
argmin
θ∈R
q
1
2
˜y −
˜
Hθ
2
2
+ λWθ
1
(20)
where W ∈ R
p×q
is an appropriate penalty matrix that is
not necessarily square (p denotes the number of constraints
associated with the unknown parameter vector θ) and needs
to be specified by the user. Of course, when p = q and W is
diagonal, the generalized LASSO reduces to the reweighted-
1
method.
IV. A N
EW MULTIPATH MITIGATION METHOD FOR GNSS
A. Problem Formulation
The proposed MP mitigation method assumes that the bias
vector m = (m
p
, m
r
)
T
resulting from (2) and (6) is sparse.
Exploiting this sparsity property, we propose to solve the
following problem
argmin
x,m
1
2
y − Hx− m
2
2
+ λWm
1
(21)
in order to detect and correct measurements affected by MP,
i.e., measurements affected by the presence of additive biases.
Note that these corrected measurements will be used as input
of the EKF presented in II-B. In order to obtain a formulation
similar to (16), it is interesting to note that the minimization
