Near Real-Time Monitoring of Large Gradient Nonlinear Subsidence in Mining Areas: A Hybrid SBAS-InSAR Method Integrating Robust Sequential Adjustment and Deep Learning

Abstract

With the increasing availability of satellite monitoring data, the demand for storage and computational resources for updating the results of monitoring the surface subsidence in a mining area continues to rise. Sequential adjustment (SA) models are considered effective for rapidly updating time series interferometry synthetic aperture radar (TS-InSAR) measurements. However, the accuracy of surface subsidence values estimated through traditional sequential adjustment is highly sensitive to abnormal observations or prior information on anomalies. Moreover, the surface subsidence associated with mining exhibits nonlinear and large gradient characteristics, making general InSAR methods challenging for obtaining reliable monitoring results. In this study, we employ the phase unwrapping network (PUNet) to obtain unwrapped values of differential interferograms. To mitigate the impact of abnormal errors in the near real-time small baseline subset InSAR (SBAS-InSAR) sequential updating process in mining areas, a robust sequential adjustment method based on M-estimation is proposed to estimate the temporal deformation parameters by using the equivalent weight model. Using a coal backfilling mining face in Shanxi, China, as the study area and the Sentinel-1 SAR dataset, we comprehensively evaluate the performance of unwrapping methods and subsidence time series estimation techniques and evaluate the effect of filling mining on surface subsidence control. The results are validated using leveling measurements within the study area. The relative error of the proposed method is less than 5%, which can meet the requirements of monitoring the surface subsidence in mining areas. The method proposed in this study not only enhances computational efficiency but also addresses the issue of underestimation encountered by InSAR methods in mining area applications. Furthermore, it also mitigates unwrapping phase anomalies on the monitoring results.

1. Introduction

Coal as a crucial primary energy source serves as a strong pillar for social development. With the increasing demand for energy, coal mining not only creates significant value for the world economy but also raises concerns about associated hazards. The subsidence and collapse resulting from large-scale and sustained mining operations not only damage surface infrastructure but also have the potential to trigger various geological disasters, posing a significant threat to the safety of people’s lives and property. Therefore, to prevent disasters and ensure safe production, it is imperative to monitor the subsidence in mining areas.

Differential Interferometric Synthetic Aperture Radar (D-InSAR) is a technology that utilizes SAR images from different time periods to obtain surface deformation information through differential interferometry. This technique, known for its cost-effectiveness and efficiency, is widely applied in monitoring the surface subsidence in mining areas [1,2]. Phase unwrapping is a crucial step in the processing of InSAR data, and the unwrapped phase directly influences the accuracy of the subsidence parameter estimation [3,4]. Owing to the periodic nature of the phase data with its inherent 2π ambiguity, phase unwrapping is employed to restore the absolute phase values, thereby yielding more accurate deformation information [5]. In mining scenarios, significant nonlinear subsidence with large gradients often arises as the mining face advances. Conventional phase unwrapping methods, such as branch–cut and minimum cost flow, tend to underestimate the subsidence with large gradients, leading to biases in the subsequent subsidence parameter estimates [6]. Therefore, it is essential to employ unwrapping methods suitable for large gradient subsidence. In recent years, scholars have applied deep learning methods to phase unwrapping, achieving considerable success [7]. Interferograms depicting large gradient subsidence exhibit complex fringe patterns and significant decorrelation noise. However, specific features or patterns within interferograms can still be extracted through deep learning to enhance the phase unwrapping results. To address the challenges of unwrapping the large gradient subsidence in mining areas, this study employs the PUNet for phase unwrapping operations. The PUNet treats phase unwrapping as a regression problem, learning a direct mapping from the interferograms to the unwrapped phases. This approach can effective handle deformation the interferograms characterized by dense fringes and decorrelation noise.

Traditional D-InSAR techniques suffer from interference factors, such as temporal and spatial decorrelation, and atmospheric effects, leading to a lower measurement accuracy [8,9]. TS-InSAR technology, by processing historical SAR images of the same area, can mitigate or eliminate the impact of the aforementioned factors, resulting in high-precision deformation parameters. SBAS-InSAR is currently a widely applied TS-InSAR method, playing a crucial role in studies related to ground subsidence, landslides, and other scenarios [10,11]. The SBAS-InSAR employs short baseline interferometric combinations for the measurements, effectively mitigating the decorrelation issues caused by spatial baselines. Additionally, since all interferograms satisfying the short baseline condition are involved in the estimation and the singular value decomposition method is employed, the method can obtain the minimum norm solution of the observation time series, thereby improving the stability and reliability of the solution. Therefore, as a post-processing technique for D-InSAR, SBAS-InSAR is more practical and reliable.

In the mining area, rapid and near real-time updates of surface subsidence monitoring values are of significant importance for guiding the mining progress and for disaster prevention and mitigation. With the development of the measurement technology, SAR satellites such as the European Space Agency’s Sentinel-1 operate continuously. This enables the stable acquisition of SAR data, providing a rich data basis for the InSAR time series analysis for near real-time monitoring in mining areas. However, a challenge arises with the increasing volume of SAR data processing, which is to dynamically update the deformation parameters without compromising the accuracy of the parameter estimation while also improving the computational efficiency [12]. General SBAS-InSAR typically involves processing a certain amount of data in a region to obtain deformation monitoring results for the accumulated time period. When new SAR images are added, each update calculation requires merging the previous interference phase and date matrices with the newly added data, thereby expanding the dimensions of the matrices, which hinders the efficient, continuous, and dynamic monitoring requirements [13,14].

To address this challenge, numerous scholars have proposed various solutions. Ansari et al. introduced an efficient InSAR time series algorithm. This method divides the data into small datasets, compresses the datasets, and avoids reprocessing the entire dataset, greatly improving the time efficiency [15]. Hu et al. explored a general pattern of surface deformation by applying a sequential adjustment method after adding two or more SAR images, validating the feasibility of the approach with simulated data [9]. Wang et al. conducted research on sequential adjustment and dynamic Kalman filtering methods in the efficient computation of SBAS-InSAR incremental data [16].

Sequential adjustment requires no storage of a large amount of historical observational data and can obtain the same optimal solution as the overall adjustment without the need for large matrix inversions. Currently, standard sequential adjustment is widely adopted for dynamically solving SBAS-InSAR deformation parameters efficiently. Existing research often assumes that the points do not change after the addition of new data, representing the ideal scenario of consistently effective monitoring points. However, in practical engineering applications, it is challenging to meet this assumption. In the process of InSAR monitoring, changes in land features may occur, leading to variations in the coherence of objects depicted in SAR images. The decreased coherence of land features may result in anomalous values in phase unwrapping, and when features lose coherence, effective monitoring values cannot be obtained [17].

Furthermore, with the increase in temporal baselines, whether using conventional unwrapping methods or the PUNet, the accuracy of phase unwrapping decreases, leading to the presence of outliers in the results. Therefore, when dynamically solving SBAS-InSAR deformation parameters, it is essential not only to consider the efficiency of computation but also to account for anomalies and invalid values caused by changes in the coherence of land features. This consideration helps avoid introducing outliers in the solution results. In summary, there are two key concerns regarding InSAR technology in the context of mining-induced subsidence: one pertains to reliable phase unwrapping methods, while the other concerns the efficient and accurate computation of temporal subsidence. This study employs the PUNet for the interferogram unwrapping to obtain the unwrapped phases of the large gradient deformation. To overcome the influence of outliers caused by the intermittent scattering of ground objects, a robust estimation theory is introduced into the SBAS-InSAR parameter estimation. Robust estimation methods are utilized to estimate the initial deformation parameters of the SBAS-InSAR, providing reliable a priori estimates. Subsequently, a combination of the least squares Bayesian estimation and robust sequential adjustment theory is employed for the subsidence parameter estimation. During the dynamic parameter estimation, only the archived data related to the pairs of new data are processed. Finally, based on the results of the InSAR monitoring, we also evaluated the effect of the filling mining method on controlling the surface subsidence in the study area. The primary innovation of this study lies in the introduction of robust estimation and sequential computation into InSAR monitoring in mining areas. Compared to slow subsidence, such as urban land subsidence, the subsidence induced by mining operations occurs at a faster rate and requires a higher monitoring frequency, necessitating more efficient monitoring and computation methods.

2. Methodology

Prior to the estimation of the deformation parameters, preprocessing steps are applied to the SAR data, involving the following procedures:

(1) Utilization of time and spatial baseline constraints to generate M interferograms from the adjacent N archived datasets. An External Digital Elevation Model (DEM) was employed for the topographic phase removal, and the flat Earth phase was corrected based on a reference ellipsoid and imaging parameters. The DEM utilized in this study comprised the data collected by unmanned aerial vehicles with a resolution of 10 m.

(2) Selection of a common location as a reference point, followed by phase unwrapping for the points meeting coherence threshold criteria. Based on the characteristics of the PUNet and the geomorphological features of the study area, the candidate points with a coherence greater than 0.2 were selected.

(3) Correction of orbit errors, turbulent atmospheric effects, and vertical stratified atmospheric effects. In the atmospheric correction, a phase-to-height linear model was first employed to correct for vertical stratified atmospheric effects, followed by the application of the Goldstein filtering method to correct for turbulent atmospheric effects [18].

2.1. Introduction to PUNet

The PUNet was first proposed in reference [19]. The design of the PUNet model was inspired by classical denoising models such as the Denoising Convolutional Neural Network (DnCNN) [20], residual networks, and dilated convolutions. The DnCNN has proven effective in image denoising, avoiding the downsampling process. As a network backbone, it can effectively handle the SAR interferograms with dense fringes and decorrelation noise. Simultaneously, dilated convolutions significantly increase the network’s receptive field, and residual modules expedite the training process.

The PUNet models the subsidence in mining areas by utilizing a distorted two-dimensional Gaussian surface. Additionally, by adjusting certain parameters, it generates different subsidence signals with varying amplitudes and patterns [19]. The two-dimensional elliptical Gaussian function can be expressed as follows:

$$f(X)=\frac{1}{2\pi {\left|\Sigma \right|}^{1/2}}\exp \left\{\left.-\frac{1}{2}{(X-u)}^T{\Sigma }^{-1}(X-u)\right\}\right.$$

(1)

Here, $X=(x_1,x_2)$ represents a two-dimensional grid of a training sample, while $u=(u_1,u_2)$ controls the position of the subsidence center. The covariance matrix is expressed as follows:

$$\Sigma =s\times U^{\prime }DU$$

(2)

It is used to control the shape and size of the subsidence region, where $D$ represents a two-dimensional random diagonal matrix and $U$ represents an orthogonal basis of another two-dimensional random matrix, while $s$ is the scaling factor [19]. As shown in Figure 1, the PUNet model structure does not include the downsampling layers, and all convolutional layers use zero padding at the edges, maintaining an output size equal to the input size. Due to its fully convolutional architecture, the PUNet can accommodate input samples of arbitrary sizes.

2.2. Initial Parameter Robust Estimation

After removing the topographic residual phase from the unwrapped phase during the phase unwrapping, the following model can be established to solve the time series deformation of the archived data [12]:

$$\underset{A_1}{\left[\begin{array}{ccccccc}-1& 1& 0& \dots & 0& 0& 0\\ -1& 0& 1& \dots & 0& 0& 0\\ .& .& .& .& .& .& .\\ 0& 0& 0& \dots & 0& -1& 1\end{array}\right]}\underset{X}{\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\\ \dots \\ {\phi }_N\end{array}\right]}=\underset{L_1}{\left[\begin{array}{c}un{w}_1\\ un{w}_2\\ \dots \\ un{w}_M\end{array}\right]}$$

(3)

In the above equation, $\phi$ represents the accumulated deformation phase and $unw$ represents the unwrapped phase from interferogram. The equation can be expressed in the following form:

$$\begin{gathered} V_1=A_{1}X-L_1 \\ {X}_1={(A_1^T{A}_1)}^{-1}{A}_1^T{L}_1 \end{gathered}$$

(4)

where $V_1$ represents the measurement noise. To mitigate the potential impact of phase outliers, we employed a robust M-estimator to estimate the final parameters [21,22]. The specific process was as follows:

a. The initial weight matrix $\stackrel{-}{P}$ was an identity matrix:

$${\stackrel{-}{P}}_1^{(l=0)}=I$$

(5)

b. There was then a weighted least squares adjustment estimation.

$$X_1^l={(A_1^T{\stackrel{-}{P}}_1^l{A}_1)}^{-1}{A}_1^T{\stackrel{-}{P}}_1^l{L}_1$$

(6)

c. The residual phase matrix $v_1^l={[v_{11}^l,\dots ,v_{1N}^l]}^T$ was computed, and the equivalent weight matrix ${\stackrel{-}{P}}_1^{(l+1)}=\mathrm{diag}\left\{p_{1i}^{(l+1)}\right\}$ was as follows [23]:

$$v_i^l=L_1-A_1{X}_1^l$$

(7)

$$p_{1i}^{(l+1)}=\left\{\begin{array}{c}{p}_{1i}^l,{|\mathrm{V}}_i|\le k_0\\ p_{1i}^l\frac{k_0}{{|\mathrm{V}}_i|}{(\frac{k_1-{|\mathrm{V}}_i|}{k_1-k_0})}^2,k_0\le {|\mathrm{V}}_i|\le k_1\\ 0,{|\mathrm{V}}_i|\ge k_1\end{array}\right.$$

(8)

In Equation (8), $|V_i|=|\frac{v_{1i}^l}{{\sigma }_{v_i^l}}|$, $k_0=1,k_1=2.5$.

d. The computation was terminated when the residual phase converged; otherwise, $l=l+1$ was set, and the iteration was continued by returning to b. Here, we set the termination condition as the cessation of the robust estimation algorithm when there is no change between consecutive iterations. Finally, ${\stackrel{-}{Q}}_{X_1}={(A_1^T{\stackrel{-}{P}}_1{A}_1)}^{-1}$.

In the aforementioned process, the robust estimation iteratively reallocates the smaller weights to the larger residual phases, making the estimation less susceptible to the influence of phase outliers.

2.3. Robust Sequential Adjustment Method

When obtaining new SAR images, unlike the general SBAS-InSAR structure, we only considered the unwrapped interferograms related to the new SAR image. Sequential adjustment was employed to dynamically update the deformation time series. The observation equation can be reformulated as follows:

$$V_2=\left[\begin{array}{cc}{A}_2& B\end{array}\right]\left[\begin{array}{c}{X}_2\\ Y\end{array}\right]-L_2$$

(9)

where $X_2,Y$ represent the cumulative deformation variables for the archived time corresponding to the new data pair and the cumulative deformation variable for the new time, respectively. Current methods typically involve computations with all archived data, which can be computationally intensive. Without compromising accuracy, we only processed the time period corresponding to the new data pair. The parameters of the prior matrix do not need to be reprocessed; only truncating the existing matrix is necessary, minimizing computational resources.

The solutions for variables $X_2$ and $Y$ based on the Bayesian least squares and sequential adjustment with varying parameters are as follows [24]:

$$\begin{gathered} \left[\begin{array}{c}{X}_2\\ Y\end{array}\right]={\left[\begin{array}{cc}{\stackrel{-}{Q}}_{X_1}^{-1}+A_2^T\stackrel{-}{P_2}{A}_2& A_2^T\stackrel{-}{P_2}{A}_2\\ B^T\stackrel{-}{P_2}{A}_2& B^T\stackrel{-}{P_2}B\end{array}\right]}^{-1}\ast \left[\begin{array}{c}{\stackrel{-}{Q}}_{X_1}^{-1}{X}_1+A_2^T\stackrel{-}{P_2}{L}_2\\ B^T\stackrel{-}{P_2}{L}_2\end{array}\right] \\ {\stackrel{-}{Q}}_{\left[X_2;Y\right]}={\left[\begin{array}{cc}{\stackrel{-}{Q}}_{X_1}^{-1}+A_2^T\stackrel{-}{P_2}{A}_2& A_2^T\stackrel{-}{P_2}{A}_2\\ B^T\stackrel{-}{P_2}{A}_2& B^T\stackrel{-}{P_2}B\end{array}\right]}^{-1} \end{gathered}$$

(10)

$\stackrel{-}{P_2}$ represents the equivalent weight matrix, and for the purpose of enhancing the computational efficiency, ${({\stackrel{-}{Q}}_{X_1}^{-1}+A_2^T\stackrel{-}{P_2}{A}_2)}^{-1},{\left[\begin{array}{cc}{\stackrel{-}{Q}}_{X_1}^{-1}+A_2^T\stackrel{-}{P_2}{A}_2& A_2^T\stackrel{-}{P_2}{A}_2\\ B^T\stackrel{-}{P_2}{A}_2& B^T\stackrel{-}{P_2}B\end{array}\right]}^{-1}$ can be simplified using the following matrix inversion formula:

$${(P_X\pm A^{T}PA)}^{-1}=P_X^{-1}\mp P_X^{-1}{A}^T{(P^{-1}\pm A{P}_X^{-1}{A}^T)}^{-1}A{P}_X^{-1}$$

(11)

The solution for the final deformation parameters can be expressed as follows:

$$\begin{gathered} \left[\begin{array}{c}{X}_2\\ Y\end{array}\right]=\left[\begin{array}{c}{X}_1+J_x(L_2-A_2{X}_1-BY)\\ {(B^T{{\stackrel{-}{Q}}_J}^{-1}B)}^{-1}{B}^T{{\stackrel{-}{Q}}_J}^{-1}(L_2-A_2{X}_1)\end{array}\right] \\ {\stackrel{-}{Q}}_{\left[X_2;Y\right]}=\left[\begin{array}{cc}{\stackrel{-}{Q}}_{X_2}& {\stackrel{-}{Q}}_{X_2,Y}\\ {\stackrel{-}{Q}}_{X_2,Y}& {\stackrel{-}{Q}}_Y\end{array}\right] \end{gathered}$$

(12)

In the above equation, $J_x$ represents the sequential adjustment gain matrix, and each element can be specified as follows:

$$\begin{gathered} {\stackrel{-}{Q}}_{X_2}={\stackrel{-}{Q}}_{X_1}-J_x{A}_2{\stackrel{-}{Q}}_{X_1}+J_{x}B{\stackrel{-}{Q}}_Y{B}^T{J_x}^T \\ {\stackrel{-}{Q}}_{X_2,Y}=-J_{x}B{\stackrel{-}{Q}}_Y \\ {\stackrel{-}{Q}}_Y={(B^T{{\stackrel{-}{Q}}_J}^{-1}B)}^{-1} \\ {J}_x={\stackrel{-}{Q}}_{X_1}{A_2}^T{{\stackrel{-}{Q}}_J}^{-1} \\ {\stackrel{-}{Q}}_J={\stackrel{-}{P_2}}^{-1}+A_2{\stackrel{-}{Q}}_{X_1}{A_2}^T \end{gathered}$$

(13)

Similar to the robust estimation method mentioned in the initial deformation parameters, the iterative computation process denoted by steps a–d was executed. However, the weights were now represented by $\stackrel{-}{P_2}$, and the unknowns to be solved were $X_2$ and $Y$. The introduction of robust estimation in the sequential update can reduce the influence of new anomalous observations. This ensures consistency in the points and epochs involved in the sequential adjustment, while also considering the variations in the scattering characteristics of some points in practical scenarios. This approach suppresses the influence of outliers and invalid values on the results. Furthermore, the proposed method in this study reduces the storage space and computational time for matrix inversion compared to the general SBAS-InSAR method. The flowchart of the above method is shown in Figure 2.

3. Experimental Test

3.1. Simulation and Real Data

In our study, we employed a baseline network as depicted in Figure 3 below to simulate the unwrapped phase of 38 interferometric pairs over 14 epochs. With reference to the methods in previous studies [14,25], the initial nine epochs were designated as a priori observations, while the subsequent five epochs were considered as additional observations. The data from the eighth and eleventh epochs were intentionally introduced with poor coherence, thereby introducing outliers with gross errors of 5–10 mm on the relevant interferometric pairs. To assess the applicability of the algorithm, we conducted simulations of the evolution of the subsidence using both a linear model and a Weibull model, which exhibits similarities to the subsidence patterns in mining areas.

This paper employs the proposed robust sequential adjustment method for the dynamic monitoring of a backfilled coal mining face in Shanxi, China. The study area is characterized by mid-elevation topography dominated by limestone and sandstone formations, with relatively gentle slopes. The working face extends 1709 m in strike length and 176 m in dip length. The working face advances in a northeast-to-southwest direction. The coal seam being mined is structurally complex, classified as a stable, moderate thickness seam with a Puri hardness index of 2.0. The thickness of the coal seam ranges from 1.66 to 2.69 m, averaging 2.25 m, and dips at angles varying between 1° and 8°, with an average of 2°. Geological investigations have revealed multiple challenges during extraction, including the presence of faults, karst sinkholes, and hydrogeological issues. A total of 46 ascending track Sentinel images were utilized, resulting in the generation of 145 interferometric pairs. The central incidence angle of the data was 38.4°, with a multilook resolution of 20 m. The satellite data coverage and location of the study area are shown in Figure 4. A measured DEM with a resolution of 10 m was employed to eliminate the terrain-induced phase and atmospheric effects. To assess the suitability of a single-track satellite for the study area, the visibility of the satellite was computed based on the DEM and satellite imaging parameters. The calculation results from Figure 4 indicate a favorable visibility within the study area, with no overlapping or shadowed regions. The study area features gentle slopes, resulting in minimal displacements along the slope direction. Consequently, the vertical subsidence induced by mining activities constitutes the primary component of the satellite line-of-sight displacement.

Utilizing a network, as illustrated in Figure 5, the data spanning from November 2021 to June 2023 were divided into four monitoring phases: November 2021 to March 2022, November 2021 to November 2022, November 2021 to March 2023, and November 2021 to June 2023.

3.2. Simulation Data Results

The results presented in Figure 6 below compare the outcomes of the traditional sequential adjustment with the proposed improved robust estimation sequential adjustment introduced in this paper. The traditional sequential adjustment exhibits significant biases, while the method proposed in this study demonstrates a better conformity with the actual results. In the linear model, the Root Mean Square Error (RMSE) for the general results is 1.34 mm, whereas the RMSE for the proposed method is 0 mm. In the Weibull model, the RMSE for general results is 2.63 mm, while the RMSE for the proposed method is 1.19 mm. These findings further validate the effectiveness of the method proposed in this paper.

3.3. PUNet Unwrapping Result

As shown in Figure 7, we compared the results of the surface deformation interferogram unwrapping using the PUNet and the minimum cost flow (MCF) method for the same coal mining face. The results indicate that, in interferograms with small deformations, the unwrapping differences between the two methods are minimal. However, in rapidly subsiding regions, severe decorrelation and dense fringes hinder the MCF method from accurately recovering the deformations within complex fringe patterns. At the moment of maximum subsidence velocity, the interferogram exhibits the most severe coherence loss, with the lowest coherence observed at the epicenter. The phase unwrapping results obtained through the MCF method show a significant underestimation phenomenon, while those from the PUNet demonstrate a notable improvement. We computed the time series subsidence above the mining face using the SA method and compared the effectiveness of the two unwrapping methods, as depicted in Figure 8. The MCF method has an obvious underestimation. In contrast, the relative error of the PUNet method is less than 5%. The unwrapping results using the PUNet exhibit a higher accuracy in capturing the deformation patterns. We analyzed the reasons for the significant differences between the two methods and identified several factors. The MCF method relies on the assumption of phase continuity during the unwrapping process. However, in areas of large gradient deformation, the absolute phase difference between adjacent pixels is likely to exceed π, violating the continuity assumption and making it difficult for the MCF method to achieve the ideal unwrapping results. Furthermore, combining coherence maps reveals a poor coherence at the center of the subsidence basin, directly leading to an ineffective observation phase acquisition, which is another factor contributing to the failure of the MCF method. Through analysis of the PUNet, we found two main advantages of this method: (1) the PUNet treats phase unwrapping as a regression problem to learn a direct mapping from the interferograms to the unwrapped phases, allowing the recovery of the unwrapped phases based on contextual information rather than solely relying on the phase difference between adjacent pixels, thus not being constrained by the assumption of phase continuity; and (2) the PUNet incorporates a denoising function that exhibits a strong robustness against noise at different levels, ensuring the quality of the unwrapping.

3.4. Dynamic Subsidence of Mining Area

The obtained results from InSAR represent the line-of-sight displacement, which is transformed into the vertical ground subsidence using the following equation, where $\theta$ denotes the incidence angle [26]:

$$d=los/\cos \theta$$

(14)

As shown in Figure 9a, the dynamic monitoring results above reveal a continual expansion in the subsidence area over time. Figure 9b was obtained by averaging the coherence across all interferometric pairs, revealing a notable consistency between the locations of the subsidence zones and the areas exhibiting a low coherence. To assess the accuracy and robustness of the monitoring outcomes, we compared the precision and robustness of the sequential adjustment and the proposed robust sequential adjustment method at two benchmark points: A and B. As shown in Figure 10, the robust sequential adjustment aligns more closely with the leveling data, exhibiting a reduced data dispersion. A statistical analysis of the accuracy and variance indicates that the RMSE for the sequential adjustment and robust sequential adjustment at point A are 35 mm and 33 mm, respectively. At point B, the RMSE values for the sequential adjustment and robust sequential adjustment are 37 mm and 35 mm, respectively. The relative error is less than 5%, which can meet the requirements of the monitoring of the surface subsidence in mining areas. More intuitively, after the ground surface has reached a state of relative stability, the SA method continues to compute the subsidence fluctuations, which deviates from the actual scenario. In contrast, the RSA method achieves a robustness, and the subsidence values tend toward 0 after the surface has stabilized. The empirical data validate the higher reliability of the proposed RSA method presented in this study.

3.5. Evaluation of the Surface Subsidence Control Effect of Filling Mining

Caving mining is a traditional mining method characterized by the natural collapse of the overlying strata after extraction, resulting in significant surface subsidence values and gradients [27]. In contrast, backfilling mining involves the timely filling of the mined-out voids with specific filling materials during the extraction process to support the overlying strata and to reduce or prevent the collapse of the strata. In order to evaluate the effectiveness of backfilling in the study area, which is a backfilled mining face, we simulated the subsidence under the scenario of the caving mining method, as illustrated in Figure 11. This study employs a probability integration model (PIM) to simulate the surface subsidence caused by caving mining, treating subsidence as the cumulative effect of numerous small unit movements, each of which can be regarded as a random variable [9]. Through probabilistic integration, the effects of these random variables are combined to predict the overall surface subsidence. The maximum subsidence obtained from the simulation of the caving mining method is 1400 mm. The maximum subsidence value obtained by the SBAS-InSAR monitoring is 1150 mm, which accounts for approximately 80% of the subsidence caused by caving mining. It is evident from both the simulation and the measurement results that the control of the surface subsidence in the study area due to backfilling mining is not satisfactory, with a limited magnitude of subsidence control. Furthermore, as depicted in a differential interferogram (Figure 12), the center of the subsidence basin shows decorrelation, indicating significant post-mining surface subsidence gradients occurring in the working face. It should be noted that the presence of two subsidence basins in the figure is due to the existence of another active working face near the studied working face. However, the influence zones of the two working faces overlap only in the peripheral areas.

4. Discussion

In order to evaluate the advantages of the near real-time estimation method proposed in this paper, we discuss the efficiency of the quantization operation from three aspects: the matrix dimension, the number of iterations, and the calculation time. In addition, we also discuss the improvement direction of the robust estimation method in InSAR monitoring.

4.1. Changes in Matrix Dimensions

The previous SBAS-InSAR sequential adjustment method involved incorporating all available prior date data into the computation. The dimensionality of the covariance matrix P1 for the a priori parameters was N × N, where N represents the total number of initial epochs considered. The covariance matrix of the coefficients had a similar N × N dimension. By selectively extracting the matrices related only to the dates overlapping with subsequent times, the dimensionality of the covariance matrix reduced to n × n, where n represents the number of epochs initially computed and the overlapping dates with subsequent updates (n < N). The covariance matrix of the coefficients also reduced to n × n. This reduction in the matrix dimensions facilitates more efficient matrix operations and storage.

4.2. Calculation Time Change

A total of 145 interferometric pairs were processed in the experiment, with a multilook factor of 4:1, resulting in the extraction of 42,798 monitoring points. The CPU utilized for the testing was an Intel Core i9-13900k running at 3.0 GHz. We conducted a statistical analysis of the computational time consumed by various methods during the data processing. As depicted in Figure 13, the iterative estimation process of the robust sequential adjustment method sacrifices some time compared to the general SBAS method, but the overall computation time remains superior. Importantly, our method also possesses the capability to suppress outliers, addressing a critical aspect of the data processing.

4.3. Iteration Number Distribution

To investigate the suitability of the algorithm proposed in this paper, we conducted a statistical analysis of the number of iterations required for each point during the robust estimation process. As shown in Figure 14, the majority of the points exhibited iteration counts equal to or less than two, with only a small subset of points requiring more than two iterations. Therefore, the method proposed in this paper does not increase the computation time too much.

4.4. Improvements in Robust Estimation Methods

The main issues in current research on the sequential adjustment updating of the InSAR data are assuming that the pixels involved from the initial to subsequent updates are the same and that all pixels are persistent scatterers and neglecting the influence of intermittent coherence points on the results. When local changes occur, the pixels that were previously highly coherent may lose coherence, making it impossible to obtain accurate monitoring values using conventional unwrapping methods. To address this problem, two approaches can be considered. First, in scenarios where deformation exhibits spatial correlation, deep learning methods such as the PUNet can be used to obtain unwrapped phases, which can then be used to update the temporal deformation data. Second, in scenarios where deformation lacks spatial correlation, an equivalent weight model considering coherence can be adopted. This model assigns a weight of 0 to unreliable unwrapped phases obtained from low-coherence regions using conventional methods, thus avoiding the introduction of outliers during sequential updating. The equivalent weight model considering coherence can be represented in the form of Equation (15), where, compared to Equation (8), it includes the coherence parameter ${\gamma }_i$. When ${\gamma }_i$ < 0.3, the weight of the observation is set to 0.

$$p_{1i}^{(l+1)}=\left\{\begin{array}{c}{p}_{1i}^l,{\gamma }_i>0.3 \& {|\mathrm{V}}_i|\le k_0\\ p_{1i}^l\frac{k_0}{{|\mathrm{V}}_i|}{(\frac{k_1-{|\mathrm{V}}_i|}{k_1-k_0})}^2,{\gamma }_i>0.3 \& k_0\le {|\mathrm{V}}_i|\le k_1\\ 0,{\gamma }_i\le 0.3{ \mathrm{or} |\mathrm{V}}_i|\ge k_1\end{array}\right.$$

(15)

4.5. The Analysis of Potential Data Leakage and Other Issues with PUNet

The PUNet network was first proposed in [19] without analyzing the potential overfitting caused by the leakage of subsidence peak values. However, the article extensively tested the reliability of the model using real data, indicating that the PUNet network is suitable for scenarios with subsidence rates below 500 cm/a. As long as the subsidence rate does not exceed this threshold, the method can achieve relatively reliable results. Additionally, the training, testing, and validation datasets used by the PUNet are all independently generated, thus minimizing the possibility of data leakage. Due to the direct use of the pre-trained model from reference [19] and the utilization of an unsupervised training paradigm, this paper achieves the cross-domain transfer of knowledge on large gradient deformation phase unwrapping without the need for additional training conditions. As indicated in reference [19], the training data for the PUNet is simulated using Gaussian functions, while the test data in our experiments is derived from real mining area scenes, with no intersection between them, thus eliminating the possibility of data leakage.

The PUNet is a phase unwrapping method aimed at large-scale deformation in mining areas, which has achieved relatively ideal results in current research. Through an analysis of the training sample generation method and network structure of this network, we identify several potential issues. (1) The PUNet uses Gaussian and distortion functions to generate training samples to simulate subsidence basins under different conditions. However, under mining conditions, the PIM is used to describe subsidence caused by mining. It may be more appropriate to use the PIM to simulate samples. (2) According to [19], the PUNet is suitable for scenarios with subsidence rates below 500 cm/a. Surface subsidence monitoring in scenarios with large mining heights and shallow burial depths may not be appropriate. This also indicates that the training data for the PUNet are not comprehensive enough. By considering actual mining scenarios and enriching training samples, it is expected that the application effectiveness of the PUNet can be improved.

5. Conclusions

In this paper, we propose an approach that combines the PUNet and RSA to estimate the near real-time large gradient nonlinear subsidence time series with increasing SAR data. The results of the proposed method are used to evaluate the effect of filling the mining face in controlling the surface subsidence.

(1) The proposed method can overcome the issues of underestimation in traditional InSAR methods and the susceptibility of low-coherence areas to outliers.

(2) The effectiveness of both the PUNet and MCF unwrapping methods is compared, and the results indicate that the PUNet method is more suitable for phase unwrapping in mining areas.

(3) The time series estimation outcomes of the RSA and general SBAS-InSAR are validated using simulations and leveling data, demonstrating that the RSA method achieves a higher accuracy. The relative error is less than 5%, which can meet the requirements of the monitoring of the surface subsidence in mining areas.

(4) The study area of this paper has large gradient subsidence, and the effect of filling mining on controlling the surface subsidence in this study area is not ideal.

The proposed approach provides an efficient means for monitoring large gradient nonlinear subsidence.