Liu Baochen and Dai Huayang developed the probabilistic integration prediction method, deriving a simplified solution based on random media theory and formulating a comprehensive, applicable approach (Liu and Dai, 2016).

(1) The semi-infinite mining prediction formula is as follows:

The insufficient settlement function is given in Equation 1, where the error function (erf) describes the integral results of the probability distribution, reflecting the cumulative effects of random medium movement. W x = W 0 2 erf π r x + 1 (1)

By differentiating Equation 1, the formula for calculating surface inclination (slope) is obtained, as shown in Equation 2. i x = W 0 r e − π 2 r 2 (2)

Further differentiation of Equation 2 yields the formula for calculating subsidence curvature along the main section, as shown in Equation 3. K x = − 2 π W 0 r 3 x e − π x 2 r 2 (3)

By differentiating Equation 3, the formula for calculating horizontal deformation along the main section is obtained, as shown in Equation 4. U x = b W 0 e − π x 2 r 2 (4)

By taking the derivative of Equation 4, the formula for calculating the horizontal deformation along the main section can be obtained, as shown in Equation 5. ε x = − 2 π b W 0 r 2 x e − π x 2 r 2 (5) Where, w ₀ −Maximum surface subsidence, mm; r−Influence radius, m; x−Horizontal position of any point; b−Horizontal movement coefficient.

(2) Prediction of Surface Movement and Deformation on the Main Inclined Cross-Section

When the coal seam is fully mined along the strike direction (parallel to the coal seam extension) with a sufficiently large mining range and complete surface subsidence, but mining in the dip direction (perpendicular to the strike) remains limited, the subsidence basin becomes asymmetric. To account for this asymmetry, the probabilistic integration method must be modified by introducing the main influence radii (r ₁ and r ₂ ) in the lower and upper mountain directions. These parameters correspond to the influence range beneath the coal seam (lower mountain) and above the coal seam (upper mountain), respectively.

The modified prediction formulas for the dip-direction main cross-section are presented in Equations 6–10: W 0 y = W y , t 1 − W y − L , t 2 (6) i 0 y = i y , t 1 − i y − L , t 2 (7) K 0 y = K y , t 1 − K y − L , t 2 (8) U 0 y = U y , t 1 − U y − L , t 2 (9) ε 0 y = ε y , t 1 − ε y − L , t 2 (10)

In the formula, t ₁ and t ₂ represent parameters for the lower and upper mountain boundaries, respectively. The main influence radii of the equivalent horizontal coal seam in these two directions are no longer identical. The radii r ₁ and r ₂ , can be determined using Equation 11: r 1 = H 2 tan ⁡ β 1 , r 2 = H 2 tan ⁡ β 2 (11)

Here, tanβ ₁ and tanβ ₂ represent the tangents of the main influence angles for the lower and upper mountain regions, respectively. The inclined working face length L is calculated using Equation 12: L = D 1 − s 1 − s 2 sin θ 0 + α sin ⁡ θ 0 (12) where D ₁ is the working face length, θ ₀ is the mining influence propagation angle, and α is the coal seam inclination angle.

The above formula applies to limited mining and does not differentiate between sufficient and insufficient mining conditions. Typically, insufficient mining values are obtained by multiplying the sufficient mining values by a mining coefficient C, as expressed in Equation 13: W 0 y = C W 0 2 erf π r x + 1 − erf π r x − D 1 + s 1 + s 2 + 1 (13)

In the probabilistic integration method, subsidence basin parameters are typically solved using mathematical techniques, with curve fitting based on actual measurement data being the most common approach. The mathematical solution process involves identifying a smooth curve that minimizes the sum of squared errors of discrete data points. The predicted parameters include the subsidence coefficient q , the tangent of the main influence angle tanβ, the mining influence propagation coefficient K, the horizontal movement coefficient b, and the inflection point offset distance S, etc.

The primary challenges in parameter estimation are as follows:

(1) High correlation between parameters: When multiple parameters exhibit strong correlation, they become indistinguishable during curve fitting. For instance, the mining coefficient C and the subsidence coefficient q are difficult to separate. In the formula Cw ₀ /2, the expression transforms into Cmqcosα/2, meaning that during fitting, it effectively scales mmm by a factor of C and q, making it difficult to distinguish between C and q.

The probability density function (PDF) method postulates that under conditions of insufficient mining, the displacement probability of discrete medium fragments gradually converges to the probability density function of a normal distribution (Guo et al., 2004). Assuming that the movement of fractured rock mass units follows a normal distribution, the subsidence model can be directly constructed using the probability density function, effectively mitigating strong correlations between parameters. The specific subsidence formula based on the probability density function method is presented in Equation 14: w x = w max ⁡ exp − π x − L 1 2 r 2 (14) i x = − 2 π w max x − L 1 r 2 exp − π x − L 1 2 r 2 (15) k x = 2 π w max r 2 2 π x − L 1 2 r 2 − 1 exp − π x − L 1 2 r 2 (16) U x = 2 π b w max x − L 1 r exp − π x − L 1 2 r 2 (17) ε x = 2 π w max r 2 2 π x − L 1 2 r 2 − 1 exp − π x − L 1 2 r 2 (18)

Equations 15–18 provide the calculation formulas for surface inclination, curvature, horizontal displacement, and horizontal deformation under the probability density function (PDF) method. These equations follow a similar structure to Equations 2–5.

In the formula: w max represents the maximum surface subsidence under insufficient mining conditions, given by w max = m q c o s , where m is the coal seam thickness, q is the subsidence rate, and α is the coal seam inclination angle; x is the horizontal distance from a point within the main cross-section to the center of the mined-out area; r is the main influence radius, expressed as r = H/tanβ, where H is the mining depth and β is the main influence angle. b is the horizontal movement coefficient. L is the working face length, where L ₁ = 1/2.

Compared to the probabilistic integration method, the PDF method directly uses the subsidence rate to predict maximum subsidence under insufficient mining conditions. This approach eliminates the need for two additional parameters—the mining coefficient and the subsidence coefficient—simplifying the prediction process.

The Xiao Ba Dang coal mine, operated by Xiao Ba Dang Mining Co., Ltd., is located in southwest Shenmu County, Yulin City, Shaanxi Province, and northeast Yuyang District. The first mining working face, 112,201, employs the large mining height longwall top coal caving method for roof management. After mining, the roof collapses and fractures, forming an “internal three zones” within the rock strata.

The 112,201 working face has a strike length of 4,660 m and a dip length of 350 m. The coal seam thickness ranges from 5.5 to 6.0 m, with an average of 5.8 m, and an inclination angle averaging 0.5°.

The surface topography of the 112,201 working face is relatively simple, with no buildings, railways, pipelines, or other obstructions. Consequently, the observation station in this study adopts a linear layout consisting of a strike observation line and a dip observation line, both perpendicular at the center of the subsidence basin.

Design of the strike observation line:

Given the thick surface loose layer of the 112,201 working face, the strike observation line is designed to cover half of the subsidence basin. Based on surface movement observations from the 20,102 working face of the Yushuwan coal mine, the strike line length is set at 900 m, with 300 m extending beyond the mining boundary of the cut eye and 600 m within it.

At one end of the strike line, three control points are placed outside the cut eye, spaced 100 m apart. Considering the average mining depth of 302 m, the point spacing is 25 m, resulting in a total of 37 monitoring points along the strike line. The monitoring point numbering follows a sequential pattern: Control points are labeled KZ01, KZ02, KZ03. Monitoring points are identified using a character and two Arabic numerals (e.g., Z01, Z02, …). “K” denotes a control point. “Z” denotes the strike observation line. The following numbers represent the sequential measuring points.

Design of the dip observation line:

To ensure complete subsidence basin coverage, the dip observation line length is set at 1,000 m, meeting observation requirements.

Each end of the dip line has two control points, totaling four, with a spacing of 100 m between control points and working points. Since the working face width (l = 350 m) is less than (1.4 H₀ = 423 m), the movement basin may not have reached full mining conditions along the dip. To accurately determine point spacing, 25 m intervals are used, with 40 working measuring points laid out. The order of the dip observation line control points is KQ01, KQ02, KQ04; the working measuring point number is composed of a character and two Arabic numerals. For example, the first monitoring point on the dip observation line is numbered as Q01, and the rest follow suit. Where “K” is the control point identification code, “Q” is the monitoring point identification code, and the following numbers are the sequential numbers of the measuring points.

The layout of strike and dip observation line monitoring points is illustrated in Figure 1.

The working face is fully mined along the strike direction, with a maximum settlement of 3,664 mm recorded at monitoring point Z18. In the dip direction, the maximum settlement reaches 3,327 mm at point Q20, and the settlement curve does not exhibit a horizontal section. A comprehensive analysis indicates that full mining conditions have not been achieved in the dip direction.

Using the nonlinear fitting function in Origin, a probabilistic integration fitting function was constructed based on Equation 13. Since the existing surface subsidence data use the left mining boundary as the origin, the independent variable x in the original formula requires a coordinate translation by s ₁ . Thus, the original coordinates are replaced with (x−s ₁ ). Incorporating the parameters of the 112,201 working face and surface observation data, the reconstructed probabilistic integration fitting function, which includes the inflection point shift, is presented in Equation 19: W 0 y = C 5800 × q 2 erf π × tan ⁡ β 302 x − s 1 + 1 − erf π × tan ⁡ β 302 x − 305 + s 2 + 1 (19)

Fitting the model using the aforementioned function reveals a correlation between C and q. Therefore, it is necessary to assign a value to C based on actual conditions. The value of C is determined as the ratio of the maximum subsidence in the dip direction to the maximum subsidence in the strike direction. After calculation, C is found to be 0.9.

Since the coal seam is nearly horizontal, it can be assumed that s ₁ and s ₂ are equal. Thus, a single variable s is used in the formula instead. The simplified probabilistic integration fitting function is presented in Equation 20: W 0 y = 2600 × q erf π × tan ⁡ β 302 x − s + 1 − erf π × tan ⁡ β 302 x − 305 + s + 1 (20)

Among them, q, tanβ and s are the predicted parameters.

The fitting curve is shown in Figure 2, and the fitting residual results are shown in Table 1:

The fitting parameters of the aforementioned curve are as follows: q = 1.02, tanβ = 2.97, and S = 82.38.

The probability density function (PDF) fitting function, constructed based on Equation 14, is presented in Equation 21: w x = − 5800 × q × exp − π x − 175 2 × t a n 2 β 302 2 (21)

The fitting curve is illustrated in Figure 3, while the fitting residuals are presented in Table 2.

The fitting parameters of the aforementioned curve are as follows: q = 0.61184, tanβ = 1.5167, and S = 82.38.

The sensitivity of q and tan β is now discussed and analyzed.

(1) q

Under the condition that tanβ remains constant at 1.5167, Figure 4 shows the fitted probability density function curves for different cases where q = 0.4, q = 0.5, and q = 0.61184.

Overall, under the constant condition of tanβ = 1.5167, an increase in the value of q leads to a reduction in the degree of surface subsidence, while the symmetrical U-shaped distribution characteristic of the subsidence remains unchanged. This trend indicates that q is a critical parameter affecting the depth of subsidence, with larger q values corresponding to shallower subsidence.

(2) tanβ

Under the condition that q remains constant at 0.61184, Figure 5 shows the fitted probability density function curves for different cases where tanβ = 1.5167, tanβ = 2, and tanβ = 3.

Under the condition that q remains constant at 0.61184, as tanβ increases from 1.5167 to two and then to 3, the following trends are observed: The maximum depth (peak value) of surface subsidence gradually increases. The range of subsidence distribution gradually narrows, with the curve becoming steeper, indicating that subsidence is more concentrated in the central region. The overall distribution retains a symmetrical U-shaped characteristic, with the basic shape unchanged by variations in tanβ.

By adjusting the expected parameters in the fitting of the predicted curve to the observation data, it is difficult to ensure that both the interior and the edges of the basin can be well fitted.

When q remains constant at 0.61184, and tanβ increases from 1.5167 to 2, and then to 3, the following trends are observed:

1. The maximum depth (peak value) of surface subsidence gradually increases.

When adjusting expected parameters to fit the predicted curve to observation data, it remains challenging to achieve an optimal fit for both the interior and edges of the subsidence basin.

Given that both the probabilistic integration method and the probability density function (PDF) method exhibit large edge fitting errors, this study proposes a regional inversion approach for the PDF method, which demonstrates a good overall fit. The objective is to reconstruct the fitting function for edge data, ensuring an accurate fit for both the edge and central regions.

For the edge data Q01-Q10 and Q31-Q40, a new fitting was performed with the fitting parameters q set to 0.02994, tanβ set to 0.5864, and the R-squared value at 0.988. The fitted curves from regional inversion are shown in Figure 6. The fitting residuals are shown in Table 3.

Based on the above analysis, it is recommended to adopt the probability density function (PDF) method, utilizing different predicted parameters for edge and central data. The inverted predicted parameters are presented in the following Table 4.

The calculations in this study were conducted using MIDAS GTS NX, a geotechnical engineering analysis software. To simulate the soil body, solid elements were employed due to their ability to accurately capture the three-dimensional stress-strain behavior of soil.

The model dimensions were set to 630 m (X) × 300 m (Y) × 351 m (Z), ensuring that boundaries were sufficiently distant from the area of interest to minimize boundary effects on simulation results. The soil material was represented using the Mohr-Coulomb model, chosen for its proven effectiveness in capturing the shear strength and failure characteristics of soil under the studied conditions. Key parameters of the Mohr-Coulomb model, including the internal friction angle (ϕ), cohesion (c), and unit weight (γ), were obtained from laboratory tests on site-collected soil samples, ensuring consistency with real-world conditions.

The boundary conditions were defined to reflect the physical constraints of the problem:

1. The upper boundary was set as a free surface, simulating an unconstrained ground surface.

To enhance reproducibility, the model assumes a homogeneous, isotropic soil body with no groundwater, neglecting fluid-solid coupling effects. The mesh was constructed using hexahedral solid elements with an average size of 5 m, balancing computational efficiency and accuracy.

The simulation was conducted in two stages:

1. Establishing an initial stress field under gravitational loading to reflect in situ conditions.

The model setup is illustrated in Figure 7, providing a visual representation of the geometry and boundary conditions.

Mining of the underground coal seam resulted in the loss of support for the overlying soil and rock mass, triggering a stress redistribution within the strata. The area adjacent to the coal pillar experienced significant pressure, leading to compressive failure of the soil and rock mass in that region (G.W. P., 2017). The simulation results are illustrated in Figure 8.

After excavation, stress redistribution induces uneven settlement within a certain range. The settlement value gradually decreases upward from the bottom of the mined-out area. On the surface, settlement decreases progressively from the center of the excavation toward both ends. The maximum surface settlement is 2.91 m.

Figures 9, 10 illustrate the fitting results of the probability density function (PDF) method and the probabilistic integration method, respectively, compared with numerical calculation results. The figures clearly depict the comparison between numerical results and fitting results.

A comparative analysis of the probabilistic integration method and the PDF method reveals that the probabilistic integration method may exhibit fitting accuracy issues within specific intervals. This phenomenon likely arises from its inherent limitations in global fitting capability. In contrast, the PDF method demonstrates a high degree of consistency with numerical simulation results in overall trends. This finding suggests that the PDF method offers a significant advantage in global fitting performance.