Deformation monitoring and prediction for residential areas in the Panji mining area based on an InSAR time series analysis and the GM-SVR model

2.1 Study area and data preparation

In this work, the village of Yang juzhuang in the Huainan mining area was selected as a case study (Figure 1). The Huainan coalfield is located in the central Anhui Province, East China,which has a mining history that covers the past 100 years. It is one of the regions that has been noted in Chinese history for the early exploitation and use of coal resources. The rich coal resources in this region make up for a shortage of energy in the eastern part of China. The Huainan coalfield plays an irreplaceable role in ensuring the sustainable and healthy development of the national economy [39]. Mining activities in the Huainan area began in the 1920s. However, due to continuous improvement in the development of coal resources in mining areas, the geography of the coal mining area has also increased sharply. As a result, the spatial extent of surface subsidence in the Huainan coalfield was only 90 km² in 2004, which grew to 150 km² in 2010 and is still increasing [40]. The urban and rural environments of the mining area are deteriorating, and the lives of the people in the mining areas have been greatly affected. All of these problems have created security issues, which have attracted the attention of the relevant local authorities.

Figure 1

A regional overview of the village of Yang Juzhuang, Huainan, China: (a, b, c) location of the study area, (d) the waterlogged area, and (e) the damaged road

Sentinel-1 is the first of the Copernicus Programme satellite constellation studies conducted by the European Space Agency. Following ERS and ENVISAT it is another C-band remote sensing satellite of great importance that collects data in every weather condition, day or night. The first satellite, Sentinel-1A, launched on April 3rd, 2014. Sentinel-1 provided a data continuity bridge from the ERS and Envisat missions, with further enhancements in terms of revisits, coverage, timeliness and reliability of service. There are wide ranges of applications for the data, which include land monitoring (agriculture, forestry) [41, 42], marine monitoring (sea-ice levels, ocean oil spills, ship activity) [43], and emergency responses (flooding, landslide and volcanic) [44]. In contrast to the traditional Strip Mape (SM) mode and ScanSAR mode, the Sentinel-1 data adopts a new Interferometric Wide Swath (IW) mode of terrain observation by employing a progressive scans (TOPS) technology. In this paper, Yang Juzhuang was selected as the study area (Figure 1). A total of 13 C-band Sentinel-1 IW images from December 27, 2016 to May 20, 2017 were collected (polarization: VV, revisit period: 12 days, orbit configuration: ascending). The acquisition dates and other relevant parameters are shown in Table 1. The corresponding Sentinel-1 Precise Orbit Ephemerides were collected.

2.2 Rationale of the SBAS-InSAR technique

An InSAR time series analysis was developed on the basis of conventional InSAR technology. Estimating the parameters of the phase statistical characteristics of the multi-temporal SAR data can effectively reduce the errors introduced

by factors such as atmospheric delays and low coherence in the interferogram. The basic principle of SBAS, which is one of the branches of the time-series InSAR technology, is to use the observation results of the conventional D-InSAR monitoring as a single observation value, and then obtain a high-precision deformation time series according to the least squares method. A detailed discussion of this approach and its implementation can be found in Lanari et al. (2004) [45]; blow is a brief summary.

Considering a set of N + 1SAR images acquired at the ordered times (t₀, t₁..., t_n) over an area of interest, then M interferograms will be generated with small temporal and spatial baselines to minimize the decorrelation phenomena. In addition, M meets the following criteria [46]:

assuming that the interferogramj generated from the time t_A, t_B (t_B > t_A). After removing the effects of the flat effect

(the reference ellipsoid) and the topographic phase, the interferometric phase of any pixel in the interferogram j can be expressed as:

where φ (t_B) and φ (t_A) are the phase values of the SAR images at time tB and tA, respectively; φ_def,j is the deformation phase along the line of sight (LOS) direction; φ_topo,j is the phase error due to the inaccuracy of the reference DEM; φ_atm,j is the atmospheric delay phase; and φ_noise,j is the noise effects phase.

The time series deformation, calculated with the SBAS-InSAR algorithm can be represented as follows:

where B is an M × N matrix and δφ is a vector which represents the interferometric phase values. When the coefficient matrix B is a full rank (M≥ N), the deformation rate can be obtained by the least squares method, and when M<N, the matrix B has a rank deficit. The singular value decomposition (SVD) is used to obtain the deformation rate. Finally, based on the time interval of the image, we can determine the deformation of the corresponding time period.

The SBAS-InSAR process includes the following key steps:

1. The master image was selected by comprehensively considering a temporal baseline, a spatial baseline and a Doppler centroid frequency baseline, and the interferogram connection diagram was generated according to the small baseline set principle.

2. For each small baseline pair, we do differential interferometry, and select the appropriate algorithm to unwrap the interferogram.

3. We artificially remove some interferometric pairs with very poor interferometry quality and select ground control points for orbit refining and reflattening.

4. By establishing a model on a pixel point with good coherence, the deformation and deformation rate of the target points are extracted by the SVD. The detailed processing steps are shown in Figure 2.

Figure 2

Processing flowchart of mine deformation monitoring and prediction in mining area based on SBAS-InSAR technology and GM-SVR model

2.3 Establishment of GM-SVR prediction model

Due to the destruction of in situ stress distributions in the mining area, the overall structure of the surrounding rock and overlying strata is affected, and under the influence of some external factors, nonlinear deformations of the ground is caused during the mining process. Since the deformation process contains many unknown factors, the deformation process is very complicated. Therefore, the traditional single linear or single nonlinear forecasting model is used for forecasting, which is not very practical and cannot achieve the desired prediction effect. In this paper, based on the GM (1, 1) model and the optimized Support Vector Machine Regression model, a combined prediction model GM-SVR is established. This model makes use of the linear prediction advantages of the grey model and the excellent generalization abilities and learning abilities of the nonlinear data of the support vector machine model to perform high-performance predictions of the ground subsidence in the mining area.

Generally, a system with poor and incomplete information is called a grey system. Before establishing the grey prediction model, all original time series need to be regarded as having grey amounts with grey features. The randomness of the initial data are weakened by the pre-preprocessing of chaotic data into ordered data. Then, differential equations are established to predict future development trends. The advantage of the model is that for systems with small data samples or inaccurate information, the orderly mathematical rules in the internal continuous change process can be excavated, and the model has a wide range of adaptations. Detailed principles can be found in Mao et al. 2006 [47]. This article only gives a brief introduction.

By accumulating the original data, a new data sequence is generated. The fluctuations and randomness of the new data series are weakened [48]:

where X⁽¹⁾ are the 1-AGO (accumulated generating operator) of the initial data. Therefore, the GM (1, 1) model can be established as Eq. (5) (a differential equation):

Using the least squares method, the values of the parameters α and ucan be obtained. Therefore, X^(1)(t+1)can be obtained by using Eq. (6).

whereX^(1)(t+1)is the predicted data. Finally, the last predicted data can be obtained byX^(1)(t+1)minusX^(1)(t).

Support vector machine (SVM) is a new method in data mining. This method is based on a general learning method and a statistical learning theory that specializes in the study of small samples. Support vector machine maps nonlinear problems in low-dimensional space to high-dimensional space and solves nonlinear problems in low-dimensional space through decision functions in high-dimensional space. Additionally, the support vector machine has a corresponding penalty mechanism to prevent the occurrence of overfitting. SVM not only considers the requirements of the approximation accuracy but also considers the complexity of the approximation function and can obtain optimal result under limited conditions. Only a brief summary is presented here. The SVM objective function can be expressed as [49, 50]:

where ω denotes the direction vector, C denotes the adjustment factor, and ξ _i and ξi⋆are slack variables. The deviation of the non-insensitive (ϵ) out-of-band training sample points is measured by introducing slack variables. Then, a Lagrange multiplier α was introduced, and samples were mapped to high-dimensional Hilbert space according to KKT (Karush–Kuhn-Tucker constraint). Finally, a linear regression is applied to the transformed results in the high-dimensional space to obtain a nonlinear regression function.

In this paper, the initial fitting of samples is carried out through the GM (1, 1) model to grasp the original internal development trend of the samples. Then, the residual value sequence is obtained by comparing the predicted value of the grey model with the initial sample. Normally, the deformation process of the surface of the mining area is very complicated, and the fluctuations in the data that are obtained by monitoring are relatively large. Therefore, the nonlinearity of the residual sequence predicted by the grey model is usually high. Based on this, an optimized SVM regression model is used to establish the residual correction model and further explore the internal rules of the residual series. The support vector machine regression model is optimized by searching the optimal penalty parameter c and the variance G in kernel function by the grid search method. Finally, by training the sample in the model, the predicted value of the grey model and the model residual correction of the SVR are obtained. The sum of the predicted value and the residual error correction value is the final predicted value of the combined model. The technical process of the deformation prediction combined with SBAS-InSAR is shown in Figure 2.

3.1 Data processing and analysis of deformation monitored by the SBAS-InSAR

In this paper, the SBAS-InSAR method was used to process Sentinel-1 data. This section describes the main parameters and processing steps and provides details of the analysis of the deformation of the study area.

First, an interference connection graph was generated with a 120-m perpendicular baseline threshold and a 60-day temporal baseline. Sixty-three interferometric pairs were obtained, and the relative relationship between the time and space position of the connection graph is shown in Figure 3. Next, according to the connection relationship of the connection diagrams, a differential interferometry set is generated, precision orbit data is used to remove the flat effect, and the external Shuttle Radar Topography Mission digital elevation model SRTM3 DEM data are used to eliminate topography phase [51]. The phase unwrapping of the small baseline interferogram was carried out by the minimum cost flow method. Subsequently, based on the stability of the amplitude and phase, the stability point in the study area was selected. Then, a third-order polynomial model was used for orbit refinement and to remove the trend phase caused by the orbital error. Finally, deformation rate inversion and geocoding were performed, the atmospheric delay phase of the coherent point is separated according to the temporal and spatial characteristics of the atmospheric phase, and the line-of-sight deformation is segmented according to the cosine of the incident angle to obtain the ground vertical displacement.

Figure 3

Temporal and spatial baseline connection diagram of interferogram pairs

Figure 4 shows a time series diagram of deformations with typical representative characteristics in the study area from December 27, 2016 to May 20, 2017. All images show the cumulative amount of deformation using December 7, 2016 as the reference time. The upper right corner of each figure shows the acquisition time of the image. Figure 5 shows the deformation rate of the study area. From Figure 4 and Figure 5, the maximum cumulative vertical subsidence caused by mining in the entire residential area during the monitoring period is 48 mm; the average cumulative displacement is −17.8 mm; the maximum deformation rate reaches -90 mm/year.

Figure 4

Deformation sequence of the study area (all maps refer to December 27, 2016): (a) January 8, 2017, (b) February 1, 2017, (c) March 9, 2017, (d) April 2, 2017, (e) April 26, 2017, and (f) May 20, 2017

Figure 5

Cumulative deformation rate of the study area and location of GPS station: (a), (b) GPS station, and (c) cumulative deformation rate

The accuracy of SBAS in monitoring deformation that is caused by mining is verified by using the double-differenced static positioning measurement of BDS/GPS combined system (hereinafter referred to as GPS). The accuracy of the elevation direction of the 1 h solution of the GPS system is better than 4 mm, which meets the requirements of the deformation monitoring accuracy of this residential area. The spatial distribution of the site is shown in Figure 5. Figure 6 shows the comparison of the observation results of the deformation sequence at the observation points between the SBAS-InSAR technology and the GPS technology. Analysis of Figure 6 reveals that the maximum absolute error, relative error and average absolute error of vertical displacement between the measurement results of SBAS-InSAR and GPS at the G1 point are 4.33 mm, 2.38mm and 2.01 mm respectively, and at the G2 point, they are 2.87 mm, 1.61 mm and 1.34 mm, respectively.

Figure 6

Comparison of the vertical deformation time series of two observation points obtained by SBAS-InSAR and GPS

3.2 Assessment of prediction model

The time series data of the observation points (spatial distribution shown in Figure 5) that were obtained by the SBAS-InSAR technology were input to the GM-SVR model for fitting and prediction. To verify the performance of the proposed GM-SVR model for deformation predictions in the mining area, two other schemes are established in this paper for comparison: scheme one, using a single GM(1, 1) model for prediction; scheme two, using a single Support Vector Machine Regression model for prediction. From December 27, 2016 to April 02, 2017, eight settlement data with a time interval of 12 days were used as training samples for modelling, and four settlement data from April 02, 2017 to May 20, 2017 were used as a test sample for the prediction

model. The detailed data and results are shown in Tables 2-5 and Figure 6.

Figure 7 shows that the deformation sequences of the two observation points are similar during the observation period and demonstrate an overall growth trend. By comparison, in Figure 7, the GM (1, 1) model, SVR model and GM-SVR model proposed in this paper have better fitting effects at the two observation points. However, in the comparison of the forecasting results, the GM (1, 1) model shows unsatisfactory prediction results and large deviations, which is obviously lower than the accuracy of the SVR model and the model proposed in this paper.

Figure 7

Comparison of the original values of the observation points and the results of the three models

Table 2 shows that during the observation period of the G1 point, the fitting effect of the three models are good and the result is relatively stable. Among the three models, the model proposed in this paper has the best fitting effect. The relative error is generally maintained below 2%. Only the data of No.3 and No.5 exceed 2%. Table 3 shows the prediction results of the three models based on the training samples at the G1 point. It can be seen that the deviations between the predicted and the measured values of the grey model are relatively large; the SVR model and the GM-SVR model of this paper have relatively good results which are relatively stable and accurate. Table 4 shows the fitting results of the three models at G2. Differing from the result at G1, the fitting effect of the GM (1, 1) model and the SVR model is relatively poor, and the accuracy and stability of the fitting result are all reduced. The maximum relative errors of the fitting results of the GM (1, 1) model and the SVR model reach 52.8% and −32.78%, respectively. However, the results of the GM-SVR model still maintain high degrees of accuracy and stability. Table 5 shows the settlement prediction results of the three models at G2. Similar to the prediction result at G1, the accuracy of the GM (1, 1) model at G2 is still not high. In contrast to the prediction result at G1 point, the prediction effect of the SVR model at G2 is not ideal, and the maximum relative error even reaches −35.5%. Only the prediction results of the GM-SVR model are similar to those predicted at G1, and still maintain high accuracy and stability. The relative error of the GM-SVR results is generally maintained below 10%. Based on this, for the deformation sequence in this study area, the GM (1, 1)

model prediction result is not ideal and has essentially deviated from the monitoring value; the prediction result of a single SVR model is not sufficiently stable, and only the GM-SVR model result shows high accuracy and stability.

To test a model’s accuracy, there are generally three methods: residual size, correlation degree, and posterior error. To evaluate the performance of the model more comprehensively and intuitively, this paper proposes to evaluate the accuracy of the model using two indicators, root mean square error (RMSE) and mean absolute error (MAE). The results are shown in Table 6.