Next Article in Journal
Follow-Up Control and Image Recognition of Neck Level for Standard Metal Gauge
Next Article in Special Issue
A General Solution for the Errors in Variables (EIV) Model with Equality and Inequality Constraints
Previous Article in Journal
A New Challenge: Path Planning for Autonomous Truck of Open-Pit Mines in The Last Transport Section
Previous Article in Special Issue
Determination of the Long-Term Ground Surface Displacements Using a PSI Technique—Case Study on Wrocław (Poland)
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Extraction of Irregularly Shaped Coal Mining Area Induced Ground Subsidence Prediction Based on Probability Integral Method

1
Lunan Geo-Engineering Exploration Institute of Shandong Province, Jining 272100, China
2
Technology Innovation Center of Restoration and Reclamation in Mining-induced Subsidence Land, Ministry of Natural Resources, Jining 272100, China
3
School of Geodesy and Geomatics, Wuhan University, Wuhan 430000, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2020, 10(18), 6623; https://doi.org/10.3390/app10186623
Submission received: 13 July 2020 / Revised: 24 August 2020 / Accepted: 21 September 2020 / Published: 22 September 2020
(This article belongs to the Special Issue Land Subsidence: Monitoring, Prediction and Modeling)

Abstract

:
Underground coal mining-induced ground subsidence (or major ground vertical settlement) is a major concern to the mining industry, government and people affected. Based on the probability integral method, this paper presents a new ground subsidence prediction method for predicting irregularly shaped coal mining area extraction-induced ground subsidence. Firstly, the Delaunay triangulation method is used to divide the irregularly shaped mining area into a series of triangular extraction elements. Then, the extraction elements within the calculation area are selected. Finally, the Monte Carlo method is used to calculate extraction element-induced ground subsidence. The proposed method was tested by two experimental data sets: the simulation data set and direct leveling-based subsidence observations. The simulation results show that the prediction error of the proposed method is proportional to mesh size and inversely proportional to the amount of generated random points within the auxiliary domain. In addition, when the mesh size is smaller than 0.5 times the minimum deviation of the inflection point of the mining area, and the amount of random points within an auxiliary domain is greater than 800 times the area of the extraction element, the difference between the proposed method-based subsidence predictions and the traditional probability integral method-based subsidence predictions is marginal. The measurement results show that the root-mean-square error of the proposed method-based subsidence predictions is smaller than 3 cm, the average of absolute deviations of the proposed method-based subsidence predictions is 2.49 cm, and the maximum absolute deviation is 4.05 cm, which is equal to 0.75% of the maximum direct leveling-based subsidence observation.

1. Introduction

Ground subsidence induced by underground coal mining activities is a very complicated process, which depends on mining thickness, mining depth, production methods, overlying strata properties and many other factors [1]. Because coal mining-induced land subsidence prediction plays an important role in underground mining under buildings, water bodies, railways and other important structures [2], a number of subsidence prediction methods have been proposed since the middle of the last century. Based on the differences between modeling approaches, these existing methods could be divided into three categories, including empirical methods, influence function methods and theoretical simulation methods. By making use of a larger number of in-situ subsidence observations, empirical methods could be used to develop prediction models (e.g., typical curves or profile functions) for ground movement and deformation, and the empirical formula which is used to calculate prediction parameters via the mining and geological parameters of the mining area [3,4,5,6,7]. Empirical methods are simple to use and yield fairly satisfactory results; however, the application of the methods is restricted to areas where the geological and mining conditions are identical, and can only be used to predict subsidence for simple two-dimensional problems of rectangular extraction. Theoretical simulating methods (e.g., the finite element method, boundary element simulation method and discrete element method) consider the material of the overlying strata, such as the cohesion less stochastic or elastic theoretical models [8,9,10,11,12,13]. Based on the model, mining-induced strata or ground movement and deformation could be calculated. Although the methods could be treated as very useful tools for investigating the mechanism of underground coal mining-induced ground subsidence, their effectivity in the accuracy of their subsidence prediction has not been etablished. Influence function methods (e.g., Konthe theory, Sann’s method, Beyer’s method or Litwiniszyn’s method) fall somewhere between empirical methods and theoretical simulating methods [14,15,16,17]. The methods typically divide the extraction area into a series of infinitesimal areas (or elementary cubes) [18,19,20]. Based on theoretical investigation, the methods could develop influence functions to describe the influence of the extraction of infinitesimal areas. The subsidence induced by the extraction of a whole mining area is equivalent to the amount of subsidence induced by the extraction of an infinitesimal area. The parameters of influence functions could be obtained by making use of in-situ subsidence observations in general. Influence methods can be mathematically validated and are applicable in various mining conditions. While the methods would become more complicated in the case of, for example, extensive irregular area mining-induced subsidence prediction, they predict symmetrical subsidence profiles about the trough center, which is not always the case. In addition, artificial neural network-based methods have also recently been used for the prediction of mining-induced subsidence [21,22,23,24,25,26]. Artificial neural network-based subsidence prediction seems to be a promising method for predicting coal mining-induced subsidence, although there are still many problems, such as the unexplained behavior of the network, the difficulty of determining the proper network structure, and the unknown reliability.
The Probability Integration Method (PIM), the theory of which was developed by Knothe and functions via Knothe’s theory, is one of the influential function methods [27,28,29,30]. The method makes use of the probability integration function as its main function. As the method is more applicable in practice and could produce a satisfactory accuracy, especially for the prediction of the extraction of rectangular mining area-induced subsidence, the method has been widely used in the prediction of underground coal mining-induced subsidence throughout the world, especially in China, and is also one of the recommended methods in the regulations of the Nation Coal Bureau [31]. However, because the difficulties of the accurate determination of the calculation area, and the complicated calculation of double integrals with an irregularly shaped integral domain, it is still a challenging problem to predict ground subsidence induced by the extraction of an irregularly shaped coal mining area based on PIM. Generally, the irregularly shaped mining area is divided into multiple rectangular sub-mining areas in order to calculate mining subsidence individually. However, this method necessitates modifying the prediction parameters for each sub-mining area, which is time-consuming. The division is subjective, and cannot approximate the in-situ mining area with the calculation area accurately, resulting in a worse agreement between the subsidence predictions and the observation ones [32]. In this paper, a new PIM-based method is proposed to predict ground subsidence induced by the extraction of irregularly shaped coal mining area, by making use of the Delaunay triangular mesh method and the Monte Carlo integration method. This method is usually applicable, since the Delaunay triangular mesh method could separate the arbitrarily shaped mining area and the calculation area with sufficient accuracy if the mesh size is small enough. The method is easily implemented because of the simplicity of the calculation of the double integral using the Monte Carlo method.
The remainder of this paper is organized as follows. In the next section, the fundamentals of the PIM-based subsidence prediction are provided. Section 3 presents the details of the proposed subsidence prediction method. Section 4 evaluates comprehensively the proposed method, using two data sets, which are the simulation data and direct leveling-based subsidence observations. Section 5 concludes this paper.

2. Fundamental of PIM-Based Ground Subsidence Prediction

PIM is a coal mining subsidence prediction method which is based on the stochastic medium theory. Because of the advantages, such as simple parameters and clear physical definition, the method has been widely used for predicting ground subsidence, tilt, curvature, horizontal displacement and deformation. Figure 1 shows a typical Cartesian coordinate system which is used for PIM-based coal mining subsidence prediction. The region of OABC shown in the figure denotes a mining area, S1 is the deviation of the inflection point in the strike direction, and S2 and S3 are the deviations of the inflection point in the lower and upper parts of the dip direction, respectively. All of those three deviations of the inflection point are related to mining depth, and can be calculated by [3]:
S i = k H i ( i = 1 , 2 , 3 )
where k is inflection point offset coefficient, and Hi is the average mining depth in the corresponding direction. Once the deviations of the inflection point are obtained, the calculation area (O’A’B’C’), which is denoted by D in the figure, can be determined from the mining area. Before the prediction of mining-induced subsidence, a coordinate system should be established. In general, the principles of coordinate axes are as follows: (a) the x axis is parallel to the strike direction of the coal seam and the right side of the strike direction is the positive direction of the x axis; (b) the y axis is parallel to the dip direction of the coal seam and the upper part of the dip direction is the positive direction of the y axis. In addition, the lower left corner of the mining area is taken as the origin of the coordinate system.
Based on the PIM and the coordinate system established above, coal mining-induced ground subsidence can be obtained by [30]:
W ( ε , η ) = W c m D 1 r 2 e π ( x ε ) 2 + ( y η ) 2 r 2 d x d y
where W(ε,η) is the coal mining-induced subsidence at the ground point with the coordinate of (ε,η); (x,y) is the planner coordinate of the mining unit, and r is the main influence radius of the mining area. This last can be related to the mining depth H by:
r = H tan β
where tanβ is the tangent of the main effect angle, and Wcm is the maximum ground subsidence at critical mining. This last can be written as:
W c m = m q cos α
where m is mining thickness; q is the subsidence factor; and α is the dip angle of the coal seam.

3. Proposed Ground Subsidence Prediction Method

In this Section, the details of the proposed subsidence prediction method are presented. The Delaunay triangulation method (DTM) is used to divide the irregularly shaped mining area into a series of triangular extraction elements at first. Then, the extraction elements within the calculation area are selected by comparing the deviation of the inflection point and the distance from the incenter of the element to the mining boundaries. Finally, the Monte Carlo method is used to calculate the extraction element-induced subsidence.

3.1. Irregularly Shaped Mining Area Segmentation Using DTM

The Delaunay triangulation method, which was invented by Boris Delaunay in 1934, is an important preprocessing technique for computational geometry and mathematics [33]. The method could divide arbitrarily shaped polygons (including convex polygons and concave polygons) into multiple triangular sub-regions. Because DTM has many advantages (e.g., optimality, uniqueness) compared with other triangular mesh methods, and can be easily implemented by programing languages, it is more applicable in the segmentation of irregularly shaped mining areas. Based on DTM and the Monte Carlo integration method, computer software (Subsidence Observation Data processing Program, SODP) has been developed by the Lunan Geo-Engineering Exploration Institute of Shandong Province and Wuhan University. The software was used to perform the following numerical calculation. Figure 2a shows the flowchart of the DTM-based segmentation algorithm of an irregularly shaped mining area. The inputs of the algorithm are the vertex coordinates of the mining area and mesh size. Based on the inputting parameters, mesh points in the mining area could be generated. There are two kinds of generated mesh points; one is the points at the boundary of the mining area and the other is the points in the interior of the mining area. The former could be generated by inserting points along with the boundary lines of the mining area (parallel to the direction of the x axis and y axis, inserting mesh points respectively), and the latter could be generated. The distance between two adjacent mesh points would be equal to the input mesh size in general. After the generation, all the mesh points are marked as “unused”. Two adjacent mesh points (e.g., p1 and p2) are selected to construct a vector, p 1 p 2 (or starting edge), as shown in Figure 2b, and at the same time these two points are marked as “used”; then, the third point, p3, is selected as follows:
(a)
Search for a point p3i from unused points. If the point satisfies
p 1 p 2 × p 1 p 3 i > 0
then add the point into the point set P. Repeat the search for each unused point;
(b)
Select a point, p3i, from P constructing two vectors p 1 p 3 i and p 2 p 3 i , then calculate the angle between these two vectors (α3i). Add the angle into an angle set Λ. Repeat the selection and calculation for each point in the point set P;
(c)
Search the maximum angle from Λ. Assuming the maximum angle is calculated based on vectors of p 1 p 3 k and p 2 p 3 k , p3k is the third point (p3) of a Delaunay triangle.
Once the third point is selected, these three points (p1, p2 and p3) could be used to construct a Delaunay triangle. Using p 1 p 3 and p 2 p 3 as the starting edge respectively and repeating the process above produces the other Delaunay triangles. When all the points are marked as “used”, the segmentation algorithm is finished, and the irregular mining area has been divided into multiple triangular extraction elements. Note that, because of the uniqueness property of the Delaunay triangle, constructing the starting edge with any two adjacent points would produce the same result. Figure 3 shows an example of the segmentation result of an irregularly shaped mining area by making use of DTM. The minimum and maximum lengths of the mining area boundary line are 26.8 m and 180.0 m, respectively, and the mesh size is set as 10.0 m. It can be seen that all the generated mesh points are evenly distributed in the mining area. Due to the optimality property of the Delaunay triangle, the lengths of the three sides are almost the same for each extraction element. This characteristic is important, because it is helpful in determining the calculation area accurately.

3.2. Selection of Extraction Elements within the Calculation Area

By making use of DTM, an irregularly shaped mining area could be divided into a series of triangular extraction elements. Because of the cantilever effect of the coal wall, only the extraction elements within the calculation area would induce ground subsidence. Although the calculation area can be easily determined though deviations of the inflection point for the regular-shaped mining area (i.e., rectangle-shaped mining area), the determination of the calculation area for an irregularly shaped mining area is a challenging problem, due to the complicated mining area boundary. Due to the fact that the calculation area is related to the deviation of the inflection point, it is possible to determine whether an extraction element is within the calculation area by comparing the deviation of the inflection point and the distance from the mining area boundary to the incenter of the extraction element. If the latter is greater than the former, the extraction element is free of the effect of the cantilever, that is, the element is within the calculation, whilst this is not true otherwise.
As shown in Figure 4, a quadrilateral mining area is constructed by four boundary lines (L1, L2, L3 and L4); P1, P2, P3 and P4 are four vertexes of the area. C is the incenter of an extraction element. The coordinates of the incenter of the element can be calculated by:
{ x 0 = a x 1 + b x 2 + c x 3 a + b + c y 0 = a y 1 + b y 2 + c y 3 a + b + c
where (x0, y0) is the coordinate of the extraction element incenter; a, b and c are the lengths of three sides of the triangular extraction element; (x1, y1), (x2, y2) and (x3, y3) are the coordinates of three vertexes of the extraction element. Boundary distance (BD), which is the distance from the extraction element to a mining boundary, is treated as the minimum distance from the incenter of the element to the boundary line. For exampe, the BD from the incenter of the element to a mining boundary line is d1 for the boundary line L1, and d2 for L2. Therefore, the amount of BDs is equals to the amount of mining area boundary lines for each extraction element. As shown in Figure 4, four BDs are termed as d1, d2, d3 and d4 respectively for the extraction element. For each BD, there is a corresponding point on the boundary line, which construct a BD with the incenter of the extraction element. These points (or boundary points, BPs) are termed as P1, P2, P3’and P4 respectively in Figure 4. Note that, because of the different relative positions between the incenter of the extraction element and the boundary line of the mining area, the determination of BPs is different. As shown in Figure 5a (i.e., Case A), the incenter of the extraction element is located between two endpoints of the mining boundary line L3; the BP of the boundary line is the projection of the incenter on the boundary line (i.e., P3). As shown in Figure 5b (i.e., Case B), the incenter is out of two endpoints of the mining boundary line L2; the BP of the boundary line is the endpoint of the boundary line which is the closet to the incenter (i.e., P3 or P2). These two cases can be determined by:
{ P i P i + 1 × P i C > 0 P i + 1 P i × P i + 1 C > 0 ( Case   A )
{ P i P i + 1 × P i C > 0 P i + 1 P i × P i + 1 C < 0 or { P i P i + 1 × P i C < 0 P i + 1 P i × P i + 1 C > 0 ( Case   B )
where Pi and Pi+1 are two endpoints of a mining boundary line. Once a BP is determined, the mining depth on the BP can be easily obtained by making use of the points for which the mining depth is known in advance and determined via the interpolation method (e.g., Kriging interpolation method). The deviation of the inflection point corresponding to the BP can be calculated by substituting the mining depth on the BP into Equation (1). If Equation (9) is true for each BP of an extraction element, the extraction element could be treated as an element within the calculation area:
S i < d i
where di is the BD of i-th mining boundary line, and Si´ is the deviation of the inflection point corresponding to the BP of i-th mining boundary line.
A flowchart of the determination of an extraction element within the calculation area is shown in Figure 6. Note that n is the amount of mining boundary lines. The corresponding algorithm can be described as follows:
(a)
Input vertex coordinates of a triangular extraction element and the vertex coordinates of the mining area;
(b)
Calculate the incenter of the extraction element by Equation (6);
(c)
Determine the relative position of the incenter and the i-th mining boundary line by Equation (7) and Equation (8); then calculate the BD (di) and BP for the i-th mining boundary line;
(d)
Calculate the deviation of the inflection point on the BP, that is, Si;
(e)
Compare di and Si. If the former is smaller than the latter, the extraction element is not in the calculation area. Otherwise, make i equal to i + 1, and repeat (c)–(e) for the (i + 1)-th mining boundary line;
(f)
If di is greater than or equal to Si’ for each mining boundary line, the extraction area is in the calculation area.
Repeating the above process for each extraction element, all the extraction elements within the calculation can be obtained. Note that, in order to guarantee the sufficient approximation to the in-situ calculation area, the area of the extraction element should be small enough. Assuming the minimum mining depth of the mining area is Hmin, the data processing results shows that there is a good agreement between the proposed method-based calculation area and the in-situ one if the mesh size is smaller than 0.5·k·Hmin. Note that k is the inflection point offset coefficient.

3.3. Extraction Element-Induced Subsidence Prediction Based on the Monte Carlo Method

In Section 3.2, the extraction elements within the calculation area have been selected. The mining of extraction element-induced ground subsidence could be obtained by making use of Equation (2). However, given the complicated calculation of the double integral function and integral domain, it is hard to calculate the extraction element induced surface subsidence using Equation (2) directly. Referring to Equation (2), the coal mining-induced surface subsidence for a given extraction element within the mining area can be written as:
W E ( ε , η ) = W c m r 2 E e π ( x ε ) 2 + ( y η ) 2 r 2 d x d y = W c m r 2 I E
where E is the domain of extraction element, or integral domain of the probability integration function, and IE is a definite integral:
I E = E f ( x , y ) d x d y = E e π ( x ε ) 2 + ( y η ) 2 r 2 d x d y
Clearly, the integrand f(x,y) satisfy the following condition over the domain E:
0 f ( x , y ) 1
Letting DE represent a three-dimensional domain {xE, yE, 0 ≤ zf(x, y)} and VE represent the volume of the domain, it can be determined from the geometric interpretation of double the integral that the value of the integral IG is equal to the volume of the domain DE. Letting DA (auxiliary domain) represent a three-dimensional domain {xE, yE, 0 ≤ z ≤ 1}, it can be easily derived that the volume of the auxiliary domain is equal to the area of the extraction element, which can be obtained by the vertex coordinates of the element:
V A = S E = 1 2 [ ( x 1 y 2 x 2 y 1 ) + ( x 2 y 3 x 3 y 2 ) + ( x 3 y 1 x 1 y 3 ) ]
where VA is the volume of the auxiliary domain, and SE is the area of the extraction element. The volume ratio of the domain DE to the domain DA can be written as:
μ = V E V A
where VE is the volume of the domain DE. If the volume ratio can be estimated, the volume of the domain DE or the integral IE can be obtained by:
V E = I E = μ V A = μ S E
Then, substituting (15) into (10), the extraction element-induced ground subsidence can be rewritten as:
W E ( ε , η ) = W c m r 2 μ S E
As shown in Figure 7, by randomly rolling a larger number of three-dimensional points over the domain DA, it can be derived via the Monte Carlo integration method [34] that the volume ratio in Equation (14) is approximately equal to the ratio of the number of points within the domain DE to the number of points within the domain DA.
μ n E n A
where nE is the number of points within the domain DE, and nA is the number of points within the domain DA (or the number of rolling points). Substituting Equation (17) into Equation (16), the extraction element-induced ground subsidence can be estimated by:
W E ( ε , η ) W c m r 2 n E n A S E
Note that the Monte Carlo integration method-based ratio estimation of nE to nA is guaranteed to converge with the in-situ ratio of VE to VA when the number of rolling points increases to infinity [35]. In addition, the accuracy of the Monte Caro method-based subsidence prediction is proportional to the number of rolling points (i.e., nA). The data processing results show that the difference between the Monte Caro method-based subsidence prediction and the in-situ one is marginal for a given extraction element if the number of rolling points is greater than or equal to 800 times the area of the extraction element.
The flowchart of the Monte Carlo method-based subsidence prediction for a given extraction element is shown in Figure 8, and the corresponding algorithm can be described as follows:
(a)
Input the vertex coordinates of a triangular extraction element;
(b)
Calculate the area of the element by Equation (18);
(c)
Generate nA random points over the domain DA;
(d)
Calculate the number of points within the domain DE. The random point (xi,yi,zi) within the domain can be determined by
z i f ( x i , y i )
(e)
Estimate the extraction element-induced ground subsidence by Equation (18).
Repeating the process described above produces ground subsidence predictions for each extraction element within the calculation area; the sum of each of the extraction element-induced ground subsidences is the final subsidence prediction induced by the extraction of the mining area.
Note that, although the summation method is comparatively simple in terms of the subsidence calculating process and the algorithm implemented [19], the data processing results show that the calculation accuracy of the summation method is inferior to that of the Monte Carlo integration method if the number of divided triangular extraction elements is equal. In addition, in order to reach the same calculation accuracy as the Monte Carlo method, the summation method needs more intensive triangular division, which is time-consuming.

4. Experimental Results

In this section, two field data sets were used to evaluate the proposed method. The first data set was obtained by simulating the extraction of the rectangular mining area. The second data set was collected from a coal mine in Shandong Province, China.

4.1. Simulation Result

Except for the typical PIM parameters and mining parameters (e.g., mining depth, mining thickness), the mesh size and the number of generated random points within the extraction element also need to be set before the prediction for the proposed method. In order to evaluate the impact of these two parameters on the accuracy of ground subsidence prediction, a rectangular mining area is simulated. As shown in Figure 9, the lengths of the mining area in the x direction and y direction are 600 m and 300 m, respectively. The dip angle of the coal seam, mining depth and mining thickness are set as 0°, −700 m and 300 cm, respectively. Typical PIM parameters of the YanZhou coalfield, China, are used to calculate ground subsidence. The PIM parameters and simulated mining area parameters are listed in Table 1. The ground subsidence predictions obtained by the traditional PIM described in Section 2 are treated as the in-situ ones. The contours of the in-situ ground subsidence are shown in Figure 9. The proposed method-based subsidence predictions on 402 ground points, which are within the subsidence contour of 1 cm, are selected to compare with the in-situ subsidence; the distribution of these ground points is also shown in Figure 9. The minimum and maximum in-situ subsidence of these selected ground points are 1.0 cm and 233.4 cm, respectively.
The approximation degree between the proposed method-based calculation area and the in-situ one is related to the ratio of mesh size to the minimum deviation of the inflection point of the mining area; the smaller ratio would produce a higher approximation between the two calculation areas, and thus more accurate subsidence predictions. Figure 10 shows a scatterplot of the RMS of the proposed method-based subsidence prediction errors versus the ratios of mesh size to deviation of the inflection point of the simulated mining area. The number of generated random points within an auxiliary domain is set as 1000 times the area of the corresponding extraction element. Note that, as the number of generated random points within the auxiliary domain is large enough, the Monte Carlo integration-induced subsidence prediction error could be ignored. It can be seen from Figure 10 that, before the ratio reaches 0.5, the RMS (root-mean-square) values of the proposed method-based subsidence prediction errors are very small, and slowly increase with the increase of the ratio. When the ratio is in the range of 0.5 to 3, the RMSE (root-mean-square error) values linearly increase with the increase of the ratio. With the increase of the ratio, the area of the extraction element and the BDs of the element also increase; once each BD of the element is greater than the corresponding deviation of the inflection point (which case corresponds to the ratio greater than 3 in Figure 10), all extraction elements would be treated as the elements within the calculation area, resulting in an overestimation of subsidence predictions for the proposed method. This could also be seen clearly in Table 2, which shows the Mean, STD (standard deviation) and RMS of the proposed method-based subsidence prediction errors under different ratios. When the ratio is greater than 3, the RMS of the proposed method-based subsidence prediction error reaches its maximum (about 8 cm) for the simulated mining area, and the RMSE tends towards stability as the calculation area does not change any more with the increase of the ratio.
In addition, it can be observed in Figure 10 and Table 2 that the proposed method-based subsidence predictions are much closer to the in-situ ones (RMSE smaller than 0.2 cm) when the ratio is less than or equal to 0.5 (or the mesh size is less than or equal to 0.5 times the deviation of the inflection point). This observation is useful because we could select a reasonable mesh size to minimize the mining area segmentation-induced prediction error.
If the mesh size is small enough, the mining area segmentation-induced subsidence prediction error could be ignored; in this circumstance, the accuracy of the proposed method-based subsidence prediction depends on the integration calculation accuracy. For a given extraction element, the integration calculation accuracy is proportional to the density of random points within the auxiliary domain, which is the ratio of the amount of generated random points within an auxiliary domain to the volume of the auxiliary domain; a larger ratio would produce a larger number of random points within the unit’s volume, resulting in an accurate convergence of the probability integral function, and thus more accurate subsidence predictions. Because the volume of the auxiliary domain is equal to the area of the corresponding extraction element, the ratio could also be expressed as a ratio of the amount of generated random points within an auxiliary to the area of the corresponding extraction element.
Taking the ratio of the amount of generated points to the area of the corresponding extraction element as the independent variable, and the RMS of the proposed method-based subsidence prediction errors as dependent variable, a scatterplot of the RMSE versus the ratio is shown in Figure 11. Note that the mesh size is set as 0.1 times the deviation of the inflection point (i.e., 2.25 m); because the mesh size is small enough, the mining area segmentation-induced subsidence prediction error could be ignored. It can be seen from Figure 11 that the RMS values of the proposed method-based subsidence prediction errors exponentially decrease with the increase of the ratio. When the ratio is greater than 800, the amount of generated random points within an extraction element is large enough, and the proposed method-based values of probability integral function are much closer to the in-situ ones for the proposed method. This could also be seen from Table 3, which shows the Mean, STD and RMS of the proposed method-based subsidence prediction errors under different ratios. This observation could be used to select a reasonable number of random points within the auxiliary domain for the purpose of minimizing Monte Carlo integration-induced prediction errors; specifically, the number of random points within an auxiliary domain should be greater than or equal to 800 times the area of the corresponding extraction element.

4.2. Validation Using Direct Leveling Based Subsidence Observations

The proposed method was validated by direct leveling-based subsidence observations, which were collected from a coal mine in Shandong Province, China. An irregular mining area is located in the northern No.1 district of the coal mine; as shown in Figure 12, the irregular mining area is constructed by eight longwall working faces. The maximum and minimum underground elevations are about −375 m and −670 m respectively, and the ground of the mining area is quite flat, with an elevation of 42.1 m before the underground extraction. The average coal seam thickness is 6.65 m. The average dip angle of the coal seam is about 10° and the dip direction and strike direction of the coal seam are as shown in Figure 12. The maximum lengths of the mining area in these two directions are about 1560 m and 1010 m, respectively. The starting mining date of the first working face (No. 1301) of the mining area was 12 March 2006, and the last working face (No.1308) was terminated on 16 May 2015. The terminating date of each working face is also shown in Figure 12. In January 2006, three ground observation lines (observation line of A, B and C) were set up, with a total of 130 observation points; the relative position of the mining area and these three observation lines is as shown in Figure 12. The observation line of A, B and C includes 46, 59 and 25 observation points respectively. Subsidence observations of these ground points were measured 22 times with fourth-order direct leveling from 25 February 2006 to 9 March 2017. However, because many points were destroyed during the observation period by human or subsidence-induced water pooling, only 46 surface points were available when the last measurement was performed. The numbers of available ground observation points are 15, 19 and 12 for the observation line of A, B and C respectively. The last subsidence observations of these available ground points are given in Figure 13.
According to the previously observed data in the southern No.1 district of the coal mine, the subsidence prediction parameters could be obtained, including the subsidence factor (0.83), the tangent of the main effect angle (2.0), and the inflection point offset coefficient (0.05). A minimum deviation of the inflection point of the mining area of 0.5 times was selected as the mesh size for the proposed method, which is equal 10.4 m, and the amount of random points within an extraction element was set as 800 times the area of the element. By making use of these parameters, the proposed method-based subsidence predictions on the ground points and the predicted subsidence contours could be obtained. The predicted subsidence contours of the mining area are also shown in Figure 12. Taking the direct leveling-based subsidence observations as the in-situ ones, the proposed method-based subsidence prediction error could be calculated by deducting the direct leveling-based subsidence observation from the predicted one. A scatter plot of the prediction error versus the distance from the ground point to the starting point of the corresponding observation line is shown in Figure 14. Table 4 shows the mean error, the error STD and the RMSE of the proposed method-based subsidence predictions. It can be seen from Figure 14 and Table 4 that there exists a significantly positive mean error, indicating that the direct leveling-based subsidences are greater than the predicted ones. This is probably because the ground subsidence induced a decline of the water table, and thus produced an additional subsidence in the ground [36]. The average of the absolute deviations of the proposed method-based subsidence predictions is 2.49 cm, and the maximum of the absolute deviation is 4.05 cm, which is equal to 0.75% of the maximum direct leveling-based subsidence observation. Basically, the results show a good agreement between the proposed method-based subsidence predictions and the direct leveling-based subsidence observations.

5. Conclusions

In order to solve the problem of subsidence prediction for an irregularly shaped coal mining area using PIM, a new DTM- and Monte Carlo integration-based method has been proposed. In Section 2, the fundamentals of PIM-based subsidence prediction are introduced. Details of the proposed method are provided in Section 3. In that section, DTM is used to firstly divide the irregularly shaped mining area into a series of triangular extraction elements. Then, the extraction elements within the calculation area are selected by comparison of the distance from the incenter of the extraction element to the mining area boundary and the deviation of the inflection point. Finally, based on the PIM described in Section 2, the Monte Carlo method is used to calculate the triangular extraction element-induced ground subsidence. The method could convert the complicated calculation of the double integral with a triangular domain into a simple multiplication operation. By summing each ground subsidence induced by triangular extraction elements within the calculation area, the final ground subsidence prediction of the irregularly shaped mining area can be obtained. The proposed method is evaluated comprehensively in Section 4 by making use of two data sets, including the simulation data and direct leveling-based subsidence measurements. The simulation results show that the prediction error of the proposed method is proportional to mesh size and inversely proportional to the amount of generated random points within the auxiliary domain. The difference in subsidence predictions between the proposed method and the traditional PIM can be ignored when the mesh size is smaller than 0.5 times the minimum deviation of the inflection point of the mining area, and the amount of random points within an auxiliary domain are greater than 800 times the area of the extraction element. The measurement results show that the root-mean-square error of the proposed method-based subsidence predictions is smaller than 3 cm, the average of the absolute deviations of the proposed method-based subsidence predictions is 2.49 cm, and the maximum of the absolute deviation is 4.05 cm, which is equal to 0.75% of the maximum direct leveling-based subsidence observation.
Future research will focus on using the proposed method to predict other coal mining-induced ground movements, including ground tilt, curvature, horizontal displacement and deformation. Predicting dynamic ground movement using the proposed method is another future research focus.

6. Patents

A software with copyright No. 2018SR057375 resulting from the work reported in this manuscript has been patented by National Copyright Administration, China.

Author Contributions

Conceptualization, Y.L. and H.B.; methodology, B.S.; software, Y.L.; validation, B.S., M.W. and G.L.; formal analysis, X.T.; data curation, X.T.; writing—original draft preparation, X.T.; visualization, Y.L.; funding acquisition, Y.L.; All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Key Laboratory of Geospace Environment and Geodesy, Ministry of Education of China, grant number 17-02-07, Key Research and Development Plan of Shandong Province, China, grant number 2019GSF109031.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Unlu, T.; Akcin, H.; Yilmaz, O. An integrated approach for the prediction of subsidence for coal mining basins. Eng. Geol. 2013, 166, 186–203. [Google Scholar] [CrossRef]
  2. Liang, L.; Kan, W.; Dawei, Z. AutoCAD-based prediction of 3D dynamic ground movement for underground coal mining. Int. J. Rock Mech. Min. Sci. 2014, 71, 194–203. [Google Scholar]
  3. Guoquan, Z.; Jixian, C. Mining under Buildings; China Coal Industry Publishing House: Beijing, China, 1983; p. 263. [Google Scholar]
  4. Jinzhuang, W.; Weimin, M.; Mengxun, N.; Xueli, F. Strata and Surface Movement in Coal Mine; China Coal Industry Publishing House: Beijing, China, 1981; p. 328. [Google Scholar]
  5. Asadi, A.; Shakhriar, K.; Goshtasbi, K. Profiling Function for Surface Subsidence Prediction in Mining Inclined Coal Seams. J. Min. Sci. 2004, 40, 142–146. [Google Scholar] [CrossRef]
  6. Bahuguna, P.P.; Srivastava, A.M.C.; Saxena, N.C. A critical review of mine subsidence prediction methods. Min. Sci. Technol. 1991, 13, 369–382. [Google Scholar] [CrossRef]
  7. Luo, Y. An improved influence function method for predicting subsidence caused by longwall mining operations in inclined coal seams. Int. J. Coal Sci. Technol. 2015, 3, 1–7. [Google Scholar] [CrossRef] [Green Version]
  8. Guo, Y.; Meng, F.; Chen, J. Numerical simulation analysis of surface subsidence under thick loose mining conditions. Coal Eng. 2014, 46, 103–105. [Google Scholar]
  9. Wei, S.; Qi, Z.; Yuanzhong, L.; Xiaoping, Z. A study of surface subsidence and coal pillar safety for strip mining in a deep mine. Environ. Earth Sci. 2018, 77, 627–639. [Google Scholar]
  10. Alejano, L.R.; Ramirez-Oyanguren, P.; Taboada, J. FDM predictive methodology for subsidence due to flat and inclined coal seam mining. Int. J. Rock Mech. Min. Sci. 1999, 36, 475–491. [Google Scholar] [CrossRef]
  11. Yaqiang, G.; Guangli, G. A Data-Intensive FLAC3D Computation Model: Application of Geospatial Big Data to Predict Mining Induced Subsidence. Comput. Modeling Eng. Sci. 2019, 5, 395–408. [Google Scholar]
  12. Abouzar, V.; John, A.; William, G. Mine-Scale Numerical Modeling of Longwall Operations. In Proceedings of the 2010 Underground Coal Operators’ Conference, Wollongong, Australia, 11–12 February 2010; pp. 115–124. [Google Scholar]
  13. Walter, K.; Ross, S.; Naj, A. Numerical Modelling of Mining Subsidence. In Proceedings of the 2006 Coal Operators’ Conference, Wollongong, Australia, 6–7 July 2006; pp. 313–326. [Google Scholar]
  14. Knothe, S. Observations of surface movements under influence of mining and their theoretical interpretation. In Proceedings of the European Conference on Ground Movement, Berlin, Germany, 19 February 1957; pp. 210–218. [Google Scholar]
  15. Keinhorst, H. Considerations on the problem of mining damage. Glukauf 1934, 70, 149–155. [Google Scholar]
  16. Kratzsch, H. Mining Subsidence Engineering; Springer: Berlin, Germany, 1983; p. 543. [Google Scholar]
  17. Whittaker, B.N.; Reddish, D.J. Subsidence: Occurrence, Prediction and Control; Elsevier: Amsterdam, The Netherlands, 1989; p. 528. [Google Scholar]
  18. Hejmanowski, R. Zur Vorausberechnung Förderbedingter Bodensenkungen über Erdöl-und Erdgaslagerstätten; TU Clausthal: Berlin, Germany, 1993. [Google Scholar]
  19. Malinowska, A.; Hejmanowski, R.; Dai, H. Ground movements modeling applying adjusted influence function. Int. J. Min. Sci. Technol. 2020, 30, 243–249. [Google Scholar] [CrossRef]
  20. Hejmanowski, R.; Malinowska, A.; Kwinta, A.; Ulmaniec, P. Modeling of land subsidence caused by salt cavern convergence applying Knothe’s theory. Markscheidewesen 2018, 25, 423–431. [Google Scholar]
  21. Garrett, J. Where and why artificial neural networks are applicable in civil engineering. Comput. Aided Civ. Infrastruct. Eng. 1994, 8, 129–130. [Google Scholar] [CrossRef]
  22. Rafie, M.; Samimi Namin, F. Prediction of subsidence risk by FMEA using artificial neural network and fuzzy inference system. Int. J. Min. Sci. Technol. 2015, 25, 655–663. [Google Scholar] [CrossRef]
  23. Pishro, M.; Khosravi, S.; Tehrani, S.M.; Mousavi, S.R. Modeling and zoning of land subsidence in the southwest of Tehran using artificial neural networks. Int. J. Hum. Cap. Urban. Manag. 2016, 1, 159–168. [Google Scholar]
  24. Tien Bui, D.; Shahabi, H.; Shirzadi, A.; Chapi, K.; Pradhan, B.; Chen, W.; Khosravi, K.; Panahi, M.; Bin Ahmad, B.; Saro, L. Land subsidence susceptibility mapping in south Korea using machine learning algorithms. Sensors 2018, 18, 2464. [Google Scholar] [CrossRef] [Green Version]
  25. Nefeslioglu, H.; Gokceoglu, C.; Sonmez, H. An assessment on the use of logistic regression and artificial neural networks with different sampling strategies for the preparation of landslide susceptibility maps. Eng. Geol. 2008, 97, 171–191. [Google Scholar] [CrossRef]
  26. Zhao, K.; Chen, S. Study on artificial neural network method for ground subsidence prediction of metal mine. Procedia Earth Planet. Sci. 2011, 2, 177–182. [Google Scholar] [CrossRef] [Green Version]
  27. Sroka, A.; Knothe, S.; Tajdus, K.; Misa, R. Underground exploitation inside safety pillar shafts considering the effective use of a coal deposit. Mineral. Resour. Manag. 2015, 31, 93–110. [Google Scholar] [CrossRef]
  28. Kwinta, A.; Gradka, R. Mining exploitation influence range. Nat. Hazards 2018, 94, 979–997. [Google Scholar] [CrossRef] [Green Version]
  29. Cheng, J.; Zhao, G.; Li, S. Predicting underground strata movements model with considering key strata effects. Geotech. Geol. Eng. 2018, 36, 621–640. [Google Scholar] [CrossRef]
  30. Baochen, L.; Guohua, L. The Basic Rule of Surface Movement of Coal Mine; Chinese Industry Press: Beijing, China, 1965. [Google Scholar]
  31. National Coal Bureau of P.R.C. Regulation of Mining and Pillar Leaving under Building, Water-Body, Railway and Main Underground Engineer; Coal Industry Press: Beijing, China, 2000; pp. 81–90. (In Chinese) [Google Scholar]
  32. Peixian, L.; Zhixiang, T.; Kazhong, D. Calculation of maximum ground movement and deformation caused by mining. Trans. Nonferrous Metals Soc. China 2011, 21, 562–569. [Google Scholar]
  33. Lee, D.T.; Schachter, B.J. Two Algorithms for Constructing a Delaunay Triangulation. Int. J. Comput. Inf. Sci. 1980, 9, 219–242. [Google Scholar] [CrossRef]
  34. Xiaoqun, W.; Kaitai, F. The effective dimension and quasi-Monte Carlo integration. J. Complex. 2003, 19, 101–124. [Google Scholar]
  35. Kong, A.; McCullagh, P.; Meng, X.L.; Nicolae, D.; Tan, Z. A Theory of Statistical Models for Monte Carlo Integration. J. R. Stat. Soc. Ser. B (Stat. Methodol.) 2003, 65, 604–618. [Google Scholar] [CrossRef]
  36. Yang, Y.; Song, X.F.; Zheng, F.D. Simulation of fully coupled finite element analysis of nonlinear hydraulic properties in land subsidence due to groundwater pumping. Environ. Earth Sci. 2015, 73, 4191–4199. [Google Scholar] [CrossRef]
Figure 1. Coordinate system used by PIM-based coal mining subsidence prediction. DD and SD denote dip direction and strike direction respectively.
Figure 1. Coordinate system used by PIM-based coal mining subsidence prediction. DD and SD denote dip direction and strike direction respectively.
Applsci 10 06623 g001
Figure 2. (a) Flowchart of DTM-based segmentation algorithm of irregularly shaped mining area. (b) Example of Delaunay triangulations generated by DTM.
Figure 2. (a) Flowchart of DTM-based segmentation algorithm of irregularly shaped mining area. (b) Example of Delaunay triangulations generated by DTM.
Applsci 10 06623 g002
Figure 3. Example of segmentation result of an irregularly shaped mining area using DTM.
Figure 3. Example of segmentation result of an irregularly shaped mining area using DTM.
Applsci 10 06623 g003
Figure 4. Geometry of mining area boundary lines and the incenter of an extraction element.
Figure 4. Geometry of mining area boundary lines and the incenter of an extraction element.
Applsci 10 06623 g004
Figure 5. Relative position of the incenter of the extraction element and mining boundary lines: (a) the incenter of the extraction element located between two endpoints of a mining boundary line; (b) the incenter of the extraction element located out of two endpoints of a mining boundary line.
Figure 5. Relative position of the incenter of the extraction element and mining boundary lines: (a) the incenter of the extraction element located between two endpoints of a mining boundary line; (b) the incenter of the extraction element located out of two endpoints of a mining boundary line.
Applsci 10 06623 g005
Figure 6. Flowchart of determination of an extraction element within the calculation area.
Figure 6. Flowchart of determination of an extraction element within the calculation area.
Applsci 10 06623 g006
Figure 7. The idea of the integral IE estimation model based on the Monte Carlo method.
Figure 7. The idea of the integral IE estimation model based on the Monte Carlo method.
Applsci 10 06623 g007
Figure 8. Flowchart of the extraction element-induced ground subsidence prediction using the Monte Carlo method.
Figure 8. Flowchart of the extraction element-induced ground subsidence prediction using the Monte Carlo method.
Applsci 10 06623 g008
Figure 9. Simulated mining area and traditional PIM predicted subsidence contours.
Figure 9. Simulated mining area and traditional PIM predicted subsidence contours.
Applsci 10 06623 g009
Figure 10. Scatterplot of RMS of the proposed method-based subsidence prediction errors versus ratio of mesh size to deviation of inflection point of the mining area.
Figure 10. Scatterplot of RMS of the proposed method-based subsidence prediction errors versus ratio of mesh size to deviation of inflection point of the mining area.
Applsci 10 06623 g010
Figure 11. Scatterplot of the RMS of the proposed method-based subsidence estimation errors versus the ratios of the amount of random points to the area of the corresponding extraction element.
Figure 11. Scatterplot of the RMS of the proposed method-based subsidence estimation errors versus the ratios of the amount of random points to the area of the corresponding extraction element.
Applsci 10 06623 g011
Figure 12. Mining area and the location of ground observation lines.
Figure 12. Mining area and the location of ground observation lines.
Applsci 10 06623 g012
Figure 13. Fourth-order direct leveling-based subsidence observations for the observation line of A, B and C.
Figure 13. Fourth-order direct leveling-based subsidence observations for the observation line of A, B and C.
Applsci 10 06623 g013
Figure 14. Scatterplot of the proposed method-based subsidence prediction error versus the distance from the surface point to the starting point of the corresponding observation line.
Figure 14. Scatterplot of the proposed method-based subsidence prediction error versus the distance from the surface point to the starting point of the corresponding observation line.
Applsci 10 06623 g014
Table 1. PIM parameters and simulated mining area parameters.
Table 1. PIM parameters and simulated mining area parameters.
Length on x Direction (m)Length on y Direction (m)α (°)H (m)M (cm)qtanβk
6003000−7003000.852.20.05
Table 2. Mean, STD and RMS of the proposed method-based subsidence estimation errors under different ratios of mesh size to deviation of inflection point.
Table 2. Mean, STD and RMS of the proposed method-based subsidence estimation errors under different ratios of mesh size to deviation of inflection point.
RatioMean (cm)STD (cm)RMSE (cm)
0.10.000.020.02
0.2−0.010.040.04
0.30.020.070.07
0.4−0.030.090.09
0.50.030.150.15
1.01.101.121.57
2.04.033.485.32
3.05.845.047.72
4.05.865.117.78
5.05.905.337.95
Table 3. Mean, STD and RMS of the proposed method-based subsidence estimations errors under different ratios of the amount of random points to the area of the corresponding extraction element.
Table 3. Mean, STD and RMS of the proposed method-based subsidence estimations errors under different ratios of the amount of random points to the area of the corresponding extraction element.
RatioMean (cm)STD (cm)RMSE (cm)
1002.4523.6823.81
200−1.7610.1710.32
300−0.535.775.79
4000.322.132.15
5000.181.011.03
600−0.150.750.76
7000.060.260.27
800−0.010.070.07
9000.040.050.06
10000.000.020.02
Table 4. Mean, STD and RMS of the proposed method-based subsidence prediction errors.
Table 4. Mean, STD and RMS of the proposed method-based subsidence prediction errors.
Observation LineMean (cm)STD (cm)RMSE (cm)
A2.381.062.61
B2.351.392.73
C2.471.142.72

Share and Cite

MDPI and ACS Style

Tan, X.; Song, B.; Bo, H.; Li, Y.; Wang, M.; Lu, G. Extraction of Irregularly Shaped Coal Mining Area Induced Ground Subsidence Prediction Based on Probability Integral Method. Appl. Sci. 2020, 10, 6623. https://doi.org/10.3390/app10186623

AMA Style

Tan X, Song B, Bo H, Li Y, Wang M, Lu G. Extraction of Irregularly Shaped Coal Mining Area Induced Ground Subsidence Prediction Based on Probability Integral Method. Applied Sciences. 2020; 10(18):6623. https://doi.org/10.3390/app10186623

Chicago/Turabian Style

Tan, Xianfeng, Bingzhong Song, Huaizhi Bo, Yunwei Li, Meng Wang, and Guohong Lu. 2020. "Extraction of Irregularly Shaped Coal Mining Area Induced Ground Subsidence Prediction Based on Probability Integral Method" Applied Sciences 10, no. 18: 6623. https://doi.org/10.3390/app10186623

APA Style

Tan, X., Song, B., Bo, H., Li, Y., Wang, M., & Lu, G. (2020). Extraction of Irregularly Shaped Coal Mining Area Induced Ground Subsidence Prediction Based on Probability Integral Method. Applied Sciences, 10(18), 6623. https://doi.org/10.3390/app10186623

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop