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

張 峰， 秦洪巖， 題正義

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

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

1 概率積分法簡介

1.1 方法特點

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

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

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

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

1.2 預計模型參數選擇

地表移動變形預計模型包含計算模型的構建和輸入參數選擇兩個方面.

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

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

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

(1)

(2)

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

(3)

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

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

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

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

(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為下沉系數；b為水平移動系數.

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

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

各主要預計參數的常規計算方法為：

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

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

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

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

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

1.3 預計方法存在的不足

由上述計算公式和輸入參數分析可知，概率積分法存在兩個方面的不足：一是輸入參數產生的誤差，二是計算模型自身產生的誤差.

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

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

1.3.2 計算模型自身產生的誤差 由于概率積分法基于隨機介質理論，計算模型中輸入參數沒有考慮待采工作面的巖性及組合特征，其計算結果中也必然會產生一定的誤差.

2 主要預計參數確定方法的建立及預計模型修正

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

2.1 主要預計參數確定方法的提出

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

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

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

W0=M-Hc，

(5)

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

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

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

(6)

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

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

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

(7)

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

2.2 預計模型修正

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

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

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

Wxz=Wyj+f(x),

(8)

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

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

傾斜值為下沉值的一階導數，得到傾斜值為：

(9)

曲率值為下沉值的二階導數，得到曲率值為：

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

(10)

由水平移動值根據與傾斜值之間的比例關系，得到水平移動值為：

(11)

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

(12)

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

(13)

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

3 現場實測驗證分析

3.1 現場實測

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

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

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

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

續表1

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

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

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

由于曲線具有微小的雙波峰、微小的雙低谷和中間平滑的特性，采用一般多項式擬合線性關系較差，因此，選用MATLAB中三角函數算法對均值曲線進行擬合，得出系數等級為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可以看出：利用多項式法修正后的計算模型預計已采工作面的下沉值與實測數據進行對比，下沉值的誤差范圍僅為0.1～0.4 m，極大的減小了由計算模型自身帶來的預計誤差，修正效果較好，能夠精準預計該井田內煤層開采帶來的地表移動變形.

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

4 結語

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

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

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