CALIBRATION OF A MULTIBEAM LASER SYSTEM BY USING A TLS-GENERATED REFERENCE

Rotating multi-beam LIDARs mounted on moving platforms have become very successful for many applications such as autonomous navigation, obstacle avoidance or mobile mapping. To obtain accurate point coordinates, a precise calibration of such a LIDAR system is required. For the determination of the corresponding parameters we propose a calibration scheme which exploits the information of 3D reference point clouds captured by a terrestrial laser scanning (TLS) device. It is assumed that the accuracy of this point clouds is considerably higher than that from the multi-beam LIDAR and that the data represent faces of man-made objects at different distances. After extracting planes in the reference data sets, the point-plane-incidences of the measured points and the reference planes are used to formulate the implicit constraints. We inspect the Velodyne HDL-64E S2 system as the best-known representative for this kind of sensor system. The usability and feasibility of the calibration procedure is demonstrated with real data sets representing building faces (walls, roof planes and ground). Beside the improvement of the point accuracy by considering the calibration results, we test the significance of the parameters related to the sensor model and consider the uncertainty of measurements w.r.t. the measured distances. The Velodyne returns two kinds of measurements—distances and encoder angles. To account for this, we perform a variance component estimation to obtain realistic standard deviations for the observations.


Motivation
The time-of-flight multi-beam LIDAR Velodyne HDL-64E S2 (Velodyne) is especially known from the DARPA Urban Challenge, where it was used for autonomous navigation.Another application is mobile mapping.In this field precision and accuracy of the measured point coordinates are important.But a first visual inspection of a Velodyne scan showed, that it is quite imprecise (see Figure 1).According to the manual (Velodyne LiDAR, Inc., 2011) the distance accuracy is below 2 cm (one sigma).In (Glennie and Lichti, 2010) a RMSE of 3.6 cm is reported for the planar misclosure w.r.t. the manufactures calibration values.After they applied their presented calibration procedure, the RMSE reduced to 1.3 cm.Existing approaches lack the usage of precise reference data for a broad range of distances.Therefore the existing calibration approach presented in (Glennie and Lichti, 2010) was extended to incorporate precise reference data.

Related Work
Several approaches exist to calibrate multi-beam LIDAR systems.Conceptually they can be divided into approaches which require data of additional sensors and calibrations based only on data constraints (especially planarity).
Approaches exploiting additional sensors An unsupervised technique is presented in (Levinson and Thrun, 2010) for moving sensors.The approach relies on the weak assumption that the 3D points lie on contiguous surfaces and that points from different beams pass the same surface area during the motion which is captured by an inertial measurement unit.The solution is obtained by minimizing an energy functional for aggregated points acquired across a series of system poses.Figure 1: Top view of a point cloud representing a wall area.Mirzaei et al. (2012) showed how to calibrate a multi-beam system with a rigidly connected camera (intrinsic and extrinsic) by observing a planar calibration board with fiducial markers.The camera supplies the position and orientation of the calibration board.Additionally, an approach is presented to obtain a closedform solution for the determination of approximate parameters.These are used as an initialization for a non-linear optimization, which also incorporates the planarity.
Approaches using data constraints Glennie and Lichti (2010) recorded a "courtyard between four identical buildings" from multiple positions with different orientations.Planes extracted in the 3D point clouds were automatically detected and used for calibration.To get a unique solution, the angular offset for one of the lasers had been fixed among others.Also the effect of the incidence and vertical angles on the residuals was analyzed, which showed raising residuals at an incidence angle of 60 • to 65 • and above and a vertical angle of −20 • and lower.
An analysis of the temporal stability of calibration values for the Velodyne system is given in (Glennie and Lichti, 2011).
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey Muhammad and Lacroix (2010) presented a simpler approach, which requires only one calibration plane measured at different distances between 4 m and 14 m.They calculated the plane parameters and minimized the point to plane distance iteratively, but also updated the plane parameters in every iteration step.

Velodyne setup
The 64 beams of the Velodyne are placed on top of each other in a spinning head.Additionally for each laser a detector is located in the head.The head spins with a rotation speed between 5 Hz and 20 Hz.
The lasers are aligned at different vertical angles and heights.They cover a vertical field of view from 2 • to −24.33 • .In the head they are split into four groups (see Figure 2

Contribution
Multiple scenes captured by a reference data scanner are used.
The reference gives the ability for exact error estimation.The scenes contain a wide range of distances, so that a large part of the possible range of the LIDAR is mapped.Also a significance test is performed to show the importance of the scale factor.Furthermore we investigate the distance depending uncertainty.A variance component estimation is performed, which leads to reasonable standard deviations for the different observation types, i.e. distances and encoder angles.

Approach
For the calibration the terrestrial laser scanner Z+F Imager 5010 supplies the reference data, which provides measurements one magnitude more accurate1 than that of the Velodyne2 .To get a mathematical manageable form, planes are extracted from the reference scans.Therefore the scenes should contain enough planes, which is the case for man-made objects like buildings.
In addition the scenes used for calibration should provide a broad range of distances.The perfect case would be to cover the whole Velodyne range from 0.9 m to 120 m.The Velodyne scanner has only a vertical field of view of 26.33 • .Therefore the scene should be recorded from more than one viewpoint and the Velodyne should be tilted at each viewpoint to cover more parts of the scene (geometric configurations) and to get additional constraints.As a starting point the calibration values supplied by the manufacturer are taken.During the calibration process the registration parameters are adjusted, too.
As pointed out by Glennie and Lichti ( 2010), there exists a strong correlation between the horizontal and vertical angles and respective offsets.Therefore the horizontal and vertical offsets are excluded from the calibration procedure and the manufacture's values are considered to be fix.
The proposed calibration process consists of the following steps: 1. Generation of reference data: (a) Record one or more scenes with the reference scanner.
(b) Extract planes from the reference point cloud.

Measurements:
(a) Record every scene from multiple viewpoints and tilting angles with the Velodyne, to get a good coverage of the planar parts of the scene.
(b) Register every Velodyne scan to the related reference cloud.

Data association and calibration:
(a) Extract all Velodyne points including raw data, which lie next to a reference plane.
(b) Use the Gauss-Helmert-Model to estimate the calibration parameters.

Generation of reference data
For most of the steps described in this section, our implementation uses the Point Cloud Library (PCL) (Rusu and Cousins, 2011) -especially for point normal estimation, region growing and registration.
The first step is to reduce the reference datasets.Therefore a subsampling is applied, which is driven by a simple one-per-gridded cell strategy, i.e. selecting the point closest to the cell center.Each grid cell has a side length of 1 cm.The reference data is very dense, especially near the scanner positions.Such dense data is not needed and slows down the data processing and visualization.Also it does not provide extra precision and accuracy.
Then the planes are extracted by a region growing method.As decision criterion for adding points the angle between the normal vector of the seed point and the normal vector of the picked point is calculated.If it is below 2 • , it will be added to the region.Afterwards a RANSAC-based plane fitting is carried out for each region.For the derivation of the required threshold we assume Gaussian-distributed point coordinates with a standard deviation of σ = 1 cm and a probability of 95% for selecting an inlier point.The threshold for the point-to-plane distance is then The scans are manually preregistered.After that a registration method based on the Normal-Distributions Transform is used, see (Magnusson, 2009).A point from a scan is included for the calibration if the point-to-plane distance is below 10 cm and the distance of the perpendicular foot of this point to the next reference point is shorter than 5 cm.The reference point density is high enough to permit only a low distance along the surface.Along the plane normal a higher threshold is needed to prevent an exclusion of valid points from the calibration.The Velodyne raw data and the corresponding plane are saved for the calibration process.

Constraints and sensor model
For the calibration we consider the incidence of observed laser points Xij captured by the laser i and corresponding planes {np, dp}, represented by their normals np and the shortest distances dp to the origin of a world coordinate system.The implicit constraints are then for the adjusted parameters and observations in world coordinate system Sw.The transformation from the sensor to the world coordinate system is for all points with the scanner's position w t and the rotation matrix w Rs.According to the Velodyne setup (Section 1.3), the 3D point is given by with the projection center Cij, the direction Dij of laser i and the measured distance sij, corrected by an offset bi and a scale factor ai.With the measured encoder angles ij , the angular horizontal corrections βi, and the horizontal and vertical offset Hi and Vi the center for each laser is and with the vertical correction angle δi the direction for each laser is Gauge constraint As pointed out by Glennie and Lichti (2010), not all 64 horizontal angular offsets can be estimated simultaneously since only 63 parameters βi are independent.To fix this gauge freedom, we could fix the parameter of one laser, say β1 := 0. Alternatively, we can introduce the centroid constraint which is linear in the unknown parameters.This constraint can be incorporated easily into the calibration process, which will be explained in the next section.
Functional model We use two different types of constraints for the observations l and the parameters x which have to be fulfilled for their estimates: conditions g( l, x) = 0 for the point-planeincidences and restrictions h( x) = 0 to fix the gauge freedom.
Stochastic model For the application of the Gauss-Helmert-Model an initial covariance matrix Σ (0) ll of the observations l is assumed to be known which is related to the true covariance matrix Σ ll with an usually unknown scale factor σ 2 0 , called variance factor.
The Velodyne returns two types of observations for each laser beam: the measured distance and the encoder angle.For the variances of these observations, initially only a rough guess exists.Furthermore, the ratio of the variances for both observation groups is unknown.Therefore, we apply a variance components estimation by dividing the observations into distance measurements l1 and encoder angles l2 and introducing corresponding covariance matrices for both groups assuming statistical independence, cf.(McGlone et al., 2004).The variance components σ 2 01 and σ 2 02 are estimated within the iterative estimation procedure, see Section 2.5.
Estimation The optimal estimates in terms of a least-squaresadjustment can be found by minimizing the Lagrangian with the Lagrangian vectors λ and µ.With the Jacobians and we obtain the linear constraints by Taylor series expansion For the residuals v and approximate observations l0 the relation Setting the partial derivatives to zero yields the neccessary conditions for a minimum: ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey 16) into ( 17)) yields where Σww = BΣ ll B T is the covariance matrix of the contradictions (15) derived by error propagation.Thus we obtain the normal equation system . ( 21) Since N has full rank, the solution is given by The covariance matrix of the estimated parameters Σ xx consists of the upper square of N −1 corresponding to x and the update for the estimated parameters is The adjusted observations are with the estimated corrections due to ( 16).The iteration is started with given initial values x0.

Variance components
For the estimation of the variance components the covariance matrix of the residuals ( 25) is needed Splitting the common variance factor into two parts, the estimates are within each iteration step, with the redundancy numbers for both observation groups in the denominator and improved covariance matrices Σii for both observation groups.

RESULTS AND DISCUSSION
After the explication of the calibration setup and the generation of reference data, we summarize and interpret the estimation results.These give reasons for a variation of the uncertainty of the measurements depending on the distances which is investigated in more detail.

Setup and generation of reference data
Three scenes of a building scanned by a Z+F Imager 5010 serve as reference data.The first scene was taken from two positions.
In contrast to the first two scenes the third is indoor.
The first scene was recorded by the Velodyne from a stairway to get different heights (three floors).Additionally, the Velodyne Sc.Pos.Rot.w/o tilted base Rot. with tilted base ht 1 1 0 Table 1: Position numbers, sensor orientation and height above ground for the collected data sets.
was turned four times at each position.The same was done with the Velodyne fixed on a tilted base of 25 • , which results overall in 24 scans of the first scene (see Table 1).Because of insufficient contradictions for registration two scans could not be incorporated in the calibration process.
The scan-position for the second scan was on the ground with azimuths of 0 • , 45 • and 90 • .The last scan inside a room was made with the tilted base on a table at four azimuth positions.ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume II-5/W2, 2013 ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey procedure described in Section 2.2, which in addition to the extracted planes are the base for the parameter estimation.The resulting distances are shown in Figure 4.

Results
For the stochastic model we assume a standard deviation of 0.02 m for the distance observations and of 0.09 • for the encoder angles, which is the standard deviation for the distance accuracy and the resolution of the encoder, both provided by the manufacturer.Table 2 summarizes selected results for the parameter estimation.For the set of 64 lasers, the occurring parameter range and the minimal and maximal estimated standard deviation for the estimated parameters are given.
In (Glennie and Lichti, 2010) the scale factor a is introduced for each laser (see ( 3)).According to  2: Shows the estimation results comprising the results for the 64 lasers.For all estimated parameters, the range, the minimal and maximal estimated standard deviations are listed.
Table 3 depicts the initially assumed and the estimated standard deviations for the measurements.The estimates reveal that the assumed standard deviation had been somewhat too optimisticprovided that the mathematical model holds.
observation type a priori a posteriori (estimated) distance 0.02 m 0.0366 m encoder angle 0.09 • 0.1426   that the points are calculated with the manufacturer calibration data and procedure.In contrast to this for estimation the new calibration data and (3) are applied.The minimal and maximal manufacturer misclosure values result from the allowed point-toplane distance of 10 cm.The standard deviation of the misclosures of the estimation is smaller compared to the manufacturer one, which can also be seen in Fig. 5.In opposite to this the range of the misclosure increased.This is caused by incorrect point-to-plane mapping, but very few points are affected.Fig. 6 and Fig. 7 show the histograms for the residuals, separated for distances and encoder angles.Obviously, the residuals feature a Gaussian scale mixture distribution.This provokes the assumption, that the standard deviations differ depending on the measured distance or the incidence angles.To check this, we grouped the distance residuals and planar misclosures into distance classes with a binning of one meter, see Fig. 8.The plot reveals that the observations are subject to large fluctuations for different distances.This gives evidence for the assumption.It is also illustrated, that the misclosures standard deviations increase according to the distance.

CONCLUSIONS AND FUTURE WORK
We presented a reference-based calibration procedure which incorporates the advantage of a broad range of distances.After the proposed calibration process the standard deviation of the planar misclosure is nearly halved from 3.2 cm to 1.7 cm.The variance component estimation as well as the standard deviation of the range residuals reveal that the manufactures standard deviation of the distance accuracy with 2 cm is a bit too optimistic.
The histograms of the planar misclosures and the residuals reveal that this quantities are not normal distributed.Our investigation of the distance depending misclosure variance change is one reason.Other sources were investigated by Glennie and Lichti (2010): the incidence angle and the vertical angle.A further possibility is the focal distance, which is different for each laser and the average is at 8 m for the lower block and at 15 m for the upper block (Velodyne LiDAR, Inc., 2011).This may introduce a distance depending-but nonlinear-variance change.Further research is needed to find the sources of these observations.Another point is the question, if the mathematical model holds, or if additions are needed.
(a)).This results in a horizontal offset from the center of the head.In addition to this a horizontal rotation angle exists for each laser.The resulting beam pattern is shown in Figure2(b): The numbers from 0 to 31 belong to the upper groups and from 32 to 63 to the lower groups (see Figure2(a)).Even numbers correspond to lasers on the left side of Figure2(a) and odd to the right side.Laser beam pattern projected on a plane in a distance of 50 m (perspective projection).

Figure 3 :Figure 4 :
Figure 3: Extracted regions of the first scene.After the scanning, a region growing for planar patches had been performed for the reference data sets, e.g.see Fig.3.Then the Velodyne data sets had been registered to the corresponding reference scan.After that, the points are extracted according to the

Figure 5 :
Figure 5: Point-to-plane misclosure with points calculated with manufacturer calibration procedure and values and calculated with (3) and estimated values plotted in the range of −0.11 m to 0.11 m.

Figure 6 :
Figure 6: Histogram of the distance residuals plotted in the range of −0.15 m to 0.15 m.

Figure 7 :Figure 8 :
Figure 7: of the angle residuals plotted in the range of −0.3 • to 0.3 • .
Table 2 the values of this parameter are close to one.Therefore we performed a test on significance by checking if one or all of these scale parameters are equal to one based on the estimates for the 64 scale factors a. −24.849108 • 1.670113 • 0.002068 • 0.003310 • β −9.106730 • 9.349247 • 0.002841 • 0.005297 •

Table 3 :
Initially assumed and estimated standard deviations for the observations (var.comp.estimation).The histogram of the point to the associated plane distances (misclosures) is shown in Fig.5.Table4shows the statistical properties of the misclosures.In both cases by manufacturer is meant,