磁窯溝煤礦黃土溝壑地貌開采沉陷規律研究

楊玉龍，宋選民

(1.山西河曲晉神磁窯溝煤業有限公司, 山西 忻州 036500；2.太原理工大學原位改性采礦教育部重點實驗室，太原 030024)

目前，我國能源結構中煤炭仍占據著主導地位，絕大部分煤炭資源來自于井工煤礦，其工作面的回采在不同程度上引起了地表沉陷變形[1]，使得煤炭資源開發與礦區環境保護之間的矛盾日益凸顯[2-3]。河東煤田是山西省重要產煤基地之一，煤炭資源賦存較淺、煤層厚度較大、上覆黃土溝壑錯綜復雜，工作面回采后對礦區地表產生較大影響，制約了礦區經濟社會的可持續發展[4-5]。

針對黃土溝壑地貌的開采沉陷問題，國內外眾多專家學者取得了大量有益研究成果。ALAEE et al[6]基于模糊理論，建立了有關地表開采沉陷的模糊隸屬度函數，并于伊朗Tabas煤礦進行了現場驗證。KOPEC et al[7]對波蘭上西西里亞煤田開采引發的地表沉陷進行了實測研究，分析了其對地表道路、房屋等建筑物的影響。魏國武[8]以大佛寺煤礦30201工作面為背景，總結得出礦山開采沉陷影響因子，構建了礦山開采沉陷危險性系數計算模型。吳文敏等[9]研究了黃土高原丘陵地貌對地表移動變形的影響，得出受采動影響溝谷兩側向下坡方向滑移，溝谷底部的地表會發生抬升，滑移分量比值與坡度存在線性遞減的關系。趙兵朝等[10]研究了生態脆弱區厚松散層薄基巖、黃土溝壑及山區丘陵下煤層開采的地表損害特征，并以基采比和基載比為關鍵參數，揭示了生態脆弱區煤層開采損害特征形成的機理。陳超等[11]根據采動地裂縫的力學成因進行分類，系統梳理了不同類型采動地裂縫形成機理研究進展，分析了采動地裂縫的影響因素與研究方法，討論了相關的關鍵科學問題。張磊等[12]運用數值模擬方法，研究了黃土覆蓋區域煤層開采引起的地表沉陷特征，認為受采動及上覆巖層自重的影響，地表沉陷變形隨采掘范圍的增加而加大。賀國偉等[13]基于支持向量機，構建了適用于黃土礦區的最大下沉預計模型，提高了黃土地表最大下沉值預計精度和可靠性。李春意等[14]采用現場實測、理論分析和離散元數值模擬相結合的方法，分析了地表沉陷的基本特征，構建了基于坐標變換的傾向地表下沉和水平移動擬合函數模型，確定了地表沉陷的角量參數和預計參數。溫建生[15]以柳灣礦61115工作面地表沉陷實測資料為基礎，分析了地表移動變形動靜態規律，得出地表下沉集中在活躍期并且下沉速度較大，基巖移動角偏大而綜合邊界角偏小。

綜上所述，對于煤層開采影響下的黃土溝壑地表沉陷預測的研究需要根據礦井實際地質生產條件綜合分析判斷。本文以磁窯溝煤礦13101工作面生產地質條件為工程背景，運用理論分析及數值模擬相結合的方法，分析其開采引起的地表移動變形規律，對該礦厚煤層綜放開采條件下的地表沉陷預判具有一定的指導意義。

1 工程背景

1.1 13101工作面工程地質條件

13101工作面位于磁窯溝礦二水平一盤區，主采13#煤層，位于井田南部，上部為10#煤層采空區，平均間距為20 m；東部為13102工作面(未形成)，以東為13#煤層未開采的實體煤；南部為山西華鹿陽坡泉煤礦，留設20 m保護煤柱，該煤礦目前未開采13#煤層;西部為13#煤層未開采的實體煤，與留設的寺也村村莊保護煤柱相接；北部為131盤區輔助運輸大巷/較大大巷/回風大巷，以北為13#煤一盤區，為未開采的實體煤。工作面寬度135.6 m，長度為1 579.4 m，煤層總厚9.40～11.10 m，平均10.25 m，煤層傾角2.8°～4.6°，平均3.7°，采用綜合機械化放頂煤工藝進行回采。

1.2 13101工作面地表典型特征

13101工作面地面標高+1 026～+1 135 m，工作面標高+889～+921 m，蓋山厚度122～215 m;松散層厚度40.4～108.0 m，平均74.2 m;基巖厚度65.5～119.0 m，平均95.0 m。工作面地表地形為山西西北部黃土高原中-低山區地形，地面植被稀少，沖刷劇烈，溝谷成“V”字形羽毛狀分布，地表大部分為新生界黃土及紅土覆蓋，且陡峭懸崖有較小規模的黃土崩塌現象。地表地形總體北高南低，山坡和溝谷分布廣泛，大部分為旱地和荒地。距切眼1 500～1 580 m地表為溝谷，大氣降水匯集后形成的積水會在短時間內流入北部的磁窯溝-洞溝一帶溝壑，其地表形貌如圖1所示。

圖1 13101工作面地表

2 概率積分法預計參數及計算分析

2.1 相關參數選取

13101工作面地表變形理論分析中，首先需要根據工作面實際生產、地質條件，依據《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規范》的相關規定，對下沉系數、水平移動系數、主要影響角正切等參數進行確定。

1)下沉系數。由于13101工作面上覆10#煤層均已采空，且兩層煤平均間距為20 m，故根據以下公式計算：

qf=(1+a)qc.

(1)

式中:qf為復采下沉系數；a為下沉活化系數，取0.2；qc為初采下沉系數，取0.55。代入磁窯溝礦相關數據可知，13101工作面上覆巖層下沉系數為0.66。

2)水平移動系數b。重復采動條件下，水平移動系數與初次采動相同，即：bf=bc，故取bf為0.25。

3)主要影響范圍角正切tanβf。重復采動時tanβf較初次采動增加0.3～0.8。對于中硬巖層可按式(2)計算：

tanβf=tanβc+0.06236lnh-0.017.

(2)

式中：tanβf為重采時主要影響范圍角正切；tanβc為初采時主要影響范圍角正切，取2.16；h為第二層煤的采深，取平均埋深164.4 m。帶入磁窯溝礦相關數據可知，13101工作面上覆巖層主要影響范圍角正切tanβf為2.48。

2.2 地表變形計算結果

根據以上分析所得相關系數，對13101工作面回采影響下地表變形進行理論預測，包括地表最大下沉值、最大水平移動值、最大傾斜變形值、最大曲率變形值及最大水平變形值，各參數的具體計算過程如下。

1)地表最大下沉值計算。

充分采動

Wcm=dqcosα.

(3)

非充分采動

Wbm=dqncosα.

(4)

式中:Wcm為充分采動條件下地表最大下沉值，mm；Wbm為非充分采動條件下地表最大下沉值，mm；q為充分采動條件下的下沉系數，取0.66；d為煤層法向開采厚度，取值10.25 m；α為煤層傾角，(°)；n為地表充分采動系數,其中

(5)

式中：n1和n3大于1時取1；k1、k3為與覆巖巖性有關的系數，堅硬巖層：k1、k3=0.7；中硬巖層：k1、k3=0.8；軟弱巖層k1、k3=0.9；l1、l3為傾向及走向工作面長度，m；h0為工作面平均開采深度，m。13101工作面寬度135.6 m，長度為1 579.4 m，平均埋深164.4 m，工作面上覆巖層以砂巖為主，為中硬巖層，代入上式計算可知13101工作面地表最大下沉值為5.48 m。

2)地表最大水平移動值計算。沿煤層走向方向上的最大水平移動

Um=bWmax.

(6)

式中，Um為最大水平移動值，mm；b為水平移動系數，取0.25；Wmax為最大下沉值，mm。代入13101工作面地表最大下沉值5.48 m計算可知，沿煤層走向的最大水平移動為1.37 m。

沿煤層傾斜方向的最大水平移動

Um=(b+0.7P)Wmax，

(7)

式中:P為覆巖綜合評價系數；h為表土層厚度，m；α為煤層傾角，(°)。磁窯溝13101工作面，上覆松散層平均厚度為69.4 m，平均埋深164.4 m，煤層平均傾角3.7°，地表最大下沉值為5.48 m。代入上式計算可知13101工作面地表沿煤層傾向方向上的最大水平移動為1.37 m。

3)最大傾斜變形值計算

(8)

4)最大曲率變形值計算

(9)

式中:kmax為最大曲率變形，mm-1。磁窯溝13101工作面地表最大下沉值為5.48 m，帶入上式計算可知，最大曲率變形1.89 mm-1。

5)最大水平變形值計算

(10)

式中:εmax為最大水平變形，mm/m。磁窯溝13101工作面地表最大下沉值為5.48 m，帶入上式計算可知，最大水平變形值0.47 mm/m。

3 黃土溝壑地表三維數值模型

3.1 數值模型及邊界條件

借助Midas-GTS/ANSYS～FLAC3D耦合建模方法,可以實現復雜地質模型的構建，模型可將各巖層、地質構造按照等高線圖真實地表現出來，大大提高了模擬預計分析的準確度和可靠性。選取13101工作面中部作為建模區域，以該區域內勘探鉆孔為巖層賦存依據，由于煤層開采所引起的地表下沉范圍大于工作面長度(136 m)，因此數值模型沿工作面方向分別向外擴展100 m，故數值模型在工作面方向的長度為430 m，工作面推進方向上模型長度為500 m，模型垂直方向按照區內實際鉆孔勘測結果進行建模，共包含巖層30層，計算單元327 495個，計算節點342 380個，如圖2所示。

圖2 三維數值模型

基于上述模型，對13101工作面進行回采，根據實際情況設置邊界條件，模型包含5個位移邊界，即模型四周平面固定對應法向變形，模型底面固定3個主方向的變形。

3.2 數值模擬本構及參數

由于本次數值模擬中所涉及材料均為巖土類，因此選用FLAC3D軟件中內嵌摩爾-庫倫本構，計算過程中需用戶賦值的材料參數包括：容重、抗拉強度、內聚力、內摩擦角、彈性模量及泊松比，由區內鉆孔數據可知，13101工作面上覆巖層大致可分為9種巖性，各巖性及參數匯總見表1所示。

表1 上覆巖層巖性及參數

3.3 數值模擬方案

首先根據前述參數對模型材料進行賦參，計算初始平衡，因區內鉆孔勘測結果中，13#煤層上部存在10-1#煤層采空區，且目前該采空區已壓實穩定，因此將10-1#煤層所對應模型設置為采空區模擬參數，故在本章數值模擬過程中不需對10-1#煤層采空區進行開采模擬。初始平衡后對13#煤層進行分步開挖，初次開挖50 m，以模擬工作面的前期開挖工況;其后，以20 m為開挖步距對剩余工作面進行回采模擬，直至工作面回采至450 m，總計開挖模擬21次。從工作面回采過程中的地表下沉量、水平移動量、剪切應變及剪應變增量4個方面出發，分別提取工作面推進不同長度時的地表各項參數模擬結果，對地表變形機理展開分析。

4 數值模擬結果及分析

數值模擬共開挖21次，由于篇幅所致，以工作面推進50 m、150 m、250 m、350 m及450 m為例，提取各參數模擬結果，對13101綜放工作面回采引起的黃土溝壑地表沉陷變形進行分析。

4.1 地表下沉量

工作面推進過程中，黃土溝壑地表下沉量分布云圖如圖3所示。由圖3可知，隨著工作面推進至50 m、150 m、250 m、350 m以及450 m時，工作面上方地表最大下沉量為1.50 m、5.08 m、6.12 m、6.41 m以及6.55 m，可見隨著工作面的推進，地表下沉量不斷增加，但下沉增速由0～50 m的30.02 mm/m，下降至350～450 m的1.38 mm/m，可見采空區覆巖運移趨于穩定。在任一推進距離，采區地表垂直方向沉陷出現“中部大兩側小”的現象，這主要是由于采空區兩側為實體煤，對兩側覆巖運移具有限制作用，導致采空區中部覆巖以豎向直接下沉為主，兩側覆巖向中部傾斜。空區地表西部黃土坡沿垂直方向下沉大于東側，推測是由于西側山坡較為陡峭，可能出現黃土陡坡失穩滑落，影響地表構筑物安全。黃土“山脊”位移較黃土“山谷”位移大，則表明黃土溝谷地表的防護重點應在邊坡山脊處。

(a)推進50 m

4.2 地表水平移動量

工作面的回采會引起上覆地表產生水平移動，將其沿工作面長度和推進兩個方向進行分解，分別提取分布云圖，分析其水平移動變形特征。

1)工作面長度方向。工作面推進過程中，黃土溝壑地表沿工作面長度方向水平移動分布云圖如圖4所示。

由圖4可知，隨著工作面推進至50 m、150 m、250 m、350 m以及450 m時，工作面上方地表東側沿工作面長度方向最大水平位移值為0.67 m、2.37 m、2.92 m、3.00 m以及3.09 m。隨著工作面的推進，地表水平位移值不斷增加，但變形增速由0～50 m的13.28 mm/m，下降至350～450 m的0.78 mm/m，可見東側覆巖運移趨于穩定。工作面上方地表西側沿工作面長度方向最大水平位移值為0.86 m、3.12 m、3.84 m、3.92 m以及3.99 m，隨著工作面的推進，地表水平位移值不斷增加，但變形增速由0～50 m的17.26 mm/m，下降至350～450 m的0.71 mm/m，可見西側覆巖運移趨于穩定。同時，在工作面長度方向上，采區地表沉陷出現“中部小兩側大”的現象，這主要是由于采空區兩側為實體煤，對兩側覆巖運移具有限制作用，導致采空區中部覆巖以豎向直接下沉為主，兩側覆巖向中部傾斜，出現沿著工作面長度方向的水平位移。采空區地表西部黃土坡沿工作面長度方向水平位移大于東側，推測是由于西側山坡較為陡峭，可能出現黃土陡坡失穩滑落，影響地表構筑物安全。黃土“山脊”位移較黃土“山谷”位移大，則表明黃土溝谷地表的防護重點應在邊坡山脊處。

(a)推進50 m

2)工作面推進方向。工作面推進過程中，黃土溝壑地表沿工作面推進方向水平移動分布云圖如圖5所示。

(a)推進50 m

由圖5可知，當工作面推進至50 m、150 m以及250 m時，工作面上方地表沿工作面推進方向最大水平位移值為0.71 m、2.54 m以及2.31 m，可見隨著工作面的推進，地表水平位移值不斷增加并趨于穩定；同時工作面上方覆巖向采空區方向運移，導致地表出現超前于工作面推進距離的變形。由于西側黃土山較高，向采空區后方偏移量較大，已出現黃土陡坡滑塌危險。當工作面推進至350 m與450 m時，采空區后側上方地表出現向采空區方向偏移趨勢，偏移量分別為0.17 m及0.26 m，同時西側陡坡位移量增加，東側緩坡亦出現少量位移。

4.3 地表最大剪切應變

工作面推進過程中，黃土溝壑地表最大剪應變分布云圖如圖6所示。由圖6可知，隨著工作面推進至50 m、150 m、250 m、350 m以及450 m時，工作面上方地表剪切應變范圍不斷增加，說明地表失穩、破壞范圍逐漸增加。如圖6所示，采空區地表中部剪切應變較小，四周剪切應變量大，說明采空區四周圍巖均出現向采空區中部傾斜，造成地表出現盆地狀變形現象；同時采空區上覆地表形貌復雜，對剪切應變分布產生一定的影響，集中體現為地表山坡頂端剪切應變加大，坡中與坡底剪切應變較小，說明黃土坡整體穩定性較差，在工作面推進時以及推進后，應當對黃土坡采取相應治理措施，防止黃土坡整體翻覆影響地表結構。

(a)推進50 m

4.4 地表最大剪應變增量

工作面推進過程中，黃土溝壑地表最大剪應變增量分布云圖如圖7所示。由圖7可知，工作面推進50 m、150 m、250 m、350 m以及450 m時，采空區上方地表最大剪應變增量范圍不斷增加，且最大剪切應變增量集中分布于采空區地表中部溝底位置，說明采空區中部溝底應變量最大，最容易發生破壞，當工作面推進時應當注重對工作面中部溝底區域地表的維護與及時處理。同時可見，兩側黃土坡剪切應變增量集中于坡頂，表明隨著工作面推進，地表兩側黃土坡易出現整體傾覆，在工作面推進時以及推進后應當注重對黃土坡進行保護，防止出現黃土坡整體滑塌。

(a)推進50 m

綜合概率積分計算結果及數值模擬結果，由于13101工作面地表為黃土溝壑地貌，因此在概率積分計算時選取了煤層平均埋深作為計算參數，所獲得的值可以認為是地表移動變形的平均值。而在數值模擬部分，通過真實地表形貌特征建立了三維數值模型，獲取了13101地表移動變形的具體數據。本文概率積分法計算結果中，13101工作面地表平均沉陷值為5.48 m，平均水平移動值為1.37 m，在回采450 m數值模型中，于最大下沉值位置取平行于工作面的剖面，計算得到該剖面地表各節點下沉的平均值為5.59 m，同樣的方法可以計算得到水平移動平均值的數值模擬結果為1.64 m。因此，概率積分計算結果與數值模擬結果相互驗證，保證了研究結果的可靠性。

5 結論

以磁窯溝煤礦13101工作面為實際地表等高線及巖層賦存情況為依據，利用概率積分及數值模擬相結合的方法，研究了其回采引起的地表移動變形規律，分析了13101工作面回采后的黃土溝壑地表穩定性。

1)隨著工作面推進至50 m、150 m、250 m、350 m以及450 m時，工作面上方地表最大下沉量為1.50 m、5.08 m、6.12 m、6.41 m以及6.55 m，最大水平移動量分別為0.86 m、3.12 m、3.84 m、3.92 m以及3.99 m。

2)概率積分法計算結果中，地表平均下沉值為5.48 m，平均水平移動值為1.37 m；而數值模擬結果中，平均下沉值為5.59 m，平均水平移動值為1.64 m，兩種方法的分析結果相互驗證，保證了研究結果的可靠性。

3)受工作面采動影響，黃土溝壑地表剪切應變增量集中于坡頂，隨著工作面推進，地表兩側黃土坡易出現整體傾覆，黃土“山脊”位移較黃土“山谷”位移大，黃土溝壑地表的防護重點應在邊坡山脊處，防止出現黃土坡整體滑塌。