喻成林 ,張宏貞 ,范洪冬
(1.徐州市自然資源和規劃局,江蘇 徐州 221008;2.中國礦業大學 環境與測繪學院,江蘇 徐州 221116)
地下資源開采后,破壞其原有平衡狀態,經過覆巖傳遞到地表,是1 個復雜的動態移動過程,涉及地質、采礦、空間等多方面因素,許多學者在開采沉陷動態變形原理、模型、精度等方面進行了研究,并取得了大量成果。開采沉陷動態預測模型主要有動態力學模型和時間函數模型2 類,目前主要采用時間函數構建開采沉陷的動態變形過程;常用的時間函數有Knothe 時間函數模型、雙曲線時間函數模型、Gompertz 時間函數模型、logistic 曲線時間函數模型及 Weibull 曲線時間函數模型等。
王悅漢等[1]將動態移動過程分為頂板初次垮落前、頂板垮落未充滿采空區、頂板垮落已充滿采空區、表土移動4 個階段,采動巖體分為連續介質、擬連續介質、非連續介質3 種介質,分階段建立了動態力學預測模型;吳侃等[2]提出了基于時間序列分析的動態預測模型;崔希民等[3]建立基于Knothe 時間函數的概率積分預計方法,根據臨界開采尺寸確定時間影響系數;張兵等[4-5]建立了1 種新的優化分段Knothe 時間函數,并提出了基于實測數據的反算時間函數對比求參法、基于概率積分參數的直接計算法的2 種優化分段Knothe 時間函數參數求取方法;劉玉成等[6-7]分析了下沉動態過程的Knothe 時間函數模型的不足,提出了冪指數的Knothe 時間函數,分析了常見地表沉陷時間函數模型的下沉、下沉速度、下沉加速度與時間關系,認為Weibull 曲線時間函數模型能完整地描述地表沉陷的動態過程;孫闖等[8]建立了考慮松散層下沉特點的雙因素Knothe 時間函數;王軍保等[9]借鑒巖石流變力學中非定常流變模型,構建了變時間影響系數的Knothe 的地表動態沉陷預測模型;許國勝等[10]回歸分析得到地表動態移動變形參數與地質及開采技術參數之間的關系;李春意等[11-12]進行了基于Logistic 時間函數、正態分布時間函數的動態沉陷預測;張永勝等[13]提出了以地表點最大下沉速度時刻為分界點,結合偏差改正、生長函數模型的1 種分段Weibull時間函數;劉東海[14]根據常村觀測資料,指出Weibull 時間函數參數不固定,兩參數變化規律為距邊界變化類偏態分布;JAROSZ 等[15]考慮開采產生的空間收斂,建立了雙參數Sroka-Schober 時間函數;KOWALSKI[16]在Knothe 和Sroka-Schober雙參數時間函數基礎上,提出了廣義時間函數;González Nicieza 等[17]在Knothe 時間函數基礎上,引入正態分布函數,建立了正態分布時間函數。
綜上,總體認為地表動態下沉量與時間的關系曲線呈“S”形曲線,下沉速度為0→Vmax→0 過程,下沉加速度為0→+amax→-amax→0 過程。
工作面上方的覆巖及地表在井下資源開采后,經歷彎曲、垮落、再次壓縮等過程,但在工作面不同位置的巖層完整厚度、位移、受力等經歷過程不同,根據工作面回采、覆巖受力與垮落情況,將工作面回采結束后的地表分為3 個部分,工作面動態分析分區圖如圖1。
圖1 工作面動態分析分區圖Fig.1 Partition diagram of dynamic analysis of working face
1)豎向整體受壓區(Ⅰ區)。處于工作面外側煤柱區域,工作面開采邊界至開采影響邊界之間。此部分區域在整個回采期間,覆巖經歷了彎曲,但沒有垮落、再次壓縮等過程,在豎向上是完整的。在整個開采期間豎向方向整體是受壓狀態。
2)中間區(Ⅱ區)。處于工作面內測靠近煤柱區域,工作面開采邊界至覆巖動態完全區邊界(覆巖垮落角確定邊界與工作面開采邊界之間)。此部分區域在整個回采期間,覆巖經歷了彎曲、垮落、非完全再次壓縮過程等部分或全部過程,主要是由于巖梁的存在導致垮落非充分及垮落巖石受力非充分,在豎向上部分破裂,且離工作面邊界(煤柱)越遠,破裂高度越大,即完整覆巖厚度減小。在整個回采期間主要經歷如下幾個階段:①受到影響到工作面推到下方前,在豎直方向上整體是受到壓縮變形,豎向受力增大;②推過后至覆巖垮落前,在豎直方向上整體是受到拉伸變形、豎向受力減小;③垮落后至移動穩定,在豎直方向上整體是受到壓縮變形,豎向受力增大至穩定。垮落巖石部分再次受到壓縮變形。
3)動態完全區(Ⅲ區)。處于工作面內測遠離煤柱區域,覆巖動態完全區(覆巖垮落角確定邊界至采空區中心)。此部分區域在整個回采期間,覆巖經歷了彎曲、垮落、完整再次壓縮過程,在豎向上部分破裂,破裂高度一致。在整個回采期間主要經歷如下幾個階段:①受到影響到工作面推到下方(垮落前),在豎直方向上整體是受到壓縮變形,豎向受力增大;②到再次受到壓縮前,在豎直方向上整體是受到拉伸變形,豎向受力減小;③到移動穩定,在豎直方向上整體是受到壓縮變形,豎向受力增大至穩定。
根據開采沉陷巖梁相關理論,地表下沉可簡化為各種類型巖梁彎曲,其彎曲下沉值與剛度(彈性模量與慣性矩之積)成反比,其中慣性矩與巖梁高度的三次方成正比。結合圖1,在同種情況下,上覆巖層剛度、慣性矩變化規律為:Ⅰ區(整體一致,最大)>Ⅱ區(與離開切眼距離成反比)>Ⅲ區(整體一致,最小)。
目前以時間函數為基礎的各種地表動態沉陷預測模型參數主要包括:最大變形值、時間函數(時間函數不同,其參數不同),各類參數與覆巖、開采等因素有關,體現其差異,但目前時間函數都采用相同參數,即在工作面不同位置(煤柱內外)的時間函數參數都一樣,這與處在工作面不同位置的受力過程、變形過程、巖層破壞情況及巖石力學性質等明顯不符。
Weibull 時間函數f(t)數學表達式為:
式中:k、m為參數,正數;t為時間。
據此建立開采沉陷過程中地表下沉時間函數、下沉速度時間函數、下沉加速度時間函數,開采沉陷動態下沉過程曲線如圖2。
圖2 Weibull 時間函數下沉、速度、加速度變化圖Fig.2 Weibull time function subsidence, velocity and acceleration variation diagram
式中:W0為該點下沉的最終值。
下沉速度時間函數v(t):
根據前面分析可知:工作面不同位置在回采期間,其覆巖移動、破壞、受力、巖梁剛度等不同,則其動態時間函數參數應不同,m、k決定了Weibull 時間函數的形態變化過程,因此m、k在工作面不同位置其存在差異,以垮落區邊界為原點、指向煤柱方向為正,根據前述巖梁彎曲特征簡述及相關巖石力學理論,建立其與地面點距離垮落區域邊界距離S的函數,工作面不同位置m、k變化情況圖如圖3。
圖3 工作面不同位置m、k 變化情況圖Fig.3 Changes of m and k at different positions of working face
在煤柱側區域(Ⅰ區)m、k為:
式中:a1、b1、c1、a2、b2為系數;S0為煤柱邊界距離;H為工作面采高。
從m、k計算公式及圖3 中可以看出:
1)在垮落區外側為負指數關系,以工作面邊界為界,此區域分為2 部分;Ⅰ區(煤柱范圍)m變化較小,基本為常數;Ⅱ區(中間區)m變化較大,在垮落區內側,m為一常數。
2)以工作面邊界、垮落區邊界為界,此區域分為3 部分;Ⅰ區(煤柱范圍)k為常數;Ⅱ區(中間區)m為線性變化;在垮落區內側,k為一常數。
由于不同位置上覆巖層經歷的受力、破壞、變形過程不同,根據開采沉陷相關研究,可采用如下方法確定垮落范圍d。① 方法1:采用覆巖力學理論,根據巖梁的相關破裂理論確定垮落邊界;② 方法2:采用概率積分法拐點偏移距;③ 方法3:采用頂板覆巖垮落角φ確定,d=H/tan φ。
某工作面采用長壁一次性采全高采煤法回采。走向為774 m,傾斜寬為220.7 m。工作面開采深度為470~515 m。根據走向觀測數據,采用上面的計算模型,獲得的反演模型系數見表1。
表1 反演模型系數Table 1 Coefficients of inversion model
根據工作面開采情況,以工作面內側90 m 為垮落邊界,以此確定各點距離邊界的S值,根據表1 反演動態時間函數 φ(t)如下:
根據工作面開采情況及參數確定情況,分別以開采結束后實測結果作為最終值、概率積分法預測值作為最終值,進行開采第523 d 的動態預測,實測結果、概率積分法預測值作為最終值的動態預測與實測對比如圖4(圖中距離以工作面邊界為原點,采空區方向為正)。
圖4 擬合與實測對比圖(523 d)Fig.4 Comparison between fitting and actual measurement (523 d)
以實測結果作為最終值:動態預測中誤差為110 mm,為最大下沉值的3.5%;以概率積分法預測結果作為最終值:動態預測中誤差為130 mm,為最大下沉值的4.2%。
1)分析了工作面上方不同位置的覆巖在開采過程中在豎向上的受力、變形、破壞等變化過程,根據煤柱、垮落情況將覆巖與地表劃分為豎向整體受壓區、中間區、動態完全區3 個區域。
2)以Weibull 時間函數動態預測模型為例,分析了m、k2 個參數的空間分布規律及特地,構建了2 個參數的分區數學模型及選取準則。
3)分別以實測結果、概率積分法預測結果作為工作面開采結束后的最終值,依據本文動態預測模型及參數確定方法,動態預測中誤差分別為110、130 mm,為最大下沉值的3.5%、4.2%,證明本模型的可靠性。