# Analysis of Mining Subsidence Characteristics and Deformation Prediction considering Size Parameters and Mechanical Parameters

Authors: Shangguo Liu; Kaixin Li; Weiping Shi; Zijun Wang; Haoran Zhang; Zhaoying Li
Journal: Geofluids
 (2022)

---

## Abstract

Crushing and caving of roof strata in goaf is an important factor affecting surface subsidence. In order to study the influence law of mining under thick loose layer, overburden size parameters, and initial mechanical parameters of collapse zone on surface subsidence, the law of surface subsidence was analyzed through measured data. On this basis, with the help of FLAC3D numerical simulation software, the surface subsidence was simulated by changing the above factors, and the evolution characteristics of the maximum surface subsidence under the influence of various factors were analyzed. Based on the grey relational analysis, the relational degree between each factor and the maximum surface subsidence value is analyzed, and through multiple nonlinear regression, the prediction model of mining subsidence deformation under thick loose layer considering the influence of seven factors is established. The results show that the mining thickness has the greatest influence on the surface subsidence, and cohesion has the least influence on surface subsidence. The rationality of the subsidence prediction model is verified by an example.


---

## Body

## 1. Introduction

Underground coal mining destroys the original mechanical balance of rock mass, and the overlying strata move and deform. When the mining area reaches a certain range, the movement and deformation will extend to the surface [1]. In order to protect the roadway, ground facilities, and water body from the harmful effects of mining, it is necessary to master the law of surface movement and deformation affected by mining and establish a reasonable prediction model of surface subsidence. For the analysis of the law of surface subsidence in coal mining, it is mainly to set up the surface movement observation station in the mining area to obtain the measured data, supplemented by physical material similarity simulation and three-dimensional numerical simulation. Many scholars have obtained rich research results on this basis. Wang et al. [2, 3] analyzed the apparent characteristics of surface subsidence under the mining conditions of thick loose layers on the basis of measurement data and summarized the intrinsic mechanism of movement and damage produced by thick loose layers and rock. Zhang et al. [4] studied the law of surface subsidence in deep mining through field measurement. Tajduś et al. [5] analyzed the ground movement and deformation in longwall face mining based on the measured data of the ground observation station in the mining area. Xu et al. [6] analyzed the law of surface movement and deformation when high-intensity mining was adopted through the measured data of surface movement observation station. Wang et al. [7] solved the rock movement parameters and dynamic surface movement parameters of inclined coal seam mining by systematically analyzing the measured data of observation stations and summarized the corresponding laws and characteristics of surface movement and deformation. Guo et al. [8] based on the measured data analyzed the appearance of surface dynamic movement characteristics under the mining condition of thick collapsible loess layer, studied the surface failure characteristics caused by mining under this geological condition, and defined what is the advanced crack distance and angle. With the rapid development of science and technology, through the continuous improvement of measuring instruments, the monitoring technology of surface subsidence continues to develop, and the monitoring efficiency is obviously improved. For example, 3D laser scanning technology [9–12], InSAR/D-InSAR [13–17], UAV deformation monitoring [18–21], and other methods have become very important surface subsidence monitoring methods.Considering the convenience, simplicity, and rapidity of numerical simulation and similar material simulation, these methods are mostly used to analyze the law of surface subsidence. Xu et al. [22] simulated and analyzed the overlying strata movement law and fissure development in fully mechanized mining face of “three soft” coal seam with the help of similar materials. Dai et al. [23] put forward the anti-inclined layered upward mining method to solve the problem of severe surface subsidence in mining extrathick steep coal seam and confirmed the mechanism of controlling strata movement by numerical simulation and similar material experiments. Zhang et al. [24] based on the geological and mining conditions of a working face in Zhaogu No.2 Coal Mine made a research by using similar material simulation experiment method and revealed the subsidence mechanism of bedrock and loose layer. Chen et al. [25] used the discrete element numerical software (UDEC) to study the law of surface movement and deformation under the condition of mining thick loose layer. Zhang et al. [26] used FLAC3D numerical experiments to study surface subsidence and ranked the importance of factors influencing surface subsidence in strip mining by designing orthogonal experiments. Wang et al. [27] rely on 3DEC numerical simulation software to analyze the sensitivity of the main control factors of mining subsidence.In surface subsidence prediction, the main methods that are widely used are the probability integral method and the profile function method. The probability integral method is a prediction method developed by Liu [28] on the basis of random medium theory. Relevant scholars have constantly revised and improved the probability integral method [29–33]. The profile function method is based on the measured typical subsidence curve characteristics, selects the appropriate fitting function, determines the constant of the function or equation, and then, applies it to the deformation prediction of the main section of the surface movement, which mainly includes hyperbolic function, error function, exponential function, and trigonometric function [34]. In addition, Chen et al. [16] used the combination of D-InSAR technology and SVR algorithm to monitor and dynamically predict the mining subsidence. Zha et al. [35] inversed the predicted parameters of probability integral method by genetic algorithm. Li et al. [36] established a parameter calculation model of probability integral method based on support vector machine, which improved the prediction accuracy. Li et al. [37] combined GIS technology and probability integral method to predict mining subsidence.Based on the geological conditions of 1301 working face in Lilou Coal Mine, on the basis of studying and analyzing the law of surface subsidence in Lilou Coal Mine, this paper analyzes the law of its influence on surface subsidence from the perspectives of the size parameters of overlying bedrock and loose layer (thickness of coal, thickness of loose layer, and thickness of bedrock) and the initial mechanical parameters of collapse zone (elastic modulus, tensile strength, cohesion, and internal friction angle) and establishes a surface subsidence prediction model to realize surface subsidence prediction before coal seam mining, and provides theoretical basis for surface subsidence control measures in advance.


## 2. Analysis of Measured Surface Subsidence of 1301 Working Face

### 2.1. Overview of Mine

The LiLou Coal Mine is located in Yuncheng County, Shandong Province, and its mine field is located at the northernmost part of Juye coal field, adjacent to Guotun mine field and Pengzhuang mine field in the south and southeast, respectively. The mine geological structure is complex and of medium type. The 1301 working face mainly mines 3 coal, which is located in the middle and lower part of Shanxi Group, with an average distance of 53.2 m from the bottom of the three ash. According to the borehole histogram at the site, the average coal thickness is 6.8 m, and the coal seam dip angle is 7°~13°. The strike length of the working face is 2265.8 m, and the inclined length is 223.4 m. The longwall fully mechanized top coal caving mining method is adopted, the backward mining is adopted, and the roof is managed by all caving method.


### 2.2. Analysis of Measured Data of Surface Subsidence

In order to fully reflect the law of surface subsidence in 1301 working face of Lilou coal mine, four observation lines were laid along the main section of ore body strike and dip. The observation line Z is arranged along the strike direction of the working face, and the observation line H, line N, and line B are arranged along the dip direction of the working face. The relative position relationship between the working face and the observation lines is shown in Figure1. Considering all kinds of constraints, this analysis is mainly based on the measured data along the observation line Z. There are 121 measuring points in the observation line Z, numbered Z01-Z121. The first observation was made in December 2016, and 30 observations were made in November 2019. In terms of data acquisition, the plane coordinates of control points are measured by GPS static relative positioning method, and the elevation is measured by bubble level, all of which are carried out according to relevant requirements.Figure 1
Relative position diagram of measuring line and working face.By analyzing the monitoring data, the surface subsidence curves of each measurement point on the Z line during the advancement of 1301 working face at different times were compiled, as shown in Figure2. As can be seen from Figure 2, as the working face continues to advance, the surface subsidence of each measuring point on the observation line Z gradually increases, and the maximum subsidence point gradually moves to the center of the goaf. The increasing process of surface subsidence conforms to the general law of mining subsidence. The maximum surface subsidence observed at the last time was 2966 mm, located at Z56; the subsidence coefficient is 0.44. On the basis of the last phase of surface subsidence data in Figure 2, the missing data are supplemented by interpolation method to obtain the final surface subsidence inclination curve as shown in Figure 3, and the surface curvature curve is shown in Figure 4. It can be seen from Figure 3 that the inclination value of the observation line is roughly antisymmetric about the center of the goaf. The maximum inclination value of the surface at the side of the open-off cut is 4.859 mm/m, which is located at Z38; the maximum inclination value of the ground surface on one side of the stop line is 5.183 mm/m, which is located at Z81. As can be seen from Figure 4, the maximum positive curvature value of the ground surface at the end of the mining observation line is 0.11 mm/m2, located at point Z84, and the maximum negative curvature value is -0.12 mm/m2, located at point Z80. Overall, the surface curvature value is small, the curvature deformation is not significant, and the curvature curve does not show obvious features.Figure 2
Strike observation subsidence curve of 1301 working face.Figure 3
Surface inclination curve of 1301 working face towards the main section.Figure 4
Surface curvature curve of 1301 working face towards main section.


## 2.1. Overview of Mine

The LiLou Coal Mine is located in Yuncheng County, Shandong Province, and its mine field is located at the northernmost part of Juye coal field, adjacent to Guotun mine field and Pengzhuang mine field in the south and southeast, respectively. The mine geological structure is complex and of medium type. The 1301 working face mainly mines 3 coal, which is located in the middle and lower part of Shanxi Group, with an average distance of 53.2 m from the bottom of the three ash. According to the borehole histogram at the site, the average coal thickness is 6.8 m, and the coal seam dip angle is 7°~13°. The strike length of the working face is 2265.8 m, and the inclined length is 223.4 m. The longwall fully mechanized top coal caving mining method is adopted, the backward mining is adopted, and the roof is managed by all caving method.


## 2.2. Analysis of Measured Data of Surface Subsidence

In order to fully reflect the law of surface subsidence in 1301 working face of Lilou coal mine, four observation lines were laid along the main section of ore body strike and dip. The observation line Z is arranged along the strike direction of the working face, and the observation line H, line N, and line B are arranged along the dip direction of the working face. The relative position relationship between the working face and the observation lines is shown in Figure1. Considering all kinds of constraints, this analysis is mainly based on the measured data along the observation line Z. There are 121 measuring points in the observation line Z, numbered Z01-Z121. The first observation was made in December 2016, and 30 observations were made in November 2019. In terms of data acquisition, the plane coordinates of control points are measured by GPS static relative positioning method, and the elevation is measured by bubble level, all of which are carried out according to relevant requirements.Figure 1
Relative position diagram of measuring line and working face.By analyzing the monitoring data, the surface subsidence curves of each measurement point on the Z line during the advancement of 1301 working face at different times were compiled, as shown in Figure2. As can be seen from Figure 2, as the working face continues to advance, the surface subsidence of each measuring point on the observation line Z gradually increases, and the maximum subsidence point gradually moves to the center of the goaf. The increasing process of surface subsidence conforms to the general law of mining subsidence. The maximum surface subsidence observed at the last time was 2966 mm, located at Z56; the subsidence coefficient is 0.44. On the basis of the last phase of surface subsidence data in Figure 2, the missing data are supplemented by interpolation method to obtain the final surface subsidence inclination curve as shown in Figure 3, and the surface curvature curve is shown in Figure 4. It can be seen from Figure 3 that the inclination value of the observation line is roughly antisymmetric about the center of the goaf. The maximum inclination value of the surface at the side of the open-off cut is 4.859 mm/m, which is located at Z38; the maximum inclination value of the ground surface on one side of the stop line is 5.183 mm/m, which is located at Z81. As can be seen from Figure 4, the maximum positive curvature value of the ground surface at the end of the mining observation line is 0.11 mm/m2, located at point Z84, and the maximum negative curvature value is -0.12 mm/m2, located at point Z80. Overall, the surface curvature value is small, the curvature deformation is not significant, and the curvature curve does not show obvious features.Figure 2
Strike observation subsidence curve of 1301 working face.Figure 3
Surface inclination curve of 1301 working face towards the main section.Figure 4
Surface curvature curve of 1301 working face towards main section.


## 3. Theoretical Analysis of Collapse Zone Height in 1301 Working Face

After coal seam mining, the collapse of overlying rock mass is the main cause of surface subsidence. The change of the mechanical properties of rock mass within the collapse zone will affect the height of the collapse zone, and it will also cause the change of the surface subsidence value. According to the borehole histogram of 1301 working face, the roof strata of 3 coal in 1301 working face mainly contain sandstone and sandy mudstone, among which the experimental value of sandstone compressive strength is 54~77.2 MPa and platts hardness coefficientf=6.3. The experimental values of compressive strength of argillaceous sandstone and sandy mudstone are 20~42 MPa, platts hardness coefficient f=3.2, and the roof rock generally belongs to medium hard rock mass. According to the formula for calculating the height of collapse zone in layered mining of medium-thick coal seam in “Regulations for setting coal pillars and pressing coal mining in building, water body, railway and main roadway” [38], it can be known that the height of collapse zone in medium-hard roof is:
(1)Hk=100∑M4.7∑M+19±2.2,where ∑M is cumulative mining thickness, m.According to the mine data and formula (1), the height of roof collapse zone in 1301 working face is within the range of 11.14~15.54 m.


## 4. Numerical Simulation of Surface Subsidence Law considering Size and Mechanical Parameters

### 4.1. Establishment of Numerical Model

Taking the 1301 working face in Lilou Coal Mine as the research object, referring to the comprehensive histogram of the working face, and considering the thickness and lithology of each rock layer between the surface and the coal seam, the numerical model is established. The working face has a strike length of 2265.8 m and a dip length of 223.4 m, while the final model design size was2400m∗1200m∗960m (length∗width∗height) in order to eliminate the influence of boundary effects. The numerical model is shown in Figure 5. A reasonable constitutive model of rock and soil mass can reflect the main deformation characteristics of rock and soil materials. The main reason for adopting Mohr-Coulomb model in this simulation is that the failure forms of overburden are mainly tensile shear and compressive shear failure in the process of underground coal seam excavation. Mohr-Coulomb model criterion can consider both shear yield and tensile failure, so this model is applicable to the shear-tensile failure characteristics of rock mass. The left, right, and bottom boundaries of the model are fixed, and the top is set as a free boundary. Meanwhile, based on the research content, it can be known from the above that the maximum height of the collapse zone of 1301 working face can reach 15.54 m Therefore, the rock strata within the collapse zone above the coal seam are merged, and the FLAC3D modeling parameters of each rock stratum are shown in Table 1. In the middle of the top of the model, survey lines are laid along the strike direction to monitor the surface subsidence.Figure 5
Numerical model diagram.Table 1
Rock mechanics parameters.
NameThickness/mDensity/kg∗m-3Modulus of elasticity/GPaTensile strength/MPaCohesion/MPaInternal friction angle/°Loose layer86618000.140.010.0217Mudstone25.9524305.63.84.129Fine sandstone6.2526009.85.45.235Siltstone2.08275011.43.53.632Mudstone13.9024305.63.84.129Fine sandstone4.9726009.85.45.235Fine sandstone15.5426009.85.45.235Coal6.814001.161.11.2525Fine sandstone3.1526009.85.45.235Mudstone4.7724305.63.84.129Fine sandstone10.5926009.85.45.235


### 4.2. Analysis of the Influence of Size Parameters on Surface Subsidence

In order to analyze the influence law of size parameters on the surface subsidence caused by coal seam mining, based on the numerical simulation model of 1301 working face in Lilou coal mine, the influence law of different thickness of loose layer, thickness of bedrock, and thickness of coal on the surface subsidence was studied. In order to make the simulation results of various schemes comparable, when one influencing factor is analyzed, other influencing factors remain unchanged.


#### 4.2.1. Influence of Thickness of Loose Layer on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of loose layer were drawn, as shown in Figure6. From the simulation results, it can be seen that when the thickness of loose layer gradually increases from 0 m to 866 m, the maximum surface subsidence values are 3184 mm, 3319 mm, 3477 mm, 3796 mm, 3543 mm, 3398 mm, 3109 mm, and 2965 mm, respectively, and the corresponding subsidence coefficients are 0.47, 0.49, 0.51, 0.56, 0.52, 0.50, 0.46, and 0.44, respectively. It can be clearly seen from this, under the condition of constant mining thickness and bedrock thickness, with the increase of loose layer thickness, the maximum subsidence value of the surface changes nonlinearly. When the thickness of loose layer changes from 0 m to 250 m, the maximum surface subsidence gradually increases from 3184 mm to 3796 mm. When the thickness of loose layer changes from 250 m to 866 m, the maximum surface subsidence gradually decreases from 3796 mm to 2965 mm. The amount of change of the subsidence coefficient is smaller, but overall, with the increase of the thickness of loose layer, it shows roughly the same evolution trend as the maximum subsidence value of the surface.Figure 6
Surface subsidence curve in strike direction with different thickness of loose layer.As the thickness of the loose layer increases, the surface subsidence caused by coal mining does not change in a linear relationship, which means, it is not the case that the thicker the loose layer is, the greater the surface subsidence is. The main reason for this is the loose characteristics of the loose layer itself. When the loose layer is in a critical thickness and the influence of mining spreads to the loose layer from the bottom up, the consolidation and compression of the soil continues to occur in the process of its downward movement due to the loose characteristics of the loose layer structure, and the value of surface subsidence increases continuously as a result. However, when the thickness of loose layer continues to increase, the load effect produced by loose layer is less than its absorption effect. Under the same thickness, the load produced by loose layer is smaller than that of rock layer, and the damage degree of bedrock under its influence is relatively small. At the same time, the lower loose layer digests and absorbs the damage of overburden, which makes it difficult for its influence to reach the upper loose layer, resulting in the gradual decrease of land subsidence.The fitted curves of the evolutionary relationship between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer are shown in Figure7. The fitting expression of the function between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer was obtained as:
(2)W=4.699H−0.011H2+6.774×10−6H3+3149.806R2=0.938,q=6.783h−1.673×10−6h2+9.895×10−10h3+0.465R2=0.993.Figure 7
Fitting curve of maximum subsidence/subsidence coefficient and thickness of loose layer.
(a)
Maximum subsidence(b)
Subsidence coefficient


#### 4.2.2. Influence of Thickness of Bedrock on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of bedrock were drawn, as shown in Figure8.Figure 8
Surface subsidence curve in strike direction with different thickness of bedrock.It can be seen from the simulation results that when the thickness of bedrock gradually increases from 50 m to 120 m, the maximum subsidence are 6734 mm, 5249 mm, 2965 mm, 2256 mm, 1415 mm, 1083 mm, 762 mm, and 399 mm, respectively, and the corresponding subsidence coefficients are 0.99, 0.77, 0.44, 0.33, 0.21, 0.16, 0.11, and 0.06, respectively. From this, it is clear that under the condition of constant thickness of coal and thickness of loose layer, the maximum subsidence is negatively correlated with thickness of bedrock, which means, with the increase of thickness of bedrock, the maximum subsidence gradually decreases. When the thickness of bedrock increases from 50 m to 90 m, the maximum subsidence decreases by 5319 mm (from 6734 mm to 1415 mm), with a large decline. When the thickness of bedrock gradually increases from 90 m to 120 m, the maximum ground settlement decreases slowly (decreases by 1016 mm). With the increase of thickness of bedrock, the subsidence coefficient presents the same evolution characteristics as the maximum subsidence. The main reason is that when the bedrock is thin, it cannot form its own equilibrium structure inside the bedrock, so under the influence of mining and the load of the overlying thick loose layer, the displacement and bending deformation of the bedrock will be more obvious, and the transmission to the surface will cause large surface subsidence, while, when the bedrock is thick, a more stable load-bearing structure will be formed in the bedrock, which can withstand the load of the overlying loose layer, and the breakage is not easy to occur. Therefore, the surface subsidence will be correspondingly reduced.The fitting curve of the evolution relationship between the maximum subsidence, subsidence coefficient, and thickness of bedrock is shown in Figure9. The function expressions of maximum subsidence, subsidence coefficient, and thickness of bedrock obtained by fitting are as follows:
(3)W=−346.567D+1.530D2+20183.363R2=0.981,q=−0.051D+2.245×10−4D2+2.964R2=0.982.Figure 9
Fitting curve of maximum subsidence/subsidence coefficient and thickness of bedrock.
(a)
Maximum subsidence(b)
Subsidence coefficient


#### 4.2.3. Influence of Thickness of Coal on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of coal were drawn, as shown in Figure10.Figure 10
Surface subsidence curve in strike direction with different thickness of coal.It can be seen from the simulation results that when the thickness of coal gradually increases from 1.8 m to 8.8 m, the maximum subsidence values are 120 mm, 275 mm, 696 mm, 1337 mm, 1997 mm, 2965 mm, 4005 mm, and 5189 mm, respectively, and the corresponding subsidence coefficients are 0.07, 0.10, 0.18, 0.28, 0.34, 0.44, 0.51, and 0.59, respectively. Therefore, the maximum subsidence is proportional to the thickness of coal, that is, with the increase of the thickness of coal, the maximum subsidence increases linearly under the same other conditions. This is mainly because the larger the thickness of coal, the larger the space of mined-out area, the larger the range of collapse zone, and the more severe the movement and destruction of overlying strata. When it is transmitted to the surface, the surface subsidence also increases correspondingly. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and thickness of coal evolution is shown in Figure11. The function expressions of the maximum subsidence, subsidence coefficient, and thickness of coal obtained by fitting are as follows:
(4)W=730.976h−1809.173R2=0.953,q=0.078h−0.098R2=0.993.Figure 11
Fitting curve of maximum subsidence/subsidence coefficient and thickness of coal.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.3. Analysis of the Influence of Mechanical Parameters on Surface Subsidence

In order to analyze the influence law of mechanical parameters of collapse zone on surface subsidence in coal seam mining, based on the above numerical simulation model of 1301 working face in Lilou coal mine, the influence law of elastic modulus, tensile strength, cohesion, and internal friction angle of different collapse zones on surface subsidence was studied. For making the simulation results of various schemes comparable, when one influencing factor is analyzed, other influencing factors remain unchanged.


#### 4.3.1. Influence of Elastic Modulus of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different elastic modulus were drawn, as shown in Figure12.Figure 12
Surface subsidence curve of different elastic modulus in strike direction.It can be seen from the simulation results that when the elastic modulus gradually increases from 4.8 GPa to 39.8 GPa, the maximum subsidence are 3116 mm, 3023 mm, 2965 mm, 2899 mm, 2852 mm, 2798 mm, and 2778 mm, respectively, and the corresponding subsidence coefficients are 0.46, 0.44, 0.44, 0.43, 0.42, 0.41, 0.41, and 0.41, respectively. From this, it is clear that, with the increase of the elastic modulus of the collapse zone, the maximum subsidence of the ground surface gradually decreases under the same other conditions. However, compared with the influence of other parameters, the total change of the maximum subsidence caused by it is not significant. When the elastic modulus increases from 4.8 GPa to 39.8 GPa, the maximum subsidence decreases by only 338 mm. Therefore, the change of elastic modulus has no obvious influence on the subsidence coefficient. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient and elastic modulus evolution is shown in Figure13. The function expressions of maximum subsidence, subsidence coefficient, and elastic modulus obtained by fitting are:
(5)W=3.423E−0.056R2=0.986,q=−0.0014E+0.459R2=0.912.Figure 13
Fitting curve of maximum subsidence/subsidence coefficient and elastic modulus.
(a)
Maximum subsidence(b)
Subsidence coefficient


#### 4.3.2. Influence of Tensile Strength of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different tensile strength were drawn, as shown in Figure14.Figure 14
Surface subsidence curve of different tensile strength in strike direction.It can be seen from the simulation results that when the tensile strength gradually increases from 0.4 MPa to 15.4 MPa, the maximum subsidence values are 3077 mm, 3050 mm, 3022 mm, 3000 mm, 2968 mm, 2900 mm, 2891 mm, and 2891 mm, respectively, and the corresponding subsidence coefficients are 0.45, 0.45, 0.44, 0.44, 0.44, 0.43, 0.43, and 0.43, respectively. From this, it is clear that, overall, under the condition of other conditions being constant, the maximum subsidence gradually decreases with the increase of tensile strength. But it has obvious segmentation characteristics. When the tensile strength increases from 0.4 MPa to 5.4 MPa, the surface subsidence decreases from 3077 mm to 2900 mm. However, when the tensile strength increases from 5.4 MPa to 15.4 MPa, the maximum subsidence basically does not change with the increase of tensile strength. Therefore, there is a critical value of the tensile strength. When the tensile strength changes within a range less than the critical value, the maximum subsidence decreases with its increase. When the tensile strength exceeds the critical value, the increase of tensile strength has almost no effect on the maximum subsidence. The change of the maximum subsidence is relatively small, so the change of subsidence coefficient is not significant. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and tensile strength evolution is shown in Figure15. The function expressions of the maximum subsidence, subsidence coefficient, and tensile strength obtained by fitting are as follows:
(6)W=237.025∗exp−T4.316+2873.316R2=0.932,q=−6.429T+0.456R2=0.890.Figure 15
Fitting curve of maximum subsidence/subsidence coefficient and tensile strength.
(a)
Maximum subsidence(b)
Subsidence coefficient


#### 4.3.3. Influence of Cohesion of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different cohesion were drawn, as shown in Figure16.Figure 16
Surface subsidence curve of different cohesion in strike direction.It can be seen from the simulation results that when the cohesion gradually increases from 1.2 MPa to 30.2 MPa, the maximum subsidence are 3702 mm, 3373 mm, 2965 mm, 2609 mm, 2287 mm, 2251 mm, and 2251 mm, respectively, and the corresponding subsidence coefficients are 0.54, 0.50, 0.44, 0.38, 0.34, 0.33, 0.33, and 0.33, respectively. As a whole, the maximum subsidence and subsidence coefficient gradually decrease with the increase of cohesion. However, it has obvious segmentation characteristics. When the cohesion increases from 1.2 MPa to 15.2 MPa, the surface subsidence value and coefficient decrease quickly, from 3702 mm to 2287 mm; however, when the cohesion increases from 15.2 MPa to 30.2 MPa, the subsidence value and subsidence coefficient decrease slowly. Therefore, there is a critical value of cohesion, and when the cohesion changes within a range less than the critical value, the maximum subsidence and subsidence coefficient change greatly. When the cohesion exceeds the critical value, the maximum subsidence and subsidence coefficient will no longer change obviously with the increase of cohesion. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and cohesion evolution is shown in Figure17. The function expressions of maximum subsidence, subsidence coefficient, and cohesion obtained by fitting are as follows:
(7)W=1854.91∗exp−c6.21+2201.23R2=0.993,q=0.270∗exp−c6.543+0.322R2=0.993.Figure 17
Fitting curve of maximum subsidence/subsidence coefficient and cohesion.
(a)
Maximum subsidence(b)
Subsidence coefficient


#### 4.3.4. Influence of Internal Friction Angle of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different internal friction angles were drawn, as shown in Figure18.Figure 18
Surface subsidence curve of different internal friction angles in strike direction.It can be seen from the simulation results that when the internal friction angle gradually increases from 10° to 45°, the maximum subsidence values are 3938 mm, 3690 mm, 3486 mm, 3314 mm, 3118 mm, 2965 mm, 2826 mm, and 2688 mm, respectively, and the corresponding subsidence coefficients are 0.58, 0.54, and 0.51, respectively. It can be seen that, with the increase of the internal friction angle in the collapse zone, the maximum subsidence and subsidence coefficient basically show a linear decreasing trend under the same other conditions. The fitting curve of the relationship between the maximum surface subsidence value and subsidence coefficient and the evolution of internal friction angle is shown in Figure19. The function expressions of maximum subsidence, subsidence coefficient, and internal friction angle obtained by fitting are:
(8)W=−35.31φ+4224.07R2=0.991,q=−0.005φ+0.618R2=0.986.Figure 19
Fitting curve of maximum subsidence/subsidence coefficient and internal friction angle.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.1. Establishment of Numerical Model

Taking the 1301 working face in Lilou Coal Mine as the research object, referring to the comprehensive histogram of the working face, and considering the thickness and lithology of each rock layer between the surface and the coal seam, the numerical model is established. The working face has a strike length of 2265.8 m and a dip length of 223.4 m, while the final model design size was2400m∗1200m∗960m (length∗width∗height) in order to eliminate the influence of boundary effects. The numerical model is shown in Figure 5. A reasonable constitutive model of rock and soil mass can reflect the main deformation characteristics of rock and soil materials. The main reason for adopting Mohr-Coulomb model in this simulation is that the failure forms of overburden are mainly tensile shear and compressive shear failure in the process of underground coal seam excavation. Mohr-Coulomb model criterion can consider both shear yield and tensile failure, so this model is applicable to the shear-tensile failure characteristics of rock mass. The left, right, and bottom boundaries of the model are fixed, and the top is set as a free boundary. Meanwhile, based on the research content, it can be known from the above that the maximum height of the collapse zone of 1301 working face can reach 15.54 m Therefore, the rock strata within the collapse zone above the coal seam are merged, and the FLAC3D modeling parameters of each rock stratum are shown in Table 1. In the middle of the top of the model, survey lines are laid along the strike direction to monitor the surface subsidence.Figure 5
Numerical model diagram.Table 1
Rock mechanics parameters.
NameThickness/mDensity/kg∗m-3Modulus of elasticity/GPaTensile strength/MPaCohesion/MPaInternal friction angle/°Loose layer86618000.140.010.0217Mudstone25.9524305.63.84.129Fine sandstone6.2526009.85.45.235Siltstone2.08275011.43.53.632Mudstone13.9024305.63.84.129Fine sandstone4.9726009.85.45.235Fine sandstone15.5426009.85.45.235Coal6.814001.161.11.2525Fine sandstone3.1526009.85.45.235Mudstone4.7724305.63.84.129Fine sandstone10.5926009.85.45.235


## 4.2. Analysis of the Influence of Size Parameters on Surface Subsidence

In order to analyze the influence law of size parameters on the surface subsidence caused by coal seam mining, based on the numerical simulation model of 1301 working face in Lilou coal mine, the influence law of different thickness of loose layer, thickness of bedrock, and thickness of coal on the surface subsidence was studied. In order to make the simulation results of various schemes comparable, when one influencing factor is analyzed, other influencing factors remain unchanged.


### 4.2.1. Influence of Thickness of Loose Layer on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of loose layer were drawn, as shown in Figure6. From the simulation results, it can be seen that when the thickness of loose layer gradually increases from 0 m to 866 m, the maximum surface subsidence values are 3184 mm, 3319 mm, 3477 mm, 3796 mm, 3543 mm, 3398 mm, 3109 mm, and 2965 mm, respectively, and the corresponding subsidence coefficients are 0.47, 0.49, 0.51, 0.56, 0.52, 0.50, 0.46, and 0.44, respectively. It can be clearly seen from this, under the condition of constant mining thickness and bedrock thickness, with the increase of loose layer thickness, the maximum subsidence value of the surface changes nonlinearly. When the thickness of loose layer changes from 0 m to 250 m, the maximum surface subsidence gradually increases from 3184 mm to 3796 mm. When the thickness of loose layer changes from 250 m to 866 m, the maximum surface subsidence gradually decreases from 3796 mm to 2965 mm. The amount of change of the subsidence coefficient is smaller, but overall, with the increase of the thickness of loose layer, it shows roughly the same evolution trend as the maximum subsidence value of the surface.Figure 6
Surface subsidence curve in strike direction with different thickness of loose layer.As the thickness of the loose layer increases, the surface subsidence caused by coal mining does not change in a linear relationship, which means, it is not the case that the thicker the loose layer is, the greater the surface subsidence is. The main reason for this is the loose characteristics of the loose layer itself. When the loose layer is in a critical thickness and the influence of mining spreads to the loose layer from the bottom up, the consolidation and compression of the soil continues to occur in the process of its downward movement due to the loose characteristics of the loose layer structure, and the value of surface subsidence increases continuously as a result. However, when the thickness of loose layer continues to increase, the load effect produced by loose layer is less than its absorption effect. Under the same thickness, the load produced by loose layer is smaller than that of rock layer, and the damage degree of bedrock under its influence is relatively small. At the same time, the lower loose layer digests and absorbs the damage of overburden, which makes it difficult for its influence to reach the upper loose layer, resulting in the gradual decrease of land subsidence.The fitted curves of the evolutionary relationship between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer are shown in Figure7. The fitting expression of the function between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer was obtained as:
(2)W=4.699H−0.011H2+6.774×10−6H3+3149.806R2=0.938,q=6.783h−1.673×10−6h2+9.895×10−10h3+0.465R2=0.993.Figure 7
Fitting curve of maximum subsidence/subsidence coefficient and thickness of loose layer.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.2.2. Influence of Thickness of Bedrock on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of bedrock were drawn, as shown in Figure8.Figure 8
Surface subsidence curve in strike direction with different thickness of bedrock.It can be seen from the simulation results that when the thickness of bedrock gradually increases from 50 m to 120 m, the maximum subsidence are 6734 mm, 5249 mm, 2965 mm, 2256 mm, 1415 mm, 1083 mm, 762 mm, and 399 mm, respectively, and the corresponding subsidence coefficients are 0.99, 0.77, 0.44, 0.33, 0.21, 0.16, 0.11, and 0.06, respectively. From this, it is clear that under the condition of constant thickness of coal and thickness of loose layer, the maximum subsidence is negatively correlated with thickness of bedrock, which means, with the increase of thickness of bedrock, the maximum subsidence gradually decreases. When the thickness of bedrock increases from 50 m to 90 m, the maximum subsidence decreases by 5319 mm (from 6734 mm to 1415 mm), with a large decline. When the thickness of bedrock gradually increases from 90 m to 120 m, the maximum ground settlement decreases slowly (decreases by 1016 mm). With the increase of thickness of bedrock, the subsidence coefficient presents the same evolution characteristics as the maximum subsidence. The main reason is that when the bedrock is thin, it cannot form its own equilibrium structure inside the bedrock, so under the influence of mining and the load of the overlying thick loose layer, the displacement and bending deformation of the bedrock will be more obvious, and the transmission to the surface will cause large surface subsidence, while, when the bedrock is thick, a more stable load-bearing structure will be formed in the bedrock, which can withstand the load of the overlying loose layer, and the breakage is not easy to occur. Therefore, the surface subsidence will be correspondingly reduced.The fitting curve of the evolution relationship between the maximum subsidence, subsidence coefficient, and thickness of bedrock is shown in Figure9. The function expressions of maximum subsidence, subsidence coefficient, and thickness of bedrock obtained by fitting are as follows:
(3)W=−346.567D+1.530D2+20183.363R2=0.981,q=−0.051D+2.245×10−4D2+2.964R2=0.982.Figure 9
Fitting curve of maximum subsidence/subsidence coefficient and thickness of bedrock.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.2.3. Influence of Thickness of Coal on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of coal were drawn, as shown in Figure10.Figure 10
Surface subsidence curve in strike direction with different thickness of coal.It can be seen from the simulation results that when the thickness of coal gradually increases from 1.8 m to 8.8 m, the maximum subsidence values are 120 mm, 275 mm, 696 mm, 1337 mm, 1997 mm, 2965 mm, 4005 mm, and 5189 mm, respectively, and the corresponding subsidence coefficients are 0.07, 0.10, 0.18, 0.28, 0.34, 0.44, 0.51, and 0.59, respectively. Therefore, the maximum subsidence is proportional to the thickness of coal, that is, with the increase of the thickness of coal, the maximum subsidence increases linearly under the same other conditions. This is mainly because the larger the thickness of coal, the larger the space of mined-out area, the larger the range of collapse zone, and the more severe the movement and destruction of overlying strata. When it is transmitted to the surface, the surface subsidence also increases correspondingly. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and thickness of coal evolution is shown in Figure11. The function expressions of the maximum subsidence, subsidence coefficient, and thickness of coal obtained by fitting are as follows:
(4)W=730.976h−1809.173R2=0.953,q=0.078h−0.098R2=0.993.Figure 11
Fitting curve of maximum subsidence/subsidence coefficient and thickness of coal.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.2.1. Influence of Thickness of Loose Layer on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of loose layer were drawn, as shown in Figure6. From the simulation results, it can be seen that when the thickness of loose layer gradually increases from 0 m to 866 m, the maximum surface subsidence values are 3184 mm, 3319 mm, 3477 mm, 3796 mm, 3543 mm, 3398 mm, 3109 mm, and 2965 mm, respectively, and the corresponding subsidence coefficients are 0.47, 0.49, 0.51, 0.56, 0.52, 0.50, 0.46, and 0.44, respectively. It can be clearly seen from this, under the condition of constant mining thickness and bedrock thickness, with the increase of loose layer thickness, the maximum subsidence value of the surface changes nonlinearly. When the thickness of loose layer changes from 0 m to 250 m, the maximum surface subsidence gradually increases from 3184 mm to 3796 mm. When the thickness of loose layer changes from 250 m to 866 m, the maximum surface subsidence gradually decreases from 3796 mm to 2965 mm. The amount of change of the subsidence coefficient is smaller, but overall, with the increase of the thickness of loose layer, it shows roughly the same evolution trend as the maximum subsidence value of the surface.Figure 6
Surface subsidence curve in strike direction with different thickness of loose layer.As the thickness of the loose layer increases, the surface subsidence caused by coal mining does not change in a linear relationship, which means, it is not the case that the thicker the loose layer is, the greater the surface subsidence is. The main reason for this is the loose characteristics of the loose layer itself. When the loose layer is in a critical thickness and the influence of mining spreads to the loose layer from the bottom up, the consolidation and compression of the soil continues to occur in the process of its downward movement due to the loose characteristics of the loose layer structure, and the value of surface subsidence increases continuously as a result. However, when the thickness of loose layer continues to increase, the load effect produced by loose layer is less than its absorption effect. Under the same thickness, the load produced by loose layer is smaller than that of rock layer, and the damage degree of bedrock under its influence is relatively small. At the same time, the lower loose layer digests and absorbs the damage of overburden, which makes it difficult for its influence to reach the upper loose layer, resulting in the gradual decrease of land subsidence.The fitted curves of the evolutionary relationship between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer are shown in Figure7. The fitting expression of the function between the maximum subsidence, the subsidence coefficient, and the thickness of loose layer was obtained as:
(2)W=4.699H−0.011H2+6.774×10−6H3+3149.806R2=0.938,q=6.783h−1.673×10−6h2+9.895×10−10h3+0.465R2=0.993.Figure 7
Fitting curve of maximum subsidence/subsidence coefficient and thickness of loose layer.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.2.2. Influence of Thickness of Bedrock on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of bedrock were drawn, as shown in Figure8.Figure 8
Surface subsidence curve in strike direction with different thickness of bedrock.It can be seen from the simulation results that when the thickness of bedrock gradually increases from 50 m to 120 m, the maximum subsidence are 6734 mm, 5249 mm, 2965 mm, 2256 mm, 1415 mm, 1083 mm, 762 mm, and 399 mm, respectively, and the corresponding subsidence coefficients are 0.99, 0.77, 0.44, 0.33, 0.21, 0.16, 0.11, and 0.06, respectively. From this, it is clear that under the condition of constant thickness of coal and thickness of loose layer, the maximum subsidence is negatively correlated with thickness of bedrock, which means, with the increase of thickness of bedrock, the maximum subsidence gradually decreases. When the thickness of bedrock increases from 50 m to 90 m, the maximum subsidence decreases by 5319 mm (from 6734 mm to 1415 mm), with a large decline. When the thickness of bedrock gradually increases from 90 m to 120 m, the maximum ground settlement decreases slowly (decreases by 1016 mm). With the increase of thickness of bedrock, the subsidence coefficient presents the same evolution characteristics as the maximum subsidence. The main reason is that when the bedrock is thin, it cannot form its own equilibrium structure inside the bedrock, so under the influence of mining and the load of the overlying thick loose layer, the displacement and bending deformation of the bedrock will be more obvious, and the transmission to the surface will cause large surface subsidence, while, when the bedrock is thick, a more stable load-bearing structure will be formed in the bedrock, which can withstand the load of the overlying loose layer, and the breakage is not easy to occur. Therefore, the surface subsidence will be correspondingly reduced.The fitting curve of the evolution relationship between the maximum subsidence, subsidence coefficient, and thickness of bedrock is shown in Figure9. The function expressions of maximum subsidence, subsidence coefficient, and thickness of bedrock obtained by fitting are as follows:
(3)W=−346.567D+1.530D2+20183.363R2=0.981,q=−0.051D+2.245×10−4D2+2.964R2=0.982.Figure 9
Fitting curve of maximum subsidence/subsidence coefficient and thickness of bedrock.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.2.3. Influence of Thickness of Coal on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different thicknesses of coal were drawn, as shown in Figure10.Figure 10
Surface subsidence curve in strike direction with different thickness of coal.It can be seen from the simulation results that when the thickness of coal gradually increases from 1.8 m to 8.8 m, the maximum subsidence values are 120 mm, 275 mm, 696 mm, 1337 mm, 1997 mm, 2965 mm, 4005 mm, and 5189 mm, respectively, and the corresponding subsidence coefficients are 0.07, 0.10, 0.18, 0.28, 0.34, 0.44, 0.51, and 0.59, respectively. Therefore, the maximum subsidence is proportional to the thickness of coal, that is, with the increase of the thickness of coal, the maximum subsidence increases linearly under the same other conditions. This is mainly because the larger the thickness of coal, the larger the space of mined-out area, the larger the range of collapse zone, and the more severe the movement and destruction of overlying strata. When it is transmitted to the surface, the surface subsidence also increases correspondingly. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and thickness of coal evolution is shown in Figure11. The function expressions of the maximum subsidence, subsidence coefficient, and thickness of coal obtained by fitting are as follows:
(4)W=730.976h−1809.173R2=0.953,q=0.078h−0.098R2=0.993.Figure 11
Fitting curve of maximum subsidence/subsidence coefficient and thickness of coal.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.3. Analysis of the Influence of Mechanical Parameters on Surface Subsidence

In order to analyze the influence law of mechanical parameters of collapse zone on surface subsidence in coal seam mining, based on the above numerical simulation model of 1301 working face in Lilou coal mine, the influence law of elastic modulus, tensile strength, cohesion, and internal friction angle of different collapse zones on surface subsidence was studied. For making the simulation results of various schemes comparable, when one influencing factor is analyzed, other influencing factors remain unchanged.


### 4.3.1. Influence of Elastic Modulus of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different elastic modulus were drawn, as shown in Figure12.Figure 12
Surface subsidence curve of different elastic modulus in strike direction.It can be seen from the simulation results that when the elastic modulus gradually increases from 4.8 GPa to 39.8 GPa, the maximum subsidence are 3116 mm, 3023 mm, 2965 mm, 2899 mm, 2852 mm, 2798 mm, and 2778 mm, respectively, and the corresponding subsidence coefficients are 0.46, 0.44, 0.44, 0.43, 0.42, 0.41, 0.41, and 0.41, respectively. From this, it is clear that, with the increase of the elastic modulus of the collapse zone, the maximum subsidence of the ground surface gradually decreases under the same other conditions. However, compared with the influence of other parameters, the total change of the maximum subsidence caused by it is not significant. When the elastic modulus increases from 4.8 GPa to 39.8 GPa, the maximum subsidence decreases by only 338 mm. Therefore, the change of elastic modulus has no obvious influence on the subsidence coefficient. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient and elastic modulus evolution is shown in Figure13. The function expressions of maximum subsidence, subsidence coefficient, and elastic modulus obtained by fitting are:
(5)W=3.423E−0.056R2=0.986,q=−0.0014E+0.459R2=0.912.Figure 13
Fitting curve of maximum subsidence/subsidence coefficient and elastic modulus.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.3.2. Influence of Tensile Strength of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different tensile strength were drawn, as shown in Figure14.Figure 14
Surface subsidence curve of different tensile strength in strike direction.It can be seen from the simulation results that when the tensile strength gradually increases from 0.4 MPa to 15.4 MPa, the maximum subsidence values are 3077 mm, 3050 mm, 3022 mm, 3000 mm, 2968 mm, 2900 mm, 2891 mm, and 2891 mm, respectively, and the corresponding subsidence coefficients are 0.45, 0.45, 0.44, 0.44, 0.44, 0.43, 0.43, and 0.43, respectively. From this, it is clear that, overall, under the condition of other conditions being constant, the maximum subsidence gradually decreases with the increase of tensile strength. But it has obvious segmentation characteristics. When the tensile strength increases from 0.4 MPa to 5.4 MPa, the surface subsidence decreases from 3077 mm to 2900 mm. However, when the tensile strength increases from 5.4 MPa to 15.4 MPa, the maximum subsidence basically does not change with the increase of tensile strength. Therefore, there is a critical value of the tensile strength. When the tensile strength changes within a range less than the critical value, the maximum subsidence decreases with its increase. When the tensile strength exceeds the critical value, the increase of tensile strength has almost no effect on the maximum subsidence. The change of the maximum subsidence is relatively small, so the change of subsidence coefficient is not significant. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and tensile strength evolution is shown in Figure15. The function expressions of the maximum subsidence, subsidence coefficient, and tensile strength obtained by fitting are as follows:
(6)W=237.025∗exp−T4.316+2873.316R2=0.932,q=−6.429T+0.456R2=0.890.Figure 15
Fitting curve of maximum subsidence/subsidence coefficient and tensile strength.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.3.3. Influence of Cohesion of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different cohesion were drawn, as shown in Figure16.Figure 16
Surface subsidence curve of different cohesion in strike direction.It can be seen from the simulation results that when the cohesion gradually increases from 1.2 MPa to 30.2 MPa, the maximum subsidence are 3702 mm, 3373 mm, 2965 mm, 2609 mm, 2287 mm, 2251 mm, and 2251 mm, respectively, and the corresponding subsidence coefficients are 0.54, 0.50, 0.44, 0.38, 0.34, 0.33, 0.33, and 0.33, respectively. As a whole, the maximum subsidence and subsidence coefficient gradually decrease with the increase of cohesion. However, it has obvious segmentation characteristics. When the cohesion increases from 1.2 MPa to 15.2 MPa, the surface subsidence value and coefficient decrease quickly, from 3702 mm to 2287 mm; however, when the cohesion increases from 15.2 MPa to 30.2 MPa, the subsidence value and subsidence coefficient decrease slowly. Therefore, there is a critical value of cohesion, and when the cohesion changes within a range less than the critical value, the maximum subsidence and subsidence coefficient change greatly. When the cohesion exceeds the critical value, the maximum subsidence and subsidence coefficient will no longer change obviously with the increase of cohesion. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and cohesion evolution is shown in Figure17. The function expressions of maximum subsidence, subsidence coefficient, and cohesion obtained by fitting are as follows:
(7)W=1854.91∗exp−c6.21+2201.23R2=0.993,q=0.270∗exp−c6.543+0.322R2=0.993.Figure 17
Fitting curve of maximum subsidence/subsidence coefficient and cohesion.
(a)
Maximum subsidence(b)
Subsidence coefficient


### 4.3.4. Influence of Internal Friction Angle of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different internal friction angles were drawn, as shown in Figure18.Figure 18
Surface subsidence curve of different internal friction angles in strike direction.It can be seen from the simulation results that when the internal friction angle gradually increases from 10° to 45°, the maximum subsidence values are 3938 mm, 3690 mm, 3486 mm, 3314 mm, 3118 mm, 2965 mm, 2826 mm, and 2688 mm, respectively, and the corresponding subsidence coefficients are 0.58, 0.54, and 0.51, respectively. It can be seen that, with the increase of the internal friction angle in the collapse zone, the maximum subsidence and subsidence coefficient basically show a linear decreasing trend under the same other conditions. The fitting curve of the relationship between the maximum surface subsidence value and subsidence coefficient and the evolution of internal friction angle is shown in Figure19. The function expressions of maximum subsidence, subsidence coefficient, and internal friction angle obtained by fitting are:
(8)W=−35.31φ+4224.07R2=0.991,q=−0.005φ+0.618R2=0.986.Figure 19
Fitting curve of maximum subsidence/subsidence coefficient and internal friction angle.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.3.1. Influence of Elastic Modulus of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different elastic modulus were drawn, as shown in Figure12.Figure 12
Surface subsidence curve of different elastic modulus in strike direction.It can be seen from the simulation results that when the elastic modulus gradually increases from 4.8 GPa to 39.8 GPa, the maximum subsidence are 3116 mm, 3023 mm, 2965 mm, 2899 mm, 2852 mm, 2798 mm, and 2778 mm, respectively, and the corresponding subsidence coefficients are 0.46, 0.44, 0.44, 0.43, 0.42, 0.41, 0.41, and 0.41, respectively. From this, it is clear that, with the increase of the elastic modulus of the collapse zone, the maximum subsidence of the ground surface gradually decreases under the same other conditions. However, compared with the influence of other parameters, the total change of the maximum subsidence caused by it is not significant. When the elastic modulus increases from 4.8 GPa to 39.8 GPa, the maximum subsidence decreases by only 338 mm. Therefore, the change of elastic modulus has no obvious influence on the subsidence coefficient. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient and elastic modulus evolution is shown in Figure13. The function expressions of maximum subsidence, subsidence coefficient, and elastic modulus obtained by fitting are:
(5)W=3.423E−0.056R2=0.986,q=−0.0014E+0.459R2=0.912.Figure 13
Fitting curve of maximum subsidence/subsidence coefficient and elastic modulus.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.3.2. Influence of Tensile Strength of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different tensile strength were drawn, as shown in Figure14.Figure 14
Surface subsidence curve of different tensile strength in strike direction.It can be seen from the simulation results that when the tensile strength gradually increases from 0.4 MPa to 15.4 MPa, the maximum subsidence values are 3077 mm, 3050 mm, 3022 mm, 3000 mm, 2968 mm, 2900 mm, 2891 mm, and 2891 mm, respectively, and the corresponding subsidence coefficients are 0.45, 0.45, 0.44, 0.44, 0.44, 0.43, 0.43, and 0.43, respectively. From this, it is clear that, overall, under the condition of other conditions being constant, the maximum subsidence gradually decreases with the increase of tensile strength. But it has obvious segmentation characteristics. When the tensile strength increases from 0.4 MPa to 5.4 MPa, the surface subsidence decreases from 3077 mm to 2900 mm. However, when the tensile strength increases from 5.4 MPa to 15.4 MPa, the maximum subsidence basically does not change with the increase of tensile strength. Therefore, there is a critical value of the tensile strength. When the tensile strength changes within a range less than the critical value, the maximum subsidence decreases with its increase. When the tensile strength exceeds the critical value, the increase of tensile strength has almost no effect on the maximum subsidence. The change of the maximum subsidence is relatively small, so the change of subsidence coefficient is not significant. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and tensile strength evolution is shown in Figure15. The function expressions of the maximum subsidence, subsidence coefficient, and tensile strength obtained by fitting are as follows:
(6)W=237.025∗exp−T4.316+2873.316R2=0.932,q=−6.429T+0.456R2=0.890.Figure 15
Fitting curve of maximum subsidence/subsidence coefficient and tensile strength.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.3.3. Influence of Cohesion of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different cohesion were drawn, as shown in Figure16.Figure 16
Surface subsidence curve of different cohesion in strike direction.It can be seen from the simulation results that when the cohesion gradually increases from 1.2 MPa to 30.2 MPa, the maximum subsidence are 3702 mm, 3373 mm, 2965 mm, 2609 mm, 2287 mm, 2251 mm, and 2251 mm, respectively, and the corresponding subsidence coefficients are 0.54, 0.50, 0.44, 0.38, 0.34, 0.33, 0.33, and 0.33, respectively. As a whole, the maximum subsidence and subsidence coefficient gradually decrease with the increase of cohesion. However, it has obvious segmentation characteristics. When the cohesion increases from 1.2 MPa to 15.2 MPa, the surface subsidence value and coefficient decrease quickly, from 3702 mm to 2287 mm; however, when the cohesion increases from 15.2 MPa to 30.2 MPa, the subsidence value and subsidence coefficient decrease slowly. Therefore, there is a critical value of cohesion, and when the cohesion changes within a range less than the critical value, the maximum subsidence and subsidence coefficient change greatly. When the cohesion exceeds the critical value, the maximum subsidence and subsidence coefficient will no longer change obviously with the increase of cohesion. The fitting curve of the relationship between the maximum subsidence, subsidence coefficient, and cohesion evolution is shown in Figure17. The function expressions of maximum subsidence, subsidence coefficient, and cohesion obtained by fitting are as follows:
(7)W=1854.91∗exp−c6.21+2201.23R2=0.993,q=0.270∗exp−c6.543+0.322R2=0.993.Figure 17
Fitting curve of maximum subsidence/subsidence coefficient and cohesion.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 4.3.4. Influence of Internal Friction Angle of Collapse Zone on Surface Subsidence

By compiling and analyzing the simulation results, the comparison curves of surface subsidence in the direction of the strike under different internal friction angles were drawn, as shown in Figure18.Figure 18
Surface subsidence curve of different internal friction angles in strike direction.It can be seen from the simulation results that when the internal friction angle gradually increases from 10° to 45°, the maximum subsidence values are 3938 mm, 3690 mm, 3486 mm, 3314 mm, 3118 mm, 2965 mm, 2826 mm, and 2688 mm, respectively, and the corresponding subsidence coefficients are 0.58, 0.54, and 0.51, respectively. It can be seen that, with the increase of the internal friction angle in the collapse zone, the maximum subsidence and subsidence coefficient basically show a linear decreasing trend under the same other conditions. The fitting curve of the relationship between the maximum surface subsidence value and subsidence coefficient and the evolution of internal friction angle is shown in Figure19. The function expressions of maximum subsidence, subsidence coefficient, and internal friction angle obtained by fitting are:
(8)W=−35.31φ+4224.07R2=0.991,q=−0.005φ+0.618R2=0.986.Figure 19
Fitting curve of maximum subsidence/subsidence coefficient and internal friction angle.
(a)
Maximum subsidence(b)
Subsidence coefficient


## 5. Grey Relational Analysis and Establishment of Prediction Model of Subsidence Deformation

In order to analyze the relation between the above factors and surface subsidence, combined with the numerical simulation results under various schemes, the grey relational degree between the maximum subsidence and these seven factors was calculated by grey relational analysis, and the main factors affecting the maximum surface subsidence value were found out. The prediction model of surface subsidence was built based on multiple regression analysis.


### 5.1. Basic Theory of Grey Relational Analysis

The basic principle of grey relational theory is to judge whether there is a close relationship between data series, which is mainly achieved by comparing the similarity of development trajectories between data series curves, and the degree of close relationship is quantified by grey relational degree [39, 40]. The larger the grey correlation degree is, the more closely the two data sequences is. On the contrary, the smaller the grey relational degree is, the smaller the closeness of the two sequences is naturally.Grey relational analysis is mainly realized through the following steps: firstly, the reference sequence is determined to reflect the behavior characteristics of the research system. Secondly, the reference and comparison series are dimensionless. Thirdly, the gray relational coefficients between the reference and comparison series are obtained based on different resolution coefficients. Finally, the grey relational coefficient obtained based on the same resolution coefficient is sorted in grey correlation order. The general relational degree proposed by Professor Deng is one of the widely used grey relational degrees [41, 42], and its basic methods can be summarized as follows.Let the reference sequencex0 with n elements be expressed as:
(9)x0=x01,x02,⋯x0n.Compare sequencesx1,x2,⋯,xj as follows:
(10)x1=x11,x12,⋯x1n,x2=x21,x22,⋯x2n,⋮xj=xj1,xj2,⋯xjn.In essence, relative analysis is to calculate the difference of geometric shapes between reference curves and comparison curves, which is formed by connecting reference series and comparison series as data points. The relational degree can be expressed by the area between the comparison curve and the reference curve, which can be seen intuitively from the geometric shape. The correlation degree can be calculated as follows:(1)
Dimensionless raw dataBecause there are many factors in the system, and each has different physical meanings and measurement units, the dimensions of the data are different, and sometimes, the magnitude of the values are far different. In order to make the analysis more convenient, and to ensure the equivalence and order of the data, before comparing various factors, it is necessary to make the original data dimensionless and normalized. The processing formula is as follows:(11)Xi=xi−ximinximax−ximin,where i∈0,j.
(2)
Calculate relational coefficientThe similarity of geometric shapes between curves can be used to analyze the relational degree between factors, so the difference between curves can be used to measure the relational degree. The relational coefficient of curvei at point k can be expressed as:
(12)γx0k,xik=iminkminx0k−xik+ρ.imaxkmaxx0k−xikkx0k−xik+ρ.imaxkmaxx0k−xik,where ρ is the resolution coefficient and ρ takes values in the range of 0 to 1 and usually takes 0.5 [43–45]. iminkminx0k−xik is the minimum difference between two poles; imaxkmaxx0k−xik is the maximum difference between two poles.
(3)
Calculate relational degreeAccording to the grey relational space theory [46], the calculation formula of relational degree between two curves is:
(13)γX0,X0=1n∑k=1nγx0k,xik,where i∈1,j.


### 5.2. Grey Relational Analysis of the Maximum Surface Subsidence in the Study Area

The general relational degree is used to analyze the correlation between the maximum surface subsidence in the strike direction of the study area and the thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of collapse zone. Taking the maximum subsidence as the reference sequence, the initial matrix of the comparison sequence is constructed, and the thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters (elastic modulus, tensile strength, cohesion, and internal friction angle) of the collapse zone are taken as the comparison sequence. The initial matrix is processed dimensionless.The results of relational analysis are shown in Table2; it can be seen that the grey relational degrees of thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of collapse zone (elastic modulus, tensile strength, cohesion, and internal friction angle) are 0.865, 0.845, 0.922, 0.843, 0.861, 0.799, and 0.876, respectively. According to the results of grey relational analysis, the influence degree of seven factors on the maximum surface subsidence under numerical simulation is as follows: thickness of coal > internal friction angle > thickness loose layer > tensile strength > thickness of bedrock > elastic modulus > cohesion. Thickness of coal has the greatest influence on surface subsidence, and the cohesion of collapse zone is the least.Table 2
Results of grey relational degree analysis.
IndexDegree of association (resolution coefficient is 0.5)Relational degree rankingThickness of coal0.9221Internal friction angle0.8762Thickness of loose layer0.8653Tensile strength0.8614Thickness of bedrock0.8455Modulus of elasticity0.8436Cohesion0.7997


### 5.3. Multivariate Nonlinear Regression Model

Multiple nonlinear regression refers to nonlinear regression involving two or more independent variables in the process of causality studies. In building nonlinear regression models, it is necessary to select reasonable model functions and initial valuation of parameters based on specific data in order to improve the speed of model convergence and avoid model misfit [47, 48]. In this paper, the establishment of regression relations was implemented by SPSS software to establish a multivariate nonlinear regression model for the maximum subsidence containing the thickness of the loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of the collapse zone (modulus of elasticity, tensile strength, cohesion, and angle of internal friction). Set the fitting relationship between the initial maximum surface subsidence value and various factors as follows:
(14)W=x1H3+x2H2+x3H+x4D2+x5D+x6h+Exγ+x8⋅expx9⋅T+x10⋅expx11⋅c+x12φ+x13,where x1,x2,⋯,x13 are the coefficients of each term, h is the thickness of loose layer, m; D is the thickness of bedrock, m; H is the thickness of coal, m; E is the elastic modulus of collapse zone, GPa; T is the tensile strength of collapse zone, MPa; c is cohesion of collapse zone, MPa; φ is the internal friction angle of the collapse zone, °.The iterative calculation was performed according to the relevant theory of nonlinear iterative algorithm. Eventually, the prediction model converged to the residual sum of squares after several iterations, indicating that the prediction model structure was reasonable and effective, and thus, the fitted relationship equation between the maximum subsidence and each factor was obtained as:(15)W=4.455×10−6H3−0.006H2+1.459H+1.870D2−396.152D+721.397h+E−0.280−35exp0.143T−10438.323exp0.004c−30.307φ+28468.524.


### 5.4. Engineering Verification of Prediction Model

In order to verify the accuracy and applicability of the prediction model, it is used to predict the maximum subsidence of some coal mines in China. According to reference [49–51], it can be known that the geological conditions of Sandaogou coal mine, Baodian coal mine, and Zhaogu No. 1 coal mine are similar to those of Lilou coal mine, so they are selected to verify the prediction model. The required parameters were collated, substituted into the prediction model for calculation, and compared and analyzed with the actual measured data from the surface observatory, and the results are detailed in Table 3. As can be seen from Table 3, the relative error between the predicted result and the measured subsidence data is less than 10% through the prediction formula (15), which shows that the prediction model is basically consistent with the measured result, which indicates that it has certain reference significance for guiding the mining in similar mines.Table 3
Comparison of measured and predicted results of maximum subsidence in some coal mines.
NameThickness of loose layer/mThickness of bedrock/mThickness of coal/mElastic modulus/GPaTensile strength/MPaCohesion/MPaInternal friction anglePredicted value/mmMeasured value/mmSandaogou coal mine70676.62.84.32.93435023346Baodian coal mine200738.55.49.02.53639743931Zhaogu no. 1 mine400543.014.58.24.23727532992


## 5.1. Basic Theory of Grey Relational Analysis

The basic principle of grey relational theory is to judge whether there is a close relationship between data series, which is mainly achieved by comparing the similarity of development trajectories between data series curves, and the degree of close relationship is quantified by grey relational degree [39, 40]. The larger the grey correlation degree is, the more closely the two data sequences is. On the contrary, the smaller the grey relational degree is, the smaller the closeness of the two sequences is naturally.Grey relational analysis is mainly realized through the following steps: firstly, the reference sequence is determined to reflect the behavior characteristics of the research system. Secondly, the reference and comparison series are dimensionless. Thirdly, the gray relational coefficients between the reference and comparison series are obtained based on different resolution coefficients. Finally, the grey relational coefficient obtained based on the same resolution coefficient is sorted in grey correlation order. The general relational degree proposed by Professor Deng is one of the widely used grey relational degrees [41, 42], and its basic methods can be summarized as follows.Let the reference sequencex0 with n elements be expressed as:
(9)x0=x01,x02,⋯x0n.Compare sequencesx1,x2,⋯,xj as follows:
(10)x1=x11,x12,⋯x1n,x2=x21,x22,⋯x2n,⋮xj=xj1,xj2,⋯xjn.In essence, relative analysis is to calculate the difference of geometric shapes between reference curves and comparison curves, which is formed by connecting reference series and comparison series as data points. The relational degree can be expressed by the area between the comparison curve and the reference curve, which can be seen intuitively from the geometric shape. The correlation degree can be calculated as follows:(1)
Dimensionless raw dataBecause there are many factors in the system, and each has different physical meanings and measurement units, the dimensions of the data are different, and sometimes, the magnitude of the values are far different. In order to make the analysis more convenient, and to ensure the equivalence and order of the data, before comparing various factors, it is necessary to make the original data dimensionless and normalized. The processing formula is as follows:(11)Xi=xi−ximinximax−ximin,where i∈0,j.
(2)
Calculate relational coefficientThe similarity of geometric shapes between curves can be used to analyze the relational degree between factors, so the difference between curves can be used to measure the relational degree. The relational coefficient of curvei at point k can be expressed as:
(12)γx0k,xik=iminkminx0k−xik+ρ.imaxkmaxx0k−xikkx0k−xik+ρ.imaxkmaxx0k−xik,where ρ is the resolution coefficient and ρ takes values in the range of 0 to 1 and usually takes 0.5 [43–45]. iminkminx0k−xik is the minimum difference between two poles; imaxkmaxx0k−xik is the maximum difference between two poles.
(3)
Calculate relational degreeAccording to the grey relational space theory [46], the calculation formula of relational degree between two curves is:
(13)γX0,X0=1n∑k=1nγx0k,xik,where i∈1,j.


## 5.2. Grey Relational Analysis of the Maximum Surface Subsidence in the Study Area

The general relational degree is used to analyze the correlation between the maximum surface subsidence in the strike direction of the study area and the thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of collapse zone. Taking the maximum subsidence as the reference sequence, the initial matrix of the comparison sequence is constructed, and the thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters (elastic modulus, tensile strength, cohesion, and internal friction angle) of the collapse zone are taken as the comparison sequence. The initial matrix is processed dimensionless.The results of relational analysis are shown in Table2; it can be seen that the grey relational degrees of thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of collapse zone (elastic modulus, tensile strength, cohesion, and internal friction angle) are 0.865, 0.845, 0.922, 0.843, 0.861, 0.799, and 0.876, respectively. According to the results of grey relational analysis, the influence degree of seven factors on the maximum surface subsidence under numerical simulation is as follows: thickness of coal > internal friction angle > thickness loose layer > tensile strength > thickness of bedrock > elastic modulus > cohesion. Thickness of coal has the greatest influence on surface subsidence, and the cohesion of collapse zone is the least.Table 2
Results of grey relational degree analysis.
IndexDegree of association (resolution coefficient is 0.5)Relational degree rankingThickness of coal0.9221Internal friction angle0.8762Thickness of loose layer0.8653Tensile strength0.8614Thickness of bedrock0.8455Modulus of elasticity0.8436Cohesion0.7997


## 5.3. Multivariate Nonlinear Regression Model

Multiple nonlinear regression refers to nonlinear regression involving two or more independent variables in the process of causality studies. In building nonlinear regression models, it is necessary to select reasonable model functions and initial valuation of parameters based on specific data in order to improve the speed of model convergence and avoid model misfit [47, 48]. In this paper, the establishment of regression relations was implemented by SPSS software to establish a multivariate nonlinear regression model for the maximum subsidence containing the thickness of the loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of the collapse zone (modulus of elasticity, tensile strength, cohesion, and angle of internal friction). Set the fitting relationship between the initial maximum surface subsidence value and various factors as follows:
(14)W=x1H3+x2H2+x3H+x4D2+x5D+x6h+Exγ+x8⋅expx9⋅T+x10⋅expx11⋅c+x12φ+x13,where x1,x2,⋯,x13 are the coefficients of each term, h is the thickness of loose layer, m; D is the thickness of bedrock, m; H is the thickness of coal, m; E is the elastic modulus of collapse zone, GPa; T is the tensile strength of collapse zone, MPa; c is cohesion of collapse zone, MPa; φ is the internal friction angle of the collapse zone, °.The iterative calculation was performed according to the relevant theory of nonlinear iterative algorithm. Eventually, the prediction model converged to the residual sum of squares after several iterations, indicating that the prediction model structure was reasonable and effective, and thus, the fitted relationship equation between the maximum subsidence and each factor was obtained as:(15)W=4.455×10−6H3−0.006H2+1.459H+1.870D2−396.152D+721.397h+E−0.280−35exp0.143T−10438.323exp0.004c−30.307φ+28468.524.


## 5.4. Engineering Verification of Prediction Model

In order to verify the accuracy and applicability of the prediction model, it is used to predict the maximum subsidence of some coal mines in China. According to reference [49–51], it can be known that the geological conditions of Sandaogou coal mine, Baodian coal mine, and Zhaogu No. 1 coal mine are similar to those of Lilou coal mine, so they are selected to verify the prediction model. The required parameters were collated, substituted into the prediction model for calculation, and compared and analyzed with the actual measured data from the surface observatory, and the results are detailed in Table 3. As can be seen from Table 3, the relative error between the predicted result and the measured subsidence data is less than 10% through the prediction formula (15), which shows that the prediction model is basically consistent with the measured result, which indicates that it has certain reference significance for guiding the mining in similar mines.Table 3
Comparison of measured and predicted results of maximum subsidence in some coal mines.
NameThickness of loose layer/mThickness of bedrock/mThickness of coal/mElastic modulus/GPaTensile strength/MPaCohesion/MPaInternal friction anglePredicted value/mmMeasured value/mmSandaogou coal mine70676.62.84.32.93435023346Baodian coal mine200738.55.49.02.53639743931Zhaogu no. 1 mine400543.014.58.24.23727532992


## 6. Conclusion

In this paper, through FLAC3D numerical simulation, the evolution law of surface subsidence under different overburden size parameters and mechanical parameters of collapse zone is analyzed, and the prediction model of surface subsidence is established. The research results are as follows:(1)
Based on the measured data of surface movement and deformation in no. 1301 working face of Lilou coal mine, the law of surface substance is analyzed(2)
The response characteristics of different thickness of loose layer, thickness of bedrock, thickness of coal, and initial mechanical parameters of collapse zone (elastic modulus, tensile strength, cohesion, and internal friction angle) in the surface subsidence are analyzed by numerical simulation(3)
Through grey relational analysis, it is shown that under the condition of numerical simulation, the order of influence degree of seven factors on the maximum subsidence is thickness of coal > internal friction angle > thickness of loose layer > tensile strength > thickness of bedrock > elastic modulus > cohesion(4)
Through multiple nonlinear regression analysis, the prediction model of surface subsidence including various factors is established:(16)W=4.455×10−6H3−0.006H2+1.459H+1.870D2−396.152D+721.397h+E−0.280−35exp0.143T−10438.323exp0.004c−30.307φ+28468.524(5)
In order to verify the rationality of the prediction model, the measured subsidence data of some coal mine surface observation stations with similar mine conditions in China are compared with the prediction results. The comparison results show that the prediction model has good applicability, which has certain guiding and reference significance for similar mine conditions



---

*Source: 5495509-2022-08-22.xml*