A Novel Analytical Model of Mining Subsidence Considering Time Effect Based on the Probability Integral Theory

: Surface subsidence caused by underground coal mining has received wide attention due to its impact on the ecological environment. To obtain ﬁrst-hand data on mining subsidence, we arranged line measurement stations in the mining area of the Weiqiang coalﬁeld and implemented surface movement observation at its ﬁrst mining working face and obtained the dynamic subsidence curve. The subsidence curves reﬂected the initial, active and declining stages of subsidence at each measurement point at a basic level, and were consistent with the general rule of surface subsidence. To accurately analyze the surface movement caused by coal mining, based on the probability integral method and considering the actual advancement progress of the working face, a new three-dimensional dynamic prediction model, namely the strip unit mining model, was proposed. The core assumption of this model was that the longwall mining is regarded as the result of superposition of many units of mining. Considering the inﬂuence of different mining times and periods after mining, the subsidence time function of two-time factors was introduced into the strip unit mining model, and the ﬁnal formulation of the dynamic subsidence prediction of the three-dimensional surface is derived by considering the time factor. Based on the analysis of the measured data of the Weiqiang coalﬁeld, the prediction parameters of surface movement were obtained, and the surface movement and deformation are predicted via the probability integral method. Finally, the measured surface subsidence was compared with the theoretically-predicted one. The good match showed that the parameters of the probability integral method, determined according to the measured data, meet the requirements of geological conditions and mining settings in the mining area, and the predicted subsidence curves agreed well with the ﬁeld measurements, demonstrating the effectiveness of the newly-developed strip unit mining model.


Introduction
The rock is a very complex medium and when the useful minerals are extracted from the rock formations, the original stress equilibrium of the rock around the mining area is disrupted and the stress is redistributed to reach a new equilibrium. In this process, continuous movement, deformation and discontinuous damage (cracking, falling, etc.) of the rock layer and the surface are produced, which is called "mining subsidence" [1]. When the rock mass is affected by various mining methods, the mining subsidence is a very complicated process in terms of both time and space [2].
The surface subsidence caused by underground mining is a global problem. The surface subsidence caused by underground mining attracted the attention of people in the mining field and related fields in the mid-19th century, and gradually formed the of strip units were considered in the process of establishing the model. Firstly, the surface movement caused by mining is studied using the random medium theory. In addition, considering the actual implementation of mining in the working face, the establishment of a strip unit mining model along the entire working face will be closer to the actual situation. It can be regarded as the result of the superposition of many units of body mining, the subsidence time function is introduced into the strip unit mining model [15], and the strip unit mining model with time factor is established. At the same time, based on the observation on the surface movement of the first working face of the Weiqiang coalfield, the measured surface subsidence curve caused by mining and the theoretically predicted curve are compared in order to analyze the selected theoretically predicted model and parameters, and the fitting situation of the actual measurement. The fitting results show that the parameters of the probability integration method, determined according to the measured data, meet the requirements under the geological mining conditions of the mining area, and the predicted results are in positive agreement with the measured results.

Surface Movement Observation in Weiqiang Coalfield
The Weiqiang coalfield is located in the north of the Yuheng mining area and the south area of the Jurassic coalfield in northern Shaanxi. As a result of the loess covering the surface of the mining area being cut, by crisscrossing the loess gullies to form a complex surface shape, the mining damage in the area is much more serious than that under standard conditions. So far, due to special conditions, such as topography, the mining area has not carried out a relatively complete observation and study of the movement law of the mining surface rock layers. The movement law of the surface caused by mining and its strata movement parameters are the basis for studying and solving the "three-under" coal mining and environmental protection problems in the mining area, and are also the key technology to ensure the safety of surface construction facilities [10]. Therefore, carrying out special research on the observation of the law of surface movement and deformation caused by mining in these areas has important theoretical research significance and application value for the rational development and utilization of coal seams, the protection of surface structures in the mining area, the protection of the ecological environment and the sustainable and coordinated development of the economy of the mining area [31]. The 1301 working face is located in the north of the Heng Shan County planning reserve area (wellfield boundary), west of the Baomao Expressway approach road and 600 m east of Majialiang Village. The workings are located in the southern part of the main road, in the west wing of the first pan area, and are the first workings of the mine. It is the first working face of the mine. It is located south of the shaft boundary, west of the unmined area, and east of the coal pillar protection boundary.
The surface elevation of the 1301 working face is 1037.7~1140.3 m, the underground elevation is 789~794.92 m, the recovery area is 634,200 m 2 , and the daily advance speed of the working face is about 8.8 m/d in the early stage of mining.

Overview of the Coal Seam
The average thickness of the coal seam at the 1301 working face is 3.2 m, and the dip angle of the coal seam is <1 • ; the depth of the coal seam is 291-325 m, and the 3-1 coal seam is long-flame coal with a simple structure and stable seam.

Overview of the Roof and Floor
The friable immediate roof of the coal seam is gray silty mudstone, 0.1~0.2 m thick, and the thickness is not constant. It is easy to delaminate and collapse with coal seam during mining, and the rock hardness is relatively soft. The immediate roof is dark gray mediumthick layered silty mudstone, 1.14~7.55 m thick, with argillaceous siltstone interbedded. The main roof is siltstone, medium-grained sandstone, siltstone mudstone and fine-grained sandstone, 17. 43-38.4

m thick.
The floor is siltstone and fine-grained sandstone, 1.62-4.77 m thick. The main floor is light grey medium-thick laminated siltstone with undulating laminations, with obvious contact with the underlying layers, and the layers and thickness are relatively stable.

Coal Mining Technique
The 1301 working face was retrieved using the longwall mining and the roof is managed using the full caving method. The length of the working face in the strike and dip are 2114 m and 300 m, respectively ( Figure 1).
angle of the coal seam is <1°; the depth of the coal seam is 291-325 m, and the 3-1 coal seam is long-flame coal with a simple structure and stable seam.

Overview of the Roof and Floor
The friable immediate roof of the coal seam is gray silty mudstone, 0.1~0.2 m thick, and the thickness is not constant. It is easy to delaminate and collapse with coal seam during mining, and the rock hardness is relatively soft. The immediate roof is dark gray medium-thick layered silty mudstone, 1.14~7.55 m thick, with argillaceous siltstone interbedded. The main roof is siltstone, medium-grained sandstone, siltstone mudstone and fine-grained sandstone, 17.43-38.4 m thick.
The floor is siltstone and fine-grained sandstone, 1.62-4.77 m thick. The main floor is light grey medium-thick laminated siltstone with undulating laminations, with obvious contact with the underlying layers, and the layers and thickness are relatively stable.

Coal Mining Technique
The 1301 working face was retrieved using the longwall mining and the roof is managed using the full caving method. The length of the working face in the strike and dip are 2114 m and 300 m, respectively ( Figure 1).

Layout of Mobile Surface Observation Stations
The comprehensive analysis of the observation area shows that the length of the mining influence radius is 161 m. Considering the randomness of the mining influence and the fact that the sliding range of the ditch slope increases, the mining influence radius is 175 m. In order to enable the observation line to effectively control the impact of mining on the half-basin, the observation station plans to arrange two observation lines: Due to the influence of the terrain, the lengths of the two survey lines were adjusted accordingly when laying the survey lines: the Z line was 500 m long, and the closest distance to the open-off cut was 475 m; the A line was 800 m long, and the distance to the open-off cut was 800 m.

Layout of Mobile Surface Observation Stations
The comprehensive analysis of the observation area shows that the length of the mining influence radius is 161 m. Considering the randomness of the mining influence and the fact that the sliding range of the ditch slope increases, the mining influence radius is 175 m. In order to enable the observation line to effectively control the impact of mining on the half-basin, the observation station plans to arrange two observation lines:

1.
Arranging a strike observation line (Z line) along the strike direction at the center of the first mining face. The strike observation line is perpendicular to the working face open-off cut, and the length of LZ should try to control the half-basin range.

2.
An observation line (A line) is arranged along the inclination, and the observation line is perpendicular to the strike observation line (Z line) and intersects with the strike observation line. According to requirements, the dip observation line shall be able to control the full basin along the inclination.
Due to the influence of the terrain, the lengths of the two survey lines were adjusted accordingly when laying the survey lines: the Z line was 500 m long, and the closest distance to the open-off cut was 475 m; the A line was 800 m long, and the distance to the open-off cut was 800 m.
According to the buried depth of the working face of the Weiqiang coalfield, the distance between measuring points is 25 m. According to the actual shape of the surface, the final survey line points are arranged as follows: the A Line has a total of 33 survey points, and 29 survey points are set up in the early stage, numbered A1~A29, then add 4 measuring points, numbered JA1~JA4; the Z line has a total of 20 measuring points, Z1~Z20. The layout of the actual measuring points is shown in Figure 2.
According to the buried depth of the working face of the Weiqiang coalfield, the distance between measuring points is 25 m. According to the actual shape of the surface, the final survey line points are arranged as follows: the A Line has a total of 33 survey points, and 29 survey points are set up in the early stage, numbered A1~A29, then add 4 measuring points, numbered JA1~JA4; the Z line has a total of 20 measuring points, Z1~Z20. The layout of the actual measuring points is shown in Figure 2.

Implementation of Surface Mobile Observation
The internal collation method of the continuous measurement results of the observatory is the same as the conventional method. Finally, the plane coordinates and elevation of the control point of the observatory need to be obtained. In the process of connecting and measuring the working face of the Weiqiang coalfield, RTK measuring technology [32] was used to measure the coordinates and the elevation of control points.
To accurately determine the spatial position of the working survey point before the start of the surface movement, two comprehensive observations were independently carried out after the connection measurement and before the start of the surface movement. The content of the comprehensive observation includes the determination of the plane position and the elevation of each measuring point, as well as the distance between each measuring point, the appropriate increase in level measurement work between the first and last comprehensive observations, and the time interval of repeated level measurement, depending on the surface subsidence rate.
By November 2016, a total of 2 comprehensive observations and 27 daily observations had been carried out in the observation area of the first mining face of the Weiqiang coalfield, and a total of 1508 sets of observation data had been obtained.

Calculation of Observed Results
In order to ensure the correctness of the observation results, the field observation results should be checked again before the collation of the internal results, and then the calculation and adjustment of various correction data should be carried out.
The internal calculation of observation data is primarily used to calculate the elevation of each measuring point and the horizontal distance between two adjacent measuring points, followed by calculating the movement deformation value and subsidence speed value of each measuring point. Attached below is the internal calculation data of some measuring points, as shown in Table 1.

Implementation of Surface Mobile Observation
The internal collation method of the continuous measurement results of the observatory is the same as the conventional method. Finally, the plane coordinates and elevation of the control point of the observatory need to be obtained. In the process of connecting and measuring the working face of the Weiqiang coalfield, RTK measuring technology [32] was used to measure the coordinates and the elevation of control points.
To accurately determine the spatial position of the working survey point before the start of the surface movement, two comprehensive observations were independently carried out after the connection measurement and before the start of the surface movement. The content of the comprehensive observation includes the determination of the plane position and the elevation of each measuring point, as well as the distance between each measuring point, the appropriate increase in level measurement work between the first and last comprehensive observations, and the time interval of repeated level measurement, depending on the surface subsidence rate.
By November 2016, a total of 2 comprehensive observations and 27 daily observations had been carried out in the observation area of the first mining face of the Weiqiang coalfield, and a total of 1508 sets of observation data had been obtained.

Calculation of Observed Results
In order to ensure the correctness of the observation results, the field observation results should be checked again before the collation of the internal results, and then the calculation and adjustment of various correction data should be carried out.
The internal calculation of observation data is primarily used to calculate the elevation of each measuring point and the horizontal distance between two adjacent measuring points, followed by calculating the movement deformation value and subsidence speed value of each measuring point. Attached below is the internal calculation data of some measuring points, as shown in Table 1. After that, the subsidence calculation of each measuring point is carried out, and the subsidence of n points in m times of observation is: where w n is the subsidence value of n point, mm; H 0 n and H m n are the elevation of point n at the first and m observations, respectively, mm.
As shown in Figure 3, the trend subsidence curve illustrates that the subsidence value of the middle measurement point is relatively smaller than that of the two end measurement points after critical mining, primarily because the slope of the middle measurement point is small, the slip amount caused by the slope is small, and therefore the subsidence of the measuring point in the middle gentle area is slightly smaller.
where w n is the subsidence value of n point, mm; H n 0 and H n m are the elevation of point n at the first and m observations, respectively, mm. As shown in Figure 3, the trend subsidence curve illustrates that the subsidence value of the middle measurement point is relatively smaller than that of the two end measurement points after critical mining, primarily because the slope of the middle measurement point is small, the slip amount caused by the slope is small, and therefore the subsidence of the measuring point in the middle gentle area is slightly smaller. It can be seen from the subsidence curve of the dip observation line in Figure 4 that the change of the surface subsidence basin is relatively gentle near the edge, while the change of its value increases near the goaf directly above, which indicates that the surface subsidence and slope degree of the goaf boundary is small. Affected by the larger buried depth, the subsidence curve has good continuity and a smooth transition, which is also a feature of surface subsidence in general deep-buried coal mines. From the comparative analysis of the relationship between the subsidence value and the time change in the above figure, it can be seen that the amount of surface subsidence was not obvious in the previous observations. For the strike observation line, the subsidence of each measuring point is mainly divided into the subsidence initiation stage, the subsidence active stage, and the subsidence decline stage. Among all the survey points so far, the last major subsidence occurred at the Z19 survey point, and the maximum subsidence value was 1.460 The strike measuring points It can be seen from the subsidence curve of the dip observation line in Figure 4 that the change of the surface subsidence basin is relatively gentle near the edge, while the change of its value increases near the goaf directly above, which indicates that the surface subsidence and slope degree of the goaf boundary is small. Affected by the larger buried depth, the subsidence curve has good continuity and a smooth transition, which is also a feature of surface subsidence in general deep-buried coal mines. From the comparative analysis of the relationship between the subsidence value and the time change in the above figure, it can be seen that the amount of surface subsidence was not obvious in the previous observations. For the strike observation line, the subsidence of each measuring point is mainly divided into the subsidence initiation stage, the subsidence active stage, and the subsidence decline stage. Among all the survey points so far, the last major subsidence occurred at the Z19 survey point, and the maximum subsidence value was 1.460 m. From the last few subsidence curves of the Z line, we can see the subsidence development trend of the goaf outside the observation area; that is, the maximum subsidence should be between 1.3 and 1.5 m.
Sustainability 2022, 14, x FOR PEER REVIEW 7 of 20 m. From the last few subsidence curves of the Z line, we can see the subsidence development trend of the goaf outside the observation area; that is, the maximum subsidence should be between 1.3 and 1.5 m. The point Z13 on the strike observation line is observed, and three stages can be observed from the start of the movement of the surface to the stop of the movement: the initiation stage, the active stage and the decay stage. The maximum sinking rate at this point is reached only after the working surface has pushed through the point for a certain distance. Table 2 shows that the sinking rate at this point is less than 10 mm/d during the start and decay phases and between 10 mm/d and 90 mm/d during the active phase, with a maximum sinking rate of 67.5 mm/d. It can also be concluded from the data in the table that the amount of sinking at this point during the active phase exceeds 80% of the total amount of sinking, thus indicating that the sinking of the ground surface primarily occurs during the sinking. This shows that the surface subsidence mainly occurs in the active period of subsidence. The change in the subsidence speed and the cumulative subsidence value of measuring point Z13 with time is shown in Figure 5.  The point Z13 on the strike observation line is observed, and three stages can be observed from the start of the movement of the surface to the stop of the movement: the initiation stage, the active stage and the decay stage. The maximum sinking rate at this point is reached only after the working surface has pushed through the point for a certain distance. Table 2 shows that the sinking rate at this point is less than 10 mm/d during the start and decay phases and between 10 mm/d and 90 mm/d during the active phase, with a maximum sinking rate of 67.5 mm/d. It can also be concluded from the data in the table that the amount of sinking at this point during the active phase exceeds 80% of the total amount of sinking, thus indicating that the sinking of the ground surface primarily occurs during the sinking. This shows that the surface subsidence mainly occurs in the active period of subsidence. The change in the subsidence speed and the cumulative subsidence value of measuring point Z13 with time is shown in Figure 5. and decay phases and between 10 mm/d and 90 mm/d during the active phase, with a maximum sinking rate of 67.5 mm/d. It can also be concluded from the data in the table that the amount of sinking at this point during the active phase exceeds 80% of the total amount of sinking, thus indicating that the sinking of the ground surface primarily occurs during the sinking. This shows that the surface subsidence mainly occurs in the active period of subsidence. The change in the subsidence speed and the cumulative subsidence value of measuring point Z13 with time is shown in Figure 5. From the overall analysis, the subsidence curves of the two observation lines reflect the initial stage, active stage, and recession stage of the subsidence of each measuring point on the surface, which is in line with the general law of surface subsidence.

Mining Subsidence Prediction Theory
For one or more working faces that are intended to be mined, the work of precalculating the surface movement and deformation affected by the mining, according to its geological coal mining conditions and the selected predicted functions and parameters, is called "surface movement mining subsidence prediction". The mining subsidence prediction is one of the core contents of mining subsidence and protection, and it is of great significance to the theoretical research and production practice of mining subsidence [4].
According to the different ways of establishing prediction methods, they can be divided into three types: empirical method, theoretical model method and influence function method. The most widely used methods in China are the probability integration method, typical curve method, profile function method, similar material simulation and numerical calculation method [33,34].

Probability Integral Method
The probability integral method regards the rock mass as a random medium. According to the random medium theory, the subsidence basin of the surface unit caused by the unit mining is normally distributed and is consistent with the distribution of the probability density [9]. Therefore, the subsidence profile equation caused by the whole mining can be expressed as the integral formula of the probability density function. I surface movement and deformation caused by the whole mining are calculated by the superposition principle [35,36].
The random particle medium model is shown in Figure 6a [37]. It takes the assumption that these particles are small spheres of the same size and uniform mass and are packed in squares of the same size and uniform arrangement. When the ball in square a1, in the first layer, is removed, one of the balls in squares b1 and b2 in the second layer will roll into this square due to gravity, and the probability that they roll in both are 1/2; if the ball in square b1 rolls into a1, then b1 will be occupied by the ball from square c1 or c2 in the third layer. Similarly, if the ball in square b2 rolls into a1, then b2 will be occupied by the ball rolling from the c2 or c3 square in the third layer. According to the multiplication and addition theorem in probability, the probability that the ball in the a1 grid is removed and the c1, c2 and c3 grids are emptied are 1/4, 1/2 and 1/4, respectively; the probability of emptying caused by the removal of a1 is written in the corresponding grid, and its probability distribution histogram can be formed, as shown in Figure 6b. If the size of the grid is very small, the histogram will approximate a smooth normal distribution probability density curve [38].
ball in square b1 rolls into a1, then b1 will be occupied by the ball from square c1 or c2 in the third layer. Similarly, if the ball in square b2 rolls into a1, then b2 will be occupied by the ball rolling from the c2 or c3 square in the third layer. According to the multiplication and addition theorem in probability, the probability that the ball in the a1 grid is removed and the c1, c2 and c3 grids are emptied are 1/4, 1/2 and 1/4, respectively; the probability of emptying caused by the removal of a1 is written in the corresponding grid, and its probability distribution histogram can be formed, as shown in Figure 6b. If the size of the grid is very small, the histogram will approximate a smooth normal distribution probability density curve [38]. Figure 6. Random granular medium model.
Since the probability integration method is based on the random medium theory, the rock formation movement is regarded as a random process that obeys statistical laws. From a statistical point of view, the entire mining area is decomposed into the mining of infinite tiny unit mining. The impact of the surface is equal to the sum of the impact of each unit of mining on the rock formation and the surface. The subsidence profile equation caused by mining can be expressed as the integral formula of the probability density function. On the basis of the subsidence profile equation, other movement and deformation equations of the measuring points are derived. Therefore, it can be said that the subsidence profile equation is the basis of the probability integral method [35].

Knothe Time Function Model
From the statistical point of view, in order to infer the formation process of the unit basin, it can be considered that the unit mining is completed in an instant, but when the unit body is excavated, a final unit subsidence basin is not immediately formed on the surface; a certain period of time and process has to be experienced. In order to describe the formation process of the subsidence basin of the unit, it can be represented by a subsidence time function.
In 1952, in order to describe the complex time-space process of surface subsidence caused by underground mining with mathematical methods, Knothe made the following assumptions: the instantaneous subsidence speed v(t) of the surface point is proportional to the difference between the instantaneous subsidence value w(t) and the final subsidence value w max of the point, the difference between the subsidence and the final subsidence value is proportional, and the following differential equation is established to describe the subsidence time course of the surface point [16,39]: Since the probability integration method is based on the random medium theory, the rock formation movement is regarded as a random process that obeys statistical laws. From a statistical point of view, the entire mining area is decomposed into the mining of infinite tiny unit mining. The impact of the surface is equal to the sum of the impact of each unit of mining on the rock formation and the surface. The subsidence profile equation caused by mining can be expressed as the integral formula of the probability density function. On the basis of the subsidence profile equation, other movement and deformation equations of the measuring points are derived. Therefore, it can be said that the subsidence profile equation is the basis of the probability integral method [35].

Knothe Time Function Model
From the statistical point of view, in order to infer the formation process of the unit basin, it can be considered that the unit mining is completed in an instant, but when the unit body is excavated, a final unit subsidence basin is not immediately formed on the surface; a certain period of time and process has to be experienced. In order to describe the formation process of the subsidence basin of the unit, it can be represented by a subsidence time function.
In 1952, in order to describe the complex time-space process of surface subsidence caused by underground mining with mathematical methods, Knothe made the following assumptions: the instantaneous subsidence speed v(t) of the surface point is proportional to the difference between the instantaneous subsidence value w(t) and the final subsidence value w max of the point, the difference between the subsidence and the final subsidence value is proportional, and the following differential equation is established to describe the subsidence time course of the surface point [16,39]: where c is the time coefficient; w max is the maximum subsidence value when the surface subsidence is stable, mm; w(t) is the instantaneous subsidence value at time t, mm. Substituting the initial condition t = 0 and w(0) = 0, the solution of Equation (2) can be obtained: Equation (3) is the final expression of the Knothe time function. The Knothe time function is the earliest "theoretical model" to study the dynamic subsidence law of the ground surface and is considered to be able to approximately describe the whole process of the subsidence of the ground point [16] (Figure 7).
Equation (3) is the final expression of the Knothe time function. The Knothe time function is the earliest "theoretical model" to study the dynamic subsidence law of the ground surface and is considered to be able to approximately describe the whole process of the subsidence of the ground point [16] (Figure 7).

Establishment of the Strip Unit Mining Model
Coal seam mining movement is always a dynamic process in space, so it is necessary to establish a three-dimensional dynamic prediction model, in order to be able to intuitively and accurately analyze the surface movement caused by coal seam mining. Based on this, this paper has made a bold attempt on the establishment of the three-dimensional dynamic prediction model of coal seam mining. The whole modeling idea mainly has the following points: (1) At present, among the theories on the surface movement caused by mining, the stochastic medium theory is the most mature and widely used. It has been recognized and applied by many scholars for its rigorous reasoning and high accuracy conclusions. It is completely feasible to study the surface movement caused by mining with stochastic medium theory. (2) There are two speeds in the advancing process of the working face: one is the speed of mining a knife along the working face; the other is the speed of the whole working face. These two speeds are different and related, but in terms of mining, the first speed is generally relatively fast, and it can be considered that during this process, the overlying rock layer and the surface have not been affected. Therefore, the first speed can be ignored and only the effect of the second velocity is considered. (3) The mining of the working face is not carried out unit by unit, but along the working face. It can be imagined that if the strip unit mining model along the entire working face is established, it will be closer to the actual situation. The strip unit mining of the working face can be regarded as the result of the superposition of many unit body mining [40]. (4) After the mining of each working face strip unit, the surface does not immediately form the final strip subsidence basin, but must go through a certain amount of time and process. The subsidence time function is introduced into the strip unit mining model to establish a strip unit mining model with time factor [41].

Establishment of the Strip Unit Mining Model
Coal seam mining movement is always a dynamic process in space, so it is necessary to establish a three-dimensional dynamic prediction model, in order to be able to intuitively and accurately analyze the surface movement caused by coal seam mining. Based on this, this paper has made a bold attempt on the establishment of the three-dimensional dynamic prediction model of coal seam mining. The whole modeling idea mainly has the following points: (1) At present, among the theories on the surface movement caused by mining, the stochastic medium theory is the most mature and widely used. It has been recognized and applied by many scholars for its rigorous reasoning and high accuracy conclusions. It is completely feasible to study the surface movement caused by mining with stochastic medium theory. (2) There are two speeds in the advancing process of the working face: one is the speed of mining a knife along the working face; the other is the speed of the whole working face. These two speeds are different and related, but in terms of mining, the first speed is generally relatively fast, and it can be considered that during this process, the overlying rock layer and the surface have not been affected. Therefore, the first speed can be ignored and only the effect of the second velocity is considered. (3) The mining of the working face is not carried out unit by unit, but along the working face. It can be imagined that if the strip unit mining model along the entire working face is established, it will be closer to the actual situation. The strip unit mining of the working face can be regarded as the result of the superposition of many unit body mining [40]. (4) After the mining of each working face strip unit, the surface does not immediately form the final strip subsidence basin, but must go through a certain amount of time and process. The subsidence time function is introduced into the strip unit mining model to establish a strip unit mining model with time factor [41]. (5) The mining time of each strip unit and the elapsed time after mining are inconsistent.
These two times are both different and related. Considering the influence of the "double time" factors, it can better reflect the subsidence of the surface. It can realize the dynamic three-dimensional prediction of the surface movement caused by coal seam mining. (6) On the time issue of the dynamic movement of the surface, the unit is generally "day (d)". Although the subsidence process is a continuous time function, if it is accurate to "hours" or "seconds", it will inevitably bring many mathematical problems to the model. Therefore, in the prediction of surface movement, the time step can be defined as day; that is, the strip unit mined per day is a unit subsidence factor The text continues here.

Strip Unit Model without Time Factor
First, considering strip unit mining per unit thickness; assuming that the incline length of the working face is D, the average buried depth of the mining coal seam is H, and selecting projection point of intersection point of strike axis line and coal wall plane in surface water plane as coordinate origin point, it can be obtained: Using MATLAB software (Matlab Version 2015a under the license of education version of Dalian University of Technology), the result of the above formula can be obtained as: where: er f (x) = 2 √ π x 0 e −t 2 dt is a built-in function in MATLAB. D is the inclination length of the working face, m; r is the major influence radius, m. Equation (4) is the expected subsidence basin expression for the mining of the strip unit per unit thickness of the working face, without considering the time factor. Since the thickness of the mining coal seam is generally small, usually within 10 m, this thickness is very small compared with the buried depth of the coal seam. Therefore, when establishing the model of the working face inclined to the full-height strip unit mining, the full-height strip unit mining of the working face is regarded as the superposition of the unit thickness strip unit mining; that is, the thickness of a mined coal seam is multiplied on the basis of the original expression, as shown in the following formula: where m is the thickness of the mined coal seam, m. Considering the fragmentation of the rock and soil in the subsidence basin of the strip unit, the space of the goaf will not be completely filled; that is, the maximum subsidence of the surface will not exceed the mining height of the coal seam m. If the subsidence coefficient of the subsidence basin of the full-height strip unit be η e , then Equation (6) is transformed into: where η e is the subsidence coefficient. Equation (7) is the final model for full-height strip unit mining, without considering the time factor. However, it is worth noting that in the influence range of the goaf and the filling rock blocks in the goaf, the force of each spatial position is also different, which makes the bulking coefficient of each strip unit different after mining; that is, the subsidence coefficient η e of each strip unit is different. However, the geological structure is relatively simple, and each rock layer is relatively horizontal mining conditions, and the subsidence coefficient η e can be approximately considered equal. After each strip unit of the working face is excavated, the surface does not immediately form the final subsidence basin of the unit. Generally, it will go through several steps of "start subsidence → maximum subsidence speed → subsidence speed gradually decreases → subsidence stabilizes". This process can be described by the dynamic subsidence time function of the surface.

Strip Unit Model Considering Time
The classic Knothe subsidence time function model is used here. The characteristic of this model is that the subsidence speed reaches the maximum at the beginning of the subsidence, which is similar to the failure law of the overlying rock and soil layer of the shallow coal seam; that is, during the first subsidence, the surface had a sudden change in subsidence, and the stratum was completely fractured and collapsed within a short period of time. During the subsequent surface movement, the subsidence speed gradually decreased and, finally, stabilized. Therefore, choosing the Knothe time function can better explain the surface movement law of shallow coal seams.
Equation (3) suggests that the time function of the time factor T experienced by strip unit mining is: where T is the mining time of the strip unit, d.
Combining Equations (6) and (7), we can obtain the subsidence prediction expression considering the time factor T of the strip unit mining experience:

Considering the Time Factor of Full-Thickness Strip Unit Advancement, t
Analyzing the actual mining process of the working face, it can be seen that the time experienced by each strip unit after mining is inconsistent with the advancing time of the working face. It is reasonable and necessary to take the mining speed into consideration in the strip unit mining model. The x axis is parallel to the advancing direction of the working face, and the w axis is perpendicular to the working face. The advancing speed of the working face along the positive direction of the x axis is v (m/d), and the advancing distance over time t is S. Subsequently, the dynamic subsidence basin within the mining distance is obtained from the integral of Equation (9): Since s = vt, ds = vdt, Equation (9) can be transformed into: where v is the speed of advancing along the positive direction of the x-axis, m/d. The time experienced by each strip unit after mining is different. The strip unit mined first has a relatively long time, and the strip unit mined later has a relatively short time. If one strip unit is used as a factor, it will inevitably bring a lot of inconvenience to the calculation. Here, the time "day (d)" is used as the unit: the strip mined every day is used as a calculation factor, and the time experienced by each factor is considered. The surface subsidence is regarded as the result of the superposition of various strip factors after different mining times. Then the dynamic subsidence model of strip unit factor mining with k = (1, 2, 3, · · · · · · ) on the k day is determined by:

The final Three-Dimensional Dynamic Subsidence Model of the Surface
For underground mining, the surface begins to subside only when the mining reaches a certain distance. The mining distance when the surface begins to subside is called the starting length, which is represented by d. Here, the corresponding time may be called the starting time, which is represented by T d . Then, from the above assumptions, we have T d = d v . Therefore, the following improvements can be made to the Knothe subsidence time function: where d is the starting length, m.
In the strip mining within the starting length d, the surface does not subside; that is, the strip mining units within the mining distance range d are synchronized in time, so the strip unit mining of this distance segment can be regarded as a whole, that is: It can be calculated by MATLAB: Substitute Equation (14) into Equation (15), we obtain: When the advancing distance of the working face exceeds the starting length, the ground surface begins to subside. The strip unit mining formula for each day thereafter is as follows: Calculating based on the day when the surface moved, the surface subsidence on the k day after the start length d should be the superposition of the previous surface subsidence and the subsidence of the strip unit on each subsequent day, which can be determined by Equation (18): Equation (18) is the final expression of the three-dimensional surface dynamic subsidence prediction, and the expression of the dynamic predicted subsidence velocity can be obtained by partial differential of time variable of Equation (18).

Determination of Prediction Parameters of Surface Movement
When using the probability integral method to predict the movement and deformation of the mining surface, five predicted parameters should first be determined: subsidence coefficient η, tangent value of main influence angle tanβ, horizontal movement coefficient b, inflection point offset d * and propagation angle of extraction θ. Through the analysis of the measured data, these parameters can be obtained by calculation [42].
The most important step for obtaining parameters is to determine the maximum subsidence value. Under the condition of critical mining (ensuring that full subsidence is achieved in both strike and dip directions, it can be called critical mining), the maximum subsidence value of each survey line is obtained through the observation data of mobile observation station.
Under the condition of full mining of nearly horizontal coal seam, the maximum subsidence value, maximum horizontal displacement and maximum inclination value in strike and dip are approximately equal. In this paper, the same value is taken in strike and dip for calculation.
(1) Subsidence coefficient η The subsidence coefficient is determined by the ratio of the maximum subsidence value on the surface to the average mining thickness under critical mining conditions: where α is the dip angle of the coal seam. The calculation results are shown in the Table 3. (2) The tangent value of main influence angle tan β The angle between the boundary point connecting the main influence range and the mining boundary line and the horizontal line is called the tangent value of main influence angle [43], and the calculation formula is shown in Equation (20): From r =H/tanβ, r = ω max /i max , we have: where r is the main influence radius, m; H is the average buried depth of coal seam, mm; ω max is the maximum subsidence value of the surface, mm; i max is the maximum tilt value of the observation point, mm/m. The calculated results are shown in Table 4. (3) Propagation angle of extraction θ.
The propagation angle of extraction represents the offset degree of the surface subsidence curve to the downhill direction, under the condition of inclined coal seam mining. The value is the angle between the inflection point of the subsidence curve and the calculation boundary of the same side goaf (considering the offset distance of the inflection point) and the horizontal line on the downhill side. Its value is primarily related to the coal seam dip angle and lithology [43]. In general, the influence propagation angle of extraction is equal to the maximum subsidence angle. It is calculated by: where c is the coefficient related to lithology, α is the dip angle of the coal seam, • . The dip angle of the working face studied in this study is 0.5 • , therefore, the propagation angle of extraction can be approximately 90 • . The horizontal movement coefficient is a coefficient that reflects the relationship between the maximum horizontal movement value and the maximum subsidence value [45]. Under critical/supercritical mining conditions, it is calculated by: where u max is the maximum value of the horizontal movement of the surface. According to the observation data of the two survey lines, the horizontal movement factor can be calculated (Table 5). Several groups of surface movement parameters in different observation periods are selected below (Table 6). It is clear that when the goaf range is small, the fluctuation of each parameter is large, indicating that the change is unstable. With the expansion of the mining range, the parameters gradually become stable, indicating that the surface movement at this time has also stabilized. From this result, it can be inferred that with the continuous increase of the goaf, the deformation of the ground surface will not change much, and the range of surface subsidence behind the goaf will not change much. The last few measurements show that the subsidence has stabilized.

Performance of the Probability Integral Theory
A large number of measured data have confirmed that the surface movement in mountainous areas has certain rules to follow. However, the surface deformation caused by mining subsidence has a positive relationship with the topography and landforms, and the regularity is very poor under general conditions. Figures 8 and 9 are the predicted parameters of the surface probability integral, based on the measured surface movement and deformation, comparing the measured surface subsidence curve caused by mining with the theoretical predicted curve, in order to analyze the selected theoretical prediction model and parameters and the measured fitting.

Performance of the Probability Integral Theory
A large number of measured data have confirmed that the surface movement in mountainous areas has certain rules to follow. However, the surface deformation caused by mining subsidence has a positive relationship with the topography and landforms, and the regularity is very poor under general conditions. Figures 8 and 9 are the predicted parameters of the surface probability integral, based on the measured surface movement and deformation, comparing the measured surface subsidence curve caused by mining with the theoretical predicted curve, in order to analyze the selected theoretical prediction model and parameters and the measured fitting.
Since the entry measurement work started after 540 m of mining at the working face, in order to better reflect the law of surface movement in the early stage of mining, the Knothe subsidence time function model was used to numerically calculate the surface movement in the early stage of mining. Figure 8 shows that the surface activity near the open-off cut is active at the initial stage of mining, and the active subsidence period of each measuring point changes with the advancement of the working face. The theoretical prediction of maximum surface subsidence reaches approximately 1.4 m, which is basically consistent with the measured value. The subsidence to the strike observation line gradually expands forward with the advancement of the working face. On the strike main section, the surface points subsidence sequentially, from near to far. When the working face advances to 400 m, the goaf will have a certain range. When the working face advances to 400 m, the subsidence of the surface point within a certain range of the goaf is stable. As the working face continues to advance, the surface points in front of the working face cycle the previous subsidence process in turn, which is consistent with the actual measurement.  Figure 9 shows that after the critical mining is achieved, the measured values are greatly affected by the terrain, and the values of local points jump. But on the whole, the subsidence trend of the theoretical value and the measured value is basically consistent. The maximum subsidence value of the surface is basically consistent with the theoretical prediction value, reaching approximately 1.4 m. Therefore, the measured surface subsidence value of the dip observation line in the influence range of the working face has a good fit with the theoretically predicted curve. The subsidence values (mm)  Since the entry measurement work started after 540 m of mining at the working face, in order to better reflect the law of surface movement in the early stage of mining, the Knothe subsidence time function model was used to numerically calculate the surface movement in the early stage of mining. Figure 8 shows that the surface activity near the open-off cut is active at the initial stage of mining, and the active subsidence period of each measuring point changes with the advancement of the working face. The theoretical prediction of maximum surface subsidence reaches approximately 1.4 m, which is basically consistent with the measured value. The subsidence to the strike observation line gradually expands forward with the advancement of the working face. On the strike main section, the surface points subsidence sequentially, from near to far. When the working face advances to 400 m, the goaf will have a certain range. When the working face advances to 400 m, the subsidence of the surface point within a certain range of the goaf is stable. As the working face continues to advance, the surface points in front of the working face cycle the previous subsidence process in turn, which is consistent with the actual measurement.
Sustainability 2022, 14, x FOR PEER REVIEW 18 of 20 Figure 9. Comparison between measured surface subsidence curve and theoretical prediction curve of dip observation line A.

Conclusions
This study conducted surface movement observations at the first mining face in the Weiqiang coalfield, obtained the dynamic subsidence curves, and proposed a three-dimensional dynamic prediction model to obtain the predicted deformation of surface movement, and compares the measured values of subsidence with the predicted values. The conclusions of this paper are as follows.
(1) According to the mining environment of the first working face of the Weiqiang coalfield, the design of surface movement observation was carried out. A total of 2 systematic observations and 27 daily observations were carried out in the observation area, and 1508 sets of reliable observation data were obtained. Through data analysis, it can be seen that the surface movement of each measuring point in the observation area has experienced the start-up stage, the active stage, and the decline stage. The duration of the three stages varies at each measuring point, and the last few observations show that the surface movement has basically stabilized. (2) A new three-dimensional dynamic prediction model is proposed: the strip unit model. Considering the influence of dual time factors, the probability integration method is used to predict the deformation of the main section of the surface mobile basin during semi-infinite mining. The estimated parameters suitable for the probability integration method are obtained. (3) The measured surface subsidence curve is fitted by the theoretical prediction formula. The fitting results show that the parameters of the probability integration method determined according to the measured data meet the requirements of the geological mining conditions in the mining area, and the predicted results are in positive agreement with the measured results.   Figure 9 shows that after the critical mining is achieved, the measured values are greatly affected by the terrain, and the values of local points jump. But on the whole, the subsidence trend of the theoretical value and the measured value is basically consistent. The maximum subsidence value of the surface is basically consistent with the theoretical prediction value, reaching approximately 1.4 m. Therefore, the measured surface subsidence value of the dip observation line in the influence range of the working face has a good fit with the theoretically predicted curve.

Conclusions
This study conducted surface movement observations at the first mining face in the Weiqiang coalfield, obtained the dynamic subsidence curves, and proposed a threedimensional dynamic prediction model to obtain the predicted deformation of surface movement, and compares the measured values of subsidence with the predicted values. The conclusions of this paper are as follows.
(1) According to the mining environment of the first working face of the Weiqiang coalfield, the design of surface movement observation was carried out. A total of 2 systematic observations and 27 daily observations were carried out in the observation area, and 1508 sets of reliable observation data were obtained. Through data analysis, it can be seen that the surface movement of each measuring point in the observation area has experienced the start-up stage, the active stage, and the decline stage. The duration of the three stages varies at each measuring point, and the last few observations show that the surface movement has basically stabilized. (2) A new three-dimensional dynamic prediction model is proposed: the strip unit model.
Considering the influence of dual time factors, the probability integration method is used to predict the deformation of the main section of the surface mobile basin during semi-infinite mining. The estimated parameters suitable for the probability integration method are obtained. (3) The measured surface subsidence curve is fitted by the theoretical prediction formula.
The fitting results show that the parameters of the probability integration method determined according to the measured data meet the requirements of the geological mining conditions in the mining area, and the predicted results are in positive agreement with the measured results.