On the role of lidar measurements in speeding up precise point positioning convergence

Global navigation satellite system (GNSS) and light detection and ranging (lidar) are well known to be complementary for vehicle positioning in urban canyons, where GNSS observations are prone to signal blockage and multi-path. As one of the most common carrier-phase-based precise positioning techniques, precise point positioning (PPP) enables single-receiver positioning as it utilizes state-space representation corrections for satellite orbits and clocks and does not require a nearby reference station. Yet PPP suffers from a long positioning convergence time. In this contribution, we propose to reduce the PPP convergence using an observation-level integration of GNSS and lidar. Lidar measurements, in the form of 3D keypoints, are generated by registering online scans to a pre-built high-definition map through deep learning and are then combined with dual-frequency PPP (DF-PPP) observations in an extended Kalman filter implementing the constant-velocity model that captures the vehicle dynamics. We realize real-time PPP (RT-PPP) in this integration using the IGS real-time service products for vehicle positioning. Comprehensive analyses are provided to evaluate different combinations of measurements and PPP corrections in both static and simulated kinematic modes using data captured by multiple receivers. Experimental results show that the integration achieves cm-level accuracy and instantaneous convergence by using redundant measurements. Accordingly, for classical PPP accuracy of 10 cm and convergence within minutes, respectively, lidar input is only required once every 10 s.


Introduction
Accurate and continuous positioning is the cornerstone of modern driving systems. Reid et al. (2019) proposed that the highly anticipated autonomous vehicles should possess the positioning accuracy of 10, 10 and 48 cm in the across-theroad, along-the-road and vertical directions, respectively, on city roads. Among numerous sensing devices used for positioning purposes, global navigation satellite system (GNSS) receivers are considered the fundamental component as they provide globally referenced solutions (Joubert et al. 2020). Unfortunately, GNSS signals suffer from interference and blockage in dense urban environments, where accurate positioning is the most demanded (Groves et al. 2013). Therefore, the integration of GNSS and other sensors such as inertial measurement unit (IMU), cameras and light detection and ranging (lidar) scanners is widely used to improve the accuracy and continuity of positioning (Wen and Hsu 2022), wherein continuity refers to the probability of maintaining a specified level of positioning accuracy throughout the operation (Teunissen and Montenbruck 2017).
Among the two most commonly used precise GNSS positioning techniques, namely real-time kinematic (RTK) and precise point positioning (PPP), the latter is gaining increasing attention in multi-sensor positioning and navigation. Compared with RTK, PPP utilizes precise state-space representation (SSR) corrections for satellite orbits and clocks and hence does not require reference GNSS receivers in a close proximity, leading to higher service availability (Zumberge et al. 1997). Despite the potential of achieving cmlevel accuracy using single-receivers, the major drawback of PPP is the long convergence time that is due to the presence of real-valued carrier phase ambiguities and unmodeled ionospheric delays prohibiting integer ambiguity resolution (Teunissen 1997;Guo et al. 2017). The convergence is even longer for kinematic receivers as compared with static ones and can take up to one hour (Wang et al. 2019;Erol et al. 2021). Moreover, high-quality PPP products, such as the final products provided by International GNSS Service (IGS), can require weeks to prepare, hindering real-time applications (Kouba and Héroux 2001).
In contrast, various real-time PPP (RT-PPP) services with lower nominal accuracy are now offered by agencies including commercial companies, Analysis Centers (AC) and IGS. For instance, Trimble CenterPoint RTX and Novatel TerraStar-C Pro have been extensively studied under different signal-reception conditions to show cm-to dm-level positioning accuracy after a convergence period of tens of minutes (Alkan et al. 2020b(Alkan et al. , 2020aİlçi and Peker 2022). However, such commercial services are not openly accessible. Meanwhile, with the contributions of ACs, IGS launched the IGS real-time service (RTS) in 2013, which provides precise satellite orbit and clock corrections continuously and is freely available over the Internet (Hadas and Bosy 2015). The stated orbital accuracy, clock accuracy and latency of IGS RTS are approximately 5 cm, 0.3 ns and 25 s, respectively (IGS 2019), which have been confirmed by numerous studies, e.g., Elsobeiey and Al-Harbi (2016) and Li et al. (2022a, b). Furthermore, it has been found that RTS, similar to the majority of other PPP services, offers most accurate solutions for the global positioning system (GPS) constellation with the product availability exceeding 95%, while delivering the positioning performance comparable with post-processed PPP (Wang et al. 2018b;Kazmierski et al. 2020;Alkan et al. 2022).
In spite of the impressive accuracy of IGS RTS corrections, RT-PPP is still restricted by the slow convergence, which prevents it from being implemented for real-time vehicle positioning. This obstacle can be overcome by the provision of complimentary PPP-RTK corrections (Teunissen and Khodabandeh 2015;Zhang et al. 2022a) and/or by the realization of multi-sensor integration (İlçi and Toth 2020;Du et al. 2021). Among the popular technologies for aiding GNSS positioning, lidar is known for its high accuracy and insensitivity to illumination (Li et al. 2023). Since the raw point clouds collected by lidar are only locally referenced for ego-motion estimations (Zhang et al. 2021b;Li et al. 2022b), high-definition (HD) maps, which supply geospatial information of the road environments, are often jointly used to transform lidar measurements to align with their GNSS counterparts for fusion . In addition to the traditional integration of GNSS and lidar aided by HD maps at the solution level (i.e., loose coupling), which mainly aims for multi-path and outlier mitigation such as Wen et al. (2019Wen et al. ( , 2020 and Zhong and Groves (2022), lidar has also been exploited for integer ambiguity resolution in RTK by improving the precision of the float solutions (Qian et al. 2020;Li et al. 2021b;Zhang et al. 2022b). Most recently, the fusion of PPP and lidar have been explored for continuous vehicle positioning by decreasing the convergence time of PPP. Li et al. (2021a) proposed P 3 -LOAM, which combines singlefrequency PPP and lidar scan matching, yet the integration is at the solution-level and only achieves m-level positioning accuracy. Observation-level integrations (i.e., tight coupling) have also been studied, e.g., Li et al. (2021c) and Li et al. (2022bLi et al. ( , 2023 exhibited dm-level accuracy using PPP, IMU and lidar. However, they rely on geometrically extracting and matching edge and planar features from lidar point clouds to derive the measurements, which can be sparse or unavailable in complex and dynamic environments. On the other hand, recent developments of deep learning in the scope of point clouds have enabled datadriven approaches for describing and registering features, which are more robust against irregular geometries Si et al. 2022). Yet, lidar complemented by HD maps and deep learning has not been leveraged to improve the convergence speed and accuracy of RT-PPP.
We propose an observation-level integration of dualfrequency PPP (DF-PPP) using IGS RTS and learningbased, map-aided lidar for real-time positioning and conduct comprehensive analyses to study the role of lidar for speeding up the convergence. The main contributions include: (1) an integrated positioning framework wherein georeferenced lidar measurements are produced through deep learning using an HD map comprising downsampled point clouds of the roads, which are then weighted by their intensity values and directly fused with their DF-PPP counterparts in an extended Kalman filter (EKF) implementing the constant-velocity model; (2) a novel weighting strategy of lidar measurements based on intensity, which is further adjusted by point distributions; (3) detailed evaluations which show that DF-PPP using highrate IGS RTS precise products can be severely affected by the temporal correlation of measurements, but is nevertheless capable of real-time vehicle positioning given the provision of lidar for accelerating the PPP convergence.
An overview of the structure of the present contribution, together with the main findings, is presented in Fig. 1. Accordingly, we commence with a discussion on the way lidar measurements are extracted through deep learning and an HD map. We then present our proposed observation-level PPP-Lidar integrated model and the measurement weighting scheme. After that, a simulation strategy and experimental details are discussed, followed by the static positioning performance comparison between tested IGS precise products. Afterward, the results of the simulated kinematic experiment are provided and discussed. Finally, we provide a summary with concluding remarks.

Lidar measurements extracted from point cloud registration
The proposed integrated positioning method combines DF-PPP observations and lidar measurements directly at the observation level in an EKF adopting the constant-velocity model. The lidar measurements are generated by registering a rover scan collected by the on-board lidar sensor with a georeferenced reference scan from the pre-built HD map, from which the 3D coordinates of the matched keypoints contribute to the estimation of the position of the vehicle. This section provides the procedure of extracting such lidar measurements using deep learning, which will then be integrated with DF-PPP.

Definitions
Two different coordinate systems are used throughout this paper. The geocentric ITRF2014 system, denoted as e-frame, is used for GNSS observations and precise PPP products, as well as the georeferenced lidar point coordinates. The l-frame is defined as the local Cartesian system centered around the lidar sensor, in which the point clouds are originally collected. It is assumed that the lidar sensor and the GNSS receiver have been aligned such that their measurements are referenced to the same point, which represents the vehicle position. In addition, we make the following definitions: 1. GNSS measurement interval ( o ): the frequency of GNSS observations in seconds. 2. Precise product sampling rate ( p ): the frequency of receiving precise products for GNSS orbital and clock corrections in seconds. In case the orbital and clock corrections have different sampling rates, p refers to the smaller value, which is often the update frequency of the clock products. 3. Lidar availability: the frequency of using lidar data in the integrated positioning method in seconds.

Learning-based keypoint extraction
In order to georeference the collected lidar point clouds so that they can be fused with their GNSS counterparts in e-frame, a globally referenced HD map of the road environment is required. In this work, we use 3D point clouds collected on the roads from a previous time in their original l-frame as the HD map, each with a georeferencing matrix to transform it to e-frame. The point clouds are also downsampled by the voxel size of 20 cm. This effectively reduces the storage size, while making the registration more challenging than using the original scans. For each rover scan collected by the lidar sensor, a nearby overlapped reference scan is selected from the HD map. Corresponding keypoints from both 3D point clouds are matched using MS-SVConv, which is a multi-scale deep neural network that shows state-of-the-art registration speed and accuracy (Horache et al. 2021). With the model pre-trained with extensive data, the network computes 32-dimensional feature vectors for randomly selected points from the two scans. The matching keypoints are found as the nearest neighbors based on the Euclidean norms of such features vectors, whose majority are unfortunately incorrect matches. Random sample consensus (RANSAC), which is an iterative estimation method for data that contains a large amount of outliers, is therefore used to identify and remove the outliers (Fischler and Bolles 1981). As output, the successfully registered keypoints have coordinates in l-frame referenced from the lidar sensor, as well as their coordinates in e-frame obtained from the HD map.

Lidar observation equation
The vehicle position can be estimated as the translational component of a transformation that registers the rover scan to the reference scan since it is equivalent to the origin of the rover scan. For keypoint i extracted at epoch t , the lidar observation equation to estimate such a transformation is given as follows: in which the vector y i (t) contains the 3D coordinates of the keypoint i in l-frame from the rover scan, with e i (t) being Overview of the paper structure (bold text) with the main findings and/or contributions (italic text) their measurement residuals. c i (t) is the vector of the coordinates of the same keypoint in e-frame obtained from the reference scan. Finally, R(t) and x r (t) denote the rotational and translational components of the estimated transformation, with the latter being also the unknown vehicle position. Since the unknown R(t) is intertwined with the measurements, (1) needs to be cast into a mixed measurement model in order to formulate the EKF (Teunissen 2003). The mathematical details of this model can be found in our previous work (Zhang et al. 2021a).

PPP-Lidar integrated model
In this section, we discuss the observation equations of the adopted DF-PPP model, together with their underlying EKF formulation, allowing one to integrate lidar and PPP at the observation level. In addition, a weighting strategy is presented to carefully consider the role of lidar keypoints in the integration.

DF-PPP observation equations
In this work, we adopt the ionosphere-unknown DF-PPP model. With precise products for satellite orbits and clocks, as well as a priori PPP corrections applied (Kouba 2009), at epoch t , let Δp r,j (t) and Δ r,j (t) , respectively, denote the observed-minus-computed code and phase observables for receiver r tracking m satellites on frequency j (j = 1, 2) , the DF-PPP observation equations read as follows (Teunissen and Montenbruck 2017) with Δx r (t) being the increment to the unknown position, whose satellite-to-receiver direction unit vectors form the m × 3 matrix G r (t) . 1 represents a matrix of ones. By defining j = 2 j ∕ 2 1 with j being the wavelength of the j th frequency, dt r (t) and r (t) denote the estimable receiver clock and ionospheric delays experienced on the first frequency containing the code biases, respectively. a r,j represents the real-valued phase ambiguities including the phase biases on frequency j . Lastly, e p,j (t) and e ,j (t) denote the code and phase measurement residuals. Apart from a r,j , which is expressed in cycles, all parameters are expressed in units of length.

EKF time update
Since (1) and (2) have the unknown position x r (t) in common, the lidar and GNSS measurements can be integrated in an EKF. The real-valued phase ambiguities a r,j in (2) are assumed constant in time, unless cycle slips occur, while the other quantities are treated as being unlinked in time. The EKF time update also applies a constant-velocity model, in which the acceleration of the vehicle motion is modeled by a zero-mean white-noise vector with three distinct spectral densities in the East-North-Up (ENU) directions. Therefore, the EKF state vector, including the vehicle position, velocity and carrier phase ambiguities at epoch t , respectively, is time-updated as follows: where the superscript (⋅) TU refers to the EKF time update, represent the aforementioned state parameters estimated at the previous epoch t − o and the transition matrix t|t− o reads as (Teunissen 2001) Accordingly, the error-variance matrix of these timeupdated parameters can be constructed as in which Q t− o is the error-variance matrix of the state vector from the previous epoch. In addition, with S = diag( 2 E , 2 N , 2 U ) formed by the ENU spectral densities of the vehicle acceleration and R e transforming S to e-frame, where diag(⋅) is the diagonal matrix operator, the system noise matrix S t is derived below where ⊗ is the Kronecker product (Henderson et al. 1983) and blkdiag(⋅) is the block diagonal matrix operator. The square matrix 0 with the dimension of 2(m − 1) indicates that a r,j are held constant over time.

EKF measurement update and measurement-weighting
With the time-updated parameters in place, they are appended to the lidar and GNSS measurements in (1) and (2), if they are available, to obtain the filtered solution in the weighted least squares (WLS) manner (Teunissen 2003). (3) Since lidar points with lower intensity are less precise, we weight the lidar measurements using the following intensity-based weight matrix by assuming n keypoints extracted at epoch t (Wujanz et al. 2017) in which i is the normalized intensity value for keypoint i , ranging between 0 and 1, with 1 representing the highest intensity. a and b are constants empirically determined depending on the lidar sensor specifications. W i denotes the weight matrix for the 3D coordinates of keypoint i , while W L is the weight matrix for all lidar measurements at epoch t used in the EKF.
In addition, due to RANSAC randomly selecting subsets of data for estimations, point cloud registrations performed by matching learning-based feature vectors, as described earlier, can occasionally be incorrect and unidentifiable. This often occurs when extracted keypoints are clustered (Zhang et al. 2021a). To address this issue, the position dilution of precision (PDOP) concept is applied on lidar data as follows (Srinara and Chiu 2022): where tr(⋅) represents the trace of a matrix and the n × 3 matrix G L (t) contains the keypoint-to-sensor direction unit vectors for the n lidar keypoints extracted at epoch t . (8) is used to scale W L to uniformly downweight lidar measurements when they have likely been derived from an incorrect registration due to clustered keypoints.
The GNSS observations are weighted in the EKF by the elevation angle of each satellite s (s = 1, … , m) . By defining p and as the zenith-referenced standard deviations of the code and phase observations, respectively, the weight matrix for all GNSS observations can be constructed as , respectively (Eueler and Goad 1991).
In order to mitigate the effects of measurement outliers and cycle clips, we apply detection, identification and adaptation (DIA) following the Chi-squared distributed test statistics and w-test statistics, whose procedure can be found in Teunissen and Montenbruck (2017). Finally, the weight matrix of the measurement update consists of the individual weight matrices of the measurements and the prediction which can be written as W = blkdiag W L , W G , Q −1 t .

Experimental setup
In order to comprehensively evaluate the proposed positioning method, GNSS observations from 20 Continuously Operating Reference Stations (CORS) distributed across Australia were collected and used as shown in Fig. 2. The GNSS receivers installed at these stations are either Septentrio PolaRx5 or Trimble NetR9, while they are equipped with a range of antennas including Javad RingAnt-DM, Leica AR25, Trimble 598,000.00 and ASH701945C_M.
Since it is not feasible to collect lidar data on such sites simultaneously or to move these fixed stations, point clouds from the KITTI dataset (Geiger et al. 2013) were simulated as if both types of data had been obtained on the same kinematic trajectory. This section presents the data collection details and the kinematic experiment simulation strategy.

GNSS data collection
The GNSS observations from the 20 selected stations were acquired from Geoscience Australia GNSS Data Archive, while two sources of IGS precise products were obtained including final products and RTS. The former was downloaded from Crustal Dynamics Data Information System (CDDIS), and the latter was recorded in real time from the single-epoch SSRA01 stream using BKG Ntrip Client (BNC). Meanwhile, 13 h of GPS observations from each station on 04 August 2022 were collected for the experiments, thereby generating 260 1 h samples in total, with 9

Simulation of kinematic experiment
For assessing the proposed PPP-Lidar integration in the simulated kinematic experiment, which will be shown later, a 1 h trajectory from the KITTI dataset is utilized with its associated lidar point clouds, which are split as rover and reference scans. The point clouds from KITTI dataset have been collected by a Velodyne HDL-64E with 64 channels with a maximum measurement range of 120 m, while the horizontal and vertical fields of view are 360° and 26.9°, respectively (Geiger et al. 2013). Figure 3 illustrates this trajectory and the time series of the ground truth acceleration of the vehicle, which can be classified as zero-mean white noise on average, agreeing with the constant-velocity assumption. The reference scans are given georeferencing matrices based on the ground truth poses of the vehicle provided in the dataset to construct the HD map. Each rover scan is kept approximately 10 m away from its corresponding reference scan to ensure a realistic amount of overlap to enable point cloud registration.
Although the GNSS observations were recorded on static points, they are transformed using the aforementioned trajectory to simulate the kinematic experiment. For each 1 h sample of GNSS data, the georeferencing matrices of the reference scans in the HD map are applied with a transformation so that the starting point of the trajectory coincides with the ground truth position of the CORS. Next, to distribute the GNSS observations on the simulated trajectory epochby-epoch, the following range modification term Δ s r,sim (t) is defined for receiver r tracking satellite s where || ⋅ || denotes the Euclidean norm. x r,gt (t) and x r,sim (t) represent the known coordinates of the CORS and the original ground truth position of the vehicle simulated at epoch t , respectively. x s (t) contains the coordinates of the satellites. Δ s r,sim (t) is added to the real-world code and phase observations so that they are simulated to have been collected by a GNSS receiver on-board the vehicle traversing the trajectory. Note that this simulation of the kinematic experiment does not consider negative effects in urban canyons that can degrade the positioning performance, such as signal blockage and multi-path (Lyu and Gao 2020).

Comparison of IGS precise products
In this section, we present the positioning results using static PPP comparing different types of IGS precise products to evaluate the accuracy of IGS RTS, since post-processed and real-time precise products can have distinct levels of performance in the case of vehicle positioning. We derive three types of precise corrections, namely Final, RT-30 s and RT-10 s. Final refers to the IGS final products, while RT-30 s and RT-10 s use identical real-time products from the SSRA01 stream sampled at different p . The orbital and clock sampling rates of these product types are summarized in Table 1. PPP positioning is performed using 260 1 h samples from the selected CORS assuming the receivers are stationary over time. Note that when precise products are missing for specific epochs due to low sampling rates, Lagrange multiplier is used to extrapolate the orbital or clock corrections. In addition, we define that PPP convergence is achieved when absolute ENU errors are all below 10 cm.
To verify the appropriateness of the chosen GNSS weighting functions, we begin by presenting the elevationdependency of the GNSS observational residuals. Figure 4 exhibits the formal values of squared residuals of the code and carrier phase observations following the inverse of (9) and their empirical counterparts by processing 3 h of (10) Δ s r,sim (t) = ||x r,sim (t) − x s (t)|| − ||x r,gt (t) − x s (t)|| GNSS data obtained from one of the used CORS, namely CUT0 in Perth, Australia, against satellite elevation angles ranging from 15° to 90°. It can be inferred that the squared residuals agree with the expectation and the GNSS measurement weight matrix is appropriate for approximating the minimum-variance property in the EKF estimation, provided that no mismodeled effect is present in the data (Teunissen 2000). Figure 5 shows the median values and the interquartile ranges of the absolute ENU errors of static PPP positioning results using Final products, with the left panel using o = 30 s and the right panel using o = 1 s. The classical static PPP performance is observed when o = 30 s that convergence is reached within 40 min (Teunissen and Montenbruck 2017). Surprisingly, when GNSS observations are retrieved more frequently at o = 1 s, the absolute ENU errors are significantly larger with higher dispersions and median convergence cannot be achieved within 1 h, as shown in the right panel. This can be attributed to the multi-path-induced temporal correlation of the observations, which is incorrectly ignored in the EKF time update that severely degrades PPP convergence, especially when the GNSS measurement interval is high (Zaminpardaz et al. 2017;Wang et al. 2018a;Erol et al. 2021). This phenomenon would be further amplified when the constant-velocity model is used instead of assuming a static receiver. However, in the case of vehicle positioning, frequently updated measurements are required (e.g., o = 1 s). In the next section, we will demonstrate the role of lidar for speeding up the convergence under such influence.
By maintaining o = 1 s, the positioning performances of using RT-30 s and RT-10 s products, respectively, are  Figure 6 exhibits the median values and the interquartile ranges of the absolute ENU errors of the static PPP positioning results based on the specified types of IGS RTS products, with the left panel using RT-30 s and the right panel using RT-10 s. Referring to the right panel of Fig. 5, despite having orbital corrections sampled at 900 s, Final products display higher precision and accuracy than RT-30 s, which is featured in the left panel of Fig. 6. The same behavior has been observed by Zhang et al. (2018), who concluded that orbital corrections can be extrapolated/ interpolated more accurately since their accuracy degrades linearly and hence have a lower requirement on the sampling frequency, whereas clock accuracy degradation varies with satellite clock types and requires higher sampling rates to enable accurate positioning. The RT-10 s positioning results in the right panel of Fig. 6 emphasize that by increasing p to 10 s, the positioning accuracy and precision are much improved as compared with the RT-30 s results and are even higher than the Final-product results, with the North and Up components converging at around 45 min. It can be concluded that the main advantage of RT-10 s is its high product sampling rate which compensates for the lower accuracy of IGS RTS compared with its final counterparts, making it suitable for real-time vehicle positioning. Therefore, RT-10 s products are used for the kinematic experiment discussed in the next section.

Simulated kinematic experiment results
To study the role of lidar for speeding up PPP convergence and realizing accurate and continuous positioning, the simulated kinematic experiment evaluates 3 positioning methods, namely PPP-only (standalone kinematic DF-PPP), Lidaronly using point cloud registration alone and the proposed PPP-Lidar integration. Moreover, the lidar availability ranges from 1 to 90 s for PPP-Lidar in order to quantify the influence of lidar on positioning and convergence performance. All the positioning methods discussed in this section utilize the constant-velocity model for the vehicle dynamics and the GNSS measurement interval of o = 1 s, while the RT-10 s products are used for those involving DF-PPP. Note that the Lidar-only results are based on 1 h of lidar data obtained from KITTI dataset, whereas PPP-only and PPP-Lidar results are summarized from 260 1 h samples of GNSS data utilizing the same lidar input.
The constants a and b in (7) are empirically determined as 0.03 and -1, respectively, by curve fitting of measurement residuals from independent data collected by the same lidar sensor model used in the experiment, which is Velodyne HDL-64E (Geiger et al. 2013). Among 3600 epochs of point cloud registration, by training the model on ETH dataset (Pomerleau et al. 2012), MS-SVConv successfully extracted 16 keypoints on average for 3105 scans obtained from the KITTI dataset, showing good transferability of the deep neural network. Keypoint extraction failed for the remaining epochs due to the HD map being downsampled to reduce the storage requirement. Moreover, the ENU spectral densities for the dynamic model in the EKF are configured as 0.05, 0.05 and 0.03 m 2 ∕s 3 , respectively, based on the mean ground truth acceleration of the vehicle (see Fig. 3).
To check the appropriateness of the intensity-based lidar weighting scheme, the formal squared lidar measurement residuals according to the inverse of (7) are compared with their empirical counterparts obtained from the Lidar-only positioning results with respect to the normalized intensity values of the keypoints. Figure 7 shows that the assumption of intensity-dependency of the lidar measurements is validated and highlights a pattern similarity with that of the GNSS weighting functions in Fig. 4.
We first present the Lidar-only (1 s) positioning results utilizing all available lidar keypoints (1 s indicates the lidar availability of 1 s, i.e., lidar input is available at every epoch). Figure 8 displays the time series of absolute ENU Lidar-only (1 s) positioning errors, with the red background highlighting the 495 epochs without lidar input, which solely rely on the EKF time update using the constant-velocity model for positioning. Despite having 43.9% of the epochs with cm-level 3D errors, the horizontal and 3D root mean squared error (RMSE) reach 0.76 and 0.80 m, respectively, as illustrated in Table 2, which shows the ENU, horizontal and 3D RMSE, as well as the median convergence time of the tested positioning methods. Large errors are observed when lidar measurements are absent for consecutive epochs due to registration failures as the accuracy of the time update degrades over time without the contribution of the measurement update.
In addition, PPP-only exhibits even poorer positioning performance. Figure 9 depicts the time series of the median absolute ENU errors of PPP-only and PPP-Lidar with different lidar availabilities, in which PPP-only has the largest positioning errors overall and does not converge within the tested period. The 3D RMSE of 1.55 m is the worst among all tested methods according to Table 2, which is insufficient  for vehicle positioning. However, it is evident from Fig. 9 that the East component is more accurate than the North component during the early stage of positioning, which can also be observed for PPP-Lidar results. This is due to the vehicle travelling mostly in the West-East direction for the first 20 min of the trajectory, creating larger variation in the satellite-receiver geometry in this dimension to improve the positioning precision over time.
In contrast, by combining the two aforementioned positioning methods, which individually showed low positioning accuracy, in the proposed integration, the PPP-Lidar (1 s) method utilizing all available GNSS and lidar measurements exhibits drastically improved positioning performance. Figure 9 depicts that PPP-Lidar (1 s) maintains the absolute ENU positioning errors well below the 0.1 m threshold for the entire tested duration. The horizontal and 3D RMSE are computed as 0.05 and 0.06 m, respectively, with immediate convergence within 2 s, according to Table 2, thereby, the convergence period of PPP is effectively eliminated. In addition, Fig. 10 displays the spatial and temporal distributions of the positioning errors using PPP, Lidar-only (1 s) and PPP-Lidar (1 s), respectively, by colorizing the median 3D errors. It is clear that PPP-only initially shows large errors over 1 m, despite slowly reducing them to decimeters. Although Lidar-only (1 s) contains the median 3D positioning errors to centimeters for most of the trajectory, they can reach m-level when lidar input is missing consecutively. In comparison, PPP-Lidar (1 s) utilizes both input in the EKF measurement update and successfully maintains cm-level accuracy throughout the 1-h duration.
It is important to note that PPP-Lidar (1 s) achieves such positioning performance thanks to the lidar measurements being available for 86.3% of the epochs. In reality, point cloud registration failure can occur more often due to outdated or inaccurate HD maps, thus decreasing the number of epochs with both types of measurements. Therefore, we   Fig. 9 and Table 2.
As the lidar availability decreases, the positioning errors increase accordingly. Table 2 suggests that from the lidar availability of 30 s onward, the ENU RMSE of PPP-Lidar reaches dm-level. However, the horizontal RMSE is always smaller than those of standalone PPPonly or Lidar-only (1 s). To maintain cm-level positioning accuracy in all directions, lidar measurements are only required once every 10 s, provided that 2 min are reserved for convergence. Figure 9 further shows that the RMSE values are mostly increased by the convergence period, which is shorter than 10 min for every lidar availability. After a brief initialization, even PPP-Lidar (90 s) can keep most of the East and North errors below 0.1 m. Note that the sudden drop of positioning errors around 9 min for all PPP-Lidar methods is caused by exceptionally large numbers of extracted keypoints comparing with other epochs, which are around 100. This reduction in PPP convergence time can also be supported by the fact that the float ambiguity solutions are considerably improved in precision by PPP-Lidar integration. For instance, in the case of CUT0 station tracking the GPS satellites G13 and G14 at the time of declared convergence, the standard deviation of the between-satellite single-differenced ambiguity solutions is improved by 5 times, namely, from 0.3 cycles to 0.06 cycles using PPP-Lidar (30 s).
Moreover, the horizontal and 3D errors of the kinematic positioning methods discussed in this section are summarized by way of cumulative distribution functions (CDF) in Fig. 11. Evidently, PPP-Lidar (1 s) outperforms all other methods since it utilizes all the available GNSS and lidar measurements. With less frequent lidar input, the integration always yields smaller positioning errors than PPP. Interestingly, the CDF curves of horizontal and 3D errors of Lidar-only (1 s) and PPP-Lidar (30 s) practically agree with each other for small error magnitudes until those of Lidar-only (1 s) escalate suddenly due to large errors caused by missing lidar input. This suggests that by fusing with PPP using RT-10 s products, the integration gives similar positioning performance as Lidaronly (1 s), which has been shown to be acceptable when lidar is present, by requiring only 1/30 of the lidar measurements. In other words, by using the proposed method for accurate vehicle positioning, the density of reference scans in the HD map can be dramatically relaxed, making it less costly to develop and deploy.
In addition, we recognize that accuracy requirements and possible lidar availabilities can vary with the types of positioning applications. Therefore, we provide an analysis of the convergence time to satisfy certain accuracy requirements using PPP-Lidar with different lidar availabilities. This is based on the same kinematic positioning results obtained above and shown in Fig. 12. For stringent accuracy conditions such as 5 cm or 10 cm in both vertical and horizontal directions, respectively, or 10 cm of 3D accuracy, the proposed method benefits more when the lidar input is more frequent. Nevertheless, the median values of the convergence time are always below 10 min. Moreover, instantaneous convergence is consistently achievable for more relaxed requirements such as 1 m vertically and 0.5 m horizontally, and 0.5 m in 3D, while a higher lidar availability results in a smaller dispersion of the convergence time, i.e., more reliable rapid convergence. These statistics highlight the much improved convergence behavior of the PPP-Lidar integration and provide further guidelines on the configurations of lidar and/or HD maps to realize accurate and continuous vehicle positioning.

Concluding remarks
In this contribution, we proposed an observation-level integration of GNSS PPP and lidar for accurate vehicle positioning with a significantly reduced convergence requirement. Specifically, lidar measurements were generated as keypoints by registering point clouds with a prebuilt HD map using deep learning, before being combined with DF-PPP observations corrected by real-time precise products. The integration was realized in an EKF employing the constant-velocity model that captures the vehicle's motion, in which lidar measurements were weighted based on their intensity and geometric distributions. Moreover, a simulation strategy to couple raw GNSS and lidar data collected from different campaigns for kinematic positioning experiments was designed, through which the proposed positioning method was evaluated, in addition to the IGS precise product comparison conducted in the static mode.
The static experiment using GNSS data collected from 20 Australian CORS showed that IGS RTS precise products offer positioning accuracy comparable with that using final products, provided that the product sampling rate is higher (e.g., p = 10 s), and are suitable for real-time positioning applications. By using the proposed PPP-Lidar method for kinematic positioning with RT-10 s, in which IGS RTS products are sampled every 10 s, PPP convergence can be reached instantaneously, while maintaining cm-level positioning accuracy throughout the experimental duration and giving the 3D RMSE of 0.06 m. In pessimistic cases where lidar availability was reduced, e.g., when lidar keypoint extraction fails due to faulty HD maps, lidar input was only required in the proposed integration once every 10 s for reliably constraining the positioning errors and convergence time to cm-level and approximately 2 min, respectively. Recommendations for lidar availabilities required for satisfying different levels positioning accuracy were also given through analyses, showing that convergence can be achieved within 10 min for all tested configurations, while more frequent lidar input leads to more reliable rapid convergence.