Coherence Difference Analysis of Sentinel-1 SAR Interferogram to Identify Earthquake-Induced Disasters in Urban Areas

This study proposes a workflow that enables the accurate identification of earthquake-induced damage zones by using coherence image pairs of the Sentinel-1 satellite before and after an earthquake event. The workflow uses interferometric synthetic aperture radar (InSAR) processing to account for coherence variations between coseismic and preseismic image pairs. The coherence difference between two image pairs is useful information to detect specific disasters in a regional-scale area after an earthquake event. To remove background effects such as the atmospheric effect and ordinal surface changes, this study employs the two-step threshold method to develop the coseismic coherence difference (CCD) map for our analyses. Thirty-four Sentinel-1 images between January 2015 and February 2016 were collected to process 30 preseismic image pairs and two coseismic image pairs for assessing multiple types of disasters in Tainan City of southwestern Taiwan, where severe damages were observed after the Meinong earthquake event. The coseismic unwrapping phases were further calculated to estimate the surface displacement in east-west and vertical directions. Results in the CCD map agree well with the observations from post-earthquake field surveys. The workflow can accurately identify earthquake-induced land subsidence and surface displacements, even for areas with insufficient geological data or for areas that had been excluded from the liquefaction potential map. In addition, the CCD details the distribution of building damages and structure failures, which might be useful information for emergency actions applied to regional-scale problems. The conversion of 2D surface displacement reveals the complex behavior of geological activities during the earthquake. In the foothill area of Tainan City, the opposite surface displacements in local areas might be influenced by the axis activities of the Kuanmiao syncline.


Introduction
The Meinong earthquake struck southwestern Taiwan on 5 February 2016 at 07:57 Coordinated Universal Time (UTC).The hypocenter of the M L 6.4 earthquake was located on 22.92 • N and 120.54 • E (Meinong district, Kaohsiung City) with a 16.7-km focal depth.The focal mechanism of this event was a thrust mechanism with two nodal planes striking in the north-south and the northwest-southeast directions (see Figure 1).This has been the most devastating earthquake in southwestern inland Taiwan after the 2010 Jiasiang earthquake.Serious damages were observed in Tainan City, which is about 35 km west of the epicenter, and the earthquake caused about 117 deaths, 551 injuries, 10 collapsed buildings and 247 seriously damaged buildings [1][2][3][4][5].The Taiwan Earthquake Model team announced that this earthquake event did reduce the seismic hazard potential in southern Taiwan, but more than about 80% of the seismic hazard potential was not yet released [6].After the earthquake event, field surveys provided detailed records of surface deformations and building damages in Tainan City.The field surveys indicated that the urban area exhibited multiple types of disasters, including collapsed buildings, damaged buildings and infrastructures, liquefaction features, surface cracks and small-scale landslides [7][8][9].However, the limited resources, large scale and multiple types of disasters made the field surveys difficult for emergency response and reconstruction activities.Moreover, the high population density in Tainan City has limited the efficiency of obtaining a relatively complete disaster map after an earthquake.These limitations motivated the use of remote sensing technologies to improve the efficiency of disaster identifications.
The radar interferometry (InSAR, interferometric synthetic aperture radar) technique has been frequently used for studies relevant to land surface deformations in Taiwan.These implementations include coseismic deformation, interseismic deformation and land subsidence due to pumping up of groundwater [10][11][12][13][14][15][16][17][18][19][20][21].Huang et al. [22] processed radar images with the InSAR technique to obtain the coseismic displacement and discussed the relationship between surface deformation and fault movements during this earthquake.However, the relationship between regional deformation and disaster type was not clearly addressed.Figure 2 demonstrates the case in Tainan City and shows that the local fringes exist for both ascending and descending orbits.The coseismic interferograms of the Meinong earthquake present many local fringes, and the diameters vary form 1-3 km at Annan, Yongkang, Xinshi, Xinhua and Guanmiao districts in Tainan City.Table 1 lists the information of the coseismic image pair.Although large-scale deformation can be obtained by using coseismic InSAR results, it is challenging to identify the affected ranges and differentiate the tectonic-induced and the regional disaster-induced displacements after an earthquake event.Conventional InSAR uses two complex SAR images acquired at the same location, but at different times, and calculates their correlation.The process is recognized as coherence analysis.When an event changes the ground characteristics and makes the coherence value reduce, this process is called decorrelation [23][24][25][26][27].The analysis of coherence variation has been successfully applied in many fields of Earth sciences, including volcanology, landslide hazard and earthquake hazard [28][29][30].Among the earthquake hazards, the decorrelation characteristic is useful for detecting the surface deformation and cracks [30,31], distinguishing the collapsed and damaged buildings [31,32] and mapping the ranges of soil liquefaction [32][33][34][35].
Few studies have analyzed the coherence difference before and after a large earthquake event in Taiwan [36].The main limitations come from the low frequency of acquired radar images in Taiwan, the complex process of obtaining images and the expensive image price in the past.The Sentinel-1A/B satellites from the European Space Agency (ESA, Paris, France) can acquire images in the same place every six days with alternate screens, and the images are free to download from the official website (https://scihub.copernicus.eu).With the improvement of satellite technologies and the associated image processing, the efficient development of a disaster map has become possible.
This study aims to develop a workflow that integrates a series of imaging processes to obtain an earthquake-induced disaster map in the urban areas.The data in Tainan City during the Meinong earthquake were collected and used for the study.We analyzed the coherence of Sentinel-1 image pairs and estimated the surface displacements in east-west and vertical directions by calculating the unwrapping phase images.The results of the earthquake-induced disaster map were then compared with the field surveys and geophysical observations to assess the reliability of the disaster map.Additionally, the interactions between disaster types and regional deformation were systematically discussed.In Section 2 of this paper, we briefly introduce the study site and relevant studies of the Meinong earthquake, the basics of the InSAR technique and the details of the coherence analysis method.Section 3 shows the results of the background effect and the comparison between coseismic coherence difference and field surveys.Section 4 discusses the comparison of the coseismic coherence difference and other geophysical data such as the liquefaction potential map and peak ground acceleration data, as well as the 2D geophysical conversion of surface displacement.Section 5 provides the conclusions and recommendations for the study.

Study Area
The study area covered 15 districts of Tainan City and two active faults (see Figure 1).The Hsinhua fault is located in the Xinhua district, and the Houchiali fault crosses the Yongkang, East and Rende districts.The Yenshui River and Tsengwen River pass through rural and partial urban areas in Tainan City.These two rivers are the boundaries of many districts.
After the Meinong earthquake, Wen et al. [37] relocated the aftershock catalog and indicated three clusters of aftershocks.The first cluster of aftershocks near the epicenter was 10-20 km in depth, then the second cluster was 5-10 km in depth and 10 km north of the hypocenter.However, the third cluster of aftershocks was 25 km away in the west of the hypocenter (between southern Hsinhua fault and eastern Houchiali fault) and not associated with new activities of geologic structures.Their study suggested that a 20-30-km depth pre-existing normal fault beneath Tainan City was triggered by the Meinong earthquake.Moreover, the peak ground acceleration (PGA) shake map of the Meinong earthquake displayed that the area that suffered a ground acceleration higher than 400 gal was not at the main shock location, but was near the aftershocks' distribution in Tainan City [4].
The continuous GPS data show that the coseismic deformation in the Longci district (west side of Guanmiao district) had the largest vertical displacement of 100 mm and west horizontal displacement of 50 mm [38].Based on the seismic data, GPS data and InSAR, Huang et al. [22] presented a two-fault model, which included a 15-20-km depth main thrust fault and a 5-10-km depth shallower thrust fault.The results indicated that the shallower fault was triggered by the Meinong earthquake under a high-background stress level and high fluid pressure.Furthermore, Kuo-Chen et al. [39] deployed 36 temporary seismic stations immediately after the earthquake to obtain the shallow crustal velocity structure.They found that the distribution of low S-wave velocity at a 0-4-km depth corresponded to the surface deformation obtained from coseismic InSAR, interseismic GPS and leveling data.The complex structure and propagation mechanisms triggered by the earthquake have never been observed in inland southwestern Taiwan [5].
These previous investigations had observed that an unknown fault beneath Tainan City was triggered by Meinong earthquake and that the unknown fault may have caused higher PGA than 400 gal.Although there is no common consensus about the geologic structure of the unknown fault, the structure was under a high fluid pressure based on the observation of low S-wave velocity at a depth of 0-4 km.This observation also implied that the surface deformation would be mainly controlled by soft materials in this area.

Interferogram Process
The C-band radar satellite Sentinel-1 is the most advanced SAR mission of ESA to implement terrain observation with progressive scans (TOPS) in azimuth mode.The instrument provides four exclusive modes with difference resolutions and land coverage.The interferometric wide swath mode (IW) is a widely-used mode in most applications.The IW mode can obtain surface information with wide area and good resolution at one time.It supports three sub-swathes, and the total swath width is more than 250 km.Moreover, the single-look resolution can reach 5 m × 20 m (range × azimuth) with the 6-day repeat cycle, and it is free to download for scientific research purposes.Therefore, it is useful to monitor the surface variation with time in our study.The high-resolution (12.5 m) digital elevation model (DEM) data were obtained from the Alaska Satellite Facility Distributed Active Archive Center (ASF DAAC), and the data were reprocessed to be 15 m in resolution for removing the topographic effect (dataset: ASF DACC 2017, ALOS-1 PALSAR_Radiometric_Terrain_Corrected_high_res; includes Material c JAXA/METI 2008.10.5067/Z97HFCNKR6VA).
We used an open-source radar interferometry software, Generic Mapping Tools Synthetic Aperture Radar (GMTSAR), which was developed by Sandwell et al. [40], for processing interferograms.GMTSAR has five main processors: (1) a preprocessor to convert the radar image and orbital information into the general format; (2) an InSAR processor to co-register two radar images, remove the topographic phase and form the complex interferogram; (3) a post-processor, based on Generic Mapping Tools (GMT) [41], to filter the interferogram and construct interferometric phase, amplitude and coherence; (4) an unwrapping processor based on SNAPHU software [42] to convert the interferometric phase into line-of-sight (LOS) displacement; and (5) a geocoding processor to switch the radar image coordinates into geographic coordinates [43].
A SAR image contains amplitude and phase information.The amplitude is the strength of scattered radar signal, and the phase is the fraction of a complete wavelength cycle, which represents the distance between the satellite antenna and the ground targets.The general terms of the single look complex (SLC) image can be written as [44]: where C(x) is a complex value, A(x) is amplitude, φ(x) is phase, e is Euler's number of exponential function and i = √ −1 is an imaginary number.Because the SLC image is composed of a regular grid, x can be replaced by (ρ, a).
C(ρ, a) = A(ρ, a)e iφ(ρ,a) , where ρ is the range and a is the azimuth.When using two SLC images acquired at different times to produce a complex interferogram, this complex interferogram consists of backscatter amplitude and phase change between observations.We can write down the terms of the complex interferogram as: where C 1 and C 2 are master and slave SLC images, A 1 and A 2 are master and slave amplitudes and φ 1 and φ 2 are master and slave phases.The notation * is the complex conjugation, and R is the real part of the complex interferogram.I represents the imaginary part of the complex interferogram.The interferometric phase is extracted in the usual way and can be written as [44]: The interferometric coherence (γ) refers to the amplitude of the complex correlation coefficient between two SAR images and can be written as [28]: where angle brackets are the statistical expectation with a rectangular filter.Although spatial resolution decreases with the rectangular filter size (also called the window size) [32], the filter is designed to further reduce the difference in radar impulse response perceived by each satellite track from the same piece of ground [25].The single-look resolution of Sentinel-1 is 5 m in range and 20 m in azimuth.To improve the quality of the amplitude image, we first averaged 6 × 2 multi-looks for a resolution cell as 30 × 40 m; we then used a window size with 5 × 5 pixels to obtain a spatial averaging coherence.The corresponding spatial resolution of the image in this study was roughly 150 × 200 m.

Coherence Difference Induced by the Earthquake
The value γ = 1 indicates no change.When surface changes alter the exact complex backscatter and cause decorrelation, the value γ would reduce.There are three major factors that contribute to coherence in repeat-pass interferometry [23]: where γ N is the thermal noise, γ S is the spatial correlation and γ T is the temporal correlation.When two radar signals are acquired by two different antennas, receivers capture the same target at the same time, and the thermal noise might influence the γ N .However, when signals are acquired by the same antenna, the variations of thermal noise can be assumed to be the same and is a trivial extension to the γ N .The sources of spatial decorrelation are mainly influenced by the spatial baseline and the sensor flight attitude.Although the spatial baseline decorrelation can be derived by applying Fourier transform to the difference in look angles for two observations' increases, using smaller baselines is relatively more efficient to increase the γ S .Under such conditions, the decorrelation of sensor flight attitude changes can be made negligible by using near-repeat precise orbits.The temporal decorrelation reflects the physical change in the surface and contains the motion of backscatters within a resolution element between observations [23,28].
When an earthquake triggers disasters of soil liquefaction and building structure damages, the paths of the radar wave for the same location are different before and after the earthquake.Therefore, scatterers such as those from building damages, the change or movement of scatters, or the variations of the sum of scatter returns can lead to temporal decorrelation for a specified location.Figure 3 demonstrates the concept that the coherence of an image pair decreased during the Meinong earthquake.Under this condition, the coherence of the coseismic image pair would be lower than that of the preseismic image pair.The coherence difference between two image pairs is useful information to detect disaster areas after the earthquake event.
Ishitsuka et al. [32] analyzed the coherence difference for mapping the soil liquefaction areas induced by the 2011 Tohoku earthquake.They assumed that the components γ N and γ S were identical and the soil liquefaction was represented by a decrease in the γ T component between preseismic and coseismic image pairs.The coherence difference (γ di f f ) influenced by the temporal decorrelation before and after the earthquake can have the formula: where γ p is the preseismic coherence and γ c is the coseismic coherence.(γ is the difference of temporal correlation between preseismic and coseismic image pairs.Following the concept in Equation ( 7), we then focused on processing the Sentinel-1 images (Copernicus Sentinel data 2016) based on two different orbits (ascending and descending) before and after the Meinong earthquake.The information of image pairs is listed in Table 1.Note that the time format of acquisition for the SAR image is YYYY/MM/DD.In this study, the perpendicular baseline (Bperp) and time difference were controlled under 150 m and 24 days, respectively.Figure 4A-C shows the preseismic coherence, coseismic coherence and the coherence difference for ascending image pairs, while Figure 4D-F shows the preseismic coherence, coseismic coherence and their coherence difference for descending image pairs.In Figure 4A,B,D,E, the warm color represents a high correlation.The correlation had values ranging from 0-1.The coherence values in western hills (i.e., along the Tsengwen River and Yueshui River) were relatively low as compared with those in the urban areas.In Figure 4C,F, the values of coherence difference were between −1 and 1.The warm color indicates a positive value, indicating that the coherence decreased from the preseismic image pair to the coseismic image pair.Such high coherence differences could provide the evidence for high variations of the land surface before and after the earthquake.Based on the processed images from either ascending or descending orbits, Figure 4C,F consistently show many hot spots with high coherence differences.These high coherence difference areas, including Xinshi, Annan, Xinhua, Yongkang and Guanmiao districts, had coherence difference values higher than 0.7.However, the associated noise was relatively high and might have influenced our analyses.Many locations exhibited uncertain information based on the temporal decorrelation processing.To extract the coherence difference solely caused by an earthquake, statistical manipulations applied to the images are essential to remove the background effects.Such background effects typically involve the atmospheric effect and ordinal surface changes [32].
The 30 Sentinel-1 images between January 2015 and January 2016 (15 ascending orbit images and 15 descending orbit images) were collected, and these images were processed for 28 coherence images (see Tables 2 and 3).The time format of acquisition for the SAR image is YYYY/MM/DD.The preseismic coherence image (the ascending orbit was 2016/01/09-2016/02/02 and the descending orbit was 2016/01/11-2016/02/04) was considered as the reference and was used to subtract all previous coherence images.The mean and standard deviation values for the previous coherence difference images were used as the threshold for selecting every pixel influenced by the earthquake.Following the concept proposed in Ishitsuka et al. [32], we used three standard deviation values to be the first threshold.In the normal distribution, the value of three standard deviations represents 99.73 percent coverage of the data.Based on this concept, the first threshold (γ thre di f f ) can be written as: where γ mean di f f (ρ, a) is the mean value at pixel (ρ, a) and γ std di f f (ρ, a) is the standard deviation at pixel (ρ, a).When the coherence difference value between preseismic and coseismic image pairs was higher than the first threshold at the same pixel, we considered this pixel to be the candidate of a disaster area.Ishitsuka et al. [32] used the ALOS PLASAR images from Japan Aerospace Exploration Agency to detect soil liquefaction.The features of the satellite such as the resolution, the radar wavelength and the repeat cycle, are quite different as compared with the Sentinel-1 satellite.ALOS provides images with relatively high spatial resolution (10 m), and the L-band with 23 cm wavelength provides better penetration for vegetation.However, the repeat cycle of ALOS is 46 days.The relatively large sampling rate of ALOS images then motivated the use of Sentinel-1 images for a particular sampling rate in our study.To improve the accuracy of selected pixels to detect the disaster areas, more thresholds might be required to filter out the candidate pixels that are directly influenced by an earthquake.
Zebker and Villasenor [23] assessed the temporal decorrelation of different surface features, including desert, forest and lava, for a time difference of 18 days.Their results showed that the values of desert correlation were constantly reaching 1, and they found that the lava had slightly less correlation than that of the desert.The correlation of forest appeared to be a linear fashion, and most correlation values dropped down to 0.5 in the study areas.Their study was the first to indicate the natural coherence decay of the unstable scatterer.In their study, the coherence difference for the forest was almost 0.5, with a time difference of 20 days.Lee et al. [45] reported that the natural temporal decay of phase coherence was not a linear behavior.However, we focus on assessing the variation of coherence decay influenced by an earthquake event in this study.The extreme coherence decay value in a constant time period could be an efficient threshold to selected the pixels of coherence difference.The coherence decay progress with time might not be significant for the analysis.
In addition, it is clear that a perfect relationship can have a correlation of 1. Correlations of 0.75 and 0.5 are considered to be strong and moderate relationships, respectively.Additionally, a weak relationship should have a value less than 0.25.When the coherence drops to a value less than 0.5, it means that the correlation is lower than a moderate relationship, and the coherence difference is higher than natural coherence decay.Therefore, in this study, we used the value of 0.5 to be the second threshold to filter out the effect of surface features.The process of the two selected thresholds was named as the two-step threshold method in this study and can be written as: where γ p−c di f f (ρ, a) is the coherence difference, which has passed the two thresholds, at pixel (ρ, a) between preseismic and coseismic images pairs.We determined the pixel selected by Equation (9) to be the coseismic coherence difference (CCD).Figure 5 summarizes the workflow of the two-step threshold method for the study.Note that the thresholds determined here might be adjusted based on site-specific problems.

Background Effects
To reduce the temporal and spatial effects, this study used 28 background image pairs (14 ascending orbit image pairs and 14 descending orbit image pairs) for the analyses.These image pairs generally had a spatial baseline of less than 200 m and had a temporal baseline of less than 24 days.Moreover, the day gap between the same pair number with different orbits was kept to 12 days.It is helpful to view the coherence difference with different look angles in a similar time period.More detailed information relevant to the image pairs is shown in Tables 2 and 3.
Figure 6 shows the results of background coherence differences for ascending orbits.The top two rows of the figure are the time series coherence difference.The lower-left quarter of Figure 6 is the mean of coherence differences, while the lower-right quarter is the standard deviation of the coherence differences.The time series of coherence differences were obtained from the preseismic coherence, which is based on processing 2016/01/09-2016/02/02 image pair and then subtracting all background coherence results.The mean and standard deviation of coherence differences were calculated by the results of 14 coherence difference maps.Figure 7 indicates the results of descending orbits based on the same processes.The reference image pairs were 2016/01/11-2016/02/04.In Figures 6 and 7, the warm color of coherence difference represents the preseismic coherence value that was higher than the background coherence.The bright color for the subfigure at the lower-right quarter shows the high variations (i.e., the standard deviation) of coherence variations.
The time series of coherence difference for ascending orbits presented consistent results.No obvious variation pattern was obtained during this time period (the top two rows in Figure 6).The positive values of coherence difference covered nearly 80% of the study area.Such observations suggest that the coherence values of the reference image pair might be generally higher than those of all background image pairs.We found that the areas with negative coherence difference approximately coexisted with positive coherence differences.However, the time series of coherence difference in descending orbit varied obviously (the top two rows in Figure 7), especially in the Anping, West Central, North, East and South districts.In these districts, the coherence differences showed negative values in image pairs Pair 02, Pair 05 and Pairs 09-12.We knew that the lagoons, canals and ponds are located near the coastal area in Anping and other southern districts in Tainan City.The regional effects might be induced by the disturbance of water vapor involved in the paths of radar waves between the reference image pairs.
In ascending and descending orbits, the mean values of coherence difference kept consistently between −0.4 and 0.6 (lower-left quarters of Figures 6 and 7).The standard deviation values were constrained under 0.2 (lower-right quarters of Figures 6 and 7), and even some variations and the burst issues might have occurred in the descending orbits (i.e., Pair 08 and Pair 13 in Figure 7).In summary, the spatial distribution of coherence values before the Meinong earthquake was relatively stable and could be used to remove the background effects in our study.

CCD Map
Following the concept described in Equation ( 9), the values of the mean and standard deviation at each pixel provided the threshold value by using three standard deviations.We then used the natural coherence decay value of 0.5 to filter out the unstable scatterers to obtain the CCD in our study.Figure 8 shows the results of the CCD that integrated the image pairs from ascending and descending orbits.Based on the concept in Figure 3, the CCD can represent the decreased coherence influenced by the earthquake.Although it is possible that the earthquake might have caused the increase of coherence in the coseismic image pair, there was no obvious negative value of coherence difference in Figure 4C.Moreover, the time series of coherence differences in Figures 6 and 7 (the top two rows) indicated a similar preseismic coherence in the previous year.The preseismic coherence was relatively more stable than the coseismic coherence.Based on the results in the study, there was no event-induced decrease of coherence in the preseismic image pair.The negative coherence difference between preseismic and coseismic image pairs might not be the issue in the study.

Comparison with Field Surveys
After the earthquake, the Central Geological Survey (CGS) of the Ministry of Economic Affairs (MOEA), Taiwan, conducted field surveys to identify the disasters in Tainan City.The disasters were classified into three different types, including destroyed buildings, ruptured structures and sand boils.The destroyed buildings were defined based on the conditions that buildings were tilted or collapsed and were difficult to reconstruct.The ruptured structures had land surface ruptures and cracks or the structure of the buildings were damaged, but the building could be reconstructed.The sand boils were defined based on the observations at sites that had obvious sand ejection with water and land surface subsidence [7].
The CCD results were identical to the results of a field survey by CGS (in Figure 8) and accurately indicated the main disasters in five districts, including Annan, Xinshin, Xinhua, Yongkang and Guanmiao districts.More specifically, the land subsidence events caused by soil liquefaction in Annan district were precisely identified by the CCD map in the study.There was no disaster near the Houchiali fault.However, two different disaster types were distributed along the Hsinhua fault.On the eastern side of the Hsinhua fault, the disaster type was mainly structural damage, while on the western side, it was mainly sand boils.The behavior of this particular distribution was the same as that shown in the 1946 Magnitude 6.3 Hsinhua earthquake [46][47][48].There was a clear boundary for two disaster types located between the plains and foothills.This phenomenon implies that the hilly areas in Tainan City might have relatively low groundwater levels in the unconfined aquifer.[7]).The orange area is the union results from the CCD map obtained from the ascending and descending image pairs.Small gray stars represent ruptured structure disasters.Large gray stars and dark gray circles indicate destroyed buildings and sand boil disasters, respectively.The pictures show the field survey records in five different areas (black rectangles on the map).
Comparing the CCD results and the CGS field surveys gives general insight into the linkage between occurrences of sand boils and variations of the coherence.However, detecting the coherence difference in the area only with the sand ejection phenomenon is dominated by the situations of surrounding environments.In special situations such as the clusters of sand boils in rice farms in Xinhua districts, the CCD map failed to identify the coherence differences accurately.Because the earthquake occurred during the early time period for paddy fields to plant rice crops, the water covered the paddy fields and provided a perfect reflector for estimating the coherence (see Area 2 of Figure 8).Even if the sand boils changed the land form of the paddy fields, the coherence differences were difficult to obtain in many locations in Xinhua district.However, many farms were abandoned in Yongkang district, and the bare lands provided a good environment for imaging processes.In these areas, the sand boils led to the change of land form and to additional water on the land.Therefore, considerable drops of coherence values in these locations were obtained after the earthquake.

Comparison with the Land Use Investigation Map
Santoro et al. [49] indicated that the coherence quality is highly dependent on the types of the land cover, especially the vegetation.To check the accuracy of CCD results, the CCD results was overlaid on the website of National Geographic Information System (http://tmap.geospatial.org.tw) and compared with the land use investigation maps (LUIM) from the National Land Surveying and Mapping Center (NLSC), Ministry of the Interior [50].
The LUIM was based on the 25 cm-resolution aerial images from the Aerial Survey Office, Forestry Bureau, Council of Agriculture, Executive Yuan, and classified into nine land use categories.We sorted the categories of LUIM into three simple classes of land covers, non-vegetation cover, vegetation cover and water cover.The non-vegetation cover mainly includes the traffic, building, public and mining land use categories; the vegetation cover includes the agricultural, forest and rest land use categories; the water cover is the hydraulic land use category, which includes rivers, canals, lakes, lagoons and ponds.
The CCD results in Annan, Xinhua, Xinshi and Yongkang districts (in Areas 1-4 of Figure 8) were validated by the three simple classes of land covers.The statistics were compiled based on the comparison between the CCD result and LUIM and are shown in Table 4.In Annan, Xinshi and Yongkang districts, more than 85% of the CCD was located on the non-vegetation cover areas of ground truth.In Xinhua, however, the CCD was partially located on the vegetation cover areas (nearly 25%), where the radar signal might have been influenced by the surrounding vegetation environments.Altogether, the CCD results are reliable for mapping the disasters in the urban area during this earthquake event.

Comparison of CCD and Liquefaction Potential and PGA Maps
In this section, we compare the CCD results with the liquefaction potential map (LPM) and the PGA map and discuss their relationship.Figure 9 shows the comparison of CCD and the LPM.The LPM was created based on data from 3D stratum, groundwater levels and the seismic hazards with a 475-year return period [51].In the study area, the development of the 3D stratum relied on standardized 133 geological and 2160 engineering boreholes.The depth of the 3D stratum was 30 m.An interpolation algorithm was used to construct the horizontal extent of physical characteristics in the 3D stratum.The mean and standard deviation of local groundwater levels were obtained based on data from the groundwater monitoring network developed by the Water Resource Agency (WRA) of the MOEA.The seismic hazards return period was based on the assumption that the dominated earthquake event was Mw 6.7 and the PGA level was between 274.4 and 424.34 gal [9,51].The liquefaction potential index was classified into three categories (i.e., high, middle and low) by the method developed by Iwasaki et al. [52].The PGA data recorded the shaking amplitudes during the earthquake event and can be related to the induced disasters.The blue dashed lines in Figure 9 are the PGA contours induced by the Meinong earthquake.This contour map was calculated by 643 stations involved in the P-alert network.The highest PGA value in this area was greater than 400 gal and was located in Xinhua district [4].The PGA values in Figure 9 decrease from the west to the south, and the values reached 200 gal in Anping, North, East, Yongkang and Xinshi districts.The results clearly show that most liquefaction phenomena investigated by CGS [7] are actually located in areas between 200 gal and 400 gal.However, the CGS report also showed that some liquefaction disasters happened in the areas with the PGA lower than 200 gal.Note that the the moment magnitude scale of the Meinong earthquake was lower than that used in the LPM, although the PGA level for the Meinong earthquake was similar to the PGA level used in the LPM.
In this study, most CCD results were mainly located in areas with middle and high liquefaction potential and in PGA levels between 200 and 400 gal, except for the clusters in Annan district.Field surveys indicated that the residential liquefaction in Annan district occurred within the old pond area, which had poor quality of underground backfill [9].Additionally, the engineering boreholes used in 3D stratum data were typically available along the highways and high-speed rail line [9,51].
The insufficiency and inequality of the data distribution might lead to a high uncertainty of the interpolated 3D stratum.These factors could lead to the underestimation of PGA level and LPM in Annan district.
In Figure 9, Guanmiao district is clearly outside of the high LPM released by Taiwan CGS.However, the Meinong earthquake event severely destroyed buildings and led to many ruptured structures.The field surveys also identified many liquefaction events after the earthquake [7].Lin [53] showed that Guanmiao district is located at the terrace deposits of the Holocene, and the main sedimental materials are composed of sand, mud and gravel.Under such terrace deposits, the formation is the Chiting Formation, which is mainly thick-bedded sandstone intercalated with thin beds of mudstone.Kuanmiao syncline is an open fold and crosses the northern part of Guanmiao district with the Chiting Formation exposed in the axial part and in both limbs.The above descriptions of the geologic structure indicate good environmental conditions for a groundwater reservoir.Moreover, our field observations showed that the groundwater level was generally one meter below the land surface in this district.Therefore, the high groundwater levels in the unconfined aquifer meet the liquefaction conditions in Guanmiao district.The strong ground shaking would rearrange the near-surface materials such as clay or fine sand even if the area were located in the foothills.

Overview of Surface Deformation
Huang et al. [22] ignored the surface displacement in the north-south direction and calculated the line-of-sight (LOS) displacements of radar images in the ascending and descending orbits for inversion of 2D (east-west direction and vertical direction) surface displacements during the Meinong earthquake.Their results agreed with the coseismic movement records obtained from the GPS.In this study, we employed the same method and estimated 2D surface displacements in CCD area.Here, the uniform grid size was 250 m × 250 m in the entire area.Figure 10 shows the conversion results and is similar to the results of Huang et al. [22].The warm color in east-west surface displacement expresses the movement toward the eastern direction (see Figure 10A), while the warm color in vertical surface displacement indicates the uplift movement (see Figure 10B).[9], and Area (II) was modified from the Central Geological Survey (CGS) field surveys [7].
The specific observations from the 2D inversion of surface displacement are summarized as follows: (1) In Annan district, the eastward displacement was close to 40 mm and the subsidence displacement was near −30 mm.(2) In Xinshin district, the westward displacement reached near −40 mm and the subsidence displacement was close to −20 mm.(3) In Xinhua district, the eastward movement was from 40-60 mm and the uplift displacement was around 10 mm after the earthquake event.(4) In Yongkang district, the eastward displacement reached 20 mm and the vertical movement varied from −5-5 mm.(5) Guanmiao district had the greatest variation of east-west movement, which varied from −100-30 mm.Additionally, the obtained high variation of vertical movement in Guanmiao district was from −20-50 mm.
In the area between Xinhua and Xinshi districts, the eastward displacement reached 50-100 mm.However, no field survey supports the values of deformation obtained in the study.To ensure the accuracy of regional-scale results, here we focus on the comparison of the 2D conversion results with filed surveys and discuss the uncertainties.Additionally, we quantify the tectonic and regional displacements during the earthquake event.The following sections provide detailed discussion for these issues.

Land Subsidence in Annan District
Tsai et al. [9] pointed out the subsidence zones in Annan district.Their study disagreed with the general results from CGS field surveys [7].However, the most serious liquefaction area was the same.The former marked more western areas to be the serious land subsidence zones, while the latter drew more land subsidence zones on the eastern side.In this study, the CCD map more clearly showed results that include these two regions (upper-left corner of Figure 10B).Most importantly, the subsidence center was the same as that obtained from field surveys.
To further address the resolution issue on a regional-scale problem, here we focus on the discussion of land subsidence in Annan district and compare the CCD map with the field measurements from the study of Tsai et al. [9].We assumed that the entire subsidence area (40,000 m 2 ) existed in one grid cell.Ground Failures (GF) 1 and 2 were defined as follows: the GF1 area was 20,000 m 2 , and the average subsidence was 5 cm; while the GF2 area was 19,500 m 2 , and the average subsidence was 17.5 cm.GF3 included five buildings with serious subsidence on Lane 161, Huian Street.The total amount of subsidence for the buildings was: 60 cm at No. 6, 90 cm at No. 8, 40 cm at No. 12 and No. 14 and 50 cm at No. 24.Assuming that the building size was 10 m × 10 m and the average subsidence in the GF3 area was 56 cm, the above information leads to an average vertical displacement of −8.28 cm in one grid cell (250 m × 250 m).In fact, the real conditions showed that the subsidence areas were located in four grid cells and the locations were nonuniform.If we consider the subsidence areas that had a uniform separated evenness, the average vertical displacement in each grid cell would be estimated as −2.07 cm, which is identical to the inversion displacement in the vertical direction.
In summary, the inversion of 2D surface displacement obtained results similar to the field measurements, but slight uncertainty still exists.The sources of the uncertainty might be the ignorance of the north-south directional movement, the low spatial resolution after the inversion and the obvious ionosphere effects [54].However, current resolution with approximately 250 m in length can monitor regional-scale disasters with reliable accuracy.

Surface Deformation in Guanmiao District
The high variation of the CCD map for Guanmiao district was mainly shown in the urban areas (see lower-right corners in Figure 10A,B).The overlapping maps of CCD and surface displacement show that the surface displacement behaved differently in the northern and southern portions of the urban area.We observed movements of 30-100 mm in the west direction and about 30-50-mm uplift in the northern portion.However, the southern portion presented movements of 10-30 mm in the east direction and exhibited 10-20-mm subsidence.The field surveys from CGS [7] reported that the northern portion in Guanmiao district was compressed westward and uplifted from a few centimeters to tens of centimeters.However, the southern portion yielded a few centimeters of subsidence based on observation of tensional fractures and liquefaction.Again, the results from CCD and the inversion of surface displacement matched well with the observations from CGS field surveys.
By integrating field measurements and InSAR estimates, Le Béon et al. [8] neglected the urban area of Guanmiao district and focused on the structurally meaningful cracks on countryside roads to explain the geologic structure activities.They argued that the extensional surface deformations were from 30-50 mm in the east-west direction in western Guanmiao district.The study concluded that the compressional displacement of 80-mm LOS in eastern Guanmiao district could be too gentle to result in the field measurements.Moreover, the Kuanmiao syncline penetrates into the northern part of the urban area in Guanmiao district [53].The axis activities of the Kuanmiao syncline might also lead to different movement mechanisms.
ElGharbawi and Tamura [35] reported that a tilted building could be a signature of a local deformation.Assuming that all buildings in the northern urban area of Guanmiao district tilted only 0.5 • toward the west direction, this would make a house 10 m in height (approximately three floors) move 87 mm in a westerly direction.In fact, the systematical building tilt is difficult to observe in most field surveys.Therefore, the field measurements of displacements were systematically lower than those from InSAR techniques.In the northern urban area in Guanmiao district, the observed displacements could be the integrated information from both tectonic movement and tilted building deformation.

Conclusions
This study proposed an efficient workflow of a two-step threshold method to remove the background effects and obtain CCD pixels in Tainan City during the Meinong earthquake.By validating with collected field survey data and sorted LUIM, we can develop a CCD map to identify multiple types of earthquake-induced disasters for regional-scale problems.Based on the results of this study, we can conclude as follows: (1) The CCD could accurately map the disaster locations in numerous districts, including Xinshi, Annan, Xinhua, Yongkang and Guanmiao.More specifically, the land subsidence induced by the soil liquefaction was clearly distinguished by the CCD results.Obvious land surface changes such as the change from sand and rocks to mud with water can result in significant drops of coherence and were also identified by the developed CCD map.(2) Except for Annan district, most high CCD values were located in middle and high potential liquefaction areas and were distributed in areas of PGA between 200 and 400 gal.Special conditions of the disasters in Annan district may result from the previous land use types and human activities.The disasters were mainly located in the so-called old pond area, where the lowlands had poor quality of backfill materials.Moreover, the low resolution of engineering boreholes in the areas might also have led to uncertainties in quantifying the index of liquefaction potential.
(3) The study found that Guanmiao district encountered severe liquefaction in many locations.
We examined the sedimentary environment and groundwater levels in Guanmiao district and found that the sand and mud deposits associated with shallow groundwater levels could be the key factors for soil liquefaction, even if the district is in the foothill area and is not included in the LPM released by Taiwan CGS.(4) The conversion of 2D surface displacement can be used to estimate the surface displacements in east-west and vertical directions.Summarizing the estimated displacements in Tainan City, Guanmiao district had high variations of east-west and vertical displacements.In other districts, such as Annan, Xinshin, Xinhua and Yongkang, the displacements could be east-west or the opposite way, depending on the location.(5) Our results showed that most differences between the estimated vertical displacements and field surveys were on the order of a few centimeters in Annan district.The inaccuracy may be influenced by neglecting surface movements in the north-south direction, the coarse spatial resolution or ionosphere effects.(6) The surface displacements in Guanmiao district behaved differently in the urban areas.The main factor could be inferred from the axis activities of the Kuanmiao syncline on the northern side.Additionally, tectonic movements and tilted buildings might also contribute to the calculation of displacements.

Figure 1 .
Figure 1.Study area.Tainan City is located in southwest Taiwan.The Roman numerals (from I-XV) are the 15 districts of Tainan City in the study area.The focal mechanism shows the epicenter location.The blue lines indicate rivers.The red lines indicate the locations of two active faults: the Hsinhua and Houchiali faults.The dashed rectangle is the study area.

Figure 2 .
Figure 2. Coseismic interferometric results in Tainan City.(A) The results of ascending image pairs.The master image 2016/02/02, and the slave image was 2016/02/14.(B) The results of descending image pairs.The master image was 2016/02/04, and the slave image was 2016/02/28.There are some small and clear fringes (diameters are from 1 km-3 km) that exist in study area.The time format of acquisition for the SAR image is provided as YYYY/MM/DD.

Figure 3 .
Figure 3. Conceptual diagram that shows the decrease of coherence before and after an earthquake.SAR: synthetic aperture radar.

Figure 4 .
Figure 4.The coherence difference obtained from the preseismic and coseismic image pairs.The subfigures show the preseismic coherence, coseismic coherence and their coherence difference obtained by using ascending image pairs (A-C) and descending image pairs (D-F).In (A,B,D,E), the value of the coherence color bar is from 0-1, and the warm color represents a high correlation value.However, the color bar for coherence difference in (C,F) is from −1-1, and the warm color indicates the decrease of a coherence value from the preseismic image pair to the coseismic image pair.

Figure 5 .
Figure 5.The workflow of the two-step threshold method.

Figure 6 .
Figure 6.The background coherence difference results for ascending orbits.The top two rows are the time series of coherence differences for the ascending background image pairs.The subfigure in the lower-left quarter displays the mean values of coherence differences for ascending background image pairs.The subfigure in the lower-right quarter shows the values of standard deviations for ascending background image pairs.

Figure 7 .
Figure 7.The background coherence difference results for descending orbits.The top two rows are the time series of coherence differences for the descending background image pairs.The subfigure in the lower-left quarter displays the mean values of coherence differences for descending background image pairs.The subfigure in the lower-right quarter shows the values of standard deviations for descending background image pairs.

Figure 8 .
Figure 8.Comparison of the CCD map and the field survey map (modified from the Central Geological Survey (CGS)[7]).The orange area is the union results from the CCD map obtained from the ascending and descending image pairs.Small gray stars represent ruptured structure disasters.Large gray stars and dark gray circles indicate destroyed buildings and sand boil disasters, respectively.The pictures show the field survey records in five different areas (black rectangles on the map).

Figure 9 .
Figure 9.The comparison of CCD with liquefaction potential and peak ground acceleration (PGA) maps.The dark gray areas show the CCD results.The liquefaction potential map was classified into three categories.The high, middle and low potential areas are marked as red, yellow and green, respectively.The blue dashed lines are the PGA contour map (modified from Wu et al. [4]).

Figure 10 .
Figure 10.Inversion of 2D surface displacements for the study area: (A) The distribution of east-west directional surface displacement, where the warm color represents eastward movement and the cold color represents westward movement.(B) The distribution of vertical surface displacement in the study area, where the warm color shows the upward movement and the cold color shows the downward movement.Red dashed lines in the enlarged figures enclose the local subsidence area.Area (I) was modified from the results in Tsai et al.[9], and Area (II) was modified from the Central Geological Survey (CGS) field surveys[7].

Table 1 .
The information of preseismic and coseismic image pairs for the study.

Table 2 .
The information of background image pairs for ascending orbits.

Table 3 .
The information of background image pairs for descending orbits.

Table 4 .
The identifiable accuracy of CCD results.