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

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

0 引言

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

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

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

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

1.1 ARIMA模型

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

式中：

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

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

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

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

ε—預(yù)測(cè)殘差；

μ—常數(shù)。

1.2 XGBoost算法

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

(4)

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

式中：

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

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

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

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

F—回歸樹(shù)的集合空間；

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

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

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

(5)

(6)

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

(7)

式中：

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

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

(8)

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

(9)

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

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

(10)

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

(11)

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

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

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

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

2.2 ARIMA模型預(yù)測(cè)

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

根據(jù)圖1中的折線變化趨勢(shì)，初步判斷3組原始序列均不符合零均值同方差的特性。3組原始序列單位根檢驗(yàn)結(jié)果，見(jiàn)表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)行差分處理。對(duì)3組原始序列一階差分后得到新序列{△Yrate}、{△Yrate}和{△YGtoll}，再次檢驗(yàn)，結(jié)果見(jiàn)表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ú)單位根，序列平穩(wěn)，已滿足了ARIMA建模條件。以序列{△Ytoll}為例對(duì)后續(xù)建模預(yù)測(cè)過(guò)程進(jìn)行說(shuō)明。

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

首先，生成序列{△Ytoll}自相關(guān)系數(shù)(Autocorrelation Function， ACF)和偏自相關(guān)系數(shù)(Partial Autocorrelation Function， PACF)，如圖2。由圖2可知，序列{△Ytoll}的自相關(guān)系數(shù)從第2項(xiàng)開(kāi)始均在虛線范圍內(nèi)波動(dòng)，因此可視為1階截尾；偏自相關(guān)系數(shù)亦是從第2項(xiàng)開(kāi)始在虛線范圍內(nèi)波動(dòng)，也為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ù)測(cè)的最優(yōu)模型。按上述過(guò)程，“百萬(wàn)噸死亡率”與“瓦斯事故死亡人數(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ù)測(cè)

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

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

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

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

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

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

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

3 結(jié)論

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

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