# Nonlinear Filtering for Hybrid GPS/GSM Mobile Terminal Tracking

Authors: Carsten Fritsche; Anja Klein
Journal: International Journal of Navigation and Observation
 (2010)

---

## Abstract

The Global Positioning System (GPS) has become one of the state-of-the-art 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 GPS-measured 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 Rao-Blackwellized 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ér-Rao 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.


---

## Body

## 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 E-911 caller within a specified accuracy [2]. This FCC mandate together with the emerging Location-Based 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 three-dimensional (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 (E-OTD) 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 GNSS-based 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 E-OTD-measured values from GSM is investigated. In [8], a hybrid method is presented that is based on the fusion of PR-measured values from GPS and round trip delay-measured 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 3GPP-LTE communication systems using an extended Kalman filter (EKF). In [12, 13], we have developed an extended Kalman filter- (EKF-) based and Rao-Blackwellized unscented Kalman filter- (RBUKF-) based MT tracking algorithm that fuses TA- and RSS-measured values from GSM- and PR-measured values from GPS.This paper deals with the combination of RSS-, TA-, and PR-measured values from GSM and GPS, as they can be easily obtained from off-the-shelf 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ér-Rao 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 GPS-measured values, and their enhanced performance compared to the GSM-based localization method is demonstrated.The remainder of this paper is organized as follows. In Section2, 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 discrete-time state-space model with additive noise:(1)ProcessModel:x(k)=f(x(k-1))+v(k-1),MeasurementModel:y(k)=h(x(k))+w(k),
where k denotes the discrete-time index, x(k)∈ℝnx denotes the state vector, y(k)∈ℝny denotes the measurement vector, and f(·) and h(·) are some known vector-valued, possibly nonlinear, mapping functions. Here, it is worth noting that the function f(·) models the deterministic relationship between x(k) and x(k-1). Similarly, the function h(·) models the deterministic relationship between the state vector x(k) and the corresponding measurements y(k) available from GPS and GSM. The process and measurement noise v(k-1) and w(k) are assumed to be mutually independent zero-mean white Gaussian noise sequences with covariances Q(k-1) and R(k), respectively.The aim in nonlinear filtering is to recursively compute estimates of the statex(k) using the sequence of all available measurements Y(k)={y(l),l=1,…,k} up to and including time k. From a Bayesian point of view, the aim is to recursively compute the posterior probability density function (pdf) p(x(k)∣Y(k)), since it provides a complete statistical description of the state x(k) at that time. The optimal Bayesian solution is given by the following recursions:Time Update:(2)p(x(k)∣Y(k-1))=∫ℝnxp(x(k)∣x(k-1))×p(x(k-1)∣Y(k-1))dx(k-1).
Measurement Update:(3)p(x(k)∣Y(k))=p(y(k)∣x(k))p(x(k)∣Y(k-1))p(y(k)∣Y(k-1)),
where p(y(k)∣Y(k-1)) is a normalizing constant given by(4)p(y(k)∣Y(k-1))=∫ℝnxp(y(k)∣x(k))p(x(k)∣Y(k-1))dx(k),
and where the pdfs p(x(k)∣x(k-1)) and p(y(k)∣x(k)) can be determined from (1). The above recursions are initiated by p(x(0)∣Y(0))=p(x(0)) [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 Rao-Blackwellized 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 locationxMT=[xMT,yMT]⊤ to be estimated and the known base station (BS) locations xBS(n)=[xBS(n),yBS(n)]⊤, n=1,…,NBS, lie in the xy-plane, where [·]⊤ denotes the transpose of a vector or matrix. The known satellite locations are given by xSAT(l)=[xSAT(l),ySAT(l),zSAT(l)]⊤, l=1,…,NSAT. For the case of 3-D 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 PR-measured 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 k 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 2-D MT location and velocity, the MT clock bias, and clock drift, that is,x=[xMT,ẋMT,yMT,ẏMT,c0·δt,c0·δṫ]⊤, where c0 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 second-order Gauss-Markov process [19, 20]. The resulting linear process model for the hybrid localization method is, thus, given by(5)x(k)=Φ·x(k-1)+Γ·v(k-1)
with(6)Φ=I3⊗[1TS01],Γ=diag(I2⊗[TS22,TS]⊤,I2·c0),
where Iq is the identity matrix of size q, ⊗ denotes the Kronecker product, and TS is the sampling time. The process noise v=[vx,vy,vδt,vδṫ]⊤ is assumed to be a zero-mean white Gaussian noise sequence with block diagonal covariance matrix Q=diag(QCV,Qδt). The covariance matrix QCV is given by QCV=diag(σx2,σy2), where σx2 and σy2 denote the noise variances in the x- and y-direction. The elements of the symmetric 2×2 matrix Qδt are given by(7)Q11=h0TS2+2h-1TS2+23π2h-2TS3,Q12=Q21=2h-1TS+π2h-2TS2,Q22=h02TS+2h-1+83π2h-2TS,
where the parameters h0,h-1, and h-2 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 ToA-measured 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 time-synchronized to the clocks of the GPS satellites, resulting in an unknown receiver clock bias δt(k) that has to be estimated. The corresponding measured biased ranges or measured pseudoranges can be obtained from multiplying the biased ToA-measured values by c0.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 yPR(k) denote the vector of NSAT-corrected PR-measured values. Then, the PR measurement model can be written as(8)yPR(k)=hPR(xMT(k),δt(k))+wPR(k)
with(9)hPR(xMT(k),δt(k))=[dSAT(1)(xMT(k)),…,dSAT(NSAT)(xMT(k))]⊤+c0·δt(k),
where dSAT(l)(xMT(k)) denotes the Euclidean distance between the MT and the lth satellite. The random variable wPR(k) describes unmodeled effects, modeling errors, and measurement errors; each PR-measured 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, line-of-sight (LOS) or non-line-of-sight (NLOS) situations. It is assumed that wPR(k) is Gaussian distributed with mean vector μPR=[μPR(1),…,μPR(NSAT)]⊤ accounting for NLOS propagation and covariance matrix RPR=diag((σPR(1))2,…,(σPR(NSAT))2), where μPR(l) and σPR(l) denote the mean and standard deviation from the PR-measured value of the lth 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 yTA(k) denote the vector of NBS TA-measured values multiplied by c0/2. Then, the TA-measurement model is given by(10)yTA(k)=hTA(xMT(k))+wTA(k)
with(11)hTA(xMT(k))=[dBS(1)(xMT(k)),…,dBS(NBS)(xMT(k))]⊤,
where dBS(n)(xMT(k)) denotes the Euclidean distance between the MT and the nth BS. The random variable wTA(k) accounts for the errors each TA-measured 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 μTA=[μTA(1),…,μTA(NBS)]⊤ accounting for NLOS propagation and covariance matrix RTA=diag((σTA(1))2,…,(σTA(NBS))2), where μTA(n) and σTA(n) denote the mean and standard deviation from the TA-measured value of the nth 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 RSS-measured values are averaged over several time-consecutive measurements, the error due to fast fading can be neglected. The model for the path loss in dB is given by(12)L(n)(xMT(k))=A(n)+10·B(n)·log10(dBS(n)(xMT(k))1km).
Reference[3], where A(n) denotes the reference path loss at a BS to MT distance of 1km and B(n) is the path loss exponent of the nth BS. Both parameters A(n) and B(n) strongly depend on the propagation conditions and BS antenna settings and can be determined either empirically or from well-known path loss models as, for example, Hata [21] or COST 231 Walfisch-Ikegami [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. LetAm(n) and φ3dB(n) denote the minimum gain and 3dB beamwidth of the BS antenna. Let further φBS(n)(xMT(k)) denote the azimuth angle between the MT and the nth 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(13)g(φBS(n)(xMT(k)))=-min{12(φBS(n)(xMT(k))φ3dB(n))2,Am(n)}.
Reference [10], where min{a,b} denotes the smallest value in the set {a,b}. Let yRSS(k) denote the vector of NBS RSS-measured values. Then, the RSS measurement model in dB scale is given by(14)yRSS(k)=hRSS(xMT(k))+wRSS(k)
with(15)hRSS(xMT(k))=[hRSS(1)(xMT(k)),…,hRSS(NBS)(xMT(k))]⊤,hRSS(n)(xMT(k))=PT(n)-{L(n)(xMT(k))-g(φBS(n)(xMT(k)))},
where PT(n) denotes the nth BS's equivalent isotropic radiated power. The random variable wRSS(k) accounts for errors, such as errors due to slow fading, quantization, and NLOS propagation. It is assumed that wRSS(k) is zero-mean Gaussian distributed with covariance matrix RRSS=diag((σRSS(1))2,…,(σRSS(NBS))2), where σRSS(n) denotes the standard deviation from the RSS-measured value of the nth 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δt(k) 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 so-called 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 δt(k)=tGPS(k)-tMTC(k), where the difference describes the offset between the GPS reference time scale tGPS(k), which is unknown to the MT, and the known MT clock time-scale tMTC(k). Here, it is worth noting that the bias is not constant over time, since the MT clock experiences errors due to clock drifts. Let yRTU(k) denote the GPS reference time uncertainty measurement. Then, the GPS reference time uncertainty measurement model is given by(16)yRTU(k)=tGPS(k)+wRTU(k),
which can be directly converted into an MT clock bias measurement model:(17)yBIAS(k)=yRTU(k)-tMTC(k)=δt(k)+wRTU(k),
where the noise wRTU(k) models the GPS reference time uncertainty which is assumed to be zero-mean Gaussian distributed with standard deviation σRTU. 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 second-order Gauss-Markov process.


#### 3.3.5. Combined

In the following, the PR, TA, RSS, and MT clock bias-measured values are concatenated into a single measurement vector, yieldingy(k)=[yPR⊤(k),yTA⊤(k),yRSS⊤(k),yBIAS(k)]⊤. Here, GPS reference time uncertainty-measured values are treated as MT clock bias-measured values according to (17). The corresponding combined nonlinear measurement model for the hybrid localization problem can be written as(18)y(k)=h(x(k))+w(k),
where(19)h(x(k))=[hPR⊤(xMT(k),δt(k)),hTA⊤(xMT(k)),hRSS⊤(xMT(k)),δt(k)]⊤,(20)w(k)=[wPR⊤(k),wTA⊤(k),wRSS⊤(k),wRTU(k)]⊤.
The random variable w(k) is Gaussian distributed with mean vector μ=[μPR⊤(k),μTA⊤(k),0]⊤, where 0 denotes the zero vector of size 1×NBS+1, and block diagonal covariance matrix R=diag(RPR,RTA,RRSS,σRTU).


## 3.1. Introduction

In the following, it is assumed that the MT locationxMT=[xMT,yMT]⊤ to be estimated and the known base station (BS) locations xBS(n)=[xBS(n),yBS(n)]⊤, n=1,…,NBS, lie in the xy-plane, where [·]⊤ denotes the transpose of a vector or matrix. The known satellite locations are given by xSAT(l)=[xSAT(l),ySAT(l),zSAT(l)]⊤, l=1,…,NSAT. For the case of 3-D 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 PR-measured 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 k 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 2-D MT location and velocity, the MT clock bias, and clock drift, that is,x=[xMT,ẋMT,yMT,ẏMT,c0·δt,c0·δṫ]⊤, where c0 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 second-order Gauss-Markov process [19, 20]. The resulting linear process model for the hybrid localization method is, thus, given by(5)x(k)=Φ·x(k-1)+Γ·v(k-1)
with(6)Φ=I3⊗[1TS01],Γ=diag(I2⊗[TS22,TS]⊤,I2·c0),
where Iq is the identity matrix of size q, ⊗ denotes the Kronecker product, and TS is the sampling time. The process noise v=[vx,vy,vδt,vδṫ]⊤ is assumed to be a zero-mean white Gaussian noise sequence with block diagonal covariance matrix Q=diag(QCV,Qδt). The covariance matrix QCV is given by QCV=diag(σx2,σy2), where σx2 and σy2 denote the noise variances in the x- and y-direction. The elements of the symmetric 2×2 matrix Qδt are given by(7)Q11=h0TS2+2h-1TS2+23π2h-2TS3,Q12=Q21=2h-1TS+π2h-2TS2,Q22=h02TS+2h-1+83π2h-2TS,
where the parameters h0,h-1, and h-2 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 ToA-measured 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 time-synchronized to the clocks of the GPS satellites, resulting in an unknown receiver clock bias δt(k) that has to be estimated. The corresponding measured biased ranges or measured pseudoranges can be obtained from multiplying the biased ToA-measured values by c0.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 yPR(k) denote the vector of NSAT-corrected PR-measured values. Then, the PR measurement model can be written as(8)yPR(k)=hPR(xMT(k),δt(k))+wPR(k)
with(9)hPR(xMT(k),δt(k))=[dSAT(1)(xMT(k)),…,dSAT(NSAT)(xMT(k))]⊤+c0·δt(k),
where dSAT(l)(xMT(k)) denotes the Euclidean distance between the MT and the lth satellite. The random variable wPR(k) describes unmodeled effects, modeling errors, and measurement errors; each PR-measured 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, line-of-sight (LOS) or non-line-of-sight (NLOS) situations. It is assumed that wPR(k) is Gaussian distributed with mean vector μPR=[μPR(1),…,μPR(NSAT)]⊤ accounting for NLOS propagation and covariance matrix RPR=diag((σPR(1))2,…,(σPR(NSAT))2), where μPR(l) and σPR(l) denote the mean and standard deviation from the PR-measured value of the lth 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 yTA(k) denote the vector of NBS TA-measured values multiplied by c0/2. Then, the TA-measurement model is given by(10)yTA(k)=hTA(xMT(k))+wTA(k)
with(11)hTA(xMT(k))=[dBS(1)(xMT(k)),…,dBS(NBS)(xMT(k))]⊤,
where dBS(n)(xMT(k)) denotes the Euclidean distance between the MT and the nth BS. The random variable wTA(k) accounts for the errors each TA-measured 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 μTA=[μTA(1),…,μTA(NBS)]⊤ accounting for NLOS propagation and covariance matrix RTA=diag((σTA(1))2,…,(σTA(NBS))2), where μTA(n) and σTA(n) denote the mean and standard deviation from the TA-measured value of the nth 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 RSS-measured values are averaged over several time-consecutive measurements, the error due to fast fading can be neglected. The model for the path loss in dB is given by(12)L(n)(xMT(k))=A(n)+10·B(n)·log10(dBS(n)(xMT(k))1km).
Reference[3], where A(n) denotes the reference path loss at a BS to MT distance of 1km and B(n) is the path loss exponent of the nth BS. Both parameters A(n) and B(n) strongly depend on the propagation conditions and BS antenna settings and can be determined either empirically or from well-known path loss models as, for example, Hata [21] or COST 231 Walfisch-Ikegami [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. LetAm(n) and φ3dB(n) denote the minimum gain and 3dB beamwidth of the BS antenna. Let further φBS(n)(xMT(k)) denote the azimuth angle between the MT and the nth 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(13)g(φBS(n)(xMT(k)))=-min{12(φBS(n)(xMT(k))φ3dB(n))2,Am(n)}.
Reference [10], where min{a,b} denotes the smallest value in the set {a,b}. Let yRSS(k) denote the vector of NBS RSS-measured values. Then, the RSS measurement model in dB scale is given by(14)yRSS(k)=hRSS(xMT(k))+wRSS(k)
with(15)hRSS(xMT(k))=[hRSS(1)(xMT(k)),…,hRSS(NBS)(xMT(k))]⊤,hRSS(n)(xMT(k))=PT(n)-{L(n)(xMT(k))-g(φBS(n)(xMT(k)))},
where PT(n) denotes the nth BS's equivalent isotropic radiated power. The random variable wRSS(k) accounts for errors, such as errors due to slow fading, quantization, and NLOS propagation. It is assumed that wRSS(k) is zero-mean Gaussian distributed with covariance matrix RRSS=diag((σRSS(1))2,…,(σRSS(NBS))2), where σRSS(n) denotes the standard deviation from the RSS-measured value of the nth 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δt(k) 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 so-called 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 δt(k)=tGPS(k)-tMTC(k), where the difference describes the offset between the GPS reference time scale tGPS(k), which is unknown to the MT, and the known MT clock time-scale tMTC(k). Here, it is worth noting that the bias is not constant over time, since the MT clock experiences errors due to clock drifts. Let yRTU(k) denote the GPS reference time uncertainty measurement. Then, the GPS reference time uncertainty measurement model is given by(16)yRTU(k)=tGPS(k)+wRTU(k),
which can be directly converted into an MT clock bias measurement model:(17)yBIAS(k)=yRTU(k)-tMTC(k)=δt(k)+wRTU(k),
where the noise wRTU(k) models the GPS reference time uncertainty which is assumed to be zero-mean Gaussian distributed with standard deviation σRTU. 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 second-order Gauss-Markov process.


### 3.3.5. Combined

In the following, the PR, TA, RSS, and MT clock bias-measured values are concatenated into a single measurement vector, yieldingy(k)=[yPR⊤(k),yTA⊤(k),yRSS⊤(k),yBIAS(k)]⊤. Here, GPS reference time uncertainty-measured values are treated as MT clock bias-measured values according to (17). The corresponding combined nonlinear measurement model for the hybrid localization problem can be written as(18)y(k)=h(x(k))+w(k),
where(19)h(x(k))=[hPR⊤(xMT(k),δt(k)),hTA⊤(xMT(k)),hRSS⊤(xMT(k)),δt(k)]⊤,(20)w(k)=[wPR⊤(k),wTA⊤(k),wRSS⊤(k),wRTU(k)]⊤.
The random variable w(k) is Gaussian distributed with mean vector μ=[μPR⊤(k),μTA⊤(k),0]⊤, where 0 denotes the zero vector of size 1×NBS+1, and block diagonal covariance matrix R=diag(RPR,RTA,RRSS,σRTU).


## 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 ToA-measured 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 time-synchronized to the clocks of the GPS satellites, resulting in an unknown receiver clock bias δt(k) that has to be estimated. The corresponding measured biased ranges or measured pseudoranges can be obtained from multiplying the biased ToA-measured values by c0.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 yPR(k) denote the vector of NSAT-corrected PR-measured values. Then, the PR measurement model can be written as(8)yPR(k)=hPR(xMT(k),δt(k))+wPR(k)
with(9)hPR(xMT(k),δt(k))=[dSAT(1)(xMT(k)),…,dSAT(NSAT)(xMT(k))]⊤+c0·δt(k),
where dSAT(l)(xMT(k)) denotes the Euclidean distance between the MT and the lth satellite. The random variable wPR(k) describes unmodeled effects, modeling errors, and measurement errors; each PR-measured 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, line-of-sight (LOS) or non-line-of-sight (NLOS) situations. It is assumed that wPR(k) is Gaussian distributed with mean vector μPR=[μPR(1),…,μPR(NSAT)]⊤ accounting for NLOS propagation and covariance matrix RPR=diag((σPR(1))2,…,(σPR(NSAT))2), where μPR(l) and σPR(l) denote the mean and standard deviation from the PR-measured value of the lth 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 yTA(k) denote the vector of NBS TA-measured values multiplied by c0/2. Then, the TA-measurement model is given by(10)yTA(k)=hTA(xMT(k))+wTA(k)
with(11)hTA(xMT(k))=[dBS(1)(xMT(k)),…,dBS(NBS)(xMT(k))]⊤,
where dBS(n)(xMT(k)) denotes the Euclidean distance between the MT and the nth BS. The random variable wTA(k) accounts for the errors each TA-measured 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 μTA=[μTA(1),…,μTA(NBS)]⊤ accounting for NLOS propagation and covariance matrix RTA=diag((σTA(1))2,…,(σTA(NBS))2), where μTA(n) and σTA(n) denote the mean and standard deviation from the TA-measured value of the nth 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 RSS-measured values are averaged over several time-consecutive measurements, the error due to fast fading can be neglected. The model for the path loss in dB is given by(12)L(n)(xMT(k))=A(n)+10·B(n)·log10(dBS(n)(xMT(k))1km).
Reference[3], where A(n) denotes the reference path loss at a BS to MT distance of 1km and B(n) is the path loss exponent of the nth BS. Both parameters A(n) and B(n) strongly depend on the propagation conditions and BS antenna settings and can be determined either empirically or from well-known path loss models as, for example, Hata [21] or COST 231 Walfisch-Ikegami [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. LetAm(n) and φ3dB(n) denote the minimum gain and 3dB beamwidth of the BS antenna. Let further φBS(n)(xMT(k)) denote the azimuth angle between the MT and the nth 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(13)g(φBS(n)(xMT(k)))=-min{12(φBS(n)(xMT(k))φ3dB(n))2,Am(n)}.
Reference [10], where min{a,b} denotes the smallest value in the set {a,b}. Let yRSS(k) denote the vector of NBS RSS-measured values. Then, the RSS measurement model in dB scale is given by(14)yRSS(k)=hRSS(xMT(k))+wRSS(k)
with(15)hRSS(xMT(k))=[hRSS(1)(xMT(k)),…,hRSS(NBS)(xMT(k))]⊤,hRSS(n)(xMT(k))=PT(n)-{L(n)(xMT(k))-g(φBS(n)(xMT(k)))},
where PT(n) denotes the nth BS's equivalent isotropic radiated power. The random variable wRSS(k) accounts for errors, such as errors due to slow fading, quantization, and NLOS propagation. It is assumed that wRSS(k) is zero-mean Gaussian distributed with covariance matrix RRSS=diag((σRSS(1))2,…,(σRSS(NBS))2), where σRSS(n) denotes the standard deviation from the RSS-measured value of the nth 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δt(k) 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 so-called 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 δt(k)=tGPS(k)-tMTC(k), where the difference describes the offset between the GPS reference time scale tGPS(k), which is unknown to the MT, and the known MT clock time-scale tMTC(k). Here, it is worth noting that the bias is not constant over time, since the MT clock experiences errors due to clock drifts. Let yRTU(k) denote the GPS reference time uncertainty measurement. Then, the GPS reference time uncertainty measurement model is given by(16)yRTU(k)=tGPS(k)+wRTU(k),
which can be directly converted into an MT clock bias measurement model:(17)yBIAS(k)=yRTU(k)-tMTC(k)=δt(k)+wRTU(k),
where the noise wRTU(k) models the GPS reference time uncertainty which is assumed to be zero-mean Gaussian distributed with standard deviation σRTU. 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 second-order Gauss-Markov process.


## 3.3.5. Combined

In the following, the PR, TA, RSS, and MT clock bias-measured values are concatenated into a single measurement vector, yieldingy(k)=[yPR⊤(k),yTA⊤(k),yRSS⊤(k),yBIAS(k)]⊤. Here, GPS reference time uncertainty-measured values are treated as MT clock bias-measured values according to (17). The corresponding combined nonlinear measurement model for the hybrid localization problem can be written as(18)y(k)=h(x(k))+w(k),
where(19)h(x(k))=[hPR⊤(xMT(k),δt(k)),hTA⊤(xMT(k)),hRSS⊤(xMT(k)),δt(k)]⊤,(20)w(k)=[wPR⊤(k),wTA⊤(k),wRSS⊤(k),wRTU(k)]⊤.
The random variable w(k) is Gaussian distributed with mean vector μ=[μPR⊤(k),μTA⊤(k),0]⊤, where 0 denotes the zero vector of size 1×NBS+1, and block diagonal covariance matrix R=diag(RPR,RTA,RRSS,σRTU).


## 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 RSS-measured 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 Rao-Blackwellized 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 functionsf(x(k-1)) and h(x(k)) are approximated with their first-order 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 well-known EKF equations, adopted to the proposed hybrid localization method, are summarized in Algorithm 1 [12].Algorithm 1:Extended Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1),P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤.Measurement UpdateEstimate the innovation covariance and cross-covariance matrixPyy(k∣k-1)=H(k)P(k∣k-1)HT(k)+RPxy(k∣k-1)=P(k∣k-1)HT(k).Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-h(x̂(k∣k-1))-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k),whereH(k) denotes the Jacobian matrix of h(x(k)) evaluated at x̂(k∣k-1).


### 4.3. Rao-Blackwellized 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 second-order of the Taylor series expansion whereas the EKF only achieves first-order 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 Rao-Blackwellization and has the advantage that the quasi-Monte Carlo variance and computational complexity can be reduced [28]. The computational complexity can be further reduced by taking into account that the posterior pdf p(x(k-1)∣Y(k-1)) 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 well-known Kalman filter update equations. The corresponding RBUKF, adapted to the hybrid localization method, is summarized in Algorithm 2 [13].Algorithm 2:Rao-Blackwellized Unscented Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1)P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤Measurement UpdateEvaluate the sigma points (i=1,…,nx)𝒳0(k∣k-1)=x̂(k∣k-1),𝒳i(k∣k-1)=x̂(k∣k-1)+((nx+λ)P(k∣k-1))i,𝒳nx+i(k∣k-1)=x̂(k∣k-1)-((nx+λ)P(k∣k-1))i.Evaluate the weights (i=1,…,2nx)W0(m)=λnx+λ,W0(c)=λnx+λ+(1-α2+β),Wi(m)=Wi(c)=12(nx+λ).Evaluate the propagated sigma points (i=0,…,2nx)𝒴i(k∣k-1)=h(𝒳i(k∣k-1)).Estimate the predicted measurementŷ(k∣k-1)=∑i=02nxWi(m)𝒴i(k∣k-1).Estimate the innovation covariance and cross-covariance matrix (i=0,…,2nx)𝒳̃i(k∣k-1)=𝒳i(k∣k-1)-x̂(k∣k-1),𝒴̃i(k∣k-1)=𝒴i(k∣k-1)-ŷ(k∣k-1),Pyy(k∣k-1)=∑i=02nxWi(c)𝒴̃i(k∣k-1)𝒴̃i⊤(k∣k-1)+R,Pxy(k∣k-1)=∑i=02nxWi(c)𝒳̃i(k∣k-1)𝒴̃i⊤(k∣k-1).Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-ŷ(k∣k-1)-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k),wherenx is the dimension of state vector, λ=α2(nx+κ)-nx is a scalingparameter,α determines the spread of the sigma points, κ is a secondary scalingparameter andβ is a weight parameter. (A)i denotes the ith column of thematrixA


### 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 p(x(k)∣Y(k-1)) and p(y(k)∣Y(k-1)) in (2) and (4) are assumed Gaussian, solving approximately the multidimensional integrals is equivalent to the evaluation of the corresponding means and covariances of p(x(k)∣Y(k-1)) andp(y(k)∣Y(k-1)). 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 nonlinearfunction×Gaussiandensity [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 sigma-point 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 α=1,β=0, and κ=0 in Algorithm 2.Algorithm 3:Modified Cubature Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1),P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤.Measurement UpdateEvaluate the cubature points (i=1,…,2nx)P(k∣k-1)=S(k∣k-1)S⊤(k∣k-1),Ξ=nx·[Inx,-Inx],𝒳i(k∣k-1)=x̂(k∣k-1)+S(k∣k-1)ξi.Evaluate the weights (i=1,…,2nx)Wi(m)=Wi(c)=1/(2nx).Evaluate the propagated cubature points (i=1,…,2nx)𝒴i(k∣k-1)=h(𝒳i(k∣k-1)).Estimate the predicted measurementŷ(k∣k-1)∑i=12nxWi(m)𝒴i(k∣k-1)Estimate the innovation covariance and cross-covariance matrixPyy(k∣k-1)=∑i=12nxWi(c)𝒴i(k∣k-1)𝒴i⊤(k∣k-1)-ŷ(k∣k-1)ŷ⊤(k∣k-1)+R,Pxy(k∣k-1)=∑i=12nxWi(c)𝒳i(k∣k-1)𝒴i⊤(k∣k-1)-x̂(k∣k-1)ŷ⊤(k∣k-1)Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-ŷ(k∣k-1)-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k).whereξi denotes the ith column vector of the matrix Ξ.


### 4.5. Posterior Cramér-Rao 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ér-Rao lower bound for the hybrid localization problem is presented that gives the best achievable performance for nonlinear filtering [16, 29]. Let x̂(k∣k) be an unbiased estimate of the state vector x(k). Then, the covariance matrix of the estimation error satisfies the inequality(21)𝔼{(x̂(k∣k)-x(k))(x̂(k∣k)-x(k))⊤}≥J-1(k),
where 𝔼{·} is the expectation with respect to x(k), J(k) denotes the filtering information matrix, and its inverse is the PCRLB matrix. The matrix inequality A≥B should be interpreted as the matrix A-B being positive semidefinite. The aim is now to calculate J(k). In [29], an elegant method is presented, where J(k) can be determined recursively. This recursion, adapted to the hybrid localization problem involving additive Gaussian noise (cf. (5) and (18)) can be written as(22)J(k)=(ΓQΓ⊤+ΦJ-1(k-1)Φ⊤)-1+𝔼{H̃⊤(k)R-1H̃(k)},
where the expectation is with respect to x(k) and H̃(k) denotes the Jacobian matrix of the nonlinear measurement function h(·) (cf. (19)) evaluated at the true value of the state x(k). Since the initial distribution p(x(0)) is assumed to be Gaussian, the recursions are initialized with the information matrix J(0)=P-1(0∣0) [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 floating-point operations (FLOPs). A FLOP is here defined as one addition, subtraction, multiplication, or division of two floating-point numbers. In Table1, 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.Table 1
Computational complexity of some common matrix operations [30].
OperationSizeMult.Add.OtherA+AA∈ℝn×m—nm—A·BA∈ℝn×m,B∈ℝm×llmn(m-1)ln—C-1C∈ℝn×nn3——CC∈ℝn×n——n3/3+2n2In 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 stepk the Jacobian matrix H(k) and the nonlinear function h(·) (cf. Algorithm 1). In the RBUKF and MCKF, one has to propagate at every time step 2nx+1 sigma points and 2nx cubature points through the nonlinear function h(·) (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 Table2, 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(23)𝒞EKF(nx,ny)=8nx3+ny3+6nx2ny+6ny2nx-13nx2-2nxny+9nx+2ny,𝒞RBUKF(nx,ny)=373nx3+ny3+6nx2ny+8ny2nx-9nx2+7nxny+4ny2+10nx+5ny,𝒞MCKF(nx,ny)=373nx3+ny3+6nx2ny+8ny2nx-11nx2+4nxny+3ny2+9nx+2ny,
where nx and ny denote the dimension of the state and measurement vector, respectively.Table 2
Computational complexity of the EKF, RBUKF, and MCKF. The dimension of the state vector is given bynx and the dimension of the measurement vector is given by ny. 𝒳(k∣k-1) denotes the matrix composed of sigma/cubature point vectors.
QuantityComplexityEKFRBUKFx̂(k∣k-1)2nx2-nx2nx2-nx2nx2-nxP(k∣k-1)8nx3-15nx2+10nx8nx3-15nx2+10nx8nx3-15nx2+10nxX(k∣k-1)—13nx3/3+2nx213nx3/3+2nx2ŷ(k∣k-1)—2nxny+2ny2nxnyPyy(k∣k-1)2ny2nx+2nx2ny-nxny4ny2nx+2nxny+4ny2+ny4ny2nx+3ny2Pxy(k∣k-1)2nx2ny-nxny4nx2ny+3nxny+2nx2+nx4nx2ny+2nxnyK(k)ny3+2ny2nx-nxnyny3+2ny2nx-nxnyny3+2ny2nx-nxnyx̂(k∣k)2nxny+2ny2nxny+2ny2nxny+2nyP(k∣k)2ny2nx+2nx2ny-nxny2ny2nx+2nx2ny-nxny2ny2nx+2nx2ny-nxny


## 4.1. Introduction

After having described the linear process model and the nonlinear relationship between the MT location and the PR-, TA-, and RSS-measured 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 Rao-Blackwellized 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 functionsf(x(k-1)) and h(x(k)) are approximated with their first-order 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 well-known EKF equations, adopted to the proposed hybrid localization method, are summarized in Algorithm 1 [12].Algorithm 1:Extended Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1),P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤.Measurement UpdateEstimate the innovation covariance and cross-covariance matrixPyy(k∣k-1)=H(k)P(k∣k-1)HT(k)+RPxy(k∣k-1)=P(k∣k-1)HT(k).Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-h(x̂(k∣k-1))-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k),whereH(k) denotes the Jacobian matrix of h(x(k)) evaluated at x̂(k∣k-1).


## 4.3. Rao-Blackwellized 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 second-order of the Taylor series expansion whereas the EKF only achieves first-order 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 Rao-Blackwellization and has the advantage that the quasi-Monte Carlo variance and computational complexity can be reduced [28]. The computational complexity can be further reduced by taking into account that the posterior pdf p(x(k-1)∣Y(k-1)) 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 well-known Kalman filter update equations. The corresponding RBUKF, adapted to the hybrid localization method, is summarized in Algorithm 2 [13].Algorithm 2:Rao-Blackwellized Unscented Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1)P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤Measurement UpdateEvaluate the sigma points (i=1,…,nx)𝒳0(k∣k-1)=x̂(k∣k-1),𝒳i(k∣k-1)=x̂(k∣k-1)+((nx+λ)P(k∣k-1))i,𝒳nx+i(k∣k-1)=x̂(k∣k-1)-((nx+λ)P(k∣k-1))i.Evaluate the weights (i=1,…,2nx)W0(m)=λnx+λ,W0(c)=λnx+λ+(1-α2+β),Wi(m)=Wi(c)=12(nx+λ).Evaluate the propagated sigma points (i=0,…,2nx)𝒴i(k∣k-1)=h(𝒳i(k∣k-1)).Estimate the predicted measurementŷ(k∣k-1)=∑i=02nxWi(m)𝒴i(k∣k-1).Estimate the innovation covariance and cross-covariance matrix (i=0,…,2nx)𝒳̃i(k∣k-1)=𝒳i(k∣k-1)-x̂(k∣k-1),𝒴̃i(k∣k-1)=𝒴i(k∣k-1)-ŷ(k∣k-1),Pyy(k∣k-1)=∑i=02nxWi(c)𝒴̃i(k∣k-1)𝒴̃i⊤(k∣k-1)+R,Pxy(k∣k-1)=∑i=02nxWi(c)𝒳̃i(k∣k-1)𝒴̃i⊤(k∣k-1).Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-ŷ(k∣k-1)-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k),wherenx is the dimension of state vector, λ=α2(nx+κ)-nx is a scalingparameter,α determines the spread of the sigma points, κ is a secondary scalingparameter andβ is a weight parameter. (A)i denotes the ith column of thematrixA


## 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 p(x(k)∣Y(k-1)) and p(y(k)∣Y(k-1)) in (2) and (4) are assumed Gaussian, solving approximately the multidimensional integrals is equivalent to the evaluation of the corresponding means and covariances of p(x(k)∣Y(k-1)) andp(y(k)∣Y(k-1)). 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 nonlinearfunction×Gaussiandensity [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 sigma-point 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 α=1,β=0, and κ=0 in Algorithm 2.Algorithm 3:Modified Cubature Kalman Filter.
Initializationx̂(0∣0)=𝔼{x(0)}P(0∣0)=𝔼{(x(0)-x̂(0∣0))(x(0)-x̂(0∣0))⊤}Time UpdateEstimate the predicted state vector and error covariance matrixx̂(k∣k-1)=Φx̂(k-1∣k-1),P(k∣k-1)=ΦP(k-1∣k-1)Φ⊤+ΓQΓ⊤.Measurement UpdateEvaluate the cubature points (i=1,…,2nx)P(k∣k-1)=S(k∣k-1)S⊤(k∣k-1),Ξ=nx·[Inx,-Inx],𝒳i(k∣k-1)=x̂(k∣k-1)+S(k∣k-1)ξi.Evaluate the weights (i=1,…,2nx)Wi(m)=Wi(c)=1/(2nx).Evaluate the propagated cubature points (i=1,…,2nx)𝒴i(k∣k-1)=h(𝒳i(k∣k-1)).Estimate the predicted measurementŷ(k∣k-1)∑i=12nxWi(m)𝒴i(k∣k-1)Estimate the innovation covariance and cross-covariance matrixPyy(k∣k-1)=∑i=12nxWi(c)𝒴i(k∣k-1)𝒴i⊤(k∣k-1)-ŷ(k∣k-1)ŷ⊤(k∣k-1)+R,Pxy(k∣k-1)=∑i=12nxWi(c)𝒳i(k∣k-1)𝒴i⊤(k∣k-1)-x̂(k∣k-1)ŷ⊤(k∣k-1)Estimate the Kalman gain, updated state vector and error covariance matrixK(k)=Pxy(k∣k-1)Pyy-1(k∣k-1),x̂(k∣k)=x̂(k∣k-1)+K(k)[y(k)-ŷ(k∣k-1)-μ],P(k∣k)=P(k∣k-1)-K(k)Pyy(k∣k-1)K⊤(k).whereξi denotes the ith column vector of the matrix Ξ.


## 4.5. Posterior Cramér-Rao 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ér-Rao lower bound for the hybrid localization problem is presented that gives the best achievable performance for nonlinear filtering [16, 29]. Let x̂(k∣k) be an unbiased estimate of the state vector x(k). Then, the covariance matrix of the estimation error satisfies the inequality(21)𝔼{(x̂(k∣k)-x(k))(x̂(k∣k)-x(k))⊤}≥J-1(k),
where 𝔼{·} is the expectation with respect to x(k), J(k) denotes the filtering information matrix, and its inverse is the PCRLB matrix. The matrix inequality A≥B should be interpreted as the matrix A-B being positive semidefinite. The aim is now to calculate J(k). In [29], an elegant method is presented, where J(k) can be determined recursively. This recursion, adapted to the hybrid localization problem involving additive Gaussian noise (cf. (5) and (18)) can be written as(22)J(k)=(ΓQΓ⊤+ΦJ-1(k-1)Φ⊤)-1+𝔼{H̃⊤(k)R-1H̃(k)},
where the expectation is with respect to x(k) and H̃(k) denotes the Jacobian matrix of the nonlinear measurement function h(·) (cf. (19)) evaluated at the true value of the state x(k). Since the initial distribution p(x(0)) is assumed to be Gaussian, the recursions are initialized with the information matrix J(0)=P-1(0∣0) [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 floating-point operations (FLOPs). A FLOP is here defined as one addition, subtraction, multiplication, or division of two floating-point numbers. In Table1, 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.Table 1
Computational complexity of some common matrix operations [30].
OperationSizeMult.Add.OtherA+AA∈ℝn×m—nm—A·BA∈ℝn×m,B∈ℝm×llmn(m-1)ln—C-1C∈ℝn×nn3——CC∈ℝn×n——n3/3+2n2In 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 stepk the Jacobian matrix H(k) and the nonlinear function h(·) (cf. Algorithm 1). In the RBUKF and MCKF, one has to propagate at every time step 2nx+1 sigma points and 2nx cubature points through the nonlinear function h(·) (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 Table2, 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(23)𝒞EKF(nx,ny)=8nx3+ny3+6nx2ny+6ny2nx-13nx2-2nxny+9nx+2ny,𝒞RBUKF(nx,ny)=373nx3+ny3+6nx2ny+8ny2nx-9nx2+7nxny+4ny2+10nx+5ny,𝒞MCKF(nx,ny)=373nx3+ny3+6nx2ny+8ny2nx-11nx2+4nxny+3ny2+9nx+2ny,
where nx and ny denote the dimension of the state and measurement vector, respectively.Table 2
Computational complexity of the EKF, RBUKF, and MCKF. The dimension of the state vector is given bynx and the dimension of the measurement vector is given by ny. 𝒳(k∣k-1) denotes the matrix composed of sigma/cubature point vectors.
QuantityComplexityEKFRBUKFx̂(k∣k-1)2nx2-nx2nx2-nx2nx2-nxP(k∣k-1)8nx3-15nx2+10nx8nx3-15nx2+10nx8nx3-15nx2+10nxX(k∣k-1)—13nx3/3+2nx213nx3/3+2nx2ŷ(k∣k-1)—2nxny+2ny2nxnyPyy(k∣k-1)2ny2nx+2nx2ny-nxny4ny2nx+2nxny+4ny2+ny4ny2nx+3ny2Pxy(k∣k-1)2nx2ny-nxny4nx2ny+3nxny+2nx2+nx4nx2ny+2nxnyK(k)ny3+2ny2nx-nxnyny3+2ny2nx-nxnyny3+2ny2nx-nxnyx̂(k∣k)2nxny+2ny2nxny+2ny2nxny+2nyP(k∣k)2ny2nx+2nx2ny-nxny2ny2nx+2nx2ny-nxny2ny2nx+2nx2ny-nxny


## 5. Simulation Results

### 5.1. Scenario I

In the first simulation scenario (Scenario I), it is assumed that a car is equipped with an MT that is capable of providing PR-measured values from GPS and TA, RSS, and GPS reference time uncertainty-measured values from GSM. The car moves with a constant speed of45km/h in a dense urban scenario of size 3km×3km as it is shown in Figure 1. The GSM network is composed of NBS=7 BSs, where each BS is equipped with a directional antenna. The BS locations as well as the BS antenna parameters are a priori known. The satellite locations are taken from the real GPS satellite constellation taking into account realistic satellite elevation masks and are assumed to be known. The parameters used in the simulations are summarized in Table 3 and are assumed to be equal for all BSs and all satellites for the sake of simplicity. The following combinations of measured values are investigated:  (i)
GSM method: one TA-measured value from the serving BS and a total of seven RSS-measured values from serving and neighbouring BS antennas,(ii)
Hybrid 1 method: measured values of GSM method and, in addition, one PR-measured value from one satellite,(iii)
Hybrid 2 method: measured values of GSM method and, in addition, two PR-measured values from two different satellites.Table 3
Simulation parameters.
ParameterValueParameterValueA in dB132.8σx in m/s210-2B in dB3.8σy in m/s210-2Pt in dBm50c0 in m/s3·108Am in dB20Ts in s0.48φ3dB in ∘60h09.4·10-20σRSS in dB8h-11.8·10-19σTA in m300h-23.8·10-21μTA in m0α10-3σPR in m15β2μPR in m0κ0Figure 1
Simulation Scenario I withNBS=7 BSs (•). The arrows (→) indicate the BS antenna boresight direction.Here, it is worth noting that the case when more than two PR-measured values are available has been investigated in [12]. For simplicity, the serving BS is assumed to be the BS located at [750m,1000m]⊤. Here, it is worth noting that the serving BS has been placed very close to the MT trajectory in order to better illustrate how the nonlinearity, inherent in the TA- and RSS-measured values, influences the performance of the different filters. The PR-, TA-, and RSS-measured values are updated every Ts=0.48s, which corresponds to the reporting period of measured values in GSM networks. The performance of the proposed algorithms for the hybrid localization method is evaluated in terms of the root mean square error (RMSE) determined from NMC=500 Monte Carlo trials [16]. For each Monte Carlo trial, the MT trajectory is generated using (5) with process noise parameters as given in Table 3 and initial state vector x(0)=[-200m,8.84m/s,-200m,8.84m/s,0m,0m]⊤. For the simulations, the initial error covariance matrix of the three filters is set to P(0∣0)=diag((200m)2, (10m/s)2, (200m)2, (10m/s)2, (300km)2/3, (10m)2/3), and the initial state vector x̂(0∣0) is obtained from random initialization [20].  The covariance matrix Q for the three filters, which is a design parameter, is chosen to be Q=diag(100·Qcv,10·Qδt), in order to account for possible MT maneuvers and receiver clock uncertainties. The measurement covariance matrix R for the simulations and the filters is assumed to be the same.


### 5.2. Scenario II

The second simulation scenario (Scenario II) investigates the tracking performance of the three different filters for a more general and realistic scenario. It is assumed that a car is equipped with an MT that is capable of providing PR-measured values from GPS and TA, RSS from GSM. Here, it is worth noting that the results for taking into account the GPS reference time uncertainty measurements are very similar to those of Scenario I and, thus, are not further elaborated. The car moves clockwise on a trapezoidal route, divided into 4 sections, in an urban scenario of size5km×5km as it is shown in Figure 2. In each section, the car moves with a different velocity in order to reflect a more realistic car movement as depicted in Figure 3. The GSM network is composed of NBS=12 BSs, where each BS is equipped with three directional antennas. The BS and satellite locations as well as the BS antenna parameters are assumed a priori known. The combination of measurements investigated and the parameters used in the simulations is the same as in Scenario I. The serving BS is assumed to be the BS antenna providing the largest RSS-measured value. The performance of the proposed algorithms is evaluated from NMC=500 Monte Carlo trials. For each Monte Carlo trial, the MT trajectory is generated based on Figure 3 and initial state vector x(0)=[500m,8m/s,1500m,0m/s,0m,0m]⊤. For the simulations, the initial state vector x̂(0∣0) is obtained from random initialization and the initial error covariance matrix of the three filters is set to P(0∣0)=diag((200m)2,(10m/s)2,(200m)2,(10m/s)2, (300km)2/3, (10m)2/3). The covariance matrix Q and the measurement covariance matrix R for the three filters are chosen as in Scenario I.Figure 2
Simulation Scenario II withNBS=12 BSs (•). Each BS is equipped with three directional antennas. The arrows (→) indicate the BS antenna boresight direction. The MT starts at xMT(0)=[500,1500]⊤ and moves clockwise on the trapezoidal route.Figure 3
MT speed versus time indexk for Scenario II.


### 5.3. Simulation Results for Scenario I without GPS Reference Time Uncertainty Measurements

In this section, the simulation results for the case when there are no GPS reference time uncertainty measurements available from the GSM network are presented. In Figure4, the simulation results for the MT location RMSE in dependence of the time index k for the GSM method are shown. From Figure 4, it can be seen that during the first 200 time steps, the performance of the three filters is approximately the same. However, when the MT is located close to the serving BS, there is a “high degree” of nonlinearity in the TA- and RSS-measured values. In this region, it can be clearly seen that the MCKF outperforms the EKF and RBUKF in terms of RMSE.Figure 4
Scenario I: MT location RMSE of the GSM method for the different filters and the corresponding PCRLB.In Figures5 and 6, the MT location RMSE for the Hybrid 1 and Hybrid 2 method for the different filters is shown. Compared to the GSM method, the MT location RMSE can only be marginally improved by the Hybrid 1 method, which additionally takes into account one PR-measured value from GPS. This can be explained by the fact that it is not possible to accurately estimate the unknown MT clock bias with the available TA- and RSS-measured values, which is depicted in Figure 7. Furthermore, when looking at the PCRLBs of the GSM and Hybrid 1 method for the MT location, one can conclude that from a theoretical point of view no RMSE improvements are possible, since the PCRLBs practically coincide with each other.Figure 5
Scenario I: MT location RMSE of the Hybrid 1 method for the different filters and the corresponding PCRLB.Figure 6
Scenario I: MT location RMSE of the Hybrid 2 method for the different filters and the corresponding PCRLB.Figure 7
Scenario I: MT clock bias RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding PCRLB.For the Hybrid 2 method (cf. Figure6) the improvements are significant. Due to the fact that two PR-measured values are available, the filters can much more accurately estimate the MT clock bias (cf. Figure 7) which has a direct impact on the achievable MT location RMSE. However, from Figure 6 it can be also seen that the performance of the three filters is approximately the same. This can be explained by the fact that the two PR-measured values have a larger influence on the MT location estimate of the filters, because these measured values are more accurate than the TA- and RSS-measured values. Since the distances between the MT and the satellites are very large, and the nonlinearity of the PR-measured values is thus “mild”, more weight is put on the PR-measured values when the MT location estimates are evaluated in the filters. As a result, the influence of the nonlinearities of the TA- and RSS-measured values is weakened and, thus, for the Hybrid 2 method, the three filters have approximately the same performance. In Table 4, the performance of the EKF is compared to the performance of the RBUKF and MCKF in terms of the RMSE averaged over the whole time period. From Table 4 it can be concluded that for this specific scenario, the MCKF slightly outperforms the RBUKF and EKF in terms of average RMSE.Table 4
Average RMSE performance of EKF, RBUKF, and MCKF for Scenario I.
AlgorithmMethodLocation in mVelocity in m/sBias in mDrift in m/sGSM73.93.6——EKF [12]Hybrid 173.03.434.21.9Hybrid 244.32.74.61.2GSM73.83.6——RBUKF [13]Hybrid 172.13.533.71.9Hybrid 244.22.74.61.2GSM71.53.6——MCKFHybrid 170.83.433.51.9Hybrid 244.02.74.61.2


### 5.4. Simulation Results for Scenario I with GPS Reference Time Uncertainty Measurements

In this section, simulation results for the case when there are GPS reference time uncertainty measurements available from the GSM network are presented. Here, the important question is investigated, what accuracy of the GPS reference time measurements is needed, in order to improve the performance of the MT location RMSE.In Figure8, the average MT location RMSE in dependence of the GPS reference time uncertainty standard deviation σRTU for the different filters are shown. From Figure 8, it can be seen that for the Hybrid 2 method no significant improvements are possible and the performance of the different filters is approximately the same. That is, when two PR-measured values are additionally available, then, the additional evaluation of the GPS reference time information from the GSM network does not help to improve the MT location RMSE of the Hybrid 2 method. However, the improvements for the Hybrid 1 method can be significant, when GPS reference time information is taken into account. For the investigated scenario, the major performance improvements for the Hybrid 1 method can be obtained for values of σRTU that lie in the range between 10-7 and 10-6 seconds. For GPS reference time uncertainty values outside this region, no significant improvements in terms of average MT location RMSE can be achieved. However, for the case when σRTU<10-7, then, the performance of the Hybrid 1 method reaches quickly a lower bound which is approximately equal to the performance of the Hybrid 2 method.Figure 8
Scenario I: Average MT location RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.When the different filters are compared with each other, one can clearly see that for large values ofσRTU, the MCKF outperforms the RBUKF and EKF. However, when the GPS reference time uncertainty decreases, then, the performance of the three filters is approximately the same. In Figures 9 and 10, the average RMSEs of the MT clock bias and drift in dependence of the GPS reference time uncertainty standard deviation σRTU for the different filters are shown. For the Hybrid 1 method, the average bias and drift RMSEs can be gradually decreased for values of σRTU smaller than 10-6. For the Hybrid 2 method, this is true for values of σRTU smaller than 10-7. From the achieved results one can conclude that the MT location RMSE can be decreased until the MT clock bias RMSE falls below a certain threshold. Beyond this point, the MT location RMSE reaches a lower bound, even though the RMSEs of the MT clock bias and drift states can be further decreased.Figure 9
Scenario I: Average MT clock bias RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.Figure 10
Scenario I: Average MT clock drift RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.In Table5, the computational complexity in terms of FLOPs for the different filters and methods is presented. From Table 5, it can be seen that for all investigated methods, the EKF has the lowest computational complexity, followed by the MCKF and RBUKF. The complexity reduction of the EKF compared to the RBUKF is about 30%. Using an MCKF rather than an RBUKF results in a complexity reduction of only 4%.Table 5
Computational complexity in FLOPs of EKF, RBUKF, and MCKF for one time step. Numbers in parentheses denote the FLOP complexity taking into account GPS reference time information.
AlgorithmMethodComplexity in FLOPsGSM3108EKFHybrid 16813(7974)Hybrid 27974(9267)GSM4534RBUKFHybrid 19708(11230)Hybrid 211230(12916)GSM4314MCKFHybrid 19360(10842)Hybrid 210842(12486)


### 5.5. Simulation Results for Scenario II

In this section, the simulation results for Scenario II are presented. In Figures11, 12, and 13, the MT location RMSE of the GSM, Hybrid 1, and Hybrid 2 method for the three different filters is shown. From these three figures it can be clearly seen that the performance of the three filters is approximately the same. Again, the GSM method provides the worst results in terms of MT location RMSE. The Hybrid 1 method only marginally improves the MT location RMSE at the cost of an increased complexity (cf. Table 5) while the best performance can be achieved by the Hybrid 2 method.Figure 11
Scenario II: MT location RMSE of the GSM method for the different filters.Figure 12
Scenario II: MT location RMSE of the Hybrid 1 method for the different filters.Figure 13
Scenario II: MT location RMSE of the Hybrid 2 method for the different filters.The equal performance of the three filters for the different methods can be explained by the fact that the distances between the BSs and the MT are large and, thus, the impact of the nonlinearities, inherent in the TA and RSS measured values, is small. Compared to the simulation results of Scenario I, there are now several distinct peaks in the MT location RMSE. These peaks that are more pronounced for the GSM and Hybrid 1 method (cf. Figures11 and 12) result from the geometric constellation of the BSs and satellites relative to the MT. However, for the Hybrid 2 method (cf. Figure 13) the magnitude of the peaks becomes smaller, because additional information from two PR-measured values is processed and more weight is put on the geometric constellation of the satellites relative to the MT.  The average RMSE performance of the EKF, RBUKF, and MCKF is summarized in Table 6.Table 6
Average RMSE performance of EKF, RBUKF, and MCKF for Scenario II.
AlgorithmMethodLocation in mVelocity in m/sBias in mDrift in m/sGSM119.14.1——EKF [12]Hybrid 1118.43.957.12.0Hybrid 269.92.84.51.1GSM118.84.1——RBUKF [13]Hybrid 1118.23.956.82.0Hybrid 270.02.84.51.1GSM118.54.1——MCKFHybrid 1118.13.956.72.0Hybrid 270.02.84.51.1


## 5.1. Scenario I

In the first simulation scenario (Scenario I), it is assumed that a car is equipped with an MT that is capable of providing PR-measured values from GPS and TA, RSS, and GPS reference time uncertainty-measured values from GSM. The car moves with a constant speed of45km/h in a dense urban scenario of size 3km×3km as it is shown in Figure 1. The GSM network is composed of NBS=7 BSs, where each BS is equipped with a directional antenna. The BS locations as well as the BS antenna parameters are a priori known. The satellite locations are taken from the real GPS satellite constellation taking into account realistic satellite elevation masks and are assumed to be known. The parameters used in the simulations are summarized in Table 3 and are assumed to be equal for all BSs and all satellites for the sake of simplicity. The following combinations of measured values are investigated:  (i)
GSM method: one TA-measured value from the serving BS and a total of seven RSS-measured values from serving and neighbouring BS antennas,(ii)
Hybrid 1 method: measured values of GSM method and, in addition, one PR-measured value from one satellite,(iii)
Hybrid 2 method: measured values of GSM method and, in addition, two PR-measured values from two different satellites.Table 3
Simulation parameters.
ParameterValueParameterValueA in dB132.8σx in m/s210-2B in dB3.8σy in m/s210-2Pt in dBm50c0 in m/s3·108Am in dB20Ts in s0.48φ3dB in ∘60h09.4·10-20σRSS in dB8h-11.8·10-19σTA in m300h-23.8·10-21μTA in m0α10-3σPR in m15β2μPR in m0κ0Figure 1
Simulation Scenario I withNBS=7 BSs (•). The arrows (→) indicate the BS antenna boresight direction.Here, it is worth noting that the case when more than two PR-measured values are available has been investigated in [12]. For simplicity, the serving BS is assumed to be the BS located at [750m,1000m]⊤. Here, it is worth noting that the serving BS has been placed very close to the MT trajectory in order to better illustrate how the nonlinearity, inherent in the TA- and RSS-measured values, influences the performance of the different filters. The PR-, TA-, and RSS-measured values are updated every Ts=0.48s, which corresponds to the reporting period of measured values in GSM networks. The performance of the proposed algorithms for the hybrid localization method is evaluated in terms of the root mean square error (RMSE) determined from NMC=500 Monte Carlo trials [16]. For each Monte Carlo trial, the MT trajectory is generated using (5) with process noise parameters as given in Table 3 and initial state vector x(0)=[-200m,8.84m/s,-200m,8.84m/s,0m,0m]⊤. For the simulations, the initial error covariance matrix of the three filters is set to P(0∣0)=diag((200m)2, (10m/s)2, (200m)2, (10m/s)2, (300km)2/3, (10m)2/3), and the initial state vector x̂(0∣0) is obtained from random initialization [20].  The covariance matrix Q for the three filters, which is a design parameter, is chosen to be Q=diag(100·Qcv,10·Qδt), in order to account for possible MT maneuvers and receiver clock uncertainties. The measurement covariance matrix R for the simulations and the filters is assumed to be the same.


## 5.2. Scenario II

The second simulation scenario (Scenario II) investigates the tracking performance of the three different filters for a more general and realistic scenario. It is assumed that a car is equipped with an MT that is capable of providing PR-measured values from GPS and TA, RSS from GSM. Here, it is worth noting that the results for taking into account the GPS reference time uncertainty measurements are very similar to those of Scenario I and, thus, are not further elaborated. The car moves clockwise on a trapezoidal route, divided into 4 sections, in an urban scenario of size5km×5km as it is shown in Figure 2. In each section, the car moves with a different velocity in order to reflect a more realistic car movement as depicted in Figure 3. The GSM network is composed of NBS=12 BSs, where each BS is equipped with three directional antennas. The BS and satellite locations as well as the BS antenna parameters are assumed a priori known. The combination of measurements investigated and the parameters used in the simulations is the same as in Scenario I. The serving BS is assumed to be the BS antenna providing the largest RSS-measured value. The performance of the proposed algorithms is evaluated from NMC=500 Monte Carlo trials. For each Monte Carlo trial, the MT trajectory is generated based on Figure 3 and initial state vector x(0)=[500m,8m/s,1500m,0m/s,0m,0m]⊤. For the simulations, the initial state vector x̂(0∣0) is obtained from random initialization and the initial error covariance matrix of the three filters is set to P(0∣0)=diag((200m)2,(10m/s)2,(200m)2,(10m/s)2, (300km)2/3, (10m)2/3). The covariance matrix Q and the measurement covariance matrix R for the three filters are chosen as in Scenario I.Figure 2
Simulation Scenario II withNBS=12 BSs (•). Each BS is equipped with three directional antennas. The arrows (→) indicate the BS antenna boresight direction. The MT starts at xMT(0)=[500,1500]⊤ and moves clockwise on the trapezoidal route.Figure 3
MT speed versus time indexk for Scenario II.


## 5.3. Simulation Results for Scenario I without GPS Reference Time Uncertainty Measurements

In this section, the simulation results for the case when there are no GPS reference time uncertainty measurements available from the GSM network are presented. In Figure4, the simulation results for the MT location RMSE in dependence of the time index k for the GSM method are shown. From Figure 4, it can be seen that during the first 200 time steps, the performance of the three filters is approximately the same. However, when the MT is located close to the serving BS, there is a “high degree” of nonlinearity in the TA- and RSS-measured values. In this region, it can be clearly seen that the MCKF outperforms the EKF and RBUKF in terms of RMSE.Figure 4
Scenario I: MT location RMSE of the GSM method for the different filters and the corresponding PCRLB.In Figures5 and 6, the MT location RMSE for the Hybrid 1 and Hybrid 2 method for the different filters is shown. Compared to the GSM method, the MT location RMSE can only be marginally improved by the Hybrid 1 method, which additionally takes into account one PR-measured value from GPS. This can be explained by the fact that it is not possible to accurately estimate the unknown MT clock bias with the available TA- and RSS-measured values, which is depicted in Figure 7. Furthermore, when looking at the PCRLBs of the GSM and Hybrid 1 method for the MT location, one can conclude that from a theoretical point of view no RMSE improvements are possible, since the PCRLBs practically coincide with each other.Figure 5
Scenario I: MT location RMSE of the Hybrid 1 method for the different filters and the corresponding PCRLB.Figure 6
Scenario I: MT location RMSE of the Hybrid 2 method for the different filters and the corresponding PCRLB.Figure 7
Scenario I: MT clock bias RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding PCRLB.For the Hybrid 2 method (cf. Figure6) the improvements are significant. Due to the fact that two PR-measured values are available, the filters can much more accurately estimate the MT clock bias (cf. Figure 7) which has a direct impact on the achievable MT location RMSE. However, from Figure 6 it can be also seen that the performance of the three filters is approximately the same. This can be explained by the fact that the two PR-measured values have a larger influence on the MT location estimate of the filters, because these measured values are more accurate than the TA- and RSS-measured values. Since the distances between the MT and the satellites are very large, and the nonlinearity of the PR-measured values is thus “mild”, more weight is put on the PR-measured values when the MT location estimates are evaluated in the filters. As a result, the influence of the nonlinearities of the TA- and RSS-measured values is weakened and, thus, for the Hybrid 2 method, the three filters have approximately the same performance. In Table 4, the performance of the EKF is compared to the performance of the RBUKF and MCKF in terms of the RMSE averaged over the whole time period. From Table 4 it can be concluded that for this specific scenario, the MCKF slightly outperforms the RBUKF and EKF in terms of average RMSE.Table 4
Average RMSE performance of EKF, RBUKF, and MCKF for Scenario I.
AlgorithmMethodLocation in mVelocity in m/sBias in mDrift in m/sGSM73.93.6——EKF [12]Hybrid 173.03.434.21.9Hybrid 244.32.74.61.2GSM73.83.6——RBUKF [13]Hybrid 172.13.533.71.9Hybrid 244.22.74.61.2GSM71.53.6——MCKFHybrid 170.83.433.51.9Hybrid 244.02.74.61.2


## 5.4. Simulation Results for Scenario I with GPS Reference Time Uncertainty Measurements

In this section, simulation results for the case when there are GPS reference time uncertainty measurements available from the GSM network are presented. Here, the important question is investigated, what accuracy of the GPS reference time measurements is needed, in order to improve the performance of the MT location RMSE.In Figure8, the average MT location RMSE in dependence of the GPS reference time uncertainty standard deviation σRTU for the different filters are shown. From Figure 8, it can be seen that for the Hybrid 2 method no significant improvements are possible and the performance of the different filters is approximately the same. That is, when two PR-measured values are additionally available, then, the additional evaluation of the GPS reference time information from the GSM network does not help to improve the MT location RMSE of the Hybrid 2 method. However, the improvements for the Hybrid 1 method can be significant, when GPS reference time information is taken into account. For the investigated scenario, the major performance improvements for the Hybrid 1 method can be obtained for values of σRTU that lie in the range between 10-7 and 10-6 seconds. For GPS reference time uncertainty values outside this region, no significant improvements in terms of average MT location RMSE can be achieved. However, for the case when σRTU<10-7, then, the performance of the Hybrid 1 method reaches quickly a lower bound which is approximately equal to the performance of the Hybrid 2 method.Figure 8
Scenario I: Average MT location RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.When the different filters are compared with each other, one can clearly see that for large values ofσRTU, the MCKF outperforms the RBUKF and EKF. However, when the GPS reference time uncertainty decreases, then, the performance of the three filters is approximately the same. In Figures 9 and 10, the average RMSEs of the MT clock bias and drift in dependence of the GPS reference time uncertainty standard deviation σRTU for the different filters are shown. For the Hybrid 1 method, the average bias and drift RMSEs can be gradually decreased for values of σRTU smaller than 10-6. For the Hybrid 2 method, this is true for values of σRTU smaller than 10-7. From the achieved results one can conclude that the MT location RMSE can be decreased until the MT clock bias RMSE falls below a certain threshold. Beyond this point, the MT location RMSE reaches a lower bound, even though the RMSEs of the MT clock bias and drift states can be further decreased.Figure 9
Scenario I: Average MT clock bias RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.Figure 10
Scenario I: Average MT clock drift RMSE of the Hybrid 1 and Hybrid 2 method for the different filters and the corresponding average PCRLB.In Table5, the computational complexity in terms of FLOPs for the different filters and methods is presented. From Table 5, it can be seen that for all investigated methods, the EKF has the lowest computational complexity, followed by the MCKF and RBUKF. The complexity reduction of the EKF compared to the RBUKF is about 30%. Using an MCKF rather than an RBUKF results in a complexity reduction of only 4%.Table 5
Computational complexity in FLOPs of EKF, RBUKF, and MCKF for one time step. Numbers in parentheses denote the FLOP complexity taking into account GPS reference time information.
AlgorithmMethodComplexity in FLOPsGSM3108EKFHybrid 16813(7974)Hybrid 27974(9267)GSM4534RBUKFHybrid 19708(11230)Hybrid 211230(12916)GSM4314MCKFHybrid 19360(10842)Hybrid 210842(12486)


## 5.5. Simulation Results for Scenario II

In this section, the simulation results for Scenario II are presented. In Figures11, 12, and 13, the MT location RMSE of the GSM, Hybrid 1, and Hybrid 2 method for the three different filters is shown. From these three figures it can be clearly seen that the performance of the three filters is approximately the same. Again, the GSM method provides the worst results in terms of MT location RMSE. The Hybrid 1 method only marginally improves the MT location RMSE at the cost of an increased complexity (cf. Table 5) while the best performance can be achieved by the Hybrid 2 method.Figure 11
Scenario II: MT location RMSE of the GSM method for the different filters.Figure 12
Scenario II: MT location RMSE of the Hybrid 1 method for the different filters.Figure 13
Scenario II: MT location RMSE of the Hybrid 2 method for the different filters.The equal performance of the three filters for the different methods can be explained by the fact that the distances between the BSs and the MT are large and, thus, the impact of the nonlinearities, inherent in the TA and RSS measured values, is small. Compared to the simulation results of Scenario I, there are now several distinct peaks in the MT location RMSE. These peaks that are more pronounced for the GSM and Hybrid 1 method (cf. Figures11 and 12) result from the geometric constellation of the BSs and satellites relative to the MT. However, for the Hybrid 2 method (cf. Figure 13) the magnitude of the peaks becomes smaller, because additional information from two PR-measured values is processed and more weight is put on the geometric constellation of the satellites relative to the MT.  The average RMSE performance of the EKF, RBUKF, and MCKF is summarized in Table 6.Table 6
Average RMSE performance of EKF, RBUKF, and MCKF for Scenario II.
AlgorithmMethodLocation in mVelocity in m/sBias in mDrift in m/sGSM119.14.1——EKF [12]Hybrid 1118.43.957.12.0Hybrid 269.92.84.51.1GSM118.84.1——RBUKF [13]Hybrid 1118.23.956.82.0Hybrid 270.02.84.51.1GSM118.54.1——MCKFHybrid 1118.13.956.72.0Hybrid 270.02.84.51.1


## 6. Experimental Results

In this section, the performance of the proposed hybrid localization method is verified for the different filters with experimental data available from a field trial. The field trial was conducted in an operating GSM network in the city center of a German city, with a test area of approximately3km×3km. During the field trial, a car equipped with a standard cellular phone is collected every Ts=480ms RXLEV- (quantized RSS) and TA-measured values from GSM. Here, it is worth noting that in GSM the RXLEV-measured values are available from the serving BS and between one and six strongest RXLEVs from the neighboring BS whereas the TA-measured value is only available from the serving BS. In addition, it should be noted that GPS reference time uncertainty measurements have not been collected during the field trial, so that this issue will not be further elaborated in this section. The GSM network is composed of NBS=13 fixed BSs with known locations. The BSs are equipped with either directional antennas or a single omnidirectional antenna. The antenna boresight directions, equivalent isotropic radiated powers, and half-power beamwidths are a priori known, and the unknown antenna gain patterns are approximated with (13). The remaining parameters of the TA and RSS measurement models (cf. (10) and (14)) are estimated from the available field trial data. Here, it is worth noting that in order to not overfit the different filters for this single trajectory the standard deviations of the TA- and RSS-measured values were chosen to be σRSS(n)≥3dB and σTA(n)≥1μs. For the path loss model (cf. (12)) the parameters are in the range of 110dB≤A(n)≤150dB and 2dB≤B(n)≤5dB.For the GPS network, PR-measured values collected from a field trial are not available, so that synthetic PR measurement data have been generated with the parameters given in Table3. The constellation of the GPS satellites during the field trial is reconstructed by taking true satellite locations from the real satellite constellation, where the satellite locations are assumed a priori known. For simplicity, it is assumed that PR-measured values are available every Ts=480ms. The satellite's visibility status during the field trial cannot be reproduced subsequently, so that it is assumed that either NSAT=1 or NSAT=2 satellites are visible to the MT. However, this assumption is only made in order to demonstrate the improvements that can be achieved by the proposed hybrid localization algorithm. In reality, the number of visible satellites changes with time, so that there will be situations where GPS (i.e., NSAT≥3) is available.In Figure14, the true MT trajectory together with the trajectories estimated by the EKF, RBUKF, and MCKF for the GSM method of the field trial is shown. Here, the true MT location was obtained from detailed maps and from GPS, where GPS was available. From Figure 14 it can be seen that all three filters can moderately track the MT. In Figure 15, the corresponding MT location error in dependence of the time index k is presented. Here, the MT location error denotes the Euclidean distance between the true and estimated MT location. From Figure 15 it can be seen that the RBUKF marginally outperforms the MCKF and EKF. This rather surprising result can be explained as follows. The MCKF has been derived in a setting where the measurement noise has to be Gaussian distributed [14]. The available field trial measurements, however, are not Gaussian distributed, so that this assumption is explicitly violated. The RBUKF has been derived in a more general setting where the involved densities have to be symmetric, but not necessarily Gaussian. As a result, the performance of the RBUKF is slightly better than the performance of the MCKF. The peak values in the MT location error can be explained by the geometric constellation of the BSs relative to the MT location and the change of the MT velocity during the field trial which results in a mismatch to the CV model that is assumed in the filter's process model.Figure 14
Field trial scenario with true MT trajectory, approximate BS locations (•), and estimated trajectories of the GSM method for the different filters.Figure 15
MT location error of the GSM method for the different filters.In Figures16 and 17, the estimated trajectories and the MT location error of the Hybrid 1 method for the different filters is presented. Again, all three filters can track the MT and the performance of the MT location error is approximately the same. However, compared to the GSM method, the MT location error of the Hybrid 1 method can be slightly improved. In Figures 18 and 19, the filter's estimated trajectories and MT location error of the Hybrid 2 method are shown. It can be clearly seen that, compared to the GSM and Hybrid 1 method, the performance can be significantly improved using the Hybrid 2 method. For the Hybrid 2 method, the MCKF yields the best performance, followed by the RBUKF and EKF.Figure 16
Field trial scenario with true MT trajectory, approximate BS locations (•), and estimated trajectories of the Hybrid 1 method for the different filters.Figure 17
MT location error of the Hybrid 1 method for the different filters.Figure 18
Field trial scenario with true MT trajectory, approximate BS locations (•), and estimated trajectories of the Hybrid 2 method for the different filters.Figure 19
MT location error of the Hybrid 2 method for the different filters.In Table7, the average MT location error of the GSM, Hybrid 1, and Hybrid 2 method for the EKF, RBUKF, and MCKF is summarized. It can be seen that for all three filters the GSM method provides the worst performance. The MT location accuracy can be marginally improved with the Hybrid 1 method and significantly improved with the Hybrid 2 method.Table 7
Average MT location error in meters of EKF, RBUKF, and MCKF for field trial scenario.
AlgorithmMethodGSMHybrid 1Hybrid 2EKF [12]59.255.543.8RBUKF55.152.342.8MCKF56.553.841.7


## 7. Conclusion

In this paper, the performance and computational complexity of three different MT tracking algorithms, namely, the EKF, RBUKF, and MCKF are investigated that combine TA- and RSS-measured values from GSM and one or two PR measurements from GPS. It has been shown by simulations and experiments that, compared to existing GSM-based MT tracking solutions, the location accuracy can be significantly improved by using hybrid GPS/GSM-based MT tracking algorithms, where all three filters have approximately the same performance, and the EKF offers the best trade-off between performance and computational complexity. When GPS reference time is additionally available from the GSM network, then, the hybrid localization method using one PR measurement can be significantly improved. It has been shown by simulations that in this case, it is sufficient to have GPS reference time information available, whose uncertainty lies in the range between10-6 to 10-7 seconds, in order to improve the hybrid localization method.Although this paper is focused on a specific scenario with measured values from GPS and GSM, the general equations and methodology presented in this paper can be easily adapted to other network configurations and measured values from other systems, for example, universal mobile telecommunications system (UMTS) or wireless local area networks. By additionally considering these measured values in the hybrid localization solution, it is expected that the performance will be further improved.



---

*Source: 149065-2010-08-17.xml*