A MARKER-FREE CALIBRATION METHOD FOR MOBILE LASER SCANNING POINT CLOUDS CORRECTION

Mobile laser scanning systems (MLS) have been widely used in collecting three-dimensional point clouds for many applications, such as 3D mapping, road facilities inventory and high definition map. Although MLS is calibrated accurately to obtain precise locations of point clouds, it is still challenging to obtain precise locations of point clouds in the areas of GPS signal denied or narrow streets with high dense buildings, resulting in uneven position deviations of point clouds between the overlapping trajectory areas. In this paper, a marker-free calibration method is proposed to solve the above problems. The proposed method firstly partitions the trajectory into segments according to the error distribution while collecting the point clouds. Secondly, the features in each overlapped area are extracted and a kind of Locally Aggregated Descriptors are built for the matching. Thirdly, a coarse-to-fine pairwise point clouds alignment is applied on the overlapping areas. Finally, the global alignment of point clouds is fulfilled with minimizing the position deviations between the overlapping areas and the adjacent segments. The proposed method has been used to correct the point clouds from several different MLSs. Experiments show that this method automatically locates and corrects the uneven position deviations in terms of good robustness and efficiencies. Besides, it proves that the proposed method is also an easy-to-use tool for the automatic correction of MLS point clouds position and boresights.


INTRODUCTION
Mobile laser scanning system can efficiently and densely obtain three-dimensional geometry and texture information of the road environment, which provides a new technical means for highresolution earth observation. Mobile laser point cloud plays a very important role in national major engineering applications such as 3D mapping, road facilities inventory, high definition map (Kukko, 2013). However, due to the comprehensive influence of GNSS positioning error, IMU attitude determination error, scanner angle and distance measurement error, multisensor installation error and so on, the position accuracy of MLS point clouds is difficult to meet the centimeter-level accuracy requirements, showing in Figure 1. Especially in high-rise urban areas, satellite signal occlusion is serious and positioning accuracy is poor, resulting in a decimeter or even meter-level deviation between mobile laser scanning system. The position inconsistency phenomenon occurs among the point clouds collected from the same area, which is difficult to meet the accuracy requirements of the above applications without making a registration pass. (Xu et al, 2015;Yang et al, 2017). More severely, it is quite difficult to automatically locate and correct the uneven position deviations from a long-distance trajectory (e.g., over 20 km). Directly align the MLS point clouds using the registration method such as Iterative Closest Point (ICP) or normal distribution transform (NDT) will not solve the inconsistency and even causing a bad matching result due to the non-rigid deformation inside the point clouds, the mimic geometry and context and the limited overlapped areas. In this paper, a series of approaches including point cloud partition, feature extraction and description, hierarchical pairwise registration and global optimization are designed to address this  Corresponding author issue and effectively produce a more consistency point clouds in the same areas. Overall, the main contributions of this paper are three-fold: 1) An adaptive partition method for MLS point clouds considering the error distribution characteristics of point clouds is proposed, which ensures the overlapping areas among the overlapping segments and improves the reasonableness of partition for MLS point clouds.
2) A pairwise registration strategy of MLS point clouds from coarse to fine is designed, which effectively solves the large deviation and non-rigid deformation of MLS point clouds.
3) The global alignment optimization of MLS point clouds is applied through optimizing the constructed objective function, which combines the distance of correspondence between overlapping sub-regions with the transformation smoothness between adjacent sub-regions.
Before giving a detailed description of the proposed approach for solving the uneven position deviations from a long-distance trajectory in Section 3, the outline of MLS position inconsistency correction and point cloud alignment is elaborated in Section 2. In Section 4 experimental results for correcting point clouds collected from different MLSs are described. The conclusion and feature work are drawn at the end of this paper.

RELATED WORKS
Alignment of point clouds is a long-standing problem in the field of computer vision and photogrammetry. The steps to deal with the inconsistency of point clouds are usually point clouds partition, local pairwise registration and global optimization.
The MLS point clouds partition strategy is summarized into two categories: time-based partition (Han et al., 2014;Yan et at., 2018) and space-based partition (Takai et al., 2013). The former usually builds a time error model to refine the accumulation error from IMU drift and GPS signals loss. One assumption is that there is a consecutive changing relationship between pose error and time stamp, thus the error could be compensated by a pairwise low order Polynomial function. Once the errors in trajectory are corrected, the MLS point clouds could be recalculated with the pose parameter and sensor observation equation. Han et al. (2014), proposed a time-variant model of vehicle trajectory and calculated the coordinate of the point clouds using GCPs. Yan et al. (2018) sliced the MLS point clouds by a fixed period of time to solve the non-rigid registration of multi-strip MLS point clouds. The registration of MLS point clouds could usually obtain high precision but once the vehicle runs slowly or even stops at a traffic light, the MLS point clouds in this area will be very dense, and time-based partition will not work at all. The latter method directly corrects MLS point clouds by slicing the MLS point clouds along the trajectory into segments, assumes that the error in each segment is very small considering the IMU's stability. The length of the segments is usually altered. Treat segment as process unit and then calculate the transformation relationship of overlapping areas' point clouds through geometric and semantic features extracted from point clouds. Takai et al. (2013) partitioned the point cloud into segments according to the key points of the trajectory. After the registration the difference among point clouds are reduced to 50mm.
The point clouds registration has been drawn considerable attention (Pomerleau et al., 2015;Cheng et al., 2018;Dong et al., 2020). Coarse-to-fine is a common strategy employed. The coarse alignment approximately finds an initial transformation contains rotation and translation matrix in a relatively fast and accurate way, thus the fine registration approach such as the NDT and its variants (Magnusson et al,.2007;Das * and Waslander †, 2012), and the ICP algorithm (Besl and McKay, 1992) and its variants (Tazir et al., 2018), are implemented to achieve the final transformation. Most of the fine registration algorithm have a bottleneck that the processing time is too long while the point clouds' number is more than 0.1 million, but the dense point clouds are widely seen in MLS. Although there are also many tricks to accelerate the calculation, when the overlapping area is small and the context is similar, the result is poor and unstable.
Recently the learning-based point cloud alignment methods developed rapidly, such as deep neural network (Choy et al.,2019;Lu et al., 2019), but the point clouds are more sophisticated and these methods have difficulties in handling this.
Based on the point cloud partition strategy and point clouds registration method, Yan et al., (2018) split point clouds into segments uniformly according to scanning time, and obtained accurate homonymous correspondence through two-step ICP algorithm, finally corrected trajectory error by using the time-variant error model. Yu et al., (2015) extracted the object of different scales in MLS point cloud and found correspondences between the objects by using multi-scale ICP, improving the quality of MLS point cloud data from coarse-to-fine. This method can solve the problem of position consistency correction of MLS point cloud in urban scale, but it depends heavily on the quality of an object and feature extraction of various types of objects and is less feasible in complex urban scenarios. Besides, Shiratori et al. (2015) detected loop closures in MLS point clouds, and used loop constraint to boost the overall quality of point clouds in the overlapping regions. Yang et al. (2018) proposed a method of urban point cloud mapping based on pose map. They combined the classifier and optimization method to eliminate the unstable local correction results, and the conditional random field is used to identify and prune the dynamic targets, preventing the negative impact from dynamic objects.
Global optimization has been a necessary step on point cloud position correction (Theiler et al. 2015) Especially in MLS point clouds registration mission, the individual point cloud segment pose could be tackled as nodes in a graph and the constraints between point clouds be formulated as edges. And to find a spatial configuration of the nodes means solving the MLS point cloud registration problem and achieving the global position consistency in overlapped areas. Theiler et al. (2015) formulated the TLS point clouds into a graph structure and conducted a global optimization to evenly spread the residual errors into all TLS point clouds after coarse-to-fine pairwise alignment.
Although the above methods have achieved promising results in specific scenarios, there are still many limitations in robustness and large deviations correction. To address the challenging issues, this paper proposes a novel method to address the position inconsistency of MLS point clouds.

METHODOLOGY
The proposed method consists of four key steps: MLS point clouds partition, feature extraction and description, pairwise registration of overlapping point clouds, and global optimization, as illustrated in Figure 2.

MLS Point Clouds Partition
Inspired by the Takai et al. (2013), This paper slices the MLS point clouds into segments based on the error distribution. Specifically, the position inconsistency deviations in the MLS point clouds often occur to the places where the speed and direction of the moving vehicle change sharply, the position deviation in these places are non-linear distribution. Besides, to improve the adaptability in large deviation, symmetry, and sparse geometric features scenarios, the overlapping areas and crossroads should stay completion for the convenience of the further process. The proposed partition strategy could be summarized as the following four steps: selection of candidate segment points, overlapping points detection, segment points filter and MLS point clouds partition.

Selection of Candidate Segment Points
Based on the trajectory of vehicle, the acceleration and angular velocity of the trajectory will be calculated and the trajectory point whose acceleration or heading angular velocity is larger than the threshold ℎ or ℎ is treated as a velocity change point or angular changing point respectively. Let , , and be the first, last velocity changing points, angular changing points of a continuous pose trajectory.  selected by extending a certain distance ℎ1 outward along the trajectory as the candidate segment point, dotted in red solid points in Figure 3(a), 3(b). Meanwhile, the intersection of trajectories is detected. Crossroad candidate segment point is found by extending a certain distance ℎ2 outward along , shown in Figure 3(c). When the interval distance between adjacent candidate partition points is larger than , the trajectory points at the midpoint are taken as new candidate segment points, dotted in hollow red point in Figure 3(d). Then the candidate segment dotted in red solid points are found and illustrated in Figure 3(e).

Overlapping Points Detection
To correct the inconsistency in the overlapping areas, for each red point as dotted in Figure 3(e), the nearest trajectory point whose accumulated distance along the trajectory is large than but also falls in a radius region are labelled as its overlapping point, as dotted in the hollow red point encircle by a rounded rectangle in Figure 3(f). And the candidate segment point who has the overlapping point will not be candidate segment point any more. The purpose of searching for overlapping segment points is to maximize the overlapping part between the overlapping segments and make subsequent registration more robust.

Segment Points Filter
The distribution of candidate segment points may be uneven, and the segment points near the intersection may be too dense. Hence, the final segment points need to be selected again from the candidate segment points. Firstly, the candidate segment point who has the overlapping point should be deleted; Secondly, candidate segment points in the middle of the intersection points should be deleted and then the trajectory is evenly segmented with a certain interval , as shown in Figure 3(g). And the non-rigid deviation in each segment might be too small to ignore, the rigid transformation between small segments should be handled.

MLS Point Clouds Partition.
According to above three steps, the MLS point cloud could be segmented by the time stamp between trajectory and point clouds. With the small segments obtained, the small segments are then merged into bigger large sub-segments according to the accumulated length , as shown in Figure 3(h).

Feature Extraction And Description
The small segments and the large segments are used as process units. A feature-based approach is applied to correct the position deviation of the overlapping areas in MLS point clouds. Binary Shape Context (BSC, ) is adopted to extract the feature points descriptors in each small segment for matching the correspondence without artificial markers. The covariance matrix is constructed based on a point and its neighboring points in a certain radius, and the eigen vector and eigen value is calculated. Once a point's curvature large enough, it is selected as a key point and then BSC feature descriptor is extracted around it. The BSC features in small segments are concatenated to form the large segments' BSC features. After local feature descriptor is extracted, a common method Vector of Locally Aggregated Descriptor (VLAD, David et al., 2006) was adopted to leverage the function of the feature descriptors. Specifically, BSC features are clustered by K-means to build a visual word dictionary for large segments and small segments.

Coarse-To-Fine Pairwise Registration
For MLS point cloud data acquired in complex urban scenarios, the general pairwise registration algorithm is of poor robustness and is difficult to apply into complex scenarios with large deviations, symmetrical structures, and sparse geometric features. When the segment length is long, the non-rigid deformation may be large, but the geometric features are rich with a high precision in pairwise registration under the above mentioned hard scenarios. When the segments' length is short, the deformation in the segment may be small, but the geometric features are relatively sparse, and the accuracy of pairwise registration is high in local region, but the robustness is relatively poor. Thus, a hierarchical coarse-to-fine method combining long segment and short segment is proposed to complement each other in pairwise registration. The approach includes three parts: visual words accelerated large segments registration; initial pose information assisted small segments registration and mismatch removal.

Large Segments Registration
When 1 between two large segments is larger than , the two large segments are considered as revisit. For the subsegments contained in the large overlapping segments, the 2 between the possible overlapping sub-segments is calculated. When 2 is larger than , the two small segments are considered as revisit. For a large segment , a set = { |IOU > IOU , j < i} , containing its' overlapped segments , is built. IOU is the overlap value between the and . And the segments in whose unique id is consecutive are grouped as sub-segments. That is, for instance, is an element in , and its adjacent large segments consist of = { | ∈ , |k − m| < ε} , shown in Figure 4. BSC features of all sub-segments contained in are grouped to get a feature set Ω , while BSC features of sub-segments in forming a feature set . Then the visual word dictionary is used to search the correspondence of , for further eliminating completely dissimilar features and improving the efficiency of feature matching. The root mean square error (RMSE) then are spread to the inside small segments. Figure 4. Pairwise registration of large segments. Seg i is the target point cloud segment, and Seg k−1 , Seg k , Seg k+1 are its overlapping segments.

Small Segments Registration
According to the pairwise registration result of large segments, the search range could be narrowed during the small segments registration. Searching radius boosts the BSC feature matching considering not only the similarity between features, but also the location constraints of feature points. However, because of the feature representation capabilities, there are mismatch segments in the registration, thus, a mismatch removal is essential.

Mismatch Removal
The mismatch removal includes: the number of the feature matching, self-rotation of each segment and the relative transformation relationship between consecutive adjacent segments. Apparently, the more the number of the homonymous features among the segments, the more reliable the registration results are. The feature matching number which is smaller than a certain threshold is pruned. The POS system carried by the mobile laser scanning system generally has high accuracy, and the attitude deviation of the acquired point clouds in the overlapping area should be small. Therefore, when the relative rotation (roll, pitch, yaw) is greater than ℎ , the pairwise registration result is considered unreliable.
T-test is utilized to eliminate pairwise registration results of abnormal relative transformation in the consecutive adjacent segments. Formulating the moving distance of bounding box center after registration as statistical variables. For segments Seg i and Seg j , the bounding box center's local coordinates are O i and O j , corresponding pose in the global coordinate system separately are T i and T j , and the relative pose obtained by pairwise registration is . The distance that Seg j moves after registration could be calculated by the equation as: Mean value, D ̅ and the standard deviation, σ, could be calculated from all overlapping segments' moving distance. Once a sampled , its corresponding t value calculated using t = Dist−D σ n−1 is greater than the critical value, its corresponding feature matching is discarded.

Global Optimization
After correcting the local position of the MLS point cloud from coarse-to-fine pairwise registration, the relative transformation relationship between the segments of the overlapping point cloud is obtained, but the transformation relationship between the segments is complex and inconsistent. Therefore, these transformation relationships need to be adjusted globally so that the position deviations are evenly distributed in the point clouds.
The small segments and correspondence between segments are

· · · · · ·
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) described as a graph structure, and are optimized by solving a non-linear least square equation through Levenberg-Marquardt method (Hartley et al., 2003). When there are many overlapping areas, the optimization will loop more iterations and consume too much time. Thus, a simple graph is firstly constructed based on the relative transformation in adjacent segments to get a finer pose, then construct the graph with additional constraint to get the optimal pose transformation.

First Stage Optimization
The overlapping small segments are constrained with small variance before and after optimization, while maintaining the transformation between the adjacent segments. Based on this idea, an object function is formulated which could be divided into data term and smooth term, as equation 1: where S and |S| = the segments and the total number of small segments Φ(S i , S j ) = the relationship between segment S i and S j = the transformation matrix corresponding to S i , −1 = the transformation matrix between segment S i and S j T S i S j ′ = the transformation matrix between segment S i and S j after global optimization d(T S i S j −1 T S i S j ′ ) = the variance of relative transformation matrix between S i and S j before and after global optimization, N(S i ) = set of adjacent segments of S i .

Φ(S i , S j )
is an indicator, that is, equal to 1 when the segment is an overlapping segment, otherwise 0. Data term forces the variance during iteration will not be too large but just a slight adjustment of the pair-wise registration. And the smooth term keeps the transform consistency of the adjacent segments. In case of multi overlapping segments, the pairs that have the most numbers of feature correspondence are remained, ensuring the high precision of correspondence. Here, the weights are set to be equal between data term and smooth term.

Second Stage Optimization
Rely on the fast convergence of the first stage optimization, each overlapping segment has a better initial pose. Then, a more powerful constraint is designed to cover all the segments with adding an inertial term, representing as equation 2. The new object function consisting of data term E data , smooth term E smooth , and inertial term E inertial .
where C = the set of the homonymous correspondence (F i , F i ′ ) = a pair of homonymous features T F i and T F i ′ = transformation matrices associated with F i and F i ′ T F i (F i ) and T F i ′ (F i ′ ) = means transform the F i and F i ′ to target coordinate system d(T F i (F i ), T F i ′ (F i ′ ) = distance residuals φ(S i, S i+1 ) and |φ(S i, S i+1 )| = homonymous points and size between adjacent segment and +1 and ′ = jth pair of homonymous points = pose of segment ( ( ), +1 ( ′ ) = Euclidean distance residuals of homonymous points in adjacent segments and +1 ( , ′ ) = variance of transformation of segment Note that, ω smooth and ω inertial are the weights of the smooth term and inertial term, ω smooth = E data / E smooth , ω inertial is very small and is generally set to 10 −3 . The data term is the sum of Euclidean distance residuals of the homonymous points inside the overlapping segment. The data term minimizes the gap in geometry space between point clouds scanned in the mostly same area. Due to the mismatch in hard scenarios, the smooth term is used to punish the abnormal transformation at a local scope of adjacent segments, controlling the distance residuals of the homonymous point and the third inertial term forces the variance of each segment limited in a small scope during the optimization process. Both smooth and inertial term compact the segments and relief the gap between segments. The graph is constructed incrementally along the vehicle trajectory. If there are overlapping between the current and previous segments, both data term and smooth term are updated. If there are no overlapping detected, only the smooth term will be updated.

Dataset
The MLS point clouds partition performance is evaluated on point clouds acquired from a RIEGL scanner in Shanghai including 1 billion points. The total trajectory is about 20km, shown in Figure 6. And small point clouds, as listed in Table 1 and shown in Figure 5, are used to demonstrate the result of the position inconsistency correction.

Point Clouds Segments
We set the segment length separately from 5m to 40m and the result is displayed in Figure 6. It can be seen that no matter how long the sub-segment is, the MLS point clouds are basically in the same position, there is adequate overlap between the overlapping segments. The MLS point clouds are completely preserved at intersections and bends. Other parameter we set is similar to Takai et al. (2013), but adjusted according to our real data. And we find there are slight effect on choosing parameters except the , which dominates the length of the segments, the , which enlarges the possibility of detecting overlapped and the which keeps the even distribution of segment points. Figure 6. Result of point cloud partition, segment length is (a) 5m, (b) 10m, (c) 20m, (d) 40m. The left column is top view of segments; the upper center shows the segments at round-trip; the down center shows the segments at crossroads, and the right column shows the segments at bend.

Position Deviation Correction
Given no ground truth provided. the manually labelled homonymous points in the overlapping point cloud are used as the ground truth to evaluate the relative position accuracy. RMSE of homonymous points between overlapping segments is used to evaluate the relative position accuracy of overlapping segments, and RMSE of virtual homonymous points between adjacent segments is used to evaluate the relative position accuracy. Shown in Table 2.
From considerable experiments, we find that 20m is almost the optimal length for the MLS point clouds in Shanghai. When the segment length is small, such as 5m, or even smaller, mismatch results in registration stage are too many to obtain a good optimization result, although the smoothing terms of adjacent segments are the largest. These mismatch makes the optimization convergence very slow, even failure, as shown in Table 2. When the segment length is 40m, or even larger, the non-rigid deformation inside the segment is relatively large, and there is still a relatively large position inconsistency after position consistency correction. But larger segments we merged from small segments stabilize the registration result, preventing many gross errors and the final result is also acceptable. When the segment length is between 10 and 20m, the registration results of pairwise registration are still reliable, with negligible non-rigid deformation inside the segments. Meanwhile, the global optimization succeeds as expecting. Figure 7 shows the location changings during the process, and the adjacent segments' consistency is qualitatively demonstrated. The position inconsistency in the same areas after processing by our algorithm could be reduced to cm level depending on the origin point cloud, for example, the beam parameter calibration value. If enough GCPs are available, the MLS point clouds could be easily transformed into absolute coordinate. Figure 8 shows more position correction result.

CONCLUSION
A marker-free MLS point clouds position inconsistency correction method is proposed in this paper. The point clouds collected from different sensors were used to verify the effectiveness and robustness of the proposed method. The adaptive hierarchical partition considering the error distribution characteristics of point clouds guarantees the overlapping segments and improves the reasonableness of MLS point clouds partition. Visual word and prior pose information acceleration matching strategy ensures the efficiencies of the matching, and a coarse-to-fine hierarchical registration method guarantees the robustness of point clouds alignment, leading to effective correction of large deviation in the overlapping MLS point clouds.
In the future, learning-based methods will be evaluated on the large urban scene point cloud alignment and pruning the moving object to relieve the mismatch.