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

: 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 veri ﬁ ed. 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, Gomper tz time function is 79.77 mm, and Weibull time function is 90.61 mm; for the average RRMSE, Logistic time function is 7.66%, Gomper tz 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.


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 Wtotal 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  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).

(a) (b)
Figure 1.Dynamic changes about surface subsidence and horizontal deformation caused by coal mining (from position 1 to position 7): (a) Dynamic changes of surface subsidence, the curve of W1 is formed when the working panel advances from open-off cut to position 1, the curve of W1+2 is formed when advancing to position 2, the others are the same, the curve of Wtotal is the final subsidence when advancing to position 7; (b) Dynamic changes of horizontal deformation.the curve of ε1 is formed when the working panel advances from open-off cut to position 1, the curve of ε1+2 is formed when advancing to position 2, the others are the same, the curve of εtotal is the final horizontal deformation when advancing to position 7.
In Figure 1, the building in this area has been affected or even damaged due to the dynamic deformation that the area underwent (ε1 → ε1+2 → … → ε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.

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 ) 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 xp 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 xp reaches its maximum, the subsidence velocity v(t) is correlated with the potential subsidence (W0(xp) -W(t)), where W0(xp) is the final subsidence of the point xp.According to Figure 2, when the subsidence velocity of the surface point xp reaches its maximum, the potential subsidence (W0(xp) − 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 ((W0(xp) − W(t))/W0(xp)).This summary is the foundation of the surface dynamic subsidence modeling.

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 W0 and the dynamic subsidence W(t) at the time t [33], shown in Equation (1): where c is related to overburden lithology.W0 is the final subsidence.
Considering W0 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 → +vmax → 0, while the subsidence acceleration is not 0 → +amax → 0 → −amax → 0. Therefore, Knothe time function cannot effectively describe the dynamic subsidence.

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 ) acceleration velocity 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).
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).
W0 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.

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].
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.

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: where g, f are the coefficients related to overburden lithology.W0 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: ) gt gt gt W g fe fe a t fe The time-subsidence curve of the Logistic time function also has a "S" shape, while when t = 0, W(t) = W0/(1 + f) ≠ 0, v(t) ≠ 0, a(t) ≠ 0.Moreover, considering a(τ) = 0, then e −gτ = 1/f, W(τ) = W0/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 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) = W0•exp(−e kh ) ≠ 0, v(t) ≠ 0, a(t) ≠ 0.Moreover, considering a(τ) = 0, then W(τ) = W0/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: 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.

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 ((W0(xp) − W(t))/W0(xp)), shown in Equation ( 13): where b is a parameter related to geological and mining technical parameters.
Rearranging Equation ( 13),we have: Integrating both sides of Equation ( 14): Decomposing the rational fraction of Equation ( 15): Integrating both sides of Equation ( 16): where C1 is a constant.let C1 = ln a, we have: Removing the logarithmic symbol of Equation ( 18): Adding 1 to both sides of Equation ( 19) and simplifying: 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: The curves of the MMF time function are shown in Figure 5.Let a(t) = 0, and then t1 = 0, t2 =  ] and decreases at the interval ( , ∞), 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 7 , 1.0 × 10 8 , 1.0 × 10 9 , 1.0 × 10 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 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.

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 H0, 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, H0, 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.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 (H0 − h)/M, advance velocity v and overburden properties.The smaller the ratio of bedrock thickness to mining thickness (i.e., the larger M/(H0 − 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/(H0 − h), v and q is shown in Figure 7.According to Figure 7, the statistical relationship of parameters a, b and M/(H0 − h), v and q is: 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, H0, 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.

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): where W 0 (x) and W 0 (y) are the final subsidence of strike main section and inclined main section, respectively.W0 is the maximum surface subsidence, W0 = 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.H1 and H2 are the mining depth in downhill and uphill directions, respectively; r, r1 and r2 are the main influence radius of strike, downhill and uphill directions, respectively; s3, s4, s1 and s2 are the inflection point offsets of strike left, right, downhill and uphill directions, respectively.D1 and D3 are mining lengths along the strike and inclined directions, respectively; θ0 is the mining influence propagation angle.
The coordinate system used in the construction of the mathematical model is shown in Figure 8a

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 v1, v2, …, vn, and the advance times of each unit are t1, t2, …, tn, respectively, then the size of each unit is v1t1, v2t2, …, vntn; after the working panel is fully extracted, the sum of surface subsidence generated by all mining units is shown in Figure 9. W1, W2, W3, …, Wn is the subsidence contributed by each mining unit, and Wt is the total subsidence affected by mining.Considering the amount of surface subsidence Wt 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 W1 in Figure 9. Accordingly, after the second mining unit is mined, it experiences the time from t1 to t, it has the contribution to the surface subsidence, as W2 in Figure 9.After the kth unit is extracted (t − t1 − t2 − … − tk-1), it has the minimal contribution to surface subsidence, like Wn in Figure 9.After all n units are extracted, the surface subsidence is Wt (Figure 9).When the coal seam is horizontal (α = 0°), the dynamic prediction model along the strike main section is where W(x − Σviti) can be calculated according to the first formula in Equation ( 24), and s3, s4 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: where Wi 0 (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.

Study Area
The study working panel is 22101 in Taiyuan, Shanxi Province, China, as shown in Figure 10.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.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.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 Wi/W0, 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.

MMF Time Function Validation with Measured Data
The fitting results for the six points are shown in Figure 13.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.

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. 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.1H0 = 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 n1 = 5, n2 = 7, n3 = 9, n4 = 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 Figures 12 and 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.where WPi is the prediction subsidence of the ith point; WMi is the measured subsidence of the ith point; Wmax 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. 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.

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.
as an example, coal was extracted from open-off cut to position 7, the surface movement basin gradually enlarged (W1 → W1+2 → … → Wtotal).When the working panel advances to position 7, the final subsidence basin shows a flat bottom (the subsidence curve about Wtotal), which is the supercritical mining Citation: Zhou, B.; Yan, Y.; Kang, J. Dynamic Prediction Model for Progressive Surface Subsidence Based on MMF Time Function.

Figure 2 .
Figure 2. Surface movement process of a point.(The initial phase is phase I.The active phase is divided into phase II and phase III based on the maximum subsidence velocity.The weakening phase is phase IV.)

Figure 3 .
Figure 3. (a) The subsidence of Knothe time function; (b) The velocity and acceleration of Knothe time function.

Figure 4 .
Figure 4. (a) The subsidence and velocity of the segmented Knothe time function; (b) The acceleration of the segmented Knothe time function.

Figure 5 .
Figure 5. Dynamic subsidence, velocity and acceleration of MMF time function.

bFigure 6 .
Figure 6.Variation of surface subsidence and subsidence velocity curves with parameters a and b based on MMF time function; (a) Variation of parameter a on MMF subsidence curve; (b) Variation of parameter a on MMF subsidence velocity curve; (c) Variation of parameter b on MMF subsidence curve; (d) Variation of parameter b on MMF subsidence velocity curve.

Figure 8 .
Figure 8.(a) Coordinate system along strike main section; (b) coordinate system along inclined main section.

Figure 10 .
Figure 10.(a) Location of Shanxi province in China; (b) Location of Taiyuan city in Shanxi province; (c) Location of the working panel 22101.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 (D1) and 140 m (D3), respectively; the dip angle is α = 1~3°, which can be considered horizontal seam mining; the average mining depth is H0 = 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 Figure11.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.

Figure 11 .
Figure 11.The Layout of working panel 22101.
100 m (52 days after the first observation) A→C: advance 130 m (67 days after the first observation) A→D: advance 150 m (84 days after the first observation) A→E: advance 190 m (114 days after the first observation) A→F: advance 310 m (149 days after the first observation) A→G: advance 410 m (178 days after the first observation)

Figure 12 .
Figure 12.Measured subsidence curves of the strike observation line.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 Wi/W0, 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 Figure13.
open-off cut (m) 15 days after the first observation 52 days after the first observation 84 days after the first observation 93 days after the first observation 100 days after the first observation 107 days after the first observation 114 days after the first observation 129 days after the first observation 136 days after the first observation 149 days after the first observation 164 days after the first observation 192 days after the first observation 265 days after the first observation 309 days after the first observation 325 days after the first observation No.1 point No.26 point The maximum subsidence 1898 mm (No.18)

Figure 14 .
Figure 14.Comparison of measured data and prediction results based on MMF time function.
open-off cut (m) 52 days after the first observation 84 days after the first observation 114 days after the first observation 178 days after the first observation 52 days after the first observation 84 days after the first observation 114 days after the first observation 178 days after the first observation measured data prediction results

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

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

Table 3 .
Accuracy comparison for prediction results among four time functions.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.57mm and 7.66%, respectively, the mean values of RMSE-Gompertz and RRMSE-Gompertz are 79.77mm and 7.73%, respectively, and the mean values of RMSE-Weibull and RRMSE-Weibull are 90.61mm 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.