Machine Learning Prediction and Interpretability Analysis of Coal and Gas Outburst

1. Introduction

As the leading global producer and consumer of coal, ensuring mining safety is critical to the sustainable development of China's economy and energy sector [1,2]. Coal and gas outbursts represent a highly complex gas-dynamic hazard in underground coal mining. These phenomena are characterized by the rapid and violent ejection of substantial volumes of coal and gas from the seam within an extremely short period, posing severe threats to mining safety [3,4]. With the initial recorded outburst in France in 1834, more than 40,000 incidents have been recorded in over 20 countries. China has experienced over 13,000 outburst incidents to date, resulting in significant economic damage and casualties [5]. Figure.1 shows statistical data on outbursts in China from 2010 to 2024. Although the number of incidents and fatalities has decreased, the average number of deaths per outburst remains high. In January 2024, a coal and gas outburst at Pingmei No. 12 Mine resulted in 16 fatalities, highlighting the ongoing challenges in preventing and controlling outburst hazards.

Figure 1. Statistical data on coal and gas outbursts in China from 2010 to 2024.

Figure 1. Statistical data on coal and gas outbursts in China from 2010 to 2024.

The precise forecasting of outbursts is pivotal to against this hazard [6]. With rapid progress in computer technology, artificial neural networks, extreme learning machines, ensemble learning methods, hybrid and evolutionary algorithms, and other artificial intelligence techniques were employed to improve prediction performance [7,8]. For example, Wang [9] developed a chaos mapping and Lévy flight-improved crow search algorithm (ICSA) to optimize CNN, and constructed an outburst prediction model. Liu et al. [10] using support vector machine (SVM) enhanced by an improved snake optimizer (ISO), in which five refinement strategies were implemented to optimize the original SO algorithm. Nevertheless, most current studies have concentrated on improving model performance, ignoring the interpretability of the results. Many machine learning models lack transparency, exhibiting poor interpretability and making their decision-making logic difficult to comprehend [11]. To improve model explainability, interpretative methods are introduced, aiming to provide technical support and reference for preventing and controlling this hazard.

2. Materials and Methods

2.1. Analysis of Outburst Influencing Factors

The Pingdingshan mining area is a critical coal production base in China, with 156 outburst incidents recorded by the end of 2024. Figure.2 (a) illustrates the mining area in central Henan Province. Its structural framework is primarily controlled by the Qinling Orogenic Belt, characterized mainly by nearly EW trending thrust nappes parallel to the Qinling-Dabie tectonic zone [12]. Subjected to multiple phases of tectonic activity [13], the coal-bearing strata have been subjected to significant tectonic compression and shearing, resulting in significant displacement and deformation. The coal structure has been severely damaged, while the high compressive stress has effectively sealed gas within the coal seams. This regional tectonic evolution has provided the material basis and essential conditions: high tectonic stress, low-strength coal body, and enriched coal seam gas.

2.1.1. Geological Structure

The spatial distribution of outburst incidents exhibits distinct zoning characteristics. Figure.2 (b) illustrates the incidents are concentrated within four major outburst zones: Zone I (NW-trending fault-fold outburst zone): Composed of the Niuzhuang syncline, Guozhuang anticline, Niuzhuang thrust fault, and the former No.11 mine thrust fault (F1). This zone has experienced 16 outbursts; Zone II (Structural intersection zone): Located at the junction of the Xindian fault and Zhangwan fault, with 15 recorded outbursts; Zone III (Compound fold zone): Consisting of the Likou syncline and Guozhuang anticline, this zone has recorded 85 outbursts; Zone IV (uper zone of the Guodishan fault): This zone has recorded 13 outbursts.

Figure 2. (a) Location of Henan province and regional geological map; (b) Distribution of coal and gas outburst accidents.

Figure 2. (a) Location of Henan province and regional geological map; (b) Distribution of coal and gas outburst accidents.

All four zones are situated within intensely deformed tectonic belts characterized by concentrated tectonic stress, well-developed structural coal, and enriched coal seam gas. These geological conditions collectively establish these areas as high-risk zones.

2.1.2. Geological Structure

These incidents include 15 cases of eruption, 115 cases of extrusion, and 26 cases of outburst. These represent 9.61%, 73.72%, and 16.667% of the total incidents, respectively. As shown in Figure.3, extrusion constitutes the predominant type.

In-situ stress represents the superposition of gravitational and tectonic stresses [14]. The measured horizontal principal stress comprises tectonic stress and the horizontal component of gravitational stress, expressed as follows:

$${\sigma }_{h1}={\displaystyle \frac{v}{1-v}}\gamma D$$

(1)

$${\sigma }_{tH}={\sigma }_H-{\displaystyle \frac{v}{1-v}}\gamma D$$

(2)

where σ_h1 is the horizontal component of gravitational stress (MPa); σ_tH is the maximum horizontal tectonic stress (MPa); v denotes the Poisson’s ratio of the overlying rock mass; γ represents its unit weight (kN/m³); D corresponds to the burial depth (m).

Figure 3. Different types of coal and gas outbursts.

Figure 3. Different types of coal and gas outbursts.

Drawing on field in-situ stress measurements from 50 monitoring points in the Pingdingshan mining area [15,16,17], σ_tH was calculated using Equations (1) and (2). The strata overlying the coal seams are predominantly composed of sandstone and sandy mudstone. Therefore, this study assumes v = 0.25 and γ = 24 kN/m³. As illustrated in Table 1, it is found that σ_tH in eastern part (including Mines No.10, No.12, No.8, No.13 and Shoushan No.1) is significantly higher than that in western (Mines No.9, No.5, No.7, No.11) and central parts (Mines No.1, No.2, No.3, No.4, No.6). Correspondingly, the number of outburst incidents in the eastern part far exceeds those in the western and central parts. This indicates that tectonic stress shows a high correlation with outbursts.

2.1.3. Coal Structure

Tectonic coal is formed through the deformation or destruction of the primary coal structure by prolonged tectonic activity, resulting in a series of characteristics distinct from those of primary structured coal [18,19,20]. Tectonic coal exhibits lower cohesion and strength, making it prone to continuous fragmentation and ejection under in-situ stress and gas pressure [21,22,23].

D, E, F, and G coal seams are the main outburst-prone seams, having experienced 40, 41, 74, and 1 outburst incidents, respectively. Among these, the F coal seam has recorded the highest number of incidents, with outbursts occurring in this seam at Mines No.4, No.5, No.8, No.9, No.10, No.12, No.13, and Shoushan No.1. The distribution characteristics of tectonic coal in the F seam are complex, showing significant variations in coal seam damage types and the development thickness of tectonic coal across different areas. Figure.4 and Table 2 illustrates the distribution pattern of tectonic coal correlates well with the characteristics of the eastern mines, which experience a higher frequency and greater intensity of outbursts.

Figure 4. Amount of coal discharged on different coal mines.

Figure 4. Amount of coal discharged on different coal mines.

2.1.4. Coal Seam Gas

The maximum gas content data from D, E, and F coal seams in Mines No.9, No.5, No.6, No.4, No.1, No.2, No.10, No.12, No.8, No.13 and Shoushan No.1. along with the corresponding number of outburst incidents, were compiled and are presented in Figure.5 and Table 3. The eastern sector exhibits generally higher gas content compared to the central and western sectors, and accordingly records the highest frequency of outburst incidents. The maximum gas content on D, E, and F coal seams follows the order: F > E > D seam, which corresponds to the number of outburst incidents and the average outburst gas emission. The distribution pattern of gas content aligns with the frequency and intensity of outburst incidents.

Figure 5. Maximum gas content on D, E and F coal seams in different coal mines.

Figure 5. Maximum gas content on D, E and F coal seams in different coal mines.

The occurrence of outbursts is materially based on tectonic coal and gas [24,25,26,27], whereas in-situ stress provides the driving force [28,29]. Geological structures play the dominant controlling role [30,31,32], dictating the characteristics of tectonic coal, the gas reservoir conditions, and the in-situ stress field [33,34,35]. This relationship can be expressed as [36]:

$$F_{outburst}={\displaystyle \int (G,M,H,\sigma )}$$

(3)

where G, M, H, and σ represents the gas state, coal properties, geological structure, and in-situ stress, respectively.

Figure. 6 provides a visual representation of this formula. Correspondingly, the risk is heightened by greater gas content, higher in-situ stress, increased structural complexity, and reduced coal strength.

Figure 6. Maximum gas content on D, E and F coal seams in different coal mines.

Figure 6. Maximum gas content on D, E and F coal seams in different coal mines.

2.2. Model Selection

The selection of appropriate machine learning models is critical for accurate prediction. To comprehensively evaluate model performance, a range of classical and contemporary machine learning algorithms were selected for comparative analysis:

(1) K-Nearest Neighbors (KNN): This method is grounded in the principle that similar data points are presumed to share identical class labels [37].

(2) Back Propagation Neural Network (BPNN): This architecture employs a multilayer feedforward design, the parameters of which are iteratively optimized via the error backpropagation learning algorithm [38].

(3) Random Forest (RF): This method operates by constructing a multitude of decision trees during training and outputs the class that is the mode of the classes from the individual trees [39].

(4) Support Vector Machine (SVM): A powerful classifier that works by finding the hyperplane in a high-dimensional space that maximizes the margin between different data classes [40].

(5) eXtreme Gradient Boosting (XGBoost): An advanced implementation of gradient boosting ensemble method that sequentially optimizes weak learners to minimize a predefined loss function [41].

2.3. Bayesian Optimization (BO)

Robust tuning methods are essential [42]. BO is an efficient algorithm for hyperparameter tuning [43]. The BO employs an intelligent search strategy. It is capable of systematically traversing the hyperparameter space. Consequently, it can rapidly pinpoint hyperparameter configurations that maximize model performance.

The execution of BO primarily involves two key steps:

(1) A Gaussian Process (GP) is employed as a non-parametric surrogate model within the BO framework. Expressed as:

$$f\left(x\right)=gp\left[\mu \left(x\right),k\left(x,x^{\prime }\right)\right]$$

(4)

where μ(x) denotes mean function; k (x, x´) denotes covariance kernel function.

(2) The Expected Improvement (EI) acquisition function governed the selection of successive evaluation points during hyperparameter optimization, navigating the trade-off between venturing into unexplored regions of the parameter space and refining solutions near the current optimum. The formula for the EI function is given by:

$$EI(x)=E\left[\max \left(f(x)-f(x_{best}),0\right)\right]$$

(5)

2.4. Interpretability

The SHapley Additive exPlanations (SHAP) method was introduced to interpret the optimal model [44]. This method decomposes and attributes predictive outputs to individual input features via SHAP values, quantifying each feature’s marginal contribution across all instances. The model output g(z) can be expressed as the sum of the base value ϕ₀ (the mean prediction over all samples) and the SHAP value ϕ_i for each feature, expressed as:

$$g(z)={\varphi }_0+{\displaystyle \sum _{i=1}^M{\varphi }_i{z}_i}$$

(6)

where M denotes the feature count. z is an indicator function (taking values 0 or 1) denoting the presence of the feature.

The SHAP value is calculated as:

$${\varphi }_i={\displaystyle \sum _{S\subseteq \left\{X_1,X_2,\dots X_p\right\}\backslash \left\{X_i\right\}=1}^M\frac{\left|S\right|!\left(p-\left|S\right|-1\right)!}{p!}}\left[f\left(S\cup \left\{X_i\right\}\right)-f\left(S\right)\right]$$

(7)

where p denotes feature count; S represents a subset of features from the complete set {X₁, X₂, …, X_p}; f(S) denotes model output; $f\left(S\cup \left\{X_i\right\}\right)$denotes the model output after incorporating feature X_i; ${\displaystyle \frac{\left|S\right|!\left(p-\left|S\right|-1\right)!}{p!}}$ represents the weighting factor, which accounts for all possible permutations of the feature subset. For a fixed feature i, the number of possible combinations of the subset S is given by $\left|S\right|!\left(p-\left|S\right|-1\right)!$, considering the total number of permutations p! of all features.

2.5. Evaluation Index

Model performance was assessed via the Receiver Operating Characteristic (ROC) curve and its corresponding Area Under the Curve (AUC) [45]. According to established conventions [46], AUC values are interpreted as follows: a value below 0.5 suggests a model with no discriminative ability; 0.7–0.8 indicates acceptable; 0.8–0.9 indicates excellent; above 0.9 indicates outstanding. The key classification metrics referenced in this analysis are defined by the following formulas:

$$Accuracy={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

(8)

$$\mathrm{Pr}ecision={\displaystyle \frac{TP}{TP+FP}}$$

(9)

$$Recall={\displaystyle \frac{TP+TN}{TP+FN}}$$

(10)

$$F1Score=2\times {\displaystyle \frac{\mathrm{Pr}ecision\times Recall}{\mathrm{Pr}ecision+Recall}}$$

(11)

$$FPR={\displaystyle \frac{TP}{TN+FP}}$$

(12)

$$AUC={\displaystyle {\int }_0^1\mathrm{Re}call(FPR)dFPR}$$

(13)

3. Results

3.1. Data Description and Pre-Processing

As shown in Table 4, the original data dataset for coal and gas outbursts consists of 60 sample data points sourced from Reference [47], which originated from Pingdingshan No. 8 Coal Mine, China. The input features comprised multiple influencing factors. To evaluate feature relevance, the Pearson correlation coefficient served as the metric for gauging the predictive linear relationship of individual features to the target. A heatmap of the resulting correlation matrix is presented in Figure 7 for visual analysis.

X₁: Coal seam depth (m); X₂: Geological structure; X₃: Change of coal thickness (m); X₄: Soft layer thickness Variation; X₅: Change of coal seam dip angle; X₆: Coal seam thickness; X₇: Coal seam soft and fallen; X₈: Coal hardiness coefficients; X₉: Absolute gas emission volume (m³/min); X₁₀: Gas volume fraction (%); X₁₁: Initial gas desorption rate (cm³/g); Q: Amount of Coal Discharged (t); Y: Outburst risk level. The risk level Y is categorized as follows: Y = 1 (no risk) when Q = 0; Y = 2 (general risk) when 0 < Q ≤50; and Y = 3 (severe risk) when Q > 50.

The heatmap visualizes the Pearson correlation matrix, revealing significant linear associations between several predictive features and the target variable. Features exhibiting low correlation with the target variable were removal to simplify the model structure and enhance computational efficiency. Consequently, a total of six features, X₂, X₃, X₄, X₅, X₇ and X₈, were selected as input variables. To mitigate the influence of disparate feature scales, the dataset underwent normalization, constraining all variables to the interval [0, 1] before training commenced. The normalization is expressed as:

$$\left\{\begin{array}{l}{x}^{\prime }=\left(x_i-x_{min}\right)/\left(x_{max}-x_{min}\right);{\rho }_{xy}>0\\ x^{\prime }=\left(x_{max}-x_i\right)/\left(x_{max}-x_{min}\right);{\rho }_{xy}<0\end{array}\right.$$

(14)

3.2. Comparative Analysis of Model Accuracy

The comparative analysis of classification models employed standard evaluation metrics, where superior predictive capability is indicated by higher numerical values. A comprehensive summary of these metrics across all tested models is provided in Table 5. As shown, XGBoost achieved the highest test AUC (0.97) and Accuracy (0.90), followed closely by SVM.

3.3. Interpretability Analysis

To provide a global explanation for the model, the SHAP framework was adopted. This method operates by distributing the prediction output for each instance among the input features, proportionate to their computed marginal contributions over the entire dataset. This analysis is visualized in Figure 8, which combines a SHAP summary plot with a feature importance bar chart for the XGBoost model. In the summary plot, feature importance is ranked vertically, while the horizontal dispersion of points reflects the magnitude and direction (positive/negative) of SHAP values. Point color corresponds to the original feature value, scaled from low to high.

The results confirm that geological structure is a top-contributing feature, a finding supported by the wide dispersion of its SHAP values in Figure 8. The model's identification of geological structure as high impact supports its outburst relevance and interpretability. This aligns with practical observations. Numerous experts and scholars used methods such as field investigations and empirical analyses, demonstrated that geological structure is the dominant factor [31,32,48,49,50]. Geological structure is a critical determinant of outburst occurrence.

In addition, features such as the coal hardiness coefficients and change of coal thickness also exhibit relatively high SHAP values, serving as important secondary indicators of outburst risk. In contrast, features including coal seam soft and fallen, change of coal seam dip angle, and soft layer thickness variation were corresponded to low SHAP values, reflecting their marginal impact on the final predictions.

4. Discussion

In recent years, China has experienced frequent coal and gas outburst disasters. Meanwhile, many complex system challenges are increasingly being addressed relies on the advancement and implementation of intelligent computing methods, including multi-factor predictive modeling. However, current research places excessive emphasis on the models themselves, often improving performance by combining multiple algorithms. This results in overly complex models, making it challenging for coal mine decision-makers to extract key insights from them. In contrast, this study introduces interpretability methods. Beyond comparing the performance of various machine learning models, this study also provides interpretability for model predictions. This approach aids in identifying key factors influencing outbursts, thereby increasing the confidence of mine decision-makers in adopting such models.

While this study presents a preliminary framework for outburst risk prediction, several limitations remain and suggest directions for future research:

(1) Future work could employ more advanced feature engineering techniques, including graph-based feature interaction and automated feature selection, to better capture nonlinear relationships among predictors and further enhance model performance.

(2) Beyond the models evaluated here, other advanced architectures, including deep learning models such as Temporal Convolutional network (TCN) and Generative Adversarial network (GAN), could be investigated to offer new methodological insights for out-burst risk prediction.

5. Conclusions

By modeling the association between outburst drivers and outburst intensity, a data-informed predictive framework was constructed and rigorously evaluated, establishing a machine learning approach for outburst risk prediction.

(1) Incorporating an excessive number of indicators can degrade model performance. Analyzing and visualizing the correlations between influencing factors and the target variable via Pearson's correlation coefficient enables the identification of key predictors.

(2) A comparative evaluation framework was implemented in which all models were optimized via BO and assessed using ROC curve analysis. Among them, XGBoost delivered the most robust predictive performance.

(3) Analysis of SHAP values revealed geological structure occupies a crucial position in coal and gas outburst relevance and interpretability. coal hardiness coefficients and change of coal thickness also exhibit strong SHAP values, which are important secondary indicators outburst risk.