Upper Silesian Coal Basin (USCB) in Poland is one of the most significant coal basins in Europe. Carboniferous hard coal underground exploitation has been conducted here since the XXI Century (Konopko 2010). Nowadays, due to significant reduction of production since 1990 there are only around forty deposits exploited by thirty active hard coal mines, where longwall mining typically operates at depths reaching more than 1000 m (Konopko 2010). At the same time, Upper Silesia is a relatively big metropolitan area occupied by 37 towns with almost 3 million inhabitants. More than a century of mining activity has caused irreversible changes on the earth’s surface over large areas. Extraction of coal is accompanied by various geohazards like ground displacements, seismic tremors, sinkholes, flooding and inundations, slope instabilities, and collapses of spoil heaps, air, soil, and water pollution by toxic waste from mining. Subsidence is one of the most significant threats related to underground mining activities. As a result of coal mining, the rock mass is disturbed, and stress changes occur in the overburden, arising from the displacements on the ground surface. These can be immediate or long-term, related to ground consolidation and dependent on the rate of pressure dissipation. Ground surface deformations in USCB can reach several meters (Borecki 1980) thus subsidence can cause severe damage to building structures and infrastructure in the area of operation. Mining subsidence in the area of USCB was the subject of many researches like: (Budryk 1947; Borecki 1980; Knothe 1957; Knothe et al. 1995; Kwiatek 1997; Popiołek 2009; Konopko 2010; Białek et al. 2014; Kowalski 2015; 2020).

Vertical terrain displacements are monitored using classical geodetic methods and occasionally using differential analyzes of high-resolution numerical terrain models obtained by LIDAR. Levelling is the most popular subsidence monitoring method used in USCB. However, the main limitation of field-based monitoring techniques is the need for frequent field observations, which will increase the monitoring costs. The rapid development of remote sensing Earth observation methods allows them to detect changes in the terrain surface. Among these methods, Synthetic Aperture Radar Interferometry (InSAR) (Massonnet and Feigl 1998) has proved to be the most useful (Carnec et al. 2000; Wegmuller et al. 2004; Herrera et al. 2010; Benecke et al. 2012; Engelbrecht et al. 2013).

Since the beginning of the technique, studies applications of InSAR for monitoring of mining inducted surface subsidence were conducted in many different extraction sites, e.g. (Herrera et al. 2007; Engelbrecht and Inggs 2013; Samsonov et al. 2013; Liu et al. 2014). Nowadays, technology is recognized as a good tool for monitoring mining influences in Europe and the world (Zheng et al. 2018; Falorni et al. 2018; Modeste et al. 2021). USCB has been the subject of similar researches since 1998 (Perski 1998; Perski and Jura 1999). Currently, several research centres in Poland are involved in the studies of InSAR use for mining-induced subsidence monitoring, mentioning Ilieva et al. 2019, Pawłuszek-Filipiak and Borkowski 2020, Sopata et al. 2020, Witkowski et al. 2021, Dwornik et al. 2021 as an example.

The presented work aimed to combine several different InSAR datasets to marge subsidence information. Scientific projects conducted at the PGI-NRI since 2005 (Graniczny et al. 2005) resulted in a collection of substantial InSAR data sets in the USCB area. Over the years, PGI-NRI has carried out a series of InSAR ground subsidence studies in the mining area (Graniczny et al. 2007, 2014a, b; Czarnogórska 2010; Del Ventisette et al. 2013; Przyłucka et al. 2015a), as well as the advantages and disadvantages of using individual SAR sensors to monitor these special geological processes (Przyłucka et al. 2014; Przyłucka and Graniczny 2015; Przyłucka et al. 2015b; Graniczny et al. 2015; Przyłucka 2017).

Among other things, it was found that the PS data from the C range made it possible to determine the boundary of the area subjected to the influence of mine activity but did not provide complete information on subsidence greater than 3 cm per year. This information can be supplemented by differential interferograms, which also allow identification of the development of subsidence troughs. However, a different way of processing and presenting the results makes the joint analysis of the data difficult.

This work aimed to combine the various data sets into one easy-to-read and understand map that shows the subsidence over twenty years in a compiled way. Although the data come from different sensors covering different time and spatial ranges, all carry information about vertical displacements of the terrain related to mining therefore an attempt was made to combine them.

As a result, a map of areas that have been exposed to subsidence within 42 hard coal exploited deposits (3045 square kilometres) was obtained. The map covers the period 1992–2012. A total of 81 interferograms (18 from the archival ERS and ENVISAT C-band scenes, 26 from the processing of ALOS L-band scenes, and 37 from the X-band scenes of TerraSAR-X satellite) were used in the study. The interferograms allowed not only to determine subsidence troughs (basins) formed from 1992 to 2012 but also to observe subsidence development over time in several dozen days long increments. The work also included five sets of PSI processing, covering different time and spatial ranges. Three sets come from ERS satellite (period 1992–2003), one from the Envisat satellite (2003–2010), and one from the TerraSAR-X satellite (2011–2012). The PSI point datasets were used to determine zones of residual subsidence. The information obtained successfully reflects the scope and scale of the underground mining influence on the surface, which occurred in 1992–2012.

The distribution of coal mining in the area of the USCB depends on geological factors: the depth of coal seams below the surface, as well as tectonic and hydrogeological factors. At the end of the 1990s, the mining development of the USCB area comprised 70 areas of mines and reserve fields. Average annual production increased from around 41 million tonnes in 1920, 104 million tonnes in 1960, reaching the maximum value of 193 million tonnes in 1980. In 1996, 65 mines were operating, extracting 130 million tons of coal. In 1997, production was 137.13 million tonnes, and 122.7 million tonnes in 1998, carried out by 54 mines. Currently, the annual exploitation is carried out at the level of 60 million tons. The depth of exploitation is now 600–1300 m, and it is decreasing by about 10 m per year (Kowalski 2020). The hard coal mining in the USCB is carried out in a longwall system, mainly with roof collapse (85%) or hydraulic backfilling (15%). The longwall lengths are 100–200 m, and their progress is on average 40 m per month. The most extended lengths of longwall – 275 m were used in 1991–1993 on Staszic Mine with the progress of up to 200 m per month (Kowalski 2020).

The contemporary distribution of the areas of active exploitation depends on the degree of depletion of the most advantageously located deposits. At present, mining is carried out: in the area of the main saddle within the latitudinal belt Bielszowice - Halemba - Katowice - Murcki - Tychy, and in the area of the southern wing of the Bytom basin in the belt: - Bytom - Piekary Śląskie and in the area of Rybnik - Rydułtowy - Wodzisław – Żory (Fig. 1).

As a result of the exploitation, large areas were permanently transformed. Among geohazards such as sinkholes, inundations, and landslides, the impact of subsidence is the most noticeable. According to the current estimates, the direct influence covers approx. 650 km² and indirect influence of the mining industry approx. 1000 km² (Kowalski 2020). The subsidence over a typical mining wall at a depth of 680 m, thickness of 2.5 m and 250 to 400 m long, reaches up to 70% of the layer height, which corresponds to the largest vertical displacement from 0.75 to 2.0 m in the central part of the basin (Kowalski 2017). Due to the use of many seams, vertical displacements in the area reach up to tens of meters in some places. Mining works are usually concentrated on 3–4 faces per year so that the impact on the surface is as low as possible (Borecki 1980).

Nevertheless, long-term subsidence reaches values of 20–25 m in many regions of USCB. For the greater part of the area, depressions reached 0.5–1.0 m, and in extreme cases, they were even 30–40 m (Kowalski 2017). The most significant increase in subsidence was recorded in the 1970s, during intensive exploitation in the Bytom protection pillar. At that time, the average daily increments of subsidence of up to 5 mm were recorded. The materials published so far show that the centuries-old history of coal mining has left vast subsidence basins with a depth of up to tens of meters (Kowalski et al. 2020).

Synthetic Aperture Radar Interferometry (InSAR) is a remote sensing technique that has revolutionized the field of ground deformation research since its first use about three decades ago (Massonnet et al. 1998). The technique is based on the processing of satellite radar scenes to detect changes on the terrain surface.

In this work, datasets processed in two different approaches were used. The Differential Synthetic Aperature Radar Interferometry (DInSAR), a classical method of creating interferograms, is based on composing images of the same area taken at different times (Ferretti et al. 2007). Interferometric fringes are obtained by calculating the phase difference of the reflected wave between the corresponding pixels from two images. Fringes on the interferogram show the change on the terrain surface in the studied area that happened between the recordings of these images (usually about a month or two). This enables deformation registration from centimetres to decimeters (depending on the wavelength used for the SAR scene).

Further processing of the interferograms produces a vertical displacement image instead of a pure interferogram. Information of the displacement is therefore presented as a raster image. The informative content of such maps is heavily dependent on the radar coherence during the acquisition period. Consequently, land cover change and the atmospheric disturbances limit the measurement precision to centimeters (Massonnet and Feigl 1998).

At the beginning of the XXI century, advanced DInSAR processing techniques relying on multi-interferogram or multi-image algorithms to study slow ground deformation have been developed to overcome some of the limitations related to conventional DInSAR approaches (Ferretti et al. 2000, 2001). Over the years, many different algorithms for processing a stack of SAR scenes have been invented, which together can be called Persistent Scatter Interferometry (PSI). Their list is presented, e.g. in (Crosetto et al. 2016). Through the processing of a large dataset (usually more than 20 images), a selected subset of pixels is used to estimate atmospheric phase screen, allowing measuring surface deformation with millimetre accuracy (Crosetto et al. 2016). As a result of the processing, a set of several hundred thousand points, persistent scatterers (PS), is kept, which reflect the pixels where radar wave has acted coherently throughout the acquisition period. For each pixel selected by the PSI method, the Line of Sight (LOS) deformation time series and the LOS velocity are obtained. The PS method can measure surface motions at a 1 mm/year level and resolve very small-scale motion not previously recognized in classical SAR interferometry (Crosetto et al. 2016).

In summary, the deformation is presented as a dense but not uniform vector set of surface points with assigned average velocity values of movement over the imaging period. This technique makes it possible to obtain information about the surface deformations, as long as they are slow. Points are not generated if the terrain changes faster than approx. 15 cm/year for C-band satellites and 25 cm/year for X-band satellites.

3.1 DInSAR datasets and processing

In this work, three sets of interferograms were gathered and used to characterize vertical, fast displacement over the whole mining area in Upper Silesia (Table 1). Part of the interferograms was already processed within previous research projects realized in PGI-NRI. Newly processed interferograms supplement existing sets to fulfil the twenty years as much as possible. The total of 81 interferograms were used:

1. (1)

   18 from the archival ERS and ENVISAT C-band scenes, from the period 30 July 1993–5 July 2010. Although the total time of acquiring the scenes is very long, it should be noted that the interferograms do not cover the entire period but only selected one-two-month intervals. The years 1993–2007 in the analysis are therefore characterized by incomplete information about the deformations at that time. Nevertheless, based on the obtained images, it was possible to assess the scale and extent of subsidence at the turn of the 20th and 21st centuries. The data’s large time gaps were caused by acquiring incoherent phase images in these periods. The SLC scenes from frames 494, 222, 143, and 415 were downloaded from the archives of the European Space Agency and processed in SARScape ENVI software.

2. (2)

   26 from the processing of ALOS L-band scenes, from the period 22 February 2007 to 13 August 2009. The L-band scenes showed great relevance for monitoring ground movement in mining areas in previous studies. L-band wavelength penetrates vegetation; therefore, interferograms are highly correlated, even in non-urbanized areas. The previously processed ALOS interferograms from the Terrafirma project were used (Przyłucka et al. 2015a), supplemented with archival scenes of the same satellite, downloaded from the NASA Earth Science Data archive. Data were processed for four overlapping frames (paths 624 and 625, frames 990 and 1000) in SARScape ENVI software.

3. (3)

   37 from the X-band scenes of TerraSAR-X satellite, orbital frame 108, from 5 to 2011 to 27 January 2013. 28 interferograms with 11-days intervals were already processed by Tele-Rilevamento Europa—T.R.E within the DORIS EU project covering the period 5 July 2011- 21 June 2012. The complete description of this processing can be found in (Przyłucka et al. 2015b). Further, 9 interferograms from 21 to 2012 to 27 January 2013 were processed in SARScape software and supplemented the information with the next half of the year of subsidence. Interferograms from X-band are characterized by high spatial resolution and good phase correlation; however, it was often observed that in the centre of the subsidence basin, where the displacement is very fast, the displacement pattern is disturbed or underestimated.

3.2 PSI datasets and processing

The work also included five sets of PSI processing, covering different time and spatial ranges.

1. (1)

   Two sets, processed within the Terrafirma project (Graniczny et al. 2005, 2014a, b), come from ERS C-band satellite and the cover period 1992–2003. Dataset covering the Northern part of Upper Silesia (the “Sosnowiec” area) was acquired from 54 scenes and resulted in 30,869 PS points. Dataset covering the Southern part of Upper Silesia (the “Rybnik” area) comes from the processing of 77 scenes and consists of 64,734 PS points.

2. (2)

   One similar ERS dataset was processed within the Doris project (Del Ventisette et al. 2013) and covered the period 1992–2000 and only the central part of Upper Silesia (the “Katowice” area). The processing included 70 ERS scenes and resulted in 71,709 PS points.

3. (3)

   One set from the Envisat satellite covered 2003–2010 and was also acquired within the Doris project. The set consists of 32,341 PS points and covers only the central part of the Upper Silesia, the “Katowice” area.

4. (4)

   One set from the TerraSAR-X satellite was based on the set of X-band scenes used for DInSAR interferograms and covered 2011–2012. Thanks to high spatial resolution, the density of PS points is significantly higher than the C-band datasets, and the processing resulted in 1.7 million PS points covering most of the Upper Silesia central part area.

The summary of all InSAR datasets is presented in Table 1 and their spatial coverage on Fig. 1.

Spatial coverage of InSAR datasets over Upper Silesian Coal Basin

The datasets significantly differ from each other. There were two different processing techniques used, but also within some groups of datasets, there are differences in spatial and temporal data cover.

Figure 2 presents a graph with all 81 interferograms divided by satellite frames of SAR images used. The graph illustrates the temporal coverage of each interferogram. It can be noted that archival C-band images are characterized by around one-month intervals and large temporal gaps, which are due to the fact that most of the archival interferograms did not provide phase images coherent enough and were disregarded. Although the ERS and Envisat SAR images usually cover the whole Upper Silesia and are easy to process, they do not provide complete information about the deformations. On the other hand, 2007–2009 are fully fulfilled by L-band data with long interferogram time intervals. Years 2011 and 2012 are covered by X-band short time intervals (usually 11 days) from X-band scenes.

In order to gather on the one map information about ground subsidence from all available raster and vector InSAR data, a methodology presented in the following paragraph was applied.

4.1 DInSAR cumulative subsidence map

In order to combine information on subsidence, the differential interferograms were elaborated as follows. Firstly, by the phase unwrapping, each acquired interferometric image was converted into a deformation surface. In this way, the images of the fringes were changed into surface deformation maps, where subsidence value was assigned to each pixel. The raster image was converted to a contour image in order to remove erroneous displacement information. At this stage, manual interpretation of the obtained subsidence was of great importance. Most interferograms contained significant errors (e.g., subsidence values less than 2–5 cm were found over the entire area). The operator focused on selecting only those places where mining-related subsidence troughs characteristic for the area of Upper Silesia were visible. The subsidence contours formed for these areas were considered in a further stage of the analysis, i.e. the re-conversion of line contours into a raster image. In this way, subsidence maps for all periods of the interferometric images were created. Then, maps of cumulative deformations for Upper Silesia were developed based on summing up successive rasters. The time interval of each selected separated deformation surface of Upper Silesia is presented in Fig. 2 (“DInSAR surfaces”). These displacement surfaces were sequentially added to each other, creating successive images of cumulated subsidence and development of troughs in mining areas. Finally, a map was prepared to contain the subsidence from all interferograms used.

Temporal cover of DInSAR interferograms used. Each horizontal line refers to the frame of SAR images used for interferograms creation. The lowest line named “DInSAR surfaces” refers to the images of subsidence over the whole USCB created based on interferograms

4.2 Mining influence zones from PSI data

The PSI data are vector point datasets containing information on small vertical displacements in mm/yr units. The maximum values of subsidence recorded in the used datasets were − 16 mm/yr (PSI ERS “Rybnik”), − 23 mm/yr (PSI “Sosnowiec”), − 28 mm/yr (Envisat “Katowice”), − 40 mm/yr (PSI ERS “Katowice”) and − 335 mm/yr (PSI TSX). Apart from the X-band data, PS points did not allow the observation of larger and faster subsidence, characteristic for the formation of troughs in the mining areas. Nevertheless, all the datasets used presented an overall image of the changes taking place in the study area. The colour scale data clearly determines the area affected by the subsidence (Fig. 3).

In order to determine the limits of the residual subsidence values, the following methodology was used. In the first stage, the PS points sets were interpolated using the kriging method. Then, from the resulting surfaces, contours with the values of − 4 and − 8 mm/yr were selected. The value of −2 mm/yr is considered the limit of the reliable processing value of PSI (Crosetto et al. 2016). It was considered that for used datasets, mainly archival C-band which are characterized by a lot of noise its multiplication would determine credible border of minor subsidence accruing in the area and named zone one. The second contour (−8 mm/yr) is the boundary of designated zone two and marks the faster displacement areas at the center of which it is possible to find individual subsidence troughs. In the next stage, the zones of mining activity influence obtained from individual data sets were combined. In this way, an image of the zones of mining influence on the surface was created.

a Used PSI datasets. Yellow-orange and red points represent subsidence values; b Interpolated PSI datasets with designated contour values of − 4 and − 8 mm/yr velocity

The map of ground subsidence caused by mining in Upper Silesia in 1992–2012 based on InSAR data

Comprehensive merging of InSAR data resulted in creating a map illustrating the extent and scale of terrain subsidence caused by mining in Upper Silesia in the years 1992–2012 (Fig. 4). Single interferograms made it possible to create maps of rapid subsidence, appearing immediately after the conducted exploitation. The preparation of the cumulative image of the vertical displacements made it possible to identify those places where new subsidence troughs appeared in the analyzed period, changing the ground in a way typical for underground exploitation areas. A total of 312.5 square kilometres of such sites have been identified, of which 18 square kilometres overlap with urbanized areas.

The designated zones of mining influence, using the PSI datasets, supplemented the information on subsidence basins. They show an area where slow, residual changes in the terrain are likely to occur. The area of the designated first zone covers as much as 521 square kilometers, while the area of the second zone, where subsidence was greater than 8 mm/yr, is 152 square kilometers. The comparison with topographic data showed that there were 98 square kilometers of compact buildings under the influence of ground movements.

The final map allows assessing the impact of mining on the land surface in Upper Silesia over 20 years at the turn of the 20th and 21st centuries. The long analysis period is undoubtedly an advantage of the work, which, together with a large area of research, gives for a first time a general, large-area picture, especially since it concerns periods in the last 30 years and therefore is impossible to study using classical geodetic methods. The map shows the enormous influence of mining on the surface. In combination with topographic data, it can provide generalized information about endangered areas for local and regional decision makers.

Cumulative subsidence based on DinSAR surfaces from February 2007 to December 2012 with a six-month interval, with the ranges of mining operation carried out at that time

Unfortunately, the data provided are incomplete. The information about subsidence from the interferograms was burdened with errors due to manual data cleaning, incoherent phase image in many places, especially those covered with vegetation, and above all, long periods of lack of data. Therefore, the resulting subsidence values should be considered inaccurate, somewhat strongly underestimated, and incomplete. The map can be used for a general analysis of mining influences. However, since the displacement values identified by InSAR are underestimated (due to the lack of information in some months or years), the obtained values should be read as “not less than”. Therefore, the red spots on the final map are areas that have settled not less than 1 m.

There are also places of subsidence identified on interferograms but not located in the zones of influence defined by the PS points, especially in the vicinity of the towns of Żory, Tychy, Bieruń, Chelmek, Libiąż, and Brzeszcze. This is due to the lack of PSI data for these areas.

In addition to the general picture of the long-term impact of coal mining on the surface, the presented elaboration of the data can be used for a detailed analysis of the development of subsidence basins in selected places. The individual surfaces of ground changes (Fig. 2, “DInSAR surfaces”) allow for the development of the final summary image of subsidence, but also intermediate images showing the increase in displacements over time. An example of such an analysis is presented in Fig. 5 for the city of Bytom. Subsequent images represent subsidence with a six-month interval from February 2007 to December 2012. The juxtaposition of images allows for tracking the increase in changes in the land surface, and thus the development of the basin threatening the city in relation to the conducted exploitation.

The article elaborates on InSAR data coming from different processing techniques, covering different time and spatial ranges. The result of the study is the map showing the impact of underground coal mining on vertical changes in the area that occurred in 1992–2012 in the Upper Silesian Coal Basin. The primary purpose of the work was to combine various data sets for comprehensive long-term and large-scale analysis. The resulting final map contains two layers: a subsidence values layer obtained based on traditional DInSAR differential interferograms, of which a total of 81 were used, and a layer of boundaries of the occurrence of residual ground displacement divided into two zones.

The final map shows the tremendous influence of mining on subsidence. In the analyzed 20-year period, on 5600 square kilometers in the USCB area, over 500 square kilometers of more or less subsidence were detected, of which almost 100 square kilometers were located in urbanized areas. Over 300 square kilometers have been identified as areas of rapid displacement seen on interferometric images. In addition, the individual images of terrain deformation allow for the analysis of the development of individual subsidence troughs. This may be of particular importance for developing a forecast of the impact of exploitation and detailed spatial urban planning.

Satellite Radar Interferometry has once again proved its rightness to use it for monitoring mining-induced ground displacement. The presented work focuses on its unique ability to analyze large areas and to study changes that occurred in the past, in this case, three decades ago. No other monitoring technique provides the ability to perform such detailed analyzes so far.