GIL: a tightly coupled GNSS PPP/INS/LiDAR method for precise vehicle navigation

Accurate positioning and navigation play a vital role in vehicle-related applications, such as autonomous driving and precision agriculture. With the rapid development of Global Navigation Satellite Systems (GNSS), Precise Point Positioning (PPP) technique, as a global positioning solution, has been widely applied due to its convenient operation. Nevertheless, the performance of PPP is severely affected by signal interference, especially in GNSS-challenged environments. Inertial Navigation System (INS) aided GNSS can significantly improve the continuity and accuracy of navigation in harsh environments, but suffers from degradation during GNSS outages. LiDAR (Laser Imaging, Detection, and Ranging)-Inertial Odometry (LIO), which has performed well in local navigation, can restrain the divergence of Inertial Measurement Units (IMU). However, in long-range navigation, error accumulation is inevitable if no external aids are applied. To improve vehicle navigation performance, we proposed a tightly coupled GNSS PPP/INS/LiDAR (GIL) integration method, which tightly integrates the raw measurements from multi-GNSS PPP, Micro-Electro-Mechanical System (MEMS)-IMU, and LiDAR to achieve high-accuracy and reliable navigation in urban environments. Several experiments were conducted to evaluate this method. The results indicate that in comparison with the multi-GNSS PPP/INS tightly coupled solution the positioning Root-Mean-Square Errors (RMSEs) of the proposed GIL method have the improvements of 63.0%, 51.3%, and 62.2% in east, north, and vertical components, respectively. The GIL method can achieve decimeter-level positioning accuracy in GNSS partly-blocked environment (i.e., the environment with GNSS signals partly-blocked) and meter-level positioning accuracy in GNSS difficult environment (i.e., the environment with GNSS hardly used). Besides, the accuracy of velocity and attitude estimation can also be enhanced with the GIL method.


Introduction
Accurate and continuous navigation is one of the fundamental prerequisites for a reliable and intelligent driving system. Nevertheless, it is often difficult for a stand-alone sensor to meet the needs of robust navigation in complex scenarios (Groves et al., 2014;Li et al., 2021). Thus, multisensor data fusion, which takes full advantage of different sensors (e.g., Global Navigation Satellite Systems (GNSS), Inertial Navigation System (INS), Laser Imaging, Detection, and Ranging (LiDAR), camera), has become a hotspot in both academic and industrial sectors. 2015). With more GNSS satellites available, the GPS (Global Positioning System) + BDS + Galileo PPP has been confirmed to have faster convergence and higher accuracy than the GPS-only PPP solution (Guo et al., 2018;Li et al., 2018).
Despite its advantages, the multi-GNSS PPP cannot be immune to degradation in poor satellite visibility or weak constellation geometry (Zhang & Li, 2012), thereby making continuous and precise navigation in urban areas an intractable problem. To address this problem, a lot of work has been done to aid GNSS with INS (Cui et al., 2019;Gao et al., 2017). As for the PPP, Roesler (2009) and Shin (2009) proved that the PPP/INS integration could achieve superior positioning accuracy and continuity in both open-sky and GNSS difficult environments (i.e., the environment with GNSS hardly used). Compared to loosely coupled approaches, tightly coupled integration has been demonstrated to be more effective and robust when the satellite availability is limited (Abd Rabbou & El-Rabbany, 2015a, 2015b. However, when GNSS signals are blocked, the navigation result degrades rapidly owing to the error accumulation, especially for Micro-Electro-Mechanical System (MEMS) Inertial Measurement Units (IMUs) (Abd Rabbou & El-Rabbany, 2015a).
Fortunately, the integration of LiDAR and IMU, which shows excellent precision and reliability in local navigation, provides another solution for vehicle navigation. Typical methods, such as Laser Odometry and Mapping (LOAM) (Zhang & Singh, 2014), Catergapher (Hess et al., 2016), and Lidar-IMU Odometry (LIO)-mapping (Ye et al., 2019) can achieve decimeter-level positioning accuracy for indoor navigation tasks. These methods commonly formulate an optimization problem to determine the best estimation of LiDAR poses. Besides, there are also substantial researches that tried to realize the fusion through filter-based practices. For instance, Zhen et al. (2017) proposed a tightly coupled INS/LiDAR method based on a Gaussian particle filter and verified its effectiveness in various challenging conditions, such as a kidnapped robot. Qin et al. (2020) utilized an iterated error-state Kalman filter to achieve real-time ego-motion estimation. Generally, the optimization-based methods require significant computational resources to obtain high accuracy through an iterative process, while the filter-based algorithms have advantages in operational efficiency with comparable accuracy (Qin et al., 2020). However, the error accumulation is unavoidable for the INS/LiDAR odometry in long-duration navigation without available global corrections.
To take advantage of both local and global navigation and achieve global drift-free localization, many scholars have conducted in-depth research on the integration of GNSS, INS, and LiDAR. One of the fusion schemes utilizes GNSS/INS results to estimate gravity's orientation and predict LiDAR poses (Hess et al., 2016). However, the scheme does not consider GPS and IMU data processing and heavily depends on LiDAR registration. Alternatively, by extending the conventional GNSS/INS scheme with LiDAR Simultaneous Location and Mapping (SLAM), the authors (Chiang et al., 2020) presented a lane-level navigation approach for moving vehicles. Other integration strategies aided the tightly coupled LiDAR/IMU odometry with optional GPS position constraints (Shan et al., 2020;Koide, 2019). These strategies take GNSS as a "black box", describe the poses in a self-defined local framework and loosely merge local odometer results with GPS positioning results. Soloviev (2008) introduced a tightly coupled solution of GPS, two-Dimensional (2D) laser, and INS to improve the performance of 2D plane positioning in urban areas. Nevertheless, this study only utilized GPS-derived horizontal positions.
In this contribution, we propose a tightly coupled multi-GNSS PPP/INS/LiDAR (GIL) algorithm to perform three-Dimensional (3D) large-scale vehicle navigation in urban environments. Raw observations from multi-GNSS PPP, MEME-IMU, and 16-line LiDAR are integrated through an Extended Kalman Filter (EKF) to enhance the navigation performance in terms of position, velocity and attitude. The experiments in the GNSS challenging environments around Wuhan University are designed to evaluate this method. The remaining paper is organized as follows. Firstly, we introduce the system state description. Then the process model is investigated, followed by the multi-GNSS PPP and LiDAR measurement models. Thirdly, the overall flow of the GIL is organized. Subsequently, the experimental platform is briefly introduced, and the experiments and corresponding analyses are presented. Finally, we summarize this study and the prospect for the follow-up research.

State description
The state vector δx consists of the error states from INS, PPP, and LiDAR, which are represented by δx ins , δx ppp and δx lidar , respectively. The specific forms can be defined as: In δx ins , δθ e b denotes the three-dimensional attitude errors. The symbols b and e indicate the body coordinate (1) frame ( b-frame) and earth-centered earth-fixed frame ( e -frame), respectively. Note that the IMU frame coincides with the b-frame in this study. The vectors δv e b and δp e b are corrections of velocity and position in e-frame. Besides, the biases of the accelerometer and gyroscope b a , b g are modeled as random walks. The corresponding bias error variations δb a and δb g are included in the state vector and estimated online.
As for δx ppp , the position of the antenna center can be related to the INS position through the lever-arm. Therefore, only the receiver clock offsets of different constellation systems δt r , zenith tropospheric wet delay δZ wet , and Ionospheric-Free (IF) ambiguity δN IF are stored in the PPP state vector (Li et al., 2015).
Furthermore, the LiDAR state can be represented by the attitude error δθ e l and position error δp e l . Several LiDAR poses maintained by sliding window mechanism are stored in the state vector for the multistate constraint Kalman Filter model (Mourikis & Roumeliotis, 2007). Thus, the LiDAR-correlated state vector at epoch k can be expressed as:

Process model
The linearized continuous dynamic INS model (Shin, 2005) can be written as: where δṗ e b , δv e b and δθ e b refer to the derivative of position, velocity, and attitude errors, respectively. ω e ie is the angular rate of the e-frame with respect to the inertial frame ( i -frame), and f b and ω b ib are the measurements of accelerometer and gyroscope, respectively. R e b denotes the rotation matrix from b-frame to e-frame.
Specifically, the system model of the GIL method used for updating EKF is expressed as:  where I is the identity matrix, n t r , n zwd , n a , n g represent the noise of t r , d zwd , b a , b g measurements.
When a new LiDAR frame is available at epoch k+1 , the initial LiDAR pose R e l k+1 ,p e l k+1 can be derived from the simultaneous IMU pose R e b k+! ,p e b k+1 through the pre-calibrated extrinsic parameters between the IMU and LiDAR R b l , p b l , which can be written as: Then, the new LiDAR state δx k+1 lidar = δθ e l k+1 , δp e l k+1 T at epoch k + 1 will be added into the system state vector, and the corresponding system covariance is also augmented (Mourikis & Roumeliotis, 2007): where δx k+1,k is the system state vector after δx k+1 lidar is added, P k,k and P k+1.k are the covariance before and after the augmentation procedure, respectively. The non-zero elements in the jacobian matrix J can be expressed as:

Multi-GNSS PPP measurement model
In PPP processing, the ionospheric-free (IF) model is used to eliminate the first-order ionospheric delay, which can be written as: where s and r refer to the satellite and the receiver. P s r,IF and L s r,IF are the ionosphere-free psedorange and carrier phase observations, respectively, ρ=|p r − p s | stands for the distance between the receiver and the satellite, c is the speed of light, t s and t r represent the offsets of satellite clock and receiver clock, T s r denotes the tropospheric delay. The wavelength and ambiguity of the IF model are represented by IF and N s r,IF . e s r,IF and ε s r,IF denote the measurement noise and multipath error of the IF pseudorange and carrier phase, respectively.
According to (Böhm et al., 2006), the tropospheric delay is composed of dry and wet components, which are computed with their zenith components Z dry , Z wet and corresponding mapping functions m dry , m wet . In this article, the dry component is corrected by a prior model, and the wet component is modeled as random walking. In the multi-constellation case of GPS (G), BDS (C), and Galileo (E), the linearized PPP observation equations can be expressed as: where p r and p s are the position vector of receiver and satellite, respectively. µ denotes the unit vector from the satellite to the receiver.
is the ambiguity correction vector, and δt r = δ t r,G δt r,C δt r,E T is the receiver clock correction vector. In the case of multi-GNSS observations, the different hardware time delays called Inter-System Biases (ISBs) should be taken into account (Li et al., 2015), where ISB G,C and ISB G,E are the ISBs of BDS and Galileo with respect to GPS, which are modeled as random walks. Thus, the receiver clock offsets δt r = δ t r,G δt r,C δt r,E et al. Satellite Navigation (2021) 2:26 can be re-parameterized as δt r = δ t r,G ISB G,C ISB G,E T .
In the integrated system, the simultaneous IMU state is used to predict the GNSS-derived position, but the antenna center has a different reference point with respect to the INS predicted position, which leads to lever-arm offsets. The lever-arm is precisely measured in advance and its correction can be expressed as: where l b denotes the lever-arm. The relationship between δp r and δp e b can be expressed from the differential of (11): The linearized GNSS observation equations can be written as: where k refers to the time epoch, δx denotes the system state vector, r P k and r L k are the residuals of psedorange and carrier measurements, n P k and n L k represent the corresponding noises. P k ins and L k ins are the INS-predicted GNSS measurements, and the designed matrix H P k , H L k can be determined directly from (9) and (12).

LiDAR measurement model
Before introducing the LiDAR measurements, the sliding window mechanism is described in advance. As shown in Fig. 1, the threshold N of LiDAR frames stored in the state vector is pre-defined, and the oldest frame will be discarded when the window is full. The newly obtained LiDAR frame will be matched with other LiDAR scans in the sliding window to establish the geometry relationships between frames. We note that the N value should not be set too large. A small window benefits computational efficiency, but the scan matching tends to be inaccurate with a large distance between frames, particularly when the vehicle moves fast on the road. Considering some unreliable observations exist in the registration, a Chi-square test (Sun et al., 2018;Zuo et al., 2019) is employed to remove outliers from the scan matching, and only measurements that pass the test will be used for EKF measurement update. When a new LiDAR frame arrives, feature extraction is first performed to select high-quality points that are approximately on edges or planes in the scan (Zhang & Singh, 2014). Then the edge or plane features are de-skew corrected (Ye et al., 2019) and projected to the old LiDAR frames to find the nearest line or plane correspondences. The spatial octree (De Berg et al., 2008) is constructed for fast indexing in this process. For more details on matching the point-line pairs and point-plane pairs refer to .
In the case of edge matching, assuming that p edge,l k+1 i is an edge feature in the newest LiDAR frame k + 1 , its projection in the frame k is p  The LiDAR observation equations can be represented by the distances in formulas (14) and (15): where r edge i , r plane i are observation residuals, and n edge and n plane represent noises of edge and plane observations, respectively. H k,k+1 edge,i and H k,k+1 plane,i are the design matrixes, which can be calculated following the chain rule: The subterms in the formula,  (14) and (15). The relationship between p edge,l k i , p plane,l k i and LiDAR states x k lidar , x k+1 lidar can be derived from the projection formulas: ,p e l k and R e l k+1 ,p e l k+1 are the poses of LiDAR frame k and k+1 , which can be obtained by referring to (5). Thus, the non-zero elements of edge or plane measurement in ∂δp l k i ∂x lidar can be expressed as:

Measurement update
When GNSS or LiDAR pre-processed observations are available at epoch k+1 , the corresponding measurement update can be performed according to (13) and (16), which follows the standard EKF routines (Qin et al., 2020) as: (20) where n k+1 is the observation noise that meets the zeromean Gaussian noise assumption, whose covariance is R n . K denotes the near-optimal Kalman gain, and the symbol ⊗ represents the compensation of the states. Meanwhile, the updated covariance matrix P k+1,k+1 can be calculated by:

GIL algorithm implementation
Based on the mathematical model introduced above, the main process of the GIL method can be organized as Fig. 2. After the system initialization, the INS dynamic alignment (Han & Wang, 2010) is performed to obtain the initial pose of the vehicle, then the bias-compensated IMU data will be adopted for INS mechanization.
The EKF update will be performed to predict the states as well as the covariance. Before constructing the IF (21)    (Fang et al., 2009). Similarly, LiDAR feature detection (Zhang et al., 2014) and de-skew (Ye et al., 2019) are carried out in LiDAR data preparation, and the initial LiDAR pose is added to the system state through the augmentation procedure. When there are pre-processed GNSS and LiDAR measurements according to the mathematical models, the corresponding measurement updates will be performed sequentially to update the filter states. Furthermore, the corrected gyroscope and accelerometer biases will be fed back to the next IMU data for restraining the INS divergence.

Experimental platform and strategies
To evaluate the performance of the proposed GIL method, experiments with different GNSS conditions were conducted in Wuhan, Hubei, on 27th October 2020. A road with tall trees and buildings on both sides was selected for the tests, which reflects the real situation of frequent signal interruptions in urbanized areas. Some typical scenarios in the experiments are shown in Fig. 3. The detailed strategies used in the GNSS PPP and LiDAR processing are listed in Tables 1 and 2, respectively. For the GNSS PPP, multiple GNSS systems, including GPS (G), BDS (C), and Galileo (E) were considered in the GNSS processing. The precise orbit and clock products used in the processing were provided by the Center for Orbit Determination in Europe (CODE). As for LiDAR, the size of the sliding window was set to 4, and the observation noise was predefined separately for the two kinds of observations based on the LiDAR resolution.
The experiment data were collected with a mobile platform shown in Fig. 4. A LiDAR (Velodyne Corporation, 2021) with 16 lines was mounted above a MEMS-IMU (Analog Devices Corporation, 2019) to ensure unimpeded scans during the experiments. A tactical IMU (StarNeto Corporation, 2014) was installed in the middle of the platform. Additionally, a GNSS antenna (NovAtel Corporation, 2015), which was linked to a GNSS receiver (Septentrio Corporation, 2019), is equipped. Time synchronization on the hardware level unifies timestamps of all collected data to GPS time. For spatial synchronization, pre-calibration was performed to identify the extrinsic parameters between LiDAR and IMU (Lv et al., 2020), and the lever arms between IMUs and the GNSS antenna were measured accurately.
Performance specifications of the two different grade IMUs are shown in Table 3. During the experiments, the sampling rate of the GNSS receiver was set to 1 Hz, and that of LiDAR was configured to 10 Hz. Both the

MEMS-IMU and
Tactical IMU worked at 100 Hz. In addition, another GNSS receiver (Septentrio PolaRx5) was fixed on the rooftop of a tall building as a base station. The reference solution calculated by commercial Inertial Explorer 8.9 software (NovAtel Corporation, 2018) was the smoothed results of multi-GNSS Post-Processing Kinematic (PPK) and INS tight integration. The overall positioning accuracy of the solution can reach the centimeter-level with a fixed rate of more than 90% and high-precision extrapolated results of tactical IMU are provided in non-fixed periods (NovAtel Corporation, 2018). On this basis, the accuracy of the reference solution is higher than that of the solutions to be analyzed by an order of magnitude, which can facilitate a fair and reliable performance comparison between different algorithms in the following evaluation.

Result analysis
In this section, we separately analyzed the vehicle navigation results in GNSS partly-blocked environment and GNSS difficult environment. In each experiment, the multi-GNSS PPP solution and satellite availability are first introduced. Then, the Position, Velocity, and Attitude (PVA) estimation of the GIL method is presented. Meanwhile, the multiple GNSS PPP/INS tightly coupled solution is used for comparison to assess the GIL performance in urban vehicle navigation.

PPP solution in GNSS partly-blocked experiment
The first experiment focuses on GNSS partly-blocked environment. Figure 5 shows the top view of the route and the multi-GNSS PPP (G + C + E) positioning result, where colors represent the corresponding error distribution intervals. It is noticed that in some road sections, the harsh environment caused frequent and slow re-convergence of PPP, resulting in positioning discontinuity and large abnormal offsets. Clearly, the multi-GNSS PPP result in this environment can hardly meet the requirement of high accuracy and continuous navigation.
In this study, the number of visible satellites, satellite elevation, continuity of satellite signal tracking, and the corresponding Position Dilution of Precision (PDOP) are utilized to verify the GNSS availability. The number of available satellites is shown in Fig. 6a. The enhancement of GNSS observability is considerable when using  (Li et al., 2015). The average number of satellites for the G + C + E combination is 13.46, which is 3.5, 2.5, and 3.3 times the number of the visible satellites for the single GPS, BDS, and Galileo systems, respectively. According to the continuity of available satellites depicted in Fig. 6b, interruptions occur frequently in some segments, particularly for the satellites with low elevation. This is mainly due to the fact that roadside houses and trees cause severe GNSS signal obstructions. As one can see from the time series of the PDOP values shown in Fig. 6c, the dramatic jumps in PDOP always accompany severe tracking loss, and the maximum value of PDOP reaches 49.25. These facts reflect that the discontinuous tracking of satellites is a challenge for GNSS positioning.

GNSS partly-blocked experiment analysis
Different from the absolute PPP, the PPP/INS integration system can produce high-frequency navigation output and attitude information due to the characteristics of INS. To further overcome the drawback of INS divergency in GNSS-constrained areas, we extend the fusion with LiDAR measurements. The solution of tightly coupled multi-GNSS PPP/INS is employed as a contrast to the GIL solution.
The position errors in three directions over time are displayed in Fig. 7, and the corresponding RMSEs, and the maximum errors during the experiment are summarized in Table 4. It can be found that the position differences in three directions have been improved noticeably with the aid of LiDAR data. The RMSEs of position are 0.81, 0.76, and 0.82 m for the PPP/INS method in east, north, and up components, while these are 0.30, 0.37, and 0.31 m for the GIL method. Additionally, the maximum positioning errors are reduced from 7.01, 3.58, and 7.55 m in east, north, and up directions to 1.39, 0.79, and 1.00 m, with improvements of 80.2%, 77.9%, and 86.8%, respectively. Besides, the region with GNSS signal blockage is usually rich in scene features, which is conducive for LiDAR registration. Thus, the LiDAR measurements can effectively assist the system to infer the vehicle movement in the environments with no GNSS signals.
To further analyze the absolute positioning performance with our method, we calculated the distribution frequency of position errors in different intervals, as plotted in Fig. 8. According to the statistical diagram, the percentage of positioning difference within 0.5 m is 19.5% for the PPP/INS solution, which increases to 41.9% with LiDAR aided. Besides, the positioning accuracy for the GIL is within 1.0 m over 98% of the time in the field test.
The results indicate that the addition of LiDAR can constrain the divergence of IMU error and contribute to the accurate positioning in urban areas. These improvements may benefit applications that require continuous and precise dynamic positioning.
For vehicle navigation, accurate velocity is another crucial factor. Figure 9 shows the velocity difference between the two methods. For the PPP/INS solutions, velocity offsets above 0.5 m/s mainly occur in low GNSS availability periods. Comparatively, after the inclusion of LiDAR data, the velocity offsets are within 0.1 m/s during most of the time. Besides, the statistics in Table 5 illustrate that the velocity RMSEs are 0.06, 0.03, and 0.04 m/s for the PPP/INS/LiDAR method in east, north, and up directions, with improvements of 60.0%, 50.0%, and 60.0% in comparison with the INS-aided PPP. Similarly, the maximum error drops from 0.66, 1.08, and 1.10 m/s in   the east, north, and up components to 0.23, 0.28, and 0.22 m/s, respectively. Note that the velocity is propagated by INS mechanization, whose performance mainly depends on the IMU. In this situation, LiDAR observations can effectively mitigate the drift problem of INS, especially for MEMS-IMU, which makes a great contribution to the velocity estimation. Similar to position and velocity, attitude is another vital information for vehicle applications. However, when the GNSS performance is degraded, the accuracy of attitude will also be affected especially for the yaw component. In this experiment, the two solutions' attitude differences over time are depicted in Fig. 10, and corresponding RMSEs of attitude are shown in Table 6. In terms of pitch and roll, the RMSEs of the two methods show comparable performance. Comparatively, the achievable accuracy of the yaw angle in both methods is much worse. This is due to the case that the yaw angle is unobservable in INS or LiDAR/INS systems (Hesch et al., 2013). Specifically, the RMSE of yaw angle for the PPP/INS solution is 3.97 degrees, and the corresponding figure for the GIL method is 1.45 degrees, which has an improvement of 63.5%. Generally, the LiDAR aiding can help the attitudecorrelated parameters to be estimated more accurately.

PPP solution in GNSS difficult experiment
To verify the universality of the GIL method, we conducted another onboard vehicle experiment. As illustrated in Fig. 11, the condition was more challenging than the previous one with dense, tall houses and trees beside the road. The track length is approximately 1700 m and the vehicle traveled about 7 min. Severe signal blockages caused frequent re-convergence and the insufficient number of the observed satellites, resulting in poor positioning and frequently discontinuous PPP solution in this experiment. According to Fig. 11, though multi-GNSS observations are used, the PPP solution can hardly maintain meter level positioning accuracy. Figure 12 shows the satellite availability over time. Compared with the previous experiment, severe interruptions occur even for the satellites with high elevation. The number of available satellites is 50% less than the previous test, and more fluctuates. The average number of GPS satellites is 2.23, while the corresponding value for the G + C combination is 4.49, which further increases to 7.74 for the G + C + E combination. We also notice that, compared to the previous test, the average PDOP variation degrades from 3.54 to 6.99, and the corresponding Standard Deviation (SD) increases from 3.40 to 5.95. Obviously, the GNSS observation condition in this experiment is much worse than the previous test, which leads to poor PPP result.

GNSS difficult experiment analysis
The position differences of the two methods with respect to the reference solution are shown in Fig. 13, and the corresponding position offsets distribution is shown in Fig. 14. The absolute RMSE of the PPP/INS positioning in this experiment reaches 4.56 m, which is about three times more than the value in the previous field test. Besides, position errors of more than 10 m can be found in both vertical and horizontal directions.
As shown in Fig. 13, the results with the inclusion of LiDAR measurements are significantly improved compared to the PPP/INS solution. Table 7    aided. Generally, the positioning performance is significantly enhanced when LiDAR data are augmented, which means the proposed fusion method can ameliorate the positioning even in bad environments. This progress is mainly due to the constraints between adjacent epochs provided by LiDAR measurements. The improved accuracy of position can benefit those position-critical tasks. For velocity estimation, outperformance is also shown in the results of the GIL method, and the statistics are summarized in Table 8. According to the velocity offsets reflected in Fig. 15, the velocity error is within 0.3 m/s during most of the periods when using the GIL method. Specifically, the velocity RMSEs are 0.85, 1.09, and 0.93 m/s for the PPP/INS method in east, north, and up directions, which are improved by 42.4%, 56.9%, and 53.8% after integrating LiDAR together. Those improvements are mainly due to the geometry restraint established by LiDAR measurements. As presented in Fig. 15, some obvious offsets occur in the PPP/INS velocity estimation. Referring to Fig. 12, it is noticed that the corresponding procedures are usually with insufficient satellite observations, and the deterioration over time occurs when GNSS cannot provide effective constraints for the PPP/INS integration. As for the GIL method, additional observations from LiDAR can deliver accurate velocity output, which helps maintain a high-accuracy velocity solution even when there are few available satellites.
The time series of attitude error is presented in Fig. 16, and the related RMSEs and maximum errors are listed in Table 9. Similarly, the results show higher accuracy in the estimation of roll and pitch angle, which may be related to the better observability of INS in the two directions. For the heading component, the RMSE of the GIL method has a significant improvement of 68.9% compared to the PPP/INS method. Slight benefits are also visible in roll and pitch angles with the improvement of 45.5% and 41.2%.

Conclusion
In recent years, the demand for stable, continuous, and accurate navigation has gradually become prominent in autonomous driving applications. In this contribution, we proposed a tightly coupled multi-GNSS PPP/INS/ LiDAR method to improve the navigation performance in heavily urbanized areas. By formulating an EKF system driven by INS mechanization, the observations from LiDAR and PPP constructed the corresponding measurement updates to restrict the rapid error growth of INS. An efficient sliding window strategy was designed to deal with LiDAR data processing. Several field experiments in typical urban areas were conducted to validate the model, and the performance of the GIL method was assessed in terms of position, velocity, and attitude. The statistics data demonstrate that compared to the GPS-only constellation, there are about three times more available satellites when three constellations of GPS, BDS, and Galileo are considered. However, the observation quality is still limited due to GNSS signal interruptions. In addition, the experimental results illustrate that our GIL method can maintain more accurate and stable navigation than the traditional tightly coupled PPP/ INS solution. In GNSS partly-blocked environments, the position RMSEs in east-north-vertical directions are dramatically reduced from 0.81, 0.76, and 0.82 m to 0.30, 0.37, and 0.31 m compared with the PPP/INS solution. When the conditions are more difficult for GNSS measurements, the GIL method also shows an excellent performance in positioning with the improvements of 50.0%, 36.5%, and 43.1% in east, north, and vertical components, respectively. It is worth mentioning that the maximum positioning error has also been significantly reduced when using the GIL method. According to the experiment in GNSS partly-blocked areas, the velocity RMSEs of the PPP/INS solution are 0.15, 0.06, and 0.10 m/s in east, north, and up directions, respectively. The addition of LiDAR to the above solution can reduce the RMSEs to 0.06, 0.03, and 0.04 m/s. The attitude-related parameters can also be accurately estimated with integration of multi-GNSS and LiDAR data. In both GNSS partlyblocked and difficult environments, the GIL method has superior performances in the determination of roll and pitch angles, with average improvements of 37.5% and 25.9%. Meanwhile, significant improvements by 63.5% and 68.9% for both experiments are found in the yaw direction compared to the multi-GNSS PPP/INS results.
In conclusion, this paper presented a novel tight fusion method, which shows great potential in position, velocity, and attitude determination. The results from the field tests demonstrate that the GIL fusion outperforms the multi-GNSS PPP/INS integration and achieves better accuracy in dynamic vehicle scenarios. The further study is to take full advantage of LiDAR measurements. A point cloud map will be established for registration and closedloop detection. Meanwhile, more efforts should also be devoted to eliminating the adverse effects of the outliers in each sensor measurements.

Availability of data and materials
The datasets used and analyzed in this study are available from the corresponding author on reasonable request.