2.1.

Multisource Data Analysis Approaches in Mapping

The integration and fusion of different remotely sensed datasets have received much attention during recent decades (see, e.g., Refs. 18 and 19) although wider discussions on the possible roles of airborne, spaceborne, and ground-based remote sensing techniques for a given application are not common. In the field of mapping, typical examples of data integration in research literature include the use of ALS and aerial image data or satellite image data for mapping buildings (see, e.g., Refs. 20 and 21) and land cover (see, e.g., Ref. 22). As mentioned above, Holland et al.¹⁴ discussed the numerous survey techniques from the viewpoint of a mapping agency. Data analysis methods were not presented in their article. The integration of ALS and MLS data has been discussed, for example, by Rutzinger et al.²³ regarding the extraction of vertical walls. Zhou and Vosselman²⁴ compared curbstone mapping results from ALS and MLS data. Zhu and Hyyppä²⁵ presented an approach for the 3-D modeling of railway environments by utilizing both ALS and MLS data. Bamber et al.²⁶ presented a digital elevation model of Greenland based on various spaceborne and airborne datasets and digitized map data.

2.2.

Mapping with Multispectral ALS Data

The feasibility of conventional single-channel ALS data for creating DTMs and mapping elevated objects, such as buildings and trees, is well known (see, e.g., Ref. 27). Several studies related to change detection and map updating have also been carried out. Change detection between multitemporal single-channel ALS datasets has been discussed, for example, by Murakami et al.,²⁸ Richter et al.,²⁹ and recently by Tran et al.³⁰ Since the emergence of the first operational multispectral ALS system, Optech Titan, in late 2014, the potential of this single-sensor technology for several applications has also been demonstrated. These include land cover classification (see, e.g., Refs. 89.10.–11, 31, and 32), tree species classification,^33,34 road mapping,³⁵ and map updating.¹⁰ In particular, the ability of the technique to map ground-level classes—such as low vegetation, roads, paved and unpaved surfaces—is a valuable addition to the conventional ALS technique.

The Optech Titan sensor has three channels: infrared 1550 nm (channel 1), near-infrared 1064 nm (channel 2), and green 532 nm (channel 3). Each channel produces a separate point cloud. The most common approach to utilizing the multispectral intensity information has been to produce multispectral intensity images in raster format. A point-based approach that exploits the nearest neighbor approach in order to merge intensity information from different channels has also been used.⁹ Radiometric correction of Optech Titan data was discussed by Yan and Shaker.³⁶ To our knowledge, change detection based on multitemporal multispectral ALS data has not yet been investigated.

2.3.

Change Detection from Optical Satellite Images

Monitoring Earth surface dynamics and tracking changes from satellite images has been an important research topic in the past, and it still is. For example, multispectral optical data from Landsat satellite series have been applied for land use mapping for decades. The opening of the Landsat archive for free use has further increased the use of Landsat data.³⁷ A review of change detection using Landsat has been presented by Zhu.³⁸ Generally, a number of change detection techniques have been developed, and it continues to increase.^39–41 When selecting a change detection method, the characteristics of the remotely sensed images used and the phenomenon under study should be considered.⁴² In a built-up area, this could include, for example, the type and density of urban structures and vegetation, and seasonal changes in the landscape.

The Copernicus Sentinel-2 mission⁴³ is a constellation of twin polar-orbiting satellites with global coverage and a revisit time of five days at the equator with the two satellites. The satellites Sentinel-2A and Sentinel-2B were launched in June 2015 and March 2017, respectively. Each of the satellites carries a multispectral optical instrument with a swath width of 290 km. They acquire information in 13 spectral bands (from 443 to 2190 nm) at different resolutions (10, 20, and 60 m). The service is mainly focused on land monitoring and security. All the data are available for scientific and commercial use, free of charge. One of the first insights into the use of Sentinel-2 data for urban areas was provided by Lefebvre et al.⁴⁴ In their study, change detection approaches were divided into image-to-image comparison and postclassification comparison. They used the fusion of the independent image classifications of Sentinel-2 and Landsat 8 to update the Copernicus High Resolution Imperviousness Layer. Pesaresi et al.⁴⁵ tested the applicability of a classification method called symbolic machine learning on Sentinel-2 imagery and also demonstrated the added value of Sentinel-2 products over Landsat-derived products for mapping built-up areas. Being a recently launched system, the real usability of Sentinel-2 data has not yet been widely and fully analyzed. There are also many commercial satellite systems providing data with higher spatial resolution, but they exceed the scope of this study.

2.4.

Mapping of Urban Objects from MLS Data

Much research has been carried out on the mapping of urban objects from MLS data including pole-like objects (e.g., Refs. 4647.48.–49) and more general object classes (see, e.g., Refs. 50 and 51). More recently, Li et al.⁵² segmented pole-like road furniture, i.e., street lights and traffic signs, into different components based on their spatial relations and compared a proposed knowledge-driven method and different machine learning classifiers to label these components. Also, the extraction of building facades has been studied^23,51 as well as the retrieval of road markings and curbstones.⁵³ The MLS classification studies reviewed above follow the object-based approach, that is, a point cloud is first segmented into objects and then the objects are classified. Also, two other approaches exist, namely the point-based approach⁵⁴ and the voxel- or supervoxel-based approach (see, e.g., Ref. 55) where either single points or (super) voxels are classified. In this study, we adopted the object-based approach because it can preserve the spatial relations between different objects.

The most common object recognition workflow starts with the removal of ground and building facade points.^23,48,50,51 Then, a connected component analysis^48,51 or some more advanced segmentation method, such as graph-cut,⁵⁰ is carried out on the above-ground points. Finally, the above-ground objects are classified into different classes, such as cars, lamp posts, and traffic signs. The classification algorithms may apply machine-learning methods^50,51 or knowledge-based rules.⁴⁸ Some approaches only concentrate on finding objects of a pole-like shape.^46,47,49

3.1.

Study Area, Reference Data, and Map Data

Our study area is located in Espoonlahti, in the city of Espoo, about 15 to 20 km from the center of Helsinki. The area is a growing suburban area with high-rise and low-rise residential areas and an industrial area. Areas that are not built-up are mainly covered by coniferous and mixed forest. In method development and testing, an “intensive study area” covered by both ALS and Sentinel-2 datasets was used. This area is illustrated in Fig. 1 (the center point is approx. 60° 9′ 18” N, 24° 38’ 24” E). The Sentinel-2 satellite images used in the study covered a larger area. On the outskirts of this larger area, there are also some agricultural fields and a golf course. These areas, however, were not used in method development or evaluation. The study area is located on the coast, and the multispectral ALS data and Sentinel-2 images covered some water areas, but these were also excluded from the analyses as our study concentrated on the changing suburban area.

Fig. 1

An illustration of the intensive study area covered by different datasets (the yellow polygon). The Sentinel-2 satellite image in the background is from June 29, 2016 (Copernicus Sentinel data [2016]). A water mask that contained data from the NLS Topographic Database 2015 was used to define the boundary of the study area.

A permanent land cover classification test field consisting of reference points with known land cover is located in the study area (see Refs. 10 and 56). In the present study, 427 “unchanged” reference points were used as training points in Sentinel-2 satellite image analyses to define training segments for developing a change detection method (with the classes “changed” and “unchanged”). The points used in the study were selected so that they cover part of the intensive study area, and the density of the points is feasible for segment-based analysis of Sentinel-2 satellite images ($1\text{point}/100\mathrm{m}\times 100\mathrm{m}$; see Fig. 1). Part of the intensive study area remained outside the training points for quality evaluation purposes (quality evaluation in this study was based on ALS classifications rather than the reference points). In addition, only those points were used that remained unchanged during the time period of interest, that is, between August 2015 and June 2016. To select unchanged points, two versions of the reference points were used. The first one is the same as that used in Ref. 10 and corresponds to the situation in August 2015. The second is a newer one corresponding to the situation in June 2016.⁵⁷ The main data sources used in the manual checking and updating of the land cover information of these point sets were multispectral ALS data acquisitions from the corresponding points of time. To select unchanged points for the satellite image analysis, the land cover was generalized into three classes: “built-up area,” “tree,” and “low vegetation.” Unchanged points considering these classes were selected.

For evaluating the quality of change detection results from multitemporal multispectral ALS data, a specific set of test points was collected. These were collected from changes that are important for topographic mapping and included 34 points on new buildings, 56 points on new roads, and 57 points on new clear-cut forest areas. The selection of these points was based on visual comparison of the two multispectral ALS datasets used in the study. Map data were used as an aid. Some of the new buildings and roads were still under construction, but they were already clearly visible in the ALS data from the second date. The clear-cut points were selected from areas where trees had been cut, but the construction of new roads and buildings had not yet begun.

Map data used in the study included information extracted from the Topographic Database maintained by the NLS.³ Buildings and roads from early 2015 were used in change detection tests with multispectral ALS data. Water areas and agricultural polygons were used to mask these areas out of analyses. Water areas were excluded from multispectral ALS and Sentinel-2 analyses. Agricultural areas were excluded from Sentinel-2 analyses. Road vectors in the database represent the center lines of the roads, and the width of the roads is not available. In change detection, we used buffered road vectors. The widths of the buffers were based on approximate information available on the widths of different types of roads. Different versions of the topographic database data, and openly available map and aerial image data from the city of Espoo (map service⁵⁸) were also used as an aid at various stages of the study. It should be noted that due to the active construction work in the area, we did not have map data that exactly corresponded to the situation at the time of data acquisitions.

3.2.

Multispectral ALS Data

Multispectral ALS data from two dates, August 21, 2015, and June 14, 2016, were used in the study. The trees were in leaf in both of these datasets. The point density in preprocessed data was about 8 to $9\text{points}/{\mathrm{m}}^2$ for each channel in the August 2015 data and about 11 to $13\text{points}/{\mathrm{m}}^2$ in the June 2016 data. The main steps in preprocessing of the August 2015 data included basic preprocessing, the range correction of intensity values, cutting out overlap points, removing some error points, and the creation of digital surface models (DSMs), DTMs, and three-channel intensity images in raster format (for more details, see Ref. 10). The same procedures, except the procedure for DTM generation, were also applied to the June 2016 data. The pixel size of the DSMs and DTM used in the study was 1 m on the ground. The pixel values in the DSMs represent the maximum height of points from all channels. The pixel size of the three-channel intensity images was 0.2 m, and the intensity values represent the average intensity of first-pulse and only-pulse laser points divided by 100.

An additional preprocessing step compared to those described in Ref. 10 was a calibration to adjust the intensity information from June 2016 so that it would be comparable with the August 2015 data. Originally, the intensity levels from the two data acquisitions were not comparable, which means that direct change detection between the intensity values is not possible without adjustment. The adjustment used 18 natural and manmade targets found in the study area. The targets included building roofs, parking places, sports fields, and beaches. The ALS points falling to each of the areas were inspected to only include target hits, not, for example, cars or reflective lines on an asphalt-covered car park. One of the targets was discarded because its behavior in the time series data was anomalous. For each channel, linear fit with scalar adjustment was found between the mean intensities from the two dates. Intensities in the later dataset were then adjusted with this linear model. To avoid extrapolation and negative intensity values, values below the minimum mean intensity value of the reference targets were not adjusted. According to visual evaluation, the intensity levels in the two datasets matched each other relatively well after the adjustment.⁵⁷

In addition to the preprocessed height and intensity data, land cover classifications carried out for the data^10,57 were used in the study. The land cover classifications were carried out using object-based random forest⁵⁹ classification. The classification results from the first date, together with the DSMs and intensity images from both dates, were used as input data in change detection tests based on multispectral ALS data. The classification results from both dates were used as reference data in evaluating change detection results based on Sentinel-2 satellite images.

3.3.

Sentinel-2 Satellite Image Data

All Sentinel-2 images from our study area were downloaded from the Copernicus Open Access Hub (previously known as the Sentinels Scientific Data Hub⁶⁰). Starting from July 15, 2015, and concluding on June 30, 2016, a total of 44 Sentinel-2 images were downloaded. Many of these images had a high percentage of cloud cover, and only nine images were finally used in the study (Table 1). The selection criterion for images was that they should be completely cloud-free over our study area.

Table 1

The acquisition dates of the remotely sensed datasets used in the study.

Sentinel-2 data at the time of image acquisition were only available in Level-1C format,⁶¹ so it was necessary to manually produce a level-2A product, which represents bottom-of-atmosphere information (i.e., atmospherically corrected data). For this purpose, the sen2cor tool provided by the ESA⁶² was used. Nowadays, users can also download level-2A data directly from the hub.

Having all the images processed to the level-2A product level, the calculation of a normalized difference vegetation index (NDVI) and normalized difference built-up index (NDBI) was done for each image. NDVI was selected because it is the most frequently used vegetation index and it easily distinguishes healthy vegetated areas, whereas NDBI is suitable for highlighting urban areas. Appropriate Sentinel-2 bands and the official software SNAP (Sentinel-2 Toolbox)⁶³ were used to calculate the indices:

3.4.

MLS Data

The MLS data for the study were collected using a Finnish Geospatial Research Institute (FGI) proprietary MLS system called Roamer. Roamer used a Riegl VUX-1HA laser scanner to collect 3-D data of the scenery. The frequency of the scanning mirror was 250 Hz, resulting in cross-track profiles; the scanner emitted 1 017 000 pulses/s. In the online signal processing, the scanner could resolve for up to seven echoes per pulse, depending on the object characteristics. Each 360 deg profile corresponded to 4068 transmitted pulses, resulting in an angular resolution of 1.5 mrad. The point spacing was $\sim 4\mathrm{mm}$ right below the scanner due to the sensor elevation of about 2.8 m from the ground. The point spacing on a plane perpendicular to the laser rays was 77 mm at a 50-m distance to the scanner. The positioning of the MLS system was carried out using the NovAtel SPAN GNSS-IMU (Global Navigation Satellite System—Inertial Measurement Unit; Flexpak6 receiver and ISA-100C IMU) positioning system for raw observables. The trajectory and orientation data were postprocessed with base station data from Trimnet virtual data service using Waypoint Inertial Explorer (NovAtel Inc., Calgary, Canada). The data were then georeferenced into a 3-D point cloud using Riegl RiWorld software (Riegl, Horn, Austria) and using previously acquired system offsets and calibration. The point cloud data were stored in LAS 1.2 format for delivery and processing.

The data used in the study were acquired on May 9, 2017, when trees were not yet in full leaf. Roamer was driven mainly on streets and, in some cases, on pedestrian routes in order to cover the area, buildings, and construction sites of interest. An excerpt of the data is presented in Fig. 2, showing the 3-D point cloud colored by the pulse return amplitude in tones of gray.

Fig. 2

An excerpt of the MLS data shows a building under construction. The data capture lots of details for interpretation.

4.1.

Change Detection from Multitemporal Multispectral ALS data

4.1.1.

Change detection method

The objective of the change detection test was to study the feasibility of direct height and intensity comparisons between multitemporal multispectral ALS datasets for map updating. A change detection process relying on two datasets was developed and mainly implemented with the eCognition software (Trimble Germany GmbH, Munich). The basic steps of the method are illustrated in Fig. 3. The main input data for the method include multispectral intensity images and DSMs from the two dates, and the results show the changes that have occurred in height or intensity. The land cover classification results from the first date are used as a starting point for the analysis. This ties the detected changes to land cover classes and also allows the use of different change detection rules for different classes. The land cover classes were generalized for the purpose of change detection (see Table 2).

Fig. 3

The main steps of the method used for change detection from multitemporal multispectral ALS data. Pseudo-NDVI values for segments were calculated from their mean values in channels 2 and 3. The height values of segments were their mean values in the maximum DSM that represents the highest elevation values for individual pixels.

Table 2

The relation between land cover classes in the suburban study area10 and change detection from multitemporal multispectral ALS data.

Class from land cover classification	Input class in change detection based on intensity	Input class in change detection based on DSM	Output classes from change detection
Building	—	Building	Building, height increase
			Building, height decrease
Tree	—	Tree	Tree, height increase
			Tree, height decrease
Asphalt	Artificial ground	Ground	Artificial ground, NDVI increase
			Ground, height increase
			Ground, height decrease
Gravel	Artificial ground	Ground	See above
Rocky area	Natural ground	Ground	Natural ground, NDVI decrease
			Ground, height increase
			Ground, height decrease
Low vegetation	Natural ground	Ground	See above

The most typical changes in our study area include the felling of trees due to the beginning of construction work and various other changes related to the construction of new buildings and roads. The felling of trees and the construction of new buildings are basically clearly visible in the height data. The same applies to changes in existing buildings, such as their expansion, demolition, and changes in height. The change detection of objects previously classified as “buildings” and “trees” was thus studied on the basis of the two DSMs. Height changes in ground objects were also analyzed because these can reveal new buildings and the extensions of previous ones.

Changes in ground-level objects, such as the construction of new roads, may also be visible in the height data, but it can be expected that intensity changes are the most important cue in this case. Changes in ground-level objects were thus also studied on the basis of the two intensity images, in practice, pseudo-NDVI values calculated from the images. A pseudo-NDVI can be calculated from Optech Titan channels 2 and 3.⁹ The aim was to detect changes from natural to artificial ground, indicating new roads, parking places, and construction sites. The opposite case (i.e., change from artificial to natural ground) was also included. This is a less typical change in practice, but it can occur, for example, when construction work is finished and newly sowed grass begins to grow on roadsides or in courtyards.

The most typical uninteresting changes detected with our approach include large vehicles, such as buses on roads, and temporary buildings or other constructions not presented on the map. Special rules exploiting the existing map data were developed to recognize these cases. The parameter values used in change detection and the reasoning for the selection of these values are shown in Table 3. Generally, the parameter values were based on general knowledge or analyses and experiments carried out in our training area. For extending the method to larger areas and different seasons, more sophisticated learning approaches could be utilized.

Table 3

The parameter values used in change detection from multitemporal multispectral ALS data.

Parameter	Value	Reasoning
Threshold values for height changes	Less than $-2.5$ m or greater than 2.5 m	The expected minimum height for a building or one floor in a building; trees are also normally higher than this
Threshold values for pseudo-NDVI changes	Less than $-0.2$ or greater than 0.2	According to histogram analysis in the training area, intensity changes in ground classes were normally within the range of $-0.2$ to 0.2
Threshold value for eliminating height changes with small object size	$<20{\mathrm{m}}^2$	Buildings are normally larger than this; the threshold value helps to avoid a large number of very small changes in the results
Threshold value for eliminating intensity changes with a small object size	$<200{\mathrm{m}}^2$	Experiments were carried out in the training area to find a suitable value; ground-level objects normally presented in the topographic database are also clearly larger than this
Threshold values for separating vehicles from the classes “ground, height increase” and “building, height decrease”: Maximum overlap with buffered road objects from the map, Object size	$>30\%$$<150{\mathrm{m}}^2$	Experiments and object size measurements were carried out in the training area; it should be noted that moving vehicles can appear as elongated objects in ALS data
The criterion for separating temporary buildings from the class “building, height decrease”: The number of overlapping building objects on the map	0	These buildings were not presented on the map; it is thus likely that they were building-like objects or temporary buildings that are not of interest for mapping

The land cover classification method and the quality of the results have been described by Matikainen et al.¹⁰ Segmentation was carried out with the multiresolution segmentation method,^65,66 and the parameter values were the same as those that were used in land cover classification. The multiresolution segmentation method is a widely used segmentation method based on bottom-up region merging. It uses an optimization procedure to minimize the average heterogeneity of segments.⁶⁶ By varying parameters and input data, users can create segments with different sizes and shapes. These segmentation results can form a hierarchical structure, where segments on an upper level consist of one or several segments on a lower level. In our ALS change detection test, we created only one segment level based on the DSM data and one based on the intensity data and these were separated from each other. The segment hierarchy, however, was utilized to create segments inside certain land cover classes and to merge neighboring segments with similar changes before analyzing the sizes of changed areas.

4.1.2.

Change detection results

The method development, including the setting of parameter values, was carried out in a $2\mathrm{km}\times 2\mathrm{km}$ subarea (about $3.9{\mathrm{km}}^2$ excluding water). Another $2\mathrm{km}\times 2\mathrm{km}$ subarea (about $3.4{\mathrm{km}}^2$ excluding water and a small land area outside the intensive study area) was used for testing and evaluating the usefulness of the detected changes. The results for this test area are shown in Fig. 4. Figure 5 shows some details of the results, together with the DSMs, intensity images, and land cover classification results from the first date.¹⁰ The figures include changes that are potentially of interest for map updating (i.e., changes classified as “vehicles” or “temporary buildings” were excluded).

Fig. 4

Change detection results based on multitemporal multispectral ALS data from August 2015 and June 2016. The land cover classes correspond to the classification of the August 2015 data.¹⁰ Buildings and road center lines from the map are based on data from the NLS Topographic Database 2015. The boundary of the study area is mainly based on a water mask that also contained data from the NLS Topographic Database 2015.

Fig. 5

Close-up views of the multitemporal multispectral ALS change detection results for (a) subarea 1 and (b) subarea 2. Left column: DSMs from August 2015 and June 2016. Middle column: intensity images from August 2015 and June 2016. Right column: land cover classification results from August 2015¹⁰ and change detection results between August 2015 and June 2016. For the legend of the change detection results, see Fig. 4. Buildings and road center lines in the change detection results are based on data from the NLS Topographic Database 2015.

4.2.

Continuous Change Monitoring with Sentinel-2 Satellite Images

4.2.1.

Change monitoring method

The method developed for continuous change monitoring from Sentinel-2 satellite images consists of preprocessing, segmentation, and change detection. An overview of the method is presented in Fig. 6.

Fig. 6

Continuous change monitoring from Sentinel-2 satellite images.

Preprocessing included atmospheric correction and the calculation of an NDVI and NDBI, as described in Sec. 3.3. Segmentation was performed in a hierarchical manner using the multiresolution segmentation algorithm^65,66 and appropriate bands and parameter values (see Table 4) in the eCognition software. The basic idea of this hierarchical segmentation method was already described in Sec. 4.1.1. In the Sentinel-2 change detection test, the first and the last image of the time series were used as input data. The first image was first segmented to obtain initial segments on a higher hierarchical level, and a lower hierarchical level of segments was then created inside the initial ones using the last image. In this way, boundaries already detected from the first image remained and new ones only detected from the last image were added. The final segments on the lower level were analyzed in change detection. The boundaries of changed areas should be visible in these segmentation results. Water and agricultural polygons from the topographic database were used as thematic layers in segmentation and excluded from further change analysis. Agricultural areas have their specific behavior considering spectral changes, and they require separate rules for change detection. In our intensive study area covered with reference points, there were only a couple of agricultural areas, and therefore, this topic was not studied in our case. Generally, most of the nonbuilt-up land in Finland is covered by forest, and the study could be later extended to agricultural areas. The parameters for segmentation were selected based on visual analysis and interpretation of segmentation results while varying their values. A total of 1686 segments (polygons) were created. The segmentation results for a subarea can be seen in Fig. 7.

Table 4

The parameters used in the segmentation of Sentinel-2 images for change detection.

Parameter	Segmentation level I	Segmentation level II
Input level	Pixels	Segmentation level I
Source image	August 17, 2015	June 29, 2016
Image bands in use	Red, green, blue, near-infrared (each with equal weight)	Red, green, blue, near-infrared (each with equal weight)
Thematic information	Water and agricultural areas	Water and agricultural areas
Scale	100	100
Shape	0.6	0.6
Compactness	0.8	0.8

Fig. 7

The segmentation of Sentinel-2 images for change detection. The Sentinel-2 satellite image in the background is from June 29, 2016 (Copernicus Sentinel data [2016]).

To develop rules for the change detection, training segments with no changes during the time frame were analyzed. There are natural changes in reflectance values during the year, such as the normal, expected changes in vegetation, and the purpose of analyzing the training segments was to obtain an understanding of these changes. The basic assumption was that changes in interest when mapping are larger than these normal seasonal changes. Training segments were selected on the basis of reference points that remained unchanged (see Sec. 3.1). Every segment that had at least one reference point within its boundaries was selected as a training segment. The total number of training segments was 100, and all of them were considered unchanged.

For each of the training segments $S_i,i\in \{\mathrm{1,2},\dots ,100\}$, the median value was calculated from NDVI and NDBI bands for each date $D_n,n\in \{\mathrm{1,2},\dots ,9\}$ to represent the segment reflectance. Median values were used instead of mean values in order to avoid the influence of outliers and averaging in the case of mixed segments. The differences of median values between every two consecutive dates were then calculated for the segments (${\text{MEDIAN}}_{S_i{D}_{n+1}}-{\text{MEDIAN}}_{S_i{D}_n}$). For each time fragment (two consecutive dates), the minimum and maximum segment differences were found, and these were the threshold values for change detection rules:

Eq. (3)

$$T_{\min D_{n+1}{D}_n}=\mathrm{MIN}({\mathrm{MEDIAN}}_{S_i{D}_{n+1}}-{\mathrm{MEDIAN}}_{S_i{D}_n}),i\in \{\mathrm{1,2},\dots ,100\},n\in \{\mathrm{1,2},\dots ,9\},$$

Eq. (4)

$$T_{\max D_{n+1}{D}_n}=\mathrm{MAX}({\mathrm{MEDIAN}}_{S_i{D}_{n+1}}-{\mathrm{MEDIAN}}_{S_i{D}_n}),i\in \{\mathrm{1,2},\dots ,100\},n\in \{\mathrm{1,2},\dots ,9\}.$$

Time series statistics, including the average, minimum, and maximum values of the median differences for each time fragment, are shown in Fig. 8. Change detection for a given segment $S_j,j\in \{\mathrm{1,2},\dots ,1686\}$ was based on the threshold values:

$$\begin{gathered} \text{if}[({\mathrm{MEDIAN}}_{S_j{D}_{n+1}}-{\mathrm{MEDIAN}}_{S_j{D}_n})<T_{\min D_{n+1}{D}_n}\text{or} \\ ({\mathrm{MEDIAN}}_{S_j{D}_{n+1}}-{\mathrm{MEDIAN}}_{S_j{D}_n})>T_{\max D_{n+1}{D}_n}] \\ \text{then Change} \\ \text{else No Change} \end{gathered}$$

Fig. 8

(a) NDVI and (b) NDBI time series statistics that were calculated for the training segments.

An additional rule applicable to an NDVI was based on the assumption that the NDVI value should decrease if a change from a vegetation area to a built-up area occurs. This was the change we were mostly interested in. In our case, this additional condition was met by only using the first inequality of the criterion above. This was possible because all the minimum threshold values ($T_{\text{min}}$) were negative, that is to say, they were related to an NDVI decrease, and all the maximum threshold values ($T_{\text{max}}$) were positive, that is to say, they were related to an NDVI increase.

4.2.2.

Change monitoring results

After having set the rules using the training segments, the median value of the NDVI and NDBI was calculated and change detection was performed for all 1686 segments. The change detection results based on the NDVI decrease rule are shown in Fig. 9. The results shown in the figure are a combination of all the changes detected within the selected time span (August 2015 to June 2016). A total of 432 polygons were detected to have changed.

Fig. 9

(a) Changed segments detected from Sentinel-2 images based on the NDVI decrease rule. The large yellow/orange polygon shows the coverage of the area that was used to select segments for the quality evaluation, and the smaller orange square shows the coverage of the subarea that is shown in (b), (c), and (d). The Sentinel-2 images in the background are from August 17, 2015 (b) and June 29, 2016 (a, c, d) (Copernicus Sentinel data [2015 and 2016]). A water mask that contained data from the NLS Topographic Database 2015 was used to define the boundary of the study area. The area of the smaller square with changes detected from multitemporal multispectral ALS data and some topographic database data can be seen in Fig. 4.

4.3.

Detailed 3-D Mapping with MLS

4.3.1.

Methodology for urban scene classification

Our method for classifying MLS data consists of three steps: data preprocessing, ground removal, and classification.⁶⁴ In the first phase, we remove noise and partition data into blocks. Then, ground points are removed. Lastly, above-ground points are clustered using connected component analysis and these clustered components are classified into different categories. The workflow of the method is shown in Fig. 10.

Fig. 10

The workflow of street scene classification from MLS data.

Data preprocessing

On occasion, the scanner produces noise, in other words, points in the air, which are mainly due to multiple reflections. Long-distance echoes, that is to say, objects at a longer range than the given measurement rate permits, are also not solved for the correct range. Presumably, some misinterpretations are caused by scattering from aerosols (dust) and other scanner online echo detection faults (e.g., direct or reflected sunlight entering the scanner optics). The noise causes mistakes in the following stages if not treated properly. Because of the low intensity of these noisy points, a low-intensity threshold was used to remove them. In this paper, this threshold was empirically defined to be 1200.

MLS datasets can be very large and the processing can be computationally expensive if we process the whole dataset in one time. Additionally, there are many urban objects that are scanned multiple times, which will cause shifts of objects. In order to solve these problems, the data are partitioned into blocks based on trajectory points and their sequenced time. Then, these generated, sequenced blocks are cropped into a regular size along the trajectory direction, 30 m along the trajectory, and 60 m across the trajectory.

Ground point removal

Compared with the points of other classes, ground points have smaller height differences locally. Therefore, ground points are extracted by filtering out points with large height differences in their local neighborhood. Here, the neighborhood is defined to include the 100 nearest neighboring points of a certain point. The threshold for a large height difference was set to be 0.15 m. In addition, we filtered out points that were too high compared with the height of a trajectory point at the same time instant. Ground points are removed and the rest of the points remain for the next phases of processing.

Clustering and classification

After ground points are removed, a connected component analysis⁶⁷ is carried out to obtain above-ground components. The distance threshold in the connected component analysis was set to be 0.5 m. Rules were defined and applied to classify the components into three classes: “buildings,” “trees,” and “pole-like road furniture.” Building and nonbuilding components are first detected in the following way. From every above-ground component, planar patches are extracted using surface growing. If there is at least one detected planar patch satisfying the following conditions, the above-ground component is classified as a building component: (1) the area of the patch should be large enough (the threshold was set to be $18{\mathrm{m}}^2$) (2) the normal of the detected planar patch should be parallel to a horizontal direction, that is, the angle between the normal of the detected plane and the horizontal plane should be less than 5 deg; and (3) the patch should mostly contain points with a prominent planarity, that is, points whose local neighborhood are situated approximately on a plane. The planarity $F_P$ for a point $p$ is calculated based on the eigenvalues ${\lambda }_1\ge {\lambda }_2\ge {\lambda }_3$ of the covariance matrix of the points in the neighborhood of $p$ using the equation $F_P=({\lambda }_2-{\lambda }_3)/{\lambda }_1$. The point $p$ is classified as prominently planar if planarity $F_P$ is larger than linearity $F_L=({\lambda }_1-{\lambda }_2)/{\lambda }_1$ and scatterness $F_S={\lambda }_3/{\lambda }_1$. In a building component, the ratio of points with a prominent planarity had to be larger than 0.65. Apart from the conditions 1 to 3 above, the horizontal distance between at least one pair of points in a detected vertical planar patch should be larger than 2.5 m in order to remove narrow nonbuilding segments.

Next, trees are extracted from the nonbuilding components. One sent laser pulse can reflect multiple times from vegetation because of gaps between small branches or leaves. In contrast, bare ground or a building facade can only give a single return for one transmitted pulse. In addition, trees are typically higher than other types of vegetation. Therefore, a nonbuilding component is classified as a tree if the following two conditions are met: (1) the percentage of multiple returns in the component is high enough ($>5\%$) and (2) the highest point of the component is 2 m or higher from ground level.

Last, pole-like road furniture are recognized from the remaining nonbuilding and nontree components. This is performed by first applying a slice analyzing method and then a cylinder mask method for each component. We first cut the component into slices in a vertical direction. If the horizontal diameter of any slice is larger than 0.5 m, the component is discarded; otherwise, it is considered as a pole candidate. Next, for each pole candidate, we calculate the shift of center points and the difference of the diameters between every two adjacent slices. If the shift of center points or difference in the diameter between any two adjacent slices is too large ($>0.05\mathrm{m}$), the slices are removed from the pole candidate. If the number of remaining slices is lower than 3, the component is removed from the pole candidates. For the remaining pole candidates, a cylinder mask⁴⁷ is used to perform a final check regarding if the candidate is pole-like road furniture. A detailed description of a cylinder mask can be found in Ref. 47; here, we use the center points of slices to define the cylinder axis. The main difference to the method of Ref. 47 is that we use unorganized point clouds, whereas in Ref. 47, the scan line information, also called profile or row and column information, is used. With the scan line information, neighborhood searches are faster, but such information is not always available, and thus, the method used in this paper is more generic.

4.3.2.

Classification results

Figure 11 shows objects of different classes detected from the MLS data. According to visual evaluation, most objects were correctly classified. Most misdetections come from connected objects such as trees connected with buildings. This is the main drawback of object-based classification. Some tree trunks with only a few branches were misclassified as pole-like objects. It is difficult to distinguish them from pole-like road furniture. Except for the classified objects, the rest of the objects remained unclassified with the method and were removed. For example, there were some trees that were visible in the original data but not in the final results because they were not recognized as trees with the set of classification rules in use.

Fig. 11

(a) Ground, (b) buildings, (c) trees, and (d) pole-like road furniture detected from the MLS data. Different objects are shown with different colors in the detection of buildings, trees, and pole-like road furniture.

5.1.

Potential of Multitemporal Multispectral ALS Data for Change Detection and Map Updating

When changes detected from the multitemporal multispectral ALS data are combined with map data, an illustrative presentation of the changes is obtained (Fig. 4). It can be seen that many changes related to construction work in the study area were detected. The most prominent of these was the felling of several forest areas before the beginning of construction work. These changes were classified as “tree, height decrease.” Several construction sites previously classified as “ground” are also visible in the results because the height of these areas has increased or their pseudo-NDVI value has decreased. All these changes are obviously interesting for map updating, although their direct linking to map objects is not straightforward. Buildings are normally presented on the map in their finalized form, not as construction sites. Forest areas are not presented as polygon objects on the present Finnish topographic maps because they are considered to belong to the background. Larger clear-cut areas, however, are presented. Many small changes are also visible in the results. The most typical of these is the felling of individual trees. Such changes are not needed for nationwide topographic mapping, but they are interesting for detailed city models that also represent trees. The short time interval between data acquisitions also affected the nature of the changes detected in this study, for example, there were many construction sites instead of finalized buildings.

To better understand the usefulness and nature of the changes detected from multitemporal multispectral ALS data, two evaluations were carried out. The first one was a numerical one and its objective was to test how well new objects interesting to topographic mapping were detected. It was based on the reference points collected on new buildings, new roads, and new clear-cut forest areas. The results are summarized in Table 5. They show that 95% of the changes were detected. Two buildings points and five road points were not detected as changes. All the clear-cut points were detected. The classification of the detected changes varied depending on the land cover class of the segments on the first date.

Table 5

A comparison of the changes detected from multitemporal multispectral ALS data to reference points.

The second evaluation was a visual one, and the objective was to understand the nature of changes that were detected by direct automatic comparison of the multitemporal datasets. The evaluation was based on visual inspection of the two multispectral ALS datasets with some help from existing map data. The information given by this evaluation could be used in further research and development. Basically, the detected differences can include changes that are interesting to mapping but also many other changes, which, from the viewpoint of map updating, could be considered as false detections. Categorical classification of the detected changes as desired and undesired ones was not included in this study. This topic is closely related to the definition of mapping rules for different objects and requires special attention if the method is further developed for operational use. A total of 10 changes were picked from each change class in the results, including vehicles and temporary buildings. If the number of changes in a class was smaller than 10, each change was evaluated. The changes were picked evenly from attribute tables listing the detected DSM and intensity changes. The results are reported in Table 6.

Table 6

A visual evaluation of changes detected from multitemporal multispectral ALS data.

A visual evaluation of the detected changes (Table 6) showed that most of them were real changes that had occurred in the area. However, the changes were not always important for mapping or the importance was difficult to determine. For example, many small changes had occurred in the industrial area. In an industrial area, there are containers, large vehicles, and other large manmade objects, and changes in such objects are visible in the change detection results. They can be difficult to separate from small buildings in aerial data, even visually. Seven out of the 96 evaluated changes were interpreted as false detections due to minor geometric differences in the data from two dates. For example, in some cases, the boundary of a large building was in slightly different positions in the two DSMs. This was probably related to the lack of laser points besides an obscuring building. In some cases, there was a real change in the area but the classification of the change was incorrect due to a misclassification of the object from the first ALS dataset. For example, parts of a building that was under construction in both datasets had been originally misclassified as tree, and these parts were now assigned to the class “tree, height increase.”

Some new buildings were present in both ALS datasets although they were still missing from the map data used in the study. These buildings were not detected as changes and they are missing from the final results, that is to say, they are presented neither as map buildings nor as changes. Some examples of these can be seen in Fig. 5(a). In a real mapping situation, it would be important that the contents of the map correspond exactly to the situation in the first dataset.

5.2.

Potential of Sentinel-2 Satellite Images for Continuous Change Monitoring

From the viewpoint of data integration in map updating, an important question related to the Sentinel-2 monitoring process is its ability to map changes that are meaningful for mapping and appear in the same areas as those detected from higher-resolution datasets. In our study, this was evaluated by comparing the results with the multispectral ALS classifications from the beginning and ending of the monitoring period. Reference changes were derived from the land cover classifications of the two multispectral ALS datasets from August 2015 and June 2016. The ALS results originally had six land cover classes, which were combined into three classes: “built-up area” (including the original classes “asphalt,” “gravel,” and “building”), “low vegetation” (originally “low vegetation” and “rocky area”), and “tree” (“tree”). The three classes were the same as those that were used to select unchanged reference points for the Sentinel analysis. Cross-classification of the August 2015 and June 2016 results was carried out, and a reference change map with changed and unchanged areas was obtained. Water areas and areas outside the ALS data coverage were excluded. Three different Sentinel-2 change detection results created with the three different change rules (NDVI change, NDBI change, and NDVI decrease) were first evaluated by using the entire ALS data coverage (i.e., also including the areas of training segments). The best result was then also evaluated with the training segments excluded.

In practice, the Sentinel-2 segment vectors were first converted into raster format with the same pixel size as the ALS classifications (0.2 m). Segments with more than 90% of their area inside the ALS data coverage were selected for further analyses. This included a total of 223 segments. For each segment, the percentage of changed pixels in the ALS change map was calculated. The segments were then divided into four change categories using 25% intervals based on this reference change percentage. From the Sentinel-2 change detection results, the number of changed segments in each of these categories was calculated and divided by the total number of segments in the reference change category. The results are presented in Table 7. Among the three different change detection rules, the one with an NDVI decrease as the change criterion gave the best correspondence with the reference results. The criterion based on NDBI change clearly gave the lowest correspondence. The results based on the NDVI decrease criterion were thus also evaluated with the training segments excluded. There were more changes in the western part of the study area outside the training points (see Fig. 9), which is reflected in these results. Overall, the correspondence between the changes detected with the NDVI decrease rule and the ALS land cover classifications was good. In the highest change category (75% to 100% of the segment changed according to the reference data), 100% of the segments were detected to have changed. In the lowest category (0% to 25%), the percentage was only 9%.

Table 7

A comparison of Sentinel-2 change detection results with reference change categories derived from two multispectral ALS classifications. The best results were obtained with the NDVI decrease rule, and they are highlighted bold in the table.

A more detailed analysis of the reference classifications showed that the most common changes inside the changed Sentinel-2 segments obtained with the NDVI decrease rule were a change from tree to built-up area (39.3% of changed pixels), a change from tree to low vegetation (34.9%), and a change from low vegetation to built-up area (18.7%). It can thus be seen that most of the changes detected from the Sentinel data were clearly related to construction. Visual analysis of the detected changes and reference changes was also carried out. The overall accuracy of the original ALS classifications was as high as 96%.⁵⁷ However, there were some false positive changes in the reference change map due to classification errors (e.g., noise in forested areas, sports fields with artificial cover, nonhomogeneous gravel areas). In Sentinel-2 images, small image shifts (1 to 2 pixels) and some special vegetation types (e.g., waterfront vegetation) caused false-positive changes. Such shifts occurred in the early stage of Sentinel-2A image acquisition, and it is expected that this effect will decrease in the future. A suburban area is a mixture of different built-up and vegetated surfaces, and for small segments, even small shifts can cause remarkable changes in reflectance. Also, the small temporal difference between the Sentinel-2 and multispectral ALS acquisitions caused errors in active construction sites.

Despite the minor issues discussed above, we can conclude that the presented Sentinel-2 change monitoring method is feasible for detecting changes in a growing built-up area. Being fully automated, this method allows an easy and affordable way of performing continuous change detection over large areas. In the future, further development of the method can be carried out, for example, to better analyze special cases, such as waterfront vegetation and to extend the analyses to agricultural areas.

5.3.

Potential of MLS Classification Results to Complement ALS Data in Mapping Urban Areas

The ability of the classified MLS point cloud to complement ALS data and provide finer details for mapping an urban area was evaluated by visual interpretation. The MLS classification test did not involve change detection. The comparison was thus based on the ALS data from the second date and concentrated on the appearance of buildings and poles in ALS and MLS data. These are important objects for detailed city models of urban areas.

Figure 12 shows a comparison of ALS and MLS data in a small area of active construction. The upper row shows a top view of the multispectral ALS point cloud and the MLS point cloud, and the lower row shows the results from building extraction. The buildings extracted from the ALS data are based on the classification of the June 2016 multispectral ALS data in raster format.⁵⁷ Figure 13 presents a side view of one building in the ALS data, classified MLS data, and combined ALS and MLS data. A comparison for poles is presented in Fig. 14.

Fig. 12

(a) ALS data and (b) MLS data from an area under active construction (top view; colors are based on elevation; black areas represent areas with no data; the MLS data only include classified points). (c) Buildings extracted from the ALS data (based on Ref. 57) and (d) buildings extracted from the MLS data (based on the present study).

Fig. 13

(a) ALS data from around a building (side view; all points in white). (b) Classified MLS points around the same building (buildings are in red, poles in blue, vegetation in green and the ground surface in orange). (c) Buildings, poles, and the ground from the MLS classification inserted into the ALS point cloud.

Fig. 14

A side view of poles from two subareas. (a) and (d) ALS data (colors based on elevation). (b) and (e) MLS data (colors based on elevation; only classified points are included). (c) and (f) Poles extracted from the MLS data and inserted into the ALS point cloud (inserted poles appear in blue color).

The comparisons show the complementary nature of the MLS and ALS point clouds. Building roofs are well presented in the ALS data, and the number of reflections from building walls is small. In MLS data, the situation is quite the opposite, and wall structures are very well presented. A combination of these two datasets would allow accurate modeling of both building roofs and facades. ALS and MLS data also provide two alternative ways for defining the outlines of 2-D buildings, ALS according to the roof outlines, and MLS according to the position of the walls. This might be useful, especially if map data present the bases of the buildings instead of roof edges, as for example the NLS Topographic Database.⁶⁸ In such cases, building outlines derived from MLS data could correspond to maps better than those derived from ALS data. However, MLS data acquired from roads also have deficiencies when considering building mapping, such as the lack of courtyards.

Poles are another interesting example, where MLS data could be used to complement ALS data. The tops of poles such as street lamps are often visible in dense ALS data, but the sides of the poles are not visible. In the MLS data, the poles are clearly visible and they can be extracted. Poles are not needed for nationwide topographic mapping, but they are interesting for detailed city mapping and modeling, street lighting design, and optimization,⁶⁹ and also for power line monitoring.⁷⁰