基于FLAC3D和DEM數(shù)據(jù)的緩傾斜煤層開采沉陷分析

甘智慧，尚 慧，杜榮軍，占惠珠

(西安科技大學(xué) 地質(zhì)與環(huán)境學(xué)院，陜西 西安 710054)

近年來，隨著地下煤炭資源的采掘活動日益劇增，從而引發(fā)了一系列采空區(qū)地面沉陷問題，嚴(yán)重影響了社會經(jīng)濟發(fā)展和人類正常的生產(chǎn)生活。

眾多學(xué)者借助現(xiàn)場實測[1]、理論分析[2-3]、數(shù)值模擬[4-5]和物理相似模擬試驗[6-7]等方法對采煤引起的地面沉陷問題進行了大量研究。許家林等[8]采用物理和數(shù)值模擬方法，就覆巖主關(guān)鍵層對地表下沉動態(tài)過程的影響進行了研究；尹光志等[9]通過相似模型試驗與數(shù)值分析相結(jié)合的方法，對大傾角煤層深部開采的采巖移動基本規(guī)律、礦山壓力分布規(guī)律和地表移動變形規(guī)律進行了研究；孫學(xué)陽等[10-11]采用數(shù)值模擬與相似模擬實驗相結(jié)合的方法研究雙煤層開采對覆巖的破壞影響，結(jié)果顯示，工作面留設(shè)煤柱寬度越大，煤層開采對其覆巖的影響越小；王永國等[12-13]采用理論分析、現(xiàn)場實測、數(shù)值模擬、水文監(jiān)測等相結(jié)合的研究方法，深入研究覆巖運移特征及其主控影響因素；黃慶享等[14-15]采用物理模擬、數(shù)值計算和現(xiàn)場實測相結(jié)合的方法，研究淺埋煤層群開采覆巖及地表裂隙演化規(guī)律及機理。

綜上可知，煤層開采引起的地面沉陷問題多采用2 種或2 種以上研究方法，但利用多時相DEM數(shù)據(jù)進行采煤沉陷的時空演化規(guī)律分析鮮見報道。筆者以寧夏石嘴山礦區(qū)為對象，建立研究區(qū)三維地質(zhì)模型，模擬實際開采情況，分析不同開采期覆巖移動變形、應(yīng)力分布及塑性區(qū)變化規(guī)律，并將地表位移模擬結(jié)果與多時相DEM 數(shù)據(jù)疊加分析結(jié)果進行相互驗證，以期為緩傾斜煤層開采導(dǎo)致地面沉陷問題的研究提供新思路，為提高礦區(qū)開采沉陷監(jiān)測效率提供新方法。

1 石嘴山礦區(qū)概況

1.1 自然地理與地質(zhì)環(huán)境條件

寧夏石嘴山礦區(qū)位于黃河中游上段，寧夏回族自治區(qū)石嘴山市惠農(nóng)區(qū)境內(nèi)，東臨黃河，西靠賀蘭山，北與內(nèi)蒙古自治區(qū)阿拉善盟、烏海市伊克昭盟接壤，屬典型的大陸半干旱氣候，常年干燥多風(fēng)，日照充足、雨量稀少，蒸發(fā)強烈。

石嘴山礦區(qū)地貌屬于石嘴山侵蝕臺地，地勢西高東低；構(gòu)造位置位于東鄂博梁背斜以東，總體呈北東向較寬緩的向斜構(gòu)造，向斜軸自北東向南西蜿蜒伸展，略呈“S”型，兩翼地層不對稱，東南翼稍緩，地層傾角18°～32°，西北翼略陡，地層傾角13°～45°；出露地層自上而下依次為：第四系、新近系、三疊系、二疊系、石炭系、寒武系、薊縣系及長城系，主要含煤地層為石炭系上統(tǒng)太原組和二疊系下統(tǒng)山西組[16]，含煤地層總厚為185.36 m；水文地質(zhì)條件中等，地下水類型屬堅硬基巖裂隙水。

1.2 采礦條件及開采現(xiàn)狀

寧夏石嘴山礦區(qū)主要包括石嘴山一礦和石嘴山二礦(圖1)，為地下開采[17]，發(fā)育煤層9 層，煤層總厚為32.09 m，平均采深285 m。石嘴山礦區(qū)已有50余年開采歷史，地表變形嚴(yán)重，總體上形成7 個較大的塌陷坑，塌陷總面積為9.1 km2，最大塌陷深度為24.39 m[16]。采空區(qū)塌陷造成礦區(qū)附近房屋產(chǎn)生裂縫、有效耕地面積減少、井水河流污染等一系列地質(zhì)環(huán)境問題[3]，因此，對該礦區(qū)采煤引發(fā)地面沉陷問題進行研究已迫在眉睫。

本文選取石嘴山一礦Ⅱ—Ⅱ′剖面作為研究對象，該剖面沿煤層傾向橫跨3 號、4 號塌陷坑，具有一定代表性(圖1)。石嘴山一礦走向長約4.35 km，傾向?qū)捈s 1.19 km，面積 5.18 km2。開采高程為+1 055～+600 m，煤層平均傾角20°，為緩傾斜煤層，其中二、三、五、六、七、九號煤屬于穩(wěn)定可開采煤層，煤層覆巖主要為砂巖和砂質(zhì)頁巖(圖2)。根據(jù)井田內(nèi)煤層賦存條件、多煤層可分組特點，將可采煤層分為上下兩組，其中二、三煤層為上組煤，五、六、七、九煤層為下組煤。煤層開采采用斜井、分別集中運輸?shù)穆?lián)合開拓方式，分為3 個水平。一水平+1 050～+974 m；二水平+974～+830 m；三水平+830～+600 m。一水平已于1968 年回采完畢，二水平已于1989 年2 月回采完畢[16]，目前正在開采三水平，亞階段+680～+600 m。

圖1 石嘴 山礦區(qū)位置Fig.1 Location of Shizuishan Mining Area

圖2 石嘴山礦區(qū)Ⅱ—Ⅱ′工程地質(zhì)剖面Fig.2 Engineering geological section at row of Ⅱ-Ⅱ′ in Shizuishan Mining Area

2 FLAC3D 數(shù)值模擬

2.1 FLAC3D 簡介

FLAC(Fast Lagrangian Analysis of Continua)是由美國Itasca 公司開發(fā)的(連續(xù)介質(zhì)力學(xué)分析軟件)顯式有限差分程序[18]，采用混合離散法和動態(tài)松弛法，能較好地模擬地質(zhì)材料漸近破壞和失穩(wěn)，適用于模擬大變形、非線性問題；在模擬材料的屈服過程中，采用混合離散化方法模擬塑性破壞與塑性流動，比有限元有效[19]。

2.2 模型建立及參數(shù)選取

為研究石嘴山礦區(qū)地面沉陷及覆巖移動變形規(guī)律，以Ⅱ—Ⅱ′剖面為對象，運用FLAC3D軟件進行數(shù)值模擬分析。依據(jù)Ⅱ-Ⅱ′剖面上鉆孔所揭露的地層情況(圖2)、工作面范圍及采空區(qū)特點建立三維數(shù)值模型(圖3)。該區(qū)域地層結(jié)構(gòu)復(fù)雜，在數(shù)值模型建立過程中，對物理力學(xué)性質(zhì)相近的巖土體進行適當(dāng)簡化，將煤層按上下兩組簡化為兩層煤，平均厚度均為10 m，兩組煤同時開采，煤層開采遵循由上向下分3 水平開采的順序，最終確定模型尺寸長(X)×寬(Y)×高(Z)設(shè)置為2 200 m×600 m×850 m，包括79 560 個節(jié)點和86 121 個單元格。結(jié)合實際開采條件，石嘴山礦區(qū)煤層開采模擬時計算步驟如下：①單元網(wǎng)格化分利用FLAC3D5.0 Extrusion 建模；② 初始應(yīng)力場求解；③初始應(yīng)力位移清零；④ 根據(jù)礦體實際開挖順序分3 水平進行開采，并在各開采水平間預(yù)留煤柱15 m。

圖3 三維數(shù)值模型Fig.3 Three dimensional numerical model

模型基本假定為各向同性連續(xù)均勻介質(zhì)，力學(xué)模型采用莫爾–庫倫彈塑性模型。模型頂部地表為自由面，底部采用固定約束，四周采用水平位移約束。礦區(qū)初始應(yīng)力以自重應(yīng)力為主，設(shè)置重力加速度為10 m/s2，方向垂直向下，巖層物理力學(xué)參數(shù)見表1。

表1 巖層物理力學(xué)參數(shù)[16]Table 1 Physical and mechanical parameters of rock strata[16]

2.3 數(shù)值計算結(jié)果分析

2.3.1 應(yīng)力場分析

采礦工程破壞原巖應(yīng)力，周圍巖體應(yīng)力重新平衡，維持采空區(qū)穩(wěn)定的部分圍巖會產(chǎn)生局部應(yīng)力集中現(xiàn)象。隨著開采深度的增加，應(yīng)力逐漸增加，當(dāng)超過其強度極限時就會發(fā)生損傷、破壞甚至采空區(qū)整體失穩(wěn)。本模型分3 水平開挖，各水平開挖后采場周邊巖層的煤層走向方向和傾向方向的應(yīng)力分布情況如圖4 所示。在沿礦體走向方向中心處(Y=300 m)設(shè)置平行礦體傾向方向剖面，該剖面各水平開采后應(yīng)力場變化如下：

圖4 各開采水平傾向及走向方向應(yīng)力云圖Fig.4 Stress cloud maps of inclination and strike direction of each mining level

①第一水平開采應(yīng)力場分布呈現(xiàn)從上往下迭代遞增的變化趨勢，最大主應(yīng)力(壓應(yīng)力)均小于22.1 MPa。煤柱出現(xiàn)明顯的應(yīng)力集中，豎向壓應(yīng)力約為4 MPa。煤層頂板出現(xiàn)拉應(yīng)力，最大拉應(yīng)力出現(xiàn)在頂板中間偏上的位置(圖4a)；

② 隨著工作面推進，第二水平開采的采空區(qū)圍巖壓應(yīng)力擾動效應(yīng)較為明顯，采空區(qū)兩側(cè)煤柱出現(xiàn)高壓應(yīng)力集中，隨著與煤柱距離的增大，壓應(yīng)力逐漸變小，但變化梯度遠小于圍巖應(yīng)力變化梯度。采空區(qū)頂板出現(xiàn)拉應(yīng)力集中，受拉破壞并產(chǎn)生斷裂破碎、垮落(圖4b)；

③第三水平開采后，采空區(qū)圍巖最大拉應(yīng)力達到3.2 MPa。頂板拉應(yīng)力也不斷增大，頂板角隅處出現(xiàn)拉應(yīng)力集中，是影響采空區(qū)穩(wěn)定性的主要因素。煤柱應(yīng)力集中較為明顯，隨著開采深度增加，壓應(yīng)力增加明顯，最大壓應(yīng)力約5 MPa。拉應(yīng)力和壓應(yīng)力共同作用形成明顯的剪切破壞，以壓應(yīng)力為主，其最大值達到19.28 MPa。上下兩層煤層距離較近，采空區(qū)圍巖出現(xiàn)應(yīng)力疊加，表現(xiàn)出“群效應(yīng)”，應(yīng)力疊加區(qū)域的巖體更容易發(fā)生破壞(圖4c)。

針對沿礦體走向方向剖面(X=600 m)在Y軸方向上的應(yīng)力場進行分析。隨著開挖進程不斷推進，采空區(qū)體積不斷增大，圍巖應(yīng)力釋放明顯。走向應(yīng)力形態(tài)基本呈現(xiàn)對稱分布，且由上到下逐漸增大。

開采第一水平，右圍巖主要表現(xiàn)為壓應(yīng)力，最大值約為17.9 MPa，在采空區(qū)中央形成應(yīng)力等值線軌跡拱，角隅處出現(xiàn)應(yīng)力集中，采空區(qū)直接頂板也表現(xiàn)為應(yīng)力，其值為0.065 MPa(圖4d)；開采第二水平后上組煤頂板壓應(yīng)力由 0.065 MPa 減小為0.042 MPa，采空區(qū)頂板處于卸壓狀態(tài)(圖4e)；開采第三水平主要表現(xiàn)為壓應(yīng)力，最大值為17.6 MPa。采空區(qū)頂板應(yīng)力最小值都超過了頂板巖體的極限抗拉強度，頂板巖層出現(xiàn)拉張破壞，引起采空區(qū)頂板垮落(圖4f)。

2.3.2 位移場分析

礦體開采必然引起圍巖的變形和位移，隨著采空區(qū)面積和體積的增大，圍巖變形也持續(xù)增大，最終出現(xiàn)頂板垮落。各開采水平傾向及走向垂直位移數(shù)值模擬結(jié)果如圖5 所示。

從傾向方向的位移云圖(圖 5a—圖 5c)可以得出：煤層開挖后，煤層頂板巖層出現(xiàn)懸空，頂板下方?jīng)]有支撐體，其頂板巖層出現(xiàn)垮落、坍塌，遠離采空區(qū)的覆巖受礦層開采影響較小。隨著采空區(qū)推進，移動現(xiàn)象延伸至地表并形成沉陷盆地，其在地表的影響范圍遠大于采空區(qū)實際尺寸。沉陷盆地不對稱，上部巖層不僅向采空區(qū)移動，而且因為受重力作用由上山方向向下山方向位移。

圖5 各開采水平傾向及走向方向位移云圖Fig.5 Displacement cloud maps of inclination and strike direction of each mining level

一水平煤層開采時，采空區(qū)頂板垂直位移最大值–1.3 m，影響范圍在開采部位正上方；地表沉降最大值約0.6 m。上組煤底板出現(xiàn)隆起現(xiàn)象，其最大位移量為0.17 m，這是由于煤層回采后導(dǎo)致采空區(qū)應(yīng)力釋放，底板巖層受到向上的力，導(dǎo)致出現(xiàn)上拱現(xiàn)象(圖5a)。

二水平開采完成后，一水平最大變形位置的垂直位移由–1.3 m 增加到–5 m；隨著開采深度加大，最大位移出現(xiàn)在上組煤第二水平頂板中心偏上山方向，約–7.4 m。采空區(qū)底板隆起達到最大值，位移量為0.73 m，采空區(qū)中間煤柱受到的上、下圍巖擠壓作用最為嚴(yán)重。地表沉降明顯增大，地表沉降最大值從第一水平的 0.6 m 增加到第二水平的3.5 m(圖5b)。

開采第三水平后最大垂直位移位置與二水平一致，從二水平–7.4 m 增加到–19.1 m，其正上方地表沉降值為12 m，對地表沉降影響顯著；而三水平采空區(qū)頂板垂直位移–2～–8 m，對地表沉降影響較小，這是因為煤層三水平采深為 320～550 m，采高10～15 m，采高比均大于30(圖5c)。

對平行煤層走向剖面(X＝600 m)在Y軸方向上的位移場進行分析。第一水平開采導(dǎo)致頂板沉降，位移最大值約1.0 m，地表沉降約為0.8 m，煤層底板稍有隆起(圖5d)；隨著工作面繼續(xù)推進，開采第二水平時，底板出現(xiàn)隆起，最大底鼓量為0.09 m，頂板位移量比底板大，最大位移出現(xiàn)在上組煤采空區(qū)頂板中間偏上山方向，約為5.18 m。采空區(qū)周圍覆巖位移最大，往外距離采空區(qū)越遠位移越小，且圍巖位移方向都指向采空區(qū)(圖5e)；開采第三水平后最大沉降量為10 m，由于煤層開采引起采空區(qū)底板由三向受力狀態(tài)變?yōu)樗蕉蚴芰顟B(tài)，垂直方向應(yīng)力釋放引起底板回彈隆起，最大達到0.19 m。整體位移分布具有左右對稱性規(guī)律，采空區(qū)圍巖位移變化峰值均處在頂板及底板中間位置(圖5f)。

通過 FLAC3D數(shù)值模擬得到地表沉降等值線(圖6)。煤層回采后在地表形成2 個巨大的沉陷盆地，其中較小的沉陷區(qū)面積約為36 520 m2，最大沉降為7 m；另一個較大的地面沉陷區(qū)域呈“類W 型”，最大沉降量為12 m，影響范圍約259 380 m2，地表垂直位移略小于采空區(qū)最大沉降量，但地面沉陷范圍遠大于采空區(qū)范圍。地表垂直位移等值線圖近似呈現(xiàn)橢圓形，采空區(qū)地表影響區(qū)分別由2 個沉降中心向四周擴散，等值線逐漸稀疏。由于研究區(qū)為傾斜煤層，因此，地面沉陷盆地不對稱，2 個沉降中心均發(fā)生在沉陷盆地中部且偏下山方向，下山方向比上山方向影響范圍更大。

圖6 最終地表位移下沉等值線Fig.6 The final surface displacement subsidence contour

2.3.3 塑性區(qū)

塑性區(qū)形態(tài)可作為開挖擾動區(qū)的參考依據(jù)，因此，可基于數(shù)值模擬結(jié)果對圍巖開挖塑性區(qū)擴展進行分析[20]?；贔LAC3D模擬得到塑性區(qū)分布狀況(圖7)，模型單元在塑性狀態(tài)下將直接呈現(xiàn)出拉張或剪切破壞，可直觀反映圍巖穩(wěn)定性。通過塑性區(qū)分析可知，塑性區(qū)集中在煤層頂板、底板及煤柱。由于煤層傾斜，煤柱及鄰近采空區(qū)圍巖基本呈剪切破壞，地表及底板局部受到拉張破壞，表明剪切破壞對采空區(qū)的穩(wěn)定性影響較大。開采第一水平頂板出現(xiàn)貫通剪切破壞，頂板垂直位移較大，應(yīng)力集中明顯(圖 7a)；隨著煤層開采深度增加，塑性區(qū)向深度擴展，開采第二水平圍巖剪切破壞擾動范圍擴大，而拉張破壞影響范圍較小，主要發(fā)生在地表及煤層底板角隅處(圖 7b)。巷道底部煤層受到剪切破壞，導(dǎo)致開采層下部卸壓煤體瓦斯流入開采層采空區(qū)，應(yīng)對巷道下部的卸壓煤體區(qū)域重點進行瓦斯抽采工作；隨著掘進的深入塑性區(qū)不斷向圍巖深部擴展，開采第三水平塑性區(qū)影響范圍擴大到采空區(qū)地表，地面沉陷量達到最大值(圖7c)。

圖7 各開采水平塑性區(qū)分布Fig.7 The horizontal plastic zone distribution of each mining level

3 不同時相DEM 疊加分析與數(shù)值模擬結(jié)果對比

數(shù)字高程模型(Digital Elevation Model，DEM)作為地形表面形態(tài)的一種離散數(shù)字表達，能夠直觀、全面地反映地表位移變化情況。為實現(xiàn)對采空區(qū)塌陷坑深度及范圍監(jiān)測，分別選取礦區(qū)20 世紀(jì)70 年代(1︰50 000)和2003 年(1︰10 000)等高線生成DEM(圖8)，且將兩時相DEM 進行疊加分析，得到20 世紀(jì)70 年代至2003 年期間地表垂直位移變化圖(圖9)。小部分高程增加(即圖9 中紅色區(qū)域)，但多數(shù)是由于城鎮(zhèn)化建設(shè)造成。其他區(qū)域高程均有所降低，沉降幅度較大處均為塌陷坑，其中Ⅱ—Ⅱ′剖面穿過的3 號、4 號塌陷坑塌陷形態(tài)與數(shù)值模擬地表位移下沉范圍和形狀(圖6)基本一致。

圖8 石嘴山礦區(qū)兩時相DEMFig.8 Bi-temporal images DEM in Shizuishan Mining Area

圖9 石嘴山礦區(qū)20 世紀(jì)70 年代—2003 年沉降變化幅度Fig.9 The subsidence variation range chart of Shizuishan Mining Area from 1970s to 2003

為更進一步了解采空區(qū)地面塌陷垂直位移的變化情況，基于DEM 自動提取兩時相地形進行疊加對比，得到Ⅱ—Ⅱ′剖面地形線變化對比圖(圖10)。靠近惠農(nóng)區(qū)城區(qū)的0～1 000 m 處地形明顯上升，大多是由于在城鎮(zhèn)化建設(shè)中新增建(構(gòu))筑物造成；1 000～3 000 m 處地形出現(xiàn)明顯不均勻沉降，分別形成3 號和4 號塌陷坑。其中3 號塌陷坑范圍為167 m，豎向位移最大值為–4.2 m；4 號塌陷坑范圍為933 m，豎向位移最大值為–12.2 m。

將Ⅱ—Ⅱ′剖面基于兩時相DEM 差值運算結(jié)果與FLAC3D數(shù)模模擬得到的地表沉降曲線進行對比分析。由于兩沉降曲線不完全重合，選取FLAC3D模型原點為基準(zhǔn)點，即Ⅱ—Ⅱ′剖面地形線變化對比圖1 100 m 位置處(圖10)。數(shù)值模擬與DEM 分析結(jié)果對比如圖11 所示，其中模擬結(jié)果的3 號塌陷中心水平位置在240 m 處且垂直位移約為–7 m，而在30～70 m 位置引起地表隆起現(xiàn)象，其最大值為1.24 m，這是由于模型邊界效應(yīng)引起。兩時相DEM疊加統(tǒng)計得到3 號塌陷坑塌陷中心在60 m 處，垂直位移為–4.2 m。對比發(fā)現(xiàn)2 種方法所得3 號塌陷坑位置及沉降量相差較大，其主要原因是一水平煤層開采深度較淺，切眼接近地表，因此，存在一定的誤差。4 號塌陷坑處2 種方法確定的塌陷中心水平位置基本相同，沉降趨勢基本吻合，且兩者變形量也基本一致，相差均不超過3.08 m，大部分點位相差在1 m 以下。數(shù)值模擬得到的地面沉陷值整體略大于兩時相 DEM 數(shù)據(jù)垂直變化量，主要原因是DEM 數(shù)據(jù)來源分別是20 世紀(jì)70 年代和2003 年，而FLAC3D數(shù)值模擬的時間范圍是從1956 年規(guī)劃開采到礦山開采初期結(jié)束，F(xiàn)LAC3D的時間數(shù)據(jù)范圍遠大于DEM，其開采深度及開采范圍也更大。綜合以上分析可知，數(shù)值模擬計算的沉降量與兩時相DEM疊加統(tǒng)計分析的變化量結(jié)果及趨勢基本一致，計算結(jié)果合理。

圖10 Ⅱ—Ⅱ′剖面地形線變化對比圖(20 世紀(jì)70 年代—2003 年)Fig.10 Contrast of section Ⅱ-Ⅱ′ topographic line change contrast diagram(1970s-2003)

圖11 兩時相DEM 疊加分析與FLAC3D 數(shù)值模擬地表位移對比Fig.11 Comparison of surface displacement by bi-temporal images DEM superposition analysis and FLAC3D numerical Simulation

4 結(jié)論

a.通過FLAC3D模擬采煤過程，煤層開采后巖層主應(yīng)力呈現(xiàn)自上而下逐漸增大的變化趨勢，采空區(qū)上山方向出現(xiàn)拉應(yīng)力，最大值約3.2 MPa，煤柱出現(xiàn)明顯的應(yīng)力集中，最大應(yīng)力值達到5 MPa。

b.地面沉陷盆地不對稱，兩個塌陷坑沉降中心均發(fā)生在沉陷盆地中部且偏下山方向，下山方向比上山方向影響范圍大，2 個沉陷盆地最大沉降量分別為7 m 和12 m。

c.采空區(qū)圍巖以剪切破壞為主，塑性區(qū)范圍隨煤層開采深度增加逐漸擴展，主要發(fā)生在地表及采空區(qū)圍巖角隅處。

d.通過礦區(qū)20 世紀(jì)70 年代與2003 年DEM數(shù)據(jù)進行疊加分析得到：3 號塌陷坑最大垂直位移為–4.2 m，4 號塌陷坑最大垂直位移為–12.2 m。而數(shù)值模擬計算得到的兩塌陷坑最大沉降量分別為–7 m 和–12 m，兩者計算結(jié)果對比發(fā)現(xiàn)，3 號塌陷坑塌陷中心水平位置差異較大，4 號塌陷坑塌陷中心水平位置與沉降趨勢基本吻合。DEM 數(shù)據(jù)可與數(shù)值模擬相互驗證，對提高礦區(qū)沉陷監(jiān)測效率有一定參考價值。