Application of DInSAR and Spatial Statistics Methods in Analysis of Surface Displacements Caused by Induced Tremors

: Induced seismicity is one of the negative phenomena caused by anthropogenic activities that include mining of minerals. This phenomenon manifests itself as sudden and unpredictable shocks of rock mass, which can cause surface deformation and damage to ground infrastructure. Until the advent of satellite radar interferometry that enables analysis of historical events, the characteristics of these unexpected surface deformations were di ﬃ cult to assess. The main aim of the research was the spatial analysis of the geometry of surface displacements caused by eight induced tremors in the Rudna copper mine (SW Poland) and the dependence of deformation characteristics (vertical displacements, extent) on the induced shock energy. For this purpose, Sentinel-1 satellite imagery, the di ﬀ erential radar satellite interferometry (DInSAR) method and geographic information systems (GIS) based spatial statistics were used. Vertical displacements were mapped on the basis of 37 calculated interferograms. Spatial statistics on the pixel-to-pixel level were performed in the GIS Map Algebra environment. In the result, descriptive and spatial statistics characterizing deformations caused by individual shocks were calculated. The average values of vertical displacements ranged from − 44 to − 119 mm. Strong, statistical correlation between the extent, maximum vertical displacement, and energy values was determined. In addition, geometries of the formed deformation areas were analyzed and presented graphically. The results obtained in this research constitute development of a knowledge base on surface displacements caused by induced tremors in underground copper mining.


Introduction
Induced seismicity is a phenomenon that is directly or indirectly caused by industrial operations, such as mining of minerals, production of geothermal energy, underground storage of liquids and gases, extraction of conventional and unconventional hydrocarbons, construction of dams, and retention reservoirs. These technological processes change the stress distribution of the Earth's crust and may cause seismic shocks. The mechanism of induced seismicity begins with anthropogenic activity in the rock mass or on its surface. Engineering works violate the natural state of equilibrium in the rock mass, which releases potential energy from the rocks. Then, part of this energy changes into seismic energy, which propagates from the focus of shock as elastic waves [1]. The resulting shocks can damage infrastructure, cause deformations on the surface of the earth, and pose danger to people. In our paper, we deal with terrain deformations that are the consequence of underground mining of minerals.
Induced seismicity and natural seismicity are the same in the form of surface vibrations. It would seem that the only differences between mining tremors and earthquakes are depth and energy of the event (magnitude). However, Zembaty [2,3] presented quantitative and qualitative differences in the context of surface effects, such as differences in size and intensity of shocks, differences in the duration of kinematic excitation, differences in the spectral properties of records of mining tremors and earthquakes, differences in peak values of ground movement, and differences between horizontal and vertical components of rock mass.
In recent years, there has been a significant increase in interest in induced seismicity by the scientific community, industry, government, as well as general society. The growing interest in this topic is due to the growing demand for energy and mineral resources [4]. To meet this demand, it is necessary to reach for (increasingly difficult to access) deposits. Complicated and complex mining conditions contribute to a greater number of induced shocks, which cause terrain deformations that, due to their usually unexpected nature, are difficult to detect with traditional surveying methods. Developments in space observation techniques make it possible, for the first time, to record the state of the surface before and after a seismic event and thus determine the resulting change. The process of the development of terrain deformations (the authors use the term as subsidence of the ground's surface) resulting from seismicity induced by underground mining and characteristics (shape) of these deformations are not yet thoroughly investigated. The state of research in this area requires development to better understand the ground surface displacement process and to describe the geometry of the resulting terrain deformations.
Thus, the primary subject of this research is focused on spatial analysis of the geometry of terrain deformation area caused by induced shocks (tremors) on the case study of underground copper mining in Poland. In addition, the relationship of the terrain deformation characteristics (vertical displacement, extent) and the energy of induced shocks has been analyzed. For this purpose, Sentinel-1A/B data for eight different seismic events, differential radar satellite interferometry (DInSAR), and geographic information systems (GIS) have been used. The DInSAR method has the advantage over other geodetic techniques e.g., precise leveling, Global Navigation Satellite Systems (GNSS), as through continuous observation of earth surface at a given interval, it enables to use data recorded prior to unexpected in nature seismic events. Whereas, spatial statistic functions in GIS allow to combine and analyze statistically numerous datasets related to a seismic event.
Further description of the research in this paper has the following structure: presentation of the state of the knowledge in induced seismicity studies with InSAR methods, description of the study area, characteristics of the data and applied methodology, presentation, and discussion of results.

State of Knowledge
Many scientific papers have been published on the study of surface displacements caused by mining using geodetic and remote sensing methods [5][6][7][8]. However, Satellite Interferometric Synthetic Aperture Radar (InSAR)has become the leading method in this field [9]. InSAR data have been utilized to analyze, among other things, long term and short term deformations caused by mining activity [10][11][12][13][14][15]. The emergence of the European Space Agency's Sentinel-1A/B satellite mission in 2014 enabled universal and open access to images from land and ocean surface observations. In the result, research using alternative to commercial spaceborne satellite radar interferometry data sources has been made possible to determine terrain deformations caused by induced shocks. Noteworthy studies in this field have been presented in Table 1. In Poland, research focuses primarily on the underground exploitation of coal and copper, as induced seismic activity is associated with two mining centers: Upper Silesia Coal Basin and Legnica-Głogów Copper District [18][19][20][21][22]. However, internationally, the main areas of induced seismicity research focus on hydraulic fracturing used to extract gas and oil, e.g., in Oklahoma in the USA [23][24][25]27], and production of geothermal energy, e.g., Germany, and Iceland [16,17]. So far, induced tremors have been studied in order to determine the magnitude and location of terrain deformations, in relation to natural and technological conditions using satellite radar interferometry. Results of these studies [16,[18][19][20][21][22][23][24]26,27] indicate that methods of satellite radar interferometry DInSAR, Small Baseline Subset (SBAS) method are suitable for detecting surface deformations measured in Line of Sight (LOS) of the satellite caused by induced seismicity. Up to now, the research on the geometry of surface deformations developed in the result of seismic events caused by underground mining and based on a large number of data has been limited. The presented studies investigate deformations caused primarily by a particular seismic event.
On the other hand, geographic information systems (GIS) provide tools for managing and processing various spatially referenced datasets and enable integration and modeling of the analyzed data. Apart from simple map overlay, spatial statistics, geodata mining and prediction methods, geographical information systems are used to augment studies of terrain deformations in mining areas. Table 2 shows notable examples of the use of GIS in research related to surface displacements caused by mining.  The publications presented above concern mainly modelling and prediction of long-term surface displacements resulting from past and active mining. Geographic information systems analytical functions have been used for: assessing the impact of mining on the surface and surface infrastructure [30][31][32][33], time-spatial analysis of risks related to deformations [34,35,37], and to support public administration [29]. In addition, attempts were also made to identify and quantify the relationship between the location of terrain deformations and factors related to mining and geological conditions [28,36]. The hybrid (GIS and other) methods proposed in some cases were aimed at preparing preliminary data, which were then processed using GIS analytical functions, e.g., spatial regression methods to obtain the final result. In addition to using GIS to study displacements, this environment has been widely applied to assess the environmental impact of mining [38,39] and mineral resource management [40,41] as a leading analytical support system facilitating studies of mining influence on its surroundings. However, apart from mapping terrain deformations [34,35], GIS has not so far been extensively used to analyze surface displacements resulting from mining induced seismicity.
Thus, in our paper we investigate the characteristics of such terrain deformations based on a large number of datasets (8 seismic events and 37 maps of vertical displacements) and map algebra spatial statistics functions.

Study Area
The Rudna copper ore mine is located in the southwestern part of Poland in the Lower Silesian Voivodeship, north of the city of Polkowice. This underground mine is one of the seven mining fields forming the Legnica-Głogów Copper District. The Rudna mining terrain covers about 78 km 2 . The mines are operated by the KGHM Polish Copper S.A. Company. The general location of the mining area is shown in Figure 1.

Geology and Tectonics
In terms of geology, the mine is located on the Fore-Sudetic Monocline (New Copper Basin), where the richest documented accumulation of copper in Poland has been documented ( Figure 2). The depth of deposition of the copper orebody in the Rudna deposit ranges from 844 m up to 1250 m below the terrain surface. The deposit is made of sandstone ores (approximately 80% of resources), carbonate ores (approximately 15%) and copper-bearing shale constituting only 5% of total deposit mass [43]. The deposit layers have an north western-south eastern (North Western-South Eastern) extent according to the course of the monocline border with the Fore-Sudetic block and lean towards NE at an angle of 1 • to 6 • [44]. The southwestern part of the mining area is intersected by faults running east and north-east belonging to the Biedrzychowa fault line, extending from the Polkowice area [45]. Fault throw to the north-west range from 40 m to 140 m. The second most important dislocation named Paulinowa is located west of the shafts in western part of the mine. Paulinowa dislocation is a trough with a depth of 20 m to 30 m, parallel to the Biedrzychowa fault. Another dislocation covers the fault zone of the "Rudna Główna" in the southern part of the deposit area. It consists of bundles of faults forming strata and ditches with amplitudes ranging from several to 30 m. In addition to fault groups, there are single discontinuous dislocations with amplitudes of up to several meters.

Mining Operation
The Rudna operation is the largest copper ore mine in Europe and one of the largest underground mines of its kind in the world. The mine extracts copper ore from the Rudna deposit, part of the Sieroszowice deposit (15%), half of the Głogów Głęboki-Przemysłowy deposit and, to a small extent, a fragment of the Lubin-Małomice deposit [43]. The Rudna mine is distinguished by the largest thickness reaching up to several meters, with the average thickness of the deposit currently at the level of over 4 m. The mining of the deposit at the Rudna mine is carried out in three mining regions: Main Rudna, Western Rudna, Northern Rudna. The orebody is accessed through 10 shafts with depths of 941 m to 1244 m-3 mining, 4 ventilation, and 3 downhill-material. The mining process is carried out using chamber-pillar operation systems. Rudna is the most at risk of seismic activity among all the copper mines in the area, primarily due to the depth of mining.

Induced Seismicity
Over the years, the frequency of induced seismicity events in the Legnica-Glogów Copper District (LGCD) area has been increasing with the development of copper ore mining. Shocks registered here achieve high energy values (a magnitude above 4.0). These higher energy values of shocks are associated with the occurrence of clearly more durable ceiling rocks formed by anhydrites and dolomites. According to publicly available European Mediterranean Seismological Centre (EMSC) database, several dozen shocks occur annually in the LGCD area. In total, from January 2016 to July 2020, there were about 355 induced shocks [46]. A significant proportion of these events are weak shocks that are not felt by people and do not cause terrain deformation. The strongest shocks reaching magnitudes above 4 occur several times a year and are felt throughout the LGCD area. These high energy shocks are believed to cause surface subsidence over an area of several square kilometers and reaching several to a dozen or so centimeters [19,20]. Since 2016, 27 shocks with a magnitude above 4.0 have been recorded, nine shocks in 2016, six shocks in 2017, again, six shocks in 2018, four shocks in 2019, and, so far, two shocks in 2020.

Materials and Methods
In the analysis of terrain deformations caused by induced seismicity, eight events with magnitudes M ≥ 2.71 (by regional seismological station of the Polish Academy of Sciences) were investigated. These shocks occurred between 29 November 2016, and 7 February 2020, i.e., during the operation of the Sentinel-1A/B satellite missions. Six shocks occurred in the winter period, and two in the summer period. The highest magnitude of 3.75 was associated with two events, one on the 26 December 2017, and the other one 29 January 2019. The lowest magnitude 2.71 was associated with the event of 8 April 2017. The research methodology involved four main stages, which are presented on a scheme in Figure 3. To determine terrain displacements caused by these events satellite imagery of the European Space Agency Sentinel-1A/B mission were used. These twin satellites deliver new images of a given location on the Earth's surface every 6 days. All of the data acquired by these satellites are available on the Alaska Satellite Facility (ASF). Data Search website [47]. For each of the analyzed seismic shocks, pairs of images were selected with a time interval of 6, 12, 18, or 24 days for ascending (73 track) and descending (22 track) orbits, with the first image recorded before the event and the second after. In addition, files specifying the precise position of the satellites in orbit and Shuttle Radar Topography Mission (SRTM) Digital Elevation Model of 30 m resolution were used.
The above input data were then used in the second stage to calculate the terrain displacements with differential satellite radar interferometry (DInSAR). The DInSAR method was developed by [48][49][50] and is widely used for detecting displacements caused by, among other things: landslides [51,52], earthquakes [53][54][55], mining exploitation [56,57], volcano eruption [58,59], and glacier motion [60,61]. This method is a development of the classic InSAR method and enables the determination of relative terrain deformations based on two radar scenes that were acquired by satellites at different times for a given place on the Earth's surface. Satellite measurement is based on sending an electromagnetic wave towards the Earth, which after reflection from the earth returns to the satellite. The shift in the wave phase recorded in subsequent satellite passes indicate displacement towards the satellite's line of sight (LOS) [62].
Calculations were performed in the GMT5SAR program [63] and began with coregistration of image pairs. For each pair, the master image being the scene registered earlier. Each master image in a given pair came from the closest available date before the shock occurred. In the next step, an interferogram was generated, i.e., the phase difference of the electromagnetic waves coiled in 2π cycles, was obtained. Then, the difference between the interferogram and SRTM in radians was calculated. Through this operation, the component responsible for the topography of the examined area was eliminated from the data. To improve the quality of the differential interferogram, data filtration (Gaussian filter) was performed to remove noise. In the final step of this stage unwrapping of the phase was performed, which allows reconstruction of the full signal from the differential interferogram. The phase unwrapping was completed in the SNAPHU program using the minimum cost flow method (MCF) [64]. Based on the results of phase unwrapping, Line of Sight (LOS) displacements of terrain were generated and given the WGS-84 reference system.
The LOS displacement raster has been resampled to a pixel dimension of 20 × 30 m. Using the GMT5SAR program, three components of the view geometry were determined (i.e., look vectors: horizontal north-south and east-west, vertical up-down) to calculate the incidence angle for each pixel representing the area of deformation. Based on the geometrical relationships between the Sentinel-1A/B satellites orbit and the measured point on the Earth's surface, the incidence angle was determined using the Formula (1) [65]: where: θ is incidence angle; look_N is horizontal north-south of look vector; look_E is horizontal east-west of look vector; look_U is vertical up-down of look vector.
The value of the incidence angle for Sentinel-1A/B satellites ranges from 29 • to 46 • . In the next step vertical displacements were calculated using the following Formula (2) [66,67]: where: θ is the incidence angle.
Once vertical displacements for each image pair were calculated, spatial statistics with the use of GIS raster map algebra's cell statistic functions were used to determine shape of the deformation area after each of the induced shocks. The functions calculate a per-cell statistic from multiple rasters, such as maximum, minimum, mean or range values. In our case, the number of input rasters (vertical displacement maps) varied from 2 to 8 depending on the event, and the shape of deformation (vertical displacement at each location) was determined as mean value of all spatially congruent cells in raster maps.
These results were used in the next stage alongside descriptive summaries of vertical displacements to reflect the shape of terrain deformations. Descriptive statistics calculated in the study included: minimum, maximum, and mean vertical displacement values. The results have been presented in the form of graphs, and Table 3 shows the Results part. In addition, the Pearson correlation coefficient test statistic was calculated to assess statistically the potential relationship between the maximum observed value of vertical displacements, extent of deformation, and energy of each event.

Results
The terrain deformations resulting from the induced shocks, which occurred in the period from 29 November 2016 to 7 February 2020, were analyzed. These shocks had energies from 2.9 × 10 6 J (registered as 2.71 magnitude by regional seismological station of the Polish Academy of Sciences) to 3.1 × 10 8 J (registered as 3.75 magnitude by regional seismological station of the Polish Academy of Sciences), location and date of these events is shown in Figure 4. Strong mining tremors in the area of LGOM are registered by the regional seismological station of the Polish Academy of Sciences, which is located several dozen kilometers to the south. In total, 68 pairs of images were processed, and 37 generated interferograms were investigated further (Appendix A). The remaining ones were characterized with low coherence values below 0.30 due to dense vegetation or snow coverage. Between 2 and 8 interferograms were obtained for each seismic shock. The time interval between pairs of images ranged from 6 to 24 days. The largest number of correct interferograms (8) was generated for the shock of 7 December 2017, and the lowest number of correct interferograms (2)  The descriptive statistics for each of the analyzed seismic shocks are presented in Table 3. The individual vertical displacement maps obtained from both ascending and descending orbits were combined into a single vertical displacement map for each seismic event. The resulting maps represent the probable shape of the deformation area. For this purpose map algebra raster functions were used where average values of spatially congruent pixels representing vertical displacements were calculated. The resulting maps of terrain deformations have been presented in Figure 5.
In order to examine more precisely the shape of deformations of the ground surface, the generalized horizontal extent of troughs that developed after each event were superimposed with their centroids as the joint center ( Figure 6). The value of −10 mm was assumed as the limit of the deformation area.
Generally, the resulting deformation areas have oval shapes (five cases), slightly elongated in the N-S direction that differs in extent depending on the energy of the seismic event that caused it. Two events produced almost circular deformation areas. Slight rotation of several deformation areas, greatest in the case of the event from 7 December 2017, can be noticed. The S-W elongation of the 26 December 2017, trough, may have been caused by the influence of the earlier shock that occurred in close vicinity on the 7 December 2017. This will beyond the scope of this paper and will be investigated in a separate study.   In the next stage of analysis, the potential relationships between the energy of tremors induced by underground mining activity and maximum vertical displacements, as well as energy and their spatial ranges in W-E and N-S directions were assessed. For this purpose Pearson statistical correlation coefficient (r) was calculated for three cases, energy vs. maximum mean vertical displacement, energy and W-E extent and energy and N-S extent of deformation areas. These relationships have been shown graphically in Figure 8a-c, respectively. The strongest relationship between the analyzed variables was determined for the pair of energy vs. W-E extent and N-S extent variables (r = 0.79). The relationship of energy vs. vertical displacements is the weakest of all three analyzed pairs but still with strong positive correlation of r = 0.74. The outliers from the trend line mainly concern the shocks: 7 December 2017, and 15 September 2018, which have been characterized with the greatest noise (weakest coherence), as was mentioned before.

Discussion
In our study, satellite radar interferometry DInSAR method and Sentinel-1A/B radar images were used to detect LOS displacements caused by tremors induced by underground mining activity. Additionally, GIS functions were used to calculate vertical displacements and analyze geometry of the resulting deformations. Altogether, 68 interferograms were processed, and 37 interferograms were obtained for time bases from 6 days to 24 days. The remaining 31 interferograms were characterized with low coherence values, which could be caused by the coverage of the area with vegetation, presence of snow cover or the influence of weather conditions on the satellite signal. The greatest number of interferograms was obtained for the months of November, December, and February, ranging from 5 to 8 interferograms for a given seismic event. The smallest numbers of interferograms were generated for the months of January, April, and September, ranging from 2 to 4.
As a result, a multiple data set was obtained for a statistical investigation of the geometry of deformations and an attempt to determine the shape of troughs resulting from tremors (generalized shape). Our study, contrary to other works [19][20][21][68][69][70], focused on determining the mean maximum values of vertical and not LOS displacements based on number of vertical displacement maps for each event.
In addition, we have analyzed a large number (8) of seismic events caused by underground mining.
Among the seismic events analyzed in our study three have been the subject of other research. The shock on 29 November 2016, was studied by [19,21,[68][69][70] using both the DInSAR and/or SBAS methods, but based on a smaller number of interferograms. Malinowska et al. [21,68] determined the value of LOS displacements equal to −90 mm. Whereas, Milczarek [69] in the his study determined this value as equal to −80 mm. In another study [19], the author calculated the cumulative LOS displacements using the SBAS method for seismic events that occurred in 2014-2018. Our results indicate value of maximum vertical displacements of −96 mm and generalized extent of the deformation area of 1.6 km × 2.25 km based on superposition of five vertical displacement maps.
For the two other previously analyzed events, i.e., shocks on 7 December 2017, and 26 December 2017, Hejmanowski et al. [20] determined maximum LOS displacements at approximately −80 and −110 mm, respectively. These events were also studied by Owczarz and Blachowski [71] who determined maximum LOS displacements at −83 mm and −94 mm, respectively. Our present results, based on the superposition of 8 and 5 vertical displacement maps indicate subsidence of −85 and −109 mm, respectively and extent of the deformation area approximately 1.4 km × 1.55 km and 1.85 km × 2.35 km, respectively.
In addition, Owczarz [72] calculated the LOS displacement of deformation area (seismic event of the 29 January 2019) using various SAR data processing software and obtained max values of displacements ranging from −92 to −118 mm.
In the result of our calculations, knowledge database of surface deformations manifesting after seismic events was created using the DInSAR and GIS methods. Processing of a large number of image pair for each seismic event allowed for verification of the results and identification of the geometry of the resulting area of deformation based on pixel-to-pixel relationship between vertical displacement maps for different time intervals and different orbits. Based on the literature [20,21] it was assumed that deformations develop suddenly in a short period of time after the seismic event and thus various temporal baselines were analyzed together, i.e., between 6 and 24 days. All of the analyzed deformation areas occurred within the zone of direct influence of underground mining. This may point to the hypothesis that the seismic events accelerate manifestation of these direct effects (surface subsidence).
We have proposed a methodology to determine vertical, in contrast to LOS, displacements caused by seismic events related to underground mining based on pixel-to-pixel relationship between a number of vertical displacement maps for different time intervals. In addition, spatial extent and shape of the subsidence through can be determined in the same way providing more insight into the phenomenon of ground subsidence after induced seismic activity.

Conclusions
The presented research concerned the determination of surface effects caused by eight shocks induced in the area of underground copper mining. For this purpose, the DInSAR method was used to calculate the resulting surface displacements and GIS to analyze the obtained results. On the basis of the Sentinel-1 A/B data 68 image pairs were processed and 37 interferograms were obtained and further analyzed, which enabled the detection of terrain deformations. This research made it possible to determine the average values of: vertical displacements, spatial ranges in the W-E and N-S directions, the geometry of troughs, as well as the relationship between energy and displacements, spatial ranges. As a result of the research, the following conclusions were formulated: 1.
The Sentinel-1 A/B mission ensures sufficient frequency of delivering new satellite imagery to detect terrain displacements due to induced seismicity in the area of underground copper mining.

2.
The DInSAR method enables the determination of displacements caused by induced shocks with energy from 2.9 × 10 6 J to 3.1 × 10 8 J (magnitude from 2.71 to 3.75).

3.
The greatest number of interferograms was calculated for the winter months of: November December, February, while the smallest number in January, April, and September. Low coherence for these months was due to the area being covered with dense vegetation or snow. 4.
GIS provides tools facilitating further processing, as well as analysis of the results and their visualization. The use of spatial statistics made it possible to develop the characteristics of the geometry of the resulting troughs. 5.
The deepest trough was created after the shock of 29 January 2019 (energy 3.1 × 10 8 J) with mean maximum vertical displacement equal to −119 mm. The shallowest trough was caused by the shock of 8 April 2017 (energy 2.9 × 10 6 J) with an mean maximum vertical displacement of −44 mm. The mean values of the spatial ranges of these troughs ranged from 500 m and 1950 m in the W-E direction and from 600 m to 2700 m in the N-S direction. 6.
Strong positive relationship has been determined between energy of induced shocks and the resulting vertical displacements of terrain, as well as between energy and spatial extents of deformation areas. The highest value of the coefficient of Pearson (r = 0.79) was determined for the energy vs. W-E extent and N-S extent relationship while the lowest (r = 0.76) for energy vs. mean maximum vertical displacements.
In the future, the same researches are planned for other induced shocks in order to create a knowledge base on deformations in the Rudna mine.