基于主成分回歸分析法的導水裂縫帶高度預測

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

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

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

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

1 方法和原理

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

1) 構建數據矩陣并標準化處理。設有n個礦區樣品，每個礦區樣品有p個主要影響因素，將原始數據寫成矩陣，見式(1)。

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

(1)

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

(2)

3) 計算相關系數矩陣，見式(3)。

(3)

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

(4)

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

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

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

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

Y=a0+a1Z1+a2Z2+…+amZm

(5)

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

Y=b0+b1X1+b2X2+…+bmXp

(6)

2 實例分析

2.1 研究區概況

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

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

2.2 主回歸方程的建立

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

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

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

第一步：利用MATLAB數據處理軟件對式(7)做標準化處理得到標準化矩陣B，見式(8)。

標準化處理之后的7個無量綱指標分別表示為

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

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

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

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

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

(7)

(8)

(9)

(10)

(11)

(12)

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

對于標準化數據，把各主成分的回歸系數還原成關于原自變量的回歸系數，經計算得到式(15)；將標準化數據還原成原始數據，經過計算得到關于原始數據的回歸系數見表4。回歸方程式見式(16)。

表4 原始數據回歸系數Table 4 Raw data regression coefficient

(15)

(16)

2.3 回歸模型檢驗分析

2.3.1 回歸顯著性檢驗

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

2.3.2 回歸系數顯著性檢驗

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

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

2.3.3 殘差分析

計算得到用于主成分回歸分析的14個工作面導水裂縫帶高度的殘差依次分別為：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可知，除個別數據距離零點較遠之外，其余數據的殘差離零點均較近，且所有殘差的置信區間均包含零點，這說明回歸模型(式(17))能較好地擬合原始數據。

圖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 蔣莊煤礦導水裂縫帶發育高度預測

3.1 主成分分析回歸模型

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

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

3.2 經驗公式法

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

(18)

(19)

式中，∑m為累計采厚。

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

(20)

(21)

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

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

3.3 結果分析

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

4 結 論

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

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

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