近水平煤層開采地表移動邊界的預測研究

孫學陽 ，袁官思龍 ，苗霖田 ，杜榮軍 ，程 正

（1.西安科技大學 地質與環境學院，陜西 西安 710054；2.陜西省煤炭綠色開發地質保障重點實驗室，陜西 西安 710054；3.煤炭資源勘查與綜合利用重點實驗室，陜西 西安 710054）

獲取采煤沉陷特征值是研究煤礦開采地表移動變形規律和揭示其形成機理的基礎，現行規范和規程[1-2]根據采煤沉陷特征值對開采沉陷盆地邊界和損毀程度做了界定。理論分析方、相似材料模擬實驗、數值模擬是目前地表移動和變形預測的3 種方法。理論分析受到相關變量的限制，相似材料模擬實驗受到實驗條件的限制，20 世紀60年代R W Clough 首次將有限元引入土石壩的穩定性分析中，從此數值模擬技術成為了巖土工程領域主流的研究和設計方法[3]。因此對于采煤沉陷特征值的獲取和地表移動邊界的預測研究中，數值模擬具有優勢。

地表損毀邊界的主要影響因子為下沉、傾斜、水平變形。胡振琪等[4]以理論分析結合現場調查，提出了以耕地損毀臨界下沉值、臨界水平變形值和臨界傾斜值中的最小值圈定耕地損毀范圍；鄧偉男[5]利用FISH 語言代碼，提取模型主斷面上的下沉、水平移動、傾斜、水平變形及曲率值來繪制曲線圖，為數值模擬中數據的處理提供了參考；文獻[6-8]采用FLAC3D軟件進行模擬開采時，通過在走向和傾向主斷面上布設監測點的方式獲得其采煤沉陷特征值；LI 等[9]從應力場和位移場的角度出發，結合地表影響范圍的放大比，分析了地表沉降影響區隨采空區不斷擴大而發生的動態變化；WANG 等[10]研究了在FLAC3D建立含復雜地層邊界的三維數值計算模型的方法，并通過MATLAB 軟件導出地表沉降和水平移動數據，進而由計算公式轉化出傾斜、水平變形和曲率的等值線圖來直觀反映地表移動變形的影響范圍；蔡音飛等[11]利用Mathematica 軟件對FLAC3D模擬結果進行再次分析得到地表移動變形量，并通過列表繪圖函數ListPlot 繪制各次試驗的變形圖對地表受損情況進行評價。

基于此，從沉陷盆地全面積的角度，運用FISH 語言提取計算采煤沉陷特征值并進行運用，圈定了霍寶干河煤礦多工作面開采后地表移動盆地最外邊界、地表危險移動邊界、地表移動裂縫邊界；研究結果不僅為開采沉陷數值模擬數據后處理提供了新思路，還為生態恢復治理規劃提供參考。

1 地表移動變形數據的提取

地下巖土體的移動變形是采煤沉陷問題的核心,是地面沉陷的原因[12]。沉陷就是移動變形的總體描述，地表沉陷狀態按移動變形性質、方向、方式的不同進行描述，可分為垂直方向的移動變形和水平方向的移動變形。工程中通常采用下沉、傾斜、曲率，水平移動、水平變形等5 個采煤沉陷特征值進行定量評價[13]。無論是實測數據分析還是形變預計，采煤沉陷特征值都是基于下沉和水平位移2 個移動分量[14]。因此先將水平移動分量和垂直分量提取出來，再進行相應計算就能得到傾斜、水平變形和曲率。

在FLAC3D軟件中，力作用在單元中心坐標上，位移產生在網格節點上。針對這一特征，利用FISH 語言對六面體網格主導的計算模型中各個網格單元的特定節點進行遍歷，進行數據的提取并保存為Surfer 能夠識別的數據格式。

除此之外，FLAC3D中提供了豐富的內置函數，用戶可以運用它們進行簡單的程序設計，并得到相應結果。針對采煤沉陷問題，有必要先對FLAC3D中模型網格的節點編號和單元編號進行一下梳理，同時也是提取數據的代碼的編寫依據。

為了便于闡述，只建立了3× 3× 1 的模型，利用內置函數gp.pos(gp)、data.label.create(v)、data.label.text(lp)可以將單元編號和每個單元的節點編號在FLAC3D中顯示出來。FLAC3D中單元和節點編號示意圖如圖1，圖中：紅色數字為每個單元的編號；藍色數字為1 號單元的節點編號。規律如下：單元編號在走向上（x方向）依次排列，在傾向上（y方向）按等差數列排列；節點編號在單個單元上遵循右手螺旋法則。

圖1 FLAC3D 中單元和節點編號示意圖Fig.1 Schematic diagrams of unit and node number in FLAC3D

在FLAC3D中，對特定范圍的指定有2 種方式，自定義分組和自動分面；遍歷的對象也有2 種，一種是對單元進行遍歷，另外一種是對節點進行遍歷。

對于下沉值和水平移動值僅對節點進行遍歷即可；傾斜值和水平變形值先對單元進行遍歷，然后在具體編號的單元內通過節點指針找到相應編號節點的坐標和位移，曲率需要同時在2 個單元內通過節點指針找到相應編號節點的坐標和位移。最后，通過FISH 語言編寫出傾斜、水平變形和曲率的函數來提取數據并進行文本文件的保存，為繪制等值線提供數據基礎。

1.1 下沉值數據提取及保存

采用自動分面命令，模型表面將自動被命名為“Top i”，這就是下沉值的遍歷范圍。垂直移動分量為負，由下沉值正負號規定可知：正值表示測點下沉，負值表示測點上升，所以要在變量前面添加負號。另外，FLAC3D的默認單位為m，下沉值的單位需要轉換成mm，具體代碼如下：

1.2 水平移動值數據提取及保存

水平移動值遍歷范圍同上，它的正負號規定：指向礦層走向斷面走向方向的移動為正，反之為負；指向傾斜斷面上山方向的移動為正，下山方向的為負。水平移動值涉及傾向和走向兩個方向的分量，在FISH 變量中分別用dx 和dy 表示。單位為mm，具體代碼如下：

1.3 傾斜值數據提取及保存

傾斜值的計算涉及同一單元的相鄰兩節點，根據每個單元的節點編號規則，對提取傾斜值的代碼設計如下：走向上傾斜值為每個單元的8 號節點和6 號節點（或7 號節點和4 號節點）的下沉值之差與它們之間沿x軸水平距離的比值；傾向上傾斜值為8 號節點和7 號節點（或6 號節點和4號節點）的下沉值之差與它們之間沿y軸水平距離的比值[15]。

傾斜值反映了地表移動盆地沿某一方向的坡度，通常以i 表示，下列代碼中ix 代表走向上的傾斜值，iy 代表傾向上的傾斜值。它的正負號規定：指向礦層走向斷面走向方向的傾斜為正，反之為負；指向傾斜斷面上山方向的傾斜為正，下山方向的為負。單位需轉換為mm/m，具體代碼如下（以走向為例）：

1.4 水平變形值數據提取及保存

水平變形簡言之就是兩點間單位長度上的水平移動差，用相鄰兩點的水平移動差與兩點間水平距離的比值來表示。如圖1（b），走向上水平變形同樣用8 號和6 號節點（或7 號和4 號節點）的x軸位移分量之差與它們之間沿x軸水平距離的比值表示；傾向上水平變形用6 號和4 號節點（或8號和7 號節點）的y軸位移分量之差與它們之間沿y軸方向水平距離的比值表示。豎直方向的變化發生傾斜，水平方向的變化發生水平變形。所以參考提取傾斜值的代碼，將內置函數gp.disp.z()在走向上改成gp.disp.x()，在傾向上改成gp.disp.y()，節點編號和坐標函數對應修改即可。它的正負號規定：單位長度線段的拉伸變形為正，壓縮變形為負[16]，單位為mm/m。

1.5 曲率值數據提取及保存

曲率反映了地表垂直斷面上觀測線的彎曲程度，它可以用相鄰線段的傾斜差與兩線段中點間的水平距離之比來表示。

如圖1，走向上的曲率為相鄰單元8 號節點和6 號節點連成線段的傾斜值之差與兩線段x軸中點間水平距離的比值；傾向上的曲率為相鄰單元的6 號節點和4 號節點連成線段的傾斜值之差與兩線段y軸中點間水平距離的比值。它的正負號規定：上凸表示地表下沉曲線為正，下凹則為負。單位需轉換成mm/m2，具體代碼如下（以傾向為例）：

其中a、b、c、d為網格單元編號，e為沿走向方向的網格單元個數，g=a-e。

2 近水平煤層開采沉陷模擬

霍寶干河煤礦為近水平煤層，第四系黃土平均厚100 m，屬于厚松散層，煤層平均厚度為4.5m，開采深度為460～500 m。

模型尺寸為2 000 m ×2 050 m ×600 m，共分為8 層，平均采深為480 m，采厚4.5 m，煤層傾角為0°，共計116 850 個單元，132 396 個節點；數值計算模型分為3 個工作面進行開采模擬，分別是1101 首采工作面，1102 工作面和1103 工作面。3 個工作面尺寸均為600 m ×210 m，工作面之間均留設10 m 寬的護巷煤柱。模型剖切圖如圖2，巖層物理力學參數見表1。

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

由于地表移動觀測站只設置在1101 首采工作面，為了保證3 個工作面全采后數值模擬結果的正確性以及用FISH 語言提取特征值的合理性，有必要先將該工作面下沉盆地主斷面的特征值進行提取，1101 工作面主斷面特征值如圖3～圖5。

圖4 1101 工作面主斷面傾斜值和水平變形值Fig.4 Inclination values and horizontal deformation values of main section of 1101 working face

圖5 1101 工作面主斷面曲率值Fig.5 Curvature values of main section of 1101 working face

將圖3～圖5 中各特征值曲線的極值進行統計，1101 工作面地表移動變形量極值見表2。

表2 1101 工作面地表移動變形量極值Table 2 Extreme values of surface movement and deformation in 1101 working face

從表2 中可以看出用FISH 語言提取的各特征值與現場實測資料相吻合，驗證了該數值模擬的正確性及用FISH 語言提取地表移動變形數據的合理性。

3 應用實例

3.1 Surfer 繪圖流程

Surfer 軟件的主要優勢是利用其插值功能進行等值線圖的繪制。用它來繪制等值線圖的前提是數據文件格式要轉換成“.grd”文件格式[17]。

由于Surfer 的插值功能，即使數據不是等間距的，它依然能夠作圖。如果不進行數據的插值，只是繪制原始數據的圖，可以在“.text”轉“.grd”的對話框中選擇距離平方反比法或克里金插值法，本實例所有等值線圖的繪制過程皆采用克里金插值法。

運行提取采煤沉陷特征值的各段FISH 代碼會得到相應數據文件，按上述步驟操作之后新建等值線圖即可，按照項目的要求進行加密、著色等設置。除此之外，Surfer 中還可以添加邊界文件，可用于工作面輪廓線的繪制，為測量地表各邊界與采空區邊界線間的距離提供了便利。

3.2 地表移動邊界成圖

地表必然受到煤炭開采的影響，提前預判這一影響的1 個重要環節就是對變形后的地表進行等值線的繪制[18]。實測法獲取地表移動盆地角值參數存在困難，建立地表移動觀測站成本高，觀測數據匱乏；類比法獲取角值參數又受到地質采礦條件的制約；統計模型法通過統計經驗公式獲取角值參數的方法又受到實測數據和地質采礦條件的雙重制約。

綜上所述，隨著1101 工作面、1102 工作面和1103 工作面的開采，將整個模型通過FISH 代碼提取采煤沉陷數據并繪制等值線圖，從數值模擬的角度預測近水平煤層多工作面開采后地表移動邊界。地表下沉等值線圖如圖6，地表移動邊界等值線圖如圖7，圖中：藍色線框代表的是工作面輪廓，紅色和紫色等值線為臨界值等值線。

圖6 地表下沉等值線圖Fig.6 Contour map of surface subsidence

圖7 地表移動邊界等值線圖Fig.7 Contour diagrams of surface movement boundary

3.2.1 地表移動盆地最外邊界

地表移動盆地最外邊界以下沉10 mm 確定，該邊界在走向方向上距離采空區邊界518 m，在傾向方向上距離采空區邊界502 m。該邊界以內的區域都屬于開采沉陷影響區，這個邊界以外的區域一般認為不受開采沉陷的影響。

3.2.2 地表危險移動邊界

地表危險移動邊界按照現行規范[1]中對磚混結構建筑物臨界變形值的規定：傾斜變形值i=±3 mm/m，水平變形值ε=±2 mm/m，曲率K=±0.2 mm/m2來確定。

圖7（a）和圖7（b）中的傾斜臨界值等值線為紅色橢圓線，圖7（c）和圖7（d）走向上曲率最大值為0.045 6，傾向上曲率最大值為0.036 4 均未達到曲率臨界值，因此未能繪制出曲率的臨界等值線，圖7（e）和圖7（f）中紫色等值線為水平變形臨界值等值線；由此可以看出走向上地表危險移動邊界與采空區邊界左側相距299 m，右側相距291 m，傾向上的地表危險移動邊界與采空區邊界上側相距308 m，下側相距316 m。上下左右間距不相等的原因是由于工作面推進方向和工作面開采順序不同造成的。地表破壞的程度和范圍在等值線圖上的直觀反映就是傾斜和曲率大的區域[19]。

3.2.3 地表移動裂縫邊界

開采沉陷引起的土體拉伸破壞臨界值與土體變形值的定量關系如式（1）[20]：

式中： εxc為 土體產生裂縫的水平變形臨界值，mm/m；h為 土體埋深，m； ρ為平均密度，t/m3；φ為摩擦角，（°）； μ為泊松比；E為彈性模量，MPa；c為黏聚力，MPa。

εxc隨土體的埋深增加而增大，土體中任意一點的水平拉伸變形值εx≥εxc時，該處土體產生拉伸破壞，否則僅產生連續的沉陷變形。地表（h=0）處土體破壞極限狀態的水平變形臨界值εh=0為最小，可通過式（2）計算。

將黃土的物理力學參數E=10 MPa、μ=0.3、c=0.02 MPa、φ=15°代入式(2)得到地表產生拉伸變形的水平變形臨界值為εh=0=2.15 mm/m。

將該臨界值的等值線繪制于圖7（e）和圖7（f）中，以紅色等值線標示出，可以看出，走向上左側裂縫邊界距離采空區邊界498 m，右側裂縫邊界距離采空區邊界484 m；傾向上上側裂縫邊界距離采空區邊界476 m，下側裂縫邊界距離采空區邊界493 m。同樣，上下左右間距不相等的原因是由于工作面推進方向和工作面開采順序不同造成的。紅色等值線不僅確定了裂縫邊界，而且圈定了產生裂縫區域的范圍，如圖7（e）和圖7（f）。

4 結 語

1）依據單元和節點的編號規則，繪制曲率等值線圖的數據需要多個循環控制語句進行提取，FISH 代碼能夠在FLAC3D中實現采煤沉陷特征值的數據提取。

2）對霍寶干河煤礦進了數值模擬，提取1101工作面的采煤沉陷特征值曲線，結果與現場實測值吻合較好，驗證了數值模擬結果的可靠性。并提取了該工作面走向上的曲率值+0.007 mm/m2和-0.014 mm/m2，傾向上的曲率值+0.01 mm/m2和-0.03 mm/m2。

3）工作面推進方向和開采順序的不同會導致地表移動盆地邊界距離采空區邊界左右間距不相等。4）通過FISH 語言編程提取3 個工作面開采后地表移動變形數據，繪制的等值線圖預測了地表移動盆地最外邊界、地表危險移動邊界、地表移動裂縫邊界和裂縫區域，從數值模擬的角度為預測地表移動盆地邊界提供了新思路。