A deep learning-based combination method of spatio-temporal prediction for regional mining surface subsidence

Table of Contents

Study area

A mining area in Heze City, Shandong Province, China was selected as the study area, its geographical location and scope were shown in Fig. 1. This area belongs to the Yellow River alluvial plain, with a flat terrain, a ground elevation of + 41 m to + 46 m, a natural terrain slope of 0.2‰, and a total area of approximately 259 km². The main coal types are fat coal, one-third coking coal, and gas coal. The depth, thickness and date of coal mining is 938 m, 6.8 m and since 2006, respectively. The overlying features include mostly farmland with common crops such as wheat, corn, and sweet potato, along with economic crops such as soybean, cotton, and vegetables. Furthermore, the industrial, construction, animal husbandry, and tertiary industries are relatively developed. The above natural conditions and the coal mining have caused serious surface subsidence problems in this study area, addressing the contradiction between underground coal mining and surface protection is an important work to achieve efficient and civilized production, and it is particularly important to conduct temporal monitoring and prediction.

In addition, the leveling monitoring surface subsidence values of 211 leveling points were compiled to compare, analyse, and verify the accuracy and reliability of the SBAS InSAR monitoring surface subsidence time series. The distribution and extent of the 211 leveling points were shown in Fig. 1c. The leveling data was measured according to the fourth-order leveling survey rules and using the American Trimble DiNi03 electronic level, and the allowable closing error was $20\sqrt L$ mm ($L$ is the length of the leveling line in km).

Accuracy analysis of SBAS InSAR

The Sentinel-1 satellite constellation is an earth observation plan of the European Space Agency Copernicus Program (Global Monitoring for Environment and Security) and is composed of Sentinel-1 A satellite and Sentinel-1 B satellite. The two satellites were successfully launched on April 3, 2014, and April 25, 2016, respectively. The two satellites can observe ground with all-weather and continuous radar imaging since they carry a C-band synthetic aperture radar and fly at an on-orbit altitude of 693 km. More importantly, Sentinel-1 A and B can work in different modes, e.g., single and dual polarization work modes, with excellent timeliness and reliability³¹. In this study, 66 Sentinel-1A SAR images from January 10, 2019, to April 5, 2021 were selected and processed using SBAS InSAR to obtain the mining surface subsidence time series of the study area. In detail, due to high accuracy and less time cost, the GAMMA software was used to calculate surface subsidence based on SBAS, with temporal baseline of 365 days, max normal baseline of 45% and unwrapping coherence threshold of 0.3 for first unwrapping and 0.2 for the second. In order to obtain a uniformly distributed time series of mining surface subsidence and minimize the surface subsidence prediction error, linear interpolation was performed when encountering subsidence values with a 24 day interval (January 10, 2019 to February 3, 2019, February 3, 2019 to February 27, 2019, August 26, 2019 to September 19, 2019). The resulting time series comprised subsidence values over 69 periods, each period as 12-days. Figure 2 shows the surface subsidence time series of six periods.

The leveling monitoring covered 19 periods from January 26, 2019, to April 4, 2021, with data from 211 leveling points. The SBAS InSAR monitoring spanned 69 periods from January 10, 2019, to April 5, 2021. To align the data, a piecewise linear interpolation was applied to the leveling results to match the SBAS InSAR data. The difference between the two monitoring methods, i.e., leveling monitoring and SBAS InSAR was then calculated for analysis. Figures 4 and 5 compare surface subsidence results from leveling and SBAS InSAR. Figure 4 displays surface subsidence curves for points H2, Q56, and Q10 over 69 dates from January 10, 2019, to April 5, 2021. Figure 5 shows histogram plots of surface subsidence values at 211 leveling points across three monitoring periods.

Figures 3 and 4 shows that leveling and SBAS InSAR could successfully monitor the continuous surface subsidence during the period from January 10, 2019, to April 5, 2021, and the surface subsidence time series in mining area exhibited obvious nonlinear characteristics.

For example, on leveling point H2 in Fig. 3a, during the first 16 monitoring periods, the subsidence values were similar and highly correlated. From the 17th to 29th monitoring periods, subsidence continued but at different rates between leveling and SBAS InSAR. During the 30th to 41st periods, leveling showed a rapid subsidence with an increased velocity, while SBAS InSAR also showed rapid subsidence, but with a lower velocity. From the 41st to 69th periods, leveling and SBAS InSAR both showed slowing down of surface subsidence, with SBAS InSAR having higher velocity than leveling. In Fig. 4, the histograms of mining surface subsidence values monitored by leveling and SBAS InSAR on 211 leveling points exhibited similar shapes. This suggested that the mining surface subsidence values of leveling and SBAS InSAR monitoring were relatively close to each other.

Overall, the analysis and comparison of the mining surface subsidence of SBAS InSAR and leveling monitoring demonstrated the consistency and accuracy of the SBAS InSAR monitoring, the spatio-temporal distribution of surface subsidence in mining areas monitored by SBAS InSAR was reliable.

Spatio-temporal distance and K-means clustering

This study used a combination strategy of spatio-temporal distance and K-means clustering algorithm (hereinafter referred to as SDK) to determine the spatial partitions of the mining surface subsidence time series. This method comprehensively considers and describes the spatio-temporal correlation and similarity of surface subsidence time series between adjacent pixels, and obtains more reasonable partitioning results.

Firstly, the spatio-temporal distance between the nearest neighboring pixels is calculated. The spatio-temporal distance refers to the weighted sum of temporal and spatial distances between the two pixels.

Assuming two pixels are $P_i$ and $P_j$, with longitude and latitude coordinates of $\left( lon_i ,lat_i \right)$ and $\left( lon_j ,lat_j \right)$, respectively, and the average values of the surface subsidence time series monitored by SBAS InSAR are $t_i$ and $t_j$, the spatio-temporal distance $d_i,j$ between them can be calculated using the following equation:

$$ d_i,j = \omega_s \sqrt {\left( \fraclon_i – lon_j \alpha \right)^2 + \left( \fraclat_i – lat_j \beta \right)^2 } + \omega_t \gamma \left| t_i – t_j \right|, $$

(1)

where $\omega_s$ and $\omega_t$ refer to the weights of temporal and spatial distances, $\alpha$ and $\beta$ refer to the coefficients for scaling latitude and longitude, $\gamma$ refers to the coefficient for scaling temporal distance. In this study, $\omega_s$ and $\omega_t$ were set to 0. 5, $\alpha$, $\beta$, and $\gamma$ were set to 1.

Then, K-means clustering is performed based on the above spatio-temporal distance. K-means clustering is widely recognized as an effective clustering method^32,33. The algorithm first randomly selects k pixels as cluster centers, where k represents the desired number of clusters and can be determined by the following elbow method^34,35,36. Subsequently, calculates the spatio-temporal distance between each pixel and each initial cluster center. Based on the shortest distance, assigns each pixel to the nearest cluster center to form initial clusters. Recalculates the cluster center of each initial cluster based on the existing pixels in that cluster, and determines the new cluster centers. Subsequently, calculates the spatio-temporal distance between each pixel and each new cluster center. Based on the shortest distance again, assigns each pixel to the nearest cluster center to form new clusters.This iterative process will continue until one of the following termination conditions is met: (1) No (or the minimum number) pixels are reassigned to different clusters, (2) no (or the minimum number) cluster centers change again, and (3) the sum of squared errors (SSE) is minimized. In this study, the third condition was selected as the iterative termination condition, and the objective function of K-means clustering was defined as³⁷:

$$ \min \left( SSE \right) = \min \left( {\sum\limits_i = 1^k \sum\limits_x \in C_i \left( D_i – x \right)^2 } \right), $$

(2)

where $C_i$ refers to the ith cluster, $D_i$ refers to the cluster center of $C_i$, $x$ refers to the pixels of $C_i$.

LSTM and GRU

LSTM and GRU are two types of recurrent neural networks^38,39. The main difference between the two lies in the different gating mechanisms. LSTM uses three gates to control the information flow, namely input gate, forget gate, and output gate. The input gate and forget gate respectively control whether new input data and previous memory are written, and the output gate controls whether the output values should be passed to the next layer⁴⁰. GRU only uses two gates, namely reset gate and update gate. The update gate controls the degree that the previous state information is retained in the current state, while the reset gate is used to determine how the current state is combined with previous memory. Compared to LSTM, GRU is simpler, fewer gating structures reduce the risk of over fitting, fewer parameters reduce the computational complexity, and improve operational efficiency. Overall, the advantages of the GRU lie in its simplicity, faster training speed and computational efficiency, and higher generalization ability. Figure 5 illustrates the basic structures of LSTM and GRU neural networks.

Snake optimization algorithm

The number of neurons, learning rate, dropout rate, and batch size of training samples are key parameters that affect the performance of the GRU network model. A novel meta-heuristic optimization algorithm called snake optimization algorithm was used to globally optimize these parameters⁴¹. It simulates the feeding and breeding behaviors of a snake to reduce the average prediction error and achieve efficient parameter combination optimization, with numerous advantages, e.g., faster compute, higher precious and robustness^42,43. It is widely used in the fields of machine learning and deep learning. its process involves the following steps⁴⁴:

(1)

Parameter definition and population initialization: Determine the parameters that need to be optimized. In this study, the parameters were the number of neurons, learning rate, dropout rate, and batch_size of GRU. Additionally, apply the SO algorithm to generate an initial set of positions (parameter combinations) for the GRU model, with each position corresponding to an individual.
(2)

Fitness calculation: Calculate the RMSE of the model’s predicted subsidence values. Lower RMSE values indicate better fitness. The fitness function can be mathematically expressed as follows^45,46,47:

$$ \textRMSE = \sqrt \frac1n\sum\limits_i = 1^n \left( \haty_i – y_i \right)^2 , $$

(3)

where n denotes the number of samples, $y_i$ denotes the actual subsidence values, and $\haty_i$ denotes the predicted values.
(3)

Iterative optimization and model training: Use SO algorithm to simulate the feeding and breeding behaviors of a snake and adjust its position to find the parameters combination with the best fitness value. Obtain the optimal combination of network model parameters, and use these optimal parameters to train the GRU prediction model.

Deep learning-based combination method of spatio-temporal prediction

To address the issues of existing models not taking into account the spatio-temporal correlation this study proposed spatio-temporal prediction method (Fig. 6) for regional mining surface subsidence can adaptively determine the optimal parameters, avoiding the tedious and random manual parameter adjustment, and ensuring that the model has strong adaptability and higher accuracy. The implementation of this method mainly involves the following three steps: Firstly, the SKD method is used to divide the surface subsidence time series into a group of partitions. Then, learn different subsidence patterns and construct local models within each partition. Finally, use the well-trained model to make short-term prediction of future regional mining surface subsidence.