葉黎明 施式亮，2教授 魯 義教授 李 賀副教授 曾明圣

(1.湖南科技大學 資源環(huán)境與安全工程學院，湖南 湘潭 411100；2.湖南科技大學煤礦安全開采技術湖南省重點實驗室，湖南 湘潭411201)

0 引言

進入21世紀，隨著我國煤礦事故死亡人數(shù)進入穩(wěn)定下降時期[1]，煤礦安全生產(chǎn)形勢逐年向好。2020年我國煤礦事故死亡人數(shù)為228人、煤礦百萬噸死亡率0.058、瓦斯事故死亡人數(shù)30人，死亡人數(shù)與1989年的歷史最高值7 625人相比下降了97.0%，百萬噸死亡率與最高時1949年的22.541相比，下降了99.7%，瓦斯事故死亡人數(shù)僅為1997年歷史峰值3 800人的0.79%[2-3]。從以上數(shù)據(jù)可以看出，我國煤礦安全生產(chǎn)形勢已經(jīng)明顯改善，但仍應認識到同世界頂尖水平之間的差距。根據(jù)美國礦山安全和健康管理局(Mine Safety and Health Administration， MSHA)公布的數(shù)據(jù)，2020年美國共發(fā)生各類煤礦事故29起，共造成5人死亡，百萬噸死亡率0.001，分別為同期我國的2%與1.7%。因此分析和預測我國煤礦安全態(tài)勢，有助于明確煤礦安全態(tài)勢的發(fā)展規(guī)律，對促進煤礦安全生產(chǎn)具有重要意義。

煤礦生產(chǎn)系統(tǒng)的安全性受到諸如煤炭資源賦存狀態(tài)、煤礦開采技術條件、安全投入以及多種致災因子等系統(tǒng)自身或外部因素影響，因此煤礦生產(chǎn)過程中事故的發(fā)生具有多因素時空耦合的復雜性，進而導致表征煤礦安全態(tài)勢各項指標呈現(xiàn)出較強的非線性特征，僅以單一線性模型刻畫重要指標隨時間的變化規(guī)律存在較大的局限性。針對各項表征煤礦安全態(tài)勢指標的預測，常用的方法有R/S分析法[4]、灰色預測法[5]、神經(jīng)網(wǎng)絡法[6]、指數(shù)平滑法[7]等非線性預測方法。時間序列預測法作為一種重要的非線性預測方法，在道路交通[8]、民航[9-10]、工業(yè)生產(chǎn)[11]等領域的安全態(tài)勢指標預測中得到廣泛應用，而運用此方法對煤礦安全態(tài)勢指標進行預測的研究較少。時間序列預測法采用差分自回歸移動平均(Autoregressive Integrated Moving Average Model， ARIMA)、Holt-Winters等模型進行預測，鑒于各項煤礦安全態(tài)勢指標中所包含的線性、非線性數(shù)據(jù)將影響單一時間序列模型的預測精度。筆者運用組合預測思想[12]，改善單一模型預測效果。

本文從煤礦安全態(tài)勢指標時間序列組合預測角度出發(fā)，運用ARIMA模型提取各態(tài)勢指標時間序列的非線性主部，在此基礎上，進一步采用極端梯度提升算法(eXtreme Gradient Boosting， XGBoost)，利用該算法預測非線性數(shù)據(jù)時不易過擬合、泛化能力強的特性對ARIMA模型殘差進行修正，從而提高預測精度?；诖?，本文建立ARIMA-XGBoost的煤礦安全態(tài)勢指標預測模型，以期為煤礦安全生產(chǎn)態(tài)勢指標的預測提供一種新思路。

1 混合模型的構(gòu)建

1.1 ARIMA模型

ARIMA模型是由Box和Jenkins[13]提出的一種時間序列預測模型，模型包含參數(shù)p、d、q，其中p為自回歸階數(shù)，d為原始序列平穩(wěn)化所需的差分階數(shù)，q為移動平均階數(shù)，亦可寫作ARIMA(p，d，q)模型。該模型可以通過歷史數(shù)據(jù)對未來進行預測[14]。其建模的基本步驟是根據(jù)單位根檢驗的結(jié)果，確定原始序列平穩(wěn)化所需的差分階數(shù)d，通過分析自相關與偏自相關函數(shù)的截尾、拖尾特征，確定參數(shù)p、q。模型中，AR是自回歸如式(1)所示，MA是移動平均如式(2)所示。

AR(p):Yt=μ+α1Yt-1+α2Yt-2+···+αpYt-p+εt

(1)

MA(q):Yt=μ+β1εt-1+β2εt-2+···+βqεt-q+εt

(2)

將公式(1)與公式(2)結(jié)合，再進行d階差分得到ARIMA(p，d，q)模型，其數(shù)學描述可簡記為：

式中：

Yt—時間序列中第t項的觀測值，t為當前項在時間序列中對應的序號；

ΔdYt—Yt經(jīng)過d階差分后得到的觀測值；

αi、βi—時間序列中第i項的帶估計參數(shù)，i是循環(huán)變量；

p、q—模型的階數(shù)；

ε—預測殘差；

μ—常數(shù)。

1.2 XGBoost算法

XGBoost算法是學者陳天奇提出的一種改進的提升決策樹(Gradient Boosting Decision Tree，GBDT)算法[15]。它對傳統(tǒng)GBDT算法的革新主要體現(xiàn)在2方面，一是XGBoost算法的最終結(jié)果是由多棵決策樹通過不斷的迭代計算之后，再累加而得，在此過程中能充分利用CPU的多線程處理性能，大幅度提升計算速度；二是XGBoost算法通過在目標函數(shù)中加入正則化項的方式，簡化模型，避免過擬合。以下是它的常規(guī)建模過程：

(4)

F={f(x)=ωz(x)}

式中：

xi—第i個數(shù)據(jù)點的特征向量；

yi—第i個數(shù)據(jù)點的真實值；

K—決策樹數(shù)量；

fk—對應決策樹空間中的一個函數(shù)；

F—回歸樹的集合空間；

z—特征向量xi映射到樹的葉子節(jié)點；

ωz(x)—單棵樹的預測值。

目標函數(shù)包含2部分：

(5)

(6)

此時目標函數(shù)可描述為：

(7)

式中：

C—為常數(shù)項。

當損失函數(shù)l為非平方誤差函數(shù)時，求解最優(yōu)解的過程復雜，為簡化計算，去掉常數(shù)項并采用泰勒展開來近似定義目標函數(shù)。

(8)

此時優(yōu)化后的目標函數(shù)：

(9)

其中，Gj=∑i∈Ijgi，Hj=∑i∈Ijhi；Ij為第j棵樹每一葉子中的樣本集合；目標函數(shù)的值越小則樹的結(jié)構(gòu)越好[16-17]。

1.3 ARIMA-XGBoost混合預測模型構(gòu)建

(10)

根據(jù)ARIMA模型的預測值和XGBoost殘差修正值得到混合模型預測值，即

(11)

2 實證計算與分析

2.1 數(shù)據(jù)來源與處理

本文數(shù)據(jù)來源于《中國煤炭工業(yè)年鑒》[2]與《中國安全生產(chǎn)年鑒》[3]。

選取1949-2020年煤礦事故死亡人數(shù)(以下統(tǒng)稱死亡人數(shù))與1966-2020年煤礦百萬噸死亡率和瓦斯事故死亡人數(shù)的統(tǒng)計數(shù)據(jù)作為研究對象，以年作為時間步長構(gòu)建“死亡人數(shù)”“百萬噸死亡率”“瓦斯事故死亡人數(shù)”原始序列，分別用{Ytoll}、{Yrate}和{YGtoll}來表示，如圖1。

2.2 ARIMA模型預測

2.2.1 差分階數(shù)的確定

根據(jù)圖1中的折線變化趨勢，初步判斷3組原始序列均不符合零均值同方差的特性。3組原始序列單位根檢驗結(jié)果，見表1。由表1可知，{Ytoll}、{Yrate}和{YGtoll}的檢驗τ統(tǒng)計量分別為-1.098 92、-1.036 28和-1.251 304，均大于各自顯著水平10%的臨界值，故序列{Ytoll}、{Yrate}和{YGtoll}皆為非平穩(wěn)序列，需進行差分處理。對3組原始序列一階差分后得到新序列{△Yrate}、{△Yrate}和{△YGtoll}，再次檢驗，結(jié)果見表2。

表1 單位根檢驗Tab.1 Unit Root Test

表2 一階差分單位根檢驗Tab.2 Unit Root Test after difference

由表2可知，3組序列一階差分后的檢驗τ統(tǒng)計量的值分別為-5.322 55、-6.766 90與-8.231 455，均小于顯著性水平為1%的臨界值，故可判定序列{△Ytoll}、{△Yrate}和{△YGtoll}皆無單位根，序列平穩(wěn)，已滿足了ARIMA建模條件。以序列{△Ytoll}為例對后續(xù)建模預測過程進行說明。

2.2.2 最優(yōu)預測模型的識別

首先，生成序列{△Ytoll}自相關系數(shù)(Autocorrelation Function， ACF)和偏自相關系數(shù)(Partial Autocorrelation Function， PACF)，如圖2。由圖2可知，序列{△Ytoll}的自相關系數(shù)從第2項開始均在虛線范圍內(nèi)波動，因此可視為1階截尾；偏自相關系數(shù)亦是從第2項開始在虛線范圍內(nèi)波動，也為1階截尾。由此ARIMA模型的階數(shù)可定為ARIMA(1，1，1)、ARIMA(1，1，0)以及ARIMA(0，1，1)。

圖2 “死亡人數(shù)”一階差分圖Fig.2 The first difference diagram of "the death toll"

表3為各模型AIC、SC檢驗值。根據(jù)最小信息量準則，選擇AIC與SC值都為最小的模型，即ARIMA(1，1，0)為“死亡人數(shù)”序列預測的最優(yōu)模型。按上述過程，“百萬噸死亡率”與“瓦斯事故死亡人數(shù)”的最優(yōu)模型為ARIMA(0，1，0)模型與ARIMA(8，1，1)模型。

表3 各模型AIC與SC檢驗值Tab.3 The AIC and SC value of each model

2.2.3 基于ARIMA模型的關鍵指標預測

2.3 基于混合模型的預測分析

根據(jù)式(10)得出殘差序列{etoll}、{erate}與{eGtoll}，在Rstudio軟件中調(diào)用forecastxgb包實現(xiàn)XGBoost建模。XGBoost算法對3組殘差序列的修正值與真實值的對比情況，根據(jù)式(11)得出混合模型預測值，如圖4?；旌夏Ｐ皖A測結(jié)果與原始序列對比情況，如圖5。

單一ARIMA模型、混合預測模型與單一Holt-Winters時間模型對3項指標預測結(jié)果的MAPE值，見表4。由表4可知，混合模型的預測結(jié)果相較于單一ARIMA模型、Holt-winters模型在預測精度上均有提升，混合模型“死亡人數(shù)”與“煤礦百萬噸死亡率”的MAPE值分別為14.121%與10.302%，混合模型對“瓦斯事故死亡人數(shù)”這一指標的預測精度提升最為明顯，從單一ARIMA模型20.496%，降至12.728%，預測精度提升了37.9%。

表4 3種模型的MAPE值Tab.4 The MAPE value of three models

進一步使用ARIMA-XGBoost混合模型對2021年3項指標的預測結(jié)果，見表5。

根據(jù)表5中的數(shù)據(jù)，2021年我國煤礦事故死亡人數(shù)約為151人與煤礦百萬噸死亡0.048 17，2項指標相較于2020年分別下降33.8%與16.9%，預測結(jié)果符合近幾年來煤礦安全生產(chǎn)形勢逐年向好的大趨勢。值得注意的是，瓦斯事故死亡人數(shù)的預測值約為59人，接近2020年的2倍，雖呈現(xiàn)較強的反撲態(tài)勢，但也符合2017年以來我國煤礦瓦斯事故死亡人數(shù)表現(xiàn)出的增減交替，波動下降的總趨勢，因此煤礦安全管理部門仍需加強對高瓦斯煤礦的監(jiān)管工作。2021年作為煤礦安全生產(chǎn)專項整治“三年行動”承上啟下的關鍵一年，進一步改善煤礦安全生產(chǎn)形勢，可借助煤礦智能化開采、淘汰不具備安全開采條件的礦井等舉措。

表5 2021年3項指標預測值Tab.5 Predictive values for three indicators in 2021

3 結(jié)論

(1)為充分挖掘煤礦安全態(tài)勢指標時間序列自身演變信息，在運用ARIMA模型提取時間序列非線性主部的基礎上，進一步利用XGBoost算法修正預測殘差，建立ARIMA-XGBoost煤礦安全態(tài)勢指標混合預測模型。使用該混合模型對3項煤礦安全態(tài)勢指標進行預測，預測結(jié)果的MAPE值均低于單一ARIMA、Holt-Winters模型，因此該混合模型更適用于對煤礦安全態(tài)勢指標的預測。

(2)根據(jù)ARIMA-XGBoost混合模型的預測結(jié)果，2021年我國煤礦事故死亡人數(shù)與煤礦百萬噸死亡率將繼續(xù)保持下降趨勢，進一步降低煤礦事故死亡人數(shù)必須依賴于智能采煤工作面等新技術的應用。對于瓦斯事故死亡人數(shù)這一指標的抬頭趨勢，煤礦安全管理部門繼續(xù)加強整治力度，減少瓦斯突出、瓦斯爆炸等事故的發(fā)生。