Strip Adjustment of Airborne Laser Scanner Data Using Area-based Surface Matching

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.

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 et al., 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 et al., 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 et al. (2000) proposed a method for reducing vertical discrepancies between overlapping strips. Habib et al. (2008) utilized linear feature matching to solve corresponding problem in the overlapping areas. Kilian et al. (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.

Fig. 1.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 et al., 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.

Fig. 2.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.

The coordinates of a corner point in the local window are computed by:

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 et al., 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.

where n and m denote the template size (i.e., n×m), ZT(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, and 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.

Fig. 3.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.

Fig. 4.3D coordinate transformation between strips

where = [X2, Y2, Z2]T, = [XT, YT, ZT]T, = S·R·[X1, Y1, Z1]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.

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.

The linearized observation equation for least-square solution is written as Equation 11.

where [dXT, dYT, dZT, dS, dω, dφ, dκ]Ti are updating increments of the parameters at i-th iteration, and (eX, eY, eZ) is error vector. The matrix form of the observation equation is shown in Equation 12.

The transformation parameters by least-squares are computed by Equation 13:

where ε is permissible tolerance, and P is weight matrix. Finally, the iteration is to be terminated if () ≤ ε ,

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).

where n and m are number of conjugate points and adjustment parameters, respectively. Dispersion indicates reliability of the parameters.

where Q = (ATPA)-1

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.

Table 1.Assigned transformation parameter values

Fig. 5.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.

Fig. 7.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.

Fig. 6.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.

Fig. 8.Extraction of conjugate points

Fig. 9Conjugate 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).

Table 2.Results of adjustment parameters

Fig. 10.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 YT and XT are 0.172m and 0.357m, respectively that are relatively larger than XT and YT are 0.052m and 0.072m, respectively. Hence, accuracy of YT and XT, 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., XT and YT. 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)

Fig. 11.RMSE of adjusted coordinates

Fig. 12.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..