A Preliminary Damage Assessment Using Dual Path Synthetic Aperture Radar Analysis for the M 6.4 Petrinja Earthquake (2020), Croatia

: On 29 December 2020, an earthquake with a magnitude of M 6.4 hit the central part of Croatia. The earthquake resulted in casualties and damaged buildings in the town of Petrinja (~6 km away from the epicenter) and surrounding areas. This study aims to characterize ground displacement and to estimate the location of damaged areas following the Petrinja earthquake using six synthetic aperture radar (SAR) images (C-band) acquired from both ascending and descending orbits of the Sentinel-1 mission. Phase information from both the ascending (Sentinel-1A) and descending (Sentinel-1B) datasets, acquired from SAR interferometry (InSAR), is used for estimation of ground displacement. For damage mapping, we use histogram information along with the RGB method to visualize the affected areas. In sparsely damaged areas, we also propose a method based on multivariate alteration detection (MAD) and naive Bayes (NB), in which pre-seismic and co-seismic coherence maps and geocoded intensity maps are the main independent variables, together with elevation and displacement maps. For training, approximately 70% of the data are employed and the rest of the data are used for validation. The results show that, despite the limitations of C-band SAR images in densely vegetated areas, the overall accuracy of MAD+NB is ~68% compared with the results from the Copernicus Emergency Management Service (CEMS).


Introduction
Earthquakes are among the most destructive natural disasters on Earth. They are mainly unpredictable in terms of time and location. The most highlighted point about earthquakes is that, even though they only shake the Earth for a few seconds, the consequences left behind potentially last for a long time. Human loss, direct damage to buildings and infrastructures, and indirect damages such as fires, liquefaction, and landslides are just a few examples of the impacts of a moderate to severe earthquake [1][2][3][4][5][6][7]. People located near seismically active areas are one of the most vulnerable elements during catastrophic earthquakes, and disaster prevention plans implemented by organizations and construction companies can increase the resiliency of a society. Resilience includes both pre-disaster and post-disaster actions. Although pre-disaster plans may reduce the risk of seismic exposure, post-disaster actions also play an important role [8][9][10][11]. Remote sensing technology is a powerful tool not only to monitor environmental changes before a disaster but also to facilitate disaster responses in order to increase post-disaster resiliency. Remote sensing technology can provide various solutions and answers in disaster situations. For example, immediately after a severe earthquake, the main question is "where should the search and rescue teams go first?". Along with innovative damage assessment techniques in terms of loss/consequences [12,13], remote sensing techniques also can provide efficient damage displaced in the towns of Petrinja and Glina [35,36]. Figure 1 shows the epicenter of the earthquake near the towns of Glina and Petrinja; the mechanism of the movement was a right lateral strike-slip. Most of the heavily damaged buildings were reported from Petrinja, Glina, and Sisak; 2902 buildings were possibly damaged or destroyed out of a total of 40,000 buildings, based on the Copernicus report and OpenStreetMap [37]. This event was initiated by a foreshock in the same region, with a magnitude of 5.2 just one day before the main shock. After the main shock, about 500 aftershocks with a NW-SE trend (with a 40 km radius) were recorded in the study areas until 21 January 2021, in which the range of the aftershocks varies from M 1 to M 4.5. Most of the aftershocks occurred northwest of the main shock (Figure 1a). The Pokupsko fault (purple line in Figure 1), 75 km in length in the NW-SE direction, is the causative fault, which is consistent with the focal mechanism of the earthquake obtained from the United States Geological Survey (USGS). The felt intensities on the European macroseismic scale MMI (modified Mercalli intensity) in Petrinja and Glina were 7.5 and 7, respectively [35]. Peak ground acceleration (PGA) and peak ground velocity (PGV) in Petrinja are 73.66% g and 49.75 cm/s, respectively. However, the PGA and PGV values near the epicenter reached 121.19% g and 89.13 cm/s, respectively. The earthquake occurred in a hilly topographic location, and its depth was 13.5 km, which is rather shallow. The last earthquake in the Pokupsko-Petrinja area occurred in 1996 (M 6) in Ston; a historical event took place in 1909, which is known as the Pokupsko earthquake, with a seismic depth of 14 km [38]. PAGER (Prompt Assessment of Global Earthquake for Response), as an initial indicator of economic and human loss, reported that the probability of human loss of between 10 and 100 people is 35% while the estimated economic losses are 0% to 3% of the total gross domestic product (GDP) of Croatia.
In this study, we used two types of SAR datasets from both ascending (track 146) and descending (track 124) orbits of C-SAR (wavelength = 5.6 cm) satellites-Sentinel-1A and Sentinel-1B. The revisit cycles of Sentinel-1A and Sentinel-1B are each 12 days, but if we consider both of them together, the revisit cycle can be reduced to 6 days. Here, we aim to use SLC (single look complex) images with IW (interferometric wide) swath, which contain altitude and orbit information of the SAR images in the slant range. The acquired images in both ascending and descending orbits are dual-polarized (VV and VH). The incidence angle of the ascending master image is 39.473 • , and the absolute differential incidence angles for the master image with two slave images-6 December 2020 and 30 December 2020-are 0 • and 0.001 • , respectively. For the descending dataset, the incidence angle of the descending master image is 39.371 • and the absolute differential incidence angles for the master image with two slave images-11 December 2020 and 4 January 2021-are 0.001 • and 0.004 • , respectively. The temporal baseline of the acquired images from both the ascending and descending orbits is 12 days; the relative post-earthquake image from the ascending orbit was taken 1 day after the earthquake, while the post-earthquake descending image was taken 5 days after the earthquake. It must be noted that, due to the inherent limitations of the C-band in vegetated areas and temporal gaps, the correlation of the pre-seismic and co-seismic SAR pairs might decrease, especially for the descending dataset. A larger temporal gap might deteriorate the quality of damage mapping, since human activities such as establishing temporal shelters and tents increase in the affected area. Detailed SAR footprints and parameters are given in Table 1 and Figure 1. Table 1. Detailed information of C-SAR sentinel-1 images used in this work. "B", "T", "m", "A", and "D" denote spatial baseline, temporal baseline, master image, ascending orbit, and descending orbit, respectively.  Interferometric analysis was performed using two SAR images (12-day) for both ascending and descending datasets. Interferometric analysis is used to calculate how much the ground has displaced in the line-of-sight (LOS) due to the earthquake. As shown in Figure 2, the co-seismic interferograms with a "butterfly" shape for ascending track 146) T 146) and descending track 124 (T 124) confirm the characteristics of strike-slip movements. The movement is extended approximately 25 km in the E-W direction. The quality of the interferograms of T 146 and T 124 is different because of temporal changes with respect to the main shock, coherence, and atmospheric noises. However, assuming that the M 6.5 event is the only significant ground movement during the two acquisitions, we estimated that the ground across the Pokupsko fault moved approximately 40 cm in the look direction of the satellite (Figure 2a,b). Despite the inherent limitations of the C-band Interferometric analysis was performed using two SAR images (12-day) for both ascending and descending datasets. Interferometric analysis is used to calculate how much the ground has displaced in the line-of-sight (LOS) due to the earthquake. As shown in Figure 2, the co-seismic interferograms with a "butterfly" shape for ascending track 146) T 146) and descending track 124 (T 124) confirm the characteristics of strike-slip movements. The movement is extended approximately 25 km in the E-W direction. The quality of the interferograms of T 146 and T 124 is different because of temporal changes with respect to the main shock, coherence, and atmospheric noises. However, assuming that the M 6.5 event is the only significant ground movement during the two acquisitions, we estimated that the ground across the Pokupsko fault moved approximately 40 cm in the look direction of the satellite (Figure 2a,b). Despite the inherent limitations of the C-band (wavelength 5.6 cm) SAR images in green areas, the interferometric fringes of the co-seismic pairs are visible in both the ascending and descending datasets (~4 fringes). (wavelength ~5.6 cm) SAR images in green areas, the interferometric fringes of the coseismic pairs are visible in both the ascending and descending datasets (~4 fringes).

Methods
This study focuses on two aspects: a) visualization of the study area; b) classification of the affected buildings. We propose a framework for damage characterization after the M 6.4 earthquake using the analyzed Sentinel-1A and Sentinel-1B intensity images and RGB color composites of differential SAR coherence values as the first step and then a classification method based on MAD and NB techniques ( Figure 3). The visualization part is a semi-automated approach for fast and effective illumination of different aspects of the changes in the earthquake-stricken area based on backscattering coefficient (σ ) coherence values [25,26]. Here, information on collapsed and damaged buildings was obtained from the Copernicus emergency management service (EMS) for Petrinja and Glina, which was combined with the building information from OpenStreetMap (OSM) as an input for damage evaluation.

Methods
This study focuses on two aspects: a) visualization of the study area; b) classification of the affected buildings. We propose a framework for damage characterization after the M 6.4 earthquake using the analyzed Sentinel-1A and Sentinel-1B intensity images and RGB color composites of differential SAR coherence values as the first step and then a classification method based on MAD and NB techniques ( Figure 3). The visualization part is a semi-automated approach for fast and effective illumination of different aspects of the changes in the earthquake-stricken area based on backscattering coefficient (σ 0 ) coherence values [25,26]. Here, information on collapsed and damaged buildings was obtained from the Copernicus emergency management service (EMS) for Petrinja and Glina, which was combined with the building information from OpenStreetMap (OSM) as an input for damage evaluation. Remote Sens. 2021, 13, x 6 of 19

Visualization
The visualization of the study area is based on auxiliary information (e.g., OSM), SAR backscattering, and SAR coherence.

SAR Backscattering
Sentinel-1 SAR images are single-look complex (SLC) images gathered by a side looking sensor, which represents the data in real and imaginary portions of the SLC data in order to preserve amplitude and phase information. For change visualization, two SAR images (or more) must have the same size and pixel number. For exact pixel matching, one of these images must be selected as a master image and the other images are the slave images [39][40][41]. Since the SAR geometry is different, the range coordinate is the slant range and the azimuth coordinate is satellite movement direction. These coordinates are not consistent with geographic coordinates. SAR images are described by considerable distortions in range direction (mainly because of topographic changes). Thus, image calibration and coordinate system conversion from the SAR coordinate system to the geographic coordinate system is necessary. The coordinate transformation is possible with a backward technique that uses a digital elevation model (DEM) to convert the positions. The Doppler frequencies, pulse transmission time of the SAR image, the position, and the velocity vectors of both sensor and backscatter elements are used for calibration as follows: where R s is the slant range, S and P are the positions of the satellite and backscatter, v s and v p are the velocity of the satellite and backscatter element, c is the speed of light, f 0 is the frequency of the carrier, and f D is the processed Doppler frequency. Once the processed Doppler frequency is achieved, during the geocoding processing, the results can be converted into the required reference systems. Finally, the derived SAR brightness values (β 0 ) can be converted to the backscattering values using a simple equation between β 0 and local slope (local incidence angle) as follows: where σ 0 is the radar backscattering coefficient, β 0 is the brightness value, and α is the local incidence angle in the SAR imaging direction. The σ 0 can be expressed as a calibrated linear value or decibel value. Accordingly, in the rest of the text, the σ 0 is σ 0 dB . In this study, we express the outcomes in decibel values as follows: The level of noise (speckle) in the produced maps was reduced with the noise-removal technique called Refined Lee Filter (RLF), which is based on the standard Lee filter [38,39]. In the RLF technique, the K-nearest neighbor (KNN) algorithm adjusts the number of pixels used for filtering in each moving window [42].

RGB Coherence Visualization and InSAR Deformation
The SAR coherence visualization is based on the estimation of interferometric phase correlation of two pairs (pre-seismic and co-seismic). The phase cross-correlation of two SAR images in different times (e.g., time 1 and time 2) for the same pixel is called "SAR coherence" in which high coherence values indicate that the features have not been changed during that period [43,44]. In contrast, low coherence values represent considerable changes during that period. Since the last decade, coherence methods have been among the preferred damage detection techniques because of the availability of the repeated SAR imagery systems and coherence sensitivity to the changes. Fielding et al. [44] showed the feasibility Remote Sens. 2021, 13, 2267 8 of 18 of the differential SAR coherence technique (pre-seismic coherence and co-seismic coherence) for building damage detection after the Bam earthquake (2003) using three Envisat ASAR images. Yun et al. [45] showed damage concepts after the Gorkha earthquake (2015) using coherence analysis as a "damage proxy map" with ALOS-2 and Cosmo-SkyMed SAR images. More recently, Karimzadeh et al. [26] improved change monitoring methods using coherence analysis of Sentinel-1 SAR images as a "Sequential SAR coherence" method for continuous monitoring of changes before and after the Kermanshah earthquake (2017). Overall, for large areas, the detection of damage using coherence is a superior method, but at least three SAR images from the same track with the same geometry are needed. The coherence method provides additional information. If the temporal gap is short, the results of the changes are highly associated with the event (e.g., earthquake), but if the temporal baseline of the two images is considered large, other changes (e.g., vegetation growth) might be reflected in the final coherence map. The interferometric phase correlation measures the coherency of the complex phase values of two SAR acquisitions. It varies between 0 and 1 and is a measure of the quality of the generated interferogram. The coherence (γ) of two zero-mean circular Gaussian variables, a and b, can be defined as follows: where a represents the relative complex values of the master image and b represents the relative complex values of the slave image for interferometric analysis, * represents the complex conjugates of the images, and E is the expectation operator. To achieve the pure phase of changes, the topographic phase is removed using a 30 m digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) and atmospheric effect noise was also reduced using atmospheric models [46][47][48]. Usually, the coherence of vegetated areas is low because of the growth of leaves and random movements of the foliage of trees due to environmental effects such as wind. In contrast, urban areas are constant and movement due to short-term effects (e.g., wind) is not possible. Thus, we always expect higher coherence values from urban areas. Here, we applied multitemporal coherence analysis using both ascending and descending datasets (VV polarization) for visualization of the affected area. We produced a differential coherence map using one pre-seismic coherence and one co-seismic coherence map for each dataset: For RGB visualizations, we defined three parameters. The first parameter is normalized "forward change", which means that the changes are related to an event. The second parameter is normalized "reverse change", which means that the changes are related to post-event activities such as human activities and vegetation growth. The third parameter is "constant", which means that there are no changes. These concepts are defined as follows [25]: Once the three parameters are calculated from Equations (7)-(9), they are assigned to "red", "green", and "blue" channels, respectively.
To extract the displacement map, first, the generated interferograms are filtered by an adaptive method based on the local scene coherence and Goldstein filter [49]. The technique filters incoherent areas more than coherent areas to increase the visibility of the deformation fringes and reduces the level of noise from temporal and spatial baselines. As Remote Sens. 2021, 13, 2267 9 of 18 shown in Figure 2, the deformation signals are repetitive color bars, which indicate that the phase of the interferograms was repeated from 0 to 2π several times. When the phase exceeds 2π, the phase starts again from 0 to form another fringe. The phase unwrapping process resolves the ambiguity of interferograms for 2π. However, the errors during the phase unwrapping process must be treated. Here, the MCF (minimum cost flow) algorithm is used to improve the quality of the unwrapped displacement maps, in which all pixels with coherence values lower than the "Unwrapping Coherence Threshold" are ignored [50]. The obtained phase values are converted to the displacement values. Each 2π color bar (fringe) corresponds to half of the wavelength of the satellite in the range direction, which should be calibrated and geocoded at the final stage.

Classification
For damage or post-disaster classifications, different studies have been carried out . As an alternative, MAD and NB techniques are deployed in this study to classify collapsed/damaged buildings from intact/moderately damaged buildings. Here, a differential change detection technique is carried out for two stacked (3 bands) images (preseismic and co-seismic) that contain 4 bands (differential backscattering map, pre-and post-seismic maps, and pre-and co-seismic coherence maps) for both ascending and descending datasets using the MAD algorithm [51]. When analyzing the changes, after partial normalization, pixel values with small differences represent little to no change and areas with large changes have large differential pixel values. The stacked images are assumed as follows: where k is number of layers in the stacked images (X and Y) and the differential maps can be calculated as follows: The MAD algorithm can produce a set of change maps from Equation (11), where k is the maximum number of bands in the first (pre-seismic) and second (co-seismic) input images. The MAD change maps have two aspects: (1) there are differences between the pair of linear combinations of bands from the first stacked image and bands from the second stacked image, chosen to maximize the correlation; (2) each MAD change map is orthogonal to the others. Thus, the MAD algorithm is a statistical method that can use different modalities and bands to produce a single multiband change map, sorted by increasing correlation. The lowest MAD correlation together with a DEM and displacement map are chosen for NB classification. NB is a probabilistic machine learning model based on Bayes' theorem. The NB model considers the presence of an object/feature in a class as unrelated to the presence of any other object. For example, a building might be considered in a "collapsed" class if it has large differential coherence, a large differential backscattering coefficient, and a large displacement value (close to the epicenter). Even if these parameters depend on each other or upon the existence of the other buildings, all of these properties independently contribute to the probability that this building is a "collapsed" building and that is why it is called "naive". Its naïve forms can be formulated as follows: where the probability of A occurring can be found, given that B has occurred. Here, B is the evidence, P(B) is the prior probability of the prediction, A is the hypothesis, and P(A) is the prior probability of the class. Despite the simple structure of the NB model, it is known to outperform some other advanced classification methods [52]. Here, two classes of buildings are defined as "not collapsed" with label "0" and as "damaged/collapsed" with label "1". Since the number of buildings in class "0" and "1" is not equal (most of the buildings are in class 0), an equal number of buildings for each category should be selected for training and validation. The model is trained by assuming that the distribution of the quantitative values of the variables (e.g., coherence values, elevation, and displacements) follows a normal distribution, which does not allow zero probabilities for the classes.

Visualization Results
Figure 4a-f shows the refined backscattering SAR images of the study area in decibel values. Darker areas show that the roughness of the surface was likely very low because of the mirror-like reflection mechanism that causes the lower backscattering values. Figure  4g shows the damage distribution in Petrinja based on auxiliary information gathered by the Copernicus emergency management service (EMS). However, the damage concept in the post-event backscattering maps is not clear by visual inspection. This is mainly because of two reasons: (a) sparse distribution of the damaged or collapsed buildings in different locations of the study area and (b) the spatial resolution of the Sentinel-1 images. Although the changes in the built-up areas are not visible in SAR backscattering maps, the histograms are helpful for the interpretation of seismic changes. A histogram of a series of SAR images for an event presents useful information such as backscattering coefficients and coherence in different polarization modes. Tabulating the frequency of occurrence of each σ 0 value within the images provides statistical information, as shown in Figure  4h,i in their symmetrical shapes. If a large number of pixels in the map have the same σ 0 values, the histogram information might not be the best way to recognize the content of the SAR intensity data. As shown in Figure 4h,i, for both ascending and descending datasets (VV polarization), the pre-event backscattering histograms (red and green) are similar in shape and height because, before the event, there was no significant change (e.g., flood, rain, earthquake, etc.) in the built-up areas. In contrast to the pre-event images, the histograms of the post-event images (blue) show that changes in the urban areas due to earthquakes can be reflected in the height of the post-event histograms. The mean values as univariate statistics for the two pre-event images of the ascending orbit are −5.2 (dB) and −5.5 (dB), respectively, while the mean value for the post-event image is −4.9 (dB). The difference in mean value between the pre-event images is almost two times smaller than the difference in mean value between one pre-event and post-event image. For the descending dataset, the situation is the same. Figure 5a,b shows pre-seismic and co-seismic coherence maps for the Sentinel-1 ascending dataset, and Figure 5c,d show pre-seismic and co-seismic coherence maps for the Sentinel-1 descending dataset (differential coherence for ascending and descending datasets are available in Figure S1 in the electronic supplement to this article). In the preseismic coherence maps, the brighter areas show higher stability, likely in the urban areas, while the darker areas are associated with lower coherence values due to vegetation and other reasons. In the co-seismic coherence maps, the stability shows some perturbations. Although the higher coherence values are still related to the urban area, due to earthquake damage, heavy rains, etc., the level of coherency is decreased considerably compared to the pre-seismic maps. A coherence decay due to the earthquake is expressed above as γ di f f , which is an initial indicator of changes that can be used for damage visualization. For pre-seismic coherence maps of ascending and descending datasets in the urban areas (red polygons in Figure 5a-d), the mean coherence values are 0.55 and 0.59, respectively. The standard deviation of the pre-seismic coherence maps for both ascending and descending datasets is 0.14. For co-seismic coherence maps of ascending and descending datasets, the mean coherence values are 0.49 and 0.45, respectively. The standard deviation of the pre-seismic coherence maps for both ascending and descending datasets is 0.16. The mean coherence decay (γ di f f ) for ascending and descending datasets is 0.06 and 0.1, respectively. The higher mean coherence decay for the descending dataset does not mean that the damage visualization in descending γ di f f will perform better than ascending γ di f f . It could be related with the post-seismic images of the descending dataset (2021.01.04), which were obtained 5 days after the main event. It is likely that the heavy rains and other human activities during these 5 days will deteriorate the coherency of the pure damages. As shown in Figure 4f, there are some water accumulated areas in the northeast of Petrinja, which confirms the precipitation in the post-seismic descending image. Figure 5e,f show histograms of the pre-seismic and co-seismic ascending and descending datasets, respectively. The shape and amplitude of the pre-seismic histograms for ascending and descending pairs (red and blue) are very similar, but for co-seismic histograms of ascending and descending pairs (green and gray), the black double arrows show that the peak coherence retreat in the descending dataset is considerably higher than that of the ascending dataset. Labeling the pre-seismic and co-seismic coherence maps of ascending and descending datasets as band 1, band 2, band 3, and band 4, details of the statistical multivariate analysis are given in Table 2.   Figure 5a,b shows pre-seismic and co-seismic coherence maps for the Sentinel-1 ascending dataset, and Figure 5c,d show pre-seismic and co-seismic coherence maps for the Sentinel-1 descending dataset (differential coherence for ascending and descending datasets are available in Figure S1 in the electronic supplement to this article). In the preseismic coherence maps, the brighter areas show higher stability, likely in the urban areas, while the darker areas are associated with lower coherence values due to vegetation and other reasons. In the co-seismic coherence maps, the stability shows some perturbations. Although the higher coherence values are still related to the urban area, due to earthquake damage, heavy rains, etc., the level of coherency is decreased considerably compared to the pre-seismic maps. A coherence decay due to the earthquake is expressed above as , which is an initial indicator of changes that can be used for damage visualization. For pre-seismic coherence maps of ascending and descending datasets in the urban areas (red polygons in Figure 5a-d), the mean coherence values are 0.55 and 0.59, respectively. The standard deviation of the pre-seismic coherence maps for both ascending and de- scending datasets, respectively. The shape and amplitude of the pre-seismic histograms for ascending and descending pairs (red and blue) are very similar, but for co-seismic histograms of ascending and descending pairs (green and gray), the black double arrows show that the peak coherence retreat in the descending dataset is considerably higher than that of the ascending dataset. Labeling the pre-seismic and co-seismic coherence maps of ascending and descending datasets as band 1, band 2, band 3, and band 4, details of the statistical multivariate analysis are given in Table 2.    Figure 6a shows the damage reference map, which is a combination of the Copernicus damage map and building inventories of OSM. Figure 6b,c are RGB damage of Petrinja and Glina from Equations (7)-(9) for ascending and descending orbits. For the descending RGB map (Figure 6c), because of the rather larger temporal gap between the event and post-event image, the red pixels are more than the RGB damage map of the ascending dataset (Figure 6b). The green areas mainly represent the vegetation growth/change, but since the built-up areas are masked by the OSM map, the number of green pixels in Figure  6b,c is very low; blue (or purple) pixels are abundant in the study area, which implies that most of the buildings are not affected by the earthquake. As shown in Figure 7a,b, mean, maximum, and minimum InSAR displacement (cm) for the entire study area from T146 and T124 are −0.006, 0.43, −0.35, 0.007, 0.40, and −0.19, respectively. A cross comparison of InSAR displacement deduced from T146 and T126 is shown in Figure 7c. The mean displacement value of intact (or slightly damaged) buildings is insignificant; for T146 and T124, it is −0.07 (cm) and +0.07 (cm), respectively, while the mean displacement values of collapsed (or damaged) buildings for T146 and T124 are −0.17 (cm) +0.18 (cm), respectively. Mean railway and road displacements for T146 and T124 are −0.036, −0.013, 0.024, and 0.009, respectively. Further statistical information are given in Table 3.

Classification Results
Once the MAD map is extracted from multivariate analysis (Figure 8a-c), the NB algorithm classifies the features in a probabilistic way. Thus, an assumed function can be expressed by the independent variables to explain the outcome in a functional state. Here, we have three independent variables (i.e., displacement map, digital elevation map, and MAD map). To create the abovementioned structure, first, the independent variables of each building were extracted from geographic coordinates. Second, the dependent variables are assigned for each building. Dependent variables are assigned in a binary mode in which the "intact" or "slightly damaged" buildings are considered as "0" and "collapsed" or "damaged" buildings are considered as "1". expressed by the independent variables to explain the outcome in a functional state. Here, we have three independent variables (i.e., displacement map, digital elevation map, and MAD map). To create the abovementioned structure, first, the independent variables of each building were extracted from geographic coordinates. Second, the dependent variables are assigned for each building. Dependent variables are assigned in a binary mode in which the "intact" or "slightly damaged" buildings are considered as "0" and "collapsed" or "damaged" buildings are considered as "1". The total number of buildings used in this study is 47,982, in which 44,603 buildings are categorized as "intact" or "slightly damaged" buildings, and 3379 buildings are categorized as "collapsed" or "damaged" buildings. However, because of the coverage of the The total number of buildings used in this study is 47,982, in which 44,603 buildings are categorized as "intact" or "slightly damaged" buildings, and 3379 buildings are categorized as "collapsed" or "damaged" buildings. However, because of the coverage of the study area and uncertainties of correct identification of the "slightly damaged" or "damaged" categories, we randomly selected 384 buildings from the "intact" and "collapsed" categories. The training set is balanced, in which 50% of the selected buildings are labeled "0" and 50% of the selected buildings are labeled as "1". The total number of the training set is 70% of the total collapsed and intact buildings in Petrinja and Glina. The overall (OA) accuracy of the building classification for the proposed method using MAD and NB is 67%, while for the collapsed buildings, the OA is 61%. We also applied only NB and support vector machine (SVM) algorithms, with only differential coherence values of ascending and descending datasets (without the elevation and displacement values) as input variables. The OA for plain NB and SVM was 64% and 63%, respectively. The OA of the collapsed category for the algorithms was 56% and 57.5%, respectively.

Discussion
In this study, deformation pattern and damage maps are presented for the Petrinja earthquake. Since the damaged area is not concentrated in a specific region, as with previous earthquakes studied in [25,26], NB and MAD techniques were used to explain the classification of the damaged buildings. MAD and NB are selected for the scattered damages case study based on empirical tests with four bands. Other indices such as principal component analysis (PCA) and optimum index factor (OIF) were also tested before the classification. The MAD was suitable for our purposes since it is a simpler technique that uses a set of several change maps to generate change maps with differences between pairs of linear combinations of bands to maximize the correlation. Classification algorithms such as SVM were also tested, along with NB, for the study area.

Conclusions and Future Directions
We conclude that SAR remote sensing products including SAR histograms and RGB coherence combination are effective for rapid damage mapping. The preliminary information presented in this study for the affected roads, along with visualization maps of the affected buildings, can be used for fast and appropriate disaster response. The major outcomes and highlights of damage visualization and classification are as follows: 1.
Detection of sparsely damaged areas is generally challenging using single path SAR imagery because the extent of damage in a specific search area is small. The proposed method utilized Sentinel-1 ascending and descending datasets together to optimize the damage classification results. We conclude that the simultaneous use of ascending and descending datasets makes the proposed method applicable to sparsely damaged areas after an earthquake occurs.

2.
The accuracy evaluation showed that the combined method together with MAD and NB are effective for sparse damage detection. Comparing the classification results of the proposed method with SVM results, we also conclude that although SVM showed good performance in densely damaged areas (e.g., the Sarpol-Zahab earthquake) [29], for sparsely damaged case studies such as the Petrinja earthquake, NB is more suitable because it is not sensitive to irrelevant pixel values in sparsely damaged areas. NB is also highly scalable, with a number of variables that produce better results with multivariate datasets. However, NB is not superior for damage classification in densely damaged areas, because it assumes that all features are independent, which is a "naive" assumption in real-world implementation of damage mapping.

3.
This study presented a binary (0-1) damage detection using dual path SAR imagery which is not fully compatible with physical damage scales such as European macroseismic scale 98 (EMS-98). In EMS-98 physical conditions of side walls, rooftops, etc., are taken into account to categorize the damage states into several classes [53]. However, due to limitations of both SAR and optical remote sensing techniques, damage expression is a bit different than those of EMS-98. In order to create a more meaningful connection between damage scales (e.g., EMS-98) and remote sensing techniques, very high-resolution images are necessary. Since the SAR imagery is inherently side-looking, the relationship between damaged buildings and backscattered signals can be explained more clearly if both ascending and descending orbits are available. For optical imagery, if nadir-looking imagery is not effective enough to the walls or other elements of buildings, pictometry might also help us to explain other damage grades which are related with walls of buildings.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13122267/s1, Figure S1: Differential coherence maps of the study area for both ascending (T 146) and descending (T 124) datasets. Yellow and red polygons are Petrinja and Glina, respectively.