1 Introduction

Over the years, the geodetic capabilities for determining land surface changes have developed significantly and also have found application in areas particularly affected by mining activities. Nowadays, leveling, gravimetry, photogrammetry, Light Detection and Ranging (LiDAR), Interferometric Synthetic Aperture Radar (InSAR), and Global Navigation Satellite Systems (GNSS) techniques are the most widely adopted methods for ground deformation monitoring. However, providing a high-accuracy spatiotemporal system for 3-D terrain motion measurements can be challenging using only one of the mentioned methods. To provide a continuous ground monitoring service for ongoing long-term movements, the GNSS and InSAR methods are commonly implemented (Shen et al. 2019; Del Soldato et al. 2021).

We performed a thorough review of the strengths and weaknesses of both techniques, which is available in Appendix A (Table 5). In summary, the InSAR is a few-day periodic, one-dimensional, high spatial resolution, remote sensing technique (no ground-based equipment required), whereas the GNSS is a continuous, three-dimensional, station-dependent resolution measurement technique (requires ground-based equipment). The InSAR one-dimensional Line of Sight (LOS) observations taken using different geometries could not be directly decomposed into all three coordinate components (North, East, Up) without introducing bias (Brouwer 2021). Based on the presented analysis, we demonstrate that the combination of GNSS and InSAR as two complementary methods allows to establish an overall system for the determination of three-dimensional displacement values and movement rates. In this study, we focus on exploiting strengths and reducing weaknesses of GNSS and InSAR techniques by providing a new methodology of integration involving Kalman filter algorithms for nonlinear strong ground displacements that prevent the use of multi-temporal InSAR (MT InSAR) techniques (e.g., Przyłucka et al. 2015).

Due to the near-polar orbit of the satellites, the limited sensitivity of the northern component is one of the main problems of conducting three-dimensional InSAR monitoring. Hu et al. (2014) provide a systematic review of the methods for mapping 3-D displacements using InSAR measurements pointing out the strengths and weaknesses of each approach. The presented by Hu et al. (2014) classification has three groups: (I) a combination of multi-pass LOS and azimuth measurements from at least two different geometries (e.g., Ng et al. 2011; Yang et al. 2017), (II) integration of GNSS and InSAR data (e.g., Bozso et al. 2021; Liu et al. 2019), and (III) prior information assisted approaches (e.g., Mohr et al. 1998; Kumar et al. 2011). The authors highlight that the methods from the first group are suitable for investigating the large displacements caused mostly by natural phenomena (e.g., earthquakes, volcanic eruptions). The methods of GNSS and InSAR integration (group II) are sufficient for both large-scale and local displacements; however, the final outcome is highly dependent on the number and distribution of the GNSS receivers. The third group, related to prior data, can be applied only in combination with external data sources, e.g., glacier or landslide movement models.

In the presented study, we refer directly to the (II) group of the mentioned methods—the integration of GNSS and InSAR techniques. Currently, the fusion process involves different approaches and assumptions, such as linear prediction for large-scale displacements or spatiotemporal interpolation for rapid nonlinear movements. For instance, Fuhrmann et al. (2015) presented an integrated method using leveling, GNSS (76 stations) and InSAR (ERS-1/2 and Envisat data from ascending and descending tracks) to determine the linear velocities over a larger area (around 32 000 km\(^2\)). The nonlinear displacements (located mostly in regions with man-induced surface deformations) were excluded from the investigation. The first step of the proposed method relied on determining the offset and trend between the InSAR data and the other two techniques. In the second step, local tangent rates were estimated using the least-squares adjustment (LSA) methodology. The precision of the obtained results was estimated as 0.30 mm/a in the horizontal and 0.13 mm/a in the vertical direction.

Another concept applies linear interpolation and prediction based on the specific spatial and temporal domains of the GNSS and InSAR data (Liu et al. 2019). The authors used a dynamic filtering method that relied on a spatiotemporal Kalman filter to predict the large-scale deformations in the area of Los Angeles, USA. The research presented in this study was based on simulated and real data. The real observations comprised 15 InSAR images from the Envisat descending orbit and data from 35 GNSS receivers. Due to the single geometry orbit, horizontal and vertical deformations were not resolved, and GNSS coordinates were projected into the LOS direction. The quality of the deformation field model was established using the leave-one InSAR image-out validation method (one InSAR image was taken as the reference and compared with the corresponding image generated from the spatiotemporal deformation model). The root mean square (RMS) between the interpolated and the reference InSAR image did not exceed 7.2 mm in the LOS direction, whereas the maximum deformation did not exceed 100 mm.

More recently, a new type of procedure to integrate Differential InSAR (DInSAR) and GNSS results was presented in the study of Bozso et al. (2021). The paper shows a method for the estimation of 3-D deformation parameters and rates using Sentinel-1 SAR data and GNSS epoch observations converted to the Sentinel-1 LOS domain. The Kalman filter algorithm was tested on a landslide area in Hungary with the assumption of linear velocities in the LOS direction. A network of geodynamic benchmarks, consisting of a pair of corner reflectors (aligned with the ascending and descending Sentinel-1 geometry) and a GNSS antenna adapter for campaign measurements, was created in the study area. To determine the 3-D deformation parameters, the GNSS observations were projected onto the radar LOS vectors. The presented approach was based only on linear velocities calculated from two GNSS campaigns. As a model of the system noise, the a priori values were assumed as 2 mm and the a posteriori errors of the NEU components were 16, 2, and 3 mm, respectively. However, this way of uncertainty assignment is not rigorous.

The Kalman filter is already widely used in the analysis of InSAR (Shirzaei and Walter 2010; Dalaison and Jolivet 2020) and GNSS (Pacione et al. 2011; Bian et al. 2014) time series; however, only a few studies (Bekaert et al. 2016; Liu et al. 2019; Bozso et al. 2021) use it for the integration of these two types of measurement. The strength of properly parameterized Kalman filtering is the ability to estimate the state of a system with high accuracy, even when the measurements are noisy, incomplete, and determined in different spatiotemporal domains. Moreover, we can conclude, based on the methods presented in the literature, that nowadays the GNSS-InSAR integration is mostly applied in areas with high-spatial resolution GNSS receivers and linear motions (Liu et al. 2019). However, for local nonlinear strong deformations, current integration techniques are not well developed (Ng et al. 2011; Bozso et al. 2021) and are not robust against different levels of observation noise. A sufficient nonlinear ground movements monitoring system can be important for early warning issues, understanding the deformation process, assessment of risk and damage, and evaluation of mitigation measures. Moreover, due to the ambiguous nature of observations, MT InSAR techniques such as Persistent Scatterer Interferometry (PSI) or Small Baseline Subset Interferometry (SBAS) have a limited capability to measure rapid nonlinear deformation phenomena, greater than several tens of centimeters per month (Przyłucka et al. 2015; Pawluszek-Filipiak and Borkowski 2020).

Given the literature review, presented above and the need to find a way to monitor local nonlinear strong deformations, we come up with the following research questions:

  • What functional and stochastic model will be required to integrate on the observation level GNSS and InSAR data to retrieve consistent 3-D deformations in the underground mining region?

  • How to overcome limitations of DInSAR measurements, unwrapping problems, varying, and unstable coherence, given non-existing MT InSAR in the area of interest?

  • How to utilize the GNSS point measurements as anchor points to constrain the DInSAR field observations with decomposition into 3-D deformation?

The study presented in this manuscript extends the Kalman filter algorithm described by Bozso et al. (2021), using GNSS observations from permanent ground stations and DInSAR deformation time series derived from Sentinel-1 satellites. Our algorithm is able to ingest the time series of GNSS topocentric coordinates with gaps of several months, and noisy time series of DInSAR ascending and descending LOS velocities subject to troposphere artifacts or SAR phase unwrapping errors. The latter errors appear in different situations but mostly either in the case of increased phase noise due to decorrelation effects (e.g., dense vegetation, snow), severe weather conditions, strong and abrupt displacement behavior, or a combination thereof. The observation uncertainties are rigorously determined from the parameter estimation process for the GNSS data and from coherence values for DInSAR velocities (see maps of mean cross-correlation coefficient in Appendix B (Fig. 7)). To verify the fusion results, a quality analysis based on independent measurements was performed. The verification procedure was implemented using the external survey data source - GNSS campaign points named also as epoch-based measurements. The proposed methodology is tested in an area of intensive underground works within the Upper Silesian Coal Basin (USCB) in Southern Poland—one of the biggest deposits in Europe. The surface deformations in the area have been monitored in recent years by different geodetic approaches, revealing a complex regime of multiple dynamic subsidence patches with a substantial range of displacement values reaching in some cases 1.5 m per year and nonlinear character of the deformations, corresponding to the distributed underground mine works along the longwalls—among the others, Del Ventisette et al. (2013), Przyłucka et al. (2015), Ilieva et al. (2019), Pawluszek-Filipiak and Borkowski (2020), Kopeć (2021). In Sect. 2, we introduce our methodology for the integration of DInSAR and GNSS measurements. The results of our algorithm are reported in Sect. 3. In Sects. 4 and 5, we present the discussion and conclusions of our work.

2 Methodology

The paper presents an original methodology for the integration of two different techniques, optimal for nonlinear motions, conducted for an area affected by underground mining works. In the approach, we processed geodetic data types namely SAR and GNSS observations, aiming to achieve comprehensive 3-D monitoring of the ongoing ground deformations. The estimation workflows are presented in Sects. 2.1 and 2.2. Further, to compare the DInSAR and GNSS time series, a procedure unification of both data sets is introduced (Sect. 2.3). Section 2.4 shows a new method for integration of DInSAR and GNSS data using a Kalman filter. Finally, to verify the obtained results, the quality analysis based on independent data set is performed (Sect. 2.5).

2.1 DInSAR processing

In the current study, a consecutive cumulative DInSAR approach for surface deformation monitoring is applied (Fig. 1a), where the second (secondary) image of each interferogram is used as a first (primary) image in the next interferometric pair, forming continues time series (Fig. 1b). The consecutive DInSAR method provides a simple and fast tool for monitoring by adding new data to the subsidence stack every time a new SAR image is released. It is also relatively easy to automate the process. A serious disadvantage is that this approach does not eliminate the contribution of the interferograms with lower quality, e.g., due to temporal or geometric decorrelation (Hanssen 2014). On the other hand, in the case of atmospheric artifacts registered in one SAR image, caused by a higher percentage of water vapor which slows down the radar signal, they will be canceled out by adding the affected image with a negative sign in the following pair. However, the atmosphere delay artifacts from the first and last image will persist in the results. The DInSAR processing in this study is accomplished with the open science Sentinel Application Platform (SNAP), provided by the European Space Agency (ESA), following the standard DInSAR procedure, outlined below.

Fig. 1
figure 1

Principle of the consecutive cumulative DInSAR technique for surface deformation monitoring applied in the current study. a An example of real data from the Sentinel-1 A and B satellites’ positions (perpendicular Baselines) within a month time span showing the successive formation of DInSAR interferograms using common images between the pairs. b Relative displacements (gray solid lines) for the time of acquisition of the primary and secondary images forming each interferogram, and the cumulative displacement along the LOS (black solid line)—an example from Sentinel-1 A and B ascending data, that correspond to the time scheme above in (a)

Full size image

Sentinel-1 data covering the area of interest acquired in Interferometric Wide swath (IW) mode with VV polarization for the time span of the GNSS measurements (Sect. 2.2) are used. Sentinel-1 is a constellation of two satellites (A and B) maintained by ESA that have mounted C-band (wavelength of \(\sim 5.55\) cm) SAR sensors and provide a 6-day revisit time within the study period. The generated SAR images are freely available on the ESA’s Sentinel hub (https://scihub.copernicus.eu). The precise Sentinel-1 orbit auxiliary files prepared by Copernicus Precise Orbit Determination (POD) Service are downloaded automatically from Copernicus Sentinels POD Data Hub (https://scihub.copernicus.eu/gnss) and are applied to the SAR data to minimize the satellite orbital error. An external Digital Elevation Model (DEM), namely the 1-sec distribution of the Shuttle Radar Topography Mission (SRTM) Height file (Jarvis et al. 2008) (http://srtm.csi.cgiar.org), is used for coregistration of the secondary image to the geometry of the primary one. During this procedure, a 21-point truncated sinc interpolation is applied. Then, the secondary image is corrected for range and azimuth shifts by implementing the Network Enhanced Spectral Diversity (NESD) approach (Fattahi et al. 2017). The step is needed due to the specific format of the Sentinel-1 data which are acquired with the Terrain Observation with Progressive Scans SAR (TOPSAR) technique by switching the SAR antenna between neighboring sub-swaths on a number of bursts. The NESD algorithm aims to correct the shifts in the overlap areas of adjacent bursts. The interferogram is computed by subtracting the flat-earth phase, while the SRTM is used again to subtract the stable topography component resulting in the wrapped interferogram. A Goldstein non-adaptive filter (Goldstein and Werner 1998) was applied to reduce speckle noise, to ease the two-dimensional spatial phase unwrapping, which solves for the 2\(\pi \) ambiguity. The Statistical-Cost, Network-Flow Algorithm for Phase Unwrapping (SNAPHU) (Chen and Zebker 2001) integrated in SNAP was applied, and the values of the resulted unwrapped interferometric phase (\(\Delta \phi _d\)) are transformed to metric units based on the relation:

$$\begin{aligned} d_{LOS} = -\frac{\lambda }{4\pi } \Delta \phi _d, \end{aligned}$$
(1)

where \(d_{LOS}\) is the relative displacement projected on the radar LOS, and \(\lambda \) is the wavelength of the transmitted radar signal mentioned above.

The quality of the processed interferometric phase is represented by the cross-correlation coefficient (\(\gamma \)) between the two SAR images of each interferogram expressed as a level of coherence at each pixel (Hanssen 2014):

$$\begin{aligned} \gamma =\frac{\left| \sum _{n=1}^{N}y_{1}^{\left( n \right) }y_{2}^{*\left( n \right) } \right| }{\sqrt{\sum _{n=1}^{N}\left| y_{1}^{\left( n\right) }\right| ^{2}\sum _{n=1}^{N}\left| y_{2}^{\left( n\right) }\right| ^{2}}}, \end{aligned}$$
(2)

where \(y_1\) and \(y_2\) are the complex values at the observed pixel of the primary (1) and secondary (2) SAR images, and N is the number of pixels in a surrounding estimation window.

The total coherence (\(\gamma _{tot}\)) depends on the influence of several decorrelation factors (Hanssen 2014):

$$\begin{aligned} \gamma _{tot}=\gamma _{geom}.\gamma _{DC}.\gamma _{vol}.\gamma _{thermal}.\gamma _{temporal}.\gamma _{processing}, \end{aligned}$$
(3)

where \(\gamma _{geom}\) presents the decorrelation due to geometry differences between the two SAR images, \(\gamma _{DC}\)—differences in the Doppler centroids between the two images, \(\gamma _{vol}\)—the effect of the scattering medium on the radar signal, \(\gamma _{thermal}\)—noise due to the characteristics of the system, \(\gamma _{temporal}\)—decorrelation due to the physical changes of the scattering characteristics of the ground over time, and \(\gamma _{processing}\)—influence of the used processing algorithms. Some of these effects are limited either by providing better orbital tracks or limiting the time span between the SAR acquisitions, while others are highly dependent on the system characteristics, e.g., used wavelength. The lower coherence in bigger areas of the interferogram impedes the proper solving of the phase ambiguities resulting in jumps in the time series (Sect. 3.1).

To determine the error values of a particular pixel, the coherence coefficient (Eq. 3) was used to estimate the phase noise (\(\sigma _{\Delta \phi }\)). The phase variance can be established using the probability density functions (PDF) depending on multilook levels (Hanssen 2014). The phase variance resulting from \(\gamma _{tot} < 1\) can be defined as:

$$\begin{aligned} \sigma _{\Delta \phi }^2=\int _{-\pi }^{+\pi }[\phi -E\{\phi \}]PDF(\phi )d\phi , \end{aligned}$$
(4)

In the presented study, the multilooking level is equal to 1 (\(L=1\)), and thus, the phase variance can be expressed as:

$$\begin{aligned} \sigma _{\Delta \phi ,L=1}^2= & {} \frac{\pi ^2}{3}-\pi \arcsin (|\gamma _{tot}|)+\arcsin ^2(|\gamma _{tot}|)\nonumber \\{} & {} -\frac{1}{2}\sum _{k=1}^{\infty }\frac{|\gamma _{tot}|^{2k}}{k^2}, \end{aligned}$$
(5)

where k is the index of summation. The final interferometric phase errors in the ascending and descending directions (\(\sigma _{d_{LOS}^A}, \sigma _{d_{LOS}^D}\)) are transformed to metric units based on the relation:

$$\begin{aligned} \sigma _{d_{LOS}}=\frac{\lambda }{4\pi }\sigma _{\Delta \phi ,L=1}. \end{aligned}$$
(6)

Since the unwrapped displacements are relative, all maps of displacement must be related to a common reference. Instead of one reference point, a set of points is used by applying the method proposed in Ilieva et al. (2019). The set of points is chosen based on stable coherence characteristic over time for both orbits—ascending and descending. These points are also away from local nonlinear strong deformation zones, and the time span for points selection is the same as the study time (here 2018–2020). A correcting surface for each map is interpolated between the pixels that preserved the highest coherence (more than 90%) over a year.

As it was mentioned before, the DInSAR displacements are calculated in LOS. In order to be able to compare the 3-D GNSS data with the DInSAR results, the heading and incidence angles are also extracted from the SAR data. The values of the incidence viewing angle (\(\theta _{inc}\)) are exported for each interferogram. The satellite’s heading angle can be treated as a constant value for each geometry, since the variation is usually within 1 deg over the extent of a SAR image and hence has a negligible influence (Fuhrmann and Garthwaite 2019). Thereby, in the current study, a mean value of the heading angles (\(\alpha \)) for the two SAR images participating in an interferogram formation is taken for further calculation. Thus, the deformation vector measured in the LOS direction can be described by three-dimensional components:

$$\begin{aligned} d_{LOS}= & {} \begin{bmatrix} \sin (\theta _{inc})\sin (\alpha ) -\sin (\theta _{inc})\cos (\alpha )&\,&\cos (\theta _{inc}) \end{bmatrix}\nonumber \\{} & {} \begin{bmatrix} n_S \\ e_S \\ u_S \end{bmatrix}, \end{aligned}$$
(7)

where \(d_{LOS}\), \(n_S\), \(e_S\), and \(u_S\) are the LOS, North, East, and Up deformations, respectively, and index S is the indication of the SAR technique.

2.2 GNSS processing

The study is focused on the ground changes above the Rydułtowy coal mine, located in the south-western part of USCB. Rydułtowy is the oldest continuously operational mine in Upper Silesia (since 1792) with a planned lifetime of exploitation up to 2040 (https://pgg.pl). The mining area covers 46 km\(^2\) with operational resources of around 80.4 million tons and a production capacity of 9000–9500 tons per day. The longwall method of coal extraction is used in this mine with the depth of the seams between 800 and 1200 m. The work in the region of interest E1 of the mine is ongoing on several walls at 4 seams, each with a thickness of about 18 m (Ćwiȩkała 2019). The investigated area covers two walls, with a spatial extent of 2 km\(^2\), and a moving, nonlinear deformation bowl with a diameter of 200 m.

The coordinates of the permanent GNSS stations were estimated in post-processing mode using a double difference technique. The solution was based on the Global Positioning System (GPS), the Russian Global Navigation Satellite System (GLONASS), and European Galileo observations. The estimation of the coordinates was performed with a daily computational interval using Bernese GNSS software v.5.2 (Dach et al. 2015). The processing scenario was based on the final precise GNSS products. The mentioned products included the orbits and clocks of the multi-GNSS satellites obtained through the Center for Orbit Determination in Europe (CODE, Dach et al. (2016)). To ensure a stable reference, 14 stations belonging to the International GNSS Service (IGS, https://igs.org) network and the European Reference (EUREF) Permanent Network (EPN, https://www.epncb.oma.be) were included in post-processing calculations.

Fig. 2
figure 2

Map of the study area including undergoing mining exploitation panels at seam 703/1 (gray rectangles), locations of permanent GNSS stations (green circles and orange triangles), and campaign survey points (purple pentagons). The map in the bottom right corner shows the locations of IGS/EPN reference stations (red squares)

Full size image

In the research area affected by mining deformations (Fig. 2), two types of permanent GNSS stations can be distinguished: geodetic and low-cost receivers deployed under the EPOS-PL project (Mutke et al. 2019, https://epos-pl.eu). One geodetic receiver was located on the building “Ignacy” Historical Mine (RES1). The site was not placed directly above the longwall panels; however, ground deformations were expected at this location as well. The second type of permanent station was equipped with five low-cost receivers (PI02, PI03, PI04, PI05, PI16). These stations were installed above the longwall of “Rydułtowy” mine currently exploited. The beginning of collecting observations was 11.02.2019; however, due to technical problems, the low-cost time series are not complete (with a maximum break of around 9 months for the PI02 station). In general, the low-cost station can be defined as an inexpensive, light, small, and low-power consumption device which is able to obtain position accuracy at sub-centimeter level in static positioning. The overall price of a low-cost GNSS instrument was variously defined as less than 500€ (Cina and Piras 2015) or below 1 000$ (Jackson et al. 2018). The most important aspect of using low-cost stations in ground deformations monitoring is to maintain the quality of position determination. Various studies conducted for different types of low-cost receivers indicated sub-millimeter or millimeter differences compared to traditional geodetic GNSS stations (Cina and Piras 2015; Hamza et al. 2020, 2021).

In addition to the mentioned permanent stations, independent epoch-based measurements were also considered in this study. The GNSS campaign network contains over 120 well-stabilized and repetitively measured points developed by the Military University of Technology in Warsaw (MUT) across the Upper Silesia region. In the area of interest, 42 campaign points were located over underground “Rydułtowy” coal seams (Fig. 2). The MUT team processed the epochs of observations together with the GNSS data from the local network of the Continuously Operating Reference Station (CORS) collected in the GNSS Data Research Infrastructure Centre (Mutke et al. 2019, http://www.gnss.wat.edu.pl). The measurements at selected epoch-based points were completed within one-to-two surveying days. Data were collected at 1 or 15 s intervals with a time span of 3.2 to 5.5 h. The processing strategy followed the EPN guidelines. It is also consistent with the MUT analysis for operational EPN products and was used successfully in the EPN Repro2 project (Araszkiewicz and Völksen 2017). To provide consistency with the ITRS/ETRS realization, the estimated results were combined with the MUT-PL solution created within the EPOS-PL project (Mutke et al. 2019). The details about GNSS equipment and the execution of the measurements are presented in Table 1.

Table 1 Classification of GNSS infrastructure deployed in Upper Silesia under the EPOS-PL project (names of the stations are consistent with Fig. 2)
Full size table

The results obtained for the permanent and campaign GNSS points were determined in the ITRF2014 reference frame Altamimi et al. (2016). In order to ensure consistency with the InSAR data, which assume unchanging positions with respect to the reference points, the ITRF2014 coordinates and uncertainties were transformed to the ETRF2000 reference frame using the 7-parameter approach (Altamimi 2018). Afterward, for a better understanding of the local displacements, the geocentric ETRF2000 coordinates were converted to the topocentric reference frame (Tao et al. 2019). An illustrative way to present the position variation of a GNSS station over time is to take the first computational epoch as the reference value. To avoid adopting an outlier as an initial value in the time series of the permanent stations, the averaged coordinates (\(N_G, E_G, U_G\)) from the first five computational observations were assumed as reference values. The reference values for the positions of the campaigned GNSS points (\(N_C, E_C, U_C\)) were determined as the first epoch that is common in time with the data from the permanent stations. To establish the uncertainties for the permanent and epoch-based topocentric coordinates, the errors conversion was defined as:

$$\begin{aligned} \begin{bmatrix} \sigma _{N}^2 &{} \sigma _{NE} &{} \sigma _{NU} \\ \sigma _{NE} &{} \sigma _{E}^2 &{} \sigma _{EU}\\ \sigma _{NU} &{} \sigma _{EU} &{} \sigma _{U}^2 \end{bmatrix}= F \begin{bmatrix} \sigma _X^2 &{} 0 &{} 0 \\ 0 &{} \sigma _Y^2 &{} 0\\ 0 &{} 0 &{} \sigma _Z^2 \end{bmatrix} F^t, \end{aligned}$$
(8)

where \(\sigma _{N}, \sigma _{E}, \sigma _{U}\) are the transformed topocentric GNSS variances, \(\sigma _{NE}, \sigma _{NU}, \sigma _{EU}\) are the covariances, and \(\sigma _X, \sigma _Y, \sigma _Z\) are the errors of the geocentric ETRF2000 coordinates. F is the rotation matrix determined as:

$$\begin{aligned} F = \begin{bmatrix} -\sin (B)cos(L) &{} -\sin (B)sin(L) &{} \cos (B) \\ -\sin (L) &{} \cos (L) &{} 0 \\ \cos (B)cos(L) &{} \cos (B)sin(L) &{} \sin (B) \end{bmatrix}, \end{aligned}$$
(9)

where B, L are the longitude and latitude in ETRF2000 of the reference point.

2.3 Unification of the DInSAR and GNSS deformations dimensions

The DInSAR measurements enable spatial continuous monitoring of large-scale subsidence areas. However, the radar technique acquires the observations in ascending and descending 1-D LOS directions. Whereas, the GNSS method provides determination of the displacements in 3-D, but only on sites where the antennas were located. Therefore, to perform a unification of both techniques, it was necessary to overlay the DInSAR rasters with the GNSS locations and to focus our study on the intersecting pixels.

Fig. 3
figure 3

Scheme demonstrating the proposed method for unification of the DInSAR and GNSS observations (top part of the graph, described in Sects. 2.1 and 2.2) in absolute temporal domain, realized through DInSAR decomposition (bottom left part—Sect. 2.3), including quality analyses of the results (bottom right part) presented in Sect. 2.5

Full size image

As it was mentioned in Sect. 2.1 (Eq. 7), theoretically the DInSAR LOS could be decomposed to its 3-D components \(n_S\), \(e_S\), and \(u_S\) if a sufficient number of observations are available. In the common case two DInSAR data sets, namely from ascending (Asc or A) and from descending (Desc or D) satellite geometry, can be used for the decomposition of LOS vector in horizontal and vertical projections (Fig. 3). Due to the nearly north–south trajectory of the SAR satellites, the system has limited sensitivity to the ground movements in this direction. That is why often the parameter \(n_S\) is neglected so the large noise driven by the north component of the LOS displacement is minimized in the decomposed in 3-D LOS vector. In this case, the measurements from the two DInSAR geometries (\(d_{LOS}^A\) and \(d_{LOS}^D\)) could be enough to define the other two unknowns—\(e_S\) and \(u_S\), though the pattern of the horizontal movement will be incomplete and especially the up component will be biased (Brouwer 2021). To minimize this problem, in this study the GNSS North displacement (\(N_G\)) was implemented instead (Samsonov and Tiampo 2006; Jun et al. 2012). Within this procedure, it was kept in mind that the two sets of observations have different temporal sampling. The GNSS technique works in an absolute mode with respect to the initial time that defines the start of the GNSS time series (Sect. 2.2). On the other hand, the DInSAR method provides only the relative deformation referring to the initial time the first (primary) image of the pair forming each interferogram. To unify the data in a common temporal domain, the GNSS \(N_G\) component was divided into intervals corresponding to the time span of the DInSAR interferometric pairs (Fig. 1). The origin for each interval was established at the beginning of the period converting the absolute temporal GNSS \(N_G\) measurement to a relative one. Then, it was introduced as the known parameter \(n_G\) in the DInSAR LOS decomposition of the DInSAR \(e_S\) and \(u_S\) estimation:

$$\begin{aligned} \begin{bmatrix} e_S\\ u_S\\ \end{bmatrix}= & {} \begin{bmatrix} -\sin (\theta ^A)\cos (\alpha ^A) &{} &{} \cos (\theta ^A)\\ -\sin (\theta ^D)\cos (\alpha ^D) &{} &{} \cos (\theta ^D) \end{bmatrix}^{-1}\nonumber \\{} & {} \left( \begin{bmatrix} d_{LOS}^A\\ d_{LOS}^D \end{bmatrix} - \begin{bmatrix} n_G \sin (\theta ^A)\sin (\alpha ^A)\\ n_G \sin (\theta ^D)\sin (\alpha ^D)\\ \end{bmatrix}\right) , \end{aligned}$$
(10)

where \(\theta ^A\) and \(\theta ^D\) are the incidence angles, and \(\alpha ^A\) and \(\alpha ^D\) are the heading angles of the ascending and descending DInSAR observations, respectively. As a final step, the relative DInSAR parameters (\(e_S,u_S\)) are converted to the cumulative, absolute temporal domain using equations:

$$\begin{aligned} E_S(t)= & {} \sum _{i=1}^te_s(i), \end{aligned}$$
(11)
$$\begin{aligned} U_S(t)= & {} \sum _{i=1}^tu_s(i), \end{aligned}$$
(12)

where t is the epoch of the DInSAR measurement. As a consequence of the relative character of the DInSAR observations, every potential outlier propagates into the subsequent epochs in the absolute time series. The maps of DInSAR cumulative East and Up deformations over the subsidence bowl located in the study area are available in Appendix C (Fig. 8). To test the impact of outliers for the DInSAR displacements (\(E_S, U_S\)), a statistical analysis of the accuracy was performed. For the validation, we used an external data source from the campaigned GNSS measurements (\(N_C, E_C, U_C\)) (Sect. 2.5).

2.4 DInSAR and GNSS data integration

The process of integrating DInSAR and GNSS observations is based on the Kalman filter (Kalman 1960) algorithm (Fig. 4). To estimate the unknown observation parameters, namely NEU displacements and rates, and the uncertainties of the system, the process is built as a recursive two-step procedure. During the first stage, termed as the time-update step, the information about the system from the previous epoch (\(t-1\)) is used to predict the parameters in the current epoch (t). In the second step, named as the measurement-update step, the modeled values, as well as the uncertainties, are corrected by the real observations and measurement noise registered in the epoch t.

Fig. 4
figure 4

Scheme of the approach for the two-step Kalman filter algorithm (bottom left part of the graph) applied for integration of DInSAR and GNSS observations (upper row of the graph), including verification of the quality based on the external GNSS campaign measurements (lower right part)

Full size image

2.4.1 Time-update segment of the data integration

The time-update part of DInSAR and GNSS integration contains the dynamic model for time-varying parameters:

$$\begin{aligned} {\widehat{x}}_{t|t-1}= & {} \Phi _{t|t-1}{\widehat{x}}_{t-1|t-1}, \end{aligned}$$
(13)
$$\begin{aligned} P_{t|t-1}= & {} \Phi _{t|t-1}P_{t-1|t-1}\Phi _{t|t-1}^T + S_t, \end{aligned}$$
(14)

where \({\widehat{x}}_{t|t-1}\) and \(P_{t|t-1}\) are the predicted state vector and error variance matrix, respectively, \({\widehat{x}}_{t-1|t-1}\) and \(P_{t-1|t-1}\) are the state vector and error variance matrix in the previous epoch, while \(S_t\) is the system noise matrix, and \(\Phi _{t|t-1}\) is the transition matrix for the current moment. The 3-D description of the ground position, common to both techniques, is considered in the following calculations by the state vector \({\widehat{x}}\), while the velocities of the ground deformations are sampled to daily intervals (\(\Delta t =\) 1 day). The processing computations are conducted in the topocentric reference frame, and the initial values of the NEU positions and rates in the state vector (\({\widehat{x}}_{0|0}\)) are equal to zero. As initial values for the errors of the parameters (\(P_{0|0}\)), the system noise values (\(S_t\)) are taken.

In our approach, the zero-mean acceleration model is introduced as the system noise matrix Teunissen (2009):

$$\begin{aligned} S_t=\sigma _0^2 \begin{bmatrix} \frac{1}{4} \Delta t^4 &{} \frac{1}{2} \Delta t^3 &{} 0 &{} 0 &{} 0 &{} 0 \\ \frac{1}{2} \Delta t^3 &{} \Delta t^2 &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{1}{4} \Delta t^4 &{} \frac{1}{2} \Delta t^3 &{} 0 &{} 0 \\ 0 &{} 0 &{} \frac{1}{2} \Delta t^3 &{} \Delta t^2 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \frac{1}{4} \Delta t^4 &{} \frac{1}{2} \Delta t^3 \\ 0 &{} 0 &{} 0 &{} 0 &{} \frac{1}{2} \Delta t^3 &{} \Delta t^2\\ \end{bmatrix}. \end{aligned}$$
(15)

The \(\sigma _0\) parameter, representing the initial standard deviation, was defined empirically, and the detailed results are presented in Sect. 3.2. The transition matrix (\(\Phi _{t|t-1}\)) of the dynamic model describes the evolution of the state parameters in time and can be obtained as a first-order differential equation:

$$\begin{aligned} \Phi _{t|t-1}= \begin{bmatrix} 1 &{} \Delta t &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 1 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 1 &{} \Delta t &{} 0 &{} 0 &{} \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} \Delta t\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} 1 \end{bmatrix}. \end{aligned}$$
(16)

2.4.2 Measurement-update segment of the data integration

In the measurement-update step, to obtain a filtered system with error variances, the prediction is combined with the actual observations in terms of minimum mean square error:

$$\begin{aligned} K_t= & {} P_{t|t-1} A_t^T (A_t P_{t|t-1} A_t^T + R_t)^{-1}, \end{aligned}$$
(17)
$$\begin{aligned} {\widehat{x}}_{t|t}= & {} {\widehat{x}}_{t|t-1}+K_t (z_t - A_t {\widehat{x}}_{t|t-1}), \end{aligned}$$
(18)
$$\begin{aligned} P_{t|t}= & {} (I-K_t A_t)P_{t|t-1}, \end{aligned}$$
(19)

where \(K_t\) is the Kalman gain matrix, \({\widehat{x}}_{t|t}\) is the system of the filtered state vector (\(N_F\), \(vN_F\), \(E_F\), \(vE_F\), \(U_F\), \(vU_F\)), and \(P_{t|t}\) is the variance measurement-updated matrix and I the identity matrix. The observations are introduced with the \(z_t\) vector (Eq. 20), while \(A_t\) is the projection matrix (Eq. 21) and \(R_t\) is the measurement noise matrix (Eq. 22):

$$\begin{aligned} z_t= \begin{bmatrix} N_G&E_G&U_G&\frac{d_{LOS}^A}{\Delta T}&\frac{d_{LOS}^D}{\Delta T} \end{bmatrix}^T. \end{aligned}$$
(20)

In the vector of the observations (Eq. 20), three elements (\(N_G\), \(E_G\), and \(U_G\)) refer to the topocentric position of the GNSS station received by the processing with Bernese software (Sect. 2.2). The two last elements (\(d_{LOS}^A and d_{LOS}^D\)) are the differential LOS displacements obtained from the DInSAR ascending and descending geometries, respectively (Sect. 2.1).

The DInSAR observations were divided by the time factor \(\Delta T\), which corresponds to the period of the interferogram time span acquisition. This interval \(\Delta T\) is equal to 6 days, or 12 days in the case of a missing Sentinel-1 image. In the projection matrix (\(A_t\)), the number of rows is equal to the number of measurements, while the number of columns determines the number of parameters:

$$\begin{aligned} A_t = \begin{bmatrix} 1 &{} 0 &{} 0 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 1 &{} 0 &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 1 &{} 0\\ 0 &{} \sin (\theta ^A)\sin (\alpha ^A) &{} 0 &{} -\sin (\theta ^A)\cos (\alpha ^A) &{} 0 &{} \cos (\theta ^A)\\ 0 &{} \sin (\theta ^D)\sin (\alpha ^D) &{} 0 &{} -\sin (\theta ^D)\cos (\alpha ^D) &{} 0 &{} \cos (\theta ^D) \end{bmatrix}.\nonumber \\ \end{aligned}$$
(21)

The two last rows are related to the conversion of NEU coordinates into LOS geometry (Eq. 7). It should be highlighted that in the proposed model, GNSS observations describe positions, while the radar observations are considered as velocity measurements. This approach allows us to preserve observations in their original domains: GNSS as an absolute estimate of coordinates and DInSAR as a relative LOS distance change over revisit time. The measurement noise matrix (\(R_t\)) contains the correlated variances of the GNSS topocentric coordinates (\(\sigma _{N_G}\), \(\sigma _{E_G}\), \(\sigma _{U_G}\), \(\sigma _{NE_G}\), \(\sigma _{NU_G}\), \(\sigma _{EU_G}\), Eq. 8), whereas estimates from DInSAR in ascending and descending geometries are assumed to be uncorrelated (\(\sigma _{d_{LOS}^A}, \sigma _{d_{LOS}^D}\), Eq. 6):

$$\begin{aligned} R_t= \begin{bmatrix} \sigma _{N_G}^2 &{} \sigma _{NE_G} &{} \sigma _{NU_G} &{} 0 &{} 0 \\ \sigma _{NE_G} &{} \sigma _{E_G}^2 &{} \sigma _{EU_G} &{} 0 &{} 0 \\ \sigma _{NU_G} &{} \sigma _{EU_G} &{} \sigma _{U_G}^2 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \sigma _{d_{LOS}^A}^2 &{} 0 \\ 0 &{} 0 &{} 0 &{} 0 &{} \sigma _{d_{LOS}^D}^2 \end{bmatrix}. \end{aligned}$$
(22)

The presented Kalman filter algorithm can be applied in a real-time mode and run anytime when a new observation is available. However, this method works only in a forward direction, so the detection and elimination of outliers that potentially occurred in the past is not possible. In order to get better estimates of the forward results, the backward Kalman filter could be applied (Fig. 4).

2.4.3 Kalman backward algorithm

The Kalman backward algorithm, named also as smoothing processing, relies on the results from the forward application of the Kalman filter, but can also be used parallel to a real-time filter Verhagen and Teunissen (2017). The backward filter runs recursively with \(t=\{N-1,N-2,\ldots ,0\}\), where N is the current epoch. The smoothed state values (\(\widehat{x^b}_{t|N}\)) and error variance–covariance matrix (\(P^b_{t|N}\)) are equal to:

$$\begin{aligned} \widehat{x^b}_{t|N}= & {} {\widehat{x}}_{t|t}+L_t(\widehat{x^b}_{t+1|N}-{\widehat{x}}_{t+1|t}), \end{aligned}$$
(23)
$$\begin{aligned} P^b_{t|N}= & {} P_{t|t}+L_t(P^b_{t+1|N}-P_{t+1|t})L_t^T, \end{aligned}$$
(24)

where

$$\begin{aligned} L_t=P_{t|t}\Phi _{t+1|t}^T P_{t+1|t}^{-1}. \end{aligned}$$
(25)

The Kalman backward filter is mainly used in the post-processing mode. The smoothing algorithm cannot be run separately—the results from the forward filter are required in every recursion. Finally, the statistical accuracy analysis was performed for the filtered state vectors of forward (\(N_F\), \(vN_F\), \(E_F\), \(vE_F\), \(U_F\), \(vU_F\)) and backward (\(N_B\), \(vN_B\), \(E_B\), \(vE_B\), \(U_B\), \(vU_B\)) Kalman filters with reference to the external data source (Sect. 2.5).

2.5 Quality analysis

To verify the obtained results, a quality analysis based on an independent data set was performed. Five GNSS campaign measurements (04.2019, 08.2019, 11.2019, 06.2020, and 01.2021) covering the time span and the study area of the DInSAR imagery and the GNSS permanent observations were used in the validation process. In order to co-locate the low-cost receivers with the nearest campaign points (see Fig. 2), a cross-reference was performed. Due to the non-identical locations between the epoch-based points and the permanent GNSS sites, the closest possible site was selected for verification. The distances between the chosen campaign points and the low-cost stations ranged from 25 to 80 m. The coordinates estimated by the Kalman backward filter for the first epoch of campaign measurements were considered as a reference and subtracted from each subsequent campaign measurement. The outcome of this procedure is a set of residual values for each of the campaign measurement epochs (except the first one).

Due to the insufficient number of samples (Chai and Draxler 2014), the calculation of any statistic based on five values was not robust. That is why another approach for validation was proposed. The differences in coordinates (\(e_i, i=1,2,\ldots ,n\)) between the campaigns at the measured sites are compared with the displacements for the corresponding intervals (n) from the analyzed data sets, namely DInSAR and GNSS (Sect. 2.3), Kalman forward and Kalman backward filter (Sect. 2.4). The results from the comparison are presented in Tables 2, 3, 4. Based on the calculated sets of differences (n) in NEU directions, the overall RMS errors for each of the methods were estimated:

$$\begin{aligned} RMS=\sqrt{\frac{1}{n}\sum _{i=1}^n e_i^2}. \end{aligned}$$
(26)

3 Results

In our study, we proposed two methods for combining DInSAR and GNSS data mentioned in Sects. 2.3 and 2.4. The effectiveness of each of them is observed by testing the quality of the results for the study case of Rydułtowy coal mine. The two first approaches based on a separate comparison are presented in Sect. 3.1—once only with the introduced N component of the displacement from the permanent GNSS measurements instead of the omitted one in the DInSAR 2-D decomposition (Sect. 2.3). The assessment of the reliability of the second method—the forward–backward Kalman filter assimilation, is reported in Sect. 3.2, and the methodology of estimation is presented in Sect. 2.4.

3.1 Unified DInSAR and GNSS time series

By applying the methodologies proposed in Sect. 2.3, we observed all permanent stations, chosen in this study (see Table 1). We compared their GNSS and DInSAR time series with the overlaid data from the GNSS epoch-based measurements (if available) in order to assess the efficacy of these two techniques. The results from the comparison for three low-cost stations: PI16 (Fig. 5.1—first row), PI03 (Fig. 5.2—second row), PI02 (Fig. 5.3—third row), located just above active extraction panels (see Fig. 2), and one geodetic permanent station RES1 (Fig. 5.4—fourth row), placed in the vicinity of the extraction panels, are presented below.

Fig. 5
figure 5

Displacement values of PI16 (1), PI03 (2), PI02 (3), and RES1 (4) stations, determined in North (a), East (b), and Up (c) directions. The graphs present the results of the GNSS permanent data (light blue dots), DInSAR data (dark blue dots), and campaign GNSS measurements (pink dots with error bars)

Full size image

The PI16 deformation, visible in Fig. 5.1, demonstrates that the epoch and permanent GNSS estimates follow a similar pattern—the North component (Fig. 5.1.a) shows displacement in the first 9 months of the observations and then stabilizes around 0.04 m. The East component (Fig. 5.1.b) is demonstrating similar 0.04 m deformation for the first 6 months of 2020 and stabilizes afterward. The Up component (Fig. 5.1.c) has a similar deformation cycle, with the first significant period lasting for 9 months in 2019 introducing 0.30 m vertical subsidence, while the second stage comprises the entire 2020 and shows rather smaller deformation reaching the level of \(-\) 0.40 m.

The GNSS observations are in disagreement with the DInSAR results for the East component (Fig. 5.1.b) with exceeding deformation, which is exhibited with a sharp inflexion in the very first months of the observations in 2019. Another discrepancy is introduced in 06.2020 shifting the East DInSAR deformation time series by 0.08 m, which can be associated with an unwrapping error. Such an error can be associated with strong decorrelation related to severe weather conditions or with greater than 1 cycle displacement. It has to be pointed out that the Up component (Fig. 5.1.c) shows much better alignment with GNSS but slightly underestimates overall deformation (0.35 m instead of 0.40 m).

A similar effect is observed in the results for station PI03 (Fig. 5.2), which is located above another extraction wall to the west of the main area of interest (Fig. 2). The estimated, from the GNSS permanent observation, deformation for the North component (Fig. 5.2.a) has two clearly defined sections—almost no deformation until 03.2020 and then a rapid change to \(-\) 0.20 m in the next year, until 03.2021. The East component (Fig. 5.2.b) shows marginal positive displacements in the first nine months of the measurements (03.2019–12.2019) up to 0.02 m, followed by a shift in the western direction (12.2019–03.2021). During the second period, some slow small variations in a positive direction are also detected, but the total site displacement at the end of the period in 03.2021 is recorded at \(-\) 0.04 m. The Up component (Fig. 5.2.c) shows, after a short 3-month non-deforming period (03.2018\(-\)06.2019), a fast pacing negative deformation accumulated to 1.00 m in 21 months (06.2019–03.2021).

The DInSAR results demonstrate a discrepancy in the East component in reference to the GNSS solution, starting from 03.2020, which coincides with the acceleration in the North component (Fig. 5.2.a), but also with the results for station PI16 (Fig. 5.1.b). The Up (Fig. 5.2.c) shows a similar DInSAR pattern in comparison with the GNSS trend—after a non-deforming period there is rapidly developing deformation accumulating up to 0.70 m sink at the end of 12.2020. It is 0.20 m less than what is visible in the GNSS data, most likely associated with the unwrapping problem in the DInSAR processing.

Station PI02 is located in the northern part of the extraction wall No.VII-E-E1 (see Fig. 2), and hence, the vertical deformation presented in Fig. 5.3.c appears in the very beginning of the period-05.2019 and is rather linearly incrementing until 09.2020 when it reaches \(-\) 0.20 m. There is an agreement between the permanent GNSS (light blue), the epoch GNSS (pink), and the DInSAR-based (dark blue) estimates. Interestingly, the East component from DInSAR (Fig. 5.3.b) demonstrates variations in the order of 0.04 m. However, the final East value is similar to the one from the GNSS observation (\(-\) 0.04 m).

Substantially different results, in comparison with the low-cost (epoch-based) GNSS observations, are derived from the geodetic permanent station RES1. As shown in Fig. 5.4, the completeness of the GNSS data is close to 100%. The coordinates variability is in the range of 0.001 and 0.002 m for the horizontal components (Fig. 5.4.a and 4.b), whereas for the vertical component it is between 0.003 and 0.004 m (Fig. 5.4.c). The accumulated displacement of North direction is 0.12 m (Fig. 5.4.a), and the Eastward movement comprises two periods of \(-\) 0.015 m and 0.020 m and finally reaches 0.0 m (Fig. 5.4.b), while the Up component drops to \(-\) 0.12 m (Fig. 5.4.c). The DInSAR results underestimate the vertical deformation by \(-\) 0.06 m and interestingly show similar variability in the Up component as in the East component (compare Fig. 5.4.c and Fig. 5.4.b), most likely again due to the unwrapping problem.

3.2 Kalman filtering application

The methodology for the application of the Kalman filter algorithm for DInSAR and GNSS integration is outlined in Sect. 2.4. However, before executing the estimation process, it was necessary to provide a priori knowledge related to the dynamics of the deformation. In the presented study, we used the zero-mean acceleration model to represent the noise of the system (Eq. 15). To provide the most optimal Kalman system noise value (\(\sigma _0\)), the RMS errors for 4 permanent stations were analyzed with respect to the GNSS epoch measurements. Among the nine tested levels (10, 5, 1, 0.5, 0.1, 0.05, 0.01, 0.005, and 0.001 mm/day\(^2\)), the lowest errors for both the Kalman forward and Kalman backward solutions were obtained for acceleration level at 0.05 mm/day\(^2\) with 89 and 91 mm RMS errors, respectively. For the Kalman forward filter, the largest RMS was calculated for \(\sigma _0=0.001\) mm/day\(^2\) (205  mm error), while for the Kalman backward filter the worst result equals the 138 mm error obtained for \(\sigma _0=10\) mm/day\(^2\).

The application of the methodology to the permanent GNSS station and DInSAR data shows significantly different results (Fig. 6) than DInSAR data and GNSS data alone (Fig. 5). It has to be also noted that the integrated results are available for all six components: N, E, U, \(v_N\), \(v_E\), \(v_U\); therefore, for each station, the results are presented in six panels—first row for the deformations and the second row for the velocities.

Fig. 6
figure 6

Displacement values (top rows) and rates (bottom rows) of the PI16 (1), PI03 (2), PI02 (3), and RES1 (4) stations, determined in North (a), East (b), and Up (c) directions. The charts present results of the GNSS permanent data (light blue dots), campaign GNSS measurements (pink dots with error bars), DInSAR data (dark blue dots), forward Kalman filter (green lines), and backward Kalman filter results (red lines)

Full size image

It is clearly visible in Fig. 6.1 that the Kalman forward model (green solid line) well traces the estimated PI16 displacements in the epoch solution of the North component (Fig. 6.1.a, top row). Similarly to Fig. 5.1.a, the positive deformation is estimated to be 0.04 m at the last period of 03.2021. It is worth to note that the GNSS permanent station demonstrates periods with noisy results in 06.2020. However, it does not have an impact on the Kalman filter, which produces stable outputs. Both East and Up components (Fig. 6.1.b and 1.c, top row panels) resulting from the forward filter show a strong linear extrapolated trend from 06.2019 until 12.2019. The velocities (Fig. 6.1, bottom row) that are taken from the DInSAR observations (dark blue dots) show very noisy output (Fig. 6.1.b and 1.c). It is interesting to observe the fast re-convergence of the Kalman estimates to the GNSS results once the GNSS data are available again (after the data gap 06.2019–12.2019). The Kalman derived displacements align well with the epoch solution except for the last height values for 12.2020 (Fig. 6.1.c, top row), which shows \(-\)0.35 m instead of \(-\) 0.40 m. The backward Kalman filter applied as a smoothing solution aligns well with the GNSS observations in all three components (Fig. 6.1, top row panels), tracing the permanent GNSS solution.

Figure 6.2 demonstrates the time series of the fast deforming station PI03. The overall Kalman forward solution (green solid line) is tracing the GNSS observations that transpire as linear extrapolation sections visible in the long GNSS data gaps (06.2019–12.2019; 04.2020–06.2020). However, a quicker return to the tracing data trend is notable in the case of the East component (Fig. 6.2.b, top row) on 10.2019, two months before the GNSS data reappear in the solution. At the same time, the velocity calculated from the DInSAR data is much less noisy than the rest of the time series (Fig. 6.2.b, bottom row). It is also remarkable that the velocities of the Up component (Fig. 6.2.c, bottom row) show significant negative DInSAR values with noise low enough to make the trend estimation viable (up to 2 mm/day). The smoothed solution (backward Kalman filter-red line) shows no significant jumps or breaks, tracing the GNSS estimates closely.

The relatively slow deformation process at station PI02 with the final accumulated vertical displacement of \(-\) 0.20 m is manifested in Fig. 6.3.c. The Kalman forward filter (green line) shows dependence on the GNSS data (light blue dots). However, the estimated trends in all components—North, East, and Up (Fig. 6.3, top row panels) indicate a high level of agreement with the epoch measurements (pink bars). It is worth to mention that the East and Up components are visible using information derived from the DInSAR observations (Fig. 6.3.b and 3.c, bottom row) that happens to be less noisy in winter (09.2019–03.2020). Similarly to the results of the previous sites, Kalman filter outputs are not sensitive to the outliers in the GNSS results (light blue dots in the top row—close to 06.2020).

Figure 6.4 presents the Kalman filter estimates for the RES1 station that is away from the major deformation zone and for which almost the entire GNSS series is completed. It is clear that the Kalman filters, both forward and backward, are tracing the GNSS observation curve, while the DInSAR observations for this particular point are not as critical (no GNSS data gaps). Moreover, evaluating the bottom panels of Fig. 6.4, it is clear to see that the breaks in 2018 and 2019 brought a lot of noisy data that were difficult to process (see Fig. 5.4). Surprisingly, the Up component (Fig. 6.4.c) demonstrates sections of different deformation velocities aligned well with the East component inflection points (03.2019 and 12.2019). The overall quality of Kalman retrieval at the RES1 station could be taken as a reference for other sites.

3.3 Intersecting results

The results presented in the previous sections demonstrate capabilities for deformation monitoring by DInSAR and GNSS data separately, or in an integrated concept based on the Kalman filter algorithm. The subsequent validation is based on the residual values for the GNSS permanent sites located in the vicinity of the GNSS campaign points. The choice of the residual analysis against statistical evaluation is dictated by the small number of low-cost stations (only four) in the vicinity of campaign sites, and only five overlapping dates within the time of the permanent and campaign points installation. The low-cost stations used in the intersect comparison are PI02, PI04, PI05, and PI16, with residuals in North (Table 2), East (Table 3), Up (Table 4) directions, calculated as a difference with the values for these three components measured at the closest points of the GNSS campaign network. What can be definitively observed in all Tables is that the RMS (Eq. 26) for the GNSS observations is similar to the Kalman forward and backward model, regardless of the station and site.

Table 2 North residuals for low-cost permanent stations in the vicinity of campaign points—see Fig. 2, for location
Full size table

In Table 2, the total deformation in the North direction is relatively small and does not exceed 0.10 m (point PI02). The overall RMS for all calculated residuals is 0.013 m regardless of the technique applied to obtain these deformations. It is worth to notice that Kalman forward and backward demonstrates availability for all reference epochs.

Table 3 East residuals for low-cost permanent stations in the vicinity of campaign points—see Fig. 2, for location
Full size table

In Table 3, the total deformation in the East direction does not exceed 0.11 m (point PI04). The discrepancies between campaign points and low-cost observations do not exceed 0.09 m; however, the uncertainty of the campaign points is 0.02 m. The largest residuals are recorded for the DInSAR result up to the maximum of \(-\) 0.052 m (PI04, 01.2021). It is also worth to mention that the Kalman filter provides the best quality solution across all reference epochs regardless of the presence or absence of GNSS data (see PI02, PI05).

Table 4 Up residuals for low-cost permanent stations in the vicinity of campaign points—see Fig. 2, for location
Full size table

In Table 4, the Up direction deformation is substantial with amounts up to \(-\) 0.388 m (point PI16), and discrepancies between campaign points and low-cost observations up to 0.16 m. Note that the uncertainty of campaign points is 0.03 m. The largest residuals are recorded again for the DInSAR results up to the maximum of \(-\) 0.142 m (PI04, 06.2020). For the PI04 station, the distance to the closest campaign point is equal to 80 m and the residuals are 142, 77, 91, and 91 mm for InSAR, GNSS, Kalman forward, and Kalman backward, respectively. It is also worth to note that the Kalman filter provides a solution across all reference epochs regardless of the completeness of the GNSS data. The residuals for PI04 are higher than half of the deformation, ranging from 0.003 m to 0.142 m, whereas the maximum recorded deformation was \(-\) 0.201 m. In the case of PI16, we may observe growing residuals with increasing deformation (maximum \(-\) 0.388 m).

4 Discussion

In our study, we developed a novel Kalman filtering approach combining high-temporal resolution displacement time series (GNSS) with a high-spatial resolution measurements (InSAR) to provide consistent 3-D coordinates and velocities for the location of the GNSS stations. It is classified as one of the most promising solutions postulated by Hu et al. (2014) for local and regional scale deformations. In the presented study, we used data from six permanent GNSS stations for a relatively small extraction field (2 \(\textrm{km}^2\)). For comparisons, in Bian et al. (2014) 19 receivers were distributed over 8000 \(\textrm{km}^2\) large-scale coal mining area in China with an average distance between the monitoring points of about 35 km, while in Liu et al. (2019), the number of GNSS sites is 35 on over 2000 \(\textrm{km}^2\). As the local nonlinear strong deformations were investigated in a relatively small area with quite large number of GNSS receivers, we can identify and correct the number of DInSAR processing issues that limit the quality of deformation estimates.

The DInSAR time series of displacements can be affected by outlier values related to the limitations of the technique, e.g., decorrelation in vegetated areas, local atmospheric effects, or other phase unwrapping problems (Osmanoğlu et al. 2016; Fattahi and Amelung 2014). The unwrapping procedure is a non-deterministic problem that has many equally correct solutions and can only be solved under certain assumptions. Nevertheless, the high subsidence rate provides a strong interferometric signal. Thereby, the fringe rate observed in interferograms increases, making it more difficult to unwrap properly (Osmanoğlu et al. 2016). In the presented study, the largest discrepancy between the GNSS and DInSAR results occurred for the PI03 station (max. disagreement equal to 0.30 m for the Up component), which corresponds directly to the most prominent observed vertical rate of deformation (\(-\) 1.00 m over 2 years of GNSS monitoring).

Another major issue related to the transformation of ground deformations measured in the DInSAR slant direction to LOS-decomposed results can introduce a significant distortion into the three-dimensional representation of the displacement depending on the applied method. As the authors presented in Wright (2004), the a posteriori Up errors decreased from 286 to 4 mm after neglecting the North component, which is the result of constraining the null space of the solution (Brouwer 2021). In our study, in the LOS decomposition into East and Up components, we applied the North displacement from the GNSS data sets (Eq. 10). To test the impact of the North component magnitude on the presented results, we introduced two separate decomposition processes (with and without North) for the PI03 station, which is characterized by the largest GNSS North deformation (\(-\) 0.20 m). Taking into account the cumulative DInSAR ascending (\(-\) 0.63 m) and descending (\(-\) 0.45 m) displacements, the Up component becomes \(-\)0.72 m after decomposition without North (North component equals 0), while with North is \(-\) 0.77 m (North component from GNSS). Therefore, the impact of the North GNSS information is equal to 6% (0.05 m). For the East deformation, the 0.106 m is obtained after decomposition without North (North component equals 0), and 0.101 with North (North component from GNSS); therefore, the impact is 5% (0.005 m). However, the final Kalman filter solution is fully constrained in all components, which results in the same accuracy for GNSS and Kalman filter solution for the North component (0.013 m) (see Table 2).

We successfully overcome the geometrical and technological limitations of InSAR, but also we had faced challenges related to GNSS processing. One of the issues with low-cost sensors is the incomplete coordinate time series recorder by the GNSS receivers. In the presented literature review for studies using permanent GNSS observations (Fuhrmann et al. 2015; Liu et al. 2019; Tao et al. 2019), there were no significant interruptions (longer than a month). However, in the current study, as we used low-cost receivers, we observed up to nine months of breaks in GNSS data flow. Such data gaps have a great impact on the continuous monitoring of nonlinear deformations. The missing GNSS time spans can be filled by DInSAR data using a Kalman filter approach (Eqs. 20–22). However, calculations provided only in the forward mode can introduce significant leaps and overestimates in the time series. Such an example is visible for the PI16 station around 12.2019 Fig. 6.1.a) in North direction or PI02 station around 04.2020 in the Up direction Fig. 6.3.c). Hence, to smooth and eliminate unexpected discontinuity of the data, we introduced a backward Kalman filter algorithm (Eqs. 23–25). The applied backward filter significantly reduced the impact of gaps in the GNSS data. In our case, we have run the backward filter only at the last epoch; however, subsequent backward filter runs are possible and would reduce the systematic error introduced by GNSS data gaps.

Finally, one of the significant problems that many deformation studies face is validation. In Bozso et al. (2021), the a posteriori errors of the Kalman filter algorithm were assessed to be 16, 2, and 3 mm in North, East, and Up directions, respectively, given the standard deviation of GNSS and InSAR \(d_{LOS}\) to be 2 mm a priori. No input noise estimation was performed, and no external data were applied for validation. The fusion of GNSS and InSAR in the paper of Liu et al. (2019) generated spatiotemporal deformation time series, however, the distribution of GNSS receivers and the sampling interval of InSAR interferograms had an impact on the integration accuracy. The quality of the deformation field model, established using the leave-one-out method, did not exceed 7.2 mm in the LOS direction. In contrast to the above, in the current paper, the verification procedure was performed using the external data source - GNSS campaign measurements. Due to the non-identical location of the campaign points and permanent GNSS stations, the closest possible locations were selected for verification. The most significant residual values are presented in Tables 2, 3, 4 and correspond to the longest distance between these points (station PI04 and the closest campaign point at 80 m distance). For the Kalman backward approach, the maximum values are -38, 34, and 91 mm for N, E, and U, respectively. The overall RMS values of all stations by the Kalman backward filter are equal to 13, 17, and 34 mm for N, E, and U, respectively. These values are well below the campaign measurement uncertainties where the maximum values reach 102 mm, 75 mm, and 98 mm for N, E, and U, respectively.

5 Conclusions

The paper presents an original methodology for the integration of two different techniques, optimal for nonlinear motions, conducted for an area affected by underground mining works. The algorithm is able to ingest the noisy GNSS topocentric coordinates with significant gaps, as well as non-corrected for the troposphere and unwrapping errors DInSAR ascending and descending velocities. The observation uncertainties are rigorously determined for GNSS during the parameter estimation process, and for DInSAR error calculated from coherence coefficient values.

The presented methods are tested and verified for the “Rydułtowy” mine (2 extraction walls with a total area of \(2~\textrm{km}^2\)), located in the south-western part of Upper Silesia in Poland, where multilevel coal exploitation induces rapid ground deformations up to 1.00 m over 2 years. The cross-reference of the data sets is presented for low-cost GNSS stations PI16, PI03, PI02, and geodetic GNSS station RES1 (Figs. 5, 6). The nearest GNSS campaign points to the permanent stations PI02, PI04, PI05, and PI16 were used for validation (Tables 2-4). The lowest overall RMS errors were reached for the Kalman backward approach (13, 17, and 34 mm for the N, E, and U, respectively), whereas the highest values were obtained for the sole use of the DInSAR technique (24 and 52 mm for the E and U, respectively). The missing North displacements data for the DInSAR are a result of the limited sensitivity of radar observations to the deformations in the along-track direction.

It can be concluded that the introduced Kalman forward and backward methods of integration provide around 1.5 times better results for the East and Up components of the deformation vector than the standard DInSAR decomposition alone. Simultaneously, the RMS uncertainties of Kalman filters in the case of the North displacements are on the same level as for GNSS data sets (13 mm). What is worth to note, the overall RMS errors were calculated based on residuals from only common epochs of all methods (e.g., the measurements from 08.2019, 11.2019, and 01.2021 have to be omitted due to technical problems of PI02 and PI05 GNSS stations, see Tables 2-4). One of the most significant advantages of the Kalman fusion, compared to the standard GNSS or DInSAR solution, is the completion of potential breaks in the time series of observations (see the GNSS data in Fig. 6).

Another benefit of integration is the preservation of the data in their original domains: GNSS as an absolute 3-D estimate of coordinates and DInSAR as a relative 1-D LOS distance change over revisit time. Moreover, the proposed measurement-update model allows us to extract separate NEU parts from ascending and descending DInSAR observations avoiding the decomposition process (see Eq. 21). Otherwise, using the standard decomposition approach, to obtain the absolute DInSAR time series, the conversion from relative to the cumulative domain is required. The incorrect identification and elimination of any outlier may bias subsequent epochs in the data set (see Fig. 5). The difficulty of removing effects of the local atmospheric conditions and other phase unwrapping problems can be reduced by thresholding of system noise level in the Kalman filter time-update segment (see Eq. 15). Subsequently, the fusion processing can be easily extended (adding additional rows in Eqs. 20, 21, 22) by data from other sources, e.g., leveling, or LiDAR, to create a framework for multi-sensor ground deformation monitoring system.

In the further research steps, the knowledge gained about radar phase discontinuities in co-located GNSS-DInSAR points can be extended to surrounding pixels. The developed methodology of integration creates the foundation which can allow nonlinear areal modeling. Currently, numerous empirical models and influence function models can be used to predict mining-induced subsidence, e.g., the Probability Integral Method, the Logistic Model, or Knothe theory. The Knothe model is a commonly applied method in the case of coal mining deformation monitoring and can be adopted in the integration approach together with improving the fusion process by other ground monitoring techniques.