Multipath modeling and mitigation by using sparse estimation in global navigation satellite system-challenged urban vehicular environments

Multipath interference has been one of the most difficult problems when using global navigation satellite system-based vehicular navigation in urban environments. In this article, we develop a multipath mitigation algorithm exploiting the sparse estimation theory that improves the absolute positioning accuracy in urban environments. The navigation observation model is established by considering the multipath bias as additive positioning errors, and the assumption for the proposed method is that global navigation satellite system signals contaminated due to multipath are the minority among the received signals, which makes the unknown bias vector sparse. We investigated an improved elastic net method to estimate the sparse multipath bias vector, and the global navigation satellite system measurements can be corrected by subtracting the estimated multipath error. The positioning performance of the proposed method is verified by analytical and experimental results.


Introduction
In recent years, with the rapid development of the economy, automobiles have gradually become the main means of transportation in daily life. However, as the number of cars is increasing, the safety and environmental problems caused by urban road congestion and frequent traffic accidents have become a concern. Vehicular navigation technologies have become a vital component in the intelligent transportation system for they provide convenience, location-related information, and safety to transportation users. 1,2 Global navigation satellite system (GNSS) is an all-weather, worldwide, continuous coverage, satellitebased radio navigation system. [3][4][5] As one of the most essential component technologies, GNSS provides precise positions of the user vehicle and neighboring vehicles in urban environments.
The urban environment presents great challenges to common commercial GNSS receivers while vehicular positioning. [6][7][8][9] This is mainly because the GNSS positioning performance can be severely degraded by the multipath effect due to shadowing behind increasing number of buildings, interference, and other undesired impairments such as trees. [10][11][12] As illustrated in Figure 1, multipath interference is generated by receiving reflected or diffracted signals, possibly signals other than line-of-sight signals, and is one of the most harmful sources of error in GNSS positioning applications. Several efforts have been devoted to mitigate the effect of multipath interference from different aspects, for example, the antenna design techniques, 13,14 the receiver-based techniques, 15 as well as the post-receiver techniques, 16 which help to improve accuracy and reliability of the GNSS positioning in urban environment. But these techniques can be expensive and bulky or require deploying additional hardware.
In this article, we proposed a multipath mitigation method exploiting the sparsity properties of channels affected by multipath to estimate multipath biases. Compared with previous works, the proposed method is based on fusing GNSS pseudoranges and pseudorange rates measurements and has an advantage that it does not require additional hardware. Giremus et al. investigated Sequential Monte Carlo methods, also referred to as particle filters, for this estimation. 17 But these approaches are computationally intensive, making a real-time implementation complicated in practical applications. In the proposed method, the computational time for solving the sparse problem is fast, which can allow real-time positioning results. The main hypothesis making this method available is that most satellites are not affected by multipath prompting the biases sparse with respect to number of received measurements. Note that the sparse estimation theory has been considered in turbine and turbofan engine theory, 18,19 however, the proposed method results from the application of a penalized least-squares approach exploiting the sparse estimation for GNSS applications. The navigation observation model is established by considering the multipath bias as additive positioning error, we can estimate the multipath by solving an improved elastic net problem, then the GNSS measurements can be corrected by removing the estimated bias vector and the positioning accuracy will be improved.
This article is organized as follows. The second section summarizes some basic principles on satellite navigation, describing how measurements (pseudoranges and pseudorange rates) are related to the state vector (position, velocity, clock offset, clock drift) and the possible multipath biases. In the third section, the improved elastic net problem is described, and the solution that can be used to estimate multipath bias vector is investigated. The fourth section discusses the experimental results, and the performance of the proposed method is evaluated. Finally, relevant conclusions are drawn in the fifth section.

State model
The navigation solution is obtained by solving these equations with observations from at least four satellites. The state equation can be written as In equation (1), the unknown state vector at time k (to be estimated) is described as where ðx k ; y k ; z k Þ T is the vehicle position in Earth-centered Earth-fixed (ECEF) coordinate, ðx 0 k ; y 0 k ; z 0 k Þ T is the vehicle velocity in ECEF coordinate, d k is the GPS receiver clock offset, and d 0 k is the GPS receiver clock offset drift.
The state transition matrix F k can be described as where Dt is the duration between t k and t kÀ1 , and the process noise vector w k is considered to be zero-mean Gaussian with known covariance matrix.

Observation model
To estimate the unknown state vector X k , we will use the measurement vector comprises the pseudoranges and pseudorange rates. Denote as s k the number of satellites visible at time instant k. The number of measurements acquired by the receiver is 2s k , namely s k pseudoranges denoted as r 1;k ; . . . ; r s k ;k and s k pseudorange rates denoted as r 0 1;k ; . . . ; r 0 s k ;k . These measurements are gathered in the vector z k ¼ ðz k;1 ; : : : ; z k;2s k Þ T 2 R 2s k whose components are z k;i ¼ r i;k and z k;iþs k ¼ r 0 i;k for i ¼ 1; : : : ; s k ð5Þ More precisely, using the notations r i;k ¼ ðx i;k ; y i;k ; z i;k Þ T and v i;k ¼ ðx 0 i;k ; y 0 i;k ; z 0 i;k Þ T for the ith satellite position and velocity at time instant k, we can obtain where r i;k ¼ ðx i;k ; y i;k ; z i;k Þ T is the ith satellite position at time k expressed in the same frame as r k , v i;k ¼ ðx 0 i;k ; y 0 i;k ; z 0 i;k Þ T is the ith satellite velocity at time instant k expressed in the same frame as v k , is the geometric distance between the user and the ith satellite, b k and b 0 k are the receiver clock bias and receiver clock drift at time instant k, respectively, and e i;k and e 0 i;k are the error terms associated with the ith propagation channel (modeling ionospheric delay, tropospheric delay, satellite clock bias, satellite position uncertainty, relativistic effects, multipath and receiver noise).
Note that b k and b 0 k are dependent of i and that e i;k and e 0 i;k summarize all the error sources affecting the s k corresponding measurements. After applying correction models for corresponding error except multipath, the measurement equations can be expressed as where m k ¼ ðm 1;k ; : : : ; m 2s k ;k Þ T 2 R 2s k is the vector accounting for the potential multipath biases, n k ¼ ðn 1;k ; : : : ; n 2s k ;k Þ T 2 R 2s k is the receiver noise vector supposed zero-mean and Gaussian with known covariance matrix R k 2 R 2s k Â2s k , and H k is the system transition matrix and can be described as follows Since the measurement equation is nonlinear, it is natural to use the extended Kalman filter (EKF) to estimate the state vector X k .

EKF for navigation
The classical EKF was used to estimate the state vector from the state equation (1) and the measurements can be summarized as followŝ Proposed GNSS multipath mitigation method The improved elastic net problem Suppose that y ¼ ðy 1 ; y 2 ; : : : ; y n Þ T be the response vector and D ¼ ðD 1 ; D 2 ; : : : ; D p Þ be the predictor matrix, where D j ¼ ðd 1j ; d 2j ; : : : ; d n j Þ T ; j ¼ 1; : : : ; p are predictors. Without loss of generality we assume that the response y is centered and the columns of X are standardized, that is, Hence the usual linear model without intercept term is given by The elastic net solves the following problem (16) can be equivalently represented aŝ where l > 0 and a 2 ð0; 1Þ are the new model parameters. akbk 1 þ ð1 À aÞkbk 2 is the elastic net penalty, 20 and is a compromise between the ridge regression penalty (a ¼ 0) and the lasso penalty (a ¼ 1). To achieve oracle properties, the adaptive elastic net 21 is developed by combining the elastic net and the adaptive lasso, that is where theŵ j ¼ jb j ðenÞj Àl for j ¼ 1; 2; : : : ; p and g is the positive constant. Li et al. investigated a so-called partly adaptive elastic net defined as follows 22 In this article, we propose an improved elastic net model with weighting matrix In addition, we give the algorithm for solving the optimization problem proposed in elastic net numerically; some popular algorithms like last angle regression (LARS) 23 and coordinate descent 24 can be used to solve the optimization problem in the proposed method efficiently. Friedman et al. compared the average CPU timings for the coordinate descent and LARS algorithm for linear regression, 25 and the results showed that the former seems to be competitive and faster than the latter. Hence, we tend to use the coordinate descent algorithm to solve the proposed sparse problem. Since the MATLAB package "glmnet" is publicly available, we prefer this algorithm for solving the elastic net problem. The run times were carried out on an Intel Core 2.80 GHz processor and 8 GB RAM memory.

Solution of multipath mitigation model
Thanks to the sparse estimation theory introduced above, our proposed positioning algorithm can estimate the multipath bias on GNSS pseudoranges and pseudorange rates measurements. This method assumes that the bias vector m k is sparse, that is, some of its components are exactly equal to 0. In order to estimate the bias vector m k and feed it to the proposed EKF, we propose to solve the following problem Regarding the weighting matrix W k , we propose to use some key parameters, somehow representative of the measurements quality, provided by most GNSS receivers, namely the carrier-to-noise density ratio C=N 0 and the satellite elevations which reflect the good or bad positions of the different satellites. 26 More precisely, we want to consider an improved elastic net method favoring satellites with large C=N 0 values and high satellite elevations, which provide reduced tracking noise and little multipath interferene. The C=N 0 weighting function can be described as follows where cn is the value of C=N 0 expressed in dB-Hz, T parameter is a threshold which determines the performance of the weighting function, a defines the bending of the curve, and F defines the value of C=N 0 weighting function and yields to A.
And the satellite elevation that is a second important parameter can be considered as where e is a given satellite elevation expressed in degrees. The multipath of low elevation satellites tends to replace the blocked direct signals and cause large position errors. It is common to reduce the effects of low elevation satellites by defining an elevation mask of 5 before positioning solution. However, due to limited view of the sky in urban environments, we propose to preserve the largest number of satellite measurements, including low elevations. As a consequence, we take different elevations of satellites into account and consider penalizing satellites whose elevations are lower than 5 .
In order to appreciate the impact of low elevation satellites considered in weighting matrix, we tested the performance of the proposed method without penalizing satellite elevations, and the cumulative distribution functions corresponding to the position error are displayed in Figure 2. It can be seen that the weighting with elevation constraint provides a more accurate position solution.
The final weight for a given satellite introduced in the improved elastic net approach is defined as the product of the two previous functions, that is where cn i and e i are the C=N 0 and elevation associated with the pseudorange and pseudorange rate from the ith satellite. In addition, it is easy to acquire C=N 0 and e of each receiver, since C=N 0 is directly estimated by the receiver and that the elevation can be computed using the actual and previous positions of each satellite.
In order to obtain a formulation similar to equation (20), it is interesting to note that the minimization of equation (21) with respect to X k for a fixed m k admits the following least-squares solution After replacing this expression of X k in equation (21), we can reexpress equation (21) as follows where P k ¼ H k ðH T k H k Þ À1 H T k is a projection matrix, and I 2s k is the 2s k Â 2s k identity matrix. Define the following notationsz The original problem (equation (21)) reduces to The proposed multipath mitigation method requires the following steps. First, transform the GNSS observation model into the sparse estimation problem (equation (30)). Based on the carefully selected weighting matrix and parameters, we can estimate the unknown sparse vector x as the solution of the improved elastic net problem. Second, estimate the multipath bias vector m k based on equation (29). Finally, correct the pseudoranges and pseudorange rates measurements by subtracting the estimated errors, and then the corrected measurements are processed using the EKF to obtain improved positioning results. Figure 3 shows the flowchart of the proposed multipath mitigation approach based on sparse estimation.

Simulation results
To appreciate the efficiency of the proposed method, this section studies several experiments conducted with simulated data. It is assumed that the synthetic data are generated using a reference position according to equation (1), and with noisy measurements resulting from equations (6) and (10). In the simulations, the sample size that we apply is K ¼ 200 and have generated artificial additive biases, modeling multipath conditions, affecting pseudoranges and pseudorange rates between time instants k ¼ 60 and k ¼ 140 for satellite channels #2, #3, and #8. These bias amplitudes have been adjusted to 60, À60, and 25 m for pseudoranges and to 40, À32, and À8 m/s for pseudorange rates. Finally, we generated C=N 0 uniformly between 47 dBHz and 50 dBHz in absence of multipath and between 25 dBHz and 28 dBHz in presence of multipath (i.e. in channels #2, #3, and #8).
Currently, the Chinese BeiDou navigation satellite system (BDS) has launched over 40 BDS satellites and begun offering services to customers in the Asia-Pacific region, which make it easy to observe eight satellites in urban environment at one time. In addition, with the joint cooperation of GPS, GLONASS, and GALILEO system, over 12 satellites can be observed to obtain more accurate positioning results. Hence we consider performance of 8 and 12 observed satellites, respectively, in the following experiments.

Synthetic data with eight observed satellites
Firstly, we consider performance of eight satellites to be observed. Figure 4 shows the estimated pseudorange and pseudorange rate biases for three representative channels. It can be seen from these two figures that additional multipath biases are capable of being extracted with accuracy between the disturbed time instant 60 and 140. The positioning root mean square errors (RMSE) of the GNSS state vectors including position, speed, clock bias, and clock drift by using different methods are displayed in Figure 5. This positioning conclusion is confirmed by the results of Figure 6 showing the impact of multipath mitigation method on localization errors. Since the parameter a has an effect on the performance of the proposed method, we compared the performance on different values of parameter a in the proposed method. In this experiment, the performance of EKF, lasso (a¼ 1), and ridge (a¼ 0) methods are compared. To summarize, the results in Figures 5 and 6 illustrate the lasso method has the best accuracy with errors less than 25 m in the four aspects of the state vectors, the ridge method works moderately well while EKF can cause the largest positioning errors as large as a few hundred meters.
Due to the limit of the sparsity assumption in practical applications, the performance of the proposed method when different number of satellites affected by multipath is further illustrated in Figure 7. As can be seen, in the case of eight observed satellites, the proposed sparse method can estimate multipath errors accurately when the number of multipath-contaminated satellites is less than five of eight. However, the performance of the proposed method deteriorates if the number of multipath-polluted satellites is more than six of eight. Because the number of satellites affected by multipath increases gradually, the sparse assumption for this scenario is not appropriate.

Synthetic data with 12 observed satellites
In the case of observations of 12 satellites, the pseudorange and pseudorange rate measurements consisting of the reference and estimated value are depicted in Figure 8.  The simulations conducted with different values of a in the proposed method and EKF methods on localization errors are depicted in Figure 9. As can be seen, all of these results show that the proposed method outperforms the EKF regardless of different values of the parameter and the proposed method with a¼ 1 provides a more accurate position among these estimation strategies. The position RMSE of the EKF, lasso (a¼ 1), and ridge (a¼ 0) methods in the east, north, and up directions are summarized in Table 1.
In the case of 8 or 12 observed satellites, we can see that the lasso method outperforms the EKF and the ridge method for both horizontal and vertical errors. In addition, the position accuracy in the case of 12 observed satellites outperforms that in the case of 8 observed satellites, which shows that multipath can be estimated comparatively accurately with 12 observed satellites, that is, the better performance obtained when the rate of satellites affected by multipath decreases.

Real data
In order to further evaluate the presented multipath mitigation method, a static scenario was considered between twoside buildings, where GPS signals subject to multipath severely causing positioning accuracy reduced. The noisy     L1 pseudoranges and pseudorange rates were used as measurements to the proposed algorithm. The raw GPS data were provided by a u-blox ZED-F9P receiver, and a laptop was used for data logging. The coordinates of the receiver are precisely known by applying the real-time kinematic (RTK) data. Figure 10 shows the number of observed GPS satellites by the receiver during the test, and Figure 11 compares the performance of the proposed multipath mitigation method when a¼ 1 with that of EKF. It can be seen that the proposed method outperforms EKF in terms of accuracy of position estimation, with the RMSE for the proposed and EKF methods of 13.1 m and 44.3 m, respectively. In addition, according to the experimental results, the computing time for the proposed method is about 0.005 s when averaged over three runs.

Conclusion
In this article, a multipath mitigation method exploiting the sparsity estimation is proposed for vehicular positioning in urban environments. We have shown via experiments conducted on synthetic and real data that the improved elastic net method can estimate the multipath bias effectively and outperforms the EKF estimation for improving the GNSSbased positioning. In addition, the computational time for solving the sparse problem of the proposed method is fast, which can allow real-time positioning results. The performance of the proposed method in different number of satellites affected by multipath is also evaluated; the analytical results show some limits for the proposed method when the number of channels affected by multipath is more than six of eight. As future work, more sophisticated solutions such as incorporating with multiple sensors are expected to be designed and evaluated to enhance positioning performance in urban environments.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.