2.1. Probability integral method

The probability integral method (PIM) is based on random medium theory and is widely used in mining subsidence prediction. It applies granular medium mechanics in discontinuous medium mechanics to study strata and surface movement. From a macro point of view, the law of strata and surface movement caused by mining is similar to that described by the granular medium model as a random medium (^{Chen et al., 2018}). The formula of surface subsidence at any point W(x, y) caused by mining is as follows:

Where q is subsidence coefficient, r is the main influencing radius, Y = H ₀ / tanb, H ₀ is the average mining depth, tanb is the tangent value of the main influence angle, (xi _, y _i ) planar coordinates for mining unit i, (x, y) is the coordinate of any surface point.

2.2. Power exponent Knothe function model

Knothe time function model was proposed by Knothe, a Polish scholar, in 1952. The model is used to describe the dynamic process of surface subsidence. During the mining of the working panel, the subsidence value is assumed to be 0. The power exponential Knothe time function model can be expressed as follows:

Where w(t) is the subsidence value at a certain time, c is the time influence coefficient related to mining conditions of the working panel and characteristics of overlying strata; k is the fitting parameter of the model.

The power exponential Knothe function has three parameters to be estimated, namely, the maximum subsidence value wm, and the constant parameters c and k. c and k mainly affect the shape of the function image (^{Yang et al., 2016}). Suppose w _m= 1000mm, as shown in Fig 1(a), when c = 0.03, k values are taken as 5, 10, 50, 100, 200 and 400, respectively. It can be seen that the larger the parameter k is, the longer the initial subsidence time is, and the longer the time to reach the maximum subsidence value is, but the shape of the function image is roughly unchanged. As shown in Fig 1(b), when k = 70, c values are 0.01, 0.03, 0.05, 0.1, 0.2 and 0.4 respectively. It can be seen that the smaller the parameter c is, the faster the surface subsidence speed is, but the change of c value has a greater impact on the shape of the function image (^{Diao et al., 2016}).

Figure 1 Effects of parameters k and c on the shape of the Knothe function 

2.3. Construction of InSAR-PEK Model

a. Subsidence values of InSAR data

^{Chen et al. (2018}) achieved large-gradient deformation monitoring by using lidar to repair the hole field of SAR data. They proposed that when the coherence coefficient was lower than 0.32, the maximum detectable deformation gradient of D-InSAR was close to 0. To obtain the interferogram with high coherence coefficient, the coherence threshold of phase unwrapping was set to 0.45 in SAR data processing of the working panel, and the interferogram with a "cavity field" was formed.

The double track method was used to obtain the LOS deformation values of the sentinel data, with the assumption that the relative subsidence value between the two scenes of the InSAR Data was ∆W _los = [W _{los 2-1} , W _{los 3-2} , W _{los 4-3} ^, L, W _{los 42-41} ], and the vertical subsidence value of the mining area extracted by the LOS three-dimensional decomposition method was W _sar =[W _{sar 2-1}, W _{sar 3-2} , W _{sar 4-3 ,} L_, W _{sar 42-41}]. Among them, the SAR data in the main subsidence period were removed due to the influence of large-gradient strain.

b. Fitting parameters

1) PIM predicted the maximum subsidence value of the working panel according to geological and mining conditions. Assuming that the initial subsidence value was 0 mm, the initial subsidence values were accumulated to obtain the data of the initial period;

2) PIM predicted the final subsidence value, and the data in residual period were obtained by accumulating the final SAR data;

3) The power exponent Knothe time function model was used for parameter inversion:

• (1) The first group of parameters c _{level_wz} and k _{level_wz} were obtained from the complete leveling observation values;

• (2) The parameters c _{level_bf} and k _{level_bf} were obtained from the leveling data excluding the data in the main settlement period;

• (3) The parameters c _sar and k _sar of SAR Data were obtained.

4) The RMSEs were calculated by comparing the predicted subsidence values of the main subsidence period and the subsidence values of the whole observation period with the leveling data.

c. Accuracy of parameters

Based on the calculated parameters, the subsidence values of each period were estimated and compared with the measured values at the same time, and the fitting accuracy of the three parameters was analyzed. The specific method flow is shown in Figure 2.

Figure 2 Flow chart of the experiment 

4.1. Introduction to the experimental area and SAR data

The experimental area is located in a coal mining panel in Huainan city, Anhui province. The mining of the working panel began on Mar. 26, 2017, and finished on Dec. 25, 2017. The hydraulic support was shipped out by the end of Jan. 2018. Actual stoping length of the working panel was about 1528 m (including 8m width of open cut) and width 251m (including transport and track paths).The mining lasted for 275 days, the average mining speed was 5.56m/d, and the coal output was 1.428 million tons. The thickness of the working panel was 1.0 ~ 4.2m, with an average of 2.8 m and an average mining height of 2.9m. The dip angle was 0 ~ 6°with an average of 3°, belonging to a nearly horizontal coal seam. The buried depth of the coal seam was 618 ~ 723 m, with an average of 668 m. The floor elevation of the coal seam was -564.1 -712.8m. The Cenozoic loose bed was about 420 m thick and belonged to the super thick loose bed. The surface above the working panel was flat with farmlands, roads and ditches, and the elevation was between 18.0m and 23.0m, with an average of 21.5m. The surface movement lasted for 597 days. The initial period of surface movement was short, about 7 days, the active period was 123 days, and the recession period was 467 days. The subsidence in the active stage accounted for 92.2% of the total subsidence.

The image data were Sentinel-1A data, C-band fully-polarized SAR image, and the data modes were SW, IW and EW. The experiment adopted IW polarization mode, with a resolution of 20 m and a standard width of 250 km. The revisit period was 12 days, and the revisit period among some images was 24 days. The subsidence values in initial period and residual settlement period were obtained from 41-scene Sentinel-1A data from May 8, 2017 to Jan. 10, 2019.

4.2. Subsidence prediction by PIM

Due to the thick loose layer in the working panel, PIM was used to predict the subsidence coefficient q =1.1, in which the horizontal movement coefficient b =0.32, which affected the propagation angle q = 86°. The tangent of the main influence angle tg = 1.70. The expected subsidence is shown in Fig. 4. The maximum subsidence value was 1882.46 mm. The black rectangular box is the working panel range, and the black line is the mining progress and the corresponding time. The mining of the working panel, was stopped on Mar. 26, 2017, while the first observation was made on May 4, 2017. The first observation points along the strike line were located near the drawdown on July 11, 2017. Six research points were selected in the research area, namely, three strike points ML03, ML06 and ML11, and three dip points MS30, MS33 and MS36, among which ML03 was at the intersection of strike and dip.

Figure 4 Subsidence prediction of the working panel by PIM 

4.3. Data pre-processing

4.3.1 SAR data extraction

The double-track method was used to process the images, and the LOS deformation values were extracted and converted into vertical settlement values by Equation (4). As shown in Fig. 5, the blank area is the empty field formed by the unwrapping errors in the processing of images. The working panel advanced from northwest to southeast, and the displacement area expanded to southeast. Due to large-gradient subsidence during the two periods, some images formed void fields. The mining of the working panel was finished on Dec. 25, 2017, and the images of February 8, 2018 solstice and February 20, 2018 still subsided near the stoping line, and then tended to be stable with slight subsidence, which was within the measurable range of InSAR data.

Figure 5 SAR images processing by dual-track method 

The first leveling observation was conducted on May 4, 2017. The first SAR data were extracted between May 8, 2017 and May 20, 2017. The last leveling observation was carried out on Jan. 7,2019, the last SAR data were extracted between Jan. 10, 2019 and Dec. 29,2018, ∆w _los = [W _{los 2-1} , W _{los 3-2} , W _{los 4-3} ^, L, W _{los n-(n-1)} ].

Assuming that the initial subsidence value was 0 mm, the subsidence values were continuously accumulated until the false value. The subsidence values predicted by PIM were taken as the final subsidence values, and they were accumulated forward from the final SAR data to the false value. Taking Point ML03 as an example, the data of this point were accumulated from May 8, 2017 to July 19, 2017, based on the default initial value of 0 mm, and then the maximum subsidence value predicted by PIM of 1674.3mm was taken as the benchmark, and then accumulated forward from Jan. 10, 2019 to Nov.4, 2017. The data were available in initial period and residual period, while the data of main period was missing.

4.3.2 Interpolation processing of measured data

The measured time of leveling is shown in Table 2. The discrete monitoring data were interpolated to extract the subsidence values consistent with the time of SAR data, W’ _level = [W’ _level-1 , W’ _level-2 , W’ _level-3 ,L, W’ _level-n ]. As shown in Fig. 6, ML03, ML06, ML11, MS30, MS33 and MS36 are interpolated data.

Table 2 Leveling schedule 

Figure 6 Interpolation of ML03, ML06, ML11, MS30, MS33 and MS36  

4.4. Dynamical prediction by the exponential Knothe function model

Two observation lines were arranged along the strike and dip of the working face, with three points on each line. There were 18 monitoring periods. As shown in Fig. 7, the data in main subsidence period were obtained by InSAR-PEK method and compared with the measured data.

Figure 7 Predicted subsidence curves by InSAR-PEK method 

Table 3 shows the comparison of RMSEs of the three groups of tests at each point in the main period and the whole observation period. The minimum RMSE obtained by InSAR-PEK method was 28.1273 mm and the maximum was 91.6846 mm in the main settlement period. The minimum RMSE was 30.9307 mm and the maximum was 58.7326 mm in the whole period. The maximum and minimum predicted RMSEs of the measured data were 21.5255 mm and 148.2012 mm in the main period, respectively. They are 28.2700mm and 87.2497mm in the whole period, respectively. The measured data were eliminated in the main period, and the minimum and maximum RMSEs were 48.6616 mm and 212.6215 mm. While they were 31.6467 mm and 120.5947 mm in the whole period, respectively.

Table 3 Comparison of dynamic predicted results of leveling data and measured data 

The accuracy values of ML03 and ML06 points were slightly lower than the calculated inverse results of the complete measured data, but better than those in main settlement period without the measured data. Among them, the RMSEs obtained by ML11 InSAR-PEK method were better than the predicted results after eliminating the measured data in the main settlement period and the complete measured data. Although the predicted accuracy values ofMS30, MS33 and MS36 were slightly lower than the calculated reverse results of the complete measured data, they were all better than those in the main settlement period excluding measured data. Among them, the prediction accuracy ofMS36 in the main subsidence period was slightly higher than that of the complete measured data, and the accuracy of MS33 point was also slightly higher.

The reason for the low prediction accuracy of ML11 point may be that there is no encryption monitoring, which can not better reflect the subsidence of the whole time period. However, data can be extracted from SAR data in the whole period, and there is still a large amount of data after eliminating false values in the main period, so the predicted accuracy is slightly higher than that of the measured data.

As shown in Figure 8, by comparing the predicted values obtained by the InSAR-PEK method and the extracted values from SAR data with measured values at each point, the extracted values of SAR had a large gap with the measured values in the main subsidence period, which were false values. However, the predicted values by the InSAR-PEK method were close to the measured values, which could meet the mining subsidence prediction requirements.

Figure 8 Comparison of subsidence values obtained by three ways