1194 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 3, MARCH 2005
ML Estimator and Hybrid Beamformer for Multipath
and Interference Mitigation in GNSS Receivers
Gonzalo Seco-Granados, Member, IEEE, Juan A. Fernández-Rubio, Member, IEEE, and
Carles Fernández-Prades, Student Member, IEEE
Abstract—This paper addresses the estimation of the code-phase
(pseudorange) and the carrier-phase of the direct signal received
from a direct-sequence spread-spectrum satellite transmitter. The
signal is received by an antenna array in a scenario with interfer-
ence and multipath propagation. These two effects are generally
the limiting error sources in most high-precision positioning ap-
plications. A new estimator of the code- and carrier-phases is de-
rived by using a simplified signal model and the maximum like-
lihood (ML) principle. The simplified model consists essentially of
gathering all signals, except for the direct one, in a component with
unknown spatial correlation. The estimator exploits the knowledge
of the direction-of-arrival of the direct signal and is much simpler
than other estimators derived under more detailed signal models.
Moreover, we present an iterative algorithm, that is adequate for a
practical implementation and explores an interesting link between
the ML estimator and a hybrid beamformer. The mean squared
error and bias of the new estimator are computed for a number
of scenarios and compared with those of other methods. The pre-
sented estimator and the hybrid beamforming outperform the ex-
isting techniques of comparable complexity and attains, in many
situations, the Cramér–Rao lower bound of the problem at hand.
Index Terms—Adaptive arrays, adaptive estimation, array
signal processing, beamforming, beam steering, bias, calibration,
Code division multiaccess, Cramér–Rao bounds, delay estima-
tion, early-late estimator, Global positioning system, GPS, GPS
positioning, GPS receiver, GPS signal, Gold codes, interference
suppression, iterative methods, low-complexity constraint, max-
imum likelihood estimation, maximum likelihood estimator, MLE,
modeling, multipath channels, multipath environment, multipath
propagation, pseudo random codes, Radio Navigation, radio re-
ceivers, reflected components, Satellite navigation systems, signal
model, simulations, single-path environment, standard deviation,
Time of arrival estimation.
I. INTRODUCTION
T
HE term Global Navigation Satellite Systems (GNSS) is
a generic expression referring to any system that enables
the calculation of the user position based on signals trans-
mitted by a constellation of satellites. At the present time, the
Global Positioning System (GPS) is the only fully operational
system. The European augmentation of GPS (EGNOS) and a
Manuscript received May 16, 2003; revised February 17, 2004. This work was
supported in part by the Spanish Science and Technology Commission under
project TIC2001-2356-C02-01 and by the Spanish Ministry of Education under
Grant FPU AP2000-3893. The associate editor coordinating the review of this
manuscript and approving it for publication was Dr. Chong-Yung Chi.
G. Seco-Granados is with the Electrical Engineering Department, Euro-
pean Space Agency (ESA), 2200 AG Noordwijk, The Netherlands (e-mail:
gonzalo@ieee.org).
J. A. Fernández-Rubio and C. Fernández-Prades are with the Signal Theory
and Communications Department, Polytechnic University of Catalonia, 08034
Barcelona, Spain (e-mail: juan@gps.tsc.upc.edu; carlos@gps.tsc.upc.edu).
Digital Object Identifier 10.1109/TSP.2004.842193
new system named GALILEO will be operational in the next
few years. All GNSS share the same operating principle: The
receiver position is computed based on the distances between
the receiving antenna and a set of satellites, and the receiver
determines these distances by measuring the propagation time
of the signals transmitted by the satellites. This propagation
time can be obtained from the delay (referred to as pseudo-
range or code-phase) of the complex envelope and from the
carrier-phase [1].
The surprising evolution of GNSS applications has led to
stringent requirements for GNSS receivers, particularly in
regard to their accuracy. Augmentations such as differential
operation help to reduce or eliminate many sources of errors
(e.g. common-mode atmospheric, orbit-, and satellite-induced
errors), but multipath remains the dominant error source in
most high-precision applications and is the limiting factor in
achieving the ultimate GNSS accuracy [2]. Due to the operating
principle of the GNSS, only the direct signal [which is also
called the line-of-sight signal (LOSS)] bears useful information
about the distance between the receiver and the satellite. Sig-
nificant research and development efforts have been devoted to
the mitigation of multipath effects, and a number of techniques
have been proposed so far. They may be classified according
to a variety of criteria, e.g., real-time versus post-processing
techniques and multiple versus single antenna techniques.
GNSS are also subject to external interferers, which have to
be cancelled in the receiver in order to make GNSS adequate
for many safety-critical applications, such as aircraft automatic
guidance and landing systems. Several methods can be used to
mitigate narrowband interferences in single-antenna receivers,
which are usually based on a linear interpolator-subtracter
structure [3]. However, in general, single-antenna methods
cannot combat wideband interferences.
Errors in the pseudorange and carrier-phase measurements
produced by the multipath propagation have been studied in
[4]–[7], among others. Only reflections correlated with the di-
rect signal, which are usually referred to as coherent multipath,
cause these errors; this is the type of reflections considered in
this paper. Their main characteristic is that their relative delays
with respect to the LOSS are on the order of or smaller than the
inverse of the signal bandwidth. For instance, in a GPS receiver
employing a delay locked loop (DLL), which is the synchro-
nization method used in the vast majority of GNSS receivers,
the multipath components may bias the pseudoranges in several
tens or even a hundred of meters, and at the same time, they
hamper the ambiguity resolution process needed for carrier-
phase ranging. The bias in the carrier phases may reach some
1053-587X/20.00 © 2005 IEEE
SECO-GRANADOS et al.: ML ESTIMATOR AND HYBRID BEAMFORMER FOR MULTIPATH AND INTERFERENCE MITIGATION 1195
centimeters (with the favorable assumption of perfect ambiguity
resolution, which is not likely the case in a multipath scenario).
The most widespread multipath-mitigation techniques are those
based on modifications of the conventional DLL, which are real-
time signal processing methods. Some of these techniques are
the Narrow Correlator DLL [7], the Multipath Estimating DLL
(MEDLL) [8], the Pulse Aperture Correlator (PAC) (patented by
NovAtel Inc.), and the Edge Correlator and (Enhanced) Strobe
Correlator [9], [10] (patented by Ashtech Inc.). These single-
sensor techniques are able to discriminate the LOSS from the
reflections only in the temporal domain, and hence, their perfor-
mance is still limited for many precise applications. A number
of post-processing techniques for multipath mitigation have also
been proposed [6]. Since most of these techniques require data
recording for several minutes or hours, they do not work in
real-time and are restricted to a small number of applications.
On the other hand, spatial filtering is probably the most ef-
fective approach to combat both interference and multipath. Un-
like in communication systems, the potential benefits of antenna
arrays in navigation systems have not been investigated thor-
oughly. The use of antenna arrays in GNSS has been centered
on interference mitigation; see, for example, [11]–[15]. In these
works, the array processor operates directly with the received
signals, and conventional array processing techniques, such as
the minimum-variance or the power-inversion beamformers, are
applied. Combating the multipath propagation with antenna ar-
rays is a much more powerful but also involved approach, and it
has hardly been studied in the context of GNSS. Some works in
this direction are [16]–[19]. One of the major difficulties in em-
ploying an antenna array in a multipath environment is that con-
ventional array processing techniques completely fail because
of the high degree of coherence existing between the LOSS
and the reflections. The maximum likelihood (ML) approach
has been used for synchronization in communications (see, e.g.,
[20] for single-antenna and [21] and [22] for multiple-antenna
systems) but never exploiting the particularities of GNSS.
The goal of this paper is to present real-time (not post- or
data-processing) array processing techniques that can be imple-
mented in a GNSS receiver in order to mitigate the effects of
interferences and any kind of multipath.
1
The ML estimator of
the pseudorange and the carrier-phase derived from a simpli-
fied model of the received signals is proposed and analyzed. The
equivalence with a hybrid beamformer is shown, which leads to
an iterative algorithm.
II. S
IGNAL MODEL
This section presents the formulation of the signals received
by the antenna array and justifies the simplified model on which
the rest of the paper is based.
A. Description of the Received Signals
Let us consider that an arbitrary
-antenna array receives the
signal transmitted by a given GNSS satellite. Assume also that,
besides the LOSS,
reflections of the GNSS signal impinge
1
In the GNSS literature, it is normal to differentiate between specular (a few
reflections produced by smooth surfaces) and diffuse (a large number of weaker
reflections) multipath.
on the array. Then, the complex baseband representation of the
signal received at the array is the
vector
2
(1)
where
is the sampling period, and are the spatial signa-
ture and the time-delay of the
th component, and models
the thermal noise and all other interference. The subscript 0
stands for the LOSS. Hence,
, when expressed in length units,
is the pseudorange, and the phases of the elements of
are
the carrier-phase observables at each antenna. The underlying
analog signal
represents the contribution of one the signals
transmitted by a GNSS satellite; it is considered to be a Direct-
Sequence Spread-Spectrum (DS-SS) signal since all present and
planned navigation systems use this modulation format [1], [23].
Therefore, the signal
can be expressed as
(2)
where
is the spreading waveform or pseudo-noise code pos-
sibly after some kind of preprocessing. The sequence of sym-
bols
forms the navigation message of the satellite, essen-
tially bearing ephemerides information. The symbol period
comprises chips, each of duration , i.e., . Note
that
is not necessarily equal to the signal transmitted by
the satellite since in general, the signals can be preprocessed
in the receiver before applying any array processing algorithm.
In our case, the need for this preprocessing arises from the fact
that the received GNSS signals are well below the noise floor. If
the array processor operated directly on the received signals, it
would not be possible to infer any spatial information about the
reflected replicas. That is to say, the array would be insensitive
to the direct and reflected signals, and it would be able to cancel
only powerful interferences, as in [12]–[15], because the con-
tribution to the spatial correlation matrix would be dominated
by the noise and interference. Since our overall objective is to
use the spatial dimension to cancel both multipath components
and interference, some kind of preprocessing that increases the
signal-to-noise ratio (SNR) of the LOSS, and its reflections, is
mandatory. When dealing with DS-SS signals, the standard ap-
proach to attain this purpose is to perform the despreading (i.e.,
simply to correlate the received signal with a local replica of the
pseudo-noise code). In any case, the advantage of the techniques
proposed in the following section is that they are independent of
the particular way
is obtained.
An assumption needed by the techniques proposed in this
paper is that
is known. This assumption does not represent
any limitation because of the following.
• The underlaying shaping pulse in
(i.e., )is
known at the receiver (it is a design parameter of the
system).
• The techniques can be applied to portions of
that
span one symbol interval, and in this case, the knowl-
edge of the transmitted symbol is not needed.
2
Vectors are denoted by lowercase bold letters, and matrices are denoted by
uppercase bold letters.
1196 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 3, MARCH 2005
• Only when the observed portion of spans more
than one symbol interval, the symbol transitions need
to be known. Even this situation is not a handicap since
GNSS transmits training sequences [1], which can be
used, for instance, to adapt the weights of the hybrid
beamformer presented in Section V. The key issue in
GNSS and in this paper is accurate synchronization
and not data detection, which usually is not a problem
in GNSS, even in the presence of strong interference
and without using training sequences, if appropriate
techniques are employed (e.g., [3], [11]–[15]). This is
also facilitated by the large redundancy present in the
system: The same parts of the message are transmitted
at different time periods and by different satellites.
The standard “narrowband assumption” used in many array
signal processing problems has been made in writing (1). This
is well justified because the bandwidth of the GNSS signals is on
the order of few megahertz, and the carrier frequency is between
1 and 2 GHz.
B. Simplified Signal Model
A number of techniques that estimate the unknown parame-
ters of the detailed (as opposed to simplified) signal model in
(1) have been developed. Many of them assume that the spa-
tial signatures are parameterized by the corresponding direc-
tions-of-arrival (DOA), as in [24] and [25]. While these methods
may exploit the full space-time structure of the signals, they in-
volve the optimization of multidimensional nonlinear cost func-
tions. Only a few cases that resort to particular array configura-
tions allow a closed-form estimation of the delays and DOAs,
e.g., [25]. To obtain simpler criteria, most methods presume that
the noise
is spatially white, which makes them incapable
of mitigating directional interferences. The large computational
load of the previous techniques can be alleviated by using an un-
structured parameterization of the spatial signatures, which are
modeled as arbitrary deterministic vectors, as in [18] and [21].
The large complexity of all the methods enumerated above,
which are based on the model (1), is mainly due to the fact that
an important effort is devoted to estimating certain parameters
that are not of interest in a GNSS receiver, such as the DOAs
and/or the time-delays of the reflections. This justifies the search
of a simpler technique allowing us to estimate only the relevant
parameters. Such a technique should, however, be capable of
using the diversity introduced by the antenna array in order to
discriminate the signals in the spatial dimension. The proposed
technique will be based on the use of a simplified signal model.
A particularity of the GNSS systems is that the receiver has
accurate estimates of its own position as well as that of the satel-
lite. Therefore, assuming that the antenna array is calibrated, it
is possible to know the spatial signature of the LOSS up to a
scaling factor, since the DOA of the direct signal can be com-
puted beforehand. We will consider that the following relation
holds:
(3)
where
is the known steering vector of the LOSS, and
is an unknown complex amplitude. This is a rather common
Fig. 1. CRB for the LOSS time-delay with and without a priori knowledge of
a
. Parameters as in Fig. 3 except for
=5
.
assumption in many GPS-related papers [11]–[14], [19]. The
attitude of the antenna array is also needed for the computation
of
, and this information is available in many applications, as
in the static receivers present in the large number of differential
reference stations. In a dynamic environment, it usually requires
the use of an attitude sensor or data from the navigation unit.
Even in this case, in aeronautical applications, the attitude is
mostly available. Furthermore, even in very adverse multipath
scenarios, the inaccuracies of the satellite and receiver positions,
which are at most on the order of a few hundred meters, result
in negligible errors in the calculation of the DOA of the LOSS
because the distance between the receiver and the satellite is
larger than 20 000 km.
Another justification for the simplified signal model comes
from the comparison of the two Cramér–Rao bounds (CRBs)
for the LOSS time-delay estimate obtained under the detailed
signal model (1) when
is known and when it is unknown.
Fig. 1 contains one particular example of these CRBs (the sce-
nario parameters are described in Section VI-A together with
the figure caption), but the conclusions are general. Both CRBs
are very close to each other and are indistinguishable in most
cases. This implies that the knowledge of
does not provide
essential information for the detailed signal model and does not
allow us to improve significantly the performance of the ML es-
timator for that model. Moreover, that ML estimator would not
result in a computational efficient implementation. Therefore,
instead of using the a priori knowledge of
together with the
detailed signal model, it is preferred to exploit this additional
information in order to simplify the signal model itself, that is
to say, to use a simple model (and hence simple estimators) that
would be unfeasible without the knowledge of
.
In the simplified signal model, the actual received signal (1)
is expressed as the addition of only two terms
(4)
together with the assumption in (3) and
being known. The
first term is the LOSS contribution, and the second one is an
equivalent noise that includes the contribution of all the unde-
sired signals, i.e., reflections, interferences, and thermal noise.
SECO-GRANADOS et al.: ML ESTIMATOR AND HYBRID BEAMFORMER FOR MULTIPATH AND INTERFERENCE MITIGATION 1197
This model is appropriate for our goals because only the desired
parameters, that is,
and , remain explicitly. This fact pro-
vides additional support for separating the contribution of the
direct signal from that of the rest of the signals.
Vector
is modeled as a complex, circularly-symmetric,
zero-mean Gaussian process that has an unknown and arbitrary
spatial correlation matrix
:
3
(5)
For simplicity, the process is assumed to be temporally white.
Matrix
is intended to model the directional or spatial charac-
teristics of both interference and multipath components. Indeed,
it is the fact that the correlation matrix is unknown and has to
be estimated that will make the estimator capable of using the
diversity introduced by the antenna array to discriminate the sig-
nals in the spatial dimension, despite the approximate modeling
of
. The simplification of the signal model comes clearly at
the expense of a certain mismatch between the model and the
actual received signals because the previous assumptions about
the equivalent noise need not be satisfied in a real scenario. For
instance, the zero-mean assumption of
(or, equivalently, the
assumption that signal and noise terms are uncorrelated) is vi-
olated when
contains the contribution of (coherent) reflec-
tions. However, this mismatch also has some positive effects on
the performance of the estimators, as further discussed in Sec-
tion VI-C. The Gaussian hypothesis for the equivalent noise is
of interest because it allows an analytical treatment of the ML
estimator and easily captures in
the spatial structure of the
multipath and interference.
All in all, the use of the simplified signal model is justified by
the fact that it allows derivation of simple estimators, whose per-
formance is excellent, as shown in Section VI. Indeed, the same
idea has been applied successfully to synchronization in mo-
bile communication systems with multiuser interference (see,
e.g., [22], [26], and [27]) and to Doppler and DOA estimation
in radar systems (see, e.g., [28]).
The
samples of (4) collected during an observation interval
can be arranged into the following
matrix:
(6)
where is formed identically to , and we have defined the
signal vector
(7)
III. ML C
ODE AND CARRIER-PHASE ESTIMATORS
The problem addressed in this section may be stated as fol-
lows: Given the collection of data
modeled by (6), the vector
, and the signal , estimate the unknown parameters , ,
and
. To this end, the ML principle is going to be applied. The
inverse of the likelihood function of the data is
4
Tr (8)
3
The transpose, conjugate, and conjugate transpose operations are designated
by
(
1
)
,
(
1
)
, and
(
1
)
, respectively.
is equal to 1 if
n
=
l
and 0 otherwise.
E
f1g
is the mathematical expectation.
4
Throughout the paper, all parameter-independent additive or multiplicative
constants in likelihood functions are neglected. Tr
f1g
denotes the trace opera-
tion.
j1j
denotes the determinant for matrices and the absolute value for scalars.
where
(9)
The ML estimates of the parameters (
, , and ) are
those values that minimize (8). The value of
that nulls the
gradient of (8) with respect to
is given by
(10)
where we have assumed that
in order for
to be invertible with probability one.
Let us define the following sample correlations:
(11)
(12)
Matrix
is an unstructured estimate of the noise cor-
relation matrix
. It is called unstructured because, unlike
, it does not use the knowledge of the spatial signa-
ture of the LOSS. Substituting (10) into (8) yields the following
concentrated inverse likelihood function:
(13)
(14)
(15)
Equation (15) stems directly from the following property of the
determinant:
, which is valid for matrices
with the appropriate dimensions. A straightforward minimiza-
tion of (15) with respect to
results in the ML estimate of the
LOSS amplitude:
(16)
Thanks to the invariance principle of the ML estimates, the ML
estimate of the carrier phase is directly the phase of
.
After plain but lengthy calculations, in which (16) is substituted
into (15) and
is expanded using the matrix inversion
lemma, the criterion in (15) can be expressed as a function of
only:
(17)
(18)
1198 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 53, NO. 3, MARCH 2005
Therefore, the ML estimate of the time delay is
(19)
The computation of the ML estimates involves only the search of
the maximum of the one-dimensional function
; there-
fore, this method is a good candidate for implementation in a
real receiver. Moreover, this ML estimator is applicable in the
same way in the presence of any type (specular or diffuse) of
multipath, which is an important advantage with respect to the
methods based on the detailed model in (1). If the term
is temporally correlated, the presented estimator is not the ML
one, but it can be applied in the same way without expecting in
general significantly worse performance because the estimator
maintains the ability to cancel interferences in the spatial do-
main. Obviously, better performance can achieved if the tem-
poral correlation is exploited, following, for instance, the ap-
proaches in [22] or [29], but this comes at the expense of a much
higher computational complexity. No general statements about
the potential improvement are possible since it depends on the
spectral shape of
and .
The previous ML estimator will be compared in Section VI
with three other methods, which are outlined below. The first
one consists in assuming that the steering vector of the LOSS
is arbitrary and unknown along with the model in (4), (5). The
minimization of (15) with respect to
is trivial, and the
ML time-delay estimate obtained using only a temporal refer-
ence is
(20)
The second method relies on the same simplified model as the
ML estimator proposed in this paper but with the additional as-
sumption that the noise is spatially white. In this method, vector
is again considered to be known. The derivation of the ML
estimates when
is replaced by is simple and yields
(21)
(22)
The last approach involves the spatial filtering of the received
signals using the classical minimum-variance or Capon beam-
former (MVB):
(23)
At first glance, this may seem to be a logical solution, and it
has been proposed for the problem under consideration in some
works, such as [12] and [14]. If the ML criterion is applied to
the output signal of the beamformer
, the re-
sulting estimates are
(24)
(25)
It is very interesting to observe that the ML criterion proposed
in (19) can be expressed as a function of the ML criterion using
only temporal information and the cost function based on the
MVB:
(26)
In Sections IV-C and VI-C2, we will dwell on the implications
of this expression.
IV. P
ERFORMANCE ANALYSIS
A complete characterization of the performance of the pro-
posed ML estimator can only be obtained by simulation since
analytical expressions of the bias and the variance under all pos-
sible circumstances (i.e., under the validity of the detailed model
(1) and under the validity of the simplified model (4)) do not yet
exist. However, it is possible to carry out an analytical study
of certain characteristics, and together with the simulation re-
sults, the study reveals important features of the ML estimator.
First, the asymptotic behavior of the cost function is addressed;
second, the CRBs for both models are presented; and last, the
effect of errors in the
a priori knowledge of
is analyzed.
A. Asymptotic Properties of the ML Cost Function
It can be shown that if the model (3)–(5) holds exactly, the
previous ML estimators of
(19), (16), and (10) are con-
sistent. The proof follows closely the proof in [28]. Next, we
analyze the asymptotic (hereinafter, in the number of samples
) expression of the cost function (19) when the actual received
signal corresponds to the detailed signal model (1). We use the
following notation:
; is a ma-
trix whose
, th element is , where is the
asymptotic autocorrelation of
and is the th element of
; ; denotes the true value of the cor-
responding parameter; and
. The
asymptotic expression of cost function
, shown in (19),
is
(27)
with probability one. This cost function is not maximum at
, in general; however, further insight can be gained by as-
suming a high SNR, which is the usual situation after prepro-
