基于FLAC3D和DEM數據的緩傾斜煤層開采沉陷分析

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

(西安科技大學 地質與環境學院，陜西 西安 710054)

近年來，隨著地下煤炭資源的采掘活動日益劇增，從而引發了一系列采空區地面沉陷問題，嚴重影響了社會經濟發展和人類正常的生產生活。

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

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

1 石嘴山礦區概況

1.1 自然地理與地質環境條件

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

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

1.2 采礦條件及開采現狀

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

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

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

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

2 FLAC3D 數值模擬

2.1 FLAC3D 簡介

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

2.2 模型建立及參數選取

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

圖3 三維數值模型Fig.3 Three dimensional numerical model

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

表1 巖層物理力學參數[16]Table 1 Physical and mechanical parameters of rock strata[16]

2.3 數值計算結果分析

2.3.1 應力場分析

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

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

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

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

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

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

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

2.3.2 位移場分析

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

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

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

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

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

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

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

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

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

2.3.3 塑性區

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

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

3 不同時相DEM 疊加分析與數值模擬結果對比

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

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

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

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

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

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

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

4 結論

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

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

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

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