Abstract

Hydraulic support plays a key role in ground control of longwall mining. The smart prediction methods of support load are important for achieving intelligent mining. In this paper, the hydraulic support load data is decomposed into trend term, cycle term, and residual term, and it is found that the data has clear trend and period features, which can be called time series data. Based on the autoregression theory and weighted moving average method, the time series model is built to analyze the load data and predict its evolution trend, and the prediction accuracy of the sliding window model, ARIMA (Autoregressive Integrated Moving Average) model, and SARIMA (Seasonal Autoregressive Integrated Moving Average) model to the hydraulic support load under different parameters are evaluated, respectively. The results of single-point and multipoint prediction test with various sliding window values indicate that the sliding window method has no advantage in predicting the trend of the support load. The ARIMA model shows a better short-term trend prediction than the sliding window model. To some extent, increasing the length of the autoregressive term can improve the long-term prediction accuracy of the model, but it also increases the sensitivity of the model to support load fluctuation, and it is still difficult to predict the load trend in one support cycle. The SARIMA model has better prediction results than the sliding window model and the ARIMA model, which reveals the load evolution trend accurately during the whole support cycle. However, there are many external factors affecting the support load, such as overburden properties, hydraulic support moving speed, and worker’s operation. The smarter model of SARIMA considering these factors should be developed to be more suitable in predicting the hydraulic support load.

1. Introduction

With the rapid development of the Internet of Things, Cloud Computing, Big Data, Artificial Intelligence, and coal mine, integrating these emerging technologies will greatly change the traditional way of coal extraction [13]. Nowadays, the longwall mining face, integrating human-machine-environmental intelligent perception technique, big data analysis and decision-making method, and smart cooperation control technology for mining equipment, will improve the safety of underground workers significantly and reduce the number of underground operators and labor intensity. The hydraulic support is the main equipment of longwall mining face, and it plays an important role to ensure the safety of the working space. Therefore, the analysis, prediction, and prewarning of the hydraulic support load are the basis of achieving smart mining in the longwall face.

The evolution laws of mining-induced stress caused by coal seam excavation is directly related to the safety of workers and equipment in the working face, which has always been the hot spot and difficulty in the research of coal mine. Although many scholars have deeply analyzed the fracture structure of overburden strata of working face, the mechanism of load transfer, and the law of energy evolution, it is still difficult to carry out accurate quantitative calculation and predict the load change of hydraulic support in advance [412]. Some scholars had tried to use a data analysis method to analyze hydraulic support load [4], but the characteristics and the suitable analysis method for the load of hydraulic support is not clear. Due to the practical needs of engineering, it is necessary to predict the load of hydraulic support at least one cycle in advance, while many data analysis methods do not meet the requirements. With the increase of coal mining depth, the hydraulic support load prediction technology becomes more and more important, but there is no effective method to predict hydraulic support load in advance.

In the past, many scholars had studied the relationship between hydraulic support and surrounding rock using theoretical modeling. Qian and colleagues [5, 6] put forward the theory of “voussoir beam” and “key stratum,” which formed the foundational mechanics models for ground pressure and strata control in the longwall face. Song and Jiang [7] proposed the transfer-beam structure model and studied the relationship between hydraulic support and surrounding rock based on the given deformation and limit-deformation conditions. Wang and Pang [8, 9] studied the process of roof stratum breakdown instability and dynamic evolution of hydraulic support load, putting forward the “stiffness-strength-stability coupling model” between hydraulic support and surrounding rock, which provides an approach for dynamic analysis and prediction of hydraulic support load on longwall mining face. Many other scholars [1019] analyzes the fracture instability process of the overburden on working face in various mining conditions, proposes a method for calculating the suitable working resistance of the hydraulic support, and reveals the relationship between surrounding rock fracture instability and hydraulic support load.

With the continuous advancement of the working face, the hydraulic support load experiences increasing resistance, fluctuate resistance, and resistance relief in one cycle, and the peak load and load curve of each cycle are also affected by the dynamic fracture of the roof. The load curve, which can be called time series data, has an obvious trend and period. The existing hydraulic support load analysis methods just count the interval pressure and capture the initiation and peak support load, sometimes the resistance forces at the end of the support cycle are also monitored. Based on these data, the roof fracture step and the pressure in surrounding rock can be estimated roughly. However, it does not fully mine the information of the hydraulic support load data, and it cannot reflect the correlation of support load at different times. Therefore, a smarter analysis method needs to be developed for the advance prediction and early warning of the instability of surrounding rock.

At present, the analysis and prediction methods of time series data include two categories, which are listed as following [2022]: (1)time series models based on statistics(2)time series modeling using intelligent algorithms such as machine learning and deep learning (RNN, LSTM) analysis

It should be noted that intelligent algorithms such as machine learning and deep learning generally perform poorly on small sample data [23]. At the initial stage of cutting in longwall mining, only a few load samples can be obtained from hydraulic support. Therefore, this paper uses a time series model based on statistics to analyze and predict the load data of hydraulic support under the condition of small samples.

At present, the prediction methods for time series model based on statistics include the sliding window prediction method, exponential smoothing method, and autoregressive prediction method, which are widely applied in predict of economic development, climate change, energy demand, and other fields [2426]. The sliding window method is similar to the single exponential smoothing method, while the double exponential smoothing method and the triple exponential smoothing method are similar to the ARIMA method and the SARIMA method. In this paper, based on the feature decomposition of the support load data, the sliding window method, ARIMA method, and SARIMA method are applied to create the time series model of the hydraulic support load, respectively. By comparing the prediction results of varies models with responding parameters, and the adaptability of different models to the prediction of hydraulic support load, it obtains a more reasonable analysis and prediction method for hydraulic support load.

2. General Methodology

2.1. Engineering Background

With the rapid development of intelligent perception and network transmission technology, hydraulic support load monitoring technology and equipment have been mature [27]; the electrohydraulic control system can implement real-time monitoring of the entire hydraulic support load and upload, which provides data source for support load analysis and prediction.

The No. 121304 ultralength longwall mining face, in Kouzidong Coal mine, China, mines the 13-1 coal seam. The thickness of the coal seam is 2.2~6.66 m, the maximum inclination length of the cutting face is 350 m, and the strike length is about 1000 m. The ZZ13000/27/60D, which is a four-column shield hydraulic support, are applied in the cutting face, and the ultimate bearing pressure of the column is 33 MPa. The KJ216 system is introduced to monitor the pressure in the front column. The data of hydraulic support (No. 90) during two working days at the middle of cutting face is extracted for analysis. During the monitoring period, the cutting face advanced in a total of 10 cycles. The preprocessed data of the hydraulic support load is shown in Figure 1. The average cycle period has 37 data points. The load value is from zero to 35.4 MPa, and it has great fluctuation during initial working and support advancing.

The load of the hydraulic support experiences resistance increasing, fluctuation resistance, and pressure relief in one operating cycle, and the period of cycle and the support load value are various. As a whole, the sample data presents a regular trend with the characteristics of time series.

2.2. Feature Decomposition of Support Load and Predictability Analysis

Considering the period trend of hydraulic support load, the load data are decomposed based on statistics and data mining methods [28]. The trend term, cycle term, and residual term are extracted from the sample data, respectively, as shown in Figure 2.

It is found from Figure 2 that the trend term presents obvious oscillating and slight trend. The cycle term presents obvious periodicity, and a total of 10 cycles are extracted. The data of each period experiences obvious rising, fluctuation, and falling, which correspond to the process of increasing resistance, fluctuation resistance, and pressure relief of hydraulic supports. The extracted 10 cycles correspond to the cutting face advanced a total of 10 cycles. It means that the extracted cycles are well matching with the support cycle operation process. The residual term floats up and down around zero, which is in line with the characteristics of white noise, but the value of the residual term is quite large, indicating that the random fluctuation of the data is obvious. Overall, the hydraulic support load has strong time series characteristics, and the time series model can be used for prediction and analysis.

3. Analysis and Prediction of Support Load Based on Sliding Window

The sliding window method uses the average or weighted average value of historical monitoring data to predict the value at the current moment or any moment in the future. The equation can be expressed as following [29], where is the predicted value of the next moment, k is sliding window range, is weight of the th monitoring value within the sliding window, and is actual monitoring value of the t-nth value in the sliding window.

Based on the monitoring sample data, the initial 400 data are used as training data, and the last 60 data are used as verification data. The Python software is used to build a support load prediction model based on sliding window method. With different values of the sliding window, single-point prediction and multipoint prediction value based on average and weighted average value are calculated. Among them, single-point prediction just predicts the value at the next moment based on historical sliding window data, while multipoint prediction refers to predicting the value at multiple consecutive times in the future based on historical sliding window data. The comparison results are shown in Figure 3.

From Figure 3, it is found that the single-point prediction results based on the average and weighted average sliding window prediction method is in line with the verification data, and the prediction curve tends to be smooth with the value of the sliding window increasing, which means the prediction value at the peak or change-point becomes inaccurate. The prediction result using the weighted average method is more accurate than the average method, while the prediction result has a lag comparing with the monitoring data, and the lag accumulates with the sliding window value increasing. In Figure 3(b), the weights value of the sliding window with a value of 5 are [0.7, 0.1, 0.1, 0.05, 0.05]; [0.5, 0.2, 0.1, 0.1, 0.05], respectively. The weights value of the sliding window with a value of 10 are [0.7, 0.05, 0.05, 0.0325, 0.025, 0.025, 0.025,0.0125, 0.0125, 0.0125]; [0.5, 0.1, 0.1, 0.1, 0.05, 0.05, 0.025,0.025, 0.025, 0.025], respectively. With the same sliding window values, the greater weight of neighboring monitoring data, the more accurate prediction result, which means the single-point prediction result is determined by the neighboring monitoring values. Different sliding window values are used to predict 37 points of the sample data (one operating cycle of the hydraulic support), as shown in Figure 3(c), the results are shown at the right part of the curve. When the sliding window value is 5, it is accurate for predicting in short-term. With the predicting period extending, the accuracy reduces significantly, so as the other sliding window values.

4. Analysis and Prediction of Support Load Based on Autoregressive Model

4.1. Stationarity Analysis of Hydraulic Support Load Data

The autoregressive prediction method can establish the correlation between the historical monitoring data and the current and future data by decomposing the historical monitoring data into the trend term, seasonal term, noise term, etc., and then, based on historical data, it can calculate the future data using the autoregressive method. In general, the hydraulic support load just has the data character of monitoring pressure of the column, and the column pressure presents a cyclic change as the hydraulic support moves and the roof fractures. It is more suitable to use the autoregressive model for single variable predictive analysis.

The unit root test, self-correlation, and partial self-correlation analysis methods are used to verify the stationarity of the hydraulic support load data, and the original data, first-order difference data, and seasonal difference data are compared and analyzed to build the time series data modeling.

The autoregressive model requires that historical and current data have a strong self-correlation (self-correlation coefficient should not be less than 0.5), and the model can be described as follows [30], where is the current value, is the constant term, is the order, is the self-correlation coefficient, is the error term, and is the historical monitoring values at time before time .

The correlation between the monitoring value at the current time and the history monitoring value can be calculated using a self-correlation function [31], where is the correlation coefficient between the monitoring value at time and the current monitoring value, is the current value, and is the historical monitoring values at time before time .

Autoregressive models require time series data to be stationary, which means the mean value and variance do not change with time. The unit root test (ADF) is used to test the stationarity of the sample data of the hydraulic support load. The calculation results are shown in Table 1.

From Table 1, it is found that the value of the preprocessed hydraulic support load sample data tends to 0, and the value is significantly less than the confidence intervals values with the 99%, 95%, and 90%, which indicates that the data is steady. The autocorrelation function (ACF) and partial autocorrelation function (PACF) are used to analyze the stationarity and self-correlation of the data, as shown in Figure 4, the shaded area covers the confidence interval.

From Figure 4, it is found that both the self-correlation coefficient and partial self-correlation coefficient of the hydraulic support load sample data decline slowly and tend to be flat (fluctuating around zero value) eventually. The self-correlation coefficient presents periodic cyclic fluctuations within the confidence interval (shaded area in Figure 4(a)), while there are many obvious lag points (points outside the shaded area in Figure 4(a)). The partial self-correlation coefficient also shows up and down during the confidence interval (shaded area in Figure 4(b)), but the periodicity is not obvious. From Figure 1, we know that the hydraulic support load in cyclic operation is affected by roof fracture behavior, cutting face advance speed, the operating state of hydraulic support, and the habits of the operators. The support load change trend in each cycle still shows a large difference, so did the length of each cycle. Therefore, the periodicity of support load is similar to the “seasonal” of the traditional time series model, while the length of the “season” and the change rule of sample value in the seasonal cycle are various. In this paper, the ARIMA model (not considering the influence of the cycle period) and the SARIMA model (considering the influence of the cycle period) are used to analyze and predict the load sample data of the hydraulic support.

Although the hydraulic support load monitoring samples have passed the stationary test from Figure 4, there are still many lag points in the ACF and PACF. The first order difference of sample data is introduced to reduce the number of lag points, so that the autocorrelation coefficient and partial autocorrelation coefficient converge to the confidence interval quickly. According to the seasonality of the sample data and the cycle periodicity data extracted in Figure 2(b), 37 monitoring values are chosen in the cycle period, and the data is subjected to “seasonal” difference to obtain stationary time series data of hydraulic support load samples.

The unit root test method is used to perform stationarity test on the preprocessed data. The calculation results are shown in Table 2. The self-correlation coefficients and partial self-correlation coefficients after data processing are shown in Figure 5.

The unit root test results show that with the first-order difference or the seasonal difference of the sample data, the test results are significantly less than the results of original data and the test value of the confidence interval; the value is also closer to zero than the original data, which indicates that the sample data become more steady with the first-order difference or seasonal difference.

From Figure 5, it is found that the lag of the self-correlation coefficient and partial self-correlation coefficient is significantly improved (the points outside the shaded area in Figure 5). In Figures 5(c) and 5(d), the partial autocorrelation coefficients at and increase significantly, and the partial autocorrelation coefficients at are greater than the value at , reflecting the periodicity of the data. It also indicates that the cyclical change characteristics of recent data have a greater influence on current and future data.

4.2. Support Load Prediction Based on ARIMA Model

Considering the variety of the hydraulic support load cycle period, the “seasonal” factor was suspended in the first model. Autoregressive moving average modeling analysis (ARIMA) is firstly applied for the hydraulic support load sample data with first-order difference. From Figures 5(a) and 5(b), it is found that the partial self-correlation of the data appears “censoring” with 5 of lags, and the self-correlation appears “tailing,” so the autoregressive value of the model is determined as . Similarly, the self-correlation of the model also appears “censored” when , and the partial self-correlation shows “tailing,” the moving average term of the model is initially determined to be 5.

According to the determined model autoregressive terms and moving average terms, the grid search method is used to create 16 models in total under the condition that the value range of the autoregressive term is , the value range of the moving average term is , and the value of the differential term is . Akaike Information Criterion (AIC) and Bayes Information Criterion (BIC) are used to optimize the model parameters [32]. The calculation results of AIC and BIC for these models are shown in Figure 6.

It can be seen that, when and , the AIC value is the smallest (2456.925), and when and , the BIC value is the smallest (2488.379), and when and , the BIC value is slightly larger than the result when and . Therefore, according to the traditional time series model parameter optimization method, the optimal parameter of the model is determined to be , , and . The residuals of the model, as shown in Figure 7, basically conform to the characteristics of white noise, which means the high reliability of the model.

Based on the above results, the training data set is used for data modeling of the hydraulic support load, and the model is used to predict the hydraulic support test samples (the next operating cycle data). The result is shown in Figure 8.

From Figure 8, it is found that the prediction result of the training sample has the similar trend and peak value with the sample data, but the prediction result of the validation data just presents the similar trend and sample value in a very short stage, and then it shows a similar horizontal straight line trend (similar to the sliding window method). Therefore, the model can and can only predict the development trend of the data well in short term, which has a slight advantage comparing with the sliding window prediction method.

The parameters of the model are , , and , that means, the previous four data before the current time are used to predict the future trend, and the previous three data before the current time are used to judge the error of the sample, the trend term of the reference sample can predict the short-term development trend of the support load, but the small number of reference samples is just suitable to predict the development trend in short term, while cannot predict the long-term development trend of the support load well. In order to make full use of historical data, the historical sample data of one operating cycle is used for data modeling and analysis, and the model parameters are modified to , , and . The results are shown in Figure 9.

From Figure 9, it is found that with the number of the historical reference data increasing, the model becomes more sensitive to data fluctuations, the predicted value change significantly with slight sample data fluctuating. However, the prediction result is more accurate. Although the model’s residual value still fluctuates around zero, the fluctuation range enlarges significantly. It may be because the value increases resulting in an increasing of the model’s cumulative error value. To further optimize the model parameters, reducing the value to 3, the residual value and the prediction result are shown in Figure 10.

From Figure 10(a), we can know that reducing the value can greatly weaken the fluctuation of the model’s residual value, as well as the fluctuation of the model’s predicted value. The prediction results have been significantly improved. The predicted value is basically consistent with the trend of the sample verification data set, while there has slight difference at the peak value. The reason is that there are many external interfering factors for the hydraulic support during the operation cycle, the fluctuation range of the residual value is -30~ 25 MPa, which leads to the deviation of the predicted value and the actual value. The prediction results can be improved with an appropriately larger value, while the exact results still do not cover the entire hydraulic support cycle. The long-term prediction results are still no reference value.

4.3. Hydraulic Support Load Prediction Based on SARIMA Model

Through the feature decomposition of the hydraulic support load data (Figure 2), it is found that the sample data has some cycle change characteristics, while the cycle period and the load value are various. Considering the periodicity, the Seasonal Autoregressive Integrated Moving Average (SARIMA) is applied to analyze and predict the load of the hydraulic support, and the seasonal difference processed data is used for modeling and analysis. The self-correlation and partial self-correlation analysis results are shown in Figures 5(c) and 5(d). It can be seen that the partial self-correlation of the sample data appears “censored” at , while the self-correlation presents “tailing.” The order of the nonseasonal autoregressive component of the model can be primarily determined to be ; and at , , the partial self-correlation is significantly enhanced, so the seasonal component order () of the autoregressive model is 2. The self-correlation of the model also appears “censored” at , and the partial self-correlation appears “tailed,” so the model’s nonseasonal moving average value is determined to 6 primarily, the self-correlation of the model enhances significantly at , so the order of the seasonal component in the moving average model is determined to 1 primarily.

Based on the primarily determined model parameters, the grid search method is used to create 630 models in total under the condition that the range of nonseasonal autoregressive terms is , the range of nonseasonal moving average terms is , and the nonseasonal difference term is , the range of seasonal autoregressive terms is , the range of seasonal moving average terms is , the seasonal difference term is , and the period length is . Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC) are also used to optimize the model parameters. When , , , , , and , the model has the smallest AIC value (1285.72), when , , , , , and , the model has the smallest BIC value (1341.07), and when , , , , , and , the BIC value (1345.62) is just slightly larger than the parameter with , , , , , and . The optimal parameters of the model can be determined as , , , , , and . The residuals of the model basically conform to the characteristics of white noise, which means the model is available. The data of the last operating cycle is used as a sample test set, and the prediction results based on the optimal parameters are shown in Figure 11.

From Figure 11, it is found that this model is also very sensitive to the fluctuation of the sample data. The predicted value fluctuates greatly with the sample data fluctuating. The residuals of the model also fluctuate around zero, but the fluctuation range is larger than that of the ARIMA model. It is inferred that the wider fluctuation range of the residual value leads to the more sensitive model to the fluctuation of the sample data. The model predicts that the load of the hydraulic support will decrease rapidly and increase rapidly in the short term, and it also predicts the load change trend in the next hydraulic support working cycle very well. The predicted trend is very close to the field trend; the predicted value is slightly lower than the monitoring value. The predicted peak value is close to the monitoring value, with a significant lag. There are many factors affecting the load of hydraulic supports, so the residual range of the data samples is wide, resulting the fluctuation of the data. Despite the low prediction precision of this model, it predicts the development trend of the hydraulic support in the next entire cycle. The prediction accuracy is affected by the roof properties, support working condition, and other external factors, which has certain randomness, so we do not take more attention on adjusting parameter to optimize the prediction result.

5. Discussion

The variance of the hydraulic support load reflects the roof fracturing process. However, the hydraulic support load is determined by multiple factors, such as roof fracture, support working condition, working face advancing speed, and the operator’s capabilities. As a result, the cycle length of the hydraulic support load varies with the condition changing, as well as the peak load, which is not good for predicting the hydraulic support load.

The sliding window prediction method has a good short-term prediction effect on the support load, but the value of the sliding window should not be too large. The weighted mean method is significantly better than the mean method. Due to the heavy dependence of the sliding window prediction method on the recent monitoring data, it is not possible to judge the trend of the data in advance, so it is difficult to predict the medium-term and long-term law of the hydraulic support load changes. Only predicting the short-term load of hydraulic supports (less than one cycle of support cycle operation) has great limitations in practical engineering applications.

By comparing the prediction results based on different model parameters in Table 3, as shown in Figure 12, it is found that both the ARIMA model and the SARIMA model can predict the load trend. However, the ARIMA model is only suitable for predicting the trend in a short period. The ARIMA model using the parameters of three groups just obtain good result at the first half of the hydraulic support cycle, while fails at the remaining period (unloading process), and the peak load of the hydraulic support cannot be calculated. The SARIMA model with the suitable parameters can obtain the support load development trend in one working cycle, and the predicted peak load is close to the field monitoring data. It indicates that the SARIMA model is more suitable for predicting the support load development trend. It should be noted that the results of the SARIMA model has a lag comparing with the monitoring time series data. More external factors need to be considered in the SARIMA model in future study.

The hydraulic support load is often affected by external factors. To identify the prediction accuracy of these models, it should be taken more attention in geological condition, recovery method, and so on, which need a new study to present it. In this paper, we mainly focus on comparing the advantage of the three models. As we know, the roof fracturing is a dynamic continuous process. In the condition of good support effect on roof and normal advancing speed of mining face, the SARIMA model can obtain a good prediction result of the hydraulic support load. In addition, the SARIMA model is more sensitive to the fluctuation of the monitoring data, and it will amplify the abnormal value slightly, as shown in Figure 11, the maximum and minimum values of the predicted result are significantly higher or smaller than the monitoring value. It indicates the SARIMA model can capture the abnormal monitoring value more effectively.

The time series model based on statistical analysis is more suitable to perform data model and prediction on small sample data. If the monitored values that are less affected by the external factors are used as training samples, the load of the hydraulic support in the subsequent cycle can be better predicted. The predicted results can also be used as the basis for judging whether the hydraulic support is affected by external factors in the next operating cycle. If it is affected, the value will deviate obviously; otherwise, the prediction and monitoring results are basic identical.

Obviously, the time series models based on statistical theory have limitations to predict the hydraulic support load. It is difficult for these models to predict large-scale pressure accidents and rock burst accidents of cutting face occurring instantaneous. Therefore, the mechanical model, numerical simulation results of the roof fracture on cutting face should be coupled in the time series model in future study, so that it improves the prediction accuracy.

6. Conclusions

(1)By decomposing the load data of the hydraulic support into trend term, cycle term, and residual term, the periodicity load in the hydraulic support cycle can be seen obviously. The trend term of the load data presents a certain rule, but the residual value is quite large, which means that the hydraulic support load fluctuates greatly with the various external factors, and it increases the difficulty to predict support load(2)The sliding window prediction method is good at predicting the next single-point. However, it fails to predict the development trend of load data. With the value of the sliding window increasing, the prediction results of the peak point or the abrupt point become inaccurate. The value of the sliding window should be in a reasonable range. The weighted average method is significantly better than the average method(3)Increasing the value of the autoregressive term of the ARIMA model can properly improve the long-term prediction result, while it also increases the sensitivity of the model to data, and the model still fails to predict the load trend in the next support cycle. Therefore, this model is not suitable to predict the load trend in a long-term mining process(4)The prediction data generated by the SARIMA model has the similar trend with the field monitoring data, which means this model is the better approach to predict the load trend in a hydraulic support cycle. The load data has the strong periodicity rather than tendency. However, the predicted value is slightly lower than the monitoring data, and it exists time lag effect

Data Availability

The No. 121304 ultralength longwall mining face, in Kouzidong Coal mine, China, mines the 13-1 coal seam. The thickness of the coal seam is 2.2~6.66 m, the maximum inclination length of the cutting face is 350 m, and the strike length is about 1000 m. The ZZ13000/27/60D, which is four-column shield hydraulic support, are applied in the cutting face, and the ultimate bearing pressure of the column is 33 MPa. The KJ216 system is introduced to monitor the pressure in the front column. The data of hydraulic support (No. 90) during two working days at the middle of cutting face is extracted for analysis. During the monitoring period, the cutting face advanced in total of 10 cycles.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The study was supported by the National Key Research and Development Program (2017YFC0603005) and by the National Natural Science Foundation of China (52004124). The authors express their gratitude to the Kouzidong Coal Mine for providing the field measurements and observations used in this work.