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

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

0 引言

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

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

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

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

1.1 ARIMA模型

ARIMA模型是由Box和Jenkins[13]提出的一種時(shí)間序列預(yù)測模型，模型包含參數(shù)p、d、q，其中p為自回歸階數(shù)，d為原始序列平穩(wěn)化所需的差分階數(shù)，q為移動平均階數(shù)，亦可寫作ARIMA(p，d，q)模型。該模型可以通過歷史數(shù)據(jù)對未來進(jìn)行預(yù)測[14]。其建模的基本步驟是根據(jù)單位根檢驗(yàn)的結(jié)果，確定原始序列平穩(wěn)化所需的差分階數(shù)d，通過分析自相關(guān)與偏自相關(guān)函數(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é)合，再進(jìn)行d階差分得到ARIMA(p，d，q)模型，其數(shù)學(xué)描述可簡記為：

式中：

Yt—時(shí)間序列中第t項(xiàng)的觀測值，t為當(dāng)前項(xiàng)在時(shí)間序列中對應(yīng)的序號；

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

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

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

ε—預(yù)測殘差；

μ—常數(shù)。

1.2 XGBoost算法

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

(4)

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

式中：

xi—第i個(gè)數(shù)據(jù)點(diǎn)的特征向量；

yi—第i個(gè)數(shù)據(jù)點(diǎn)的真實(shí)值；

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

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

F—回歸樹的集合空間；

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

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

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

(5)

(6)

此時(shí)目標(biāo)函數(shù)可描述為：

(7)

式中：

C—為常數(shù)項(xiàng)。

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

(8)

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

(9)

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

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

(10)

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

(11)

2 實(shí)證計(jì)算與分析

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

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

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

2.2 ARIMA模型預(yù)測

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

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

表1 單位根檢驗(yàn)Tab.1 Unit Root Test

表2 一階差分單位根檢驗(yàn)Tab.2 Unit Root Test after difference

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

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

首先，生成序列{△Ytoll}自相關(guān)系數(shù)(Autocorrelation Function， ACF)和偏自相關(guān)系數(shù)(Partial Autocorrelation Function， PACF)，如圖2。由圖2可知，序列{△Ytoll}的自相關(guān)系數(shù)從第2項(xiàng)開始均在虛線范圍內(nèi)波動，因此可視為1階截尾；偏自相關(guān)系數(shù)亦是從第2項(xiàng)開始在虛線范圍內(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檢驗(yàn)值。根據(jù)最小信息量準(zhǔn)則，選擇AIC與SC值都為最小的模型，即ARIMA(1，1，0)為“死亡人數(shù)”序列預(yù)測的最優(yōu)模型。按上述過程，“百萬噸死亡率”與“瓦斯事故死亡人數(shù)”的最優(yōu)模型為ARIMA(0，1，0)模型與ARIMA(8，1，1)模型。

表3 各模型AIC與SC檢驗(yàn)值Tab.3 The AIC and SC value of each model

2.2.3 基于ARIMA模型的關(guān)鍵指標(biāo)預(yù)測

2.3 基于混合模型的預(yù)測分析

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

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

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

進(jìn)一步使用ARIMA-XGBoost混合模型對2021年3項(xiàng)指標(biāo)的預(yù)測結(jié)果，見表5。

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

表5 2021年3項(xiàng)指標(biāo)預(yù)測值Tab.5 Predictive values for three indicators in 2021

3 結(jié)論

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

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