Accuracy Improvement of Airborne Lidar Strip Adjustment by Using Height Data and Surface Feature Strength Information Derived from the Tensor Voting Algorithm

Abstract

Light detection and ranging (Lidar) spatial coordinates, especially height data, and the intensity data of point clouds are often used for strip adjustment in airborne Lidar. However, inconsistency in the intensity data and then intensity gradient data because of the variations in the incidence and reflection angles in the scanning direction and sunlight incident in the same areas of different strips may cause problems in the Lidar strip adjustment process. Instead of the Lidar intensity, a new type of data, termed surface feature strength data derived by using the tensor voting method, were introduced into the strip adjustment process using the partial least squares method in this study. These data are consistent in the same regions of different strips, especially on the roofs of buildings. Our experimental results indicated a significant improvement in the accuracy of strip adjustment results when both height data and surface feature strength data were used.

1. Introduction

Currently, light detection and ranging (Lidar) systems are widely used for rapidly acquiring various types of three-dimensional (3D) point cloud data to reduce the working time effectively in many fields, especially in surveying, engineering, and 3D smart city construction. Digital terrain models have been constructed from Lidar point clouds by using interpolation and filtering methods [1,2]. Lidar data have been applied for extracting building outlines and subsequently reconstructing 3D buildings [3,4]. In [5], a neural network was employed to classify different plant species according to the Lidar intensity. Moreover, hyperspectral data and Lidar data have been combined for improving classification confidence [6,7]. Because Lidar data can be updated quickly, Lidar technology has been adopted for detecting and monitoring urban changes [8,9]. Moreover, by using a waveform airborne Lidar, one can classify urban and forestry areas and analyze opaque solids and structures [10,11,12,13].

In general, errors in Lidar data can be ascribed to factors such as distance measurement errors, orientation determination errors of the inertial measurement unit (IMU), Global Navigation Satellite System (GNSS) positioning errors, timing errors, mounting errors, and errors of the geoidal model (if using the orthometric height system) [14,15]. The accuracy of spatial positioning with Lidar may be affected by these errors. The performance enhancement achieved when using Lidar, especially the geometric performance enhancement, depends on the quality and accuracy of Lidar datasets. Therefore, Lidar data should be first calibrated or corrected. Individual instruments such as GNSS receivers and IMUs as well as their integration with the system should be calibrated. Often, a Lidar system can be effectively calibrated using the so-called calibration field [16,17,18,19,20]. Another calibration method involves the use of planar patches derived from the photogrammetric bundle adjustment as control surfaces instead of a calibrated field [21]. After calibration, the dm accuracy of the Lidar point coordinates can be usually achieved.

However, errors still may exist in Lidar data because of imperfect calibration. For reducing these errors, strip adjustment is frequently adopted to refine Lidar data. A few studies have introduced polynomial models or have matched flat and smooth planes in the strip adjustment process to eliminate displacements of the same objects in overlapping strips, especially height differences [22,23,24,25].

In addition to the use of height data in Lidar strip adjustment, other data, such as aerial image and intensity data, have been employed to acquire highly accurate Lidar data [26,27,28]. For instance, Burman [26] and Maas [27] have corrected height discrepancies by using the height and intensity data of Lidar points in strip adjustment. Zhang et al. [28] used aero-triangulated images to enhance Lidar strip adjustment by matching conjugated points and conjugated building corner features in intensity and aerial images. For correctly extracting and matching surface features, Filin [29] presented an error-recovery model for Lidar systems and established tie surfaces in overlapping strips. Lee et al. [30] corrected height discrepancies between overlapping strips with conjugated linear features by using the adjustment method. In summary, calibration and strip adjustment are necessary for acquiring highly accurate Lidar data.

Although Lidar intensity data may be applied for distinguishing objects according to the reflectivity of different materials and for strip adjustment, these data are easily affected by laser spreading loss, the incidence angles of sunlight and laser beams, scan angles, atmospheric attenuation, surface roughness, materials of objects, and other factors. Therefore, intensity data must be usually corrected in advance by following the so-called intensity normalization procedure. Even after subjecting Lidar intensity data to the intensity normalization procedure, guaranteeing their consistency on the same surfaces of different strips is not easy because of unpredictable factors such as unstable emission power, marginal variance in atmosphere conditions, and shadow regions (lack of texture) [31,32,33]. If the Intensity data are consistent, they are useful data for Lidar strip adjustment. Some researchers have successfully adopted such data for strip adjustment [26,27,28]. If they are inconsistent, it is difficult to use them in the strip adjustment.

Consistency in height and intensity between strips is crucial for obtaining high-quality Lidar strip adjustment results. The present study focuses on the Lidar strip adjustment by using consistent data, namely height data and the surface feature strength of objects. The surface feature strength data derived using the tensor voting method (TVM) provided useful geometric information [3,34,35].

In this study, the feasibility of using height data with surface feature strength data for conducting Lidar strip adjustment computation according to a method proposed in preliminary work studies [36,37] is further discussed and analyzed. The properties of surface feature strength data are introduced in Section 2. Given that a large number of parameters of the Lidar strip adjustment system are unknown, the partial least squares (PLS) method introduced in Section 3 is adopted instead of the ordinary least squares (OLS) method to improve computational effectiveness. The experiments conducted in this study are discussed in Section 4. The concluding remarks are presented in Section 5.

2. Geometric Feature Information

The surface feature strength of Lidar data is a type of geometric feature data derived from the TVM, and these data can be used to detect and extract geometric features of buildings [3]. You and Lin [35] used the TVM with topographic maps to reconstruct three-dimensional buildings from airborne Lidar data and achieved favorable results. We used the geometric feature extraction results obtained with the TVM, especially the surface feature strength data, for strip adjustment.

2.1. Geometric Feature Strength

You and Lin [3,35] demonstrated that the TVM algorithm is suitable for extracting geometric features from Lidar point clouds. By using the TVM, the geometric features of each Lidar point can be represented with a second-order tensor as follows [34]:

$$T=\left({\lambda }_1-{\lambda }_2\right)v_1{v}_1^T+\left({\lambda }_2-{\lambda }_3\right)(v_1{v}_1^T+v_2{v}_2^T)+{\lambda }_3(v_1{v}_1^T+v_2{v}_2^T+v_3{v}_3^T)$$

(1)

where λ₁, λ₂, and λ₃ are the eigenvalues of the tensor T and v₁, v₂, and v₃ are the corresponding eigenvectors. The geometric features of a Lidar point can be captured according to the following rules [34]: point features (${\lambda }_1\approx {\lambda }_2\approx {\lambda }_3)$ can be expressed in terms of λ₃, linear features (${\lambda }_1\approx {\lambda }_2\gg {\lambda }_3)$ can be expressed in terms of (λ₂ − λ₃), and surface feature (${\lambda }_1\gg {\lambda }_2\approx {\lambda }_3$) can be expressed in terms of (λ₁ − λ₂). Detailed decompositions of the tensors and the relationship between the tensors and the geometric features of a point can be obtained by referring to [3,34].

The geometric relationship between a desired point and its adjacent points within a specified region can be determined from the tensor T in Equation (1) when using the TVM algorithm. After TVM processing, the main geometric feature of a point can be determined using the eigenvalues λ₁, λ₂, and λ₃ of the tensor T according to the aforementioned capture rules. The geometric feature strengths may be affected by the number of adjacent points. To overcome this influence, You and Lin [3] suggested normalizing the values of the feature strengths as follows:

Surface feature strength:   C^s = (λ₁ − λ₂)/λ₁,

Line feature strength:   C^l = (λ₂ − λ₃)/λ₁,

Point feature strength:   C^p = λ₃/λ₁,

where $0\le C^s,C^l,C^p\le 1$ and $C^s+C^l+C^p=1.$ The surface feature strength ($C^s$) data and height data of Lidar points were used for strip adjustment in our study.

2.2. Properties of the Surface Feature Strength

In this section, we first review the important surface feature strength $C^s$ properties of Lidar points on a plane, as mentioned in [3]. The surface feature strengths of points on planes with the same slope are almost equal and approach approximately 1. By contrast, the strengths may vary at the junction of two adjacent planes of different slopes, for example, on roofs. The greater the slope difference between the two adjacent planes, the more significant is the $C^s$ variation, as depicted in Figure 1.

As an example, we illustrate the height and $C^s$ variation of a building consisting of several gable roofs by using real airborne Lidar data that were calibrated in advance by using the calibration field method. Figure 2a–c show the aerial photos, height data, and $C^s$ values, respectively, of the multi-gable roofs of the building. The height values (as z-data in this study) of the section line $\overline{PQ}$ (Figure 2a) are depicted in Figure 2c, and the $C^s$ values of $\overline{PQ}$ in different strips are illustrated in Figure 2f,g. Although the Lidar data were calibrated in advance, height discrepancies remained, as illustrated in Figure 2e.

The $C^s$ values are nearly the same on planes, for example, flat planes (e.g., the ground) or roof planes of buildings, inclined or flat roof planes, or ground planes. However, the $C^s$ values of the section line $\overline{PQ}$ vary at the roof ridges of adjacent roof planes with differential slopes (Figure 2c,d,f). The $C^s$ values in the neighborhood of the roof ridge are evidently different from those in the other parts of roof planes; however, their patterns in different strips of the same building are extremely similar, as illustrated in Figure 2f,g. The surface feature strength values of the same objects in different strips are consistent.

This consistency of $C^s$ values can be used for strip adjustment. Therefore, in addition to the height data of Lidar points, the $C^s$ values of Lidar points were used as complementary data when performing strip adjustment in this study. The flat roof planes or inclined roof planes and flat ground planes, for which the $C^s$ values vary marginally, are termed as the areas that exhibit slight variation in $C^s$ values (SVC areas) in this study. These planes are indicated in brown in Figure 2d. The neighborhood of the roof ridges of adjacent roof planes with different slopes are the areas with high variation in $C^s$ values (hereafter, HVC areas), as indicated by the blue parts in Figure 2d. The green regions in Figure 2d are the $C^s$ images of trees or the underbrush, and their $C^s$ values vary considerably. However, such messy distributed point clouds are unsuitable for strip adjustment and are excluded. In this study, $C^s$ data pertaining to the HVC areas were used in the strip adjustment to study the effects of surface feature strength data.

3. Strip Adjustment

Inaccuracy in Lidar data may exist even after the Lidar system has been calibrated. To correct such inaccuracy, strip adjustment is a powerful tool. Strip adjustment usually involves using the height discrepancies of the same object in the overlapping regions and occasionally involves using the Lidar intensity data if consistent intensity data are available. The three-, seven-, nine-, and 12-parameter transformation models as well as polynomial models, bilinear interpolation models using grids, linear interpolation models using a triangulated irregular network (TIN), corrections based on the generation function of Lidar data, and other methods have been adopted for Lidar strip adjustment [22,23,24,25,26,27,28].

Because each intensity-normalized method is associated with different restricted items [31,32,33], maintaining a consistent intensity for the same object in different strips is difficult even after normalization. Therefore, intensity data are excluded from strip adjustment in this study. Instead of intensity data, surface feature strength ($C^s$) values are used. The z-values (height values) and $C^s$ values of objects in the overlapping regions of different strips serve as the main data in our strip adjustment scheme. In addition, PLS is adopted for the strip adjustment because it offers advantages in terms of computational efficiency. The conjugated points for the adjustment are the grid points of the overlapping regions. An interpolation method is also required for adjustment. Details regarding the grids and bilinear interpolation method used in this study can be obtained from [26].

3.1. Mathematical Models and Grid Setting

Every Lidar strip is assigned to a seven-parameter similarity transformation model, and these parameters are used to absorb the displacement errors between overlapping strips. The model is expressed as follows:

$${\left(\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right)}_i^j={\kappa }^j·\mathrm{R}\left({\mathsf{\alpha }}^j,{\mathsf{\beta }}^j,{\mathsf{\gamma }}^j\right)·{\left(\begin{array}{c}x\\ y\\ z\end{array}\right)}_i^j+\left(\begin{array}{c}{t}_x^j\\ t_y^j\\ t_z^j\end{array}\right)$$

(2)

where the elements of the rotational matrix $\mathrm{R}\left({\mathsf{\alpha }}^j,{\mathsf{\beta }}^j,{\mathsf{\gamma }}^j\right)$ are as follows:

$$\begin{array}{l}\mathrm{R}\left(1,1\right)={\cos \mathsf{\beta }}^j·{\cos \mathsf{\gamma }}^j,\\ \mathrm{R}\left(1,2\right)={\cos \mathsf{\alpha }}^j·{\sin \mathsf{\gamma }}^j+{\sin \mathsf{\alpha }}^j·{\sin \mathsf{\beta }}^j·{\cos \mathsf{\gamma }}^j,\\ \mathrm{R}\left(1,3\right)={\sin \mathsf{\alpha }}^j·{\sin \mathsf{\gamma }}^j-{\cos \mathsf{\alpha }}^j·{\sin \mathsf{\beta }}^j·{\cos \mathsf{\gamma }}^j,\\ \mathrm{R}\left(2,1\right)=-{\cos \mathsf{\beta }}^j·{\sin \mathsf{\gamma }}^j,\\ \mathrm{R}\left(2,2\right)={\cos \mathsf{\alpha }}^j·{\cos \mathsf{\gamma }}^j-{\sin \mathsf{\alpha }}^j·{\sin \mathsf{\beta }}^j·{\sin \mathsf{\gamma }}^j,\\ \mathrm{R}\left(2,3\right)={\sin \mathsf{\alpha }}^j·{\cos \mathsf{\gamma }}^j+{\cos \mathsf{\alpha }}^j·{\sin \mathsf{\beta }}^j·{\sin \mathsf{\gamma }}^j,\\ \mathrm{R}\left(3,1\right)={\sin \mathsf{\beta }}^j,\\ \mathrm{R}\left(3,2\right)=-{\sin \mathsf{\alpha }}^j·{\cos \mathsf{\beta }}^j,\\ \mathrm{R}\left(3,3\right)={\cos \mathsf{\alpha }}^j·{\cos \mathsf{\beta }}^j.\end{array}$$

(3)

The superscript j represents the strip number; subscript i represents the Lidar point i; ${\mathsf{\alpha }}^j,{\mathsf{\beta }}^j,{\mathrm{and}\mathsf{\gamma }}^j$ represent the rotational angles of the jth strip, ${\kappa }^j$ is the scale factor, $t_x^j,t_y^j,\mathrm{and} t_z^j$ are the shift parameters; and ($x_i^j,y_i^j,z_i^j$) and (${x^{\prime }}_i^j,{y^{\prime }}_i^j,{z^{\prime }}_i^j$) are the coordinates of a laser point before and after transformation, respectively.

A series of square grids are first set in the flight regions to connect different strips. The (x, y) coordinates of every grid point are set to be known and are calculated as follows:

$$\begin{array}{ll}{X}_s=M_0+G·s,& \begin{array}{cc}& \begin{array}{ccc}& & s=0\cdots m,\end{array}\end{array}\\ Y_t=N_0+G·t,& \begin{array}{cc}& \begin{array}{ccc}& & t=0\cdots n,\end{array}\end{array}\end{array}$$

(4)

where ($M_0,N_0$) are the known coordinates of the original point of the continuous grid and G is the given length of a square grid. The height ($Z_{s,t})$ and the surface feature strength $\left(C_{s,t}\right)$ of every grid point are functions of the coordinates ($X_s,Y_t)$, that is, $Z_{s,t}$ = $f_Z$($X_s,Y_t)$ and $C_{s,t}$ = $f_C$($X_s,Y_t)$. The height and the surface feature strength data of each Lidar point can be derived from those of the grid points by using bilinear interpolation [26]:

$$\begin{array}{cc}\hfill {z^{\prime }}_i^j& =\frac{\left(X_{s+1,t+1}-{x^{\prime }}_i^j\right)\left(Y_{s+1,t+1}-{y^{\prime }}_i^j\right)}{G^2}{Z}_{s,t}+\frac{\left({x^{\prime }}_i^j-X_{s,t}\right)\left(Y_{s,t+1}-{y^{\prime }}_i^j\right)}{G^2}{Z}_{s+1,t}\hfill \\ & +\frac{\left(X_{s+1,t}-{x^{\prime }}_i^j\right)\left({y^{\prime }}_i^j-Y_{s,t}\right)}{G^2}{Z}_{s,t+1}+\frac{\left({x^{\prime }}_i^j-X_{s,t}\right)\left({y^{\prime }}_i^j-Y_{s,t}\right)}{G^2}{Z}_{s+1,t+1},\hfill \end{array}$$

(5)

$$\begin{array}{c}{C^{\prime }}_i^j=\frac{\left(X_{s+1,t+1}-{x^{\prime }}_i^j\right)\left(Y_{s+1,t+1}-{y^{\prime }}_i^j\right)}{G^2}{C}_{s,t}+\frac{\left({x^{\prime }}_i^j-X_{s,t}\right)\left(Y_{s,t+1}-{y^{\prime }}_i^j\right)}{G^2}{C}_{s+1,t}\\ +\frac{\left(X_{s+1,t}-{x^{\prime }}_i^j\right)\left({y^{\prime }}_i^j-Y_{s,t}\right)}{G^2}{C}_{s,t+1}+\frac{\left({x^{\prime }}_i^j-X_{s,t}\right)({y^{\prime }}_i^j-Y_{s,t})}{G^2}{C}_{s+1,t+1},\end{array}$$

(6)

By introducing Equation (2) into Equations (5) and (6) as well as the “total errors” v_z and v_C, the following complex observation equations are adopted for our strip adjustment scheme:

$$\begin{array}{cc}\hfill 0+v_Z=& z_i^j-{z^{\prime }}_i^j=F_Z{(}^j,{\mathsf{\alpha }}^j,{\mathsf{\beta }}^j,{\mathsf{\gamma }}^j,t_x^j,t_y^j,t_z^j,Z_{s,t},Z_{s+1,t},Z_{s,t+1},Z_{s+1,t+1},\hfill \\ & x_i^j,y_i^j,z_i^j,X_{s,t},Y_{s,t},X_{s+1,t},Y_{s+1,t},X_{s,t+1},Y_{s,t+1},X_{s+1,t+1},Y_{s+1,t+1}),\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill 0+v_C=& C_i^j-{C^{\prime }}_i^j=F_C{(}^j,{\mathsf{\alpha }}^j,{\mathsf{\beta }}^j,{\mathsf{\gamma }}^j,t_x^j,t_y^j,t_z^j,C_{s,t},C_{s+1,t},C_{s,t+1},C_{s+1,t+1},\hfill \\ & x_i^j,y_i^j,X_{s,t},Y_{s,t},X_{s+1,t},Y_{s+1,t},X_{s,t+1},Y_{s,t+1},X_{s+1,t+1},Y_{s+1,t+1}),\hfill \end{array}$$

(8)

where $x_i^j,y_i^j,z_i^j$ and $C_i^j$ denote the measurement values of the coordinates and the surface feature strength value at a Lidar point, respectively. The remaining variables in the functions F_z and F_C are unknown, including the transformation parameters of each strip and the heights and the surface feature strength values of grid points.

Compared with the results obtained using height observation data alone, it may be expected that the reliability of the results of strip adjustment would be higher when using both height and surface feature strength data. Notably, the parameter $t_z^j$ could not be solved when using $C^s$ observation data alone, which is similar to using only intensity data in the adjustment [26].

3.2. PLS Method

Because the number of Lidar points and the associated observation equations are large in number and Equations (5) and (6) are nonlinear, a linearization and iteration procedure is necessary for strip adjustment. Because of the demand for computational effectiveness, the PLS method is selected for strip adjustment. Helland [38,39] and Young [40] have suggested the use of the PLS method to reduce the dimensionality of unknown parameters for increasing operational efficiency. The analysis results of the OLS and PLS program source codes indicated that the time complexity and operational time of PLS are higher than those of OLS [41]. In addition, the PLS method may reduce the correlation among parameters, which may lead to a more reliable solution [42].

The PLS method also uses the least squares principle to solve the unknowns; however, it does not solve all unknown parameters simultaneously. Instead, it divides the unknown parameters into two parts. A variable transformation may be used for reducing the correlation between parameters [42]. In principle, at the beginning of the adjustment, the parameters in one of the parts (named the first part hereafter) are treated as knowns by assigning initial values to them and the parameters in the second part are regarded as unknowns that are determined by the least squares principle. Once the parameters are estimated, the parameters belonging to the second part are treated as knowns in the second adjustment step, whereas the parameters in the first part are now considered unknowns. In the third adjustment step, the parameters in the first part estimated in the previous adjustment step are treated as knowns and the parameters in second part are now considered unknowns again. The adjustment process continues with such a sequential exchange of knowns and unknowns until the estimation is convergent [43].

In this study, the unknown parameters are divided into two parts: the heights and $C^s$ values of the grid points are the parameters in the first part and the datum transformation parameters are the parameters in the second part. A flowchart of the strip adjustment procedure performed in this study with the PLS method is illustrated in Figure 3. Before strip adjustment by using the PLS method, the heights and surface feature strength of the grids are assigned initial values according to the closest points. The initial value of the scale parameter is set to 1 and those of the other datum parameters are set to 0. At the beginning of strip adjustment, the parameters in the first part, namely the initial heights and $C^s$ values of the aforementioned grid points, are considered known information and the parameters in the second part, namely the transformation parameters of individual Lidar strips, are estimated with the OLS method. After each step of strip adjustment, new coordinates of laser points are updated according to the new estimated transformation parameters. The final solutions were obtained when the PLS procedure ceased until the changes in the unknown parameters (including grid values and transformation parameters) became lower than the corresponding thresholds.

If Lidar data (including the z-data and $C^s$ data) contain outliers, the results of strip adjustment may deteriorate. Therefore, an outlier detection procedure should be considered. In this study, the outliers in Lidar data were detected using the robust method [44]. Moreover, observational data of different types, namely height data and $C^s$ data in this case, should be assigned correct weight values in the adjustment process. When the observational data lack correct priori information about weights, variance component estimation should be performed for the observation data of different types. We employed Helmert variance component estimation to obtain accurate variances of the observation data for achieving the optimal estimation of unknown parameters [45].

4. Experiments and Discussion

The study region covered with four Lidar strips is located on the campus of National Cheng-Kung University in Tainan City, Taiwan (Figure 4). The area of the study region is approximately 360,000 m². The width of flights A and B was approximately 350 m, and the length of flights A and B was 1,070 m. The width of flights C and D was approximately 250 m, and the length of flights A and B was 570 m. The flight height was approximately 500 m above ground level. Lidar point clouds were captured using an Optech ALTM 30/70 scanner, and the Lidar density was found to be approximately 6 pts/m². Although the Lidar point clouds in the study region were calibrated in advance, discrepancies were readily apparent, as depicted in Figure 2e. Therefore, the strip adjustment and grid interpolation methods mentioned in Section 3 were used to correct such discrepancies.

First, the study region was set into grids with a width of 1 m in the experiment. The regions O1, O3, O4, O5, and O7 in Figure 4 were overlapped by two strips, and the regions O2 and O6 were overlapped by three strips. Second, the observational equations (Equations (7) and (8)) of all Lidar points, excluding those of trees, the underbrush, the ground, and some selected check planes in the overlapping areas, were adjusted using the PLS method. A few inclined or flat roof planes and ground planes were selected as the check planes. Observations of Lidar points on these planes were not used in the adjustment. After the adjustment was performed, the adjusted height and surface feature strength of each Lidar point were calculated using the transformation parameters and the interpolation method expressed in Equations (2), (5), and (6).

In this study, the accuracy of the adjustment results was evaluated using two quality indicators, namely relative deviation (RED) and absolute deviation (ABD). The RED was calculated using all the residuals of observations in the check planes as follows:

$$\mathrm{RED}=\sqrt{\frac{{\displaystyle \sum }\left(V_Z{}^T{V}_Z\right)}{n}},$$

(9)

where n is the number of Lidar points. The RED can be used to reveal the quality of matching between the strips after adjustment are made. The residual $v_z$ in Cases 2–4 was equal to the z-observational value minus the adjusted z value, whereas in Case 1, the residual $v_z$ was calculated from the differences in the calibrated z-data of the same object in adjacent strips. Cases 1–4 are defined in the following paragraph.

Measurements of 20 flat check planes were conducted using three to five points through the real-time kinematic GPS (RTK-GPS) positioning by using a Leica SR530 GPS receiver. Such check planes were used for external examination of the adjustment results (Figure 4). The ABD indicator can be expressed as follows:

$$\mathrm{ABD}=\sqrt{\frac{{\displaystyle \sum }(\left(d_i^w{)}^T·d_i^w\right)}{n}},$$

(10)

where n represents the sum of the numbers of check points on all such check planes, $d_i^w={z^{\prime }}_i^w-{{\displaystyle \sum }}_{k=1}^m{z}_{k,GPS }^w$/m (m = 3~5) represents the height residual of the ith Lidar point on the wth flat plane, and $z_{k, GPS}^w$ represents the kth GPS measurement height on the wth flat plane. The GPS heights and the heights of Lidar points were obtained by referring to the same height system in this study.

To analyses the effects of the surface feature strength on our strip adjustment scheme, we designed the following four cases to evaluate the quality of the adjustment results with or without the surface feature strength values. In Case 1, the initial height data (z-data) after calibration were used without the subsequent PLS strip adjustment. In Case 2, the z-data of the Lidar points in the HVC areas obtained after strip adjustment were used. In Case 3, the z-data of the Lidar points in the SVC areas only, as used commonly in airborne Lidar adjustment, were used. In Case 4, the z-data of the Lidar points in the SVC areas and the $C^s$ data in the HVC areas were used. In the aforementioned cases, the $C^s$ data and z-data of the same locations were not used simultaneously in the adjustment. The z-data used in Cases 2–4 were the same as those used in Case 1.

Because the $C^s$ data in the SVC areas, such as in the plane roofs, were nearly identical, some unknown parameters could not be solved in our experiments when only the $C^s$ data were used in the SVC areas. By contrast, z-data are usually useful in the SVC areas. Therefore, the use of $C^s$ data for strip adjustment in the SVC areas is not discussed in this paper.

The adjustment results are summarized in

Table 1. Summation of RED and ABD.

Case	RED (m)	ABD (m)
1	0.203	0.111
2	0.406	0.261
3	0.067	0.088
4	0.057	0.073

. The experimental results indicate that in Case 1, in which only the calibrated z-data were used, and without using strip adjustment, the RED accuracy was 20 cm and the ABD accuracy was 11 cm. Figure 5 illustrates the discrepancies between the strips in section lines $\overline{GH}, \overline{KL}, \overline{PQ}$, $\overline{RS}, \overline{TU},\mathrm{and} \overline{VW}$ in Cases 1, 2, 3, and 4. The abscissa axis indicates the direction along the section line and the vertical axis represents the elevation in meters. These discrepancies negatively affect the working of highly accurate applications using Lidar data. The adjustment results of Case 2 were worse than those of Case 1. This indicates that conducting Lidar strip adjustment is not feasible with the use of the z-data in the HVC areas. By using z-data in the SVC areas (Case 3), one may obtain better results, as summarized in Table 1 and illustrated in Figure 5. The same conclusion was arrived at in [26].

According to statistical hypothesis testing at the 95% significance level, the adjustment results obtained in Cases 3 and 4 were significantly superior to those obtained in Case 1. Compared with Case 1, the RED and ABD results improved by 67 and 21% in Case 3, respectively, and by 72 and 34% in Case 4, respectively. The results of Case 4 were superior to and more accurate than those of Case 3. In Case 4, the RED and ABD improved significantly by approximately 15 and 17%, respectively, compared with the corresponding values in Case 3. The quality of the strip adjustment results obtained using both z-data and $C^s$ data was superior to that obtained using z-data alone. In general, the experiments indicate that the surface feature strength ($C^s$) is useful for improving the results of strip adjustment significantly.

To analyses the levels of improvement in different areas when using $C^s$ data, the adjustment results obtained in regions O1–O7 are discussed in the following text. The related RED and ABD values are illustrated in Figure 6. The accuracy of the results obtained in regions O2, O4, and O6 in Case 1 was higher than that of the results obtained in the other regions. The errors may have resulted from the large discrepancies between strips A and B despite the fact that the Lidar source data used in the process were calibrated in advance. After strip adjustment by using Case 3 data, the RED of all regions except that of O5 reached the cm-level. As illustrated in Figure 4, this result was expected because fewer suitable z-data points were available from region O5 than from other regions for use in the adjustment. Moreover, the data available from region O5 had poor distribution. However, the accuracy was improved when Case 4 data were used (Figure 6). The use of strip adjustment to correct discrepancies between strips, especially discrepancies in height, may not always be possible even when the z-data are distributed uniformly and over a sufficient SVC terrain. Strip adjustment may yield poor results if z-data alone are used.

However, the use of the z-data of Lidar points in the SVC areas and the $C^s$ data of the points in the HVC areas can improve the quality and accuracy of strip adjustment, as illustrated in Figure 5 and Figure 6. Notably, the z-data of SVC areas and the $C^s$ data of HVC areas can be used as complementary data for Lidar strip adjustment.

5. Conclusions

In this study, a new Lidar data type, namely the surface feature strength derived from the TVM, was introduced into Lidar strip adjustment. Compared with the use of height data alone for up-to-date airborne Lidar strip adjustment, the use of the surface feature strength data of Lidar points in the HVC areas as auxiliary data can improve the strip adjustment accuracy by 15–17% according to our experimental results. The new data, namely the surface feature strength, are useful auxiliary data for strip adjustment, especially when elevation data in SVC areas are insufficient.