王文才 ，吳周康 ，王 鵬 ，蘇保山 ，楊少晨

（1.內(nèi)蒙古科技大學 礦業(yè)與煤炭學院，內(nèi)蒙古 包頭 014010；2.國能烏海能源五虎山礦業(yè)有限責任公司，內(nèi)蒙古 烏海 016000）

煤炭開采過程中，采區(qū)巖體內(nèi)部原有的力學平衡狀態(tài)被破壞，引起了地表的移動變形，并進一步容易誘發(fā)多種地質(zhì)災害[1]。工作面開采過程中采空區(qū)的形成使得巖體內(nèi)部形成空腔，原巖應力轉(zhuǎn)變?yōu)閿_動應力，且上覆巖層在自重以及其他巖層應力的疊加作用下，使得覆巖發(fā)生運移破壞。這種破壞隨著工作面的推進，由起初的覆巖移動發(fā)展為垮落，進一步發(fā)展至關鍵層形成破斷再進而發(fā)展至地表形成地表的沉陷[17]。地表沉陷的危害性尤其嚴重[2-3]。因此開展地表沉陷預測成為采礦工程中重要的一項工作[4-6]。針對這一問題，郭惟嘉等[7]通過現(xiàn)場與理論相結(jié)合的分析手段，對翟鎮(zhèn)煤礦采動下地表移動規(guī)律進行了分析，得出了地表非連續(xù)移動的誘發(fā)因素是覆巖的結(jié)構(gòu)與采深與采高的比值；劉寶琛等[8]引入概率積分法并進行改進，為我國地表移動規(guī)律研究指出了更加明朗的研究方向；劉詩杰等[9]將數(shù)值模擬與概率積分法配合使用，對由于工作面開采而造成的地表沉陷進行了準確的預測，有效保證了地表建筑物的安全；張世良[10]使用地表監(jiān)測點實測數(shù)據(jù)對概率積分法地表沉陷預測參數(shù)進行反演，從而做到對地表沉陷的二次預測；郭瑞瑞[11]通過對布爾臺礦區(qū)工作面開采沉陷規(guī)律的研究，得到地表下沉的特點為啟動速度快，穩(wěn)定態(tài)緩慢；韓亞鵬等[12]以觀測站實測數(shù)據(jù)為基礎，通過建立概率積分法時間與坐標的預測函數(shù)，對開采任意過程進行了地表沉陷的預測。

上述研究多集中于靜態(tài)的下沉預測，然而地表沉陷是1 個時間和空間上的動態(tài)過程，因此針對這一特征需做出動態(tài)下沉預測[13]。張兵等[14]、崔希民等[15]通過運用一系列的時間函數(shù)從不同的角度，對地表沉陷動態(tài)預測的理論與方法進行了改進與提高。然而針對不同的工況條件，需要對其做出不同的改進。為此，在以往學者的基礎上[16]，通過將 Knothe 時間函數(shù)、概率積分法以及 Usher 函數(shù)相結(jié)合，構(gòu)建出動態(tài)沉陷預測模型，并結(jié)合數(shù)值模擬與現(xiàn)場實測驗證了模型公式的實用性與準確性，可為其他工作面開采地表沉陷預測提供參考。

1 地表預測時間函數(shù)與概率積分交叉

1.1 Knothe 時間函數(shù)

Knothe 在分析工作面開采形成的地表下沉與時間上的關系時，對一假定點在時間為t時刻監(jiān)測到的下沉速度為dW（t）/dt與此點的最大地表沉陷值Wmax與動態(tài)沉陷函數(shù)W（t）的差值為反比，則有：

式中：A為與上覆巖層性質(zhì)有關的影響系數(shù)；Wmax為地表最大沉陷值，mm；t為沉陷發(fā)展時間，d。

對式（1）積分，并對初始時間取0，初始沉陷值取0，可得:

令時間函數(shù)f(t)=1-e-At，則式（2）可表示為：

1.2 時間函數(shù)結(jié)合概率積分法

走向主斷面沉陷原理分析圖如圖1。

圖1 走向主斷面沉陷原理分析圖Fig.1 Principle analysis diagram of strike main section settlement

假設工作面開采預計時間為t，即工作面從開切眼推進之日起，到工作面開采結(jié)束的時間差值；l為工作面開采長度，單位m；v為工作面推進速度，單位m/d。若有tv≥l，則表示工作面內(nèi)所有單元都參與地表下沉過程；若有tv

由概率積分法預測[18]可知，開采沉陷為階段累積變化過程，對于如圖1 的單一長壁式工作面開采的過程，可以將工作面按照推進過程劃分為n個單元：則對于第1 個階段，假設工作面推進速度為v1，工作面單元沉陷預計時間為t1，即有：l1=v1t1；對于第2 個階段則有l(wèi)2=v2t2；可知第n個開采階段，工作面的推進距離有l(wèi)n=vntn，而各個階段之間的持續(xù)時間為相鄰之間的差值，則第1 個單元一直持續(xù)到第n個單元地表下沉的持續(xù)時間為：

則開采第1 個單元體的動態(tài)下沉值為：

第2 個開采單元的動態(tài)下沉值為：

則開采至第n個單元的動態(tài)下沉值為:

對于預計時刻，地表走向累積動態(tài)下沉W(x,t)為：

W（x）概率積分法計算公式為:

式中：r為走向開采影響半徑，m。

同理可知，傾向累積動態(tài)下沉W(y,t)為：

W(y)概率積分法計算公式為：

式中：R為傾向開采影響半徑。

2 基于Usher 函數(shù)對動態(tài)下沉模型的構(gòu)建

Usher 模型為描述增長信息隨時間變化的模型[19]，其微分方程表達式為:

式中：z為模型函數(shù)；a為變化增長因子；b為形狀因子（取值范圍為0～1）；zm為極限值，又稱最值。

式（12）經(jīng)過變形積分后得到用于表示采空區(qū)地表動態(tài)沉陷預測的函數(shù)，其表達式為：

式中：W(t)為t時刻地表的沉陷量；c為沉陷參數(shù)；a為時間系數(shù)；b為下沉曲線形狀控制參數(shù)。

對式（14）進行求導得：

而對于單點的下沉速度即為2 點下沉值和時間差值的比值，則有：

將式（16）代入式（15）可得下沉速度與時間的函數(shù)v(t)為：

而由《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采指南》一書中可知，地表動態(tài)下沉時間函數(shù)積分W(x,y)t為[20]：

則引申出地表走向動態(tài)下沉函數(shù)W(x,t)為：

引申出地表傾向動態(tài)下沉函數(shù)W(y,t)為：

將式（8）、式（17）代入式（19）得：

同理將式（11）、式（17）代入式（20）得：

3 工作面及觀測站概況

1）工作面概況。研究礦井的工作面開采4#煤，埋深為234 m，采高為4.0 m，傾角為2°，煤層產(chǎn)狀由終采線向切眼方向傾斜，煤層為南北走向，設計傾向長為200 m，走向長為1 200 m?；夭山Y(jié)束后采空區(qū)實際傾向長為195 m，走向長為1 157 m。

2）觀測站布置情況。研究中采用十字交叉法布置測線，觀測線布置圖如圖2。共布設2 條觀測線，1 條為沿工作面傾向觀測線（C-D），長度為850 m；1 條為沿工作面走向觀測線（A-B），長度為1 050 m，該觀測線從采空區(qū)的外側(cè)開始，直至工作面內(nèi)部。同時在觀測線上按一定的方法要求布設若干個觀測點，測點之間間距為20 m。

在確定坐標的前提下，以多次觀測得出的數(shù)據(jù)進行計算分析。從2021 年2 月觀測站建成起，至2021 年12 月共進行14 次監(jiān)測，觀測間隔平均約為21 d。

4 數(shù)值模擬

4.1 模型構(gòu)建

為了獲得由于工作面開采而導致得地表移動下沉的清晰展示，使用FLAC3D數(shù)值模擬軟件，以某煤礦工程地質(zhì)條件為模擬基礎，建立立體，模型整體圖，模型如圖3。設計模型尺寸為500 m×1 000 m×300 m 模型采用Moore-Coulomb criterion模型。模型除頂部外其他7 個面固定。模型上部施加0.16 MPa 模型各巖性參數(shù)見表1。

表1 巖性參數(shù)表Table 1 Table of lithology parameters

圖3 模型整體圖Fig.3 Overall model

4.2 走向下沉

走向下沉模擬圖如圖4。

圖4 走向下沉模擬圖Fig.4 Simulation of subsidence

由圖4 可以看出，工作面推進200 m 時地表表現(xiàn)為隆起，且隆起值為649 mm，當推進至400 m 時由于工作面上覆巖層向采空區(qū)的垮落而表現(xiàn)為沉陷，此時地表隆起值減小至500 mm，于是得出隨著工作面的推進，地表下沉值沿著走向不斷增大，且在開采初期表現(xiàn)為地表隆起現(xiàn)象。由“砌體梁”理論可知，工作面開采后采空區(qū)內(nèi)側(cè)隨工作面的推進而向下破斷形成“倒三角”下沉堆疊，由于應力的集中，使得采空區(qū)外側(cè)對應地表位置發(fā)生局部隆起[21]。在工作面開采中期即工作面推進至600 m 范圍內(nèi)地表沉陷值急劇增大，由原本的隆起500 mm 左右，變?yōu)橄鲁? 840 mm，并且達到最大地表沉陷，此時工作面達到充分采動。之后的開采過程中工作面下沉速度減緩，直到工作面推進至950 m 即開采結(jié)束，地表仍舊發(fā)生緩慢下沉，表現(xiàn)為殘余下沉。

非充分開采階段走向下沉形狀表現(xiàn)為“V”形下窄上寬，且“V”形底隨工作面走向推進移動。充分開采階段走向下沉形狀表現(xiàn)為“U”形下平上寬，且“U”形底位置不隨工作面走向推進而發(fā)生移動。

4.3 傾向下沉

傾向下沉模擬圖如圖5。

圖5 傾向下沉模擬圖Fig.5 Tendency subsidence simulation diagram

由圖5 可以看出，傾向下沉在工作面推進400 m 以內(nèi)表現(xiàn)出與走向下沉相同的性質(zhì)，但地表隆起值比走向大80 mm 左右，當工作面推進至600 m 時地表表現(xiàn)為下沉，且下沉值為2 035 mm。下沉曲線表現(xiàn)出“V”形，隨工作面推進至950 m時地表傾向下沉達到最大值2 507 mm，其傾向下沉規(guī)律表現(xiàn)為：工作面推進200～400 m 時地表隆起降低，工作面推進400～600 m 時地表由隆起轉(zhuǎn)為沉陷且下沉速度激增，工作面推進600～950 m時地表沉陷值處于緩慢增長。

5 地表移動動態(tài)下沉值預測

5.1 走向動態(tài)下沉預測

根據(jù)觀測站觀測數(shù)據(jù)，確定走向Wmax=2 848.094 3 mm，r=450.675 3 m。代入式（21）進行擬合求得參數(shù)a=0.075，b=0.859，c=2 040 000。故走向動態(tài)預測下沉函數(shù)為：

根據(jù)式（23）由Origin 作出的走向下沉預測曲線如圖6。

圖6 走向動態(tài)下沉預測圖Fig.6 Strike dynamic subsidence prediction

5.2 傾向動態(tài)下沉預測

根據(jù)觀測站觀測數(shù)據(jù)，確定傾向Wmax=2 664.753 mm，R=390.245 3 m。代入式（22）進行擬合求得參數(shù)a=0.065，b=0.653，c=2 032 000。故傾向動態(tài)預測下沉函數(shù)為：

根據(jù)式（24）由Origin 作出的傾向下沉預測曲線如圖7。

圖7 傾向動態(tài)下沉預測圖Fig.7 Prediction of tendency dynamic subsidence

5.3 結(jié)果分析

采用數(shù)值模擬結(jié)果繪制的下沉曲線與實際曲線有一定的區(qū)別，從圖中可明顯發(fā)現(xiàn)，模擬下沉圖較實際測量曲線圖值偏低。說明模擬過程中地表移動條件對真實的地層反映具有不完善性。但是在驗證預測函數(shù)的準確性時具有一定的參考性。

經(jīng)過對走向、傾向預測函數(shù)進行的圖繪分析，結(jié)果表明該預測函數(shù)能夠準確地對走向以及傾向動態(tài)下沉進行精準定位，并且能夠有效地對下沉過程進行動態(tài)反映。

6 結(jié) 語

1）以Knothe 時間函數(shù)為基準，結(jié)合概率積分法做過程演化，最后聯(lián)立Usher 函數(shù)進行推導，分別得到了走向、傾向動態(tài)下沉預測函數(shù)模型。

2）通過使用FLAC3D數(shù)值模擬軟件構(gòu)建模型，得到了工作面開采動態(tài)下沉規(guī)律：開采初期，地表表現(xiàn)為隆起，且隆起量為 600 mm 左右；開采中期，地表表現(xiàn)為下沉，且在工作面推進距離為 600 m 時達到充分采動，此時下沉值為 2 850 mm，地表沉陷形狀表現(xiàn)為“U”型。

3）對預測函數(shù)進行參數(shù)擬合，得到了走向、傾向動態(tài)下沉值預測公式。并通過對公式進行繪圖與地表下沉模擬曲線、實測曲線對比分析結(jié)果顯示相似度很高，驗證了走向、傾向動態(tài)下沉預測公式的準確性。

4）通過對該公式的使用，能夠有效地對由于工作面充分開采而引發(fā)的地表下沉做出1 個動態(tài)的結(jié)果預測，進而有助于對地表建筑物等采取防護措施，對工作面是否需要采用充填開采方案做出判定。