Application of a Small Baseline Subset Time Series Method with Atmospheric Correction in Monitoring Results of Mining Activity on Ground Surface and in Detecting Induced Seismic Events

The article presents the results of applying a Small Baseline Subset (SBAS) method to the monitoring of the influence of underground mining exploitation on terrain surface. Calculations were performed on the basis of two large SAR data sets. Calculations of time series additionally employed the Tymofyeyeva and Fialko empirical method for estimating atmospheric delay. A series of test calculations allowed verifying the potential of the method. The SBAS calculations were based on theoretical data prepared with the use of the Mogi model. This approach allowed simulating both underground mining exploitation and induced tremors. The results indicate that the model is useful in detecting sudden ground surface deformations. In such a situation, the image of a tremor can be observed in the atmospheric component. The actual data used in the analysis covered a region of underground copper ore extraction in southwest Poland. The calculations were based on SAR data obtained from the Sentinel 1A/B satellites, for the period between the end of 2014 and May 2018, from paths 22 and 73. In both cases, the total number of images exceeded 120. They served to calculate over 420 interferograms for each set. The results of ground surface displacements allowed precise identification of areas affected by underground mining activity. Over the acquisition period of six days, sufficient amount of SAR data was obtained to precisely monitor surface changes in the analyzed area. The employed empirical atmospheric delay reduction method was demonstrated to enable the detection of sudden ground surface changes due to mining tremors. This feature makes the method useful in detecting induced seismic events in areas characterized by intensive mining activity. It allows determining tremor locations in cases when acquisition dates are greater than the optimal value of 6–12 days (in the case of Sentinel 1A/B).


Introduction
Ground surface deformations due to underground mining activity belong to the most important results of some human interventions in the rock mass. Depending on the geological conditions, they may have a continuous or non-continuous character [1]. Ground surface displacements may be particularly dangerous in urban areas, where they may pose a threat to both the ground-level and underground infrastructure. Therefore, areas affected by mining operations are constantly monitored in order to provide up-to-date information on the scale and range of the observed deformations. The measurement results also serve to prepare and verify forecasts of how underground operations influence changes on the ground surface. The methods most frequently used to measure ground surface deformations in mining areas include: precise leveling lines [2,3], GNSS measurements [4],

Methods
A general equation of the interferometric phase ∆ϕ calculated from two SAR acquisitions (t i and t j ) may be shown as follows: where: λ represents the carrier wavelength of the radar, ∆r ij disp is the ground displacement in radar LOS direction, ∆r ij geom repre sents the geometrical range difference from radar to the target caused by a nonzero spatial baseline between the two orbits, and ∆r ij atm is the atmospheric delay, respectively. As is the case with other methods for measuring deformations of ground surfaces or engineering structures, the currently used InSAR methods have significant limitations [20]. Apart from the unknown deformation component, a representative interferogram also includes other components which significantly influence the resulting value of the interferometric phase. The geometric part of the signal ∆r ij geom , which includes mainly orbit errors and topographic errors, may be relatively easy to separate. On the other hand, the atmospheric component (∆r ij atm ), composed of ionospheric delay (ionospheric screen) and tropospheric delay, is one of the main areas in current research focused on increasing the accuracy of ground surface deformation calculations performed with the use of InSAR methods [21,22]. In order to determine the values and subsequently to eliminate the influence of errors related to the passing of the signal through the ionosphere and the troposphere, an iterative method was used, as proposed by Tymofyeyev and Fialko [19]. The main principle behind the method is that two (or more) interferograms having the same scene contain an identical (with respect to the value) atmospheric component of the phase due to wave propagation in the troposphere and in the ionosphere. Such an approach for determining the atmospheric component is empirical. The theoretical foundations of the methods are presented below.
Let us assume that we have two interferograms calculated on the basis of three images (in such case the acquisition times are t 1 , t 2 and t 3 , respectively- Figure 1). In the case of each interferogram, the equations will have the following form, respectively: where α i -atmospheric component due to signal propagation among others in the atmosphere and in the ionosphere, in time t i and ∆τ ij is the phase component related to the observed ground surface deformations in time t ij . Figure 1. Graphic representation of the basic assumptions behind the method. Part A shows two interferometric pairs having one common image (number 2, denoted in red color). Part B represents the idea behind the method: the more pairs of SAR images with the same recurring image (in red), the more accurate the estimation of atmospheric correction for the image (date). Part C shows the baseline plot for the analyzed SAR data from path 73. Part D, on the other hand, shows the actual implementation of the model. Parts A and B are based on [19].
At this stage, an additional assumption should be formulated that the ∆φ 12 − ∆φ 23 difference is independent of the deformation component on condition that the time base is constant and deformation increments are constant in time. This assumption is of key importance for the results presented in this paper. It is also similar to the properties of other empirical methods employed to eliminate atmospheric delays in InSAR calculations [23].
The accuracy of the determined atmospheric component of the phase increases together with the number of interferograms for larger aperture of the averaging stencil N. The general equation describing atmospheric component has the following form: The above-described method has been implemented among others in the SBAS time series method. It was successfully employed in seismically active regions of California and Mexico [19,24].
The method was verified with three test models. The data were prepared on the basis of the Mogi model [25], which was developed to test ground deformations in the vicinity of volcanoes.
However, the universal character of the model allows a number of further applications [26][27][28][29]. The Mogi model ( Figure 2) is an analytical model, which describes the magnitude and the direction of ground displacements due to changed volume or pressure of a point in a uniform elastic half space. The displacements of any point on the surface can be calculated from the following equation: The main purpose of using the Mogi model, as well as the test calculations was to check how the empirical model of atmospheric delay estimation [19] will affect the results of displacement calculations in the case of an induced seismic event. Taking into consideration one of the assumptions of the model stating that displacements are constant in time, it was checked whether in the case of a mining tremor, it will be possible to extract the signal caused by this tremor. If so, the model proposed by Tymofyeyev and Fialko could be used to determine the impact of mining tremors on the surface indirectly. Detection of the impact of mining tremors in real conditions based on SAR data is often limited. It arises from the fact that intensive and large scale underground mining (e.g., the LGOM area) causes numerous local areas of subsidence-the so-called subsidence basins. In the majority of known cases, the impact of mining tremor is observed within the limits of existing subsidence basins. This is an additional difficulty in determining the impact of mining tremors on the surface.

Area of Study
Copper ore is mined in the southwest part of Poland. The mining areas having total surface exceeding 750 km 2 are located in the Fore Sudetic monocline, whose structure was shaped during the Cimmerian and Laramie movements. The copper ore deposit in this region is located in the roof part of the Rotliegend and in the floor part of the Zechstein. The ore is currently extracted at a depth of over 1100 m below the terrain surface. The LGCB comprises seven mining areas ( Figure 3) operated by three mining plants. Copper ore is mined using the room and pillar system. In the case of the LGCB, the level of continuous ground surface deformations is consistent with the range of a particular deposit within an area of the copper belt. The scale of deformations is mostly influenced by the thickness of the mined rock mass and by the method applied to the liquidation of the mined-out areas. Induced tremors observed in the region are related to the type of rock in the roof layers of the mined deposit. Every year, the region suffers from several seismic events exceeding magnitude Mw 4.0. The simultaneous occurrence of mining-related deformations, mining tremors and deformations due to rock mass drainage makes the LGCB region a valuable area of research into the influence of mining activity on the surface.

Data Time Series Processing
The calculations of the satellite radar data were based on the SBAS method [30][31][32] and with simultaneous allowance for the method of reducing the influence of the atmosphere [19]. As the results could not be compared with the results from other measurement methods (GNSS, precise leveling), the calculations were performed independently for both paths. The analyzed period spanned the years from 2014 to 2018 (Table 1).

Test Calculations: Application of the Mogi Model
The aim of the test calculations was to verify the empirical solution for the reduction of atmospheric disturbances in the phase signal of the interferogram, as suggested by Fialko. In order to simulate ground surface subsidence due to underground mining activity, the Mogi model was employed, as it allowed simulations of mining at a set depth. 24 time periods were assumed (t 1 . . . t 24 ) between which time-constant ground surface subsidences of −40 mm occurred due to mining activity at a depth of 500 m. The time interval between the periods was constant. Figure 4 shows baseline plot for the test data. The spatial range of the test data was 10 × 10 km. Based on the defined assumptions, three independent calculation cases were prepared: mining activity without a mining tremor, (II) mining activity with a mining tremor recorded between t 18 and t 19 , the tremor having occurred exactly in the center of the already existing subsidence trough, (III) mining activity with a mining tremor recorded between t 18 and t 19 , but unlike in case (II) the tremor occurred outside the mining area influenced by mining operations and was located at a distance of approximately 3 km from the center of the already existing subsidence trough.

Calculations for the Lgcb Region
The calculations of the interferograms were based on two-pass processing with the use of GMTSAR ver. 5.4.4 [70]. Phase unwrapping was performed the with Snaphu ver. 1.4.2 application [71]. Wave phase correction in relation to the ground surface was performed with the data from Shuttle Radar Topography Mission (SRTM) version 1 [72]. The enhanced spectral diversity (ESD) algorithm was not used at the preprocessing stage. This algorithm may remove some part of the atmospheric component in the burst boundaries [24]. Basic information on the dataset for the analyzed paths is provided in Table 2. Figure 5 shows the calculated interferometric pairs.

Test Calculations-Application of the Mogi Model
In the first case, without a mining tremor, maximum ground subsidence was −969 mm. Figure 6 shows subsidence increments along a selected profile and allows an observation that increments are constant between subsequent periods. As is expected, no significant values of atmospheric component are observed.
In the second analyzed case (Figure 7), a tremor was assumed in the period between t 18 and t 19 . In this scenario, total subsidence was −1150 mm. The graph (Figure 7, right-hand side, in the middle), shows the values of the determined atmospheric corrections. Corrections having maximum values between −90 and +75 mm were observed in the periods directly before and after the tremor (t 18 and  The third case also involved a mining tremor between periods t 18 and t 19 . While in the previous case the tremor was located exactly in the already existing subsidence trough, the tremor in this scenario occurred outside the surface area influenced by mining activity (Figure 8). The tests proved that allowance for atmospheric delay reduction in SBAS calculations does not affect or distort the calculated displacement results. However, the test calculations demonstrated that in the case of a sudden event (e.g. a seismic event), the image of the event is to some an extent recorded in the atmospheric component of the signal.

Calculations for the Lgcb Region
SBAS calculations of the data obtained from the Sentinel 1A/B satellites allowed identification of 27 regions in which underground mining activity influences the ground surface. Total displacements for paths 22 and 73 ( Figure 9) without and with allowance for atmospheric correction are insufficient to estimate the influence of the model on the final results. Therefore, a more detailed analysis of the obtained results is offered in the Discussion section.
Graphs representing ground subsidence were also plotted for a selected area ( Figure 10). The results for the complete period subjected to analysis clearly demonstrate that the influence of underground mining on the surface is constant. This fact is of key importance in the context of the method for atmospheric delay reduction used in the calculations. Analysis of the displacements indicates regions in which successive curves intersect. In some locations, especially those in which the observed subsidences are small-scale, subsidence increments are inconsistent. In the majority of the graphs, however, successive curves have similar shapes, with a certain offset. This observation confirms the quality of the results obtained for the analyzed area. It also proves that SAR data and the method are useful in monitoring this type of phenomena on the ground surface.

Discussion
Section 5.2 presents the calculation results for ground displacements in mining areas. The results allow an observation that SAR data (in particular data obtained from new satellites) may successfully supplement information gathered with the help of methods commonly used to measure ground surface deformations. First, however, it was necessary to determine the degree to which the Fialko method influences final results in the SBAS calculations. Figure 11 provides the results of LOS displacements for selected points with and without allowance for atmospheric correction. Based on the graphs, a conclusion may be offered that the obtained displacement results (without allowance for atmospheric correction) show limited disturbances. The results for the points located in selected subsidence troughs demonstrate a clear trend downwards, due to the progress of mining works. No points were observed to have outlying values. As a result, the values of displacements with allowance for atmospheric correction are similar to the original values. The differences between the results reach more than ten millimeters, but are seldom greater than several tens of millimeters. Significantly higher values are observed in areas not affected by mining activities (points 5 and 6, Figure 11).
In order to compare displacements from path 22 and from path 73, 10 regions were selected to plot graphs of total displacements ( Figure 12).
Analysis of the LOS displacements from the two paths for 10 representative profiles indicates that the majority of the results are consistent (Figure 13). The results are significantly different only for two paths on profiles 3 and 10. In the remaining cases, the displacements have similar curves and maximum values.
The results of the test calculations are also presented. These calculations were intended to verify whether the model for reducing the influence of atmospheric delay on final results may be considered useful. This empirical method has been successfully used in areas demonstrating displacements which are constant in time [24]. Like some other empirical methods [23], this method does not work well in the case of sudden events (outlying events), such as mining tremors. Figure 14 shows an example of how the method is applied in an area, in which time-constant displacements are accompanied by sporadic sudden events. The results here presented relate to the test models described above. Time-constant ground displacements determined with the SBAS method without the model are correct (i.e., they show as time-constant subsidence increments). The situation is different in the case when an event results in greater displacements observed on the surface. In such case the model reduces the displacement values to approximate the level from the periods when surface displacements were constant. As a result, displacement analysis cannot be used to identify tremor location. In such situation, the results are not fully reliable. The time of the tremor may be identified only by analyzing the atmospheric component. Interestingly, upon comparison, the calculated total displacement values show to be practically identical for cases in which tremors are analyzed with and without allowance for the reduction model. Figure 11. Comparison of the results of LOS calculations without and with allowance for atmospheric correction (black dots and red dots, respectively). Locations of selected points are shown in Figure 9. This fact allowed a hypothesis that the method in question can be applied to detecting induced seismic events. The model seems applicable in the case of such regions as the LGCB. Intensive, albeit time-constant underground mining activity results in the creation of subsidence troughs, while specific geological conditions may contribute to mining tremors. Tremors of this type are located in the mining area, i.e., in the area affected by ground surface subsidence. With all of these conditions overlapping, using InSAR against the background of mining deformations to detect a potential mining tremor may prove difficult or impossible. The test calculations demonstrated that the model allows the detection of a mining tremor in a separated atmospheric signal. Detailed analysis covered 23 induced seismic events, which occurred between December 2014 and June 2018 in the LGCB region. Their distribution in time with respect to SAR data is shown in Figure 5. Tremors selected for analysis had magnitude greater than Mw 4.0. The selection of tremors for analysis was based on the data obtained from Helmholtz Centre Potsdam. Each tremor was matched with two pairs of SAR images obtained from the Sentinel 1A/B satellites, from paths 22 and 73 (Table 3). Time consistency was ensured by selecting images dated closest to a particular tremor (before and after, at a maximum of 12 days).
The results (Table 3, Figure 15) indicate the possibility to observe an image of a tremor in the calculated atmospheric component. In the case of the LGCB region, the model for reducing atmospheric delays allows detecting the majority of induced seismic events having magnitude Mw greater than 4.1.
The demonstrated property of SBAS calculations that takes into account the empirical model of atmospheric delay estimation, can be used to detect the impact of mining tremors on the surface. Analysis of the results showed that detection of mining tremors based on SAR data obtained for a period longer than 6-12 days is limited due to influence mining operation observed in the phase of the signal. Figure 15 shows the results for three selected induced seismic events. The first two lines concern the same event. The first row shows the results for SAR data with a temporal base of 12 days. The results obtained in this case allow to clearly locate the mining tremor. The second row presents displacement results for longer temporal bases, corresponding to: 36 days (for track 22) and 24 days (for track 73) respectively. In the case of track 22 where the temporal base was over one month, additional subsidence basins can be observed. These are the result of the impact of mining on the surface. The rows 3 and 4 present the results for selected mining tremors that were determined based on SAR data with an interval of 6 days.

Conclusions
Application of satellite interferometry in the monitoring of mining and post-mining areas (secondary deformations) may prove an important aid to the methods currently used for measuring ground surface deformations. This paper presents an analysis of the results of surface deformations due to underground mining of copper ore. The analyses are based on calculations of a large number of SAR data obtained from satellites Sentinel 1A/B. The calculations employed the SBAS method and the empirical method for the elimination of atmospheric delay, as proposed by Tymofyeyev and Fialko [19]. Prior to calculating actual data, a series of test calculations was performed. These calculations employed the Mogi model in order to simulate the influence of underground mining activity on ground surface and on mining tremors. The test calculations performed with the use of the atmospheric delay reduction model resulted in two observations which are of significance from the perspective of monitoring surfaces of mining areas. Firstly, if constant displacements occur in the analyzed area, then the application of the model does not have a significant effect on the scale of calculated displacements. Secondly, if a tremor occurs in a considered area, then the applied model reduces displacements due to the seismic event. Some values of these displacements are observed in the calculated atmospheric correction, related to images forming interferometric pairs. If a tremor occurs outside the area affected by mining activity, the displacement value is practically completely reduced. The image of such an event is recorded in the part of the signal corresponding to the atmospheric correction. Therefore, the method should not be applied in areas where ground surface deformations are not time-constant. The method may be used, however, in such regions as the LGCB, where the employed mining system results in relatively time-constant displacements. The results indicate that the method is effective at detecting changes resulting from induced seismic events. It allows determining tremor locations in cases when acquisition dates are greater than the optimal value of 6-12 days (in the case of Sentinel 1A/B). Importantly, the calculation results presented here have some weaknesses. One such flaw is lack of comparison with the results from other measurement methods (e.g., GNSS, precise leveling). Unfortunately, access to the GNSS and precise leveling data is currently impossible due to the policy implemented by the mines of the LGCB region.