POST-DISASTER DAMAGE ASSESSMENT THROUGH COHERENT CHANGE DETECTION ON SAR IMAGERY

Damage assessment is a fundamental step to support emergency response and recovery activities in a post-earthquake scenario. In recent years, UAVs and satellite optical imagery was applied to assess major structural damages before technicians could reach the areas affected by the earthquake. However, bad weather conditions may harm the quality of these optical assessments, thus limiting the practical applicability of these techniques. In this paper, the application of Synthetic Aperture Radar (SAR) imagery is investigated and a novel approach to SAR-based damage assessment is presented. Coherent Change Detection (CCD) algorithms on multiple interferometrically pre-processed SAR images of the area affected by the seismic event are exploited to automatically detect potential damages to buildings and other physical structures. As a case study, the 2016 Central Italy earthquake involving the cities of Amatrice and Accumoli was selected. The main contribution of the research outlined above is the integration of a complex process, requiring the coordination of a variety of methods and tools, into a unitary framework, which allows end-to-end application of the approach from SAR data pre-processing to result visualization in a Geographic Information System (GIS). A prototype of this pipeline was implemented, and the outcomes of this methodology were validated through an extended comparison with traditional damage assessment maps, created through photo-interpretation of high resolution aerial imagery. The results indicate that the proposed methodology is able to perform damage detection with a good level of accuracy, as most of the detected points of change are concentrated around highly damaged buildings.


INTRODUCTION
Damage assessment is a fundamental step to support emergency response and recovery activities in a post-earthquake scenario.This process is typically supported by damage surveys carried out by specialists who visit the earthquake area, inspect buildings and other physical structures, and classify the damages occurred.In recent years, UAVs and satellite optical imagery, often made available with open licenses to the public (Harris and Baumann, 2015), was applied to assess major structural damages in a faster way, before than technicians could reach the areas affected by the earthquake.However, the detection of changes by photointerpretation is not straightforward, while bad weather conditions may harm the quality of these optical assessments, thus limiting the practical applicability of these techniques.
In this paper, the application of Synthetic Aperture Radar (SAR) imagery is investigated and a novel approach to SAR-based damage assessment is presented.While damage assessment techniques based on optical images typically require a manual interpretation activity (Milisavljevic et al., 2015), SAR-based techniques are more sensible to very small variations and may be used together with automatic processing tools.In this context, Coherent Change Detection (CCD) algorithms on multiple interferometrically pre-processed SAR images of the area 1 * Corresponding author affected by the seismic event can be exploited to automatically detect potential damages to buildings and other physical structures.
As a case study, the 2016 Central Italy earthquake involving the cities of Amatrice and Accumoli was selected.Relevant SAR data were captured by Sentinel-1A, a European radar imaging satellite launched as part of the European Union's Copernicus program in partnership with the European Space Agency (ESA).
The main contribution of the research outlined above is the integration of a complex process, requiring the coordination of a variety of methods and tools, into a unitary framework, which allows end-to-end application of the approach from SAR data collection to result visualization.A prototype of this pipeline was implemented, and the outcomes of this methodology were validated through an extended comparison with traditional damage assessment maps, created through photo-interpretation of high resolution aerial imagery.The paper is organized as follows: section 2 contains presents the research methodology that was adopted; section 3 describes how the processing pipeline was implemented in terms of software components and their interactions; section 4 presents the 2016 Central Italy earthquake case study; finally, section 5 includes conclusions and future directions for the work.

METHODOLOGY
The procedure presented in this paper and implemented in the software pipeline is mainly based on the application of Coherent Change Detection (CCD) techniques (Wright et al., 2005) on a multi-temporal stack of SAR images pre-processed by the means of interferometry.When this approach is applied on a pair of interferometric SAR images, it uses the magnitude of their complex cross correlation to quantify changes in the observed amplitude and phase of the image pixels (Bamler and Hartl, 1998).An alternative to CCD is Incoherent Change Detection (ICD) that detects changes dependent on the mean backscatter power of the scene.The intensity of the backscatter depends on the structural and dielectric properties of the objects in the scene which proves useful at detecting changes on vegetation (Oliver and Quegan, 1998), which is not directly useful for damage assessment or urban monitoring but may find use in separating the scenes of low and high coherence.Though both CCD and ICD measure different properties of the scene, they provide complementary characterizations of scene changes (Preiss and Stacy, 2006).Therefore, a technique that exploits a two-stage change detection (Cha et al., 2015) results in a better robustness at the cost of greatly increasing complexity of the algorithm.
A CCD algorithm, such as the one implemented, can detect very subtle scene changes, making it prone to false alarms and extremely sensitive to the relative spatial geometry of the scene.For example, sometimes CCD detects undesired changes such as motions of leaves and changes in the vegetation.For this reason, images related to multiple SAR passes on the same area might be combined to cancel out disturbances and limit the number of false alarms (Barber and Kogon, 2012) in order to detect only major scene changes, such as the ones due to natural and man-made disasters.Supported by the fact that by multiple pass SAR acquisitions we can recognize the regions of low coherence, such as that of vegetation and landscape, in our implementation the algorithm is limited to detecting changes on buildings.Regarding the changes in the spatial geometry of the acquisition, in our implemented pipeline we rely on the inSAR processing provided by Sentinel Application Platform SNAP to rectify and align each acquisition of the passing Sentinel-1 satellite.
To investigate a change for an area of interest P, a stack of Ni time-sorted repeated-pass interferometric SAR acquisitions, that may differ in the geometry position of their sensor, is needed.Within that stack of Ni images, we differentiate between the group before and after the timestamp of the occurring event we want to analyze.Furthermore, for each acquisition of the stack Ni centered around the area of interest P, there are Ns samples, containing both amplitude and phase, forming the matrix X [Ni,Ns].The procedure relies onto the calculation of the probability density functions for change and no-change upon the information of the columns of the matrix X under the assumption that after the event Nc no further changes occur.To decide upon the occurrence of a change we use the log likelihood ratio Λ(Nc) for coherence matrices C0 and CNc, which are normalized covariances of the cases of no change and change respectively; C0 has a diagonal of ones and constant γ everywhere else, CNc should continue the behavior of C 0 except that the fields representing cross correlation for the acquisitions before and after event Nc should be zero resulting in Figure 1.= estimate of the data covariance Following the (Fan et al., 2001) Generalized Likelihood Ratio Test (GLRT) for evaluating the significance of changes, they can be valued as detectable by the local minima in the L to which the GLRT is applied.Looking at Figure 1 it can be noticed that the entries are averaged together, separately for those pre-change and those post change resulting in: Nc × (Nc -1)/2, Nc × Nq , and Nq × (Nq -1)/2 where Nq is the of acquisitions before event Nc equating to Nq = Ni -Nc.This helps the robustness of the estimation, by reducing bias as proven by (Novak, 2005).In the end, the resulting GLRT percentile is the estimated probability which inevitably helps to assign the geographical location of the destructive change detected by the algorithm.

An integrated pipeline
One of the main goals of this research was to wrap the process into an integrated pipeline, that may be used to automatically preprocess SAR data for a specific disaster, perform change detection, and visualize the results.This is made possible by means of QGIS, one of the most popular Free and Open Source Software (FOSS) cross-platform GIS.Therefore, for this work the Sentinel Application Platform (SNAP) is used and its integration with QGIS is investigated in order to facilitate the replication of the analysis as well as to promote its usage among a wider and less specialized community.

QGIS plugin
QGIS allows users to create custom plugins as well as complex data processing pipelines using Python, which can integrate both external libraries and custom algorithms with native geospatial data manipulation and analysis tools.Moreover, being QGIS a FOSS, developed source code can be shared among the community of users thus enabling further revisions and improvements of the presented procedure.Nevertheless, QGIS does not currently provide extensive support for SAR data processing.For this reason, a QGIS plugin was implemented in order to hide the complexity of the processing activities included in the pipeline behind an easy-to-use component fully integrated within the QGIS user interface.
The plugin was implemented leveraging the power of the QGIS Application Programming Interface (API) which allows developers to access functionalities of the GIS software and to integrate their modules directly inside the toolbars of QGIS.The plugin was implemented in Python because, when compared to other programming languages which may be used to develop QGIS plugins (C++), Python is typically easier to understand, but also quicker to write and maintain.The development process of the plugin was also supported by the usage of the QGIS Plugin Builder, a template which may be used to speed up the plugin creation process, and by Qt Designer, a tool for designing and building Graphical User Interfaces.The plugin was implemented with a focus on the usage of FOSS solutions, and on crossplatform compatibility: the software components were tested on machines running Windows, macOS and Linux.
The two main modules which were integrated in the QGIS plugin are: 1. SAR data pre-processing module 2. Coherent Change Detection module Further details about each software component and its implementation are provided in following paragraphs.

SAR data pre-processing module
SAR images captured in Interferometric Wide (IW) swath mode (data with a 250 swath at 5 m by 20 m spatial resolution) need to be pre-processed by the means of terrain correction, coregistration and interferometric analysis before running Coherent Change Detection algorithms on them.
In more detail, the pre-processing activities are carried out by using the SNAP platform and specifically the Sentinel-1 Toolbox (S1TBX), a collection of tools which may be used to support processing activities such as data conversion, calibration, noise filtering, coregistration, orthorectification, mosaicking and interferometry on SAR data.Using these tools, for each SAR product, the sub-swath containing the area of interest, is coregistered with respect to a master, obtaining a pixel-to-pixel mapping.The two images, the master and the slave, are then debursted.At this stage, the area of interest can be selected.The two images need to be phase calibrated, prior to perform the analysis, to remove most of the systematic phase contributions.The topographic phase is estimated by adding the contribution of flat terrain and that of the DEM: this phase is removed from each of the slave images.All these steps are performed both on the VV and the VH polarization, giving two phase compensated stacks of images.
From a technical standpoint, the processing activity described above was implemented as a Python script which calls externally SNAP to execute functional graphs created by means of SNAP Graph Builder.More specifically, as previously mentioned, SNAP tools were used to perform a set of operations: 1. Terrain Observation with Progressive Scans in azimuth coregistration (reading master and slave images, performing TOPSAR splitting, applying orbit files and back-geocoding); 2. TOPSAR debursting and interferometry; 3. results writing (with Topographic Phase Removal on a subset).
Finally, geocoding is performed in order to create a regularly gridded file in latitude/longitude that contains as x-y indexes azimuth and range of the master image in the SAR reference.This procedure is carried out by reading the results written by the previous workflow, by performing terrain correction and finally by writing the geoTIFF images that will be provided as input to the Coherent Change Detection module.

Coherent Change Detection module
Since SAR is a coherent imaging system, in principle coherent and incoherent change detection approaches may be adopted.Considering that coherent change detection allows to detect very subtle changes that incoherent techniques cannot identify, the former approach is selected and a Coherent Change Detection (CCD) algorithm based on normalized log-likelihood ratio test is used to perform change detection over the aforementioned multi-temporal stack of SAR images.A small sample of data is used to identify a proper change detection threshold, thus minimizing both the probability to miss a detection (false negatives) and the probability of detecting a nonexistent change (false positives).
The pre-processed stacks are given as inputs to the CCD algorithm, which returns latitude and longitude coordinates for each detected change in the area of interest.The output of CCD is then geocoded to identify changes at ground level.The result of this process is a CSV file containing, for each detected change, the latitude and longitude coordinates, the value of change statistics, and date.
From a technical standpoint, the aforementioned CCD algorithm was implemented as a collection of scripts written in GNU Octave, a scientific programming language largely compatible with Matlab which has the advantage of being part of the GNU Project, thus being a Free and Open Source Software released under the terms of a GNU General Public License.Moreover, the two additional packages, namely image and mapping, from Octave-Forge, a collection of packages providing extra functionality for GNU Octave, were used for their image processing, mapping and raster reading functions.Access to GNU Octave from Python code was enabled by the usage of three Python packages: 1. oct2py, a package which allows developers to seamlessly call Octave functions and M-files within a Python script; 2. NumPy, a package for scientific computing; 3. SciPy, a package which provides many efficient numerical routines.

Data and experiment setup
As a case study, the 2016 Central Italy earthquake involving the cities of Amatrice and Accumoli was selected.Relevant SAR data were captured by Sentinel-1A, a European polar-orbiting, all-weather, day-and-night radar imaging satellite launched as part of the European Union's Copernicus program in partnership with the European Space Agency (ESA), capable of taking copolarized V/V and cross-polarized V/H images simultaneously.
A dataset of 32 images including the area of interest and collected between February and October 2016 was selected and downloaded from the ESA Copernicus Open Access Hub website.More specifically, such dataset is made by a stack of 15 descending-pass images with relative orbit number 22 and a stack of 17 ascending-pass images with relative orbit number 117.
In Tables 1 and 2 the sensing dates of the images of the two stacks are listed.Since the earthquake struck on August 24th, 2016, the descending pass stack has 10 pre-disaster images (Nc=10) and 5 post-disaster images, whereas the ascending pass stack has 14 pre-disaster images (Nc=14) and 3 post-disaster images.

Image number
Moreover, it is important to underline that the two stacks have been independently processed, but both refer to the same area of interest covering the cities of Amatrice and Accumoli.In particular, the pre-processing module of the pipeline was run by giving as input a rectangular bounding box, whose latitude and longitude in decimal degrees respectively range from 42.60 N to 42.71 N and from 13.24 E to 13.33 E, as shown in Figure 4.

Experimental Results and Validation
The integrated pipeline described in the previous section was applied to the two aforementioned stacks of SAR images.Therefore, two different CSV files, containing for each detected point of change the latitude and longitude coordinates, were obtained and then graphically visualized in QGIS as an additional layer over the map imported from OpenStreetMap, as shown in Figures 5 and 6.Furthermore, the effectiveness of the obtained outcomes in terms of damage identification was assessed by means of an extended comparison with data taken from independent and traditional damage assessment campaign in the considered area, which realized (10 cm) aerial imagery (Copernicus Emergency Management Service, 2016).More specifically, this kind of damage assessment allows to determine a damage grade scale for buildings and roads with 5 levels of damage (Destroyed, Highly Damaged, Moderately Damaged, Negligible to slight damage and Not Affected).Such preliminary assessment of the performance of our methodology is presented in Figure 7, which has been obtained by overlapping in QGIS the layer of the points of change (red dots) and the Copernicus EMS buildings damage grading over aerial imagery of the urban area of Amatrice.By visually inspecting the map, it is possible to notice that the points of change are mostly concentrated around the buildings marked as Destroyed (in red), Highly Damaged and Moderately Damaged (in orange), demonstrating that the presented pipeline is very sensitive in detecting changes occurred to the buildings which have been consistently affected the earthquake.

CONCLUSION
Since environmental hazards happen without any warning in an unexpected area, it is hard to aid such area within a strict time frame using CCD algorithm if all steps to its utilization are performed manually.In that aspect, the project succeeds in creating a seamless early damage evaluation pipeline which may be used with little effort, except for the retrieval of the imagery from the Copernicus platform, which implies a manual selection of the proper pass acquisitions related to the disaster area and might require a significant amount of time due to the size of the images to retrieve.In the end, this project is giving a chance to provide an early damage evaluation to aid some more detailed evaluation later in the course of disaster management.
As suggested from the results from the validation, the CCD pipeline takes to a reliable semi-automated change detection given an adequate amount of post-event acquisitions.However, as the detection algorithm relies on at least two acquisitions to eliminate low coherence such as vegetation and landscape, it is crucial that those acquisitions occur immediately after the analyzed event.For the Sentinel-1 imagery, that is a major drawback as the repeated passes of the satellite cannot be finely tuned to perform multiple passes over an area in a short time.If the automated detection provided by the pipeline is crucial for a very early damage evaluation, it may fail under a stricter deadline.In essence, since SAR technology is superior in bad conditions compared to optical imagery, the information provided is considered extremely useful if such conditions prolong trough an extended period.Because of the aforementioned behavior of our implementation, the pipeline may become essential in an event of hazard happens during a prolonged bad weather conditions.In case of normal weather conditions, it is still a great guideline to start the damage mapping and offer some basis of analysis.
The next step of the development of the project is to test the pipeline against a larger variety of scenarios and locations, and consider the possibility of its use in environments that may have structures that differ in geometry than the ones in Amatrice.Furthermore, aware of the shortcomings of the project, we plan on a future work that may involve crowdsourced data to provide a very fast detection validation of results.In more details, the microblogging platform of Twitter provides a source of information that is really fast in responding to environmental disasters.Using such data may inevitably change the performance of the evaluation to rely on fewer acquisition.In addition, a more ambitious vision for such project is to allow for full automation even for image retrieval from the Copernicus platform and an improved per-scenario performance, that may better help in the case of evaluating a wider range of disaster scenarios.

Figure 1
Figure 1.coherency matrixIn the end we formulate the Log-Likelihood Ratio as in equation (1), where the absolute value of the covariance matrix is assumed here to remove the random phase field superposed at each image.

Figure 2 .
Figure 2. Component diagram of the processing pipeline The output of the Coherent Change Detection module, a CSV file containing, for each detected change, the WGS84 coordinates, the value of change statistics, and the date, is automatically imported in QGIS as a Delimited Text Layer using the QgsVectorLayer class of the QGIS API.result of the pipeline is a graphical visualization of the change detection map, as it will be observed in Section 4.

Figure 4 .
Figure 4. Area of interest and bounding box, basemap by OpenStreetMap

Figure 5 .
Figure 5. CCD results visualization for the ascending-pass stack (area of Amatrice), basemap by OpenStreetMap

Figure 7 .
Figure 7. CCD results visualization for the descending-pass stack over aerial imagery and grading map of the urban area of Amatrice produced by Copernicus Emergency Management Service and e-GEOS (Situation as of 25/08/2016).

Table 2 .
Stack of ascending-pass images