Dynamic Prediction Model for Progressive Surface Subsidence Based on MMF Time Function

Abstract

It is imperative to timely and accurately predict the progressive surface subsidence caused by coal mining in the context of precision coal mining. However, the existing dynamic prediction methods that use time functions still have limitations, especially in the description of the moments of initiation and maximum subsidence velocity, which hinder their wide application. In this study, we proposed the MMF (Morgan–Mercer–Flodin) time function for predicting progressive surface subsidence based on the model assumptions and formula derivations. MMF time function can resolve the limitations in the description of the moments of initiation and maximum subsidence velocity perfectly. Afterward, we established the dynamic prediction model by combining the probability integral method with the MMF time function. Finally, using the measured subsidence data of working panel 22101 as an example, the accuracy and reliability of the dynamic prediction model was verified. The average RMSE and average relative RMSE (RRMSE) of prediction progressive subsidence using MMF time function are 46.65 mm and 4.63%, respectively. The accuracy is optimal compared with other time functions (for the average RMSE, Logistic time function is 80.57 mm, Gompertz time function is 79.77 mm, and Weibull time function is 90.61 mm; for the average RRMSE, Logistic time function is 7.66%, Gompertz time function is 7.73%, and Weibull time function is 8.62%). The results show that the method proposed in this paper can fully meet the requirements of practical engineering applications, achieve accurate dynamic prediction during the coal mining process, and provide good guidance for surface deformation and building protection.

1. Introduction

China leads the world in coal production and consumes a huge amount of coal every year. The mining of coal resources will inevitably result in overburden movement, damage to buildings (structures), surface subsidence, etc. It will cause a series of social and ecological problems such as life and property safety and sanding of land [1,2,3]. It is urgent to recognize the surface movement and deformation law and make accurate predictions to solve these problems well.

Conventional surface subsidence prediction results often describe the final static surface movement and deformation after the coal mining finished [4,5,6,7], such as the final subsidence W_total in Figure 1a and the final horizontal deformation ε_total in Figure 1b. In practice, it is not hard to consider the final static surface movement and deformation results only. Taking the situation shown in Figure 1 as an example, coal was extracted from open-off cut to position 7, the surface movement basin gradually enlarged (W₁ → W₁₊₂ → … → W_total). When the working panel advances to position 7, the final subsidence basin shows a flat bottom (the subsidence curve about W_total), which is the supercritical mining condition (the maximum surface settlement value does not increase with the size of the mining area and the number of surface maximum subsidence points is more than one). The amount of subsidence in the flat bottom is essentially the same and the horizontal deformation in this region is minimal (the values of ε_total are around the building in Figure 1b).

In Figure 1, the building in this area has been affected or even damaged due to the dynamic deformation that the area underwent (ε₁ → ε₁₊₂ → … → ε_total). The horizontal deformation around the building has dynamic values which include 0, +ε (tensile deformation), −ε (compression deformation) and ε_total. Although this dynamic deformation does not act for a long time, there is a risk that the building may be affected by the deformation and rendered unusable. It is therefore important to determine the amount of ground movement and deformation over different time periods in order to take the necessary protective measures. From the analysis of the above example, it is clear that the actual problem cannot be adequately solved based only on the final state of the surface movement [4,7]. Therefore, the dynamic prediction of mining subsidence has important practical engineering applications [8,9,10]. Moreover, in the context of precision coal mining [11], the prediction subsidence amount at a certain time under water, buildings, railways or roads has higher accuracy requirements during the process of coal mining. Accurate dynamic prediction is crucially important.

Most existing dynamic prediction models are modelled using static prediction model with different time functions [12,13,14]. In China, the static prediction model for mining subsidence is the probability integral method (PIM) [4,15,16]. The key to improving the accuracy of dynamic prediction is to determine the optimal time function [12,17]. Many scholars have conducted relevant studies. In 1952, Knothe proposed the Knothe time function, which describes the dynamic subsidence in the form of a differential equation [18], but it has drawbacks in describing surface subsidence velocity and acceleration [10,14,19]. Considering the flaws of the Knothe time function, Sroka proposed a Sroka-Schober’s function considering the delay of the overburden movement and surface subsidence [20,21,22], but it is difficult to obtain the parameters. Chang and Wang proposed the segmented Knothe time function [23], which does not match the actual values at the maximum and segmented points. Liu presented the power exponent of Knothe time function [24]. Zhang also proposed an improved Knothe time function [19] where the prediction accuracy was also improved. However, the physical meanings of the parameters at the exponential position of two-time functions are not clear. Zhang and Cui developed an optimized segmented Knothe time function [25], but the moment at which the maximum subsidence velocity occurred, it needed to be determined empirically. In addition to Knothe time function and its improved forms, there are other time functions such as Logistic time function [9,13,26], Gompertz time function [27,28], Weibull time function [10,29] and so on.

When synthesizing the above-mentioned studies related to the mining subsidence dynamic prediction based on time functions, it can be found that they almost focus only on a few single points locally, rather than the progressive surface subsidence. As a result, the reliability and practicality of these models are somewhat lacking and do not serve the fundamental purpose of accurately predicting mining subsidence.

In this study, the features of the ideal time function are described, the characteristics and shortcomings of the existing time functions are summarized, and through model analysis and formula derivation it is concluded that the MMF function can accurately describe the dynamic process of surface subsidence. Then, the MMF time function is combined with the PIM to construct a dynamic prediction model for progressive surface subsidence. Finally, an engineering example is applied to verify the accuracy and validity of the proposed model, and the validation is good. This work can be used for the dynamic prediction to provide scientific reference for practical engineering, ensure the safe and effective mining of coal resources and build protection.

2. Models and Methods

2.1. Dynamic Process Analysis of Mining Subsidence

According to movement of a surface point P shown in Figure 1, the surface point P goes through a complicated process. The evolution of its subsidence velocity can be divided into three phases [16,30,31], as shown in Figure 2.

In Figure 2, the initial phase is from the start of surface subsidence until the surface subsidence velocity reaches 1.67 mm/d. The active phase occurs when the velocity exceeds 1.67 mm/d. The weakening phase ends when the cumulative surface subsidence is less than 30 mm for six consecutive months.

In general, the dynamic process of surface subsidence is a limited growth process [32]. The ideal time function used to describe the dynamic process must have the following characteristics:

(1)

When t = 0, the subsidence (w), subsidence velocity (v) and subsidence acceleration (a) of a point on the surface are all zero (w(t) = 0, v(t) = 0, a(t) = 0);

(2)

During the initial phase, the subsidence, subsidence velocity and subsidence acceleration are all small, positive and gradually increasing;

(3)

During the active phase, the subsidence velocity and acceleration are larger, the subsidence acceleration changes gradually to 0 after increasing to the peak, then becomes negative, and gradually decreases after the negative value reaches the peak (0 → v_max → v, 0 → +a_max → 0 → −a_max → −a);

(4)

During the weakening phase, the subsidence velocity gradually becomes smaller and gradually tends to 0, the subsidence acceleration gradually increases from negative and gradually tends to 0 (v → 0, −a → 0).

As shown in Figure 2, a point x_p on the surface is considered. Its subsidence and subsidence velocity at the time t are W(t) and v(t), respectively. It is obvious that the subsidence velocity v(t) is correlated with the average subsidence velocity (W(t)/t). Moreover, after the subsidence velocity of the surface point x_p reaches its maximum, the subsidence velocity v(t) is correlated with the potential subsidence (W₀(x_p) – W(t)), where W₀(x_p) is the final subsidence of the point x_p. According to Figure 2, when the subsidence velocity of the surface point x_p reaches its maximum, the potential subsidence (W₀(x_p) − W(t)) decreases, and correspondingly, the surface subsidence velocity also decreases. Therefore, the subsidence velocity at the time t is correlated with the average subsidence velocity and the potential subsidence ratio is ((W₀(x_p) − W(t))/W₀(x_p)). This summary is the foundation of the surface dynamic subsidence modeling.

2.2. The Existing Time Functions and Their Flaws

2.2.1. Knothe Time Function

In the early 1950s, Polish scholar Knothe proposed the hypothesis of surface subsidence coefficient based on the results of soil compaction experiments. He assumed that the rate of subsidence is proportional to the difference between the final subsidence W₀ and the dynamic subsidence W(t) at the time t [33], shown in Equation (1):

$$\stackrel{•}{W}(t)=\frac{dW(t)}{dt}=c\cdot [W_0-W(t)]$$

(1)

where c is related to overburden lithology. W₀ is the final subsidence.

According to Equation (1) and considering W(0) = 0, W(t) is:

$$W(t)=W_0(1-e^{-ct})$$

(2)

The expressions of v(t) and a(t) are shown in Equations (3) and (4), respectively.

$$v(t)=\frac{dW(t)}{dt}=W_0\cdot e^{-ct}$$

(3)

$$a(t)=\frac{d^{2}W(t)}{d{t}^2}=-W_0\cdot c^2{e}^{-ct}$$

(4)

Considering W₀ is 1500 mm, the influence coefficient c is 0.04 and the curves of the Knothe time function are plotted in Figure 3a,b, respectively.

As can be seen from Figure 3a, the shape of subsidence-time curve is not the “S” shape and cannot reflect the three phases of surface subsidence completely. According to Figure 3b, it can be seen that the surface subsidence velocity is not 0 → +v_max → 0, while the subsidence acceleration is not 0 → +a_max → 0 → −a_max → 0. Therefore, Knothe time function cannot effectively describe the dynamic subsidence.

2.2.2. Segmented Knothe Time Function

Knothe time function does not correspond to the law of dynamic subsidence [12,14]. In 2003, Chang and Wang proposed a segmented Knothe time function which assumes that the maximum subsidence velocity occurs at half the total time and that the subsidence at that moment is half the total subsidence [23], as shown in Equation (5).

$$W(t)=\{\begin{array}{l}0.5W_0[e^{-c(\tau -t)}-e^{-c\tau }], 0<t\le \tau \\ 0.5W_0[2-e^{-c(t-\tau )}-e^{-c\tau }], \tau <t<T\end{array}$$

(5)

where τ is the moment of maximum subsidence velocity, T is the total subsidence time.

The expression of subsidence velocity v(t) is shown in Equation (6).

$$v(t)=\{\begin{array}{c}0.5W_0[c{e}^{-c(\tau -t)}],0<t\le \tau \\ 0.5W_0[c{e}^{-c(t-\tau )}],\tau <t<T\end{array}$$

(6)

W₀ is 1500 mm, the influence coefficient c is 0.1, the total subsidence time T is 150 days and τ is 35 days (τ ≈ T/4) [25]. The subsidence, velocity and acceleration are plotted in Figure 4a,b, respectively.

According to Figure 4a,b, the subsidence has the shape of “S” and can improve the prediction accuracy. While the process of acceleration and deceleration of subsidence is symmetrical, it does not reflect the feature of rapidity during the active phase and the feature of long-term during the weakening phase.

2.2.3. Optimized Segmented Knothe Time Function

To overcome the theoretical deficiency of the segmented Knothe time function, Zhang and Cui developed the optimized segmented Knothe time function [25,34].

$$W(t)=\{\begin{array}{l}\frac{1}{2}{W}_0[\frac{t-\tau (1-e^{ct})}{\tau e^{c\tau }}],0<t\le \tau \\ W_0[1-\frac{1}{2}{e}^{c(\tau -t)}],\tau <t<T\end{array}$$

(7)

where τ and T are the same as those in Equation (5).

In fact, the time-subsidence curve is asymmetrical, with a short acceleration period and long deceleration period, and the moment when maximum subsidence velocity occurs is always before half of the total subsidence time T [7,20,28,32]. The optimized segmented Knothe time function still makes the first half (acceleration period) and second half (deceleration period) of the time-subsidence curve strictly symmetrical. It is not quite consistent with the actual observed data.

2.2.4. Other Time Functions

In addition to the Knothe time function and its improved forms, there are other time functions: Logistic time function, Gompertz time function, Sigmoidal Function [35] and so on [22,36,37].

(1)

Logistic time function

The Logistic time function has been widely used in ecology, demography and other fields. It reflects the process by which things occur, develop, mature and tend to saturation [13,26,38,39]. It is expressed as:

$$W(t)=\frac{W_0}{1+f{e}^{-gt}}$$

(8)

where g, f are the coefficients related to overburden lithology. W₀ and W(t) are the same as above.

The expressions of v(t) and a(t) for Logistic time function are shown in Equations (9) and (10), respectively:

$$v(t)=\frac{W_{0}gf{e}^{-gt}}{{(1+f{e}^{-gt})}^2}$$

(9)

$$a(t)=\frac{W_0{g}^{2}f{e}^{-gt}(f{e}^{-gt}-1)}{{(1+f{e}^{-gt})}^2}$$

(10)

The time-subsidence curve of the Logistic time function also has a “S” shape, while when t = 0, W(t) = W₀/(1 + f) ≠ 0, v(t) ≠ 0, a(t) ≠ 0. Moreover, considering a(τ) = 0, then e^−gτ = 1/f, W(τ) = W₀/2, it means the subsidence amount of the Logistic time function at the inflection point is constant, which is equal to 1/2 of the final subsidence. These cases do not accord with the measured data.

(2)

Gompertz time function

The Gompertz time function is proposed by the British mathematician B. Gomepertz. It is applied to the law of the extinction of animals and plants, and is now mostly used for the settlement of embankments or foundations under the influence of non-mining [27,28,40], which is expressed as

$$W(t)=W_0\exp (-e^{-k(t-h)})$$

(11)

where k and h are parameters that determine the shape and the position, respectively.

Like the Logistic time function, although the time-subsidence curve of the Gompertz time function also has a “S” shape, when t = 0, W(t) = W₀·exp(−e^kh) ≠ 0, v(t) ≠ 0, a(t) ≠ 0. Moreover, considering a(τ) = 0, then W(τ) = W₀/e, it means the subsidence amount at the inflection point is constant, which is equal to 1/e of the final subsidence. Obviously, these cases do not accord with the measured data.

(3)

Weibull time function

The Weibull time function is the famous S-shaped growth curve model proposed by Swedish mathematician Weibull in 1951 and is particularly suitable for mining conditions, which have thicker mining depths and thicker overburden [10,29,41], which is expressed as:

$$W(t)=W_0(1-e^{-n{t}^p})$$

(12)

where n, p are parameters related to the properties of overlying strata.

Weibull time function has the following characteristics:

(a)

Compared to the Knothe time function (Equation (2)) the Weibull time function (Equation (12)), is an improvement of the Knothe time function (n = c, p= 1).

(b)

The convergence rate of this function is fast, and the tail properties of subsidence during the weakening period is insufficient for description [42].

(c)

The physical meaning of the parameters is not clear and difficult for promotion [42].

Based on a large amount of measured data and the analysis above, it is demonstrated that the time functions analyzed in this paper and other common time functions do not exactly match the actual dynamic process of surface subsidence [23,39,42]. Although these functions can describe the dynamic subsidence to some extent, these functions either do not reflect the change law of velocity and acceleration, or the physical meaning of the parameters in the functions is unclear and does not reflect the change of each parameter with geological mining conditions.

Therefore, there is an urgent need to determine an optimal time function that can accurately describe the dynamic progressive subsidence in mining areas in a timely and accurate manner.

2.3. Establishment of MMF Time Function

Based on the analysis of dynamic surface movement shown in Figure 2, we found that the subsidence velocity at time t is correlated with the average subsidence velocity and the potential subsidence ratio ((W₀(x_p) – W(t))/W₀(x_p)), shown in Equation (13):

$$\stackrel{•}{W}(t)=\frac{dW(t)}{dt}=b\cdot \frac{W(t)-W(0)}{t-0}\cdot [\frac{W_0(x_p)-W(t)}{W_0(x_p)}]=b\cdot \frac{W(t)}{t}\cdot [\frac{W_0(x_p)-W(t)}{W_0(x_p)}]$$

(13)

where b is a parameter related to geological and mining technical parameters.

Rearranging Equation (13),we have:

$$\frac{W_0(x_p)}{W(t)\cdot \left[W_0(x_p)-W(t)\right]}dW(t)=\frac{b}{t}dt$$

(14)

Integrating both sides of Equation (14):

$${\displaystyle \int \frac{W_0(x_p)}{W(t)\cdot \left[W_0(x_p)-W(t)\right]}}dW={\displaystyle \int \frac{b}{t}}dt$$

(15)

Decomposing the rational fraction of Equation (15):

$${\displaystyle \int [\frac{1}{W(t)}+\frac{1}{W_0(x_p)-W(t)}]}dW={\displaystyle \int \frac{b}{t}}dt$$

(16)

Integrating both sides of Equation (16):

$$\ln W(t)-\ln [W_0(x_p)-W(t)]+C_1=b\ln t$$

(17)

where C₁ is a constant.

Let C₁ = ln a, we have:

$$\ln W(t)-\ln [W_0(x_p)-W(t)]+\ln a=\ln t^b$$

(18)

Removing the logarithmic symbol of Equation (18):

$$\frac{W(t)\cdot a}{W_0(x_p)-W(t)}=t^b$$

(19)

Adding 1 to both sides of Equation (19) and simplifying:

$$W(t)=\frac{W_0(x_p)\cdot t^b}{a+t^b}=W_0(x_p)\cdot \phi (t)$$

(20)

The form of W(t) in Equation (20) is similar to the Morgan–Mercer–Flodin growth model [43]. φ(t) in Equation (20) is the MMF time function for the surface dynamic subsidence prediction proposed in this paper.

The subsidence velocity and acceleration expressions of MMF time function are shown in Equations (21) and (22), respectively:

$$v(t)=\frac{W_0(x_p)\cdot a\cdot b\cdot t^{b-1}}{{(a+t^b)}^2}$$

(21)

$$a(t)=\frac{W_0(x_p)\cdot a\cdot b\cdot t^{b-2}}{{(a+t^b)}^3}[(b-1)(a+t^b)-2b{t}^b]$$

(22)

The curves of the MMF time function are shown in Figure 5.

Let a(t) = 0, and then t₁ = 0, t₂ = $\sqrt[b]{a(b-1)/(b+1)}$ and t₃ = ∞. When 0 ≤ t ≤ $\sqrt[b]{a(b-1)/(b+1)}$, a(t) > 0; when $\sqrt[b]{a(b-1)/(b+1)}$ < t < ∞, a(t) < 0. It means the subsidence velocity increases at the interval [0, $\sqrt[b]{a(b-1)/(b+1)}$] and decreases at the interval ($\sqrt[b]{a(b-1)/(b+1)}$, ∞), where the subsidence velocity reaches the maximum. It is proven that there is an inflection point in the curve of the MMF time function. Moreover, the inflection point is not fixed and is correlated with the parameters a and b, this characteristic is more accord with the practical situation of subsidence engineering.

Let b = 5.5, a = 1.0 × 10⁷, 1.0 × 10⁸, 1.0 × 10⁹, 1.0 × 10¹⁰, respectively; the effect of the variation of parameter a is on the curves of subsidence and velocity is analyzed, as shown in Figure 6a,b. Let a = 1.0 × 10⁸, b = 4.5, 5.0, 5.5, 6.0, respectively; the variation of parameter b is on the curves of subsidence and subsidence velocity is analyzed, as shown in Figure 6c,d.

As shown in Figure 6, the variation of parameters a and b has the opposite effect on the curves. For a surface point, as the parameter a increases, the subsidence velocity decreases and the time required from the beginning of subsidence to stabilization increases significantly. As the parameter b increases, the subsidence velocity of increases significantly, and the time required from the beginning of subsidence to the stabilization is significantly shortened. Therefore, it can be considered that parameters a and b reflect the subsidence velocity together, and the values of these two parameters depend on the nature of the overlying rock layers in the mining area.

With suitable parameters a and b, the MMF time function can provide a more realistic description of the dynamic process of surface subsidence.

2.4. Determination of parameters for MMF Time Function

There is a close relationship between the parameters in dynamic prediction models and the geological and mining technical conditions, and the dynamic subsidence is different under different geological and mining technical conditions [8,12,14,28]. Therefore, it is necessary to analyze the statistical relationship between model parameters and geological and mining technical parameters in order to guide the actual modeling work.

The geological and mining parameters include mining thickness M, average mining depth H₀, loose layer thickness h, advance velocity of a working panel v and subsidence coefficient q.

The parameters of MMF time function are a and b. In order to obtain the statistical relationship between the two parameters (a, b) and geological and mining technical parameters (M, H₀, h, v), the geological and mining technical parameters of sixteen working panels were selected for research. MMF time function was also applied to measured data of the maximum subsidence point (MSP) on the sixteen strike main sections of the working panels [8,14]. The method for determining parameters a and b is nonlinear least squares algorithm based on MATLAB curve fitting toolbox. The relevant MMF time function parameters of the maximum subsidence points and geological and mining technical parameters of each working panel are shown in

Table 1. Parameters of MMF time function and geologic and mining conditions.

Name of Working Panels	Mining Height M, m	Average Mining Depth H0, m	Thickness of Loose Layer h, m	Advance Velocity v, m/d	Subsidence Coefficient q	ln a	b
RL7222	2.2	390	285	1.64	1.08	18.02	3.83
BZ1203	2.3	255	42.5	7.19	1.1	19.25	5.24
BZ4303	4.6	217	52.4	2.30	0.91	23.61	5.87
BZ4410	5.78	314	50	5.06	0.81	32.31	8.40
LJ1203	4.19	327	45	2.30	0.77	19.74	4.57
LJ2201	3.93	311	40	2.60	0.8	12.25	3.58
LJ4412	6.53	439	35	2.98	0.55	18.56	4.83
LJ4108	6.48	387	38	6.35	0.71	30.16	7.94
WL11206	2	285	26	3.37	0.88	17.78	4.17
WL1201	1.4	66.5	24	2.03	1.09	22.30	5.75
WL10103	2.03	282	30	1.56	0.95	16.62	3.60
XL5306	7.25	420	183	2.31	0.81	26.70	5.48
JZ3702	1.6	330	50.5	2.50	0.35	12.03	2.37
MM2311D	2	190	30	0.91	0.64	17.52	4.07
MM3212	0.93	557	34	2.25	0.44	10.45	2.09
ZZ6210	2.8	396	60	2.57	0.65	20.49	4.27

. Considering the magnitude variance between parameter a and parameter b, the relationship between the logarithm value of parameter a (ln a) and geological mining conditions is statistically analyzed.

For the exploitation cases in Table 1, where q is greater than 1, which can be caused by reactivation of goafs. The reactivation of goafs is related to the ratio of the thickness of the loose layer to the thickness of the bedrock; the larger the ratio, the more likely the goafs are activated. Furthermore, when the thickness of the loose layer exceeds 100 m, the soil has additional settlement (settlement caused by consolidation and compression of the soil and water table drop), which also results in q greater than 1.

According to Table 1, we found that the MMF model parameters (a, b) are mainly related to the ratio of bedrock thickness to mining thickness (H₀ – h)/M, advance velocity v and overburden properties. The smaller the ratio of bedrock thickness to mining thickness (i.e., the larger M/(H₀ – h)), the larger the advance velocity; the smaller the average solidity coefficient of the overlying rock layers, the more intense the surface movement and deformation, and accordingly the larger the values of a and b of MMF time function. The relationship between parameters a and b and M/(H₀ – h), v and q is shown in Figure 7.

According to Figure 7, the statistical relationship of parameters a, b and M/(H₀ – h), v and q is:

$$\{\begin{array}{l}\ln a=157.7452\frac{M}{H_0-h}vq+17.6518\\ b=48.6275\frac{M}{H_0-h}vq+4.3393\end{array}$$

(23)

In Equation (23), the subsidence coefficient q reflects the nature of the overburden. If the overburden is softer, q is larger, while the overburden is harder, q is smaller; other parameters are the same above.

In practice, if the geological mining parameters (i.e., M, H₀, h, v, q) are outside the range of values in Table 1, it is advisable to obtain the model parameters based on the actual data of similar geological and mining technical conditions.

2.5. The Final State Prediction Method

In China, the most popular final state prediction method is the probability integral method (PIM), which is based on the stochastic medium theory [4]. Consider horizontal or gently inclined coal seam mining, if it is a rectangular working panel, the subsidence prediction model based on coordinate transformation can be constructed along the strike or inclined main section of the working panel, as shown in Equation (24):

$$\{\begin{array}{l}{W}^0(x)=\frac{W_0}{2}\left\{erf\left[\frac{\sqrt{\pi }}{r}(x-s_3)\right]-erf\left[\frac{\sqrt{\pi }}{r}(x-D_3+s_4)\right]\right\}\\ W^0(y)=\frac{W_0}{2}\left\{erf\left[\frac{\sqrt{\pi }}{r_1}(y+H_1\cot {\theta }_0-s_1\cos \alpha )\right]-erf\left[\frac{\sqrt{\pi }}{r_2}(y-(D_1-s_2)\cos \alpha +H_2\cot {\theta }_0)\right]\right\}\end{array}$$

(24)

where W⁰(x) and W⁰(y) are the final subsidence of strike main section and inclined main section, respectively. W₀ is the maximum surface subsidence, W₀ = Mqcos α, M is the mining thickness, q is the subsidence coefficient, or the ratio of maximum surface subsidence value to mining thickness under critical mining conditions and α is the coal seam dip angle. H₁ and H₂ are the mining depth in downhill and uphill directions, respectively; r, r₁ and r₂ are the main influence radius of strike, downhill and uphill directions, respectively; s₃, s₄, s₁ and s₂ are the inflection point offsets of strike left, right, downhill and uphill directions, respectively. D₁ and D₃ are mining lengths along the strike and inclined directions, respectively; θ₀ is the mining influence propagation angle.

The coordinate system used in the construction of the mathematical model is shown in Figure 8a,b.

2.6. Dynamic Prediction Model for Progressive Surface Subsidence

The dynamic prediction model based on MMF function can be determined by multiplying the surface subsidence of each point by the MMF time function. During the process of calculating, the mining area is divided into n independent mining units according to the time series. First, the surface subsidence generated by each unit is calculated. Then, the surface subsidence generated by all units is superimposed to obtain the total surface subsidence.

The start of mining is defined as the start time, and the advance velocity of each unit is v₁, v₂, …, v_n, and the advance times of each unit are t₁, t₂, …, t_n, respectively, then the size of each unit is v₁t₁, v₂t₂, …, v_nt_n; after the working panel is fully extracted, the sum of surface subsidence generated by all mining units is shown in Figure 9. W₁, W₂, W₃, …, W_n is the subsidence contributed by each mining unit, and W_t is the total subsidence affected by mining. Considering the amount of surface subsidence W_t at time t, when the first mining unit is mined, it experiences the longest time from 0 to t, and has the maximum contribution to the surface subsidence, as W₁ in Figure 9. Accordingly, after the second mining unit is mined, it experiences the time from t₁ to t, it has the contribution to the surface subsidence, as W₂ in Figure 9. After the kth unit is extracted (t − t₁ − t₂ − … − t_k−1), it has the minimal contribution to surface subsidence, like W_n in Figure 9. After all n units are extracted, the surface subsidence is W_t (Figure 9).

When the coal seam is horizontal (α = 0°), the dynamic prediction model along the strike main section is

$$\begin{array}{ll}W(x,t)& =\phi (t)W_1^0(x)+\phi (t-t_1)W_2^0(x)+\cdots +\phi (t-t_1-t_2-\cdots t_n)W_n^0(x)\\ & =\phi (t)[W(x-s_3)-W(x-v_1{t}_1)]+\phi (t-t_1)[W(x-v_1{t}_1)-W(x-v_1{t}_1-v_2{t}_2)]\\ & +\phi (t-t_1-t_2)[W(x-v_1{t}_1-v_2{t}_2)-W(x-v_1{t}_1-v_2{t}_2-v_3{t}_3)]\\ & +\cdots +\phi (t-t_1-t_2-\cdots t_{n-1})\{W(x-v_1{t}_1-\cdots -v_{n-1}{t}_{n-1})-W[(x-(v_1{t}_1+v_2{t}_2+\cdots +v_n{t}_n-s_4)]\}\end{array}$$

(25)

where W(x – Σv_it_i) can be calculated according to the first formula in Equation (24), and s₃, s₄ are the inflection point offsets of strike left and right directions. φ(t) is MMF time function in Equation (20).

When the inclined coal seam is mined, the dynamic prediction model along the inclined main section is:

$$W(y,t)=\phi (t)W_1^0(y)+\phi (t-t_1)W_2^0(y)+\cdots +\phi (t-t_1-t_2-\cdots t_n)W_n^0(y)$$

(26)

where W_i⁰(y) is the second formula in Equation (24) and φ(t) is MMF time function in Equation (20).

Equation (25) and (26) are the progressive dynamic prediction model along strike and inclined main sections, respectively.

3. Case Study

3.1. Study Area

The study working panel is 22101 in Taiyuan, Shanxi Province, China, as shown in Figure 10.

The overlying rocks in the working panel 22101 are the Shihezi Group and Shanxi Group strata, mainly composed of fine sandstone, sandy mudstone and sand shale interbedded, medium-hard lithology, with a thickness of about 140 m. The Quaternary loess has a thickness of about 70 m. The soil has a large sand content, is not strongly bonded and has vertical joints; thus, the tensile and shear strength is low.

The working panel adopted an all-collapse method. The mining height is M = 3.18 m; the lengths of the working panel are 440 m (D₁) and 140 m (D₃), respectively; the dip angle is α = 1~3°, which can be considered horizontal seam mining; the average mining depth is H₀ = 220 m; the average thickness of loose layers h = 80 m; the average advance velocity is v = 2.2 m/d. According to Equation (23), the parameters of MMF time function which include ln a and b, are 22.38 and 5.80, respectively.

Two observation lines were laid and there were 45 monitoring points, as shown in Figure 11. The traverse measurement and the trigonometric elevation measurement were performed to gain the plane position (x, y) and the elevation (h) of each point. A Leica T402 total station was used for measurement, the angular nominal accuracy was 2” and the distance accuracy was 3 mm + 2 ppm.

3.2. MMF Time Function Validation with Measured Data

Data observation of mining-induced surface subsidence included a traverse survey and a trigonometric elevation survey. Fifteen days after the first observation with a total station, the workings began to advance. The final coordinate survey was carried out 325 days after the first observation. Fifteen observations were conducted at 26 points on the strike observation line, as shown in Figure 12. Five observations were made on the tendency observation line, however, there are a lack of data, so the dynamic prediction model proposed in this paper is validated and analyzed mainly based on the data on the strike observation line. As demonstrated in Figure 12, from 309 days after the first observation, the surface movement stabilized and the maximum subsidence reached 1898 mm, which is located at No. 18.

To illustrate the feasibility and applicability of the MMF time function for predicting the progressive surface subsidence, the MMF time function was fitted to the measured settlement data. The measured data of six observation points on the strike line of the working panel 22101 were selected, including No. 10, No. 13, No. 16, No. 19, No. 22 and No. 25.

In order to eliminate differences of subsidence values at different points and to highlight the common feature of subsidence over time, we transformed the subsidence of each point to a subsidence ratio, which is equal to W_i/W₀, and took the values in the range 0–1. The parameters a and b of MMF time function are solved by nonlinear least squares algorithm through loop iteration. Considering the large discrepancy in order of magnitude between parameter a and parameter b, we replaced parameter a with logarithm a.

The fitting results for the six points are shown In Figure 13.

Each MMF fitting result is consistent with the corresponding time-subsidence ratio data, with R² greater than 0.97 and even greater than 0.99 for the majority of R². This suggests that the MMF time function is well suited for predicting dynamic subsidence at surface points caused by coal mining.

Furthermore, the surface subsidence when the maximum subsidence velocity occurs is not always about half of the final subsidence, like the points near the boundary of the working panel (No. 10, No. 13, No. 16). In addition, regardless of the time function and coordinate system used, the parameters at different locations are different.

3.3. Dynamic Prediction Results

Based on the progressive surface subsidence parameters calculation method and dynamic prediction model above, surface subsidence was predicted at 52, 84, 114 and 178 d after the first observation for the working panel 22101. Considering that the surface subsidence had basically stabilized at 178 days after the first observation, the measured data at 178 d were interpolated from the data at 164d and the data at 192 d.

The parameters of the probability integral method (PIM) were obtained based on the final subsidence curve (the subsidence curve of 325 days after the first observation) [16], as shown in

Table 2. Prediction parameters of the PIM and number of exploitation times for working panel 22101.

q	tan β	θ0/°	s1/m	s2/m	s3/m	s4/m
0.59	2.19	83.19	−2.43	−2.43	7.57	7.57
The number of exploitation times:			First exploitation

.

During the process of dynamic prediction for four observation subsidence curves, each unit length along the strike main section direction can be 0.1H₀ = 22 m [44]. The number of units for the prediction of four subsidence curves (A → B (52 days, advance 100 m), A → D (84 days, advance 150 m), A → E (114 days, advance 190 m), A → G (178 days, advance 410 m)) can be n₁ = 5, n₂ = 7, n₃ = 9, n₄ = 19, respectively. Equation (25) was applied for predicting the dynamic subsidence for the selected four subsidence. The prediction results are shown in Figure 14.

According to Figure 12 and Figure 14, with the advance of the working panel, the range of subsidence gradually increases, and eventually a flat bottom appears. During this process, the surface deformation is also changing, so the dynamic prediction is of great value for grasping the surface deformation in real time, enhancing the understanding of the dynamic movement and deformation, reasonably making the mining plan and protecting the surface buildings and structures.

3.4. Model Validation

The accuracy of the model proposed in this paper was examined by comparing the prediction results with the measured data along the strike observation line at four prediction times for the working panel 22101. In addition, the prediction results using the time functions in Section 2.2.4 at the same prediction times were also compared in order to provide a clearer description of the prediction accuracy and advantages of the MMF time function. The comparison curves are shown in Figure 15.

A comparison of the results predicted by four time functions in Figure 15 shows that the difference between the progressive surface subsidence predicted by MMF time function and measured subsidence is much smaller, particularly at the extreme values of subsidence for each time period and near the boundaries.

The prediction accuracy was evaluated in order to compare the prediction effectiveness of each time function more intuitively. The root mean square error (RMSE) (Equation (27)) and relative root mean square error (RRMSE) (Equation (28)) of the progressive surface subsidence for each time period were calculated.

$$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum _{i=1}^m(W_{Pi}^{}-W_{Mi}^{})}}{m}}$$

(27)

$$\mathrm{RRMSE}=\frac{\mathrm{RMSE}}{W_{\max }^{}}$$

(28)

where W_Pi is the prediction subsidence of the ith point; W_Mi is the measured subsidence of the ith point; W_max represents the maximum value of the measured subsidence for each time period; m is the number of observation points (m = 26).

The accuracy of the progressive surface subsidence predicted by four time functions is shown in

Table 3. Accuracy comparison for prediction results among four time functions.

Time Period/d	Advance Distance/m	RMSE/mm				RRMSE
		Logistic	Gompertz	Weibull	MMF	Logistic	Gompertz	Weibull	MMF
52	100	3.55	4.73	3.76	2.79	7.25%	9.66%	7.67%	5.70%
84	150	92.40	57.38	107.75	36.10	8.83%	5.48%	10.29%	3.45%
114	190	114.58	85.10	150.42	67.16	8.32%	6.18%	10.92%	4.88%
178	410	111.75	171.88	100.49	80.53	6.24%	9.59%	5.61%	4.49%
Average		80.57	79.77	90.61	46.65	7.66%	7.73%	8.62%	4.63%
Increased accuracy		33.92	33.12	43.96		3.03%	3.10%	3.99%

.

With surface subsidence increases, the RMSE and RRMSE gradually increase. The MMF time function outperforms other time functions in all time periods in terms of prediction accuracy. The mean values of RMSE-Logistic and RRMSE-Logistic are 80.57 mm and 7.66%, respectively, the mean values of RMSE-Gompertz and RRMSE-Gompertz are 79.77 mm and 7.73%, respectively, and the mean values of RMSE-Weibull and RRMSE-Weibull are 90.61 mm and 8.62%, respectively. While the mean values of RMSE-MMF and RRMSE-MMF are 46.65 mm and 4.63%, respectively, the accuracy of RMSE is improved by 33.92 mm, 33.12 mm, 43.96 mm, respectively, and the accuracy of RRMSE is improved by 3.03%, 3.10%, 3.99%, respectively.

Combined with the analysis of the subsidence data above, the average RMSE of the progressive surface subsidence predicted by MMF time function is 46.65 mm and an average RRMSE of 4.63%. This meets the requirements of practical engineering applications. It can be used fully for dynamic prediction of mining subsidence.

4. Conclusions

An accurate and effective dynamic prediction of progressive surface subsidence is important and indispensable for mine production and planning and for the protection of surface structures. To this end, this paper investigates the dynamic prediction model based on MMF time function and the conclusions are summarized here.

The MMF time function for progressive prediction was proposed based on the model assumption and formula derivation, and the characteristics of the MMF time function parameters for describing surface dynamic subsidence were clarified. Subsequently, the variation relationship between the MMF time function parameters and the geological mining technical conditions was analyzed and fitting formulas were given. The MMF time function was then combined with the probability integral method (PIM), which is applied for predicting the final state surface subsidence to develop a dynamic progressive prediction model.

For the proposed model validation, the measured data of working panel 22101 were used as an engineering example. The average RMSE and average RRMSE, which used the MMF time function, were 46.65 mm and 4.63%, respectively. Compared with Logistic time function, Gompertz time function and Weibull time function, the accuracy of RMSE is improved by 33.92 mm, 33.12 mm, 43.96 mm, respectively, and the accuracy of RRMSE is improved by 3.03%, 3.10%, 3.99%, respectively. The MMF time function outperforms other time functions in terms of prediction accuracy at all time periods.

The results show that the accuracy of the model proposed in this paper can meet the requirements of practical engineering applications. This model can assess surface building damage and determine repair time so that appropriate repair measures can be given to reduce economic losses and alleviate the conflict between mining companies and residents. Therefore, the method is expected to be widely used.

Nevertheless, it should be pointed out that the prediction model mentioned in this paper only considers the mining case of a single coal seam and does not consider the mining case of multiple coal seams. However, this is a common case in actual production and will affect dynamic surface subsidence prediction. Therefore, predicting progressive subsidence under multiple coal seams is a future endeavor for our country.