高強度開采地表損傷程度分類判別與控制研究

李曉斌，李全生,2,3,4，韓鵬華，許旭輝，皇甫龍

(1.中國礦業(yè)大學(北京) 能源與礦業(yè)學院，北京 100083；2.煤炭開采水資源保護與利用國家重點實驗室，北京 102209；3.國家能源投資集團有限責任公司，北京 100011；4.北京低碳清潔能源研究院，北京 102211)

神東礦區(qū)是我國西部礦區(qū)高強度開采的典型代表，工作面具有埋深小、采高與開采尺寸大、推進速度快、開采設備機械化程度高等特征，覆巖破壞與地表損傷比較嚴重，容易產生大量覆巖裂隙與地裂縫[1-2]。高強度開采引起地下水資源系統(tǒng)發(fā)生劇烈變化，影響地表土壤水分、有機質與礦物質含量，加之西部地區(qū)自身生態(tài)環(huán)境脆弱，造成煤炭高效開采與生態(tài)保護的矛盾日益凸顯[3-4]。因此，從工作面、覆巖與地表3個角度出發(fā)，研究高強度采動損傷定量評價與減損技術迫在眉睫。

采動損傷是指煤層開采后改變圍巖應力狀態(tài)，形成高地應力區(qū)域，導致覆巖發(fā)生變形、破壞與運移，最終傳導至地表，產生開采沉陷和生態(tài)損傷的過程[5]。由此可見，采動損傷傳導涉及煤炭開采的全生命周期。我國采礦科技工作者針對高強度采動損傷也開展了一系列研究，李全生[6]等從開采源頭與采動損傷能量角度出發(fā)，建立高強度開采損傷傳導模型，并在上灣煤礦12401工作面得到應用；何祥[7]等基于定義覆巖損傷度，探究了工作面長度與推進速度對覆巖損傷程度的影響；胡振琪[8]等基于現(xiàn)場實測數(shù)據(jù)，給出了風積沙區(qū)工作面埋深、推進速度與地裂縫發(fā)育周期的函數(shù)模型；此外還有眾多采礦學者[9-11]以工作面開采參數(shù)作為指標，采用隨機森林算法、層次分析法、遺傳算法等數(shù)學方法對覆巖損傷高度與開采沉陷進行預測。上述研究主要集中在工作面開采參數(shù)對覆巖或地表損傷的影響，未將工作面開采參數(shù)、覆巖與地表損傷三者緊密聯(lián)系起來，且鮮有研究涉及采動損傷程度定量評價與分類判別。目前，主要通過改變開采方法(限厚開采、條帶開采、充填開采、采-充-留協(xié)調開采)、留設煤柱、覆巖離層注漿技術等控制采動損傷[12-15]。整體治理手段比較單一，部分方法實施成本較高、效率低，且損傷程度影響因素眾多，不同開采條件下?lián)p傷程度與生態(tài)脆弱度不盡相同。在實際應用過程中，需針對具體地質條件，綜合考慮各種采礦與生態(tài)環(huán)境因素，選取合適的控制方法，才能降低損傷程度，達到保護生態(tài)環(huán)境的目標。

鑒于此，筆者從工作面開采參數(shù)、覆巖硬度系數(shù)與地表松散層等方面建立PCA-DDA數(shù)學評價模型，借助主成分分析法獲得高強度采動損傷類型，進一步通過距離判別分析法驗證模型可靠性。根據(jù)采動損傷類型與現(xiàn)場采礦地質條件，在上灣煤礦12401工作面提出“大工作面快速推進法、覆巖局部注漿加固法與裂縫分區(qū)治理法”三步法，實現(xiàn)了源頭減損(工作面)、傳導過程(覆巖)、末端治理(地表)的開采全周期采動損傷控制。

1 神東礦區(qū)開采強度特征

神東礦區(qū)是我國最大的煤炭地下開采礦區(qū)，各煤礦主要集中分布在山西、陜西與內蒙，礦區(qū)分布面積廣，東西方向寬度達55 km，南北方向長度可至90 km，已探明煤炭儲量超過2 000億t，同時具有上灣煤礦、補連塔煤礦、布爾臺煤礦、大柳塔煤礦等千萬噸大型礦井，整體產能在2 億t以上。礦區(qū)煤層地質條件簡單，埋藏較淺，工作面開采尺寸與采高大，推進速度較快，地表容易產生連續(xù)與非連續(xù)變形，整體開采強度比較大，神東礦區(qū)開采強度分布情況如圖1[16]所示。

圖1 神東礦區(qū)開采強度分布[16] Fig.1 Distribution of mining intensity in Shendong mining area[16]

2 主成分分析法與距離判別分析法原理

2.1 主成分分析法

由于高強度采動損傷程度影響因素眾多，包括工作面開采參數(shù)，覆巖綜合硬度與地表松散層厚度等，故引入主成分分析法(PCA)，它是將多個變量轉化為幾個綜合指標的一種方法。在多變量分析過程中，經常由于變量多且變量間存在相關性，造成數(shù)據(jù)信息發(fā)生重疊。在維數(shù)較多的情況下，變量分析就變得更加復雜，且部分數(shù)據(jù)處理系統(tǒng)對指標數(shù)量有一定限制。主成分分析采用一種降低維度的辦法[17]，通過求解幾個主成分因子來最大程度地反映原有變量信息，且綜合因子之間沒有相關性，實現(xiàn)簡化變量的目標，具體求解過程包括下述幾個步驟。

2.1.1 相關系數(shù)矩陣

式中，R為相關系數(shù)矩陣；ri為原指標變量xi與yi之間的相關系數(shù)，其中，xi，yi為樣本值。

ri的計算方法為

2.1.2 特征值和特征向量

相關系數(shù)矩陣R的特征方程為

式中，μ為特征值，取值為μ1，μ2，…，μp，且滿足μ1≥μ2≥…≥μp；E為單位矩陣。

特征向量計算方法為

式中，α為特征向量。

根據(jù)式(4)求出特征值(μ1，μ2，…，μp)對應的特征向量值為α1，α2，…，αp。

2.1.3 主成分貢獻率和累計貢獻率

主成分貢獻率計算方法為

式中，V為主成分貢獻率；μi為第i個特征值；為p個特征值之和。

主成分累計貢獻率計算方法為

式中，V′為主成分累計貢獻率；為前i個特征值之和。

通常認為主成分累計貢獻率超過85%，就可認定前m個主成分因子包括原有變量的絕大部分信息。μ1，μ2，…，μm對應第1，第2，…，第m個主成分因子，其中m≤p。

2.1.4 主成分因子載荷矩陣與主成分得分

主成分因子載荷計算方法為

式中，li為第i個主成分因子載荷；zi為第i個得分情況；αi為第i個特征向量值。

得到主成分因子載荷后，求解主成分得分矩陣G為

根據(jù)特征值與主成分因子得分，可求得綜合得分，具體計算方法為

式中，Z為綜合得分；G1，G2，…，Gm對應第1，第2，…，第m個主成分因子得分。

2.2 距離判別分析法

距離判別分析法(DDA)是根據(jù)多變量確定已有樣本類別，用來判別未知樣本類別的一種分析方法，通常用誤判率來表示判別的可靠性[18]。假設Q1與Q2為2個不同總體，X1，X2，…，Xm是來自Q1的m個訓練樣本，Y1，Y2，…，Yn是來自Q2的n個訓練樣本，如果屬于Q1的樣本被誤判為屬于Q2的個數(shù)為u個，屬于Q2的樣本被誤判為屬于Q1的個數(shù)為v個，則誤判率η可表示為

式中，u，v為誤判樣本數(shù)。

3 高強度開采損傷評價指標

煤層高強度開采后從直接頂開始上覆巖層逐次發(fā)生彎曲破斷，采動損傷逐層由下向上傳導直至地表，引起地表產生一系列連續(xù)與非連續(xù)變形，威脅西部脆弱的生態(tài)環(huán)境。鑒于此，高強度開采損傷評價指標體系的確定要包含工作面、覆巖與地表整個開采全周期。從工作面角度考慮，選取工作面寬度、工作面推進長度、采高、推進速度、埋深與寬深比等指標，覆巖方面引入覆巖硬度系數(shù)，地表方面包括松散層厚度與沙基比(松散層厚度與基巖厚度之比)。關于工作面開采參數(shù)對地表損傷的影響的研究較多，在此不再贅述。筆者重點研究松散層厚度、覆巖硬度系數(shù)與地表下沉系數(shù)之間的內在聯(lián)系。

3.1 松散層厚度與地表下沉系數(shù)的關系

由于松散層具有力學強度與硬度小，容易跟隨巖層發(fā)生流動的特征，在開采下沉過程中，松散層不僅能夠減小覆巖綜合硬度，還會增加覆巖載荷，促使采空區(qū)上方巖層處在高載荷應力中[19]。由于松散層的特性，在開采引起覆巖斷裂時，會很快傳導至地表，產生地裂縫等采動損害。如圖2所示，基于高強度開采工作面部分數(shù)據(jù)擬合，結果表明松散層厚度與埋深比值與地表下沉系數(shù)呈二次函數(shù)關系，且在松散層厚度與埋深比值為0.25附近出現(xiàn)轉折。

圖2 松散層厚度與埋深比值與地表下沉系數(shù)的關系 Fig.2 Relationship between the ratio of loose layer thickness to buried depth and surface subsidence coefficient

這是由于當比值小于0.25時，沒有形成松散層拱結構，松散層以載荷形式作用在基巖上面，隨著松散層厚度增加，地表下沉系數(shù)增大；當比值大于0.25時，在厚度較大的松散層中形成拱結構，對地表下沉起到一定控制作用，地表下沉系數(shù)減小。

3.2 覆巖硬度系數(shù)與地表下沉系數(shù)的關系

覆巖硬度與巖石中所含礦物成分、巖石類型及單軸抗壓強度等息息相關。在高強度開采條件下，堅硬頂板懸露面積大，不易發(fā)生垮落，當?shù)竭_巖層的破壞強度時，覆巖發(fā)生切落式破壞，引起地表突然陷落。在覆巖硬度較小時，隨著工作面的開采，地表逐漸出現(xiàn)緩慢式下沉。

覆巖硬度系數(shù)[20]計算公式為

根據(jù)式(11)可得K∈[0,1]，當K=0時，上覆全部為松散層；當K=1時，上覆全部為巖層；當K∈（0,1）時，煤層上方存在巖層與松散層互層。

圖3為覆巖硬度系數(shù)與地表下沉系數(shù)的關系。由圖3可知，在覆巖硬度系數(shù)較小時，地表下沉系數(shù)呈上升趨勢；在覆巖硬度系數(shù)較大時，起到相當于關鍵層的作用，從而控制地表下沉的效果，地表下沉系數(shù)逐漸減小。

圖3 覆巖硬度系數(shù)與地表下沉系數(shù)的關系 Fig.3 Relationship between overburden hardness coefficient and surface subsidence coefficient

4 數(shù)學評價模型的建立

通過收集神東礦區(qū)25個高強度開采工作面相關數(shù)據(jù)，對工作面開采參數(shù)、覆巖硬度系數(shù)、地表松散層厚度等數(shù)據(jù)進行了匯總，見表1。

4.1 Pearson相關系數(shù)矩陣的確定

由于各個指標量綱不同，在進行主成分分析之前，需要對原始數(shù)據(jù)進行標準化處理，將收集的數(shù)據(jù)導入DPS數(shù)據(jù)處理系統(tǒng)得到相關系數(shù)矩陣，見表2。

表1 不同采礦地質條件下高強度開采工作面數(shù)據(jù) Table 1 Data of high-intensity mining face under different mining geological conditions

表2 Pearson相關系數(shù)矩陣 Table 2 Pearson correlation coefficient matrix

由表2可以看出，寬深比(X6)與工作面寬度(X1)相關系數(shù)為0.887，松散層厚度(X7)與沙基比(X8)、覆巖硬度系數(shù)(X9)相關系數(shù)分別為0.880，-0.865，覆巖硬度系數(shù)(X9)與沙基比(X8)相關系數(shù)為-0.853，這些變量之間存在顯著性較高的相關關系，進一步反映原始指標彼此之間存在信息重疊的現(xiàn)象，因此很有必要運用主成分分析法對評價指標進行維度降低處理，獲取新的互不相關指標。

4.2 主成分因子分析

運用Data Processing System(DPS)軟件對數(shù)據(jù)進行主成分分析，得到KMO(Kaiser-Meyer-Olkin)檢驗值為0.526，大于0.500，Bartlett球形檢驗p=0，表明采用主成分分析法是可行的。特征值在一定情況下能夠反映主成分對指標的影響程度。根據(jù)主成分分析法，選取累計百分率超過85%即能概括原變量的大部分信息，DPS系統(tǒng)默認取90%。特征值分布見表3。

表3 特征值分布 Table 3 Distribution of eigenvalues

由表3可知，前5個主成分特征值分別為3.757，2.051，1.403，0.800，0.615，方差累計百分率為95.857%。由圖4可知，當主成分數(shù)超過5，特征值逐漸趨向收斂狀態(tài)。故選擇前5個主成分能夠較好地反映原有變量信息，即該5個主成分能解釋原變量95.857%的內容。基于此，獲得這5個主成分的因子載荷矩陣，見表4。

該矩陣揭示了主成分因子與原有變量之間的相關程度，矩陣中數(shù)值絕對值越大，表示主成分因子與原變量相關程度越高。第1個主成分因子與X9最相關，在很大程度上解釋了覆巖硬度系數(shù)，從第2，3，4，5個主成分因子來看，X7(0.724)，X6(0.778)，X3(-0.737)，X4(-0.502)的絕對值最大，分別反映了與松散層厚度、寬深比、采高和推進速度的相關程度，其中推進速度的相關程度較小。

圖4 主成分分析結果 Fig.4 Results of principal component analysis

表4 因子載荷矩陣 Table 4 Factor load matrix

然而，從各個主成分因子反映的原變量來看，用單一主成分來確定高強度開采沉陷程度的大小是比較困難的，因此采用前5個主成分作為綜合評估指標，以各個主成分方差貢獻率為權重建立綜合評價數(shù)學模型[21]為

基于綜合評價數(shù)學模型，獲得表1中25個工作面高強度采動損傷程度的排名情況，見表5。

表5 高強度開采損傷程度綜合得分與排名 Table 5 Comprehensive score and ranking of damage degree of high-intensity mining

續(xù) 表

根據(jù)排名與綜合得分情況得到圖5，將高強度采動損傷分為Ⅰ型、Ⅱ型、Ⅲ型、Ⅳ型，分別對應極高、高、中、低4個損傷程度，其中Ⅰ型包括4，5，24，18，22樣本組，綜合得分＞0.5；Ⅱ型包括17，25，21，1，20，2，13樣本組，綜合得分0～0.5；Ⅲ型包括6，16，11，3，8，7，12，14樣本組，綜合得分-0.5～0；Ⅳ型包括10，15，9，23，19樣本組，綜合得分＜-0.5。此方法可實現(xiàn)高強度開采損傷程度的定量判別，為地表損傷治理與生態(tài)修復提供科學指導。

圖5 綜合得分與分類 Fig.5 Comprehensive score and classification chart

4.3 數(shù)學模型驗證

采用DDA對定量分類結果進一步驗證，以PCA定量評價得到的Ⅰ型、Ⅱ型、Ⅲ型、Ⅳ型作為4個不同總體，學習樣本選取圖中25組散點數(shù)據(jù)，同時為降低訓練樣本與測試樣本個數(shù)對分類結果的影響，取測試樣本占學習樣本的20%，選擇5組情況進行學習，訓練樣本與測試樣本比例分別為24∶1，23∶2，22∶3，21∶4，20∶5，見表6。當訓練樣本與測試樣本比例為23∶2，誤判概率最大為4.35%，可能與樣本數(shù)量少、指標選取不全面、模型泛化能力不強等因素有關，但該模型整體判對概率超過95%，表明在高強度開采損傷程度評價方面是可行的。

表6 不同樣本比例的5組學習情況 Table 6 Learning of 5 groups with different sample proportions

5 現(xiàn)場應用

5.1 工程背景

以神東礦區(qū)上灣煤礦12401工作面為工程背景，工作面寬度300 m，推進長度5 429 m，采高8.6 m，平均推進速度13.6 m/d，埋深184 m，寬深比1.63，松散層厚度62 m，沙基比0.51，覆巖硬度系數(shù)0.663，采用PCA-DDA數(shù)學模型判定采動損傷類型為Ⅱ型，相應損傷程度為高。由于開采強度較大，地表損傷控制層發(fā)生破斷，傳導至地表產生大量地裂縫及塌陷槽，圖6[22]為地裂縫分布平面圖，不同顏色代表不同裂縫寬度。

圖6 地裂縫分布平面[22] Fig.6 Ground fissures distribution plan[22]

5.2 高強度采動損傷控制

針對上灣煤礦12401工作面地表損傷程度較高的特點，綜合采用大工作面快速推進法、覆巖局部注漿加固法與裂縫分區(qū)治理法，實現(xiàn)了源頭減損(采場)、傳導過程(覆巖)、末端治理(地表)的開采全周期采動損傷控制。

5.2.1 大工作面快速推進法

大工作面采動影響范圍內地表變形很小，該工作面在開采初期處于設備調試階段，平均推進速度＜4 m/d，覆巖及地表損傷嚴重，后期平均推進速度達到13.8 m/d，覆巖卸壓不充分，基本頂來壓步距增加，在時間和空間上損傷傳導表現(xiàn)出一定滯后性，損傷程度有所減小。

圖7[7]為不同推進速度下地表中部裂縫發(fā)育情況，表明加快推進速度可使地表損傷程度更為緩和。

圖7 不同推進速度下地裂縫發(fā)育特征[7] Fig.7 Development characteristics of ground fissures at different advancing speeds[7]

5.2.2 覆巖局部注漿加固法

覆巖局部注漿加固法是通過在采空區(qū)垮落帶與離層裂隙帶注漿實現(xiàn)的，垮落帶注漿能夠維持基本頂穩(wěn)定性，減小頂板下沉量與垮落帶上方巖層移動變形程度。離層裂隙帶注漿可以支撐關鍵層的下沉，進一步加固關鍵層，將采動損傷阻隔在覆巖中，降低覆巖破壞高度，促進“三帶”模式的出現(xiàn)，減弱地表損傷程度。

圖8為注漿前后裂縫寬度隨時間的變化曲線，由圖8可知注漿后裂縫發(fā)育更加平緩，裂縫寬度顯著減小，能有效實現(xiàn)采動損傷控制。

圖8 注漿前后裂縫寬度變化特征 Fig.8 Variation characteristics of fracture width before and after grouting

5.2.3 地裂縫分區(qū)治理法

針對開采邊界永久裂縫采取源頭減損與加固關鍵層的控制方法，對于中部動態(tài)裂縫，由于關鍵層破斷巖塊回轉，能閉合的裂縫讓其自修復，對于不能閉合的裂縫采取人工填埋的方法。

6 結 論

(1) 基于工作面開采參數(shù)、覆巖硬度系數(shù)與地表松散層厚度等9個評價指標，提出一種高強度采動損傷定量評價數(shù)學模型——PCA-DDA模型，該模型可以在保留原變量信息的情況下，轉化5個主成分因子分析，實現(xiàn)高強度采動損傷定量評價與分類判別，模型判對概率高達95%以上。

(2) 基于PCA-DDA模型，得到高強度采動損傷程度的排名與綜合得分情況，在此基礎上將高強度采動損傷分為Ⅰ型、Ⅱ型、Ⅲ型、Ⅳ型，分別對應極高、高、中、低4個損傷程度，此分類結果為采動損傷治理與生態(tài)修復研究提供科學指導。

(3) 將PCA-DDA模型應用于上灣煤礦12401工作面，結果表明該工作面采動損傷程度較高，綜合采用大工作面快速推進法、覆巖局部注漿加固法與裂縫分區(qū)治理法，實現(xiàn)了源頭減損(工作面)、傳導過程(覆巖)、末端治理(地表)的開采全周期采動損傷控制。

需要說明的是，由于數(shù)據(jù)來源有限，筆者只針對部分高強度開采工作面及評價指標初步建立數(shù)學分類判別模型，沒有考慮其他因素的影響，如導水裂隙帶發(fā)育高度、關鍵層位置及數(shù)量、地裂縫參數(shù)等，其將在后續(xù)工作中做進一步研究。