COHERENT CHANGE DETECTION FOR REPEATED-PASS INTERFEROMETRIC SAR IMAGES : AN APPLICATION TO EARTHQUAKE DAMAGE ASSESSMENT ON BUILDINGS

During disaster response, the availability of relevant information, delivered in a proper format enabling its use among the different actors involved in response efforts, is key to lessen the impact of the disaster itself. Focusing on the contribution of geospatial information, meaningful advances have been achieved through the adoption of satellite earth observations within emergency management practices. Among these technologies, the Synthetic Aperture Radar (SAR) imaging has been extensively employed for large-scale applications such as flood areas delineation and terrain deformation analysis after earthquakes. However, the emerging availability of higher spatial and temporal resolution data has uncovered the potential contribution of SAR to applications at a finer scale. This paper proposes an approach to enable pixel-wise earthquake damage assessments based on Coherent Change Detection methods applied to a stack of repeated-pass interferometric SAR images. A preliminary performance assessment of the procedure is provided by processing Sentinel-1 data stack related to the 2016 central Italy earthquake for the towns of Ametrine and Accumoli. Damage assessment maps from photo-interpretation of high-resolution airborne imagery, produced in the framework of Copernicus EMS (Emergency Management Service European Commission) and cross-checked with field survey, is used as ground truth for the performance assessment. Results show the ability of the proposed approach to automatically identify changes at an almost individual building level, thus enabling the possibility to empower traditional damage assessment procedures from optical imagery with the centimetric change detection sensitivity characterizing SAR. The possibility of disseminating outputs in a GIS-like format represents an asset for an effective and cross-cutting information sharing among decision makers and analysts.


INTRODUCTION
According to the general definition, disasters always imply a serious disruption of the functioning of a community involving widespread human, material, economic, and/or environmental losses (UNISDR, 2009).Focusing on material losses, timely availability and accessibility of spatially contextualized information, enabling disaster responders to identify affected areas as well as to assess the damage grade on infrastructures, is critical for the implementation of any effective Disaster Management (DM) practices.In the recent past, geospatial data and technologies have been extensively adopted to fulfil this informational requirement.Indeed, Geographic Information (GI) and Geographic Information Systems (GIS) are nowadays recognized as keys within all the phases of the DM cycle (Thomas et al., 2007).Among the available technologies, spaceborne Earth Observation (EO) has terrifically changed the perspective for the operational disaster response by reducing both effort and uncertainties connected to the traditional ground-based information collection (Voigt et al., 2016).
During the last two decades, national and international organizations have set up extensive investment plans targeting EO technologies.This has resulted in a growing availability of sensors and platforms delivering high-resolution and up-to-date EO data, which is more often released with open licenses to the public (Harris and Baumann, 2015).One of the most relevant examples of the above is represented by the Copernicus Earth Observation Programme of the European Union (http://copernicus.eu).
In the context of disasters, this emerging EO information has paved the way to a number of applications addressing heterogenous topics, including floods (e.g.Kwak et al., 2016), tsunamis (e.g.Yamada, 2015), and earthquakes (e.g.Koyama et al., 2016).From the operational point of view, optical and Synthetic Aperture Radar (SAR) satellite images are nowadays the most popular data sources providing emergency services and decision makers with essential information during disaster response (Joyce et al., 2009).Optical images are the most common EO instrument having applications in different areas such as agriculture, land-cover mapping, damage assessment and urban planning.Optical images, however, are limited to cloudfree conditions and daytime operation which, in the case of disasters, might represent a serious limitation to the use of this data for response and recovery purposes.SAR provides instead day-night and weather independent images.Moreover, the precise repetition of the acquisition geometry combined with the coherent imaging makes of SAR a valuable sensor for automatic detection of Earth surface changes -either at a small or at a large-scale -related to disasters.Being based on interferometry, change detection with SAR is sensible to centimetric variations and generally less biased than any techniques based on optical images classification, which usually requires manual interpretation (Milisavljevic et al., 2015).Multi-temporal SAR observations are currently employed to perform change detection at a large-scale, such as terrain deformation analysis after earthquakes and flood areas delineation.However, for particular tasks such as the damage assessment on buildings and infrastructures after a disaster, a higher spatial resolution is required.Most of the operative spaceborne SAR platforms have reached a spatial resolution for data acquisition which potentially allows change detection at a building scale (Plank, 2014).This is followed by a decreased revisit time due to the presence of larger satellite constellations.Nevertheless, both in the literature as much as in practical applications, there is still a lack of examples considering change detection techniques to cope with the increasing spatial and temporal resolution of the available SAR observations.
In this paper, we propose a procedure based on Coherent Change Detection (CCD) techniques (Wright et al., 2005) optimized for the application to multi-temporal stacks of interferometric SAR images, allowing a more robust pixel-wise change detection than the traditional two-images based approach.This asset is exploited in the assessment of damages to single buildings and infrastructures after disasters such as earthquakes.The procedure is tested by taking advantage of the Sentinel-1 (S1) open SAR data provided by the Copernicus Programme.The paper is organized as follows: Section 2 contains a summary of the methodology adopted.Section 3 introduces the potential mutual benefit for the combination of the proposed methodology with the S1 SAR data.In Section 4, the experimental software pipeline implementing the methodology is outlined.Section 5 describes preliminary outputs obtained for the 2016 central Italy earthquake, which was selected here as the case study.Conclusions and future directions for the work are finally discussed in Section 6.

METHODOLOGY OUTLINE
We approach here the identification of a possible change occurring for the target of interest P between the images Nc and Nc+1 from a stack of Ni repeated-pass SAR images sorted by acquisition date.The reference geometry for the CCD estimates is shown in Figure 1.We define X as a [Ni , Ns] matrix with Ni repeated-pass observations, each made by Ns samples, taken from a window centred on the target of interest P. The detection of changes is performed considering both amplitude and phase measurements by evaluating the ratio between the probability of no-change and that of a change occurred after epoch Nc.This is the ratio of the two zero-mean normal Probability Density Functions (PDFs) of all the columns of X under the assumption of no-change after Nc and the assumption of change after Nc.This is expressed in (1) as the Likelihood Ratio Λ(Nc).In order to compute the Λ(Nc), we assume C0 and CNc to be the coherence matrices (normalized covariance) for the two cases of no-change and change after Nc images.C0 is constant everywhere and 1 on the diagonal.CNc should be instead defined according to some coherence change model.Here we assume that a target has constant coherence γ in all the epochs pre-event and the same for all the epochs postchange, whereas we assume that there is no correlation between and after the change (Figure 2).The likelihood of a change after Nc is: Significance for the detected changes is inferred by means of the Generalized Likelihood Ratio Test (GLRT) (Fan et al., 2001).
Changes are detectable by the local minima in the L to which GLRT is applied.Notice that the coherence pre and post-change is positively weighted, thus favouring the detectability, whereas the one prepost change has a negative weight, that would reduce the detectabilityunless coherence is null.The block structure of the coherence matrix shown in Figure 2 implies that all the entries pre-change are averaged together, and the same happens for the post-change ones.It is like coherence is estimated by averaging Nc × (Nc -1)/2, Nc × Nq , and Nq × (Nq -1)/2 where Nq = Ni -Nc is the number of data pre-event.Such averaging improves the quality of the estimation, by reducing the bias, and leads to a much more robust approach than the two-images one, that is the desired goal.A proof of this latter, together with a comprehensive theoretical discussion and the mathematical formulation of the methodology, is provided in (Novak, 2005).
Outcomes of the CCD are assigned back to their targets of reference, together with the indication of the change epoch.This measure, i.e. the GLRT percentile, is an estimation of the probably of change which can be later exploited for damage assessment by knowing the location of each target in the geographical space

IMPLEMENTATION BY SENTINEL-1 SAR
The detection method theoretically proposed by (Novak, 2005) in the multi-pass case has been implemented for the twopass high resolution airborne (Preiss et al., 2006;Wahl et al., 2016), and extended to the repeat-pass spaceborne SAR by means of a proper phase calibration (Monti Guarnieri et al., 2018).This is of particular interest for the Copernicus S1 SAR.
The constellation is made by two C-band SAR, while other two will be launched after 2020, that systematically acquires data over land masses, worldwide, and they made available nearly in real time and free on Copernicus scientific data hub.The orbit repeat cycle is twelve days for each sensor, that gives and interferometric repeat of six days when two satellites are combined.Then, a revisit shorter than three days can be achieved can combining interferometric series acquired by different view angles and complementary ascending and descending geometries.This makes of S1 a precious system to precise and rapid mapping of changes.

EXPERIMENTAL SOFTWARE PIPELINE
The methodology outlined in Section 2 is wrapped into an experimental software pipeline to allow its application and testing (Figure 3).The pipeline combines existing SAR data preprocessing algorithms, provided by the Sentinel Application Platform (SNAP: http://step.esa.int/main/toolboxes/snap),together with custom scripts which perform the CCD on the preprocessed SAR stacks by converting the outputs into a geospatial layer to allow visualization and post-processing operations within a GIS environment.
The procedure requires a set of repeated-pass SAR images covering the area impacted by the disaster event and acquired both before and after the time when the disaster struck.The number of pre and post-event images is not a priori defined.This choice influences the averaging procedure for coherence estimation as mentioned in the previous section.Theoretically, the larger is the number of images both before and after the disaster event, the better the estimation of coherence changes.Nevertheless, the application to damage assessment after disasters might not consider long series of post-event images, due to the need of time-effective estimations which are expected from automatic procedures such as the one here presented.With this in mind, it is important to feed the procedure with a rich set of pre-event images and at least one post-event image, eventually performing additional CCD estimations at any time new postevent images become available, in order to improve the results reliability.
As a mere example, for the case study reported in the following section, we processed a set of ~ 20 SAR images with ~ ⅔ of preevent images.Two independent sets of images, both having the aforementioned characteristics, were collected to account for ascending and descending pass orbits.
The required pre-processing steps are shown in figure 3.These include: images coregistration, orbit correction, topographic phase compensation using a Digital Elevation Model (DEM), and phase calibration.These steps must be performed for creating suitable image stacks from which the coherence matrices are computed.SNAP provides users with separate functionalities to perform each step.The pre-processing procedure is automatized for multiple images by taking advantage of the SNAP Graph Builder, which enables to assemble and run chains of user-selected functionalities.The pre-processed stacks are then passed as input to custom scripts performing the CCD estimations by exploiting the three types of coherent combinations, i.e. the different polarizations: VV, VH, and joint VV+VH as well as different windows shapes and sizes.Outputs are geocoded and provided in a GIS-like format, thus enabling post-analysis and visualization within any GIS software.These consist of layers of points enriched with the computed GLRT percentile and the corresponding time of change.The whole pipeline is under development.This ongoing work aims at a future deployment of a stable software tool implementing the complete procedure.
Figure 3. Experimental software pipeline, adopted for implementing the proposed procedure

CASE STUDY: THE 2016 CENTRAL ITALY EARTHQUAKE
The proposed procedure is applied to the damages assessment on buildings after central Italy earthquake, that struck on August 24th, 2016.The analysis is performed for the town centres of Amatrice and Accumoli which were highly impacted by the disaster event.A preliminary assessment of the ability of the procedure in identifying damaged buildings is achieved by means The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3/W4, 2018 GeoInformation For Disaster Management (Gi4DM), 18-21 March 2018, Istanbul, Turkey of comparison with data from an independent damage assessment campaign in the area.

Data collection and processing
Only The two image stacks were processed separately using the procedure described in the previous section by testing different windows size (i.e.3,4,5 pixels) and shape, different polarization (i.e.VV, VH, VV+VH) and different options for phase calibration (space domain, or joint space and time domain).The detected changes were geocoded to a grid size of 10 x 10 m.Significant detected changes, expressed as GLRT percentiles, were associated with their reference targets.This information was then stored in a layer of points depicting the CCD targets (i.e.pixel centres) and containing their geographical coordinates, the GLRT percentile and the epochs of change (Figure 4).
The independent damage assessment data consists of a shapefile of Amatrice and Accumoli buildings enriched with damage grade (Figure 4).Damage grade scale provides five different damage classes namely: Completely Destroyed, Highly Damaged, Moderately Damaged, Negligible to slightly damaged, Not Affected.These grades were assessed by means of photointerpretation of 10 cm high-resolution airborne imagery, produced in the framework of Copernicus Emergency Management Service (EMS: http://emergency.copernicus.eu)and cross-checked with field survey (Sandu et al., 2017).Data is limited to the town centres and it was used as ground truth for the preliminary assessment of the damage detection performance of the proposed procedure.
Figure 4. Grading map of buildings and targets with significant post-event detected changes (descendent stack, Nc=13) for the town of Amatrice At a first visual inspection of the resulting maps, it is possible to observe higher concentration of targets, showing significant changes at Nc, overlapping urbanized areas.This proves the sensibility of the proposed methodology in detecting significant changes caused by the earthquake, which likely occurred in buildup areas.

Preliminary performance assessment
In order to assess damage detection performance of the procedure, we select for both layers of points computed from the descending and ascending stacks only those points for which detected changes happened at the Nc of the source stack.This is done for all the different CCD tests performed by varying both the windows size and the coherent combinations.Nine different configurations for CCD test were tested on both stacks.Details about the configurations are reported in Table 1.
Using the selected target subsets, we perform multiple spatial joins between selected points and the building polygons of the ground truth shapefile.The spatial join is based on the intersection between each building polygon and the 20 m buffer of each point.The buffer radius is selected assuming a positional uncertainty for the CCD targets equal to the maximum dimension of the original SAR image pixel.The number of significant targets results to be generally higher for the descending stack.The total number of buildings included in the ground truth shapefile is 627.The number of buildings included in the multiple spatial join are 560 for the descending stack and 513 for the ascending one.This operation has two purposes: a) to assign a GLRT percentile to each building involved in the spatial join by averaging the GLRT percentiles of all its spatially linked CCD targets, and b) to assign a damage class to each target by looking at their spatially linked buildings.By doing so, we assume an equal probability for each target of being related to a building and vice versa, providing that both elements are less than 20 m far from each other.Moreover, to a single target might be assigned more than one damage class.For this preliminary performance assessment, we accepted these approximations.

ID Test Code Window Size Polarizations
We select suitable CCD test configurations by testing if the lowest average GLRT percentile computed for the classes of Damage is higher than the highest average GLRT percentile computed for the classes of No-Damage, thus using this letter as a threshold for damaged buildings classification (Figure 5).The thresholds are compared with the average GLRT percentiles computed for each building by means of the multiple spatial join with CCD targets.The damage detection performance is assessed by means of binary confusion matrices, classifying damaged / not-damaged buildings according to the comparison between their average GLRT percentiles and the thresholds.The ground truth data is used for assessing the classification performances.A summary of the results is included in Table 2.  Preliminary results from the case study show generally higher performance for the descended stack, using windows size of 4 pixels and the joint coherent combination VV + VH.The maximum obtainable accuracy computed for these test settles around 60 % which can be considered satisfactory for this preliminary application.Nevertheless, the whole processing as well as performance assessment procedures require further improvements to remove possible bias generated by the approximations and artefacts we introduced in this first analysis.

CASE STUDY: THE 2016 CENTRAL ITALY EARTHQUAKE
In this paper, we presented the application of CCD techniques on repeated-pass interferometric SAR images to address automatic damage assessment on buildings after an earthquake.The theoretical methodology is outlined and the benefits than of the traditional two-images based approach are discussed.An experimental implementation of the procedure exploiting software solutions is also proposed.Results from the employment of the procedure for the central Italy earthquake in 2016 are reported together with a preliminary assessment of the achievable damage detection performances.
The proposed procedure as presented here requires additional advances to be viable for performing real automatic damage assessments.Besides the selection of the optimal configuration for the CCD tests, which might vary from case to case due to different settlement and landscape characteristics of the study area, the introduced damage / no-damage thresholds are here estimated through the use of existing damage information.A priori definition of such thresholds is required for an automatic damage estimation.The best setting of thresholds should be adapted on each single target by exploring its historical series prior to the event.Future investigations will address the two points above by means of extensive tests on additional case studies.These will aim at the parametrization of the procedure by looking for underlined links between real damages on buildings and the GLRT percentiles for the targets computed from the SAR image stacks.Nevertheless, preliminary results demonstrate the capability of the procedure to automatically locate damages on buildings with a satisfying precision, justifying the research interest in this topic.
The implementation of the procedure into a stable software tool is also planned.This will be addressed by an extensive use of Free and Open Source Software (FOSS) with a particular consideration of FOSS GIS solutions.The possibility of presenting results by means of simple maps while enabling postprocessing operations within a GIS environment is key to interoperability and future applications in the operative DM practices.

Figure 1 .
Figure 1.Left: geometry of a repeated-pass interferometric SAR acquisitions, showing the sensor position at different time, and the target under test.Right: model of the stack of Ni images, where a change occurred to the target P creates two subsets of Nc and Nq images

Figure 2 .
Figure 2. Model of the coherence and the weight matrix, for the case of a change occurring at Nc = 13 and coherence γ = 0.7 for all the targets before and after changes repeated-pass S1-A, were available in 2016, in the usual Interferometric Wide Swath (IWS) mode, double polarization (VV + VH) and ground resolution of 20 m (along track) x 5 m (across track).The data analysed included a) Ni = 18 descending pass images collected from February 11th up to October 20th, 2016 (Nc = 13), and b) Ni = 16 ascending pass images collected from March 12th up to September 26th, 2016 (Nc = 12).
After the definition of the spatial relationship between targets and buildings and in turn the identification of targets falling in each damage class, we compute the average GLRT percentile for each damage class for the different CCD tests.Different settings of the windows size and the coherent combinations produce different results.In order to identify tests providing adequate performances, the five damage classes are aggregate in two classes a) Damage, containing the former classes: Completely Destroyed, Highly Damaged, Moderately Damaged, and b) No-Damage, containing the former classes: Negligible to slightly damaged, Not Affected.

Figure 5 .
Figure 5. Average GLRT percentile by damage classes with damage / no-damage thresholds extracted for the selected CCD test configurations.Data refers to the descending pass stack.

Table 2 .
Summary results from the damage performance assessment of the procedure for the case study