基于多項式法的地表移動變形計算模型修正

張 峰， 秦洪巖， 題正義

(1.山西工程技術(shù)學院采礦工程系， 山西 陽泉 045000；2.遼寧工程技術(shù)大學礦業(yè)學院， 遼寧 阜新 123000； 3.華北科技學院安全工程學院， 北京 101601)

煤層開采后，采場上覆巖層結(jié)構(gòu)經(jīng)歷了從移動變形到變形破壞的過程，巖層的變形破壞自下而上波及至地表，使得地表的地形地貌出現(xiàn)不同程度變形破壞，影響礦區(qū)自然生態(tài)結(jié)構(gòu)和構(gòu)筑物的安全.由于不同區(qū)域煤層上覆巖層結(jié)構(gòu)不同、巖石的性質(zhì)千差萬別，精確預(yù)計煤層開采后上覆巖層結(jié)構(gòu)和地表的移動變形過程一直是困擾學者研究地表構(gòu)筑物下煤層安全開采評價的主要難題，多數(shù)學者專家對此進行了一系列的研究.李世保[1]結(jié)合word和CAD研發(fā)了地表移動變形數(shù)據(jù)處理與預(yù)計一體化系統(tǒng)，克服了傳統(tǒng)方法只能查看，不能定量分析的不足.李青峰[2]對比分析數(shù)值模擬和概率積分法預(yù)計地表移動變形值的結(jié)果，認為數(shù)值模擬得到的結(jié)果較為準確.周婷婷[3]分析了地質(zhì)采礦條件對地表移動參數(shù)的影響機理，得到了各參數(shù)與地質(zhì)采礦條件之間的對應(yīng)關(guān)系，對該礦區(qū)地表移動變形預(yù)計提供了較高的實用價值.賈林剛[4]利用相似模擬方法研究了軟巖煤層組的采動力學演化特性、頂板破斷發(fā)展軌跡和地表移動變形規(guī)律，揭示了圍巖移動變形過程中的影響因素.韓永斌[5]利用現(xiàn)場監(jiān)測方法研究了厚煤層采動后地表移動變形大小和變形破壞規(guī)律，有效控制地表出現(xiàn)非連續(xù)變形.季鯤鵬[6-8]通過數(shù)值模擬和相似材料模擬研究了煤層開采后地表移動變形與預(yù)計參數(shù)之間的關(guān)系，對預(yù)計公式進行修正.綜上，大多數(shù)學者對于概率積分法預(yù)計地表移動變形存在誤差較大的問題具有共鳴之處，常采用數(shù)值模擬和參數(shù)修正等方法，沒有從不同地質(zhì)條件和計算公式推算過程上去開展相關(guān)研究工作，為此，依據(jù)地表下沉值與其它移動變形值之間的計算原理關(guān)系，結(jié)合煤層開采后地表移動變形計算結(jié)果和實測結(jié)果，本文提出了利用多項式法對計算模型進行修正的方法，構(gòu)建基于概率積分法和多項式法相疊加的地表移動變形值計算模型，對于構(gòu)筑物下煤層開采的安全性評價具有較高的實用價值.

1 概率積分法簡介

1.1 方法特點

概率積分法的基礎(chǔ)是隨機介質(zhì)理論[9-10]，按照其計算原理，共包含4個方面的基本假設(shè)：巖體為均質(zhì)、非連續(xù)介質(zhì)；具備線性疊加原理；彎曲帶內(nèi)巖體只發(fā)生形變，不發(fā)生體積的變化；開采結(jié)束后，地表的下沉體積等同于采區(qū)礦物的體積.

使用概率積分法對地表移動變形進行預(yù)計的過程為：將巖體看成特別小的顆粒介質(zhì)，把顆粒按埋深高度不同分層劃分，抽出下部顆粒，預(yù)測上部各分層中顆粒的移動規(guī)律，如圖1.

圖1 隨機介質(zhì)力學模型Fig.1 Stochastic medium mechanics model

假設(shè)顆粒之間完全失去聯(lián)系，可以相互運動，當某一層的顆粒介質(zhì)移動后，其上部的顆粒介質(zhì)有向下移動的趨勢，并伴有一定的概率，顆粒介質(zhì)的移動將形成一系列反應(yīng)，被影響到的顆粒介質(zhì)都會產(chǎn)生相應(yīng)的移動，最終將會影響到地表，引起地表的移動變形.

1.2 預(yù)計模型參數(shù)選擇

地表移動變形預(yù)計模型包含計算模型的構(gòu)建和輸入?yún)?shù)選擇兩個方面.

1.2.1 預(yù)計模型 將隨機介質(zhì)力學模型引入到煤層開采后的地表移動變形預(yù)計中，是將預(yù)計工作面的實際參數(shù)帶入到最大下沉值公式中，預(yù)計煤層開采條件下地表各點的下沉值.如圖2，以O(shè)1作為煤壁與采空區(qū)的分界線，當工作面推至煤壁O1點處時，采空區(qū)側(cè)頂板冒落穩(wěn)定后剩余自由下沉空間高度為W0(頂板彎曲下沉空間高度)，測得的地表最大下沉值也為W0，說明煤層開采達到充分采動.

圖2 走向主斷面地表下沉和水平移動Fig.2 Towards the main section surface subsidence and horizontal movement

根據(jù)隨機介質(zhì)理論，在X軸方向上任意點的地表下沉值和水平移動值表達式分別為：

(1)

(2)

式中,We(x)為X軸方向上任意點的地表下沉值，m；Ue(x)為X軸方向上任意點的水平移動值，m；r為主要影響半徑，m；B為比例常數(shù).

(3)

式中，H為走向主斷面采深，m；tanβ為走向主斷面主要影響角正切.

根據(jù)式(1)和式(2)，可繪制煤層開采后地表的下沉變化曲線和其對應(yīng)的地表水平移動值變化曲線，如圖3所示.

圖3 地表下沉曲線和水平移動曲線Fig.3 Surface subsidence and horizontal movement curve

對走向主斷面內(nèi)坐標為x的地表任意點A，其水平變形值、傾斜值和曲率值由地表水平移動值和下沉值表達式在[0，+∞]區(qū)間的積分及其一階或二階導(dǎo)數(shù)進行計算，預(yù)計公式為：

(4)

式中，W(x)為計算點x的下沉值，mm；i(x)為計算點x的傾斜值，mm·m-1；K(x)為計算點x的曲率值，mm·m-2；U(x)為計算點x的水平移動值，mm；ε(x)為計算點x的水平變形值，mm·m-1；W0為最大下沉值，mm；M為采高，m；α為煤層傾角，°；q為下沉系數(shù)；b為水平移動系數(shù).

1.2.2 輸入?yún)?shù) 根據(jù)公式(4)計算地表各點移動變形值所需的輸入?yún)?shù)，可將其分為待采工作面的地質(zhì)采礦基礎(chǔ)參數(shù)和預(yù)計參數(shù).其中工作面采高、煤層傾角為基礎(chǔ)參數(shù)，工作面的最大下沉值、最大下沉系數(shù)、水平移動系數(shù)、開采影響傳播角、主要影響角正切值、拐點偏移距和主要影響半徑為預(yù)計參數(shù).

基礎(chǔ)參數(shù)選擇依據(jù)待采工作面的開采設(shè)計可以得到，預(yù)計參數(shù)通常采用相鄰工作面直接引用、同一礦區(qū)多個工作面均值處理后引用、不同礦區(qū)相同地質(zhì)采礦條件下直接引用等工程類比法得到.

各主要預(yù)計參數(shù)的常規(guī)計算方法為：

1) 最大下沉值：最大下沉值依據(jù)該井田內(nèi)相鄰工作面的平均下沉系數(shù)，結(jié)合煤層傾角進行估算.

2) 水平移動系數(shù)：水平移動系數(shù)是煤層開采后地表最大水平移動值與最大下沉值之比，其值一般根據(jù)工程類比法得到，取值為0.3左右.

3) 開采影響傳播角：開采影響傳播角表示地表移動盆地向下山方向的偏移程度，主要用來確定下沉盆地拐點位置，與覆巖巖性、煤層傾角等因素有關(guān).與煤層傾角的表達關(guān)系式為θ=90°-kα，其中α為煤層傾角，k與覆巖巖性相關(guān)(且為小于1的常數(shù))，當工作面頂板覆巖巖性為堅硬時k為0.7～0.8，中硬時k為0.6～0.7，軟弱時k為0.5～0.6.

4) 主要影響角正切值：主要影響角正切值是主斷面邊界采深與主要影響半徑之比，與頂板上覆巖層的物理性質(zhì)、采區(qū)尺寸、煤層傾角有關(guān).

5) 拐點偏移距：拐點偏移距是自下沉曲線的拐點沿開采影響傳播角作直線，與煤層相交的交點沿煤層方向到采空區(qū)邊界的距離.

1.3 預(yù)計方法存在的不足

由上述計算公式和輸入?yún)?shù)分析可知，概率積分法存在兩個方面的不足：一是輸入?yún)?shù)產(chǎn)生的誤差，二是計算模型自身產(chǎn)生的誤差.

1.3.1 輸入?yún)?shù)產(chǎn)生的誤差 對于地表移動變形預(yù)計中需要輸入的計算參數(shù)為采高、埋深、煤層傾角，開采影響傳播角、主要影響角正切值、拐點偏移距、最大下沉系數(shù)、水平移動系數(shù)和任意點的坐標位置.其中采高、埋深、煤層傾角和兩個參數(shù)均為工作面已知定值，任意點的坐標是預(yù)計不同位置所需輸入的數(shù)值，也為定值，對于模型的計算結(jié)果不會產(chǎn)生誤差.

待采工作面的開采影響傳播角、拐點偏移距和主要影響角正切值均不能直接測量，主要由最大下沉值和水平移動值計算所得，而水平移動值的大小也與最大下沉值存在一定關(guān)系.最大下沉值是由下沉系數(shù)、采高和煤層傾角計算所得.對于待采工作面，其下沉系數(shù)需要通過工程類比得到，如果輸入的計算參數(shù)不夠準確，必然產(chǎn)生計算誤差.

1.3.2 計算模型自身產(chǎn)生的誤差 由于概率積分法基于隨機介質(zhì)理論，計算模型中輸入?yún)?shù)沒有考慮待采工作面的巖性及組合特征，其計算結(jié)果中也必然會產(chǎn)生一定的誤差.

2 主要預(yù)計參數(shù)確定方法的建立及預(yù)計模型修正

針對預(yù)計參數(shù)難以精準獲取和計算模型自身存在誤差問題[11]，通過實測結(jié)果與預(yù)計結(jié)果的誤差原因分析，提出相對準確的預(yù)計參數(shù)確定方法，并對預(yù)計模型進行修正，為地表移動變形規(guī)律分析提供可靠的手段.

2.1 主要預(yù)計參數(shù)確定方法的提出

根據(jù)基于概率積分法的地表移動變形預(yù)計計算公式可知，最大下沉值選擇的合理與否，直接影響著其它移動變形值的預(yù)計精度.為此，提出了覆巖剩余自由空間高度等于地表最大下沉值的觀點，并基于剩余自由空間高度與冒高、采高和冒落帶殘余碎脹高度的關(guān)系，構(gòu)建待采工作面最大下沉值的計算模型.

在采空區(qū)內(nèi)，煤層開采形成的空間首先由頂板垮落巖層所充滿[12]，然后冒落巖石在其上的裂隙帶和彎曲下沉帶的巖層重力作用下逐漸壓實、高度逐漸縮小.冒落巖石初始高度與壓實高度之差稱之為覆巖剩余自由空間高度，即為地表最大下沉值.因此，地表最大下沉值的實質(zhì)是地表在覆巖重力作用下能夠向下移動的垂直距離，其大小等于覆巖剩余自由空間高度.

根據(jù)上述定義，只要計算出殘余碎脹高度，根據(jù)采高，即可計算出地表最大下沉值.

W0=M-Hc，

(5)

式中，Hc為殘余碎脹高度，m.

冒落巖體的殘余碎脹高度計算公式為：

Hc=kiHm-Hm=(ki-1)Hm，

(6)

式中，ki為冒落巖體的殘余碎脹系數(shù)，Hm為冒落帶高度，m.

將式(6)代入式(5)中，得：

W0=M-Hc=M-(ki-1)Hm.

(7)

由于各巖層巖體的碎脹性系數(shù)和殘余碎脹性系數(shù)不相同，為了得到煤層開采后冒落巖體被壓實后的高度，首先要對冒落巖層發(fā)育的最高層位和最高位置進行判定，然后根據(jù)各個冒落巖層的高度和各個巖層巖體的殘余碎脹性系數(shù)，才能求得冒落巖體被壓實后的高度.

2.2 預(yù)計模型修正

為了削減計算誤差，根據(jù)相鄰工作面的實測數(shù)據(jù)，設(shè)想從最大下沉值的計算模型入手，在最大下沉值計算模型之后增加一修正函數(shù)項f(x)，對各移動變形值的預(yù)計進行修正，即采取多項式法進行修正[15-16].

具體方法是根據(jù)實際測量多個工作面的下沉值與預(yù)計結(jié)果的下沉值之間的差值，將多個差值的計算結(jié)果取平均值，對最終的均值曲線采用MATLAB軟件擬合成對應(yīng)的函數(shù)公式，添加到基礎(chǔ)參數(shù)修正后的下沉值預(yù)計結(jié)果中，完成概率積分法的下沉值計算模型修正.

下沉值的計算模型修正后如式(8).

Wxz=Wyj+f(x),

(8)

式中，Wxz為修正后的預(yù)計下沉值；Wyj為概率積分法基礎(chǔ)參數(shù)修正后的預(yù)計下沉值；f(x)為修正項.

根據(jù)下沉值的多項式修正模型在[0，+∞]區(qū)間的積分及其一階或二階導(dǎo)數(shù)進行計算，得到其它變形值的多項式修正計算模型.

傾斜值為下沉值的一階導(dǎo)數(shù)，得到傾斜值為：

(9)

曲率值為下沉值的二階導(dǎo)數(shù)，得到曲率值為：

K(x)=W″xz=W″yj+f″(x)=

(10)

由水平移動值根據(jù)與傾斜值之間的比例關(guān)系，得到水平移動值為：

(11)

由水平變形值與曲率之間的比例關(guān)系，得到水平變形值為：

(12)

根據(jù)各移動變形值的計算公式和修正后的計算公式、以及其一階導(dǎo)函數(shù)f′(x)、二階導(dǎo)函數(shù)f″(x)，結(jié)合地表最大下沉值的計算公式，得到修正后的各變形值計算公式[17-19].

(13)

將上述公式中的二維模型轉(zhuǎn)化為三維地表移動變形值計算模型，即可求取地表受采動影響的各變形值.

3 現(xiàn)場實測驗證分析

3.1 現(xiàn)場實測

以某礦S2S2工作面測點布置為例，沿工作面推進方向中線位置布置一條觀測線，測線距進回風巷道的垂直距離均為113.5 m，共計34個測點，間距為30 m，其中1號測點距切眼2 m，如圖4所示.

圖4 工作面測點布置示意圖Fig.4 Diagram of measuring point arrangement in working face

當工作面回采結(jié)束后，工作面上覆地表中各測點連續(xù)6個月內(nèi)累計下沉值小于30 mm時，即認為地表移動變形結(jié)束，工作面達到充分采動[20].測點1至測點34的移動變形值如表1所示.

表1 各測點的地表移動變形實測值Tab.1 Actual measured values of surface movement and deformation at each measuring point

續(xù)表1

3.2 多項式修正項f(x)計算

f(x)是由各測點實測下沉值與預(yù)計下沉值之間的差值求取，以該礦S2S2、N1S3、S2N1、N1S1共計4個綜放工作面為例，將預(yù)計結(jié)果與實測做差值，各差值的結(jié)果求取均值，將均值曲線利用MATLAB軟件回歸擬合，差值曲線圖如圖5.

圖5 4個工作面各測點實測下沉值與預(yù)計值之間差值曲線圖Fig.5 Curve of the difference between the actual measured subsidence value and the predicted value at each measuring point of four working faces

由于曲線具有微小的雙波峰、微小的雙低谷和中間平滑的特性，采用一般多項式擬合線性關(guān)系較差，因此，選用MATLAB中三角函數(shù)算法對均值曲線進行擬合，得出系數(shù)等級為7次，R2為0.996 8，擬合精度較高.擬合的曲線如圖6所示.

圖6 下沉值差值的均值擬合曲線圖Fig.6 The mean fitting curve of subsidence difference

f(x)=1.899sin(0.08489x+3.179)+0.1936sin(0.7251x-2.12)+0.2275sin(0.4684x-0.6443)+0.9147sin(0.3432x+1.774)+0.966sin(0.3165x-0.9018)+0.8867sin(1.186x-0.9285)+0.8808sin(1.203x+1.907).(14)

按照上述計算方法計算得到4個綜放工作面的下沉值，繪制下沉曲線對比圖(圖7).

由圖7可以看出：利用多項式法修正后的計算模型預(yù)計已采工作面的下沉值與實測數(shù)據(jù)進行對比，下沉值的誤差范圍僅為0.1～0.4 m，極大的減小了由計算模型自身帶來的預(yù)計誤差，修正效果較好，能夠精準預(yù)計該井田內(nèi)煤層開采帶來的地表移動變形.

圖7 綜放工作面實測下沉值與修正預(yù)計下沉值對比圖Fig.7 Contrast chart between actual measured subsidence value and revised predicted subsidence value of each mechanized caving face

4 結(jié)語

1) 由基于概率積分法的地表移動變形計算模型推算過程，確定預(yù)計方法存在的不足之處有兩個方面：輸入?yún)?shù)引起的誤差和模型自身產(chǎn)生的誤差.

2) 提出了利用多項式法對計算模型進行修正的方法，構(gòu)建基于概率積分法和多項式法相疊加的地表移動變形值修正計算模型.

3) 以下沉值為例，計算模型修正后得到的地表下沉值與實測結(jié)果之間的誤差范圍僅為0.1～0.4 m，極大的減小了計算模型的誤差范圍.