RIGOROUS LIDAR STRIP ADJUSTMENT WITH TRIANGULATED AERIAL IMAGERY

This paper proposes a POS aided LiDAR strip adjustment method. Firstly, aero-triangulation of the simultaneously obtained aerial images is conducted with a few photogrammetry-specific ground control points. Secondly, LiDAR intensity images are generated from the reflectance signals of laser foot points, and conjugate points are automatically matched between the LiDAR intensity image and the aero-triangulated aerial image. Control points used in LiDAR strip adjustment are derived from these conjugate points. Finally, LiDAR strip adjustment of real data is conducted with the POS aided LiDAR strip adjustment method proposed in this paper, and comparison experiment using three-dimensional similarity transformation method is also performed. The results indicate that the POS aided LiDAR strip adjustment method can significantly correct the planimetric and vertical errors of LiDAR strips. The planimetric correction accuracy is higher than average point distance while the vertical correction accuracy is comparable to that of the result of aero-triangulation. Moreover, the proposed method is obliviously superior to the traditional three-dimensional similarity transformation method.


INTRODUCTION
Airborne Light Detecting and Ranging (Airborne LiDAR) is a commonly used system for the fast acquisition of topographic information, which is also known as Airborne Laser Scanning (ALS).The three components of LiDAR are Laser Scanner, Global Positioning System (GPS), and Inertial Measurement Unit (IMU).The GPS and IMU are highly integrated to compose the Positioning and Orientation System (POS).The LiDAR data collection is usually conducted in strip-wise pattern, and the LiDAR data obtained during one flight strip is called as a LiDAR strip.Due to its character as a combined system, the coordinates of points collected by LiDAR suffer from many systematic errors and some random errors.According to previous researches, the major systematic errors include range errors, mirror angle scale errors, mounting errors, and POS systematic errors (Baltsavias, 1999;Filin, 2003;Habib et. al., 2009;Habib et. al., 2010;Kumari et. al., 2011).There will be discrepancies among overlapping LiDAR strips and between LiDAR strips and ground truth if these systematic errors haven't been eliminated.
Various LiDAR Strip Adjustment (LSA) methods have been conducted for the aim of eliminating these discrepancies.Existing approaches can be classified into two main categories: data-driven methods and sensor system-driven (Filin, 2003;Lee et. al., 2007;Willers et. al., 2008;Habib et. al., 2009).Datadriven methods apply translation and rotation models in strip adjustment, e.g. a six-parameter rigid-body transformation (three shifts and three rotations) (Lee et. al., 2007;Vosselman, 2008;Habib et. al., 2009;Habib et. al., 2010;Rentsch and Krzystek, 2012) or even a simple vertical shift (Crombaghs et. al., 2000;Willers et. al., 2008).In order to provide both planimetric and vertical constrains for LiDAR strips, conjugate features such as lines and surfaces are matched and then colinearity or coplanarity relationship between these features are adopted to calculate the translation and rotation parameters.
Due to the sophisticated affection from various systematic errors to the LiDAR point coordinates, the consequential distortion of LiDAR strips can not be simply assumed as linear (Maas, 2002;Filin, 2003).So the adoption of linear models such as translation and rotation cannot eliminate the discrepancies completely, and rigorous approaches belonging to the sensor system-driven category are required.Sensor systemdriven methods adopt LiDAR geo-location equations as adjustment model.These methods require the original observations (GPS, IMU and the laser measurements) or at least the trajectory and the time-tagged point clouds (Filin, 2003;Kumari et. al., 2011).Among most sensor system-driven methods, only systematic errors of mounting parameters and ranging are calibrated (Burman, 2000;Filin, 2003;Skaloud and Lichti, 2006;Favalli et. al., 2009;Kumari et. al., 2011;Hebel and Stilla, 2012).They believed that the state-of-art GPS/IMU on board has a sufficient precision, so the POS system errors are negligible.Nevertheless, owing to the high price and some other reasons, a lot of LiDAR system used in China are equipped with low precision POS.If other negative factors are considered, such as weak GPS constellation geometry, inadequate reference GPS stations, and even the limited professional ability for GPS/IMU data post-processing, the POS system errors can be very large in magnitude and will be the major error sources of LiDAR strip distortion.So the calibration of POS systematic errors is very essential for strip adjustment.
Both sensor system-driven and data-driven methods need control elements or tie elements.Current used control or tie elements are mainly lines and surfaces (Maas, 2002;Filin, 2003;Kager, 2004;Skaloud and Lichti, 2006;Lee et. al., 2007;Habib et. al., 2009;Habib et. al., 2010;Kumari et. al., 2011;Hebel and Stilla, 2012).Considering many LiDAR systems are equipped with high resolution digital cameras, which can obtain high resolution digital imagery simultaneously during the collection of LiDAR strips.If a few photogrammetry-specific ground control points are measured, POS aided bundle block adjustment can be conducted to obtain the precise exterior orientation elements of aerial images (Zhang et. al., 2012).So the aero-triangulated aerial images can provide control information for the LSA by matching conjugate features between LiDAR data and aerial imagery.
Based on the above analysis, this paper proposed a rigorous POS aided LSA method.A flowchart of this paper is given in Fig. 1 The coordinate of laser foot point P in local reference frame can be derived by the LiDAR geo-location equation (Filin, 2003;Skaloud and Lichti, 2006;Favalli, et. al., 2009;Habib et. al., 2009;Habib, et. al., 2010;Kumari, et. al., 2011;Hebel and Stilla, 2012) as follow: Where [X P(t) , Y P(t) , Z P(t) ] is the coordinate of point P in local reference frame at time t, [X G(t) , Y G(t) , Z G(t) ] is the coordinate of the POS system central in local reference frame at time t, R I(t) is the rotation matrix composed by the attitude angles of the POS system in local reference frame at time t, [X L , Y L , Z L ] is the lever-arm bias of the laser scanner in POS body frame, R MIS_L is the rotation matrix of the bore-sight angles of the laser scanner in POS body frame, R R(t) is the rotation matrix of laser beam angles in laser scanner frame at time t, d(t) is the range between laser fire point and foot point.
Since the POS equipment and the laser scanner are highly integrated and there would be regular calibration on their mounting errors, during one flight mission the lever-arm bias [X L , Y L , Z L ] and the rotation matrix R MIS_L of the bore-sight angles in eq. ( 1) are fixed value and the lever-arm bias errors can be assumed as very small.After manufacture's calibration, the measurement error for range and laser beam angles of laser scanner could be negligible.Additionally, there have been researches revealing that the POS angular errors, mounting angular errors and scan angle of laser scanning are highly correlated and cannot be introduced as unknowns simultaneously in the estimation process (Kumari et. al., 2011).Therefore, only POS position and attitude angular errors are introduced in the LSA.Based on these analyses, eq. ( 1) can be simplified as follows: where, [X P(t) , Y P(t) , Z P(t) ], [X G(t) , Y G(t) , Z G(t) ] and R I(t) share the same definition in eq.( 1), [X LPOS(t) , Y LPOS(t) , Z LPOS(t) ] is the coordinate of LiDAR point P in the POS body frame at time t, which is given in eq. ( 3 Owing to the POS systematic errors, there exists discrepancy between the observation value (X S(t) , Y S(t) , Z S(t) , ω S(t) , φ S(t) , κ S(t) ) and the true value (X ).The systematic errors of the position observations and the angular observations of the POS system are time dependent (Zhang et. al., 2012).And Kager (2004) uses time dependent polynomials to correct POS observation values.In this paper second order polynomials are adopted to fit the systematic errors of the POS observations.Such a relationship is described mathematically in Where t is the time when the laser point is obtained, (a is respectively the constant coefficient, the 1st order coefficient and the 2nd coefficient of the POS systematic error model, and t 0 is the reference time (usually the average time of strip).
For a certain LiDAR point obtained at time t, its coordinates [X P(t) , Y P(t) , Z P(t) ] in local reference frame can be found in LiDAR strip file.The POS observations at time t can be interpolated from trajectory.Then the point's coordinates [X LPOS (t) , Y LPOS(t) , Z LPOS(t) ] in POS body frame can be calculated using eq.( 2), which are used as original information of this point.If its true ground coordinates are known, this point can be used as a control point to calculate the systematic errors of the POS observations.If the planimetric coordinates and vertical coordinate are accurate, then this point can be viewed as a planimetric-vertical control point.If only planimetric coordinate is accurate, it can be viewed as planimetric control point.
Taken the true coordinates of the control points [X P(t) , Y P(t) , Z P(t) ] in local reference frame as observation value, their coordinates [X LPOS (t) , Y LPOS(t) , Z LPOS(t) ] in POS body frame are known numbers while the constant term, the 1st order term and the 2nd term of the POS systematic error model are unknowns.Integrate eq.4 with eq.2.As the angular related unknowns can not be separated form each other, the equations are linearized and the error equation is shown as follow: which is the unknowns to be calculated in strip adjustment, A is the coefficient matrix of the unknowns, L is the constant term of the error equation, P is the weight matrix of the ground control points.Eq. ( 5) includes three equations corresponding with X, Y, Z coordinates respectively.Precisely speaking, the unknowns r are not equal to real POS systematic error model coefficients, because they may be affected by mounting angular errors.
All the control points can help to build the error equations and find the solution to the unknowns using least square method.
After get the coefficients of the POS systematic error model for each LiDAR strip, the trajectory can be rectified by eliminating the POS systematic errors, and the rectified coordinates of LiDAR points in local frame can be recalculated using eq.( 2).

Automatic Matching of Conjugate Points among Airborne LiDAR Strips and Aerial Images
LiDAR intensity images are generated from reflectance signals of laser foot points, with the TIN (Triangulated Irregular Network) interpolation method (the sampling distance is usually about 1/3 of the point distance).And conjugate points are matched automatically among LiDAR intensity images and aerial images.Then, control points for LSA can be obtained from the successfully matched conjugate points.
In this paper, before the automatic conjugate points matching, the first procedure is to automatically match the tie points on the aerial images aided by POS data, and POS aided bundle block adjustment of the aerial images are complemented using a few photogrammetry-specific ground control points.Then, down sample aerial images to 3 times (the down sampled images are called 1:3 aerial images here) and extract Harris feature points on the 1:3 images.Each Harris feature point is matched with the LiDAR intensity images to acquire the conjugate point.During the matching process, an affine transformation model is adopted to correct the distortion between aerial image and LiDAR intensity image.If successfully matched, automatically match the conjugate points on the overlapped 1:3 aerial images.After calculate the coordinates of the feature point (X I , Y I , Z I ) in local reference frame using space intersection, use the coordinates as a reference.The coordinates (X D , Y D , Z D ) in local reference frame of the conjugate point in LiDAR data can be calculated from its coordinates in LiDAR intensity image.The point-pair (X I , Y I , Z I ) and (X D , Y D , Z D ) is used as control point in LSA.The principal to judge the attribute of one control point is: Draw a circle around the point with a certain radius in LiDAR strip, get all the elevation information of the points fall in the circle and calculate the difference between the maximum and the minimum.If the difference is less than the threshold (0.1m), the control point is located on a flat area and it can be viewed as a planimetric-vertical control point.Otherwise the interpolated elevation of the feature point might consist of large error therefore merely qualified as a planimetric control point.

Experimental Data
To verify the correctness of the presented POS aided LSA method, the airborne LiDAR point cloud data of three strips in a certain region (named Las1, Las2, and Las3 respectively) and simultaneously obtained digital aerial images are used in this paper to test the proposed method.The set of LIDAR data points was captured at a mean relative flight height of 700 m with a mean point density around 4.8 points/m 2 and the average point distance is 0.50m.The aerial image has a frame of 8984×6732 pixels, the CCD pixel size is about 0.006mm, the focus of the camera is 51.695mm and the ground resolution is around 0.09m.The frame rate of camera is about 5.5 seconds.

Results of Aero-triangulation
28 photogrammetry-specific ground feature points are surveyed using GPS on the test field, and their coordinate accuracy is around 0.05m.After manually measured their image point coordinates on the aerial images, select five ground points located on the four corners and the center as control points and leave the rest 23 points as check points in POS aided bundle block adjustment.The Root Mean Square Error (RMSE) of the planimetric and vertical coordinate residues of the control points after POS aided bundle block adjustment is 0.05m and 0.07m respectively, the RMSE of the planimetric and vertical coordinate residues of the check points is 0.17m and 0.11m respectively.The results show that the POS aided bundle block adjustment achieves good accuracy.

Results of Automatic Matching among Aerial Imagery and LiDAR Strips
In the experiment, the LiDAR intensity images are generated from LiDAR point cloud at a sampling interval of 0. There is certain scale and angular discrepancy between LiDAR intensity images and the 1:3 aerial images.However, as shown in Fig. 2, adopting the proposed conjugate points matching method can well extract conjugate points like zebra strips' feature points and building roof corner points.Because the 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 initial matching position can be predicted very precisely using the exterior orientation parameters of aerial images, the number of incorrect matching is small.Moreover, the incorrect matching can be removed in the LSA procedure using iteration method with variable weights.To analyze the planimetric coordinate accuracy of each LiDAR strip, 37 pairs of conjugate points are manually measured accurately among the original aerial images and the LiDAR intensity images as planimetric check points, mainly the obvious feature points like zebra strips' endpoints and roof corner points.There are 10, 14, and 13 planimetric check points on Las1, Las2, and Las3, respectively.According to the corresponding image point coordinates, the coordinates in local reference frame of the planimetric check points are calculated through space intersection using the image exterior orientation parameters obtained in aerotriangulation.Since the accuracy of aerotriangulation is very good, the coordinate in local frame can be regarded as the ground truth and used as a reference.Then, the planimetric coordinates of the check points in LiDAR strip are calculated using its image coordinates in LiDAR intensity images.The differences between the reference planimetric coordinates and the planimetric coordinates in LiDAR strip, also known as the check point residuals, are computed., dX and dY represent the X coordinate residual and the Y coordinate residual respectively) of the planimetric check points before and after adjustment.As shown in the table, the RMSE of planimetric coordinate residuals of the planimetric check points on three LiDAR strips are 1.57m, 0.71m and 1.68m respectively before adjustment.Except for Las2, the point clouds strips all have big planimetric coordinate error.After TDST LSA, the RMSE of planimetric coordinate residuals of the planimetric check points has been reduced to 0.70m, 0.35m and 0.96m respectively.Only the planimetric coordinate error of Las2 has been well corrected, the other two strips still have big errors.After adopted proposed method in this paper, the three RMSEs reduced to 0.36m, 0.25m and 0.22m respectively.Considering that the average point distance of the LiDAR point cloud is around 0.5m, the planimetric coordinate error correction effects of the three strips are generally ideal.
3.2.3.2Vertical Accuracy Analysis before and after Adjustment (1) Vertical discrepancies between overlapping LiDAR strips To evaluate the vertical discrepancies between the overlapping LiDAR strips, several positions between the overlapping LiDAR strips are selected to compute their vertical discrepancies between overlapping strips.Specifically, select a point location every 5m along the flight direction in the overlapping area and search the elevation of all the LiDAR points within certain radius around the point (the radius is set to be 1.2m in this paper).If the difference between the maximum and the minimum elevation is less than the threshold (0.1m), the point is located on a flat area and it can be viewed as a vertical check point.For every vertical check point, its elevation can be calculated through interpolation on the two overlapped LiDAR strips.The discrepancies of the elevation value between the two strips can be used to evaluate the vertical discrepancies between the LiDAR strips.Fig. 3 and Fig. 4 show the discrepancies between Las1 and Las2, Las2 and Las3 respectively before and after performing LSA.As shown in Fig. 3, after POS aided LSA, the vertical discrepancies have obviously reduced to the magnitude between -0.20m and 0.05m.This is approximately comparable to the accuracy of aero-triangulation.From Fig. 4, after POS aided LSA, the vertical discrepancies have obviously reduced to the magnitude between -0.22m and 0.08m, and most of them are around 0.0m.Seen from the experiment results of Fig. 3 and Fig. 4, adopting the proposed POS aided LSA method can well improve the vertical discrepancies between the LiDAR strips.
(2) Analysis of the Absolute Vertical Accuracy of the LiDAR Strips before and after Adjustment To further analysis the absolute vertical accuracy of each LiDAR strips, 72 tie points located in flat region is precisely measured on the aerial images.Calculate their coordinates in local frame through the image exterior orientation parameters by space intersection as ground truth and take them along with the above 28 photogrammetry-specific ground points as the vertical check points of the LiDAR strip.There are respectively 20, 25, 28 qualified points within the coverage of Las1, Las2 and Las3.The point clouds vertical accuracy can be measured through calculating the discrepancies between the ground truth of the check point and its interpolated vertical value in the LiDAR strip.The statistics of the mean, RMSE and the maximum of the vertical coordinate residuals of the vertical check points on three LiDAR strips before and after LSA are shown in Table 3.According to the corrected results of the three LiDAR strips, the vertical precision of the point cloud after adopting the POS aided LSA is up to 0.1m.The largest vertical residual is smaller than 0.2m, which shares the same level of accuracy with aerotriangulation.Obviously, the POS aided correction is better than the TDST LSA.data are performed using the proposed POS aided LSA method, and comparision experiment using traditional TDST LSA method is also conducted.Results show that our proposed method achieves very good improvement of both planimetric and vertical accuracy and is superior to the traditional TDST LSA method.
the coefficients of the POS systematic error model, 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 2m and 9331 pair of conjugate points are matched between the 1:3 aerial images and LiDAR intensity images.Among them, 2544 pairs are on Las1 including 918 pairs of planimetric-vertical control points and 1626 pairs of planimetric control points, 3724 pairs are on Las2 including 1096 pairs of planimetricvertical control points and 2628 pairs of planimetric control points, 3063 pairs are on Las3 including 906 pairs of planimetric-vertical control points and 2157 pairs of planimetric control points.Part of the matching result is shown in Fig. 2.
Part of conjugate points matched between the 1:3 aerial images and LiDAR intensity images (There are two sets of matching results (a), (b), the upper half part of each set is the 1:3 aerial image feature point and the bottom half is the conjugate point matched on LiDAR intensity, the cross indicates their locations) 3.2.3Results of LSA 3.2.3.1 Analysis of the LiDAR Strip Planimetric Coordinate Accuracy before and after Adjustment

Figure 3 .
Figure3.Discrepancies between Las1 and Las2 before and after LSA (In Fig.3and Fig.4, the horizontal axis means the vertical check point selected along the flight direction on the LiDAR strips.The vertical axis presents the vertical discrepancies of the vertical check points between overlapping strips, and the unit is meter.The results demonstrated in blue, pink and green are the results before LSA, after using POS aided LSA and after using TDST LSA)

Figure
Figure 4. Discrepancies between Las2 and Las3 before and after LSA

2. METHODOLOGY 2.1 POS Aided LSA
. Aerotriangulation of aerial images is firstly conducted with a few Photogrammetry-Specific Ground Control Points (P-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 SGCP) (usually five P-SGCPs in the whole surveying field), and the precise exterior orientation parameters of aerial images are obtained.LiDAR intensity images are generated from the reflectance signals of laser foot points, and conjugate points are matched between the LiDAR intensity images and the aerotriangulated aerial images.Control points used in LSA are derived from these conjugate points.LSA of real data is conducted with the proposed POS aided LSA method, and comparison experiment using Three-Dimensional Similarity Transformation (TDST) LSA method is also conducted.

Table 1 .
Each strip consists of 28 images, so the flight time of each strip is about 2.6 minutes.And the laser foot point size is smaller than 0.35m.The raw LiDAR data is obtained with Trimble Harrier 56 LiDAR system, which is equipped with POS AV 510.Systematic and random errors added to the POS data of 3 LiDAR strips (unit: rad)

Table 2 .
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 The planimetric coordinate residuals of check points before and after correction (units：m)

Table 3 .
4. Discrepancies between Las2 and Las3 before and after LSA The vertical coordinate residuals of vertical check point before and after LSA (units：m) LiDAR strips accuracy can be affected by many systematic errors, such as range errors, mirror angle scale errors, mounting errors, and POS systematic errors.However, apart from POS systematic errors, the other systematic errors can be calibrated in laboratory by the manufactures.Moreover, many LiDAR systems used in China are equipped with low precision POS for the economical and other reasons.Thus the POS systematic errors single out to be the main error sources for postprocessing.This paper proposed a rigorous POS aided LSA method.Aerotriangulation of aerial images are conducted with a few photogrammetry-specific ground control points.Then the aero-triangulated aerial imagery are used as control data for LSA by implementing feature point matching between aerial imagery and LiDAR intensity imagery.Strip adjustment of real ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, VolumeII-5/W2, 2013  ISPRS Workshop Laser Scanning 2013, 11 -13 November 2013, Antalya, Turkey