3-DIMENSIONAL PRE-AND POST-FIRE COMPARISON OF FOREST AREAS

Forest agencies give special attention to forest fires where post-disaster loss can rarely be gauged in a quick and economic way unless an appropriate technology is adopted. Determination of the planimetric and volumetric changes between preand post-fire is in high demand. The FORSAT system (a satellite processing platform for high resolution forest assessment) was developed to meet the relevant demands. It has the capability of extracting of 3D geometric information from the very-high resolution (VHR) imagery from satellite optical sensors and automatic 3D change detection. FORSAT includes two main units. The first one is dedicated to the geometric and radiometric processing of satellite optical imagery and 2D/3D information extraction. This includes: image radiometric pre-processing, image and ground point measurement, improvement of geometric sensor orientation, quasi-epipolar image generation for stereo measurements, digital surface model (DSM) extraction by using a precise and robust image matching approach specially designed for VHR satellite imagery, generation of orthoimages, and 3D measurements in single images using mono-plotting and in stereo images. FORSAT supports most of the VHR optically imagery commonly used for civil applications: IKONOS, OrbView – 3, SPOT – 5 HRS, SPOT – 5 HRG, QuickBird, GeoEye-1, WorldView-1/2, Pléiades 1A/1B, SPOT 6/7. The second unit of FORSAT is dedicated to 3D surface comparison for change detection. It allows users to import DSMs, to co-register them using an advanced 3D surface matching approach and to calculate the planimetric and volumetric changes between epochs. The capacity and benefits of FORSAT have been tested in two real cases, where are burned areas located in Cyprus and Austria. The geometric characteristics of burned forest areas have been identified both in 2D plane and 3D volume dimensions, using preand post-fire optical data from different sensors. FORSAT is a single source and flexible forest information solution, allowing expert and non-expert remote sensing users to monitor forests in three and four dimensions from VHR optical imagery for many forest information needs. * Corresponding author


INTRODUCTION
Deforestation is one of the major sources of carbon emission which threatens the global climate targets.Several satellites and sensors have been used to monitor deforestation areas.The Advanced Very High-Resolution Radiometer (AVHRR) on board the NOAA-series satellites (Di Maio Mantovani and Setzer, 1997), JERS-1 (Almeida-Filho et al., 2005), MODIS onboard the NASA EOS satellites (Anderson et al., 2005), ASTER (Haboudane and Bahri, 2007), Formosat-2 (Baillarin et al., 2008), PALSAR on-board ALOS (Isoguchi et al., 2009) constitutes a small set of examples.Low to medium resolution level optical images present drawbacks for operation in the moist tropics and in all weather conditions, synthetic aperture radar (SAR) data might be seen as an alternative (Santos et al., 2008;Solberg et al., 2013).Applications of SAR data to map deforestation are generally based on the assumption that undisturbed forests consistently exhibit higher radar backscatter than deforested areas.Depending on the stage of the deforestation process (slashing, burning and terrain clearing), this assumption is not always valid, and deforested areas may display a stronger radar return backscatter than primary forest (Almeida-Filho et al., 2007).Especially, new deforested areas are not unequivocally detected in some cases (Almeida-Filho et al., 2009).With over 30 years of directly comparable satellite observations now freely available, and new imagery being added to the archive every day, Landsat time series legacy affords novel opportunities for ecosystem mapping, environmental monitoring and comparative ecology (Pasquarella et al., 2016).It has been used to understand the space-time dynamics of deforestation over large areas, but at moderate resolution (Alves, 2002;Ichii et al., 2003;Bodart et al., 2011;Souza et al., 2013).
Measuring the aerial extent of deforestation for other than localized areas requires the use of fine resolution satellite data.An accurate determination of deforestation is very difficult to achieve by a random sampling analysis of Landsat or similar resolution data unless a very high percentage of the area to be studied is sampled (Tucker and Townshend, 2000).The Very High Resolution (VHR) satellite imagery is an alternative and effective solution (Mora et al., 2013).The list and metric capability of the VHR satellite imagery is given in Remondino (2011) and Sefercik et al. (2013).
Traditionally, forest monitoring and information gathering is done through expensive and time-consuming field studies in combination with aerial optical imagery and/or LiDAR surveys (Eva et al., 2010;Koch, 2010).This is followed by digital data analysis in image processing suites and geographic information systems (GIS) to derive basic parameters such as area coverage, species, height, volume, health, damage, change, deforestation, etc. Subsequently the map and statistical data is used for management decisions, implementations of regulations or simply for industry operations.Although the traditional forest inventory methods are the most accurate, they are neither agile nor economic.Alternative management strategies are required (Mondal et al., 2010).The Global Forest Watch (GFW) is a worldwide response (http://www.globalforestwatch.org)as an open-source web application to monitor global forests in near real-time.It is an initiative of the World Resources Institute (WRI) with partners including Google, Esri and many other academic, non-profit, public and private organizations.
Deforestation is caused by ever-increasing activities of the growing human population (Pahari and Murai, 1999), its density (Svancara et al., 2009) and agricultural colonisation (Millington et al., 2003).Their effects are seen in long terms.On the other hand, the forest fire, which is another main cause of deforestation, is a rapid event whose effects are seen in very short terms (Lee, 2008).A rapid, effective and economic way of change detection is required in order to understand the pre-and post-fire changes of forest areas.The VHR satellite imagery which allows stereo and triplet acquisitions at very fine spatial resolutions offers less expensive, faster and more agile remote sensing capacities than the alternative technologies, thereby providing a complementary and single source solution for such change detection tasks.The other advantages are no overflight permissions needed, an optimum ground coverage capacity versus spatial resolution, and repetition of image acquisitions on a certain area of interest until cloud coverage is free.We have identified this special requirement in the forestry sector, and carried a project to develop a new software platform and associated services specifically forest monitoring and management tasks.The project is entitled as "a satellite processing platform for high resolution forest assessment" (FORSAT), which is a research & development project cofunded by the national funding authorities of the participating partners under Eurostars and the European Commission.
The aim of the FORSAT project (http://www.forsat.eu) is to develop a standalone satellite-based monitoring capacity for 3D forest cover mapping and change detection applications.The existing precise processing capabilities of terrestrial and airborne techniques are transferred to VHR optical satellite imagery level, thereby a single source forest information solution is obtained.The FORSAT software platform, described in Section 2, enables high resolution thematic mapping of forest areas along with 3D canopy height model (CHM), thereby allows derivation of forest volume.Based on new image acquisitions as well as historical data, the system allows automatic 3D change detection of forest and non-forest areas along with change in area and volume.The FORSAT capability and applicability to forest fires were tested by executing case studies located in Cyprus and Austria.The results have generated the necessary facts and needed trust of public bodies and private organisations to use VHR satellite data for many forest information needs, which do not require centimetre level accuracy.In Section 3 the results achieved on test cases are reported and discussed, followed by final conclusions.

TECHNICAL APROACH
The FORSAT software allows analysing raw satellite imagery and extracting meaningful and quantitative information about the world's forests, such as area, volume and measurements of deforestation or even regeneration of a forest.The FORSAT software has a workflow which is graphically represented in Figure 1.The overall architecture comprises four building blocks (pre-processing, geo-referencing, digital surface models (DSMs) generation and 3D comparison), being each block tightly coupled in a software suite framework.This workflow is a combination of two main tasks.The first one is dedicated to the geometric and radiometric processing of satellite imagery and 2D/3D information extraction, that is: radiometric pre-processing, image and ground point measurement and sensor orientation, quasi-epipolar image derivation, DSMs extraction, orthorectification, 3D vector measurement.This unit supports most of the VHR optical imagery commonly used for civil applications, for example IKONOS, GeoEye-1, WorldView-1/2, SPOT-5/6/7, Pleiades-1A/1B, and it can be easily updated to similar images from future missions.The second task is dedicated to 3D surface comparison for change detection.It allows the users to import DSMs and digital terrain models (DTMs), to align them using an advanced 3D surface matching technique and to calculate the 3D volume differences.The technical approach comprises a bunch of interconnected algorithms whose details are given in Poli ( 2005), Zhang (2005) and Akca (2007) all of which are the doctoral studies performed at the photogrammetry group of ETH Zurich.

Metadata files
The VHR satellite images are provided together with their metadata by the vendors.Before applying the algorithms for the geo-referencing of the images, some operations are required in order to prepare the input data.The pre-processing includes both the analysis of the metadata files for the extraction of the required information and the radiometric improvement of the images in order to facilitate the point measurements (Poli, 2007).

Radiometric pre-processing
The performance of the image matching and feature extraction procedures depend on the quality and quantity of information carried out by the images.Compared to the traditional scanned 8-bit/pixel images, digital imagery from linear array sensors has better radiometric performance e.g. higher dynamic range and signal-to-noise ratio.Most of the linear array sensors have the ability to provide high quality digital images.However, some radiometric problems still have to be considered: poor image contrast, the image blur problems mainly caused by CCD line jitter, kappa jitter and motion blur and deficiencies of the lens system, image noise, and radiometric problems caused by the variations in the sensor view angle, the sun angle, shadowing, and the seasonal and the atmospheric conditions.These problems are usually beyond the control of the users.However, they have to be restored as much as possible.In order to reduce the effects of such radiometric problems and optimise the images for subsequent feature extraction and image matching, the image pre-processing methods have to be employed (Zhang, 2005).Although the gamma correction, contrast enhancement, histogram equalisation are trivial applications and can be found in many standard image processing software, the Wallis filter is a specific and powerful algorithm (Wallis, 1976).The filter forces the mean and standard deviation of an image to given target values.The FORSAT software includes the Wallis image pre-processing method.

Image pyramids
An image pyramid is a multi-resolution representation of the original image.It is used to speed-up the image matching computation while at same time keeping the finest spatial resolution of the final DSM output.The image pyramids are used for coarse-to-fine hierarchical image matching methods.At present, many of the modern matching algorithms are based on the image pyramids.With coarse-to-fine hierarchical strategy based on image pyramid representation, the matches obtained at a coarse resolution are used to guide and limit the search space for the matching of finer-resolution features.In this strategy, the usual way is to start matching at a low-resolution pyramid level, where the influence of image noise is reduced and coarse approximate values are sufficient to stay within the pull-in range of the matching procedure.In addition, in the lowresolution images, the regions of interest for correspondence analysis in levels of higher resolution can be found at low cost because irrelevant details are no longer available there.The computations are usually performed successively on each level of the hierarchy using the results of the higher level as approximate information (Ackermann and Hahn, 1991;Zhang and Gruen, 2006).FORSAT software generates the image pyramids after the pre-processing step, starting from the original resolution images.Each pyramid level is generated by multiplying a generating kernel and reduces the resolution by factor 3.

Geo-referencing
The FORSAT software uses the rational function model (RFM), a well-known non-rigorous (generalised) orientation method based on the rational polynomials functions.A RFM is the ratio of two polynomials derived from the rigorous sensor model and the corresponding terrain information, which does not reveal the sensor parameters explicitly.The VHR satellite images are supplied with only rational polynomials coefficients (RPCs) instead of rigorous sensor model parameters.In the RFM, the image pixel coordinates are expressed as the ratios of polynomials of object coordinates, which usually correspond to latitude, longitude and height values.The RFM is computed based on a rigorous sensor model (Figure 2).
With the given parameters of the rigorous model and by projecting evenly distributed image points into the object space, multiple-layer 3D object points can be computed and used as virtual (fictitious) control points.The control points are created based on the full extent of the image and the range of elevation variation in the object space.The entire range of elevation variation is sliced into several layers.Then, the RPCs are calculated by a least squares adjustment with these virtual control points (Tao and Hu, 2001).Grodecki and Dial (2003) proposed a block adjustment method for the VHR satellite imagery where the geo-referencing accuracy of the RFM is improved by use of a few numbers of control points.This RPC block adjustment method was implemented in the FORSAT software.It can run for stereo, triplet and block image configurations.

DSM generation
The automated DSM generation was performed using a modified version of the multiple primitive multi-image (MPM) matching method introduced by Zhang and Gruen (2004), Zhang (2005) and Zhang and Gruen (2006).In order to achieve successful and reliable results, the method matches a dense pattern of features with an appropriate matching strategy, making use of all available and explicit knowledge, concerning sensor model, network structure, image content and geometrical constraints such as the epipolar geometry constraint.The approach combines area-based matching (ABM) and featurebased matching (FBM), matching parameter self-tuning, generation of more redundant matches and a coarse-to-fine hierarchical matching strategy (Zhang et al., 2006;Baltsavias et al., 2007).The workflow is given schematically in Figure 3.
After the pre-processing of the original images and production of the image pyramids, the area based and the feature based matching methods are run in parallel.Starting from the lowdensity features on the images with the low resolution, the matching procedure progressively approaches finally on the original resolution images.Since all the matching procedures are based on the concept of multi-image matching (two-fold and three-fold images) guided from the object space, any number of images could be processed simultaneously.The triangulated irregular network (TIN) from DSM is reconstructed from the matched features on each level of the pyramid using the Delaunay triangulation method, which in turn is used in the subsequent pyramid level for the approximations and adaptively computation of the matching parameters.Finally, least squares matching methods are used to achieve more precise matches for all the features and to identify some false matches.
Figure 3. Automated DSM generation in the FORSAT software.
The entire system consists of three mutually connected subsystems: the image pre-processing module, the multiple primitive multi-image matching (MPIM) module and the refined matching module.The image pre-processing module is used to reduce the effects of the radiometric problems and optimise the images for subsequent feature extraction and image matching procedure.A combined matching process (point matching, edge matching and relational matching processes) goes through all the image pyramid levels in the MPIM module and generates good enough approximations for the refined matching module.In the final refined matching module, the least squares matching methods are performed only on the original resolution images to sub-pixel accuracy for all matched features obtained in the MPIM module (Zhang, 2005).

3D comparison and analysis
The co-registration is crucially needed wherever spatially related data sets, described as surfaces, has to be transformed to each other.Examples can be found in medicine, computer graphics, animation, cartography, virtual reality, industrial inspection and quality control, change detection, spatial data fusion, cultural heritage, photogrammetry, etc.Since DSMs represent the object surface, the problem can be defined as a surface co-registration problem.There have been some studies on the co-registration of DSMs for control information and for change detection tasks.This work is known as DEM matching (Ebner et al., 1988;Rosenholm and Torlegard, 1988;Mitchell and Chadwick, 1999).This method basically estimates the 3D similarity transformation parameters between two DEM patches, minimising the sum of the squares of the elevation differences (1D along the z-axes).The 1D elevation differences may not truly represent the surface-to-surface distance where terrain is complex with undulations.
For quality evaluation of DSMs, often a reference DSM is interpolated in the DSM to be checked.This approach is suboptimal (Gruen et al., 2004), since: • at surface discontinuities surface modelling errors may lead to large height differences although the measurements are correct (Figure 4a), and if the reference frames of the two DSMs differ (e.g.shifts and tilts), then again large differences occur, especially at discontinuities although the heights may be correct (Figure 4b).
These shortcomings can be overcome by employing the approach where the shortest 3D (Euclidean) distance between each reference point and the produced DSM is used (Gruen and Akca, 2005;Akca, 2010;Akca et al., 2010).Although the registration of 3D point clouds and DSMs is a very actively working area in many disciplines, we notice that a contribution that responds favourably to the following aspects is needed: (1) co-registration capability with higher order spatial transformation models, (2) co-registration and comparison of full 3D surfaces (as opposed to 2.5D), (3) a rigorous mathematical formulation for high accuracy demands, (4) a flexible model for further algorithmic extensions, (5) mechanisms and statistical tools for internal quality control, and (6) capability of matching of data sets in different quality and resolution.As a consequence, a fully satisfying general solution was developed and implemented in the FORSAT software.
We opted for the least squares 3D surface matching (LS3D) method (Gruen and Akca, 2005;Akca, 2007;Akca, 2010).The LS3D method is a rigorous algorithm for the matching of overlapping 3D surfaces and/or point clouds.It estimates the transformation parameters of one or more fully 3D surfaces with respect to a template surface, using the generalised Gauss-Markov model, minimising the sum of the squares of the Euclidean distances between the surfaces.This formulation gives the opportunity to match arbitrarily oriented 3D surfaces, without using explicit tie points.It is a powerful method whose accuracy and precision potential is directly dependent on the quality of the input data.Details of the procedure can be found in Akca and Gruen (2005;2007).Several applications ranging from 3D modelling (Akca, 2012) to geomorphology (Akca and Seybold, 2016) are available.The 3D co-registration and comparison module of the FORSAT software is a specialised implementation of the LS3D method.

User interface
FORSAT Software suite is not a monolithic system which means every core module works independently but related to each other at the same time.FORSAT Software suite uses a graphical user interface (GUI) which eases users to work with the software (Figure 5).
Figure 5. GUI of the FORSAT software.

Cyprus test site
The climate of Cyprus is characterised by mild and rain-laden winters as well as dry and hot summers.Due to its climate conditions it is predestined for the breakout of forest fires.The high risk for forest ecosystems is also boosted through human interventions.The forest fire of Saittas, raged in 2007, was defined as a demonstration site for the generation of the pilot application.A reason for this decision was the disastrous impact of the fire to the local vegetation which is still visible despite the long period of years between the event and image acquisition.The area of interest (AOI) delimited by a red frame covers an area of about 45 km² and includes the city of Pelentri in the Limassol district (Figure 6).There was an Ikonos stereo pair from October 2001 available which was used for the calculation of the DSM before the forest fire.For the actual date no VHR stereo coverage was available so a new Pléiades stereo acquisition was ordered for the specific area.The Pléiades coverage was acquired in July 2014.
Figure 6.The AOI is located in the centre of Cyprus.
Based on the archived Ikonos stereo pair and the new Pléiades stereo pair, a historic DSM (Figure 7a) and an actual DSM (Figure 7b) were generated using the FORSAT software.The resulting DSMs were generated with a resolution of 2.0 meters.On the basis of the generated DSMs, the change analysis was performed with "3D Comparison and Analysis" module of the FORSAT software.The effect of the forest fire of year 2007 can be seen in the west of the change map (see the large blue ellipse in Figure 8 and Figure 9a).Furthermore, there are smaller burned areas near Pelentri in the north-east of the AOI (see the blue circle in the north-east of Figure 8 and Figure 9b).The decreased areas in the south-east of the AOI are man-made changes (Figure 9c).  Figure 10 shows that more than half (54.6 %) of the broadleaved forest area, 47.5 % of conifer forest and more than a third (36.0 %) of the undefined (forest) area were affected by the decrease in height.We noticed that the forest has still not recovered from the fire in 2007.Figure 11 visualizes the change in volume in fire-affected vegetation cover based on height differences in DSMs between 2001 and 2014.It illustrates that there is a large decrease (approx.13 million m 3 ) of conifer forest.However, this class has also the highest volume increase (365'892 m 3 ), but it must be noted that conifer forest cover almost 70% of the whole fire-affected area.Other land cover classes with a volume decrease of nearly 200'000 m 3 and more are undefined and broadleaved forest, tree cultivations as well as bush and shrubs (maqui).Besides conifer forest only tree cultivations have a volume increase of more than 100'000 m 3 .An interesting fact is that the land cover class annual cultivations decreases the area, but increases the volume in total.
Figure 10.Change in fire-affected vegetation cover.
Figure 11.Change in volume in fire-affected vegetation cover.

Austria test site
Monitoring of the forest areas in terms of quantifying changes over time is a crucial objective in Austria.High winds cannot only cause direct damages, but also lead to threats concerning the break out and the spreading of forest fires.Fires are great dangers for forests and subsequently for humans as the forests have a protection function for inhabited areas regarding other natural hazards like avalanches, landslides or rock falls.Due to a recent forest fire event near Absam in Tyrol (Austria), which started on 20th of March 2014 and raged about two days, the corresponding burned area was selected as a test site.The area of interest (AOI) covers about 62 km² and includes the burned area from the aforementioned forest fire as well as the eastern suburbs of Innsbruck.Hence, the test site includes a variety of landscapes like mountains, dense urbanised areas surrounded by intensive agriculture as well as rural areas (Figure 12).The occasional clear-cuts and storm damages in the north and the vegetation growth in the south are dominating patterns.The enlarged view of the forest fire, delineated with the blue rectangle, is given in Figure 15a.The dark green region in the centre of the false colour composite image (Figure 15b) is the fire affected area, which is shown in orange colour in its change DSM counterpart (Figure 15a).The fire burned area was dominated by shrubs (proportionally shorter than trees) according to the reports of the Forest Department.This is the reason why the burned areas are not as dominant as the deforestation areas due to clear-cuts or storms such as the red plots in the upper left and lower right parts of Figure 15a.There are three classes in the test area: timber forest, shrubs and clear-cuts.The forest map showing the boundaries of these three classes is obtained in vector format and overlaid with the change results given in Figure 14.Thus, changes for each individual class are computed (Figure 16).In spite of the forest fire, changes of cover of shrubs are minor, and 99.1% area of the shrubs coverage remained unchanged.17.5% area of the timber forest class increases in height, which is clearly visible as growth of the vegetation in the southern part of the test area (below green areas in Figure 14).The results of the volumetric comparison (Figure 17) show that the trend is in the same direction with the area comparison, but the magnitude is overwhelmingly nonlinear.The timber forest class gained 17.7 million m 3 volume of stock whereas it lost 7.2 million m 3 due to deforestation.A detailed analysis of the burnt area, visualised in Figure 15, shows a loss of volume caused by the fire in March 2014 of 41'800 m 3 shrubs and 40'100 m 3 timber forest.

CONCLUSIONS
We developed a software suite, so called FORSAT which is capable of providing spatial information for rapid monitoring of forest areas.The basic input data is the VHR satellite imagery.
The software consists of mutually linked modules which are pre-processing, geo-referencing, DSM generation, and 3D comparison and analysis.The pilot application study demonstrates the capability of the FORSAT software especially for change analysis of the forest burnt areas.The results of FORSAT provide single source, flexible forest information solutions with a very competitive price vs. quality ratio, allowing for new market entry in the forest sector.

Figure 1 .
Figure 1.Main components of the FORSAT software.

Figure 12 .
Figure 12.The test site is located near Innsbruck in the Austrian province of Tyrol.The pre-fire data, which dates back to 2006/2007, is a 1-meter resolution DSM derived from an airborne LiDAR flight (Figure13a).A set of Pléiades triplet images was acquired in June 2014 and used to represent the post-fire situation.The Pléiades triplet was processed using the FORSAT software and a 1-meter resampled DSM was generated (Figure13b).
Figure 14.3D Change analysis.The red colour indicates deforestation and the green colour for growth of the vegetation.
Figure 15.(a) Enlarged view of the fire area which is shown as a blue rectangle in Figure 14, (b) Pléiades (in June 2014) false colour composite of the same extend.

Figure 16 .
Figure 16.Change of each class coverage.

Figure 17 .
Figure 17.Change in volume of each class.