開采沉陷動(dòng)態(tài)預(yù)測時(shí)間函數(shù)參數(shù)變化規(guī)律研究
——以Weibull 為例

喻成林 ，張宏貞 ，范洪冬

（1.徐州市自然資源和規(guī)劃局，江蘇 徐州 221008；2.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116）

地下資源開采后，破壞其原有平衡狀態(tài)，經(jīng)過覆巖傳遞到地表，是1 個(gè)復(fù)雜的動(dòng)態(tài)移動(dòng)過程，涉及地質(zhì)、采礦、空間等多方面因素，許多學(xué)者在開采沉陷動(dòng)態(tài)變形原理、模型、精度等方面進(jìn)行了研究，并取得了大量成果。開采沉陷動(dòng)態(tài)預(yù)測模型主要有動(dòng)態(tài)力學(xué)模型和時(shí)間函數(shù)模型2 類，目前主要采用時(shí)間函數(shù)構(gòu)建開采沉陷的動(dòng)態(tài)變形過程；常用的時(shí)間函數(shù)有Knothe 時(shí)間函數(shù)模型、雙曲線時(shí)間函數(shù)模型、Gompertz 時(shí)間函數(shù)模型、logistic 曲線時(shí)間函數(shù)模型及 Weibull 曲線時(shí)間函數(shù)模型等。

王悅漢等[1]將動(dòng)態(tài)移動(dòng)過程分為頂板初次垮落前、頂板垮落未充滿采空區(qū)、頂板垮落已充滿采空區(qū)、表土移動(dòng)4 個(gè)階段，采動(dòng)巖體分為連續(xù)介質(zhì)、擬連續(xù)介質(zhì)、非連續(xù)介質(zhì)3 種介質(zhì)，分階段建立了動(dòng)態(tài)力學(xué)預(yù)測模型；吳侃等[2]提出了基于時(shí)間序列分析的動(dòng)態(tài)預(yù)測模型；崔希民等[3]建立基于Knothe 時(shí)間函數(shù)的概率積分預(yù)計(jì)方法，根據(jù)臨界開采尺寸確定時(shí)間影響系數(shù)；張兵等[4-5]建立了1 種新的優(yōu)化分段Knothe 時(shí)間函數(shù)，并提出了基于實(shí)測數(shù)據(jù)的反算時(shí)間函數(shù)對比求參法、基于概率積分參數(shù)的直接計(jì)算法的2 種優(yōu)化分段Knothe 時(shí)間函數(shù)參數(shù)求取方法；劉玉成等[6-7]分析了下沉動(dòng)態(tài)過程的Knothe 時(shí)間函數(shù)模型的不足，提出了冪指數(shù)的Knothe 時(shí)間函數(shù)，分析了常見地表沉陷時(shí)間函數(shù)模型的下沉、下沉速度、下沉加速度與時(shí)間關(guān)系，認(rèn)為Weibull 曲線時(shí)間函數(shù)模型能完整地描述地表沉陷的動(dòng)態(tài)過程；孫闖等[8]建立了考慮松散層下沉特點(diǎn)的雙因素Knothe 時(shí)間函數(shù)；王軍保等[9]借鑒巖石流變力學(xué)中非定常流變模型，構(gòu)建了變時(shí)間影響系數(shù)的Knothe 的地表動(dòng)態(tài)沉陷預(yù)測模型；許國勝等[10]回歸分析得到地表動(dòng)態(tài)移動(dòng)變形參數(shù)與地質(zhì)及開采技術(shù)參數(shù)之間的關(guān)系；李春意等[11-12]進(jìn)行了基于Logistic 時(shí)間函數(shù)、正態(tài)分布時(shí)間函數(shù)的動(dòng)態(tài)沉陷預(yù)測；張永勝等[13]提出了以地表點(diǎn)最大下沉速度時(shí)刻為分界點(diǎn)，結(jié)合偏差改正、生長函數(shù)模型的1 種分段Weibull時(shí)間函數(shù)；劉東海[14]根據(jù)常村觀測資料，指出Weibull 時(shí)間函數(shù)參數(shù)不固定，兩參數(shù)變化規(guī)律為距邊界變化類偏態(tài)分布；JAROSZ 等[15]考慮開采產(chǎn)生的空間收斂，建立了雙參數(shù)Sroka-Schober 時(shí)間函數(shù)；KOWALSKI[16]在Knothe 和Sroka-Schober雙參數(shù)時(shí)間函數(shù)基礎(chǔ)上，提出了廣義時(shí)間函數(shù)；González Nicieza 等[17]在Knothe 時(shí)間函數(shù)基礎(chǔ)上，引入正態(tài)分布函數(shù)，建立了正態(tài)分布時(shí)間函數(shù)。

綜上，總體認(rèn)為地表動(dòng)態(tài)下沉量與時(shí)間的關(guān)系曲線呈“S”形曲線，下沉速度為0→Vmax→0 過程，下沉加速度為0→+amax→-amax→0 過程。

1 地表動(dòng)態(tài)移動(dòng)過程覆巖變化特征

工作面上方的覆巖及地表在井下資源開采后，經(jīng)歷彎曲、垮落、再次壓縮等過程，但在工作面不同位置的巖層完整厚度、位移、受力等經(jīng)歷過程不同，根據(jù)工作面回采、覆巖受力與垮落情況，將工作面回采結(jié)束后的地表分為3 個(gè)部分，工作面動(dòng)態(tài)分析分區(qū)圖如圖1。

圖1 工作面動(dòng)態(tài)分析分區(qū)圖Fig.1 Partition diagram of dynamic analysis of working face

1）豎向整體受壓區(qū)（Ⅰ區(qū)）。處于工作面外側(cè)煤柱區(qū)域，工作面開采邊界至開采影響邊界之間。此部分區(qū)域在整個(gè)回采期間，覆巖經(jīng)歷了彎曲，但沒有垮落、再次壓縮等過程，在豎向上是完整的。在整個(gè)開采期間豎向方向整體是受壓狀態(tài)。

2）中間區(qū)（Ⅱ區(qū)）。處于工作面內(nèi)測靠近煤柱區(qū)域，工作面開采邊界至覆巖動(dòng)態(tài)完全區(qū)邊界（覆巖垮落角確定邊界與工作面開采邊界之間）。此部分區(qū)域在整個(gè)回采期間，覆巖經(jīng)歷了彎曲、垮落、非完全再次壓縮過程等部分或全部過程，主要是由于巖梁的存在導(dǎo)致垮落非充分及垮落巖石受力非充分，在豎向上部分破裂，且離工作面邊界（煤柱）越遠(yuǎn)，破裂高度越大，即完整覆巖厚度減小。在整個(gè)回采期間主要經(jīng)歷如下幾個(gè)階段：①受到影響到工作面推到下方前，在豎直方向上整體是受到壓縮變形，豎向受力增大；②推過后至覆巖垮落前，在豎直方向上整體是受到拉伸變形、豎向受力減?。虎劭迓浜笾烈苿?dòng)穩(wěn)定，在豎直方向上整體是受到壓縮變形，豎向受力增大至穩(wěn)定。垮落巖石部分再次受到壓縮變形。

3）動(dòng)態(tài)完全區(qū)（Ⅲ區(qū)）。處于工作面內(nèi)測遠(yuǎn)離煤柱區(qū)域，覆巖動(dòng)態(tài)完全區(qū)（覆巖垮落角確定邊界至采空區(qū)中心）。此部分區(qū)域在整個(gè)回采期間，覆巖經(jīng)歷了彎曲、垮落、完整再次壓縮過程，在豎向上部分破裂，破裂高度一致。在整個(gè)回采期間主要經(jīng)歷如下幾個(gè)階段：①受到影響到工作面推到下方（垮落前），在豎直方向上整體是受到壓縮變形，豎向受力增大；②到再次受到壓縮前，在豎直方向上整體是受到拉伸變形，豎向受力減?。虎鄣揭苿?dòng)穩(wěn)定，在豎直方向上整體是受到壓縮變形，豎向受力增大至穩(wěn)定。

根據(jù)開采沉陷巖梁相關(guān)理論，地表下沉可簡化為各種類型巖梁彎曲，其彎曲下沉值與剛度（彈性模量與慣性矩之積）成反比，其中慣性矩與巖梁高度的三次方成正比。結(jié)合圖1，在同種情況下，上覆巖層剛度、慣性矩變化規(guī)律為：Ⅰ區(qū)（整體一致，最大）>Ⅱ區(qū)（與離開切眼距離成反比）>Ⅲ區(qū)（整體一致，最小）。

目前以時(shí)間函數(shù)為基礎(chǔ)的各種地表動(dòng)態(tài)沉陷預(yù)測模型參數(shù)主要包括：最大變形值、時(shí)間函數(shù)（時(shí)間函數(shù)不同，其參數(shù)不同），各類參數(shù)與覆巖、開采等因素有關(guān)，體現(xiàn)其差異，但目前時(shí)間函數(shù)都采用相同參數(shù)，即在工作面不同位置（煤柱內(nèi)外）的時(shí)間函數(shù)參數(shù)都一樣，這與處在工作面不同位置的受力過程、變形過程、巖層破壞情況及巖石力學(xué)性質(zhì)等明顯不符。

2 Weibull 時(shí)間函數(shù)參數(shù)變化特點(diǎn)

Weibull 時(shí)間函數(shù)f(t)數(shù)學(xué)表達(dá)式為：

式中：k、m為參數(shù)，正數(shù)；t為時(shí)間。

據(jù)此建立開采沉陷過程中地表下沉?xí)r間函數(shù)、下沉速度時(shí)間函數(shù)、下沉加速度時(shí)間函數(shù)，開采沉陷動(dòng)態(tài)下沉過程曲線如圖2。

圖2 Weibull 時(shí)間函數(shù)下沉、速度、加速度變化圖Fig.2 Weibull time function subsidence, velocity and acceleration variation diagram

式中：W0為該點(diǎn)下沉的最終值。

下沉速度時(shí)間函數(shù)v(t)：

根據(jù)前面分析可知：工作面不同位置在回采期間，其覆巖移動(dòng)、破壞、受力、巖梁剛度等不同，則其動(dòng)態(tài)時(shí)間函數(shù)參數(shù)應(yīng)不同，m、k決定了Weibull 時(shí)間函數(shù)的形態(tài)變化過程，因此m、k在工作面不同位置其存在差異，以垮落區(qū)邊界為原點(diǎn)、指向煤柱方向?yàn)檎鶕?jù)前述巖梁彎曲特征簡述及相關(guān)巖石力學(xué)理論，建立其與地面點(diǎn)距離垮落區(qū)域邊界距離S的函數(shù)，工作面不同位置m、k變化情況圖如圖3。

圖3 工作面不同位置m、k 變化情況圖Fig.3 Changes of m and k at different positions of working face

在煤柱側(cè)區(qū)域（Ⅰ區(qū)）m、k為：

式中：a1、b1、c1、a2、b2為系數(shù)；S0為煤柱邊界距離；H為工作面采高。

從m、k計(jì)算公式及圖3 中可以看出：

1）在垮落區(qū)外側(cè)為負(fù)指數(shù)關(guān)系，以工作面邊界為界，此區(qū)域分為2 部分；Ⅰ區(qū)（煤柱范圍）m變化較小，基本為常數(shù)；Ⅱ區(qū)（中間區(qū)）m變化較大，在垮落區(qū)內(nèi)側(cè)，m為一常數(shù)。

2）以工作面邊界、垮落區(qū)邊界為界，此區(qū)域分為3 部分；Ⅰ區(qū)（煤柱范圍）k為常數(shù)；Ⅱ區(qū)（中間區(qū)）m為線性變化；在垮落區(qū)內(nèi)側(cè)，k為一常數(shù)。

由于不同位置上覆巖層經(jīng)歷的受力、破壞、變形過程不同，根據(jù)開采沉陷相關(guān)研究，可采用如下方法確定垮落范圍d。① 方法1：采用覆巖力學(xué)理論，根據(jù)巖梁的相關(guān)破裂理論確定垮落邊界；② 方法2：采用概率積分法拐點(diǎn)偏移距；③ 方法3：采用頂板覆巖垮落角φ確定，d=H/tan φ。

3 案例分析

某工作面采用長壁一次性采全高采煤法回采。走向?yàn)?74 m，傾斜寬為220.7 m。工作面開采深度為470～515 m。根據(jù)走向觀測數(shù)據(jù)，采用上面的計(jì)算模型，獲得的反演模型系數(shù)見表1。

表1 反演模型系數(shù)Table 1 Coefficients of inversion model

根據(jù)工作面開采情況，以工作面內(nèi)側(cè)90 m 為垮落邊界，以此確定各點(diǎn)距離邊界的S值，根據(jù)表1 反演動(dòng)態(tài)時(shí)間函數(shù) φ(t)如下：

根據(jù)工作面開采情況及參數(shù)確定情況，分別以開采結(jié)束后實(shí)測結(jié)果作為最終值、概率積分法預(yù)測值作為最終值，進(jìn)行開采第523 d 的動(dòng)態(tài)預(yù)測，實(shí)測結(jié)果、概率積分法預(yù)測值作為最終值的動(dòng)態(tài)預(yù)測與實(shí)測對比如圖4(圖中距離以工作面邊界為原點(diǎn)，采空區(qū)方向?yàn)檎?。

圖4 擬合與實(shí)測對比圖（523 d）Fig.4 Comparison between fitting and actual measurement (523 d)

以實(shí)測結(jié)果作為最終值：動(dòng)態(tài)預(yù)測中誤差為110 mm，為最大下沉值的3.5%；以概率積分法預(yù)測結(jié)果作為最終值：動(dòng)態(tài)預(yù)測中誤差為130 mm，為最大下沉值的4.2%。

4 結(jié) 語

1）分析了工作面上方不同位置的覆巖在開采過程中在豎向上的受力、變形、破壞等變化過程，根據(jù)煤柱、垮落情況將覆巖與地表劃分為豎向整體受壓區(qū)、中間區(qū)、動(dòng)態(tài)完全區(qū)3 個(gè)區(qū)域。

2）以Weibull 時(shí)間函數(shù)動(dòng)態(tài)預(yù)測模型為例，分析了m、k2 個(gè)參數(shù)的空間分布規(guī)律及特地，構(gòu)建了2 個(gè)參數(shù)的分區(qū)數(shù)學(xué)模型及選取準(zhǔn)則。

3）分別以實(shí)測結(jié)果、概率積分法預(yù)測結(jié)果作為工作面開采結(jié)束后的最終值，依據(jù)本文動(dòng)態(tài)預(yù)測模型及參數(shù)確定方法，動(dòng)態(tài)預(yù)測中誤差分別為110、130 mm，為最大下沉值的3.5%、4.2%，證明本模型的可靠性。