Mine subsidence monitoring and prediction integrating SBAS-InSAR technology and BO-Prophet model

Data sources

This study collected ascending synthetic aperature radar (SAR) data from June 23, 2023 to August 4, 2024 and descending SAR data spanning June 21, 2023 to August 2, 2024. The descending data is used to verify the reliability of ascending monitoring results. The ascending dataset comprises 34 scenes of Sentinel-1A IW-mode images with 20 m spatial resolution. Its 12-day revisit cycle ensures temporal consistency for deformation monitoring. Descending data contains 17 scenes of RADARSAT-2 XF-mode images achieving 5 m resolution. The 24-day revisit cycle provides complementary observation angles. In addition, Precise orbit files were collected to eliminate orbital errors through ESA’s Precise Orbit Determination (POD) service. A total of 30 m digital elevation model (DEM) data were obtained to remove topographic phase effects. Table 1 shows the comparison between the two types of data.

SBAS-InSAR

SBAS-InSAR (Small Baseline Subset InSAR) is a time-series InSAR analysis method using multiple master images, Berardino et al. (2002) proposed this technology in 2002. Figure 2 shows the SBAS-InSAR workflow. The core principle involves organizing numerous SAR datasets into interferometric subsets. Each subset contains multiple master images. All interferometric pairs in a subset maintain baseline lengths below critical baseline thresholds. The time baselines are also minimized. This approach resolves decorrelation issues in both temporal and spatial dimensions. The method then applies minimum norm criteria for deformation velocity estimation. Singular value decomposition (SVD) helps obtain time-series deformation of coherent targets. It also generates average deformation rate maps (Li, Xu & Li, 2022). SBAS-InSAR achieves reliable monitoring results with relatively small datasets. This technology proves particularly effective for monitoring surface deformation in mining areas.

BO-Prophet prediction model

This study builds BO-Prophet and BO-LSTM models using PyCharm and Python 3.8, primarily employing the Pandas, NumPy, Kats and Ax libraries. The BO-Prophet model combines Bayesian optimization and the Prophet forecasting framework. Facebook released Prophet as an open-source time series prediction tool in 2017 (Taylor & Letham, 2018).

The Prophet model uses an additive regression framework. It consists of three core components: trend, seasonality, and holiday effects. This structure enables Prophet to handle diverse forecasting problems with unique characteristics (Chang, Yang & Zhang, 2025; Chérif et al., 2023). It effectively captures trend changes, seasonal patterns, holidays, and irregular events in time series data. The basic formula of the Prophet model is:

(1) $$y\left(t\right)=g\left(t\right)+s\left(t\right)+h\left(t\right)+{\epsilon }_t.$$

In the formula, $g\left(t\right)$ represents the trend component, modeling non-periodic changes using piecewise linear or logistic regression; $s\left(t\right)$ denotes the seasonality component, capturing periodic fluctuations like weekly or yearly cycles; $h\left(t\right)$ accounts for holiday effects, addressing impacts from special events; ${\epsilon }_t$ is the error term, and it is used to represent unpredicted random fluctuations.

The logistic growth model employs a sigmoid function to describe growth patterns. The logistic function (sigmoid) is expressed as:

(2) $$g\left(t\right)={\displaystyle \frac{C\left(t\right)}{1+exp\left(-\left(k+a{\left(t\right)}^T\delta \right)\left(t-\left(m+a{\left(t\right)}^T\gamma \right)\right)\right)}.}$$

In the formula, $C\left(t\right)$ indicates the trend’s saturation capacity; $k$ controls the growth rate; $m$ represents the base inflection point parameter of the trend; $a\left(t\right)$ serves as the trend’s changepoint indicator matrix; $\delta$ represents the trend slope adjustment coefficient; $\gamma$ represents the trend inflection point time adjustment coefficient.

Prophet incorporates Fourier series to model periodic variations. This approach represents seasonality through linear combinations of sine and cosine functions:

(3) $$s\left(t\right)=\sum _{n=1}^N\left(a_{n}cos\left({\displaystyle \frac{2\pi nt}{P}}\right)+b_{n}sin\left({\displaystyle \frac{2\pi nt}{P}}\right)\right)=x\left(t\right)\beta .$$

In the formula, $a_n$, $b_n$ represent the Fourier coefficients; $P$ indicates the period length; $x\left(t\right)$ serves as the Fourier basis matrix; $\beta$ denotes parameters following a normal distribution: β ~ N (0, σ²); A larger σ amplifies seasonal effects, while a smaller σ reduces them.

Hyperparameters are parameters set before training machine learning models. Their values cannot be learned automatically during training. Researchers select these values based on experience or experimentation (Bischl et al., 2023). Manual hyperparameter adjustment is time-consuming and expertise-dependent. BO is an algorithm for finding optimal solutions to functions. It efficiently explores parameter spaces. This method identifies strong parameter combinations with minimal trials and time. It is widely used for hyperparameter tuning in machine learning (Snoek, Larochelle & Adams, 2012). BO builds a probabilistic model of the objective function. This model guides the search process. The goal is to find parameter configurations that optimize the objective function. The algorithm handles complex objective functions with multiple peaks or non-convex shapes. It performs effectively in many practical applications (Shahriari et al., 2015).

The BO-Prophet framework applies BO to tune Prophet hyperparameters, as shown in Fig. 3. The process begins by randomly sampling an initial set of hyperparameters. These hyperparameters initialize the model configuration. The framework then evaluates the model training outcomes using predefined metrics. A Gaussian process model integrates historical hyperparameter evaluations and their performance results. An acquisition function analyzes this probabilistic model. It selects the next hyperparameter set with the highest potential performance improvement. The system iteratively repeats hyperparameter selection, model retraining, and performance evaluation. This loop continues until meeting predefined termination criteria (e.g., maximum iterations or convergence thresholds).

LSTM prediction model

This study compares the LSTM neural network prediction model with the Prophet model. LSTM represents a specialized type of recurrent neural network (RNN) (Yu et al., 2019). Compared to conventional RNNs, LSTM demonstrates superior performance in processing time-series data with extended intervals between critical events. The LSTM architecture incorporates memory cells and three gating units. Memory cells store historical information systematically. The input gate regulates the integration of current input data into memory cells. The forget gate manages the elimination of irrelevant information from memory cells (Smagulova & James, 2019). The output gate controls the transmission of stored information to subsequent layers. These gating mechanisms collectively preserve crucial temporal patterns in sequences. They effectively address long-term dependency challenges in time-series analysis. Additionally, this structure resolves gradient-related instability problems common in standard RNNs. Figure 4 illustrates the memory unit configuration of the LSTM neural network.

The memory unit structure of LSTM neural networks can be mathematically represented by the following equations:

(4) $$f_t=\sigma \left(W_f\left[h_{t-1},x_t\right]+b_f\right)$$

(5) $$i_t=\sigma \left(W_i\left[h_{t-1},x_t\right]+b_i\right)$$

(6) $${\tilde{C}}_t=\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(W_C\left[h_{t-1},x_t\right]+b_C\right)$$

(7) $$C_t=f_t\ast C_{t-1}+i_t\ast {\tilde{C}}_t$$

(8) $$o_t=\sigma \left(W_o\left[h_{t-1},x_t\right]+b_o\right)$$

(9) $$h_t=o_t\ast \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(C_t\right).$$

In the formulas, the notation $\left[h_{t-1},x_t\right]$ represents the concatenated matrix of the previous hidden state $h_{t-1}$ and current input $x_t$. The variables $f_t$, $i_t$, and $o_t$ denote the forget gate, input gate, and output gate values at time step t. The operator $\sigma$ indicates the Sigmoid activation function. $W_f$, $W_i$ and $W_o$ are the weight matrices of the forget gate, input gate and output gate, respectively; $W_C$ is the weight matrix for the candidate memory cell ${\tilde{C}}_t$, acting on the concatenated matrix $\left[h_{t-1},x_t\right]$. $b_f$, $b_i$ and $b_o$ are the bias terms of the forget gate, input gate and output gate, respectively; $b_C$ is the bias term for the candidate memory cell ${\tilde{C}}_t$. $C_t$ represents the memory cell state at time step t, which retains long-term information by integrating the previous cell state $C_{t-1}$ and the candidate cell state ${\tilde{C}}_t$.

Mine surface subsidence monitoring and identification

This study processed 34 Sentinel-1A images covering the study area (June 23, 2023–August 4, 2024) using SBAS-InSAR technology. To effectively reduce temporal and spatial decorrelation, the temporal baseline threshold was set to 120 days. The spatial baseline threshold was controlled within 3% of the critical baseline. This configuration generated 134 interferometric pairs (Fig. 5).

During data processing, a coherence threshold of 0.3 was applied. This threshold selected high-coherence points to ensure stability for subsequent analysis. A 30-m resolution DEM of the study area was introduced. This DEM removed topographic phase contributions. Goldstein filtering optimized the interferogram quality. Phase unwrapping was completed using the minimum cost flow (MCF) method.

Atmospheric delay and residual topographic errors were addressed through two SBAS-InSAR inversion iterations. The first inversion separated atmospheric phase from deformation signals. The second inversion further corrected residual errors. This iterative approach effectively enhanced deformation extraction accuracy.

Vertical deformation dominates surface movements in mining areas. This dominance justifies converting line-of-sight (LOS) displacements to vertical deformations. The processing yields annual average vertical deformation rate and cumulative settlement time series for Yinying Coal Mine.

Figure 6 displays the annual deformation rate map of the study area from June 23, 2023 to August 4, 2024. The figure reveals annual deformation rate ranging from −115 to 28 mm/yr. The maximum subsidence rate reaches −115 mm/yr. Four distinct subsidence areas are identified in the study area. Subsidence areas A and B are located within the mining boundaries of Yinying Coal Mine. Subsidence areas C and D occur in adjacent mining areas. Area A covers 1.04 km², and Area B spans 0.64 km². Field investigations confirm both areas are active coal mining zones. The observed subsidence patterns align with on-site geological conditions.

Reliability analysis of monitoring results

To validate the reliability of ascending-orbit SBAS-InSAR monitoring results, this study employed dual validation methods: cross-validation using ascending/descending data and comparison with Global Navigation Satellite System (GNSS) measurements (Tran et al., 2020; Chen et al., 2025). First, 17 descending-orbit RADARSAT-2 images (June 21, 2023 to August 2, 2024) were processed using SBAS-InSAR with consistent parameters (120-day temporal baseline, 3% spatial baseline threshold, 0.3 coherence threshold), generating 36 interferometric pairs with the December 30, 2023 image as the master image. Processing steps followed the same workflow as “Mine Surface Subsidence Monitoring and Identification”, including DEM-based flat-earth phase removal, Goldstein filtering, MCF unwrapping, and two inversion iterations for error correction.

Comparative analysis extracts 3,000 high-coherence points from subsidence areas. Figure 7 demonstrates agreement between ascending (Sentinel-1A) and descending (RADARSAT-2) orbit SBAS-InSAR results. The Pearson correlation coefficient reaches 0.946. This high consistency confirms the high consistency between both datasets.

Furthermore, we collected GNSS monitoring data from four stations. These stations were distributed across the subsidence area in the eastern part of the mining area and residential area west of Region B. We compared GNSS-derived cumulative subsidence time series against SBAS-InSAR results. Figure 8 shows the comparisons: The maximum absolute errors between both monitoring results were 4.684, 3.328, 3.194, and 2.462 mm for points G1, G2, G3, and G4, respectively. GNSS and InSAR deformation trends showed good agreement.

Both validation approaches confirm SBAS-InSAR result reliability. This technique is suitable for mining subsidence monitoring applications.

Characteristics of subsidence basin

To further analyze the deformation characteristics and patterns of the subsidence basin, this study selects Area B as the primary research region. Area B exhibits moderate cumulative subsidence and average subsidence rates. Although Area D shows larger surface subsidence values, mining activities in Area D ceased around March 2024. Within this study’s monitoring period (June 23, 2023 to August 4, 2024), mining only occurred during the first 8 months in Area D. The subsequent 5 months represented a residual subsidence phase. Consequently, Area D provides an insufficiently long effective subsidence time series. Capturing the complete characteristics of surface subsidence during active mining is difficult with such short data.

In contrast, Area B is the core mining area of the Yinying Coal Mine from 2023 to 2025. Area B offers a sufficiently long time series within the monitoring period. Furthermore, the mining advance information for the 150,302 working face within Area B is complete. This complete data enables better comparison with SBAS-InSAR monitoring results. Such comparisons better reflect the relationship between subsidence and mining progress.

Moreover, the subsidence characteristics in Area B closely align with the overall subsidence pattern of the mining region. This alignment effectively represents typical subsidence behavior under standard mining conditions. Selecting Area B for prediction enables a more accurate linkage between mining activity and the subsidence process. It also provides a reliable scenario for validating the prediction model’s practicality. Consequently, results obtained from Area B offer greater reference value for subsidence early warning. Figure 9 illustrates the time-series subsidence plot of this subsidence basin.

The plot demonstrates progressive expansion of subsidence from northeast to southwest. Subsidence area B exhibits continuous downward movement between June 23, 2023 and August 4, 2024. A profile line is established along the mining face orientation. Twenty-five sampling points are selected at 32-m intervals from southwest to northeast along this profile. The SBAS-InSAR monitoring results provide cumulative subsidence time series at these points. These time-series data generate the temporal cumulative subsidence curve along the profile line. Figure 10A displays the time-series cumulative subsidence curve of the profile line derived from this analysis. Furthermore, six characteristic points (P1–P6) were selected at uniform 75-m intervals along the profile line. These points comprehensively cover the subsidence basin’s gradient zone from the center to the margin. P1 and P2 are located at the subsidence center. This area experiences the maximum mining impact, reflecting core extraction zone deformation. The surface cover consists of open woodland. This cover lacks dense root consolidation, making it highly susceptible to subsidence from coal seam extraction. P3, P4, and P5 lie within the subsidence transition zone. Moderate subsidence here illustrates the subsidence funnel expansion process. The surface cover is shrubland. Dense shrub roots provide topsoil consolidation, mitigating coal mining impacts. P6 is positioned at the subsidence margin. Minor subsidence here reflects indirect mining disturbance. The surface cover is farmland. Loose soil structure and severe erosion make this farmland vulnerable to mining-induced subsidence. These six points form a uniform distribution along the working face’s strike-direction profile line. They collectively represent the subsidence center, transition zone, and margin area. This distribution effectively captures Area B’s subsidence patterns.

Figures 9 and 10 reveal the evolution of surface subsidence in area B corresponding to the advance of the 150,302 mining face. The mining face advanced from northeast to southwest starting in June 2023. By July 17, 2023, SBAS-InSAR monitoring detected initial subsidence near the northeastern end of the face. This subsidence formed a small funnel-shaped area with a maximum cumulative subsidence of approximately −20 mm. The surface response was initial at this stage due to the short advance distance. Subsequently, the mining face continued advancing. By October 9, 2023, SBAS-InSAR monitoring showed the subsidence area expanding southwestward, covering the extracted portion of the face. The maximum cumulative subsidence increased to −36 mm, and the subsidence funnel elongated along the face strike. Simultaneously, slight subsidence appeared in surrounding areas affected by indirect mining disturbance.

By March 13, 2024, the mining face had advanced to its central section. SBAS-InSAR monitoring indicated continued subsidence across the affected area. The maximum cumulative subsidence reached −63 mm. However, the spatial extent of the subsidence area showed no significant increase compared to February 6, 2024. From March 13, 2024, to August 2024, the mining face advanced further towards the southwestern end. SBAS-InSAR monitoring then showed the subsidence area fully covering the mining face and a surrounding zone within 300 m. The subsidence center moved approximately 80 m southwestward. The average subsidence value of the subsidence center is −78 mm. Subsidence continued within the existing surrounding disturbance zone, but no significant new subsidence areas emerged.

Furthermore, Fig. 10A shows higher subsidence rates occurred in the northeastern part of the face between June 2023 and February 6, 2024. After the mining advance reached the central section on March 13, 2024, the subsidence rate decreased in the northeastern part but increased in the southwestern part. This differential subsidence rate further demonstrates the direct link between surface movement and mining face activity.

In summary, the direction of subsidence area expansion and the migration path of the subsidence center in Area B, monitored by SBAS-InSAR, synchronized precisely with the mining face advance. SBAS-InSAR also clearly captured the subsidence response in areas indirectly disturbed by mining. These observations confirm the spatial and temporal correlation between mining activity and surface deformation.

Bayesian optimization-based hyperparameter tuning

To optimize hyperparameters for Prophet and LSTM models, this study randomly selects 10 high-coherence subsidence points in the study area. SBAS-InSAR-derived time-series subsidence data are acquired for these points. The first 27 steps form the training set, while the last seven steps serve as the testing set. Mean absolute error (MAE) is selected as the evaluation metric for prediction performance.

The parameter search spaces for Prophet and LSTM are defined as specified in Table 2. The Bayesian optimization framework yields the following optimal configurations: Prophet (seasonality prior scale = 0.008, weekly seasonality = True, changepoint prior scale = 0.69, changepoint range = 0.5, seasonality mode = ‘additive’); LSTM (hidden size = 17, time window = 18, num epochs = 61, num layers = 2).

Model prediction and results analysis

Characteristic point prediction and model performance comparison

This study extracted time-series cumulative subsidence data of six characteristic points in subsidence area B using SBAS-InSAR technology. The first 27 steps of data for each characteristic point served as the training set, and the subsequent seven steps as the validation set. The Prophet model and LSTM model with hyperparameter optimization were employed to predict subsidence for the subsequent 12 steps, The first seven steps of the prediction results were compared with the validation set, while the last five steps were used to demonstrate the future deformation trends. Two accuracy metrics were introduced to evaluate model performance: MAE and root mean square error (RMSE). The formulas are defined as follows:

(10) $$MAE={\displaystyle \frac{1}{n}\sum _{i=1}^n\left|y_i-{\hat{y}}_i\right|}$$

(11) $$RMSE=\sqrt{{\displaystyle \frac{1}{n}\sum _{i=1}^n{(y_i-{\hat{y}}_i)}^2}}.$$

In these formulas, $y_i$ represents measured values, ${\hat{y}}_i$ denotes predicted values. Using MAE and RMSE, this study systematically compared the Prophet and LSTM models. The analysis revealed performance differences in subsidence prediction, providing a scientific basis for model selection.

Figure 11 illustrates the comparative performance of the prediction models. The BO-Prophet model exhibits specific error distribution characteristics across seven prediction periods for six characteristic points. The maximum absolute error values for each characteristic point are 4.35, 5.38, 4.45, 3.72, 2.33, and 2.63 mm. Notably, the maximum absolute error for all characteristic points (except P4) concentrates in the third prediction step. This trend aligns with significant subsidence observed in Area B from June 17 to June 29, 2024.

The minimum absolute errors of the BO-Prophet model at six points are 0.59, 0.09, 0.79, 0.15, 0.27, and 0.48 mm. These values confirm high consistency between the BO-Prophet predictions and SBAS-InSAR monitoring data. Visual comparisons demonstrate strong agreement between predicted results and monitoring data.

According to the error indicators presented in Table 3, the BO-Prophet model achieves MAE values of 2.18, 2.69, 2.46, 1.90, 1.25, and 1.63 mm across the six characteristic points. The RMSE values are 2.52, 3.26, 2.85, 2.31, 1.45, and 1.79 mm. Compared to the BO-LSTM model, the BO-Prophet model shows lower errors in both MAE and RMSE. These results validate the higher accuracy and stability of the BO-Prophet model in predicting surface subsidence.

Table 3:

Comparison table of the accuracy of each prediction model.

Point	BO-LSTM		BO-Prophet
	MAE/mm	RMSE/mm	MAE/mm	RMSE/mm
P1	3.08	4.16	2.18	2.52
P2	2.52	3.65	2.69	3.26
P3	2.96	3.81	2.46	2.85
P4	2.59	3.44	1.90	2.31
P5	2.34	3.52	1.25	1.45
P6	2.69	3.41	1.63	1.79

DOI: 10.7717/peerj-cs.3327/table-3

To comprehensively assess prediction accuracy, the BO-Prophet model generated subsidence forecasts across area B. These predictions were compared with SBAS-InSAR measurements through correlation analysis. Figure 12 demonstrates strong consistency between predicted values and SBAS-InSAR values across multiple time spans (12d, 24d, 48d, 72d, 84d). All correlation coefficients exceed 0.96, confirming the model’s effectiveness in large-scale subsidence prediction. Further analysis examines correlation coefficient variations under different time spans. As shown in Fig. 12F, the coefficients gradually decrease over time. This trend indicates diminishing prediction accuracy with increasing forecast duration.

Distribution of errors

The study investigates spatial error patterns by overlaying BO-Prophet predictions with SBAS-InSAR measurements in area B. This process absolute error distribution map. Figure 13 reveals two critical trends. First, prediction errors progressively increase over time. Second, errors predominantly cluster in the northwestern sector of area B. This area corresponds to the marginal slow subsidence zone of the depression basin. These findings demonstrate BO-Prophet’s superior short-term forecasting capability. However, error accumulation becomes evident with extended prediction horizons. The spatial error concentration suggests systematic biases in specific geological subregions.