AUTOMATIC DETECTION OF STORM DAMAGES USING HIGH-ALTITUDE PHOTOGRAMMETRIC IMAGING

The risks of storms that cause damage in forests are increasing due to climate change. Quickly detecting fallen trees, assessing the amount of fallen trees and efficiently collecting them are of great importance for economic and environmental reasons. Visually detecting and delineating storm damage is a laborious and error-prone process; thus, it is important to develop cost-efficient and highly automated methods. Objective of our research project is to investigate and develop a reliable and efficient method for automatic storm damage detection, which is based on airborne imagery that is collected after a storm. The requirements for the method are the before-storm and after-storm surface models. A difference surface is calculated using two DSMs and the locations where significant changes have appeared are automatically detected. In our previous research we used four-year old airborne laser scanning surface model as the before-storm surface. The after-storm DSM was provided from the photogrammetric images using the Next Generation Automatic Terrain Extraction (NGATE) algorithm of Socet Set software. We obtained 100% accuracy in detection of major storm damages. In this investigation we will further evaluate the sensitivity of the storm-damage detection process. We will investigate the potential of national airborne photography, that is collected at no-leaf season, to automatically produce a before-storm DSM using image matching. We will also compare impact of the terrain extraction algorithm to the results. Our results will also promote the potential of national open source data sets in the management of natural disasters. * Corresponding author.


INTRODUCTION
The photogrammetric imaging is a cost-efficient method for collecting data for mapping purposes.Recent developments in digital surface model (DSM) generation from stereoscopic imaging by image matching techniques have made the method serious competitor for laser scanning in applications which utilize the object surface information (Leberl et al 2010).The imagery also provides visual, stereoscopic support for various forest storm damage management tasks.
Airborne laser scanning (ALS) provides accurate information about the canopy height, structure and the underlying terrain (Hyyppä et al, 2012a(Hyyppä et al, ,2012b) ) and is widely considered to be the most promising method for various forest inventory and management tasks.ALS data was used for damage detection in forest (Rahman 2011, Vastaranta et al 2011), both studies used the height difference between the ALS measurements.The availability of ALS data that is recorded right after storm is unlikely, because the data is costly and therefore the flights are well planned.
Storm damage detection process is characterized by the fact that storms can occur at any season, which also influences the properties and quality of the data used in remote sensing tasks.The most challenging conditions for matching forest surfaces occur during the seasons when the deciduous trees do not have leaves, which can lead to a decrease in the matching quality (Baltsavias et al 2008).Low solar altitudes during autumn and winter seasons can be challenging, especially due to shadows, which can deteriorate the quality of image matching (Honkavaara et al 2012).We expect that the methods that are the least sensitive to these natural variations have the highest automation potential.Methods that are based on comparisons of DSMs collected at two points in time (bi-temporal) could be efficient for automating storm damage detection.Other promising methods for automatic DSM generation include laser scanning (Hyyppä et al, 2012a(Hyyppä et al, , 2012b)), radargrammetry using SAR images (Karjalainen et al 2012), and the automatic matching of images (Leberl et al 2010, Baltsavias et al 2008, Hirshmüller 2011, Järnstedt et al 2012).
In our previous research (Honkavaara et al 2013) we used fouryear old airborne laser scanning surface model as the beforestorm surface.The after-storm DSM was provided from the photogrammetric images using the Next Generation Automatic Terrain Extraction (NGATE) algorithm of Socet Set software.We obtained 100% accuracy in detection of major storm damages.We identified several topics to be considered in our future investigations.The before-storm DSM is a critical requirement in the method.One of the most interesting sources for before-storm DSM is to produce the model with image matching from national high-altitude photogrammetric imagery.These data sets are available much more widely and with higher temporal resolution than airborne laser scanning based surface models.Another important issue was the influence of the terrain extraction algorithm in the results.We used area based approach, but very interesting new developments are methods, that are based on pixel wise cost functions, such as the the semiglobal matching (sgm) (Hirschmüller, 2011) expected to provide even denser DSMs.We will study these two topics in this investigation.We will describe the data and methods in Section 2, we give the results in Section 3 and conclusions are given in Section 4.

DATA AND METHODS
In this study, an area of size 21 km by 5 km was used.The area under study has a flat topography with maximum of 50 m height differences.The main land cover consists of coniferous forests; mixed forests and transitional woodland are also present.

Data
We used as the pre-storm data the national ALS data and the national photogrammetric imagery.As the after-storm data we used imagery collected after the storm.We have used in our studies a clipping of size 21 km by 5 km from one image strip of the after-storm flight

After-storm imagery
To collect after-storm imagery, airborne photogrammetric image flights were carried out using a Microsoft UltraCamXp largeformat mapping camera on January 8th, 2012, right after the first snow had fallen.The average flying height was 5370 m above ground level (AGL), which yielded a ground sample distance (GSD) of 32 cm.The image block was with16 image strips and approximately 30 images per strip; the forward overlaps were 65% and the side overlaps were 30%; the distances of the image strips were approximately 3900 m.The atmosphere was clear and the solar elevation was only 5-7°.
Total area covered during one day flight was 1620 km2.

Pre-storm ALS data
The ALS data were from the Finnish national land survey (NLS).The minimum point density of the NLS ALS data is half a point per square meter, and the elevation accuracy of the points in well-defined surfaces is 15 cm; the horizontal accuracy is 60 cm.The point cloud has automatic ground classification.
The ALS data used in this study was collected in 2008 during the spring.LasTools software (Isenburg 2013) was used to merge the map sheets of NLS laser data that covered the area and to make a DSM of the point cloud with 1 m grid spacing.The first and only pulses were used for DSM.

Pre-storm image data
Pre-storm photogrammetric imaging was done in 2010 at 28 th of April, when deciduous trees are gaining leaves.The flying height was approximately 7485 m above mean sea level giving a GSD of 45 cm.The camera used was a Vexcel Imaging UltraCam Xp, and pansharpened level 3 images were used.The forward overlap was 60% and the side overlap was 32%.
In the case of after-storm data, the orientations were calculated for a sub-block with five image strips.26 ground control points were taken from the ALS point cloud and the ALS intensity image.Automatic tie points were then extracted and the orientations were determined using a bundle block adjustment.
The registration quality of the ALS and photogrammetric data was estimated to be approximately 1 m.Details of the processing are given by Honkavaara et al.,(2013).
Orientations of the before-storm images were determined for a small block with three image strips and 9-13 images per strip.
For this image block we used altogether 11 ALS XYZ GCPs with a priori standard deviations of 1 m in all coordinates.We extracted 145 tie points per image by automatic image matching.We did not have a priori orientation information for this image block.The standard error of unit weight was 0.202 pixels and the root-mean-square residuals were 0.440, 0.678 and 0.768 m for X, Y and Z, respectively.
SOCET SET Next Generation Automatic Terrain Extraction software (NGATE) (DeVenecia et al 2007, Socet Set 2009) was used to generate DSMs.The matching strategy employs the edge and correlation matching methods and applies image pyramids and back matching.Different predefined and selftuned strategies were tested.Finally, one of the default NGATE strategies, "ngate_urban_canyon.strategy," was selected as the method to be used in this investigation, because it appeared to provide the highest success rates in matching in forest areas.DSMs were generated in an irregular point-network (TIN) format of the NGATE software, with a point interval of 1 m and with a no-thinning option.The process is explained in more detail in (Honkavaara et al. 2013).This methodology was used to process after-storm and before-storm imagery into surfaces that are referred to as NGATE surfaces in the following parts of the study.
Another after-storm surface was computed by Blom Kartta Oy using a dense matching algorithm of the Microsoft UltraMap software.The UltraMap software uses a dense matcher algorithm to provide point clouds (or DSM) with the same resolution as the image data (Wiechert et al., 2012).The UltraMap point cloud has point spacing of 32 cm according to the imagery GSD.To be able to compare the surfaces from the two photogrammetric point clouds, we have thinned the points with LasTools using the thin_with_grid filter with grid size of 1m before computing the surface grid for the point cloud.The surfaces generated with this method are referred to as UMD (UltraMap dense) surfaces in the following parts of the study.

Computation of the surface grids
The LasTools blast2dem (Isenburg, 2013) function was used to make a grid DSM of the point cloud with a 1 m grid spacing.All of the point clouds from different sources were processed with the same function and the same size, grid cell size, and limiting lower left coordinates of the resulting DSM were given.Further registration was not performed.

DIFFERENCE SURFACES IN DAMAGE DETECTION
Four different difference surfaces were computed.2 and 3.An intact area in forest is marked with 1.In the area marked with 2, the pre-storm NGATE surface is high and the others low, this could be a place where windfall has occurred.In the area marked with 3, thinning has probably taken place between the ALS and the pre-storm NGATE, and the storm has finished most of the trees that were left after thinning.The area marked with 4 shows how the UMD surface forms steeper edges than the NGATE surface.
In Figure 4, a profile of the example area in Figures 2 and 3 is shown.The profile goes from bottom to top in the forested area to the left in Figure 2. In Figure 4, the ALS data has the largest variability in height.The UMD surface has more height variability than the NGATE surfaces, which tend to keep to the top part of the canopy.

Classification
We have developed a classification method for the difference surface between ALS and NGATE surfaces in (Honkavaara et al 2013).Using this method, we computed pixel wise classifications for each of the difference surfaces in Figure 1.
The method classifies the surface into three classes: 'no change',' sparse change' and 'large change'.We studied the difference surfaces together with classified surfaces in 70 squares of 100 x 100 m that are located in forested areas in the ALS data.We had a visual classification of the squares from the after-storm images.

Performance assessment
The classifications of Figure 5 were used in Table 1 so In Table 1 and Figure 5, it can be seen, that the performance of the classification method is far from optimal for the UMD surfaces a) and c).The number of changed pixels is too large for usability.The differences in the logged areas in Table 1 are caused by the fact that three of the eight areas marked as logged using the after-storm imagery had been logged before the prestorm imagery was taken, thus there is no change in the height between the pre-and after-storm surfaces c) and d) in these areas, and the 63% result in the Table 1 is the best that can be achieved.

Height distribution
The height distribution in the four different difference surfaces was compared by computing the percentage of points above height 3 m in the difference surface in the test squares.The test squares were grouped to five different groups according to the level of damage that was visually apparent on the square.The damage level groups and the number of squares in them are given in Table 1, leftmost column.It is evident from Figure 6, that with large number of fallen trees, the number of high height values in the difference surface increases and this is the base for the automatic storm damage detection process.The intact area 'no damage' in Figure 6 has both lower mean and standard deviation in the percentage of heights above 3 meters in the test areas.These test squares are clearly distinct from the squares, where either windfall or logging has occurred.

CONCLUSION
The results of this study show, that the storm damages are visible on the difference surfaces regardless of the type of the surface or surface generation technique.Because of this, the methods that are based on comparisons of surface models are very promising for the automating the storm damage detection process.The main difference of results between the studied photogrammetric surface generation methods Socet Set NGATE and Ultra Map dense matching (UMD) was that the UMD surface made steeper height changes than the NGATE.This was visible in the forest areas.The difference in the surfaces did not affect the visibility of the 'large change' areas (>10 fallen trees in ha) in the difference surfaces.The UMD appeared to provide slightly better results in cases where there were less changes in the forests (less than 10 fallen trees in ha).
The used classification method was designed with a difference surface ALS -NGATE.The method performed surprisingly well with the tested materials, but we expect that better results could be obtained by selecting the parameters based on the data set.Automating the parameter selection process will be of importance.Our future study in this topic will further develop the classification method.So far, a simple approach which mainly demonstrates the usability of the surface difference has been used.With a more sophisticated method, we expect better results.
Important further evaluations include the impact of DSM generation method to the result as well as further tuning of parameters of the change detection method.However, the data was with 44 -45 cm GSD and collected from a flight altitude of 7400 m above the terrain level so this is very extreme data to be used.Furthermore, the forward overlaps of the data were 60%, which is not ideal for dense surface model generation.
As the final conclusion we would like to emphasize that the high-altitude photogrammetric imagery and automated digital surface model generation are promising and cost-efficient tool for operational, automatic storm damage detection process.The further advantage of the imagery is that it can be used for training the automatic damage detection method and for quality control purposes.
Figure 1.Difference surfaces from left: a)ALSafter-storm UMD, b) ALSafter-storm NGATE, c) pre-storm NGATEafter-storm UMD, d) pre-storm NGATEafter-storm NGATE.The colour bar at the bottom shows the height colour coding.Different nature of ALS and photogrammetric surfaces is useful in the interpretation of the difference surfaces of Figure 1.In Figure 1 a) and b), the forested areas are visible, and this makes the difference surface usable without map.In Figure 1 c) and d), both of the surfaces are photogrammetric, and only the forest edges are visiblethe division between field and forest patches is not as clear as in the Figure 1 a) and b).

Figure 2
Figure2illustrates a small area with a variety of different items: fallen trees, logging, field, and forest are shown in the afterstorm and pre-storm images.In Figure2, the poor imaging conditions in the after-storm image are visible.The shadows are very long due to the extremely low solar elevations (5-7 o ) and the ground is covered with snow.The images show a logged area on the right side of the picture.

Figure 3 .
Figure 3. Difference surfaces from the area of Figure 2. a) ALS after-storm UMD b) ALSafter-storm NGATE c) pre-storm NGATEafter-storm UMD d) pre-storm NGATEafter-storm NGATE.The colour bar in the middle shows the height colour coding.The essential information on the height change from damage view-point is the same in all of the comparisons between prestorm and after-storm surfaces in Figure3.There are changes in the detail and especially in the visibility of small changes.Surfaces in Figures3c) and 3d) that show comparison between photogrammetric models show the edges of the forest areas visible as lines in both UMD and NGATE surfaces.This is probably a sum of multiple factors including different viewing angles and heights, registration errors in the imagery and different lighting conditions.The ALS and the photogrammetric surfaces are different in the penetration to the forest ground.In ALS data, there are more hits from the forest floor inside the forest than in photogrammetric data.When the NGATE and UMD surfaces are compared, the UMD surface is closer to the ALSit shows sharper tree and forest edges than the NGATE surface.

Figure 5 .
Figure 5. Classification of the surfaces: blue -'no change', green -'sparse change', dark red -'large change'.Surfaces from left: a) ALS-after-storm UMD, b) ALSafter-storm NGATE, c) pre-storm NGATEafter-storm UMD, d) prestorm NGATEafter-storm NGATE The main steps of the classification are enhancement of pixels with height change (large in earlier data and small in later data), averaging and thresholding.Two different sets of classifications with different parameters are made and then combined to find both sparse-and large changes.The classified surfaces are , that both 'sparse change' and 'large change' were combined into 'change' and the classification of a test square was 'change' if over 10% of the pixels inside were classified to 'change'.The upper figures in Table 1 are for class 'change', and the lower ones for 'large change' with the same 10% of the square criterion as for the change class.Table 1.Percentage of test squares classified as changed with different surface combinations.a) ALSafter-storm UMD, b) ALSafter-storm NGATE, c) pre-storm NGATEafter-storm UMD, d) pre-storm NGATEafter-storm NGATE.The upper figures are results for classification into any change class, lower bold figures apply for classifications to the class 'large change'.
Figure 6.Mean percentage of height values above 3m in test squares of different damage groups.The test squares are coloured with different colour for each surface.
, which are