EVALUATING THE POTENTIAL OF MULTISPECTRAL AIRBORNE LIDAR FOR TOPOGRAPHIC MAPPING AND LAND COVER CLASSIFICATION

Recently multispectral LiDAR became a promising research field for enhanced LiDAR classification workflows and e.g. the assessment of vegetation health. Current analyses on multispectral LiDAR are mainly based on experimental setups, which are often limited transferable to operational tasks. In late 2014 Optech Inc. announced the first commercially available multispectral LiDAR system for airborne topographic mapping. The combined system makes synchronic multispectral LiDAR measurements possible, solving time shift problems of experimental acquisitions. This paper presents an explorative analysis of the first airborne collected data with focus on class specific spectral signatures. Spectral patterns are used for a classification approach, which is evaluated in comparison to a manual reference classification. Typical spectral patterns comparable to optical imagery could be observed for homogeneous and planar surfaces. For rough and volumetric objects such as trees, the spectral signature becomes biased by signal modification due to multi return effects. However, we show that this first flight data set is suitable for conventional geometrical classification and mapping procedures. Additional classes such as sealed and unsealed ground can be separated with high classification accuracies. For vegetation classification the distinction of species and health classes is possible.


INTRODUCTION
Commercial airborne laser scanning (ALS) systems are commonly monochromatic systems recording either discrete echoes or the full waveform (FWF) of the reflected laser beam.Accordingly to multispectral and hyperspectral remote sensing using passive optical sensors there have already been several studies conducted making use of multiple wavelengths for topographic LiDAR mapping..Recent studies present experimental multispectral or hyperspectral LiDAR systems, which are used for specific tasks in laboratories or as terrestrial laser scanning (TLS) platforms.
Many studies use systems for vegetation analysis.Danson et al. (2014) present the capabilities of the Salford Advanced Laser Canopy Analyser (SALCA), which is a hemispherical dualwavelength full-waveform TLS.The scanner operates in nearand middle infrared (1063 and 1545 nm) and is designed for field experiments i.e. outdoor applications.SALCA is tested for forest understory mapping by comparing raw and normalized intensities.A derived normalized ratio index is calculated from both wavelengths and shows promising results for applications such as analysis of biophysical variables of forests and water stress monitoring.Hakala et al. (2012) differentiate tree trunks from tree tops of norway spruce (Picea abies) with a TLS FWF hyperspectral LiDAR operating in the wavelength domain from 480 nm to 1,000 nm.In a laboratory test they compare the reflectance behavior to passive spectrometer measurements showing good correlation of backscatter reflectance curves between sensors.Additionally high separability between trunk and tree crown is given.Douglas et al. (2015) apply the so called dual wavelength Echidna LiDAR (DWEL), which is a TLS scanner operating with 1,064 nm and 1,548 nm wavelengths.They use the system for an outdoor experiment separating canopy (leafs) from tree trunks calculating a normalized difference of channel intensity.Shi et al. (2015) develop an intensity calibration strategy for multiwavelength LiDAR operating with four wavelengths (556, 670, 700, and 780 nm) considering incidence angle and surface roughness improving classification using vegetation indices.The approach is a laboratory setup for scanning single leafs.Puttonen et al. (2015) present a study on an outdoor scanning experiment with different high vegetation species.They calculate several vegetation indices relating them to leaf area index, chlorophyll in leafs and biomass estimates.The hyperspectral LiDAR emits white laser pulses in the wavelength domain from 420 to 2,100 nm with a channel width of 20 nm.In a time series experiment of 26h, selected channels (545.5, 641.2, 675.0, 711.0, 741.5, 778.4, 978.0) are used to investigate the influence of sunlight on spectral and geometrical properties and the normalized difference vegetation index (NDVI).
First operational systems using different wavelengths in ALS were developed for coastal areas and surface mapping in shallow water.One example is the Scanning Hydrographic Operational Airborne LiDAR Survey (SHOALS) System (e.g.Irish and Lillycrop, 1999) combining lasers operating in infrared 1,064 nm and blue-green 532 nm.The combination of a red and green laser with different beam divergence allows the measurement of the water surface and the penetration of the water surface to a certain depth.The potential of multispectral full waveform ALS for forest structure analysis has been investigated in a detailed study by Morsdorf et al. (2009) simulating waveforms, NDVI, backscatter energy distributions and comparing these to canopy volume distributions along tree canopies making use of virtual 3D tree models.A comprehensive overview on the concept of multiwavelength ALS has been presented by Pfennigbauer and Ullrich (2011).Based on this study, one of the first multiwavelength LiDAR ALS experiments for topographic land surface and object differentiation has been conducted by Briese et al. (2012Briese et al. ( , 2013)).They assemble FWF ALS point clouds from three independent flight missions using RIEGL scanners operating in different wavelengths (532 nm, 1,064 nm, 1,550 nm).After radiometric calibration composites of amplitudes are generated.In Briese et al. (2013) multiwavelength LiDAR amplitudes are used for archeological prospection on grass vegetated terrain.A classification study of land cover (roads and gravel, bare soil, low vegetation, high vegetation, roofs, and water bodies) is presented by Wang et al. (2014).They combine surveys from an Optech ALTM Pegasus HD400 (1,064 nm) and a Riegl LMS-Q680i (1,550 nm).The features amplitude, echo width and surface height are the input for a support vector machine classification.
Within this study, we use a data set from a prototype of the Optech Titan ALTM system, the first commercial system combining three wavelengths (532, 1,064 and 1,550 nm) in a single sensor system.After a short introduction to the data set, we describe a first approach to process such data and present results on spectral pattern analysis and point cloud classification opportunities for topographic mapping and land cover classification.
The main aims of this work are: The exploration of characteristic surface class spectral signatures and LiDAR specific implications.
The evaluation of the suitability of the data set for state-of-theart geometrical LiDAR classifications.
The evaluation of spectral signature based classifications as improvements to conventional, only geometry based LiDAR classifications.

DATASET
A single flight strip of 3 km 2 (4 km in length and 600 m wide) covering the city of Oshawa (Ontario, Canada) was captured with a prototype of the airborne multispectral LiDAR System Optech Titan in September 2014.The system works with three independent active imaging channels at 1,550 (C1), 1,064 (C2) and 532 nm (C3).Laser pulses of the three channels are sent to the same oscillating mirror with individual deflecting angles for each channel.This leads to separate optical light paths for each channel: C2 0° nadir looking, C1 3.5° forward looking, and C3 7° forward looking.The system has a maximum laser measurement rate of 300 kHz per channel (900 kHz in total) and provides three independent point clouds, one for each channel.
The data was acquired during leave-on conditions to optimize the geometrical properties of the system, operating at 980 m above ground with a field of view (FOV) of +/-20 degree.The point clouds are geometrically registered, but provide only raw intensity values.The three point clouds have the following properties: -C1: 11,573,302 points with an average point density of 4 pts/m 2 (stddev.4 pts/m 2 ) and an average point distance of 0.32 m (stddev.0.27 m) -C2: 12,262,554 points with an average point density of 5 pts/m 2 (stddev.4 pts/m 2 ) and an average point distance of 0.32 m (stddev.0.29 m) -C3: 9,966,192 points with an average point density of 4 pts/m 2 (stddev.4 pts/m 2 ) and an average point distance of 0.29 m (stddev.0.14 m) The point cloud which resulted from merging the three channels has 33,802,048 points with an average point density of 11 pts/m 2 (stddev.8 pts/m 2 ) and an average point distance of 0.20 m (stddev.0.21 m).

METHODS
The methodical approach is adapted to meet the limited information available for this prototype data set.As no sufficient meta data (e.g.flight trajectory or atmospheric conditions) and ground truth data are available, no detailed analysis on beam incidence angle effects as presented by Jutzi and Gross (2010) or radiometric calibration approaches (e.g.Briese et al., 2012) could be conducted.All processing is done with the System for Automated Geoscientific Analysis (SAGA, Conrad et al., 2015) and the extension Laserdata Information System (LIS, Laserdata, 2015, Rieg et al., 2014).
In order to avoid gridding and to preserve the 3-dimensionality of the input data the three independent point clouds of each channel are merged into a single point cloud by a nearest neighbor approach.As the laser beams of channel C1 and C3 are tilted by +3.5 and +7 degree from the nadir direction, it is not possible to use the GPS time to assign the intensity values.Therefore a tool was developed, which merges the channels geometrically by assigning to each point the intensity value of the nearest neighbor point found in each other channel point cloud.This results in a point cloud containing all points of the independent point clouds, with each point having the three channel intensities attached (its own intensity and two assigned intensities).The nearest neighbor search is limited by a maximum distance of 1 m to a neighboring point in order to prevent the grouping of points which are located on different objects.This should help to avoid mixed signal effects.In case there is no point found within 1 m distance, we assume a dropout in this channel and assign an intensity value of zero (complete absorption).
Besides point cloud merging and intensity grouping, the tool performs additional point feature calculations and attaches the results as attributes to each point.This includes the channel number of the point itself and a binary encoding of the availability of each channel within the search distance.The latter allows analyzing the spatial occurrence of drop-outs later.Additionally, the geometric distances between the points used for each channel combination are calculated.For this calculation also points outside the maximum nearest neighbor search distance of 1 m, and thus classified as drop-outs, are included.These distances are used to calculate a mean distance, which is used as quality measure of the channel grouping.
Due to the limited availability of meta data, we apply basic image processing strategies (e.g.contrast stretching, Lillesand and Kiefer, 1999) to the raw intensity values.For each channel, the 99.99th percentile is computed and used as intensity cutoff threshold.This way, outliers from specular reflections are eliminated.Visual inspection of these points revealed that such outliers occur on moving objects like cars and other man-made objects like building facades.Finally, a histogram stretch to twice the standard deviation is computed for each channel, resulting in value ranges from 0 to 255.Finally, the point cloud is classified into the main classes 'unsealed ground', 'sealed ground', 'buildings', 'mid vegetation', 'high vegetation', 'water surface' and 'water body / bottom' and for some small training areas into the subclasses 'green grass', 'dried up grass', 'sand / bare soil', 'wetlands', 'darker asphalt' and 'lighter asphalt'.The classification is done by a two stage approach consisting of automatic and semiautomatic classification.First the point cloud is automatically classified into 'ground', 'buildings', and 'mid vegetation' and 'high vegetation' classes by a "geometrical classification" using information about x,y,z coordinates and neighborhood exclusively.Ground points are classified using a hybrid approach of progressive TIN densification (Axelsson, 2000) and RANSAC-based point cloud segmentation.The latter is also used (together with other point cloud features such as eigenvalue based omnivariance, Mallet et al., 2011) to differentiate the remaining non-ground points into buildings and vegetation.For the accuracy assessment of the automatic classification, two test areas, each 600 x 600 m, were selected and manually revised.In a second semi-automatic step, the ground class is further subclassified by introducing three additional classes: (i) 'sealed ground' (e.g.roads), (ii) 'water surface', and (iii) 'water body / bottom'.For further analysis of small training samples the heterogeneous classes of 'sealed' and 'unsealed ground' were semi-automatically classified into 'green grass', 'dried up grass', 'sand / bare soil', 'wetlands', 'darker asphalt' and 'lighter asphalt'.
For each channel, histograms showing the intensity value distribution for the seven classes are calculated and analyzed in order to evaluate the potential of extending the geometrical classification approaches currently available for airborne LiDAR mapping.
Because of their heterogeneity, the ground classes 'unsealed' and 'sealed' are examined in greater detail.Therefore, training areas composed of different materials are manually selected and again histograms are calculated for each sub-class.The channel peak values of each sub-class are then used for a supervised classification by pattern matching based on the Mahalanobis distance.

RESULTS
For the merged channel datasets we computed false color composite point clouds as shown in Fig. 1.Based on the manual classification of the datasets, we looked at the 8 bit scaled channel histograms and the pseudo NDVI for the main classes.
As the majority of the data showed only a spread of 201, 167 and 133 integer values for the channels C1, C2 and C3, the 8 bit scaled histograms contain data gaps (Fig. 2,3 and 5,6).

Spectral patterns of main classes
Looking at the two subareas, a set of separated peaks is visible in the histograms of the three channels.This indicates heterogeneity of spectral patterns within one object class (Fig. 2 and 3).

Spectral patterns of ground sub-class
For manually selected small training areas within the sealed and unsealed ground class we created also channel histograms for the subclasses 'green grass', 'dried up grass','sand / bare soil', 'wetlands', 'darker asphalt', and 'lighter asphalt' (Fig. 5 and 6).This was supported by the inspection of ortho images from web mapping services.
For the sealed classes 'lighter asphalt' and 'darker asphalt' a clear separation is possible due to intensity differences in all channels.The 'darker asphalt' (~50, 25, 70) shows lower intensities than the 'lighter asphalt' (~110, 50, 140).Looking at single channels, the discrimination from unsealed classes is difficult.
In C1 'darker asphalt' shows the same peak as 'wetlands' and 'lighter asphalt' is very similar to 'dried up grass'.In C2 'lighter asphalt' shows the same peak as 'wetlands' and 'dried up grass'.In C1 'darker asphalt' shows the same peak as 'dried up grass' and 'lighter asphalt' is very similar to 'sand / bare soil'.
The unsealed group of classes shows fluent transitions from 'wetlands' (~50,50,40) to 'green grass ' (170,200,110), to 'dried up grass' (110,60,75) and 'sand / bare soil' (210,130,140).This is visible in the histogram of 'green grass' including different species and all states of health.In theory, the moisture and vegetation density is increasing with channel intensity.Table 2. Confusion matrix of the automatic geometrical classification

Geometrical Classification
Looking at the potential of multispectral LiDAR processing, it is important to evaluate the geometrical suitability of the data set for conventional processing workflows.Due to the setup of the technical system it can be observed that e.g. the standard deviation of the point densities doubles with the merge of the three individual channel point clouds (see section 2).Therefore we tested, whether the geometrical representation meets the quality of state-of-the-art single wavelength data sets.
Evaluating the confusion matrix (Table 2) of the automatic geometrical classification shows sufficient results, compared to single wavelength systems.The overall accuracy is 0.99.Higher numbers of wrong classifications (User's Accuracy of 0.88, Table 3) can be observed for the 'undefined' class.
Here, 3,335 points of 76,086 reference points are falsely labelled as 'high vegetation', this is mainly due to cabels that are not addressed in detail within the applied classification approach.3,121 points are wrongly labelled as ground, mainly due to moving objects and other smaller structures.A remarkable error is the wrong classification of 2,070 'mid vegetation' points as buildings, which is due to the planar appearance of hedges in the data set.

DISCUSSION
The class specific spectral patterns observed in the dataset show similarities to those found in optical imagery.One example is given by the low reflectance values of moist surfaces and increasing reflectance values with increasing dryness of surfaces.However, the representation as a point cloud and the acquisition technique include additional signatures interfering with the spectral ones.This can be explained especially for the classes 'high' and 'mid vegetation'.

Multi return effects
As known from optical imagery, lower reflectance values for the 532 and 1,550 nm wavelength are expected for vegetation cover, while higher reflectance values are expected for the wavelength of 1,064 nm (Lillesand and Kiefer, 1999).This is the case for the class 'green grass' showing good agreement with this assumption.However, as a ground class, it shows a high planarity and homogeneity.With increasing roughness of the surfaces and decreasing size of single targets (getting smaller than the laser beam cross-section), the geometrical effect is intensified.For 'high vegetation' and 'mid vegetation' the classical spectral pattern of vegetation is altered by splitting the emitted signal into multiple echoes and the convolution of the original signal at rough surfaces.Thereby the signal intensity becomes distributed over a varying number of returns, leading to lowered intensities for each single point of a multi return shot.In high vegetation, the single points represent only small targets such as leaf clumps and thus are supporting the effect of lowered intensities.In combination, the spectral and geometrical properties decrease the intensity values of high vegetation in all channels, making the separation from e.g.buildings difficult.However, as all channels are affected by return modification effects, the normalized ratio of the NDVI shows a good potential for the discrimination of vegetation and sealed surfaces.
In comparison to 'high vegetation' the signature of 'mid vegetation' shows slightly higher reflectance values.This can also be explained by geometrical effects.As most of the 'mid vegetation' in the given residential area are well trimmed hedges showing dense and almost planar surfaces, the reflected signals are less effected by signal split-up and modification.
Besides, using e.g. the NDVI for the classification of these features can lead to improved results compared to pure geometrical classifications, where a planar surface might be interpreted as man-made object (see section 4.3).
Multi echoes also affect the distinction of ground classes.Ground points below vegetation cover or close to buildings share their corresponding shot intensities with additional points in the canopy or on a building.Thus, the expected range of intensity values cannot be reached for these points and leads to altered spectral patterns.The most obvious example is the false classification of multi echoes from 'green grass' into the rarely reflecting 'wetlands' class.To manage the effects of signal modification, a strategy for return number normalization should be developed.

Drop-out effects
A typical feature introduced during the geometrical channel merging approach are points for which no corresponding point in another channel can be found.Here, we call these missing points a 'drop-out' in this channel.Drop-outs can have several reasons (Höfle et al., 2009), including i) complete absorption, ii) too low reflection values because of strong scattering, iii) specular reflection away from the sensor and iv) no laser shot of this channel hitting a target within the search distance.The most obvious example is the channel pattern for 'water body / bottom', showing strong agreement with the principles used for bathymetric LiDAR.Here, the point distribution of C3, which is able to penetrate shallow water bodies, allows detecting water points in the data set.Drop-outs can also be observed for 'mid vegetation' and 'high vegetation'.Besides a complete absorption because of the material properties, the complex geometrical distribution of targets and the target sizes can result in strong scattering.The proportion of light reflected by smaller targets might become too small to be detected by the sensor.C3 shows a higher probability for drop-outs, which can be explained by wavelength and the larger footprint compared to the other channels.The emitted energy becomes distributed over a larger cross-sectional area.Thus, the ratio of target area and cross-sectional area becomes smaller, leading to even smaller proportions reflected from a small scatterer.Because of the tilt (looking forward) of channel C1 and C3, the lines of sight of the three emitted laser beams are different.Larger and smaller targets are hit by chance by the three lasers, causing random drop-outs on each channel.

Implications for classification approaches
As shown in section 4, a pattern based classification of ground classes, such as 'green grass', 'wetlands' or 'asphalt' provide promising results.Combined with conventional geometrical classification workflows it can improve the information extraction from LiDAR data.However, multi echo effects have to be taken into account in order to stabilize classification results.Without consideration, points with higher return numbers are prone to classification errors.As these points contain underestimated reflectance values, they tend to be grouped to classes with low reflectance values such as 'wetlands' or 'darker asphalt'.Drop-outs are an advantage for classification approaches, giving additional information onto object heterogeneities and the occurrence of water surfaces.

CONCLUSION
In the presented paper an explorative analysis of class specific spectral signatures is conducted for the first commercial multispectral LiDAR system Optech Titan.Spectral patterns are used for a classification approach and evaluated.Typical spectral patterns comparable to optical imagery could be observed for homogeneous and planar surfaces.LiDAR specific signatures introduced through multiple echo detection and the separate optical light paths of the channels, lead to underestimations of channel intensities and channel drop-outs.While drop-outs might improve the possibility to detect water and vegetation, intensity underestimations are confusing class signatures and classification results.Here, a strategy for return number normalization of the intensities has to be developed in the future.
We also looked at the geometrical suitability of the data set for conventional LiDAR classification workflows.Although the data acquisition technique of the system has to manage the timeand orientation-shifted scan pattern, the geometrical characteristics of the given dataset are sufficient for geometrybased processing workflows.Thus, while additionally providing spectral information, the system maintains the geometrical characteristics of single wavelength LiDAR systems.
As the test data set does not provide a complete set of meta data, all the described processing is done on raw intensity data.In future, capturing a dataset including the full range of input meta data and ground truth data sets, would allow for a complete radiometric calibration.In practise, most of the information necessary for a radiometric calibration will not be available in most cases.Therefore approaches using raw intensity values in a stable manner would be welcome.
The simultaneous recording of geometrical and spectral information with an active sensor has advantages compared to classical monochromatic laser scanning systems and optical imagery: -almost no time gap between the acquisition of geometrical information per channel (potential for the classification of e.g.moving objects).
-almost no time gap between geometrical and spectral information, which is an important improvement to experimental approaches, such as the merging of separate flight missions (e.g.Briese et al., 2012Briese et al., , 2013)), or the merging of ortho photos to the point cloud acquired at different time steps.Compared to optical imagery, it is also influenced by comparable effects such as mixed pixels (geometrical channel grouping) or sun illumination (surface roughness and laser beam interaction).

Figure 2 .
Figure 2. Class specific 8 bit scaled channel histograms for 1,550 nm and 1,064 nm wavelength laser signals (counts/bin)The 'unsealed ground' includes a variety of sub-classes representing ground vegetation of varying health and density.This is most obvious in C2 showing multiple peaks.C2 shows the best contrast to the 'sealed ground' class, while C1 and C3 show distributions very similar to the 'sealed ground' distributions.The pseudo NDVI also shows a good discrimination between 'sealed ground' (low NDVI) and 'unsealed ground' (high NDVI).The suitability of the different channels for the separation of 'sealed ground' and 'unsealed ground' is shown in Fig.4.

Figure 3 .
Figure 3. Class specific 8 bit scaled channel histograms for 532 nm wavelength laser signals and pseudo NDVI (counts/bin)Combining all three channel intensities for a pattern matching classification (Fig.7) leads to a good agreement with the visual inspection obtained from the false colour coded point cloud.Besides, the pseudo NDVI as a two channel ratio shows a good contrast between sealed and unsealed surface classes.

Figure 4 .
Figure 4. Contrast between sealed and unsealed classes: a) low contrast for channel C1; b) high contrast for channel C2; c) moderate contrast for channel C3 (inverse to C1 and C2); d) high contrast for the pseudo NDVI

Table 1 .
Confusion matrix of the pattern based classification of sealed and unsealed Table 1).