Enhanced trajectory estimation of mobile laser scanners using aerial images

Multipath effects and signal obstruction by buildings in urban canyons can lead to inaccurate GNSS measurements and therefore errors in the estimated trajectory of Mobile Laser Scanning (MLS) systems; consequently, derived point clouds are distorted and lose spatial consistency. We obtain decimetre-level trajectory accuracy making use of corresponding points between the MLS data and aerial images with accurate exterior orientations instead of using ground control points. The MLS trajectory is estimated based on observation equations resulting from these corresponding points, the original IMU observations, and soft constraints on the pitch and yaw rotations of the vehicle. We analyse the quality of the trajectory enhancement under several conditions where the experiments were designed to test the influence of the number and quality of corresponding points and to test different settings for a B-spline representation of the vehicle trajectory. The method was tested on two independently acquired MLS datasets in Rotterdam by enhancing the trajectories and evaluating them using checkpoints. The RMSE values of the original GNSS/IMU based Kalman filter results at the checkpoints were 0.26 m, 0.30 m, and 0.47 m for the X, Y-and Z-coordinates in the first dataset and 1.10 m, 1.51 m, and 1.81 m in the second dataset. The latter dataset was recorded with a lower quality IMU in an area with taller buildings. After trajectory adjustment these RMSE values were reduced to 0.09 m, 0.11 m, and 0.16 m for the first dataset and 0.12 m, 0.14 m, and 0.18 m for the second dataset. The results confirmed that, if sufficient tie points between the point cloud and aerial imagery are available, the method supports geo-referencing of MLS point clouds in urban canyons with a near-decimetre accuracy.


Introduction
Multipath effects and signal obstruction by buildings in urban canyons can lead to inaccurate measurements with Global Navigation Satellite Systems (GNSS) and therefore errors in the estimated trajectory of Mobile Laser Scanning (MLS) systems. Kukko (2013) demonstrated that GNSS measurement accuracy can degrade to>50 cm during an outage of GNSS signals. In this case, acquired point cloud quality suffers from an inaccurate trajectory and the 3D data becomes less useful for mapping applications. Under ideal circumstances, without any GNSS signal outage or multipath, a state-of-the-art Mobile Mapping (MM) platform can achieve 2-3 cm accuracy, estimated by Haala et al. (2008); Kaartinen et al. (2012). However, this is not possible in urban canyons. The task of trajectory correction is crucial to ensure the quality of the mobile laser scanning point clouds. Commercial software first tries to correct the trajectory by automatic registration of multiple passes, similar to the techniques reported by Levinson et al. (2007), Ding et al. (2007); Zhao (2011). In the same vein, Hunter et al. (2006) and Bornaz et al. (2003) proposed a consecutive strip adjustment to improve the misalignments in MLS data sets. Similarly, Bosse and Zlot (2009) described a scanmatching method based on iterative closest points (ICP) to recover an accurate trajectory. However, during long-term GNSS signal outages, the achieved accuracy is still metre-level, as expressed by Chiang and Huang (2008). Furthermore, all of these techniques can only be used to increase relative accuracy and require data that has multiple scans/ passes of the same scene. To minimise costs it is, however, desirable to scan an area only once.
Normally, after automatic registration, Ground Control Points (GCPs) are collected in the target area and then their correspondences in the MLS point cloud are manually selected. Although some software does provide assistance to automatically detect such landmarks, the final decision and effort remain with a human operator. Finally, the erroneous trajectory is adjusted based on the manually established correspondences and the MLS point cloud is then regenerated. However, manual acquisition of the GCPs is very labour intensive, costly and errorprone. Therefore, an automated method is desirable to replace manual correction, particularly to improve the MLS platform's trajectory.
We make use of well-oriented aerial imagery as the georeferencing source for the MLS data. In many countries such imagery is available at a national level, hence, costs are only incurred for the joint processing of the MLS data and aerial imagery, not for data acquisition. The whole workflow for trajectory improvement is comprised of two major steps: The first step is registration of the MLS point cloud with aerial images to determine corresponding points. This has been addressed in our earlier work (Hussnain et al., 2019), briefly summarised in the related work section. The second step, the main contribution of this article, is the automatic 6 degree of freedom (DOF) trajectory adjustment in which we utilize three types of observation equations: (1) those resulting from the established correspondences, (2) Inertial Measurement Unit (IMU) observations, and (3) soft constraints to the vehicle pitch and yaw rotations. We demonstrate that we nearly reach a decimetre accuracy which was set as a goal in discussion with four mobile mapping companies operating in the Netherlands. This accuracy was considered sufficient for most of the mapping applications of their clients. We furthermore analyse how the trajectory accuracy decreases with reduced numbers of correspondences to aerial imagery or a temporary unavailability of such correspondences during which the trajectory is estimated on the basis of IMU data only. We model the six pose parameters of the trajectory with B-spline functions. For modelling with sufficient accuracy, the optimal settings of the curve order and knot interval are determined.
After a description of related literature in Section 2, Section 3 describes our method of the trajectory adjustment by introducing the various observation equations. Section 4 explains the design of the experiments and Section 5 discusses the results obtained with datasets acquired with two different mobile laser scanners and demonstrates a near-decimetre accuracy can be obtained, even by a mobile laser scanner with a relatively low quality IMU sensor.

Related work
Several studies have been performed to improve the georeferencing of mobile mapping sensors by matching the recorded data to another already properly georeferenced dataset. Most relevant to our research is the work presented by Gao et al. (2015) in which the trajectory of a mobile laser scanner is adjusted based on correspondences between the point cloud and imagery. Gao et al. (2015) employ a shape-preserving piecewise cubic Hermite interpolation to estimate a time-dependent offset to be applied to the original MLS trajectory and improve RMSE values in X-, Y-and Z-coordinates from respectively 0.13 m, 0.21 m, and 0.47 m to 0.09 m, 0.06 m, and 0.10 m. High resolution UAV imagery is required for this purpose. No use is made of IMU measurements and sensor attitude estimates are not updated. Hence, the shape of the trajectory in between subsequent tie points to the aerial imagery is still dependent on the original trajectory solution. Wolcott and Eustice (2014) developed an image-based navigation for self-driving systems. The mobile mapping camera images were registered with previously acquired 3D lidar data by maximizing the normalized mutual information method. Their developed approach achieved an RMS error of 0.19-0.45 m in longitudinal direction and 0.14-0.21 m in lateral direction. Kümmerle et al. (2011) developed a graph-based SLAM procedure in which correspondences between laser range scans and aerial images were formulated as constraints. Edges extracted from the aerial images were matched against point clouds captured on a robot by a 2D laser scanner which was slightly tilted up and down with respect to the horizontal plane. The scanner was primarily used to support the 2D navigation as it was continuously capturing vertical structures for the SLAM. Compared to a 2D cadastral map an RMSE of 0.2 m was obtained on the six building corners used as checkpoints.
Recently, Javanmardi et al. (2017Javanmardi et al. ( , 2018 proposed a technique for MLS platform localization based on the purpose-made HD maps generated from an accurate point cloud recorded before the MLS survey. Their approach constrained the registration problem to 2DOF, reported accuracies ranged from 0.1 to 0.5 m.
Surprisingly, none of the above approaches takes advantage of IMU data to improve the trajectory estimate and to bridge trajectory parts without correspondences to the data source used for georeferencing.
In our work we make use of B-splines to model the six pose parameters of the MLS over time. B-splines are frequently used for trajectory modelling as they allow local updates of optimal trajectories if additional control points become available and are able to represent the trajectory as a continuous function over time. Usenko et al. (2017) used B-splines to re-plan trajectories for micro aerial vehicles when previously unknown obstacles were detected. The main purpose of this study was to show that the developed system could accommodate real-world dynamics into the trajectory planning and did not assess achieved positioning accuracy. Patron-Perez et al. (2015) exploited benefits of B-splines for egomotion estimation with unsynchronised sensors including rolling shutter cameras which involved individual time stamps for each pixel. Results were only presented for simulated datasets and the accuracy and reliability of the developed system was not discussed. Furgale et al. (2012) as well as Ovrén and Forssén (2018) used Bsplines to model trajectories for SLAM based on visual-inertial odometry. Both noted that first and second B-spline derivatives of pose parameters can be elegantly linked respectively to the angular velocities and accelerations measured by IMUs. Unaware of their work at that time, we formulated similar relationships in our previous work (Hussnain et al., 2018). The same IMU observation equations were also used to support the estimate of the trajectory of our backpack laser scanning system in the case of insufficient matches between points seen by the laser scanners and the SLAM map (Karam et al., 2019).
In our method for trajectory adjustment we make use of corresponding points between the MLS point cloud and aerial imagery. Our matching procedure as well as a review of other work on matching point clouds to images has been described in (Hussnain et al., 2019). To make this article self-contained, we briefly review our procedure to establish correspondences between the point cloud and aerial imagery. Based on the trajectory information we split the point cloud into tiles. The orientations of the aerial images are used to select point cloud tiles showing the same piece of the terrain and to project the points of those tiles back into the aerial images. The reflectance strength of the MLS points are used to create a reflectance strength image in the coordinate system of the aerial image. Key points are extracted with the Harris operator (Harris and Stephens, 1998) and matched on the basis of the LATCH descriptor (Levi and Hassner, 2016). The same key points are also matched between overlapping aerial images. For each key point in the point cloud with corresponding key points in multiple aerial images we triangulate the aerial image key points to obtain a 3D point in the terrain coordinate system. The objective of the trajectory adjustment will be to minimize the distances between these triangulated points and the corresponding key points in the MLS point clouds. We will refer to these points as 3D tie points in the remainder of this article. For further details on the matching procedure, please refer to (Hussnain et al., 2019).

B-spline based 6DOF trajectory adjustment
In this paper, we further develop and analyse the trajectory adjustment procedure initially formulated in Hussnain et al. (2018). The workflow of the adjustment procedure is presented in Fig. 1. The developed workflow employs 3D tie points between the point cloud and aerial images, IMU measurements, and soft constraints on the trajectory. The trajectory is parameterized with six B-splines modelling the six pose parameters of the MLS as a function of time. The GNSS/IMU based Kalman filtering solution 1 is used as the approximation of the trajectory to be improved in the adjustment. Although the Kalman filter solution will be affected by multi-path and occlusion of GNSS signals, errors in the trajectory are typically below 1-2 m and hence the Kalman filter solution is accurate enough to obtain initial values of all spline coefficients. The workflow incrementally improves the trajectory estimates until it converges.
In the remainder of this section we first discuss the required parameters for the representation of trajectory with B-splines. We then describe the observation equations obtained from 3D tie points, IMU measurements of angular velocities and accelerations and soft constraints.

B-spline order and knot interval optimization
To model the MLS trajectory with B-splines of sufficient accuracy we need to estimate the optimal values of the knot interval and the curve order. The B-spline with optimal parameters only needs a minimum number of coefficients, otherwise, for every additional knot interval an extra spline coefficient needs to be estimated for each of the six pose parameters. Because we want to achieve decimetre-level accuracy in the improved MLS point cloud we should ensure that the errors introduced by modelling the trajectory with splines are significantly smaller. We set the maximum position error introduced by the spline modelling to 4 cm.
For the rotational error the criterion is also based on the maximum positional error allowed. Suppose a mobile mapping car observes a 3D point 20 m away (based on an average road width in the datasets), then the maximum error allowed in any angle is θ ε = 0.12 • if the effect of the error on the coordinate calculation should remain within 4 cm as shown in Fig. 2.  1 While we refer to this solution as the GNSS/IMU based Kalman filtering solution, mobile mapping companies also typically apply a backward pass with the Rauch-Tung-Striebel smoother to use all measurements to estimate all state vectors. In addition constraints are imposed wherever the trajectory intersects itself.
In Section 5 we analyse different combinations of spline order and knot intervals to find the splines with the minimum number of parameters which still fulfil the accuracy requirement.

3D tie point observations
The correspondences between aerial images and the MLS point cloud are the most important input to our adjustment procedure because the aerial images are the georeferenced source.
In order to linearize the 3D tie point observation for the B-spline based adjustment we start with a point in the recorded MLS point cloud. A laser scanner observes a 3D point X C PC in the car coordinate system, which can be rotated to the world coordinate system with a rotation matrix R W C representing the attitude of the car in the world coordinate system. The point X C PC can then be translated to the world coordinate system by adding the location of the car T W C in the world coordinate system. This relation is represented in Eq. (1). The upper indices C and W indicate the coordinate system of the car and world respectively. The lower indexing specifies the source of the coordinate vector, e.g. AI for aerial image and PC for the point cloud.
Here X W PC is the point cloud point as observed in the world coordinate system. R W C (t) and T W C (t) represent rotation and translation which change over time t.
Because the GNSS data was unreliable, we want to re-estimate the rotation R W C (t)and translation T W C (t) based on the 3D tie point X W AI obtained from the aerial images and the corresponding point in the point cloud (Fig. 3). This relation is represented in Eq. (2).
The original GNSS and IMU data have been processed by a Kalman filter to estimate the rotationR W C,Kalman (t) and translation T W C,Kalman (t) between the car and the world coordinate system. The Kalman filter results have been used to reconstruct the original point cloud points X C PC in the car coordinate system as in Eq. (3).
The result of the Kalman filtering is used to obtain the approximate spline coefficients of the six pose parameters over time. The six pose parameters comprise three angles ω(t), φ(t), κ(t) and a translation vector described by translations along the three axes, T X (t), T Y (t), T Z (t) The modelling of e.g. ω(t) by a B-spline function is given in Eq. (4).
Where B i is i'th B-spline function and α ω,i its coefficient to be estimated. Eq. (2) is linearized with the upper index 0 denoting the approximate value. In the first iteration, the output of the Kalman filter is used to obtain the spline coefficients for the approximate rotation and translation. The time dependency (t)is omitted below to shorten the expression. The values of Euler angles are unwrapped to handle the jumps from − 180 to 180 degree, which is necessary for the estimation of smooth spline functions.
At the left-hand side is the observed misclosure between the 3D tie point obtained from the aerial images and the 3D point obtained in the point cloud expressed in the world coordinate system. On the right-hand side are the increments to the spline coefficients multiplied with the elements of the Jacobian.
The 3D-3D correspondences are not available consistently over the whole trajectory. Instead, their availability is dependent on the presence of the road markings. Moreover, even if available, tie points are sparse and cannot cover every spline interval. Therefore, the positioning has to be supported by inertial navigation.

Acceleration observation
The IMU observes accelerations for the car's position in the sensor coordinate system S. These are denoted Ẍ S IMU . After rotation to the car coordinate system by the rotation matrix R C S and then to the world coordinate system by R W C , and subtraction of the gravity vector, these accelerations should correspond to the second derivatives of the car's location in the world coordinate system, i.e.T W C . Hence, The angles of the rotation matrix R C S are calibrated and remain constant during the laser scanning operation.
To model the second derivative of the car locationT W C , it is straightforward to take the second derivative of spline as it is a polynomial function, like for T X (t) : The B-spline coefficients α TX ,i to be estimated remain exactly the same. It's only the B-splines functions themselves that need to be differentiated twice. Also, note that the differentiation only applies to the B-splines of the translation and not to those of the rotation angles inR W C .
With this the linearized equation becomes:

Angular velocity observation
The IMU also observe angular velocities in the sensor coordinate The first derivatives of the angles describing the rotation from the car to world coordinate To determine the relationship between the observed angular velocities, we first need to define the order and direction of rotationR W C from the car coordinate system to the world coordinate system, where components of the rotation matrix R W C are defined as; The rotation from the world coordinate system to the car coordinate system is therefore defined by: The ω W C is the first rotation applied when rotating from the car to the world coordinate system, therefore the angular velocity ω W C of the car in the world coordinate system is only observed around the X-axis of the IMU. The first derivative of car's rotation around the Y-axis in the world coordinate system does not directly correspond to the rotational velocity around Y-axis in car coordinate system because Y-axis of the car has been rotated by − ω W C around the X-axis. Therefore, the derivative of φ W C should be rotated to the car coordinate system by the rotation matrixR 1 ( . Similarly, the derivative of κ W C needs to be rotated around the Y-axis by R 2 (− φ W C ) and then X-axis by R 1 before it corresponds to the angular velocity vector in the car coordinate system.
which can be written to define the angular velocity observation equation, ⎛ where ω is a vector of angular velocities around X-, Y-and Zcoordinates.
In order to model the angular velocities, we use the first derivatives of the B-splines function, e.g. for ω(t) Linearization of Eq. (9) leads to where the partial derivate of S C W w.r.t. ω and φ are

Soft constraint observation
The direction of the car is described by the heading κ W C and pitch φ W C , but these two rotations can also be inferred from the first derivatives of the car's positions T W C as shown in Fig. 4. Therefore, we can add two constraints to ensure that the rotations are consistent with the trajectory, effectively reducing the degrees of freedom in the transformation from the car to the world coordinate system from six to four.
The heading can be inferred from the first derivatives of the X-and Ycoordinates of the car trajectory. As, the first rotation from the world coordinate system to the car coordinate system was defined -κ. Hence, the soft constraint observation related to the heading of car is, Fig. 4. The relationship between positions and direction of the car.
To shorten the notation below, first derivatives are denoted where Ṫ 0 X and Ṫ 0 Y are the first derivatives of the X-and Y-coordinates derived from the Kalman filtered trajectory for the first iteration of the adjustment process.
Similarly, the pitch angle of the car in the world coordinate system can be inferred from

Design of the experiments
In this section we outline the design of experiments aiming to assess the trajectory after the adjustment. As multiple inputs take part in the adjustment process various experiments are needed to analyse the effect of individual input on the accuracy. This section proposes several experiments to determine the minimum amount and quality of the tie points and the maximum time without tie points which still allows an accurate estimation of the trajectory. The proposed set of experiments are summarised in Table 1. The description of these experiments is provided in the next sections. Experiment 0 is the original trajectory as obtained from the GNSS/IMU based Kalman filtering which is used as a baseline for the adjusted trajectories making use of the different types of observation equations.

Datasets
For the experiments we use two independently acquired trajectory datasets; Trajectory-I and Trajectory-II as plotted in Fig. 5 and Fig. 6 respectively. Both systems are based on different laser scanning and position estimation systems. The Trajectory-I dataset is acquired by Topcon® IP-S3 mobile mapping system, which is a single 360-degree lidar sensor. The poses are estimated by an industrial grade IMU (KVH® CG-5100). The length of Trajectory-I on the road surface is 13  Z. Hussnain et al. km with the total scanning time of 56 min. Only considering the points from the single pass, the mean density of points on the road surface is 3.2 points/cm 2 . Trajectory-II is acquired by Topcon® IP-S2 system that is based on two 360-degree laser scanners mounted in a cross-configuration and a Honeywell HG1700 AG58 IMU. The length of the trajectory on the road surface is 10 km with the scanning operation lasting 45 min. Only considering the points from the single pass and single scanner, the mean density of points on the road surface is 0.74 points/cm 2 . Both used mapping systems employed a comparable dual frequency GNSS receiver.
The technical specification of the IMUs of both systems are summarised in Table 2. The second dataset is more challenging for various reasons: the IMU performance specifications are clearly lower, but the IMU used for trajectory II was also known to be degraded because of frequent use. Furthermore, the buildings next to the trajectory are taller than those along trajectory I, leading to more occlusions of GNSS signals.
Fifteen aerial images of 20010x13080 pixels were captured with the UltraCam® camera on a manned flight at the approx. height of 4500 m. The quality of the orientation of the aerial images has been assessed by the Pix4Dmapper software using checkpoints and shown to be accurate to 10 cm corresponding to 1 GSD. The aerial images have 60% forward overlap and 40% across track overlap. The orientation parameters of both aerial images and MLS datasets are known in the Dutch reference system RD-NAP.
A total of 50 corners of road marks or potholes were surveyed with RTK with a few cm accuracy. As the point clouds had a point spacing of around 1 cm, the checkpoints could also be accurately located in the point clouds colour-coded by the pulse reflectance strength (e.g. Fig. 8) and hence be used to assess the accuracy of the estimated trajectories.

Quantitative analysis for the tie point observation
Experiments 1, 2, and 3 were used to determine the minimum number of tie points needed to accurately correct the trajectory. This was achieved by performing the adjustment with gradually decreasing quantities of the tie points and then evaluating the accuracy of each adjustment with checkpoints which have been surveyed in the terrain by GNSS and corresponding locations are hand-picked in the MLS point cloud.
The 2nd experiment reduced the tie points set by removing tie points from randomly selected locations. The reduction could therefore also eliminate the tie points in areas where they are already scarce. Thereby, the experiment analysed the situation where some segments provided a more than sufficient amount of tie points while the others provided none or just a few. This lead to a situation where some locations in the trajectory may not be improved as the trajectory reconstruction in these parts will completely depend on the IMU.
In the 3rd experiment we removed tie points such that the remaining points were equally distributed as much as possible along the trajectory. For this purpose, we first calculated the distribution of the tie points along the trajectory over time. Then the tie points were grouped in bins w.r.t. the interval size that is optimized in Section 5.2 and then the bins containing more tie points were reduced first to get a more even distribution. This strategy ensured that the reduced tie points remained well distributed along the trajectory.

Qualitative analysis for the tie point observation
Experiments 4a-4d were conducted to analyse the relationship between the quality of the tie points and that of the reconstructed trajectory. For each tie point we measure the misclosure between the tie point as obtained by triangulation of the key points in the aerial images and the corresponding point in the MLS point cloud (Eq. (5)). The accuracy of this observed misclosure primarily depends on the measurements in the image as the pixel size of 10 cm is much larger than the point spacing in the MLS point cloud. Even when road marks are more located towards to road side the average point spacing is below 2 cm. We therefore assess the quality of each tie point by error propagation in the triangulation of  the aerial image key points. By only using the tie points with a certain range of standard deviations we can analyse the relationship between the quality of the tie points and the accuracy achieved in the trajectory after the enhancement. The experiments 4a-4d required a trajectory containing tie point sets of varying standard deviations. Based on the standard deviation we can divide tie point sets into low, medium and high quality. Normally, the high-quality tie points are computed from more than two aerial images and entail small back-projection error. Moreover, the high-quality tie points are the result of triangulation of subpixel accurate 2D key points which are the best estimation of 2D corner features in the aerial image whereas the low-quality tie points have the opposite characteristics. As indicated in Table 1, we separately enhanced the trajectory using tie points of high-quality (4c) with standard deviation ≤ 5 cm, mediumquality (4b) with standard deviation > 5 cm and ≤ 8 cm and lowquality (4a) with standard deviation > 8 cm, and then compared the results to establish the minimum level of quality needed. We used all tie points with weights based on their standard deviations in experiment 4d.

Assessment of trajectory constructed only by IMU and soft constraints observations
The IMU observations provide inertial navigation whereas the 3D tie point observations provide a geo-referencing relative to the aerial imagery. However, in locations where 3D tie points are not available, the IMU-based positioning suffers from the accumulation of noise over time and drift errors in the inertial measurements. Even when 3D tie points are available, in the construction of the entire trajectory, there are small segments of various lengths that are reconstructed based on IMU observations only. In the 5th experiment, we prepared segments of stepwise increasing lengths from the same trajectory covered in 30, 60 and 90 s. Moreover, we fixed the start and end pose of the trajectory segments because not fixing the start end pose results in a singular equation system. The fixed pose observation is explained in Hussnain et al. (2018). An illustration of the start and end poses of three trajectory segments and their unified formation is shown in Fig. 7. Between each start and end pose the trajectory is constructed only based on the IMU and soft constraint observations. Fixing the pose at both ends is important because it forces the trajectory to pass through the start and end poses, this situation is compatible with the real world scenario where an interval of no tie points is always adjoined at both ends with the segments containing tie points. Otherwise, if we only fix the pose at the start, it will increase the drift error near the end pose.

Impact of soft constraints
Experiment 6 has been used to distinguish the improvement ascribed to the soft constraints by not applying these constraints and then comparing the result with the same trajectory which has already been improved by the soft-constraints as in experiment 1. All of the above defined experiments employ the soft constraints because it is not necessary to remove them while analysing the influence of the other observations on the adjustment. Therefore, we also use soft constraints for the IMU based trajectory in Section 4.3 because it focuses on exploiting other available observations except those of the 3D tie points.

Results
In this section we present the results of the experiments described in the previous section starting with the criteria of the trajectory segment selection followed by the optimization of the B-spline knot interval and order. Next, the procedure to evaluate an updated trajectory using checkpoints is described. Lastly, the main evaluations of the experiments are presented where each experiment is associated with the analysis of individual observations.

Trajectory selection criteria
To conduct the experiments, we needed to select a part of a trajectory where most subparts are supplied with the road marks and thereby the checkpoints are attainable. The trajectory can be considered as a collection of subparts or segments where each segment is related to a point cloud tile and its associated tie points. The extraction of tie points related to a point cloud tile and its registration to aerial images is described in Hussnain et al. (2019). When all consecutive segments of a trajectory contain tie points it is feasible to evaluate and ensure that a certain accuracy persists in the whole trajectory. Moreover, it provides reliable and accurate start and end poses, which are also needed for the IMU-based observation. For the experiments in this section we select parts of the trajectories which fulfil these criteria.

Optimal knot interval and B-spline order
We needed to understand how well a trajectory could be modelled by B-splines. Therefore, we took a trajectory as processed by the Kalman filter, sampled it densely by extracting points at every 10 ms and used those points to reconstruct the trajectory with B-splines. These reconstructed trajectories were then compared against the original trajectory to analyse the errors introduced by the approximation of the trajectory by B-splines. For the optimization, various combinations of knot interval and curve order parameters were used. Based on preliminary experiments on the spline fitting, testable ranges of knot interval and curve order were selected. As a result, three curve orders 2, 3 and 4 in combination with five knot intervals of respectively 0.1, 0.2, 0.5, 1.0 and 2.0 s were constructed. Each combination of parameters sets was used  separately to fit B-splines to the six parameters X, Y, Z, ω, φ and ϰ. The goal was to identify the combination of largest possible knot interval with a smallest possible curve order which satisfies accuracy criteria, 4 cm for the position and 0.12 degrees for the rotation angles as mentioned in Section 3.1.
The dataset of Trajectory-I was used for this experiment, results are provided in Table 3. RMSE values not meeting the requirements are highlighted in bold. Those for the knot intervals 0.1 and 0.2 s have been omitted as these were clearly too short. The results in Table 3 show that the criteria for the positional splines can be more easily met than those for the rotation angles. With a knot interval of 2.0 s none of the tested orders could meet the 0.12 degrees. Results for the knot interval of 1.0 s showed that the spline order did not have a large influence on the RMSE of the angle splines. The worst results were obtained for the φ spline. For all three orders the RMSE was around 0.12 degrees, be it slightly above for orders 2 and 3. To stay on the safe side, order 4 in combination with a knot interval of 1.0 s was selected for all experiments.

Regeneration of point cloud and 3D tie points
The main motivation for the trajectory adjustment is to reconstruct an accurate MLS point cloud. The point cloud in which the 3D tie points have been measured has originally been computed with the trajectory resulting from the Kalman filtering. Hence, a point in the world coordinate system X W PC was computed based on the observed point X C PC by the laser scanner in the car coordinate system with Eq. (16).
where the rotation matrix R W C,Kalman and translation vector T W C,Kalman are the inaccurate pose parameters from the GNSS/IMU based Kalman filter solution. So, the measured point X W PC is the inaccurate point cloud point, which needs to be recalculated based on the updated pose parameters after the trajectory enhancement.
Assuming that the trajectory pose after the adjustment corresponds where X W ' PC is a 3D point of the improved point cloud in the world coordinate system. This is also how the improved position of the 3D tie points will be computed. Note that the local laser scanner observation X C PC remains the same, comparing Eq. (16) with Eq. (17). It is only the pose parameters R W ' C (t)and T W ' C (t) that are updated. For every experiment we first select a trajectory and then adjust it using the combination of observations stated in Section 4. Then we verify the accuracy by regenerating the tie points in the MLS point cloud from the adjusted trajectory and measuring the improvement compared to the checkpoints for which two sets of corresponding checkpoints are acquired and labelled. Set A contains 38 checkpoints and set B contains 12 checkpoints. As an evaluation example the adjustment of Trajectory-I using all observation is conducted and the residuals between the regenerated tie points and checkpoint A47 and checkpoint A50 are provided in Table 4. Plots showing the checkpoints before and after the adjustment are presented in Fig. 8 and Fig. 9 for the checkpoint A47 and A50 respectively. The residual of A47 is comparatively larger than A50 on XY plane, however, considering the Z coordinate there is more error in A50 than in A47. Moreover, calculating the norm e.g. 23 cm for A47 and 18 cm for A50 shows that the trajectory has been adjusted comparatively better for A50 than for the A47.

Experiments and evaluation
We performed the experiments according to the guidelines prepared in Table 1. We conducted the experiment on Trajectory-I and Trajectory-II datasets.
We selected two subparts Trajectory-IA and Trajectory-IB from the Trajectory-I which are plotted in Fig. 10 and Fig. 11 respectively. These subparts fulfilled the criteria for the experiments as mention in Section 5.1. The first subpart was used for experiments 1, 2, 3 and 4 because these experiments could be conducted with infrequent checkpoints. As Table 3 Results for B-spline fitting to Trajectory-I dataset using combinations of knot interval and curve order. Results not meeting the criteria are shown in bold.

B-spline order
Knot interval (sec)  Table 4 Two examples of residuals measured using checkpoints on the point cloud regenerated from the trajectory enhanced using all observations.  Z. Hussnain et al. these experiments were based on the tie points we needed to select the trajectory part which had road marks at some places. For experiment 5 we needed to select a trajectory with abundant and consistent availability of checkpoints regardless of the road marks. For this purpose, the second subpart was feasible because it had more checkpoints to test the IMU-based trajectory over small intervals without the involvement of the tie points.
To conduct the experiment, first, we performed the adjustment of a particular trajectory using the observations mentioned in Table 1. Then the evaluation was performed using the regenerated point clouds around the locations of the checkpoints. Table 5 presents the RMSEs at the checkpoints as well as the minimum and maximum error.
The experiment 0 in Table 5 was conducted to analyse the accuracy of the Kalman trajectory before any observations were used for the adjustment. This experiment provided an assessment of the current state of the Trajectory-I dataset. Before analysing the relation between the quantity of the tie points and the accuracy of the trajectory, experiment 1 applied all tie points and other observations to provide the reference for further analysis. It results in the accuracy of the trajectory without reducing the number of tie points.
In experiment 2 and 3 in Table 5, we tested the impact of a reduced number of 3D tie points. In experiment 2 we removed the tie points situated in consecutive order. This operation lengthened the trajectory intervals without tie points where tie points are scarce. Dissimilarly in experiment 3, the total number of tie points were reduced from where they are already in majority. This removal technique circumvents the wide gaps of no tie points, therefore, we expected to obtain better results than for experiment 2. For evaluation of the trajectory produced by experiments 2 and 3, we used 14 checkpoints. By comparing the results of experiments 2 and 3 we notice that while reducing the tie points, the accuracy considerably suffers in experiment 2. As with the wide gaps of no tie points the IMU is the only source to construct the trajectory. However, it is not directly comparable to experiment 5 where we only used IMU observation without any tie points to test the reliability of IMU for small intervals.
In experiment 3, the accuracy does not deteriorate when applying the 80% tie points because the remaining tie points will nonetheless provide observation for approximately all locations. However, when we use only 50% tie points in experiment 3 the accuracy approaches an undesirable level.
In the 4th experiment it is clearly noticeable that a decimetre-level accuracy can only be obtained with the high-quality tie points, either by only using those tie points (experiment 4a) or by giving lower weights to the low accuracy tie points (experiment 4d). Results of these two experiments are quite similar. Conversely, only using the low-quality tie point leads to large errors in the trajectory. Yet, the utilization of medium or even low-quality tie points is better than using no tie points at all which is apparent from comparing the results of experiment 4 with experiment 5. Even for a small interval of the 120 s without tie points the accuracy is poorer than using the low-quality tie points alone. Note that the low-quality tie points are used for the whole trajectory sustaining 45 min but still, the accuracy is comparatively not as poor as the IMU-based trajectory exhibit for the small intervals. The high-quality tie points are only 55% of the total set yet they are alone enough to correct the trajectory to the maximum accuracy. The result of experiment 3 with 80% tie points and high-quality tie points are comparable but the accuracy is little decreased because 20% removed tie points set comprise some highquality tie points. Moreover, the remaining 80% tie points can comprise low-quality tie points, therefore, the accuracy achieved in experiment 3 with 80% tie points is still less than in experiment 4c. Furthermore, a comparison of the result of experiment 4c with experiment 1 confirms that the use of low-quality tie points introduces error in the trajectory. It is therefore advisable to only perform the trajectory adjustment with high-quality tie points.
In experiment 5, the trajectory was constructed for the interval of 30, 60 and 90 s which is scanning distance of 396, 792 and 1101 m respectively. We only checked the trajectory for the interval of up to 90 s because the resulted accuracy was already degraded more than the desired level. The main reason to conduct experiment 5 was to determine the maximum tolerable time of having no tie points. The largest interval of no tie points in the Trajectory-I lasts for 48 s, which is roughly 556 m of distance on the road. This interval arrives between the range of 30 and 60 s so we expect the accuracy to be better than indicated in 60 s interval, while the accuracy of 60 s interval is a little above the decimetre level. Therefore, we expect the accuracy to be below or near the decimetre level during the 48 s interval of no tie points.
In experiment 6, we evaluated the improvement in the trajectory without using the soft constraint observation. Comparing this result with experiment 1 confirms that the soft constraint observation is valuable. Without using them the RMSE values increased by 2, 1 and 5 cm for the X-, Y-and Z-coordinates respectively.
The Trajectory-II is in the area with the tallest buildings in Rotterdam and contained 21 checkpoints as plotted in Fig. 12. We repeated the six experiments for this trajectory as well, the results are provided in Table 6.
Overall the Kalman filter estimate of Trajectory-II is less accurate than that of Trajectory-I (cf. Table 5 and Table 6), likely caused by the  lower quality IMU and the presence of taller buildings along the trajectory. However, after using the tie points in the adjustment, the quality of the Trajectory-II is similar to that of Trajectory-I. RMSE values quickly rise when fewer tie points are used (experiments 2 and 3) or when less accurate tie points are used (experiment 4). This shows that the trajectory accuracy stronger depends on the availability of tie points than for Trajectory-I where the superior IMU is able to bridge longer stretches without tie points. In experiment 5 it is noticeable that the car speed of Trajectory-II is lower than Trajectory-I. Ignoring the distance travelled by the laser scanning car and only considering the scanning time confirms that the IMU system used for Trajectory-II is comparatively less accurate than the IMU used for Trajectory-I. Even compared to experiment 4b, in which 20% tie points of a medium quality were used, the RMSE values for the check points in the 90 s trajectory reconstructed with IMU data only increase by 21, 7 and 16 cm for the X-, Y-and Z-coordinates respectively. Similar to our experience with Trajectory-I, when comparing the result of experiment 6 conducted without soft constraints with experiment 1, it is also clear that soft constraints further improved 5 cm in the RMSE in both X and Z coordinates.

Conclusions
We developed and described an automatic method for the enhancement of 6DOF MLS platform trajectories. For this purpose, we described multiple observations related to automatically acquired 3D tie points matched to corresponding points extracted from aerial photographs, IMU measurements and soft constraints which were used collectively for the trajectory adjustment. We proposed to model the trajectory parameters with B-Splines and also estimated the optimal values of the B-spline knot interval and curve order. To evaluate the developed method we designed several experiments with each experiment focussing on a specific aspect of the correction.
To evaluate the accuracy of the enhanced trajectory the developed method was tested on two independent datasets. The adjustment with all observations brought a large improvement when results were compared with the original Kalman filter solution. The Kalman filter solution of Trajectory-I had RMSE values of 26 cm, 36 cm and 47 cm for the X-, Yand Z-coordinates respectively. With the adjustment, using all introduced observations, these RMSE values were reduced to 9 cm, 14 cm and 17 cm. When only high-quality tie points (55% with standard deviation ≤ 5 cm) were used, RMSE values were further reduced to 9 cm, 11 cm and 16 cm.
For Trajectory-II, the residuals in the Kalman filter trajectory were 1.10 m, 1.51 m, 1.81 m for the X-, Y-and Z-coordinates respectively. After enhancement of the Trajectory-II, using all observations, the remaining errors in the trajectory amount to 12 cm, 15 cm, and 20 cm. After the adjustment with all observations Trajectory-II shows a larger improvement than Trajectory-I because Trajectory-II was more erroneous due to the lower quality IMU. Similar to the first dataset the use of only high-quality tie points (60%) leads to a further improvement in trajectory with RMSE values at 12 cm, 14 cm and 18 cm. The evaluation of the developed method on both data sets confirms that our method can significantly improve the trajectory for applications demanding high accuracy in urban canyons. It also shows that the significantly lower quality of the IMU in the second dataset has little to no effect on the quality of the reconstructed trajectory if sufficient tie points to aerial photographs are available.
It is preferred to use tie points with weights or only high-quality tie points for Trajectory-I and -II respectively because both tie point sets provide comparable results.. Moreover, the results obtained by 80% or 90% inconsistently located tie points lead to RMSE values of 12 cm, 19 cm and 30 cm for Trajectory-I and 15 cm, 21 cm and 40 cm for Trajectory-II. However, when removed consistently, the same amount of tie points can achieve RMSE values of 10 cm, 14 cm and 17 for Trajectory-I and 12 cm, 14 cm, and 20 cm for Trajectory-II dataset. Therefore, we conclude that the consistent availability of the tie point is crucial for the developed method. The consistent availability of the tie points is dependent on the existence of the road-marks in the dataset. The earlier developed image matching procedure (Hussnain et al. (2019)) is already designed to be robust enough to handle the difficult matching situations, since we cannot afford long stretches of trajectory without tie points.
The quality of the tie points depends on the accuracy of the exterior and interior orientation parameters of the aerial images as well as the accuracy of matching in the aerial images. The tie point quality is assessed by the standard deviation of the multi-view triangulation. Moreover, the high aerial image resolution provides precise information required to reliably detect and estimate the subpixel-level features which can improve the estimation of high-quality tie points. In this way the use of aerial images can strongly improve the MLS trajectory estimation and thereby reduce the need for surveys on the ground.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.