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

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.


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.

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, S 1 is the deviation of the inflection point in the strike direction, and S 2 and S 3 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]: where k is inflection point offset coefficient, and H i 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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 18 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.

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]: 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]: 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: where tanβ is the tangent of the main effect angle, and W cm is the maximum ground subsidence at critical mining. This last can be written as: where m is mining thickness; q is the subsidence factor; and α is the dip angle of the coal seam.

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.

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., p 1 and p 2 ) 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, p 3 , is selected as follows:  Once the third point is selected, these three points (p1, p2 and p3) could be used to construct a Delaunay triangle. Using 1 3 p p   and 2 3 p p  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. Once the third point is selected, these three points (p 1 , p 2 and p 3 ) 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.

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: 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

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 (L 1 , L 2 , L 3 and L 4 ); P 1 , P 2 , P 3 and P 4 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: where (x 0 , y 0 ) is the coordinate of the extraction element incenter; a, b and c are the lengths of three sides of the triangular extraction element; (x 1 , y 1 ), (x 2 , y 2 ) and (x 3 , y 3 ) 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 d 1 for the boundary line L 1 , and d 2 for L 2 . 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 d 1 , d 2 , d 3 and d 4 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 P 1 ', P 2 ', P 3 'and P 4 ' 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 L 3 ; the BP of the boundary line is the projection of the incenter on the boundary line (i.e., P 3 '). As shown in Figure 5b (i.e., Case B), the incenter is out of two endpoints of the mining boundary line L 2 ; the BP of the boundary line is the endpoint of the boundary line which is the closet to the incenter (i.e., P 3 or P 2 '). These two cases can be determined by: where P i and P i+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: where d i is the BD of i-th mining boundary line, and S i´i s the deviation of the inflection point corresponding to the BP of i-th mining boundary line.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 18 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: 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: 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.   Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 18 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: 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: 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 (d i ) and BP for the i-th mining boundary line; (d) Calculate the deviation of the inflection point on the BP, that is, S i '; (e) Compare d i and S i '. 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 d i is greater than or equal to S i ' for each mining boundary line, the extraction area is in the calculation area.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 18 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.

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: 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 H min , 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·H min . Note that k is the inflection point offset coefficient.

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)  Referring to Equation (2), the coal mining-induced surface subsidence for a given extraction element within the mining area can be written as: where E is the domain of extraction element, or integral domain of the probability integration function, and I E is a definite integral: Clearly, the integrand f (x,y) satisfy the following condition over the domain E: Letting D E represent a three-dimensional domain {x ∈ E, y ∈ E, 0 ≤ z ≤ f (x, y)} and V E represent the volume of the domain, it can be determined from the geometric interpretation of double the integral that the value of the integral I G is equal to the volume of the domain D E . Letting D A (auxiliary domain) represent a three-dimensional domain {x ∈ E, y ∈ E, 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: where V A is the volume of the auxiliary domain, and S E is the area of the extraction element. The volume ratio of the domain D E to the domain D A can be written as: where V E is the volume of the domain D E . If the volume ratio can be estimated, the volume of the domain D E or the integral I E can be obtained by: Then, substituting (15) into (10), the extraction element-induced ground subsidence can be rewritten as: As shown in Figure 7, by randomly rolling a larger number of three-dimensional points over the domain D A , 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 D E to the number of points within the domain D A .
µ ≈ n E n A where n E is the number of points within the domain D E , and n A is the number of points within the domain D A (or the number of rolling points). Substituting Equation (17) into Equation (16), the extraction element-induced ground subsidence can be estimated by: 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: (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 Note that the Monte Carlo integration method-based ratio estimation of n E to n A is guaranteed to converge with the in-situ ratio of V E to V A 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., n A ). 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: (e) Estimate the extraction element-induced ground subsidence by Equation (18). 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: (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 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.

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.

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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 18 reach the same calculation accuracy as the Monte Carlo method, the summation method needs more intensive triangular division, which is time-consuming.

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.

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  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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 18 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.    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. 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. 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.

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

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.

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

Conflicts of Interest:
The authors declare no conflict of interest.