基于隨機森林回歸的圍巖應力插值方法

王蘇健，賈澎濤，金聲堯

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

0 引 言

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

圍巖應力監測數據是典型的時間序列數據，可以采用時間序列數據插值方法進行插補。常用的時間序列插值方法有傳統的數理統計方法和機器學習方法。傳統的數理統計方法有均值插值、中值插值、線性插值、最鄰近插值、拉格朗日插值、樣條插值等[8]?；跈C器學習的插值方法有基于聚類的插值[9-10]、基于回歸分析的插值[11]、基于神經網絡的插值[12-14]、基于集成學習的插值方法[15]等。

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

綜合來看，以上學者多采用的是改進的傳統插值方法，對基于機器學習的插值方法在圍巖應力數據領域的研究還較少。

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

1 基于隨機森林回歸的圍巖應力數據補全

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

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

1.1 圍巖應力回歸決策樹

隨機森林回歸預測算法的基礎是回歸決策樹CART算法。下面介紹圍巖應力回歸決策樹建模過程。

(1)

(2)

(3)

(4)

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

1.2 IRFR模型構建

CART決策樹能在一定程度上有效表示原始訓練樣本中的潛在統計關系，但往往較為粗糙，且不穩定。因此，應用集成學習的思想，在單棵CART樹的基礎上，構建基于隨機森林回歸方法的圍巖應力插值模型——IRFR模型，則可有效彌補單棵CART樹的不足。IRFR模型結構如圖1所示。

圖1 IRFR模型結構Fig.1 Structure diagram of IRFR model

從圖1可以看出，IRFR模型構建步驟如下

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

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

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

步驟4：根據數據的屬性特征，生成的B棵決策樹按照以下規則生成IRFR預測模型

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

(5)

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

2 實驗數據與實驗設計

2.1 實驗數據

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

圖2 數據缺失頻度情況分布Fig.2 Distribution of data missing frequency

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

1)連續缺失值1～6個，對應缺失值較少的情況。

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

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

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

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

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

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

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

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

2.2 數據預處理

對實驗數據集做歸一化處理，將數值規范化到[0，1]之間。設圍巖應力時間序列為R={x1,x2,…,xn}，標準化公式見式(6)。

(6)

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

圖3 實驗數據Fig.3 Experimental data

2.3 對比插值方法

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

2.4 IRFR參數確定

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

預測步長L設定為缺失數據個數；訓練窗口的長度n一般由經驗給出，這里設定為2倍的最大缺失值個數，即360。

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

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

(7)

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

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

(8)

式中ε為預先給定的精度值。當δ<ε(ε>0)時，B為所求的最優棵數。

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

2.5 實驗結果評價

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

3 結果分析

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

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

3.1 缺失值較少的情況

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

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

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

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

3.2 缺失值較多的情況

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

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

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

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

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

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

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

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

4 結 論

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

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

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