基于主成分回歸分析法的導(dǎo)水裂縫帶高度預(yù)測

曹始友，成文舉，張歷峰，徐德寶，王 松，陳大林，王 鵬，尹會永

(1.棗莊礦業(yè)(集團)有限責(zé)任公司，山東 棗莊 277100；2.山東科技大學(xué)地球科學(xué)與工程學(xué)院山東省沉積成礦作用與沉積礦產(chǎn)重點實驗室，山東 青島 266590；3.棗莊礦業(yè)(集團)有限責(zé)任公司濱湖煤礦，山東 棗莊 277599；4.棗莊礦業(yè)集團新安煤業(yè)有限公司，山東 濟寧 277607；5.棗莊礦業(yè)(集團)付村煤業(yè)有限公司，山東 濟寧 277605；6.棗莊礦業(yè)(集團)有限責(zé)任公司蔣莊煤礦，山東 棗莊 277519；7.山東省三河口礦業(yè)有限責(zé)任公司，山東 濟寧 277600)

煤層開采后，圍巖的原始應(yīng)力平衡被打破，頂板上覆圍巖發(fā)生移動變形，極易形成裂縫，進而發(fā)育成為地表水及地下水運移的通道[1-3]。依據(jù)現(xiàn)有的“上三帶”理論[4-5]，冒落帶及裂隙帶共同組成的導(dǎo)水裂縫帶成為煤層頂板突水危險性的決定性因素，即原始狀態(tài)隔水層的隔水性是否在開采過程中遭到采動的破壞。影響導(dǎo)水裂縫帶發(fā)育高度的因素有很多，大部分可進行定量描述，但有些只能進行定性說明[6-8]。目前,研究認(rèn)為其主要影響因素有上覆巖層性質(zhì)、煤層厚度、開采深度、工作面傾向長度、頂板厚度及煤層傾角；另外，開采時間、重復(fù)采動等也會對導(dǎo)水裂縫帶高度產(chǎn)生影響。如何高精度地預(yù)測頂板導(dǎo)水裂縫帶發(fā)育高度，一直是許多科技工作者和現(xiàn)場工程技術(shù)人員工作研究的重點問題之一。

目前對于導(dǎo)水裂縫帶高度預(yù)測的方法仍以《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)范》[9]所提供的經(jīng)驗公式法為主，但隨著開采方式的變化及開采難度的增大，原有的經(jīng)驗方法已不能完全滿足生產(chǎn)的需求。目前國內(nèi)的學(xué)者及科技工作者對于導(dǎo)水裂縫帶高度的預(yù)測方法也做了一定的探索和研究，提出了一些具有實用價值的新方法，如：量綱分析法、PSO-RBF神經(jīng)網(wǎng)絡(luò)模型、偏最小二乘回歸法、基于最小二乘法的PLS-BP神經(jīng)網(wǎng)絡(luò)模型、支持向量機模型、模糊聚類支持向量機模型等[10-15]，另外還有SAS數(shù)學(xué)軟件[8]及多方法綜合預(yù)測[16]等。相較于傳統(tǒng)的預(yù)測方法，新的模擬模型及綜合方法在特定條件下具有更高的準(zhǔn)確性和針對性，但模型方法較為復(fù)雜，參數(shù)選取主觀性較大；多方法綜合確定需要較大的工作量，仍具有一定的不足。主成分分析，也稱主分量分析，是通過降維的方式將多指標(biāo)通過一定的重新組合轉(zhuǎn)化為少數(shù)幾個綜合指標(biāo)，即所謂主成分，用以解釋資料的綜合性指標(biāo)[17-18]。本文基于蔣莊煤礦及其周邊多個煤礦的礦井資料，選取了影響導(dǎo)水裂縫帶發(fā)育的主要影響因子，利用主成分分析法及MATLAB數(shù)據(jù)處理軟件，結(jié)合回歸模型檢驗[19]，建立蔣莊煤礦導(dǎo)水裂縫帶發(fā)育高度的回歸模型，對煤礦導(dǎo)水裂縫帶高度的預(yù)測方法的研究和應(yīng)用具有一定的借鑒意義。

1 方法和原理

主成分分析是在給定的一組相關(guān)變量的基礎(chǔ)上，通過一定的線性關(guān)系，將其轉(zhuǎn)換成另一組按照遞減方差的順序排列的不相關(guān)變量，確定各變量的主成分，在此基礎(chǔ)上建立起與實際相關(guān)的回歸模型。主成分分析法在水質(zhì)評價方面具有較廣的應(yīng)用，參照類比的思想，通過分析影響導(dǎo)水裂縫帶高度的主要因素，將其運用到導(dǎo)水裂縫帶高度的預(yù)測中，通過確定幾個最主要的影響因素，來替代原有非必要指標(biāo)，從而達到減少任務(wù)量的目的。其主要的計算步驟如下所述。

1) 構(gòu)建數(shù)據(jù)矩陣并標(biāo)準(zhǔn)化處理。設(shè)有n個礦區(qū)樣品，每個礦區(qū)樣品有p個主要影響因素，將原始數(shù)據(jù)寫成矩陣，見式(1)。

[X1,X2,…,Xp-1,Xp]

(1)

2) 標(biāo)準(zhǔn)化處理。由于不同指標(biāo)的量綱不同，在計算方差會產(chǎn)生很大的誤差，對判斷主成分具有較大的影響，因此需要將選用指標(biāo)進行無量綱化處理，即標(biāo)準(zhǔn)化，其公式見式(2)。

(2)

3) 計算相關(guān)系數(shù)矩陣，見式(3)。

(3)

式中，rij(i，j=1，2，…，p)為原變量xi與xj的相關(guān)系數(shù)，rij=rji，其計算公式為式(4)。

(4)

4) 計算特征值與特征向量。

①解特征方程|λI-R|=0，常用雅可比法(Jacobi)求出特征值，并使其按大小順序排列λ1≥λ2≥…≥λp≥0。

在實際工作中，主成分個數(shù)的多少取決于能夠反映原來變量85%以上的信息量為依據(jù)，即選擇累計貢獻率不小于85%的主成分個數(shù)。

6) 建立回歸模型。在確定主成分指標(biāo)后，通過選取的主成分建立回歸模型(式(5))。

Y=a0+a1Z1+a2Z2+…+amZm

(5)

通過上述計算得到了關(guān)于Z1，Z2，…，Zm的主成分回歸，將主成分式子代入得到關(guān)于原指標(biāo)X1，X2，X3，…，Xp的回歸方程(式(6))。

Y=b0+b1X1+b2X2+…+bmXp

(6)

2 實例分析

2.1 研究區(qū)概況

蔣莊煤礦位于山東省西南部的滕南礦區(qū)，北與程樓斷層相接，南接李集斷層，西鄰高廟斷層、徐莊斷層，東連尹家洼斷層，形成一斷裂構(gòu)造發(fā)育，以地塹、地壘為主要特點的寬緩褶皺區(qū)，構(gòu)造發(fā)育，斷裂構(gòu)造眾多，構(gòu)造中等偏復(fù)雜。地層屬華北型沉積，廣為第四系覆蓋，未見基巖露頭。井田主要含水層為二疊系山西組3煤層頂板砂巖含水層、石炭系太原組第三層石灰?guī)r含水層、石炭系太原組第十層灰?guī)r含水層、第四系砂礫層、上侏羅統(tǒng)礫巖含水層、二疊系石盒子組頂部砂巖含水層、二疊系石盒子組底部砂巖含水層、石炭系本溪組第十四層石灰?guī)r含水層及奧陶系石灰?guī)r含水層，富水性不均一[20]。

本文所采用數(shù)據(jù)來源于蔣莊煤礦及其四周緊鄰的多個煤礦，如崔莊煤礦、高莊煤礦、付村煤礦等，所選取煤礦均位于滕州礦區(qū)，在巖體力學(xué)參數(shù)、地質(zhì)構(gòu)造及地層分布上有很大的相似之處，因此在對導(dǎo)水裂縫帶高度預(yù)測的研究過程中，可以利用這些礦井長期生產(chǎn)實踐中得到的影響導(dǎo)水裂縫帶高度因素的數(shù)據(jù)加以定量分析，既對蔣莊煤礦導(dǎo)水裂縫帶高度的研究具有很高的實用價值，也為其他相似煤礦的導(dǎo)水裂縫帶高度預(yù)測工作提供了一種可供借鑒的方法。

2.2 主回歸方程的建立

在研究蔣莊煤礦導(dǎo)水裂縫帶高度發(fā)育程度的過程中，參照了蔣莊煤礦及其周邊相鄰?fù)粎^(qū)域的12個生產(chǎn)礦井的相關(guān)資料，由于這些礦井具有基本相同的頂板管理方法，因此對于導(dǎo)水裂縫帶具有主要影響的因子確定為以下6個主要因素：采深、煤層厚度、工作面傾向長度、巖性參數(shù)(單軸抗壓強度)、頂板厚度及傾角。具體的導(dǎo)水裂縫帶高度影響因素量化數(shù)值見表1。

表1 導(dǎo)水裂縫帶高度主要影響因素量化表Table 1 Quantification table of main factors affecting the height of water-conducting fracture zone

對于上述14個工作面所選取的7個指標(biāo)分別依次確定為x1、x2、x3、x4、x5、x6、x7，由于所選指標(biāo)的量綱不同，因此需要對各指標(biāo)的量綱進行消除，使數(shù)據(jù)具有可比性，即對數(shù)據(jù)進行標(biāo)準(zhǔn)化處理。為了研究主成分回歸分析對于導(dǎo)水裂縫帶高度的預(yù)測可靠性，將前12個工作面的數(shù)據(jù)作為分析數(shù)據(jù)，將蔣莊煤礦2個工作面數(shù)據(jù)作為預(yù)測結(jié)果的驗證數(shù)據(jù)。分析數(shù)據(jù)構(gòu)成的矩陣為矩陣A，見式(7)。

第一步：利用MATLAB數(shù)據(jù)處理軟件對式(7)做標(biāo)準(zhǔn)化處理得到標(biāo)準(zhǔn)化矩陣B，見式(8)。

標(biāo)準(zhǔn)化處理之后的7個無量綱指標(biāo)分別表示為

第二步：經(jīng)過標(biāo)準(zhǔn)化處理后數(shù)據(jù)的相關(guān)系數(shù)矩陣用Rp×p=(rij)p×p表示。由于該矩陣為對稱矩陣，因此可以用矩陣B計算得到的相關(guān)系數(shù)矩陣C，見式(9)。

第三步：求出相關(guān)系數(shù)矩陣的特征向量和特征值，利用MATLAB軟件計算得到具體數(shù)值，即矩陣D，見式(10)。

第四步：確定主成分，見式(11)。

通過累計貢獻率的大小來選擇m(m

表2 主成分特征值和貢獻率Table 2 Eigenvalue and contribution rate value ofprincipal component

(7)

(8)

(9)

(10)

(11)

(12)

通過上述四個主成分的基本線性組合可知：第一主成分F1可看作x2的變量，因為其在式中的系數(shù)絕對值大于其他變量的系數(shù)，即反映了煤層厚度對導(dǎo)水裂縫帶高度的影響；第二主成分F2可看作x4和x6的綜合變量，反映了巖性、傾角對導(dǎo)水裂縫帶高度的影響；F3可看作x5的變量，反映頂板厚度對導(dǎo)水裂縫帶高度的影響；F4可看作x1和x3的變量，反映采深、傾向長度對導(dǎo)水裂縫帶高度的影響。四個主成分的得分情況見表3。

表3 主成分分布表Table 3 Principal component distribution table

取前四個主成分建立回歸方程(式(13))。

y′=a0+a1F1+a2F2+a3F3+a4F4

(13)

采用最小二乘估計可得式(14)。

y′=-0.538 8F1-0.371 1F2+

0.06F3+0.445 4F4

(14)

對于標(biāo)準(zhǔn)化數(shù)據(jù)，把各主成分的回歸系數(shù)還原成關(guān)于原自變量的回歸系數(shù)，經(jīng)計算得到式(15)；將標(biāo)準(zhǔn)化數(shù)據(jù)還原成原始數(shù)據(jù)，經(jīng)過計算得到關(guān)于原始數(shù)據(jù)的回歸系數(shù)見表4?；貧w方程式見式(16)。

表4 原始數(shù)據(jù)回歸系數(shù)Table 4 Raw data regression coefficient

(15)

(16)

2.3 回歸模型檢驗分析

2.3.1 回歸顯著性檢驗

將顯著水平α的值定為0.05，通過MATLAB軟件計算可以得到F的值為37.646 9，可知F>F1-α，表示拒絕H0，因此認(rèn)為此差別不大，可能僅由抽樣誤差所致，存在實驗因素不同造成誤差的可能，故在統(tǒng)計上成立，回歸方程顯著；同時計算相關(guān)系數(shù)r2為0.978 3，數(shù)值非常接近1，說明回歸方程顯著；與F對應(yīng)的概率p為0.000 5，p<α，回歸模型成立。

2.3.2 回歸系數(shù)顯著性檢驗

通過比較各個變量的t值及其相應(yīng)的概率，進行t檢驗，如果相應(yīng)的概率小于給定的顯著水平，則Xi對Y顯著的線性作用顯示為0，否則為1。采用MATLAB軟件得到的計算結(jié)果見表5。由表5可知，選取的自變量對因變量有顯著的線性作用。

表5 Xi對Y的線性作用Table 5 Linear effect of Xi on Y

2.3.3 殘差分析

計算得到用于主成分回歸分析的14個工作面導(dǎo)水裂縫帶高度的殘差依次分別為：2.155 6，-2.709 3，4.650 8，-2.033 0，-3.498 9，0.018 4，0.513 1，-1.694 2，-1.356 1，-0.278 2，-0.306 7，4.538 7。采用MATLAB軟件做出殘差圖，如圖2所示。由圖2可知，除個別數(shù)據(jù)距離零點較遠(yuǎn)之外，其余數(shù)據(jù)的殘差離零點均較近，且所有殘差的置信區(qū)間均包含零點，這說明回歸模型(式(17))能較好地擬合原始數(shù)據(jù)。

圖2 回歸方程殘差圖Fig.2 Regression equation residual plot

y=76.370 7-0.023 3x1+6.46x2-0.067 8x3-0.739 3x4-0.495 1x5+0.066x6

(17)

3 蔣莊煤礦導(dǎo)水裂縫帶發(fā)育高度預(yù)測

3.1 主成分分析回歸模型

通過上述對蔣莊煤礦周邊生產(chǎn)礦井影響導(dǎo)水裂縫帶高度的因素進行主成分回歸分析，以及回歸模型檢驗分析，得出了回歸顯著性較強、相對準(zhǔn)確合理的回歸方程。將已有的蔣莊煤礦3下803工作面、3上603工作面導(dǎo)水裂縫帶的量化指標(biāo)代入回歸方程，可以得到導(dǎo)水裂縫帶高度，見表6。

表6 實測值與預(yù)測值對比Table 6 Comparison of measured value and predicted value

3.2 經(jīng)驗公式法

由于蔣莊煤礦工作面頂板多為粗、中、細(xì)砂巖、粉砂巖和泥巖，屬中硬偏堅硬巖性，因此工作面頂板定為中硬巖層情況考慮，按照國家煤礦安監(jiān)局制定的《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)范》的規(guī)定，本次蔣莊煤礦兩工作面的冒落帶及導(dǎo)水裂縫帶高度計算公式選取中硬巖性公式[9](式(18)和式(19))。

(18)

(19)

式中，∑m為累計采厚。

已知蔣莊煤礦3上工作面煤層采厚3.8 m，3下工作面煤層采厚亦為3.8 m，按式(18)和式(19)計算其冒落帶及導(dǎo)水裂隙帶高度見式(20)和式(21)。

(20)

(21)

通過計算可以得出，按照中硬巖層考慮預(yù)計采3上(3下)工作面煤層時，兩帶高度分別為H冒=8.1～12.5 m，H裂=33.7～44.9 m，實際上工作面頂板為中偏硬巖層，為了使觀測方案安全可靠，裂隙帶探測高度定為60.0 m，見表7。

表7 實測值與預(yù)測值對比Table 7 Comparison of measured value and predicted value

3.3 結(jié)果分析

主成分回歸分析法可消除回歸分析中出現(xiàn)的指標(biāo)間多重共線關(guān)系影響，使所建回歸模型更符合實際情況。通過實測數(shù)據(jù)可知蔣莊煤礦3下803工作面和3上603工作面導(dǎo)水裂縫帶高度分別為54.60 m、47.75 m，基于山東滕州礦區(qū)12個煤礦的實測數(shù)據(jù)建立的主成分回歸模型，對兩工作面導(dǎo)水裂縫帶高度預(yù)測分別為56.239 1 m、49.102 1 m，預(yù)測相對誤差分別為3.00%、2.83%；運用經(jīng)驗公式法對兩工作面導(dǎo)水裂縫帶高度預(yù)測結(jié)果均定為60.00 m，預(yù)測相對誤差分別為-9.89%、-25.65%。由以上預(yù)測結(jié)果可知，主成分回歸分析模型方法預(yù)測回采工作面的導(dǎo)水裂縫帶發(fā)育高度較經(jīng)驗公式法具有較高的準(zhǔn)確性，可以更好地預(yù)測導(dǎo)水裂縫帶發(fā)育高度。

4 結(jié) 論

1) 主成分回歸分析法對于數(shù)據(jù)量較大的實際問題具有較好的數(shù)據(jù)處理能力，通過降維的處理思想，可以通過將原始具有重疊關(guān)系的數(shù)據(jù)重新組合優(yōu)化，通過貢獻率的大小選定多個主成分，基本代表原有數(shù)據(jù)的所有信息。

2) 基于山東滕州礦區(qū)12個煤礦的實測數(shù)據(jù)建立的主成分回歸模型，對蔣莊煤礦3下803工作面和3上603工作面導(dǎo)水裂縫帶高度預(yù)測結(jié)果分別為56.239 1 m、49.102 1 m，預(yù)測相對誤差分別為3.00%、2.83%；運用經(jīng)驗公式法對兩工作面導(dǎo)水裂縫帶高度預(yù)測結(jié)果均定為60.00 m，預(yù)測相對誤差分別為-9.89%、-25.65%。

3) 本文受驗證數(shù)據(jù)的數(shù)量限制，未能夠排除數(shù)據(jù)量較少帶來的誤差，但在已有數(shù)據(jù)的基礎(chǔ)上表明，主成分回歸分析模型對于蔣莊煤礦的導(dǎo)水裂縫帶高度預(yù)測符合實際，具有較好的預(yù)測精度，可以作為一種導(dǎo)水裂縫帶高度預(yù)測的輔助方法。