基于多函數交叉的地表動態沉陷預測模型及應用

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

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

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

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

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

1.1 Knothe 時間函數

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

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

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

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

1.2 時間函數結合概率積分法

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

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

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

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

則開采第1 個單元體的動態下沉值為：

第2 個開采單元的動態下沉值為：

則開采至第n個單元的動態下沉值為:

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

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

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

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

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

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

2 基于Usher 函數對動態下沉模型的構建

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

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

式（12）經過變形積分后得到用于表示采空區地表動態沉陷預測的函數，其表達式為：

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

對式（14）進行求導得：

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

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

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

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

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

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

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

3 工作面及觀測站概況

1）工作面概況。研究礦井的工作面開采4#煤，埋深為234 m，采高為4.0 m，傾角為2°，煤層產狀由終采線向切眼方向傾斜，煤層為南北走向，設計傾向長為200 m，走向長為1 200 m。回采結束后采空區實際傾向長為195 m，走向長為1 157 m。

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

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

4 數值模擬

4.1 模型構建

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

表1 巖性參數表Table 1 Table of lithology parameters

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

4.2 走向下沉

走向下沉模擬圖如圖4。

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

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

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

4.3 傾向下沉

傾向下沉模擬圖如圖5。

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

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

5 地表移動動態下沉值預測

5.1 走向動態下沉預測

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

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

圖6 走向動態下沉預測圖Fig.6 Strike dynamic subsidence prediction

5.2 傾向動態下沉預測

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

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

圖7 傾向動態下沉預測圖Fig.7 Prediction of tendency dynamic subsidence

5.3 結果分析

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

經過對走向、傾向預測函數進行的圖繪分析，結果表明該預測函數能夠準確地對走向以及傾向動態下沉進行精準定位，并且能夠有效地對下沉過程進行動態反映。

6 結 語

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

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

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

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