Estimation of Mining-Induced Horizontal Strain Tensor of Land Surface Applying InSAR

Abstract

Horizontal strains related to mining-induced subsidence may endanger infrastructure and surface users’ safety. While directional horizontal strains should be well determined, appropriate solutions for a complete assessment of the terrain surface deformation field are still required. As a result, the presented study examined a new method for calculating horizontal strain tensor based on the decomposition of satellite radar interferometry (InSAR) observations into vertical and azimuth look direction (ALD) displacements. Based on a geometric integral model, we tested our method on experimental data before applying it to an underground copper ore mine in Poland. In the case study, the displacement field was determined using the Multi-Temporal InSAR method on Sentinel-1 data. The model data relative error did not exceed 0.02 at σ = ±0.003. For the case study, land subsidence of up to −167 mm and ALD displacements ranging from −110 mm to +62 mm was obtained, whereas the extreme values of horizontal strains ranged from −0.52 mm/m to +0.36 mm/m at σ = ±0.050 mm/m. Our results demonstrate the high accuracy of the method in determining the horizontal strain tensor. As a result, the approach can broaden the assessment of the environmental impact of land subsidence worldwide.

1. Introduction

Land subsidence can be caused by both natural Earth factors and anthropogenic activities [1,2,3,4,5,6,7,8,9]. This type of phenomenon has been observed in various parts of the world for decades. Land subsidence harms structures, both surface and subsurface, and may endanger the safety of surface users [10,11,12]. As a result, it is critical to constantly seek solutions that allow for the phenomenon’s negative environmental effects to be minimised [1].

Land subsidence due to mining is primarily caused by the displacement of elements of the rock mass above the seam and post-extraction void migration to the terrain surface [13]. These displacements cause continuous surface deformations. Subsidence, tilt, curvature, horizontal displacement, and horizontal strain are common parameters used to describe continuous surface deformations [13,14]. However, strain is the most important subsidence parameter for analysing natural and built surface features. While several well-established methods exist for forecasting mine subsidence [15,16,17,18,19,20,21,22,23,24,25,26], predicting horizontal strain is significantly more difficult since the strain is sensitive to variations in surface features and the presence of near-surface landforms [27,28].

Early strain prediction methods relied on simple empirical relationships with mining geometry or other predicted subsidence parameters [27,29,30]. When the surface and overburden are relatively flat and there are no significant near-surface geological components, these empirical methods provide reasonable predictions of regular ground movements. However, due to irregular movements, the prediction of horizontal strain acquired using these methods is frequently exceeded at distinct points.

Recently, statistical methods have been applied to estimate strain distribution based on regular and irregular movements [27,28,31,32]. These approaches assess the likelihood of exceeding strain limits based on previously recorded ground monitoring data. However, these statistical methods are generally limited to regions where the mining geometry and overburden lithology are comparable to the mining sites where the monitoring data was obtained.

Importantly, whereas directional horizontal strains may be accurately estimated using, for example, precise levelling and Global Navigation Satellite Systems (GNSS) measurements [33,34,35,36], the variance in horizontal stress magnitudes is poorly understood. As a result, new strategies for assessing the terrain surface deformation field are still necessary. A better understanding of how the magnitudes fluctuate would be beneficial in formulating mine design solutions to address horizontal stress-related ground control issues.

In recent years, InSAR can effectively measure ground surface displacement. This technique has a high spatial resolution and accuracy and is not dependent on weather conditions [37]. Therefore, InSAR is now commonly used to assess ground movements caused by natural and anthropogenic phenomena such as seismicity [38], volcano eruptions [39], landslides [40], glacier motions [41], mining activity [42], and tunnelling [43]. The displacements of the land surface using InSAR data are first calculated in the line-of-sight (LOS) direction [44,45]. Determining the values of such movements in different directions necessitates the use of additional satellite data acquisition geometry [46] or physics assumptions about the observed phenomenon [47].

The estimation of a 3D displacement field from InSAR data can be carried out by applying three types of methods [48]. The first group uses a single-geometry InSAR, and thus LOS observation, in conjunction with an established deformation model of the studied phenomenon [49]. This approach has been used for building damage assessment [47]. The authors of the mentioned approach utilised InSAR observations and a prediction model, which required additional information about the deformation process, e.g., geological condition, mining situation, and parametrisation of the theoretical model. Unfortunately, such information is often difficult to obtain with sufficient reliability. Aside from single-geometry InSAR data sets, Multiple-Aperture InSAR (MAI) [50] or Offset-Tracking (OT) [51] methods can be used to estimate a 3D displacement field. However, the MAI and OT methods typically yield less accurate results than the single-geometry InSAR data set approach [48]. The second set of techniques is based on multi-pass LOS displacements combined with MAI or OT methods. However, their use in slow surface deformation regions is limited due to the low accuracy of the multi-pass LOS technique [52]. Finally, the third group of methods includes GNSS observations. Due to the difficulty in obtaining archival or current GNSS data sets and relatively sparse distribution of GNSS observation points in many parts of the world, their use, however, in determining a 3D displacement field in conjunction with InSAR data is highly limited [48].

When observing natural or anthropogenic phenomena, it is necessary to understand a 3D displacement field related to their occurrence and deformation field components such as the horizontal strain tensor and the rotation tensor [42]. Since surface and underground infrastructure is particularly vulnerable to high horizontal strain tensor values, it is critical to recognise their extreme values for risk assessment and mitigation plan implementation in prone areas [47,53]. Furthermore, the impact of high horizontal strain values on infrastructure and surface user safety is particularly important when underground mining is carried out in densely populated urban areas [54]. As a result, numerous approaches have been proposed to address this research problem. One of them used a strain model in conjunction with the MAI technique to determine the displacement field caused by underground mining applying InSAR double geometry [52]. The size of the calculation window determines the accuracy of the component determination in this solution. The greater the window size, the lower the spatial resolution of the data and, as a result, the accuracy of the modelling results. InSAR observations and theoretical model calculations were also used to determine the horizontal strain tensor for a subsidence trough over a coal longwall [47]. Similarly, an InSAR-based damage assessment method was proposed that utilises horizontal strain tensor obtained from InSAR measurements [55]. The presented results, however, are limited to the ground inclination (slope). As a result, there is a need to develop a high accuracy solution for the horizontal strain tensor determination not constrained by strict geometrical limitations or additional model assumptions.

The presented study fills a gap in the current knowledge by proposing a new method for determining horizontal strain tensors based on double-geometry InSAR. The third component of the displacement field, which is needed to calculate the strain tensor, is calculated by assuming that the horizontal displacement field is proportional to the first derivative of the vertical displacement field. This assumption is also used in the theoretical modelling of mining subsidence displacement fields [56]. Moreover, a similar solution is used in gravity-driven phenomena such as landslides or ice velocity, where the first derivative of the land surface elevation is utilised [57,58]. The proposed method is entirely based on the interpretation of InSAR LOS movements. Additionally, geostatistical tools were used to filter InSAR observations for noise reduction [59]. As a result, the derivative of a displacement field can be reliably computed, and therefore, the horizontal strain tensor is accurately determined.

2. Study Area and Input Data

The presented study was conducted on experimental data first and then applied to real-world data.

The research was carried out in underground copper ore mining areas in Legnica-Głogów Copper District (LGOM), Poland. The study focused on a 55 km² area that included sections of the mining sites of “Sieroszowice”, “Głogów Głęboki-Przemysłowy”, and “Gaworzyce” (Figure 1).

The geological and tectonic conditions in the study area are complicated (Figure 1B). The LGOM is located in the Fore-Sudetic Monocline. The Proterozoic crystalline rocks and Carboniferous sandstones form the foundation of this formation [60]. These strata are overlain by Paleozoic sediments that dip down to the northeast at up to 6°. A deep series of Tertiary sediments represent the Rotliegend and Zechstein strata. Sediments from the Neogene and Paleogene periods formed the youngest rock strata [61]. There are numerous fault zones in the study area (Figure 1B). These zones run northwestern to the southeast. The main faults range in inclination from 30° to 90° [61,62,63].

Since the 1980s, copper ore extraction has been carried out in the study area using the pillar-and-chamber method with bending of the roof layer, partially with surface protection [64,65]. As a result, maximum values of land subsidence recorded in the study area reached a few meters. Notably, the entire subsidence trough usually occurs within 5–9 years of the end of mining operations [15,66]. This implies that the expected dynamics of land subsidence is relatively fast, with a maximum of about 30 mm/month [36]. Furthermore, previous mining operations have revealed the risk of high-energy mining tremors occurring in the study area. Every year, at least a few earthquakes with energies greater than 10E7 J are recorded [61,67,68,69,70,71]. Moreover, mining-induced earthquakes in the study area caused much greater velocity and quantity of ground displacements than land subsidence caused directly by the extraction of copper ore deposits [67,72,73,74,75,76,77,78,79,80,81,82,83].

In a computational experiment, a geometric-integral model based on the Knothe-Budryk theory was used [15,84,85]. The displacement field theoretical values for a rectangular disturbance located at a depth of 800 m with dimensions of 1500 m × 900 m were calculated to reflect local mining and geological conditions. The computations assumed a layer thickness of 2.8 m and two model parameters: 0.70 for the influence factor and 1.5 for the layer strength parameter. These parameters, according to a geometric-integral model, allowed to determine the influence radius $r$. The numerical calculations were carried out on a grid of 50 m × 50 m, with the use of SubCom software (v1.0., manufactured by Wojciech Witkowski and Ryszard Hejmanowski, AGH University of Science and Technology, Cracow, Poland) [86] (Figure 2).

Ground movements were analysed in the presented study using data from Sentinel-1 satellites [87]. The Sentinel-1 acquisitions were obtained between November 2017 and April 2018. The dataset contained 110 ascending geometry interferograms and 86 descending geometry interferograms. The Multi-Temporal InSAR (MT-InSAR) algorithm based on the Small-Baseline (SB) approach was used to process the data (Figure 3) [88]. The MT-InSAR analysis is a type of advanced InSAR data processing technique developed to overcome the limitations of the classical differential InSAR method, such as temporal and geometrical decorrelation, as well as to compensate for errors caused by atmospheric distortions, inaccurate terrain models, and uncertain satellite orbits [89,90,91]. This method’s ability to detect and investigate deformation phenomena has already been demonstrated in various applications [92,93,94,95,96,97,98,99,100]. Furthermore, since our study area is mostly covered by arable fields and forests (Figure 1C), we could retrieve reliable results of ground movements by using the MT-InSAR approach, which can be successfully applied for agricultural and non-urban areas [91].

Interferometric processing was done in Sentinel Application Platform (SNAP v8.0.0, designed by European Space Agency, Paris, France) and StaMPS (v4.1b1, designed by Andy Hooper, University of Leeds, Leeds, UK), whereas modelling and postprocessing of obtained results were carried out in Matlab (vR2021a, designed by Cleve Moler) and ArcGIS (v10.7.1, designed by ESRI, West Redlans, CA, USA). Since the primary goal of the presented research was not to process InSAR data, a more detailed description of MT-InSAR processing was omitted from the article. However, detailed descriptions of these techniques can be found in numerous research papers [88,101,102].

3. Research Methodology

The developed methodology retrieves the results of ascending and descending double geometry InSAR acquisitions (Figure 4) [45]. Based on these observations, the displacement field is decomposed according to Equation (1):

$$\left[\begin{array}{c}{d}_{los}^{asc}\\ d_{los}^{desc}\end{array}\right]=\left[\begin{array}{cc}\cos {\theta }^{asc}& \frac{\sin {\theta }^{asc}}{cos \Delta \alpha }\\ \cos {\theta }^{desc}& \sin {\theta }^{desc}\end{array}\right]\left[\begin{array}{c}{d}_v\\ d_{ALD}\end{array}\right],$$

(1)

where $d_{los}^{asc}$ and $d_{los}^{asc}$ are LOS displacements obtained from ascending and descending geometry mode of InSAR, respectively; ${\theta }^{asc}$ and ${\theta }^{desc}$ are incident angles of ascending and descending InSAR acquisitions, respectively; $\Delta \alpha$ is the satellite heading difference between ascending and descending acquisitions of InSAR data; $d_v$is the vertical displacement of the deformation field and $d_{ALD}$ is the horizontal displacement of the deformation field in the Azimuth Look Direction.

The horizontal components of the displacement field were determined using the concept of gravity-driven processes. Horizontal displacements are assumed to be proportional to the first derivative of vertical displacements in this principle. According to Equation (2), $d_{ALD}$ values in our solution are functionally related to $\frac{\partial d_v}{\partial ALD}$ by parameter $B$:

$$d_{ALD}=B\frac{\partial d_v}{\partial ALD}.$$

(2)

Using the determined parameter $B$, displacements in the directions of the $x$ and $y$ axes can be calculated using Equations (3) and (4), respectively:

$$d_x=B\frac{\partial d_v}{\partial x},$$

(3)

$$d_y=B\frac{\partial d_v}{\partial y}.$$

(4)

The gradient of the displacement field, however, can be used to compute a horizontal deformation tensor $F$. Because of the information about the plane state of displacements, the considerations can be limited to four components of the matrix $F$, as shown in Equation (5):

$$F=T_{\epsilon }+R_{\epsilon }=\left[\begin{array}{cc}\frac{\partial d_x}{\partial x}& \frac{1}{2}\left(\frac{\partial d_x}{\partial y}+\frac{\partial d_y}{\partial x}\right)\\ \frac{1}{2}\left(\frac{\partial d_y}{\partial x}+\frac{\partial d_x}{\partial y}\right)& \frac{\partial d_y}{\partial y}\end{array}\right]+\left[\begin{array}{cc}0& \frac{1}{2}\left(\frac{\partial d_x}{\partial y}-\frac{\partial d_y}{\partial x}\right)\\ \frac{1}{2}\left(\frac{\partial d_y}{\partial x}-\frac{\partial d_x}{\partial y}\right)& 0\end{array}\right],$$

(5)

where $T_{\epsilon }$ is the matrix of the strain tensor and $R_{\epsilon }$ is the antisymmetric matrix of the rigid body rotation tensor. Since rotations have low importance in the displacement field, the analyses were limited to strain tensor $T_{\epsilon }$.

In addition, horizontal displacements $d_{ALD}$ derived from LOS observations are used in the calculations. By determining the directional derivative, the directional horizontal strain value ${\epsilon }_{ALD}$ is obtained. After modifying Equation (5) by entering the symbols ${\epsilon }_{xx}=\frac{\partial d_x}{\partial x}$, ${\epsilon }_{yy}=\frac{\partial d_y}{\partial y}$, ${\gamma }_{xy}=\frac{1}{2}\left(\frac{\partial d_x}{\partial y}+\frac{\partial d_y}{\partial x}\right)$, the tensor component ${\epsilon }_{ALD}$ is determined using Equation (6):

$${\epsilon }_{ALD}={\epsilon }_{xx}co{s}^2\left(ALD\right)+2{\gamma }_{xy}sin\left(ALD\right)cos\left(ALD\right)+{\epsilon }_{yy}si{n}^2\left(ALD\right).$$

(6)

Equation (7) can be used to calculate the extreme values of the horizontal strain tensor after determining the components of the $T_{\epsilon }$ matrix:

$$\left.\begin{array}{c}{\mathsf{\epsilon }}_{\max }\\ {\mathsf{\epsilon }}_{\min }\end{array}\right\}=\frac{{\mathsf{\epsilon }}_{\mathrm{xx}}+{\mathsf{\epsilon }}_{\mathrm{yy}}}{2}\pm \frac{1}{2}\sqrt{{\left({\mathsf{\epsilon }}_{\mathrm{xx}}+{\mathsf{\epsilon }}_{\mathrm{yy}}\right)}^2+4{\mathsf{\gamma }}_{\mathrm{xy}}^2}.$$

(7)

Empirical observations are used to calculate derivatives of the displacement field. As a result, the information obtained is susceptible to observational noise. Data filtering tools were included in the method presented to reduce the effect of noise on empirical data. The geostatistical filtering method based on ordinary kriging was used. In the calculations, the semivariance function $\gamma \left(h\right)$ expressed by Equation (8) was used. The value differences ($z_i$, $z_{i+h}$) between the points in the space are counted by this function. All pairs of points $N\left(h\right)$ separated by the vector $h$ are analysed in the calculations. Semivariance analysis enables the solution of kriging equations and the reduction of system variance. As a result, it is possible to obtain quasi-smooth displacement field distributions and data standard deviation information.

$$\gamma \left(h\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^{N\left(h\right)}{\left(z_i-z_{i+h}\right)}^2.$$

(8)

4. Results and Discussion

4.1. Experimental Data

For the experimental data (see Section 2), the calculated maximum land subsidence is 1.892 m, which occurs in the centre of the analysed deformation field (Figure 5A). Land subsidence causes horizontal strains ranging from −2.00 mm/m to +1.82 mm/m. The asymmetric distribution of these values results from the accumulation of compression zones in the centre of the analysed field. By examining the spatial distribution of the calculated horizontal strains, we can conclude that tension zones are located outside the deformation field boundary, whereas compression zones are located within the deformation field. Furthermore, the extreme horizontal strain values occur near the boundary vertices of the field (Figure 5B).

The proposed methodology is based on two parameters derived from the decomposition vector of the displacement field. As a result, the model values of the vertical axis movements $d_v$ and the calculated directional horizontal displacement $d_{ALD}$ were used to determine the components of the strain tensor $T_{\epsilon }$. These components were obtained at each point of the computational grid. Finally, the extreme horizontal strain values (Figure 5C) were calculated and compared to the model values (Figure 5D). The average difference between the model and the calculated extreme deformation values was equal to 0.00 mm/m, with σ = ±0.007 mm/m. For the model calculations, the relative error did not exceed 0.02, with σ = ±0.003.

The experiment results demonstrated that the proposed method is correct and accurate. However, the accuracy of the results depends on the quality of the InSAR data processing and the real-life geological and mining conditions in subsidence areas (e.g., anomaly in terrain ground deformation). As a result, if the above assumptions are met, our approach can successfully assess the land surface deformation fields associated with processes characterised by gravity-driven mechanisms.

4.2. Real-World Case Study

Based on InSAR observations, six subsidence troughs were detected in the case study area. The following analyses were carried out for centroids determined for common PS points from two geometries (ascending and descending). As a result, we obtained the density of the observations of approximately 232 points/km². The maximum vertical displacements reached −167 mm with σ = ±$d_{ALD}$ ranged from −110 mm to +62 mm with σ = ±8 mm. The detected land subsidence has a symmetrical distribution, with the highest values located in the centre of each subsidence trough (Figure 6A). The least visible trough is in the study area’s $d_{ALD}$ displacements (Figure 6B). The negative horizontal displacement values reflect westward movement while the positive values show movement into the east. Thus, the eastern parts of each trough have negative horizontal displacement values (displacements towards the centre of land subsidence trough), whereas the western parts have positive displacement values (displacements also towards the centre of land subsidence trough). This demonstrates that the land subsidence in the case study is a gravity-driven process.

Based on the InSAR-retrieved ground movements in the case study area, Equation (2) was used to determine the relationship between the horizontal directional displacement and the first derivative of the vertical displacements (Figure 7). As a result, the proportionality coefficient $B$ equal to −308.9 m was obtained. The outcome of the calculations is characterised by a strong correlation coefficient R² of 0.805. In general, the region with the greatest absolute horizontal displacements corresponds to the highest values of the first derivative of vertical displacement (Figure 7).

The remaining distributions for $d_x, d_y,\frac{\partial d_x}{\partial y},\frac{\partial d_y}{\partial y}$ and $\frac{\partial d_y}{\partial x}$ values were determined using the $B$ parameter. Geostatistical filtering was used to denoise the data before the derivative calculation stage. The most reliable theoretical model for semivariogram function was the Gaussian model with nugget effect in all analysed data. The established relationship between the directional strain and the horizontal tensor strain components allowed the displacement value in $d_{ALD}$ to be used and the last strain tensor element $\frac{\partial d_x}{\partial x}$ related to the $x$ axis to be calculated. As a result, the horizontal strain tensor for each site point determined by the MT-InSAR method was obtained. Moreover, the spatial resolution of interferometric analysis data was preserved.

The final stage of the calculations was to determine the maximum and minimum values of horizontal strains in the case study area (Figure 8). Based on these findings, it is possible to conclude that the centres of land subsidence troughs have extreme compressive strains of −0.52 mm/m. In addition, extreme tensile strains of up to +0.36 mm/m occur at the trough’s edge. The asymmetry of these values is related to the accuracy of LOS observations and the physics of the observed phenomenon.

The accuracy of the results was calculated in the potentially stable region, outside the subsidence troughs, considering the variability of these observations. The −10 mm contour was used to establish the border of the land subsidence troughs initially. This value is defined by the vertical displacement accuracy at the level of 2σ of determined vertical displacements. For the following calculations, points from outside the observed troughs were chosen (Figure 8). The average horizontal strain value for these selected points was +0.01 mm/m, with σ = ±0.050 mm/m. Therefore, assuming a 95% confidence level, observations with an accuracy greater than 2σ, namely ±0.100 mm/m, are considered reliable.

5. Conclusions

Land subsidence is a common occurrence worldwide, causing infrastructure damage and posing a direct threat to residents in vulnerable areas. However, since the rock mass and terrain surface deformation process is a complicated research subject reliant on numerous geological, hydrogeological, and geomechanical variables, the spatial extent and values of the phenomenon, and thus its environmental impact, vary significantly.

There are numerous methods for estimating ground surface movements in the case of mining-induced land subsidence. Nevertheless, the horizontal strains that appear in mining areas may also cause significant damage to surface objects. So far, the accurate determination of these values remains a difficult scientific problem. As a result, solutions that allow for the effective determination of the full deformation field of the terrain surface are required. Therefore, the goal of the presented study was to propose a novel method for the reliable estimation of horizontal strains related to the occurrence of land subsidence. Based on double-geometry InSAR and the assumption that the horizontal displacement field is proportional to the first derivative of the vertical displacement field, this method was developed.

It should be noted that in the real-world case study presented in this study, the land surface displacements determined using MT-InSAR were relatively small, ranging from a few milimetres to a few centimetres. As a result, the deformation values of the terrain surface were also small, rarely exceeding 1 mm/m. Nonetheless, the computation results for the real-world case study confirm that the proposed algorithm can determine reliable modelling results with high accuracy even when horizontal strains are small.

Furthermore, since InSAR observations cover large parts of the Earth’s surface and modern satellites, including Sentinel-1, and are characterised by a short revisit time, the method can be successfully applied to gravity-driven phenomena worldwide. As a result, the research presented here may contribute to increased safety in areas of land subsidence caused by various natural and anthropogenic phenomena around the world.