Airborne LiDAR Intensity Correction Based on a New Method for Incidence Angle Correction for Improving Land-Cover Classiﬁcation

: Light detection and range (LiDAR) intensity is an important feature describing the char-acteristics of a target. The direct use of original intensity values has limitations for users, because the same objects may have different spectra, while different objects may have similar spectra in the overlapping regions of airborne LiDAR intensity data. The incidence angle and range constitute the geometric conﬁguration of the airborne measurement system, which has an important inﬂuence on the LiDAR intensity. Considering positional shift and rotation angle deviation of the laser scanner and the inertial measurement unit (IMU), a new method for calculating the incident angle is presented based on the rigorous geometric measurement model for airborne LiDAR. The improved approach was applied to experimental intensity data of two forms from a RIEGL laser scanner system mounted on a manned aerial platform. The results showed that the variation coefﬁcient of the intensity values after correction in homogeneous regions is lower than that obtained before correction. The overall classiﬁcation accuracy of the corrected intensity data of the ﬁrst form (amplitude) is signiﬁcantly improved by 30.01%, and the overall classiﬁcation accuracy of the corrected intensity data of second form (reﬂectance) increased by 18.21%. The results suggest that the correction method is applicable to other airborne LiDAR systems. Corrected intensity values can be better used for classiﬁcation, especially in more reﬁned target recognition scenarios, such as road mark extraction and forest monitoring. This study provides useful guidance for the development of future LiDAR data processing systems.


Introduction
Airborne LiDAR sensors can obtain high-precision three-dimensional (3D) coordinates of objects, and it is well known that height information (Z value) of LiDAR data can be used for ground object recognition and extraction [1]. The potential of LiDAR information is made even greater with the inclusion of intensity data [2][3][4][5]. In addition to acquiring geometric information, LiDAR sensors can also record the intensity information of the surface targets. The intensity values record the amplitude of the return signal from the illuminated object, which usually can be the analog electrical signal output from the photodetector or the digitized waveform [1]. Commercial LiDAR sensors generally use Nd:YAG lasers with a wavelength of 1064 nm. The spectral reflectance of different materials has separability in the near-infrared spectral range [6]. The use of LiDAR intensity data relationship between the intensity data with AGC-on and AGC-off in the same area, and found that the model can be useful for many of the natural targets with a reflectance around 20-30%, while some problems were encountered with low-reflectance targets. Korpela [36] compensated AGC from the Leica ALS50 sensors, and found that it can be optimized for maximum classification performance. Yan and Shaker [13] proposed a sub-histogram matching technique based on Gaussian mixture model to align the overlapping intensity data among different LiDAR data strips and minimized the effects of AGC for intensity.
Range and the incidence angle constitute the geometric configuration of the airborne measurement system, which has an important influence on the laser intensity. Greater incidence angles (θ) typically result in less incident laser energy being backscattered in the direction of the receiver, thereby affecting the intensity. The emitted and received laser energy decays as greater range (R), thus influence the intensity (see Figure 1). The incidence angle and range were usually considered in the correction process [37][38][39]. Range between the sensor and individual returns can be obtained relatively easily, which is contained within the trajectory file of the plane [40]. Some researchers have calculated the incident angle by different methods to correct intensity. One is to use scan angle instead of the incident angle for intensity correction. The scan angle was feasible for the flat ground, but there were errors in the case of rugged terrain [22,33]. Yan et al. [22] improved the estimation of this angle by calculating the scan angle and surface slope and aspect, and found that this incidence angle would lead to the effects of overcorrection for steep slopes [41]. Yan and Shaker [13] further proposed an approach by using a slope threshold to select either the scan angle or the incidence angle in the radar range equation for intensity correction. Yi et al. [42] improved the calculation of this angle by using the object tilt angle and the scan angle through determination of the laser scanning direction. Another method calculated approximately the incident angle between the surface normal and each laser pulse vector from the instantaneous 3D coordinates of the aircraft, and the 3D coordinates of point cloud. The instantaneous 3D coordinates of the aircraft was from position and orientation system (POS) integrated by global positioning system (GPS) and inertial measurement unit (IMU). The laser scanner and IMU have position shift and rotation angle deviation. However, based on the rigorous geometric measurement model of airborne LiDAR, each laser pulse vector should be calculated from the origin of the laser scanner center, not from the center of the aircraft.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 22 normalized LiDAR intensity by a liner function with range and AGC, and showed that the model can improve the separability of lichens from other surface materials. Vain et al. [35] built linear empirical models based on the statistical relationship between the intensity data with AGC-on and AGC-off in the same area, and found that the model can be useful for many of the natural targets with a reflectance around 20-30%, while some problems were encountered with low-reflectance targets. Korpela [36] compensated AGC from the Leica ALS50 sensors, and found that it can be optimized for maximum classification performance. Yan and Shaker [13] proposed a sub-histogram matching technique based on Gaussian mixture model to align the overlapping intensity data among different Li-DAR data strips and minimized the effects of AGC for intensity. Range and the incidence angle constitute the geometric configuration of the airborne measurement system, which has an important influence on the laser intensity. Greater incidence angles (θ) typically result in less incident laser energy being backscattered in the direction of the receiver, thereby affecting the intensity. The emitted and received laser energy decays as greater range (R), thus influence the intensity (see Figure 1). The incidence angle and range were usually considered in the correction process [37][38][39]. Range between the sensor and individual returns can be obtained relatively easily, which is contained within the trajectory file of the plane [40]. Some researchers have calculated the incident angle by different methods to correct intensity. One is to use scan angle instead of the incident angle for intensity correction. The scan angle was feasible for the flat ground, but there were errors in the case of rugged terrain [22,33]. Yan et al. [22] improved the estimation of this angle by calculating the scan angle and surface slope and aspect, and found that this incidence angle would lead to the effects of overcorrection for steep slopes [41]. Yan and Shaker [13] further proposed an approach by using a slope threshold to select either the scan angle or the incidence angle in the radar range equation for intensity correction. Yi et al. [42] improved the calculation of this angle by using the object tilt angle and the scan angle through determination of the laser scanning direction. Another method calculated approximately the incident angle between the surface normal and each laser pulse vector from the instantaneous 3D coordinates of the aircraft, and the 3D coordinates of point cloud. The instantaneous 3D coordinates of the aircraft was from position and orientation system (POS) integrated by global positioning system (GPS) and inertial measurement unit (IMU). The laser scanner and IMU have position shift and rotation angle deviation. However, based on the rigorous geometric measurement model of airborne LiDAR, each laser pulse vector should be calculated from the origin of the laser scanner center, not from the center of the aircraft. There exists a gap in calculating the incident angle in the previous studies. This paper proposes a new method for calculating the incident angle and correcting intensity values for improving land cover classification. The paper evaluates the consistency of Remote Sens. 2021, 13, 511 4 of 21 intensities before and after correction in a homogeneous surface, and further addresses the classification accuracy based on airborne LiDAR intensity before and after correction. The study aims to enhance the utilization value of LiDAR intensity data and provide potential information for classification, especially in more refined target recognition scenarios, such as road mark extraction and forest monitoring. The study expands the knowledge base of correction and application of airborne LiDAR intensity data.
Following this introduction, Section 2 presents an approach to correcting the effects of incidence angle, and Section 3 introduces experiments and datasets. The results, discussion, and recommendations for future work are presented in Section 4. Conclusions are presented in Section 5.

Radar (Range) Equation
A laser scanner sends laser pulses, which return to the laser receiver. This physical process can be described quantitatively by the radar (range) equation [43]. The received laser energy (P r ) can be expressed as a product function related to transmitted signal power of optical transmission characteristics, laser atmospheric transmission system, target characteristics, and laser receiver parameters.
where P r is the received laser energy, P E is the transmitted laser energy, G is the gain factor of the antenna, R is the range between target and sensor, σ is the effective target cross-section, D is the receiver aperture diameter, η sys is the system transmission factor, and η atm is the atmospheric transmission factor. LiDAR intensity represents the energy and power of the received signal, amplitude, or reflectance in different laser scanner.
Because some parameters are difficult to obtain, Equation (1) is further simplified in previous studies [22,25]. Finally, the derived formula is radar range Equation (2) with the assumption of extended lambert target.
where ρ is the spectral reflectance of the surface. P E , D, η sys and η atm are assumed to be constant C during the same flight. This leads to Equation (3).
The LiDAR intensity is directly proportional to the amount of received laser energy Pr, as reported in Höfle and Pfeifer [25] and Yan et al. [22]. Therefore, LiDAR intensity is proportional to the spectral reflectance ρ, which is related to the surface properties. LiDAR intensity is also proportional to cos θ, and inversely proportional to squared range R. In theory, the intensity of the same target should be the same in different flight lines of the same scanning task. However, the range and the incident angle θ are different in different flight lines, resulting in difference intensity values of the same target. It is necessary to reduce or eliminate the effects of the range and incidence angle for extracting intensity value proportional to the reflectance of the target. The following section describes how to calculate the laser incidence angle.

A New Calculation Method for Incident Angle
The laser incidence angle (direction of reflection) is the angle between the instantaneous laser beam and the surface normal. Here the laser pulse vector is calculated by converting the laser original coordinates (X laser , Y laser , Z laser ) into projection coordinates (X, Y, Z). The proposed calculation method for incident angle is described below.
Firstly, the laser original coordinates (X laser , Y laser , Z laser ) are transformed into the coordinate of the IMU system: where (X I MU , Y I MU , Z I MU ) is the coordinate of the IMU system. The laser coordinate system and the IMU coordinate system are adjusted to be substantially parallel by matrix transformation of R 0 .
where R 1 is a rotation matrix, which represents the rotation transformation relationship between the three coordinate axes of the laser scanner and the three coordinate axes of the IMU. Ω, Φ, K are attitude angles between laser scanner and IMU. Then, the coordinates of the IMU system (X I MU , Y I MU , Z I MU ) are transformed into the projection coordinates (X, Y, Z): cos γ cos α + sin β sin γ sin α cos β sin α − sin γ cos α + sin β cos γ sin α − cos γ sin α + sin β sin γ cos α cos β cos α sin γ sin α + sin β cos γ cos α cos β sin γ − sin β cos γ cos β   where R 2 is a rotation matrix, which represents the rotation transformation relationship between the three coordinate axes of the IMU and the three coordinate axes of the geographic coordinates. γ is the roll angle by the IMU, β is the pitch angle by the IMU, and is the difference angle between the heading angle ψ and the meridian convergence angle ε [44]. Through the above transformation, the laser original coordinate system is parallel to the projection coordinates system. The laser pulse vector from each laser footprint is (−X, −Y, −Z), and the cosine of the laser incidence angle θ is calculated: where (n 1 , n 2 , n 3 ) is the surface normal. The following section describes how to calculate the surface normal.
In the above progress, each instantaneous laser beam is calculated from the origin of the laser scanner center, not the center of the aircraft calculated by the previous researches. This could be an improvement in applying rigorous geometric theory for LiDAR intensity correction so it will fill in the previous researches. Because the instantaneous 3D coordinates of the aircraft was from position and orientation system (POS) integrated by GPS and IMU. The laser scanner and IMU have position shift and rotation angle deviation.

The Surface Normal
In order to obtain the surface normal corresponding to each laser footprint, the surrounding points are taken for plane fitting, and the normal vector of the corresponding point is calculated by the fitting equation. The calculation of the surface normal is as follows: Firstly, the K-Nearest Neighbor (KNN) method is used to find the nearest N points for each point, and then the height differences between the current point and the neighboring points are calculated. A threshold value of 0.4 m is set to remove the neighboring points higher than the threshold value, and the set of points meeting the threshold condition is used as the initial set of points. The setting of the elevation threshold is mainly to reduce the calculation error of the surface normal of the building edge points and the ground points nearby. If the number of points in the initial set of points is less than 3, the surface normal of the point is (0, 0, 0). Then the initial set of points is fitted to a plane, and the normal vector of the plane at this point is calculated by the covariance matrix, singular value decomposition, and the eigenvector corresponding to the minimum singular value. The outliers for the cosine of the incident angle of the laser points are post-processed and replaced with the average of the neighboring points. Figure 2 shows the workflow of the radiometric correction of the intensity data for assessment of land cover classification. of the aircraft was from position and orientation system (POS) integrated by GPS and IMU. The laser scanner and IMU have position shift and rotation angle deviation.

The Surface Normal
In order to obtain the surface normal corresponding to each laser footprint, the surrounding points are taken for plane fitting, and the normal vector of the corresponding point is calculated by the fitting equation. The calculation of the surface normal is as follows: Firstly, the K-Nearest Neighbor (KNN) method is used to find the nearest N points for each point, and then the height differences between the current point and the neighboring points are calculated. A threshold value of 0.4 m is set to remove the neighboring points higher than the threshold value, and the set of points meeting the threshold condition is used as the initial set of points. The setting of the elevation threshold is mainly to reduce the calculation error of the surface normal of the building edge points and the ground points nearby. If the number of points in the initial set of points is less than 3, the surface normal of the point is (0, 0, 0). Then the initial set of points is fitted to a plane, and the normal vector of the plane at this point is calculated by the covariance matrix, singular value decomposition, and the eigenvector corresponding to the minimum singular value. The outliers for the cosine of the incident angle of the laser points are post-processed and replaced with the average of the neighboring points. Figure 2 shows the workflow of the radiometric correction of the intensity data for assessment of land cover classification.

RIEGL Laser Scanning Equipment
In this study, we correct the intensity of RIEGL laser scanning equipment. RIEGL provides two types of intensity data: calibrated amplitude and relative reflectance. The calibrated amplitude is given relative to the amplitude of an echo signal at the detection threshold of the instrument. The value of the amplitude reading is a ratio, given in the units of decibel (dB). The formula of the amplitude:  In this study, we correct the intensity of RIEGL laser scanning equipment. RIEGL provides two types of intensity data: calibrated amplitude and relative reflectance. The calibrated amplitude is given relative to the amplitude of an echo signal at the detection Remote Sens. 2021, 13, 511 7 of 21 threshold of the instrument. The value of the amplitude reading is a ratio, given in the units of decibel (dB). The formula of the amplitude: where A dB is the amplitude in decibel, P echo is the optical input power, P DL is the minimum detectable input power. The amplitude of the echo signal reaching the laser scanner depends on a number of parameters, including system parameters like the emitted laser pulse peak power and the receiver aperture, but also including target parameters like the target's reflectance and range. A dB suffers from scan range dependence, so the relative reflectance as an additional attribute is provided by RIEGL. The relative reflectance is a ratio of the actual amplitude of that target to the amplitude of a white flat target at the same range, orientated orthonormal to the beam axis, and with a size in excess of the laser footprint. The value of the reflectance is also a ratio, given in the units of decibel: where ρ rel is relative reflectance in dB, A dB is the calibrated amplitude, A dB,Re f (R) is amplitude of a reference target at range R. The reflectance is nearly independent from the scan range [45]. According to Equations (11) and (12), the amplitude and reflectance of RIEGL laser scanning equipment are in dB units, which are converted as follows: where A dB is the calibrated amplitude, I amp is the converted amplitude of target.
where ρ rel is the relative reflectance, I re f is the converted reflectance of target.
Finally, we correct the intensity of first type (the calibrated amplitude) by eliminating the effects of incident angle and range. The equation is: where I Camp is the intensity of first type after correction, R is the range between the sensor and the target, R s is the reference distance, cos θ is the cosine of incident angle (direction of reflection). The range is calculated by laser original coordinates (X laser , Y laser , Z laser ), which are decoded by the RIEGL RXP file. The intensity of second type (the relative reflectance) is corrected by eliminating the effects of laser incident angle. The equation is: where I Cre f is the intensity of second type after correction, and cos θ is the cosine of incident angle (direction of reflection).

Study Area and Dataset
The study area is located in Xuchang City, Henan Province, China. The flight experiment was carried out on 16 March 2016. LiDAR raw data were acquired by RIEGL laser scanner system mounted on a manned aerial platform. Three strips of experimental area were shown in Figure 3. Flight strip 1 was scanned from south to north, flight strip 2 was scanned from north to south, and flight strip 3 was scanned from east to west. Three strips  Through coordinate transformation, LiDAR raw data ( , , , and ) decoded by RXP file data were fused into discrete LiDAR point cloud data of ( , , , ) and ( , , , ) in the coordinate system of WGS-84 UTM Zone 49N. We corrected the intensity of the two types by Equations (15) and (16), and obtained the corrected LiDAR point clouds ( , , , ) and ( , , , ). The LiDAR intensity data of two types before correction and after correction were converted into intensity raster images with 0.5 m resolution.

Evaluation Method of Intensity Correction
We used statistical analysis of samples to assess the intensity values in homogeneous areas before and after correction. Samples from homogeneous road surfaces in flight strips 1 and 2 were collected. The sample data of different scan ranges and scan angles from nadir to edges cover the entire sample area. Since the roads are flat, the incidence angle is nearly equal to the scan angle. A total of 30 samples from small rectangular areas with different scan ranges and scan angles were selected from the road surfaces. Each small rectangular area included 25-40 points. The standard deviation of the scan ranges for these points is less than 0.4 m, and the standard deviation of the scan angles for these points is less than 3°, similar to the study by Oh [46]. Figures 4 and 5 show the sampling areas in flight strips 1 and 2, respectively, where the red boxes represent the rectangular sample areas and the blue line represents the flight path. Through coordinate transformation, LiDAR raw data (X laser , Y laser , Z laser , A dB and ρ rel ) decoded by RXP file data were fused into discrete LiDAR point cloud data of X, Y, Z, I amp and X, Y, Z, I re f in the coordinate system of WGS-84 UTM Zone 49N. We corrected the intensity of the two types by Equations (15) and (16), and obtained the corrected LiDAR point clouds X, Y, Z, I Camp and X, Y, Z, I Cre f . The LiDAR intensity data of two types before correction and after correction were converted into intensity raster images with 0.5 m resolution.

Evaluation Method of Intensity Correction
We used statistical analysis of samples to assess the intensity values in homogeneous areas before and after correction. Samples from homogeneous road surfaces in flight strips 1 and 2 were collected. The sample data of different scan ranges and scan angles from nadir to edges cover the entire sample area. Since the roads are flat, the incidence angle is nearly equal to the scan angle. A total of 30 samples from small rectangular areas with different scan ranges and scan angles were selected from the road surfaces. Each small rectangular area included 25-40 points. The standard deviation of the scan ranges for these points is less than 0.4 m, and the standard deviation of the scan angles for these points is less than 3 • , similar to the study by Oh [46]. nearly equal to the scan angle. A total of 30 samples from small rectangular areas w different scan ranges and scan angles were selected from the road surfaces. Each sm rectangular area included 25-40 points. The standard deviation of the scan ranges these points is less than 0.4 m, and the standard deviation of the scan angles for th points is less than 3°, similar to the study by Oh [46]. Figures 4 and 5 show the sampli areas in flight strips 1 and 2, respectively, where the red boxes represent the rectangu sample areas and the blue line represents the flight path.  Because of some artificial errors in the sampling process of natural objects, it is not possible that all sampling points have unique intensity values on uniform surfaces. Meanwhile, LiDAR intensity values of the two types have been stretched after correction, leading to different dimensions, so we used the median and the coefficient of variation (CV) to assess the correction effects. The formula of the coefficient of variation is: where CV denotes the dispersion degree of the sample data, σ is the standard deviation of the sample value, and μ is the mean of the sample values.
We assessed the classification accuracies for land cover types based on uncorrected and corrected intensity data of two types by using the random forest method. Six classes were differentiated: building, pavement/road, tree, cropland, grass, and bare soil. Building, pavement/road, and tree are predominant land-cover types, while cropland, grass, and bare soil are minor classes in the scenes of the study area (see Figure 3). Google map image data with a resolution of 0.25 m was used as a reference to randomly select samples of six classes (see Table 1). Table 1 shows total number of samples for six land cover types and independent-sample Kruskal-Wallis test for differences. Moreover, 70% of samples were used as training data and 30% of samples were used for verification. There is no overlap between the training samples and the verification samples. Overall accuracy were used to evaluate the classification accuracy of the whole data set. Producer's accuracy (PA) and user's accuracy (UA) were used to assess the classification accuracy of each class. We added overall class accuracy (CA) for a single class defined by Singh et al. [12], which takes into account the misclassification error and omission error of each class.  Because of some artificial errors in the sampling process of natural objects, it is not possible that all sampling points have unique intensity values on uniform surfaces. Meanwhile, LiDAR intensity values of the two types have been stretched after correction, leading to different dimensions, so we used the median and the coefficient of variation (CV) to assess the correction effects. The formula of the coefficient of variation is: where CV denotes the dispersion degree of the sample data, σ is the standard deviation of the sample value, and µ is the mean of the sample values.
We assessed the classification accuracies for land cover types based on uncorrected and corrected intensity data of two types by using the random forest method. Six classes were differentiated: building, pavement/road, tree, cropland, grass, and bare soil. Building, pavement/road, and tree are predominant land-cover types, while cropland, grass, and bare soil are minor classes in the scenes of the study area (see Figure 3). Google map image data with a resolution of 0.25 m was used as a reference to randomly select samples of six classes (see Table 1). Table 1 shows total number of samples for six land cover types and independent-sample Kruskal-Wallis test for differences. Moreover, 70% of samples were used as training data and 30% of samples were used for verification. There is no overlap between the training samples and the verification samples. Overall accuracy were used to evaluate the classification accuracy of the whole data set. Producer's accuracy (PA) and user's accuracy (UA) were used to assess the classification accuracy of each class. We added overall class accuracy (CA) for a single class defined by Singh et al. [12], which takes into account the misclassification error and omission error of each class.

Assessment of Homogeneous Areas
LiDAR intensity data of the two types before and after correction are shown in Figures 6 and 7. Visual analyses suggest that the intensity values at the nadir are significantly higher than those at the edges before correction, which also proved that intensity correction is meaningful. This can be explained by the increase in the range at the edges, and the decrease of received laser energy. Figure 6a shows that the intensity values of the overlapping regions of strips 1, 2, and 3 were not consistent. Theoretically, when scanning the same target, the intensity values of the same target were close to each other in different flight strips. In the paper, in order to highlight the intensity correction effects of proposed incident angle, we only collected the samples from roads, not other classes. Different scan ranges and incident angles from nadir to edges cover the entire sample area in flight strips 1 and 2.

Assessment of Homogeneous Areas
LiDAR intensity data of the two types before and after correction are shown in Figures 6 and 7. Visual analyses suggest that the intensity values at the nadir are significantly higher than those at the edges before correction, which also proved that intensity correction is meaningful. This can be explained by the increase in the range at the edges, and the decrease of received laser energy. Figure 6a shows that the intensity values of the overlapping regions of strips 1, 2, and 3 were not consistent. Theoretically, when scanning the same target, the intensity values of the same target were close to each other in different flight strips. In the paper, in order to highlight the intensity correction effects of proposed incident angle, we only collected the samples from roads, not other classes. Different scan ranges and incident angles from nadir to edges cover the entire sample area in flight strips 1 and 2.  In order to further evaluate the results of intensity correction, we compared the intensity statistics of the road before and after correction. In flight strip 1, the median values of the intensity of the first type for each sample were shown in Figure 8a,b. The intensity values before correction I amp were between 2 and 8 while the intensity values after correction I Camp were from 5 to 8, which indicated that the intensity values of the first type after correction I Camp were more concentrated. The median values of intensity of the second type for each sample were shown in Figure 9a,b. The intensity values of the second type before correction I re f were between 1200 and 2600, while the intensity values after correction I Cre f were between 1500 and 2600. The results show that the intensity values of the second type after correction I Cre f were more concentrated, so were the corrected intensity values of the first type I Camp . Based on the results in Table 2, the CV of the intensity of the first type before correction I amp was 0.36, the CV of the intensity of the first type after correction I Camp was 0.14, and the ratio of variation coefficient after correction to that before correction was 0.39 (less than 1), indicating that the theoretical model using the proposed correction method was effective. The CV of the intensity of the second type before correction I re f was 0.20; the CV of the intensity of the second type after correction I Cre f was 0.14, the ratio of variation coefficient after correction to that before correction was 0.7, which also showed that the theoretical correction using the proposed incidence angle was effective. Moreover, it was found that the CV of the intensity values of the second type I re f was less than that of the first type I amp , the relative reflectance before correction is more useful than the calibrated amplitude provided by the RIEGL system. In order to further evaluate the results of intensity correction, we compared the intensity statistics of the road before and after correction. In flight strip 1, the median values of the intensity of the first type for each sample were shown in Figure 8a,b. The intensity values before correction were between 2 and 8 while the intensity values after correction were from 5 to 8, which indicated that the intensity values of the first type after correction were more concentrated. The median values of intensity of the second type for each sample were shown in Figure 9a,b. The intensity values of the second type before correction were between 1200 and 2600, while the intensity values after correction were between 1500 and 2600. The results show that the intensity values of the second type after correction were more concentrated, so were the corrected intensity values of the first type . Based on the results in Table 2, the CV of the intensity of the first type before correction was 0.36, the CV of the intensity of the first type after correction was 0.14, and the ratio of variation coefficient after correction to that before correction was 0.39 (less than 1), indicating that the theoretical model using the proposed correction method was effective. The CV of the intensity of the second type before correction was 0.20; the CV of the intensity of the second type after correction was 0.14, the ratio of variation coefficient after correction to that before correction was 0.7, which also showed that the theoretical correction using the proposed incidence angle was effective. Moreover, it was found that the CV of the intensity values of the second type was less than that of the first type , the relative reflectance before correction is more useful than the calibrated amplitude provided by the RIEGL system.    In flight strip 2, the intensity of the first type before correction and after correction were 2 to 8 and 5 to 8, respectively. The intensity of the second type before correction and after correction were 1200 to 2600 and 1600 to 2600, respectively (see Figures 10 and 11). The results show that the intensity values of the two types ( and ) were more concentrated after correction, which was in line with the conclusions obtained in flight strip 1, indicating that the correction was effective. We further compared the CV of the intensity before and after correction in Table 3, the CV of the intensity of the first type before correction and after correction was 0.30 and 0.13, the ratio of variation coefficient after correction to that before correction was 0.43. The CV of the intensity of the second type before correction and after correction was 0.18 and 0.12, the ratio of variation coefficient after correction to that before correction was 0.67 (less than 1). These results are also in line with the results from flight strip 1, the smaller the CV and ratio of coefficient variation, the better the correction results, which indicated that the improved theoretical correction by using proposed incidence angle is effective.  In flight strip 2, the intensity of the first type before correction I amp and after correction I Camp were 2 to 8 and 5 to 8, respectively. The intensity of the second type before correction I re f and after correction I Cre f were 1200 to 2600 and 1600 to 2600, respectively (see Figures 10 and 11). The results show that the intensity values of the two types (I Camp and I Cre f ) were more concentrated after correction, which was in line with the conclusions obtained in flight strip 1, indicating that the correction was effective. We further compared the CV of the intensity before and after correction in Table 3, the CV of the intensity of the first type before correction I amp and after correction I Camp was 0.30 and 0.13, the ratio of variation coefficient after correction to that before correction was 0.43. The CV of the intensity of the second type before correction I re f and after correction I Cre f was 0.18 and 0.12, the ratio of variation coefficient after correction to that before correction was 0.67 (less than 1). These results are also in line with the results from flight strip 1, the smaller the CV and ratio of coefficient variation, the better the correction results, which indicated that the improved theoretical correction by using proposed incidence angle is effective.   According to the median intensity values of the sample data before and after correction in flight strip 1 and 2 in Figures 8-11, the intensity variation of the first type before correction was larger than the intensity variation of the second type , which partly explains the range correction for the intensity of the first type is necessary in addition to incidence angle correction. The intensity of the first type after range and incidence angle correction is a value directly proportional to the spectral reflectance. Although the intensity of the second type has been calibrated by the RIEGL manufacturer, it is still necessary to eliminate the influence of incident angle, obtaining the intensity value proportional to the spectral reflectance.

Comparison of Land Cover Classification Performance Based on Airborne LiDAR Intensity before and after Correction
Using the random forest method, which is one of the best classification methods in machine learning [47], the classification accuracies for land cover based on the uncorrected and corrected intensity data of two types were compared. The four scenarios include the  According to the median intensity values of the sample data before and after correction in flight strip 1 and 2 in Figures 8-11, the intensity variation of the first type before correction I amp was larger than the intensity variation of the second type I re f , which partly explains the range correction for the intensity of the first type I amp is necessary in addition to incidence angle correction. The intensity of the first type after range and incidence angle correction is a value directly proportional to the spectral reflectance. Although the intensity of the second type has been calibrated by the RIEGL manufacturer, it is still necessary to eliminate the influence of incident angle, obtaining the intensity value proportional to the spectral reflectance.

Comparison of Land Cover Classification Performance Based on Airborne LiDAR Intensity before and after Correction
Using the random forest method, which is one of the best classification methods in machine learning [47], the classification accuracies for land cover based on the uncorrected and corrected intensity data of two types were compared. The four scenarios include the intensity of the first type (calibrated amplitude) before correction Intensity amp , the intensity of the first type (calibrated amplitude) after correction Intensity Camp , the intensity of the second type (reflectance) before correction Intensity re f , and the intensity of the second type (reflectance) after correction Intensity Cre f .
From overall classification accuracies shown in Figure 12, the overall accuracy based on the intensity of the first type before correction Intensity amp and after correction Intensity Camp was 44.10% and 74.11%, respectively. The overall accuracy based on the intensity of the second type before correction Intensity re f and after correction Intensity Cre f was 61.38% and 79.59%. After correction by the improved theoretical model, the overall classification accuracy based on the intensity of first type (Intensity Camp ) was significantly improved by 30.01%, the overall classification accuracy based on the intensity of second type (Intensity Cre f ) increased by 18.21%. The overall classification accuracy based on the intensity of two types after correction was higher than that obtained before correction. The overall classification accuracy based on the uncorrected intensity of the first type Intensity amp was Remote Sens. 2021, 13, 511 14 of 21 less than that the uncorrected intensity of the second type Intensity re f . For both uncorrected and corrected intensity data, the intensity of the second type (Intensity re f and Intensity Cre f ) was more valuable than the intensity of first type (Intensity amp and Intensity Camp ) for land cover classification. was 61.38% and 79.59%. After correction by the improved theoretical model, the overall classification accuracy based on the intensity of first type ( ) was significantly improved by 30.01%, the overall classification accuracy based on the intensity of second type ( ) increased by 18.21%. The overall classification accuracy based on the intensity of two types after correction was higher than that obtained before correction. The overall classification accuracy based on the uncorrected intensity of the first type was less than that the uncorrected intensity of the second type . For both uncorrected and corrected intensity data, the intensity of the second type ( and ) was more valuable than the intensity of first type ( and ) for land cover classification.  Table 4 shows the comparison of classification accuracies for the six land-cover types based on airborne LiDAR intensity before and after correction. The corrected airborne Li-DAR intensity data has good performance in identifying tree types, with the highest CA (78.14% and 80.19%) among all ground objects. The second highest classification accuracy was buildings with a CA of 61.97% and 68.70%, followed by roads with a CV of 48.33% and 58.98%, respectively. Bare soil, grass, and cropland classes have the worst classification accuracy. The reason is that when the laser pulses hit grasses and croplands, the returned signal may come from bare soil. Although not all land cover types can be well classified, the proposed intensity correction method can enhance the separability of the land cover features.  Table 4 shows the comparison of classification accuracies for the six land-cover types based on airborne LiDAR intensity before and after correction. The corrected airborne LiDAR intensity data has good performance in identifying tree types, with the highest CA (78.14% and 80.19%) among all ground objects. The second highest classification accuracy was buildings with a CA of 61.97% and 68.70%, followed by roads with a CV of 48.33% and 58.98%, respectively. Bare soil, grass, and cropland classes have the worst classification accuracy. The reason is that when the laser pulses hit grasses and croplands, the returned signal may come from bare soil. Although not all land cover types can be well classified, the proposed intensity correction method can enhance the separability of the land cover features.
By analyzing the classification accuracies of the uncorrected intensity of second types Intensity re f and the first Intensity amp presented in Figure 13, the classification accuracy of all classes showed a net increase. The highest accuracy improvements (28.21% and 28.94%) in PA and UA were achieved for the roads. Net increases of buildings in PA and UA were 26.21% and 13.20%, whereas the net increase in PA of cropland, grass, bare soil, and tree types was all less than 5% (0.51%, 1.69%, 3.13%, and 2.08%, respectively). The net increase in UA of cropland, grass, bare soil, and tree types were larger than 5% (9.77%, 10.46%, 6.84%, and 12.00%, respectively). Single class of all types with CV showed a net increase, as well as the PA and UA. The road achieved the highest CA of the net increases with 23.43%, followed by building (19.27%), the net increase of tree was 16.54%, the net growth of bare soil, cropland, grass type were significantly lower than the above types, bare soil types was 2.32%, the grass type was 1.89%, the cropland type was 1.50%. One should note that using the intensity of the second type Intensity re f was more beneficial than the intensity of the first type for land-cover classification, especially distinguish road, building, and tree classes. It demonstrated that the Intensity re f data calibrated by RIEGL had a positive impact on the accuracy of land cover classification. building, and tree classes. It demonstrated that the data calibrated by RIEGL had a positive impact on the accuracy of land cover classification.  Figure 14 illustrated net changes in the classification accuracies of single class based on the corrected intensity of the first type vs. the uncorrected intensity of the first type . Net increases in PA and UA were observed for all classes. The classification accuracy for bare soil and grass improved more 40%, the bare soil and grass were barely identified based on the uncorrected intensity of the first type. The net increase in PA and UA of cropland and road was more than 30%. The building and tree showed minor changes in PA and UA, they produced more than 60% in PA. All classes resulted in net gains in CA, especially the bare soil produced the highest net gains (50.39%). The classification accuracy of the intensity after correction was significantly improved compared with the uncorrected intensity of the first type , especially in distinguishing the bare soil and grass. The correction can help improve the classification of some land cover features. Before correction, it was explained that overlap between intensity values of grass and soil and intensity values of soil and cropland. After correction, the intensity  The classification accuracy for bare soil and grass improved more 40%, the bare soil and grass were barely identified based on the uncorrected intensity of the first type. The net increase in PA and UA of cropland and road was more than 30%. The building and tree showed minor changes in PA and UA, they produced more than 60% in PA. All classes resulted in net gains in CA, especially the bare soil produced the highest net gains (50.39%). The classification accuracy of the intensity after correction Intensity Camp was significantly improved compared with the uncorrected intensity of the first type Intensity amp , especially in distinguishing the bare soil and grass. The correction can help improve the classification of some land cover features. Before correction, it was explained that overlap between intensity values of grass and soil and intensity values of soil and cropland. After correction, the intensity values of these land cover features showed some degree of separability.
first type . Net increases in PA and UA were observed for all classes. Th sification accuracy for bare soil and grass improved more 40%, the bare soil and gras barely identified based on the uncorrected intensity of the first type. The net increase and UA of cropland and road was more than 30%. The building and tree showed changes in PA and UA, they produced more than 60% in PA. All classes resulted in ne in CA, especially the bare soil produced the highest net gains (50.39%). The classif accuracy of the intensity after correction was significantly improved pared with the uncorrected intensity of the first type , especially in guishing the bare soil and grass. The correction can help improve the classification o land cover features. Before correction, it was explained that overlap between intensi ues of grass and soil and intensity values of soil and cropland. After correction, the in values of these land cover features showed some degree of separability.  Figure 15 showed the net changes in the classification accuracies of single class based on the corrected intensity of the second type Intensity Cre f vs. the uncorrected intensity of the second type Intensity re f . While PA decreased 2.59% for building, there are net increases in PA and UA for all classes. Bare soil and grass showed the highest net increase (73.44% and 60.16%, 61.91% and 57.45%) in PA and UA. The classification results of bare soil and grass were poor based on the intensity of the second type before correction Intensity re f , which was in line with the conclusion based on the intensity of the first type before correction Intensity amp . All classes resulted in net gains in CA, especially the bare soil and grass produced the highest net gains (54.38% and 49.54%), this is in line with the results from the intensity of first type after correction Intensity Camp vs. before correction Intensity amp . It also indicates that the identification of bare soil and grass can be effectively improved based on the corrected intensity. The net increase in the classification accuracy of cropland was 47.15%, and the corrected intensity was significant in the identification of cropland. The net increase in the classification accuracy of buildings and tree classes was the lowest but greater than 10%. The separability amongst the intensity values of grass, cropland, and soil were significantly increased after correction. Figure 16 shows net changes in the classification accuracies of single class based on the corrected intensity of the second type Intensity Cre f vs. the corrected intensity of the first type Intensity Camp . While Cropland and grass show a net decrease of 0.18% and 0.72% in UA, the net increase in PA were the highest (19.37% and 18.20%), improvement in CA were 12.44% and 12.07%. Road produced net gains in PA(12.50%), 4.42% for UA, 10.65% for CA. Building shows a net increases in UA(9.92%) and CA(6.72%), while building slightly decreased 1.8% in PA. Improvement in PA for bare soil is 9.38%, 1.82% for UA, 6.31% for CA. The PA, UA, and CA of tree also showed a net increase, but only increased by 0.04%, 2.50%, and 2.05%. When identifying tree types, their PA and UA were both greater than 85%, so the intensity of two types after correction in identifying tree types was basically equivalent. A significant improvement in separability between the intensity values of tree and other classes was shown after correction. This can be explained by the fact that vegetation has a strong reflective peak in the near-infrared band. The separability of trees were higher than other vegetation, e.g., grasses and croplands. The reason is that when the laser pulses hit grasses and croplands, the returned signal may partly come from bare soil.
In general, the intensity of the second type is more useful than the intensity of first type for land cover classification.
soil and grass produced the highest net gains (54.38% and 49.54%), this is in line with the results from the intensity of first type after correction vs. before correction . It also indicates that the identification of bare soil and grass can be effectively improved based on the corrected intensity. The net increase in the classification accuracy of cropland was 47.15%, and the corrected intensity was significant in the identification of cropland. The net increase in the classification accuracy of buildings and tree classes was the lowest but greater than 10%. The separability amongst the intensity values of grass, cropland, and soil were significantly increased after correction.  . While Cropland and grass show a net decrease of 0.18% and 0.72% in UA, the net increase in PA were the highest (19.37% and 18.20%), improvement in CA were 12.44% and 12.07%. Road produced net gains in PA(12.50%), 4.42% for UA, 10.65% for CA. Building shows a net increases in UA(9.92%) and CA(6.72%), while building slightly decreased 1.8% in PA. Improvement in PA for bare soil is 9.38%, 1.82% for UA, 6.31% for CA. The PA, UA, and CA of tree also showed a net increase, but only increased by 0.04%, 2.50%, and 2.05%. When identifying tree types, their PA and UA were both greater than 85%, so the intensity of two types after correction in identifying tree types was basically equivalent. A significant improvement in separability between the intensity values of tree and other classes was shown after correction. This can be explained by the fact that vegetation has a strong reflective peak in the near-infrared band. The separability of trees were higher than other vegetation, e.g., grasses and croplands. The reason is that when the laser pulses hit grasses and croplands, the returned signal may partly come from bare soil. In general, the intensity of the second type is more useful than the intensity of first type for land cover classification.

Discussion
As can be seen from the above results, the intensity of the two types before and after correction were compared and assessed. The corrected intensity data have been proven to be effective by sample statistical analysis and classification accuracy assessment. The CV and ratio of coefficient variation of the two corrected intensity are smaller. This finding is consistent with Höfle and Pfeifer [25], where the variation coefficient of the intensity value by using the theoretical model in the homogeneous region is lower than that before correction. As shown in Figures 7-10, when the incident angle is greater than 15 degrees, the influence of the incident angle on the intensity value is significantly enhanced. Kukko et al. [48] found that the effect of incident angle on the intensity was obvious, which is the consistent with this study.
To evaluate the performance of the proposed method, the method is compared with the method proposed by Yan and Shaker [13]. In the method proposed by Yan and Shaker [13], a surface slope threshold was used to determine whether the scan angle or the incident angle should be used for intensity correction. When the slope was less than or equivalent to 40°, the incidence angle was used in the radar (range) equation; when the slope exceeded 40°, the scan angle was used instead. Table 5 shown comparisons of classification accuracies for six land cover types obtained from the proposed method (Proposed correction I and Proposed correction II) and the method of Yan and Shaker [13] (Method I and Method II) during the same region. The incident angle is referenced and calculated by Yan and Shaker [13] in Method Ⅰ and Ⅱ. Then, in Method I, we corrected the intensity

Discussion
As can be seen from the above results, the intensity of the two types before and after correction were compared and assessed. The corrected intensity data have been proven to be effective by sample statistical analysis and classification accuracy assessment. The CV and ratio of coefficient variation of the two corrected intensity are smaller. This finding is consistent with Höfle and Pfeifer [25], where the variation coefficient of the intensity value by using the theoretical model in the homogeneous region is lower than that before correction. As shown in Figures 7-10, when the incident angle is greater than 15 degrees, the influence of the incident angle on the intensity value is significantly enhanced. Kukko et al. [48] found that the effect of incident angle on the intensity was obvious, which is the consistent with this study.
To evaluate the performance of the proposed method, the method is compared with the method proposed by Yan and Shaker [13]. In the method proposed by Yan and Shaker [13], a surface slope threshold was used to determine whether the scan angle or the incident angle should be used for intensity correction. When the slope was less than or equivalent to 40 • , the incidence angle was used in the radar (range) equation; when the slope exceeded 40 • , the scan angle was used instead. Table 5 shown comparisons of classification accuracies for six land cover types obtained from the proposed method (Proposed correction I and Proposed correction II) and the method of Yan and Shaker [13] (Method I and Method II) during the same region. The incident angle is referenced and calculated by Yan and Shaker [13] in Method I and II. Then, in Method I, we corrected the intensity of first form (the calibrated amplitude) by the Equation (15). It eliminated the effects of incident angle and range. In Method II, we corrected the intensity data of the second form (the relative reflectance) by the Equation (16). It only eliminated the effects of incident angle. The classification accuracy of the results obtained by the proposed method were almost the highest. The results suggest that the proposed method had better performance than the method of Yan and Shaker [13] in this area. The recognition of building and tree were superior to other types in all methods. The reason is that when the laser pulses hit low vegetation (grass and cropland), the returned signal may come from bare soil. The proposed intensity correction method in identifying building and tree achieved better results than Method I and II. Cropland, grass, and soil were the worst of classification accuracy in Method I and II, which were almost unrecognized. These may be explained by considering that the accuracy of the incidence angle derived on from the TIN model is not guaranteed and limited. Moreover, slope threshold setting is highly dependent on the terrain and cover of the study area, which may affect the correction results and the separability of the land cover features. The intensity values after correction were also more useful than the intensity before correction for land cover classification. Although the intensity of reflectance has been calibrated by the RIEGL manufacturer, our study suggests that it is still necessary to eliminate the influence of incident angle using intensity for identification and classification. The overall classification accuracy of the corrected intensity of reflectance was improved by 18.21%. Calders et al. [39] evaluated the intensity of RIEGL VZ-400, and concluded that the reflectance data provided by RIEGL VZ-400 has certain errors in users' analysis that need radiation correction. The classification accuracy based on corrected intensity has been improved in identifying six land-cover types, especially tree types, with the highest CA (78.14% and 80.19%) among all ground objects. This is different from the conclusions of the study by Yan and Shaker [13] in which intensity correction was not applicable for tree canopies.
However, the corrected intensity values are only proportional to the ground reflectance and do not reflect the true reflectance values. The radar (range) equation was simplified under the assumption that the target was an extended Lambertian. The effect of incidence angle was corrected by using Lambert's cosines law. In future research, we will consider a realistic illumination model by taking into account the superposition of specular reflected light, diffuse reflected light and ambient light to eliminate the incident effect and improve the actual precision of intensity values. We will perform rigorous radiometric correction and calibration to extract the reflectance values, using the brightness measurement.
From the transmitted signal to the received signal, the system, environment, and target effects may lead to errors of intensity values. In addition to the influence of the incidence angle and scanning range, other parameters such as transmission power and AGC and transmission errors of the system may affect the intensity. However, we assumed that instrument parameters were stable, without considering the influence of the instrument system parameters, therefore instrument effects need to be further analyzed. Because the aircraft flew at a low altitude of 300 m, we assumed that the effect of atmosphere attenuation on LiDAR intensity can be ignored, but, in fact, the laser energy in the process of transmission is affected by temperature, humidity, and air pressure. The atmospheric attenuation need to be corrected in future research.

Conclusions
We presented a theoretical correction using the proposed incidence angle to solve the problem of inconsistent in the overlapping region of airborne LiDAR intensity data. The corrected intensity is directly proportional to the ground reflectance for land cover classification. The intensity of two types (amplitude and reflectance information) from the RIEGL scanner were corrected, and the effectiveness of the improved method was evaluated by the sample statistical method and classification accuracy assessment. Compared with the intensity values before correction, the intensity values after correction are more concentrated and the coefficient of variation is smaller, which proves that the method is effective. It is found that when the incident angle is greater than 15 degrees, the influence of the incident angle on the intensity value is significantly enhanced. It also indicates that the effect of incidence angle cannot be ignored in the intensity correction of airborne LiDAR.
The classification accuracies for land cover based on the uncorrected and corrected intensity data of two types were compared. Six classes are differentiated: building, pavement/road, tree, cropland, grass, and bare soil. In the case of uncorrected data, the intensity of the second type (reflectance) is more valuable than the intensity of first type (amplitude) for land cover classification, especially for distinguishing road, building, and tree classes. After correction the overall accuracy of the intensity of the first type (amplitude) is significantly improved by 30.01%, the overall accuracy of the intensity of second form (reflectance) is increased by 18.21%. The classification results of bare soil and grass were poor based on the intensity of the two type before correction. Distinguishing bare soil and grass after correction show large improvements. Moreover, the corrected intensity data has good performance in identifying trees. In the case of correction, the intensity of the second type (reflectance) is more useful than the intensity of first type (amplitude) for land cover classification. We recommend that radiometric correction of airborne LiDAR intensity data be conducted to improve classification accuracy. This correction approach not only avoids manual intervention, but also has theoretical bases. The study is significant for improving the application of multispectral and hyper-spectral LiDAR intensity data.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to collecting data by the team and the partner, the partner demanded strictly not sharing data.