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

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

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

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

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

1.1 Knothe 時(shí)間函數(shù)

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

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

對(duì)式（1）積分，并對(duì)初始時(shí)間取0，初始沉陷值取0，可得:

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

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

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

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

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

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

則開(kāi)采第1 個(gè)單元體的動(dòng)態(tài)下沉值為：

第2 個(gè)開(kāi)采單元的動(dòng)態(tài)下沉值為：

則開(kāi)采至第n個(gè)單元的動(dòng)態(tài)下沉值為:

對(duì)于預(yù)計(jì)時(shí)刻，地表走向累積動(dòng)態(tài)下沉W(x,t)為：

W（x）概率積分法計(jì)算公式為:

式中：r為走向開(kāi)采影響半徑，m。

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

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

式中：R為傾向開(kāi)采影響半徑。

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

Usher 模型為描述增長(zhǎng)信息隨時(shí)間變化的模型[19]，其微分方程表達(dá)式為:

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

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

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

對(duì)式（14）進(jìn)行求導(dǎo)得：

而對(duì)于單點(diǎn)的下沉速度即為2 點(diǎn)下沉值和時(shí)間差值的比值，則有：

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

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

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

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

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

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

3 工作面及觀測(cè)站概況

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

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

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

4 數(shù)值模擬

4.1 模型構(gòu)建

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

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

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

4.2 走向下沉

走向下沉模擬圖如圖4。

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

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

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

4.3 傾向下沉

傾向下沉模擬圖如圖5。

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

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

5 地表移動(dòng)動(dòng)態(tài)下沉值預(yù)測(cè)

5.1 走向動(dòng)態(tài)下沉預(yù)測(cè)

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

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

圖6 走向動(dòng)態(tài)下沉預(yù)測(cè)圖Fig.6 Strike dynamic subsidence prediction

5.2 傾向動(dòng)態(tài)下沉預(yù)測(cè)

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

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

圖7 傾向動(dòng)態(tài)下沉預(yù)測(cè)圖Fig.7 Prediction of tendency dynamic subsidence

5.3 結(jié)果分析

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

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

6 結(jié) 語(yǔ)

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

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

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

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