Long-Term Remote Monitoring of Ground Deformation Using Sentinel-1 Interferometric Synthetic Aperture Radar (InSAR): Applications and Insights into Geotechnical Engineering Practices

: Development of synthetic aperture radar (SAR) technology and the dedicated suite of processing tools have aided the evolution of remote sensing techniques for various Earth Observation (EO) applications. Interferometric SAR (InSAR) is a relatively new geodetic technique which provides high-speed and reliable geographic, geologic, and hazards information allowing the prognosis of future environmental and urban planning. In this study, we explored the applicability of two di ﬀ erential interferometry techniques, conventional and advanced di ﬀ erential InSAR (A-DInSAR), for topographic mapping and long-term geotechnical monitoring by exploiting satellite data, particularly Sentinel-1 SAR data, which is publicly shared. We speciﬁcally used the open-source tools of SeNtinel Application Platform (SNAP) and Stanford Method for Persistent Scatterers (StaMPS) for interferometric data processing to implement A-DInSAR. This study presents various applications, which include generation of a digital elevation model (DEM), mapping of seismically induced displacement and associated damages, and detection and long-term monitoring of tunneling-induced ground deformation and rainfall-induced landslide. Geometric and temporal decorrelations posed challenges and limitations in the successful implementation of Sentinel-1 SAR interferometry speciﬁcally in vegetated areas. The presented results proved the validity and reliability of the exploited SAR data and InSAR techniques for addressing geotechnical engineering related problems. Advanced di ﬀ erential interferometric synthetic aperture radar (A-DInSAR) (based on Permanent Scatterer InSAR (PS-InSAR) technique) processing workﬂow through integration of SNAP and StaMPS. The workﬂow was implemented for tunneling-induced ground deformation monitoring in Site C and rainfall-induced landslide analysis in Site D.


Introduction
High precision geodetic measurements with coarse-to-fine spatial and temporal resolutions significantly complement geological, seismic, and geotechnical investigations and characterizations of geodynamic processes related to natural and human-induced activities. The continual development of synthetic aperture radar (SAR) technology and the dedicated suite of processing tools over the past years have introduced numerous opportunities in the fields of remote sensing and civil engineering. Interferometric SAR (InSAR) and advanced differential InSAR (A-DInSAR) are relatively new geodetic techniques in these fields which have proven their capability of providing centimeter-to millimeter-level change records of terrain topography and surface deformations. SAR pixel-based derived classifications and measurements promote it from the sparse point Global Positioning System (GPS) technology, while processing both SAR amplitude and phase information make it unique and advantageous over conventional imaging techniques. For DEM generation, Seogwipo, the second-largest city located in the southern part of Jeju Island, was chosen. This volcanic island houses a volcanic mountain, called Mt. Halla, at the center of the island (Figure 1), which has a height of 1950 m above sea level with an inclination of approximately 3-5 • from the peak [40]. For the earthquake-induced ground deformation analysis, the city of Pohang was selected. This city stands on the alluvium deposit located in North Gyeongsang Province and serves as a main seaport in the southeast coast of South Korea. An earthquake with the moment magnitude M w = 5.4 struck the city on 15 November 2017, one of the strongest and most damaging earthquakes in the Korean Peninsula. It was a moderate-sized earthquake with an estimated focal depth of 5 km [35]. It has an oblique strike-slip mechanism with a reverse faulting. Over 1000 aftershocks followed the mainshock, and among them, the largest event was with the magnitude 4.6 on 11 February 2018. Table A1 shows the epicenters of the major seismic events of the Pohang earthquake [41]. Moreover, a local magnitude 3.1 earthquake was relocated during the hydraulic stimulation operations at an Enhanced Geothermal System (EGS) site in Pohang [42]. The mainshock caused widespread liquefaction, and heavily damaged the town of Heunghae and many neighboring villages which in turn resulted in over 80 fatalities and left thousands of people homeless.

Site C: South Chungcheong Province
We chose a tunnel construction site in South Chungcheong Province for the tunneling-induced ground deformation analysis. The site has notable geologic characteristics, as follows. (a) This AOI is a reclamation area, which had been subjected to dramatic land development, transforming itself into a coastal industrial zone which houses several manufacturing factories and steel mills. (b) The site is located seaside close to a shore. (c) The geologic strata are composed of a Songak gneiss rock formation and a Pyeongtaek migmatitic gneiss rock formation with overlaying Quaternary alluvial deposits. In the chosen site, excavation of a vertical shaft for the tunnel construction using a shield tunnel boring machine (TBM) began in 2017. The excavation caused excessive inflow of groundwater into the excavated vertical shaft and called pumping-out of the inflowing groundwater with intermittent suspensions of excavation works. After 13 months, the horizontal tunneling of the TBM started. However, when the TBM advanced approximately 95 m after 3 months, it was again temporarily terminated, as damages in neighboring structures and ground cracks along the roads near the vertical shaft were reported, as shown in Figure 2. Table A2 briefly outlines the construction sequence of the tunnel construction schedule. Photos were taken in June 2019. Road pavement cracks (photo (e)) and masonry fences cracks (photos (c,d,g,h)) were observed during the field investigation. Path sections near the vertical shaft were either depressed or raised (photo (f)).

Site D: Gyeongju, North Gyeongsang Province
For rainfall-induced landslide analysis, the AOI covers a cut slope along a motorway and its uphill Yangbuk Mountain area in the city of Gyeongju. On 7 October 2018, more than 300 mm of rain fell under the influence of Typhoon Kongrey. This heavy rainfall caused failure of the reinforced slope of 20 m wide and 250 m long, as shown in Figure 3. Slope stabilization works which involved pile reinforcement and groundwater level control were implemented to improve safety. On 30 June Photos were taken in June 2019. Road pavement cracks (photo (e)) and masonry fences cracks (photos (c,d,g,h)) were observed during the field investigation. Path sections near the vertical shaft were either depressed or raised (photo (f)). For rainfall-induced landslide analysis, the AOI covers a cut slope along a motorway and its uphill Yangbuk Mountain area in the city of Gyeongju. On 7 October 2018, more than 300 mm of rain fell under the influence of Typhoon Kongrey. This heavy rainfall caused failure of the reinforced slope of 20 m wide and 250 m long, as shown in Figure 3. Slope stabilization works which involved pile reinforcement and groundwater level control were implemented to improve safety. On 30 June 2019, 0.70 km of the Janghang Intersection was opened to relieve the transportation inconvenience. Ramp sections connecting the Janghang Intersection to the building complex were opened on 5 July 2019. Photos were taken in June 2019. Road pavement cracks (photo (e)) and masonry fences cracks (photos (c,d,g,h)) were observed during the field investigation. Path sections near the vertical shaft were either depressed or raised (photo (f)).

Site D: Gyeongju, North Gyeongsang Province
For rainfall-induced landslide analysis, the AOI covers a cut slope along a motorway and its uphill Yangbuk Mountain area in the city of Gyeongju. On 7 October 2018, more than 300 mm of rain fell under the influence of Typhoon Kongrey. This heavy rainfall caused failure of the reinforced slope of 20 m wide and 250 m long, as shown in Figure 3. Slope stabilization works which involved pile reinforcement and groundwater level control were implemented to improve safety. On 30 June 2019, 0.70 km of the Janghang Intersection was opened to relieve the transportation inconvenience. Ramp sections connecting the Janghang Intersection to the building complex were opened on 5 July 2019.

Datasets Used
This study used Level 1 Single-Look-Complex (SLC) SAR dataset acquired by Sentinel-1 radar system equipped with C-band. This C-band has a radar wavelength of 5.6 cm. The images were acquired using the Interferometric Wide (IW) swath mode with co-vertical polarization (VV). Sentinel-1, a constellation of two satellites, is in a near-polar, sun-synchronous orbit. Therefore, it passes over the AOI in both ascending and descending flights with a 12-day repeat cycle alone. The SAR data are freely shared in public. The Sentinel-1A was launched in 2014 and the Sentinel-1B in 2016 through the COPERNICUS Programme of the European Space Agency (ESA). ESA's Sentinel-1 missions operate in Terrain Observation by Progressive Scans SAR (TOPSAR) mode which provides a SAR image with wide area coverage of around 250 × 200 km. The IW swath mode captures three swaths comprising a series of bursts. The SAR image has the moderate ground resolution of approximately 14 × 4 m in azimuth (parallel to the flight direction) and range (perpendicular to the flight direction) directions, respectively. Orbit state vectors in the metadata of the acquired SAR products are updated using the Sentinel-1 Precise Orbit Ephemerides (POE) in order to provide refined satellite position and velocity information. Aside from SLC dataset, Ground Range Detected (GRD) and Ocean (OCN) data can also be searched and ordered from the COPERNICUS Sentinel Open Access Hub or Alaska SAR Facility. Table A3 lists the characteristics of Sentinel-1 missions.
A total of 114 Sentinel-1 SAR scenes were processed: two images for Seogwipo (Site A, descending track 134), three images for Pohang (Site B, descending track 61), 78 images for South Chungcheong province (Site C, descending track 134), and 31 images for Gyeongju (Site D, descending track 61). All images were acquired at night to ensure high signal-to-noise ratio (SNR) and high coherence on the AOIs and mitigate atmospheric delay in the interferometric phase [43][44][45]. Figure 4 shows the Sentinel-1 images connection graph for the A-DInSAR analysis performed in Site C and Site D.
(GRD) and Ocean (OCN) data can also be searched and ordered from the COPERNICUS Sentinel Open Access Hub or Alaska SAR Facility. Table A3 lists the characteristics of Sentinel-1 missions.
A total of 114 Sentinel-1 SAR scenes were processed: two images for Seogwipo (Site A, descending track 134), three images for Pohang (Site B, descending track 61), 78 images for South Chungcheong province (Site C, descending track 134), and 31 images for Gyeongju (Site D, descending track 61). All images were acquired at night to ensure high signal-to-noise ratio (SNR) and high coherence on the AOIs and mitigate atmospheric delay in the interferometric phase [43][44][45]. Figure 4 shows the Sentinel-1 images connection graph for the A-DInSAR analysis performed in Site C and Site D. rainfall-induced landslide in Site D. The time interval between successive SAR acquisition is equal to or less than 24 days. All perpendicular baselines for each interferometric pair were less than 150 m. Each point in the graphs represents a SAR image and each line represents an interferometric pair. Horizontal and vertical axes define the temporal (BTemp) and perpendicular (BPerp) baselines, respectively. Table 1 lists the summary of applications of interferometry in various geotechnical engineering related problems investigated in South Korea.   Table 1 lists the summary of applications of interferometry in various geotechnical engineering related problems investigated in South Korea. 2.3. Sentinel-1 SAR Interferometry

Basic Concept
SAR is an active imaging system where the sensor emits microwave radiation to illuminate the Earth's surface and receives the reflected signal, called backscatterer. It works any time of the day and in any cloud-cover condition towards the surface. The backscattered signals are expressed as an image of which each pixel contains both amplitude (A) and phase (ϕ) information, recorded as complex numbers (c) [46]: where i is the imaginary number. By pairing two SAR images obtained at distant periods with disparate look angles over the same location, one image is set as master (m) and the other as slave (s), an interferogram (I) is generated. An interferogram is computed by multiplying the amplitude of both images and calculating the difference in phase (Equation (2)). The InSAR technique exploits the phase difference between ϕ m and ϕ s , named interferometric phase (ϕ ip ), expressed in Equation (3).

of 20
where ϕ elev is the topographic phase, ϕ defo is the phase contribution due to surface deformation, ϕ flat is the flat-earth phase, ϕ atm is the atmospheric phase screening, and ϕ noise is the phase noise related to some processing-induced factors. Through interferometry, the desired phase contributor, which is typically the elevation or the surface deformation, can be isolated from the other phase components. However, the final quantitative measurements of elevation or surface deformation from interferometry will greatly depend of the quality of the generated interferogram. While the SAR technology creates satellite images that cover almost the entire planet, not every resolution cell contains valuable information. The main condition for the InSAR application is the degree of coherence. When ground targets or objects change between acquisition periods, the coherent sum of backscattered signals varies resulting from coherence degradation due to decorrelation. To quantify the impact of decorrelation for our InSAR application, the coherence (γ) is calculated: where N is the coherence window size. The calculated coherence indicates how identical each pixel is between master and slave scenes from a scale of 0 (completely out of phase) to 1 (completely in phase). Any ground deformation associated with natural hazard generates different radar signal paths for the same ground target between two observation periods. Hence, radar signals received by the sensor from damaged structures such as buildings and road pavements can lead to coherence loss. In principle, under ideal condition co-event image pair coherence (γ co ) would be lower than that of the pre-event image pair coherence (γ pre ). The coherence difference (γ diff ) between those two pairs can be used to delineate damaged areas after the disastrous event. In this study, the coherence difference is calculated as follows [17,18]:

Standard Repeat-Pass Interferometry (Single Pair)
For the DEM generation in Site A (Seogwipo) and earthquake-induced ground deformation in Site B (Pohang), a standard repeat-pass interferometry was conducted by using SNAP (Sentinel-1 Toolbox) developed and freely distributed by ESA. Figure 5 presents the processing chain for the single pair standard repeat-pass InSAR.

Standard Repeat-Pass Interferometry (Single Pair)
For the DEM generation in Site A (Seogwipo) and earthquake-induced ground deformation in Site B (Pohang), a standard repeat-pass interferometry was conducted by using SNAP (Sentinel-1 Toolbox) developed and freely distributed by ESA. Figure 5 presents the processing chain for the single pair standard repeat-pass InSAR. First, coregistration of two SAR images employing the Enhanced Spectral Diversity (ESD) method was performed after accurate orbit data implementation to create the interferogram of AOI and after subtraction of flat-earth phase. During the interferogram formation, coherence was also estimated. The output is then debursted to smoothly stitch all burst data into a single image. For the DEM generation, topographic phase was not removed as it is the contributor of interest to be estimated. Moreover, multilooking was not performed as it reduces the spatial resolution of the First, coregistration of two SAR images employing the Enhanced Spectral Diversity (ESD) method was performed after accurate orbit data implementation to create the interferogram of AOI and after subtraction of flat-earth phase. During the interferogram formation, coherence was also estimated. The output is then debursted to smoothly stitch all burst data into a single image. For the DEM generation, topographic phase was not removed as it is the contributor of interest to be estimated. Moreover, multilooking was not performed as it reduces the spatial resolution of the resulting DEM product. For ground deformation estimation, the topographic phase was removed using precise orbit data and Shuttle Radar Topography Mission (SRTM) DEM 1 arcsec data. Phase filtering was then performed to increase the SNR. Phase unwrapping was performed using the SNAPHU program. Thereafter, the unwrapped phase was converted into either elevation or displacement. The converted product which was still in the radar geometry was then geocoded into the World Geodetic System 1984 (WGS84) and exported for post analysis and visualization in Google Earth and QGIS.

A-DInSAR Analysis with SNAP and StaMPS Integration
The A-DInSAR analysis based on the PS-InSAR technique was performed using the acquired Sentinel-1 free dataset. The PS-InSAR technique has emerged as a matured geodetic monitoring tool characterized by millimeter accuracy. The processing is split into two independent workflows [47][48][49] and presented in Figure 6: (1) single master DInSAR processing through SNAP, and (2)  In SNAP processing, the master image was selected using the InSAR Stack Overview tool. The selected master image was approximately in the barycenter of the dataset considering both the temporal and spatial baselines to ensure high coherence in all interferometric pairs generated over the whole stack. Image splitting and orbit file correction were applied, and slave images were coregistered to the geometry of the master image. Interferograms were generated and topographic phase was removed using SRTM DEM 1 arcsec data. Then, debursting was applied to both the original coregistered images and interferograms. Finally, the processed data was exported in a format that is compatible with StaMPS readers. The snap2stamps scripts were used to automate the process except for the master image selection [50]. It is worth noting that interferograms need to be checked before proceeding to the next part to avoid processing errors during the PS-InSAR analysis.
The second part is StaMPS processing. The first step is the preparation of the SNAP exports into StaMPS using the mt_prep_snap script included in the distribution. After running the script, Matlab was launched, and StaMPS processing was subsequently carried out from Step 1 (Load data) through Step 8 (Atmospheric filtering). Additionally, to mitigate the tropospheric phase screens, the aps_linear approach was used as part of the integrated Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [43]. Finally, the displacement time series and mean annual velocity were exported for post analysis and visualization in StaMPS-Visualizer [51] and QGIS.

DEM Generation over Site A-Seogwipo City
An interferometric pair with the shortest temporal baseline and the highest perpendicular baseline assessed the DEM based on Sentinel-1 SAR data over Seogwipo. Figure 7 shows the In SNAP processing, the master image was selected using the InSAR Stack Overview tool. The selected master image was approximately in the barycenter of the dataset considering both the temporal and spatial baselines to ensure high coherence in all interferometric pairs generated over the whole stack. Image splitting and orbit file correction were applied, and slave images were coregistered to the geometry of the master image. Interferograms were generated and topographic phase was removed using SRTM DEM 1 arcsec data. Then, debursting was applied to both the original coregistered images and interferograms. Finally, the processed data was exported in a format that is compatible with StaMPS readers. The snap2stamps scripts were used to automate the process except for the master image selection [50]. It is worth noting that interferograms need to be checked before proceeding to the next part to avoid processing errors during the PS-InSAR analysis.
The second part is StaMPS processing. The first step is the preparation of the SNAP exports into StaMPS using the mt_prep_snap script included in the distribution. After running the script, Matlab was launched, and StaMPS processing was subsequently carried out from Step 1 (Load data) through Step 8 (Atmospheric filtering). Additionally, to mitigate the tropospheric phase screens, the aps_linear approach was used as part of the integrated Toolbox for Reducing Atmospheric InSAR Noise (TRAIN) [43]. Finally, the displacement time series and mean annual velocity were exported for post analysis and visualization in StaMPS-Visualizer [51] and QGIS.

DEM Generation over Site A-Seogwipo City
An interferometric pair with the shortest temporal baseline and the highest perpendicular baseline assessed the DEM based on Sentinel-1 SAR data over Seogwipo. Figure 7 shows the interferometric phase generated and the corresponding coherence to check the quality of phase information. The interferometric fringes are clearly shown in the lowlands towards the coast, which is a good indication of topographic information, but are not visible in the vicinity of Mt. Halla. The dominant grassy and vegetated nature of the AOI in the vicinity of Mt. Halla caused considerable loss of coherence due to strong temporal decorrelation (i.e., surface change between the 12-day time interval). The mean coherence of the whole AOI is below 0.30. Figure 8 shows the InSAR-derived DEM.   Figure 9 shows the comparison to the reference model (SRTM DEM). Figure 10 shows the elevation profiles along sections A-A' and B-B'. Although the InSAR-derived DEM and the SRTM DEM appear to be almost similar to each other, the elevation profiles derived from InSAR show more fluctuation compared to the smooth SRTM elevation profiles. More importantly, the InSAR result shows that better elevation measurements can be obtained from plain areas compared to hilly areas and vegetated areas. It underestimates the elevation value towards the mountain peak; the absolute elevation difference can be as high as 250 m. For the lowlands, the absolute elevation difference between InSAR-DEM and SRTM DEM can be as low as 1 m.   Figure 9 shows the comparison to the reference model (SRTM DEM). Figure 10 shows the elevation profiles along sections A-A' and B-B'. Although the InSAR-derived DEM and the SRTM DEM appear to be almost similar to each other, the elevation profiles derived from InSAR show more fluctuation compared to the smooth SRTM elevation profiles. More importantly, the InSAR result shows that better elevation measurements can be obtained from plain areas compared to hilly areas and vegetated areas. It underestimates the elevation value towards the mountain peak; the absolute elevation difference can be as high as 250 m. For the lowlands, the absolute elevation difference between InSAR-DEM and SRTM DEM can be as low as 1 m. Figure 8. InSAR-derived DEM using Sentinel-1 SAR data processed through SNAP. Figure 9 shows the comparison to the reference model (SRTM DEM). Figure 10 shows the elevation profiles along sections A-A' and B-B'. Although the InSAR-derived DEM and the SRTM DEM appear to be almost similar to each other, the elevation profiles derived from InSAR show more fluctuation compared to the smooth SRTM elevation profiles. More importantly, the InSAR result shows that better elevation measurements can be obtained from plain areas compared to hilly areas and vegetated areas. It underestimates the elevation value towards the mountain peak; the absolute elevation difference can be as high as 250 m. For the lowlands, the absolute elevation difference between InSAR-DEM and SRTM DEM can be as low as 1 m. DEM appear to be almost similar to each other, the elevation profiles derived from InSAR show more fluctuation compared to the smooth SRTM elevation profiles. More importantly, the InSAR result shows that better elevation measurements can be obtained from plain areas compared to hilly areas and vegetated areas. It underestimates the elevation value towards the mountain peak; the absolute elevation difference can be as high as 250 m. For the lowlands, the absolute elevation difference between InSAR-DEM and SRTM DEM can be as low as 1 m.

Earthquake-Induced Ground Deformation Mapping over Site B-Pohang City
Pre-and co-seismic interferometric pairs were processed to detect the surface changes related to the 15 November 2017 earthquake. Figure 11 presents the pre-and co-seismic coherence maps. It is observed that the overall co-seismic coherence is less than the pre-seismic coherence in most urban areas located in the west, southeast, and southern regions of the epicenter. This coherence degradation indicates the damaged buildings and the inherent differential settlement caused by seismic ground motion and associated liquefaction. Meanwhile, liquefaction phenomena have been reported in the northern and western parts of the epicenter, which is consistent with the observed coherence reduction.   [41,52,53]. The scale of coherence difference is between −1 and 1. The negative coherence difference (warm color) shows areas that were modified (i.e., by natural processes) before the event but were not modified during the co-seismic interval. The positive coherence difference (cold color) represents areas that had surface changes during the co-seismic interval but little to no change before the event. In the map, positive coherence difference is found

Earthquake-Induced Ground Deformation Mapping over Site B-Pohang City
Pre-and co-seismic interferometric pairs were processed to detect the surface changes related to the 15 November 2017 earthquake. Figure 11 presents the pre-and co-seismic coherence maps. It is observed that the overall co-seismic coherence is less than the pre-seismic coherence in most urban areas located in the west, southeast, and southern regions of the epicenter. This coherence degradation indicates the damaged buildings and the inherent differential settlement caused by seismic ground motion and associated liquefaction. Meanwhile, liquefaction phenomena have been reported in the northern and western parts of the epicenter, which is consistent with the observed coherence reduction.

Earthquake-Induced Ground Deformation Mapping over Site B-Pohang City
Pre-and co-seismic interferometric pairs were processed to detect the surface changes related to the 15 November 2017 earthquake. Figure 11 presents the pre-and co-seismic coherence maps. It is observed that the overall co-seismic coherence is less than the pre-seismic coherence in most urban areas located in the west, southeast, and southern regions of the epicenter. This coherence degradation indicates the damaged buildings and the inherent differential settlement caused by seismic ground motion and associated liquefaction. Meanwhile, liquefaction phenomena have been reported in the northern and western parts of the epicenter, which is consistent with the observed coherence reduction.   [41,52,53]. The scale of coherence difference is between −1 and 1. The negative coherence difference (warm color) shows areas that were modified (i.e., by natural processes) before the event but were not modified during the co-seismic interval. The positive coherence difference (cold color) represents areas that had surface changes during the co-seismic interval but little to no change before the event. In the map, positive coherence difference is found   [41,52,53]. The scale of coherence difference is between −1 and 1. The negative coherence difference (warm color) shows areas that were modified (i.e., by natural processes) before the event but were not modified during the co-seismic interval. The positive coherence difference (cold color) represents areas that had surface changes during the co-seismic interval but little to no change before the event. In the map, positive coherence difference is found mostly in the southeastern part of the AOI, which is associated with damaged buildings in urban areas. Whereas, negative coherence difference is placed dominantly in the northern region. Most areas in this northern part from the epicenter are arable land where little to no liquefaction features were reported [54].  Figure 13 shows the co-seismic differential interferogram and the corresponding ground deformation, which were obtained by using the descending track SAR data. The interferometric fringe delineates the surface displacement in the direction of the radar's line-of-sight (LOS) in a unit of the radar half-wavelength (i.e., one fringe cycle equals 2.8 cm for Sentinel-1). The interferogram clearly reveals visible co-seismic deformation along the radar LOS. In the displacement map, uplift and subsidence areas were distinguished which caused damage to neighboring buildings. Since Sentinel-1B acquires images in a right-looking mode, positive values indicate uplift and easting, and the negative values suggest subsidence and westing. The reference area, which was assumed not to be affected by the seismic ground motion, is marked on the map and has a coherence value close to 1. The areas (or pixels) with the coherence value less than 0.30 were masked out from the map. Displacement values of more than 5 cm (uplift) from the vicinity of the epicenter were reported.  Figure 13 shows the co-seismic differential interferogram and the corresponding ground deformation, which were obtained by using the descending track SAR data. The interferometric fringe delineates the surface displacement in the direction of the radar's line-of-sight (LOS) in a unit of the radar half-wavelength (i.e., one fringe cycle equals 2.8 cm for Sentinel-1). The interferogram clearly reveals visible co-seismic deformation along the radar LOS. In the displacement map, uplift and subsidence areas were distinguished which caused damage to neighboring buildings. Since Sentinel-1B acquires images in a right-looking mode, positive values indicate uplift and easting, and the negative values suggest subsidence and westing. The reference area, which was assumed not to be affected by the seismic ground motion, is marked on the map and has a coherence value close to 1. The areas (or pixels) with the coherence value less than 0.30 were masked out from the map. Displacement values of more than 5 cm (uplift) from the vicinity of the epicenter were reported.  Figure 14 shows the LOS mean velocity for a period of 1009 days estimated by using the PS-InSAR technique from Sentinel-1B data. In the map, the color scale is between −50 and +30 mm/yr. Positive (in blue color) and negative (in red color) values denote movement towards and away from the radar LOS, respectively. Spatial distribution of the PS points was relatively dispersed, showing the complete absence of a rebound deformation trough. Due to the presence of buildings in the project site, the map shows a dense distribution of PS points except in the NW and SW corner regions, where the regions are mostly cultivation lands. A remarkable subsidence trough was observed in the industrial zone where the tunnel construction has commenced and was almost nonexistent in the SE corner region. The subsidence trough consists of yellow, orange, and red color variants depending on the magnitude of deformation.   Figure 14 shows the LOS mean velocity for a period of 1009 days estimated by using the PS-InSAR technique from Sentinel-1B data. In the map, the color scale is between −50 and +30 mm/yr. Positive (in blue color) and negative (in red color) values denote movement towards and away from the radar LOS, respectively. Spatial distribution of the PS points was relatively dispersed, showing the complete absence of a rebound deformation trough. Due to the presence of buildings in the project site, the map shows a dense distribution of PS points except in the NW and SW corner regions, where the regions are mostly cultivation lands. A remarkable subsidence trough was observed in the industrial zone where the tunnel construction has commenced and was almost nonexistent in the SE corner region. The subsidence trough consists of yellow, orange, and red color variants depending on the magnitude of deformation.  Figure 14 shows the LOS mean velocity for a period of 1009 days estimated by using the PS-InSAR technique from Sentinel-1B data. In the map, the color scale is between −50 and +30 mm/yr. Positive (in blue color) and negative (in red color) values denote movement towards and away from the radar LOS, respectively. Spatial distribution of the PS points was relatively dispersed, showing the complete absence of a rebound deformation trough. Due to the presence of buildings in the project site, the map shows a dense distribution of PS points except in the NW and SW corner regions, where the regions are mostly cultivation lands. A remarkable subsidence trough was observed in the industrial zone where the tunnel construction has commenced and was almost nonexistent in the SE corner region. The subsidence trough consists of yellow, orange, and red color variants depending on the magnitude of deformation.   Figure 15 shows the time histories of InSAR-derived displacements of the selected PS points located along the two road segments (A and B) passing through the vertical shaft. All displacements are estimated in the LOS direction. The tunnel construction phases are highlighted with grey color and marked with letters A-E (Table A2). It is observed that ground deformations near the vertical shaft accelerated during Phase C. Prior to this, no appreciable displacements were recorded along the two roads. Moreover, subsidence increased as the PS distance from the shaft decreased. The settlement started to slow down gradually at the end of Phase D after the shield tunneling was temporarily stopped. No extreme ground deformation was then observed at the end of Phase D that continued towards the end of Phase E. PS-A01, B04, and B05 show the cumulative displacements of more than 5 cm towards the end of the monitoring period.

Tunneling-Induced Ground Deformation Detection and Monitoring over Site C-Chungcheong
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 22 Figure 15 shows the time histories of InSAR-derived displacements of the selected PS points located along the two road segments (A and B) passing through the vertical shaft. All displacements are estimated in the LOS direction. The tunnel construction phases are highlighted with grey color and marked with letters A-E (Table A2). It is observed that ground deformations near the vertical shaft accelerated during Phase C. Prior to this, no appreciable displacements were recorded along the two roads. Moreover, subsidence increased as the PS distance from the shaft decreased. The settlement started to slow down gradually at the end of Phase D after the shield tunneling was temporarily stopped. No extreme ground deformation was then observed at the end of Phase D that continued towards the end of Phase E. PS-A01, B04, and B05 show the cumulative displacements of more than 5 cm towards the end of the monitoring period.  Figure 16 shows the LOS mean velocity for a period of 361 days prior to the landslide as well as the selected PS points over the target AOI. In the map, the color scale is between −18 and +16 mm/yr. The PS points were densely placed in a building complex located on the upper part of the figure and along the National Route 4 of the Janghang Intersection. These PS points showed various motions related to natural factors and seasonal changes (i.e., contraction and expansion due to temperature change).  Figure 16 shows the LOS mean velocity for a period of 361 days prior to the landslide as well as the selected PS points over the target AOI. In the map, the color scale is between −18 and +16 mm/yr. The PS points were densely placed in a building complex located on the upper part of the figure and along the National Route 4 of the Janghang Intersection. These PS points showed various motions related to natural factors and seasonal changes (i.e., contraction and expansion due to temperature change).

Rainfall-Induced Landslide Detection over Site D-Gyeongju
On the other hand, limited measurement points were present within Yangbuk Mountain due to low coherence related to significant environmental changes between the image acquisitions. The available revisit period was longer than 6 days for the target AOI. Dense vegetation over Yangbuk Mountain led to temporal decorrelation which resulted in bad interferogram quality. Figure 17 shows the temporal changes in displacement of the selected PS points located within the failed slope and its periphery (Figure 16b). There existed instability patterns throughout the analysis period. The initial plunge in the displacement history plot, which occurred between day 252 and day 288, represents a considerable surface displacement. This time interval in Gyeongju, measured with respect to the first SAR image, has the greatest number of rainfall events that last for about 13 days and typically aggregates to almost 200 mm of precipitation. This also coexists with the occurrence of the monsoon season in the region and with the continuous rain events which starts in late June until early August, displacements are triggered. The mean displacement time series plots of the two groups of measurement points are also shown ( Figure 17). On the other hand, limited measurement points were present within Yangbuk Mountain due to low coherence related to significant environmental changes between the image acquisitions. The available revisit period was longer than 6 days for the target AOI. Dense vegetation over Yangbuk Mountain led to temporal decorrelation which resulted in bad interferogram quality. Figure 17 shows the temporal changes in displacement of the selected PS points located within the failed slope and its periphery (Figure 16b). There existed instability patterns throughout the analysis period. The initial plunge in the displacement history plot, which occurred between day 252 and day 288, represents a considerable surface displacement. This time interval in Gyeongju, measured with respect to the first SAR image, has the greatest number of rainfall events that last for about 13 days and typically aggregates to almost 200 mm of precipitation. This also coexists with the

Advantages and Limitations of InSAR
In this study, our results obtained by using standard repeat-pass InSAR and A-DInSAR (based on PS-InSAR) techniques with Sentinel-1 SAR data and open-source SNAP and StaMPS software packages show its capability in generating DEM and estimating changes in deformation or damages in various geotechnical engineering practices. SAR, as an active radar imaging technique, operates day and night in all weather conditions which has greater advantage over optical remote sensing techniques. InSAR and A-DInSAR techniques have matured and become viable solutions complimenting conventional geodetic, optical, and LiDAR techniques and integrated in multidisciplinary fields. Moreover, the improved temporal and spatial resolutions, wider area coverage, and free access of Sentinel missions allow a weekly repeat observation for ground anomaly detection, and near real-time monitoring of various geohazards.
Meanwhile, the following limiting factors should be accounted for prior to practice. Geometric and temporal decorrelations are the major limiting factors in the successful implementation of Sentinel-1 SAR interferometry [55]. Geometric decorrelation occurs when the geometries of the master and slave images poorly match each other in association with the changes in a look angle and non-parallel orbits. This also happens when the baseline between the two SAR acquisitions deviates, as a result of the satellite's orbit separation in space, from zero [56]. Non-zero baseline tends to have higher geometric decorrelation. Temporal decorrelation on the other hand occurs as the time separation between SAR observation is large owing to the dynamic surface changes of the Earth. Shorter temporal baseline ensures that these changes can be captured related to the SAR phase, hence, minimizing the temporal decorrelation.
For deformation detection, mapping, and monitoring at the rest of the test sites (Sites B, C, and D), a shorter perpendicular baseline is favored to reduce the phase signal due to topography. Meanwhile, the temporal baseline depends on the magnitude and type of deformation being investigated. Thanks to the short revisit times and tight orbital navigation of Sentinel-1 constellation (Figure 4 and Table 1), these decorrelation effects were minimized. In addition to these decorrelations, atmospheric perturbations also affect the propagation of the SAR signals on standard repeat-pass InSAR that was used for the Pohang earthquake at Site B.
The proportion of Sentinel-1 s short wavelength to the size of water molecules as it comes closer to the edge of the atmospheric window makes it sensitive to atmospheric effects. The troposphere and ionosphere tend to delay and accelerate, respectively, the electromagnetic wave, which yields incorrect elevation and deformation calculations. Geometric distortions due to layover, foreshortening, and shadowing in mountainous regions could result in miscalculations especially for landslide applications (i.e., in Gyeongju) [29]. Layover occurs when the top of the mountain appears before its bottom in the radar LOS. Foreshortening exists when the slope angle is less than the incidence angle of the observation. Whereas, shadowing happens when the mountain blocks the view of observation.

Possible Improvement of InSAR-Based Detection and Monitoring
The C-band Sentinel-1 mission currently provides free, open-access, medium-resolution SAR image data with worldwide coverage, which we used in this study for various geotechnical engineering related phenomena. The datasets were processed using SNAP and StaMPS which are also open-source software packages. The potential of using shared data and tools can be further utilized for recurring geohazards caused by natural phenomena and human activities around the world. This is possible by (1) combining both Sentinel-1A/B SAR data to achieve an effective 6-day revisiting period which is shorter compared to the 12-day repeat cycle for one constellation alone (if available), (2) combining both ascending and descending data from the same orbit to obtain the true N-S and E-W ground motions, (3) acquiring short-time span and high-resolution SAR products (i.e., from COSMO-SkyMed and TerraSAR-X), and (4) using SAR data acquired with longer radar wavelength (i.e., from ALOS PALSAR 2). Furthermore, (5) no means least is the integration of InSAR techniques with other surveying methods. Use of precise leveling for monitoring of ground movements can complement the satellite-based remote sensing techniques after InSAR detects a potential threat due to ground deformations caused by natural factors and human activities. These further improvements ensure more efficient temporal decorrelation consideration and more PS detection capability even in sub-urban and/or rural areas. As one of the possible examples, these enhancements will allow to inspect and manage the tunnel construction sequence before proceeding to the next stage when untoward incidents happen. Infrastructures and people's lives can also be safeguarded against landslides if an early detection method of precursory slope movement is in place.

Conclusions
This work highlighted the limitations, performance, and the great contribution of using medium-resolution SAR images for various geotechnical engineering related phenomena. The study demonstrated the applicability of using shared Sentinel-1 SAR data provided by ESA under the COPERNICUS Programme processed with open-source software packages, SNAP and StaMPS, employing the conventional and advanced interferometric processing. Two well-defined processing workflows were presented and tested in four geotechnical engineering applications namely DEM generation, earthquake-and tunneling-induced ground deformations, and rainfall-induced landslide in the Korean Peninsula. The major limitations of SAR interferometry were also discussed, which will serve as important pointers in future projects involving continuous deformation monitoring and landslide early detection with improved satellites features.
Overall, the challenges that the authors have encountered on using open data and open-source software tools have created and promoted a conducive and scientific learning environment between remote sensing and geotechnical engineering fields. The experience gained through this work opens new opportunities to explore in addressing timely social and environmental impacts of various geohazards.  Acknowledgments: The authors would like to thank the European Space Agency (ESA) for freely supplying Sentinel-1 SAR data under the COPERNICUS Programme (https://scihub.copernicus.eu). The authors would also like to thank the United States Geological Survey (USGS) in cooperation with the National Aeronautics and Space Administration (NASA) for providing the SRTM DEM (https://earthexplorer.usgs.gov). SNAP (by ESA), StaMPS (by Andrew Hooper), and TRAIN (by David Bekaert) software packages are freely distributed under GNU-GPL license.

Conflicts of Interest:
The authors declare no conflict of interest.