LAND VEHICLE NAVIGATION USING LOW-COST INTEGRATED SMARTPHONE GNSS MEMS AND MAP MATCHING TECHNIQUE

. The demand for smartphone positioning has grown rapidly due to increased positioning accuracy applications, such as land vehicle navigation systems used for vehicle tracking, emergency assistance, and intelligent transportation systems. The integration between navigation systems is necessary to maintain a reliable solution. High-end inertial sensors are not preferred due to their high cost. Smartphone microelectromechanical systems (MEMS) are attractive due to their small size and low cost; however, they suffer from long-term drift, which highlights the need for additional aiding solutions using road network that can perform efficiently for longer periods. In this research, the performance of the Xiaomi MI 8 smartphone's single-frequency precise point positioning was tested in kinematic mode using the between-satellite single-difference (BSSD) technique. A Kalman filter algorithm was used to integrate BSSD and inertial navigation system (INS)-based smartphone MEMS. Map matching technique was proposed to assist navigation systems in global navigation satellite system (GNSS)-denied environments, based on the integration of BSSD–INS and road network models applying hidden Marcov model and Viterbi algorithm. The results showed that BSSD–INS– map performed consistently better than BSSD solution and BSSD–INS integration, irrespective of whether simulated outages were added or not. The root mean square error (RMSE) values for 2D horizontal position accuracy when applying BSSD–INS–map integration improved by 29% and 22%, compared to BSSD and BSSD–INS navigation solutions, respectively, with no simulated outages added. The overall average improvement of proposed BSSD–INS–map integration was 91%, 96%, and 98% in 2D horizontal positioning accuracy, compared to BSSD–INS algorithm for six GNSS simulated signal outages with duration of 10, 20, and 30 s, respectively.


INTRODUCTION
The precise point positioning (PPP) technique has been widely used for the positioning of the global navigation satellite system (GNSS).The benefits of PPP are its low cost, where there is no need for the base receiver, no limitation of distance relative to the differential technique, and providing up to decimeter positioning accuracy in kinematic mode (Hamed et al., 2019).Google announced in May 2016 that raw GNSS measurements from the smartphone Android 7 became available (European, 2017).Many smartphone applications give developers access to code and carrier measurements, leading to the development of GNSS processing techniques using lowcost sensors.These techniques were previously only available to expensive GNSS receivers (GSA, 2020).In May 2018, Xiaomi launched the first dual-frequency smartphone, the Xiaomi MI 8.It has a GNSS chip with the availability of observing codes and carrier phase measurements of L1 frequency on global positioning system (GPS) of the United States and L1 & L5 frequency on the global navigation satellite system of the Russian Federation (GLONASS), and Galileo E1 & E5a (El-Mezayen et al., 2019).
Smartphone's single frequency of GNSS measurements was tested in real-time kinematic (RTK) and static modes (Dabove et al., 2019).The PPP performance of the Xiaomi MI 8 smartphone was evaluated in kinematic mode.The positioning accuracy was at the level of a few meters, and the positioning performance had an offset of more than 4 m from the real track when compared to the geodetic receiver (Wu et al.,2019).
Increased positioning accuracy can support many mass consumer applications such as land vehicle navigation systems (LVNSs).LVNS is used for vehicle tracking, collision avoidance, emergency assistance, and several mobile applications that require increasing location accuracy.Most vehicles travel in urban areas with difficult navigational conditions (tunnels, forests, in between skyscrapers, etc.).GNSS systems are often combined with inertial sensors to maintain a reliable and continuous navigation solution (Karamat, 2014).High-end inertial sensors are not preferred due to their high cost.Smartphone microelectromechanical systems (MEMS) are attractive due to their small size, light weight, and low cost, but they suffer from severe drifts because of inherent biases (Angrisano, 2010).Map matching is a basic operation for improving positioning accuracy by integrating GNSS tracking data with spatial road network data (roadway centerlines).It is utilized to identify the correct road link on which a vehicle is traveling and determine the location of a vehicle on a network road (Attia, 2013).Map matching is useful for intelligent transportation systems and urban traffic modeling by optimizing the traffic flow and avoiding traffic jams (Trough et al., 2020).
LVNS technology based on GNSS and inertial navigation system (INS) is a subject of interest due to its potential for consumers and commercial vehicle markets.The performance of LVNS has been improved over the years by a lot of approaches.The first approach is to improve the stochastic modeling of the MEMS inertial sensors and enhance the quality of the MEMS data.Niu et al. (2010) added constraint update information to the MEMS sensors to improve Kalman filter errors, especially during GNSS signal blockages.The second approach is to utilize the road network data as an additional constraint.Attia (2013) provided several geospatial data models for both indoor and outdoor environments.He introduced a map matching algorithm to match and project navigation positions on the geospatial map links to increase positioning accuracy.
The third method is to use aid sensors and apply a reduced inertial sensor integration.Moussa et al. (2019) introduced an approach for land vehicle navigation in GNSS-denied environments by aiding INS with a low-cost ultrasonic sensor.They used the Kalman filter to bound the INS drift during GNSS blockage based on heading change update for enhancing the reduced integration estimation.The fourth approach is to introduce priori navigation system information (vehicle dynamic models and kinematic restrictions), which can reduce estimation uncertainty.Zhang et al. (2021) developed a velocity model composed ofspeed sensors and motion constraints used as direct auxiliary information to the GNSS/INS integrated system.This model was developed to improve the land vehicle navigation accuracy in a blocked GNSS signal environment.
The main purpose of this research is to improve the smartphone's positioning accuracy for land vehicle navigation applications in environments where GPS accuracy has either deteriorated or its signal has been blocked.This can be accomplished by (1) developing an undifferenced and between-satellite single-difference (BSSD) PPP technique using single-frequency GNSS observations; (2) evaluating the model results in kinematic mode using a Xiaomi MI 8 smartphone and a high-end geodetic receiver; (3) providing an integrated positioning model based on combining the BSSD PPP results with the smartphone inertial sensors output; (4) developing a rectified road network mapping for the area of interest through a new map matching model based on hidden Marcov model (HMM) and Viterbi algorithms (Luo et al., 2017); (5) combining the positioning of the integrated system with a map matching model to obtain a restricted land vehicle navigation solution; and (6) assessment of the restricted integrated system in comparison with the reference RTK relative positioning solution.

BSSD PPP Technique for GPS and GLONASS
The standard single-frequency GNSS PPP technique of code and carrier phases, considering GPS time as a reference time system, can be expressed by the following (Abd Rabbou and El-Rabbany, 2015b) equations: To completely remove the receiver-related biases from both the pseudorange (P) and phase (∅) observations and intersystem biases, a satellite from each system is selected as reference.GPS (m) and GLONASS (n) satellites are taken as reference satellites to form the BSSD singlefrequency PPP technique (Abd Rabbou and El-Rabbany, 2015b): For equations (5-8), the term ∆ for GPS system (G) is the observation difference between any satellite and the reference satellite m, and for GLONASS system (L) is the observation difference of any satellite referenced to the satellite n, where: dt r -clock error for the receiver    ,    -clock errors for GPS and GLONASS satellites systems, respectively; ,  the true geometric range from receiver to satellite for GPS and GLONASS systems, respectively; d r , d s frequency-dependent code hardware delay for receiver and satellite, respectively; δ r , δ s frequency-dependent phase hardware delay for receiver and satellite, respectively; e  , e  code relevant system noise for GPS, GLONASS systems, respectively; ε  , ε  phase relevant system noise for GPS, GLONASS systems, respectively; T G , T Lionosphere errors for GPS and GLONASS systems, respectively; I G , I L troposphere errors for GPS and GLONASS systems, respectively; ISB L the intersystem bias for GLONASS satellites; N undifferenced ambiguity parameter.

INS mechanization
In the process of converting the output of an inertial measurement unit (IMU) into position, velocity, and attitude, the outputs include (  ,   ,   ) measured by gyroscopes and (  ,   ,   ) measured by accelerometers with respect to the body frame.The mechanization is executed in navigation frame (n-frame) using equations (9-12) (Noureldin et al., 2013): where:   = [  ℎ]  the position vector in the -frame, in latitude, longitude, and elevation;   = [      ] the velocity vector in the -frame, east, north, and up components; M the meridian radius of the ellipsoid and N the normal radius of the ellipsoid; ̇ the kinematic acceleration vector in the -frame; 2   the Coriolis acceleration vector and   is the gravity vector;    the skew-symmetric matrix of rotation rate vector of the earth (expressed in n-frame);    the skew-symmetric matrix of rotation rate vector of the n-frame with respect to e-frame;   the specific force vector in the -frame, measured by accelerometers;    the skew-symmetric matrix of the vector    , measured by gyroscopes;    the skew-symmetric matrix of the angular velocity vector    of the -frame w.r.t -frame.

Kalman filter algorithm
Recursive algorithm usually used in navigational application utilizes a set of equations to acquire an optimum estimate of the state vector which has the least variance and the Kalman Filter equations come under two steps: (1) prediction and (2) correction.The first step (equations 13 and 14) includes the prediction equations to predict the state vector and the associated covariance matrix of the system, depending on the current state (Abd Rabbou and El-Rabbany, 2015a): where: the superscript "^" indicates the estimated quantity; the superscript "-" represents a predicted quantity; the superscript "+" represents the corrected quantity.
The second step (equations 15 and 16) is to update the state vector and the associated covariance matrix with the available measurements regarding measurement model: refers to the Kalman gain at time   .Concerning the correction step, it depends on the Kalman gain, which is computed to produce a minimum error variance and is defined as in following equation: where: the state vector of the process at time   ;   the covariance matrix associated with state vector at epoch; Φ −1.thestate transition matrix from epoch  −1 to   ;   the system process noise covariance matrix;   the measurement (observation) vector, at time   ;   the observation matrix of the system at epoch   ;   the measurement noise covariance matrix.

Loosely coupled integration of BSSD with INS
The BSSD PPP output was loosely integrated with the INS using the Kalman filter, as shown in Figure 1.First, the raw IMU measurements obtained through mechanization equations using the quaternion method estimation were processed.Second, the resulting position and velocity were differenced from GPS measurement to enter the Kalman filter to get the error estimate of position and velocity.Finally, the resulting error estimate was used to update INS estimation to get the full state vector (position, velocity, and attitude).

Map Matching Technique
Map matching is a basic operation for improving the positioning accuracy by integrating GNSS tracking data with spatial road network data (roadway centerlines).It is used to identify the correct road link on which a vehicle is traveling and determine the location of a vehicle on the road network as shown in Figure 2 (Trough et al., 2020).Map matching algorithms can be either global or incremental.Global algorithms have batch processing for the entire input trajectory before generating the solution.Incremental algorithms employ localizing strategies that divide the input trajectory into smaller segments and process them sequentially without the knowledge of future inputs (Goh et al., 2012).
There are different types of map matching algorithms such as geometric, topological, and probabilistic.The probabilistic map matching requires applying the HMM algorithm, which is a recursive algorithm.It is used to find the unknown or hidden positions (candidates) on the road network with the help of the known measured positions (observations) of the trajectory.The mathematical model performs the following tasks sequentially to perform map matching:  The maximum probabilities are selected using the Viterbi algorithm as described later.

Calculation of the observation and transition probabilities
Each GNSS-INS point is assigned an observation probability value with candidates, as shown in Figure 3, which refers to how likely this candidate represents an observation on the road network.For the calculation of the observation probability, it is assumed that a position on the road network with certain distribution positions can be determined by GNSS.To simplify matters, it is postulated that the measurements from GNSS-INS point to the candidates are normally distributed (blue arrows in Figure 3).The observation probability is calculated with the following formula (Jung, 2019): where:   - . distance between the observation points and the candidate; σ z standard deviation.
Transition probability gives the probability of a vehicle moving between the road candidates.Some transitions will be very unlikely and require complicated maneuvers.Practically, transitions whose driving distance is about the same as the distance between the measurements were preferred.As shown in Figure 4, the P GNSS−INS −1 point was matched to the road network.
For matching P GNSS−INS  , the model forms transitions between the previous matched point and all selected candidates of the observation (C1, C2, and C3).A probability value is calculated for each transition.This probability is calculated using two terms.The first is based on the distance between GNSS-INS observations and the distance of the previous matched point to each of the selected candidates.The second is based on the orientation angles of these two distances, as shown in the following formula (Jung, 2019):

Application of the Viterbi algorithm and determination of the candidates that best represent the given trajectory
The HMM is applied in real-time mode as soon as the candidates have been determined and all probabilities have been calculated for each epoch.By applying the Viterbi algorithm to the HMM, the sequence of candidates is determined for which total probability is maximized.The next equations is to explain some of Viterbi algorithm basics.
For state space of the candidates of road links, n the number of observations and the time in seconds ranges from 1 to n.The objective is to find the candidates' sequence ( 1, ,   ), that has the maximum probability given the observations.Using the probability theory (Luo et al., 2017), we can conclude the following: In the HMM, the candidates at time t depend only on the state of the system as time  − 1 and the observation at time  depends only on the state of the system as time t.By substituting equations (25 and 26) in equation ( 24), the computable formula can be obtained as in equation ( 27).For each epoch, different paths from the previous matched point to all selected candidates are used to calculate the transition probabilities, as described in Figure 5.By merging the obtained candidates, the final matching path on the network is obtained and represents the most probable trajectory that, in this case, passes with candidates  1−1 ,  2−3 ,  3−2 , and  4−2 .

DATA COLLECTION
A real-life kinematic trajectory was provided to collect GNSS raw observations and inertial sensor data, where the two geodetic receivers of the Leica GS15 with the capability of tracking GNSS satellites were used to perform an RTK solution.The base receiver was at the point with known coordinates and the smartphone was held on a vehicle as shown in Figure 6.The kinematic trajectory's length was about 9.20 km, and it was performed in Madinaty city, Cairo.The layout plan of the trajectory is shown in Figure 7. GNSS data were observed for about 28 min.The reference positioning solution for the smartphone was the RTK relative positioning of the geodetic receivers.
Inertial sensor data and GNSS observation epochs were collected every second with 1674 used epochs.The RTK solution was also performed with 1 s interval, 100% fixed solutions were obtained, and the trajectory was observed as a loop and not a straight line.Epoch synchronization of different solutions was important for correlating and comparing the results to the reference RTK, where unsynchronized epochs were removed.At the beginning, the trajectory was in an open sky low-rise area (average height = 9 m); after 4 km, the trajectory led through high-rise areas (average height = 18 m) and the streets were so wide (up to eight lanes and traffic islands).
There are many Android applications capable of logging raw GNSS and sensor measurements.During the experiment, Geo++ RINEX Logger was used, which was released in 2017 by the Geo++ company, because it provides observation measurements for GPS, Galileo, and GLONASS systems, records the raw data in RINEX 3.03 format, allows processing by different software, and is available online for free.Also, "Androsensor" application was used to collect smartphone MEMS data to be utilized in INS mechanization.

METHODOLOGY
In the kinematic trajectory, the coordinate system used was the Universal Transverse Mercator (UTM) for both the smartphone and the geodetic rover receiver.The RINEX files were processed using the BSSD PPP technique, which made use of GPS/GLONASS satellite systems as well as orbital and clock corrections from the International GNSS Service (IGS) stations.INS data was mechanized by the quaternion method and the results were loosely integrated with BSSD outcome using the Kalman filter.Trajectory road network was downloaded from the open street map.The trajectory roads were buffered by their width.the points outside the buffer were matched and returned to the network centerlines by a probabilistic map matching technique, which utilized the HMM and Viterbi algorithm, as described before, at each epoch.As shown in Figure 8, three navigational solutions were obtained: BSSD PPP results (BSSD solution), the result of GNSS and INS integration (BSSD-INS solution), and the result of BSSD-INS-map matching (BSSD-INS-MM solution).For using the map matching results as a constraint in the navigation solution and simulating the signal loss problems, six outages along the trajectory were added.As shown in Figure 6, the outages lasted 10, 20, and 30 s for outages 1-2, 3-4, and 5-6, respectively.The accuracy of the resulting navigation solution was represented for the trajectory with and without added simulated outages.

ANALYSIS AND RESULTS
The results of the navigation solutions were correlated and plotted with reference to the relative RTK positioning using the UTM coordinate system.The results showed that the proposed BSSD-INS-map matching integration enhanced the navigation solution, irrespective of whether simulated outages were added or not.The RMSE for easting and northing positioning accuracy when applying BSSD-INS-based smartphone MEMS integration improved by 9%, compared to applying BSSD.The RMSE for easting and northing position accuracy when applying BSSD-INS-map matching integration improved by 30% and 28%, respectively, compared to the BSSD solution, in case of no outages inserted.The mean values for improving the proposed BSSD-INS-map matching integration were 91%, 96%, and 98% in 2D horizontal position accuracy compared to the BSSD-INS algorithm for six GNSS added signal outages with 10, 20, and 30 s periods, respectively.

Figure 1 .
Figure 1.Block diagram of loosely coupled GNSS-INS integration process

Figure 2 .
Figure 2. Map matching process illustration (left) and focusing at specific position (right) candidate point is a possible position of an observation on the network within a search distance.The road network is downloaded using open street maps (Export, 2020) and divided into equal line segments according to the data sampling rate and the speed of the vehicle.The beginning and endpoints of all road segments constitute a grid of points on the network.The calculation of candidates affects the computation time.A higher search distance can result in more candidates per observation and more time.

Figure 3 .
Figure 3. Selection of the candidate points of the trajectory between GNSS-INS points of successive epochs; D c the distance between the previous matched point and each of selected candidates; A  the orientation angle of GNSS-INS points' distance; A D c the orientation angle of the D c candidate distance.

Figure 5 .
Figure 5. Determination of the final matching trajectory

Figure 8 .
Figure 8. Block diagram of GPS-INS-map matching navigation solution Figures 9 and 10 represent easting and northing direction solutions of BSSD, BSSD-INS integration, and BSSD-INS-map matching integration tracks, respectively.The two charts are clarified with focused charts at 1 min along the track.The root mean square error (RMSE) and the maximum error for the three navigation solutions are shown in Figure 11.The RMSE improved by about 9% in easting, northing, and 2D horizontal when applying BSSD-INS compared to the BSSD navigation solution.The RMSE improved by about 30%, 28%, and 29% in easting, northing, and 2D horizontal, respectively, when applying BSSD-INS-MM compared to BSSD, and it also improved by about 22%, 21%, and 22% in easting, northing, and 2D horizontal, respectively, when applying BSSD-INS-MM compared to BSSD-INS.Figure 12 shows the navigation positioning solutions with the relative positioning RTK track for a part of the layout trajectory (the left) and at outage No. 1 layout (the right).

Table 1 .
Maximum error for BSSD-INS and BSSD-INS-MM tracks at the inserted outages