Coal and gas outbursts are influenced by four main factors: geological conditions, coal seam characteristics, gas-related factors, and operational practices. Based on field observations, these factors are further subdivided into 14 specific elements (He et al., 2010) (Figure 1). The risk level (L) of coal and gas outbursts is classified into five categories based on the amount of ejected coal (Table 1).

In this study, we used coal and gas outburst data collected from a coal mine in Pingdingshan over a period spanning from 1984 to 2009 (Xie et al., 2018). The first 50 data points were selected as the training set, while the final 10 data points were reserved for testing purposes. The coal and gas outburst data from the Pingdingshan Coal Mine are summarized in Table 2. The geographical location of Pingdingshan is shown in Figure 2.

Before feeding the training data into the model, it is crucial to carefully preprocess the dataset by identifying and addressing any anomalies.

Step 1: The Pauta criterion was applied to establish the outlier threshold. Data points that deviated from the mean by more than three standard deviations were considered outliers. These outliers were then removed, and their corresponding entries were set to null values.

Step 2: Missing values were then filled by extracting the five data points preceding and following each missing entry. In this study, Lagrange interpolation was used to estimate the missing data, as shown in Equations 1, 2.

where l_i(x) represents the interpolation basis function; L_n(x) denotes the interpolated value of the missing data; y_i is a known (non-missing) value; x is the index corresponding to the missing value; x_i denotes the index of the known data point; and x_j is the interpolation node.

The primary indicator for determining the severity of a coal and gas outburst was the quantity of coal thrown (Table 1). As shown in Figure 3, Pearson's correlation analysis was performed to rigorously investigate both the interrelationships among the various influencing factors and their individual correlations with the quantity of coal thrown, as defined by Equation 3.

where cov(X, Y) is covariance and σ_X and σ_Y are the standard deviations of X and Y, respectively.

Following the strong and weak correlation partitions established by Zhang et al. (2022), the correlation results presented in Figure 3 were subsequently classified. The outcomes of this classification are comprehensively detailed in Table 3. As demonstrated in Table 3, this study meticulously selected six factors exhibiting medium to high correlation levels as input variables for the model, specifically including A5, B3, D2, A4, B2, and A3.

XGBoost is a prominent and highly efficient boosting ensemble learning algorithm, representing an advanced evolution of the Gradient Boosted Decision Tree (GBDT) model (Yao et al., 2022; Xiong et al., 2024; Utkarsh, 2024). The predictive output of the XGBoost model is formulated as shown in Equation 4:

where ŷ_i denotes the predicted value for the i-th sample; K represents the number of trees; F signifies the function space of the tree; x_i is the feature vector of the i-th data point; and f_k refers to the function learned by the k-th tree, which is characterized by its structure q and leaf weights w.

The loss function of the XGBoost model comprises two components, as shown in Equation (5):

where the first term represents the training error between the predicted value ŷ_i and the true target value y_i; the second term denotes the sum of tree complexities, which serves as a regularization term to control the model's complexity, as presented in Equation 6:

where γ and λ are the penalty coefficients.

During the minimization process of the objective function defined in Equation 5, the incremental function f_t(x_i) is added at each iteration to reduce the loss function. The objective function at the t-th iteration is presented in Equation 7:

For Equation 7, the The sample set is defined in each leaf of t objective function is approximated using a second-order Taylor expansion. The j-th tree as I_j = {i|q(x_i = j)}. Here, gi=∂y^i(t-1)l(yi,y^i(t-1)) and hi=∂y^i(t-1)2l(yi,y^i(t-1)) represent the first and second derivatives of the loss function, respectively. From these definitions, Equation 8 can be derived as follows:

Defining Gj=∑i∈Ijgi and Hj=∑i∈Ijhi leads to Equation 9:

The partial derivative with respect to w_j yields Equation 10:

Substituting the weights into the objective function yields Equation 11:

A smaller loss function signifies enhanced model performance. A greedy algorithm is used to partition the subtree by enumerating feasible split points: a new split is added to existing leaves at each step, and the maximum gain is computed accordingly. The gain is then calculated as shown in Equation 12.

where the first and second terms represent the gains achieved by splitting the left and right subtrees, respectively, whereas the third term corresponds to the gain obtained without any split.

The attention mechanism receives the output from the Stacking model as its input and adaptively assigns weights to its input features, thereby emphasizing the most relevant ones and suppressing less important features, consequently facilitating more accurate feature selection (Zhu et al., 2021; Wankhade et al., 2023). The structure of the attention mechanism is shown in Figure 4.

In Figure 4, x₁, x₂, …, x_t represents the inputs from multi-source fusion data; h₁, h₂, …, h_t corresponds to the output state values generated by the ensemble model, while a_ti signifies the adaptive attention weight assigned to each output; s_t denotes the final output.

By calculating the correlation between h₁, h₂, ⋯ , h_t and the current decoding time, the e_{t, i} for each influencing factor is obtained. These updated values are presented in Equation 13.

According to the computed probability e_{t, i} of each influencing factor within the population, this value is then used to compute the attention weight for each output of the ensemble model. The updated output is presented in Equation 14.

The hidden states h₁, h₂, ⋯ , h_t are weighted by their corresponding attention values and then linearly combined. The updated output is shown in Equation 15.

S_t represents the final output derived through the attention mechanism, as shown in Equation 16.

where V, W, and U represent the trainable parameters, which are continuously updated during model training.

The performance of the Stacking model is directly influenced by the number of base models used. Using too few base models may not provide adequate diversity to effectively support the meta-model, while using too many could result in redundancy, higher computational costs, and a more intricate tuning process. Typically, 3–5 base models are recommended (Kumar et al., 2024; Zhao et al., 2024b).

Based on the predictive capabilities of various base learners, this article selected high-performing models as the first-layer base learners in the Stacking model. This selection is driven by the fact that base models with strong learning abilities contribute to improving the overall predictive accuracy of the ensemble. Specifically, RF, which uses the bagging technique, is preferred for its robust learning capacity and well-established theoretical foundation, making it applicable across a wide range of domains. SVM is selected for its unique strengths in handling small datasets, non-linear relationships, and high-dimensional regression problems. KNN is included due to its solid theoretical background and efficient training process, delivering strong practical performance. For the second layer, models with robust generalization capabilities are chosen to aggregate and correct biases from the multiple base learners in the training set while mitigating overfitting through ensemble strategies. Consequently, the Stacking ensemble model incorporates RF, KNN, and SVM as the first-layer base learners, with Attention-XGBoost serving as the meta-learner in the second layer. The overall architecture is shown in Figure 5.

It is important to note that the training set for the meta-learner is derived from the outputs of the base learners. Directly utilizing the base learners' training data to form the meta-learner's training set could result in significant overfitting. To prevent the data from being redundantly learned by both layers and to avoid overfitting, an appropriate data usage strategy must be implemented. The dataset is first split into training and testing sets using cross-validation, with the three base learners making independent predictions. For each base learner, the original training dataset is partitioned into six mutually exclusive subsets, ensuring that no data IDs are repeated across subsets. For each base learner, one data subset is reserved as the validation set, while the remaining five subsets serve as the training set. Each base learner produces prediction results on its own validation subset. These predictions from the three base learners are then combined to form a new dataset, equal in size to the original dataset, as illustrated in Figure 6. This approach facilitates a comprehensive feature transformation from the original input features to the meta-learner's input features. Since each base learner's predicted data subset was excluded from its own training, this method guarantees that every data point is used only once during training, effectively preventing overfitting.

The training and prediction process of the AXGBoost-Stacking model is shown in Figure 6, and the detailed training procedure is outlined as follows:

Step 1: The coal and gas outburst dataset is defined as presented in Equation 17.

where x_n represents the feature vector of the n-th sample, y_n is the corresponding target (predicted) value, and p is the number of features, meaning that each feature vector can be expressed as (x₁, x₂,..., x_p). Next, the dataset is partitioned into Z equally sized subsets: S₁, S₂,..., S_z. The cross-validation between datasets as presented in Equation 18.

where S_z denotes the z-th test set and S¯z represents the corresponding training set.

Step 2: The training set S̄z is fed into the first layer of the XGBoost-Stacking ensemble model, where three base learners are trained to obtain the base model L. Simultaneously, each sample x_n in the cross-validation test set S_z is passed through the trained base model L to generate the corresponding predictions.

Step 3: The output predictions from the three base learners are concatenated to form a new data sample, which is then used as the input for the second layer of the Stacking model. At this stage, a prediction algorithm that integrates the attention mechanism with XGBoost is used to aggregate these outputs and finalize the prediction of coal ejection volume.

In this study, the AXGBoost-Stacking model is implemented using the scikit-learn library in Python. A detailed description of the algorithm is provided in Table 4.

The multi-model Stacking prediction framework proposed in this study adopts the AXGBoost-Stacking model, which uses SVM, RF, and KNN as base learners and an attention-enhanced XGBoost model as the meta-learner. The Stacking ensemble learning algorithm enables a two-layer fusion of the SVM, RF, KNN, and XGBoost models. In addition to AXGBoost-Stacking model, three alternative Stacking models can also be constructed for comparative analysis:

To evaluate the predictive performance of the AXGBoost-Stacking model and compare it with the individual predictive capabilities of the other three Stacking models, this study uses mean squared error (MSE), mean error (ME), and the F1-score [as defined in reference Xie et al. (2018)] as evaluation metrics. The formulas for calculating MSE and ME are provided in Equations 19, 20.

where y_i is true data, and ŷ_i is prediction data.

The calculation formula for the F1-score is presented in Equation 21.

where Precision indicates the model's accuracy, and R_ecall represents the model's recall rate.