Article Direct Georeferencing of Stationary LiDAR

Unlike mobile survey systems, stationary survey systems are given very little direct georeferencing attention. Direct Georeferencing is currently being used in several mobile applications, especially in terrestrial and airborne LiDAR systems. Georeferencing of stationary terrestrial LiDAR scanning data, however, is currently performed indirectly through using control points in the scanning site. The indirect georeferencing procedure is often troublesome; the availability of control stations within the scanning range is not always possible. Also, field procedure can be laborious and involve extra equipment and target setups. In addition, the conventional method allows for possible human error due to target information bookkeeping. Additionally, the accuracy of this procedure varies according to the quality of the control used. By adding a dual GPS antenna apparatus to the scanner setup, thereby supplanting the use of multiple ground control points scattered throughout the scanning site, we mitigate not only the problems associated with indirect georeferencing but also induce a more efficient set up procedure while maintaining sufficient precision. In this paper, we describe a new method for determining the 3D absolute orientation of LiDAR point cloud using GPS measurements from two antennae firmly mounted on the optical head of a stationary LiDAR system. In this paper, the general case is derived where the orientation angles are not small; this case completes the theory of stationary LiDAR direct georeferencing. Simulation and real world field experimentation of the prototype implementation suggest a precision of about 0.05 degrees (~1 milli-radian) for the three orientation angles.


Introduction
In stationary applications, especially in the case of terrestrial LiDAR systems, direct georeferencing is still in its infancy.Direct georeferencing of mobile mapping systems, on the other hand, has gained wide acceptance in wide range of applications.In its simplest form, a direct georeferencing system is composed of a GPS-aided inertial navigation system mounted rigidly with the mapping sensor.The georeferencing system and the mapping system are both treated as rigid body where the equations for rigid body motion apply.Mobile LiDAR systems are being developed for different applications, terrestrial and marine [1].
The typical approach for georeferencing LiDAR point clouds from a stationary system is through the use of reflective targets placed at ground control points with known coordinates.The transformation (geo-referencing) parameters between a Scanner Own Coordinate System (SOCS) and a Global Coordinate System (GLCS) are obtained by using the known coordinates of the targets.All other points in the point cloud can then be transformed from the SOCS coordinate system to the GLCS coordinate system using the transformation parameters, hence geo-registering the point cloud in the ground coordinate system.
In addition to the coordinate transformation method described above, there are other methods for georeferencing LiDAR data such as back-sighting, sensor-driven and data-driven approaches, see [2] for a comparison of these approaches and [3,17,18] for an application of a data-driven approach.However, these methods rely on the identification of control targets in the point cloud and in the project area.Some of the disadvantages associated with these methods are: the need for additional, sometimes cumbersome equipment; extra time required to set up the targets; and precision surveying of target control points.Obtaining control often requires personnel work in the scene of interest.This may not always be possible due to physical obstacles and/or safety concerns, e.g., in glaciers or on top of mountains; see Section 2 that describes our motivation to such development.One can increase the speed and efficiency of data acquisition georeferencing and circumvent limitations due to the nature of the scanned area by using a direct georeferencing method that doesn't rely on scanner targets.
Determination of the angular attitude of earth-bound objects using satellites was considered even prior to the launch of the first GPS satellite [6].Obtaining attitude information from GPS multi-antenna systems has since proven to be a feasible application [7].Most applications of multi-antenna attitude determination are performed using high precision GPS; attitude determination by GPS is often part of an integrated system using both inertial and GPS measurements [8,9].Nevertheless, it has been shown that low-cost systems that provide practical attitude precision are viable [10][11][12], and navigation systems using only GPS for attitude determination have been developed with successful results [11].One of the systems described in [11] uses four GPS antennas to determine azimuth, pitch, and roll, providing an accuracy of approximately 1 mrad (about 0.05 degree) using an antenna baseline distance of 3-5 meters.The two antenna GPS-compass system described in [13] reports results of 4 mrad precision using two off-the-shelf GPS receivers.A four-antenna GPS attitude determination system, with baselines of 11 and 42 m, has been shown to produce precisions of 0.3, 0.6 and 0.9 mrad for pitch, heading and roll respectively as reported in [13].
In this paper, we present a direct georeferencing method for terrestrial LiDAR point clouds via a dual antenna attitude determination scheme.To circumvent the singularity of the dual-antenna system to recover the three orientation angles, we allow the dual-antenna vector to rotate about the LiDAR standing axis to form a multi-vector plane.The formation of the plane allows for the solution of the fundamental rotation problem, where the dual-antenna vectors are defined in the GLCS and in the SOCS at every stop of the system; finding the rotation matrix between SOCS and GLCS defines the rotation matrix.Through simulation and real-world experimentation, we show that our method can provide solutions with a useful level of precision for the three orientation angles.

Motivation
The authors were inspired by work done by Distinguished Professor Marte S. Guitierrez of the Colorado School of Mines in the Philippines in 2006 [14].Figure 1 shows a survey crew installing reflex targets manually on the main slip surface to indirectly geo-reference the LiDAR system.According to [15], "The initial survey required two days as it took about four hours to hike and climb the slide surface, and another four hours to return."The authors of this paper, in collaboration with Riegl USA, are developing a method to make such georeferencing safer and more efficient.

Methodology
In this section, we describe the mathematical model and its linearization technique used in the development of the Dual-Antennae System (DAS).Direct Georeferencing is defined as the computation of the transformation parameters between an arbitrary input coordinate system (represented in this application by the scanner own coordinate system) and a global or mapping coordinate system.The mathematical model used is fundamentally based on the 3D conformal transformation between two orthogonal systems; see [15] for more details while some equations are repeated here for convenience.The 3D conformal transformation from arbitrary scanner coordinates to global coordinates has the form [16]: where: G p X is the position vector of a point in the ground system, s is a scaling factor, G S R is the rotation matrix from the scanner coordinate system to the ground coordinate system, G p X is the position vector of the point in the Scanner coordinate system, and, T is the translation matrix from the scanner coordinate system to the ground coordinate system.
With redundant observations, a least-squares solution can be used to compute s, G S R , and G S T .The traditional approach uses scanned points with known coordinates in the scanner system and in the global system for the two vectors X S and X G respectively.In the proposed method, however, the coordinates (and ultimately the differential vectors) of the phase centers of scanner-mounted GPS antennae are used instead, as will be described later.
The number of angular stops of the antenna-mounted apparatus has to be larger than three to achieve a solution.In order to exploit the higher precision of differential GPS, a two antenna configuration is used.The formulation of this case is as follows: at each stop, two equations for the two antennae are formed as in Equation 2: Subtracting both antennae positions in Equation 2 and assuming a unit scale between the scanner coordinate system and the global coordinate system yields: By assuming small angle rotations in the rotation matrix R and re-arranging, we can rewrite Equation 3 as follows: where ω, φ, κ are the three Euler rotations about the x, y, and z axes of the two coordinate systems, respectively; the small angle rotation assumption is practically possible by leveling and initializing the LiDAR sensor in the north direction.
By applying the identity of skew symmetric matrices and their vector equivalents found in Equation 4 [17], it becomes Equation 5 is linear in the rotation vector (• , • , • ) T .Knowing the measured vectors of the dual-antenna system S X ∆ and G X ∆ , in the scanner and global coordinate systems, respectively , the solution of the angular values • , • , and • can be found using least squares.
Without losing generality, it can also be shown that the small-angle condition above is not necessary to reach Equation 5 from Equation 3. We will proof that through decomposing the rotation matrix into two matrices, one known at initial approximate rotations ) , , ( and another one with unknown small rotations to correct the initial rotations, G S R ∆ .In this case, Equation 3 can be rewritten as: where is defined as in Equation 6but for small correction angles (d• , d• , d• ) T .Using the development as above, Equation 6 can be linearized for the small angle corrections as follows: Where the full rotation angles can be recovered after solving Equation 7, as follows: Equations 5 and 8 are both linear in the unknown parameter vector (• , • , • ), but include observations on both sides of the observation equation.A proper rigorous approach to solving this system of equations requires the general least squares formulation [18].
It should be clear in the reader's mind that the solution we achieve through this algorithm is not equivalent to the epoch-by-epoch solution one achieves from mere computation of orientation angles of the GPS vector.Besides, the rotation of the GPS vector of this algorithm enables overcoming the rank deficiency of the 2D system to calculate 3D orientation angles.The solution we achieve in this algorithm improves the single observation precision by a factor of square root of the number of stops as we will show in the experimental results.

Configuration and Calibration of the Georeferencing System
The dual-antennae georeferencing system was designed for and implemented on a Reigl Z390i laser scanner.A detailed description of the system configuration and its calibration was the subject of [15]; interested readers are directed to the mentioned publication for details while a picture of the experimental setup is depicted in Figure 2

Testing, Results, and Analysis
In this section, we discuss results obtained from simulation and field testing.

Simulation
In order to assert the validity of the proposed georeferencing method, several simulations were performed.The implementation of the simulation allowed for variable parameters which led to better understanding of the problem and avoidance of potential pitfalls of the formulation.The output of the simulation, using parameters that are quite feasible for terrestrial LiDAR scanning, results in solutions with precisions useful for practical applications.Standard deviations of orientation angles of about 0.05 degrees (~1 mrad) were calculated using the dual antenna system.The results of some special-case simulations are discussed in this section.

Single Antenna Case
The simulations for the single antenna case consisted of inputting orientation angles with magnitudes of about 2 degrees and generating antenna positions based on a bar length of about 1 m (for simplicity) and varying numbers of stops.The solution for a 3D conformal transformation using the points was computed via least-squares (Equation 2). Figure 3 is a plot of the number of stops (the number of times the sensor head paused at specific angles, abscissa) versus the precision of the single antenna solution in degrees (ordinate) for multiple simulations with the standard deviation of GPS vectors set at 0.015 m (a realistic estimate for absolute positioning of a single point).The figure displays the asymptotic nature of the precision of the solution with respect to the number of stops.This indicates that after a certain number of stops, the standard deviations of the solution don't get considerably smaller with any more stops.The orientation angle standard deviations for a 30 stop setup were about 0.4 degrees.This level of precision is about ten times worse than what is considered useful.

Dual Antenna Case
The dual antenna simulation was performed by taking user input consisting of the number of stops; the incremental alpha angle (the change in sensor head angle between stops); observable rotation angles omega, phi, and kappa; the length of the baseline between GPS antennas; and the standard deviation of the differential GPS solution.The antenna baseline was assumed to be perfectly aligned with respect to the sensor coordinate system.In other words, calibration errors were considered to be zero in the simulation.The scanner coordinate antenna baseline vectors for each stop were obtained by rotating a vector with magnitude equal to the input length by the input incremental alpha angle at each stop.Next, the GPS coordinate differential vectors were created by applying rotation to the input angles and then adding random Gaussian noise based on the input GPS standard deviation.The solution was then calculated using a least-squares adjustment with output of the standard deviations of the solved angle values (Equation 10).By running multiple simulations with different parameters, we were able to observe how they affect the precision of the orientation solution.

Effect of Number of Stops on Attitude Precision
The standard deviation of the solved orientation angles is greatly affected by the number of stops.By increasing the number of stops, the increased redundancy compensates for the errors caused by GPS vector uncertainty and uncertainty stemming from the use of the approximated skew-symmetric rotation matrix model.The antenna configuration for each simulation was kept symmetric by using an incremental alpha angle equal to 180/(number of stops).As the number of stops approaches infinity, the estimated standard deviations for each angle will approach zero.Because of the symmetry of the configurations, the standard deviations of omega and phi can be considered identical.Figure 3 is a plot of the number of stops (abscissa) versus the precision of the solution in degrees (ordinate) for multiple simulations with bar length of 1 m and standard deviation of GPS vectors at 0.001 m.
As in the single antenna simulation, Figure 4 shows the asymptotic nature of the solution precision with respect to the number of stops.Although a very high number of stops defeats one of the main goals of this project, we can see that precision sharply increases between 2 and 10 stops with a gradual increase between 10 and 50 stops.At 10 stops the standard deviation for the angles is about 0.04 degrees.The plot levels off at about 0.01 as the number of stops approaches 150.It is also notable that kappa is consistently more precise than omega and phi, an inherent property of the design.By comparison with the plot for the precision versus the length of baseline, we can see that the redundancy increase from the increased number of stops has more effect on the precision than the baseline length.This may be due to the fact that the inexactness of the computed angles caused by the skew-symmetric matrix approximation is better compensated for by redundancy than the angle inexactness caused by GPS error is compensated for by increased baseline length.

Effect of Orientation Angle Magnitude on Attitude Precision
Because we are using a skew-symmetric matrix to solve for orientation, the observed angles must be small for this solution to be valid.Thus simulations were run for several angles of increasing magnitude, with the resulting standard deviations as output.Figures 5 and 6 are plots of orientation angles in degrees (abscissa) and resulting standard deviations in degrees (ordinate).For these simulations, the length of the bar was set at 1, the number of stops at 10, and the incremental alpha angle at 18 degrees.The GPS position was assumed to be perfect in order to isolate

Standard Deviation (degrees) Angle Magnitude (degrees)
Omega Phi Kappa errors associated with using the skew-symmetric rotation matrix.Figures 5 and 6 indicate that we should not expect to solve for orientation angles, with our desired uncertainty (0.05 degrees), when the angle magnitude is greater than about 2 degrees.It is also good to note that again, due to the symmetric configuration, estimates for omega and phi are identical.

Effect of Baseline Length on Attitude Precision
In order to study their effect on the precision of the solution, simulations were run using varying baseline lengths.The standard deviation of the GPS position was again set to 0.001 and the magnitude of the observed angles omega, phi and kappa were set at 2 degrees.Figures 7 and 8 are plots of the standard deviations with respect to the baseline length for symmetric configurations of 2 and 10 stops respectively.These plots show that although baseline length plays a large role in the precision of the solution, when the baseline length is greater than 4m, the precision levels off.Thus after this point, the limiting error in the solution of the angles is probably due to uncertainty caused by using the skewsymmetric matrix rotation model.This can be compensated for by increased redundant observations.

Performance Analysis of the Dual Antenna System (DAS)
It was concluded from the simulation results that the single antenna case will not provide sufficient accuracy for terrestrial stationary LiDAR georeferencing.A field test of the dual antenna system was carried out to validate the results obtained from the simulation.The scanner was oriented near north and leveled in order to satisfy the small-angle requirement.Differential GPS observations were collected at 24 stops of the rotating scanner head at 15 degree increments.After removing two obvious blunders, the GPS solutions were then used to solve for the orientation angles via the formulas described in the Method section.
Several different configurations were used to find the solution.Solutions for every symmetric configuration at varying numbers of stops were used to find the angle solution and adjustment standard deviations.For instance, at the two-stop level, there were 15 possible symmetric configurations with 90 degrees between stops: one with stop angles of 15 and 105 degrees; one with stop angles of 30 and 120 degrees; one at 45 and 135 degrees and so on.The solution results are listed in   The solution for the angles has differing values based on the number of stops.Figure 9 shows that the 2 stop configurations result in a large difference between solutions for the different configurations, and as the number of stops is increased, the difference in solution between configurations is decreased.We see the same effect in Figures 10 and 11, the standard deviations plots.As the number of stops increases, the difference in standard deviations between different configurations is reduced.The range of magnitudes of the standard deviations may be due to the geometric configuration of the apparatus with respect to the datum axes and the varying precision of the GPS solution at each stop.In other words, there are situations where fewer numbers of stops result in more accurate precisions due to the above factors.Figure 11, a plot of the precision of omega against the number of stops, illustrates that as the number of stops is increased, the standard deviation converges on a viable value close to our goal of 0.05 degrees.Comparison between these standard deviation plots and the results shown in Figure 4, simulated angle precision versus number of stops, indicates that our field experiment solutions have behaved as the simulations predicted.That is, we are able to observe the predicted asymptotic trend in precision versus number of stops.Also, as indicated in simulations and verified in this real-world experiment, there is a practical number of stops after which the precision stops improving.This is reflected in asymptotic nature of Figures 10,11,and 4. We have also shown, through these experiments, that we can achieve the small angle requirement of our method through leveling and orienting the scanner with north.

External Reference Comparison
In an attempt to characterize the errors of the dual antenna system to an external reference, we mounted a tactical-grade IMU on the LiDAR moving scanner head.The IMU axes were aligned and leveled to the dual antenna bar as much as possible.So, both the IMU and the dual antenna system should experience the same attitude change, except for a small shift.It should be noted that this was not a test of our overall method.This test was used to see how our dual GPS antenna system performs measuring orientation with respect to a datum in terms of tilt/roll and azimuth/heading.
Figure 12 shows the difference between the roll and heading as measured by the IMU and those calculated from the GPS trajectory.It should be noted that some of this error is due to the fact that both systems were not perfectly aligned.It can also be noted that some blunders could not be avoided in the setup and operation of the IMU and the dual antenna system.Angle (Degrees)

Number of Stops
Although the IMU is of medium accuracy and should perform better that a dual antenna system of 0.88 m span separation, Figure 12 shows that there is clear evidence that the heading estimate of the IMU was drifting as a result of a poor alignment.Therefore, the resultant RMSE of the heading error is partially due to the IMU and partially due to the GPS dual antenna system.The calculated RMSE of the tilt error is 0.1 deg with an average error of 0.25 deg, which we consider as the expected accuracy of the dual antenna system in single measurement mode.This on one hand proves the feasibility of the dual antenna system to achieve reasonable accuracy for geo-referencing.Comparing the results from the single measurement RMSE to the achieved accuracy from our method, one can see clearly the advantage of using the new method to estimate the geo-referencing parameters.While our method is capable of achieving better than 0.05 degrees RMSE in estimating the different orientation parameters, the single measurement accuracy is hardly half as good in RMSE with the possibility of committing 0.25 degree error.

Comparison to Conventional Approach
Table 1 show the orientation and shift parameters of the SOCS with respect to GLCS when conventional control targets were used.The orientation results show accuracy better than half a milliradian while shift results show accuracy better than one centimeter.The scale factor of non-unity was the result of using non-metric units for the GLCS coordinates of the control stations, but nevertheless, show accuracy of better than 80 ppm.The solution achieved from the DAS system achieved one half of the conventional accuracy for 12 stops as shown in Figure 9 above.

Summary and Conclusions
Current techniques for georeferencing stationary LiDAR data involve the direct measurement of objects in the point cloud.The collection of external control can be laborious and in some cases impossible due to physical, safety, and other limitations.We have presented an autonomous way of determining the orientation parameters for transforming the point clouds from the Scanner Own Coordinate System to a Global Coordinate System.By using a dual-antenna apparatus mounted on a moving scanner head, we have shown that orientation angle precision of about 0.05 degrees (~1 mrad) can be achieved for the three rotation angles by using a feasible number of angular stops of the scanner head.This is approximately one half of the accuracy achieved with the conventional technique.We presented the theoretical development of a small angle rotation case and supplanted it with the general case, without losing generality.The length of the apparatus bar did not seem to have affected the attained accuracy substantially leading to the conclusion that a reasonable bar length will suffice in most applications.
The proposed apparatus and algorithm are simple and easy to implement.Currently, tests indicate that at small angles the linear (proposed) and nonlinear solution differences are minimal.Also, solving the model using the general least-squares technique does not improve the accuracy significantly.We are planning more rigorous tests of the system.We will include multiple setups and data from a camera to increase the redundancy and the internal consistency with the conventional technique.This work paves the way for augmentations of the presented dual antenna method with other techniques including Kalman filtering.We also plan to use measurements from the scanner's horizontal rotation controller as updates to the Kalman filter.

Figure 1 .
Figure 1.LiDAR survey crew (left) and manual survey crew installing the reflex targets on the slip surface of the February 17, 2006 landslide in Guinsaugon, Southern Leyte, Philippines, photos courtesy Professor Guitierrez and Riegl USA. below.

Figure 2 .
Figure 2. Experimental Setup of the Dual-Antenna System (DAS): Two GPS antennae fixed on a metal bar mounted on top of the LiDAR rotating sensor head.The Inertial Measuring Unit (IMU) in the middle of the setup is for check purpose only and is not part of the DAS system.

Figure 3 .
Figure 3. Simulation Results of the Single-Antenna Case showing precision of calculated orientation angles vs. number of device stops; the curve practically starts to dampen at about 20 stops, however the achieved accuracy is not close to what is required for stationary LiDAR.

Figure 4 .
Figure 4. Simulation Results of the Dual-Antenna Case showing precision of calculated orientation angles vs. number of device stops; the curve's actual dampening happens at close to 100 stops; the achieved accuracy is twice as high as what is required for stationary LiDAR.

Figure 5 .
Figure 5. Simulation Results of the Dual-Antenna Case showing precision of calculated orientation angles vs. their magnitude; as the orientation angle gets larger than 5 degrees, the assumption of small angle rotation gets violated and small-angle model becomes invalid.

Figure 6 .
Figure 6.Simulation Results of the Dual-Antenna Case showing precision of calculated orientation angles vs. their magnitude over the first 5 degrees; the 1 mrad precision is possible with orientation angles smaller than 3 degrees.

Figure 7 .Figure 8 .
Figure 7. Simulation Results of the Dual-Antenna Case showing precision of calculated orientation angles vs. mounting bar length for systematic two device stops; the curve's actual dampening happens at about 1 m bar length which show the practicality of the method; it can also be noticed that the achieved accuracy at the 1 m bar length value is close to what is required for stationary LiDAR but not enough with two device stops.

Figure 9 .
Figure 9. Field Test Results of the DAS showing values of calculated orientation angles vs. number of device stops for different symmetrical configurations of the device stops; the chart shows that configuration plays a role in determining the value of the calculated angle.

Figure 10 .
Figure 10.Field Test Results of the DAS showing precision of calculated orientation angles vs. number of device stops for different symmetrical configurations of the device stops; the chart shows the asymptotic nature of the precision curve as the number of device stops approaches 10 stops.

Figure 11 .
Figure 11.Field Test Results of the DAS showing precision of calculated orientation angle Omega vs. number of device stops for different symmetrical configurations of the device stops; the chart shows the asymptotic nature of the precision curve as the number of device stops approaches 10 stops.

Figure 12 .
Figure 12.GPS Tilt and Azimuth Errors Compared to calculated instantaneous orientation angles using the mounted Inertial Navigation System (INS); one milli-radian accuracy is not achievable from a single GPS DAS observation justifying the need for multiple setups and the need for the algorithm proposed.

Table 1 .
Results of a Total Station survey of control targets and corresponding orientation parameters of the LiDAR sensor head compared with respect to the global coordinate system.