Object-based analysis of multispectral airborne laser scanner data for land cover classification and map updating

ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Contents lists available at ScienceDirect

ISPRS Journal of Photogrammetry and Remote Sensing journal homepage: www.elsevier.com/locate/isprsjprs

Object-based analysis of multispectral airborne laser scanner data for land cover classification and map updating Leena Matikainen ⇑, Kirsi Karila, Juha Hyyppä, Paula Litkey, Eetu Puttonen, Eero Ahokas Finnish Geospatial Research Institute FGI at the National Land Survey of Finland, Centre of Excellence in Laser Scanning Research, Geodeetinrinne 2, FI-02430 Masala, Finland

a r t i c l e

i n f o

Article history: Received 21 November 2016 Received in revised form 7 April 2017 Accepted 10 April 2017

Keywords: Laser scanning Lidar Multispectral Land cover Updating Change detection

a b s t r a c t During the last 20 years, airborne laser scanning (ALS), often combined with passive multispectral information from aerial images, has shown its high feasibility for automated mapping processes. The main benefits have been achieved in the mapping of elevated objects such as buildings and trees. Recently, the first multispectral airborne laser scanners have been launched, and active multispectral information is for the first time available for 3D ALS point clouds from a single sensor. This article discusses the potential of this new technology in map updating, especially in automated object-based land cover classification and change detection in a suburban area. For our study, Optech Titan multispectral ALS data over a suburban area in Finland were acquired. Results from an object-based random forests analysis suggest that the multispectral ALS data are very useful for land cover classification, considering both elevated classes and ground-level classes. The overall accuracy of the land cover classification results with six classes was 96% compared with validation points. The classes under study included building, tree, asphalt, gravel, rocky area and low vegetation. Compared to classification of single-channel data, the main improvements were achieved for ground-level classes. According to feature importance analyses, multispectral intensity features based on several channels were more useful than those based on one channel. Automatic change detection for buildings and roads was also demonstrated by utilising the new multispectral ALS data in combination with old map vectors. In change detection of buildings, an old digital surface model (DSM) based on single-channel ALS data was also used. Overall, our analyses suggest that the new data have high potential for further increasing the automation level in mapping. Unlike passive aerial imaging commonly used in mapping, the multispectral ALS technology is independent of external illumination conditions, and there are no shadows on intensity images produced from the data. These are significant advantages in developing automated classification and change detection procedures. Ó 2017 The Authors. Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction 1.1. Background on multispectral laser scanning applications Airborne laser scanning (ALS) has shown its high feasibility for automated mapping processes. During the last 20 years, the accurate 3D data, often combined with spectral information from digital aerial images, have allowed the development of methods for automated object extraction and change detection. In particular, important advancements have been achieved in the mapping of elevated objects such as buildings and trees (see, for example, ⇑ Corresponding author. E-mail addresses: [email protected] (L. Matikainen), [email protected] (K. Karila), [email protected] (J. Hyyppä), [email protected] (P. Litkey), eetu. [email protected] (E. Puttonen), [email protected] (E. Ahokas).

Hug, 1997; Rottensteiner et al., 2007; Brenner, 2010; Guo et al., 2011). Some researchers have concentrated on roads (e.g., Hu et al., 2004; Clode et al., 2007), and numerous studies related to more general land cover mapping have also been carried out (Yan et al., 2015). The role of laser intensity has been relatively small in the development. Intensity information, increasingly also from fullwaveform data, has been used in many classification studies (e.g., Hug, 1997; Guo et al., 2011), but generally the geometric information of the laser scanner data has been more important. In 2008, EuroSDR (European Spatial Data Research) initiated a project ‘Radiometric Calibration of ALS Intensity’ led by the Finnish Geospatial Research Institute FGI and TU Wien in order to increase the awareness of intensity calibration. It was expected that someday multispectral ALS would be available and then the intensity of

http://dx.doi.org/10.1016/j.isprsjprs.2017.04.005 0924-2716/Ó 2017 The Authors. Published by Elsevier B.V. on behalf of International Society for Photogrammetry and Remote Sensing, Inc. (ISPRS). This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

ALS will be of high value. Thus, earlier intensity studies (e.g., Luzum et al., 2004; Kaasalainen et al., 2005; Ahokas et al., 2006; Wagner et al., 2006; Höfle and Pfeifer, 2007; Wagner, 2010) showing the concepts for radiometric calibration of ALS intensity have been preparations for multispectral ALS. Laser scanning systems providing geometric and multispectral information simultaneously provide interesting possibilities for further development of automated object detection and change detection methods. Experiences on this field have been obtained, for example, from studies with a hyperspectral laser scanner developed at the FGI (Chen et al., 2010; Hakala et al., 2012). The system is a terrestrial one and has been successfully applied for vegetation analyses and separation of man-made targets from natural ones (e.g., Puttonen et al., 2015). Multiple wavelengths also allow detection of detailed spatiotemporal changes in vegetation (e.g., Hakala et al., 2015; Puttonen et al., 2016). Pfennigbauer and Ullrich (2011) discussed the development and potential of multi-wavelength ALS systems. Briese et al. (2013) realised such an approach by using three separate ALS systems. The study concentrated on radiometric calibration of intensity data from an archaeological study area. Wang et al. (2014) used two separate ALS systems to acquire dual-wavelength data for classifying land cover. They found that the use of dualwavelength data can substantially improve classification accuracy compared to single-wavelength data. More extensive discussions and reference lists on these earlier developments of multispectral laser scanning can be found, for example, in Vauhkonen et al. (2013), Puttonen et al. (2015), Wichmann et al. (2015) and Eitel et al. (2016). The first operational multispectral ALS system was launched by Teledyne Optech (Ontario, Canada) in late 2014 with the product name Titan. With this scanner, active multispectral information is for the first time available for 3D ALS point clouds from a single sensor. The channels of the sensor are: infrared 1550 nm (channel 1, Ch1), near-infrared 1064 nm (Ch2) and green 532 nm (Ch3). Each of the channels produces a separate point cloud. The look directions of the channels are: Ch1: 3.5° forward, Ch2: nadir, Ch3: 7° forward. Due to the separate point clouds, multispectral intensity values are not originally available for single points. Some preprocessing is thus needed before the multispectral information can be utilised. The first studies based on Optech Titan data have been published, most of them in conference proceedings and based on Optech sample data acquired from Ontario in 2014. Wichmann et al. (2015) evaluated the potential of Titan sample data for topographic mapping and land cover classification. They analysed spectral patterns of land cover classes and carried out a classification test. The analysis was point-based and used a single point cloud that was obtained by merging three separate point clouds from the three channels by a nearest neighbour approach. Wichmann et al. (2015) concluded that the data are suitable for conventional geometrical classification and that additional classes such as sealed and unsealed ground can be separated with high classification accuracies using the intensity information. They also discussed multi-return and drop-out effects, drop-outs referring to cases where returns are not obtained in all channels. For example, the spectral signature of objects such as trees becomes biased due to multi-return effects. Bakuła (2015) created digital surface models (DSMs) and digital terrain models (DTMs) separately from the different channels and analysed them. Differences occurred in vegetated areas and building edges due to various distributions of points. For DTMs, the differences were small. The potential of the data for land cover mapping was also discussed. A further land cover classification experiment was presented in Bakuła et al. (2016). An overall accuracy of about 91% was achieved in classifying six classes: water, sand and gravel, concrete and asphalt, low

299

vegetation, trees, buildings. This classification was raster-based and used spectral data, elevation, and textural data. Morsy et al. (2016) applied three normalized difference feature indices calculated from the Titan channels for land–water and vegetation– built-up separation. Classification experiments with data from the same area were also reported by Miller et al. (2016). Zou et al. (2016) used Titan data acquired from Ontario in 2015. They applied an object-based classification approach. Segmentation was based on the multispectral intensity data and classification rules were constructed by using four indices based on intensity, return count and elevation information. The overall accuracy achieved in classification with nine classes was about 92%. The classes included water bodies, bare soil, lawn, road, building, low vegetation, medium vegetation, high vegetation and power line. Hopkinson et al. (2016) investigated the characteristics of Titan data and multisensor data acquired with the same wavelengths in a forested test site in Canada. Differences in proportions of returns at ground level, vertical foliage distributions and gap probability were observed across the wavelengths. A multispectral classification test with rasterised intensity data was also carried out. Eight classes were included in the Titan data classification: corn, gravel, hay, asphalt, larch, hardwood, pine/spruce and sand. The highest overall accuracy from training site accuracy statistics was about 88%. The overall accuracies achieved with the multisensor data were lower, although the number of classes in that case was only five. Hopkinson et al. (2016) also found that the overlap in single-channel intensity values between different land cover classes was so high that it prevents accurate classification from one channel data. Fernandez-Diaz et al. (2016) discussed the performance of a Titan system on the basis of test campaigns during two years. Among many other topics, results of a land cover classification test from a university campus area in Houston, Texas, were reported briefly. Intensity images and structural images derived from the elevation data and return counts were used in supervised classification. The resolution of the rasterised data was 2 m, and five classes were included in the study: grass/lawn, road/parking, trees and short vegetation, commercial buildings, and residential buildings. The best overall accuracy was about 90%, and it was obtained when structural images and two intensity channels (Ch2 and Ch3) were used. In most of the Titan studies discussed above, radiometric calibration of intensity data has not been reported. Morsy et al. (2016) mentioned an intensity correction mainly for the range and the scan angle in land–water mapping. Fernandez-Diaz et al. (2016) made an intensity correction by normalising raw intensities by range. Some other recent studies with multiple wavelength ALS data can also be found. Leigh and Magruder (2016) investigated the use of dual-wavelength, full-waveform data for surface classification. The data were acquired with Leica Chiroptera sensor (Leica Geosystems AG, Switzerland), which has been planned for shallow water and coastal mapping. The sensor has green and near-infrared channels. Leigh and Magruder (2016) used a voxelisation approach to utilise the full-waveform information and to extract features for random forest classification. The best overall accuracy was about 81%, and it was achieved when both wavelengths were used. The classes under study included pavement, cleared paths, grass, scrub (two classes), rock, pine, and cypress. 1.2. Objectives and contribution of our study Our article discusses the potential of the new single-sensor multispectral ALS technology in land cover classification and map updating. Our main interest is in nationwide mapping and updating of databases in the future. In addition to 3D modelling of objects, an important application area of ALS data has been automated change detection, especially for buildings. Despite promis-

300

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

ing results achieved with ALS and/or aerial image data in research projects (e.g., Holland et al., 2008; Champion et al., 2009), the current updating processes in mapping organisations are still largely based on visual analysis and manual digitising. There is thus much potential for further developments. Methods developed for change detection of buildings from ALS data can be roughly divided into two groups based on the input data available. The first group uses recently acquired ALS data, possibly combined with new image data, and aims to detect changes compared to an existing building map (e.g., Matikainen et al., 2003, 2010; Vosselman et al., 2005; Malpica et al., 2013). The second approach assumes that both old and new ALS datasets are available and the change detection method can rely on the comparison of these (e.g., Murakami et al., 1999; Richter et al., 2013; Teo and Shih, 2013). In the first group of methods, the methods are directly linked to the existing building objects on a map, and change information is normally obtained for these objects. In the second group, the link between the change detection method and existing building maps or models is often weaker. For our study, Optech Titan multispectral ALS data over a large suburban area in Finland were acquired. The data were acquired with an operational system, i.e., we did not use sample data provided by Optech. The objective of our study is twofold. Firstly, we studied how the novel multispectral information can benefit land cover mapping in a suburban area. An object-based classification approach and the random forests method (Breiman, 2001) were used. The importance of various features for classification was investigated and the accuracy of the results was analysed. For comparison, a classification test with single-channel data was also carried out. The motivation for this study is to provide basic information for development of automated mapping procedures in the future. Secondly, our objective is to demonstrate the potential of the new data for such procedures with two examples: buildings and roads. The test with buildings is an extension to the previous change detection research and illustrates an approach that uses the new multispectral ALS data in combination with old single-channel ALS data and existing building vectors. The test with roads exploits the specific capability of the multispectral ALS data in discriminating ground-level classes. The preliminary results of our study were presented in conference articles (Matikainen et al., 2016; Ahokas et al., 2016). An analysis concentrating on the potential of the multispectral ALS data for road classification can be found in Karila et al. (2017). The present article includes several extensions and improvements such as a larger study area, larger number of multispectral features analysed and quality analysis of the land cover classification results with independent validation data. According to our knowledge, it is the first study in international journals focusing on land cover classification and map updating with multispectral ALS data in a suburban area.

2. Study area and data 2.1. Study area and reference data The study area is a suburban area located in Espoonlahti in the city of Espoo, near Helsinki on the southern coast of Finland (centre point approx. 60°90 1800 N, 24°380 2400 E). It comprises about 10 km2 of land area (Fig. 1). In this study we concentrated on mapping of the land area, and sea was excluded by using a water mask derived from topographic map data. The change detection applications were demonstrated in a smaller subarea where a new residential area has been built during the last 15 years. The study area includes a land cover classification test field established previously (see Matikainen and Karila, 2011). The location and land cover of all test field points were now checked, and if

needed, corrected to correspond to the current situation. This updating of the points was based on visual interpretation of the Titan data and open datasets available from the City of Espoo (http://kartat.espoo.fi/ims), the National Land Survey of Finland (NLS) (http://www.maanmittauslaitos.fi/aineistot-palvelut/latauspalvelut/avoimien-aineistojen-tiedostopalvelu), Google Maps (Google Inc., Mountain View, CA, USA) (https://maps.google.com) and Google Street View. Field checks were used to support the analysis and to check uncertain cases. The updated test field points were then used as reference data in the study. A set of 332 points was used for training and a separate set of 566 points was used for validation of the classification results. Fig. 1 illustrates the study area and distribution of the reference points. The training points were originally collected by using a 50 m  50 m grid and the validation points with a 100 m  100 m grid so that one point is located within each grid cell. The points are located in homogeneous places, i.e., not on borders between different classes. The attributes of the reference points basically include object types and surface types. For the present study, this information was generalised to the classes presented in Table 1. Due to the use of separate training and validation areas, there were minor differences in the contents of the land cover classes in the training and validation sets. For example, there were vegetable garden points only in the training area and most of the forest floor points were located in the validation area. However, due to the use of generalised classes such as low vegetation, this is not likely to be a significant problem in the classification. The use of separate training and validation areas can give more realistic information on the performance of the classification in a new area. Grass and meadows were combined into the low vegetation class because they are very difficult to differentiate strictly in our suburban study area even in field checks. In evaluating change detection results, map data from 2015 and current aerial imagery were used. These were downloaded from the NLS link given above. 2.2. Multispectral ALS data Multispectral Optech Titan data over the study area were acquired on 21 August 2015 in cooperation with TerraTec Oy (Helsinki, Finland). The data were acquired from a fixed-wing aircraft, and the flying altitude was approx. 650 m above sea level. The dataset includes the three Titan channels available as separate point clouds. The point densities in our preprocessed dataset over land areas are approximately the following: Ch1: 9 points/m2; Ch2: 9 points/m2; Ch3: 8 points/m2. The first part of preprocessing was carried out by TerraTec Oy and included some basic processing steps such as geoid correction and strip matching. Further processing was carried out at the FGI. Radiometric calibration of ALS intensity is an important prerequisite for successful classification (Ahokas et al., 2006; Höfle and Pfeifer, 2007; Kaasalainen et al., 2011). In this study, we used relative radiometric calibration. As expected, our preliminary analyses with uncorrected data showed that intensity values were highest in the middle of a flight line and decreased as the distance from the scanner increased. A range correction was applied to correct this effect:

icorr ¼ i 

R2i R2ref

;

ð1Þ

where icorr is the corrected intensity, i is the original intensity, Ri is the distance from the scanner to the scanned point, and Rref is the flying altitude (650 m). After the radiometric correction, functions of TerraScan software (Terrasolid Ltd., Helsinki, Finland) were used to cut overlap points and remove erroneously high or low points.

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

301

Fig. 1. Espoonlahti study area with Titan intensity image in the background (Red: Ch1; Green: Ch2; Blue: Ch3). The black rectangle shows the subarea used in change detection tests. The water mask contains data from the National Land Survey of Finland Topographic Database 2015.

Finally, three-channel intensity images and DSMs in raster format were created. We selected object-based analysis of raster format data for this study because it is well suited for producing land cover classifications of large areas. In addition, image format is a natural way to present and analyse intensity information, which was the main interest with multispectral ALS data in our study. Objectbased analysis of intensity data diminishes possible problems related to individual intensity values and helps to give an overall understanding of the capabilities of the new data. In the intensity images, pixel values in each channel represent the average intensity value of first pulse and only pulse laser points divided by 100. For creating the DSMs, a combined point cloud from all channels was used. Maximum and minimum DSMs, presenting the highest or lowest height within the pixel, respectively, were created. Interpolation was used to obtain intensity and height values for pixels without laser points. The processing was carried out with the TerraScan software. The pixel size of the intensity images used in the present study was 0.2 m. The pixel size of the DSMs was 1 m. For the intensity images, the small pixel size was preferred to better represent small details such as narrow roads.

For the DSMs, the larger pixel size appeared advantageous to avoid some noise in building edges and to remove most of the vegetation from the minimum DSMs. Fig. 2 shows a three-channel intensity image for a part of our study area. In the figure, the intensity image can be compared with a conventional aerial ortho image. A DTM was also produced from last pulse ALS points (included only and last of many pulses). The ground classification of the TerraScan software was used to detect ground points. This method is originally based on a filtering method developed by Axelsson (2000). A raster DTM with the pixel size of 1 m was then created. The process was fully automatic, i.e., manual corrections, which are typically used in operational work, were not used in our case. This meant, for example, that bridges were included in the ground surface in the DTM.

2.3. Old map data and old DSM Old map data were used to demonstrate change detection for buildings and roads. In the case of buildings, a building map repre-

302

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Table 1 Reference points used for training and validation of the land cover classification. Class High objects Building

Number of training points

Number of validation points

Comments

88

130

Different types of buildings (e.g., high-rise residential, low-rise residential, industrial, public) with different roof materials and colours Coniferous and deciduous trees in forest, in gardens and along roads

Tree Low objects Asphalt Gravel Rocky area Low vegetation

83

142

62 15 16 68

158 16 44 76

Total

332

566

In addition to 213 asphalt points, 7 road or parking place points with tile surface were included. Soft, non-vegetated surfaces with different grain sizes (e.g., small roads, sports fields, beaches) Rocky areas with bare or slightly vegetated surface (typically some moss or patchy grass) Includes grass (42 points in total), meadow (66), forest floor (20), vegetable gardens (9), and low bushes (7)

Fig. 2. Comparison of an aerial ortho image (left) and Titan intensity image (right) in a 300 m  300 m subarea. Ortho image Ó National Land Survey of Finland 2013. It should be noted that the acquisition time of the images was different. In the Titan data, deciduous trees are in full leaf, while the aerial ortho image is from leaf-off conditions in early spring.

senting the situation in 2005 was used as an old map. The data were originally obtained from the City of Espoo and edited at the FGI to correspond to the situation in 2005 (Matikainen et al., 2010). Old roads were obtained from an old version of the NLS Topographic Database from 2000. An old DSM based on Optech Airborne Laser Terrain Mapper (ALTM) 3100 ALS data acquired in July 2005 was also used as input data in the change detection of buildings. The DSM represented the minimum height values for 0.3 m  0.3 m pixels (Matikainen et al., 2010). It should be noted that the height values in the old DSM were in the old Finnish N60 height system, while the values in the new DSMs were in the newer N2000 height system. The height difference between these systems is about 25 cm in our study area. Transformation was not considered necessary because the change detection method used only a rough comparison between the DSMs and the systematic difference could be taken into account in selecting threshold values for the comparison. 3. Methods 3.1. Land cover classification method 3.1.1. Segmentation and calculation of features The DSM and intensity image data were first segmented into homogeneous regions, and various features were calculated for

the segments. These steps were carried out by using the eCognition software (Trimble Germany GmbH, Munich) and its multiresolution segmentation algorithm (Baatz and Schäpe, 2000). Segmentation was first carried out based on the maximum DSM. This DSM shows tree canopies in their full size, and it basically represents information corresponding to the intensity image that was created from first pulse and only pulse points. Segments were divided into high and low segments based on their mean height (difference between mean values in the maximum DSM and DTM). The threshold value was 2.5 m. Low segments were merged to each other and resegmented based on the three-channel intensity image. Fig. 3 illustrates the segmentation results. Parameters for segmentation were selected on the basis of previous experience and visual analysis of different results. They are presented in Table 2. ‘‘Scale” indirectly defines the size of the segments, and ‘‘Shape” and ‘‘Compactness” are related to the composition of the heterogeneity criterion used when merging neighbouring segments. When ‘‘Shape” is 0, the segmentation is purely based on ‘‘colour” information of the image layers, without aiming to form compact or smooth segments. In our DSM segmentation, ‘‘colour” information means height values in practice. It is important to note that suitable values for the parameters are dependent on the data, for example, range of values in use. Thus, some changes are likely to be needed for other datasets, unless the properties of these datasets are similar to ours. The features calculated for the segments

303

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Fig. 3. Segmentation levels based on the DSM (left) and intensity image (right). In the intensity level, low segments from the DSM level have been further segmented on the basis of the intensity image. High segments are the same in both cases. The width of the area is about 170 m.

Table 2 Segmentation parameters used for land cover classification. High and low segments refer to segments with their mean heights >2.5 m and 2.5 m above the ground level, respectively. Segmentation level

Image layers used

Scale

Shape

Compactness

DSM, high segments Intensity image, low segments

Maximum DSM Intensity image Ch1, Ch2 and Ch3 (each with equal weight)

15 2

0 0.01

– 0.5

and used in further analyses and classification are presented in Table 3.

Table 3 Features used in land cover classification. Feature number

3.1.2. Random forests and histogram analyses Our main interest in this study was to investigate the usefulness of multispectral ALS features in land cover classification. The reference points were used to determine training segments for different land cover classes. Training segments for buildings and trees were picked from the high segments, and training segments for asphalt, gravel, rocky areas and low vegetation were picked from the low segments. If there was a reference point inside a segment, this segment became a training segment of the corresponding class. The potential of the features for separating the classes was then investigated by using histograms and machine learning analyses with the random forests method. Functions of the Matlab software (The Mathworks, Inc., Natick, MA, USA) and its Statistics and Machine Learning Toolbox were used for this purpose. The random forests method (Breiman, 2001; Breiman and Cutler, 2004) is a further development of classification trees (Breiman et al., 1984), and it is an ensemble classification method. The basic classification tree algorithm creates one classification tree automatically based on training data. In the random forests method, a large number of classification trees is created and the classification result is obtained as a voting result. Each classification tree is constructed by using a part of the training cases. The remaining cases are called out-of-bag data. If a training set has N cases, N cases are selected at random, but with replacement, for creating a classification tree. At each node of a classification tree, m variables (features) out of M (the total number of variables) are selected at random and used to find the best split. The classification trees are not pruned. This approach is known to lead to an accurate classifier, and it also provides the possibility to calculate feature importance and an estimate of the classification error

Feature name

Intensity features 1 Brightness (mean value of the mean intensity values in different channels) 2, 3, 4 Mean intensity in Ch1, Ch2, Ch3 5, 6, 7 Intensity quantile 10% in Ch1, Ch2, Ch3 8, 9, 10 Intensity quantile 25% in Ch1, Ch2, Ch3 11, 12, 13 Intensity quantile 50% in Ch1, Ch2, Ch3 14, 15, 16 Intensity quantile 75% in Ch1, Ch2, Ch3 17, 18, 19 Intensity quantile 90% in Ch1, Ch2, Ch3 20, 21, 22 Difference between intensity quantiles 90% and 10% in Ch1, Ch2, Ch3 23, 24, 25 Intensity ratio in Ch1, Ch2, Ch3 (ratio is calculated by dividing the mean intensity in one channel by the sum of the mean intensity values in all channels) 26 Intensity ratio Ch1/Ch2 (calculated from the mean values) 27 Intensity ratio Ch1/Ch3 (calculated from the mean values) 28 Intensity ratio Ch2/Ch3 (calculated from the mean values) 29 PseudoNDVI (normalized difference vegetation index) = (Mean Ch2 – Mean Ch 3)/(Mean Ch2 + Mean Ch3) (Wichmann et al., 2015) 30 PseudoNDBI (normalized difference built-up index) = (Mean Ch1  Mean Ch2)/(Mean Ch1 + Mean Ch2) (index for Titan data modified from Morsy et al. (2016); NDBI originally presented by Zha et al. (2003)) 31, 32, 33 Standard deviation of intensity in Ch1, Ch2, Ch3 34, 35, 36 Grey-level co-occurrence matrix (GLCM) homogeneity of intensity in Ch1, Ch2, Ch3 (texture feature originally presented by Haralick et al. (1973)) DSM features 37 38 39 40 41

Maximum DSM – minimum DSM (difference between mean values) Standard deviation of the maximum DSM Standard deviation of the minimum DSM GLCM homogeneity of the maximum DSM GLCM homogeneity of the minimum DSM

304

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

based on the out-of-bag data. The method does not overfit. (Breiman and Cutler, 2004.) Both the classification tree and the random forests methods are well suited for classifying data with a large number of features. They do not require separate feature selection or normalisation of feature values. Both methods have also been successfully used for classifying ALS data (for application of random forests, see, for example, Guo et al., 2011). In our study, we used the Matlab function ‘fitensemble’ with default parameters for bagging and classification trees to construct ensembles of classification trees. The number of classification trees was 1000 in each tested classification case (these will be listed in Section 4.1). The functions ‘oobLoss’ and ‘predictorImportance’ were used to calculate the out-of-bag classification error and feature importance, respectively. The classification error was calculated as the fraction of misclassified data. The feature importance analysis is based on changes in the risk due to splits on every feature (Mathworks, 2016). 3.1.3. Classification workflow and accuracy estimation Finally, actual land cover classification for the whole study area was carried out. This included segmentation as described above (eCognition) and random forests classification separately for high and low objects (Matlab). Classification ensembles created based on the training data were used. High objects were classified as buildings and trees using intensity and DSM features. Low objects were classified as asphalt, gravel, rocky areas and low vegetation using intensity features only. As a postprocessing step, very small building objects (< 20 m2) were eliminated (Matlab). In practice, these were reclassified as trees, because it is most probably the correct class for high objects with a small area, and this significantly improves the visual appearance of the results considering buildings. The accuracy of the final results was evaluated by using the validation points. For comparison purposes, we also carried out the land cover classification by using intensity information from one channel of the Titan data only (for details, see Section 4.2). 3.2. Change detection methods 3.2.1. Change detection of buildings In the case of buildings, we tested a ‘‘second-generation” ALS-based map updating approach where two ALS datasets are available: an old one based on conventional single-channel ALS and a new one from multispectral ALS. This presents a realistic scenario for many mapping organisations in future years. The idea is that the change detection approach exploits the possibility to compare the old and new datasets, but it also makes use of the multispectral information and it is directly linked to the existing building vectors. The change detection process uses old building vectors, old minimum DSM, new minimum DSM, new maximum DSM, new multispectral intensity image and DTM as input data. It was mainly implemented with eCognition, but the classification of new buildings was carried out with Matlab, with the same classification ensemble that was used for classifying high objects in the land cover classification. The change detection procedure includes the following main steps: 1. Analysis of building segments derived from the old map. This is carried out by creating segments directly from the building vectors (one segment corresponds to one building) and classifying these using the following rules:  IF h > 2.5 m & abs_dsm_dif  3 m THEN Old building unchanged,  IF h > 2.5 m & abs_dsm_dif > 3 m THEN Height change in building,  IF h  2.5 m & dsm_dif < -2 m THEN Demolished building,  ELSE Old building unchanged,

where h is the height of the segment, i.e., the difference between the mean heights in the new minimum DSM and DTM, abs_dsm_dif is the absolute value of the difference between the mean heights in the new and old minimum DSMs, and dsm_dif is the same difference with sign. The threshold value 2.5 m is the same that was used to classify high objects, i.e. building and trees, in the land cover classification. The other threshold values were defined after experiments. The main idea with the threshold 3 m is that a real height change in a building, typically an increase, is likely to correspond to the height of at least one floor. In addition, some tolerance was needed because of the different height systems and because the shape of the buildings was slightly different in the old and new DSMs with different pixel sizes. For demolished buildings, a threshold value of -2 m was selected to also detect demolitions of the lowest buildings. The purpose of this threshold value is also to ensure that very low objects presented on the building map are not classified as demolished simply based on the 2.5 m threshold value. 2. Detection of new buildings from the new data. This is carried out by segmenting the minimum DSM outside the old buildings, separating high objects from ground with a threshold value of 2.5 m and using the random forest classification ensemble to detect buildings. Finally, new buildings smaller than 20 m2 were removed. In change detection we used the minimum DSM as a basis for segmentation because it is advantageous in the analysis of buildings. The scale parameter in segmentation was 20, and the shape parameter was 0. In the final results, new building parts (>20 m2) connected to old buildings are also visible. It is common that the shape of a building in the DSM is slightly different from the same building on the map, which results in small building parts classified as new. Further processing, such as analysing the relative size of the parts and defining a single class for each building would also be possible (Matikainen et al., 2016), but was not included in this test. It depends on further usage of the results whether such postprocessing is desirable and which are the exact criteria that should be used. Here, our aim is to demonstrate automated change detection with the new data, but not to present universally applicable classification rules. The application of the method in another area or larger area may require modifications and additional rules. 3.2.2. Change detection of roads For change detection the road centre line vectors were extracted from the obsolete version of the Topographic Database (corresponding to year 2000). A test point set was derived from the vectors with circa 2 m interval. Each point was compared to the land cover classification result from the multispectral ALS data using geographic information system (GIS) software QGIS (QGIS, 2016) and classified using the following rules:  IF class = (asphalt OR gravel) THEN label point as Old road unchanged.  ELSE label point as Old road not found. Potential new roads were detected by extracting asphalt and gravel classes from the classification result. 4. Results 4.1. Analysis of class separability and feature importance The importance of various features (see Table 3) for separating different classes was tested by analysing the six classification cases shown in Table 4. In each case, 1000 classification trees were created on the basis of the training data, and the out-of-bag

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

305

Table 4 Results of the random forests analysis. In each case, 1000 classification trees were created based on training data. Cases marked with asterisk (*) are visualised in Figs. 4–7. Classification case: Classes to separate

Features in use

Out-of-bag classification error

Five most important features (feature numbers from Table 3 in brackets)

1: Buildings/Trees 2: Buildings/Trees

Intensity DSM

0.006 0.018

3: Buildings/Trees *)

Intensity and DSM Intensity

0.000

GLCM hom. Ch2 (35), Ratio Ch2 (24), PseudoNDVI (29), Ch2/Ch3 (28), PseudoNDBI (30) MaxDSM–MinDSM (37), Std MaxDSM (38), Std MinDSM (39), GLCM hom. MaxDSM (40), GLCM hom. MinDSM (41) Ratio Ch2 (24), GLCM hom. Ch2 (35), PseudoNDVI (29), Ch2/Ch3 (28), Ch1/Ch2 (26)

0.031

Ch1/Ch2 (26), PseudoNDBI (30), Ratio Ch3 (25), Ch1/Ch3 (27), Ratio Ch2 (24)

Intensity

0.019

PseudoNDVI (29), Ch2/Ch3 (28), PseudoNDBI (30), Ch1/Ch2 (26), Ratio Ch2 (24)

Intensity

0.026

Ch1/Ch3 (27), Ratio Ch3 (25), Ratio Ch1 (23), Quantile 75% Ch1 (14), Quantile 90% Ch1 (17)

4: Asphalt/Gravel/Rocky areas/Low vegetation *) 5: Road-like surfaces (included asphalt and gravel)/Rocky areas/Low vegetation *) 6: Asphalt/Gravel *)

classification error and feature importance were calculated. The results are shown in Table 4 and Figs. 4–7. Our main interest was to investigate the intensity features and ground-level classes that are challenging to distinguish based on single-channel ALS data. For buildings and trees, different feature combinations (intensity or DSM) were tested. In this case, it is well known that DSM data provide useful features, such as texture (e.g., Hug, 1997). Here, we were interested in studying what is the additional contribution of the multispectral intensity features. Histograms of five most important intensity features are shown for three classification cases (3, 5, 6) in Fig. 8. 4.2. Land cover classification results To obtain the final land cover classification results for the entire study area, high objects were classified into buildings and trees using the classification ensemble from Classification case 3 (see Table 4). Low objects were classified into asphalt, gravel, rocky areas and low vegetation using the classification ensemble from Classification case 4. The final results after postprocessing, i.e. elimination of the smallest buildings, are shown in Fig. 9. Fig. 10 shows three

close-up views of the results. The quality of the results compared with the validation points is presented in Table 5. Two of the validation points were located under the water mask, and were thus excluded from the analysis. For comparison, we also carried out the land cover classification by using intensity information from one channel of the Titan data only. Conventional single-channel ALS systems normally work in the infrared or near-infrared region (cf., Titan channels 1 and 2). We selected channel 2 for this additional test because it appeared more useful for land cover classification than channel 1 on the basis of the feature importance analyses (Section 4.1). However, to keep the effect of point density constant, the DSMs and DTM used in this test were still based on points from all three Titan channels. The segmentation parameters were the same as presented in Table 2, but only Ch2 was used in the intensity image segmentation stage. The features used in classification included those features from Table 3 that can be calculated from Ch2 only or the DSMs. In practice, 9 intensity features and 5 DSM features were available. Brightness, ratio features, PseudoNDVI, PseudoNDBI, and all intensity features based on Ch1 or Ch3 were unavailable. Similar to the multispectral classification, high objects

Fig. 4. Importance of features in separating buildings and trees from each other. Both intensity and DSM features were available (Classification case 3).

306

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Fig. 5. Importance of features in separating asphalt, gravel, rocky areas and low vegetation from each other. Intensity features were available (Classification case 4).

Fig. 6. Importance of features in separating road-like surfaces, rocky areas and low vegetation from each other. Road-like surfaces included asphalt and gravel. Intensity features were available (Classification case 5).

were classified using the intensity and DSM features and low objects were classified using the intensity features only. The quality of the classification results was evaluated by comparing them with the validation points, and the results are presented in Table 6.

4.3. Change detection results Fig. 11 shows the change detection results. Building changes and road changes have been combined in the figure and overlaid on the Titan intensity image to illustrate how automatically

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

307

Fig. 7. Importance of features in separating asphalt and gravel from each other. Intensity features were available (Classification case 6).

derived changes might be presented for an operator who is updating a map database. However, it should be noted that in practice our old building and road maps were not from the same date. The results for buildings are shown for a subarea in more detail in Fig. 12. The quality of the change detection results for buildings was evaluated visually by comparing them with the old and new building maps, the old and new DSMs, and the Titan intensity image. The quality was generally satisfactory:  Except for a few multi-storey car parks mentioned in the next list and one low-rise residential building, all new major buildings were detected, at least partly.  In some cases the results revealed very new changes that were not yet presented in the reference map from 2015.  A few buildings were already presented on the old map, but they had apparently been under construction when the old DSM data were acquired. There was a clear height change between the two DSMs, and these buildings were now classified as changed. For example, one such residential building can be seen in Fig. 12 (highlighted with a black circle). The main problems detected in the results were the following:  A few multi-storey car parks were not detected as buildings, although they were presented as buildings on the map. The problem here is that the car parks are located on a hill slope and their upper levels are partly on the ground level or near it. They were thus classified as ground.  One high-rise building was falsely classified as changed. The shape of the buildings was rougher in the new DSM with 1 m pixel size than in the old DSM with 0.3 m size. In the case of the misclassified building, the difference in shape led to mean height difference larger than the threshold value.  Tree cover caused some misclassifications of small buildings, but also uncertainty in evaluating the results for these buildings.

The road change detection results are best suited for visual inspection, since roads are combined in the same class with other gravel and asphalt surfaces. The result shows plenty of new roads that have been constructed in the area over the last 16 years (see Fig. 11). Driveways not included in the road database, parking lots and construction sites are also visible. The result also clearly shows that several roads have been demolished or moved to a new location. Tree cover was the main reason for not detecting roads that have not changed, therefore leaf-off data would probably give more accurate results. Road classes in the study area varied from express ways to narrow driveways and cycle ways. Most of the problems in road detection were caused by the trees occluding narrow roads. Minor problems were caused by road markings with high reflectance and reflections from vehicles. These could be avoided to some extent by applying simple postprocessing methods (e.g. finding vehicle height objects on roads).

5. Discussion Visual analysis of multispectral intensity images produced from the Optech Titan point clouds suggests that multispectral ALS data are very promising for mapping applications, especially considering object-based automated analyses. Unlike conventional passive aerial imaging, the multispectral ALS technology is independent of external illumination conditions, and there are no shadows on intensity images produced from the data (see Fig. 2). This means a significant advantage when the development of automated classification and change detection procedures is considered. For example, roads near trees or buildings are better visible on the Titan intensity image from leaf-on conditions than on an aerial ortho image from leaf-off conditions. The direct availability of both geometric and multispectral information from a single sensor also simplifies preprocessing and data analysis workflows. Results from the random forests analysis clearly support the hypothesis that the multispectral intensity information is useful for land cover mapping. Out-of-bag classification errors calculated

308

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Fig. 8. Histograms of training segments in different classes. Left column: buildings and trees (Classification case 3). Middle column: road-like surfaces (includes asphalt and gravel), rocky areas and low vegetation (Classification case 5). Right column: asphalt and gravel (Classification case 6). For each case, the five most important features according to the feature importance analysis (Table 4, Figs. 4–7) are shown.

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

309

Fig. 9. Land cover classification results for the entire study area. The black rectangles show the subareas of Fig. 10. The boundary of the study area is mainly based on the water mask, which contains data from the National Land Survey of Finland Topographic Database 2015.

based on the training data were remarkably small in all tested cases (0.00–0.03; Table 4). According to the feature importance analysis, features derived using multiple channels, i.e., intensity ratios (feature numbers 23–28), PseudoNDVI (29) and PseudoNDBI (30), were very useful. Some of these always appeared among the most important features. Histograms of the most important features (Fig. 8) further illustrate their high potential for separating classes. Comparison of the classification results to separate validation points gave an overall accuracy of 96%, which is a good result in classification of a large study area with six classes. The additional classification test with intensity information from Ch2 alone gave a slightly lower overall accuracy of about 93%. Considering individual classes, the results suggest that multispectral ALS data have high potential for separating ground surface classes from each other. Promising results for ground surface classes have also been reported in other studies with Titan data (Wichmann et al., 2015; Bakuła et al., 2016; Zou et al., 2016). In our case, the low classes included asphalt, gravel, rocky areas and low vegetation. The good results for these classes are very important from the practical point of view because they cannot

be separated on the basis of height data alone, and results based on one channel intensity data are not very reliable. As shown in Tables 5 and 6, the difference in quality between the multispectral and single-channel classifications was more pronounced when the ground surface classes are considered. In particular, the completeness and correctness values of classes gravel and rocky area were clearly lower when intensity from Ch2 alone was used. The correctness of low vegetation also decreased, meaning that more points from other classes were misclassified as low vegetation. The benefits of the multispectral information can be seen in the histograms of Fig. 8. Asphalt and gravel as one class (roads, etc.) were well separated from rocky areas and low vegetation using the features PseudoNDVI (29) and Ch2/Ch3 (28). The histograms of rocky areas and low vegetation, however, had much overlap in the case of these features. In the case of PseudoNDBI (30) and Ch1/Ch2 (26), rocky areas appeared more like roads and were well separated from low vegetation. As already discussed in Section 4.3, minor problems in road classification were caused by road markings and vehicles. Examples of such classification errors can be seen in Fig. 10 (left column).

310

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

Fig. 10. Land cover classification results for three subareas (330 m  330 m each). The upper row shows the Titan intensity image for the same areas. For the legend of the classification results, see Fig. 9.

Table 5 Confusion matrix and accuracy estimates for the land cover classification results. Multispectral Titan data were used in classification. Classification result

Reference data (validation points) Building

Building 125 Tree 3 Asphalt 1 Gravel 1 Rocky area 0 Low veg. 0 Total 130 Completeness 96.2% Overall accuracy: 95.9%, Kappa coefficient: 0.95

Tree

Asphalt

Gravel

Rocky area

Low veg.

Total

Correctness

0 142 0 0 0 0 142 100%

0 0 149 5 1 2 157 94.9%

0 0 3 12 0 0 15 80.0%

0 0 0 0 40 4 44 90.9%

0 0 0 0 3 73 76 96.1%

125 145 153 18 44 79 564

100% 97.9% 97.4% 66.7% 90.9% 92.4%

Table 6 Confusion matrix and accuracy estimates for the additional land cover classification test, where intensity information from Titan channel 2 only was used in classification. Classification result

Reference data (validation points) Building

Building 124 Tree 4 Asphalt 0 Gravel 1 Rocky area 0 Low veg. 1 Total 130 Completeness 95.4% Overall accuracy: 92.9%, Kappa coefficient: 0.91

Tree

Asphalt

Gravel

Rocky area

Low veg.

Total

Correctness

0 142 0 0 0 0 142 100%

0 0 147 2 6 2 157 93.6%

0 0 3 7 4 1 15 46.7%

0 0 1 4 30 9 44 68.2%

0 0 0 1 1 74 76 97.4%

124 146 151 15 41 87 564

100% 97.3% 97.4% 46.7% 73.2% 85.1%

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

311

Fig. 11. Change detection results for buildings and roads overlaid on the Titan intensity image. The black rectangle shows the subarea of Fig. 12. Original map data for old buildings Ó the City of Espoo. Original map data for old roads is based on the National Land Survey of Finland Topographic Database 2000.

Fig. 12. Change detection of buildings in a subarea. Left: old DSM; Middle: new DSM and old building vectors; Right: change detection results (for legend, see Fig. 11). The building inside the black circle is an example of a building with a clear height change between the two DSMs (mean height difference about 6 m). Original map data for old buildings Ó the City of Espoo.

In the comparison to validation points, the class with the lowest completeness and correctness values in the multispectral classification was gravel. Confusion occurred mainly with asphalt. Since the number of gravel points was relatively small and gravel is a

class with significant in-class variation, it is possible that the training and validation points for gravel had somewhat different characteristics. The potential of the Titan data for detecting roads and classifying road surface into gravel and paved was further studied

312

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313

and discussed in Karila et al. (2017). In that study, comparison to a large reference dataset based on map data showed accuracy of 90% and 95% in classifying detected roads into gravel and paved, respectively. Buildings and trees can be separated from each other by using DSM data and geometric features. However, also in this case the multispectral feature Ratio Ch2 (24) was the most important among the tested intensity and DSM features. Another useful feature with almost equal importance was the textural feature GLCM homogeneity Ch2 (35). Some other multispectral features also appeared more important than the DSM features. Buildings and trees had very high completeness and correctness values in the classification (96–100%). The numerical accuracy values remained almost at the same level when intensity from Ch2 alone was used. The variation of roof materials and colours was large in our study area, and some errors occurred in the classifications, such as misclassification of a bright roof partly as tree (Fig. 10, middle column). Classification of bridges as buildings, which occurred in our preliminary study (Matikainen et al., 2016), was now mostly avoided as bridges became part of the ground surface in our new DTM created automatically from the Titan data. The results for buildings and trees suggest that multispectral ALS data further improve the possibility to develop automated mapping methods for these objects with an accuracy level acceptable for operational work. In the future, more land cover classes, such as different types of vegetation, could be included in the analysis. Our Titan data acquired from a flying altitude of 650 m also show some potential for detecting smaller objects than tested in this study. For example, many poles, traffic signs and buses were now classified as buildings and trees. Some examples of misclassified traffic signs and light poles above the road can be seen in Fig. 10 (left column). Especially for classifying small objects and objects and surfaces under trees, classification of the data in point-based format would be useful. There is also further potential in analysing the intensity and geometric properties of the different channels in more detail. From the practical point of view, an interesting topic is the effect of data acquisition date on the classification accuracy. Unlike the present study, data used operationally for mapping is normally acquired early in spring in leaf-off conditions. This has clear advantages for the visibility of many objects, but the spectral properties might be less useful for automated classification. If the intention is for the data to be used operationally for national mapping, there are also other questions such as the flying altitude, which should be optimised for effective work. Change detection of buildings and roads was demonstrated in a new residential area with plenty of changes. The change detection results clearly illustrate the development of the area. This type of results could be helpful in operational map updating work to show where changes are needed in the database. The final decision on changes and outlining of objects could be done by an operator. Further automated steps such as boundary extraction and 3D modelling of the new buildings are also possible. Considering the change detection methods, the availability of two DSMs allowed straightforward height comparisons for existing buildings. Multispectral information was exploited in the detection of new buildings. The method detected major changes such as new residential buildings well. For difficult cases such as checking small buildings in the forest, more sophisticated methods should be developed. In the case of roads, the special capability of the multispectral ALS data for ground surface classification was exploited. The change detection was simply based on the comparison of the asphalt and gravel surfaces from the land cover classification with old road vectors. Using more advanced road network extraction algorithms, the change detection could be more automatic in the future. The development of automated change detection and map updating

procedures could also be extended to other types of objects including, for example, forests, agricultural fields and land cover of the ground surface. Based on our study, we believe that multispectral ALS data provide a good basis for this work. The first methods can utilise multispectral data together with information from existing map databases and single-channel ALS datasets. Later, further improvements are expected when comparison approaches based on two or more multispectral ALS datasets can be developed. 6. Conclusions With the Optech Titan sensor, multispectral information is for the first time available for 3D ALS point clouds from a single sensor. Unlike passive aerial imaging, the technology is independent of external illumination conditions, and there are no shadows on intensity images produced from the data. These are significant advantages when the development of automated analysis methods is considered. We studied the potential of this new single-sensor technology in map updating, especially in object-based automated land cover classification and change detection in a suburban area. The study is the first one on this topic in international journals. Results from a random forests analysis suggest that the multispectral intensity information is very useful for land cover classification, especially when considering low objects such as roads and land cover classes of the ground surface. It also further enhances the possibility to automate building and tree detection with high reliability. The overall accuracy of the land cover classification results with six classes was 96% compared with validation points. The main improvements compared to classification of singlechannel data were achieved for ground-level classes. According to a feature importance analysis, multispectral features based on several channels were generally more useful than those based on one channel. Automatic change detection for buildings and roads was also demonstrated by utilising the new multispectral ALS data in combination with old map vectors. In change detection of buildings, an old DSM based on single-channel ALS data was also used. In the future, the development of automated change detection and map updating procedures can be extended to other types of objects and land cover classes. In addition to methodological developments, important topics for further study include the optimal data acquisition time for automated mapping applications and the use of multitemporal multispectral data. Acknowledgements The authors wish to thank Juha Kareinen for providing information on the ground classification process used in the production work of the National Land Survey (NLS), TerraTec Oy for cooperation in Optech Titan data acquisition, the City of Espoo for map data, and the NLS Department of Production for aerial ortho images and map data. The financial support from the Academy of Finland (projects ‘Integration of Large Multisource Point Cloud and Image Datasets for Adaptive Map Updating’ (project decision number 295047) and ‘Centre of Excellence in Laser Scanning Research’ (272195)) is gratefully acknowledged. References Ahokas, E., Kaasalainen, S., Hyyppä, J., Suomalainen, J., 2006. Calibration of the Optech ALTM 3100 laser scanner intensity data using brightness targets. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXVI (Part 1). 7 p. Ahokas, E., Hyyppä, J., Yu, X., Liang, X., Matikainen, L., Karila, K., Litkey, P., Kukko, A., Jaakkola, A., Kaartinen, H., Holopainen, M., Vastaranta, M., 2016. Towards automatic single-sensor mapping by multispectral airborne laser scanning. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XLI-B3, 155–162. Axelsson, P., 2000. DEM generation from laser scanner data using adaptive TIN models. Int. Archiv. Photogram. Rem. Sens. XXXIII (Part B4), 110–117.

L. Matikainen et al. / ISPRS Journal of Photogrammetry and Remote Sensing 128 (2017) 298–313 Baatz, M., Schäpe, A., 2000. Multiresolution segmentation – an optimization approach for high quality multi-scale image segmentation. In: Strobl, J., Blaschke, T., Griesebner, G. (Eds.), Angewandte Geographische Informationsverarbeitung XII: Beiträge Zum AGIT-Symposium Salzburg 2000. Wichmann, Heidelberg, pp. 12–23. Bakuła, K., 2015. Multispectral airborne laser scanning – a new trend in the development of LiDAR technology. Archiwum Fotogrametrii, Kartografii i Teledetekcji 27, 25–44. Bakuła, K., Kupidura, P., Jełowicki, Ł., 2016. Testing of land cover classification from multispectral airborne laser scanning data. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XLI-B7, 161–169. Breiman, L., 2001. Random forests. Mach. Learn. 45 (1), 5–32. Breiman, L., Cutler, A., 2004. Random Forests. (accessed 15 November, 2016). Breiman, L., Friedman, J.H., Olshen, R.A., Stone, C.J., 1984. Classification and Regression Trees. Wadsworth Inc, Belmont, CA, USA. Brenner, C., 2010. Building extraction. In: Vosselman, G., Maas, H.-G. (Eds.), Airborne and Terrestrial Laser Scanning. Whittles Publishing, Dunbeath, Scotland, UK, pp. 169–212. Briese, C., Pfennigbauer, M., Ullrich, A., Doneus, M., 2013. Multi-wavelength airborne laser scanning for archaeological prospection. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XL-5/W2, 119–124. Champion, N., Rottensteiner, F., Matikainen, L., Liang, X., Hyyppä, J., Olsen, B.P., 2009. A test of automatic building change detection approaches. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXVIII (Part 3/W4), 145–150. Chen, Y., Räikkönen, E., Kaasalainen, S., Suomalainen, J., Hakala, T., Hyyppä, J., Chen, R., 2010. Two-channel hyperspectral LiDAR with a supercontinuum laser source. Sensors 10, 7057–7066. Clode, S., Rottensteiner, F., Kootsookos, P., Zelniker, E., 2007. Detection and vectorization of roads from lidar data. Photogramm. Eng. Remote Sens. 73 (5), 517–535. Eitel, J.U.H., Höfle, B., Vierling, L.A., Abellán, A., Asner, G.P., Deems, J.S., Glennie, C.L., Joerg, P.C., LeWinter, A.L., Magney, T.S., Mandlburger, G., Morton, D.C., Müller, J., Vierling, K.T., 2016. Beyond 3-D: The new spectrum of lidar applications for earth and ecological sciences. Remote Sens. Environ. 186, 372–392. Fernandez-Diaz, J.C., Carter, W.E., Glennie, C., Shrestha, R.L., Pan, Z., Ekhtari, N., Singhania, A., Hauser, D., Sartori, M., 2016. Capability assessment and performance metrics for the Titan multispectral mapping lidar. Rem. Sens. 8 (11), 936. Guo, L., Chehata, N., Mallet, C., Boukir, S., 2011. Relevance of airborne lidar and multispectral image data for urban scene classification using Random Forests. ISPRS J. Photogram. Rem. Sens. 66, 56–66. Hakala, T., Suomalainen, J., Kaasalainen, S., 2012. Full waveform active hyperspectral lidar. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXIX-B7, 459–462. Hakala, T., Nevalainen, O., Kaasalainen, S., Mäkipää, R., 2015. Technical note: multispectral lidar time series of pine canopy chlorophyll content. Biogeosciences 12 (5), 1629–1634. Haralick, R.M., Shanmugam, K., Dinstein, I., 1973. Textural features for image classification. IEEE Trans. Syst. Man Cybernet. 3 (6), 610–621. Höfle, B., Pfeifer, N., 2007. Correction of laser scanning intensity data: data and model-driven approaches. ISPRS J. Photogram. Rem. Sens. 62 (6), 415–433. Holland, D.A., Sanchez-Hernandez, C., Gladstone, C., 2008. Detecting changes to topographic features using high resolution imagery. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXVII (Part B4), 1153–1158. Hopkinson, C., Chasmer, L., Gynan, C., Mahoney, C., Sitar, M., 2016. Multisensor and multispectral LiDAR characterization and classification of a forest environment. Can. J. Remote. Sens. 42 (5), 501–520. Hu, X., Tao, C.V., Hu, Y., 2004. Automatic road extraction from dense urban area by integrated processing of high resolution imagery and lidar data. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXV (Part B3), 320–324. Hug, C., 1997. Extracting artificial surface objects from airborne laser scanner data. In: Gruen, A., Baltsavias, E.P., Henricsson, O. (Eds.), Automatic Extraction of Man-Made Objects From Aerial and Space Images (II). Birkhäuser Verlag, Basel, Switzerland, pp. 203–212. Kaasalainen, S., Ahokas, E., Hyyppä, J., Suomalainen, J., 2005. Study of surface brightness from backscattered laser intensity: calibration of laser data. IEEE Geosci. Remote Sens. Lett. 2 (3), 255–259. Kaasalainen, S., Pyysalo, U., Krooks, A., Vain, A., Kukko, A., Hyyppä, J., Kaasalainen, M., 2011. Absolute radiometric calibration of ALS intensity data: effects on accuracy and target classification. Sensors 11 (11), 10586–10602. Karila, K., Matikainen, L., Puttonen, E., Hyyppä, J., 2017. Feasibility of multispectral airborne laser scanning data for road mapping. IEEE Geosci. Remote Sens. Lett. 14 (3), 294–298. Leigh, H.W., Magruder, L.A., 2016. Using dual-wavelength, full-waveform airborne lidar for surface classification and vegetation characterization. J. Appl. Rem. Sens. 10 (4), 045001.

313

Luzum, B., Starek, M., Slatton, K.C., 2004. Normalizing ALSM intensities. Geosensing Engineering and Mapping (GEM) Center Report No. Rep_2004-07-001, Civil and Coastal Engineering Department, University of Florida, Gainesville, FL, USA. Malpica, J.A., Alonso, M.C., Papí, F., Arozarena, A., Martínez De Agirre, A., 2013. Change detection of buildings from satellite imagery and lidar data. Int. J. Remote Sens. 34 (5), 1652–1675. Mathworks, 2016. Online Documentation for Statistics and Machine Learning Toolbox, Version R2016a. Matikainen, L., Karila, K., 2011. Segment-based land cover mapping of a suburban area – comparison of high-resolution remotely sensed datasets using classification trees and test field points. Rem. Sens. 3 (8), 1777–1804. Matikainen, L., Hyyppä, J., Hyyppä, H., 2003. Automatic detection of buildings from laser scanner data for map updating. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XXXIV (Part 3/W13), 218–224. Matikainen, L., Hyyppä, J., Ahokas, E., Markelin, L., Kaartinen, H., 2010. Automatic detection of buildings and changes in buildings for updating of maps. Rem. Sens. 2 (5), 1217–1248. Matikainen, L., Hyyppä, J., Litkey, P., 2016. Multispectral airborne laser scanning for automated map updating. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XLI-B3, 323–330. Miller, C.I., Thomas, J.J., Kim, A.M., Metcalf, J.P., Olsen, R.C., 2016. Application of image classification techniques to multispectral lidar point cloud data. In: Proc. SPIE 9832 (Laser Radar Tehnology and Applications XXI), 98320X, 12 p. Morsy, S., Shaker, A., El-Rabbany, A., LaRocque, P.E., 2016. Airborne multispectral lidar data for land-cover classification and land/water mapping using different spectral indexes. ISPRS Ann. Photogram. Rem. Sens. Spatial Inform. Sci. III-3, 217–224. Murakami, H., Nakagawa, K., Hasegawa, H., Shibata, T., Iwanami, E., 1999. Change detection of buildings using an airborne laser scanner. ISPRS J. Photogram. Rem. Sens. 54, 148–152. Pfennigbauer, M., Ullrich, A., 2011. Multi-wavelength airborne laser scanning. In: Proc. 2011 International Lidar Mapping Forum, ILMF, New Orleans, 10 p. Puttonen, E., Hakala, T., Nevalainen, O., Kaasalainen, S., Krooks, A., Karjalainen, M., Anttila, K., 2015. Artificial target detection with a hyperspectral LiDAR over 26h measurement. Opt. Eng. 54 (1), 013105. Puttonen, E., Briese, C., Mandlburger, G., Wieser, M., Pfennigbauer, M., Zlinszky, A., Pfeifer, N., 2016. Quantification of overnight movement of birch (betula pendula) branches and foliage with short interval terrestrial laser scanning. Front. Plant Sci. 7, 222. QGIS Development Team, 2016. QGIS Geographic Information System. Open Source Geospatial Foundation. (accessed 17 February, 2017). Richter, R., Kyprianidis, J.E., Döllner, J., 2013. Out-of-core GPU-based change detection in massive 3D point clouds. Trans. GIS 17 (5), 724–741. Rottensteiner, F., Trinder, J., Clode, S., Kubik, K., 2007. Building detection by fusion of airborne laser scanner data and multi-spectral images: Performance evaluation and sensitivity analysis. ISPRS J. Photogram. Rem. Sens. 62, 135–149. Teo, T.-A., Shih, T.-Y., 2013. Lidar-based change detection and change-type determination in urban areas. Int. J. Remote Sens. 34 (3), 968–981. Vauhkonen, J., Hakala, T., Suomalainen, J., Kaasalainen, S., Nevalainen, O., Vastaranta, M., Holopainen, M., Hyyppä, J., 2013. Classification of spruce and pine trees using active hyperspectral LiDAR. IEEE Geosci. Remote Sens. Lett. 10 (5), 1138–1141. Vosselman, G., Kessels, P., Gorte, B., 2005. The utilisation of airborne laser scanning for mapping. Int. J. Appl. Earth Obs. Geoinf. 6, 177–186. Wagner, W., 2010. Radiometric calibration of small-footprint full-waveform airborne laser scanner measurements: basic physical concepts. ISPRS J. Photogram. Rem. Sens. 65 (6), 505–513. Wagner, W., Ullrich, A., Ducic, V., Melzer, T., Studnicka, N., 2006. Gaussian decomposition and calibration of a novel small-footprint full-waveform digitising airborne laser scanner. ISPRS J. Photogram. Rem. Sens. 60 (2), 100– 112. Wang, C.-K., Tseng, Y.-H., Chu, H.-J., 2014. Airborne dual-wavelength LiDAR data for classifying land cover. Rem. Sens. 6 (1), 700–715. Wichmann, V., Bremer, M., Lindenberger, J., Rutzinger, M., Georges, C., PetriniMonteferri, F., 2015. Evaluating the potential of multispectral airborne lidar for topographic mapping and land cover classification. ISPRS Ann. Photogram. Rem. Sens. Spatial Inform. Sci. II-3/W5, 113–119. Yan, W.Y., Shaker, A., El-Ashmawy, N., 2015. Urban land cover classification using airborne LiDAR data: a review. Remote Sens. Environ. 158, 295–310. Zha, Y., Gao, J., Ni, S., 2003. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. Int. J. Remote Sens. 24 (3), 583–594. Zou, X., Zhao, G., Li, J., Yang, Y., Fang, Y., 2016. 3D land cover classification based on multispectral lidar point clouds. Int. Archiv. Photogram. Rem. Sens. Spatial Inform. Sci. XLI-B1, 741–747.