王蘇健，賈澎濤，金聲堯

(1.陜西煤業(yè)化工技術(shù)研究院有限責(zé)任公司，陜西 西安 710100；2.西安科技大學(xué) 計算機科學(xué)與技術(shù)學(xué)院，陜西 西安 710054)

0 引 言

隨著我國礦井開采深度的增加，因圍巖應(yīng)力增高而造成的事故時有發(fā)生，嚴(yán)重威脅著煤礦生產(chǎn)安全[1]。為了有效控制圍巖應(yīng)力，需要采取必要的監(jiān)測預(yù)警手段，實現(xiàn)對礦井巷道和采掘工作面的實時監(jiān)控預(yù)警[2-4]。在圍巖應(yīng)力安全監(jiān)測中，采集數(shù)據(jù)的準(zhǔn)確與否對礦山安全預(yù)警結(jié)果有著至關(guān)重要的影響[5]。但是在監(jiān)測數(shù)據(jù)采集過程中，由于受井下惡劣的現(xiàn)場環(huán)境、設(shè)備和線路故障、人為因素等條件限制，采集到的圍巖應(yīng)力數(shù)據(jù)經(jīng)常會存在缺失值的情況，嚴(yán)重影響了后續(xù)的數(shù)據(jù)分析、安全預(yù)警等工作。因此，對圍巖應(yīng)力監(jiān)測數(shù)據(jù)的缺失值進行有效插補是亟待解決的問題[6-7]。

圍巖應(yīng)力監(jiān)測數(shù)據(jù)是典型的時間序列數(shù)據(jù)，可以采用時間序列數(shù)據(jù)插值方法進行插補。常用的時間序列插值方法有傳統(tǒng)的數(shù)理統(tǒng)計方法和機器學(xué)習(xí)方法。傳統(tǒng)的數(shù)理統(tǒng)計方法有均值插值、中值插值、線性插值、最鄰近插值、拉格朗日插值、樣條插值等[8]。基于機器學(xué)習(xí)的插值方法有基于聚類的插值[9-10]、基于回歸分析的插值[11]、基于神經(jīng)網(wǎng)絡(luò)的插值[12-14]、基于集成學(xué)習(xí)的插值方法[15]等。

在圍巖應(yīng)力數(shù)據(jù)插值領(lǐng)域，陸振裕等采用線性插值法，對于鉆屑法獲得圍巖壓力系列離散點數(shù)據(jù)進行插值，為后期沖擊地壓危險預(yù)測提供了有效的數(shù)據(jù)[16]。李凱拓同樣采用線性插值法，對王莊煤礦三維地應(yīng)力分量做了平滑處理，最終繪制出了各應(yīng)力量分布云圖[17]。宋道柱等采用等值線插值對單體液壓支柱壓力的誤差數(shù)據(jù)進行修正，采用雙線性內(nèi)插法對缺失值數(shù)據(jù)進行插值，為支撐受力變化趨勢分析提供有效數(shù)據(jù)[18]。尹時雨采用3次樣條插值以及高階多元非線性回歸方程，對工作面各支架測點工作阻力數(shù)據(jù)進行插值擬合，為研究礦山壓力顯現(xiàn)規(guī)律提供理論依據(jù)[19]。徐恩虎等提出了一種連續(xù)插值模型，可以得到等步長與非等步長數(shù)例的任意內(nèi)插值與外插值[20]。陳輝等提出了一種基于粒子群的3次樣條插值算法，取得了比最鄰近插值、拉格朗日插值、線性插值方法更好的插值效果[21]。馮俊軍等采用克里金插值法得到的應(yīng)力云圖能夠反映試驗工作面巷幫煤體應(yīng)力場分布規(guī)律，為巷道超前支護、頂板穩(wěn)定性控制、圍巖穩(wěn)定性分析、沖擊地壓預(yù)防等提供理論依據(jù)[22]。

綜合來看，以上學(xué)者多采用的是改進的傳統(tǒng)插值方法，對基于機器學(xué)習(xí)的插值方法在圍巖應(yīng)力數(shù)據(jù)領(lǐng)域的研究還較少。

由于圍巖應(yīng)力數(shù)據(jù)缺失問題可以看作預(yù)測問題，用缺失值之前的數(shù)據(jù)預(yù)測缺失位置的數(shù)據(jù)，從而實現(xiàn)插值。因此，文中嘗試采用機器學(xué)習(xí)中有監(jiān)督集成學(xué)習(xí)方法——隨機森林回歸方法(random forest regression，RFR)，進行缺失數(shù)據(jù)的補全處理。首先，在集成學(xué)習(xí)理論的基礎(chǔ)上，提出基于隨機森林回歸的圍巖應(yīng)力插值方法(interpolation based on random forest regression，IRFR)。然后，以無缺失的圍巖應(yīng)力時間序列數(shù)據(jù)為樣本集，構(gòu)建不同缺失情況的數(shù)據(jù)集，作為實驗用數(shù)據(jù)。最后，在不同缺失值情況下，研究不同插值方法對不同缺失情況的圍巖應(yīng)力數(shù)據(jù)的影響，并就其插值效果進行比較分析。實驗結(jié)果驗證了IRFR的正確性與有效性。

1 基于隨機森林回歸的圍巖應(yīng)力數(shù)據(jù)補全

隨機森林算法是由Brieman于2001年提出的一種集成學(xué)習(xí)算法[23]，用于解決高維非線性數(shù)據(jù)的分類預(yù)測、回歸預(yù)測與特征選擇問題。隨機森林回歸預(yù)測算法是bagging算法的改進算法，它用K個分類回歸決策樹(classification and regression tree，CART)作為基學(xué)習(xí)器[24]，以K個基學(xué)習(xí)器預(yù)測值的平均值作為最終結(jié)果。

基于隨機森林回歸方法的圍巖應(yīng)力插值方法(IRFR)的思路是輸入圍巖應(yīng)力監(jiān)測數(shù)據(jù)，用當(dāng)前缺失值之前的指定窗口長度的歷史數(shù)據(jù)訓(xùn)練隨機森林回歸模型，并輸出預(yù)測值，該預(yù)測值就為插值的結(jié)果。如果是連續(xù)缺失的情況，采用遞推方法進行預(yù)測補全。

1.1 圍巖應(yīng)力回歸決策樹

隨機森林回歸預(yù)測算法的基礎(chǔ)是回歸決策樹CART算法。下面介紹圍巖應(yīng)力回歸決策樹建模過程。

(1)

(2)

(3)

(4)

按照上述分割方法，分別將Rl和Rr作為父節(jié)點，遞歸進行分割，直至當(dāng)前父節(jié)點中樣本的y值方差小于給定方差閾值。條件滿足時，停止遞歸并將當(dāng)前父節(jié)點設(shè)置為葉子節(jié)點。至此，單棵圍巖應(yīng)力CART樹就建立起來了。

1.2 IRFR模型構(gòu)建

CART決策樹能在一定程度上有效表示原始訓(xùn)練樣本中的潛在統(tǒng)計關(guān)系，但往往較為粗糙，且不穩(wěn)定。因此，應(yīng)用集成學(xué)習(xí)的思想，在單棵CART樹的基礎(chǔ)上，構(gòu)建基于隨機森林回歸方法的圍巖應(yīng)力插值模型——IRFR模型，則可有效彌補單棵CART樹的不足。IRFR模型結(jié)構(gòu)如圖1所示。

圖1 IRFR模型結(jié)構(gòu)Fig.1 Structure diagram of IRFR model

從圖1可以看出，IRFR模型構(gòu)建步驟如下

步驟1：采用Bootstrap抽樣技術(shù)，從原始樣本集R中有放回的抽取B個樣本集，構(gòu)建B個CART樹進行學(xué)習(xí)訓(xùn)練，剩余未被抽取的樣本作為袋外(out of bag，OOB)數(shù)據(jù)，形成模型的測試樣本數(shù)據(jù)。

步驟2：設(shè)原始數(shù)據(jù)集變量個數(shù)為p，在每棵決策樹模型的內(nèi)部節(jié)點隨機抽取k(k≤p)個變量作為備選分枝變量，按照單棵決策樹構(gòu)建過程尋找最佳分枝。

步驟3：每棵決策樹自頂向下遞歸分枝，直至當(dāng)前父節(jié)點中樣本的y值方差小于給定方差閾值。條件滿足時，停止遞歸并將當(dāng)前父節(jié)點設(shè)置為葉子節(jié)點。

步驟4：根據(jù)數(shù)據(jù)的屬性特征，生成的B棵決策樹按照以下規(guī)則生成IRFR預(yù)測模型

yIRFR=ave(y(X,Tb)),b=1,2,…,B

(5)

式中Tb為第b棵回歸樹；y(X,Tb)為第b棵回歸樹的預(yù)測值。yIRFR的預(yù)測值是B棵樹預(yù)測值的平均值。

2 實驗數(shù)據(jù)與實驗設(shè)計

2.1 實驗數(shù)據(jù)

在某煤礦綜采工作面進行圍巖應(yīng)力數(shù)據(jù)采樣，從2019年5月7日10：00開始至2019年11月6日7：30分結(jié)束。按照采樣間隔10 min，應(yīng)采樣26 338個數(shù)據(jù)，實際采樣14 306個數(shù)據(jù)，缺失12 032個數(shù)據(jù)，缺失情況較為嚴(yán)重。連續(xù)缺失數(shù)據(jù)的頻度情況如圖2所示。

圖2 數(shù)據(jù)缺失頻度情況分布Fig.2 Distribution of data missing frequency

從圖2可以看出，數(shù)據(jù)連續(xù)缺失1～30個的情況較多，連續(xù)缺失30個值以上的情況較少，連續(xù)缺失180個以上的情況只有2次。因此，為了方便實驗，按照數(shù)據(jù)缺失的分布情況，把原始數(shù)據(jù)的缺失情況歸納為8種情況，這8種情況涵蓋了大多數(shù)的數(shù)據(jù)缺失情況。

1)連續(xù)缺失值1～6個，對應(yīng)缺失值較少的情況。

2)連續(xù)缺失值12個(2個小時)；

3)連續(xù)缺失值30個(5個小時)。

4)連續(xù)缺失值60個(10個小時)。

5)連續(xù)缺失值90個(15個小時)。

6)連續(xù)缺失值120個(20個小時)。

7)連續(xù)缺失值150個(25個小時)。

8)連續(xù)缺失值180個(30個小時)。

為了驗證模型效果，選取無缺失值的一段序列，即2019年7月8日00：00至2019年7月19日23：50的1 728個圍巖應(yīng)力傳感數(shù)據(jù)，作為實驗數(shù)據(jù)。實驗數(shù)據(jù)均值為24.317 7 MPa，標(biāo)準(zhǔn)差為6.508 4 MPa，最小值0.4 MPa，最大值45.8 MPa。按照以上8種缺失值情況，人為設(shè)置其缺失值，便于后續(xù)實驗比較。

2.2 數(shù)據(jù)預(yù)處理

對實驗數(shù)據(jù)集做歸一化處理，將數(shù)值規(guī)范化到[0，1]之間。設(shè)圍巖應(yīng)力時間序列為R={x1,x2,…,xn}，標(biāo)準(zhǔn)化公式見式(6)。

(6)

式中i=1，2，…，n。歸一化后的數(shù)據(jù)如圖3所示。

圖3 實驗數(shù)據(jù)Fig.3 Experimental data

2.3 對比插值方法

選擇均值插值、中值插值、線性插值、最鄰近插值(簡稱最鄰近)、Zero階梯插值(簡稱Zero)、3次B樣條插值(簡稱Cube)、拉格朗日3次多項式插值(簡稱Lagrange)7種插值方法作為實驗對比插值方法。這些插值方法均為較常用的插值方法，這里對其原理不再贅述。

2.4 IRFR參數(shù)確定

針對不同的數(shù)據(jù)缺失情況，IRFR模型需確定預(yù)測步長L、訓(xùn)練窗口長度n和決策樹棵樹B。

預(yù)測步長L設(shè)定為缺失數(shù)據(jù)個數(shù)；訓(xùn)練窗口的長度n一般由經(jīng)驗給出，這里設(shè)定為2倍的最大缺失值個數(shù)，即360。

IRFR模型中，決策樹的棵數(shù)對預(yù)測結(jié)果的準(zhǔn)確率和性能有較大的影響。因此針對某種缺失情況，通過比較不同決策樹棵樹下訓(xùn)練集的均方誤差(mean square error，MSE)，從而確定決策樹的棵數(shù)。

設(shè)圍巖應(yīng)力時間序列為R={x1,x2,…,xn,xn+1,…,xn+s}，某一時刻圍巖應(yīng)力的真實值為xi，擬合的插值為yi，預(yù)測步長為s，則預(yù)測值和真實值的MSE公式如式(7)所示。

(7)

設(shè)E={e1,e2,…,eB}為決策樹棵樹為1，2，…，B時的MSE誤差，相鄰誤差的差值δ定義為

δ=ei+1-ei(i=1,2,…,B-1)

(8)

式中ε為預(yù)先給定的精度值。當(dāng)δ<ε(ε>0)時，B為所求的最優(yōu)棵數(shù)。

以缺失值為120個的情況為例，不同棵數(shù)決策樹情況下的MSE如圖4所示。從圖4可以看出，隨著決策樹數(shù)量的增多，誤差呈下降趨勢。在棵樹為500時，相鄰誤差差值小于指定的ε=0.000 1，因此選擇決策樹棵樹為500。

2.5 實驗結(jié)果評價

采用均方誤差(MSE)作為評價指標(biāo)，對不同插值方法的結(jié)果進行比較。

3 結(jié)果分析

采用IRFR方法和2.3小節(jié)中的對比插值方法，在8種缺失值情況(2.1小節(jié)羅列)的數(shù)據(jù)集上進行試驗，將擬合值和真實值進行對比。

圖4 不同決策樹數(shù)量情況下的MSE誤差Fig.4 MSE errors of different number decision trees

3.1 缺失值較少的情況

針對缺失值較少的情況(1～6個)，各插值方法的誤差比較見表1(誤差值取小數(shù)點后4位)。限于篇幅，僅列出連續(xù)缺失6個值的插值效果圖，如圖5所示(均值插值和中值插值效果較差，未顯示)。

表1 缺失值較少情況下不同插值方法的誤差比較Table 1 Errors comparison of different interpolation methods with fewer missing values

圖5 連續(xù)缺失6個值情況下的插值效果對比Fig.5 Comparison of interpolation effects in the case of continuous missing 6 values

從表1可以看出，均值插值和中值插值效果較差；拉格朗日插值在連續(xù)缺失值1個和2個的情況下取得了最好的插值效果；線性插值在連續(xù)缺失值3～4個情況下取得了最好的插值效果；IRFR在連續(xù)缺失值5個情況下取得了最好的插值效果；從表1和圖5可以看出，在連續(xù)缺失6個值情況下，Cube方法取得了最好的插值效果，IRFR次之。因此，在連續(xù)缺失值較少的情況下，線性插值、最鄰近插值、Zero階梯插值、3次B樣條插值、拉格朗日插值和IRFR插值效果相當(dāng)。

3.2 缺失值較多的情況

針對缺失值較多的情況(連續(xù)缺失值個數(shù)為12、30、60、90、120、150和180)，各插值方法的插值效果見表2(誤差值取小數(shù)點后4位)。限于篇幅，此處只展示連續(xù)缺失12個值、30個值和180個值的插值效果圖，如圖6、圖7、圖8所示(圖7、圖8中，因為Lagrange插值效果太差，未顯示)。

從表1和表2可以看出，均值插值、中值插值在缺失值較少和較多的情況下，誤差變化情況不大，但是相較于除Lagrange的其他插值方法，效果一般；線性插值在缺失值少的情況下，取得了較好的插值效果，但隨著缺失值的增加，插值效果逐漸變差；Zero階梯插值隨著缺失值的增加，插值效果逐漸變差；Lagrange插值隨著缺失值數(shù)量增多，均方誤差增加較大，表現(xiàn)最差。

從表2和圖6至圖8可以看出，缺失值較多情況下，所有插值方法中IRFR方法取得了最好的插值效果，且隨著缺失值的增加，誤差也沒有明顯增大。

表2 缺失值較多情況下不同插值方法的誤差比較Table 2 Errors comparison of different interpolation methods with more missing values

圖6 連續(xù)缺失12個值情況下的插值效果對比Fig.6 Comparison of interpolation effects in the case of continuous missing 12 values

圖7 連續(xù)缺失30個值情況下的插值效果對比Fig.7 Comparison of interpolation effects in the case of continuous missing 30 values

圖8 連續(xù)缺失180個值情況下的插值效果對比Fig.8 Comparison of interpolation effects in the case of continuous missing 180 values

實驗結(jié)果表明，針對圍巖應(yīng)力數(shù)據(jù)，相較于對比的均值插值、中值插值、線性插值、最鄰近插值、Zero階梯插值、3次B樣條插值、拉格朗日3次多項式插值插值方法，IRFR方法為最佳插值方法。

4 結(jié) 論

1)提出了一種基于隨機森林回歸預(yù)測方法的圍巖應(yīng)力插值方法，該方法利用歷史數(shù)據(jù)訓(xùn)練隨機森林回歸模型，通過對缺失值進行預(yù)測，實現(xiàn)了針對圍巖應(yīng)力時間序列缺失值的插補。

2)對不同缺失值情況下，不同插值方法在圍巖應(yīng)力數(shù)據(jù)集上的插值效果進行了比較分析。

3)實驗結(jié)果表明，均值差值、中值插值方法效果較差；拉格朗日插值方法僅適用于缺失值較少的情況，在缺失值較多的情況下，均方誤差有隨著缺失值數(shù)量增多而增大的趨勢，效果最差；線性插值、最鄰近插值、Zero階梯插值、3次B樣條插值同樣適用于缺失值較少的情況，在缺失值較多的情況下效果較差；IRFR方法在不同缺失值情況下，均取得了較好的插值效果，且隨著缺失值數(shù)量的增加，這種優(yōu)勢尤為明顯。因此，IRFR方法適用于圍巖應(yīng)力插值。