Research Article  Open Access
Carsten Fritsche, Anja Klein, "Nonlinear Filtering for Hybrid GPS/GSM Mobile Terminal Tracking", International Journal of Navigation and Observation, vol. 2010, Article ID 149065, 17 pages, 2010. https://doi.org/10.1155/2010/149065
Nonlinear Filtering for Hybrid GPS/GSM Mobile Terminal Tracking
Abstract
The Global Positioning System (GPS) has become one of the stateoftheart location systems that offers reliable mobile terminal (MT) location estimates. However, there exist situations where GPS is not available, for example, when the MT is used indoors or when the MT is located close to high buildings. In these scenarios, a promising approach is to combine the GPSmeasured values with measured values from the Global System for Mobile Communication (GSM), which is known as hybrid localization method. In this paper, three nonlinear filters, namely, an extended Kalman filter, a RaoBlackwellized unscented Kalman filter, and a modified version of the recently proposed cubature Kalman filter, are proposed that combine pseudoranges from GPS with timing advance and received signal strengths from GSM. The three filters are compared with each other in terms of performance and computational complexity. Posterior CramérRao lower bounds are evaluated in order to assess the theoretical performance. Furthermore, it is investigated how additional GPS reference time information available from GSM influences the performance of the hybrid localization method. Simulation and experimental results show that the proposed hybrid method outperforms the GSM method.
1. Introduction
In the past few years, there is an increased interest in wireless location systems offering reliable mobile terminal (MT) location estimates. On the one hand, this is due to upcoming and already available commercial services (aka Location Based Services) such as intelligent transport systems, fraud detection, yellow page services, location sensitive billing, and other promising services that rely on accurate MT location estimates [1]. On the other hand, the United States Federal Communications Commission (FCC) issued an order, in which all wireless service providers are required to report the location of an E911 caller within a specified accuracy [2]. This FCC mandate together with the emerging LocationBased Services has pushed further the research and standardization activities in the field of MT localization.
Until now, several localization methods have been proposed to solve the problem of locating an MT in a wireless network [3, 4]. The global navigation satellite systems (GNSSs), such as the Global Positioning System (GPS) and the prospective European counterpart Galileo, are promising candidates to fulfill the FCC requirements [5]. In the GNSS, the MT location is estimated from the propagation time; the satellite (SAT) signals need to propagate to the MT, which is known as time of arrival (ToA) method. If the MT receives satellite signals from at least four different satellites, a threedimensional (3D) MT location estimate can be found, where the fourth satellite signal is needed to resolve the unknown bias between the MT and satellite clock [5]. In a similar manner, one can obtain a 2D MT location estimate if the MT receives signals from at least three different satellites. However, there exist situations where the GNSS signals are blocked, for example, when the MT is located indoors or in urban canyons. In these scenarios, the number of satellites in view is often not sufficient to obtain a 3D or even 2D MT location estimate.
An alternative to the GNSS is the exploitation of communication signals of the cellular radio network, in order to obtain MT location estimates. In the Global System for Mobile Communication (GSM), for example, measurements such as the received signal strength (RSS), timing advance (TA), angle of arrival (AoA), or enhanced observed time difference (EOTD) exist that give information on the MT location. An appealing advantage of these measurements is that they are almost everywhere available. However, the corresponding localization methods that are based on these measurements cannot offer the same accuracy as their GNSSbased counterpart. The combination of measured values from the GNSS and the cellular radio network is, thus, a promising approach in order to obtain MT location estimates even if less than four or three satellites are in view [6–13]. The resulting hybrid localization methods are expected to improve the accuracy and availability of MT location estimates.
In [6, 7], a hybrid localization method combining pseudorange (PR) measured values from GPS and EOTDmeasured values from GSM is investigated. In [8], a hybrid method is presented that is based on the fusion of PRmeasured values from GPS and round trip delaymeasured values from a cellular radio network that is, perfectly synchronized to GPS time. However, [6–8] only provide general descriptions of their hybrid methods and no algorithms or theoretical performance bounds are given. In [9], a hybrid method based on the combination of PR measured values from GPS and time difference of arrival (TDoA) measured values from a cellular radio network using a least squares approach is introduced. In [11], a hybrid data fusion approach is presented that combines pseudoranges from the GNSS with TDoA measurements from future 3GPPLTE communication systems using an extended Kalman filter (EKF). In [12, 13], we have developed an extended Kalman filter (EKF) based and RaoBlackwellized unscented Kalman filter (RBUKF) based MT tracking algorithm that fuses TA and RSSmeasured values from GSM and PRmeasured values from GPS.
This paper deals with the combination of RSS, TA, and PRmeasured values from GSM and GPS, as they can be easily obtained from offtheshelf mobile handsets and conventional GPS receivers. The underlying hybrid mobile terminal tracking problem is then solved using the EKF, RBUKF, and the recently proposed cubature Kalman filter (CKF) [14]. Here, a novel extension of the CKF is introduced, accounting for the linear process model structure, which is called the modified cubature Kalman filter (MCKF). Furthermore, it is investigated in this paper how GPS reference time information from the GSM network, which is available from the Radio Resource Location Services Protocol (RRLP) [15], can help to improve the performance of the hybrid localization methods. The different filtering approaches are then compared to each other, the expected computational complexity is evaluated, and their achievable performance is compared in a realistic simulation study with the posterior CramérRao lower bound (PCRLB). The PCRLB gives the theoretical best achievable performance of nonlinear filters [16] and serves here as an important tool for the design of a hybrid MT tracking system. Finally, the three different algorithms are tested on “real world” GSM measurements together with synthetic GPSmeasured values, and their enhanced performance compared to the GSMbased localization method is demonstrated.
The remainder of this paper is organized as follows. In Section 2, the hybrid localization problem is formulated as a nonlinear filtering problem where the optimal solution is given, at least conceptually, by the Bayesian filter. In Section 3, the MT process model and the measurement models for the PR, TA, RSS, and GPS reference time uncertainty are presented that are required to use the different filters. In Section 4, three different nonlinear filters, namely, the EKF, RBUKF, and MCKF are introduced for the hybrid localization problem as well as the PCRLB. The main differences between the different filters are highlighted and the computational complexity is analyzed. In Sections 5 and 6, simulation and experimental results are presented, where the proposed algorithms are compared to each other, and where the advantage of the proposed hybrid method is demonstrated. Finally, Section 7 concludes the work.
2. Problem Statement
In this paper, the MT tracking problem is formulated as a nonlinear filtering problem, where a sequence of measurements available from GSM and GPS is used to estimate the actual state of the MT. Consider the following discretetime statespace model with additive noise: where denotes the discretetime index, denotes the state vector, denotes the measurement vector, and and are some known vectorvalued, possibly nonlinear, mapping functions. Here, it is worth noting that the function models the deterministic relationship between and . Similarly, the function models the deterministic relationship between the state vector and the corresponding measurements available from GPS and GSM. The process and measurement noise and are assumed to be mutually independent zeromean white Gaussian noise sequences with covariances and , respectively.
The aim in nonlinear filtering is to recursively compute estimates of the state using the sequence of all available measurements up to and including time . From a Bayesian point of view, the aim is to recursively compute the posterior probability density function (pdf) , since it provides a complete statistical description of the state at that time. The optimal Bayesian solution is given by the following recursions:
Time Update: Measurement Update: where is a normalizing constant given by and where the pdfs and can be determined from (1). The above recursions are initiated by [17]. It is well known that the nonlinear recursive filtering problem only allows analytical solutions in a few special cases, for example, for linear Gaussian models, where the Kalman filter provides the optimal solution [18]. However, for the general model (1), an analytical solution to the above recursions is intractable and, thus, one has to resort to suboptimal algorithms.
For the hybrid localization method three suboptimal nonlinear filters are investigated, namely, the extended Kalman filter, the RaoBlackwellized unscented Kalman filter, and the modified cubature Kalman filter. But before these filters will be explained in more detail, it will be first shown how the process model and measurement models are chosen for the hybrid localization method.
3. Process and Measurement Model
3.1. Introduction
In the following, it is assumed that the MT location to be estimated and the known base station (BS) locations , , lie in the plane, where denotes the transpose of a vector or matrix. The known satellite locations are given by , . For the case of 3D MT and BS locations, the process and measurement models can be obtained in a similar way. The measurements that are used for the hybrid localization method are the PRmeasured values from GPS and TA, RSS, and GPS reference time uncertainty measured values from GSM. Here, it is worth noting that the hybridization takes place by combining different types of measurements from GPS and GSM rather than location estimates from GPS and GSM. That is, one first collects at every time step all the measurements from GPS and GSM and then these measurements are processed jointly in the filter in order to estimate the MT location. With this strategy, it is possible to obtain MT location estimates even if less than three satellites are visible to the MT.
3.2. Process Model
For the hybrid localization method, the states of the process model include the 2D MT location and velocity, the MT clock bias, and clock drift, that is, , where is the speed of light. The movement of the MT is approximated with a nearly constant velocity (CV) model and the receiver clock bias is modeled by a secondorder GaussMarkov process [19, 20]. The resulting linear process model for the hybrid localization method is, thus, given by with where is the identity matrix of size , denotes the Kronecker product, and is the sampling time. The process noise is assumed to be a zeromean white Gaussian noise sequence with block diagonal covariance matrix . The covariance matrix is given by , where and denote the noise variances in the  and direction. The elements of the symmetric matrix are given by where the parameters and correspond to values of a typical quartz standard [19].
3.3. Measurement Model
3.3.1. Pseudorange
In GPS, the MT is measuring the time the satellite signal requires to travel from the satellite to the MT, which is known as ToA principle [5]. The corresponding ToAmeasured values are affected by delays due to the transmission of the satellite signal through the ionosphere and the troposphere and due to other errors, for example, receiver noise or multipath propagation [5]. In addition to that, the MT's clock is generally not timesynchronized to the clocks of the GPS satellites, resulting in an unknown receiver clock bias that has to be estimated. The corresponding measured biased ranges or measured pseudoranges can be obtained from multiplying the biased ToAmeasured values by .
In the following, it is assumed that each measured pseudorange is corrected for the known errors that are available using parameter values in the navigation message from the satellite [5]. Let denote the vector of corrected PRmeasured values. Then, the PR measurement model can be written as with where denotes the Euclidean distance between the MT and the th satellite. The random variable describes unmodeled effects, modeling errors, and measurement errors; each PRmeasured value is affected by, for example, delays as the signal propagates through the atmosphere, receiver noise, as well as errors due to changing propagation conditions, that is, lineofsight (LOS) or nonlineofsight (NLOS) situations. It is assumed that is Gaussian distributed with mean vector accounting for NLOS propagation and covariance matrix , where and denote the mean and standard deviation from the PRmeasured value of the th satellite.
3.3.2. Timing Advance
In GSM, the Timing Advance (TA) is a parameter that is used to synchronize the transmitted bursts of the MTs to the frame of the receiving BS [1]. In principle, the TA is a quantized value of the round trip time, that is, the time the radio signal requires to propagate from the BS to the MT and back. Let denote the vector of TAmeasured values multiplied by . Then, the TAmeasurement model is given by with where denotes the Euclidean distance between the MT and the th BS. The random variable accounts for the errors each TAmeasured value is affected by, such as quantization, changing propagation conditions—LOS or NLOS situation—and measurement noise. These errors are assumed to be Gaussian distributed with mean vector accounting for NLOS propagation and covariance matrix , where and denote the mean and standard deviation from the TAmeasured value of the th BS.
3.3.3. Received Signal Strength
In GSM, the RSS value is an averaged value of the strength of a radio signal received by the MT. The attenuation of the signal strength through a mobile radio channel is caused by three factors, namely, fast fading, slow fading, and path loss. Since, in GSM, the RSSmeasured values are averaged over several timeconsecutive measurements, the error due to fast fading can be neglected. The model for the path loss in dB is given by Reference[3], where denotes the reference path loss at a BS to MT distance of and is the path loss exponent of the th BS. Both parameters and strongly depend on the propagation conditions and BS antenna settings and can be determined either empirically or from wellknown path loss models as, for example, Hata [21] or COST 231 WalfischIkegami [22].
In real systems, the BSs may be equipped with directional antennas in order to increase the cell's capacity. However, the employment of directional antennas at the BSs should be directly taken into account in the model for the RSS measured value, because otherwise the performance of the tracking algorithms will considerably degrade. In the following, it is assumed that antenna gain models are a priori available. Let and denote the minimum gain and beamwidth of the BS antenna. Let further denote the azimuth angle between the MT and the th BS antenna, counted counterclockwise from the boresight direction of the BS antenna. Then, a model for the normalized antenna gain in dB scale is given by Reference [10], where denotes the smallest value in the set . Let denote the vector of RSSmeasured values. Then, the RSS measurement model in dB scale is given by with where denotes the th BS's equivalent isotropic radiated power. The random variable accounts for errors, such as errors due to slow fading, quantization, and NLOS propagation. It is assumed that is zeromean Gaussian distributed with covariance matrix , where denotes the standard deviation from the RSSmeasured value of the th BS.
3.3.4. GPS Reference Time Uncertainty
In GSM, there exists the possibility to obtain GPS reference time information that can be used to estimate the unknown clock bias in the pseudorange equations (cf. (9)) according to the available RRLP [15]. However, in [15] it is stated that this reference time can be provided only with a specified accuracy which is expressed by the socalled GPS reference time uncertainty. In the following, a model connecting the GPS reference time uncertainty to the unknown MT clock bias will be derived.
Since the satellite clocks can be assumed to be mutually synchronized [5], the MT clock bias can be written as , where the difference describes the offset between the GPS reference time scale , which is unknown to the MT, and the known MT clock timescale . Here, it is worth noting that the bias is not constant over time, since the MT clock experiences errors due to clock drifts. Let denote the GPS reference time uncertainty measurement. Then, the GPS reference time uncertainty measurement model is given by which can be directly converted into an MT clock bias measurement model: where the noise models the GPS reference time uncertainty which is assumed to be zeromean Gaussian distributed with standard deviation . Since we now have related the GPS reference time uncertainty measurement to the clock bias, the uncertainty of the MT clock is implicitly modelled with the process model (cf. (5)) where the MT clock bias evolves according to a secondorder GaussMarkov process.
3.3.5. Combined
In the following, the PR, TA, RSS, and MT clock biasmeasured values are concatenated into a single measurement vector, yielding . Here, GPS reference time uncertaintymeasured values are treated as MT clock biasmeasured values according to (17). The corresponding combined nonlinear measurement model for the hybrid localization problem can be written as where The random variable is Gaussian distributed with mean vector , where denotes the zero vector of size , and block diagonal covariance matrix
4. Nonlinear Filters for Hybrid Localization
4.1. Introduction
After having described the linear process model and the nonlinear relationship between the MT location and the PR, TA, and RSSmeasured values, the problem at hand is how one can efficiently sequentially estimate the MT state from these measured values. The optimal Bayesian solution given by ((2), (3), and (4)) provides a unified approach for nonlinear filtering problems. However, due to the fact that the measurement model is nonlinear (cf. (18)) the multidimensional integral involved in (4) is intractable and, thus, one has to resort to suboptimal algorithms [14, 16, 23–25]. In this paper, three different suboptimal algorithms, namely, the extended Kalman filter, the RaoBlackwellized unscented Kalman filter, and the modified cubature Kalman filter, are proposed in order to solve the underlying hybrid localization problem. These filters belong to the class of approaches where all densities in ((2), (3), and (4)) are assumed to be Gaussian. An appealing advantage of this approximation is that the functional recursion in ((2), (3), and (4)) reduces to an algebraic recursion, where only means and covariances have to be calculated.
4.2. Extended Kalman Filter
In the EKF, the nonlinear functions and are approximated with their firstorder Taylor series expansion, so that an analytical solution of (2) and (4) is possible. This approach, however, leads to several shortcomings. On the one hand, the EKF may have suboptimal performance or even will diverge, if we have a high degree of nonlinearities in the measurement function. On the other hand, the linearization of the measurement model implies the evaluation of Jacobian matrices, which in some cases may become difficult; for example, consider the case when antenna gain models (cf. (13)) are available only from measurements and, consequently, no closed form expressions for these models exist. In these cases, it is much easier to approximate the models using interpolation than trying to evaluate the corresponding Jacobian matrices. The wellknown EKF equations, adopted to the proposed hybrid localization method, are summarized in Algorithm 1 [12].

4.3. RaoBlackwellized Unscented Kalman Filter
While the EKF is based on a simple linear approximation of the nonlinear measurement equation, the unscented Kalman filter (UKF) approximates the multidimensional integrals in (2) and (4) using the (scaled) unscented transformation [26, 27]. In the (scaled) unscented transformation, the multidimensional integrals are approximated using a deterministic sampling procedure. The sampling scheme consists of deterministically choosing a symmetric set of sigma points and weights. These sigma points are then propagated through the true nonlinearity and the corresponding mean and covariances are approximated using a weighted sample mean and covariance. Compared to the EKF, the advantage of the UKF is that no Jacobian matrices have to be evaluated, since the sigma points are transformed through the true nonlinearity. Furthermore, it can be shown that the nonlinear transformed samples capture the mean and covariance accurately to at least the secondorder of the Taylor series expansion whereas the EKF only achieves firstorder accuracy [27]. Since the noise in the process model and measurement model is assumed to be additive and Gaussian distributed, the dimension of the vector, from which the sigma points are sampled, can be reduced. This technique is also known as RaoBlackwellization and has the advantage that the quasiMonte Carlo variance and computational complexity can be reduced [28]. The computational complexity can be further reduced by taking into account that the posterior pdf is assumed Gaussian and the process model is linear Gaussian. In this case, the multidimensional integral, given in (2), can be evaluated in closed form, resulting in the wellknown Kalman filter update equations. The corresponding RBUKF, adapted to the hybrid localization method, is summarized in Algorithm 2 [13].

4.4. Modified Cubature Kalman Filter
In the recently proposed CKF [14], the multidimensional integrals in (2) and (4) are approximated in a different way. Since the conditional pdfs and in (2) and (4) are assumed Gaussian, solving approximately the multidimensional integrals is equivalent to the evaluation of the corresponding means and covariances of and. It can be shown that the evaluation of the mean and covariance leads again to multidimensional integrals, but whose integrands are now all of the form [14]. These integrals are then solved using highly efficient numerical integration methods known as cubature rules. As a result, one obtains a set of cubature points and weights, from which the corresponding mean and covariance can be computed without evaluating Jacobian matrices. Due to the fact that (2) can be evaluated in closed form, the computational complexity of the CKF can be further reduced. The modified version of the CKF that is used for the hybrid localization problem is summarized in Algorithm 3. Although the CKF and UKF seem to be very similar, the authors in [14] claim that the cubature approach is more accurate and more principled in mathematical terms than the sigmapoint approach used in the UKF. Here, it is worth noting that the RBUKF reduces to the MCKF, when the unknown parameters in the scaled unscented transformation are chosen as , and in Algorithm 2.

4.5. Posterior CramérRao Lower Bound
After having introduced the different filters for the hybrid localization problem, their performance should be compared to a theoretical bound. In the following, the posterior CramérRao lower bound for the hybrid localization problem is presented that gives the best achievable performance for nonlinear filtering [16, 29]. Let be an unbiased estimate of the state vector . Then, the covariance matrix of the estimation error satisfies the inequality where is the expectation with respect to , denotes the filtering information matrix, and its inverse is the PCRLB matrix. The matrix inequality should be interpreted as the matrix being positive semidefinite. The aim is now to calculate . In [29], an elegant method is presented, where can be determined recursively. This recursion, adapted to the hybrid localization problem involving additive Gaussian noise (cf. (5) and (18)) can be written as where the expectation is with respect to and denotes the Jacobian matrix of the nonlinear measurement function (cf. (19)) evaluated at the true value of the state . Since the initial distribution is assumed to be Gaussian, the recursions are initialized with the information matrix [16].
4.6. Computational Complexity
In this section, the computational complexity of the EKF, RBUKF, and MCKF for the hybrid localization method is investigated in terms of floatingpoint operations (FLOPs). A FLOP is here defined as one addition, subtraction, multiplication, or division of two floatingpoint numbers. In Table 1, the computational complexity of some common matrix operations is summarized. Here, it is worth noting that the matrix square root, which is needed to evaluate the set of cubature and sigma points, is computed using Cholesky decomposition, whose complexity grows cubically.

In the EKF as well as in the RBUKF and MCKF, there are certain steps that cannot be measured in FLOPs. In the EKF, for example, one has to evaluate at every time step the Jacobian matrix and the nonlinear function (cf. Algorithm 1). In the RBUKF and MCKF, one has to propagate at every time step sigma points and cubature points through the nonlinear function (cf. Algorithms 2 and 3). In the following, the cost of evaluating a certain nonlinear function and Jacobian matrix is neglected. Furthermore, the computation of the weights in the RBUKF, and MCKF as well as the initialization of all three filters can be neglected, since these steps are done only once.
In Table 2, the computational complexity of the different quantities that have to be evaluated in the EKF, RBUKF, and MCKF is presented. Summing up the computational complexity of the different quantities results in the total FLOP complexity of the EKF, RBUKF, and MCKF for one time step which is given by where and denote the dimension of the state and measurement vector, respectively.
