Evolution Assessment of Mining Subsidence Characteristics Using SBAS and PS Interferometry in Sanshandao Gold Mine, China

Ground subsidence is a common geological phenomenon occurring in mining areas. As an important Chinese gold mine, Sanshandao Gold Mine has a mining history of 25 years, with remarkable ground subsidence deformation. Mining development, life security, property security and ecological protection all require comprehension of the ground subsidence characteristics and evolution in the mining area. In this study, the mining subsidence phenomenon of the Sanshandao Gold Mine was investigated and analyzed based on Persistent Scatterer Interferometry (PSI) and small baseline subset (SBAS). The SAR (synthetic aperture radar) images covering the study area were acquired by the Sentinel-1A satellite between 2018 and 2021; 54 images (between 22 February 2018 and 25 May 2021) were processed using the PSI technique and 24 images (between 11 April 2018 and 12 July 2021) were processed using the SBAS technique. In addition, GACOS (generic atmospheric correction online service) data were adopted to eliminate the atmospheric error in both kinds of data processing. The interferometric synthetic aperture radar (InSAR) results showed a basically consistent subsidence area and a similar subsidence pattern. Both InSAR results indicated that the maximum LOS (line of sight) subsidence velocity is about 49 mm/year. The main subsidence zone is situated in the main mining area, extending in the northwest and southeast directions. According to the subsidence displacement of several representative sites in the mining area, we found that the PSI result has a higher subsidence displacement value compared to the SBAS result. Mining activities were accompanied by ground subsidence in the mining area: the ground subsidence phenomenon is exacerbated by the increasing mining quantity. Temporally, the mining subsidence lags behind the increase in mining quantity by about three months. In summary, the mining area has varying degrees of ground subsidence, monitored by two reliable time-series InSAR techniques. Further study of the subsidence mechanism is necessary to forecast ground subsidence and instruct mining activities.


Introduction
In recent years, ground subsidence has caused growing ecological problems and social stability concerns worldwide. As a common geological phenomenon, it is driven by the soil type, hydrogeological conditions, bedrock configuration, tectonic faults and human activities [1][2][3][4][5]. Moreover, the subsidence mechanism is complicated. Ground subsidence induced by the groundwater table lowering is related to the external loading that varies with time. Different types of soil (sand, silt, clay, gravel, etc.) have distinct compressibility and hydraulic conductivity [6][7][8][9]; soil thickness also affects settlement. Under the action of external loading, soils with different properties will have different mechanical responses. Furthermore, plenty of factors that trigger ground subsidence exist [10][11][12][13][14]; underground mining activities are a vital contributor [15][16][17][18][19][20][21][22]. The deformation and movement of rock mass induced by mining activities is a complicated mechanical problem. During mining, the formation of the goaf causes stress redistribution, deformation and destruction of rock mass until a new balance is reached. Due to the self-weight of the rock mass and the existence of the free surface, the rock mass movement eventually leads to ground subsidence. When rock mass has been repeatedly disturbed, the ground subsidence shape varies from its original and normal profile [23]. In addition, the steep-lying mining of minerals also has an important influence on the ground subsidence shape [24]. Most mining areas are adjacent to densely populated towns, where ground subsidence causes flooding, the degradation of farmland and soil erosion [22], which, in turn, poses a threat to the population and infrastructure as well as raising conflicts between residents and mining operators. In the past few years, the Interferometric Synthetic Aperture Radar (InSAR) technique has been broadly applied to various ground deformation detection purposes, from monitoring subsidence [25][26][27], landslides [28][29][30], tectonics [31], the coseismic deformation fields of earthquakes [32], to volcanology [33][34][35]. In addition, InSAR has the ability to detect the ground deformation that results from the underground storage of CO 2 [36] and cryospheric processes [37]. Compared with traditional monitoring methods, such as optical remote sensing, GPS, levelling, and total station, InSAR technology has high temporal and spatial coverage [32,38]. Moreover, it has millimeter-level line-of-sight (LOS) deformation precision and the ability to uninterruptedly acquire earth surface information regardless of the weather conditions [39][40][41][42][43]. InSAR techniques for detecting and monitoring ground subsidence have been demonstrated to be very reliable on a global level [44,45]. However, traditional InSAR technology is significantly affected by spatiotemporal coherence and atmospheric delay effects [46]. During the past few decades, advanced InSAR technology was developed on the basis of traditional InSAR technology: the time-series InSAR technology, such as the Small Baseline Subset (SBAS) technique and Persistent Scatterer Interferometry (PSI) technique, was proposed and widely used [47][48][49][50][51][52]. Moreover, time-series InSAR technology presents a period of recording the changeable radar signal return phase, which represents the LOS distance between the ground surface and the radar platform [53].
Based on advanced technology development, a large number of ground subsidence studies [54][55][56][57][58][59] have focused on coal mining. Nevertheless, studies on metal mining using InSAR technology are scarce [60][61][62][63], especially in cases of a complicated geological structure; the mining of vein metal deposits with faults as the roof in Sanshandao is one of only a few examples worldwide. Due to its great industrial importance, human geographic background and tourism-related economic value to Sanshandao Town, since it includes Sanshandao Gold Mine and gold coast resorts, the ground subsidence within the mining area needs to be investigated urgently. Moreover, due to the outbreak of novel coronavirus pneumonia, field subsidence monitoring is unachievable.
In this study, we investigated the mining-induced subsidence with time-series In-SAR analysis using Sentinel-1A images, examined the accuracy of the analysis by crossvalidating the results of the PSI and SBAS techniques and analyzed the patterns of the subsidence deformation within the Sanshandao Gold Mine over the period between 2018 and 2021. Hereafter, further study on the Sanshandao ground subsidence pattern will be summarized in our next study, which is meaningful for forecasting potential hazards and guiding mining activities.

Region of Interest
The area of interest, Sanshandao Gold Mine, is located in Sanshandao Town, Laizhou City, bordering the Bohai Sea in the northwest (Figure 1). It has a temperate monsoon continental climate, according to the monitoring data of the Laizhou Meteorological Bureau. In addition, its average temperature in the past 40 years is 12.5 • C and the average annual precipitation reaches about 600 mm. The study area is between longitudes 119.94 • and 119.96 • , and latitudes 37.39 • and 37.41 • . Overall, the terrain of the entire mining area is flat over a distance of 5 km (between 1 and 6 m in elevation) [64]. There are around 20,000 full-time residents in the study area [65]. Moreover, Sanshandao Town (a coastal tourist town) has a large floating population, with provincial highways 218 and 304 passing through the territory.
annual precipitation reaches about 600 mm. The study area is between longitudes 119.94°and 119.96°, and latitudes 37.39°and 37.41°. Overall, the terrain of the entire mining area is flat over a distance of 5 km (between 1 and 6 m in elevation) [64]. There are around 20,000 full-time residents in the study area [65]. Moreover, Sanshandao Town (a coastal tourist town) has a large floating population, with provincial highways 218 and 304 passing through the territory.

Geological Settings
The geotectonic location of the Sanshandao Gold Mine is in the Jiaobei uplift area of the Jiaobei terrane on the southern margin of the North China Platform, bound by the Yishu fault zone, the most active fault zone in eastern China, to the west [66], and by the Jiaolai depression to the south. The Sanshandao gold deposit is controlled by the northnortheast (NNE) to the northeast (NE)-trending Sanshandao-Cangshang fault that belongs to the Yishu fault zone [67] (Figure 2). There are three controlling faults, F1, F2, and F3 ( Figure 2), as well as some typical secondary joints in the mining area: Fault F1 (dipping southeastward) is the ore-bearing fault (reverse fault) with an extension depth of about 1000 m, which is located in the northeastern section of the Sanshandao-Cangshang fault zone [68]. The gold ore body is altered rock type gold deposits controlled by the F1 fault. The attitude of F1 is NE 35-40°/SE 38° ("NE 35-40°" refers to a fault strike in the 35-40° range, "SE" refers to a fault dip, "38°" refers to the dip angle of a fault, the following expressions have the same meanings) with slight variations in the strike and dip directions [69]. Fault F2 is a small-scale normal fault with an attitude of about NE 15°/NW 60-80°. The strike-slip fault F3 has the attitude of NW 300°/SW 90°, with an extension depth of about 600 m, extending westward to Laizhou Bay. Quaternary sediment (about 25 m) covers the mining area. In terms of hydrogeological features, quaternary sediment is divided into the following four units from top to bottom: (a) the first aquifer (medium to coarse sand) ranges in thickness from 4 to 7 m; (b) the first aquiclude (sandy clay) ranges in thickness from 7 to 8 m; (c) the second aquifer (coarse sand and gravel) has an average thickness of 4 m; and (d) the second aquiclude (red clay), which has an average thickness of 4 m, lies above the weathered crust ( Figure 3). The bedrock consists of sericitic and silicic granite, pyrite sericite, gneiss, cataclasite, etc. The in-situ stress field of the Sanshandao Gold Mine is dominated by horizontal stress: The maximum principal stress lies in a nearly horizontal direction (NWW), and its value is larger than the self-weight stress of the rock mass. The bedrock is subjected to strong weathering and tectonic effects, and fractures are

Geological Settings
The geotectonic location of the Sanshandao Gold Mine is in the Jiaobei uplift area of the Jiaobei terrane on the southern margin of the North China Platform, bound by the Yishu fault zone, the most active fault zone in eastern China, to the west [66], and by the Jiaolai depression to the south. The Sanshandao gold deposit is controlled by the north-northeast (NNE) to the northeast (NE)-trending Sanshandao-Cangshang fault that belongs to the Yishu fault zone [67] (Figure 2). There are three controlling faults, F1, F2, and F3 (Figure 2), as well as some typical secondary joints in the mining area: Fault F1 (dipping southeastward) is the ore-bearing fault (reverse fault) with an extension depth of about 1000 m, which is located in the northeastern section of the Sanshandao-Cangshang fault zone [68]. The gold ore body is altered rock type gold deposits controlled by the F1 fault. The attitude of F1 is NE 35-40 • /SE 38 • ("NE 35-40 • " refers to a fault strike in the 35-40 • range, "SE" refers to a fault dip, "38 • " refers to the dip angle of a fault, the following expressions have the same meanings) with slight variations in the strike and dip directions [69]. Fault F2 is a small-scale normal fault with an attitude of about NE 15 • /NW 60-80 • . The strike-slip fault F3 has the attitude of NW 300 • /SW 90 • , with an extension depth of about 600 m, extending westward to Laizhou Bay. Quaternary sediment (about 25 m) covers the mining area. In terms of hydrogeological features, quaternary sediment is divided into the following four units from top to bottom: (a) the first aquifer (medium to coarse sand) ranges in thickness from 4 to 7 m; (b) the first aquiclude (sandy clay) ranges in thickness from 7 to 8 m; (c) the second aquifer (coarse sand and gravel) has an average thickness of 4 m; and (d) the second aquiclude (red clay), which has an average thickness of 4 m, lies above the weathered crust ( Figure 3). The bedrock consists of sericitic and silicic granite, pyrite sericite, gneiss, cataclasite, etc. The in-situ stress field of the Sanshandao Gold Mine is dominated by horizontal stress: The maximum principal stress lies in a nearly horizontal direction (NWW), and its value is larger than the selfweight stress of the rock mass. The bedrock is subjected to strong weathering and tectonic effects, and fractures are relatively developed, forming a densely developed fracture zone of fault F1. According to the data sharing infrastructure of the National Earthquake Data Center (http://data.earthquake.cn, accessed on 8 December 2021), there have been four earthquakes with a depth of less than 7 km and a magnitude below 5 in Laizhou city thus far ( Ms 3.0; Location: 120.02 • E, 37

Mining Settings
The gold mine began its basic construction in August 1984. At the end of 1990, the production capacity of mining reached its designed scale, with an expected production capacity of 1500 t/d. The exploitation activities mainly developed along the strike direction of the F1, from southwest to northeast. The mechanical upward horizontal slice stoppingfilling method with pointed pillars is adopted in this mining area. A stope with good ventilation and drainage is positioned along the strike, which has a length of 100 m, and the width is the thickness of the ore body. At present, the exploitation length of the −375 m middle section is the longest, about 2.54 km along the strike direction of the F1. The exploitation length of the deep middle section is gradually increasing, and the deepest exploiting middle section has reached −825 m. At present, the exploitation of the shallow

Mining Settings
The gold mine began its basic construction in August 1984. At the end of 1990, the production capacity of mining reached its designed scale, with an expected production capacity of 1500 t/d. The exploitation activities mainly developed along the strike direction of the F1, from southwest to northeast. The mechanical upward horizontal slice stoppingfilling method with pointed pillars is adopted in this mining area. A stope with good ventilation and drainage is positioned along the strike, which has a length of 100 m, and the width is the thickness of the ore body. At present, the exploitation length of the −375 m middle section is the longest, about 2.54 km along the strike direction of the F1. The exploitation length of the deep middle section is gradually increasing, and the deepest exploiting middle section has reached −825 m. At present, the exploitation of the shallow

Mining Settings
The gold mine began its basic construction in August 1984. At the end of 1990, the production capacity of mining reached its designed scale, with an expected production capacity of 1500 t/d. The exploitation activities mainly developed along the strike direction of the F1, from southwest to northeast. The mechanical upward horizontal slice stoppingfilling method with pointed pillars is adopted in this mining area. A stope with good ventilation and drainage is positioned along the strike, which has a length of 100 m, and the width is the thickness of the ore body. At present, the exploitation length of the −375 m middle section is the longest, about 2.54 km along the strike direction of the F1. The exploitation length of the deep middle section is gradually increasing, and the deepest exploiting middle section has reached −825 m. At present, the exploitation of the shallow gold ore body has basically ended, and the depth of the mining activities is close to −1000 m. During the mining period, the gold ore body was exploited at different depths.

SAR Images
The SAR images used in this study were acquired by Sentinel-1A, which is the Earth observation satellite in the European Space Agency (ESA)'s "Copernicus" plan, following ERS and ENVISAT SAR missions [49]. The advantage of the Sentinel-1A satellite is that it provides C-band (5.6 cm wavelength) SAR data that are considerably enhanced in aspects of coverage, repeat cycle, etc. In addition, the Sentinel-1A satellite with a total swath width of 250 km (pixel size: 5 m × 20 m) has the potential to map the global landmasses once every 12 days according to the acquisition plan of various sites [70,71]. Fifty-seven radar images of the study area were obtained using the C-band sensor on board the Sentinel-1A satellite. The images used for InSAR analysis in this study were collected between 22 February 2018 and 12 July 2021. In addition, the images were obtained in the Interferometric Wide swath (IW) scanning mode (angle of view θ = 38.9 • ; swath width = 250 km) with VV polarization and along the ascending orbit (track 171). Specific information on the SAR images is shown in Table 1. Twenty-four of these fifty-seven images were used for SBAS analysis, which were collected between 11 April 2018 and 12 July 2021. Fifty-four of these fifty-seven images were used for time series InSAR analysis using the PSI technique. The Precise Orbit Determination data (obtained from the European Space Agency (ESA) https://scihub.copernicus.eu/gnss/#/home, accessed on 8 December 2021) and the digital elevation model from the Shuttle Radar Topography Mission [72] with a resolution of 30 m were used to remove the orbit error and topography phase. The corresponding unit is displayed in parentheses.

Data Processing of Time-Series InSAR Technology
The SBAS technique was developed by the CNR-IREA group in Italy [48]. It was improved to deal with the full resolution interferograms [73] and applies multi-looking [74]. Its critical characteristics are the utilization of abundant SAR images distributed in small baseline subsets and the singular value decomposition (SVD) method [48,75]. In terms of analyzing distributed scatterers, the SBAS technique performed better than the PSI technique in capturing strongly nonlinear displacement rates [76]. The SBAS method enables us to generate multi-master interferograms instead of single-master interferograms. In this study, the SBAS processing was implemented through isce [77] and mintpy (https://github.com/insarlab/MintPy, accessed on 8 December 2021) software. The stack included 24 SAR images from the Sentinel-1A satellite (polarization: VV, orbit configuration: ascending), dating from 11 April 2018 to 12 July 2021.
In recent years, many studies have focused on the following two aspects of the PSI technique: detecting permanent scatterers [78][79][80] and modeling [48,49,79,81]. The first concept of the PSI technique is based on the identification of Persistent Scatterers (PS) whose phase has stable amplitude values and consistently low noise levels over a long time [48,78]. The key steps are as follows: establish a model based on the PSs, solve the deformation phase, and separate the atmospheric delay phase in order to obtain the deformation information.
In this study, PSI processing was performed with the support of SNAP (Sentinel Application Platform) software (ESA, 2017) [82] and the StaMPS (Stanford Method for Persistent Scatterers) algorithm. Forming interferograms from single-look complex (SLC) images is the preprocessing target of SNAP [5,34,[83][84][85][86]. It does this by applying precise orbit correction to co-register the stack. choosing the master image. generating interferograms and removing the topographic phase, and exporting the processed stack into Gamma format to enable PSI processing with the StaMPS algorithm (version 3.3b1). Moreover, the Goldstein adaptative filter (window size: 32; value of α: 0.8) [87] was applied during the phase unwrapping process to reduce the noise phase. The flowchart of the whole PSI processing chain is shown in Figure 4.
formation information.
In this study, PSI processing was performed with the support of SNAP (Sentinel Application Platform) software (ESA, 2017) [82] and the StaMPS (Stanford Method for Persistent Scatterers) algorithm. Forming interferograms from single-look complex (SLC) images is the preprocessing target of SNAP [5,34,[83][84][85][86]. It does this by ① applying precise orbit correction to co-register the stack. ② choosing the master image. ③ generating interferograms and removing the topographic phase, and ④ exporting the processed stack into Gamma format to enable PSI processing with the StaMPS algorithm (version 3.3b1). Moreover, the Goldstein adaptative filter (window size: 32; value of α: 0.8) [87] was applied during the phase unwrapping process to reduce the noise phase. The flowchart of the whole PSI processing chain is shown in Figure 4.
As mentioned before, the PSI analysis included 54 radar images (ascending orbit) and covered the period from 22 February 2018 to 25 May 2021. The date 20,200,804 was used as a master scene for track 171 in PSI processing. The Shuttle Radar Topography Mission (SRTM) 30 m data were used to remove the topographic phase contribution. In addition, it is a characteristic of the PSI technique that the retained pixels of the subsidence velocity maps are located mainly over artificial structures. Furthermore, we applied GACOS data that are broadly used by scholars to eliminate atmospheric errors in both kinds of data processing [88][89][90]. The reason why we do not use all available images from the monitoring period is as follows. Some images of the study area are corrupted and unavailable, and our computer configuration is not powerful enough. In consideration of computational efficiency and image date regularity, we only select a part of the images for analysis. This is because the PSI technique and SBAS technique have different key mechanisms, and these mechanisms need a different number of images. For instance, the SBAS technique is an analysis method for the multi-master-image InSAR time series that uses interferometric pairs of small baselines to obtain the surface deformation, whereas, the PSI technique is an analysis method for the mono-master-image. In addition, the perpendicular baselines and temporal baselines of the images are important for forming interferometric image pairs. Considering the As mentioned before, the PSI analysis included 54 radar images (ascending orbit) and covered the period from 22 February 2018 to 25 May 2021. The date 20,200,804 was used as a master scene for track 171 in PSI processing. The Shuttle Radar Topography Mission (SRTM) 30 m data were used to remove the topographic phase contribution. In addition, it is a characteristic of the PSI technique that the retained pixels of the subsidence velocity maps are located mainly over artificial structures. Furthermore, we applied GACOS data that are broadly used by scholars to eliminate atmospheric errors in both kinds of data processing [88][89][90].
The reason why we do not use all available images from the monitoring period is as follows. Some images of the study area are corrupted and unavailable, and our computer configuration is not powerful enough. In consideration of computational efficiency and image date regularity, we only select a part of the images for analysis. This is because the PSI technique and SBAS technique have different key mechanisms, and these mechanisms need a different number of images. For instance, the SBAS technique is an analysis method for the multi-master-image InSAR time series that uses interferometric pairs of small baselines to obtain the surface deformation, whereas, the PSI technique is an analysis method for the mono-master-image. In addition, the perpendicular baselines and temporal baselines of the images are important for forming interferometric image pairs. Considering the problems of image date regularity and the temporal baseline, we used a dataset with a slight difference.

Spatial Characteristics of LOS Velocity Fields
The PSI processing identified 7218 PS within the selected area, which means an average density of 150 PS/km 2 . The mean LOS velocity maps of various archives are calculated by estimating the linear regression of displacement over time. As shown in Figure 5, the maximum deformation rate was about −49 mm/year, furthermore, the standard deviation of the mean LOS velocity result is shown in Figure 6, which is an indication of the result's quality. In addition, the LOS velocity value produced by PSI is a relative value, this explains why positive values exist. In addition, the outcome of the SBAS analysis demonstrated that the mining area shows a maximum subsidence velocity in LOS of about 49 mm/year (Figure 7). the two subsidence funnels obtained using the PSI and SBAS techniques share the same scope, between 119.94° E and 119.96° E, and 37.39° N and 37.41° N. Figure 8 shows the important construction areas in the study area. As displayed in Figure 9, a general trend of negative LOS velocity is distributed in the important construction areas. The heavily subsided area is lower than the surrounding area. In the rainy summer, rainwater gathers due to poor drainage (Figure 10), affecting transport in Sanshandao Town. However, the specular reflection of radar waves by water surface results in missing the deformation data in the fish farm area. There is a gradual change in LOS velocity from negative values to positive values in areas away from the important constructions. It turns out that the subsidence may be bounded by the F3 fault ( Figure 2). As Figure 9 shows, the ground deformation extends along the northwest and southeast direction, similar to an oblique egg in shape (the same part is outlined by a solid, yellow, rectangular line in Figure 1, using the PSI and SBAS techniques).      According to the LOS deformation results, we discovered a similar deformation pattern through two time-series InSAR techniques: the shape of the subsidence funnel resembles an egg. At the same time, they all possess a stratified structure. The external part has a lower LOS velocity value, and the inner part has a higher LOS velocity value. Moreover, the two subsidence funnels obtained using the PSI and SBAS techniques share the same scope, between 119.94 • E and 119.96 • E, and 37.39 • N and 37.41 • N. Figure 8 shows the important construction areas in the study area. As displayed in Figure 9, a general trend of negative LOS velocity is distributed in the important construction areas. The heavily subsided area is lower than the surrounding area. In the rainy summer, rainwater gathers due to poor drainage (Figure 10), affecting transport in Sanshandao Town. However, the specular reflection of radar waves by water surface results in missing the deformation data in the fish farm area. There is a gradual change in LOS velocity from negative values to positive values in areas away from the important constructions. It turns out that the subsidence may be bounded by the F3 fault ( Figure 2). As Figure 9 shows, the ground deformation extends along the northwest and southeast direction, similar to an oblique egg in shape (the same part is outlined by a solid, yellow, rectangular line in Figure 1, using the PSI and SBAS techniques).

Comparison between PSI and SBAS Results
As mentioned above, we attained similar deformation results using the PSI and SBAS techniques applied to the same study area. In order to cross-validate the specific deformation results and advance the detailed study, we selected six sites represented by the letters R1, R2, S1, S2, SC1, and SC2 (The Google Earth images of six critical sites are displayed in Figure 11). Due to the fact that PSI and SBAS data are both relative results, one stable zero-deformation point (37.372° N, 119.934° E) was selected to be the reference in

Comparison between PSI and SBAS Results
As mentioned above, we attained similar deformation results using the PSI and SBAS techniques applied to the same study area. In order to cross-validate the specific deformation results and advance the detailed study, we selected six sites represented by the letters R1, R2, S1, S2, SC1, and SC2 (The Google Earth images of six critical sites are displayed in Figure 11). Due to the fact that PSI and SBAS data are both relative results, one stable zero-deformation point (37.372° N, 119.934° E) was selected to be the reference in

Comparison between PSI and SBAS Results
As mentioned above, we attained similar deformation results using the PSI and SBAS techniques applied to the same study area. In order to cross-validate the specific deformation results and advance the detailed study, we selected six sites represented by the letters R1, R2, S1, S2, SC1, and SC2 (The Google Earth images of six critical sites are displayed in Figure 11). Due to the fact that PSI and SBAS data are both relative results, one stable zero-deformation point (37.372 • N, 119.934 • E) was selected to be the reference in PSI and SBAS analysis processing, and it is far away from the mining subsidence funnel with low LOS velocity and standard deviation value (as shown by the black dot in Figure 7). Consistent time series behavior was discovered after contrasting several representative sites in the constructions area. As displayed in Figure 12, the LOS displacements in the Residential Area zone are relatively consistent, as are those shown in Figures 13 and 14. Figure 13 shows the LOS displacements results in the sea view community zone;, and Figure 14 shows the LOS displacements of various sites in the area behind the school. The results indicated the high consistency of subsidence monitoring by time-series InSAR techniques. In the area of import construction (Figures 8 and 11), the sea view community area (represented by S1, S2), the residential area (represented by R1, R2) and the area behind the school (represented by SC1, SC2) have different levels of severity. Most of the buildings in these areas have wall cracking. In general, damages to buildings due to slope and curvature are relatively slight, and the main cause of building damage is the tensile and compressible deformation, with the former being more harmful to buildings. We believe that both PSI and SBAS results are plausible; the degree of building damage is obviously different in the different construction areas, as shown in Figures 15-17. The buildings in the sea view community which were made of reinforced concrete were less damaged. PSI and SBAS analysis processing, and it is far away from the mining subsidence funnel with low LOS velocity and standard deviation value (as shown by the black dot in Figure  7). Consistent time series behavior was discovered after contrasting several representative sites in the constructions area. As displayed in Figure 12, the LOS displacements in the Residential Area zone are relatively consistent, as are those shown in Figures 13 and 14. Figure 13 shows the LOS displacements results in the sea view community zone;, and Figure 14 shows the LOS displacements of various sites in the area behind the school. The results indicated the high consistency of subsidence monitoring by time-series InSAR techniques. In the area of import construction (Figures 8 and 11), the sea view community area (represented by S1, S2), the residential area (represented by R1, R2) and the area behind the school (represented by SC1, SC2) have different levels of severity. Most of the buildings in these areas have wall cracking. In general, damages to buildings due to slope and curvature are relatively slight, and the main cause of building damage is the tensile and compressible deformation, with the former being more harmful to buildings. We believe that both PSI and SBAS results are plausible; the degree of building damage is obviously different in the different construction areas, as shown in Figures 15, 16 and 17. The buildings in the sea view community which were made of reinforced concrete were less damaged.

Spatiotemporal Correlation between InSAR Results and Mining Activities
As mentioned before, we performed PSI processing with the 54 radar images obtained from the Sentinel-1A satellite. The cumulative LOS displacements in several representative sites were provided in the previous section. Moreover, we attained more relatively continuous subsidence displacement values compared with SBAS processing. Therefore, we chose the PSI results as the following discussion object. Firstly, we broke the cumulative LOS displacement results into two parts for analysis, using the date 20,200,801 as the cut-off point. Secondly, we performed linear fitting on the two parts of the data. Table 2 shows that the first part has a low subsidence velocity and the second part has a higher subsidence velocity. Multiple potential mechanisms that trigger the subsidence phenomenon could exist, as demonstrated by the distinct temporal behavior of the LOS displacement value. For instance, this phenomenon can be explained by changes in mining quantity.

Spatiotemporal Correlation between InSAR Results and Mining Activities
As mentioned before, we performed PSI processing with the 54 radar images obtained from the Sentinel-1A satellite. The cumulative LOS displacements in several representative sites were provided in the previous section. Moreover, we attained more relatively continuous subsidence displacement values compared with SBAS processing. Therefore, we chose the PSI results as the following discussion object. Firstly, we broke the cumulative LOS displacement results into two parts for analysis, using the date 20,200,801 as the cut-off point. Secondly, we performed linear fitting on the two parts of the data. Table 2 shows that the first part has a low subsidence velocity and the second part has a higher subsidence velocity. Multiple potential mechanisms that trigger the subsidence phenomenon could exist, as demonstrated by the distinct temporal behavior of the LOS displacement value. For instance, this phenomenon can be explained by changes in mining quantity.

Spatiotemporal Correlation between InSAR Results and Mining Activities
As mentioned before, we performed PSI processing with the 54 radar images obtained from the Sentinel-1A satellite. The cumulative LOS displacements in several representative sites were provided in the previous section. Moreover, we attained more relatively continuous subsidence displacement values compared with SBAS processing. Therefore, we chose the PSI results as the following discussion object. Firstly, we broke the cumulative LOS displacement results into two parts for analysis, using the date 20,200,801 as the cut-off point. Secondly, we performed linear fitting on the two parts of the data. Table 2 shows that the first part has a low subsidence velocity and the second part has a higher subsidence velocity. Multiple potential mechanisms that trigger the subsidence phenomenon could exist, as demonstrated by the distinct temporal behavior of the LOS displacement value. For instance, this phenomenon can be explained by changes in mining quantity. Dates are in the format yyyy mm dd; The unit of velocity is mm/year; R 2 is the goodness of fit.

Discussion
The water supplies for production and daily life in Sanshandao Town come from the urban centralized water supply measure. There is no need to exploit the groundwater of the study area for water supply, which means that the effective stress of the underlying soil does not increase significantly as the water table falls. Therefore, ground subsidence induced by groundwater over-exploitation, as a result of an external loading varying with time, is almost impossible. Moreover, the structure of the overlying soil layer in the study area is relatively uniform. In other words, the deformation of the overlying soil is related to the deformation of the underlying rock mass. For these reasons, the ground subsidence changes the original stress state of construction, causing construction deformation damage.
However, both of the LOS displacement results in the same site are not totally coincidental, and the common character of the displacement results indicated that the PSI displacement value is bigger than the SBAS displacement value. We selected several sites (as shown in Figure 11: R1-R4, S1-S4, SC1-SC4) for error analysis. As shown in Figure 18, the scatter plot with regular dispersion shows a low standard deviation. We used IBM SPSS (version 22) for error analysis: the standard deviation is about 6 mm, and the deviations do not have normal error distribution (Kolmogorov-Smirnov consistency test [91,92]), which implies that the results are reliable. The absolute deviation result is displayed in the lower right corner of Figure 18. We calculated the root mean square error (RMSE) between the results of the PS and SBAS. The RMSE value is about 8 mm. In general, according to the results of the previous chapters and the error analysis, the results obtained using the PSI and SBAS techniques are in good agreement from both the qualitative and quantitative perspectives. The potential reasons for this phenomenon are described below: both time-series InSAR techniques have different principles and algorithms, as well as work platforms. For instance, the SBAS technique is an analysis method for the multi-master-image InSAR time series that uses interferometric pairs of small baselines to obtain the surface deformation; the PSI technique is an analysis method for the mono-master-image. SBAS regards better nonlinear movements, and the SBAS technique reduces the decorrelation by exploiting distributed targets; PSI enables a precise characterization of the linear deformation affecting urban areas, and PSI can detect more PS in urban areas. Therefore, the results are not in perfect accord.
As displayed in Figure 9, the subsidence deformation may be bounded by the F3 fault, as the ground deformation extends along the northwest and southeast direction, similar in shape to an oblique egg. We have selected several sites in the dashed box (black) in Figure 11 from the northwest to the southeast for LOS velocity analysis. It turns out that the subsidence rate is different on both sides of the fault. The south area has a bigger LOS velocity value than the north area ( Figure 19). From a qualitative and quantitative perspective, our results suggested that the fault F3 may restrict and affect the subsidence development. The ground subsidence shape is affected by the discontinuous interface (in this case, the fault) which is the product of internal and external dynamic geological processes. In the case of mining, the heterogeneous rock mass with geological discontinuities affects stress redistribution. Therefore, the ground subsidence can be influenced [5,93,94]. In addition, there have been no earthquakes in Laizhou City since 2018 according to the earthquake data. The subsidence does not appear to have been caused by strong fault movement. Unfortunately, we cannot discriminate the mechanism between ground subsidence and fault attributes due to a lack of enough underground mining information and field tests. Further studies on this aspect will be summarized in our next study. As displayed in Figure 9, the subsidence deformation may be bounded by the F3 fault, as the ground deformation extends along the northwest and southeast direction, similar in shape to an oblique egg. We have selected several sites in the dashed box (black) in Figure 11 from the northwest to the southeast for LOS velocity analysis. It turns out that the subsidence rate is different on both sides of the fault. The south area has a bigger LOS velocity value than the north area ( Figure 19). From a qualitative and quantitative perspective, our results suggested that the fault F3 may restrict and affect the subsidence development. The ground subsidence shape is affected by the discontinuous interface (in this case, the fault) which is the product of internal and external dynamic geological processes. In the case of mining, the heterogeneous rock mass with geological discontinuities affects stress redistribution. Therefore, the ground subsidence can be influenced [5,93,94]. In addition, there have been no earthquakes in Laizhou City since 2018 according to the earthquake data. The subsidence does not appear to have been caused by strong fault movement. Unfortunately, we cannot discriminate the mechanism between ground subsidence and fault attributes due to a lack of enough underground mining information and field tests. Further studies on this aspect will be summarized in our next study. The underground rock mass is heterogeneous due to the existence of geological discontinuities. Before mining, the rock mass under the effect of the in-situ stress field is in a balanced state. After the local ore body is exploited, the goaf is formed inside the rock mass, which leads to changes in the stress state of the surrounding rock mass, causing stress redistribution. This change in the stress states induces a deformation that could reach the failure threshold. Throughout this process, the rock mass deforms, becomes broken, and generates numerous fractures. Finally, the underlying bedrock deforms and the overlying quaternary soil settles down. Many scholars have studied the coal mining subsidence [43,[54][55][56][57][58][59], but studies on steeply dipping metal deposits are scarce. Coalfields are more regular than metal deposits and thus form a more regular subsidence area. The underground rock mass is heterogeneous due to the existence of geological discontinuities. Before mining, the rock mass under the effect of the in-situ stress field is in a balanced state. After the local ore body is exploited, the goaf is formed inside the rock mass, which leads to changes in the stress state of the surrounding rock mass, causing stress redistribution. This change in the stress states induces a deformation that could reach the failure threshold. Throughout this process, the rock mass deforms, becomes broken, and generates numerous fractures. Finally, the underlying bedrock deforms and the overlying quaternary soil settles down. Many scholars have studied the coal mining subsidence [43,[54][55][56][57][58][59], but studies on steeply dipping metal deposits are scarce. Coalfields are more regular than metal deposits and thus form a more regular subsidence area.
According to the schematic diagram of ground subsidence (Figure 20), the settlement shape is affected by the shape of the ore body [24]: when the ore body is horizontal and the goaf is located in the center of the ore body, the subsidence boundary is symmetrical and the location of the subsidence center corresponds to the position of the goaf; when the ore body is steep-lying and the goaf is located in the center of the ore body, the sedimentary boundary (asymmetrical) and subsidence center moves in the opposite direction to the ore body inclination. According to the detailed deformation figures (Figure 9), we discovered that the displacement value is asymmetric. We circled the settlement center with a smaller dotted ellipse; and we circled one LOS velocity contour with a bigger dotted ellipse. The subsidence center is roughly biased in the uphill direction (northwest). In other words, the displacement gradient (between two ellipses) in the northwest of the subsidence center is larger. These results imply the possibility that the subsidence of the study area is related to the geometric configuration of the ore body. Further studies on this aspect will be summarized in our next study. To analyze the relationship between the subsidence velocity and mining activities, we collected the mining quantity information (between January 2018 and December 2020) of the mining area above −600 m and selected R3, R2, S1 and S2 for analysis. Since June 2020, the mining quantity (above −600 m) has doubled compared with previous months, and has continued to grow exponentially. Consequently, the cumulative LOS subsidence displacement also changed with the mining quantity ( Figure 21). Remarkably, the results indicate that the occurrence of subsidence lags behind mining and excavation activities. This result may be due to the fact that the gold ore body was exploited at different depths during the mining period. During mining, the formation of the goaf causes stress redistribution, deformation and destruction of the rock mass until a new balance is reached. It takes time for the underground rock mass and the overlying soil to deform. As shown in Figure 21, a sudden increase in mining quantity does not significantly change the amount of subsidence, and this is because the two different measuring units are not comparable; the rock has a greater density (2.78 t/m 3 ) and the total goaf volume has a small change. According to the administrative staff of the Sanshandao Gold Mine, the total mining quantity has remained almost stable in recent years. However, the mining quantity of the mining area (above −600 m) has increased dramatically since June 2020. To analyze the relationship between the subsidence velocity and mining activities, we collected the mining quantity information (between January 2018 and December 2020) of the mining area above −600 m and selected R3, R2, S1 and S2 for analysis. Since June 2020, the mining quantity (above −600 m) has doubled compared with previous months, and has continued to grow exponentially. Consequently, the cumulative LOS subsidence displacement also changed with the mining quantity ( Figure 21). Remarkably, the results indicate that the occurrence of subsidence lags behind mining and excavation activities. This result may be due to the fact that the gold ore body was exploited at different depths during the mining period. During mining, the formation of the goaf causes stress redistribution, deformation and destruction of the rock mass until a new balance is reached. It takes time for the underground rock mass and the overlying soil to deform. As shown in Figure 21, a sudden increase in mining quantity does not significantly change the amount of subsidence, and this is because the two different measuring units are not comparable; the rock has a greater density (2.78 t/m 3 ) and the total goaf volume has a small change. According to the administrative staff of the Sanshandao Gold Mine, the total mining quantity has remained almost stable in recent years. However, the mining quantity of the mining area (above −600 m) has increased dramatically since June 2020.
Remote Sens. 2022, 14, x FOR PEER REVIEW 21 of 25 Figure 21. The cumulative LOS subsidence displacement and mining quantity change (the turning points of subsidence displacement and mining quantity are circled by yellow ellipses).

Conclusions
In this study, we monitored the ground subsidence over the Sanshandao Gold Mine area between 2018 and 2021 using PSI and SBAS techniques. The subsidence results derived using the time-series InSAR techniques have high consistency, indicating that the mining area has a similar deformation pattern. The main subsidence deformation area is located between 119.94° E and 119.96° E, and 37.39° N and 37.41° N, corresponding to the main mining area. In addition, the outcome of the PSI and SBAS analysis indicated that the maximum deformation rate is about −49 mm/year. Spatially, the subsidence deformation fits with the F3 fault, extending along the northwest and southeast direction, similar in shape to an oblique egg. Temporally, the subsidence velocity significantly accelerated after July 2020, providing insights into the time of changes in mining activities; the subsidence velocity changed over the mining boom period, with a three-month delay. In general, our results help to better characterize the spatiotemporal evolution of ground subsidence. The findings of this study suggest that mining quantity is a vital factor affecting ground subsidence in the Sanshandao Gold Mine area. Due to a lack of detailed mining material, this study did not address the specific mechanism of ground subsidence in the Sanshandao Gold Mine area. Therefore, the local government should pay more attention to residential areas with severe subsidence. Measures to displace the population can be taken when necessary. Due to the complexity of the geological environment, more intensive field investigations of ground subsidence and the collection of adequate mining data are necessary for in-depth researches in the future.

Conclusions
In this study, we monitored the ground subsidence over the Sanshandao Gold Mine area between 2018 and 2021 using PSI and SBAS techniques. The subsidence results derived using the time-series InSAR techniques have high consistency, indicating that the mining area has a similar deformation pattern. The main subsidence deformation area is located between 119.94 • E and 119.96 • E, and 37.39 • N and 37.41 • N, corresponding to the main mining area. In addition, the outcome of the PSI and SBAS analysis indicated that the maximum deformation rate is about −49 mm/year. Spatially, the subsidence deformation fits with the F3 fault, extending along the northwest and southeast direction, similar in shape to an oblique egg. Temporally, the subsidence velocity significantly accelerated after July 2020, providing insights into the time of changes in mining activities; the subsidence velocity changed over the mining boom period, with a three-month delay. In general, our results help to better characterize the spatiotemporal evolution of ground subsidence. The findings of this study suggest that mining quantity is a vital factor affecting ground subsidence in the Sanshandao Gold Mine area. Due to a lack of detailed mining material, this study did not address the specific mechanism of ground subsidence in the Sanshandao Gold Mine area. Therefore, the local government should pay more attention to residential areas with severe subsidence. Measures to displace the population can be taken when necessary. Due to the complexity of the geological environment, more intensive field investigations of ground subsidence and the collection of adequate mining data are necessary for in-depth researches in the future.