CORRECTION OF OVERLAPPING MULTISPECTRAL LIDAR INTENSITY DATA : POLYNOMIAL APPROXIMATION OF RANGE AND ANGLE EFFECTS

Recent development of radiometric calibration, correction and normalization approaches have facilitated the use of monochromatic LiDAR intensity and waveform data for land surface analysis and classification. Despite the recent successful attempts, the majority of existing approaches are mainly tailor made for monochromatic LiDAR toward specific land surface scenario. In view of the latest development of multispectral LiDAR sensor, such as the Optech Titan manufactured by Teledyne Optech, a more generic approach should be developed so that the radiometric correction model is able to handle and compensate the laser energy loss with respect to different wavelengths. In this study, we propose a semi-physical approach that aims to utilize high order polynomial functions to model the distortion effects due to the range and the angle. To estimate the parameters of the respect polynomial functions for the range and angle, our approach first locates a pair of closest points within the overlapping LiDAR data strips and subsequently uses a non-linear least squares adjustment to retrieve the polynomial parameters based on the Levenberg-Marquardt algorithm. The approach was tested on a multispectral airborne LiDAR dataset collected by the Optech Titan for the Petawawa Research Forest located in Ontario, Canada. The experimental results demonstrated that the coefficient of variation of the intensity of channel 1 (1550 nm), channel 2 (1064 nm) and channel 3 (532 nm) were reduced by 0.1% to 39%, 10% to 45% and 12% to 54%, respectively. The striping noises, no matter found within single strip and overlapping strips, were significantly reduced after implementing the proposed radiometric correction.


INTRODUCTION
Remote sensing community has evidenced a boom of sensor development, particularly the Light Detection And Ranging (Li-DAR) sensor, for Earth observation and planetary exploration.Among which the ICESat (Schutz et al., 2005), Mars Orbiter Laser Altimeter (Smith et al., 2001), and Lunar Orbiter Laser Altimeter (Smith et al., 2010) were launched for the exploration of the Earth, Mars and Lunar, respectively, during the last decade.Due to the operation of monochromatic laser, the LiDAR data being collected are mainly used for topographic mapping or geometric measurement, instead of studying the Earth and planetary surface properties.In view of the forthcoming space borne LiDAR mission, it is foreseeable that future space mission will adopt the use of multi-wavelength lasers in order to study the Earth and planetary surface properties, through analyzing the backscattered laser signal strength.
Teledyne Optech has officially manufactured the world's first commercial multispectral airborne LiDAR sensor, Optech Titan, where the sensor is capable of collecting laser point clouds at three different wavelengths (532 nm, 1064 nm and 1550 nm) simultaneously.Such a break-through in LiDAR technology compensates the drawback of high resolution satellite remote sensing imagery in various ways and opens many doors in topographic and environmental applications.Recent attempts can be found in different applications using multispectral LiDAR data, including land cover classification (Wang et al., 2014;Matikainen et al., 2017;Morsy et al., 2017), tree species classification (Yu et al., 2017), road extraction (Karila et al., 2017), water mapping (Morsy et al., 2016).One can foresee the multispectral LiDAR sensors will be deployed in different platforms in near future.
With various radiometric pre-processing techniques being proposed and developed (Kashani et al., 2015), LiDAR data users reap the benefits of using the monochromatic laser intensity data for surface classification and object recognition (Yan et al., 2015).Most of the existing radiometric calibration and correction models are built upon the radar (range) equation that considers the system-and environment induced distortions, including the range (Kaasalainen et al., 2011), incidence angle (You et al., 2017), atmospheric attenuation (Errington and Daku, 2017) and other system parameters (Vain et al., 2010).To minimize the discrepancy in overlapping strips, certain radiometric normalization approaches have been proposed (Yan andShaker, 2014, 2016) which have proven to reduce the striping noise found within the overlapping region.Indeed, all these radiometric pre-processing techniques not only improve the intensity data quality, but also maximize the benefits of using the intensity data in different application domains.
Despite the above-mentioned attempts, existing radiometric calibration and correction approaches mainly consider the the square of range and the inverse cosine correction of angle.Such a setting can only handle simple land cover scenario and stable system settings for monochromatic LiDAR system.Recent researches further explore the use of different bidirectional reflectance distribution functions (BDRF), such as Phong model (Jutzi and Gross, 2010;Ding et al., 2013), Oren-Nayar reflectance model (Carrea et al., 2016;Tan et al., 2016) .In view of the latest development of multispectral LiDAR sensor, a more generic approach should be developed so that the radiometric correction model is capable of handling and compensating the laser energy loss with respect to different wavelengths and land cover scenarios.This thus inspires us to develop a semi-physical model that built upon the use of radar (range) equation, where the range and angle factors are represented by a respective higher order polynomial function.

PROPOSED RADIOMETRIC CORRECTION
The majority of existing radiometric correction model for LiDAR intensity data is built upon the radar (range) equation (Jelalian, 1992) that describes the relationship of the backscattered laser energy with respect to various system and environmental parameters as shown in the following equation: where Pr refers to the received laser energy, Pt refers to the transmitted laser energy, Dr is the diameter aperture, R is the range, ηsys and ηatm are the system and atmospheric attenuation, respectively, and finally σ describes the characteristics of the laser pulse being backscattered from the ground objects: where ρs refers to the spectral reflectance of the surface and θ refers to the incidence angle.Since most of the existing studies assume Lambertian reflectance occurred in the ground surface; therefore cosine of incidence angle is being used to describe the diffuse reflectance.With several parameters assumed as constant and unchanged during the flight survey, the radiometric correction model is commonly simplified as follows: where Ic refers to the corrected intensity data and I refers to the original intensity data.Although the above Eq. 3 is commonly used for radiometric calibration and radiometric correction by normalizing the data with either 1/R 2 minimum or 1/R 2 median , the use of cosine correction for incidence angle seems to be inefficient for laser pulses collected by different wavelengths and different land cover scenarios.Our previous study also pointed out the use of cosine incidence angle may induce overcorrection in certain scenarios (Yan and Shaker, 2014).Therefore, we propose to use polynomial approximation to model the laser attenuation effects due to the range and angle loss.In this case, Eq. 3 becomes To estimate the polynomial coefficients ai and bi that are used to correct the effects of range and incidence angle, we reply on the use of overlapping data strips to estimate these coefficients.Similar to our previous work (Yan and Shaker, 2016), assume there exist a set of multispectral LiDAR data having two overlapping strips, i.e., X (A j ) and X (B j ) , where j refers to the data channel, the process begins by looking for a set of closest points between the two overlapping data strips in each of the data channels within a specific threshold.Once a set of data point is being paired up from the two data strips, i.e., x (A j ) and x (B j ) , we assume their corrected intensity (spectral reflectance) would be identical, i.e.
. Therefore, Eq. 4 can be expressed as: Since the above Eq.5 is a non-linear equation, non-linear least squares adjustment should be conducted in order to estimate the parameter of ai and bi with respect to the original intensity data, range and incidence angle provided by the matched data points from the two strips, i.e.I (A j ) and I (B j ) , R (A j ) and R (B j ) , and θ (A j ) and θ (B j ) , respectively.To yield a robust estimation of the non-linear least squares adjustment, the Levenberg-Marquardt algorithm is implemented to estimate the polynomial coefficients.
Once the coefficients are solved, they are applied back to the Eq. 4 in order to perform the radiometric correction for the corresponding channel j of the two LiDAR data strips, and the abovementioned process is implemented for all the available channels independently.

Multispectral LiDAR Dataset
The multispectral LiDAR datasets were collected for the Laurentian Hills (approximately 350 m above mean sea level), situated at the northwest of the town of Petawawa (45

Data Processing
The LiDAR data files were provided in both las format and ASCII format, where the following data fields were extracted for computation: x-coordinate, y-coordinate, z-coordinate, intensity, range, scan angle, return number, total number of returns, and scan direction flag.Prior to the radiometric correction, we observed a mild intensity striping noise occurred in each of the individual LiDAR data strip in channel 1 and channel 2. The reason is due to the effect of "intensity banding", where a portion of the laser beam may wander off the detector edge in a specific scan direction.Therefore, a pre-processing was first conducted to split each of the individual data strip into positive scan direction and negative scan direction based on the "scan direction flag" in the las file.After that, we implemented our previously proposed radiometric normalization technique (Yan and Shaker, 2016) to normalize the dataset of positive scan direction with respect to the dataset of negative scan direction based on a second order polynomial model.After the normalization, the LiDAR dataset with both positive and negative scans were combined and the observed "intensity banding" effects were significantly reduced.Subsequently, the proposed radiometric correction was implemented.
As mentioned in section 2, the process began by pairing up the closest points within the overlapping region of the three strips.Therefore, kd-tree was built on the LiDAR data strips, regardless of the data channels, in order to speed up the searching of closest points in each pair of the overlapping LiDAR data strips.Half of the meaning point spacing was used as a threshold to pair up the closest points within the overlapping region.The range and angle polynomial models, as noted in Eq. 4, were set to be in the third order based on the experimental trials.To initialize the non-linear least squares solution, only a2 and b1 were set as 1 and the rest of the parameters were assigned as zero.Such an initialization settings deemed to be similar as the traditional radiometric correction being commonly used (Yan et al., 2012).Since the study area covers mainly the forest canopies, the use of incidence angle for radiometric correction may induce certain overcorrection effects as reported in our previous study (Yan and Shaker, 2014).Therefore, we only used the scan angle instead of incidence angle for the radiometric correction.

Experimental Evaluation
Similar to our previous work (Yan andShaker, 2014, 2016), several land cover samples were collected in order to compute the mean and standard deviation of the intensity data before and after radiometric correction.Four types of land cover samples, including bare soil, tree canopies, grass cover and road, were selected within the overlapping LiDAR data strips as well as the individual strips in order to assess the coefficient of variation (cv), which can be used to assess the homogeneity of the selected land cover feature (ω).
Such an evaluation can help to quantitatively assess the level of noise within the intensity data for all the three data channels before and after radiometric correction.If the value of cv is reduced after radiometric correction, one can conclude that the radiometric correction process can help to reduce the intensity noise.

RESULTS AND ANALYSIS
Fig. 2 shows the cv of the four land cover samples computed for channel 1, channel 2 and channel 3 before and after radiometric correction.One can easily observe a reduction of cv in the majority of the land cover samples after performing the radiometric normalization for "intensity banding" and the proposed radiometric correction.In channel 1 (Fig. 2(a)), the road data samples recorded a reduction of cv by 39% (from 0.407 to 0.369 after radiometric correction).The soil and grass samples had a cv 0.198 and 0.407 before radiometric correction, and their corresponding cv recorded a reduction by 15% and 9% after implementing the proposed correction.Among all the data samples, tree canopies had the least reduction of cv (only 0.1%).Figure 2. Coefficient of variation (cv) of four selected land cover samples before and after radiometric correction.
Among the three LiDAR data channels, channel 2 had the serious striping noises appeared in each of the single strip as well as the overlapping strips.Therefore, significant improvement of intensity discrepancy was found after running the normalization (issue regarding the "intensity banding") and radiometric correction.Unlike channel 1, the tree samples recorded a mild reduction The proposed radiometric correction performed least well among the three data channels.One of the reasons can be explained by the relative low reflectance of land covers in the green wavelength (532 nm), and thus the striping noise may not appear obvious.
In addition, the LiDAR sensor of channel 3 did not suffer from the problem of "intensity banding".Both soil and grass samples had a similar cv value at around 0.196.After radiometric correction, both cv values decreased to 0.172 (soil) and 0.150 (grass) by 12% and 23%, respectively.The road samples performed the best among the four samples with a reduction of cv by 54%.However, the absolute value of cv in channel 3 was comparatively low than the other two channels.Surprisingly, the tree samples recorded an increase of cv from 0.2 to 0.212 after radiometric correction.Despite these, the impact of radiometric correction on channel 3 was not significant at all.As shown in both Figs. 3(d) and 3(i), both multispectral LiDAR intensity images suffered from systematic line striping noise that are obvious on the ground level (both bare soil, road and grass cover).Such striping noises were not found visually obvious on the tree canopies, but they were still noticeable.After running the proposed data correction scheme, the line striping noises were significantly reduced.Among the three LiDAR channels, channel 2 had the highest level of striping noise (see Figs. 3(b) and 3(j)) in the original intensity data.Therefore, the the corrected intensity of channel 2 had an obvious reduction of noise (see Figs. 3(f) and 3(n)).Channels 1 and 3 did not suffer from the high level of striping noise as channel 2, and thus the visual impact of radiometric correction was not as noticeable as that of channel 2. Regardless of the LiDAR data channels, there still existed a low level of intensity noise due to the the occurrence of shadows and occlusions on tree crowns.In addition, significant energy loss was found along the edge of the scan, and this thus caused a high intensity discrepancy comparing to those at the nadir.

CONCLUSIONS
Recent development of multispectral airborne LiDAR has deemed to be a break-through in the remote sensing community, and such a cutting edge technology can improve the analytical capability and accuracy in different remote sensing applications.
In order to maximize the benefits of using the multispectral Li-DAR intensity data for thematic analysis, radiometric correction should be first conducted to remove or reduce the system-and environmental-induced distortions.In spite of the existing radiometric calibration and correction techniques for monochromatic LiDAR, there is a need for further refinement of the radiometric correction model with respect to the various wavelengths.We thus propose a semi-physical correction model that is built upon the radar (range) equation, where the effects of range and angle are modelled by a respective polynomial equation.Such a setting provides a flexible and approximate description of different types of reflection through modelling the intensity discrepancy within the overlapping data strips.Our experimental work demonstrated an improvement of radiometric quality after implementing the proposed radiometric correction model by up to 39% to 54% in the respective three LiDAR data channels collected by Optech Titan for the Petawawa Research Forest, Ontario, Canada.The line striping noises were obviously reduced after implementing the proposed correction model.Further effort should be placed toward incorporating such polynomial based approximation of range and angle with the considerations of both specular and diffuse components, and fusing with different BRDFs.

Figure 1 .
Figure 1.Study area and multispectral LiDAR data strips.
Figure 3. Multispectral LiDAR intensity image before and after radiometric correction.

Fig. 3
Fig.3shows the multispectral LiDAR intensity image combined from the three channels with band combination: red band = chan-

Table 1 .
• 53' N, 77 • 16' W), Ontario, Canada (see Fig.1).A subset of the data was clipped, in which it covers south of the Petawawa Research Forest.The multispectral LiDAR dataset contains three overlapping strips with more than 50% overlap, in which they have an approximate extent of 1.5 km (across track) by 21.5 km (along track).The flight mission was accomplished on July 20 th , 2016.The Optech Titan sensor was operated with scan angle ±24 • , pulse repetition frequency 300 kHz per channel and flying attitude 1,100 m above sea level.With these settings, the mean point density yields approximately 4.36 points/m 2 , 4.66 points/m 2 , and 1.61 points/m 2 for the channel 1, channel 2 and channel 3, respectively.Table1summarizes the LiDAR system settings and data specification.LiDAR system settings and data specification.