An Automated Rectification Method for Unmanned Aerial Vehicle LiDAR Point Cloud Data Based on Laser Intensity

Abstract

Point cloud rectification is an efficient approach to improve the quality of laser point cloud data. Conventional rectification methods mostly relied on ground control points (GCPs), typical artificial ground objects, and raw measurements of the laser scanner which impede automation and adaptability in practice. This paper proposed an automated rectification method for the point cloud data that are acquired by an unmanned aerial vehicle LiDAR system based on laser intensity, with the goal to reduce the dependency of ancillary data and improve the automated level of the rectification process. First, laser intensity images were produced by interpolating the intensity data of all the LiDAR scanning strips. Second, a scale-invariant feature transform algorithm was conducted to extract two dimensional (2D) tie points from the intensity images; the pseudo tie points were removed by using a random sample consensus algorithm. Next, all the 2D tie points were transformed to three dimensional (3D) point cloud to derive 3D tie point sets. After that, the observation error equations were created with the condition of coplanar constraints. Finally, a nonlinear least square algorithm was applied to solve the boresight angular error parameters, which were subsequently used to correct the laser point cloud data. A case study in Shehezi, Xinjiang, China was implemented with our proposed method and the results indicate that our method is efficient to estimate the boresight angular error between the laser scanner and inertial measurement unit. After applying the results of the boresight angular error solution to rectify the laser point cloud, the planar root mean square error (RMSE) is 5.7 cm and decreased by 1.1 cm in average; the elevation RMSE is 1.4 cm and decreased by 0.8 cm in average. Comparing with the stepwise geometric method, our proposed method achieved similar horizontal accuracy and outperformed it in vertical accuracy of registration.

1. Introduction

Unmanned aerial vehicle (UAV) Light Detection and Ranging (LiDAR) is a new technology in the field of survey and mapping that is equipped with low-altitude UAV platform for LiDAR data acquisition and composed of three core components, including laser measurement, differential Global Navigation Satellite System (GNSS), and inertial navigation unit (IMU) [1,2]. Comparing with conventional aerial photogrammetry techniques, UAV LiDAR bears a number of advantages such as being less impacted by flying conditions (e.g., cloud cover and flexible ground control), high-level automation, more precision and density data, and high flexibility. Thus, it has been widely used in acquiring digital elevation model [3,4,5], disaster monitoring [6], heritage protection [7], forestry survey [8,9,10,11], and 3D modeling [12]. In these applications, attention was often given to the surveying accuracy of height or altitude, which is important for most situations. Among those factors that degrade the surveying accuracy of UAV LiDAR data [13,14,15], boresight angular error is a systematic error and has significant influence on the geolocation of laser points [15,16]. It usually happens during the payload installation of a UAV LiDAR system and the unstable UAV flying process, and consequently the boresight angular error cannot be considered to be negligible and is hard to be directly measured [17,18]. Moreover, it can lead to systematic positional deviations in all LiDAR scanning strips and misalignment of the same objects in overlapping areas in different LiDAR scanning strips [19,20]. Therefore, it is necessary to estimate and reduce down the errors in LiDAR point data processing. Obviously the accuracy of laser points has significant impact on the quality of mapping products, especially for applications such as matching and fusion of multisource remote sensing data [15,21]. Thus, it is currently a key step to explore reliable, intelligent, and efficient methods for error rectification to improve the quality of laser point cloud data that are captured by UAV LiDAR systems [2].

The mainstream idea for the error rectification of point cloud data captured by a LiDAR system is based on the LiDAR georeferencing equations, in which, all the possible error sources of a LiDAR system are considered to build and resolve error equations [18,19,22,23]. The core procedure is the construction of the observation error equation that is usually achieved by tie points, and the error parameters can be solved by a matrix operation. Zhang and Forsberg [24] proposed a simple boresight angular error rectification method, termed stepwise geometric correction. In this method the boresight angular parameters was examined by analyzing the relationship between the positioning displacement and the selected regular shape of the ground objects such as horizontal plane surface, which was used to build a model for estimating the geometric errors. Considering a calibration field or sufficient ground control points are necessary in the method; consequently the UAV flight lines need to be specially designed in the implementation. However, due to the fact that the modeling is still not accurate enough, the method is only applicable for rough estimation [24]. To overcome the difficulty in identifying tie points from overlapping LiDAR scanning strips, the complanate features were extracted interactively and the coplanar condition was used as a constraint to solve the boresight angular parameters [18]. However this method showed a low-level automation and heavily depended on typical artificial objects in the scanned area. In addition to the flight path data of the UAV, it also requires raw data of the laser scanner which is not accessible for most users in most cases because most producers of the UAV LiDAR systems do not release the raw data format in the software package. Considering the complexity of the rectification model and systematic error sources, a new boresight angular error rectification method that is not affected by GNSS and IMU observation errors was proposed by Le Scouarnec et al. [25]. This method requires scanning horizontal and vertical planes as reference objects, and the laser scanner must be kept static in the scanning process. The original observation information of the LiDAR scanner is a must when solving the boresight angular error, and as a result, it is more suitable for terrestrial laser scanning system. Due to the fact that the end user has difficulty in accessing raw observations of the laser scanner after a flight, a rigorous rectification model that does not need the raw observations of the laser scanning system has been proposed by Bang et al. [19], which considers almost all error sources besides boresight angular error, but still needs many ground control points and ground features of objects to resolve the model. In addition to the above-mentioned methods, Zhang et al. [15] noticed that previous studies paid less attention to the rectification of relatively low accurate position and orientation system (POS) data, and proposed an aerotriangulation-aided adjustment rectification model for LiDAR scanning strips to eliminate positioning and angular errors caused by boresight angular errors and POS errors. In the model proposed by Zhang et al. [15], the error rectification model was established by combining time-independent boresight angular error and time-dependent POS data error in order to correct the LiDAR strip data. This method requires sufficient accurate GCPs to ensure accurate aerial triangulation, and consequently the model is complicated with relatively low automation. In order to improve the automation of rectification of boresight angular error, a novel rectification model was proposed using the coplanar constraints, in which an automated approach for retrieving the rooftop facet and building walls was designed based on a region growth and a Random Sample Consensus (RANSAC) segmentation algorithm, and the tie planes were extracted from different LiDAR scanning strips and used as matching elements [16]. However, this method was limited to urban areas where dense rooftops exist, and the LiDAR data must be acquired by an oblique forward-looking full-waveform laser scanner to guarantee the acquisition of laser point cloud data of the building rooftop and building facades [16].

Generally speaking, the existing methods for boresight angular error rectification in previous studies mostly depend on precise ground control points or manually selected features of objects, and consequently they were implemented in a low degree of automation. Some researchers tried to create automated rectification models independent of ground control points, but to solve the models requires GNSS trajectory data as well as raw observations of the laser scanner(e.g., laser ranges and scan mirror angles), which were not often available for the end users in many cases. Thus, some studies added some restrictions to the experimental conditions to simplify the solving of the models, but the generality of the methods was reduced down. Most importantly, previous studies paid little attention to the laser intensity in building rectification models for boresight angular error, and the values of the laser point cloud have not been fully explored and utilized. In addition, a low-altitude UAV LiDAR system is prone to the boresight angular error during the pre-fly and flying stages because of its much lower installation precision and stability than a manned airborne platform. Removal of boresight angular error can ameliorate the quality of LiDAR point cloud data, which is the key step to realize accurate matching of different laser scanning strips. It is significant to develop a method for boresight angular error rectification with high degree of automation independent of ground control points, feature objects, and raw observations.

As mentioned above, the geolocation error induced by boresight angle is a main systematic error in a low-altitude UAV LiDAR system compared to other ones, and has a non-negligible impact on the laser point cloud data. Thus, this study is focused on the boresight angular error rectification problem, and presents an automated boresight angular error rectification method for a low-altitude UAV LiDAR system based on the laser intensity information. A case study in the Shihezi area, Xinjiang, China has demonstrated that our proposed method could reduce down the geolocation error of the LiDAR point cloud data without support of any ground control points, feature objects and raw observations of the scanner. In Section 2 we talk about data acquisition and the method; the results and analysis will be presented in Section 3. In Section 4 we discuss the influence of some factors on the accuracy of error rectification, and finally we conclude this research in Section 5.

2. Data and Methods

2.1. Study Area and Data Acquisition

The study area is located in a farming area in the northwest of Shihezi city, Xinjiang, China (44°24′06″N, 85°53′41″E). The UAV LiDAR payload contains the OXTS xNAV550 GNSS/IMU dual-GNSS antenna navigation system and the RIEGL VUX-1UAV laser scanner, which were integrated on a Swiss-made ScoutB1-100 unmanned helicopter [26]. Detailed information of the payload is listed in

Table 1. The core specifications of the UAV Light and Detection Ranging (LiDAR) system.

Laser Scanner1	Specifications	GNSS/IMU2	Specifications
Minimum Range	5 m	Positioning Mode	RTK
Pulse Repetition Rate	550 KHz	Data Frequency	100 Hz
Measurement Accuracy	0.015 m	Position Accuracy(CEP)	H:0.02 m; V:0.03 m
Scanning Speed	200 scan/s	Speed Accuracy	0.1 km/h
Angle Resolution	0.001°	Roll Accuracy (1σ)	0.05°
Field of View	330°	Pitch Accuracy (1σ)	0.05°
Echo Signal Intensity	16 bit	Heading Accuracy (1σ)	0.10°

¹http://www.riegl.com/products/newriegl-vux-1-series/newriegl-vux-1lr/. ²https://www.oxts.com/products/xnav/.

. Time synchronization is based on GNSS provided by the GNSS/IMU unit, and the UAV LiDAR system flew above the experimental area on 29 July, 2017 with the proper flight configuration of planned flight lines (Figure 1). The laser scanner was set to a field angle of 110° with a scanning rate at 550 kHz. The UAV flew at a low altitude above ground level of 30 m and a cruising speed of 5 m/s, and captured two flight strips with opposite directions in the experiment. The point density for each scanning line is ~400–600 pts/m² (Figure 2), the ground scanning width is ~85.7 m, and the spacing between the two flight trajectories is ~10 m. Each laser point attribute includes X, Y, and Z coordinates, the number of echoes at that point, scanning time, and laser intensity. In addition, the POS data of the strips can be derived from the inertial navigation system, including time, latitude, longitude, elevation, roll angle, pitch angle, and heading angle. These two data sets are the basic data sources for the following boresight angular error rectification.

2.2. Rectification of UAV LiDAR System Errors

2.2.1. LiDAR Georeferencing Equations

The geolocation error sources of a UAV LiDAR system usually include laser ranging error, GNSS positioning error, orientation angular error, lever arm error (displacements between laser scanner, IMU, and GNSS antenna), boresight angular error, and others. Among these errors, the boresight angular error and lever arm error belong to systematic errors, and the lever arm error can be measured by surveying instrument and usually has little effect on the accuracy of LiDAR point cloud data. However, the boresight angular error is usually caused by artificial factors or instable platform during the UAV LiDAR system installment stage prior to the flight. Although boresight angular error is not easy to be measured by surveying instruments, it can be estimated through a rigorous mathematic model. Based on the previous studies [20,22,27,28], the LiDAR georeferencing equation can be expressed as

$$\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]=R_W{R}_N\left[R_I{R}_L\left[\begin{array}{c}0\\ 0\\ \rho \end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_G\\ Y_W\\ Z_G\end{array}\right]$$

(1)

where [X_W Y_W Z_W]^T refers to the mapping frame coordinate of laser points; ρ denotes the distance measured by the laser scanner; [dx dy dz]^T refers to the direction vector from the laser scanner center to GNSS antenna phase center, which can be measured by surveying instruments or obtained from the configuration file; [X_W Y_G Z_G]^T refers to the mapping frame coordinate of GNSS antenna phase center; R_L is a 3*3 rotation matrix that transforms the instantaneous scanning coordinate system into the scanning reference coordinate system, i.e., the current direction of laser pulse emission; R_I is a 3*3 rotation matrix that transforms the scanner reference coordinate system into the IMU coordinate system; R_N is a 3*3 rotation matrix that transforms the IMU coordinate system into the navigation coordinate system; and R_W is a 3*3 rotation matrix that transforms the navigation system into the mapping frame coordinate system. Let

$$\left[\begin{array}{c}{X}_I\\ Y_I\\ Z_I\end{array}\right]=R_I{R}_L\left[\begin{array}{c}0\\ 0\\ \rho \end{array}\right]$$

(2)

Then Equation (1) can be rewritten as

$$\left[\begin{array}{c}{X}_I\\ Y_I\\ Z_I\end{array}\right]=R_N^{-1}{R}_W^{-1}\left(\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]-\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]\right)-\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]$$

2.2.2. Boresight Alignment Model

The principle of boresight angular error rectification is based on the LiDAR georeferencing model and the laser scanning reference coordinate system with axes deviation error can be rectified by rotating to the real laser scanning reference coordinate system with three boresight angular corrections. If the boresight angular error parameters of the LiDAR system in three directions of rolling, pitching, and heading are ω, ϕ and κ, respectively, the rotation matrix derived from the boresight angular error in the three directions can be described as [29]

$$R_M=\left[\begin{array}{ccc}\cos \kappa & -\sin \kappa & 0\\ \sin \kappa & \cos \kappa & 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}\cos \phi & 0& \sin \phi \\ 0& 1& 0\\ -\sin \phi & 0& \cos \phi \end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \omega & -\sin \omega \\ 0& \sin \omega & \cos \omega \end{array}\right]$$

(3)

where, R_M is the rotation matrix of the boresight angular error parameters. Considering that the boresight angular error parameter is usually small in number, Equation (3) can be approximated:

$$R_M\approx \left[\begin{array}{ccc}1& -\kappa & \phi \\ \kappa & 1& -\omega \\ -\phi & \omega & 1\end{array}\right]$$

(4)

If we only consider the boresight angular error and ignore the other error sources, the LiDAR georeferencing equation can be described as

$$\left[\begin{array}{c}{X}_W\\ Y_W\\ Z_W\end{array}\right]=R_W{R}_N\left[R_M\left[\begin{array}{c}{X}_I\\ Y_I\\ Z_I\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]$$

(5)

where [X_W Y_W Z_W]^T denotes the coordinate of laser points in the mapping frame and R_W and R_N refer to the matrix that transforms navigation system into the mapping frame coordinate system and the matrix that transforms the IMU coordinate system into the navigation coordinate system, respectively. [X_I Y_I Z_I]^T denotes the coordinate of laser points in IMU reference frame and [X_G Y_G Z_G]^T denotes the coordinate of GNSS antenna phase center in the mapping frame. [dx dy dz]^T refers to the direction vector from the laser scanner center to GNSS antenna phase center.

When P₁ and P₂ are two tie points, they should satisfy the following conditions.

$$\begin{array}{l}{P}_1:\left[\begin{array}{c}{X}_{W1}\\ Y_{W1}\\ Z_{W1}\end{array}\right]=R_W{R}_{N1}\left[R_M\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G1}\\ Y_{G1}\\ Z_{G1}\end{array}\right]\\ P_2:\left[\begin{array}{c}{X}_{W2}\\ Y_{W2}\\ Z_{W2}\end{array}\right]=R_W{R}_{N2}\left[R_M\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G2}\\ Y_{G2}\\ Z_{G2}\end{array}\right]\end{array}\}$$

(6)

where, [X_W1 Y_W1 Z_W1]^T and [X_W2 Y_W2 Z_W2]^T are the coordinate of tie points P₁ and P₂ in the mapping frame, respectively; [X_I1 Y_I1 Z_I1]^T and [X_I2 Y_I2 Z_I2]^T are the coordinate of tie points P₁ and P₂ in IMU reference frame, respectively; [X_G1 Y_G1 Z_G1]^T and [X_G2 Y_G2 Z_G2]^T are the coordinate of GNSS antenna phase center of tie points P₁ and P₂ in the mapping frame, respectively.

After applying the rectification model, conceptually, P₁ and P₂ would be the same point, which means

$$R_W{R}_{N1}\left[R_M\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G1}\\ Y_{G1}\\ Z_{G1}\end{array}\right]=R_W{R}_{N2}\left[R_M\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G2}\\ Y_{G2}\\ Z_{G2}\end{array}\right]=\left[\begin{array}{c}{X}_{true}\\ Y_{true}\\ Z_{true}\end{array}\right]$$

(7)

where R_N1 and R_N2 are matrixes that transform from the IMU coordinate system to the navigation coordinate system for tie points P₁ and P₂, respectively, and R_M is the rotation matrix of the boresight angular error parameters. [X_true Y_true Z_true]^T denotes the actual coordinate of tie points P₁ and P₂. Let

$$\begin{array}{l}{F}_1=R_W{R}_{N1}\left[R_M\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G1}\\ Y_{G1}\\ Z_{G1}\end{array}\right]\\ F_2=R_W{R}_{N2}\left[R_M\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]+\left[\begin{array}{c}dx\\ dy\\ dz\end{array}\right]\right]+\left[\begin{array}{c}{X}_{G2}\\ Y_{G2}\\ Z_{G2}\end{array}\right]\end{array}$$

Then

$$F=F_1-F_2=0$$

(8)

In Equation (8), F is a nonlinear function and can be linearized by using the Taylor series expansion. F₁ and F₂ are the rectified coordinates of P₁ and P₂ in the mapping frame, respectively. If only the first-order term and the constant term remain, Equation (8) can be rewritten as

$$F\approx F_0+\frac{\partial F}{\partial \omega }\Delta \omega +\frac{\partial F}{\partial \phi }\Delta \phi +\frac{\partial F}{\partial \kappa }\Delta \kappa \approx 0$$

(9)

where, F₀ represents the constant term that is an approximate value estimated using Equation (8) when R_M is set with the initial value of boresight angular error parameters. Δω, Δϕ, and Δκ are the first-order terms.

Therefore, the observation error equation can be expressed as

$$V=\frac{\partial F}{\partial \omega }\Delta \omega +\frac{\partial F}{\partial \phi }\Delta \phi +\frac{\partial F}{\partial \kappa }\Delta \kappa -(-F_0)$$

(10)

where, V denotes the residual error matrix. If there are n pairs of tie points, then the observation error equation is

$$\underset{3n\times 1}{{\displaystyle V}}=\underset{3n\times 3}{{\displaystyle A}}\underset{3\times 1}{{\displaystyle X}}-\underset{3n\times 1}{{\displaystyle L}}$$

(11)

Equation (11) is given with the following notions.

$$\begin{array}{l}A=\left[\begin{array}{ccc}\frac{\partial F}{\partial \omega }& \frac{\partial F}{\partial \phi }& \frac{\partial F}{\partial \kappa }\end{array}\right]\\ X={\left[\begin{array}{ccc}\Delta \omega & \Delta \phi & \Delta \kappa \end{array}\right]}^{\mathrm{T}}\end{array}$$

$$\begin{array}{c}\frac{\partial F}{\partial \omega }=\frac{\partial F_1}{\partial \omega }-\frac{\partial F_2}{\partial \omega }=R_W{R}_{N1}\frac{\partial R_M}{\partial \omega }\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]-R_W{R}_{N2}\frac{\partial R_M}{\partial \omega }\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]\\ \frac{\partial F}{\partial \phi }=\frac{\partial F_1}{\partial \phi }-\frac{\partial F_2}{\partial \phi }=R_W{R}_{N1}\frac{\partial R_M}{\partial \phi }\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]-R_W{R}_{N2}\frac{\partial R_M}{\partial \phi }\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]\\ \frac{\partial F}{\partial \kappa }=\frac{\partial F_1}{\partial \kappa }-\frac{\partial F_2}{\partial \kappa }=R_W{R}_{N1}\frac{\partial R_M}{\partial \kappa }\left[\begin{array}{c}{X}_{I1}\\ Y_{I1}\\ Z_{I1}\end{array}\right]-R_W{R}_{N2}\frac{\partial R_M}{\partial \kappa }\left[\begin{array}{c}{X}_{I2}\\ Y_{I2}\\ Z_{I2}\end{array}\right]\end{array}\}$$

Applying the least squares algorithm proposed by Marchant et al. [30], the solution of the Equation (11) is

$$X={(A{A}^{\mathrm{T}})}^{-1}AL$$

(12)

Thus we can get the approximate solution of three unknown parameters: ∆ω, ∆φ, and ∆κ. Since only the first-order term of the Taylor expansion in Equation (9) is considered, we can solve the observation error equations in an iterative way, where the coefficient values and constant terms are modified successively till the observation error converges to a presetting threshold. Finally, the value of each unknown boresight angular error parameter can be obtained as follows

$$\begin{array}{l}\omega ={\omega }_0+\Delta {\omega }_1+\Delta {\omega }_2+\cdots \\ \phi ={\phi }_0+\Delta {\phi }_1+\Delta {\phi }_2+\cdots \\ \kappa ={\kappa }_0+\Delta {\kappa }_1+\Delta {\kappa }_2+\cdots \end{array}$$

where ω, ϕ, and κ are the solution of the boresight angular error parameters, and ω₀, ϕ₀, and κ₀ are the initial values of the boresight angular error parameters. Δω_i, Δϕ_i, and Δκ_i represent the increment of the boresight angular error parameters at each iteration.

2.3. Automated Rectification Based on the Laser Intensity

When a UAV LiDAR system is used to scan terrain or targets over the scanned area, it can acquire the geospatial information and also record the reflection intensity of the scanned terrain and targets. The Rigel VUX-1LR laser scanner emits laser light pulses at the near-infrared wavelength centered at 1550nm. Due to the differences in the reflectivity of near-infrared laser light among different ground targets, the laser intensity may be helpful for feature points extraction and matching of the LiDAR flight strips. The Scale-Invariant Feature Transform (SIFT) algorithm [31] is more robust and ascendant in feature extraction and matching, but is only applicable for 2D image data. Therefore, we proposed a new approach, and the technical steps can be described as follows. First, the intensity of the point cloud data are interpolated to produce intensity images and the tie points are then retrieved based on the 2D adjacent intensity images, and based on the tie points in 2D intensity images, the tie points in 3D point cloud data can be also determined using a 2D-to-3D mapping strategy. After that, an observation error equation will be built based on the LiDAR georeferencing equation. The boresight angular error can be finally estimated by resolving the error equation using the extracted 3D tie point sets.

Because the distribution of laser points is spatially irregular and discrete, it is difficult to guarantee that the same ground point can be scanned in different scanning strips. Consequently, the tie points extracted from different LiDAR strips using SIFT algorithm are not real laser footprint points. This situation makes it difficult to obtain the corresponding observation information such as scanning time and POS data and, as a result, it becomes difficult or impossible to construct the LiDAR georeferencing equation. Considering there is a spatial constrained relationship between the tie points and their surrounded laser points, the K-nearest neighbor points of the tie points can be used as matching unit and an error equation can be constructed for every matching unit. The boresight angular error can then be solved iteratively and different flight lines can be aligned by applying the boresight angular error corrections.

If the boresight angular error is not well rectified, the geolocations of the laser points might have spatial displacements, which cause tie points in adjacent flight strips to not be spatially coincident. As shown in Figure 3, S₁ and S₂ are the same area in the two adjacent UAV scanning strips, but the boresight angular error caused geolocation displacement, and the centroids in dark black color do not coincide spatially each other.

2.4. Workflow of Our Proposed Method

2.4.1. Generation of the Intensity Images

The laser point intensity data is to be transformed onto a horizontal plane and then rasterized, and each grid cell is assigned mean intensity value of all the laser points that fall into current grid cell. The resolution of the grid cell can be set as the average point spacing, as for some grid cells that have null value (i.e., hollow), the strategy is to find the nearest grid cell that is not more than two pixels apart and to assign the nearest grid cell value to this null grid cell, which can enable the intensity image to be more smoothly and homogeneously. Finally, the intensity image of each strip can be obtained after the processing abovementioned.

2.4.2. Tie Point Extraction in 2D Space

SIFT is robust for changes in illumination or viewing angle and shows a good strong potential in antinoise and is widely used in target tracking, image mosaic, etc. [32]. The SIFT algorithm was adopted in our study to extract the key points from the intensity image of each strip and obtain initial tie points by matching key points from the adjacent images. It is inevitable that there will be some pairs of pseudo matching points in the initial tie points, and thus, the RANSAC algorithm is used to optimize the initial tie points [33]. The local affine transformation invariance is used as a constraint between the adjacent strips to eliminate pseudo matching points.

2.4.3. Refining Tie Point Sets in 3D Space

The tie points retrieved from the intensity images are mapped to 3D point cloud space, and a K-neighborhood search algorithm was applied to find the K-nearest neighboring points in the 3D point cloud space. The resultant K-nearest points are considered as a tie point set or matching units among the different LiDAR scanning strips. The value of K needs to be determined according to terrain conditions and laser point spacing, so that the K-nearest neighboring points can be located on the same surface as much as possible to avoid 3D leap. Since we only consider the radiation characteristics of targets in the extraction of tie point sets and do not take into account the geometric characteristics inside the neighborhood of the tie points, as well as the influence of observation noise, some tie point sets can be found in areas with obvious radiation characteristics such as buildings edges and vegetation canopies. However, these tie point sets with unstable geometric features are more likely to cause matching errors, and need further optimization therefore. In our study, a 3D plane is fitted for each tie point set firstly, and then the normal angle between the pair of planes derived from the tie point set is calculated. The pairs of the tie points with a normal angle greater than a certain threshold δ_n and a height difference greater than a certain threshold δ_h will be eliminated and the optimized tie point sets can be achieved finally.

2.4.4. Estimation of Boresight Angular Error Parameters

Step 1: Assigning initial values

The initial values of boresight angular parameters are given based on prior knowledge, since the number of the boresight angular error parameters is relatively small, therefore the initial values of the three boresight angular parameters can be all set to 0.

Step 2: Construction of the observation error equation

For the tie point set in the n-th LiDAR scanning strip, if there is a point that is named $P_i=(x_W^i\hspace{1em}{y}_W^i\hspace{1em}{z}_W^i)$ in it, then the coordinate of the point in IMU coordinate system can be derived from Equation (1). Thus its corrected coordinate in the mapping frame coordinate system can be derived from Equation (5) with known boresight angular parameter values, and the centroid coordinate of the tie point set can be derived with following equation.

$$\overline{F_n}=\left[\begin{array}{c}\overline{X_n}\\ \overline{Y_n}\\ \overline{Z_n}\end{array}\right]=\left[\begin{array}{c}\frac{1}{K}{\displaystyle \sum _{i=1}^K{X}_{Wn}^i}\\ \frac{1}{K}{\displaystyle \sum _{i=1}^K{Y}_{Wn}^i}\\ \frac{1}{K}{\displaystyle \sum _{i=1}^K{Z}_{Wn}^i}\end{array}\right](n=1,2)$$

(13)

where n denotes the number of the UAV flightlines. The residual of the observation error equation is

$$L=\overline{F_1}-\overline{F_2}$$

(14)

We first calculate the unknown coefficient for each point of the tie point set, and then average the coefficients as the coefficient matrix of the normal equations with the following two equations.

$$\begin{array}{c}\frac{\partial \overline{F_n}}{\partial \omega }=\frac{1}{K}{\displaystyle \sum _{i=1}^K\frac{\partial F_n^i}{\partial \omega }}\\ \frac{\partial \overline{F_n}}{\partial \phi }=\frac{1}{K}{\displaystyle \sum _{i=1}^K\frac{\partial F_n^i}{\partial \phi }}\\ \frac{\partial \overline{F_n}}{\partial \kappa }=\frac{1}{K}{\displaystyle \sum _{i=1}^K\frac{\partial F_n^i}{\partial \kappa }}\end{array}\}$$

(15)

$$A=\left[\begin{array}{ccc}\frac{\partial \overline{F_1}}{\partial \omega }-\frac{\partial \overline{F_2}}{\partial \omega }& \frac{\partial \overline{F_1}}{\partial \phi }-\frac{\partial \overline{F_2}}{\partial \phi }& \frac{\partial \overline{F_1}}{\partial \kappa }-\frac{\partial \overline{F_2}}{\partial \kappa }\end{array}\right]$$

(16)

If there are only two adjacent strips, only three observation error equations can be derived for each pair of tie points. For three or more LiDAR strips, three observation error equations can be listed between any two strips and then a total of 3*C²_n equations will be derived.

Step 3: Solving observation error equations

There are three unknown parameters in the unknown matrix, since each pair of tie points can list three equations, so the boresight angular error can be solved using at least one pair of tie point set. Many pairs of tie point sets will cause redundant observations and form a complete observation error equation. The approximate solution to the boresight angular error is obtained by Equation (12) in this section, and the three boresight angular parameters are updated as follows

$$\begin{array}{c}\omega =\omega +\Delta \omega \\ \phi =\phi +\Delta \phi \\ \kappa =\kappa +\Delta \kappa \end{array}\}$$

(17)

Step 4: Determining the termination of iterative running

Whether the residual is convergent and whether the parameter value keeps almost unchanged is used as a termination condition of the iterative running. If both are satisfied, the iterative calculation will be stopped and the estimated values of the boresight angular parameters will be output. If both are not satisfied, then the current estimated parameter values will be set as new values of the boresight angular parameters of the next iteration, and go back to step 2 to continue the iterative running. The flowchart for the proposed method is illustrated in Figure 4.

In the implementation of the technical flowchart in our study, the 3D tie points were classified into two groups, the training set and testing set. The training samples were extracted from the total samples of 3D tie points in a random way, and the rest of the total samples were used as the testing samples. The training samples were applied to estimate the parameters for the boresight angular error while the testing ones were prepared for the verification of the model, and the disparity between the tie point pairs was adopted to be an index for measuring the model fitting effect. The random sample dividing operation were performed 100 times in the experiment to guarantee the best estimate of the model parameters, and finally the model parameters of the best verification were chosen to achieve robust and reliable estimation of the boresight angular error.

3. Results

3.1. 2D Tie Points Extraction in the Case Study

The two LiDAR strips’ data was captured with the RIEGL VUX-1UAV laser scanner onboard on a Swiss-made ScoutB1-100 unmanned helicopter on 29 July, 2017 were de-noised with an anomaly detector such as Gaussian distribution statistics model to remove the elevation anomaly points in a local neighborhood. After that, the intensity images were generated by interpolating the laser intensity data, and the SIFT algorithm was applied to detect feature points on the two intensity images. A 128-dimensional feature descriptor was generated using the SIFT algorithm to describe the local features of the key points such as position, scale, and rotation. Next, the Euclidean distance between the two descriptors was calculated and used as a similarity metric of the two-image matching points. The matching strategy is as follows. Take a given key point in the query image and find the key point that is closest to the target one based on the ranking of the similarity. A certain number of pairs of tie points of the two images were able to be extracted out from the 2D intensity images. The RANSAC algorithm was applied to remove pseudo matching pairs with the condition of affine invariance. Finally a total of 63 pairs of tie points were extracted with the SIFT algorithm and 23 pairs of refined tie points were remained after removing pseudo matching points. It can be seen that the tie points are mainly located near the flat ground and the edge of the buildings. The spatial distribution is relatively uniform and the link lines between all pairs of refined tie points are basically parallel and show a good consistency.

3.2. 3D Tie Point Sets Construction with the Two Flight Lines

Considering the flat topography in the study area and the spatial distribution of the tie point pairs in 2D space, the thresholding values of the normal vector angle parameter δ_n of the tie point set and the elevation range parameter δ_h of each point set were set to 30° and 0.5 m, respectively. In the transformation of the tie points from 2D to 3D, the value of the core parameter K has a great influence on the final results (for details see Section 4.1). After a series of experiments and comparisons in the selection of optimal tie points in 3D space, the best performance can be derived when K is set to 200, the discrepancy of the tie point set is the best and the standard deviation is the smallest at the same time. Therefore, the value of K is set to 200 to extract the tie point set in 3D space finally, and a total of 12 pairs of optimal 3D tie points were retrieved in the study (Figure 5).

3.3. Correction of the Boresight Angular Error of the Two Strips

The boresight angular errors of the Scount-100/RIEGL VUX-1UAV system were resolved with an automation workflow, in which only the POS data and laser point cloud data needed to be inputs. The operation for randomly creating the training set and testing set from the total 3D tie point samples were conducted 100 times in the experiment of this study. Each training set was used to train the model and estimate the parameters, and each testing one was used to verify the parameters. The training and testing sets that achieved the best estimation of the model parameters were adopted and the resultant estimation of the boresight angular error were used to correct the laser point cloud data of the two flight lines. To validate our proposed method, the data obtained by the commonly used stepwise geometric method [15] were used as a comparison. It can be seen in

Table 2. Estimated boresight angle error parameters.

Method	Stepwise Geometric Method			Our Proposed Method
Parameter	ω	ϕ	κ	ω	ϕ	κ
Estimated value	−1.050°	−0.2580°	−0.7980°	−0.7384°	−0.2245°	−0.7219°

that the values of the three boresight angular parameters are relatively small in number, almost within 1°, which is consistent with the results in previous studies.

The resultant boresight angular error parameters were then substituted into Equation (5) to rectify the point cloud data of the two LiDAR scanning strips captured in Shihezi, Xinjiang, China on 29 July, 2017, and adjusted laser point data were achieved after the rectification using the parameters. A visual check of the corrections was illustrated in Figure 6, and it can be seen that the horizontal offset of the building facades have been effectively corrected after applying the rectification parameters.

Specifically it showed that our proposed method has a better performance than the stepwise geometric rectification method (e.g., 1# in Figure 6), but the stepwise geometric rectification method showed a slight better result in vertical direction (e.g., 2#). The displacement of the rooftop surfaces of the buildings were well aligned in horizontal direction and the rectification performance of the two methods showed few differences (e.g., 3#). The laser points of the vehicles have been well adjusted and became more coincident in the overlapping areas, and our proposed method slightly outperformed the stepwise geometric method (e.g., 4# and 5#). Because the diameter of the electricity poles is small and the spacing distance between the poles is relatively large, consequently, the rectification performance is not very distinct. However, the positional shift has been improved after the rectification (e.g., 6# and 7#). Generally speaking, our proposed method outperformed the stepwise geometric rectification method in most cases in this study, and has better automation and less ancillary data requirements such as GCPs, and thus can effectively improve the quality of laser point cloud data and enable good matching between adjacent flight lines.

3.4. Accuracy Assessment

The effect of the boresight angular error rectification on the position of the tie point set was analyzed quantitatively by using the absolute coordinate deviation of the tie point set before and after rectification processing. In Figure 7, the X axis indicates the sequence number of the tie point set, and the Y axis represents the discrepancy of the centroid coordinates of the tie point set. On the whole, the discrepancy of the XYZ coordinate before the rectification is basically between 0.04 m and 0.35 m, and the discrepancy of the rectified XYZ coordinate is basically less than 0.1 m, which shows that after the rectification the offset error between the tie points was effectively corrected. In addition, the discrepancy of Z axis is significantly smaller than that of the XY axes, indicating the rectification performance in vertical direction is better than in horizontal directions.

Due to the lack of accurate ground control points in study area, the RMSE statistic is chosen to evaluate the rectification accuracy of all the laser point cloud data. The planar RMSE of the point cloud data was calculated by projecting vertical wall points into horizontal plane to form a series of discrete points, based on which, the fitted residual error was calculated. The elevation RMSE was also calculated by projecting the roof points to vertical plane to form a series of discrete points, and the standard deviation of elevation of these points was calculated. The rooftop and vertical walls of the buildings in the study area were selected as planar and vertical reference planes. Totally six vertical reference planes and six horizontal rooftop planes were selected in each flight strip to perform the accuracy assessment (

Table 3. Accuracy assessment based on all the laser points.

Error	Planar RMSE/m			Elevation RMSE/m
Method	Raw data	Stepwise geometric method	Proposed method	Raw data	Step-wise geometric method	Proposed method
Strip 1	0.060	0.049	0.050	0.024	0.015	0.014
Strip 2	0.075	0.059	0.064	0.020	0.015	0.014
Average	0.068	0.054	0.057	0.022	0.015	0.014

).

It can be seen that after the rectification of boresight angular error, the accuracy in planar and elevation improved and the planar RMSE is 5.7 cm and decreases by 1.0 cm to 2.0 cm, while the elevation RMSE is ~1.4 cm and decreases by 0.5 cm to 1.0 cm. The RMSE reduction in elevation is slightly smaller than that in X–Y plane. In addition, our rectification method can achieve better correction in elevation than the stepwise geometric rectification method in the case study. As for the two LiDAR strips, the rectification of strip 1 is better than that of strip 2 in general.

4. Discussion

It can be seen in the case study that the proposed method can remove most of the boresight angular error caused by the unstable low-altitude UAV LiDAR system, and can achieve good matching of the two adjacent flight lines based on the laser intensity information. Compared to the stepwise geometric correction method, our method does not require any ground control points, feature objects, or raw observations of the laser scanner, and thus has larger degree of automation. However, the parameterization in our method may have influence on the final result, so we will discuss about it in the next subsections.

4.1. The Influence of Parameter K

The K-nearest neighbor has a great impact on the spatial distribution of 3D tie points, and the proposed method requires extremely good quality tie point sets for accurately solving the error model. Thus, it is necessary to understand how K value influences on the final result. To empirically explore the influence of K, we selected a series of K values that fall between 0 and 500, and tried to resolve the parameters for the boresight angular error for each K value. Some results of the K values are shown in

Table 4. The boresight angle errors of different K values.

K value	Match count	ω/°	ϕ/°	κ/°
20	13	-0.7302	-0.3353	-2.3283
40	12	-0.7366	-0.2530	-1.1801
60	12	-0.7384	-0.2507	-1.1319
80	12	-0.7373	-0.2451	-1.0975
100	12	-0.7395	-0.2382	-0.9498
200	12	-0.7384	-0.2245	-0.7219
300	11	-0.7492	-0.2951	-1.6186
400	11	-0.7500	-0.2955	-1.6275
500	10	-0.7498	-0.2490	-1.6450

. It can be seen that the most sensitive parameter of the boresight angular error is κ, second ϕ, and third ω. Applying each solution under different K values to make a rectification for the tie point sets and the disparity statistical result is shown in Figure 8. It is found that regardless of the value of K, the minimum value of disparity tends to be 0 and shows little change. However, the maximum value of disparity is not stable, and shows an obvious trend of first increasing and then decreasing. The mean disparity seems to be stable and shows a general trend, decreasing first and followed by an increase, and reaches its minimum value when the K value is 200. The standard deviation of the disparity shows similar trends and reaches the minimum value when the K is 200. Empirically, the optimal value of K is determined by comparing the mean value and the standard deviation of disparity. Therefore, the rectification was processed with the optimal K value of 200.

4.2. The Influence of Calibration Parameter on Geolocation Error

Besides the K value, it is necessary to examine the influence of the boresight angular error on the positioning error of laser point cloud data quantitatively. Assuming that the remaining errors are all 0, according to the basic principle of laser point cloud geolocation, the influence of the boresight angular error parameters can be derived from the difference between Equation (5) and Equation (1), and is expressed in Equation (18).

$$\begin{array}{l}\left[\begin{array}{c}{e}_X\\ e_Y\\ e_Z\end{array}\right]=R_W{R}_N(R_W-I)R_I{R}_L\left[\begin{array}{c}0\\ 0\\ \rho \end{array}\right]\\ =R_W{R}_N\left[\begin{array}{ccc}0& -\kappa & \phi \\ \kappa & 0& -\omega \\ -\phi & \omega & 0\end{array}\right]·\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]·\left[\begin{array}{ccc}0& 0& 0\\ 0& \cos \theta & -\sin \theta \\ 0& \sin \theta & \cos \theta \end{array}\right]·\left[\begin{array}{c}0\\ 0\\ \rho \end{array}\right]\\ =R_W{R}_N\left[\begin{array}{c}\kappa \rho \sin \theta +\phi \rho \cos \theta \\ -\omega \rho \cos \theta \\ -\omega \rho \sin \theta \end{array}\right]=R_W{R}_{N}H\left[\begin{array}{c}\kappa \tan \theta +\phi \\ -\omega \\ -\omega \tan \theta \end{array}\right]\end{array}$$

(18)

where θ is the laser scanning angle; H is the altitude of the flight; ω, ϕ, and κ are boresight angular error parameters in three directions of rolling, pitching, and heading, respectively; and ρ denotes the distance measured by the laser scanner. It can be understood in Equation (18) that the geolocation error of a LiDAR system is proportional to the altitude H, and the higher the altitude is, the larger the planar and elevation error. When the laser scanning angle θ is constant, the pitching direction φ and heading direction κ together affect the X-direction error of the point cloud data, and the error increases with the increase of these two parameters. The rolling direction ω mainly affects the Y-direction and Z-direction errors of the point cloud data, and both are proportional to ω, and as the error increases, ω increases. It is worth noting that the scanning angle is also an important factor influencing the positioning error of the LiDAR system, and has only a direct effect on the X and Z direction errors of the laser point cloud data.

4.3. Influence of the Initial Values on Model Convergence

In addition to the parameters abovementioned, the convergence speed also has a significant impact on the practicability and robustness of the model for estimating the boresight angular error. Considering the assumption that the boresight angular error is relatively small in number, thus we preset the initial values of w, ϕ, and κ to zero, and the parameters of the boresight angular error can be solved after three iterations. If this is not true and the parameters are big in number, does the initial values of w, ϕ, and κ affect the convergence speed of the model? A series of initial values was preset and used to test the convergence speed (

Table 5. The relationship between the convergence speed and the initial values.

Initial Value			Iteration Count	Converges to Same Value
ω/°	ϕ/°	κ/°
0	0	0	3	—
10	0	0	4	yes
0	10	0	3	yes
0	0	10	3	yes
10	10	0	3	yes
10	0	10	3	yes
0	10	10	4	yes
10	10	10	4	yes
30	30	30	4	yes
60	60	60	6	yes

). It can be seen in Table 5 that the initial values of the boresight angular error parameters have slight influence on the convergence speed, and three or four iterations can achieve the convergence in our experiment. Even if the initial values were intentionally set to a very large value, i.e., 60°, the convergence can also be achieved only after six iterations. Therefore, empirically speaking, the initial values of the boresight angular error parameters have slight influence on the convergence speed. In other words, our proposed method is not sensitive to the initial values and possesses robustness and stability.

5. Conclusions

This paper presents a new method for the boresight angular error rectification of a UAV LiDAR system based on the laser intensity information. Our proposed method has been verified with tens of millions of laser point cloud data acquired by the Scount-100/RIEGL VUX-1UAV LiDAR system in a farmland located in the northwestern of Shihezi city, Xinjiang, China. A comparison with conventional stepwise geometric rectification method [15] was also conducted with the same data sets. It can be concluded that the boresight angular error is one of the main error sources leading to the positional error between different scanning strips of UAV LiDAR data. The boresight angular errors in the UAV LiDAR system used in our study were estimated by our proposed method, and the angular parameters are ω=-0.7385°, φ=-0.2245°, and κ=-0.7219°. After the rectification, the planar RMSE is 5.7 cm and decreased by 1.0 to 2.0 cm, and the elevation RMSE is 1.4 cm and decreased by 0.5 to 1.0 cm. It is also found that K value in K-nearest analysis has great influence on the estimation of the boresight angular error parameters. Due to the difficulty in the determination of a best K value theoretically, an empirical solution has been proposed in our study, simply put, it can be determined through comparing the rectification performance under a series of different K values. In this experiment, the best K value of 200 was adopted for the implementation of our proposed method, and can achieve best adjustment of the two LiDAR scanning strips. The sensitivity analysis of the angular parameters shows that ω is the weakest, followed by φ, and κ is the strongest. It should be noted that our proposed method is not sensitive to the initial values of ω, φ, and κ. In spite of different initial values, the resolving of the boresight angular error model can be accomplished in less than ten iterations. Thus our proposed method shows good robustness and stability in practice.

One limitation of this study is that the testing was just based on the two UAV LiDAR strips, and the experimental area is a relatively flat area and the main ground object types are vegetation and road pavement. Consequently, there is a potential uncertainty in the extraction of tie points. In the future, more validation work needs to be done on a larger dataset with more UAV flight lines. For example, testing work can be expanded to mining areas, built-up areas, and other regions with various topographical conditions and ground object types. In addition, in this paper we only considered the boresight angular error of the point cloud data, lacking of the analysis and comparison of other potential error sources. Future work will be focused on the combination of other error sources to design robust algorithms and make more improvement in the practicability of our proposed method.