Advanced
Strip Adjustment of Airborne Laser Scanner Data Using Area-based Surface Matching
Strip Adjustment of Airborne Laser Scanner Data Using Area-based Surface Matching
Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography. 2014. Dec, 32(6): 625-635
Copyright © 2014, Korean Society of Surveying, Geodesy, Photogrammetry and Cartography
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.
  • Received : November 11, 2014
  • Accepted : December 12, 2014
  • Published : December 31, 2014
Download
PDF
e-PUB
PubReader
PPT
Export by style
Share
Article
Author
Metrics
Cited by
TagCloud
About the Authors
Dae Geon Lee
Department of Geoinformation Engineering, Sejong University, Korea (E-mail:dglee@sju.ac.kr)
Eun Jin Yoo
Member, Department of Geoinformation Engineering, Sejong University, Korea (E-mail:ejyoo@sju.ac.kr)
Jae-Hong Yom
Member, Department of Geoinformation Engineering, Sejong University, Korea (E-mail:jhyom@sejong.ac.kr)
Dong-Cheon Lee
Corresponding Author, Member, Department of Geoinformation Engineering, Sejong University, Korea (E-mail:dclee@sejong.ac.kr)
Abstract
Multiple strips are required for large area mapping using ALS (Airborne Laser Scanner) system. LiDAR (Light Detection And Ranging) data collected from the ALS system has discrepancies between strips due to systematic errors of on-board laser scanner and GPS/INS, inaccurate processing of the system calibration as well as boresight misalignments. Such discrepancies deteriorate the overall geometric quality of the end products such as DEM (Digital Elevation Model), building models, and digital maps. Therefore, strip adjustment for minimizing discrepancies between overlapping strips is one of the most essential tasks to create seamless point cloud data. This study implemented area-based matching (ABM) to determine conjugate features for computing 3D transformation parameters. ABM is a well-known method and easily implemented for this purpose. It is obvious that the exact same LiDAR points do not exist in the overlapping strips. Therefore, the term “conjugate point” means that the location of occurring maximum similarity within the overlapping strips. Coordinates of the conjugate locations were determined with sub-pixel accuracy. The major drawbacks of the ABM are sensitive to scale change and rotation. However, there is almost no scale change and the rotation angles are quite small between adjacent strips to apply AMB. Experimental results from this study using both simulated and real datasets demonstrate validity of the proposed scheme.
Keywords
1. Introduction
Over the last two decades LiDAR data has been one of the major sources for establishing geospatial information infrastructure including DEM, topographic maps, urban model, GIS database, biomass measurement, and so on. LiDAR data is collected in strip-wise with multiple strips to scan entire areas. In order to obtain directly geo-referenced point cloud data, laser scanner, GPS and INS are integrated. However, systematic errors of the individual devices and inaccurate integration of the multi-sensor ( i.e ., insufficient system calibration or mounting bias) lead discrepancies in the overlapping areas between adjacent strips. Therefore, strip adjustment to remove or minimize the discrepancies is required to create seamless and compatible data. Precise calibration requires original observations from GPS/INS and laser range measurements, which are not available to the end users in most cases ( Bang ., 2009 ; Kersting and Habib, 2012 ; Rentsch and Krzystek, 2012 ). Instead, 3D coordinates of the LiDAR point data are provided. Therefore, the strip adjustment is based on the 3D coordinate transformation between overlapping strips. The LiDAR script adjustment is similar to independent model adjustment of the aerial triangulation. Therefore, the adjustment parameters are basically composed of scale, translations, and rotations between adjacent strips.
Strip adjustment can be also utilized as internal quality evaluation of the LiDAR data associated with quality of the system calibration even though there is no unique and standard strip adjustment method yet ( Habib ., 2008 ; Pfeifer and Bӧhm, 2008 ; Schenk, 2001 ). The central issue of the strip adjustment is to define corresponding features and to find locations of the features in the overlapping strips. The corresponding features could be points, lines or planes ( e.g. , corners, edges or roof facets of buildings). Several studies have carried out to improve LiDAR data quality through the strip adjustment. Crombaghs . (2000) proposed a method for reducing vertical discrepancies between overlapping strips. Habib . (2008) utilized linear feature matching to solve corresponding problem in the overlapping areas. Kilian . (1996) introduced a procedure similar to the aerial photo strip adjustment for detecting discrepancies and improving the compatibility. Maas (2000) suggested a method to establish correspondence by minimizing the normal distances between discrete points in one LiDAR strip and TIN patches in the other strip.
In this paper, the conjugate point features were identified using correlation template matching with sub-pixel accuracy. Template matching takes characteristics of the surface into account. The method finds similarity of the surface patches in the overlapping strips in terms of cross-correlation instead of point-to-point correspondence. The template matching is quite straight forward and easy to implement. The strip adjustment was performed with 3D similarity transformation based on the identified conjugate templates. Finally, the accuracy of the strip adjustment was evaluated by recovery of the discrepancy and RMSE (Root Mean Square Error) of the adjustment parameters, and RMSE of LiDAR point coordinates as well as visual inspection of the profiles after strip adjustment. The results show that the proposed procedure is suitable for airborne LiDAR strip adjustment.
2. Surface Matching for Conjugate Features
- 2.1 Overview of study
LiDAR strip adjustment procedure consists of three parts:
  • (1) Determination of the correspondence between overlapping strips
  • (2) Computation of the adjustment parameters and coordinate transformation
  • (3) Quality assessment of the parameters and transformed coordinates after adjustment
Correspondence problem is associated with conjugate feature matching. Strip adjustment reduces discrepancies by 3D coordinate transformation. Finally, quality of the strip adjustment is evaluated by checking compatibility of the LiDAR point clouds. Fig. 1 shows the overview of the proposed method.
PPT Slide
Lager Image
Work flow of strip adjustment
To find the conjugate points might not be suitable for LiDAR data-derived surfaces because it is almost impossible to identify conjugate points in the overlapping areas. Identification of the distinct points is quite difficult and not reliable due to the irregular nature of the LiDAR footprints. Therefore, it is obvious that most probably the exact same LiDAR points do not exist in the overlapping strips ( Habib ., 2008 ). In other words, point-to-point correspondence could not be valid as shown in Fig. 2(a) . Therefore, the term corresponding feature should reflect topographic characteristics of the surface of the interest regions. Template matching depicted in Fig. 2(b) is one of the traditional methods to find similarity between features.
PPT Slide
Lager Image
LiDAR point distribution and matching conjugate features
Advantages of the template matching are straight forward, simple and practical to implement while the major drawback is unreliable under scale change, rotations and geometric distortions. In addition, irregularly distributed points are interpolated to create regular grid data. The templates of the reference and target surfaces implicitly represent surface characteristics, and the similarity between surfaces is computed by cross-correlation coefficients. Therefore, template matching is basically not based on point-to-point correspondence but surface-to-surface correspondence. Some sophisticated matching methods utilize LiDAR dataderived secondary products such as lines or planes that require further processes including extraction of the line features, segmentation of the LiDAR points, grouping and modeling. In consequence, quality of the correspondence would be affected by accuracy of the secondary products.
- 2.2 Automatic extraction of conjugate features
It is not necessary that the matching features are to be specific objects in template matching. However, to find reliable matching features such as well-defined points in the homogeneous areas is difficult. Therefore, to improve matching accuracy and efficiency, selection of the reliable conjugate feature is recommendable. Uniquely identifiable features would be reliable matching candidates to avoid ambiguous or multiple matching. Reliable matching candidates are distinct model key points such as corners of the buildings or other man-made structures. In this study, Canny edge operator ( Canny, 1986 ) and Harrison corner detector ( Harris and Stephens, 1988 ) described from Equation 1 to Equation 4, were used to extract conjugate features. This procedure aims at automatic detection of the most possible conjugate features. Locations of the extracted conjugate feature are to be reference for determining search surfaces.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
The coordinates of a corner point in the local window are computed by:
PPT Slide
Lager Image
where M is a gradient matrix computed from partial derivatives with respect to X and Y. Z ( X,Y ) denotes Z-coordinate at ( X,Y ). w is local window size. K and T denote threshold values. Rij is R value at ( i,j ) in the local windows.
- 2.3 Determination of conjugate features
Incorrectly determined conjugate features deteriorate strip adjustment ( Pfeifer ., 2005 ). Template size and locations of the search surface are significant factors in the template matching. The size and locations depend on the point density (or average ground distance) and tolerance of the GPS/INS accuracy. The centers of the templates are not necessarily distinct points ( e.g. , building corners). However, it is recommended to avoid locating templates in homogeneous areas because to determine the maximum correlation is ambiguous in such area. The cross-correlation in the overlapping strips is computed by Equation 5.
PPT Slide
Lager Image
where n and m denote the template size ( i.e., n×m ), Z T ( Xi,Yi ) is Z-coordinate at ( Xi,Yi ) in template, ZM ( Xi,Yi ) is Z-coordinate at ( Xi,Yi ) in the matching windows inside of the search surface,
PPT Slide
Lager Image
and
PPT Slide
Lager Image
are averages of the Z-coordinates in template and matching window, respectively.
The correlation surfaces were formed by fitting mathematical functions. To choose an appropriate mathematical function is crucial since shape of the surfaces depend on the fitting functions. Fig. 3 shows two correlation surfaces. Fig. 3(a) and Fig. 3(b) are results from 3rd-order polynomial and bi-harmonic fitting in a search window, respectively. Since a conjugate feature with sub-pixel accuracy is to be located near maximum correlation, polynomial fitting using Equation 6 might not be appropriate representation of the correlation surfaces. Instead, the bi-harmonic function according to Equation 7 renders the correlation surface reasonably compared with polynomial surface.
PPT Slide
Lager Image
PPT Slide
Lager Image
PPT Slide
Lager Image
Correlation surfaces
where n and m are order of the polynomial. aij is a coefficient of the polynomials. ∇ 2 denotes Laplacian operator that performs sum of the second derivatives of the correlation surface with respect to each variable, X and Y .
Location of the conjugate feature in each search surface occurs at the maximum correlation. However, the location with maximum correlation might not guarantee to be a conjugate point. Therefore, sub-pixel accuracy should be considered by analyzing correlation surfaces to find accurate matching location. The sub-pixel coordinates of the matching locations could be computed by partial derivatives of the fitting surfaces. i.e. , sub-pixel locations of the maximum correlation were estimated by partial derivatives with respect to X and Y where ∂ρ/∂ X = 0 and ∂ρ/∂ Y = 0.
3. LiDAR Strip Adjustment
- 3.1 Computation of adjustment parameters
As shown in Fig. 4 , relationship between adjacent strips could be determined by 3D similarity transformation having three translations ( XT, YT, ZT ), three rotation angles ( ω, φ, κ ) and a scale ( S ) based on Equation 8.
PPT Slide
Lager Image
3D coordinate transformation between strips
PPT Slide
Lager Image
where
PPT Slide
Lager Image
= [X 2 , Y 2 , Z 2 ] T ,
PPT Slide
Lager Image
= [X T , Y T , Z T ] T ,
PPT Slide
Lager Image
= S·R·[X 1 , Y 1 , Z 1 ] T with S is scale, and R denotes 3D rotation matrix defined by ω, φ, κ. ( X1, Y1, Z1 ) and ( X2, Y2, Z2 ) are LiDAR points in strip 1 and strip 2, respectively. Therefore, the coordinate transformation can be performed through the 3D similarity transformation. It is noted that actually there is no scale change between strips, i.e. , S=1 ( Filin and Vosselman, 2004 ).
In order to solve the transformation parameters, initial approximations of the parameters, linearization of the transformation, and iteration of the computation until reaching required accuracy are necessary. Linearized form of the 3D similarity transformation is written as Equation 9.
PPT Slide
Lager Image
where subscript i indicates iteration step. ( Xo, Yo, Zo ) denote transformed coordinates that are computed using approximate parameters in previous iteration ( i.e. i-1 ) The linearized rotation matrix is given by equation 10.
PPT Slide
Lager Image
The linearized observation equation for least-square solution is written as Equation 11.
PPT Slide
Lager Image
where [ dXT, dYT, dZT, dS, dω, dφ, dκ ] T i are updating increments of the parameters at i -th iteration, and (e X , e Y , e Z ) is error vector. The matrix form of the observation equation is shown in Equation 12.
PPT Slide
Lager Image
The transformation parameters by least-squares are computed by Equation 13:
PPT Slide
Lager Image
where ε is permissible tolerance, and P is weight matrix. Finally, the iteration is to be terminated if (
PPT Slide
Lager Image
) ≤ ε ,
- 3.2 Accuracy assessment
Accuracy of the results was analyzed based on the “variance component” (Equation 14), “dispersion” of the parameters (or RMSE of the parameters) (Equation 15), and “RMSE” of the adjusted coordinates (Equation 16).
PPT Slide
Lager Image
where
PPT Slide
Lager Image
n and m are number of conjugate points and adjustment parameters, respectively. Dispersion indicates reliability of the parameters.
PPT Slide
Lager Image
where Q = (A T PA) -1
PPT Slide
Lager Image
where N is number of LiDAR points.
4. Experimental Results and Analysis
- 4.1 Description of test data
Fig. 5 depicts how the simulated data was created for the experiments. Actually the test area was scanned with one strip. The data was intentionally divided into two strips with allowing overlapping between strips, and strip 1 serves as reference strip. Therefore, coordinates of the LiDAR points in strip 2 were transformed into strip 1 using the conjugate features in the overlapping area between adjacent strips. Table 1 shows the transformation parameters of each dataset.
Assigned transformation parameter values
PPT Slide
Lager Image
Assigned transformation parameter values
PPT Slide
Lager Image
Creation of test data strips
It is noted that the real datasets (Dataset B-1 and B-2) were formed with two different strips for E-W ( Fig. 7(a) ) and N-S ( Fig. 7(b) flight directions to analyze effect of the strip formation direction. Eventually the experiments aim to evaluate accuracy of the strip adjustment. Therefore, reference data is required to compute residuals of (X, Y, Z) coordinates between the reference and transformed data after strip adjustment. The original data of strip 2, ( i.e. , before applying transformation parameters given in Table 1 ) is regarded as reference data. Strip 2 datasets were created by Equation 17.
PPT Slide
Lager Image
PPT Slide
Lager Image
Real LiDAR data
where [ X’2, Y’2, X’2 ] T is strip 2 data after applying transformation using parameters listed in Table 1 . [ X2, Y2, X2 ] T is original strip 2 data ( i.e. , reference data) before transformation. Therefore, the transformation parameters and coordinates were computed through the adjustment procedure, and accuracies of the results were analyzed by comparing with the reference data.
Following LiDAR datasets were used for the experiments:
  • (1) Simulated data (Dataset A) includes buildings with various roof shapes with 0.25m average GSD (Ground Sampling Distance) after resampling and ±0.10m random noise (Fig. 6).
  • (2) Real data (Dataset B-1, B-2) of urban area provided by Saehan Aero Survey Co. Ltd. using Optech ALTM Gemini with 70Hz/sec scan frequency, 0.1m accuracy, and 0.6m average GSD after resampling (Fig. 7). Nearest neighbor interpolation was applied to both simulated and real datasets for resampling.
PPT Slide
Lager Image
Simulated LiDAR data
- 4.2 Results and analysis
Automatically extracted conjugate candidates in the overlapping strips are displayed in Fig. 8 , and conjugate feature locations after matching were derived as shown in Fig. 9 . Sizes of the template and search area are 13×13 pixels and 29×29 pixels, respectively in this experiment. The corner points were detected on the edges and the edges represent building boundaries that were extracted by Canny operator. In Fig. 9 , intersections of the cross lines indicate the sub-pixel locations of the conjugate points where the maximum correlation value occur in each search surface. It is obvious that the correlation values increase with subpixel accuracy as shown in Fig. 9 .
PPT Slide
Lager Image
Extraction of conjugate points
PPT Slide
Lager Image
Conjugate locations with sub-pixel accuracy
Table 2 shows adjustment parameters, discrepancies, dispersions of the parameters for each dataset. Discrepancies are the differences of adjustment parameters, and dispersions of the parameters (or RMSE of the parameters) indicate accuracy of the parameters. As for the simulated data “A”, the accuracy of the parameters are quite high based on the dispersion values. However, dispersion of ω (=0.057º) is relatively lager than other rotation parameters (φ=0.006º, and κ=0.005º) because the flying direction is E-W. In other words, rotation angle with respect to X-axis is obviously not stable for E-W strip direction. This fact is also found in the real dataset “B-1” and “B-2”. For dataset “B-1”, dispersion of ω (=0.030º) is larger than dispersions of other rotation parameters (φ=0.001º, and κ=0.003º), and for dataset “B-2” dispersion of φ (=0.136º) is larger than dispersions of other rotation parameters (ω =0.005º, and κ=0.006º). Graphical representation would be helpful for better analysis of the results. Fig. 10 shows the discrepancy and dispersion of the parameters (or RMSE of the parameters).
Results of adjustment parameters
PPT Slide
Lager Image
Results of adjustment parameters
PPT Slide
Lager Image
Discrepancy and RMSE of adjustment parameters
Therefore, it is conclusive that the accuracy of the rotation parameter depends upon the strip formation direction. In consequence, accuracy of the rotation parameters also affects accuracy of the planimetric translation parameters. For dataset “B-1” and “B-2”, dispersion of Y T and X T are 0.172m and 0.357m, respectively that are relatively larger than X T and Y T are 0.052m and 0.072m, respectively. Hence, accuracy of Y T and X T , are influenced by accuracy of ω and φ, respectively. In Table 2 , gray shaded cells indicate influence of the rotation parameters to the translation parameters. In the simulated data, the conjugate points are well distributed of even though the strip is elongated in E-W direction. Therefore, the ω does not affect to the planimetric translation parameters, i.e. , X T and Y T . Finally, there is no significant scale change between strips. Therefore, scale parameter could be ignored in the LiDAR strip adjustment.
Comparison of the RMSE of (X, Y, Z) coordinates before and after adjustment is important to evaluate quality of the adjustment procedure. RMSEs for all LiDAR points were computed as shown in Table 3 and Fig. 11 . Total LiDAR points are 85,272; 148,104; 216,270 for dataset A, B-1 and B-2, respectively. It is noticeable that there is significant accuracy improvement in Z-coordinate, especially. RMSE of each dataset was reduced by 0.504m, 2.461m, and 1.677m. This fact is quite meaningful through the LiDAR strip adjustment because discrepancy of Z-coordinate should be minimized in order for LiDAR data to be consistent and seamless. Fig. 12 shows parts of the profiles for before and after strip adjustment. It is visually verified that the discrepancy of two profiles were reduced after strip adjustment.
Table 3. RMSE of before and adjusted coordinates, and differences (unit: m)
PPT Slide
Lager Image
Table 3. RMSE of before and adjusted coordinates, and differences (unit: m)
PPT Slide
Lager Image
RMSE of adjusted coordinates
PPT Slide
Lager Image
Profi les of before and after adjustment
5. Conclusions
ALS system reaches a degree of the required accuracy in surveying, mapping and various applications. In order to produce high quality spatial products, calibration and strip adjustment of the LiDAR data are important tasks to create seamless and compatible point cloud data. In this study, strip adjustment was performed by 3D transformation between overlapping strips by utilizing conjugate features. The crucial task is automatic identification of the conjugate feature with accuracy requirement.
In order to increase accuracy and efficiency for determining conjugate features, it is recommendable to select well-defined points in both strips. Cross-correlation matching with sub-pixel accuracy was implemented on building corners that were extracted by corner detector. The proposed scheme was applied to the simulated and real LiDAR data. Following concluding remarks and future works were drawn based on the experimental results and analysis:
  • (1) In general, cross-correlation matching method is sensitive to rotation and scale, however, influence of the rotation and scale are not significant under the LiDAR data acquisition condition.
  • (2) Matching for determining conjugate points with subpixel reached up to half pixel accuracy.
  • (3) Overall accuracies of the adjustment parameters for each dataset were estimated as follows in terms of dispersions:
  • • Rotations = 0.06º, Planimetric shift = less than 0.01m, Vertical shift = less than 0.05m for simulated data (A)
  • • Rotations = 0.03º, Planimetric shift = less than 0.30m, Vertical shift = less than 0.40m for real data (B-1)
  • • Rotations = 0.14º, Planimetric shift = less than 0.40m, Vertical shift = less than 1.00m for real data (B-2)
  • (4) Accuracies of the strip coordinates after adjustment were improved as follows in terms of absolute RMSE difference:
  • • △RMSEX= 0.049m, △RMSEY= 0.421m, △RMSEZ= 0.504m for simulated data (A)
  • • △RMSEX= 0.521m, △RMSEY= 0.304m, △RMSEZ= 2.461m for real data (B-1)
  • • △RMSEX= 0.981m, △RMSEY= 0.741m, △RMSEZ= 1.677m for real data (B-2)
  • (5) Cross-correlation matching requires regular grid DEM, therefore, there might be errors due to resampling of the LiDAR data.
  • (6) Accuracies of the rotation parameters are sensitive to flight direction.
  • (7) It is recommendable to divide entire overlapping strips into several regions, and to perform strip adjustment independently in order to reflect local characteristics of the terrain features..
Acknowledgements
This research was supported a grant from geospatial information workforce development program funded by the Ministry of Land, Infrastructure and Transport of Korean Government (2014-05-03).
References
Bang K. , Habib A. , Müller M. 2009 LiDAR system calibration using overlapping strips Bol. Ciênc. Geod. 15 (5) 725 - 742
Canny J. 1986 A computational approach to edge detection IEEE Transactions on Pattern Analysis and Machine Intelligence 8 (6) 679 - 698
Crombaghs M. , De Min E. , Bruegelmann R. 2000 On the adjustment of overlapping strips of laser altimeter height data International Archives of Photogrammetry and Remote Sensing 33 (B3/1) 230 - 237
Filin S. , Vosselman G. 2004 Adjustment of airborne laser altimetry strips International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences 35 (B3) 285 - 289
Habib A. , Kersting A. , Ruifanga Z. , Al-Durgham M. , Kim C. , Lee D.C. 2008 Lidar strip adjustment using conjugate linear features in overlapping strips International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences XXXVII (B1) 385 - 390
Harris C. , Stephens M. 1988 A combined corner and edge detector Proceedings of Fourth Alvey Vision Conference Manchester, UK 31 August - 2 September, 1998 147 - 151
Kersting A. , Habib A. 2012 A comparative analysis between rigorous and approximate approaches for LiDAR system calibration Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography 30 (6-2) 593 - 605
Kilian J. , Haala N. , Englich M. 1996 Capture of elevation of airborne laser scanner data International Archives of Photogrammetry and Remote Sensing 31 (B3) 383 - 388
Lee D. , Yoo E. , Yom J. H. , Lee D. C. 2014 3D surface matching for LiDAR strip adjustment Proceedings of Korean Society of Surveying, Geodesy, Photogrammetry and Cartography Seoul, Korea 24-25 April (in Korean with English abstract) 147 - 149
Maas H. G. 2000 Least-squares matching with airborne laserscanning data in a TIN structure International Archives of Photogrammetry and Remote Sensing 33 (B3/1) 548 - 555
Pfeifer N. , Bӧhm J. , Li Z. , Chen J. , Baltsavias E. 2008 Advanced in Photogrammetry, Remote Sensing and Spatial Information Sciences: 2008 ISPRS Congress Book CRC Press/Balkema Leiden, The Netherlands Early stage of LiDAR data processing 169 - 184
Pfeifer N. , Oude Elberink S. , Filin S. 2005 Automatic tie elements detection for laser scanner strip adjustment International Archives of Photogrammetry and Remote Sensing 36 (3/W3) 1682 - 1750
Rentsch M. , Krzystek P. 2012 Lidar strip adjustment with automatically reconstructed roof shapes The Photogrammetric Record 27 (139) 272 - 292
Schenk T. 2001 Modeling and Analyzing Systematic Errors of Airborne Laser Scanners, Technical Notes in Photogrammetry, No. 19 Department of Civil and Environmental Engineering and Geodetic Science, The Ohio State University Columbus, OH. 40 -