礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測

黃 悅，幸福堂，2，石癸鑫，吳孟龍

(1.武漢科技大學(xué) 資源與環(huán)境工程學(xué)院，湖北 武漢 430081；2.武漢科技大學(xué) 冶金礦產(chǎn)資源高效利用與造塊重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430081)

近年來，隨著采礦業(yè)的發(fā)展，礦業(yè)安全生產(chǎn)事故數(shù)量和傷亡人數(shù)居高不下，給社會(huì)發(fā)展帶來了負(fù)面影響。據(jù)近年來我國礦山事故統(tǒng)計(jì)顯示，2005年至2019年，全國共發(fā)生礦業(yè)安全生產(chǎn)事故3 826起[1]，死亡人數(shù)達(dá)到14 251人[1]，事故起數(shù)和死亡人數(shù)超過發(fā)達(dá)國家，安全問題嚴(yán)重制約了我國礦山企業(yè)的發(fā)展。因此，準(zhǔn)確預(yù)測礦業(yè)安全生產(chǎn)事故對制訂合理的礦業(yè)安全決策，以及提高礦業(yè)安全管理水平具有重要的意義。

礦業(yè)安全系統(tǒng)是一個(gè)復(fù)雜的循環(huán)系統(tǒng)[2]，礦業(yè)安全生產(chǎn)受自然條件、作業(yè)環(huán)境等諸多因素影響[3]。目前，國內(nèi)外眾多學(xué)者對礦業(yè)安全生產(chǎn)事故預(yù)測的方法展開了研究[4]。劉素兵等[5]將GM(1,1)與支持向量機(jī)融合，建立了灰色支持向量機(jī)組合預(yù)測模型；李懷俊等[6]通過參數(shù)尋優(yōu)改進(jìn)了ARMA模型，提高了ARMA模型的預(yù)測精度；Barman等[7]通過螢火蟲算法(Firefly Algorithm，F(xiàn)A)優(yōu)化SVM的參數(shù)，提高了SVM模型的預(yù)測精度；劉俊娥等[8]提出了一種模糊信息?；cSVM相結(jié)合的方法，對瓦斯涌出量進(jìn)行了預(yù)測；JANG Y J等[9]使用長短期記憶(LSTM)遞歸神經(jīng)網(wǎng)絡(luò)(RNN)的深度學(xué)習(xí)模型來預(yù)測業(yè)務(wù)失??；WU Menglong等[10]對經(jīng)驗(yàn)?zāi)B(tài)分解后的序列單獨(dú)分析,建立預(yù)測模型來預(yù)測礦山安全生產(chǎn)形勢。然而上述方法預(yù)測效果的魯棒性較差，并且只能得到礦業(yè)安全生產(chǎn)事故的點(diǎn)預(yù)測結(jié)果，目前關(guān)于礦業(yè)安全區(qū)間預(yù)測的研究還較少。區(qū)間預(yù)測是一類既含有確定性又含有不確定性的預(yù)測[11]。為此，考慮將模糊信息?；氲降V業(yè)安全生產(chǎn)事故預(yù)測中，得到礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測結(jié)果。區(qū)間預(yù)測能夠提供對應(yīng)某一期望概率的區(qū)間預(yù)測結(jié)果，獲取比確定性預(yù)測結(jié)果更多的信息，有利于決策者評估礦業(yè)安全生產(chǎn)事故的不確定性和風(fēng)險(xiǎn)因素，做出更為合理的決策[12]。因此，準(zhǔn)確預(yù)測礦業(yè)安全生產(chǎn)態(tài)勢的變化趨勢和變化空間，在實(shí)際應(yīng)用中具有重要的參考價(jià)值和指導(dǎo)意義。

筆者先將礦業(yè)安全生產(chǎn)態(tài)勢時(shí)間序列進(jìn)行小波變換，然后采用樣本熵將子序列重組，將重組的分量按照一定窗口進(jìn)行模糊信息?；?，得到能夠反映原始數(shù)據(jù)特點(diǎn)的信息粒；最后基于FCM對窗口化的 3個(gè)參數(shù)序列進(jìn)行預(yù)測，以實(shí)現(xiàn)礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測。

1 模型的建立

通過小波變換和樣本熵(SE)算法將礦業(yè)安全生產(chǎn)事故的時(shí)間序列分解并重組，然后采用模糊信息?；姆椒ǎ瑢⒌V業(yè)安全生產(chǎn)事故時(shí)間序列模糊?；玫絃ow、R、Up共3個(gè)?；瘏?shù)時(shí)間序列；最后依據(jù)模糊C均值聚類算法對窗口化的時(shí)間序列進(jìn)行預(yù)測，建立礦業(yè)安全生產(chǎn)事故時(shí)間序列預(yù)測模型。

1.1 時(shí)間序列的分解與重組

原始時(shí)間序列x(t)的小波變換可以表示為[13]：

(1)

式中：s為伸縮因子；t為時(shí)間；τ為平移因子；φ*為小波基函數(shù)φs,τ的復(fù)共軛。

小波基基于以下2個(gè)原則確定：

1)能—熵比準(zhǔn)則

將能量與熵結(jié)合起來形成一個(gè)新準(zhǔn)則，即能—熵比準(zhǔn)則[14]，見下式：

R(j)=Eenergy(j)/Eentropy(j)

(2)

式中：Eenergy(j)為原始時(shí)間序列在第j分解水平時(shí)小波系數(shù)序列包含的能量；Eentropy(j)為原始時(shí)間序列在第j分解水平時(shí)小波系數(shù)序列的Shannon熵。

由式(2)可以看出，能量越大，熵越小，則能—熵比越大，小波基[15]越優(yōu)。

2)相關(guān)性準(zhǔn)則

在處理和分析時(shí)間序列時(shí)，相關(guān)性分析[16]是一種在時(shí)域中對序列特性進(jìn)行描述的重要方法。

礦業(yè)安全生產(chǎn)事故時(shí)間序列的細(xì)節(jié)和近似部分可通過低通和高通濾波器獲得，通過對第(j-1)分解水平的近似系數(shù)與低通濾波器系數(shù)進(jìn)行卷積，可以獲得第j分解水平的近似系數(shù)：

(3)

式中：aj,k為近似系數(shù)；dj,k為細(xì)節(jié)系數(shù)；h、g分別為所選取的小波函數(shù)對應(yīng)的低通和高通濾波器的抽頭系數(shù)序列；m為維度；k=1,2,…,m。

經(jīng)過小波變換后的時(shí)間序列具有一定重復(fù)性和復(fù)雜性，故采用樣本熵算法對各分量時(shí)間序列進(jìn)行重組。

樣本熵的計(jì)算過程如下[17]：

1)將小波變換后的序列按序號形成n維的向量組Xm(1),…,Xm(N-m+1),其中Xm(i)={x(i),x(i+1),…,x(i-m+1)}，1≤i≤N-m+1。

(4)

3)定義B(m)(r)為：

(5)

4)將維數(shù)m增加至m+1維，對m+1維向量重復(fù)步驟1)～3)得到B(m+1)(r)。

5)該序列的樣本熵值為：

(6)

式中:m為重構(gòu)維數(shù),一般取m=2；r為相似容限，一般為原始數(shù)據(jù)標(biāo)準(zhǔn)差的0.10～0.25；n為數(shù)據(jù)長度。

1.2 重組序列的模糊信息?；?/h3>

根據(jù)一定的窗口w將原始時(shí)間序列分割成一些子序列。原始時(shí)間序列窗口以1個(gè)季度為1個(gè)單元，然后進(jìn)行模糊信息粒化，得到窗口化后的3個(gè)?；瘏?shù)時(shí)間序列[18]。模糊信息粒化后模糊粒子的隸屬函數(shù)為[19]：

(7)

式中：x為論域中的變量；a、m、b為參數(shù)，對應(yīng)窗口模糊化后得到的3個(gè)變量Low、R、Up；Low為窗口變化的最小值，R為窗口變化的平均值，Up為窗口變化的最大值。

1.3 基于FCM的窗口化序列模糊預(yù)測

基于FCM對窗口化的時(shí)間序列進(jìn)行預(yù)測，F(xiàn)CM的建模過程如下[20]：

1)確定聚類中心數(shù)。設(shè)原始礦業(yè)安全生產(chǎn)事故時(shí)間序列為Xt，t=1，2，…，N；待預(yù)測值為XN+1。首先對原始序列進(jìn)行聚類，通過下式確定原始序列的聚類中心數(shù)：

(8)

式中：c為聚類中心數(shù)；Xmax、Xmin為原始序列中的最大值和最小值；[]為取整運(yùn)算。

2)初始化聚類中心。根據(jù)式(8)計(jì)算得到聚類中心數(shù)，隨機(jī)選定初始聚類中心Ck。Ck是一個(gè)長度為c的向量，k=1，2…，c。

3)計(jì)算目標(biāo)化函數(shù)。FCM的目標(biāo)化函數(shù)為：

(9)

式中：Ck為模糊組i的聚類中心；uij為屬于第i個(gè)聚類的隸屬度，介于0～1；dij為第i個(gè)聚類中心與第j個(gè)數(shù)據(jù)點(diǎn)間的歐幾里得距離。

2 模型應(yīng)用

以2005—2019年礦業(yè)單位從業(yè)人員事故的萬人死亡率為研究對象,統(tǒng)計(jì)得到全國礦業(yè)事故萬人死亡率時(shí)間序列,如圖1所示。

圖1 2005—2019年全國礦業(yè)事故萬人死亡率時(shí)間序列

2.1 小波變換

采用能—熵比準(zhǔn)則及相關(guān)性準(zhǔn)則選取最優(yōu)小波基，并依據(jù)相關(guān)性準(zhǔn)則確定臨界分解水平。

筆者采用較為常用的dbN(1～10)、coifN(1～5)、symN(2～10)3種小波系選取最優(yōu)小波基。

2.1.1 能—熵比準(zhǔn)則的計(jì)算結(jié)果

基于能—熵比準(zhǔn)則的計(jì)算結(jié)果如表1所示。

表1 基于能—熵比準(zhǔn)則R(j)的計(jì)算結(jié)果

將表1中數(shù)據(jù)依據(jù)時(shí)間序列比值高低排序，結(jié)果如表2所示。

表2 基于能—熵比準(zhǔn)則R(j)對小波基的排序

2.1.2 相關(guān)性準(zhǔn)則的計(jì)算結(jié)果

基于相關(guān)性準(zhǔn)則的計(jì)算結(jié)果如表3所示。

表3 基于相關(guān)性準(zhǔn)則的計(jì)算結(jié)果

將表3中數(shù)據(jù)依據(jù)相關(guān)性系數(shù)大小進(jìn)行排序，排序結(jié)果如表4所示。

表4 基于相關(guān)性準(zhǔn)則對小波基的排序結(jié)果

2.1.3 小波基優(yōu)選結(jié)果

計(jì)算表2、表4中小波基的得分，根據(jù)其對應(yīng)的分值進(jìn)行排序，排序結(jié)果如表5所示。

表5 小波基排序結(jié)果

由表5可知，小波基coif5得分最高，coif5為最優(yōu)小波基。

2.1.4 臨界分解水平的確定

基于臨界分解水平判別準(zhǔn)則，對曲線的臨界分解水平進(jìn)行分析判斷,得到各分解水平趨勢序列與原始礦業(yè)安全生產(chǎn)事故時(shí)間序列的相關(guān)性系數(shù),結(jié)果如圖2所示。

圖2 各分解水平的相關(guān)性系數(shù)

由圖2可知，礦業(yè)安全生產(chǎn)事故時(shí)間序列曲線分解水平大于5時(shí)，相關(guān)性系數(shù)發(fā)生顯著變化，故確定臨界分解水平為5。

2.1.5 時(shí)間序列的分解

在確定最優(yōu)小波基和臨界分解水平后，利用MATLAB小波工具箱對其進(jìn)行小波變換。礦業(yè)安全生產(chǎn)事故時(shí)間序列5分解水平小波變換后的結(jié)果如圖3所示(圖中S代表原始序列，ai代表近似部分，di代表細(xì)節(jié)部分)。

(a)細(xì)節(jié)部分 (b)近似部分

由圖3可知，隨著分解水平的增加，細(xì)節(jié)部分逐漸被提取出來。經(jīng)過5分解水平小波變換后，可以滿足礦業(yè)安全生產(chǎn)事故時(shí)間序列分解的要求。

2.2 分量重組

基于礦業(yè)安全生產(chǎn)事故時(shí)間序列進(jìn)行小波變換后的子序列，計(jì)算給定序列與子序列的樣本熵值，各分量的樣本熵如圖4所示。

圖4 各分量的樣本熵

根據(jù)子序列樣本熵的分布情況，取λ=0.2。d1為隨機(jī)分量；d2～d4重組為細(xì)節(jié)分量s；d5、a5重組為趨勢分量r。重組后各分量結(jié)果如圖5所示。

圖5 重組后各分量

2.3 重組分量的信息?；?/h3>

由圖5可知，樣本熵重組后的序列具有隨機(jī)性、波動(dòng)性等特點(diǎn)。為此，先將原始數(shù)據(jù)按照一定窗口進(jìn)行模糊信息?；?，將礦業(yè)安全生產(chǎn)態(tài)勢時(shí)間序列的萬人死亡率模糊信息粒化為Low、R、Up共3個(gè)參數(shù)，3個(gè)子序列的模糊信息?；Y(jié)果如圖6所示。該例中Low具體為窗口化子序列每月變化的最小值，R具體為礦業(yè)安全生產(chǎn)態(tài)勢萬人死亡率每月變化的平均值，Up具體為礦業(yè)安全生產(chǎn)態(tài)勢萬人死亡率每月變化的最大值。

(a)d2序列粒化

(b)d3序列?；?/p>

(c)d4序列?；?/p>

2.4 基于FCM的窗口化子序列模糊預(yù)測

對于窗口化的時(shí)間序列，按照窗口模糊信息?；Y(jié)果將其分成6個(gè)數(shù)據(jù)集，根據(jù)每年數(shù)據(jù)情況確定討論區(qū)間的上下界計(jì)算模型的論域。基于樣本數(shù)據(jù)求得的隸屬度值見表6，樣本模糊化后的結(jié)果如表7 所示。

表6 樣本隸屬度

表7 樣本模糊關(guān)系

由表6和表7并結(jié)合式(9)計(jì)算，求得測試集預(yù)測結(jié)果，如表8所示。

表8 測試集預(yù)測結(jié)果

合并窗口化后，礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測結(jié)果如圖7所示。

圖7 窗口化區(qū)間預(yù)測結(jié)果

3 模型的驗(yàn)證

利用訓(xùn)練集150組樣本建立的模型，對2018年7月至2019年10月測試集的21組數(shù)據(jù)進(jìn)行預(yù)測，并計(jì)算區(qū)間上下邊界值與中值的預(yù)測值，結(jié)果如表9 所示。

表9 測試集預(yù)測誤差對比

由表9可知：R值的平均相對誤差為17.737 1%，Up值的平均相對誤差為8.771 2%。Up值的預(yù)測精度明顯高于R值的預(yù)測精度，表明模糊信息?；蟮纳舷逻吔缰的軌蜉^為準(zhǔn)確地描述礦業(yè)安全生產(chǎn)事故的變化范圍。

將礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測模型與小波-SVM(支持向量機(jī))、小波-ARMA(自回歸滑動(dòng)平均模型)和小波神經(jīng)網(wǎng)絡(luò)的預(yù)測效果進(jìn)行對比，各模型的預(yù)測誤差分析如表10所示。其中，MAE為預(yù)測值和實(shí)際值之間的平均絕對誤差；MRE為絕對誤差與真實(shí)值的比率，即平均相對誤差；RMSE為均方根誤差，其對異常值更為敏感。

表10 各模型預(yù)測誤差對比分析結(jié)果

由表10可知，本文預(yù)測模型的MAE、MRE、RMSE分別為0.445 7%、13.254 1%、0.372 3%，各項(xiàng)指標(biāo)均低于對比模型，表明本文模型預(yù)測效果較為穩(wěn)定，精度較高。

4 結(jié)論

1)提出一種礦業(yè)安全生產(chǎn)事故時(shí)間序列的模糊區(qū)間預(yù)測模型，能夠?qū)⒌V業(yè)系統(tǒng)中的隨機(jī)因素等考慮進(jìn)去，相較于定量值的點(diǎn)預(yù)測，區(qū)間預(yù)測的結(jié)果更為可靠，更具有實(shí)用價(jià)值。

2)R、Up值時(shí)間序列預(yù)測的平均相對誤差分別為17.737 1%、8.771 2%。邊界值的預(yù)測精度較高，表明模型的區(qū)間范圍較為可靠，也進(jìn)一步地說明了模糊信息粒化后的邊界值結(jié)果的合理性。

3)相對于小波-支持向量機(jī)、小波-自回歸滑動(dòng)平均模型和小波神經(jīng)網(wǎng)絡(luò)的預(yù)測結(jié)果，本文的預(yù)測模型平均絕對誤差(MAE)、平均相對誤差(MRE)、均方根誤差(RMSE)3個(gè)指標(biāo)均為最小，表明訓(xùn)練得到的模型具有較高的可靠性，對礦業(yè)安全決策具有一定的指導(dǎo)意義。