# Enhanced Map-Aided GPS/3D RISS Combined Positioning Strategy in Urban Canyons

Authors: Xiang Song; Chunxiao Ren; Huilin Jiang; Liping Li; Wei Wu; Ling Li; Shun Yan; Bingyu Zhang; Jiaen Wu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

To realize the effective positioning in urban canyons, an enhanced map-aided Global Positioning System (GPS)/three-dimensional (3D) reduced inertial sensor system (RISS) tightly combined positioning strategy is proposed. First, the 3D RISS is only based on the built-in controller area network (CAN). CAN bus sensor without additional sensors is first constructed to lower the cost. Then, a simple but effective enhanced map is created to assist positioning. Based on the map, a Kalman filtering (KF) tightly coupled method is proposed to fuse the 3D RISS with GPS information and to achieve the preliminary positioning. In KF-based preliminary positioning method, a simply observation noise variance optimization algorithm based on 2D enhanced map is proposed to improve KF method. In this algorithm, the value of the observation noise variance matrix is determined only according to the building plane information which is contained in the enhanced map. Further, a multiweight map matching algorithm is proposed for optimizing the initial positioning results. In this algorithm, factors such as distance, direction, road network topology, and lane change are considered and applied to map matching to further increase the positioning performance and form the final positioning results. Finally, the effectiveness of the strategy is proved by field test. The results show that this method has better accuracy and reliability than the conventional method.


---

## Body

## 1. Introduction

Accurate location information is the premise for ground vehicles to complete guidance- and safety-related tasks such as route planning, self-driving, and so on [1]. The most widely used technology for vehicle positioning is the Global Navigation Satellite System (GNSS) [2]. As a typical GNSS used in vehicle, Global Positioning System (GPS) provides all-weather and globally referenced positioning in open environments [3]. However, in sheltered environments such as urban canyons, the satellite signals are often blocked by tall buildings or thick tree cover and make GPS impractical [4]. In contrast, an inertial navigation system (INS) [5] based on a full inertial measurement unit (IMU) can realize high-frequency active positioning and has little dependence on the external environment. However, the positioning performance of INS will decline over time due to error accumulation [6]. Therefore, combining the respective advantages of GPS and INS and integrating them to form GPS/INS is an effective vehicle positioning solution [7].Although micro-electro-mechanical system (MEMS)-based IMU greatly reduces the cost of INS, its price is still expensive for ordinary vehicles. To further lower the cost, a reduced inertial sensor system (RISS) [8] involving one single-axis gyroscope, and the vehicle odometer was proposed. However, RISS only provides a two-dimensional (2D) positioning solution. To obtain three-dimensional (3D) solution, two additional accelerometers in the forward direction and transverse direction of the vehicle are used to determine the pitch angle and roll angle [9]. The main advantage of 3D RISS [10] is that it can obtain a full navigation solution with fewer sensors. In addition, the vehicle’s forward speed comes from the speedometer or odometer in 3D RISS, which reduces the accumulated error of the dual integration process.As an inertial system, 3D RISS still has the same drawbacks as INS especially the use of low-cost inertial sensors [11]. During GPS interruption, GPS/3D RISS integrated system will still cause significant drift of positioning solution. Especially in urban canyons, due to the satellite occlusion or reflection caused by tall building, the performance of GPS/3D RISS integrated positioning system decreases dramatically. Therefore, other auxiliary sensors are used to suppress this drift and improve positioning performance, such as motion constraints [12], speed sensors [13], camera [14], LiDAR sensors [15], and so on. However, most of these methods have their limitations. For example, the camera is easily affected by weather and light, and its real-time performance and reliability are poor; LiDAR is expensive and cannot be widely used in a short period of time.Recently, several 3D maps have been constructed and applied for map-matching or map-aiding in urban canyon [16–18]. The map-aiding technology mainly uses the height of buildings and other information to identify blocked satellites, which require accurate and rich 3D building information.The combined localization technology based on 3D map has effectively improved the localization performance with low cost and has attracted a lot of attention from researchers. However, the performance of these methods strongly depends on the accuracy of the 3D map. The current 3D-map making methods are difficult to achieve a balance between accuracy and storage efficiency. For example, in [19], the polygonal line is used to represent the road, and the road is simplified to the center line of the road. This method is simple and convenient, but the map accuracy is not high, and the data are not sufficient. The method proposed in [20] can achieve lane level accuracy but take up too much data storage space, which is difficult to meet the requirements of digital map storage efficiency. In addition, some methods use more expensive sensors, but it is difficult to achieve on the premise of low-cost [21].Overall, the 3D maps with high accuracy and rich information are not easily accessible and very expensive. For example, information such as the height of 3D buildings must be obtained from the local Planning Bureau. At the same time, 3D maps have large storage capacity due to complex structure, which is difficult to be suitable for real vehicle applications.Meanwhile, although RISS greatly reduces the number of sensors and thus reduces costs, automotive companies are very sensitive to hardware cost, and they more prefer to further reduce costs. On the other hand, as more vehicles are equipped with safety control systems [22], it can provide wheel speed, acceleration and other information via the vehicle’s controller area network (CAN) bus without any additional cost. In other words, all the sensors needed to build 3D RISS can be obtained from CAN bus without any additional cost, which gives us a possible solution.To address above problems, we propose an enhanced map-aided GPS/3D RISS tightly combined positioning strategy to improve positioning accuracy in urban canyons. In the proposed strategy, the 3D RISS only based on the built-in CAN bus sensor is firstly constructed. Then, a simple but effective enhanced map is created. Based on the map, an Kalman filtering (KF) method with the adaptive adjustment of observation noise variance is proposed to fuse the 3D RISS with GPS information and to achieve the preliminary positioning. Finally, a multiweight map matching method is proposed to correct the initial positioning results and obtain the final positioning results. The main contributions can be listed as follows:(1)
The 3D RISS relies only on the on-board CAN bus sensors for implementation without additional sensors, which significantly reduces the cost of in-vehicle applications.(2)
Based on the previous studies of our research group [23–26], an enhanced map containing rich information such as intersection details was designed and applied to assist fusion positioning. It should be noted that this map is a 2D map and only includes building plan information, which can be easily obtained from Open Street Map (OSM), so as to effectively reduce data volume while ensuring precision.(3)
In KF-based preliminary positioning method, a simply observation noise variance optimization algorithm based on 2D enhanced map is proposed to improve KF method. In this algorithm, the value of the observation noise variance matrix is determined only according to the building plane information which contained in the enhanced map without 3D building information. The algorithm solves the problem that the observation noise of each satellite may not be the same due to building occlusion, thus selecting available satellites, eliminating signals from satellites with large errors, and reducing the amount of computation. Compared with the map-aiding technology based on 3D map, it greatly reduces the difficulty of map accessible and the map-aiding complexity of the algorithm.(4)
A multiweight map matching (MWMM) algorithm is proposed for optimizing the initial positioning results. In this algorithm, factors such as distance, direction, road network topology, and lane change are considered and applied to map matching to further improve the positioning accuracy and form the final positioning results.The remainder of the article proceeds as follows. The overview of proposed strategy is given in Section2. Section 3 presents the CAN bus-based 3D RISS construction method and the conventional KF–based GPS/3D RISS combined positioning method. Section 4 is devoted to the enhanced map building method and 2D map-based observation noise variance adjustment algorithm of KF and multi-weight map matching algorithm. Section 5 shows the experimental results and analysis. Section 6 gives the conclusion.


## 2. Overview of the Proposed Positioning Strategy

The proposed strategy is shown as Figure1. There are three main components, that is, the environment perception part, the data processing part and information fusion part.Figure 1
Proposed vehicle positioning strategy.The environment perception part includes low-cost GPS receiver, in-vehicle sensors, and high-precision integrated navigation and positioning system. The in-vehicle sensors contain a wheel speed sensor for obtaining the vehicle’s travel speed, two accelerometers for obtaining transverse and forward acceleration, and a gyroscope for obtaining yaw rate. All these on-board sensors are included in the CAN bus without additional sensors. The low-cost GPS is used to obtain pseudorange with noise. It should be noted that the high-precision integrated navigation and positioning system is only used to obtain accurate vehicle trajectory information to make enhanced digital map and is not used in the process of fusion positioning.In the data processing part, the information from CAN bus information is processed to obtain the vehicle attitude and speed, and then the vehicle position estimation is obtained through the derivation of 3D RISS motion equations; the pseudorange, satellite ephemeris, ionospheric error, and other information output by the low-cost GPS receiver are solved to obtain the corrected pseudorange information; the accurate vehicle trajectory data obtained by the high-precision integrated navigation and positioning system is used for road modeling, intersection modeling and data screening to generate an enhanced map.The information fusion part is the key of the whole positioning strategy, which is divided into two parts: preliminary fusion positioning and final positioning. First, KF method is designed to fuse GPS and 3D RISS. Then, a simply observation noise variance optimization algorithm based on 2D enhanced map is proposed to improve KF method to obtain preliminary location estimation. In this algorithm, the value of the observation noise variance matrix can adaptive adjustment according to the building plane information which contained in the enhanced map. Finally, a multiweight map matching algorithm is proposed to further optimize the preliminary positioning results and output the final vehicle positioning results.


## 3. Conventional KF-Based GPS/3D RISS Tightly Combined Positioning Method

### 3.1. 3D Reduced Inertial Sensor System Based on CAN Bus Sensors

Two accelerometers mounted on the forward and transverse directions of the vehicle frame and a vertically aligned gyroscope in addition to the wheel speed sensors [11] are used in 3D RISS. In this article, all measurement data for 3D RISS are provided by CAN bus without additional sensors, which significantly reduces the cost.In this study, the state vector of RISS isX=ϕλhvEvNvUprA′, where φ, λ, and h are longitude, latitude, and altitude of vehicle; vE, vN, and vU represent the velocity components along the east, north, and up directions; indicate the pitch angle; r indicate the roll angle; A indicate the azimuth angle. The input vector of RISS is U=vfaffxfyωz′, where vf is the speed obtained from the wheel speed sensors; af is the acceleration derived from the wheel speed sensors; fx is the transversal accelerometer measurement; fy is the forward accelerometer reading; ωz is the angular rate obtained from the vertically aligned gyroscope. The nonlinear motion model for a 3D RISS involving the position, velocity, and attitude states is presented below.


#### 3.1.1. Pitch and Roll Calculation

The pitch angle is calculated as follows:(1)p=sin−1fy−afg,where g is the acceleration of gravity.Similarly, the roll angle is calculated as follows:(2)r=−sin−1fx+vfωzgcosp.


#### 3.1.2. Azimuth Calculation

The heading angle can be expressed as follows:(3)Ak=Ak−1−ωkzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt,where the subscript k indicates at time k, Δt is the sampling time, ωe denotes the magnitude of the Earth’s rotation rate, and RN is the normal radius.


#### 3.1.3. 3D Position and Velocity Calculations

According to the coordinate transformation relationship, the equation for converting the vehicle speed obtained by the wheel speed sensor from the body frame to the local-level frame is the following:(4)vkEvkNvkU=Rb,kl0vkf0=vkfsinAkcospkvkfcosAkcospkvkfsinpk,where Rb,kl is direction cosine matrix which transforms a vector from the body Frame to the Local-Level Frame.According to the north speedvkN, the latitude at time k can be expressed as follows:(5)φk=φk−1+dφdt|kΔt=φk−1+vkNRM+hkΔt=φk−1+vkfcosAkcospkRM+hkΔt,where RM is the meridian radius of curvature of the Earth.Similarly, according to the east speedvkE, the longitude is(6)λk=λk−1+dλdt|kΔt=λk−1+vkERN+hkvkUΔt=λk−1+vkfsinAkcospkRN+hkcosφkΔt,where RN is the normal radius of curvature of the Earth.According to the up speedvkU, the altitude is (7)hk=hk−1+dhdt|kΔt=hk−1+vkUΔt=hk−1+vkfsinpkΔt.The overall motion model is represented as follows:(8)Xk=φkλkhkvkEvkNvkUpkrkAk=fXk−1,Uk=φk−1+vkNRM+hkΔtλk−1+vkERN+hkcosφkΔthk−1+vkUΔtvkfsinAkcospkvkfcosAkcospkvkfsinpksin−1fky−akfg−sin−1fky+vkfωkzgcospkAk−1−ωzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt.


### 3.2. KF for Tightly Coupled 3D RISS/GPS Integration

To effectively reduce the influence of satellite signal occlusion on positioning results, the error of navigation parameters is used as the state vector, and GPS and 3D RISS are fused at the level of satellite pseudo range. It can be seen from the above that the motion model of 3D RISS is nonlinear, so it is linearized first to obtain the error model that can be used as the KF system model.It should be noted that this section only constructs the conventional GPS and 3D RISS fusion based on conventional KF, while the improvement of observation noise variance of KF filter method is designed in Section4.2. The preliminary positioning results of vehicles can be obtained by combining the methods in this section and Section 4.2.


#### 3.2.1. System State Model

The error-state system model [9] for KF tightly coupled positioning system can be written as follows:(9)δXk=Φk−1δXk−1+Gk−1ωk−1,where k is the discrete-time step; δXk indicates the state vector, Φk−1 is the state transition matrix, Gk−1 is the noise distribution matrix, and ωk−1 is the unit-variance white Gaussian noise which has a zero-mean and covariance Q. The error-state vector is δXk=δφk,δλk,δhk,δvkE,δvkN,δvkU,δAk,δFkw,δωkz,δbkGPS,δdkGPS′, where δφk, δλk, and δhk represent the position error in latitude, longitude, and altitude, respectively; δvkE, δvkN, and δvkU indicate the error in the velocity component along east direction, north direction, and up direction, respectively; δAk is the error in the azimuth angle, δFkw is the scale factor error of the wheel speed-derived speed, δωkz is the stochastic gyroscope drift, δbkGPS and δdkGPS are the bias of the GPS receiver’s clock and its drift.These equations are linearized by applying a Taylor series expansion and ignoring the higher-order terms. The corresponding linearized error-state system model is(10)δXk=δrkδvkδekδlk=I3×3F1O3×3O3×2O3×3I3×3F2O3×2O3×3O3×3F3O3×2O2×3O2×3O2×3F4δrk−1δvk−1δek−1δlk−1+O3×1O3×1σekσlk,where δrk=δφkδλkδhk; δvk=δvkEδvkNδvkU; δek=δAkδFkwδωkz; δlk=δbkGPSδdkGPS; σek=02gfσf2Δt2bzσz2Δt; σlk=σbΔtσdΔt;(11)F1=01RM+hk−1Δt01RN+hk−1cosφk−1Δt0000Δt;F2=vk−1fcosAk−1cospk−1vk−1fsinAk−1cospk−10−vk−1fsinAk−1cospk−1vk−1fcosAk−1cospk−100vk−1fsinpk−10;F3=10−Δt01−γfΔt0001−βzΔt;F4=1Δt01.The stochastic errors associated with the gyroscope and the wheel speed sensor are modeled by a Gauss–Markov model whereγf is the time-related error of wheel speed sensor, σf2 is the variance of the noise associated with it; βz is the reciprocal of the autocorrelation time for the gyroscope’s stochastic drift, and σz2 is the variance of the noise associated with it. σb and σd are the standard deviation of white noise for the clock bias and clock drift.


#### 3.2.2. Measurement Model

The linearized measurement model for KF [9] in tightly coupled integration is the following:(12)δZk=HkδXk+εk,where δZk is the measurement vector; Hk is the measurement design matrix; εk is a vector of measurement random noise which has a zero-mean and covariance R.In this article, the difference between the pseudo range of GPS and 3D RISS is used as the observation:(13)δZk=δZkρ,1δZkρ,2⋮δZkρ,m−1δZkρ,m=ρkRISS,1−ρkGPS,1ρkRISS,2−ρkGPS,2⋮ρkRISS,m−1−ρkGPS,m−1ρkRISS,m−ρkGPS,m,where m is the number of visible satellites to the receive at time k; ρkRISS,m is the RISS estimated pseudorange between the mth satellite and the receiver; ρkGPS,m is the GPS measured pseudorange between the mth satellite and the receiver.The pseudorange for themth satellite can therefore be modeled by(14)ρkGPS,m=rk−rkm+δbkGPS+εkρ,m,where rk is the receiver’s position in earth-centered earth-fix (ECEF) rectangular coordinates; rkm is position of the mth satellite in ECEF rectangular coordinates; δbkGPS is the receiver’s clock bias; εkρ,m is the total effect of residual errors due to atmospheric delays (after correcting for these errors using models inside the receiver itself), receiver noise, etc.The estimated pseudorange from the output of the RISS navigation system is defined as follows:(15)ρkRISS,m=rkRISS−rkm,where rkRISS is the position of the vehicle calculated from the position output of the RISS mechanization originally in ECEF geodetic coordinates and transformed to ECEF rectangular coordinates.The measurement design matrixHk=GkLk, and the calculation formula of Gk and Lk are [9, 13]:(16)Gk=lk1lk2⋮lkm−1lkm=lx,k1ly,k1lz,k1lx,k2ly,k2lz,k2⋮⋮⋮lx,km−1ly,km−1lz,km−1lx,kmly,kmlz,km,Lk=−RN+hk−1sinφk−1cosλk−1−RN+hk−1cosφk−1sinλk−1cosφk−1cosλk−1−RN+hk−1sinφk−1sinλk−1RN+hk−1cosφk−1cosλk−1cosφk−1sinλk−1RN1−e2+hk−1cosλk−10sinφk−1,where lkm is the line of sight unit vector from the mth satellite to the position of the receiver based on the output of RISS mechanization.


#### 3.2.3. Conventional KF for Fusion

Since the state equation and measurement equation have been linearized, the conventional KF can be used to fuse GPS and 3D RISS. The KF process consists of the following two phases:Time update:(17)δX^k,k−1=Φk,k−1δX^k−1,Pk,k−1=Φk,k−1Pk−1Φk,k−1′+Qk−1.Measurement update:(18)Kk=Pk,k−1HkHkPk,k−1Hk′+Rk−1,δX^k=δX^k,k−1+KkδZk−HkδX^k,k−1,Pk=I−KkHkPk,k−1,where I is an identity matrix.In general, the observation noise variance matrixR of conventional KF usually set as a fixed value. However, the observation noise of each satellite may not be the same under actual conditions. To solve this problem, an observation noise variance optimization algorithm based on 2D enhanced map will be designed in Section 4.2. By using the proposed KF method with the adaptive adjustment of observation noise variance, GPS and 3D RISS are fused to realize preliminary positioning.


## 3.1. 3D Reduced Inertial Sensor System Based on CAN Bus Sensors

Two accelerometers mounted on the forward and transverse directions of the vehicle frame and a vertically aligned gyroscope in addition to the wheel speed sensors [11] are used in 3D RISS. In this article, all measurement data for 3D RISS are provided by CAN bus without additional sensors, which significantly reduces the cost.In this study, the state vector of RISS isX=ϕλhvEvNvUprA′, where φ, λ, and h are longitude, latitude, and altitude of vehicle; vE, vN, and vU represent the velocity components along the east, north, and up directions; indicate the pitch angle; r indicate the roll angle; A indicate the azimuth angle. The input vector of RISS is U=vfaffxfyωz′, where vf is the speed obtained from the wheel speed sensors; af is the acceleration derived from the wheel speed sensors; fx is the transversal accelerometer measurement; fy is the forward accelerometer reading; ωz is the angular rate obtained from the vertically aligned gyroscope. The nonlinear motion model for a 3D RISS involving the position, velocity, and attitude states is presented below.


### 3.1.1. Pitch and Roll Calculation

The pitch angle is calculated as follows:(1)p=sin−1fy−afg,where g is the acceleration of gravity.Similarly, the roll angle is calculated as follows:(2)r=−sin−1fx+vfωzgcosp.


### 3.1.2. Azimuth Calculation

The heading angle can be expressed as follows:(3)Ak=Ak−1−ωkzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt,where the subscript k indicates at time k, Δt is the sampling time, ωe denotes the magnitude of the Earth’s rotation rate, and RN is the normal radius.


### 3.1.3. 3D Position and Velocity Calculations

According to the coordinate transformation relationship, the equation for converting the vehicle speed obtained by the wheel speed sensor from the body frame to the local-level frame is the following:(4)vkEvkNvkU=Rb,kl0vkf0=vkfsinAkcospkvkfcosAkcospkvkfsinpk,where Rb,kl is direction cosine matrix which transforms a vector from the body Frame to the Local-Level Frame.According to the north speedvkN, the latitude at time k can be expressed as follows:(5)φk=φk−1+dφdt|kΔt=φk−1+vkNRM+hkΔt=φk−1+vkfcosAkcospkRM+hkΔt,where RM is the meridian radius of curvature of the Earth.Similarly, according to the east speedvkE, the longitude is(6)λk=λk−1+dλdt|kΔt=λk−1+vkERN+hkvkUΔt=λk−1+vkfsinAkcospkRN+hkcosφkΔt,where RN is the normal radius of curvature of the Earth.According to the up speedvkU, the altitude is (7)hk=hk−1+dhdt|kΔt=hk−1+vkUΔt=hk−1+vkfsinpkΔt.The overall motion model is represented as follows:(8)Xk=φkλkhkvkEvkNvkUpkrkAk=fXk−1,Uk=φk−1+vkNRM+hkΔtλk−1+vkERN+hkcosφkΔthk−1+vkUΔtvkfsinAkcospkvkfcosAkcospkvkfsinpksin−1fky−akfg−sin−1fky+vkfωkzgcospkAk−1−ωzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt.


## 3.1.1. Pitch and Roll Calculation

The pitch angle is calculated as follows:(1)p=sin−1fy−afg,where g is the acceleration of gravity.Similarly, the roll angle is calculated as follows:(2)r=−sin−1fx+vfωzgcosp.


## 3.1.2. Azimuth Calculation

The heading angle can be expressed as follows:(3)Ak=Ak−1−ωkzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt,where the subscript k indicates at time k, Δt is the sampling time, ωe denotes the magnitude of the Earth’s rotation rate, and RN is the normal radius.


## 3.1.3. 3D Position and Velocity Calculations

According to the coordinate transformation relationship, the equation for converting the vehicle speed obtained by the wheel speed sensor from the body frame to the local-level frame is the following:(4)vkEvkNvkU=Rb,kl0vkf0=vkfsinAkcospkvkfcosAkcospkvkfsinpk,where Rb,kl is direction cosine matrix which transforms a vector from the body Frame to the Local-Level Frame.According to the north speedvkN, the latitude at time k can be expressed as follows:(5)φk=φk−1+dφdt|kΔt=φk−1+vkNRM+hkΔt=φk−1+vkfcosAkcospkRM+hkΔt,where RM is the meridian radius of curvature of the Earth.Similarly, according to the east speedvkE, the longitude is(6)λk=λk−1+dλdt|kΔt=λk−1+vkERN+hkvkUΔt=λk−1+vkfsinAkcospkRN+hkcosφkΔt,where RN is the normal radius of curvature of the Earth.According to the up speedvkU, the altitude is (7)hk=hk−1+dhdt|kΔt=hk−1+vkUΔt=hk−1+vkfsinpkΔt.The overall motion model is represented as follows:(8)Xk=φkλkhkvkEvkNvkUpkrkAk=fXk−1,Uk=φk−1+vkNRM+hkΔtλk−1+vkERN+hkcosφkΔthk−1+vkUΔtvkfsinAkcospkvkfcosAkcospkvkfsinpksin−1fky−akfg−sin−1fky+vkfωkzgcospkAk−1−ωzΔt+ωesinφk−1Δt+vk−1Etanφk−1RN+hk−1Δt.


## 3.2. KF for Tightly Coupled 3D RISS/GPS Integration

To effectively reduce the influence of satellite signal occlusion on positioning results, the error of navigation parameters is used as the state vector, and GPS and 3D RISS are fused at the level of satellite pseudo range. It can be seen from the above that the motion model of 3D RISS is nonlinear, so it is linearized first to obtain the error model that can be used as the KF system model.It should be noted that this section only constructs the conventional GPS and 3D RISS fusion based on conventional KF, while the improvement of observation noise variance of KF filter method is designed in Section4.2. The preliminary positioning results of vehicles can be obtained by combining the methods in this section and Section 4.2.


### 3.2.1. System State Model

The error-state system model [9] for KF tightly coupled positioning system can be written as follows:(9)δXk=Φk−1δXk−1+Gk−1ωk−1,where k is the discrete-time step; δXk indicates the state vector, Φk−1 is the state transition matrix, Gk−1 is the noise distribution matrix, and ωk−1 is the unit-variance white Gaussian noise which has a zero-mean and covariance Q. The error-state vector is δXk=δφk,δλk,δhk,δvkE,δvkN,δvkU,δAk,δFkw,δωkz,δbkGPS,δdkGPS′, where δφk, δλk, and δhk represent the position error in latitude, longitude, and altitude, respectively; δvkE, δvkN, and δvkU indicate the error in the velocity component along east direction, north direction, and up direction, respectively; δAk is the error in the azimuth angle, δFkw is the scale factor error of the wheel speed-derived speed, δωkz is the stochastic gyroscope drift, δbkGPS and δdkGPS are the bias of the GPS receiver’s clock and its drift.These equations are linearized by applying a Taylor series expansion and ignoring the higher-order terms. The corresponding linearized error-state system model is(10)δXk=δrkδvkδekδlk=I3×3F1O3×3O3×2O3×3I3×3F2O3×2O3×3O3×3F3O3×2O2×3O2×3O2×3F4δrk−1δvk−1δek−1δlk−1+O3×1O3×1σekσlk,where δrk=δφkδλkδhk; δvk=δvkEδvkNδvkU; δek=δAkδFkwδωkz; δlk=δbkGPSδdkGPS; σek=02gfσf2Δt2bzσz2Δt; σlk=σbΔtσdΔt;(11)F1=01RM+hk−1Δt01RN+hk−1cosφk−1Δt0000Δt;F2=vk−1fcosAk−1cospk−1vk−1fsinAk−1cospk−10−vk−1fsinAk−1cospk−1vk−1fcosAk−1cospk−100vk−1fsinpk−10;F3=10−Δt01−γfΔt0001−βzΔt;F4=1Δt01.The stochastic errors associated with the gyroscope and the wheel speed sensor are modeled by a Gauss–Markov model whereγf is the time-related error of wheel speed sensor, σf2 is the variance of the noise associated with it; βz is the reciprocal of the autocorrelation time for the gyroscope’s stochastic drift, and σz2 is the variance of the noise associated with it. σb and σd are the standard deviation of white noise for the clock bias and clock drift.


### 3.2.2. Measurement Model

The linearized measurement model for KF [9] in tightly coupled integration is the following:(12)δZk=HkδXk+εk,where δZk is the measurement vector; Hk is the measurement design matrix; εk is a vector of measurement random noise which has a zero-mean and covariance R.In this article, the difference between the pseudo range of GPS and 3D RISS is used as the observation:(13)δZk=δZkρ,1δZkρ,2⋮δZkρ,m−1δZkρ,m=ρkRISS,1−ρkGPS,1ρkRISS,2−ρkGPS,2⋮ρkRISS,m−1−ρkGPS,m−1ρkRISS,m−ρkGPS,m,where m is the number of visible satellites to the receive at time k; ρkRISS,m is the RISS estimated pseudorange between the mth satellite and the receiver; ρkGPS,m is the GPS measured pseudorange between the mth satellite and the receiver.The pseudorange for themth satellite can therefore be modeled by(14)ρkGPS,m=rk−rkm+δbkGPS+εkρ,m,where rk is the receiver’s position in earth-centered earth-fix (ECEF) rectangular coordinates; rkm is position of the mth satellite in ECEF rectangular coordinates; δbkGPS is the receiver’s clock bias; εkρ,m is the total effect of residual errors due to atmospheric delays (after correcting for these errors using models inside the receiver itself), receiver noise, etc.The estimated pseudorange from the output of the RISS navigation system is defined as follows:(15)ρkRISS,m=rkRISS−rkm,where rkRISS is the position of the vehicle calculated from the position output of the RISS mechanization originally in ECEF geodetic coordinates and transformed to ECEF rectangular coordinates.The measurement design matrixHk=GkLk, and the calculation formula of Gk and Lk are [9, 13]:(16)Gk=lk1lk2⋮lkm−1lkm=lx,k1ly,k1lz,k1lx,k2ly,k2lz,k2⋮⋮⋮lx,km−1ly,km−1lz,km−1lx,kmly,kmlz,km,Lk=−RN+hk−1sinφk−1cosλk−1−RN+hk−1cosφk−1sinλk−1cosφk−1cosλk−1−RN+hk−1sinφk−1sinλk−1RN+hk−1cosφk−1cosλk−1cosφk−1sinλk−1RN1−e2+hk−1cosλk−10sinφk−1,where lkm is the line of sight unit vector from the mth satellite to the position of the receiver based on the output of RISS mechanization.


### 3.2.3. Conventional KF for Fusion

Since the state equation and measurement equation have been linearized, the conventional KF can be used to fuse GPS and 3D RISS. The KF process consists of the following two phases:Time update:(17)δX^k,k−1=Φk,k−1δX^k−1,Pk,k−1=Φk,k−1Pk−1Φk,k−1′+Qk−1.Measurement update:(18)Kk=Pk,k−1HkHkPk,k−1Hk′+Rk−1,δX^k=δX^k,k−1+KkδZk−HkδX^k,k−1,Pk=I−KkHkPk,k−1,where I is an identity matrix.In general, the observation noise variance matrixR of conventional KF usually set as a fixed value. However, the observation noise of each satellite may not be the same under actual conditions. To solve this problem, an observation noise variance optimization algorithm based on 2D enhanced map will be designed in Section 4.2. By using the proposed KF method with the adaptive adjustment of observation noise variance, GPS and 3D RISS are fused to realize preliminary positioning.


## 3.2.1. System State Model

The error-state system model [9] for KF tightly coupled positioning system can be written as follows:(9)δXk=Φk−1δXk−1+Gk−1ωk−1,where k is the discrete-time step; δXk indicates the state vector, Φk−1 is the state transition matrix, Gk−1 is the noise distribution matrix, and ωk−1 is the unit-variance white Gaussian noise which has a zero-mean and covariance Q. The error-state vector is δXk=δφk,δλk,δhk,δvkE,δvkN,δvkU,δAk,δFkw,δωkz,δbkGPS,δdkGPS′, where δφk, δλk, and δhk represent the position error in latitude, longitude, and altitude, respectively; δvkE, δvkN, and δvkU indicate the error in the velocity component along east direction, north direction, and up direction, respectively; δAk is the error in the azimuth angle, δFkw is the scale factor error of the wheel speed-derived speed, δωkz is the stochastic gyroscope drift, δbkGPS and δdkGPS are the bias of the GPS receiver’s clock and its drift.These equations are linearized by applying a Taylor series expansion and ignoring the higher-order terms. The corresponding linearized error-state system model is(10)δXk=δrkδvkδekδlk=I3×3F1O3×3O3×2O3×3I3×3F2O3×2O3×3O3×3F3O3×2O2×3O2×3O2×3F4δrk−1δvk−1δek−1δlk−1+O3×1O3×1σekσlk,where δrk=δφkδλkδhk; δvk=δvkEδvkNδvkU; δek=δAkδFkwδωkz; δlk=δbkGPSδdkGPS; σek=02gfσf2Δt2bzσz2Δt; σlk=σbΔtσdΔt;(11)F1=01RM+hk−1Δt01RN+hk−1cosφk−1Δt0000Δt;F2=vk−1fcosAk−1cospk−1vk−1fsinAk−1cospk−10−vk−1fsinAk−1cospk−1vk−1fcosAk−1cospk−100vk−1fsinpk−10;F3=10−Δt01−γfΔt0001−βzΔt;F4=1Δt01.The stochastic errors associated with the gyroscope and the wheel speed sensor are modeled by a Gauss–Markov model whereγf is the time-related error of wheel speed sensor, σf2 is the variance of the noise associated with it; βz is the reciprocal of the autocorrelation time for the gyroscope’s stochastic drift, and σz2 is the variance of the noise associated with it. σb and σd are the standard deviation of white noise for the clock bias and clock drift.


## 3.2.2. Measurement Model

The linearized measurement model for KF [9] in tightly coupled integration is the following:(12)δZk=HkδXk+εk,where δZk is the measurement vector; Hk is the measurement design matrix; εk is a vector of measurement random noise which has a zero-mean and covariance R.In this article, the difference between the pseudo range of GPS and 3D RISS is used as the observation:(13)δZk=δZkρ,1δZkρ,2⋮δZkρ,m−1δZkρ,m=ρkRISS,1−ρkGPS,1ρkRISS,2−ρkGPS,2⋮ρkRISS,m−1−ρkGPS,m−1ρkRISS,m−ρkGPS,m,where m is the number of visible satellites to the receive at time k; ρkRISS,m is the RISS estimated pseudorange between the mth satellite and the receiver; ρkGPS,m is the GPS measured pseudorange between the mth satellite and the receiver.The pseudorange for themth satellite can therefore be modeled by(14)ρkGPS,m=rk−rkm+δbkGPS+εkρ,m,where rk is the receiver’s position in earth-centered earth-fix (ECEF) rectangular coordinates; rkm is position of the mth satellite in ECEF rectangular coordinates; δbkGPS is the receiver’s clock bias; εkρ,m is the total effect of residual errors due to atmospheric delays (after correcting for these errors using models inside the receiver itself), receiver noise, etc.The estimated pseudorange from the output of the RISS navigation system is defined as follows:(15)ρkRISS,m=rkRISS−rkm,where rkRISS is the position of the vehicle calculated from the position output of the RISS mechanization originally in ECEF geodetic coordinates and transformed to ECEF rectangular coordinates.The measurement design matrixHk=GkLk, and the calculation formula of Gk and Lk are [9, 13]:(16)Gk=lk1lk2⋮lkm−1lkm=lx,k1ly,k1lz,k1lx,k2ly,k2lz,k2⋮⋮⋮lx,km−1ly,km−1lz,km−1lx,kmly,kmlz,km,Lk=−RN+hk−1sinφk−1cosλk−1−RN+hk−1cosφk−1sinλk−1cosφk−1cosλk−1−RN+hk−1sinφk−1sinλk−1RN+hk−1cosφk−1cosλk−1cosφk−1sinλk−1RN1−e2+hk−1cosλk−10sinφk−1,where lkm is the line of sight unit vector from the mth satellite to the position of the receiver based on the output of RISS mechanization.


## 3.2.3. Conventional KF for Fusion

Since the state equation and measurement equation have been linearized, the conventional KF can be used to fuse GPS and 3D RISS. The KF process consists of the following two phases:Time update:(17)δX^k,k−1=Φk,k−1δX^k−1,Pk,k−1=Φk,k−1Pk−1Φk,k−1′+Qk−1.Measurement update:(18)Kk=Pk,k−1HkHkPk,k−1Hk′+Rk−1,δX^k=δX^k,k−1+KkδZk−HkδX^k,k−1,Pk=I−KkHkPk,k−1,where I is an identity matrix.In general, the observation noise variance matrixR of conventional KF usually set as a fixed value. However, the observation noise of each satellite may not be the same under actual conditions. To solve this problem, an observation noise variance optimization algorithm based on 2D enhanced map will be designed in Section 4.2. By using the proposed KF method with the adaptive adjustment of observation noise variance, GPS and 3D RISS are fused to realize preliminary positioning.


## 4. Enhanced Map and Map-Aided Positioning Algorithms

In this section, the construction of enhanced map is introduced to assist in the implementation of two map-based positioning algorithms, that is, the observation noise variance optimization algorithm for improving KF and the MWMM algorithm.


### 4.1. The Enhanced Map

To improve the positioning accuracy, the enhanced map in this article is constructed. In the proposed enhanced map, not only the accuracy has been improved to the sub meter level but also rich attribute information such as road geographic information, lane geographic information, and the topological relationship between roads and lanes has been added. At the same time, considering that satellite signals are easily blocked by dense buildings in the urban environment, the building plane information on both sides of the road is also added to the enhanced map. To reduce the map capacity and acquisition difficulty, the 3D information such as height of buildings is not added to the map. Moreover, intersections and their nearby areas are high incidence areas of satellite occlusion due to their more complex traffic environment (turning, merging, etc.) and the clustering of surrounding tall buildings, so the detailed information about the intersections is also added to the enhanced map in this article.For the related content of enhanced map construction involved in this article, the research group has carried out a lot of related studies. For example, [22] proposed a digital map making method based on cardinal spline to realize an optimal balance between the efficiency and reliability requirements of road models, [23] proposed a modeling method for urban complex intersections, [20] proposed a simple 3D map modeling method including the height of roadside buildings, which can be directly applied to this research.Since the purpose of this article is to use enhanced maps combined with low-cost sensors for high-precision positioning of vehicles, the accuracy of digital maps is required to be able to reach sub-meter level, but it is not necessary to include road attribute information such as traffic signs and traffic lights. Therefore, only Novatel SPAN-CPT high-precision combined navigation and positioning system is used to obtain the precise latitude and longitude information of the target road section. The ArcGIS is used as an editing tool to generate enhanced maps in this article.The overall scheme of map construction in this paper is shown in Figure2.Figure 2
Overall scheme of map construction.The specific map production method can be found in literature [22, 23]. The key aspects are:(1)
In data preprocessing part, the raw data is segmented according to the real-world road network topology, and the possible outliers are eliminated.(2)
In the part of general road modeling and intersection modeling, the Cardinal spline is chosen to establish an initial road model, which is specified by a series of control points and tension parameters. It overcomes the shortcomings of B-spline, such as requiring more nodes, complicated calculation, not suitable for practical application, and local adjust-ability is not enough to adapt to the large change of urban road curvature.(3)
In general road modeling part, based on the initial road model, a gradual optimization algorithm as [22] is introduced to determine the final general road model. According to the degree of the change of road curvature, the reasonable control points and optimal tension parameters can be determined.(4)
In intersection modeling part, based on the initial road model, an intersection model optimization algorithm as [23] is introduced to determine the final intersection model. According to different traffic conditions (right turn, left turn and straight travel), the optimal control points and optimal tension parameters of special intersections can be determined.(5)
In the part of general road modeling and intersection modeling, the building plane information on both sides of the road is added to the road network model. Only the 2D geographic information system (GIS) data which can be extracted from Open Street Map (OSM) including the layouts and positions of every building on both sides of the roads, is required and imported to ArcGIS. In addition, in the actual urban traffic environment, many buildings on both sides of the target road section have a long “span” in the direction of the road section. Therefore, the coordinates of the two endpoints of the building plane geometry in the direction of the target road section are used as an effective evaluation of the building span and introduced into the additional layers.(6)
Because the enhanced map contains specific lane information, it has a richer road network structure than the ordinary digital map. At the same time, it also greatly increases the amount of map data, which affect the efficiency of map access and query. In data point screening part, the data points of general roads and intersections are screened respectively according to the direction information of the center line of roads or virtual lanes. Compared with the data points evenly distributed by traditional methods, the data volume of the map is effectively reduced while ensuring the accuracy of the enhanced map.


### 4.2. The Observation Noise Variance Optimization Algorithm Based on Enhanced Map

In general, the observation noise variance matrixR of KF observation model is usually set as a fixed value. However, the observation noise of each satellite may not be the same under actual conditions. In this article, the difference between GPS and 3D RISS satellite pseudorange is used as the observation measurement. Since the ionospheric error, tropospheric error, and clock error have been corrected by satellite navigation messages, the nonline of sight (NLOS) interference and multipath effect caused by the roadside tall buildings in the urban environment have become the main factors affecting the pseudorange accuracy.To minimize the influence of roadside buildings on the satellite signal accuracy, many researchers use 3D maps to identify blocked satellites. However, 3D maps are difficult to obtain and have large capacity, which limits the application of real vehicles. Some researchers have proposed to evaluate the pseudorange accuracy and then fuzzy determine the observation noise variance by comprehensively considering the factors such as satellite altitude angle, the angle between satellite azimuth angle and road direction based on 2D map. However, it still exits some problems:(1)
Due to the presence of different building heights and the effect of reflections from taller buildings, it is difficult to evaluate the pseudorange accuracy based on the unified altitude angle threshold.(2)
The accuracy of the satellite that deviates from the road direction is not necessarily poor (e.g., the scenario that tall buildings only on one side of the road, and low on the other side), and the satellite along the road direction may also be blocked (e.g., the scenario that the curved road in the urban canyon).To solve the above problems, a simply adaptive adjustment method for the observation noise varianceR based on 2D map is proposed, as shown in Figure 3. In this method, the value of the observation noise variance matrix R is determined only according to the building plane information which contained in the 2D enhanced map building geographic information layer.Step 1: The pointeri of the building sequence is initialized, i = 1.Step 2: The shortest distanceDisi from the i-th building to the vehicle is calculated, and Disi is compared with the distance threshold Disth to determine whether Disi is less than Disth. If Disi<Disth, it is considered that the i-th building will affect the positioning accuracy of the vehicle, and then Step 3 is executed. Otherwise, skip to Step 7. (The value of Disth is related to the height of the building, in this paper, Disth is taken as 5 m).Step 3: The azimuth interval (AMin, AMax) of the i-th building relative to the vehicle is calculated. That is to say, the satellite from direction AMin to direction AMax of the vehicle will be affected by the i-th building.Step 4: The pointerj of the satellite sequence is initialized, j = 1.Step 5: IfAMin<Aj<AMax, R=RMax; else if AR−15∘<Aj<AR+15∘, R=RMin; else R=RMid (where Aj is the azimuth angle of the j-th satellite; AR is the direction angle of road).Step 6: Ifj<SA (where SA is the total number of visible satellites at current time), perform Step 7; else let j = j + 1, and skip to Step 5.Step 7: Ifi<BU (where BU is the total number of buildings that may block the satellite at the current time), the process is terminated and R is determined; else let j = j+1, and skip to Step 2.Figure 3
The observation noise variance optimization algorithm.The proposed algorithm in this section is applied to the KF algorithm in Section3.2 to adaptively adjust the observation noise variance matrix, and the fusion of GPS and SD RISS is realized to output the preliminary positioning results.


### 4.3. The MWMM Algorithm Based on Enhanced Map

The preliminary combined positioning method based on GPS and 3D RISS can solve the vehicle positioning problem in most situations, but the positioning accuracy is limited in the complex urban traffic environment. In order to further improve the positioning accuracy, the enhanced map is used to optimize the preliminary positioning results.In this article, a MWMM algorithm is proposed. The basic premise of MWMM is that the vehicle must drive on the road. Assuming that vehicle is driving near the centerline of a lane at the current time, the lane where the vehicle is located at the current time can be estimated by comprehensively considering the factors such as distance, direction, road network topology, lane change, and so on. Then, the preliminary estimation of the vehicle position obtained by GPS/3D RISS fusion is projected onto the corresponding lane centerline, and the final precise positioning of the vehicle can be realized.The MWMM algorithm proposed in this study includes the following aspects.


#### 4.3.1. The Matching Weight Based on Distance

The roads in vector maps are formed by orderly connection of a series of points. The simplest matching method is to calculate the distance between the current position of the vehicle and each point, and then match the vehicle to the point with the minimum distance from it. However, there are two problems with this method: (1) as shown in Figure4(a), point Pa is closest to lane L1. However, due to the closest Euclidean distance between Pa and point Sa, Pa is wrongly matched to Lane L2; (2) as shown in Figure 4(b), Pt is the actual position of the vehicle. Although point Pb is matched to the correct lane L2, it causes a large longitudinal matching error due to point-to-point matching.Figure 4
Distance-based map matching.
(a)(b)(c)To solve the above problems, map matching can be performed by calculating the distance between the vehicle position and the centerline of each lane. However, the matching only according to the distance relation between pointPb and each lane will be greatly affected by outliers. Compared with the point-to-point and point-to-line map matching method, the line-to-line map matching method based on the current vehicle position and the previous driving track can effectively reduce the influence of abnormal points on the matching result and reduce the mismatching. However, the track points in the past can only be used for reference. The earlier the track points are, the lower the reference value of the current matching. If the distance of track points in the past period is calculated with the same weight, it will cause the matching delay problem as shown in Figure 4(c). In Figure 4(c), the vehicle changes lanes from point P4, but the matching result obtained by the conventional least square method takes several seconds to get the correct matching result. Therefore, the weighted least square method is used for lane estimation in this article.In this article, the current position of the vehicle and the track point of the previous 4 s (a point in each second) are used for weighted least squares estimation. For any lane, the matching weight based on distance is calculated as follows:(19)Wdis=WdisMaxdi=0,WdisMin+∑i=15aidL−di3dL0<di<dL,WdisMindi≥dL,where Wdis is the matching weight based on distance; WdisMax=110 represents the maximum value of Wdis; WdisMin=10 represents the minimum value of Wdis; dL is the lane width of the road; di indicates the distance between data point i and the lane; ai is the matching coefficient, which is respectively taken as follows: a1 = 100, a2 = 80, a3 = 60, a4 = 40, a5 = 20 in this article.


#### 4.3.2. The Matching Weight Based on Road Direction

Generally speaking, different lanes on the same road have the same direction, but intersections are more complex. The possible driving path of vehicles within the intersection range is generally represented by the virtual lane center line in the enhanced map in this article. The exit lane and the entry lane of the virtual lane center line are not the same, the direction of the lane is also very different, and many virtual lane center lines are very close or even cross, so the direction of the lane will play a great role in judging the location of the vehicle. The enhanced map in this study has direction information whether it is the virtual lane of intersection or the actual lane line of ordinary road, which can be used in map matching algorithm.In this study, the angleα between the vehicle driving direction and the alternative matching road direction is calculated to determine the matching weight of each lane. Under the premise of observing traffic rules, when the vehicle is driving on the road, the driving direction of the vehicle shall be consistent with the road direction or the angle α is very small. As shown in Figure 5, when the vehicle is driving on road 1, the angle between the vehicle driving direction and the directions of roads 1 and 3 should be 0° or very small, and the angle with roads 2 and 4 will be very large. If the vehicle travels to road 2 at the next time, the angle between the vehicle driving direction and road 2 will become 0° or very small, and the angle with roads 1, 3, and 4 will become larger. Therefore, the matching weight based on road direction is calculated as follows:(20)Wdir=WdirMaxα≤15o,WdirMinα>15o,where Wdir denotes the direction weight of lane; WdirMax is constant 50; WdirMin is constant 0.Figure 5
Road direction-based map matching.


#### 4.3.3. The Matching Weight Based on Road Network Topology

If a lane cannot be reached from the lane in which the vehicle is located at the previous time through lane changing and other normal driving behaviors, that is, the two lanes are not topological connection relationship, then for the lane, the matching weight based on road network topology structureWtop is 0; otherwise, Wtop is 25; in particular, for the lane where the vehicle is at the previous moment, Wtop is 50.As shown in Figure6, the topological relationship diagram describes the connectivity between roads. Assuming that the vehicle is driving on Road 1 at one time, limited by the actual driving speed, the vehicle may drive on roads 1, 2, 3, and 4 at the next time, and will not drive on roads 5, 6, and 7. Therefore, roads 1, 2, 3, and 4 are selected as alternative matching roads. Among them, it is more likely to drive on Road 1 at the next moment than on roads 2, 3 and 4. Based on the above analysis, we set the weight of road 1 as 50, the weight of roads 2, 3, and 4 as 25, and the weight of other roads as 0.Figure 6
Road topological relationship.


#### 4.3.4. The Matching Weight Based on Lane Change

The enhanced map used in this study contains the direction information of the road. In the process of vehicle lane changing as shown in Figure7, the vehicle heading will shift and then return to normal, but the direction of the road does not change. Therefore, the difference between the road direction and the vehicle heading angle will show a crest or trough due to vehicle lane change, as shown in Figure 8. In Figure 8, the x-axis represents the vehicle running time and the y-axis represents the direction angle. In this article, the matching weight Wcl based on lane change is calculated by using the above characteristics of the difference between the road direction and heading angle during the vehicle change.(21)Wcl=WclMaxGdif>12,WclMax28<Gdif<12,WclMinGdif<8,where Gdif is the absolute value of the difference between the road direction and heading angle; WclMax is the constant 50; WclMin is the constant 0.Figure 7
Vehicle lane change.Figure 8
Difference between road direction and heading angle when vehicles change lanes.


#### 4.3.5. The Overall Matching Weight

Combined with the above factors, for vehicles running on the road, the overall matching weight of each lane is(22)Wtotal=Wdis+Wdir+Wtop+Wcl.The lane with the largest weightWtotal is considered to be the lane where the vehicle is traveling. Then the preliminary estimation of the vehicle position is projected vertically to the center line of the corresponding lane, and the final vehicle positioning result is obtained.


## 4.1. The Enhanced Map

To improve the positioning accuracy, the enhanced map in this article is constructed. In the proposed enhanced map, not only the accuracy has been improved to the sub meter level but also rich attribute information such as road geographic information, lane geographic information, and the topological relationship between roads and lanes has been added. At the same time, considering that satellite signals are easily blocked by dense buildings in the urban environment, the building plane information on both sides of the road is also added to the enhanced map. To reduce the map capacity and acquisition difficulty, the 3D information such as height of buildings is not added to the map. Moreover, intersections and their nearby areas are high incidence areas of satellite occlusion due to their more complex traffic environment (turning, merging, etc.) and the clustering of surrounding tall buildings, so the detailed information about the intersections is also added to the enhanced map in this article.For the related content of enhanced map construction involved in this article, the research group has carried out a lot of related studies. For example, [22] proposed a digital map making method based on cardinal spline to realize an optimal balance between the efficiency and reliability requirements of road models, [23] proposed a modeling method for urban complex intersections, [20] proposed a simple 3D map modeling method including the height of roadside buildings, which can be directly applied to this research.Since the purpose of this article is to use enhanced maps combined with low-cost sensors for high-precision positioning of vehicles, the accuracy of digital maps is required to be able to reach sub-meter level, but it is not necessary to include road attribute information such as traffic signs and traffic lights. Therefore, only Novatel SPAN-CPT high-precision combined navigation and positioning system is used to obtain the precise latitude and longitude information of the target road section. The ArcGIS is used as an editing tool to generate enhanced maps in this article.The overall scheme of map construction in this paper is shown in Figure2.Figure 2
Overall scheme of map construction.The specific map production method can be found in literature [22, 23]. The key aspects are:(1)
In data preprocessing part, the raw data is segmented according to the real-world road network topology, and the possible outliers are eliminated.(2)
In the part of general road modeling and intersection modeling, the Cardinal spline is chosen to establish an initial road model, which is specified by a series of control points and tension parameters. It overcomes the shortcomings of B-spline, such as requiring more nodes, complicated calculation, not suitable for practical application, and local adjust-ability is not enough to adapt to the large change of urban road curvature.(3)
In general road modeling part, based on the initial road model, a gradual optimization algorithm as [22] is introduced to determine the final general road model. According to the degree of the change of road curvature, the reasonable control points and optimal tension parameters can be determined.(4)
In intersection modeling part, based on the initial road model, an intersection model optimization algorithm as [23] is introduced to determine the final intersection model. According to different traffic conditions (right turn, left turn and straight travel), the optimal control points and optimal tension parameters of special intersections can be determined.(5)
In the part of general road modeling and intersection modeling, the building plane information on both sides of the road is added to the road network model. Only the 2D geographic information system (GIS) data which can be extracted from Open Street Map (OSM) including the layouts and positions of every building on both sides of the roads, is required and imported to ArcGIS. In addition, in the actual urban traffic environment, many buildings on both sides of the target road section have a long “span” in the direction of the road section. Therefore, the coordinates of the two endpoints of the building plane geometry in the direction of the target road section are used as an effective evaluation of the building span and introduced into the additional layers.(6)
Because the enhanced map contains specific lane information, it has a richer road network structure than the ordinary digital map. At the same time, it also greatly increases the amount of map data, which affect the efficiency of map access and query. In data point screening part, the data points of general roads and intersections are screened respectively according to the direction information of the center line of roads or virtual lanes. Compared with the data points evenly distributed by traditional methods, the data volume of the map is effectively reduced while ensuring the accuracy of the enhanced map.


## 4.2. The Observation Noise Variance Optimization Algorithm Based on Enhanced Map

In general, the observation noise variance matrixR of KF observation model is usually set as a fixed value. However, the observation noise of each satellite may not be the same under actual conditions. In this article, the difference between GPS and 3D RISS satellite pseudorange is used as the observation measurement. Since the ionospheric error, tropospheric error, and clock error have been corrected by satellite navigation messages, the nonline of sight (NLOS) interference and multipath effect caused by the roadside tall buildings in the urban environment have become the main factors affecting the pseudorange accuracy.To minimize the influence of roadside buildings on the satellite signal accuracy, many researchers use 3D maps to identify blocked satellites. However, 3D maps are difficult to obtain and have large capacity, which limits the application of real vehicles. Some researchers have proposed to evaluate the pseudorange accuracy and then fuzzy determine the observation noise variance by comprehensively considering the factors such as satellite altitude angle, the angle between satellite azimuth angle and road direction based on 2D map. However, it still exits some problems:(1)
Due to the presence of different building heights and the effect of reflections from taller buildings, it is difficult to evaluate the pseudorange accuracy based on the unified altitude angle threshold.(2)
The accuracy of the satellite that deviates from the road direction is not necessarily poor (e.g., the scenario that tall buildings only on one side of the road, and low on the other side), and the satellite along the road direction may also be blocked (e.g., the scenario that the curved road in the urban canyon).To solve the above problems, a simply adaptive adjustment method for the observation noise varianceR based on 2D map is proposed, as shown in Figure 3. In this method, the value of the observation noise variance matrix R is determined only according to the building plane information which contained in the 2D enhanced map building geographic information layer.Step 1: The pointeri of the building sequence is initialized, i = 1.Step 2: The shortest distanceDisi from the i-th building to the vehicle is calculated, and Disi is compared with the distance threshold Disth to determine whether Disi is less than Disth. If Disi<Disth, it is considered that the i-th building will affect the positioning accuracy of the vehicle, and then Step 3 is executed. Otherwise, skip to Step 7. (The value of Disth is related to the height of the building, in this paper, Disth is taken as 5 m).Step 3: The azimuth interval (AMin, AMax) of the i-th building relative to the vehicle is calculated. That is to say, the satellite from direction AMin to direction AMax of the vehicle will be affected by the i-th building.Step 4: The pointerj of the satellite sequence is initialized, j = 1.Step 5: IfAMin<Aj<AMax, R=RMax; else if AR−15∘<Aj<AR+15∘, R=RMin; else R=RMid (where Aj is the azimuth angle of the j-th satellite; AR is the direction angle of road).Step 6: Ifj<SA (where SA is the total number of visible satellites at current time), perform Step 7; else let j = j + 1, and skip to Step 5.Step 7: Ifi<BU (where BU is the total number of buildings that may block the satellite at the current time), the process is terminated and R is determined; else let j = j+1, and skip to Step 2.Figure 3
The observation noise variance optimization algorithm.The proposed algorithm in this section is applied to the KF algorithm in Section3.2 to adaptively adjust the observation noise variance matrix, and the fusion of GPS and SD RISS is realized to output the preliminary positioning results.


## 4.3. The MWMM Algorithm Based on Enhanced Map

The preliminary combined positioning method based on GPS and 3D RISS can solve the vehicle positioning problem in most situations, but the positioning accuracy is limited in the complex urban traffic environment. In order to further improve the positioning accuracy, the enhanced map is used to optimize the preliminary positioning results.In this article, a MWMM algorithm is proposed. The basic premise of MWMM is that the vehicle must drive on the road. Assuming that vehicle is driving near the centerline of a lane at the current time, the lane where the vehicle is located at the current time can be estimated by comprehensively considering the factors such as distance, direction, road network topology, lane change, and so on. Then, the preliminary estimation of the vehicle position obtained by GPS/3D RISS fusion is projected onto the corresponding lane centerline, and the final precise positioning of the vehicle can be realized.The MWMM algorithm proposed in this study includes the following aspects.


### 4.3.1. The Matching Weight Based on Distance

The roads in vector maps are formed by orderly connection of a series of points. The simplest matching method is to calculate the distance between the current position of the vehicle and each point, and then match the vehicle to the point with the minimum distance from it. However, there are two problems with this method: (1) as shown in Figure4(a), point Pa is closest to lane L1. However, due to the closest Euclidean distance between Pa and point Sa, Pa is wrongly matched to Lane L2; (2) as shown in Figure 4(b), Pt is the actual position of the vehicle. Although point Pb is matched to the correct lane L2, it causes a large longitudinal matching error due to point-to-point matching.Figure 4
Distance-based map matching.
(a)(b)(c)To solve the above problems, map matching can be performed by calculating the distance between the vehicle position and the centerline of each lane. However, the matching only according to the distance relation between pointPb and each lane will be greatly affected by outliers. Compared with the point-to-point and point-to-line map matching method, the line-to-line map matching method based on the current vehicle position and the previous driving track can effectively reduce the influence of abnormal points on the matching result and reduce the mismatching. However, the track points in the past can only be used for reference. The earlier the track points are, the lower the reference value of the current matching. If the distance of track points in the past period is calculated with the same weight, it will cause the matching delay problem as shown in Figure 4(c). In Figure 4(c), the vehicle changes lanes from point P4, but the matching result obtained by the conventional least square method takes several seconds to get the correct matching result. Therefore, the weighted least square method is used for lane estimation in this article.In this article, the current position of the vehicle and the track point of the previous 4 s (a point in each second) are used for weighted least squares estimation. For any lane, the matching weight based on distance is calculated as follows:(19)Wdis=WdisMaxdi=0,WdisMin+∑i=15aidL−di3dL0<di<dL,WdisMindi≥dL,where Wdis is the matching weight based on distance; WdisMax=110 represents the maximum value of Wdis; WdisMin=10 represents the minimum value of Wdis; dL is the lane width of the road; di indicates the distance between data point i and the lane; ai is the matching coefficient, which is respectively taken as follows: a1 = 100, a2 = 80, a3 = 60, a4 = 40, a5 = 20 in this article.


### 4.3.2. The Matching Weight Based on Road Direction

Generally speaking, different lanes on the same road have the same direction, but intersections are more complex. The possible driving path of vehicles within the intersection range is generally represented by the virtual lane center line in the enhanced map in this article. The exit lane and the entry lane of the virtual lane center line are not the same, the direction of the lane is also very different, and many virtual lane center lines are very close or even cross, so the direction of the lane will play a great role in judging the location of the vehicle. The enhanced map in this study has direction information whether it is the virtual lane of intersection or the actual lane line of ordinary road, which can be used in map matching algorithm.In this study, the angleα between the vehicle driving direction and the alternative matching road direction is calculated to determine the matching weight of each lane. Under the premise of observing traffic rules, when the vehicle is driving on the road, the driving direction of the vehicle shall be consistent with the road direction or the angle α is very small. As shown in Figure 5, when the vehicle is driving on road 1, the angle between the vehicle driving direction and the directions of roads 1 and 3 should be 0° or very small, and the angle with roads 2 and 4 will be very large. If the vehicle travels to road 2 at the next time, the angle between the vehicle driving direction and road 2 will become 0° or very small, and the angle with roads 1, 3, and 4 will become larger. Therefore, the matching weight based on road direction is calculated as follows:(20)Wdir=WdirMaxα≤15o,WdirMinα>15o,where Wdir denotes the direction weight of lane; WdirMax is constant 50; WdirMin is constant 0.Figure 5
Road direction-based map matching.


### 4.3.3. The Matching Weight Based on Road Network Topology

If a lane cannot be reached from the lane in which the vehicle is located at the previous time through lane changing and other normal driving behaviors, that is, the two lanes are not topological connection relationship, then for the lane, the matching weight based on road network topology structureWtop is 0; otherwise, Wtop is 25; in particular, for the lane where the vehicle is at the previous moment, Wtop is 50.As shown in Figure6, the topological relationship diagram describes the connectivity between roads. Assuming that the vehicle is driving on Road 1 at one time, limited by the actual driving speed, the vehicle may drive on roads 1, 2, 3, and 4 at the next time, and will not drive on roads 5, 6, and 7. Therefore, roads 1, 2, 3, and 4 are selected as alternative matching roads. Among them, it is more likely to drive on Road 1 at the next moment than on roads 2, 3 and 4. Based on the above analysis, we set the weight of road 1 as 50, the weight of roads 2, 3, and 4 as 25, and the weight of other roads as 0.Figure 6
Road topological relationship.


### 4.3.4. The Matching Weight Based on Lane Change

The enhanced map used in this study contains the direction information of the road. In the process of vehicle lane changing as shown in Figure7, the vehicle heading will shift and then return to normal, but the direction of the road does not change. Therefore, the difference between the road direction and the vehicle heading angle will show a crest or trough due to vehicle lane change, as shown in Figure 8. In Figure 8, the x-axis represents the vehicle running time and the y-axis represents the direction angle. In this article, the matching weight Wcl based on lane change is calculated by using the above characteristics of the difference between the road direction and heading angle during the vehicle change.(21)Wcl=WclMaxGdif>12,WclMax28<Gdif<12,WclMinGdif<8,where Gdif is the absolute value of the difference between the road direction and heading angle; WclMax is the constant 50; WclMin is the constant 0.Figure 7
Vehicle lane change.Figure 8
Difference between road direction and heading angle when vehicles change lanes.


### 4.3.5. The Overall Matching Weight

Combined with the above factors, for vehicles running on the road, the overall matching weight of each lane is(22)Wtotal=Wdis+Wdir+Wtop+Wcl.The lane with the largest weightWtotal is considered to be the lane where the vehicle is traveling. Then the preliminary estimation of the vehicle position is projected vertically to the center line of the corresponding lane, and the final vehicle positioning result is obtained.


## 4.3.1. The Matching Weight Based on Distance

The roads in vector maps are formed by orderly connection of a series of points. The simplest matching method is to calculate the distance between the current position of the vehicle and each point, and then match the vehicle to the point with the minimum distance from it. However, there are two problems with this method: (1) as shown in Figure4(a), point Pa is closest to lane L1. However, due to the closest Euclidean distance between Pa and point Sa, Pa is wrongly matched to Lane L2; (2) as shown in Figure 4(b), Pt is the actual position of the vehicle. Although point Pb is matched to the correct lane L2, it causes a large longitudinal matching error due to point-to-point matching.Figure 4
Distance-based map matching.
(a)(b)(c)To solve the above problems, map matching can be performed by calculating the distance between the vehicle position and the centerline of each lane. However, the matching only according to the distance relation between pointPb and each lane will be greatly affected by outliers. Compared with the point-to-point and point-to-line map matching method, the line-to-line map matching method based on the current vehicle position and the previous driving track can effectively reduce the influence of abnormal points on the matching result and reduce the mismatching. However, the track points in the past can only be used for reference. The earlier the track points are, the lower the reference value of the current matching. If the distance of track points in the past period is calculated with the same weight, it will cause the matching delay problem as shown in Figure 4(c). In Figure 4(c), the vehicle changes lanes from point P4, but the matching result obtained by the conventional least square method takes several seconds to get the correct matching result. Therefore, the weighted least square method is used for lane estimation in this article.In this article, the current position of the vehicle and the track point of the previous 4 s (a point in each second) are used for weighted least squares estimation. For any lane, the matching weight based on distance is calculated as follows:(19)Wdis=WdisMaxdi=0,WdisMin+∑i=15aidL−di3dL0<di<dL,WdisMindi≥dL,where Wdis is the matching weight based on distance; WdisMax=110 represents the maximum value of Wdis; WdisMin=10 represents the minimum value of Wdis; dL is the lane width of the road; di indicates the distance between data point i and the lane; ai is the matching coefficient, which is respectively taken as follows: a1 = 100, a2 = 80, a3 = 60, a4 = 40, a5 = 20 in this article.


## 4.3.2. The Matching Weight Based on Road Direction

Generally speaking, different lanes on the same road have the same direction, but intersections are more complex. The possible driving path of vehicles within the intersection range is generally represented by the virtual lane center line in the enhanced map in this article. The exit lane and the entry lane of the virtual lane center line are not the same, the direction of the lane is also very different, and many virtual lane center lines are very close or even cross, so the direction of the lane will play a great role in judging the location of the vehicle. The enhanced map in this study has direction information whether it is the virtual lane of intersection or the actual lane line of ordinary road, which can be used in map matching algorithm.In this study, the angleα between the vehicle driving direction and the alternative matching road direction is calculated to determine the matching weight of each lane. Under the premise of observing traffic rules, when the vehicle is driving on the road, the driving direction of the vehicle shall be consistent with the road direction or the angle α is very small. As shown in Figure 5, when the vehicle is driving on road 1, the angle between the vehicle driving direction and the directions of roads 1 and 3 should be 0° or very small, and the angle with roads 2 and 4 will be very large. If the vehicle travels to road 2 at the next time, the angle between the vehicle driving direction and road 2 will become 0° or very small, and the angle with roads 1, 3, and 4 will become larger. Therefore, the matching weight based on road direction is calculated as follows:(20)Wdir=WdirMaxα≤15o,WdirMinα>15o,where Wdir denotes the direction weight of lane; WdirMax is constant 50; WdirMin is constant 0.Figure 5
Road direction-based map matching.


## 4.3.3. The Matching Weight Based on Road Network Topology

If a lane cannot be reached from the lane in which the vehicle is located at the previous time through lane changing and other normal driving behaviors, that is, the two lanes are not topological connection relationship, then for the lane, the matching weight based on road network topology structureWtop is 0; otherwise, Wtop is 25; in particular, for the lane where the vehicle is at the previous moment, Wtop is 50.As shown in Figure6, the topological relationship diagram describes the connectivity between roads. Assuming that the vehicle is driving on Road 1 at one time, limited by the actual driving speed, the vehicle may drive on roads 1, 2, 3, and 4 at the next time, and will not drive on roads 5, 6, and 7. Therefore, roads 1, 2, 3, and 4 are selected as alternative matching roads. Among them, it is more likely to drive on Road 1 at the next moment than on roads 2, 3 and 4. Based on the above analysis, we set the weight of road 1 as 50, the weight of roads 2, 3, and 4 as 25, and the weight of other roads as 0.Figure 6
Road topological relationship.


## 4.3.4. The Matching Weight Based on Lane Change

The enhanced map used in this study contains the direction information of the road. In the process of vehicle lane changing as shown in Figure7, the vehicle heading will shift and then return to normal, but the direction of the road does not change. Therefore, the difference between the road direction and the vehicle heading angle will show a crest or trough due to vehicle lane change, as shown in Figure 8. In Figure 8, the x-axis represents the vehicle running time and the y-axis represents the direction angle. In this article, the matching weight Wcl based on lane change is calculated by using the above characteristics of the difference between the road direction and heading angle during the vehicle change.(21)Wcl=WclMaxGdif>12,WclMax28<Gdif<12,WclMinGdif<8,where Gdif is the absolute value of the difference between the road direction and heading angle; WclMax is the constant 50; WclMin is the constant 0.Figure 7
Vehicle lane change.Figure 8
Difference between road direction and heading angle when vehicles change lanes.


## 4.3.5. The Overall Matching Weight

Combined with the above factors, for vehicles running on the road, the overall matching weight of each lane is(22)Wtotal=Wdis+Wdir+Wtop+Wcl.The lane with the largest weightWtotal is considered to be the lane where the vehicle is traveling. Then the preliminary estimation of the vehicle position is projected vertically to the center line of the corresponding lane, and the final vehicle positioning result is obtained.


## 5. Experiments and Results

To evaluate the proposed strategy, several experiments were conducted by using a Chery TIGGO5 sport utility vehicle on real road. The positioning sensors used in the vehicle only include low-cost GPS and Controller Area Network (CAN) bus data acquisition card. The low-cost GPS adopted the C230-AT provided by Beijing BDStar Navigation Company Ltd. The CAN bus data acquisition card adopted the PCIe-9221 provided by of Zhiyuan Electronics Co., Ltd., which is used to collect the on-board sensor information from the CAN bus. In addition, the high-precision integrated navigation system Novatel SPAN-CPT was used to make the enhanced digital map and also used as a reference to verify the accuracy of the proposed strategy.The area near Nanjing Olympic Sports Center was selected as the experimental site in this study, and its satellite aerial photos are shown in Figure9. In the test, the vehicle starts from the yellow sign and drives clockwise. The blue closed line segment represents the test route, with a total length of about 4.1 km and a time of 459 s. The test route includes a variety of typical sections such as straight lines, curves, and intersections. There are several tall buildings on the east route and dense trees on the west route.Figure 9
Experimental site.


### 5.1. The Enhanced Map

The enhanced digital map generated in this study is only used to optimize the vehicle positioning results on the test road and does not involve the surrounding geographical environment. Therefore, this article only focuses on the production of enhanced digital map of the test road and its related sections.


#### 5.1.1. Road Modeling

In this article, the original trajectory data points collected by the experimental vehicle are used for road modeling based on Cardinal spline. First, an initial road model is established by using Cardinal spline, and then the model is optimized to determine its reasonable control points and optimal tension parameters. Finally, the efficiency and reliability of this method and B-spline road modeling method are compared and analyzed. Figure10 shows road modeling process based on Cardinal spline.Figure 10
The road modeling process. (a) Original track data points; (b) results of the selection of control points; (c) results of tension parameter adjustment; (d) the road model based on Cardinal spline.
(a)(b)(c)(d)The 378 raw trajectory data points from the detection vehicles were used as input to the curvature analysis-based road model optimization algorithm. After determining reasonable control points and optimal tension parameters, an optimized road model was obtained, which contains 193 control points and 192 Cardinal spline curve segments. To evaluate the performance of the proposed road modeling method, we compared the method with the B-spline-based road modeling method proposed in the literature [20]. The B-spline based road model has been shown to outperform various previous road modeling methods such as those using polygons, gyratory curves, and partial spline curves. The detailed comparison results are shown in Tables 1 and 2.Table 1
Comparison of the efficiency of two road modeling methods.
B-splineCardinal splineReliabilityHighHighGlobal position error RMS (m)1.251.15EfficiencyMiddleHighData storage1142578FlexibilityMiddleHighTable 2
Comparison of the reliability of two road modeling methods.
B-splineCardinal splineReliabilityLowHighGlobal position error RMS (m)2.831.15EfficiencyHighHighData storage578578FlexibilityMiddleHighFrom the above comparison results, it can be seen that the road modeling method based on B-spline uses 1142 floating-point numbers to achieve the global positioning accuracy of 1.25 m and 578 floating-point numbers to achieve the global positioning accuracy of 2.83 m. At the same time, the proposed road modeling method based on Cardinal spline uses 578 floating-point numbers to achieve a global accuracy of 1.15 m. Figure11(a) shows the efficiency comparison of the two road modeling methods. The results show that compared with the road modeling method based on B-spline, the road modeling method based on Cardinal spline uses fewer data points and can meet the desired global accuracy. Figure 11(b) shows the reliability comparison of the two road modeling methods. The results show that when the same number of data points is used, the road modeling method based on Cardinal spline has higher global accuracy than the road modeling method based on B-spline. To sum up, the road modeling method based on Cardinal spline is better than the road modeling method based on B-spline, that is, this method can better achieve the balance between efficiency and reliability.Figure 11
Comparison diagram of two road modeling methods. (a) Sematic diagram of efficiency comparison; (b) sematic diagram of reliability comparison.
(a)(b)


#### 5.1.2. Map Visualization

According to the overall production scheme of enhanced map proposed in Section4.1 of this article, the two-dimensional enhanced digital map of the experimental site is produced. Figure 12 shows the overall effect diagram of the enhanced map of the test site. (a) is the vector map, the red line segment represents the lane, the yellow area represents the intersection, and the light cyan block represents the plane structure of the building. (b) is the effect of the lane and intersection in the vector digital map in Google Earth.Figure 12
Overall effect diagram of the enhanced map of the test site.
(a)(b)Figure13 shows the enhanced map near the intersection. (a) is the vector map effect diagram, (b) is the vector map effect in Google Earth. In Figure 13, the red line is the lane center line of the conventional road, and the yellow line is the virtual lane center line of the intersection.Figure 13
The enhanced map near the intersection.
(a)(b)


### 5.2. Performance of Fusion Positioning

According to the proposed enhanced map-based vehicle fusion positioning scheme, first, low-cost GPS and RISS are fused through KF with adaptive adjustment of observation noise variance to obtain the preliminary positioning results, and then the preliminary positioning results are further optimized and improved by using the produced enhanced map to obtain the final positioning results.In order to more intuitively show the experimental effect of the combined positioning scheme, the vehicle position is uniformly transformed from the geocentric geostationary coordinate system to the local tangent plane coordinate system. The positioning results of the two stages are shown in Figure14. The experiment takes 436 s (excluding the seconds when the vehicle is stationary at the beginning and end of the experiment). The red curve in the figure represents the reference value of the vehicle's horizontal position (output by the high-precision integrated navigation and positioning system), and the blue curve represents the preliminary positioning results of adaptive setting of observation noise, the cyan curve represents the final positioning result obtained by the multi-weight map matching algorithm. In Figure 14, three local areas are enlarged, including two straight roads and an intersection.Figure 14
Comparison diagram of positioning results.To separately analyze the effect of observation noise variance optimization algorithm and multiweight matching algorithm based on enhanced map. Figure15 shows the comparison between KF positioning results with fixed observation noise variance (KF1) and KF positioning results with adaptive observation noise variance setting (KF2), Figure 16 shows the comparison between KF positioning results with adaptive observation noise variance setting (KF2, i.e., the preliminary positioning) and multiweight map matching positioning results (MMWM, i.e., the final positioning), and Table 3 gives the gives the statistics of Euclidean distance errors (i.e., horizontal position errors) contain the maximum value and the root mean square (RMS) value of the three positioning methods.Figure 15
Comparison between KF1 and KF2.Figure 16
Comparison between KF2 and MWMM.Table 3
Comparison of three positioning methods.
MethodStatistics of Euclidean distance errors (unit: m)MaxRMSKF117.913.664KF212.563.114MWMM6.2402.499From Figures15 and 16 and Table 3, it can be seen that the KF1 has the worst accuracy, that is, both the RMS and the maximum values of its Euclidean distance error are the largest. It can be attributed to the blockage of satellites by roadside buildings or trees. The preliminary positioning (KF2) can achieve higher accuracy than KF1. This is because the adaptive adjustment of observation noise variance can alleviate the pollution of invalid satellite observation to positioning results. However, the max value of its Euclidean distance error is large. The main reason is that in the area with serious satellite occlusion, the number of visible satellites after removing invalid satellites is insufficient, resulting in serious positioning solution error.From Table3, it is obviously that the positioning accuracy of final fusion positioning (MWMM) is improved obviously when compared with KF1 and KF2. For instance, the RMS value of MWMM error is decreased to 2.499 m from the value 3.114 m of preliminary positioning and the value 3.664 m of KF1. The improvement of maximum error is more obvious, that is, achieve 65% accuracy improvement compared to KF1 and 50% improvement compared to conventional KF2. It can be attributed that the accurate maps and MWMM algorithms can effectively correct outliers.To further verify the effect of the proposed algorithm, the results of another experiment are shown in Table4 and Figure 17. The experiment took 400 seconds. The satellite occlusion on the experimental route was less than that in the previous experiment, and it was mainly covered by trees or buildings on one side. From Figure 17 and Table 4, it can be seen that the RMS of KF2 is lower than that of KF1, but it cannot alleviate the maximum error. The RMS of MWMM is only 2.236 m, which is much smaller than the first two cases, and has a significant inhibitory effect on the maximum deviation.Table 4
Comparison of three positioning methods in another experiment.
MethodStatistics of Euclidean distance errors (unit: m)MaxRMSKF17.9452.932KF27.5452.696MWMM5.6592.236Figure 17
Comparison of the three methods in another experiment.Many groups of experiments have been carried out, and the results show that the localization result of KF2 is better than that of KF1. The final positioning result obtained by MWMM is greatly improved compared with both KF1 and KF2, especially for the points with large error.


## 5.1. The Enhanced Map

The enhanced digital map generated in this study is only used to optimize the vehicle positioning results on the test road and does not involve the surrounding geographical environment. Therefore, this article only focuses on the production of enhanced digital map of the test road and its related sections.


### 5.1.1. Road Modeling

In this article, the original trajectory data points collected by the experimental vehicle are used for road modeling based on Cardinal spline. First, an initial road model is established by using Cardinal spline, and then the model is optimized to determine its reasonable control points and optimal tension parameters. Finally, the efficiency and reliability of this method and B-spline road modeling method are compared and analyzed. Figure10 shows road modeling process based on Cardinal spline.Figure 10
The road modeling process. (a) Original track data points; (b) results of the selection of control points; (c) results of tension parameter adjustment; (d) the road model based on Cardinal spline.
(a)(b)(c)(d)The 378 raw trajectory data points from the detection vehicles were used as input to the curvature analysis-based road model optimization algorithm. After determining reasonable control points and optimal tension parameters, an optimized road model was obtained, which contains 193 control points and 192 Cardinal spline curve segments. To evaluate the performance of the proposed road modeling method, we compared the method with the B-spline-based road modeling method proposed in the literature [20]. The B-spline based road model has been shown to outperform various previous road modeling methods such as those using polygons, gyratory curves, and partial spline curves. The detailed comparison results are shown in Tables 1 and 2.Table 1
Comparison of the efficiency of two road modeling methods.
B-splineCardinal splineReliabilityHighHighGlobal position error RMS (m)1.251.15EfficiencyMiddleHighData storage1142578FlexibilityMiddleHighTable 2
Comparison of the reliability of two road modeling methods.
B-splineCardinal splineReliabilityLowHighGlobal position error RMS (m)2.831.15EfficiencyHighHighData storage578578FlexibilityMiddleHighFrom the above comparison results, it can be seen that the road modeling method based on B-spline uses 1142 floating-point numbers to achieve the global positioning accuracy of 1.25 m and 578 floating-point numbers to achieve the global positioning accuracy of 2.83 m. At the same time, the proposed road modeling method based on Cardinal spline uses 578 floating-point numbers to achieve a global accuracy of 1.15 m. Figure11(a) shows the efficiency comparison of the two road modeling methods. The results show that compared with the road modeling method based on B-spline, the road modeling method based on Cardinal spline uses fewer data points and can meet the desired global accuracy. Figure 11(b) shows the reliability comparison of the two road modeling methods. The results show that when the same number of data points is used, the road modeling method based on Cardinal spline has higher global accuracy than the road modeling method based on B-spline. To sum up, the road modeling method based on Cardinal spline is better than the road modeling method based on B-spline, that is, this method can better achieve the balance between efficiency and reliability.Figure 11
Comparison diagram of two road modeling methods. (a) Sematic diagram of efficiency comparison; (b) sematic diagram of reliability comparison.
(a)(b)


### 5.1.2. Map Visualization

According to the overall production scheme of enhanced map proposed in Section4.1 of this article, the two-dimensional enhanced digital map of the experimental site is produced. Figure 12 shows the overall effect diagram of the enhanced map of the test site. (a) is the vector map, the red line segment represents the lane, the yellow area represents the intersection, and the light cyan block represents the plane structure of the building. (b) is the effect of the lane and intersection in the vector digital map in Google Earth.Figure 12
Overall effect diagram of the enhanced map of the test site.
(a)(b)Figure13 shows the enhanced map near the intersection. (a) is the vector map effect diagram, (b) is the vector map effect in Google Earth. In Figure 13, the red line is the lane center line of the conventional road, and the yellow line is the virtual lane center line of the intersection.Figure 13
The enhanced map near the intersection.
(a)(b)


## 5.1.1. Road Modeling

In this article, the original trajectory data points collected by the experimental vehicle are used for road modeling based on Cardinal spline. First, an initial road model is established by using Cardinal spline, and then the model is optimized to determine its reasonable control points and optimal tension parameters. Finally, the efficiency and reliability of this method and B-spline road modeling method are compared and analyzed. Figure10 shows road modeling process based on Cardinal spline.Figure 10
The road modeling process. (a) Original track data points; (b) results of the selection of control points; (c) results of tension parameter adjustment; (d) the road model based on Cardinal spline.
(a)(b)(c)(d)The 378 raw trajectory data points from the detection vehicles were used as input to the curvature analysis-based road model optimization algorithm. After determining reasonable control points and optimal tension parameters, an optimized road model was obtained, which contains 193 control points and 192 Cardinal spline curve segments. To evaluate the performance of the proposed road modeling method, we compared the method with the B-spline-based road modeling method proposed in the literature [20]. The B-spline based road model has been shown to outperform various previous road modeling methods such as those using polygons, gyratory curves, and partial spline curves. The detailed comparison results are shown in Tables 1 and 2.Table 1
Comparison of the efficiency of two road modeling methods.
B-splineCardinal splineReliabilityHighHighGlobal position error RMS (m)1.251.15EfficiencyMiddleHighData storage1142578FlexibilityMiddleHighTable 2
Comparison of the reliability of two road modeling methods.
B-splineCardinal splineReliabilityLowHighGlobal position error RMS (m)2.831.15EfficiencyHighHighData storage578578FlexibilityMiddleHighFrom the above comparison results, it can be seen that the road modeling method based on B-spline uses 1142 floating-point numbers to achieve the global positioning accuracy of 1.25 m and 578 floating-point numbers to achieve the global positioning accuracy of 2.83 m. At the same time, the proposed road modeling method based on Cardinal spline uses 578 floating-point numbers to achieve a global accuracy of 1.15 m. Figure11(a) shows the efficiency comparison of the two road modeling methods. The results show that compared with the road modeling method based on B-spline, the road modeling method based on Cardinal spline uses fewer data points and can meet the desired global accuracy. Figure 11(b) shows the reliability comparison of the two road modeling methods. The results show that when the same number of data points is used, the road modeling method based on Cardinal spline has higher global accuracy than the road modeling method based on B-spline. To sum up, the road modeling method based on Cardinal spline is better than the road modeling method based on B-spline, that is, this method can better achieve the balance between efficiency and reliability.Figure 11
Comparison diagram of two road modeling methods. (a) Sematic diagram of efficiency comparison; (b) sematic diagram of reliability comparison.
(a)(b)


## 5.1.2. Map Visualization

According to the overall production scheme of enhanced map proposed in Section4.1 of this article, the two-dimensional enhanced digital map of the experimental site is produced. Figure 12 shows the overall effect diagram of the enhanced map of the test site. (a) is the vector map, the red line segment represents the lane, the yellow area represents the intersection, and the light cyan block represents the plane structure of the building. (b) is the effect of the lane and intersection in the vector digital map in Google Earth.Figure 12
Overall effect diagram of the enhanced map of the test site.
(a)(b)Figure13 shows the enhanced map near the intersection. (a) is the vector map effect diagram, (b) is the vector map effect in Google Earth. In Figure 13, the red line is the lane center line of the conventional road, and the yellow line is the virtual lane center line of the intersection.Figure 13
The enhanced map near the intersection.
(a)(b)


## 5.2. Performance of Fusion Positioning

According to the proposed enhanced map-based vehicle fusion positioning scheme, first, low-cost GPS and RISS are fused through KF with adaptive adjustment of observation noise variance to obtain the preliminary positioning results, and then the preliminary positioning results are further optimized and improved by using the produced enhanced map to obtain the final positioning results.In order to more intuitively show the experimental effect of the combined positioning scheme, the vehicle position is uniformly transformed from the geocentric geostationary coordinate system to the local tangent plane coordinate system. The positioning results of the two stages are shown in Figure14. The experiment takes 436 s (excluding the seconds when the vehicle is stationary at the beginning and end of the experiment). The red curve in the figure represents the reference value of the vehicle's horizontal position (output by the high-precision integrated navigation and positioning system), and the blue curve represents the preliminary positioning results of adaptive setting of observation noise, the cyan curve represents the final positioning result obtained by the multi-weight map matching algorithm. In Figure 14, three local areas are enlarged, including two straight roads and an intersection.Figure 14
Comparison diagram of positioning results.To separately analyze the effect of observation noise variance optimization algorithm and multiweight matching algorithm based on enhanced map. Figure15 shows the comparison between KF positioning results with fixed observation noise variance (KF1) and KF positioning results with adaptive observation noise variance setting (KF2), Figure 16 shows the comparison between KF positioning results with adaptive observation noise variance setting (KF2, i.e., the preliminary positioning) and multiweight map matching positioning results (MMWM, i.e., the final positioning), and Table 3 gives the gives the statistics of Euclidean distance errors (i.e., horizontal position errors) contain the maximum value and the root mean square (RMS) value of the three positioning methods.Figure 15
Comparison between KF1 and KF2.Figure 16
Comparison between KF2 and MWMM.Table 3
Comparison of three positioning methods.
MethodStatistics of Euclidean distance errors (unit: m)MaxRMSKF117.913.664KF212.563.114MWMM6.2402.499From Figures15 and 16 and Table 3, it can be seen that the KF1 has the worst accuracy, that is, both the RMS and the maximum values of its Euclidean distance error are the largest. It can be attributed to the blockage of satellites by roadside buildings or trees. The preliminary positioning (KF2) can achieve higher accuracy than KF1. This is because the adaptive adjustment of observation noise variance can alleviate the pollution of invalid satellite observation to positioning results. However, the max value of its Euclidean distance error is large. The main reason is that in the area with serious satellite occlusion, the number of visible satellites after removing invalid satellites is insufficient, resulting in serious positioning solution error.From Table3, it is obviously that the positioning accuracy of final fusion positioning (MWMM) is improved obviously when compared with KF1 and KF2. For instance, the RMS value of MWMM error is decreased to 2.499 m from the value 3.114 m of preliminary positioning and the value 3.664 m of KF1. The improvement of maximum error is more obvious, that is, achieve 65% accuracy improvement compared to KF1 and 50% improvement compared to conventional KF2. It can be attributed that the accurate maps and MWMM algorithms can effectively correct outliers.To further verify the effect of the proposed algorithm, the results of another experiment are shown in Table4 and Figure 17. The experiment took 400 seconds. The satellite occlusion on the experimental route was less than that in the previous experiment, and it was mainly covered by trees or buildings on one side. From Figure 17 and Table 4, it can be seen that the RMS of KF2 is lower than that of KF1, but it cannot alleviate the maximum error. The RMS of MWMM is only 2.236 m, which is much smaller than the first two cases, and has a significant inhibitory effect on the maximum deviation.Table 4
Comparison of three positioning methods in another experiment.
MethodStatistics of Euclidean distance errors (unit: m)MaxRMSKF17.9452.932KF27.5452.696MWMM5.6592.236Figure 17
Comparison of the three methods in another experiment.Many groups of experiments have been carried out, and the results show that the localization result of KF2 is better than that of KF1. The final positioning result obtained by MWMM is greatly improved compared with both KF1 and KF2, especially for the points with large error.


## 6. Conclusion

To achieve accurate and reliable positioning in urban canyons, an enhanced map-aided GPS/3D RISS tightly combined positioning strategy is proposed. In the proposed strategy, the 3D RISS only based on the built-in CAN bus sensor without additional sensors is first constructed. Then, a simple but effective enhanced map is created. Based on the map, a KF method with the adaptive adjustment of observation noise variance is proposed to fuse the 3D RISS with GPS information and to achieve the preliminary positioning. In KF-based preliminary positioning method, the value of the observation noise variance matrix is determined only according to the building plane information which is contained in the enhanced map. Finally, a multiweight map matching method is proposed to correct the initial positioning results and obtain the final positioning results. The experimental results show that there are more obvious advantages of the proposed strategy than other methods in urban canyon, that is, the RMS achieve 32.8% accuracy improvement compared to conventional KF.Further research includes using sensors such as fisheye camera and machine learning methods to further effectively identify failed satellites, replacing KF algorithm with more effective algorithm, and carrying out more comprehensive experiments, especially in deep urban canyon environment.



---

*Source: 7650435-2022-04-14.xml*