Damage Detection from SAR Imagery: Application to the 2003 Algeria and 2007 Peru Earthquakes

This paper is focused on the improvement and further validation of a recently proposed approach for the joint use of radar satellite imagery of an area affected by a major disaster and ancillary data. The study was carried out at different sites on imagery of two different earthquakes occurred one in the Mediterranean coast of Algeria on May 21st, 2003, which severely affected the city of Boumerdes, and one in the Pacific Coast of Peru on August, 15th, 2007. The combination of different radar-extracted features results in very fuzzy classification of the damage patterns, far less detailed than what available using optical imagery. However, focused results using the above-mentioned ancillary data provide enough detail and precision to be comparable with them. In particular, quantized damage level at the block level is achieved at enough detail using ALOS/PALSAR data and thus validates the original idea.


INTRODUCTION
One of the most important issues in disaster damage detection is time, and map timeliness is as important as precision. As a consequence of that, every kind of data available is used in order to provide information to emergency operators [1]. To this aim, remotely sensed imagery can be instrumental, as well as geographical information system (GIS) layers, printed maps, and historical datasets. However, while the use of this kind of imagery has been constantly growing in the past few years, image interpretation tools, though fast and efficient even if not highly accurate, are still not used in many applications.
Damage assessment is actually a big challenge, being impossible to get the right data at the right time. For this reason, the organizations such as the International Charter on "Space and Major Disasters" are forced to consider a wide range of sources. The need of rapid damage pattern estimate requires information extraction about damages using quick and possibly efficient approaches suited to the decided spatial scale and the available data. So, the scale at which the analysis can be carried out is determined by the spatial resolution of the available data; for example, while SAR data can be useful to extract damage information only at a parcel level, it should be possible to recognize the single collapsed buildings from HR images.
The aim of this work is to understand the viability of radar satellite approaches to damage patterns by analyzing many different disasters all around the world, looking also at different scales of work. SAR data are becoming widely available with more and more fine spatial resolution, and thus with larger usability for urban area management. This improvement actually allows a better match between the growing requests for focused analysis in these areas (due to the concentration of population) with the enhanced availability of dataset. Still, approaches to disaster management in urban areas using SAR data are very limited, due to the problems in data interpretation and the lack of automated or semiautomated tools.

DAMAGE PATTERN ESTIMATE FROM SAR DATA
In recent technical literature, some works have already suggested that multitemporal SAR data may provide, at 2 International Journal of Navigation and Observation    a proper temporal and spatial scale, interesting information about disaster like earthquakes and floods. Most of these works concerning earthquakes need data coming from ground surveys to validate but also to initiate the process of information extraction. For this reason, these strategies are very useful in order to correlate damage patterns with ground displacements and soil properties [2,3], or to provide precise 3D changes of the earth crusts [4], but offer very poor results in terms of damage assessment and rapid damage mapping of an affected area. However, it is important to say that classification and change detection methodologies solely used cannot provide immediately usable results to the final user. Instead, these methods integrated with some kind of ancillary data allow obtaining more precise and understandable results. Moreover, damage analysis is almost uniquely required in urban areas and human settlements in general, where it is often easily feasible to collect layers of Geographic Information System (GIS) data.
In this work, we apply a technique recently proposed in [5] for the two test cases of the Bam (Iran) [6] and of the Golcük (Turkey) earthquakes. The first aim of this work is indeed to show that the proposed approach is valid in other situations and produces useful results fro damage assessment in different areas in the world. As a second objective, this paper reports the results of an investigations about the robustness of the approach to the lack of some of the features originally used in the cited papers.
The overall methodologyof the data analysis is proposed in Figure 1, even if, for sake of brevity, we do not recall here all the details of the algorithm. The procedure involves first of all the extraction of a suitable set of features from the original multitemporal dataset comprising pre-and postevent SAR imagery of the area under test. This feature set is then input to a multiband supervised classifier, whose output is a multiple class change detection map. Finally, a postclassification fusion step is performed at the end of the procedure involving the use of the above mentioned ancillary data.
G. Trianni and P. Gamba  The set of feature extracted from the multitemporal SAR images of the area under test depends on the trivial assumption that radar returns in damaged areas are quite different than in the original "undisturbed" configuration of the buildings. There are studies showing that urban areas show a remarkably strong coherence in complex return values and correlation in amplitude/intensity values during time. It is interesting therefore to use as hint to damages the change in the complex coherence and in the intensity correlation. In particular, intensity correlation has been tested in technical literature for the Hyogoken-nambu (Japan) and Bam (Iran) earthquakes. Following the paper where this was originally proposed [2], each (complex) SAR image of the available data sequence X n , n = 1, 2, . . ., is prefiltered with a Lee filter, and intensity correlation r i between the ith image and the previous one in the temporal sequence is computed according to the formula where x i = X i (recall that SAR data are complex values), the W (·) notation means that computation is done for each image element in a window W = N ×X around it, and finally the mean value x i is similarly computed in W. Along with intensity correlation, another valuable input feature is the difference between the logarithmic value of the mean prefiltered data intensity d i = 10 log 10 (x i ) −  10 log 10 (x i−1 ), and the feature set is completed by the original pre-and postevent pointwise intensities.
Of course, according to the dimension of the window W and thus to the geographic area of computation of the spatial features, multiple scales of analysis of the data can be enhanced. A convenient value for the window size, according to our past experience, depends on the ground spatial resolution of the data and the mean dimension of a meaningful block of buildings in the human settlement under test. In all the considered cases, a value of N between 15 and 21 is equally valuable, when the SAR images have spatial resolution in the 10 meters' range.
The second step of the approach, as shown in Figure 1, is a multiband classification, performed in this work comparing two different approaches: a neuro-fuzzy per-pixel Fuzzy ARTMAP (FA) classifier [7], and a contextual classifier based on the assumption of a Markov random field (MRF) spatial model [8]. Generally speaking, the neuro-fuzzy classifier has been chosen because of its proven capability to provide good 4 International Journal of Navigation and Observation  results when performing a multiband per pixel classification, while the MRF approach allows a spatially joint analysis. Both algorithms are supervised, and as such a minimal knowledge about the damages on the ground, their locations, and level is required.
After classification, and due to the complex interactions between radar waves and the urban environment, either damaged or undamaged, it is very likely that the damage classification map has a "blurred" or "fuzzy" appearance. To improve the results, and to meaningfully focus the damage map at a spatial scale of interest to the final user, a fusion step between the map results and ancillary GIS data is performed. To this aim, the best results between the two classification methods are used to make a decision, by means of a decision fusion process, about each of the areas detailed by the GIS ancillary information. This processing step involves a data fusion procedure which has been detailed in [5] and in this work will be degraded to the simplest situation, that is, majority voting. This basically means that the most voted damage class in each block individuated by the GIS layer is considered as representative of the whole block.
Following this procedure, in next section two different testcases will be considered, referring to very different countries in different parts of the world. Moreover, different SAR sensors are considered, and thus different scale of analysis and data availability. With the results of the following section we want to stress how much the simple procedure presented here can be helpful in real situations, and compare how much and how well different choices of the input feature set can highlight the damage patterns in the area. The goal of having more test cases and comparing with the originally studied Bam (Iran) and Golcük (Turkey) cases is also a way to check for the best combination (if any) for all of them.

APPLICATIVE TEST CASES
In order to test the proposed methodology, the aim of our tests was to consider different SAR datasets, coming from sensors on board of different satellites. The combination of different bands of work, spatial resolution, polarization information, availability or not of phase information was meant to provide a method to test the robustness of the approach and find where it has to be adapted. Moreover, the very diverse damage patterns, connected to the original spatial urban patterns and the effect of the earthquake, make the two test site analyzed in the following. (in addition to the 2003 Bam test case in [5] and the 1999 Turkey case in [6] a definitely valid series of applicative results.)

First test site: Boumerdes (2003 Algeria earthquake)
The first results refer to the magnitude 6.

Second test site: Pisco (2007 Peru earthquake)
The second example refers to the test case of Peru, whose consist of a GIS layer depicting the borders of the parcels in the urban area, and were obtained by manual digitalization of the information in [9], and validated by comparison with the same SPOT pre-event image used in that paper. From the same paper, also the information related to damaged areas obtained by in situ measurements was extracted (see Figure 5). Since the data have been provided as amplitude images, no phase information could have been considered. This prevents us from using the bands which were considered as G. Trianni and P. Gamba 7 the best choice in [5] for the problem of damage detection, that is, pre-and postevent intensity, pre-post coherence and difference between pre-post and pre-pre coherence. Instead, from the available images only some intensity features have been extracted, in particular the intensity correlation r and the backscattering coefficient d and taking into account the different spatial resolution of the ALOS/PALSAR scene than the ERS and JERS data used in all the works on the same subject so far. Other considered features are the pre-and postevent intensities, computed from the original data after despeckling with a 5 × 5 Gamma filter. As in the first test case the first classification by means of the above neuro-fuzzy perpixel classifier or the context-aware MRF classifier is followed by a fusion with the GIS layer, where each parcel of the GIS is assigned to the class to which the majority of mapped pixels belong.
For this test case, a wider range of classification maps is proposed in Figures 6-9 to allow a better comparison of the combinations of features and classifiers, as well as to appreciate the improvement in understanding the results by using the data fusion final step. The need for this extended result analysis is connected to the different spatial resolution of ALOS/PALSAR with respect to the ERS/JERS data used so far. This makes the classification map more precise at the per-pixel level and allows defining spatial units smaller than in the previous test site (compare Figures 5(b) with 3(b)). Moreover, the lack of original complex data does not allow computing the phase coherence, one of the most important bands for the multiband/multitemporal damage assessment classification step according to [5]. We thus intend to analyze which, among the computed features, are the most viable ones to get a rapid damage map in this situation.
Finally, Table 1 reports the overall accuracy for the maps in the rightmost column in Figures 6-9 and allows to improve the visual comparison with a quantitative assessment.
A first comment to the results is that the information fusion step is really mandatory to achieve results not only with a decent mapping accuracy, but also understandable to anyone looking to the map. A second comment is that the accuracy values are in the same range as the one reported for the first test case, the Algeria earthquake. Although the ground truth in the present case is more detailed, the higher spatial resolution of the SAR data allow matching the accuracy values obtained from Boumerdes' images.
According to the maps and the accuracy values, the best result is obtained by using a combination of the amplitude pre-event image (both polarizations HH and HV) and the amplitude postevent image (again with both polarizations), and the second best approach is the use of the backscattering and the correlation computed between the pre-event image (HH polarization) and the postevent image (HV polarization), and the postevent image (with both polarizations).
Since the ALOS/PALSAR image pre-and postevent image pairs are alternate polarization (AP) images, it was also possible to compare the effect of polarization with respect to damage mapping task using the proposed approach. However, it was found that no particular choice can be made, and the problems in mapping the damage to the right extent in some portion of the area are equally in place suing one or the other of the two polarizations, or even a combination of both.
A last comment is driven by the fact that the damage maps in [9] involve four classes, instead of the three used in our validation. In Figure 9, in fact, medium and high damages were considered as a single class, in orange. The reason is that previous trials attempting to obtain maps with high level of damage discrimination did fail. The highest accuracy values, as reported in Table 1 for sake of completeness, were around 45% at their best. The corresponding maps, where a high level of damage is depicted in red, are reported in the last row of Figures 6 and 8. Apparently, thereare limits inherent to the structure of the artificial elements of the landscape, the typology of earthquake, and the spatial resolution that prevent the 7 m ALOS/PALSAR to be enough for this precise recognition task.

CONCLUSIONS
In this paper, a rapid damage mapping approach is proposed, based on the exploitation and interpretation of satellite SAR data. The approach proves to be robust and useful to detect the damage patterns. It is however imprecise with respect to accuracy and needs improvements.
More specifically, although semiautomatic SAR data interpretation in urban areas is an open research issue, this paper shows that a combination (fusion) of remotely sensed data and geographical databases may lead to a real improvement in this interpretation, making the data more useful for the end-user. Accuracy and robustness of the procedure, together with its affordability, were proved by the analysis of extreme events (notably, earthquakes) in many different parts of the world.
There are some commonalities among the choices of input features used in the cases presented in this work, and the combination of the pre-and postevent intensity data with the intensity correlation and the difference of the logarithmic means achieves always better results. The possibility to incorporate some phase information by means of coherence, not exploited in this work but proposed in the original paper [5], leads apparently to the best among these multiple possible choices.
Very interesting and still open issues are those connected to the analysis of additional spatial feature and the correlation between the maps and the features itself to the actual, on site generated, damage map.

ACKNOWLEDGMENTS
The first part of this work has been carried out during the stay of G. Trianni at the John A. Blume Earthquake Engineering Center, Stanford University, Calif, USA, as a visiting researcher. The authors would like to thank Professor Anne Kiremidjian and Dr. Pooya Sarabandi for their hospitality and M. Matsuoka and Professor F. Yamazaki for providing the SAR dataset and the damage map for the Boumerdes test site. The PALSAR satellite data over Peru were provided by