Point Cloud Intensity Correction for 2D LiDAR Mobile Laser Scanning

,


Introduction
Light detection and ranging (LiDAR), a type of noncontact active remote sensing sensor, can quickly scan highresolution and high-precision 3D point cloud data on the target surface. LiDAR has been successfully applied in fields, such as global climate change research, smart cities, forest resource surveys, environmental monitoring, and basic mapping [1,2]. LiDAR can not only obtain the 3D geometric information of the target surface, but it can also record the intensity information. The intensity corresponds to the coordinate information one-to-one without registration and has the feature of pixel-level fusion. It represents the reflection spectral characteristics of the object target to the laser and can be used as an important feature of target classification [3][4][5]. However, due to the influence of various factors, such as scanner characteristics, atmospheric transmission characteristics, target surface parameters, and data acquisition parameters, there is a large deviation between the intensity and the spectral reflection characteristics of the object target. The phenomena of the same object target with different spectra and the same spectrum of different object targets are apt to occur. Therefore, it is necessary to eliminate the influence of various factors through correction [6][7][8].
Intensity correction methods can be divided into theoretical correction and empirical correction [9]. The theoretical correction method focuses on the analysis of the relationship between multiple factors causing the intensity change, and the regression model of each influencing factor and intensity is established according to the LiDAR ranging principle. Song et al. [10] systematically discussed and analyzed the radiation characteristics, as well as the key influencing factors of intensity data from the perspective of the LiDAR radiation transmission mechanism. They also eliminated these factors through a theoretical correction model. Bolkas [11] studied the intensity correction method of the incident angle and evaluated the intensity correction effect of four kinds of light reflection models under different colors and glosses. Cheng et al. [12] proposed a method to eliminate the distance effect of laser intensity by using a piecewise polynomial model, which effectively eliminated the intensity deviation caused by distance. Fang et al. [13], based on the principle of laser imaging, deduced an intensity theoretical model that considered the defocus effect of the receiving optical system and applied the model to the intensity correction of fresco. The theoretical correction model is applicable to a wide range of scenarios, but the process of the established model is more complex. The parameters in the theoretical correction model are always set as the measured value under the ideal state, which produces large errors in the actual intensity correction. For example, the measurement of the atmospheric attenuation coefficient and transmission coefficient of a LiDAR optical system usually has fluctuation errors. The reflection characteristics of the scanned object's surface are very complex, and it is difficult to regard it as an ideal diffuse reflection target.
The empirical correction method does not rely on any theoretical model. However, it only establishes the intensity correction model in the form of elementary functions (such as polynomial, cosine, or exponential functions). It estimates the function parameters from the actual measured data, which is suitable for situations where the structure and parameters of LiDAR are poorly known. Li and Cheng [14] established two data-driven models to correct the effects of the distance and incident angle on intensity and then applied the corrected intensity to the damage detection of historic buildings. Coren and Sterzai [15] adopted the empirical correction method to establish an exponential function model between the intensity and atmospheric attenuation coefficient, which reduced the fluctuation of the intensity data of asphalt pavement. Vain et al. [16] and Korpela et al. [17] used the empirical correction method to compensate for the intensity change caused by the LiDAR internal automatic gain control system. For the intensity data of specific scenes, the empirical correction method has a higher accuracy than the theoretical correction method, but it is only applicable to this specific application scene; otherwise, it needs to be remodeled.
LiDAR measurements are often classified into three types: terrestrial laser scanning (TLS), mobile laser scanning (MLS), and airborne laser scanning (ALS). TLS can directly obtain high-precision 3D point cloud data from the surrounding natural environment by setting up fixed sites. MLS, on the other hand, can quickly scan 3D point cloud data on both sides of the road or on one side by carrying LiDAR on mobile devices, such as cars [18]. Compared with TLS, MLS is flexible and expands the scanning range of LiDAR. MLS also has mature applications in digital city 3D modeling, urban environment monitoring, and urban resource surveys [19,20]. ALS can obtain the surface spatial information within a large area in a short time with high working efficiency by carrying LiDAR onto the aircraft to realize scanning. However, compared with MLS, point cloud data obtained via ALS measurements are less accurate. Moreover, MLS has an advantage over ALS in obtaining vertical point cloud data. Yang et al. [21] showed that, compared with ALS, MLS has considerable advantages in street tree identification at the scale of a single tree, as well as a strong data acquisition ability to penetrate the inner canopy and trunk.
LiDAR sensors are usually divided into two kinds: 3D and 2D (also known as multithread and single-thread sensors). The 3D LiDAR can accurately obtain 3D geometric information of the surrounding environment, and it is mainly used in the field of unmanned driving, which requires high accuracy. However, 3D LiDAR is usually expensive. In contrast, 2D LiDAR can obtain highprecision 3D information about the environment throughout the process of moving through fan-shaped scanning frame by frame, which is advantageous because of the fast scanning speed, the small size, low power consumption, and low manufacturing cost. In addition, 2D LiDAR's acquisition of point cloud data has low redundancy and simple data fusion. Thus, it can be directly indexed according to the frame number and in-frame number of the measured points, avoiding accuracy loss caused by grid partitioning [22]. 2D LiDAR using the MLS measurement method is widely used in the field of tree parameter extraction [23,24] and urban block ground object classification [25,26]. The point cloud data collection process by 2D LiDAR is shown in Figure 1. Xu et al. [27] and Nan et al. [28] designed a real-time automatic target spray system based on MLS 2D LiDAR detection and compared the performance with that of infrared, ultrasonic sensors, and 3D LiDAR, proving that 2D LiDAR has great advantages in the identification accuracy and rapid measurement of distance.
To the author's knowledge, most previous intensity correction studies have focused on 3D LiDAR sensors for ALS and TLS, whereas 2D LiDAR sensors for MLS are still relatively rare. Based on the advantages of MLS over ALS and TLS analyzed above in application value, the important contribution of this work proposes a high-precision intensity correction method suitable for 2D MLS. In recent years, there have been some studies on intensity correction for TLS similar to this paper. Tan and Cheng [29] studied TLS intensity correction based on the laser ranging formula. The relationship between intensity and new variables was established by combining the cosine of the incident angle and the square of the distance. In the intensity correction for different areas of the white wall, the intensity varianceto-mean ratio ε was between 0.6 and 0.78, indicating that the improvements in intensity consistency ranged from 22% to 40%. Subsequently, Tan and Cheng [30] investigated the effects of incidence angle and distance on intensity data and corrected the intensity data of Faro Focus3D 120. For four reference targets with reflectances of 20%, 40%, 60%, and 80%, the linear interpolation method was used to fit the relationship between incidence angle versus intensity and distance versus intensity. A total of 20 small regions with a size of approximately 15 cm × 15 cm in the white lime wall surface were randomly sampled to verify the intensity correction effect. The intensity variance-to-mean ratio ε  [31] proposed a distance intensity correction workflow for MLS road point clouds using a data-driven approach. The relationship between distance and intensity was fitted by a piecewise polynomial. By comparing the differences between the scanners and the stripe intensities before and after correction, the improvements in intensity consistency ranged from 47% to 56%. In contrast from the general data-driven intensity correction model, which directly uses various models to establish the approximate relationship between intensity and its influencing factors, lacking the necessary theoretical basis, a new point cloud intensity correction method for 2D MLS is proposed by combining theoretical derivation with empirical correction. In addition, according to the scan pattern, a 2D scan grid is constructed to organize the MLS intensity of the wall, and a new method of a spherical neighborhood search fitting plane is proposed to accurately calculate the cosine of the incident angle. In particular, the relationship between the neighborhood radius and the distance is discussed in the process of plane fitting. The accuracy of incident angle measurement is improved by selecting an appropriate neighborhood radius. Two groups of verification experiments are carried out on single sites and multiple sites to test the effect of the intensity correction model. A single-site experiment shows that the ε value of the five same area regions within the range of the distance and incident angle [0.52 m-1.55 m, 0°-74°] is approximately 0.3, indicating that the consistency of intensity has been improved by 70% after correction. Multisite experiments concluded that the ε values of sites A, B, C, and D were 0.073, 0.079, 0.233, and 0.280, respectively, within the range of distance and incident angle [1.52 m-5.34 m, 0°-62°]. This means that the improvements in intensity consistency range from 72% to 92.7% after correction. Overall, this approach is superior to the latest study mentioned above.
The proposed method of the intensity correction model is introduced in Section 2. The process of obtaining the intensity correction model is described in Section 3. To demonstrate the practical effectiveness of the intensity correction model, two different scenes of wall intensity correction are presented and analyzed in Section 4. Section 5 concludes with the findings. Figure 2 is the flow chart of intensity correction for the MLS 2D LiDAR point cloud. In this paper, the intensity correction process consists of two steps. The first step is for model establishment, and the second step is for model validation. Model establishment is when the intensity correction model is obtained. First, the influencing factors of MLS intensity (including target reflectivity, distance, and incident angle) are analyzed for the acquisition of the intensity correction model. Second, the intensity multiplicative model is established, and the intensity correction formula for the distance and incident angle is deduced. Finally, the parameters of the intensity correction model are estimated, including the order and coefficient of the polynomial. Model validation uses the intensity correction model to correct the intensity data of the wall point cloud. For the wall point cloud intensity data scanned by MLS 2D LiDAR, the distance intensity and incident angle intensity are extracted. Then, the obtained intensity correction model is used to correct the intensity data of the wall point cloud.

Influencing Factors of MLS Intensity.
Assuming that the measured target surface is a Lambert surface (the surface with ideal diffuse reflection characteristics), according to the laser ranging formula [32], the echo power received by the laser detector is

Wireless Communications and Mobile Computing
where the received laser power P r is a function of the transmitting laser power P t , the aperture of the receiving detector is denoted by D r , the atmospheric one-way extinction coefficient is η atm , the transmission coefficient of the optical system is η sys , the target reflectivity is λ, the distance R is between the measuring point and the receiving end of the LiDAR laser (referred to as the distance), and the cosine of the laser incident angle is denoted by cos θ. For the MLS system with a short distance, it can be considered that P t , D r , η atm , and η sys are constants. Tan and Cheng [29] showed that intensity is the laser echo energy represented digitally. Additionally, there is a certain linear relationship between the received laser power and intensity. According to the above analysis, the main influencing factors of MLS intensity are target reflectivity λ, distance R, and incident angle θ.
2.3. The Multiplicative Model of Intensity by MLS. The objective of intensity correction in this paper is to remove the influence of the distance and incident angle on MLS intensity. Then, the intensity value is only related to the reflectivity of the target. Since the received laser power is nonlinearly processed internally by LiDAR, it is not possible to directly use Formula (1) for theoretical correction. In addition, considering that the effects of the reflectivity, distance, and incident angle are independent in theory, they can each be solved separately. The multiplicative model of intensity can be established as Iðλ, R, θÞ is the original intensity, and f λ ðλÞ, f R ðRÞ, and f θ ðcos θÞ represent functions of the independent influence of the target reflectivity, distance, and incident angle as Iðλ , R, θÞ. To eliminate the influence of the distance and incident angle on intensity, Iðλ, R, θÞ at any distance and incident angle should be transformed to the corrected intensity I c at reference distance R 0 and reference incident angle θ 0 for a Lambert body with the same target reflectivity.
Formula (3) shows that the transformation relationship between I c and Iðλ, R, θÞ is established by f R ðRÞ and f θ ðcos θÞ.
After satisfying the conditions of a target with reflectivity, λ is scanned at different distances R and at a constant incident angle θ 0 ; then, the intensity after the distance correction I rc can be written as Under the above premise conditions, R min ≤ R ≤ R max (R min and R max are the minimum value and maximum value of the distance measuring range, respectively). According to the Weierstrass theorem [33], a continuous function on a closed interval can be uniformly approximated by a polynomial series. According to the definition, f R ðRÞ describes the relationship between distance R and Iðλ, R, θÞ in the distance-intensity data. The intensity regression value Iðλ, R, θ 0 Þ under reflectivity λ and reference incident angle θ 0 can be obtained by polynomial regression fitting. Tan and Cheng [30] showed that due to the short-distance effect of the LiDAR optical system, the intensity increases with an increasing distance when the distance is relatively close and decreases with an increasing distance when the distance is relatively far. According to the variable of distance R, piecewise polynomial modeling can be adopted here to illustrate the relationship between R and Iðλ, R, θ 0 Þ as where a k and b l are the coefficients of the distance polynomial, K and L are the order of the distance polynomial, and    (4), I rc can be expressed as follows: Similarly, satisfying the conditions of a reference target with reflectivity, λ is scanned at different incident angles θ and at a constant distance R 0 , deduced from Formula (3). The intensity after the incident angle correction I θc can be written as In the above incident angle-intensity data, 0 ≤ cos θ ≤ 1. According to the Weierstrass theorem mentioned above, the intensity regression value Iðλ, R 0 , θÞ under the reference reflectivity λ 0 and the reference incident angle θ 0 can also be obtained by polynomial regression fitting. A cosine polynomial is used for modeling f θ ðcos θÞ, as shown in the following formula: where c m is the coefficient of the incident angle cosine polynomial and M is the order of the incident angle cosine polynomial. By substituting the obtained f θ ðcos θÞ into Formula (7), I θc is established as follows:

Parameter Estimation of the MLS Intensity Correction
Model. In this paper, the elbow rule method and the least square method are used to determine the polynomial order and the fitting polynomial coefficients. For example, in the polynomial in the short-distance segment (R ≤ R t ) of Formula (5), the values of K and a k need to be determined in a way that enables the error sum of squares Sða 0 , a 1 ,⋯a K Þ between Iðλ, R, θ 0 Þ and Iðλ, R, θÞ to be minimized. Sða 0 , a 1 ,⋯,a K Þ is the cost function. The root mean square error (RMSE) of the intensity sample is defined as where N is the number of Iðλ, R, θÞ and Iðλ, R, θ 0 Þ and n is the serial number of the intensity sample. The smaller the RMSE value is, the lower the fitting error is. In general, the higher the polynomial order (the larger the K value), the smaller the RMSE value. However, it should be noted that if the fitting order K value is selected too high, it is easy to overfit when the distance correction model is applied to the actual scene for intensity correction and the correction effect is poor. If the K value is too low, the distance intensity correction cannot be carried out effectively. To avoid overfitting, according to the change curve of different polynomial order RMSE values, the elbow rule in machine learning is used to determine that the elbow position is the best polynomial order.
To obtain the polynomial coefficient a k under Sða 0 , a 1 , ⋯a K Þ the minimum value, the polynomial in the shortdistance segment (R ≤ R t ) of Formula (5) is brought into Formula (10). Then, the partial derivative of each polynomial coefficient a k is transformed into the problem of finding the extremum; thus,

Wireless Communications and Mobile Computing
It can be further obtained that The linear formulas can be obtained as follows: : ð13Þ In this way, the polynomial coefficients ½a 0 , a 1 ,⋯,a k of f R ðRÞ can be solved by the Gaussian elimination method. Similarly, bringing the polynomial in the long-distance segment (R > R t ) of Formula (5) and the cosine polynomial of Formula (8) into Formula (10), based on the elbow method and the least square method, the corresponding polynomial coefficient parameters ½b 0 , b 1 ,⋯,b l and ½c 0 , c 1 ,⋯c m are obtained.

Evaluation Indices of the MLS Intensity Correction
Model. To evaluate the effect of distance and incident angle correction models on MLS intensity correction, the coefficient of variation CV is used to represent the degree of discretization before and after intensity correction. The coefficient of variation CV is defined by where STD is the standard deviation of intensity and Mean is the mean value of intensity. CV is positively correlated with the dispersion of intensity data. The smaller its value, the smaller the data dispersion. In addition, the evaluation index variance-to-mean ratio ε of the intensity correction model is defined as follows: CV ori and CV cor represent the variation coefficient of the original intensity and the variation coefficient of the corrected intensity, which can be calculated by Formula (14). When ε is less than 1, it means that the variability of the intensity value after correction is less than that before correction, and the intensity correction model is effective. The smaller the value of ε is, the better the intensity consistency of the model after correction. The reference target is a standard diffuse reflectance plate with a size of 50 cm × 50 cm and a reflectance of 50% (that is, the reflectance of the reference target is λ 0 =0.5), as shown in Figure 4. The diffuse reflector plate is composed of highly diffuse reflector materials. When irradiated by a laser beam with a wavelength of 905 nm, the diffuse reflector plate can be regarded as a standard Lambert body.

Experiment and Data Acquisition
The UTM-30LX 2D LiDAR was placed vertically against the diffuse reflectance plate in an indoor environment. Since

Distance
Correction Model Acquisition Experiment. The reference incident angle was set to θ 0 = 0°, and the distance-intensity data ðR i , Iðλ 0 , R i , θ 0 ÞÞ were obtained by setting different distance sites to scan the diffuse reflectance plate. Since the intensity measurement error was large when the distance was long, the intensity data of the diffuse reflectance plate were collected within the range of 0.   Figure 5.
The 2D LiDAR was fixed on the ground and placed horizontally opposite the diffuse reflector plate. To ensure that the distance-intensity data ðR i , Iðλ 0 , R i , θ 0 ÞÞ were obtained at the reference angle θ 0 = 0°, the laser beam at the intermediate point of the LiDAR at each site (frame 541) was perpendicular to the diffuse reflector plate. To improve the measurement accuracy of distance-intensity data, 10 frames were collected continuously at each site, and then, the average intensity value was calculated. Setting the reference incident angle θ 0 = 0°, within the range of 0.1 m-14.4 m, the relationship between different distances and the corresponding intensity is shown in Figure 6.
In the range of 0.1-0.7 m, the intensity increases rapidly with an increasing distance value, and after 0.7 m, the intensity decreases slowly with an increasing distance value. Therefore, the polynomial function can be fitted by Formula (5) with 0.7 m as the distance cutoff point.

Incident Angle Correction Model Acquisition Experiment.
The reference distance was set to R 0 = 1:2 m, and the intensities of different incident angles within the range of 0°-80°w ere collected. To ensure that the incident angle-intensity data ðθ i , Iðλ 0 , R 0 , θ i ÞÞ were measured at R 0 = 1:2 m, the laser beam at frame 541 of LiDAR was positioned directly against the central axis of the diffuse reflector plate. The position of LiDAR was fixed, and the diffuse reflector plate central axis was taken to change the position within the range of 0°-80°a t intervals of 10°. The experiment is shown in Figure 7. Ten frames were collected by continuous scanning at each angle site (frame 541), the incident angle-intensity data of the laser beam position were obtained, and the average intensity value was calculated. The experimental result is shown in Figure 8.

Determining the Order and Coefficient of the Polynomial
Model. The least square method described in Section 2.4 is used to determine the order and coefficient of the polynomial model of distance and incident angle. To avoid overfitting the intensity data, according to the elbow rule in machine learning, the best order of the polynomial is determined to be the elbow position. To determine the ideal polynomial order, the distance correction model acquisition experiment in Section 3.2 and the incident angle correction model acquisition experiment in Section 3.3 were repeated 5 times to expand the dataset. Table 1 lists RMSE for different polynomial order K, L, and M values. The process of determining the optimal values of K, L, and M based on RMSE is shown in Figure 9. It is proven that K = 4 and L = 3 are the ideal distance piecewise polynomial orders in the short-distance and long-distance segments, respectively,   and M = 1 is the ideal incident angle cosine polynomial order. As shown in Figure 9(a), the RMSE decreases substantially with increasing K and L values in the distance correction experiment, and K ≥ 4 and L ≥ 3 tend to be flat. According to the principle of the elbow rule in machine learning, the order of the distance polynomial model is set to K = 4 and L = 3. The blue and orange circles in Figure 9(a) represent the optimal values of K and L, respectively. In Figure 9(b), compared with K and L, the different RMSE values and fluctuation ranges corresponding to the order M of the incident angle cosine polynomial are relatively small. To avoid overfitting, the order of the incident angle polynomial model is set to M = 1. The gray circle in Figure 9(b) represents the optimal value of M.
The function expression in Table 2 is only applicable to the point cloud intensity data obtained by UTM-30LX 2D LiDAR scanning the same target reflectivity Lambert body and can only fit the intensity data within the range of 0.1 m-14.4 m and incident angle value of 0°-80°. The fitting results of the piecewise distance and incident angle cosine polynomial are shown in Figure 10.
R 0 =1.2 m and θ 0 = 0°are taken and substituted by the distance-intensity data ðR i , Iðλ 0 , R i , θ 0 ÞÞ and the incident angle-intensity data ðθ i , Iðλ 0 , R 0 , θ i ÞÞ mentioned above into Formulas (6) and (9) for intensity correction. The intensity after correction at different sites is shown in Figure 11. The intensity distribution of each site is approximately 3400, although there are some model errors.

MLS 2D LiDAR Point Cloud Data Acquisition System.
To verify the validity of the intensity correction model for the distance and incident angle, a flat wall was selected as the experimental object, which was roughly considered a Lambert surface. To visually observe the intensity point cloud, the intensity value was converted into RGB in Cloud-Compare, a software point cloud development tool. The 3D point cloud RGB intensity map was also converted into a 2D pseudocolor map for display. The actual scene diagram of the experimental scanning wall is shown in Figure 12. The MLS 2D LiDAR measurement system emitted laser pulses in all directions through internal rotating optical components to form a 2D fan-shaped scanning surface. The moving platform carried the laser pulses along the direction perpendicular to the scanning surface to realize the 3D     Figure 13.
The x-axis is the moving direction of LiDAR, the y-axis is the scanning direction of LiDAR, the z-axis is vertical to the ground, αði, jÞ is the scanning angle of the i th measuring point in the j th frame, and R is the distance between the LiDAR laser receiver and the measuring where i is the in-frame number, j is the frame number of the measured points, and xði, jÞ, yði, jÞ, and zði, jÞ represent the X, Y, and Z coordinates of the i th measuring point in the j th frame. v is the speed of the moving slide platform in the direction of the motion, and Δt is the scan cycle of the 2D LiDAR. The geometric relationship between the moving slide with 2D LiDAR and the wall is shown in Figure 14. The coordinate of LiDAR center point O at different scanning moments was set as ðxði, jÞ, 0, 0Þ; the 3D coordinate of the i th measurement point in frame j th was pði, jÞ =(xði, jÞ ,yði, jÞ,zði, jÞ), incident laser vector lði, jÞ = ð0, yði, jÞ, zði, jÞÞ, and the normal vector of scanning point nði, jÞ = ðn1, n2, n3Þ.

4.2.
Obtaining R and cos θ. Generally, the density of point cloud data obtained by 3D LiDAR is relatively high. Usually, to reduce data redundancy, a support vector machine (SVM) is used to integrate different feature weights of the point cloud into the classifier for training [35][36][37][38]. In addition, to improve the efficiency of neighborhood searching, the k-D tree algorithm is usually used for downsampling 3D LiDAR    [39]. Compared with 3D LiDAR point clouds, the data structure of 2D LiDAR point clouds is relatively simple. It only needs to be indexed by the grid in the established coordinate system, and the k-nearest neighbor algorithm [40] can be used to establish the point cloud neighborhood set. Referring to our previous research [41], a 3D spherical neighborhood Sði, jÞ is defined as a set of the nearest adjacent points within a sphere centered at the measurement point pði, jÞ with a radius of δ. According to the idea of the least square method, the plane is fitted in the neighborhood set Sði, jÞ, and the normal vector nði, jÞ corresponding to each measuring point pði, jÞ is obtained. Then, cos θði, jÞ of the LiDAR laser receiver to each measuring point can be calculated according to Formula (17). Obviously, Rði, jÞ can be obtained directly.
Since lði, jÞ does not change, the accuracy of cos θði, jÞ depends on nði, jÞ. In general, the smaller the value of δ, the smaller the fitting plane, and the more accurate the value of nði, jÞ. However, if the value of δ is too small, the point cloud data in the two adjacent frames cannot be included in the neighborhood; then, the plane cannot be fitted either. In the MLS 2D LiDAR measurement system, the value of δ is closely related to distance Rði, jÞ and the resolution of αði, jÞ. The resolution of the moving direction and the scanning direction of LiDAR is defined as Δx and Δs, respectively: Δα is the 2D LiDAR angular resolution, and its value is constant. The larger the Rði, jÞ is, the greater the value of Δs. Therefore, if the value of δ is less than Δs, in the LiDAR scanning direction, all of the in-frame data points adjacent to the measurement point pði, jÞ cannot be included in Sði, jÞ. Then, nði, jÞ cannot be obtained by plane fitting. According to Formula (18), the necessary conditions for the measurement point pði, jÞ to have As for UTM-30LX, the value of Δα is 0.25°. In the wall intensity correction experiment, the maximum distance of LiDAR from the wall is not more than 5.5 m. According to Formula (19), Rði, jÞ max sinΔα is 0.024. Considering that the range resolution of the LiDAR sensor within 10 m is ±30 mm, the value of δ is set to 0.03 m.

Correction of the Wall Point Cloud Intensity Data.
According to the calculation method of the incident angle and distance, the incident angle θ and the distance value R of each measurement point in the wall point cloud intensity data were calculated. The study of Tan and Cheng [42] shows that since the impact of the incident angle and distance on intensity is independent, the sequence of the incident angle correction and distance correction does not affect the correction results. According to the established intensity correction model, the reference incident angle and reference distance are set to be θ 0 = 0°and R 0 =1.2 m, respectively. The incident angle and distance are corrected for the wall point cloud intensity data. Figure 15 shows the intensity 2D pseudocolor map and intensity distribution histogram of the red rectangular wall area described in Figure 12 before and after intensity correction.
As seen from the RGB diagram of the original intensity in Figure 15(a), the intensity value at different distances  Figure 17: Intensity box diagram of different areas on the same wall before and after correction: (a) intensity box diagram of area S1-S5 before correction; (b) intensity box diagram of area S1-S5 after correction. and incident angles is considerably different before correction under the same target reflectivity. The closer the distance or the smaller the incident angle, the higher the intensity value will be. In contrast, the farther the distance or the larger the incident angle, the lower the intensity value will be. The RGB of the intensity after the correction of the incident angle in Figure 15(b) shows that the uniformity of the wall point cloud intensity data has been improved. The RGB diagram of the distance correction intensity after the incident angle correction is shown in Figure 15(c). The influence of the distance and incident angle on intensity is eliminated, and the corrected intensity value is consistent and similar. The histogram of the intensity distribution in Figure 15(d) shows that the original intensity data of the wall point cloud under the same target reflectivity have no obvious distribution pattern. From Figure 15(f), it is clear that after the correction of the incident angle and distance, the intensity value of the wall is concentrated and presents a Gaussian distribution.

Using an Evaluation
Index ε to Verify the Intensity Correction Model 4.4.1. Single-Site Multiregion Verification Experiment. In the wall point cloud intensity data obtained from the above experiment, 5 different areas of 20 cm × 20 cm were randomly selected, as shown in Figure 16. According to the established intensity correction model, the intensity data in different areas of S1-S5 were also corrected with reference distance R 0 =1.2 m and incident angle θ 0 = 0°. In addition, the intensity value before and after correction was statistically analyzed. The intensity box chart is shown in Figure 17. The distribution of intensity maximum value (MAX), minimum value (MIN), mean, STD, CV, and the evaluation index ε before and after correction is shown in Table 3. As shown in Figure 17(a), the intensity value of different regions is substantially different and fluctuates in a large range within each region before intensity correction. From Figure 17(b), it is evident that the difference in the intensity value between different regions decreases, and the fluctuation range of the intensity within the region also becomes considerably smaller after the correction of the incident angle and distance. It is clear that the consistency of intensity between different areas and within different areas has been substantially improved by comparing the box graphs of intensity data from different areas of walls S1-S5 before and after intensity correction under the same reflectivity of the target. Table 3 also shows a conclusion similar to that in Figure 17. Due to the small area of the S1-S5 region, the dif-ference between the internal distance and incident angle of each region is not large. The STD values of the original intensity data are all below 21. In addition, due to the large difference in distance and incident angle value between different regions, their STD values are also relatively large, with a minimum value of 15.9 and a maximum value of 20.44. After intensity correction, the STD values of the intensity data in each region are all below 10, and the difference between them is small, indicating that the uniformity of the intensity distribution in each region has been substantially improved. The ε value fluctuates approximately 0.3 in S1-S5 different regions, indicating that the intensity consistency of the five same area regions with the range of distance and incident angle [0.52 m-1.55 m, 0°-74°] has been massively improved by 70% after applying the proposed correction method.   Similar to the single-site multiregion verification experiment in this paper, Tan and Cheng [29] utilized Faro Focus3D TLS to obtain intensity point cloud data of ordinary white wall surfaces. In the range of distance and incident angle [6.70 m-14.76 m, 0°-80°], intensity correction experiments were carried out for the A-F region with the same area. The ε value corresponding to the intensity correction of the white wall obtained by the Tan method and the method proposed in this paper is shown in Table 4 below.
The ε value of the Tan method has a large fluctuation range, approximately 0.6, while the ε value of the method proposed in this paper has a small fluctuation range, approximately 0.3. In the Tan experiment, especially in regions E and F, the ε value is greater than 1, indicating that the consistency of the corrected intensity is not as good as the original intensity, and there is an overfitting phenomenon. Moreover, the Tan method only conducts intensity correction research at distances where the intensity value is relatively stable, ignoring the LiDAR short-distance effect, and does not correct the intensity at distances less than 1 m. Overall, we can say that the proposed method in this paper is better than the Tan method and provides a higher accuracy of intensity correction.

Multisite
Scanning of the Whole Wall for the Verification Experiment. In the previous single-site experiment, the area of the S1-S5 region was small, and the inci-dent angle and distance value of the point clouds in each region were not substantially different. The STD of the original intensity was also relatively small and could not effectively reflect the distribution characteristics of the original intensity data of the whole wall. A multisite experiment was devoted to the intensity correction of the same red rectangular area on the wall, as shown in Figure 12, under multiple sites. Four sites were set up: A, B, C, and D. The vertical distances between the sliding platform and the wall were 1.5 m, 2.5 m, 3.5 m, and 4.5 m, respectively. The rectangular area of the wall was scanned to obtain the point cloud intensity data, as shown in Figure 18. According to the established intensity correction model, the intensity data at different sites were corrected at reference distance R 0 =1.2 m and incident angle θ 0 = 0°. Table 5 lists the distance and incident angle range values of different sites A, B, C, and D. Figure 19 shows the intensity pseudocolor map of the red rectangular area on the wall described in Figure 12 before and after intensity correction at sites A, B, C and D. The histogram of the intensity distribution before and after intensity correction is shown in Figure 20.
It is evident from Figures 19(a)-19(d) that as the distance increases, the density of point clouds on walls with the same area decreases, while the intensity value shows great differences within the whole wall of each site. It can be clearly seen from Figures 19(e)-19(h) that the variation of wall intensity of the four sites has been considerably Figure 19: The intensity 2D pseudocolor map of the red rectangular area on the wall described in Figure 12     Meanwhile, the values of CV ori at A, B, C, and D are 0.0797, 0.0419, 0.0270, and 0.0215, respectively. The value of CV ori decreases gradually with increasing distance between the sites and finally stabilizes at approximately 5 m. In combination with the distance and incident angle value in Table 5, it can be seen that the closer the site is to the wall, the larger the range of fluctuation of distance and incident angle in the same area will be, leading to the larger difference of original intensity value. Therefore, the values of CV ori at sites A and B are slightly higher than those at sites C and D. After correction, the values of CV cor at the four sites are 0.0058, 0.0033, 0.0063, and 0.0061, respectively. It has been shown that the intensity variability after correction is  smaller and the intensity correction effect is obvious. Among them, the CV cor values of sites A, C, and D are approximately 0.006, while the CV cor value of site B is only 0.0033, indicating that the corrected intensity consistency of site B is substantially better than that of other sites. The same conclusion can be obtained from Figure 19(f). The reason is that the original intensity distribution of site B is more concentrated than that of sites A, C, and D under the joint influence of distance and incident angle, as shown in Figure 20(b). Consequently, the intensity variability after correction at site B is lower than that at sites A, C, and D. Meanwhile, the ε values of sites A, B, C, and D are 0.073, 0.079, 0.233, and 0.280, respectively. This means that the intensity consistency of the four sites increased by 92.7%, 92.1%, 76.7%, and 72%, respectively. From Figure 19, we can intuitively see that the intensity correction effect of the shortdistance sites A and B is better than that of the longdistance sites C and D. Similar to the research topic in this paper, Tan and Cheng [30] adopted a linear interpolation method to fit the relationship between incident angle versus intensity and distance versus intensity. In the range of distances and incident angles [1 m-29 m, 0°-80°], four reference targets with reflectance of 20%, 40%, 60%, and 80% were established for intensity correction models. Then, a total of 20 small regions with a size of approximately 15 cm × 15 cm in the white lime wall surface were randomly sampled to verify the intensity correction effect over the whole range of distances and incident angles. After using the above four intensity correction models, the intensity variance-to-mean ratio ε was 0.26, 0.14, 0.19, and 0.21, meaning that the intensity consistency was improved by 74%, 86%, 81%, and 79%.
Compared with Tan's work, an intensity correction model based on a reference target with 50% reflectivity in the range of distance and incident angle [0.1-14.4 m, 0°-80°] is proposed. The intensity correction effect in this paper is better than Tan's in eliminating the factors of the distance and incident angle. As mentioned above, the improvement rates of intensity consistency before and after correction for site A and site B are 92.7% and 92.1%, respectively. Obviously, regardless of what kind of reference target reflectance Tan's intensity correction model is based on, the intensity consistency of the corrected white wall at sites A and B in this paper is higher than that of his research work. The reasons are as follows: first, to eliminate the distance factor, distance-intensity data with a distance value less than 1 m are considered in this paper. Therefore, compared with Tan's method, the method presented in this paper performs better in short-distance intensity correction. Then, a new method of spherical neighborhood search fitting plane is proposed to accurately calculate the cosine of the incident angle. In particular, the relationship between neighborhood radius and distance is discussed in the process of plane fitting. The accuracy of incident angle measurement is improved by selecting an appropriate neighborhood radius. Finally, unlike Tan's interpolation method, this paper adopts a piecewise distance polynomial and an incident angle cosine polynomial to fit distance-intensity and incident angle-intensity data. In conclusion, compared with Tan's research, the intensity correction model established in this paper has a higher reliability and accuracy.

Conclusion
In this paper, a new point cloud intensity correction method for 2D MLS based on theoretical derivation and empirical correction is proposed to solve the problem that intensity information cannot be directly used for target recognition. Based on the diffuse reflection Lambert body of the same target reflectance, the intensity correction model of the piecewise distance polynomial and incident angle cosine polynomial is adopted, and the model parameters are calculated by experiments. The effectiveness of the intensity correction method is verified by single-site and multisite experiments on a white wall using MLS 2D LiDAR. The experimental results show that the intensity consistency is substantially improved by 70% to 92.7% after correction within the range of the distance and incident angle [0.52 m-5.34 m, 0°-74°]. Compared with the latest research, the intensity correction model proposed in this paper has a higher fitting accuracy and can effectively eliminate the MLS intensity deviation caused by distance and incident angle.
However, the intensity correction model is only suitable for the specific UTM-30LX 2D LiDAR and applicable to targets similar to the standard Lambert surface. The intensity deviation still exists after the correction, indicating that further research is needed, especially to reduce the model error. In addition, it is found in the incident angle correction experiment that the correction of the incident angle is not thorough in place with a large angle, which requires further study of a more rigorous correction model. Furthermore, considering the long running time of LiDAR, the influence of internal temperature rise on intensity can further improve the accuracy of model correction, which is also a research direction for the future.

Data Availability
No data were used to support this study.

Conflicts of Interest
The authors declare that there is no conflict of interests regarding the publication of this paper.