Sequential SAR Coherence Method for the Monitoring of Buildings in Sarpole-Zahab, Iran

: In this study, we used ﬁfty-six synthetic aperture radar (SAR) images acquired from the Sentinel-1 C-band satellite with a regular period of 12 days (except for one image) to produce sequential phase correlation (sequential coherence) maps for the town of Sarpole-Zahab in western Iran, which experienced a magnitude 7.3 earthquake on 12 November 2017. The preseismic condition of the buildings in the town was assessed based on a long sequential SAR coherence (LSSC) method, in which we considered 55 of the 56 images to produce a coherence decay model with climatic and temporal parameters. The coseismic condition of the buildings was assessed with 3 later images and normalized RGB visualization using the short sequential SAR coherence (SSSC) method. Discriminant analysis between the completely collapsed and uncollapsed buildings was also performed for approximately 700 randomly selected buildings (for each category) by considering the heights of the buildings and the SSSC results. Finally, the area and volume of debris were calculated based on a fusion of a discriminant map and a 3D vector map of the town. of speckle noise in SLC images and increase the reliability of the coherence maps, we considered the multilook factor of range 4 and azimuth 1 to obtain square pixels in the final product of the Sentinel-1 data. Once coherence values over an urban area are calculated for at least one preseismic pair, the concept of damage can be developed by calculating the coseismic coherence values of one post-event image and one of the pre-event images that were used for the preseismic coherence analysis. A decorrelation between the pre-event and post-event data manifests itself as a decreased coherence map. To understand the extent of the damage, simple differencing and normalized differencing can be considered as follows:


Introduction
One of the necessities of disaster management before or during a natural or anthropogenic event is the inexpensive monitoring of buildings in urban areas. Earthquakes are disastrous natural events that must be constantly monitored to improve the damage response and to prevent numerous deaths. Several spaceborne remote sensing datasets and methods have been proposed to achieve these goals [1][2][3][4][5][6][7][8][9][10][11]. Because earthquakes are significant threats to cities, for successful post-disaster management and rapid damage response, the remote sensing datasets and their derived products (e.g., the extracted footprints of buildings, building heights and corresponding pixel values) should interact with a geographic information system (GIS) that incorporates different characteristics of the city and the remote sensing data. The two main remote sensing methods for monitoring earthquake-induced building damage are based on optical and synthetic aperture radar (SAR) data. One example is the application of optical imagery to building damage assessments after the Gujarat earthquake (2001) using pre-event and post-event spaceborne Landsat-panchromatic images [2]. However, optical imagery is not always successful, especially when the weather is not clear. SAR imagery has the advantage of all-weather imaging capability during both the day and the night. The intensity information (e.g., difference and correlation of the backscattering coefficient) of 5 pre-event and 4 post-event SAR images acquired by the ERS-1 satellite were used to assess damaged and undamaged buildings after the 1995 Kobe earthquake [3]. This method was also used with the ENIVSAT satellite's SAR images and the Bam earthquake (2003) [4]. The methodology has also been successfully applied in other case studies; however, due to the lack of sufficient ground truth data from around the world, the model sometimes does not provide accurate damage estimation, especially if the backscattering reflections from the rooftops of buildings are different. To evaluate the building damage after the 2006 Mid Java earthquake, three images were used to calculate preseismic and coseismic correlation maps of the study area using the ALOS-1 satellite's PALSAR-1 data [5]. In addition to intensity information, phase information plays an important role in the detection of earthquake-induced building damage. The Kobe and Bam earthquakes were evaluated based on the phase correlation (coherence) of SAR images from the JERS-1 satellite and ENVISAT satellite, respectively [6,7]. Coherence is the spatial correlation of the interferometric phases of two SAR images and is controlled by the position and stability of individual scatterers. Thus, the method requires at least three images to create two coherence maps to assess the pre-event and coseismic conditions of the affected areas. The coherence method has become popular and is accurate because it minimizes seasonal surface changes and signal noise, especially if repeated observations over a specific area are available. However, the retired SAR missions that were mentioned previously (i.e., ERS, ENVISAT, JERS-1, ALOS-1) had longer revisit intervals over specific areas, which made accurate and timely damage assessment difficult. In addition, their spatial resolution was not sufficient for building-by-building damage assessments. New SAR missions such as the TerraSAR-X, COSMO-SkyMed and Sentinel-1 missions, which have considerably shorter revisit intervals and better spatial resolution, provide valuable information about the Earth's surface. Sentinel-1 s coverage is better than those of the other missions. It is the first satellite with an interferometric wide swath (IW) mode that uses operationally the terrain observation with progressive scan (TOPS) technique. The Sentinel-1A satellite, which is part of the Sentinel-1 mission, was launched on 3 April 2014 to collect images of the Earth using a C-band SAR sensor with a 12-day revisit interval. The revisit interval of the mission was reduced to 6 days after the launch of the Sentinel-1B satellite almost two years after Sentinel-1A. From the start of the mission to the end of 2016, a total of 884,628 Copernicus Sentinel-1 products were made available on the Copernicus Open Access portal, which have a total data volume of 0.77 petabyte (PB) [12]. In the era of SAR dataset abundance, the monitoring of stationary (e.g., seasonal effects) and nonstationary (e.g., earthquake-related) changes on the Earth is now easier than it was in previous years.
Recently, the impacts of coherence changes during the Iran-Iraq earthquake (2017) on the town of Sarpole-Zahab were evaluated using fifteen Sentinel-1 images with temporal baselines less than 100 days and perpendicular baselines less than 150 m [10]. The method used one image as the master image and its coregistration to subsequent slaves. The pre-event SAR interferometry (InSAR) analysis exhibited a coherence loss of 37.8% in the town, whereas the coherence loss for the coseismic analysis was 70.3% and was related to the destruction of roads and buildings. In this study, we analyze fifty-six SAR images acquired from Sentinel-1 in the dual-polarized mode (VV and VH) from January 2016 to November 2017 at a regular interval of 12 days to fill the gap between the previous studies and show the approximately 2-year changes in the town of Sarpole-Zahab using InSAR before and during the earthquake. The method used in this study is an urban monitoring technique that uses the chain-like coherence of SAR images with respect to each previous image to reveal the preseismic and coseismic conditions of the town. The chain-like coherence method is composed of two parts: The long sequential SAR coherence (LSSC) monitors the pre-event conditions and reveals the seasonal effects, and the short sequential SAR coherence (SSSC) monitors the coseismic changes. We model the coherence values with respect to the elapsed days for each pair and the amount of precipitation in Sarpole-Zahab to constrain the decrease of the coherence. The proposed model addresses the coherence decays related to precipitation and the earthquake. It also predicts the 1-year coherence changes in the absence of human activities under the post-event conditions. The earthquake damage assessment is also performed by (1) normalized RGB color composition of the two latter coherence maps and (2) logistic discriminant analysis of uncollapsed and collapsed buildings considering both the coherence values and the building heights. Finally, the volume of debris and the 3D perspective of the town are calculated based on the discriminant cut-off values, the building footprints and their corresponding y values.

Earthquake and Study Area
Iran is located in a seismically active region between two tectonic plates to the north (Eurasia) and the south (Arabia). In recent decades, Iran has experienced several deadly earthquakes that caused considerable casualties, which makes a rapid observation system necessary to increase disaster response and improve emergency management in the aftermath of an earthquake. Recent major earthquakes in different parts of Iran include the 2003 Bam earthquake (M 6.3) in central Iran, which caused approximately 40,000 causalities, the twin earthquakes in Ahar-Varzaghan in NW Iran (M 6.4 and 6.4) in 2012, which caused more than 300 casualties, and the 2013 Makran earthquake (M 7.3) in the SE subduction zone [13][14][15]. More recently, at 21:48 Iran Standard Time on 12 November 2017, a magnitude 7.3 earthquake occurred in the Zagros seismotectonic zone in western Iran at a depth of 18 km and particularly affected Kermanshah Province [16,17]. This event was the deadliest earthquake in 2017; it caused more than 600 casualties, injured 7000 people and left 15,000 people homeless, and it was not only felt in the neighboring provinces of Iran but was also recorded by seismological stations in neighboring countries, such as Iraq, Turkey and Azerbaijan [16]. According to the Iranian seismic network, the earthquake occurred in Iranian territory in Kermanshah Province, near the Iraqi border. The event was called the Kermanshah or Sarpole-Zahab earthquake in Iran; in Iraq, it was called the Halabja earthquake. In the early hours after the earthquake, no structural damage was reported in towns located north of the epicenter (e.g., Paveh and Nowsud). In contrast, the main damage reported by visual inspections conducted by experts and local residents was located in towns south of the epicenter, such as Salase-Babajani, Sarpole-Zahab and Qasre-Shirin [16,17]. As shown in Figure 1, the epicenter was located 28, 51 and 56 km from Salase-Babajani, Sarpole-Zahab and Qasre-Shirin, respectively. The highest peak ground acceleration (PGA) record of the Iranian Strong Motion Network, which is operated by the Building and Housing Research Center (BHRC), was at the Sarpole-Zahab station, where the longitudinal (L), traverse (T) and vertical (V) components were 681 cm/s 2 , 562 cm/s 2 and 404 cm/s 2 , respectively. The event's PGA records at the Nowsud station, which is located 37 km NW of the epicenter, were 48 cm/s 2 , 45 cm/s 2 and 44 cm/s 2 for the L, T and V components, respectively, which indicates that the seismic waves propagated from north to south [16]. Because the Sarpole-Zahab station is located in a forward direction, its high PGA record is reasonable and implies that the shorter distance between the site and the epicenter cannot exclusively reflect the nature of damage to buildings. If we assume that the quality and construction methods of the structures in the towns north and south of the epicenter are similar, the geological/topographical conditions and shear wave velocity of the upper 30 m (Vs30) are the most important factors that affected the heavy damage in the towns south of the epicenter. A shaded relief map of the study area indicates that Sarpole-Zahab and Qasre-Shirin are located in relatively flat regions, whereas the majority of the towns to the north of the epicenter, such as Paveh, are located in topographically rugged regions, which may reduce the impacts of local seismic amplification ( Figure 1). The average Vs30 values for Nowsud and Sarpole-Zahab are approximately 760 m/s and 290 m/s, respectively. These Vs30 values were extracted from Iran's national Vs30 map, which was developed from topographic and geological data as a proxy for seismic amplification and indicates that sites with lower Vs30 values have greater potential of being damaged during earthquakes [18]. Although the purpose of this study is to demonstrate the feasibility of using SAR data for rapid monitoring of urban regions, large-scale damage monitoring of the stricken area using SAR data is not currently possible because auxiliary data about building heights from commercial sources, such as digital surface models (DSMs), are not available. Thus, we selected an area of approximately 7 km 2 around Sarpole-Zahab for a pilot study due to the heavy damage in some parts of the town and the availability of test data for the collapsed buildings from other sources. The bright pixels in Figure 2 are buildings in Sarpole-Zahab overlaid on a backscatter map (30 October 2017); 75% of the buildings in the town are 1-or 2-story buildings, and the remaining buildings are 3-to 6-story buildings. The bright pixels that are extended NW-SE are not associated with buildings; they correspond to a layover distortion of sharp topographic features (e.g., mountains), in which the foreslope angle (α = 49 • ) of the features exceeds the incidence angle of the Sentinel-1 image (θ = 42 • ). The layover distortion cannot be removed from the image because the radar receivers receive the signals simultaneously, regardless of how far apart the locations are. Therefore, they are not separated in the superimposed image [19]. The dark pixels in the western area with low intensity values are vegetated areas, cropland and filled lands, which are classified as soft soil in most building construction codes.
The sequential SAR coherence method used in this study is derived from the concept of repetitive SAR data and the conventional InSAR method. The techniques used in InSAR rely on phase differences of different epochs (i.e., ground displacements). Interferometric methods can provide information such as the coherence, which is an index used to measure the degree of interferometric correlation between two coregistered SAR images as follows: where and are the master and slave images, respectively; the brackets denote averaging over a finite number of signal measurements, and and ̅ are the complex conjugates of the master and slave images, respectively. There is a trade-off between the window size and parameters such as the ground resolution and coherence. As the window size increases during phase processing, the resolution of the results will decrease, whereas small windows will result in overestimated coherence values. To decrease the level of speckle noise in SLC images and increase the reliability of the coherence maps, we considered the multilook factor of range 4 and azimuth 1 to obtain square pixels in the final product of the Sentinel-1 data. Once coherence values over an urban area are calculated for at least one preseismic pair, the concept of damage can be developed by calculating the coseismic coherence values of one post-event image and one of the pre-event images that were used for the preseismic coherence analysis. A decorrelation between the pre-event and post-event data manifests itself as a decreased coherence map. To understand the extent of the damage, simple differencing and normalized differencing can be considered as follows:
The sequential SAR coherence method used in this study is derived from the concept of repetitive SAR data and the conventional InSAR method. The techniques used in InSAR rely on phase differences of different epochs (i.e., ground displacements). Interferometric methods can provide information such as the coherence, which is an index used to measure the degree of interferometric correlation between two coregistered SAR images as follows: where M and S are the master and slave images, respectively; the brackets denote averaging over a finite number of signal measurements, and M and S are the complex conjugates of the master and slave images, respectively. There is a trade-off between the window size and parameters such as the ground resolution and coherence. As the window size increases during phase processing, the resolution of the results will decrease, whereas small windows will result in overestimated coherence values. To decrease the level of speckle noise in SLC images and increase the reliability of the coherence maps, we considered the multilook factor of range 4 and azimuth 1 to obtain square pixels in the final product of the Sentinel-1 data. Once coherence values over an urban area are calculated for at least one preseismic pair, the concept of damage can be developed by calculating the coseismic coherence values of one post-event image and one of the pre-event images that were used for the preseismic coherence analysis. A decorrelation between the pre-event and post-event data manifests itself as a decreased coherence map. To understand the extent of the damage, simple differencing and normalized differencing can be considered as follows: For local urban studies, simple differencing, as calculated in Equation (2), results in satisfactory accuracy; however, RGB visualizations are better applied in the normalized approach for potential global approaches in future studies [9]. Although our study is local, we follow both simple and normalized approaches; we use the former to evaluate the coherence values of the pairs with other effective parameters and the latter for mere RGB visualization. Generally, the use of pre-event images with shorter temporal gaps in coseismic InSAR analysis will increase the reliability and quality of the coherence proxy maps [8]. In this study, we use this idea not only for the coseismic coherence but also for all of the datasets acquired from Sentinel-1. As shown in Figure 3a, regardless of the relative spatial baseline of the images, we focus only on sequential coherence analysis because the images are regularly obtained over time. The long sequential SAR coherence (LSSC) and short sequential SAR coherence (SSSC) can be defined as follows: Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 17 For local urban studies, simple differencing, as calculated in Equation (2), results in satisfactory accuracy; however, RGB visualizations are better applied in the normalized approach for potential global approaches in future studies [9]. Although our study is local, we follow both simple and normalized approaches; we use the former to evaluate the coherence values of the pairs with other effective parameters and the latter for mere RGB visualization. Generally, the use of pre-event images with shorter temporal gaps in coseismic InSAR analysis will increase the reliability and quality of the coherence proxy maps [8]. In this study, we use this idea not only for the coseismic coherence but also for all of the datasets acquired from Sentinel-1. As shown in Figure 3a, regardless of the relative spatial baseline of the images, we focus only on sequential coherence analysis because the images are regularly obtained over time. The long sequential SAR coherence (LSSC) and short sequential SAR coherence (SSSC) can be defined as follows: In Equations (4) and (5), represents the label of the images, whereas the label of the sequential coherence is − 1. The spatiotemporal distribution of the baselines for a timespan of two years (9 In Equations (4) and (5), i represents the label of the images, whereas the label of the sequential coherence is i − 1. The spatiotemporal distribution of the baselines for a timespan of two years (9 January 2016 to 17 November 2017) of Sentinel-1A and B (Figure 3a) shows that for a pool of 57 interferometric wide swath (IW) images with a 5 × 20 m spatial resolution from an ascending orbit (72), 56 coherence maps with VV polarization are produced. Depending on the preseismic or coseismic conditions, 54 and 2 coherence maps are considered for the LSSC and SSSC, respectively. All of the epochs have a regular 12-day latency except for one 24-day epoch (between 13 June 2016 and 7 July 2016). Because this 24-day epoch occurred during the dry season in the study area, its impact on the general trend of the mean coherence is negligible. Accordingly, the mean coherence of the LSSC can be used as an indicator of stationary changes, such as seasonal effects and precipitation. As shown in Figure 3a Although the mean coherence value decreases during these winter-to-spring periods, histograms of the coherence in Sarpole-Zahab show that the high coherence values are influenced less by seasonal effects than by the earthquake (Figure 3c). Figure 3c shows that after the earthquake, the peak coherence shifted considerably toward lower values, whereas the seasonal effects flattened the coherence peaks without any considerable shift. Note that we did not exclude some small vegetated areas, such as parks and promenades in Sarpole-Zahab, during the mean coherence calculation. The most important characteristic of the decrease in coherence of the buildings in Sarpole-Zahab in the winter-to-spring periods is its indirect relationship with the height of buildings. This relationship was extracted from a comparison between the heights of all of the buildings in Sarpole-Zahab (9230 buildings) and the corresponding centroid coherence values (Figure 3d). The correlation coefficient between the coherence and the building height is 0.21, which implies that as the building height increases, the InSAR results become more coherent, which might be related to the type of materials used in rooftops. This correlation is not only found for the height of buildings; other height-derived parameters, such as the volume of buildings, have positive correlations with the coherence. Most new buildings in Sarpole-Zahab are relatively tall (between 3 and 6 stories), whereas the old buildings are usually 1 to 2 stories tall. During a field survey in December 2017, we observed that the materials used in the rooftops of old buildings are more natural than the synthetic materials used in new buildings ( Figure 4). Natural building materials created from soil, water and fibrous products such as straw are dominant in old buildings. Usually, these cob-made rooftops absorb water during rain, which changes the volume of the roof; thus, the coherence in this type of building would be relatively low. In contrast, as shown in Figure 4d, the walls and rooftops of new buildings are composed of bituminous materials to protect the insides of the buildings from rain and moisture. Bituminous waterproofing sheets do not change in size during rain and exhibit higher coherence values. Although a correlation between the building height and coherence exists, it is not a high correlation. One local reason might be the replacement of cob materials in the old buildings with bituminous materials in recent years ( Figure 4a). Old buildings were built in different locations with dark rooftops, whereas new buildings have been constructed in the western part of the town and can be distinguished by silver rooftops (Figures 2 and 4). Unfortunately, because of our incomplete building inventory, discrimination of all of the repaired old buildings was not possible. Because we focus on LSSC and SSSC for the coherence monitoring and building damage assessment in this study, we did not present the backscattering coefficients of the buildings. The backscattering behavior is complex; during rainy periods, it would be high in old buildings because of the electrical properties of air and water, whereas in the dry seasons, it would be lower in old buildings because some of the radar energy will be absorbed by the dry soil.

Coherence Modeling
The dynamic state of SAR data includes both stationary and nonstationary parts. In Sarpole-Zahab, from the point of view of stochastic processes, the stationary elements, such as climatic and seasonal conditions (i.e., LSSC), can be represented by coherence decay models, whereas nonstationary features, such as the impacts of human activities after disasters, are more difficult to model. One reason why we cannot focus on the coherence modeling of human changes (i.e., SSSC) is that these changes are highly dependent on local parameters, which are more complex. For example, the post-earthquake modeling of the coherence changes in Sarpole-Zahab requires extensive information about the post-earthquake construction activities, debris removal tasks, and other activities. Here, we focus only on the phase stability of Sarpole-Zahab according to stationary parameters in which the temporal decorrelation is reduced with regular Sentinel-1 acquisitions. We use the mean coherence value of all of the Sarpole-Zahab buildings for which the timespan of the results of 660 days contains two winter-to-spring periods. The mean coherence shows that in a normal stationary condition, the minimum and maximum mean coherence values are 0.86 and 0.75, respectively. A stationary process as a cause of coherence changes can be represented by a polynomial temporal-climatic (TC) model that includes nonlinear (precipitation) and linear (time difference) parts as follows: where ( , ) is the modeled coherence of the sequential method, is precipitation, is the time differential between image − 1 and image , which is 12 days (or more) for Sentinel-1A products, and − are coefficients estimated from least squares and iteration methods. We impose a constraint on the nonlinear part of the TC model such that should be nonzero. This constraint does not change the outcomes significantly, but it forces the TC model to produce logical coherence values that are not greater than one [25]. Thus, the TC model is useful for separating the seasonal effects from the earthquake effects. Figure 5 shows that the TC model fits the calculated results well, and the overall RMSE is 0.02. The RMSE in the winter-to-spring period of 2016 (i.e., days 0-156) was 0.007, whereas the RMSE in the winter-to-spring period of 2017 (i.e., days 348-492) was 0.03. The TC model in the summer to autumn of 2016 (i.e., days 157-336) tends to overestimate the mean coherence

Coherence Modeling
The dynamic state of SAR data includes both stationary and nonstationary parts. In Sarpole-Zahab, from the point of view of stochastic processes, the stationary elements, such as climatic and seasonal conditions (i.e., LSSC), can be represented by coherence decay models, whereas nonstationary features, such as the impacts of human activities after disasters, are more difficult to model. One reason why we cannot focus on the coherence modeling of human changes (i.e., SSSC) is that these changes are highly dependent on local parameters, which are more complex. For example, the post-earthquake modeling of the coherence changes in Sarpole-Zahab requires extensive information about the post-earthquake construction activities, debris removal tasks, and other activities. Here, we focus only on the phase stability of Sarpole-Zahab according to stationary parameters in which the temporal decorrelation is reduced with regular Sentinel-1 acquisitions. We use the mean coherence value of all of the Sarpole-Zahab buildings for which the timespan of the results of 660 days contains two winter-to-spring periods. The mean coherence shows that in a normal stationary condition, the minimum and maximum mean coherence values are 0.86 and 0.75, respectively. A stationary process as a cause of coherence changes can be represented by a polynomial temporal-climatic (TC) model that includes nonlinear (precipitation) and linear (time difference) parts as follows: where γ d (t i−1 , t i ) is the modeled coherence of the sequential method, p is precipitation, t is the time differential between image i − 1 and image i, which is 12 days (or more) for Sentinel-1A products, and a − e are coefficients estimated from least squares and iteration methods. We impose a constraint on the nonlinear part of the TC model such that a should be nonzero. This constraint does not change the outcomes significantly, but it forces the TC model to produce logical coherence values that are not greater than one [25]. Thus, the TC model is useful for separating the seasonal effects from the earthquake effects. Figure 5 shows that the TC model fits the calculated results well, and the overall RMSE is 0.02. The RMSE in the winter-to-spring period of 2016 (i.e., days 0-156) was 0.007, whereas the RMSE in the winter-to-spring period of 2017 (i.e., days 348-492) was 0.03. The TC model in the summer to autumn of 2016 (i.e., days 157-336) tends to overestimate the mean coherence values, whereas during the 2017 summer-to-autumn period (i.e., days 494-660), the TC model underestimates the mean coherence value until the sequential coherence chain is broken by the earthquake between days 660 and 678. Table 1 shows the summary of the statistical parameters, the estimated coefficients and the sensitivity of the TC model. values, whereas during the 2017 summer-to-autumn period (i.e., days 494-660), the TC model underestimates the mean coherence value until the sequential coherence chain is broken by the earthquake between days 660 and 678. Table 1 shows the summary of the statistical parameters, the estimated coefficients and the sensitivity of the TC model.

Damage Mapping and Discriminant Analysis of SSSC
A visual comparison of the collapsed buildings from UNITAR's (United Nations Institute for Training and Research) report with the results of the simple differential coherence (∆ ) method shows that the method enables the detection of collapsed buildings. As shown in Figure 6a,b, most of the collapsed buildings are located in the NW, SW and central parts of the town, but the range of ∆ cannot be quantified by simple inspections. The normalized RGB color composition map of the study area provides a way to more realistically visualize the extent of the damage, in which the forward value of ∆ (C54-C55) in the red band shows collapsed buildings, the backward value of ∆ (C55-C54) in the green band shows human activities, such as debris removal that occurred immediately after the event, and the mean values of the SSSC in the blue band show the urban pixels without distinct changes (Figure 6c). The normalized RGB map of Sarpole-Zahab is a good example of a qualitative damage map that presents both earthquake-related and human activities on the same map. However, such mapping does not always produce acceptable results depending on the density of the vegetated and built up areas and the wavelength of the microwave sensors used [9].

Damage Mapping and Discriminant Analysis of SSSC
A visual comparison of the collapsed buildings from UNITAR's (United Nations Institute for Training and Research) report with the results of the simple differential coherence (∆γ) method shows that the method enables the detection of collapsed buildings. As shown in Figure 6a,b, most of the collapsed buildings are located in the NW, SW and central parts of the town, but the range of ∆γ cannot be quantified by simple inspections. The normalized RGB color composition map of the study area provides a way to more realistically visualize the extent of the damage, in which the forward value of ∆γ (C54-C55) in the red band shows collapsed buildings, the backward value of ∆γ (C55-C54) in the green band shows human activities, such as debris removal that occurred immediately after the event, and the mean values of the SSSC in the blue band show the urban pixels without distinct changes (Figure 6c). The normalized RGB map of Sarpole-Zahab is a good example of a qualitative damage map that presents both earthquake-related and human activities on the same map. However, such mapping does not always produce acceptable results depending on the density of the vegetated and built up areas and the wavelength of the microwave sensors used [9]. To further investigate the extent of damage, we apply discriminant analysis to create a damage proxy map of Sarpole-Zahab. This method was used in previous studies with the SAR correlation coefficient and the difference in backscattering values [3,4] and, more recently, with the difference in coherence and the difference in backscattering [8,9]. Here, we apply discriminant analysis based on ∆ estimated from the SSSC and the height of buildings (H). A 3D vector map of the buildings in Sarpole-Zahab with an acceptable vertical and horizontal accuracy of ~0.5 m is created from seven VHR optical images from the latest acquisitions of DigitalGlobe's constellation between 2015 and 2017 [26][27][28]. The 3D vector map of the town over an area of ~15 km 2 includes both the footprints and heights of the buildings, which are used to compare the centroid coherence values of each individual building with the corresponding heights in Figure 3d. The discriminant analysis considers two types of buildings (two dependent variables) to be discriminated: (1) collapsed buildings, which are those in which the walls or rooftop have been significantly deformed (e.g., Grades 4 and 5 in EMS-98); and (2) uncollapsed buildings (e.g., Grades 1 to 3 in EMS-98), which are those that are either intact or slightly deformed with damage not visible from satellite imagery [29]. These variables are assigned a value of 1 if severe damage characteristics are present and 0 if severe damage characteristics are not present. The locations of the collapsed buildings in Sarpole-Zahab from the high-resolution optical imagery prepared by UNITAR are used as calibration data in our discriminant analysis; thus, the aptitude tests of ∆ (C54-C55) and the building heights are two representatives of the discriminant function [30]. According to UNITAR's report, 683 collapsed buildings were identified. To prevent a violated discriminant analysis and to maintain a normal distribution between the collapsed and uncollapsed buildings, we also randomly select 683 uncollapsed buildings (Figure 6a). Figure 6a shows the locations of 683 collapsed buildings (red circles) from UNITAR and 683 randomly selected To further investigate the extent of damage, we apply discriminant analysis to create a damage proxy map of Sarpole-Zahab. This method was used in previous studies with the SAR correlation coefficient and the difference in backscattering values [3,4] and, more recently, with the difference in coherence and the difference in backscattering [8,9]. Here, we apply discriminant analysis based on ∆γ estimated from the SSSC and the height of buildings (H). A 3D vector map of the buildings in Sarpole-Zahab with an acceptable vertical and horizontal accuracy of~0.5 m is created from seven VHR optical images from the latest acquisitions of DigitalGlobe's constellation between 2015 and 2017 [26][27][28]. The 3D vector map of the town over an area of~15 km 2 includes both the footprints and heights of the buildings, which are used to compare the centroid coherence values of each individual building with the corresponding heights in Figure 3d. The discriminant analysis considers two types of buildings (two dependent variables) to be discriminated: (1) collapsed buildings, which are those in which the walls or rooftop have been significantly deformed (e.g., Grades 4 and 5 in EMS-98); and (2) uncollapsed buildings (e.g., Grades 1 to 3 in EMS-98), which are those that are either intact or slightly deformed with damage not visible from satellite imagery [29]. These variables are assigned a value of 1 if severe damage characteristics are present and 0 if severe damage characteristics are not present. The locations of the collapsed buildings in Sarpole-Zahab from the high-resolution optical imagery prepared by UNITAR are used as calibration data in our discriminant analysis; thus, the aptitude tests of ∆γ (C54-C55) and the building heights are two representatives of the discriminant function [30]. According to UNITAR's report, 683 collapsed buildings were identified. To prevent a violated discriminant analysis and to maintain a normal distribution between the collapsed and uncollapsed buildings, we also randomly select 683 uncollapsed buildings (Figure 6a). Figure 6a shows the locations of 683 collapsed buildings (red circles) from UNITAR and 683 randomly selected uncollapsed buildings (blue circles), and Figure 6d shows that the ∆γ value of the collapsed buildings (0.24) is nearly twice that of the uncollapsed buildings (0.11). Because some of the collapsed buildings in Sarpole-Zahab had areas smaller than 100 m 2 , the higher standard deviation of the collapsed buildings is an indication of a limitation of the coarse resolution SAR imagery, in which certain small buildings with severe damage are not identified by the SSSC method. In the discriminant analysis, we try to predict the probability of the buildings falling into categories 1 or 0 using linear regression. However, the interpretation of the outputs of this method is different from that applied for ordinary least squares. In the ordinary least squares method, it is difficult to force the predicted value to remain between 0 and 1. Thus, we apply the logistic regression method to account for the possibility of being 0 and 1 and fix the violation of dummy values with constraints. In the ordinary method for the two-category case, the dependent variables are defined as follows: where y is the categorical dependent variable and an indicator of the degree of building damage. If y is 1, the observation falls into category 1. If it does not fall into category 1, y equals 0, β 0 is a constant, and β 1 through β i are regression coefficients. The parameters x 1 through x i are independent variables inferred from each individual building. For the logistic regression, Equation (7) can be rewritten as follows: where p is the probability that a collapsed building is present. In the logistic method we defined, the probability of the logarithmic value of the odds of the buildings is a linear combination of the independent variables (H and ∆γ) as follows: when an x variable changes by 1 unit, and all the other factors remain unchanged, the odds will increase by e b i . Once y(p) is calculated based on an interaction between the dependent and independent variables, we define a constraint on the sum of all of the possibilities to calculate the intercept (β 0 ) and the coefficients of the independent variables (β i ). Here, we solve a new function to maximize the sum of the logarithmic likelihoods of the building damage type (i.e., uncollapsed and collapsed) and find the related parameters as follows: To create a damage proxy map from the logistic equation (Equation (11)), we calculated the threshold between the uncollapsed and collapsed groups as follows: where T is the threshold, n 0 and n 1 are the numbers of uncollapsed and collapsed buildings, respectively, and y 0 and y 1 are the average values of the averaged discriminant scores for uncollapsed and collapsed buildings, respectively. The calculated threshold value is 0.07, which means that buildings with a discriminant score (y) lower than 0.07 are categorized in group 0 and that those with a discriminant score (y) higher than 0.07 are categorized in group 1. To compare the results of the discriminant score with the calibration data from UNITAR, we randomly selected uncollapsed buildings to calculate the overall accuracy of the logistic discriminant analysis. The overall accuracy of the building classification in Sarpole-Zahab is 68%, in which 932 of the 1366 buildings from the training data are correctly classified. This result indicated slightly higher accuracy than that achieved through the simple discriminant analysis. The number of misclassified buildings in categories 0 and 1 were 144 and 290, respectively. Details of the building classification with respect to logistic regression are given in Table 2. Figure 7 shows the damage proxy map of Sarpole-Zahab with scatter plots of ∆γ and H, in which the majority of the collapsed buildings were two-story buildings. This finding does not mean that the 3-to 6-story buildings were intact or slightly damaged. The scatterplot of the collapsed (red circles) buildings in Figure 7 indicates that almost all of the buildings with heights more than 20 m in Sarpole-Zahab collapsed. The higher number of misclassifications among collapsed buildings shown by the thick dashed line in Figure 7 can be interpreted as a weakness of the discriminant analysis, mainly because of the lower accuracy of the calibration data, which were estimated from a space-based method instead of a field investigation, the limitations of the low-resolution SAR data and the heteroscedasticity between ∆γ and H. to logistic regression are given in Table 2. Figure 7 shows the damage proxy map of Sarpole-Zahab with scatter plots of ∆ and H, in which the majority of the collapsed buildings were two-story buildings. This finding does not mean that the 3-to 6-story buildings were intact or slightly damaged. The scatterplot of the collapsed (red circles) buildings in Figure 7 indicates that almost all of the buildings with heights more than 20 m in Sarpole-Zahab collapsed. The higher number of misclassifications among collapsed buildings shown by the thick dashed line in Figure 7 can be interpreted as a weakness of the discriminant analysis, mainly because of the lower accuracy of the calibration data, which were estimated from a space-based method instead of a field investigation, the limitations of the low-resolution SAR data and the heteroscedasticity between ∆ and H.

Estimation of Debris
The definition of debris generally includes natural (e.g., trees) and man-made (e.g., cement blocks) materials generated by a disaster. Several SAR-based methods are available for debris handling and volume and height estimation. If airborne SAR data are available, very high-resolution stereo SAR data with different polarizations can be used to define the base, which is a reference point on the nondebris-covered ground surface [31]. If airborne SAR data are not available, very high-resolution spaceborne data provide an effective means of estimating the building footprint and height, and the modeling of backscattering patterns can contribute to estimating the extent of damage and the debris volume [32]. In most cases, very high-resolution SAR images, either airborne or spaceborne, are not available. Here, we focused only on debris consisting of man-made materials (i.e., buildings), for which indirect relationships with the building heights and the discriminant threshold values are considered. The volumes of buildings, construction materials and debris are defined by Hashemi and Alesheikh [33] as follows: where Vb i and Ab i are the volume and area of building i, respectively, which are extracted from 3D vector data, H i is the height of building i, Vcm i and Vd are the volume of the construction materials of building i and the total debris volume, respectively, and r i is the percentage of damage to the building. In this paper, r i is unknown. As discussed above, the SSSC method does not explain the degree of damage because some buildings might be deformed but remain standing. However, we were able to calculate the discriminant threshold value based on the logistic method to separate the range of probability of collapsed buildings and force the buildings to be considered as either "collapsed" or "uncollapsed." This approach will not produce accurate debris estimations because assuming that debris will only be produced from G4 and G5 buildings, we should consider that according to the EMS-98 code, G5 buildings will have more debris and waste materials. Thus, we assume that the percentage of damage in G4 buildings is 50% less than that in G5 buildings. We then divide the buildings with values higher than 0.07 into two groups, which are likely to be similar to the G4 and G5 classes, using the "natural breaks" optimization. The method seeks to classify the buildings by minimizing the average deviation from the class mean while maximizing the deviation from the means of the other class, which enables it to decrease the variance within classes and maximize the variance between classes [34]. This is a reliable classification method when field survey results are incomplete or lacking. The break value in the "collapsed" class (range from 0.07 to 4.82) is 1.37. Values larger than the break value are likely to be in group G5, whereas values between 0.07 and 1.37 are in group G4 (Figure 8). The amount of debris is based on the initial results of the earthquake between two latter coherence maps; therefore, post-earthquake activities and the effects of human activities, such as the removal of debris related to shear walls of uncollapsed buildings, are not considered. Most of the debris is produced from heavily damaged buildings (G4) that have significant failures of their walls while their frames remain standing. Of the 799 suspicious buildings identified by the logistic and natural breaks methods, 652 are in group G4, and 147 are in group G5. The amount of debris for the G4 group is 95,454 m 3 , which is approximately three times larger than that of the G5 buildings. This finding implies that the total failure of buildings is not the main cause of debris in the town. The average heights of the G4 and G5 buildings are 6.2 m and 5.8 m, respectively, which indicates that in both damage categories, 1-and 2-story buildings are the most vulnerable structures. Details of the EMS-98 damage grades and the estimated volumes of debris in Sarpole-Zahab are given in Table 3. maximize the variance between classes [34]. This is a reliable classification method when field survey results are incomplete or lacking. The break value in the "collapsed" class (range from 0.07 to 4.82) is 1.37. Values larger than the break value are likely to be in group G5, whereas values between 0.07 and 1.37 are in group G4 (Figure 8). The amount of debris is based on the initial results of the earthquake between two latter coherence maps; therefore, post-earthquake activities and the effects of human activities, such as the removal of debris related to shear walls of uncollapsed buildings, are not considered. Most of the debris is produced from heavily damaged buildings (G4) that have significant failures of their walls while their frames remain standing. Of the 799 suspicious buildings identified by the logistic and natural breaks methods, 652 are in group G4, and 147 are in group G5. The amount of debris for the G4 group is 95,454 m 3 , which is approximately three times larger than that of the G5 buildings. This finding implies that the total failure of buildings is not the main cause of debris in the town. The average heights of the G4 and G5 buildings are 6.2 m and 5.8 m, respectively, which indicates that in both damage categories, 1-and 2-story buildings are the most vulnerable structures. Details of the EMS-98 damage grades and the estimated volumes of debris in Sarpole-Zahab are given in Table 3.

Discussion
The mean coherence over the town with an area approximately 7 km 2 decreases by up to 0.1 of the coherence (~10% of the coherence) due to the rooftop deformations of old buildings after heavy precipitation events. Washaya et al. also showed that a coherence loss of approximately 38% occurred immediately before the earthquake [10]. Using LSSC, the physical condition of each individual building, especially old and weak buildings, can be monitored if they are influenced by long-term stationary changes such as heavy rain. Additionally, the mean coherence decay caused by the earthquake was 0.24, which implies that stationary and nonstationary events can be easily distinguished with continuous urban monitoring using this method. The TC model is used to estimate the sensitivity of the mean coherence during the spring-to-winter and summer-to-autumn periods and to predict the mean coherence values in the absence of an earthquake. The model fits the data well, but it is not guaranteed to be reliable if the temporal parameter in Equation (6) is large. Accordingly, the impacts of different temporal baselines must be quantified in future studies. Because the relationship between the building height and the coherence is not sufficiently strong, the overall accuracy of the outputs of the logistic discriminant analysis is low. In addition, the number of misclassified buildings in the collapsed category was higher than that of the uncollapsed category. This was probably caused by two reasons: (1) because of the coarse resolution of the Sentinel-1 dataset, the pixel size is larger than the size of the buildings, so most buildings near collapsed buildings cannot be distinguished; and (2) for large buildings that cover more than one coherence pixel, we used only the pixel value that covers the centroid of the buildings, which increases the uncertainty of the building classification. Different pixel averaging methodologies or different window sizes may reduce the number of misclassified buildings and increase the overall accuracy.
Sarpole-Zahab was not sufficiently large to demonstrate all the capabilities of repeated SAR observations, but the model applied in Sarpole-Zahab can be used in other cities around the world with their relative coherence values (LSSC), relative urban footprints extracted by different methods and additional information about the cities in an online GIS geodatabase. Because the availability of digital surface models (DSMs) for the extraction of building heights is limited, the priority should be given to seismically vulnerable cities with extensive geodatabases. Using the large amounts of SAR data provided by Sentinel-1, new portals and tools with second-order SAR products (e.g., SAR coherence maps) have already been established [35,36]. Thus, for further progress in Earth observations, the model proposed in this paper can also be used to create new portals with third-order SAR products, such as RGB coherence maps, and for monitoring industrial cities before and during natural or anthropogenic disasters. Although vulnerability is a complex concept that consists of variable components, the development of a mean coherence index for building damage in each city or area with a specified boundary can provide a better and faster understanding of large changes, such as damage caused by earthquakes.

Conclusions
This paper presented the application of a sequential interferometric pair for the assessment of Sarpole-Zahab before and during the 2017 earthquake using nearly two years of Sentinel-1 data. Typically, interferograms with larger temporal baselines generate lower coherence. This effect is reduced by selecting regular 12-day intervals followed by the application of the sequential coherence method for the constant screening of Sarpole-Zahab, whereas seasonal coherence changes are explored during the rainy and snowy seasons.
The binomial logistic regression method was used to discriminate two types of buildings. This approach provides a more realistic method when two groups of buildings (i.e., collapsed and uncollapsed) are present. Because of the inherent limitations of SAR imagery for building damage assessments, it was difficult to consider more sophisticated logistic regressions (i.e., multinomial logistic regression), which apply to situations that can have three or more possible outcomes. In addition, we combined the results of the logistic regression analysis and the natural breaks classification to estimate the amount of debris in G4 and G5 buildings, assuming that the amount of debris is based on the immediate results of the SSSC method; therefore, the debris generated after the event cannot be directly addressed with this method. Funding: This study was funded by the Japan Society for the Promotion of Science (JSPS) Grants-in-Aid for scientific research (KAKENHI) number 16F16380 and 17H02050.