Terrestrial Laser Scanning Intensity Correction by Piecewise Fitting and Overlap-Driven Adjustment

Terrestrial laser scanning sensors deliver not only three-dimensional geometric information of the scanned objects but also the intensity data of returned laser pulse. Recent studies have demonstrated potential applications of intensity data from Terrestrial Laser Scanning (TLS). However, the distance and incident angle effects distort the TLS raw intensity data. To overcome the distortions, a new intensity correction method by combining the piecewise fitting and overlap-driven adjustment approaches was proposed in this study. The distance effect is eliminated by the piecewise fitting approach. The incident angle effect is eliminated by overlap-driven adjustment using the Oren–Nayar model that employs the surface roughness parameter of the scanned object. The surface roughness parameter at a certain point in an overlapped region of the multi-station scans is estimated by using the raw intensity data from two different stations at the point rather than estimated by averaging the surface roughness at other positions for each kind of object, which eliminates the estimation deviation. Experimental results obtained by using a TLS sensor (Riegl VZ-400i) demonstrate that the proposed method is valid and the deviations of the retrieved reflectance values from those measured by a spectrometer are all less than 3%.


Introduction
During the last two decades, terrestrial laser scanning (TLS), an active remote sensing technique, is one of the most significant means to acquire three-dimensional (3D) point cloud (X, Y, Z) containing high-accuracy and high-density surface topography, which has been widely used in a variety of applications [1,2]. A large number of researchers have focused on data processing algorithms for 3D point clouds, such as the segmentation, classification, and extraction approaches from a mass of point cloud data based on the geometric property of the scanned objects [3][4][5]. However, the TLS sensor not only delivers the 3D geometric information of the scanned objects, but also can receive the power of the light backscattered from the scanned object. The received power is transformed and recorded as a digital number called "intensity", which involves the physical characteristics of the scanned object at that point [6]. Owing to the advantage of intensity that it is insensitive to the ambient light and shadowing, unlike the red-green-blue (RGB) values of the target images, it can be used for point cloud data processing as a significant complement [7].
Recently, many studies have demonstrated the potential of intensity data, which is widely used in many applications, such as cultural heritage [8], identification of different rock and soil layers [9], structural damage detection [10], water content extraction [11], road traffic marking identification [12], glacier and snow distribution [13][14][15]. In previous studies, the TLS intensity data was directly used without any radiometric processing [16,17]. Even without radiometric correction, it was proved that the intensity data was a better selection for classifying the tree species and species mixtures [18]. However, the TLS intensity data are influenced by some factors. Some are related to the TLS sensor itself such as the laser wavelength, laser divergence angle and laser power, which can be collectively called as the system transmission factor and considered as a constant. The others are related to the scanned objects including distance, incident angle and physical characteristics of the scanned surface [8,19], which have significant effects on the intensity data and are variable for different scanned objects. The atmospheric effect is negligible and can be ignored for short-range TLS. Therefore, the system transmission factor, distance effect and incident angle effect on TLS intensity data should be eliminated so that the intensity data is only related to the physical characteristics of the scanned objects.
There are some studies related to TLS intensity correction, which have successfully reduced the distance and incident angle effects [6,7,[20][21][22][23]. These methods can be mainly divided into model-driven [20,21] and data-driven methods [6,21]. The model-driven approach uses a physical mathematic model based on the radar range equation to convert the intensity data to a backscattering cross-section or backscattering coefficient with physical meaning. In contrast, the data-driven approach uses an empirical mathematic model derived from the dataset itself. However, these existing approaches have their own limitations as follows. Firstly, the scanned objects were supposed to meet various kinds of reflectance models, such as Minnaert model for the reflectance of the surface of the moon [24], Henyey-Greenstein model [25], and Lambertian model [26] in some previous studies. The Lambertian reflectance model describes a perfect diffuse reflector, the cosine law of which is widely used to correct the incident angle effect. However, the perfectly diffuse reflections hardly exist in real applications and the reflectance characteristics of the scanned objects are very complicated in the real world [6,27]. Furthermore, it is not sufficient to correct the incident angle effect for TLS intensity data using the cosθ relation due to the surface roughness or grain size of scanned objects that makes the reflectance characteristics of the scanned objects deviate significantly from the Lambertian model [23]. Secondly, as the system transmission factor includes the brightness reduction [27], automatic gain control (AGC) [28], and receiver optics' defocusing [20,29], the distance effect correction of TLS intensity data is more complex than that of intensity data of airborne laser scanning and does not simply follow the R −2 relation within the entire distance range, especially when the distance is less than 15 m [6]. Thirdly, some homogeneous points need to be selected manually to estimate the surface roughness parameter of scanned object and eliminate the incident angle effect on the intensity data [7]. It needs much more human interactions and consumes more time during the intensity correction. Therefore, a new method is highly desired to improve the correction accuracy of intensity data.
The objective of this research is to develop a new intensity correction method to improve the correction accuracy of TLS intensity data. The contributions of the research include (1) a new intensity correction method is proposed by combining the piecewise fitting and overlap-driven adjustment approaches; and (2) the proposed method is used to retrieve the surface roughness parameter and reflectance of the scanned object, which can be further used for building classification, structural damage detection and road traffic marking identification. The feasibility of the proposed method is validated by experiments and the results show that the reflectance values of all kinds of objects retrieved from the proposed method are very close to those measured by a commercial spectrometer.

Fundamentals of Terrestrial Laser Scanning
The power of light backscattered from the scanned objects in terrestrial laser scanning mainly depends on the system transmission factor of the sensor, the atmospheric effect, distance effect, incident angle effect and the physical characteristic of the scanned objects [8]. As the electromagnetic waves of Remote Sens. 2017, 9, 1090 3 of 16 laser and radar both follow the same principles, the radar range equation can be applied to describe these effects on the received power. The radar range equation can be written as [30,31] where P t and P r are the transmitted power and received power of the sensor, respectively. D r is the receiver aperture diameter of the TLS sensor, R is the distance from the sensor to the scanned object, β t is the transmitter beam width, η sys is the system transmission factor of the sensor, η atm is the atmospheric attenuation factor and σ cross is the backscatter cross section. Assume that the homogeneous target is larger than the laser footprint and its surface is a perfect diffuse reflector that meets the Lambertian model, as shown in Figure 1. The cross section σ cross can be expressed as where θ is the incident angle, which is the angle between the incident light and the surface normal, and ρ λ is the surface reflectance of the scanned object at a certain laser wavelength λ.
where Pt and Pr are the transmitted power and received power of the sensor, respectively. Dr is the receiver aperture diameter of the TLS sensor, R is the distance from the sensor to the scanned object, βt is the transmitter beam width, ηsys is the system transmission factor of the sensor, ηatm is the atmospheric attenuation factor and σcross is the backscatter cross section. Assume that the homogeneous target is larger than the laser footprint and its surface is a perfect diffuse reflector that meets the Lambertian model, as shown in Figure 1. The cross section σcross can be expressed as where θ is the incident angle, which is the angle between the incident light and the surface normal, and ρλ is the surface reflectance of the scanned object at a certain laser wavelength λ.
where the parameters t P, 2 r D and sys η can be considered as three constants, which depend on the sensor only. The atmospheric attenuation factor ηatm can be also considered as a constant during the process of terrestrial laser scanning [32]. Thus, Equation (3) can be further simplified as where C is a constant that depends on the system transmission factor and the atmospheric effects. For a homogeneous object, the received power Pr is proportional to cos θ and R −2 . As the intensity data recorded by the TLS sensor is generally proportional to the received power Pr, the intensity should be also proportional to cosθ and R −2 for the homogeneous object that meets the Lambertian model. However, in many applications, the recorded intensity data is not in a linear relationship with R −2 within the entire distance range, especially at a short distance. Furthermore, it is not sufficient to correct the incident angle effect for TLS intensity data using the cos θ relation because a large surface roughness or grain size will make the scanned object deviate significantly from the Lambertian model. Substituting Equation (2) into Equation (1), the radar range equation can be simplified as where the parameters P t , D 2 r and η sys can be considered as three constants, which depend on the sensor only. The atmospheric attenuation factor η atm can be also considered as a constant during the process of terrestrial laser scanning [32]. Thus, Equation (3) can be further simplified as where C is a constant that depends on the system transmission factor and the atmospheric effects. For a homogeneous object, the received power P r is proportional to cos θ and R −2 . As the intensity data recorded by the TLS sensor is generally proportional to the received power P r , the intensity should be also proportional to cosθ and R −2 for the homogeneous object that meets the Lambertian model. However, in many applications, the recorded intensity data is not in a linear relationship with R −2 within the entire distance range, especially at a short distance. Furthermore, it is not sufficient to correct the incident angle effect for TLS intensity data using the cos θ relation because a large surface roughness or grain size will make the scanned object deviate significantly from the Lambertian model.
Theoretically, the distance and incident angle effects are independent from each other and can be eliminated separately [22,33]. Therefore, the TLS intensity data I can be expressed as where G 1 (R), G 2 (θ) and G 3 (ρ λ ) are the distance effect function, the incident angle effect function and the reflectance effect function, respectively. In this paper, RIEGL VZ-400i (provided by Five-star Electronic Technology Company Limited, Beijing, China), with the wavelength of 1550 nm and range accuracy of 5 mm was used to obtain the 3D point cloud data and the intensity data of the scanned objects [34]. The returned intensity I dB is given in the unit of decibel (dB) for the instrument, which can be expressed as [35] I dB = 10 log(I) = 10 log( P r P th ), where P r is the received power and P th is the detection threshold power of the sensor, which is a constant. Thus, substituting Equation (5) into Equation (6), we obtained where F 1 (R), F 2 (θ) and F 3 (ρ λ ) are decibel forms of G 1 (R), G 2 (θ) and G 3 (ρ λ ), respectively. Therefore, the corrected intensity I c can be described as where I c = F 3 (ρ λ ) depends on the reflectance of the scanned object only. F 1 (R) is determined by the piecewise fitting approach. F 2 (θ) is determined by the Oren-Nayar model using the surface roughness parameter of the scanned object [36]. The Oren-Nayar model is a physical reflectivity model to describe the surface roughness using a single parameter: the standard deviation of the slope angle of facets. Details on solving F 1 (R) and F 2 (θ) are provided below. Some discrete points in a certain distance interval are acquired by the TLS sensor using a reference target. The distance and the incident angle can be calculated using the 3D coordinates of the footprint on the surface of the scanned object, as shown in Figure 2. Based on the geometric principle, the distance R and the incident angle θ can be described as where → n s is the normal vector of the scanned surface, O s is the center of the TLS sensor, and D is a laser footprint on the scanned surface, which can be considered as a point due to its small size for TLS. Thus, the line O s D represents the incident light from the TLS sensor. The normal vector → n s is estimated from the best-fitting plane of a small homogenous area surrounding the point of interest in this study. During the distance experiment, the incident angle θ remains the same to eliminate the incident angle effect on the intensity data, which is simply set as θ = 0 • .

Distance Effect Correction
As detailed system parameters of the TLS sensor are not provided by the manufacturer, it is hard to correct the intensity data directly using the radar range equation. Moreover, the distance effect does not follow the R −2 relation within the entire distance range. Due to the brightness reduction and the receiver optics' defocusing for near distance, the radar range equation cannot be met any more, especially when the distance is less than around 15 m [6]. As the distance effect mainly depends on the instrument, the piecewise fitting approach was applied in this research, and 1 ( ) F R can be expressed as two segments where two segments are split by the separation point (Rsep), n represents the polynomial order of 11 ( ) F R , and 0 1 , , , n a a a  are the coefficients of the polynomial function 11 ( ) F R , respectively. 0 b is the coefficient of 12 ( ) F R , which can be determined by the continuity of 1 ( ) F R . The constraint for the continuity should be set as 11 12 ( The coefficients of 11 ( ) F R are estimated through least squares fitting. The root mean square error (RMSE) of the metric is applied to evaluate the fitting results. 12 ( ) F R is derived from Equation (12).
After these coefficients of the two segments are determined, 1 ( ) F R is substituted into Equation (9) to eliminate the distance effect on the raw intensity data.

Incident Angle Effect Correction
In the past decades, correction of TLS intensity data has been widely studied, where perfect diffuse scattering from the object was commonly assumed and the Lambertian model was applied. However, some scanned objects with rough surfaces deviate significantly from the Lambertian reflectors. Several reflectance models are used to optimize Lambertian model in recent studies, such as Lommel-Seeliger model [6,37], Oren-Nayar model [32,36] and Phong model [38,39]. Lommel-Seeliger model is a scattering model particularly adapted to the dark surfaces [37]. Phong model is an empirical model, which is not sufficiently well adapted to correct intensity data [38]. Oren-Nayar is a physical reflectivity model to describe the surface roughness and more realistic in considering local roughness and micro structure [36]. Compared with other models, Oren-Nayar model is more

Distance Effect Correction
As detailed system parameters of the TLS sensor are not provided by the manufacturer, it is hard to correct the intensity data directly using the radar range equation. Moreover, the distance effect does not follow the R −2 relation within the entire distance range. Due to the brightness reduction and the receiver optics' defocusing for near distance, the radar range equation cannot be met any more, especially when the distance is less than around 15 m [6]. As the distance effect mainly depends on the instrument, the piecewise fitting approach was applied in this research, and F 1 (R) can be expressed as two segments where two segments are split by the separation point (R sep ), n represents the polynomial order of F 11 (R), and a 0 , a 1 , · · · , a n are the coefficients of the polynomial function F 11 (R), respectively. b 0 is the coefficient of F 12 (R), which can be determined by the continuity of F 1 (R). The constraint for the continuity should be set as The coefficients of F 11 (R) are estimated through least squares fitting. The root mean square error (RMSE) of the metric is applied to evaluate the fitting results. F 12 (R) is derived from Equation (12). After these coefficients of the two segments are determined, F 1 (R) is substituted into Equation (9) to eliminate the distance effect on the raw intensity data.

Incident Angle Effect Correction
In the past decades, correction of TLS intensity data has been widely studied, where perfect diffuse scattering from the object was commonly assumed and the Lambertian model was applied. However, some scanned objects with rough surfaces deviate significantly from the Lambertian reflectors. Several reflectance models are used to optimize Lambertian model in recent studies, such as Lommel-Seeliger model [6,37], Oren-Nayar model [32,36] and Phong model [38,39]. Lommel-Seeliger model is a scattering model particularly adapted to the dark surfaces [37]. Phong model is an empirical model, which is not sufficiently well adapted to correct intensity data [38]. Oren-Nayar is a physical reflectivity model to describe the surface roughness and more realistic in considering local roughness and micro structure [36]. Compared with other models, Oren-Nayar model is more consistent with the real laser reflection conditions of natural surfaces and can accurately and quantitatively simulate the luminance of light backscattered from natural surfaces. Therefore, Oren-Nayar model was used to eliminate the incident angle effect in this research. In TLS, the incident and returned lights are coincident and the Oren-Nayar model can be simplified as [32,36] where L in and L out are the transmitted intensity on the scanned surface and returned intensity backscattered from the scanned surface, respectively. The surface roughness parameter ε ∈ (0, π/2) is the standard deviation of the slope angle distribution in radians. According to the Oren-Nayar model, F 2 (θ) can be written as where F 2 (θ) is used as a normalized factor to eliminate the incident angle effect, and the surface roughness parameter ε of F 2 (θ) will be estimated in the next section.

Estimation of Surface Roughness Parameter
The estimation of the surface roughness parameter is important for the intensity correction. As the roughness parameter may differ greatly from object to object and from point to point even for the same kind of object, it should be estimated individually for each point of the object. Therefore, a new method to estimate the surface roughness parameter of the scanned object at each point individually is proposed in this paper, as shown in Figure 3.
Remote Sens. 2017, 9,1090 6 of 16 consistent with the real laser reflection conditions of natural surfaces and can accurately and quantitatively simulate the luminance of light backscattered from natural surfaces. Therefore, Oren-Nayar model was used to eliminate the incident angle effect in this research. In TLS, the incident and returned lights are coincident and the Oren-Nayar model can be simplified as [32,36] cos ( sin tan ) where 2 2 2 2 in L and out L are the transmitted intensity on the scanned surface and returned intensity backscattered from the scanned surface, respectively. The surface roughness parameter ε ∈ (0, π/2) is the standard deviation of the slope angle distribution in radians. According to the Oren-Nayar model, 2 ( ) F θ can be written as where 2 ( ) F θ is used as a normalized factor to eliminate the incident angle effect, and the surface roughness parameter ε of 2 ( ) F θ will be estimated in the next section.

Estimation of Surface Roughness Parameter
The estimation of the surface roughness parameter is important for the intensity correction. As the roughness parameter may differ greatly from object to object and from point to point even for the same kind of object, it should be estimated individually for each point of the object. Therefore, a new method to estimate the surface roughness parameter of the scanned object at each point individually is proposed in this paper, as shown in Figure 3. Firstly, all the distance corrected intensity data sets  . Schematic diagram of the proposed method to estimate the optimal value of the surface roughness parameter using the intensity data from TLS stations p and q in a small homogenous area Ω. R p and R q are the distances from laser points to TLS stations p and q, respectively. θ p and θ q are the incident angles between the surface normal and the incident lights from TLS station p and q, respectively.
Firstly, all the distance corrected intensity data sets q ) (i = 1, 2, . . . , M) were selected from the stations p and q, respectively, where M is the total number of laser points from the two stations in the same small homogenous area Ω. Secondly, to obtain the optimal value of the surface roughness parameter, all the distance corrected intensities from stations p and q are further corrected according to Equation (15) with the surface roughness parameter from 0 • to 90 • in a step interval of 1 • . In the overlapped region, the corrected intensities of the homologous laser points from stations p and q should be the same and be expressed as where q ) are the corrected intensities of the ith laser point from stations p and q in the same small homogenous area Ω, respectively. Owing to the impacts of the background noise and the measurement error of the TLS system, the corrected intensities q ) may slightly differ from each other. The difference can be written as where i and M are the serial and total numbers of the laser points in the small homogenous area Ω. Then, the surface roughness parameter can be obtained by using the objective function O(ε), which is where the objective function O(ε) is the root mean square value of the difference δ(i) (i = 1, 2, . . . , M).
The results of the objective function are calculated using the values of surface roughness parameter from 0 • to 90 • in a step interval of 1 • . Thus, we can get 91 groups of results of the objective function, and the surface roughness parameter is selected as the value that minimizes the result of the objective function.
Once the surface roughness parameter is determined, it is substituted in Equation (9) to correct the intensity data for the entire dataset.

Retrieving Reflectance
After the distance and incident effects are eliminated, the intensity data mainly depends on the surface reflectance of the scanned object. The corrected intensity I c has a logarithmic relationship with the reflectance ρ λ as I c = 10 log(ρ λ ).
Thus, the reflectance ρ λ can be obtained by where the corrected intensity I c can be obtained from Equation (9) and the reflectance ρ λ is expressed as a percentage.
To further validate the proposed method for TLS intensity correction, we used an FieldSpec ® 4 Hi-Res spectrometer (provided by Analytical Spectral Devices (ASD), Inc., Longmont, CO, USA) with the spectral resolution of 10 nm and wavelength accuracy of 0.5 nm [40] to measure the reflectance of the three kinds of objects. The reflectance values of the three kinds of objects measured by the spectrometer at the wavelength of 1550 nm were treated as the reference values. Each reference value was the mean of the reflectance values measured by the spectrometer at 10 equally spaced sample points on a homogeneous area of the corresponding kind of object.

Experiments and Data Acquisition
In this section, two sets of experiments were carried out to verify the proposed method. The first set is on intensity correction to distance effect only using three reference targets with various reflectance values. The second set is on intensity correction to both distance and incident angle effects using multi-station scans.
To obtain the coefficients in Equation (11), three reference targets with a size of 60 cm × 60 cm, Teflon material and reflectance of 15%, 30% and 60% were scanned by the TLS sensor at a series of distances from 5 m to 50 m. The TLS sensor utilized in this research is Riegl VZ-400i. The wavelength of Riegl VZ-400i is 1550 nm, the range accuracy is 5 mm, the angular resolution is better than 0.007 • , and the diameter of the laser beam footprint is 35 mm @100 m [34]. Experimental setup and scene of the distance experiment are shown as Figure 4. To eliminate the incident angle effect on the intensity data, the incident angle was set to 0 • during the scanning process. Each of the three reference targets was firstly scanned at the distance of 5 m, then moved forward at a step interval of 0.3 m until the distance reached 12 m. From the distance of 12 m, the reference targets were moved forward with a step interval of 1.2 m until the distance reached 50 m. For each reference target, intensity data of 50 laser points around the center of the reference target were obtained at each distance, and the means of the intensity data of the laser points were applied to eliminate the distance effect.

Experiments and Data Acquisition
In this section, two sets of experiments were carried out to verify the proposed method. The first set is on intensity correction to distance effect only using three reference targets with various reflectance values. The second set is on intensity correction to both distance and incident angle effects using multi-station scans.
To obtain the coefficients in Equation (11), three reference targets with a size of 60 cm × 60cm, Teflon material and reflectance of 15%, 30% and 60% were scanned by the TLS sensor at a series of distances from 5 m to 50 m. The TLS sensor utilized in this research is Riegl VZ-400i. The wavelength of Riegl VZ-400i is 1550 nm, the range accuracy is 5 mm, the angular resolution is better than 0.007°, and the diameter of the laser beam footprint is 35 mm @100 m [34]. Experimental setup and scene of the distance experiment are shown as Figure 4. To eliminate the incident angle effect on the intensity data, the incident angle was set to 0° during the scanning process. Each of the three reference targets was firstly scanned at the distance of 5 m, then moved forward at a step interval of 0.3 m until the distance reached 12 m. From the distance of 12 m, the reference targets were moved forward with a step interval of 1.2 m until the distance reached 50 m. For each reference target, intensity data of 50 laser points around the center of the reference target were obtained at each distance, and the means of the intensity data of the laser points were applied to eliminate the distance effect. In order to estimate the surface roughness parameter and validate the feasibility of the proposed method, the second set of experiments was conducted to obtain multi-station scans. Before intensity correction, the multi-station scans were preprocessed and registered using the RiSCAN PRO software (v2.42, RIEGL Laser Measurement Systems GmbH, Vienna, Austria). In the registration process, some corresponding points from multi-station scans were selected manually to make coarse registration. Then the fine registration was completed using the Iterative Closest Point (ICP) algorithm [41]. Figure  5 is the point cloud of the multi-station scans after registration, which is colored with raw intensity data and displayed by the CloudCompare software (v2.9, D. Girardeau-Montaut, EDF Research and Development, Telecom Paris Tech, Paris, France) [42]. It consists of three kinds of objects including paving brick, concrete road and marking line, which are numbered as 1, 2 and 3, respectively. In order to estimate the surface roughness parameter and validate the feasibility of the proposed method, the second set of experiments was conducted to obtain multi-station scans. Before intensity correction, the multi-station scans were preprocessed and registered using the RiSCAN PRO software (v2.42, RIEGL Laser Measurement Systems GmbH, Vienna, Austria). In the registration process, some corresponding points from multi-station scans were selected manually to make coarse registration. Then the fine registration was completed using the Iterative Closest Point (ICP) algorithm [41]. Figure 5 is the point cloud of the multi-station scans after registration, which is colored with raw intensity data and displayed by the CloudCompare software (v2.9, D. Girardeau-Montaut, EDF Research and Development, Telecom Paris Tech, Paris, France) [42]. It consists of three kinds of objects including paving brick, concrete road and marking line, which are numbered as 1, 2 and 3, respectively.

Results of Distance Effect Correction
To implement the intensity correction, the distance effect should be first corrected. The variations of raw intensity data with distance for the three reference targets are shown in Figure 6. We can see that the raw intensity data increases firstly with increase of the distance up to around 12

Results of Distance Effect Correction
To implement the intensity correction, the distance effect should be first corrected. The variations of raw intensity data with distance for the three reference targets are shown in Figure 6. We can see that the raw intensity data increases firstly with increase of the distance up to around 12 m for all the three reference targets. Then, the intensity decreases with increase of the distance until 50 m. When the distance is shorter than around 20 m, the system transmission factor, such as the brightness reducer or the receiver optics' defocusing effect, has a significant effect on the raw intensity data, and it makes that the intensity data dissatisfy the radar range equation (refer to Equation (4)). When the distance is longer than around 20 m, the distance effect is the main factor that influences the raw intensity data and other effects can be ignored or remain unchanged with distance. Thus, we can see that the intensity data of each of the three reference targets is proportional to 10 × log(b 0 /R 2 ) when the distance is longer than 20 m, which meets the radar range equation. Therefore, the separation point is selected as R sep = 20 m in this research.

Results of Distance Effect Correction
To implement the intensity correction, the distance effect should be first corrected. The variations of raw intensity data with distance for the three reference targets are shown in Figure 6. We can see that the raw intensity data increases firstly with increase of the distance up to around 12 m for all the three reference targets. Then, the intensity decreases with increase of the distance until 50 m. When the distance is shorter than around 20 m, the system transmission factor, such as the brightness reducer or the receiver optics' defocusing effect, has a significant effect on the raw intensity data, and it makes that the intensity data dissatisfy the radar range equation (refer to Equation (4)). When the distance is longer than around 20 m, the distance effect is the main factor that influences the raw intensity data and other effects can be ignored or remain unchanged with distance. Thus, we can see that the intensity data of each of the three reference targets is proportional to 10 × log(b0/R 2 ) when the distance is longer than 20 m, which meets the radar range equation. Therefore, the separation point is selected as Rsep = 20 m in this research. Figure 6. Variation of raw intensity data with distance for all of the three reference targets with the reflectance of 15%, 30% and 60% at the incident angle of 0°.
The intensity data within a short distance range, less than 20 m, were fit by the nonlinear least square estimation, and the fitting results relative to the polynomial order (n) of the function are shown in Figure 7. The intensity data within a short distance range, less than 20 m, were fit by the nonlinear least square estimation, and the fitting results relative to the polynomial order (n) of the function are shown in Figure 7.
In the polynomial fitting process, the RMS values of fitting errors and coefficient of determination (R-square) are used to evaluate the influence of polynomial order (n) on the fitting results in this paper. From Figure 7, we know that when the value of polynomial order is larger than 3, the R-square and RMS values of fitting errors remain almost unchanged. Moreover, the R-square is 0.962 and very close to 1, and the RMS value is only 0.261 dB at the polynomial order of 3. Therefore, the optimal order of the polynomial function F 11 (R) was selected as n = 3. According to the continuity of the function, the coefficient b 0 of the function F 12 (R) can be calculated as b 0 = 3.218 × 10 5 using Equations (11) and (12). After the fitting process, the two segments of F 1 (R) were obtained and expressed as F 1 (R) = F 11 (R) = 1.623 × 10 −3 R 3 − 9.287 × 10 −2 R 2 + 1.367 R + 25.88, R < 20 F 12 (R) = 10 log(3.218 × 10 5 /R 2 ), R ≥ 20 .
After calculating for all the intensity data, the RMS value of fitting errors of the polynomial model F 1 (R) is 0.296 dB, and the R-square is 0.984, which is very close to 1. Therefore, the polynomial model has a good fitting result. In the polynomial fitting process, the RMS values of fitting errors and coefficient of determination (R-square) are used to evaluate the influence of polynomial order (n) on the fitting results in this paper. From Figure 7, we know that when the value of polynomial order is larger than 3, the R-square and RMS values of fitting errors remain almost unchanged. Moreover, the R-square is 0.962 and very close to 1, and the RMS value is only 0.261 dB at the polynomial order of 3. Therefore, the optimal order of the polynomial function 11 ( ) F R was selected as n = 3. According to the continuity of the function, the coefficient b0 of the function 12 ( ) F R can be calculated as b0 = 3.218 × 10 5 using Equations (11) and (12). After the fitting process, the two segments of 1 ( ) F R were obtained and expressed as After calculating for all the intensity data, the RMS value of fitting errors of the polynomial model 1 ( ) F R is 0.296 dB, and the R-square is 0.984, which is very close to 1. Therefore, the polynomial model has a good fitting result.

Results of Incident Angle Effect Correction
Assume that Lin = 1, the simulated results of the Oren-Nayar model for different surface roughness parameters are shown in Figure 8. For flat surfaces without facet variation, i.e., ε = 0°, the Oren-Nayar model is equal to the Lambertian model, shown in the blue curve of Figure 8. From the results of Figure 8a, we know that the normalized intensity data decreases with the increase of the incident angle, and it also decreases with the increase of ε at the same time for incident angle less than about 50°. However, for incident angle larger than about 50°, it increases with the increase of ε at the same incident angle. Furthermore, the absolute deviations of  Figure   8b. It shows that the deviation is as large as 0.4 for the surface roughness parameter ε of 60°, even for the small surface roughness parameter ε of 15°, and the deviation also reaches up to 0.2. Thus, the surface roughness parameter has a large effect on the intensity data. Moreover, the rougher the surface of the scanned object, the more significant the deviation of the Oren-Nayar model from the Lambertian model. Only if the incident angle is around in the interval [35°, 50°], the deviation is small, less than 0.05, as shown in the red box of Figure 8b.
As mentioned, we know that the surface roughness parameter is crucial for intensity correction, and different surface roughness parameter will generate different correction results for the same object. When the Lambertian model is used to correct the intensity data, the surface roughness

Results of Incident Angle Effect Correction
Assume that L in = 1, the simulated results of the Oren-Nayar model for different surface roughness parameters are shown in Figure 8. For flat surfaces without facet variation, i.e., ε = 0 • , the Oren-Nayar model is equal to the Lambertian model, shown in the blue curve of Figure 8. From the results of Figure 8a, we know that the normalized intensity data decreases with the increase of the incident angle, and it also decreases with the increase of ε at the same time for incident angle less than about 50 • . However, for incident angle larger than about 50 • , it increases with the increase of ε at the same incident angle. Furthermore, the absolute deviations of L ε=15 • out , L ε=30 • out , L ε=45 • out , L ε=60 • out from L ε=0 • out decrease firstly and then increase with the increase of the incident angle, as shown in Figure 8b. It shows that the deviation is as large as 0.4 for the surface roughness parameter ε of 60 • , even for the small surface roughness parameter ε of 15 • , and the deviation also reaches up to 0.2. Thus, the surface roughness parameter has a large effect on the intensity data. Moreover, the rougher the surface of the scanned object, the more significant the deviation of the Oren-Nayar model from the Lambertian model. Only if the incident angle is around in the interval [35 • , 50 • ], the deviation is small, less than 0.05, as shown in the red box of Figure 8b.
As mentioned, we know that the surface roughness parameter is crucial for intensity correction, and different surface roughness parameter will generate different correction results for the same object. When the Lambertian model is used to correct the intensity data, the surface roughness parameter cannot be involved, resulting in correction results deviating greatly from the real intensity data that depend on the physical characteristic of the scanned objects. Only when the surface roughness parameter is considered can the corrected intensity data better approach the real intensity data. Therefore, the surface roughness parameter should be estimated for each kind of object before being used for intensity correction.
To eliminate the incident angle effect on the TLS intensity data, the surface roughness parameter of the three kinds of objects should be determined according to Equations (17) and (18). The results are shown in Table 1. It can be seen that the three kinds of objects have different surface roughness parameter, the means of which are 20.6 • , 17.9 • and 20.8 • , respectively. The mean and standard deviation (Std.) of surface roughness parameter of the paving brick, concrete road and marking line are large. The divergence of the surface roughness parameter of each of the three kinds of objects indicates that the surface roughness parameter differs from position to position even for the same kind of object. Therefore, the surface roughness parameter of the same kind of object should be estimated individually from point to point.
Remote Sens. 2017, 9, 1090 11 of 16 parameter cannot be involved, resulting in correction results deviating greatly from the real intensity data that depend on the physical characteristic of the scanned objects. Only when the surface roughness parameter is considered can the corrected intensity data better approach the real intensity data. Therefore, the surface roughness parameter should be estimated for each kind of object before being used for intensity correction.
To eliminate the incident angle effect on the TLS intensity data, the surface roughness parameter of the three kinds of objects should be determined according to Equations (17) and (18). The results are shown in Table 1. It can be seen that the three kinds of objects have different surface roughness parameter, the means of which are 20.6°, 17.9° and 20.8°, respectively. The mean and standard deviation (Std.) of surface roughness parameter of the paving brick, concrete road and marking line are large. The divergence of the surface roughness parameter of each of the three kinds of objects indicates that the surface roughness parameter differs from position to position even for the same kind of object. Therefore, the surface roughness parameter of the same kind of object should be estimated individually from point to point. In order to eliminate the distance and incident angle effects on intensity data, three methods were used to correct the raw intensity data in this research, including the proposed method, the Tan method [7] and the Lambertian model. In the Tan method, the distance effect is eliminated by reference targets not using the distance effect function and the surface roughness parameter is estimated by using the average surface roughness parameter of 20 manually selected points for each kind of object. However, in the proposed method, the local surface roughness is estimated from point to point by using the intensity data of corresponding points in the overlapped regions of two stations, rather than being estimated by averaging the surface roughness at 20 manually selected points for each kind of object. The distributions of raw and corrected intensity data for the three kinds of objects from TLS stations 1 and 2 are shown in Figure 9.
Firstly, it can be seen that, as the incident angle and distance in station 1 are different from those in station 2, the raw intensity data of all the three kinds of objects in station 1 are significantly different from those in station 2. The differences of the distribution of raw intensity between stations 1 and 2 are as large as 7 dB for the three kinds of objects, and the distribution widths are from 4 dB to 7 dB for the three kinds of objects in stations 1 and 2. However, the distribution of corrected intensity   In order to eliminate the distance and incident angle effects on intensity data, three methods were used to correct the raw intensity data in this research, including the proposed method, the Tan method [7] and the Lambertian model. In the Tan method, the distance effect is eliminated by reference targets not using the distance effect function and the surface roughness parameter is estimated by using the average surface roughness parameter of 20 manually selected points for each kind of object. However, in the proposed method, the local surface roughness is estimated from point to point by using the intensity data of corresponding points in the overlapped regions of two stations, rather than being estimated by averaging the surface roughness at 20 manually selected points for each kind of object. The distributions of raw and corrected intensity data for the three kinds of objects from TLS stations 1 and 2 are shown in Figure 9.
Firstly, it can be seen that, as the incident angle and distance in station 1 are different from those in station 2, the raw intensity data of all the three kinds of objects in station 1 are significantly different from those in station 2. The differences of the distribution of raw intensity between stations 1 and 2 are as large as 7 dB for the three kinds of objects, and the distribution widths are from 4 dB to 7 dB for the three kinds of objects in stations 1 and 2. However, the distribution of corrected intensity between stations 1 and 2 are almost the same, and the distribution widths are less than 3 dB for the three kinds of objects. Therefore, the raw intensity data should be corrected before used.
Secondly, after being corrected by the Lambertian model, the distributions of corrected intensity data of all the three kinds of objects from station 1 still differ greatly from those from station 2. Especially for the marking line, the difference of distributions of the corrected intensity data between stations 1 and 2 is as large as 6 dB. After being corrected by the Tan method, the distributions of intensity data are better than those corrected by the Lambertian model. However, the difference still exists between the distributions of intensity data from stations 1 and 2.
Compared with the distributions of intensity data corrected by the other two methods, the distributions of intensity data corrected by the proposed method from station 1 are almost the same as those from station 2 for all the three kinds of objects. Therefore, the correction accuracy using the proposed method is better than those using the other two methods.
Especially for the marking line, the difference of distributions of the corrected intensity data between stations 1 and 2 is as large as 6 dB. After being corrected by the Tan method, the distributions of intensity data are better than those corrected by the Lambertian model. However, the difference still exists between the distributions of intensity data from stations 1 and 2.
Compared with the distributions of intensity data corrected by the other two methods, the distributions of intensity data corrected by the proposed method from station 1 are almost the same as those from station 2 for all the three kinds of objects. Therefore, the correction accuracy using the proposed method is better than those using the other two methods. Figure 9. Distributions of raw and corrected intensity data for the three kinds of objects from TLS stations 1 and 2. (a,c,e) are the distributions of raw intensity data for paving brick, concrete road and marking line, respectively; (b,d,f) are the distributions of intensity data corrected by the three methods for paving brick, concrete road and marking line, respectively. (For the box chart, the three long horizontal lines represent the lower quartile, median and upper quartile from down to up, respectively. The small circle represents the mean value. The two short horizontal lines are the lower whisker and upper whisker of the box, respectively. The two signs 'x' at both ends represent the 1% and 99% range of the box [43]).

Reflectance Retrieved from Corrected Intensity Data
The means and standard deviations of the reflectance values for the three kinds of objects obtained by the spectrometer, the proposed method, the Tan method and the Lambertian model are shown in Table 2. From the results, it can be seen that the deviations of the mean reflectance values obtained by Lambertian model from those obtained by the spectrometer are larger than 37%. The deviations of the mean reflectance values obtained by the Tan method from those obtained by the spectrometer are 6.2%, 6.7% and −6.0% for the three kinds of objects, respectively. The deviations of the mean reflectance values obtained by the proposed method from those obtained by the spectrometer are 1.0%, 2.5% and −0.8% for the three kinds of objects, respectively, and the standard deviations of the reflectance values obtained by the proposed method are close to those obtained by the spectrometer. The mean reflectance values retrieved from the proposed method for all three kinds of objects are closer to the reference values than those retrieved from the Tan method and Lambertian model. Overall, we can say that the proposed method is better than the Tan method and provides higher accuracy reflectance of the scanned object.

Discussion
In this study, the intensity correction combining the piecewise fitting and overlap-driven adjustment approaches has been carried out for TLS. Results of the distance experiment performed with three different reflectance targets demonstrated that the intensity values recorded by the used TLS sensor cannot follow the 1/R 2 relation at near distance, which is consistent with some previous studies [6,20]. However, at far distance (about 20 m) the intensity values are proportional to the range (1/R 2 ), which overcomes the technical issue of defocusing on TLS instruments. Thus, the distance effect function is established using the piecewise fitting approach. The incident angle effect function is established using the Oren-Nayar model that employs the surface roughness parameter of the scanned object. The distance effect and incidence angle effect on intensity data are corrected independently as previous studies suggested [7,33].
Results of incident angle effect correction show that as the surface roughness parameter of the three kinds of objects are large, their surface characteristics do not meet the Lambertian model, and intensity correction using the Lambertian model generates large deviations. Therefore, the Lambertian model is not suitable for the intensity correction of the scanned object with a large surface roughness parameter, and not as good as the Tan method. For the marking line, the surface roughness parameter and the difference of incident angles between stations 1 and 2 are both large so that the intensity correction using the Lambertian model generates large deviation. If the surface of scanned object is homogenous and has no inhomogeneous points, such as crack points, contaminants and pits, the Tan method is suitable and can eliminate the random measurement error of the TLS sensor through the 20 manually selected points for each kind of object. However, the scanned object often has structural damage or contaminants on its surface in the real condition. When the manually selected points include inhomogeneous points, the estimation of surface roughness parameter will generate a large deviation that can further affect the following intensity correction. As the local surface roughness parameter of the scanned object is estimated using the intensity data of a small area around the corresponding point in the overlapped region of the multi-station scans in the proposed method, the estimation of the local surface roughness parameter at each point is not affected by the inhomogeneous points at other positions of the scanned object like the Tan method, and hence more accurate. Therefore, the proposed method can improve the correction accuracy of intensity data and obtain more accurate reflectance for the three kinds of objects than that derived from the Tan method.
The advantage of the proposed method is that the local surface roughness parameter is considered in TLS intensity correction and estimated by using the intensity data of the corresponding points in the overlapped regions of multi-station scans, where there is no need to manually select the homogenous points as the Tan method described therein. Moreover, because the points selected manually in the Tan method can include inhomogeneous points, the estimation of surface roughness parameter will generate a larger deviation that can further affect the following intensity correction. Thus, the proposed method for TLS intensity correction is better than the Tan method. However, the proposed method for the intensity correction seems to have some limitations. First, the proposed method cannot be used to well correct the intensity data for the objects with glossy surfaces, especially when the incident angle is 0 • . Because the Oren-Nayar is a physical model, which only considers the diffuse reflection, without the specular reflection. Although mostly natural objects meet the Oren-Nayar model, the objects with smooth or glossy surfaces, such as the glasses and glazed tile, cannot meet the Oren-Nayar model. Second, the proposed method is not suitable for a long-range TLS. As the size of laser footprint is large for the long-range TLS, the small homogenous area between two stations hardly exists and even the laser footprint is larger than the scanned object.

Conclusions
In this study, a new method was proposed to correct the intensity data for TLS by combining the piecewise fitting and overlap-driven adjustment approaches. In the method, the local surface roughness parameter is estimated from point to point by using the intensity data of corresponding points in the overlapped regions of multi-station scans, rather than being estimated by averaging the surface roughness parameter at 20 manually selected points for each kind of object in the Tan method. Our method estimating the local surface roughness parameter is not influenced by the inhomogeneous points, whereas the Tan method is affected by the inhomogeneous points due to selecting the points at different positions on the scanned object. Thus, the estimation using the proposed method is more accurate than that using the Tan method.
The experimental results demonstrate that the distribution of corrected intensity data between stations 1 and 2 are almost the same, and the distribution widths of corrected intensity data are reduced to 3 dB compared with the distribution width of 7 dB for raw intensity data. The reflectance values derived from the proposed method are much closer to the reference values obtained by the spectrometer than those derived from the Tan method, and the reflectance deviations are all less than 3% for the three kinds of objects. It proves that the proposed method can improve the correction accuracy of intensity data, and hence be better applied for building classification, structural damage detection and road traffic marking identification.